WSTester updated to work plus hopefully all the other changes that need to go into...

[jabaws.git] / binaries / src / ViennaRNA / lib / mm.c

diff --git a/binaries/src/ViennaRNA/lib/mm.c b/binaries/src/ViennaRNA/lib/mm.c
new file mode 100644 (file)
index 0000000..248dff4
--- /dev/null
+++ b/binaries/src/ViennaRNA/lib/mm.c
@@ -0,0 +1,99 @@
+/*
+      Implementation of Nussinov Maximum Matching
+      Ronny Lorenz
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <ctype.h>
+#include <string.h>
+#include "utils.h"
+#include "energy_par.h"
+#include "fold_vars.h"
+#include "pair_mat.h"
+#include "params.h"
+
+/* the encoded string MUST have the length of the sequence at position 0!!! */
+PUBLIC unsigned int maximumMatching(const char *string){
+  unsigned int i, j, l, length, max = 0;
+  unsigned int *mm;            /* holds maximum matching on subsequence [i,j] */
+  short *encodedString = encode_sequence(string, 0);
+  int *iindx = get_iindx((unsigned) encodedString[0]);
+  make_pair_matrix();
+  length = (unsigned int)encodedString[0];
+  mm = (unsigned int *) space(sizeof(unsigned int)*((length*(length+1))/2+2));
+  for(j = 1; j<=length; j++)
+    for(i=(j>TURN?(j-TURN):1); i<j; i++)
+      mm[iindx[i]-j] = 0;
+  for(i=length-TURN-1;i>0; i--)
+    for(j=i+TURN+1; j<= length; j++){
+      max = mm[iindx[i]-j+1];
+      for(l=j-TURN-1; l>=i; l--)
+        if(pair[encodedString[l]][encodedString[j]])
+          max = MAX2(max, ((l>i) ? mm[iindx[i]-l+1] : 0) + 1 + mm[iindx[l+1]-j+1]);
+       mm[iindx[i]-j] = max;
+    }
+  max = mm[iindx[1]-length];
+  free(mm);
+  free(iindx);
+  free(encodedString);
+  return max;
+}
+
+/* the encoded string MUST have the length of the sequence at position 0!!! */
+PUBLIC unsigned int *maximumMatchingConstraint(const char *string, short *ptable){
+  unsigned int i, j, l, length, max = 0;
+  unsigned int *mm;            /* holds maximum matching on subsequence [i,j] */
+  short *encodedString = encode_sequence(string, 0);
+  int *iindx = get_iindx((unsigned) encodedString[0]);
+  make_pair_matrix();
+  length = (unsigned int)encodedString[0];
+  mm = (unsigned int *) space(sizeof(unsigned int)*((length*(length+1))/2+2));
+  for(j = 1; j<=length; j++)
+    for(i=(j>TURN?(j-TURN):1); i<j; i++)
+      mm[iindx[i]-j] = 0;
+  for(i=length-TURN-1;i>0; i--)
+    for(j=i+TURN+1; j<= length; j++){
+      max = mm[iindx[i]-j+1];
+      for(l=j-TURN-1; l>=i; l--){
+        if(pair[encodedString[l]][encodedString[j]]){
+          if(ptable[l] != j)
+            max = MAX2(max, ((l>i) ? mm[iindx[i]-l+1] : 0) + 1 + mm[iindx[l+1]-j+1]);
+        }
+      }
+      mm[iindx[i]-j] = max;
+    }
+  free(iindx);
+  free(encodedString);
+  return mm;
+}
+
+/* the encoded string MUST have the length of the sequence at position 0!!! */
+PUBLIC unsigned int *maximumMatching2Constraint(const char *string, short *ptable, short *ptable2){
+  unsigned int i, j, l, length, max = 0;
+  unsigned int *mm;            /* holds maximum matching on subsequence [i,j] */
+  short *encodedString = encode_sequence(string, 0);
+  int *iindx = get_iindx((unsigned) encodedString[0]);
+  make_pair_matrix();
+  length = (unsigned int)encodedString[0];
+  mm = (unsigned int *) space(sizeof(unsigned int)*((length*(length+1))/2+2));
+  for(j = 1; j<=length; j++)
+    for(i=(j>TURN?(j-TURN):1); i<j; i++)
+      mm[iindx[i]-j] = 0;
+  for(i=length-TURN-1;i>0; i--)
+    for(j=i+TURN+1; j<= length; j++){
+      max = mm[iindx[i]-j+1];
+      for(l=j-TURN-1; l>=i; l--){
+        if(pair[encodedString[l]][encodedString[j]]){
+          if(ptable[l] != j && ptable2[l] != j)
+            max = MAX2(max, ((l>i) ? mm[iindx[i]-l+1] : 0) + 1 + mm[iindx[l+1]-j+1]);
+        }
+      }
+      mm[iindx[i]-j] = max;
+    }
+  free(iindx);
+  free(encodedString);
+  return mm;
+}
+