Merge remote branch 'prashant/rewrite'

gtfold-mfe/include/constraints.h

@@ -18,9 +18,15 @@ void print_constraints(int length) ;
 #ifdef __cplusplus
 extern "C" {
 #endif
-int ssOK(int i, int j); 
-int baseOK(int i);
-int pairOK(int i, int j);
+int is_ss(int i, int j); 
+int prohibit_base(int i) ;
+int check_base(int i) ;
+int force_pair(int i, int j) ;
+int force_pair1(int i, int j) ;
+int check_iloop(int i, int j, int p, int q) ;
+int check_pair(int i, int j) ;
+int check_stack(int i, int j) ;
+int check_hairpin(int i, int j) ;
 #ifdef __cplusplus
 }
 #endif

gtfold-mfe/include/options.h

@@ -16,6 +16,7 @@ extern bool LIMIT_DISTANCE;
 extern bool BPP_ENABLED;
 extern bool SUBOPT_ENABLED;
 extern bool CONS_ENABLED;
+extern bool VERBOSE;
 
 extern string seqfile;
 extern string constraintsFile;

gtfold-mfe/include/traceback.h

@@ -26,7 +26,7 @@
 #ifdef __cplusplus
 extern "C" {
 #endif
-	void trace(int len);
+	void trace(int len, int vbose);
 #ifdef __cplusplus
 }
 #endif

gtfold-mfe/src/algorithms.c

@@ -58,10 +58,8 @@ int calculate(int len, int nThreads) {
 
 			if (canPair(RNA[i], RNA[j])) {
 				flag = 1;
-
-				int eh = eH(i, j); //hair pin
-				int es = eS(i, j) +  V(i+1,j-1); // stack
-
+				int eh = check_hairpin(i,j)?INFINITY_:eH(i, j); //hair pin
+				int es = check_stack(i,j)?INFINITY_:(eS(i, j) + V(i+1,j-1)); // stack
 				if (j-i > 6) {  // Internal Loop BEGIN
 					int p=0, q=0;
 					int VBIij = INFINITY_;
@@ -71,12 +69,11 @@ int calculate(int len, int nThreads) {
 						if (minq < p+1+TURN) minq = p+1+TURN;
 						for (q = minq; q < j; q++) {
 							if (!canPair(RNA[p], RNA[q])) continue;
-							if (!ssOK(i,p) || !ssOK(q,j)) continue;
-
+							if (check_iloop(i,j,p,q)) continue;
 							VBIij = MIN(eL(i, j, p, q) + V(p,q), VBIij);
 						}
 					}
-					VBI(i,j) = pairOK(i,j)?VBIij:INFINITY_;
+					VBI(i,j) = check_pair(i,j)?INFINITY_:VBIij;
 				} 	// Internal Loop END
 
 				if (j-i > 10) {	 // Multi Loop BEGIN
@@ -95,57 +92,63 @@ int calculate(int len, int nThreads) {
 					int d3 = Ed3(i,j,j-1);
 					int d5 = Ed5(i,j,i+1);
 
-					VMij = MIN(VMij, baseOK(i+1)?(VMidj + d5 +Ec):INFINITY_);
-					VMij = MIN(VMij, baseOK(j-1)?(VMijd + d3 +Ec):INFINITY_);
-					VMij = MIN(VMij, (baseOK(i+1)&&baseOK(j-1))?(VMidjd + d5 +  d3+ 2*Ec):INFINITY_);
+					VMij = MIN(VMij, check_base(i+1)?(VMidj + d5 +Ec):INFINITY_) ;
+					VMij = MIN(VMij, check_base(j-1)?(VMijd + d3 +Ec):INFINITY_);
+					VMij = MIN(VMij, (check_base(i+1)&&check_base(j-1))?(VMidjd + d5 +  d3+ 2*Ec):INFINITY_);
 					VMij = VMij + Ea + Eb + auPenalty(i,j);
-
-					VM(i,j) = pairOK(i,j)?VMij:INFINITY_;
+					VM(i,j) = check_pair(i,j)?INFINITY_:VMij;
 				} // Multi Loop END
 
-				V(i,j) = pairOK(i,j)?MIN4(eh,es,VBI(i,j),VM(i,j)):INFINITY_;
+				V(i,j) = check_pair(i,j)?INFINITY_: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; 
 				if (!flag) {
 					for (h = i+TURN+1 ; h <= j-TURN-1; h++) {
 						newWM = MIN(newWM, WMU(i,h-1) + WML(h,j));
 					}
 				}
 
-				newWM = MIN(V(i,j) + auPenalty(i,j) + Eb, newWM);  //V(i,j)
-				newWM = baseOK(i)?MIN(V(i+1,j) + Ed3(j,i+1,i) + auPenalty(i+1,j) + Eb + Ec, newWM): INFINITY_; //V(i+1,j) 
-				newWM = baseOK(j)?MIN(V(i,j-1) + Ed5(j-1,i,j) + auPenalty(i,j-1) + Eb + Ec, newWM): INFINITY_; //V(i,j-1)
-				newWM = (baseOK(i)&&baseOK(j))?MIN(V(i+1,j-1) + Ed3(j-1,i+1,i) + Ed5(j-1,i+1,j) + auPenalty(i+1,j-1) + Eb + 2*Ec, newWM): INFINITY_ ;
+				newWM = MIN(V(i,j) + auPenalty(i,j) + Eb, newWM); 
+				newWM = check_base(i)?MIN(V(i+1,j) + Ed3(j,i+1,i) + auPenalty(i+1,j) + Eb + Ec, newWM): INFINITY_; 
+				newWM = check_base(j)?MIN(V(i,j-1) + Ed5(j-1,i,j) + auPenalty(i,j-1) + Eb + Ec, newWM) : INFINITY_; 
+				newWM = (check_base(i)&&check_base(j))?MIN(V(i+1,j-1) + Ed3(j-1,i+1,i) + Ed5(j-1,i+1,j) + auPenalty(i+1,j-1) + Eb + 2*Ec, newWM): INFINITY_;
 
-				newWM = baseOK(i)?MIN(WMU(i+1,j) + Ec, newWM):INFINITY_;
-				newWM = baseOK(j)?MIN(WML(i,j-1) + Ec, newWM):INFINITY_;
+				newWM = check_base(i)?MIN(WMU(i+1,j) + Ec, newWM):INFINITY_;
+				newWM = check_base(j)?MIN(WML(i,j-1) + Ec, newWM):INFINITY_;
 
 				WMU(i,j) = WML(i,j) = newWM;
-				// WM END
-			}
+			} // WM END
 		}
-
 	}
-	
+
 	for (j = TURN+2; j <= len; j++)	{
-		int i, Wj, Widjd, Wijd, Widj, Wij, Wim1;
+		int i, branch=0, Wj, Widjd, Wijd, Widj, Wij, Wim1;
 		Wj = INFINITY_;
 		for (i = 1; i < j-TURN; i++) {
 			Wij = Widjd = Wijd = Widj = INFINITY_;
 			Wim1 = MIN(0, W[i-1]); 
 			Wij = V(i, j) + auPenalty(i, j) + Wim1;
-			Widjd = V(i+1,j-1) + auPenalty(i+1,j-1) + Ed3(j-1,i + 1,i) + Ed5(j-1,i+1,j) + Wim1;
-			Wijd = V(i,j-1) + auPenalty(i,j- 1) + Ed5(j-1,i,j) + Wim1;
-			Widj = V(i+1, j) + auPenalty(i+1,j) + Ed3(j,i + 1,i) + Wim1;
-			Wj = MIN(MIN(MIN(Wij, Widjd), MIN(Wijd, Widj)), Wj); 
+			Widjd = (check_base(i)&&check_base(j))?(V(i+1,j-1) + auPenalty(i+1,j-1) + Ed3(j-1,i + 1,i) + Ed5(j-1,i+1,j) + Wim1):INFINITY_;
+			Wijd = check_base(j)?(V(i,j-1) + auPenalty(i,j-1) + Ed5(j-1,i,j) + Wim1):INFINITY_;
+			Widj =  (check_base(i))?(V(i+1, j) + auPenalty(i+1,j) + Ed3(j,i + 1,i) + Wim1):INFINITY_;
+			Wj = MIN(MIN4(Wij, Widjd, Wijd, Widj), Wj); 
+
+			if (Wj<INFINITY_) {
+				if (Wj==Wij && force_pair1(i,j))
+				   branch = 1;
+				else if	(Wj==Widjd && force_pair1(i+1,j-1))
+				   branch = 1;
+				else if	(Wj==Wijd && force_pair1(i,j-1))
+				  branch = 1;	
+				else if	(force_pair1(i+1,j))
+					branch = 1;
+			}
 		}
-
-		W[j] = MIN(Wj, W[j-1]);
+		W[j] = branch?Wj:MIN(Wj, W[j-1]);
 	}
 
 	return W[len];

gtfold-mfe/src/constraints.cc

@@ -3,6 +3,8 @@
 #include <stdio.h>
 #include <iostream>
 #include <fstream>
+#include <vector>
+#include <algorithm>
 
 #include "global.h"
 #include "options.h"
@@ -15,10 +17,15 @@ int** FBP;
 int nPBP;
 int nFBP;
 
+bool compare_bp(const std::pair<int,int>& o1, 
+			   	const std::pair<int,int>& o2) {
+	return o1.first < o2.first;
+}
+
+
 static int load_constraints(const char* constr_file, int verbose=0) {
 	std::ifstream cfcons;
     fprintf(stdout, "- Running with constraints\n");
-    //fprintf(stdout, "Opening constraint file: %s\n", constr_file);
 
     cfcons.open(constr_file, std::ios::in);
     if (cfcons == 0) {
@@ -77,20 +84,21 @@ static int load_constraints(const char* constr_file, int verbose=0) {
         }
     }
 
-	if (verbose == 1) {
-		fprintf(stdout, "Forced base pairs: ");
-		for(it=0; it< nFBP; it++) {
-			for(int k=1;k<= FBP[it][2];k++)
-				fprintf(stdout, "(%d,%d) ", FBP[it][0]+k-1, FBP[it][1]-k+1);
-		}
-		fprintf(stdout, "\nProhibited base pairs: ");
-		for(it=0; it< nPBP; it++) {
-			for(int k=1;k<= PBP[it][2];k++)
-				fprintf(stdout, "(%d,%d) ", PBP[it][0]+k-1, PBP[it][1]-k+1);
+	std::vector<std::pair<int,int> > v_fbp;
+	for(it=0; it< nFBP; it++) {
+		for(int k=1;k<= FBP[it][2];k++)
+			v_fbp.push_back(std::pair<int,int>(FBP[it][0]+k-1, FBP[it][1]-k+1));
+	}
+	std::sort(v_fbp.begin(), v_fbp.end(), compare_bp);
+	for (size_t ii = 0; ii < v_fbp.size() -1 ; ++ii) {
+		if (v_fbp[ii].second <= v_fbp[ii+1].second) {
+			fprintf(stderr, "\nConstraints create pseudoknots, exiting !!!\n");
+			exit(-1);
 		}
-		fprintf(stdout, "\n\n");
+			
 	}
 
+
     return 0;
 }
 
@@ -110,6 +118,10 @@ int init_constraints(const char* constr_file,int length) {
 
 	if (nPBP != 0) {
 		for (it = 0; it < nPBP; it++) {
+			if (PBP[it][2] < 1) {
+				printf("Invalid entry (%d %d %d)\n", PBP[it][0], PBP[it][1],PBP[it][2]);
+				continue;
+			}	
 			for(k= 1; k <= PBP[it][2];k++) {
 				BP[PBP[it][0]+k-1] = -1;
 				if(PBP[it][1]!=0) {
@@ -120,11 +132,21 @@ int init_constraints(const char* constr_file,int length) {
 	}
 	if (nFBP != 0) {
 		for (it = 0; it < nFBP; it++) {
+			if (FBP[it][2] < 1) {
+				printf("Invalid entry (%d %d %d)\n", FBP[it][0], FBP[it][1], FBP[it][2]);
+				continue;
+			}
 			for(k=1; k <= FBP[it][2];k++) {
+				int i1 = FBP[it][0]+k-1;
+				int j1 = FBP[it][1]-k+1;
 				if (!canPair(RNA[FBP[it][0]+k-1], RNA[FBP[it][1]-k+1])) {
 					printf("Can't constrain (%d,%d)\n", FBP[it][0]+k-1, FBP[it][1]-k+1);
 					continue;
 				}
+				if (j1-i1 < TURN) {
+					printf("Can't constrain (%d,%d)\n", i1, j1);
+					continue;
+				}	
 				BP[FBP[it][0]+k-1] = FBP[it][1]+1-k;
 				BP[FBP[it][1]+1-k] = FBP[it][0]+k-1;
 			}
@@ -153,32 +175,64 @@ void print_constraints(int len) {
     printf("\n");
 }
 
-int ssOK(int i, int j) {
+int is_ss(int i, int j) {
 	if (CONS_ENABLED) {
 		int it;
 		for (it = i + 1; it < j; it++) {
-			if (BP[it] > 0) return 0;
+			if (BP[it] > 0) return 1;
 		}
-		return 1;
+		return 0;
 	}
 	else
-		return 1;
+		return 0;
+}
+
+int prohibit_base(int i) {
+		return (BP[i] == -1);
 }
 
-int baseOK(int i) {
-	if (CONS_ENABLED)
-		return (BP[i] <=0);
+int check_base(int i) {
+	if (CONS_ENABLED) 
+		return (BP[i] <= 0);
 	else
 		return 1;
 }
 
-int pairOK(int i, int j) {
-	if (CONS_ENABLED)
-	return !(BP[i] < 0 || BP[j] < 0 ||  
-			(BP[i] > 0 && j != BP[i]) ||
-			(BP[j] > 0 && i != BP[j]));  			
+int force_pair(int i, int j) {
+		return (BP[i] > 0 && j != BP[i]);
+}
+
+int force_pair1(int i, int j) {
+	if (CONS_ENABLED) 
+		return (BP[i]==j);
 	else
-		return 1;
+		return 0;
+}
+
+int check_iloop(int i, int j, int p, int q) {
+	if (CONS_ENABLED) 
+		return is_ss(i,p) || is_ss(q,j);
+	else 
+		return 0;
+}
+
+int check_pair(int i, int j) {
+	if (CONS_ENABLED) 
+		return prohibit_base(i) || prohibit_base(j) || force_pair(i,j) || force_pair(j,i);
+	else
+		return 0;
 }
 
+int check_stack(int i, int j) {
+	if (CONS_ENABLED) 
+		return force_pair(i,j) || force_pair(j,i);
+	else
+		return 0;
+}
 
+int check_hairpin(int i, int j) {
+	if (CONS_ENABLED) 
+		return is_ss(i,j) || force_pair(i,j) || force_pair(j,i);
+	else
+		return 0;
+}

gtfold-mfe/src/main.cc

@@ -220,7 +220,7 @@ int main(int argc, char** argv) {
 	}
 
 	t1 = get_seconds();
-	trace(seq.length());
+	trace(seq.length(), VERBOSE);
 	t1 = get_seconds() - t1;
 
 	printf("\n");

gtfold-mfe/src/options.cc

@@ -11,6 +11,7 @@ bool LIMIT_DISTANCE;
 bool BPP_ENABLED;
 bool SUBOPT_ENABLED;
 bool CONS_ENABLED = false;
+bool VERBOSE = false;
 
 string seqfile = "";
 string constraintsFile = "";
@@ -87,7 +88,10 @@ void parse_options(int argc, char** argv) {
 					nThreads = atoi(argv[++i]);
 				else
 					help();	
-			} else if(strcmp(argv[i], "--bpp") == 0) {
+			} else if (strcmp(argv[i], "--verbose") == 0 || strcmp(argv[i], "-v") == 0) {
+                VERBOSE = true;
+			}
+			else if(strcmp(argv[i], "--bpp") == 0) {
 				BPP_ENABLED = true;
 			} else if(strcmp(argv[i], "--subopt") == 0) {
 				SUBOPT_ENABLED = true;
@@ -112,6 +116,12 @@ void parse_options(int argc, char** argv) {
 
 		// base it off the input file
 		outputFile += seqfile;
+
+		size_t pos;
+		// extract file name from the path
+		if ((pos=outputFile.find_last_of('/')) > 0) {
+			outputFile = outputFile.substr(pos+1);
+		}
 
 		// and if an extension exists, remove it ...
 		if(outputFile.find(".") != string::npos)
@@ -132,16 +142,6 @@ void printRunConfiguration(string seq) {
 
 	printf("Run Configuration:\n");
 
-//#ifdef _OPENMP
-//	if(nThreads == -1)
-//		printf("- [OMP] thread count: %d\n", omp_get_num_threads());
-//	else
-//		printf("- [OMP] thread count: %d\n", nThreads);
-//#else
-//	printf("- thread count: 1\n");
-//#endif
-//
-
 	if (NOISOLATE == true) {
 		printf("- preventing isolated base pairs\n");
 		standardRun = false;