Removing old version of mafft to replace with the new
[jabaws.git] / binaries / src / mafft / extensions / mxscarna_src / McCaskill.cpp
diff --git a/binaries/src/mafft/extensions/mxscarna_src/McCaskill.cpp b/binaries/src/mafft/extensions/mxscarna_src/McCaskill.cpp
deleted file mode 100644 (file)
index 39497fb..0000000
+++ /dev/null
@@ -1,839 +0,0 @@
-#include <iostream>
-#include "McCaskill.hpp"
-//#include "energy_param3.hpp"
-#include "Util.hpp"
-#include <cstring>
-
-namespace MXSCARNA {
-energy_param  McCaskill::energyParam;
-
-float *McCaskill::exphairpin;
-float McCaskill::expninio[5][32];
-float McCaskill::expdangle5[6][4];
-float McCaskill::expdangle3[6][4];
-float McCaskill::expinternalLoop[31];
-float McCaskill::expbulge[31];
-char   McCaskill::exptriLoop[2][6];
-float McCaskill::exptriLoopEnergy[2];
-char   McCaskill::exptetraLoop[30][7];
-float McCaskill::exptetraLoopEnergy[30];
-float McCaskill::expStack[6][6];
-float McCaskill::expTstackH[6][16];
-float McCaskill::expTstackI[6][16];
-float McCaskill::expint11[6][6][16];
-float McCaskill::expint21[6][6][16][4];
-float McCaskill::expint22[6][6][16][16];
-float McCaskill::expMLclosing;
-float McCaskill::expMLintern[8];
-float McCaskill::expTermAU;
-float McCaskill::expMLbase[31];
-
-const int McCaskill::TURN         = 3;
-const float McCaskill::GASCONST   = 1.98717;
-const float McCaskill::T          = 37 + 273.15;
-const int McCaskill::MAXLOOP      = 30;
-const int McCaskill::TETRA_ENTH37 = -400;
-const int McCaskill::NBPAIRS      = 7;
-const int McCaskill::SCALE        = 10;
-
-
-void 
-McCaskill::
-calcPartitionFunction()
-{ 
-    initParameter();
-    Inside();
-    Outside();
-    convertProbability();
-
-/*
-    for(int i = 0; i < n_seq; i++) {
-       for(int j = i; j < n_seq; j++) {
-           cout << getProb(i, j) << " ";
-       }
-       cout << endl;
-    }
-*/
-}
-
-void
-McCaskill::
-convertProbability()
-{
-    float *pPointer  = p.getPointer(0, 0);
-    float *abPointer = ab.getPointer(0,0);
-    for(int i = 0; i < n_seq; i++) {
-       for(int j = i; j < n_seq; j++) {
-           *pPointer += *abPointer++;
-           *pPointer++ = EXP(*pPointer);
-       }
-    }
-}
-
-void 
-McCaskill::
-Inside()
-{
-
-    for (int i = n_seq - TURN - 2; i >= 0; i--) {
-       float *abPointer   = ab.getPointer(i, i + TURN + 1);
-       float *am1Pointer  = am1.getPointer(i, i + TURN + 1);
-       float *amPointer   = am.getPointer(i, i + TURN + 1);
-       float *q1Pointer   = q1.getPointer(i, i + TURN + 1);
-       float *aPointer    = a.getPointer(i, i + TURN + 1);
-       int   *typePointer = typeMat.getPointer(i, i+TURN+1);
-       for (int j = i + TURN + 1; j < n_seq; j++) {
-           int tmpType   = *typePointer++;
-                           *abPointer++  = compQb(i, j, tmpType);
-           am1v.ref(i,j) = *am1Pointer++ = compQm1(i,j, tmpType);
-           amv.ref(i,j)  = *amPointer++  = compQm(i,j);
-           q1v.ref(i,j)  = *q1Pointer++  = compQ1(i,j, tmpType);
-                           *aPointer++   = compQ(i,j);   
-       }
-    }
-}
-
-inline float
-McCaskill::
-compQb(int i, int j, int tmpType) 
-{
-
-    float tmpAb;
-    int type  = tmpType;
-    int u     = j - i - 1;
-    if(Beta::isCanonicalReducedPairCode(type)) {
-       // hairpin energy 
-       assert(u >= 3);
-       tmpAb = expHairpinEnergy(type, u, i + 1, j - 1);
-               
-       // base pair, bulge, interior loop energy
-       for(int h = i + 1; h <= MIN(i + MAXLOOP + 1, j - TURN - 2); h++) {
-           int u1 = h-i-1;
-           int max = MAX(h + TURN + 1, j-1-MAXLOOP+u1);
-           float *abPointer     = ab.getPointer(h, max - 1);
-           const int *typeMatPointer = typeMat.getPointer(h, max);
-
-           for(int l = max; l < j; l++) {
-               int type2 = *typeMatPointer++;
-               abPointer++;
-               if(!Beta::isCanonicalReducedPairCode(type2)) continue;
-               
-               assert(h >= 0 && h < n_seq && l >= 0 && l < n_seq);
-               type2 = Beta::flipReducedPairCode(type2);
-               assert(h-i-1 >= 0); assert(j-l-1 >= 0);
-               float loopE = *abPointer;
-               loopE += expLoopEnergy(u1, j-l-1, tmpType, type2, i+1, j-1, h-1, l+1); 
-               tmpAb = LOG_ADD(tmpAb, loopE);
-           }
-       }
-
-       // multi loop
-       float tmpQm = IMPOSSIBLE;
-       float *amPointer  = am.getPointer(i + 1, j - TURN - 3);
-       float *am1Pointer = am1v.getPointer(j-TURN-2, j-1);
-       for(int h = j - TURN - 2; h >= i + TURN + 3; h--) {
-           assert(h >= 0 && h < n_seq);
-           float tmpScore = *amPointer--;
-           tmpScore += *am1Pointer--;
-           tmpQm = LOG_ADD(tmpQm, tmpScore);
-       }
-       tmpQm += expMLclosing + expMLintern[type];
-       tmpQm += endStemScore(i, j);
-       tmpAb = LOG_ADD(tmpAb, tmpQm);
-    }
-    else {
-       tmpAb = IMPOSSIBLE;
-    }
-    return tmpAb;
-}
-
-//F = a:ML_closing + b:ML_intern*k + c:ML_BASE*u
-
-inline float
-McCaskill::
-compQm1(int i, int j, int tmpType)
-{
-    float tmpQm1 = IMPOSSIBLE;
-
-    int l = j;
-    if (i + TURN + 1 <= l) {
-       int type = typeMat.ref(i,l);
-       if(Beta::isCanonicalReducedPairCode(type)) {
-           float tmpScore = ab.ref(i,l);
-           tmpScore += beginStemScore(i, l);
-           //tmpScore += expMLbase[1]*(j-l) + expMLintern[type];
-           tmpScore += expMLintern[tmpType];
-           tmpQm1 = LOG_ADD(tmpQm1, tmpScore);
-       }
-    }
-    if(i < j) {
-       tmpQm1 = LOG_ADD(tmpQm1, am1.ref(i,j-1));
-    }
-
-    return tmpQm1;
-}
-
-inline float
-McCaskill::
-compQm(int i, int j)
-{
-    float tmpQm = IMPOSSIBLE;
-    float *amPointer  = am.getPointer(i,j-TURN-2);
-    float *am1Pointer = am1v.getPointer(j-TURN-1, j);
-    for(int h = j - TURN - 1; h >= i ; h--) {
-       float tmpScore = 0;
-       float tmpAm1 = *am1Pointer--;
-       
-       tmpScore += tmpAm1;
-       tmpQm = LOG_ADD(tmpQm, tmpScore);
-       tmpScore = *amPointer--;
-       tmpScore += tmpAm1;
-       tmpQm = LOG_ADD(tmpQm, tmpScore);
-    }
-
-    return tmpQm;
-}
-
-inline float
-McCaskill::
-compQ1(int i, int j, int tmpType)
-{
-    float tmpQ1 = IMPOSSIBLE;
-
-    if(Beta::isCanonicalReducedPairCode(tmpType)) {
-       float tmpScore = ab.ref(i, j);
-       tmpScore += beginStemScore(i, j);
-       tmpQ1 = LOG_ADD(tmpQ1, tmpScore);
-    }
-    tmpQ1 = LOG_ADD(tmpQ1, q1.ref(i, j - 1));
-
-    return tmpQ1;
-}
-
-inline float
-McCaskill::
-compQ(int i, int j)
-{
-    float tmpQ = 0;
-    tmpQ = LOG_ADD(tmpQ, q1.ref(i,j));
-    
-    float *aPointer  = a.getPointer(i,j-TURN-2);
-    float *q1Pointer = q1v.getPointer(j-TURN-1, j);
-    for(int h = j - TURN - 1; h >= i + 1; h--) {
-       float tmpScore  = *aPointer--;
-       tmpScore       += *q1Pointer--;
-       tmpQ            = LOG_ADD(tmpQ, tmpScore);
-    }
-
-    return tmpQ;
-}
-
-inline float
-McCaskill::
-beginStemScore(const int i, const int j) const
-{
-    float temp = 0;
-    int type = typeMat.ref(i,j);
-
-    if(0 < i)                   { temp += expdangle5[type][numSeq[i-1]]; }
-    if(j < n_seq-1)             { temp += expdangle3[type][numSeq[j+1]]; }
-    if(type != Beta::REDUCED_CG_CODE && type != Beta::REDUCED_GC_CODE)  { temp += expTermAU; }
-    return temp;
-}
-
-inline float
-McCaskill::
-endStemScore(const int i, const int j) const
-{
-    float temp = 0;
-    int type = typeMat.ref(i,j);
-
-    type = Beta::flipReducedPairCode(type);
-
-    if(i < n_seq-1)             { temp += expdangle3[type][numSeq[i+1]]; }
-    if(j > 0)                   { temp += expdangle5[type][numSeq[j-1]]; }
-    if(type != Beta::REDUCED_CG_CODE && type != Beta::REDUCED_GC_CODE)  { temp += expTermAU; }
-    return temp;
-}
-
-inline float
-McCaskill::
-compP(int h, int l, int tmpType)
-{
-    float prob = IMPOSSIBLE;
-           
-    int type = tmpType;
-    if(Beta::isCanonicalReducedPairCode(type)) {
-               
-       /* exterior loop */
-       float tmp_p = 0;
-       tmp_p -= a.ref(0,n_seq-1);
-       if(0 < h) { 
-           tmp_p += a.ref(0,h-1);
-       }
-       if(l < n_seq-1) {
-           tmp_p += a.ref(l+1, n_seq-1);
-       }
-       tmp_p += beginStemScore(h, l);
-       prob = LOG_ADD(prob, tmp_p);
-
-       assert(IMPOSSIBLE <= prob && prob <= 0);
-
-       /* internal loop */
-       tmp_p = IMPOSSIBLE;
-       int tt = Beta::flipReducedPairCode(tmpType);
-       int max = MAX(0,h-MAXLOOP-1);
-       for(int i = max; i <= h - 1; i++) {
-           float min = MIN(l+MAXLOOP-h+i+2, n_seq-1);
-           int   *typeMatPointer = typeMat.getPointer(i,l+1);
-           float *pPointer       = p.getPointer(i,l);
-           for(int j = l + 1; j <= min; j++) {
-               int type2    = *typeMatPointer++;
-               pPointer++;
-               if(!Beta::isCanonicalReducedPairCode(type2)) continue;
-               assert(i >= 0 && i < n_seq && j >= 0 && j < n_seq);
-
-               float tmpScore  = *pPointer;
-               tmpScore       += expLoopEnergy(h-i-1, j-l-1, type2, tt, i+1, j-1, h-1, l+1);
-               tmp_p = LOG_ADD(tmp_p, tmpScore);
-           }
-       }
-       prob = LOG_ADD(prob, tmp_p); 
-       assert(IMPOSSIBLE <= prob && prob <= 0);
-
-       /* multi loop */
-       tmp_p = IMPOSSIBLE;
-       float tmp_begin   = beginStemScore(h, l);
-       float *q1Pointer  = q1v.getPointer(0, l);
-       float *am1Pointer = am1v.getPointer(0, l);
-       float *amPointer  = amv.getPointer(1,h-1);
-       for(int i = 0; i <= h-TURN-1; i++) {
-           float tmpq1    = *q1Pointer++;
-           float tmpam    = *amPointer++;
-           float tmpScore = *am1Pointer++;
-
-           tmpScore += tmpam;
-           tmpScore += tmp_begin;
-           tmpScore += expMLclosing + expMLintern[tt];
-           tmp_p = LOG_ADD(tmp_p, tmpScore);
-
-           tmpScore  = tmpq1;
-           tmpScore += tmpam;
-           tmpScore += tmp_begin;
-           tmpScore += expMLclosing + expMLintern[tt];
-           tmp_p = LOG_ADD(tmp_p, tmpScore);
-
-           tmpScore  = tmpq1;
-           tmpScore += tmp_begin;
-           tmpScore += expMLclosing + expMLintern[tt];
-           tmp_p = LOG_ADD(tmp_p, tmpScore);
-       }
-               
-       assert(IMPOSSIBLE <= tmp_p && tmp_p <= 0);
-       prob = LOG_ADD(prob, tmp_p); 
-       
-       tmp_p = IMPOSSIBLE;
-       for(int i = h-TURN; i <= h-1; i++) {
-           if(i >= 0) {
-               float tmpScore  = q1.ref(i,l);
-               tmpScore       += tmp_begin;
-               tmpScore       += expMLclosing + expMLintern[tt];
-               tmp_p           = LOG_ADD(tmp_p, tmpScore);
-           }
-       }
-       assert(IMPOSSIBLE <= tmp_p && tmp_p <= 0); 
-       prob = LOG_ADD(prob, tmp_p);
-    }
-    else {
-       prob = IMPOSSIBLE;
-    }
-
-    return prob;
-}
-
-inline float
-McCaskill::
-compPm(int i, int l)
-{
-  float tmpPm  = IMPOSSIBLE;
-
-  int *typeMatPointer   = typeMat.getPointer(i,n_seq-1);
-  float *pPointer       = p.getPointer(i,n_seq);
-  float *amPointer      = am.getPointer(l+1,n_seq-1);
-  float *abPointer      = ab.getPointer(i, n_seq);
-  for(int j = n_seq - 1; j >= l + TURN + 1; j--) {
-      int type = *typeMatPointer--;
-      pPointer--;
-      amPointer--;
-      abPointer--;
-      if(Beta::isCanonicalReducedPairCode(type)) {
-         float tmp  = *pPointer;
-         tmp       += *amPointer;
-         tmp       += endStemScore(i, j);
-         tmpPm = LOG_ADD(tmpPm, tmp);
-      }
-  }
-  tmpPm += expMLintern[1];
-
-  return tmpPm;
-}
-
-inline float
-McCaskill::
-compPm1(int i, int l)
-{
-  float tmpPm1 = IMPOSSIBLE;
-
-  int j =  l + 1;
-  if(j <= n_seq-1) {
-      int type = typeMat.ref(i,j);
-      if(Beta::isCanonicalReducedPairCode(type)) {
-         float tmp = p.ref(i,j);
-         tmp += endStemScore(i, j);
-         tmpPm1 = tmp;
-      }
-      tmpPm1 += expMLintern[1];
-  }
-  if(l+1 <= n_seq - 1) {
-      tmpPm1 = LOG_ADD(tmpPm1, am1.ref(i, l+1));
-  }
-
-  return tmpPm1;
-}
-
-void 
-McCaskill::
-Outside()
-{
-    for(int h = 0; h <= n_seq - TURN - 2; h++) {
-       float *pPointer    = p.getPointer(h, n_seq-1);
-       float *q1Pointer   = q1.getPointer(h, n_seq-1);
-       float *am1Pointer  = am1.getPointer(h, n_seq-1);
-       int   *typePointer = typeMat.getPointer(h, n_seq-1);
-       for(int l = n_seq-1; l >= h + TURN + 1; l--) {
-           int tmpType    = *typePointer--;
-           pv.ref(h,l)    = *pPointer--   = compP(h,l,tmpType);
-           q1v.ref(h,l)   = *q1Pointer--  = compPm(h,l);
-           am1v.ref(h,l)  = *am1Pointer-- = compPm1(h,l);
-
-           assert(p.ref(h,l) <= 0);
-       }
-    }
-}
-
-void
-McCaskill::
-printProbMat()
-{
-    int m = 0;
-    for(int i = 0; i < n_seq; i++) cout << " " << seq[i];
-    cout << endl;
-    for(int i = 0; i < n_seq; i++) {
-       if(m < n_seq) {
-           cout << seq[m];
-       }
-       for(int j = 0; j <= i-1; j++) {
-           if(j != i-1) cout << "  ";
-           else         cout << " ";
-       }
-       if(i != 0 && i != n_seq-1) {
-           cout << "\\";
-       }
-
-       for(int j = i; j < n_seq; j++) {
-           if(p.ref(i,j) > 0.01) {
-
-               int type = Beta::getReducedPairCode(numSeq[i], numSeq[j]);
-               
-               if(!Beta::isCanonicalReducedPairCode(type)) {
-                   cout << "\n" << seq[i] << " " << seq[j] << " " << exp(p.ref(i,j)) << endl;
-               }
-
-               if(j != n_seq-1) {
-                   cout << "* ";
-               }
-               else {
-                   cout << "*";
-               }
-
-           }
-           else {
-
-               if(j != n_seq-1) {
-                   cout << "  ";
-               }
-               else {
-                   cout << " ";
-               }
-
-           }
-
-       }
-       if(m < n_seq) {
-           cout << seq[m++] << endl;
-       }
-       if(i == n_seq - 1) cout << endl;
-    }
-    for(int i = 0; i < n_seq; i++) cout << " " << seq[i];
-    cout << endl;
-}
-
-void
-McCaskill::
-initParameter()
-{
-    float GT;
-    float RT_KCAL_MOL = McCaskill::T*McCaskill::GASCONST;
-    int len = 31;
-
-    for(int i = 0; i < len; i++) {
-       GT = energyParam.getHairpin(i);
-       exphairpin[i] = -GT*10/RT_KCAL_MOL;
-    }
-    
-    for (int i = 0; i < len; i++) {
-       GT = energyParam.getBulge(i);
-       expbulge[i] = -GT*10/RT_KCAL_MOL;
-       GT = energyParam.getInternalLoop(i);
-       expinternalLoop[i] = -GT*10/RT_KCAL_MOL;
-    }
-    expinternalLoop[2] = -80*10/RT_KCAL_MOL; /* special case of size 2 interior loops (single mismatch) */
-    
-    // float lxc = energy_param3::lxc37;
-    for (int i = 31; i < n_seq; i++) {
-       GT = energyParam.getHairpin(30) + (107.856*LOG((float)i/30));
-       exphairpin[i] = -GT*10/RT_KCAL_MOL;
-    }
-
-    for(int i = 0; i < 5; i++) {
-       GT = energyParam.getNinio(i);
-       for(int j = 0; j <= MAXLOOP; j++) {
-           expninio[i][j] = -MIN(energyParam.getMaxNinio(), j*GT)*10/RT_KCAL_MOL;
-       }
-    }
-
-    for(int i = 0; i < 30; i++) {
-      GT = energyParam.getTetraLoopEnergy(i);
-      exptetraLoopEnergy[i] = -GT*10/RT_KCAL_MOL;
-    }
-    
-    /*no parameters for Triloop*/
-    for(int i = 0; i < 2; i++) {
-       GT = 0;
-       exptriLoopEnergy[i] = -GT*10/RT_KCAL_MOL;
-    }
-
-    GT = energyParam.getMLclosing();
-    expMLclosing = -GT*10/RT_KCAL_MOL;
-
-    for(int i = 0; i <= NBPAIRS; i++) {
-      GT = energyParam.getMLintern();
-      expMLintern[i] = -GT*10/RT_KCAL_MOL;
-    }
-
-    expTermAU = -energyParam.getTerminalAU()*10/RT_KCAL_MOL;
-    
-    GT = energyParam.getMLBase();
-    for(int i = 0; i < len; i++) {
-       expMLbase[i] = -GT*10*(float)i/RT_KCAL_MOL;
-    }
-
-    /*
-       if danlges = 0 just set their energy to 0,
-       don't let dangle energyies become > 0 (at large temps)
-       but make sure go smoothly to 0
-    */
-    for(int i = 0; i < 6; i++) {
-       for(int j =0; j < 4; j++) {
-           GT = energyParam.getDangle5(i,j);
-           expdangle5[i][j] = -GT*10/RT_KCAL_MOL;
-           GT = energyParam.getDangle3(i,j);
-           expdangle3[i][j] = -GT*10/RT_KCAL_MOL;
-       }
-    }
-
-    /* stacking energies  */
-    for(int i = 0; i < 6; i++) {
-       for(int j = 0; j < 6; j++) {
-           GT = energyParam.getStack(i,j);
-           expStack[i][j] = -GT*10/RT_KCAL_MOL;
-       }
-    }
-
-    /* mismatch energies */
-    for (int i = 0; i < 6; i++) {
-       for (int j = 0; j < 16; j++) {
-         GT = energyParam.getTstackI(i, j);
-         //  cout << i << " " << " " << j << " " << GT << endl;
-         expTstackI[i][j] = -GT*10/RT_KCAL_MOL;
-         GT = energyParam.getTstackH(i, j);
-         expTstackH[i][j] = -GT*10/RT_KCAL_MOL;
-       }
-    }
-    
-    /* interior loops of length 2*/
-    for(int i = 0; i < 6; i++) {
-       for(int j = 0; j < 6; j++) {
-           for(int k = 0; k < 16; k++) {
-             GT = energyParam.getInt11(i, j, k);
-             expint11[i][j][k] = -GT*10/RT_KCAL_MOL;
-           }
-       }
-    }
-
-    /* interior 2*1 loops */
-    for(int i = 0; i < 6; i++) {
-       for(int j =0; j < 6; j++) {
-           for(int k = 0; k < 16; k++) {
-               for(int l = 0; l < 4; l++) {
-                 GT = energyParam.getInt21(i,j,k,l);
-                 expint21[i][j][k][l] = -GT*10/RT_KCAL_MOL;
-               }
-           }
-       }
-    }
-
-    /* interior 2*2 loops */
-    for (int i = 0; i < 6; i++) {
-       for(int j = 0; j < 6; j++) {
-           for(int k = 0; k < 16; k++) {
-               for(int l = 0; l < 16; l++) {
-                 GT = energyParam.getInt22(i,j,k,l);
-                 expint22[i][j][k][l] = -GT*10/RT_KCAL_MOL;
-               }
-           }
-       }
-    }
-}
-
-
-inline float 
-McCaskill::
-expHairpinEnergy(const int type, const int l, const int i, const int j)
-{
-  float q;
-  int    k;
-  
-//  assert(l >= 0);
-  q = exphairpin[l];
-
-  if(l == 4) {
-      char temp_seq[7];
-
-      for(int iter = i - 1; iter < i + 5; iter++) {
-       temp_seq[iter - (i-1)] = seq[iter];
-      }
-      temp_seq[6] = '\0';
-
-      for(k = 0; k < 30; k++) {
-       if(strcmp(temp_seq, energyParam.getTetraLoop(k)) == 0) break;
-      }
-      if(k != 30) {
-       q += exptetraLoopEnergy[k];
-      }
-  }
-  if(l == 3) {
-
-      /* no triloop bonus
-    char temp_seq[6];
-    
-    for(int iter = i - 1; iter < i + 4; iter++) {
-       temp_seq[iter - (i-1)] = seq[iter];
-    }
-    temp_seq[6] = '\0';
-    for(k = 0; k < 2; k++) {
-      if(strcmp(temp_seq, energyParam.getTriLoop(k)) == 0) break;
-    }
-    if(k != 2) {
-      q *= exptriLoopEnergy[k];
-    }
-      */
-
-    if(type != Beta::REDUCED_CG_CODE && type != Beta::REDUCED_GC_CODE) q += expTermAU;
-  }
-  else {
-      int type2 = Beta::getPairCode(numSeq[i], numSeq[j]);
-      q += expTstackH[type][type2];
-  }
-
-  return q;
-}
-
-
-inline float 
-McCaskill::
-expLoopEnergy(int u1, int u2, int type, int type2, 
-             int si1, int sj1, int sp1, int sq1)
-{
-  float z = 0;
-
-  if((u1 == 0) && (u2 == 0)) { z = expStack[type][type2]; }
-  else {
-    if((u1 == 0) || (u2 == 0)) {
-      int u;
-      if(u1 == 0) { u = u2; }
-      else        { u = u1; }
-      z = expbulge[u];
-      if(u1 + u2 == 1) z += expStack[type][type2];
-      else {
-       if (type  != Beta::REDUCED_CG_CODE && type  != Beta::REDUCED_GC_CODE) z += expTermAU;
-       if (type2 != Beta::REDUCED_CG_CODE && type2 != Beta::REDUCED_GC_CODE) z += expTermAU;
-      }
-    }
-    else {
-      if(u1 + u2 == 2) {
-         z = expint11[type][type2][Beta::getPairCode(numSeq[si1], numSeq[sj1])];
-      }
-      else if((u1 == 1) && (u2 == 2)) {
-         z = expint21[type][type2][Beta::getPairCode(numSeq[si1], numSeq[sj1])][numSeq[sq1]];
-      }
-      else if((u1 == 2) && (u2 == 1)) {
-         z = expint21[type2][type][Beta::getPairCode(numSeq[sq1], numSeq[sp1])][numSeq[si1]];
-      }
-      else if((u1 == 2) && (u2 == 2)) {
-       z = expint22[type][type2][Beta::getPairCode(numSeq[si1], numSeq[sj1])][Beta::getPairCode(numSeq[sp1], numSeq[sq1])];
-      }
-      else {
-       z = expinternalLoop[u1 + u2] +
-           expTstackI[type][Beta::getPairCode(numSeq[si1], numSeq[sj1])]
-            + expTstackI[type2][Beta::getPairCode(numSeq[sq1], numSeq[sp1])];
-       z += expninio[2][abs(u1-u2)];
-      }
-    }
-  }
-
-  return z;
-}
-
-void
-McCaskill::
-printExpEnergy()
-{
-    cout << "exphairpin:" << endl;
-    for(int i = 0; i < 31; i++) {
-       cout << exphairpin[i] << endl;
-    }
-    
-    cout << "expninio[5][32]:" << endl;
-    for(int i = 0; i < 5; i++) {
-       for(int j = 0; j < 32; j++) {
-           cout << expninio[i][j] << " ";
-       }
-       cout << endl;
-    }
-
-    cout << "expdangle5[6][4]:" << endl;
-    for(int i = 0; i < 6; i++) {
-       for(int j = 0; j < 4; j++) {
-           cout << expdangle5[i][j] << " ";
-       }
-       cout << endl;
-    }
-
-    cout << "expdangle3[6][4]:" << endl;
-    for(int i = 0; i < 6; i++) {
-       for(int j = 0; j < 4; j++) {
-           cout << expdangle3[i][j] << " ";
-       }
-       cout << endl;
-    }
-
-    cout << "expinternalLoop[31]:" << endl;
-    for(int i = 0; i < 31; i++) {
-       cout << i << ":" << expinternalLoop[i] << endl;
-    }
-    cout << "expbulge[31]:" << endl;
-    for(int i = 0; i < 31; i++) {
-       cout << i << ":" << expbulge[i] << endl;
-    }
-    
-    cout << "exptriLoopEnergy[2]:" << endl;
-    for(int i = 0; i < 2; i++) {
-       cout << i << ":" << exptriLoopEnergy[i] << endl;
-    }
-
-    cout << "exptetraLoopEnergy[15]" << endl;
-    for(int i = 0; i < 15; i++) {
-       cout << i << ":" << exptetraLoopEnergy[i] << endl;
-    }
-    
-    cout << "expStack[6][6]:" << endl;
-    for(int i = 0; i < 6; i++) {
-       for(int j = 0; j < 6; j++) {
-           cout << expStack[i][j] << " ";
-       }
-       cout << endl;
-    }
-
-    cout << "expTstackH[6][16]:" << endl;
-    for(int i = 0; i < 6; i++) {
-       for(int j = 0; j < 16; j++) {
-           cout << expTstackH[i][j] << " ";
-       }
-       cout << endl;
-    }
-
-    cout << "expTstackI[6][16]:" << endl;
-    for(int i = 0; i < 6; i++) {
-       for(int j = 0; j < 16; j++) {
-           cout << expTstackI[i][j] << " ";
-       }
-       cout << endl;
-    }
-    
-    cout << "expMLclosing=" << expMLclosing << endl;
-    cout << "expMLintern:" << endl;
-    for(int i = 0; i < 8; i++) {
-       cout << expMLintern[i] << " ";
-    }
-    cout << endl;
-
-    cout << "expMLbase[31]:";
-    for(int i = 0; i < 31; i++) {
-       cout << i << ":" << expMLbase[i] << endl;
-    }
-
-    cout << "expint11[6][6][16]:";
-    for(int i = 0; i < 6; i++) {
-       for(int j = 0; j < 6; j++) {
-           for(int k = 0; k < 16; k++) {
-               cout << expint11[i][j][k] << " ";
-           }
-           cout << endl;
-       }
-       cout << endl;
-    }
-
-    cout << "expint21[6][6][16][4]:" << endl;
-    for(int i = 0; i < 6; i++) {
-       for(int j = 0; j < 6; j++) {
-           for(int k = 0; k < 16; k++) {
-               for(int l = 0; l < 4; l++) {
-                   cout << expint21[i][j][k][l] << " ";
-               }
-               cout << endl;
-           }
-           cout << endl;
-       }
-       cout << endl;
-    }
-
-    
-    cout << "expint22[6][6][16][16]:" << endl;
-    for(int i = 0; i < 6; i++) {
-       for(int j = 0; j < 6; j++) {
-           for(int k = 0; k < 16; k++) {
-               for(int l = 0; l < 16; l++) {
-                   cout << expint22[i][j][k][l] << " ";
-               }
-               cout << endl;
-           }
-           cout << endl;
-       }
-       cout << endl;
-    }
-
-}
-
-}