diff --git a/.github/workflows/ci-publish.yml b/.github/workflows/ci-publish.yml index b10f004f..06c6b812 100644 --- a/.github/workflows/ci-publish.yml +++ b/.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 }} diff --git a/beast-base/pom.xml b/beast-base/pom.xml index c010deb0..5f02f9ea 100644 --- a/beast-base/pom.xml +++ b/beast-base/pom.xml @@ -91,11 +91,7 @@ io.github.compevol beagle 1.0 - - - io.github.compevol - colt - 1.0 + true diff --git a/beast-base/src/main/java/beast/base/evolution/substitutionmodel/ColtEigenSystem.java b/beast-base/src/main/java/beast/base/evolution/substitutionmodel/ColtEigenSystem.java index b209721c..7bc1a1e6 100644 --- a/beast-base/src/main/java/beast/base/evolution/substitutionmodel/ColtEigenSystem.java +++ b/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; diff --git a/beast-base/src/main/java/beast/base/evolution/substitutionmodel/ComplexColtEigenSystem.java b/beast-base/src/main/java/beast/base/evolution/substitutionmodel/ComplexColtEigenSystem.java index 58745cdb..c63c87dd 100644 --- a/beast-base/src/main/java/beast/base/evolution/substitutionmodel/ComplexColtEigenSystem.java +++ b/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"); } diff --git a/beast-base/src/main/java/beast/base/math/matrixalgebra/RobustEigenDecomposition.java b/beast-base/src/main/java/beast/base/math/matrixalgebra/RobustEigenDecomposition.java index 89859110..d3177e4a 100644 --- a/beast-base/src/main/java/beast/base/math/matrixalgebra/RobustEigenDecomposition.java +++ b/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(); diff --git a/beast-base/src/main/java/beast/base/math/matrixalgebra/RobustSingularValueDecomposition.java b/beast-base/src/main/java/beast/base/math/matrixalgebra/RobustSingularValueDecomposition.java index 4cd71d83..dbd653ce 100644 --- a/beast-base/src/main/java/beast/base/math/matrixalgebra/RobustSingularValueDecomposition.java +++ b/beast-base/src/main/java/beast/base/math/matrixalgebra/RobustSingularValueDecomposition.java @@ -1,9 +1,6 @@ package beast.base.math.matrixalgebra; import beast.base.math.MathUtils; -import cern.colt.matrix.DoubleFactory2D; -import cern.colt.matrix.DoubleMatrix2D; -import cern.colt.matrix.linalg.Property; public class RobustSingularValueDecomposition implements java.io.Serializable { @@ -26,8 +23,6 @@ public class RobustSingularValueDecomposition implements java.io.Serializable { */ private int m, n; -// private int maxIterations; - static private int maxIterationsDefault = 100000; private static final String ERROR_STRING = "SVD is not converged."; @@ -35,23 +30,23 @@ public class RobustSingularValueDecomposition implements java.io.Serializable { /** Constructs and returns a new singular value decomposition object; The decomposed matrices can be retrieved via instance methods of the returned decomposition object. - @param Arg A rectangular matrix. Returns a decomposition object to access {@code U}, {@code S} and {@code V}. - @throws IllegalArgumentException if {@code A.rows() < A.columns()}. + @param Arg A rectangular matrix as a 2D array. + @throws IllegalArgumentException if {@code Arg} is empty. */ - public RobustSingularValueDecomposition(DoubleMatrix2D Arg) throws ArithmeticException { + public RobustSingularValueDecomposition(double[][] Arg) throws ArithmeticException { this(Arg,maxIterationsDefault); } - public RobustSingularValueDecomposition(DoubleMatrix2D Arg, int maxIterations) throws ArithmeticException { - Property.DEFAULT.checkRectangular(Arg); - -// this.maxIterations = maxIterations; + public RobustSingularValueDecomposition(double[][] Arg, int maxIterations) throws ArithmeticException { // Derived from LINPACK code. // Initialize. - double[][] A = Arg.toArray(); - m = Arg.rows(); - n = Arg.columns(); + double[][] A = new double[Arg.length][]; + for (int i = 0; i < Arg.length; i++) { + A[i] = Arg[i].clone(); + } + m = Arg.length; + n = Arg[0].length; int nu = Math.min(m,n); s = new double [Math.min(m+1,n)]; U = new double [m][nu]; @@ -470,7 +465,7 @@ public double cond() { Returns the diagonal matrix of singular values. @return S */ - public DoubleMatrix2D getS() { + public double[][] getS() { double[][] S = new double[n][n]; for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { @@ -478,7 +473,7 @@ public DoubleMatrix2D getS() { } S[i][i] = this.s[i]; } - return DoubleFactory2D.dense.make(S); + return S; } /** Returns the diagonal of {@code S}, which is a one-dimensional array of singular values @@ -491,16 +486,23 @@ public double[] getSingularValues() { Returns the left singular vectors {@code U}. @return {@code U} */ - public DoubleMatrix2D getU() { - //return new DoubleMatrix2D(U,m,Math.min(m+1,n)); - return DoubleFactory2D.dense.make(U).viewPart(0,0,m,Math.min(m+1,n)); + public double[][] getU() { + double[][] result = new double[m][Math.min(m+1,n)]; + for (int i = 0; i < m; i++) { + System.arraycopy(U[i], 0, result[i], 0, Math.min(m+1,n)); + } + return result; } /** Returns the right singular vectors {@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; } /** Returns the two norm, which is {@code max(S)}. @@ -552,20 +554,17 @@ public String toString() { catch (IllegalArgumentException exc) { buf.append(unknown+exc.getMessage()); } buf.append("\n\nU = "); - try { buf.append(String.valueOf(this.getU()));} + try { buf.append(java.util.Arrays.deepToString(this.getU()));} catch (IllegalArgumentException exc) { buf.append(unknown+exc.getMessage()); } buf.append("\n\nS = "); - try { buf.append(String.valueOf(this.getS()));} + try { buf.append(java.util.Arrays.deepToString(this.getS()));} 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(); } } - - - diff --git a/beast-base/src/main/java/module-info.java b/beast-base/src/main/java/module-info.java index 981af86f..22876638 100644 --- a/beast-base/src/main/java/module-info.java +++ b/beast-base/src/main/java/module-info.java @@ -20,7 +20,6 @@ // Automatic modules requires org.antlr.antlr4.runtime; requires beagle; - requires colt; // Export all packages exports beast.base; diff --git a/beast-base/src/test/java/beast/base/math/matrixalgebra/RobustEigenDecompositionTest.java b/beast-base/src/test/java/beast/base/math/matrixalgebra/RobustEigenDecompositionTest.java index 22c6218e..642d5a6a 100644 --- a/beast-base/src/test/java/beast/base/math/matrixalgebra/RobustEigenDecompositionTest.java +++ b/beast-base/src/test/java/beast/base/math/matrixalgebra/RobustEigenDecompositionTest.java @@ -1,8 +1,5 @@ package beast.base.math.matrixalgebra; -import cern.colt.matrix.DoubleFactory2D; -import cern.colt.matrix.DoubleMatrix2D; -import cern.colt.matrix.linalg.Algebra; import org.junit.jupiter.api.Test; import java.util.Arrays; @@ -19,20 +16,18 @@ */ public class RobustEigenDecompositionTest { - private static final Algebra ALG = new Algebra(); - // ---- Symmetric matrices with known eigenvalues ---- @Test public void testSymmetric2x2() { // [[3, 1], [1, 3]] has eigenvalues 2, 4 double[][] A = {{3, 1}, {1, 3}}; - RobustEigenDecomposition eig = decompose(A); + RobustEigenDecomposition eig = new RobustEigenDecomposition(A); - double[] lambda = sorted(eig.getRealEigenvalues().toArray()); + double[] lambda = sorted(eig.getRealEigenvalues()); assertEquals(2.0, lambda[0], 1e-12); assertEquals(4.0, lambda[1], 1e-12); - assertAllZero(eig.getImagEigenvalues().toArray(), 1e-12); + assertAllZero(eig.getImagEigenvalues(), 1e-12); assertEigenProperty(A, eig, 1e-12); } @@ -40,13 +35,13 @@ public void testSymmetric2x2() { public void testSymmetric3x3Tridiagonal() { // Eigenvalues: 2 - sqrt(2), 2, 2 + sqrt(2) double[][] A = {{2, -1, 0}, {-1, 2, -1}, {0, -1, 2}}; - RobustEigenDecomposition eig = decompose(A); + RobustEigenDecomposition eig = new RobustEigenDecomposition(A); - double[] lambda = sorted(eig.getRealEigenvalues().toArray()); + double[] lambda = sorted(eig.getRealEigenvalues()); assertEquals(2 - Math.sqrt(2), lambda[0], 1e-12); assertEquals(2.0, lambda[1], 1e-12); assertEquals(2 + Math.sqrt(2), lambda[2], 1e-12); - assertAllZero(eig.getImagEigenvalues().toArray(), 1e-12); + assertAllZero(eig.getImagEigenvalues(), 1e-12); assertEigenProperty(A, eig, 1e-12); } @@ -54,10 +49,10 @@ public void testSymmetric3x3Tridiagonal() { public void testSymmetricOrthogonalEigenvectors() { // For symmetric matrices, V should be orthogonal: V^T * V = I double[][] A = {{4, 1, 2}, {1, 3, 0}, {2, 0, 5}}; - RobustEigenDecomposition eig = decompose(A); + RobustEigenDecomposition eig = new RobustEigenDecomposition(A); - DoubleMatrix2D V = eig.getV(); - DoubleMatrix2D VtV = ALG.mult(ALG.transpose(V), V); + double[][] V = eig.getV(); + double[][] VtV = multiply(transpose(V), V); assertIdentity(VtV, 1e-12); assertEigenProperty(A, eig, 1e-12); } @@ -75,7 +70,7 @@ public void testNonSymmetric4x4RateMatrix() { 2.0, 0.25, 0.25 }; double[][] Q = buildNormalizedRateMatrix(f, r); - RobustEigenDecomposition eig = decompose(Q); + RobustEigenDecomposition eig = new RobustEigenDecomposition(Q); assertEigenProperty(Q, eig, 1e-10); } @@ -90,7 +85,7 @@ public void testExtremeRates4x4() { 1E6, 0.0, 0.0 }; double[][] Q = buildNormalizedRateMatrix(f, r); - RobustEigenDecomposition eig = decompose(Q); + RobustEigenDecomposition eig = new RobustEigenDecomposition(Q); assertEigenProperty(Q, eig, 1e-6); } @@ -162,16 +157,12 @@ public void testNonSymmetric20x20RateMatrix() { 0.009942, 0.023276, 0.003999, 0.010515, 0.013447, 0.022690, 0.012800, 0.019909, 0.038235 }; double[][] Q = buildNormalizedRateMatrix(f, r); - RobustEigenDecomposition eig = decompose(Q); + RobustEigenDecomposition eig = new RobustEigenDecomposition(Q); assertEigenProperty(Q, eig, 1e-8); } // ---- Helpers ---- - private static RobustEigenDecomposition decompose(double[][] data) { - return new RobustEigenDecomposition(DoubleFactory2D.dense.make(data)); - } - private static double[] sorted(double[] values) { double[] copy = values.clone(); Arrays.sort(copy); @@ -181,17 +172,16 @@ private static double[] sorted(double[] values) { /** * Verify A * V = V * D (the fundamental eigendecomposition property). */ - private static void assertEigenProperty(double[][] data, RobustEigenDecomposition eig, double tol) { - DoubleMatrix2D A = DoubleFactory2D.dense.make(data); - DoubleMatrix2D V = eig.getV(); - DoubleMatrix2D D = eig.getD(); - DoubleMatrix2D AV = ALG.mult(A, V); - DoubleMatrix2D VD = ALG.mult(V, D); - - int n = A.rows(); + private static void assertEigenProperty(double[][] A, RobustEigenDecomposition eig, double tol) { + double[][] V = eig.getV(); + double[][] D = eig.getD(); + double[][] AV = multiply(A, V); + double[][] VD = multiply(V, D); + + int n = A.length; for (int i = 0; i < n; i++) for (int j = 0; j < n; j++) - assertEquals(AV.get(i, j), VD.get(i, j), tol, + assertEquals(AV[i][j], VD[i][j], tol, "A*V != V*D at [" + i + "][" + j + "]"); } @@ -200,14 +190,33 @@ private static void assertAllZero(double[] values, double tol) { assertEquals(0.0, values[i], tol, "Expected zero imaginary part at index " + i); } - private static void assertIdentity(DoubleMatrix2D M, double tol) { - int n = M.rows(); + private static void assertIdentity(double[][] M, double tol) { + int n = M.length; for (int i = 0; i < n; i++) for (int j = 0; j < n; j++) - assertEquals(i == j ? 1.0 : 0.0, M.get(i, j), tol, + assertEquals(i == j ? 1.0 : 0.0, M[i][j], tol, "Not identity at [" + i + "][" + j + "]"); } + static double[][] multiply(double[][] A, double[][] B) { + int m = A.length, n = B[0].length, p = B.length; + double[][] C = new double[m][n]; + for (int i = 0; i < m; i++) + for (int j = 0; j < n; j++) + for (int k = 0; k < p; k++) + C[i][j] += A[i][k] * B[k][j]; + return C; + } + + static double[][] transpose(double[][] A) { + int m = A.length, n = A[0].length; + double[][] T = new double[n][m]; + for (int i = 0; i < m; i++) + for (int j = 0; j < n; j++) + T[j][i] = A[i][j]; + return T; + } + /** * Build a normalized rate matrix following transprob.R conventions. * Rates are ordered: for each row i, first columns j < i, then j > i. diff --git a/beast-base/src/test/java/beast/base/math/matrixalgebra/RobustSingularValueDecompositionTest.java b/beast-base/src/test/java/beast/base/math/matrixalgebra/RobustSingularValueDecompositionTest.java index 3fd3464c..201faf6c 100644 --- a/beast-base/src/test/java/beast/base/math/matrixalgebra/RobustSingularValueDecompositionTest.java +++ b/beast-base/src/test/java/beast/base/math/matrixalgebra/RobustSingularValueDecompositionTest.java @@ -1,10 +1,8 @@ package beast.base.math.matrixalgebra; -import cern.colt.matrix.DoubleFactory2D; -import cern.colt.matrix.DoubleMatrix2D; -import cern.colt.matrix.linalg.Algebra; import org.junit.jupiter.api.Test; +import static beast.base.math.matrixalgebra.RobustEigenDecompositionTest.*; import static org.junit.jupiter.api.Assertions.*; /** @@ -14,13 +12,11 @@ */ public class RobustSingularValueDecompositionTest { - private static final Algebra ALG = new Algebra(); - @Test public void testDiagonal3x3() { // Singular values of a diagonal matrix are the absolute values of diagonal entries double[][] A = {{3, 0, 0}, {0, -5, 0}, {0, 0, 2}}; - RobustSingularValueDecomposition svd = decompose(A); + RobustSingularValueDecomposition svd = new RobustSingularValueDecomposition(A); double[] sv = svd.getSingularValues(); assertEquals(5.0, sv[0], 1e-12); @@ -34,7 +30,7 @@ public void testKnownSingularValues2x2() { // A^T*A = [[25, -15], [-15, 25]], eigenvalues 40, 10 // Singular values: sqrt(40), sqrt(10) double[][] A = {{4, 0}, {3, -5}}; - RobustSingularValueDecomposition svd = decompose(A); + RobustSingularValueDecomposition svd = new RobustSingularValueDecomposition(A); double[] sv = svd.getSingularValues(); assertEquals(Math.sqrt(40), sv[0], 1e-10); @@ -50,8 +46,8 @@ public void testReconstructionSquare() { {0.4, 6.1, 2.2, 3.7}, {5.5, 1.0, 3.3, 0.6} }; - RobustSingularValueDecomposition svd = decompose(A); - assertReconstruction(DoubleFactory2D.dense.make(A), svd, 1e-10); + RobustSingularValueDecomposition svd = new RobustSingularValueDecomposition(A); + assertReconstruction(A, svd, 1e-10); } @Test @@ -63,25 +59,25 @@ public void testReconstructionRectangular() { {7, 8, 10}, {2, 1, 4} }; - RobustSingularValueDecomposition svd = decompose(A); - assertReconstruction(DoubleFactory2D.dense.make(A), svd, 1e-10); + RobustSingularValueDecomposition svd = new RobustSingularValueDecomposition(A); + assertReconstruction(A, svd, 1e-10); } @Test public void testOrthogonalV() { // V^T * V = I double[][] A = {{1, 2, 3}, {4, 5, 6}, {7, 8, 10}}; - RobustSingularValueDecomposition svd = decompose(A); + RobustSingularValueDecomposition svd = new RobustSingularValueDecomposition(A); - DoubleMatrix2D V = svd.getV(); - DoubleMatrix2D VtV = ALG.mult(ALG.transpose(V), V); + double[][] V = svd.getV(); + double[][] VtV = multiply(transpose(V), V); assertIdentity(VtV, 1e-12); } @Test public void testRankFullRank() { double[][] A = {{1, 0, 0}, {0, 2, 0}, {0, 0, 3}}; - RobustSingularValueDecomposition svd = decompose(A); + RobustSingularValueDecomposition svd = new RobustSingularValueDecomposition(A); assertEquals(3, svd.rank()); } @@ -89,7 +85,7 @@ public void testRankFullRank() { public void testRankDeficient() { // All rows are multiples of [1, 2, 3] — rank 1 double[][] A = {{1, 2, 3}, {2, 4, 6}, {3, 6, 9}}; - RobustSingularValueDecomposition svd = decompose(A); + RobustSingularValueDecomposition svd = new RobustSingularValueDecomposition(A); assertEquals(1, svd.rank()); } @@ -101,7 +97,7 @@ public void testSingularValuesNonNegativeDescending() { {2.4, 1.7, 5.0, 0.3}, {1.9, 3.3, 0.6, 2.8} }; - RobustSingularValueDecomposition svd = decompose(A); + RobustSingularValueDecomposition svd = new RobustSingularValueDecomposition(A); double[] sv = svd.getSingularValues(); for (int i = 0; i < sv.length; i++) { @@ -116,7 +112,7 @@ public void testSingularValuesNonNegativeDescending() { public void testConditionNumber() { // cond = max(S) / min(S) double[][] A = {{3, 0, 0}, {0, 6, 0}, {0, 0, 2}}; - RobustSingularValueDecomposition svd = decompose(A); + RobustSingularValueDecomposition svd = new RobustSingularValueDecomposition(A); assertEquals(6.0 / 2.0, svd.cond(), 1e-12); } @@ -124,13 +120,12 @@ public void testConditionNumber() { public void testNorm2() { // norm2 = max singular value double[][] A = {{3, 0, 0}, {0, 6, 0}, {0, 0, 2}}; - RobustSingularValueDecomposition svd = decompose(A); + RobustSingularValueDecomposition svd = new RobustSingularValueDecomposition(A); assertEquals(6.0, svd.norm2(), 1e-12); } @Test public void test4x4RateMatrix() { - // Use the same rate matrix as RobustEigenDecompositionTest double[] f = {0.25, 0.25, 0.25, 0.25}; double[] r = { 0.0, 1.0, 1.0, @@ -138,39 +133,33 @@ public void test4x4RateMatrix() { 0.25, 0.75, 1.0, 2.0, 0.25, 0.25 }; - double[][] Q = RobustEigenDecompositionTest.buildNormalizedRateMatrix(f, r); - DoubleMatrix2D A = DoubleFactory2D.dense.make(Q); - RobustSingularValueDecomposition svd = new RobustSingularValueDecomposition(A); - assertReconstruction(A, svd, 1e-10); + double[][] Q = buildNormalizedRateMatrix(f, r); + RobustSingularValueDecomposition svd = new RobustSingularValueDecomposition(Q); + assertReconstruction(Q, svd, 1e-10); } // ---- Helpers ---- - private static RobustSingularValueDecomposition decompose(double[][] data) { - return new RobustSingularValueDecomposition(DoubleFactory2D.dense.make(data)); - } - /** * Verify A = U * S * V^T. */ - private static void assertReconstruction(DoubleMatrix2D A, RobustSingularValueDecomposition svd, double tol) { - DoubleMatrix2D U = svd.getU(); - DoubleMatrix2D S = svd.getS(); - DoubleMatrix2D V = svd.getV(); - DoubleMatrix2D US = ALG.mult(U, S); - DoubleMatrix2D reconstructed = ALG.mult(US, ALG.transpose(V)); - - for (int i = 0; i < A.rows(); i++) - for (int j = 0; j < A.columns(); j++) - assertEquals(A.get(i, j), reconstructed.get(i, j), tol, + private static void assertReconstruction(double[][] A, RobustSingularValueDecomposition svd, double tol) { + double[][] U = svd.getU(); + double[][] S = svd.getS(); + double[][] V = svd.getV(); + double[][] reconstructed = multiply(multiply(U, S), transpose(V)); + + for (int i = 0; i < A.length; i++) + for (int j = 0; j < A[0].length; j++) + assertEquals(A[i][j], reconstructed[i][j], tol, "A != U*S*V^T at [" + i + "][" + j + "]"); } - private static void assertIdentity(DoubleMatrix2D M, double tol) { - int n = M.rows(); + private static void assertIdentity(double[][] M, double tol) { + int n = M.length; for (int i = 0; i < n; i++) for (int j = 0; j < n; j++) - assertEquals(i == j ? 1.0 : 0.0, M.get(i, j), tol, + assertEquals(i == j ? 1.0 : 0.0, M[i][j], tol, "Not identity at [" + i + "][" + j + "]"); } } diff --git a/beast-base/src/test/java/test/beast/evolution/substmodel/BinaryCovarionModelTest.java b/beast-base/src/test/java/test/beast/evolution/substmodel/BinaryCovarionModelTest.java index b24f9211..fe9c24ff 100644 --- a/beast-base/src/test/java/test/beast/evolution/substmodel/BinaryCovarionModelTest.java +++ b/beast-base/src/test/java/test/beast/evolution/substmodel/BinaryCovarionModelTest.java @@ -12,7 +12,7 @@ import beast.base.spec.type.RealVector; import beast.base.spec.type.Simplex; import beast.base.util.Randomizer; -import cern.colt.Arrays; +import java.util.Arrays; import static org.junit.jupiter.api.Assertions.assertEquals; public class BinaryCovarionModelTest { diff --git a/beast-fx/pom.xml b/beast-fx/pom.xml index b938df31..37847658 100644 --- a/beast-fx/pom.xml +++ b/beast-fx/pom.xml @@ -52,11 +52,7 @@ io.github.compevol beagle 1.0 - - - io.github.compevol - colt - 1.0 + true diff --git a/beast-fx/src/main/java/beastfx/app/treeannotator/KernelDensityEstimator2D.java b/beast-fx/src/main/java/beastfx/app/treeannotator/KernelDensityEstimator2D.java index 2efd7b60..00a009d2 100644 --- a/beast-fx/src/main/java/beastfx/app/treeannotator/KernelDensityEstimator2D.java +++ b/beast-fx/src/main/java/beastfx/app/treeannotator/KernelDensityEstimator2D.java @@ -3,8 +3,6 @@ import java.util.Arrays; import beast.base.util.DiscreteStatistics; -import cern.colt.list.DoubleArrayList; -import cern.jet.stat.Descriptive; /** @@ -238,18 +236,29 @@ private void setupH() { // } public double bandwidthNRD(double[] in) { - DoubleArrayList inList = new DoubleArrayList(in.length); - for (double d : in) - inList.add(d); - inList.sort(); + double[] sorted = in.clone(); + Arrays.sort(sorted); - final double h = (Descriptive.quantile(inList, 0.75) - Descriptive.quantile(inList, 0.25)) / 1.34; + final double h = (quantile(sorted, 0.75) - quantile(sorted, 0.25)) / 1.34; return 4 * 1.06 * Math.min(Math.sqrt(DiscreteStatistics.variance(in)), h) * Math.pow(in.length, -0.2); } + /** + * Returns the quantile of a sorted array using linear interpolation. + */ + private static double quantile(double[] sorted, double phi) { + int n = sorted.length; + double index = phi * (n - 1); + int lo = (int) index; + int hi = lo + 1; + double frac = index - lo; + if (hi >= n) return sorted[n - 1]; + return sorted[lo] + frac * (sorted[hi] - sorted[lo]); + } + public static void main(String[] arg) { double[] x = {3.4, 1.2, 5.6, 2.2, 3.1}; diff --git a/beast-fx/src/main/java/module-info.java b/beast-fx/src/main/java/module-info.java index d8d034f1..770bd6e8 100644 --- a/beast-fx/src/main/java/module-info.java +++ b/beast-fx/src/main/java/module-info.java @@ -18,7 +18,6 @@ // Automatic modules requires beagle; - requires colt; // Export all beastfx packages exports beastfx.app.beast; diff --git a/lib/colt.jar b/lib/colt.jar deleted file mode 100644 index 8ad0541f..00000000 Binary files a/lib/colt.jar and /dev/null differ diff --git a/lib/colt/module-info.java b/lib/colt/module-info.java deleted file mode 100644 index f14b993a..00000000 --- a/lib/colt/module-info.java +++ /dev/null @@ -1,26 +0,0 @@ -module colt { - exports cern.clhep; - exports cern.colt; - exports cern.colt.bitvector; - exports cern.colt.buffer; - exports cern.colt.function; - exports cern.colt.list; - exports cern.colt.list.adapter; - exports cern.colt.map; - exports cern.colt.matrix; - exports cern.colt.matrix.bench; - exports cern.colt.matrix.doublealgo; - exports cern.colt.matrix.impl; - exports cern.colt.matrix.linalg; - exports cern.colt.matrix.objectalgo; - exports cern.jet.math; - exports cern.jet.random; - exports cern.jet.random.engine; - exports cern.jet.random.sampling; - exports cern.jet.stat; - exports cern.jet.stat.quantile; - exports corejava; - exports hep.aida; - exports hep.aida.bin; - exports hep.aida.ref; -}