Search in sources :

Example 6 with DenseMatrix64F

use of org.ejml.data.DenseMatrix64F in project GDSC-SMLM by aherbert.

the class EJMLLinearSolver method solve.

/**
	 * Solves (one) linear equation, a x = b
	 * <p>
	 * On output b replaced by x. Matrix a may be modified.
	 * 
	 * @return False if the equation is singular (no solution)
	 */
private boolean solve(LinearSolver<DenseMatrix64F> solver, DenseMatrix64F A, DenseMatrix64F B) {
    boolean copy = (errorChecking || solver.modifiesB());
    DenseMatrix64F x = (copy) ? getX() : B;
    if (!solve(solver, A, B, x))
        return false;
    // Copy back result if necessary
    if (copy)
        System.arraycopy(x.data, 0, B.data, 0, B.numRows);
    return true;
}
Also used : DenseMatrix64F(org.ejml.data.DenseMatrix64F)

Example 7 with DenseMatrix64F

use of org.ejml.data.DenseMatrix64F in project GDSC-SMLM by aherbert.

the class EJMLLinearSolverTest method runInversionSpeedTest.

private void runInversionSpeedTest(int flags) {
    final Gaussian2DFunction f0 = GaussianFunctionFactory.create2D(1, 10, 10, flags, null);
    int n = f0.size();
    final double[] y = new double[n];
    final TurboList<DenseMatrix64F> aList = new TurboList<DenseMatrix64F>();
    double[] testbackground = new double[] { 0.2, 0.7 };
    double[] testsignal1 = new double[] { 30, 100, 300 };
    double[] testcx1 = new double[] { 4.9, 5.3 };
    double[] testcy1 = new double[] { 4.8, 5.2 };
    double[] testw1 = new double[] { 1.1, 1.2, 1.5 };
    int np = f0.getNumberOfGradients();
    GradientCalculator calc = GradientCalculatorFactory.newCalculator(np);
    final RandomDataGenerator rdg = new RandomDataGenerator(new Well19937c(30051977));
    //double lambda = 10;
    for (double background : testbackground) // Peak 1
    for (double signal1 : testsignal1) for (double cx1 : testcx1) for (double cy1 : testcy1) for (double w1 : testw1) {
        double[] p = new double[] { background, signal1, 0, cx1, cy1, w1, w1 };
        f0.initialise(p);
        f0.forEach(new ValueProcedure() {

            int i = 0;

            public void execute(double value) {
                // Poisson data 
                y[i++] = rdg.nextPoisson(value);
            }
        });
        double[][] alpha = new double[np][np];
        double[] beta = new double[np];
        //double ss = 
        calc.findLinearised(n, y, p, alpha, beta, f0);
        //System.out.printf("SS = %f\n", ss);
        // As per the LVM algorithm
        //for (int i = 0; i < np; i++)
        //	alpha[i][i] *= lambda;
        aList.add(EJMLLinearSolver.toA(alpha));
    }
    DenseMatrix64F[] a = aList.toArray(new DenseMatrix64F[aList.size()]);
    boolean[] ignore = new boolean[a.length];
    double[][] answer = new double[a.length][];
    int runs = 100000 / a.length;
    TimingService ts = new TimingService(runs);
    TurboList<InversionTimingTask> tasks = new TurboList<InversionTimingTask>();
    tasks.add(new PseudoInverseInversionTimingTask(a, ignore, answer));
    tasks.add(new CholeskyInversionTimingTask(a, ignore, answer));
    tasks.add(new LinearInversionTimingTask(a, ignore, answer));
    tasks.add(new CholeskyLDLTInversionTimingTask(a, ignore, answer));
    tasks.add(new DirectInversionInversionTimingTask(a, ignore, answer));
    tasks.add(new DiagonalDirectInversionInversionTimingTask(a, ignore, answer));
    for (InversionTimingTask task : tasks) if (!task.badSolver)
        ts.execute(task);
    ts.repeat();
    ts.report();
}
Also used : ValueProcedure(gdsc.smlm.function.ValueProcedure) TurboList(gdsc.core.utils.TurboList) RandomDataGenerator(org.apache.commons.math3.random.RandomDataGenerator) Well19937c(org.apache.commons.math3.random.Well19937c) DenseMatrix64F(org.ejml.data.DenseMatrix64F) Gaussian2DFunction(gdsc.smlm.function.gaussian.Gaussian2DFunction) GradientCalculator(gdsc.smlm.fitting.nonlinear.gradient.GradientCalculator) TimingService(gdsc.core.test.TimingService)

Example 8 with DenseMatrix64F

use of org.ejml.data.DenseMatrix64F in project GDSC-SMLM by aherbert.

the class ChiSquaredDistributionTableTest method canComputeChiSquared.

@Test
public void canComputeChiSquared() {
    // We have to use the transpose of the table
    DenseMatrix64F m = new DenseMatrix64F(chi2);
    CommonOps.transpose(m);
    int max = m.numCols;
    double[] et = m.data;
    for (int i = 0, j = 0; i < p.length; i++) {
        ChiSquaredDistributionTable upperTable = ChiSquaredDistributionTable.createUpperTailed(p[i], max);
        // Use 1-p as the significance level to get the same critical values
        ChiSquaredDistributionTable lowerTable = ChiSquaredDistributionTable.createLowerTailed(1 - p[i], max);
        for (int df = 1; df <= max; df++) {
            double o = upperTable.getCrititalValue(df);
            double e = et[j++];
            //System.out.printf("p=%.3f,df=%d = %f\n", p[i], df, o);
            Assert.assertEquals(e, o, 1e-2);
            // The test only stores 2 decimal places so use the computed value to set upper/lower
            double upper = o * 1.01;
            double lower = o * 0.99;
            Assert.assertTrue("Upper did not reject higher", upperTable.reject(upper, df));
            Assert.assertFalse("Upper did not reject actual value", upperTable.reject(o, df));
            Assert.assertFalse("Upper did not accept lower", upperTable.reject(lower, df));
            Assert.assertTrue("Lower did not reject lower", lowerTable.reject(lower, df));
            Assert.assertFalse("Lower did not accept actual value", lowerTable.reject(o, df));
            Assert.assertFalse("Lower did not accept higher", lowerTable.reject(upper, df));
        }
    }
}
Also used : DenseMatrix64F(org.ejml.data.DenseMatrix64F) ChiSquareTest(org.apache.commons.math3.stat.inference.ChiSquareTest) Test(org.junit.Test)

Example 9 with DenseMatrix64F

use of org.ejml.data.DenseMatrix64F in project GDSC-SMLM by aherbert.

the class LSQLVMGradientProcedureTest method gradientProcedureComputesSameAsGradientCalculator.

private void gradientProcedureComputesSameAsGradientCalculator(int nparams, BaseLSQLVMGradientProcedureFactory factory) {
    int iter = 10;
    rdg = new RandomDataGenerator(new Well19937c(30051977));
    double[][] alpha = new double[nparams][nparams];
    double[] beta = new double[nparams];
    ArrayList<double[]> paramsList = new ArrayList<double[]>(iter);
    ArrayList<double[]> yList = new ArrayList<double[]>(iter);
    int[] x = createFakeData(nparams, iter, paramsList, yList);
    int n = x.length;
    FakeGradientFunction func = new FakeGradientFunction(blockWidth, nparams);
    GradientCalculator calc = GradientCalculatorFactory.newCalculator(nparams, false);
    String name = factory.getClass().getSimpleName();
    for (int i = 0; i < paramsList.size(); i++) {
        BaseLSQLVMGradientProcedure p = factory.createProcedure(yList.get(i), func);
        p.gradient(paramsList.get(i));
        double s = p.value;
        double s2 = calc.findLinearised(n, yList.get(i), paramsList.get(i), alpha, beta, func);
        // Exactly the same ...
        Assert.assertEquals(name + " Result: Not same @ " + i, s, s2, 0);
        Assert.assertArrayEquals(name + " Observations: Not same beta @ " + i, p.beta, beta, 0);
        double[] al = p.getAlphaLinear();
        Assert.assertArrayEquals(name + " Observations: Not same alpha @ " + i, al, new DenseMatrix64F(alpha).data, 0);
        double[][] am = p.getAlphaMatrix();
        for (int j = 0; j < nparams; j++) Assert.assertArrayEquals(name + " Observations: Not same alpha @ " + i, am[j], alpha[j], 0);
    }
}
Also used : RandomDataGenerator(org.apache.commons.math3.random.RandomDataGenerator) ArrayList(java.util.ArrayList) Well19937c(org.apache.commons.math3.random.Well19937c) DenseMatrix64F(org.ejml.data.DenseMatrix64F) FakeGradientFunction(gdsc.smlm.function.FakeGradientFunction)

Example 10 with DenseMatrix64F

use of org.ejml.data.DenseMatrix64F in project GDSC-SMLM by aherbert.

the class LVMGradientProcedureTest method gradientProcedureComputesSameAsGradientCalculator.

private void gradientProcedureComputesSameAsGradientCalculator(int nparams, boolean mle) {
    int iter = 10;
    rdg = new RandomDataGenerator(new Well19937c(30051977));
    double[][] alpha = new double[nparams][nparams];
    double[] beta = new double[nparams];
    ArrayList<double[]> paramsList = new ArrayList<double[]>(iter);
    ArrayList<double[]> yList = new ArrayList<double[]>(iter);
    int[] x = createFakeData(nparams, iter, paramsList, yList);
    int n = x.length;
    FakeGradientFunction func = new FakeGradientFunction(blockWidth, nparams);
    GradientCalculator calc = GradientCalculatorFactory.newCalculator(nparams, mle);
    String name = String.format("[%d] %b", nparams, mle);
    for (int i = 0; i < paramsList.size(); i++) {
        LVMGradientProcedure p = LVMGradientProcedureFactory.create(yList.get(i), func, (mle) ? Type.MLE : Type.LSQ);
        p.gradient(paramsList.get(i));
        double s = p.value;
        double s2 = calc.findLinearised(n, yList.get(i), paramsList.get(i), alpha, beta, func);
        // Exactly the same ...
        Assert.assertEquals(name + " Result: Not same @ " + i, s, s2, 0);
        Assert.assertArrayEquals(name + " Observations: Not same beta @ " + i, p.beta, beta, 0);
        double[] al = p.getAlphaLinear();
        Assert.assertArrayEquals(name + " Observations: Not same alpha @ " + i, al, new DenseMatrix64F(alpha).data, 0);
        double[][] am = p.getAlphaMatrix();
        for (int j = 0; j < nparams; j++) Assert.assertArrayEquals(name + " Observations: Not same alpha @ " + i, am[j], alpha[j], 0);
    }
}
Also used : RandomDataGenerator(org.apache.commons.math3.random.RandomDataGenerator) ArrayList(java.util.ArrayList) Well19937c(org.apache.commons.math3.random.Well19937c) DenseMatrix64F(org.ejml.data.DenseMatrix64F) FakeGradientFunction(gdsc.smlm.function.FakeGradientFunction)

Aggregations

DenseMatrix64F (org.ejml.data.DenseMatrix64F)20 RandomDataGenerator (org.apache.commons.math3.random.RandomDataGenerator)6 Well19937c (org.apache.commons.math3.random.Well19937c)6 ArrayList (java.util.ArrayList)4 FakeGradientFunction (gdsc.smlm.function.FakeGradientFunction)3 Test (org.junit.Test)3 TimingService (gdsc.core.test.TimingService)2 TurboList (gdsc.core.utils.TurboList)2 GradientCalculator (gdsc.smlm.fitting.nonlinear.gradient.GradientCalculator)2 ExtendedGradient2Procedure (gdsc.smlm.function.ExtendedGradient2Procedure)2 ValueProcedure (gdsc.smlm.function.ValueProcedure)2 Gaussian2DFunction (gdsc.smlm.function.gaussian.Gaussian2DFunction)2 Gaussian2DFunctionTest (gdsc.smlm.function.gaussian.Gaussian2DFunctionTest)2 DoubleEquality (gdsc.core.utils.DoubleEquality)1 ChiSquareTest (org.apache.commons.math3.stat.inference.ChiSquareTest)1