1 /*
2 * NAVIEW -- A program to make a modified radial drawing of an RNA
3 * secondary structure.
4 *
5 * Copyright (c) 1988 Robert E. Bruccoleri
6 * Copying of this software, in whole or in part, is permitted
7 * provided that the copies are not made for commercial purposes,
8 * appropriate credit for the use of the software is given, this
9 * copyright notice appears, and notice is given that the copying
10 * is by permission of Robert E. Bruccoleri. Any other copying
11 * requires specific permission.
12 *
13 * See R. Bruccoleri and G. Heinrich, Computer Applications in the
14 * Biosciences, 4, 167-173 (1988) for a full description.
15 *
16 * In November 1997, Michael Zuker made a number of changes to bring
17 * naview up to modern standards. All functions defined in naview are
18 * now declared before main() with arguments and argument types.
19 * When functions are defined, their argument types are declared
20 * with the function and these definitions are removed after the '{'.
21 * The 'void' declaration was used as necessary for functions.
22 *
23 * The troublesome na_scanf function was deleted and replaced by
24 * scanf. Finally, there is now no default for the minimum separation
25 * of bases. A floating point number must be entered. However, as
26 * before an entry < 0 will be moved up to 0 and an entry > 0.5
27 * will be reduced to 0.5.
28 *
29 * Adapted for use as a subroutine in the Vienna RNA Package
30 * by Ivo Hofacker, May 1998:
31 * deleted output routines, replaced main() by naview_xy_coordinates(),
32 * which fills the X and Y arrays used by PS_rna_plot() etc.
33 * added ansi prototypes and fixed memory leaks.
34 */
35
36 #ifdef HAVE_CONFIG_H
37 #include "config.h"
38 #endif
39
40 #include <stdio.h>
41 #include <stdlib.h>
42 #include <string.h>
43 #include <math.h>
44
45 #include "ViennaRNA/utils/basic.h"
46 #include "ViennaRNA/plotting/naview.h"
47
48 typedef int LOGICAL;
49 #define logical LOGICAL
50
51 #define true 1
52 #define false 0
53 #define FATAL_ERROR 1
54 #define SUCCESS 0
55
56 #define type_alloc(type) (type *) vrna_alloc(sizeof(type))
57
58 #define struct_alloc(structure_name) type_alloc(struct structure_name)
59
60 #define add_double_list(head,tail,newp) {\
61 (newp)->next = (newp)->prev = NULL; \
62 if ((head) == NULL) (head) = (tail) = (newp); \
63 else { \
64 (tail)->next = (newp); \
65 (newp)->prev = (tail); \
66 (tail) = (newp); \
67 } \
68 }
69
70 static double pi = 3.141592653589793;
71 static double anum = 9999.0;
72
73
74
75 /*
76 * Function data type definitions
77 */
78
79 #define minf2(x1, x2) ((x1)<(x2))?(x1):(x2)
80 #define maxf2(x1, x2) ((x1)>(x2))?(x1):(x2)
81
82 static struct base {
83 int mate;
84 double x,y;
85 logical extracted;
86 struct region *region;
87 } *bases;
88
89 struct region {
90 int start1,end1,start2,end2;
91 };
92
93 struct loop {
94 int nconnection;
95 struct connection **connections;
96 int number;
97 int depth;
98 logical mark;
99 double x,y,radius;
100 };
101
102 struct connection {
103 struct loop *loop;
104 struct region *region;
105 int start,end; /* Start and end form the 1st base pair of the region. */
106 double xrad,yrad,angle;
107 logical extruded; /* True if segment between this connection and
108 the next must be extruded out of the circle */
109 logical broken; /* True if the extruded segment must be drawn long. */
110 };
111
112 static int nbase, nregion, loop_count;
113
114 static struct loop *root, *loops;
115
116 static struct region *regions;
117
118 static struct loop *construct_loop(int ibase);
119
120 struct radloop {
121 double radius;
122 int loopnumber;
123 struct radloop *next, *prev;
124 };
125
126 static struct radloop *rlphead;
127
128 static double lencut;
129
130 static logical debug = false;
131
132 static void read_in_bases(const short *pair_table);
133 static void find_regions(void);
134 static void dump_loops(void);
135 static void find_central_loop(void);
136 static void determine_depths(void);
137 static void traverse_loop(struct loop *lp,struct connection *anchor_connection);
138 static void determine_radius(struct loop *lp,double lencut);
139 static void generate_region(struct connection *cp);
140 static void construct_extruded_segment(struct connection *cp,struct connection *cpnext);
141 static void find_center_for_arc(int n,double b,double *hp,double *thetap);
142 static int depth(struct loop *lp);
143
144 static logical connected_connection(struct connection *cp, struct connection *cpnext);
145 static int find_ic_middle(int icstart, int icend, struct connection *anchor_connection, struct connection *acp, struct loop *lp);
146
147
148
149 int
vrna_plot_coords_naview(const char * structure,float ** x,float ** y)150 vrna_plot_coords_naview(const char *structure,
151 float **x,
152 float **y)
153 {
154 if (structure) {
155 int ret = 0;
156 short *pt = vrna_ptable(structure);
157
158 ret = vrna_plot_coords_naview_pt(pt, x, y);
159
160 free(pt);
161
162 return ret;
163 }
164
165 if (x)
166 *x = NULL;
167
168 if (y)
169 *y = NULL;
170
171 return 0;
172 }
173
174
175 int
vrna_plot_coords_naview_pt(const short * pt,float ** x,float ** y)176 vrna_plot_coords_naview_pt(const short *pt,
177 float **x,
178 float **y)
179 {
180 int i;
181
182 if ((pt) && (x) && (y)) {
183 nbase = pt[0]; /* length */
184 *x = (float *)vrna_alloc(sizeof(float) * (nbase + 1));
185 *y = (float *)vrna_alloc(sizeof(float) * (nbase + 1));
186 bases = (struct base *) vrna_alloc(sizeof(struct base)*(nbase+1));
187 regions = (struct region *) vrna_alloc(sizeof(struct region)*(nbase+1));
188 loops = (struct loop *) vrna_alloc(sizeof(struct loop) * (nbase + 1));
189 lencut = 0.5;
190 rlphead = NULL;
191 loop_count = 0;
192
193 read_in_bases(pt);
194
195 find_regions();
196
197 construct_loop(0);
198
199 find_central_loop();
200
201 if (debug)
202 dump_loops();
203
204 traverse_loop(root, NULL);
205
206 for (i = 0; i < nbase; i++) {
207 (*x)[i] = 100 + 15 * bases[i + 1].x;
208 (*y)[i] = 100 + 15 * bases[i + 1].y;
209 }
210 free(bases);
211 free(regions);
212 free(loops);
213
214 return nbase;
215 }
216
217 if (x)
218 *x = NULL;
219
220 if (y)
221 *y = NULL;
222
223 return 0;
224 }
225
226
naview_xy_coordinates(short * pair_table,float * X,float * Y)227 int naview_xy_coordinates(short *pair_table, float *X, float *Y) {
228 int i;
229
230 nbase = pair_table[0]; /* length */
231 bases = (struct base *) vrna_alloc(sizeof(struct base)*(nbase+1));
232 regions = (struct region *) vrna_alloc(sizeof(struct region)*(nbase+1));
233 read_in_bases(pair_table);
234 lencut = 0.5;
235 rlphead = NULL;
236 find_regions();
237 loop_count = 0;
238 loops = (struct loop *) vrna_alloc(sizeof(struct loop)*(nbase+1));
239 construct_loop(0);
240 find_central_loop();
241 if (debug) dump_loops();
242
243 traverse_loop(root,NULL);
244 for (i=0; i<nbase; i++) {
245 X[i] = 100 + 15*bases[i+1].x;
246 Y[i] = 100 + 15*bases[i+1].y;
247 }
248 free(bases);
249 free(regions);
250 free(loops);
251 return nbase;
252 }
253
254
read_in_bases(const short * pair_table)255 static void read_in_bases(const short *pair_table)
256 {
257 int i,npairs;
258
259 /* Set up an origin. */
260 bases[0].mate = 0;
261 bases[0].extracted = false;
262 bases[0].x = anum;
263 bases[0].y = anum;
264
265 for (npairs=0,i=1; i<=nbase; i++) {
266 bases[i].extracted = false;
267 bases[i].x = anum;
268 bases[i].y = anum;
269 bases[i].mate = pair_table[i];
270 if ((int) pair_table[i]>i) npairs++;
271 }
272 if (npairs==0) { /* must have at least 1 pair to avoid segfault */
273 bases[1].mate=nbase;
274 bases[nbase].mate=1;
275 }
276 }
277
278
find_regions(void)279 static void find_regions(void)
280 /*
281 * Identifies the regions in the structure.
282 */
283
284 {
285 int i,mate,nb1;
286 logical *mark;
287
288 nb1 = nbase + 1;
289 mark = (int *) vrna_alloc(sizeof(int)*nb1);
290 for (i = 0; i < nb1; i++) mark[i] = false;
291 nregion = 0;
292 for (i=0; i<=nbase; i++) {
293 if ( (mate = bases[i].mate) && !mark[i]) {
294 regions[nregion].start1 = i;
295 regions[nregion].end2 = mate;
296 mark[i] = true;
297 mark[mate] = true;
298 bases[i].region = bases[mate].region = ®ions[nregion];
299 for (i++,mate--;
300 i<mate && bases[i].mate == mate;
301 i++,mate--) {
302 mark[i] = mark[mate] = true;
303 bases[i].region = bases[mate].region = ®ions[nregion];
304 }
305 regions[nregion].end1 = --i;
306 regions[nregion].start2 = mate+1;
307 if (debug) {
308 if (nregion == 0) printf("\nRegions are:\n");
309 printf("Region %d is %d-%d and %d-%d with gap of %d.\n",
310 nregion+1,regions[nregion].start1,regions[nregion].end1,
311 regions[nregion].start2,regions[nregion].end2,
312 regions[nregion].start2-regions[nregion].end1+1);
313 }
314 nregion++;
315 }
316 }
317 free(mark);
318 }
319
320
construct_loop(int ibase)321 static struct loop *construct_loop(int ibase)
322 /*
323 * Starting at residue ibase, recursively constructs the loop containing
324 * said base and all deeper bases.
325 */
326
327 {
328 int i,mate;
329 struct loop *retloop,*lp;
330 struct connection *cp;
331 struct region *rp;
332 struct radloop *rlp;
333
334 retloop = &loops[loop_count++];
335 retloop->nconnection = 0;
336 retloop->connections = (struct connection **) vrna_alloc(sizeof(struct connection *));
337 retloop->depth = 0;
338 retloop->number = loop_count;
339 retloop->radius = 0.0;
340 for (rlp = rlphead; rlp; rlp = rlp->next)
341 if (rlp->loopnumber == loop_count) retloop->radius = rlp->radius;
342 i = ibase;
343 do {
344 if ((mate = bases[i].mate) != 0) {
345 rp = bases[i].region;
346 if (!bases[rp->start1].extracted) {
347 if (i == rp->start1) {
348 bases[rp->start1].extracted = true;
349 bases[rp->end1].extracted = true;
350 bases[rp->start2].extracted = true;
351 bases[rp->end2].extracted = true;
352 lp = construct_loop(rp->end1 < nbase ? rp->end1+1 : 0);
353 }
354 else if (i == rp->start2){
355 bases[rp->start2].extracted = true;
356 bases[rp->end2].extracted = true;
357 bases[rp->start1].extracted = true;
358 bases[rp->end1].extracted = true;
359 lp = construct_loop(rp->end2 < nbase ? rp->end2+1 : 0);
360 }
361 else {
362 vrna_message_error("naview: Error detected in construct_loop. i = %d not found in region table.",i);
363 exit(FATAL_ERROR);
364 }
365 retloop->connections = (struct connection **)
366 realloc(retloop->connections,
367 (++retloop->nconnection+1) *
368 sizeof(struct connection *));
369 retloop->connections[retloop->nconnection-1] = cp =
370 struct_alloc(connection);
371 retloop->connections[retloop->nconnection] = NULL;
372 cp->loop = lp;
373 cp->region = rp;
374 if (i == rp->start1) {
375 cp->start = rp->start1;
376 cp->end = rp->end2;
377 }
378 else {
379 cp->start = rp->start2;
380 cp->end = rp->end1;
381 }
382 cp->extruded = false;
383 cp->broken = false;
384 lp->connections = (struct connection **)
385 realloc(lp->connections,
386 (++lp->nconnection+1) *
387 sizeof(struct connection *));
388 lp->connections[lp->nconnection-1] = cp =
389 struct_alloc(connection);
390 lp->connections[lp->nconnection] = NULL;
391 cp->loop = retloop;
392 cp->region = rp;
393 if (i == rp->start1) {
394 cp->start = rp->start2;
395 cp->end = rp->end1;
396 }
397 else {
398 cp->start = rp->start1;
399 cp->end = rp->end2;
400 }
401 cp->extruded = false;
402 cp->broken = false;
403 }
404 i = mate;
405 }
406 if (++i > nbase) i = 0;
407 }
408 while (i != ibase);
409 return retloop;
410 }
411
412
dump_loops(void)413 static void dump_loops(void)
414 /*
415 * Displays all the loops.
416 */
417
418 {
419 int il,ilp,irp;
420 struct loop *lp;
421 struct connection *cp,**cpp;
422
423 printf("\nRoot loop is #%d\n",(int)(root-loops)+1);
424 for (il=0; il < loop_count; il++) {
425 lp = &loops[il];
426 printf("Loop %d has %d connections:\n",il+1,lp->nconnection);
427 for (cpp = lp->connections; (cp = *cpp); cpp++) {
428 ilp = (cp->loop - loops) + 1;
429 irp = (cp->region - regions) + 1;
430 printf(" Loop %d Region %d (%d-%d)\n",
431 ilp,irp,cp->start,cp->end);
432 }
433 }
434 }
435
436
find_central_loop(void)437 static void find_central_loop(void)
438 /*
439 * Find node of greatest branching that is deepest.
440 */
441
442 {
443 struct loop *lp;
444 int maxconn,maxdepth,i;
445
446 determine_depths();
447 maxconn = 0;
448 maxdepth = -1;
449
450 for (i=0; i<loop_count; i++) {
451 lp = &loops[i];
452 if (lp->nconnection > maxconn) {
453 maxdepth = lp->depth;
454 maxconn = lp->nconnection;
455 root = lp;
456 }
457 else if (lp->depth > maxdepth && lp->nconnection == maxconn) {
458 maxdepth = lp->depth;
459 root = lp;
460 }
461 }
462 }
463
464
determine_depths(void)465 static void determine_depths(void)
466 /*
467 * Determine the depth of all loops.
468 */
469
470 {
471 struct loop *lp;
472 int i,j;
473
474 for (i=0; i<loop_count; i++) {
475 lp = &loops[i];
476 for (j=0; j<loop_count; j++) loops[j].mark = false;
477 lp->depth = depth(lp);
478 }
479 }
480
481
482
depth(struct loop * lp)483 static int depth(struct loop *lp)
484 /*
485 * Determines the depth of loop, lp. Depth is defined as the minimum
486 * distance to a leaf loop where a leaf loop is one that has only one
487 * or no connections.
488 */
489
490 {
491 struct connection *cp,**cpp;
492 int count,ret,d;
493
494 if (lp->nconnection <= 1) return 0;
495 if (lp->mark) return -1;
496 lp->mark = true;
497 count = 0;
498 ret = 0;
499 for (cpp=lp->connections; (cp = *cpp); cpp++) {
500 d = depth(cp->loop);
501 if (d >= 0) {
502 if (++count == 1) ret = d;
503 else if (ret > d) ret = d;
504 }
505 }
506 lp->mark = false;
507 return ret+1;
508 }
509
510
traverse_loop(struct loop * lp,struct connection * anchor_connection)511 static void traverse_loop(struct loop *lp, struct connection *anchor_connection)
512 /*
513 * This is the workhorse of the display program. The algorithm is
514 * recursive based on processing individual loops. Each base pairing
515 * region is displayed using the direction given by the circle diagram,
516 * and the connections between the regions is drawn by equally spaced
517 * points. The radius of the loop is set to minimize the square error
518 * for lengths between sequential bases in the loops. The "correct"
519 * length for base links is 1. If the least squares fitting of the
520 * radius results in loops being less than 1/2 unit apart, then that
521 * segment is extruded.
522 *
523 * The variable, anchor_connection, gives the connection to the loop
524 * processed in an previous level of recursion.
525 */
526
527 {
528 double xs,ys,xe,ye,xn,yn,angleinc,r;
529 double radius,xc,yc,xo,yo,astart,aend,a;
530 struct connection *cp,*cpnext,**cpp,*acp,*cpprev;
531 int i,j,n,ic;
532 double da,maxang;
533 int count,icstart,icend,icmiddle,icroot;
534 logical done,done_all_connections,rooted;
535 int sign;
536 double midx,midy,nrx,nry,mx,my,vx,vy,dotmv,nmidx,nmidy;
537 int icstart1,icup,icdown,icnext,direction;
538 double dan,dx,dy,rr;
539 double cpx,cpy,cpnextx,cpnexty,cnx,cny,rcn,rc,lnx,lny,rl,ac,acn,sx,sy,dcp;
540 int imaxloop;
541
542 angleinc = 2 * pi / (nbase+1);
543 acp = NULL;
544 icroot = -1;
545 for (cpp=lp->connections, ic = 0; (cp = *cpp); cpp++, ic++) {
546 /* xs = cos(angleinc*cp->start);
547 ys = sin(angleinc*cp->start);
548 xe = cos(angleinc*cp->end);
549 ye = sin(angleinc*cp->end); */
550 xs = -sin(angleinc*cp->start);
551 ys = cos(angleinc*cp->start);
552 xe = -sin(angleinc*cp->end);
553 ye = cos(angleinc*cp->end);
554 xn = ye-ys;
555 yn = xs-xe;
556 r = sqrt(xn*xn + yn*yn);
557 cp->xrad = xn/r;
558 cp->yrad = yn/r;
559 cp->angle = atan2(yn,xn);
560 if (cp->angle < 0.0) cp->angle += 2*pi;
561 if (anchor_connection != NULL &&
562 anchor_connection->region == cp->region) {
563 acp = cp;
564 icroot = ic;
565 }
566 }
567
568 set_radius:
569 determine_radius(lp,lencut);
570 radius = lp->radius;
571 if (anchor_connection == NULL) xc = yc = 0.0;
572 else {
573 xo = (bases[acp->start].x+bases[acp->end].x) / 2.0;
574 yo = (bases[acp->start].y+bases[acp->end].y) / 2.0;
575 xc = xo - radius * acp->xrad;
576 yc = yo - radius * acp->yrad;
577 }
578
579 /*
580 * The construction of the connectors will proceed in blocks of
581 * connected connectors, where a connected connector pairs means
582 * two connectors that are forced out of the drawn circle because they
583 * are too close together in angle.
584 */
585
586 /*
587 * First, find the start of a block of connected connectors
588 */
589
590 if (icroot == -1)
591 icstart = 0;
592 else icstart = icroot;
593 cp = lp->connections[icstart];
594 count = 0;
595 if (debug) printf("Now processing loop %d\n",lp->number);
596 done = false;
597 do {
598 j = icstart - 1;
599 if (j < 0) j = lp->nconnection - 1;
600 cpprev = lp->connections[j];
601 if (!connected_connection(cpprev,cp)) {
602 done = true;
603 }
604 else {
605 icstart = j;
606 cp = cpprev;
607 }
608 if (++count > lp->nconnection) {
609 /*
610 * Here everything is connected. Break on maximum angular separation
611 * between connections.
612 */
613 maxang = -1.0;
614 for (ic = 0; ic < lp->nconnection; ic++) {
615 j = ic + 1;
616 if (j >= lp->nconnection) j = 0;
617 cp = lp->connections[ic];
618 cpnext = lp->connections[j];
619 ac = cpnext->angle - cp->angle;
620 if (ac < 0.0) ac += 2*pi;
621 if (ac > maxang) {
622 maxang = ac;
623 imaxloop = ic;
624 }
625 }
626 icend = imaxloop;
627 icstart = imaxloop + 1;
628 if (icstart >= lp->nconnection) icstart = 0;
629 cp = lp->connections[icend];
630 cp->broken = true;
631 done = true;
632 }
633 } while (!done);
634 done_all_connections = false;
635 icstart1 = icstart;
636 if (debug) printf("Icstart1 = %d\n",icstart1);
637 while (!done_all_connections) {
638 count = 0;
639 done = false;
640 icend = icstart;
641 rooted = false;
642 while (!done) {
643 cp = lp->connections[icend];
644 if (icend == icroot) rooted = true;
645 j = icend + 1;
646 if (j >= lp->nconnection) {
647 j = 0;
648 }
649 cpnext = lp->connections[j];
650 if (connected_connection(cp,cpnext)) {
651 if (++count >= lp->nconnection) break;
652 icend = j;
653 }
654 else {
655 done = true;
656 }
657 }
658 icmiddle = find_ic_middle(icstart,icend,anchor_connection,acp,lp);
659 ic = icup = icdown = icmiddle;
660 if (debug)
661 printf("IC start = %d middle = %d end = %d\n",
662 icstart,icmiddle,icend);
663 done = false;
664 direction = 0;
665 while (!done) {
666 if (direction < 0) {
667 ic = icup;
668 }
669 else if (direction == 0) {
670 ic = icmiddle;
671 }
672 else {
673 ic = icdown;
674 }
675 if (ic >= 0) {
676 cp = lp->connections[ic];
677 if (anchor_connection == NULL || acp != cp) {
678 if (direction == 0) {
679 astart = cp->angle - asin(1.0/2.0/radius);
680 aend = cp->angle + asin(1.0/2.0/radius);
681 bases[cp->start].x = xc + radius*cos(astart);
682 bases[cp->start].y = yc + radius*sin(astart);
683 bases[cp->end].x = xc + radius*cos(aend);
684 bases[cp->end].y = yc + radius*sin(aend);
685 }
686 else if (direction < 0) {
687 j = ic + 1;
688 if (j >= lp->nconnection) j = 0;
689 cp = lp->connections[ic];
690 cpnext = lp->connections[j];
691 cpx = cp->xrad;
692 cpy = cp->yrad;
693 ac = (cp->angle + cpnext->angle) / 2.0;
694 if (cp->angle > cpnext->angle) ac -= pi;
695 cnx = cos(ac);
696 cny = sin(ac);
697 lnx = cny;
698 lny = -cnx;
699 da = cpnext->angle - cp->angle;
700 if (da < 0.0) da += 2*pi;
701 if (cp->extruded) {
702 if (da <= pi/2) rl = 2.0;
703 else rl = 1.5;
704 }
705 else {
706 rl = 1.0;
707 }
708 bases[cp->end].x = bases[cpnext->start].x + rl*lnx;
709 bases[cp->end].y = bases[cpnext->start].y + rl*lny;
710 bases[cp->start].x = bases[cp->end].x + cpy;
711 bases[cp->start].y = bases[cp->end].y - cpx;
712 }
713 else {
714 j = ic - 1;
715 if (j < 0) j = lp->nconnection - 1;
716 cp = lp->connections[j];
717 cpnext = lp->connections[ic];
718 cpnextx = cpnext->xrad;
719 cpnexty = cpnext->yrad;
720 ac = (cp->angle + cpnext->angle) / 2.0;
721 if (cp->angle > cpnext->angle) ac -= pi;
722 cnx = cos(ac);
723 cny = sin(ac);
724 lnx = -cny;
725 lny = cnx;
726 da = cpnext->angle - cp->angle;
727 if (da < 0.0) da += 2*pi;
728 if (cp->extruded) {
729 if (da <= pi/2) rl = 2.0;
730 else rl = 1.5;
731 }
732 else {
733 rl = 1.0;
734 }
735 bases[cpnext->start].x = bases[cp->end].x + rl*lnx;
736 bases[cpnext->start].y = bases[cp->end].y + rl*lny;
737 bases[cpnext->end].x = bases[cpnext->start].x - cpnexty;
738 bases[cpnext->end].y = bases[cpnext->start].y + cpnextx;
739 }
740 }
741 }
742 if (direction < 0) {
743 if (icdown == icend) {
744 icdown = -1;
745 }
746 else if (icdown >= 0) {
747 if (++icdown >= lp->nconnection) {
748 icdown = 0;
749 }
750 }
751 direction = 1;
752 }
753 else {
754 if (icup == icstart) icup = -1;
755 else if (icup >= 0) {
756 if (--icup < 0) {
757 icup = lp->nconnection - 1;
758 }
759 }
760 direction = -1;
761 }
762 done = icup == -1 && icdown == -1;
763 }
764 icnext = icend + 1;
765 if (icnext >= lp->nconnection) icnext = 0;
766 if (icend != icstart && (! (icstart == icstart1 && icnext == icstart1))) {
767 /*
768 * Move the bases just constructed (or the radius) so
769 * that the bisector of the end points is radius distance
770 * away from the loop center.
771 */
772 cp = lp->connections[icstart];
773 cpnext = lp->connections[icend];
774 dx = bases[cpnext->end].x - bases[cp->start].x;
775 dy = bases[cpnext->end].y - bases[cp->start].y;
776 midx = bases[cp->start].x + dx/2.0;
777 midy = bases[cp->start].y + dy/2.0;
778 rr = sqrt(dx*dx + dy*dy);
779 mx = dx / rr;
780 my = dy / rr;
781 vx = xc - midx;
782 vy = yc - midy;
783 rr = sqrt(dx*dx + dy*dy);
784 vx /= rr;
785 vy /= rr;
786 dotmv = vx*mx + vy*my;
787 nrx = dotmv*mx - vx;
788 nry = dotmv*my - vy;
789 rr = sqrt(nrx*nrx + nry*nry);
790 nrx /= rr;
791 nry /= rr;
792 /*
793 * Determine which side of the bisector the center should be.
794 */
795 dx = bases[cp->start].x - xc;
796 dy = bases[cp->start].y - yc;
797 ac = atan2(dy,dx);
798 if (ac < 0.0) ac += 2*pi;
799 dx = bases[cpnext->end].x - xc;
800 dy = bases[cpnext->end].y - yc;
801 acn = atan2(dy,dx);
802 if (acn < 0.0) acn += 2*pi;
803 if (acn < ac) acn += 2*pi;
804 if (acn - ac > pi) sign = -1;
805 else sign = 1;
806 nmidx = xc + sign*radius*nrx;
807 nmidy = yc + sign*radius*nry;
808 if (rooted) {
809 xc -= nmidx - midx;
810 yc -= nmidy - midy;
811 }
812 else {
813 for (ic=icstart; ; ++ic >= lp->nconnection ? (ic = 0) : 0) {
814 cp = lp->connections[ic];
815 i = cp->start;
816 bases[i].x += nmidx - midx;
817 bases[i].y += nmidy - midy;
818 i = cp->end;
819 bases[i].x += nmidx - midx;
820 bases[i].y += nmidy - midy;
821 if (ic == icend) break;
822 }
823 }
824 }
825 icstart = icnext;
826 done_all_connections = icstart == icstart1;
827 }
828 for (ic=0; ic < lp->nconnection; ic++) {
829 cp = lp->connections[ic];
830 j = ic + 1;
831 if (j >= lp->nconnection) j = 0;
832 cpnext = lp->connections[j];
833 dx = bases[cp->end].x - xc;
834 dy = bases[cp->end].y - yc;
835 rc = sqrt(dx*dx + dy*dy);
836 ac = atan2(dy,dx);
837 if (ac < 0.0) ac += 2*pi;
838 dx = bases[cpnext->start].x - xc;
839 dy = bases[cpnext->start].y - yc;
840 rcn = sqrt(dx*dx + dy*dy);
841 acn = atan2(dy,dx);
842 if (acn < 0.0) acn += 2*pi;
843 if (acn < ac) acn += 2*pi;
844 dan = acn - ac;
845 dcp = cpnext->angle - cp->angle;
846 if (dcp <= 0.0) dcp += 2*pi;
847 if (fabs(dan-dcp) > pi) {
848 if (cp->extruded) {
849 vrna_message_warning("...from traverse_loop. Loop %d has crossed regions",
850 lp->number);
851 }
852 else if ((cpnext->start - cp->end) != 1) {
853 cp->extruded = true;
854 goto set_radius; /* Forever shamed */
855 }
856 }
857 if (cp->extruded) {
858 construct_extruded_segment(cp,cpnext);
859 }
860 else {
861 n = cpnext->start - cp->end;
862 if (n < 0) n += nbase + 1;
863 angleinc = dan / n;
864 for (j = 1; j < n; j++) {
865 i = cp->end + j;
866 if (i > nbase) i -= nbase + 1;
867 a = ac + j*angleinc;
868 rr = rc + (rcn-rc)*(a-ac)/dan;
869 bases[i].x = xc + rr*cos(a);
870 bases[i].y = yc + rr*sin(a);
871 }
872 }
873 }
874 for (ic=0; ic < lp->nconnection; ic++) {
875 if (icroot != ic) {
876 cp = lp->connections[ic];
877 generate_region(cp);
878 traverse_loop(cp->loop,cp);
879 }
880 }
881 n = 0;
882 sx = 0.0;
883 sy = 0.0;
884 for (ic = 0; ic < lp->nconnection; ic++) {
885 j = ic + 1;
886 if (j >= lp->nconnection) j = 0;
887 cp = lp->connections[ic];
888 cpnext = lp->connections[j];
889 n += 2;
890 sx += bases[cp->start].x + bases[cp->end].x;
891 sy += bases[cp->start].y + bases[cp->end].y;
892 if (!cp->extruded) {
893 for (j = cp->end + 1; j != cpnext->start; j++) {
894 if (j > nbase) j -= nbase + 1;
895 n++;
896 sx += bases[j].x;
897 sy += bases[j].y;
898 }
899 }
900 }
901 lp->x = sx / n;
902 lp->y = sy / n;
903
904 /* free connections (ih) */
905 for (ic = 0; ic < lp->nconnection; ic++)
906 free(lp->connections[ic]);
907 free(lp->connections);
908 }
909
910
911
determine_radius(struct loop * lp,double lencut)912 static void determine_radius(struct loop *lp,double lencut)
913 /*
914 * For the loop pointed to by lp, determine the radius of
915 * the loop that will ensure that each base around the loop will
916 * have a separation of at least lencut around the circle.
917 * If a segment joining two connectors will not support this separation,
918 * then the flag, extruded, will be set in the first of these
919 * two indicators. The radius is set in lp.
920 *
921 * The radius is selected by a least squares procedure where the sum of the
922 * squares of the deviations of length from the ideal value of 1 is used
923 * as the error function.
924 */
925
926 {
927 double mindit,ci,dt,sumn,sumd,radius,dit;
928 int i,j,end,start,imindit;
929 struct connection *cp,*cpnext;
930 static double rt2_2 = 0.7071068;
931
932 do {
933 mindit = 1.0e10;
934 for (sumd=0.0, sumn=0.0, i=0;
935 i < lp->nconnection;
936 i++) {
937 cp = lp->connections[i];
938 j = i + 1;
939 if (j >= lp->nconnection) j = 0;
940 cpnext = lp->connections[j];
941 end = cp->end;
942 start = cpnext->start;
943 if (start < end) start += nbase + 1;
944 dt = cpnext->angle - cp->angle;
945 if (dt <= 0.0) dt += 2*pi;
946 if (!cp->extruded)
947 ci = start - end;
948 else {
949 if (dt <= pi/2) ci = 2.0;
950 else ci = 1.5;
951 }
952 sumn += dt * (1.0/ci + 1.0);
953 sumd += dt * dt / ci;
954 dit = dt/ci;
955 if (dit < mindit && !cp->extruded && ci > 1.0) {
956 mindit = dit;
957 imindit = i;
958 }
959 }
960 radius = sumn/sumd;
961 if (radius < rt2_2) radius = rt2_2;
962 if (mindit*radius < lencut) {
963 lp->connections[imindit]->extruded = true;
964 }
965 } while (mindit*radius < lencut);
966 if (lp->radius > 0.0)
967 radius = lp->radius;
968 else lp->radius = radius;
969 }
970
971
connected_connection(struct connection * cp,struct connection * cpnext)972 static logical connected_connection(struct connection *cp, struct connection *cpnext)
973 /*
974 * Determines if the connections cp and cpnext are connected
975 */
976
977 {
978
979 if (cp->extruded) {
980 return true;
981 }
982 else if (cp->end+1 == cpnext->start) {
983 return true;
984 }
985 else {
986 return false;
987 }
988 }
989
990
find_ic_middle(int icstart,int icend,struct connection * anchor_connection,struct connection * acp,struct loop * lp)991 static int find_ic_middle(int icstart, int icend, struct connection *anchor_connection, struct connection *acp, struct loop *lp)
992 /*
993 * Finds the middle of a set of connected connectors. This is normally
994 * the middle connection in the sequence except if one of the connections
995 * is the anchor, in which case that connection will be used.
996 */
997
998 {
999 int count,ret,ic,i;
1000 logical done;
1001
1002 count = 0;
1003 ret = -1;
1004 ic = icstart;
1005 done = false;
1006 while (!done) {
1007 if (count++ > lp->nconnection * 2) {
1008 printf("Infinite loop detected in find_ic_middle\n");
1009 exit(FATAL_ERROR);
1010 }
1011 if (anchor_connection != NULL && lp->connections[ic] == acp) {
1012 ret = ic;
1013 }
1014 done = ic == icend;
1015 if (++ic >= lp->nconnection) {
1016 ic = 0;
1017 }
1018 }
1019 if (ret == -1) {
1020 for (i=1, ic=icstart; i<(count+1)/2; i++) {
1021 if (++ic >= lp->nconnection) ic = 0;
1022 }
1023 ret = ic;
1024 }
1025 return ret;
1026 }
1027
1028
generate_region(struct connection * cp)1029 static void generate_region(struct connection *cp)
1030 /*
1031 * Generates the coordinates for the base pairing region of a connection
1032 * given the position of the starting base pair.
1033 */
1034
1035 {
1036 int l,start,end,i,mate;
1037 struct region *rp;
1038
1039 rp = cp->region;
1040 l = 0;
1041 if (cp->start == rp->start1) {
1042 start = rp->start1;
1043 end = rp->end1;
1044 }
1045 else {
1046 start = rp->start2;
1047 end = rp->end2;
1048 }
1049 if (bases[cp->start].x > anum - 100.0 ||
1050 bases[cp->end].x > anum - 100.0) {
1051 printf("Bad region passed to generate_region. Coordinates not defined.\n");
1052 exit(FATAL_ERROR);
1053 }
1054 for (i=start+1; i<=end; i++) {
1055 l++;
1056 bases[i].x = bases[cp->start].x + l*cp->xrad;
1057 bases[i].y = bases[cp->start].y + l*cp->yrad;
1058 mate = bases[i].mate;
1059 bases[mate].x = bases[cp->end].x + l*cp->xrad;
1060 bases[mate].y = bases[cp->end].y + l*cp->yrad;
1061 }
1062 }
1063
1064
construct_circle_segment(int start,int end)1065 static void construct_circle_segment(int start, int end)
1066 /*
1067 * Draws the segment of residue between the bases numbered start
1068 * through end, where start and end are presumed to be part of a base
1069 * pairing region. They are drawn as a circle which has a chord given
1070 * by the ends of two base pairing regions defined by the connections.
1071 */
1072
1073 {
1074 double dx,dy,rr,h,angleinc,midx,midy,xn,yn,nrx,nry,mx,my,a;
1075 int l,j,i;
1076
1077 dx = bases[end].x - bases[start].x;
1078 dy = bases[end].y - bases[start].y;
1079 rr = sqrt(dx*dx + dy*dy);
1080 l = end - start;
1081 if (l < 0) l += nbase + 1;
1082 if (rr >= l) {
1083 dx /= rr;
1084 dy /= rr;
1085 for (j = 1; j < l; j++) {
1086 i = start + j;
1087 if (i > nbase) i -= nbase + 1;
1088 bases[i].x = bases[start].x + dx*(double)j/(double)l;
1089 bases[i].y = bases[start].y + dy*(double)j/(double)l;
1090 }
1091 }
1092 else {
1093 find_center_for_arc(l-1,rr,&h,&angleinc);
1094 dx /= rr;
1095 dy /= rr;
1096 midx = bases[start].x + dx*rr/2.0;
1097 midy = bases[start].y + dy*rr/2.0;
1098 xn = dy;
1099 yn = -dx;
1100 nrx = midx + h*xn;
1101 nry = midy + h*yn;
1102 mx = bases[start].x - nrx;
1103 my = bases[start].y - nry;
1104 rr = sqrt(mx*mx + my*my);
1105 a = atan2(my,mx);
1106 for (j = 1; j < l; j++) {
1107 i = start + j;
1108 if (i > nbase) i -= nbase + 1;
1109 bases[i].x = nrx + rr*cos(a+j*angleinc);
1110 bases[i].y = nry + rr*sin(a+j*angleinc);
1111 }
1112 }
1113 }
1114
1115
construct_extruded_segment(struct connection * cp,struct connection * cpnext)1116 static void construct_extruded_segment(struct connection *cp, struct connection *cpnext)
1117 /*
1118 * Constructs the segment between cp and cpnext as a circle if possible.
1119 * However, if the segment is too large, the lines are drawn between
1120 * the two connecting regions, and bases are placed there until the
1121 * connecting circle will fit.
1122 */
1123
1124 {
1125 double astart,aend1,aend2,aave,dx,dy,a1,a2,ac,rr,da,dac;
1126 int start,end,n,nstart,nend;
1127 logical collision;
1128
1129 astart = cp->angle;
1130 aend2 = aend1 = cpnext->angle;
1131 if (aend2 < astart) aend2 += 2*pi;
1132 aave = (astart + aend2) / 2.0;
1133 start = cp->end;
1134 end = cpnext->start;
1135 n = end - start;
1136 if (n < 0) n += nbase + 1;
1137 da = cpnext->angle - cp->angle;
1138 if (da < 0.0) {
1139 da += 2*pi;
1140 }
1141 if (n == 2) construct_circle_segment(start,end);
1142 else {
1143 dx = bases[end].x - bases[start].x;
1144 dy = bases[end].y - bases[start].y;
1145 rr = sqrt(dx*dx + dy*dy);
1146 dx /= rr;
1147 dy /= rr;
1148 if (rr >= 1.5 && da <= pi/2) {
1149 nstart = start + 1;
1150 if (nstart > nbase) nstart -= nbase + 1;
1151 nend = end - 1;
1152 if (nend < 0) nend += nbase + 1;
1153 bases[nstart].x = bases[start].x + 0.5*dx;
1154 bases[nstart].y = bases[start].y + 0.5*dy;
1155 bases[nend].x = bases[end].x - 0.5*dx;
1156 bases[nend].y = bases[end].y - 0.5*dy;
1157 start = nstart;
1158 end = nend;
1159 }
1160 do {
1161 collision = false;
1162 construct_circle_segment(start,end);
1163 nstart = start + 1;
1164 if (nstart > nbase) nstart -= nbase + 1;
1165 dx = bases[nstart].x - bases[start].x;
1166 dy = bases[nstart].y - bases[start].y;
1167 a1 = atan2(dy,dx);
1168 if (a1 < 0.0) a1 += 2*pi;
1169 dac = a1 - astart;
1170 if (dac < 0.0) dac += 2*pi;
1171 if (dac > pi) collision = true;
1172 nend = end - 1;
1173 if (nend < 0) nend += nbase + 1;
1174 dx = bases[nend].x - bases[end].x;
1175 dy = bases[nend].y - bases[end].y;
1176 a2 = atan2(dy,dx);
1177 if (a2 < 0.0) a2 += 2*pi;
1178 dac = aend1 - a2;
1179 if (dac < 0.0) dac += 2*pi;
1180 if (dac > pi) collision = true;
1181 if (collision) {
1182 ac = minf2(aave,astart + 0.5);
1183 bases[nstart].x = bases[start].x + cos(ac);
1184 bases[nstart].y = bases[start].y + sin(ac);
1185 start = nstart;
1186 ac = maxf2(aave,aend2 - 0.5);
1187 bases[nend].x = bases[end].x + cos(ac);
1188 bases[nend].y = bases[end].y + sin(ac);
1189 end = nend;
1190 n -= 2;
1191 }
1192 } while (collision && n > 1);
1193 }
1194 }
1195
1196
find_center_for_arc(int n,double b,double * hp,double * thetap)1197 static void find_center_for_arc(int n,double b,double *hp,double *thetap)
1198 /*
1199 * Given n points to be placed equidistantly and equiangularly on a
1200 * polygon which has a chord of length, b, find the distance, h, from the
1201 * midpoint of the chord for the center of polygon. Positive values
1202 * mean the center is within the polygon and the chord, whereas
1203 * negative values mean the center is outside the chord. Also, the
1204 * radial angle for each polygon side is returned in theta.
1205 *
1206 * The procedure uses a bisection algorithm to find the correct
1207 * value for the center. Two equations are solved, the angles
1208 * around the center must add to 2*pi, and the sides of the polygon
1209 * excluding the chord must have a length of 1.
1210 */
1211
1212 {
1213 double h,hhi,hlow,r,disc,theta,e,phi;
1214 int iter;
1215 #define maxiter 500
1216
1217 hhi = (n+1) / pi;
1218 hlow = - hhi - b/(n+1.000001-b); /* changed to prevent div by zero if (ih) */
1219 if (b<1) hlow = 0; /* otherwise we might fail below (ih) */
1220 iter = 0;
1221 do {
1222 h = (hhi + hlow) / 2.0;
1223 r = sqrt(h*h + b*b/4.0);
1224 /* if (r<0.5) {r = 0.5; h = 0.5*sqrt(1-b*b);} */
1225 disc = 1.0 - 0.5/(r*r);
1226 if (fabs(disc) > 1.0) {
1227 vrna_message_error("Unexpected large magnitude discriminant = %g %g", disc,r);
1228 exit(FATAL_ERROR);
1229 }
1230 theta = acos(disc);
1231 /* theta = 2*acos(sqrt(1-1/(4*r*r))); */
1232 phi = acos(h/r);
1233 e = theta * (n+1) + 2*phi - 2*pi;
1234 if (e > 0.0) {
1235 hlow = h;
1236 }
1237 else {
1238 hhi = h;
1239 }
1240 } while (fabs(e) > 0.0001 && ++iter < maxiter);
1241 if (iter >= maxiter) {
1242 vrna_message_warning("Iteration failed in find_center_for_arc");
1243 h = 0.0;
1244 theta = 0.0;
1245 }
1246 *hp = h;
1247 *thetap = theta;
1248 }
1249
1250