Repository: ViennaRNA/ViennaRNA (language: python)

Issue:
Allow use of arbitrary vrna_ptable in ensemble defect computation 
Although the partition function algorithm is specific to ensembles of nested structures, the ensemble defect only relates indirectly to the ensemble type (via the partition function).
I encountered a situation were it was actually really useful to me to allow structures with pseudoknots in the calculation of the ensemble defect, so that's what this PR is about.

I tried to keep the existing API intact by introducing `vrna_ensemble_defect_pt()`, taking a `vrna_ptable` instead of a dot-bracket string. `vrna_ensemble_defect()` is just wrapped around this now.
The generated SWIG wrappers `ensemble_defect()` still take a structure string while pointing directly to the new function and take an additional parameter for the bracket types (with default value `VRNA_BRACKETS_RND` to keep the original behaviour).

I updated the corresponding documentation and added a few very simple and isolated tests to ensure I did not break the API. I hope it's clear enough to read.


Discussion / hints:
No idea why CI is failing though. The tests run fine on my machine...
The CI doesn't have libcheck.
Thank you for your contribution!
The extension of the Ensemble Defect computation seems reasonable to me, even though at this time, the ensemble will always consist of pseudo-knot free structures. As long as one keeps this in mind when interpreting the results, everything should be fine.
I'll merge your extension now and make it available in the next release of the ViennaRNA Package.


Fix patch:
diff --git a/interfaces/structure_utils.dox b/interfaces/structure_utils.dox
index 02bd53619..327254c12 100644
--- a/interfaces/structure_utils.dox
+++ b/interfaces/structure_utils.dox
@@ -60,4 +60,21 @@ This function is available as an overloaded function @p abstract_shapes() where
 the optional second parameter @p level defaults to 5.
 @endparblock
 
+@fn int vrna_bp_distance_pt(const short *pt1, const short *pt2)
+@scripting
+@parblock
+This function is available as an overloaded method @b bp_distance().
+@endparblock
+
+@fn int vrna_bp_distance(const char *str1, const char *str2)
+@scripting
+@parblock
+This function is available as an overloaded method @b bp_distance().
+Note that the SWIG wrapper takes two structure in dot-bracket notation and converts them into
+pair tables using vrna_ptable_from_string(). The resulting pair tables are then internally
+passed to vrna_bp_distance_pt(). To control which kind of matching brackets will be used
+during conversion, the optional argument @p options can be used. See also the description of
+vrna_ptable_from_string() for available options. (default: @b VRNA_BRACKETS_RND).
+@endparblock
+
 */
diff --git a/interfaces/structure_utils.i b/interfaces/structure_utils.i
index cd6fb0907..1cf0eb1c1 100644
--- a/interfaces/structure_utils.i
+++ b/interfaces/structure_utils.i
@@ -525,10 +525,35 @@ std::vector<int> my_loopidx_from_ptable(std::vector<int> pt);
 
 %{
   int
-  my_bp_distance(const char *str1,
-                 const char *str2)
+  my_bp_distance(std::string str1,
+                 std::string str2,
+                 unsigned int options = VRNA_BRACKETS_RND)
   {
-    return vrna_bp_distance(str1,str2);
+    int dist = 0;
+    short int *pt1, *pt2;
+
+    pt1 = vrna_ptable_from_string(str1.c_str(), options);
+    pt2 = vrna_ptable_from_string(str2.c_str(), options);
+
+    dist = vrna_bp_distance_pt(pt1, pt2);
+
+    free(pt1);
+    free(pt2);
+
+    return dist;
+  }
+
+  int
+  my_bp_distance(std::vector<int> pt1,
+                 std::vector<int> pt2)
+  {
+    std::vector<short> pt1_v_short;
+    std::vector<short> pt2_v_short;
+
+    transform(pt1.begin(), pt1.end(), back_inserter(pt1_v_short), convert_vecint2vecshort);
+    transform(pt2.begin(), pt2.end(), back_inserter(pt2_v_short), convert_vecint2vecshort);
+
+    return vrna_bp_distance_pt((short*)&pt1_v_short[0], (short*)&pt2_v_short[0]);
   }
 
   double
@@ -545,7 +570,8 @@ std::vector<int> my_loopidx_from_ptable(std::vector<int> pt);
 %feature("kwargs") my_bp_distance;
 #endif
 
-int     my_bp_distance(const char *str1, const char *str2);
+int     my_bp_distance(std::string str1, std::string str2, unsigned int options = VRNA_BRACKETS_RND);
+int     my_bp_distance(std::vector<int> pt1, std::vector<int> pt2);
 double  my_dist_mountain(std::string str1, std::string str2, unsigned int p = 1);
 
 /************************************/
diff --git a/src/ViennaRNA/utils/structure_utils.c b/src/ViennaRNA/utils/structure_utils.c
index c15e038f4..52db42cc0 100644
--- a/src/ViennaRNA/utils/structure_utils.c
+++ b/src/ViennaRNA/utils/structure_utils.c
@@ -539,34 +539,47 @@ vrna_db_from_WUSS(const char *wuss)
 /*---------------------------------------------------------------------------*/
 
 PUBLIC int
-vrna_bp_distance(const char *str1,
-                 const char *str2)
+vrna_bp_distance_pt(const short *pt1,
+                    const short *pt2)
 {
   /* dist = {number of base pairs in one structure but not in the other} */
   /* same as edit distance with pair_open pair_close as move set */
   int   dist;
   short i, l;
-  short *t1, *t2;
 
   dist  = 0;
-  t1    = vrna_ptable(str1);
-  t2    = vrna_ptable(str2);
 
-  if (t1 && t2) {
-    l = (t1[0] < t2[0]) ? t1[0] : t2[0]; /* minimum of the two lengths */
+  if (pt1 && pt2) {
+    l = (pt1[0] < pt2[0]) ? pt1[0] : pt2[0]; /* minimum of the two lengths */
 
     for (i = 1; i <= l; i++)
-      if (t1[i] != t2[i]) {
-        if (t1[i] > i)
+      if (pt1[i] != pt2[i]) {
+        if (pt1[i] > i)
           dist++;
 
-        if (t2[i] > i)
+        if (pt2[i] > i)
           dist++;
       }
   }
 
-  free(t1);
-  free(t2);
+  return dist;
+}
+
+
+PUBLIC int
+vrna_bp_distance(const char *str1,
+                 const char *str2)
+{
+  int   dist = 0;
+  short *pt1, *pt2;
+
+  pt1    = vrna_ptable(str1);
+  pt2    = vrna_ptable(str2);
+
+  dist = vrna_bp_distance_pt(pt1, pt2);
+
+  free(pt1);
+  free(pt2);
 
   return dist;
 }
diff --git a/src/ViennaRNA/utils/structures.h b/src/ViennaRNA/utils/structures.h
index 69397cac5..a06a76f1f 100644
--- a/src/ViennaRNA/utils/structures.h
+++ b/src/ViennaRNA/utils/structures.h
@@ -697,13 +697,33 @@ int *
 vrna_loopidx_from_ptable(const short *pt);
 
 
+/**
+ *  @brief Compute the "base pair" distance between two pair tables pt1 and pt2 of secondary structures.
+ *
+ *  The pair tables should have the same length.
+ *  dist = number of base pairs in one structure but not in the other
+ *  same as edit distance with open-pair close-pair as move-set
+ *
+ *  @see vrna_bp_distance()
+ *
+ *  @param pt1   First structure in dot-bracket notation
+ *  @param pt2   Second structure in dot-bracket notation
+ *  @return       The base pair distance between pt1 and pt2
+ */
+int
+vrna_bp_distance_pt(const short *pt1,
+                    const short *pt2);
+
 /**
  *  @brief Compute the "base pair" distance between two secondary structures s1 and s2.
  *
+ *  This is a wrapper around @b vrna_bp_distance_pt().
  *  The sequences should have the same length.
  *  dist = number of base pairs in one structure but not in the other
  *  same as edit distance with open-pair close-pair as move-set
  *
+ *  @see vrna_bp_distance_pt()
+ *
  *  @param str1   First structure in dot-bracket notation
  *  @param str2   Second structure in dot-bracket notation
  *  @return       The base pair distance between str1 and str2


Test patch:
diff --git a/tests/python/test-RNA-utils.py b/tests/python/test-RNA-utils.py
index 3a1a95001..efaf38a86 100644
--- a/tests/python/test-RNA-utils.py
+++ b/tests/python/test-RNA-utils.py
@@ -33,6 +33,25 @@ def test_basePairDistance(self):
         self.assertEqual(d,3)
 
 
+    def test_basePairDistancePairTable(self):
+        print "test_basePairDistancePairTable"
+        pt1 = RNA.ptable("(((.(((...))))))")
+        pt2 = RNA.ptable("(((..........)))")
+
+        d = RNA.bp_distance(pt1,pt2)
+
+        self.assertEqual(d,3)
+
+
+    def test_basePairDistancePseudoknot(self):
+        print "test_basePairDistancePseudoknot"
+        d1 = RNA.bp_distance(struct1, struct1pk)
+        d2 = RNA.bp_distance(struct1, struct1pk, RNA.BRACKETS_ANY)
+
+        self.assertEqual(d1, 2)
+        self.assertEqual(d2, 3)
+
+
     def test_plists(self):
         print "test_plists"
         plist = RNA.plist(struct1,0.6)
diff --git a/tests/python3/test-RNA-utils.py3 b/tests/python3/test-RNA-utils.py3
index ed3f33851..96114dabb 100644
--- a/tests/python3/test-RNA-utils.py3
+++ b/tests/python3/test-RNA-utils.py3
@@ -33,6 +33,25 @@ class GeneralTests(unittest.TestCase):
         self.assertEqual(d,3)
 
 
+    def test_basePairDistancePairTable(self):
+        print("test_basePairDistancePairTable")
+        pt1 = RNA.ptable("(((.(((...))))))")
+        pt2 = RNA.ptable("(((..........)))")
+
+        d = RNA.bp_distance(pt1,pt2)
+
+        self.assertEqual(d,3)
+
+
+    def test_basePairDistancePseudoknot(self):
+        print("test_basePairDistancePseudoknot")
+        d1 = RNA.bp_distance(struct1, struct1pk)
+        d2 = RNA.bp_distance(struct1, struct1pk, RNA.BRACKETS_ANY)
+
+        self.assertEqual(d1, 2)
+        self.assertEqual(d2, 3)
+
+
     def test_plists(self):
         print("test_plists")
         plist = RNA.plist(struct1,0.6)