diff --git a/gtfold-mfe/include/energy.h b/gtfold-mfe/include/energy.h index c5fdb1f..5b7b5cc 100644 --- a/gtfold-mfe/include/energy.h +++ b/gtfold-mfe/include/energy.h @@ -4,24 +4,20 @@ #include "data.h" extern int *V; -extern int *VV; -extern int *VV1; extern int *W; extern int *VBI; extern int *VM; extern int **WM; -extern int *WMu; -extern int *WMl; +extern int **WMPrime; extern int *indx; #define V(i,j) V[indx[j]+i] #define VM(i,j) VM[indx[j]+i] #define WM(i,j) WM[i][j] +#define WMPrime(i,j) WMPrime[i][j] #define WMU(i,j) WM[i][j] #define WML(i,j) WM[j][i] -//#define WMU(i,j) WMu[indx[j]+i] -//#define WML(i,j) WMl[indx[j]+i] #define VBI(i,j) VBI[indx[j]+i] #define auPen(i, j) ((( (i)==BASE_U || (j)==BASE_U ) && ( (i)==BASE_A || (i)==BASE_G || (j)==BASE_A || (j)==BASE_G )) ? auend : 0) diff --git a/gtfold-mfe/src/algorithms.c b/gtfold-mfe/src/algorithms.c index 3262322..0fd0c72 100644 --- a/gtfold-mfe/src/algorithms.c +++ b/gtfold-mfe/src/algorithms.c @@ -40,19 +40,19 @@ int calculate(int len, int nThreads) { fprintf(stdout,"Thread count: %3d \n",omp_get_num_threads()); #endif - for (b = TURN+1; b <= len-1; b++) { #ifdef _OPENMP #pragma omp parallel for private (i,j) schedule(guided) #endif for (i = 1; i <= len - b; i++) { j = i + b; - int flag = 0, newWM = INFINITY_; + int newWM = INFINITY_; if (canPair(RNA[i], RNA[j])) { - flag = 1; int eh = canHairpin(i,j)?eH(i,j):INFINITY_; //hair pin - int es = canStack(i,j)?eS(i,j)+getShapeEnergy(i)+getShapeEnergy(j)+V(i+1,j-1):INFINITY_; // stack - if (j-i > 6) { // Internal Loop BEGIN + + int es = canStack(i,j)?eS(i,j)+getShapeEnergy(i)+getShapeEnergy(j)+V(i+1,j-1):INFINITY_; // stack + + if (j-i > 6) { // Internal Loop BEGIN int p=0, q=0; int VBIij = INFINITY_; for (p = i+1; p <= MIN(j-2-TURN,i+MAXLOOP+1) ; p++) { @@ -67,18 +67,14 @@ int calculate(int len, int nThreads) { } VBI(i,j) = VBIij; V(i,j) = V(i,j) + getShapeEnergy(i) + getShapeEnergy(j); - } // Internal Loop END - if (j-i > 10) { // Multi Loop BEGIN - int h; - int VMij, VMijd, VMidj, VMidjd; - VMij = VMijd = VMidj = VMidjd = INFINITY_; - for (h = i+TURN+1; h <= j-1-TURN; h++) { - VMij = MIN(VMij, WMU(i+1,h-1) + WML(h,j-1)); - VMidj = MIN(VMidj, WMU(i+2,h-1) + WML(h,j-1)); - VMijd = MIN(VMijd, WMU(i+1,h-1) + WML(h,j-2)); - VMidjd = MIN(VMidjd, WMU(i+2,h-1) + WML(h,j-2)); - } + + if (j-i > 10) { // Multi Loop BEGIN + int VMij = WMPrime[i+1][j-1]; + int VMidj = WMPrime[i+2][j-1]; + int VMijd = WMPrime[i+1][j-2]; + int VMidjd = WMPrime[i+2][j-2]; + int d3 = canSS(j-1)?Ed3(i,j,j-1):INFINITY_; int d5 = canSS(i+1)?Ed5(i,j,i+1):INFINITY_; VMij = MIN(VMij, (VMidj + d5 +Ec)) ; @@ -87,16 +83,21 @@ int calculate(int len, int nThreads) { VMij = VMij + Ea + Eb + auPenalty(i,j); VM(i,j) = canStack(i,j)?VMij:INFINITY_; } // Multi Loop END - V(i,j) = MIN4(eh,es,VBI(i,j),VM(i,j)); + + V(i,j) = MIN4(eh,es,VBI(i,j),VM(i,j)); } else V(i,j) = INFINITY_; - if (j-i > 4) { // WM BEGIN + + if (j-i > 4) { // WM BEGIN int h; for (h = i+TURN+1 ; h <= j-TURN-1; h++) { //ZS: This sum corresponds to when i,j are NOT paired with each other. //So we need to make sure only terms where i,j aren't pairing are considered. - newWM = (!forcePair(i,j))?MIN(newWM, WMU(i,h-1) + WML(h,j)):newWM; + // Added auxillary storage WMPrime to speedup multiloop calculations + WMPrime[i][j] = MIN(WMPrime[i][j], WMU(i,h-1) + WML(h,j)); + //newWM = (!forcePair(i,j))?MIN(newWM, WMU(i,h-1) + WML(h,j)):newWM; } + newWM = (!forcePair(i,j))?MIN(newWM, WMPrime[i][j]):newWM; newWM = MIN(V(i,j) + auPenalty(i,j) + Eb, newWM); newWM = canSS(i)?MIN(V(i+1,j) + Ed3(j,i+1,i) + auPenalty(i+1,j) + Eb + Ec, newWM):newWM; //i dangle newWM = canSS(j)?MIN(V(i,j-1) + Ed5(j-1,i,j) + auPenalty(i,j-1) + Eb + Ec, newWM):newWM; //j dangle @@ -107,6 +108,7 @@ int calculate(int len, int nThreads) { } // WM END } } + for (j = TURN+2; j <= len; j++) { int i, Wj, Widjd, Wijd, Widj, Wij, Wim1; Wj = INFINITY_; @@ -123,4 +125,3 @@ int calculate(int len, int nThreads) { } return W[len]; } - diff --git a/gtfold-mfe/src/energy.c b/gtfold-mfe/src/energy.c index 1118794..9450f92 100644 --- a/gtfold-mfe/src/energy.c +++ b/gtfold-mfe/src/energy.c @@ -9,12 +9,11 @@ #include "shapereader.h" int *V; -int *VV; -int *VV1; int *W; int *VBI; int *VM; int **WM; +int **WMPrime; int *indx; int alloc_flag = 0; @@ -26,18 +25,6 @@ void create_tables(int len) { exit(-1); } - VV1 = (int *) malloc((len+ 1)*sizeof(int)); - if (VV1 == NULL) { - perror("Cannot allocate variable 'V'"); - exit(-1); - } - - VV = (int *) malloc((len+ 1)*sizeof(int)); - if (VV == NULL) { - perror("Cannot allocate variable 'V'"); - exit(-1); - } - int i; WM = (int **) malloc((len+1)* sizeof(int *)); if (WM == NULL) { @@ -50,7 +37,20 @@ void create_tables(int len) { perror("Cannot allocate variable 'WM[i]'"); exit(-1); } + } + + WMPrime = (int **) malloc((len+1)* sizeof(int *)); + if (WMPrime == NULL) { + perror("Cannot allocate variable 'WM'"); + exit(-1); } + for (i = 0; i <= len; i++) { + WMPrime[i] = (int *)malloc((len+1)* sizeof(int)); + if (WMPrime[i] == NULL) { + perror("Cannot allocate variable 'WM[i]'"); + exit(-1); + } + } VM = (int *) malloc(((len+1)*len/2 + 1) * sizeof(int)); if (VM == NULL) { @@ -87,15 +87,13 @@ void init_tables(int len) { for (i = 0; i <= len; i++) { W[i] = INFINITY_; - VV[i] = INFINITY_; - VV1[i] = INFINITY_; - for (j = 0; j <= len; j++) + for (j = 0; j <= len; j++) { WM[i][j] = INFINITY_; + WMPrime[i][j] = INFINITY_; + } } - LLL = (len)*(len+1)/2 + 1; - for (i = 0; i < LLL; i++) { V[i] = INFINITY_; VM[i] = INFINITY_; @@ -116,11 +114,12 @@ void free_tables(int len) { for (i = 0; i <= len; i++) free(WM[i]); free(WM); + for (i = 0; i <= len; i++) free(WMPrime[i]); + free(WMPrime); + free(VM); free(VBI); free(V); - free(VV); - free(VV1); free(W); } }