/**
 * Title:        ProAlign<p>
 * Description:  <p>
 * Copyright:    Copyright (c) Ari Loytynoja<p>
 * License:      GNU GENERAL PUBLIC LICENSE<p>
 * @see          http://www.gnu.org/copyleft/gpl.html
 * Company:      ULB<p>
 * @author       Ari Loytynoja
 * @version      1.0
 */
package proalign;

import java.util.HashMap;

public class ParameterEstimates {

    String[][] pairs;
    HashMap seqs;
    PwAlignment pwa;

    ParameterEstimates(ResultWindow rw,SetParameters sp) { 

	this.seqs = rw.seqs;

	PwSubstitutionMatrix psm = new PwSubstitutionMatrix();
	String pwAlphabet;
	int[][] pwSubst;
	int gOpen, gExt;

	if(ProAlign.isDna){
	    pwAlphabet = psm.dnaAlphabet;
	    pwSubst = psm.swdna;
	    gOpen = -1*rw.pa.pwDnaOpen;
	    gExt = -1*rw.pa.pwDnaExt;

	} else {

	    pwAlphabet = psm.protAlphabet;
	    if(rw.pa.pwProtMatrix.equals("pam60")) {
		pwSubst = psm.pam60;
	    } else if(rw.pa.pwProtMatrix.equals("pam160")) {
		pwSubst = psm.pam160;
	    } else if(rw.pa.pwProtMatrix.equals("pam250")) {
		pwSubst = psm.pam250;
	    } else {
		pwSubst = psm.pam120;
	    }
	    gOpen = -1*rw.pa.pwProtOpen;
	    gExt = -1*rw.pa.pwProtExt;
	}
	
	pwa = new PwAlignment(pwSubst,gOpen,gExt,pwAlphabet,ProAlign.isDna);

	pairs = rw.root.getTerminalPairNames();

	this.estimate();

	sp.textDelta.setText((""+ProAlign.modelDelta+"      ").substring(0,5));
	sp.textEpsil.setText((""+ProAlign.modelEpsilon+"      ").substring(0,5));
	
    }

    ParameterEstimates(RunCommandLine rcl) {

	ProAlign.log("ParameterEstimates");

	AlignmentNode root = rcl.root;

	seqs = rcl.pa.seqs;
	pwa = rcl.pwa;

	pairs = root.getTerminalPairNames();

	this.estimate();
    }

    void estimate() {

	double gapFreq = 0d;
	double sumDelta = 0d;
	double sumEpsilon = 0d;

	for(int i=0; i<pairs.length; i++) {

	    String s0 = (String)seqs.get(pairs[i][0]);
	    String s1 = (String)seqs.get(pairs[i][1]);

	    String[] revSeq =  pwa.revAligned(s0,s1);
	    
	    int all = 0;
	    int ms = 0; // match sum
	    int mn = 1; // match numebr
	    int ml = 0; // match length
	    int gs = 0;
	    int gn = 1;
	    int gl = 0;
	    double mm = 5d; // match mean
	    double mg = 5d;
	    boolean isM = false;
	    
	    for(int k=0; k<revSeq[0].length(); k++) {
		if(k==0) {
		    if(revSeq[0].charAt(k)=='-' || revSeq[1].charAt(k)=='-') {
			gn++;
		    } else {
			mn++;
		    }
		}
		
		if(revSeq[0].charAt(k)=='-' || revSeq[1].charAt(k)=='-') {
		    gs++;
		    gl++;
		    if(isM) {
			if(k>0) {
			    if(mn>1) {
				mm = ((double)mm*(mn-1)+ml)/mn;
			    } else {
				mm = (mm+ml)/2;
			    }
			    ml=0;
			    gn++;
			}
			isM = false;				
		    }
		} else {
		    ms++;
		    ml++;
		    if(!isM) {		    
			if(k>0) {
			    if(gn>1) {
				mg = ((double)mg*(gn-1)+gl)/gn;
			    } else {
				mg = (mg+gl)/2;
			    }
			    gl=0;
			    mn++;
			}
			isM = true;				
		    }
		}
		
		if(k+1==revSeq[0].length()) {
		    if(isM) {
			if(mn>1) {
			    mm = ((double)mm*(mn-1)+ml)/mn;
			} else {
			    mm = ml;
			}
		    } else {
			if(gn>1) {
			    mg = ((double)mg*(gn-1)+gl)/gn;
			} else {
			    mg = gl;
			}
		    }
		}
		all++;
	    }
	    gapFreq += (double)gs/(2*all);
	    sumDelta += mm;
	    sumEpsilon += mg;
	}

	gapFreq /= pairs.length;
	sumDelta /= pairs.length;
	sumEpsilon /= pairs.length;

	if(ProAlign.estimateDelta) {
	    ProAlign.modelDelta = 0.5d/(sumDelta+1);
	    ProAlign.log.println(" HMM delta estimate: "+ProAlign.modelDelta);
//	    System.out.println("modelDelta: "+ProAlign.modelDelta);	    
	}
	if(ProAlign.estimateEpsilon) {
	    ProAlign.modelEpsilon = (1d-1d/(sumEpsilon+1));
	    ProAlign.log.println(" HMM epsilon estimate: "+ProAlign.modelEpsilon);
//	    System.out.println("modelEpsilon: "+ProAlign.modelEpsilon);
	}
	if(ProAlign.estimateGapFreq) {
	    ProAlign.gapFreq = gapFreq;
	    ProAlign.log.println(" gap frequency estimate: "+ProAlign.gapFreq);
//	    System.out.println("gapFreq: "+ProAlign.gapFreq);
	}
	if(ProAlign.estimateGapProb) {
	    if(ProAlign.isDna) {
		ProAlign.gapProb = gapFreq/4d;
	    } else {
		ProAlign.gapProb = gapFreq/20d;
	    }
	    ProAlign.log.println(" gap probability estimate: "+ProAlign.gapProb);
//	    System.out.println("gapProb: "+ProAlign.gapProb);
	}

    }	

}