JPRED-2 Add sources of all binaries (except alscript) to Git
[jpred.git] / sources / multicoil / PairCoilDiffer_cuts.c
diff --git a/sources/multicoil/PairCoilDiffer_cuts.c b/sources/multicoil/PairCoilDiffer_cuts.c
new file mode 100644 (file)
index 0000000..c8f7994
--- /dev/null
@@ -0,0 +1,257 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+#include "scio.h"
+#include "switches.h"
+#include "likelihood.h"
+#include "sc.h"
+#include "scconst.h"
+#include "options.h"
+#include "compute_like.h"
+#include "interface.h"
+#include "stats.h"
+
+void average_score_over_coils2(Sequence seq, 
+                               double scores[MAXSEQLEN][POSNUM+1],
+                              double PreprocLike[MAXSEQLEN][POSNUM+1],
+                              int CoilDiffMethod, int offset_to_use,
+                              char offsets[MAXSEQLEN],
+                              int start_at_reg_shift)
+{
+  int i=0, j,k;
+  int start_coil[POSNUM+1];
+  double total_coil_score[POSNUM+1];
+  double avg_coil_score[POSNUM+1];
+  int regist;
+
+  double total_coil_length[POSNUM +1];  /* A weighted length. */
+
+  for (j=0; j<POSNUM+1; j++) {
+    start_coil[j] = -1;
+    total_coil_score[j]=0;
+    avg_coil_score[j]=0;
+    total_coil_length[j]=0;
+  }
+
+  if (offset_to_use == -1 ) offset_to_use = 7;  /* Adjust "all" offset */
+                                              /* to max offset locally. */
+
+  while (i<=seq.seqlen) {
+    for (k=0; k<=POSNUM; k++) {
+      if (k!=7) regist = (i+k)%POSNUM;
+      else regist =7;
+
+      /*  Consider < 1/10 percent to be zero */
+      if ( ( PreprocLike[i][regist] < .001) || (i==seq.seqlen) ||
+          ( (i>0) && (offsets[i] != offsets[i-1])  && start_at_reg_shift && 
+           (PreprocLike[i-1][7]>= .001) ) )  {       
+      /** If not in a coil (value of -HUGE_VAL signals not in coil).  */
+        if (start_coil[k] != -1)  {  /* Ended coil at i-1. */
+          avg_coil_score[k] = total_coil_score[k]/total_coil_length[k];
+          for  (j=start_coil[k]; j<i; j++) {
+            if (k!=7)
+              scores[j][(j+k)%POSNUM]= avg_coil_score[k];
+            else
+              scores[j][7]= avg_coil_score[k];
+
+          }
+        }
+        total_coil_score[k] =0;
+void average_score_over_coils2(Sequence seq, 
+                               double scores[MAXSEQLEN][POSNUM+1],
+                              double PreprocLike[MAXSEQLEN][POSNUM+1],
+                              int CoilDiffMethod, int offset_to_use,
+                              char offsets[MAXSEQLEN],
+                              int start_at_reg_shift)
+{
+  int i=0, j,k;
+  int start_coil[POSNUM+1];
+  double total_coil_score[POSNUM+1];
+  double avg_coil_score[POSNUM+1];
+  int regist;
+
+  double total_coil_length[POSNUM +1];  /* A weighted length. */
+
+  for (j=0; j<POSNUM+1; j++) {
+    start_coil[j] = -1;
+    total_coil_score[j]=0;
+    avg_coil_score[j]=0;
+    total_coil_length[j]=0;
+  }
+
+  if (offset_to_use == -1 ) offset_to_use = 7;  /* Adjust "all" offset */
+                                              /* to max offset locally. */
+
+  while (i<=seq.seqlen) {
+    for (k=0; k<=POSNUM; k++) {
+      if (k!=7) regist = (i+k)%POSNUM;
+      else regist =7;
+
+      /*  Consider < 1/10 percent to be zero */
+      if ( ( PreprocLike[i][regist] < .001) || (i==seq.seqlen) ||
+          ( (i>0) && (offsets[i] != offsets[i-1])  && start_at_reg_shift && 
+           (PreprocLike[i-1][7]>= .001) ) )  {       
+      /** If not in a coil (value of -HUGE_VAL signals not in coil).  */
+        if (start_coil[k] != -1)  {  /* Ended coil at i-1. */
+          avg_coil_score[k] = total_coil_score[k]/total_coil_length[k];
+          for  (j=start_coil[k]; j<i; j++) {
+            if (k!=7)
+              scores[j][(j+k)%POSNUM]= avg_coil_score[k];
+            else
+              scores[j][7]= avg_coil_score[k];
+
+          }
+        }
+        total_coil_score[k] =0;
+        if (i != seq.seqlen)   /* Consider < 1/10 percent to be 0 */
+          if  ( PreprocLike[i][regist] < .001)  {  /* Not in coil anymore. */
+            scores[i][regist]=0;
+            start_coil[k] = -1;
+            total_coil_length[k] =0;
+          }
+          else {  /* Changing registers. */
+            start_coil[k]= i;
+            if (CoilDiffMethod==0) {
+              total_coil_length[k] = 1;
+              total_coil_score[k] = scores[i][regist];
+            }
+            else if (CoilDiffMethod ==1) {
+              total_coil_length[k]= PreprocLike[i][regist];
+              total_coil_score[k] = scores[i][regist]*
+                PreprocLike[i][regist];
+            }
+          }
+      }   /** Ends if ending a coil.  **/
+
+      else {    /* If in a coil.  */
+        if (CoilDiffMethod ==0) {   /* Average scores over coil residues. */
+          total_coil_score[k] += scores[i][regist];
+          
+          total_coil_length[k] += 1;
+        }
+        else if (CoilDiffMethod == 1) {  /* Weighted average scores.  */
+          total_coil_score[k] += scores[i][regist]*
+                 PreprocLike[i][regist];
+          total_coil_length[k] += PreprocLike[i][regist];
+        }
+        if (start_coil[k] == -1) 
+          start_coil[k] =i;
+      }
+    }
+
+    i++;
+  }
+
+  /*  Need the following hack to get the correct register in max off method.*/
+  if ( (offset_to_use == 7) && (!start_at_reg_shift)) 
+    for (i=0; i<seq.seqlen; i++) {
+      for (j=0; j<POSNUM; j++)
+        scores[i][j]= -HUGE_VAL;
+      scores[i][(i+offsets[i]) %POSNUM] = scores[i][7];
+    }
+}
+
+
+
+/*******************************************************************/
+
+void PairCoilDiffer(int mode, Sequence sequence, char lib[MAXFUNCTNUM], 
+                     double *maxscore,char offsets[MAXSEQLEN], 
+                     double preproc_like[MAXSEQLEN][POSNUM+1], 
+                     int offset_to_use, int table, 
+                     double scores[MAXSEQLEN][POSNUM+1],
+                     int Avg_Coil_Score, int CoilDiffMethod,
+                     int start_coil_at_reg_shift)
+
+{  
+  int i,j,k;
+  double sc2scores[MAXSEQLEN][POSNUM];
+  char diff_offsets[MAXSEQLEN];
+  double total_coilscore[POSNUM+1];
+  int coil_length[POSNUM +1];
+  int like2or3;
+  
+  extern long pfreqs[MAX_TABLE_NUMBER][AANUM][POSNUM]; 
+  extern long ptotals[MAX_TABLE_NUMBER][POSNUM];
+  extern long pfreqp[MAX_TABLE_NUMBER][AANUM][AANUM][POSNUM][POSNUM]; 
+  extern long ptotalp[MAX_TABLE_NUMBER][POSNUM][POSNUM];
+
+
+  int ProlineFreeWindow= mode & PROLINE_FREE_MODE;
+  
+  if (table ==0) like2or3= 2;
+  else like2or3=3;
+
+  for (i=0; i<=POSNUM; i++) {
+    coil_length[i] = 0;
+    total_coilscore[i] =0;
+  }
+
+
+
+  /* score the sequence using sc2seq */
+  scseqadj(lib,
+          sequence.seq, sequence.seqlen,
+          sc2scores,
+          diff_offsets,
+          maxscore, mode, 1);  /* Get differentiator score.*/
+
+
+
+
+
+/*  *maxscore= exp(*maxscore);   */
+
+  *maxscore =0;
+
+  for (i=0; i<sequence.seqlen; ++i) {
+    /*  Do the most probabable register if offset_to_use is -1 or 7. */
+    if ((offset_to_use == -1) || (offset_to_use ==7)) 
+      if (like2or3 ==2)  
+       scores[i][7] = 1 - exp(sc2scores[i][offsets[i]]); 
+      else 
+       scores[i][7] =  exp(sc2scores[i][offsets[i]]); 
+     
+    else    /*  Offset to use is 0-6. */
+      if (like2or3 ==  2) 
+       scores[i][7] = 1 - exp(sc2scores[i][offset_to_use]);
+      else
+       scores[i][7] = exp(sc2scores[i][offset_to_use]);
+                /*****Trimer like is 1 - dimerlike ***/
+
+    if (preproc_like[i][7] <.001)   /* Consider < 1/10 percent to be 0 */
+      scores[i][7]=0;  /* Don't score non-coiled regions. */
+
+
+
+    for (j=0; j<7; ++j)  { /** sc2scores[i][j] is the score for position
+                            i+j%POSNUM  **/
+      if (preproc_like[i][(i+j)%POSNUM]> 0) { 
+       if (like2or3==2) 
+         scores[i][(i+j)%POSNUM] = 1 - exp(sc2scores[i][j]);
+       else 
+         scores[i][(i+j)%POSNUM] = exp(sc2scores[i][j]);
+      }
+      else 
+       scores[i][(i+j)%POSNUM]= 0;
+    }
+    
+  }
+       
+
+  if (Avg_Coil_Score)
+    average_score_over_coils2(sequence, scores, preproc_like,
+                             CoilDiffMethod, offset_to_use,offsets,
+                             start_coil_at_reg_shift);
+
+
+
+  for (i=0; i<sequence.seqlen; i++)
+    if (scores[i][7] > *maxscore) *maxscore = scores[i][7];
+
+
+}
+
+
+