#!/usr/bin/perl

# implement forward algorithm on a specific pair of sequences

$seq1 = "CGTC";
$seq2 = "TGTC";


#model

$tau = 0.3;
$delta = 0.1;
$epsilon = 0.5;

#housekeeping

if(length($seq2) > length($seq1)) {
$temp = $seq2;
$seq2 = $seq1;
$seq1 = $temp;
}


#begin state 

$seq1 = "0".$seq1;
$seq2 = "0".$seq2;


#initialization

$fm[0][0] = 1; 
$fx[0][0] = 0;
$fy[0][0] = 0;

for($x = 1; $x < length($seq1); $x++) {
$fx[$x][0] = $delta + $epsilon*$x;
$fy[$x][0] = 0;
$fm[$x][0] = 0;
}
for($y = 1; $y < length($seq2); $y++) {
$fx[0][$y] = 0;
$fy[0][$y] = $delta + $epsilon*$y;
$fm[0][$y] = 0;
}

for($x = 1; $x < length($seq1); $x++) {
	for($y = 1; $y < $x; $y++) {

		if(substr($seq1,$x,1) eq substr($seq2,$y,1)) {
		$match_emission = 0.1;
		}
		else {
		$match_emission = 0.05;
		}

	$emission_x = 0.25;
	$emission_y = 0.25;
	
	
$fm[$x][$y] = $match_emission * ((1-2*$delta-$tau)*$fm[$x-1][$y-1] + (1-$epsilon-$tau)*($fx[$x-1][$y-1]+$fy[$x-1][$y-1]));
print "res = ", (1-$epsilon-$tau)*($fx[$x-1][$y-1]+$fy[$x-1][$y-1]),"\n";
$fx[$x][$y] = $emission_x*($delta*$fm[$x-1][$y]+$epsilon*$fx[$x-1][$y]);
$fy[$x][$y] = $emission_y*($delta*$fm[$x][$y-1]+$epsilon*$fy[$x][$y-1]);
}

	for($q = 1; $q <= $x; $q++) {
		if(substr($seq1,$q,1) eq substr($seq2,$x,1)) {
		$match_emission = 0.1;
		}
		else {
		$match_emission = 0.05;
		}
	
	$emission_x = 0.25;
	$emission_y = 0.25;
	
	$fm[$q][$x] = $match_emission * ((1-2*$delta-$tau)*$fm[$q-1][$x-1] + (1-$epsilon-$tau)*($fx[$q-1][$x-1]+$fy[$q-1][$x-1]));
	$fx[$q][$x] = $emission_x*($delta*$fm[$q-1][$x]+$epsilon*$fx[$q-1][$x]);
	$fy[$q][$x] = $emission_y*($delta*$fm[$q][$x-1]+$epsilon*$fy[$q][$x-1]);
	}
}

print "$fm[length($seq1)-1][length($seq2)-1] $fx[length($seq1)-1][length($seq2)-1] $fy[length($seq1)-1][length($seq2)-1] \n";

print $tau*($fm[length($seq1)-1][length($seq2)-1] + $fx[length($seq1)-1][length($seq2)-1] + $fy[length($seq1)-1][length($seq2)-1]),"\n";