#558 WWZ: replace Gauss–Jordan matinv with closed-form 3×3 symmetric inverse

- Faster, no per-call double[][] allocation; add useLegacyMatinv for Gauss–Jordan benchmark
- Add WWZMatinvBenchmarkTest to compare implementations

Made-with: Cursor

Author: dbenn

src/org/aavso/tools/vstar/util/period/wwz/WeightedWaveletZTransform.java

@@ -55,6 +55,9 @@
  */
 public class WeightedWaveletZTransform implements IAlgorithm {
 
+	/** For benchmark only: use legacy Gauss-Jordan matinv instead of closed-form. Package visibility for tests. */
+	static boolean useLegacyMatinv = false;
+
 	// Observations to be analysed.
 	private List<ValidObservation> obs;
 
@@ -399,25 +402,63 @@ public void make_freqs_from_period_range(double minPer, double maxPer,
 	}
 
 	/**
-	 * Invert the matrix of the wwz equations...
+	 * Invert the 3×3 symmetric matrix in dmat (in-place).
+	 * Uses a closed-form cofactor formula. Faster than Gauss-Jordan for fixed 3×3,
+	 * no pivoting logic, and avoids the per-(tau,freq) loop overhead.
 	 */
 	private void matinv() throws InterruptedException {
-		double dsol[][] = new double[3][3];// (0:2,0:2);
-		double dfac;
+		if (useLegacyMatinv) {
+			matinvGaussJordan();
+			return;
+		}
+		// Symmetric 3×3: a=dmat[0][0], b=dmat[0][1], c=dmat[0][2], d=dmat[1][1], e=dmat[1][2], f=dmat[2][2]
+		double a = dmat[0][0], b = dmat[0][1], c = dmat[0][2];
+		double d = dmat[1][1], e = dmat[1][2], f = dmat[2][2];
+
+		// det = a(df - e²) - b(bf - ce) + c(be - cd)
+		double det = a * (d * f - e * e) - b * (b * f - c * e) + c * (b * e - c * d);
+		if (Math.abs(det) <= 1.0e-14) {
+			return; // singular, leave dmat unchanged (matches previous behaviour)
+		}
+		double invDet = 1.0 / det;
 
-		int ndim = 2;
+		if (interrupted) {
+			throw new InterruptedException();
+		}
 
+		// Adjugate (symmetric for symmetric input), then inv = adj/det
+		double a00 = (d * f - e * e) * invDet;
+		double a01 = (c * e - b * f) * invDet;
+		double a02 = (b * e - c * d) * invDet;
+		double a11 = (a * f - c * c) * invDet;
+		double a12 = (c * b - a * e) * invDet;
+		double a22 = (a * d - b * b) * invDet;
+
+		dmat[0][0] = a00;
+		dmat[0][1] = a01;
+		dmat[0][2] = a02;
+		dmat[1][0] = a01;
+		dmat[1][1] = a11;
+		dmat[1][2] = a12;
+		dmat[2][0] = a02;
+		dmat[2][1] = a12;
+		dmat[2][2] = a22;
+	}
+
+	/** Legacy Gauss-Jordan 3×3 in-place inverse; used only when useLegacyMatinv is true (benchmark). */
+	private void matinvGaussJordan() throws InterruptedException {
+		double dsol[][] = new double[3][3];
+		double dfac;
+		int ndim = 2;
 		for (int i = 0; i <= 2; i++) {
 			for (int j = 0; j <= 2; j++) {
 				dsol[i][j] = 0.0;
 			}
 			dsol[i][i] = 1.0;
-			
 			if (interrupted) {
 				throw new InterruptedException();
 			}
 		}
-
 		for (int i = 0; i <= ndim; i++) {
 			if (dmat[i][i] == 0.0) {
 				if (i == ndim)
@@ -426,16 +467,14 @@ private void matinv() throws InterruptedException {
 					if (dmat[j][i] != 0.0) {
 						for (int k = 0; k <= ndim; k++) {
 							dmat[i][k] = dmat[i][k] + dmat[j][k];
-							dsol[i][j] = dsol[i][j] + dsol[j][k];
+							dsol[i][k] = dsol[i][k] + dsol[j][k];
 						}
 					}
 				}
-				
 				if (interrupted) {
 					throw new InterruptedException();
 				}
 			}
-			
 			dfac = dmat[i][i];
 			for (int j = 0; j <= ndim; j++) {
 				dmat[i][j] = dmat[i][j] / dfac;
@@ -448,7 +487,6 @@ private void matinv() throws InterruptedException {
 						dmat[j][k] = dmat[j][k] - (dmat[i][k] * dfac);
 						dsol[j][k] = dsol[j][k] - (dsol[i][k] * dfac);
 					}
-					
 					if (interrupted) {
 						throw new InterruptedException();
 					}
@@ -459,7 +497,6 @@ private void matinv() throws InterruptedException {
 			for (int j = 0; j <= ndim; j++) {
 				dmat[i][j] = dsol[i][j];
 			}
-			
 			if (interrupted) {
 				throw new InterruptedException();
 			}

test/org/aavso/tools/vstar/util/period/wwz/WWZMatinvBenchmarkTest.java

@@ -0,0 +1,86 @@
+/**
+ * VStar: a statistical analysis tool for variable star data.
+ * Copyright (C) 2010  AAVSO (http://www.aavso.org/)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Affero General Public License as
+ * published by the Free Software Foundation, either version 3 of the
+ * License, or (at your option) any later version.
+ *
+ * This program 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 Affero General Public License for more details.
+ *
+ * You should have received a copy of the GNU Affero General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+package org.aavso.tools.vstar.util.period.wwz;
+
+import org.aavso.tools.vstar.util.period.dcdft.DataTestBase;
+
+/**
+ * Benchmark to quantify speedup of closed-form 3×3 matinv over legacy Gauss-Jordan.
+ * Same scenario as WWZTUmi2420000To2425000Test; runs multiple iterations and prints times.
+ */
+public class WWZMatinvBenchmarkTest extends DataTestBase {
+
+    private static final int WARMUP = 1;
+    private static final int ITERATIONS = 5;
+
+    public WWZMatinvBenchmarkTest() {
+        super("WWZ matinv benchmark", TUmi2420000To2425000Data.data);
+    }
+
+    /**
+     * Run WWZ with new (closed-form) matinv then with legacy (Gauss-Jordan), report total times and speedup.
+     */
+    public void testMatinvBenchmark() throws Exception {
+        double minFreq = 0.01;
+        double maxFreq = 0.02;
+        double deltaFreq = 0.001;
+        double decay = 0.01;
+        double timeDivisions = 50.0;
+
+        // Warmup and time: closed-form (new) matinv
+        WeightedWaveletZTransform.useLegacyMatinv = false;
+        runWarmup(obs, decay, timeDivisions, minFreq, maxFreq, deltaFreq);
+        long newNs = runIterations(obs, decay, timeDivisions, minFreq, maxFreq, deltaFreq, ITERATIONS);
+
+        // Warmup and time: legacy Gauss-Jordan matinv
+        WeightedWaveletZTransform.useLegacyMatinv = true;
+        runWarmup(obs, decay, timeDivisions, minFreq, maxFreq, deltaFreq);
+        long legacyNs = runIterations(obs, decay, timeDivisions, minFreq, maxFreq, deltaFreq, ITERATIONS);
+
+        WeightedWaveletZTransform.useLegacyMatinv = false;
+
+        double newMs = newNs / 1_000_000.0;
+        double legacyMs = legacyNs / 1_000_000.0;
+        double speedup = (double) legacyNs / (double) newNs;
+
+        System.out.println("WWZ matinv benchmark (TUMi 2420000–2425000, " + ITERATIONS + " runs each):");
+        System.out.println("  Closed-form (new) matinv: " + String.format("%.2f", newMs) + " ms total");
+        System.out.println("  Legacy Gauss-Jordan     : " + String.format("%.2f", legacyMs) + " ms total");
+        System.out.println("  Speedup (legacy/new): " + String.format("%.2fx", speedup));
+    }
+
+    private void runWarmup(java.util.List<org.aavso.tools.vstar.data.ValidObservation> obs,
+            double decay, double timeDivisions, double minFreq, double maxFreq, double deltaFreq) throws Exception {
+        for (int i = 0; i < WARMUP; i++) {
+            WeightedWaveletZTransform wwt = new WeightedWaveletZTransform(obs, decay, timeDivisions);
+            wwt.make_freqs_from_freq_range(minFreq, maxFreq, deltaFreq);
+            wwt.execute();
+        }
+    }
+
+    private long runIterations(java.util.List<org.aavso.tools.vstar.data.ValidObservation> obs,
+            double decay, double timeDivisions, double minFreq, double maxFreq, double deltaFreq, int n) throws Exception {
+        long start = System.nanoTime();
+        for (int i = 0; i < n; i++) {
+            WeightedWaveletZTransform wwt = new WeightedWaveletZTransform(obs, decay, timeDivisions);
+            wwt.make_freqs_from_freq_range(minFreq, maxFreq, deltaFreq);
+            wwt.execute();
+        }
+        return System.nanoTime() - start;
+    }
+}