--- /dev/null
+#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 */
+/*Fri Feb 18 08:27:45 CET 2011 - Revision 596. */
+/*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*******************************/
git://source.jalview.org / jabaws.git / blobdiff