2 VARNA is a tool for the automated drawing, visualization and annotation of the secondary structure of RNA, designed as a companion software for web servers and databases.
3 Copyright (C) 2008 Kevin Darty, Alain Denise and Yann Ponty.
4 electronic mail : Yann.Ponty@lri.fr
5 paper mail : LRI, bat 490 Université Paris-Sud 91405 Orsay Cedex France
7 This file is part of VARNA version 3.1.
8 VARNA version 3.1 is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License
9 as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
11 VARNA version 3.1 is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
12 without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 See the GNU General Public License for more details.
15 You should have received a copy of the GNU General Public License along with VARNA version 3.1.
16 If not, see http://www.gnu.org/licenses.
18 package fr.orsay.lri.varna.models.naView;
21 import java.util.ArrayList;
23 import fr.orsay.lri.varna.exceptions.ExceptionNAViewAlgorithm;
24 import fr.orsay.lri.varna.interfaces.InterfaceVARNAListener;
25 import fr.orsay.lri.varna.interfaces.InterfaceVARNAObservable;
28 private final double ANUM = 9999.0;
29 private final int MAXITER = 500;
31 private ArrayList<Base> bases;
32 private int nbase, nregion, loop_count;
34 private Loop root = new Loop();
35 private ArrayList<Loop> loops;
37 private ArrayList<Region> regions;
39 private Radloop rlphead = new Radloop();
41 private double lencut=0.8;
42 private final double RADIUS_REDUCTION_FACTOR = 1.4;
44 // show algorithm step by step
45 private boolean debug = false;
47 private double angleinc;
51 private ArrayList<InterfaceVARNAListener> _listeVARNAListener = new ArrayList<InterfaceVARNAListener>();
53 private boolean noIterationFailureYet = true;
55 double HELIX_FACTOR = 0.6;
56 double BACKBONE_DISTANCE = 27;
58 public int naview_xy_coordinates(ArrayList<Short> pair_table2,
59 ArrayList<Double> x, ArrayList<Double> y)
60 throws ExceptionNAViewAlgorithm {
62 System.out.println("naview_xy_coordinates");
63 if (pair_table2.size() == 0)
66 ArrayList<Integer> pair_table = new ArrayList<Integer>(pair_table2
68 pair_table.add(pair_table2.size());
70 for (int j = 0; j < pair_table2.size(); j++) {
71 pair_table.add(pair_table2.get(j) + 1);
75 infoStructure(pair_table);
78 nbase = pair_table.get(0);
79 bases = new ArrayList<Base>(nbase + 1);
81 for (int index = 0; index < bases.size(); index++) {
82 bases.add(new Base());
85 regions = new ArrayList<Region>();
86 for (int index = 0; index < nbase + 1; index++) {
87 regions.add(new Region());
90 read_in_bases(pair_table);
103 loops = new ArrayList<Loop>(nbase + 1);
104 for (int index = 0; index < nbase + 1; index++) {
105 loops.add(new Loop());
111 infoBasesExtracted();
121 traverse_loop(root, null);
123 for (i = 0; i < nbase; i++) {
124 x.add(100 + BACKBONE_DISTANCE * bases.get(i + 1).getX());
125 y.add(100 + BACKBONE_DISTANCE * bases.get(i + 1).getY());
131 private void infoStructure(ArrayList<Integer> pair_table) {
132 System.out.println("structure:");
133 for (int j = 0; j < pair_table.size(); j++) {
134 System.out.print("#" + j + ":" + pair_table.get(j) + "\t");
136 System.out.println();
138 System.out.println();
141 private void infoBasesMate() {
142 System.out.println("Bases mate:");
143 for (int index = 0; index < bases.size(); index++) {
144 System.out.print("#" + index + ":" + bases.get(index).getMate()
147 System.out.println();
149 System.out.println();
152 private void infoRegions() {
153 System.out.println("regions:");
154 for (int index = 0; index < regions.size(); index++) {
155 System.out.print("(" + regions.get(index).getStart1() + ","
156 + regions.get(index).getStart2() + ";"
157 + regions.get(index).getEnd1() + ","
158 + regions.get(index).getEnd2() + ")\t\t");
160 System.out.println();
162 System.out.println();
165 private void infoBasesExtracted() {
166 System.out.println("Bases extracted:");
167 for (int index = 0; index < bases.size(); index++) {
168 System.out.print("i=" + index + ":"
169 + bases.get(index).isExtracted() + "\t");
171 System.out.println();
173 System.out.println();
176 private void infoRoot() {
177 System.out.println("root" + root.getNconnection() + ";"
179 System.out.println("\troot : ");
180 System.out.println("\tdepth=" + root.getDepth());
181 System.out.println("\tmark=" + root.isMark());
182 System.out.println("\tnumber=" + root.getNumber());
183 System.out.println("\tradius=" + root.getRadius());
184 System.out.println("\tx=" + root.getX());
185 System.out.println("\ty=" + root.getY());
186 System.out.println("\tnconnection=" + root.getNconnection());
189 private void read_in_bases(ArrayList<Integer> pair_table) {
191 System.out.println("read_in_bases");
196 bases.add(new Base());
197 bases.get(0).setMate(0);
198 bases.get(0).setExtracted(false);
199 bases.get(0).setX(ANUM);
200 bases.get(0).setY(ANUM);
202 for (npairs = 0, i = 1; i <= nbase; i++) {
203 bases.add(new Base());
204 bases.get(i).setExtracted(false);
205 bases.get(i).setX(ANUM);
206 bases.get(i).setY(ANUM);
207 bases.get(i).setMate(pair_table.get(i));
208 if ((int) pair_table.get(i) > i)
211 // must have at least 1 pair to avoid segfault
213 bases.get(1).setMate(nbase);
214 bases.get(nbase).setMate(1);
219 * Identifies the regions in the structure.
221 private void find_regions()
225 System.out.println("find_regions");
228 ArrayList<Boolean> mark = new ArrayList<Boolean>(nb1);
229 for (i = 0; i < nb1; i++)
232 for (i = 0; i <= nbase; i++) {
233 if ((mate = bases.get(i).getMate()) != 0 && !mark.get(i)) {
234 regions.get(nregion).setStart1(i);
235 regions.get(nregion).setEnd2(mate);
237 mark.set(mate, true);
238 bases.get(i).setRegion(regions.get(nregion));
239 bases.get(mate).setRegion(regions.get(nregion));
240 for (i++, mate--; i < mate && bases.get(i).getMate() == mate; i++, mate--) {
241 mark.set(mate, true);
243 bases.get(i).setRegion(regions.get(nregion));
244 bases.get(mate).setRegion(regions.get(nregion));
246 regions.get(nregion).setEnd1(--i);
247 regions.get(nregion).setStart2(mate + 1);
250 System.out.printf("\nRegions are:\n");
252 "Region %d is %d-%d and %d-%d with gap of %d.\n",
253 nregion + 1, regions.get(nregion).getStart1(),
254 regions.get(nregion).getEnd1(), regions
255 .get(nregion).getStart2(), regions.get(
256 nregion).getEnd2(), regions.get(nregion)
258 - regions.get(nregion).getEnd1() + 1);
266 * Starting at residue ibase, recursively constructs the loop containing
267 * said base and all deeper bases.
269 * @throws ExceptionNAViewAlgorithm
271 private Loop construct_loop(int ibase) throws ExceptionNAViewAlgorithm {
273 System.out.println("construct_loop");
275 Loop retloop = new Loop(), lp = new Loop();
276 Connection cp = new Connection();
277 Region rp = new Region();
278 Radloop rlp = new Radloop();
279 retloop = loops.get(loop_count++);
280 retloop.setNconnection(0);
281 //System.out.println(""+ibase+" "+nbase);
282 //ArrayList<Connection> a = new ArrayList<Connection>(nbase + 1);
283 //retloop.setConnections(a);
284 // for (int index = 0; index < nbase + 1; index++)
285 // retloop.getConnections().add(new Connection());
286 //for (int index = 0; index < nbase + 1; index++)
287 // retloop.addConnection(index,new Connection());
289 retloop.setNumber(loop_count);
290 retloop.setRadius(0.0);
291 for (rlp = rlphead; rlp != null; rlp = rlp.getNext())
292 if (rlp.getLoopnumber() == loop_count)
293 retloop.setRadius(rlp.getRadius());
296 if ((mate = bases.get(i).getMate()) != 0) {
297 rp = bases.get(i).getRegion();
298 if (!bases.get(rp.getStart1()).isExtracted()) {
299 if (i == rp.getStart1()) {
300 bases.get(rp.getStart1()).setExtracted(true);
301 bases.get(rp.getEnd1()).setExtracted(true);
302 bases.get(rp.getStart2()).setExtracted(true);
303 bases.get(rp.getEnd2()).setExtracted(true);
304 lp = construct_loop(rp.getEnd1() < nbase ? rp.getEnd1() + 1
306 } else if (i == rp.getStart2()) {
307 bases.get(rp.getStart2()).setExtracted(true);
308 bases.get(rp.getEnd2()).setExtracted(true);
309 bases.get(rp.getStart1()).setExtracted(true);
310 bases.get(rp.getEnd1()).setExtracted(true);
311 lp = construct_loop(rp.getEnd2() < nbase ? rp.getEnd2() + 1
314 throw new ExceptionNAViewAlgorithm(
315 "naview:Error detected in construct_loop. i = "
316 + i + " not found in region table.\n");
318 retloop.setNconnection(retloop.getNconnection() + 1);
319 cp = new Connection();
320 retloop.setConnection(retloop.getNconnection() - 1, cp);
321 retloop.setConnection(retloop.getNconnection(), null);
324 if (i == rp.getStart1()) {
325 cp.setStart(rp.getStart1());
326 cp.setEnd(rp.getEnd2());
328 cp.setStart(rp.getStart2());
329 cp.setEnd(rp.getEnd1());
331 cp.setExtruded(false);
333 lp.setNconnection(lp.getNconnection() + 1);
334 cp = new Connection();
335 lp.setConnection(lp.getNconnection() - 1, cp);
336 lp.setConnection(lp.getNconnection(), null);
339 if (i == rp.getStart1()) {
340 cp.setStart(rp.getStart2());
341 cp.setEnd(rp.getEnd1());
343 cp.setStart(rp.getStart1());
344 cp.setEnd(rp.getEnd2());
346 cp.setExtruded(false);
353 } while (i != ibase);
358 * Displays all the loops.
360 private void dump_loops() {
361 System.out.println("dump_loops");
366 System.out.printf("\nRoot loop is #%d\n", loops.indexOf(root) + 1);
367 for (il = 0; il < loop_count; il++) {
369 System.out.printf("Loop %d has %d connections:\n", il + 1, lp
371 for (int i = 0; (cp = lp.getConnection(i)) != null; i++) {
372 ilp = (loops.indexOf(cp.getLoop())) + 1;
373 irp = (regions.indexOf(cp.getRegion())) + 1;
374 System.out.printf(" Loop %d Region %d (%d-%d)\n", ilp, irp, cp
375 .getStart(), cp.getEnd());
381 * Find node of greatest branching that is deepest.
383 private void find_central_loop() {
385 System.out.println("find_central_loop");
386 Loop lp = new Loop();
387 int maxconn, maxdepth, i;
392 for (i = 0; i < loop_count; i++) {
394 if (lp.getNconnection() > maxconn) {
395 maxdepth = lp.getDepth();
396 maxconn = lp.getNconnection();
398 } else if (lp.getDepth() > maxdepth
399 && lp.getNconnection() == maxconn) {
400 maxdepth = lp.getDepth();
407 * Determine the depth of all loops.
409 private void determine_depths() {
411 System.out.println("determine_depths");
412 Loop lp = new Loop();
415 for (i = 0; i < loop_count; i++) {
417 for (j = 0; j < loop_count; j++)
418 loops.get(j).setMark(false);
419 lp.setDepth(depth(lp));
424 * Determines the depth of loop, lp. Depth is defined as the minimum
425 * distance to a leaf loop where a leaf loop is one that has only one or no
428 private int depth(Loop lp) {
430 System.out.println("depth");
433 if (lp.getNconnection() <= 1)
440 for (int i = 0; lp.getConnection(i) != null; i++) {
441 d = depth(lp.getConnection(i).getLoop());
454 * This is the workhorse of the display program. The algorithm is recursive
455 * based on processing individual loops. Each base pairing region is
456 * displayed using the direction given by the circle diagram, and the
457 * connections between the regions is drawn by equally spaced points. The
458 * radius of the loop is set to minimize the square error for lengths
459 * between sequential bases in the loops. The "correct" length for base
460 * links is 1. If the least squares fitting of the radius results in loops
461 * being less than 1/2 unit apart, then that segment is extruded.
463 * The variable, anchor_connection, gives the connection to the loop
464 * processed in an previous level of recursion.
466 * @throws ExceptionNAViewAlgorithm
468 private void traverse_loop(Loop lp, Connection anchor_connection)
469 throws ExceptionNAViewAlgorithm {
471 System.out.println(" traverse_loop");
472 double xs, ys, xe, ye, xn, yn, angleinc, r;
473 double radius, xc, yc, xo, yo, astart, aend, a;
474 Connection cp, cpnext, acp, cpprev;
477 int count, icstart, icend, icmiddle, icroot;
478 boolean done, done_all_connections, rooted;
480 double midx, midy, nrx, nry, mx, my, vx, vy, dotmv, nmidx, nmidy;
481 int icstart1, icup, icdown, icnext, direction;
482 double dan, dx, dy, rr;
483 double cpx, cpy, cpnextx, cpnexty, cnx, cny, rcn, rc, lnx, lny, rl, ac, acn, sx, sy, dcp;
486 angleinc = 2 * Math.PI / (nbase + 1);
491 for (ic = 0; (cp = lp.getConnection(indice)) != null; indice++, ic++) {
492 // xs = cos(angleinc*cp.setStart(); ys = sin(angleinc*cp.setStart();
494 // cos(angleinc*cp.setEnd()); ye = sin(angleinc*cp.setEnd());
495 xs = -Math.sin(angleinc * cp.getStart());
496 ys = Math.cos(angleinc * cp.getStart());
497 xe = -Math.sin(angleinc * cp.getEnd());
498 ye = Math.cos(angleinc * cp.getEnd());
501 r = Math.sqrt(xn * xn + yn * yn);
504 cp.setAngle(Math.atan2(yn, xn));
505 if (cp.getAngle() < 0.0)
506 cp.setAngle(cp.getAngle() + 2 * Math.PI);
507 if (anchor_connection != null
508 && anchor_connection.getRegion() == cp.getRegion()) {
513 // remplacement d'une etiquette de goto
514 set_radius: while (true) {
515 determine_radius(lp, lencut);
516 radius = lp.getRadius()/RADIUS_REDUCTION_FACTOR;
517 if (anchor_connection == null)
520 xo = (bases.get(acp.getStart()).getX() + bases
521 .get(acp.getEnd()).getX()) / 2.0;
522 yo = (bases.get(acp.getStart()).getY() + bases
523 .get(acp.getEnd()).getY()) / 2.0;
524 xc = xo - radius * acp.getXrad();
525 yc = yo - radius * acp.getYrad();
528 // The construction of the connectors will proceed in blocks of
529 // connected connectors, where a connected connector pairs means two
530 // connectors that are forced out of the drawn circle because they
531 // are too close together in angle.
533 // First, find the start of a block of connected connectors
539 cp = lp.getConnection(icstart);
543 System.out.printf("Now processing loop %d\n", lp.getNumber());
544 System.out.println(" "+lp);
550 j = lp.getNconnection() - 1;
551 cpprev = lp.getConnection(j);
552 if (!connected_connection(cpprev, cp)) {
558 if (++count > lp.getNconnection()) {
559 // Here everything is connected. Break on maximum angular
560 // separation between connections.
563 for (ic = 0; ic < lp.getNconnection(); ic++) {
565 if (j >= lp.getNconnection())
567 cp = lp.getConnection(ic);
568 cpnext = lp.getConnection(j);
569 ac = cpnext.getAngle() - cp.getAngle();
578 icstart = imaxloop + 1;
579 if (icstart >= lp.getNconnection())
581 cp = lp.getConnection(icend);
586 done_all_connections = false;
589 System.out.printf(" Icstart1 = %d\n", icstart1);
590 while (!done_all_connections) {
596 cp = lp.getConnection(icend);
600 if (j >= lp.getNconnection()) {
603 cpnext = lp.getConnection(j);
604 if (connected_connection(cp, cpnext)) {
605 if (++count >= lp.getNconnection())
612 icmiddle = find_ic_middle(icstart, icend, anchor_connection,
614 ic = icup = icdown = icmiddle;
616 System.out.printf(" IC start = %d middle = %d end = %d\n",
617 icstart, icmiddle, icend);
623 } else if (direction == 0) {
629 cp = lp.getConnection(ic);
630 if (anchor_connection == null || acp != cp) {
631 if (direction == 0) {
632 astart = cp.getAngle()
633 - Math.asin(1.0 / 2.0 / radius);
635 + Math.asin(1.0 / 2.0 / radius);
636 bases.get(cp.getStart()).setX(
637 xc + radius * Math.cos(astart));
638 bases.get(cp.getStart()).setY(
639 yc + radius * Math.sin(astart));
640 bases.get(cp.getEnd()).setX(
641 xc + radius * Math.cos(aend));
642 bases.get(cp.getEnd()).setY(
643 yc + radius * Math.sin(aend));
644 } else if (direction < 0) {
646 if (j >= lp.getNconnection())
648 cp = lp.getConnection(ic);
649 cpnext = lp.getConnection(j);
652 ac = (cp.getAngle() + cpnext.getAngle()) / 2.0;
653 if (cp.getAngle() > cpnext.getAngle())
659 da = cpnext.getAngle() - cp.getAngle();
662 if (cp.isExtruded()) {
663 if (da <= Math.PI / 2)
670 bases.get(cp.getEnd()).setX(
671 bases.get(cpnext.getStart()).getX()
673 bases.get(cp.getEnd()).setY(
674 bases.get(cpnext.getStart()).getY()
676 bases.get(cp.getStart()).setX(
677 bases.get(cp.getEnd()).getX() + cpy);
678 bases.get(cp.getStart()).setY(
679 bases.get(cp.getEnd()).getY() - cpx);
683 j = lp.getNconnection() - 1;
684 cp = lp.getConnection(j);
685 cpnext = lp.getConnection(ic);
686 cpnextx = cpnext.getXrad();
687 cpnexty = cpnext.getYrad();
688 ac = (cp.getAngle() + cpnext.getAngle()) / 2.0;
689 if (cp.getAngle() > cpnext.getAngle())
695 da = cpnext.getAngle() - cp.getAngle();
698 if (cp.isExtruded()) {
699 if (da <= Math.PI / 2)
706 bases.get(cpnext.getStart()).setX(
707 bases.get(cp.getEnd()).getX() + rl
709 bases.get(cpnext.getStart()).setY(
710 bases.get(cp.getEnd()).getY() + rl
712 bases.get(cpnext.getEnd()).setX(
713 bases.get(cpnext.getStart()).getX()
715 bases.get(cpnext.getEnd()).setY(
716 bases.get(cpnext.getStart()).getY()
722 if (icdown == icend) {
724 } else if (icdown >= 0) {
725 if (++icdown >= lp.getNconnection()) {
733 else if (icup >= 0) {
735 icup = lp.getNconnection() - 1;
740 done = icup == -1 && icdown == -1;
743 if (icnext >= lp.getNconnection())
746 && (!(icstart == icstart1 && icnext == icstart1))) {
748 // Move the bases just constructed (or the radius) so that
749 // the bisector of the end points is radius distance away
750 // from the loop center.
752 cp = lp.getConnection(icstart);
753 cpnext = lp.getConnection(icend);
754 dx = bases.get(cpnext.getEnd()).getX()
755 - bases.get(cp.getStart()).getX();
756 dy = bases.get(cpnext.getEnd()).getY()
757 - bases.get(cp.getStart()).getY();
758 midx = bases.get(cp.getStart()).getX() + dx / 2.0;
759 midy = bases.get(cp.getStart()).getY() + dy / 2.0;
760 rr = Math.sqrt(dx * dx + dy * dy);
765 rr = Math.sqrt(dx * dx + dy * dy);
768 dotmv = vx * mx + vy * my;
769 nrx = dotmv * mx - vx;
770 nry = dotmv * my - vy;
771 rr = Math.sqrt(nrx * nrx + nry * nry);
775 // Determine which side of the bisector the center should
778 dx = bases.get(cp.getStart()).getX() - xc;
779 dy = bases.get(cp.getStart()).getY() - yc;
780 ac = Math.atan2(dy, dx);
783 dx = bases.get(cpnext.getEnd()).getX() - xc;
784 dy = bases.get(cpnext.getEnd()).getY() - yc;
785 acn = Math.atan2(dy, dx);
790 if (acn - ac > Math.PI)
794 nmidx = xc + sign * radius * nrx;
795 nmidy = yc + sign * radius * nry;
800 for (ic = icstart;;) {
801 cp = lp.getConnection(ic);
804 bases.get(i).getX() + nmidx - midx);
806 bases.get(i).getY() + nmidy - midy);
809 bases.get(i).getX() + nmidx - midx);
811 bases.get(i).getY() + nmidy - midy);
814 if (++ic >= lp.getNconnection())
820 done_all_connections = icstart == icstart1;
822 for (ic = 0; ic < lp.getNconnection(); ic++) {
823 cp = lp.getConnection(ic);
825 if (j >= lp.getNconnection())
827 cpnext = lp.getConnection(j);
828 dx = bases.get(cp.getEnd()).getX() - xc;
829 dy = bases.get(cp.getEnd()).getY() - yc;
830 rc = Math.sqrt(dx * dx + dy * dy);
831 ac = Math.atan2(dy, dx);
834 dx = bases.get(cpnext.getStart()).getX() - xc;
835 dy = bases.get(cpnext.getStart()).getY() - yc;
836 rcn = Math.sqrt(dx * dx + dy * dy);
837 acn = Math.atan2(dy, dx);
843 dcp = cpnext.getAngle() - cp.getAngle();
846 if (Math.abs(dan - dcp) > Math.PI) {
847 if (cp.isExtruded()) {
848 warningEmition("Warning from traverse_loop. Loop "
849 + lp.getNumber() + " has crossed regions\n");
850 } else if ((cpnext.getStart() - cp.getEnd()) != 1) {
851 cp.setExtruded(true);
852 continue set_radius; // remplacement du goto
855 if (cp.isExtruded()) {
856 construct_extruded_segment(cp, cpnext);
858 n = cpnext.getStart() - cp.getEnd();
862 for (j = 1; j < n; j++) {
866 a = ac + j * angleinc;
867 rr = rc + (rcn - rc) * (a - ac) / dan;
868 bases.get(i).setX(xc + rr * Math.cos(a));
869 bases.get(i).setY(yc + rr * Math.sin(a));
875 for (ic = 0; ic < lp.getNconnection(); ic++) {
877 cp = lp.getConnection(ic);
879 traverse_loop(cp.getLoop(), cp);
885 for (ic = 0; ic < lp.getNconnection(); ic++) {
887 if (j >= lp.getNconnection())
889 cp = lp.getConnection(ic);
890 cpnext = lp.getConnection(j);
892 sx += bases.get(cp.getStart()).getX()
893 + bases.get(cp.getEnd()).getX();
894 sy += bases.get(cp.getStart()).getY()
895 + bases.get(cp.getEnd()).getY();
896 if (!cp.isExtruded()) {
897 for (j = cp.getEnd() + 1; j != cpnext.getStart(); j++) {
901 sx += bases.get(j).getX();
902 sy += bases.get(j).getY();
911 * For the loop pointed to by lp, determine the radius of the loop that will
912 * ensure that each base around the loop will have a separation of at least
913 * lencut around the circle. If a segment joining two connectors will not
914 * support this separation, then the flag, extruded, will be set in the
915 * first of these two indicators. The radius is set in lp.
917 * The radius is selected by a least squares procedure where the sum of the
918 * squares of the deviations of length from the ideal value of 1 is used as
919 * the error function.
921 private void determine_radius(Loop lp, double lencut) {
923 System.out.println(" Determine_radius");
924 double mindit, ci, dt, sumn, sumd, radius, dit;
925 int i, j, end, start, imindit = 0;
926 Connection cp = new Connection(), cpnext = new Connection();
927 double rt2_2 = 0.7071068;
931 for (sumd = 0.0, sumn = 0.0, i = 0; i < lp.getNconnection(); i++) {
932 cp = lp.getConnection(i);
934 if (j >= lp.getNconnection())
936 cpnext = lp.getConnection(j);
938 start = cpnext.getStart();
941 dt = cpnext.getAngle() - cp.getAngle();
944 if (!cp.isExtruded())
947 if (dt <= Math.PI / 2)
952 sumn += dt * (1.0 / ci + 1.0);
953 sumd += dt * dt / ci;
955 if (dit < mindit && !cp.isExtruded() && ci > 1.0) {
960 radius = sumn / sumd;
963 if (mindit * radius < lencut) {
964 lp.getConnection(imindit).setExtruded(true);
966 } while (mindit * radius < lencut);
967 if (lp.getRadius() > 0.0)
968 radius = lp.getRadius();
970 lp.setRadius(radius);
974 * Determines if the connections cp and cpnext are connected
976 private boolean connected_connection(Connection cp, Connection cpnext) {
978 System.out.println(" Connected_connection");
979 if (cp.isExtruded()) {
981 } else if (cp.getEnd() + 1 == cpnext.getStart()) {
989 * Finds the middle of a set of connected connectors. This is normally the
990 * middle connection in the sequence except if one of the connections is the
991 * anchor, in which case that connection will be used.
993 * @throws ExceptionNAViewAlgorithm
995 private int find_ic_middle(int icstart, int icend,
996 Connection anchor_connection, Connection acp, Loop lp)
997 throws ExceptionNAViewAlgorithm {
999 System.out.println(" Find_ic_middle");
1000 int count, ret, ic, i;
1008 if (count++ > lp.getNconnection() * 2) {
1009 throw new ExceptionNAViewAlgorithm(
1010 "Infinite loop detected in find_ic_middle");
1012 if (anchor_connection != null && lp.getConnection(ic) == acp) {
1016 if (++ic >= lp.getNconnection()) {
1021 for (i = 1, ic = icstart; i < (count + 1) / 2; i++) {
1022 if (++ic >= lp.getNconnection())
1031 * Generates the coordinates for the base pairing region of a connection
1032 * given the position of the starting base pair.
1034 * @throws ExceptionNAViewAlgorithm
1036 private void generate_region(Connection cp) throws ExceptionNAViewAlgorithm {
1038 System.out.println(" Generate_region");
1039 int l, start, end, i, mate;
1042 rp = cp.getRegion();
1044 if (cp.getStart() == rp.getStart1()) {
1045 start = rp.getStart1();
1048 start = rp.getStart2();
1051 if (bases.get(cp.getStart()).getX() > ANUM - 100.0
1052 || bases.get(cp.getEnd()).getX() > ANUM - 100.0) {
1053 throw new ExceptionNAViewAlgorithm(
1054 "Bad region passed to generate_region. Coordinates not defined.");
1056 for (i = start + 1; i <= end; i++) {
1059 bases.get(cp.getStart()).getX() + HELIX_FACTOR * l
1062 bases.get(cp.getStart()).getY() + HELIX_FACTOR * l
1064 mate = bases.get(i).getMate();
1065 bases.get(mate).setX(
1066 bases.get(cp.getEnd()).getX() + HELIX_FACTOR * l
1068 bases.get(mate).setY(
1069 bases.get(cp.getEnd()).getY() + HELIX_FACTOR * l
1076 * Draws the segment of residue between the bases numbered start through
1077 * end, where start and end are presumed to be part of a base pairing
1078 * region. They are drawn as a circle which has a chord given by the ends of
1079 * two base pairing regions defined by the connections.
1081 * @throws ExceptionNAViewAlgorithm
1083 private void construct_circle_segment(int start, int end)
1084 throws ExceptionNAViewAlgorithm {
1086 System.out.println(" Construct_circle_segment");
1087 double dx, dy, rr, midx, midy, xn, yn, nrx, nry, mx, my, a;
1090 dx = bases.get(end).getX() - bases.get(start).getX();
1091 dy = bases.get(end).getY() - bases.get(start).getY();
1092 rr = Math.sqrt(dx * dx + dy * dy);
1099 for (j = 1; j < l; j++) {
1104 bases.get(start).getX() + dx * (double) j / (double) l);
1106 bases.get(start).getY() + dy * (double) j / (double) l);
1109 find_center_for_arc((l - 1), rr);
1112 midx = bases.get(start).getX() + dx * rr / 2.0;
1113 midy = bases.get(start).getY() + dy * rr / 2.0;
1116 nrx = midx + _h * xn;
1117 nry = midy + _h * yn;
1118 mx = bases.get(start).getX() - nrx;
1119 my = bases.get(start).getY() - nry;
1120 rr = Math.sqrt(mx * mx + my * my);
1121 a = Math.atan2(my, mx);
1122 for (j = 1; j < l; j++) {
1126 bases.get(i).setX(nrx + rr * Math.cos(a + j * angleinc));
1127 bases.get(i).setY(nry + rr * Math.sin(a + j * angleinc));
1133 * Constructs the segment between cp and cpnext as a circle if possible.
1134 * However, if the segment is too large, the lines are drawn between the two
1135 * connecting regions, and bases are placed there until the connecting
1138 * @throws ExceptionNAViewAlgorithm
1140 private void construct_extruded_segment(Connection cp, Connection cpnext)
1141 throws ExceptionNAViewAlgorithm {
1143 System.out.println(" Construct_extruded_segment");
1144 double astart, aend1, aend2, aave, dx, dy, a1, a2, ac, rr, da, dac;
1145 int start, end, n, nstart, nend;
1148 astart = cp.getAngle();
1149 aend2 = aend1 = cpnext.getAngle();
1151 aend2 += 2 * Math.PI;
1152 aave = (astart + aend2) / 2.0;
1153 start = cp.getEnd();
1154 end = cpnext.getStart();
1158 da = cpnext.getAngle() - cp.getAngle();
1163 construct_circle_segment(start, end);
1165 dx = bases.get(end).getX() - bases.get(start).getX();
1166 dy = bases.get(end).getY() - bases.get(start).getY();
1167 rr = Math.sqrt(dx * dx + dy * dy);
1170 if (rr >= 1.5 && da <= Math.PI / 2) {
1173 nstart -= nbase + 1;
1177 bases.get(nstart).setX(bases.get(start).getX() + 0.5 * dx);
1178 bases.get(nstart).setY(bases.get(start).getY() + 0.5 * dy);
1179 bases.get(nend).setX(bases.get(end).getX() - 0.5 * dx);
1180 bases.get(nend).setY(bases.get(end).getY() - 0.5 * dy);
1186 construct_circle_segment(start, end);
1189 nstart -= nbase + 1;
1190 dx = bases.get(nstart).getX() - bases.get(start).getX();
1191 dy = bases.get(nstart).getY() - bases.get(start).getY();
1192 a1 = Math.atan2(dy, dx);
1203 dx = bases.get(nend).getX() - bases.get(end).getX();
1204 dy = bases.get(nend).getY() - bases.get(end).getY();
1205 a2 = Math.atan2(dy, dx);
1214 ac = minf2(aave, astart + 0.5);
1215 bases.get(nstart).setX(
1216 bases.get(start).getX() + Math.cos(ac));
1217 bases.get(nstart).setY(
1218 bases.get(start).getY() + Math.sin(ac));
1220 ac = maxf2(aave, aend2 - 0.5);
1221 bases.get(nend).setX(bases.get(end).getX() + Math.cos(ac));
1222 bases.get(nend).setY(bases.get(end).getY() + Math.sin(ac));
1226 } while (collision && n > 1);
1231 * Given n points to be placed equidistantly and equiangularly on a polygon
1232 * which has a chord of length, b, find the distance, h, from the midpoint
1233 * of the chord for the center of polygon. Positive values mean the center
1234 * is within the polygon and the chord, whereas negative values mean the
1235 * center is outside the chord. Also, the radial angle for each polygon side
1236 * is returned in theta.
1238 * The procedure uses a bisection algorithm to find the correct value for
1239 * the center. Two equations are solved, the angles around the center must
1240 * add to 2*Math.PI, and the sides of the polygon excluding the chord must
1241 * have a length of 1.
1243 * @throws ExceptionNAViewAlgorithm
1245 private void find_center_for_arc(double n, double b)
1246 throws ExceptionNAViewAlgorithm {
1248 System.out.println(" Find_center_for_arc");
1249 double h, hhi, hlow, r, disc, theta, e, phi;
1252 hhi = (n + 1.0) / Math.PI;
1253 // changed to prevent div by zero if (ih)
1254 hlow = -hhi - b / (n + 1.000001 - b);
1256 // otherwise we might fail below (ih)
1260 h = (hhi + hlow) / 2.0;
1261 r = Math.sqrt(h * h + b * b / 4.0);
1262 // if (r<0.5) {r = 0.5; h = 0.5*Math.sqrt(1-b*b);}
1263 disc = 1.0 - 0.5 / (r * r);
1264 if (Math.abs(disc) > 1.0) {
1265 throw new ExceptionNAViewAlgorithm(
1266 "Unexpected large magnitude discriminant = " + disc
1269 theta = Math.acos(disc);
1270 // theta = 2*Math.acos(Math.sqrt(1-1/(4*r*r)));
1271 phi = Math.acos(h / r);
1272 e = theta * (n + 1) + 2 * phi - 2 * Math.PI;
1278 } while (Math.abs(e) > 0.0001 && ++iter < MAXITER);
1279 if (iter >= MAXITER) {
1280 if (noIterationFailureYet) {
1281 warningEmition("Iteration failed in find_center_for_arc");
1282 noIterationFailureYet = false;
1291 private double minf2(double x1, double x2) {
1292 return ((x1) < (x2)) ? (x1) : (x2);
1295 private double maxf2(double x1, double x2) {
1296 return ((x1) > (x2)) ? (x1) : (x2);
1299 public void warningEmition(String warningMessage) throws ExceptionNAViewAlgorithm {
1300 throw (new ExceptionNAViewAlgorithm(warningMessage));