package fmph.features.supporting.phylogenetic;

import calhoun.util.Assert;

import fmph.features.util.GeneticCodes;

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

public class GoldmanYang94 extends CodonSubstitutionModel {
    private static final long serialVersionUID = -1555350959092722736L;
    private static final Log log = LogFactory.getLog(GoldmanYang94.class);
    // Model for nucleotide substitution proposed by Kimura in 1980.
    // Models the overall rate of nucleotide substitution and the difference
    // between rate of transitions and transversions.

    double[] piCodon = new double[64];
    double kappa;
    double omega;
    int code;


    public GoldmanYang94(double[] parms) {
        setParameters(parms);
    }


    @Override
    public String getEvolutionaryModelName() {
        return "GoldmanYang94_codon_evolution_model";
    }

    @Override
    public double[] getParameters() {
        double[] ret = new double[6];

        return ret;
    }


    @Override
    public void setParameters(double[] parms) {
        Assert.a(parms.length == 66);
        System.arraycopy(parms, 0, piCodon, 0, 63);

        double sum = 0;
        for (int i = 0; i < 63; i++) {
            sum += parms[i];
        }
        piCodon[63] = 1 - sum;

        kappa = parms[63];
        omega = parms[64];
        code = (int)parms[65];

        String codeTable = GeneticCodes.getCode(code)[1];

        double[][] X = new double[64][64];

        int T = 0, C = 1, A = 2, G = 3;
        for (int i = 0; i < 64; i++) {
            int ns1 = i / 16;
            int ns2 = (i % 16) / 4;
            int ns3 = (i % 4);
            char aas = codeTable.charAt(i);
            for (int j = 0; j < 64; j++) {
                if (i == j)
                    continue;
                int nd1 = j / 16;
                int nd2 = (j % 16) / 4;
                int nd3 = (j % 4);
                char aad = codeTable.charAt(j);
                /* TCAG */
                boolean synonymous = aas == aad;
                int ns = -1;
                int nd = -1;

                int mismatches = 0;
                if (ns1 != nd1) {
                    ns = ns1;
                    nd = nd1;
                    mismatches++;

                }
                if (ns2 != nd2) {
                    ns = ns2;
                    nd = nd2;
                    mismatches++;
                }
                if (ns3 != nd3) {
                    ns = ns3;
                    nd = nd3;
                    mismatches++;
                }

                if (mismatches == 1) {
                    if (ns + nd == 1 || ns + nd == 5) { //transition
                        if (synonymous) {
                            X[i][j] = kappa * piCodon[j];
                        } else {
                            X[i][j] = omega * kappa * piCodon[j];
                        }
                    } else { //transversion
                        if (synonymous) {
                            X[i][j] = piCodon[j];
                        } else {
                            X[i][j] = omega * piCodon[j];
                        }
                    }
                } else {
                    X[i][j] = 0;
                }

            }
        }

        for (int i = 0; i < 64; i++) {
            for (int j = 0; j < 64; j++) {
                if (j != i) {
                    X[i][i] -= X[i][j];
                }
            }
        }

        R.assign(X);
    }
    
    

    @Override
    public void summarize() {
        log.debug("GoldmanYang94 codon substitution rate matrix; kappa=" + kappa + "  omega=" + omega + "  code=" + code);
        
        StringBuffer buff = new StringBuffer();
        for (int i = 0; i < 64; i++) {
           buff.append("\t"+GeneticCodes.getCodon(i));
        }
        buff.append("\n");
        
        
        double[][] X = R.toArray();
        
        for (int i = 0; i < 64; i++) {
            buff.append(GeneticCodes.getCodon(i));
            for (int j = 0; j < 64; j++) {
                buff.append("\t"+X[i][j]);
            }
            buff.append("\n");
        }
        log.debug(buff);
    }
}
