#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdarg.h>
#include <string.h>

#include "io_lib_header.h"
#include "util_lib_header.h"
#include "define_header.h"
#include "dp_lib_header.h"

#define hmm_add(x,y) ((x==UNDEFINED || y==UNDEFINED)?UNDEFINED:(x+y))
#define MAX_EMISSION 256

/*********************************************************************************/
/*                                                                               */
/*                                                                               */
/*                       Procons dp                                              */
/*                                                                               */
/*                                                                               */
/*********************************************************************************/
char alphabetDefault[] = "ARNDCQEGHILKMFPSTWYV";
double emitPairsDefault[20][20] = {
  {0.02373072f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  {0.00244502f, 0.01775118f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  {0.00210228f, 0.00207782f, 0.01281864f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  {0.00223549f, 0.00161657f, 0.00353540f, 0.01911178f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  {0.00145515f, 0.00044701f, 0.00042479f, 0.00036798f, 0.01013470f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  {0.00219102f, 0.00253532f, 0.00158223f, 0.00176784f, 0.00032102f, 0.00756604f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  {0.00332218f, 0.00268865f, 0.00224738f, 0.00496800f, 0.00037956f, 0.00345128f, 0.01676565f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  {0.00597898f, 0.00194865f, 0.00288882f, 0.00235249f, 0.00071206f, 0.00142432f, 0.00214860f, 0.04062876f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  {0.00114353f, 0.00132105f, 0.00141205f, 0.00097077f, 0.00026421f, 0.00113901f, 0.00131767f, 0.00103704f, 0.00867996f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  {0.00318853f, 0.00138145f, 0.00104273f, 0.00105355f, 0.00094040f, 0.00100883f, 0.00124207f, 0.00142520f, 0.00059716f, 0.01778263f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  {0.00449576f, 0.00246811f, 0.00160275f, 0.00161966f, 0.00138494f, 0.00180553f, 0.00222063f, 0.00212853f, 0.00111754f, 0.01071834f, 0.03583921f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  {0.00331693f, 0.00595650f, 0.00257310f, 0.00252518f, 0.00046951f, 0.00312308f, 0.00428420f, 0.00259311f, 0.00121376f, 0.00157852f, 0.00259626f, 0.01612228f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  {0.00148878f, 0.00076734f, 0.00063401f, 0.00047808f, 0.00037421f, 0.00075546f, 0.00076105f, 0.00066504f, 0.00042237f, 0.00224097f, 0.00461939f, 0.00096120f, 0.00409522f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  {0.00165004f, 0.00090768f, 0.00084658f, 0.00069041f, 0.00052274f, 0.00059248f, 0.00078814f, 0.00115204f, 0.00072545f, 0.00279948f, 0.00533369f, 0.00087222f, 0.00116111f, 0.01661038f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  {0.00230618f, 0.00106268f, 0.00100282f, 0.00125381f, 0.00034766f, 0.00090111f, 0.00151550f, 0.00155601f, 0.00049078f, 0.00103767f, 0.00157310f, 0.00154836f, 0.00046718f, 0.00060701f, 0.01846071f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
  {0.00631752f, 0.00224540f, 0.00301397f, 0.00285226f, 0.00094867f, 0.00191155f, 0.00293898f, 0.00381962f, 0.00116422f, 0.00173565f, 0.00250962f, 0.00312633f, 0.00087787f, 0.00119036f, 0.00180037f, 0.01346609f, 0.0f, 0.0f, 0.0f, 0.0f},
  {0.00389995f, 0.00186053f, 0.00220144f, 0.00180488f, 0.00073798f, 0.00154526f, 0.00216760f, 0.00214841f, 0.00077747f, 0.00248968f, 0.00302273f, 0.00250862f, 0.00093371f, 0.00107595f, 0.00147982f, 0.00487295f, 0.01299436f, 0.0f, 0.0f, 0.0f},
  {0.00039119f, 0.00029139f, 0.00021006f, 0.00016015f, 0.00010666f, 0.00020592f, 0.00023815f, 0.00038786f, 0.00019097f, 0.00039549f, 0.00076736f, 0.00028448f, 0.00016253f, 0.00085751f, 0.00015674f, 0.00026525f, 0.00024961f, 0.00563625f, 0.0f, 0.0f},
  {0.00131840f, 0.00099430f, 0.00074960f, 0.00066005f, 0.00036626f, 0.00070192f, 0.00092548f, 0.00089301f, 0.00131038f, 0.00127857f, 0.00219713f, 0.00100817f, 0.00054105f, 0.00368739f, 0.00047608f, 0.00102648f, 0.00094759f, 0.00069226f, 0.00999315f, 0.0f},
  {0.00533241f, 0.00169359f, 0.00136609f, 0.00127915f, 0.00119152f, 0.00132844f, 0.00178697f, 0.00194579f, 0.00071553f, 0.01117956f, 0.00914460f, 0.00210897f, 0.00197461f, 0.00256159f, 0.00135781f, 0.00241601f, 0.00343452f, 0.00038538f, 0.00148001f, 0.02075171f}
};

static void DisplayMatState ( MatState *S, char *s);






void check_viterbiL ( Alignment *A,int *ns, int **ls, Constraint_list *CL);
MatState *viterbi2path2 ( double ***Sc, int ***St, Hmm *H, MatState *S, MatState *E);
void testfunc ( MatState *S, char *s);

#ifdef IN_PGROGRESS
/*********************************************************************************/
/*                                                                               */
/*                                                                               */
/*                     MSA Analyzer                                              */
/*                                                                               */
/*                                                                               */
/*********************************************************************************/
Alignment * analyze_alignment ( Alignment *A)
{
  evaluate_alignment (A);
  H=define_msa_model (-100);
  M=seq_viterbi_hmm (A->seq_al[0], H);
  path=seq_viterbi2path ( seq, H, M);
}
    

Hmm* define_msa_model(double penalty)
{
  Hmm *H;
  double freeT=0;
  int n=0;
  HmmState *S;
  
  
  H=declare_hmm(2);
  H->freeT=freeT=0;
 
  H->forbiden=FORBIDEN;
  H->start=START_STATE;
  H->end=END_STATE;
 
  /*define START*/
  S=H->S[n];  
  sprintf (S->name, "START"); S->state=n;
  
  S->DI=0;
  S->DJ=0;
  S->em=freeT;
 
  sprintf ( (S->T[S->nT])->name, "C") ;(S->T[S->nT])->tr=freeT  ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "W");(S->T[S->nT])->tr=freeT  ;S->nT++;
  n++;
  /*define END*/
  S=H->S[n];  
  sprintf (S->name, "END"); S->state=n;
  S->DI=0;
  S->DJ=0;
  S->em=freeT;
  n++;

  /*define Correct*/
  S=H->S[n];  
  sprintf (S->name, "C"); S->state=n;
  S->DI=1;
  S->DJ=0;
  S->em=H->forbiden;
  S->em_func=em_correct_msa;
  
  sprintf ( (S->T[S->nT])->name, "C") ;(S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "W");(S->T[S->nT])->tr=penalty  ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "END");   (S->T[S->nT])->tr=freeT;S->nT++;
  n++;
  
  /*define Wrong*/
  S=H->S[n];  
  sprintf (S->name, "INSERT"); S->state=n;
  S->DI=1;
  S->DJ=0;
  S->em=H->forbiden;
  S->em_func=em_wrong_msa;
  sprintf ( (S->T[S->nT])->name, "C") ; (S->T[S->nT])->tr=penalty;S->nT++;
  sprintf ( (S->T[S->nT])->name, "W"); (S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "END");    (S->T[S->nT])->tr=-gop;S->nT++;
  n++;
  
  /*define LInsert*/
  S=H->S[n];  
  sprintf (S->name, "LINSERT"); S->state=n;
  S->DI=1;
  S->DJ=0;
  S->em=lgep;
  
  sprintf ( (S->T[S->nT])->name, "INSERT")  ;(S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "LINSERT");(S->T[S->nT])->tr=freeT;S->nT++;
  n++;
  
  H=bound_hmm ( H);
  return H;
}
#endif
/*********************************************************************************/
/*                                                                               */
/*                                                                               */
/*                     simple HMM: Viterbi                                       */
/*                                                                               */
/*                                                                               */
/*********************************************************************************/
double pavie_em_func (Hmm*H, HmmState *S, int v);
Hmm* define_full_model(int nstate,char **state_list, char *model_name,Generic_em_func evaluation_func );
char **produce_state_name (int nstate,char **list, char *model_name, Hmm* H);
double** seq_viterbi_hmm (char *seq, Hmm *H);
int * seq_viterbi2path (char *s, Hmm *H, double **M);
double analyze_sequence ( char *seq, Hmm*H);

double pavie_emission (Hmm*H, HmmState *S, int v)
{
  char *n;


  n=S->name;
  
  if ( v==n[0] || ( v=='*' && n[0]=='E')) return H->freeT;
  return H->forbiden;
}				     
Hmm* define_full_model(int nstate, char **list, char *model_name, Generic_em_func emission_function)
{
  /*list: a list of the state names: does not include START or END*/
  /*model_name: a string that will be appended to the names in list*/

  Hmm *H;
  int a,n;
  HmmState *S;
  
  
  H=declare_hmm(nstate+2);
  H->freeT=0;
 
  H->forbiden=FORBIDEN;
  H->start=START_STATE;
  H->end=END_STATE;
 
  list=produce_state_name (nstate,list, model_name, H);
  nstate+=2;
  
  for (n=0; n<nstate; n++)
    {
      
      S=H->S[n];
      S->state=n;
      sprintf ( S->name, "%s", list[n]);
      S->em_func2=emission_function;
      if (n==H->end || n==H->start){S->DI=0;S->DJ=0;}
      else S->DI=1;S->DJ=0;
      
      /*Emmissions*/
      S->em_func2=emission_function;
      for (a=0; a< MAX_EMISSION; a++)S->em2[a]=H->freeT;
      
      for (a=0; a<nstate && n!=H->end; a++)
	{
	  if (a!=H->start && !(n==H->start && a==H->end) )
	  {
	    sprintf ( (S->T[S->nT])->name, "%s", list[a]);
	    (S->T[S->nT])->tr=H->freeT;
	    S->nT++;
	  }
	}
    }
  return H;
}

char **produce_state_name (int nstate,char **list, char *model_name, Hmm* H)
{
  int a,b,c;
  char **new_list;
  nstate+=2;
  
  new_list=declare_char ( nstate, 100);
  for ( a=0, b=0, c=0; a< nstate; a++)
    {
      if ( a==H->start)sprintf ( new_list[a], "START");
      else if ( a==H->end)sprintf ( new_list[a], "END");
      else if ( list==NULL){sprintf ( new_list[a], "%c%s", 'a'+b, (model_name)?model_name:"");b++;}
      else {sprintf ( new_list[a], "%c%s", list[a][c], (model_name)?model_name:"");c++;}
    }
  return new_list;
}

int seq_viterbi_pair_wise (Alignment *A,int*ns, int **ls,Constraint_list *CL)
{
  ungap(A->seq_al[0]);
  analyze_sequence (A->seq_al[0], NULL);
  myexit (EXIT_FAILURE);
  return 1;
}
double analyze_sequence ( char *seq, Hmm *H)
{

  double **M;
  int *path;

  if ( H==NULL)
    {
      H=define_full_model(5, NULL,"_first", pavie_emission);
      H=bound_hmm(H);
      DisplayHmm (H);
    }
  M=seq_viterbi_hmm (seq, H);
  path=seq_viterbi2path (seq, H, M);
  return M[H->end][strlen (seq)];
}
   

double** seq_viterbi_hmm (char *seq, Hmm *H)
{
  /*Given a model H and a sequence seq*/
  double **M;
  double e, v, max;
  int i,pi, bestk, s, k, l1;
  HmmState *S1, *S2;
  
  
  l1=strlen (seq);
  M=declare_double (H->nS*2,l1+2);
  
  /*Handle the start*/
  M[H->start][0]=0;
  for ( i=0; i<=l1; i++)
    {
      for ( s=0; s< H->nS; s++)
	{
	  S1=H->S[s];
	  pi=i-S1->DI;
	  max=H->forbiden;
	  bestk=H->forbiden;
	  if ( pi<0){M[s][i]=H->forbiden;}/*Boundary*/
	  else 
	    {
	      if (pi==0) {max=H->T[(int)H->start][s];bestk=H->start;}/*Start*/
	      else
		{
		  for (k=1; k<=H->fromM[S1->state][0]; k++)
		    {
		      S2=H->S[H->fromM[s][k]];
		      if ( S2->state==H->start || S2->state==H->end)continue;
		      v=hmm_add((M[S2->state][pi]),(H->T[S2->state][S1->state]));
		      if ( v!=H->forbiden && (max==H->forbiden || v>max)){max=v;bestk=S2->state;}
		    }
		} 	
	      if (S1->em2)e=S1->em2[(int)seq[pi]];
	      else e=S1->em_func2(H,S1, (int)seq[pi]);
	      
	      e=hmm_add (e,max);
	      
	      M[s][i]=e;
	      M[s+H->nS][i]=bestk;
	      }
	  }
      }   
  /*Terminate viterbi: connect the path to the END state*/
  max=UNDEFINED;
  bestk=UNDEFINED;
  for (k=0; k< H->nS; k++)
    {
      if (k==H->start || k==H->end);
      else
	{
	  v=(M[k][l1]==H->forbiden || H->T[k][H->end]==H->forbiden)?H->forbiden:M[k][l1]+H->T[k][H->end];
	  if ( max==H->forbiden || v>max){bestk=k;max=v;}
	}
    }
  M[H->end][l1]=max;
  M[H->nS+H->end][l1]=bestk; 
  return M;
}
  
int * seq_viterbi2path (char *s, Hmm *H, double **M)
{
  int i,l,l1;
  int *path;
  HmmState *S1;
  int cs;

  l1=strlen (s);
  path=vcalloc (l1+1, sizeof (int));
  i=l1;
  l=0;
  cs=M[H->nS+H->end][i];
  
  while (i>0)
    {
      
      S1=H->S[cs];
      path[l++]=cs;

      cs=M[H->nS+cs][i];
      i-=S1->DI;
      /*fprintf ( stderr, "%d", cs);*/
    }
  invert_list_int (path, l);
  path[l++]=H->forbiden;
  
  return path;
}

/*********************************************************************************/
/*                                                                               */
/*                                                                               */
/*                     pairHMM: Viterbi                                              */
/*                                                                               */
/*                                                                               */
/*********************************************************************************/
Hmm* define_mnm_model(Constraint_list *CL);
int viterbi_pair_wise_OLD (Alignment *A,int*ns, int **ls,Constraint_list *CL)
{
  int l1, l2, a;
  double ***M;
  int *path;
  Hmm * H;
  
  A->pos=aln2pos_simple( A, -1, ns, ls);
  
  //  H=define_mnm_model (CL);
  H=define_two_mat_model (CL);
  
  l1=strlen (A->seq_al[ls[0][0]]);
  l2=strlen (A->seq_al[ls[1][0]]);
  M=viterbi_hmm (A, ns, ls, H, CL);
  path=viterbi2path (l1,l2, H,M);
  A=viterbipath2aln (A,ns,ls,path, H); 
  A->score=A->score_aln=M[H->end][l1][l2];
  for ( a=0; a< H->nS*2; a++)free_double (M[a], -1);
  vfree (M);
  free_int (A->pos, -1);
  A->pos=NULL;

  free_Hmm (H);
  vfree (path);
  
  return A->score_aln;
}

Alignment * viterbipath2aln (Alignment *A, int *ns,int **ls,int *tb, Hmm *H)
{
  char **aln;
  char *char_buf;
  int a, b, c, len, ch;
  HmmState *S;
  int l[2];
  
  len=0;while (tb[len]!=H->forbiden)len++;
  
  if ( A->declared_len<=len)A=realloc_aln2  ( A,A->max_n_seq,2*len);
  aln=A->seq_al;
  
  char_buf=vcalloc (len+1, sizeof (char));
  l[0]=strlen ( A->seq_al[ls[0][0]]);
  l[1]=strlen ( A->seq_al[ls[1][0]]);
  
  for ( c=0; c< 2; c++)
    for ( a=0; a< ns[c]; a++) 
      {
	for (ch=0, b=0; b<len; b++)
	  {
	    S=H->S[tb[b]];
	    if ( (c==0 && S->DI)|| (c==1 && S->DJ) )
	      char_buf[b]=aln[ls[c][a]][ch++];
	    else 
	      char_buf[b]='-';
	  }
	char_buf[b]='\0';
	sprintf (aln[ls[c][a]],"%s", char_buf);
	if ( l[c]!=ch){fprintf (stderr, "\nERROR: Wrong Size Of Alignmnent (Real %d, Observed %d)[FATAL:%s]",l[c], ch, PROGRAM);}
      }
  A->len_aln=len;
  A->nseq=ns[0]+ns[1];
  
  vfree(char_buf);
  return A;
}

double*** viterbi_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL)
{
  double ***M;
  double e, v, max;
  int a, i,pi, bestk,j,pj, s, k, l1, l2;
  HmmState *S1, *S2;
  
  l1=strlen (A->seq_al[ls[0][0]]);
  l2=strlen (A->seq_al[ls[1][0]]);
  
  M=vcalloc (H->nS*2, sizeof (double**));
  for ( a=0; a<H->nS*2; a++)M[a]=declare_double (l1+2, l2+2); 
  
  /*Handle the start*/
  
  M[H->start][0][0]=0;
  for ( i=0; i<=l1; i++)
    for ( j=0; j<=l2; j++)
      {
	for ( s=0; s< H->nS; s++)
	  {
	    S1=H->S[s];
	    pi=i-S1->DI;
	    pj=j-S1->DJ;
	    max=H->forbiden;
	    bestk=H->forbiden;
	    if ( pi<0 ||pj<0){M[s][i][j]=H->forbiden;}/*Boundary*/
	    else 
	      {
		if (pi+pj==0) {max=H->T[H->start][s];bestk=H->start;}/*Start*/
		else
		  {
		    for (k=1; k<=H->fromM[S1->state][0]; k++)
		      {
			S2=H->S[H->fromM[s][k]];
			if ( S2->state==H->start || S2->state==H->end)continue;
			v=(M[S2->state][pi][pj]==H->forbiden)?H->forbiden:(M[S2->state][pi][pj]+H->T[S2->state][S1->state]);
			if ( v!=H->forbiden && (max==H->forbiden || v>max)){max=v;bestk=S2->state;}
		      }
		  }
	
		e=(S1->em==H->forbiden)?S1->em_func (A, A->pos, ns[0], ls[0],i-1, A->pos,ns[1], ls[1], j-1, CL):S1->em;
		e=(max==H->forbiden || e==H->forbiden)?H->forbiden:e+max;
		
		M[s][i][j]=e;
		M[s+H->nS][i][j]=bestk;
	      }
	  }
      }
  
  /*Terminate viterbi: connect the path to the END state*/
  max=UNDEFINED;
  bestk=UNDEFINED;
  for (k=0; k< H->nS; k++)
    {
      if (k==H->start || k==H->end);
      else
	{
	  v=(M[k][l1][l2]==H->forbiden || H->T[k][H->end]==H->forbiden)?H->forbiden:M[k][l1][l2]+H->T[k][H->end];
	  if ( max==H->forbiden || v>max){bestk=k;max=v;}
	}
    }
  M[H->end][l1][l2]=max;
  M[H->nS+H->end][l1][l2]=bestk; 
    
  return M;
}

/*********************************************************************************/
/*                                                                               */
/*                                                                               */
/*                      HMM: Decode/Traceback                                       */
/*                                                                               */
/*                                                                               */
/*********************************************************************************/
int * traceback (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL,MatState *S, MatState *E, int **seg_list)

{
  int *path;
  int l=0;
  MatState *N;
  int l1, l2;
  
  l1=strlen (A->seq_al[ls[0][0]]);
  l2=strlen (A->seq_al[ls[1][0]]);
  path=vcalloc ( l1+l2+1, sizeof (int));
  
  while ( S->st!=H->end)
    {
      DisplayMatState (S, "\n\tTraceback");
      N=S->n;
      if ( N && S && (((N->i-S->i)>1) ||((N->j-S->j)>1)))
	{
	  RviterbiD_hmm (A,ns,ls,H,CL,S,N,seg_list);
	  N=S->n;
	}

      path[l++]=S->st;
      ManageMatState (FREE,S);
      S=N;
    }
  
  path[l]=H->forbiden;
  return path;
}
	
int * viterbi2path (int l1,int l2, Hmm *H, double ***M)
{
  int i, j,l;
  int *path;
  HmmState *S1;
  int cs;

  l=0;
  path=vcalloc (l1+l2+1, sizeof (int));
  i=l1;j=l2;
  l=0;
  cs=M[H->nS+H->end][i][j];
  
  while (i>0|| j>0)
    {
      
      S1=H->S[cs];
      path[l++]=cs;

      cs=M[H->nS+cs][i][j];
      i-=S1->DI;
      j-=S1->DJ;
      /*fprintf ( stderr, "%d", cs);*/
    }
  invert_list_int (path, l);
  path[l++]=H->forbiden;
  
  return path;
}

/*********************************************************************************/
/*                                                                               */
/*                                                                               */
/*                      HMM Viterbi Linear                                       */
/*                                                                               */
/*                                                                               */
/*********************************************************************************/


int viterbiL_pair_wise (Alignment *A,int*ns, int **ls,Constraint_list *CL)
{
  int l1, l2;
  int *path;
  Hmm * H;
  MatState *Start;
  MatState *End;

  A->pos=aln2pos_simple( A, -1, ns, ls);
  Start=ManageMatState ( DECLARE, NULL);
  End=ManageMatState ( DECLARE, NULL);
  H=define_simple_model (CL);
  l1=strlen (A->seq_al[ls[0][0]]);
  l2=strlen (A->seq_al[ls[1][0]]);

  

  Start->i=0 ;Start->j=0 ; Start->st=H->start;Start->sc=0;
  End->i  =l1;  End->j=l2; End  ->st=H->end;
  Start=RviterbiL_hmm (A, ns, ls, H, CL, Start,End);
  path=traceback (A, ns, ls, H, CL, Start,NULL, NULL);
  
  A=viterbipath2aln (A,ns,ls,path, H); 
  
  free_Hmm (H);
  free_int (A->pos, -1);
  A->pos=NULL;
  
  return A->score_aln;
}


MatState* RviterbiL_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL,MatState *S, MatState *E)
{
  MatState *MS, *ME;
  MS=S;
  ME=E;
  
  viterbiL_hmm (A,ns, ls,H, CL, S, E);
  
  
  if ( S->n==E)return S;
  if ( E->sc==H->forbiden)
	{
	  DisplayHmm (H);
	  fprintf ( stderr, "\nERROR: The Requested Model (Cf Above) Cannot Produce the Pair-Alignment\nYou must allow extra possible transitions\n[FATAL:%s]", PROGRAM);
	  myexit ( EXIT_FAILURE);
	}
  E=S->n;
  
  while (S!=ME)
    {
      int d1, d2, align;
      d1=MinDeltaMatState(S,E);
      d2=MaxDeltaMatState(S,E);
      align=((d1==1 && d2==1) || ( d1==0))?0:1;
      if (align)RviterbiL_hmm (A,ns, ls,H, CL,S,E);
      S=E;
      E=S->n;
    }
  return MS;
}

#define Dim(i,j) (i*LenJ+j)
MatState* viterbiL_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL, MatState *S,MatState *E)
{
  int current, previous,row, prow;
  double v;
  int a,i,j,pi,pj, s, k;
  int start_i, start_j, end_i, end_j, l1, l2;
  HmmState *S1, *S2;
  MatState *CC, *PCC,*tS, *tE, *mark=NULL;
  int midpoint;


  static MatState ***M;


  static int LenJ, LenI;
  int MaxDelta=50, DeltaI, DeltaJ;

  

  DisplayMatState (S, "\n\tS");
  DisplayMatState (E, "\n\tE");
  
  
  if ( A==NULL)
    {
      for ( a=0; a<2; a++)memset(M[a],0,LenJ*LenI*sizeof (MatState*));
      free_arrayN((void **)M, 3);M=NULL;
      ManageMatState ( FREE_ALL, NULL);
      return NULL;
    }

  
  if ( MatStateAreIdentical ( S, E))return NULL;
  l1=strlen (A->seq_al[ls[0][0]]);l2=strlen (A->seq_al[ls[1][0]]);
  
  midpoint=S->i+((E->i-S->i)/2);
  DeltaI=E->i-S->i;
  DeltaJ=E->j-S->j;
  
  start_i=S->i;end_i=E->i;start_j=S->j;end_j=E->j;
  current=0;previous=1;
 
  
  if ( !M)
    {
      LenI=l2+1;
      LenJ=H->nS;
      M=declare_arrayN(3, sizeof ( MatState),2,LenI*LenJ,0);
    }
  
  
  /*MAKE THE VITERBI FROM S(tart) to E(nd)*/
  mark=ManageMatState ( MARK, mark);
  for (i=start_i; i<=end_i; i++)
    {
      row=current;
      for ( j=start_j; j<=end_j; j++)
	{
	  DeltaJ=((FABS(j-i))<MaxDelta)?1:MaxDelta+1; /*Keep Pointers on a DealtaMax Band around the main diagonal*/
	  for ( s=H->nS-1;s>=0; s--)
	    {
	      S1=H->S[s];pi=i-S1->DI;prow=S1->DI;pj=j-S1->DJ;
	      	      
	      CC=M[row][Dim(j,s)]=CopyMatState(NULL, M[row][Dim(j,s)]);
	      CC->i=i; CC->j=j; CC->st=s;PCC=NULL;
	      	      
	      if (i==start_i && j==start_j && s==S->st){CC=CopyMatState(S,CC);}
	      else if ( i==end_i && j==end_j && E->st!=H->end && s!=E->st)CC->sc=H->forbiden;
	      else if ( pi<start_i || pj<start_j)        {CC->sc=H->forbiden;}
	      else 
		{
		  for (k=1; k<=H->fromM[S1->state][0]; k++)
		    {
		      S2=H->S[H->fromM[s][k]];
		      PCC=M[prow][Dim((j-S1->DJ),(S2->state))];
		      		    
		      if ( !PCC)PCC=NULL;
		      else if      ( pi+pj!=0 && S2->state==H->start);
		      else if ( !(pi==l1 && pj==l2) && s==H->end);
		      else
			{
			  
			  v=hmm_add(CC->sc,H->T[PCC->st][CC->st]);
			  
			  v=lu_RviterbiD_hmm(A,ns, ls, H, CL,PCC,CC, NULL);
			  if ( v!=H->forbiden && (CC->sc==H->forbiden || v> CC->sc)){CC->sc=v; CC->pst=S2->state;CC->p=PCC;}
			}
		    }
		}
	      if (CC->sc==H->forbiden);
	      else if (i==midpoint || DeltaI<=MaxDelta||DeltaJ<=MaxDelta ||(i==start_i && j==start_j && s==S->st) )
		{
		  CC->m=(CC->p)?(CC->p)->m:NULL;
		  PCC=CopyMatState(CC,NULL);
		  PCC->m=CC->m;CC->m=PCC;
		}
	      else CC->m=(CC->p)?(CC->p)->m:NULL;
	    }
	}
      prow=previous;
      for ( j=start_j; j<=end_j && i!=end_i; j++)
	{
	  for ( s=H->nS-1;s>=0; s--)
	    {
	      
	      CC=(M[prow][Dim(j,s)]);M[prow][Dim(j,s)]=M[row][Dim(j,s)];M[row][Dim(j,s)]=CC;
	      if (M[prow][Dim(j,s)]) M[row ][Dim(j,s)]=CopyMatState ( M[prow][Dim(j,s)], M[row][Dim(j,s)]);
	      
	    }
	} 
      
    }  
  
  mark=ManageMatState ( MARK,mark);
  row=current;
  
  
  if ( E->st==H->end || E->st==H->forbiden){E=CopyMatState ((M[row][Dim(end_j,E->st)]),E);}

  
  

  PCC=CopyMatState (M[row][Dim(end_j,E->st)], NULL);
  
  if ( MatStateAreIdentical(PCC,PCC->m))PCC=PCC->m;
    tS=tE=PCC;
  while (PCC->m)
    {
      tS=CopyMatState (PCC->m,NULL); tS->n=PCC; PCC->p=tS;PCC=tS;
    }
    
  if (tS==tE);
  else
    {
      S->n=tS->n; (S->n)->p=S;
      E->p=tE->p; (E->p)->n=E;
    }
  for ( a=0; a<2; a++)memset(M[a],0,LenJ*LenI*sizeof (MatState*));
  ManageMatState ( FREE_MARK,mark);
  
    
  while (S && S->p!=E){S->m=NULL;S=S->n;}/*Clean the memory of the rturned Cells*/
  return NULL;
}

/*********************************************************************************/
/*                                                                               */
/*                                                                               */
/*                      HMM Viterbi Diagonals                                    */
/*                                                                               */
/*                                                                               */
/*********************************************************************************/
int viterbiD_pair_wise (Alignment *A,int*ns, int **ls,Constraint_list *CL)
{
  int l1, l2;
  int *path;
  Hmm * H;
  MatState *Start;
  MatState *End;
  int **seg_list;
  int a, b, c;
  int main_i;
  int main_j;

  
  A->pos=aln2pos_simple( A, -1, ns, ls);
  
  Start=ManageMatState ( DECLARE, NULL);
  End=ManageMatState ( DECLARE, NULL);
  H=define_simple_model (CL);
  l1=strlen (A->seq_al[ls[0][0]]);
  l2=strlen (A->seq_al[ls[1][0]]);

  main_i=MAX(1,(l2-l1)+1);
  main_j=MAX(1,(l1-l2)+1);
  
  seg_list=declare_arrayN(2, sizeof (int), l1+l2+3, 3);
  seg_list[0][0]=DIAGONALS;

  
  c=1;
  for ( b=1,a=l1; a>= 1; a--)
    {
      if (a<50 || (b==main_i && a==main_j))
	{
	seg_list[c][0]=a;
	seg_list[c][1]=b;
	seg_list[c][2]=MIN((l1-a), (l2-b));
	c++;
	}
      }
  
  
  for ( b=2,a=1; b<= l2; b++, c++)
    {
      if (b<50 || (b==main_i && a==main_j))
	{
	  seg_list[c][0]=a;
	  seg_list[c][1]=b;
	  seg_list[c][2]=MIN((l1-a), (l2-b));
	}
      }


  seg_list[c][0]=FORBIDEN;
  
  Start->i=0 ;Start->j=0 ; Start->st=H->start;Start->sc=0;
  End->i  =l1;  End->j=l2; End  ->st=H->end;
  Start=RviterbiD_hmm (A, ns, ls, H, CL, Start,End,seg_list);
  
    
  path=traceback (A, ns, ls, H, CL, Start,NULL, NULL);

  

  A=viterbipath2aln (A,ns,ls,path, H); 
  
  viterbiD_hmm (NULL, ns, ls, H, CL, Start,End, seg_list);
  free_Hmm (H);
  free_int (A->pos, -1);
  free_arrayN((void **)seg_list, 2);
    
  A->pos=NULL;
  return A->score_aln;
}


double lu_RviterbiD_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL,MatState *S, MatState *E, int **seg_list)
{
  HmmState *S1;
  double sc, sc2,e, t;
  static MatState *cS=NULL, *cE=NULL;
  double min, max;
  max=MAX((E->i-S->i), (E->j-S->j));
  min=MIN((E->i-S->i), (E->j-S->j));
  
  
  if ( S->sc==H->forbiden) return H->forbiden;
  else if (min==0)
    {
      e=hmm_add(S->sc,H->T[S->st][E->st]);
      if (  H->T[E->st][E->st]!=H->forbiden)e=hmm_add(e, (max-1)*H->T[E->st][E->st]);
      if ( (H->S[E->st])->em!=H->forbiden)  e=hmm_add(e, max    *(H->S[E->st])->em );
      return e;
    }
  else if ( min>0 && max>1)
    { 
      
      fprintf ( stderr, "\nWarning: Disjoined Diagonals");
      DisplayMatState (S, "\n\tS");
      DisplayMatState (E, "\n\tE");
      
      
      cS=CopyMatState ( S,cS);
      cE=CopyMatState ( E,cE);
      cE->sc=H->forbiden;
      viterbiD_hmm (A,ns,ls, H,CL,cS, cE, NULL);
      sc2=cE->sc;
     
      return sc2;
    }
  else
    {
      S1=H->S[E->st];
      t=H->T[S->st][E->st];
      e=(S1->em==H->forbiden)?S1->em_func (A, A->pos, ns[0], ls[0],E->i-1, A->pos,ns[1], ls[1], E->j-1, CL):S1->em;
      sc=hmm_add(S->sc,t);
      sc=hmm_add(sc,e);
      return sc;
    }
  return H->forbiden;
}


MatState* RviterbiD_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL,MatState *S, MatState *E, int **seg_list)
{
  MatState *MS, *ME;
  MS=S;
  ME=E;
  
  viterbiD_hmm (A,ns, ls,H, CL, S, E, seg_list);
  
  
  if ( S->n==E)return S;
  if ( E->sc==H->forbiden)
	{
	  DisplayHmm (H);
	  fprintf ( stderr, "\nERROR: The Requested Model (Cf Above) Cannot Produce the Pair-Alignment\nYou must allow extra possible transitions\n[FATAL:%s]", PROGRAM);
	  myexit ( EXIT_FAILURE);
	}
  E=S->n;
  
  while (S!=ME)
    {
      int d1, d2, align;
      d1=MinDeltaMatState(S,E);
      d2=MaxDeltaMatState(S,E);
      align=((d1==1 && d2==1) || ( d1==0))?0:1;
      if (align)RviterbiD_hmm (A,ns, ls,H, CL,S,E, seg_list);
      S=E;
      E=S->n;
    }
  return MS;
}

#define Dim(i,j) (i*LenJ+j)
MatState* viterbiD_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL, MatState *S,MatState *E, int **seg_list)
{
  int current, previous,row, prow;
  double v;
  int a,b,i,j,pi,pj, s, k;
  int start_i, start_j, end_i, end_j, l1, l2;
  HmmState *S1, *S2;
  MatState *CC, *PCC,*tS, *tE, *mark=NULL;
  int midpoint;
  
  int dj;
  int dc;
  int *jlist=NULL;
  static int **main_jlist;
  static MatState ***M;
  static int *toclean;
  int ntoclean;
  static int LenJ, LenI;
  int MaxDelta=50, DeltaI, DeltaJ;
  int mode;

  DisplayMatState (S, "\n\tS");
  DisplayMatState (E, "\n\tE");
  
  
  if ( A==NULL)
    {
      free_arrayN((void **)main_jlist, 2);main_jlist=NULL;

      for ( a=0; a<2; a++)memset(M[a],0,LenJ*LenI*sizeof (MatState*));
      free_arrayN((void **)M, 3);M=NULL;
      vfree (toclean);
      ManageMatState ( FREE_ALL, NULL);
      return NULL;
    }

  
  if ( MatStateAreIdentical ( S, E))return NULL;
  l1=strlen (A->seq_al[ls[0][0]]);l2=strlen (A->seq_al[ls[1][0]]);
  
  midpoint=S->i+((E->i-S->i)/2);
  DeltaI=E->i-S->i;
  
  
  start_i=S->i;end_i=E->i;start_j=S->j;end_j=E->j;
  current=0;previous=1;
 
  
  if ( !M)
    {
      LenI=l2+1;
      LenJ=H->nS;
      M=declare_arrayN(3, sizeof ( MatState),2,LenI*LenJ,0);
      toclean=vcalloc ( LenI*LenJ, sizeof (int));
    }
  
  if ( !main_jlist)main_jlist= seglist2table(seg_list, l1, l2);
  

  /*MAKE THE VITERBI FROM S(tart) to E(nd)*/
  mark=ManageMatState ( MARK, mark);
  mode=(!seg_list)?ALL:seg_list[0][0];
  
  for (ntoclean=0,i=start_i; i<=end_i; i++)
    {
      row=current;
      
      if ( mode==ALL)jlist=main_jlist[0];
      else if ( mode==DIAGONALS)jlist=(i==0)?main_jlist[0]:main_jlist[1];
      else if ( mode==SEGMENTS)  jlist=main_jlist[i+2];
      
     
      for ( dj=1; dj<=jlist[0]; dj++)
	{
	  DeltaJ=((FABS(dj-i))<MaxDelta)?1:MaxDelta+1; /*Keep Pointers on a 2*DeltaMax Band around the main diagonal*/

	  dc=(mode==DIAGONALS && dj!=1)?i:0;/*Make sure the diagonal Mode uses the sides of the array*/
	  j=jlist[dj]+dc;
	  
	  
	  if ( j<start_j || j>end_j)continue;
	  for ( s=H->nS-1;s>=0; s--)
	    {
	      S1=H->S[s];pi=i-S1->DI;prow=S1->DI;
	      
	      if ( S1->DI && S1->DJ){pj=j-S1->DJ;}
	      else if ( !S1->DJ)pj=j;
	      else if ( dj>1)pj=jlist[dj-S1->DJ]+dc;
	      else pj=-1;
	      
	      if (!M[row][Dim(j,s)])toclean[ntoclean]=Dim(j,s);
	      
	      CC=M[row][Dim(j,s)]=CopyMatState(NULL, M[row][Dim(j,s)]);
	      CC->i=i; CC->j=j; CC->st=s;PCC=NULL;
	      	      
	      if (i==start_i && j==start_j && s==S->st){CC=CopyMatState(S,CC);}
	      else if ( i==end_i && j==end_j && E->st!=H->end && s!=E->st)CC->sc=H->forbiden;
	      else if ( pi<start_i || pj<start_j)        {CC->sc=H->forbiden;}
	      else 
		{
		  for (k=1; k<=H->fromM[S1->state][0]; k++)
		    {
		      S2=H->S[H->fromM[s][k]];
		      
		      if ( S1->DI && S1->DJ)PCC=M[prow][Dim((j-S1->DJ),(S2->state))];
		      else PCC=M[prow][Dim((jlist[dj-S1->DJ]+dc),(S2->state))];
		    
		      if ( !PCC)PCC=NULL;
		      else if      ( pi+pj!=0 && S2->state==H->start);
		      else if ( !(pi==l1 && pj==l2) && s==H->end);
		      else
			{
			  v=lu_RviterbiD_hmm(A,ns, ls, H, CL,PCC,CC, NULL);
			  if ( v!=H->forbiden && (CC->sc==H->forbiden || v> CC->sc)){CC->sc=v; CC->pst=S2->state;CC->p=PCC;}
			}
		    }
		}
	      if (CC->sc==H->forbiden);
	      else if (i==midpoint || DeltaI<=MaxDelta||DeltaJ<=MaxDelta ||(i==start_i && j==start_j && s==S->st) )
		{
		  CC->m=(CC->p)?(CC->p)->m:NULL;
		  PCC=CopyMatState(CC,NULL);
		  PCC->m=CC->m;CC->m=PCC;
		}
	      else CC->m=(CC->p)?(CC->p)->m:NULL;
	    }
	}
      prow=previous;
      for ( dj=1; dj<=jlist[0] && i!=end_i; dj++)
	{
	  dc=(mode==DIAGONALS && dj!=1)?i:0;
	  j=jlist[dj]+dc;
	  if ( j<start_j || j>end_j)continue;
	  
	  for ( s=H->nS-1;s>=0; s--)
	    {
	      
	      CC=(M[prow][Dim(j,s)]);M[prow][Dim(j,s)]=M[row][Dim(j,s)];M[row][Dim(j,s)]=CC;
	      if (!M[row][Dim(j,s)])toclean[ntoclean++]=Dim(j,s);
	      if (M[prow][Dim(j,s)]) M[row ][Dim(j,s)]=CopyMatState ( M[prow][Dim(j,s)], M[row][Dim(j,s)]);
	      
	    }
	} 
      
    }  
  
  mark=ManageMatState ( MARK,mark);
  row=current;
  
  
  if ( E->st==H->end || E->st==H->forbiden){E=CopyMatState ((M[row][Dim(end_j,E->st)]),E);}

  
  

  PCC=CopyMatState (M[row][Dim(end_j,E->st)], NULL);
  
  if ( MatStateAreIdentical(PCC,PCC->m))PCC=PCC->m;
    tS=tE=PCC;
  while (PCC->m)
    {
      tS=CopyMatState (PCC->m,NULL); tS->n=PCC; PCC->p=tS;PCC=tS;
    }
    
  if (tS==tE);
  else
    {
      S->n=tS->n; (S->n)->p=S;
      E->p=tE->p; (E->p)->n=E;
    }

  ManageMatState ( FREE_MARK,mark);
  
  
  for ( a=0; a<ntoclean; a++)
    {
      for ( b=0; b< 2; b++){M[b][toclean[a]]=NULL;}
      toclean[a]=0;
    }
  
  while (S && S->p!=E){S->m=NULL;S=S->n;}/*Clean the memory of the rturned Cells*/
  return NULL;
}

/*********************************************************************************/
/*                                                                               */
/*                                                                               */
/*                      HMM Viterbi Diagonals  GLOBAL/LOCAL                                  */
/*                                                                               */
/*                                                                               */
/*********************************************************************************/

int viterbiDGL_pair_wise (Alignment *A,int*ns, int **ls,Constraint_list *CL)
{
  int l1, l2;
  int *path;
  Hmm * H;
  MatState *Start;
  MatState *End;
  int **seg_list;
  int a, b, c;
  int main_i;
  int main_j;

  
  A->pos=aln2pos_simple( A, -1, ns, ls);
  
  Start=ManageMatState ( DECLARE, NULL);
  End=ManageMatState ( DECLARE, NULL);
  H=define_simple_model (CL);
  l1=strlen (A->seq_al[ls[0][0]]);
  l2=strlen (A->seq_al[ls[1][0]]);

  main_i=MAX(1,(l2-l1)+1);
  main_j=MAX(1,(l1-l2)+1);
  
  seg_list=declare_arrayN(2, sizeof (int), l1+l2+3, 3);
  seg_list[0][0]=DIAGONALS;

  
  c=1;
  for ( b=1,a=l1; a>= 1; a--)
    {
      if (a<50 || (b==main_i && a==main_j))
	{
	seg_list[c][0]=a;
	seg_list[c][1]=b;
	seg_list[c][2]=MIN((l1-a), (l2-b));
	c++;
	}
      }
  
  
  for ( b=2,a=1; b<= l2; b++, c++)
    {
      if (b<50 || (b==main_i && a==main_j))
	{
	  seg_list[c][0]=a;
	  seg_list[c][1]=b;
	  seg_list[c][2]=MIN((l1-a), (l2-b));
	}
      }


  seg_list[c][0]=FORBIDEN;
  
  Start->i=0 ;Start->j=0 ; Start->st=H->start;Start->sc=0;
  End->i  =l1;  End->j=l2; End  ->st=H->end;
  Start=RviterbiDGL_hmm (A, ns, ls, H, CL, Start,End,seg_list);
  
    
  path=traceback (A, ns, ls, H, CL, Start,NULL, NULL);

  

  A=viterbipath2aln (A,ns,ls,path, H); 
  
  viterbiD_hmm (NULL, ns, ls, H, CL, Start,End, seg_list);
  free_Hmm (H);
  free_int (A->pos, -1);
  free_arrayN((void **)seg_list, 2);
    
  A->pos=NULL;
  return A->score_aln;
}


double lu_RviterbiDGL_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL,MatState *S, MatState *E, int **seg_list)
{
  HmmState *S1;
  double sc, sc2,e, t;
  static MatState *cS=NULL, *cE=NULL;
  double min, max;
  max=MAX((E->i-S->i), (E->j-S->j));
  min=MIN((E->i-S->i), (E->j-S->j));
  
  
  
  
  if ( S==NULL || E==NULL || S->sc==H->forbiden) return H->forbiden;
  else if ( S->st==H->start) return 0;
  else if ( E->st==H->end) return S->sc;
  else if (min==0)
    {
      e=hmm_add(S->sc,H->T[S->st][E->st]);
      if (  H->T[E->st][E->st]!=H->forbiden)e=hmm_add(e, (max-1)*H->T[E->st][E->st]);
      if ( (H->S[E->st])->em!=H->forbiden)  e=hmm_add(e, max    *(H->S[E->st])->em );
      return e;
    }
  else if ( min>0 && max>1)
    { 
      
      fprintf ( stderr, "\nWarning: Disjoined Diagonals");
      DisplayMatState (S, "\n\tS");
      DisplayMatState (E, "\n\tE");
      
      
      cS=CopyMatState ( S,cS);
      cE=CopyMatState ( E,cE);
      cE->sc=H->forbiden;
      viterbiD_hmm (A,ns,ls, H,CL,cS, cE, NULL);
      sc2=cE->sc;
     
      return sc2;
    }
  else
    {
      S1=H->S[E->st];
      t=H->T[S->st][E->st];
      e=(S1->em==H->forbiden)?S1->em_func (A, A->pos, ns[0], ls[0],E->i-1, A->pos,ns[1], ls[1], E->j-1, CL):S1->em;
      sc=hmm_add(S->sc,t);
      sc=hmm_add(sc,e);
      return sc;
    }
  return H->forbiden;
}


MatState* RviterbiDGL_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL,MatState *S, MatState *E, int **seg_list)
{
  MatState *MS, *ME;
  MS=S;
  ME=E;
  
  viterbiDGL_hmm (A,ns, ls,H, CL, S, E, seg_list);
  
  
  if ( S->n==E)return S;
  if ( E->sc==H->forbiden)
	{
	  DisplayHmm (H);
	  fprintf ( stderr, "\nERROR: The Requested Model (Cf Above) Cannot Produce the Pair-Alignment\nYou must allow extra possible transitions\n[FATAL:%s]", PROGRAM);
	  myexit ( EXIT_FAILURE);
	}
  E=S->n;
  
  while (S!=ME)
    {
      int d1, d2, align;
      d1=MinDeltaMatState(S,E);
      d2=MaxDeltaMatState(S,E);
      align=((d1==1 && d2==1) || ( d1==0))?0:1;
      if (align)RviterbiDGL_hmm (A,ns, ls,H, CL,S,E, seg_list);
      S=E;
      E=S->n;
    }
  return MS;
}


#define Dim(i,j) (i*LenJ+j)
MatState* viterbiDGL_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL, MatState *S,MatState *E, int **seg_list)
{
  int current, previous,row, prow;
  double v;
  int a,i,j,pi,pj, s, k;
  int start_i, start_j, end_i, end_j, l1, l2;
  HmmState *S1, *S2;
  MatState *CC, *PCC,*tS, *tE,*bestE,*bestS, *mark=NULL;
  int midpoint;
  
  int dj;
  int dc;
  int *jlist=NULL;
  static int **main_jlist;
  static MatState ***M;
  static int *toclean;
  int ntoclean;
  static int LenJ, LenI;
  int MaxDelta=50, DeltaI, DeltaJ;
  int mode;

  

  DisplayMatState (S, "\n\tS");
  DisplayMatState (E, "\n\tE");
  
  
  if ( A==NULL)
    {
      free_arrayN((void **)main_jlist, 2);main_jlist=NULL;

      for ( a=0; a<2; a++)memset(M[a],0,LenJ*LenI*sizeof (MatState*));
      free_arrayN((void **)M, 3);M=NULL;
      vfree (toclean);
      ManageMatState ( FREE_ALL, NULL);
      return NULL;
    }

  
  if ( MatStateAreIdentical ( S, E))return NULL;
  l1=strlen (A->seq_al[ls[0][0]]);l2=strlen (A->seq_al[ls[1][0]]);
  
  midpoint=S->i+((E->i-S->i)/2);
  DeltaI=E->i-S->i;
  
  
  start_i=S->i;end_i=E->i;start_j=S->j;end_j=E->j;
  current=0;previous=1;
 
  
  if ( !M)
    {
      LenI=l2+1;
      LenJ=H->nS;
      M=declare_arrayN(3, sizeof ( MatState),2,LenI*LenJ,0);
      toclean=vcalloc ( LenI*LenJ, sizeof (int));
    }
  
  if ( !main_jlist)main_jlist= seglist2table(seg_list, l1, l2);
  

  /*MAKE THE VITERBI FROM S(tart) to E(nd)*/
  mark=ManageMatState ( MARK, mark);
  mode=(!seg_list)?ALL:seg_list[0][0];
  bestE=CopyMatState (E, NULL);
  bestS=CopyMatState (NULL, NULL);
  for (ntoclean=0,i=start_i; i<=end_i; i++)
    {
      row=current;
      
      if ( mode==ALL)jlist=main_jlist[0];
      else if ( mode==DIAGONALS)jlist=(i==0)?main_jlist[0]:main_jlist[1];
      else if ( mode==SEGMENTS)  jlist=main_jlist[i+2];
      
     
      for ( dj=1; dj<=jlist[0]; dj++)
	{
	  DeltaJ=(FABS(dj-i)<MaxDelta)?1:MaxDelta+1; /*Keep Pointers on a 2*DeltaMax Band around the main diagonal*/

	  dc=(mode==DIAGONALS && dj!=1)?i:0;/*Make sure the diagonal Mode uses the sides of the array*/
	  j=jlist[dj]+dc;
	  
	  
	  if ( j<start_j || j>end_j)continue;
	  for ( s=H->nS-1;s>=0; s--)
	    {
	      if ( s==S->st)continue;
	      S1=H->S[s];pi=i-S1->DI;prow=S1->DI;
	      
	      if ( S1->DI && S1->DJ){pj=j-S1->DJ;}
	      else if ( !S1->DJ)pj=j;
	      else if ( dj>1)pj=jlist[dj-S1->DJ]+dc;
	      else pj=-1;
	      
	      if (!M[row][Dim(j,s)])toclean[ntoclean]=Dim(j,s);
	      
	      CC=M[row][Dim(j,s)]=CopyMatState(NULL, M[row][Dim(j,s)]);
	      CC->i=i; CC->j=j; CC->st=s;PCC=NULL;
	      	      
	      if (i==start_i && j==start_j && s==S->st){CC=CopyMatState(S,CC);}
	      else if ( s==S->st);
	      else if ( i==end_i && j==end_j && E->st!=H->end && s!=E->st)CC->sc=H->forbiden;
	      
	      else if ( pi<start_i || pj<start_j)        {CC->sc=H->forbiden;}
	      else 
		{
		  for (k=1; k<=H->fromM[S1->state][0]; k++)
		    {
		      S2=H->S[H->fromM[s][k]];
		      
		      if ( S1->DI && S1->DJ)PCC=M[prow][Dim((j-S1->DJ),(S2->state))];
		      else PCC=M[prow][Dim((jlist[dj-S1->DJ]+dc),(S2->state))];
		    
		      if ( S2->state==H->start){PCC=bestS;PCC->st=0;PCC->sc=0;PCC->m=PCC->n=PCC->p=NULL;}
		      
		      v=lu_RviterbiDGL_hmm(A,ns, ls, H, CL,PCC,CC, NULL);
		      if ( v!=H->forbiden && (CC->sc==H->forbiden || v> CC->sc)){CC->sc=v; CC->pst=S2->state;CC->p=PCC;}
		    }
		}
	      if ( CC->sc==H->forbiden);
	      else if ( bestE->sc==H->forbiden || bestE->sc>CC->sc) 
		{
		  bestE=CopyMatState(CC, bestE);
		  bestE->m=(CC->p)->m;
		}
	      else if (CC->p && (CC->p)->st==H->start)
		{
		  CC->m=CopyMatState (CC->p, NULL);		  
		}
	      else if (i==midpoint || DeltaI<=MaxDelta||DeltaJ<=MaxDelta ||(i==start_i && j==start_j && s==S->st) )
		{
		  CC->m=(CC->p)?(CC->p)->m:NULL;
		  PCC=CopyMatState(CC,NULL);
		  PCC->m=CC->m;CC->m=PCC;
		}
	      else CC->m=(CC->p)?(CC->p)->m:NULL;
	    }
	}
      prow=previous;
      for ( dj=1; dj<=jlist[0] && i!=end_i; dj++)
	{
	  dc=(mode==DIAGONALS && dj!=1)?i:0;
	  j=jlist[dj]+dc;
	  if ( j<start_j || j>end_j)continue;
	  
	  for ( s=H->nS-1;s>=0; s--)
	    {
	      
	      CC=(M[prow][Dim(j,s)]);M[prow][Dim(j,s)]=M[row][Dim(j,s)];M[row][Dim(j,s)]=CC;
	      /*if (!M[row][Dim(j,s)])toclean[ntoclean++]=Dim(j,s);*/
	      if (M[prow][Dim(j,s)]) M[row ][Dim(j,s)]=CopyMatState ( M[prow][Dim(j,s)], M[row][Dim(j,s)]);
	      
	    }
	} 
      
    }  
  
  mark=ManageMatState ( MARK,mark);
  row=current;
  
  
  if ( E->st==H->end || E->st==H->forbiden){E=CopyMatState ((M[row][Dim(end_j,E->st)]),E);}
  PCC=CopyMatState (bestE, NULL);
  
  if ( MatStateAreIdentical(PCC,PCC->m))PCC=PCC->m;
    tS=tE=PCC;
  while (PCC->m)
    {
      tS=CopyMatState (PCC->m,NULL); tS->n=PCC; PCC->p=tS;PCC=tS;
    }
    
  if (tS==tE);
  else
    {
      CopyMatState ( tS, S);
      CopyMatState ( tE, E);
    }
  ManageMatState ( FREE_MARK,mark);
  
  for ( a=0; a<2; a++)memset(M[a],0,LenJ*LenI*sizeof (MatState*));
  
  while (S && S->p!=E){S->m=NULL;S=S->n;}/*Clean the memory of the rturned Cells*/
  return NULL;
}


/*********************************************************************************/
/*                                                                               */
/*                                                                               */
/*                      HMM Viterbi Diagonals  PROCESSING                        */
/*                                                                               */
/*                                                                               */
/*********************************************************************************/
int **seglist2table ( int **seglist,int l1, int l2)
  {
    int **valuesT;
    int *bvalues;
    int line, a,si, sj, ei, j, c;
    
    /*All: 0*/
    valuesT=vcalloc ((l1+2)+3, sizeof (int*));
    valuesT[0]=vcalloc (l2+2, sizeof (int));
    for (a=0; a<=l2; a++)valuesT[0][++valuesT[0][0]]=a;

    if ( !seglist) return valuesT;
    /*Diagonals: 1*/
    valuesT[1]=vcalloc (l1+l2+2, sizeof (int));
    bvalues=vcalloc (l1+l2+2, sizeof (int));
    c=1;
    while (seglist[c][0]!=FORBIDEN)
      {
	
	si=seglist[c][0];
	sj=seglist[c][1];
	
	bvalues[(sj-si)+l1]=1;
	c++;
      }
    valuesT[1][++valuesT[1][0]]=0;
    for (a=0; a<=(l1+l2); a++)
      {
	if (bvalues[a])
	  {
	    valuesT[1][++valuesT[1][0]]=a-l1;
	  }
	
      }
    vfree (bvalues);

    /*Segments: 2*/
    valuesT[2]=vcalloc (l2+2, sizeof (int));
    for (a=0; a<=l2; a++)valuesT[2][++valuesT[2][0]]=a;
    
    bvalues=vcalloc (l2+2, sizeof (int));
    for ( line=1; line<=l1; line++)
      {
	bvalues[0]=c=0;
	bvalues[++bvalues[0]]=0;
	while (seglist[c][0]!=FORBIDEN)
	  {
		si=seglist[c][0];
		ei=si+seglist[c][2];
		sj=seglist[c][1];
		j=sj+(line-si);
		if ( line<si || line>ei);
		else if (j>=0 && j<=l2 && seglist[c][2])
		  {
		    bvalues[++bvalues[0]]=j;
		  }
		c++;
	  }
	valuesT[line+2]=vcalloc (bvalues[0]+1, sizeof (int));
	for ( a=0; a<=bvalues[0]; a++) valuesT[line+2][a]=bvalues[a];
      }
    vfree (bvalues);
    return valuesT;
  }


  
/*********************************************************************************/
/*                                                                               */
/*                                                                               */
/*                      HMM modeling                                             */
/*                                                                               */
/*                                                                               */
/*********************************************************************************/

Hmm* declare_hmm(int n)
  {
    Hmm *H;
    int a, b;
    
    H=vcalloc (1, sizeof (Hmm));
    H->nS=n;
    H->S=vcalloc (H->nS, sizeof (HmmState*));
    for (a=0; a<H->nS; a++)
      {
	H->S[a]=vcalloc (1, sizeof (HmmState));
	(H->S[a])->em2=vcalloc (MAX_EMISSION, sizeof (double));
	
	(H->S[a])->T=vcalloc ( H->nS, sizeof (StateTrans*));
	for ( b=0; b< H->nS; b++)
	  (H->S[a])->T[b]=vcalloc (1, sizeof (StateTrans));
      }
    return H;
  }

Hmm* free_Hmm(Hmm*H)
  {
    int a, b;
    
    H=vcalloc (1, sizeof (Hmm));
    free_double (H->T, -1);
    free_int ( H->fromM, -1);
    free_int ( H->toM, -1);
    
    for (a=0; a< H->nS; a++)
      {
	
	for ( b=0; b< H->nS; b++)
	  {
	    vfree ((H->S[a])->em2);
	    vfree((H->S[a])->T[b]);
	  }
	vfree((H->S[a])->T);
	vfree(H->S[a]);
      }
    vfree (H->S);
    vfree (H);
    return NULL;
  }

void DisplayHmm ( Hmm *H)
{
  int a, b;
  HmmState *S1, *S2;
  
  for ( a=0; a< H->nS; a++)
    {
      S1=H->S[a];
      fprintf ( stderr, "\nSTATE %d: %s\n",S1->state,S1->name);
      fprintf ( stderr, "\n\tDI %d", S1->DI);
      fprintf ( stderr, "\n\tDJ %d", S1->DJ);
      fprintf ( stderr, "\n\tE  %f", (float)S1->em);
            
      fprintf ( stderr, "\nReached FROM: ");
      for ( b=1; b<=H->fromM[a][0]; b++)
	{
	  S2=H->S[H->fromM[a][b]];
	  fprintf ( stderr, "[ %s %f ] ", S2->name, H->T[S2->state][S1->state]);
	}
      fprintf ( stderr, "\nGoes TO: ");
      for ( b=1; b<=H->toM[a][0]; b++)
	{
	  S2=H->S[H->toM[a][b]];
	  fprintf ( stderr, "[ %s %f ] ", S2->name, H->T[S1->state][S2->state]);
	}
    }
  return;
}
Hmm * bound_hmm ( Hmm *H)
{
  int a, b, c;
  char **name;
  HmmState *S;
  
  name=declare_char(H->nS, 100);
  H->T=declare_double ( H->nS, H->nS);
  
  for ( a=0; a< H->nS; a++)
    {
      sprintf ( name[a], "%s", (H->S[a])->name);
      H->order=MAX(H->order, (H->S[a])->DI);
      H->order=MAX(H->order, (H->S[a])->DJ);
    }
  
  for ( a=0; a< H->nS; a++)for (b=0; b<H->nS; b++)H->T[a][b]=H->forbiden;
  for (a=0; a< H->nS; a++)
    {
      S=H->S[a];
      for ( b=0; b< S->nT; b++)
	{
	  c=name_is_in_list ((S->T[b])->name, name, H->nS, 100);
	  if ( c!=-1)H->T[a][c]=(S->T[b])->tr;
	}
    }
  
  /*Bound the model:
    bM[state][0]=n_allowed transitions
    bM[state][1]=first allowed transition
  */

  H->toM=declare_int ( H->nS, H->nS);
  H->fromM=declare_int ( H->nS, H->nS);
  
  for ( a=0; a< H->nS; a++)
    for ( b=0; b< H->nS; b++)
      {
	if ( H->T[a][b]!=H->forbiden )
	  {
	    {H->fromM[b][0]++; H->fromM[b][H->fromM[b][0]]=a;}
	    {H->toM[a][0]++; H->toM[a][H->toM[a][0]]=b;}
	  }
      }
  for ( a=0; a< H->nS; a++)
    {
      if (( H->S[a])->em!=H->forbiden)( H->S[a])->em*=SCORE_K;
      for ( b=0; b< H->nS; b++)
	if ( H->T[a][b]!=H->forbiden)H->T[a][b]*=SCORE_K;
    }
  free_arrayN((void**)name, 2);
  return H;
}

/*********************************************************************************/
/*                                                                               */
/*                                                                               */
/*                      Memory Management                                        */
/*                                                                               */
/*                                                                               */
/*********************************************************************************/

MatState * ManageMatState(int Mode, MatState *C)
{
  static MatState *Fheap;
  static MatState *Aheap;
  MatState *Cmark, *Pmark;
  static int alloc, free;
  if (!Fheap || Fheap->Hp==NULL)
    {
      int c=0;
      int extension=1000;
      if (!Fheap){Fheap=vcalloc (1, sizeof (MatState));Fheap->free=1;free++;}
      if (!Aheap)Aheap=vcalloc (1, sizeof (MatState));
      while ( c!=extension)
	{
	  C=vcalloc ( 1, sizeof (MatState));
	  C->free=1;Fheap->Hn=C;C->Hp=Fheap;
	  Fheap=C;
	  c++;
	  free++;
	}
    }
  
  if ( Mode==DECLARE)
    {
     
      C=Fheap;
      Fheap=Fheap->Hp;
      C->Hn=C->Hp=NULL;
      if ( Aheap){Aheap->Hn=C;C->Hp=Aheap;Aheap=C;}
      else Aheap=C;
      alloc++;
      free--;
      C->free=0;
      C=CopyMatState(NULL, C);
      return C;
    }
  else if ( Mode==FREE)
    {
      if ( !C || C->free==1);
      else
	{
	  C=CopyMatState(NULL, C);
	  C->free=1;
	  if (C->Hp==NULL && C==Aheap)crash ("");
	  if (C==Aheap)Aheap=C->Hp;
	  if (C->Hn){(C->Hn)->Hp=C->Hp;}
	  if (C->Hp){(C->Hp)->Hn=C->Hn;}
	  C->Hp=C->Hn=NULL;
	  Fheap->Hn=C;C->Hp=Fheap;
	  Fheap=C;
	  alloc--;
	  free++;
	}
      return NULL;
    }
  else if ( Mode==FREE_ALL)
    {
      while ( Aheap)
	{
	  C=Aheap->Hp;
	  vfree (Aheap);
	  Aheap=C;
	}
      while ( Fheap)
	{
	  C=Fheap->Hp;
	  vfree (Fheap);
	  Fheap=C;
	}
    }
  else if ( Mode==INFO)
    {
      fprintf ( stderr, "\nAllocated: %d Free %d", alloc, free);
    }
  else if ( Mode==MARK)
    {
      
      if (C==NULL);
      else {C->Mn=Aheap;Aheap->Mp=C;}
      
      return Aheap;
    }
  else if ( Mode==UNMARK)
    {
      Pmark=Cmark=NULL;
    }
  else if ( Mode == FREE_MARK)
    {
      Cmark=C;
      Pmark=C->Mp;
      
      if ( Cmark==Pmark)return NULL;
      else if ( Cmark==Aheap)
	{Aheap=Pmark;C=Pmark->Hn;Pmark->Hn=NULL;}
      else
	{
	  (Cmark->Hn)->Hp=Pmark;
	  C=Pmark->Hn;
	  Pmark->Hn=Cmark->Hn;
	}
      
      Fheap->Hn=C;
      C->Hp=Fheap;
      Fheap=Cmark;
      Fheap->Hn=NULL;
  
      C=Fheap;
      while (C && !C->free)
	{
	  free++;alloc--;
	  C->free=1;
	  C=C->Hp;
	}
        
    }
  return NULL;
}


MatState* CopyMatState ( MatState*I, MatState*O)
{
  if (O==NULL || O->free==1) O=ManageMatState(DECLARE, NULL);
  if (I==NULL || I->free==1)I=NULL;
  O->i  =(I)?I->i:0;
  O->j  =(I)?I->j:0;
  O->st =(I)?I->st:FORBIDEN;
  O->pst=(I)?I->pst:FORBIDEN;
  O->sc =(I)?I->sc:FORBIDEN;
  O->n  =(I)?I->n:NULL;
  O->p  =(I)?I->p:NULL;
  O->m  =(I)?I->m:NULL;
  O->s  =(I)?I->m:NULL;
  
  return O;
}

/*********************************************************************************/
/*                                                                               */
/*                                                                               */
/*                      Comparisons                                              */
/*                                                                               */
/*                                                                               */
/*********************************************************************************/
int MaxDeltaMatState (MatState*S, MatState*E)
{
  if ( !S || !E) return -1;
  else return MAX((E->i-S->i),(E->j-S->j));
}
int MinDeltaMatState (MatState*S, MatState*E)
{
  if ( !S || !E) return -1;
  else return MIN((E->i-S->i),(E->j-S->j));
}
int MatStateAreIdentical (MatState*I, MatState*O)
{
  if ( !I || !O)return 0;

  if ( I->i!=O->i)return 0;
  if ( I->j!=O->j)return 0;
  if ( I->st!=O->st)return 0;
  return 1;
}
  


Hmm* define_probcons_model(Constraint_list *CL)
{
  Hmm *H;
  double gop=-10;
  double gep=-1;
  double lgop=-100;
  double lgep=-100;
  double freeT=0;
  int n=0;
  HmmState *S;
  
  
  H=declare_hmm(7);
  H->freeT=freeT=0;
 
  H->forbiden=FORBIDEN;
  H->start=START_STATE;
  H->end=END_STATE;
 
  /*define START*/
  S=H->S[n];  
  sprintf (S->name, "START"); S->state=n;
  
  S->DI=0;
  S->DJ=0;
  S->em=freeT;
 
  sprintf ( (S->T[S->nT])->name, "MATCH") ;(S->T[S->nT])->tr=freeT  ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=freeT  ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=freeT  ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "END")   ;(S->T[S->nT])->tr=freeT  ;S->nT++;
  
  n++;
  /*define END*/
  S=H->S[n];  
  sprintf (S->name, "END"); S->state=n;
  S->DI=0;
  S->DJ=0;
  S->em=freeT;
  n++;

  /*define Match*/
  S=H->S[n];  
  sprintf (S->name, "MATCH"); S->state=n;
  S->DI=1;
  S->DJ=1;
  S->em=H->forbiden;
  S->em_func=CL->get_dp_cost;
  
  sprintf ( (S->T[S->nT])->name, "MATCH") ;(S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=gop  ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=gop  ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "END");   (S->T[S->nT])->tr=freeT;S->nT++;
  
  n++;
  
  /*define Insert*/
  S=H->S[n];  
  sprintf (S->name, "INSERT"); S->state=n;
  S->DI=1;
  S->DJ=0;
  S->em=gep;
  sprintf ( (S->T[S->nT])->name, "MATCH") ; (S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "INSERT"); (S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "LINSERT");(S->T[S->nT])->tr=lgop ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "END");    (S->T[S->nT])->tr=-gop;S->nT++;
  
  n++;
  
  /*define LInsert*/
  S=H->S[n];  
  sprintf (S->name, "LINSERT"); S->state=n;
  S->DI=1;
  S->DJ=0;
  S->em=lgep;
  
  sprintf ( (S->T[S->nT])->name, "INSERT")  ;(S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "LINSERT");(S->T[S->nT])->tr=freeT;S->nT++;
   
  n++;
  
  
  /*define Delete*/
  S=H->S[n];  
  sprintf (S->name, "DELETE"); S->state=n;
  S->DI=0;
  S->DJ=1;
  S->em=gep;
  
  sprintf ( (S->T[S->nT])->name, "MATCH")   ;(S->T[S->nT])->tr=freeT;S->nT++;
   sprintf ( (S->T[S->nT])->name, "DELETE")  ;(S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "LDELETE") ;(S->T[S->nT])->tr=lgop ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "END")     ;(S->T[S->nT])->tr=-gop;S->nT++;
  
  n++;
  
  /*define LDelete*/
  S=H->S[n];  
  sprintf (S->name, "LDELETE"); S->state=n;
  S->DI=0;
  S->DJ=1;
  S->em=lgep;
  sprintf ( (S->T[S->nT])->name, "DELETE")  ;(S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "LDELETE");(S->T[S->nT])->tr=freeT;S->nT++;
    
  n++;
  
  
  H=bound_hmm ( H);
  return H;
}

Hmm* define_mnm_model(Constraint_list *CL)
{
  Hmm *H;
  double gop=20;



  double freeT=0;
  int n=0;
  HmmState *S;
  
  
  H=declare_hmm(6);
  H->freeT=freeT=0;
 
  H->forbiden=FORBIDEN;
  H->start=START_STATE;
  H->end=END_STATE;
 
  /*define START*/
  S=H->S[n];  
  sprintf (S->name, "START"); S->state=n;
  
  S->DI=0;
  S->DJ=0;
  S->em=freeT;
 
  sprintf ( (S->T[S->nT])->name, "MATCH") ;(S->T[S->nT])->tr=freeT  ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=freeT  ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=freeT  ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "NOMATCH");(S->T[S->nT])->tr=freeT  ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "END")   ;(S->T[S->nT])->tr=freeT  ;S->nT++;
  
  n++;
  /*define END*/
  S=H->S[n];  
  sprintf (S->name, "END"); S->state=n;
  S->DI=0;
  S->DJ=0;
  S->em=freeT;
  n++;

  /*define Match*/
  S=H->S[n];  
  sprintf (S->name, "MATCH"); S->state=n;
  S->DI=1;
  S->DJ=1;
  S->em=H->forbiden;
  S->em_func=CL->get_dp_cost;
  
  sprintf ( (S->T[S->nT])->name, "MATCH")  ;(S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "NOMATCH");(S->T[S->nT])->tr=gop  ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "END");   (S->T[S->nT])->tr=freeT;S->nT++;
  
  n++;
  
  /*define NOMatch*/
  S=H->S[n];  
  sprintf (S->name, "NOMATCH"); S->state=n;
  S->DI=1;
  S->DJ=1;
  S->em=freeT;
  S->em_func=NULL;
  
  sprintf ( (S->T[S->nT])->name, "NOMATCH")  ;(S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "MATCH")  ;(S->T[S->nT])->tr=freeT;S->nT++;  
  sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=freeT  ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=freeT  ;S->nT++;  
  sprintf ( (S->T[S->nT])->name, "END");   (S->T[S->nT])->tr=freeT;S->nT++;
  
  n++;
  /*define Insert*/
  S=H->S[n];  
  sprintf (S->name, "INSERT"); S->state=n;
  S->DI=1;
  S->DJ=0;
  S->em=freeT;
  sprintf ( (S->T[S->nT])->name, "NOMATCH") ; (S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "INSERT")  ; (S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "END");    (S->T[S->nT])->tr=freeT;S->nT++;
  
  n++;
  
  /*define Delete*/
  S=H->S[n];  
  sprintf (S->name, "DELETE"); S->state=n;
  S->DI=0;
  S->DJ=1;
  S->em=freeT;
  
  sprintf ( (S->T[S->nT])->name, "NOMATCH")   ;(S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "DELETE")  ;(S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "END")     ;(S->T[S->nT])->tr=freeT;S->nT++;
  
  n++;
  
  
  H=bound_hmm ( H);
  return H;
}

Hmm* define_simple_model(Constraint_list *CL)
{
  Hmm *H;
  double gop=-10;
  double gep=-1;
  double freeT=0;
  int n=0;
  HmmState *S;
  
  
  H=declare_hmm(5);
  H->freeT=freeT=0;
 
  H->forbiden=FORBIDEN;
  H->start=START_STATE;
  H->end=END_STATE;
 
  /*define START*/
  S=H->S[n];  
  sprintf (S->name, "START"); S->state=n;
  
  S->DI=0;
  S->DJ=0;
  S->em=freeT;
 
  sprintf ( (S->T[S->nT])->name, "MATCH") ;(S->T[S->nT])->tr=freeT  ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=freeT  ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=freeT  ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "END")   ;(S->T[S->nT])->tr=freeT  ;S->nT++;
  
  n++;
  /*define END*/
  S=H->S[n];  
  sprintf (S->name, "END"); S->state=n;
  S->DI=0;
  S->DJ=0;
  S->em=freeT;
  n++;

  /*define Match*/
  S=H->S[n];  
  sprintf (S->name, "MATCH"); S->state=n;
  S->DI=1;
  S->DJ=1;
  S->em=H->forbiden;
  S->em_func=CL->get_dp_cost;
  
  sprintf ( (S->T[S->nT])->name, "MATCH") ;(S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=gop  ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=gop  ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "END");   (S->T[S->nT])->tr=freeT;S->nT++;
  
  n++;
  
  /*define Insert*/
  S=H->S[n];  
  sprintf (S->name, "INSERT"); S->state=n;
  S->DI=1;
  S->DJ=0;
  S->em=gep;
  sprintf ( (S->T[S->nT])->name, "MATCH") ; (S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "INSERT"); (S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "DELETE"); (S->T[S->nT])->tr=freeT;S->nT++;
  
  sprintf ( (S->T[S->nT])->name, "END");    (S->T[S->nT])->tr=-gop;S->nT++;
  
  n++;
  
    
  /*define Delete*/
  S=H->S[n];  
  sprintf (S->name, "DELETE"); S->state=n;
  S->DI=0;
  S->DJ=1;
  S->em=gep;
  
  sprintf ( (S->T[S->nT])->name, "MATCH")   ;(S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "DELETE")  ;(S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "INSERT"); (S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "END")     ;(S->T[S->nT])->tr=-gop;S->nT++;
  
  n++;
  
  H=bound_hmm ( H);
  return H;
}

Hmm* define_two_mat_model(Constraint_list *CL)
{
  Hmm *H;
  double gop=-15;
  double gep=-2;
  double lgop=-6;
  double lgep=-1;
  double freeT=0;
  int n=0;
  HmmState *S;
  
  
  H=declare_hmm(8);
  H->freeT=freeT=0;
 
  H->forbiden=FORBIDEN;
  H->start=START_STATE;
  H->end=END_STATE;
 
  /*define START*/
  S=H->S[n];  
  sprintf (S->name, "START"); S->state=n;
  
  S->DI=0;
  S->DJ=0;
  S->em=freeT;
 
  sprintf ( (S->T[S->nT])->name, "MATCH1") ;(S->T[S->nT])->tr=freeT  ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "MATCH2") ;(S->T[S->nT])->tr=freeT  ;S->nT++;
  
  sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=freeT  ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=freeT  ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "END")   ;(S->T[S->nT])->tr=freeT  ;S->nT++;
  
  n++;
  /*define END*/
  S=H->S[n];  
  sprintf (S->name, "END"); S->state=n;
  S->DI=0;
  S->DJ=0;
  S->em=freeT;
  n++;

  /*define Match*/
  S=H->S[n];  
  sprintf (S->name, "MATCH1"); S->state=n;
  S->DI=1;
  S->DJ=1;
  S->em=H->forbiden;
  S->em_func=get_dp_cost_pam_matrix;
  
  sprintf ( (S->T[S->nT])->name, "MATCH1") ;(S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=gop  ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=gop  ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "END");   (S->T[S->nT])->tr=freeT;S->nT++;
  
  n++;
  
  /*define Match*/
  S=H->S[n];  
  sprintf (S->name, "MATCH2"); S->state=n;
  S->DI=1;
  S->DJ=1;
  S->em=H->forbiden;
  S->em_func=get_dp_cost_blosum_matrix;
  
  sprintf ( (S->T[S->nT])->name, "MATCH2") ;(S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=gop  ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=gop  ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "END");   (S->T[S->nT])->tr=freeT;S->nT++;
  
  n++;

  /*define Insert*/
  S=H->S[n];  
  sprintf (S->name, "INSERT"); S->state=n;
  S->DI=1;
  S->DJ=0;
  S->em=gep;
  sprintf ( (S->T[S->nT])->name, "MATCH2") ; (S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "MATCH1") ; (S->T[S->nT])->tr=freeT;S->nT++;
  
  sprintf ( (S->T[S->nT])->name, "INSERT"); (S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "LINSERT");(S->T[S->nT])->tr=lgop ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "END");    (S->T[S->nT])->tr=-gop;S->nT++;
  
  n++;
  
  /*define LInsert*/
  S=H->S[n];  
  sprintf (S->name, "LINSERT"); S->state=n;
  S->DI=1;
  S->DJ=0;
  S->em=lgep;
  
  sprintf ( (S->T[S->nT])->name, "INSERT")  ;(S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "LINSERT");(S->T[S->nT])->tr=freeT;S->nT++;
    
  n++;
  
  
  /*define Delete*/
  S=H->S[n];  
  sprintf (S->name, "DELETE"); S->state=n;
  S->DI=0;
  S->DJ=1;
  S->em=gep;
  
  sprintf ( (S->T[S->nT])->name, "MATCH2")   ;(S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "MATCH1")   ;(S->T[S->nT])->tr=freeT;S->nT++;
  
  sprintf ( (S->T[S->nT])->name, "DELETE")  ;(S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "LDELETE") ;(S->T[S->nT])->tr=lgop ;S->nT++;
  sprintf ( (S->T[S->nT])->name, "END")     ;(S->T[S->nT])->tr=-gop;S->nT++;
  n++;
  
  /*define LDelete*/
  S=H->S[n];  
  sprintf (S->name, "LDELETE"); S->state=n;
  S->DI=0;
  S->DJ=1;
  S->em=lgep;
  sprintf ( (S->T[S->nT])->name, "DELETE")  ;(S->T[S->nT])->tr=freeT;S->nT++;
  sprintf ( (S->T[S->nT])->name, "LDELETE");(S->T[S->nT])->tr=freeT;S->nT++;
  n++;
  
  if ( n!=H->nS)
    {
      fprintf ( stderr, "\nERROR in HMM definition [FATAL:%s]", PROGRAM);
      myexit (EXIT_FAILURE);
    }
  
  H=bound_hmm ( H);
  return H;
} 
void DisplayMatState ( MatState *S, char *s)
{
  if ( S==NULL)fprintf ( stderr, "%s: Cell is undefined", s);
  else  fprintf ( stderr, "%s: i=%d j=%d st=%d pst=%d sc=%d Free %d", s, S->i, S->j, S->st, S->pst, (int)S->sc, S->free);
}
void testfunc ( MatState *S, char *s)
{
  if ( S==NULL)return;
  fprintf ( stderr, "\n#### %s ", s);
  while ( S){DisplayMatState ( S,"\n\t");S=S->n;}
  fprintf ( stderr, "\n");
}
  
































#ifdef BACKHERE
	  
	  if ( i>0 && j>0)
	  m=emit_pair_default[alphabetDefault[seq1[a]]][alphabetDefault[seq2[a]]];
	  /*Match*/
	  F[M][i][j]=F[M][i-1][j-1];
	  
	  
	
	M[Match][i][j]=m+log_add3( M[Match][i-step_i][j-step_j],M[I][i-step_i][j],M[D][i][j-step_j]);		
	M[D    ][i][j]=log_add3(gep,M[Match][i       ][j-step_j]+gop,M[D][i       ][j-step_j]);
	M[I    ][i][j]=log_add3(gep,M[Match][i-step_i][j       ]+gop,M[I][i-step_i][j       ]);
	
	/*Long gaps
	M[Match][i][j]=log_add3(M[Match][i][j], M[LI][i-step_i][j],M[LD][i][j-step_j]);
	M[LI    ][i][j]=log_add3(lgep, M[I][i-step_i][j       ]+lgop,M[LI][i-step_i][j      ]);	
	M[LD    ][i][j]=log_add3(lgep, M[D][i       ][j-step_j]+lgop,M[LD][i       ][j-step_j]);
	*/
	
	}
    }
  retun M;
MatState* RviterbiL_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL,MatState *S, MatState *E)
{
  MatState *Mid;
  Mid=viterbiL_hmm (A,ns, ls,H, CL, S, E);
  
  if (!Mid)
    {
      return S;
    }
  else if ( Mid->n)
    {
      return Mid;
    }
      
  else
    {
      Mid->p=S;S->n=Mid;
      Mid->n=E;E->p=Mid;
      RviterbiL_hmm (A,ns, ls,H, CL,S,   Mid);
      RviterbiL_hmm (A,ns, ls,H, CL,Mid, E);
      return S;
    }
}

MatState* viterbiL_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL, MatState *S,MatState *E)
{
  int current,memory, dim;
  double e, v,t;
  int i,j,pi,pj, s, k;
  int start_i, start_j, end_i, end_j, l1, l2;
  HmmState *S1, *S2;
  static MatState ****M;
  static int maxl;
  MatState *Mid=NULL;

  MatState *CC, *PCC;
  int midpoint;
  int Delta;
  

  if ( MatStateAreIdentical ( S, E))return NULL;
  
  l1=strlen (A->seq_al[ls[0][0]]);
  l2=strlen (A->seq_al[ls[1][0]]);
    
  midpoint=S->i+(E->i-S->i)/2;
  Delta=E->i-S->i;

  start_i=S->i;end_i=E->i;
  start_j=S->j;end_j=E->j;
    
  dim=H->order+2;current=0;memory=H->order+1;
  if (!M || (l2+1)>maxl)
    { free_arrayN((void **)M, 4);
      M=declare_arrayN(4, sizeof ( MatState), dim, maxl=(l2+1), H->nS,1);
    }
      
  /*MAKE THE VITERBI FROM S(tart) to E(nd)*/
  for ( i=start_i; i<=end_i; i++)
    {
      M= (MatState****)recycle  ( (void **)M,H->order+1,1);
      for ( j=start_j; j<=end_j; j++)
	{
	  for ( s=H->nS-1;s>=0; s--)
	    {

	      S1=H->S[s];
	      pi=i-S1->DI;pj=j-S1->DJ;
	      CC=M[current][j][s];
	      CC->i=i; CC->j=j; CC->st=s;CC->sc=H->forbiden;CC->p=CC->n=CC->m=NULL;CC->sc=H->forbiden; 
	      if (i==start_i && j==start_j && s==S->st)  {CopyMatState(S,CC);}
	      else if ( i==end_i && j==end_j && s==E->st && s!=H->end)
		{
		  S2=H->S[E->pst];
		  CopyMatState(E,CC);
		  CC->p=M[S1->DI][j-S1->DJ][S2->state];
		}
	      else if ( pi<start_i || pj<start_j)        {CC->sc=H->forbiden;}
	      else 
		{
		  for (k=1; k<=H->fromM[S1->state][0]; k++)
		    {
		      S2=H->S[H->fromM[s][k]];
		      PCC=M[S1->DI][j-S1->DJ][S2->state];
	
		      if      ( pi+pj!=0 && S2->state==H->start) {t=H->forbiden;}
		      else if ( !(pi==l1 && pj==l2) && s==H->end){t=H->forbiden;}
		      else   t=H->T[S2->state][S1->state];

		      v=hmm_add(t,PCC->sc);
		      if ( v!=H->forbiden && (CC->sc==H->forbiden || v> CC->sc)){CC->sc=v; CC->pst=S2->state;CC->p=PCC;}
		    }
		  
		  e=(S1->em==H->forbiden)?S1->em_func (A, A->pos, ns[0], ls[0],i-1, A->pos,ns[1], ls[1], j-1, CL):S1->em;
		  CC->sc=hmm_add(CC->sc,e);
		}
	      
	      if (i==midpoint)CC->m=CopyMatState(CC, M[memory][j][s]);
	      else if (i>midpoint && CC->sc!=H->forbiden) CC->m=(M[S1->DI][j-S1->DJ][CC->pst])->m;
	    }
	}
    }  
  
  if ( E->st==H->end)CopyMatState ((M[current][end_j][E->st]),E);
  
  if ( Delta>1)
    {
      Mid=CopyMatState ((M[current][end_j][E->st])->m,NULL);
    }
  else if ( Delta==1)
    { 
      CC=M[current][E->j][E->st];
      Mid=E;
      while (!MatStateAreIdentical (CC->p, S) )
	{
	  Mid->p=CopyMatState(CC->p,NULL);
	  (Mid->p)->n=Mid;
	  Mid=Mid->p;CC=CC->p;
	}
      Mid->p=S;
      S->n=Mid;
      Mid=S;
    }

  return Mid;
}
#endif
/*********************************COPYRIGHT NOTICE**********************************/
/*© Centre National de la Recherche Scientifique (CNRS) */
/*and */
/*Please Cite: Notredame*/
/*Mon May 17 20:15:35 MDT 2004. */
/*All rights reserved.*/
/*NOTICE:                                                                                                                                     |*/
/*  This file is an integral part of the */
/*  ALIGN_TWO_SEQ Software. */
/*  Its content is protected and all */
/*  the conditions mentioned in the licensing */
/*  agreement of the software apply to this file.*/
/*...............................................                                                                                      |*/
/*  If you need some more information, or if you */
/*  wish to obtain a full license, please contact: */
/*  cedric.notredame@europe.com*/
/*...............................................                                                                                                                                     |*/
/**/
/**/
/*	*/
/*********************************COPYRIGHT NOTICE**********************************/
/*********************************COPYRIGHT NOTICE**********************************/
/*© Centro de Regulacio Genomica */
/*and */
/*Cedric Notredame */
/*Tue Oct 27 10:12:26 WEST 2009. */
/*All rights reserved.*/
/*This file is part of T-COFFEE.*/
/**/
/*    T-COFFEE is free software; you can redistribute it and/or modify*/
/*    it under the terms of the GNU General Public License as published by*/
/*    the Free Software Foundation; either version 2 of the License, or*/
/*    (at your option) any later version.*/
/**/
/*    T-COFFEE is distributed in the hope that it will be useful,*/
/*    but WITHOUT ANY WARRANTY; without even the implied warranty of*/
/*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the*/
/*    GNU General Public License for more details.*/
/**/
/*    You should have received a copy of the GNU General Public License*/
/*    along with Foobar; if not, write to the Free Software*/
/*    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*/
/*...............................................                                                                                      |*/
/*  If you need some more information*/
/*  cedric.notredame@europe.com*/
/*...............................................                                                                                                                                     |*/
/**/
/**/
/*	*/
/*********************************COPYRIGHT NOTICE**********************************/