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