add grad wrt mean for GP and clean some warnings

beast-dev · Nov 30, 2023 · c66bbea · c66bbea
1 parent 3e13a4e
commit c66bbea
Show file tree

Hide file tree

Showing 4 changed files with 37 additions and 9 deletions.
diff --git a/src/dr/evomodelxml/substmodel/GlmSubstitutionModelParser.java b/src/dr/evomodelxml/substmodel/GlmSubstitutionModelParser.java
@@ -99,7 +99,7 @@ public XMLSyntaxRule[] getSyntaxRules() {
         return rules;
     }
 
-    private XMLSyntaxRule[] rules = new XMLSyntaxRule[]{
+    private final XMLSyntaxRule[] rules = new XMLSyntaxRule[]{
             new XORRule(
                     new StringAttributeRule(DataType.DATA_TYPE, "The type of sequence data",
                             DataType.getRegisteredDataTypeNames(), false),

diff --git a/src/dr/inferencexml/distribution/GaussianProcessFieldParser.java b/src/dr/inferencexml/distribution/GaussianProcessFieldParser.java
@@ -37,7 +37,6 @@
 
 import static dr.math.distributions.gp.AdditiveGaussianProcessDistribution.BasisDimension;
 import static dr.inferencexml.distribution.RandomFieldParser.WEIGHTS_RULE;
-import static dr.inferencexml.distribution.RandomFieldParser.parseWeightProvider;
 
 public class GaussianProcessFieldParser extends AbstractXMLObjectParser {
 

diff --git a/src/dr/math/distributions/GaussianMarkovRandomField.java b/src/dr/math/distributions/GaussianMarkovRandomField.java
@@ -496,8 +496,7 @@ static class SymmetricTriDiagonalMatrix {
         double[] offDiagonal;
 
         SymmetricTriDiagonalMatrix(int dim) {
-            this.diagonal = new double[dim];
-            this.offDiagonal = new double[dim - 1];
+            this(new double[dim], new double[dim - 1]);
         }
 
         SymmetricTriDiagonalMatrix(double[] diagonal, double[] offDiagonal) {

diff --git a/src/dr/math/distributions/gp/AdditiveGaussianProcessDistribution.java b/src/dr/math/distributions/gp/AdditiveGaussianProcessDistribution.java
@@ -226,11 +226,6 @@ public double[] getMean() {
         return mean;
     }
 
-    @Override
-    public GradientProvider getGradientWrt(Parameter parameter) {
-        throw new RuntimeException("Unknown parameter");
-    }
-
     @Override
     public String getType() {
         return TYPE;
@@ -313,6 +308,41 @@ protected void restoreState() {
     @Override
     protected void acceptState() { }
 
+    @Override
+    public GradientProvider getGradientWrt(Parameter parameter) {
+        if (parameter == meanParameter) {
+            return new GradientProvider() {
+                @Override
+                public int getDimension() {
+                    return meanParameter.getDimension();
+                }
+
+                @Override
+                public double[] getGradientLogDensity(Object x) {
+
+                    double[] gradient = gradLogPdf((double[]) x, getMean(), getPrecision());
+
+                    if (meanParameter.getDimension() == dim) {
+                        for (int i = 0; i < dim; ++i) {
+                            gradient[i] *= -1;
+                        }
+                        return gradient;
+                    } else if (meanParameter.getDimension() == 1) {
+                        double sum = 0.0;
+                        for (int i = 0; i < dim; ++i) {
+                            sum += gradient[i];
+                        }
+                        return new double[]{sum}; // TODO should this be -sum?
+                    }
+
+                    throw new IllegalArgumentException("Unknown mean parameter structure");
+                }
+            };
+        } else {
+            throw new RuntimeException("Unknown parameter");
+        }
+    }
+
     public static class BasisDimension {
 
         private final GaussianProcessKernel kernel;