/*  
 *  src/org/metaqtl/util/OffspringManager.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.util;


import org.metaqtl.bio.IBioCrossTypes;
import org.metaqtl.numrec.NumericalUtilities;

/**
 * 
 * Class Description Here
 * 
 * @author Jean-Baptiste Veyrieras
 *
 */
public final class OffspringManager {
	
	public OffspringManager() {}

	/**
	 * Returns the inverse of the Fisher Information depending on
	 * the type of cross for a given recombination rate and a given
	 * population size.
	 * @param r the recombination rate estimate.
	 * @param size
	 * @param type
	 * @return the variance for the recombination rate estimate.
	 */
	public static double fisherVariance(double r, int size, int type, int gen) {
		int g;
		double R,tmp,var=0.0;
		
		if (type == IBioCrossTypes.BC) {
			var = r * (1.0 - r) / (double)size;
		} else if (type == IBioCrossTypes.RILSi) {
			var = r*(1.0+6.0*r)*(1.0+6.0*r)*(1.0+2.0*r)/(4.0*(double)size);
		} else if (type == IBioCrossTypes.RILSe) {
			var = r*(1.0+2.0*r)*(1.0+2.0*r)/(2.0*(double)size);
		} else if (type == IBioCrossTypes.SF) {
			if (gen == 2) { /* F2 */
				var = r * (1.0 - r) * (1.0 - 2.0 * r + 2.0 * r * r);
				var /= (2.0 * (double)size * (1.0 - 3.0 * r + 3.0 * r * r));			
			} else if (gen > 2) {
				/* Use the Haldane - Waddington recurrence formula */
				HaldaneWaddingtonProportion p = new OffspringManager().new HaldaneWaddingtonProportion();
				g=1;
				do {
					selfingRecurrence(p, r);
					g++;
				} while(g < gen);
				var = selfingIndividualInformation(p);
				var = 1.0/(var*(double)size);
			}
		} else if (type == IBioCrossTypes.RF) {
			if (gen == 2) { /* F2 */
				var = r * (1.0 - r) * (1.0 - 2.0 * r + 2.0 * r * r);
				var /= (2.0 * (double)size * (1.0 - 3.0 * r + 3.0 * r * r));				
			} else if (gen > 2) {
				/*
				 * See Liu, S.C. et al., Genome-Wide High-Resolution Mapping by
				 * Recurrent Intermating Using Arabidopsis thaliana as a Model,
				 * 1996, Genetics 142:247-258. 
				 */
				tmp  = Math.pow(1.0-r, 4*(gen-1));
				var  = (1.0-Math.pow(1-2.0*r, 4)*tmp);
				var  /= Math.pow(1.0-r, 2*(gen-1)-2);
				var  /= Math.pow(2.0*(1.0-r)+(gen-1)*(1.0-2.0*r), 2);
				tmp   = Math.pow(1.0-r, 2*(gen-1));
				var  /= (1.0+3.0*(1-2.0*r)*(1-2.0*r)*tmp);
				var  /= (double)size;
			}
		} else if (type == IBioCrossTypes.IRISe) {
			IRISEFunction func1 = new OffspringManager().new IRISEFunction();
			/* Set the number of random mating generations */
			func1.t = gen;
			/* Get the propertion of recombinant genotypes from r */
			R    = r2R(r, type, gen);
			var  = fisherVariance(R, size, IBioCrossTypes.BC, gen);
			/* Compute the derivate of r = f(R) in R */
			tmp   = Math.min(R, 0.5-R)/10000.0;
			NumericalUtilities.dfridr(func1, R, tmp);
			tmp  = func1.getDF();
			/* Here we assume a first order approximation */
			var  = var*(tmp*tmp);
		} else if (type == IBioCrossTypes.IRISi) {
			IRISIFunction func1 = new OffspringManager().new IRISIFunction();
			/* Set the number of random mating generations */
			func1.t = gen;
			/* Get the propertion of recombinant genotypes from r */
			R    = r2R(r, type, gen);
			var  = fisherVariance(R, size, IBioCrossTypes.BC, gen);
			/* Compute the derivate of r = f(R) in R */
			tmp   = Math.min(R, 0.5-R)/10000.0;
			NumericalUtilities.dfridr(func1, R, tmp);
			tmp  = func1.getDF();
			/* Here we assume a first order approximation */
			var  = var*(tmp*tmp); 
		} else if (type == IBioCrossTypes.UNKNOWN) {
			var = fisherVariance(r, size, IBioCrossTypes.BC, gen);
		}
		return var;
	}
	/**
	 * Given the current proportion C,D,E,F and G and the recombination
	 * rate r this routine computes the next generation proportions following 
	 * the recurrence formula of Haldane and Waddington,
	 * 1931, Inbreeding and linkage, Genetics 16:357-374
	 * (see section self-fertilization p358).
	 * 
	 * Cn+1 = Cn + 0.5*En + 0.25*(1-r)^2*Fn + 0.25*r^2*Gn
	 * Dn+1 = p.D + 0.5*En + 0.25*r^2Fn + 0.25*(1-r)^2Gn
	 * En+1 = 0.5*En + 0.5*r*(1-r)*(Fn + Gn)
	 * Fn+1 = 0.5*(1-r)^2Fn + 0.5*r^2Gn
	 * Gn+1 = 0.5*r^2Fn + 0.5*(1-r)^2*Gn
	 * 
	 * This also computes the next generation derivates dC,dD,dE,dF and dG which can
	 * be easily obtained from the above equations.
	 *  
	 * @param proportion
	 */
	public static void selfingRecurrence(HaldaneWaddingtonProportion  p, double r) {
		 double C,D,E,F,G,dC,dD,dE,dF,dG;
		
		 C = p.C + 0.5*p.E + 0.25*(1.0-r)*(1.0-r)*p.F + 0.25*r*r*p.G;
		 D = p.D + 0.5*p.E + 0.25*r*r*p.F + 0.25*(1.0-r)*(1.0-r)*p.G;
		 E = 0.5*p.E + 0.5*r*(1.0-r)*(p.F + p.G);
		 F = 0.5*(1.0-r)*(1.0-r)*p.F + 0.5*r*r*p.G;
		 G = 0.5*r*r*p.F + 0.5*(1.0-r)*(1.0-r)*p.G;
		 
		 dC = p.dC + 0.5*p.dE - 0.5*(1.0-r)*p.F + 0.25*(1.0-r)*(1.0-r)*p.dF + 0.5*r*p.G + 0.25*r*r*p.dG;
		 dD = p.dD + 0.5*p.dE + 0.5*r*p.F + 0.25*r*r*p.dF - 0.5*(1.0-r)*p.G + 0.25*(1.0-r)*(1.0-r)*p.dG ;
		 dE = 0.5*p.dE + 0.5*(1.0-2.0*r)*(p.F + p.G) + 0.5*r*(1.0-r)*(p.dF + p.dG);
		 dF = (r-1.0)*p.F + 0.5*(1.0-r)*(1.0-r)*p.dF + r*p.G + 0.5*r*r*p.dG;
		 dG = (r-1.0)*p.G + 0.5*(1.0-r)*(1.0-r)*p.dG + r*p.F + 0.5*r*r*p.dF;
		 
		 p.C = C;
		 p.D = D;
		 p.E = E;
		 p.F = F;
		 p.G = G;
		 p.dC = dC;
		 p.dD = dD;
		 p.dE = dE;
		 p.dF = dF;
		 p.dG = dG;		 
	}
	/**
	 * This routine returns the individual means amount of information
	 * for the given generation propertions for a self-fertilization
	 * design. This assumes that the 4 zygotic classes can be discriminated
	 * without uncertainty - i.e the case where markers are codominants.
	 * @param p
	 */
	public static double selfingIndividualInformation(HaldaneWaddingtonProportion  p) {
		double m1,m2,m3,m4,dm1,dm2,dm3,dm4,ir=0.0;
		
		m1 = p.C/HaldaneWaddingtonProportion.C_coeff;
		m2 = p.D/HaldaneWaddingtonProportion.D_coeff;
		m3 = p.E/HaldaneWaddingtonProportion.E_coeff;
		m4 = p.F/HaldaneWaddingtonProportion.F_coeff;
		m4 += p.G/HaldaneWaddingtonProportion.G_coeff;
		dm1 = p.dC;
		if (m1==0.0) {dm1=0.0;m1=1.0;}
		dm2 = p.dD;
		if (m2==0.0) {dm2=0.0;m2=1.0;}
		dm3 = p.dE;
		if (m3==0.0) {dm3=0.0;m3=1.0;}
		dm4 = p.dF + p.dG;
		if (m4==0.0) {dm4=0.0;m4=1.0;}
		
		ir  = (dm1*dm1)/m1;
		ir += (dm2*dm2)/m2;
		ir += (dm3*dm3)/m3;
		ir += (dm4*dm4)/m4;
		
		return ir;
	}
	
	/**
	 * This routine returns the recombination rate r for a given 
	 * value of the propertion of recombibant genotypes R according
	 * to the cross design.
	 * 
	 * type :
	 * 	
	 * IBioCrossTypes.RILSe	Recombinant Inbred Lines via selfing
	 * 
	 * 			R = 2r/(1+2r) (Haldane and Waddington 1931)
	 * 
	 * 	IBioCrossTypes.RILSi Recombinant Inbred Lines via sib-mating
	 * 
	 * 			R = 4r/(1+6r) (Haldane and Waddington 1931)
	 * 	
	 * IBioCrossTypes.IRISe High Recombinant Inbred Lines via selfing
	 * 
	 * 			R = 0.5(1-((1-2r)/(1+2r))(1-r)^t) (Winkler et al. 2003)
	 * 
	 * IBioCrossTypes.IRISi High Recombinant Inbred Lines via sib-mating
	 * 
	 * 			R = 0.5(1-((1-2r)/(1+6r))(1-r)^t) (Winkler et al. 2003)
	 * 
	 * @param R the proportion of recombinant genotypes. 
	 * @param size the population size.
	 * @param type the type of pedigree.
	 * @param gen the number of random mating generations (only for IRISi and IRISe).
	 * @return the recombination rate r.
	 */
	public static double R2r(double R, int type, int gen) {
		double r = 0.0;
		
		if (R==0.0) return 0.0;
		if (R==0.5) return 0.5;
		
		switch(type) {
			case IBioCrossTypes.RILSe :
				r = R/(2.0*(1.0-R));
				break;
			case IBioCrossTypes.RILSi :
				r = R/(4.0-6.0*R);
				break;
			case IBioCrossTypes.RF :
				/* Following Liu et al. we find the root z in [1/2,1]
				 * of the following polynome :
				 * 
				 * 2z^{t+1} - z^t + (2R-1) = 0
				 * 
				 * where z = 1-r
				 * 
				 */
				IRFunction func0 = new OffspringManager().new IRFunction();
				func0.t = gen;
				r = func0.func(R);
				break;
			case IBioCrossTypes.IRISe :
				/* Following Winkler et al. we find the root z in [1/2,1]
				 * of the following polynome :
				 * 
				 * 2z^{t+1} - z^t + 2(1-2R)z + 3(2R-1) = 0
				 * 
				 * where z = 1-r
				 * 
				 */
				IRISEFunction func1 = new OffspringManager().new IRISEFunction();
				func1.t = gen;
				r = func1.func(R);
				break;
			case IBioCrossTypes.IRISi :
				/* Following Winkler et al. we find the root z in [1/2,1]
				 * of the following polynome :
				 * 
				 * 2z^{t+1} - z^t + 6(1-2R)z + 7(2R-1) = 0
				 * 
				 * where z = 1-r
				 * 
				 */
				IRISIFunction func2 = new OffspringManager().new IRISIFunction();
				func2.t = gen;
				r = func2.func(R);
				break;
			default :
				r = R;
				break;
		}
		
		return r;
	}
	
	/**
	 * This routine returns the recombination rate r for a given 
	 * value of the propertion of recombibant genotypes R according
	 * to the cross design.
	 * 
	 * type :
	 * 	
	 * IBioCrossTypes.RILSe	Recombinant Inbred Lines via selfing
	 * 
	 * 			R = 2r/(1+2r) (Haldane and Waddington 1931)
	 * 
	 * 	IBioCrossTypes.RILSi Recombinant Inbred Lines via sib-mating
	 * 
	 * 			R = 4r/(1+6r) (Haldane and Waddington 1931)
	 * 	
	 * IBioCrossTypes.IRISe High Recombinant Inbred Lines via selfing
	 * 
	 * 			R = 0.5(1-((1-2r)/(1+2r))(1-r)^t) (Winkler et al. 2003)
	 * 
	 * IBioCrossTypes.IRISi High Recombinant Inbred Lines via sib-mating
	 * 
	 * 			R = 0.5(1-((1-2r)/(1+6r))(1-r)^t) (Winkler et al. 2003)
	 * 
	 * @param R the proportion of recombinant genotypes. 
	 * @param size the population size.
	 * @param type the type of pedigree.
	 * @param gen the number of random mating generations (only for IRISi and IRISe).
	 * @return the recombination rate r.
	 */
	public static double r2R(double r, int type, int gen) {
		double R = 0.0;
		
		if (r==0.0) return 0.0;
		if (r==0.5) return 0.5;
		
		switch(type) {
			case IBioCrossTypes.RILSe :
				R = 2.0*r/(1.0+2.0*r);
				break;
			case IBioCrossTypes.RILSi :
				R = 4.0*r/(1.0+6.0*r);
				break;
			case IBioCrossTypes.RF :
				R = 0.5*(1-Math.pow(1.0-r, gen)*(1.0-2.0*r));
			case IBioCrossTypes.IRISe :
				R = 0.5*(1.0-(1.0-2.0*r)/(1.0+2.0*r)*Math.pow(1.0-r, gen)); 
				break;
			case IBioCrossTypes.IRISi :
				R = 0.5*(1.0-(1.0-2.0*r)/(1.0+6.0*r)*Math.pow(1.0-r, gen)); 
				break;
			default :
				R = r;
				break;
		}
		
		return R;
	}
	
	private class IRFunction implements NumericalUtilities.Function {
		
		int t;
		
		double df;
		
		IRPolynome func1 = new OffspringManager().new IRPolynome();
		
		/* (non-Javadoc)
		 * @see org.thalia.metaqtl.numrec.NumericalUtilities.Function#func(double)
		 */
		public double func(double x) {
			double r;
			func1.R = x;
			func1.t = t;
			r = NumericalUtilities.rtsafe(func1, 0.5, 1.0, 1.e-8);
			return  1.0-r;
		}

		/* (non-Javadoc)
		 * @see org.thalia.metaqtl.numrec.NumericalUtilities.Function#getDF()
		 */
		public double getDF() {
			return df;
		}

		/* (non-Javadoc)
		 * @see org.thalia.metaqtl.numrec.NumericalUtilities.Function#setDF(double)
		 */
		public void setDF(double df) {
			this.df = df;
		}
		
	
	}
	
	
	private class IRISEFunction implements NumericalUtilities.Function {
		
		int t;
		
		double df;
		
		IRISEPolynome func1 = new OffspringManager().new IRISEPolynome();
		
		/* (non-Javadoc)
		 * @see org.thalia.metaqtl.numrec.NumericalUtilities.Function#func(double)
		 */
		public double func(double x) {
			double r;
			func1.R = x;
			func1.t = t;
			r = NumericalUtilities.rtsafe(func1, 0.5, 1.0, 1.e-8);
			return  1.0-r;
		}

		/* (non-Javadoc)
		 * @see org.thalia.metaqtl.numrec.NumericalUtilities.Function#getDF()
		 */
		public double getDF() {
			return df;
		}

		/* (non-Javadoc)
		 * @see org.thalia.metaqtl.numrec.NumericalUtilities.Function#setDF(double)
		 */
		public void setDF(double df) {
			this.df = df;
		}
		
	
	}
	
	
	private class IRISIFunction implements NumericalUtilities.Function {
		
		int t;
		
		double df;
		
		IRISIPolynome func1 = new OffspringManager().new IRISIPolynome();
		
		/* (non-Javadoc)
		 * @see org.thalia.metaqtl.numrec.NumericalUtilities.Function#func(double)
		 */
		public double func(double x) {
			double r;
			func1.R = x;
			func1.t = t;
			r = NumericalUtilities.rtsafe(func1, 0.5, 1.0, 1.e-8);
			return  1.0-r;
		}

		/* (non-Javadoc)
		 * @see org.thalia.metaqtl.numrec.NumericalUtilities.Function#getDF()
		 */
		public double getDF() {
			return df;
		}

		/* (non-Javadoc)
		 * @see org.thalia.metaqtl.numrec.NumericalUtilities.Function#setDF(double)
		 */
		public void setDF(double df) {
			this.df = df;
		}
		
	
	}
	
	
	
	private class IRPolynome implements NumericalUtilities.NRFunction {
		
		public double f;
		
		public double df;
		
		public int t=0;
		
		public double R=0.0;
		
		/* (non-Javadoc)
		 * @see org.thalia.metaqtl.numrec.NumericalUtilities.NRFunction#funcd(double)
		 */
		public void funcd(double x) {
			/* Function value */
			f = 2.0*Math.pow(x, t+1) - Math.pow(x, t) + (2.0*R-1.0);
			/* Derivate value */
			df = 2.0*(t+1)*Math.pow(x, t) - t*Math.pow(x, t-1);
		}
		
		public double getF() {return f;}
		
		public double getDF() {return df;}
		
	}
	/**
	 * 
	 * @author jveyrier
	 *
	 *  
	 */
	private class IRISEPolynome implements NumericalUtilities.NRFunction {
		
		public double f;
		
		public double df;
		
		public int t=0;
		
		public double R=0.0;
		
		/* (non-Javadoc)
		 * @see org.thalia.metaqtl.numrec.NumericalUtilities.NRFunction#funcd(double)
		 */
		public void funcd(double x) {
			/* Function value */
			f = 2.0*Math.pow(x, t+1) - Math.pow(x, t) + (2.0-4.0*R)*x + 3.0*(2.0*R-1.0);
			/* Derivate value */
			df = 2.0*(t+1)*Math.pow(x, t) - t*Math.pow(x, t-1) + (2.0-4.0*R);
		}
		
		public double getF() {return f;}
		
		public double getDF() {return df;}
		
	}
	/**
	 * 
	 * @author jveyrier
	 *
	 *  
	 */
	private class IRISIPolynome implements NumericalUtilities.NRFunction {
		
		public double f;
		
		public double df;
		
		public int t=0;
		
		public double R=0.0;
		
		/* (non-Javadoc)
		 * @see org.thalia.metaqtl.numrec.NumericalUtilities.NRFunction#funcd(double)
		 */
		public void funcd(double x) {
			/* Function value */
			f = 2.0*Math.pow(x, t+1) - Math.pow(x, t) + 6.0*(1.0-2.0*R)*x + 7.0*(2.0*R-1.0);
			/* Derivate value */
			df = 2.0*(t+1)*Math.pow(x, t) - t*Math.pow(x, t-1) + 6.0*(1.0-2.0*R);
		}
		
		public double getF() {return f;}
		
		public double getDF() {return df;}
		
	}
	
	
	public class HaldaneWaddingtonProportion {
		
		public final static double C_coeff = 1.0;
		
		public final static double D_coeff = 1.0;
		
		public final static double E_coeff = 2.0;
		
		public final static double F_coeff = 0.5;
		
		public final static double G_coeff = 0.5;
		
		public double C;
		
		public double D;
		
		public double E;
		
		public double F;
		
		public double G;
		
		public double dC;
		
		public double dD;
		
		public double dE;
		
		public double dF;
		
		public double dG;
		
		public HaldaneWaddingtonProportion() {
			C=D=E=G=0.0;
			F=2.0;
			dC=dD=dE=dF=dG=0.0;
		}
		
	}
}
