ReversibleNucleotide.java

package statalign.model.subst.plugins;

import java.io.IOException;

import statalign.base.Utils;

/**
 * This class implements the general time reversible parameter model.
 * 
 * @author lyngsoe
 */
public class ReversibleNucleotide extends NucleotideModel{

	public static final String menuName = "General Reversible (Nucl.)";

	//static final double span = 0.1; // largest parameter change when resampling

	/**
	 * This constructor uses numerical methods to diagonalise the rate
	 * matrix and to set diagonalising matrices.
	 * @throws IOException
	 */
	public ReversibleNucleotide() throws IOException{
		/* Initialise NucleotideModel part of object */
		super();
		/* Initialise Reversible specific part of object */
		params = new double[8];
		oldparams = new double[8];
		/* Set all frequencies to 1/4 */
		params[0] = params[1] = params[2] = 0.25;
		/* Set all rate ratios to 1 */
		params[3] = params[4] = params[5] = params[6] = params[7] = 1.0;
		/* Set diagonalising matrices and equilibrium frequencies */
		setDiagonal();
	}

	void setDiagonal(){
		/* Set equilibrium frequencies */
		e[alphabet.length - 1] = 1.0;
		for (int i = 0; i < alphabet.length - 1; i++){
			e[i] = params[i];
			e[alphabet.length - 1] -= params[i];
		}
		/* Set diagonalising matrices */
		// Start by constructing rate matrix
		double Q[][] = new double[alphabet.length][alphabet.length];
		for (int i = 0; i < alphabet.length; i++)
			Q[i][i] = 0.0;
		for (int i = 1; i < alphabet.length; i++)
			for (int j = 0; j < i; j++){
				if (j < 2){
					Q[i][j] = e[j] * params[2 + 2 * j + i];
					Q[j][i] = e[i] * params[2 + 2 * j + i];
				}
				else{
					Q[i][j] = e[j];
					Q[j][i] = e[i];
				}
				Q[i][i] -= Q[i][j];
				Q[j][j] -= Q[j][i];
			}
		NucleotideModelUtils.numericalDiagonalisation(this, Q);
	}

	/**
	 * Prints the parameters of the model
	 */
	@Override
	public String print(){
		String s = "";
		for(int i = 3; i < params.length; i++){
			s+="\trho_"+(i-2)+"\t"+params[i];
		}
		return "A\t"+e[0]+"\tC\t"+e[1]+"\tG\t"+e[2]+"\tT\t"+e[3]+s;
	}

	/**
	 * This function implements a proposal for new parameter values.
	 * Returns with the logarithm of the Metropolis-Hastings ratio.
	 */
	@Override
	public double sampleParameter() {
		double x;

		if (Utils.generator.nextInt(8) < 3){
			/* Resample frequency parameter */
			x = NucleotideModelUtils.sampleFrequencies(e, span / 2.0);
			/* Save current parameters as old parameters, store new
			 * frequency parameters, and ensure frequencies sum to one.
			 */
			e[alphabet.length - 1] = 1.0;
			for (int i = 0; i < params.length; i++){
				oldparams[i] = params[i];
				if (i < alphabet.length - 1){
					params[i] = e[i];
					e[alphabet.length - 1] -= e[i];
				}
			}
		}
		else{
			/* Resample rate ratio parameter */
			double r[] = new double[5];
			for (int i = 0; i < 5; i++)
				r[i] = params[3 + i];
			x = NucleotideModelUtils.sampleRates(r, span / 2.0);
			/* save all params (bugfix), we thank William Majoros for reporting */
            for (int i = 0; i < params.length; i++)
            	oldparams[i] = params[i];
            for (int i = 0; i < 5; i++)
				params[3 + i] = r[i];

		}
		setDiagonal();
		return x;
	}
}