package fmph.features.util;

import calhoun.analysis.crf.Conrad;
import calhoun.analysis.crf.features.supporting.phylogenetic.EvolutionaryModel;
import calhoun.analysis.crf.features.supporting.phylogenetic.HKY85Model;
import calhoun.analysis.crf.features.supporting.phylogenetic.Kimura80Model;
import calhoun.analysis.crf.features.supporting.phylogenetic.PhylogeneticTreeFelsensteinOrder;
import calhoun.analysis.crf.io.FilenameMapper;
import calhoun.analysis.crf.io.InputSequence;
import calhoun.analysis.crf.io.MultipleAlignmentInputSequence;
import calhoun.analysis.crf.io.TrainingSequence;
import calhoun.analysis.crf.io.MultipleAlignmentInputSequence.MultipleAlignmentColumn;

import calhoun.analysis.crf.io.NameInputSequence;

import calhoun.util.Assert;

import fmph.features.supporting.phylogenetic.CodonEvolutionaryModel;

import fmph.features.supporting.phylogenetic.GoldmanYang94;

import java.io.FileNotFoundException;
import java.io.FileOutputStream;

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

import java.util.ArrayList;
import java.util.List;

public class PhyloViewer {
    EvolutionaryModel emodelIntergenic,emodelIntronic;
    ArrayList<EvolutionaryModel> emodelExonic;
    String nameComponent, multipleAlnComponent, sourceDir;
    
    
    
    CodonEvolutionaryModel emodelCodonExonic,emodelCodonNoncoding;
    
   
    double[] emodelExonicParams = { 
            //codon freq
    0.06860156,0.01163551,0.06247777,0.00667686,
    0.03160051,0.00535196,0.02873777,0.00307114,
    0.03568897,0.00606133,0.00010000,0.00010000,
    0.02109063,0.00357198,0.01918000,0.00204972,
    0.02191554,0.00371168,0.01993017,0.00212989,
    0.01008044,0.00170726,0.00916724,0.00097968,
    0.01141654,0.00193354,0.01038230,0.00110953,
    0.00672783,0.00113945,0.00611835,0.00065385,
    0.06803812,0.01152315,0.06187443,0.00661238,
    0.03129535,0.00530028,0.02846025,0.00304148,
    0.03544335,0.00600280,0.03223248,0.00344461,
    0.02088696,0.00353748,0.01899478,0.00202993,
    0.04945072,0.00837513,0.04497090,0.00480594,
    0.02274574,0.00385229,0.02068517,0.00221058,
    0.02576055,0.00436289,0.02342686,0.00250358,
    0.01518083,0.00257107,0.01380558,//0.00147537,
    0.65454,//kappa
    0.08654,//omega
    4 //icode
            };
    
    double[] emodelCodonNoncodingParams = { 
            //codon freq
            0.04849882,0.01136236,0.05020109,0.01369843,
            0.01142529,0.00267673,0.01182630,0.00322706,
            0.04754752,0.01113183,0.00100000,0.00100000,
            0.01617676,0.00378991,0.01674456,0.00456911,
            0.01487535,0.00348501,0.01539747,0.00420153,
            0.00350432,0.00082099,0.00362731,0.00098979,
            0.01458664,0.00341737,0.01509862,0.00411998,
            0.00496167,0.00116243,0.00513582,0.00140142,
            0.07856974,0.01840741,0.08132748,0.02219193,
            0.01850935,0.00433639,0.01915902,0.00522794,
            0.07704480,0.01805014,0.07974901,0.02176121,
            0.02620691,0.00613978,0.02712675,0.00740211,
            0.02138010,0.00500895,0.02213053,0.00603879,
            0.00503670,0.00118000,0.00521348,0.00142261,
            0.02096514,0.00491174,0.02170100,0.00592158,
            0.00713133,0.00167074,0.00738163,//0.00201423,
            1,//kappa
            1,//omega
            4 //icode
            };
    
    public static void main(String[] args) throws Exception {
        new PhyloViewer(args);
    }
    
    public PhyloViewer(String[] args) throws IOException{
        nameComponent = args[2]; 
        multipleAlnComponent = args[3]; 
        sourceDir = args[1]; 
        
        
        Conrad crf = new Conrad(args[0]);
        List<? extends TrainingSequence<?>> t = crf.getInputHandler().readTrainingData(args[1], false);
        PhylogeneticTreeFelsensteinOrder felsOrder = ((MultipleAlignmentColumn)t.get(0).getComponent(args[3]).getX(0)).getMultipleAlignment().getFelsensteinOrder();
        PhylogeneticTreeFelsensteinOrder felsOrder2 = ((MultipleAlignmentColumn)t.get(0).getComponent("aln2").getX(0)).getMultipleAlignment().getFelsensteinOrder();
        
        //emodelIntergenic = new EvolutionaryModel(felsOrder, new double[]{0.44448760556490713,0.05551239443509287,0.05551239443509287,0.44448760556490713}, new Kimura80Model(new double[]{0.23740500177751464,0.6536810457617929}));
        
        emodelIntergenic = new EvolutionaryModel(felsOrder, new double[]{0.44448760556490713,0.06946919299912885,0.06946919299912885,0.44448760556490713}, new HKY85Model(new double[]{9.027941760821149,4.033718455425241,0.44448760556490713,0.06946919299912885,0.05551239443509287}));
        //emodelIntronic = new EvolutionaryModel(felsOrder, new double[]{0.43914528597477975,0.08378919793512597,0.1197832395921617,0.35728227649793254}, new Kimura80Model(new double[]{0.6492450873908504,0.9403504300614604}));
        emodelIntronic = new EvolutionaryModel(felsOrder, new double[]{0.43914528597477975,0.08378919793512597,0.1197832395921617,0.35728227649793254}, new HKY85Model(new double[]{7.1637616103302575,3.975049018330243,0.43914528597477975,0.08378919793512597,0.1197832395921617}));
        emodelExonic = new ArrayList<EvolutionaryModel>();
        //emodelExonic.add(new EvolutionaryModel(felsOrder, new double[]{0.337173765122421,0.09340437555501736,0.2211996447888916,0.3482222145336701}, new Kimura80Model(new double[]{0.29820919354158326,0.20188192420999146})));
        emodelExonic.add(new EvolutionaryModel(felsOrder, new double[]{0.337173765122421,0.09340437555501736,0.2211996447888916,0.3482222145336701}, new HKY85Model(new double[]{2.3728814207126065,1.425772860742017,0.337173765122421,0.09340437555501736,0.2211996447888916})));
        //emodelExonic.add(new EvolutionaryModel(felsOrder, new double[]{0.24048632499336725,0.19082719082719082,0.12942519984773507,0.43926128433170686}, new Kimura80Model(new double[]{0.22162973326046825,0.1656976351280791})));
        emodelExonic.add(new EvolutionaryModel(felsOrder, new double[]{0.24048632499336725,0.19082719082719082,0.12942519984773507,0.43926128433170686}, new HKY85Model(new double[]{1.7390188901730215,1.2445736320384442,0.24048632499336725,0.19082719082719082,0.12942519984773507})));
        //emodelExonic.add(new EvolutionaryModel(felsOrder, new double[]{0.5036354823073195,0.05671352399418323,0.04398356531172818,0.39566742838676916}, new Kimura80Model(new double[]{0.25051188519131423,0.40170737974641324})));
        emodelExonic.add(new EvolutionaryModel(felsOrder, new double[]{0.5036354823073195,0.05671352399418323,0.04398356531172818,0.39566742838676916}, new HKY85Model(new double[]{6.682136848355302,2.7266919623589505,0.5036354823073195,0.05671352399418323,0.04398356531172818})));
        
        double[] codonNoncodingPi = new double[64];
        double[] codonExonicPi = new double[64];
        System.arraycopy(emodelExonicParams,0,codonExonicPi,0,63);
        System.arraycopy(emodelCodonNoncodingParams,0,codonNoncodingPi,0,63);
        
        double sum = 0;
        for(int i=0;i<63;i++){
            sum += emodelExonicParams[i];
        }
        codonExonicPi[63] = 1 - sum;
        
        sum = 0;
        for(int i=0;i<63;i++){
            sum += emodelCodonNoncodingParams[i];
        }
        codonNoncodingPi[63] = 1 - sum;
        
        emodelCodonExonic = new CodonEvolutionaryModel(felsOrder2, codonExonicPi, new GoldmanYang94(emodelExonicParams));
        emodelCodonNoncoding = new CodonEvolutionaryModel(felsOrder2, codonNoncodingPi, new GoldmanYang94(emodelCodonNoncodingParams));
        //emodelCodonExonic.summarize();
        //emodelIntronic.summarize();
        //System.out.println(emodelCodonExonic.logprob(null,null,null,false));
        //printValues(t,args);
        printCodonValues3(t);
    }
    
    private void printValues(List<? extends TrainingSequence<?>> t, String[] args) throws FileNotFoundException {
        //PrintWriter pwInter = new PrintWriter(new FileOutputStream(intergenicOutputFile));
        //PrintWriter pwIntron = new PrintWriter(new FileOutputStream(intronOutputFile));
        PrintWriter pwExon = new PrintWriter(new FileOutputStream(sourceDir+"/comparative.values"));
        
        for(TrainingSequence<?> seq:t){
            String name = ((NameInputSequence)seq.getInputSequence().getComponent(nameComponent)).getName();
            //pwInter.println(">"+name);
            //pwIntron.println(">"+name);
            pwExon.println(">"+name);
            InputSequence<? extends MultipleAlignmentColumn> aligSeq = (InputSequence<? extends MultipleAlignmentColumn>)seq.getComponent(multipleAlnComponent);
            for(int pos=0;pos<aligSeq.length();pos++){
                MultipleAlignmentColumn col = aligSeq.getX(pos);
                int state = seq.getY(pos);
                int ephase;
                int featureOffset = Integer.MIN_VALUE;
                switch (state) {
                case 0:
                    //pwInter.println(emodelIntergenic.logprob(col, true));
                    //pwIntron.println(emodelIntronic.logprob(col,true));
                    //pwExon.println(emodelExonic.get(1).logprob(col,true));
                    double intronicValue = emodelIntronic.logprob(col,true);
                    pwExon.println(
                    Math.max(emodelExonic.get(0).logprob(col,true)-intronicValue,
                        Math.max( emodelExonic.get(1).logprob(col,true)-intronicValue,
                                  emodelExonic.get(2).logprob(col,true)-intronicValue)));
                    break;
                case 1:
                case 2:
                case 3:
                    ephase = ((pos - state + 1) % 3 + 3) % 3; //((pos-(state-1))%3 +3)%3;
                    //val = emodelExonic.get(0).logprob(col,true);
                    //pwInter.println(emodelIntergenic.logprob(col, true));
                    //pwIntron.println(emodelIntronic.logprob(col,true));
                    //pwExon.println(emodelExonic.get(ephase).logprob(col,true));
                    pwExon.println(emodelExonic.get(ephase).logprob(col,true)-emodelIntronic.logprob(col,true));
                    break;
                case 4:
                case 5:
                case 6:
                    //pwInter.println(emodelIntergenic.logprob(col, true));
                    //pwIntron.println(emodelIntronic.logprob(col,true));
                    //pwExon.println(emodelExonic.get(1).logprob(col,true));
                    intronicValue = emodelIntronic.logprob(col,true);
                    pwExon.println(
                    Math.max(emodelExonic.get(0).logprob(col,true)-intronicValue,
                        Math.max( emodelExonic.get(1).logprob(col,true)-intronicValue,
                                  emodelExonic.get(2).logprob(col,true)-intronicValue)));
                    break;
                case 7:
                case 8:
                case 9:
                    ephase = ((-pos + state + 1) % 3 + 3) % 3; // ((-pos+2+(state-7))%3 +3)%3;
                    //pwInter.println(emodelIntergenic.logprob(col, true));
                    //pwIntron.println(emodelIntronic.logprobRC(col,true));
                    //pwExon.println(emodelExonic.get(ephase).logprobRC(col,true));
                    pwExon.println(emodelExonic.get(ephase).logprobRC(col,true)-emodelIntronic.logprobRC(col,true));
                    break;
                case 10:
                case 11:
                case 12:
                    //pwInter.println(emodelIntergenic.logprob(col, true));
                    //pwIntron.println(emodelIntronic.logprobRC(col,true));
                    //pwExon.println(emodelExonic.get(1).logprobRC(col,true));
                    intronicValue = emodelIntronic.logprobRC(col,true);
                    pwExon.println(
                    Math.max( 
                        emodelExonic.get(0).logprobRC(col,true)-intronicValue,
                        Math.max( emodelExonic.get(1).logprobRC(col,true)-intronicValue,
                                  emodelExonic.get(2).logprobRC(col,true)-intronicValue)));
                    break;
                default:
                    Assert.a(false);
                }
            }
        }
        
        //pwInter.close();
        //pwIntron.close();
        pwExon.close();
        
    }
    
    private void printCodonValues(List<? extends TrainingSequence<?>> t) throws FileNotFoundException {
        //PrintWriter pwInter = new PrintWriter(new FileOutputStream(intergenicOutputFile));
        //PrintWriter pwIntron = new PrintWriter(new FileOutputStream(intronOutputFile));
        PrintWriter pwExon = new PrintWriter(new FileOutputStream(sourceDir+"/codon.values"));
        //PrintWriter pwExon2 = new PrintWriter(new FileOutputStream(exonOutputFile+".codon2"));
        int ix = 0;
        for(TrainingSequence<?> seq:t){
            System.out.println(++ix);
            String name = ((NameInputSequence)seq.getInputSequence().getComponent(nameComponent)).getName();
            //pwInter.println(">"+name);
            //pwIntron.println(">"+name);
            pwExon.println(">"+name);
            //pwExon2.println(">"+name);
            InputSequence<? extends MultipleAlignmentColumn> aligSeq = (InputSequence<? extends MultipleAlignmentColumn>)seq.getComponent(multipleAlnComponent);
            for(int pos=0;pos<aligSeq.length();pos++){
                
                int state = seq.getY(pos);
                int ephase;
                int featureOffset = Integer.MIN_VALUE;
                MultipleAlignmentColumn colm2 = null;
                MultipleAlignmentColumn colm1 = null;
                MultipleAlignmentColumn col0 = null;
                MultipleAlignmentColumn col1 = null;
                MultipleAlignmentColumn col2 = null;
                switch (state) {
                case 1:
                case 2:
                case 3:
                    ephase = ((pos - state + 1) % 3 + 3) % 3; //((pos-(state-1))%3 +3)%3;
                    if(ephase != 0) break;
                    col0 = aligSeq.getX(pos-ephase);
                    col1 = aligSeq.getX(pos-ephase+1);
                    col2 = aligSeq.getX(pos-ephase+2);
                    //pwExon.println(emodelCodonExonic.logprob(col0,col1,col2,true));
                    //val = emodelExonic.get(0).logprob(col,true);
                    //pwInter.println(emodelIntergenic.logprob(col, true));
                    //pwInter.println(emodelCodonNoncoding.logprob(col0,col1,col2, true));
                    //pwIntron.println(emodelCodonNoncoding.logprob(col0,col1,col2,true));
                    double value = emodelCodonExonic.logprob(col0,col1,col2,true) - (emodelIntronic.logprob(col0,true) + 
                                                                                    emodelIntronic.logprob(col1,true) + 
                                                                                    emodelIntronic.logprob(col2,true));
                    pwExon.println(value);
                    pwExon.println(value);
                    pwExon.println(value);
                    //pwExon2.println((emodelCodonExonic.logprob(col0,col1,col2,false) - emodelCodonNoncoding.logprob(col0,col1,col2,false))/3);
                    /*pwExon.println(emodelCodonExonic.logprob(col0,col1,col2,false)*-2.499858e-03 + ((col0==null?0:emodelIntronic.logprob(col0,false)) + 
                                                                                                        (col1==null?0:emodelIntronic.logprob(col1,false)) + 
                                                                                                        (col2==null?0:emodelIntronic.logprob(col2,false)))*4.687246e-04);
                    
                     */break;
                case 0:
                case 4:
                case 5:
                case 6:
                    colm2 = pos>1?aligSeq.getX(pos-2):null;
                    colm1 = pos>0?aligSeq.getX(pos-1):null;
                    col0 = aligSeq.getX(pos);
                    col1 = pos<aligSeq.length()-1?aligSeq.getX(pos+1):null;
                    col2 = pos<aligSeq.length()-2?aligSeq.getX(pos+2):null;
                    //pwExon.println(emodelCodonExonic.logprob(col0,col1,col2,true));
                    //pwInter.println(emodelCodonNoncoding.logprob(col0,col1,col2, true));
                    //pwIntron.println(emodelCodonNoncoding.logprob(col0,col1,col2,true));
                    //pwExon.println(emodelCodonExonic.logprob(col0,col1,col2,true));
                    pwExon.println(
                    Math.max( 
                        emodelCodonExonic.logprob(colm2,colm1,col0,true)- ((colm2==null?0:emodelIntronic.logprob(colm2,true)) + 
                                                                           (colm1==null?0:emodelIntronic.logprob(colm1,true)) + 
                                                                           (col0==null?0:emodelIntronic.logprob(col0,true))),
                        Math.max( emodelCodonExonic.logprob(colm1,col0,col1,true)-((colm1==null?0:emodelIntronic.logprob(colm1,true)) + 
                                                                           (col0==null?0:emodelIntronic.logprob(col0,true)) + 
                                                                           (col1==null?0:emodelIntronic.logprob(col1,true))),
                                emodelCodonExonic.logprob(col0,col1,col2,true)-((col0==null?0:emodelIntronic.logprob(col0,true)) + 
                                                                           (col1==null?0:emodelIntronic.logprob(col1,true)) + 
                                                                           (col2==null?0:emodelIntronic.logprob(col2,true))))));
                    /*pwExon.println((emodelCodonExonic.logprob(col0,col1,col2,true)-((col0==null?0:emodelIntronic.logprob(col0,true)) + 
                                                               (col1==null?0:emodelIntronic.logprob(col1,true)) + 
                                                               (col2==null?0:emodelIntronic.logprob(col2,true))))/3);*/
                    //pwExon2.println((emodelCodonExonic.logprob(col0,col1,col2,false) - emodelCodonNoncoding.logprob(col0,col1,col2,false))/3);
                    /*pwExon.println(emodelCodonExonic.logprob(col0,col1,col2,false)*-2.499858e-03 + ((col0==null?0:emodelIntronic.logprob(col0,false)) + 
                                                                                                        (col1==null?0:emodelIntronic.logprob(col1,false)) + 
                                                                                                        (col2==null?0:emodelIntronic.logprob(col2,false)))*4.687246e-04);
                    
                    */break;
                case 7:
                case 8:
                case 9:
                    ephase = ((-pos + state + 1) % 3 + 3) % 3; //((pos-(state-1))%3 +3)%3;
                    if(ephase != 0) break;
                    col0 = aligSeq.getX(pos+ephase-2);
                    col1 = aligSeq.getX(pos+ephase-1);
                    col2 = aligSeq.getX(pos+ephase);
                    //pwExon.println(emodelCodonExonic.logprobRC(col0,col1,col2,true));
                    //val = emodelExonic.get(0).logprob(col,true);
                    //pwInter.println(emodelIntergenic.logprob(col, true));
                    //pwInter.println(emodelCodonNoncoding.logprobRC(col0,col1,col2, true));
                    //pwIntron.println(emodelCodonNoncoding.logprobRC(col0,col1,col2,true));
                    value = emodelCodonExonic.logprobRC(col0,col1,col2,true) - (emodelIntronic.logprobRC(col0,true) + 
                                                                                       emodelIntronic.logprobRC(col1,true) + 
                                                                                       emodelIntronic.logprobRC(col2,true));
                    pwExon.println(value);
                    pwExon.println(value);
                    pwExon.println(value);
                    //pwExon2.println((emodelCodonExonic.logprob(col0,col1,col2,false) - emodelCodonNoncoding.logprobRC(col0,col1,col2,false))/3);
                    /*pwExon.println(emodelCodonExonic.logprobRC(col0,col1,col2,false)*-2.499858e-03 + ((col0==null?0:emodelIntronic.logprobRC(col0,false)) + 
                                                                                                        (col1==null?0:emodelIntronic.logprobRC(col1,false)) + 
                                                                                                        (col2==null?0:emodelIntronic.logprobRC(col2,false)))*4.687246e-04);
                    */break;
                case 10:
                case 11:
                case 12:
                    colm2 = pos>1?aligSeq.getX(pos-2):null;
                    colm1 = pos>0?aligSeq.getX(pos-1):null;
                    col0 = aligSeq.getX(pos);
                    col1 = pos<aligSeq.length()-1?aligSeq.getX(pos+1):null;
                    col2 = pos<aligSeq.length()-2?aligSeq.getX(pos+2):null;
                    //pwInter.println(emodelCodonNoncoding.logprob(col0,col1,col2, true));
                    //pwIntron.println(emodelCodonNoncoding.logprob(col0,col1,col2,true));
                    //pwExon.println(emodelCodonExonic.logprobRC(col0,col1,col2,true));
                    pwExon.println(
                    Math.max( 
                        emodelCodonExonic.logprobRC(colm2,colm1,col0,true)- ((colm2==null?0:emodelIntronic.logprobRC(colm2,true)) + 
                                                                           (colm1==null?0:emodelIntronic.logprobRC(colm1,true)) + 
                                                                           (col0==null?0:emodelIntronic.logprobRC(col0,true))),
                        Math.max( emodelCodonExonic.logprobRC(colm1,col0,col1,true)-((colm1==null?0:emodelIntronic.logprobRC(colm1,true)) + 
                                                                           (col0==null?0:emodelIntronic.logprobRC(col0,true)) + 
                                                                           (col1==null?0:emodelIntronic.logprobRC(col1,true))),
                                emodelCodonExonic.logprobRC(col0,col1,col2,true)-(col0==null?0:emodelIntronic.logprobRC(col0,true) + 
                                                                           (col1==null?0:emodelIntronic.logprobRC(col1,true)) + 
                                                                           (col2==null?0:emodelIntronic.logprobRC(col2,true))))));
                    /*pwExon.println((emodelCodonExonic.logprobRC(col0,col1,col2,true)-(col0==null?0:emodelIntronic.logprobRC(col0,true) + 
                                                               (col1==null?0:emodelIntronic.logprobRC(col1,true)) + 
                                                               (col2==null?0:emodelIntronic.logprobRC(col2,true))))/3);*/
                    //pwExon2.println((emodelCodonExonic.logprob(col0,col1,col2,false) - emodelCodonNoncoding.logprobRC(col0,col1,col2,false))/3);
                    /*pwExon.println(emodelCodonExonic.logprobRC(col0,col1,col2,false)*-2.499858e-03 + ((col0==null?0:emodelIntronic.logprobRC(col0,false)) + 
                                                                                                        (col1==null?0:emodelIntronic.logprobRC(col1,false)) + 
                                                                                                        (col2==null?0:emodelIntronic.logprobRC(col2,false)))*4.687246e-04);
                    */break;
                default:
                    Assert.a(false);
                }
            }
        }
        
        //pwInter.close();
        //pwIntron.close();
        pwExon.close();
        //pwExon2.close();
        
    }
    
    private void printCodonValues2(List<? extends TrainingSequence<?>> t) throws FileNotFoundException {
        //PrintWriter pwInter = new PrintWriter(new FileOutputStream(intergenicOutputFile));
        //PrintWriter pwIntron = new PrintWriter(new FileOutputStream(intronOutputFile));
        PrintWriter pwExon = new PrintWriter(new FileOutputStream(sourceDir+"/codon2.values"));
        //PrintWriter pwExon2 = new PrintWriter(new FileOutputStream(exonOutputFile+".codon2"));
        int ix = 0;
        for(TrainingSequence<?> seq:t){
            System.out.println(++ix);
            String name = ((NameInputSequence)seq.getInputSequence().getComponent(nameComponent)).getName();
            //pwInter.println(">"+name);
            //pwIntron.println(">"+name);
            pwExon.println(">"+name);
            //pwExon2.println(">"+name);
            InputSequence<? extends MultipleAlignmentColumn> aligSeq = (InputSequence<? extends MultipleAlignmentColumn>)seq.getComponent(multipleAlnComponent);
            for(int pos=0;pos<aligSeq.length();pos++){
                
                int state = seq.getY(pos);
                int ephase;
                int featureOffset = Integer.MIN_VALUE;
                MultipleAlignmentColumn colm2 = null;
                MultipleAlignmentColumn colm1 = null;
                MultipleAlignmentColumn col0 = null;
                MultipleAlignmentColumn col1 = null;
                MultipleAlignmentColumn col2 = null;
                switch (state) {
                case 1:
                case 2:
                case 3:
                    ephase = ((pos - state + 1) % 3 + 3) % 3; //((pos-(state-1))%3 +3)%3;
                    if(ephase != 0) break;
                    col0 = aligSeq.getX(pos-ephase);
                    col1 = aligSeq.getX(pos-ephase+1);
                    col2 = aligSeq.getX(pos-ephase+2);
                    //pwExon.println(emodelCodonExonic.logprob(col0,col1,col2,true));
                    //val = emodelExonic.get(0).logprob(col,true);
                    //pwInter.println(emodelIntergenic.logprob(col, true));
                    //pwInter.println(emodelCodonNoncoding.logprob(col0,col1,col2, true));
                    //pwIntron.println(emodelCodonNoncoding.logprob(col0,col1,col2,true));
                    String val = "[" + emodelCodonExonic.logprob(col0,col1,col2,true) + ";" + emodelIntronic.logprob(col0,true) + ";" + emodelIntronic.logprob(col1,true) + ";" + emodelIntronic.logprob(col2,true) + "]";
                    /*double value = emodelCodonExonic.logprob(col0,col1,col2,true) - (emodelIntronic.logprob(col0,true) + 
                                                                                    emodelIntronic.logprob(col1,true) + 
                                                                                    emodelIntronic.logprob(col2,true));*/
                    pwExon.println(val);
                    pwExon.println(val);
                    pwExon.println(val);
                    
                    //pwExon2.println((emodelCodonExonic.logprob(col0,col1,col2,false) - emodelCodonNoncoding.logprob(col0,col1,col2,false))/3);
                    /*pwExon.println(emodelCodonExonic.logprob(col0,col1,col2,false)*-2.499858e-03 + ((col0==null?0:emodelIntronic.logprob(col0,false)) + 
                                                                                                        (col1==null?0:emodelIntronic.logprob(col1,false)) + 
                                                                                                        (col2==null?0:emodelIntronic.logprob(col2,false)))*4.687246e-04);
                    
                     */break;
                case 0:
                case 4:
                case 5:
                case 6:
                    colm2 = pos>1?aligSeq.getX(pos-2):null;
                    colm1 = pos>0?aligSeq.getX(pos-1):null;
                    col0 = aligSeq.getX(pos);
                    col1 = pos<aligSeq.length()-1?aligSeq.getX(pos+1):null;
                    col2 = pos<aligSeq.length()-2?aligSeq.getX(pos+2):null;
                    //pwExon.println(emodelCodonExonic.logprob(col0,col1,col2,true));
                    //pwInter.println(emodelCodonNoncoding.logprob(col0,col1,col2, true));
                    //pwIntron.println(emodelCodonNoncoding.logprob(col0,col1,col2,true));
                    //pwExon.println(emodelCodonExonic.logprob(col0,col1,col2,true));
                                /*pwExon.println((emodelCodonExonic.logprob(col0,col1,col2,true)-((col0==null?0:emodelIntronic.logprob(col0,true)) + 
                                                               (col1==null?0:emodelIntronic.logprob(col1,true)) + 
                                                               (col2==null?0:emodelIntronic.logprob(col2,true))))/3);*/
                    //pwExon2.println((emodelCodonExonic.logprob(col0,col1,col2,false) - emodelCodonNoncoding.logprob(col0,col1,col2,false))/3);
                    /*pwExon.println(emodelCodonExonic.logprob(col0,col1,col2,false)*-2.499858e-03 + ((col0==null?0:emodelIntronic.logprob(col0,false)) + 
                                                                                                        (col1==null?0:emodelIntronic.logprob(col1,false)) + 
                                                                                                        (col2==null?0:emodelIntronic.logprob(col2,false)))*4.687246e-04);
                    
                    */break;
                case 7:
                case 8:
                case 9:
                    ephase = ((-pos + state + 1) % 3 + 3) % 3; //((pos-(state-1))%3 +3)%3;
                    if(ephase != 0) break;
                    col0 = aligSeq.getX(pos+ephase-2);
                    col1 = aligSeq.getX(pos+ephase-1);
                    col2 = aligSeq.getX(pos+ephase);
                    //pwExon.println(emodelCodonExonic.logprobRC(col0,col1,col2,true));
                    //val = emodelExonic.get(0).logprob(col,true);
                    //pwInter.println(emodelIntergenic.logprob(col, true));
                    //pwInter.println(emodelCodonNoncoding.logprobRC(col0,col1,col2, true));
                    //pwIntron.println(emodelCodonNoncoding.logprobRC(col0,col1,col2,true));
                        //pwExon2.println((emodelCodonExonic.logprob(col0,col1,col2,false) - emodelCodonNoncoding.logprobRC(col0,col1,col2,false))/3);
                    /*pwExon.println(emodelCodonExonic.logprobRC(col0,col1,col2,false)*-2.499858e-03 + ((col0==null?0:emodelIntronic.logprobRC(col0,false)) + 
                                                                                                        (col1==null?0:emodelIntronic.logprobRC(col1,false)) + 
                                                                                                        (col2==null?0:emodelIntronic.logprobRC(col2,false)))*4.687246e-04);
                    */break;
                case 10:
                case 11:
                case 12:
                    colm2 = pos>1?aligSeq.getX(pos-2):null;
                    colm1 = pos>0?aligSeq.getX(pos-1):null;
                    col0 = aligSeq.getX(pos);
                    col1 = pos<aligSeq.length()-1?aligSeq.getX(pos+1):null;
                    col2 = pos<aligSeq.length()-2?aligSeq.getX(pos+2):null;
                    //pwInter.println(emodelCodonNoncoding.logprob(col0,col1,col2, true));
                    //pwIntron.println(emodelCodonNoncoding.logprob(col0,col1,col2,true));
                    //pwExon.println(emodelCodonExonic.logprobRC(col0,col1,col2,true));
                    pwExon.println(
                    Math.max( 
                        emodelCodonExonic.logprobRC(colm2,colm1,col0,true)- ((colm2==null?0:emodelIntronic.logprobRC(colm2,true)) + 
                                                                           (colm1==null?0:emodelIntronic.logprobRC(colm1,true)) + 
                                                                           (col0==null?0:emodelIntronic.logprobRC(col0,true))),
                        Math.max( emodelCodonExonic.logprobRC(colm1,col0,col1,true)-((colm1==null?0:emodelIntronic.logprobRC(colm1,true)) + 
                                                                           (col0==null?0:emodelIntronic.logprobRC(col0,true)) + 
                                                                           (col1==null?0:emodelIntronic.logprobRC(col1,true))),
                                emodelCodonExonic.logprobRC(col0,col1,col2,true)-(col0==null?0:emodelIntronic.logprobRC(col0,true) + 
                                                                           (col1==null?0:emodelIntronic.logprobRC(col1,true)) + 
                                                                           (col2==null?0:emodelIntronic.logprobRC(col2,true))))));
                    /*pwExon.println((emodelCodonExonic.logprobRC(col0,col1,col2,true)-(col0==null?0:emodelIntronic.logprobRC(col0,true) + 
                                                               (col1==null?0:emodelIntronic.logprobRC(col1,true)) + 
                                                               (col2==null?0:emodelIntronic.logprobRC(col2,true))))/3);*/
                    //pwExon2.println((emodelCodonExonic.logprob(col0,col1,col2,false) - emodelCodonNoncoding.logprobRC(col0,col1,col2,false))/3);
                    /*pwExon.println(emodelCodonExonic.logprobRC(col0,col1,col2,false)*-2.499858e-03 + ((col0==null?0:emodelIntronic.logprobRC(col0,false)) + 
                                                                                                        (col1==null?0:emodelIntronic.logprobRC(col1,false)) + 
                                                                                                        (col2==null?0:emodelIntronic.logprobRC(col2,false)))*4.687246e-04);
                    */break;
                default:
                    Assert.a(false);
                }
            }
        }
        
        //pwInter.close();
        //pwIntron.close();
        pwExon.close();
        //pwExon2.close();
        
    }
    
    private void printCodonValues3(List<? extends TrainingSequence<?>> t) throws FileNotFoundException {
        //PrintWriter pwInter = new PrintWriter(new FileOutputStream(intergenicOutputFile));
        //PrintWriter pwIntron = new PrintWriter(new FileOutputStream(intronOutputFile));
        PrintWriter pwExon = new PrintWriter(new FileOutputStream(sourceDir+"/codon3.values"));
        //PrintWriter pwExon2 = new PrintWriter(new FileOutputStream(exonOutputFile+".codon2"));
        int ix = 0;
        for(TrainingSequence<?> seq:t){
            System.out.println(++ix);
            String name = ((NameInputSequence)seq.getInputSequence().getComponent(nameComponent)).getName();
            pwExon.println(">"+name);
            InputSequence<? extends MultipleAlignmentColumn> aligSeq = (InputSequence<? extends MultipleAlignmentColumn>)seq.getComponent(multipleAlnComponent);
            for(int pos=0;pos<aligSeq.length();pos++){
                
                int state = seq.getY(pos);
                int ephase;
                
                MultipleAlignmentColumn col0 = null;
                MultipleAlignmentColumn col1 = null;
                MultipleAlignmentColumn col2 = null;
                switch (state) {
                case 1:
                case 2:
                case 3:
                    ephase = ((pos - state + 1) % 3 + 3) % 3; //((pos-(state-1))%3 +3)%3;
                    
                    col0 = aligSeq.getX(pos-ephase);
                    col1 = aligSeq.getX(pos-ephase+1);
                    col2 = aligSeq.getX(pos-ephase+2);
                    //pwExon.println(emodelCodonExonic.logprob(col0,col1,col2,true));
                    //val = emodelExonic.get(0).logprob(col,true);
                    //pwInter.println(emodelIntergenic.logprob(col, true));
                    //pwInter.println(emodelCodonNoncoding.logprob(col0,col1,col2, true));
                    //pwIntron.println(emodelCodonNoncoding.logprob(col0,col1,col2,true));
                    //String val = "[" + emodelCodonExonic.logprob(col0,col1,col2,true) + ";" + emodelIntronic.logprob(col0,true) + ";" + emodelIntronic.logprob(col1,true) + ";" + emodelIntronic.logprob(col2,true) + "]";
                    double value = emodelCodonExonic.logprob(col0,col1,col2,true) - emodelCodonNoncoding.logprob(col0,col1,col2,true);
                    pwExon.println(value);
                    //pwExon.println(val);
                    //pwExon.println(val);
                    
                    //pwExon2.println((emodelCodonExonic.logprob(col0,col1,col2,false) - emodelCodonNoncoding.logprob(col0,col1,col2,false))/3);
                    /*pwExon.println(emodelCodonExonic.logprob(col0,col1,col2,false)*-2.499858e-03 + ((col0==null?0:emodelIntronic.logprob(col0,false)) + 
                                                                                                        (col1==null?0:emodelIntronic.logprob(col1,false)) + 
                                                                                                        (col2==null?0:emodelIntronic.logprob(col2,false)))*4.687246e-04);
                    
                     */break;
                case 0:
                case 4:
                case 5:
                case 6:
                    
                    col0 = aligSeq.getX(pos);
                    col1 = pos<aligSeq.length()-1?aligSeq.getX(pos+1):null;
                    col2 = pos<aligSeq.length()-2?aligSeq.getX(pos+2):null;
                    value = emodelCodonExonic.logprob(col0,col1,col2,true) - emodelCodonNoncoding.logprob(col0,col1,col2,true);
                    pwExon.println(value);
                    break;
                case 7:
                case 8:
                case 9:
                    ephase = ((-pos + state + 1) % 3 + 3) % 3; //((pos-(state-1))%3 +3)%3;
                    col0 = aligSeq.getX(pos+ephase-2);
                    col1 = aligSeq.getX(pos+ephase-1);
                    col2 = aligSeq.getX(pos+ephase);
                    value = emodelCodonExonic.logprobRC(col0,col1,col2,true) - emodelCodonNoncoding.logprobRC(col0,col1,col2,true);
                    pwExon.println(value);
                    break;
                case 10:
                case 11:
                case 12:
                    
                    col0 = aligSeq.getX(pos);
                    col1 = pos<aligSeq.length()-1?aligSeq.getX(pos+1):null;
                    col2 = pos<aligSeq.length()-2?aligSeq.getX(pos+2):null;
                    value = emodelCodonExonic.logprobRC(col0,col1,col2,true) - emodelCodonNoncoding.logprobRC(col0,col1,col2,true);
                    pwExon.println(value);
                    break;
                default:
                    Assert.a(false);
                }
            }
        }
        
        //pwInter.close();
        //pwIntron.close();
        pwExon.close();
        //pwExon2.close();
        
    }
                            
}
