14 #include "energy_param.hpp"
15 //#include "energy_param3.hpp"
17 #include "ScoreType.hpp"
21 using namespace ProbCons;
28 static energy_param energyParam;
30 Trimat<float> a, q1, ab, am, am1, p;
31 TriVertMat<float> q1v, abv, amv, am1v, pv;
34 void printExpEnergy();
38 void convertProbability();
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);
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;
55 static const float GASCONST;
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];
85 McCaskill(int n, const char *mySeq) {
86 seq = new char[n + 1];
87 numSeq = new int[n + 1];
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++; }
106 am1.Allocator(n_seq);
108 q1v.Allocator(n_seq);
109 abv.Allocator(n_seq);
110 amv.Allocator(n_seq);
111 am1v.Allocator(n_seq);
113 typeMat.Allocator(n_seq);
116 exphairpin = new float[n_seq + 1];
119 exphairpin = new float[31];
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;
129 for(int i = 0; i < n_seq; i++) {
131 q1.ref(i,i) = IMPOSSIBLE;
132 q1v.ref(i,i) = IMPOSSIBLE;
135 for(int i = 0; i < n_seq-1; i++) {
137 q1.ref(i,i+1) = IMPOSSIBLE;
138 q1v.ref(i,i+1) = IMPOSSIBLE;
141 for(int i = 0; i < n_seq-2; i++) {
143 q1.ref(i,i+2) = IMPOSSIBLE;
144 q1v.ref(i,i+2) = IMPOSSIBLE;
147 for(int i = 0; i < n_seq-3; i++) {
149 q1.ref(i,i+3) = IMPOSSIBLE;
150 q1v.ref(i,i+3) = IMPOSSIBLE;
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;
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]);
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));
193 void calcPartitionFunction();
196 inline float getProb(const int i, const int j) const {
197 // 0 origin : 0..(n-1)
202 #endif // MCCASKILL_H