binaries/src/ViennaRNA/lib/naview.c

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