7 #include "io_lib_header.h"
8 #include "util_lib_header.h"
9 #include "define_header.h"
10 #include "dp_lib_header.h"
13 int cl2pair_list_ecl ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
17 /*******************************************************************************/
18 /* idscore_pairseq: measure the % id without delivering thze aln*/
20 /* makes DP between the the ns[0] sequences and the ns[1] */
22 /* for MODE, see the function get_dp_cost */
23 /*******************************************************************************/
24 int idscore_pairseq (char *s1, char *s2, int gop, int gep, int **m, char *comp_mode)
26 int **I, **D, **M, *P;
27 int i, j,l1, l2, l,score, id, igop,match;
30 l1=strlen (s1); l2=strlen (s2);
31 lower_string (s1); lower_string (s2);
33 I=declare_int (6,l2+1);D=declare_int (6,l2+1);M=declare_int (6,l2+1);
36 D[0][j]=gep*j;M[0][j]=2*gep*j;D[4][j]=0;
47 score=m[s1[i-1]-'a'][s2[j-1]-'a'];
48 id=(s1[i-1]==s2[j-1])?1:0;
50 igop=(i==l1 || j==l2)?0:gop;
52 if ((D[0][j]+gep)>(M[0][j]+igop+gep)) {D[1][j]=D[0][j]+gep; D[3][j]=D[2][j]; D[5][j]=D[4][j];}
53 else {D[1][j]=M[0][j]+igop+gep; D[3][j]=M[2][j]; D[5][j]=M[4][j];}
55 if ( (I[1][j-1]+gep)>(M[1][j-1]+igop+gep)){I[1][j]=I[1][j-1]+gep; I[3][j]=I[3][j-1]; I[5][j]=I[5][j-1];}
56 else {I[1][j]=M[1][j-1]+igop+gep; I[3][j]=M[3][j-1]; I[5][j]=M[5][j-1];}
58 match=M[0][j-1]+score;
59 if (I[1][j]>match && I[1][j]>D[1][j]) {M[1][j]=I[1][j] ; M[3][j]=I[3][j]; M[5][j]=I[5][j];}
60 else if (D[1][j]>match) {M[1][j]=D[1][j] ; M[3][j]=D[3][j]; M[5][j]=D[5][j];}
61 else {M[1][j]=match ; M[3][j]=M[2][j-1]+id; M[5][j]=M[4][j-1]+1;}
63 P=I[0]; I[0]=I[1]; I[1]=P;
64 P=I[2]; I[2]=I[3]; I[3]=P;
65 P=I[4]; I[4]=I[5]; I[5]=P;
67 P=D[0]; D[0]=D[1]; D[1]=P;
68 P=D[2]; D[2]=D[3]; D[3]=P;
69 P=D[4]; D[4]=D[5]; D[5]=P;
71 P=M[0]; M[0]=M[1]; M[1]=P;
72 P=M[2]; M[2]=M[3]; M[3]=P;
73 P=M[4]; M[4]=M[5]; M[5]=P;
79 if ( strstr (comp_mode, "sim2"))
82 score=(l==0)?0:(M[2][l2]*100)/l;
84 else if ( strstr (comp_mode, "sim3"))
87 score=(l==0)?0:(M[2][l2]*100)/l;
89 else if ( strstr (comp_mode, "cov"))
92 score=(l==0)?0:((M[4][l2]*100)/l);
98 score=(l==0)?0:(M[2][l2]*100)/l;
108 int test_pair_wise (Alignment *A, int *ns, int **l_s, Constraint_list *CL)
111 char buf[VERY_LONG_STRING];
114 l0=strlen (A->seq_al[l_s[0][0]]);
115 l1=strlen (A->seq_al[l_s[1][0]]);
118 gap=generate_null(l1-n);
119 for (a=0;a<ns[0]; a++)
121 seq=A->seq_al[l_s[0][a]];
122 sprintf (buf, "%s%s",seq, gap);
123 sprintf (seq, "%s", buf);
126 gap=generate_null(l0-n);
128 for (a=0;a<ns[1]; a++)
130 seq=A->seq_al[l_s[1][a]];
131 sprintf (buf, "%s%s",seq, gap);
132 sprintf (seq, "%s", buf);
137 A->len_aln=strlen (A->seq_al[l_s[0][0]]);
138 A->score=A->score_aln=100;
142 int idscore_pair_wise (Alignment *A,int*ns, int **l_s,Constraint_list *CL)
145 A->score_aln=A->score=idscore_pairseq (A->seq_al[l_s[0][0]], A->seq_al[l_s[1][0]], CL->gop, CL->gep,CL->M, "sim3");
148 int dp_max (int *trace, int n, ...);
149 int dp_max (int *trace, int n, ...)
152 int a, v, t, best_v=0;
177 int is_tied (int *trace, int n, ...);
178 int is_tied(int *trace, int n, ...)
181 int a, v, t, best_v=0;
210 if (v==best_v && trace[0]!=t)
217 void display_mat (int **M, int l1, int l2, char *title);
218 void display_mat (int **M, int l1, int l2, char *title)
222 fprintf ( stdout, "\n\nTitle %s\n", title);
223 for ( a=0; a<=l1; a++)
225 fprintf ( stdout, "\n");
226 for ( b=0; b<=l2; b++)
227 fprintf ( stdout, "%3d ", M[a][b]);
230 int glocal_pair_wise ( Alignment *A, int *ns, int **l_s, Constraint_list *CL)
233 int i,j, l1, l2, n, sub, trace,ntrace, a, b, c, score;
234 int gop,rgop,tgop, gep, unmatch;
235 int M1, M2, I1, D1, LEN;
236 char **al, *char_buf, **aln;
240 l1=strlen (A->seq_al[l_s[0][0]]);
241 l2=strlen (A->seq_al[l_s[1][0]]);
244 M1=n++;D1=n++;I1=n++;M2=n++;
245 t=declare_arrayN(3, sizeof (int),n, l1+1, l2+1);
246 m=declare_arrayN(3, sizeof (int),n, l1+1, l2+1);
254 pos0=aln2pos_simple ( A,-1, ns, l_s);
257 for (j=1; j<=l2; j++)
266 for (i=1; i<=l1; i++)
272 for ( j=1; j<=l2; j++)
274 rgop=(i==l1 || j==1)?0:gop;
276 sub=(CL->get_dp_cost) (A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],j-1,CL);
277 m[M1][i][j]=dp_max (&trace,4,M1,m[M1][i-1][j-1],I1, m[I1][i-1][j-1],D1,m[D1][i-1][j-1],M2,m[M2][i-1][j-1])+sub;
280 m[D1][i][j]=dp_max (&trace,3, M1,m[M1][i][j-1]+rgop,D1, m[D1][i][j-1]+gep, M2, m[M2][i][j-1]);
283 m[I1][i][j]=dp_max (&trace,3, M1,m[M1][i-1][j]+rgop, I1, m[I1][i-1][j]+gep, M2, m[M2][i-1][j]);
286 m[M2][i][j]=dp_max (&trace,4,M1,m[M1][i-1][j-1]+tgop,I1, m[I1][i-1][j-1]+tgop,D1,m[D1][i-1][j-1]+tgop,M2,m[M2][i-1][j-1])+unmatch;
292 score=dp_max (&trace,4, M1,m[M1][l1][l2],D1,m[D1][l1][l2],I1, m[I1][l1][l2],M2,m[M2][l1][l2]);
294 al=declare_char (2, l1+l2+1);
297 trace=t[trace][i][j];
298 while (!(i==0 &&j==0))
301 ntrace=t[trace][i][j];
343 else if ( trace == I1)
354 invert_list_char ( al[0], LEN);
355 invert_list_char ( al[1], LEN);
356 if ( A->declared_len<=LEN)A=realloc_aln2 ( A,A->max_n_seq, 2*LEN);
359 char_buf= vcalloc (LEN+1, sizeof (char));
360 for ( c=0; c< 2; c++)
362 for ( a=0; a< ns[c]; a++)
365 for ( b=0; b< LEN; b++)
368 char_buf[b]=aln[l_s[c][a]][ch++];
373 sprintf (aln[l_s[c][a]],"%s", char_buf);
380 free_arrayN((void *)m, 3);
381 free_arrayN((void *)t, 3);
387 int ** aln2local_penalties (Alignment *A, int n, int *ls, Constraint_list *CL, int **lg);
388 int ** aln2local_penalties (Alignment *A, int n, int *ls, Constraint_list *CL, int **lg)
390 //adapted from gap_count in MAFFT V 5.5
393 int open=3, close=4, gap=5;
398 l=strlen (A->seq_al[ls[0]]);
402 lg=declare_int (6, l);
405 if ( read_array_size_new (lg[0])<l)
408 lg=declare_int (6, l);
416 c2=A->seq_al[ls[s]][p];
418 if (c1!='-' && c2=='-')lg[open][p]++;
419 if (c1=='-' && c2!='-')lg[close][p]++;
420 if ( c1=='-')lg[gap][p]++;
433 lg[GOP][p]=0.5*(1-(go/nn))*gop;
434 lg[GCP][p]=0.5*(1-(gc/nn))*gop;
435 //Checked locacal gep => gives low quality results
436 lg[GEP][p]=gep;//(1-((float)lg[gap][p]/(float)n))*gep;
437 lg[open][p]=lg[close][p]=lg[gap][p]=0;
443 int free_gotoh_pair_wise_lgp()
445 return gotoh_pair_wise_lgp (NULL, NULL, NULL, NULL);
447 int gotoh_pair_wise_lgp ( Alignment *A, int *ns, int **l_s, Constraint_list *CL)
449 int i,j, li, lj, n, sub, trace,ntrace, a, b, c, score;
452 char **al, *char_buf, **aln;
456 int gop[2], gcp[2], gep[2];
457 static int ***gpl, ***t, ***m;
458 static int max_li, max_lj;
462 //gotoh_pair_wise ( A, ns, l_s,CL);
463 //ungap_sub_aln (A, ns[0], l_s[0]);
464 //ungap_sub_aln (A, ns[1], l_s[1]);
468 free_arrayN((void**)gpl, 3);
469 free_arrayN((void**)t, 3);
470 free_arrayN((void**)m, 3);
478 li=strlen (A->seq_al[l_s[I][0]]);
479 lj=strlen (A->seq_al[l_s[J][0]]);
481 if ( !gpl)gpl=vcalloc ( 2, sizeof (int**));
482 gpl[I]=aln2local_penalties (A,ns[I], l_s[I], CL,gpl[I]);
483 gpl[J]=aln2local_penalties (A,ns[J], l_s[J], CL,gpl[J]);
487 M1=n++;D1=n++;I1=n++;
489 if ( li>max_li ||lj>max_lj )
491 free_arrayN((void**)t, 3);
492 free_arrayN((void**)m, 3);
497 t=declare_arrayN(3, sizeof (int),n, max_li+1, max_lj+1);
498 m=declare_arrayN(3, sizeof (int),n, max_li+1, max_lj+1);
501 pos0=aln2pos_simple ( A,-1, ns, l_s);
503 //Compatibility with Macro
507 for (j=1; j<=lj; j++)
509 gep[J]=gpl[J][GEP][j-1];
510 m[D1][0][j]=gep[J]*j;
511 m[I1][0][j]=m[D1][0][j]-1;
512 m[M1][0][j]=m[D1][0][j]-1;
515 //D1: gap in sequence I
516 //I1: gap in sequence J
519 for (i=1; i<=li; i++)
521 gep[I]=gpl[I][GEP][i-1];
522 gop[I]=gpl[I][GOP][i-1];
523 gcp[I]=gpl[I][GCP][i-1];
525 m[I1][i][0]=i*gep[I];
526 m[D1][i][0]= m[I1][i][0]-1;
527 m[M1][i][0]= m[I1][i][0]-1;
531 gop[I]=(i==1 || i==li )?0:gop[I];
532 gcp[I]=(i==1 || i==li )?0:gcp[I];
535 for ( j=1; j<=lj; j++)
538 gep[J]=gpl[J][GEP][j-1];
539 gop[J]=gpl[J][GOP][j-1];
540 gcp[J]=gpl[J][GCP][j-1];
542 //gep[J]=gep[I]=(gep[J]+gep[I])/2;
543 //gop[J]=gop[I]=(gop[J]+gop[I])/2;
544 //gcp[J]=gcp[I]=(gcp[J]+gcp[I])/2;
547 gop[J]=(j==1 || j==lj )?0:gop[J];
548 gcp[J]=(j==1 || j==lj )?0:gcp[J];
551 //sub=(CL->get_dp_cost) (A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],j-1,CL);
552 sub=TC_SCORE((i-1), (j-1));
554 m[M1][i][j]=dp_max (&trace,3,M1,m[M1][i-1][j-1],I1, m[I1][i-1][j-1]+gcp[I],D1,m[D1][i-1][j-1]+gcp[J])+sub;
558 m[D1][i][j]=dp_max (&trace,2, M1,m[M1][i][j-1]+gop[J]+gep[J],D1, m[D1][i][j-1]+gep[J]);
562 m[I1][i][j]=dp_max (&trace,2, M1,m[M1][i-1][j]+gop[I]+gep[I],I1, m[I1][i-1][j]+gep[I]);
568 score=dp_max (&trace,3, M1,m[M1][li][lj],D1,m[D1][li][lj],I1, m[I1][li][lj]);
571 al=declare_char (2, li+lj);
574 trace=t[trace][i][j];
575 while (!(i==0 &&j==0))
578 ntrace=t[trace][i][j];
609 else if ( trace == I1)
620 invert_list_char ( al[0], LEN);
621 invert_list_char ( al[1], LEN);
622 if ( A->declared_len<=LEN)A=realloc_aln ( A,2*LEN+1);
625 char_buf= vcalloc (LEN+1, sizeof (char));
626 for ( c=0; c< 2; c++)
628 for ( a=0; a< ns[c]; a++)
631 for ( b=0; b< LEN; b++)
634 char_buf[b]=aln[l_s[c][a]][ch++];
639 sprintf (aln[l_s[c][a]],"%s", char_buf);
651 /*******************************************************************************/
654 /* makes DP between the the ns[0] sequences and the ns[1] */
656 /* for MODE, see the function get_dp_cost */
657 /*******************************************************************************/
658 int glocal2_pair_wise (Alignment *IN,int*ns, int **ls,Constraint_list *CL)
664 buf=vcalloc (1000, sizeof (char));
665 seq=vcalloc (1000, sizeof (char));
667 A=copy_aln (IN,NULL);
668 L=copy_aln (IN,NULL);
669 R=copy_aln (IN,NULL);
671 gotoh_pair_wise_sw (A, ns, ls, CL);
676 for (b=0; b<ns[a]; b++)
679 sprintf ( seq,"%s", IN->seq_al[s]);
681 seq[A->order[s][2]]='\0';
682 sprintf (L->seq_al[s], "%s", seq);
683 sprintf (R->seq_al[s], "%s", seq+A->order[s][3]+1);
687 print_sub_aln (A, ns, ls);
688 gotoh_pair_wise(L, ns, ls, CL);
689 print_sub_aln (L, ns, ls);
690 gotoh_pair_wise(R, ns, ls, CL);
691 print_sub_aln (R, ns, ls);
693 IN=realloc_aln (IN, A->len_aln+L->len_aln+R->len_aln+1);
696 for (b=0; b<ns[a]; b++)
699 sprintf (IN->seq_al[s], "%s%s%s",L->seq_al[s], A->seq_al[s], R->seq_al[s]);
702 IN->len_aln=strlen (IN->seq_al[s]);
704 print_sub_aln (IN, ns, ls);
705 vfree (seq); vfree (buf);
706 free_aln (A); free_aln (L);free_aln (R);
707 return IN->score_aln;
711 int gotoh_pair_wise (Alignment *A,int*ns, int **l_s,Constraint_list *CL)
713 /*******************************************************************************/
714 /* NEEDLEMAN AND WUNSCH (GOTOH) */
716 /* makes DP between the the ns[0] sequences and the ns[1] */
718 /* for MODE, see the function get_dp_cost */
719 /*******************************************************************************/
722 /*TREATMENT OF THE TERMINAL GAP PENALTIES*/
723 /*TG_MODE=0---> gop and gep*/
724 /*TG_MODE=1---> --- gep*/
725 /*TG_MODE=2---> --- ---*/
732 /*VARIANLES FOR THE MULTIPLE SEQUENCE ALIGNMENT*/
750 /*trace back variables */
751 FILE *long_trace=NULL;
752 TRACE_TYPE *buf_trace=NULL;
753 static TRACE_TYPE **trace;
756 int long_trace_flag=0;
758 /********Prepare penalties*******/
762 maximise=CL->maximise;
765 /********************************/
766 /*CLEAN UP AFTER USE*/
776 /*DO MEMORY ALLOCATION FOR DP*/
778 lenal[0]=strlen (A->seq_al[l_s[0][0]]);
779 lenal[1]=strlen (A->seq_al[l_s[1][0]]);
780 len= MAX(lenal[0],lenal[1])+1;
781 buf_trace=vcalloc ( len, sizeof (TRACE_TYPE));
782 buffer=vcalloc ( 2*len, sizeof (char));
783 al=declare_char (2, 2*len);
785 char_buf= vcalloc (2*len, sizeof (char));
788 dd = vcalloc (len, sizeof (int));
791 cc = vcalloc (len, sizeof (int));
792 ddg=vcalloc (len, sizeof (int));
796 if ( len>=MAX_LEN_FOR_DP)
799 long_trace=vtmpfile();
804 dim=(trace==NULL)?0:read_size_int ( trace,sizeof (int*));
805 trace =realloc_int ( trace,dim,dim,MAX(0,len-dim), MAX(0,len-dim));
808 /*END OF MEMORY ALLOCATION*/
823 pos0=aln2pos_simple ( A,-1, ns, l_s);
827 tr=(long_trace_flag)?buf_trace:trace[0];
829 for ( j=1; j<=lenal[1]; j++)tr[j]=(TRACE_TYPE)-1;
830 if (long_trace_flag)fwrite (buf_trace, sizeof ( TRACE_TYPE),lenal[1]+1, long_trace);
833 t=(TG_MODE==0)?gop:0;
836 for (cc[0]=0,j=1; j<=lenal[1]; j++)
839 l_gop=(TG_MODE==0)?gop:0;
840 l_gep=(TG_MODE==2)?0:gep;
846 t=(TG_MODE==0)?gop:0;
848 for (i=1; i<=lenal[0];i++)
850 tr=(long_trace_flag)?buf_trace:trace[i];
853 l_gop=(TG_MODE==0)?gop:0;
854 l_gep=(TG_MODE==2)?0:gep;
864 for (eg=0,j=1; j<=lenal[1];j++)
867 sub=(CL->get_dp_cost) (A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],j-1,CL);
869 /*get the best Insertion*/
870 l_gop=(i==lenal[0] || i==1 )?((TG_MODE==0)?gop:0):gop;
871 l_gep=(i==lenal[0] || i==1)?((TG_MODE==2)?0:gep):gep;
874 if ( a_better_than_b ( e,c+l_gop, maximise))eg++;
876 e=best_of_a_b (e, c+l_gop, maximise)+l_gep;
878 /*Get the best deletion*/
879 l_gop=(j==lenal[1] || j==1)?((TG_MODE==0)?gop:0):gop;
880 l_gep=(j==lenal[1] || j==1)?((TG_MODE==2)?0:gep):gep;
883 if ( a_better_than_b ( dd[j], cc[j]+l_gop, maximise))ddg[j]++;
885 dd[j]=best_of_a_b( dd[j], cc[j]+l_gop,maximise)+l_gep;
889 c=best_int(3,maximise,&fop, e, s+sub,dd[j]);
890 /*Chose Substitution for tie breaking*/
891 if ( fop==0 && (s+sub)==e)fop=1;
892 else if ( fop==2 && (s+sub)==dd[j])fop=1;
893 /*Chose Deletion for tie breaking*/
894 else if ( fop==2 && e==dd[j])fop=1;
902 {tr[j]=(TRACE_TYPE)fop*eg;
905 {tr[j]=(TRACE_TYPE)fop*ddg[j];
908 {tr[j]=(TRACE_TYPE)0;
914 fwrite ( buf_trace, sizeof (TRACE_TYPE), lenal[1]+1, long_trace);
925 while (i>=0 && j>=0 && ((i+j)!=0))
931 else if ( j==0 && i==0)
937 fseek ( long_trace, sizeof (TRACE_TYPE)*((lenal[1]+1)*(i)+j),SEEK_SET);
938 fread ( &k, sizeof (TRACE_TYPE), 1, long_trace);
959 for ( a=0; a< k; a++)
983 invert_list_char ( al[0], LEN);
984 invert_list_char ( al[1], LEN);
985 if ( A->declared_len<=LEN)A=realloc_aln2 ( A,A->max_n_seq, 2*LEN);
988 for ( c=0; c< 2; c++)
990 for ( a=0; a< ns[c]; a++)
993 for ( b=0; b< LEN; b++)
996 char_buf[b]=aln[l_s[c][a]][ch++];
1001 sprintf (aln[l_s[c][a]],"%s", char_buf);
1007 A->nseq=ns[0]+ns[1];
1016 free_char ( al, -1);
1017 free_int (pos0, -1);
1018 if ( long_trace_flag)fclose (long_trace);
1026 int get_transition_cost (Alignment *A, int **posi, int ni, int *li, int i, int **posj, int nj, int *lj, int j,Constraint_list *CL);
1027 int gotoh_pair_wise_lgp_sticky ( Alignment *A, int *ns, int **l_s, Constraint_list *CL)
1029 int i,j, li, lj, n, sub, trace,ntrace, a, b, c, score;
1031 int M1, I1, D1, LEN;
1032 char **al, *char_buf, **aln;
1035 int gop[2], gcp[2], gep[2];
1036 static int ***gpl, ***t, ***m;
1037 static int max_li, max_lj;
1041 //gotoh_pair_wise ( A, ns, l_s,CL);
1042 //ungap_sub_aln (A, ns[0], l_s[0]);
1043 //ungap_sub_aln (A, ns[1], l_s[1]);
1048 li=strlen (A->seq_al[l_s[I][0]]);
1049 lj=strlen (A->seq_al[l_s[J][0]]);
1051 if ( !gpl)gpl=vcalloc ( 2, sizeof (int**));
1052 gpl[I]=aln2local_penalties (A,ns[I], l_s[I], CL,gpl[I]);
1053 gpl[J]=aln2local_penalties (A,ns[J], l_s[J], CL,gpl[J]);
1057 M1=n++;D1=n++;I1=n++;
1059 if ( li>max_li ||lj>max_lj )
1061 free_arrayN((void**)t, 3);
1062 free_arrayN((void**)m, 3);
1067 t=declare_arrayN(3, sizeof (int),n, max_li+1, max_lj+1);
1068 m=declare_arrayN(3, sizeof (int),n, max_li+1, max_lj+1);
1071 pos0=aln2pos_simple ( A,-1, ns, l_s);
1074 for (j=1; j<=lj; j++)
1076 gep[J]=gpl[J][GEP][j-1];
1077 m[D1][0][j]=gep[J]*j;
1078 m[I1][0][j]=m[D1][0][j]-1;
1079 m[M1][0][j]=m[D1][0][j]-1;
1082 //D1: gap in sequence I
1083 //I1: gap in sequence J
1086 for (i=1; i<=li; i++)
1088 gep[I]=gpl[I][GEP][i-1];
1089 gop[I]=gpl[I][GOP][i-1];
1090 gcp[I]=gpl[I][GCP][i-1];
1092 m[I1][i][0]=i*gep[I];
1093 m[D1][i][0]= m[I1][i][0]-1;
1094 m[M1][i][0]= m[I1][i][0]-1;
1098 gop[I]=(i==1 || i==li )?0:gop[I];
1099 gcp[I]=(i==1 || i==li )?0:gcp[I];
1102 for ( j=1; j<=lj; j++)
1106 gep[J]=gpl[J][GEP][j-1];
1107 gop[J]=gpl[J][GOP][j-1];
1108 gcp[J]=gpl[J][GCP][j-1];
1110 //gep[J]=gep[I]=(gep[J]+gep[I])/2;
1111 //gop[J]=gop[I]=(gop[J]+gop[I])/2;
1112 //gcp[J]=gcp[I]=(gcp[J]+gcp[I])/2;
1115 gop[J]=(j==1 || j==lj )?0:gop[J];
1116 gcp[J]=(j==1 || j==lj )?0:gcp[J];
1119 sub=(CL->get_dp_cost) (A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],j-1,CL);
1120 transition=get_transition_cost (A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],j-1,CL);
1122 m[M1][i][j]=dp_max (&trace,3,M1,m[M1][i-1][j-1]+transition,I1, m[I1][i-1][j-1]+gcp[I],D1,m[D1][i-1][j-1]+gcp[J])+sub;
1126 m[D1][i][j]=dp_max (&trace,2, M1,m[M1][i][j-1]+gop[J]+gep[J],D1, m[D1][i][j-1]+gep[J]);
1130 m[I1][i][j]=dp_max (&trace,2, M1,m[M1][i-1][j]+gop[I]+gep[I],I1, m[I1][i-1][j]+gep[I]);
1136 score=dp_max (&trace,3, M1,m[M1][li][lj],D1,m[D1][li][lj],I1, m[I1][li][lj]);
1139 al=declare_char (2, li+lj);
1142 trace=t[trace][i][j];
1143 while (!(i==0 &&j==0))
1146 ntrace=t[trace][i][j];
1163 else if ( trace==M1)
1170 else if ( trace==D1)
1177 else if ( trace == I1)
1188 invert_list_char ( al[0], LEN);
1189 invert_list_char ( al[1], LEN);
1190 if ( A->declared_len<=LEN)A=realloc_aln ( A,2*LEN+1);
1193 char_buf= vcalloc (LEN+1, sizeof (char));
1194 for ( c=0; c< 2; c++)
1196 for ( a=0; a< ns[c]; a++)
1199 for ( b=0; b< LEN; b++)
1202 char_buf[b]=aln[l_s[c][a]][ch++];
1207 sprintf (aln[l_s[c][a]],"%s", char_buf);
1213 A->nseq=ns[0]+ns[1];
1216 free_int (pos0, -1);
1219 int get_transition_cost (Alignment *A, int **posi, int ni, int *li, int i, int **posj, int nj, int *lj, int j,Constraint_list *CL)
1221 /*counts the number of identical transitions between position i-1, i and j-1..j*/
1226 if (i==0 || j==0)return 0;
1228 for (a=0; a<ni; a++)
1231 if (posi[s][i]<0 || posi[s][i-1]<0)continue;
1232 if (S->seq[li[a]][i-1]==S->seq[li[a]][i-1])t++;
1235 for (a=0; a<nj; a++)
1238 if (posj[s][j]<0 || posj[s][j-1]<0)continue;
1239 if (S->seq[li[a]][j-1]==S->seq[li[a]][j-1])t++;
1242 t=(t*10)/(float)(ni+nj);
1245 /*******************************************************************************/
1246 /* idscore_pairseq: measure the % id without delivering thze aln*/
1248 /* makes DP between the the ns[0] sequences and the ns[1] */
1250 /* for MODE, see the function get_dp_cost */
1251 /*******************************************************************************/
1252 int cl2pair_list ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in, int mode, int ndiag);
1253 int cl2pair_list_ref ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
1254 int cl2pair_list_ecf ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
1255 int cl2pair_list_diag ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in, int add);
1256 int cl2list_borders (Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
1257 int cl2diag_cap (Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in); //add one element at the end of each segment so that they can be joined
1258 int** cl2sorted_diagonals ( Alignment *A, int *ns, int **ls, Constraint_list *CL);
1259 int** cl2sorted_diagonals_mat ( Alignment *A, int *ns, int **ls, Constraint_list *CL);
1260 int** cl2sorted_diagonals_cs ( Alignment *A, int *ns, int **ls, Constraint_list *CL);
1261 int list2nodup_list (Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
1262 int fill_matrix ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
1263 int cl2pair_list ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in, int mode, int ndiag)
1268 free_int (list_in[0], -1); list_in[0]=NULL;
1272 cl2list_borders(A, ns, ls, CL, list_in, n_in);
1275 v=cl2pair_list_ref (A, ns, ls, CL, list_in, n_in);
1277 v=cl2pair_list_ecl (A, ns, ls, CL, list_in, n_in);
1279 v=cl2pair_list_diag (A, ns, ls, CL, list_in, n_in,ndiag); //add diagonals
1281 cl2diag_cap (A, ns, ls, CL, list_in, n_in);
1282 //fill_matrix (A, ns, ls, CL, list_in, n_in);//Fill matrix with 0s
1283 sort_list_int (list_in[0],7, 1, 0, n_in[0]-1);
1284 list2nodup_list (A, ns, ls, CL, list_in, n_in);
1288 int fill_matrix( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list, int *n)
1290 int a, b, l1, l2, n2=0;
1295 pos=aln2pos_simple ( A,-1, ns, ls);
1296 l1=strlen (A->seq_al[ls[0][0]]);
1297 l2=strlen (A->seq_al[ls[1][0]]);
1299 max_n=read_array_size (list[0], sizeof (int));
1301 for (a=0; a<=l1; a++)
1302 for (b=0; b<=l2; b++)
1305 score=(a==0 || b==0)?0:slow_get_dp_cost (A, pos, ns[0], ls[0],a-1, pos, ns[1], ls[1], b-1, CL);
1306 if ( score>0 && a!=0 && b!=0 && a!=l1 && b!=l2)
1308 if (n[0]==max_n){max_n+=1000;list[0]=vrealloc (list[0], max_n*sizeof (int*));}
1309 if (!list[0][n[0]])list[0][n[0]]=vcalloc (7, sizeof (int));
1312 list[0][n[0]][3]=(l1-a)+b;
1313 list[0][n[0]][2]=score;
1314 if ( a!=0 && b!=0 && a!=l1 && b!=l2)
1325 int list2nodup_list ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
1335 for (b=a=1; a<n; a++)
1337 if (list[a][0]==list[b-1][0] && list[a][1]==list[b-1][1])
1339 //HERE ("Duplicate");
1340 list[b-1][2]=MAX(list[b-1][2],list[a][2]);
1344 for (c=0; c<4; c++)list[b][c]=list[a][c];
1353 int** cl2sorted_diagonals ( Alignment *A, int *ns, int **ls, Constraint_list *CL)
1355 if ( CL && CL->L)return cl2sorted_diagonals_cs ( A, ns, ls, CL);
1356 else return cl2sorted_diagonals_mat ( A, ns, ls, CL);
1360 static char **warray;
1361 int cmp_word ( const int**a, const int**b);
1362 int ** seq2index_list ( Sequence *S, int k);
1363 int** cl2sorted_diagonals_mat ( Alignment *A, int *ns, int **ls, Constraint_list *CL)
1366 int a,b,c,d, comp, k, l1, l2, ndiag;
1370 static char **alp, alps=5;
1375 alp=make_group_aa_upgma ("blosum62mt",alps);
1379 l1=strlen (A->seq_al[ls[0][0]]);
1380 buf1=vcalloc ( l1+1, sizeof (char));
1381 l2=strlen (A->seq_al[ls[1][0]]);
1382 buf2=vcalloc ( l2+1, sizeof (char));
1385 diag=declare_int (ndiag+3,2);
1386 for (a=0; a<=ndiag; a++) diag[a][0]=a;
1387 vfree (diag[ndiag+1]);
1390 for ( a=0; a<ns[0]; a++)
1392 sprintf (buf1, A->seq_al[ls[0][a]]);
1393 lower_string (buf1);
1394 string_convert (buf1, alps, alp);
1395 for (b=0; b<ns[1]; b++)
1397 sprintf (buf2, A->seq_al[ls[1][b]]);
1398 lower_string (buf2);
1399 string_convert (buf2, alps, alp);
1400 for (c=0; c<l1-k; c++)
1402 if ( strnrchr (buf1+c, '-', k))continue;
1403 for (d=0; d<l2-k; d++)
1405 if ( strnrchr (buf2+d, '-', k))continue;
1406 comp=strncmp (buf1+c,buf2+d, k);
1407 diag[l1-c+d][1]+=(comp==0)?1:0;
1414 max_len=MAX(l1, l2);
1415 for (a=1; a<ndiag; a++)
1417 int l, d, p1_0, p1_1, p2_0, p2_1;
1420 if (d<=l1)l=MIN(d,l2);
1421 else l=MIN(((l1+l2)-d),l1);
1423 diag[a][1]=(diag[a][1]*1000)/max_len;
1424 diag[a][1]=(float)((float)diag[a][1]*(((float)max_len)/(float)l));
1427 sort_int_inv (diag, 2, 1,0,ndiag);
1438 int ** seq2index_list ( Sequence *S, int k)
1440 int **list,**mlist=NULL;
1441 int a, b,c,ml, n, l, e, s, max=0, nm=0;
1445 for (ml=0,a=0; a<S->nseq; a++)ml+=strlen (S->seq[a]);
1446 list=declare_int (ml+1, 2);
1448 for (n=0,a=0; a<S->nseq; a++)
1450 l=strlen (S->seq[a])-k;
1451 for ( b=0; b<l; b++, n++)
1461 qsort (list, n, sizeof (long**), (int(*)(const void*,const void*))(cmp_word));
1466 for (a=0; a<=n; a++)
1469 if (!cw ||a==n|| strncmp (warray[list[a][0]]+list[a][1],cw, k)!=0)
1471 if (a<n)cw=warray[list[a][0]]+list[a][1];
1472 for (b=s; b<a-1; b++)
1473 for (c=b+1; c<a; c++)
1476 if (list[b][0]<list[c][0])
1490 if (s1==s2)continue;
1493 if (nm>=max){max+=1000; mlist=vrealloc (mlist, max*sizeof (int*));}
1494 mlist[nm]=vcalloc (4, sizeof (int));
1506 if (nm>=max){max+=1000;mlist=vrealloc (mlist, max*sizeof (int));}
1507 sort_list_int ( mlist,4,1, 0, nm-1);
1510 int cmp_word ( const int**a, const int**b)
1514 c=strncmp (warray[a[0][0]]+a[0][1], warray[b[0][0]]+b[0][1], kword);
1522 if ( a[0][c]>b[0][c])return 1;
1523 else if (a[0][c]<b[0][c])return -1;
1529 int** cl2sorted_diagonals_cs ( Alignment *A, int *ns, int **ls, Constraint_list *CL)
1531 int p1, p2, s1, s2, r1, r2;
1539 pos=aln2pos_simple ( A,-1, ns, ls);
1541 l1=strlen (A->seq_al[ls[0][0]]);
1542 l2=strlen (A->seq_al[ls[1][0]]);
1546 CL=index_res_constraint_list (CL, CL->weight_field);
1548 diag=declare_int (ndiag+3, 2);
1550 for (a=1; a<=ndiag; a++)diag[a][0]=a;
1552 for (p1=1; p1<=l1; p1++)
1554 for (p2=1; p2<=l2; p2++)
1556 for (a=0; a<ns[0]; a++)
1561 for (b=0; b<ns[1]; b++)
1568 diag[diag_i][1]+=residue_pair_extended_list_raw (CL,s1, r1-1,s2, r2-1);
1574 sort_int_inv (diag, 2, 1,0,ndiag);
1576 vfree (diag[ndiag+1]);
1581 int** cl2sorted_diagonals_cs_old2 ( Alignment *A, int *ns, int **ls, Constraint_list *CL)
1583 int p1, p2, si, s, r, t_s, t_r;
1585 int *sl2, **pos,**inv_pos;
1591 pos=aln2pos_simple ( A,-1, ns, ls);
1592 inv_pos=vcalloc ((CL->S)->nseq, sizeof (int*));
1593 for (a=0; a<ns[1]; a++)inv_pos[ls[1][a]] =seq2inv_pos(A->seq_al[ls[1][a]]);
1595 l1=strlen (A->seq_al[ls[0][0]]);
1596 l2=strlen (A->seq_al[ls[1][0]]);
1597 sl2=vcalloc ((CL->S)->nseq, sizeof (int));
1599 for (a=0;a<ns[1]; a++)sl2[ls[1][a]]=1;
1600 CL=index_res_constraint_list (CL, CL->weight_field);
1602 diag=declare_int (ndiag+3, 2);
1604 for (a=1; a<=ndiag; a++)diag[a][0]=a;
1605 for (p1=0; p1<=l1; p1++)
1607 for (si=0;p1>0 && si<ns[0]; si++)
1612 for (a=1; r>0 && a<CL->residue_index[s][r][0];a+=3)
1615 t_s=CL->residue_index[s][r][a];
1616 t_r=CL->residue_index[s][r][a+1];
1620 p2=inv_pos[t_s][t_r];
1622 diag[diag_i][1]+=(CL->get_dp_cost) (A, pos, ns[0], ls[0], p1-1, pos, ns[1], ls[1],p2-1,CL);
1627 max_len=MAX(l1, l2);
1628 for (a=1; a<ndiag; a++)
1633 if (d<=l1)l=MIN(d,l2);
1634 else l=MIN(((l1+l2)-d),l1);
1636 diag[a][1]/=max_len;
1637 diag[a][1]=(float)((float)diag[a][1]*(((float)max_len)/(float)l));
1639 sort_int_inv (diag, 2, 1,0,ndiag);
1641 vfree (diag[ndiag+1]);
1644 free_int (inv_pos, -1);
1649 int cl2list_borders (Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
1651 int maxlen, a,n, p1, p2, l1, l2;
1656 // printf("ICH BIN HIER\n");
1657 if ( list_in[0] && list_in[0][0]==0)
1658 return read_array_size (list, sizeof (int*));
1660 // printf("UND HIER IRGENDWIE AUCH\n");
1665 else maxlen=read_array_size (list, sizeof (int*));
1668 l1=strlen (A->seq_al[ls[0][0]]);
1669 l2=strlen (A->seq_al[ls[1][0]]);
1670 // pos=aln2pos_simple ( A,-1, ns, ls);
1672 for (p1=0; p1<=l1; p1++)
1674 if (p1==0 || p1==l1)
1676 for (p2=0; p2<=l2; p2++)
1678 if (n==maxlen){maxlen+=1000;list=vrealloc (list,maxlen*sizeof (int*));}
1679 if (!list[n])list[n]=vcalloc (7, sizeof (int));
1682 list[n][3]=(l1-(p1))+(p2);
1683 //list[n][2]=(p1==0||p2==0)?0:(CL->get_dp_cost) (A, pos, ns[0], ls[0], p1-1, pos, ns[1], ls[1],p2-1,CL);;
1684 list[n][2]=(CL->gep)*SCORE_K*p2;
1693 if (n==maxlen){maxlen+=1000;list=vrealloc (list,maxlen*sizeof (int*));}
1694 if (!list[n])list[n]=vcalloc (7, sizeof (int));
1697 list[n][3]=(l1-(p1))+(p2);
1698 //list[n][2]=(p1==0||p2==0)?0:(CL->get_dp_cost) (A, pos, ns[0], ls[0], p1-1, pos, ns[1], ls[1],p2-1,CL);;
1699 list[n][2]=(CL->gep)*SCORE_K*p1;
1704 // free_int (pos, -1);
1707 return read_array_size (list, sizeof (int*));
1710 int cl2diag_cap (Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
1713 int n, in, a, b, al1, al2;
1719 al1=strlen (A->seq_al[ls[0][0]]);
1720 al2=strlen (A->seq_al[ls[1][0]]);
1724 max_n=read_array_size (list, sizeof (int*));
1729 for (a=0; a< n; a++)
1732 list[a][3]=list[a][0];
1736 sort_list_int (list, 4, 1, 0, n-1);
1737 for (a=0; a< n; a++)
1740 list[a][3]=list[a][0];
1747 for (a=0; a<in; a++)
1749 int i, j, pi, pj, ni, nj;
1750 if (list[a][2]==0)continue;
1754 if (a==0){pi=-10;pj=-10;}
1755 else {pi=list[a-1][0];pj=list[a-1][1];}
1757 if (a==in-1){ni=-10; nj=-10;}
1758 else {ni=list[a+1][0]; nj=list[a+1][1];}
1762 else if ( i==pi || j==pj);
1763 else if ( i-pi!=1 || j-pj!=1)
1766 if (n==max_n){max_n+=1000;list=vrealloc (list, max_n*sizeof (int*));}
1767 if (!list[n])list[n]=vcalloc (7, sizeof (int));
1771 list[n][3]=list[a][3];
1777 if (i==al1 || j==al2);
1778 else if ( i==ni || j==nj);
1779 else if ( ni-i!=1 || nj-j!=1)
1782 if (n==max_n){max_n+=1000;list=vrealloc (list, max_n*sizeof (int*));}
1783 if (!list[n])list[n]=vcalloc (7, sizeof (int));
1787 list[n][3]=list[a][3];
1798 int cl2pair_list_diag_mat ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in, int add );
1799 int cl2pair_list_diag_cl ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in, int add );
1800 int cl2pair_list_diag ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in, int add )
1802 if (CL->L)return cl2pair_list_diag_cl (A, ns, ls, CL, list_in, n_in, add);
1803 else return cl2pair_list_diag_mat (A, ns, ls, CL, list_in, n_in, add);
1805 int cl2pair_list_diag_mat ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in, int add )
1817 free_int (pos, -1);pos=NULL;
1818 free_int (diag, -1);diag=NULL;
1819 //free_int (list_in[0], -1); list_in[0]=NULL;
1825 pos=aln2pos_simple ( A,-1, ns, ls);
1826 diag=cl2sorted_diagonals (A,ns,ls,CL);
1831 max_n=read_array_size (list, sizeof (int**));
1835 l1=strlen (A->seq_al[ls[0][0]]);
1836 l2=strlen (A->seq_al[ls[1][0]]);
1841 while ( diag[d] && diag[d][1]==-1)d++;
1847 while (diag[add++]);
1850 HERE ("Add %d diagonals, starts %d N=%d", add, d, n);
1852 for (d=0; d<add && diag[d]; d++)
1858 HERE ("\t S=%d", diag[d][1]);
1860 p1_0=MAX(0,l1-diag[d][0]);
1861 p2_0=MAX(0,diag[d][0]-l1);
1864 for (p1=p1_0, p2=p2_0; p1<=l1 && p2<=l2; p1++,p2++)
1866 if (!BORDER(p1,l1, p2,l2) )
1869 if (n==max_n){max_n+=1000;list=vrealloc (list, max_n*sizeof (int*));}
1870 if (!list[n])list[n]=vcalloc (7, sizeof (int));
1874 list[n][3]=(l1-(p1))+(p2);
1875 list[n][2]=(CL->get_dp_cost) (A, pos, ns[0], ls[0], p1-1, pos, ns[1], ls[1],p2-1,CL);
1880 HERE ("Addition Finished n=%d", n);
1882 sort_list_int (list,4, 1, 0, n-1);
1886 HERE ("\nN=%d r=%.3f [l1=%d l2=%d]", n, (float)n/(float)(l1*l2), l1, l2);
1890 int cl2pair_list_diag_cl ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in, int add )
1902 free_int (pos, -1);pos=NULL;
1903 free_int (diag, -1);diag=NULL;
1904 //free_int (list_in[0], -1); list_in[0]=NULL;
1910 pos=aln2pos_simple ( A,-1, ns, ls);
1911 diag=cl2sorted_diagonals (A,ns,ls,CL);
1916 max_n=read_array_size (list, sizeof (int**));
1920 l1=strlen (A->seq_al[ls[0][0]]);
1921 l2=strlen (A->seq_al[ls[1][0]]);
1922 if ( add==0)add=l1+l2;
1924 while ( diag[d] && diag[d][1]==-1)d++;
1925 HERE ("Add %d diagonals, starts %d N=%d", add, d, n);
1927 for (d; d<add && diag[d]; d++)
1930 if (diag[d][1]==0)continue;
1933 HERE ("\t S=%d", diag[d][1]);
1935 p1_0=MAX(0,l1-diag[d][0]);
1936 p2_0=MAX(0,diag[d][0]-l1);
1939 for (p1=p1_0, p2=p2_0; p1<=l1 && p2<=l2; p1++,p2++)
1941 if (!BORDER(p1,l1, p2,l2) && (score=(CL->get_dp_cost) (A, pos, ns[0], ls[0], p1-1, pos, ns[1], ls[1],p2-1,CL))!=0)
1944 if (n==max_n){max_n+=1000;list=vrealloc (list, max_n*sizeof (int*));}
1945 if (!list[n])list[n]=vcalloc (7, sizeof (int));
1949 list[n][3]=(l1-(p1))+(p2);
1955 HERE ("Addition Finished n=%d", n);
1957 sort_list_int (list,4, 1, 0, n-1);
1961 HERE ("\nN=%d r=%.3f [l1=%d l2=%d]", n, (float)n/(float)(l1*l2), l1, l2);
1965 int cl2pair_list_ecl_norm ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
1966 int cl2pair_list_ecl_raw ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
1967 int cl2pair_list_ecl_pc ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
1968 int cl2pair_list_ecl_noext_raw ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
1969 int cl2pair_list_ecl_rna2 ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
1971 int cl2pair_list_ecl_rawquad ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
1972 int cl2pair_list_ecl ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
1976 if ( mode==1)return cl2pair_list_ecl_norm (A, ns, ls, CL, list_in, n_in);
1977 else if ( mode==2)return cl2pair_list_ecl_raw (A, ns, ls, CL, list_in, n_in);
1978 else if ( mode==3)return cl2pair_list_ecl_rawquad (A, ns, ls, CL, list_in, n_in);
1979 else if ( mode==4)return cl2pair_list_ecl_noext_raw (A, ns, ls, CL, list_in, n_in);
1980 else if ( mode==5)return cl2pair_list_ecl_pc (A, ns, ls, CL, list_in, n_in);
1982 int cl2pair_list_ecl_noext_raw ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
1984 int p1, p2, si, s, r, t_s2, t_r2, t_w2, n,n2;
1994 int *used_list, **used;
1995 int *sl2, **inv_pos;
2004 max_n=read_array_size (list, sizeof (int*));
2008 pos=aln2pos_simple ( A,-1, ns, ls);
2009 inv_pos=vcalloc ((CL->S)->nseq, sizeof (int*));
2010 for (a=0; a<ns[1]; a++)inv_pos[ls[1][a]] =seq2inv_pos(A->seq_al[ls[1][a]]);
2012 l1=strlen (A->seq_al[ls[0][0]]);
2013 l2=strlen (A->seq_al[ls[1][0]]);
2014 sl2=vcalloc ((CL->S)->nseq, sizeof (int));
2016 for (a=0;a<ns[1]; a++)sl2[ls[1][a]]=1;
2018 CL=index_res_constraint_list (CL, CL->weight_field);
2020 used=declare_int (l2+1,2);
2021 used_list=vcalloc (l2+1, sizeof (int));
2024 for (p1=0; p1<=l1; p1++)
2026 for (nused=0,si=0;p1>0 && si<ns[0]; si++)
2028 s=ls [0][si];r=pos[s][p1-1];
2029 for (a=1; r>0 && a<CL->residue_index[s][r][0];a+=3)
2031 t_s2=CL->residue_index[s][r][a];t_r2=CL->residue_index[s][r][a+1];t_w2=CL->residue_index[s][r][a+2];
2034 p2=inv_pos[t_s2][t_r2];
2036 if (!used[p2][1] && score>0)
2038 used_list[nused++]=p2;
2045 for (a=0; a<nused; a++)
2050 score=used[p2][0]*SCORE_K;
2051 used[p2][0]=used[p2][1]=0;
2053 if (n==max_n){max_n+=1000;list=vrealloc (list, max_n*sizeof (int*));}
2054 if (score>filter && p1!=0 && p2!=0 && p1!=l1 && p2!=l2)
2056 if (!list[n])list[n]=vcalloc (7, sizeof (int));
2060 list[n][3]=(l1-(p1))+(p2);
2070 free_int (inv_pos, -1);
2079 int cl2pair_list_ecl_raw ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
2081 int p1, p2, si, s, r, t_s, t_r,t_w, t_s2, t_r2, t_w2, n,tot;
2087 int set, raw_max,nscore, score, nused;
2088 int *used_list, **used;
2089 int *sl2, **inv_pos;
2091 int filter1=0, filter2=0, max=0;
2093 long tot_score=0, avg;
2099 max_n=read_array_size (list, sizeof (int*));
2102 pos=aln2pos_simple ( A,-1, ns, ls);
2103 inv_pos=vcalloc ((CL->S)->nseq, sizeof (int*));
2104 for (a=0; a<ns[1]; a++)inv_pos[ls[1][a]] =seq2inv_pos(A->seq_al[ls[1][a]]);
2106 l1=strlen (A->seq_al[ls[0][0]]);
2107 l2=strlen (A->seq_al[ls[1][0]]);
2108 sl2=vcalloc ((CL->S)->nseq, sizeof (int));
2109 nr=declare_int (2, MAX(l1,l2)+1);
2112 for (b=0; b<ns[0]; b++)
2113 if (!is_gap(A->seq_al[ls[0][b]][a]))nr[0][a+1]++;
2115 for (b=0; b<ns[1]; b++)
2116 if (!is_gap(A->seq_al[ls[1][b]][a]))nr[1][a+1]++;
2118 for (a=0;a<ns[1]; a++)sl2[ls[1][a]]=1;
2120 CL=index_res_constraint_list (CL, CL->weight_field);
2122 used=declare_int (l2+1,2);
2123 used_list=vcalloc (l2+1, sizeof (int));
2126 for (raw_max=0,p1=0; p1<=l1; p1++)
2128 for (tot=0,nused=0,si=0;p1>0 && si<ns[0]; si++)
2130 s=ls [0][si];r=pos[s][p1-1];
2131 for (a=1; r>0 && a<CL->residue_index[s][r][0];a+=3)
2133 t_s=CL->residue_index[s][r][a];t_r=CL->residue_index[s][r][a+1];t_w=CL->residue_index[s][r][a+2];
2134 for (b=0; b<CL->residue_index[t_s][t_r][0];)
2136 if (b==0){t_s2=t_s;t_r2=t_r;t_w2=t_w;b++;}
2139 t_s2=CL->residue_index[t_s][t_r][b];t_r2=CL->residue_index[t_s][t_r][b+1];t_w2=CL->residue_index[t_s][t_r][b+2];b+=3;
2144 p2=inv_pos[t_s2][t_r2];
2145 score=MIN(t_w,t_w2);
2146 if (score<filter1)score=0;
2147 //if ( score<3) continue;
2148 if (!used[p2][1] && score>0)
2150 used_list[nused++]=p2;
2159 //set the threshold to 1/2 of the best normalised score
2161 for (filter2=0,set=0,a=0; a<nused; a++)
2167 nscore=(score*100)/tot;
2168 if (set==0){filter2=nscore;set=1;}
2169 filter2=MAX(nscore,filter2);
2170 if ( score<0)HERE ("*********** %d", score);
2174 for (a=0; a<nused; a++)
2179 //score=used[p2][0];
2180 nscore=(used[p2][0]*100)/tot; //Normalized score used for filtering
2182 raw_max=MAX(score, raw_max);
2184 used[p2][0]=used[p2][1]=0;
2186 if (n==max_n){max_n+=1000;list=vrealloc (list, max_n*sizeof (int*));}
2187 if (nscore>=filter2 && p1!=0 && p2!=0 && p1!=l1 && p2!=l2)
2189 if (!list[n])list[n]=vcalloc (7, sizeof (int));
2193 list[n][3]=(l1-(p1))+(p2);
2201 avg=tot_score/new_n;
2203 //CL->gop=-1*avg*3;CL->gep=0;
2204 HERE ("FILTER: %d->%d [THR=%d]", max, n-n_in[0], filter2);
2208 free_int (inv_pos, -1);
2219 * Calculates scores for diagonal segments.
2221 * \param Alignment The sequences.
2222 * \param ns Number of sequences in each group
2223 * \param ls sequences in in groups (ls[0][x] sequences in group 1, ls[1][x] squences in group 2).
2224 * \param CL the constraint list
2225 * \param list_in the diagonals
2226 * \param n_in number of sequences?
2228 int cl2pair_list_ecl_pc ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
2230 int p1, p2, si, s, r, t_s, t_r,t_w, t_s2, t_r2, t_w2, n;
2237 int *sl2, **inv_pos;
2242 float nscore, score, tot, filter, avg=0, new=0;
2248 max_n=read_array_size (list, sizeof (int*));
2250 pos=aln2pos_simple ( A,-1, ns, ls);
2251 inv_pos=vcalloc ((CL->S)->nseq, sizeof (int*));
2252 for (a=0; a<ns[1]; a++)inv_pos[ls[1][a]] =seq2inv_pos(A->seq_al[ls[1][a]]);
2254 l1=strlen (A->seq_al[ls[0][0]]);
2255 l2=strlen (A->seq_al[ls[1][0]]);
2256 sl2=vcalloc ((CL->S)->nseq, sizeof (int));
2257 nr=declare_int (2, MAX(l1,l2)+1);
2260 for (b=0; b<ns[0]; b++)
2261 if (!is_gap(A->seq_al[ls[0][b]][a]))nr[0][a+1]++;
2263 for (b=0; b<ns[1]; b++)
2264 if (!is_gap(A->seq_al[ls[1][b]][a]))nr[1][a+1]++;
2266 for (a=0;a<ns[1]; a++)sl2[ls[1][a]]=1;
2268 CL=index_res_constraint_list (CL, CL->weight_field);
2270 used=declare_float (l2+1,2);
2271 used_list=vcalloc (l2+1, sizeof (int));
2274 for (p1=0; p1<=l1; p1++)
2277 for (tot=0,nused=0,si=0;p1>0 && si<ns[0]; si++)
2279 s=ls [0][si];r=pos[s][p1-1];
2280 for (a=1; r>0 && a<CL->residue_index[s][r][0];a+=3)
2282 t_s=CL->residue_index[s][r][a];t_r=CL->residue_index[s][r][a+1];t_w=CL->residue_index[s][r][a+2];
2283 for (b=0; b<CL->residue_index[t_s][t_r][0];)
2285 if (b==0){t_s2=t_s;t_r2=t_r;t_w2=t_w;b++;}
2288 t_s2=CL->residue_index[t_s][t_r][b];t_r2=CL->residue_index[t_s][t_r][b+1];t_w2=CL->residue_index[t_s][t_r][b+2];b+=3;
2293 p2=inv_pos[t_s2][t_r2];
2294 //score=((float)t_w/(float)NORM_F)*((float)t_w2/(float)NORM_F);
2295 score=MIN(((float)t_w/(float)NORM_F),((float)t_w2/(float)NORM_F));
2297 if (!used[p2][1] && score>0)
2299 used_list[nused++]=p2;
2309 //FILTER: Keep in the graph the edges where (p1->p2/(Sum (P1->x))>0.01
2312 for (a=0; a<nused; a++)
2316 nscore=used[p2][0]/tot; //Normalized score used for filtering
2319 used[p2][0]=used[p2][1]=0;
2322 max_n+=10000;list=vrealloc (list, max_n*sizeof (int*));
2324 if (nscore>filter && p1!=0 && p2!=0 && p1!=l1 && p2!=l2)
2327 list[n]=vcalloc (7, sizeof (int));
2331 list[n][3]=(l1-(p1))+(p2);
2332 score/=(float)((CL->S)->nseq*nr[0][p1]*nr[1][p2]);
2333 list[n][2]=(int)((float)score*(float)NORM_F);
2334 avg+=(int)((float)score*(float)NORM_F);
2340 free_float (used, -1);
2342 free_int (inv_pos, -1);
2354 int cl2pair_list_ecl_rawquad ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
2356 int p1, p2, si, s, r, t_s, t_r,t_w, t_s2, t_r2, t_w2,t_s3, t_r3, t_w3, n,n2;
2366 int *used_list, **used;
2367 int *sl2, **inv_pos;
2374 max_n=read_array_size (list, sizeof (int*));
2378 pos=aln2pos_simple ( A,-1, ns, ls);
2379 inv_pos=vcalloc ((CL->S)->nseq, sizeof (int*));
2380 for (a=0; a<ns[1]; a++)inv_pos[ls[1][a]] =seq2inv_pos(A->seq_al[ls[1][a]]);
2382 l1=strlen (A->seq_al[ls[0][0]]);
2383 l2=strlen (A->seq_al[ls[1][0]]);
2384 sl2=vcalloc ((CL->S)->nseq, sizeof (int));
2386 for (a=0;a<ns[1]; a++)sl2[ls[1][a]]=1;
2388 CL=index_res_constraint_list (CL, CL->weight_field);
2390 used=declare_int (l2+1,2);
2391 used_list=vcalloc (l2+1, sizeof (int));
2393 nseq2=(CL->S)->nseq*(CL->S)->nseq;
2395 for (p1=0; p1<=l1; p1++)
2397 for (nused=0,si=0;p1>0 && si<ns[0]; si++)
2399 s=ls [0][si];r=pos[s][p1-1];
2400 for (a=1; r>0 && a<CL->residue_index[s][r][0];a+=3)
2402 t_s=CL->residue_index[s][r][a];t_r=CL->residue_index[s][r][a+1];t_w=CL->residue_index[s][r][a+2];
2403 for (b=0; b<CL->residue_index[t_s][t_r][0];)
2405 if (b==0){t_s2=t_s;t_r2=t_r;t_w2=t_w;b++;}
2408 t_s2=CL->residue_index[t_s][t_r][b];t_r2=CL->residue_index[t_s][t_r][b+1];t_w2=CL->residue_index[t_s][t_r][b+2];b+=3;
2412 for (c=0; c<CL->residue_index[t_s2][t_r2][0];)
2414 if (c==0){t_s3=t_s2;t_r3=t_r2;t_w3=t_w2;c++;}
2417 t_s3=CL->residue_index[t_s2][t_r2][c];t_r3=CL->residue_index[t_s2][t_r2][c+1];t_w3=CL->residue_index[t_s2][t_r2][c+2];c+=3;
2422 p2=inv_pos[t_s3][t_r3];
2423 score=MIN(t_w,t_w2);
2424 score=MIN(score,t_w3);
2425 if (!used[p2][1] && score>0)
2427 used_list[nused++]=p2;
2437 for (a=0; a<nused; a++)
2443 used[p2][0]=used[p2][1]=0;
2445 if (n==max_n){max_n+=1000;list=vrealloc (list, max_n*sizeof (int*));}
2446 if (score>filter && p1!=0 && p2!=0 && p1!=l1 && p2!=l2)
2448 if (!list[n])list[n]=vcalloc (7, sizeof (int));
2452 list[n][3]=(l1-(p1))+(p2);
2462 free_int (inv_pos, -1);
2470 int cl2pair_list_ecl_norm ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
2472 int p1, p2, si, s, r, t_s, t_r, n,n2;
2482 int *used_list, *used;
2483 int *sl2, **inv_pos;
2490 max_n=read_array_size (list, sizeof (int*));
2494 pos=aln2pos_simple ( A,-1, ns, ls);
2495 inv_pos=vcalloc ((CL->S)->nseq, sizeof (int*));
2496 for (a=0; a<ns[1]; a++)inv_pos[ls[1][a]] =seq2inv_pos(A->seq_al[ls[1][a]]);
2498 l1=strlen (A->seq_al[ls[0][0]]);
2499 l2=strlen (A->seq_al[ls[1][0]]);
2500 sl2=vcalloc ((CL->S)->nseq, sizeof (int));
2502 for (a=0;a<ns[1]; a++)sl2[ls[1][a]]=1;
2504 CL=index_res_constraint_list (CL, CL->weight_field);
2506 used=vcalloc (l2+1, sizeof (int));
2507 used_list=vcalloc (l2+1, sizeof (int));
2512 for (p1=0; p1<=l1; p1++)
2514 for (si=0;p1>0 && si<ns[0]; si++)
2519 for (a=1; r>0 && a<CL->residue_index[s][r][0];a+=3)
2522 t_s=CL->residue_index[s][r][a];
2523 t_r=CL->residue_index[s][r][a+1];
2525 for (b=0; b<CL->residue_index[t_s][t_r][0];)
2528 if (b==0){t_s2=t_s;t_r2=t_r;b++;}
2529 else { t_s2=CL->residue_index[t_s][t_r][b];t_r2=CL->residue_index[t_s][t_r][b+1];b+=3;}
2533 p2=inv_pos[t_s2][t_r2];
2534 if (!used[p2]){used[p2]=1;used_list[nused++]=p2;}
2541 if (p1==0 || p1==l1)
2543 for (nused=0,p2=0; p2<=l2; p2++)used_list[nused++]=p2;
2547 if (!used[0])used_list[nused++]=0;
2548 if (!used[l2])used_list[nused++]=l2;
2550 for (a=0; a<nused; a++)
2553 if (p2==0 || p1==0)score=0;
2555 else score=(CL->get_dp_cost) (A, pos, ns[0], ls[0], p1-1, pos, ns[1], ls[1],p2-1,CL);
2556 if (score>filter && p1!=0 && p2!=0 && p1!=l1 && p2!=l2)
2558 if (n==max_n){max_n+=1000;list=vrealloc (list, max_n*sizeof (int*));}
2559 if (!list[n])list[n]=vcalloc (7, sizeof (int));
2562 list[n][3]=(l1-(p1))+(p2);
2565 if (p1!=0 && p2!=0 && p1!=l1 && p2!=l2)n2++;
2575 free_int (inv_pos, -1);
2588 int cl2pair_list_ref( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list, int *n)
2590 int a, b, l1, l2, n2=0;
2595 pos=aln2pos_simple ( A,-1, ns, ls);
2596 l1=strlen (A->seq_al[ls[0][0]]);
2597 l2=strlen (A->seq_al[ls[1][0]]);
2599 max_n=read_array_size (list[0], sizeof (int));
2601 for (a=0; a<=l1; a++)
2602 for (b=0; b<=l2; b++)
2604 score=(a==0 || b==0)?0:slow_get_dp_cost_pc(A, pos, ns[0], ls[0],a-1, pos, ns[1], ls[1], b-1, CL);
2607 if ( score>0 && a!=0 && b!=0 && a!=l1 && b!=l2)
2609 if (n[0]==max_n){max_n+=1000;list[0]=vrealloc (list[0], max_n*sizeof (int*));}
2610 if (!list[0][n[0]])list[0][n[0]]=vcalloc (7, sizeof (int));
2613 list[0][n[0]][3]=(l1-a)+b;
2614 list[0][n[0]][2]=score;
2615 if ( a!=0 && b!=0 && a!=l1 && b!=l2)
2627 int list2linked_pair_wise ( Alignment *A, int *ns, int **l_s, Constraint_list *CL, int **list, int n);
2628 int two_pass_linked_pair_wise ( Alignment *A, int *ns, int **l_s, Constraint_list *CL)
2630 int n=0, **list=NULL;
2635 cl2pair_list (A,ns, l_s, CL, &list, &n, mode, 10);
2636 nscore=list2linked_pair_wise (A, ns, l_s, CL, list, n);
2637 id=sub_aln2sim (A, ns, l_s, "idmat_sim");
2639 if (id>50)return nscore;
2640 ungap_sub_aln ( A, ns[0], l_s[0]);
2641 ungap_sub_aln ( A, ns[1], l_s[1]);
2642 cl2pair_list (A,ns, l_s, CL, &list, &n,mode,0);
2643 nscore=list2linked_pair_wise (A, ns, l_s, CL, list, n);
2644 cl2pair_list (NULL,ns, l_s, CL, &list, &n, mode, 0);
2647 int clinked_pair_wise ( Alignment *A, int *ns, int **l_s, Constraint_list *CL)
2649 int n=0, **list=NULL;
2650 int nscore, pscore=0;
2654 cl2pair_list (A,ns, l_s, CL, &list, &n, mode, 1000);
2655 nscore=list2linked_pair_wise (A, ns, l_s, CL, list, n);
2656 HERE ("***********First: %d", nscore);
2659 while (nscore>pscore)
2662 ungap_sub_aln ( A, ns[0], l_s[0]);
2663 ungap_sub_aln ( A, ns[1], l_s[1]);
2664 cl2pair_list (A,ns, l_s, CL, &list, &n, mode, 10);
2665 nscore=list2linked_pair_wise (A, ns, l_s, CL, list, n);
2666 HERE ("****************New: %d", nscore);
2669 cl2pair_list (NULL,ns, l_s, CL, &list, &n, mode, 0);
2672 int linked_pair_wise ( Alignment *A, int *nsi, int **lsi, Constraint_list *CL)
2675 static int **list=NULL;
2678 int mode=1;//1:ecl, 0:ref
2680 ns=vcalloc (2, sizeof (int));
2681 ns[0]=nsi[1]; ns[1]=nsi[0];
2683 ls=declare_int (2, ns[0]+ns[1]);
2684 for (a=0; a<ns[1]; a++)
2686 for (a=0; a<ns[0]; a++)
2690 cl2pair_list (A,ns, ls, CL, &list, &n, mode, 0);
2692 score=list2linked_pair_wise (A, ns, ls, CL, list, n);
2693 cl2pair_list (NULL,ns, ls, CL, &list, &n, mode, 0);
2698 int list2linked_pair_wise_nolgp ( Alignment *A, int *ns, int **l_s, Constraint_list *CL, int **list, int n);
2699 int list2linked_pair_wise ( Alignment *A, int *ns, int **l_s, Constraint_list *CL, int **list, int n)
2702 if (mode==0)return list2linked_pair_wise_nolgp (A, ns, l_s, CL, list, n);
2707 #define LIN(a,b,c) a[b*5+c]
2708 int list2linked_pair_wise_nolgp( Alignment *A, int *ns, int **l_s, Constraint_list *CL, int **list, int n)
2710 int a,b,c, i, j, LEN=0, start_trace;
2711 int pi, pj,ij, delta_i, delta_j, prev_i, prev_j;
2713 static long *MI, *MJ, *MM,*MT2;
2714 static int *sortseq;
2715 static int max_size;
2716 int gop, gep, igop, igep;
2719 char **aln,*char_buf;
2724 l1=strlen (A->seq_al[l_s[0][0]]);
2725 l2=strlen (A->seq_al[l_s[1][0]]);
2726 al=declare_char (2,l1+l2+1);
2729 //Penalties: max score is NORM_F
2730 //Penalties must be negative
2738 vfree (MI);vfree (MJ); vfree (MM);
2739 free_int (slist, -1);
2741 slist=declare_int (n,3);
2743 MI=vcalloc (5*n, sizeof (long));
2744 MJ=vcalloc (5*n, sizeof (long));
2745 MM=vcalloc (5*n, sizeof (long));
2751 for (b=0; b<5; b++)LIN(MI,a,b)=LIN(MJ,a,b)=LIN(MJ,a,b)=-1000000;
2755 if (!sortseq) sortseq=vcalloc( 7, sizeof (int));
2756 sortseq[0]=0; sortseq[1]=1;sortseq[2]=-1;
2757 sort_list_int2 (list, sortseq,7, 0, n-1);
2764 sortseq[0]=1; sortseq[1]=0;sortseq[2]=-1;
2765 sort_list_int2 (list, sortseq,7, 0, n-1);
2768 slist[a][1]=list[a][4];
2772 sortseq[0]=3; sortseq[1]=0;sortseq[2]=1;sortseq[3]=-1;
2773 sort_list_int2 (list, sortseq,7, 0, n-1);
2776 slist[a][2]=list[a][4];
2780 sortseq[0]=0; sortseq[1]=1;sortseq[2]=-1;
2781 sort_list_int2 (list, sortseq,7, 0, n-1);
2798 if (i==l1 || j==l2)gop=0;
2801 if (i==l1 && j==l2)start_trace=a;
2802 else if ( i==0 || j==0)
2804 LIN(MM,a,0)=-1000000;
2808 LIN(MJ,a,0)=-10000000;
2815 LIN(MI,a,0)=-10000000;
2820 LIN(MI,a,1)=LIN(MJ,a,1)=-1;
2821 LIN(MI,a,2)=LIN(MJ,a,2)=i;
2822 LIN(MI,a,3)=LIN(MJ,a,3)=j;
2846 delta_i=list[a][0]-list[pi][0];
2847 delta_j=list[a][1]-list[pj][1];
2850 LIN(MI,a,0)=MAX(LIN(MI,pi,0),(LIN(MM,pi,0)+gop))+delta_i*gep;
2852 LIN(MI,a,2)=delta_i;
2854 LIN(MI,a,4)=(LIN(MI,pi,0)>=(LIN(MM,pi,0)+gop))?'i':'m';
2857 LIN(MJ,a,0)=MAX(LIN(MJ,pj,0),(LIN(MM,pj,0)+gop))+delta_j*gep;
2860 LIN(MJ,a,3)=delta_j;
2862 LIN(MJ,a,4)=(LIN(MJ,pj,0)>=LIN(MM,pj,0)+gop)?'j':'m';
2865 if (a>1 && (ls=list[a][0]-list[ij][0])==(list[a][1]-list[ij][1]))
2867 LIN(MM,a,0)=MAX3(LIN(MM,ij,0),LIN(MI,ij,0),LIN(MJ,ij,0))+list[a][2]-(ls*CL->nomatch);
2872 if ( LIN(MM,ij,0)>=LIN(MI,ij,0) && LIN(MM,ij,0)>=LIN(MJ,ij,0))LIN(MM,a,4)='m';
2873 else if ( LIN(MI,ij,0) >= LIN(MJ,ij,0))LIN(MM,a,4)='i';
2874 else LIN(MM,a,4)='j';
2879 LIN(MM,a,0)=UNDEFINED;
2885 if (LIN(MM,a,0)>=LIN(MI,a,0) && LIN(MM,a,0) >=LIN(MJ,a,0))MT2=MM;
2886 else if ( LIN(MI,a,0)>=LIN(MJ,a,0))MT2=MI;
2889 score=MAX3(LIN(MM,a,0), LIN(MI,a,0), LIN(MJ,a,0));
2895 while (!(i==0 &&j==0))
2898 l=MAX(LIN(MT2,a,2),LIN(MT2,a,3));
2899 // HERE ("%c from %c %d %d SCORE=%d [%d %d] [%2d %2d]", T2[a][5],T2[a][4], T2[a][2], T2[a][3], T2[a][0], gop, gep, i, j);
2919 else if (l==0) {HERE ("L=0 i=%d j=%d",l, i, j);exit (0);}
2922 for (b=0; b<l; b++, LEN++)
2924 if (LIN(MT2,a,2)){al[0][LEN]=1;i--;ni++;}
2927 if (LIN(MT2,a,3)){al[1][LEN]=1;j--;nj++;}
2931 next_a=LIN(MT2,a,1);
2932 if (LIN(MT2,a,4)=='m')MT2=MM;
2933 else if (LIN(MT2,a,4)=='i')MT2=MI;
2934 else if (LIN(MT2,a,4)=='j')MT2=MJ;
2941 invert_list_char ( al[0], LEN);
2942 invert_list_char ( al[1], LEN);
2945 if ( A->declared_len<=LEN)A=realloc_aln ( A,2*LEN+1);
2947 char_buf= vcalloc (LEN+1, sizeof (char));
2949 for ( c=0; c< 2; c++)
2951 for ( a=0; a< ns[c]; a++)
2954 for ( b=0; b< LEN; b++)
2957 char_buf[b]=aln[l_s[c][a]][ch++];
2962 sprintf (aln[l_s[c][a]],"%s", char_buf);
2967 A->nseq=ns[0]+ns[1];
2974 int ** aln2local_penalties4link (Alignment *A, int n, int *ls, Constraint_list *CL, int **lg);
2975 int ** aln2local_penalties4link (Alignment *A, int n, int *ls, Constraint_list *CL, int **lg)
2977 //adapted from gap_count in MAFFT V 5.5
2980 int open=3, close=4, gap=5;
2985 l=strlen (A->seq_al[ls[0]]);
2989 lg=declare_int (6, l);
2992 if ( read_array_size_new (lg[0])<l)
2995 lg=declare_int (6, l);
2998 for( s=0; s<n; s++ )
3003 c2=A->seq_al[ls[s]][p];
3005 if (c1!='-' && c2=='-')lg[open][p]++;
3006 if (c1=='-' && c2!='-')lg[close][p]++;
3007 if ( c1=='-')lg[gap][p]++;
3020 lg[GOP][p]=0.5*(1-(go/nn))*gop;
3021 lg[GCP][p]=0.5*(1-(gc/nn))*gop;
3022 //Checked locacal gep => gives low quality results
3023 lg[GEP][p]=gep;//(1-((float)lg[gap][p]/(float)n))*gep;
3024 lg[open][p]=lg[close][p]=lg[gap][p]=0;
3031 /*********************************COPYRIGHT NOTICE**********************************/
3032 /*© Centro de Regulacio Genomica */
3034 /*Cedric Notredame */
3035 /*Tue Oct 27 10:12:26 WEST 2009. */
3036 /*All rights reserved.*/
3037 /*This file is part of T-COFFEE.*/
3039 /* T-COFFEE is free software; you can redistribute it and/or modify*/
3040 /* it under the terms of the GNU General Public License as published by*/
3041 /* the Free Software Foundation; either version 2 of the License, or*/
3042 /* (at your option) any later version.*/
3044 /* T-COFFEE is distributed in the hope that it will be useful,*/
3045 /* but WITHOUT ANY WARRANTY; without even the implied warranty of*/
3046 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the*/
3047 /* GNU General Public License for more details.*/
3049 /* You should have received a copy of the GNU General Public License*/
3050 /* along with Foobar; if not, write to the Free Software*/
3051 /* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA*/
3052 /*............................................... |*/
3053 /* If you need some more information*/
3054 /* cedric.notredame@europe.com*/
3055 /*............................................... |*/
3059 /*********************************COPYRIGHT NOTICE**********************************/