JPRED-2 Add sources of all binaries (except alscript) to Git
[jpred.git] / sources / multicoil / PairCoilDiffer.c
diff --git a/sources/multicoil/PairCoilDiffer.c b/sources/multicoil/PairCoilDiffer.c
new file mode 100644 (file)
index 0000000..f0661af
--- /dev/null
@@ -0,0 +1,179 @@
+#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 PairCoilDifferDimension(int mode, Sequence sequence, 
+                            char lib[MAXFUNCTNUM], double *maxscore,
+                            char preproc_offsets[MAXSEQLEN],
+                            double preproc_scores[MAXSEQLEN][POSNUM+1],
+                 double all_scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
+                 char all_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN], int dimension,
+                            int offset_to_use,
+                            int number_differ_score_dim, double bound)
+{
+  int i,j, pass,pos;
+  double sc2scores[MAXSEQLEN][POSNUM];
+  char diff_offsets[MAXSEQLEN];
+  extern int which_differentiator;
+  extern int differ_window_length;
+  int number_passes=1;
+  int start_coil =-1;
+
+  double coil_score[POSNUM];
+
+  /** which_differentiator = 1 means do window score starting at residue**/
+  /** If it is 2 means find both max and min window.  ******/
+  /** Passing 1 into scseqadj makes it return residue scores so the **/
+  /** residue gets the window score for the window starting at that **/
+  /** residue.  2 returns the max window score containing the residue **/
+  /** and 3 returns the min window score. **/
+
+  WINDOW = differ_window_length;  /** TEMPORARY.  **/
+
+  if (which_differentiator ==2) number_passes=2;
+  else number_passes=1;
+
+  for ( pass=0; pass<number_passes; pass++) {
+    if (which_differentiator ==3 || which_differentiator==4) { 
+                                /** get coiled regions from preproc. **/
+      i=0;
+      while(i<=sequence.seqlen) {
+       if ( (start_coil >= 0) && ( (i==sequence.seqlen) || 
+          (preproc_scores[i][7]<=bound) ||  /* End coil */
+          ( (i>0) && 
+           (preproc_offsets[i] != preproc_offsets[i-1])) )) { /*regshift*/
+/*       WINDOW= i-start_coil; */
+/*******
+         scseqadj(lib,
+                  &(sequence.seq[start_coil]), i-start_coil,
+                  &(sc2scores[start_coil][0]),
+                  &(diff_offsets[start_coil]),
+                  maxscore, mode, 1);  
+*******/
+/**          score_differentiator_window(lib,
+              &(sequence.seq[start_coil]), i-start_coil,
+              coil_score,mode, 1);
+**/
+          score_differentiator_window(lib,
+              &(sequence.seq[start_coil]), i-start_coil,
+              coil_score,mode, 1);
+
+/*      fprintf(stderr,"\n start_coil = %d, end_coil =%d",start_coil, i-1); */
+
+         for(j=start_coil; j<i; j++)
+           for(pos=0;pos<POSNUM;pos++)
+             if (which_differentiator==3) 
+               sc2scores[j][pos]= coil_score[pos];
+             else if (which_differentiator==4) 
+               sc2scores[j][pos] = coil_score[pos]/(i-start_coil);
+
+         start_coil = -1;
+         
+       }
+       
+       if ( (preproc_scores[i][7] > bound) && (start_coil ==-1))
+         start_coil =i;
+       i++;
+      }
+    }
+    else
+      scseqadj(lib,
+              sequence.seq, sequence.seqlen,
+              sc2scores,
+              diff_offsets,
+              maxscore, mode, pass + which_differentiator); 
+                               /* Get differentiator score.*/
+                              /* pass + which_diff = 2 is max_window_score **/
+                              /*                     3 is min_window_score **/
+
+    for (i=0; i<sequence.seqlen; ++i) {
+      /* Use the offset given by the preproc if offset_to_use is -1 or 7 */
+      if ((offset_to_use == -1) || (offset_to_use ==7))  {
+       all_scores[dimension+ pass*number_differ_score_dim/2][i][7] = 
+         sc2scores[i][preproc_offsets[i]];
+       all_offsets[dimension+ pass*number_differ_score_dim/2][i] = 
+         preproc_offsets[i];
+      }
+      else {    /*  Offset to use is 0-6. */
+       all_scores[dimension+ pass*number_differ_score_dim/2][i][7] = 
+         sc2scores[i][offset_to_use];
+       all_offsets[dimension+ pass*number_differ_score_dim/2][i]= 
+         offset_to_use;
+      }
+    
+
+      /* do other registers */
+      for (j=0; j<7; ++j)  { /** sc2scores[i][j] is the score for position
+                            i+j%POSNUM  **/
+       all_scores[dimension+ pass*number_differ_score_dim/2][i][(i+j)%POSNUM]=
+         sc2scores[i][j];
+      }
+    }
+  }
+  
+}
+
+
+
+
+/******************************************************************/
+void zero_out_non_coils( 
+                    double all_scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
+                    double preproc_like[MAXSEQLEN][POSNUM+1],
+                    int seqlen, double bound)
+{
+
+  int res, reg, olig;
+
+  for(olig=0; olig<2; olig++) 
+    for (res=0; res<seqlen; res++)
+      for (reg=0;reg<=POSNUM;reg++)
+       if (preproc_like[res][7] > bound)
+         all_scores[olig][res][reg] = all_scores[olig+3][res][reg];
+       else
+         all_scores[olig][res][reg] =0;
+}
+
+
+void zero_out_bad_windows(int differ_window_length, 
+                    double all_scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
+                    double preproc_like[MAXSEQLEN][POSNUM+1],
+                    int seqlen, double bound)
+{
+  int i,j;
+  int last_bad_window_start=-1;
+  int olig_state, pos;
+/*** For now just look at score of first and last residue in window to **/
+/**  see if they are above bound.  To really do it right should look at */
+/** min score in window (using min_over_windows on preproc_like), but   */
+/** the format expected by min_over_windows would require the POSNUM+1  */
+/** indice of preproc_like be in the first, not second position....     */
+
+  for (i=0; i<seqlen; i++)
+    if (preproc_like[i][7] < bound)  { /** Zero out windows. **/
+      for (j=MAX(last_bad_window_start+1,i-differ_window_length); j<=i; j++)
+       for(olig_state=0; olig_state<2; olig_state++)
+         for(pos=0; pos<=POSNUM; pos++) 
+           all_scores[olig_state][j][pos]=0;
+
+      last_bad_window_start=i;
+    }
+    else   /* Do this in case bound has changed....using stored value in +3 */
+      for(olig_state=0; olig_state<2; olig_state++)
+       for(pos=0; pos<=POSNUM; pos++)
+         all_scores[olig_state][i][pos]=all_scores[olig_state+3][i][pos];
+}
+
+
+
+