Examples with SimpleRegression org.apache.commons.math3.stat.regression.SimpleRegression used on opensource projects

Example 26 with SimpleRegression

use of org.apache.commons.math3.stat.regression.SimpleRegression in project GDSC-SMLM by aherbert.

the class DiffusionRateTest method run.

@Override
public void run(String arg) {
    SmlmUsageTracker.recordPlugin(this.getClass(), arg);
    pluginSettings = Settings.load();
    pluginSettings.save();
    if (IJ.controlKeyDown()) {
        simpleTest();
        return;
    }
    extraOptions = ImageJUtils.isExtraOptions();
    if (!showDialog()) {
        return;
    }
    lastSimulation.set(null);
    final int totalSteps = (int) Math.ceil(settings.getSeconds() * settings.getStepsPerSecond());
    conversionFactor = 1000000.0 / (settings.getPixelPitch() * settings.getPixelPitch());
    // Diffusion rate is um^2/sec. Convert to pixels per simulation frame.
    final double diffusionRateInPixelsPerSecond = settings.getDiffusionRate() * conversionFactor;
    final double diffusionRateInPixelsPerStep = diffusionRateInPixelsPerSecond / settings.getStepsPerSecond();
    final double precisionInPixels = myPrecision / settings.getPixelPitch();
    final boolean addError = myPrecision != 0;
    ImageJUtils.log(TITLE + " : D = %s um^2/sec, Precision = %s nm", MathUtils.rounded(settings.getDiffusionRate(), 4), MathUtils.rounded(myPrecision, 4));
    ImageJUtils.log("Mean-displacement per dimension = %s nm/sec", MathUtils.rounded(1e3 * ImageModel.getRandomMoveDistance(settings.getDiffusionRate()), 4));
    if (extraOptions) {
        ImageJUtils.log("Step size = %s, precision = %s", MathUtils.rounded(ImageModel.getRandomMoveDistance(diffusionRateInPixelsPerStep)), MathUtils.rounded(precisionInPixels));
    }
    // Convert diffusion co-efficient into the standard deviation for the random walk
    final DiffusionType diffusionType = CreateDataSettingsHelper.getDiffusionType(settings.getDiffusionType());
    final double diffusionSigma = ImageModel.getRandomMoveDistance(diffusionRateInPixelsPerStep);
    ImageJUtils.log("Simulation step-size = %s nm", MathUtils.rounded(settings.getPixelPitch() * diffusionSigma, 4));
    // Move the molecules and get the diffusion rate
    IJ.showStatus("Simulating ...");
    final long start = System.nanoTime();
    final UniformRandomProvider random = UniformRandomProviders.create();
    final Statistics[] stats2D = new Statistics[totalSteps];
    final Statistics[] stats3D = new Statistics[totalSteps];
    final StoredDataStatistics jumpDistances2D = new StoredDataStatistics(totalSteps);
    final StoredDataStatistics jumpDistances3D = new StoredDataStatistics(totalSteps);
    for (int j = 0; j < totalSteps; j++) {
        stats2D[j] = new Statistics();
        stats3D[j] = new Statistics();
    }
    final SphericalDistribution dist = new SphericalDistribution(settings.getConfinementRadius() / settings.getPixelPitch());
    final Statistics asymptote = new Statistics();
    // Save results to memory
    final MemoryPeakResults results = new MemoryPeakResults(totalSteps);
    results.setCalibration(CalibrationHelper.create(settings.getPixelPitch(), 1, 1000.0 / settings.getStepsPerSecond()));
    results.setName(TITLE);
    results.setPsf(PsfHelper.create(PSFType.CUSTOM));
    int peak = 0;
    // Store raw coordinates
    final ArrayList<Point> points = new ArrayList<>(totalSteps);
    final StoredData totalJumpDistances1D = new StoredData(settings.getParticles());
    final StoredData totalJumpDistances2D = new StoredData(settings.getParticles());
    final StoredData totalJumpDistances3D = new StoredData(settings.getParticles());
    final NormalizedGaussianSampler gauss = SamplerUtils.createNormalizedGaussianSampler(random);
    for (int i = 0; i < settings.getParticles(); i++) {
        if (i % 16 == 0) {
            IJ.showProgress(i, settings.getParticles());
            if (ImageJUtils.isInterrupted()) {
                return;
            }
        }
        // Increment the frame so that tracing analysis can distinguish traces
        peak++;
        double[] origin = new double[3];
        final int id = i + 1;
        final MoleculeModel m = new MoleculeModel(id, origin.clone());
        if (addError) {
            origin = addError(origin, precisionInPixels, gauss);
        }
        if (pluginSettings.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 = (diffusionType == DiffusionType.LINEAR_WALK) ? nextVector(gauss) : null;
            for (int j = 0; j < totalSteps; j++) {
                double[] xyz = m.getCoordinates();
                final double[] originalXyz = xyz.clone();
                for (int n = pluginSettings.confinementAttempts; n-- > 0; ) {
                    if (diffusionType == DiffusionType.GRID_WALK) {
                        m.walk(diffusionSigma, random);
                    } else if (diffusionType == DiffusionType.LINEAR_WALK) {
                        m.slide(diffusionSigma, axis, random);
                    } 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, gauss);
                }
                peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
            }
            asymptote.add(distance(m.getCoordinates()));
        } else if (diffusionType == 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, gauss);
                }
                peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
            }
        } else if (diffusionType == DiffusionType.LINEAR_WALK) {
            final double[] axis = nextVector(gauss);
            for (int j = 0; j < totalSteps; j++) {
                m.slide(diffusionSigma, axis, random);
                double[] xyz = m.getCoordinates();
                points.add(new Point(id, xyz));
                if (addError) {
                    xyz = addError(xyz, precisionInPixels, gauss);
                }
                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, gauss);
                }
                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 long nanoseconds = System.nanoTime() - start;
    IJ.showProgress(1);
    MemoryPeakResults.addResults(results);
    simulation = new SimulationData(results.getName(), myPrecision);
    // Convert pixels^2/step to um^2/sec
    final double msd2D = (jumpDistances2D.getMean() / conversionFactor) / (results.getCalibrationReader().getExposureTime() / 1000);
    final double msd3D = (jumpDistances3D.getMean() / conversionFactor) / (results.getCalibrationReader().getExposureTime() / 1000);
    ImageJUtils.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", MathUtils.rounded(settings.getDiffusionRate()), MathUtils.rounded(myPrecision), jumpDistances2D.getN(), MathUtils.rounded(results.getCalibrationReader().getExposureTime() / 1000), MathUtils.rounded(jumpDistances2D.getMean() / conversionFactor), MathUtils.rounded(msd2D), MathUtils.rounded(jumpDistances3D.getMean() / conversionFactor), MathUtils.rounded(msd3D));
    aggregateIntoFrames(points, addError, precisionInPixels, gauss);
    IJ.showStatus("Analysing results ...");
    if (pluginSettings.showDiffusionExample) {
        showExample(totalSteps, diffusionSigma, random);
    }
    // Plot a graph of mean squared distance
    final double[] xValues = new double[stats2D.length];
    final double[] yValues2D = new double[stats2D.length];
    final double[] yValues3D = new double[stats3D.length];
    final double[] upper2D = new double[stats2D.length];
    final double[] lower2D = new double[stats2D.length];
    final double[] upper3D = new double[stats3D.length];
    final double[] lower3D = new double[stats3D.length];
    final SimpleRegression r2D = new SimpleRegression(false);
    final SimpleRegression r3D = new SimpleRegression(false);
    final int firstN = (pluginSettings.useConfinement) ? pluginSettings.fitN : totalSteps;
    for (int j = 0; j < totalSteps; j++) {
        // Convert steps to seconds
        xValues[j] = (j + 1) / settings.getStepsPerSecond();
        // 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;
    final PolynomialFunction 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 diffCoeff = best2D[1] / 4.0;
        final String msg = "2D Diffusion rate = " + MathUtils.rounded(diffCoeff, 4) + " um^2 / sec (" + TextUtils.nanosToString(nanoseconds) + ")";
        IJ.showStatus(msg);
        ImageJUtils.log(msg);
        diffCoeff = best3D[1] / 6.0;
        ImageJUtils.log("3D Diffusion rate = " + MathUtils.rounded(diffCoeff, 4) + " um^2 / sec (" + TextUtils.nanosToString(nanoseconds) + ")");
    } 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);
    // Show the total jump length for debugging
    // plotJumpDistances(TITLE + " total", totalJumpDistances1D, 1, totalSteps);
    // plotJumpDistances(TITLE + " total", totalJumpDistances2D, 2, totalSteps);
    // plotJumpDistances(TITLE + " total", totalJumpDistances3D, 3, totalSteps);
    windowOrganiser.tile();
    if (pluginSettings.useConfinement) {
        ImageJUtils.log("3D asymptote distance = %s nm (expected %.2f)", MathUtils.rounded(asymptote.getMean() * settings.getPixelPitch(), 4), 3 * settings.getConfinementRadius() / 4);
    }
}

Example 27 with SimpleRegression

use of org.apache.commons.math3.stat.regression.SimpleRegression in project GDSC-SMLM by aherbert.

the class BenchmarkSpotFit method summariseResults.

private void summariseResults(BenchmarkSpotFitResult spotFitResults, long runTime, final PreprocessedPeakResult[] preprocessedPeakResults, int uniqueIdCount, CandidateData candidateData, TIntObjectHashMap<List<Coordinate>> actualCoordinates) {
    // Summarise the fitting results. N fits, N failures.
    // Optimal match statistics if filtering is perfect (since fitting is not perfect).
    final StoredDataStatistics distanceStats = new StoredDataStatistics();
    final StoredDataStatistics depthStats = new StoredDataStatistics();
    // Get stats for all fitted results and those that match
    // Signal, SNR, Width, xShift, yShift, Precision
    createFilterCriteria();
    final StoredDataStatistics[][] stats = new StoredDataStatistics[3][filterCriteria.length];
    for (int i = 0; i < stats.length; i++) {
        for (int j = 0; j < stats[i].length; j++) {
            stats[i][j] = new StoredDataStatistics();
        }
    }
    final double nmPerPixel = simulationParameters.pixelPitch;
    double tp = 0;
    double fp = 0;
    int failCtp = 0;
    int failCfp = 0;
    int ctp = 0;
    int cfp = 0;
    final int[] singleStatus = new int[FitStatus.values().length];
    final int[] multiStatus = new int[singleStatus.length];
    final int[] doubletStatus = new int[singleStatus.length];
    final int[] multiDoubletStatus = new int[singleStatus.length];
    // Easier to materialise the values since we have a lot of non final variables to manipulate
    final TIntObjectHashMap<FilterCandidates> fitResults = spotFitResults.fitResults;
    final int[] frames = new int[fitResults.size()];
    final FilterCandidates[] candidates = new FilterCandidates[fitResults.size()];
    final int[] counter = new int[1];
    fitResults.forEachEntry((frame, candidate) -> {
        frames[counter[0]] = frame;
        candidates[counter[0]] = candidate;
        counter[0]++;
        return true;
    });
    for (final FilterCandidates result : candidates) {
        // Count the number of fit results that matched (tp) and did not match (fp)
        tp += result.tp;
        fp += result.fp;
        for (int i = 0; i < result.fitResult.length; i++) {
            if (result.spots[i].match) {
                ctp++;
            } else {
                cfp++;
            }
            final MultiPathFitResult fitResult = result.fitResult[i];
            if (singleStatus != null && result.spots[i].match) {
                // Debugging reasons for fit failure
                addStatus(singleStatus, fitResult.getSingleFitResult());
                addStatus(multiStatus, fitResult.getMultiFitResult());
                addStatus(doubletStatus, fitResult.getDoubletFitResult());
                addStatus(multiDoubletStatus, fitResult.getMultiDoubletFitResult());
            }
            if (noMatch(fitResult)) {
                if (result.spots[i].match) {
                    failCtp++;
                } else {
                    failCfp++;
                }
            }
            // We have multi-path results.
            // We want statistics for:
            // [0] all fitted spots
            // [1] fitted spots that match a result
            // [2] fitted spots that do not match a result
            addToStats(fitResult.getSingleFitResult(), stats);
            addToStats(fitResult.getMultiFitResult(), stats);
            addToStats(fitResult.getDoubletFitResult(), stats);
            addToStats(fitResult.getMultiDoubletFitResult(), stats);
        }
        // Statistics on spots that fit an actual result
        for (int i = 0; i < result.match.length; i++) {
            if (!result.match[i].isFitResult()) {
                // For now just ignore the candidates that matched
                continue;
            }
            final FitMatch fitMatch = (FitMatch) result.match[i];
            distanceStats.add(fitMatch.distance * nmPerPixel);
            depthStats.add(fitMatch.zdepth * nmPerPixel);
        }
    }
    if (tp == 0) {
        IJ.error(TITLE, "No fit results matched the simulation actual results");
        return;
    }
    // Store data for computing correlation
    final double[] i1 = new double[depthStats.getN()];
    final double[] i2 = new double[i1.length];
    final double[] is = new double[i1.length];
    int ci = 0;
    for (final FilterCandidates result : candidates) {
        for (int i = 0; i < result.match.length; i++) {
            if (!result.match[i].isFitResult()) {
                // For now just ignore the candidates that matched
                continue;
            }
            final FitMatch fitMatch = (FitMatch) result.match[i];
            final ScoredSpot spot = result.spots[fitMatch.index];
            i1[ci] = fitMatch.predictedSignal;
            i2[ci] = fitMatch.actualSignal;
            is[ci] = spot.spot.intensity;
            ci++;
        }
    }
    // We want to compute the Jaccard against the spot metric
    // Filter the results using the multi-path filter
    final ArrayList<MultiPathFitResults> multiPathResults = new ArrayList<>(fitResults.size());
    for (int i = 0; i < frames.length; i++) {
        final int frame = frames[i];
        final MultiPathFitResult[] multiPathFitResults = candidates[i].fitResult;
        final int totalCandidates = candidates[i].spots.length;
        final List<Coordinate> list = actualCoordinates.get(frame);
        final int nActual = (list == null) ? 0 : list.size();
        multiPathResults.add(new MultiPathFitResults(frame, multiPathFitResults, totalCandidates, nActual));
    }
    // Score the results and count the number returned
    final List<FractionalAssignment[]> assignments = new ArrayList<>();
    final TIntHashSet set = new TIntHashSet(uniqueIdCount);
    final FractionScoreStore scoreStore = set::add;
    final MultiPathFitResults[] multiResults = multiPathResults.toArray(new MultiPathFitResults[0]);
    // Filter with no filter
    final MultiPathFilter mpf = new MultiPathFilter(new SignalFilter(0), null, multiFilter.residualsThreshold);
    mpf.fractionScoreSubset(multiResults, NullFailCounter.INSTANCE, this.results.size(), assignments, scoreStore, CoordinateStoreFactory.create(0, 0, imp.getWidth(), imp.getHeight(), config.convertUsingHwhMax(config.getDuplicateDistanceParameter())));
    final double[][] matchScores = new double[set.size()][];
    int count = 0;
    for (int i = 0; i < assignments.size(); i++) {
        final FractionalAssignment[] a = assignments.get(i);
        if (a == null) {
            continue;
        }
        for (int j = 0; j < a.length; j++) {
            final PreprocessedPeakResult r = ((PeakFractionalAssignment) a[j]).peakResult;
            set.remove(r.getUniqueId());
            final double precision = Math.sqrt(r.getLocationVariance());
            final double signal = r.getSignal();
            final double snr = r.getSnr();
            final double width = r.getXSdFactor();
            final double xShift = r.getXRelativeShift2();
            final double yShift = r.getYRelativeShift2();
            // Since these two are combined for filtering and the max is what matters.
            final double shift = (xShift > yShift) ? Math.sqrt(xShift) : Math.sqrt(yShift);
            final double eshift = Math.sqrt(xShift + yShift);
            final double[] score = new double[8];
            score[FILTER_SIGNAL] = signal;
            score[FILTER_SNR] = snr;
            score[FILTER_MIN_WIDTH] = width;
            score[FILTER_MAX_WIDTH] = width;
            score[FILTER_SHIFT] = shift;
            score[FILTER_ESHIFT] = eshift;
            score[FILTER_PRECISION] = precision;
            score[FILTER_PRECISION + 1] = a[j].getScore();
            matchScores[count++] = score;
        }
    }
    // Add the rest
    set.forEach(new CustomTIntProcedure(count) {

        @Override
        public boolean execute(int uniqueId) {
            // This should not be null or something has gone wrong
            final PreprocessedPeakResult r = preprocessedPeakResults[uniqueId];
            if (r == null) {
                throw new IllegalArgumentException("Missing result: " + uniqueId);
            }
            final double precision = Math.sqrt(r.getLocationVariance());
            final double signal = r.getSignal();
            final double snr = r.getSnr();
            final double width = r.getXSdFactor();
            final double xShift = r.getXRelativeShift2();
            final double yShift = r.getYRelativeShift2();
            // Since these two are combined for filtering and the max is what matters.
            final double shift = (xShift > yShift) ? Math.sqrt(xShift) : Math.sqrt(yShift);
            final double eshift = Math.sqrt(xShift + yShift);
            final double[] score = new double[8];
            score[FILTER_SIGNAL] = signal;
            score[FILTER_SNR] = snr;
            score[FILTER_MIN_WIDTH] = width;
            score[FILTER_MAX_WIDTH] = width;
            score[FILTER_SHIFT] = shift;
            score[FILTER_ESHIFT] = eshift;
            score[FILTER_PRECISION] = precision;
            matchScores[count++] = score;
            return true;
        }
    });
    final FitConfiguration fitConfig = config.getFitConfiguration();
    // Debug the reasons the fit failed
    if (singleStatus != null) {
        String name = PeakFit.getSolverName(fitConfig);
        if (fitConfig.getFitSolver() == FitSolver.MLE && fitConfig.isModelCamera()) {
            name += " Camera";
        }
        IJ.log("Failure counts: " + name);
        printFailures("Single", singleStatus);
        printFailures("Multi", multiStatus);
        printFailures("Doublet", doubletStatus);
        printFailures("Multi doublet", multiDoubletStatus);
    }
    final StringBuilder sb = new StringBuilder(300);
    // Add information about the simulation
    final double signal = simulationParameters.averageSignal;
    final int n = results.size();
    sb.append(imp.getStackSize()).append('\t');
    final int w = imp.getWidth();
    final int h = imp.getHeight();
    sb.append(w).append('\t');
    sb.append(h).append('\t');
    sb.append(n).append('\t');
    final double density = ((double) n / imp.getStackSize()) / (w * h) / (simulationParameters.pixelPitch * simulationParameters.pixelPitch / 1e6);
    sb.append(MathUtils.rounded(density)).append('\t');
    sb.append(MathUtils.rounded(signal)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.sd)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.pixelPitch)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.depth)).append('\t');
    sb.append(simulationParameters.fixedDepth).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.gain)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.readNoise)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.background)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.noise)).append('\t');
    if (simulationParameters.fullSimulation) {
    // The total signal is spread over frames
    }
    sb.append(MathUtils.rounded(signal / simulationParameters.noise)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.sd / simulationParameters.pixelPitch)).append('\t');
    sb.append(spotFilter.getDescription());
    // nP and nN is the fractional score of the spot candidates
    addCount(sb, (double) candidateData.countPositive + candidateData.countNegative);
    addCount(sb, candidateData.countPositive);
    addCount(sb, candidateData.countNegative);
    addCount(sb, candidateData.fractionPositive);
    addCount(sb, candidateData.fractionNegative);
    String name = PeakFit.getSolverName(fitConfig);
    if (fitConfig.getFitSolver() == FitSolver.MLE && fitConfig.isModelCamera()) {
        name += " Camera";
    }
    add(sb, name);
    add(sb, config.getFitting());
    spotFitResults.resultPrefix = sb.toString();
    // Q. Should I add other fit configuration here?
    // The fraction of positive and negative candidates that were included
    add(sb, (100.0 * ctp) / candidateData.countPositive);
    add(sb, (100.0 * cfp) / candidateData.countNegative);
    // Score the fitting results compared to the original simulation.
    // Score the candidate selection:
    add(sb, ctp + cfp);
    add(sb, ctp);
    add(sb, cfp);
    // TP are all candidates that can be matched to a spot
    // FP are all candidates that cannot be matched to a spot
    // FN = The number of missed spots
    FractionClassificationResult match = new FractionClassificationResult(ctp, cfp, 0, simulationParameters.molecules - ctp);
    add(sb, match.getRecall());
    add(sb, match.getPrecision());
    add(sb, match.getF1Score());
    add(sb, match.getJaccard());
    // Score the fitting results:
    add(sb, failCtp);
    add(sb, failCfp);
    // TP are all fit results that can be matched to a spot
    // FP are all fit results that cannot be matched to a spot
    // FN = The number of missed spots
    add(sb, tp);
    add(sb, fp);
    match = new FractionClassificationResult(tp, fp, 0, simulationParameters.molecules - tp);
    add(sb, match.getRecall());
    add(sb, match.getPrecision());
    add(sb, match.getF1Score());
    add(sb, match.getJaccard());
    // Do it again but pretend we can perfectly filter all the false positives
    // add(sb, tp);
    match = new FractionClassificationResult(tp, 0, 0, simulationParameters.molecules - tp);
    // Recall is unchanged
    // Precision will be 100%
    add(sb, match.getF1Score());
    add(sb, match.getJaccard());
    // The mean may be subject to extreme outliers so use the median
    double median = distanceStats.getMedian();
    add(sb, median);
    final WindowOrganiser wo = new WindowOrganiser();
    String label = String.format("Recall = %s. n = %d. Median = %s nm. SD = %s nm", MathUtils.rounded(match.getRecall()), distanceStats.getN(), MathUtils.rounded(median), MathUtils.rounded(distanceStats.getStandardDeviation()));
    new HistogramPlotBuilder(TITLE, distanceStats, "Match Distance (nm)").setPlotLabel(label).show(wo);
    median = depthStats.getMedian();
    add(sb, median);
    // Sort by spot intensity and produce correlation
    double[] correlation = null;
    double[] rankCorrelation = null;
    double[] rank = null;
    final FastCorrelator fastCorrelator = new FastCorrelator();
    final ArrayList<Ranking> pc1 = new ArrayList<>();
    final ArrayList<Ranking> pc2 = new ArrayList<>();
    ci = 0;
    if (settings.showCorrelation) {
        final int[] indices = SimpleArrayUtils.natural(i1.length);
        SortUtils.sortData(indices, is, settings.rankByIntensity, true);
        correlation = new double[i1.length];
        rankCorrelation = new double[i1.length];
        rank = new double[i1.length];
        for (final int ci2 : indices) {
            fastCorrelator.add(Math.round(i1[ci2]), Math.round(i2[ci2]));
            pc1.add(new Ranking(i1[ci2], ci));
            pc2.add(new Ranking(i2[ci2], ci));
            correlation[ci] = fastCorrelator.getCorrelation();
            rankCorrelation[ci] = Correlator.correlation(rank(pc1), rank(pc2));
            if (settings.rankByIntensity) {
                rank[ci] = is[0] - is[ci];
            } else {
                rank[ci] = ci;
            }
            ci++;
        }
    } else {
        for (int i = 0; i < i1.length; i++) {
            fastCorrelator.add(Math.round(i1[i]), Math.round(i2[i]));
            pc1.add(new Ranking(i1[i], i));
            pc2.add(new Ranking(i2[i], i));
        }
    }
    final double pearsonCorr = fastCorrelator.getCorrelation();
    final double rankedCorr = Correlator.correlation(rank(pc1), rank(pc2));
    // Get the regression
    final SimpleRegression regression = new SimpleRegression(false);
    for (int i = 0; i < pc1.size(); i++) {
        regression.addData(pc1.get(i).value, pc2.get(i).value);
    }
    // final double intercept = regression.getIntercept();
    final double slope = regression.getSlope();
    if (settings.showCorrelation) {
        String title = TITLE + " Intensity";
        Plot plot = new Plot(title, "Candidate", "Spot");
        final double[] limits1 = MathUtils.limits(i1);
        final double[] limits2 = MathUtils.limits(i2);
        plot.setLimits(limits1[0], limits1[1], limits2[0], limits2[1]);
        label = String.format("Correlation=%s; Ranked=%s; Slope=%s", MathUtils.rounded(pearsonCorr), MathUtils.rounded(rankedCorr), MathUtils.rounded(slope));
        plot.addLabel(0, 0, label);
        plot.setColor(Color.red);
        plot.addPoints(i1, i2, Plot.DOT);
        if (slope > 1) {
            plot.drawLine(limits1[0], limits1[0] * slope, limits1[1], limits1[1] * slope);
        } else {
            plot.drawLine(limits2[0] / slope, limits2[0], limits2[1] / slope, limits2[1]);
        }
        ImageJUtils.display(title, plot, wo);
        title = TITLE + " Correlation";
        plot = new Plot(title, "Spot Rank", "Correlation");
        final double[] xlimits = MathUtils.limits(rank);
        double[] ylimits = MathUtils.limits(correlation);
        ylimits = MathUtils.limits(ylimits, rankCorrelation);
        plot.setLimits(xlimits[0], xlimits[1], ylimits[0], ylimits[1]);
        plot.setColor(Color.red);
        plot.addPoints(rank, correlation, Plot.LINE);
        plot.setColor(Color.blue);
        plot.addPoints(rank, rankCorrelation, Plot.LINE);
        plot.setColor(Color.black);
        plot.addLabel(0, 0, label);
        ImageJUtils.display(title, plot, wo);
    }
    add(sb, pearsonCorr);
    add(sb, rankedCorr);
    add(sb, slope);
    label = String.format("n = %d. Median = %s nm", depthStats.getN(), MathUtils.rounded(median));
    new HistogramPlotBuilder(TITLE, depthStats, "Match Depth (nm)").setRemoveOutliersOption(1).setPlotLabel(label).show(wo);
    // Plot histograms of the stats on the same window
    final double[] lower = new double[filterCriteria.length];
    final double[] upper = new double[lower.length];
    final double[] min = new double[lower.length];
    final double[] max = new double[lower.length];
    for (int i = 0; i < stats[0].length; i++) {
        final double[] limits = showDoubleHistogram(stats, i, wo, matchScores);
        lower[i] = limits[0];
        upper[i] = limits[1];
        min[i] = limits[2];
        max[i] = limits[3];
    }
    // Reconfigure some of the range limits
    // Make this a bit bigger
    upper[FILTER_SIGNAL] *= 2;
    // Make this a bit bigger
    upper[FILTER_SNR] *= 2;
    final double factor = 0.25;
    if (lower[FILTER_MIN_WIDTH] != 0) {
        // (assuming lower is less than 1)
        upper[FILTER_MIN_WIDTH] = 1 - Math.max(0, factor * (1 - lower[FILTER_MIN_WIDTH]));
    }
    if (upper[FILTER_MIN_WIDTH] != 0) {
        // (assuming upper is more than 1)
        lower[FILTER_MAX_WIDTH] = 1 + Math.max(0, factor * (upper[FILTER_MAX_WIDTH] - 1));
    }
    // Round the ranges
    final double[] interval = new double[stats[0].length];
    interval[FILTER_SIGNAL] = SignalFilter.DEFAULT_INCREMENT;
    interval[FILTER_SNR] = SnrFilter.DEFAULT_INCREMENT;
    interval[FILTER_MIN_WIDTH] = WidthFilter2.DEFAULT_MIN_INCREMENT;
    interval[FILTER_MAX_WIDTH] = WidthFilter.DEFAULT_INCREMENT;
    interval[FILTER_SHIFT] = ShiftFilter.DEFAULT_INCREMENT;
    interval[FILTER_ESHIFT] = EShiftFilter.DEFAULT_INCREMENT;
    interval[FILTER_PRECISION] = PrecisionFilter.DEFAULT_INCREMENT;
    interval[FILTER_ITERATIONS] = 0.1;
    interval[FILTER_EVALUATIONS] = 0.1;
    // Create a range increment
    final double[] increment = new double[lower.length];
    for (int i = 0; i < increment.length; i++) {
        lower[i] = MathUtils.floor(lower[i], interval[i]);
        upper[i] = MathUtils.ceil(upper[i], interval[i]);
        final double range = upper[i] - lower[i];
        // Allow clipping if the range is small compared to the min increment
        double multiples = range / interval[i];
        // Use 8 multiples for the equivalent of +/- 4 steps around the centre
        if (multiples < 8) {
            multiples = Math.ceil(multiples);
        } else {
            multiples = 8;
        }
        increment[i] = MathUtils.ceil(range / multiples, interval[i]);
        if (i == FILTER_MIN_WIDTH) {
            // Requires clipping based on the upper limit
            lower[i] = upper[i] - increment[i] * multiples;
        } else {
            upper[i] = lower[i] + increment[i] * multiples;
        }
    }
    for (int i = 0; i < stats[0].length; i++) {
        lower[i] = MathUtils.round(lower[i]);
        upper[i] = MathUtils.round(upper[i]);
        min[i] = MathUtils.round(min[i]);
        max[i] = MathUtils.round(max[i]);
        increment[i] = MathUtils.round(increment[i]);
        sb.append('\t').append(min[i]).append(':').append(lower[i]).append('-').append(upper[i]).append(':').append(max[i]);
    }
    // Disable some filters
    increment[FILTER_SIGNAL] = Double.POSITIVE_INFINITY;
    // increment[FILTER_SHIFT] = Double.POSITIVE_INFINITY;
    increment[FILTER_ESHIFT] = Double.POSITIVE_INFINITY;
    wo.tile();
    sb.append('\t').append(TextUtils.nanosToString(runTime));
    createTable().append(sb.toString());
    if (settings.saveFilterRange) {
        GUIFilterSettings filterSettings = SettingsManager.readGuiFilterSettings(0);
        String filename = (silent) ? filterSettings.getFilterSetFilename() : ImageJUtils.getFilename("Filter_range_file", filterSettings.getFilterSetFilename());
        if (filename == null) {
            return;
        }
        // Remove extension to store the filename
        filename = FileUtils.replaceExtension(filename, ".xml");
        filterSettings = filterSettings.toBuilder().setFilterSetFilename(filename).build();
        // Create a filter set using the ranges
        final ArrayList<Filter> filters = new ArrayList<>(4);
        // Create the multi-filter using the same precision type as that used during fitting.
        // Currently no support for z-filter as 3D astigmatism fitting is experimental.
        final PrecisionMethod precisionMethod = getPrecisionMethod((DirectFilter) multiFilter.getFilter());
        Function<double[], Filter> generator;
        if (precisionMethod == PrecisionMethod.POISSON_CRLB) {
            generator = parameters -> new MultiFilterCrlb(parameters[FILTER_SIGNAL], (float) parameters[FILTER_SNR], parameters[FILTER_MIN_WIDTH], parameters[FILTER_MAX_WIDTH], parameters[FILTER_SHIFT], parameters[FILTER_ESHIFT], parameters[FILTER_PRECISION], 0f, 0f);
        } else if (precisionMethod == PrecisionMethod.MORTENSEN) {
            generator = parameters -> new MultiFilter(parameters[FILTER_SIGNAL], (float) parameters[FILTER_SNR], parameters[FILTER_MIN_WIDTH], parameters[FILTER_MAX_WIDTH], parameters[FILTER_SHIFT], parameters[FILTER_ESHIFT], parameters[FILTER_PRECISION], 0f, 0f);
        } else {
            // Default
            generator = parameters -> new MultiFilter2(parameters[FILTER_SIGNAL], (float) parameters[FILTER_SNR], parameters[FILTER_MIN_WIDTH], parameters[FILTER_MAX_WIDTH], parameters[FILTER_SHIFT], parameters[FILTER_ESHIFT], parameters[FILTER_PRECISION], 0f, 0f);
        }
        filters.add(generator.apply(lower));
        filters.add(generator.apply(upper));
        filters.add(generator.apply(increment));
        if (saveFilters(filename, filters)) {
            SettingsManager.writeSettings(filterSettings);
        }
        // Create a filter set using the min/max and the initial bounds.
        // Set sensible limits
        min[FILTER_SIGNAL] = Math.max(min[FILTER_SIGNAL], 30);
        max[FILTER_SNR] = Math.min(max[FILTER_SNR], 10000);
        max[FILTER_PRECISION] = Math.min(max[FILTER_PRECISION], 100);
        // Make the 4-set filters the same as the 3-set filters.
        filters.clear();
        filters.add(generator.apply(min));
        filters.add(generator.apply(lower));
        filters.add(generator.apply(upper));
        filters.add(generator.apply(max));
        saveFilters(FileUtils.replaceExtension(filename, ".4.xml"), filters);
    }
    spotFitResults.min = min;
    spotFitResults.max = max;
}

Example 28 with SimpleRegression

use of org.apache.commons.math3.stat.regression.SimpleRegression in project GDSC-SMLM by aherbert.

the class BenchmarkFilterAnalysis method reportResults.

private ComplexFilterScore reportResults(boolean newResults, List<ComplexFilterScore> filters) {
    if (filters.isEmpty()) {
        IJ.log("Warning: No filters pass the criteria");
        return null;
    }
    getCoordinateStore();
    Collections.sort(filters);
    FractionClassificationResult topFilterClassificationResult = null;
    ArrayList<FractionalAssignment[]> topFilterResults = null;
    String topFilterSummary = null;
    if (settings.showSummaryTable || settings.saveTemplate) {
        final Consumer<String> summaryWindow = createSummaryWindow();
        int count = 0;
        final double range = (settings.summaryDepth / simulationParameters.pixelPitch) * 0.5;
        final int[] np = { 0 };
        fitResultData.depthStats.forEach(depth -> {
            if (Math.abs(depth) < range) {
                np[0]++;
            }
        });
        for (final ComplexFilterScore fs : filters) {
            final ArrayList<FractionalAssignment[]> list = new ArrayList<>(fitResultData.resultsList.length);
            final FractionClassificationResult r = scoreFilter(fs.getFilter(), defaultMinimalFilter, fitResultData.resultsList, list, coordinateStore);
            final StringBuilder sb = createResult(fs.getFilter(), r);
            if (topFilterResults == null) {
                topFilterResults = list;
                topFilterClassificationResult = r;
            }
            // Show the recall at the specified depth. Sum the distance and signal factor of all scored
            // spots.
            int scored = 0;
            double tp = 0;
            double distance = 0;
            double sf = 0;
            double rmsd = 0;
            final SimpleRegression regression = new SimpleRegression(false);
            for (final FractionalAssignment[] assignments : list) {
                if (assignments == null) {
                    continue;
                }
                for (int i = 0; i < assignments.length; i++) {
                    final CustomFractionalAssignment c = (CustomFractionalAssignment) assignments[i];
                    if (Math.abs(c.peak.getZPosition()) <= range) {
                        tp += c.getScore();
                    }
                    distance += c.distToTarget;
                    sf += c.getSignalFactor();
                    rmsd += c.distToTarget * c.distToTarget;
                    regression.addData(c.peakResult.getSignal(), c.peak.getIntensity());
                }
                scored += assignments.length;
            }
            final double slope = regression.getSlope();
            sb.append('\t');
            sb.append(MathUtils.rounded(tp / np[0])).append('\t');
            sb.append(MathUtils.rounded(distance / scored)).append('\t');
            sb.append(MathUtils.rounded(sf / scored)).append('\t');
            // RMSD to be the root mean square deviation in a single dimension so divide by 2.
            // (This assumes 2D Euclidean distances.)
            sb.append(MathUtils.rounded(Math.sqrt(MathUtils.div0(rmsd / 2, scored)))).append('\t');
            sb.append(MathUtils.rounded(slope)).append('\t');
            if (fs.atLimit() != null) {
                sb.append(fs.atLimit());
            }
            String text = sb.toString();
            if (topFilterSummary == null) {
                topFilterSummary = text;
                if (!settings.showSummaryTable) {
                    break;
                }
            }
            if (fs.time != 0) {
                sb.append('\t');
                sb.append(fs.algorithm);
                sb.append('\t');
                sb.append(org.apache.commons.lang3.time.DurationFormatUtils.formatDurationHMS(fs.time));
            } else {
                sb.append("\t\t");
            }
            if (fs.paramTime != 0) {
                sb.append('\t');
                sb.append(fs.getParamAlgorithm());
                sb.append('\t');
                sb.append(org.apache.commons.lang3.time.DurationFormatUtils.formatDurationHMS(fs.paramTime));
            } else {
                sb.append("\t\t");
            }
            text = sb.toString();
            summaryWindow.accept(text);
            count++;
            if (settings.summaryTopN > 0 && count >= settings.summaryTopN) {
                break;
            }
        }
        // Add a spacer to the summary table if we have multiple results
        if (count > 1 && settings.showSummaryTable) {
            summaryWindow.accept("");
        }
    }
    final DirectFilter localBestFilter = filters.get(0).getFilter();
    if (settings.saveBestFilter) {
        saveFilter(localBestFilter);
    }
    if (topFilterClassificationResult == null) {
        topFilterResults = new ArrayList<>(fitResultData.resultsList.length);
        scoreFilter(localBestFilter, defaultMinimalFilter, fitResultData.resultsList, topFilterResults, coordinateStore);
    }
    if (newResults || filterAnalysisResult.scores.isEmpty()) {
        filterAnalysisResult.scores.add(new FilterResult(settings.failCount, residualsThreshold, settings.duplicateDistance, settings.duplicateDistanceAbsolute, filters.get(0)));
    }
    if (settings.saveTemplate) {
        saveTemplate(topFilterSummary);
    }
    showPlots();
    calculateSensitivity();
    topFilterResults = depthAnalysis(topFilterResults, localBestFilter);
    topFilterResults = scoreAnalysis(topFilterResults, localBestFilter);
    componentAnalysis(filters.get(0));
    PreprocessedPeakResult[] filterResults = null;
    if (isShowOverlay()) {
        filterResults = showOverlay(topFilterResults, localBestFilter);
    }
    saveResults(filterResults, localBestFilter);
    wo.tile();
    return filters.get(0);
}

Example 29 with SimpleRegression

use of org.apache.commons.math3.stat.regression.SimpleRegression in project GDSC-SMLM by aherbert.

the class Gaussian2DPeakResultHelperTest method canComputePixelAmplitude.

@Test
void canComputePixelAmplitude() {
    final float[] x = new float[] { 0f, 0.1f, 0.3f, 0.5f, 0.7f, 1f };
    final float[] s = new float[] { 0.8f, 1f, 1.5f, 2.2f };
    final float[] paramsf = new float[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
    paramsf[Gaussian2DFunction.BACKGROUND] = 0;
    paramsf[Gaussian2DFunction.SIGNAL] = 105;
    final Gaussian2DFunction f = GaussianFunctionFactory.create2D(1, 1, 1, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    final SimpleRegression r = new SimpleRegression(false);
    for (final float tx : x) {
        for (final float ty : x) {
            for (final float sx : s) {
                for (final float sy : s) {
                    paramsf[Gaussian2DFunction.X_POSITION] = tx;
                    paramsf[Gaussian2DFunction.Y_POSITION] = ty;
                    paramsf[Gaussian2DFunction.X_SD] = sx;
                    paramsf[Gaussian2DFunction.Y_SD] = sy;
                    // Get the answer using a single pixel image
                    // Note the Gaussian2D functions set the centre of the pixel as 0,0 so offset
                    final double[] params = SimpleArrayUtils.toDouble(paramsf);
                    params[Gaussian2DFunction.X_POSITION] -= 0.5;
                    params[Gaussian2DFunction.Y_POSITION] -= 0.5;
                    f.initialise0(params);
                    final double e = f.eval(0);
                    final PSF psf = PsfHelper.create(PSFType.TWO_AXIS_GAUSSIAN_2D);
                    final CalibrationWriter calibration = new CalibrationWriter();
                    calibration.setCountPerPhoton(1);
                    calibration.setIntensityUnit(IntensityUnit.PHOTON);
                    calibration.setNmPerPixel(1);
                    calibration.setDistanceUnit(DistanceUnit.PIXEL);
                    final Gaussian2DPeakResultCalculator calc = Gaussian2DPeakResultHelper.create(psf, calibration, Gaussian2DPeakResultHelper.AMPLITUDE | Gaussian2DPeakResultHelper.PIXEL_AMPLITUDE);
                    final double o1 = calc.getAmplitude(paramsf);
                    final double o2 = calc.getPixelAmplitude(paramsf);
                    // logger.fine(FunctionUtils.getSupplier("e=%f, o1=%f, o2=%f", e, o1, o2));
                    Assertions.assertEquals(e, o2, 1e-3);
                    r.addData(e, o1);
                }
            }
        }
    }
    // logger.fine(FunctionUtils.getSupplier("Regression: pixel amplitude vs amplitude = %f,
    // slope=%f, n=%d", r.getR(), r.getSlope(),
    // r.getN()));
    // The simple amplitude over estimates the actual pixel amplitude
    Assertions.assertTrue(r.getSlope() > 1);
}

Example 30 with SimpleRegression

use of org.apache.commons.math3.stat.regression.SimpleRegression in project ta4j by ta4j.

the class SimpleLinearRegressionIndicatorTest method calculateLinearRegression.

@Test
public void calculateLinearRegression() {
    double[] values = new double[] { 1, 2, 1.3, 3.75, 2.25 };
    ClosePriceIndicator indicator = new ClosePriceIndicator(new MockTimeSeries(values));
    SimpleLinearRegressionIndicator reg = new SimpleLinearRegressionIndicator(indicator, 5);
    SimpleRegression origReg = buildSimpleRegression(values);
    assertDecimalEquals(reg.getValue(4), origReg.predict(4));
}