Examples with GammaSiteRateModel dr.evomodel.siteratemodel.GammaSiteRateModel used on opensource projects

Example 1 with GammaSiteRateModel

use of dr.evomodel.siteratemodel.GammaSiteRateModel in project beast-mcmc by beast-dev.

the class BeagleTreeLikelihood method main.

public static void main(String[] args) {
    try {
        MathUtils.setSeed(666);
        System.out.println("Test case 1: simulateOnePartition");
        int sequenceLength = 1000;
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        // create tree
        NewickImporter importer = new NewickImporter("(SimSeq1:73.7468,(SimSeq2:25.256989999999995,SimSeq3:45.256989999999995):18.48981);");
        Tree tree = importer.importTree(null);
        TreeModel treeModel = new TreeModel(tree);
        // create Frequency Model
        Parameter freqs = new Parameter.Default(new double[] { 0.25, 0.25, 0.25, 0.25 });
        FrequencyModel freqModel = new FrequencyModel(Nucleotides.INSTANCE, freqs);
        // create branch model
        Parameter kappa1 = new Parameter.Default(1, 1);
        Parameter kappa2 = new Parameter.Default(1, 1);
        HKY hky1 = new HKY(kappa1, freqModel);
        HKY hky2 = new HKY(kappa2, freqModel);
        HomogeneousBranchModel homogenousBranchSubstitutionModel = new HomogeneousBranchModel(hky1);
        List<SubstitutionModel> substitutionModels = new ArrayList<SubstitutionModel>();
        substitutionModels.add(hky1);
        substitutionModels.add(hky2);
        List<FrequencyModel> freqModels = new ArrayList<FrequencyModel>();
        freqModels.add(freqModel);
        Parameter epochTimes = new Parameter.Default(1, 20);
        // create branch rate model
        Parameter rate = new Parameter.Default(1, 0.001);
        BranchRateModel branchRateModel = new StrictClockBranchRates(rate);
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        BranchModel homogeneousBranchModel = new HomogeneousBranchModel(hky1);
        BranchModel epochBranchModel = new EpochBranchModel(treeModel, substitutionModels, epochTimes);
        // create partition
        Partition partition1 = new //
        Partition(//
        treeModel, //
        homogenousBranchSubstitutionModel, //
        siteRateModel, //
        branchRateModel, //
        freqModel, // from
        0, // to
        sequenceLength - 1, // every
        1);
        partitionsList.add(partition1);
        // feed to sequence simulator and generate data
        BeagleSequenceSimulator simulator = new BeagleSequenceSimulator(partitionsList);
        Alignment alignment = simulator.simulate(false, false);
        BeagleTreeLikelihood nbtl = new BeagleTreeLikelihood(alignment, treeModel, homogeneousBranchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, false);
        System.out.println("nBTL(homogeneous) = " + nbtl.getLogLikelihood());
        nbtl = new BeagleTreeLikelihood(alignment, treeModel, epochBranchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, false);
        System.out.println("nBTL(epoch) = " + nbtl.getLogLikelihood());
    } catch (Exception e) {
        e.printStackTrace();
        System.exit(-1);
    }
// END: try-catch block
}

Example 2 with GammaSiteRateModel

use of dr.evomodel.siteratemodel.GammaSiteRateModel in project beast-mcmc by beast-dev.

the class CompleteHistorySimulatorParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    int nReplications = xo.getIntegerAttribute(REPLICATIONS);
    Tree tree = (Tree) xo.getChild(Tree.class);
    GammaSiteRateModel siteModel = (GammaSiteRateModel) xo.getChild(GammaSiteRateModel.class);
    BranchRateModel rateModel = (BranchRateModel) xo.getChild(BranchRateModel.class);
    if (rateModel == null)
        rateModel = new DefaultBranchRateModel();
    DataType dataType = siteModel.getSubstitutionModel().getDataType();
    String jumpTag = xo.getAttribute(JUMP_TAG_NAME, JUMP_TAG);
    boolean sumAcrossSites = xo.getAttribute(SUM_SITES, false);
    Parameter branchSpecificParameter = null;
    Parameter variableValueParameter = null;
    if (xo.hasChildNamed(BRANCH_SPECIFIC_SPECIFICATION)) {
        XMLObject cxo = xo.getChild(BRANCH_SPECIFIC_SPECIFICATION);
        branchSpecificParameter = (Parameter) cxo.getChild(BRANCH_VARIABLE_PARAMETER).getChild(Parameter.class);
        variableValueParameter = (Parameter) cxo.getChild(VARIABLE_VALUE_PARAMETER).getChild(Parameter.class);
    }
    CompleteHistorySimulator history = new CompleteHistorySimulator(tree, siteModel, rateModel, nReplications, sumAcrossSites, branchSpecificParameter, variableValueParameter);
    XMLObject cxo = xo.getChild(COUNTS);
    if (cxo != null) {
        MarkovJumpsTreeLikelihoodParser.parseAllChildren(cxo, history, dataType.getStateCount(), jumpTag, MarkovJumpsType.COUNTS, false);
    }
    cxo = xo.getChild(REWARDS);
    if (cxo != null) {
        MarkovJumpsTreeLikelihoodParser.parseAllChildren(cxo, history, dataType.getStateCount(), jumpTag, MarkovJumpsType.REWARDS, false);
    }
    if (dataType instanceof Codons) {
        Codons codons = (Codons) dataType;
        if (xo.getAttribute(SYN_JUMPS, false)) {
            // use base 61
            double[] synRegMatrix = CodonLabeling.getRegisterMatrix(CodonLabeling.SYN, codons, false);
            Parameter registerParameter = new Parameter.Default(synRegMatrix);
            registerParameter.setId("S");
            history.addRegister(registerParameter, MarkovJumpsType.COUNTS, false);
        }
        if (xo.getAttribute(NON_SYN_JUMPS, false)) {
            // use base 61
            double[] nonSynRegMatrix = CodonLabeling.getRegisterMatrix(CodonLabeling.NON_SYN, codons, false);
            Parameter registerParameter = new Parameter.Default(nonSynRegMatrix);
            registerParameter.setId("N");
            history.addRegister(registerParameter, MarkovJumpsType.COUNTS, false);
        }
    }
    if (xo.getAttribute(ANNOTATE_WITH_ALIGNMENT, false)) {
        history.addAlignmentTrait();
    }
    boolean alignmentOnly = xo.getAttribute(ALIGNMENT_ONLY, false);
    if (dataType instanceof Codons && !alignmentOnly) {
        System.out.println("Codon models give exception when count statistics are done on them. " + "You can supress this by setting alignmentOnly to true.");
    }
    if (alignmentOnly) {
        history.setAlignmentOnly();
    }
    history.simulate();
    return history;
}

Example 3 with GammaSiteRateModel

use of dr.evomodel.siteratemodel.GammaSiteRateModel in project beast-mcmc by beast-dev.

the class AncestralSequenceAnnotator method loadSiteModel.

/*
     * This method is equivalent to the SubstitutionModelLoader without having to
     * be object orientated and can be much more flexible.
     */
private GammaSiteRateModel loadSiteModel(Tree tree) {
    String modelType = (String) tree.getAttribute(SUBST_MODEL);
    /* Identify the datatype and substitution model. Load the model */
    SubstitutionModelLoader sml = null;
    String substModelName = modelType.replaceFirst("\\*.+", "").replaceFirst("\\+.+", "").trim();
    System.out.println("Basic Substitution Model is " + substModelName);
    for (int i = 0; i < GENERAL_MODELS_LIST.length; i++) {
        if (substModelName.matches(GENERAL_MODELS_LIST[i])) {
            sml = new GeneralSubstitutionModelLoader(tree, modelType);
        }
    }
    for (int i = 0; i < NUCLEOTIDE_MODELS_LIST.length; i++) {
        if (substModelName.matches(NUCLEOTIDE_MODELS_LIST[i])) {
            sml = new NucleotideSubstitutionModelLoader(tree, modelType);
        }
    }
    for (int i = 0; i < AMINO_ACID_MODELS_LIST.length; i++) {
        if (substModelName.matches(AMINO_ACID_MODELS_LIST[i])) {
            sml = new AminoAcidSubstitutionModelLoader(tree, modelType);
        }
    }
    for (int i = 0; i < TRIPLET_MODELS_LIST.length; i++) {
        if (substModelName.matches(TRIPLET_MODELS_LIST[i])) {
            sml = new TripletSubstitutionModelLoader(tree, modelType);
        }
    }
    for (int i = 0; i < CODON_MODELS_LIST.length; i++) {
        if (substModelName.matches(CODON_MODELS_LIST[i])) {
            sml = new CodonSubstitutionModelLoader(tree, modelType);
        }
    }
    if (sml.getSubstitutionModel() == null) {
        System.err.println("Substitution model type '" + modelType + "' not implemented");
        System.exit(-1);
    }
    //SiteModel siteModel = new GammaSiteModel(sml.getSubstitutionModel(), new Parameter.Default(1.0), null, 0, null);
    //SiteModel siteModel = new GammaSiteModel(sml.getSubstitutionModel(), null, null, 0, null);
    String siteRatesModels = modelType.substring(modelType.indexOf("+") + 1, modelType.length());
    //String[] siteRatesModels = siteRatesParameters.split(" + ");
    System.out.println("Site rate models: " + siteRatesModels);
    if (sml.getSubstitutionModel().getDataType().getClass().equals(Codons.class) && siteRatesModels.length() > 0) {
        if (siteRatesModels.indexOf("+M2") >= 0) {
            /* M2 */
            System.out.println("Site model - M2 Codon site model used");
            Parameter m2FrequencyAAInv = new Parameter.Default(Double.parseDouble(tree.getAttribute("M2_f[AA INV]").toString()));
            Parameter m2FrequencyNeutral = new Parameter.Default(Double.parseDouble(tree.getAttribute("M2_f[Neutral]").toString()));
            Parameter m2FrequencySelected = new Parameter.Default(Double.parseDouble(tree.getAttribute("M2_f[Selected]").toString()));
            Parameter m2Omega = new Parameter.Default(Double.parseDouble(tree.getAttribute("M2_omega").toString()));
            //Parameter m2FrequencyAAInv = new Parameter.Default((Double) tree.getAttribute("M2_f[AA INV]"));
            //Parameter m2FrequencyNeutral = new Parameter.Default((Double) tree.getAttribute("M2_f[Neutral]"));
            //Parameter m2FrequencySelected = new Parameter.Default((Double) tree.getAttribute("M2_f[Selected]"));
            //Parameter m2Omega = new Parameter.Default((Double) tree.getAttribute("M2_omega"));
            System.err.println("Sorry, M2 substitution model is not yet implemented in BEAST");
            System.exit(0);
        } else if (siteRatesModels.indexOf("+M3") >= 0) {
            /* M3 */
            System.out.println("Site model - M3 Codon site model used");
            int numberOfBins = Integer.parseInt(siteRatesModels.replaceFirst(".+M3\\[", "").replaceFirst("\\].+", ""));
            System.out.println(" + M3 n value: " + numberOfBins);
            Parameter[] m3Frequencies = new Parameter[numberOfBins];
            Parameter[] m3Omegas = new Parameter[numberOfBins];
            for (int i = 1; i <= numberOfBins; i++) {
                m3Frequencies[i - 1] = new Parameter.Default(Double.parseDouble(tree.getAttribute("M3_f" + i).toString()));
                m3Omegas[i - 1] = new Parameter.Default(Double.parseDouble(tree.getAttribute("M3_omega" + i).toString()));
            //m3Frequencies[i-1] = new Parameter.Default((Double) tree.getAttribute("M3_f"+i));
            //m3Omegas[i-1] = new Parameter.Default((Double) tree.getAttribute("M3_omega"+i));
            }
            System.err.println("Sorry, M3 substitution model is not yet implemented in BEAST");
            System.exit(0);
        } else if (siteRatesModels.indexOf("+M0_omega~Beta(") >= 0) {
            /* M7 */
            System.out.println("Site model - M7 Codon site model used");
            int numberOfBins = Integer.parseInt(siteRatesModels.replaceFirst("M0_omega~Beta\\(", "").replaceFirst("\\)", ""));
            System.out.println(" + M7 n value: " + numberOfBins);
            Parameter m7BetaMu = new Parameter.Default(Double.parseDouble(tree.getAttribute("beta_mu").toString()));
            Parameter m7BetaVarMu = new Parameter.Default(Double.parseDouble(tree.getAttribute("beta_Var/mu").toString()));
            //Parameter m7BetaMu = new Parameter.Default((Double) tree.getAttribute("beta_mu"));
            //Parameter m7BetaVarMu = new Parameter.Default((Double) tree.getAttribute("beta_Var/mu"));
            System.err.println("Sorry, M7 substitution model is not yet implemented in BEAST");
            System.exit(0);
        }
    } else if (siteRatesModels.length() > 0) {
        /* i.e. for other data types. */
        /* Do gamma/lognormal + pinv */
        Parameter pInvParameter = null;
        int categories = -1;
        Parameter alphaParameter = null;
        if (siteRatesModels.indexOf("+INV") >= 0) {
            System.out.println("Site model -  proportion of invariable sites used");
            //pInvParameter = new Parameter.Default(((Double) tree.getAttribute("INV_p")).doubleValue());
            pInvParameter = new Parameter.Default(Double.parseDouble(tree.getAttribute("INV_p").toString()));
        }
        if (siteRatesModels.indexOf("+rate~Gamma(") >= 0) {
            System.out.println("Site model - gamma site rate heterogeneity used");
            categories = Integer.parseInt(siteRatesModels.replaceFirst(".+rate~Gamma\\(", "").replaceFirst("\\).*", ""));
            //double sigmaMu = (Double) tree.getAttribute("gamma_sigma/mu");
            double sigmaMu = Double.parseDouble(tree.getAttribute("gamma_sigma/mu").toString());
            sigmaMu = (1.0 / sigmaMu) * (1.0 / sigmaMu);
            /* BAli-Phy is parameterised by sigma/mu instead of alpha */
            alphaParameter = new Parameter.Default(sigmaMu);
        } else if (siteRatesModels.indexOf("+rate~LogNormal(") >= 0) {
            // TODO implement lognormal site model
            System.out.println("Site model - lognormal site rate heterogeneity used");
            System.err.println("Sorry, lognormal site rates are not yet implemented in BEAST");
            System.exit(0);
            categories = Integer.parseInt(siteRatesModels.replaceFirst(".+rate~LogNormal\\(", "").replaceFirst("\\).*", ""));
            //double sigmaMu = (Double) tree.getAttribute("log-normal_sigma/mu");
            double sigmaMu = Double.parseDouble(tree.getAttribute("log-normal_sigma/mu").toString());
            sigmaMu = (1.0 / sigmaMu) * (1.0 / sigmaMu);
            /* BAli-Phy is parameterised by sigma/mu instead of alpha */
            alphaParameter = new Parameter.Default(sigmaMu);
        } else if (siteRatesModels.indexOf("+GAMMA(") >= 0) {
            /* For BEAST output */
            System.out.println("Site model - gamma site rate heterogeneity used");
            categories = Integer.parseInt(siteRatesModels.replaceFirst(".+GAMMA\\(", "").replaceFirst("\\).*", ""));
            //double sigmaMu = (Double) tree.getAttribute("gamma_sigma/mu");
            double alpha = Double.parseDouble(tree.getAttribute("alpha").toString());
            alphaParameter = new Parameter.Default(alpha);
        }
        //System.out.println("alpha and pinv parameters: " + alphaParameter.getParameterValue(0) + "\t" + pInvParameter.getParameterValue(0));
        //GammaSiteRateModel siteModel = new GammaSiteRateModel(sml.getSubstitutionModel(), new Parameter.Default(1.0), alphaParameter, categories, pInvParameter);
        GammaSiteRateModel siteModel = new GammaSiteRateModel(GammaSiteModelParser.SITE_MODEL, new Parameter.Default(1.0), alphaParameter, categories, pInvParameter);
        siteModel.setSubstitutionModel(sml.getSubstitutionModel());
        //SiteModel siteModel = new GammaSiteModel(sml.getSubstitutionModel(), null, null, 0, null);
        return siteModel;
    }
    /* Default with no gamma or pinv */
    //SiteRateModel siteModel = new GammaSiteRateModel(sml.getSubstitutionModel());
    GammaSiteRateModel siteModel = new GammaSiteRateModel(GammaSiteModelParser.SITE_MODEL);
    siteModel.setSubstitutionModel(sml.getSubstitutionModel());
    return siteModel;
}

Example 4 with GammaSiteRateModel

use of dr.evomodel.siteratemodel.GammaSiteRateModel in project beast-mcmc by beast-dev.

the class AncestralSequenceAnnotator method processTree.

private Tree processTree(Tree tree) {
    // Remake tree to fix node ordering - Marc
    GammaSiteRateModel siteModel = loadSiteModel(tree);
    SimpleAlignment alignment = new SimpleAlignment();
    alignment.setDataType(siteModel.getSubstitutionModel().getDataType());
    if (siteModel.getSubstitutionModel().getDataType().getClass().equals(Codons.class)) {
        //System.out.println("trololo");
        alignment.setDataType(Nucleotides.INSTANCE);
    }
    //System.out.println("BOO BOO " + siteModel.getSubstitutionModel().getDataType().getClass().getName()+"\t" + Codons.UNIVERSAL.getClass().getName() + "\t" + alignment.getDataType().getClass().getName());
    // Get sequences
    String[] sequence = new String[tree.getNodeCount()];
    for (int i = 0; i < tree.getNodeCount(); i++) {
        NodeRef node = tree.getNode(i);
        sequence[i] = (String) tree.getNodeAttribute(node, SEQ_STRING);
        if (tree.isExternal(node)) {
            Taxon taxon = tree.getNodeTaxon(node);
            alignment.addSequence(new Sequence(taxon, sequence[i]));
        //System.out.println("seq " + sequence[i]);
        }
    }
    // Make evolutionary model
    BranchRateModel rateModel = new StrictClockBranchRates(new Parameter.Default(1.0));
    FlexibleTree flexTree;
    if (siteModel.getSubstitutionModel().getDataType().getClass().equals(Codons.class)) {
        ConvertAlignment convertAlignment = new ConvertAlignment(siteModel.getSubstitutionModel().getDataType(), ((Codons) siteModel.getSubstitutionModel().getDataType()).getGeneticCode(), alignment);
        flexTree = sampleTree(tree, convertAlignment, siteModel, rateModel);
    //flexTree = sampleTree(tree, alignment, siteModel, rateModel);
    } else {
        flexTree = sampleTree(tree, alignment, siteModel, rateModel);
    }
    introduceGaps(flexTree, tree);
    return flexTree;
}

Example 5 with GammaSiteRateModel

use of dr.evomodel.siteratemodel.GammaSiteRateModel in project beast-mcmc by beast-dev.

the class DataLikelihoodTester2 method main.

public static void main(String[] args) {
    // turn off logging to avoid screen noise...
    Logger logger = Logger.getLogger("dr");
    logger.setUseParentHandlers(false);
    SimpleAlignment alignment = createAlignment(sequences, Nucleotides.INSTANCE);
    TreeModel treeModel;
    try {
        treeModel = createSpecifiedTree("((human:0.1,chimp:0.1):0.1,gorilla:0.2)");
    } catch (Exception e) {
        throw new RuntimeException("Unable to parse Newick tree");
    }
    System.out.print("\nTest BeagleTreeLikelihood (kappa = 1): ");
    //substitutionModel
    Parameter freqs = new Parameter.Default(new double[] { 0.25, 0.25, 0.25, 0.25 });
    Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 1.0, 0, 100);
    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);
    //siteModel
    double alpha = 0.5;
    GammaSiteRateModel siteRateModel = new GammaSiteRateModel("gammaModel", alpha, 4);
    //        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteRateModel");
    siteRateModel.setSubstitutionModel(hky);
    Parameter mu = new Parameter.Default(GammaSiteModelParser.SUBSTITUTION_RATE, 1.0, 0, Double.POSITIVE_INFINITY);
    siteRateModel.setRelativeRateParameter(mu);
    FrequencyModel f2 = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    Parameter kappa2 = new Parameter.Default(HKYParser.KAPPA, 10.0, 0, 100);
    HKY hky2 = new HKY(kappa2, f2);
    GammaSiteRateModel siteRateModel2 = new GammaSiteRateModel("gammaModel", alpha, 4);
    siteRateModel2.setSubstitutionModel(hky2);
    siteRateModel2.setRelativeRateParameter(mu);
    //treeLikelihood
    SitePatterns patterns = new SitePatterns(alignment, null, 0, -1, 1, true);
    BranchModel branchModel = new HomogeneousBranchModel(siteRateModel.getSubstitutionModel(), siteRateModel.getSubstitutionModel().getFrequencyModel());
    BranchModel branchModel2 = new HomogeneousBranchModel(siteRateModel2.getSubstitutionModel(), siteRateModel2.getSubstitutionModel().getFrequencyModel());
    BranchRateModel branchRateModel = new DefaultBranchRateModel();
    BeagleTreeLikelihood treeLikelihood = new BeagleTreeLikelihood(patterns, treeModel, branchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.AUTO, true);
    double logLikelihood = treeLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood);
    System.out.print("\nTest BeagleDataLikelihoodDelegate (kappa = 1): ");
    BeagleDataLikelihoodDelegate dataLikelihoodDelegate = new BeagleDataLikelihoodDelegate(treeModel, patterns, branchModel, siteRateModel, false, PartialsRescalingScheme.NONE, false);
    TreeDataLikelihood treeDataLikelihood = new TreeDataLikelihood(dataLikelihoodDelegate, treeModel, branchRateModel);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood);
    hky.setKappa(5.0);
    System.out.print("\nTest BeagleDataLikelihoodDelegate (kappa = 5): ");
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood);
    System.out.print("\nTest BeagleDataLikelihoodDelegate (kappa = 10): ");
    dataLikelihoodDelegate = new BeagleDataLikelihoodDelegate(treeModel, patterns, branchModel2, siteRateModel2, false, PartialsRescalingScheme.NONE, false);
    treeDataLikelihood = new TreeDataLikelihood(dataLikelihoodDelegate, treeModel, branchRateModel);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood);
    hky2.setKappa(11.0);
    System.out.print("\nTest BeagleDataLikelihoodDelegate (kappa = 11): ");
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood);
    hky.setKappa(1.0);
    hky2.setKappa(10.0);
    MultiPartitionDataLikelihoodDelegate multiPartitionDataLikelihoodDelegate;
    System.out.print("\nTest MultiPartitionDataLikelihoodDelegate 1 partition (kappa = 1):");
    multiPartitionDataLikelihoodDelegate = new MultiPartitionDataLikelihoodDelegate(treeModel, Collections.singletonList((PatternList) patterns), Collections.singletonList((BranchModel) branchModel), Collections.singletonList((SiteRateModel) siteRateModel), true, PartialsRescalingScheme.NONE, false);
    treeDataLikelihood = new TreeDataLikelihood(multiPartitionDataLikelihoodDelegate, treeModel, branchRateModel);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood);
    hky.setKappa(5.0);
    System.out.print("\nTest MultiPartitionDataLikelihoodDelegate 1 partition (kappa = 5):");
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood);
    hky.setKappa(1.0);
    System.out.print("\nTest MultiPartitionDataLikelihoodDelegate 1 partition (kappa = 10):");
    multiPartitionDataLikelihoodDelegate = new MultiPartitionDataLikelihoodDelegate(treeModel, Collections.singletonList((PatternList) patterns), Collections.singletonList((BranchModel) branchModel2), Collections.singletonList((SiteRateModel) siteRateModel2), true, PartialsRescalingScheme.NONE, false);
    treeDataLikelihood = new TreeDataLikelihood(multiPartitionDataLikelihoodDelegate, treeModel, branchRateModel);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood);
    System.out.print("\nTest MultiPartitionDataLikelihoodDelegate 2 partitions (kappa = 1, 10): ");
    List<PatternList> patternLists = new ArrayList<PatternList>();
    patternLists.add(patterns);
    patternLists.add(patterns);
    List<SiteRateModel> siteRateModels = new ArrayList<SiteRateModel>();
    siteRateModels.add(siteRateModel);
    siteRateModels.add(siteRateModel2);
    List<BranchModel> branchModels = new ArrayList<BranchModel>();
    branchModels.add(branchModel);
    branchModels.add(branchModel2);
    multiPartitionDataLikelihoodDelegate = new MultiPartitionDataLikelihoodDelegate(treeModel, patternLists, branchModels, siteRateModels, true, PartialsRescalingScheme.NONE, false);
    treeDataLikelihood = new TreeDataLikelihood(multiPartitionDataLikelihoodDelegate, treeModel, branchRateModel);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood + " (NOT OK: this is 2x the logLikelihood of the 2nd partition)\n\n");
    System.exit(0);
    //START ADDITIONAL TEST #1 - Guy Baele
    System.out.println("-- Test #1 SiteRateModels -- ");
    //alpha in partition 1 reject followed by alpha in partition 2 reject
    System.out.print("Adjust alpha in partition 1: ");
    siteRateModel.setAlpha(0.4);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood);
    System.out.print("Return alpha in partition 1 to original value: ");
    siteRateModel.setAlpha(0.5);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood + " (i.e. reject: OK)\n");
    System.out.print("Adjust alpha in partition 2: ");
    siteRateModel2.setAlpha(0.35);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood);
    System.out.print("Return alpha in partition 2 to original value: ");
    siteRateModel2.setAlpha(0.5);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood + " (i.e. reject: OK)\n");
    //alpha in partition 1 accept followed by alpha in partition 2 accept
    System.out.print("Adjust alpha in partition 1: ");
    siteRateModel.setAlpha(0.4);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood);
    System.out.print("Adjust alpha in partition 2: ");
    siteRateModel2.setAlpha(0.35);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood + " (NOT OK: same logLikelihood as only setting alpha in partition 2)");
    System.out.print("Return alpha in partition 1 to original value: ");
    siteRateModel.setAlpha(0.5);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood + " (NOT OK: alpha in partition 2 has not been returned to original value yet)");
    System.out.print("Return alpha in partition 2 to original value: ");
    siteRateModel2.setAlpha(0.5);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood + "\n");
    //adjusting alphas in both partitions without explicitly calling getLogLikelihood() in between
    System.out.print("Adjust both alphas in partitions 1 and 2: ");
    siteRateModel.setAlpha(0.4);
    siteRateModel2.setAlpha(0.35);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood);
    System.out.print("Return alpha in partition 2 to original value: ");
    siteRateModel2.setAlpha(0.5);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood + " (NOT OK: alpha in partition 1 has not been returned to original value yet)");
    System.out.print("Return alpha in partition 1 to original value: ");
    siteRateModel.setAlpha(0.5);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood + "\n\n");
    //END ADDITIONAL TEST - Guy Baele
    //START ADDITIONAL TEST #2 - Guy Baele
    System.out.println("-- Test #2 SiteRateModels -- ");
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood);
    //1 siteRateModel shared across 2 partitions
    siteRateModels = new ArrayList<SiteRateModel>();
    siteRateModels.add(siteRateModel);
    multiPartitionDataLikelihoodDelegate = new MultiPartitionDataLikelihoodDelegate(treeModel, patternLists, branchModels, siteRateModels, true, PartialsRescalingScheme.NONE, false);
    treeDataLikelihood = new TreeDataLikelihood(multiPartitionDataLikelihoodDelegate, treeModel, branchRateModel);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood + "\n");
    System.out.print("Adjust alpha in shared siteRateModel: ");
    siteRateModel.setAlpha(0.4);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood + " (NOT OK: same logLikelihood as only adjusted alpha for partition 1)");
    siteRateModel.setAlpha(0.5);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood + "\n\n");
    //END ADDITIONAL TEST - Guy Baele
    //START ADDITIONAL TEST #3 - Guy Baele
    System.out.println("-- Test #3 SiteRateModels -- ");
    siteRateModel = new GammaSiteRateModel("gammaModel");
    siteRateModel.setSubstitutionModel(hky);
    siteRateModel.setRelativeRateParameter(mu);
    siteRateModel2 = new GammaSiteRateModel("gammaModel2");
    siteRateModel2.setSubstitutionModel(hky2);
    siteRateModel2.setRelativeRateParameter(mu);
    siteRateModels = new ArrayList<SiteRateModel>();
    siteRateModels.add(siteRateModel);
    siteRateModels.add(siteRateModel2);
    multiPartitionDataLikelihoodDelegate = new MultiPartitionDataLikelihoodDelegate(treeModel, patternLists, branchModels, siteRateModels, true, PartialsRescalingScheme.NONE, false);
    treeDataLikelihood = new TreeDataLikelihood(multiPartitionDataLikelihoodDelegate, treeModel, branchRateModel);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood + "\n");
    System.out.print("Adjust kappa in partition 1: ");
    hky.setKappa(5.0);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood + " (NOT OK: logLikelihood has not changed?)");
    System.out.print("Return kappa in partition 1 to original value: ");
    hky.setKappa(1.0);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood + "\n");
    System.out.print("Adjust kappa in partition 2: ");
    hky2.setKappa(11.0);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood);
    System.out.print("Return kappa in partition 2 to original value: ");
    hky2.setKappa(10.0);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.println("logLikelihood = " + logLikelihood + " (i.e. reject: OK)\n\n");
    //END ADDITIONAL TEST - Guy Baele
    //START ADDITIONAL TEST #4 - Guy Baele
    System.out.println("-- Test #4 SiteRateModels -- ");
    SimpleAlignment secondAlignment = createAlignment(moreSequences, Nucleotides.INSTANCE);
    SitePatterns morePatterns = new SitePatterns(secondAlignment, null, 0, -1, 1, true);
    BeagleDataLikelihoodDelegate dataLikelihoodDelegateOne = new BeagleDataLikelihoodDelegate(treeModel, patterns, branchModel, siteRateModel, false, PartialsRescalingScheme.NONE, false);
    TreeDataLikelihood treeDataLikelihoodOne = new TreeDataLikelihood(dataLikelihoodDelegateOne, treeModel, branchRateModel);
    logLikelihood = treeDataLikelihoodOne.getLogLikelihood();
    System.out.println("\nBeagleDataLikelihoodDelegate logLikelihood partition 1 (kappa = 1) = " + logLikelihood);
    hky.setKappa(10.0);
    logLikelihood = treeDataLikelihoodOne.getLogLikelihood();
    System.out.println("BeagleDataLikelihoodDelegate logLikelihood partition 1 (kappa = 10) = " + logLikelihood);
    hky.setKappa(1.0);
    BeagleDataLikelihoodDelegate dataLikelihoodDelegateTwo = new BeagleDataLikelihoodDelegate(treeModel, morePatterns, branchModel2, siteRateModel2, false, PartialsRescalingScheme.NONE, false);
    TreeDataLikelihood treeDataLikelihoodTwo = new TreeDataLikelihood(dataLikelihoodDelegateTwo, treeModel, branchRateModel);
    logLikelihood = treeDataLikelihoodTwo.getLogLikelihood();
    System.out.println("BeagleDataLikelihoodDelegate logLikelihood partition 2 (kappa = 10) = " + logLikelihood + "\n");
    multiPartitionDataLikelihoodDelegate = new MultiPartitionDataLikelihoodDelegate(treeModel, Collections.singletonList((PatternList) patterns), Collections.singletonList((BranchModel) branchModel), Collections.singletonList((SiteRateModel) siteRateModel), true, PartialsRescalingScheme.NONE, false);
    treeDataLikelihood = new TreeDataLikelihood(multiPartitionDataLikelihoodDelegate, treeModel, branchRateModel);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.print("Test MultiPartitionDataLikelihoodDelegate 1st partition (kappa = 1):");
    System.out.println("logLikelihood = " + logLikelihood);
    hky.setKappa(10.0);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.print("Test MultiPartitionDataLikelihoodDelegate 1st partition (kappa = 10):");
    System.out.println("logLikelihood = " + logLikelihood);
    hky.setKappa(1.0);
    multiPartitionDataLikelihoodDelegate = new MultiPartitionDataLikelihoodDelegate(treeModel, Collections.singletonList((PatternList) morePatterns), Collections.singletonList((BranchModel) branchModel2), Collections.singletonList((SiteRateModel) siteRateModel2), true, PartialsRescalingScheme.NONE, false);
    treeDataLikelihood = new TreeDataLikelihood(multiPartitionDataLikelihoodDelegate, treeModel, branchRateModel);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.print("Test MultiPartitionDataLikelihoodDelegate 2nd partition (kappa = 10):");
    System.out.println("logLikelihood = " + logLikelihood + "\n");
    patternLists = new ArrayList<PatternList>();
    patternLists.add(patterns);
    patternLists.add(morePatterns);
    multiPartitionDataLikelihoodDelegate = new MultiPartitionDataLikelihoodDelegate(treeModel, patternLists, branchModels, siteRateModels, true, PartialsRescalingScheme.NONE, false);
    treeDataLikelihood = new TreeDataLikelihood(multiPartitionDataLikelihoodDelegate, treeModel, branchRateModel);
    logLikelihood = treeDataLikelihood.getLogLikelihood();
    System.out.print("Test MultiPartitionDataLikelihoodDelegate 2 partitions (kappa = 1, 10): ");
    System.out.println("logLikelihood = " + logLikelihood + " (NOT OK: should be the sum of both separate logLikelihoods)\nKappa value of partition 2 is used to compute logLikelihood for both partitions?");
//END ADDITIONAL TEST - Guy Baele
}