Examples with DiagonalMatrix org.apache.commons.math3.linear.DiagonalMatrix used on opensource projects

Example 21 with DiagonalMatrix

use of org.apache.commons.math3.linear.DiagonalMatrix in project GDSC-SMLM by aherbert.

the class BlinkEstimator method fit.

/**
 * Fit the dark time to counts of molecules curve. Only use the first n fitted points.
 *
 * <p>Calculates:<br> N = The number of photoblinking molecules in the sample<br> nBlink = The
 * average number of blinks per flourophore<br> tOff = The off-time
 *
 * @param td The dark time
 * @param ntd The counts of molecules
 * @param numberOfFittedPointsSetting the number of fitted points
 * @param log Write the fitting results to the ImageJ log window
 * @return The fitted parameters [N, nBlink, tOff], or null if no fit was possible
 */
@Nullable
public double[] fit(double[] td, double[] ntd, int numberOfFittedPointsSetting, boolean log) {
    blinkingModel = new BlinkingFunction();
    blinkingModel.setLogging(true);
    for (int i = 0; i < numberOfFittedPointsSetting; i++) {
        blinkingModel.addPoint(td[i], ntd[i]);
    }
    final LevenbergMarquardtOptimizer optimiser = new LevenbergMarquardtOptimizer(INITIAL_STEP_BOUND_FACTOR, COST_RELATIVE_TOLERANCE, PAR_RELATIVE_TOLERANCE, ORTHO_TOLERANCE, THRESHOLD);
    try {
        final double[] obs = blinkingModel.getY();
        // @formatter:off
        final LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(1000).start(new double[] { ntd[0], 0.1, td[1] }).target(obs).weight(new DiagonalMatrix(blinkingModel.getWeights())).model(blinkingModel, blinkingModel::jacobian).build();
        // @formatter:on
        blinkingModel.setLogging(false);
        final Optimum optimum = optimiser.optimize(problem);
        final double[] parameters = optimum.getPoint().toArray();
        double mean = 0;
        for (final double d : obs) {
            mean += d;
        }
        mean /= obs.length;
        double ssResiduals = 0;
        double ssTotal = 0;
        for (int i = 0; i < obs.length; i++) {
            ssTotal += (obs[i] - mean) * (obs[i] - mean);
        }
        // This is true if the weights are 1
        ssResiduals = optimum.getResiduals().dotProduct(optimum.getResiduals());
        r2 = 1 - ssResiduals / ssTotal;
        adjustedR2 = getAdjustedCoefficientOfDetermination(ssResiduals, ssTotal, obs.length, parameters.length);
        if (log) {
            ImageJUtils.log("  Fit %d points. R^2 = %s. Adjusted R^2 = %s", obs.length, MathUtils.rounded(r2, 4), MathUtils.rounded(adjustedR2, 4));
            ImageJUtils.log("  N=%s, nBlink=%s, tOff=%s (%s frames)", MathUtils.rounded(parameters[0], 4), MathUtils.rounded(parameters[1], 4), MathUtils.rounded(parameters[2], 4), MathUtils.rounded(parameters[2] / msPerFrame, 4));
        }
        return parameters;
    } catch (final TooManyIterationsException ex) {
        if (log) {
            ImageJUtils.log("  Failed to fit %d points: Too many iterations: (%s)", blinkingModel.size(), ex.getMessage());
        }
    } catch (final ConvergenceException ex) {
        if (log) {
            ImageJUtils.log("  Failed to fit %d points", blinkingModel.size());
        }
    }
    return null;
}

Example 22 with DiagonalMatrix

use of org.apache.commons.math3.linear.DiagonalMatrix in project GDSC-SMLM by aherbert.

the class PcPalmFitting method fitRandomModel.

/**
 * Fits the correlation curve with r>0 to the random 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 estimate protein density
 * @param resultColour the result colour
 * @return The fitted parameters [precision, density]
 */
@Nullable
private double[] fitRandomModel(double[][] gr, double sigmaS, double proteinDensity, String resultColour) {
    final RandomModelFunction function = new RandomModelFunction();
    randomModel = function;
    ImageJUtils.log("Fitting %s: Estimated precision = %f nm, estimated protein density = %g um^-2", randomModel.getName(), sigmaS, proteinDensity * 1e6);
    randomModel.setLogging(true);
    for (int i = offset; i < gr[0].length; i++) {
        // error)
        if (gr[0][i] > sigmaS * settings.fitAboveEstimatedPrecision) {
            randomModel.addPoint(gr[0][i], gr[1][i]);
        }
    }
    final LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
    Optimum optimum;
    try {
        // @formatter:off
        final LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(3000).start(new double[] { sigmaS, proteinDensity }).target(function.getY()).weight(new DiagonalMatrix(function.getWeights())).model(function, function::jacobian).build();
        // @formatter:on
        optimum = optimizer.optimize(problem);
    } catch (final TooManyIterationsException ex) {
        ImageJUtils.log("Failed to fit %s: Too many iterations (%s)", randomModel.getName(), ex.getMessage());
        return null;
    } catch (final ConvergenceException ex) {
        ImageJUtils.log("Failed to fit %s: %s", randomModel.getName(), ex.getMessage());
        return null;
    }
    randomModel.setLogging(false);
    final double[] parameters = optimum.getPoint().toArray();
    // Ensure the width is positive
    parameters[0] = Math.abs(parameters[0]);
    final double ss = optimum.getResiduals().dotProduct(optimum.getResiduals());
    final double totalSumSquares = MathUtils.getTotalSumOfSquares(randomModel.getY());
    randomModelAdjustedR2 = MathUtils.getAdjustedCoefficientOfDetermination(ss, totalSumSquares, randomModel.size(), parameters.length);
    final double fitSigmaS = parameters[0];
    final double fitProteinDensity = parameters[1];
    // Check the fitted parameters are within tolerance of the initial estimates
    final double e1 = parameterDrift(sigmaS, fitSigmaS);
    final double e2 = parameterDrift(proteinDensity, fitProteinDensity);
    ImageJUtils.log("  %s fit: SS = %f. Adj.R^2 = %f. %d evaluations", randomModel.getName(), ss, randomModelAdjustedR2, optimum.getEvaluations());
    ImageJUtils.log("  %s parameters:", randomModel.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));
    valid1 = 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%%)", randomModel.getName(), MathUtils.rounded(settings.fittingTolerance, 4), fitSigmaS, MathUtils.rounded(e1, 4), fitProteinDensity * 1e6, MathUtils.rounded(e2, 4));
        valid1 = false;
    }
    if (valid1) {
        // ---------
        // TODO - My data does not comply with this criteria.
        // This could be due to the PC-PALM Molecule code limiting the nmPerPixel to fit the images in
        // memory thus removing correlations at small r.
        // It could also be due to the nature of the random simulations being 3D not 2D membranes
        // as per the PC-PALM paper.
        // ---------
        // Evaluate g(r)protein where:
        // g(r)peaks = g(r)protein + g(r)stoch
        // g(r)peaks ~ 1 + g(r)stoch
        // Verify g(r)protein should be <1.5 for all r>0
        final double[] grStoch = randomModel.value(parameters);
        final double[] grPeaks = randomModel.getY();
        final double[] radius = randomModel.getX();
        for (int i = 0; i < grPeaks.length; i++) {
            // Only evaluate above the fitted average precision
            if (radius[i] < fitSigmaS) {
                continue;
            }
            // Note the RandomModelFunction evaluates g(r)stoch + 1
            final double gr_protein_i = grPeaks[i] - (grStoch[i] - 1);
            if (gr_protein_i > settings.grProteinThreshold) {
                // Failed fit
                ImageJUtils.log("  Failed to fit %s: g(r)protein %s > %s @ r=%s", randomModel.getName(), MathUtils.rounded(gr_protein_i, 4), MathUtils.rounded(settings.grProteinThreshold, 4), MathUtils.rounded(radius[i], 4));
                valid1 = false;
            }
        }
    }
    addResult(randomModel.getName(), resultColour, valid1, fitSigmaS, fitProteinDensity, 0, 0, 0, 0, randomModelAdjustedR2);
    return parameters;
}