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");
}
}
}
}
}
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");
}
}
}
}
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));
}
}
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");
}
}
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));
}
}
Aggregations