X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=srcjar%2Ffr%2Forsay%2Flri%2Fvarna%2Fcontrolers%2FControleurInterpolator.java;fp=srcjar%2Ffr%2Forsay%2Flri%2Fvarna%2Fcontrolers%2FControleurInterpolator.java;h=a77bac6e5be581e195411a1059d6c33b812931c1;hb=4898f0ae429e0c61ddba72ca46be89b34bb4df8b;hp=0000000000000000000000000000000000000000;hpb=5a6ac5b535856903629234ad43a71319a91ebee5;p=jalview.git

diff --git a/srcjar/fr/orsay/lri/varna/controlers/ControleurInterpolator.java b/srcjar/fr/orsay/lri/varna/controlers/ControleurInterpolator.java
new file mode 100644
index 0000000..a77bac6
--- /dev/null
+++ b/srcjar/fr/orsay/lri/varna/controlers/ControleurInterpolator.java
@@ -0,0 +1,660 @@
+/*
+ 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();
+		// }
+	}
+
+}