binaries/src/ViennaRNA/lib/plot_layouts.c

   1 /**
   2  * This file is a container for all plotting layout algorithms
   3  *
   4  *  c Ronny Lorenz
   5  *    The ViennaRNA Package
   6  */
   7
   8 #include <config.h>
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <math.h>
  12 #include "utils.h"
  13 #include "plot_layouts.h"
  14
  15 #ifdef _OPENMP
  16 #include <omp.h>
  17 #endif
  18
  19 #define PUBLIC
  20 #define PRIVATE static
  21
  22 PUBLIC  int     rna_plot_type = 1;  /* 0 = simple, 1 = naview, 2 = circular plot */
  23
  24 PRIVATE float   *angle;
  25 PRIVATE int     *loop_size, *stack_size;
  26 PRIVATE int     lp, stk;
  27
  28 PRIVATE void  loop(int i, int j, short *pair_table);
  29
  30 #ifdef _OPENMP
  31 /* NOTE: all threadprivate variables are uninitialized when entering a thread! */
  32 #pragma omp threadprivate(angle, loop_size, stack_size, lp, stk)
  33 #endif
  34
  35 /*---------------------------------------------------------------------------*/
  36
  37 PUBLIC int simple_xy_coordinates(short *pair_table, float *x, float *y)
  38 {
  39   float INIT_ANGLE=0.;     /* initial bending angle */
  40   float INIT_X = 100.;     /* coordinate of first digit */
  41   float INIT_Y = 100.;     /* see above */
  42   float RADIUS =  15.;
  43
  44   int i, length;
  45   float  alpha;
  46
  47   length = pair_table[0];
  48   angle =      (float*) space( (length+5)*sizeof(float) );
  49   loop_size  =   (int*) space( 16+(length/5)*sizeof(int) );
  50   stack_size =   (int*) space( 16+(length/5)*sizeof(int) );
  51   lp = stk = 0;
  52   loop(0, length+1, pair_table);
  53   loop_size[lp] -= 2;     /* correct for cheating with function loop */
  54
  55   alpha = INIT_ANGLE;
  56   x[0]  = INIT_X;
  57   y[0]  = INIT_Y;
  58
  59   for (i = 1; i <= length; i++) {
  60     x[i] = x[i-1]+RADIUS*cos(alpha);
  61     y[i] = y[i-1]+RADIUS*sin(alpha);
  62     alpha += PI-angle[i+1];
  63   }
  64   free(angle);
  65   free(loop_size);
  66   free(stack_size);
  67
  68   return length;
  69
  70 }
  71
  72 /*---------------------------------------------------------------------------*/
  73
  74 PRIVATE void loop(int i, int j, short *pair_table)
  75              /* i, j are the positions AFTER the last pair of a stack; i.e
  76                 i-1 and j+1 are paired. */
  77 {
  78   int    count = 2;   /* counts the VERTICES of a loop polygon; that's
  79                            NOT necessarily the number of unpaired bases!
  80                            Upon entry the loop has already 2 vertices, namely
  81                            the pair i-1/j+1.  */
  82
  83   int    r = 0, bubble = 0; /* bubble counts the unpaired digits in loops */
  84
  85   int    i_old, partner, k, l, start_k, start_l, fill, ladder;
  86   int    begin, v, diff;
  87   float  polygon;
  88
  89   short *remember;
  90
  91   remember = (short *) space((1+(j-i)/5)*2*sizeof(short));
  92
  93   i_old = i-1, j++;         /* j has now been set to the partner of the
  94                                previous pair for correct while-loop
  95                                termination.  */
  96   while (i != j) {
  97     partner = pair_table[i];
  98     if ((!partner) || (i==0))
  99       i++, count++, bubble++;
 100     else {
 101       count += 2;
 102       k = i, l = partner;    /* beginning of stack */
 103       remember[++r] = k;
 104       remember[++r] = l;
 105       i = partner+1;         /* next i for the current loop */
 106
 107       start_k = k, start_l = l;
 108       ladder = 0;
 109       do {
 110         k++, l--, ladder++;        /* go along the stack region */
 111       }
 112       while (pair_table[k] == l);
 113
 114       fill = ladder-2;
 115       if (ladder >= 2) {
 116         angle[start_k+1+fill] += PIHALF;   /*  Loop entries and    */
 117         angle[start_l-1-fill] += PIHALF;   /*  exits get an        */
 118         angle[start_k]        += PIHALF;   /*  additional PI/2.    */
 119         angle[start_l]        += PIHALF;   /*  Why ? (exercise)    */
 120         if (ladder > 2) {
 121           for (; fill >= 1; fill--) {
 122             angle[start_k+fill] = PI;    /*  fill in the angles  */
 123             angle[start_l-fill] = PI;    /*  for the backbone    */
 124           }
 125         }
 126       }
 127       stack_size[++stk] = ladder;
 128       loop(k, l, pair_table);
 129     }
 130   }
 131   polygon = PI*(count-2)/(float)count; /* bending angle in loop polygon */
 132   remember[++r] = j;
 133   begin = i_old < 0 ? 0 : i_old;
 134   for (v = 1; v <= r; v++) {
 135     diff  = remember[v]-begin;
 136     for (fill = 0; fill <= diff; fill++)
 137       angle[begin+fill] += polygon;
 138     if (v > r)
 139       break;
 140     begin = remember[++v];
 141   }
 142   loop_size[++lp] = bubble;
 143   free(remember);
 144 }
 145
 146 /*---------------------------------------------------------------------------*/
 147
 148 PUBLIC int simple_circplot_coordinates(short *pair_table, float *x, float *y){
 149   unsigned int  length = (unsigned int) pair_table[0];
 150   unsigned int  i;
 151   float         d = 2*PI/length;
 152   for(i=0; i < length; i++){
 153     x[i] = cos(i * d - PI/2);
 154     y[i] = sin(i * d - PI/2);
 155   }
 156   return length;
 157 }