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

# shuffle an array
sub shuffle {
    my $array = shift;
    my $i = @$array;
    while ( --$i ) {
        my $j = int rand( $i+1 );
        @$array[$i,$j] = @$array[$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,$seeds_string) = @_;
    my $result = 0;
    my @seeds = split(/[,;]/,$seeds_string);
    for (my $s=0; $s <= $#seeds;$s++) {
        my @seed = split(//,$seeds[$s]);
         for (my $i = 0; $i < (length $alignment) - $#seed; $i++) {
             for (my $j = 0; $j <= $#seed ; $j++) {
                 if ((substr($alignment,$i+$j,1) eq "0") && ($seed[$j] eq "#")) {
                     goto no_match;
                 }
             }
             $result++;
           no_match:
         }
    }
    return $result;
}

# compute the multihit for one/several seeds
sub coveragehit_seed_matches_alignment($$){
    my ($alignment,$seeds_string) = @_;
    my @seeds = split(/[,;]/,$seeds_string);
    my @coverage = ();
    @coverage = (@coverage, ('0') x (length $alignment));

    for (my $s=0; $s <= $#seeds;$s++) {
        my @seed = split(//,$seeds[$s]);
         for (my $i = 0; $i < (length $alignment) - $#seed; $i++) {
             for (my $j = 0; $j <= $#seed ; $j++) {
                 if ((substr($alignment,$i+$j,1) eq "0") && ($seed[$j] eq "#")) {
                     goto no_match;
                 }
             }
             for (my $j = 0; $j <= $#seed ; $j++) {
                 if ($seed[$j] eq "#") {
                     $coverage[$i+$j] = 1;
                 }
             }
           no_match:
         }
    }
    #print join(",", @coverage);
    my $result = eval join('+', @coverage);
    return $result;
}


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

# foreach seed
open S,$seedfile;
while (<S>) {
    chomp $_;
    if ($_ eq ""){
        next;
    }
    my $s = "$_";
    # foreach percentage 
    for (my $p = 0.2 ; $p <= 1.0; $p += 0.05) {
        for (my $trial = 0; $trial < 1000 ; $trial++) {
            my $a = &generate_alignment($l,int($p*$l));
            print "${l}\t${s}\t${p}\t";
            print &multihit_seed_matches_alignment($a,$s);
            print "\t";
            print &coveragehit_seed_matches_alignment($a,$s);
            print "\n";
        }
    }
}
close S;