A=reorder_aln (A, (CL->S)->name,(CL->S)->nseq);
pos=aln2pos_simple ( A,A->nseq);
- CL=index_res_constraint_list (CL, CL->weight_field);
+
for (p1=0; p1<A->len_aln; p1++)
{
A=reorder_aln (A, (CL->S)->name,(CL->S)->nseq);
pos=aln2pos_simple ( A,A->nseq);
- CL=index_res_constraint_list (CL, CL->weight_field);
+
for (p1=0; p1<A->len_aln; p1++)
{
free_int (pos, -1);
score=(int)(((float)score)/(A->len_aln*SCORE_K));
- score=(int)(CL->L && CL->normalise)?((score*MAXID)/(CL->normalise)):(score);
+ score=(int)(CL->residue_index && CL->normalise)?((score*MAXID)/(CL->normalise)):(score);
return (int)score;
}
int sub_aln2sub_aln_raw_score ( Alignment *A,Constraint_list *CL, const char *mode, int *ns, int **ls)
Alignment *TopS=NULL, *LastS=NULL, *CurrentS=NULL;
-
+
if ( IN->A)IN=IN->A;
while (IN)
{
else if ( CL->ne==0) return matrix_evaluate_output (IN, CL);
}
else if ( strm (mode, "no"))return NULL;
- else if ( strm4 ( mode, "t_coffee_fast","tcoffee_fast","fast_tcoffee", "fast_t_coffee"))
+ else if ( strstr( mode, "triplet"))
+ {
+ return triplet_coffee_evaluate_output ( IN,CL);
+ }
+ else if ( strstr ( mode, "fast"))
{
return fast_coffee_evaluate_output ( IN,CL);
}
- else if ( strm4 ( mode, "t_coffee_slow","tcoffee_slow","slow_tcoffee","slow_t_coffee" ))
+ else if ( strstr ( mode, "slow"))
{
return slow_coffee_evaluate_output ( IN,CL);
}
- else if ( strm4 ( mode, "tcoffee_non_extended","t_coffee_non_extended","non_extended_tcoffee","non_extended_t_coffee"))
+ else if ( strstr ( mode, "non_extended"))
{
return non_extended_t_coffee_evaluate_output ( IN,CL);
}
- else if ( strm5 ( mode, "tcoffee_heuristic","t_coffee_heuristic","heuristic_tcoffee","heuristic_t_coffee", "dali"))
+
+ else if ( strm (mode, "sequences"))
{
- return heuristic_coffee_evaluate_output ( IN,CL);
+ return coffee_seq_evaluate_output ( IN,CL);
}
else
{
vfree(aa);
return OUT;
}
-
+Alignment * fork_triplet_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL);
+Alignment * nfork_triplet_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL);
+Alignment * triplet_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL)
+{
+
+ if (!IN || !CL || !CL->residue_index) return IN;
+
+
+ if ( get_nproc()==1)return nfork_triplet_coffee_evaluate_output (IN,CL);
+ else if (strstr ( CL->multi_thread, "evaluate"))return fork_triplet_coffee_evaluate_output (IN,CL);
+ else return nfork_triplet_coffee_evaluate_output (IN,CL);
+}
+Alignment * fork_triplet_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL)
+ {
+ Alignment *OUT=NULL;
+ int **pos;
+
+ double score_col=0, score_aln=0, score_res=0, score=0;
+ double max_col=0, max_aln=0, max_res=0;
+ double *max_seq, *score_seq;
+
+ int a,b, x, y,res;
+ int s1,r1,s2,r2,w2,s3,r3,w3;
+ int **lu;
+
+ //multi-threading
+ FILE *fp;
+ char **pid_tmpfile;
+ int sjobs, njobs;
+ int **sl;
+ int j;
+
+ OUT=copy_aln (IN, OUT);
+ pos=aln2pos_simple(IN, IN->nseq);
+ sprintf ( OUT->name[IN->nseq], "cons");
+
+ max_seq=vcalloc ( IN->nseq, sizeof (double));
+ score_seq=vcalloc ( IN->nseq, sizeof (double));
+ lu=declare_int (IN->nseq, IN->len_aln+1);
+
+ //multi Threading stuff
+ njobs=get_nproc();
+ sl=n2splits (njobs,IN->len_aln);
+ pid_tmpfile=vcalloc (njobs, sizeof (char*));
+
+ for (sjobs=0,j=0; sjobs<njobs; j++)
+ {
+ pid_tmpfile[j]=vtmpnam(NULL);
+ if (vvfork (NULL)==0)
+ {
+ initiate_vtmpnam(NULL);
+ fp=vfopen (pid_tmpfile[j], "w");
+ score_aln=max_aln=0;
+ for (a=sl[j][0]; a<sl[j][1]; a++)
+ {
+ if (j==0)output_completion (stderr,a,sl[0][1],1, "Final Evaluation");
+ score_col=max_col=0;
+ for (b=0; b<IN->nseq; b++)
+ {
+ s1=IN->order[b][0];
+ r1=pos[b][a];
+ if (r1>=0)lu[s1][r1]=1;
+ }
+ for (b=0; b<IN->nseq; b++)
+ {
+ score_res=max_res=NORM_F;
+ res=NO_COLOR_RESIDUE;
+ s1=IN->order[b][0];
+ r1=pos[b][a];
+ if (r1<=0)
+ {
+ fprintf (fp, "%d ", res);
+ continue;
+ }
+
+ for (x=1;x<CL->residue_index[s1][r1][0];x+=ICHUNK)
+ {
+ s2=CL->residue_index[s1][r1][x+SEQ2];
+ r2=CL->residue_index[s1][r1][x+R2];
+ w2=CL->residue_index[s1][r1][x+WE];
+ for (y=1; y<CL->residue_index[s2][r2][0];y+=ICHUNK)
+ {
+ s3=CL->residue_index[s2][r2][y+SEQ2];
+ r3=CL->residue_index[s2][r2][y+R2];
+ w3=CL->residue_index[s2][r2][y+WE];
+
+ score=MIN(w2,w3);
+
+ max_res+=score;
+ max_col+=score;
+ max_seq[b]+=score;
+ max_aln+=score;
+
+ if (lu[s3][r3])
+ {
+ score_res+=score;
+ score_col+=score;
+ score_seq[b]+=score;
+ score_aln+=score;
+ }
+ }
+ }
+ res=(max_res==0)?NO_COLOR_RESIDUE:((score_res*10)/max_res);
+ res=(res==NO_COLOR_RESIDUE)?res:(MIN(res,9));
+ fprintf ( fp, "%d ", res);
+ }
+ for (b=0; b<IN->nseq; b++)
+ {
+ s1=IN->order[b][0];
+ r1=pos[b][a];
+ if (r1>0)lu[s1][r1]=0;
+ }
+
+ res=(max_col==0)?NO_COLOR_RESIDUE:((score_col*10)/max_col);
+ res=(res==NO_COLOR_RESIDUE)?res:(MIN(res,9));
+ fprintf (fp, "%d ", res);
+ for (b=0; b<IN->nseq; b++)fprintf (fp, "%f %f ", score_seq[b], max_seq[b]);
+ }
+ fprintf (fp, "%f %f ", score_aln, max_aln);
+ 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++)
+ {
+ float sseq, mseq;
+ fp=vfopen (pid_tmpfile[j], "r");
+ for (a=sl[j][0];a<sl[j][1]; a++)
+ {
+ for (b=0; b<=IN->nseq; b++)//don't forget the consensus
+ {
+ fscanf (fp, "%d ", &res);
+ OUT->seq_al[b][a]=res;
+ }
+ for (b=0; b<IN->nseq; b++)
+ {
+ fscanf (fp, "%f %f", &sseq, &mseq);
+ score_seq[b]+=sseq;
+ max_seq[b]+=mseq;
+ }
+ }
+ fscanf (fp, "%f %f", &sseq, &mseq);
+ score_aln+=sseq;
+ max_aln+=mseq;
+ vfclose (fp);
+ remove (pid_tmpfile[j]);
+ }
+ fprintf ( stderr, "\n");
+
+ IN->score_aln=OUT->score_aln=(max_aln==0)?0:((score_aln*100)/max_aln);
+ for ( a=0; a< OUT->nseq; a++)
+ {
+ OUT->score_seq[a]=(max_seq[a]==0)?0:((score_seq[a]*100)/max_seq[a]);
+ }
+
+ free_int (lu,-1);
+ free_int (pos , -1);
+ vfree (pid_tmpfile);
+ free_int (sl, -1);
+ vfree ( score_seq);
+ vfree ( max_seq);
+ return OUT;
+ }
+
+Alignment * nfork_triplet_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL)
+ {
+ Alignment *OUT=NULL;
+ int **pos;
+
+ double score_col=0, score_aln=0, score_res=0, score=0;
+ double max_col=0, max_aln=0, max_res=0;
+ double *max_seq, *score_seq;
+
+ int a,b, x, y,res;
+ int s1,r1,s2,r2,w2,s3,r3,w3;
+ int **lu;
+
+
+
+ OUT=copy_aln (IN, OUT);
+ pos=aln2pos_simple(IN, IN->nseq);
+ sprintf ( OUT->name[IN->nseq], "cons");
+
+ max_seq=vcalloc ( IN->nseq, sizeof (double));
+ score_seq=vcalloc ( IN->nseq, sizeof (double));
+ lu=declare_int (IN->nseq, IN->len_aln+1);
+
+
+
+
+ score_aln=max_aln=0;
+ for (a=0; a<IN->len_aln; a++)
+ {
+ output_completion (stderr,a,IN->len_aln,1, "Final Evaluation");
+ score_col=max_col=0;
+ for (b=0; b<IN->nseq; b++)
+ {
+ s1=IN->order[b][0];
+ r1=pos[b][a];
+ if (r1>=0)lu[s1][r1]=1;
+ }
+ for (b=0; b<IN->nseq; b++)
+ {
+ score_res=max_res=NORM_F;
+ OUT->seq_al[b][a]=NO_COLOR_RESIDUE;
+ s1=IN->order[b][0];
+ r1=pos[b][a];
+ if (r1<=0)continue;
+ for (x=1;x<CL->residue_index[s1][r1][0];x+=ICHUNK)
+ {
+ s2=CL->residue_index[s1][r1][x+SEQ2];
+ r2=CL->residue_index[s1][r1][x+R2];
+ w2=CL->residue_index[s1][r1][x+WE];
+ for (y=1; y<CL->residue_index[s2][r2][0];y+=ICHUNK)
+ {
+ s3=CL->residue_index[s2][r2][y+SEQ2];
+ r3=CL->residue_index[s2][r2][y+R2];
+ w3=CL->residue_index[s2][r2][y+WE];
+
+ score=MIN(w2,w3);
+
+ max_res+=score;
+ max_col+=score;
+ max_seq[b]+=score;
+ max_aln+=score;
+
+ if (lu[s3][r3])
+ {
+ score_res+=score;
+ score_col+=score;
+ score_seq[b]+=score;
+ score_aln+=score;
+ }
+ }
+ }
+ res=(max_res==0)?NO_COLOR_RESIDUE:((score_res*10)/max_res);
+ res=(res==NO_COLOR_RESIDUE)?res:(MIN(res,9));
+ OUT->seq_al[b][a]=res;
+ }
+ for (b=0; b<IN->nseq; b++)
+ {
+ s1=IN->order[b][0];
+ r1=pos[b][a];
+ if (r1>0)lu[s1][r1]=0;
+ }
+
+ res=(max_col==0)?NO_COLOR_RESIDUE:((score_col*10)/max_col);
+ res=(res==NO_COLOR_RESIDUE)?res:(MIN(res,9));
+ OUT->seq_al[IN->nseq][a]=res;
+ }
+ fprintf ( stderr, "\n");
+
+ IN->score_aln=OUT->score_aln=(max_aln==0)?0:((score_aln*100)/max_aln);
+ for ( a=0; a< OUT->nseq; a++)
+ {
+ OUT->score_seq[a]=(max_seq[a]==0)?0:((score_seq[a]*100)/max_seq[a]);
+ }
+
+ free_int (lu,-1);
+ free_int (pos , -1);
+ vfree ( score_seq);
+ vfree ( max_seq);
+ return OUT;
+ }
Alignment * fast_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL)
{
int a,b, c, m,res, s, s1, s2, r1, r2;
vfree ( score_seq);
vfree ( max_seq);
return OUT;
- }
-
-
+ }
Alignment * slow_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL)
{
return OUT;
}
-
-Alignment * heuristic_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL)
+double genepred2sd (Sequence *S);
+double genepred2avg (Sequence *S);
+double genepred2zsum2 (Sequence *S);
+int * genepred2orf_len (Sequence *S);
+double genepred2acc (Sequence *S, Sequence *PS);
+double seq_genepred2acc (Sequence *S, Sequence *PS, char *name);
+Alignment *coffee_seq_evaluate_output ( Alignment *IN, Constraint_list *CL)
+{
+ int **w, a,b,c,l,i;
+ int avg, min_avg, best_cycle;
+ char **pred_list,**weight_list;
+ FILE *fp;
+ Sequence *S, *PS;
+ int ncycle=20;
+ float nsd, best_score,score, acc;
+ int min_ne;
+ Alignment *A;
+ float *count;
+ min_ne=CL->ne/10; //avoid side effects due to very relaxed libraries
+ w=NULL;
+ pred_list=declare_char (ncycle, 100);
+ weight_list=declare_char(ncycle,100);
+ S=CL->S;
+ fprintf ( stderr, "\nSCORE_MODE: %s\n", CL->genepred_score);
+ for (a=0; a<ncycle && CL->ne>min_ne; a++)
{
- int a,b, c,res, s, s1, s2, r1, r2;
- Alignment *OUT=NULL;
- int **pos;
- int max_score_r, score_r;
- double score_col=0, score_aln=0;
- int max_score_col, max_score_aln;
- double *max_score_seq, *score_seq;
- int **tot_extended_weight;
- int **res_extended_weight;
- int n_res_in_col;
+ if (w);free_int (w, -1);
+ w=list2residue_total_weight (CL);
+ PS=dnaseq2geneseq (S, w);
+ nsd=(float)(genepred2sd(PS)/genepred2avg(PS));
+ if ( strm (CL->genepred_score, "nsd"))score=nsd;
+ else score=1-seq_genepred2acc (S, PS, CL->genepred_score);
+ fprintf ( stderr, "Cycle: %2d LIB: %6d AVG LENGTH: %6.2f NSD: %8.2f SCORE: %.2f ",a+1,CL->ne,(float) genepred2avg(PS),nsd,score);
+
+ vfree (display_accuracy (genepred_seq2accuracy_counts (S, PS, NULL),stderr));
+ pred_list[a]=vtmpnam(NULL);
+ weight_list[a]=vtmpnam(NULL);
+ output_fasta_seqS(pred_list[a],PS);
+ output_array_int (weight_list[a],w);
+ free_sequence (PS, -1);
+
+ if (a==0 || score<best_score)
+ {
+ best_score=score;
+ best_cycle=a;
+ }
+ CL=relax_constraint_list_4gp(CL);
+ }
+ if (w)free_int(w,-1);
+ S=get_fasta_sequence (pred_list[best_cycle], NULL);
- /*
- Residue x: sum of observed extended X.. /sum of possible X..
- */
+ fprintf ( stderr, "\nSelected Cycle: %d (SCORE=%.4f)----> ", best_cycle+1, best_score);
+ vfree (display_accuracy (genepred_seq2accuracy_counts ((CL->S), S, NULL),stderr));
- if ( !CL->evaluate_residue_pair){fprintf ( stderr, "\nWARNING: CL->evaluate_residue_pair Not set\nSet to: extend_residue_pair\n");CL->evaluate_residue_pair= extend_residue_pair; }
-
- OUT=copy_aln (IN, OUT);
- pos=aln2pos_simple(IN, IN->nseq);
+ genepred_seq2accuracy_counts4all((CL->S), S);
+
+ IN=seq2aln (S,NULL, 1);
+ IN->score_res=input_array_int (weight_list[best_cycle]);
+ return IN;
+}
+double seq_genepred2acc (Sequence *S, Sequence *PS, char *name)
+{
+ float *count;
+ float *acc;
+ float r;
+ int ref, target;
- max_score_seq=vcalloc ( IN->nseq, sizeof (double));
- score_seq=vcalloc ( IN->nseq, sizeof (double));
+ if ( strm (name, "best"))return genepred2acc (S, PS);
+
+ ref=name_is_in_list (name, S->name, S->nseq, 100);
+ target=name_is_in_list (name, PS->name, PS->nseq, 100);
+
+ if ( target==-1 || ref==-1)
+ {
+ printf_exit (EXIT_FAILURE,stderr, "\nERROR: %s is not a valid sequence", name);
+ }
+ count=genepred2accuracy_counts (S->seq[ref],PS->seq[target],NULL);
+ acc=counts2accuracy (count);
+
+ r=acc[3];
+ vfree (acc);
+ vfree (count);
+ return r;
+}
+
+
- tot_extended_weight=list2residue_partial_extended_weight(CL);
- res_extended_weight=declare_int ((CL->S)->nseq, (CL->S)->max_len+1);
-
- for (a=0; a< IN->len_aln; a++)
- {
- for ( b=0; b< IN->nseq-1; b++)
- {
- s1=IN->order[b][0];
- r1=pos[b][a];
- for ( c=b+1; c< IN->nseq; c++)
- {
- s2=IN->order[c][0];
- r2=pos[c][a];
- if ( s1==s2 && !CL->do_self)continue;
- else if ( r1<=0 || r2<=0) continue;
- else
- {
- if ( s1< s2)s=(CL->evaluate_residue_pair)( CL, s1, r1, s2, r2);
- else s=(CL->evaluate_residue_pair)( CL, s2, r2, s1, r1);
- res_extended_weight[s1][r1]+=s;
- res_extended_weight[s2][r2]+=s;
- }
- }
- }
- }
-
+double genepred2acc (Sequence *S, Sequence *PS)
+{
+ float *count;
+ float *acc;
+ float r;
+ count=genepred_seq2accuracy_counts (S, PS, NULL);
+ acc=counts2accuracy (count);
- sprintf ( OUT->name[IN->nseq], "cons");
- for ( max_score_aln=0,score_aln=0,a=0; a< IN->len_aln; a++)
+ r=acc[3];
+ vfree (acc);
+ vfree (count);
+ return r;
+}
+
+double genepred2sd (Sequence *S)
+{
+ double sum=0, sum2=0;
+ int a, b;
+ int *len;
+ len=genepred2orf_len (S);
+
+ for (a=0; a<S->nseq; a++)
+ {
+ sum+=(double)len[a];
+ sum2+=(double)len[a]*(double)len[a];
+ }
+ sum/=(double)S->nseq;
+ sum2/=S->nseq;
+ vfree (len);
+ return sqrt (sum2-(sum*sum));
+}
+double genepred2avg (Sequence *S)
+{
+ int a, b;
+ double avg=0;
+ int *len;
+ len=genepred2orf_len (S);
+
+ for (a=0; a<S->nseq; a++)avg+=len[a];
+ vfree (len);
+ return avg/(double)S->nseq;
+}
+double genepred2zsum2 (Sequence *S)
+{
+ double sum=0, sum2=0, zscore=0, zsum2=0;
+ int a, b;
+ int *len;
+ double sd, avg;
+ sd=genepred2sd(S);
+ avg=genepred2avg (S);
+
+ len=genepred2orf_len (S);
+ for (a=0; a<S->nseq; a++)
+ {
+ zscore=((double)len[a]-avg)/sd;
+ zsum2+=FABS(zscore);
+ }
+ zsum2/=(float)S->nseq;
+ vfree (len);
+ return zsum2;
+}
+int *genepred2orf_len (Sequence *S)
+{
+ int a,b, *len;
+ len=vcalloc (S->nseq, sizeof (int));
+ for (a=0; a<S->nseq; a++)
+ for (b=0; b<S->len[a]; b++)
+ len[a]+=(isupper(S->seq[a][b]));
+ return len;
+}
+
+Alignment *coffee_seq_evaluate_output_old2 ( Alignment *IN, Constraint_list *CL)
+{
+ int **w, a,b,c,l,i;
+ int avg, min_avg, best_cycle;
+ char **pred_list; FILE *fp;
+ Sequence *S, *PS;
+ int ncycle=100;
+ w=NULL;
+ pred_list=declare_char (ncycle, 100);
+ S=CL->S;
+
+ //CL=expand_constraint_list_4gp(CL);
+ min_avg=constraint_list2avg(CL);
+
+ for (a=1; a<ncycle && CL->ne>0; a++)
+ {
+ int t;
+ if (w);free_int (w, -1);
+ w=list2residue_total_weight (CL);
+ CL=relax_constraint_list_4gp (CL);
+ fprintf (stderr,"\nRELAX CYCLE: %2d AVG: %5d [%10d] ", a, avg=constraint_list2avg(CL), CL->ne);
+ for (b=0; b<S->nseq; b++)for (c=0; c<S->len[b]; c++)w[b][c]-=1;
+
+ //rescale nuclotide coding weights
+
+ PS=dnaseq2geneseq (S, w);
+ vfree (display_accuracy (genepred_seq2accuracy_counts (S, PS, NULL),stderr));
+
+
+ free_sequence (PS, PS->nseq);
+ if ( avg<min_avg)
{
- OUT->seq_al[IN->nseq][a]=NO_COLOR_RESIDUE;
- for ( n_res_in_col=0,b=0; b<IN->nseq; b++){n_res_in_col+=(pos[b][a]>0)?1:0;}
- for ( max_score_col=0, score_col=0,b=0; b< IN->nseq; b++)
- {
- OUT->seq_al[b][a]=NO_COLOR_RESIDUE;
- s1=IN->order[b][0];
- r1=pos[b][a];
- if (r1<=0)continue;
- else
- {
- max_score_r =tot_extended_weight[s1][r1];
- score_r=res_extended_weight[s1][r1];
- res=(max_score_r==0 ||n_res_in_col<2 )?NO_COLOR_RESIDUE:((score_r*10)/max_score_r);
- (OUT)->seq_al[b][a]=(res==NO_COLOR_RESIDUE)?res:(MIN(res, 9));
- max_score_col+=max_score_r;
- score_col+=score_r;
- max_score_seq[b]+=max_score_r;
- score_seq[b]+=score_r;
- max_score_aln+=max_score_r;
- score_aln+=score_r;
- }
- res=(max_score_col==0 || n_res_in_col<2)?NO_COLOR_RESIDUE:((score_col*10)/max_score_col);
- OUT->seq_al[IN->nseq][a]=(res==NO_COLOR_RESIDUE)?res:(MIN(res,9));
- }
+ min_avg=avg;
+ best_cycle=a;
}
- IN->score_aln=OUT->score_aln=MIN(100,((max_score_aln==0)?0:((score_aln*100)/max_score_aln)));
- for ( a=0; a< OUT->nseq; a++)
- {
- OUT->score_seq[a]=MIN(100,((max_score_seq[a]==0)?0:((score_seq[a]*100)/max_score_seq[a])));
- }
-
-
- vfree ( score_seq);
- vfree ( max_score_seq);
-
- free_int (tot_extended_weight, -1);
- free_int (res_extended_weight, -1);
- free_int (pos, -1);
-
- return OUT;
}
+
+ S=get_fasta_sequence (pred_list[best_cycle], NULL);
+ vfree (display_accuracy (genepred_seq2accuracy_counts ((CL->S), S, NULL),stderr));myexit (0);
+ for (a=0; a<S->nseq; a++)
+ HERE (">%s\n%s", S->name[a], S->seq[a]);
+ myexit (0);
+ return seq2aln (S,NULL, 1);
+}
+
+
+
Alignment * non_extended_t_coffee_evaluate_output ( Alignment *IN,Constraint_list *CL)
{
int a,b, c,res, s1, s2, r1, r2;
int p;
int max_score=0;
- entry=vcalloc (CL->entry_len, CL->el_size);
+ entry=vcalloc (CL->entry_len+1, CL->el_size);
if ( !CL->evaluate_residue_pair){fprintf ( stderr, "\nWARNING: CL->evaluate_residue_pair Not set\nSet to: extend_residue_pair\n");CL->evaluate_residue_pair= extend_residue_pair; }
OUT=copy_aln (IN, OUT);
/* FUNCTIONS FOR GETING THE COST : (Sequences) ->evaluate_residue_pair */
/* */
/*********************************************************************************************/
+int initialize_scoring_scheme (Constraint_list *CL)
+{
+ if (!CL) return 0;
+ else if (!CL->evaluate_residue_pair)return 0;
+ else if ( !CL->S) return 0;
+ else if ( (CL->S)->nseq<2) return 0;
+ else if ( strlen ((CL->S)->seq[0])==0)return 0;
+ else if ( strlen ((CL->S)->seq[1])==0)return 0;
+ else
+ {
+ //CL->evaluate_residue_pair (CL,13,1,16,1);
+ CL->evaluate_residue_pair (CL,0,1,1,1);
+
+ }
+ return 1;
+}
+
int evaluate_blast_profile_score (Constraint_list *CL, int s1, int r1, int s2, int r2)
{
Alignment *PRF1;
fp=vfopen ( name, "r");
while ( fscanf (fp, "%s %d %c %f\n",b1, &a, &c,&f )==4)
{
- if ( c!=(CL->S)->seq[s1][n]){HERE ("ERROR: %c %c", c,(CL->S)->seq[s1][n] );exit (0);}
+ if ( c!=(CL->S)->seq[s1][n]){HERE ("ERROR: %c %c", c,(CL->S)->seq[s1][n] );myexit (0);}
else array[n++]=(int)(float)100*(float)f;
}
vfclose (fp);
int evaluate_tm_matrix_score ( Constraint_list *CL, int s1, int r1, int s2, int r2)
{
static char **st;
- static int **tmat;
+ float F, RF;
-
- if (!tmat)
+ if (!st)
{
- tmat=read_matrice ("blosum62mt");
+ st= vcalloc ((CL->S)->nseq, sizeof (char*));
+ RF=atoigetenv ("TM_FACTOR_4_TCOFFEE");
+ if ( !RF)RF=10;
}
- if (!st) st= vcalloc ((CL->S)->nseq, sizeof (char*));
if (!st[s1])st[s1]=seq2T_template_string((CL->S),s1);
if (!st[s2])st[s2]=seq2T_template_string((CL->S),s2);
-
-
-
- if (r1>0 && r2>0)
- {
- int p1, p2;
-
- r1--;
- r2--;
- p1=p2=-1;
-
- if (st[s1])p1=tolower (st[s1][r1]);
- if (st[s2])p2=tolower (st[s2][r2]);
-
- if ( p1=='h' && p2=='h')return tmat[(CL->S)->seq[s1][r1]-'A'][(CL->S)->seq[s2][r2]-'A']*SCORE_K;
- //else if (p1=='h' || p2=='h' ) return -100*SCORE_K;
- else return CL->M[(CL->S)->seq[s1][r1]-'A'][(CL->S)->seq[s2][r2]-'A']*SCORE_K;
-
- }
- else
- {
- return 0;
- }
-}
-int evaluate_ssp_matrix_score ( Constraint_list *CL, int s1, int r1, int s2, int r2)
-{
- static char **st;
- static int **alpha, **beta, **coil;
-
-
- if (!alpha)
- {
- beta=read_matrice ("beta_mat");
- alpha=read_matrice ("alpha_mat");
- coil=read_matrice ("coil_mat");
- }
-
- if (!st) st= vcalloc ((CL->S)->nseq, sizeof (char*));
- if (!st[s1])st[s1]=seq2E_template_string((CL->S),s1);
- if (!st[s2])st[s2]=seq2E_template_string((CL->S),s2);
-
-
- if ( !st[s1])HERE ("1******");
- if ( !st[s2])HERE ("2******");
-
- if (r1>0 && r2>0)
- {
- char p1, p2;
- float F;
-
- int score=0;
- p1=p2=-1;
- r1--;
- r2--;
- if (st[s1])p1=tolower (st[s1][r1]);
- if (st[s2])p2=tolower (st[s2][r2]);
-
- if (p1!=-1 && p1==p2)F=1.3;
- else F=1;
-
+ if (r1<=0 || r2<=0)return 0;
+ else if (st[s1] && st[s2] && (st[s1][r1-1]==st[s2][r2-1]))F=RF;
+ else if (st[s1] && st[s2] && (st[s1][r1-1]!=st[s2][r2-1]))F=0;
+ else F=0;
- score= CL->M[(CL->S)->seq[s1][r1]-'A'][(CL->S)->seq[s2][r2]-'A']*F*SCORE_K;
-
-
- return score;
- }
+ return (int)(evaluate_matrix_score (CL, s1, r1, s2, r2))+F;
+}
+int evaluate_ssp_matrix_score ( Constraint_list *CL, int s1, int r1, int s2, int r2)
+{
+ static char **st;
+ float F, RF;
- else
+ if (!st)
{
- return 0;
+ st= vcalloc ((CL->S)->nseq, sizeof (char*));
+ RF=atoigetenv ("SSP_FACTOR_4_TCOFFEE");
+ if ( !RF)RF=10;
}
+
+ if (!st[s1])st[s1]=seq2T_template_string((CL->S),s1);
+ if (!st[s2])st[s2]=seq2T_template_string((CL->S),s2);
+
+ if (r1<=0 || r2<=0)return 0;
+ else if (st[s1] && st[s2] && (st[s1][r1-1]==st[s2][r2-1]))F=RF;
+ else if (st[s1] && st[s2] && (st[s1][r1-1]!=st[s2][r2-1]))F=0;
+ else F=0;
+
+ return (int)(evaluate_matrix_score (CL, s1, r1, s2, r2))+F;
}
for ( a=0; a< (CL->S)->nseq; a++)hasch[a]=vcalloc ( (CL->S)->len[a]+1, sizeof (int));
}
- CL=index_res_constraint_list ( CL, field);
+
hasch[s1][r1]=100000;
- for (a=1; a< CL->residue_index[s1][r1][0]; a+=3)
+ for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s1][r1][a];
- t_r=CL->residue_index[s1][r1][a+1];
- t_w=CL->residue_index[s1][r1][a+2];
+ 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];
if ( CL->seq_for_quadruplet[t_s])
{
- for ( b=1; b<CL->residue_index[t_s][t_r][0]; b+=3)
+ for ( b=1; b<CL->residue_index[t_s][t_r][0]; b+=ICHUNK)
{
- q_s=CL->residue_index[t_s][t_r][b];
- q_r=CL->residue_index[t_s][t_r][b+1];
- q_w=CL->residue_index[t_s][t_r][b+2];
+ q_s=CL->residue_index[t_s][t_r][b+SEQ2];
+ q_r=CL->residue_index[t_s][t_r][b+R2];
+ q_w=CL->residue_index[t_s][t_r][b+WE];
if (CL-> seq_for_quadruplet[q_s])
- hasch[q_s][q_r]=MIN(q_w,t_w);
+ hasch[q_s][q_r]=MIN(q_w,t_w);
}
}
}
- for (a=1; a< CL->residue_index[s2][r2][0]; a+=3)
+ for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s2][r2][a];
- t_r=CL->residue_index[s2][r2][a+1];
- t_w=CL->residue_index[s2][r2][a+2];
+ t_s=CL->residue_index[s2][r2][a+SEQ2];
+ t_r=CL->residue_index[s2][r2][a+R2];
+ t_w=CL->residue_index[s2][r2][a+WE];
if ( CL->seq_for_quadruplet[t_s])
{
- for ( b=1; b<CL->residue_index[t_s][t_r][0]; b+=3)
+ for ( b=1; b<CL->residue_index[t_s][t_r][0]; b+=ICHUNK)
{
- q_s=CL->residue_index[t_s][t_r][b];
- q_r=CL->residue_index[t_s][t_r][b+1];
- q_w=CL->residue_index[t_s][t_r][b+2];
+ q_s=CL->residue_index[t_s][t_r][b+SEQ2];
+ q_r=CL->residue_index[t_s][t_r][b+R2];
+ q_w=CL->residue_index[t_s][t_r][b+WE];
if (hasch[q_s][q_r] && CL->seq_for_quadruplet[q_s])
score+=MIN(hasch[q_s][q_r],MIN(q_w,t_w));
}
score=(CL->normalise)?((score*CL->normalise)/CL->max_ext_value):score;
- for (a=1; a< CL->residue_index[s1][r1][0]; a+=3)
+ for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s1][r1][a];
- t_r=CL->residue_index[s1][r1][a+1];
- t_w=CL->residue_index[s1][r1][a+2];
+ 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];
if ( CL->seq_for_quadruplet[t_s])
{
- for ( b=1; b<CL->residue_index[t_s][t_r][0]; b+=3)
+ for ( b=1; b<CL->residue_index[t_s][t_r][0]; b+=ICHUNK)
{
- q_s=CL->residue_index[t_s][t_r][b];
- q_r=CL->residue_index[t_s][t_r][b+1];
+ q_s=CL->residue_index[t_s][t_r][b+SEQ2];
+ q_r=CL->residue_index[t_s][t_r][b+R2];
hasch[q_s][q_r]=0;
}
}
if (!R)
{
-
+
R=read_rna_lib (CL->S, CL->rna_lib);
rna_lib_extension (CL,R);
n1=n2=0;
list1[n1++]=r1;
- for (a=1; a<R->residue_index[s1][r1][0]; a+=3)
+ for (a=1; a<R->residue_index[s1][r1][0]; a+=ICHUNK)
{
- list1[n1++]=R->residue_index[s1][r1][a+1];
+ list1[n1++]=R->residue_index[s1][r1][a+R2];
}
list2[n2++]=r2;
- for (a=1; a<R->residue_index[s2][r2][0]; a+=3)
+ for (a=1; a<R->residue_index[s2][r2][0]; a+=ICHUNK)
{
- list2[n2++]=R->residue_index[s2][r2][a+1];
+ list2[n2++]=R->residue_index[s2][r2][a+R2];
}
score=residue_pair_extended_list ( CL, s1,list1[0], s2,list2[0]);
-
+
for (score2=0,a=1; a<n1; a++)
for ( b=1; b<n2; b++)
score2=MAX(residue_pair_extended_list ( CL, s1,list1[a], s2,list2[b]), score2);
-
+
if ( n1==1 || n2==1);
else score=MAX(score,score2);
{
R=read_constraint_list_file (R, list[a]);
}
- R=index_res_constraint_list (R, CL->weight_field);
}
n1=n2=0;
list1[n1++]=r1;
- for (a=1; a<R->residue_index[s1][r1][0]; a+=3)
+ for (a=1; a<R->residue_index[s1][r1][0]; a+=ICHUNK)
{
- list1[n1++]=R->residue_index[s1][r1][a+1];
+ list1[n1++]=R->residue_index[s1][r1][a+R2];
}
list2[n2++]=r2;
- for (a=1; a<R->residue_index[s2][r2][0]; a+=3)
+ for (a=1; a<R->residue_index[s2][r2][0]; a+=ICHUNK)
{
- list2[n2++]=R->residue_index[s2][r2][a+1];
+ list2[n2++]=R->residue_index[s2][r2][a+R2];
}
*/
field=CL->weight_field;
-
+ field=WE;
+
if ( r1<=0 || r2<=0)return 0;
if ( !hasch || max_len!=(CL->S)->max_len)
{
hasch=declare_int ( (CL->S)->nseq, (CL->S)->max_len+1);
}
- CL=index_res_constraint_list ( CL, field);
+
/* Check matches for R1 in the indexed lib*/
hasch[s1][r1]=FORBIDEN;
- for (a=1; a< CL->residue_index[s1][r1][0]; a+=3)
+ for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s1][r1][a];
- t_r=CL->residue_index[s1][r1][a+1];
- hasch[t_s][t_r]=CL->residue_index[s1][r1][a+2];
+ t_s=CL->residue_index[s1][r1][a+SEQ2];
+ t_r=CL->residue_index[s1][r1][a+R2];
+ hasch[t_s][t_r]=CL->residue_index[s1][r1][a+field];
}
/*Check Matches for r1 <-> r2 in the indexed lib */
- for (a=1; a< CL->residue_index[s2][r2][0]; a+=3)
+ for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s2][r2][a];
- t_r=CL->residue_index[s2][r2][a+1];
+ t_s=CL->residue_index[s2][r2][a+SEQ2];
+ t_r=CL->residue_index[s2][r2][a+R2];
if (hasch[t_s][t_r])
{
if (hasch[t_s][t_r]==FORBIDEN)
{
- score+=CL->residue_index[s2][r2][a+2];
+ score+=CL->residue_index[s2][r2][a+field];
}
else
{
- delta=MIN(hasch[t_s][t_r],CL->residue_index[s2][r2][a+2]);
+ delta=MIN(hasch[t_s][t_r],CL->residue_index[s2][r2][a+field]);
score+=delta;
}
}
clean_residue_pair_hasch ( s1, r1,s2, r2, hasch, CL);
return score;
}
-int residue_pair_extended_list_pc ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
+int residue_pair_extended_list_4gp ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
{
double score=0;
*/
field=CL->weight_field;
+ field=WE;
if ( r1<=0 || r2<=0)return 0;
if ( !hasch || max_len!=(CL->S)->max_len)
hasch=declare_int ( (CL->S)->nseq, (CL->S)->max_len+1);
}
- CL=index_res_constraint_list ( CL, field);
+
/* Check matches for R1 in the indexed lib*/
hasch[s1][r1]=FORBIDEN;
- for (a=1; a< CL->residue_index[s1][r1][0]; a+=3)
+ for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s1][r1][a];
- t_r=CL->residue_index[s1][r1][a+1];
- hasch[t_s][t_r]=CL->residue_index[s1][r1][a+2];
+ t_s=CL->residue_index[s1][r1][a+SEQ2];
+ t_r=CL->residue_index[s1][r1][a+R2];
+ hasch[t_s][t_r]=CL->residue_index[s1][r1][a+field];
}
/*Check Matches for r1 <-> r2 in the indexed lib */
- for (a=1; a< CL->residue_index[s2][r2][0]; a+=3)
+ for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s2][r2][a];
- t_r=CL->residue_index[s2][r2][a+1];
+ t_s=CL->residue_index[s2][r2][a+SEQ2];
+ t_r=CL->residue_index[s2][r2][a+R2];
if (hasch[t_s][t_r])
}
else
{
- //delta=((float)hasch[t_s][t_r]/NORM_F)*((float)CL->residue_index[s2][r2][a+2]/NORM_F);
- delta=MIN((((float)hasch[t_s][t_r]/NORM_F)),(((float)CL->residue_index[s2][r2][a+2]/NORM_F)));
+ delta=MAX((((float)hasch[t_s][t_r]/NORM_F)),(((float)CL->residue_index[s2][r2][a+field]/NORM_F)));
score+=delta;
}
}
clean_residue_pair_hasch ( s1, r1,s2, r2, hasch, CL);
score/=(CL->S)->nseq;
- return score*NORM_F;
+ score *=NORM_F;
+
+ return score;
}
-
-int residue_pair_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
+int residue_pair_extended_list_pc ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
{
- double score=0;
- double max_score=0;
- double max_val=0;
-
- int a, t_s, t_r;
- static int **hasch;
- static int max_len;
- int field;
-
- /*
- new function: self normalized
+ //old norm does not take into account the number of effective intermediate sequences
+
+ double score=0;
+ int a, t_s, t_r;
+ static int **hasch;
+ int nclean;
+ static int max_len;
+ int field;
+ double delta;
+ float norm1=0;
+ float norm2=0;
+
+ /*
+
function documentation: start
-
+
int residue_pair_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2);
Computes the extended score for aligning residue seq1(r1) Vs seq2(r2)
*/
field=CL->weight_field;
-
+ field=WE;
+
if ( r1<=0 || r2<=0)return 0;
if ( !hasch || max_len!=(CL->S)->max_len)
{
max_len=(CL->S)->max_len;
if ( hasch) free_int ( hasch, -1);
hasch=declare_int ( (CL->S)->nseq, (CL->S)->max_len+1);
+
}
-
- CL=index_res_constraint_list ( CL, field);
+
/* Check matches for R1 in the indexed lib*/
hasch[s1][r1]=FORBIDEN;
- for (a=1; a< CL->residue_index[s1][r1][0]; a+=3)
+ for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s1][r1][a];
- t_r=CL->residue_index[s1][r1][a+1];
- hasch[t_s][t_r]=CL->residue_index[s1][r1][a+2];
- max_score+=CL->residue_index[s1][r1][a+2];
+ t_s=CL->residue_index[s1][r1][a+SEQ2];
+ t_r=CL->residue_index[s1][r1][a+R2];
+ hasch[t_s][t_r]=CL->residue_index[s1][r1][a+field];
+ norm1++;
+
}
-
+
+
/*Check Matches for r1 <-> r2 in the indexed lib */
- for (a=1; a< CL->residue_index[s2][r2][0]; a+=3)
+ for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s2][r2][a];
- t_r=CL->residue_index[s2][r2][a+1];
-
+ t_s=CL->residue_index[s2][r2][a+SEQ2];
+ t_r=CL->residue_index[s2][r2][a+R2];
+ norm2++;
if (hasch[t_s][t_r])
{
if (hasch[t_s][t_r]==FORBIDEN)
{
- score+=CL->residue_index[s2][r2][a+2];
- max_score+=CL->residue_index[s2][r2][a+2];
+ score+=((float)CL->residue_index[s2][r2][a+field]/NORM_F);
}
else
{
- double delta;
- delta=MIN(hasch[t_s][t_r],CL->residue_index[s2][r2][a+2]);
-
+ delta=MIN((((float)hasch[t_s][t_r]/NORM_F)),(((float)CL->residue_index[s2][r2][a+field]/NORM_F)));
score+=delta;
- max_score-=hasch[t_s][t_r];
- max_score+=delta;
- max_val=MAX(max_val,delta);
}
}
- else
- {
- max_score+=CL->residue_index[s2][r2][a+2];
- }
}
- max_score-=hasch[s2][r2];
- clean_residue_pair_hasch ( s1, r1,s2, r2, hasch, CL);
+ clean_residue_pair_hasch ( s1, r1,s2, r2, hasch, CL);
+
+ //New Normalization
+ norm1=MIN(norm1,norm2);
+
+ //Old Normailzation: on the number of sequences, useless when not doiang an all against all
+ //norm1=(CL->S)->nseq;
+
+ score=(norm1)?score/norm1:0;
+
+ return score*NORM_F;
+ }
+int residue_pair_extended_list_pc_old ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
+ {
+ //old norm does not take into account the number of effective intermediate sequences
+
+ double score=0;
+ int a, t_s, t_r;
+ static int **hasch;
+
+ static int max_len;
+ int field;
+ double delta;
+
+ /*
+
+ function documentation: start
+
+ int residue_pair_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2);
+
+ Computes the extended score for aligning residue seq1(r1) Vs seq2(r2)
+ Computes: matrix_score
+ non extended score
+ extended score
- if ( max_score==0)score=0;
- else if ( CL->normalise)
+ The extended score depends on the function index_res_constraint_list.
+ This function can compare a sequence with itself.
+
+ Associated functions: See util constraint list, list extention functions.
+
+ function documentation: end
+ */
+
+ field=CL->weight_field;
+ field=WE;
+
+ if ( r1<=0 || r2<=0)return 0;
+ if ( !hasch || max_len!=(CL->S)->max_len)
+ {
+ max_len=(CL->S)->max_len;
+ if ( hasch) free_int ( hasch, -1);
+ hasch=declare_int ( (CL->S)->nseq, (CL->S)->max_len+1);
+ }
+
+
+ /* Check matches for R1 in the indexed lib*/
+ hasch[s1][r1]=FORBIDEN;
+ for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
{
- score=((score*CL->normalise)/max_score)*SCORE_K;
- if (max_val> CL->normalise)
+ t_s=CL->residue_index[s1][r1][a+SEQ2];
+ t_r=CL->residue_index[s1][r1][a+R2];
+ hasch[t_s][t_r]=CL->residue_index[s1][r1][a+field];
+ }
+
+
+ /*Check Matches for r1 <-> r2 in the indexed lib */
+ for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
+ {
+ t_s=CL->residue_index[s2][r2][a+SEQ2];
+ t_r=CL->residue_index[s2][r2][a+R2];
+
+
+ if (hasch[t_s][t_r])
{
- score*=max_val/(double)CL->normalise;
+ if (hasch[t_s][t_r]==FORBIDEN)
+ {
+ score+=((float)CL->residue_index[s2][r2][a+2]/NORM_F);
+ }
+ else
+ {
+ delta=MIN((((float)hasch[t_s][t_r]/NORM_F)),(((float)CL->residue_index[s2][r2][a+field]/NORM_F)));
+ score+=delta;
+ }
}
}
- return (int) score;
+
+ clean_residue_pair_hasch ( s1, r1,s2, r2, hasch, CL);
+ score/=(CL->S)->nseq;
+ return score*NORM_F;
}
-int ** clean_residue_pair_hasch (int s1, int r1, int s2, int r2,int **hasch, Constraint_list *CL)
- {
- int a, t_s, t_r;
- if ( !hasch) return hasch;
-
- for (a=1; a< CL->residue_index[s1][r1][0]; a+=3)
- {
- t_s=CL->residue_index[s1][r1][a];
- t_r=CL->residue_index[s1][r1][a+1];
- hasch[t_s][t_r]=0;
- }
- hasch[s1][r1]=hasch[s2][r2]=0;
- return hasch;
- }
-int residue_pair_extended_list_old ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
+
+
+int residue_pair_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
{
double score=0;
-
+ double max_score=0;
+ double max_val=0;
+
int a, t_s, t_r;
static int **hasch;
static int max_len;
-
int field;
+
/*
+ new function: self normalized
function documentation: start
int residue_pair_extended_list ( Constraint_list *CL, int s1, int r1, int s2, int r2);
function documentation: end
*/
-
-
-
+
field=CL->weight_field;
-
+ field=WE;
+
if ( r1<=0 || r2<=0)return 0;
if ( !hasch || max_len!=(CL->S)->max_len)
{
hasch=declare_int ( (CL->S)->nseq, (CL->S)->max_len+1);
}
- CL=index_res_constraint_list ( CL, field);
-
- hasch[s1][r1]=100000;
- for (a=1; a< CL->residue_index[s1][r1][0]; a+=3)
+
+
+ /* Check matches for R1 in the indexed lib*/
+ hasch[s1][r1]=FORBIDEN;
+ for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s1][r1][a];
- t_r=CL->residue_index[s1][r1][a+1];
- hasch[t_s][t_r]=CL->residue_index[s1][r1][a+2];
+ t_s=CL->residue_index[s1][r1][a+SEQ2];
+ t_r=CL->residue_index[s1][r1][a+R2];
+ hasch[t_s][t_r]=CL->residue_index[s1][r1][a+field];
+ max_score+=CL->residue_index[s1][r1][a+field];
}
-
-
- for (a=1; a< CL->residue_index[s2][r2][0]; a+=3)
+
+ /*Check Matches for r1 <-> r2 in the indexed lib */
+ for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s2][r2][a];
- t_r=CL->residue_index[s2][r2][a+1];
+ t_s=CL->residue_index[s2][r2][a+SEQ2];
+ t_r=CL->residue_index[s2][r2][a+R2];
+
+
if (hasch[t_s][t_r])
{
- if (field==WE)
+ if (hasch[t_s][t_r]==FORBIDEN)
+ {
+ score+=CL->residue_index[s2][r2][a+field];
+ max_score+=CL->residue_index[s2][r2][a+field];
+ }
+ else
{
+ double delta;
+ delta=MIN(hasch[t_s][t_r],CL->residue_index[s2][r2][a+field]);
- score+=MIN(hasch[t_s][t_r],CL->residue_index[s2][r2][a+2]);
+ score+=delta;
+ max_score-=hasch[t_s][t_r];
+ max_score+=delta;
+ max_val=MAX(max_val,delta);
}
}
+ else
+ {
+ max_score+=CL->residue_index[s2][r2][a+field];
+ }
}
-
- score=(CL->normalise)?((score*CL->normalise)/CL->max_ext_value):score;
-
-
- for (a=1; a< CL->residue_index[s1][r1][0]; a+=3)
+ max_score-=hasch[s2][r2];
+ clean_residue_pair_hasch ( s1, r1,s2, r2, hasch, CL);
+
+
+ if ( max_score==0)score=0;
+ else if ( CL->normalise)
{
- t_s=CL->residue_index[s1][r1][a];
- t_r=CL->residue_index[s1][r1][a+1];
- hasch[t_s][t_r]=0;
+ score=((score*CL->normalise)/max_score)*SCORE_K;
+ if (max_val> CL->normalise)
+ {
+ score*=max_val/(double)CL->normalise;
+ }
}
- hasch[s1][r1]=hasch[s2][r2]=0;
-
-
- return (int)(score*SCORE_K);
+ return (int) score;
}
+int ** clean_residue_pair_hasch (int s1, int r1, int s2, int r2,int **hasch, Constraint_list *CL)
+ {
+ int a, t_s, t_r;
+ if ( !hasch) return hasch;
+
+ for (a=1; 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];
+ hasch[t_s][t_r]=0;
+ }
+ hasch[s1][r1]=hasch[s2][r2]=0;
+ return hasch;
+ }
+
int residue_pair_test_function ( Constraint_list *CL, int s1, int r1, int s2, int r2 )
{
double score=0;
hasch=declare_int ( (CL->S)->nseq, (CL->S)->max_len+1);
}
- CL=index_res_constraint_list ( CL, field);
+
hasch[s1][r1]=1000;
- for (a=1; a< CL->residue_index[s1][r1][0]; a+=3)
+ for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s1][r1][a];
- t_r=CL->residue_index[s1][r1][a+1];
- hasch[t_s][t_r]=CL->residue_index[s1][r1][a+2];
+ t_s=CL->residue_index[s1][r1][a+SEQ2];
+ t_r=CL->residue_index[s1][r1][a+R2];
+ hasch[t_s][t_r]=CL->residue_index[s1][r1][a+field];
}
- for (a=1; a< CL->residue_index[s2][r2][0]; a+=3)
+ for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s2][r2][a];
- t_r=CL->residue_index[s2][r2][a+1];
+ t_s=CL->residue_index[s2][r2][a+SEQ2];
+ t_r=CL->residue_index[s2][r2][a+R2];
if (hasch[t_s][t_r])
{
cons1=hasch[t_s][t_r];
- cons2=CL->residue_index[s2][r2][a+2];
+ cons2=CL->residue_index[s2][r2][a+field];
score +=MIN(cons1,cons2);
}
}
score=(CL->normalise)?((score*CL->normalise)/CL->max_ext_value):score;
- for (a=1; a< CL->residue_index[s1][r1][0]; a+=3)
+ for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s1][r1][a];
- t_r=CL->residue_index[s1][r1][a+1];
+ t_s=CL->residue_index[s1][r1][a+SEQ2];
+ t_r=CL->residue_index[s1][r1][a+R2];
hasch[t_s][t_r]=0;
}
hasch[s1][r1]=hasch[s2][r2]=0;
-
+
return (int)(score*SCORE_K);
}
field=CL->weight_field;
-
+ field=WE;
+
if ( r1<=0 || r2<=0)return 0;
if ( !hasch || max_len!=(CL->S)->max_len)
{
hasch=declare_int ( (CL->S)->nseq, (CL->S)->max_len+1);
}
- CL=index_res_constraint_list ( CL, field);
+
hasch[s1][r1]=100000;
- for (a=1; a< CL->residue_index[s1][r1][0]; a+=3)
+ for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s1][r1][a];
- t_r=CL->residue_index[s1][r1][a+1];
- hasch[t_s][t_r]=CL->residue_index[s1][r1][a+2];
- total_score+=CL->residue_index[s1][r1][a+2];
+ t_s=CL->residue_index[s1][r1][a+SEQ2];
+ t_r=CL->residue_index[s1][r1][a+R2];
+ hasch[t_s][t_r]=CL->residue_index[s1][r1][a+field];
+ total_score+=CL->residue_index[s1][r1][a+field];
}
- for (a=1; a< CL->residue_index[s2][r2][0]; a+=3)
+ for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s2][r2][a];
- t_r=CL->residue_index[s2][r2][a+1];
- total_score+=CL->residue_index[s1][r1][a+2];
+ t_s=CL->residue_index[s2][r2][a+SEQ2];
+ t_r=CL->residue_index[s2][r2][a+R2];
+ total_score+=CL->residue_index[s1][r1][a+field];
if (hasch[t_s][t_r])
{
- if (field==WE){score+=2*MIN(hasch[t_s][t_r],CL->residue_index[s2][r2][a+2]);}
+ if (field==WE){score+=2*MIN(hasch[t_s][t_r],CL->residue_index[s2][r2][a+field]);}
}
}
score=((CL->normalise*score)/total_score);
- for (a=1; a< CL->residue_index[s1][r1][0]; a+=3)
+ for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s1][r1][a];
- t_r=CL->residue_index[s1][r1][a+1];
+ t_s=CL->residue_index[s1][r1][a+SEQ2];
+ t_r=CL->residue_index[s1][r1][a+R2];
hasch[t_s][t_r]=0;
}
hasch[s1][r1]=hasch[s2][r2]=0;
for ( a=0; a< (CL->S)->nseq; a++)hasch[a]=vcalloc ( (CL->S)->len[a]+1, sizeof (int));
}
- CL=index_res_constraint_list ( CL, field);
+
hasch[s1][r1]=100000;
- for (a=1; a< CL->residue_index[s1][r1][0]; a+=3)
+ for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s1][r1][a];
- t_r=CL->residue_index[s1][r1][a+1];
- t_w=CL->residue_index[s1][r1][a+2];
+ 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+field];
if ( CL->seq_for_quadruplet[t_s])
{
- for ( b=1; b<CL->residue_index[t_s][t_r][0]; b+=3)
+ for ( b=1; b<CL->residue_index[t_s][t_r][0]; b+=ICHUNK)
{
- q_s=CL->residue_index[t_s][t_r][b];
- q_r=CL->residue_index[t_s][t_r][b+1];
+ q_s=CL->residue_index[t_s][t_r][b+SEQ2];
+ q_r=CL->residue_index[t_s][t_r][b+R2];
if (CL-> seq_for_quadruplet[q_s])
{
}
- for (s=0,score=0,a=1; a< CL->residue_index[s2][r2][0]; a+=3)
+ for (s=0,score=0,a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s2][r2][a];
- t_r=CL->residue_index[s2][r2][a+1];
- t_w=CL->residue_index[s2][r2][a+2];
+ t_s=CL->residue_index[s2][r2][a+SEQ2];
+ t_r=CL->residue_index[s2][r2][a+R2];
+ t_w=CL->residue_index[s2][r2][a+field];
if ( CL->seq_for_quadruplet[t_s])
{
- for ( b=1; b<CL->residue_index[t_s][t_r][0]; b+=3)
+ for ( b=1; b<CL->residue_index[t_s][t_r][0]; b+=ICHUNK)
{
- q_s=CL->residue_index[t_s][t_r][b];
- q_r=CL->residue_index[t_s][t_r][b+1];
- q_w=CL->residue_index[t_s][t_r][b+2];
+ q_s=CL->residue_index[t_s][t_r][b+SEQ2];
+ q_r=CL->residue_index[t_s][t_r][b+R2];
+ q_w=CL->residue_index[t_s][t_r][b+field];
if (hasch[q_s][q_r] && CL->seq_for_quadruplet[q_s])
s=MIN(hasch[q_s][q_r],MIN(CL->residue_index[t_s][t_r][b+2],q_w));
score=MAX(score, s);
}
}
- for (a=1; a< CL->residue_index[s1][r1][0]; a+=3)
+ for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s1][r1][a];
- t_r=CL->residue_index[s1][r1][a+1];
- t_w=CL->residue_index[s1][r1][a+2];
+ 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+field];
if ( CL->seq_for_quadruplet[t_s])
{
- for ( b=1; b<CL->residue_index[t_s][t_r][0]; b+=3)
+ for ( b=1; b<CL->residue_index[t_s][t_r][0]; b+=ICHUNK)
{
- q_s=CL->residue_index[t_s][t_r][b];
- q_r=CL->residue_index[t_s][t_r][b+1];
+ q_s=CL->residue_index[t_s][t_r][b+SEQ2];
+ q_r=CL->residue_index[t_s][t_r][b+R2];
hasch[q_s][q_r]=0;
}
}
for ( a=0; a< (CL->S)->nseq; a++)hasch[a]=vcalloc ( (CL->S)->len[a]+1, sizeof (int));
}
- CL=index_res_constraint_list ( CL, field);
+
hasch[s1][r1]=100000;
- for (a=1; a< CL->residue_index[s1][r1][0]; a+=3)
+ for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s1][r1][a];
- t_r=CL->residue_index[s1][r1][a+1];
- hasch[t_s][t_r]=CL->residue_index[s1][r1][a+2];
+ t_s=CL->residue_index[s1][r1][a+SEQ2];
+ t_r=CL->residue_index[s1][r1][a+R2];
+ hasch[t_s][t_r]=CL->residue_index[s1][r1][a+field];
}
- for (s=0, score=0,a=1; a< CL->residue_index[s2][r2][0]; a+=3)
+ for (s=0, score=0,a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s2][r2][a];
- t_r=CL->residue_index[s2][r2][a+1];
+ t_s=CL->residue_index[s2][r2][a+SEQ2];
+ t_r=CL->residue_index[s2][r2][a+R2];
if (hasch[t_s][t_r])
{
if (field==WE)
- {s=MIN(hasch[t_s][t_r],CL->residue_index[s2][r2][a+2]);
+ {s=MIN(hasch[t_s][t_r],CL->residue_index[s2][r2][a+WE]);
score=MAX(s,score);
}
}
}
- for (a=1; a< CL->residue_index[s1][r1][0]; a+=3)
+ for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s1][r1][a];
- t_r=CL->residue_index[s1][r1][a+1];
+ t_s=CL->residue_index[s1][r1][a+SEQ2];
+ t_r=CL->residue_index[s1][r1][a+R2];
hasch[t_s][t_r]=0;
}
hasch[s1][r1]=hasch[s2][r2]=0;
field=CL->weight_field;
if ( r1<=0 || r2<=0)return 0;
- else if ( !CL->L && CL->M)
+ else if ( !CL->residue_index && CL->M)
{
return evaluate_matrix_score (CL, s1,r1, s2, r2);
}
for ( a=0; a< (CL->S)->nseq; a++)hasch[a]=vcalloc ( (CL->S)->len[a]+1, sizeof (int));
}
- CL=index_res_constraint_list ( CL, field);
+
hasch[s1][r1]=100000;
- for (a=1; a< CL->residue_index[s1][r1][0]; a+=3)
+ for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s1][r1][a];
- t_r=CL->residue_index[s1][r1][a+1];
- hasch[t_s][t_r]=CL->residue_index[s1][r1][a+2];
+ t_s=CL->residue_index[s1][r1][a+SEQ2];
+ t_r=CL->residue_index[s1][r1][a+R2];
+ hasch[t_s][t_r]=CL->residue_index[s1][r1][a+WE];
}
- for (a=1; a< CL->residue_index[s2][r2][0]; a+=3)
+ for (a=1; a< CL->residue_index[s2][r2][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s2][r2][a];
- t_r=CL->residue_index[s2][r2][a+1];
+ t_s=CL->residue_index[s2][r2][a+SEQ2];
+ t_r=CL->residue_index[s2][r2][a+R2];
if (hasch[t_s][t_r])
{
- if (field==WE)score+=MIN(hasch[t_s][t_r],CL->residue_index[s2][r2][a+2] );
+ if (field==WE)score+=MIN(hasch[t_s][t_r],CL->residue_index[s2][r2][a+WE] );
}
}
score=(CL->normalise)?((score*CL->normalise)/CL->max_ext_value):score;
- for (a=1; a< CL->residue_index[s1][r1][0]; a+=3)
+ for (a=1; a< CL->residue_index[s1][r1][0]; a+=ICHUNK)
{
- t_s=CL->residue_index[s1][r1][a];
- t_r=CL->residue_index[s1][r1][a+1];
+ t_s=CL->residue_index[s1][r1][a+SEQ2];
+ t_r=CL->residue_index[s1][r1][a+R2];
hasch[t_s][t_r]=0;
}
hasch[s1][r1]=hasch[s2][r2]=0;
if (A==NULL)return 0;
- if (MODE!=2 || MODE==0 || (!CL->L && CL->M) || (!CL->L && CL->T)|| ns1==1 || ns2==1)
+ if (MODE!=2 || MODE==0 || (!CL->residue_index && CL->M) || (!CL->residue_index && CL->T)|| ns1==1 || ns2==1)
score=slow_get_dp_cost ( A, pos1, ns1, list1,col1, pos2, ns2, list2, col2, CL);
else if (MODE==1 || MODE==2)
score=fast_get_dp_cost ( A, pos1, ns1, list1,col1, pos2, ns2, list2, col2, CL);
col2=buf_col;
}
- CL=index_res_constraint_list ( CL, WE);
+
if ( !hasch1)
{
else
{
is_in_col1[s1][r1]=1;
- for (b=1; b< CL->residue_index[s1][r1][0]; b+=3)
+ for (b=1; b< CL->residue_index[s1][r1][0]; b+=ICHUNK)
{
- t_s=CL->residue_index[s1][r1][b];
- t_r=CL->residue_index[s1][r1][b+1];
- t_w=CL->residue_index[s1][r1][b+2];
- for ( c=1; c< CL->residue_index[t_s][t_r][0]; c+=3)
+ t_s=CL->residue_index[s1][r1][b+SEQ2];
+ t_r=CL->residue_index[s1][r1][b+R2];
+ t_w=CL->residue_index[s1][r1][b+WE];
+ for ( c=1; c< CL->residue_index[t_s][t_r][0]; c+=ICHUNK)
{
- q_s=CL->residue_index[t_s][t_r][c];
- q_r=CL->residue_index[t_s][t_r][c+1];
- q_w=CL->residue_index[t_s][t_r][c+2];
+ q_s=CL->residue_index[t_s][t_r][c+SEQ2];
+ q_r=CL->residue_index[t_s][t_r][c+R2];
+ q_w=CL->residue_index[t_s][t_r][c+WE];
hasch1[q_s][q_r]+=MIN(q_w, t_w);
n_hasch1[q_s][q_r]++;
}
else
{
is_in_col2[s2][r2]=1;
- for (b=1; b< CL->residue_index[s2][r2][0]; b+=3)
+ for (b=1; b< CL->residue_index[s2][r2][0]; b+=ICHUNK)
{
- t_s=CL->residue_index[s2][r2][b];
- t_r=CL->residue_index[s2][r2][b+1];
- t_w=CL->residue_index[s2][r2][b+2];
- for ( c=1; c< CL->residue_index[t_s][t_r][0]; c+=3)
+ t_s=CL->residue_index[s2][r2][b+SEQ2];
+ t_r=CL->residue_index[s2][r2][b+R2];
+ t_w=CL->residue_index[s2][r2][b+WE];
+ for ( c=1; c< CL->residue_index[t_s][t_r][0]; c+=ICHUNK)
{
- q_s=CL->residue_index[t_s][t_r][c];
- q_r=CL->residue_index[t_s][t_r][c+1];
- q_w=CL->residue_index[t_s][t_r][c+2];
+ q_s=CL->residue_index[t_s][t_r][c+SEQ2];
+ q_r=CL->residue_index[t_s][t_r][c+R2];
+ q_w=CL->residue_index[t_s][t_r][c+WE];
hasch2[q_s][q_r]+=MIN(t_w, q_w);
n_hasch2[q_s][q_r]++;
}
if (r2<0);
else
{
- for (b=1; b< CL->residue_index[s2][r2][0]; b+=3)
+ for (b=1; b< CL->residue_index[s2][r2][0]; b+=ICHUNK)
{
- t_s=CL->residue_index[s2][r2][b];
- t_r=CL->residue_index[s2][r2][b+1];
+ t_s=CL->residue_index[s2][r2][b+SEQ2];
+ t_r=CL->residue_index[s2][r2][b+R2];
- for ( c=1; c< CL->residue_index[t_s][t_r][0]; c+=3)
+ for ( c=1; c< CL->residue_index[t_s][t_r][0]; c+=ICHUNK)
{
- q_s=CL->residue_index[t_s][t_r][c];
- q_r=CL->residue_index[t_s][t_r][c+1];
+ q_s=CL->residue_index[t_s][t_r][c+SEQ2];
+ q_r=CL->residue_index[t_s][t_r][c+R2];
if ( hasch2[q_s][q_r] && hasch1[q_s][q_r]&& !(is_in_col1[q_s][q_r] || is_in_col2[q_s][q_r]))
{
score+=MIN(hasch2[q_s][q_r]*(n_hasch1[q_s][q_r]),hasch1[q_s][q_r]*(n_hasch2[q_s][q_r]));
{
is_in_col1[s1][r1]=0;
hasch1[s1][r1]=0;
- for (b=1; b< CL->residue_index[s1][r1][0]; b+=3)
+ for (b=1; b< CL->residue_index[s1][r1][0]; b+=ICHUNK)
{
- t_s=CL->residue_index[s1][r1][b];
- t_r=CL->residue_index[s1][r1][b+1];
- for ( c=1; c< CL->residue_index[t_s][t_r][0]; c+=3)
+ t_s=CL->residue_index[s1][r1][b+SEQ2];
+ t_r=CL->residue_index[s1][r1][b+R2];
+ for ( c=1; c< CL->residue_index[t_s][t_r][0]; c+=ICHUNK)
{
- q_s=CL->residue_index[t_s][t_r][c];
- q_r=CL->residue_index[t_s][t_r][c+1];
+ q_s=CL->residue_index[t_s][t_r][c+SEQ2];
+ q_r=CL->residue_index[t_s][t_r][c+R2];
hasch1[q_s][q_r]=0;
n_hasch1[q_s][q_r]=0;
}
- CL=index_res_constraint_list ( CL, WE);
+
if ( !hasch1)
{
else
{
is_in_col1[s1][r1]=1;
- for (b=1; b< CL->residue_index[s1][r1][0]; b+=3)
+ for (b=1; b< CL->residue_index[s1][r1][0]; b+=ICHUNK)
{
- t_s=CL->residue_index[s1][r1][b];
- t_r=CL->residue_index[s1][r1][b+1];
- hasch1[t_s][t_r]+=CL->residue_index[s1][r1][b+2];
+ t_s=CL->residue_index[s1][r1][b+SEQ2];
+ t_r=CL->residue_index[s1][r1][b+R2];
+ hasch1[t_s][t_r]+=CL->residue_index[s1][r1][b+WE];
n_hasch1[t_s][t_r]++;
}
}
else
{
is_in_col2[s2][r2]=1;
- for (b=1; b< CL->residue_index[s2][r2][0]; b+=3)
+ for (b=1; b< CL->residue_index[s2][r2][0]; b+=ICHUNK)
{
- t_s=CL->residue_index[s2][r2][b];
- t_r=CL->residue_index[s2][r2][b+1];
+ t_s=CL->residue_index[s2][r2][b+SEQ2];
+ t_r=CL->residue_index[s2][r2][b+R2];
- hasch2[t_s][t_r]+=CL->residue_index[s2][r2][b+2];
+ hasch2[t_s][t_r]+=CL->residue_index[s2][r2][b+WE];
n_hasch2[t_s][t_r]++;
}
}
if (r2<0);
else
{
- for (b=1; b< CL->residue_index[s2][r2][0]; b+=3)
+ for (b=1; b< CL->residue_index[s2][r2][0]; b+=ICHUNK)
{
- t_s=CL->residue_index[s2][r2][b];
- t_r=CL->residue_index[s2][r2][b+1];
+ t_s=CL->residue_index[s2][r2][b+SEQ2];
+ t_r=CL->residue_index[s2][r2][b+R2];
if ( hasch2[t_s][t_r] && hasch1[t_s][t_r]&& !(is_in_col1[t_s][t_r] || is_in_col2[t_s][t_r]))
{
{
is_in_col1[s1][r1]=0;
hasch1[s1][r1]=0;
- for (b=1; b< CL->residue_index[s1][r1][0]; b+=3)
+ for (b=1; b< CL->residue_index[s1][r1][0]; b+=ICHUNK)
{
- t_s=CL->residue_index[s1][r1][b];
- t_r=CL->residue_index[s1][r1][b+1];
+ t_s=CL->residue_index[s1][r1][b+SEQ2];
+ t_r=CL->residue_index[s1][r1][b+R2];
hasch1[t_s][t_r]=0;
n_hasch1[t_s][t_r]=0;
if (r1<0);
else
{
- for (b=1; b< CL->residue_index[s1][r1][0]; b+=3)
+ for (b=1; b< CL->residue_index[s1][r1][0]; b+=ICHUNK)
{
- t_s=CL->residue_index[s1][r1][b];
- t_r=CL->residue_index[s1][r1][b+1];
+ t_s=CL->residue_index[s1][r1][b+SEQ2];
+ t_r=CL->residue_index[s1][r1][b+R2];
if ( hasch1[t_s][t_r] && hasch2[t_s][t_r]&& !(is_in_col2[t_s][t_r] || is_in_col1[t_s][t_r]))
{
{
is_in_col2[s2][r2]=0;
hasch1[s2][r2]=0;
- for (b=1; b< CL->residue_index[s2][r2][0]; b+=3)
+ for (b=1; b< CL->residue_index[s2][r2][0]; b+=ICHUNK)
{
- t_s=CL->residue_index[s2][r2][b];
- t_r=CL->residue_index[s2][r2][b+1];
+ t_s=CL->residue_index[s2][r2][b+SEQ2];
+ t_r=CL->residue_index[s2][r2][b+R2];
hasch2[t_s][t_r]=0;
n_hasch2[t_s][t_r]=0;
return (int)score;
}
-int fast_get_dp_cost_3 ( Alignment *A, int**pos1, int ns1, int*list1, int col1, int**pos2, int ns2, int*list2, int col2, Constraint_list *CL)
-{
- static int last_tag;
- static Constraint_list *NCL;
- int score;
-
- if ( ns1==1 && ns2==1)
- {
- return slow_get_dp_cost( A,pos1, ns1,list1, col1, pos2, ns2, list2, col2,CL);
- }
-
- if ( last_tag !=A->random_tag)
- {
- int *ns, **ls;
-
- last_tag=A->random_tag;
- ns=vcalloc (2, sizeof (int));ns[0]=ns1; ns[1]=ns2;
- ls=vcalloc (2, sizeof (int*));ls[0]=list1; ls[1]=list2;
-
- NCL=progressive_index_res_constraint_list ( A, ns, ls, CL);
- vfree (ls); vfree (ns);
- }
- score=residue_pair_extended_list ( NCL,list1[0],col1, list2[0], col2);
- score=(CL->normalise)?((score*CL->normalise)/CL->max_ext_value):score;
- score=(score-SCORE_K*CL->nomatch);
- return score;
-}
}
Constraint_list * choose_extension_mode ( char *extend_mode, Constraint_list *CL)
{
+ //evaluation_functions: residues start at 1, sequences at 0;
+
if ( !CL)
{
fprintf ( stderr, "\nWarning: CL was not set");
CL->evaluate_residue_pair=residue_pair_extended_list;
CL->get_dp_cost =fast_get_dp_cost;
}
- else if ( strm ( extend_mode, "test_triplet") && !CL->M)
- {
- CL->evaluate_residue_pair=residue_pair_extended_list;
- CL->get_dp_cost =fast_get_dp_cost_3;
- }
+
else if ( strm ( extend_mode, "very_fast_triplet") && !CL->M)
{
CL->evaluate_residue_pair=residue_pair_extended_list;
fprintf ( stderr, "\nERROR: %s is an unknown extend_mode[FATAL:%s]\n", extend_mode, PROGRAM);
myexit (EXIT_FAILURE);
}
+
return CL;
}
int ** combine_two_matrices ( int **mat1, int **mat2)
return -1;
}
-/*********************************COPYRIGHT NOTICE**********************************/
+/******************************COPYRIGHT NOTICE*******************************/
/*© Centro de Regulacio Genomica */
/*and */
/*Cedric Notredame */
-/*Tue Oct 27 10:12:26 WEST 2009. */
+/*Fri Feb 18 08:27:45 CET 2011 - Revision 596. */
/*All rights reserved.*/
/*This file is part of T-COFFEE.*/
/**/
/**/
/**/
/* */
-/*********************************COPYRIGHT NOTICE**********************************/
+/******************************COPYRIGHT NOTICE*******************************/
git://source.jalview.org / jabaws.git / blobdiff