binaries/src/ViennaRNA/Utils/popt.c

   1 /*
   2    extract Zuker's p-optimal folds from subopt output.
   3    prints p-optimal structures followed by the list of pairs
   4    with respect to which they are optimal
   5
   6    input must be sorted by energy!
   7
   8    Ivo Hofacker, Walter Fontana
   9 */
  10
  11 #include <stdio.h>
  12 #include <string.h>
  13 #include <stdlib.h>
  14 #include "utils.h"
  15 #define PRIVATE static
  16 #define PUBLIC
  17
  18
  19 PRIVATE void
  20 help ()
  21 {
  22   printf ("p-optimal filter to subopt output. usage:\n"
  23           "RNAsubopt -s < seq | popt \n");
  24   exit(EXIT_FAILURE);
  25 }
  26
  27 /*-------------------------------------------------------------------------*/
  28
  29 int main (int argc, char *argv[])
  30 {
  31   char *line, *sequence, *structure;
  32   int i, j, length;
  33   float energy;
  34
  35   short *ptable;
  36   char *pair_seen;
  37   int  *iindx;
  38
  39   if (argc>1) help();
  40
  41   line = get_line (stdin);
  42   if (*line=='>') {
  43     free (line);
  44     line = get_line (stdin);
  45   }
  46
  47   /* read the sequence */
  48   sequence = (char *) space (sizeof (char) * (strlen(line)+1));
  49   (void) sscanf(line, "%s", sequence);
  50   free (line);
  51   length = (int) strlen(sequence);
  52
  53   iindx = (int *) space(sizeof(int)*(length+1));
  54   for (i=1; i<=length; i++)
  55     iindx[i] = ((length+1-i)*(length-i))/2 +length+1;
  56   pair_seen = (char *) space(((length+1)*(length+2))/2 * sizeof(char));
  57
  58   structure = (char *) space (sizeof (char) * (length+1));
  59
  60   /* get of suboptimal structures */
  61
  62   while ((line = get_line (stdin))) {
  63     int r, popt;
  64
  65     r = sscanf(line, "%s %f", structure, &energy);
  66     free(line);
  67     if (r==0) continue;
  68     ptable = make_pair_table(structure);
  69
  70     popt = 0;
  71     for (i = 0; i < length; i++) {
  72       if ((j=ptable[i])>i)
  73         if (pair_seen[iindx[i]-j]==0) {
  74           pair_seen[iindx[i]-j]=(char) 1;
  75           if (popt==0) {
  76             printf("%s %6.2f", structure, energy);
  77             popt = 1;
  78           }
  79           printf(" (%d,%d)", i,j);
  80         }
  81     }
  82     if (popt==1) printf("\n");
  83     free(ptable);
  84   }
  85   free(sequence); free(structure); free(iindx); free(pair_seen);
  86
  87   return 0;
  88 }