binaries/src/ViennaRNA/lib/mm.c

   1 /*
   2       Implementation of Nussinov Maximum Matching
   3       Ronny Lorenz
   4 */
   5
   6 #include <stdio.h>
   7 #include <stdlib.h>
   8 #include <math.h>
   9 #include <ctype.h>
  10 #include <string.h>
  11 #include "utils.h"
  12 #include "energy_par.h"
  13 #include "fold_vars.h"
  14 #include "pair_mat.h"
  15 #include "params.h"
  16
  17 /* the encoded string MUST have the length of the sequence at position 0!!! */
  18 PUBLIC unsigned int maximumMatching(const char *string){
  19   unsigned int i, j, l, length, max = 0;
  20   unsigned int *mm;            /* holds maximum matching on subsequence [i,j] */
  21   short *encodedString = encode_sequence(string, 0);
  22   int *iindx = get_iindx((unsigned) encodedString[0]);
  23   make_pair_matrix();
  24   length = (unsigned int)encodedString[0];
  25   mm = (unsigned int *) space(sizeof(unsigned int)*((length*(length+1))/2+2));
  26   for(j = 1; j<=length; j++)
  27     for(i=(j>TURN?(j-TURN):1); i<j; i++)
  28       mm[iindx[i]-j] = 0;
  29   for(i=length-TURN-1;i>0; i--)
  30     for(j=i+TURN+1; j<= length; j++){
  31       max = mm[iindx[i]-j+1];
  32       for(l=j-TURN-1; l>=i; l--)
  33         if(pair[encodedString[l]][encodedString[j]])
  34           max = MAX2(max, ((l>i) ? mm[iindx[i]-l+1] : 0) + 1 + mm[iindx[l+1]-j+1]);
  35        mm[iindx[i]-j] = max;
  36     }
  37   max = mm[iindx[1]-length];
  38   free(mm);
  39   free(iindx);
  40   free(encodedString);
  41   return max;
  42 }
  43
  44 /* the encoded string MUST have the length of the sequence at position 0!!! */
  45 PUBLIC unsigned int *maximumMatchingConstraint(const char *string, short *ptable){
  46   unsigned int i, j, l, length, max = 0;
  47   unsigned int *mm;            /* holds maximum matching on subsequence [i,j] */
  48   short *encodedString = encode_sequence(string, 0);
  49   int *iindx = get_iindx((unsigned) encodedString[0]);
  50   make_pair_matrix();
  51   length = (unsigned int)encodedString[0];
  52   mm = (unsigned int *) space(sizeof(unsigned int)*((length*(length+1))/2+2));
  53   for(j = 1; j<=length; j++)
  54     for(i=(j>TURN?(j-TURN):1); i<j; i++)
  55       mm[iindx[i]-j] = 0;
  56   for(i=length-TURN-1;i>0; i--)
  57     for(j=i+TURN+1; j<= length; j++){
  58       max = mm[iindx[i]-j+1];
  59       for(l=j-TURN-1; l>=i; l--){
  60         if(pair[encodedString[l]][encodedString[j]]){
  61           if(ptable[l] != j)
  62             max = MAX2(max, ((l>i) ? mm[iindx[i]-l+1] : 0) + 1 + mm[iindx[l+1]-j+1]);
  63         }
  64       }
  65       mm[iindx[i]-j] = max;
  66     }
  67   free(iindx);
  68   free(encodedString);
  69   return mm;
  70 }
  71
  72 /* the encoded string MUST have the length of the sequence at position 0!!! */
  73 PUBLIC unsigned int *maximumMatching2Constraint(const char *string, short *ptable, short *ptable2){
  74   unsigned int i, j, l, length, max = 0;
  75   unsigned int *mm;            /* holds maximum matching on subsequence [i,j] */
  76   short *encodedString = encode_sequence(string, 0);
  77   int *iindx = get_iindx((unsigned) encodedString[0]);
  78   make_pair_matrix();
  79   length = (unsigned int)encodedString[0];
  80   mm = (unsigned int *) space(sizeof(unsigned int)*((length*(length+1))/2+2));
  81   for(j = 1; j<=length; j++)
  82     for(i=(j>TURN?(j-TURN):1); i<j; i++)
  83       mm[iindx[i]-j] = 0;
  84   for(i=length-TURN-1;i>0; i--)
  85     for(j=i+TURN+1; j<= length; j++){
  86       max = mm[iindx[i]-j+1];
  87       for(l=j-TURN-1; l>=i; l--){
  88         if(pair[encodedString[l]][encodedString[j]]){
  89           if(ptable[l] != j && ptable2[l] != j)
  90             max = MAX2(max, ((l>i) ? mm[iindx[i]-l+1] : 0) + 1 + mm[iindx[l+1]-j+1]);
  91         }
  92       }
  93       mm[iindx[i]-j] = max;
  94     }
  95   free(iindx);
  96   free(encodedString);
  97   return mm;
  98 }
  99