matching "<= " and "="

beast-dev · Jul 26, 2023 · a329583 · a329583
1 parent 1702300
commit a329583
Show file tree

Hide file tree

Showing 2 changed files with 27 additions and 12 deletions.
diff --git a/src/dr/evomodel/speciation/MasBirthDeathSerialSamplingModel.java b/src/dr/evomodel/speciation/MasBirthDeathSerialSamplingModel.java
@@ -59,7 +59,7 @@ final void accumulateGradientForInterval(final double[] gradient, final int curr
         }
     }
 
-    final void accumulateGradientForSampling(double[] gradient, int currentModelSegment, double term1,
+    final void accumulateGradientForSerialSampling(double[] gradient, int currentModelSegment, double term1,
                                              double[] intermediate) {
 
         for (int k = 0; k <= currentModelSegment; k++) {
@@ -70,6 +70,17 @@ final void accumulateGradientForSampling(double[] gradient, int currentModelSegm
 
     }
 
+    final void accumulateGradientForIntensiveSampling(double[] gradient, int currentModelSegment, double term1,
+                                                   double[] intermediate) {
+
+        for (int k = 0; k < currentModelSegment; k++) {
+            for (int p = 0; p < 4; ++p) {
+                gradient[k * 5 + p] += term1 * intermediate[k * 4 + p];
+            }
+        }
+
+    }
+
     final void dPCompute(int model, double t, double intervalStart, double eAt, double[] dP, double[] dG2) {
 
         double G1 = g1(eAt);

diff --git a/src/dr/evomodel/speciation/NewBirthDeathSerialSamplingModel.java b/src/dr/evomodel/speciation/NewBirthDeathSerialSamplingModel.java
@@ -868,14 +868,13 @@ void accumulateGradientForInterval(double[] gradient, int currentModelSegment, i
         for (int p = 0; p < 4; ++p) {
             if (gradientFlags[p]) {
                 for (int k = 0; k <= currentModelSegment; k++) {
-                    gradient[k * 5 + p] += nLineages * (partialQ_all_old[k * 4 + p] / Q_Old
-                            - partialQ_all_young[k * 4 + p] / Q_young);
+                    gradient[k * 5 + p] += nLineages * (partialQ_all_old[k * 4 + p] / Q_Old - partialQ_all_young[k * 4 + p] / Q_young);
                 }
             }
         }
     }
 
-    void accumulateGradientForSampling(double[] gradient, int currentModelSegment, double term1,
+    void accumulateGradientForSerialSampling(double[] gradient, int currentModelSegment, double term1,
                                        double[] intermediate) {
         for (int p = 0; p < 4; p++) {
             if (gradientFlags[p]) {
@@ -886,6 +885,17 @@ void accumulateGradientForSampling(double[] gradient, int currentModelSegment, d
         }
     }
 
+    void accumulateGradientForIntensiveSampling(double[] gradient, int currentModelSegment, double term1,
+                                             double[] intermediate) {
+        for (int p = 0; p < 4; p++) {
+            if (gradientFlags[p]) {
+                for (int k = 0; k < currentModelSegment; k++) {
+                    gradient[k * 5 + p] += term1 * intermediate[k * 4 + p];
+                }
+            }
+        }
+    }
+
     @Override
     public void processGradientSampling(double[] gradient, int currentModelSegment, double intervalEnd) {
 
@@ -935,7 +945,7 @@ public void processGradientSampling(double[] gradient, int currentModelSegment,
 //                }
 //            }
 
-            accumulateGradientForSampling(gradient, currentModelSegment, term1, dPIntervalEnd);
+            accumulateGradientForSerialSampling(gradient, currentModelSegment, term1, dPIntervalEnd);
         }
 
         if (sampleIsAtEventTime && currentModelSegment > 0) {
@@ -944,12 +954,6 @@ public void processGradientSampling(double[] gradient, int currentModelSegment,
 
             double term1 = (1 - r) / ((1 - r) * previousP + r);
 
-//            for (int k = 0; k < currentModelSegment; k++) {
-//                gradient[birthIndex(k, numIntervals)] += term1 * dPModelEnd_prev[k * 4 + 0];
-//                gradient[deathIndex(k, numIntervals)] += term1 * dPModelEnd_prev[k * 4 + 1];
-//                gradient[samplingIndex(k, numIntervals)] += term1 * dPModelEnd_prev[k * 4 + 2];
-//                gradient[fractionIndex(k, numIntervals)] += term1 * dPModelEnd_prev[k * 4 + 3];
-//            }
 
 //            for (int p = 0; p < 4; p++) {
 //                if (gradientFlags[p]) {
@@ -959,7 +963,7 @@ public void processGradientSampling(double[] gradient, int currentModelSegment,
 //                }
 //            }
 
-            accumulateGradientForSampling(gradient, currentModelSegment, term1, dPModelEnd_prev);
+            accumulateGradientForIntensiveSampling(gradient, currentModelSegment, term1, dPModelEnd_prev);
         }
     }