+++ /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 addE(i,j,d,s,list,n)\
- if (list)\
- {\
- int max_n; \
- Memcontrol *ppp; \
- if (!list[0])\
- {\
- list[0]=vcalloc ( 1000, sizeof (int*));\
- }\
- \
- ppp=(Memcontrol*)list[0]; \
- ppp-=2; \
- max_n=ppp[0].size/sizeof(int*);\
- if (n[0]>=max_n){max_n+=1000;list[0]=vrealloc (list[0], max_n*sizeof (int*));} \
- if (!list[0][n[0]])list[0][n[0]]=vcalloc (7, sizeof (int)); \
- list[0][n[0]][0]=i; \
- list[0][n[0]][1]=j; \
- list[0][n[0]][3]=d; \
- list[0][n[0]][2]=s; \
- n[0]++; \
- } \
-
-int cl2pair_list_ecl ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
-
-
-
-/*******************************************************************************/
-/* idscore_pairseq: measure the % id without delivering thze aln*/
-/* */
-/* makes DP between the the ns[0] sequences and the ns[1] */
-/* */
-/* for MODE, see the function get_dp_cost */
-/*******************************************************************************/
-int idscore_pairseq (char *s1, char *s2, int gop, int gep, int **m, char *comp_mode)
-{
- int **I, **D, **M, *P;
- int i, j,l1, l2, l,score, id, igop,match;
-
-
- l1=strlen (s1); l2=strlen (s2);
- lower_string (s1); lower_string (s2);
-
- I=declare_int (6,l2+1);D=declare_int (6,l2+1);M=declare_int (6,l2+1);
- for (j=0; j<=l2; j++)
- {
- D[0][j]=gep*j;M[0][j]=2*gep*j;D[4][j]=0;
- }
-
- for (i=1; i<=l1; i++)
- {
-
- I[1][0]=i*gep;
- M[1][0]=2*i*gep;
-
- for (j=1; j<=l2; j++)
- {
- score=m[s1[i-1]-'a'][s2[j-1]-'a'];
- id=(s1[i-1]==s2[j-1])?1:0;
-
- igop=(i==l1 || j==l2)?0:gop;
-
- 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];}
- else {D[1][j]=M[0][j]+igop+gep; D[3][j]=M[2][j]; D[5][j]=M[4][j];}
-
- 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];}
- else {I[1][j]=M[1][j-1]+igop+gep; I[3][j]=M[3][j-1]; I[5][j]=M[5][j-1];}
-
- match=M[0][j-1]+score;
- 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];}
- else if (D[1][j]>match) {M[1][j]=D[1][j] ; M[3][j]=D[3][j]; M[5][j]=D[5][j];}
- else {M[1][j]=match ; M[3][j]=M[2][j-1]+id; M[5][j]=M[4][j-1]+1;}
- }
- P=I[0]; I[0]=I[1]; I[1]=P;
- P=I[2]; I[2]=I[3]; I[3]=P;
- P=I[4]; I[4]=I[5]; I[5]=P;
-
- P=D[0]; D[0]=D[1]; D[1]=P;
- P=D[2]; D[2]=D[3]; D[3]=P;
- P=D[4]; D[4]=D[5]; D[5]=P;
-
- P=M[0]; M[0]=M[1]; M[1]=P;
- P=M[2]; M[2]=M[3]; M[3]=P;
- P=M[4]; M[4]=M[5]; M[5]=P;
- }
-
-
-
-
- if ( strstr (comp_mode, "sim2"))
- {
- l=MIN(l1,l2);
- score=(l==0)?0:(M[2][l2]*100)/l;
- }
- else if ( strstr (comp_mode, "sim3"))
- {
- l=MAX(l1,l2);
- score=(l==0)?0:(M[2][l2]*100)/l;
- }
- else if ( strstr (comp_mode, "cov"))
- {
- l=MAX(l1,l2);
- score=(l==0)?0:((M[4][l2]*100)/l);
- }
- else
- {
- //default: simple sim
- l=M[4][l2];
- score=(l==0)?0:(M[2][l2]*100)/l;
- }
-
- free_int (I, -1);
- free_int (D, -1);
- free_int (M, -1);
-
- return score;
-}
-
-int test_pair_wise (Alignment *A, int *ns, int **l_s, Constraint_list *CL)
-{
- int a,l0, l1, n;
- char buf[VERY_LONG_STRING];
- char *gap, *seq;
-
- l0=strlen (A->seq_al[l_s[0][0]]);
- l1=strlen (A->seq_al[l_s[1][0]]);
-
- n=(l0<5)?l0/2:5;
- gap=generate_null(l1-n);
- for (a=0;a<ns[0]; a++)
- {
- seq=A->seq_al[l_s[0][a]];
- sprintf (buf, "%s%s",seq, gap);
- sprintf (seq, "%s", buf);
- }
- vfree (gap);
- gap=generate_null(l0-n);
-
- for (a=0;a<ns[1]; a++)
- {
- seq=A->seq_al[l_s[1][a]];
- sprintf (buf, "%s%s",seq, gap);
- sprintf (seq, "%s", buf);
- }
- vfree(gap);
-
-
- A->len_aln=strlen (A->seq_al[l_s[0][0]]);
- A->score=A->score_aln=100;
- return 100;
-}
-
-int idscore_pair_wise (Alignment *A,int*ns, int **l_s,Constraint_list *CL)
-{
-
- 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");
- return A->score_aln;
-}
-int dp_max (int *trace, int n, ...);
-int dp_max (int *trace, int n, ...)
-{
- va_list ap;
- int a, v, t, best_v=0;
-
- va_start (ap, n);
- for (a=0; a< n; a++)
- {
- t=va_arg (ap, int);
- v=va_arg (ap, int);
-
- if (a==0)
- {
- best_v=v;
- trace[0]=t;
- }
- else
- {
- if (v>best_v)
- {
- best_v=v;
- trace[0]=t;
- }
- }
- }
-
- return best_v;
-}
-int is_tied (int *trace, int n, ...);
-int is_tied(int *trace, int n, ...)
-{
- va_list ap;
- int a, v, t, best_v=0;
- int nties=0;
-
- va_start (ap, n);
- for (a=0; a< n; a++)
- {
- t=va_arg (ap, int);
- v=va_arg (ap, int);
-
- if (a==0)
- {
- best_v=v;
- trace[0]=t;
- }
- else
- {
- if (v>best_v)
- {
- best_v=v;
- trace[0]=t;
- }
- }
- }
- va_end(ap);
- va_start (ap,n);
- for (a=0; a<n; a++)
- {
- t=va_arg (ap, int);
- v=va_arg (ap, int);
- if (v==best_v && trace[0]!=t)
- nties++;
- }
- va_end (ap);
- return nties;
-}
-
-void display_mat (int **M, int l1, int l2, char *title);
-void display_mat (int **M, int l1, int l2, char *title)
-{
- int a, b;
-
- fprintf ( stdout, "\n\nTitle %s\n", title);
- for ( a=0; a<=l1; a++)
- {
- fprintf ( stdout, "\n");
- for ( b=0; b<=l2; b++)
- fprintf ( stdout, "%3d ", M[a][b]);
- }
-}
-int glocal_pair_wise ( Alignment *A, int *ns, int **l_s, Constraint_list *CL)
-{
- int ***t, ***m;
- int i,j, l1, l2, n, sub, trace,ntrace, a, b, c, score;
- int gop,rgop,tgop, gep, unmatch;
- int M1, M2, I1, D1, LEN;
- char **al, *char_buf, **aln;
- int **pos0;
-
-
- l1=strlen (A->seq_al[l_s[0][0]]);
- l2=strlen (A->seq_al[l_s[1][0]]);
-
- n=1;
- M1=n++;D1=n++;I1=n++;M2=n++;
- t=declare_arrayN(3, sizeof (int),n, l1+1, l2+1);
- m=declare_arrayN(3, sizeof (int),n, l1+1, l2+1);
-
-
- gop=CL->gop*SCORE_K;
- gep=CL->gep*SCORE_K;
- tgop=gop;
- unmatch=gep;
-
- pos0=aln2pos_simple ( A,-1, ns, l_s);
-
-
- for (j=1; j<=l2; j++)
- {
- m[D1][0][j]=gep*j;
-
- m[M1][0][j]=2*gep*j;
- m[M2][0][j]=4*gep*j;
- }
-
-
- for (i=1; i<=l1; i++)
- {
- m[I1][i][0]=i*gep;
- m[M2][i][0]=4*i*gep;
- m[M1][i][0]=2*i*gep;
-
- for ( j=1; j<=l2; j++)
- {
- rgop=(i==l1 || j==1)?0:gop;
- rgop=gop;
- sub=(CL->get_dp_cost) (A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],j-1,CL);
- 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;
- t[M1][i][j]=trace;
-
- 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]);
- t[D1][i][j]=trace;
-
- 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]);
- t[I1][i][j]=trace;
-
- 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;
- t[M2][i][j]=trace;
-
- }
-
- }
- 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]);
- LEN=0;i=l1;j=l2;
- al=declare_char (2, l1+l2+1);
-
-
- trace=t[trace][i][j];
- while (!(i==0 &&j==0))
- {
-
- ntrace=t[trace][i][j];
- if (i==0)
- {
- al[0][LEN]=0;
- al[1][LEN]=1;
- j--;
- LEN++;
- }
- else if ( j==0)
- {
- al[0][LEN]=1;
- al[1][LEN]=0;
- i--;
- LEN++;
- }
- else if ( trace==M1)
- {
- al[0][LEN]=1;
- al[1][LEN]=1;
- i--; j--;
- LEN++;
- }
- else if ( trace==M2)
- {
- al[0][LEN]=1;
- al[1][LEN]=0;
- LEN++;
-
- al[0][LEN]=0;
- al[1][LEN]=1;
- LEN++;
-
- i--; j--;
-
- }
- else if ( trace==D1)
- {
- al[0][LEN]=0;
- al[1][LEN]=1;
- j--;
- LEN++;
- }
- else if ( trace == I1)
- {
- al[0][LEN]=1;
- al[1][LEN]=0;
- i--;
- LEN++;
- }
- trace=ntrace;
-
- }
-
- invert_list_char ( al[0], LEN);
- invert_list_char ( al[1], 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));
- for ( c=0; c< 2; c++)
- {
- for ( a=0; a< ns[c]; a++)
- {
- int ch=0;
- for ( b=0; b< LEN; b++)
- {
- if (al[c][b]==1)
- char_buf[b]=aln[l_s[c][a]][ch++];
- else
- char_buf[b]='-';
- }
- char_buf[b]='\0';
- sprintf (aln[l_s[c][a]],"%s", char_buf);
- }
- }
-
-
- A->len_aln=LEN;
- A->nseq=ns[0]+ns[1];
- free_arrayN((void *)m, 3);
- free_arrayN((void *)t, 3);
- vfree (char_buf);
- free_char (al, -1);
- return score;
-}
-
-int ** aln2local_penalties (Alignment *A, int n, int *ls, Constraint_list *CL, int **lg);
-int ** aln2local_penalties (Alignment *A, int n, int *ls, Constraint_list *CL, int **lg)
-{
- //adapted from gap_count in MAFFT V 5.5
- int p,s,l, c1, c2;
- int gep,gop;
- int open=3, close=4, gap=5;
-
- gop=CL->gop*SCORE_K;
- gep=CL->gep*SCORE_K;
-
- l=strlen (A->seq_al[ls[0]]);
-
- if (!lg)
- {
- lg=declare_int (6, l);
- }
-
- if ( read_array_size_new (lg[0])<l)
- {
- free_int (lg, -1);
- lg=declare_int (6, l);
- }
-
- for( s=0; s<n; s++ )
- {
- c1='x';
- for (p=0; p<l; p++)
- {
- c2=A->seq_al[ls[s]][p];
-
- if (c1!='-' && c2=='-')lg[open][p]++;
- if (c1=='-' && c2!='-')lg[close][p]++;
- if ( c1=='-')lg[gap][p]++;
- c1=c2;
- }
- }
-
- for (p=0; p<l; p++)
- {
- float go, gc, nn;
- nn=n;
- go=lg[open ][p];
- gc=lg[close][p];
-
-
- lg[GOP][p]=0.5*(1-(go/nn))*gop;
- lg[GCP][p]=0.5*(1-(gc/nn))*gop;
- //Checked locacal gep => gives low quality results
- lg[GEP][p]=gep;//(1-((float)lg[gap][p]/(float)n))*gep;
- lg[open][p]=lg[close][p]=lg[gap][p]=0;
-
- }
-
- return lg;
-}
-int free_gotoh_pair_wise_lgp()
-{
- return gotoh_pair_wise_lgp (NULL, NULL, NULL, NULL);
-}
-int gotoh_pair_wise_lgp ( Alignment *A, int *ns, int **l_s, Constraint_list *CL)
-{
- int i,j, li, lj, n, sub, trace,ntrace, a, b, c, score;
- int I, J;
- int M1, I1, D1, LEN;
- char **al, *char_buf, **aln;
- int **pos0, **pos;
- Alignment *Aln;
-
- int gop[2], gcp[2], gep[2];
- static int ***gpl, ***t, ***m;
- static int max_li, max_lj;
-
-
-
- //gotoh_pair_wise ( A, ns, l_s,CL);
- //ungap_sub_aln (A, ns[0], l_s[0]);
- //ungap_sub_aln (A, ns[1], l_s[1]);
-
- if (!A)
- {
- free_arrayN((void**)gpl, 3);
- free_arrayN((void**)t, 3);
- free_arrayN((void**)m, 3);
- max_li=max_lj=0;
- return 0;
- }
-
- I=0;J=1;
-
-
- li=strlen (A->seq_al[l_s[I][0]]);
- lj=strlen (A->seq_al[l_s[J][0]]);
-
- if ( !gpl)gpl=vcalloc ( 2, sizeof (int**));
- gpl[I]=aln2local_penalties (A,ns[I], l_s[I], CL,gpl[I]);
- gpl[J]=aln2local_penalties (A,ns[J], l_s[J], CL,gpl[J]);
-
-
- n=1;
- M1=n++;D1=n++;I1=n++;
-
- if ( li>max_li ||lj>max_lj )
- {
- free_arrayN((void**)t, 3);
- free_arrayN((void**)m, 3);
-
-
- max_li=li;
- max_lj=lj;
- t=declare_arrayN(3, sizeof (int),n, max_li+1, max_lj+1);
- m=declare_arrayN(3, sizeof (int),n, max_li+1, max_lj+1);
-
- }
- pos0=aln2pos_simple ( A,-1, ns, l_s);
-
- //Compatibility with Macro
- Aln=A;
- pos=pos0;
-
- for (j=1; j<=lj; j++)
- {
- gep[J]=gpl[J][GEP][j-1];
- m[D1][0][j]=gep[J]*j;
- m[I1][0][j]=m[D1][0][j]-1;
- m[M1][0][j]=m[D1][0][j]-1;
- }
-
- //D1: gap in sequence I
- //I1: gap in sequence J
-
-
- for (i=1; i<=li; i++)
- {
- gep[I]=gpl[I][GEP][i-1];
- gop[I]=gpl[I][GOP][i-1];
- gcp[I]=gpl[I][GCP][i-1];
-
- m[I1][i][0]=i*gep[I];
- m[D1][i][0]= m[I1][i][0]-1;
- m[M1][i][0]= m[I1][i][0]-1;
-
-
-
- gop[I]=(i==1 || i==li )?0:gop[I];
- gcp[I]=(i==1 || i==li )?0:gcp[I];
-
-
- for ( j=1; j<=lj; j++)
- {
-
- gep[J]=gpl[J][GEP][j-1];
- gop[J]=gpl[J][GOP][j-1];
- gcp[J]=gpl[J][GCP][j-1];
-
- //gep[J]=gep[I]=(gep[J]+gep[I])/2;
- //gop[J]=gop[I]=(gop[J]+gop[I])/2;
- //gcp[J]=gcp[I]=(gcp[J]+gcp[I])/2;
-
-
- gop[J]=(j==1 || j==lj )?0:gop[J];
- gcp[J]=(j==1 || j==lj )?0:gcp[J];
-
-
- //sub=(CL->get_dp_cost) (A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],j-1,CL);
- sub=TC_SCORE((i-1), (j-1));
-
- 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;
- t[M1][i][j]=trace;
-
-
- 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]);
- t[D1][i][j]=trace;
-
-
- 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]);
- t[I1][i][j]=trace;
-
- }
-
- }
- score=dp_max (&trace,3, M1,m[M1][li][lj],D1,m[D1][li][lj],I1, m[I1][li][lj]);
-
- LEN=0;i=li;j=lj;
- al=declare_char (2, li+lj);
-
-
- trace=t[trace][i][j];
- while (!(i==0 &&j==0))
- {
-
- ntrace=t[trace][i][j];
-
-
- if (i==0)
- {
- al[0][LEN]=0;
- al[1][LEN]=1;
- j--;
- LEN++;
- }
- else if ( j==0)
- {
- al[0][LEN]=1;
- al[1][LEN]=0;
- i--;
- LEN++;
- }
- else if ( trace==M1)
- {
- al[0][LEN]=1;
- al[1][LEN]=1;
- i--; j--;
- LEN++;
- }
- else if ( trace==D1)
- {
- al[0][LEN]=0;
- al[1][LEN]=1;
- j--;
- LEN++;
- }
- else if ( trace == I1)
- {
- al[0][LEN]=1;
- al[1][LEN]=0;
- i--;
- LEN++;
- }
- trace=ntrace;
-
- }
-
- invert_list_char ( al[0], LEN);
- invert_list_char ( al[1], LEN);
- if ( A->declared_len<=LEN)A=realloc_aln ( A,2*LEN+1);
-
- aln=A->seq_al;
- char_buf= vcalloc (LEN+1, sizeof (char));
- for ( c=0; c< 2; c++)
- {
- for ( a=0; a< ns[c]; a++)
- {
- int ch=0;
- for ( b=0; b< LEN; b++)
- {
- if (al[c][b]==1)
- char_buf[b]=aln[l_s[c][a]][ch++];
- else
- char_buf[b]='-';
- }
- char_buf[b]='\0';
- sprintf (aln[l_s[c][a]],"%s", char_buf);
- }
- }
-
-
- A->len_aln=LEN;
- A->nseq=ns[0]+ns[1];
- vfree (char_buf);
- free_char (al, -1);
- free_int (pos0, -1);
- return score;
-}
-/*******************************************************************************/
-/* GLOCAL 2 */
-/* */
-/* makes DP between the the ns[0] sequences and the ns[1] */
-/* */
-/* for MODE, see the function get_dp_cost */
-/*******************************************************************************/
-int glocal2_pair_wise (Alignment *IN,int*ns, int **ls,Constraint_list *CL)
-{
- int a, b, s=0;
- Alignment *A, *R,*L;
- char *seq, *buf;
-
- buf=vcalloc (1000, sizeof (char));
- seq=vcalloc (1000, sizeof (char));
-
- A=copy_aln (IN,NULL);
- L=copy_aln (IN,NULL);
- R=copy_aln (IN,NULL);
-
- gotoh_pair_wise_sw (A, ns, ls, CL);
-
- HERE ("1");
- for (a=0; a<2; a++)
- {
- for (b=0; b<ns[a]; b++)
- {
- s=ls[a][b];
- sprintf ( seq,"%s", IN->seq_al[s]);
-
- seq[A->order[s][2]]='\0';
- sprintf (L->seq_al[s], "%s", seq);
- sprintf (R->seq_al[s], "%s", seq+A->order[s][3]+1);
- }
- }
- HERE ("2");
- print_sub_aln (A, ns, ls);
- gotoh_pair_wise(L, ns, ls, CL);
- print_sub_aln (L, ns, ls);
- gotoh_pair_wise(R, ns, ls, CL);
- print_sub_aln (R, ns, ls);
-
- IN=realloc_aln (IN, A->len_aln+L->len_aln+R->len_aln+1);
- for (a=0; a<2; a++)
- {
- for (b=0; b<ns[a]; b++)
- {
- s=ls[a][b];
- sprintf (IN->seq_al[s], "%s%s%s",L->seq_al[s], A->seq_al[s], R->seq_al[s]);
- }
- }
- IN->len_aln=strlen (IN->seq_al[s]);
-
- print_sub_aln (IN, ns, ls);
- vfree (seq); vfree (buf);
- free_aln (A); free_aln (L);free_aln (R);
- return IN->score_aln;
-}
-
-
-int gotoh_pair_wise (Alignment *A,int*ns, int **l_s,Constraint_list *CL)
- {
-/*******************************************************************************/
-/* NEEDLEMAN AND WUNSCH (GOTOH) */
-/* */
-/* makes DP between the the ns[0] sequences and the ns[1] */
-/* */
-/* for MODE, see the function get_dp_cost */
-/*******************************************************************************/
-
-
-/*TREATMENT OF THE TERMINAL GAP PENALTIES*/
-/*TG_MODE=0---> gop and gep*/
-/*TG_MODE=1---> --- gep*/
-/*TG_MODE=2---> --- ---*/
-
-
- int TG_MODE;
- int l_gop, l_gep;
- int gop, gep;
- int maximise;
-/*VARIANLES FOR THE MULTIPLE SEQUENCE ALIGNMENT*/
- int a, b, i, j;
-
- int *cc;
- int *dd,*ddg;
- int e, eg;
-
- int lenal[2], len;
- int t, c=0,s, ch;
- int sub;
- int fop;
- int score=0;
- int **pos0;
- static char **al;
- char **aln;
- int ala, alb,LEN;
- char *buffer;
- char *char_buf;
-/*trace back variables */
- static int **trace;
- static int **bit;
- int *bi;
- int *tr;
- int dim;
- int ibit=0;
- int k;
- int sample=0;//road==0, random tie; road=1: upper road; road=2 lower road;
- /********Prepare penalties*******/
- gop=CL->gop*SCORE_K;
- gep=CL->gep*SCORE_K;
- TG_MODE=CL->TG_MODE;
- maximise=CL->maximise;
-
-
-/********************************/
-/*CLEAN UP AFTER USE*/
- if ( A==NULL)
- {
- free_int (trace,-1);
- free_int (bit, -1);
- trace=NULL;
- bit=NULL;
-
- free_char (al,-1);
- al=NULL;
- return 0;
- }
-
-/*DO MEMORY ALLOCATION FOR DP*/
-
-
- sample=atoigetenv ("SAMPLE_DP_4_TCOFFEE");
-
- lenal[0]=strlen (A->seq_al[l_s[0][0]]);
- lenal[1]=strlen (A->seq_al[l_s[1][0]]);
- len= MAX(lenal[0],lenal[1])+1;
-
- buffer=vcalloc ( 2*len, sizeof (char));
- al=declare_char (2, 2*len);
-
- char_buf= vcalloc (2*len, sizeof (char));
-
-
- dd = vcalloc (len, sizeof (int));
-
-
- cc = vcalloc (len, sizeof (int));
- ddg=vcalloc (len, sizeof (int));
-
-
-
-
-
- dim=(trace==NULL)?0:read_size_int ( trace,sizeof (int*));
- trace =realloc_int ( trace,dim,dim,MAX(0,len-dim), MAX(0,len-dim));
- bit =realloc_int ( bit,dim,dim,MAX(0,len-dim), MAX(0,len-dim));
-
-/*END OF MEMORY ALLOCATION*/
-
-
- /*
- 0(s) +(dd)
- \ |
- \ |
- \ |
- \ |
- \ |
- \ |
- \|
- -(e)----O
- */
-
- pos0=aln2pos_simple ( A,-1, ns, l_s);
-
-
- cc[0]=0;
- tr=trace[0];
- bi=bit[0];
- tr[0]=1;
- for ( j=1; j<=lenal[1]; j++)tr[j]=-1;
-
- t=(TG_MODE==0)?gop:0;
-
-
- for (cc[0]=0,j=1; j<=lenal[1]; j++)
- {
-
- l_gop=(TG_MODE==0)?gop:0;
- l_gep=(TG_MODE==2)?0:gep;
-
- cc[j]=t=t+l_gep;
- dd[j]= t+ gop;
- }
-
- t=(TG_MODE==0)?gop:0;
-
- for (i=1; i<=lenal[0];i++)
- {
- tr=trace[i];
- bi=bit[i];
- s=cc[0];
-
- l_gop=(TG_MODE==0)?gop:0;
- l_gep=(TG_MODE==2)?0:gep;
-
-
-
- cc[0]=c=t=t+l_gep;
- e=t+ gop;
- tr[0]=1;
-
-
-
- for (eg=0,j=1; j<=lenal[1];j++)
- {
-
- sub=(CL->get_dp_cost) (A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],j-1,CL);
-
- /*get the best Insertion*/
- l_gop=(i==lenal[0] || i==1 )?((TG_MODE==0)?gop:0):gop;
- l_gep=(i==lenal[0] || i==1)?((TG_MODE==2)?0:gep):gep;
-
-
- if ( a_better_than_b ( e,c+l_gop, maximise))eg++;
- else eg=1;
- e=best_of_a_b (e, c+l_gop, maximise)+l_gep;
-
- /*Get the best deletion*/
- l_gop=(j==lenal[1] || j==1)?((TG_MODE==0)?gop:0):gop;
- l_gep=(j==lenal[1] || j==1)?((TG_MODE==2)?0:gep):gep;
-
-
- if ( a_better_than_b ( dd[j], cc[j]+l_gop, maximise))ddg[j]++;
- else ddg[j]=1;
- dd[j]=best_of_a_b( dd[j], cc[j]+l_gop,maximise)+l_gep;
-
-
-
- c=best_int(3,maximise,&fop, e, s+sub,dd[j]);
-
-
- if (sample==1)
- {
- int rr[3];
- int nn=0;
- int fop2;
- int ind;
- if (c==e)rr[nn++]=0;
- if (c==(s+sub))rr[nn++]=1;
- if (c==dd[j])rr[nn++]=2;
- ind=rand()%(nn);
- fop=rr[ind];
- if (nn>1)
- {
- // HERE ("NN=%d index=%d",nn, ind);
- //HERE ("%d ->%d", fop, fop2);
- //HERE ("%d %d %d", e, s+sub,dd[j]);
- ;
- }
- }
- else if (sample==0)
- {
- /*Chose Substitution for tie breaking*/
- if ( fop==0 && (s+sub)==e)fop=1;
- else if ( fop==2 && (s+sub)==dd[j])fop=1;
- /*Chose Deletion for tie breaking*/
- else if ( fop==2 && e==dd[j])fop=2;
- }
- else if (sample==-1)
- {
-
- if ( fop==0 && (s+sub)==e)fop=1;
- else if ( fop==1 && (s+sub)==dd[j])fop=2;
- /*Chose Deletion for tie breaking*/
- else if ( fop==2 && e==dd[j])fop=1;
- }
- bi[j]=0;
- if (c==e){bi[j]++;}
- if (c==(s+sub)){bi[j]++;}
- if (c==(dd[j])){bi[j]++;}
- //bi[j]--;
-
-
- fop-=1;
- s=cc[j];
- cc[j]=c;
-
-
-
-
-
- if ( fop<0)
- {tr[j]=(TRACE_TYPE)fop*eg;
- }
- else if ( fop>0)
- {tr[j]=(TRACE_TYPE)fop*ddg[j];
- }
- else if (fop==0)
- {tr[j]=(TRACE_TYPE)0;
- }
- fop= -2;
- }
-
- }
-
- score=c;
-
- i=lenal[0];
- j=lenal[1];
- ala=alb=0;
-
- if (!A->ibit)A->ibit=1; //set the bit counter on
- while (i>=0 && j>=0 && ((i+j)!=0))
- {
- if ( i==0)
- k=-1;
- else if ( j==0)
- k=1;
- else if ( j==0 && i==0)
- k=1;
- else
- {
- k=trace[i][j];
- A->ibit*=bit[i][j];
- }
-
-
- if (k==0)
- {
-
- al[0][ala++]=1;
- al[1][alb++]=1;
- i--;
- j--;
- }
- else if (k>0)
- {
-
- for ( a=0; a< k; a++)
- {
- al[0][ala++]=1;
- al[1][alb++]=0;
- i--;
- }
- }
- else if (k<0)
- {
-
- for ( a=0; a>k; a--)
- {
- al[0][ala++]=0;
- al[1][alb++]=1;
- j--;
- }
- }
- }
-
- LEN=ala;
- c=LEN-1;
-
-
-
- invert_list_char ( al[0], LEN);
- invert_list_char ( al[1], LEN);
- if ( A->declared_len<=LEN)A=realloc_aln2 ( A,A->max_n_seq, 2*LEN);
- aln=A->seq_al;
-
- for ( c=0; c< 2; c++)
- {
- for ( a=0; a< ns[c]; a++)
- {
- ch=0;
- for ( b=0; b< LEN; b++)
- {
- if (al[c][b]==1)
- char_buf[b]=aln[l_s[c][a]][ch++];
- else
- char_buf[b]='-';
- }
- char_buf[b]='\0';
- sprintf (aln[l_s[c][a]],"%s", char_buf);
- }
- }
-
-
- A->len_aln=LEN;
- A->nseq=ns[0]+ns[1];
-
-
- vfree ( cc);
- vfree (dd);
- vfree (ddg);
- vfree (buffer);
- vfree (char_buf);
-
- free_char ( al, -1);
- free_int (pos0, -1);
-
-
-
-
- return score;
- }
-
-
-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);
-int gotoh_pair_wise_lgp_sticky ( Alignment *A, int *ns, int **l_s, Constraint_list *CL)
-{
- int i,j, li, lj, n, sub, trace,ntrace, a, b, c, score;
- int I, J;
- int M1, I1, D1, LEN;
- char **al, *char_buf, **aln;
- int **pos0;
-
- int gop[2], gcp[2], gep[2];
- static int ***gpl, ***t, ***m;
- static int max_li, max_lj;
-
-
-
- //gotoh_pair_wise ( A, ns, l_s,CL);
- //ungap_sub_aln (A, ns[0], l_s[0]);
- //ungap_sub_aln (A, ns[1], l_s[1]);
-
- I=0;J=1;
-
-
- li=strlen (A->seq_al[l_s[I][0]]);
- lj=strlen (A->seq_al[l_s[J][0]]);
-
- if ( !gpl)gpl=vcalloc ( 2, sizeof (int**));
- gpl[I]=aln2local_penalties (A,ns[I], l_s[I], CL,gpl[I]);
- gpl[J]=aln2local_penalties (A,ns[J], l_s[J], CL,gpl[J]);
-
-
- n=1;
- M1=n++;D1=n++;I1=n++;
-
- if ( li>max_li ||lj>max_lj )
- {
- free_arrayN((void**)t, 3);
- free_arrayN((void**)m, 3);
-
-
- max_li=li;
- max_lj=lj;
- t=declare_arrayN(3, sizeof (int),n, max_li+1, max_lj+1);
- m=declare_arrayN(3, sizeof (int),n, max_li+1, max_lj+1);
-
- }
- pos0=aln2pos_simple ( A,-1, ns, l_s);
-
-
- for (j=1; j<=lj; j++)
- {
- gep[J]=gpl[J][GEP][j-1];
- m[D1][0][j]=gep[J]*j;
- m[I1][0][j]=m[D1][0][j]-1;
- m[M1][0][j]=m[D1][0][j]-1;
- }
-
- //D1: gap in sequence I
- //I1: gap in sequence J
-
-
- for (i=1; i<=li; i++)
- {
- gep[I]=gpl[I][GEP][i-1];
- gop[I]=gpl[I][GOP][i-1];
- gcp[I]=gpl[I][GCP][i-1];
-
- m[I1][i][0]=i*gep[I];
- m[D1][i][0]= m[I1][i][0]-1;
- m[M1][i][0]= m[I1][i][0]-1;
-
-
-
- gop[I]=(i==1 || i==li )?0:gop[I];
- gcp[I]=(i==1 || i==li )?0:gcp[I];
-
-
- for ( j=1; j<=lj; j++)
- {
- int transition;
-
- gep[J]=gpl[J][GEP][j-1];
- gop[J]=gpl[J][GOP][j-1];
- gcp[J]=gpl[J][GCP][j-1];
-
- //gep[J]=gep[I]=(gep[J]+gep[I])/2;
- //gop[J]=gop[I]=(gop[J]+gop[I])/2;
- //gcp[J]=gcp[I]=(gcp[J]+gcp[I])/2;
-
-
- gop[J]=(j==1 || j==lj )?0:gop[J];
- gcp[J]=(j==1 || j==lj )?0:gcp[J];
-
-
- sub=(CL->get_dp_cost) (A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],j-1,CL);
- transition=get_transition_cost (A, pos0, ns[0], l_s[0], i-1, pos0, ns[1], l_s[1],j-1,CL);
-
- 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;
- t[M1][i][j]=trace;
-
-
- 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]);
- t[D1][i][j]=trace;
-
-
- 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]);
- t[I1][i][j]=trace;
-
- }
-
- }
- score=dp_max (&trace,3, M1,m[M1][li][lj],D1,m[D1][li][lj],I1, m[I1][li][lj]);
-
- LEN=0;i=li;j=lj;
- al=declare_char (2, li+lj);
-
-
- trace=t[trace][i][j];
- while (!(i==0 &&j==0))
- {
-
- ntrace=t[trace][i][j];
-
-
- if (i==0)
- {
- al[0][LEN]=0;
- al[1][LEN]=1;
- j--;
- LEN++;
- }
- else if ( j==0)
- {
- al[0][LEN]=1;
- al[1][LEN]=0;
- i--;
- LEN++;
- }
- else if ( trace==M1)
- {
- al[0][LEN]=1;
- al[1][LEN]=1;
- i--; j--;
- LEN++;
- }
- else if ( trace==D1)
- {
- al[0][LEN]=0;
- al[1][LEN]=1;
- j--;
- LEN++;
- }
- else if ( trace == I1)
- {
- al[0][LEN]=1;
- al[1][LEN]=0;
- i--;
- LEN++;
- }
- trace=ntrace;
-
- }
-
- invert_list_char ( al[0], LEN);
- invert_list_char ( al[1], LEN);
- if ( A->declared_len<=LEN)A=realloc_aln ( A,2*LEN+1);
-
- aln=A->seq_al;
- char_buf= vcalloc (LEN+1, sizeof (char));
- for ( c=0; c< 2; c++)
- {
- for ( a=0; a< ns[c]; a++)
- {
- int ch=0;
- for ( b=0; b< LEN; b++)
- {
- if (al[c][b]==1)
- char_buf[b]=aln[l_s[c][a]][ch++];
- else
- char_buf[b]='-';
- }
- char_buf[b]='\0';
- sprintf (aln[l_s[c][a]],"%s", char_buf);
- }
- }
-
-
- A->len_aln=LEN;
- A->nseq=ns[0]+ns[1];
- vfree (char_buf);
- free_char (al, -1);
- free_int (pos0, -1);
- return score;
-}
-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)
-{
- /*counts the number of identical transitions between position i-1, i and j-1..j*/
- float t=0;
- int a,s;
- Sequence *S;
-
- if (i==0 || j==0)return 0;
-
- for (a=0; a<ni; a++)
- {
- s=li[a];
- if (posi[s][i]<0 || posi[s][i-1]<0)continue;
- if (S->seq[li[a]][i-1]==S->seq[li[a]][i-1])t++;
- }
-
- for (a=0; a<nj; a++)
- {
- s=lj[a];
- if (posj[s][j]<0 || posj[s][j-1]<0)continue;
- if (S->seq[li[a]][j-1]==S->seq[li[a]][j-1])t++;
- }
-
- t=(t*10)/(float)(ni+nj);
- return t;
-}
-/*******************************************************************************/
-/* idscore_pairseq: measure the % id without delivering thze aln*/
-/* */
-/* makes DP between the the ns[0] sequences and the ns[1] */
-/* */
-/* for MODE, see the function get_dp_cost */
-/*******************************************************************************/
-
-int cl2pair_list ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in, int mode, int ndiag);
-int cl2pair_list_ref ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
-int cl2pair_list_ecf ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
-int cl2pair_list_diag ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in, int add);
-int cl2list_borders (Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
-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
-int** cl2sorted_diagonals ( Alignment *A, int *ns, int **ls, Constraint_list *CL);
-int** cl2sorted_diagonals_mat ( Alignment *A, int *ns, int **ls, Constraint_list *CL);
-int** cl2sorted_diagonals_cs ( Alignment *A, int *ns, int **ls, Constraint_list *CL);
-int list2nodup_list (Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
-int fill_matrix ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
-
-int list2nodup_list ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
-{
- int **list;
- int n, a, b, c;
-
- list=list_in[0];
- n=n_in[0];
-
- if ( !A)return 0;
-
-
- sort_list_int (list,7, 1, 0, n-1);
- for (b=a=1; a<n; a++)
- {
- if (list[a][0]==list[b-1][0] && list[a][1]==list[b-1][1])
- {
- //HERE ("Duplicate");
- list[b-1][2]=MAX(list[b-1][2],list[a][2]);
- }
- else
- {
- for (c=0; c<4; c++)list[b][c]=list[a][c];
- b++;
- }
-
- }
- n_in[0]=b;
- list_in[0]=list;
- return b;
-}
-
-int cl2list_borders (Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
-{
- int a,n, p1, p2, l1, l2;
- int **list;
- int **pos;
- if (!A)return 0;
-
-
-
-
- l1=strlen (A->seq_al[ls[0][0]]);
- l2=strlen (A->seq_al[ls[1][0]]);
-
- for (p1=0; p1<=l1; p1++)
- {
- if (p1==0 || p1==l1)
- {
- for (p2=0; p2<=l2; p2++)
- {
- addE(p1,p2,((l1-(p1))+(p2)),((CL->gep)*SCORE_K*p2), list_in, n_in);
- }
- }
- else
- {
- for (a=0; a<2; a++)
- {
- p2=(a==0)?0:l2;
- addE(p1,p2,((l1-(p1))+(p2)),((CL->gep)*SCORE_K*p1), list_in, n_in);
- }
- }
- }
-
- return read_array_size (list_in[0], sizeof (int*));
-}
-
-int cl2diag_cap (Alignment *A, int *nns, int **ls, Constraint_list *CL, int ***list, int *n)
-{
- int *sortseq;
-
- int in, a, b, al1, al2;
- int max_n;
- int cap=0;
- int k=0;
-
- static int **ll;
- static int max_ll;
- int nll=0;
-
- int ns=0;
- int nt=0;
- int i,j,si,sj,ti,tj;
-
- if ( !A)vfree (ll);max_ll=0;
-
- al1=strlen (A->seq_al[ls[0][0]]);
- al2=strlen (A->seq_al[ls[1][0]]);
-
- sortseq=vcalloc (7, sizeof (int));
- sortseq[0]=3;sortseq[1]=0;sortseq[2]=-1;
- sort_list_int2 (list[0], sortseq,4, 0, n[0]-1);
- vfree(sortseq);
- in=n[0];
-
-
- if (!ll){max_ll=100;ll=vcalloc(max_ll,sizeof(int*));}
-
- for (a=0; a<in; a++)
- {
- int i, j, pi, pj, ni, nj;
- if (list[0][a][2]==0)continue;//this is where borders are avoided
- i=list[0][a][0];
- j=list[0][a][1];
-
- if (a==0){pi=-10;pj=-10;}
- else {pi=list[0][a-1][0];pj=list[0][a-1][1];}
-
- if (a==in-1){ni=-10; nj=-10;}
- else {ni=list[0][a+1][0]; nj=list[0][a+1][1];}
-
-
- if ((i==0 || j==0));
- else if ( i==pi || j==pj);
- else if ( i-pi!=1 || j-pj!=1)
- {
- if (nll>=max_ll){max_ll+=1000;ll=vrealloc (ll, max_ll*sizeof (int*));}
- ll[nll++]=list[0][a];
- list[0][a][6]=_START;
- }
-
- if (i==al1 || j==al2);
- else if ( i==ni || j==nj);
- else if ( ni-i!=1 || nj-j!=1)
- {
- if (nll>=max_ll){max_ll+=1000;ll=vrealloc (ll, max_ll*sizeof (int*));}
- ll[nll++]=list[0][a];
- list[0][a][6]=_TERM;
- }
- }
-
- sortseq=vcalloc (7, sizeof (int));
- sortseq[0]=0;sortseq[1]=1;sortseq[2]=-1;
- sort_list_int2 (ll, sortseq,4, 0,nll-1);
- vfree (sortseq);
-
- for (a=0; a<nll; a++)
- {
- int ci, nl,max_nl,best_d,d,best_s;
- max_nl=100;
- if (ll[a][6]!=_TERM)continue;
-
- ti=ll[a][0];
- tj=ll[a][1];
- ci=ti;
-
- for (nl=0,best_d=-1,b=a+1;b<nll && nl<max_nl; b++)
- {
- if (ll[b][6]!=_START)continue;
-
- si=ll[b][0];
- sj=ll[b][1];
-
- if (si>ci){nl++;ci=si;}
- d=MIN((si-ti), (sj-tj));
- if (d<=0);
- else if (best_d==-1 || best_d>d){best_d=d; best_s=b;}
- }
- if (best_d==-1)continue;
-
- si=ll[best_s][0];
- sj=ll[best_s][1];
-
- for (i=ti, j=tj; (i<=si && j<=sj); i++, j++)//extend the top diagonal
- {
- addE(i,j,(al1-i+j),cap, list,n);
- }
-
- for (i=si, j=sj; (i>=ti && j>=tj); i--, j--)//extend the bottom diagonal
- {
- addE(i,j,(al1-i+j),cap, list,n);
- }
- }
-
- for (a=0; a<nll; a++)ll[a][6]=0;
-
- return n[0];
-}
-
-int cl2pair_list_ecl_pc ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
-
-/**
- * Calculates scores for diagonal segments.
- *
- * \param Alignment The sequences.
- * \param ns Number of sequences in each group
- * \param ls sequences in in groups (ls[0][x] sequences in group 1, ls[1][x] squences in group 2).
- * \param CL the constraint list
- * \param list_in the diagonals
- * \param n_in number of sequences?
- */
-int fork_cl2pair_list_ecl_pc ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
-int nfork_cl2pair_list_ecl_pc ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
-int cl2pair_list_ecl_pc ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
-{
-
- if (!CL || !CL->S || !CL->residue_index) return 0;
-
-
- if ( get_nproc()==1)return nfork_cl2pair_list_ecl_pc(A,ns,ls,CL,list_in,n_in);
- else if (strstr ( CL->multi_thread, "pairwise"))return fork_cl2pair_list_ecl_pc(A,ns,ls,CL,list_in,n_in);
- else return nfork_cl2pair_list_ecl_pc(A,ns,ls,CL,list_in,n_in);
-}
-
-
-int fork_cl2pair_list_ecl_pc ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
-{
- int p1, p2,diag, si, s, r, t_s, t_r,t_w, t_s2, t_r2, t_w2;
- int a, b, l1, l2;
- int **pos;
-
- int nused;
- int *used_list;
- int *sl2,*sl1, **inv_pos;
-
-
-
- float nscore, score, tot, filter, avg=0, new=0;
- float **used;
- float *norm;
-
- //variables for fork
- FILE *fp;
- char **pid_tmpfile;
- int sjobs, njobs,j;
- int **sl;
-
-
- if ( !A) return 0;
-
-
-
- pos=aln2pos_simple ( A,-1, ns, ls);
- inv_pos=vcalloc ((CL->S)->nseq, sizeof (int*));
- for (a=0; a<ns[1]; a++)inv_pos[ls[1][a]] =seq2inv_pos(A->seq_al[ls[1][a]]);
-
- l1=strlen (A->seq_al[ls[0][0]]);
- l2=strlen (A->seq_al[ls[1][0]]);
- sl1=vcalloc ((CL->S)->nseq, sizeof (int));
- sl2=vcalloc ((CL->S)->nseq, sizeof (int));
-
- for (a=0;a<ns[0]; a++)sl1[ls[0][a]]=1;
- for (a=0;a<ns[1]; a++)sl2[ls[1][a]]=1;
- norm=vcalloc ( l1+1, sizeof (float));
-
- njobs=get_nproc();
- sl=n2splits (njobs,l1+1);
- pid_tmpfile=vcalloc (njobs, sizeof (char*));
-
- used=declare_float (l2+1,2);
- used_list=vcalloc (l2+1, sizeof (int));
- nused=0;
-
- for (sjobs=0, j=0; j<njobs; j++)
- {
- pid_tmpfile[j]=vtmpnam(NULL);
- if (vvfork (NULL)==0)
- {
- initiate_vtmpnam(NULL);
- fp=vfopen (pid_tmpfile[j], "w");
- for (p1=sl[j][0]; p1<sl[j][1]; p1++)
- {
- for (tot=0,nused=0,si=0;p1>0 && si<ns[0]; si++)
- {
- s=ls [0][si];r=pos[s][p1-1];
- for (a=1; r>0 && a<CL->residue_index[s][r][0];a+=ICHUNK)
- {
- t_s=CL->residue_index[s][r][a+SEQ2];
- t_r=CL->residue_index[s][r][a+R2];
- t_w=CL->residue_index[s][r][a+WE];
- if (sl1[t_s])continue;//do not extend within a profile
-
- norm[p1]++;
- for (b=0; b<CL->residue_index[t_s][t_r][0];)
- {
- if (b==0){t_s2=t_s;t_r2=t_r;t_w2=t_w;b++;}
- else
- {
- t_s2=CL->residue_index[t_s][t_r][b+SEQ2];
- t_r2=CL->residue_index[t_s][t_r][b+R2];
- t_w2=CL->residue_index[t_s][t_r][b+WE];
- b+=ICHUNK;
- }
- if (sl2[t_s2])
- {
- p2=inv_pos[t_s2][t_r2];
- score=MIN(((float)t_w/(float)NORM_F),((float)t_w2/(float)NORM_F));
-
- if (!used[p2][1] && score>0)
- {
- used_list[nused++]=p2;
- }
-
- tot+=score;
- used[p2][0]+=score;
- used[p2][1]++;
- }
- }
- }
- }
- filter=0.01;
- for (a=0; a<nused; a++)
- {
-
- p2=used_list[a];
- nscore=used[p2][0]/tot; //Normalized score used for filtering
- score =used[p2][0];
- used[p2][0]=used[p2][1]=0;
-
- if (nscore>filter && p1!=0 && p2!=0 && p1!=l1 && p2!=l2)
- {
- score=((norm[p1]>0)?score/norm[p1]:0)*NORM_F;
- fprintf (fp, "%d %d %d %f ", p1, p2, ((l1-(p1))+(p2)), score);
- }
- }
- }
- vfclose (fp);
- myexit (EXIT_SUCCESS);
- }
- else
- {
- sjobs++;
- }
- }
- while (sjobs>=0){vwait(NULL); sjobs--;}//wait for all jobs to complete
- for (j=0; j<njobs; j++)
- {
- fp=vfopen (pid_tmpfile[j], "r");
- while ((fscanf(fp, "%d %d %d %f ", &p1,&p2, &diag, &score))==4)
- addE (p1,p2,((l1-(p1))+(p2)),score,list_in, n_in);
- vfclose (fp);
- remove (pid_tmpfile[j]);
- }
-
- free_float (used, -1);
- vfree (used_list);
- free_int (inv_pos, -1);
- free_int (pos, -1);
- vfree (sl2);vfree (sl1);
- vfree(norm);
- return n_in[0];
-}
-
-
-
-
-int nfork_cl2pair_list_ecl_pc ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
-{
- int p1, p2, si, s, r, t_s, t_r,t_w, t_s2, t_r2, t_w2;
- int a, b, l1, l2;
- int **pos;
-
- int nused;
- int *used_list;
- int *sl2,*sl1, **inv_pos;
-
-
- float nscore, score, tot, filter, avg=0, new=0;
- float **used;
- float *norm;
-
- if ( !A) return 0;
-
- pos=aln2pos_simple ( A,-1, ns, ls);
- inv_pos=vcalloc ((CL->S)->nseq, sizeof (int*));
- for (a=0; a<ns[1]; a++)inv_pos[ls[1][a]] =seq2inv_pos(A->seq_al[ls[1][a]]);
-
- l1=strlen (A->seq_al[ls[0][0]]);
- l2=strlen (A->seq_al[ls[1][0]]);
- sl1=vcalloc ((CL->S)->nseq, sizeof (int));
- sl2=vcalloc ((CL->S)->nseq, sizeof (int));
-
- norm=vcalloc ( l1+1, sizeof (float));
-
-
- for (a=0;a<ns[0]; a++)sl1[ls[0][a]]=1;
- for (a=0;a<ns[1]; a++)sl2[ls[1][a]]=1;
-
-
-
- used=declare_float (l2+1,2);
- used_list=vcalloc (l2+1, sizeof (int));
- nused=0;
-
- for (p1=0; p1<=l1; p1++)
- {
-
- for (tot=0,nused=0,si=0;p1>0 && si<ns[0]; si++)
- {
- s=ls [0][si];r=pos[s][p1-1];
- for (a=1; r>0 && a<CL->residue_index[s][r][0];a+=ICHUNK)
- {
- t_s=CL->residue_index[s][r][a+SEQ2];
- t_r=CL->residue_index[s][r][a+R2];
- t_w=CL->residue_index[s][r][a+WE];
- if (sl1[t_s])continue;//do not extend within a profile
-
- norm[p1]++;
- for (b=0; b<CL->residue_index[t_s][t_r][0];)
- {
- if (b==0){t_s2=t_s;t_r2=t_r;t_w2=t_w;b++;}
- else
- {
- t_s2=CL->residue_index[t_s][t_r][b+SEQ2];
- t_r2=CL->residue_index[t_s][t_r][b+R2];
- t_w2=CL->residue_index[t_s][t_r][b+WE];
- b+=ICHUNK;
- }
-
- if (sl2[t_s2])
- {
- p2=inv_pos[t_s2][t_r2];
- score=MIN(((float)t_w/(float)NORM_F),((float)t_w2/(float)NORM_F));
-
- if (!used[p2][1] && score>0)
- {
- used_list[nused++]=p2;
- }
-
- tot+=score;
- used[p2][0]+=score;
- used[p2][1]++;
- }
- }
- }
- }
- //FILTER: Keep in the graph the edges where (p1->p2/(Sum (P1->x))>0.01
- filter=0.01;
-
- for (a=0; a<nused; a++)
- {
-
- p2=used_list[a];
- nscore=used[p2][0]/tot; //Normalized score used for filtering
- score =used[p2][0];
- used[p2][0]=used[p2][1]=0;
-
- if (nscore>filter && p1!=0 && p2!=0 && p1!=l1 && p2!=l2)
- {
- score=((norm[p1]>0)?score/norm[p1]:0)*NORM_F;
- addE (p1,p2,((l1-(p1))+(p2)),score,list_in, n_in);
- }
- }
- }
- free_float (used, -1);
- vfree (used_list);
- free_int (inv_pos, -1);
- free_int (pos, -1);
- vfree (sl2);vfree (sl1);
- vfree(norm);
- return n_in[0];
-}
-
-
-
-
-int list2linked_pair_wise ( Alignment *A, int *ns, int **l_s, Constraint_list *CL, int **list, int n, char ***al, int *len);
-int linked_pair_wise ( Alignment *A, int *ns, int **ls, Constraint_list *CL)
-{
- int n=0;
- static int **list=NULL;
- int score, a;
- char **al;
- int len=0;
- int invert=0;
- int tr0,tr1;
-
- if ( !A)free_int (list, -1);
- if ( !CL->residue_index)return myers_miller_pair_wise (A, ns,ls,CL);
-
-
- tr0=ns[0]*strlen (A->seq_al[ls[0][0]]);
- tr1=ns[1]*strlen (A->seq_al[ls[1][0]]);
-
- if (tr0>tr1)
- {
- int *ins;
- int **ils;
- int a,b,c;
- invert=1;
- ins=vcalloc (2, sizeof(int));
- ils=declare_int (2, (CL->S)->nseq);
-
- for ( a=0; a<2; a++)
- {
- ins[a]=ns[a];
- for (b=0; b<ns[a]; b++)ils[a][b]=ls[a][b];
- }
-
- for (c=1,a=0; a<2; a++,c--)
- {
- ns[c]=ins[a];
- for (b=0; b<ins[a]; b++)
- ls[c][b]=ils[a][b];
- }
-
- vfree (ins);
- free_int (ils, -1);
- }
-
-
-
-
-
- /*Prepare the list*/
-
-
- cl2pair_list_ecl_pc (A, ns, ls, CL, &list, &n);
- cl2diag_cap (A, ns, ls, CL, &list, &n);
- cl2list_borders (A, ns, ls, CL, &list, &n);
- list2nodup_list (A, ns, ls, CL, &list, &n);
- /*Do the DP*/
- score=list2linked_pair_wise (A, ns, ls, CL, list, n, &al,&len);
- free_char (al, -1);
-
- if (invert)
- {
- int *ins;
- int **ils;
- int a,b,c;
-
- ins=vcalloc (2, sizeof(int));
- ils=declare_int (2, (CL->S)->nseq);
-
- for ( a=0; a<2; a++)
- {
- ins[a]=ns[a];
- for (b=0; b<ns[a]; b++)ils[a][b]=ls[a][b];
- }
-
- for (c=1,a=0; a<2; a++,c--)
- {
- ns[c]=ins[a];
- for (b=0; b<ins[a]; b++)
- ls[c][b]=ils[a][b];
- }
-
- vfree (ins);
- free_int (ils, -1);
- }
-
- /*Free the list*/
- return score;
-}
-
-#define LIN(a,b,c) a[b*5+c]
-int list2linked_pair_wise( Alignment *A, int *ns, int **l_s, Constraint_list *CL, int **list, int n, char ***tb, int *len)
-{
- int a,b,c, i, j, LEN=0, start_trace;
- int pi, pj,ij, delta_i, delta_j, prev_i, prev_j;
- static int **slist;
- static long *MI, *MJ, *MM,*MT2;
- static int *sortseq;
- static int max_size;
- int gop, gep, igop, igep;
- int l1, l2, l, ls;
- char **al;
- char **aln,*char_buf;
- int ni=0, nj=0;
- long score;
- int nomatch;
-
- l1=strlen (A->seq_al[l_s[0][0]]);
- l2=strlen (A->seq_al[l_s[1][0]]);
- al=declare_char (2,l1+l2+1);
- tb[0]=al;
-
-
- //Penalties: max score is NORM_F
- //Penalties must be negative
- igop=CL->gop;
- gep=igep=CL->gep;
-
- if (n>max_size)
- {
- max_size=n;
-
- vfree (MI);vfree (MJ); vfree (MM);
- free_int (slist, -1);
-
- slist=declare_int (n,3);
-
- MI=vcalloc (5*n, sizeof (long));
- MJ=vcalloc (5*n, sizeof (long));
- MM=vcalloc (5*n, sizeof (long));
-
- }
- else
- {
- for (a=0; a<n; a++)
- for (b=0; b<5; b++)LIN(MI,a,b)=LIN(MJ,a,b)=LIN(MJ,a,b)=-1000000;
- }
-
- /*New Bit: Start*/
- if (!sortseq) sortseq=vcalloc( 7, sizeof (int));
- sortseq[0]=0; sortseq[1]=1;sortseq[2]=-1;
- sort_list_int2 (list, sortseq,7, 0, n-1);
-
- for (a=0; a<n; a++)
- {
-
- slist[a][0]=a;
- list[a][4]=a;
- }
-
- sortseq[0]=1; sortseq[1]=0;sortseq[2]=-1;
- sort_list_int2 (list, sortseq,7, 0, n-1);
- for (a=0; a<n; a++)
- {
- slist[a][1]=list[a][4];
- list[a][5]=a;
- }
-
- sortseq[0]=3; sortseq[1]=0;sortseq[2]=1;sortseq[3]=-1;
- sort_list_int2 (list, sortseq,7, 0, n-1);
- for (a=0; a<n; a++)
- {
- slist[a][2]=list[a][4];
- list[a][6]=a;
- }
-
- sortseq[0]=0; sortseq[1]=1;sortseq[2]=-1;
- sort_list_int2 (list, sortseq,7, 0, n-1);
-
- /*New Bit: EnD*/
-
-
-
-
-
-
- for (a=0; a<n; a++)
- {
-
-
- i=list[a][0];
- j=list[a][1];
-
-
- if (i==l1 || j==l2)gop=0;
- else gop=igop;
-
- if (i==l1 && j==l2)start_trace=a;
- else if ( i==0 || j==0)
- {
- LIN(MM,a,0)=-1000000;
- if (j==0)
- {
-
- LIN(MJ,a,0)=-10000000;
- LIN(MI,a,0)=gep*i;
-
- }
- else if (i==0)
- {
-
- LIN(MI,a,0)=-10000000;
- LIN(MJ,a,0)=gep*j;
-
- }
-
- LIN(MI,a,1)=LIN(MJ,a,1)=-1;
- LIN(MI,a,2)=LIN(MJ,a,2)=i;
- LIN(MI,a,3)=LIN(MJ,a,3)=j;
- continue;
- }
-
- pi=list[a][5];
- pi=slist[pi-1][1];
-
- pj=list[a][4];
- pj=slist[pj-1][0];
-
- ij=list[a][6];
- ij=slist[ij-1][2];
-
-
- ij=list[a][6];
- ij=slist[ij-1][2];
-
-
-
-
-
- prev_i=list[pi][0];
- prev_j=list[pj][1];
-
- delta_i=list[a][0]-list[pi][0];
- delta_j=list[a][1]-list[pj][1];
-
- /*Linear Notation*/
- LIN(MI,a,0)=MAX(LIN(MI,pi,0),(LIN(MM,pi,0)+gop))+delta_i*gep;
- LIN(MI,a,1)=pi;
- LIN(MI,a,2)=delta_i;
- LIN(MI,a,3)=0;
- LIN(MI,a,4)=(LIN(MI,pi,0)>=(LIN(MM,pi,0)+gop))?'i':'m';
-
-
- LIN(MJ,a,0)=MAX(LIN(MJ,pj,0),(LIN(MM,pj,0)+gop))+delta_j*gep;
- LIN(MJ,a,1)=pj;
- LIN(MJ,a,2)=0;
- LIN(MJ,a,3)=delta_j;
-
- LIN(MJ,a,4)=(LIN(MJ,pj,0)>=LIN(MM,pj,0)+gop)?'j':'m';
-
-
-
- if (a>1 && (ls=list[a][0]-list[ij][0])==(list[a][1]-list[ij][1]))
- {
- LIN(MM,a,0)=MAX3(LIN(MM,ij,0),LIN(MI,ij,0),LIN(MJ,ij,0))+list[a][2]-(ls*CL->nomatch);
-
- LIN(MM,a,1)=ij;
- LIN(MM,a,2)=ls;
- LIN(MM,a,3)=ls;
- if ( LIN(MM,ij,0)>=LIN(MI,ij,0) && LIN(MM,ij,0)>=LIN(MJ,ij,0))LIN(MM,a,4)='m';
- else if ( LIN(MI,ij,0) >= LIN(MJ,ij,0))LIN(MM,a,4)='i';
- else LIN(MM,a,4)='j';
-
- }
- else
- {
- LIN(MM,a,0)=UNDEFINED;
- LIN(MM,a,1)=-1;
- }
- }
-
- a=start_trace;
- if (LIN(MM,a,0)>=LIN(MI,a,0) && LIN(MM,a,0) >=LIN(MJ,a,0))MT2=MM;
- else if ( LIN(MI,a,0)>=LIN(MJ,a,0))MT2=MI;
- else MT2=MJ;
-
- score=MAX3(LIN(MM,a,0), LIN(MI,a,0), LIN(MJ,a,0));
-
- i=l1;
- j=l2;
-
-
- while (!(i==0 &&j==0))
- {
- int next_a;
- l=MAX(LIN(MT2,a,2),LIN(MT2,a,3));
- // 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);
- if (i==0)
- {
- while ( j>0)
- {
- al[0][LEN]=0;
- al[1][LEN]=1;
- j--; LEN++;
- }
- }
- else if (j==0)
- {
- while ( i>0)
- {
- al[0][LEN]=1;
- al[1][LEN]=0;
- i--; LEN++;
- }
- }
-
- else if (l==0) {HERE ("L=0 i=%d j=%d",l, i, j);exit (0);}
- else
- {
- for (b=0; b<l; b++, LEN++)
- {
- if (LIN(MT2,a,2)){al[0][LEN]=1;i--;ni++;}
- else al[0][LEN]=0;
-
- if (LIN(MT2,a,3)){al[1][LEN]=1;j--;nj++;}
- else al[1][LEN]=0;
- }
-
- next_a=LIN(MT2,a,1);
- if (LIN(MT2,a,4)=='m')MT2=MM;
- else if (LIN(MT2,a,4)=='i')MT2=MI;
- else if (LIN(MT2,a,4)=='j')MT2=MJ;
- a=next_a;
- }
- }
-
-
-
- invert_list_char ( al[0], LEN);
- invert_list_char ( al[1], LEN);
-
-
- if ( A->declared_len<=LEN)A=realloc_aln ( A,2*LEN+1);
- aln=A->seq_al;
- char_buf= vcalloc (LEN+1, sizeof (char));
-
- for ( c=0; c< 2; c++)
- {
- for ( a=0; a< ns[c]; a++)
- {
- int ch=0;
- for ( b=0; b< LEN; b++)
- {
- if (al[c][b]==1)
- char_buf[b]=aln[l_s[c][a]][ch++];
- else
- char_buf[b]='-';
- }
- char_buf[b]='\0';
- sprintf (aln[l_s[c][a]],"%s", char_buf);
- }
- }
-
- A->len_aln=LEN;
- A->nseq=ns[0]+ns[1];
-
- vfree (char_buf);
- len[0]=LEN;
- return score;
-}
-
-
-
-
-//linked_pair_wise_collapse
-//Collapses the CL as it proceeds during the progressive alignment
-//Cannot be parralelized
-
-
-void display_ns (Alignment *A,Constraint_list *CL, int *ns, int **ls, char *txt);
-int cl2pair_list_collapse ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in);
-Constraint_list* collapse_list (Alignment *A,int *ns, int **ls, char**al, int len, Constraint_list *CL);
-int ns2s (int *ns, int **ls, int *is1, int *is2, int *is);
-int linked_pair_wise_collapse ( Alignment *A, int *ns, int **ls, Constraint_list *CL)
-{
- int n=0;
- static int **list=NULL;
- int score, a;
- char **al;
- int len=0;
-
-
- if ( !A)free_int (list, -1);
- if ( !CL->residue_index)return myers_miller_pair_wise (A, ns,ls,CL);
-
- /*Prepare the list*/
-
-
- cl2pair_list_collapse (A, ns, ls, CL, &list, &n);
-
- cl2diag_cap (A, ns, ls, CL, &list, &n);
- cl2list_borders (A, ns, ls, CL, &list, &n);
- list2nodup_list (A, ns, ls, CL, &list, &n);
-
- /*Do the DP*/
- score=list2linked_pair_wise (A, ns, ls, CL, list, n, &al,&len);
- CL=collapse_list (A,ns, ls, al, len, CL);
- free_char (al, -1);
- /*Free the list*/
- return score;
-}
-
-
-Constraint_list* collapse_list (Alignment *A,int *ns, int **ls, char **al, int len, Constraint_list *CL)
-{
- int s1, s2,s, cs1, cs2, cr1, cr2,l,ll;
- int **lu;
- int a,b,c,d;
- static char *add;
- static int *p;
- FILE *fp;
-
- if (!add)
- {
- add=vtmpnam (NULL);
- p=vcalloc ( 100, sizeof (int));
- }
-
- lu=declare_int (2, len+1);
- for (a=0; a<2; a++)
- for (c=0,b=0; b<len; b++)if (al[a][b]){lu[a][++c]=b+1;}
-
-
- ns2s (ns, ls, &s1, &s2, &s);
-
-
- s1=name_is_in_list (A->name[s1], (CL->S)->name, (CL->S)->nseq, 100);
- s2=name_is_in_list (A->name[s2], (CL->S)->name, (CL->S)->nseq, 100);
- s =name_is_in_list (A->name[s ], (CL->S)->name, (CL->S)->nseq, 100);
-
-
- CL->residue_index[s]=vrealloc (CL->residue_index[s], (len+2)*sizeof (int*));
- for (a=0; a<=len; a++)
- {
- if (!CL->residue_index[s][a])
- {
- CL->residue_index[s][a]=vcalloc (1, sizeof (int));
- CL->residue_index[s][a][0]=1;
- }
- }
-
- fp=vfopen (add, "w");
- CL->ne=0;
- for (cs1=0; cs1<(CL->S)->nseq; cs1++)
- {
- cr1=1;
- while (CL->residue_index[cs1][cr1])
- {
- for (ll=l=1; l<CL->residue_index[cs1][cr1][0]; l+=ICHUNK)
- {
- cs2=CL->residue_index[cs1][cr1][l+SEQ2];
-
- if (cs1==s1 || cs1==s2 || cs2==s1 || cs2==s2)
- {
- p[SEQ1]=cs1;
- p[SEQ2]=CL->residue_index[cs1][cr1][l+SEQ2];
- p[R1] =cr1;
- p[R2] =CL->residue_index[cs1][cr1][l+R2];
- p[CONS]=CL->residue_index[cs1][cr1][l+CONS];
- p[MISC]=CL->residue_index[cs1][cr1][l+MISC];
- p[WE]=CL->residue_index[cs1][cr1][l+WE];
- if (cs1==s1)
- {
- p[SEQ1]=s;
- p[R1]=lu[0][p[R1]];
- }
- else if (cs1==s2)
- {
- p[SEQ1]=s;
- p[R1]=lu[1][p[R1]];
- }
-
- if (cs2==s1)
- {
- p[SEQ2]=s;
- p[R2]=lu[0][p[R2]];
- }
- else if (cs2==s2)
- {
-
- p[SEQ2]=s;
- p[R2]=lu[1][p[R2]];
- }
-
- if (p[SEQ1]==p[SEQ2]);
- else for (d=0; d<CL->entry_len; d++)fprintf (fp, "%d ", p[d]);
- }
- else
- {
- for (d=0; d<ICHUNK; d++) CL->residue_index[cs1][cr1][ll++]=CL->residue_index[cs1][cr1][d+l];
- CL->ne++;
- }
- }
- CL->residue_index[cs1][cr1][0]=ll;
- cr1++;
- }
- }
- vfclose (fp);
- CL=undump_constraint_list (CL,add);
- return CL;
-}
-
-
-
-
-int cl2pair_list_collapse ( Alignment *A, int *ns, int **ls, Constraint_list *CL, int ***list_in, int *n_in)
-{
- int si, r1,r2,t_s, t_r,t_w, t_s2, t_r2, t_w2, s1, s2;
- int a, b, l1, l2;
-
- int nused;
- int *used_list;
-
- float nscore, score, tot, filter, avg=0, new=0;
- float **used;
- int *norm;
-
-
-
- if ( !A) return 0;
-
- ns2s (ns, ls, &s1, &s2,NULL);
-
- l1=strlen (A->seq_al[s1]);
- l2=strlen (A->seq_al[s2]);
- used=declare_float (l2+1,2); used_list=vcalloc (l2+1, sizeof (int));
- nused=0;
- norm=vcalloc (l1+2, sizeof(int));
-
- s1=name_is_in_list (A->name[s1], (CL->S)->name, (CL->S)->nseq, 100);
- s2=name_is_in_list (A->name[s2], (CL->S)->name, (CL->S)->nseq, 100);
-
- for (r1=1; r1<=l1; r1++)
- {
- tot=0; nused=0;
- for (a=1; r1>0 && a<CL->residue_index[s1][r1][0];a+=ICHUNK)
- {
- t_s=CL->residue_index[s1][r1][a+SEQ2];
- t_r=CL->residue_index[s1][r1][a+R2];
- t_w=CL->residue_index[s1][r1][a+WE];
- norm[r1]++;
- for (b=0; b<CL->residue_index[t_s][t_r][0];)
- {
- if (b==0){t_s2=t_s;t_r2=t_r;t_w2=t_w;b++;}
- else
- {
- t_s2=CL->residue_index[t_s][t_r][b+SEQ2];
- t_r2=CL->residue_index[t_s][t_r][b+R2];
- t_w2=CL->residue_index[t_s][t_r][b+WE];
- b+=ICHUNK;
- }
-
- if (t_s2==s2)
- {
- score=MIN(((float)t_w/(float)NORM_F),((float)t_w2/(float)NORM_F));
-
- if (!used[t_r2][1] && score>0)
- {
- used_list[nused++]=t_r2;
- }
-
- tot+=score;
- used[t_r2][0]+=score;
- used[t_r2][1]++;
- }
- }
- }
-
- //FILTER: Keep in the graph the edges where (p1->p2/(Sum (P1->x))>0.01
- filter=0.01;
-
- for (a=0; a<nused; a++)
- {
-
- r2=used_list[a];
- nscore=used[r2][0]/tot; //Normalized score used for filtering
- score =used[r2][0];
-
- used[r2][0]=used[r2][1]=0;
-
- if (nscore>filter && r1!=0 && r2!=0 && r1!=l1 && r2!=l2)
- {
- score=((norm[r1]>0)?score/norm[r1]:0)*NORM_F;
- addE (r1,r2,((l1-(r1))+(r2)),score,list_in, n_in);
- }
- }
- }
-
- free_float (used, -1);
- vfree (used_list);
- vfree (norm);
- return n_in[0];
-}
-int ns2s (int *ns, int **ls, int *is1, int *is2, int *is)
-{
- int a, b;
- int s1, s2, s;
-
- s1=s2=s=-1;
-
- for (a=0; a< 2; a++)
- for (b=0; b<ns[a]; b++)
- {
- if (a==0)s1=MAX(s1,(ls[a][b]));
- if (a==1)s2=MAX(s2,(ls[a][b]));
- }
- s=MAX((s1),(s2));
- if (is1)is1[0]=s1;
- if (is2)is2[0]=s2;
- if (is) is [0]=s;
-
- return s;
-}
-
-/******************************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