/*  
 *  src/org/metaqtl/algo/MetaMapAlgorithm.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 org.metaqtl.ChromCluster;
import org.metaqtl.Chromosome;
import org.metaqtl.IMetaQtlConstants;
import org.metaqtl.MapCluster;
import org.metaqtl.MetaChrom;
import org.metaqtl.MetaDico;
import org.metaqtl.MetaMap;
import org.metaqtl.MultiFitResult;
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 is dedicated to the construction of a consensus
 * genetic map from a set of genetic maps based on a Weigthed
 * Least Square approach.  
 */
public class MetaMapAlgorithm extends MetaAlgorithm {
	/**
	 * The minimal number of times a marker must be seen in linkage groups. 
	 */
	private int mrkThresh = 1;
	/**
	 * The cluster of chromosomes on which the projection would be performed.
	 */
	private MapCluster mapCluster=null;
	/**
	 * 
	 */
	private MetaMap consMap;
	
	/**
	 * 
	 * @param maps
	 * @param skeleton
	 * @param dubious
	 * @param mrkDico
	 */
	public MetaMapAlgorithm(IBioGenome[] maps, IBioGenome skeleton, IBioLocus[] dubious, MetaDico mrkDico) {
		/*
		 * First, clusterize the linkage groups
		 */
		mapCluster = new MapCluster();
		mapCluster.setDubiousMarker(dubious);
		mapCluster.setMarkerDico(mrkDico);
		for(int i=0;i<maps.length;i++) {
			mapCluster.addMap(maps[i], false);
		}
		if (skeleton!=null) mapCluster.addMap(skeleton, true);
	}
	/**
	 * 
	 * @return
	 */
	public MetaMap getResult() {
		return consMap;
	}
	/**
	 * This method builds a consensus genetic map from the given
	 * array of <code>maps</code>. The consensus map is returned
	 * as a {@link MetaMap}. If <code>skeleton</code> is not null
	 * then the positions of the markers on this map are assumed to
	 * be known and are used as anchor points in the construction
	 * of the consensus map. If <code>dubious<code> is no null or not
	 * empty then the loci of this array won't be included in the consensus
	 * map.
	 *   
	 * @param maps the array of genetic maps.
	 * @param skeleton a genetic map used as skeleton.
	 * @param dubious an array of dubious loci to not include in the construction.
	 * @return the consensus genetic map.
	 */
	public void run() {
		 /*
		 * Then get the chromosme clusters
		 */
		ChromCluster[] chrClusters = mapCluster.getClusters();
		
		if (chrClusters != null) {
			
			consMap = new MetaMap(chrClusters.length);
			
			for(int i=0;i<chrClusters.length;i++) {
				
				consMap.chromosomes[i] = buildMetaChrom(chrClusters[i]);
				
			}
		}
	}
	/**
	 * For a cluster of chromosomes <code>cluster</code> this method
	 * builds the consensus chromosome using a Weighted Least Square
	 * strategy and returns the result as a {@link MetaChrom} object.
	 * 
	 * @param cluster the cluster of chromosomes.
	 * @return the consensus chromosome.
	 */
	protected MetaChrom buildMetaChrom(ChromCluster cluster) {
		int i,j,k,l,xi,xj1,xj2,nc,nm,ni,size,type,gen;
		double r,d,w,x1,x2;
		double[]   distances;
		double[]   weights;
		double[][] X;
		Chromosome[] chromosomes;
		MultiFitResult result;		
		
		if (cluster==null) return null;
		
		cluster.setMrkThresh(mrkThresh);
		/* Then fixe the locus name index */
		cluster.fixCluster();
		/*double[] ratios = cluster.getCMrkTotalNumbers();
		for(v=0.0,i=0;i<ratios.length;i++)
			v+=ratios[i];
		nc=cluster.getMarkerNumber();
		System.out.println(cluster.name + " " + 2*v/(nc*(nc-1)) + " " + nc);*/
		/* Get the members of the cluster  */
		chromosomes = cluster.getClusterMembers();
		
		if (chromosomes == null) return null;
		
		if (cluster.isConnected(2)) {
			/* The number of chromosomes */
			nc = chromosomes.length;
			/* The number of uniq markers     */
			nm = cluster.getMarkerNumber();
			/* The cumulated number of marker intervals 
			 * which is also the number of equations to
			 * solve */
			ni = getMarkerIntervalNumber(chromosomes);
			/*
			 * Allocate memory for the numrical workspace
			 */
			distances = new double[ni];
			weights   = new double[ni];
			X         = new double[ni][nm];
			
			/*
			 * Loop on chromosomes
			 */
			for(xi=i=0;i<nc;i++) {
				/*
				 * Ignore chromosme if skeleton is true. Not that this may
				 * be just true one time... see further
				 */
				if (!chromosomes[i].skeleton) {
					
					/* Get the cross properties for this chromosome */
					type = chromosomes[i].ct;
					size = chromosomes[i].cs;
					gen  = chromosomes[i].cg;
					
					/*
					 * Loop on the marker intervals of the
					 * chromosome.
					 */
					for(j=0;j<chromosomes[i].nm-1;j++) {
						
						/* Get the recombination rate for this interval */
						r = chromosomes[i].mr[j];
						/* Convert it to genetic distance using additive Haldane
						 * distance
						 */
						d = MappingFunction.distance(
								chromosomes[i].mr[j],
								IMetaQtlConstants.MAPPING_FUNCTION_HALDANE,
								IMetaQtlConstants.MAPPING_UNIT_M
								);
						/* Set the indices of the flanking markers 
						 * 
						 * xj1 = left marker
						 * xj2 = right marker
						 * 
						 * */
						xj1 = j;
						xj2 = j+1;
						
						/* 
						 * First is the recombination rate is null
						 * for this interval ?
						 * 
						 * If it is the case, we have to look for the first
						 * marker at the right or at the left for which the
						 * distance to the current one is not too small...
						 * */
						if (d<=IMetaQtlConstants.TINY_MAPPING_DISTANCE) {
							/* Look at the left if we can find a not null
							 * distance. Otherwise look at the right
							 */
							for(l=j-1;l>=0;l--) {
								d += MappingFunction.distance(
											chromosomes[i].mr[l],
											IMetaQtlConstants.MAPPING_FUNCTION_HALDANE,
											IMetaQtlConstants.MAPPING_UNIT_M
											);
								if (d > IMetaQtlConstants.TINY_MAPPING_DISTANCE)
									break;
							} /* end left look */
							if (l>=0) {
								/* This means that we have find a marker at the left */
								xj1 = l;
							} else {
								for(l=j+1;l<chromosomes[i].nm-1;l++) {
									d += MappingFunction.distance(
												chromosomes[i].mr[l],
												IMetaQtlConstants.MAPPING_FUNCTION_HALDANE,
												IMetaQtlConstants.MAPPING_UNIT_M
												);
									if (d > IMetaQtlConstants.TINY_MAPPING_DISTANCE)
										break;
								} /* end right look */
								if (l<chromosomes[i].nm-1) {
									/* This means that we have find a marker at the right */
									xj2 = l+1;									
								} else {
									/* this must never occured, otherwise it means that we
									 * have a map with null distances, i.e of length zero.
									 */									
								}
							}
						} /* end minimal distance condition */
						/* 
						 * Convert the original marker indices to the
						 * absolute marker indices
						 * */
						x1 = NumericalUtilities.SIGN(1.0,(double)xj1-xj2);
						x2 = NumericalUtilities.SIGN(1.0,(double)xj2-xj1);
						xj1 = cluster.getMarkerIndexByName(chromosomes[i].mrkNames[xj1]);
						xj2 = cluster.getMarkerIndexByName(chromosomes[i].mrkNames[xj2]);
						/*
						 * Then, to understand what follows, recall that we try
						 * to solve the given equation :
						 * 
						 *  d = |xj2 - xj1| + e  avec e ~ N(0,(1/weights)^2)
						 * 
						 * This can be also written as follows,
						 * 
						 *  d = SIGN(xj2-xj1)*xj2 - SIGN(xj2-xj1)xj1 + e
						 *  
						 */
						X[xi][xj1]		= x1;
						X[xi][xj2] 		= x2;
						/* The observed distance */
						distances[xi]   = d;
						/* And the weight which is computed as the inverse
						 * variance of the distance using the Fisher Information
						 * assuming that the marker are fully informative. Of course
						 * this variance depends on the offspring population under
						 * consideration.
						 */
						r               = MappingFunction.recombination(d, IMetaQtlConstants.MAPPING_FUNCTION_HALDANE, IMetaQtlConstants.MAPPING_UNIT_M);
						weights[xi]		= OffspringManager.fisherVariance(r, size, type, gen);
						/* Then convert the variance of the recombination rate into that
						 * of distance
						 */ 
						weights[xi]     = MappingFunction.varianceDistance(
													r,
													weights[xi],
													IMetaQtlConstants.MAPPING_FUNCTION_HALDANE,
													IMetaQtlConstants.MAPPING_UNIT_M);
						/* try to avoid overflow */
						if (weights[xi] > IMetaQtlConstants.TINY_MAPPING_VARIANCE) {
							weights[xi] = 1.0/weights[xi];
						} else {
							weights[xi] = 0.0;
						}
						/* Increment the system row counter */
						xi++;
						
					} /* end loop on marker intervals */
					
				} /* end test if skeleton */
				
			} /* end loop on chromosomes */
			/* 
			 * Look if there is a skeleton chromosome in the cluster
			 * 
			 * Here we assume that there is only one skeleton in the cluster,
			 * that's why we stop to the first one.
			 * 
			 * See skeleton meaning for more details...
			 * 
			 */
			for(i=0;i<nc;i++) {
				
				if (chromosomes[i].skeleton) {
					
					d = 0.0; /* Set the distance to zero */
					/*
					 * Loop on skeleton markers.
					 */
					for(j=0;j<=chromosomes[i].nm;j++) {
						/*
						 * Get the absolute index for this marker
						 */
						xj1 = cluster.getMarkerIndexByName(chromosomes[i].mrkNames[j]);
						/*
						 * Move the coresponding column of the design matrix to
						 * the end and update the vector of distances
						 */
						for(l=xj1;l<nm-1;l++) {
							for(k=0;k<ni;k++) {
								distances[k]+=X[k][l]*d;
								X[k][l]=X[k][l+1];
							}
						}
						/*
						 * Cumulate distance
						 */
						d += MappingFunction.distance(
								chromosomes[i].mr[j-1],
								IMetaQtlConstants.MAPPING_FUNCTION_HALDANE,
								IMetaQtlConstants.MAPPING_UNIT_M
								);
					}
					
					nm -= chromosomes[i].nm;
					
					break;
				}
			}
			if (i==nc) {
				/* This means that any skeleton has been found. So we have to
				 * set one marker to zero to avoid overparametrization. This
				 * is done by decreasing the number of the columns of the 
				 * design matrix by one.
				 */
				nm--;
			}
			/* Check if there are weights to zero and then set them to the
			 * mean
			 */
			for(w=0.0,i=0;i<ni;i++) {
				weights[i] = Math.sqrt(weights[i]);
				w += weights[i];
			}
			w /= (double)ni;
			for(i=0;i<ni;i++) {
				//System.out.print("\n" + distances[i]+ " " + weights[i]);
				if (weights[i]==0.0) weights[i] = w;
				/*for(j=0;j<nm;j++) {
					System.out.print(" "+X[i][j]);
				}*/
			}
			/*
			 * Chi-square/WLS minization 
			 */
			result = SVDAlgorithm.SVDMFit(X, distances, weights, ni, nm);
			
			return new MetaChrom(result, cluster);
			
		} else {
			
			getLogger().println("[ WARNING ] Unable to build a consensus linkage groups for chromosome " + cluster.name);
			getLogger().println("            => Linkage groups are deconnected");
			
			return null;
			
		}
	}
	/**
	 * @param chromosomes
	 * @return
	 */
	private int getMarkerIntervalNumber(Chromosome[] chromosomes) {
		int ni=0;
		for(int i=0;i<chromosomes.length;i++)
			ni += (chromosomes[i].nm-1);
		return ni;
	}
	/**
	 * @param mrkDico
	 */
	public void setMarkerDico(MetaDico mrkDico) {
		this.mapCluster.setMarkerDico(mrkDico);
	}
	/**
	 * @param mrkThresh
	 */
	public void setMrkThreshold(int mrkThresh) {
		this.mrkThresh = mrkThresh;
	}
}
