package fmph.features.mitochondrion1;

import calhoun.analysis.crf.AbstractFeatureManager;
import calhoun.analysis.crf.CacheStrategySpec;
import calhoun.analysis.crf.CacheStrategySpec.CacheStrategy;
import calhoun.analysis.crf.FeatureList;
import calhoun.analysis.crf.FeatureManagerEdge;
import calhoun.analysis.crf.ModelManager;
import calhoun.analysis.crf.features.supporting.phylogenetic.PhylogeneticTreeFelsensteinOrder;
import calhoun.analysis.crf.io.InputSequence;
import calhoun.analysis.crf.io.MultipleAlignmentInputSequence;
import calhoun.analysis.crf.io.MultipleAlignmentInputSequence.MultipleAlignmentColumn;
import calhoun.analysis.crf.io.TrainingSequence;

import fmph.features.supporting.phylogenetic.CodonEvolutionaryModel;
import fmph.features.supporting.phylogenetic.GoldmanYang94;

import java.util.List;

import org.apache.commons.logging.Log;
import org.apache.commons.logging.LogFactory;

public class CodonBoundaryLogprobMitochondrion1 extends AbstractFeatureManager<MultipleAlignmentColumn> implements FeatureManagerEdge<MultipleAlignmentColumn> {
    private static final long serialVersionUID = -7659288739348604129L;
    private static final Log log = LogFactory.getLog(CodonBoundaryLogprobMitochondrion1.class);


    int startIx; // The index of the first feature managed by this FeatureManager
    ModelManager model;
    
    boolean[] intergenicState;
    boolean[] plusExonState;
    boolean[] minusExonState;
    boolean[] plusIntronState;
    boolean[] minusIntronState;     

    // one model for a column of aligned sequence in intergenic region
    CodonEvolutionaryModel emodelExonic;
    
    double codonEdgeCoef = 1;

    
    double[] emodelExonicParams = {
        //codon freq
    0.06870156,0.01163551,0.06247777,0.00667686,
    0.03160051,0.00535196,0.02873777,0.00307114,
    0.03577897,0.00605133,0.00001000,0.00001000,
    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.65390,//kappa
    0.08634,//omega
    4 //icode
        } ;

    //static KmerHasher hforward = new KmerHasher(KmerHasher.ACGTother, 1); // a character hasher for forward strand
    //static KmerHasher hbackward = new KmerHasher(KmerHasher.ACGTotherRC, 1); // a character hasher for reverse strand

    ///////////////////////////////////////////////////////////////////////////////

    public CodonBoundaryLogprobMitochondrion1() {
    } // a constructor with no arguments

    public int getNumFeatures() { // there is exactly one feature
        return 2;
    }

    public String getFeatureName(int featureIndex) {
        String[] names = {"Exonic boundary codon phylogeny 2 ambigious columns","Exonic boundary codon phylogeny 1 ambigious column" };
        return names[featureIndex - startIx];
    }

    public void evaluateEdge(InputSequence<? extends MultipleAlignmentColumn> seq, int pos, int prevState, int state, FeatureList result) {
        if(pos < 2 || pos > seq.length()-3) return;
        // Evauluates boundary codons in exon around an intron
        MultipleAlignmentInputSequence.MultipleAlignmentColumn col0 = null;
        MultipleAlignmentInputSequence.MultipleAlignmentColumn col1 = null;
        MultipleAlignmentInputSequence.MultipleAlignmentColumn col2 = null;


        if (plusIntronState[prevState] && plusExonState[state] ) { 
            switch (prevState) {
            case 4:
                break; // splitting occured between codons
            case 5:
                col0 = null;
                col1 = null;
                col2 = seq.getX(pos);
                result.addFeature(startIx, emodelExonic.logprob(col0, col1, col2, true));
                break;
            case 6:
                col0 = null;
                col1 = seq.getX(pos);
                col2 = seq.getX(pos + 1);
                result.addFeature(startIx+1, emodelExonic.logprob(col0, col1, col2, true));
                break;
            }
        }
        
        if (minusExonState[prevState] && minusIntronState[state]) {
            switch (state) {
            case 10:
                break; // splitting occured between codons
            case 11:
                col0 = seq.getX(pos - 1);
                col1 = null;
                col2 = null;
                result.addFeature(startIx, emodelExonic.logprobRC(col0, col1, col2, true));
                break;
            case 12:
                col0 = seq.getX(pos - 2);
                col1 = seq.getX(pos - 1);
                col2 = null;
                result.addFeature(startIx+1, emodelExonic.logprobRC(col0, col1, col2, true));
                break;
            }
        }
        
        if (plusExonState[prevState] && plusIntronState[state]) {
            switch (state) {
            case 4:
                break; // splitting occured between codons
            case 5:
                col0 = seq.getX(pos-2);
                col1 = seq.getX(pos-1);
                col2 = null;
                result.addFeature(startIx+1, emodelExonic.logprob(col0, col1, col2, true));
                break;
            case 6:
                col0 = seq.getX(pos-1);
                col1 = null;
                col2 = null;
                result.addFeature(startIx, emodelExonic.logprob(col0, col1, col2, true));
                break;
            }
        }
        if (plusIntronState[prevState] && plusExonState[state]) { 
            switch (prevState) {
            case 10:
                break; // splitting occured between codons
            case 11:
                col0 = seq.getX(pos);
                col1 = seq.getX(pos+1);
                col2 = null;
                result.addFeature(startIx+1, emodelExonic.logprobRC(col0, col1, col2, true));
                break;
            case 12:
                col0 = seq.getX(pos);
                col1 = null;
                col2 = null;
                result.addFeature(startIx, emodelExonic.logprobRC(col0, col1, col2, true));
                break;
            }
        }
    }


    public void train(int startingIndex, ModelManager modelInfo, final List<? extends TrainingSequence<? extends MultipleAlignmentColumn>> data) {
        startIx = startingIndex;
        model = modelInfo;

        final PhylogeneticTreeFelsensteinOrder felsOrder =
            ((MultipleAlignmentInputSequence.MultipleAlignmentColumn)data.get(0).getX(0)).getMultipleAlignment().getFelsensteinOrder();
        double[] codonExonicPi = new double[64];
        System.arraycopy(emodelExonicParams, 0, codonExonicPi, 0, 63);
        
        double sum = 0;
        for(int i=0;i<63;i++){
            sum += emodelExonicParams[i];
        }
        codonExonicPi[63] = 1 - sum;

        emodelExonic = new CodonEvolutionaryModel(felsOrder, codonExonicPi, new GoldmanYang94(emodelExonicParams), codonEdgeCoef);
        log.debug("Trained codon evolutionary bundary feature");
        emodelExonic.summarize();
        
        int nStates = model.getNumStates();
        
        plusExonState = new boolean[nStates];
        plusExonState[model.getStateIndex("exon0")] = true;
        plusExonState[model.getStateIndex("exon1")] = true;
        plusExonState[model.getStateIndex("exon2")] = true;
        
        minusExonState = new boolean[nStates]; 
        minusExonState[model.getStateIndex("exon0m")] = true;
        minusExonState[model.getStateIndex("exon1m")] = true;
        minusExonState[model.getStateIndex("exon2m")] = true;           
        
        plusIntronState = new boolean[nStates];
        plusIntronState[model.getStateIndex("intron0")] = true;
        plusIntronState[model.getStateIndex("intron1")] = true;
        plusIntronState[model.getStateIndex("intron2")] = true;
        
        minusIntronState = new boolean[nStates]; 
        minusIntronState[model.getStateIndex("intron0m")] = true;
        minusIntronState[model.getStateIndex("intron1m")] = true;
        minusIntronState[model.getStateIndex("intron2m")] = true;       

        intergenicState = new boolean[nStates];
        intergenicState[model.getStateIndex("intergenic")] = true;
    }


    @Override
    public CacheStrategySpec getCacheStrategy() {
        CacheStrategySpec css = new CacheStrategySpec(CacheStrategy.UNSPECIFIED);
        return css;
    }

    public void setCodonEdgeCoef(double codonEdgeCoef) {
        this.codonEdgeCoef = codonEdgeCoef;
    }

    public double getCodonEdgeCoef() {
        return codonEdgeCoef;
    }
}
