JWS-112 Bumping version of ClustalO (src, binaries and windows) to version 1.2.4.
[jabaws.git] / binaries / src / clustalo / src / hhalign / hhhit-C.h
index cd200f9..7f2296c 100644 (file)
@@ -15,7 +15,7 @@
  ********************************************************************/
 
 /*
- * RCS $Id: hhhit-C.h 245 2011-06-15 12:38:53Z fabian $
+ * RCS $Id: hhhit-C.h 316 2016-12-16 16:14:39Z fabian $
  */
 
 // hhhit.C
@@ -47,6 +47,7 @@ using std::ofstream;
 #include "hhalignment.h" // class Alignment
 #include "hhhitlist.h"   // class HitList
 #endif
+//#include "new_new.h" /* memory tracking */
 
 #define CALCULATE_MAX6(max, var1, var2, var3, var4, var5, var6, varb) \
 if (var1>var2) { max=var1; varb=STOP;} \
@@ -159,10 +160,10 @@ Hit::Delete()
       delete[] dbfile; dbfile = NULL;
       for (int k=0; k<n_display; k++) 
           {
-              //delete[] sname[k]; sname[k] = NULL;
+              delete[] sname[k]; sname[k] = NULL;  /* this seems to be the long lost piece of memory, FS, 2016-06-09 */
               delete[] seq[k]; seq[k] = NULL;
           }
-      //delete[] sname; sname = NULL;
+      delete[] sname; sname = NULL;  /* this seems to be the long lost piece of memory, FS, 2016-06-09 */
       delete[] seq; seq = NULL;
     }
 }
@@ -176,20 +177,20 @@ void
 Hit::AllocateBacktraceMatrix(int Nq, int Nt)
 {
   int i;
-  bMM=new(char*[Nq]);
-  bMI=new(char*[Nq]);
-  bIM=new(char*[Nq]);
-  bDG=new(char*[Nq]);
-  bGD=new(char*[Nq]);
-  cell_off=new(char*[Nq]);
+  bMM=new char*[Nq];
+  bMI=new char*[Nq];
+  bIM=new char*[Nq];
+  bDG=new char*[Nq];
+  bGD=new char*[Nq];
+  cell_off=new char*[Nq];
   for (i=0; i<Nq; i++) 
     {
-      bMM[i]=new(char[Nt]);
-      bMI[i]=new(char[Nt]);
-      bIM[i]=new(char[Nt]);
-      bGD[i]=new(char[Nt]);
-      bDG[i]=new(char[Nt]);
-      cell_off[i]=new(char[Nt]);
+      bMM[i]=new char[Nt];
+      bMI[i]=new char[Nt];
+      bIM[i]=new char[Nt];
+      bGD[i]=new char[Nt];
+      bDG[i]=new char[Nt];
+      cell_off[i]=new char[Nt] ;
       if (!bMM[i] || !bMI[i] || !bIM[i] || !bGD[i] || !bDG[i] || !cell_off[i]) 
        {
          fprintf(stderr,"Error: out of memory while allocating row %i (out of %i) for dynamic programming matrices \n",i+1,Nq);
@@ -209,21 +210,23 @@ void
 Hit::DeleteBacktraceMatrix(int Nq)
 {
   int i;
-  for (i=0; i<Nq; i++) 
-    {
-      delete[] bMM[i]; bMM[i] = NULL;
-      delete[] bMI[i]; bMI[i] = NULL;
-      delete[] bIM[i]; bIM[i] = NULL;
-      delete[] bGD[i]; bGD[i] = NULL;
-      delete[] bDG[i]; bDG[i] = NULL;
-      delete[] cell_off[i]; cell_off[i] = NULL;
-    }
-  delete[] bMM; bMM = NULL;
-  delete[] bMI; bMI = NULL;
-  delete[] bIM; bIM = NULL;
-  delete[] bDG; bDG = NULL;
-  delete[] bGD; bGD = NULL;
-  delete[] cell_off; cell_off = NULL;
+
+  if (NULL != bMM){ /* FS, r259 -> r260 */
+      for (i=0; i<Nq; i++)  {
+          delete[] bMM[i]; bMM[i] = NULL;
+          delete[] bMI[i]; bMI[i] = NULL;
+          delete[] bIM[i]; bIM[i] = NULL;
+          delete[] bGD[i]; bGD[i] = NULL;
+          delete[] bDG[i]; bDG[i] = NULL;
+          delete[] cell_off[i]; cell_off[i] = NULL;
+      }
+      delete[] bMM; bMM = NULL;
+      delete[] bMI; bMI = NULL;
+      delete[] bIM; bIM = NULL;
+      delete[] bDG; bDG = NULL;
+      delete[] bGD; bGD = NULL;
+      delete[] cell_off; cell_off = NULL;
+  }
 }
 
 
@@ -234,19 +237,19 @@ Hit::DeleteBacktraceMatrix(int Nq)
 void 
 Hit::AllocateForwardMatrix(int Nq, int Nt)
 {
-  F_MM=new(double*[Nq]);
-  F_MI=new(double*[Nq]);
-  F_DG=new(double*[Nq]);
-  F_IM=new(double*[Nq]);
-  F_GD=new(double*[Nq]);
-  scale=new(double[Nq+1]); // need Nq+3?
+  F_MM=new double*[Nq];
+  F_MI=new double*[Nq];
+  F_DG=new double*[Nq];
+  F_IM=new double*[Nq];
+  F_GD=new double*[Nq];
+  scale=new double[Nq+1] ; // need Nq+3?
   for (int i=0; i<Nq; i++) 
     {
-      F_MM[i] = new(double[Nt]);
-      F_MI[i] = new(double[Nt]);
-      F_DG[i] = new(double[Nt]);
-      F_IM[i] = new(double[Nt]);
-      F_GD[i] = new(double[Nt]);
+      F_MM[i] = new double[Nt];
+      F_MI[i] = new double[Nt];
+      F_DG[i] = new double[Nt];
+      F_IM[i] = new double[Nt];
+      F_GD[i] = new double[Nt];
       if (!F_MM[i] || !F_MI[i] || !F_IM[i] || !F_GD[i] || !F_DG[i]) 
        {
          fprintf(stderr,"Error: out of memory while allocating row %i (out of %i) for dynamic programming matrices \n",i+1,Nq);
@@ -266,20 +269,22 @@ Hit::AllocateForwardMatrix(int Nq, int Nt)
 void 
 Hit::DeleteForwardMatrix(int Nq)
 {
-  for (int i=0; i<Nq; i++) 
-    {
-      delete[] F_MM[i]; F_MM[i] = NULL;
-      delete[] F_MI[i]; F_MI[i] = NULL;
-      delete[] F_IM[i]; F_IM[i] = NULL;
-      delete[] F_GD[i]; F_GD[i] = NULL;
-      delete[] F_DG[i]; F_DG[i] = NULL;
+
+    if (NULL != F_MM){ /* FS, r259 -> r260 */
+        for (int i=0; i<Nq; i++) {
+            delete[] F_MM[i]; F_MM[i] = NULL;
+            delete[] F_MI[i]; F_MI[i] = NULL;
+            delete[] F_IM[i]; F_IM[i] = NULL;
+            delete[] F_GD[i]; F_GD[i] = NULL;
+            delete[] F_DG[i]; F_DG[i] = NULL;
+        }
+        delete[] F_MM; F_MM = NULL;
+        delete[] F_MI; F_MI = NULL;
+        delete[] F_IM; F_IM = NULL;
+        delete[] F_DG; F_DG = NULL;
+        delete[] F_GD; F_GD = NULL;
+        delete[] scale; scale = NULL;
     }
-  delete[] F_MM; F_MM = NULL;
-  delete[] F_MI; F_MI = NULL;
-  delete[] F_IM; F_IM = NULL;
-  delete[] F_DG; F_DG = NULL;
-  delete[] F_GD; F_GD = NULL;
-  delete[] scale; scale = NULL;
 }
 
 /////////////////////////////////////////////////////////////////////////////////////
@@ -289,14 +294,14 @@ Hit::DeleteForwardMatrix(int Nq)
 void 
 Hit::AllocateBackwardMatrix(int Nq, int Nt)
 {
-  B_MM=new(double*[Nq]);
+  B_MM=new double*[Nq];
   B_MI=F_MI; 
   B_DG=F_DG; 
   B_IM=F_IM; 
   B_GD=F_GD; 
   for (int i=0; i<Nq; i++) 
     {
-      B_MM[i] = new(double[Nt]);
+      B_MM[i] = new double[Nt];
       if (!B_MM[i]) 
        {
          fprintf(stderr,"Error: out of memory while allocating row %i (out of %i) for dynamic programming matrices \n",i+1,Nq);
@@ -311,12 +316,14 @@ Hit::AllocateBackwardMatrix(int Nq, int Nt)
 
 void Hit::DeleteBackwardMatrix(int Nq)
 {
-  for (int i=0; i<Nq; i++) 
-    {
-      delete[] B_MM[i]; B_MM[i] = NULL;  /* is this all? FS */
+
+    if (NULL != B_MM){ /* FS, r259 -> r260 */
+        for (int i=0; i<Nq; i++) {
+            delete[] B_MM[i]; B_MM[i] = NULL;  /* is this all? FS */
+        }
+        delete[] B_MM; B_MM = NULL;
+        B_MM=B_MI=B_IM=B_DG=B_GD=NULL;
     }
-  delete[] B_MM; B_MM = NULL;
-  B_MM=B_MI=B_IM=B_DG=B_GD=NULL;
 }
 
 
@@ -367,16 +374,16 @@ Hit::Viterbi(HMM& q, HMM& t, float** Sstruc)
     //float sDG[MAXRES];          // sDG[i][j] = score of best alignment up to indices (i,j) ending in (Delete,Gap)
     //float sIM[MAXRES];          // sIM[i][j] = score of best alignment up to indices (i,j) ending in (Ins,Match)
     //float sMI[MAXRES];          // sMI[i][j] = score of best alignment up to indices (i,j) ending in (Match,Ins) 
-    float *sMM = new(float[par.maxResLen]);   // sMM[i][j] = score of best alignment up to indices (i,j) ending in (Match,Match) 
-    float *sGD = new(float[par.maxResLen]);   // sGD[i][j] = score of best alignment up to indices (i,j) ending in (Gap,Delete) 
-    float *sDG = new(float[par.maxResLen]);   // sDG[i][j] = score of best alignment up to indices (i,j) ending in (Delete,Gap)
-    float *sIM = new(float[par.maxResLen]);   // sIM[i][j] = score of best alignment up to indices (i,j) ending in (Ins,Match)
-    float *sMI = new(float[par.maxResLen]);   // sMI[i][j] = score of best alignment up to indices (i,j) ending in (Match,Ins) 
-  float smin=(par.loc? 0:-FLT_MAX);  //used to distinguish between SW and NW algorithms in maximization         
-  int i,j;      //query and template match state indices
-  float sMM_i_j=0,sMI_i_j,sIM_i_j,sGD_i_j,sDG_i_j;
-  float sMM_i_1_j_1,sMI_i_1_j_1,sIM_i_1_j_1,sGD_i_1_j_1,sDG_i_1_j_1;
-  int jmin,jmax;
+    float *sMM = new float[par.maxResLen];   // sMM[i][j] = score of best alignment up to indices (i,j) ending in (Match,Match) 
+    float *sGD = new float[par.maxResLen];   // sGD[i][j] = score of best alignment up to indices (i,j) ending in (Gap,Delete) 
+    float *sDG = new float[par.maxResLen];   // sDG[i][j] = score of best alignment up to indices (i,j) ending in (Delete,Gap)
+    float *sIM = new float[par.maxResLen];   // sIM[i][j] = score of best alignment up to indices (i,j) ending in (Ins,Match)
+    float *sMI = new float[par.maxResLen];   // sMI[i][j] = score of best alignment up to indices (i,j) ending in (Match,Ins) 
+    float smin=(par.loc? 0:-FLT_MAX);  //used to distinguish between SW and NW algorithms in maximization         
+    int i=0,j=0;      //query and template match state indices
+    float sMM_i_j=0, sMI_i_j=0, sIM_i_j=0, sGD_i_j=0, sDG_i_j=0;
+    float sMM_i_1_j_1=0, sMI_i_1_j_1=0, sIM_i_1_j_1=0, sGD_i_1_j_1=0, sDG_i_1_j_1=0;
+    int jmin=0, jmax=0;
 
   // Reset crossed out cells?
   if(irep==1) InitializeForAlignment(q,t);
@@ -690,8 +697,8 @@ Hit::Forward(HMM& q, HMM& t, float** Pstruc)
         }
     
     // Calculate log2(P_forward)
-    score = log2(Pforward)-10.0f;
-    for (i=1; i<=q.L+1; i++) score -= log2(scale[i]);
+    score = Log2(Pforward)-10.0f;
+    for (i=1; i<=q.L+1; i++) score -= Log2(scale[i]);
     //   state = MM;
     
     if (par.loc) 
@@ -732,6 +739,12 @@ Hit::Forward(HMM& q, HMM& t, float** Pstruc)
                 __FUNCTION__, __FILE__, __LINE__, score, Pforward);
         return FAILURE;
     }
+    /* alignment might fail if no useful characters in sequence, FS, r259 -> r260 */
+    if ( (q.L <= 0) || (t.L <= 0) ){
+        fprintf(stderr, "%s:%s:%d: length(s) of profile(s) invalid (q.L=%d/t.L=%d)\n",
+                __FUNCTION__, __FILE__, __LINE__, q.L, t.L);
+        return FAILURE;
+    }
     i = q.L-1; j = t.L-1; /* FS, r241 -> r243 */
     if (isinf(F_MM[i][j]+F_MI[i][j]+F_IM[i][j]+F_DG[i][j]+F_GD[i][j])){
         fprintf(stderr, "%s:%s:%d: F_MM[i][j]=%g, F_IM[i][j]=%g, F_MI[i][j]=%g, F_DG[i][j]=%g, F_GD[i][j]=%g (i=%d,j=%d)\n", 
@@ -1111,9 +1124,9 @@ Hit::Backtrace(HMM& q, HMM& t)
   nsteps=step; 
   
   // Allocate new space for alignment scores
-  if (t.Xcons) Xcons = new( char[q.L+2]); // for template consensus sequence aligned to query
-  S    = new( float[nsteps+1]);
-  S_ss = new( float[nsteps+1]);
+  if (t.Xcons) Xcons = new char[q.L+2]; // for template consensus sequence aligned to query
+  S    = new float[nsteps+1];
+  S_ss = new float[nsteps+1];
   if (!S_ss) MemoryError("space for HMM-HMM alignments");
 
   // Add contribution from secondary structure score, record score along alignment,
@@ -1213,7 +1226,7 @@ Hit::StochasticBacktrace(HMM& q, HMM& t, char maximize)
   int i,j;         // query and template match state indices
 //  float pmin=(par.loc? 1.0: 0.0);    // used to distinguish between SW and NW algorithms in maximization         
   const float pmin=0;
-  double* scale_cum = new(double[q.L+2]);
+  double* scale_cum = new double[q.L+2];
   
 
   scale_cum[0]=1;
@@ -1232,7 +1245,7 @@ Hit::StochasticBacktrace(HMM& q, HMM& t, char maximize)
   else 
     {
 //      float sumF[q.L+t.L];
-      double* sumF=new(double[q.L+t.L]);
+      double* sumF=new double[q.L+t.L];
       sumF[0]=0.0;
       for (j=1; j<=t.L; j++)        sumF[j] = sumF[j-1] + F_MM[q.L][j]/scale_cum[q.L];;
       for (j=t.L+1; j<t.L+q.L; j++) sumF[j] = sumF[j-1] + F_MM[j-t.L][t.L]/scale_cum[j-t.L];;
@@ -1356,9 +1369,9 @@ Hit::StochasticBacktrace(HMM& q, HMM& t, char maximize)
   nsteps=step; 
 
   // Allocate new space for alignment scores
-  if (t.Xcons) Xcons = new( char[q.L+2]); // for template consensus sequence aligned to query
-  S    = new( float[nsteps+1]);
-  S_ss = new( float[nsteps+1]);
+  if (t.Xcons) Xcons = new char[q.L+2]; // for template consensus sequence aligned to query
+  S    = new float[nsteps+1];
+  S_ss = new float[nsteps+1];
   if (!S_ss) MemoryError("space for HMM-HMM alignments");
 
   // Add contribution from secondary structure score, record score along alignment,
@@ -1460,10 +1473,10 @@ Hit::BacktraceMAC(HMM& q, HMM& t)
   nsteps=step; 
     
   // Allocate new space for alignment scores
-  if (t.Xcons) Xcons = new( char[q.L+2]); // for template consensus sequence aligned to query
-  S    = new( float[nsteps+1]);
-  S_ss = new( float[nsteps+1]);
-  P_posterior = new( float[nsteps+1]);
+  if (t.Xcons) Xcons = new char[q.L+2]; // for template consensus sequence aligned to query
+  S    = new float[nsteps+1];
+  S_ss = new float[nsteps+1];
+  P_posterior = new float[nsteps+1];
   if (!P_posterior) MemoryError("space for HMM-HMM alignments");
 
   // Add contribution from secondary structure score, record score along alignment,
@@ -1665,59 +1678,62 @@ Hit::InitializeForAlignment(HMM& q, HMM& t)
 void 
 Hit::InitializeBacktrace(HMM& q, HMM& t)
 {
-  if (irep==1) //if this is the first single repeat repeat hit with this template
-    {
-      //Copy information about template profile to hit and reset template pointers to avoid destruction
-      longname=new(char[strlen(t.longname)+1]);
-      name    =new(char[strlen(t.name)+1]);
-      file    =new(char[strlen(t.file)+1]);
-      if (!file) MemoryError("space for alignments with database HMMs. \nNote that all alignments have to be kept in memory");
-      strcpy(longname,t.longname);
-      strcpy(name,t.name);
-      strcpy(fam ,t.fam);
-      strcpy(sfam ,t.sfam);
-      strcpy(fold ,t.fold);
-      strcpy(cl ,t.cl);
-      strcpy(file,t.file);
-      sname=new(char*[t.n_display]);   // Call Compare only once with irep=1
-      seq  =new(char*[t.n_display]);   // Call Compare only once with irep=1
-      if (!sname || !seq) 
-       MemoryError("space for alignments with database HMMs.\nNote that all sequences for display have to be kept in memory");
-
-      for (int k=0; k<t.n_display; k++)        {
-          if (NULL != t.sname){
-              sname[k]=t.sname[k]; t.sname[k]=NULL;
-          }
-          else {
-              sname[k]=NULL;
-          }
-          seq[k]  =t.seq[k];   t.seq[k]=NULL;
-      }
+    if (irep==1) //if this is the first single repeat repeat hit with this template
+        {
+            //Copy information about template profile to hit and reset template pointers to avoid destruction
+            longname=new char[strlen(t.longname)+1]();
+            name    =new char[strlen(t.name)+1]();
+            file    =new char[strlen(t.file)+1]();
+            if (!file) {
+                MemoryError("space for alignments with database HMMs. \nNote that all alignments have to be kept in memory");
+            }
+            strcpy(longname,t.longname);
+            strcpy(name,t.name);
+            strcpy(fam ,t.fam);
+            strcpy(sfam ,t.sfam);
+            strcpy(fold ,t.fold);
+            strcpy(cl ,t.cl);
+            strcpy(file,t.file);
+            sname=new char*[t.n_display]();   // Call Compare only once with irep=1
+            seq  =new char*[t.n_display]();   // Call Compare only once with irep=1
+            if (!sname || !seq) {
+                MemoryError("space for alignments with database HMMs.\nNote that all sequences for display have to be kept in memory");
+            }
 
-      n_display=t.n_display; t.n_display=0;
-      ncons  = t.ncons;
-      nfirst = t.nfirst;
-      nss_dssp = t.nss_dssp;
-      nsa_dssp = t.nsa_dssp;
-      nss_pred = t.nss_pred;
-      nss_conf = t.nss_conf;
-      L = t.L;
-      Neff_HMM = t.Neff_HMM;
-      Eval   = 1.0;
-      Pval   = 1.0;
-      Pvalt  = 1.0;
-      logPval = 0.0;
-      logPvalt= 0.0;
-      Probab = 1.0;
-    }    
-
-  // Allocate new space
-  this->i = new( int[i2+j2+2]);
-  this->j = new( int[i2+j2+2]);
-  states  = new( char[i2+j2+2]);
-  S = S_ss = P_posterior = NULL; // set to NULL to avoid deleting data from irep=1 when hit with irep=2 is removed 
-  Xcons = NULL;
-}
+            for (int k=0; k<t.n_display; k++)  {
+                if (NULL != t.sname){
+                    sname[k]=t.sname[k]; t.sname[k]=NULL;
+                }
+                else {
+                    sname[k]=NULL;
+                }
+                seq[k]  =t.seq[k];   t.seq[k]=NULL;
+            }
+            
+            n_display=t.n_display; t.n_display=0;
+            ncons  = t.ncons;
+            nfirst = t.nfirst;
+            nss_dssp = t.nss_dssp;
+            nsa_dssp = t.nsa_dssp;
+            nss_pred = t.nss_pred;
+            nss_conf = t.nss_conf;
+            L = t.L;
+            Neff_HMM = t.Neff_HMM;
+            Eval   = 1.0;
+            Pval   = 1.0;
+            Pvalt  = 1.0;
+            logPval = 0.0;
+            logPvalt= 0.0;
+            Probab = 1.0;
+        }    
+    
+    // Allocate new space
+    this->i = new int[i2+j2+2]();
+    this->j = new int[i2+j2+2]();
+    states  = new char[i2+j2+2]();
+    S = S_ss = P_posterior = NULL; // set to NULL to avoid deleting data from irep=1 when hit with irep=2 is removed 
+    Xcons = NULL;
+} /* this is the end of Hit::InitializeBacktrace() */
 
 /////////////////////////////////////////////////////////////////////////////////////
 // Some score functions