Examples with TooManyIterationsException org.apache.commons.math3.exception.TooManyIterationsException used on opensource projects

Example 16 with TooManyIterationsException

use of org.apache.commons.math3.exception.TooManyIterationsException 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 the x
 * @param y the y
 * @param title the title
 * @param plot the plot
 * @return [D, precision]
 */
private double[] fitMsd(double[] x, double[] y, String title, Plot 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)
    // D
    double diffCoeff = 0;
    double intercept = 0;
    double precision = 0;
    final LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
    Optimum lvmSolution;
    double ic = Double.NaN;
    // Fit with no intercept
    try {
        final LinearFunction function = new LinearFunction(x, y, clusteringSettings.getFitLength());
        final double[] parameters = new double[] { function.guess() };
        // @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, function::jacobian).build();
        // @formatter:on
        lvmSolution = optimizer.optimize(problem);
        final 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 = getAkaikeInformationCriterionFromResiduals(ss, function.getY().length, 1);
        final double gradient = lvmSolution.getPoint().getEntry(0);
        diffCoeff = gradient / 4;
        ImageJUtils.log("Linear fit (%d points) : Gradient = %s, D = %s um^2/s, SS = %s, " + "IC = %s (%d evaluations)", function.getY().length, MathUtils.rounded(gradient, 4), MathUtils.rounded(diffCoeff, 4), MathUtils.rounded(ss), MathUtils.rounded(ic), lvmSolution.getEvaluations());
    } catch (final TooManyIterationsException ex) {
        ImageJUtils.log("Failed to fit : Too many iterations (%s)", ex.getMessage());
    } catch (final ConvergenceException ex) {
        ImageJUtils.log("Failed to fit : %s", ex.getMessage());
    }
    // Fit with intercept.
    // Optionally include the intercept (which is the estimated precision).
    final boolean fitIntercept = true;
    try {
        final LinearFunctionWithIntercept function = new LinearFunctionWithIntercept(x, y, clusteringSettings.getFitLength(), fitIntercept);
        // @formatter:off
        final LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(3000).start(function.guess()).target(function.getY()).weight(new DiagonalMatrix(function.getWeights())).model(function, function::jacobian).build();
        // @formatter:on
        lvmSolution = optimizer.optimize(problem);
        final RealVector residuals = lvmSolution.getResiduals();
        final double ss = residuals.dotProduct(residuals);
        // 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]);
        final double ic2 = getAkaikeInformationCriterionFromResiduals(ss, function.getY().length, 2);
        final double gradient = lvmSolution.getPoint().getEntry(0);
        final double s = lvmSolution.getPoint().getEntry(1);
        final double intercept2 = 4 * s * s;
        if (ic2 < ic || Double.isNaN(ic)) {
            if (settings.debugFitting) {
                // Convert fitted precision in um to nm
                ImageJUtils.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, MathUtils.rounded(gradient, 4), MathUtils.rounded(intercept2, 4), MathUtils.rounded(gradient / 4, 4), MathUtils.rounded(s * 1000, 4), MathUtils.rounded(ss), MathUtils.rounded(ic2), lvmSolution.getEvaluations());
            }
            intercept = intercept2;
            diffCoeff = gradient / 4;
            precision = s;
        }
    } catch (final TooManyIterationsException ex) {
        ImageJUtils.log("Failed to fit with intercept : Too many iterations (%s)", ex.getMessage());
    } catch (final ConvergenceException ex) {
        ImageJUtils.log("Failed to fit with intercept : %s", ex.getMessage());
    }
    if (clusteringSettings.getMsdCorrection()) {
        // 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
            final 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, clusteringSettings.getFitLength(), fitIntercept);
            // @formatter:off
            final LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(3000).start(function.guess()).target(function.getY()).weight(new DiagonalMatrix(function.getWeights())).model(function, function::jacobian).build();
            // @formatter:on
            lvmSolution = optimizer.optimize(problem);
            final RealVector residuals = lvmSolution.getResiduals();
            final double ss = residuals.dotProduct(residuals);
            // 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]);
            final double ic2 = getAkaikeInformationCriterionFromResiduals(ss, function.getY().length, 2);
            double gradient = lvmSolution.getPoint().getEntry(0);
            final double s = lvmSolution.getPoint().getEntry(1);
            final 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 || Double.isNaN(ic)) {
                if (settings.debugFitting) {
                    // Convert fitted precision in um to nm
                    ImageJUtils.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, MathUtils.rounded(gradient, 4), MathUtils.rounded(intercept2, 4), MathUtils.rounded(gradient / 4, 4), MathUtils.rounded(s * 1000, 4), MathUtils.rounded(ss), MathUtils.rounded(ic2), lvmSolution.getEvaluations());
                }
                intercept = intercept2;
                diffCoeff = gradient / 4;
                precision = s;
            }
        } catch (final TooManyIterationsException ex) {
            ImageJUtils.log("Failed to fit with intercept : Too many iterations (%s)", ex.getMessage());
        } catch (final ConvergenceException ex) {
            ImageJUtils.log("Failed to fit with intercept : %s", ex.getMessage());
        }
    }
    // Add the fit to the plot
    if (diffCoeff > 0) {
        plot.setColor(Color.magenta);
        plot.drawLine(0, intercept, x[x.length - 1], 4 * diffCoeff * x[x.length - 1] + intercept);
        display(title, plot);
        checkTraceSettings(diffCoeff);
    }
    return new double[] { diffCoeff, precision };
}

Example 17 with TooManyIterationsException

use of org.apache.commons.math3.exception.TooManyIterationsException in project GDSC-SMLM by aherbert.

the class TrackPopulationAnalysis method extractTrackData.

/**
 * Extract the track data. This extracts different descriptors of the track using a rolling local
 * window.
 *
 * <p>Distances are converted to {@code unit} using the provided converter and time units are
 * converted from frame to seconds (s). The diffusion coefficients is in unit^2/s.
 *
 * <p>If categories are added they are remapped to be a natural sequence starting from 0.
 *
 * @param tracks the tracks
 * @param distanceConverter the distance converter
 * @param deltaT the time step of each frame in seconds (delta T)
 * @param hasCategory if true add the category from the result to the end of the data
 * @return the track data (track lengths and descriptors)
 */
private Pair<int[], double[][]> extractTrackData(List<Trace> tracks, TypeConverter<DistanceUnit> distanceConverter, double deltaT, boolean hasCategory) {
    final List<double[]> data = new LocalList<>(tracks.size());
    double[] x = new double[0];
    double[] y = new double[0];
    final int wm1 = settings.window - 1;
    final int valuesLength = hasCategory ? 5 : 4;
    final int mid = settings.window / 2;
    final MsdFitter msdFitter = new MsdFitter(settings, deltaT);
    final double significance = settings.significance;
    final int[] fitResult = new int[4];
    // Factor for the diffusion coefficient: 1/N * 1/(2*dimensions*deltaT) = 1 / 4Nt
    // with N the number of points to average.
    final double diffusionCoefficientFactor = 1.0 / (4 * wm1 * deltaT);
    // Used for the standard deviations
    final Statistics statsX = new Statistics();
    final Statistics statsY = new Statistics();
    final Ticker ticker = ImageJUtils.createTicker(tracks.size(), 1, "Computing track features...");
    // Collect the track categories. Later these are renumbered to a natural sequence from 0.
    final TIntHashSet catSet = new TIntHashSet();
    // final StoredDataStatistics statsAlpha = new StoredDataStatistics(tracks.size());
    // Process each track
    final TIntArrayList lengths = new TIntArrayList(tracks.size());
    for (final Trace track : tracks) {
        // Get xy coordinates
        final int size = track.size();
        if (x.length < size) {
            x = new double[size];
            y = new double[size];
        }
        for (int i = 0; i < size; i++) {
            final PeakResult peak = track.get(i);
            x[i] = distanceConverter.convert(peak.getXPosition());
            y[i] = distanceConverter.convert(peak.getYPosition());
        }
        final int smwm1 = size - wm1;
        final int previousSize = data.size();
        for (int k = 0; k < smwm1; k++) {
            final double[] values = new double[valuesLength];
            data.add(values);
            // middle of even sized windows is between two localisations.
            if (hasCategory) {
                final int cat = track.get(k + mid).getCategory();
                values[4] = cat;
                catSet.add(cat);
            }
            // First point in window = k
            // Last point in window = k + w - 1 = k + wm1
            final int end = k + wm1;
            // 1. Anomalous exponent.
            msdFitter.setData(x, y);
            try {
                msdFitter.fit(k, null);
                // statsAlpha.add(msdFitter.alpha);
                int fitIndex = msdFitter.positiveSlope ? 0 : 2;
                // If better then this is anomalous diffusion
                final double alpha = msdFitter.lvmSolution2.getPoint().getEntry(2);
                values[0] = alpha;
                if (msdFitter.pValue > significance) {
                    fitIndex++;
                }
                fitResult[fitIndex]++;
            // values[0] = msdFitter.alpha;
            // // Debug
            // if (
            // // msdFitter.pValue < 0.2
            // msdFitter.pValue < 0.2 && values[0] < 0
            // // !msdFitter.positiveSlope
            // ) {
            // final RealVector p = msdFitter.lvmSolution2.getPoint();
            // final String title = "anomalous exponent";
            // final Plot plot = new Plot(title, "time (s)", "MSD (um^2)");
            // final double[] t = SimpleArrayUtils.newArray(msdFitter.s.length, deltaT, deltaT);
            // plot.addLabel(0, 0, msdFitter.lvmSolution2.getPoint().toString() + " p="
            // + msdFitter.pValue + ". " + msdFitter.lvmSolution1.getPoint().toString());
            // plot.addPoints(t, msdFitter.s, Plot.CROSS);
            // plot.addPoints(t, msdFitter.model2.value(p).getFirst().toArray(), Plot.LINE);
            // plot.setColor(Color.BLUE);
            // plot.addPoints(t,
            // msdFitter.model1.value(msdFitter.lvmSolution1.getPoint()).getFirst().toArray(),
            // Plot.LINE);
            // plot.setColor(Color.RED);
            // final double[] yy = Arrays.stream(t).map(msdFitter.reg::predict).toArray();
            // plot.addPoints(t, yy, Plot.CIRCLE);
            // System.out.printf("%s : %s", msdFitter.lvmSolution2.getPoint(), values[0]);
            // ImageJUtils.display(title, plot, ImageJUtils.NO_TO_FRONT);
            // System.out.println();
            // }
            } catch (TooManyIterationsException | ConvergenceException ex) {
                if (settings.debug) {
                    ImageJUtils.log("Failed to fit anomalous exponent: " + ex.getMessage());
                }
                // Ignore this and leave as Brownian motion
                values[0] = 1.0;
            }
            // Referenced papers:
            // Hozé, N. H., D. (2017) Statistical methods for large ensembles of super-resolution
            // stochastic single particle trajectories in cell biology.
            // Annual Review of Statistics and Its Application 4, 189-223
            // 
            // Amitai, A., Seeber, A., Gasser, S. M. & Holcman, D. (2017) Visualization of Chromatin
            // Decompaction and Break Site Extrusion as Predicted by Statistical Polymer
            // Modeling of Single-Locus Trajectories. Cell reports 18, 1200-1214
            // 2. Effective diffusion coefficient (Hozé, eq 10).
            // This is the average squared jump distance between successive points
            // divided by 1 / (2 * dimensions * deltaT), i.e. 1 / 4t.
            double sum = 0;
            for (int i = k; i < end; i++) {
                sum += MathUtils.distance2(x[i], y[i], x[i + 1], y[i + 1]);
            }
            values[1] = sum * diffusionCoefficientFactor;
            // 3. Length of confinement (Amitai et al, eq 1).
            // Compute the average of the standard deviation of the position in each dimension.
            statsX.reset();
            statsY.reset();
            for (int i = k; i <= end; i++) {
                statsX.add(x[i]);
                statsY.add(y[i]);
            }
            values[2] = (statsX.getStandardDeviation() + statsY.getStandardDeviation()) / 2;
            // 4. Magnitude of drift vector (Hozé, eq 9).
            // Note: The drift field is given as the expected distance between successive points, i.e.
            // the average step. Since all track windows are the same length this is the same
            // as the distance between the first and last point divided by the number of points.
            // The drift field in each dimension is combined to create a drift norm, i.e. Euclidean
            // distance.
            values[3] = MathUtils.distance(x[k], y[k], x[end], y[end]) / wm1;
        }
        lengths.add(data.size() - previousSize);
        ticker.tick();
    }
    ImageJUtils.finished();
    if (settings.debug) {
        ImageJUtils.log("  +Slope, significant:   %d", fitResult[0]);
        ImageJUtils.log("  +Slope, insignificant: %d", fitResult[1]);
        ImageJUtils.log("  -Slope, significant:   %d", fitResult[2]);
        ImageJUtils.log("  -Slope, insignificant: %d", fitResult[3]);
    }
    ImageJUtils.log("Insignificant anomalous exponents: %d / %d", fitResult[1] + fitResult[3], data.size());
    // System.out.println(statsAlpha.getStatistics().toString());
    final double[][] trackData = data.toArray(new double[0][0]);
    if (hasCategory) {
        final int[] categories = catSet.toArray();
        Arrays.sort(categories);
        // Only remap if non-compact (i.e. not 0 to n)
        final int max = categories[categories.length - 1];
        if (categories[0] != 0 || max + 1 != categories.length) {
            final int[] remap = new int[max + 1];
            for (int i = 0; i < categories.length; i++) {
                remap[categories[i]] = i;
            }
            final int end = valuesLength - 1;
            for (final double[] values : trackData) {
                values[end] = remap[(int) values[end]];
            }
        }
    }
    return Pair.create(lengths.toArray(), trackData);
}

Example 18 with TooManyIterationsException

use of org.apache.commons.math3.exception.TooManyIterationsException in project gdsc by aherbert.

the class FrapAnalysis_PlugIn method fitRecovery.

/**
 * Fit the bleaching curve.
 *
 * <p>The returned fit parameters depend on the model chosen to best fit the data.
 *
 * @param region the region
 * @param y the data
 * @param tau the initial estimate for the general image bleaching rate
 * @param bleachingEvent the index in y for the bleaching event (low point of curve)
 * @param size the size of the region
 * @return fit function and result
 */
private Pair<FrapFunction, Optimum> fitRecovery(int region, float[] y, double tau, int bleachingEvent, int size) {
    final boolean nested = settings.nestedModels;
    // Initial estimates
    // @formatter:off
    // 
    // ---+
    // |
    // |         ---------    A
    // |     ---/
    // |   -/
    // |  /
    // | /
    // |/
    // i0
    // 
    // ------------------------ B
    // 
    // B is background after al bleaching
    // A is the magnitude of the recovery.
    // koff is the recovery rate.
    // i0 is the baseline for the intensity.
    // tau is the bleaching rate over time.
    // @formatter:on
    final double after = y[bleachingEvent];
    double a = y[y.length - 1] - after;
    double i0 = after;
    // Find point where the recovery is half.
    // Use a simple rolling average of 3 to smooth the curve.
    final double half = i0 + a / 2;
    int t = bleachingEvent + 1;
    while (t + 1 < y.length) {
        final double mean = ((double) y[t - 1] + y[t] + y[t + 1]) / 3.0;
        if (mean > half) {
            break;
        }
        t++;
    }
    // half-life (median of an exponential) = ln(2) / koff
    // koff = ln(2) / half-life
    double koff = LN2 / (t - bleachingEvent);
    // Extract curve to fit
    final double[] yy = new double[y.length - bleachingEvent];
    for (int i = 0; i < yy.length; i++) {
        yy[i] = y[i + bleachingEvent];
    }
    final LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
    final RealVector observed = new ArrayRealVector(yy, false);
    final ConvergenceChecker<Evaluation> checker = (iteration, previous, current) -> DoubleEquality.relativeError(previous.getCost(), current.getCost()) < 1e-6;
    // Fit the simple version too (which ignores residual bleaching): y0 + A(1 - exp(-koff * t))
    final FrapFunction model1 = new ReactionLimitedRecoveryFunction(yy.length, timeScale);
    final RealVector start1 = new ArrayRealVector(new double[] { i0, a, koff }, false);
    final ParameterValidator paramValidator = point -> {
        // Do not use MIN_VALUE here to avoid sub-normal numbers
        for (int i = point.getDimension(); i-- > 0; ) {
            if (point.getEntry(i) < Double.MIN_NORMAL) {
                point.setEntry(i, Double.MIN_NORMAL);
            }
        }
        return point;
    };
    final RealMatrix weightMatrix = new DiagonalMatrix(SimpleArrayUtils.newDoubleArray(yy.length, 1.0), false);
    final int maxEvaluations = Integer.MAX_VALUE;
    final int maxIterations = 3000;
    final boolean lazyEvaluation = false;
    final LeastSquaresProblem problem1 = LeastSquaresFactory.create(model1, observed, start1, weightMatrix, checker, maxEvaluations, maxIterations, lazyEvaluation, paramValidator);
    Optimum best1 = null;
    FrapFunction fun1 = model1;
    try {
        best1 = optimizer.optimize(problem1);
        final double[] fit = best1.getPoint().toArray();
        ImageJUtils.log("  Region [%d] reaction limited recovery (ss=%s): f(t) = %s + %s(1 - exp(-%s t)); " + "Half-life = %s", region, MathUtils.rounded(getResidualSumOfSquares(best1)), MathUtils.rounded(fit[0]), MathUtils.rounded(fit[1]), MathUtils.rounded(fit[2]), MathUtils.rounded(LN2 / fit[2]));
    } catch (TooManyIterationsException | ConvergenceException ex) {
        ImageJUtils.log("Failed to fit reaction limited recovery curve: ", ex.getMessage());
    }
    if (best1 != null && nested) {
        // Use the fit as the start point for nested models
        i0 = best1.getPoint().getEntry(0);
        a = best1.getPoint().getEntry(1);
        koff = best1.getPoint().getEntry(2);
        // Fit the simple version but with a general decay of the recovered signal.
        // B can be initialised to the camera offset. Here use a small value.
        final double b = i0 / 100;
        i0 -= b;
        final FrapFunction model2 = new ReactionLimitedRecoveryFunctionB(yy.length, timeScale);
        final RealVector start2 = new ArrayRealVector(new double[] { i0, a, koff, b, tau }, false);
        final LeastSquaresProblem problem2 = LeastSquaresFactory.create(model2, observed, start2, weightMatrix, checker, maxEvaluations, maxIterations, lazyEvaluation, paramValidator);
        try {
            final Optimum lvmSolution = optimizer.optimize(problem2);
            // Check for model improvement
            final double rss1 = getResidualSumOfSquares(best1);
            final double rss2 = getResidualSumOfSquares(lvmSolution);
            double pvalue;
            double f;
            if (rss1 < rss2) {
                f = 0;
                pvalue = 1;
            } else {
                f = RegressionUtils.residualsFStatistic(rss1, start1.getDimension(), rss2, start2.getDimension(), yy.length);
                pvalue = RegressionUtils.residualsFTest(rss1, start1.getDimension(), rss2, start2.getDimension(), yy.length);
            }
            // Optionally log no improvement here...
            final double[] fit = lvmSolution.getPoint().toArray();
            ImageJUtils.log("  Region [%d] reaction limited recovery (ss=%s): f(t) = %s + " + "(%s + %s(1 - exp(-%s t))) * exp(-%s t); Half-life1 = %s; Half-life2 = %s", region, MathUtils.rounded(rss2), MathUtils.rounded(fit[3]), MathUtils.rounded(fit[0]), MathUtils.rounded(fit[1]), MathUtils.rounded(fit[2]), MathUtils.rounded(fit[4]), MathUtils.rounded(LN2 / fit[2]), MathUtils.rounded(LN2 / fit[4]));
            ImageJUtils.log("  Region [%d] : rss1=%s, rss2=%s, p(F-Test=%s) = %s; ", region, MathUtils.rounded(rss1), MathUtils.rounded(rss2), MathUtils.rounded(f), MathUtils.rounded(pvalue));
            if (pvalue < 0.01) {
                // reject null hypothesis that model 2 is not better
                fun1 = model2;
                best1 = lvmSolution;
            }
        } catch (TooManyIterationsException | ConvergenceException ex) {
            ImageJUtils.log("Failed to fit reaction limited recovery curve with decay: ", ex.getMessage());
        }
    }
    // Fit a diffusion limited model.
    // Estimate characteristic diffusion time using the region area and given diffusion coefficient.
    // tD = w^2 / 4D
    // D = w^2 / 4tD
    // Assume the region is a circle
    final double w2 = distanceScale * distanceScale * size / Math.PI;
    // Estimate td using the current fit for i0 and a
    double td = estimateTd(yy, i0, a);
    final FrapFunction model3 = new DiffusionLimitedRecoveryFunction(yy.length, timeScale);
    final RealVector start3 = new ArrayRealVector(new double[] { i0, a, td }, false);
    final LeastSquaresProblem problem3 = LeastSquaresFactory.create(model3, observed, start3, weightMatrix, checker, maxEvaluations, maxIterations, lazyEvaluation, paramValidator);
    Optimum best2 = null;
    FrapFunction fun2 = model3;
    try {
        best2 = optimizer.optimize(problem3);
        final double[] fit = best2.getPoint().toArray();
        final String dT = MathUtils.rounded(fit[2]);
        ImageJUtils.log("  Region [%d] diffusion limited recovery (ss=%s): f(t) = %s + " + "%s(exp(-2*%s/t) * (I0(2*%s/t) + I1(2*%s/t)); D = %s", region, MathUtils.rounded(getResidualSumOfSquares(best2)), MathUtils.rounded(fit[0]), MathUtils.rounded(fit[1]), dT, dT, dT, MathUtils.round(w2 / (4 * timeScale * fit[2])));
    } catch (TooManyIterationsException | ConvergenceException ex) {
        ImageJUtils.log("Failed to fit diffusion limited recovery curve: ", ex.getMessage());
    }
    if (best2 != null && nested) {
        // Use the fit as the start point for nested models
        i0 = best2.getPoint().getEntry(0);
        a = best2.getPoint().getEntry(1);
        td = best2.getPoint().getEntry(2);
        // Fit the simple version but with a general decay of the recovered signal.
        // B can be initialised to the camera offset. Here use a small value.
        final double b = i0 / 100;
        i0 -= b;
        final FrapFunction model4 = new DiffusionLimitedRecoveryFunctionB(yy.length, timeScale);
        final RealVector start4 = new ArrayRealVector(new double[] { i0, a, td, b, tau }, false);
        final LeastSquaresProblem problem4 = LeastSquaresFactory.create(model4, observed, start4, weightMatrix, checker, maxEvaluations, maxIterations, lazyEvaluation, paramValidator);
        try {
            final Optimum lvmSolution = optimizer.optimize(problem4);
            // Check for model improvement
            final double rss1 = getResidualSumOfSquares(best2);
            final double rss2 = getResidualSumOfSquares(lvmSolution);
            double pvalue;
            double f;
            if (rss1 < rss2) {
                f = 0;
                pvalue = 1;
            } else {
                f = RegressionUtils.residualsFStatistic(rss1, start3.getDimension(), rss2, start4.getDimension(), yy.length);
                pvalue = RegressionUtils.residualsFTest(rss1, start3.getDimension(), rss2, start4.getDimension(), yy.length);
            }
            // Optionally log no improvement here...
            final double[] fit = lvmSolution.getPoint().toArray();
            final String dT = MathUtils.rounded(fit[2]);
            ImageJUtils.log("  Region [%d] diffusion limited recovery (ss=%s): f(t) = %s + " + "(%s + %s(exp(-2*%s/t) * (I0(2*%s/t) + I1(2*%s/t))) * exp(-%s t); D = %s; " + "Half-life2 = %s", region, MathUtils.rounded(rss2), MathUtils.rounded(fit[3]), MathUtils.rounded(fit[0]), MathUtils.rounded(fit[1]), dT, dT, dT, MathUtils.rounded(fit[4]), MathUtils.round(w2 / (4 * timeScale * fit[2])), MathUtils.rounded(LN2 / fit[4]));
            ImageJUtils.log("  Region [%d] : rss1=%s, rss2=%s, p(F-Test=%s) = %s; ", region, MathUtils.rounded(rss1), MathUtils.rounded(rss2), MathUtils.rounded(f), MathUtils.rounded(pvalue));
            if (pvalue < 0.01) {
                // reject null hypothesis that model 2 is not better
                fun2 = model4;
                best2 = lvmSolution;
            }
        } catch (TooManyIterationsException | ConvergenceException ex) {
            ImageJUtils.log("Failed to fit diffusion limited recovery curve with decay: ", ex.getMessage());
        }
    }
    if (best1 != null) {
        if (best2 != null) {
            // Return the best. These are not nested models so we cannot use an F-test.
            final double rss1 = getResidualSumOfSquares(best1);
            final double rss2 = getResidualSumOfSquares(best2);
            if (rss2 < rss1) {
                fun1 = fun2;
                best1 = best2;
            }
        }
        return Pair.create(fun1, best1);
    } else if (best2 != null) {
        return Pair.create(fun2, best2);
    }
    // No model
    return null;
}

Example 19 with TooManyIterationsException

use of org.apache.commons.math3.exception.TooManyIterationsException in project gdsc by aherbert.

the class FrapAnalysis_PlugIn method fitBleaching.

/**
 * Fit the bleaching curve {@code f(t) = y0 + B exp(-koff * t)}.
 *
 * @param y the data
 * @param interval the time interval
 * @return {y0, B, koff}
 */
private static double[] fitBleaching(float[] y, double interval) {
    // Initial estimates
    final float[] limits = MathUtils.limits(y);
    final double y0 = limits[0];
    // Find point where the decay is half
    final float half = limits[0] + (limits[1] - limits[0]) / 2;
    int t = 1;
    while (t < y.length && y[t] > half) {
        t++;
    }
    // half-life (median of an exponential) = ln(2) / koff
    // koff = ln(2) / half-life
    final double koff = LN2 / t;
    // Solve for B
    final double b = (half - y0) / Math.exp(-koff * t);
    final LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
    final RealVector observed = new ArrayRealVector(SimpleArrayUtils.toDouble(y), false);
    final ConvergenceChecker<Evaluation> checker = (iteration, previous, current) -> DoubleEquality.relativeError(previous.getCost(), current.getCost()) < 1e-6;
    final MultivariateJacobianFunction model1 = new DecayFunction(y.length, interval);
    final RealVector start1 = new ArrayRealVector(new double[] { y0, b, koff }, false);
    final ParameterValidator paramValidator1 = point -> {
        // Do not use MIN_VALUE here to avoid sub-normal numbers
        for (int i = 0; i < 3; i++) {
            if (point.getEntry(i) < Double.MIN_NORMAL) {
                point.setEntry(i, Double.MIN_NORMAL);
            }
        }
        return point;
    };
    final RealMatrix weightMatrix = new DiagonalMatrix(SimpleArrayUtils.newDoubleArray(y.length, 1.0), false);
    final int maxEvaluations = Integer.MAX_VALUE;
    final int maxIterations = 3000;
    final boolean lazyEvaluation = false;
    final LeastSquaresProblem problem1 = LeastSquaresFactory.create(model1, observed, start1, weightMatrix, checker, maxEvaluations, maxIterations, lazyEvaluation, paramValidator1);
    try {
        final Optimum lvmSolution1 = optimizer.optimize(problem1);
        final RealVector fit1 = lvmSolution1.getPoint();
        return fit1.toArray();
    } catch (TooManyIterationsException | ConvergenceException ex) {
        ImageJUtils.log("Failed to fit bleaching curve: ", ex.getMessage());
        return null;
    }
}

Example 20 with TooManyIterationsException

use of org.apache.commons.math3.exception.TooManyIterationsException 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 nFittedPoints
	 * @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
	 */
public double[] fit(double[] td, double[] ntd, int nFittedPoints, boolean log) {
    blinkingModel = new BlinkingFunction();
    blinkingModel.setLogging(true);
    for (int i = 0; i < nFittedPoints; i++) blinkingModel.addPoint(td[i], ntd[i]);
    // Different convergence thresholds seem to have no effect on the resulting fit, only the number of
    // iterations for convergence
    double initialStepBoundFactor = 100;
    double costRelativeTolerance = 1e-6;
    double parRelativeTolerance = 1e-6;
    double orthoTolerance = 1e-6;
    double threshold = Precision.SAFE_MIN;
    LevenbergMarquardtOptimizer optimiser = new LevenbergMarquardtOptimizer(initialStepBoundFactor, costRelativeTolerance, parRelativeTolerance, orthoTolerance, threshold);
    try {
        double[] obs = blinkingModel.getY();
        //@formatter:off
        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, new MultivariateMatrixFunction() {

            public double[][] value(double[] point) throws IllegalArgumentException {
                return blinkingModel.jacobian(point);
            }
        }).build();
        //@formatter:on
        blinkingModel.setLogging(false);
        Optimum optimum = optimiser.optimize(problem);
        double[] parameters = optimum.getPoint().toArray();
        //double[] exp = blinkingModel.value(parameters);
        double mean = 0;
        for (double d : obs) mean += d;
        mean /= obs.length;
        double ssResiduals = 0, ssTotal = 0;
        for (int i = 0; i < obs.length; i++) {
            //ssResiduals += (obs[i] - exp[i]) * (obs[i] - exp[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) {
            Utils.log("  Fit %d points. R^2 = %s. Adjusted R^2 = %s", obs.length, Utils.rounded(r2, 4), Utils.rounded(adjustedR2, 4));
            Utils.log("  N=%s, nBlink=%s, tOff=%s (%s frames)", Utils.rounded(parameters[0], 4), Utils.rounded(parameters[1], 4), Utils.rounded(parameters[2], 4), Utils.rounded(parameters[2] / msPerFrame, 4));
        }
        return parameters;
    } catch (TooManyIterationsException e) {
        if (log)
            Utils.log("  Failed to fit %d points: Too many iterations: (%s)", blinkingModel.size(), e.getMessage());
        return null;
    } catch (ConvergenceException e) {
        if (log)
            Utils.log("  Failed to fit %d points", blinkingModel.size());
        return null;
    }
}