#!/usr/bin/perl -w
use strict;


# get the next array
sub nextgen($) {
  my $array_ref = shift;
  my $i = (@$array_ref) - 1;
  #print join("", @$array_ref); print "\n";
  while ( ($i > 0) && !(@$array_ref[$i-1] eq '1' && (@$array_ref[$i] eq '0'))) {$i--;}
  if ($i == 0) {
    return 0;
  }
  @$array_ref[$i-1,$i] = @$array_ref[$i,$i-1];
  my $u = $i+1;
  my $v = (@$array_ref)-1;
  while ($u<$v) {
    if (@$array_ref[$v] eq '0') {
      $v--;
    } elsif (@$array_ref[$u] eq '1') {
      $u++;
    } else {
      @$array_ref[$v,$u] = @$array_ref[$u,$v];
    }
  }
  return 1;
}

# shuffle an array
sub shuffle($) {
    my $array_ref = shift;
    my $i = @$array_ref;
    while ( --$i ) {
        my $j = int rand( $i+1 );
        @$array_ref[$i,$j] = @$array_ref[$j,$i];
    }
}

# generate an alignment of length l with i '1' and (l-i) '0'
sub generate_alignment($$) {
    my ($l,$i) = @_;
    my $r = ('1' x $i).('0' x ($l-$i));
    my @r_a = split(//, $r);
    shuffle(\@r_a);
    $r = join('', @r_a);
    return $r;
}

# compute the multihit for one/several seeds
sub multihit_seed_matches_alignment($$){
  my ($alignment_ref,$seeds_ref) = @_;
  my $result = 0;
  for (my $s=0; $s < (@$seeds_ref); $s++) {
    my $seed_ref = @$seeds_ref[$s];
    for (my $i = 0; $i < (@$alignment_ref) - (@$seed_ref) + 1; $i++) {
      for (my $j = 0; $j < (@$seed_ref) ; $j++) {
        if ((@$alignment_ref[$i+$j] eq '0') && (@$seed_ref[$j] eq '#')) {
          goto no_match;
        }
      }
      $result++;
    no_match:
    }
  }
  return $result;
}

# compute the multihit for one/several seeds
sub coveragehit_seed_matches_alignment($$) {
  my ($alignment_ref,$seeds_ref) = @_;
  my @coverage = ();
  @coverage = (@coverage, ('0') x (@$alignment_ref));
  for (my $s=0; $s < (@$seeds_ref);$s++) {
    my $seed_ref = @$seeds_ref[$s];
    for (my $i = 0; $i < (@$alignment_ref) - (@$seed_ref) + 1; $i++) {
      for (my $j = 0; $j < (@$seed_ref) ; $j++) {
        if ((@$alignment_ref[$i+$j] eq '0') && (@$seed_ref[$j] eq '#')) {
          goto no_match;
        }
      }
      for (my $j = 0; $j < (@$seed_ref) ; $j++) {
        if (@$seed_ref[$j] eq "#") {
          $coverage[$i+$j] = 1;
        }
      }
    no_match:
    }
  }
  #print join(",", @coverage);
  my $result = eval join('+', @coverage);
  return $result;
}





# main
my $length = 32;
my $seedfile = "seeds.txt";

# foreach seed
open S,$seedfile;
while (<S>) {
  chomp $_;
  if ($_ eq ""){
    next;
  }
  my $seeds = "$_";

  #>> correlation variables
  my $x  = 0;
  my $y  = 0;
  my $z  = 0;
  my $x2 = 0;
  my $y2 = 0;
  my $z2 = 0;
  my $xz = 0;
  my $yz = 0;
  my $n  = 0;
  #<<


  print "\"${seeds}\"\t";

  # set the seeds as ref of ref table
  my @seeds_table = split(/[,;]/,$seeds);
  my @seeds_2table = ();
  for (my $s =0 ; $s<@seeds_table; $s++) {
    my @seed = split(//,$seeds_table[$s]);
    push @seeds_2table, \@seed;
  }

  # foreach percentage 
  for (my $p = 0.2 ; $p <= 1.0; $p += 0.05) {
    my $alignment = &generate_alignment($length,int($p*$length));
    my @alignment_table = split(//, $alignment);
    # foreach trial
    for (my $trial = 0; $trial < 1000 ; $trial++) {
      my $multihit    = &multihit_seed_matches_alignment(\@alignment_table,\@seeds_2table);
      my $coveragehit = &coveragehit_seed_matches_alignment(\@alignment_table,\@seeds_2table);
      $x  += $multihit;
      $x2 += $multihit*$multihit;
      $y  += $coveragehit;
      $y2 += $coveragehit*$coveragehit;
      $z  += $p;
      $z2 += $p*$p;
      $xz += $multihit*$p;
      $yz += $coveragehit*$p;
      $n  += 1;
      shuffle(\@alignment_table);
    }
  }
  #>> correlation computation
  my $covxz  =      ($xz*$n) - ($x*$z);
  my $covyz  =      ($yz*$n) - ($y*$z);
  my $sigx   = sqrt(($x2*$n) - ($x*$x));
  my $sigy   = sqrt(($y2*$n) - ($y*$y));
  my $sigz   = sqrt(($z2*$n) - ($z*$z)); 
  my $corrxz = $covxz / ($sigx * $sigz);
  my $corryz = $covyz / ($sigy * $sigz);
  print "$corrxz\t$corryz\n";
  #<<
}
close S;
