[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.2, Fri Jun 12 16:03:03 2009 UTC revision 1.3, Thu Jun 18 14:42:48 2009 UTC
# Line 47  Line 47 
47  use strict;  use strict;
48    
49  my($N,$sample1,$sample2);  my($N,$sample1,$sample2);
50  my $usage = "usage: dp [-k n] Sample1 Sample2";  my $usage = "usage: dp [-k n] [-j] Sample1 Sample2";
51    
52    my $jaccard = 0;
53  $N = 4;  $N = 4;
54  while ( @ARGV && $ARGV[0] =~ s/^-// )  while ( @ARGV && $ARGV[0] =~ s/^-// )
55  {  {
56      $_ = shift;      $_ = shift;
57      if ( s/^k// ) { $N    = $_ || shift; next }      if ( s/^k// ) { $N    = $_ || shift; next }
58        elsif   ( s/^j// ) { $jaccard = 1 }
59      else      else
60      {      {
61          print STDERR "Bad flag '$_'\n", $usage;          print STDERR "Bad flag '$_'\n", $usage;
# Line 71  Line 73 
73  my $v1       = &load_sample($sample1,$role_num,$N);  my $v1       = &load_sample($sample1,$role_num,$N);
74  my $v2       = &load_sample($sample2,$role_num,$N);  my $v2       = &load_sample($sample2,$role_num,$N);
75    
76  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
77                                      # by normalizing v1 and v2                                      # by normalizing v1 and v2
78    
79  print "$dot_prod\n";  print "$dot_prod\n";
# Line 151  Line 153 
153      return $uv;      return $uv;
154  }  }
155    
156  sub dot_prod {  sub similarity {
157      my($v1,$v2) = @_;      my($v1,$v2,$jaccard) = @_;
158    
159      $v1 = &unit_vector($v1);      $v1 = &unit_vector($v1);
160      $v2 = &unit_vector($v2);      $v2 = &unit_vector($v2);
161    
162        if ($jaccard)
163        {
164            my $dp = &dot_prod($v1,$v2);
165            my @v3 = @$v1;
166            my $i;
167            for ($i=0; ($i < @$v2); $i++)
168            {
169                if ($v3[$i] && $v2->[$i])
170                {
171                    $v3[$i] += $v2->[$i];
172                }
173                elsif ($v2->[$i])
174                {
175                    $v3[$i] = $v2->[$i];
176                }
177            }
178            my $sq = &dot_prod(\@v3,\@v3);
179            return $dp / $sq;
180        }
181        else
182        {
183            return &dot_prod($v1,$v2);
184        }
185    }
186    
187    sub dot_prod {
188        my($v1,$v2) = @_;
189    
190    
191      my $tot = 0;      my $tot = 0;
192      my $i;      my $i;
193      for ($i=0; ($i < @$v1); $i++)      for ($i=0; ($i < @$v1); $i++)
194      {      {
195          if (defined($v1->[$i]) && defined($v2->[$i]))          if ($v1->[$i] && $v2->[$i])
196          {          {
197              $tot += $v1->[$i] * $v2->[$i];              $tot += $v1->[$i] * $v2->[$i];
198          }          }

Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3