Phase 1: audit and regression-backed correctness fixes

ChangeLog

@@ -12,6 +12,13 @@
 
 # IntaRNA
 
+- BUGFIX : heuristic cell initialization and noLP ensemble recurrences
+- BUGFIX : ensemble filtering and seed-extension partition reset
+- BUGFIX : target accessibility limits and seed-free outMinPu filtering
+- BUGFIX : accessibility range splitting and intramolecular ensemble energies
+- BUGFIX : interaction assignment, seed comparison and output handling
+- reject Zall/Eall-dependent output in window mode
+
 # IntaRNA_plotRegions.R
 
 - fix deprecation warning
@@ -23,6 +30,46 @@
 ################################################################################  
 ################################################################################
 
+260919 Alexander Mitrofanov
+ * IntaRNA/PredictorMfe2dHeuristic :
+ * IntaRNA/PredictorMfe2dHeuristicSeed :
+   * fillHybridE() : BUGFIX : reset best energy for each cell
+ * IntaRNA/PredictorMfeEns2dHeuristic :
+   * fillHybridZ() :
+     * BUGFIX : reset best energy for each cell
+     * BUGFIX : noLP double counting and missing stacked/bulged extensions
+ * IntaRNA/PredictorMfeEns :
+   * updateZ() : BUGFIX : apply noGUend and maxED filters
+ * IntaRNA/PredictorMfeEns2dSeedExtension :
+ * IntaRNA/PredictorMfeEns2dHeuristicSeedExtension :
+   * predict() : BUGFIX : clear partition functions when no seed is found
+ * IntaRNA/PredictorMfeEns2dSeedExtension :
+   * fillHybridZ_left() : BUGFIX : multiply stacked partition factors
+ * bin/CommandLineParsing :
+   * parse() : reject Zall/Eall-dependent output in window mode
+   * getTargetAccessibility() : BUGFIX : use target-specific length limits
+   * getQueryRanges(), getTargetRanges() : BUGFIX : outMinPu without seeds
+ * IntaRNA/Accessibility :
+   * decomposeByMaxED() : BUGFIX : exclude forbidden positions and fix range length
+ * IntaRNA/InteractionEnergyBasePair :
+   * computeES() : BUGFIX : exclude empty structure
+ * IntaRNA/NussinovHandler :
+   * getQb() : BUGFIX : zero weight for out-of-range paired intervals
+ * IntaRNA/InteractionEnergyVrna :
+   * computeES() : BUGFIX : heed maximal base pair span
+   * computeES(), computeIntraEall() : BUGFIX : free ViennaRNA data on exit
+ * IntaRNA/Interaction :
+   * operator=() : BUGFIX : handle self-assignment
+   * operator==() : BUGFIX : handle missing seeds
+   * Seed::operator<() : BUGFIX : compare all seed boundaries on energy ties
+ * IntaRNA/OutputHandlerHub :
+   * add() : BUGFIX : fix output forwarding signature
+   * reported() : BUGFIX : report maximal child count
+ * IntaRNA/OutputHandlerInteractionList :
+   * add() : BUGFIX : handle zero storage capacity
+ + tests : API and CLI regressions
+ + doc/refactor/1-current-state.md : architecture and correctness audit
+
 250219 Martin Raden
  * python/CopomuS.py
  * python/copomus/IntaRNA.py

src/IntaRNA/Accessibility.cpp

@@ -213,8 +213,8 @@ decomposeByMaxED( IndexRangeList & ranges, const E_type maxED, const size_t minR
 		for (size_t i= range->from; i <= range->to; i++) {
 			if (E_isINF(getED(i,i)) || (getED(i,i) > maxED && !E_equal(getED(i,i),maxED))) {
 				// check if end of range found and to be stored
-				if (lastStart < i && minRangeLength <= (i +1 - lastStart)) {
-					out.push_back(IndexRange(lastStart,i));
+				if (lastStart < i && minRangeLength <= (i - lastStart)) {
+					out.push_back(IndexRange(lastStart,i - 1));
 				}
 				lastStart = range->to +1;
 			} else {

src/IntaRNA/Interaction.cpp

@@ -151,6 +151,9 @@ Interaction &
 Interaction::
 operator= ( const Interaction & toCopy )
 {
+	if (this == &toCopy) {
+		return *this;
+	}
 #if INTARNA_IN_DEBUG_MODE
 	if (!toCopy.isValid())
 		throw std::runtime_error("Interaction::=("+toString(toCopy)+") not valid!");
@@ -227,7 +230,8 @@ operator == ( const Interaction &i ) const
 			&& s2 == i.s2
 			&& E_equal( energy, i.energy )
 			&& basePairs == i.basePairs
-			&& (seed == i.seed || *seed == *(i.seed))
+			&& (seed == i.seed
+					|| (seed != NULL && i.seed != NULL && *seed == *(i.seed)))
 			;
 }
 

src/IntaRNA/Interaction.h

@@ -175,16 +175,26 @@ class Interaction {
 		E_type energy;
 
 		/**
-		 * order definition: first by increasing energy using increasing seq1
-		 * index as tie breaker.
+		 * order definition: first by increasing energy, then lexicographically
+		 * by all seed boundaries as tie breakers.
 		 * @param s the seed to compare to
 		 * @return true if this seed is considered smaller than s
 		 */
 		const bool
 		operator <  ( const Seed &s ) const {
-			return ( energy < s.energy
-						|| (E_equal(energy,s.energy) && (bp_i.first < s.bp_i.first))
-						);
+			if (! E_equal(energy, s.energy )) {
+				return energy < s.energy;
+			}
+			if (bp_i.first != s.bp_i.first) {
+				return bp_i.first < s.bp_i.first;
+			}
+			if (bp_i.second != s.bp_i.second) {
+				return bp_i.second < s.bp_i.second;
+			}
+			if (bp_j.first != s.bp_j.first) {
+				return bp_j.first < s.bp_j.first;
+			}
+			return bp_j.second < s.bp_j.second;
 		}
 
 		/**

src/IntaRNA/InteractionEnergyBasePair.cpp

@@ -26,7 +26,10 @@ void InteractionEnergyBasePair::computeES(const RnaSequence &seq,
       if (Z_equal(q_val, 1.0)) {
         logQ(i, j) = E_INF;
       } else {
-        logQ(i, j) = getE(q_val);
+        // getES* covers only structures containing at least one base pair.
+        // The full monomer partition Q also contains the empty structure with
+        // unit weight, which has to be removed here.
+        logQ(i, j) = getE(q_val - Z_type(1.0));
       }
     }
   }

src/IntaRNA/InteractionEnergyVrna.cpp

@@ -3,6 +3,8 @@
 #include "IntaRNA/AccessibilityVrna.h"
 
 #include <cassert>
+#include <cstdlib>
+#include <memory>
 #include <set>
 
 // ES computation
@@ -16,6 +18,26 @@ extern "C" {
 }
 
 
+namespace {
+
+struct VrnaAllocatedDeleter {
+	void operator()( char * data ) const {
+		free(data);
+	}
+};
+
+struct VrnaFoldCompoundDeleter {
+	void operator()( vrna_fold_compound_t * foldCompound ) const {
+		vrna_fold_compound_free(foldCompound);
+	}
+};
+
+typedef std::unique_ptr<char, VrnaAllocatedDeleter> VrnaAllocatedPtr;
+typedef std::unique_ptr<vrna_fold_compound_t, VrnaFoldCompoundDeleter> VrnaFoldCompoundPtr;
+
+}
+
+
 
 namespace IntaRNA {
 
@@ -94,8 +116,10 @@ computeES( const Accessibility & acc, InteractionEnergyVrna::EsMatrix & esToFill
 	const Z_type RT = getRT();
 
 	// VRNA compatible data structures
-	char * sequence = (char *) vrna_alloc(sizeof(char) * (seqLength + 1));
-	char * structureConstraint = (char *) vrna_alloc(sizeof(char) * (seqLength + 1));
+	VrnaAllocatedPtr sequenceOwner( (char *) vrna_alloc(sizeof(char) * (seqLength + 1)) );
+	VrnaAllocatedPtr structureConstraintOwner( (char *) vrna_alloc(sizeof(char) * (seqLength + 1)) );
+	char * const sequence = sequenceOwner.get();
+	char * const structureConstraint = structureConstraintOwner.get();
 	for (int i=0; i<seqLength; i++) {
 		// copy sequence
 		sequence[i] = acc.getSequence().asString().at(i);
@@ -112,7 +136,8 @@ computeES( const Accessibility & acc, InteractionEnergyVrna::EsMatrix & esToFill
 		curModel.max_bp_span = -1;
 	}
 	// TODO check if VRNA_OPTION_WINDOW reasonable to speedup
-	vrna_fold_compound_t * foldData = vrna_fold_compound( sequence, &foldModel, VRNA_OPTION_PF);
+	VrnaFoldCompoundPtr foldDataOwner( vrna_fold_compound( sequence, &curModel, VRNA_OPTION_PF) );
+	vrna_fold_compound_t * const foldData = foldDataOwner.get();
 
 	// Adding hard constraints from pseudo dot-bracket
 	unsigned int constraint_options = VRNA_CONSTRAINT_DB_DEFAULT;
@@ -157,11 +182,6 @@ computeES( const Accessibility & acc, InteractionEnergyVrna::EsMatrix & esToFill
 			}
 		}
 	}
-	// garbage collection
-	vrna_fold_compound_free(foldData);
-	free(structureConstraint);
-	free(sequence);
-
 }
 
 ////////////////////////////////////////////////////////////////////////////
@@ -178,14 +198,16 @@ computeIntraEall( const Accessibility & acc ) const
 	const int length = acc.getSequence().size();
 
 	// copy sequence into C data structure
-	char * sequence = (char *) vrna_alloc(sizeof(char) * (length + 1));
+	VrnaAllocatedPtr sequenceOwner( (char *) vrna_alloc(sizeof(char) * (length + 1)) );
+	char * const sequence = sequenceOwner.get();
 	for (int i=0; i<length; i++) {
 		sequence[i] = acc.getSequence().asString().at(i);
 	}
 	sequence[length] = '\0';
 
-    // setup folding data
-    vrna_fold_compound_t * fold_compound = vrna_fold_compound( sequence, &curModel, VRNA_OPTION_DEFAULT );
+	// setup folding data
+	VrnaFoldCompoundPtr foldCompoundOwner( vrna_fold_compound( sequence, &curModel, VRNA_OPTION_DEFAULT ) );
+	vrna_fold_compound_t * const fold_compound = foldCompoundOwner.get();
 
     // add accessibility constraints
     AccessibilityVrna::addConstraints( *fold_compound, acc );

src/IntaRNA/NussinovHandler.cpp

@@ -94,7 +94,7 @@ NussinovHandler::getQb(const size_t i, const size_t j, const RnaSequence &seq,
     const Z_type bpWeight, const size_t minLoopLength,
     NussinovHandler::Z2dMatrix &Q, NussinovHandler::Z2dMatrix &Qb) {
   if (j >= seq.size()) {
-    return 1.0;
+    return 0.0;
   }
   if (i + minLoopLength >= j) {
     return 0.0;

src/IntaRNA/OutputHandlerHub.h

@@ -204,11 +204,11 @@ addOutputHandler( OutputHandler * handler )
 inline
 void
 OutputHandlerHub::
-add( const Interaction & inter, const OutputConstraint & outConstraint )
+add( const Interaction & inter )
 {
 	// forward to all in list
 	for (auto it=outList.begin(); it!=outList.end(); it++) {
-		(*it)->add(inter,outConstraint);
+		(*it)->add(inter);
 	}
 }
 
@@ -222,7 +222,7 @@ reported() const
 	size_t maxReported = 0;
 	// get maximal number of reports among all handlers
 	for (auto it=outList.begin(); it!=outList.end(); it++) {
-		maxReported = std::min( maxReported, (*it)->reported() );
+		maxReported = std::max( maxReported, (*it)->reported() );
 	}
 	// return maximum
 	return maxReported;

src/IntaRNA/OutputHandlerInteractionList.cpp

@@ -52,7 +52,9 @@ add( const Interaction & interaction )
 	{
 		// count interaction
 		reportedInteractions++;
-		if (storage.size() < maxToStore || lessThan_StorageContainer( &interaction, *(storage.rbegin()) )) {
+		if (maxToStore > 0
+				&& (storage.size() < maxToStore
+						|| lessThan_StorageContainer( &interaction, *(storage.rbegin()) ))) {
 			// find where to insert this interaction
 			StorageContainer::iterator insertPos = std::lower_bound( storage.begin(), storage.end(), &interaction, lessThan_StorageContainer );
 			// check if interaction is NOT already part of the list

src/IntaRNA/PredictorMfe2dHeuristic.cpp

@@ -99,6 +99,7 @@ fillHybridE()
 
 			// init as invalid boundary
 			*curCell = BestInteractionE(E_INF, RnaSequence::lastPos, RnaSequence::lastPos);
+			curCellEtotal = E_INF;
 
 			// check if positions can form interaction
 			if (	energy.isAccessible1(i1)

src/IntaRNA/PredictorMfe2dHeuristicSeed.cpp

@@ -116,6 +116,7 @@ fillHybridE()
 			// init as invalid boundary
 			*curCell = BestInteractionE(E_INF, RnaSequence::lastPos, RnaSequence::lastPos);
 			*curCellSeed = BestInteractionE(E_INF, RnaSequence::lastPos, RnaSequence::lastPos);
+			curCellEtotal = E_INF;
 
 			// check if positions can form interaction
 			if (	energy.isAccessible1(i1)

src/IntaRNA/PredictorMfeEns.cpp

@@ -46,6 +46,22 @@ updateZ( const size_t i1, const size_t j1
 	// check if something to be done
 	if (Z_equal(partZ,0) || Z_isINF(Zall))
 		return;
+
+	// Apply the same site filters used for MFE candidates before changing
+	// either the global or boundary-specific partition.
+	const OutputConstraint & outConstraint = output.getOutputConstraint();
+	if (outConstraint.noGUend
+			&& (energy.isGU(i1,i2) || energy.isGU(j1,j2)))
+	{
+		return;
+	}
+	if (outConstraint.maxED < Accessibility::ED_UPPER_BOUND
+			&& (energy.getED1(i1,j1) > outConstraint.maxED
+					|| energy.getED2(i2,j2) > outConstraint.maxED))
+	{
+		return;
+	}
+
 	// handle whether or not partZ includes ED values or not
 	Z_type partZ_withED = 0, partZ_noED = 0;
 	if (isHybridZ) {

src/IntaRNA/PredictorMfeEns2dHeuristic.cpp

@@ -101,6 +101,7 @@ fillHybridZ()
 
 			// init as invalid boundary
 			*curCell = BestInteractionZ(0.0, RnaSequence::lastPos, RnaSequence::lastPos);
+			curCellEtotal = E_INF;
 
 			// check if positions can form interaction
 			if (	energy.isAccessible1(i1)
@@ -134,24 +135,20 @@ fillHybridZ()
 					// update overall partition function information for initial bps only
 					updateZ( i1,curCell->j1, i2,curCell->j2, curCell->val, true );
 
-					if(outConstraint.noLP) {
-						/////////////////////////////////////////
-						// check direct extension to the right of the noLP stacking
-						/////////////////////////////////////////
+				}
 
-						// direct cell access (const)
-						rightExt = &(hybridZ(i1+noLpShift,i2+noLpShift));
-						// check if right side can pair
-						if (Z_equal(rightExt->val, 0.0)) {
-							continue;
-						}
-						// check if interaction length is within boundary
-						if ( (rightExt->j1 +1 -i1) > energy.getAccessibility1().getMaxLength()
-							|| (rightExt->j2 +1 -i2) > energy.getAccessibility2().getMaxLength() )
-						{
-							continue;
-						}
+				if(outConstraint.noLP) {
+					/////////////////////////////////////////
+					// check direct extension to the right of the noLP stacking
+					/////////////////////////////////////////
 
+					// direct cell access (const)
+					rightExt = &(hybridZ(i1+noLpShift,i2+noLpShift));
+					// check if right side can pair and interaction length is within boundary
+					if (!Z_equal(rightExt->val, 0.0)
+						&& (rightExt->j1 +1 -i1) <= energy.getAccessibility1().getMaxLength()
+						&& (rightExt->j2 +1 -i2) <= energy.getAccessibility2().getMaxLength() )
+					{
 						// compute Z for direct extension with stacking
 						curZ = iStackZ * rightExt->val;
 
@@ -179,6 +176,12 @@ fillHybridZ()
 				// iterate over all loop sizes w1 (seq1) and w2 (seq2) (minus 1)
 				for (w1=1; w1-1 <= energy.getMaxInternalLoopSize1() && i1+w1+noLpShift<hybridZ.size1(); w1++) {
 				for (w2=1; w2-1 <= energy.getMaxInternalLoopSize2() && i2+w2+noLpShift<hybridZ.size2(); w2++) {
+					// For noLP, the adjacent continuation is already represented
+					// by the direct extension above. Counting it again as the
+					// (w1,w2)=(1,1) loop duplicates the same interaction paths.
+					if (noLpShift != 0 && w1 == 1 && w2 == 1) {
+						continue;
+					}
 					// direct cell access (const)
 					rightExt = &(hybridZ(i1+noLpShift+w1,i2+noLpShift+w2));
 					// check if right side can pair

src/IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.cpp

@@ -61,6 +61,7 @@ predict( const IndexRange & r1, const IndexRange & r2 )
 	if (seedHandler.fillSeed( 0, range_size1-1, 0, range_size2-1 ) == 0) {
 		// trigger empty interaction reporting
 		initOptima();
+		initZ();
 		reportOptima();
 		// stop computation
 		return;

src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp

@@ -62,6 +62,7 @@ predict( const IndexRange & r1, const IndexRange & r2 )
 	if (seedHandler.fillSeed( 0, range_size1-1, 0, range_size2-1 ) == 0) {
 		// trigger empty interaction reporting
 		initOptima();
+		initZ();
 		reportOptima();
 		// stop computation
 		return;
@@ -229,7 +230,7 @@ fillHybridZ_left( const size_t si1, const size_t si2 )
 					// get stacking energy to avoid recomputation in recursion below
 					iStackZ = energy.getBoltzmannWeight(energy.getE_interLeft(i1,i1+noLpShift,i2,i2+noLpShift));
 					// check just stacked
-					curZ += iStackZ + hybridZ_left(l1-noLpShift,l2-noLpShift);
+					curZ += iStackZ * hybridZ_left(l1-noLpShift,l2-noLpShift);
 				}
 
 				// check all combinations of decompositions into (i1,i2)..(k1,k2)-(j1,j2)