Examples with DoubleEquality gdsc.core.utils.DoubleEquality used on opensource projects

Example 6 with DoubleEquality

use of gdsc.core.utils.DoubleEquality in project GDSC-SMLM by aherbert.

the class LSQLVMGradientProcedureTest method gradientProcedureComputesGradient.

private void gradientProcedureComputesGradient(ErfGaussian2DFunction func) {
    int nparams = func.getNumberOfGradients();
    int[] indices = func.gradientIndices();
    int iter = 100;
    rdg = new RandomDataGenerator(new Well19937c(30051977));
    ArrayList<double[]> paramsList = new ArrayList<double[]>(iter);
    ArrayList<double[]> yList = new ArrayList<double[]>(iter);
    createData(1, iter, paramsList, yList, true);
    double delta = 1e-3;
    DoubleEquality eq = new DoubleEquality(1e-3, 1e-3);
    for (int i = 0; i < paramsList.size(); i++) {
        double[] y = yList.get(i);
        double[] a = paramsList.get(i);
        double[] a2 = a.clone();
        BaseLSQLVMGradientProcedure p = LSQLVMGradientProcedureFactory.create(y, func);
        p.gradient(a);
        //double s = p.ssx;
        double[] beta = p.beta.clone();
        for (int j = 0; j < nparams; j++) {
            int k = indices[j];
            double d = Precision.representableDelta(a[k], (a[k] == 0) ? 1e-3 : a[k] * delta);
            a2[k] = a[k] + d;
            p.value(a2);
            double s1 = p.value;
            a2[k] = a[k] - d;
            p.value(a2);
            double s2 = p.value;
            a2[k] = a[k];
            // Apply a factor of -2 to compute the actual gradients:
            // See Numerical Recipes in C++, 2nd Ed. Equation 15.5.6 for Nonlinear Models 
            beta[j] *= -2;
            double gradient = (s1 - s2) / (2 * d);
            //System.out.printf("[%d,%d] %f  (%s %f+/-%f)  %f  ?=  %f\n", i, k, s, func.getName(k), a[k], d, beta[j],
            //		gradient);
            Assert.assertTrue("Not same gradient @ " + j, eq.almostEqualRelativeOrAbsolute(beta[j], gradient));
        }
    }
}

Example 7 with DoubleEquality

use of gdsc.core.utils.DoubleEquality in project GDSC-SMLM by aherbert.

the class LVMGradientProcedureTest method gradientProcedureSupportsPrecomputed.

private void gradientProcedureSupportsPrecomputed(final Type type) {
    int iter = 10;
    rdg = new RandomDataGenerator(new Well19937c(30051977));
    ArrayList<double[]> paramsList = new ArrayList<double[]>(iter);
    ArrayList<double[]> yList = new ArrayList<double[]>(iter);
    // 3 peaks
    createData(3, iter, paramsList, yList, true);
    for (int i = 0; i < paramsList.size(); i++) {
        final double[] y = yList.get(i);
        // Add Gaussian read noise so we have negatives
        double min = Maths.min(y);
        for (int j = 0; j < y.length; j++) y[j] = y[i] - min + rdg.nextGaussian(0, Noise);
    }
    // We want to know that:
    // y|peak1+peak2+peak3 == y|peak1+peak2+peak3(precomputed)
    // We want to know when:
    // y|peak1+peak2+peak3 != y-peak3|peak1+peak2
    // i.e. we cannot subtract a precomputed peak from the data, it must be included in the fit
    // E.G. LSQ - subtraction is OK, MLE/WLSQ - subtraction is not allowed
    Gaussian2DFunction f123 = GaussianFunctionFactory.create2D(3, blockWidth, blockWidth, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    Gaussian2DFunction f12 = GaussianFunctionFactory.create2D(2, blockWidth, blockWidth, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    Gaussian2DFunction f3 = GaussianFunctionFactory.create2D(1, blockWidth, blockWidth, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    int nparams = f12.getNumberOfGradients();
    int[] indices = f12.gradientIndices();
    final double[] b = new double[f12.size()];
    double delta = 1e-3;
    DoubleEquality eq = new DoubleEquality(5e-3, 1e-6);
    double[] a1peaks = new double[7];
    final double[] y_b = new double[b.length];
    for (int i = 0; i < paramsList.size(); i++) {
        final double[] y = yList.get(i);
        double[] a3peaks = paramsList.get(i);
        double[] a2peaks = Arrays.copyOf(a3peaks, 1 + 2 * 6);
        double[] a2peaks2 = a2peaks.clone();
        for (int j = 1; j < 7; j++) a1peaks[j] = a3peaks[j + 2 * 6];
        // Evaluate peak 3 to get the background and subtract it from the data to get the new data
        f3.initialise0(a1peaks);
        f3.forEach(new ValueProcedure() {

            int k = 0;

            public void execute(double value) {
                b[k] = value;
                // Remove negatives for MLE
                if (type == Type.MLE) {
                    y[k] = Math.max(0, y[k]);
                    y_b[k] = Math.max(0, y[k] - value);
                } else {
                    y_b[k] = y[k] - value;
                }
                k++;
            }
        });
        // These should be the same
        LVMGradientProcedure p123 = LVMGradientProcedureFactory.create(y, f123, type);
        LVMGradientProcedure p12b3 = LVMGradientProcedureFactory.create(y, PrecomputedGradient1Function.wrapGradient1Function(f12, b), type);
        // This may be different
        LVMGradientProcedure p12m3 = LVMGradientProcedureFactory.create(y_b, f12, type);
        // Check they are the same
        p123.gradient(a3peaks);
        double[][] m123 = p123.getAlphaMatrix();
        p12b3.gradient(a2peaks);
        double s = p12b3.value;
        double[] beta = p12b3.beta.clone();
        double[][] alpha = p12b3.getAlphaMatrix();
        System.out.printf("MLE=%b [%d] p12b3  %f  %f\n", type, i, p123.value, s);
        Assert.assertTrue("p12b3 Not same value @ " + i, eq.almostEqualRelativeOrAbsolute(p123.value, s));
        Assert.assertTrue("p12b3 Not same gradient @ " + i, eq.almostEqualRelativeOrAbsolute(beta, p123.beta));
        for (int j = 0; j < alpha.length; j++) Assert.assertTrue("p12b3 Not same alpha @ " + j, eq.almostEqualRelativeOrAbsolute(alpha[j], m123[j]));
        // Check actual gradients are correct
        for (int j = 0; j < nparams; j++) {
            int k = indices[j];
            double d = Precision.representableDelta(a2peaks[k], (a2peaks[k] == 0) ? 1e-3 : a2peaks[k] * delta);
            a2peaks2[k] = a2peaks[k] + d;
            p12b3.value(a2peaks2);
            double s1 = p12b3.value;
            a2peaks2[k] = a2peaks[k] - d;
            p12b3.value(a2peaks2);
            double s2 = p12b3.value;
            a2peaks2[k] = a2peaks[k];
            // Apply a factor of -2 to compute the actual gradients:
            // See Numerical Recipes in C++, 2nd Ed. Equation 15.5.6 for Nonlinear Models
            beta[j] *= -2;
            double gradient = (s1 - s2) / (2 * d);
            System.out.printf("[%d,%d] %f  (%s %f+/-%f)  %f  ?=  %f  (%f)\n", i, k, s, f12.getName(k), a2peaks[k], d, beta[j], gradient, DoubleEquality.relativeError(gradient, beta[j]));
            Assert.assertTrue("Not same gradient @ " + j, eq.almostEqualRelativeOrAbsolute(beta[j], gradient));
        }
        // Check these may be different
        p12m3.gradient(a2peaks);
        s = p12m3.value;
        beta = p12m3.beta.clone();
        alpha = p12m3.getAlphaMatrix();
        System.out.printf("%s [%d] p12m3  %f  %f\n", type, i, p123.value, s);
        if (type != Type.LSQ) {
            Assert.assertFalse("p12b3 Same value @ " + i, eq.almostEqualRelativeOrAbsolute(p123.value, s));
            Assert.assertFalse("p12b3 Same gradient @ " + i, eq.almostEqualRelativeOrAbsolute(beta, p123.beta));
            for (int j = 0; j < alpha.length; j++) {
                //System.out.printf("%s != %s\n", Arrays.toString(alpha[j]), Arrays.toString(m123[j]));
                Assert.assertFalse("p12b3 Same alpha @ " + j, eq.almostEqualRelativeOrAbsolute(alpha[j], m123[j]));
            }
        } else {
            Assert.assertTrue("p12b3 Not same value @ " + i, eq.almostEqualRelativeOrAbsolute(p123.value, s));
            Assert.assertTrue("p12b3 Not same gradient @ " + i, eq.almostEqualRelativeOrAbsolute(beta, p123.beta));
            for (int j = 0; j < alpha.length; j++) Assert.assertTrue("p12b3 Not same alpha @ " + j, eq.almostEqualRelativeOrAbsolute(alpha[j], m123[j]));
        }
        // Check actual gradients are correct
        for (int j = 0; j < nparams; j++) {
            int k = indices[j];
            double d = Precision.representableDelta(a2peaks[k], (a2peaks[k] == 0) ? 1e-3 : a2peaks[k] * delta);
            a2peaks2[k] = a2peaks[k] + d;
            p12m3.value(a2peaks2);
            double s1 = p12m3.value;
            a2peaks2[k] = a2peaks[k] - d;
            p12m3.value(a2peaks2);
            double s2 = p12m3.value;
            a2peaks2[k] = a2peaks[k];
            // Apply a factor of -2 to compute the actual gradients:
            // See Numerical Recipes in C++, 2nd Ed. Equation 15.5.6 for Nonlinear Models
            beta[j] *= -2;
            double gradient = (s1 - s2) / (2 * d);
            System.out.printf("[%d,%d] %f  (%s %f+/-%f)  %f  ?=  %f  (%f)\n", i, k, s, f12.getName(k), a2peaks[k], d, beta[j], gradient, DoubleEquality.relativeError(gradient, beta[j]));
            Assert.assertTrue("Not same gradient @ " + j, eq.almostEqualRelativeOrAbsolute(beta[j], gradient));
        }
    }
}

Example 8 with DoubleEquality

use of gdsc.core.utils.DoubleEquality in project GDSC-SMLM by aherbert.

the class SCMOSLikelihoodWrapperTest method fitNBEllipticalComputesGradient.

@Test
public void fitNBEllipticalComputesGradient() {
    // The elliptical function gradient evaluation is worse 
    DoubleEquality tmp = eq;
    eq = eqPerDatum;
    functionComputesGradient(GaussianFunctionFactory.FIT_SIMPLE_NB_ELLIPTICAL);
    eq = tmp;
}

Example 9 with DoubleEquality

use of gdsc.core.utils.DoubleEquality in project GDSC-SMLM by aherbert.

the class SpeedTest method f1ComputesSameAsf2.

void f1ComputesSameAsf2(int npeaks, int flags1, int flags2) {
    DoubleEquality eq = new DoubleEquality(1e-2, 1e-10);
    int iter = 2000;
    ArrayList<double[]> paramsList2 = (npeaks == 1) ? copyList(paramsListSinglePeak, iter) : copyList(paramsListDoublePeak, iter);
    Gaussian2DFunction f1 = GaussianFunctionFactory.create2D(1, blockWidth, blockWidth, flags1, null);
    Gaussian2DFunction f2 = GaussianFunctionFactory.create2D(1, blockWidth, blockWidth, flags2, null);
    double[] dyda1 = new double[1 + npeaks * 6];
    double[] dyda2 = new double[1 + npeaks * 6];
    int[] gradientIndices = f1.gradientIndices();
    int[] g1 = new int[gradientIndices.length];
    int[] g2 = new int[gradientIndices.length];
    int nparams = 0;
    for (int i = 0; i < gradientIndices.length; i++) {
        int index1 = f1.findGradientIndex(g1[i]);
        int index2 = f2.findGradientIndex(g2[i]);
        if (index1 >= 0 && index2 >= 0) {
            g1[nparams] = index1;
            g2[nparams] = index2;
            nparams++;
        }
    }
    for (int i = 0; i < paramsList2.size(); i++) {
        f1.initialise(paramsList2.get(i));
        f2.initialise(paramsList2.get(i));
        for (int j = 0; j < x.length; j++) {
            double y1 = f1.eval(x[j], dyda1);
            double y2 = f2.eval(x[j], dyda2);
            Assert.assertTrue("Not same y[" + j + "] @ " + i + " " + y1 + " != " + y2, eq.almostEqualRelativeOrAbsolute(y1, y2));
            for (int ii = 0; ii < nparams; ii++) Assert.assertTrue("Not same dyda[" + j + "] @ " + gradientIndices[g1[ii]] + ": " + dyda1[g1[ii]] + " != " + dyda2[g2[ii]], eq.almostEqualRelativeOrAbsolute(dyda1[g1[ii]], dyda2[g2[ii]]));
        }
    }
}

Example 10 with DoubleEquality

use of gdsc.core.utils.DoubleEquality in project GDSC-SMLM by aherbert.

the class FastMLEJacobianGradient2ProcedureTest method gradientCalculatorComputesGradient.

private void gradientCalculatorComputesGradient(int nPeaks, ErfGaussian2DFunction func) {
    // Check the first and second derivatives
    int nparams = func.getNumberOfGradients();
    int[] indices = func.gradientIndices();
    int iter = 100;
    rdg = new RandomDataGenerator(new Well19937c(30051977));
    ArrayList<double[]> paramsList = new ArrayList<double[]>(iter);
    ArrayList<double[]> yList = new ArrayList<double[]>(iter);
    createData(nPeaks, iter, paramsList, yList, true);
    double delta = 1e-5;
    DoubleEquality eq = new DoubleEquality(1e-2, 1e-3);
    for (int i = 0; i < paramsList.size(); i++) {
        double[] y = yList.get(i);
        double[] a = paramsList.get(i);
        double[] a2 = a.clone();
        FastMLEJacobianGradient2Procedure p = new FastMLEJacobianGradient2Procedure(y, (ExtendedGradient2Function) func);
        //double ll = p.computeLogLikelihood(a);
        p.computeJacobian(a);
        double[] d1 = p.d1.clone();
        double[] d2 = p.d2.clone();
        DenseMatrix64F J = DenseMatrix64F.wrap(nparams, nparams, p.getJacobianLinear());
        for (int j = 0; j < nparams; j++) {
            int k = indices[j];
            double d = Precision.representableDelta(a[k], (a[k] == 0) ? delta : a[k] * delta);
            a2[k] = a[k] + d;
            double llh = p.computeLogLikelihood(a2);
            p.computeFirstDerivative(a2);
            double[] d1h = p.d1.clone();
            a2[k] = a[k] - d;
            double lll = p.computeLogLikelihood(a2);
            p.computeFirstDerivative(a2);
            double[] d1l = p.d1.clone();
            a2[k] = a[k];
            double gradient1 = (llh - lll) / (2 * d);
            double gradient2 = (d1h[j] - d1l[j]) / (2 * d);
            //System.out.printf("[%d,%d] ll - %f  (%s %f+/-%f) d1 %f ?= %f : d2 %f ?= %f\n", i, k, ll, func.getName(k), a[k], d, 
            //		gradient1, d1[j], gradient2, d2[j]);
            Assert.assertTrue("Not same gradient1 @ " + j, eq.almostEqualRelativeOrAbsolute(gradient1, d1[j]));
            Assert.assertTrue("Not same gradient2 @ " + j, eq.almostEqualRelativeOrAbsolute(gradient2, d2[j]));
            for (int jj = 0; jj < nparams; jj++) {
                if (j == jj) {
                // This is done above
                // Check it anyway to ensure the Jacobian is correct
                //continue;
                }
                int kk = indices[jj];
                double dd = Precision.representableDelta(a[kk], (a[kk] == 0) ? delta : a[kk] * delta);
                a2[kk] = a[kk] + dd;
                p.computeFirstDerivative(a2);
                d1h = p.d1.clone();
                a2[kk] = a[kk] - dd;
                p.computeFirstDerivative(a2);
                d1l = p.d1.clone();
                a2[kk] = a[kk];
                // Use index j even though we adjusted index jj
                gradient2 = (d1h[j] - d1l[j]) / (2 * dd);
                boolean ok = eq.almostEqualRelativeOrAbsolute(gradient2, J.get(j, jj));
                //		a[k], func.getName(kk), a[kk], dd, gradient2, J.get(j, jj), ok);
                if (!ok) {
                    Assert.fail(String.format("Not same gradientJ @ [%d,%d]", j, jj));
                }
            }
        }
    }
}