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

   1
   2 #ifndef MCCAKILL_H
   3 #define MCCAKILL_H
   4
   5 //#define NDEBUG
   6
   7 #include <string>
   8 #include <iostream>
   9 #include "nrutil.h"
  10 #include <cassert>
  11 #include "Util.hpp"
  12 #include "Beta.hpp"
  13
  14 #include "energy_param.hpp"
  15 //#include "energy_param3.hpp"
  16
  17 #include "ScoreType.hpp"
  18
  19
  20 using namespace std;
  21 using namespace ProbCons;
  22
  23 namespace MXSCARNA {
  24 class McCaskill {
  25     char *seq;
  26     int  *numSeq;
  27     int  n_seq;
  28     static energy_param  energyParam;
  29
  30     Trimat<float> a, q1, ab, am, am1, p;
  31     TriVertMat<float> q1v, abv, amv, am1v, pv;
  32     Trimat<int> typeMat;
  33
  34     void printExpEnergy();
  35     void initParameter();
  36     void Inside();
  37     void Outside();
  38     void convertProbability();
  39
  40     inline float expHairpinEnergy(const int type, const int l, const int i, const int j);
  41     inline float expLoopEnergy(int u1, int u2, int type, int type2,
  42                          int si1, int sj1, int sp1, int sq1);
  43
  44     inline float compQb(int i, int j, int tmpType);
  45     inline float compQm1(int i, int j, int tmpType);
  46     inline float compQm(int i, int j);
  47     inline float compQ1(int i, int j, int tmpType);
  48     inline float compQ(int i, int j);
  49     inline float compP(int h, int l, int tmpType);
  50     inline float compPm(int i, int l);
  51     inline float compPm1(int i, int l);
  52     inline float beginStemScore(const int i, const int j) const;
  53     inline float endStemScore(const int i, const int j) const;
  54
  55     static const float GASCONST;
  56     static const float T;
  57     static const int    MAXLOOP;
  58     static const int    TETRA_ENTH37;
  59     static const int    NBPAIRS;
  60     static const int    SCALE;
  61     static const int    TURN;
  62     static float *exphairpin;
  63     static float expninio[5][32];
  64     static float expdangle5[6][4];
  65     static float expdangle3[6][4];
  66     static float expinternalLoop[31];
  67     static float expbulge[31];
  68     static char   exptriLoop[2][6];
  69     static float exptriLoopEnergy[2];
  70     static char   exptetraLoop[30][7];
  71     static float exptetraLoopEnergy[30];
  72     static float expStack[6][6];
  73     static float expTstackH[6][16];
  74     static float expTstackI[6][16];
  75     static float expint11[6][6][16];
  76     static float expint21[6][6][16][4];
  77     static float expint22[6][6][16][16];
  78     static float expMLclosing;
  79     static float expMLintern[8];
  80     static float expTermAU;
  81     static float expMLbase[31];
  82
  83  public:
  84
  85     McCaskill(int n, const char *mySeq) {
  86         seq    = new char[n + 1];
  87         numSeq = new int[n + 1];
  88         n_seq = 0;
  89
  90
  91         for(int i = 0; i < n; i++) {
  92             if     (mySeq[i] == 'a' || mySeq[i] == 'A') { seq[n_seq] = 'A'; numSeq[n_seq] = Beta::A_CODE; n_seq++; }
  93             else if(mySeq[i] == 't' || mySeq[i] == 'T' ||
  94                     mySeq[i] == 'u' || mySeq[i] == 'U') { seq[n_seq] = 'U'; numSeq[n_seq] = Beta::U_CODE; n_seq++; }
  95             else if(mySeq[i] == 'g' || mySeq[i] == 'G') { seq[n_seq] = 'G'; numSeq[n_seq] = Beta::G_CODE; n_seq++; }
  96             else if(mySeq[i] == 'c' || mySeq[i] == 'C') { seq[n_seq] = 'C'; numSeq[n_seq] = Beta::C_CODE; n_seq++; }
  97             else if(mySeq[i] == 'n' || mySeq[i] == 'N') { seq[n_seq] = 'N'; numSeq[n_seq] = Beta::N_CODE; n_seq++; }
  98             else if(mySeq[i] == '.' || mySeq[i] == '-') { seq[n_seq] = '-'; numSeq[n_seq] = Beta::GAP_CODE; n_seq++; }
  99             else { seq[n_seq] = mySeq[i]; numSeq[n_seq] = Beta::INVALID_CODE; n_seq++; }
 100         }
 101         seq[n_seq] = '\0';
 102         a.Allocator(n_seq);
 103         q1.Allocator(n_seq);
 104         ab.Allocator(n_seq);
 105         am.Allocator(n_seq);
 106         am1.Allocator(n_seq);
 107         p.Allocator(n_seq);
 108         q1v.Allocator(n_seq);
 109         abv.Allocator(n_seq);
 110         amv.Allocator(n_seq);
 111         am1v.Allocator(n_seq);
 112         pv.Allocator(n_seq);
 113         typeMat.Allocator(n_seq);
 114
 115         if(n_seq > 31) {
 116           exphairpin = new float[n_seq + 1];
 117         }
 118         else {
 119           exphairpin = new float[31];
 120         }
 121
 122         for(int i = 0; i < n_seq; i++) {
 123             for(int j = i; j < n_seq; j++) {
 124                 a.ref(i,j)  = q1.ref(i,j) = IMPOSSIBLE;
 125                 q1v.ref(i,j) = IMPOSSIBLE;
 126             }
 127         }
 128
 129         for(int i = 0; i < n_seq; i++) {
 130             a.ref(i,i)   = 0.0;
 131             q1.ref(i,i)  = IMPOSSIBLE;
 132             q1v.ref(i,i) = IMPOSSIBLE;
 133         }
 134
 135         for(int i = 0; i < n_seq-1; i++) {
 136             a.ref(i,i+1)   = 0.0;
 137             q1.ref(i,i+1)  = IMPOSSIBLE;
 138             q1v.ref(i,i+1) = IMPOSSIBLE;
 139         }
 140
 141         for(int i = 0; i < n_seq-2; i++) {
 142             a.ref(i,i+2)   = 0.0;
 143             q1.ref(i,i+2)  = IMPOSSIBLE;
 144             q1v.ref(i,i+2) = IMPOSSIBLE;
 145         }
 146
 147         for(int i = 0; i < n_seq-3; i++) {
 148             a.ref(i,i+3)   = 0.0;
 149             q1.ref(i,i+3)  = IMPOSSIBLE;
 150             q1v.ref(i,i+3) = IMPOSSIBLE;
 151
 152         }
 153
 154         for(int i = 0; i < n_seq; i++) {
 155             for(int j = i; j < n_seq; j++) {
 156                 ab.ref(i,j)  = am.ref(i,j)  = am1.ref(i,j)  = p.ref(i,j)  = IMPOSSIBLE;
 157                 abv.ref(i,j) = amv.ref(i,j) = am1v.ref(i,j) = pv.ref(i,j) = IMPOSSIBLE;
 158             }
 159         }
 160
 161         /* the type of base pair */
 162         /* C <-> G : type 1      */
 163         /* G <-> C : type 2      */
 164         /* G <-> U : type 3      */
 165         /* U <-> G : type 5      */
 166         /* A <-> U : type 0      */
 167         /* U <-> A : type 4      */
 168         /* ? <-> ? : type 6      */
 169         for(int i = 0; i < n_seq; i++) {
 170             for(int j = i; j < n_seq; j++) {
 171                 typeMat.ref(i,j) = Beta::getReducedPairCode(numSeq[i], numSeq[j]);
 172             }
 173         }
 174
 175     }
 176
 177     /*------------------------------------------------------------------------*/
 178     /* dangling ends should never be destabilizing, i.e. expdangle>=1         */
 179     /* specific heat needs smooth function (2nd derivative)                   */
 180     /* we use a*(sin(x+b)+1)^2, with a=2/(3*sqrt(3)), b=Pi/6-sqrt(3)/2,       */
 181     /* in the interval b<x<sqrt(3)/2                                          */
 182     float SMOOTH(float X) {
 183       return ((X)/SCALE<-1.2283697)?0:(((X)/SCALE>0.8660254)?(X):
 184                                   SCALE*0.38490018*(sin((X)/SCALE-0.34242663)+1)*(sin((X)/SCALE-0.34242663)+1));
 185     }
 186
 187     ~McCaskill() {
 188         delete[] seq;
 189         delete[] numSeq;
 190         delete[] exphairpin;
 191     }
 192
 193     void calcPartitionFunction();
 194     void printProbMat();
 195
 196     inline float getProb(const int i, const int j) const {
 197         // 0 origin : 0..(n-1)
 198         return p.ref(i, j);
 199     }
 200 };
 201 }
 202 #endif // MCCASKILL_H