Examples with DoubleDoubleBiPredicate uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate used on opensource projects

Example 16 with DoubleDoubleBiPredicate

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

the class PoissonGammaFunctionTest method canComputePoissonGammaGradient.

@SuppressWarnings("unused")
private static void canComputePoissonGammaGradient(final double gain, final double mu, boolean nonInteger) {
    final double o = mu;
    // * o;
    final double delta = 1e-3;
    final double uo = o + delta;
    final double lo = o - delta;
    final double diff = uo - lo;
    // The numerical gradient is poor around the switch between the use of the
    // Bessel function and the approximation. So just count the errors.
    int fail = 0;
    double sum = 0;
    final int[] range = PoissonGaussianFunctionTest.getRange(gain, mu, 0);
    final int min = Math.max(0, range[0]);
    final int max = range[1];
    final double[] dp_dt = new double[1];
    final double[] dp_dt2 = new double[1];
    final double step = (nonInteger) ? 0.5 : 1;
    // When using the approximation the gradients are not as accurate
    final boolean approx = (2 * Math.sqrt(max * o / gain) > 709);
    final double tol = approx ? 0.05 : 1e-3;
    final TDoubleArrayList list = new TDoubleArrayList();
    if (min != 0) {
        list.add(0);
    }
    for (double x = min; x <= max; x += step) {
        list.add(x);
    }
    final DoubleDoubleBiPredicate predicate = TestHelper.doublesAreClose(1e-8, 0);
    for (final double x : list.toArray()) {
        final double p1 = PoissonGammaFunction.poissonGamma(x, o, gain);
        final double p2 = PoissonGammaFunction.poissonGamma(x, o, gain, dp_dt);
        Assertions.assertEquals(p1, p2);
        // Check partial gradient matches
        double p3 = PoissonGammaFunction.poissonGammaPartial(x, o, gain, dp_dt2);
        Assertions.assertEquals(p1, p3);
        TestAssertions.assertTest(dp_dt[0] + p1, dp_dt2[0], predicate);
        // Check no dirac gradient matches
        p3 = PoissonGammaFunction.poissonGammaN(x, o, gain, dp_dt2);
        if (x == 0) {
            final double dirac = PoissonGammaFunction.dirac(o);
            // Add the dirac contribution
            p3 += dirac;
            dp_dt2[0] -= dirac;
            TestAssertions.assertTest(p1, p3, predicate);
            TestAssertions.assertTest(dp_dt[0], dp_dt2[0], predicate);
        } else {
            Assertions.assertEquals(p1, p3);
            TestAssertions.assertTest(dp_dt[0], dp_dt2[0], predicate);
        }
        final double up = PoissonGammaFunction.poissonGamma(x, uo, gain);
        final double lp = PoissonGammaFunction.poissonGamma(x, lo, gain);
        final double eg = dp_dt[0];
        final double g = (up - lp) / diff;
        final double error = DoubleEquality.relativeError(g, eg);
        final double ox = x / gain;
        if (error > tol) {
            fail++;
            sum += error;
        }
    }
    final double f = (double) fail / list.size();
    logger.log(TestLogUtils.getRecord(Level.INFO, "g=%g, mu=%g, failures=%g, mean=%f", gain, mu, f, MathUtils.div0(sum, fail)));
    if (approx) {
        Assertions.assertTrue(f < 0.2);
    } else {
        Assertions.assertTrue(f < 0.01);
    }
}

Example 17 with DoubleDoubleBiPredicate

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

the class PoissonFunctionTest method probabilityMatchesPoissonWithNoGain.

private static void probabilityMatchesPoissonWithNoGain(final double mu) {
    final double o = mu;
    final PoissonFunction f = new PoissonFunction(1.0);
    final PoissonDistribution pd = new PoissonDistribution(mu);
    final double p = 0;
    final int[] range = getRange(1, mu);
    final int min = range[0];
    final int max = range[1];
    final DoubleDoubleBiPredicate predicate = TestHelper.doublesAreClose(1e-8, 0);
    for (int x = min; x <= max; x++) {
        final double v1 = f.likelihood(x, o);
        final double v2 = pd.probability(x);
        TestAssertions.assertTest(v1, v2, predicate, FunctionUtils.getSupplier("g=%f, mu=%f, x=%d", gain, mu, x));
    }
}

Example 18 with DoubleDoubleBiPredicate

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

the class PoissonGammaGaussianConvolutionFunctionTest method probabilityMatchesLogProbability.

private static void probabilityMatchesLogProbability(final double gain, double mu, double sd) {
    final PoissonGammaGaussianConvolutionFunction f = PoissonGammaGaussianConvolutionFunction.createWithStandardDeviation(1.0 / gain, sd);
    // Evaluate an initial range.
    // Gaussian should have >99% within +/- s
    // Poisson will have mean mu with a variance mu.
    // At large mu it is approximately normal so use 3 sqrt(mu) for the range added to the mean
    // Note: The input s parameter is after-gain so adjust.
    final int[] range = PoissonGaussianFunctionTest.getRange(gain, mu, sd / gain);
    final int min = range[0];
    final int max = range[1];
    // Note: The input mu parameter is pre-gain.
    final double e = mu;
    final Supplier<String> msg = () -> String.format("g=%f, mu=%f, s=%f", gain, mu, sd);
    final DoubleDoubleBiPredicate predicate = TestHelper.doublesAreClose(1e-3, 0);
    for (int x = min; x <= max; x++) {
        final double p = f.likelihood(x, e);
        if (p == 0) {
            continue;
        }
        final double logP = f.logLikelihood(x, e);
        TestAssertions.assertTest(Math.log(p), logP, predicate, msg);
    }
}

Example 19 with DoubleDoubleBiPredicate

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

the class PoissonGaussianFunctionTest method probabilityMatchesLogProbability.

private static void probabilityMatchesLogProbability(final double gain, final double mu, final double sd, final boolean usePicard) {
    // Note: The input mu & s parameters are pre-gain.
    final PoissonGaussianFunction f = PoissonGaussianFunction.createWithStandardDeviation(1.0 / gain, mu, sd * gain);
    f.setUsePicardApproximation(usePicard);
    // Evaluate an initial range.
    // Gaussian should have >99% within +/- s
    // Poisson will have mean mu with a variance mu.
    // At large mu it is approximately normal so use 3 sqrt(mu) for the range added to the mean
    final int[] range = getRange(gain, mu, sd);
    final int min = range[0];
    final int max = range[1];
    final Supplier<String> msg = () -> String.format("g=%f, mu=%f, s=%f", gain, mu, sd);
    final DoubleDoubleBiPredicate predicate = TestHelper.doublesAreClose(1e-3, 0);
    for (int x = min; x <= max; x++) {
        final double p = f.probability(x);
        if (p == 0) {
            continue;
        }
        final double logP = f.logProbability(x);
        TestAssertions.assertTest(Math.log(p), logP, predicate, msg);
    }
}

Example 20 with DoubleDoubleBiPredicate

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

the class ScmosLikelihoodWrapperTest method instanceLikelihoodMatches.

private static void instanceLikelihoodMatches(final double mu, boolean test) {
    // Determine upper limit for a Poisson
    final int limit = new PoissonDistribution(mu).inverseCumulativeProbability(P_LIMIT);
    // Map to observed values using the gain and offset
    final double max = limit * G;
    final double step = 0.1;
    final int n = (int) Math.ceil(max / step);
    // Evaluate all values from (zero+offset) to large n
    final double[] k = SimpleArrayUtils.newArray(n, O, step);
    final double[] a = new double[0];
    final double[] gradient = new double[0];
    final float[] var = newArray(n, VAR);
    final float[] g = newArray(n, G);
    final float[] o = newArray(n, O);
    final NonLinearFunction nlf = new NonLinearFunction() {

        @Override
        public void initialise(double[] a) {
        // Ignore
        }

        @Override
        public int[] gradientIndices() {
            return new int[0];
        }

        @Override
        public double evalw(int x, double[] dyda, double[] weight) {
            return 0;
        }

        @Override
        public double evalw(int x, double[] weight) {
            return 0;
        }

        @Override
        public double eval(int x) {
            return mu;
        }

        @Override
        public double eval(int x, double[] dyda) {
            return mu;
        }

        @Override
        public boolean canComputeWeights() {
            return false;
        }

        @Override
        public int getNumberOfGradients() {
            return 0;
        }
    };
    ScmosLikelihoodWrapper func = new ScmosLikelihoodWrapper(nlf, a, k, n, var, g, o);
    final IntArrayFormatSupplier msg1 = new IntArrayFormatSupplier("computeLikelihood @ %d", 1);
    final IntArrayFormatSupplier msg2 = new IntArrayFormatSupplier("computeLikelihood+gradient @ %d", 1);
    double total = 0;
    double pvalue = 0;
    double maxp = 0;
    int maxi = 0;
    final DoubleDoubleBiPredicate predicate = TestHelper.doublesAreClose(1e-10, 0);
    for (int i = 0; i < n; i++) {
        final double nll = func.computeLikelihood(i);
        final double nll2 = func.computeLikelihood(gradient, i);
        final double nll3 = ScmosLikelihoodWrapper.negativeLogLikelihood(mu, var[i], g[i], o[i], k[i]);
        total += nll;
        TestAssertions.assertTest(nll3, nll, predicate, msg1.set(0, i));
        TestAssertions.assertTest(nll3, nll2, predicate, msg2.set(0, i));
        final double pp = StdMath.exp(-nll);
        if (maxp < pp) {
            maxp = pp;
            maxi = i;
        // TestLog.fine(logger,"mu=%f, e=%f, k=%f, pp=%f", mu, mu * G + O, k[i], pp);
        }
        pvalue += pp * step;
    }
    // Expected max of the distribution is the mode of the Poisson distribution.
    // This has two modes for integer input counts. We take the mean of those.
    // https://en.wikipedia.org/wiki/Poisson_distribution
    // Note that the shift of VAR/(G*G) is a constant applied to both the expected and
    // observed values and consequently cancels when predicting the max, i.e. we add
    // a constant count to the observed values and shift the distribution by the same
    // constant. We can thus compute the mode for the unshifted distribution.
    final double lambda = mu;
    final double mode1 = Math.floor(lambda);
    final double mode2 = Math.ceil(lambda) - 1;
    // Scale to observed values
    final double kmax = ((mode1 + mode2) * 0.5) * G + O;
    // TestLog.fine(logger,"mu=%f, p=%f, maxp=%f @ %f (expected=%f %f)", mu, p, maxp, k[maxi], kmax,
    // kmax - k[maxi]);
    TestAssertions.assertTest(kmax, k[maxi], TestHelper.doublesAreClose(1e-3, 0), "k-max");
    if (test) {
        Assertions.assertEquals(P_LIMIT, pvalue, 0.02, () -> "mu=" + mu);
    }
    // Check the function can compute the same total
    double sum;
    double sum2;
    sum = func.computeLikelihood();
    sum2 = func.computeLikelihood(gradient);
    TestAssertions.assertTest(total, sum, predicate, "computeLikelihood");
    TestAssertions.assertTest(total, sum2, predicate, "computeLikelihood with gradient");
    // Check the function can compute the same total after duplication
    func = func.build(nlf, a);
    sum = func.computeLikelihood();
    sum2 = func.computeLikelihood(gradient);
    TestAssertions.assertTest(total, sum, predicate, "computeLikelihood");
    TestAssertions.assertTest(total, sum2, predicate, "computeLikelihood with gradient");
}