/*
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.
Copyright (C) 2008 Kevin Darty, Alain Denise and Yann Ponty.
electronic mail : Yann.Ponty@lri.fr
paper mail : LRI, bat 490 University Paris-Sud 91405 Orsay Cedex France
This file is part of VARNA version 3.1.
VARNA version 3.1 is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License
as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
VARNA version 3.1 is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with VARNA version 3.1.
If not, see http://www.gnu.org/licenses.
*/
package fr.orsay.lri.varna.models.naView;
import java.util.ArrayList;
import fr.orsay.lri.varna.exceptions.ExceptionNAViewAlgorithm;
import fr.orsay.lri.varna.interfaces.InterfaceVARNAListener;
import fr.orsay.lri.varna.interfaces.InterfaceVARNAObservable;
public class NAView {
private final double ANUM = 9999.0;
private final int MAXITER = 500;
private ArrayList<Base> bases;
private int nbase, nregion, loop_count;
private Loop root = new Loop();
private ArrayList<Loop> loops;
private ArrayList<Region> regions;
private Radloop rlphead = new Radloop();
private double lencut=0.8;
private final double RADIUS_REDUCTION_FACTOR = 1.4;
// show algorithm step by step
private boolean debug = false;
private double angleinc;
private double _h;
private ArrayList<InterfaceVARNAListener> _listeVARNAListener = new ArrayList<InterfaceVARNAListener>();
private boolean noIterationFailureYet = true;
double HELIX_FACTOR = 0.6;
double BACKBONE_DISTANCE = 27;
public int naview_xy_coordinates(ArrayList<Short> pair_table2,
ArrayList<Double> x, ArrayList<Double> y)
throws ExceptionNAViewAlgorithm {
if (debug)
System.out.println("naview_xy_coordinates");
if (pair_table2.size() == 0)
return 0;
int i;
ArrayList<Integer> pair_table = new ArrayList<Integer>(pair_table2
.size() + 1);
pair_table.add(pair_table2.size());
for (int j = 0; j < pair_table2.size(); j++) {
pair_table.add(pair_table2.get(j) + 1);
}
if (debug) {
infoStructure(pair_table);
}
// length
nbase = pair_table.get(0);
bases = new ArrayList<Base>(nbase + 1);
for (int index = 0; index < bases.size(); index++) {
bases.add(new Base());
}
regions = new ArrayList<Region>();
for (int index = 0; index < nbase + 1; index++) {
regions.add(new Region());
}
read_in_bases(pair_table);
if (debug)
infoBasesMate();
rlphead = null;
find_regions();
if (debug)
infoRegions();
loop_count = 0;
loops = new ArrayList<Loop>(nbase + 1);
for (int index = 0; index < nbase + 1; index++) {
loops.add(new Loop());
}
construct_loop(0);
if (debug)
infoBasesExtracted();
find_central_loop();
if (debug)
infoRoot();
if (debug)
dump_loops();
traverse_loop(root, null);
for (i = 0; i < nbase; i++) {
x.add(100 + BACKBONE_DISTANCE * bases.get(i + 1).getX());
y.add(100 + BACKBONE_DISTANCE * bases.get(i + 1).getY());
}
return nbase;
}
private void infoStructure(ArrayList<Integer> pair_table) {
System.out.println("structure:");
for (int j = 0; j < pair_table.size(); j++) {
System.out.print("#" + j + ":" + pair_table.get(j) + "\t");
if (j % 10 == 0)
System.out.println();
}
System.out.println();
}
private void infoBasesMate() {
System.out.println("Bases mate:");
for (int index = 0; index < bases.size(); index++) {
System.out.print("#" + index + ":" + bases.get(index).getMate()
+ "\t");
if (index % 10 == 0)
System.out.println();
}
System.out.println();
}
private void infoRegions() {
System.out.println("regions:");
for (int index = 0; index < regions.size(); index++) {
System.out.print("(" + regions.get(index).getStart1() + ","
+ regions.get(index).getStart2() + ";"
+ regions.get(index).getEnd1() + ","
+ regions.get(index).getEnd2() + ")\t\t");
if (index % 5 == 0)
System.out.println();
}
System.out.println();
}
private void infoBasesExtracted() {
System.out.println("Bases extracted:");
for (int index = 0; index < bases.size(); index++) {
System.out.print("i=" + index + ":"
+ bases.get(index).isExtracted() + "\t");
if (index % 5 == 0)
System.out.println();
}
System.out.println();
}
private void infoRoot() {
System.out.println("root" + root.getNconnection() + ";"
+ root.getNumber());
System.out.println("\troot : ");
System.out.println("\tdepth=" + root.getDepth());
System.out.println("\tmark=" + root.isMark());
System.out.println("\tnumber=" + root.getNumber());
System.out.println("\tradius=" + root.getRadius());
System.out.println("\tx=" + root.getX());
System.out.println("\ty=" + root.getY());
System.out.println("\tnconnection=" + root.getNconnection());
}
private void read_in_bases(ArrayList<Integer> pair_table) {
if (debug)
System.out.println("read_in_bases");
int i, npairs;
// Set up an origin.
bases.add(new Base());
bases.get(0).setMate(0);
bases.get(0).setExtracted(false);
bases.get(0).setX(ANUM);
bases.get(0).setY(ANUM);
for (npairs = 0, i = 1; i <= nbase; i++) {
bases.add(new Base());
bases.get(i).setExtracted(false);
bases.get(i).setX(ANUM);
bases.get(i).setY(ANUM);
bases.get(i).setMate(pair_table.get(i));
if ((int) pair_table.get(i) > i)
npairs++;
}
// must have at least 1 pair to avoid segfault
if (npairs == 0) {
bases.get(1).setMate(nbase);
bases.get(nbase).setMate(1);
}
}
/**
* Identifies the regions in the structure.
*/
private void find_regions()
{
if (debug)
System.out.println("find_regions");
int i, mate, nb1;
nb1 = nbase + 1;
ArrayList<Boolean> mark = new ArrayList<Boolean>(nb1);
for (i = 0; i < nb1; i++)
mark.add(false);
nregion = 0;
for (i = 0; i <= nbase; i++) {
if ((mate = bases.get(i).getMate()) != 0 && !mark.get(i)) {
regions.get(nregion).setStart1(i);
regions.get(nregion).setEnd2(mate);
mark.set(i, true);
mark.set(mate, true);
bases.get(i).setRegion(regions.get(nregion));
bases.get(mate).setRegion(regions.get(nregion));
for (i++, mate--; i < mate && bases.get(i).getMate() == mate; i++, mate--) {
mark.set(mate, true);
mark.set(i, true);
bases.get(i).setRegion(regions.get(nregion));
bases.get(mate).setRegion(regions.get(nregion));
}
regions.get(nregion).setEnd1(--i);
regions.get(nregion).setStart2(mate + 1);
if (debug) {
if (nregion == 0)
System.out.printf("\nRegions are:\n");
System.out.printf(
"Region %d is %d-%d and %d-%d with gap of %d.\n",
nregion + 1, regions.get(nregion).getStart1(),
regions.get(nregion).getEnd1(), regions
.get(nregion).getStart2(), regions.get(
nregion).getEnd2(), regions.get(nregion)
.getStart2()
- regions.get(nregion).getEnd1() + 1);
}
nregion++;
}
}
}
/**
* Starting at residue ibase, recursively constructs the loop containing
* said base and all deeper bases.
*
* @throws ExceptionNAViewAlgorithm
*/
private Loop construct_loop(int ibase) throws ExceptionNAViewAlgorithm {
if (debug)
System.out.println("construct_loop");
int i, mate;
Loop retloop = new Loop(), lp = new Loop();
Connection cp = new Connection();
Region rp = new Region();
Radloop rlp = new Radloop();
retloop = loops.get(loop_count++);
retloop.setNconnection(0);
//System.out.println(""+ibase+" "+nbase);
//ArrayList<Connection> a = new ArrayList<Connection>(nbase + 1);
//retloop.setConnections(a);
// for (int index = 0; index < nbase + 1; index++)
// retloop.getConnections().add(new Connection());
//for (int index = 0; index < nbase + 1; index++)
// retloop.addConnection(index,new Connection());
retloop.setDepth(0);
retloop.setNumber(loop_count);
retloop.setRadius(0.0);
for (rlp = rlphead; rlp != null; rlp = rlp.getNext())
if (rlp.getLoopnumber() == loop_count)
retloop.setRadius(rlp.getRadius());
i = ibase;
do {
if ((mate = bases.get(i).getMate()) != 0) {
rp = bases.get(i).getRegion();
if (!bases.get(rp.getStart1()).isExtracted()) {
if (i == rp.getStart1()) {
bases.get(rp.getStart1()).setExtracted(true);
bases.get(rp.getEnd1()).setExtracted(true);
bases.get(rp.getStart2()).setExtracted(true);
bases.get(rp.getEnd2()).setExtracted(true);
lp = construct_loop(rp.getEnd1() < nbase ? rp.getEnd1() + 1
: 0);
} else if (i == rp.getStart2()) {
bases.get(rp.getStart2()).setExtracted(true);
bases.get(rp.getEnd2()).setExtracted(true);
bases.get(rp.getStart1()).setExtracted(true);
bases.get(rp.getEnd1()).setExtracted(true);
lp = construct_loop(rp.getEnd2() < nbase ? rp.getEnd2() + 1
: 0);
} else {
throw new ExceptionNAViewAlgorithm(
"naview:Error detected in construct_loop. i = "
+ i + " not found in region table.\n");
}
retloop.setNconnection(retloop.getNconnection() + 1);
cp = new Connection();
retloop.setConnection(retloop.getNconnection() - 1, cp);
retloop.setConnection(retloop.getNconnection(), null);
cp.setLoop(lp);
cp.setRegion(rp);
if (i == rp.getStart1()) {
cp.setStart(rp.getStart1());
cp.setEnd(rp.getEnd2());
} else {
cp.setStart(rp.getStart2());
cp.setEnd(rp.getEnd1());
}
cp.setExtruded(false);
cp.setBroken(false);
lp.setNconnection(lp.getNconnection() + 1);
cp = new Connection();
lp.setConnection(lp.getNconnection() - 1, cp);
lp.setConnection(lp.getNconnection(), null);
cp.setLoop(retloop);
cp.setRegion(rp);
if (i == rp.getStart1()) {
cp.setStart(rp.getStart2());
cp.setEnd(rp.getEnd1());
} else {
cp.setStart(rp.getStart1());
cp.setEnd(rp.getEnd2());
}
cp.setExtruded(false);
cp.setBroken(false);
}
i = mate;
}
if (++i > nbase)
i = 0;
} while (i != ibase);
return retloop;
}
/**
* Displays all the loops.
*/
private void dump_loops() {
System.out.println("dump_loops");
int il, ilp, irp;
Loop lp;
Connection cp;
System.out.printf("\nRoot loop is #%d\n", loops.indexOf(root) + 1);
for (il = 0; il < loop_count; il++) {
lp = loops.get(il);
System.out.printf("Loop %d has %d connections:\n", il + 1, lp
.getNconnection());
for (int i = 0; (cp = lp.getConnection(i)) != null; i++) {
ilp = (loops.indexOf(cp.getLoop())) + 1;
irp = (regions.indexOf(cp.getRegion())) + 1;
System.out.printf(" Loop %d Region %d (%d-%d)\n", ilp, irp, cp
.getStart(), cp.getEnd());
}
}
}
/**
* Find node of greatest branching that is deepest.
*/
private void find_central_loop() {
if (debug)
System.out.println("find_central_loop");
Loop lp = new Loop();
int maxconn, maxdepth, i;
determine_depths();
maxconn = 0;
maxdepth = -1;
for (i = 0; i < loop_count; i++) {
lp = loops.get(i);
if (lp.getNconnection() > maxconn) {
maxdepth = lp.getDepth();
maxconn = lp.getNconnection();
root = lp;
} else if (lp.getDepth() > maxdepth
&& lp.getNconnection() == maxconn) {
maxdepth = lp.getDepth();
root = lp;
}
}
}
/**
* Determine the depth of all loops.
*/
private void determine_depths() {
if (debug)
System.out.println("determine_depths");
Loop lp = new Loop();
int i, j;
for (i = 0; i < loop_count; i++) {
lp = loops.get(i);
for (j = 0; j < loop_count; j++)
loops.get(j).setMark(false);
lp.setDepth(depth(lp));
}
}
/**
* Determines the depth of loop, lp. Depth is defined as the minimum
* distance to a leaf loop where a leaf loop is one that has only one or no
* connections.
*/
private int depth(Loop lp) {
if (debug)
System.out.println("depth");
int count, ret, d;
if (lp.getNconnection() <= 1)
return 0;
if (lp.isMark())
return -1;
lp.setMark(true);
count = 0;
ret = 0;
for (int i = 0; lp.getConnection(i) != null; i++) {
d = depth(lp.getConnection(i).getLoop());
if (d >= 0) {
if (++count == 1)
ret = d;
else if (ret > d)
ret = d;
}
}
lp.setMark(false);
return ret + 1;
}
/**
* This is the workhorse of the display program. The algorithm is recursive
* based on processing individual loops. Each base pairing region is
* displayed using the direction given by the circle diagram, and the
* connections between the regions is drawn by equally spaced points. The
* radius of the loop is set to minimize the square error for lengths
* between sequential bases in the loops. The "correct" length for base
* links is 1. If the least squares fitting of the radius results in loops
* being less than 1/2 unit apart, then that segment is extruded.
*
* The variable, anchor_connection, gives the connection to the loop
* processed in an previous level of recursion.
*
* @throws ExceptionNAViewAlgorithm
*/
private void traverse_loop(Loop lp, Connection anchor_connection)
throws ExceptionNAViewAlgorithm {
if (debug)
System.out.println(" traverse_loop");
double xs, ys, xe, ye, xn, yn, angleinc, r;
double radius, xc, yc, xo, yo, astart, aend, a;
Connection cp, cpnext, acp, cpprev;
int i, j, n, ic;
double da, maxang;
int count, icstart, icend, icmiddle, icroot;
boolean done, done_all_connections, rooted;
int sign;
double midx, midy, nrx, nry, mx, my, vx, vy, dotmv, nmidx, nmidy;
int icstart1, icup, icdown, icnext, direction;
double dan, dx, dy, rr;
double cpx, cpy, cpnextx, cpnexty, cnx, cny, rcn, rc, lnx, lny, rl, ac, acn, sx, sy, dcp;
int imaxloop = 0;
angleinc = 2 * Math.PI / (nbase + 1);
acp = null;
icroot = -1;
int indice = 0;
for (ic = 0; (cp = lp.getConnection(indice)) != null; indice++, ic++) {
// xs = cos(angleinc*cp.setStart(); ys = sin(angleinc*cp.setStart();
// xe =
// cos(angleinc*cp.setEnd()); ye = sin(angleinc*cp.setEnd());
xs = -Math.sin(angleinc * cp.getStart());
ys = Math.cos(angleinc * cp.getStart());
xe = -Math.sin(angleinc * cp.getEnd());
ye = Math.cos(angleinc * cp.getEnd());
xn = ye - ys;
yn = xs - xe;
r = Math.sqrt(xn * xn + yn * yn);
cp.setXrad(xn / r);
cp.setYrad(yn / r);
cp.setAngle(Math.atan2(yn, xn));
if (cp.getAngle() < 0.0)
cp.setAngle(cp.getAngle() + 2 * Math.PI);
if (anchor_connection != null
&& anchor_connection.getRegion() == cp.getRegion()) {
acp = cp;
icroot = ic;
}
}
// remplacement d'une etiquette de goto
set_radius: while (true) {
determine_radius(lp, lencut);
radius = lp.getRadius()/RADIUS_REDUCTION_FACTOR;
if (anchor_connection == null)
xc = yc = 0.0;
else {
xo = (bases.get(acp.getStart()).getX() + bases
.get(acp.getEnd()).getX()) / 2.0;
yo = (bases.get(acp.getStart()).getY() + bases
.get(acp.getEnd()).getY()) / 2.0;
xc = xo - radius * acp.getXrad();
yc = yo - radius * acp.getYrad();
}
// The construction of the connectors will proceed in blocks of
// connected connectors, where a connected connector pairs means two
// connectors that are forced out of the drawn circle because they
// are too close together in angle.
// First, find the start of a block of connected connectors
if (icroot == -1)
icstart = 0;
else
icstart = icroot;
cp = lp.getConnection(icstart);
count = 0;
if (debug)
{
System.out.printf("Now processing loop %d\n", lp.getNumber());
System.out.println(" "+lp);
}
done = false;
do {
j = icstart - 1;
if (j < 0)
j = lp.getNconnection() - 1;
cpprev = lp.getConnection(j);
if (!connected_connection(cpprev, cp)) {
done = true;
} else {
icstart = j;
cp = cpprev;
}
if (++count > lp.getNconnection()) {
// Here everything is connected. Break on maximum angular
// separation between connections.
maxang = -1.0;
for (ic = 0; ic < lp.getNconnection(); ic++) {
j = ic + 1;
if (j >= lp.getNconnection())
j = 0;
cp = lp.getConnection(ic);
cpnext = lp.getConnection(j);
ac = cpnext.getAngle() - cp.getAngle();
if (ac < 0.0)
ac += 2 * Math.PI;
if (ac > maxang) {
maxang = ac;
imaxloop = ic;
}
}
icend = imaxloop;
icstart = imaxloop + 1;
if (icstart >= lp.getNconnection())
icstart = 0;
cp = lp.getConnection(icend);
cp.setBroken(true);
done = true;
}
} while (!done);
done_all_connections = false;
icstart1 = icstart;
if (debug)
System.out.printf(" Icstart1 = %d\n", icstart1);
while (!done_all_connections) {
count = 0;
done = false;
icend = icstart;
rooted = false;
while (!done) {
cp = lp.getConnection(icend);
if (icend == icroot)
rooted = true;
j = icend + 1;
if (j >= lp.getNconnection()) {
j = 0;
}
cpnext = lp.getConnection(j);
if (connected_connection(cp, cpnext)) {
if (++count >= lp.getNconnection())
break;
icend = j;
} else {
done = true;
}
}
icmiddle = find_ic_middle(icstart, icend, anchor_connection,
acp, lp);
ic = icup = icdown = icmiddle;
if (debug)
System.out.printf(" IC start = %d middle = %d end = %d\n",
icstart, icmiddle, icend);
done = false;
direction = 0;
while (!done) {
if (direction < 0) {
ic = icup;
} else if (direction == 0) {
ic = icmiddle;
} else {
ic = icdown;
}
if (ic >= 0) {
cp = lp.getConnection(ic);
if (anchor_connection == null || acp != cp) {
if (direction == 0) {
astart = cp.getAngle()
- Math.asin(1.0 / 2.0 / radius);
aend = cp.getAngle()
+ Math.asin(1.0 / 2.0 / radius);
bases.get(cp.getStart()).setX(
xc + radius * Math.cos(astart));
bases.get(cp.getStart()).setY(
yc + radius * Math.sin(astart));
bases.get(cp.getEnd()).setX(
xc + radius * Math.cos(aend));
bases.get(cp.getEnd()).setY(
yc + radius * Math.sin(aend));
} else if (direction < 0) {
j = ic + 1;
if (j >= lp.getNconnection())
j = 0;
cp = lp.getConnection(ic);
cpnext = lp.getConnection(j);
cpx = cp.getXrad();
cpy = cp.getYrad();
ac = (cp.getAngle() + cpnext.getAngle()) / 2.0;
if (cp.getAngle() > cpnext.getAngle())
ac -= Math.PI;
cnx = Math.cos(ac);
cny = Math.sin(ac);
lnx = cny;
lny = -cnx;
da = cpnext.getAngle() - cp.getAngle();
if (da < 0.0)
da += 2 * Math.PI;
if (cp.isExtruded()) {
if (da <= Math.PI / 2)
rl = 2.0;
else
rl = 1.5;
} else {
rl = 1.0;
}
bases.get(cp.getEnd()).setX(
bases.get(cpnext.getStart()).getX()
+ rl * lnx);
bases.get(cp.getEnd()).setY(
bases.get(cpnext.getStart()).getY()
+ rl * lny);
bases.get(cp.getStart()).setX(
bases.get(cp.getEnd()).getX() + cpy);
bases.get(cp.getStart()).setY(
bases.get(cp.getEnd()).getY() - cpx);
} else {
j = ic - 1;
if (j < 0)
j = lp.getNconnection() - 1;
cp = lp.getConnection(j);
cpnext = lp.getConnection(ic);
cpnextx = cpnext.getXrad();
cpnexty = cpnext.getYrad();
ac = (cp.getAngle() + cpnext.getAngle()) / 2.0;
if (cp.getAngle() > cpnext.getAngle())
ac -= Math.PI;
cnx = Math.cos(ac);
cny = Math.sin(ac);
lnx = -cny;
lny = cnx;
da = cpnext.getAngle() - cp.getAngle();
if (da < 0.0)
da += 2 * Math.PI;
if (cp.isExtruded()) {
if (da <= Math.PI / 2)
rl = 2.0;
else
rl = 1.5;
} else {
rl = 1.0;
}
bases.get(cpnext.getStart()).setX(
bases.get(cp.getEnd()).getX() + rl
* lnx);
bases.get(cpnext.getStart()).setY(
bases.get(cp.getEnd()).getY() + rl
* lny);
bases.get(cpnext.getEnd()).setX(
bases.get(cpnext.getStart()).getX()
- cpnexty);
bases.get(cpnext.getEnd()).setY(
bases.get(cpnext.getStart()).getY()
+ cpnextx);
}
}
}
if (direction < 0) {
if (icdown == icend) {
icdown = -1;
} else if (icdown >= 0) {
if (++icdown >= lp.getNconnection()) {
icdown = 0;
}
}
direction = 1;
} else {
if (icup == icstart)
icup = -1;
else if (icup >= 0) {
if (--icup < 0) {
icup = lp.getNconnection() - 1;
}
}
direction = -1;
}
done = icup == -1 && icdown == -1;
}
icnext = icend + 1;
if (icnext >= lp.getNconnection())
icnext = 0;
if (icend != icstart
&& (!(icstart == icstart1 && icnext == icstart1))) {
// Move the bases just constructed (or the radius) so that
// the bisector of the end points is radius distance away
// from the loop center.
cp = lp.getConnection(icstart);
cpnext = lp.getConnection(icend);
dx = bases.get(cpnext.getEnd()).getX()
- bases.get(cp.getStart()).getX();
dy = bases.get(cpnext.getEnd()).getY()
- bases.get(cp.getStart()).getY();
midx = bases.get(cp.getStart()).getX() + dx / 2.0;
midy = bases.get(cp.getStart()).getY() + dy / 2.0;
rr = Math.sqrt(dx * dx + dy * dy);
mx = dx / rr;
my = dy / rr;
vx = xc - midx;
vy = yc - midy;
rr = Math.sqrt(dx * dx + dy * dy);
vx /= rr;
vy /= rr;
dotmv = vx * mx + vy * my;
nrx = dotmv * mx - vx;
nry = dotmv * my - vy;
rr = Math.sqrt(nrx * nrx + nry * nry);
nrx /= rr;
nry /= rr;
// Determine which side of the bisector the center should
// be.
dx = bases.get(cp.getStart()).getX() - xc;
dy = bases.get(cp.getStart()).getY() - yc;
ac = Math.atan2(dy, dx);
if (ac < 0.0)
ac += 2 * Math.PI;
dx = bases.get(cpnext.getEnd()).getX() - xc;
dy = bases.get(cpnext.getEnd()).getY() - yc;
acn = Math.atan2(dy, dx);
if (acn < 0.0)
acn += 2 * Math.PI;
if (acn < ac)
acn += 2 * Math.PI;
if (acn - ac > Math.PI)
sign = -1;
else
sign = 1;
nmidx = xc + sign * radius * nrx;
nmidy = yc + sign * radius * nry;
if (rooted) {
xc -= nmidx - midx;
yc -= nmidy - midy;
} else {
for (ic = icstart;;) {
cp = lp.getConnection(ic);
i = cp.getStart();
bases.get(i).setX(
bases.get(i).getX() + nmidx - midx);
bases.get(i).setY(
bases.get(i).getY() + nmidy - midy);
i = cp.getEnd();
bases.get(i).setX(
bases.get(i).getX() + nmidx - midx);
bases.get(i).setY(
bases.get(i).getY() + nmidy - midy);
if (ic == icend)
break;
if (++ic >= lp.getNconnection())
ic = 0;
}
}
}
icstart = icnext;
done_all_connections = icstart == icstart1;
}
for (ic = 0; ic < lp.getNconnection(); ic++) {
cp = lp.getConnection(ic);
j = ic + 1;
if (j >= lp.getNconnection())
j = 0;
cpnext = lp.getConnection(j);
dx = bases.get(cp.getEnd()).getX() - xc;
dy = bases.get(cp.getEnd()).getY() - yc;
rc = Math.sqrt(dx * dx + dy * dy);
ac = Math.atan2(dy, dx);
if (ac < 0.0)
ac += 2 * Math.PI;
dx = bases.get(cpnext.getStart()).getX() - xc;
dy = bases.get(cpnext.getStart()).getY() - yc;
rcn = Math.sqrt(dx * dx + dy * dy);
acn = Math.atan2(dy, dx);
if (acn < 0.0)
acn += 2 * Math.PI;
if (acn < ac)
acn += 2 * Math.PI;
dan = acn - ac;
dcp = cpnext.getAngle() - cp.getAngle();
if (dcp <= 0.0)
dcp += 2 * Math.PI;
if (Math.abs(dan - dcp) > Math.PI) {
if (cp.isExtruded()) {
warningEmition("Warning from traverse_loop. Loop "
+ lp.getNumber() + " has crossed regions\n");
} else if ((cpnext.getStart() - cp.getEnd()) != 1) {
cp.setExtruded(true);
continue set_radius; // remplacement du goto
}
}
if (cp.isExtruded()) {
construct_extruded_segment(cp, cpnext);
} else {
n = cpnext.getStart() - cp.getEnd();
if (n < 0)
n += nbase + 1;
angleinc = dan / n;
for (j = 1; j < n; j++) {
i = cp.getEnd() + j;
if (i > nbase)
i -= nbase + 1;
a = ac + j * angleinc;
rr = rc + (rcn - rc) * (a - ac) / dan;
bases.get(i).setX(xc + rr * Math.cos(a));
bases.get(i).setY(yc + rr * Math.sin(a));
}
}
}
break;
}
for (ic = 0; ic < lp.getNconnection(); ic++) {
if (icroot != ic) {
cp = lp.getConnection(ic);
generate_region(cp);
traverse_loop(cp.getLoop(), cp);
}
}
n = 0;
sx = 0.0;
sy = 0.0;
for (ic = 0; ic < lp.getNconnection(); ic++) {
j = ic + 1;
if (j >= lp.getNconnection())
j = 0;
cp = lp.getConnection(ic);
cpnext = lp.getConnection(j);
n += 2;
sx += bases.get(cp.getStart()).getX()
+ bases.get(cp.getEnd()).getX();
sy += bases.get(cp.getStart()).getY()
+ bases.get(cp.getEnd()).getY();
if (!cp.isExtruded()) {
for (j = cp.getEnd() + 1; j != cpnext.getStart(); j++) {
if (j > nbase)
j -= nbase + 1;
n++;
sx += bases.get(j).getX();
sy += bases.get(j).getY();
}
}
}
lp.setX(sx / n);
lp.setY(sy / n);
}
/**
* For the loop pointed to by lp, determine the radius of the loop that will
* ensure that each base around the loop will have a separation of at least
* lencut around the circle. If a segment joining two connectors will not
* support this separation, then the flag, extruded, will be set in the
* first of these two indicators. The radius is set in lp.
*
* The radius is selected by a least squares procedure where the sum of the
* squares of the deviations of length from the ideal value of 1 is used as
* the error function.
*/
private void determine_radius(Loop lp, double lencut) {
if (debug)
System.out.println(" Determine_radius");
double mindit, ci, dt, sumn, sumd, radius, dit;
int i, j, end, start, imindit = 0;
Connection cp = new Connection(), cpnext = new Connection();
double rt2_2 = 0.7071068;
do {
mindit = 1.0e10;
for (sumd = 0.0, sumn = 0.0, i = 0; i < lp.getNconnection(); i++) {
cp = lp.getConnection(i);
j = i + 1;
if (j >= lp.getNconnection())
j = 0;
cpnext = lp.getConnection(j);
end = cp.getEnd();
start = cpnext.getStart();
if (start < end)
start += nbase + 1;
dt = cpnext.getAngle() - cp.getAngle();
if (dt <= 0.0)
dt += 2 * Math.PI;
if (!cp.isExtruded())
ci = start - end;
else {
if (dt <= Math.PI / 2)
ci = 2.0;
else
ci = 1.5;
}
sumn += dt * (1.0 / ci + 1.0);
sumd += dt * dt / ci;
dit = dt / ci;
if (dit < mindit && !cp.isExtruded() && ci > 1.0) {
mindit = dit;
imindit = i;
}
}
radius = sumn / sumd;
if (radius < rt2_2)
radius = rt2_2;
if (mindit * radius < lencut) {
lp.getConnection(imindit).setExtruded(true);
}
} while (mindit * radius < lencut);
if (lp.getRadius() > 0.0)
radius = lp.getRadius();
else
lp.setRadius(radius);
}
/**
* Determines if the connections cp and cpnext are connected
*/
private boolean connected_connection(Connection cp, Connection cpnext) {
if (debug)
System.out.println(" Connected_connection");
if (cp.isExtruded()) {
return true;
} else if (cp.getEnd() + 1 == cpnext.getStart()) {
return true;
} else {
return false;
}
}
/**
* Finds the middle of a set of connected connectors. This is normally the
* middle connection in the sequence except if one of the connections is the
* anchor, in which case that connection will be used.
*
* @throws ExceptionNAViewAlgorithm
*/
private int find_ic_middle(int icstart, int icend,
Connection anchor_connection, Connection acp, Loop lp)
throws ExceptionNAViewAlgorithm {
if (debug)
System.out.println(" Find_ic_middle");
int count, ret, ic, i;
boolean done;
count = 0;
ret = -1;
ic = icstart;
done = false;
while (!done) {
if (count++ > lp.getNconnection() * 2) {
throw new ExceptionNAViewAlgorithm(
"Infinite loop detected in find_ic_middle");
}
if (anchor_connection != null && lp.getConnection(ic) == acp) {
ret = ic;
}
done = ic == icend;
if (++ic >= lp.getNconnection()) {
ic = 0;
}
}
if (ret == -1) {
for (i = 1, ic = icstart; i < (count + 1) / 2; i++) {
if (++ic >= lp.getNconnection())
ic = 0;
}
ret = ic;
}
return ret;
}
/**
* Generates the coordinates for the base pairing region of a connection
* given the position of the starting base pair.
*
* @throws ExceptionNAViewAlgorithm
*/
private void generate_region(Connection cp) throws ExceptionNAViewAlgorithm {
if (debug)
System.out.println(" Generate_region");
int l, start, end, i, mate;
Region rp;
rp = cp.getRegion();
l = 0;
if (cp.getStart() == rp.getStart1()) {
start = rp.getStart1();
end = rp.getEnd1();
} else {
start = rp.getStart2();
end = rp.getEnd2();
}
if (bases.get(cp.getStart()).getX() > ANUM - 100.0
|| bases.get(cp.getEnd()).getX() > ANUM - 100.0) {
throw new ExceptionNAViewAlgorithm(
"Bad region passed to generate_region. Coordinates not defined.");
}
for (i = start + 1; i <= end; i++) {
l++;
bases.get(i).setX(
bases.get(cp.getStart()).getX() + HELIX_FACTOR * l
* cp.getXrad());
bases.get(i).setY(
bases.get(cp.getStart()).getY() + HELIX_FACTOR * l
* cp.getYrad());
mate = bases.get(i).getMate();
bases.get(mate).setX(
bases.get(cp.getEnd()).getX() + HELIX_FACTOR * l
* cp.getXrad());
bases.get(mate).setY(
bases.get(cp.getEnd()).getY() + HELIX_FACTOR * l
* cp.getYrad());
}
}
/**
* Draws the segment of residue between the bases numbered start through
* end, where start and end are presumed to be part of a base pairing
* region. They are drawn as a circle which has a chord given by the ends of
* two base pairing regions defined by the connections.
*
* @throws ExceptionNAViewAlgorithm
*/
private void construct_circle_segment(int start, int end)
throws ExceptionNAViewAlgorithm {
if (debug)
System.out.println(" Construct_circle_segment");
double dx, dy, rr, midx, midy, xn, yn, nrx, nry, mx, my, a;
int l, j, i;
dx = bases.get(end).getX() - bases.get(start).getX();
dy = bases.get(end).getY() - bases.get(start).getY();
rr = Math.sqrt(dx * dx + dy * dy);
l = end - start;
if (l < 0)
l += nbase + 1;
if (rr >= l) {
dx /= rr;
dy /= rr;
for (j = 1; j < l; j++) {
i = start + j;
if (i > nbase)
i -= nbase + 1;
bases.get(i).setX(
bases.get(start).getX() + dx * (double) j / (double) l);
bases.get(i).setY(
bases.get(start).getY() + dy * (double) j / (double) l);
}
} else {
find_center_for_arc((l - 1), rr);
dx /= rr;
dy /= rr;
midx = bases.get(start).getX() + dx * rr / 2.0;
midy = bases.get(start).getY() + dy * rr / 2.0;
xn = dy;
yn = -dx;
nrx = midx + _h * xn;
nry = midy + _h * yn;
mx = bases.get(start).getX() - nrx;
my = bases.get(start).getY() - nry;
rr = Math.sqrt(mx * mx + my * my);
a = Math.atan2(my, mx);
for (j = 1; j < l; j++) {
i = start + j;
if (i > nbase)
i -= nbase + 1;
bases.get(i).setX(nrx + rr * Math.cos(a + j * angleinc));
bases.get(i).setY(nry + rr * Math.sin(a + j * angleinc));
}
}
}
/**
* Constructs the segment between cp and cpnext as a circle if possible.
* However, if the segment is too large, the lines are drawn between the two
* connecting regions, and bases are placed there until the connecting
* circle will fit.
*
* @throws ExceptionNAViewAlgorithm
*/
private void construct_extruded_segment(Connection cp, Connection cpnext)
throws ExceptionNAViewAlgorithm {
if (debug)
System.out.println(" Construct_extruded_segment");
double astart, aend1, aend2, aave, dx, dy, a1, a2, ac, rr, da, dac;
int start, end, n, nstart, nend;
boolean collision;
astart = cp.getAngle();
aend2 = aend1 = cpnext.getAngle();
if (aend2 < astart)
aend2 += 2 * Math.PI;
aave = (astart + aend2) / 2.0;
start = cp.getEnd();
end = cpnext.getStart();
n = end - start;
if (n < 0)
n += nbase + 1;
da = cpnext.getAngle() - cp.getAngle();
if (da < 0.0) {
da += 2 * Math.PI;
}
if (n == 2)
construct_circle_segment(start, end);
else {
dx = bases.get(end).getX() - bases.get(start).getX();
dy = bases.get(end).getY() - bases.get(start).getY();
rr = Math.sqrt(dx * dx + dy * dy);
dx /= rr;
dy /= rr;
if (rr >= 1.5 && da <= Math.PI / 2) {
nstart = start + 1;
if (nstart > nbase)
nstart -= nbase + 1;
nend = end - 1;
if (nend < 0)
nend += nbase + 1;
bases.get(nstart).setX(bases.get(start).getX() + 0.5 * dx);
bases.get(nstart).setY(bases.get(start).getY() + 0.5 * dy);
bases.get(nend).setX(bases.get(end).getX() - 0.5 * dx);
bases.get(nend).setY(bases.get(end).getY() - 0.5 * dy);
start = nstart;
end = nend;
}
do {
collision = false;
construct_circle_segment(start, end);
nstart = start + 1;
if (nstart > nbase)
nstart -= nbase + 1;
dx = bases.get(nstart).getX() - bases.get(start).getX();
dy = bases.get(nstart).getY() - bases.get(start).getY();
a1 = Math.atan2(dy, dx);
if (a1 < 0.0)
a1 += 2 * Math.PI;
dac = a1 - astart;
if (dac < 0.0)
dac += 2 * Math.PI;
if (dac > Math.PI)
collision = true;
nend = end - 1;
if (nend < 0)
nend += nbase + 1;
dx = bases.get(nend).getX() - bases.get(end).getX();
dy = bases.get(nend).getY() - bases.get(end).getY();
a2 = Math.atan2(dy, dx);
if (a2 < 0.0)
a2 += 2 * Math.PI;
dac = aend1 - a2;
if (dac < 0.0)
dac += 2 * Math.PI;
if (dac > Math.PI)
collision = true;
if (collision) {
ac = minf2(aave, astart + 0.5);
bases.get(nstart).setX(
bases.get(start).getX() + Math.cos(ac));
bases.get(nstart).setY(
bases.get(start).getY() + Math.sin(ac));
start = nstart;
ac = maxf2(aave, aend2 - 0.5);
bases.get(nend).setX(bases.get(end).getX() + Math.cos(ac));
bases.get(nend).setY(bases.get(end).getY() + Math.sin(ac));
end = nend;
n -= 2;
}
} while (collision && n > 1);
}
}
/**
* Given n points to be placed equidistantly and equiangularly on a polygon
* which has a chord of length, b, find the distance, h, from the midpoint
* of the chord for the center of polygon. Positive values mean the center
* is within the polygon and the chord, whereas negative values mean the
* center is outside the chord. Also, the radial angle for each polygon side
* is returned in theta.
*
* The procedure uses a bisection algorithm to find the correct value for
* the center. Two equations are solved, the angles around the center must
* add to 2*Math.PI, and the sides of the polygon excluding the chord must
* have a length of 1.
*
* @throws ExceptionNAViewAlgorithm
*/
private void find_center_for_arc(double n, double b)
throws ExceptionNAViewAlgorithm {
if (debug)
System.out.println(" Find_center_for_arc");
double h, hhi, hlow, r, disc, theta, e, phi;
int iter;
hhi = (n + 1.0) / Math.PI;
// changed to prevent div by zero if (ih)
hlow = -hhi - b / (n + 1.000001 - b);
if (b < 1)
// otherwise we might fail below (ih)
hlow = 0;
iter = 0;
do {
h = (hhi + hlow) / 2.0;
r = Math.sqrt(h * h + b * b / 4.0);
// if (r<0.5) {r = 0.5; h = 0.5*Math.sqrt(1-b*b);}
disc = 1.0 - 0.5 / (r * r);
if (Math.abs(disc) > 1.0) {
throw new ExceptionNAViewAlgorithm(
"Unexpected large magnitude discriminant = " + disc
+ " " + r);
}
theta = Math.acos(disc);
// theta = 2*Math.acos(Math.sqrt(1-1/(4*r*r)));
phi = Math.acos(h / r);
e = theta * (n + 1) + 2 * phi - 2 * Math.PI;
if (e > 0.0) {
hlow = h;
} else {
hhi = h;
}
} while (Math.abs(e) > 0.0001 && ++iter < MAXITER);
if (iter >= MAXITER) {
if (noIterationFailureYet) {
warningEmition("Iteration failed in find_center_for_arc");
noIterationFailureYet = false;
}
h = 0.0;
theta = 0.0;
}
_h = h;
angleinc = theta;
}
private double minf2(double x1, double x2) {
return ((x1) < (x2)) ? (x1) : (x2);
}
private double maxf2(double x1, double x2) {
return ((x1) > (x2)) ? (x1) : (x2);
}
public void warningEmition(String warningMessage) throws ExceptionNAViewAlgorithm {
throw (new ExceptionNAViewAlgorithm(warningMessage));
}
}