/*
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.controlers;
import java.awt.geom.Point2D;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.LinkedList;
import java.util.Vector;
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 {
VARNAPanel _vpn;
private int _numSteps = 25;
private long _timeDelay = 15;
private boolean _running = false;
Targets _d = new Targets();
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;
}
}
/**
* 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);
}
/**
* 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(RNA _target, VARNAConfig _conf, Mapping _mapping, boolean moveTarget)
{
try {
RNA source = _vpn.getRNA();
RNA current = source;
if (moveTarget) moveNearOtherRNA(source, _target, _mapping);
ArrayList<ModeleBase> currBases = current.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
Point2D.Double[] initPosSource = new Point2D.Double[currBases
.size()];
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
*/
Point2D.Double[] finalPosSource = new Point2D.Double[initPosSource.length];
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]);
}
}
boolean firstHalf = true;
for (int i = 0; i < _numSteps; i++) {
if (i == _numSteps / 2) {
_vpn.showRNA(_target);
current = _target;
currBases = current.get_listeBases();
firstHalf = false;
if (_conf!=null)
{_vpn.setConfig(_conf);}
for (int j = 0; j < initPosSource.length; j++) {
source.setCoord(j, initPosSource[j]);
}
}
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();
sleep(_timeDelay);
}
} catch (InterruptedException e) {
e.printStackTrace();
} catch (MappingException e) {
e.printStackTrace();
}catch (Exception e) {
e.printStackTrace();
}
_vpn.showRNA(_target);
_vpn.repaint();
}
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 synchronized void add(TargetsHolder d)
{
_d.addLast(d);
notify();
}
public synchronized TargetsHolder get()
{
while(_d.size()==0)
{
try {
wait();
} catch (InterruptedException e) {
e.printStackTrace();
}
}
TargetsHolder x = _d.getFirst();
_d.removeFirst();
return x;
}
}
}