Examples with Precision org.apache.commons.math3.util.Precision used on opensource projects

Example 6 with Precision

use of org.apache.commons.math3.util.Precision in project GDSC-SMLM by aherbert.

the class DiffusionRateTest method run.

/*
	 * (non-Javadoc)
	 * 
	 * @see ij.plugin.PlugIn#run(java.lang.String)
	 */
public void run(String arg) {
    SMLMUsageTracker.recordPlugin(this.getClass(), arg);
    if (IJ.controlKeyDown()) {
        simpleTest();
        return;
    }
    extraOptions = Utils.isExtraOptions();
    if (!showDialog())
        return;
    lastSimulatedDataset[0] = lastSimulatedDataset[1] = "";
    lastSimulatedPrecision = 0;
    final int totalSteps = (int) Math.ceil(settings.seconds * settings.stepsPerSecond);
    conversionFactor = 1000000.0 / (settings.pixelPitch * settings.pixelPitch);
    // Diffusion rate is um^2/sec. Convert to pixels per simulation frame.
    final double diffusionRateInPixelsPerSecond = settings.diffusionRate * conversionFactor;
    final double diffusionRateInPixelsPerStep = diffusionRateInPixelsPerSecond / settings.stepsPerSecond;
    final double precisionInPixels = myPrecision / settings.pixelPitch;
    final boolean addError = myPrecision != 0;
    Utils.log(TITLE + " : D = %s um^2/sec, Precision = %s nm", Utils.rounded(settings.diffusionRate, 4), Utils.rounded(myPrecision, 4));
    Utils.log("Mean-displacement per dimension = %s nm/sec", Utils.rounded(1e3 * ImageModel.getRandomMoveDistance(settings.diffusionRate), 4));
    if (extraOptions)
        Utils.log("Step size = %s, precision = %s", Utils.rounded(ImageModel.getRandomMoveDistance(diffusionRateInPixelsPerStep)), Utils.rounded(precisionInPixels));
    // Convert diffusion co-efficient into the standard deviation for the random walk
    final double diffusionSigma = (settings.getDiffusionType() == DiffusionType.LINEAR_WALK) ? // Q. What should this be? At the moment just do 1D diffusion on a random vector
    ImageModel.getRandomMoveDistance(diffusionRateInPixelsPerStep) : ImageModel.getRandomMoveDistance(diffusionRateInPixelsPerStep);
    Utils.log("Simulation step-size = %s nm", Utils.rounded(settings.pixelPitch * diffusionSigma, 4));
    // Move the molecules and get the diffusion rate
    IJ.showStatus("Simulating ...");
    final long start = System.nanoTime();
    final long seed = System.currentTimeMillis() + System.identityHashCode(this);
    RandomGenerator[] random = new RandomGenerator[3];
    RandomGenerator[] random2 = new RandomGenerator[3];
    for (int i = 0; i < 3; i++) {
        random[i] = new Well19937c(seed + i * 12436);
        random2[i] = new Well19937c(seed + i * 678678 + 3);
    }
    Statistics[] stats2D = new Statistics[totalSteps];
    Statistics[] stats3D = new Statistics[totalSteps];
    StoredDataStatistics jumpDistances2D = new StoredDataStatistics(totalSteps);
    StoredDataStatistics jumpDistances3D = new StoredDataStatistics(totalSteps);
    for (int j = 0; j < totalSteps; j++) {
        stats2D[j] = new Statistics();
        stats3D[j] = new Statistics();
    }
    SphericalDistribution dist = new SphericalDistribution(settings.confinementRadius / settings.pixelPitch);
    Statistics asymptote = new Statistics();
    // Save results to memory
    MemoryPeakResults results = new MemoryPeakResults(totalSteps);
    Calibration cal = new Calibration(settings.pixelPitch, 1, 1000.0 / settings.stepsPerSecond);
    results.setCalibration(cal);
    results.setName(TITLE);
    int peak = 0;
    // Store raw coordinates
    ArrayList<Point> points = new ArrayList<Point>(totalSteps);
    StoredData totalJumpDistances1D = new StoredData(settings.particles);
    StoredData totalJumpDistances2D = new StoredData(settings.particles);
    StoredData totalJumpDistances3D = new StoredData(settings.particles);
    for (int i = 0; i < settings.particles; i++) {
        if (i % 16 == 0) {
            IJ.showProgress(i, settings.particles);
            if (Utils.isInterrupted())
                return;
        }
        // Increment the frame so that tracing analysis can distinguish traces
        peak++;
        double[] origin = new double[3];
        final int id = i + 1;
        MoleculeModel m = new MoleculeModel(id, origin.clone());
        if (addError)
            origin = addError(origin, precisionInPixels, random);
        if (useConfinement) {
            // Note: When using confinement the average displacement should asymptote
            // at the average distance of a point from the centre of a ball. This is 3r/4.
            // See: http://answers.yahoo.com/question/index?qid=20090131162630AAMTUfM
            // The equivalent in 2D is 2r/3. However although we are plotting 2D distance
            // this is a projection of the 3D position onto the plane and so the particles
            // will not be evenly spread (there will be clustering at centre caused by the
            // poles)
            final double[] axis = (settings.getDiffusionType() == DiffusionType.LINEAR_WALK) ? nextVector() : null;
            for (int j = 0; j < totalSteps; j++) {
                double[] xyz = m.getCoordinates();
                double[] originalXyz = xyz.clone();
                for (int n = confinementAttempts; n-- > 0; ) {
                    if (settings.getDiffusionType() == DiffusionType.GRID_WALK)
                        m.walk(diffusionSigma, random);
                    else if (settings.getDiffusionType() == DiffusionType.LINEAR_WALK)
                        m.slide(diffusionSigma, axis, random[0]);
                    else
                        m.move(diffusionSigma, random);
                    if (!dist.isWithin(m.getCoordinates())) {
                        // Reset position
                        for (int k = 0; k < 3; k++) xyz[k] = originalXyz[k];
                    } else {
                        // The move was allowed
                        break;
                    }
                }
                points.add(new Point(id, xyz));
                if (addError)
                    xyz = addError(xyz, precisionInPixels, random2);
                peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
            }
            asymptote.add(distance(m.getCoordinates()));
        } else {
            if (settings.getDiffusionType() == DiffusionType.GRID_WALK) {
                for (int j = 0; j < totalSteps; j++) {
                    m.walk(diffusionSigma, random);
                    double[] xyz = m.getCoordinates();
                    points.add(new Point(id, xyz));
                    if (addError)
                        xyz = addError(xyz, precisionInPixels, random2);
                    peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
                }
            } else if (settings.getDiffusionType() == DiffusionType.LINEAR_WALK) {
                final double[] axis = nextVector();
                for (int j = 0; j < totalSteps; j++) {
                    m.slide(diffusionSigma, axis, random[0]);
                    double[] xyz = m.getCoordinates();
                    points.add(new Point(id, xyz));
                    if (addError)
                        xyz = addError(xyz, precisionInPixels, random2);
                    peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
                }
            } else {
                for (int j = 0; j < totalSteps; j++) {
                    m.move(diffusionSigma, random);
                    double[] xyz = m.getCoordinates();
                    points.add(new Point(id, xyz));
                    if (addError)
                        xyz = addError(xyz, precisionInPixels, random2);
                    peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
                }
            }
        }
        // Debug: record all the particles so they can be analysed
        // System.out.printf("%f %f %f\n", m.getX(), m.getY(), m.getZ());
        final double[] xyz = m.getCoordinates();
        double d2 = 0;
        totalJumpDistances1D.add(d2 = xyz[0] * xyz[0]);
        totalJumpDistances2D.add(d2 += xyz[1] * xyz[1]);
        totalJumpDistances3D.add(d2 += xyz[2] * xyz[2]);
    }
    final double time = (System.nanoTime() - start) / 1000000.0;
    IJ.showProgress(1);
    MemoryPeakResults.addResults(results);
    lastSimulatedDataset[0] = results.getName();
    lastSimulatedPrecision = myPrecision;
    // Convert pixels^2/step to um^2/sec
    final double msd2D = (jumpDistances2D.getMean() / conversionFactor) / (results.getCalibration().getExposureTime() / 1000);
    final double msd3D = (jumpDistances3D.getMean() / conversionFactor) / (results.getCalibration().getExposureTime() / 1000);
    Utils.log("Raw data D=%s um^2/s, Precision = %s nm, N=%d, step=%s s, mean2D=%s um^2, MSD 2D = %s um^2/s, mean3D=%s um^2, MSD 3D = %s um^2/s", Utils.rounded(settings.diffusionRate), Utils.rounded(myPrecision), jumpDistances2D.getN(), Utils.rounded(results.getCalibration().getExposureTime() / 1000), Utils.rounded(jumpDistances2D.getMean() / conversionFactor), Utils.rounded(msd2D), Utils.rounded(jumpDistances3D.getMean() / conversionFactor), Utils.rounded(msd3D));
    aggregateIntoFrames(points, addError, precisionInPixels, random2);
    IJ.showStatus("Analysing results ...");
    if (showDiffusionExample) {
        showExample(totalSteps, diffusionSigma, random);
    }
    // Plot a graph of mean squared distance
    double[] xValues = new double[stats2D.length];
    double[] yValues2D = new double[stats2D.length];
    double[] yValues3D = new double[stats3D.length];
    double[] upper2D = new double[stats2D.length];
    double[] lower2D = new double[stats2D.length];
    double[] upper3D = new double[stats3D.length];
    double[] lower3D = new double[stats3D.length];
    SimpleRegression r2D = new SimpleRegression(false);
    SimpleRegression r3D = new SimpleRegression(false);
    final int firstN = (useConfinement) ? fitN : totalSteps;
    for (int j = 0; j < totalSteps; j++) {
        // Convert steps to seconds
        xValues[j] = (double) (j + 1) / settings.stepsPerSecond;
        // Convert values in pixels^2 to um^2
        final double mean2D = stats2D[j].getMean() / conversionFactor;
        final double mean3D = stats3D[j].getMean() / conversionFactor;
        final double sd2D = stats2D[j].getStandardDeviation() / conversionFactor;
        final double sd3D = stats3D[j].getStandardDeviation() / conversionFactor;
        yValues2D[j] = mean2D;
        yValues3D[j] = mean3D;
        upper2D[j] = mean2D + sd2D;
        lower2D[j] = mean2D - sd2D;
        upper3D[j] = mean3D + sd3D;
        lower3D[j] = mean3D - sd3D;
        if (j < firstN) {
            r2D.addData(xValues[j], yValues2D[j]);
            r3D.addData(xValues[j], yValues3D[j]);
        }
    }
    // TODO - Fit using the equation for 2D confined diffusion:
    // MSD = 4s^2 + R^2 (1 - 0.99e^(-1.84^2 Dt / R^2)
    // s = localisation precision
    // R = confinement radius
    // D = 2D diffusion coefficient
    // t = time
    final PolynomialFunction fitted2D, fitted3D;
    if (r2D.getN() > 0) {
        // Do linear regression to get diffusion rate
        final double[] best2D = new double[] { r2D.getIntercept(), r2D.getSlope() };
        fitted2D = new PolynomialFunction(best2D);
        final double[] best3D = new double[] { r3D.getIntercept(), r3D.getSlope() };
        fitted3D = new PolynomialFunction(best3D);
        // For 2D diffusion: d^2 = 4D
        // where: d^2 = mean-square displacement
        double D = best2D[1] / 4.0;
        String msg = "2D Diffusion rate = " + Utils.rounded(D, 4) + " um^2 / sec (" + Utils.timeToString(time) + ")";
        IJ.showStatus(msg);
        Utils.log(msg);
        D = best3D[1] / 6.0;
        Utils.log("3D Diffusion rate = " + Utils.rounded(D, 4) + " um^2 / sec (" + Utils.timeToString(time) + ")");
    } else {
        fitted2D = fitted3D = null;
    }
    // Create plots
    plotMSD(totalSteps, xValues, yValues2D, lower2D, upper2D, fitted2D, 2);
    plotMSD(totalSteps, xValues, yValues3D, lower3D, upper3D, fitted3D, 3);
    plotJumpDistances(TITLE, jumpDistances2D, 2, 1);
    plotJumpDistances(TITLE, jumpDistances3D, 3, 1);
    if (idCount > 0)
        new WindowOrganiser().tileWindows(idList);
    if (useConfinement)
        Utils.log("3D asymptote distance = %s nm (expected %.2f)", Utils.rounded(asymptote.getMean() * settings.pixelPitch, 4), 3 * settings.confinementRadius / 4);
}

Example 7 with Precision

use of org.apache.commons.math3.util.Precision in project GDSC-SMLM by aherbert.

the class TraceDiffusion method fitMSD.

/**
	 * Fit the MSD using a linear fit that must pass through 0,0.
	 * <p>
	 * Update the plot by adding the fit line.
	 * 
	 * @param x
	 * @param y
	 * @param title
	 * @param plot
	 * @return [D, precision]
	 */
private double[] fitMSD(double[] x, double[] y, String title, Plot2 plot) {
    // The Weimann paper (Plos One e64287) fits:
    // MSD(n dt) = 4D n dt + 4s^2
    // n = number of jumps
    // dt = time difference between frames
    // s = localisation precision
    // Thus we should fit an intercept as well.
    // From the fit D = gradient / (4*exposureTime)
    double D = 0;
    double intercept = 0;
    double precision = 0;
    LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
    Optimum lvmSolution;
    double ic = 0;
    // Fit with no intercept
    try {
        final LinearFunction function = new LinearFunction(x, y, settings.fitLength);
        double[] parameters = new double[] { function.guess() };
        //@formatter:off
        LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(3000).start(parameters).target(function.getY()).weight(new DiagonalMatrix(function.getWeights())).model(function, new MultivariateMatrixFunction() {

            public double[][] value(double[] point) throws IllegalArgumentException {
                return function.jacobian(point);
            }
        }).build();
        //@formatter:on
        lvmSolution = optimizer.optimize(problem);
        double ss = lvmSolution.getResiduals().dotProduct(lvmSolution.getResiduals());
        //double ss = 0;
        //double[] obs = function.getY();
        //double[] exp = lvmSolution.getValue();
        //for (int i = 0; i < obs.length; i++)
        //	ss += (obs[i] - exp[i]) * (obs[i] - exp[i]);
        ic = Maths.getAkaikeInformationCriterionFromResiduals(ss, function.getY().length, 1);
        double gradient = lvmSolution.getPoint().getEntry(0);
        D = gradient / 4;
        Utils.log("Linear fit (%d points) : Gradient = %s, D = %s um^2/s, SS = %s, IC = %s (%d evaluations)", function.getY().length, Utils.rounded(gradient, 4), Utils.rounded(D, 4), Utils.rounded(ss), Utils.rounded(ic), lvmSolution.getEvaluations());
    } catch (TooManyIterationsException e) {
        Utils.log("Failed to fit : Too many iterations (%s)", e.getMessage());
    } catch (ConvergenceException e) {
        Utils.log("Failed to fit : %s", e.getMessage());
    }
    // Fit with intercept.
    // Optionally include the intercept (which is the estimated precision).
    boolean fitIntercept = true;
    try {
        final LinearFunctionWithIntercept function = new LinearFunctionWithIntercept(x, y, settings.fitLength, fitIntercept);
        //@formatter:off
        LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(3000).start(function.guess()).target(function.getY()).weight(new DiagonalMatrix(function.getWeights())).model(function, new MultivariateMatrixFunction() {

            public double[][] value(double[] point) throws IllegalArgumentException {
                return function.jacobian(point);
            }
        }).build();
        //@formatter:on
        lvmSolution = optimizer.optimize(problem);
        double ss = lvmSolution.getResiduals().dotProduct(lvmSolution.getResiduals());
        //double ss = 0;
        //double[] obs = function.getY();
        //double[] exp = lvmSolution.getValue();
        //for (int i = 0; i < obs.length; i++)
        //	ss += (obs[i] - exp[i]) * (obs[i] - exp[i]);
        double ic2 = Maths.getAkaikeInformationCriterionFromResiduals(ss, function.getY().length, 2);
        double gradient = lvmSolution.getPoint().getEntry(0);
        final double s = lvmSolution.getPoint().getEntry(1);
        double intercept2 = 4 * s * s;
        if (ic2 < ic || debugFitting) {
            // Convert fitted precision in um to nm
            Utils.log("Linear fit with intercept (%d points) : Gradient = %s, Intercept = %s, D = %s um^2/s, precision = %s nm, SS = %s, IC = %s (%d evaluations)", function.getY().length, Utils.rounded(gradient, 4), Utils.rounded(intercept2, 4), Utils.rounded(gradient / 4, 4), Utils.rounded(s * 1000, 4), Utils.rounded(ss), Utils.rounded(ic2), lvmSolution.getEvaluations());
        }
        if (lvmSolution == null || ic2 < ic) {
            intercept = intercept2;
            D = gradient / 4;
            precision = s;
        }
    } catch (TooManyIterationsException e) {
        Utils.log("Failed to fit with intercept : Too many iterations (%s)", e.getMessage());
    } catch (ConvergenceException e) {
        Utils.log("Failed to fit with intercept : %s", e.getMessage());
    }
    if (settings.msdCorrection) {
        // i.e. the intercept is allowed to be a small negative.
        try {
            // This function fits the jump distance (n) not the time (nt) so update x
            double[] x2 = new double[x.length];
            for (int i = 0; i < x2.length; i++) x2[i] = x[i] / exposureTime;
            final LinearFunctionWithMSDCorrectedIntercept function = new LinearFunctionWithMSDCorrectedIntercept(x2, y, settings.fitLength, fitIntercept);
            //@formatter:off
            LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(3000).start(function.guess()).target(function.getY()).weight(new DiagonalMatrix(function.getWeights())).model(function, new MultivariateMatrixFunction() {

                public double[][] value(double[] point) throws IllegalArgumentException {
                    return function.jacobian(point);
                }
            }).build();
            //@formatter:on
            lvmSolution = optimizer.optimize(problem);
            double ss = lvmSolution.getResiduals().dotProduct(lvmSolution.getResiduals());
            //double ss = 0;
            //double[] obs = function.getY();
            //double[] exp = lvmSolution.getValue();
            //for (int i = 0; i < obs.length; i++)
            //	ss += (obs[i] - exp[i]) * (obs[i] - exp[i]);
            double ic2 = Maths.getAkaikeInformationCriterionFromResiduals(ss, function.getY().length, 2);
            double gradient = lvmSolution.getPoint().getEntry(0);
            final double s = lvmSolution.getPoint().getEntry(1);
            double intercept2 = 4 * s * s - gradient / 3;
            // Q. Is this working?
            // Try fixed precision fitting. Is the gradient correct?
            // Revisit all the equations to see if they are wrong.
            // Try adding the x[0] datapoint using the precision.
            // Change the formula to not be linear at x[0] and to just fit the precision, i.e. the intercept2 = 4 * s * s - gradient / 3 is wrong as the 
            // equation is not linear below n=1.
            // Incorporate the exposure time into the gradient to allow comparison to other fits 
            gradient /= exposureTime;
            if (ic2 < ic || debugFitting) {
                // Convert fitted precision in um to nm
                Utils.log("Linear fit with MSD corrected intercept (%d points) : Gradient = %s, Intercept = %s, D = %s um^2/s, precision = %s nm, SS = %s, IC = %s (%d evaluations)", function.getY().length, Utils.rounded(gradient, 4), Utils.rounded(intercept2, 4), Utils.rounded(gradient / 4, 4), Utils.rounded(s * 1000, 4), Utils.rounded(ss), Utils.rounded(ic2), lvmSolution.getEvaluations());
            }
            if (lvmSolution == null || ic2 < ic) {
                intercept = intercept2;
                D = gradient / 4;
                precision = s;
            }
        } catch (TooManyIterationsException e) {
            Utils.log("Failed to fit with intercept : Too many iterations (%s)", e.getMessage());
        } catch (ConvergenceException e) {
            Utils.log("Failed to fit with intercept : %s", e.getMessage());
        }
    }
    // Add the fit to the plot
    if (D > 0) {
        plot.setColor(Color.magenta);
        plot.drawLine(0, intercept, x[x.length - 1], 4 * D * x[x.length - 1] + intercept);
        display(title, plot);
        checkTraceDistance(D);
    }
    return new double[] { D, precision };
}

Example 8 with Precision

use of org.apache.commons.math3.util.Precision in project GDSC-SMLM by aherbert.

the class BaseFunctionSolverTest method canFitSingleGaussianBetter.

void canFitSingleGaussianBetter(FunctionSolver solver, boolean applyBounds, FunctionSolver solver2, boolean applyBounds2, String name, String name2, NoiseModel noiseModel) {
    double[] noise = getNoise(noiseModel);
    if (solver.isWeighted())
        solver.setWeights(getWeights(noiseModel));
    int LOOPS = 5;
    randomGenerator.setSeed(seed);
    StoredDataStatistics[] stats = new StoredDataStatistics[6];
    String[] statName = { "Signal", "X", "Y" };
    int[] betterPrecision = new int[3];
    int[] totalPrecision = new int[3];
    int[] betterAccuracy = new int[3];
    int[] totalAccuracy = new int[3];
    int i1 = 0, i2 = 0;
    for (double s : signal) {
        double[] expected = createParams(1, s, 0, 0, 1);
        double[] lower = null, upper = null;
        if (applyBounds || applyBounds2) {
            lower = createParams(0, s * 0.5, -0.2, -0.2, 0.8);
            upper = createParams(3, s * 2, 0.2, 0.2, 1.2);
        }
        if (applyBounds)
            solver.setBounds(lower, upper);
        if (applyBounds2)
            solver2.setBounds(lower, upper);
        for (int loop = LOOPS; loop-- > 0; ) {
            double[] data = drawGaussian(expected, noise, noiseModel);
            for (int i = 0; i < stats.length; i++) stats[i] = new StoredDataStatistics();
            for (double db : base) for (double dx : shift) for (double dy : shift) for (double dsx : factor) {
                double[] p = createParams(db, s, dx, dy, dsx);
                double[] fp = fitGaussian(solver, data, p, expected);
                i1 += solver.getEvaluations();
                double[] fp2 = fitGaussian(solver2, data, p, expected);
                i2 += solver2.getEvaluations();
                // Get the mean and sd (the fit precision)
                compare(fp, expected, fp2, expected, Gaussian2DFunction.SIGNAL, stats[0], stats[1]);
                compare(fp, expected, fp2, expected, Gaussian2DFunction.X_POSITION, stats[2], stats[3]);
                compare(fp, expected, fp2, expected, Gaussian2DFunction.Y_POSITION, stats[4], stats[5]);
            // Use the distance
            //stats[2].add(distance(fp, expected));
            //stats[3].add(distance(fp2, expected2));
            }
            // two sided
            double alpha = 0.05;
            for (int i = 0; i < stats.length; i += 2) {
                double u1 = stats[i].getMean();
                double u2 = stats[i + 1].getMean();
                double sd1 = stats[i].getStandardDeviation();
                double sd2 = stats[i + 1].getStandardDeviation();
                TTest tt = new TTest();
                boolean diff = tt.tTest(stats[i].getValues(), stats[i + 1].getValues(), alpha);
                int index = i / 2;
                String msg = String.format("%s vs %s : %.1f (%s) %s %f +/- %f vs %f +/- %f  (N=%d) %b", name2, name, s, noiseModel, statName[index], u2, sd2, u1, sd1, stats[i].getN(), diff);
                if (diff) {
                    // Different means. Check they are roughly the same
                    if (DoubleEquality.almostEqualRelativeOrAbsolute(u1, u2, 0.1, 0)) {
                        // Basically the same. Check which is more precise
                        if (!DoubleEquality.almostEqualRelativeOrAbsolute(sd1, sd2, 0.05, 0)) {
                            if (sd2 < sd1) {
                                betterPrecision[index]++;
                                println(msg + " P*");
                            } else
                                println(msg + " P");
                            totalPrecision[index]++;
                        }
                    } else {
                        // Check which is more accurate (closer to zero)
                        u1 = Math.abs(u1);
                        u2 = Math.abs(u2);
                        if (u2 < u1) {
                            betterAccuracy[index]++;
                            println(msg + " A*");
                        } else
                            println(msg + " A");
                        totalAccuracy[index]++;
                    }
                } else {
                    // The same means. Check that it is more precise
                    if (!DoubleEquality.almostEqualRelativeOrAbsolute(sd1, sd2, 0.05, 0)) {
                        if (sd2 < sd1) {
                            betterPrecision[index]++;
                            println(msg + " P*");
                        } else
                            println(msg + " P");
                        totalPrecision[index]++;
                    }
                }
            }
        }
    }
    int better = 0, total = 0;
    for (int index = 0; index < statName.length; index++) {
        better += betterPrecision[index] + betterAccuracy[index];
        total += totalPrecision[index] + totalAccuracy[index];
        test(name2, name, statName[index] + " P", betterPrecision[index], totalPrecision[index], printBetterDetails);
        test(name2, name, statName[index] + " A", betterAccuracy[index], totalAccuracy[index], printBetterDetails);
    }
    test(name2, name, String.format("All (eval [%d] [%d]) : ", i2, i1), better, total, true);
}

Example 9 with Precision

use of org.apache.commons.math3.util.Precision in project GDSC-SMLM by aherbert.

the class PulseActivationAnalysis method simulateActivations.

private void simulateActivations(RandomDataGenerator rdg, float[][][] molecules, int c, MemoryPeakResults[] results) {
    int n = molecules[c].length;
    if (n == 0)
        return;
    // Compute desired number per frame
    double umPerPixel = sim_nmPerPixel / 1000;
    double um2PerPixel = umPerPixel * umPerPixel;
    double area = sim_size * sim_size * um2PerPixel;
    double nPerFrame = area * sim_activationDensity;
    // Compute the activation probability (but set an upper limit so not all are on in every frame)
    double p = Math.min(0.5, nPerFrame / n);
    // Determine the other channels activation probability using crosstalk
    double[] p0 = { p, p, p };
    int index1, index2, c1, c2;
    switch(c) {
        case 0:
            index1 = C12;
            index2 = C13;
            c1 = 1;
            c2 = 2;
            break;
        case 1:
            index1 = C21;
            index2 = C23;
            c1 = 0;
            c2 = 2;
            break;
        case 2:
        default:
            index1 = C31;
            index2 = C32;
            c1 = 0;
            c2 = 1;
            break;
    }
    p0[c1] *= ct[index1];
    p0[c2] *= ct[index2];
    // Assume 10 frames after each channel pulse => 30 frames per cycle
    double precision = sim_precision[c] / sim_nmPerPixel;
    int id = c + 1;
    RandomGenerator rand = rdg.getRandomGenerator();
    BinomialDistribution[] bd = new BinomialDistribution[4];
    for (int i = 0; i < 3; i++) bd[i] = createBinomialDistribution(rand, n, p0[i]);
    int[] frames = new int[27];
    for (int i = 1, j = 0; i <= 30; i++) {
        if (i % 10 == 1)
            // Skip specific activation frames 
            continue;
        frames[j++] = i;
    }
    bd[3] = createBinomialDistribution(rand, n, p * sim_nonSpecificFrequency);
    // Count the actual cross talk
    int[] count = new int[3];
    for (int i = 0, t = 1; i < sim_cycles; i++, t += 30) {
        count[0] += simulateActivations(rdg, bd[0], molecules[c], results[c], t, precision, id);
        count[1] += simulateActivations(rdg, bd[1], molecules[c], results[c], t + 10, precision, id);
        count[2] += simulateActivations(rdg, bd[2], molecules[c], results[c], t + 20, precision, id);
        // Add non-specific activations
        if (bd[3] != null) {
            for (int t2 : frames) simulateActivations(rdg, bd[3], molecules[c], results[c], t2, precision, id);
        }
    }
    // Report simulated cross talk
    double[] crosstalk = computeCrosstalk(count, c);
    Utils.log("Simulated crosstalk C%s  %s=>%s, C%s  %s=>%s", ctNames[index1], Utils.rounded(ct[index1]), Utils.rounded(crosstalk[c1]), ctNames[index2], Utils.rounded(ct[index2]), Utils.rounded(crosstalk[c2]));
}

Example 10 with Precision

use of org.apache.commons.math3.util.Precision in project GDSC-SMLM by aherbert.

the class PCPALMMolecules method performDistanceAnalysis.

private void performDistanceAnalysis(double[][] intraHist, int p99) {
    // We want to know the fraction of distances between molecules at the 99th percentile
    // that are intra- rather than inter-molecule.
    // Do single linkage clustering of closest pair at this distance and count the number of 
    // links that are inter and intra.
    // Convert molecules for clustering
    ArrayList<ClusterPoint> points = new ArrayList<ClusterPoint>(molecules.size());
    for (Molecule m : molecules) // Precision was used to store the molecule ID
    points.add(ClusterPoint.newClusterPoint((int) m.precision, m.x, m.y, m.photons));
    ClusteringEngine engine = new ClusteringEngine(Prefs.getThreads(), ClusteringAlgorithm.PARTICLE_SINGLE_LINKAGE, new IJTrackProgress());
    IJ.showStatus("Clustering to check inter-molecule distances");
    engine.setTrackJoins(true);
    ArrayList<Cluster> clusters = engine.findClusters(points, intraHist[0][p99]);
    IJ.showStatus("");
    if (clusters != null) {
        double[] intraIdDistances = engine.getIntraIdDistances();
        double[] interIdDistances = engine.getInterIdDistances();
        int all = interIdDistances.length + intraIdDistances.length;
        log("  * Fraction of inter-molecule particle linkage @ %s nm = %s %%", Utils.rounded(intraHist[0][p99], 4), (all > 0) ? Utils.rounded(100.0 * interIdDistances.length / all, 4) : "0");
        // Show a double cumulative histogram plot
        double[][] intraIdHist = Maths.cumulativeHistogram(intraIdDistances, false);
        double[][] interIdHist = Maths.cumulativeHistogram(interIdDistances, false);
        // Plot
        String title = TITLE + " molecule linkage distance";
        Plot2 plot = new Plot2(title, "Distance", "Frequency", intraIdHist[0], intraIdHist[1]);
        double max = (intraIdHist[1].length > 0) ? intraIdHist[1][intraIdHist[1].length - 1] : 0;
        if (interIdHist[1].length > 0)
            max = FastMath.max(max, interIdHist[1][interIdHist[1].length - 1]);
        plot.setLimits(0, intraIdHist[0][intraIdHist[0].length - 1], 0, max);
        plot.setColor(Color.blue);
        plot.addPoints(interIdHist[0], interIdHist[1], Plot2.LINE);
        plot.setColor(Color.black);
        Utils.display(title, plot);
    } else {
        log("Aborted clustering to check inter-molecule distances");
    }
}