Examples with Nullable uk.ac.sussex.gdsc.core.annotation.Nullable used on opensource projects

Example 71 with Nullable

use of uk.ac.sussex.gdsc.core.annotation.Nullable in project gdsc-smlm by aherbert.

the class PcPalmFitting method fitClusteredModel.

/**
 * Fits the correlation curve with r>0 to the clustered model using the estimated density and
 * precision. Parameters must be fit within a tolerance of the starting values.
 *
 * @param gr the correlation curve
 * @param sigmaS The estimated precision
 * @param proteinDensity The estimated protein density
 * @param resultColour the result colour
 * @return The fitted parameters [precision, density, clusterRadius, clusterDensity]
 */
@Nullable
private double[] fitClusteredModel(double[][] gr, double sigmaS, double proteinDensity, String resultColour) {
    final ClusteredModelFunctionGradient function = new ClusteredModelFunctionGradient();
    clusteredModel = function;
    ImageJUtils.log("Fitting %s: Estimated precision = %f nm, estimated protein density = %g um^-2", clusteredModel.getName(), sigmaS, proteinDensity * 1e6);
    clusteredModel.setLogging(true);
    for (int i = offset; i < gr[0].length; i++) {
        // error)
        if (gr[0][i] > sigmaS * settings.fitAboveEstimatedPrecision) {
            clusteredModel.addPoint(gr[0][i], gr[1][i]);
        }
    }
    // The model is: sigma, density, range, amplitude
    final double[] initialSolution = { sigmaS, proteinDensity, sigmaS * 5, 1 };
    // Constrain the fitting to be close to the estimated precision (sigmaS) and protein density.
    // LVM fitting does not support constrained fitting so use a bounded optimiser.
    final SumOfSquaresModelFunction clusteredModelMulti = new SumOfSquaresModelFunction(clusteredModel);
    final double[] x = clusteredModelMulti.x;
    // Put some bounds around the initial guess. Use the fitting tolerance (in %) if provided.
    final double limit = (settings.fittingTolerance > 0) ? 1 + settings.fittingTolerance / 100 : 2;
    final double[] lB = { initialSolution[0] / limit, initialSolution[1] / limit, 0, 0 };
    // The amplitude and range should not extend beyond the limits of the g(r) curve.
    final double[] uB = { initialSolution[0] * limit, initialSolution[1] * limit, MathUtils.max(x), MathUtils.max(gr[1]) };
    ImageJUtils.log("Fitting %s using a bounded search: %s < precision < %s & %s < density < %s", clusteredModel.getName(), MathUtils.rounded(lB[0], 4), MathUtils.rounded(uB[0], 4), MathUtils.rounded(lB[1] * 1e6, 4), MathUtils.rounded(uB[1] * 1e6, 4));
    final PointValuePair constrainedSolution = runBoundedOptimiser(initialSolution, lB, uB, clusteredModelMulti);
    if (constrainedSolution == null) {
        return null;
    }
    double[] parameters = constrainedSolution.getPointRef();
    int evaluations = boundedEvaluations;
    // Refit using a LVM
    if (settings.refitWithGradients) {
        ImageJUtils.log("Re-fitting %s using a gradient optimisation", clusteredModel.getName());
        final LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
        Optimum lvmSolution;
        try {
            // @formatter:off
            final LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(3000).start(parameters).target(function.getY()).weight(new DiagonalMatrix(function.getWeights())).model(function, new MultivariateMatrixFunction() {

                @Override
                public double[][] value(double[] point) {
                    return function.jacobian(point);
                }
            }).build();
            // @formatter:on
            lvmSolution = optimizer.optimize(problem);
            evaluations += lvmSolution.getEvaluations();
            final double ss = lvmSolution.getResiduals().dotProduct(lvmSolution.getResiduals());
            if (ss < constrainedSolution.getValue()) {
                ImageJUtils.log("Re-fitting %s improved the SS from %s to %s (-%s%%)", clusteredModel.getName(), MathUtils.rounded(constrainedSolution.getValue(), 4), MathUtils.rounded(ss, 4), MathUtils.rounded(100 * (constrainedSolution.getValue() - ss) / constrainedSolution.getValue(), 4));
                parameters = lvmSolution.getPoint().toArray();
            }
        } catch (final TooManyIterationsException ex) {
            ImageJUtils.log("Failed to re-fit %s: Too many iterations (%s)", clusteredModel.getName(), ex.getMessage());
        } catch (final ConvergenceException ex) {
            ImageJUtils.log("Failed to re-fit %s: %s", clusteredModel.getName(), ex.getMessage());
        }
    }
    clusteredModel.setLogging(false);
    // Ensure the width is positive
    parameters[0] = Math.abs(parameters[0]);
    double ss = 0;
    final double[] obs = clusteredModel.getY();
    final double[] exp = clusteredModel.value(parameters);
    for (int i = 0; i < obs.length; i++) {
        ss += (obs[i] - exp[i]) * (obs[i] - exp[i]);
    }
    final double totalSumSquares = MathUtils.getTotalSumOfSquares(clusteredModel.getY());
    final double adjustedR2 = MathUtils.getAdjustedCoefficientOfDetermination(ss, totalSumSquares, clusteredModel.size(), parameters.length);
    final double fitSigmaS = parameters[0];
    final double fitProteinDensity = parameters[1];
    // The radius of the cluster domain
    final double domainRadius = parameters[2];
    // The density of the cluster domain
    final double domainDensity = parameters[3];
    // This is from the PC-PALM paper. However that paper fits the g(r)protein exponential convolved
    // in 2D with the g(r)PSF. In this method we have just fit the exponential
    final double nCluster = 2 * domainDensity * Math.PI * domainRadius * domainRadius * fitProteinDensity;
    final double e1 = parameterDrift(sigmaS, fitSigmaS);
    final double e2 = parameterDrift(proteinDensity, fitProteinDensity);
    ImageJUtils.log("  %s fit: SS = %f. Adj.R^2 = %f. %d evaluations", clusteredModel.getName(), ss, adjustedR2, evaluations);
    ImageJUtils.log("  %s parameters:", clusteredModel.getName());
    ImageJUtils.log("    Average precision = %s nm (%s%%)", MathUtils.rounded(fitSigmaS, 4), MathUtils.rounded(e1, 4));
    ImageJUtils.log("    Average protein density = %s um^-2 (%s%%)", MathUtils.rounded(fitProteinDensity * 1e6, 4), MathUtils.rounded(e2, 4));
    ImageJUtils.log("    Domain radius = %s nm", MathUtils.rounded(domainRadius, 4));
    ImageJUtils.log("    Domain density = %s", MathUtils.rounded(domainDensity, 4));
    ImageJUtils.log("    nCluster = %s", MathUtils.rounded(nCluster, 4));
    // Check the fitted parameters are within tolerance of the initial estimates
    valid2 = true;
    if (settings.fittingTolerance > 0 && (Math.abs(e1) > settings.fittingTolerance || Math.abs(e2) > settings.fittingTolerance)) {
        ImageJUtils.log("  Failed to fit %s within tolerance (%s%%): Average precision = %f nm (%s%%)," + " average protein density = %g um^-2 (%s%%)", clusteredModel.getName(), MathUtils.rounded(settings.fittingTolerance, 4), fitSigmaS, MathUtils.rounded(e1, 4), fitProteinDensity * 1e6, MathUtils.rounded(e2, 4));
        valid2 = false;
    }
    // positive
    if (domainRadius < fitSigmaS) {
        ImageJUtils.log("  Failed to fit %s: Domain radius is smaller than the average precision (%s < %s)", clusteredModel.getName(), MathUtils.rounded(domainRadius, 4), MathUtils.rounded(fitSigmaS, 4));
        valid2 = false;
    }
    if (domainDensity < 0) {
        ImageJUtils.log("  Failed to fit %s: Domain density is negative (%s)", clusteredModel.getName(), MathUtils.rounded(domainDensity, 4));
        valid2 = false;
    }
    if (adjustedR2 <= randomModelAdjustedR2) {
        ImageJUtils.log("  Failed to fit %s - Adjusted r^2 has decreased %s%%", clusteredModel.getName(), MathUtils.rounded((100 * (randomModelAdjustedR2 - adjustedR2) / randomModelAdjustedR2), 4));
        valid2 = false;
    }
    addResult(clusteredModel.getName(), resultColour, valid2, fitSigmaS, fitProteinDensity, domainRadius, domainDensity, nCluster, -1, adjustedR2);
    return parameters;
}

Example 72 with Nullable

use of uk.ac.sussex.gdsc.core.annotation.Nullable in project gdsc-smlm by aherbert.

the class PcPalmMolecules method fitGaussian.

@Nullable
private static double[] fitGaussian(float[] x, float[] y) {
    final WeightedObservedPoints obs = new WeightedObservedPoints();
    for (int i = 0; i < x.length; i++) {
        obs.add(x[i], y[i]);
    }
    final Collection<WeightedObservedPoint> observations = obs.toList();
    final GaussianCurveFitter fitter = GaussianCurveFitter.create().withMaxIterations(2000);
    final GaussianCurveFitter.ParameterGuesser guess = new GaussianCurveFitter.ParameterGuesser(observations);
    double[] initialGuess = null;
    try {
        initialGuess = guess.guess();
        return fitter.withStartPoint(initialGuess).fit(observations);
    } catch (final TooManyEvaluationsException ex) {
        // Use the initial estimate
        return initialGuess;
    } catch (final Exception ex) {
        // Just in case there is another exception type, or the initial estimate failed
        return null;
    }
}

Example 73 with Nullable

use of uk.ac.sussex.gdsc.core.annotation.Nullable in project GDSC-SMLM by aherbert.

the class MultiPathFilter method filter.

/**
 * Create a subset of multi-path results, i.e. all those that pass the filter.
 *
 * <p>The results are processed in order. Results are only processed if the fail counter
 * {@link FailCounter#isOk() isOK}.
 *
 * <p>If the subset flag is set to true the candidate Id will be used to determine the number of
 * failed fits before the current candidate, assuming candidates start at zero and increment.
 *
 * @param multiPathResults the multi path results
 * @param failCounter the counter to track the failures to allow per frame before all peaks are
 *        rejected
 * @param setup Set to true to run the {@link #setup()} method
 * @param subset True if a subset (the candidate Id will be used to determine the number of failed
 *        fits before the current candidate)
 * @return the filtered results
 */
@Nullable
private MultiPathFitResult[] filter(final IMultiPathFitResults multiPathResults, final FailCounter failCounter, boolean setup, boolean subset) {
    if (setup) {
        setup();
    }
    failCounter.reset();
    int lastId = -1;
    int size = 0;
    final MultiPathFitResult[] newMultiPathResults = new MultiPathFitResult[multiPathResults.getNumberOfResults()];
    final SimpleSelectedResultStore store = new SimpleSelectedResultStore(multiPathResults.getTotalCandidates());
    // System.out.println("Debug");
    for (int c = 0; c < newMultiPathResults.length; c++) {
        final MultiPathFitResult multiPathResult = multiPathResults.getResult(c);
        // Include the number of failures before this result from the larger set
        if (subset) {
            incrementFailures(failCounter, lastId, multiPathResult);
            lastId = multiPathResult.getCandidateId();
        }
        final boolean evaluateFit = failCounter.isOk();
        if (evaluateFit || store.isValid(multiPathResult.getCandidateId())) {
            // Evaluate the result.
            // This allows storing more estimates in the store even if we are past the failures limit.
            final PreprocessedPeakResult[] result = accept(multiPathResult, false, store);
            // Note: Even if the actual result failed, the candidate may have passed and so
            // the entire multi-path result should be retained.
            // Also note that depending on the filter, different results can be selected and pushed
            // through
            // the store to set them valid. So we must push everything through the store to ensure
            // nothing
            // is removed that could be used.
            checkIsValid(multiPathResult.getSingleFitResult(), store);
            checkIsValid(multiPathResult.getMultiFitResult(), store);
            setupFilter(FilterValidationOption.NO_SHIFT);
            checkIsValid(multiPathResult.getDoubletFitResult(), store);
            // Fix to only disable shift filtering for the doublet results...
            final FitResult multiDoubletFitResult = multiPathResult.getMultiDoubletFitResult();
            if (multiDoubletFitResult != null && multiDoubletFitResult.getResults() != null) {
                // Note: Only disable shift for the doublet results.
                // doublets = len(multi-doublet) - len(multi) + 1
                final PreprocessedPeakResult[] results = multiDoubletFitResult.getResults();
                final int nDoublets = results.length - multiPathResult.getMultiFitResult().getResults().length + 1;
                checkIsValid(results, store, 0, nDoublets);
                restoreFilterState();
                checkIsValid(results, store, nDoublets, results.length);
            } else {
                restoreFilterState();
            }
            // This has valid results so add to the output subset
            newMultiPathResults[size++] = multiPathResult;
            if (evaluateFit) {
                if (isNewResult(result)) {
                    // More results were accepted so reset the fail count
                    failCounter.pass();
                } else {
                    // Nothing was accepted, increment fail count
                    failCounter.fail();
                }
            }
        } else {
            // This was rejected, increment fail count
            failCounter.fail();
        }
        multiPathResults.complete(c);
    }
    if (size != 0) {
        return Arrays.copyOf(newMultiPathResults, size);
    }
    return null;
}

Example 74 with Nullable

use of uk.ac.sussex.gdsc.core.annotation.Nullable in project GDSC-SMLM by aherbert.

the class PeakResultsReader method createNStormResult.

/**
 * Creates the NStorm result.
 *
 * @param line the line
 * @param readPhotons Set to {@code true} if there is a Photons field
 * @return the peak result
 */
@Nullable
private static PeakResult createNStormResult(String line, boolean readPhotons) {
    // Ywc
    try (Scanner scanner = new Scanner(line)) {
        scanner.useDelimiter(tabPattern);
        scanner.useLocale(Locale.US);
        // channelName
        scanner.next();
        // x
        scanner.nextFloat();
        // y
        scanner.nextFloat();
        final float xc = scanner.nextFloat();
        final float yc = scanner.nextFloat();
        final float height = scanner.nextFloat();
        final float area = scanner.nextFloat();
        final float width = scanner.nextFloat();
        // phi
        scanner.nextFloat();
        final float ax = scanner.nextFloat();
        final float bg = scanner.nextFloat();
        // intensity
        scanner.nextFloat();
        final int frame = scanner.nextInt();
        final int length = scanner.nextInt();
        double photons = 0;
        if (readPhotons) {
            // These are not needed
            // Link
            scanner.next();
            // Valid
            scanner.next();
            // Z
            scanner.next();
            // Zc
            scanner.next();
            photons = scanner.nextDouble();
        }
        // The coordinates are in nm
        // The values are in ADUs. The area value is the signal.
        // The following relationship holds when length == 1:
        // Area = Height * 2 * pi * (Width / (pixel_pitch*2) )^2
        // => Pixel_pitch = 0.5 * Width / sqrt(Area / (Height * 2 * pi))
        final float[] params = new float[SIZE_TWO_AXIS];
        params[PeakResult.BACKGROUND] = bg;
        params[PeakResult.INTENSITY] = area;
        params[PeakResult.X] = xc;
        params[PeakResult.Y] = yc;
        // Convert width (2*SD) to SD
        final float sd = width / 2f;
        // Convert to separate XY widths using the axial ratio
        if (ax == 1) {
            params[INDEX_SX] = sd;
            params[INDEX_SY] = sd;
        } else {
            // Ensure the axial ratio is long/short
            final double a = Math.sqrt((ax < 1) ? 1.0 / ax : ax);
            params[INDEX_SX] = (float) (sd * a);
            params[INDEX_SY] = (float) (sd / a);
        }
        // Store the photons in the error value
        return new ExtendedPeakResult(frame, (int) xc, (int) yc, height, photons, 0.0f, 0, params, null, frame + length - 1, 0);
    } catch (final NoSuchElementException ignored) {
    // Ignore
    }
    return null;
}

Example 75 with Nullable

use of uk.ac.sussex.gdsc.core.annotation.Nullable in project GDSC-SMLM by aherbert.

the class BenchmarkFilterAnalysis method findOptimum.

@Nullable
@Override
public SearchResult<FilterScore> findOptimum(double[][] points) {
    gaIteration++;
    SimpleFilterScore max = filterScoreOptimum;
    final FilterScoreResult[] scoreResults = scoreFilters(setStrength(new FilterSet(searchSpaceToFilters(points))), false);
    if (scoreResults == null) {
        return null;
    }
    for (final FilterScoreResult scoreResult : scoreResults) {
        final SimpleFilterScore result = new SimpleFilterScore(scoreResult, true, scoreResult.criteria >= minCriteria);
        if (result.compareTo(max) < 0) {
            max = result;
        }
    }
    filterScoreOptimum = max;
    // Add the best filter to the table
    // This filter may not have been part of the scored subset so use the entire results set for
    // reporting
    final DirectFilter filter = max.result.filter;
    final FractionClassificationResult r = scoreFilter(filter, DEFAULT_MINIMUM_FILTER, gaResultsList, coordinateStore);
    final StringBuilder text = createResult(filter, r);
    add(text, gaIteration);
    gaWindow.accept(text.toString());
    return new SearchResult<>(filter.getParameters(), max);
}