binaries/src/mafft/extensions/mxscarna_src/McCaskill.cpp

   1 #include <iostream>
   2 #include "McCaskill.hpp"
   3 //#include "energy_param3.hpp"
   4 #include "Util.hpp"
   5 #include <cstring>
   6
   7 namespace MXSCARNA {
   8 energy_param  McCaskill::energyParam;
   9
  10 float *McCaskill::exphairpin;
  11 float McCaskill::expninio[5][32];
  12 float McCaskill::expdangle5[6][4];
  13 float McCaskill::expdangle3[6][4];
  14 float McCaskill::expinternalLoop[31];
  15 float McCaskill::expbulge[31];
  16 char   McCaskill::exptriLoop[2][6];
  17 float McCaskill::exptriLoopEnergy[2];
  18 char   McCaskill::exptetraLoop[30][7];
  19 float McCaskill::exptetraLoopEnergy[30];
  20 float McCaskill::expStack[6][6];
  21 float McCaskill::expTstackH[6][16];
  22 float McCaskill::expTstackI[6][16];
  23 float McCaskill::expint11[6][6][16];
  24 float McCaskill::expint21[6][6][16][4];
  25 float McCaskill::expint22[6][6][16][16];
  26 float McCaskill::expMLclosing;
  27 float McCaskill::expMLintern[8];
  28 float McCaskill::expTermAU;
  29 float McCaskill::expMLbase[31];
  30
  31 const int McCaskill::TURN         = 3;
  32 const float McCaskill::GASCONST   = 1.98717;
  33 const float McCaskill::T          = 37 + 273.15;
  34 const int McCaskill::MAXLOOP      = 30;
  35 const int McCaskill::TETRA_ENTH37 = -400;
  36 const int McCaskill::NBPAIRS      = 7;
  37 const int McCaskill::SCALE        = 10;
  38
  39
  40 void
  41 McCaskill::
  42 calcPartitionFunction()
  43 {
  44     initParameter();
  45     Inside();
  46     Outside();
  47     convertProbability();
  48
  49 /*
  50     for(int i = 0; i < n_seq; i++) {
  51         for(int j = i; j < n_seq; j++) {
  52             cout << getProb(i, j) << " ";
  53         }
  54         cout << endl;
  55     }
  56 */
  57 }
  58
  59 void
  60 McCaskill::
  61 convertProbability()
  62 {
  63     float *pPointer  = p.getPointer(0, 0);
  64     float *abPointer = ab.getPointer(0,0);
  65     for(int i = 0; i < n_seq; i++) {
  66         for(int j = i; j < n_seq; j++) {
  67             *pPointer += *abPointer++;
  68             *pPointer++ = EXP(*pPointer);
  69         }
  70     }
  71 }
  72
  73 void
  74 McCaskill::
  75 Inside()
  76 {
  77
  78     for (int i = n_seq - TURN - 2; i >= 0; i--) {
  79         float *abPointer   = ab.getPointer(i, i + TURN + 1);
  80         float *am1Pointer  = am1.getPointer(i, i + TURN + 1);
  81         float *amPointer   = am.getPointer(i, i + TURN + 1);
  82         float *q1Pointer   = q1.getPointer(i, i + TURN + 1);
  83         float *aPointer    = a.getPointer(i, i + TURN + 1);
  84         int   *typePointer = typeMat.getPointer(i, i+TURN+1);
  85         for (int j = i + TURN + 1; j < n_seq; j++) {
  86             int tmpType   = *typePointer++;
  87                             *abPointer++  = compQb(i, j, tmpType);
  88             am1v.ref(i,j) = *am1Pointer++ = compQm1(i,j, tmpType);
  89             amv.ref(i,j)  = *amPointer++  = compQm(i,j);
  90             q1v.ref(i,j)  = *q1Pointer++  = compQ1(i,j, tmpType);
  91                             *aPointer++   = compQ(i,j);
  92         }
  93     }
  94 }
  95
  96 inline float
  97 McCaskill::
  98 compQb(int i, int j, int tmpType)
  99 {
 100
 101     float tmpAb;
 102     int type  = tmpType;
 103     int u     = j - i - 1;
 104
 105     if(Beta::isCanonicalReducedPairCode(type)) {
 106         // hairpin energy
 107         assert(u >= 3);
 108         tmpAb = expHairpinEnergy(type, u, i + 1, j - 1);
 109
 110         // base pair, bulge, interior loop energy
 111         for(int h = i + 1; h <= MIN(i + MAXLOOP + 1, j - TURN - 2); h++) {
 112             int u1 = h-i-1;
 113             int max = MAX(h + TURN + 1, j-1-MAXLOOP+u1);
 114             float *abPointer     = ab.getPointer(h, max - 1);
 115             const int *typeMatPointer = typeMat.getPointer(h, max);
 116
 117             for(int l = max; l < j; l++) {
 118                 int type2 = *typeMatPointer++;
 119                 abPointer++;
 120                 if(!Beta::isCanonicalReducedPairCode(type2)) continue;
 121
 122                 assert(h >= 0 && h < n_seq && l >= 0 && l < n_seq);
 123                 type2 = Beta::flipReducedPairCode(type2);
 124                 assert(h-i-1 >= 0); assert(j-l-1 >= 0);
 125                 float loopE = *abPointer;
 126                 loopE += expLoopEnergy(u1, j-l-1, tmpType, type2, i+1, j-1, h-1, l+1);
 127                 tmpAb = LOG_ADD(tmpAb, loopE);
 128             }
 129         }
 130
 131         // multi loop
 132         float tmpQm = IMPOSSIBLE;
 133         float *amPointer  = am.getPointer(i + 1, j - TURN - 3);
 134         float *am1Pointer = am1v.getPointer(j-TURN-2, j-1);
 135         for(int h = j - TURN - 2; h >= i + TURN + 3; h--) {
 136             assert(h >= 0 && h < n_seq);
 137             float tmpScore = *amPointer--;
 138             tmpScore += *am1Pointer--;
 139             tmpQm = LOG_ADD(tmpQm, tmpScore);
 140         }
 141         tmpQm += expMLclosing + expMLintern[type];
 142         tmpQm += endStemScore(i, j);
 143         tmpAb = LOG_ADD(tmpAb, tmpQm);
 144     }
 145     else {
 146         tmpAb = IMPOSSIBLE;
 147     }
 148     return tmpAb;
 149 }
 150
 151 //F = a:ML_closing + b:ML_intern*k + c:ML_BASE*u
 152
 153 inline float
 154 McCaskill::
 155 compQm1(int i, int j, int tmpType)
 156 {
 157     float tmpQm1 = IMPOSSIBLE;
 158
 159     int l = j;
 160     if (i + TURN + 1 <= l) {
 161         int type = typeMat.ref(i,l);
 162         if(Beta::isCanonicalReducedPairCode(type)) {
 163             float tmpScore = ab.ref(i,l);
 164             tmpScore += beginStemScore(i, l);
 165             //tmpScore += expMLbase[1]*(j-l) + expMLintern[type];
 166             tmpScore += expMLintern[tmpType];
 167             tmpQm1 = LOG_ADD(tmpQm1, tmpScore);
 168         }
 169     }
 170     if(i < j) {
 171         tmpQm1 = LOG_ADD(tmpQm1, am1.ref(i,j-1));
 172     }
 173
 174     return tmpQm1;
 175 }
 176
 177 inline float
 178 McCaskill::
 179 compQm(int i, int j)
 180 {
 181     float tmpQm = IMPOSSIBLE;
 182     float *amPointer  = am.getPointer(i,j-TURN-2);
 183     float *am1Pointer = am1v.getPointer(j-TURN-1, j);
 184     for(int h = j - TURN - 1; h >= i ; h--) {
 185         float tmpScore = 0;
 186         float tmpAm1 = *am1Pointer--;
 187
 188         tmpScore += tmpAm1;
 189         tmpQm = LOG_ADD(tmpQm, tmpScore);
 190         tmpScore = *amPointer--;
 191         tmpScore += tmpAm1;
 192         tmpQm = LOG_ADD(tmpQm, tmpScore);
 193     }
 194
 195     return tmpQm;
 196 }
 197
 198 inline float
 199 McCaskill::
 200 compQ1(int i, int j, int tmpType)
 201 {
 202     float tmpQ1 = IMPOSSIBLE;
 203
 204     if(Beta::isCanonicalReducedPairCode(tmpType)) {
 205         float tmpScore = ab.ref(i, j);
 206         tmpScore += beginStemScore(i, j);
 207         tmpQ1 = LOG_ADD(tmpQ1, tmpScore);
 208     }
 209     tmpQ1 = LOG_ADD(tmpQ1, q1.ref(i, j - 1));
 210
 211     return tmpQ1;
 212 }
 213
 214 inline float
 215 McCaskill::
 216 compQ(int i, int j)
 217 {
 218     float tmpQ = 0;
 219     tmpQ = LOG_ADD(tmpQ, q1.ref(i,j));
 220
 221     float *aPointer  = a.getPointer(i,j-TURN-2);
 222     float *q1Pointer = q1v.getPointer(j-TURN-1, j);
 223     for(int h = j - TURN - 1; h >= i + 1; h--) {
 224         float tmpScore  = *aPointer--;
 225         tmpScore       += *q1Pointer--;
 226         tmpQ            = LOG_ADD(tmpQ, tmpScore);
 227     }
 228
 229     return tmpQ;
 230 }
 231
 232 inline float
 233 McCaskill::
 234 beginStemScore(const int i, const int j) const
 235 {
 236     float temp = 0;
 237     int type = typeMat.ref(i,j);
 238
 239     if(0 < i)                   { temp += expdangle5[type][numSeq[i-1]]; }
 240     if(j < n_seq-1)             { temp += expdangle3[type][numSeq[j+1]]; }
 241     if(type != Beta::REDUCED_CG_CODE && type != Beta::REDUCED_GC_CODE)  { temp += expTermAU; }
 242     return temp;
 243 }
 244
 245 inline float
 246 McCaskill::
 247 endStemScore(const int i, const int j) const
 248 {
 249     float temp = 0;
 250     int type = typeMat.ref(i,j);
 251
 252     type = Beta::flipReducedPairCode(type);
 253
 254     if(i < n_seq-1)             { temp += expdangle3[type][numSeq[i+1]]; }
 255     if(j > 0)                   { temp += expdangle5[type][numSeq[j-1]]; }
 256     if(type != Beta::REDUCED_CG_CODE && type != Beta::REDUCED_GC_CODE)  { temp += expTermAU; }
 257     return temp;
 258 }
 259
 260 inline float
 261 McCaskill::
 262 compP(int h, int l, int tmpType)
 263 {
 264     float prob = IMPOSSIBLE;
 265
 266     int type = tmpType;
 267     if(Beta::isCanonicalReducedPairCode(type)) {
 268
 269         /* exterior loop */
 270         float tmp_p = 0;
 271         tmp_p -= a.ref(0,n_seq-1);
 272         if(0 < h) {
 273             tmp_p += a.ref(0,h-1);
 274         }
 275         if(l < n_seq-1) {
 276             tmp_p += a.ref(l+1, n_seq-1);
 277         }
 278         tmp_p += beginStemScore(h, l);
 279         prob = LOG_ADD(prob, tmp_p);
 280
 281         assert(IMPOSSIBLE <= prob && prob <= 0);
 282
 283         /* internal loop */
 284         tmp_p = IMPOSSIBLE;
 285         int tt = Beta::flipReducedPairCode(tmpType);
 286         int max = MAX(0,h-MAXLOOP-1);
 287         for(int i = max; i <= h - 1; i++) {
 288             float min = MIN(l+MAXLOOP-h+i+2, n_seq-1);
 289             int   *typeMatPointer = typeMat.getPointer(i,l+1);
 290             float *pPointer       = p.getPointer(i,l);
 291             for(int j = l + 1; j <= min; j++) {
 292                 int type2    = *typeMatPointer++;
 293                 pPointer++;
 294                 if(!Beta::isCanonicalReducedPairCode(type2)) continue;
 295                 assert(i >= 0 && i < n_seq && j >= 0 && j < n_seq);
 296
 297                 float tmpScore  = *pPointer;
 298                 tmpScore       += expLoopEnergy(h-i-1, j-l-1, type2, tt, i+1, j-1, h-1, l+1);
 299                 tmp_p = LOG_ADD(tmp_p, tmpScore);
 300             }
 301         }
 302         prob = LOG_ADD(prob, tmp_p);
 303         assert(IMPOSSIBLE <= prob && prob <= 0);
 304
 305         /* multi loop */
 306         tmp_p = IMPOSSIBLE;
 307         float tmp_begin   = beginStemScore(h, l);
 308         float *q1Pointer  = q1v.getPointer(0, l);
 309         float *am1Pointer = am1v.getPointer(0, l);
 310         float *amPointer  = amv.getPointer(1,h-1);
 311         for(int i = 0; i <= h-TURN-1; i++) {
 312             float tmpq1    = *q1Pointer++;
 313             float tmpam    = *amPointer++;
 314             float tmpScore = *am1Pointer++;
 315
 316             tmpScore += tmpam;
 317             tmpScore += tmp_begin;
 318             tmpScore += expMLclosing + expMLintern[tt];
 319             tmp_p = LOG_ADD(tmp_p, tmpScore);
 320
 321             tmpScore  = tmpq1;
 322             tmpScore += tmpam;
 323             tmpScore += tmp_begin;
 324             tmpScore += expMLclosing + expMLintern[tt];
 325             tmp_p = LOG_ADD(tmp_p, tmpScore);
 326
 327             tmpScore  = tmpq1;
 328             tmpScore += tmp_begin;
 329             tmpScore += expMLclosing + expMLintern[tt];
 330             tmp_p = LOG_ADD(tmp_p, tmpScore);
 331         }
 332
 333         assert(IMPOSSIBLE <= tmp_p && tmp_p <= 0);
 334         prob = LOG_ADD(prob, tmp_p);
 335
 336         tmp_p = IMPOSSIBLE;
 337         for(int i = h-TURN; i <= h-1; i++) {
 338             if(i >= 0) {
 339                 float tmpScore  = q1.ref(i,l);
 340                 tmpScore       += tmp_begin;
 341                 tmpScore       += expMLclosing + expMLintern[tt];
 342                 tmp_p           = LOG_ADD(tmp_p, tmpScore);
 343             }
 344         }
 345         assert(IMPOSSIBLE <= tmp_p && tmp_p <= 0);
 346         prob = LOG_ADD(prob, tmp_p);
 347     }
 348     else {
 349         prob = IMPOSSIBLE;
 350     }
 351
 352     return prob;
 353 }
 354
 355 inline float
 356 McCaskill::
 357 compPm(int i, int l)
 358 {
 359   float tmpPm  = IMPOSSIBLE;
 360
 361   int *typeMatPointer   = typeMat.getPointer(i,n_seq-1);
 362   float *pPointer       = p.getPointer(i,n_seq);
 363   float *amPointer      = am.getPointer(l+1,n_seq-1);
 364   float *abPointer      = ab.getPointer(i, n_seq);
 365   for(int j = n_seq - 1; j >= l + TURN + 1; j--) {
 366       int type = *typeMatPointer--;
 367       pPointer--;
 368       amPointer--;
 369       abPointer--;
 370       if(Beta::isCanonicalReducedPairCode(type)) {
 371           float tmp  = *pPointer;
 372           tmp       += *amPointer;
 373           tmp       += endStemScore(i, j);
 374           tmpPm = LOG_ADD(tmpPm, tmp);
 375       }
 376   }
 377   tmpPm += expMLintern[1];
 378
 379   return tmpPm;
 380 }
 381
 382 inline float
 383 McCaskill::
 384 compPm1(int i, int l)
 385 {
 386   float tmpPm1 = IMPOSSIBLE;
 387
 388   int j =  l + 1;
 389   if(j <= n_seq-1) {
 390       int type = typeMat.ref(i,j);
 391       if(Beta::isCanonicalReducedPairCode(type)) {
 392           float tmp = p.ref(i,j);
 393           tmp += endStemScore(i, j);
 394           tmpPm1 = tmp;
 395       }
 396       tmpPm1 += expMLintern[1];
 397   }
 398   if(l+1 <= n_seq - 1) {
 399       tmpPm1 = LOG_ADD(tmpPm1, am1.ref(i, l+1));
 400   }
 401
 402   return tmpPm1;
 403 }
 404
 405 void
 406 McCaskill::
 407 Outside()
 408 {
 409     for(int h = 0; h <= n_seq - TURN - 2; h++) {
 410         float *pPointer    = p.getPointer(h, n_seq-1);
 411         float *q1Pointer   = q1.getPointer(h, n_seq-1);
 412         float *am1Pointer  = am1.getPointer(h, n_seq-1);
 413         int   *typePointer = typeMat.getPointer(h, n_seq-1);
 414         for(int l = n_seq-1; l >= h + TURN + 1; l--) {
 415             int tmpType    = *typePointer--;
 416             pv.ref(h,l)    = *pPointer--   = compP(h,l,tmpType);
 417             q1v.ref(h,l)   = *q1Pointer--  = compPm(h,l);
 418             am1v.ref(h,l)  = *am1Pointer-- = compPm1(h,l);
 419
 420             assert(p.ref(h,l) <= 0);
 421         }
 422     }
 423 }
 424
 425 void
 426 McCaskill::
 427 printProbMat()
 428 {
 429     int m = 0;
 430     for(int i = 0; i < n_seq; i++) cout << " " << seq[i];
 431     cout << endl;
 432     for(int i = 0; i < n_seq; i++) {
 433         if(m < n_seq) {
 434             cout << seq[m];
 435         }
 436         for(int j = 0; j <= i-1; j++) {
 437             if(j != i-1) cout << "  ";
 438             else         cout << " ";
 439         }
 440         if(i != 0 && i != n_seq-1) {
 441             cout << "\\";
 442         }
 443
 444         for(int j = i; j < n_seq; j++) {
 445             if(p.ref(i,j) > 0.01) {
 446
 447                 int type = Beta::getReducedPairCode(numSeq[i], numSeq[j]);
 448
 449                 if(!Beta::isCanonicalReducedPairCode(type)) {
 450                     cout << "\n" << seq[i] << " " << seq[j] << " " << exp(p.ref(i,j)) << endl;
 451                 }
 452
 453                 if(j != n_seq-1) {
 454                     cout << "* ";
 455                 }
 456                 else {
 457                     cout << "*";
 458                 }
 459
 460             }
 461             else {
 462
 463                 if(j != n_seq-1) {
 464                     cout << "  ";
 465                 }
 466                 else {
 467                     cout << " ";
 468                 }
 469
 470             }
 471
 472         }
 473         if(m < n_seq) {
 474             cout << seq[m++] << endl;
 475         }
 476         if(i == n_seq - 1) cout << endl;
 477     }
 478     for(int i = 0; i < n_seq; i++) cout << " " << seq[i];
 479     cout << endl;
 480 }
 481
 482 void
 483 McCaskill::
 484 initParameter()
 485 {
 486     float GT;
 487     float RT_KCAL_MOL = McCaskill::T*McCaskill::GASCONST;
 488     int len = 31;
 489
 490     for(int i = 0; i < len; i++) {
 491         GT = energyParam.getHairpin(i);
 492         exphairpin[i] = -GT*10/RT_KCAL_MOL;
 493     }
 494
 495     for (int i = 0; i < len; i++) {
 496         GT = energyParam.getBulge(i);
 497         expbulge[i] = -GT*10/RT_KCAL_MOL;
 498         GT = energyParam.getInternalLoop(i);
 499         expinternalLoop[i] = -GT*10/RT_KCAL_MOL;
 500     }
 501     expinternalLoop[2] = -80*10/RT_KCAL_MOL; /* special case of size 2 interior loops (single mismatch) */
 502
 503     // float lxc = energy_param3::lxc37;
 504     for (int i = 31; i < n_seq; i++) {
 505         GT = energyParam.getHairpin(30) + (107.856*LOG((float)i/30));
 506         exphairpin[i] = -GT*10/RT_KCAL_MOL;
 507     }
 508
 509     for(int i = 0; i < 5; i++) {
 510         GT = energyParam.getNinio(i);
 511         for(int j = 0; j <= MAXLOOP; j++) {
 512             expninio[i][j] = -MIN(energyParam.getMaxNinio(), j*GT)*10/RT_KCAL_MOL;
 513         }
 514     }
 515
 516     for(int i = 0; i < 30; i++) {
 517       GT = energyParam.getTetraLoopEnergy(i);
 518       exptetraLoopEnergy[i] = -GT*10/RT_KCAL_MOL;
 519     }
 520
 521     /*no parameters for Triloop*/
 522     for(int i = 0; i < 2; i++) {
 523         GT = 0;
 524         exptriLoopEnergy[i] = -GT*10/RT_KCAL_MOL;
 525     }
 526
 527     GT = energyParam.getMLclosing();
 528     expMLclosing = -GT*10/RT_KCAL_MOL;
 529
 530     for(int i = 0; i <= NBPAIRS; i++) {
 531       GT = energyParam.getMLintern();
 532       expMLintern[i] = -GT*10/RT_KCAL_MOL;
 533     }
 534
 535     expTermAU = -energyParam.getTerminalAU()*10/RT_KCAL_MOL;
 536
 537     GT = energyParam.getMLBase();
 538     for(int i = 0; i < len; i++) {
 539         expMLbase[i] = -GT*10*(float)i/RT_KCAL_MOL;
 540     }
 541
 542     /*
 543        if danlges = 0 just set their energy to 0,
 544        don't let dangle energyies become > 0 (at large temps)
 545        but make sure go smoothly to 0
 546     */
 547     for(int i = 0; i < 6; i++) {
 548         for(int j =0; j < 4; j++) {
 549             GT = energyParam.getDangle5(i,j);
 550             expdangle5[i][j] = -GT*10/RT_KCAL_MOL;
 551             GT = energyParam.getDangle3(i,j);
 552             expdangle3[i][j] = -GT*10/RT_KCAL_MOL;
 553         }
 554     }
 555
 556     /* stacking energies  */
 557     for(int i = 0; i < 6; i++) {
 558         for(int j = 0; j < 6; j++) {
 559             GT = energyParam.getStack(i,j);
 560             expStack[i][j] = -GT*10/RT_KCAL_MOL;
 561         }
 562     }
 563
 564     /* mismatch energies */
 565     for (int i = 0; i < 6; i++) {
 566         for (int j = 0; j < 16; j++) {
 567           GT = energyParam.getTstackI(i, j);
 568           //  cout << i << " " << " " << j << " " << GT << endl;
 569           expTstackI[i][j] = -GT*10/RT_KCAL_MOL;
 570           GT = energyParam.getTstackH(i, j);
 571           expTstackH[i][j] = -GT*10/RT_KCAL_MOL;
 572         }
 573     }
 574
 575     /* interior loops of length 2*/
 576     for(int i = 0; i < 6; i++) {
 577         for(int j = 0; j < 6; j++) {
 578             for(int k = 0; k < 16; k++) {
 579               GT = energyParam.getInt11(i, j, k);
 580               expint11[i][j][k] = -GT*10/RT_KCAL_MOL;
 581             }
 582         }
 583     }
 584
 585     /* interior 2*1 loops */
 586     for(int i = 0; i < 6; i++) {
 587         for(int j =0; j < 6; j++) {
 588             for(int k = 0; k < 16; k++) {
 589                 for(int l = 0; l < 4; l++) {
 590                   GT = energyParam.getInt21(i,j,k,l);
 591                   expint21[i][j][k][l] = -GT*10/RT_KCAL_MOL;
 592                 }
 593             }
 594         }
 595     }
 596
 597     /* interior 2*2 loops */
 598     for (int i = 0; i < 6; i++) {
 599         for(int j = 0; j < 6; j++) {
 600             for(int k = 0; k < 16; k++) {
 601                 for(int l = 0; l < 16; l++) {
 602                   GT = energyParam.getInt22(i,j,k,l);
 603                   expint22[i][j][k][l] = -GT*10/RT_KCAL_MOL;
 604                 }
 605             }
 606         }
 607     }
 608 }
 609
 610
 611 inline float
 612 McCaskill::
 613 expHairpinEnergy(const int type, const int l, const int i, const int j)
 614 {
 615   float q;
 616   int    k;
 617
 618 //  assert(l >= 0);
 619   q = exphairpin[l];
 620
 621   if(l == 4) {
 622       char temp_seq[7];
 623
 624       for(int iter = i - 1; iter < i + 5; iter++) {
 625         temp_seq[iter - (i-1)] = seq[iter];
 626       }
 627       temp_seq[6] = '\0';
 628
 629       for(k = 0; k < 30; k++) {
 630         if(strcmp(temp_seq, energyParam.getTetraLoop(k)) == 0) break;
 631       }
 632       if(k != 30) {
 633         q += exptetraLoopEnergy[k];
 634       }
 635   }
 636   if(l == 3) {
 637
 638       /* no triloop bonus
 639     char temp_seq[6];
 640
 641     for(int iter = i - 1; iter < i + 4; iter++) {
 642         temp_seq[iter - (i-1)] = seq[iter];
 643     }
 644     temp_seq[6] = '\0';
 645     for(k = 0; k < 2; k++) {
 646       if(strcmp(temp_seq, energyParam.getTriLoop(k)) == 0) break;
 647     }
 648     if(k != 2) {
 649       q *= exptriLoopEnergy[k];
 650     }
 651       */
 652
 653     if(type != Beta::REDUCED_CG_CODE && type != Beta::REDUCED_GC_CODE) q += expTermAU;
 654   }
 655   else {
 656       int type2 = Beta::getPairCode(numSeq[i], numSeq[j]);
 657       q += expTstackH[type][type2];
 658   }
 659
 660   return q;
 661 }
 662
 663
 664 inline float
 665 McCaskill::
 666 expLoopEnergy(int u1, int u2, int type, int type2,
 667               int si1, int sj1, int sp1, int sq1)
 668 {
 669   float z = 0;
 670
 671   if((u1 == 0) && (u2 == 0)) { z = expStack[type][type2]; }
 672   else {
 673     if((u1 == 0) || (u2 == 0)) {
 674       int u;
 675       if(u1 == 0) { u = u2; }
 676       else        { u = u1; }
 677       z = expbulge[u];
 678       if(u1 + u2 == 1) z += expStack[type][type2];
 679       else {
 680         if (type  != Beta::REDUCED_CG_CODE && type  != Beta::REDUCED_GC_CODE) z += expTermAU;
 681         if (type2 != Beta::REDUCED_CG_CODE && type2 != Beta::REDUCED_GC_CODE) z += expTermAU;
 682       }
 683     }
 684     else {
 685       if(u1 + u2 == 2) {
 686           z = expint11[type][type2][Beta::getPairCode(numSeq[si1], numSeq[sj1])];
 687       }
 688       else if((u1 == 1) && (u2 == 2)) {
 689           z = expint21[type][type2][Beta::getPairCode(numSeq[si1], numSeq[sj1])][numSeq[sq1]];
 690       }
 691       else if((u1 == 2) && (u2 == 1)) {
 692           z = expint21[type2][type][Beta::getPairCode(numSeq[sq1], numSeq[sp1])][numSeq[si1]];
 693       }
 694       else if((u1 == 2) && (u2 == 2)) {
 695         z = expint22[type][type2][Beta::getPairCode(numSeq[si1], numSeq[sj1])][Beta::getPairCode(numSeq[sp1], numSeq[sq1])];
 696       }
 697       else {
 698         z = expinternalLoop[u1 + u2] +
 699             expTstackI[type][Beta::getPairCode(numSeq[si1], numSeq[sj1])]
 700              + expTstackI[type2][Beta::getPairCode(numSeq[sq1], numSeq[sp1])];
 701         z += expninio[2][abs(u1-u2)];
 702       }
 703     }
 704   }
 705
 706   return z;
 707 }
 708
 709 void
 710 McCaskill::
 711 printExpEnergy()
 712 {
 713     cout << "exphairpin:" << endl;
 714     for(int i = 0; i < 31; i++) {
 715         cout << exphairpin[i] << endl;
 716     }
 717
 718     cout << "expninio[5][32]:" << endl;
 719     for(int i = 0; i < 5; i++) {
 720         for(int j = 0; j < 32; j++) {
 721             cout << expninio[i][j] << " ";
 722         }
 723         cout << endl;
 724     }
 725
 726     cout << "expdangle5[6][4]:" << endl;
 727     for(int i = 0; i < 6; i++) {
 728         for(int j = 0; j < 4; j++) {
 729             cout << expdangle5[i][j] << " ";
 730         }
 731         cout << endl;
 732     }
 733
 734     cout << "expdangle3[6][4]:" << endl;
 735     for(int i = 0; i < 6; i++) {
 736         for(int j = 0; j < 4; j++) {
 737             cout << expdangle3[i][j] << " ";
 738         }
 739         cout << endl;
 740     }
 741
 742     cout << "expinternalLoop[31]:" << endl;
 743     for(int i = 0; i < 31; i++) {
 744         cout << i << ":" << expinternalLoop[i] << endl;
 745     }
 746     cout << "expbulge[31]:" << endl;
 747     for(int i = 0; i < 31; i++) {
 748         cout << i << ":" << expbulge[i] << endl;
 749     }
 750
 751     cout << "exptriLoopEnergy[2]:" << endl;
 752     for(int i = 0; i < 2; i++) {
 753         cout << i << ":" << exptriLoopEnergy[i] << endl;
 754     }
 755
 756     cout << "exptetraLoopEnergy[15]" << endl;
 757     for(int i = 0; i < 15; i++) {
 758         cout << i << ":" << exptetraLoopEnergy[i] << endl;
 759     }
 760
 761     cout << "expStack[6][6]:" << endl;
 762     for(int i = 0; i < 6; i++) {
 763         for(int j = 0; j < 6; j++) {
 764             cout << expStack[i][j] << " ";
 765         }
 766         cout << endl;
 767     }
 768
 769     cout << "expTstackH[6][16]:" << endl;
 770     for(int i = 0; i < 6; i++) {
 771         for(int j = 0; j < 16; j++) {
 772             cout << expTstackH[i][j] << " ";
 773         }
 774         cout << endl;
 775     }
 776
 777     cout << "expTstackI[6][16]:" << endl;
 778     for(int i = 0; i < 6; i++) {
 779         for(int j = 0; j < 16; j++) {
 780             cout << expTstackI[i][j] << " ";
 781         }
 782         cout << endl;
 783     }
 784
 785     cout << "expMLclosing=" << expMLclosing << endl;
 786     cout << "expMLintern:" << endl;
 787     for(int i = 0; i < 8; i++) {
 788         cout << expMLintern[i] << " ";
 789     }
 790     cout << endl;
 791
 792     cout << "expMLbase[31]:";
 793     for(int i = 0; i < 31; i++) {
 794         cout << i << ":" << expMLbase[i] << endl;
 795     }
 796
 797     cout << "expint11[6][6][16]:";
 798     for(int i = 0; i < 6; i++) {
 799         for(int j = 0; j < 6; j++) {
 800             for(int k = 0; k < 16; k++) {
 801                 cout << expint11[i][j][k] << " ";
 802             }
 803             cout << endl;
 804         }
 805         cout << endl;
 806     }
 807
 808     cout << "expint21[6][6][16][4]:" << endl;
 809     for(int i = 0; i < 6; i++) {
 810         for(int j = 0; j < 6; j++) {
 811             for(int k = 0; k < 16; k++) {
 812                 for(int l = 0; l < 4; l++) {
 813                     cout << expint21[i][j][k][l] << " ";
 814                 }
 815                 cout << endl;
 816             }
 817             cout << endl;
 818         }
 819         cout << endl;
 820     }
 821
 822
 823     cout << "expint22[6][6][16][16]:" << endl;
 824     for(int i = 0; i < 6; i++) {
 825         for(int j = 0; j < 6; j++) {
 826             for(int k = 0; k < 16; k++) {
 827                 for(int l = 0; l < 16; l++) {
 828                     cout << expint22[i][j][k][l] << " ";
 829                 }
 830                 cout << endl;
 831             }
 832             cout << endl;
 833         }
 834         cout << endl;
 835     }
 836
 837 }
 838
 839 }