/*  
 *  src/org/metaqtl/main/QTLClustInfo.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.main;

import java.io.FileWriter;
import java.io.IOException;
import java.io.PrintWriter;

import org.metaqtl.EMCriteria;
import org.metaqtl.EMResult;
import org.metaqtl.IMetaQtlConstants;
import org.metaqtl.MetaQtlAnalysis;
import org.metaqtl.MetaQtlResult;
import org.metaqtl.bio.util.NumberFormat;
import org.metaqtl.numrec.NumericalUtilities;

/**
 * 
 */
public class QTLClustInfo extends MetaMain {
	
	private static final String VERSION = IMetaQtlConstants.VERSION; 
	
	private static final String syntax  =  "Syntaxe: MetaQtlStat "+
										   "[{-r, --result}]] "+
										   "[{-o, --outfile}]] "+
										   "[{-c, --chr}]] "+
										   "[{-t, --trait}]] "+
										   "[{-b, --best}]] "+
										   "[{--kmin}]] "+
										   "[{--kmax}]] "+										   
										   MetaMain.generalUsage();
	
	public void printUsage() {
        System.err.println(syntax);
    }
	
	public void printHelp() {
		MetaMain.printLicense("MetaQtlStat", VERSION);
		System.out.println();
		System.out.println(syntax);
        System.out.println();
        MetaMain.generalHelp();
        System.out.println();
        System.out.println("-r, --clust   : the clustering result file");
        System.out.println("-o, --output  : the output file");
        System.out.println("-c, --chr     : the chromosome name");
        System.out.println("-t, --trait   : the trait name");
        System.out.println("-b, --best    : the best model (number of QTL)");
        System.out.println("--kmin        : the minimal number of QTL");
        System.out.println("--kmax        : the maximal number of QTL");
        
    }
	
	public static void main(String[] args) {
		
		QTLClustInfo program = new QTLClustInfo();
		
		program.initCmdLineParser();
		
		CmdLineParser parser = program.parser;
		
        CmdLineParser.Option aresFile   = parser.addStringOption('r', "clust");
        CmdLineParser.Option aoutFile   = parser.addStringOption('o', "output");
        CmdLineParser.Option achrName   = parser.addStringOption('c', "chr");
        CmdLineParser.Option atraitName = parser.addStringOption('t', "trait");
        CmdLineParser.Option abestK     = parser.addIntegerOption('b', "best");
        CmdLineParser.Option aminK      = parser.addIntegerOption("kmin");
        CmdLineParser.Option amaxK      = parser.addIntegerOption("kmax");
        
        program.parseCmdLine(args);
        
        // Mandatory
        String resFile  	=  (String)parser.getOptionValue(aresFile);
        String outFile  	=  (String)parser.getOptionValue(aoutFile);
        String chrName  	=  (String)parser.getOptionValue(achrName);
        String traitName  	=  (String)parser.getOptionValue(atraitName);
        Integer bestK		= (Integer)parser.getOptionValue(abestK);
        Integer minK		= (Integer)parser.getOptionValue(aminK);
        Integer maxK		= (Integer)parser.getOptionValue(amaxK);
        
        if (resFile == null) {
        	System.err.println("[ ERROR ] No result file defined");
        	System.exit(2);
        }
        if (outFile == null) {
        	System.err.println("[ ERROR ] No output file defined");
        	System.exit(2);
        }
        if (bestK == null) {
        	System.err.println("[ ERROR ] No best model defined");
        	System.exit(2);
        }
        if (chrName == null) {
        	System.err.println("[ ERROR ] No chromosome defined");
        	System.exit(2);
        }
        
        MetaQtlAnalysis analysis	= null;
		
        try {
        	analysis   = getResult(resFile);
       } catch (IOException e) {
			System.err.println(e.getMessage());
			System.exit(3);			
		}
		if (analysis == null) {
			System.err.println("[ ERROR ] Unable to load result from " + resFile);
			System.exit(4);
		}
		
		int cidx             = analysis.getChromIdx(chrName);
		
		MetaQtlResult result = analysis.getResult(chrName, traitName);
		
		if (result != null) {
			
			EMResult best = result.getClustering(bestK.intValue());
			// get the unconditional standard deviation estimates 
			double[] usd  = result.getUSD(best, minK.intValue(), maxK.intValue(), EMCriteria.AIC_NAME);
			// get the predicted values according to the best model
			double[] xp   = best.getXPred();
			// get the mahalanobis distance between components.
			double[] d    = best.getMahalanobis(result.sd);
			
			double ci_b   = NumericalUtilities.sd2ci(best.getSD(0));
			double ci_all = NumericalUtilities.sd2ci(usd[0]);
			
			double[] tmp = new double[xp.length+1];
			int[] idx    = new int[xp.length+1];
			for(int i=0;i<xp.length;i++) tmp[i+1]=xp[i];
			NumericalUtilities.indexx(xp.length, tmp, idx);
			
			
			try {
			
				FileWriter writer = new FileWriter(outFile);
				PrintWriter out   = new PrintWriter(writer);
				
				out.println("##");
				out.println("# Meta-QTL Table");
				out.println("##");
				out.println("--");
				out.println("QTL\tPosition\tWeight\tDistance\tCI(95%)\tUCI(95%)");
				out.println("--");
				for(int i=0;i<best.getK();i++) {
					ci_b   = NumericalUtilities.sd2ci(best.getSD(i));
					ci_all = NumericalUtilities.sd2ci(usd[i]);
					out.print((i+1)
							     +"\t"
								 +NumberFormat.formatDouble(best.getMu(i))
								 +"\t"
								 +NumberFormat.formatDouble(best.getPi(i))
								 +"\t");
					if (i<best.getK()-1)
						out.print(NumberFormat.formatDouble(d[i]));
					else
						out.print("-");
					out.println("\t"
								 +NumberFormat.formatDouble(ci_b)
								 +"\t"
								 +NumberFormat.formatDouble(ci_all));
				}
				out.println("--");
				out.println("##");
				out.println("# Observed QTL Table");
				out.println("##");
				out.println("--");
				out.print("QTL   Position   CI(95%)   Predicted   Memberships 1");
				if (bestK.intValue()>1) out.println(" --> "+bestK.intValue());
				else  out.println();
				out.println("--");
				for(int i=1;i<=xp.length;i++) {
					String qtlName = analysis.qtlByChrom[cidx][result.qtlIdx[idx[i]-1]].getName();
					out.print(qtlName
							  +"\t"
							  +NumberFormat.formatDouble(result.x[idx[i]-1])
							  +"\t"
							  +NumberFormat.formatDouble(NumericalUtilities.sd2ci(result.sd[idx[i]-1]))
							  +"\t"
							  +NumberFormat.formatDouble(xp[idx[i]-1]));
					double[] z=best.getZ(idx[i]-1);
					for(int j=0;j<z.length;j++)
						out.print("\t"
								  +NumberFormat.formatDouble(z[j]));
					out.println();
				}
				out.println("--");
				
				writer.close();
				
			} catch(IOException e) {
				
			}
		}
	 	
		System.exit(0);
	 	
	 }
	
	

}
