source: trunk/GDE/MAFFT/mafft-7.055-with-extensions/extensions/mxscarna_src/McCaskill.cpp

Last change on this file was 10371, checked in by aboeckma, 11 years ago

updated mafft version. Added extensions (no svn ignore, yet)
File size: 19.7 KB
Line 
1#include <iostream>
2#include "McCaskill.hpp"
3//#include "energy_param3.hpp"
4#include "Util.hpp"
5#include <cstring>
6
7namespace MXSCARNA {
8energy_param  McCaskill::energyParam;
9
10float *McCaskill::exphairpin;
11float McCaskill::expninio[5][32];
12float McCaskill::expdangle5[6][4];
13float McCaskill::expdangle3[6][4];
14float McCaskill::expinternalLoop[31];
15float McCaskill::expbulge[31];
16char   McCaskill::exptriLoop[2][6];
17float McCaskill::exptriLoopEnergy[2];
18char   McCaskill::exptetraLoop[30][7];
19float McCaskill::exptetraLoopEnergy[30];
20float McCaskill::expStack[6][6];
21float McCaskill::expTstackH[6][16];
22float McCaskill::expTstackI[6][16];
23float McCaskill::expint11[6][6][16];
24float McCaskill::expint21[6][6][16][4];
25float McCaskill::expint22[6][6][16][16];
26float McCaskill::expMLclosing;
27float McCaskill::expMLintern[8];
28float McCaskill::expTermAU;
29float McCaskill::expMLbase[31];
30
31const int McCaskill::TURN         = 3;
32const float McCaskill::GASCONST   = 1.98717;
33const float McCaskill::T          = 37 + 273.15;
34const int McCaskill::MAXLOOP      = 30;
35const int McCaskill::TETRA_ENTH37 = -400;
36const int McCaskill::NBPAIRS      = 7;
37const int McCaskill::SCALE        = 10;
38
39
40void 
41McCaskill::
42calcPartitionFunction()
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
59void
60McCaskill::
61convertProbability()
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
73void 
74McCaskill::
75Inside()
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
96inline float
97McCaskill::
98compQb(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
153inline float
154McCaskill::
155compQm1(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
177inline float
178McCaskill::
179compQm(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
198inline float
199McCaskill::
200compQ1(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
214inline float
215McCaskill::
216compQ(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
232inline float
233McCaskill::
234beginStemScore(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
245inline float
246McCaskill::
247endStemScore(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
260inline float
261McCaskill::
262compP(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
355inline float
356McCaskill::
357compPm(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
382inline float
383McCaskill::
384compPm1(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
405void 
406McCaskill::
407Outside()
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
425void
426McCaskill::
427printProbMat()
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
482void
483McCaskill::
484initParameter()
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
611inline float 
612McCaskill::
613expHairpinEnergy(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
664inline float 
665McCaskill::
666expLoopEnergy(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
709void
710McCaskill::
711printExpEnergy()
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}
Note: See TracBrowser for help on using the repository browser.