7 #include "io_lib_header.h"
8 #include "util_lib_header.h"
9 #include "define_header.h"
10 #include "dp_lib_header.h"
12 #define hmm_add(x,y) ((x==UNDEFINED || y==UNDEFINED)?UNDEFINED:(x+y))
13 #define MAX_EMISSION 256
15 /*********************************************************************************/
21 /*********************************************************************************/
22 char alphabetDefault[] = "ARNDCQEGHILKMFPSTWYV";
23 double emitPairsDefault[20][20] = {
24 {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},
25 {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},
26 {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},
27 {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},
28 {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},
29 {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},
30 {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},
31 {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},
32 {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},
33 {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},
34 {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},
35 {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},
36 {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},
37 {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},
38 {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},
39 {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},
40 {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},
41 {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},
42 {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},
43 {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}
46 static void DisplayMatState ( MatState *S, char *s);
53 void check_viterbiL ( Alignment *A,int *ns, int **ls, Constraint_list *CL);
54 MatState *viterbi2path2 ( double ***Sc, int ***St, Hmm *H, MatState *S, MatState *E);
55 void testfunc ( MatState *S, char *s);
58 /*********************************************************************************/
64 /*********************************************************************************/
65 Alignment * analyze_alignment ( Alignment *A)
67 evaluate_alignment (A);
68 H=define_msa_model (-100);
69 M=seq_viterbi_hmm (A->seq_al[0], H);
70 path=seq_viterbi2path ( seq, H, M);
74 Hmm* define_msa_model(double penalty)
91 sprintf (S->name, "START"); S->state=n;
97 sprintf ( (S->T[S->nT])->name, "C") ;(S->T[S->nT])->tr=freeT ;S->nT++;
98 sprintf ( (S->T[S->nT])->name, "W");(S->T[S->nT])->tr=freeT ;S->nT++;
102 sprintf (S->name, "END"); S->state=n;
110 sprintf (S->name, "C"); S->state=n;
114 S->em_func=em_correct_msa;
116 sprintf ( (S->T[S->nT])->name, "C") ;(S->T[S->nT])->tr=freeT;S->nT++;
117 sprintf ( (S->T[S->nT])->name, "W");(S->T[S->nT])->tr=penalty ;S->nT++;
118 sprintf ( (S->T[S->nT])->name, "END"); (S->T[S->nT])->tr=freeT;S->nT++;
123 sprintf (S->name, "INSERT"); S->state=n;
127 S->em_func=em_wrong_msa;
128 sprintf ( (S->T[S->nT])->name, "C") ; (S->T[S->nT])->tr=penalty;S->nT++;
129 sprintf ( (S->T[S->nT])->name, "W"); (S->T[S->nT])->tr=freeT;S->nT++;
130 sprintf ( (S->T[S->nT])->name, "END"); (S->T[S->nT])->tr=-gop;S->nT++;
135 sprintf (S->name, "LINSERT"); S->state=n;
140 sprintf ( (S->T[S->nT])->name, "INSERT") ;(S->T[S->nT])->tr=freeT;S->nT++;
141 sprintf ( (S->T[S->nT])->name, "LINSERT");(S->T[S->nT])->tr=freeT;S->nT++;
148 /*********************************************************************************/
151 /* simple HMM: Viterbi */
154 /*********************************************************************************/
155 double pavie_em_func (Hmm*H, HmmState *S, int v);
156 Hmm* define_full_model(int nstate,char **state_list, char *model_name,Generic_em_func evaluation_func );
157 char **produce_state_name (int nstate,char **list, char *model_name, Hmm* H);
158 double** seq_viterbi_hmm (char *seq, Hmm *H);
159 int * seq_viterbi2path (char *s, Hmm *H, double **M);
160 double analyze_sequence ( char *seq, Hmm*H);
162 double pavie_emission (Hmm*H, HmmState *S, int v)
169 if ( v==n[0] || ( v=='*' && n[0]=='E')) return H->freeT;
172 Hmm* define_full_model(int nstate, char **list, char *model_name, Generic_em_func emission_function)
174 /*list: a list of the state names: does not include START or END*/
175 /*model_name: a string that will be appended to the names in list*/
182 H=declare_hmm(nstate+2);
185 H->forbiden=FORBIDEN;
186 H->start=START_STATE;
189 list=produce_state_name (nstate,list, model_name, H);
192 for (n=0; n<nstate; n++)
197 sprintf ( S->name, "%s", list[n]);
198 S->em_func2=emission_function;
199 if (n==H->end || n==H->start){S->DI=0;S->DJ=0;}
200 else S->DI=1;S->DJ=0;
203 S->em_func2=emission_function;
204 for (a=0; a< MAX_EMISSION; a++)S->em2[a]=H->freeT;
206 for (a=0; a<nstate && n!=H->end; a++)
208 if (a!=H->start && !(n==H->start && a==H->end) )
210 sprintf ( (S->T[S->nT])->name, "%s", list[a]);
211 (S->T[S->nT])->tr=H->freeT;
219 char **produce_state_name (int nstate,char **list, char *model_name, Hmm* H)
225 new_list=declare_char ( nstate, 100);
226 for ( a=0, b=0, c=0; a< nstate; a++)
228 if ( a==H->start)sprintf ( new_list[a], "START");
229 else if ( a==H->end)sprintf ( new_list[a], "END");
230 else if ( list==NULL){sprintf ( new_list[a], "%c%s", 'a'+b, (model_name)?model_name:"");b++;}
231 else {sprintf ( new_list[a], "%c%s", list[a][c], (model_name)?model_name:"");c++;}
236 int seq_viterbi_pair_wise (Alignment *A,int*ns, int **ls,Constraint_list *CL)
239 analyze_sequence (A->seq_al[0], NULL);
240 myexit (EXIT_FAILURE);
243 double analyze_sequence ( char *seq, Hmm *H)
251 H=define_full_model(5, NULL,"_first", pavie_emission);
255 M=seq_viterbi_hmm (seq, H);
256 path=seq_viterbi2path (seq, H, M);
257 return M[H->end][strlen (seq)];
261 double** seq_viterbi_hmm (char *seq, Hmm *H)
263 /*Given a model H and a sequence seq*/
266 int i,pi, bestk, s, k, l1;
271 M=declare_double (H->nS*2,l1+2);
275 for ( i=0; i<=l1; i++)
277 for ( s=0; s< H->nS; s++)
283 if ( pi<0){M[s][i]=H->forbiden;}/*Boundary*/
286 if (pi==0) {max=H->T[(int)H->start][s];bestk=H->start;}/*Start*/
289 for (k=1; k<=H->fromM[S1->state][0]; k++)
291 S2=H->S[H->fromM[s][k]];
292 if ( S2->state==H->start || S2->state==H->end)continue;
293 v=hmm_add((M[S2->state][pi]),(H->T[S2->state][S1->state]));
294 if ( v!=H->forbiden && (max==H->forbiden || v>max)){max=v;bestk=S2->state;}
297 if (S1->em2)e=S1->em2[(int)seq[pi]];
298 else e=S1->em_func2(H,S1, (int)seq[pi]);
307 /*Terminate viterbi: connect the path to the END state*/
310 for (k=0; k< H->nS; k++)
312 if (k==H->start || k==H->end);
315 v=(M[k][l1]==H->forbiden || H->T[k][H->end]==H->forbiden)?H->forbiden:M[k][l1]+H->T[k][H->end];
316 if ( max==H->forbiden || v>max){bestk=k;max=v;}
320 M[H->nS+H->end][l1]=bestk;
324 int * seq_viterbi2path (char *s, Hmm *H, double **M)
332 path=vcalloc (l1+1, sizeof (int));
335 cs=M[H->nS+H->end][i];
345 /*fprintf ( stderr, "%d", cs);*/
347 invert_list_int (path, l);
348 path[l++]=H->forbiden;
353 /*********************************************************************************/
356 /* pairHMM: Viterbi */
359 /*********************************************************************************/
360 Hmm* define_mnm_model(Constraint_list *CL);
361 int viterbi_pair_wise_OLD (Alignment *A,int*ns, int **ls,Constraint_list *CL)
368 A->pos=aln2pos_simple( A, -1, ns, ls);
370 // H=define_mnm_model (CL);
371 H=define_two_mat_model (CL);
373 l1=strlen (A->seq_al[ls[0][0]]);
374 l2=strlen (A->seq_al[ls[1][0]]);
375 M=viterbi_hmm (A, ns, ls, H, CL);
376 path=viterbi2path (l1,l2, H,M);
377 A=viterbipath2aln (A,ns,ls,path, H);
378 A->score=A->score_aln=M[H->end][l1][l2];
379 for ( a=0; a< H->nS*2; a++)free_double (M[a], -1);
381 free_int (A->pos, -1);
390 Alignment * viterbipath2aln (Alignment *A, int *ns,int **ls,int *tb, Hmm *H)
394 int a, b, c, len, ch;
398 len=0;while (tb[len]!=H->forbiden)len++;
400 if ( A->declared_len<=len)A=realloc_aln2 ( A,A->max_n_seq,2*len);
403 char_buf=vcalloc (len+1, sizeof (char));
404 l[0]=strlen ( A->seq_al[ls[0][0]]);
405 l[1]=strlen ( A->seq_al[ls[1][0]]);
407 for ( c=0; c< 2; c++)
408 for ( a=0; a< ns[c]; a++)
410 for (ch=0, b=0; b<len; b++)
413 if ( (c==0 && S->DI)|| (c==1 && S->DJ) )
414 char_buf[b]=aln[ls[c][a]][ch++];
419 sprintf (aln[ls[c][a]],"%s", char_buf);
420 if ( l[c]!=ch){fprintf (stderr, "\nERROR: Wrong Size Of Alignmnent (Real %d, Observed %d)[FATAL:%s]",l[c], ch, PROGRAM);}
429 double*** viterbi_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL)
433 int a, i,pi, bestk,j,pj, s, k, l1, l2;
436 l1=strlen (A->seq_al[ls[0][0]]);
437 l2=strlen (A->seq_al[ls[1][0]]);
439 M=vcalloc (H->nS*2, sizeof (double**));
440 for ( a=0; a<H->nS*2; a++)M[a]=declare_double (l1+2, l2+2);
445 for ( i=0; i<=l1; i++)
446 for ( j=0; j<=l2; j++)
448 for ( s=0; s< H->nS; s++)
455 if ( pi<0 ||pj<0){M[s][i][j]=H->forbiden;}/*Boundary*/
458 if (pi+pj==0) {max=H->T[H->start][s];bestk=H->start;}/*Start*/
461 for (k=1; k<=H->fromM[S1->state][0]; k++)
463 S2=H->S[H->fromM[s][k]];
464 if ( S2->state==H->start || S2->state==H->end)continue;
465 v=(M[S2->state][pi][pj]==H->forbiden)?H->forbiden:(M[S2->state][pi][pj]+H->T[S2->state][S1->state]);
466 if ( v!=H->forbiden && (max==H->forbiden || v>max)){max=v;bestk=S2->state;}
470 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;
471 e=(max==H->forbiden || e==H->forbiden)?H->forbiden:e+max;
474 M[s+H->nS][i][j]=bestk;
479 /*Terminate viterbi: connect the path to the END state*/
482 for (k=0; k< H->nS; k++)
484 if (k==H->start || k==H->end);
487 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];
488 if ( max==H->forbiden || v>max){bestk=k;max=v;}
491 M[H->end][l1][l2]=max;
492 M[H->nS+H->end][l1][l2]=bestk;
497 /*********************************************************************************/
500 /* HMM: Decode/Traceback */
503 /*********************************************************************************/
504 int * traceback (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL,MatState *S, MatState *E, int **seg_list)
512 l1=strlen (A->seq_al[ls[0][0]]);
513 l2=strlen (A->seq_al[ls[1][0]]);
514 path=vcalloc ( l1+l2+1, sizeof (int));
516 while ( S->st!=H->end)
518 DisplayMatState (S, "\n\tTraceback");
520 if ( N && S && (((N->i-S->i)>1) ||((N->j-S->j)>1)))
522 RviterbiD_hmm (A,ns,ls,H,CL,S,N,seg_list);
527 ManageMatState (FREE,S);
535 int * viterbi2path (int l1,int l2, Hmm *H, double ***M)
543 path=vcalloc (l1+l2+1, sizeof (int));
546 cs=M[H->nS+H->end][i][j];
554 cs=M[H->nS+cs][i][j];
557 /*fprintf ( stderr, "%d", cs);*/
559 invert_list_int (path, l);
560 path[l++]=H->forbiden;
565 /*********************************************************************************/
568 /* HMM Viterbi Linear */
571 /*********************************************************************************/
574 int viterbiL_pair_wise (Alignment *A,int*ns, int **ls,Constraint_list *CL)
582 A->pos=aln2pos_simple( A, -1, ns, ls);
583 Start=ManageMatState ( DECLARE, NULL);
584 End=ManageMatState ( DECLARE, NULL);
585 H=define_simple_model (CL);
586 l1=strlen (A->seq_al[ls[0][0]]);
587 l2=strlen (A->seq_al[ls[1][0]]);
591 Start->i=0 ;Start->j=0 ; Start->st=H->start;Start->sc=0;
592 End->i =l1; End->j=l2; End ->st=H->end;
593 Start=RviterbiL_hmm (A, ns, ls, H, CL, Start,End);
594 path=traceback (A, ns, ls, H, CL, Start,NULL, NULL);
596 A=viterbipath2aln (A,ns,ls,path, H);
599 free_int (A->pos, -1);
606 MatState* RviterbiL_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL,MatState *S, MatState *E)
612 viterbiL_hmm (A,ns, ls,H, CL, S, E);
615 if ( S->n==E)return S;
616 if ( E->sc==H->forbiden)
619 fprintf ( stderr, "\nERROR: The Requested Model (Cf Above) Cannot Produce the Pair-Alignment\nYou must allow extra possible transitions\n[FATAL:%s]", PROGRAM);
620 myexit ( EXIT_FAILURE);
627 d1=MinDeltaMatState(S,E);
628 d2=MaxDeltaMatState(S,E);
629 align=((d1==1 && d2==1) || ( d1==0))?0:1;
630 if (align)RviterbiL_hmm (A,ns, ls,H, CL,S,E);
637 #define Dim(i,j) (i*LenJ+j)
638 MatState* viterbiL_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL, MatState *S,MatState *E)
640 int current, previous,row, prow;
642 int a,i,j,pi,pj, s, k;
643 int start_i, start_j, end_i, end_j, l1, l2;
645 MatState *CC, *PCC,*tS, *tE, *mark=NULL;
649 static MatState ***M;
652 static int LenJ, LenI;
653 int MaxDelta=50, DeltaI, DeltaJ;
657 DisplayMatState (S, "\n\tS");
658 DisplayMatState (E, "\n\tE");
663 for ( a=0; a<2; a++)memset(M[a],0,LenJ*LenI*sizeof (MatState*));
664 free_arrayN((void **)M, 3);M=NULL;
665 ManageMatState ( FREE_ALL, NULL);
670 if ( MatStateAreIdentical ( S, E))return NULL;
671 l1=strlen (A->seq_al[ls[0][0]]);l2=strlen (A->seq_al[ls[1][0]]);
673 midpoint=S->i+((E->i-S->i)/2);
677 start_i=S->i;end_i=E->i;start_j=S->j;end_j=E->j;
678 current=0;previous=1;
685 M=declare_arrayN(3, sizeof ( MatState),2,LenI*LenJ,0);
689 /*MAKE THE VITERBI FROM S(tart) to E(nd)*/
690 mark=ManageMatState ( MARK, mark);
691 for (i=start_i; i<=end_i; i++)
694 for ( j=start_j; j<=end_j; j++)
696 DeltaJ=((FABS(j-i))<MaxDelta)?1:MaxDelta+1; /*Keep Pointers on a DealtaMax Band around the main diagonal*/
697 for ( s=H->nS-1;s>=0; s--)
699 S1=H->S[s];pi=i-S1->DI;prow=S1->DI;pj=j-S1->DJ;
701 CC=M[row][Dim(j,s)]=CopyMatState(NULL, M[row][Dim(j,s)]);
702 CC->i=i; CC->j=j; CC->st=s;PCC=NULL;
704 if (i==start_i && j==start_j && s==S->st){CC=CopyMatState(S,CC);}
705 else if ( i==end_i && j==end_j && E->st!=H->end && s!=E->st)CC->sc=H->forbiden;
706 else if ( pi<start_i || pj<start_j) {CC->sc=H->forbiden;}
709 for (k=1; k<=H->fromM[S1->state][0]; k++)
711 S2=H->S[H->fromM[s][k]];
712 PCC=M[prow][Dim((j-S1->DJ),(S2->state))];
715 else if ( pi+pj!=0 && S2->state==H->start);
716 else if ( !(pi==l1 && pj==l2) && s==H->end);
720 v=hmm_add(CC->sc,H->T[PCC->st][CC->st]);
722 v=lu_RviterbiD_hmm(A,ns, ls, H, CL,PCC,CC, NULL);
723 if ( v!=H->forbiden && (CC->sc==H->forbiden || v> CC->sc)){CC->sc=v; CC->pst=S2->state;CC->p=PCC;}
727 if (CC->sc==H->forbiden);
728 else if (i==midpoint || DeltaI<=MaxDelta||DeltaJ<=MaxDelta ||(i==start_i && j==start_j && s==S->st) )
730 CC->m=(CC->p)?(CC->p)->m:NULL;
731 PCC=CopyMatState(CC,NULL);
732 PCC->m=CC->m;CC->m=PCC;
734 else CC->m=(CC->p)?(CC->p)->m:NULL;
738 for ( j=start_j; j<=end_j && i!=end_i; j++)
740 for ( s=H->nS-1;s>=0; s--)
743 CC=(M[prow][Dim(j,s)]);M[prow][Dim(j,s)]=M[row][Dim(j,s)];M[row][Dim(j,s)]=CC;
744 if (M[prow][Dim(j,s)]) M[row ][Dim(j,s)]=CopyMatState ( M[prow][Dim(j,s)], M[row][Dim(j,s)]);
751 mark=ManageMatState ( MARK,mark);
755 if ( E->st==H->end || E->st==H->forbiden){E=CopyMatState ((M[row][Dim(end_j,E->st)]),E);}
760 PCC=CopyMatState (M[row][Dim(end_j,E->st)], NULL);
762 if ( MatStateAreIdentical(PCC,PCC->m))PCC=PCC->m;
766 tS=CopyMatState (PCC->m,NULL); tS->n=PCC; PCC->p=tS;PCC=tS;
772 S->n=tS->n; (S->n)->p=S;
773 E->p=tE->p; (E->p)->n=E;
775 for ( a=0; a<2; a++)memset(M[a],0,LenJ*LenI*sizeof (MatState*));
776 ManageMatState ( FREE_MARK,mark);
779 while (S && S->p!=E){S->m=NULL;S=S->n;}/*Clean the memory of the rturned Cells*/
783 /*********************************************************************************/
786 /* HMM Viterbi Diagonals */
789 /*********************************************************************************/
790 int viterbiD_pair_wise (Alignment *A,int*ns, int **ls,Constraint_list *CL)
803 A->pos=aln2pos_simple( A, -1, ns, ls);
805 Start=ManageMatState ( DECLARE, NULL);
806 End=ManageMatState ( DECLARE, NULL);
807 H=define_simple_model (CL);
808 l1=strlen (A->seq_al[ls[0][0]]);
809 l2=strlen (A->seq_al[ls[1][0]]);
811 main_i=MAX(1,(l2-l1)+1);
812 main_j=MAX(1,(l1-l2)+1);
814 seg_list=declare_arrayN(2, sizeof (int), l1+l2+3, 3);
815 seg_list[0][0]=DIAGONALS;
819 for ( b=1,a=l1; a>= 1; a--)
821 if (a<50 || (b==main_i && a==main_j))
825 seg_list[c][2]=MIN((l1-a), (l2-b));
831 for ( b=2,a=1; b<= l2; b++, c++)
833 if (b<50 || (b==main_i && a==main_j))
837 seg_list[c][2]=MIN((l1-a), (l2-b));
842 seg_list[c][0]=FORBIDEN;
844 Start->i=0 ;Start->j=0 ; Start->st=H->start;Start->sc=0;
845 End->i =l1; End->j=l2; End ->st=H->end;
846 Start=RviterbiD_hmm (A, ns, ls, H, CL, Start,End,seg_list);
849 path=traceback (A, ns, ls, H, CL, Start,NULL, NULL);
853 A=viterbipath2aln (A,ns,ls,path, H);
855 viterbiD_hmm (NULL, ns, ls, H, CL, Start,End, seg_list);
857 free_int (A->pos, -1);
858 free_arrayN((void **)seg_list, 2);
865 double lu_RviterbiD_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL,MatState *S, MatState *E, int **seg_list)
869 static MatState *cS=NULL, *cE=NULL;
871 max=MAX((E->i-S->i), (E->j-S->j));
872 min=MIN((E->i-S->i), (E->j-S->j));
875 if ( S->sc==H->forbiden) return H->forbiden;
878 e=hmm_add(S->sc,H->T[S->st][E->st]);
879 if ( H->T[E->st][E->st]!=H->forbiden)e=hmm_add(e, (max-1)*H->T[E->st][E->st]);
880 if ( (H->S[E->st])->em!=H->forbiden) e=hmm_add(e, max *(H->S[E->st])->em );
883 else if ( min>0 && max>1)
886 fprintf ( stderr, "\nWarning: Disjoined Diagonals");
887 DisplayMatState (S, "\n\tS");
888 DisplayMatState (E, "\n\tE");
891 cS=CopyMatState ( S,cS);
892 cE=CopyMatState ( E,cE);
894 viterbiD_hmm (A,ns,ls, H,CL,cS, cE, NULL);
902 t=H->T[S->st][E->st];
903 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;
912 MatState* RviterbiD_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL,MatState *S, MatState *E, int **seg_list)
918 viterbiD_hmm (A,ns, ls,H, CL, S, E, seg_list);
921 if ( S->n==E)return S;
922 if ( E->sc==H->forbiden)
925 fprintf ( stderr, "\nERROR: The Requested Model (Cf Above) Cannot Produce the Pair-Alignment\nYou must allow extra possible transitions\n[FATAL:%s]", PROGRAM);
926 myexit ( EXIT_FAILURE);
933 d1=MinDeltaMatState(S,E);
934 d2=MaxDeltaMatState(S,E);
935 align=((d1==1 && d2==1) || ( d1==0))?0:1;
936 if (align)RviterbiD_hmm (A,ns, ls,H, CL,S,E, seg_list);
943 #define Dim(i,j) (i*LenJ+j)
944 MatState* viterbiD_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL, MatState *S,MatState *E, int **seg_list)
946 int current, previous,row, prow;
948 int a,b,i,j,pi,pj, s, k;
949 int start_i, start_j, end_i, end_j, l1, l2;
951 MatState *CC, *PCC,*tS, *tE, *mark=NULL;
957 static int **main_jlist;
958 static MatState ***M;
961 static int LenJ, LenI;
962 int MaxDelta=50, DeltaI, DeltaJ;
965 DisplayMatState (S, "\n\tS");
966 DisplayMatState (E, "\n\tE");
971 free_arrayN((void **)main_jlist, 2);main_jlist=NULL;
973 for ( a=0; a<2; a++)memset(M[a],0,LenJ*LenI*sizeof (MatState*));
974 free_arrayN((void **)M, 3);M=NULL;
976 ManageMatState ( FREE_ALL, NULL);
981 if ( MatStateAreIdentical ( S, E))return NULL;
982 l1=strlen (A->seq_al[ls[0][0]]);l2=strlen (A->seq_al[ls[1][0]]);
984 midpoint=S->i+((E->i-S->i)/2);
988 start_i=S->i;end_i=E->i;start_j=S->j;end_j=E->j;
989 current=0;previous=1;
996 M=declare_arrayN(3, sizeof ( MatState),2,LenI*LenJ,0);
997 toclean=vcalloc ( LenI*LenJ, sizeof (int));
1000 if ( !main_jlist)main_jlist= seglist2table(seg_list, l1, l2);
1003 /*MAKE THE VITERBI FROM S(tart) to E(nd)*/
1004 mark=ManageMatState ( MARK, mark);
1005 mode=(!seg_list)?ALL:seg_list[0][0];
1007 for (ntoclean=0,i=start_i; i<=end_i; i++)
1011 if ( mode==ALL)jlist=main_jlist[0];
1012 else if ( mode==DIAGONALS)jlist=(i==0)?main_jlist[0]:main_jlist[1];
1013 else if ( mode==SEGMENTS) jlist=main_jlist[i+2];
1016 for ( dj=1; dj<=jlist[0]; dj++)
1018 DeltaJ=((FABS(dj-i))<MaxDelta)?1:MaxDelta+1; /*Keep Pointers on a 2*DeltaMax Band around the main diagonal*/
1020 dc=(mode==DIAGONALS && dj!=1)?i:0;/*Make sure the diagonal Mode uses the sides of the array*/
1024 if ( j<start_j || j>end_j)continue;
1025 for ( s=H->nS-1;s>=0; s--)
1027 S1=H->S[s];pi=i-S1->DI;prow=S1->DI;
1029 if ( S1->DI && S1->DJ){pj=j-S1->DJ;}
1030 else if ( !S1->DJ)pj=j;
1031 else if ( dj>1)pj=jlist[dj-S1->DJ]+dc;
1034 if (!M[row][Dim(j,s)])toclean[ntoclean]=Dim(j,s);
1036 CC=M[row][Dim(j,s)]=CopyMatState(NULL, M[row][Dim(j,s)]);
1037 CC->i=i; CC->j=j; CC->st=s;PCC=NULL;
1039 if (i==start_i && j==start_j && s==S->st){CC=CopyMatState(S,CC);}
1040 else if ( i==end_i && j==end_j && E->st!=H->end && s!=E->st)CC->sc=H->forbiden;
1041 else if ( pi<start_i || pj<start_j) {CC->sc=H->forbiden;}
1044 for (k=1; k<=H->fromM[S1->state][0]; k++)
1046 S2=H->S[H->fromM[s][k]];
1048 if ( S1->DI && S1->DJ)PCC=M[prow][Dim((j-S1->DJ),(S2->state))];
1049 else PCC=M[prow][Dim((jlist[dj-S1->DJ]+dc),(S2->state))];
1052 else if ( pi+pj!=0 && S2->state==H->start);
1053 else if ( !(pi==l1 && pj==l2) && s==H->end);
1056 v=lu_RviterbiD_hmm(A,ns, ls, H, CL,PCC,CC, NULL);
1057 if ( v!=H->forbiden && (CC->sc==H->forbiden || v> CC->sc)){CC->sc=v; CC->pst=S2->state;CC->p=PCC;}
1061 if (CC->sc==H->forbiden);
1062 else if (i==midpoint || DeltaI<=MaxDelta||DeltaJ<=MaxDelta ||(i==start_i && j==start_j && s==S->st) )
1064 CC->m=(CC->p)?(CC->p)->m:NULL;
1065 PCC=CopyMatState(CC,NULL);
1066 PCC->m=CC->m;CC->m=PCC;
1068 else CC->m=(CC->p)?(CC->p)->m:NULL;
1072 for ( dj=1; dj<=jlist[0] && i!=end_i; dj++)
1074 dc=(mode==DIAGONALS && dj!=1)?i:0;
1076 if ( j<start_j || j>end_j)continue;
1078 for ( s=H->nS-1;s>=0; s--)
1081 CC=(M[prow][Dim(j,s)]);M[prow][Dim(j,s)]=M[row][Dim(j,s)];M[row][Dim(j,s)]=CC;
1082 if (!M[row][Dim(j,s)])toclean[ntoclean++]=Dim(j,s);
1083 if (M[prow][Dim(j,s)]) M[row ][Dim(j,s)]=CopyMatState ( M[prow][Dim(j,s)], M[row][Dim(j,s)]);
1090 mark=ManageMatState ( MARK,mark);
1094 if ( E->st==H->end || E->st==H->forbiden){E=CopyMatState ((M[row][Dim(end_j,E->st)]),E);}
1099 PCC=CopyMatState (M[row][Dim(end_j,E->st)], NULL);
1101 if ( MatStateAreIdentical(PCC,PCC->m))PCC=PCC->m;
1105 tS=CopyMatState (PCC->m,NULL); tS->n=PCC; PCC->p=tS;PCC=tS;
1111 S->n=tS->n; (S->n)->p=S;
1112 E->p=tE->p; (E->p)->n=E;
1115 ManageMatState ( FREE_MARK,mark);
1118 for ( a=0; a<ntoclean; a++)
1120 for ( b=0; b< 2; b++){M[b][toclean[a]]=NULL;}
1124 while (S && S->p!=E){S->m=NULL;S=S->n;}/*Clean the memory of the rturned Cells*/
1128 /*********************************************************************************/
1131 /* HMM Viterbi Diagonals GLOBAL/LOCAL */
1134 /*********************************************************************************/
1136 int viterbiDGL_pair_wise (Alignment *A,int*ns, int **ls,Constraint_list *CL)
1149 A->pos=aln2pos_simple( A, -1, ns, ls);
1151 Start=ManageMatState ( DECLARE, NULL);
1152 End=ManageMatState ( DECLARE, NULL);
1153 H=define_simple_model (CL);
1154 l1=strlen (A->seq_al[ls[0][0]]);
1155 l2=strlen (A->seq_al[ls[1][0]]);
1157 main_i=MAX(1,(l2-l1)+1);
1158 main_j=MAX(1,(l1-l2)+1);
1160 seg_list=declare_arrayN(2, sizeof (int), l1+l2+3, 3);
1161 seg_list[0][0]=DIAGONALS;
1165 for ( b=1,a=l1; a>= 1; a--)
1167 if (a<50 || (b==main_i && a==main_j))
1171 seg_list[c][2]=MIN((l1-a), (l2-b));
1177 for ( b=2,a=1; b<= l2; b++, c++)
1179 if (b<50 || (b==main_i && a==main_j))
1183 seg_list[c][2]=MIN((l1-a), (l2-b));
1188 seg_list[c][0]=FORBIDEN;
1190 Start->i=0 ;Start->j=0 ; Start->st=H->start;Start->sc=0;
1191 End->i =l1; End->j=l2; End ->st=H->end;
1192 Start=RviterbiDGL_hmm (A, ns, ls, H, CL, Start,End,seg_list);
1195 path=traceback (A, ns, ls, H, CL, Start,NULL, NULL);
1199 A=viterbipath2aln (A,ns,ls,path, H);
1201 viterbiD_hmm (NULL, ns, ls, H, CL, Start,End, seg_list);
1203 free_int (A->pos, -1);
1204 free_arrayN((void **)seg_list, 2);
1207 return A->score_aln;
1211 double lu_RviterbiDGL_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL,MatState *S, MatState *E, int **seg_list)
1214 double sc, sc2,e, t;
1215 static MatState *cS=NULL, *cE=NULL;
1217 max=MAX((E->i-S->i), (E->j-S->j));
1218 min=MIN((E->i-S->i), (E->j-S->j));
1223 if ( S==NULL || E==NULL || S->sc==H->forbiden) return H->forbiden;
1224 else if ( S->st==H->start) return 0;
1225 else if ( E->st==H->end) return S->sc;
1228 e=hmm_add(S->sc,H->T[S->st][E->st]);
1229 if ( H->T[E->st][E->st]!=H->forbiden)e=hmm_add(e, (max-1)*H->T[E->st][E->st]);
1230 if ( (H->S[E->st])->em!=H->forbiden) e=hmm_add(e, max *(H->S[E->st])->em );
1233 else if ( min>0 && max>1)
1236 fprintf ( stderr, "\nWarning: Disjoined Diagonals");
1237 DisplayMatState (S, "\n\tS");
1238 DisplayMatState (E, "\n\tE");
1241 cS=CopyMatState ( S,cS);
1242 cE=CopyMatState ( E,cE);
1244 viterbiD_hmm (A,ns,ls, H,CL,cS, cE, NULL);
1252 t=H->T[S->st][E->st];
1253 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;
1254 sc=hmm_add(S->sc,t);
1262 MatState* RviterbiDGL_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL,MatState *S, MatState *E, int **seg_list)
1268 viterbiDGL_hmm (A,ns, ls,H, CL, S, E, seg_list);
1271 if ( S->n==E)return S;
1272 if ( E->sc==H->forbiden)
1275 fprintf ( stderr, "\nERROR: The Requested Model (Cf Above) Cannot Produce the Pair-Alignment\nYou must allow extra possible transitions\n[FATAL:%s]", PROGRAM);
1276 myexit ( EXIT_FAILURE);
1283 d1=MinDeltaMatState(S,E);
1284 d2=MaxDeltaMatState(S,E);
1285 align=((d1==1 && d2==1) || ( d1==0))?0:1;
1286 if (align)RviterbiDGL_hmm (A,ns, ls,H, CL,S,E, seg_list);
1294 #define Dim(i,j) (i*LenJ+j)
1295 MatState* viterbiDGL_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL, MatState *S,MatState *E, int **seg_list)
1297 int current, previous,row, prow;
1299 int a,i,j,pi,pj, s, k;
1300 int start_i, start_j, end_i, end_j, l1, l2;
1302 MatState *CC, *PCC,*tS, *tE,*bestE,*bestS, *mark=NULL;
1308 static int **main_jlist;
1309 static MatState ***M;
1310 static int *toclean;
1312 static int LenJ, LenI;
1313 int MaxDelta=50, DeltaI, DeltaJ;
1318 DisplayMatState (S, "\n\tS");
1319 DisplayMatState (E, "\n\tE");
1324 free_arrayN((void **)main_jlist, 2);main_jlist=NULL;
1326 for ( a=0; a<2; a++)memset(M[a],0,LenJ*LenI*sizeof (MatState*));
1327 free_arrayN((void **)M, 3);M=NULL;
1329 ManageMatState ( FREE_ALL, NULL);
1334 if ( MatStateAreIdentical ( S, E))return NULL;
1335 l1=strlen (A->seq_al[ls[0][0]]);l2=strlen (A->seq_al[ls[1][0]]);
1337 midpoint=S->i+((E->i-S->i)/2);
1341 start_i=S->i;end_i=E->i;start_j=S->j;end_j=E->j;
1342 current=0;previous=1;
1349 M=declare_arrayN(3, sizeof ( MatState),2,LenI*LenJ,0);
1350 toclean=vcalloc ( LenI*LenJ, sizeof (int));
1353 if ( !main_jlist)main_jlist= seglist2table(seg_list, l1, l2);
1356 /*MAKE THE VITERBI FROM S(tart) to E(nd)*/
1357 mark=ManageMatState ( MARK, mark);
1358 mode=(!seg_list)?ALL:seg_list[0][0];
1359 bestE=CopyMatState (E, NULL);
1360 bestS=CopyMatState (NULL, NULL);
1361 for (ntoclean=0,i=start_i; i<=end_i; i++)
1365 if ( mode==ALL)jlist=main_jlist[0];
1366 else if ( mode==DIAGONALS)jlist=(i==0)?main_jlist[0]:main_jlist[1];
1367 else if ( mode==SEGMENTS) jlist=main_jlist[i+2];
1370 for ( dj=1; dj<=jlist[0]; dj++)
1372 DeltaJ=(FABS(dj-i)<MaxDelta)?1:MaxDelta+1; /*Keep Pointers on a 2*DeltaMax Band around the main diagonal*/
1374 dc=(mode==DIAGONALS && dj!=1)?i:0;/*Make sure the diagonal Mode uses the sides of the array*/
1378 if ( j<start_j || j>end_j)continue;
1379 for ( s=H->nS-1;s>=0; s--)
1381 if ( s==S->st)continue;
1382 S1=H->S[s];pi=i-S1->DI;prow=S1->DI;
1384 if ( S1->DI && S1->DJ){pj=j-S1->DJ;}
1385 else if ( !S1->DJ)pj=j;
1386 else if ( dj>1)pj=jlist[dj-S1->DJ]+dc;
1389 if (!M[row][Dim(j,s)])toclean[ntoclean]=Dim(j,s);
1391 CC=M[row][Dim(j,s)]=CopyMatState(NULL, M[row][Dim(j,s)]);
1392 CC->i=i; CC->j=j; CC->st=s;PCC=NULL;
1394 if (i==start_i && j==start_j && s==S->st){CC=CopyMatState(S,CC);}
1395 else if ( s==S->st);
1396 else if ( i==end_i && j==end_j && E->st!=H->end && s!=E->st)CC->sc=H->forbiden;
1398 else if ( pi<start_i || pj<start_j) {CC->sc=H->forbiden;}
1401 for (k=1; k<=H->fromM[S1->state][0]; k++)
1403 S2=H->S[H->fromM[s][k]];
1405 if ( S1->DI && S1->DJ)PCC=M[prow][Dim((j-S1->DJ),(S2->state))];
1406 else PCC=M[prow][Dim((jlist[dj-S1->DJ]+dc),(S2->state))];
1408 if ( S2->state==H->start){PCC=bestS;PCC->st=0;PCC->sc=0;PCC->m=PCC->n=PCC->p=NULL;}
1410 v=lu_RviterbiDGL_hmm(A,ns, ls, H, CL,PCC,CC, NULL);
1411 if ( v!=H->forbiden && (CC->sc==H->forbiden || v> CC->sc)){CC->sc=v; CC->pst=S2->state;CC->p=PCC;}
1414 if ( CC->sc==H->forbiden);
1415 else if ( bestE->sc==H->forbiden || bestE->sc>CC->sc)
1417 bestE=CopyMatState(CC, bestE);
1418 bestE->m=(CC->p)->m;
1420 else if (CC->p && (CC->p)->st==H->start)
1422 CC->m=CopyMatState (CC->p, NULL);
1424 else if (i==midpoint || DeltaI<=MaxDelta||DeltaJ<=MaxDelta ||(i==start_i && j==start_j && s==S->st) )
1426 CC->m=(CC->p)?(CC->p)->m:NULL;
1427 PCC=CopyMatState(CC,NULL);
1428 PCC->m=CC->m;CC->m=PCC;
1430 else CC->m=(CC->p)?(CC->p)->m:NULL;
1434 for ( dj=1; dj<=jlist[0] && i!=end_i; dj++)
1436 dc=(mode==DIAGONALS && dj!=1)?i:0;
1438 if ( j<start_j || j>end_j)continue;
1440 for ( s=H->nS-1;s>=0; s--)
1443 CC=(M[prow][Dim(j,s)]);M[prow][Dim(j,s)]=M[row][Dim(j,s)];M[row][Dim(j,s)]=CC;
1444 /*if (!M[row][Dim(j,s)])toclean[ntoclean++]=Dim(j,s);*/
1445 if (M[prow][Dim(j,s)]) M[row ][Dim(j,s)]=CopyMatState ( M[prow][Dim(j,s)], M[row][Dim(j,s)]);
1452 mark=ManageMatState ( MARK,mark);
1456 if ( E->st==H->end || E->st==H->forbiden){E=CopyMatState ((M[row][Dim(end_j,E->st)]),E);}
1457 PCC=CopyMatState (bestE, NULL);
1459 if ( MatStateAreIdentical(PCC,PCC->m))PCC=PCC->m;
1463 tS=CopyMatState (PCC->m,NULL); tS->n=PCC; PCC->p=tS;PCC=tS;
1469 CopyMatState ( tS, S);
1470 CopyMatState ( tE, E);
1472 ManageMatState ( FREE_MARK,mark);
1474 for ( a=0; a<2; a++)memset(M[a],0,LenJ*LenI*sizeof (MatState*));
1476 while (S && S->p!=E){S->m=NULL;S=S->n;}/*Clean the memory of the rturned Cells*/
1481 /*********************************************************************************/
1484 /* HMM Viterbi Diagonals PROCESSING */
1487 /*********************************************************************************/
1488 int **seglist2table ( int **seglist,int l1, int l2)
1492 int line, a,si, sj, ei, j, c;
1495 valuesT=vcalloc ((l1+2)+3, sizeof (int*));
1496 valuesT[0]=vcalloc (l2+2, sizeof (int));
1497 for (a=0; a<=l2; a++)valuesT[0][++valuesT[0][0]]=a;
1499 if ( !seglist) return valuesT;
1501 valuesT[1]=vcalloc (l1+l2+2, sizeof (int));
1502 bvalues=vcalloc (l1+l2+2, sizeof (int));
1504 while (seglist[c][0]!=FORBIDEN)
1510 bvalues[(sj-si)+l1]=1;
1513 valuesT[1][++valuesT[1][0]]=0;
1514 for (a=0; a<=(l1+l2); a++)
1518 valuesT[1][++valuesT[1][0]]=a-l1;
1525 valuesT[2]=vcalloc (l2+2, sizeof (int));
1526 for (a=0; a<=l2; a++)valuesT[2][++valuesT[2][0]]=a;
1528 bvalues=vcalloc (l2+2, sizeof (int));
1529 for ( line=1; line<=l1; line++)
1532 bvalues[++bvalues[0]]=0;
1533 while (seglist[c][0]!=FORBIDEN)
1536 ei=si+seglist[c][2];
1539 if ( line<si || line>ei);
1540 else if (j>=0 && j<=l2 && seglist[c][2])
1542 bvalues[++bvalues[0]]=j;
1546 valuesT[line+2]=vcalloc (bvalues[0]+1, sizeof (int));
1547 for ( a=0; a<=bvalues[0]; a++) valuesT[line+2][a]=bvalues[a];
1555 /*********************************************************************************/
1561 /*********************************************************************************/
1563 Hmm* declare_hmm(int n)
1568 H=vcalloc (1, sizeof (Hmm));
1570 H->S=vcalloc (H->nS, sizeof (HmmState*));
1571 for (a=0; a<H->nS; a++)
1573 H->S[a]=vcalloc (1, sizeof (HmmState));
1574 (H->S[a])->em2=vcalloc (MAX_EMISSION, sizeof (double));
1576 (H->S[a])->T=vcalloc ( H->nS, sizeof (StateTrans*));
1577 for ( b=0; b< H->nS; b++)
1578 (H->S[a])->T[b]=vcalloc (1, sizeof (StateTrans));
1583 Hmm* free_Hmm(Hmm*H)
1587 H=vcalloc (1, sizeof (Hmm));
1588 free_double (H->T, -1);
1589 free_int ( H->fromM, -1);
1590 free_int ( H->toM, -1);
1592 for (a=0; a< H->nS; a++)
1595 for ( b=0; b< H->nS; b++)
1597 vfree ((H->S[a])->em2);
1598 vfree((H->S[a])->T[b]);
1600 vfree((H->S[a])->T);
1608 void DisplayHmm ( Hmm *H)
1613 for ( a=0; a< H->nS; a++)
1616 fprintf ( stderr, "\nSTATE %d: %s\n",S1->state,S1->name);
1617 fprintf ( stderr, "\n\tDI %d", S1->DI);
1618 fprintf ( stderr, "\n\tDJ %d", S1->DJ);
1619 fprintf ( stderr, "\n\tE %f", (float)S1->em);
1621 fprintf ( stderr, "\nReached FROM: ");
1622 for ( b=1; b<=H->fromM[a][0]; b++)
1624 S2=H->S[H->fromM[a][b]];
1625 fprintf ( stderr, "[ %s %f ] ", S2->name, H->T[S2->state][S1->state]);
1627 fprintf ( stderr, "\nGoes TO: ");
1628 for ( b=1; b<=H->toM[a][0]; b++)
1630 S2=H->S[H->toM[a][b]];
1631 fprintf ( stderr, "[ %s %f ] ", S2->name, H->T[S1->state][S2->state]);
1636 Hmm * bound_hmm ( Hmm *H)
1642 name=declare_char(H->nS, 100);
1643 H->T=declare_double ( H->nS, H->nS);
1645 for ( a=0; a< H->nS; a++)
1647 sprintf ( name[a], "%s", (H->S[a])->name);
1648 H->order=MAX(H->order, (H->S[a])->DI);
1649 H->order=MAX(H->order, (H->S[a])->DJ);
1652 for ( a=0; a< H->nS; a++)for (b=0; b<H->nS; b++)H->T[a][b]=H->forbiden;
1653 for (a=0; a< H->nS; a++)
1656 for ( b=0; b< S->nT; b++)
1658 c=name_is_in_list ((S->T[b])->name, name, H->nS, 100);
1659 if ( c!=-1)H->T[a][c]=(S->T[b])->tr;
1664 bM[state][0]=n_allowed transitions
1665 bM[state][1]=first allowed transition
1668 H->toM=declare_int ( H->nS, H->nS);
1669 H->fromM=declare_int ( H->nS, H->nS);
1671 for ( a=0; a< H->nS; a++)
1672 for ( b=0; b< H->nS; b++)
1674 if ( H->T[a][b]!=H->forbiden )
1676 {H->fromM[b][0]++; H->fromM[b][H->fromM[b][0]]=a;}
1677 {H->toM[a][0]++; H->toM[a][H->toM[a][0]]=b;}
1680 for ( a=0; a< H->nS; a++)
1682 if (( H->S[a])->em!=H->forbiden)( H->S[a])->em*=SCORE_K;
1683 for ( b=0; b< H->nS; b++)
1684 if ( H->T[a][b]!=H->forbiden)H->T[a][b]*=SCORE_K;
1686 free_arrayN((void**)name, 2);
1690 /*********************************************************************************/
1693 /* Memory Management */
1696 /*********************************************************************************/
1698 MatState * ManageMatState(int Mode, MatState *C)
1700 static MatState *Fheap;
1701 static MatState *Aheap;
1702 MatState *Cmark, *Pmark;
1703 static int alloc, free;
1704 if (!Fheap || Fheap->Hp==NULL)
1708 if (!Fheap){Fheap=vcalloc (1, sizeof (MatState));Fheap->free=1;free++;}
1709 if (!Aheap)Aheap=vcalloc (1, sizeof (MatState));
1710 while ( c!=extension)
1712 C=vcalloc ( 1, sizeof (MatState));
1713 C->free=1;Fheap->Hn=C;C->Hp=Fheap;
1726 if ( Aheap){Aheap->Hn=C;C->Hp=Aheap;Aheap=C;}
1731 C=CopyMatState(NULL, C);
1734 else if ( Mode==FREE)
1736 if ( !C || C->free==1);
1739 C=CopyMatState(NULL, C);
1741 if (C->Hp==NULL && C==Aheap)crash ("");
1742 if (C==Aheap)Aheap=C->Hp;
1743 if (C->Hn){(C->Hn)->Hp=C->Hp;}
1744 if (C->Hp){(C->Hp)->Hn=C->Hn;}
1746 Fheap->Hn=C;C->Hp=Fheap;
1753 else if ( Mode==FREE_ALL)
1768 else if ( Mode==INFO)
1770 fprintf ( stderr, "\nAllocated: %d Free %d", alloc, free);
1772 else if ( Mode==MARK)
1776 else {C->Mn=Aheap;Aheap->Mp=C;}
1780 else if ( Mode==UNMARK)
1784 else if ( Mode == FREE_MARK)
1789 if ( Cmark==Pmark)return NULL;
1790 else if ( Cmark==Aheap)
1791 {Aheap=Pmark;C=Pmark->Hn;Pmark->Hn=NULL;}
1794 (Cmark->Hn)->Hp=Pmark;
1796 Pmark->Hn=Cmark->Hn;
1805 while (C && !C->free)
1817 MatState* CopyMatState ( MatState*I, MatState*O)
1819 if (O==NULL || O->free==1) O=ManageMatState(DECLARE, NULL);
1820 if (I==NULL || I->free==1)I=NULL;
1823 O->st =(I)?I->st:FORBIDEN;
1824 O->pst=(I)?I->pst:FORBIDEN;
1825 O->sc =(I)?I->sc:FORBIDEN;
1826 O->n =(I)?I->n:NULL;
1827 O->p =(I)?I->p:NULL;
1828 O->m =(I)?I->m:NULL;
1829 O->s =(I)?I->m:NULL;
1834 /*********************************************************************************/
1840 /*********************************************************************************/
1841 int MaxDeltaMatState (MatState*S, MatState*E)
1843 if ( !S || !E) return -1;
1844 else return MAX((E->i-S->i),(E->j-S->j));
1846 int MinDeltaMatState (MatState*S, MatState*E)
1848 if ( !S || !E) return -1;
1849 else return MIN((E->i-S->i),(E->j-S->j));
1851 int MatStateAreIdentical (MatState*I, MatState*O)
1853 if ( !I || !O)return 0;
1855 if ( I->i!=O->i)return 0;
1856 if ( I->j!=O->j)return 0;
1857 if ( I->st!=O->st)return 0;
1863 Hmm* define_probcons_model(Constraint_list *CL)
1878 H->forbiden=FORBIDEN;
1879 H->start=START_STATE;
1884 sprintf (S->name, "START"); S->state=n;
1890 sprintf ( (S->T[S->nT])->name, "MATCH") ;(S->T[S->nT])->tr=freeT ;S->nT++;
1891 sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=freeT ;S->nT++;
1892 sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=freeT ;S->nT++;
1893 sprintf ( (S->T[S->nT])->name, "END") ;(S->T[S->nT])->tr=freeT ;S->nT++;
1898 sprintf (S->name, "END"); S->state=n;
1906 sprintf (S->name, "MATCH"); S->state=n;
1910 S->em_func=CL->get_dp_cost;
1912 sprintf ( (S->T[S->nT])->name, "MATCH") ;(S->T[S->nT])->tr=freeT;S->nT++;
1913 sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=gop ;S->nT++;
1914 sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=gop ;S->nT++;
1915 sprintf ( (S->T[S->nT])->name, "END"); (S->T[S->nT])->tr=freeT;S->nT++;
1921 sprintf (S->name, "INSERT"); S->state=n;
1925 sprintf ( (S->T[S->nT])->name, "MATCH") ; (S->T[S->nT])->tr=freeT;S->nT++;
1926 sprintf ( (S->T[S->nT])->name, "INSERT"); (S->T[S->nT])->tr=freeT;S->nT++;
1927 sprintf ( (S->T[S->nT])->name, "LINSERT");(S->T[S->nT])->tr=lgop ;S->nT++;
1928 sprintf ( (S->T[S->nT])->name, "END"); (S->T[S->nT])->tr=-gop;S->nT++;
1934 sprintf (S->name, "LINSERT"); S->state=n;
1939 sprintf ( (S->T[S->nT])->name, "INSERT") ;(S->T[S->nT])->tr=freeT;S->nT++;
1940 sprintf ( (S->T[S->nT])->name, "LINSERT");(S->T[S->nT])->tr=freeT;S->nT++;
1947 sprintf (S->name, "DELETE"); S->state=n;
1952 sprintf ( (S->T[S->nT])->name, "MATCH") ;(S->T[S->nT])->tr=freeT;S->nT++;
1953 sprintf ( (S->T[S->nT])->name, "DELETE") ;(S->T[S->nT])->tr=freeT;S->nT++;
1954 sprintf ( (S->T[S->nT])->name, "LDELETE") ;(S->T[S->nT])->tr=lgop ;S->nT++;
1955 sprintf ( (S->T[S->nT])->name, "END") ;(S->T[S->nT])->tr=-gop;S->nT++;
1961 sprintf (S->name, "LDELETE"); S->state=n;
1965 sprintf ( (S->T[S->nT])->name, "DELETE") ;(S->T[S->nT])->tr=freeT;S->nT++;
1966 sprintf ( (S->T[S->nT])->name, "LDELETE");(S->T[S->nT])->tr=freeT;S->nT++;
1975 Hmm* define_mnm_model(Constraint_list *CL)
1990 H->forbiden=FORBIDEN;
1991 H->start=START_STATE;
1996 sprintf (S->name, "START"); S->state=n;
2002 sprintf ( (S->T[S->nT])->name, "MATCH") ;(S->T[S->nT])->tr=freeT ;S->nT++;
2003 sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=freeT ;S->nT++;
2004 sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=freeT ;S->nT++;
2005 sprintf ( (S->T[S->nT])->name, "NOMATCH");(S->T[S->nT])->tr=freeT ;S->nT++;
2006 sprintf ( (S->T[S->nT])->name, "END") ;(S->T[S->nT])->tr=freeT ;S->nT++;
2011 sprintf (S->name, "END"); S->state=n;
2019 sprintf (S->name, "MATCH"); S->state=n;
2023 S->em_func=CL->get_dp_cost;
2025 sprintf ( (S->T[S->nT])->name, "MATCH") ;(S->T[S->nT])->tr=freeT;S->nT++;
2026 sprintf ( (S->T[S->nT])->name, "NOMATCH");(S->T[S->nT])->tr=gop ;S->nT++;
2027 sprintf ( (S->T[S->nT])->name, "END"); (S->T[S->nT])->tr=freeT;S->nT++;
2033 sprintf (S->name, "NOMATCH"); S->state=n;
2039 sprintf ( (S->T[S->nT])->name, "NOMATCH") ;(S->T[S->nT])->tr=freeT;S->nT++;
2040 sprintf ( (S->T[S->nT])->name, "MATCH") ;(S->T[S->nT])->tr=freeT;S->nT++;
2041 sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=freeT ;S->nT++;
2042 sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=freeT ;S->nT++;
2043 sprintf ( (S->T[S->nT])->name, "END"); (S->T[S->nT])->tr=freeT;S->nT++;
2048 sprintf (S->name, "INSERT"); S->state=n;
2052 sprintf ( (S->T[S->nT])->name, "NOMATCH") ; (S->T[S->nT])->tr=freeT;S->nT++;
2053 sprintf ( (S->T[S->nT])->name, "INSERT") ; (S->T[S->nT])->tr=freeT;S->nT++;
2054 sprintf ( (S->T[S->nT])->name, "END"); (S->T[S->nT])->tr=freeT;S->nT++;
2060 sprintf (S->name, "DELETE"); S->state=n;
2065 sprintf ( (S->T[S->nT])->name, "NOMATCH") ;(S->T[S->nT])->tr=freeT;S->nT++;
2066 sprintf ( (S->T[S->nT])->name, "DELETE") ;(S->T[S->nT])->tr=freeT;S->nT++;
2067 sprintf ( (S->T[S->nT])->name, "END") ;(S->T[S->nT])->tr=freeT;S->nT++;
2076 Hmm* define_simple_model(Constraint_list *CL)
2089 H->forbiden=FORBIDEN;
2090 H->start=START_STATE;
2095 sprintf (S->name, "START"); S->state=n;
2101 sprintf ( (S->T[S->nT])->name, "MATCH") ;(S->T[S->nT])->tr=freeT ;S->nT++;
2102 sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=freeT ;S->nT++;
2103 sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=freeT ;S->nT++;
2104 sprintf ( (S->T[S->nT])->name, "END") ;(S->T[S->nT])->tr=freeT ;S->nT++;
2109 sprintf (S->name, "END"); S->state=n;
2117 sprintf (S->name, "MATCH"); S->state=n;
2121 S->em_func=CL->get_dp_cost;
2123 sprintf ( (S->T[S->nT])->name, "MATCH") ;(S->T[S->nT])->tr=freeT;S->nT++;
2124 sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=gop ;S->nT++;
2125 sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=gop ;S->nT++;
2126 sprintf ( (S->T[S->nT])->name, "END"); (S->T[S->nT])->tr=freeT;S->nT++;
2132 sprintf (S->name, "INSERT"); S->state=n;
2136 sprintf ( (S->T[S->nT])->name, "MATCH") ; (S->T[S->nT])->tr=freeT;S->nT++;
2137 sprintf ( (S->T[S->nT])->name, "INSERT"); (S->T[S->nT])->tr=freeT;S->nT++;
2138 sprintf ( (S->T[S->nT])->name, "DELETE"); (S->T[S->nT])->tr=freeT;S->nT++;
2140 sprintf ( (S->T[S->nT])->name, "END"); (S->T[S->nT])->tr=-gop;S->nT++;
2147 sprintf (S->name, "DELETE"); S->state=n;
2152 sprintf ( (S->T[S->nT])->name, "MATCH") ;(S->T[S->nT])->tr=freeT;S->nT++;
2153 sprintf ( (S->T[S->nT])->name, "DELETE") ;(S->T[S->nT])->tr=freeT;S->nT++;
2154 sprintf ( (S->T[S->nT])->name, "INSERT"); (S->T[S->nT])->tr=freeT;S->nT++;
2155 sprintf ( (S->T[S->nT])->name, "END") ;(S->T[S->nT])->tr=-gop;S->nT++;
2163 Hmm* define_two_mat_model(Constraint_list *CL)
2178 H->forbiden=FORBIDEN;
2179 H->start=START_STATE;
2184 sprintf (S->name, "START"); S->state=n;
2190 sprintf ( (S->T[S->nT])->name, "MATCH1") ;(S->T[S->nT])->tr=freeT ;S->nT++;
2191 sprintf ( (S->T[S->nT])->name, "MATCH2") ;(S->T[S->nT])->tr=freeT ;S->nT++;
2193 sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=freeT ;S->nT++;
2194 sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=freeT ;S->nT++;
2195 sprintf ( (S->T[S->nT])->name, "END") ;(S->T[S->nT])->tr=freeT ;S->nT++;
2200 sprintf (S->name, "END"); S->state=n;
2208 sprintf (S->name, "MATCH1"); S->state=n;
2212 S->em_func=get_dp_cost_pam_matrix;
2214 sprintf ( (S->T[S->nT])->name, "MATCH1") ;(S->T[S->nT])->tr=freeT;S->nT++;
2215 sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=gop ;S->nT++;
2216 sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=gop ;S->nT++;
2217 sprintf ( (S->T[S->nT])->name, "END"); (S->T[S->nT])->tr=freeT;S->nT++;
2223 sprintf (S->name, "MATCH2"); S->state=n;
2227 S->em_func=get_dp_cost_blosum_matrix;
2229 sprintf ( (S->T[S->nT])->name, "MATCH2") ;(S->T[S->nT])->tr=freeT;S->nT++;
2230 sprintf ( (S->T[S->nT])->name, "INSERT");(S->T[S->nT])->tr=gop ;S->nT++;
2231 sprintf ( (S->T[S->nT])->name, "DELETE");(S->T[S->nT])->tr=gop ;S->nT++;
2232 sprintf ( (S->T[S->nT])->name, "END"); (S->T[S->nT])->tr=freeT;S->nT++;
2238 sprintf (S->name, "INSERT"); S->state=n;
2242 sprintf ( (S->T[S->nT])->name, "MATCH2") ; (S->T[S->nT])->tr=freeT;S->nT++;
2243 sprintf ( (S->T[S->nT])->name, "MATCH1") ; (S->T[S->nT])->tr=freeT;S->nT++;
2245 sprintf ( (S->T[S->nT])->name, "INSERT"); (S->T[S->nT])->tr=freeT;S->nT++;
2246 sprintf ( (S->T[S->nT])->name, "LINSERT");(S->T[S->nT])->tr=lgop ;S->nT++;
2247 sprintf ( (S->T[S->nT])->name, "END"); (S->T[S->nT])->tr=-gop;S->nT++;
2253 sprintf (S->name, "LINSERT"); S->state=n;
2258 sprintf ( (S->T[S->nT])->name, "INSERT") ;(S->T[S->nT])->tr=freeT;S->nT++;
2259 sprintf ( (S->T[S->nT])->name, "LINSERT");(S->T[S->nT])->tr=freeT;S->nT++;
2266 sprintf (S->name, "DELETE"); S->state=n;
2271 sprintf ( (S->T[S->nT])->name, "MATCH2") ;(S->T[S->nT])->tr=freeT;S->nT++;
2272 sprintf ( (S->T[S->nT])->name, "MATCH1") ;(S->T[S->nT])->tr=freeT;S->nT++;
2274 sprintf ( (S->T[S->nT])->name, "DELETE") ;(S->T[S->nT])->tr=freeT;S->nT++;
2275 sprintf ( (S->T[S->nT])->name, "LDELETE") ;(S->T[S->nT])->tr=lgop ;S->nT++;
2276 sprintf ( (S->T[S->nT])->name, "END") ;(S->T[S->nT])->tr=-gop;S->nT++;
2281 sprintf (S->name, "LDELETE"); S->state=n;
2285 sprintf ( (S->T[S->nT])->name, "DELETE") ;(S->T[S->nT])->tr=freeT;S->nT++;
2286 sprintf ( (S->T[S->nT])->name, "LDELETE");(S->T[S->nT])->tr=freeT;S->nT++;
2291 fprintf ( stderr, "\nERROR in HMM definition [FATAL:%s]", PROGRAM);
2292 myexit (EXIT_FAILURE);
2298 void DisplayMatState ( MatState *S, char *s)
2300 if ( S==NULL)fprintf ( stderr, "%s: Cell is undefined", s);
2301 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);
2303 void testfunc ( MatState *S, char *s)
2305 if ( S==NULL)return;
2306 fprintf ( stderr, "\n#### %s ", s);
2307 while ( S){DisplayMatState ( S,"\n\t");S=S->n;}
2308 fprintf ( stderr, "\n");
2346 m=emit_pair_default[alphabetDefault[seq1[a]]][alphabetDefault[seq2[a]]];
2348 F[M][i][j]=F[M][i-1][j-1];
2352 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]);
2353 M[D ][i][j]=log_add3(gep,M[Match][i ][j-step_j]+gop,M[D][i ][j-step_j]);
2354 M[I ][i][j]=log_add3(gep,M[Match][i-step_i][j ]+gop,M[I][i-step_i][j ]);
2357 M[Match][i][j]=log_add3(M[Match][i][j], M[LI][i-step_i][j],M[LD][i][j-step_j]);
2358 M[LI ][i][j]=log_add3(lgep, M[I][i-step_i][j ]+lgop,M[LI][i-step_i][j ]);
2359 M[LD ][i][j]=log_add3(lgep, M[D][i ][j-step_j]+lgop,M[LD][i ][j-step_j]);
2365 MatState* RviterbiL_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL,MatState *S, MatState *E)
2368 Mid=viterbiL_hmm (A,ns, ls,H, CL, S, E);
2383 RviterbiL_hmm (A,ns, ls,H, CL,S, Mid);
2384 RviterbiL_hmm (A,ns, ls,H, CL,Mid, E);
2389 MatState* viterbiL_hmm (Alignment *A,int *ns, int **ls, Hmm *H, Constraint_list *CL, MatState *S,MatState *E)
2391 int current,memory, dim;
2393 int i,j,pi,pj, s, k;
2394 int start_i, start_j, end_i, end_j, l1, l2;
2396 static MatState ****M;
2405 if ( MatStateAreIdentical ( S, E))return NULL;
2407 l1=strlen (A->seq_al[ls[0][0]]);
2408 l2=strlen (A->seq_al[ls[1][0]]);
2410 midpoint=S->i+(E->i-S->i)/2;
2413 start_i=S->i;end_i=E->i;
2414 start_j=S->j;end_j=E->j;
2416 dim=H->order+2;current=0;memory=H->order+1;
2417 if (!M || (l2+1)>maxl)
2418 { free_arrayN((void **)M, 4);
2419 M=declare_arrayN(4, sizeof ( MatState), dim, maxl=(l2+1), H->nS,1);
2422 /*MAKE THE VITERBI FROM S(tart) to E(nd)*/
2423 for ( i=start_i; i<=end_i; i++)
2425 M= (MatState****)recycle ( (void **)M,H->order+1,1);
2426 for ( j=start_j; j<=end_j; j++)
2428 for ( s=H->nS-1;s>=0; s--)
2432 pi=i-S1->DI;pj=j-S1->DJ;
2433 CC=M[current][j][s];
2434 CC->i=i; CC->j=j; CC->st=s;CC->sc=H->forbiden;CC->p=CC->n=CC->m=NULL;CC->sc=H->forbiden;
2435 if (i==start_i && j==start_j && s==S->st) {CopyMatState(S,CC);}
2436 else if ( i==end_i && j==end_j && s==E->st && s!=H->end)
2440 CC->p=M[S1->DI][j-S1->DJ][S2->state];
2442 else if ( pi<start_i || pj<start_j) {CC->sc=H->forbiden;}
2445 for (k=1; k<=H->fromM[S1->state][0]; k++)
2447 S2=H->S[H->fromM[s][k]];
2448 PCC=M[S1->DI][j-S1->DJ][S2->state];
2450 if ( pi+pj!=0 && S2->state==H->start) {t=H->forbiden;}
2451 else if ( !(pi==l1 && pj==l2) && s==H->end){t=H->forbiden;}
2452 else t=H->T[S2->state][S1->state];
2454 v=hmm_add(t,PCC->sc);
2455 if ( v!=H->forbiden && (CC->sc==H->forbiden || v> CC->sc)){CC->sc=v; CC->pst=S2->state;CC->p=PCC;}
2458 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;
2459 CC->sc=hmm_add(CC->sc,e);
2462 if (i==midpoint)CC->m=CopyMatState(CC, M[memory][j][s]);
2463 else if (i>midpoint && CC->sc!=H->forbiden) CC->m=(M[S1->DI][j-S1->DJ][CC->pst])->m;
2468 if ( E->st==H->end)CopyMatState ((M[current][end_j][E->st]),E);
2472 Mid=CopyMatState ((M[current][end_j][E->st])->m,NULL);
2476 CC=M[current][E->j][E->st];
2478 while (!MatStateAreIdentical (CC->p, S) )
2480 Mid->p=CopyMatState(CC->p,NULL);
2482 Mid=Mid->p;CC=CC->p;
2492 /*********************************COPYRIGHT NOTICE**********************************/
2493 /*© Centre National de la Recherche Scientifique (CNRS) */
2495 /*Please Cite: Notredame*/
2496 /*Mon May 17 20:15:35 MDT 2004. */
2497 /*All rights reserved.*/
2499 /* This file is an integral part of the */
2500 /* ALIGN_TWO_SEQ Software. */
2501 /* Its content is protected and all */
2502 /* the conditions mentioned in the licensing */
2503 /* agreement of the software apply to this file.*/
2504 /*............................................... |*/
2505 /* If you need some more information, or if you */
2506 /* wish to obtain a full license, please contact: */
2507 /* cedric.notredame@europe.com*/
2508 /*............................................... |*/
2512 /*********************************COPYRIGHT NOTICE**********************************/
2513 /*********************************COPYRIGHT NOTICE**********************************/
2514 /*© Centro de Regulacio Genomica */
2516 /*Cedric Notredame */
2517 /*Tue Oct 27 10:12:26 WEST 2009. */
2518 /*All rights reserved.*/
2519 /*This file is part of T-COFFEE.*/
2521 /* T-COFFEE is free software; you can redistribute it and/or modify*/
2522 /* it under the terms of the GNU General Public License as published by*/
2523 /* the Free Software Foundation; either version 2 of the License, or*/
2524 /* (at your option) any later version.*/
2526 /* T-COFFEE is distributed in the hope that it will be useful,*/
2527 /* but WITHOUT ANY WARRANTY; without even the implied warranty of*/
2528 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the*/
2529 /* GNU General Public License for more details.*/
2531 /* You should have received a copy of the GNU General Public License*/
2532 /* along with Foobar; if not, write to the Free Software*/
2533 /* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA*/
2534 /*............................................... |*/
2535 /* If you need some more information*/
2536 /* cedric.notredame@europe.com*/
2537 /*............................................... |*/
2541 /*********************************COPYRIGHT NOTICE**********************************/