/*  
 *  src/org/metaqtl/algo/QtlProjAlgorithm.java
 * 
 *  $Author: veyrieras $, $Date: 2006-10-17 15:41:55 $, $Version$
 *
 *  MetaQTL : a Java package for QTL Meta-Analysis
 * 
 *  Copyright (C) 2003-2006  Jean-Baptiste Veyrieras, INRA, France.
 *
 *  This program 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 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program 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 this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA * 
 */
package org.metaqtl.algo;

import java.util.ArrayList;

import org.metaqtl.CMarkerSequence;
import org.metaqtl.ChromCluster;
import org.metaqtl.Chromosome;
import org.metaqtl.IMetaQtlConstants;
import org.metaqtl.MapCluster;
import org.metaqtl.Qtl;
import org.metaqtl.Qtl.QTLProj;
import org.metaqtl.bio.IBioGenome;
import org.metaqtl.bio.IBioLocus;
import org.metaqtl.numrec.NumericalUtilities;
import org.metaqtl.util.MappingFunction;
import org.metaqtl.util.OffspringManager;
/**
 * This class defines some methods to project QTL positions from 
 * a genetic map to another.
 */
public class QtlProjAlgorithm extends MetaAlgorithm {
	
	private int _qidx;
	
	private int _chridx;
	
	private int _clustidx;
	/**
	 * The cluster of chromosomes on which the projection would be performed.
	 */
	private MapCluster mapCluster=null;
	/**
	 * 
	 */
	private Qtl[][] qtlMatrix;
	/**
	 * 
	 */
	private Chromosome[] skeletons;
	/**
	 * The total number of Qtl flanking marker intervals which are
	 * conserved between the skeleton chromosomes and the other ones.
	 */
	private int nshare;
	
	private double alpha=0.05;
	
	private double ratio=0.5;
	/**
	 * The mode of compuation of
	 * the standard deviation of the QTL positions.
	 */
	public int sd_mode = IMetaQtlConstants.CI_MODE_MAXIMAL;
	/**
	 * Creates a new instance of <code>QtlProjAlgorithm</code> for
	 * the given configuration of projection.
	 * @param maps the maps from which the Qtl would to be projected.
	 * @param skeleton the map on which the Qtl would be projected.
	 * @param dubious an optional list of dubious marker which would not
	 *        be taken into account in the projection.
	 */
	public QtlProjAlgorithm(IBioGenome[] maps, IBioGenome skeleton, IBioLocus[] dubious) {
		/* First, clusterize the linkage groups */
		mapCluster = new MapCluster();
		mapCluster.setDubiousMarker(dubious);
		for(int i=0;i<maps.length;i++)
			mapCluster.addMap(maps[i], false);
		if (skeleton!=null)
			mapCluster.addMap(skeleton, true);
	}
	/**
	 * 
	 * @return
	 */
	public Chromosome[] getResult() {
		ChromCluster[] chrClusters = mapCluster.getClusters();
		Chromosome[] out 		   = new Chromosome[chrClusters.length];
		// Now for each skeleton chromosome add the 
		// qtl which have been projected.
		while(hasNextCluster(chrClusters)) {
			ChromCluster cluster = nextCluster(chrClusters);
			Chromosome skeleton	 = getSkeleton();
			if (skeleton.hasQTL()) {
				int i;
				int np = getProjQtlNumber();
				Qtl[] tmp = new Qtl[skeleton.getQtlNumber()+np];
				for(i=0;i<skeleton.getQtlNumber();i++) {
					tmp[i]=skeleton.qtls[i];
				}
				for(i=0;i<np;i++) {
					tmp[i+skeleton.getQtlNumber()]=getProjQtl(i);
				}
				skeleton.qtls = tmp;
			} else {
				skeleton.qtls = getProjQtls();
			}
			out[_clustidx-1]=skeleton;
		}
		return out;
	}
	/**
	 * This routine performs the projection of the QTL from the maps onto the 
	 * skeleton and returns the chromosomes on which the projection has been
	 * done.
	 
	 * @return the array of <code>Chromosome</code> on which the projection has been
	 * done
	 */
	public void run() {
		
		if (mapCluster != null) {
			
			ChromCluster[] chrClusters = mapCluster.getClusters();
			
			//System.out.println("Init Workspace");
			initWorkSpace(chrClusters);
			//System.out.println("Init Projection");
			initProjection(chrClusters);
			//System.out.println("Check Projection");
			checkProjection(chrClusters);
			//System.out.println("Do Projection");
			doProjection(chrClusters);
			
		}
	}
	
	/**
	 * @return
	 */
	private Qtl[] getProjQtls() {
		Qtl[] qtls = null;
		ArrayList list = new ArrayList();
		
		if (qtlMatrix[_clustidx-1] == null) return null;
		
		for(int i=0;i<qtlMatrix[_clustidx-1].length;i++)
			if (!qtlMatrix[_clustidx-1][i].isIgnore())
				list.add(qtlMatrix[_clustidx-1][i]);
			
		if (list.size()>0) {
		
			qtls = new Qtl[list.size()];
			list.toArray(qtls);
			
		}
		
		
		return qtls;
	}
	/**
	 * @param i
	 * @return
	 */
	private Qtl getProjQtl(int i) {
		int np=0;
		for(int j=0;j<qtlMatrix[_clustidx-1].length;j++) {
			if (!qtlMatrix[_clustidx-1][j].isIgnore()) {
				if (np == i) return qtlMatrix[_clustidx-1][j]; 
				np++;
			}
		}
		return null;
	}
	/**
	 * @return
	 */
	private int getProjQtlNumber() {
		int np=0;
		for(int i=0;i<qtlMatrix[_clustidx-1].length;i++) {
			if (!qtlMatrix[_clustidx-1][i].isIgnore())
				np++;
		}
		return np;
	}
	/**
	 * 
	 * @param chromClusters
	 */
	private void initWorkSpace(ChromCluster[] chromClusters) {
		
		_qidx     = 0;
		_chridx   = 0;
		_clustidx = 0;
		nshare 	  = 0;
		qtlMatrix = new Qtl[chromClusters.length][];
		skeletons = new Chromosome[chromClusters.length];
		
	}
	
	/**
	 * Initialize the projection.
	 * @param chromClusters
	 */
	private void initProjection(ChromCluster[] chromClusters) {
		
		while(hasNextCluster(chromClusters)) {
			
			ChromCluster cluster = nextCluster(chromClusters);
			
			initCluster(cluster);
			
		}
			
		
	}
	/**
	 * For a given cluster of chromosomes this method project which
	 * must contain a skeleton chromosome this method project the QTL
	 * of the other chromosomes on the skeleton one. 
	 * 
	 * @param cluster a cluster of chromosome which must contain a skeleton.
	 * @return the modified skeleton chromosome.
	 * 
	 * @see Chromosome
	 */
	private void initCluster(ChromCluster cluster) {
		int nqtl;
		Chromosome[] chromosomes;
		Chromosome skeleton;
		
		if (cluster==null) return;
		
		// Get the chromosomes of the cluster.  
		chromosomes = cluster.getClusterMembers();
		// init the skeleton
		initSkeleton(chromosomes);
		// get the skeleton
		skeleton = getSkeleton();
		// Count the number of QTL which have to be
		// projected on the skeleton 
		nqtl = getQtlNumber(chromosomes,skeleton);
		if (nqtl == 0) return;
		// Initialize the Qtl matrix
		initQtlMatrix(nqtl);
		
		while(hasNextChromosome(chromosomes)) { // Loop on chromosomes
			// Get the next chromosome
			Chromosome chromosome = nextChromosome(chromosomes);
			// Get the QTL to project for this chromosome
			Qtl[] qtls = initChromosome(chromosome,skeleton);
			// Store these QTL into the QTL Matrix.
			setQtlMatrix(qtls);
		}
	}
	
	/** 
	 * @param skeleton
	 * @param chromosome
	 * @param qtls
	 */
	private Qtl[] initChromosome(Chromosome chromosome, Chromosome skeleton)
	{
		int j,k,l;
		int tmp,ch_mrkIdxL,sk_mrkIdxL,ch_mrkIdxR,sk_mrkIdxR;
		boolean swap;
		CMarkerSequence cms;
		
		// First of all compute the standard deviations of the 
		// Qtl positions for the given mode.
		chromosome.computeQtlSD(sd_mode);
		// Get a copy of the Qtls of this chromosome
		Qtl[] qtls = chromosome.getQtls();
		// Get the common marker sequences between the skeleton and
		// the chromosome
		cms = CMSAlgorithm.run(chromosome, skeleton);
		
		if (cms==null) {
			if (isLoggerEnable()) {
				getLogger().println("[ WARNING ] No common markers for " + chromosome.getMapName() + " on chromosome " + chromosome.getName());
			}
			for(l=0;l<qtls.length;l++) {
				qtls[l].setIgnore(true);				
			}
			return qtls;
		}
		// First get the first and the last common markers on
		// the chromosome
		ch_mrkIdxL=cms.mcidx[0][0];
		ch_mrkIdxR=cms.mcidx[0][cms.nmc-1];
		// Remove the QTL which could be before the first marker and/or
		// after the last marker.
		for(l=0;l<qtls.length;l++) {
			if (ch_mrkIdxL > 0 && chromosome.areMarkersAroundQTL(0, ch_mrkIdxL, l))
				qtls[l].setIgnore(true);
			if (ch_mrkIdxR < chromosome.nm-1 && chromosome.areMarkersAroundQTL(ch_mrkIdxR, chromosome.nm-1, l))
				qtls[l].setIgnore(true);
			if (qtls[l].isIgnore() && isLoggerEnable()) {
				getLogger().println("[ WARNING ] Unable to project QTL " + qtls[l].getName() + " on chromosome " + chromosome.getName());
			}
		}
		
		for(j=0;j<cms.ncs;j++) { // loop on common sequences
			
			for(k=0;k<cms.css[j]-1;k++) {  // loop on marker intervals inside a commone sequence
				
				sk_mrkIdxL = cms.idx2[j][k];    /* skeleton left flanking marker indice */
				sk_mrkIdxR = cms.idx2[j][k+1];  /* skeleton right flanking marker indice */
				ch_mrkIdxL = cms.idx1[j][k];    /* chromosome left flanking marker indice */
				ch_mrkIdxR = cms.idx1[j][k+1];  /* chromosome right flanking marker indice */
				
				swap = false;
				
				/* Swap indices if required */
				if (sk_mrkIdxL > sk_mrkIdxR) {
					tmp        = sk_mrkIdxR; 
					sk_mrkIdxR = sk_mrkIdxL;
					sk_mrkIdxL = tmp;
					swap       = true;
				}
				/* Swap indices if required */
				if (ch_mrkIdxL > ch_mrkIdxR) {
					tmp        = ch_mrkIdxR; 
					ch_mrkIdxR = ch_mrkIdxL;
					ch_mrkIdxL = tmp;
					swap   = !swap && true;
				}
				
				/*
				 * Get the weigths of the interval for each QTL 
				 */
				for(l=0;l<chromosome.qtls.length;l++) {
					
					if (qtls[l].isIgnore()) continue;
					
					if (qtls[l].proj == null) {
						qtls[l].proj           = qtls[l].new QTLProj();
						qtls[l].proj.mapOrigin = chromosome.getMapName();
					}
					
					updateQtlScale(qtls[l], chromosome, ch_mrkIdxL, ch_mrkIdxR, skeleton, sk_mrkIdxL, sk_mrkIdxR);
					/*
					 * Test if the given marker configuration surrounds
					 * the l^th QTL. Logically we expect at most a single positive 
					 * response...
					 */
					if (!qtls[l].proj.shareFlanking) {
						
						qtls[l].proj.shareFlanking					
							= chromosome.areMarkersAroundQTL(ch_mrkIdxL, ch_mrkIdxR, l);
						
						if (qtls[l].proj.shareFlanking) {
							
							setQtlFlanking(qtls[l],chromosome, ch_mrkIdxL, ch_mrkIdxR, skeleton, sk_mrkIdxL, sk_mrkIdxR,swap);
							
						}
					}							
				}
				
		 	} /* end of loop on inside common marker sequence */
		 	
		 } /* end of loop on common marker sequence */
		
		return qtls;
	}
	/**
	 * 
	 * @param chromClusters
	 */
	public void checkProjection(ChromCluster[] chromClusters) {
		
		while(hasNextCluster(chromClusters)) {
			
			ChromCluster cluster = nextCluster(chromClusters);
			
			checkCluster(cluster);
			
		}	
		
	}
	/**
	 * 
	 * @param cluster
	 */
	private void checkCluster(ChromCluster cluster) {
		
		if (cluster==null) return;
		
		// Get the members of the cluster 
		Chromosome[] chromosomes = cluster.getClusterMembers();
		// Get the skeleton
		Chromosome skeleton 	 = getSkeleton();
		
		while(hasNextChromosome(chromosomes)) { // Loop on chromosomes
			// Get the next chromosome
			Chromosome chromosome = nextChromosome(chromosomes);
			
			// Check the projection for this chromosome
			checkChromosome(chromosome,skeleton);
		}
		
	}
	/**
	 * 
	 * @param chromosome
	 * @param skeleton
	 */
	private void checkChromosome(Chromosome chromosome, Chromosome skeleton)
	{
		
		Qtl[] qtls = getQtlToProject(chromosome);
		
		for(int i=0;i<qtls.length;i++) {
			
			if (qtls[i].isIgnore()) continue;
			
			QTLProj projPar = qtls[i].proj;
			
			if (qtls[i].proj != null) {
				// If the flanking marker interval around the Qtl has not
				// been found in any common sequences in the previous round
				// then try to find another flanking interval.
				// If the method fails to find another flanking marker interval
				// then the projection is said to be impossible.
				if (!projPar.shareFlanking) {
					//System.out.println("Find Falnking");
					qtls[i].setIgnore(!findFlanking(chromosome,skeleton,qtls[i]));
				}
			} else {
				// this means that any common marker sequence have been found at
				// the previous step. Then the projection is not possible.
				qtls[i].setIgnore(true);
				if (isLoggerEnable()) {
					getLogger().println("[ WARNING ] : Unable to project QTL " + qtls[i].getName() +  " on chromosome " + chromosome.getName());
					getLogger().flush();
				}
				
			}
			// If the previous steps have lead to a valid QTL projection
			// context then check this context in terms of dilatation
			// parameters and if not valid try to propose another context.
			if (!qtls[i].isIgnore()) {
			
				checkFlanking(chromosome,skeleton,qtls[i]);
				
			} else {
				
				if (isLoggerEnable()) {
					getLogger().println("[ WARNING ] : Unable to project QTL " + qtls[i].getName() +  " on chromosome " + chromosome.getName());
					getLogger().flush();
				}
				
			}
			
		}
		
	}
	/**
	 * 
	 * @return
	 */
	private boolean findFlanking(Chromosome chromosome, Chromosome skeleton, Qtl qtl) {
		int tmp,ch_mrkIdxL,sk_mrkIdxL,ch_mrkIdxR,sk_mrkIdxR;
		boolean frame,swap=false;
		CMarkerSequence cms;
		boolean ok=false;
		
		/*
		 * At this step, this means that the previous loop on
		 * common marker sequences has found any common sequence
		 * that includes the QTL. Therefore we have to look if there
		 * are two continguous sequences which upper and lower boundaries
		 * are flanking the QTL.
		 */
		
		cms = CMSAlgorithm.run(chromosome, skeleton);
		
		ch_mrkIdxL = cms.getMarkerIdx(0, 0, 0);
		
		if (ch_mrkIdxL >= qtl.m2) {
			// Test if there is one marker at the left which is
			// not in a common sequence
			int j,l;
			for(j=0;j<cms.nmc;j++)
				if (cms.mcidx[0][j] <= qtl.m1)
					break;
				
			for(l=0;l<cms.ncs;l++) {
				for(int m=0;m<cms.css[l];m++) {
					if (cms.mcidx[1][j]==cms.getMarkerIdx(0, l, m)) break;
				}
			}
			if (l == cms.ncs) {
				ch_mrkIdxR = ch_mrkIdxL;
				ch_mrkIdxL = j;
				sk_mrkIdxL = cms.mcidx[1][j];        /* skeleton left flanking marker indice */
				sk_mrkIdxR = cms.getMarkerIdx(1,0,0);/* skeleton right flanking marker indice */
				/* Swap indices if required */
				if (sk_mrkIdxL > sk_mrkIdxR) {
					tmp        = sk_mrkIdxR; 
					sk_mrkIdxR = sk_mrkIdxL;
					sk_mrkIdxL = tmp;
					swap       = true;
				}
				updateQtlScale(qtl,chromosome, ch_mrkIdxL, ch_mrkIdxR, skeleton, sk_mrkIdxL, sk_mrkIdxR);
				setQtlFlanking(qtl,chromosome, ch_mrkIdxL, ch_mrkIdxR, skeleton, sk_mrkIdxL, sk_mrkIdxR, swap);
				return true;
			}
		}
		
		ch_mrkIdxR = cms.getMarkerIdx(0, cms.ncs-1, cms.css[cms.ncs-1]-1);
		
		if (ch_mrkIdxR <= qtl.m1) {
			// Test if there is one marker at the left which is
			// not in a common sequence
			int j,l;
			for(j=0;j<cms.nmc;j++)
				if (cms.mcidx[0][j] >= qtl.m2)
					break;
			for(l=0;l<cms.ncs;l++) {
				for(int m=0;m<cms.css[l];m++) {
					if (cms.mcidx[0][j]==cms.getMarkerIdx(0, l, m)) break;
				}
			}
			if (l == cms.ncs) {
				ch_mrkIdxL = ch_mrkIdxR;
				ch_mrkIdxR = cms.mcidx[0][j];
				sk_mrkIdxL = cms.getMarkerIdx(1,cms.ncs-1, cms.css[cms.ncs-1]-1);        /* skeleton left flanking marker indice */
				sk_mrkIdxR = cms.mcidx[1][j];/* skeleton right flanking marker indice */
				/* Swap indices if required */
				if (sk_mrkIdxL > sk_mrkIdxR) {
					tmp        = sk_mrkIdxR; 
					sk_mrkIdxR = sk_mrkIdxL;
					sk_mrkIdxL = tmp;
					swap       = true;
				}
				updateQtlScale(qtl,chromosome, ch_mrkIdxL, ch_mrkIdxR, skeleton, sk_mrkIdxL, sk_mrkIdxR);
				setQtlFlanking(qtl,chromosome, ch_mrkIdxL, ch_mrkIdxR, skeleton, sk_mrkIdxL, sk_mrkIdxR, swap);
				return true;
			}
		}
		
		for(int j=0;j<cms.ncs-1;j++) { // loop on common sequences
			// Get the upper marker idx on the chromosome.
			ch_mrkIdxL = cms.idx1[j][cms.css[j]-1];
			// Get the lower marker idx on the chromosome.
			ch_mrkIdxR = cms.idx1[j+1][0];
			// does the frame of the common sequences are the sames ?
			frame = true;//((cms.frames[j] && cms.frames[j+1]) || (!cms.frames[j] && !cms.frames[j+1]));
			
			// Does this marker interval includes the QTL ?
			if (frame && (ch_mrkIdxL <= qtl.m1 &&  ch_mrkIdxR >= qtl.m2)) {
				ok = true;
				
				sk_mrkIdxL = cms.idx2[j][cms.css[j]-1];    /* skeleton left flanking marker indice */
				sk_mrkIdxR = cms.idx2[j+1][0]; 			   /* skeleton right flanking marker indice */
				/* Swap indices if required */
				if (sk_mrkIdxL > sk_mrkIdxR) {
					tmp        = sk_mrkIdxR; 
					sk_mrkIdxR = sk_mrkIdxL;
					sk_mrkIdxL = tmp;
					swap       = true;
				}
				
				updateQtlScale(qtl,chromosome, ch_mrkIdxL, ch_mrkIdxR, skeleton, sk_mrkIdxL, sk_mrkIdxR);
				setQtlFlanking(qtl,chromosome, ch_mrkIdxL, ch_mrkIdxR, skeleton, sk_mrkIdxL, sk_mrkIdxR,swap);

				// exit ok
				break;
			} else if (ch_mrkIdxL > qtl.m1)
				break;
		}
		
		return ok;
	}
	
	private void updateQtlScale(Qtl qtl, Chromosome chromosome, int m1, int m2, Chromosome skeleton, int s1, int s2) {
		// Compute the ratio between the two intervals.
		double rho = getRatio(chromosome,m1,m2,skeleton,s1,s2);
		// get the weight of this interval.
		double w   = getQtlProba(qtl,chromosome,m1,m2);
		qtl.proj.mapScale  += w * rho;
		qtl.proj.mapWeight += w;
	}
	
	private double getQtlProba(Qtl qtl, Chromosome chromosome, int m1, int m2) {
		if (m1 <= m2)
			return chromosome.getQTLInMrkIntProb(m1, m2, qtl);
		else
			return chromosome.getQTLInMrkIntProb(m2, m1, qtl);
		
	}
	
	private void setQtlFlanking(Qtl qtl, Chromosome chromosome, int m1, int m2, Chromosome skeleton, int s1, int s2, boolean swap) {
		qtl.proj.o_mrkIdxL = m1;
		qtl.proj.o_mrkIdxR = m2;
		qtl.proj.s_mrkIdxL = s1;
		qtl.proj.s_mrkIdxR = s2;
		qtl.proj.swap      = swap;
		qtl.proj.intScale  = getRatio(chromosome,m1,m2,skeleton,s1,s2);;
		// increment the total number of shared flanking marker intervals.
		nshare++;
	}
	/**
	 * 
	 * @param chromosome
	 * @param skeleton
	 * @param qtl
	 */
	private void checkFlanking(Chromosome chromosome, Chromosome skeleton, Qtl qtl) {
		int m1,m2,s1,s2;
		double ratio,pvalue;
		
		QTLProj projPar = qtl.proj;
		
		m1=projPar.o_mrkIdxL;
		m2=projPar.o_mrkIdxR;
		s1=projPar.s_mrkIdxL;
		s2=projPar.s_mrkIdxR;
		
		pvalue = getPvalue(chromosome,m1,m2,skeleton,s1,s2);
		ratio  = getRatio(chromosome,m1,m2,skeleton,s1,s2);
		
		// This is totaly empiric : we try to get another interval for
		// the QTL which ratio is beyond 1.0 (compression) and the pvalue
		// is beyond a corrected threshold alpha (multiple test).
		
		if (ratio < this.ratio && pvalue < this.alpha) {
			
			qtl.setIgnore(!updateFlanking(chromosome, skeleton, qtl, 1.0, alpha));
			
			if (qtl.isIgnore()) {
				if (isLoggerEnable()) {
					//getLogger().println(chromosome.getMarkerName(m1)+" < > "+chromosome.getMarkerName(m2));
					//getLogger().println(ratio+" < > "+pvalue);
					getLogger().println("[ WARNING ] : Unable to project QTL " + qtl.getName() +  " on chromosome " + chromosome.getName());
					getLogger().flush();
				}
			}
			
		}
		
	}
	/**
	 * 
	 * @param chromosome
	 * @param skeleton
	 * @param qtl
	 * @param ratio
	 * @param pvalue
	 * @return
	 */
	private boolean updateFlanking(Chromosome chromosome, Chromosome skeleton, Qtl qtl, double ratio, double pvalue) {
		boolean ok = false;
		int iter,tmp,c1,s1,c2,s2,m1,m2,cs1,cs2,m1x,m2x;
		int[] idx;
		double r,p;
		boolean leftEnd,rightEnd,swap=false;
		CMarkerSequence cms;
		
		QTLProj projPar = qtl.proj;
		
		cs1=cs2=m1x=m2x=0;
		c1 = m1 = projPar.o_mrkIdxL;
		c2 = m2 = projPar.o_mrkIdxR;
		s1 = projPar.s_mrkIdxL;
		s2 = projPar.s_mrkIdxR;
		
		cms = CMSAlgorithm.run(chromosome, skeleton);
		
		// Get the common sequence indices on which
		// the flanking marker are.
		leftEnd = rightEnd = false;
		
		idx = cms.getCSIdx(0, m1);
		if (idx == null) {
			leftEnd = true;
		} else {
			cs1 = idx[0];
			m1x = idx[1];
		}
		idx = cms.getCSIdx(0, m2);
		if (idx == null) {
			rightEnd = true;
		} else {
			cs2 = idx[0];
			m2x = idx[1];
		}
		
		iter=0;
		
		while(!ok) {
			
			iter++;
			
			if (iter%2==0 && !leftEnd) {
				// Go to left
				idx=cms.getLeftCSMarker(cs1,m1x);
				if (idx[0]==cs1 && idx[1]==m1x) leftEnd=true;
				cs1 = idx[0];
				m1x = idx[1];
				c1 = cms.getMarkerIdx(0, cs1, m1x);
				s1 = cms.getMarkerIdx(1, cs1, m1x);    /* skeleton left flanking marker indice */				
			} else if (!rightEnd){
				// Go to rigth
				idx=cms.getRightCSMarker(cs2,m2x);
				if (idx[0]==cs2 && idx[1]==m2x) rightEnd=true;
				cs2 = idx[0];
				m2x = idx[1];
				c2 = cms.getMarkerIdx(0, cs2, m2x);
				s2 = cms.getMarkerIdx(1, cs2, m2x);    /* skeleton right flanking marker indice */
			}
			
			if (leftEnd && rightEnd) break;
			
			/* Swap indices if required */
			if (c1 > c2) {
				tmp = c2; 
				c2  = c1;
				c1  = tmp;
				swap = true;
			}
			/* Swap indices if required */
			if (s1 > s2) {
				tmp = s2; 
				s2  = s1;
				s1  = tmp;
				swap = !swap && true;
			}
			
			//does the frame of the common sequnces are the sames ?
			//if (!cms.haveSameFrame(cs1,cs2)) {
				//continue;
			//}
			//Check if the new configuration fill
			// the required conditions
			// Compute the ratio between the two intervals.
			r = getRatio(chromosome,c1,c2,skeleton,s1,s2);
			p = getPvalue(chromosome,c1,c2,skeleton,s1,s2);
			
			if (r >= ratio || p > pvalue) ok=true;
			//System.out.println(iter);
			
		}
		if (ok) {
			// Set the flanking marker interval
			setQtlFlanking(qtl,chromosome,c1,c2,skeleton,s1,s2,swap);
			//Recompute the projection scale parameters.
			projPar.mapScale  = 0.0;
			projPar.mapWeight = 0.0;
			updateQtlScale(qtl,chromosome,c1,c2,skeleton,s1,s2);
			// Loop on the cms at the left of the flanking interval
			do {
				idx = cms.getLeftCSMarker(cs1,m1x);
				if (idx[0]==cs1 && idx[1]==m1x) break;
				c1 = cms.getMarkerIdx(0, cs1, m1x);
				c2 = cms.getMarkerIdx(0, idx[0], idx[1]);
				s1 = cms.getMarkerIdx(1, cs1, m1x);
				s2 = cms.getMarkerIdx(1, idx[0], idx[1]);
				updateQtlScale(qtl,chromosome,c1,c2,skeleton,s1,s2);
				cs1 = idx[0];
				m1x = idx[1];
			} while(0==0);
			// Loop on the cms at the right of the flanking interval
			do {
				idx = cms.getRightCSMarker(cs2,m2x);
				if (idx[0]==cs2 && idx[1]==m2x) break;
				c1 = cms.getMarkerIdx(0, cs2, m2x);
				c2 = cms.getMarkerIdx(0, idx[0], idx[1]);
				s1 = cms.getMarkerIdx(1, cs2, m2x);
				s2 = cms.getMarkerIdx(1, idx[0], idx[1]);
				updateQtlScale(qtl,chromosome,c1,c2,skeleton,s1,s2);
				cs2 = idx[0];
				m2x = idx[1];
			} while(0==0);			
		}
		return ok;
	}
	/**
	 * 
	 * @param chromClusters
	 */
	public void doProjection(ChromCluster[] chromClusters) {
		
		while(hasNextCluster(chromClusters)) {
			
			ChromCluster cluster = nextCluster(chromClusters);
			
			doCluster(cluster);
			
		}	
		
	}
	/**
	 * 
	 * @param cluster
	 */
	private void doCluster(ChromCluster cluster) {
		
		if (cluster==null) return;
		
		// Get the members of the cluster 
		Chromosome[] chromosomes = cluster.getClusterMembers();
		// Get the skeleton
		Chromosome skeleton 	 = getSkeleton();
		
		while(hasNextChromosome(chromosomes)) { // Loop on chromosomes
			// Get the next chromosome
			Chromosome chromosome = nextChromosome(chromosomes);
			
			// Check the projection for this chromosome
			doChromosome(chromosome,skeleton);
		}
		
	}
	
	/**
	 * This methods project the QTL of the chromosome <code>chromosome</code>
	 * on the skeleton chromosome <code>skeleton</code>.  
	 * 
	 * @param skeleton the skeleton chromosome
	 * @param chromosome the chromosome.
	 * @param m the current first QTL indice.
	 */
	private void doChromosome(Chromosome chromosome, Chromosome skeleton) {
		int i;
		int ch_mrkIdxL,sk_mrkIdxL,ch_mrkIdxR,sk_mrkIdxR;
		double d,rho,sk_d,x_l,tmp,dic1,dic2;
		boolean swap;
		
		Qtl[] qtls = getQtlToProject(chromosome);
		
		for(i=0;i<qtls.length;i++) {
			
			if (qtls[i].isIgnore()) continue;
				
			QTLProj projPar = qtls[i].proj;
			ch_mrkIdxL = projPar.o_mrkIdxL;
			ch_mrkIdxR = projPar.o_mrkIdxR;
			/* Get the corresponding marker indices on the skeleton */
			sk_mrkIdxL = projPar.s_mrkIdxL;
			sk_mrkIdxR = projPar.s_mrkIdxR;
			swap       = projPar.swap;
			/*
			 * Then, it is the most simple case. We just have to 
			 * compute the ratio between the skeleton flanking marker
			 * interval and this of the original map, namely t;
			 * Then the position of the QTL on the skeleton is given
			 * by :
			 * 		 sk_xl = ch_xl*sk_d/ch_d   
			 * where :
			 * 
			 * 	sk_xl is the distance to the left flanking skeleton marker.
			 *  ch_xl is the distance to the left flanking chromosome marker.
			 *  sk_d  is the distance of the skeleton flanking marker interval.
			 *  ch_d is the distance of the chromosome flanking marker interval.
			 * 
			 * Note that if the flanking marker intervals are in reversed frame
			 * between the skeleton and the chromosome marker interval then ch_xl is
			 * the distance between the QTL position and the right flanking marker
			 * on the chromosome (see above).  
			 */
			/* First get the chromosome left flanking marker indice */
			if (swap) {
				/* Compute the distance from the flagged right marker indice
				 * to the QTL.*/
				d = chromosome.getDistance(ch_mrkIdxR) - qtls[i].position;
				if (d < 0.0) {
					System.out.println("SWAP::ERROR "+ qtls[i].name+" "+ch_mrkIdxL + " " + ch_mrkIdxR + " " + chromosome.getName());
					System.out.println(qtls[i].m1 + " " + qtls[i].m2);
					System.out.println(qtls[i].position + " " + chromosome.getDistance(ch_mrkIdxR));
					continue;
				}
				x_l = MappingFunction.recombination(d, chromosome.mfc, chromosome.mut);
				x_l = MappingFunction.distance(x_l, IMetaQtlConstants.MAPPING_FUNCTION_HALDANE, IMetaQtlConstants.MAPPING_UNIT_M);
			} else {
				/* Compute the distance from the flagged left marker indice
				 * to the QTL.*/
				d   = qtls[i].position - chromosome.getDistance(ch_mrkIdxL);
				if (d < 0.0) {
					System.out.println("ERROR "+ qtls[i].name+" "+ch_mrkIdxR + " " + ch_mrkIdxL + " " + chromosome.nm);
					continue;
				}
				x_l = MappingFunction.recombination(d, chromosome.mfc, chromosome.mut);
				x_l = MappingFunction.distance(x_l, IMetaQtlConstants.MAPPING_FUNCTION_HALDANE, IMetaQtlConstants.MAPPING_UNIT_M);
			}
			/* Compute the haldane distance between ch_mrkIdxL and ch_mrkIdxR */
			rho  = projPar.intScale;
			// Get the distance of the interval
			sk_d = skeleton.getDistanceBetween(sk_mrkIdxL, sk_mrkIdxR, IMetaQtlConstants.MAPPING_FUNCTION_HALDANE, IMetaQtlConstants.MAPPING_UNIT_M);
			/* Set the flanking marker indices */
			qtls[i].m1 = sk_mrkIdxL;
			qtls[i].m2 = sk_mrkIdxR;
			/* Set the flanking marker distances */
			qtls[i].c1 =  x_l * rho;
			qtls[i].c2 =  sk_d - qtls[i].c1;
			sk_d = skeleton.getDistance(sk_mrkIdxL, IMetaQtlConstants.MAPPING_FUNCTION_HALDANE, IMetaQtlConstants.MAPPING_UNIT_M);
			qtls[i].position = sk_d + qtls[i].c1;
			/* Now convert the position for the skeleton mapping context */
			qtls[i].position = MappingFunction.recombination(qtls[i].position, IMetaQtlConstants.MAPPING_FUNCTION_HALDANE, IMetaQtlConstants.MAPPING_UNIT_M);
			qtls[i].position = MappingFunction.distance(qtls[i].position, skeleton.mfc, skeleton.mut);
			/* Convert c1 and c2 into recombination fractions */
			qtls[i].c1 = MappingFunction.recombination(
											qtls[i].c1,
											IMetaQtlConstants.MAPPING_FUNCTION_HALDANE,
											IMetaQtlConstants.MAPPING_UNIT_M);
			qtls[i].c2 = MappingFunction.recombination(
											qtls[i].c2,
											IMetaQtlConstants.MAPPING_FUNCTION_HALDANE,
											IMetaQtlConstants.MAPPING_UNIT_M);
			// Compute the map scale projection parameter 
			if (qtls[i].proj.mapWeight > 0.0)
				rho = qtls[i].proj.mapScale/qtls[i].proj.mapWeight;
			else rho = 1.0;
			qtls[i].proj.mapScale = rho;
			
			// If defined scale the confidence interval
			if (qtls[i].ci != null) {
				/* If a CI is defined then apply the average
				 * dilatation to it */
				dic1 = MappingFunction.distance(
												qtls[i].ci.ic1,
												skeleton.mfc,
												skeleton.mut);
				dic2 = MappingFunction.distance(
												qtls[i].ci.ic2,
												skeleton.mfc,
												skeleton.mut);
				if (swap) {
					tmp = dic1;
					dic1 = dic2;
					dic2 = tmp;
				}
				/* Get the absolute position of the QTL */
				tmp = qtls[i].position;
				/* Get the total distance in Morgan of the skeleton chromosome */
				/* Check if dic2 is not over the chromosome range */ 
				if (tmp + dic2 > skeleton.totd) {
					dic2 = skeleton.totd - tmp;
				}
				
				qtls[i].ci.from = Math.max(0,tmp - dic1);
				qtls[i].ci.to   = tmp + dic2;
				
				qtls[i].ci.ic1  = MappingFunction.recombination(
														dic1,
														skeleton.mfc,
														skeleton.mut);
				qtls[i].ci.ic2 = MappingFunction.recombination(
														dic2,
														skeleton.mfc,
														skeleton.mut);
			}
			
		}
		
	}
	
	
	/**
	 * @param skeleton
	 * @param s1
	 * @param s2
	 * @param chromosome
	 * @param m1
	 * @param m2
	 * @return
	 */
	private static double getPvalue(Chromosome chromosome, int m1, int m2, Chromosome skeleton, int s1, int s2) {
		double d1,d2,r2,vr2,vd2,z = 0.0;
		
		d1  = skeleton.getDistanceBetween(s1, s2, IMetaQtlConstants.MAPPING_FUNCTION_HALDANE, IMetaQtlConstants.MAPPING_UNIT_M);
		d2  = chromosome.getDistanceBetween(m1, m2, IMetaQtlConstants.MAPPING_FUNCTION_HALDANE, IMetaQtlConstants.MAPPING_UNIT_M);
		r2  = MappingFunction.recombination(d2, IMetaQtlConstants.MAPPING_FUNCTION_HALDANE, IMetaQtlConstants.MAPPING_UNIT_M);
		vr2 = OffspringManager.fisherVariance(r2, chromosome.cs, chromosome.ct, chromosome.cg);
		vd2 = MappingFunction.varianceDistance(r2, vr2, IMetaQtlConstants.MAPPING_FUNCTION_HALDANE, IMetaQtlConstants.MAPPING_UNIT_M);
		
		// This is the standardized residual assuming that the 
		// true distance is the distance of the interval on the
		// skeleton map.
		
		z=(d1-d2)*(d1-d2)/vd2;
		
		return 1.0-NumericalUtilities.gammp(0.5,0.5*z);
	}
	
	/**
	 * @param skeleton
	 * @param s1
	 * @param s2
	 * @param chromosome
	 * @param m1
	 * @param m2
	 * @return
	 */
	private static double getRatio(Chromosome chromosome, int m1, int m2, Chromosome skeleton, int s1, int s2) {
		double d1,d2,r;
		
		if (s1 <= s2)
			d1  = skeleton.getDistanceBetween(s1, s2, IMetaQtlConstants.MAPPING_FUNCTION_HALDANE, IMetaQtlConstants.MAPPING_UNIT_M);
		else
			d1  = skeleton.getDistanceBetween(s2, s1, IMetaQtlConstants.MAPPING_FUNCTION_HALDANE, IMetaQtlConstants.MAPPING_UNIT_M);
		if (m1 <= m2)
			d2  = chromosome.getDistanceBetween(m1, m2, IMetaQtlConstants.MAPPING_FUNCTION_HALDANE, IMetaQtlConstants.MAPPING_UNIT_M);
		else
			d2  = chromosome.getDistanceBetween(m2, m1, IMetaQtlConstants.MAPPING_FUNCTION_HALDANE, IMetaQtlConstants.MAPPING_UNIT_M);
		
		r=d1/d2;
		
		if (Double.isInfinite(r) || Double.isNaN(r)) r = 0.0;
		
		return r;
	}
	/**
	 * 
	 * @param clusters
	 * @return
	 */
	private boolean hasNextCluster(ChromCluster[] clusters) {
		if (_clustidx == clusters.length) {
			_clustidx=0;
			return false;
		}
		boolean next;
		next  = (_clustidx < clusters.length);
		_qidx = 0;		
		return next;
	}
	/**
	 * 
	 * @param clusters
	 * @return
	 */
	private ChromCluster nextCluster(ChromCluster[] clusters) {
		return clusters[_clustidx++];
	}
	
	private boolean hasNextChromosome(Chromosome[] chromosomes) {
		
		if (_chridx == chromosomes.length)	{
			_chridx=0;
			return false;
		}
		
		boolean next = (_chridx < chromosomes.length);
		
		if (next) {
			
			if (_chridx > 0) 
				
				_qidx += chromosomes[_chridx - 1].getQtlNumber();
			
			while (next && (!chromosomes[_chridx].hasQTL() || isSkeleton(chromosomes[_chridx], skeletons[_clustidx-1]))) {
				
				_chridx++;
				
				next = (_chridx < chromosomes.length);
				
			}
			
			if (_chridx == chromosomes.length)	{
				_chridx=0;
				next = false;
			}
			
		}	
		
		return next;
	}
	
	private Chromosome nextChromosome(Chromosome[] chromosomes) {
		return chromosomes[_chridx++];
	}
	
	private void initQtlMatrix(int nqtl) {
		
		qtlMatrix[_clustidx-1] = new Qtl[nqtl];
		
	}
	
	private void setQtlMatrix(Qtl[] qtls) {
		
		for(int i=0;i<qtls.length;i++) {
			
			qtlMatrix[_clustidx-1][_qidx + i] = qtls[i];
			
		}
	}
	/**
	 * @return
	 */
	private Qtl[] getQtlToProject(Chromosome chromosome) {
		Qtl[] qtls = new Qtl[chromosome.getQtlNumber()];
		
		for(int i=0;i<qtls.length;i++)
			qtls[i]=qtlMatrix[_clustidx-1][_qidx + i];
		
		return qtls;
	}
	
	/**
	 * @return
	 */
	private Chromosome getSkeleton() {
		return skeletons[_clustidx-1];
	}
	/**
	 * @param chromosomes
	 */
	private void initSkeleton(Chromosome[] chromosomes) {
		// minus one because the method is called after
		// nextCluster() have been called.
		skeletons[_clustidx-1] = getSkeleton(chromosomes);
	}
	/**
	 * @param chromosome
	 * @param skeleton
	 * @return
	 */
	private boolean isSkeleton(Chromosome chromosome, Chromosome skeleton) {
		if (skeleton==null) return false;
		return (chromosome.hashCode()==skeleton.hashCode());
	}
	/**
	 * @param chromosomes
	 * @return
	 */
	private int getQtlNumber(Chromosome[] chromosomes, Chromosome skeleton) {
		int nqtl=0;
		for(int i=0;i<chromosomes.length;i++) {
			if (!isSkeleton(chromosomes[i],skeleton))
				nqtl += chromosomes[i].getQtlNumber();
		}
		return nqtl;
	}
	/**
	 * @param chromosomes
	 * @return
	 */
	private Chromosome getSkeleton(Chromosome[] chromosomes) {
		int i;
		/* First get the skeleton member */
		for(i=0;i<chromosomes.length;i++)
			if (chromosomes[i].skeleton)
				break;
		if (i==chromosomes.length) return null;
		/* Set the skeleton */
		return chromosomes[i];
	}
	/**
	 * @param d
	 */
	public void setRatio(double r) {
		this.ratio = r;
	}
	
	public void setPvalue(double p) {
		this.alpha = p;
	}
}
