Search in sources :

Example 1 with DoubleDoubleBiPredicate

use of uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate in project GDSC-SMLM by aherbert.

the class TrackPopulationAnalysisTest method canComputeFbmDiffusionFunction.

@Test
void canComputeFbmDiffusionFunction() {
    final int size = 10;
    final double delta = 1e-6;
    final DoubleDoubleBiPredicate test = TestHelper.doublesAreClose(1e-5);
    for (final double t : new double[] { 0.8, 1, 1.2 }) {
        final MultivariateJacobianFunction f = new FbmDiffusionFunction(size, t);
        for (final double d : new double[] { 0.8, 0.9, 1, 1.1, 1.2 }) {
            for (final double s : new double[] { 2, 20 }) {
                for (final double alpha : new double[] { 1e-2, 0.8, 0.9, 1, 1.1, 1.2, 2.0 }) {
                    // Check the value and Jacobian
                    final RealVector point = new ArrayRealVector(new double[] { d, s, alpha }, false);
                    final Pair<RealVector, RealMatrix> p = f.value(point);
                    final double[] value = p.getFirst().toArray();
                    Assertions.assertEquals(size, value.length);
                    final double a1 = alpha + 1;
                    final double a2 = alpha + 2;
                    final double ta = Math.pow(t, alpha);
                    for (int n = 1; n <= size; n++) {
                        // MSD = [4Dt^a / (a+2)(a+1)] * [(n+1)^(a+2) + (n-1)^(a+2) - 2n^(a+2)]
                        // - [8Dt^a / (a+2)(a+1)] + 4s^2
                        // assume t=1.
                        final double msd = (4 * d * ta / (a2 * a1)) * (Math.pow(n + 1, a2) + Math.pow(n - 1, a2) - 2 * Math.pow(n, a2)) - 8 * d * ta / (a2 * a1) + 4 * s * s;
                        TestAssertions.assertTest(msd, value[n - 1], test, "value");
                    }
                    // Columns of the Jacobian
                    final double[] dfda1 = p.getSecond().getColumn(0);
                    final double[] dfdb1 = p.getSecond().getColumn(1);
                    final double[] dfdc1 = p.getSecond().getColumn(2);
                    point.setEntry(0, d - delta);
                    RealVector v1 = f.value(point).getFirst();
                    point.setEntry(0, d + delta);
                    RealVector v2 = f.value(point).getFirst();
                    final double[] dfda = v2.subtract(v1).mapDivide(2 * delta).toArray();
                    point.setEntry(0, d);
                    point.setEntry(1, s - delta);
                    v1 = f.value(point).getFirst();
                    point.setEntry(1, s + delta);
                    v2 = f.value(point).getFirst();
                    final double[] dfdb = v2.subtract(v1).mapDivide(2 * delta).toArray();
                    point.setEntry(1, s);
                    point.setEntry(2, alpha - delta);
                    v1 = f.value(point).getFirst();
                    point.setEntry(2, alpha + delta);
                    v2 = f.value(point).getFirst();
                    final double[] dfdc = v2.subtract(v1).mapDivide(2 * delta).toArray();
                    // Element-by-element relative error
                    TestAssertions.assertArrayTest(dfda, dfda1, test, "jacobian dfda");
                    TestAssertions.assertArrayTest(dfdb, dfdb1, test, "jacobian dfdb");
                    TestAssertions.assertArrayTest(dfdc, dfdc1, test, "jacobian dfdc");
                }
            }
        }
    }
}
Also used : DoubleDoubleBiPredicate(uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate) RealMatrix(org.apache.commons.math3.linear.RealMatrix) RealVector(org.apache.commons.math3.linear.RealVector) ArrayRealVector(org.apache.commons.math3.linear.ArrayRealVector) FbmDiffusionFunction(uk.ac.sussex.gdsc.smlm.ij.plugins.TrackPopulationAnalysis.FbmDiffusionFunction) ArrayRealVector(org.apache.commons.math3.linear.ArrayRealVector) MultivariateJacobianFunction(org.apache.commons.math3.fitting.leastsquares.MultivariateJacobianFunction) Test(org.junit.jupiter.api.Test)

Example 2 with DoubleDoubleBiPredicate

use of uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate in project GDSC-SMLM by aherbert.

the class TrackPopulationAnalysisTest method canComputeBrownianModelUsingFbmFunction.

@Test
void canComputeBrownianModelUsingFbmFunction() {
    final int size = 10;
    final DoubleDoubleBiPredicate test = TestHelper.doublesAreClose(1e-5);
    final double alpha = 1.0;
    for (final double t : new double[] { 0.8, 1, 1.2 }) {
        final BrownianDiffusionFunction f1 = new BrownianDiffusionFunction(size, t);
        final FbmDiffusionFunction f2 = new FbmDiffusionFunction(size, t);
        for (final double d : new double[] { 0.8, 0.9, 1, 1.1, 1.2 }) {
            for (final double s : new double[] { 2, 20 }) {
                // Check the value and Jacobian
                final RealVector point1 = new ArrayRealVector(new double[] { d, s, alpha }, false);
                final Pair<RealVector, RealMatrix> p1 = f1.value(point1);
                final double[] value1 = p1.getFirst().toArray();
                final RealVector point2 = new ArrayRealVector(new double[] { d, s, alpha }, false);
                final Pair<RealVector, RealMatrix> p2 = f2.value(point2);
                final double[] value2 = p2.getFirst().toArray();
                TestAssertions.assertArrayTest(value1, value2, test, "value");
                final double[] dfda1 = p1.getSecond().getColumn(0);
                final double[] dfdb1 = p1.getSecond().getColumn(1);
                final double[] dfda2 = p2.getSecond().getColumn(0);
                final double[] dfdb2 = p2.getSecond().getColumn(1);
                TestAssertions.assertArrayTest(dfda1, dfda2, test, "jacobian dfda");
                TestAssertions.assertArrayTest(dfdb1, dfdb2, test, "jacobian dfdb");
            }
        }
    }
}
Also used : BrownianDiffusionFunction(uk.ac.sussex.gdsc.smlm.ij.plugins.TrackPopulationAnalysis.BrownianDiffusionFunction) DoubleDoubleBiPredicate(uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate) RealMatrix(org.apache.commons.math3.linear.RealMatrix) RealVector(org.apache.commons.math3.linear.RealVector) ArrayRealVector(org.apache.commons.math3.linear.ArrayRealVector) FbmDiffusionFunction(uk.ac.sussex.gdsc.smlm.ij.plugins.TrackPopulationAnalysis.FbmDiffusionFunction) ArrayRealVector(org.apache.commons.math3.linear.ArrayRealVector) Test(org.junit.jupiter.api.Test)

Example 3 with DoubleDoubleBiPredicate

use of uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate in project GDSC-SMLM by aherbert.

the class LvmGradientProcedureTest method gradientProcedureComputesSameAsGradientCalculator.

private void gradientProcedureComputesSameAsGradientCalculator(RandomSeed seed, int nparams, Type type, double error) {
    final int iter = 10;
    final double[][] alpha = new double[nparams][nparams];
    final double[] beta = new double[nparams];
    final ArrayList<double[]> paramsList = new ArrayList<>(iter);
    final ArrayList<double[]> yList = new ArrayList<>(iter);
    final int[] x = createFakeData(RngUtils.create(seed.getSeed()), nparams, iter, paramsList, yList);
    final int n = x.length;
    final FakeGradientFunction func = new FakeGradientFunction(blockWidth, nparams);
    final boolean mle = type != Type.LSQ;
    final FastLog fastLog = (type == Type.FAST_LOG_MLE) ? getFastLog() : null;
    final GradientCalculator calc = GradientCalculatorUtils.newCalculator(nparams, mle);
    final String name = String.format("[%d] %b", nparams, mle);
    // Create messages
    final IndexSupplier msgR = new IndexSupplier(1, name + "Result: Not same ", null);
    final IndexSupplier msgOb = new IndexSupplier(1, name + "Observations: Not same beta ", null);
    final IndexSupplier msgOal = new IndexSupplier(1, name + "Observations: Not same alpha linear ", null);
    final IndexSupplier msgOam = new IndexSupplier(1, name + "Observations: Not same alpha matrix ", null);
    final DoubleDoubleBiPredicate predicate = (error == 0) ? TestHelper.doublesEqual() : TestHelper.doublesAreClose(error, 0);
    for (int i = 0; i < paramsList.size(); i++) {
        // Reference implementation
        final double s = calc.findLinearised(n, yList.get(i), paramsList.get(i), alpha, beta, func);
        // Procedure
        final LvmGradientProcedure p = LvmGradientProcedureUtils.create(yList.get(i), func, type, fastLog);
        p.gradient(paramsList.get(i));
        final double s2 = p.value;
        // Value may be different depending on log implementation
        msgR.set(0, i);
        TestAssertions.assertTest(s, s2, predicate, msgR);
        // Exactly the same ...
        Assertions.assertArrayEquals(p.beta, beta, msgOb.set(0, i));
        final double[] al = p.getAlphaLinear();
        Assertions.assertArrayEquals(al, new DenseMatrix64F(alpha).data, msgOal.set(0, i));
        final double[][] am = p.getAlphaMatrix();
        Assertions.assertArrayEquals(am, alpha, msgOam.set(0, i));
    }
}
Also used : DoubleDoubleBiPredicate(uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate) ArrayList(java.util.ArrayList) FastLog(uk.ac.sussex.gdsc.smlm.function.FastLog) DenseMatrix64F(org.ejml.data.DenseMatrix64F) IndexSupplier(uk.ac.sussex.gdsc.test.utils.functions.IndexSupplier) FakeGradientFunction(uk.ac.sussex.gdsc.smlm.function.FakeGradientFunction)

Example 4 with DoubleDoubleBiPredicate

use of uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate in project GDSC-SMLM by aherbert.

the class BaseFunctionSolverTest method fitAndComputeValueMatch.

/**
 * Check the fit and compute values match. The first solver will be used to do the fit. This is
 * initialised from the solution so the convergence criteria can be set to accept the first step.
 * The second solver is used to compute values (thus is not initialised for fitting).
 *
 * @param seed the seed
 * @param solver1 the solver
 * @param solver2 the solver 2
 * @param noiseModel the noise model
 * @param useWeights the use weights
 */
void fitAndComputeValueMatch(RandomSeed seed, BaseFunctionSolver solver1, BaseFunctionSolver solver2, NoiseModel noiseModel, boolean useWeights) {
    final double[] noise = getNoise(seed, noiseModel);
    if (solver1.isWeighted() && useWeights) {
        solver1.setWeights(getWeights(seed, noiseModel));
        solver2.setWeights(getWeights(seed, noiseModel));
    }
    // Draw target data
    final UniformRandomProvider rg = RngUtils.create(seed.getSeed());
    final double[] data = drawGaussian(p12, noise, noiseModel, rg);
    // fit with 2 peaks using the known params.
    final Gaussian2DFunction f2 = GaussianFunctionFactory.create2D(2, size, size, flags, null);
    solver1.setGradientFunction(f2);
    solver2.setGradientFunction(f2);
    double[] params = p12.clone();
    solver1.fit(data, null, params, null);
    solver2.computeValue(data, null, params);
    final DoubleDoubleBiPredicate predicate = TestHelper.doublesAreClose(1e-10, 0);
    double v1 = solver1.getValue();
    double v2 = solver2.getValue();
    TestAssertions.assertTest(v1, v2, predicate, "Fit 2 peaks and computeValue");
    final double[] o1 = new double[f2.size()];
    final double[] o2 = new double[o1.length];
    solver1.fit(data, o1, params, null);
    solver2.computeValue(data, o2, params);
    v1 = solver1.getValue();
    v2 = solver2.getValue();
    TestAssertions.assertTest(v1, v2, predicate, "Fit 2 peaks and computeValue with yFit");
    final StandardValueProcedure p = new StandardValueProcedure();
    double[] expected = p.getValues(f2, params);
    Assertions.assertArrayEquals(expected, o1, 1e-8, "Fit 2 peaks yFit");
    Assertions.assertArrayEquals(expected, o2, 1e-8, "computeValue 2 peaks yFit");
    if (solver1 instanceof SteppingFunctionSolver) {
        // fit with 1 peak + 1 precomputed using the known params.
        // compare to 2 peak computation.
        final ErfGaussian2DFunction f1 = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(1, size, size, flags, null);
        final Gradient2Function pf1 = OffsetGradient2Function.wrapGradient2Function(f1, p2v);
        solver1.setGradientFunction(pf1);
        solver2.setGradientFunction(pf1);
        params = p1.clone();
        solver1.fit(data, null, params, null);
        solver2.computeValue(data, null, params);
        v1 = solver1.getValue();
        v2 = solver2.getValue();
        TestAssertions.assertTest(v1, v2, predicate, "Fit 1 peak + 1 precomputed and computeValue");
        Arrays.fill(o1, 0);
        Arrays.fill(o2, 0);
        solver1.fit(data, o1, params, null);
        solver2.computeValue(data, o2, params);
        v1 = solver1.getValue();
        v2 = solver2.getValue();
        TestAssertions.assertTest(v1, v2, predicate, "Fit 1 peak + 1 precomputed and computeValue with yFit");
        expected = p.getValues(pf1, params);
        Assertions.assertArrayEquals(expected, o1, 1e-8, "Fit 1 peak + 1 precomputed yFit");
        Assertions.assertArrayEquals(expected, o2, 1e-8, "computeValue 1 peak + 1 precomputed yFit");
    }
}
Also used : ErfGaussian2DFunction(uk.ac.sussex.gdsc.smlm.function.gaussian.erf.ErfGaussian2DFunction) DoubleDoubleBiPredicate(uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate) ErfGaussian2DFunction(uk.ac.sussex.gdsc.smlm.function.gaussian.erf.ErfGaussian2DFunction) Gaussian2DFunction(uk.ac.sussex.gdsc.smlm.function.gaussian.Gaussian2DFunction) Gradient2Function(uk.ac.sussex.gdsc.smlm.function.Gradient2Function) OffsetGradient2Function(uk.ac.sussex.gdsc.smlm.function.OffsetGradient2Function) StandardValueProcedure(uk.ac.sussex.gdsc.smlm.function.StandardValueProcedure) UniformRandomProvider(org.apache.commons.rng.UniformRandomProvider)

Example 5 with DoubleDoubleBiPredicate

use of uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate in project GDSC-SMLM by aherbert.

the class SolverSpeedTest method solveLinearAndGaussJordanReturnSameSolutionAndInversionResult.

@SeededTest
void solveLinearAndGaussJordanReturnSameSolutionAndInversionResult(RandomSeed seed) {
    final int iter = 100;
    final SolverSpeedTestData data = ensureData(seed, iter);
    final ArrayList<double[][]> adata = copyAdouble(data.adata, iter);
    final ArrayList<double[]> bdata = copyBdouble(data.bdata, iter);
    final ArrayList<double[][]> adata2 = copyAdouble(data.adata, iter);
    final ArrayList<double[]> bdata2 = copyBdouble(data.bdata, iter);
    final GaussJordan solver = new GaussJordan();
    final EjmlLinearSolver solver2 = new EjmlLinearSolver();
    final int failureLimit = TestCounter.computeFailureLimit(iter, 0.1);
    final TestCounter failCounter = new TestCounter(failureLimit, 2);
    final DoubleDoubleBiPredicate predicate = TestHelper.doublesAreClose(1e-2, 0);
    int fail = 0;
    for (int i = 0; i < adata.size(); i++) {
        final double[][] a1 = adata.get(i);
        final double[] b1 = bdata.get(i);
        final double[][] a2 = adata2.get(i);
        final double[] b2 = bdata2.get(i);
        final boolean r1 = solver.solve(a1, b1);
        final boolean r2 = solver2.solveLinear(a2, b2);
        solver2.invertLastA(a2);
        // Assertions.assertTrue("Different solve result @ " + i, r1 == r2);
        if (r1 && r2) {
            failCounter.run(0, () -> {
                TestAssertions.assertArrayTest(b1, b2, predicate, "Different b result");
            });
            failCounter.run(0, () -> {
                TestAssertions.assertArrayTest(a1, a2, predicate, "Different a result");
            });
        } else {
            fail++;
        }
    }
    if (fail > iter / 2) {
        Assertions.fail(String.format("Failed to solve %d / %d", fail, iter));
    }
}
Also used : TestCounter(uk.ac.sussex.gdsc.test.utils.TestCounter) DoubleDoubleBiPredicate(uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate) SeededTest(uk.ac.sussex.gdsc.test.junit5.SeededTest)

Aggregations

DoubleDoubleBiPredicate (uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate)63 SeededTest (uk.ac.sussex.gdsc.test.junit5.SeededTest)27 Test (org.junit.jupiter.api.Test)22 UniformRandomProvider (org.apache.commons.rng.UniformRandomProvider)12 PoissonDistribution (org.apache.commons.math3.distribution.PoissonDistribution)6 TDoubleArrayList (gnu.trove.list.array.TDoubleArrayList)4 ArrayList (java.util.ArrayList)4 MixtureMultivariateGaussianDistribution (uk.ac.sussex.gdsc.smlm.math3.distribution.fitting.MultivariateGaussianMixtureExpectationMaximization.MixtureMultivariateGaussianDistribution)4 MultivariateGaussianDistribution (uk.ac.sussex.gdsc.smlm.math3.distribution.fitting.MultivariateGaussianMixtureExpectationMaximization.MixtureMultivariateGaussianDistribution.MultivariateGaussianDistribution)4 UnivariateFunction (org.apache.commons.math3.analysis.UnivariateFunction)3 SimpsonIntegrator (org.apache.commons.math3.analysis.integration.SimpsonIntegrator)3 ArrayRealVector (org.apache.commons.math3.linear.ArrayRealVector)3 RealMatrix (org.apache.commons.math3.linear.RealMatrix)3 RealVector (org.apache.commons.math3.linear.RealVector)3 ParameterizedTest (org.junit.jupiter.params.ParameterizedTest)3 BigDecimal (java.math.BigDecimal)2 MathContext (java.math.MathContext)2 MixtureMultivariateNormalDistribution (org.apache.commons.math3.distribution.MixtureMultivariateNormalDistribution)2 MultivariateNormalDistribution (org.apache.commons.math3.distribution.MultivariateNormalDistribution)2 MultivariateJacobianFunction (org.apache.commons.math3.fitting.leastsquares.MultivariateJacobianFunction)2