1#ifndef RNAPUZZLER_RESOLVE_EXT_CHILD_INTERSECT
2#define RNAPUZZLER_RESOLVE_EXT_CHILD_INTERSECT
3
4#include <stdlib.h>
5
6#include "ViennaRNA/utils/basic.h"
7
8#include "definitions.inc"
9#include "../headers/configtree_struct.h"
10#include "../headers/tBaseInformation_struct.h"
11#include "configtree.inc"
12#include "boundingBoxes.inc"
13#include "intersectLevelTreeNodes.inc"
14
15
16PRIVATE void
17resolveExteriorChildrenIntersectionXY(treeNode            *exteriorNode,
18                                      short const *const  pair_table,
19                                      const double        unpaired,
20                                      const short         allowFlipping,
21                                      double              *myX,
22                                      double              *myY);
23
24
25PRIVATE void
26resolveExteriorChildrenIntersectionAffin(treeNode                 *exteriorNode,
27                                         short const *const       pair_table,
28                                         tBaseInformation *const  baseInformation,
29                                         const double             unpaired,
30                                         const short              allowFlipping);
31
32
33PRIVATE void
34resolveExteriorChildIntersections(treeNode                *exteriorNode,
35                                  short const *const      pair_table,
36                                  tBaseInformation *const baseInformation,
37                                  const double            unpaired,
38                                  const short             allowFlipping);
39
40
41/* ----------------------------------------------------------------------------- */
42
43PRIVATE void
44getSimpleBoundingBox(treeNode       *node,
45                     double *const  bounds,
46                     const int      recursionDepth)
47{
48  double  loopMin = node->lBox->c[0] - node->lBox->r;
49  double  loopMax = node->lBox->c[0] + node->lBox->r;
50
51  if (recursionDepth == 0) {
52    bounds[0] = loopMin;
53    bounds[1] = loopMax;
54  }
55
56  for (int i = 0; i < node->childCount; i++) {
57    treeNode *child = getChild(node, i);
58    getSimpleBoundingBox(child, bounds, recursionDepth + 1);
59  }
60
61  if (loopMin < bounds[0])
62    bounds[0] = loopMin;
63
64  if (loopMax > bounds[1])
65    bounds[1] = loopMax;
66
67  for (int i = 0; i < node->sBox->bulgeCount; i++) {
68    double  pPrev[2];
69    double  pThis[2];
70    double  pNext[2];
71    getBulgeCoordinates(node->sBox, i, pPrev, pThis, pNext);
72
73    if (pThis[0] < bounds[0])
74      bounds[0] = pThis[0];
75
76    if (pThis[0] > bounds[1])
77      bounds[1] = pThis[0];
78  }
79}
80
81
82/**
83 * Resolve the intersections of the children of the exterior loop
84 */
85PRIVATE void
86resolveExteriorChildrenIntersectionXY(treeNode            *exteriorNode,
87                                      short const *const  pair_table,
88                                      const double        unpaired,
89                                      const short         allowFlipping,
90                                      double              *myX,
91                                      double              *myY)
92{
93  /* number of subtrees */
94  int subtreeCount = exteriorNode->childCount;
95
96  if (subtreeCount < 2)
97    return;
98
99  /* for each exterior child: get first node */
100  treeNode **childTreeNode = (treeNode **)vrna_alloc(subtreeCount * sizeof(treeNode *));
101  for (int subtree = 0; subtree < subtreeCount; subtree++)
102    childTreeNode[subtree] = getChild(exteriorNode, subtree);
103
104#if 0
105  /* for each exterior child: prepare bounding box */
106  double **bounds = (double **)vrna_alloc(subtreeCount * sizeof(double *));
107  for (int subtree = 0; subtree < subtreeCount; subtree++) {
108    bounds[subtree]     = (double *)vrna_alloc(2 * sizeof(double));
109    bounds[subtree][0]  = 0.0;
110    bounds[subtree][1]  = 0.0;
111  }
112
113  /* get bounding box of first child */
114  getSimpleBoundingBox(childTreeNode[0], bounds[0], 0);
115#endif
116
117  /*
118   * for each subtree
119   * - compute number of its first non-exterior base
120   * - compute number of nucleotides before the subtree
121   * - distance between nucleotides before the subtree
122   */
123  int     *firstBase  = (int *)vrna_alloc(subtreeCount * sizeof(int));
124  int     *backbone   = (int *)vrna_alloc(subtreeCount * sizeof(int));
125  double  *distance   = (double *)vrna_alloc(subtreeCount * sizeof(double));
126  for (int subtree = 0; subtree < subtreeCount; subtree++) {
127    backbone[subtree] = 0;
128    distance[subtree] = 0.0;
129  }
130  int     subtree = 0;
131  int     base    = 1;
132  while (base < pair_table[0] && subtree < subtreeCount) {
133    if (pair_table[base] > base) {
134      firstBase[subtree] = base;
135      subtree++;
136      base = pair_table[base];
137    } else {
138      base++;
139      backbone[subtree]++;
140    }
141  }
142
143  /* store upper and lower subtrees */
144  int *upper  = (int *)vrna_alloc((subtreeCount + 1) * sizeof(int));
145  int *lower  = (int *)vrna_alloc((subtreeCount + 1) * sizeof(int));
146  upper[0]  = 0;
147  lower[0]  = 0;
148
149  /* set first subtree to upper side */
150  upper[0]++;
151  upper[upper[0]] = 0;
152
153  /* accumulated offset of children */
154  double  offset = 0.0;
155  /* accumulated translation of children */
156  double  accumulatedTranslation = 0.0;
157
158  /* for all subtrees */
159  for (int subtree = 1; subtree < subtreeCount; subtree++) {
160    /* translate current subtree by accumulated offset */
161    if (offset > 0.0) {
162      double translate[2] = {
163        offset, 0.0
164      };
165      translateBoundingBoxes(childTreeNode[subtree], translate);
166    }
167
168    /* getSimpleBoundingBox(childTreeNode[subtree], bounds[subtree], 0); */
169
170    /* as long as the current child gets translated */
171    short   changed         = 1;
172    short   intersectUpper  = 0;
173    short   intersectLower  = 0;
174    double  fixOverlap      = 0.0;
175    while (changed) {
176      /* printf("Handling subtree %d\n", subtree); */
177      changed         = 0;
178      intersectUpper  = 0;
179      intersectLower  = 0;
180
181      /* check intersection of current subtree with previous upper subtrees */
182      for (int u = 1; u <= upper[0]; u++) {
183        int upperStem = upper[u];
184        /* printf("Check intersection between %d and %d\n", subtree, upperStem); */
185        intersectUpper = intersectTrees(childTreeNode[subtree], childTreeNode[upperStem]);
186        if (intersectUpper)
187          /* printf("Intersection between %d and %d\n", subtree, upperStem); */
188          break;
189
190        /* printf("%d vs %d: upperOverlap:%f (boundsOverlap:%f)\n", subtree, upperStem, upperOverlap, boundsOverlap); */
191      }
192
193      if (allowFlipping) {
194        /*
195         * if flipping is allowed:
196         * check intersection of current subtree with previous lower subtrees
197         */
198        for (int l = 1; l <= lower[0]; l++) {
199          int lowerStem = lower[l];
200          intersectLower = intersectTrees(childTreeNode[subtree], childTreeNode[lowerStem]);
201          if (intersectLower)
202            break;
203
204          /* printf("%d vs %d: lowerOverlap:%f (boundsOverlap:%f)\n", subtree, lowerStem, lowerOverlap, boundsOverlap); */
205        }
206      }
207
208      if ((!allowFlipping && intersectUpper) ||
209          (allowFlipping && intersectUpper && intersectLower)) {
210        /*
211         * if intersections can not be resolved by flipping
212         * increase distance by constant amount per exterior base
213         */
214        distance[subtree] += unpaired;         /* minOverlap / backbone[subtree]; */
215        fixOverlap        = unpaired * backbone[subtree];
216        /* printf("Increase distance by %12.8lf\n", fixOverlap); */
217
218        double translate[2] = {
219          fixOverlap, 0.0
220        };
221        translateBoundingBoxes(childTreeNode[subtree], translate);
222
223        /*
224         * bounds[subtree][0] += fixOverlap;
225         * bounds[subtree][1] += fixOverlap;
226         */
227
228        offset += fixOverlap;
229        /* printf("Total offset: %12.8lf\n", offset); */
230
231        changed = 1;
232      } else {
233        if (allowFlipping && intersectUpper) {
234          /*
235           * if flipping is allowed and sufficient for resolving the intersection:
236           * (intersection is on the upper side)
237           */
238          lower[0]++;
239          lower[lower[0]] = subtree;
240        } else {
241          upper[0]++;
242          upper[upper[0]] = subtree;
243        }
244      }
245    }     /* end while(changed) */
246
247    /* translate exterior bases between previous and current subtree */
248    int currentBase = 1;
249    for (int base = pair_table[firstBase[subtree - 1]];
250         base < firstBase[subtree];
251         base++, ++currentBase)
252      myX[base] += currentBase * distance[subtree] + accumulatedTranslation;
253    accumulatedTranslation += distance[subtree] * backbone[subtree];
254  }
255
256  /* Last part of the exterior loop */
257  for (int base = pair_table[firstBase[subtreeCount - 1]];
258       base < pair_table[0];
259       ++base)
260    myX[base] += accumulatedTranslation;
261
262  /* modify x- and y-coordinates for all subtrees */
263  int     currentLower  = 1;
264  double  translation   = 0.0;
265  for (int subtree = 1; subtree < subtreeCount; subtree++) {
266    /* translate all bases of current subtree */
267    translation += distance[subtree] * backbone[subtree];
268    /* printf("Translate subtree %d by %12.8lf\n", subtree, translation); */
269    for (int base = firstBase[subtree];
270         base < pair_table[firstBase[subtree]];
271         ++base)
272      myX[base] += translation;
273
274    if (subtree == lower[currentLower]) {
275      /* flip subtrees */
276      double exteriorY = myY[1];
277      for (int base = firstBase[subtree];
278           base < pair_table[firstBase[subtree]];
279           ++base)
280        myY[base] = 2 * exteriorY - myY[base];
281      ++currentLower;
282    }
283  }
284  /* processing end */
285
286  free(upper);
287  free(lower);
288  free(backbone);
289  free(distance);
290  free(firstBase);
291  /*
292   * for (int subtree = 0; subtree < subtreeCount; subtree++) {
293   *  free(bounds[subtree]);
294   * }
295   * free(bounds);
296   */
297  free(childTreeNode);
298}
299
300
301/**
302 * Resolve the intersections of the children of the exterior loop
303 */
304PRIVATE void
305resolveExteriorChildrenIntersectionAffin(treeNode                 *exteriorNode,
306                                         short const *const       pair_table,
307                                         tBaseInformation *const  baseInformation,
308                                         const double             unpaired,
309                                         const short              allowFlipping)
310{
311  /* number of subtrees */
312  int subtreeCount = exteriorNode->childCount;
313
314  if (subtreeCount < 2)
315    return;
316
317  /* for each exterior child: get first node */
318  treeNode **childTreeNode = (treeNode **)vrna_alloc(subtreeCount * sizeof(treeNode *));
319  for (int subtree = 0; subtree < subtreeCount; subtree++)
320    childTreeNode[subtree] = getChild(exteriorNode, subtree);
321
322#if 0
323  /* for each exterior child: prepare bounding box */
324  double **bounds = (double **)vrna_alloc(subtreeCount * sizeof(double *));
325  for (int subtree = 0; subtree < subtreeCount; subtree++) {
326    bounds[subtree]     = (double *)vrna_alloc(2 * sizeof(double));
327    bounds[subtree][0]  = 0.0;
328    bounds[subtree][1]  = 0.0;
329  }
330
331  /* get bounding box of first child */
332  getSimpleBoundingBox(childTreeNode[0], bounds[0], 0);
333#endif
334
335  /*
336   * for each subtree
337   * - compute number of its first non-exterior base
338   * - compute number of nucleotides before the subtree
339   */
340  int *firstBase  = (int *)vrna_alloc(subtreeCount * sizeof(int));
341  int *backbone   = (int *)vrna_alloc(subtreeCount * sizeof(int));
342  for (int subtree = 0; subtree < subtreeCount; subtree++)
343    backbone[subtree] = 0;
344  int subtree = 0;
345  int base    = 1;
346  while (base < pair_table[0] && subtree < subtreeCount) {
347    if (pair_table[base] > base) {
348      firstBase[subtree] = base;
349      subtree++;
350      base = pair_table[base];
351    } else {
352      base++;
353      backbone[subtree]++;
354    }
355  }
356
357  /* store upper and lower stems */
358  int *upper  = (int *)vrna_alloc((subtreeCount + 1) * sizeof(int));
359  int *lower  = (int *)vrna_alloc((subtreeCount + 1) * sizeof(int));
360  upper[0]  = 0;
361  lower[0]  = 0;
362
363  /* set first subtree to upper side */
364  upper[0]++;
365  upper[upper[0]] = 0;
366
367  /* accumulated offset of children */
368  double offset = 0.0;
369
370  /* for all stems */
371  for (int subtree = 1; subtree < subtreeCount; subtree++) {
372    /* translate current subtree by accumulated offset */
373    if (offset > 0.0) {
374      double translate[2] = {
375        offset, 0.0
376      };
377      translateBoundingBoxes(childTreeNode[subtree], translate);
378    }
379
380    /* getSimpleBoundingBox(childTreeNode[subtree], bounds[subtree], 0); */
381
382    /* as long as the current child gets translated */
383    short   changed         = 1;
384    short   intersectUpper  = 0;
385    short   intersectLower  = 0;
386    double  fixOverlap      = 0.0;
387    while (changed) {
388      changed         = 0;
389      intersectUpper  = 0;
390      intersectLower  = 0;
391
392      /* check intersection of current subtree with previous upper stems */
393      for (int u = 1; u <= upper[0]; u++) {
394        int upperStem = upper[u];
395        intersectUpper = intersectTrees(childTreeNode[subtree], childTreeNode[upperStem]);
396        if (intersectUpper)
397          break;
398
399        /* printf("%d vs %d: upperOverlap:%f (boundsOverlap:%f)\n", subtree, upperStem, upperOverlap, boundsOverlap); */
400      }
401
402      if (allowFlipping) {
403        /*
404         * if flipping is allowed:
405         * check intersection of current subtree with previous lower stems
406         */
407        for (int l = 1; l <= lower[0]; l++) {
408          int lowerStem = lower[l];
409          intersectLower = intersectTrees(childTreeNode[subtree], childTreeNode[lowerStem]);
410          if (intersectLower)
411            break;
412
413          /* printf("%d vs %d: lowerOverlap:%f (boundsOverlap:%f)\n", subtree, lowerStem, lowerOverlap, boundsOverlap); */
414        }
415      }
416
417      if ((!allowFlipping && intersectUpper) ||
418          (allowFlipping && intersectUpper && intersectLower)) {
419        /*
420         * if intersections can not be resolved by flipping
421         * increase distance by constant amount per exterior base
422         */
423        double deltaPerDistance = unpaired;         /* minOverlap / backbone[subtree]; */
424        fixOverlap = deltaPerDistance * backbone[subtree];
425
426        for (int base = pair_table[firstBase[subtree - 1]]; base < firstBase[subtree]; base++)
427          baseInformation[base].distance += deltaPerDistance;
428
429        double translate[2] = {
430          fixOverlap, 0.0
431        };
432        translateBoundingBoxes(childTreeNode[subtree], translate);
433
434        /*
435         * bounds[subtree][0] += fixOverlap;
436         * bounds[subtree][1] += fixOverlap;
437         */
438
439        offset += fixOverlap;
440
441        changed = 1;
442      } else {
443        if (allowFlipping && intersectUpper) {
444          /*
445           * if flipping is allowed and sufficient for resolving the intersection:
446           * (intersection is on the upper side)
447           */
448          for (int base = firstBase[subtree] + 1; base <= pair_table[firstBase[subtree]] + 1;
449               base++) {
450            if (base > pair_table[0])
451              break;
452
453            baseInformation[base].angle *= -1;
454          }
455
456          lower[0]++;
457          lower[lower[0]] = subtree;
458        } else {
459          upper[0]++;
460          upper[upper[0]] = subtree;
461        }
462      }
463    }     /* end while(changed) */
464  }
465  /* processing end */
466
467  free(upper);
468  free(lower);
469  free(backbone);
470  free(firstBase);
471  /*
472   * for (int subtree = 0; subtree < subtreeCount; subtree++) {
473   *  free(bounds[subtree]);
474   * }
475   * free(bounds);
476   */
477  free(childTreeNode);
478}
479
480
481PRIVATE void
482resolveExteriorChildIntersections(treeNode                *exteriorNode,
483                                  short const *const      pair_table,
484                                  tBaseInformation *const baseInformation,
485                                  const double            unpaired,
486                                  const short             allowFlipping)
487{
488  int stemCount = exteriorNode->childCount;
489
490  if (stemCount < 2)
491    return;
492
493  /* for each exterior child: get first node */
494  treeNode  **node = (treeNode **)vrna_alloc(stemCount * sizeof(treeNode *));
495  for (int stem = 0; stem < stemCount; stem++)
496    node[stem] = getChild(exteriorNode, stem);
497
498  /* for each exterior child: prepare bounding box */
499  double    **bounds = (double **)vrna_alloc(stemCount * sizeof(double *));
500  for (int stem = 0; stem < stemCount; stem++) {
501    bounds[stem]    = (double *)vrna_alloc(2 * sizeof(double));
502    bounds[stem][0] = 0.0;
503    bounds[stem][1] = 0.0;
504  }
505  getSimpleBoundingBox(node[0], bounds[0], 0);
506
507  /*
508   * for each stem
509   * - compute number of its first base
510   * - compute number of nucleotides before the stem
511   */
512  int *firstBase  = (int *)vrna_alloc(stemCount * sizeof(int));
513  int *backbone   = (int *)vrna_alloc(stemCount * sizeof(int));
514  for (int stem = 0; stem < stemCount; stem++)
515    backbone[stem] = 0;
516  int stem  = 0;
517  int base  = 1;
518  while (base < pair_table[0] && stem < stemCount) {
519    if (pair_table[base] > base) {
520      firstBase[stem] = base;
521      stem++;
522      base = pair_table[base];
523    } else {
524      base++;
525      backbone[stem]++;
526    }
527  }
528
529  /* store upper and lower stems */
530  int *upper  = (int *)vrna_alloc((stemCount + 1) * sizeof(int));
531  int *lower  = (int *)vrna_alloc((stemCount + 1) * sizeof(int));
532  upper[0]  = 0;
533  lower[0]  = 0;
534
535  /* set first stem to upper side */
536  upper[0]++;
537  upper[upper[0]] = 0;
538
539  /* accumulated offset of children */
540  double offset = 0.0;
541
542  /* processing start */
543  for (int stem = 1; stem < stemCount; stem++) {
544    /* initial setting */
545    if (offset > 0.0) {
546      double translate[2] = {
547        offset, 0.0
548      };
549      translateBoundingBoxes(node[stem], translate);
550    }
551
552    getSimpleBoundingBox(node[stem], bounds[stem], 0);
553
554    /* as long as the current child gets translated */
555    short changed = 1;
556    while (changed) {
557      changed = 0;
558
559      /* get overlap with upper and lower side */
560      double upperOverlap = 0.0;
561      for (int u = 1; u <= upper[0]; u++) {
562        int     upperStem     = upper[u];
563        double  boundsOverlap = (bounds[upperStem][1] + unpaired) - bounds[stem][0];
564        if (boundsOverlap > upperOverlap && intersectTrees(node[stem], node[upperStem]))
565          upperOverlap = boundsOverlap;
566
567        /* printf("%d vs %d: upperOverlap:%f (boundsOverlap:%f)\n", stem, upperStem, upperOverlap, boundsOverlap); */
568      }
569      double lowerOverlap = 0.0;
570      for (int l = 1; l <= lower[0]; l++) {
571        int     lowerStem     = lower[l];
572        double  boundsOverlap = (bounds[lowerStem][1] + unpaired) - bounds[stem][0];
573        if (boundsOverlap > lowerOverlap && intersectTrees(node[stem], node[lowerStem]))
574          lowerOverlap = boundsOverlap;
575
576        /* printf("%d vs %d: lowerOverlap:%f (boundsOverlap:%f)\n", stem, lowerStem, lowerOverlap, boundsOverlap); */
577      }
578
579      /* fix minimum of both overlaps */
580      double minOverlap =
581        (allowFlipping && (lowerOverlap < upperOverlap)) ? lowerOverlap : upperOverlap;
582      if (minOverlap > 0.0) {
583        double  deltaPerDistance = unpaired;        /* minOverlap / backbone[stem]; */
584        minOverlap = deltaPerDistance * backbone[stem];
585        for (int base = pair_table[firstBase[stem - 1]]; base < firstBase[stem]; base++)
586          baseInformation[base].distance += deltaPerDistance;
587
588        double  translate[2] = {
589          minOverlap, 0.0
590        };
591        translateBoundingBoxes(node[stem], translate);
592
593        bounds[stem][0] += minOverlap;
594        bounds[stem][1] += minOverlap;
595
596        offset += minOverlap;
597
598        changed = 1;
599      } else {
600        /* flip if necessary */
601        if (lowerOverlap < upperOverlap) {
602          for (int base = firstBase[stem] + 1; base <= pair_table[firstBase[stem]] + 1; base++) {
603            if (base > pair_table[0])
604              break;
605
606            baseInformation[base].angle *= -1;
607          }
608
609          lower[0]++;
610          lower[lower[0]] = stem;
611        } else {
612          upper[0]++;
613          upper[upper[0]] = stem;
614        }
615      }
616    }     /* end while(changed) */
617  }
618  /* processing end */
619
620  free(upper);
621  free(lower);
622  free(backbone);
623  free(firstBase);
624  for (int stem = 0; stem < stemCount; stem++)
625    free(bounds[stem]);
626  free(bounds);
627  free(node);
628}
629
630
631#endif
632