some code duplication to fix dQ bug when GLM subsitution has X entry …

…negative which is turned to zero by SetupQMatrix in ComplexSubstitutionModel

src/dr/evomodel/substmodel/DifferentiableSubstitutionModelUtil.java

@@ -138,7 +138,7 @@ public static WrappedMatrix getInfinitesimalDifferentialMatrix(WrtParameter wrt,
         ((DifferentiableSubstitutionModel) substitutionModel).setupDifferentialFrequency(wrt, differentialFrequencies);
 
         double[][] differentialMassMatrix = new double[stateCount][stateCount];
-        substitutionModel.setupQMatrix(differentialRates, differentialFrequencies, differentialMassMatrix);
+        setupQDerivative(substitutionModel, differentialRates, differentialFrequencies, differentialMassMatrix);
         substitutionModel.makeValid(differentialMassMatrix, stateCount);
 
         final double weightedNormalizationGradient
@@ -160,6 +160,28 @@ public static WrappedMatrix getInfinitesimalDifferentialMatrix(WrtParameter wrt,
         return differential;
     }
 
+    private static void setupQDerivative(BaseSubstitutionModel substitutionModel, double[] differentialRates,
+                                         double[] differentialFrequencies, double[][] differentialMassMatrix) {
+        if (substitutionModel instanceof ComplexSubstitutionModel) {
+            int i, j, k = 0;
+            final int stateCount = differentialFrequencies.length;
+            for (i = 0; i < stateCount; i++) {
+                for (j = i + 1; j < stateCount; j++) {
+                    final double thisRate = differentialRates[k++];
+                    differentialMassMatrix[i][j] = thisRate * differentialFrequencies[j];
+                }
+            }
+            for (j = 0; j < stateCount; j++) {
+                for (i = j + 1; i < stateCount; i++) {
+                    final double thisRate = differentialRates[k++];
+                    differentialMassMatrix[i][j] = thisRate * differentialFrequencies[j];
+                }
+            }
+        } else {
+            substitutionModel.setupQMatrix(differentialRates, differentialFrequencies, differentialMassMatrix);
+        }
+    }
+
     private static final boolean CHECK_COMMUTABILITY = false;
     private static final double COMMUTABILITY_CHECK_THRESHOLD = 0.01;