#!/usr/bin/perl -w
use strict;
#
# Programa para hacer un algoritmo MEME
#
# Definimos todas las variables
my $longitud_motivo;
my $nomfitxersequencies;
my $secuencia;
my @vectorsecuencia;
my @matrizsecuencias;
my %hashM;
my %matrizP;
my @matrizZ;
my @salidaM;
my $diferencia;
my $vers_previa = -99999999999999;
my $versemblanca = -9999999999;
my $vueltas;
my $i = 0;

if (scalar(@ARGV) < 2) {
    print STDERR "programaambfuncions.pl <longitudmotivo> <ficherosecuencias>\n";
    exit(1);
}
########## Leemos los input: el archivo con las secuencias y la longitud del motivo. ####################

$longitud_motivo = $ARGV[0];                                   # Argumento 1 para determinar la longitud 
                                                               # del motivo que queremos encontrar
$nomfitxersequencies = $ARGV[1];
if (!open(SECUENCIA,"< $nomfitxersequencies")) {               # Filehandle que nos permite recorrer
                                                               # las secuencias de 1 en 1
    print STDERR "programaambfuncions.pl: impossible obrir $nomfitxersequencies\n";
    exit(1);
}
while (<SECUENCIA>)                                     # aqui empezamos a leer las secuencias, fila fila
    
{
    chomp;
    
    $secuencia = $_;                                              # definimos esta variable como la linea
                                                                  # que el programa lee en cada momento
    
    @vectorsecuencia = split(//,uc($secuencia));     # creeamos un vector para cada secuencia del dataset
    
    $matrizsecuencias[$i] = [@vectorsecuencia];               # creamos una matriz con todos los vectores
    
    $i++;
}
close(SECUENCIA);

############################# CUERPO PRINCIPAL DEL PROGRAMA #############################################
######### Ahora tenemos que crear una matriz P0 aleatoria a partir de las secuencias introducidas #######

%matrizP = &random (\@vectorsecuencia, \$longitud_motivo, \@matrizsecuencias);
$diferencia = 10;
$vueltas = 0;
while ($diferencia > 0.0000001 && $vueltas < 10000)  # El programa compilara hasta que la verosimilitud
{                                                    # no aumente o de 1000 vueltas
    $vers_previa = $versemblanca;
    @matrizZ = &pasoE(\%matrizP, \@vectorsecuencia, \$longitud_motivo, \@matrizsecuencias);
    
    @salidaM = &pasoM(\@matrizZ, \@vectorsecuencia, \$longitud_motivo, \@matrizsecuencias);
    
    %matrizP = %{$salidaM[0]};
    
    %hashM = %{$salidaM[1]};
    
    $versemblanca = &calculovers(\%matrizP, \%hashM, \$longitud_motivo, \@vectorsecuencia);
    
    $diferencia = &diferencia(\$versemblanca, \$vers_previa);
    
    $vueltas++;
}
print "he hecho $vueltas vueltas\n";

########################################## OUTPUT #######################################################
#########################################################################################################

my @iniciosmotivo = &output(\@matrizZ, $longitud_motivo, \@vectorsecuencia);
$i = 0;
while ($i < scalar(@matrizZ))
{
    print "el inicio de motivo de la secuencia $i es $iniciosmotivo[$i]\n";
                                             # este print no es necesario porque lo que tenemos 
    $i++;                                    # como output son las secuencias con los motivos resaltados
}                                            # y la verosimilitud.

          # esto nos permite printar la secuencia con el background separado del motivo #

my $fila= 0;
while ($fila < scalar(@matrizsecuencias))                                 
{
    my $k = 0;
    
    while ($k < scalar(@vectorsecuencia))
	
    {
	if ($k < $iniciosmotivo[$fila])
	    
	{
	    print "$matrizsecuencias[$fila][$k]";
	}
	
	if ($k == $iniciosmotivo[$fila])
	{
	    print " ";
	}
	
	if ($k >= $iniciosmotivo[$fila] && $k < $iniciosmotivo[$fila]+ $longitud_motivo)
	{
	    print "$matrizsecuencias[$fila][$k]";
	}
	
	if ($k == $iniciosmotivo[$fila]+ $longitud_motivo -1)
	{
	    print " ";
	}
	
	if ($k>= $iniciosmotivo[$fila]+ $longitud_motivo)

	{
	    print "$matrizsecuencias[$fila][$k]";
	}
	
	$k++;
    }
    print " El valor de Z para este inicio de motivo es $matrizZ[$fila][$iniciosmotivo[$fila]]\n";
    print "\n";
    
    $fila = $fila + 1;
}
print "La verosimilitud de este motivo es de $versemblanca.\n";

 ######################
 ##                  ##
 ##    FUNCIONES     ##
 ##                  ##
 ######################
sub random {
    my $sumatotal =0;
    my @numeroaleatorio;
    my @vectorsecuencia = @{$_[0]};
    my $longitud_motivo = ${$_[1]};
    my @matrizsecuencias = @{$_[2]};
    my $i = 0;
    my $posicion;
    my $sumaA = 0;
    my $sumaC = 0;
    my $sumaG = 0;
    my $sumaT = 0;
    my @nucleotidos;
    my %matrizP;
    my $posicionmotivo = 0;
    
    while ($i < scalar(@matrizsecuencias))                              # leemos las secuencias una a una
	
    {
	
	$numeroaleatorio[$i] = int (rand (scalar(@vectorsecuencia)-$longitud_motivo));
	
	print "El numero aleatorio de la secuencia $i es $numeroaleatorio[$i].\n";
	
	$posicion =0;
	
	
	while ($posicion < scalar(@vectorsecuencia))                      # leemos cada nucleotidos 1 a 1
	    
	{
	    
	    if ($posicion < $numeroaleatorio[$i] || $posicion >= $numeroaleatorio[$i] + $longitud_motivo)
	    {
		if ($vectorsecuencia[$posicion] eq 'A')
		{
		    $sumaA = $sumaA + 1;                   # cada una de estas sumas es para contabilizar
		}                                         # cuantas letras de cada tipo hay en background
		
		if ($vectorsecuencia[$posicion] eq 'C')
		{
		    $sumaC = $sumaC + 1;
		}
		
		if ($vectorsecuencia[$posicion] eq 'G')
		{
		    $sumaG = $sumaG + 1;
		}
		if ($vectorsecuencia[$posicion] eq 'T')
		{
		    $sumaT = $sumaT + 1;
		}
	    }
	    
	    $posicion++;
	}
	
	$i++;
    }
    
    $sumatotal = $sumaA + $sumaG + $sumaG + $sumaT;   # esta variable la usamos para nomalizar la primera
                                                      # columna del hash (background)
    
    $nucleotidos[0] ='A';
    $nucleotidos[1] ='C';
    $nucleotidos[2] ='G';
    $nucleotidos[3] ='T';
                                                              # aqui llenamos la primera columna del hash
                                                                      # usando cada nucleotido como clave
    $matrizP{$nucleotidos[0]}[$longitud_motivo]= ($sumaA+1)/($sumatotal+4);
    $matrizP{$nucleotidos[1]}[$longitud_motivo]= ($sumaC+1)/($sumatotal+4); 
    $matrizP{$nucleotidos[2]}[$longitud_motivo]= ($sumaG+1)/($sumatotal+4);
    $matrizP{$nucleotidos[3]}[$longitud_motivo]= ($sumaT+1)/($sumatotal+4);

    $sumaA =0;
    $sumaC =0;
    $sumaG =0;
    $sumaT =0;
    
    $posicion = 0;

##### este es el paso en el que llenamos el resto de columnas del hash 
### #(tantas como posiciones tenga el motivo) ###########
    
    while ($posicionmotivo < $longitud_motivo) # como el hash se llena por columnas teniendo encuenta la
                                               # la posicion dentro del motivo, empezamos recorriendo
    {                                          # las posiciones
	$i=0;
	$sumatotal = 0;
	$sumaA = 0;
	$sumaC = 0;
	$sumaG = 0;
	$sumaT = 0;
	while ($i < scalar(@matrizsecuencias))   # aqui recorremos todas las filas con la posicion fijada

	{
	    
	    if ($matrizsecuencias[$i][$numeroaleatorio[$i]+$posicionmotivo] eq 'A')
	    {
		$sumaA = $sumaA + 1;
	    }
	    
	    if ($matrizsecuencias[$i][$numeroaleatorio[$i]+$posicionmotivo] eq 'C')
	    {
		$sumaC = $sumaC + 1;
	    }
	    
	    if ($matrizsecuencias[$i][$numeroaleatorio[$i]+$posicionmotivo] eq 'G')
	    {
		$sumaG = $sumaG + 1;
	    }
	    if ($matrizsecuencias[$i][$numeroaleatorio[$i]+$posicionmotivo] eq 'T')
	    {
		$sumaT = $sumaT + 1;
	    }
	    
	    $i++;
	    
	}
	
	$sumatotal = $sumaA + $sumaC + $sumaG + $sumaT;
	
	$matrizP{$nucleotidos[0]}[$posicionmotivo]= ($sumaA+1)/($sumatotal+4);
	$matrizP{$nucleotidos[1]}[$posicionmotivo]= ($sumaC+1)/($sumatotal+4);
	$matrizP{$nucleotidos[2]}[$posicionmotivo]= ($sumaG+1)/($sumatotal+4);
	$matrizP{$nucleotidos[3]}[$posicionmotivo]= ($sumaT+1)/($sumatotal+4);
	
$posicionmotivo++;
    }
    
    return %matrizP;
}
sub pasoE {
    my $i = 0;
    my $contadorposicion = 0;
    my $multiplica;
    my $contadorback;
    my @matrizsecuencias = @{$_[3]};
    my %matrizP = %{$_[0]};
    my @vectorsecuencia = @{$_[1]};
    my $longitud_motivo = ${$_[2]};
    my $k =0;
    my $contadortrasero;
    my $fila = 0;
    my $logsuma;
    my $exp = 0;
    my @matrizZ;
    my $columna;
    while ($i < scalar(@matrizsecuencias))
    {
	$contadorposicion=0;                  # este marcador sirve para identificar que posicion estamos
                                              # analizando dentro de la secuencia y se corresponde con el
                                              # inicio del motivo
	
	@vectorsecuencia = @{$matrizsecuencias[$i]};          # definimos esta variable como la linea que
                                                              # el programa lee en cada momento
	
	while ($contadorposicion < scalar(@vectorsecuencia) - $longitud_motivo +1) 

                                                                  # este bucle nos permite ir recorriendo
                                                # el vector que contiene la secuencia posicion a posicion
	    
	{
	    $multiplica = 0;           # en este numero almacenamos cada multiplicacion de posibilidades
                                       # (una casilla de la matriz Zij
	    
	    
	    $contadorback = 0;                  # variable para contabilizar las posiciones de background
                                                # anteriores al motivo
	    
	    while ($contadorback < $contadorposicion)            # esta iteracion nos permite calcular el
                                                                 # valor del background que se encuentra
                                             # anterior a la posicion que ocupa el motivo en cada momento
		
	    {
		
		$multiplica = ($multiplica) + log($matrizP{$vectorsecuencia[$contadorback]}[$longitud_motivo]);
		
		$contadorback++;
	    }
	    
	    $k = 0;
	    
	    while ($k < $longitud_motivo)          # con este bucle recorremos unicamente las posiciones
                                                   # que corresponden al motivo
		
	    {
		
		$multiplica = ($multiplica) + log($matrizP{$vectorsecuencia[$contadorposicion+$k]}[$k]);
		
		
		$k++;
	    }
	    
	    $contadortrasero = $contadorposicion + $longitud_motivo; # variable que recorre las posicione
                                # de background posteriores a la localización del motivo en cada momento
	    
	    $k = $contadortrasero;
	    
	    while ($contadortrasero < scalar(@vectorsecuencia))      # bucle que contabiliza el background
                                # existente detras de la posición del motivo
		
	    {
		
		$multiplica = $multiplica + log($matrizP{$vectorsecuencia[$contadortrasero]}[$longitud_motivo]);
		
		$contadortrasero++;
		
	    }
	    
	    $matrizZ[$i][$contadorposicion] = $multiplica; # vamos llenando cada posicion de la matriz Z
                                                           # con el valor almacenado en $multiplica
	    
               	    $contadorposicion++;
	    
	}
	
	$i++;
	
    }
    $fila = 0;                             # esta variable se usa como marcador de las filas de la matriz

    $columna = 0;
    
    
  ################### Hasta aqui la matriz Z no normalizada #############################################

 ################### Sumatorio para usar como denominador en la normalizacion ########################
    
    $exp = 0;                    # variable para acumular los exponentes de los valores de cada secuencia

    $fila = 0;

    while ($fila < scalar(@matrizZ))            # iteracion para recorrer la matriz Z posicion a posicion
	
    {
	$columna = 0;

	$exp = 0;
	
	while ($columna < scalar(@{$matrizZ[$fila]}))
	    
	{
	    
	    $exp = $exp + exp($matrizZ[$fila][$columna]);    # aqui sumamos los exponentes de los valores
                                                             # de cada posicion de una secuencia
	    $columna++;
	    
	}

	$logsuma = log($exp);                        # aqui hacemos el logaritmo de la suma de exponentes
	
	$columna = 0;

	while ($columna < scalar(@{$matrizZ[$fila]}))
	    
	{
	    
	    $matrizZ[$fila][$columna] = exp($matrizZ[$fila][$columna] - $logsuma);
                                                                    # este es el calculo DE NORMALIZACION
	    
	    $columna++;
	}
	
	
	$fila++;
    }

    return @matrizZ;
}

sub pasoM {                                        # esta funcion corresponde con el paso M del algoritmo
                                     # aqui calculamos y normalizamos la matriz P a partir de la matriz Z
    my @salidaM;
    my $sumaA;
    my $sumaC;
    my $sumaG;
    my $sumaT;
    my $posicionmotivo = 0;
    my %hashM;
    my @nucleotidos;
    my $i;
    my $secuencia;
    my $fila;
    my @matrizZ = @{$_[0]};
    my @vectorsecuencia = @{$_[1]};
    my $longitud_motivo = ${$_[2]};
    my @matrizsecuencias = @{$_[3]};
    my $posicion;
    my $sumaposicion = 0;
    my @valortotal;
    my $totalA;
    my $totalC;
    my $totalG;
    my $totalT;
    my $backA;
    my $backC;
    my $backG;
    my $backT;
    my $backnormalizador;
    
    $nucleotidos[0] = 'A';
    $nucleotidos[1] = 'C';
    $nucleotidos[2] = 'G';
    $nucleotidos[3] = 'T';
    
    
    while ($posicionmotivo < $longitud_motivo) # este paso calcula el hash del motivo base es decir suma
                                               # todas las zetas de una letra y una posicion del motivo
    {
	
	$sumaA = 0;
	$sumaC = 0;
	$sumaG = 0;
	$sumaT = 0;
	
	$fila = 0;
	
	while ($fila < scalar(@matrizsecuencias))
	    
	{
	    
	    $secuencia = $matrizsecuencias[$fila];
	    
	    $posicion = 0;
	    
	    while ($posicion <= scalar(@vectorsecuencia) - $longitud_motivo)
		
	    {
		
		if ($vectorsecuencia[$posicion + $posicionmotivo] eq 'A')
		    
		{
		    
		    $sumaA = $sumaA + $matrizZ[$fila][$posicion];
		    
		}
		
		if ($vectorsecuencia[$posicion + $posicionmotivo] eq 'C')
		    
		{
		    
		    $sumaC = $sumaC + $matrizZ[$fila][$posicion];
		    
		}
		
		if ($vectorsecuencia[$posicion + $posicionmotivo] eq 'G')
		    
		{
		    
		    $sumaG = $sumaG + $matrizZ[$fila][$posicion];
		    
		}
		
		if ($vectorsecuencia[$posicion + $posicionmotivo] eq 'T')
		    
		{
		    
		    $sumaT = $sumaT + $matrizZ[$fila][$posicion];
		    
		}
		
		$posicion++;
		
	    }
	    
	    
	    $fila++;    # este contador no esta en un while porque usamos el filehandle para realizarlo,
                        # no obstante necesitamos que el marcador vaya subiendo para ir llenando el hash
                        # base. (que consta solo de las 'Z' sumadas).
	}
	
	
	
	
	$hashM{$nucleotidos[0]}[$posicionmotivo] = $sumaA; # aqui definimos cada posicion del hash, los 4
	$hashM{$nucleotidos[1]}[$posicionmotivo] = $sumaC; # nucleotidos como clave y tantas columnas
	$hashM{$nucleotidos[2]}[$posicionmotivo] = $sumaG; # como posiciones tenga el motivo
	$hashM{$nucleotidos[3]}[$posicionmotivo] = $sumaT;
	
	
	$posicionmotivo++;
    }
    
    $posicion=0;
    
    while ($posicion < $longitud_motivo)     # en este paso sumamos todos los valores de una columna para
                                             # despues poder normalizarlas
    {
	
	$fila= 0;                            # recorremos una fila (nucleotido)
	
	$sumaposicion = 0;                   # igualamos a cero para reiniciar la variable que acumula la
                                             # suma de cada columna
	
	while ($fila < scalar(@nucleotidos))
	    
	{
	    
	    $sumaposicion = $sumaposicion + $hashM{$nucleotidos[$fila]}[$posicion];
	    
                                     # aqui vamos sumando en la variable todos los valores de una columna
	    $fila++;
	    
	}
	
	$valortotal[$posicion] = $sumaposicion;
                                     # aqui pasamos el valor de la suma de una columna a un vector
	
	$posicion++;
	
    }

# Aqui hemos creado un vector con la suma de cada columna del hash BASE. a continuacion debemos
# normalizar el hash base cogiendolo y dividiendo cada posicion de una columna por la suma
# de los valores de esta

    $posicion = 0;
    
    while ($posicion < $longitud_motivo)
	
    {
	$fila = 0;
	
	while ($fila < scalar(@nucleotidos))
	    
	{
	    
	    $matrizP{$nucleotidos[$fila]}[$posicion] = ($hashM{$nucleotidos[$fila]}[$posicion]+1)/($valortotal[$posicion]+4);
	    
# este es el paso de normalizacion del motivo ya que dividimos cada posicion por la suma de la
# de la columna a la que pertenece

	    $fila++;
}
	
	$posicion++;
    }

# A partir de aqui generamos la posicion de background. Lo primero que hacemos es sumar todas las
# probabilidades de una letra en todas las posiciones de todas las secuencias

    $sumaA = 0;
    $sumaC = 0;
    $sumaG = 0;
    $sumaT = 0;
    $posicionmotivo = 0;
    
    $nucleotidos[0] = 'A';
    $nucleotidos[1] = 'C';
    $nucleotidos[2] = 'G';
    $nucleotidos[3] = 'T';

    while ($posicionmotivo < $longitud_motivo)
                                           # como el hash va por posiciones la variable principal es esta
                                           # ya que necesitamos sumar todas la probabilidades de cada
                                           # nucleotido en cada posicion y secuencia a secuencia
    {
	
	$i = 0;

	while ($i < scalar(@matrizsecuencias))
	    
	{                                       # este paso es un poco liado ya que tenemos que sumar las
                                                # probabilidades de cada nt en cada posicion, para saber
                                                # que proporcion de estos es parte del motivo
	    $posicion = 0;
	    
	    while ($posicion <= scalar(@vectorsecuencia) - $longitud_motivo)
		
	    {
		
		if ($vectorsecuencia[$posicion + $posicionmotivo] eq 'A')
		    
		{
		    
		    $sumaA = $sumaA + $matrizZ[$i][$posicion];

		}
		
		if ($vectorsecuencia[$posicion + $posicionmotivo] eq 'C')
		    
		{
		    
		    $sumaC = $sumaC + $matrizZ[$i][$posicion];
		    
		}
		
		if ($vectorsecuencia[$posicion + $posicionmotivo] eq 'G')
		    
		{
		    
		    $sumaG = $sumaG + $matrizZ[$i][$posicion];

		}
		
		if ($vectorsecuencia[$posicion + $posicionmotivo] eq 'T')
    
		{
		    
		    $sumaT = $sumaT + $matrizZ[$i][$posicion];
		    
		}
		
		$posicion++;
		
	    }                                 # Al acabar este paso obtendremos las posibilidades de cada
                                              # nucleotido en TODO el motivo. A continuacion debemos
	    $i++;                             # contar cuantos nucleotidos de cada tipo hay entre todas
	}                                     # las secuencias y para acabar, restar uno del otro para
                                              # saber cuantos de los nts totales son del background
	
	$posicionmotivo++;
    }
    $totalA = 0;
    $totalC = 0;
    $totalG = 0;
    $totalT = 0;
    
    $i = 0;
    while ($i < scalar(@matrizsecuencias))      # este es el paso que sirve para contar el numero de nts
                                                # de cada tipo en TODAS las secuencias juntas
    {
	$secuencia = $matrizsecuencias[$i];
	
	$posicion = 0;
	
	while ($posicion < scalar(@vectorsecuencia))
	    
	{
	    if ($vectorsecuencia[$posicion] eq 'A')
                                                  # comprobamos que letra hay en cada posicion del vector
                                                  # y las vamos contando
	    {
		
		$totalA = $totalA + 1;
		
	    }
	    
	    if ($vectorsecuencia[$posicion] eq 'C')
		
	    {
		
		$totalC = $totalC + 1;
		
	    }
	    
	    if ($vectorsecuencia[$posicion] eq 'G')
		
	    {
		
		$totalG = $totalG + 1;
		
	    }
	    
	    if ($vectorsecuencia[$posicion] eq 'T')
		
	    {
		
		$totalT = $totalT + 1;
		
	    }
	    
	    $posicion++;
	    
	}
	
	$i++;
    }
    $backA = ($totalA - $sumaA); # columna de background SIN normalizar. La obtenemos restando al numero
    $backC = ($totalC - $sumaC); # de nts total la proporcion de nts que son parte del motivo, de esta
    $backG = ($totalG - $sumaG); # manera nos queda la parte de los nts que son parte del background
    $backT = ($totalT - $sumaT);

    $hashM{$nucleotidos[0]}[$longitud_motivo] = $backA;       # aqui llenamos la posicion del background 
                                                              # del hashM, que es la matriz P
    $hashM{$nucleotidos[1]}[$longitud_motivo] = $backC;       # sin normalizar, lo necesitamos para 
                                                              # calcular la razon de
    $hashM{$nucleotidos[2]}[$longitud_motivo] = $backG;       # verosimilitud

    $hashM{$nucleotidos[3]}[$longitud_motivo] = $backT;

# ahora debemos normalizar el background obtenido, de la misma forma que hemos hecho con las columnas
# correspondientes a cada posicion del motivo

    $backnormalizador = $backA + $backC + $backG + $backT;   # variable para normalizar la columna
                                                             # de background, que corresponde con la suma
                                                             # de todos los valores de la columna
                                                             # sumamos cuatro para compensar que hemos
                                                             # SUMADO 1 las sumas de cada nucleotido
    $matrizP{$nucleotidos[0]}[$longitud_motivo]= ($backA+1)/($backnormalizador+4);
                                                             
    $matrizP{$nucleotidos[1]}[$longitud_motivo]= ($backC+1)/($backnormalizador+4);
                                                             
    $matrizP{$nucleotidos[2]}[$longitud_motivo]= ($backG+1)/($backnormalizador+4);

    $matrizP{$nucleotidos[3]}[$longitud_motivo]= ($backT+1)/($backnormalizador+4);

    @salidaM = (\%matrizP, \%hashM);                    # con esto conseguimos que utilice la matriz P
                                                        # calculada a partir de la Z, en vez de la random
    return @salidaM;
}

sub calculovers {                             # esta funcion sirve para calcular la verosimilitud de cada
                                              # ciclo del algoritmo, y es lo que nos permite comprobar si
    my $i = 0;                                # el motivo que encuentra nuestro programa es robusto o no
    my %matrizP = %{$_[0]};
    my %hashM = %{$_[1]};
    my $longitud_motivo = ${$_[2]};
    my @vectorsecuencia = @{$_[3]};
    my $sumaA = 0;
    my $sumaC = 0;
    my $sumaG = 0;
    my $sumaT = 0;
    my $sumabackground = 0;
    my @nucleotidos;
    while ($i <= $longitud_motivo)
    {
	$nucleotidos[0] = 'A';
	$nucleotidos[1] = 'C';
	$nucleotidos[2] = 'G';
	$nucleotidos[3] = 'T';
	$sumaA = $sumaA + ($hashM{$nucleotidos[0]}[$i] * log($matrizP{$nucleotidos[0]}[$i]));
	$sumaC = $sumaC + ($hashM{$nucleotidos[1]}[$i] * log($matrizP{$nucleotidos[1]}[$i]));
	$sumaG = $sumaG + ($hashM{$nucleotidos[2]}[$i] * log($matrizP{$nucleotidos[2]}[$i]));
	$sumaT = $sumaT + ($hashM{$nucleotidos[3]}[$i] * log($matrizP{$nucleotidos[3]}[$i]));
	$i++;
    }
    $sumabackground = $sumabackground + ($hashM{$nucleotidos[0]}[$longitud_motivo] * log($matrizP{$nucleotidos[0]}[$longitud_motivo])*(scalar(@vectorsecuencia) - $longitud_motivo));
    $sumabackground = $sumabackground + ($hashM{$nucleotidos[1]}[$longitud_motivo] * log($matrizP{$nucleotidos[1]}[$longitud_motivo])*(scalar(@vectorsecuencia) - $longitud_motivo));
    $sumabackground = $sumabackground + ($hashM{$nucleotidos[2]}[$longitud_motivo] * log($matrizP{$nucleotidos[2]}[$longitud_motivo])*(scalar(@vectorsecuencia) - $longitud_motivo));
    $sumabackground = $sumabackground + ($hashM{$nucleotidos[3]}[$longitud_motivo] * log($matrizP{$nucleotidos[3]}[$longitud_motivo])*(scalar(@vectorsecuencia) - $longitud_motivo));
    $versemblanca = $sumaA + $sumaC + $sumaG + $sumaT + $sumabackground;
    return $versemblanca;
}

sub diferencia {                                  # esta funcion sirve para calcular la diferencia entre
                                                  # la verosimilitud de cada ciclo del algoritmo. Cuando
    my $versemblanca = ${$_[0]};                  # se estanca y deja de crecer, el algoritmo para porque
    my $vers_previa = ${$_[1]};                   # ya ha encontrado la mejor matriz que define al motivo
    my $diferencia = 0;
    if ($versemblanca < $vers_previa)
	
    {
	$diferencia = $vers_previa - $versemblanca;
    }
    
    else
	
    {
	$diferencia = $versemblanca - $vers_previa;
    }
    return $diferencia;
}

sub output {                                  # esta funcion sirve para formar un vector que nos almacene
                                              # las posiciones de cada secuencia en las que empieza el
    my $imax = 0;                             # motivo que hemos encontrado, definido por la matriz P
    my $max = 0;
    my $fila = 0;
    my @iniciosmotivo;
    my @matrizZ = @{$_[0]};
    my $longitud_motivo = $_[1];
    my @vectorscuencia = @{$_[2]};
    while ($fila < scalar(@matrizZ))
	
    {
	$imax = 0;
	$i = 0;
	$max = 0;
	
	while ($i < scalar(@vectorsecuencia)-$longitud_motivo+1)
	    
	{
	    if ($matrizZ[$fila][$i] > $max)
		
	    {
		$max = $matrizZ[$fila][$i];
		$imax = $i;
	    }
	    
	    $i++;
	}
	$iniciosmotivo[$fila]= $imax;
	
	$fila++;
    }
    return @iniciosmotivo;
}

# FIN del PROGRAMA(ciao) #
# implementado por : Irene Andretta
                   # Pablo Martinez
                   # Diana Reyes