Merge branch 'hmc-clock' of github.com:beast-dev/beast-mcmc into hmc-…

…clock

src/dr/app/beast/development_parsers.properties

@@ -95,6 +95,7 @@ dr.evomodelxml.coalescent.operators.GaussianProcessSkytrackBlockUpdateOperatorPa
 dr.evomodelxml.coalescent.operators.GaussianProcessSkytrackTreeOperatorParser
 dr.inferencexml.distribution.GeneralizedAdditiveGaussianProcessModelParser
 dr.inferencexml.distribution.GaussianProcessBasisApproximationParser
+dr.evomodelxml.coalescent.MultilocusNonParametricCoalescentLikelihoodParser
 
 
 # TREE SUMMARY STATISTICS

src/dr/evomodel/coalescent/MultilocusNonparametricCoalescentLikelihood.java

@@ -123,7 +123,7 @@ protected void handleModelChangedEvent(Model model, Object object, int index) {
 
     @Override
     protected void handleVariableChangedEvent(Variable variable, int index, Parameter.ChangeType type) {
-
+        likelihoodKnown = false;
     }
 
     private int moveToNextTimeIndex(int treeIndex, int lastTimeIndex, double[] times) {

src/dr/evomodel/substmodel/InfinitesimalRatesLogger.java

@@ -30,13 +30,18 @@
 import dr.inference.loggers.LogColumn;
 import dr.inference.loggers.Loggable;
 import dr.inference.loggers.NumberColumn;
+import dr.util.Transform;
 
 public class InfinitesimalRatesLogger implements Loggable {
 
-    public InfinitesimalRatesLogger(SubstitutionModel substitutionModel) {
+    public InfinitesimalRatesLogger(SubstitutionModel substitutionModel, boolean diagonalElements, Transform transform) {
         this.substitutionModel = substitutionModel;
+        this.diagonalElements = diagonalElements;
+        this.transform = transform;
     }
 
+
+
     @Override
     public LogColumn[] getColumns() {
         int stateCount = substitutionModel.getDataType().getStateCount();
@@ -45,26 +50,54 @@ public LogColumn[] getColumns() {
             generator = new double[stateCount * stateCount];
         }
 
-        LogColumn[] columns = new LogColumn[stateCount * stateCount];
+        int nOutputs = stateCount * stateCount;
+        if (!diagonalElements) nOutputs -= stateCount;
+        LogColumn[] columns = new LogColumn[nOutputs];
 
+        int index = 0;
         for (int i = 0; i < stateCount; ++i) {
             for (int j = 0; j < stateCount; ++j) {
-                final int k = i * stateCount + j;
+                if (!diagonalElements && i == j) {
+                    continue;
+                }
+                final int row = i;
+                final int col = j;
+                final int k = index++;  // Use index to store column number
                 columns[k] = new NumberColumn(substitutionModel.getId() + "." + (i + 1) + "." + (j + 1)) {
                     @Override
                     public double getDoubleValue() {
                         if (k == 0) { // Refresh at first-element read
                             substitutionModel.getInfinitesimalMatrix(generator);
                         }
-                        return generator[k];
+                        if (transform != null) {
+                            return transform.transform(generator[row * stateCount + col]);
+                        } else {
+                            return generator[row * stateCount + col];
+                        }
                     }
                 };
             }
         }
 
+//        for (int i = 0; i < stateCount; ++i) {
+//            for (int j = 0; j < stateCount; ++j) {
+//                final int k = i * stateCount + j;
+//                columns[k] = new NumberColumn(substitutionModel.getId() + "." + (i + 1) + "." + (j + 1)) {
+//                    @Override
+//                    public double getDoubleValue() {
+//                        if (k == 0) { // Refresh at first-element read
+//                            substitutionModel.getInfinitesimalMatrix(generator);
+//                        }
+//                        return generator[k];
+//                    }
+//                };
+//            }
+//        }
+
         return columns;
     }
-
+    private final Transform transform;
     private final SubstitutionModel substitutionModel;
+    private final Boolean diagonalElements;
     private double[] generator;
 }

src/dr/evomodelxml/coalescent/MultilocusNonParametricCoalescentLikelihoodParser.java

@@ -0,0 +1,154 @@
+/*
+ * GMRFSkyrideLikelihoodParser.java
+ *
+ * Copyright (c) 2002-2024 Alexei Drummond, Andrew Rambaut and Marc Suchard
+ *
+ * This file is part of BEAST.
+ * 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 dr.evomodelxml.coalescent;
+
+import dr.evolution.coalescent.TreeIntervals;
+import dr.evomodel.coalescent.MultilocusNonparametricCoalescentLikelihood;
+import dr.evomodel.coalescent.smooth.SkyGlideLikelihood;
+import dr.evomodel.tree.TreeModel;
+import dr.evomodelxml.coalescent.smooth.SmoothSkygridLikelihoodParser;
+import dr.inference.model.Parameter;
+import dr.xml.*;
+
+import java.util.ArrayList;
+import java.util.List;
+
+public class MultilocusNonParametricCoalescentLikelihoodParser extends AbstractXMLObjectParser {
+
+    private static final String PARSER_NAME = "multiLocusNPCoalescentLikelihood";
+    private static final String POPULATION_TREE = GMRFSkyrideLikelihoodParser.POPULATION_TREE;
+    private static final String POPULATION_PARAMETER = GMRFSkyrideLikelihoodParser.POPULATION_PARAMETER;
+    //    public static final String PRECISION_PARAMETER = "precisionParameter";
+
+    public static final String PLOIDY = "ploidy";
+
+    private static final String GRID_POINTS = GMRFSkyrideLikelihoodParser.GRID_POINTS;
+    private static final String NUM_GRID_POINTS = GMRFSkyrideLikelihoodParser.NUM_GRID_POINTS;
+    private static final String CUT_OFF = GMRFSkyrideLikelihoodParser.CUT_OFF;
+
+    @Override
+    public Object parseXMLObject(XMLObject xo) throws XMLParseException {
+
+        Parameter logPopSizes = (Parameter) xo.getElementFirstChild(POPULATION_PARAMETER);
+        Parameter gridPoints = SmoothSkygridLikelihoodParser.getGridPoints(xo); //TODO should I call this parser?
+        int nGridPoints = gridPoints.getDimension();
+
+        //        cxo = xo.getChild(PRECISION_PARAMETER);
+        //        Parameter precParameter = (Parameter) cxo.getChild(Parameter.class);
+
+        List<TreeModel> trees = new ArrayList<>();
+        XMLObject cxo = xo.getChild(POPULATION_TREE);
+        for (int i = 0; i < cxo.getChildCount(); i++) {
+            trees.add((TreeModel) cxo.getChild(i));
+        }
+
+        List<TreeIntervals> intervalLists = new ArrayList<>();
+        for (int i = 0; i < trees.size(); i++) {
+            TreeIntervals treeIntervals = new TreeIntervals(trees.get(i));
+            intervalLists.add(treeIntervals);
+        }
+        //        List<BigFastTreeIntervals> intervalLists = new ArrayList<>();
+        //        for (int i = 0; i < trees.size(); i++) {
+        //            BigFastTreeIntervals treeIntervals = new BigFastTreeIntervals(trees.get(i));
+        //            intervalLists.add((BigFastTreeIntervals) treeIntervals);
+        //        }
+
+        Parameter ploidyFactors = parsePloidyFactors(xo, trees, nGridPoints);
+
+        MultilocusNonparametricCoalescentLikelihood likelihood = new MultilocusNonparametricCoalescentLikelihood(intervalLists, logPopSizes, gridPoints, ploidyFactors);
+        return likelihood;
+    }
+
+    protected Parameter parsePloidyFactors(XMLObject xo, List<TreeModel> trees, int nGridPoints) {
+        Parameter ploidyFactors = null;
+        if (xo.getChild(PLOIDY) != null) {
+            XMLObject cxo = xo.getChild(PLOIDY);
+            ploidyFactors = (Parameter) cxo.getChild(Parameter.class);
+        } else {
+            ploidyFactors = new Parameter.Default(PLOIDY, trees.size());
+            for (int i = 0; i < trees.size(); i++) {
+                ploidyFactors.setParameterValue(i, 1.0);
+            }
+//            if (nGridPoints != 0) {
+//                ploidyFactors = new Parameter.Default(PLOIDY, nGridPoints);
+//                for(int i = 0; i < nGridPoints; i++){
+//                    ploidyFactors.setParameterValue(i, 1.0);
+//                }
+//            } else {
+//                ploidyFactors = new Parameter.Default(PLOIDY, trees.size());
+//                for (int i = 0; i < trees.size(); i++) {
+//                    ploidyFactors.setParameterValue(i, 1.0);
+//                }
+//            }
+        }
+        return ploidyFactors;
+    }
+
+    @Override
+    public XMLSyntaxRule[] getSyntaxRules() {
+        return new XMLSyntaxRule[0];
+    }
+
+    private final XMLSyntaxRule[] rules = {
+            new ElementRule(POPULATION_TREE, new XMLSyntaxRule[]{
+                    new ElementRule(TreeModel.class, 1, Integer.MAX_VALUE)
+            }),
+
+            new ElementRule(POPULATION_PARAMETER, new XMLSyntaxRule[]{
+                    new ElementRule(Parameter.class)
+            }),
+
+            new XORRule(
+                    new ElementRule(GRID_POINTS, new XMLSyntaxRule[]{
+                            new ElementRule(Parameter.class)
+                    }),
+                    new AndRule(
+                            new ElementRule(CUT_OFF, new XMLSyntaxRule[]{
+                                    new ElementRule(Parameter.class)
+                            }),
+                            new ElementRule(NUM_GRID_POINTS, new XMLSyntaxRule[]{
+                                    new ElementRule(Parameter.class)
+                            })
+                    )
+            ),
+
+    };
+
+    @Override
+    public String getParserDescription() {
+        return null;
+    }
+
+    @Override
+    public Class getReturnType() {
+        return MultilocusNonparametricCoalescentLikelihood.class;
+    }
+
+    @Override
+    public String getParserName() {
+        return PARSER_NAME;
+    }
+}

src/dr/evomodelxml/substmodel/InfinitesimalRatesLoggerParser.java

@@ -29,25 +29,43 @@
 
 import dr.evomodel.substmodel.InfinitesimalRatesLogger;
 import dr.evomodel.substmodel.SubstitutionModel;
+import dr.util.Transform;
 import dr.xml.*;
 
 public class InfinitesimalRatesLoggerParser extends AbstractXMLObjectParser {
 
     private static final String NAME = "infinitesimalRatesLogger";
+    private static final String DIAGONAL_ELEMENTS = "diagonalElements";
 
     @Override
     public Object parseXMLObject(XMLObject xo) throws XMLParseException {
         SubstitutionModel substitutionModel = (SubstitutionModel) xo.getChild(SubstitutionModel.class);
-        return new InfinitesimalRatesLogger(substitutionModel);
+
+        boolean diagonalElements = true;
+        if (xo.hasAttribute(DIAGONAL_ELEMENTS)) {
+            diagonalElements = xo.getAttribute(DIAGONAL_ELEMENTS, true);
+        }
+
+        Transform.ParsedTransform pt = (Transform.ParsedTransform) xo.getChild(Transform.ParsedTransform.class);
+        Transform transform = null;
+        if (pt != null) {
+            transform = pt.transform;
+        }
+
+        return new InfinitesimalRatesLogger(substitutionModel, diagonalElements, transform);
     }
 
     @Override
     public XMLSyntaxRule[] getSyntaxRules() {
-        return new XMLSyntaxRule[] {
-                new ElementRule(SubstitutionModel.class),
-        };
+        return rules;
     }
 
+    private final XMLSyntaxRule[] rules = new XMLSyntaxRule[]{
+            AttributeRule.newBooleanRule(DIAGONAL_ELEMENTS, true),
+            new ElementRule(SubstitutionModel.class),
+            new ElementRule(Transform.ParsedTransform.class, true)
+    };
+
     @Override
     public String getParserDescription() {
         return "Logger to report infinitesimal rates of a substitution model";

src/dr/math/distributions/gp/AdditiveGaussianProcessDistribution.java

@@ -304,13 +304,6 @@ protected void handleModelChangedEvent(Model model, Object object, int index) {
           throw new IllegalArgumentException("Unknown model");
         }
     }
-//
-//        if (model == bases){
-//            precisionAndDeterminantKnown = false;
-//        } else {
-//            throw new IllegalArgumentException("Unknown model");
-//        }
-
 
     private boolean containsKernel(Model model) {
         for (BasisDimension basis : bases) {