1#ifndef RNAPUZZLER_CONFIG_H
2#define RNAPUZZLER_CONFIG_H
3
4/*
5 *      RNApuzzler template config
6 *
7 *      c  Daniel Wiegreffe, Daniel Alexander, Dirk Zeckzer
8 *      ViennaRNA package
9 */
10#include <stdlib.h>
11#include <math.h>
12
13#include "ViennaRNA/utils/basic.h"
14
15#include "definitions.inc"
16#include "vector_math.inc"
17#include "../headers/configtree_struct.h"
18#include "../headers/tBaseInformation_struct.h"
19#include "ViennaRNA/plotting/RNApuzzler/RNApuzzler.h"
20
21PRIVATE double
22getArcAngle(const config  *cfg,
23            const int     currentArc);
24
25
26PRIVATE double
27getArcAngleDegree(const config  *cfg,
28                  const int     currentArc);
29
30
31/**
32 * @brief cfgPrintConfig
33 *      - prints the given config to the terminal
34 * @param config
35 *      - config to print
36 */
37PRIVATE void
38cfgPrintConfig(config *config);
39
40
41/**
42 * @brief cfgGenerateConfig
43 *      - generates configurations for each loop
44 * @param pair_table
45 *      - the RNA's pairing information
46 * @param baseInformation
47 *      - array of tBaseInformation annotations (for saving configs)
48 * @param unpaired
49 *      - default length of backbone lines
50 * @param paired
51 *      - default distance between paired bases
52 */
53PRIVATE void
54cfgGenerateConfig(const short *const  pair_table,
55                  tBaseInformation    *baseInformation,
56                  const double        unpaired,
57                  const double        paired);
58
59
60/**
61 * @brief cfgCloneConfig
62 *      - prints the given config to the terminal
63 * @param cfg
64 *      - config to clone
65 * @return clone of cfg
66 */
67PRIVATE config *
68cfgCloneConfig(const config *cfg);
69
70
71/**
72 * @brief cfgFreeConfig
73 *      - prints the given config to the terminal
74 * @param cfg
75 *      - config to free
76 */
77PRIVATE void
78cfgFreeConfig(config *cfg);
79
80
81/**
82 * Function to apply a set of config changes.
83 *
84 * By passing -1.0 to radiusNew you will apply the minimum possible radius for this loop while not allowing to shrink the loop.
85 * In case it would actually shrink a default relative increase is applied to the radius.
86 *
87 * By passing 0.0 to radiusNew you set the minimum possible radius no matter if there would be some shrinking or not.
88 * @param loopName the node to change
89 * @param deltaCfg array of angles to change (in degree)
90 * @param radiusNew desired radius to set while applying deltas
91 * @param puzzler
92 */
93PRIVATE double
94cfgApplyChanges(config                            *cfg,
95                const char                        loopName,
96                const double                      *deltaCfg,
97                const double                      radiusNew,
98                const vrna_plot_options_puzzler_t *puzzler);
99
100
101/**
102 * @brief cfgIsValid
103 *      - check if config is valid
104 * @param config
105 *      - the config to check
106 * @param deltaCfg
107 *      - array of config changes.
108 *        contains diff-values for each config angle.
109 *        in degree format
110 * @return true iff config is valid
111 */
112PRIVATE short
113cfgIsValid(config       *config,
114           const double *deltaCfg);
115
116
117/**
118 * @brief intToMotiv
119 *      - converts a given number into a readable format
120 * @param _int
121 *      - motiv name as number (int)
122 * @return
123 *      - motiv name as char
124 */
125PRIVATE char
126intToMotiv(const int _int);
127
128
129/*
130 *--------------------------------------------------------------------------
131 *---   create, copy, and free config   ------------------------------------
132 *--------------------------------------------------------------------------
133 */
134
135/**
136 * @brief cfgCreateConfig
137 *      - constructor-like method for creating a config
138 * @param radius
139 *      - radius used for drawing that loop
140 * @return
141 *      - an initialized config struct
142 */
143PRIVATE config *
144cfgCreateConfig(const double radius)
145{
146  config *cfg = (config *)vrna_alloc(1 * sizeof(config));
147
148  cfg->radius         = radius;
149  cfg->minRadius      = radius;
150  cfg->defaultRadius  = radius;
151
152  cfg->cfgArcs      = NULL;
153  cfg->numberOfArcs = 0;
154
155  return cfg;
156}
157
158
159/**
160 * @brief cfgCreateConfigArc
161 *      - constructor-like method for adding a new config entry to a given config
162 * @param angle
163 *      - angle (radiant) that can be found between stems 'from' and 'to' later on
164 * @param numberOfArcSegments
165 *      - number of arc segments between stems 'from' and 'to'
166 * @return
167 *      - an initialized configArc struct
168 */
169PRIVATE configArc
170cfgCreateConfigArc(const double angle,
171                   const int    numberOfArcSegments)
172{
173  configArc newConfigArc;
174
175  newConfigArc.numberOfArcSegments = numberOfArcSegments;
176
177  newConfigArc.arcAngle = angle;
178
179  return newConfigArc;
180}
181
182
183PRIVATE config *
184cfgCloneConfig(const config *cfg)
185{
186  config *clonedCfg = (config *)vrna_alloc(sizeof(config));
187
188  clonedCfg->radius         = cfg->radius;
189  clonedCfg->minRadius      = cfg->minRadius;
190  clonedCfg->defaultRadius  = cfg->defaultRadius;
191  clonedCfg->numberOfArcs   = cfg->numberOfArcs;
192
193  const int numberOfArcs = cfg->numberOfArcs;
194  clonedCfg->cfgArcs = (configArc *)vrna_alloc(numberOfArcs * sizeof(configArc));
195
196  for (int currentArc = 0; currentArc < numberOfArcs; ++currentArc) {
197    clonedCfg->cfgArcs[currentArc].numberOfArcSegments =
198      cfg->cfgArcs[currentArc].numberOfArcSegments;
199    clonedCfg->cfgArcs[currentArc].arcAngle = cfg->cfgArcs[currentArc].arcAngle;
200  }
201
202  return clonedCfg;
203}
204
205
206PRIVATE void
207cfgFreeConfig(config *cfg)
208{
209  free(cfg->cfgArcs);
210  free(cfg);
211}
212
213
214/* documentation at header file */
215PRIVATE void
216cfgPrintConfig(config *cfg)
217{
218  short useDegree = 1;
219
220  for (int currentArc = 0; currentArc < cfg->numberOfArcs; ++currentArc) {
221    double angle = getArcAngle(cfg, currentArc);
222
223    if (useDegree)
224      angle = toDegree(angle);
225  }
226}
227
228
229/*
230 *--------------------------------------------------------------------------
231 *---   access to config elements   ----------------------------------------
232 *--------------------------------------------------------------------------
233 */
234PRIVATE double
235getArcAngle(const config  *cfg,
236            const int     currentArc)
237{
238  return cfg->cfgArcs[currentArc].arcAngle;
239}
240
241
242PRIVATE double
243getArcAngleDegree(const config  *cfg,
244                  const int     currentArc)
245{
246  return toDegree(getArcAngle(cfg, currentArc));
247}
248
249
250/*
251 *--------------------------------------------------------------------------
252 *---   radius computation   -----------------------------------------------
253 *--------------------------------------------------------------------------
254 */
255
256/**
257 * Approximate the radius of a circle required to draw m base pairs
258 * with a distance of 'a' to each other, and n (unpaired) consecutive
259 * nucleotides with a distance of b over a specified angle.
260 *
261 * Uses a Newton iteration to approximate solution
262 *
263 * @param a paired
264 * @param b unpaired
265 * @param m #stems
266 * @param n #backbones
267 * @param angle angle
268 */
269PRIVATE double
270approximateConfigArcRadius(const double a,
271                           const double b,
272                           const short  m,
273                           const short  n,
274                           double       angle)
275{
276  const short MAX_ITERATIONS = 1000;
277
278  /*
279   * calculation:
280   *
281   * be s the length of a line at the circle (paired or unpaired / a or b)
282   * the angle over such a single line be alpha
283   * alpha = angle / ( m + n )
284   *
285   * for such a single line the following equation holds (where r is the radius of the circle)
286   * sin( alpha / 2 ) = ( s / 2 ) / r
287   * r = ( s / 2 ) / sin( alpha / 2 )
288   * r = ( s / 2 ) / sin( ( angle / ( m + n ) ) / 2 )
289   *
290   * now we replace s with a or b to get the upper or lower bound for the radius interval
291   */
292  double  lowerBound  = (b / 2) / sin((angle / (m + n)) / 2);
293  double  upperBound  = (a / 2) / sin((angle / (m + n)) / 2);
294  /* printf("new: lower %f upper %f\n", lowerBound, upperBound); */
295  double  rtn = 0.5 * (lowerBound + upperBound);
296
297  /*
298   * there is a minimum valid radius!
299   * if rtn is smaller than 0.5*a or 0.5*b than the result will become nan
300   */
301  rtn = fmax(rtn, 0.5 * a);
302  rtn = fmax(rtn, 0.5 * b);
303
304  /*
305   * printf("lower %f upper %f %f\n", lowerBound, upperBound, rtn);
306   * printf("Angle %f\n", angle);
307   */
308
309  /* short isERROR = 0; */
310
311  int j = 0;
312
313  for (j = 0; j < MAX_ITERATIONS; j++) {
314    double dx = 2 * (m * asin(a / (2 * rtn)) + n * asin(b / (2 * rtn)) - (angle / 2))
315                / (-(a * m / (rtn * sqrt(rtn * rtn - a * a / 4)) + b * n /
316                     (rtn * sqrt(rtn * rtn - b * b / 4))));
317    rtn -= dx;
318
319    if (fabs(dx) < EPSILON_3)
320
321      break;
322  }
323
324  if (rtn < lowerBound)
325    rtn = lowerBound;
326  else if (rtn > upperBound)
327    rtn = upperBound;
328
329  return rtn;
330}
331
332
333PRIVATE double
334approximateConfigRadius(const config  *cfg,
335                        const double  unpaired,
336                        const double  paired)
337{
338  /*
339   * calculate a fitting radius for each arc without compressing or stretching arc segments
340   * return the maximum of those values as the radius fitting for the loop
341   */
342  double r = 0;
343
344  for (int currentArc = 0; currentArc < cfg->numberOfArcs; ++currentArc) {
345    int     stems               = 1;
346    int     numberOfArcSegments = (cfg->cfgArcs[currentArc]).numberOfArcSegments;
347    double  angle               = getArcAngle(cfg, currentArc);
348
349    double  tempR = approximateConfigArcRadius(paired, unpaired, stems, numberOfArcSegments, angle);
350
351    if (tempR > r)
352
353      r = tempR;
354  }
355
356  return r;
357}
358
359
360/* -------------------------------------------------------------------------------------------------------------------------- */
361
362/**
363 * @brief cfgGenerateDefaultConfig
364 *      - generates a config that resembles a drawing without any
365 *        constraints given by config input for the given loop
366 * @param pair_table
367 *      - the RNA's pairing information
368 * @param start
369 *      - index of the loop's first base
370 * @param unpaired
371 *      - default length of backbone lines
372 * @param paired
373 *      - default distance between paired bases
374 * @param radius
375 *      - radius for the given loop
376 * @return
377 *      - complete config for that loop
378 */
379PRIVATE config *
380cfgGenerateDefaultConfig(const short *const pair_table,
381                         const int          start,
382                         const int          unpaired,
383                         const int          paired,
384                         const double       radius)
385{
386  /* create loop configuration */
387  config  *cfg = cfgCreateConfig(radius);
388
389  /* compute angles for paired and unpaired bases */
390  double  anglePaired   = 2 * asin(paired / (2 * radius));        /* angle over paired */
391  double  angleUnpaired = 2 * asin(unpaired / (2 * radius));      /* angle over unpaired */
392
393  /* initialize values for first arc */
394  int     arcUnpaired = 0;
395  double  angleArc;  /* alpha + numBackbones * beta */
396
397  /* start with first base after parent stem */
398  int     i = start + 1;
399
400  while (i <= pair_table[start]) {
401    /* until last base at parent stem */
402    if (pair_table[i] == 0) {
403      /* arc */
404      i++;
405    } else {
406      /* increment number of arcs */
407      ++(cfg->numberOfArcs);
408
409      if (i != pair_table[start])
410        /* skip subtree at stem */
411        i = pair_table[i] + 1;
412      else
413        /* parent stem -> finish */
414        break;
415    }
416  }
417
418  cfg->cfgArcs = (configArc *)vrna_alloc(cfg->numberOfArcs * sizeof(configArc));
419
420  /* start with first base after parent stem */
421  i = start + 1;
422  int currentArc          = 0;
423  int numberOfArcSegments = 0;
424
425  while (i <= pair_table[start]) {
426    /* until last base at parent stem */
427    if (pair_table[i] == 0) {
428      /* arc */
429      arcUnpaired++;
430      i++;
431    } else {
432      /* stem: create arc */
433      angleArc                  = anglePaired + (arcUnpaired + 1) * angleUnpaired;
434      numberOfArcSegments       = arcUnpaired + 1;
435      cfg->cfgArcs[currentArc]  = cfgCreateConfigArc(angleArc, numberOfArcSegments);
436      ++currentArc;
437
438      if (i != pair_table[start]) {
439        /* initialize values for next arc */
440        arcUnpaired = 0;
441
442        /* skip subtree at stem */
443        i = pair_table[i] + 1;
444      } else {
445        /* parent stem -> finish */
446        break;
447      }
448    }
449  }
450
451  return cfg;
452}
453
454
455PRIVATE void
456cfgGenHandleStem(int                baseNr,
457                 const short *const pair_table,
458                 tBaseInformation   *baseInformation,
459                 const double       unpaired,
460                 const double       paired);
461
462
463/**
464 * @brief cfgGenHandleLoop
465 *      - recursively iterates through the RNA and generates default configurations.
466 *        Alternates with corresponding handleStem method.
467 * @param baseNr
468 *      - index of the loop's first base
469 * @param pair_table
470 *      - the RNA's pairing information
471 * @param baseInformation
472 *      - array of tBaseInformation annotations (to save config)
473 */
474PRIVATE void
475cfgGenHandleLoop(int                baseNr,
476                 const short *const pair_table,
477                 tBaseInformation   *baseInformation,
478                 const double       unpaired,
479                 const double       paired)
480{
481  int start = baseNr;
482  int end   = pair_table[baseNr];
483
484  int unpairedCount = 0;
485  int stemCount     = 1;
486
487  /* count stems and unpaired bases to use for bulge detection */
488  int i = start + 1;
489
490  while (i < end) {
491    if (pair_table[i] == 0) {
492      /* unpaired base */
493      unpairedCount++;
494      i++;
495    } else if (pair_table[i] > i) {
496      /* found new stem */
497      stemCount++;
498      i = pair_table[i];
499    } else {
500      /* returned from stem */
501      i++;
502    }
503  }
504
505  short isBulge = (stemCount == 2 && unpairedCount == 1);
506  if (isBulge) {
507    if (pair_table[start + 1] == 0)
508      /* unpaired on left strand */
509      cfgGenHandleStem(start + 2, pair_table, baseInformation, unpaired, paired);
510    else
511      /* unpaired on the right strand */
512      cfgGenHandleStem(start + 1, pair_table, baseInformation, unpaired, paired);
513  } else {
514    int     m             = stemCount;                  /* compare RNApuzzler.c -> f_handle_loop */
515    int     n             = unpairedCount + stemCount;  /* compare RNApuzzler.c -> f_handle_loop */
516    double  defaultRadius = approximateConfigArcRadius(paired, unpaired, m, n, MATH_TWO_PI);
517    config  *cfgLoop      = cfgGenerateDefaultConfig(pair_table,
518                                                     start,
519                                                     unpaired,
520                                                     paired,
521                                                     defaultRadius);
522    baseInformation[start].config = cfgLoop;
523
524    int i = start + 1;
525    while (i < end) {
526      if (pair_table[i] == 0) {
527        /* unpaired base */
528        i++;
529      } else if (pair_table[i] > i) {
530        /* found new stem */
531        cfgGenHandleStem(i, pair_table, baseInformation, unpaired, paired);
532        i = pair_table[i];
533      } else {
534        /* returned from stem */
535        i++;
536      }
537    }
538  }
539}
540
541
542/**
543 * @brief cfgGenHandleStem
544 *      - recursively iterates through the RNA and generates default configurations.
545 *        Alternates with corresponding handleLoop method.
546 * @param baseNr
547 *      - index of the stem's first base
548 * @param pair_table
549 *      - the RNA's pairing information
550 * @param baseInformation
551 *      - array of tBaseInformation annotations (to save config)
552 */
553PRIVATE void
554cfgGenHandleStem(int                baseNr,
555                 const short *const pair_table,
556                 tBaseInformation   *baseInformation,
557                 const double       unpaired,
558                 const double       paired)
559{
560  /* iterating over the stem and calling cfgGenHandleLoop as soon as a loop is found */
561  short continueStem  = 1;
562  int   i             = baseNr;
563
564  while (continueStem) {
565    if (pair_table[i + 1] == pair_table[i] - 1) {
566      i++;
567    } else {
568      /* found unpaired above stem */
569      cfgGenHandleLoop(i, pair_table, baseInformation, unpaired, paired);
570      continueStem = 0;
571    }
572  }
573}
574
575
576/* documentation at header file */
577PRIVATE void
578cfgGenerateConfig(const short *const  pair_table,
579                  tBaseInformation    *baseInformation,
580                  const double        unpaired,
581                  const double        paired)
582{
583  /*
584   * iterate over RNA
585   * for each loop generate a default config
586   */
587  int length  = pair_table[0];
588  int i       = 1;
589
590  while (i < length) {
591    if (pair_table[i] == 0) {
592      /* unpaired at exterior loop */
593      i++;
594    } else if (pair_table[i] > i) {
595      /* found stem */
596      cfgGenHandleStem(i, pair_table, baseInformation, unpaired, paired);
597      i = pair_table[i];
598    } else {
599      /* returned from stem */
600      i++;
601    }
602  }
603}
604
605
606/*
607 *--------------------------------------------------------------------------
608 *---   set and update config elements   -----------------------------------
609 *--------------------------------------------------------------------------
610 */
611
612/* documentation at header file */
613/**
614 * @brief cfgSetRadius
615 *      - changes the value of radius for a config to the given value
616 * @param config
617 *      - config that is being altered
618 * @param radius
619 *      - new radius
620 */
621PRIVATE void
622cfgSetRadius(config       *cfg,
623             const double radius)
624{
625  cfg->radius = radius;
626}
627
628
629/**
630 * @brief cfgUpdateMinRadius
631 *      - updates the minimum possible radius for the given config
632 * @param config
633 * @param unpaired
634 * @param paired
635 */
636PRIVATE void
637cfgUpdateMinRadius(config       *cfg,
638                   const double unpaired,
639                   const double paired)
640{
641  double minRadius = approximateConfigRadius(cfg, unpaired, paired);
642
643  cfg->minRadius = minRadius;
644}
645
646
647PRIVATE double
648cfgApplyChanges(config                            *cfg,
649                const char                        loopName,
650                const double                      *deltaCfg,
651                const double                      radiusNew,
652                const vrna_plot_options_puzzler_t *puzzler)
653{
654  /* start with adjusting config angles; if applicable */
655  if (deltaCfg != NULL) {
656    for (int currentArc = 0; currentArc < cfg->numberOfArcs; currentArc++)
657
658      (cfg->cfgArcs[currentArc]).arcAngle += deltaCfg[currentArc];
659  }
660
661  /* then, adjust config radius */
662  double  oldRadius = cfg->radius;
663  double  newRadius = -1.0;
664  if (radiusNew > 0.0) {
665    /*
666     * in case the input is a positive value
667     * we set the minimum of valid and input as new radius
668     */
669    cfgUpdateMinRadius(cfg, puzzler->unpaired, puzzler->paired);
670    newRadius = fmax(radiusNew, cfg->minRadius);
671    cfgSetRadius(cfg, newRadius);
672  } else if (radiusNew == 0.0) {
673    /*
674     * set the minRadius as new value
675     * (this allows to shrink a loop)
676     */
677    cfgUpdateMinRadius(cfg, puzzler->unpaired, puzzler->paired);
678    newRadius = cfg->minRadius;
679    cfgSetRadius(cfg, newRadius);
680  } else if (radiusNew == -1.0) {
681    /*
682     * set the minRadius as new value
683     * (this forbidds to shrink a loop)
684     */
685    cfgUpdateMinRadius(cfg, puzzler->unpaired, puzzler->paired);
686
687    if (cfg->minRadius - EPSILON_0 > oldRadius) {
688      newRadius = cfg->minRadius;
689    } else {
690      double defaultIncrease = 1.05;
691      newRadius = oldRadius * defaultIncrease;
692    }
693
694    cfgSetRadius(cfg, newRadius);
695  } else {
696    /* all unhandled inputs result in errors */
697    newRadius = -1.0;
698  }
699
700  return newRadius;
701}
702
703
704PRIVATE short
705cfgIsValid(config       *cfg,
706           const double *deltaCfg)
707{
708  if (deltaCfg == NULL)
709
710    return 0;
711
712  double  sumAngles         = 0.0;
713  short   validSingleAngles = 1;
714
715  for (int currentArc = 0; currentArc < cfg->numberOfArcs; currentArc++) {
716    double  angle = getArcAngle(cfg, currentArc) + deltaCfg[currentArc];
717    sumAngles += angle;
718
719    short   validAngle = 0.0 < angle && angle < MATH_TWO_PI;
720    validSingleAngles = validSingleAngles && validAngle;
721  }
722
723  short validSumAngles = (fabs(sumAngles - MATH_TWO_PI) < EPSILON_3);
724
725  return validSingleAngles && validSumAngles;
726}
727
728
729#endif
730