Examples with IntArrayFormatSupplier uk.ac.sussex.gdsc.test.utils.functions.IntArrayFormatSupplier used on opensource projects

Example 11 with IntArrayFormatSupplier

use of uk.ac.sussex.gdsc.test.utils.functions.IntArrayFormatSupplier in project GDSC-SMLM by aherbert.

the class FastMleGradient2ProcedureTest method gradientProcedureUnrolledComputesSameAsGradientProcedure.

private void gradientProcedureUnrolledComputesSameAsGradientProcedure(RandomSeed seed, int nparams) {
    final int iter = 10;
    final ArrayList<double[]> paramsList = new ArrayList<>(iter);
    final ArrayList<double[]> yList = new ArrayList<>(iter);
    createFakeData(RngUtils.create(seed.getSeed()), nparams, iter, paramsList, yList);
    final FakeGradientFunction func = new FakeGradientFunction(blockWidth, nparams);
    // Create messages
    final IntArrayFormatSupplier msgLl = getMessage(nparams, "[%d] LL: Not same @ %d");
    final IntArrayFormatSupplier msg1Dv = getMessage(nparams, "[%d] first derivative: Not same value @ %d");
    final IntArrayFormatSupplier msg1Dd1 = getMessage(nparams, "[%d] first derivative: Not same d1 @ %d");
    final IntArrayFormatSupplier msg2Dv = getMessage(nparams, "[%d] second derivative: Not same value @ %d");
    final IntArrayFormatSupplier msg2Dd1 = getMessage(nparams, "[%d] second derivative: Not same d1 @ %d");
    final IntArrayFormatSupplier msg2Dd2 = getMessage(nparams, "[%d] second derivative: Not same d2 @ %d");
    for (int i = 0; i < paramsList.size(); i++) {
        FastMleGradient2Procedure p1 = new FastMleGradient2Procedure(yList.get(i), func);
        FastMleGradient2Procedure p2 = FastMleGradient2ProcedureUtils.createUnrolled(yList.get(i), func);
        final double[] a = paramsList.get(i);
        final double ll1 = p1.computeLogLikelihood(a);
        final double ll2 = p2.computeLogLikelihood(a);
        Assertions.assertEquals(ll1, ll2, msgLl.set(1, i));
        p1 = new FastMleGradient2Procedure(yList.get(i), func);
        p2 = FastMleGradient2ProcedureUtils.createUnrolled(yList.get(i), func);
        p1.computeFirstDerivative(a);
        p2.computeFirstDerivative(a);
        Assertions.assertArrayEquals(p1.u, p2.u, msg1Dv.set(1, i));
        Assertions.assertArrayEquals(p1.d1, p2.d1, msg1Dd1.set(1, i));
        p1 = new FastMleGradient2Procedure(yList.get(i), func);
        p2 = FastMleGradient2ProcedureUtils.createUnrolled(yList.get(i), func);
        p1.computeSecondDerivative(a);
        p2.computeSecondDerivative(a);
        Assertions.assertArrayEquals(p1.u, p2.u, msg2Dv.set(1, i));
        Assertions.assertArrayEquals(p1.d1, p2.d1, msg2Dd1.set(1, i));
        Assertions.assertArrayEquals(p1.d2, p2.d2, msg2Dd2.set(1, i));
    }
}

Example 12 with IntArrayFormatSupplier

use of uk.ac.sussex.gdsc.test.utils.functions.IntArrayFormatSupplier 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");
}

Example 13 with IntArrayFormatSupplier

use of uk.ac.sussex.gdsc.test.utils.functions.IntArrayFormatSupplier in project GDSC-SMLM by aherbert.

the class LvmGradientProcedureTest method gradientProcedureIsFasterUnrolledThanGradientProcedure.

private void gradientProcedureIsFasterUnrolledThanGradientProcedure(RandomSeed seed, final int nparams, final Type type, final boolean precomputed) {
    Assumptions.assumeTrue(TestSettings.allow(TestComplexity.MEDIUM));
    final int iter = 100;
    final ArrayList<double[]> paramsList = new ArrayList<>(iter);
    final ArrayList<double[]> yList = new ArrayList<>(iter);
    createData(RngUtils.create(seed.getSeed()), 1, iter, paramsList, yList);
    // Remove the timing of the function call by creating a dummy function
    final FakeGradientFunction fgf = new FakeGradientFunction(blockWidth, nparams);
    final Gradient1Function func;
    if (precomputed) {
        final double[] b = SimpleArrayUtils.newArray(fgf.size(), 0.1, 1.3);
        func = OffsetGradient1Function.wrapGradient1Function(fgf, b);
    } else {
        func = fgf;
    }
    final FastLog fastLog = type == Type.FAST_LOG_MLE ? getFastLog() : null;
    final IntArrayFormatSupplier msgA = new IntArrayFormatSupplier("A [%d]", 1);
    final IntArrayFormatSupplier msgB = new IntArrayFormatSupplier("B [%d]", 1);
    for (int i = 0; i < paramsList.size(); i++) {
        final LvmGradientProcedure p1 = createProcedure(type, yList.get(i), func, fastLog);
        p1.gradient(paramsList.get(i));
        p1.gradient(paramsList.get(i));
        final LvmGradientProcedure p2 = LvmGradientProcedureUtils.create(yList.get(i), func, type, fastLog);
        p2.gradient(paramsList.get(i));
        p2.gradient(paramsList.get(i));
        // Check they are the same
        Assertions.assertArrayEquals(p1.getAlphaLinear(), p2.getAlphaLinear(), msgA.set(0, i));
        Assertions.assertArrayEquals(p1.beta, p2.beta, msgB.set(0, i));
    }
    // Realistic loops for an optimisation
    final int loops = 15;
    // Run till stable timing
    final Timer t1 = new Timer() {

        @Override
        void run() {
            for (int i = 0, k = 0; i < paramsList.size(); i++) {
                final LvmGradientProcedure p1 = createProcedure(type, yList.get(i), func, fastLog);
                for (int j = loops; j-- > 0; ) {
                    p1.gradient(paramsList.get(k++ % iter));
                }
            }
        }
    };
    final long time1 = t1.getTime();
    final Timer t2 = new Timer(t1.loops) {

        @Override
        void run() {
            for (int i = 0, k = 0; i < paramsList.size(); i++) {
                final LvmGradientProcedure p2 = LvmGradientProcedureUtils.create(yList.get(i), func, type, fastLog);
                for (int j = loops; j-- > 0; ) {
                    p2.gradient(paramsList.get(k++ % iter));
                }
            }
        }
    };
    final long time2 = t2.getTime();
    logger.log(TestLogUtils.getTimingRecord(new TimingResult(() -> String.format("%s, Precomputed=%b : Standard", type, precomputed), time1), new TimingResult(() -> String.format("Unrolled %d", nparams), time2)));
}

Example 14 with IntArrayFormatSupplier

use of uk.ac.sussex.gdsc.test.utils.functions.IntArrayFormatSupplier in project GDSC-SMLM by aherbert.

the class WPoissonGradientProcedureTest method getMessage.

private static IntArrayFormatSupplier getMessage(int nparams, String format) {
    final IntArrayFormatSupplier msg = new IntArrayFormatSupplier(format, 2);
    msg.set(0, nparams);
    return msg;
}

Example 15 with IntArrayFormatSupplier

use of uk.ac.sussex.gdsc.test.utils.functions.IntArrayFormatSupplier in project GDSC-SMLM by aherbert.

the class WPoissonGradientProcedureTest method poissonGradientProcedureComputesSameAsWLsqGradientProcedure.

private void poissonGradientProcedureComputesSameAsWLsqGradientProcedure(RandomSeed seed, int nparams) {
    final double[] var = dataCache.computeIfAbsent(seed, WPoissonGradientProcedureTest::createData);
    final int iter = 10;
    final ArrayList<double[]> paramsList = new ArrayList<>(iter);
    final UniformRandomProvider rng = RngUtils.create(seed.getSeed());
    createFakeParams(rng, nparams, iter, paramsList);
    final FakeGradientFunction func = new FakeGradientFunction(blockWidth, nparams);
    final IntArrayFormatSupplier msgOa = getMessage(nparams, "[%d] Observations: Not same alpha @ %d");
    final IntArrayFormatSupplier msgOal = getMessage(nparams, "[%d] Observations: Not same alpha linear @ %d");
    for (int i = 0; i < paramsList.size(); i++) {
        final double[] y = createFakeData(rng);
        final WPoissonGradientProcedure p1 = WPoissonGradientProcedureUtils.create(y, var, func);
        p1.computeFisherInformation(paramsList.get(i));
        final WLsqLvmGradientProcedure p2 = new WLsqLvmGradientProcedure(y, var, func);
        p2.gradient(paramsList.get(i));
        // Exactly the same ...
        Assertions.assertArrayEquals(p1.data, p2.alpha, msgOa.set(1, i));
        Assertions.assertArrayEquals(p1.getLinear(), p2.getAlphaLinear(), msgOal.set(1, i));
    }
}