<commit_before>/*
* File GeneralSubstitutionModel.java
*
* Copyright (C) 2010 <NAME> <EMAIL>
*
* This file is not copyright Remco! It is copied from BEAST 1.
*
* This file is part of BEAST2.
* See the NOTICE file distributed with this work for additional
* information regarding copyright ownership and licensing.
*
* BEAST is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as
* published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
*
*  BEAST 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 Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with BEAST; if not, write to the
* Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
* Boston, MA  02110-1301  USA
*/
package beast.evolution.substitutionmodel;


import beast.core.Description;
import beast.core.Input;
import beast.core.Valuable;
import beast.core.Input.Validate;
import beast.evolution.tree.Node;

@Description("Specifies transition probability matrix with no restrictions on the rates other " +
        "than that one of the is equal to one and the others are specified relative to " +
        "this unit rate. Works for any number of states.")
public class GeneralSubstitutionModel extends SubstitutionModel.Base {
    public Input<Valuable> m_rates =
            new Input<Valuable>("rates", "Rate parameter which defines the transition rate matrix. " +
            		"Only the off-diagonal entries need to be specified (diagonal makes row sum to zero in a " +
            		"rate matrix). Entry i specifies the rate of from i%n to floor(i/(n-1)) where " +
            		"n is the number of states.", Validate.REQUIRED);

    /** number of states **/
    int m_nStates;
    /** a square m_nStates x m_nStates matrix containing current rates  **/
    double [][] m_rateMatrix;
    
    @Override
    public void initAndValidate() throws Exception {
        updateMatrix = true;
        m_nStates = frequencies.get().getFreqs().length;
        if (m_rates.get().getDimension() != m_nStates * (m_nStates-1)) {
        	throw new Exception("Dimension of input 'rates' is " + m_rates.get().getDimension() + " but a " +
        			"rate matrix of dimension " + m_nStates + "x" + (m_nStates -1) + "=" + m_nStates * (m_nStates -1) + " was " +
        			"expected");
        }
        eigenSystem = new DefaultEigenSystem(m_nStates);
        m_rateMatrix = new double[m_nStates][m_nStates];
        relativeRates = new double[m_rates.get().getDimension()];
        storedRelativeRates = new double[m_rates.get().getDimension()];
        //q = new double[m_nStates][m_nStates];
    } // initAndValidate

    protected double[] relativeRates;
    protected double[] storedRelativeRates;

    protected EigenSystem eigenSystem;
    
    protected EigenDecomposition eigenDecomposition;
    private EigenDecomposition storedEigenDecomposition;

    protected boolean updateMatrix = true;
    private boolean storedUpdateMatrix = true;

    @Override
    public void getTransitionProbabilities(Node node, double fStartTime, double fEndTime, double fRate, double[] matrix) {
    	double distance = (fStartTime - fEndTime) * fRate;
    	
        int i, j, k;
        double temp;

        // this must be synchronized to avoid being called simultaneously by
        // two different likelihood threads - AJD
        synchronized (this) {
            if (updateMatrix) {
            	setupRelativeRates();
            	setupRateMatrix();
            	eigenDecomposition = eigenSystem.decomposeMatrix(m_rateMatrix);
            	updateMatrix = false;
            }
        }

        // is the following really necessary?
        // implemented a pool of iexp matrices to support multiple threads
        // without creating a new matrix each call. - AJD
        // a quick timing experiment shows no difference - RRB
        double[] iexp = new double[m_nStates * m_nStates];
        // Eigen vectors
        double[] Evec = eigenDecomposition.getEigenVectors();
        // inverse Eigen vectors
        double[] Ievc = eigenDecomposition.getInverseEigenVectors();
        // Eigen values
        double[] Eval = eigenDecomposition.getEigenValues();
        for (i = 0; i < m_nStates; i++) {
            temp = Math.exp(distance * Eval[i]);
            for (j = 0; j < m_nStates; j++) {
                iexp[i * m_nStates + j] = Ievc[i * m_nStates + j] * temp;
            }
        }

        int u = 0;
        for (i = 0; i < m_nStates; i++) {
            for (j = 0; j < m_nStates; j++) {
                temp = 0.0;
                for (k = 0; k < m_nStates; k++) {
                    temp += Evec[i * m_nStates + k] * iexp[k * m_nStates + j];
                }

                matrix[u] = Math.abs(temp);
                u++;
            }
        }
    } // getTransitionProbabilities

    
    protected void setupRelativeRates() {
    	Valuable rates = m_rates.get();
    	for (int i = 0; i < rates.getDimension(); i++) {
    		relativeRates[i] = rates.getArrayValue(i);
    	}
    }

    /** sets up rate matrix **/
    protected void setupRateMatrix() {
    	double [] fFreqs = frequencies.get().getFreqs();<commit_msg>Make rate matrix readable for general substitution model<commit_after>/*
* File GeneralSubstitutionModel.java
*
* Copyright (C) 2010 <NAME> <EMAIL>
*
* This file is not copyright Remco! It is copied from BEAST 1.
*
* This file is part of BEAST2.
* See the NOTICE file distributed with this work for additional
* information regarding copyright ownership and licensing.
*
* BEAST is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as
* published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
*
*  BEAST 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 Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with BEAST; if not, write to the
* Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
* Boston, MA  02110-1301  USA
*/
package beast.evolution.substitutionmodel;


import beast.core.Description;
import beast.core.Input;
import beast.core.Valuable;
import beast.core.Input.Validate;
import beast.evolution.tree.Node;

@Description("Specifies transition probability matrix with no restrictions on the rates other " +
        "than that one of the is equal to one and the others are specified relative to " +
        "this unit rate. Works for any number of states.")
public class GeneralSubstitutionModel extends SubstitutionModel.Base {
    public Input<Valuable> m_rates =
            new Input<Valuable>("rates", "Rate parameter which defines the transition rate matrix. " +
            		"Only the off-diagonal entries need to be specified (diagonal makes row sum to zero in a " +
            		"rate matrix). Entry i specifies the rate of from i%n to floor(i/(n-1)) where " +
            		"n is the number of states.", Validate.REQUIRED);

    /** number of states **/
    int m_nStates;
    /** a square m_nStates x m_nStates matrix containing current rates  **/
    double [][] m_rateMatrix;
    
    @Override
    public void initAndValidate() throws Exception {
        updateMatrix = true;
        m_nStates = frequencies.get().getFreqs().length;
        if (m_rates.get().getDimension() != m_nStates * (m_nStates-1)) {
        	throw new Exception("Dimension of input 'rates' is " + m_rates.get().getDimension() + " but a " +
        			"rate matrix of dimension " + m_nStates + "x" + (m_nStates -1) + "=" + m_nStates * (m_nStates -1) + " was " +
        			"expected");
        }
        eigenSystem = new DefaultEigenSystem(m_nStates);
        m_rateMatrix = new double[m_nStates][m_nStates];
        relativeRates = new double[m_rates.get().getDimension()];
        storedRelativeRates = new double[m_rates.get().getDimension()];
        //q = new double[m_nStates][m_nStates];
    } // initAndValidate

    protected double[] relativeRates;
    protected double[] storedRelativeRates;

    protected EigenSystem eigenSystem;
    
    protected EigenDecomposition eigenDecomposition;
    private EigenDecomposition storedEigenDecomposition;

    protected boolean updateMatrix = true;
    private boolean storedUpdateMatrix = true;

    @Override
    public void getTransitionProbabilities(Node node, double fStartTime, double fEndTime, double fRate, double[] matrix) {
    	double distance = (fStartTime - fEndTime) * fRate;
    	
        int i, j, k;
        double temp;

        // this must be synchronized to avoid being called simultaneously by
        // two different likelihood threads - AJD
        synchronized (this) {
            if (updateMatrix) {
            	setupRelativeRates();
            	setupRateMatrix();
            	eigenDecomposition = eigenSystem.decomposeMatrix(m_rateMatrix);
            	updateMatrix = false;
            }
        }

        // is the following really necessary?
        // implemented a pool of iexp matrices to support multiple threads
        // without creating a new matrix each call. - AJD
        // a quick timing experiment shows no difference - RRB
        double[] iexp = new double[m_nStates * m_nStates];
        // Eigen vectors
        double[] Evec = eigenDecomposition.getEigenVectors();
        // inverse Eigen vectors
        double[] Ievc = eigenDecomposition.getInverseEigenVectors();
        // Eigen values
        double[] Eval = eigenDecomposition.getEigenValues();
        for (i = 0; i < m_nStates; i++) {
            temp = Math.exp(distance * Eval[i]);
            for (j = 0; j < m_nStates; j++) {
                iexp[i * m_nStates + j] = Ievc[i * m_nStates + j] * temp;
            }
        }

        int u = 0;
        for (i = 0; i < m_nStates; i++) {
            for (j = 0; j < m_nStates; j++) {
                temp = 0.0;
                for (k = 0; k < m_nStates; k++) {
                    temp += Evec[i * m_nStates + k] * iexp[k * m_nStates + j];
                }

                matrix[u] = Math.abs(temp);
                u++;
            }
        }
    } // getTransitionProbabilities

    /** access to (copy of) rate matrix **/
    protected double [][] getRateMatrix() {return m_rateMatrix.clone();}
    
    protected void setupRelativeRates() {
    	Valuable rates = m_rates.get();
    	for (int i = 0; i < rates.getDimension(); i++) {
    		relativeRates[i] = rates.getArrayValue(i);
    	}
    }

    /** sets up rate matrix **/
    protected void setupRateMatrix() {
    	double [] fFreqs = frequencies.get().getFreqs();
{'fraction_non_alphanumeric': 0.09574660633484162, 'fraction_numerical': 0.005610859728506787, 'mean_word_length': 2.5385846942042907, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4154337', 'n_tokens_mistral': 3201, 'n_tokens_neox': 3074, 'n_words': 1079}