[Bio] / FigKernelScripts / dp.pl Repository:
ViewVC logotype

Diff of /FigKernelScripts/dp.pl

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1, Thu Jun 11 19:35:46 2009 UTC revision 1.4, Fri Jul 24 01:56:49 2009 UTC
# Line 7  Line 7 
7  #  #
8  #  Here the "dot product" will return a value from 0 to 1, where the higher the value  #  Here the "dot product" will return a value from 0 to 1, where the higher the value
9  #  the more similar the samples.  #  the more similar the samples.
10    #
11    #  The command will work if either file contains a single line containing a fabric.
12    #
13    
14  use Carp;  use Carp;
15  use Data::Dumper;  use Data::Dumper;
16  use strict;  use strict;
17    
18  my($N,$sample1,$sample2);  my($N,$sample1,$sample2);
19  my $usage = "usage: dp [-k n] Sample1 Sample2";  my $usage = "usage: dp [-k n] [-j] Sample1 Sample2";
20    
21    my $jaccard = 0;
22  $N = 4;  $N = 4;
23  while ( @ARGV && $ARGV[0] =~ s/^-// )  while ( @ARGV && $ARGV[0] =~ s/^-// )
24  {  {
25      $_ = shift;      $_ = shift;
26      if ( s/^k// ) { $N    = $_ || shift; next }      if ( s/^k// ) { $N    = $_ || shift; next }
27        elsif   ( s/^j// ) { $jaccard = 1 }
28      else      else
29      {      {
30          print STDERR "Bad flag '$_'\n", $usage;          print STDERR "Bad flag '$_'\n", $usage;
# Line 37  Line 42 
42  my $v1       = &load_sample($sample1,$role_num,$N);  my $v1       = &load_sample($sample1,$role_num,$N);
43  my $v2       = &load_sample($sample2,$role_num,$N);  my $v2       = &load_sample($sample2,$role_num,$N);
44    
45  my $dot_prod = &dot_prod($v1,$v2);  # note we compute the dot prod of two unit vectors created  my $dot_prod = &similarity($v1,$v2,$jaccard);  # note we compute the dot prod of two unit vectors created
46                                      # by normalizing v1 and v2                                      # by normalizing v1 and v2
47    
48  print "$dot_prod\n";  print "$dot_prod\n";
# Line 67  Line 72 
72      {      {
73          chomp;          chomp;
74          my @flds = split(/\t/,$_);          my @flds = split(/\t/,$_);
75            if (@flds > 8000)
76            {
77                @$role_counts = @flds;
78            }
79            else
80            {
81          my $func = $flds[$N-1];          my $func = $flds[$N-1];
82          if ($func)          if ($func)
83          {          {
# Line 81  Line 92 
92              }              }
93          }          }
94      }      }
95        }
96        close(IN);
97      return $role_counts;      return $role_counts;
98  }  }
99    
# Line 109  Line 122 
122      return $uv;      return $uv;
123  }  }
124    
125  sub dot_prod {  sub similarity {
126      my($v1,$v2) = @_;      my($v1,$v2,$jaccard) = @_;
127    
128      $v1 = &unit_vector($v1);      $v1 = &unit_vector($v1);
129      $v2 = &unit_vector($v2);      $v2 = &unit_vector($v2);
130    
131        if ($jaccard)
132        {
133            my $dp = &dot_prod($v1,$v2);
134            my @v3 = @$v1;
135            my $i;
136            for ($i=0; ($i < @$v2); $i++)
137            {
138                if ($v3[$i] && $v2->[$i])
139                {
140                    $v3[$i] += $v2->[$i];
141                }
142                elsif ($v2->[$i])
143                {
144                    $v3[$i] = $v2->[$i];
145                }
146            }
147            my $sq = &dot_prod(\@v3,\@v3);
148            return $dp / $sq;
149        }
150        else
151        {
152            return &dot_prod($v1,$v2);
153        }
154    }
155    
156    sub dot_prod {
157        my($v1,$v2) = @_;
158    
159    
160      my $tot = 0;      my $tot = 0;
161      my $i;      my $i;
162      for ($i=0; ($i < @$v1); $i++)      for ($i=0; ($i < @$v1); $i++)
163      {      {
164          if (defined($v1->[$i]) && defined($v2->[$i]))          if ($v1->[$i] && $v2->[$i])
165          {          {
166              $tot += $v1->[$i] * $v2->[$i];              $tot += $v1->[$i] * $v2->[$i];
167          }          }

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.4

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3