/**
 * File written by Raphael Champeimont
 * UMR 7238 Genomique des Microorganismes
 */
package fr.orsay.lri.varna.models.rna;

import java.awt.geom.AffineTransform;
import java.awt.geom.Line2D;
import java.awt.geom.Point2D;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.HashMap;
import java.util.HashSet;
import java.util.Iterator;
import java.util.Map;
import java.util.Set;

import fr.orsay.lri.varna.exceptions.ExceptionInvalidRNATemplate;
import fr.orsay.lri.varna.models.VARNAConfig;
import fr.orsay.lri.varna.models.geom.ComputeArcCenter;
import fr.orsay.lri.varna.models.geom.ComputeEllipseAxis;
import fr.orsay.lri.varna.models.geom.CubicBezierCurve;
import fr.orsay.lri.varna.models.geom.HalfEllipse;
import fr.orsay.lri.varna.models.geom.LinesIntersect;
import fr.orsay.lri.varna.models.geom.MiscGeom;
import fr.orsay.lri.varna.models.templates.DrawRNATemplateCurveMethod;
import fr.orsay.lri.varna.models.templates.DrawRNATemplateMethod;
import fr.orsay.lri.varna.models.templates.RNANodeValueTemplate;
import fr.orsay.lri.varna.models.templates.RNANodeValueTemplateBasePair;
import fr.orsay.lri.varna.models.templates.RNATemplate;
import fr.orsay.lri.varna.models.templates.RNATemplateAlign;
import fr.orsay.lri.varna.models.templates.RNATemplateDrawingAlgorithmException;
import fr.orsay.lri.varna.models.templates.RNATemplateMapping;
import fr.orsay.lri.varna.models.templates.RNATemplate.EdgeEndPointPosition;
import fr.orsay.lri.varna.models.templates.RNATemplate.In1Is;
import fr.orsay.lri.varna.models.templates.RNATemplate.RNATemplateElement;
import fr.orsay.lri.varna.models.templates.RNATemplate.RNATemplateHelix;
import fr.orsay.lri.varna.models.templates.RNATemplate.RNATemplateUnpairedSequence;
import fr.orsay.lri.varna.models.templates.RNATemplate.RNATemplateElement.EdgeEndPoint;
import fr.orsay.lri.varna.models.treealign.Tree;

public class DrawRNATemplate {
	private RNA rna;
	private RNATemplateMapping mapping;
	private ArrayList<ModeleBase> _listeBases;
	
	
	public DrawRNATemplate(RNA rna) {
		this.rna = rna;
		this._listeBases = rna.getListeBases();
	}
	
	public RNATemplateMapping getMapping() {
		return mapping;
	}
	
	
	
	/**
	 * Draw this RNA like the given template.
	 * The helixLengthAdjustmentMethod argument tells what to do in case
	 * some helices are of a different length in the template and the
	 * actual helix. See class DrawRNATemplateMethod above for possible values.
	 * @param straightBulges 
	 */
	public void drawRNATemplate(
			RNATemplate template,
			VARNAConfig conf,
			DrawRNATemplateMethod helixLengthAdjustmentMethod,
			DrawRNATemplateCurveMethod curveMethod, boolean straightBulges)
	throws RNATemplateDrawingAlgorithmException {
		
		// debug
//		try {
//			RNA perfectMatchingRNA = template.toRNA();
//			System.out.println("An RNA that would perfectly match this template would be:");
//			System.out.println(perfectMatchingRNA.getStructDBN());
//		} catch (ExceptionInvalidRNATemplate e) {
//			e.printStackTrace();
//		}
		
		mapping = RNATemplateAlign.mapRNAWithTemplate(rna, template);
		//System.out.println(mapping.showCompact(this));
		
		// debug
//			RNATemplateAlign.printMapping(mapping, template, getSeq());
//			try {
//				TreeGraphviz.treeToGraphvizPostscript(alignment, "alignment_graphviz.ps");
//			} catch (IOException e) {
//				e.printStackTrace();
//			}
		
		

		Iterator<RNATemplateElement> iter;
		double globalIncreaseFactor = 1;
		Map<RNATemplateHelix, Point2D.Double> translateVectors = null;
		if (helixLengthAdjustmentMethod == DrawRNATemplateMethod.MAXSCALINGFACTOR) {
			// Compute increase factors for helices.
			Map<RNATemplateHelix, Double> lengthIncreaseFactor = new HashMap<RNATemplateHelix, Double>();
			double maxLengthIncreaseFactor = Double.NEGATIVE_INFINITY;
			//RNATemplateHelix maxIncreaseHelix = null;
			iter = template.rnaIterator();
			while (iter.hasNext()) {
				RNATemplateElement element = iter.next();
				if (element instanceof RNATemplateHelix
						&& mapping.getAncestor(element) != null
						&& !lengthIncreaseFactor.containsKey(element)) {
					RNATemplateHelix helix = (RNATemplateHelix) element;
					int[] basesInHelixArray = RNATemplateAlign.intArrayFromList(mapping.getAncestor(helix));
					Arrays.sort(basesInHelixArray);
					double l = computeLengthIncreaseFactor(basesInHelixArray, helix, straightBulges);
					lengthIncreaseFactor.put(helix, l);
					if (l > maxLengthIncreaseFactor) {
						maxLengthIncreaseFactor = l;
						//maxIncreaseHelix = helix;
					}
				}
			}
			
			// debug
			//System.out.println("Max helix length increase factor = " + maxLengthIncreaseFactor + " reached with helix " + maxIncreaseHelix);;
			
			globalIncreaseFactor = Math.max(1, maxLengthIncreaseFactor);
			
		} else if (helixLengthAdjustmentMethod == DrawRNATemplateMethod.HELIXTRANSLATE) {
			try {
				// Now we need to propagate this helices translations
				Tree<RNANodeValueTemplate> templateAsTree = template.toTree();
				translateVectors = computeHelixTranslations(templateAsTree, mapping, straightBulges);
			
			} catch (ExceptionInvalidRNATemplate e) {
				throw (new RNATemplateDrawingAlgorithmException("ExceptionInvalidRNATemplate: " + e.getMessage()));
			}
		}
		
		// Allocate the coords and centers arrays
		// We create Point2D.Double objects in it but the algorithms
		// we use may choose to create new Point2D.Double objects or to
		// modify those created here.
		Point2D.Double[] coords = new Point2D.Double[_listeBases.size()];
		Point2D.Double[] centers = new Point2D.Double[_listeBases.size()];
		double[] angles = new double[_listeBases.size()];
		for (int i = 0; i < _listeBases.size(); i++) {
			coords[i] = new Point2D.Double(0, 0);
			centers[i] = new Point2D.Double(0, 0);
		}
		
		boolean computeCoords = true;
		while (computeCoords) {
			computeCoords = false;
			// Compute coords and centers
			Set<RNATemplateHelix> alreadyDrawnHelixes = new HashSet<RNATemplateHelix>();
			RNATemplateHelix lastMappedHelix = null;
			EdgeEndPoint howWeGotOutOfLastHelix = null;
			int howWeGotOutOfLastHelixBaseIndex = 0;
			iter = template.rnaIterator();
			RNATemplateElement element = null;
			while (iter.hasNext()) {
				element = iter.next();
				if (element instanceof RNATemplateHelix
						&& mapping.getAncestor(element) != null) {
					// We have a mapping between an helix in the RNA sequence
					// and an helix in the template.
					
					RNATemplateHelix helix = (RNATemplateHelix) element;
					boolean firstTimeWeMeetThisHelix;
					int[] basesInHelixArray = RNATemplateAlign.intArrayFromList(mapping.getAncestor(helix));
					Arrays.sort(basesInHelixArray);
					
					// Draw this helix if it has not already been done
					if (!alreadyDrawnHelixes.contains(helix)) {
						firstTimeWeMeetThisHelix = true;
						drawHelixLikeTemplateHelix(basesInHelixArray, helix, coords, centers,angles, globalIncreaseFactor, translateVectors, straightBulges);
						alreadyDrawnHelixes.add(helix);
					} else {
						firstTimeWeMeetThisHelix = false;
					}
					
					EdgeEndPoint howWeGetInCurrentHelix;
					if (firstTimeWeMeetThisHelix) {
						if (helix.getIn1Is() == In1Is.IN1_IS_5PRIME) {
							howWeGetInCurrentHelix = helix.getIn1();
						} else {
							howWeGetInCurrentHelix = helix.getIn2();
						}
					} else {
						if (helix.getIn1Is() == In1Is.IN1_IS_5PRIME) {
							howWeGetInCurrentHelix = helix.getIn2();
						} else {
							howWeGetInCurrentHelix = helix.getIn1();
						}
					}
					
					Point2D.Double P0 = new Point2D.Double();
					Point2D.Double P3 = new Point2D.Double();
					
					if (lastMappedHelix != null) {
						// Now draw the RNA sequence (possibly containing helixes)
						// between the last template drawn helix and this one.
	
						if (lastMappedHelix == helix) {
							// Last helix is the same as the current one so
							// nothing matched (or at best a single
							// non-paired sequence) so we will just
							// use the Radiate algorithm
							
							Point2D.Double helixVector = new Point2D.Double();
							computeHelixEndPointDirections(howWeGotOutOfLastHelix, helixVector, new Point2D.Double());
							
							double angle = MiscGeom.angleFromVector(helixVector);
							int b1 = basesInHelixArray[basesInHelixArray.length/2 - 1];
							P0.setLocation(coords[b1]);
							int b2 = basesInHelixArray[basesInHelixArray.length/2];
							P3.setLocation(coords[b2]);
							Point2D.Double loopCenter = new Point2D.Double((P0.x + P3.x)/2, (P0.y + P3.y)/2);
							rna.drawLoop(b1,
									 b2,
									 loopCenter.x,
									 loopCenter.y,
									 angle,
									 coords,
									 centers,
									 angles,
									 straightBulges);
							// If the helix is flipped, we need to compute the symmetric
							// of the whole loop.
							if (helix.isFlipped()) {
								symmetric(loopCenter, helixVector, coords, b1, b2);
								symmetric(loopCenter, helixVector, centers, b1, b2);
							}
						} else {
							// No helices matched between the last helix and
							// the current one, so we draw what is between
							// using the radiate algorithm but on the Bezier curve.
							
							int b1 = howWeGotOutOfLastHelixBaseIndex;
							int b2 = firstTimeWeMeetThisHelix ? basesInHelixArray[0] : basesInHelixArray[basesInHelixArray.length/2];
							P0.setLocation(coords[b1]);
							P3.setLocation(coords[b2]);
							
							Point2D.Double P1, P2;
							
							if (howWeGotOutOfLastHelix.getOtherElement() instanceof RNATemplateUnpairedSequence
									&& howWeGetInCurrentHelix.getOtherElement() instanceof RNATemplateUnpairedSequence) {
								// We will draw the bases on a Bezier curve
								P1 = new Point2D.Double();
								computeBezierTangentVectorTarget(howWeGotOutOfLastHelix, P0, P1);
								
								P2 = new Point2D.Double();
								computeBezierTangentVectorTarget(howWeGetInCurrentHelix, P3, P2);
							} else {
								// We will draw the bases on a straight line between P0 and P3
								P1 = null;
								P2 = null;
							}
							
							drawAlongCurve(b1, b2, P0, P1, P2, P3, coords, centers, angles, curveMethod, lastMappedHelix.isFlipped(), straightBulges);
						}
						
					} else if (basesInHelixArray[0] > 0) {
						// Here we draw what is before the first mapped helix.
	
						RNATemplateUnpairedSequence templateSequence;
						// Try to find our template sequence as the mapped element of base 0
						RNATemplateElement templateSequenceCandidate = mapping.getPartner(0);
						if (templateSequenceCandidate != null
								&& templateSequenceCandidate instanceof RNATemplateUnpairedSequence) {
							templateSequence = (RNATemplateUnpairedSequence) templateSequenceCandidate;
						} else {
							// Try other idea: first template element if it is a sequence
							templateSequenceCandidate = template.getFirst();
							if (templateSequenceCandidate != null
									&& templateSequenceCandidate instanceof RNATemplateUnpairedSequence) {
								templateSequence = (RNATemplateUnpairedSequence) templateSequenceCandidate;
							} else {
								// We don't know where to start
								templateSequence = null;
							}
						}
						
						int b1 = 0;
						int b2 = firstTimeWeMeetThisHelix ? basesInHelixArray[0] : basesInHelixArray[basesInHelixArray.length/2];
						P3.setLocation(coords[b2]);
						
						if (templateSequence != null) {
							coords[0].setLocation(templateSequence.getVertex5());
							coords[0].x *= globalIncreaseFactor;
							coords[0].y *= globalIncreaseFactor;
						} else {
							// Put b1 at an ideal distance from b2, using the "flat exterior loop" method
							double idealLength = computeStraightLineIdealLength(b1, b2);
							Point2D.Double j = new Point2D.Double();
							if (howWeGetInCurrentHelix != null) {
								computeHelixEndPointDirections(howWeGetInCurrentHelix, new Point2D.Double(), j);
							} else {
								j.setLocation(1, 0);
							}
							coords[b1].setLocation(coords[b2].x + j.x*idealLength, coords[b2].y + j.y*idealLength);
						}
						P0.setLocation(coords[0]);
						
						Point2D.Double P1, P2;
						
						if (howWeGetInCurrentHelix.getOtherElement() instanceof RNATemplateUnpairedSequence
								&& templateSequence != null) {
							// We will draw the bases on a Bezier curve
							P1 = new Point2D.Double();
							computeBezierTangentVectorTarget(templateSequence.getIn(), P0, P1);
							
							P2 = new Point2D.Double();
							computeBezierTangentVectorTarget(howWeGetInCurrentHelix, P3, P2);
						} else {
							// We will draw the bases on a straight line between P0 and P3
							P1 = null;
							P2 = null;
						}
						
						drawAlongCurve(b1, b2, P0, P1, P2, P3, coords, centers, angles, curveMethod, false, straightBulges);
					}
					
					lastMappedHelix = helix;
					howWeGotOutOfLastHelix = howWeGetInCurrentHelix.getNextEndPoint();
					if (firstTimeWeMeetThisHelix) {
						howWeGotOutOfLastHelixBaseIndex = basesInHelixArray[basesInHelixArray.length/2-1];
					} else {
						howWeGotOutOfLastHelixBaseIndex = basesInHelixArray[basesInHelixArray.length-1];
					}
				}
			} // end template iteration
			
			
			// Now we need to draw what is after the last mapped helix.
			if (howWeGotOutOfLastHelixBaseIndex < coords.length-1
					&& element != null
					&& coords.length > 1) {
				
				RNATemplateUnpairedSequence beginTemplateSequence = null;
				if (lastMappedHelix == null) {
					// No helix at all matched between the template and RNA!
					// So the sequence we want to draw is the full RNA.
					
					// Try to find our template sequence as the mapped element of base 0
					RNATemplateElement templateSequenceCandidate = mapping.getPartner(0);
					if (templateSequenceCandidate != null
							&& templateSequenceCandidate instanceof RNATemplateUnpairedSequence) {
						beginTemplateSequence = (RNATemplateUnpairedSequence) templateSequenceCandidate;
					} else {
						// Try other idea: first template element if it is a sequence
						templateSequenceCandidate = template.getFirst();
						if (templateSequenceCandidate != null
								&& templateSequenceCandidate instanceof RNATemplateUnpairedSequence) {
							beginTemplateSequence = (RNATemplateUnpairedSequence) templateSequenceCandidate;
						} else {
							// We don't know where to start
							beginTemplateSequence = null;
						}
					}
					
					if (beginTemplateSequence != null) {
						coords[0].setLocation(beginTemplateSequence.getVertex5());
						coords[0].x *= globalIncreaseFactor;
						coords[0].y *= globalIncreaseFactor;
					}
					
				}
				
				RNATemplateUnpairedSequence endTemplateSequence;
				// Try to find our template sequence as the mapped element of last base
				RNATemplateElement templateSequenceCandidate = mapping.getPartner(coords.length-1);
				if (templateSequenceCandidate != null
						&& templateSequenceCandidate instanceof RNATemplateUnpairedSequence) {
					endTemplateSequence = (RNATemplateUnpairedSequence) templateSequenceCandidate;
				} else {
					// Try other idea: last template element if it is a sequence
					templateSequenceCandidate = element;
					if (templateSequenceCandidate != null
							&& templateSequenceCandidate instanceof RNATemplateUnpairedSequence) {
						endTemplateSequence = (RNATemplateUnpairedSequence) templateSequenceCandidate;
					} else {
						// We don't know where to end
						endTemplateSequence = null;
					}
				}
				
				int b1 = howWeGotOutOfLastHelixBaseIndex;
				int b2 = coords.length - 1;
				
				if (endTemplateSequence != null) {
					coords[b2].setLocation(endTemplateSequence.getVertex3());
					coords[b2].x *= globalIncreaseFactor;
					coords[b2].y *= globalIncreaseFactor;
				} else {
					// Put b2 at an ideal distance from b1, using the "flat exterior loop" method
					double idealLength = computeStraightLineIdealLength(b1, b2);
					Point2D.Double j = new Point2D.Double();
					if (howWeGotOutOfLastHelix != null) {
						computeHelixEndPointDirections(howWeGotOutOfLastHelix, new Point2D.Double(), j);
					} else {
						j.setLocation(1, 0);
					}
					coords[b2].setLocation(coords[b1].x + j.x*idealLength, coords[b1].y + j.y*idealLength);
				}

				
				Point2D.Double P0 = new Point2D.Double();
				Point2D.Double P3 = new Point2D.Double();
				
				P0.setLocation(coords[b1]);
				P3.setLocation(coords[b2]);
				
				Point2D.Double P1, P2;
				
				if (howWeGotOutOfLastHelix != null
						&& howWeGotOutOfLastHelix.getOtherElement() instanceof RNATemplateUnpairedSequence
						&& endTemplateSequence != null) {
					// We will draw the bases on a Bezier curve
					P1 = new Point2D.Double();
					computeBezierTangentVectorTarget(howWeGotOutOfLastHelix, P0, P1);
					
					P2 = new Point2D.Double();
					computeBezierTangentVectorTarget(endTemplateSequence.getOut(), P3, P2);
				} else if (lastMappedHelix == null
						&& beginTemplateSequence != null
						&& endTemplateSequence != null) {
					// We will draw the bases on a Bezier curve
					P1 = new Point2D.Double();
					computeBezierTangentVectorTarget(beginTemplateSequence.getIn(), P0, P1);
					
					P2 = new Point2D.Double();
					computeBezierTangentVectorTarget(endTemplateSequence.getOut(), P3, P2);
				} else {
					// We will draw the bases on a straight line between P0 and P3
					P1 = null;
					P2 = null;
				}
				
				drawAlongCurve(b1, b2, P0, P1, P2, P3, coords, centers, angles, curveMethod, lastMappedHelix != null ? lastMappedHelix.isFlipped() : false, straightBulges);
			
			}
			
			
			if (helixLengthAdjustmentMethod == DrawRNATemplateMethod.NOINTERSECT && coords.length > 3) {
				// Are we happy with this value of globalIncreaseFactor?
				Line2D.Double[] lines = new Line2D.Double[coords.length-1];
				for (int i=0; i<coords.length-1; i++) {
					lines[i] = new Line2D.Double(coords[i], coords[i+1]);
				}
				int intersectLines = 0;
				for (int i=0; i<lines.length; i++) {
					for (int j=i+2; j<lines.length; j++) {
						if (LinesIntersect.linesIntersect(lines[i], lines[j])) {
							intersectLines++;
						}
					}
				}
				// If no intersection we keep this globalIncreaseFactor value
				if (intersectLines > 0) {
					// Don't increase more than a maximum value
					if (globalIncreaseFactor < 3) {
						globalIncreaseFactor += 0.1;
						//System.out.println("globalIncreaseFactor increased to " + globalIncreaseFactor);
						// Compute the drawing again
						computeCoords = true;
					}
				}
			}
			
		}
		
		// debug
		if (helixLengthAdjustmentMethod == DrawRNATemplateMethod.MAXSCALINGFACTOR
				|| helixLengthAdjustmentMethod == DrawRNATemplateMethod.NOINTERSECT) {
			//System.out.println("globalIncreaseFactor = " + globalIncreaseFactor);
		}
		
		// Now we actually move the bases, according to arrays coords and centers
		// and taking in account the space between bases parameter.
		for (int i = 0; i < _listeBases.size(); i++) {
			_listeBases.get(i).setCoords(
					new Point2D.Double(coords[i].x * conf._spaceBetweenBases,
							coords[i].y * conf._spaceBetweenBases));
			_listeBases.get(i).setCenter(
					new Point2D.Double(centers[i].x * conf._spaceBetweenBases,
							centers[i].y * conf._spaceBetweenBases));
		}
	
	}
	
	
	
	
	
	
	/**
	 * IN: Argument helixEndPoint is an IN argument (will be read),
	 * and must contain an helix edge endpoint.
	 * 
	 * The other arguments are OUT arguments
	 * (must be existing objects, content will be overwritten).
	 * 
	 * OUT: The i argument will contain a vector colinear to the vector
	 * from the helix startPosition to endPosition or the opposite
	 * depending on there the endpoint is (the endpoint will be on the
	 * destination side of the vector). ||i|| = 1
	 * 
	 * OUT: The j vector will contain an vector that is colinear
	 * to the last/first base pair connection on the side of this endpoint.
	 * The vector will be oriented to the side of the given endpoint.
	 * ||j|| = 1
	 */
	private void computeHelixEndPointDirections(
			EdgeEndPoint helixEndPoint, // IN
			Point2D.Double i, // OUT
			Point2D.Double j // OUT
			) {
		RNATemplateHelix helix = (RNATemplateHelix) helixEndPoint.getElement();
		Point2D.Double startpos = helix.getStartPosition();
		Point2D.Double endpos = helix.getEndPosition();
		Point2D.Double helixVector = new Point2D.Double();
		switch (helixEndPoint.getPosition()) {
		case IN1:
		case OUT2:
			helixVector.x = startpos.x - endpos.x;
			helixVector.y = startpos.y - endpos.y;
			break;
		case IN2:
		case OUT1:
			helixVector.x = endpos.x - startpos.x;
			helixVector.y = endpos.y - startpos.y;
			break;
		}
		double helixVectorLength = Math.hypot(helixVector.x, helixVector.y);
		// i is the vector which is colinear to helixVector and such that ||i|| = 1
		i.x = helixVector.x / helixVectorLength;
		i.y = helixVector.y / helixVectorLength;
		// Find j such that it is orthogonal to i, ||j|| = 1
		// and j goes to the side where the sequence will be connected
		switch (helixEndPoint.getPosition()) {
		case IN1:
		case IN2:
			// rotation of +pi/2
			j.x = - i.y;
			j.y =   i.x;
			break;
		case OUT1:
		case OUT2:
			// rotation of -pi/2
			j.x =   i.y;
			j.y = - i.x;
			break;
		}
		if (helix.isFlipped()) {
			j.x = - j.x;
			j.y = - j.y;
		}

	}
	
	/**
	 * A cubic Bezier curve can be defined by 4 points,
	 * see http://en.wikipedia.org/wiki/Bezier_curve#Cubic_B.C3.A9zier_curves
	 * For each of the curve end points, there is the last/first point of the
	 * curve and a point which gives the direction and length of the tangent
	 * vector on that side. This two points are respectively curveEndPoint
	 * and curveVectorOtherPoint.
	 * IN:  Argument helixVector is the vector formed by the helix,
	 *      in the right direction for our sequence.
	 * IN:  Argument curveEndPoint is the position of the endpoint on the helix.
	 * OUT: Argument curveVectorOtherPoint must be allocated
	 *      and the values will be modified.
	 */
	private void computeBezierTangentVectorTarget(
			EdgeEndPoint endPoint,
			Point2D.Double curveEndPoint,
			Point2D.Double curveVectorOtherPoint)
			throws RNATemplateDrawingAlgorithmException {
		
		boolean sequenceEndPointIsIn;
		RNATemplateUnpairedSequence sequence;
		
		if (endPoint.getElement() instanceof RNATemplateHelix) {
			sequence = (RNATemplateUnpairedSequence) endPoint.getOtherElement();
			EdgeEndPointPosition endPointPositionOnHelix = endPoint.getPosition();
			switch (endPointPositionOnHelix) {
			case IN1:
			case IN2:
				sequenceEndPointIsIn = false;
				break;
			default:
				sequenceEndPointIsIn = true;
			}
			
			EdgeEndPoint endPointOnHelix =
				sequenceEndPointIsIn ?
						sequence.getIn().getOtherEndPoint() :
						sequence.getOut().getOtherEndPoint();
			if (endPointOnHelix == null) {
				throw (new RNATemplateDrawingAlgorithmException("Sequence is not connected to an helix."));
			}
		} else {
			// The endpoint is on an unpaired sequence.
			sequence = (RNATemplateUnpairedSequence) endPoint.getElement();
			if (endPoint == sequence.getIn()) {
				// endpoint is 5'
				sequenceEndPointIsIn = true;
			} else {
				sequenceEndPointIsIn = false;
			}
		}

		double l =
			sequenceEndPointIsIn ?
				sequence.getInTangentVectorLength() :
				sequence.getOutTangentVectorLength();
		
		// Compute the absolute angle our line makes to the helix
		double theta =
			sequenceEndPointIsIn ?
				sequence.getInTangentVectorAngle() :
				sequence.getOutTangentVectorAngle();
		
		// Compute v, the tangent vector of the Bezier curve
		Point2D.Double v = new Point2D.Double();
		v.x = l * Math.cos(theta);
		v.y = l * Math.sin(theta);
		curveVectorOtherPoint.x = curveEndPoint.x + v.x;
		curveVectorOtherPoint.y = curveEndPoint.y + v.y;
	}
	

	/**
	 * Compute (actual helix length / helix length in template).
	 * @param straightBulges 
	 */
	private double computeLengthIncreaseFactor(
			int[] basesInHelixArray,  // IN
			RNATemplateHelix helix,    // IN
			boolean straightBulges
			) {
		double templateLength = computeHelixTemplateLength(helix);
		double realLength = computeHelixRealLength(basesInHelixArray,straightBulges);
		return realLength / templateLength;
	}
	
	/**
	 * Compute (actual helix vector - helix vector in template).
	 */
	private Point2D.Double computeLengthIncreaseDelta(
			int[] basesInHelixArray,  // IN
			RNATemplateHelix helix,   // IN
			boolean straightBulges
			) {
		double templateLength = computeHelixTemplateLength(helix);
		double realLength = computeHelixRealLength(basesInHelixArray,straightBulges);
		Point2D.Double i = new Point2D.Double();
		computeTemplateHelixVectors(helix, null, i, null);
		return new Point2D.Double(i.x*(realLength-templateLength), i.y*(realLength-templateLength));
	}
	
	/**
	 * Compute helix interesting vectors from template helix.
	 * @param helix The template helix you want to compute the vectors from.
	 * @param o This point coordinates will be set the origin of the helix (or not if null),
	 *          ie. the point in the middle of the base pair with the two most extreme bases.
	 * @param i Will be set to the normalized helix vector. (nothing done if null)
	 * @param j Will be set to the normalized helix base pair vector (5' -> 3'). (nothing done if null)
	 */
	private void computeTemplateHelixVectors(
			RNATemplateHelix helix,  // IN
			Point2D.Double o,        // OUT
			Point2D.Double i,        // OUT
			Point2D.Double j         // OUT
			) {
		Point2D.Double startpos, endpos;
		if (helix.getIn1Is() == In1Is.IN1_IS_5PRIME) {
			startpos = helix.getStartPosition();
			endpos = helix.getEndPosition();
		} else {
			endpos = helix.getStartPosition();
			startpos = helix.getEndPosition();
		}
		if (o != null) {
			o.x = startpos.x;
			o.y = startpos.y;
		}
		if (i != null || j != null) {
			// (i_x,i_y) is the vector between two consecutive bases of the same side of an helix
			if (i == null)
				i = new Point2D.Double();
			i.x = (endpos.x - startpos.x);
			i.y = (endpos.y - startpos.y);
			double i_original_norm = Math.hypot(i.x, i.y);
			// change its norm to 1
			i.x = i.x / i_original_norm;
			i.y = i.y / i_original_norm;
			if (j != null) {
				j.x = - i.y;
				j.y =   i.x;
				if (helix.isFlipped()) {
					j.x = - j.x;
					j.y = - j.y;
				}
				double j_original_norm = Math.hypot(j.x, j.y);
				// change (j_x,j_y) so that its norm is 1
				j.x = j.x / j_original_norm;
				j.y = j.y / j_original_norm;
			}
		}
	}
	
	
	/**
	 * Estimate bulge arc length.
	 */
	private double estimateBulgeArcLength(int firstBase, int lastBase) {
		if (firstBase + 1 == lastBase)
			return RNA.LOOP_DISTANCE; // there is actually no bulge
		double len = 0.0;
		int k = firstBase;
		while (k < lastBase) {
			int l = _listeBases.get(k).getElementStructure();
			if (k < l && l < lastBase) {
				len += RNA.BASE_PAIR_DISTANCE;
				k = l;
			} else {
				len += RNA.LOOP_DISTANCE;
				k++;
			}
		}
		return len;
	}
	
	
	/**
	 * Estimate bulge width, the given first and last bases must be those in the helix.
	 */
	private double estimateBulgeWidth(int firstBase, int lastBase) {
		double len = estimateBulgeArcLength(firstBase, lastBase);
		return 2 * (len / Math.PI);
	}
	
	
	/**
	 * Get helix length in template.
	 */
	private double computeHelixTemplateLength(RNATemplateHelix helix) {
		return Math.hypot(helix.getStartPosition().x - helix.getEndPosition().x,
				helix.getStartPosition().y - helix.getEndPosition().y);
	}
	
	
	/**
	 * Compute helix actual length (as drawHelixLikeTemplateHelix() would draw it).
	 */
	private double computeHelixRealLength(int[] basesInHelixArray, boolean straightBulges) {
		return drawHelixLikeTemplateHelix(basesInHelixArray, null, null, null, null, 0, null,straightBulges);
	}
	
	
	/**
	 * Result type of countUnpairedLine(), see below.
	 */
	private static class UnpairedLineCounts {
		public int nBP, nLD, total;
	}
	/**
	 * If we are drawing an unpaired region that may contains helices,
	 * and we are drawing it on a line (curve or straight, doesn't matter),
	 * how many intervals should have a base-pair length (start of an helix)
	 * and how many should have a consecutive unpaired bases length?
	 * Returns an array with three elements:
	 * - answer to first question
	 * - answer to second question
	 * - sum of both, ie. total number of intervals
	 *   (this is NOT lastBase-firstBase because bases deep in helices do not count)
	 */
	private UnpairedLineCounts countUnpairedLine(int firstBase, int lastBase) {
		UnpairedLineCounts counts = new UnpairedLineCounts();
		int nBP = 0;
		int nLD = 0;
		{
			int b = firstBase;
			while (b < lastBase) {
				int l = _listeBases.get(b).getElementStructure();
				if (b < l && l < lastBase) {
					nBP++;
					b = l;
				} else {
					nLD++;
					b++;
				}
			}
		}
		counts.nBP = nBP;
		counts.nLD = nLD;
		counts.total = nBP + nLD;
		return counts;
	}
	
	
	/**
	 * Draw the given helix (given as a *SORTED* array of indexes)
	 * like defined in the given template helix.
	 * OUT: The bases positions are not changed in fact,
	 *      instead the coords and centers arrays are modified.
	 * IN:  The helix origin position is multiplied by scaleHelixOrigin
	 *      and translateVectors.get(helix) is added.
	 * RETURN VALUE:
	 *      The length of the drawn helix.
	 * @param straightBulges 
	 * 
	 */
	private double drawHelixLikeTemplateHelix(
			int[] basesInHelixArray,  // IN
			RNATemplateHelix helix,   // IN  (optional, ie. may be null)
			Point2D.Double[] coords,  // OUT (optional, ie. may be null)
			Point2D.Double[] centers, // OUT (optional, ie. may be null)
			double[] angles,  // OUT 
			double scaleHelixOrigin,  // IN
			Map<RNATemplateHelix, Point2D.Double> translateVectors // IN (optional, ie. may be null)
, boolean straightBulges
			) {
		int n = basesInHelixArray.length / 2;
		if (n == 0)
			return 0;
		 // Default values when not template helix is provided:
		Point2D.Double o = new Point2D.Double(0, 0);
		Point2D.Double i = new Point2D.Double(1, 0);
		Point2D.Double j = new Point2D.Double(0, 1);
		boolean flipped = false;
		if (helix != null) {
			computeTemplateHelixVectors(helix, o, i, j);
			flipped = helix.isFlipped();
		}
		Point2D.Double li = new Point2D.Double(i.x*RNA.LOOP_DISTANCE, i.y*RNA.LOOP_DISTANCE);
		// We want o to be the point where the first base (5' end) is
		o.x = (o.x - j.x * RNA.BASE_PAIR_DISTANCE / 2) * scaleHelixOrigin;
		o.y = (o.y - j.y * RNA.BASE_PAIR_DISTANCE / 2) * scaleHelixOrigin;
		if (translateVectors != null && translateVectors.containsKey(helix)) {
			Point2D.Double v = translateVectors.get(helix);
			o.x = o.x + v.x;
			o.y = o.y + v.y;
		}
		
		// We need this array so that we can store positions even if coords == null
		Point2D.Double[] helixBasesPositions = new Point2D.Double[basesInHelixArray.length];
		for (int k=0; k<helixBasesPositions.length; k++) {
			helixBasesPositions[k] = new Point2D.Double();	
		}
		Point2D.Double accDelta = new Point2D.Double(0, 0);
		for (int k=0; k<n; k++) {
			int kp = 2*n-k-1;
			Point2D.Double p1 = helixBasesPositions[k]; // we assign the point *reference*
			Point2D.Double p2 = helixBasesPositions[kp];
			// Do we have a bulge between previous base pair and this one?
			boolean bulge = k >= 1 && (basesInHelixArray[k] != basesInHelixArray[k-1] + 1
					                   || basesInHelixArray[kp+1] != basesInHelixArray[kp] + 1);
			if (k >= 1) {
				if (basesInHelixArray[k] < basesInHelixArray[k-1]
				    || basesInHelixArray[kp+1] < basesInHelixArray[kp]) {
					throw new Error("Internal bug: basesInHelixArray must be sorted");
				}
				if (bulge) {
					// Estimate a good distance (delta) between the previous base pair and this one
					double delta1 = estimateBulgeWidth(basesInHelixArray[k-1], basesInHelixArray[k]);
					double delta2 = estimateBulgeWidth(basesInHelixArray[kp], basesInHelixArray[kp+1]);
					// The biggest bulge defines the width
					double delta = Math.max(delta1, delta2);

					if (coords != null) {
						// Now, where do we put the bases that are part of the bulge?
						for (int side=0; side<2; side++) {
							Point2D.Double pstart = new Point2D.Double();
							Point2D.Double pend = new Point2D.Double();
							Point2D.Double bisectVect = new Point2D.Double();
							Point2D.Double is = new Point2D.Double();
							int firstBase, lastBase;
							double alphasign = flipped ? -1 : 1;
							if (side == 0) {
								firstBase = basesInHelixArray[k-1];
								lastBase  = basesInHelixArray[k];
								pstart.setLocation(o.x + accDelta.x,
								                   o.y + accDelta.y);
								pend.setLocation(o.x + accDelta.x + i.x*delta,
								                 o.y + accDelta.y + i.y*delta);
								bisectVect.setLocation(-j.x, -j.y);
								is.setLocation(i);
							} else {
								firstBase = basesInHelixArray[kp];
								lastBase  = basesInHelixArray[kp+1];
								pstart.setLocation(o.x + accDelta.x + i.x*delta + j.x*RNA.BASE_PAIR_DISTANCE,
						                           o.y + accDelta.y + i.y*delta + j.y*RNA.BASE_PAIR_DISTANCE);
								pend.setLocation(o.x + accDelta.x + j.x*RNA.BASE_PAIR_DISTANCE,
								                 o.y + accDelta.y + j.y*RNA.BASE_PAIR_DISTANCE);

								bisectVect.setLocation(j);
								is.setLocation(-i.x, -i.y);
							}
							double arclen = estimateBulgeArcLength(firstBase, lastBase);
							double centerOnBisect = ComputeArcCenter.computeArcCenter(delta, arclen);
							
							// Should we draw the base on an arc or simply use a line?
							if (centerOnBisect > -1000) {
								Point2D.Double center = new Point2D.Double(pstart.x + is.x*delta/2 + bisectVect.x*centerOnBisect,
								                                           pstart.y + is.y*delta/2 + bisectVect.y*centerOnBisect);
								int b = firstBase;
								double len = 0;
								double r = Math.hypot(pstart.x - center.x, pstart.y - center.y);
								double alpha0 = MiscGeom.angleFromVector(pstart.x - center.x, pstart.y - center.y);
								while (b < lastBase) {
									int l = _listeBases.get(b).getElementStructure();
									if (b < l && l < lastBase) {
										len += RNA.BASE_PAIR_DISTANCE;
										b = l;
									} else {
										len += RNA.LOOP_DISTANCE;
										b++;
									}
									if (b < lastBase) {
										coords[b].x = center.x + r*Math.cos(alpha0 + alphasign*len/r);
										coords[b].y = center.y + r*Math.sin(alpha0 + alphasign*len/r);
									}
								}
							} else {
								// Draw on a line
								
								// Distance between paired bases cannot be changed
								// (imposed by helix width) but distance between other
								// bases can be adjusted.
								UnpairedLineCounts counts = countUnpairedLine(firstBase, lastBase);
								double LD = Math.max((delta - counts.nBP*RNA.BASE_PAIR_DISTANCE) / (double) counts.nLD, 0);
								//System.out.println("nBP=" + nBP + " nLD=" + nLD);
								double len = 0;
								{
									int b = firstBase;
									while (b < lastBase) {
										int l = _listeBases.get(b).getElementStructure();
										if (b < l && l < lastBase) {
											len += RNA.BASE_PAIR_DISTANCE;
											b = l;
										} else {
											len += LD;
											b++;
										}
										if (b < lastBase) {
											coords[b].x = pstart.x + is.x*len;
											coords[b].y = pstart.y + is.y*len;
										}
									}
								}
								//System.out.println("len=" + len + " delta=" + delta + " d(pstart,pend)=" + Math.hypot(pend.x-pstart.x, pend.y-pstart.y));
							}
							
							// Does the bulge contain an helix?
							// If so, use drawLoop() to draw it.
							{
								int b = firstBase;
								while (b < lastBase) {
									int l = _listeBases.get(b).getElementStructure();
									if (b < l && l < lastBase) {
										// Helix present in bulge
										Point2D.Double b1pos = coords[b];
										Point2D.Double b2pos = coords[l];
										double beta = MiscGeom.angleFromVector(b2pos.x - b1pos.x, b2pos.y - b1pos.y) - Math.PI / 2 + (flipped ? Math.PI : 0);
										Point2D.Double loopCenter = new Point2D.Double((b1pos.x + b2pos.x)/2, (b1pos.y + b2pos.y)/2);
										rna.drawLoop(b,
												 l,
												 loopCenter.x,
												 loopCenter.y,
												 beta,
												 coords,
												 centers,
												 angles,
												 straightBulges);
										// If the helix is flipped, we need to compute the symmetric
										// of the whole loop.
										if (helix.isFlipped()) {
											Point2D.Double v = new Point2D.Double(Math.cos(beta), Math.sin(beta));
											symmetric(loopCenter, v, coords, b, l);
											symmetric(loopCenter, v, centers, b, l);
										}
										// Continue
										b = l;
									} else {
										b++;
									}
								}
							}
						}
					}
					
					accDelta.x += i.x*delta;
					accDelta.y += i.y*delta;
					p1.x = o.x + accDelta.x;
					p1.y = o.y + accDelta.y;
					p2.x = p1.x + j.x*RNA.BASE_PAIR_DISTANCE;
					p2.y = p1.y + j.y*RNA.BASE_PAIR_DISTANCE;
					
				} else {
					accDelta.x += li.x;
					accDelta.y += li.y;
					p1.x = o.x + accDelta.x;
					p1.y = o.y + accDelta.y;
					p2.x = p1.x + j.x*RNA.BASE_PAIR_DISTANCE;
					p2.y = p1.y + j.y*RNA.BASE_PAIR_DISTANCE;
				}
			} else {
				// First base pair
				p1.x = o.x;
				p1.y = o.y;
				p2.x = p1.x + j.x*RNA.BASE_PAIR_DISTANCE;
				p2.y = p1.y + j.y*RNA.BASE_PAIR_DISTANCE;
			}
		}
		
		Point2D.Double p1 = helixBasesPositions[0];
		Point2D.Double p2 = helixBasesPositions[n-1];
		
		if (coords != null) {
			for (int k=0; k<helixBasesPositions.length; k++) {
				coords[basesInHelixArray[k]] = helixBasesPositions[k];
			}
		}
		
		return Math.hypot(p2.x-p1.x, p2.y-p1.y);
	}
	
	
	

	/**
	 * Return ideal length to draw unpaired region from firstBase to lastBase.
	 */
	private double computeStraightLineIdealLength(int firstBase, int lastBase) {
		UnpairedLineCounts counts = countUnpairedLine(firstBase, lastBase);
		return RNA.BASE_PAIR_DISTANCE*counts.nBP + RNA.LOOP_DISTANCE*counts.nLD;
	}
	
	/**
	 * Like drawOnBezierCurve(), but on a straight line.
	 */
	private boolean drawOnStraightLine(
			int firstBase,
			int lastBase,
			Point2D.Double P0,
			Point2D.Double P3,
			Point2D.Double[] coords,
			Point2D.Double[] centers,
			boolean cancelIfNotGood) {
		
		Point2D.Double vector = new Point2D.Double(P3.x - P0.x, P3.y - P0.y); 
		double vectorNorm = Math.hypot(vector.x, vector.y);

		UnpairedLineCounts counts = countUnpairedLine(firstBase, lastBase);
		double LD = Math.max((vectorNorm - counts.nBP*RNA.BASE_PAIR_DISTANCE) / (double) counts.nLD, 0);
		
		if (cancelIfNotGood && LD < RNA.LOOP_DISTANCE*0.75) {
			// Bases would be drawn "too" near (0.75 is heuristic)
			return false;
		}
		
		double len = 0;
		{
			int b = firstBase;
			while (b <= lastBase) {
				coords[b].x = P0.x + vector.x*len/vectorNorm;
				coords[b].y = P0.y + vector.y*len/vectorNorm;
				int l = _listeBases.get(b).getElementStructure();
				if (b < l && l < lastBase) {
					len += RNA.BASE_PAIR_DISTANCE;
					b = l;
				} else {
					len += LD;
					b++;
				}
			}
		}
		
		return true;
	}
	
	


	/**
	 * A Bezier curve can be defined by four points,
	 * see http://en.wikipedia.org/wiki/Bezier_curve#Cubic_B.C3.A9zier_curves
	 * Here we give this four points and an array of bases indexes
	 * (which must be indexes in this RNA sequence) which will be moved
	 * to be on the Bezier curve.
	 * The bases positions are not changed in fact, instead the coords and
	 * centers arrays are modified.
	 * If cancelIfNotGood, draw only if it would look good,
	 * otherwise just return false.
	 */
	private boolean drawOnBezierCurve(
			int firstBase,
			int lastBase,
			Point2D.Double P0,
			Point2D.Double P1,
			Point2D.Double P2,
			Point2D.Double P3,
			Point2D.Double[] coords,
			Point2D.Double[] centers,
			boolean cancelIfNotGood) {
		
		UnpairedLineCounts counts = countUnpairedLine(firstBase, lastBase);
		
		// Draw the bases of the sequence along a Bezier curve
		int n = counts.total;
		
		// We choose to approximate the Bezier curve by 10*n straight lines.
		CubicBezierCurve bezier = new CubicBezierCurve(P0, P1, P2, P3, 10*n);
		double curveLength = bezier.getApproxCurveLength();
		
		
		double LD = Math.max((curveLength - counts.nBP*RNA.BASE_PAIR_DISTANCE) / (double) counts.nLD, 0);
		
		if (cancelIfNotGood && LD < RNA.LOOP_DISTANCE*0.75) {
			// Bases would be drawn "too" near (0.75 is heuristic)
			return false;
		}
		
		double[] t = new double[n+1];

		{
			double len = 0;
			int k = 0;
			int b = firstBase;
			while (b <= lastBase) {
				t[k] = len;
				k++;
				
				int l = _listeBases.get(b).getElementStructure();
				if (b < l && l < lastBase) {
					len += RNA.BASE_PAIR_DISTANCE;
					b = l;
				} else {
					len += LD;
					b++;
				}
			}
		}
		
		Point2D.Double[] sequenceBasesCoords = bezier.uniformParam(t);
		
		{
			int k = 0;
			int b = firstBase;
			while (b <= lastBase) {
				coords[b] = sequenceBasesCoords[k];
				int l = _listeBases.get(b).getElementStructure();
				if (b < l && l < lastBase) {
					b = l;
				} else {
					b++;
				}
				k++;
			}
		}

		return true;
	}
	
	
	
	private void drawOnEllipse(
			int firstBase,
			int lastBase,
			Point2D.Double P0,
			Point2D.Double P3,
			Point2D.Double[] coords,
			Point2D.Double[] centers,
			boolean reverse) {		

		UnpairedLineCounts counts = countUnpairedLine(firstBase, lastBase);
		double halfEllipseLength = RNA.BASE_PAIR_DISTANCE*counts.nBP + RNA.LOOP_DISTANCE*counts.nLD;
		double fullEllipseLength = halfEllipseLength*2;
		
		double axisA = P0.distance(P3)/2;
		double axisB = ComputeEllipseAxis.computeEllipseAxis(axisA, fullEllipseLength);
		
		if (axisB == 0) {
			// Ellipse is in fact a line, and it is much faster to draw it as such.
			drawOnStraightLine(firstBase, lastBase, P0, P3, coords, centers, false);
			return;
		}
		
		//System.out.println("Ellipse a=" + axisA + " b=" + axisB);
		
		int n = counts.total;
		
		// We choose to approximate the Bezier curve by 10*n straight lines.
		HalfEllipse curve = new HalfEllipse(axisA, axisB, 10*n);
		double curveLength = curve.getApproxCurveLength();
		
		
		double LD = Math.max((curveLength - counts.nBP*RNA.BASE_PAIR_DISTANCE) / (double) counts.nLD, 0);
		
		double[] t = new double[n+1];

		{
			double len = 0;
			int k = 0;
			int b = firstBase;
			while (b <= lastBase) {
				t[k] = len;
				k++;
				
				int l = _listeBases.get(b).getElementStructure();
				if (b < l && l < lastBase) {
					len += RNA.BASE_PAIR_DISTANCE;
					b = l;
				} else {
					len += LD;
					b++;
				}
			}
		}
		
		// Ellipse with axis = X,Y
		Point2D.Double[] sequenceBasesCoords = curve.uniformParam(t);
		
		// Compute the symmetric half-ellipse
		if (reverse) {
			AffineTransform tranform1 = new AffineTransform();
			tranform1.scale(1, -1);
			tranform1.transform(sequenceBasesCoords, 0, sequenceBasesCoords, 0, sequenceBasesCoords.length);
		}

		// Translate + rotate such that ellipse is at the right place
		AffineTransform tranform = HalfEllipse.matchAxisA(P0, P3);
		tranform.transform(sequenceBasesCoords, 0, sequenceBasesCoords, 0, sequenceBasesCoords.length);
		
		{
			int k = 0;
			int b = firstBase;
			while (b <= lastBase) {
				coords[b] = sequenceBasesCoords[k];
				//coords[b] = curve.standardParam((double) k/n);
				//System.out.println(b + " " + coords[b]);
				int l = _listeBases.get(b).getElementStructure();
				if (b < l && l < lastBase) {
					b = l;
				} else {
					b++;
				}
				k++;
			}
		}

	}
	
	/**
	 * This functions draws the RNA sequence between (including)
	 * firstBase and lastBase along a curve.
	 * The sequence may contain helices.
	 * 
	 * Bezier curve:
	 * A Bezier curve can be defined by four points,
	 * see http://en.wikipedia.org/wiki/Bezier_curve#Cubic_B.C3.A9zier_curves
	 * 
	 * Straight line:
	 * If P1 and P2 are null, the bases are drawn on a straight line.
	 * 
	 * OUT: The bases positions are not changed in fact,
	 * instead the coords and centers arrays are modified.
	 * 
	 * Argument reverse says if we should reverse the direction of inserted helices.
	 */
	private void drawAlongCurve(
			int firstBase,
			int lastBase,
			Point2D.Double P0,
			Point2D.Double P1,
			Point2D.Double P2,
			Point2D.Double P3,
			Point2D.Double[] coords,
			Point2D.Double[] centers,
			double[] angles,
			DrawRNATemplateCurveMethod curveMethod,
			boolean reverse,
			boolean straightBulges) {
		
		// First we find the bases which are directly on the Bezier curve
		ArrayList<Integer> alongBezierCurve = new ArrayList<Integer>();
		for (int depth=0, i=firstBase; i<=lastBase; i++) {
			int k = _listeBases.get(i).getElementStructure();
			if (k < 0 || k > lastBase || k < firstBase) {
				if (depth == 0) {
					alongBezierCurve.add(i);
				}
			} else {
				if (i < k) {
					if (depth == 0) {
						alongBezierCurve.add(i);
						alongBezierCurve.add(k);
					}
					depth++;
				} else {
					depth--;
				}
			}
		}
		// number of bases along the Bezier curve
		int n = alongBezierCurve.size();
		int[] alongBezierCurveArray = RNATemplateAlign.intArrayFromList(alongBezierCurve);
		if (n > 0) {
			if (curveMethod == DrawRNATemplateCurveMethod.ALWAYS_REPLACE_BY_ELLIPSES) {
				drawOnEllipse(firstBase, lastBase, P0, P3, coords, centers, reverse);
			} else if (curveMethod == DrawRNATemplateCurveMethod.SMART) {
				boolean passed;
				if (P1 != null && P2 != null) {
					passed = drawOnBezierCurve(firstBase, lastBase, P0, P1, P2, P3, coords, centers, true);
				} else {
					passed = drawOnStraightLine(firstBase, lastBase, P0, P3, coords, centers, true);
				}
				if (!passed) {
					drawOnEllipse(firstBase, lastBase, P0, P3, coords, centers, reverse);
				}
			} else {
				if (P1 != null && P2 != null) {
					drawOnBezierCurve(firstBase, lastBase, P0, P1, P2, P3, coords, centers, false);
				} else {
					drawOnStraightLine(firstBase, lastBase, P0, P3, coords, centers, false);
				}
			}
		}
		// Now use the radiate algorithm to draw the helixes
		for (int k=0; k<n-1; k++) {
			int b1 = alongBezierCurveArray[k];
			int b2 = alongBezierCurveArray[k+1];
			if (_listeBases.get(b1).getElementStructure() == b2) {
				Point2D.Double b1pos = coords[b1];
				Point2D.Double b2pos = coords[b2];
				double alpha = MiscGeom.angleFromVector(b2pos.x - b1pos.x, b2pos.y - b1pos.y);
				rna.drawLoop(b1,
						 b2,
						 (b1pos.x + b2pos.x)/2,
						 (b1pos.y + b2pos.y)/2,
						 alpha - Math.PI / 2,
						 coords,
						 centers,
						 angles,
						 straightBulges);
				if (reverse) {
					Point2D.Double symAxisVect = new Point2D.Double(coords[b2].x - coords[b1].x, coords[b2].y - coords[b1].y); 
					symmetric(coords[b1], symAxisVect, coords, b1, b2);
					symmetric(coords[b1], symAxisVect, centers, b1, b2);
				}
			}
		}	
	}
	
	
	/**
	 * Compute the symmetric of all the points from first to last in the points array
	 * relative to the line that goes through p and has director vector v.
	 * The array is modified in-place.
	 */
	private static void symmetric(
			Point2D.Double p,
			Point2D.Double v,
			Point2D.Double[] points,
			int first,
			int last) {
		// ||v||^2
		double lv = v.x*v.x + v.y*v.y;
		for (int i=first; i<=last; i++) {
			// A is the coordinates of points[i] after moving the origin at p
			Point2D.Double A = new Point2D.Double(points[i].x - p.x, points[i].y - p.y);
			// Symmetric of A
			Point2D.Double B = new Point2D.Double(
					-(A.x*v.y*v.y - 2*A.y*v.x*v.y - A.x*v.x*v.x) / lv,
					 (A.y*v.y*v.y + 2*A.x*v.x*v.y - A.y*v.x*v.x) / lv);
			// Change the origin back
			points[i].x = B.x + p.x;
			points[i].y = B.y + p.y;
		}
	}
	
	private void computeHelixTranslations(
			Tree<RNANodeValueTemplate> tree, // IN
			Map<RNATemplateHelix, Point2D.Double> translateVectors, // OUT (must be initialized)
			RNATemplateMapping mapping,      // IN
			Point2D.Double parentDeltaVector, // IN
			boolean straightBulges
	) {
		RNANodeValueTemplate nvt = tree.getValue();
		Point2D.Double newDeltaVector = parentDeltaVector;
		if (nvt instanceof RNANodeValueTemplateBasePair) {
			RNATemplateHelix helix = ((RNANodeValueTemplateBasePair) nvt).getHelix();
			if (! translateVectors.containsKey(helix)) {
				translateVectors.put(helix, parentDeltaVector);
				int[] basesInHelixArray;
				if (mapping.getAncestor(helix) != null) {
					basesInHelixArray = RNATemplateAlign.intArrayFromList(mapping.getAncestor(helix));
					Arrays.sort(basesInHelixArray);
				} else {
					basesInHelixArray = new int[0];
				}
				Point2D.Double helixDeltaVector = computeLengthIncreaseDelta(basesInHelixArray, helix, straightBulges);
				newDeltaVector = new Point2D.Double(parentDeltaVector.x+helixDeltaVector.x, parentDeltaVector.y+helixDeltaVector.y);
			} 
		}
		for (Tree<RNANodeValueTemplate> subtree: tree.getChildren()) {
			computeHelixTranslations(subtree, translateVectors, mapping, newDeltaVector, straightBulges);
		}
	}
	
	private Map<RNATemplateHelix, Point2D.Double> computeHelixTranslations(
			Tree<RNANodeValueTemplate> tree, // IN
			RNATemplateMapping mapping,       // IN
			boolean straightBulges
	) {
		Map<RNATemplateHelix, Point2D.Double> translateVectors = new HashMap<RNATemplateHelix, Point2D.Double>();
		computeHelixTranslations(tree, translateVectors, mapping, new Point2D.Double(0,0), straightBulges);
		return translateVectors;
	}
}