/*  
 *  src/org/metaqtl/algo/CMSAlgorithm.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.CMarkerSequence;
import org.metaqtl.Chromosome;
import org.metaqtl.numrec.NumericalUtilities;
/**
 * Use this class to find common marker sequences between
 * two chromosomes.
 * 
 * @see org.metaqtl.Chromosome
 * @see org.metaqtl.CMarkerSequence
 */
public final class CMSAlgorithm {

	/**
	 * This methods find the common marker sequences between the 
	 * two given chromosomes.
	 * 
	 * @param chrom1 the first <code>Chromosome</code>.
	 * @param chrom2 the second <code>Chromosome</code>.
	 * @return a <code>CMarkerSequence</code> 
	 */
	public static CMarkerSequence run(Chromosome chrom1, Chromosome chrom2)
	{
		boolean frame,first,ok,swap;
		int i,j,nmc,id1,id2,lastIdx;
		double d0,d1;
		int[] idx1;
		int[][] mcidx;
		double[] tmp1;
		CMarkerSequence cms = null;
		/* 
		 * First get only the common markers
		 * between the two chromosomes and store
		 * the indices into mcidx[] array.
		 */
		mcidx = new int[2][Math.min(chrom1.nm, chrom2.nm)];
		tmp1  = new double[Math.min(chrom1.nm, chrom2.nm)+1];
		
		for(nmc=i=0;i<chrom1.nm;i++) {
			for(j=0;j<chrom2.nm;j++) {
				if (chrom1.mrkNames[i].equals(chrom2.mrkNames[j])) {
					id1=i;
					id2=j;
					if (nmc > 0) {
						d0 = chrom1.getDistanceBetween(mcidx[0][nmc-1], i);
						if (j > mcidx[1][nmc-1]) d1 = chrom2.getDistanceBetween(mcidx[1][nmc-1], j);
						else d1 = chrom2.getDistanceBetween(j, mcidx[1][nmc-1]);
						ok=(d0 != 0.0 && d1 != 0.0);
						swap=(d0==0.0 && j < mcidx[1][nmc]);
					} else {
						ok=true;
						swap=false;
					}
					if (ok) {
						mcidx[0][nmc]=id1;
						mcidx[1][nmc]=id2;
						tmp1[nmc+1]=(double)id2;
						if (swap) {
							NumericalUtilities.ISWAP(mcidx[0], nmc-1, nmc);
							NumericalUtilities.ISWAP(mcidx[1], nmc-1, nmc);
							NumericalUtilities.DSWAP(tmp1, nmc, nmc+1);
						}
						nmc++;
					}
				}
			}
		}
		
		if (nmc>=1)
			cms = new CMarkerSequence(mcidx, nmc);
		
		if (nmc > 2) {
			/* First re-index the common marker array */
			idx1  = new int[nmc+1];
			NumericalUtilities.indexx(nmc,tmp1,idx1);
			
			lastIdx=-1;
			
			i=0;
			
			id1 = mcidx[1][0];//idx1[1]-1;
			for(j=1;j<=nmc;j++)
				if (idx1[j]==i+1) break;
			id1 = j-1;	
			
			first=frame=true; /* the frame of the marker order */
			
			for(i=1;i<nmc;i++) {
				
				id2 = mcidx[1][i];//idx1[i+1]-1;
				for(j=1;j<=nmc;j++)
					if (idx1[j]==i+1) break;
				id2 = j-1;
				
				if (first) { 
					if (id2 > id1) frame=true;
					else frame=false;
				}
				
				ok=isInCMS(id1,id2,frame);
				
				if (ok) {
					
					if (lastIdx == id1) {
						/* Add this marker to the previous common sequence */
						j=++cms.css[cms.ncs-1]; /* increment the number of marker of the current cs */
					} else {
						/* Create a new common sequence */
						cms.ncs++;
						cms.frames[cms.ncs-1]=frame;
						j=cms.css[cms.ncs-1]=2;
						cms.idx1[cms.ncs-1][0]=mcidx[0][i-1];//mcidx[0][id1];
						cms.idx2[cms.ncs-1][0]=mcidx[1][i-1];//mcidx[1][id1];
						//cms.incms[id1]=true;
						cms.incms[i-1]=true;
					}
					cms.idx1[cms.ncs-1][j-1]=mcidx[0][i];//mcidx[0][id2];
					cms.idx2[cms.ncs-1][j-1]=mcidx[1][i];//mcidx[1][id2];
					//cms.incms[id2]=true;
					cms.incms[i]=true;
					lastIdx=id2;
					first=false;
					
				} else {
					
					first=true;
					
				}
				
				id1=id2;
				
			}
		
		} else if (nmc == 2 && chrom1.mrkNames[mcidx[0][0]].equals(chrom2.mrkNames[mcidx[1][0]])) {
			cms.ncs=1;
			j=cms.css[0]=2;
			cms.frames[0]=true;
			cms.idx1[0][0]=mcidx[0][0];
			cms.idx2[0][0]=mcidx[1][0];
			cms.idx1[0][1]=mcidx[0][1];
			cms.idx2[0][1]=mcidx[1][1];
			cms.incms[0]=cms.incms[1]=true;
		} 
	    
	    return cms;
	}
	/**
	 * 
	 * @param id1
	 * @param id2
	 * @param frame
	 * @param d1
	 * @param d2
	 * @return
	 */
	private static boolean isInCMS(int id1, int id2, boolean frame) {
		return ((id2 == id1+1 && frame) || (id2 == id1-1 && !frame));		
	}	
}
