+++ /dev/null
-
-#ifndef MCCAKILL_H
-#define MCCAKILL_H
-
-//#define NDEBUG
-
-#include <string>
-#include <iostream>
-#include "nrutil.h"
-#include <cassert>
-#include "Util.hpp"
-#include "Beta.hpp"
-
-#include "energy_param.hpp"
-//#include "energy_param3.hpp"
-
-#include "ScoreType.hpp"
-
-
-using namespace std;
-using namespace ProbCons;
-
-namespace MXSCARNA {
-class McCaskill {
- char *seq;
- int *numSeq;
- int n_seq;
- static energy_param energyParam;
-
- Trimat<float> a, q1, ab, am, am1, p;
- TriVertMat<float> q1v, abv, amv, am1v, pv;
- Trimat<int> typeMat;
-
- void printExpEnergy();
- void initParameter();
- void Inside();
- void Outside();
- void convertProbability();
-
- inline float expHairpinEnergy(const int type, const int l, const int i, const int j);
- inline float expLoopEnergy(int u1, int u2, int type, int type2,
- int si1, int sj1, int sp1, int sq1);
-
- inline float compQb(int i, int j, int tmpType);
- inline float compQm1(int i, int j, int tmpType);
- inline float compQm(int i, int j);
- inline float compQ1(int i, int j, int tmpType);
- inline float compQ(int i, int j);
- inline float compP(int h, int l, int tmpType);
- inline float compPm(int i, int l);
- inline float compPm1(int i, int l);
- inline float beginStemScore(const int i, const int j) const;
- inline float endStemScore(const int i, const int j) const;
-
- static const float GASCONST;
- static const float T;
- static const int MAXLOOP;
- static const int TETRA_ENTH37;
- static const int NBPAIRS;
- static const int SCALE;
- static const int TURN;
- static float *exphairpin;
- static float expninio[5][32];
- static float expdangle5[6][4];
- static float expdangle3[6][4];
- static float expinternalLoop[31];
- static float expbulge[31];
- static char exptriLoop[2][6];
- static float exptriLoopEnergy[2];
- static char exptetraLoop[30][7];
- static float exptetraLoopEnergy[30];
- static float expStack[6][6];
- static float expTstackH[6][16];
- static float expTstackI[6][16];
- static float expint11[6][6][16];
- static float expint21[6][6][16][4];
- static float expint22[6][6][16][16];
- static float expMLclosing;
- static float expMLintern[8];
- static float expTermAU;
- static float expMLbase[31];
-
- public:
-
- McCaskill(int n, const char *mySeq) {
- seq = new char[n + 1];
- numSeq = new int[n + 1];
- n_seq = 0;
-
-
- for(int i = 0; i < n; i++) {
- if (mySeq[i] == 'a' || mySeq[i] == 'A') { seq[n_seq] = 'A'; numSeq[n_seq] = Beta::A_CODE; n_seq++; }
- else if(mySeq[i] == 't' || mySeq[i] == 'T' ||
- mySeq[i] == 'u' || mySeq[i] == 'U') { seq[n_seq] = 'U'; numSeq[n_seq] = Beta::U_CODE; n_seq++; }
- else if(mySeq[i] == 'g' || mySeq[i] == 'G') { seq[n_seq] = 'G'; numSeq[n_seq] = Beta::G_CODE; n_seq++; }
- else if(mySeq[i] == 'c' || mySeq[i] == 'C') { seq[n_seq] = 'C'; numSeq[n_seq] = Beta::C_CODE; n_seq++; }
- else if(mySeq[i] == 'n' || mySeq[i] == 'N') { seq[n_seq] = 'N'; numSeq[n_seq] = Beta::N_CODE; n_seq++; }
- else if(mySeq[i] == '.' || mySeq[i] == '-') { seq[n_seq] = '-'; numSeq[n_seq] = Beta::GAP_CODE; n_seq++; }
- else { seq[n_seq] = mySeq[i]; numSeq[n_seq] = Beta::INVALID_CODE; n_seq++; }
- }
- seq[n_seq] = '\0';
- a.Allocator(n_seq);
- q1.Allocator(n_seq);
- ab.Allocator(n_seq);
- am.Allocator(n_seq);
- am1.Allocator(n_seq);
- p.Allocator(n_seq);
- q1v.Allocator(n_seq);
- abv.Allocator(n_seq);
- amv.Allocator(n_seq);
- am1v.Allocator(n_seq);
- pv.Allocator(n_seq);
- typeMat.Allocator(n_seq);
-
- if(n_seq > 31) {
- exphairpin = new float[n_seq + 1];
- }
- else {
- exphairpin = new float[31];
- }
-
- for(int i = 0; i < n_seq; i++) {
- for(int j = i; j < n_seq; j++) {
- a.ref(i,j) = q1.ref(i,j) = IMPOSSIBLE;
- q1v.ref(i,j) = IMPOSSIBLE;
- }
- }
-
- for(int i = 0; i < n_seq; i++) {
- a.ref(i,i) = 0.0;
- q1.ref(i,i) = IMPOSSIBLE;
- q1v.ref(i,i) = IMPOSSIBLE;
- }
-
- for(int i = 0; i < n_seq-1; i++) {
- a.ref(i,i+1) = 0.0;
- q1.ref(i,i+1) = IMPOSSIBLE;
- q1v.ref(i,i+1) = IMPOSSIBLE;
- }
-
- for(int i = 0; i < n_seq-2; i++) {
- a.ref(i,i+2) = 0.0;
- q1.ref(i,i+2) = IMPOSSIBLE;
- q1v.ref(i,i+2) = IMPOSSIBLE;
- }
-
- for(int i = 0; i < n_seq-3; i++) {
- a.ref(i,i+3) = 0.0;
- q1.ref(i,i+3) = IMPOSSIBLE;
- q1v.ref(i,i+3) = IMPOSSIBLE;
-
- }
-
- for(int i = 0; i < n_seq; i++) {
- for(int j = i; j < n_seq; j++) {
- ab.ref(i,j) = am.ref(i,j) = am1.ref(i,j) = p.ref(i,j) = IMPOSSIBLE;
- abv.ref(i,j) = amv.ref(i,j) = am1v.ref(i,j) = pv.ref(i,j) = IMPOSSIBLE;
- }
- }
-
- /* the type of base pair */
- /* C <-> G : type 1 */
- /* G <-> C : type 2 */
- /* G <-> U : type 3 */
- /* U <-> G : type 5 */
- /* A <-> U : type 0 */
- /* U <-> A : type 4 */
- /* ? <-> ? : type 6 */
- for(int i = 0; i < n_seq; i++) {
- for(int j = i; j < n_seq; j++) {
- typeMat.ref(i,j) = Beta::getReducedPairCode(numSeq[i], numSeq[j]);
- }
- }
-
- }
-
- /*------------------------------------------------------------------------*/
- /* dangling ends should never be destabilizing, i.e. expdangle>=1 */
- /* specific heat needs smooth function (2nd derivative) */
- /* we use a*(sin(x+b)+1)^2, with a=2/(3*sqrt(3)), b=Pi/6-sqrt(3)/2, */
- /* in the interval b<x<sqrt(3)/2 */
- float SMOOTH(float X) {
- return ((X)/SCALE<-1.2283697)?0:(((X)/SCALE>0.8660254)?(X):
- SCALE*0.38490018*(sin((X)/SCALE-0.34242663)+1)*(sin((X)/SCALE-0.34242663)+1));
- }
-
- ~McCaskill() {
- delete[] seq;
- delete[] numSeq;
- delete[] exphairpin;
- }
-
- void calcPartitionFunction();
- void printProbMat();
-
- inline float getProb(const int i, const int j) const {
- // 0 origin : 0..(n-1)
- return p.ref(i, j);
- }
-};
-}
-#endif // MCCASKILL_H
git://source.jalview.org / jabaws.git / blobdiff