+/*
+ 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 Université 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.controlers;
+
+import java.awt.event.ActionEvent;
+import java.awt.event.ActionListener;
+import java.awt.geom.Point2D;
+import java.util.ArrayList;
+import java.util.Arrays;
+import java.util.LinkedList;
+import java.util.Vector;
+
+import javax.swing.Timer;
+
+import fr.orsay.lri.varna.VARNAPanel;
+import fr.orsay.lri.varna.exceptions.MappingException;
+import fr.orsay.lri.varna.models.VARNAConfig;
+import fr.orsay.lri.varna.models.rna.Mapping;
+import fr.orsay.lri.varna.models.rna.ModeleBase;
+import fr.orsay.lri.varna.models.rna.RNA;
+
+public class ControleurInterpolator extends Thread implements ActionListener {
+
+ private static final int START = 0;
+ private static final int LOOP = 1;
+ private static final int END = 2;
+
+ VARNAPanel _vpn;
+ protected int _numSteps = 25; // BH SwingJS
+ private long _timeDelay = /** @j2sNative 2 || */15;
+ private boolean _running = false;
+ Targets _d = new Targets();
+ protected int _step; // BH SwingJS
+ protected boolean _firstHalf; // BH SwingJS
+
+ public ControleurInterpolator(VARNAPanel vpn) {
+ _vpn = vpn;
+ }
+
+ public synchronized void addTarget(RNA target, Mapping mapping) {
+ addTarget(target, null, mapping);
+ }
+
+
+ public synchronized void addTarget(RNA target, VARNAConfig conf, Mapping mapping) {
+ _d.add(new TargetsHolder(target,conf,mapping));
+ }
+
+ public synchronized void addTarget(RNA target) {
+ addTarget(target, null, Mapping.DefaultOutermostMapping(_vpn.getRNA()
+ .get_listeBases().size(), target.get_listeBases().size()));
+ }
+
+ public boolean isInterpolationInProgress()
+ {
+ return _running;
+ }
+
+ private Point2D.Double computeDestination(Point2D pli, Point2D pri,
+ Point2D pi, int i, int n, Point2D plf, Point2D prf) {
+ Point2D.Double plm = new Point2D.Double(
+ (pli.getX() + plf.getX()) / 2.0,
+ (pli.getY() + plf.getY()) / 2.0);
+ Point2D.Double prm = new Point2D.Double(
+ (pri.getX() + prf.getX()) / 2.0,
+ (pri.getY() + prf.getY()) / 2.0);
+ Point2D.Double pm = new Point2D.Double(((n - i) * plm.getX() + i
+ * prm.getX())
+ / n, ((n - i) * plm.getY() + i * prm.getY()) / n);
+ Point2D.Double v = new Point2D.Double(pm.getX() - pi.getX(), pm.getY()
+ - pi.getY());
+ Point2D.Double pf = new Point2D.Double(pi.getX() + 2.0 * v.getX(), pi
+ .getY()
+ + 2.0 * v.getY());
+ return pf;
+ }
+
+ private Vector<Vector<Integer>> clusterIndices(int numIndices,
+ int[] mappedIndices) throws MappingException {
+ int[] indices = new int[numIndices];
+ for (int i = 0; i < numIndices; i++) {
+ indices[i] = i;
+ }
+ return clusterIndices(indices, mappedIndices);
+ }
+
+ /**
+ * Builds and returns an interval array, alternating unmatched regions and
+ * matched indices. Namely, returns a vector of vector such that the vectors
+ * found at odd positions contain those indices that are NOT associated with
+ * any other base in the current mapping. On the other hand, vectors found
+ * at even positions contain only one mapped index. <br/>
+ * Ex: If indices=<code>[1,2,3,4,5]</code> and mappedIndices=
+ * <code>[2,5]</code> then the function will return
+ * <code>[[1],[2],[3,4],[5],[]]</code>.
+ *
+ * @param indices
+ * The total list of indices
+ * @param mappedIndices
+ * Matched indices, should be a subset of <code>indices</code>
+ * @return A clustered array
+ * @throws MappingException
+ * If one of the parameters is an empty array
+ */
+ private Vector<Vector<Integer>> clusterIndices(int[] indices,
+ int[] mappedIndices) throws MappingException {
+ if ((mappedIndices.length == 0) || (indices.length == 0)) {
+ throw new MappingException(
+ "Mapping Error: Cannot cluster indices in an empty mapping");
+ }
+ Vector<Vector<Integer>> res = new Vector<Vector<Integer>>();
+
+ Arrays.sort(indices);
+ Arrays.sort(mappedIndices);
+ int i, j = 0, k;
+ Vector<Integer> tmp = new Vector<Integer>();
+ for (i = 0; (i < indices.length) && (j < mappedIndices.length); i++) {
+ if (indices[i] == mappedIndices[j]) {
+ res.add(tmp);
+ tmp = new Vector<Integer>();
+ tmp.add(indices[i]);
+ res.add(tmp);
+ tmp = new Vector<Integer>();
+ j++;
+ } else {
+ tmp.add(indices[i]);
+ }
+ }
+ k = i;
+ for (i = k; (i < indices.length); i++) {
+ tmp.add(indices[i]);
+ }
+ res.add(tmp);
+ return res;
+ }
+
+
+
+ public void run() {
+ while (true)
+ {
+ TargetsHolder d = _d.get();
+ _running = true;
+ try{
+ nextTarget(d.target,d.conf,d.mapping);
+ }
+ catch(Exception e)
+ {
+ System.err.println(e);
+ e.printStackTrace();
+ }
+ _running = false;
+ /**
+ * @j2sNative
+ *
+ * break;
+ */
+ }
+
+ }
+
+ /**
+ * Compute the centroid of the RNA bases that have their indexes
+ * in the given array.
+ */
+ private static Point2D.Double computeCentroid(ArrayList<ModeleBase> rnaBases, int[] indexes) {
+ double centroidX = 0, centroidY = 0;
+ for (int i=0; i<indexes.length; i++) {
+ int index = indexes[i];
+ Point2D coords = rnaBases.get(index).getCoords();
+ centroidX += coords.getX();
+ centroidY += coords.getY();
+ }
+ centroidX /= indexes.length;
+ centroidY /= indexes.length;
+
+ return new Point2D.Double(centroidX, centroidY);
+ }
+
+ /**
+ * The purpose of this class is to compute the rotation that minimizes the
+ * RMSD between the bases of the first RNA and the bases of the second RNA.
+ *
+ * @author Raphael Champeimont
+ */
+ private static class MinimizeRMSD {
+ private double[] X1, X2, Y1, Y2;
+
+ public MinimizeRMSD(double[] X1, double[] Y1, double[] X2, double[] Y2) {
+ this.X1 = X1;
+ this.X2 = X2;
+ this.Y1 = Y1;
+ this.Y2 = Y2;
+ }
+
+ /**
+ * A function such that minimizing it is equivalent to
+ * minimize the RMSD between between the two arrays of vectors,
+ * supposing we rotate the points in [X2,Y2] by theta.
+ */
+ private double f(double theta) {
+ double cos_theta = Math.cos(theta);
+ double sin_theta = Math.sin(theta);
+ int n = X1.length;
+ double sum = 0;
+ for (int i=0; i<n; i++) {
+ double dsx = X2[i]*cos_theta - Y2[i]*sin_theta - X1[i];
+ double dsy = X2[i]*sin_theta + Y2[i]*cos_theta - Y1[i];
+ sum = sum + dsx*dsx + dsy*dsy;
+ }
+ return sum;
+ }
+
+ /**
+ * d/dtheta f(theta)
+ */
+ private double fprime(double theta) {
+ double cos_theta = Math.cos(theta);
+ double sin_theta = Math.sin(theta);
+ int n = X1.length;
+ double sum = 0;
+ for (int i=0; i<n; i++) {
+ sum = sum
+ + (X1[i]*X2[i] + Y1[i]*Y2[i]) * sin_theta
+ + (X1[i]*Y2[i] - X2[i]*Y1[i]) * cos_theta;
+ }
+ return sum;
+ }
+
+ /**
+ * d^2/dtheta^2 f(theta)
+ */
+ private double fsecond(double theta) {
+ double cos_theta = Math.cos(theta);
+ double sin_theta = Math.sin(theta);
+ int n = X1.length;
+ double sum = 0;
+ for (int i=0; i<n; i++) {
+ sum = sum
+ + (X1[i]*X2[i] + Y1[i]*Y2[i]) * cos_theta
+ + (X2[i]*Y1[i] - X1[i]*Y2[i]) * sin_theta;
+ }
+ return sum;
+ }
+
+
+ /**
+ * Find the theta that minimizes f(theta).
+ * We use Newton's method.
+ */
+ public double computeOptimalTheta() {
+ final double epsilon = 1E-5;
+ double x_n = 0;
+ int numsteps = 0;
+ // In practice the number of steps to reach precision 1E-5 is smaller that
+ // 10 most of the time.
+ final int maxnumsteps = 100;
+ double x_n_plus_1;
+ double result;
+ while (true) {
+ numsteps = numsteps + 1;
+ double d = fsecond(x_n);
+ if (d == 0) {
+ // if f''(x_n) is 0 we cannot divide by it,
+ // so we move a little.
+ x_n_plus_1 = x_n + epsilon * (Math.random() - 0.5);
+ } else {
+ x_n_plus_1 = x_n - fprime(x_n)/fsecond(x_n);
+ if (Math.abs(x_n_plus_1 - x_n) < epsilon) {
+ result = x_n_plus_1;
+ break;
+ }
+ }
+ if (numsteps >= maxnumsteps) {
+ // If things go bad (for example f''(x) keeps being 0)
+ // we need to give up after what we consider to be too many steps.
+ // In practice this can happen only in pathological cases
+ // like f being constant, which is very unlikely.
+ result = x_n_plus_1;
+ break;
+ }
+ x_n = x_n_plus_1;
+ }
+
+ // We now have either found the min or the max at x = result.
+ // If we have the max at x we know the min is at x+pi.
+ if (f(result + Math.PI) < f(result)) {
+ result = result + Math.PI;
+ }
+
+ return result;
+ }
+
+
+ }
+
+
+ /**
+ * We suppose we have two lists of points. The coordinates of the first
+ * list of points are X1 and Y1, and X2 and Y2 for the second list.
+ * We suppose that the centroids of [X1,Y1] and [X2,Y2] are both (0,0).
+ * This function computes the rotation to apply to the second set of
+ * points such that the RMSD (Root mean square deviation) between them
+ * is minimum. The returned angle is in radians.
+ */
+ private static double minimizeRotateRMSD(double[] X1, double[] Y1, double[] X2, double[] Y2) {
+ MinimizeRMSD minimizer = new MinimizeRMSD(X1, Y1, X2, Y2);
+ return minimizer.computeOptimalTheta();
+ }
+
+
+ /**
+ * Move rna2 using a rotation so that it would move as little as possible
+ * (in the sense that we minimize the RMSD) when transforming
+ * it to rna1 using the given mapping.
+ */
+ public static void moveNearOtherRNA(RNA rna1, RNA rna2, Mapping mapping) {
+ int[] rna1MappedElems = mapping.getSourceElems();
+ int[] rna2MappedElems = mapping.getTargetElems();
+ ArrayList<ModeleBase> rna1Bases = rna1.get_listeBases();
+ ArrayList<ModeleBase> rna2Bases = rna2.get_listeBases();
+ int n = rna1MappedElems.length;
+
+ // If there is less than 2 points, it is useless to rotate the RNA.
+ if (n < 2) return;
+
+ // We can now assume that n >= 2.
+
+ // Compute the centroids of both RNAs
+ Point2D.Double rna1MappedElemsCentroid = computeCentroid(rna1Bases, rna1MappedElems);
+ Point2D.Double rna2MappedElemsCentroid = computeCentroid(rna2Bases, rna2MappedElems);
+
+ // Compute the optimal rotation
+ // We first compute coordinates for both RNAs, changing the origins
+ // to be the centroids.
+ double[] X1 = new double[rna1MappedElems.length];
+ double[] Y1 = new double[rna1MappedElems.length];
+ double[] X2 = new double[rna2MappedElems.length];
+ double[] Y2 = new double[rna2MappedElems.length];
+ for (int i=0; i<rna1MappedElems.length; i++) {
+ int base1Index = rna1MappedElems[i];
+ Point2D.Double coords1 = rna1Bases.get(base1Index).getCoords();
+ X1[i] = coords1.x - rna1MappedElemsCentroid.x;
+ Y1[i] = coords1.y - rna1MappedElemsCentroid.y;
+ Point2D.Double coords2 = rna2Bases.get(mapping.getPartner(base1Index)).getCoords();
+ X2[i] = coords2.x - rna2MappedElemsCentroid.x;
+ Y2[i] = coords2.y - rna2MappedElemsCentroid.y;
+ }
+
+ // Compute the optimal rotation angle theta
+ double theta = minimizeRotateRMSD(X1, Y1, X2, Y2);
+
+ // Rotate RNA2
+ rna2.globalRotation(theta * 180.0 / Math.PI);
+ }
+
+ public void nextTarget(RNA _target, VARNAConfig _conf, Mapping _mapping) {
+ nextTarget(_target, _conf, _mapping, false);
+ }
+
+ Runnable _loop, _end;
+
+ /**
+ * The argument moveTarget specifies whether the RNA _target should
+ * be rotated so that bases move as little as possible when switching
+ * from the current RNA to _target using the animation.
+ * Note that this will modify the _target object directly.
+ */
+ public void nextTarget(final RNA _target, final VARNAConfig _conf, Mapping _mapping, boolean moveTarget)
+ {
+ _end = new Runnable() {
+
+ @Override
+ public void run() {
+ _vpn.showRNA(_target);
+ _vpn.repaint();
+ }
+
+ };
+
+ try {
+ final RNA source = _vpn.getRNA();
+
+ if (moveTarget) moveNearOtherRNA(source, _target, _mapping);
+
+ if (source.getSize()!=0&&_target.getSize()!=0)
+ {
+ ArrayList<ModeleBase> currBases = source.get_listeBases();
+ ArrayList<ModeleBase> destBases = _target.get_listeBases();
+ Vector<Vector<Integer>> intArrSource = new Vector<Vector<Integer>>();
+ Vector<Vector<Integer>> intArrTarget = new Vector<Vector<Integer>>();
+ // Building interval arrays
+ intArrSource = clusterIndices(currBases.size(), _mapping
+ .getSourceElems());
+ intArrTarget = clusterIndices(destBases.size(), _mapping
+ .getTargetElems());
+
+ // Duplicating source and target coordinates
+ final Point2D.Double[] initPosSource = new Point2D.Double[currBases
+ .size()];
+ final Point2D.Double[] finalPosTarget = new Point2D.Double[destBases
+ .size()];
+
+ for (int i = 0; i < currBases.size(); i++) {
+ Point2D tmp = currBases.get(i).getCoords();
+ initPosSource[i] = new Point2D.Double(tmp.getX(), tmp.getY());
+ }
+ for (int i = 0; i < destBases.size(); i++) {
+ Point2D tmp = destBases.get(i).getCoords();
+ finalPosTarget[i] = new Point2D.Double(tmp.getX(), tmp.getY());
+ }
+
+ /**
+ * Assigning final (Source) and initial (Target) coordinates
+ */
+ final Point2D.Double[] finalPosSource = new Point2D.Double[initPosSource.length];
+ final Point2D.Double[] initPosTarget = new Point2D.Double[finalPosTarget.length];
+ // Final position of source model
+ for (int i = 0; i < finalPosSource.length; i++) {
+ if (_mapping.getPartner(i) != Mapping.UNKNOWN) {
+ Point2D dest;
+ dest = finalPosTarget[_mapping.getPartner(i)];
+ finalPosSource[i] = new Point2D.Double(dest.getX(), dest
+ .getY());
+ }
+ }
+
+ for (int i = 0; i < intArrSource.size(); i += 2) {
+ int matchedNeighborLeft, matchedNeighborRight;
+ if (i == 0) {
+ matchedNeighborLeft = intArrSource.get(1).get(0);
+ matchedNeighborRight = intArrSource.get(1).get(0);
+ } else if (i == intArrSource.size() - 1) {
+ matchedNeighborLeft = intArrSource.get(
+ intArrSource.size() - 2).get(0);
+ matchedNeighborRight = intArrSource.get(
+ intArrSource.size() - 2).get(0);
+ } else {
+ matchedNeighborLeft = intArrSource.get(i - 1).get(0);
+ matchedNeighborRight = intArrSource.get(i + 1).get(0);
+ }
+ Vector<Integer> v = intArrSource.get(i);
+ for (int j = 0; j < v.size(); j++) {
+ int index = v.get(j);
+ finalPosSource[index] = computeDestination(
+ initPosSource[matchedNeighborLeft],
+ initPosSource[matchedNeighborRight],
+ initPosSource[index], j + 1, v.size() + 1,
+ finalPosSource[matchedNeighborLeft],
+ finalPosSource[matchedNeighborRight]);
+ }
+ }
+ for (int i = 0; i < initPosTarget.length; i++) {
+ if (_mapping.getAncestor(i) != Mapping.UNKNOWN) {
+ Point2D dest;
+ dest = initPosSource[_mapping.getAncestor(i)];
+ initPosTarget[i] = new Point2D.Double(dest.getX(), dest.getY());
+ }
+ }
+ for (int i = 0; i < intArrTarget.size(); i += 2) {
+ int matchedNeighborLeft, matchedNeighborRight;
+ if (i == 0) {
+ matchedNeighborLeft = intArrTarget.get(1).get(0);
+ matchedNeighborRight = intArrTarget.get(1).get(0);
+ } else if (i == intArrTarget.size() - 1) {
+ matchedNeighborLeft = intArrTarget.get(
+ intArrTarget.size() - 2).get(0);
+ matchedNeighborRight = intArrTarget.get(
+ intArrTarget.size() - 2).get(0);
+ } else {
+ matchedNeighborLeft = intArrTarget.get(i - 1).get(0);
+ matchedNeighborRight = intArrTarget.get(i + 1).get(0);
+ }
+ Vector<Integer> v = intArrTarget.get(i);
+ for (int j = 0; j < v.size(); j++) {
+ int index = v.get(j);
+ initPosTarget[index] = computeDestination(
+ finalPosTarget[matchedNeighborLeft],
+ finalPosTarget[matchedNeighborRight],
+ finalPosTarget[index], j + 1, v.size() + 1,
+ initPosTarget[matchedNeighborLeft],
+ initPosTarget[matchedNeighborRight]);
+ }
+ }
+
+ mode = START;
+ _loop = new Runnable(){
+
+ @Override
+ public void run() {
+ int i = _step;
+ RNA current = (_firstHalf ? source : _target);
+ if (i == _numSteps / 2) {
+ _vpn.showRNA(_target);
+ current = _target;
+ _firstHalf = false;
+ if (_conf!=null)
+ {_vpn.setConfig(_conf);}
+
+ for (int j = 0; j < initPosSource.length; j++) {
+ source.setCoord(j, initPosSource[j]);
+ }
+ }
+ ArrayList<ModeleBase> currBases = current.get_listeBases();
+ for (int j = 0; j < currBases.size(); j++) {
+ ModeleBase m = currBases.get(j);
+ Point2D mpc, mnc;
+ if (_firstHalf) {
+ mpc = initPosSource[j];
+ mnc = finalPosSource[j];
+ } else {
+ mpc = initPosTarget[j];
+ mnc = finalPosTarget[j];
+ }
+ m.setCoords(new Point2D.Double(((_numSteps - 1 - i)
+ * mpc.getX() + (i) * mnc.getX())
+ / (_numSteps - 1), ((_numSteps - 1 - i)
+ * mpc.getY() + i * mnc.getY())
+ / (_numSteps - 1)));
+ }
+ _vpn.repaint();
+ }
+
+ };
+ actionPerformed(null);
+ } else {
+ _end.run();
+ }
+ } catch (MappingException e) {
+ e.printStackTrace();
+ _end.run();
+ }catch (Exception e) {
+ e.printStackTrace();
+ _end.run();
+ }
+
+
+ }
+
+ private class TargetsHolder
+ {
+ public RNA target;
+ public VARNAConfig conf;
+ public Mapping mapping;
+ public TargetsHolder(RNA t, VARNAConfig c, Mapping m)
+ {
+ target = t;
+ conf = c;
+ mapping = m;
+ }
+ }
+
+ private class Targets
+ {
+ LinkedList<TargetsHolder> _d = new LinkedList<TargetsHolder>();
+ public Targets() {
+ // BH j2s SwingJS added only to remove Eclipse warning
+ }
+ public synchronized void add(TargetsHolder d)
+ {
+ _d.addLast(d);
+
+ @SuppressWarnings("unused")
+ Runnable interpolator = ControleurInterpolator.this;
+ /**
+ * BH SwingJS no notify()
+ * @j2sNative
+ *
+ * interpolator.run$();
+ *
+ */
+ {
+ notify();
+ }
+ }
+
+ public synchronized TargetsHolder get() {
+
+ /**
+ * BH SwingJS no wait()
+ *
+ * @j2sNative
+ *
+ *
+ */
+ {
+ while (_d.size() == 0) {
+ try {
+ wait();
+ } catch (InterruptedException e) {
+ e.printStackTrace();
+ }
+ }
+ }
+ TargetsHolder x = _d.getLast();
+ _d.clear();
+ return x;
+ }
+ }
+
+
+ @Override
+ public void actionPerformed(ActionEvent e) {
+ runAnimation();
+ }
+
+ private int mode;
+
+ private void runAnimation() {
+ switch (mode) {
+ case START:
+ _firstHalf = true;
+ _step = 0;
+ mode = LOOP;
+ // Fall through
+ case LOOP:
+ if (_step < _numSteps) {
+ _loop.run();
+ ++_step;
+ break;
+ }
+ mode = END;
+ // Fall through
+ case END:
+ _end.run();
+ return;
+ }
+ Timer t = new Timer((int) _timeDelay, this);
+ t.setRepeats(false);
+ t.start();
+ // try {
+ // for (int i = 0; i < _numSteps; i++) {
+ // _step = i;
+ // loop.run();
+ //
+ // sleep(_timeDelay);
+ // }
+ // end.run();
+ // } catch (InterruptedException e) {
+ // e.printStackTrace();
+ // end.run();
+ // }
+ }
+
+}