<commit_before>package beast.evolution.substitutionmodel;

import beast.core.Description;
import beast.core.Input;
import beast.core.Input.Validate;
import beast.core.parameter.RealParameter;
import beast.core.util.Log;
import beast.evolution.datatype.DataType;
import beast.evolution.datatype.TwoStateCovarion;


/**
 * <p/>
 * a	the rate of the slow rate class
 * 1	the rate of the fast rate class
 * p0	the equilibrium frequency of zero states
 * p1	1 - p0, the equilibrium frequency of one states
 * f0	the equilibrium frequency of slow rate class
 * f1	1 - f0, the equilibrium frequency of fast rate class
 * s, s1, s2 the rate of switching
 * <p/>
 * then the (unnormalized) instantaneous rate matrix (unnormalized Q) should be (depending on mode)
 * <p/>
 * 
 * mode = BEAST -- using classic BEAST implementation, reversible iff hfrequencies = (0.5, 0.5)
 * FLAGS: reversible = false, TSParameterisation = false
 * 
 * [ -(a*p1)-s ,   a*p1    ,    s   ,   0   ]
 * [   a*p0    , -(a*p0)-s ,    0   ,   s   ]
 * [    s      ,     0     ,  -p1-s ,  p1   ]
 * [    0      ,     s     ,    p0  , -p0-s ]
 *
 * equilibrium frequencies
 * [ p0 * f0, p1, * f0, p0 * f1, p1, * f1 ]
 *
 * mode = REVERSIBLE -- brings in hfrequencies in rate matrix
 * reversible = true, TSParameterisation = false
 * [ - , a , s , 0 ]
 * [ a , - , 0 , s ]
 * [ s , 0 , - , 1 ]
 * [ 0 , s , 1 , - ]
 * 
 * which with frequencies becomes
 * 
 * [ -(a*p1*f0)-s*f0 ,   a*p1*f0       ,    s*f0      ,   0         ]
 * [   a*p0*f0       , -(a*p0*f0)-s*f0 ,    0         ,   s*f0      ]
 * [    s*f1         ,     0           ,  -p1*f1-s*f1 ,  p1*f1      ]
 * [    0            ,     s*f1        ,    p0*f1     , -p0*f1-s*f1 ]
 * 
 * equilibrium frequencies
 * [ p0 * f0, p1, * f0, p0 * f1, p1, * f1 ]
 * 
 * mode = TUFFLEYSTEEL uses alternative parameterisation: hfrequencies is ignored, and switch parameter is set to dimension = 2
 * [ -(a*p1)-s1 ,   a*p1     ,    s1   ,   0    ]
 * [   a*p0     , -(a*p0)-s1 ,    0    ,   s1   ]
 * [    s2      ,     0      ,  -p1-s2 ,  p1    ]
 * [    0       ,     s2     ,    p0   , -p0-s2 ]
 *
 * equilibrium frequencies
 * [ f0 * s2/(s1+s2), f1, * s2/(s1+s2), f0 * s1/(s1+s2), f1, * s1/(s1+s2) ]
 *
 *
 * Note: to use Tuffley & Steel's methods, set a = 0.
 */
@Description("Covarion model for Binary data")
public class BinaryCovarion extends GeneralSubstitutionModel {
    public Input<RealParameter> alphaInput = new Input<RealParameter>("alpha", "the rate of evolution in slow mode", Validate.REQUIRED);
    public Input<RealParameter> switchRateInput = new Input<RealParameter>("switchRate", "the rate of flipping between slow and fast modes", Validate.REQUIRED);
    public Input<RealParameter> frequenciesInput = new Input<RealParameter>("vfrequencies", "the frequencies of the visible states", Validate.REQUIRED);
    public Input<RealParameter> hfrequenciesInput = new Input<RealParameter>("hfrequencies", "the frequencies of the hidden rates");

    public enum MODE {BEAST, REVERSIBLE, TUFFLEYSTEEL};
	public Input<MODE> modeInput = new Input<>("mode","one of BEAST, REVERSIBLE, TUFFLESTEEL "
			+ "BEAST = implementation as in BEAST 1 "
			+ "REVERSIBLE = like BEAST 1 implementation, but using frequencies to make it reversible "
			+ "TUFFLEYSTEEL = Tuffley & Steel (1996) impementation (no rates for ", MODE.BEAST,MODE.values());

    private RealParameter alpha;
    private RealParameter switchRate;
    private RealParameter frequencies;
    private RealParameter hiddenFrequencies;

    protected double[][] unnormalizedQ;
    protected double[][] storedUnnormalizedQ;
    int stateCount;
    MODE mode  = modeInput.get();

    public BinaryCovarion() {
        ratesInput.setRule(Validate.OPTIONAL);
        frequenciesInput.setRule(Validate.OPTIONAL);
    }

    @Override
    public void initAndValidate() throws Exception {
        alpha = alphaInput.get();
        switchRate = switchRateInput.get();
        frequencies = frequenciesInput.get();
        hiddenFrequencies = hfrequenciesInput.get();

        
        if (mode.equals(MODE.BEAST) || mode.equals(MODE.REVERSIBLE)) {
        	if (alpha.getDimension() != 1) {
        		throw new Exception("alpha should have dimension 1");
        	}
        } else {
        	if (alpha.getDimension() != 2) {
        		throw new Exception("alpha should have dimension 2");
        	}
        }
        if (switchRate.getDimension() != 1) {
            throw new Exception("switchRate should have dimension 1");
        }
        if (frequencies.getDimension() != 2) {
            throw new Exception("frequencies should have dimension 2");
        }
        if (mode.equals(MODE.BEAST) || mode.equals(MODE.REVERSIBLE)) {
        	if (hfrequenciesInput.get() == null) {
        		throw new Exception("hiddenFrequenciesshould should be specified");
        	}
            if (hiddenFrequencies.getDimension() != 2) {
                throw new Exception("hiddenFrequenciesshould have dimension 2");
            }
        } else {
        	if (hfrequenciesInput.get() != null) {
        		Log.warning.println("WARNING: hfrequencies is specified, but the BinaryCovarion model ignores it.");
        	}
        }<commit_msg>issue #398 BinaryCovarion sanity check for Tuffely & Steel mode fixed<commit_after>package beast.evolution.substitutionmodel;

import beast.core.Description;
import beast.core.Input;
import beast.core.Input.Validate;
import beast.core.parameter.RealParameter;
import beast.core.util.Log;
import beast.evolution.datatype.DataType;
import beast.evolution.datatype.TwoStateCovarion;


/**
 * <p/>
 * a	the rate of the slow rate class
 * 1	the rate of the fast rate class
 * p0	the equilibrium frequency of zero states
 * p1	1 - p0, the equilibrium frequency of one states
 * f0	the equilibrium frequency of slow rate class
 * f1	1 - f0, the equilibrium frequency of fast rate class
 * s, s1, s2 the rate of switching
 * <p/>
 * then the (unnormalized) instantaneous rate matrix (unnormalized Q) should be (depending on mode)
 * <p/>
 * 
 * mode = BEAST -- using classic BEAST implementation, reversible iff hfrequencies = (0.5, 0.5)
 * FLAGS: reversible = false, TSParameterisation = false
 * 
 * [ -(a*p1)-s ,   a*p1    ,    s   ,   0   ]
 * [   a*p0    , -(a*p0)-s ,    0   ,   s   ]
 * [    s      ,     0     ,  -p1-s ,  p1   ]
 * [    0      ,     s     ,    p0  , -p0-s ]
 *
 * equilibrium frequencies
 * [ p0 * f0, p1, * f0, p0 * f1, p1, * f1 ]
 *
 * mode = REVERSIBLE -- brings in hfrequencies in rate matrix
 * reversible = true, TSParameterisation = false
 * [ - , a , s , 0 ]
 * [ a , - , 0 , s ]
 * [ s , 0 , - , 1 ]
 * [ 0 , s , 1 , - ]
 * 
 * which with frequencies becomes
 * 
 * [ -(a*p1*f0)-s*f0 ,   a*p1*f0       ,    s*f0      ,   0         ]
 * [   a*p0*f0       , -(a*p0*f0)-s*f0 ,    0         ,   s*f0      ]
 * [    s*f1         ,     0           ,  -p1*f1-s*f1 ,  p1*f1      ]
 * [    0            ,     s*f1        ,    p0*f1     , -p0*f1-s*f1 ]
 * 
 * equilibrium frequencies
 * [ p0 * f0, p1, * f0, p0 * f1, p1, * f1 ]
 * 
 * mode = TUFFLEYSTEEL uses alternative parameterisation: hfrequencies is ignored, and switch parameter is set to dimension = 2
 * [ -(a*p1)-s1 ,   a*p1     ,    s1   ,   0    ]
 * [   a*p0     , -(a*p0)-s1 ,    0    ,   s1   ]
 * [    s2      ,     0      ,  -p1-s2 ,  p1    ]
 * [    0       ,     s2     ,    p0   , -p0-s2 ]
 *
 * equilibrium frequencies
 * [ f0 * s2/(s1+s2), f1, * s2/(s1+s2), f0 * s1/(s1+s2), f1, * s1/(s1+s2) ]
 *
 *
 * Note: to use Tuffley & Steel's methods, set a = 0.
 */
@Description("Covarion model for Binary data")
public class BinaryCovarion extends GeneralSubstitutionModel {
    public Input<RealParameter> alphaInput = new Input<RealParameter>("alpha", "the rate of evolution in slow mode", Validate.REQUIRED);
    public Input<RealParameter> switchRateInput = new Input<RealParameter>("switchRate", "the rate of flipping between slow and fast modes", Validate.REQUIRED);
    public Input<RealParameter> frequenciesInput = new Input<RealParameter>("vfrequencies", "the frequencies of the visible states", Validate.REQUIRED);
    public Input<RealParameter> hfrequenciesInput = new Input<RealParameter>("hfrequencies", "the frequencies of the hidden rates");

    public enum MODE {BEAST, REVERSIBLE, TUFFLEYSTEEL};
	public Input<MODE> modeInput = new Input<>("mode","one of BEAST, REVERSIBLE, TUFFLESTEEL "
			+ "BEAST = implementation as in BEAST 1 "
			+ "REVERSIBLE = like BEAST 1 implementation, but using frequencies to make it reversible "
			+ "TUFFLEYSTEEL = Tuffley & Steel (1996) impementation (no rates for ", MODE.BEAST,MODE.values());

    private RealParameter alpha;
    private RealParameter switchRate;
    private RealParameter frequencies;
    private RealParameter hiddenFrequencies;

    protected double[][] unnormalizedQ;
    protected double[][] storedUnnormalizedQ;
    int stateCount;
    MODE mode  = modeInput.get();

    public BinaryCovarion() {
        ratesInput.setRule(Validate.OPTIONAL);
        frequenciesInput.setRule(Validate.OPTIONAL);
    }

    @Override
    public void initAndValidate() throws Exception {
        alpha = alphaInput.get();
        switchRate = switchRateInput.get();
        frequencies = frequenciesInput.get();
        hiddenFrequencies = hfrequenciesInput.get();

        
        if (mode.equals(MODE.BEAST) || mode.equals(MODE.REVERSIBLE)) {
        	if (switchRate.getDimension() != 1) {
        		throw new Exception("switchRate should have dimension 1");
        	}
        } else {
        	if (switchRate.getDimension() != 2) {
        		throw new Exception("switchRate should have dimension 2");
        	}
        }
        if (alpha.getDimension() != 1) {
            throw new Exception("alpha should have dimension 1");
        }
        if (frequencies.getDimension() != 2) {
            throw new Exception("frequencies should have dimension 2");
        }
        if (mode.equals(MODE.BEAST) || mode.equals(MODE.REVERSIBLE)) {
        	if (hfrequenciesInput.get() == null) {
        		throw new Exception("hiddenFrequenciesshould should be specified");
        	}
            if (hiddenFrequencies.getDimension() != 2) {
                throw new Exception("hiddenFrequenciesshould have dimension 2");
            }
        } else {
        	if (hfrequenciesInput.get() != null) {
        		Log.warning.println("WARNING: hfrequencies is specified, but the BinaryCovarion model ignores it.");
        	}
        }
{'fraction_non_alphanumeric': 0.1215648854961832, 'fraction_numerical': 0.025286259541984733, 'mean_word_length': 2.4762852404643447, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 29, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '977815', 'n_tokens_mistral': 3538, 'n_tokens_neox': 3414, 'n_words': 927}