#!/usr/bin/perl -w

use strict;

open(FIT2,"< $ARGV[1]");
my $i=0;
my @v;
my $l;
my $j;
my $linea =<FIT2>;
$linea=<FIT2>;
while (<FIT2>){
     chomp($linea);
     my @tmp = split(//,$linea);
     $l=scalar(@tmp);
     $j=0;
     my $f=scalar(@v);
     while ($j<$l){
	  $v[$f+$j]=$tmp[$j];
	  $j=$j+1;
     }
     $linea=<FIT2>;
}
my @sec=@v;

my %ocu=("A"=>[0],
	 "C"=>[0],
	 "G"=>[0],
	 "T"=>[0]);

my @k=keys(%ocu);

$i=0;
my $N=scalar(@v);
my $c;
while ($i<$N){
     $c=0;
     while ($c<scalar(@k)){
	  if ($v[$i] eq $k[$c]){
		$ocu{$k[$c]}[0]=$ocu{$k[$c]}[0]+1;
	  }
	  $c=$c+1;
     }
     $i=$i+1;
 }

close (FIT2);



open(FIT,"< $ARGV[0]");

$j=1;

while ($j<=13){
    @v=@sec;
    $linea=<FIT>;
    my $nombre;
    if ($linea=~m/FA/){
	my @lll = split/\s+/, $linea;
	$nombre = "$lll[1] $lll[2]";
	print "\n$nombre\n";
	next;
    }
    elsif ($linea =~m/\A[\D]+/){
	next;
    }

    #matrix
    my %h = ("A" => [0, 0, 0, 0, 0, 0, 0, 0, 0],
	     "C" => [0, 0, 0, 0, 0, 0, 0, 0, 0],
	     "G" => [0, 0, 0, 0, 0, 0, 0, 0, 0],
	     "T" => [0, 0, 0, 0, 0, 0, 0, 0, 0]);

    my $lm=0;
    while ($linea =~ /[\d+][\s]+(\d+)[\s]+(\d+)[\s]+(\d+)[\s]+(\d+)/){
	chomp($linea);

	$h{A}[$lm]=$1;
	$h{C}[$lm]=$2;
	$h{G}[$lm]=$3;
	$h{T}[$lm]=$4;
	
	$lm=$lm+1;
	$linea=<FIT>;
    }
    
    #print de la matrix
    @k = keys(%h);
    $c = 0;
    #print "martix $j\n";
    while ($c < scalar(@k)) {
	my $d = 0;
	print "$k[$c]:";
	while ($d < $lm) {
	    print "\t$h{$k[$c]}[$d]";
	    $d = $d + 1;
	}
	print "\n";
	$c = $c + 1;
    }

    #PWM
    my $p=0;
    my %sumnuc = ("A" => [0],
		  "C" => [0],
		  "G" => [0],
		  "T" => [0]);
    $c=0;
    while ($c < scalar(@k)){
	$p = 0;
	while ($p < $lm){
	    $sumnuc{$k[$c]}[0]=$sumnuc{$k[$c]}[0]+$h{$k[$c]}[$p];
	    $p=$p+1;
	}
	$c=$c+1;
    }

    $p=0;
    while ($p < $lm){
	$c=0;
	while ($c<scalar(@k)){
	    my $primo=($h{$k[$c]}[$p])/($sumnuc{$k[$c]}[0]);
	    my $secundo;
	    if ($primo != 0){
		$secundo=($ocu{$k[$c]}[0])/$N;
		$h{$k[$c]}[$p]=log($primo)-log($secundo);
	    }
	    else {
		$h{$k[$c]}[$p]=-999;
            }
	
	    $c=$c+1;
	}
	$p=$p+1;
    }


    #print de la matrix
    #@k = keys(%h);
    #$c = 0;
    #print "martix $j\n";
    #while ($c < scalar(@k)) {
	#my $d = 0;
	#print "$k[$c]:";
	#while ($d < $lm) {
	    #print "\t$h{$k[$c]}[$d]";
	    #$d = $d + 1;
	#}
	#print "\n";
	#$c = $c + 1;
    #}

    #PUNTO 3
    my $score=0;
    my $w=0;
    my $posin=0;
    while ($w<$lm){
	$score=$score+$h{$v[$w]}[$w];
	$w=$w+1;
    }
    my $y=1;
    while ($y<=$N-$lm){
	my $s=0;
	$w=0;
	while ($w<$lm){
	    $s=$s+$h{$v[$y+$w]}[$w];
	    $w=$w+1;
	}
	if ($s>$score){
	    $score=$s;
	    $posin=$y;
	}
	$y=$y+1;
    }
    print "\n el score mas alto es $score\n";
    print "\n se encuentra en la posicion $posin\n";
    
    #PUNTO 4
    my $n=1;
    my $count=0;
    while ($n<=100){
	$N=scalar(@v);
	$i = $N - 1;
	while ($i >= 0) {
	    my $r = int(rand($i+1));
	    if ($i != $r) {
		my $tmp = $v[$i];

		$v[$i] = $v[$r];
		$v[$r] = $tmp;
	    }
	    $i = $i - 1;
	}

	$y=0;
	my $scoretmp=-10000;
	my $posintmp=0;

	while ($y<=$N-$lm){
	    my $s=0;
	    my $w=0;
	    while ($w<$lm){
		$s=$s+$h{$v[$y+$w]}[$w];
		$w=$w+1;
	    }
	    if ($s>$scoretmp){
		$scoretmp=$s;
		$posintmp=$y;
	    }
	    $y=$y+1;
	}
	if ($scoretmp>=$score){
	    $count=$count+1;
	}
	$n=$n+1;
    }
    my $pvalue=$count/100;
    print " \npvalue $pvalue\n\n";

    $j=$j+1;
}