/*  
 *  src/org/metaqtl/Chromosome.java
 * 
 *  $Author: veyrieras $, $Date: 2006-10-17 15:41:54 $, $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;

import java.util.ArrayList;

import org.metaqtl.adapter.QtlAdapter;
import org.metaqtl.bio.IBioCrossTypes;
import org.metaqtl.bio.IBioLocus;
import org.metaqtl.numrec.NumericalUtilities;
import org.metaqtl.util.MappingFunction;
import org.metaqtl.util.OffspringManager;

/**
 * On a Chromosome there are two types of locus :
 * <ul>
 * <li> Marker </li>
 * <li> QTL (Quantitative Trait Loci) </li>
 * </ul>
 * inside a mapping context which is defined by both the mapping framework
 * (function and unit) and the population cross design used to build this
 * chromosome.
 */
public class Chromosome {
	/**
	 * The number of markers.
	 */
	public int nm;
	/**
	 * The size of the offspring population.
	 */
	public int cs=0;
	/**
	 * The type of the offspring population.
	 */
	public int ct=IBioCrossTypes.UNKNOWN;
	/**
	 * The number of generations of the offspring population.
	 */
	public int cg=0;
	/**
	 * The mapping function for this chromosome
	 */
	public int mfc;
	/**
	 * The mapping unit for this chromosome
	 */
	public int mut;
	/**
	 * The skeleton status for this member
	 */
	public boolean skeleton;
	/**
	 * The total distance of the chromosome.
	 */
	public double totd;
	/**
	 *  The adjacent recombination rate between 
	 * ordered markers */
	public double[] mr=null;
	/**
	 * The ordered marker names
	 */
	public String[] mrkNames=null;
	/**
	 * The name of the chromosome.
	 */
	public String name=null;
	/**
	 * The name of the map in which this chromosome
	 * is defined.
	 */
	public String mapName=null;
	/**
	 * The QTLs located on this chromosome.
	 */
	public Qtl[] qtls=null;
	/**
	 * The meta Qtls on this chromosome.
	 */
	public Qtl[] mqtls=null;
	/**
	 * The results of a meta-analysis of QTL.
	 */
	public MetaQtlResult[] metaQtlAnalysis=null;
	/**
	 * Probabilities defined along the chromosome.
	 */
	public double[][] proba=null;
	
	/**
	 * This methods returns the probability that the QTL be inside 
	 * the marker interval with the given right and left marker indices
	 * and the given QTL indice.
	 * The method assumes that <code>mrkIdxL &lt; mrkIdxR</code>
	 *
	 * @param mrkIdxL the left marker indice. 
	 * @param mrkIdxR the right marker indice.
	 * @param qtlIdx  the QTL indice. 
	 * @return the probability that the QTL be inside the interval.
	 */
	public double getQTLInMrkIntProb(int mrkIdxL, int mrkIdxR, int qtlIdx) {
		Qtl qtl;
		
		if (qtls==null || qtls.length <= qtlIdx) return 0.0;
		
		qtl = qtls[qtlIdx];
		
		return getQTLInMrkIntProb(mrkIdxL,mrkIdxR,qtl);
		
	}
	/**
	 * This methods returns the probability that the QTL be inside 
	 * the marker interval with the given right and left marker indices
	 * and the given QTL.
	 * The method assumes that <code>mrkIdxL &lt; mrkIdxR</code>
	 *
	 * @param mrkIdxL the left marker indice. 
	 * @param mrkIdxR the right marker indice.
	 * @param qtl     the QTL. 
	 * @return the probability that the QTL be inside the interval.
	 */
	public double getQTLInMrkIntProb(int mrkIdxL, int mrkIdxR, Qtl qtl) {
		int i;
		double mu,sigma,d,h,y1,y2,w;
		
		mu    = qtl.position;
		sigma = qtl.sd_position;
		
		if (sigma!=0.0) {
			
			for(d=0.0,i=0;i<mrkIdxL;i++)
				d+=MappingFunction.distance(mr[i], mfc, mut);
			
			for(w=0.0;i<mrkIdxR;i++) {
				h = 0.0;
				/* The probability of the current point */ 
				y1 = Math.exp(-(mu-d)*(mu-d)/(2.0*sigma*sigma))/(Math.sqrt(2.0*Math.PI)*sigma);
				if (y1 < Double.MIN_VALUE) y1 = 0.0;
				h -= d;
				d += MappingFunction.distance(mr[i], mfc, mut);
				h += d;
				y2 = Math.exp(-(mu-d)*(mu-d)/(2.0*sigma*sigma))/(Math.sqrt(2.0*Math.PI)*sigma);
				if (y2 < Double.MIN_VALUE) y2 = 0.0;
				/* Trapeze area */
				w += 0.5*(y1+y2)*h;				
			}
			
			if (w < Double.MIN_VALUE) {
				w = 0.0;
			}
			
			return w;
			
		} else if (mrkIdxL <= qtl.m1 && mrkIdxR >= qtl.m2)
			return 1.0;
		else
			return 0.0;
		
	}
	/**
	 * Returns true if the QTL position is inside the marker
	 * interval which markers indices are <code>mrkIdxL</code>
	 * and <code>mrkIdxR</code>.
	 * This assumes that <code>mrkIdxL &lt; mrkIdxR</code> 
	 * 
	 * @param mrkIdxL the left marker indice.
	 * @param mrkIdxR the right marker indice.
	 * @param qtlIdx the qtl indice.
	 * @return true if the QTL position is inside the marker interval.
	 */
	public boolean areMarkersAroundQTL(int mrkIdxL, int mrkIdxR, int qtlIdx) {
		return (mrkIdxL <= qtls[qtlIdx].m1 && mrkIdxR >= qtls[qtlIdx].m2);
	}
	/**
	 * This method tests if the chromosome have a 
	 * set of QTL defined.
	 * 
	 * @return true if there are QTL.
	 */
	public boolean hasQTL() {
		return (qtls != null && qtls.length > 0);
	}
	
	public String toString() {
		int i;
		double d;
		StringBuffer buff = new StringBuffer();
		
		buff.append(this.mapName);
		buff.append(System.getProperty("line.separator"));
		
		for(d=0.0,i=0;i<this.nm-1;i++) {
			
			buff.append(this.mrkNames[i]+" "+d);
			buff.append(System.getProperty("line.separator"));
			d+=MappingFunction.distance(this.mr[i],IMetaQtlConstants.MAPPING_FUNCTION_HALDANE,IMetaQtlConstants.MAPPING_UNIT_CM);
			
		}
		buff.append(this.mrkNames[i]+" "+d);
		buff.append(System.getProperty("line.separator"));
		
		return buff.toString();
	}
	/**
	 * This method looks if it exists a marker on the chromosome
	 * with the given name and returns its indice. Otherwise it
	 * returns -1.
	 * 
	 * @param string the name of the marker.
	 * @return the indice of the marker if exists, -1 otherwise.
	 */
	public int getMarkerIdxWithName(String string) {
		int i;
		
		if (mrkNames==null) return -1;
		
		for(i=0;i<mrkNames.length;i++)
			if (mrkNames[i].equals(string))
				return i;
		
		return -1;
	}
	/**
	 * This methods returns the distance between the markers with
	 * the given indices assuming that <code>mrkIdxL &lt; mrkIdxR</code>.
	 * The distance is computed using the given mapping context.
	 * 
	 * @param mrkIdxL the left marker indice.
	 * @param mrkIdxR the right marker indice.
	 * @param mapping_function the mapping function.
	 * @param mapping_unit the mapping unit.
	 * @return the distance between the two markers. 
	 */
	public double getDistanceBetween(int mrkIdxL, int mrkIdxR, int mapping_function, int mapping_unit) {
		int l;
		double d;
		
		if (mrkIdxL > mrkIdxR) return 0.0;
		
		for(d=0.0,l=mrkIdxL;l<mrkIdxR;l++)
			d+=MappingFunction.distance(mr[l], mapping_function, mapping_unit);
		
		return d;
	}
	/**
	 * 
	 * @param mrkIdxL
	 * @param mrkIdxR
	 * @return
	 */
	public double getDistanceBetween(int mrkIdxL, int mrkIdxR) {
		return getDistanceBetween(mrkIdxL, mrkIdxR, mfc, mut);
	}
	/**
	 * This methods returns the distance between the first marker and
	 * the marker with the given indice using the given mapping context. 
	 * 
	 * @param mrkIdx the marker indice.
	 * @param mapping_function the mapping function.
	 * @param mapping_unit the mapping unit.
	 * @return the distance between the first marker and the given one.
	 */
	public double getDistance(int mrkIdx, int mapping_function, int mapping_unit) {
		int l;
		double d;
		
		if (mrkIdx < 0) return 0.0;
		if (mrkIdx >= nm) return 0.0;
		
		if (mrkIdx == nm-1) {
			d = MappingFunction.recombination(totd, this.mfc, this.mut);
			d = MappingFunction.distance(d, mapping_function,mapping_unit);
		} else {
			for(d=0.0,l=0;l<mrkIdx;l++)
				d+=MappingFunction.distance(mr[l],mapping_function,mapping_unit);
		}
		
		return d;
	}
	/**
	 * Updates the mapping context of the chromosome. 
	 * 
	 * @param mapping_function the new mapping function.
	 * @param mapping_unit the new mapping unit.
	 */
	public void updateMappingContext(int mapping_function, int mapping_unit) {
		int j;
		boolean update;
		
		update = ((mapping_unit != mut) && (mapping_function != mfc));
		
		if (update) {
			if (qtls != null) {
				for(j=0;j<qtls.length;j++) {
					qtls[j].position=MappingFunction.recombination(qtls[j].position,this.mfc, this.mut);
					qtls[j].position=MappingFunction.distance(qtls[j].position,mapping_function,mapping_unit);
					if (qtls[j].ci != null) {
						qtls[j].ci.from = qtls[j].position - MappingFunction.distance(qtls[j].ci.ic1,mapping_function,mapping_unit);
						qtls[j].ci.to   = qtls[j].position + MappingFunction.distance(qtls[j].ci.ic2,mapping_function,mapping_unit);
					}
				}
			}
			this.mfc = mapping_function;
			this.mut     = mapping_unit;
		}
	}
	/**
	 * This methods returns the distance between the first marker and
	 * the marker with the given indice using the default mapping
	 * context. 
	 * @param i the indice of the marker.
	 * @return the distance betwee the first marker and the given one. 
	 */
	public double getDistance(int i) {
		return getDistance(i, mfc, mut);
	}
	
	public void r2R() {
		r2R(this.ct, this.cg);
	}
	/**
	 * Re-scale the chromosome assuming that the distances
	 * and recombination fractions are based on the frequency of
	 * recombinant genotypes R instead of the real recombination
	 * fraction r.
	 */
	public void r2R(int type, int gen) {
		int i;
		double tmp;
		
		totd = 0.0;
		
		for(i=0;i<this.nm-1;i++) {
			mr[i] = OffspringManager.r2R(mr[i], type, gen);
			totd += MappingFunction.distance(mr[i], mfc, mut);
		}
		
		if (qtls != null) {
			/*
			 * We also have to re-scale the QTL positions
			 * and CI absolute boundaries.
			 */
			for(i=0;i<qtls.length;i++) {
				qtls[i].c1 = OffspringManager.r2R(qtls[i].c1, type, gen);
				qtls[i].c2 = OffspringManager.r2R(qtls[i].c2, type, gen);
				tmp =  MappingFunction.distance(qtls[i].c1, mfc, mut);
				tmp += getDistance(qtls[i].m1, mfc, mut);
				qtls[i].position=tmp;
				if (qtls[i].ci != null) {
					//System.out.println("<QTL "+ i + " " + (qtls[i].ci.to-qtls[i].ci.from));
					qtls[i].ci.ic1  = OffspringManager.r2R(qtls[i].ci.ic1, type, gen);
					qtls[i].ci.ic2  = OffspringManager.r2R(qtls[i].ci.ic2, type, gen);
					qtls[i].ci.from = qtls[i].position - MappingFunction.distance(qtls[i].ci.ic1, mfc, mut);
					qtls[i].ci.to   = qtls[i].position + MappingFunction.distance(qtls[i].ci.ic2, mfc, mut);
					qtls[i].ci.from = Math.max(0.0, qtls[i].ci.from);
					qtls[i].ci.to   = Math.min(qtls[i].ci.to, totd);
					//System.out.println(">QTL "+ i + " " + (qtls[i].ci.to-qtls[i].ci.from));
				}
			}
		}
	}
	
	public void R2r() {
		R2r(this.ct, this.cg);
	}
	/**
	 * Re-scale the chromosome assuming that the distances
	 * and recombination fractions are based on the frequency of
	 * recombinant genotypes R instead of the real recombination
	 * fraction r.
	 */
	public void R2r(int type, int gen) {
		int i;
		double tmp;
		
		totd = 0.0;
		for(i=0;i<this.nm-1;i++) {
			//System.out.print(mr[i]+" "+cs+" "+ct+" "+cg+" ==> ");
			//mr[i] = OffspringManager.r2R(mr[i], cs, IBioCrossTypes.RILSe, 1);
			mr[i] = OffspringManager.R2r(mr[i], type, gen);
			//System.out.println(mr[i]);
			totd += MappingFunction.distance(mr[i], mfc, mut);
		}
		//System.out.println(name+" TOTD "+totd);
		
		if (qtls != null) {
			/*
			 * We also have to re-scale the QTL positions
			 * and CI absolute boundaries.
			 */
			for(i=0;i<qtls.length;i++) {
				qtls[i].c1 = OffspringManager.R2r(qtls[i].c1, type, gen);
				qtls[i].c2 = OffspringManager.R2r(qtls[i].c2, type, gen);
				tmp =  MappingFunction.distance(qtls[i].c1, mfc, mut);
				tmp += getDistance(qtls[i].m1, mfc, mut);
				qtls[i].position=tmp;
				if (qtls[i].ci != null) {
					//System.out.println("<QTL "+ i + " " + (qtls[i].ci.to-qtls[i].ci.from));
					qtls[i].ci.ic1  = OffspringManager.R2r(qtls[i].ci.ic1, type, gen);
					qtls[i].ci.ic2  = OffspringManager.R2r(qtls[i].ci.ic2, type, gen);
					qtls[i].ci.from = qtls[i].position - MappingFunction.distance(qtls[i].ci.ic1, mfc, mut);
					qtls[i].ci.to   = qtls[i].position + MappingFunction.distance(qtls[i].ci.ic2, mfc, mut);
					qtls[i].ci.from = Math.max(0.0, qtls[i].ci.from);
					qtls[i].ci.to   = Math.min(qtls[i].ci.to, totd);
					//System.out.println(">QTL "+ i + " " + (qtls[i].ci.to-qtls[i].ci.from));
				}
			}
		}
	}
	/**
	 * According to the given standard deviation computation
	 * mode this method computes for each QTL the value of the standard
	 * deviation of the QTL position. If no standard deviation can be
	 * computed then the standard deviation is set to 0.0;  
	 * 
	 * @see IMetaQtlConstants
	 * 
	 * @param mode the standard deviation computation mode.
	 */
	public void computeQtlSD(int mode) {
		
		if (qtls==null) return;
		
		int i;
		double[] sdci=null,sdr2=null;
		boolean max = (mode == IMetaQtlConstants.CI_MODE_MAXIMAL);
		boolean ci  = (mode == IMetaQtlConstants.CI_MODE_LOD || mode == IMetaQtlConstants.CI_MODE_AVAILABLE);
		boolean r2  = (mode == IMetaQtlConstants.CI_MODE_RSQUARE || mode == IMetaQtlConstants.CI_MODE_AVAILABLE);
		ci |= max;
		
		for(i=0;i<qtls.length;i++) 
			qtls[i].sd_position = 0.0;
		
		if (max) {
			sdci = new double[qtls.length];
			sdr2 = new double[qtls.length];
		}
		
		if (ci) {
			for(i=0;i<qtls.length;i++) {
				if (qtls[i].ci != null) {
					qtls[i].sd_position = Qtl.getQTLSD(qtls[i].ci, mfc, mut);
					if (max) {
						sdci[i]=qtls[i].sd_position;
					}
				}
			}
		}
		
		r2    |= max;
		
		if (r2){
			for(i=0;i<qtls.length;i++) {
				if (qtls[i].stats != null) {
					if (mode == IMetaQtlConstants.CI_MODE_AVAILABLE && qtls[i].sd_position!=0.0) continue;
					qtls[i].sd_position = Qtl.getQTLSD(qtls[i].stats.rsquare, qtls[i].cross, mut);
					//System.out.println(qtls[i].name+" "+qtls[i].stats.rsquare+" "+qtls[i].sd_position);
					// If the qtl has been projected then re-scale the standard-deviation
					// by the projection scale parameter.
					if (qtls[i].proj != null) {
						
						qtls[i].sd_position *= qtls[i].proj.mapScale;
					}
					if (max) {
						sdr2[i]=qtls[i].sd_position;
					}
				}
			}
		}
		if (max) {
			for(i=0;i<qtls.length;i++) {
				qtls[i].sd_position=Math.max(sdci[i],sdr2[i]);
			}
		}
	}
	/**
	 * 
	 * @return
	 */
	public Qtl[] getMetaQtls() {
		return mqtls;
	}
	/**
	 * 
	 * @param metaQtl
	 */
	public void setMetaQtls(Qtl[] metaQtl) {
		mqtls = metaQtl;
	}
	/**
	 * This method returns a copy of the QTL located
	 * on thic chromosome.
	 * 
	 * @return a copy of the array of QTLlocated on this chromosome.
	 */
	public Qtl[] getQtls() {
		if (qtls!=null) {
			Qtl[] copy = new Qtl[qtls.length];
			for(int i=0;i<copy.length;i++) {
				copy[i] = qtls[i].getCopy();
			}
			return copy;
		}
		return null;
	}
	/**
	 * @return
	 */
	public int getQtlNumber() {
		return (qtls==null) ? 0 : qtls.length;
	}
	/**
	 * @param mqtls2
	 */
	public void attachQtl(Qtl[] qtls) {
		int i,j,s;
	 	double d,dd,qtlpos;
		
		if (qtls == null) return;
		s=0;
		if (this.qtls != null) {
	 		Qtl[] tmp = new Qtl[qtls.length + this.qtls.length];
	 		for(i=0;i<this.qtls.length;i++)
	 			tmp[i]=this.qtls[i];
	 		s = this.qtls.length; 
	 		this.qtls = tmp;
	 	} else {
	 		this.qtls = new Qtl[qtls.length];
	 	}
	 	
	 	for(i=s;i<this.qtls.length;i++) {
	 		this.qtls[i] = qtls[i-s].getCopy();
	 		/* Then as the QTL position on the chromosome is
	 		 * absolute we have to find the marker interval
	 		 * in which the QTL belongsand computes the left
	 		 * and rigth distance of the QTL from the flanking
	 		 * markers.
	 		 */
	 		qtlpos = this.qtls[i].position;
	 		
	 		d = 0.0;
	 		
	 		for(j=0;j<nm-1;j++) {
	 			
	 			d += MappingFunction.distance(mr[j], mfc, mut);
	 			
	 			if (d >= qtlpos) {
	 				/* Get the flanking marker indices and the recombination
	 				 * rates between the QTL position and them.
	 				 */ 
	 				this.qtls[i].m1       = j;
	 				this.qtls[i].m2       = j+1;
	 				this.qtls[i].c2       = d-qtlpos;
	 				dd	    			  = MappingFunction.distance(mr[j], mfc, mut);
	 				this.qtls[i].c1       = dd-this.qtls[i].c2;
	 				this.qtls[i].c1       = MappingFunction.recombination(this.qtls[i].c1, mfc, mut);
	 				this.qtls[i].c2       = MappingFunction.recombination(this.qtls[i].c2, mfc, mut);
	 				break;
	 			}
	 			
	 		} /* end loop on markers */
	 		
	 		/* Now, is there a CI defined for this QTL ? */
	 		if (this.qtls[i].ci != null) {
	 			/* Convert the distance in recombination rate */
	 			this.qtls[i].ci.ic1 = qtlpos - this.qtls[i].ci.from;
	 			if (this.qtls[i].ci.ic1 < 0) {
	 				System.err.println("[ WARNING ]");
	 				System.err.print("  In map " + mapName + " on chrom "+name+" : ");
	 				System.err.println("  BAD LOD SUPPORT FOR QTL " + this.qtls[i].name);
	 				this.qtls[i].ci = null;
	 			} else {
		 			this.qtls[i].ci.ic1 = MappingFunction.recombination(this.qtls[i].ci.ic1, mfc, mut);
		 			if (this.qtls[i].ci.to == 0.0) {
		 				/* This means that any upper bound has been defined
		 				 * then by default use the end of the chromosome as
		 				 * lower bound
		 				 */
		 				this.qtls[i].ci.to = totd;
		 				System.err.println("[ WARNING ]");
		 				System.err.print("  In map " + mapName + " on chrom "+name+" : ");
		 				System.err.println("  NO UPPER LOD SUPPORT BOUNDARY FOR QTL " + this.qtls[i].name);
		 			}
		 		    this.qtls[i].ci.ic2 = this.qtls[i].ci.to - qtlpos;
		 		    if (this.qtls[i].ci.ic2 < 0) {
		 				System.err.println("[ WARNING ]");
		 				System.err.print("  In map " + mapName + " on chrom "+name+" : ");
		 				System.err.println("  BAD LOD SUPPORT FOR QTL " + this.qtls[i].name);
		 				this.qtls[i].ci = null;
		 		    } else {
		 		    	this.qtls[i].ci.ic2 = MappingFunction.recombination(this.qtls[i].ci.ic2, mfc, mut);
		 		    }
	 			}
	 		}
	 		
	 		if (this.qtls[i].cross == null) {
	 			
	 			this.qtls[i].cross = this.qtls[i].new QTLCross();
	 			this.qtls[i].cross.gen  = this.cg;
	 			this.qtls[i].cross.type = this.ct;
	 			this.qtls[i].cross.size = this.cs;
	 			
	 		}
	 		
	 	} /* end loop on QTL */
		
	}
	
	/**
	 * Add  
	 * @param loci
	 */
	 public void attachQTL(IBioLocus[] qtls) {
	 	int i;
	 	
	 	if (qtls == null) return;
	 	
	 	Qtl[] tmp = new Qtl[qtls.length]; 
	 	for(i=0;i<qtls.length;i++) {
	 		tmp[i] = QtlAdapter.toQTL(qtls[i]);
	 	}
	 	
	 	attachQtl(tmp);
	 	
	 }
	/**
	 * Returns the names of the marker at the i^th position
	 * on the chromosome.
	 * @param i the position of the marker.
	 * @return the names of the marker.
	 */
	public String getMarkerName(int i) {
		if (i>=0 && i<mrkNames.length)
			return mrkNames[i];
		else
			return null;
	}
	/**
	 * Returns the name of the chromosome.
	 * 
	 * @return the name of teh chromosome.
	 */
	public String getName() {
		return name;
	}
	/**
	 * 
	 * @param ncm
	 * @param mrkIdx
	 * @param chrom
	 * @param ref
	 * @return
	 */
	public static double[] getRecHeterogeneity(int ncm, int[][] mrkIdx, Chromosome chrom, Chromosome ref)
	{
		if (ncm <= 1) return null;
		
		int mc1,mc2,mr1,mr2,mf,mu,df;
		double d1,d2,r,r1,r2,v1,v2;
		double[] residuals = new double[ncm-1];
		
		mf=IMetaQtlConstants.MAPPING_FUNCTION_HALDANE;
		mu=IMetaQtlConstants.MAPPING_UNIT_CM;
		
		for(int i=0;i<ncm-1;i++) {
			mc1=mrkIdx[0][i];
			mc2=mrkIdx[0][i+1];
			mr1=mrkIdx[1][i];
			mr2=mrkIdx[1][i+1];
			df  = 0;
			d1  = chrom.getDistanceBetween(mc1,mc2,mf,mu);
			r1  = MappingFunction.recombination(d1,mf,mu);
			v1  = 0.0;
			if (chrom.ct!=IBioCrossTypes.UNKNOWN) {
				v1 = OffspringManager.fisherVariance(r1,chrom.cs,chrom.ct,chrom.cg);
				df++;
			}
			d2  = ref.getDistanceBetween(mr1,mr2,mf,mu);
			r2  = MappingFunction.recombination(d2,mf,mu);
			v2  = 0.0;
			if (ref.ct!=IBioCrossTypes.UNKNOWN) {
				v2 = OffspringManager.fisherVariance(r2,ref.cs,ref.ct,ref.cg);
				df++;
			}
			if (df==1) {
				residuals[i]=(v1>0.0) ? (r1-r2)*(r1-r2)/v1 : (r1-r2)*(r1-r2)/v2;
				residuals[i]=NumericalUtilities.gammp(0.5,0.5*residuals[i]);
			} else if (df==2) {
				r=(r1/v1+r2/v2)/(1.0/v1+1.0/v2);
				residuals[i]=((r1-r)/v1+(r2-r)/v2)/(1.0/v1+1.0/v2);
				residuals[i]=NumericalUtilities.gammp(0.5,0.5*residuals[i]);
			} else residuals[i]=0.0;
		}
		
		return residuals;
	}
	
	public static double[] getDistHeterogeneity(int ncm, int[][] mrkIdx, Chromosome chrom, Chromosome ref)
	{
		if (ncm <= 1) return null;
		
		int mc1,mc2,mr1,mr2,mf,mu,df;
		double d1,d2,r,r1,r2,v1,v2;
		double[] residuals = new double[ncm-1];
		
		mf=IMetaQtlConstants.MAPPING_FUNCTION_HALDANE;
		mu=IMetaQtlConstants.MAPPING_UNIT_CM;
		
		for(int i=0;i<ncm-1;i++) {
			mc1=mrkIdx[0][i];
			mc2=mrkIdx[0][i+1];
			mr1=mrkIdx[1][i];
			mr2=mrkIdx[1][i+1];
			df  = 0;
			d1  = chrom.getDistanceBetween(mc1,mc2,mf,mu);
			r1  = MappingFunction.recombination(d1,mf,mu);
			v1  = 0.0;
			if (chrom.ct!=IBioCrossTypes.UNKNOWN) {
				v1 = OffspringManager.fisherVariance(r1,chrom.cs,chrom.ct,chrom.cg);
				v1 = MappingFunction.varianceDistance(r1, v1, mf, mu);
				df++;
			}
			d2  = ref.getDistanceBetween(mr1,mr2,mf,mu);
			r2  = MappingFunction.recombination(d2,mf,mu);
			v2  = 0.0;
			if (ref.ct!=IBioCrossTypes.UNKNOWN) {
				v2 = OffspringManager.fisherVariance(r2,ref.cs,ref.ct,ref.cg);
				v2 = MappingFunction.varianceDistance(r2, v2, mf, mu);
				df++;
			}
			if (df==1 && (v1 != 0.0 || v2 != 0.0)) {
				residuals[i]=(v1>0.0) ? (d1-d2)/Math.sqrt(v1) : (d1-d2)/Math.sqrt(v2);
				residuals[i]*=residuals[i];
				residuals[i]=NumericalUtilities.gammp(0.5,0.5*residuals[i]);
			} else if (df==2 && (v1 != 0.0 && v2 != 0.0)) {
				r=(d1/v1+d2/v2)/(1.0/v1+1.0/v2);
				residuals[i]=((d1-r)/v1+(d2-r)/v2)/(1.0/v1+1.0/v2);
				residuals[i]=NumericalUtilities.gammp(0.5,0.5*residuals[i]);
			} else residuals[i]=0.0;
		}
		
		return residuals;
	}
	
	/**
	 * @return
	 */
	public double[][] getProba() {
		return proba;
	}
	/**
	 * @param proba
	 */
	public void setProba(double[][] proba) {
		this.proba  = proba;
	}
	/**
	 * @param chromosome
	 * @param chromosome2
	 * @return
	 */
	public static int getFirstCommonMrkIdx(Chromosome chrom, Chromosome chrom2) {
		for(int i=0;i<chrom2.nm;i++) {
			for(int j=0;j<chrom.nm;j++) {
				if (chrom.mrkNames[j].equals(chrom2.mrkNames[i]))
					return j;
			}
		}
		return -1;
	}
	/**
	 * @param strings
	 * @return
	 */
	public Chromosome getSubChromosome(String[] markerNames) {
		
		if (markerNames == null) return null;
		
		Chromosome sub = null;
		
		ArrayList  mrk = new ArrayList();
		
		for(int i=0;i<this.nm;i++) {
			for(int j=0;j<markerNames.length;j++) {
				if (this.mrkNames[i].equals(markerNames[j])) {
					mrk.add(this.mrkNames[i]);
				}
			}
		}
		
		if (mrk.size() > 1) {
			
			sub 		   = new Chromosome();
			sub.nm         = mrk.size();
			sub.cg         = this.cg;
			sub.cs         = this.cs;
			sub.ct         = this.ct;
			sub.mfc        = this.mfc;
			sub.mut        = this.mut;
			sub.skeleton   = this.skeleton;
			sub.mapName    = new String(this.mapName);
			sub.name       = new String(this.name);
			
			sub.mrkNames   = new String[mrk.size()];
			mrk.toArray(sub.mrkNames);
			sub.mr         = new double[mrk.size()-1];
			sub.totd       = 0.0;
			for(int i=0;i<sub.nm-1;i++) {
				//System.out.println(" + " + sub.mrkNames[i]);
				int m1=this.getMarkerIdxWithName(sub.mrkNames[i]);
				int m2=this.getMarkerIdxWithName(sub.mrkNames[i+1]);
				sub.mr[i]=this.getDistanceBetween(m1, m2, sub.mfc, sub.mut);
				sub.totd+=sub.mr[i];
				sub.mr[i]=MappingFunction.recombination(sub.mr[i], sub.mfc, sub.mut);
			}
			if (qtls!=null) {
				ArrayList qtl2Attach = new ArrayList();
				for(int i=0;i<this.qtls.length;i++) {
					if (qtls[i].m1 >= getMarkerIdxWithName(sub.mrkNames[0])
					    && qtls[i].m2 <= getMarkerIdxWithName(sub.mrkNames[sub.nm-1])) {
						Qtl qtl      = qtls[i].getCopy();
						updateQtlPos(qtl, getMarkerIdxWithName(sub.mrkNames[0]), sub.totd);
						qtl2Attach.add(qtl);
					}
				}
				if (qtl2Attach.size()>0) {
					sub.qtls = new Qtl[qtl2Attach.size()];
					qtl2Attach.toArray(sub.qtls);
				}
			}
		}
		
		return sub;
	}
	/**
	 * @param qtl
	 * @param markerIdxWithName
	 */
	private void updateQtlPos(Qtl qtl, int midx, double length) {
		double shift  = this.getDistance(midx);
		qtl.position -= shift;
		if (qtl.ci != null) {
			qtl.ci.from = qtl.position - MappingFunction.distance(qtl.ci.ic1, mfc, mut);
			qtl.ci.from = (qtl.ci.from >= 0.0) ? qtl.ci.from : 0.0;
			qtl.ci.to   = qtl.position + MappingFunction.distance(qtl.ci.ic2, mfc, mut);
			qtl.ci.to   = (qtl.ci.to >= length) ? qtl.ci.to : length;
		}
	}
	/**
	 * @return
	 */
	public String getMapName() {
		return mapName;
	}
}