Remove Colt dependency, replace with commons-math4-legacy

.github/workflows/ci-publish.yml

@@ -3,6 +3,7 @@ name: CI & Publish
 on:
   push:
     branches: [ master ]
+    tags: [ 'v*' ]
   pull_request:
     branches: [ master ]
 
@@ -18,7 +19,8 @@ jobs:
       - name: Checkout repository
         uses: actions/checkout@v4
 
-      - name: Set up JDK 25 (Azul Zulu)
+      - name: Set up JDK 25 (Azul Zulu) for GitHub Packages
+        if: "!startsWith(github.ref, 'refs/tags/v')"
         uses: actions/setup-java@v4
         with:
           distribution: zulu
@@ -28,6 +30,19 @@ jobs:
           server-username: GITHUB_ACTOR
           server-password: GITHUB_TOKEN
 
+      - name: Set up JDK 25 (Azul Zulu) for Maven Central
+        if: startsWith(github.ref, 'refs/tags/v')
+        uses: actions/setup-java@v4
+        with:
+          distribution: zulu
+          java-version: '25'
+          cache: maven
+          server-id: central
+          server-username: CENTRAL_USERNAME
+          server-password: CENTRAL_PASSWORD
+          gpg-private-key: ${{ secrets.GPG_PRIVATE_KEY }}
+          gpg-passphrase: GPG_PASSPHRASE
+
       - name: Install beagle.jar to local Maven repo
         run: >
           mvn install:install-file
@@ -37,15 +52,6 @@ jobs:
           -Dversion=1.0
           -Dpackaging=jar
 
-      - name: Install colt.jar to local Maven repo
-        run: >
-          mvn install:install-file
-          -Dfile=lib/colt.jar
-          -DgroupId=io.github.compevol
-          -DartifactId=colt
-          -Dversion=1.0
-          -Dpackaging=jar
-
       - name: Build and test
         run: mvn verify
 
@@ -64,23 +70,16 @@ jobs:
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
 
-      - name: Deploy colt.jar to GitHub Packages
-        if: github.event_name == 'push' && github.ref == 'refs/heads/master'
-        continue-on-error: true  # 409 Conflict if already published
-        run: >
-          mvn deploy:deploy-file
-          -Dfile=lib/colt.jar
-          -DgroupId=io.github.compevol
-          -DartifactId=colt
-          -Dversion=1.0
-          -Dpackaging=jar
-          -DrepositoryId=github
-          -Durl=https://maven.pkg.github.com/CompEvol/beast3
-        env:
-          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
-
       - name: Deploy all modules to GitHub Packages
         if: github.event_name == 'push' && github.ref == 'refs/heads/master'
         run: mvn deploy -DskipTests
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+
+      - name: Publish to Maven Central
+        if: startsWith(github.ref, 'refs/tags/v')
+        run: mvn deploy -P release -DskipTests
+        env:
+          CENTRAL_USERNAME: ${{ secrets.CENTRAL_USERNAME }}
+          CENTRAL_PASSWORD: ${{ secrets.CENTRAL_PASSWORD }}
+          GPG_PASSPHRASE: ${{ secrets.GPG_PASSPHRASE }}

beast-base/pom.xml

@@ -91,11 +91,7 @@
             <groupId>io.github.compevol</groupId>
             <artifactId>beagle</artifactId>
             <version>1.0</version>
-        </dependency>
-        <dependency>
-            <groupId>io.github.compevol</groupId>
-            <artifactId>colt</artifactId>
-            <version>1.0</version>
+            <optional>true</optional>
         </dependency>
 
         <!-- Testing -->

beast-base/src/main/java/beast/base/evolution/substitutionmodel/ColtEigenSystem.java

@@ -26,10 +26,9 @@
  */
 
 
-import cern.colt.matrix.DoubleMatrix2D;
-import cern.colt.matrix.impl.DenseDoubleMatrix2D;
-import cern.colt.matrix.linalg.Algebra;
-import cern.colt.matrix.linalg.Property;
+import org.apache.commons.math4.legacy.linear.Array2DRowRealMatrix;
+import org.apache.commons.math4.legacy.linear.LUDecomposition;
+import org.apache.commons.math4.legacy.linear.RealMatrix;
 
 import java.util.Arrays;
 
@@ -56,11 +55,10 @@ public EigenDecomposition decomposeMatrix(double[][] matrix) {
 
         final int stateCount = matrix.length;
 
-        RobustEigenDecomposition eigenDecomp = new RobustEigenDecomposition(
-                new DenseDoubleMatrix2D(matrix), maxIterations);
+        RobustEigenDecomposition eigenDecomp = new RobustEigenDecomposition(matrix, maxIterations);
 
-        DoubleMatrix2D eigenV = eigenDecomp.getV();
-        DoubleMatrix2D eigenVInv;
+        double[][] eigenV = eigenDecomp.getV();
+        double[][] eigenVInv;
 
         if (checkConditioning) {
             RobustSingularValueDecomposition svd;
@@ -76,13 +74,12 @@ public EigenDecomposition decomposeMatrix(double[][] matrix) {
         }
 
         try {
-            eigenVInv = alegbra.inverse(eigenV);
-        } catch (IllegalArgumentException e) {
+            RealMatrix eigenVMatrix = new Array2DRowRealMatrix(eigenV);
+            eigenVInv = new LUDecomposition(eigenVMatrix).getSolver().getInverse().getData();
+        } catch (Exception e) {
             return getEmptyDecomposition(stateCount);
         }
 
-        double[][] Evec = eigenV.toArray();
-        double[][] Ievc = eigenVInv.toArray();
         double[] Eval = getAllEigenValues(eigenDecomp);
 
         if (checkConditioning) {
@@ -99,9 +96,9 @@ public EigenDecomposition decomposeMatrix(double[][] matrix) {
         double[] flatEvec = new double[stateCount * stateCount];
         double[] flatIevc = new double[stateCount * stateCount];
 
-        for (int i = 0; i < Evec.length; i++) {
-            System.arraycopy(Evec[i], 0, flatEvec, i * stateCount, stateCount);
-            System.arraycopy(Ievc[i], 0, flatIevc, i * stateCount, stateCount);
+        for (int i = 0; i < stateCount; i++) {
+            System.arraycopy(eigenV[i], 0, flatEvec, i * stateCount, stateCount);
+            System.arraycopy(eigenVInv[i], 0, flatIevc, i * stateCount, stateCount);
         }
 
         return new EigenDecomposition(flatEvec, flatIevc, Eval);
@@ -193,7 +190,7 @@ public void computeExponential(EigenDecomposition eigen, double distance, double
     }
 
     protected double[] getAllEigenValues(RobustEigenDecomposition decomposition) {
-        return decomposition.getRealEigenvalues().toArray();
+        return decomposition.getRealEigenvalues();
     }
 
     protected double[] getEmptyAllEigenValues(int dim) {
@@ -214,8 +211,7 @@ protected EigenDecomposition getEmptyDecomposition(int dim) {
 
     protected final int stateCount;
 
-    private static final double minProb = Property.DEFAULT.tolerance();
-    private static final Algebra alegbra = new Algebra(minProb);
+    private static final double minProb = 1.0E-12;
 
     public static final boolean defaultCheckConditioning = true;
     public static final int defaultMaxConditionNumber = 1000000;

beast-base/src/main/java/beast/base/evolution/substitutionmodel/ComplexColtEigenSystem.java

@@ -26,8 +26,6 @@
  */
 
 
-import cern.colt.matrix.DoubleMatrix2D;
-
 import java.util.Arrays;
 
 import beast.base.math.matrixalgebra.RobustEigenDecomposition;
@@ -46,8 +44,8 @@ public ComplexColtEigenSystem(int stateCount, boolean checkConditioning, int max
     }
 
     protected double[] getAllEigenValues(RobustEigenDecomposition decomposition) {
-        double[] realEval = decomposition.getRealEigenvalues().toArray();
-        double[] imagEval = decomposition.getImagEigenvalues().toArray();
+        double[] realEval = decomposition.getRealEigenvalues();
+        double[] imagEval = decomposition.getImagEigenvalues();
 
         final int dim = realEval.length;
         double[] merge = new double[2 * dim];
@@ -60,10 +58,6 @@ protected double[] getEmptyAllEigenValues(int dim) {
         return new double[2 * dim];
     }
 
-    protected boolean validDecomposition(DoubleMatrix2D eigenV) {
-        return true;
-    }
-
     public double computeExponential(EigenDecomposition eigen, double distance, int i, int j) {
         throw new RuntimeException("Not yet implemented");
     }

beast-base/src/main/java/beast/base/math/matrixalgebra/RobustEigenDecomposition.java

@@ -2,11 +2,6 @@
 
 
 import beast.base.math.MathUtils;
-import cern.colt.matrix.DoubleFactory1D;
-import cern.colt.matrix.DoubleFactory2D;
-import cern.colt.matrix.DoubleMatrix1D;
-import cern.colt.matrix.DoubleMatrix2D;
-import cern.colt.matrix.linalg.Property;
 
 /**
  * Copyright ? 1999 CERN - European Organization for Nuclear Research.
@@ -63,31 +58,35 @@ public class RobustEigenDecomposition implements java.io.Serializable {
 Constructs and returns a new eigenvalue decomposition object;
 The decomposed matrices can be retrieved via instance methods of the returned decomposition object.
 Checks for symmetry, then constructs the eigenvalue decomposition.
-@param A    A square matrix. Returns a decomposition object to access {@code D} and {@code V}.
+@param A    A square matrix as a 2D array. Returns a decomposition object to access {@code D} and {@code V}.
 @throws IllegalArgumentException if {@code A} is not square.
 */
-public RobustEigenDecomposition(DoubleMatrix2D A) throws ArithmeticException {
+public RobustEigenDecomposition(double[][] A) throws ArithmeticException {
     this(A,maxIterationsDefault);
 }
 
 
-public RobustEigenDecomposition(DoubleMatrix2D A, int maxIterations) throws ArithmeticException {
+public RobustEigenDecomposition(double[][] A, int maxIterations) throws ArithmeticException {
 
-	Property.DEFAULT.checkSquare(A);
+	int rows = A.length;
+	int cols = A[0].length;
+	if (rows != cols) {
+		throw new IllegalArgumentException("Matrix must be square: " + rows + " x " + cols);
+	}
 
     this.maxIterations = maxIterations;
 
-	n = A.columns();
+	n = cols;
 	V = new double[n][n];
 	d = new double[n];
 	e = new double[n];
 
-	issymmetric = Property.DEFAULT.isSymmetric(A);
+	issymmetric = isSymmetric(A);
 
 	if (issymmetric) {
 		for (int i = 0; i < n; i++) {
 			for (int j = 0; j < n; j++) {
-				V[i][j] = A.getQuick(i,j);
+				V[i][j] = A[i][j];
 			}
 		}
 
@@ -104,7 +103,7 @@ public RobustEigenDecomposition(DoubleMatrix2D A, int maxIterations) throws Arit
 
 		for (int j = 0; j < n; j++) {
 			for (int i = 0; i < n; i++) {
-				H[i][j] = A.getQuick(i,j);
+				H[i][j] = A[i][j];
 			}
 		}
 
@@ -115,6 +114,23 @@ public RobustEigenDecomposition(DoubleMatrix2D A, int maxIterations) throws Arit
 		hqr2();
 	}
 }
+
+/**
+ * Check if a matrix is symmetric within default tolerance.
+ */
+private static boolean isSymmetric(double[][] A) {
+	double tolerance = 1.0E-12;
+	int n = A.length;
+	for (int i = 0; i < n; i++) {
+		for (int j = i + 1; j < n; j++) {
+			if (Math.abs(A[i][j] - A[j][i]) > tolerance) {
+				return false;
+			}
+		}
+	}
+	return true;
+}
+
 private void cdiv(double xr, double xi, double yr, double yi) {
 	double r,d;
 	if (Math.abs(yr) > Math.abs(yi)) {
@@ -134,7 +150,7 @@ private void cdiv(double xr, double xi, double yr, double yi) {
 Returns the block diagonal eigenvalue matrix, {@code D}.
 @return     {@code D}
 */
-public DoubleMatrix2D getD() {
+public double[][] getD() {
 	double[][] D = new double[n][n];
 	for (int i = 0; i < n; i++) {
 		for (int j = 0; j < n; j++) {
@@ -148,28 +164,32 @@ else if (e[i] < 0) {
 			D[i][i-1] = e[i];
 		}
 	}
-	return DoubleFactory2D.dense.make(D);
+	return D;
 }
 /**
 Returns the imaginary parts of the eigenvalues.
 @return     imag(diag(D))
 */
-public DoubleMatrix1D getImagEigenvalues () {
-	return DoubleFactory1D.dense.make(e);
+public double[] getImagEigenvalues () {
+	return e.clone();
 }
 /**
 Returns the real parts of the eigenvalues.
 @return     real(diag(D))
 */
-public DoubleMatrix1D getRealEigenvalues () {
-	return DoubleFactory1D.dense.make(d);
+public double[] getRealEigenvalues () {
+	return d.clone();
 }
 /**
 Returns the eigenvector matrix, {@code V}
 @return     {@code V}
 */
-public DoubleMatrix2D getV () {
-	return DoubleFactory2D.dense.make(V);
+public double[][] getV () {
+	double[][] result = new double[n][n];
+	for (int i = 0; i < n; i++) {
+		System.arraycopy(V[i], 0, result[i], 0, n);
+	}
+	return result;
 }
 /**
 Nonsymmetric reduction from Hessenberg to real Schur form.
@@ -724,19 +744,19 @@ public String toString() {
 	buf.append("---------------------------------------------------------------------\n");
 
 	buf.append("realEigenvalues = ");
-	try { buf.append(String.valueOf(this.getRealEigenvalues()));}
+	try { buf.append(java.util.Arrays.toString(this.getRealEigenvalues()));}
 	catch (IllegalArgumentException exc) { buf.append(unknown+exc.getMessage()); }
 
 	buf.append("\nimagEigenvalues = ");
-	try { buf.append(String.valueOf(this.getImagEigenvalues()));}
+	try { buf.append(java.util.Arrays.toString(this.getImagEigenvalues()));}
 	catch (IllegalArgumentException exc) { buf.append(unknown+exc.getMessage()); }
 
 	buf.append("\n\nD = ");
-	try { buf.append(String.valueOf(this.getD()));}
+	try { buf.append(java.util.Arrays.deepToString(this.getD()));}
 	catch (IllegalArgumentException exc) { buf.append(unknown+exc.getMessage()); }
 
 	buf.append("\n\nV = ");
-	try { buf.append(String.valueOf(this.getV()));}
+	try { buf.append(java.util.Arrays.deepToString(this.getV()));}
 	catch (IllegalArgumentException exc) { buf.append(unknown+exc.getMessage()); }
 
 	return buf.toString();