Examples with PointValuePair org.apache.commons.math3.optim.PointValuePair used on opensource projects

Example 1 with PointValuePair

use of org.apache.commons.math3.optim.PointValuePair in project GDSC-SMLM by aherbert.

the class JumpDistanceAnalysis method doFitJumpDistanceHistogram.

/**
	 * Fit the jump distance histogram using a cumulative sum with the given number of species.
	 * <p>
	 * Results are sorted by the diffusion coefficient ascending.
	 * 
	 * @param jdHistogram
	 *            The cumulative jump distance histogram. X-axis is um^2, Y-axis is cumulative probability. Must be
	 *            monototic ascending.
	 * @param estimatedD
	 *            The estimated diffusion coefficient
	 * @param n
	 *            The number of species in the mixed population
	 * @return Array containing: { D (um^2), Fractions }. Can be null if no fit was made.
	 */
private double[][] doFitJumpDistanceHistogram(double[][] jdHistogram, double estimatedD, int n) {
    calibrated = isCalibrated();
    if (n == 1) {
        // Fit using a single population model
        LevenbergMarquardtOptimizer lvmOptimizer = new LevenbergMarquardtOptimizer();
        try {
            final JumpDistanceCumulFunction function = new JumpDistanceCumulFunction(jdHistogram[0], jdHistogram[1], estimatedD);
            //@formatter:off
            LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(3000).start(function.guess()).target(function.getY()).weight(new DiagonalMatrix(function.getWeights())).model(function, new MultivariateMatrixFunction() {

                public double[][] value(double[] point) throws IllegalArgumentException {
                    return function.jacobian(point);
                }
            }).build();
            //@formatter:on
            Optimum lvmSolution = lvmOptimizer.optimize(problem);
            double[] fitParams = lvmSolution.getPoint().toArray();
            // True for an unweighted fit
            ss = lvmSolution.getResiduals().dotProduct(lvmSolution.getResiduals());
            //ss = calculateSumOfSquares(function.getY(), function.value(fitParams));
            lastIC = ic = Maths.getAkaikeInformationCriterionFromResiduals(ss, function.x.length, 1);
            double[] coefficients = fitParams;
            double[] fractions = new double[] { 1 };
            logger.info("Fit Jump distance (N=1) : %s, SS = %s, IC = %s (%d evaluations)", formatD(fitParams[0]), Maths.rounded(ss, 4), Maths.rounded(ic, 4), lvmSolution.getEvaluations());
            return new double[][] { coefficients, fractions };
        } catch (TooManyIterationsException e) {
            logger.info("LVM optimiser failed to fit (N=1) : Too many iterations : %s", e.getMessage());
        } catch (ConvergenceException e) {
            logger.info("LVM optimiser failed to fit (N=1) : %s", e.getMessage());
        }
    }
    // Uses a weighted sum of n exponential functions, each function models a fraction of the particles.
    // An LVM fit cannot restrict the parameters so the fractions do not go below zero.
    // Use the CustomPowell/CMEASOptimizer which supports bounded fitting.
    MixedJumpDistanceCumulFunctionMultivariate function = new MixedJumpDistanceCumulFunctionMultivariate(jdHistogram[0], jdHistogram[1], estimatedD, n);
    double[] lB = function.getLowerBounds();
    int evaluations = 0;
    PointValuePair constrainedSolution = null;
    MaxEval maxEval = new MaxEval(20000);
    CustomPowellOptimizer powellOptimizer = createCustomPowellOptimizer();
    try {
        // The Powell algorithm can use more general bounds: 0 - Infinity
        constrainedSolution = powellOptimizer.optimize(maxEval, new ObjectiveFunction(function), new InitialGuess(function.guess()), new SimpleBounds(lB, function.getUpperBounds(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY)), new CustomPowellOptimizer.BasisStep(function.step()), GoalType.MINIMIZE);
        evaluations = powellOptimizer.getEvaluations();
        logger.debug("Powell optimiser fit (N=%d) : SS = %f (%d evaluations)", n, constrainedSolution.getValue(), evaluations);
    } catch (TooManyEvaluationsException e) {
        logger.info("Powell optimiser failed to fit (N=%d) : Too many evaluations (%d)", n, powellOptimizer.getEvaluations());
    } catch (TooManyIterationsException e) {
        logger.info("Powell optimiser failed to fit (N=%d) : Too many iterations (%d)", n, powellOptimizer.getIterations());
    } catch (ConvergenceException e) {
        logger.info("Powell optimiser failed to fit (N=%d) : %s", n, e.getMessage());
    }
    if (constrainedSolution == null) {
        logger.info("Trying CMAES optimiser with restarts ...");
        double[] uB = function.getUpperBounds();
        SimpleBounds bounds = new SimpleBounds(lB, uB);
        // The sigma determines the search range for the variables. It should be 1/3 of the initial search region.
        double[] s = new double[lB.length];
        for (int i = 0; i < s.length; i++) s[i] = (uB[i] - lB[i]) / 3;
        OptimizationData sigma = new CMAESOptimizer.Sigma(s);
        OptimizationData popSize = new CMAESOptimizer.PopulationSize((int) (4 + Math.floor(3 * Math.log(function.x.length))));
        // Iterate this for stability in the initial guess
        CMAESOptimizer cmaesOptimizer = createCMAESOptimizer();
        for (int i = 0; i <= fitRestarts; i++) {
            // Try from the initial guess
            try {
                PointValuePair solution = cmaesOptimizer.optimize(new InitialGuess(function.guess()), new ObjectiveFunction(function), GoalType.MINIMIZE, bounds, sigma, popSize, maxEval);
                if (constrainedSolution == null || solution.getValue() < constrainedSolution.getValue()) {
                    evaluations = cmaesOptimizer.getEvaluations();
                    constrainedSolution = solution;
                    logger.debug("CMAES optimiser [%da] fit (N=%d) : SS = %f (%d evaluations)", i, n, solution.getValue(), evaluations);
                }
            } catch (TooManyEvaluationsException e) {
            }
            if (constrainedSolution == null)
                continue;
            // Try from the current optimum
            try {
                PointValuePair solution = cmaesOptimizer.optimize(new InitialGuess(constrainedSolution.getPointRef()), new ObjectiveFunction(function), GoalType.MINIMIZE, bounds, sigma, popSize, maxEval);
                if (solution.getValue() < constrainedSolution.getValue()) {
                    evaluations = cmaesOptimizer.getEvaluations();
                    constrainedSolution = solution;
                    logger.debug("CMAES optimiser [%db] fit (N=%d) : SS = %f (%d evaluations)", i, n, solution.getValue(), evaluations);
                }
            } catch (TooManyEvaluationsException e) {
            }
        }
        if (constrainedSolution != null) {
            // Re-optimise with Powell?
            try {
                PointValuePair solution = powellOptimizer.optimize(maxEval, new ObjectiveFunction(function), new InitialGuess(constrainedSolution.getPointRef()), new SimpleBounds(lB, function.getUpperBounds(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY)), new CustomPowellOptimizer.BasisStep(function.step()), GoalType.MINIMIZE);
                if (solution.getValue() < constrainedSolution.getValue()) {
                    evaluations = cmaesOptimizer.getEvaluations();
                    constrainedSolution = solution;
                    logger.info("Powell optimiser re-fit (N=%d) : SS = %f (%d evaluations)", n, constrainedSolution.getValue(), evaluations);
                }
            } catch (TooManyEvaluationsException e) {
            } catch (TooManyIterationsException e) {
            } catch (ConvergenceException e) {
            }
        }
    }
    if (constrainedSolution == null) {
        logger.info("Failed to fit N=%d", n);
        return null;
    }
    double[] fitParams = constrainedSolution.getPointRef();
    ss = constrainedSolution.getValue();
    // TODO - Try a bounded BFGS optimiser
    // Try and improve using a LVM fit
    final MixedJumpDistanceCumulFunctionGradient functionGradient = new MixedJumpDistanceCumulFunctionGradient(jdHistogram[0], jdHistogram[1], estimatedD, n);
    Optimum lvmSolution;
    LevenbergMarquardtOptimizer lvmOptimizer = new LevenbergMarquardtOptimizer();
    try {
        //@formatter:off
        LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(3000).start(fitParams).target(functionGradient.getY()).weight(new DiagonalMatrix(functionGradient.getWeights())).model(functionGradient, new MultivariateMatrixFunction() {

            public double[][] value(double[] point) throws IllegalArgumentException {
                return functionGradient.jacobian(point);
            }
        }).build();
        //@formatter:on
        lvmSolution = lvmOptimizer.optimize(problem);
        // True for an unweighted fit
        double ss = lvmSolution.getResiduals().dotProduct(lvmSolution.getResiduals());
        // All fitted parameters must be above zero
        if (ss < this.ss && Maths.min(lvmSolution.getPoint().toArray()) > 0) {
            logger.info("  Re-fitting improved the SS from %s to %s (-%s%%)", Maths.rounded(this.ss, 4), Maths.rounded(ss, 4), Maths.rounded(100 * (this.ss - ss) / this.ss, 4));
            fitParams = lvmSolution.getPoint().toArray();
            this.ss = ss;
            evaluations += lvmSolution.getEvaluations();
        }
    } catch (TooManyIterationsException e) {
        logger.error("Failed to re-fit : Too many iterations : %s", e.getMessage());
    } catch (ConvergenceException e) {
        logger.error("Failed to re-fit : %s", e.getMessage());
    }
    // Since the fractions must sum to one we subtract 1 degree of freedom from the number of parameters
    ic = Maths.getAkaikeInformationCriterionFromResiduals(ss, function.x.length, fitParams.length - 1);
    double[] d = new double[n];
    double[] f = new double[n];
    double sum = 0;
    for (int i = 0; i < d.length; i++) {
        f[i] = fitParams[i * 2];
        sum += f[i];
        d[i] = fitParams[i * 2 + 1];
    }
    for (int i = 0; i < f.length; i++) f[i] /= sum;
    // Sort by coefficient size
    sort(d, f);
    double[] coefficients = d;
    double[] fractions = f;
    logger.info("Fit Jump distance (N=%d) : %s (%s), SS = %s, IC = %s (%d evaluations)", n, formatD(d), format(f), Maths.rounded(ss, 4), Maths.rounded(ic, 4), evaluations);
    if (isValid(d, f)) {
        lastIC = ic;
        return new double[][] { coefficients, fractions };
    }
    return null;
}

Example 2 with PointValuePair

use of org.apache.commons.math3.optim.PointValuePair in project GDSC-SMLM by aherbert.

the class JumpDistanceAnalysis method doFitJumpDistancesMLE.

/**
	 * Fit the jump distances using a maximum likelihood estimation with the given number of species.
	 * | *
	 * <p>
	 * Results are sorted by the diffusion coefficient ascending.
	 * 
	 * @param jumpDistances
	 *            The jump distances (in um^2)
	 * @param estimatedD
	 *            The estimated diffusion coefficient
	 * @param n
	 *            The number of species in the mixed population
	 * @return Array containing: { D (um^2), Fractions }. Can be null if no fit was made.
	 */
private double[][] doFitJumpDistancesMLE(double[] jumpDistances, double estimatedD, int n) {
    MaxEval maxEval = new MaxEval(20000);
    CustomPowellOptimizer powellOptimizer = createCustomPowellOptimizer();
    calibrated = isCalibrated();
    if (n == 1) {
        try {
            final JumpDistanceFunction function = new JumpDistanceFunction(jumpDistances, estimatedD);
            // The Powell algorithm can use more general bounds: 0 - Infinity
            SimpleBounds bounds = new SimpleBounds(function.getLowerBounds(), function.getUpperBounds(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY));
            PointValuePair solution = powellOptimizer.optimize(maxEval, new ObjectiveFunction(function), new InitialGuess(function.guess()), bounds, new CustomPowellOptimizer.BasisStep(function.step()), GoalType.MAXIMIZE);
            double[] fitParams = solution.getPointRef();
            ll = solution.getValue();
            lastIC = ic = Maths.getAkaikeInformationCriterion(ll, jumpDistances.length, 1);
            double[] coefficients = fitParams;
            double[] fractions = new double[] { 1 };
            logger.info("Fit Jump distance (N=1) : %s, MLE = %s, IC = %s (%d evaluations)", formatD(fitParams[0]), Maths.rounded(ll, 4), Maths.rounded(ic, 4), powellOptimizer.getEvaluations());
            return new double[][] { coefficients, fractions };
        } catch (TooManyEvaluationsException e) {
            logger.info("Powell optimiser failed to fit (N=1) : Too many evaluation (%d)", powellOptimizer.getEvaluations());
        } catch (TooManyIterationsException e) {
            logger.info("Powell optimiser failed to fit (N=1) : Too many iterations (%d)", powellOptimizer.getIterations());
        } catch (ConvergenceException e) {
            logger.info("Powell optimiser failed to fit (N=1) : %s", e.getMessage());
        }
        return null;
    }
    MixedJumpDistanceFunction function = new MixedJumpDistanceFunction(jumpDistances, estimatedD, n);
    double[] lB = function.getLowerBounds();
    int evaluations = 0;
    PointValuePair constrainedSolution = null;
    try {
        // The Powell algorithm can use more general bounds: 0 - Infinity
        constrainedSolution = powellOptimizer.optimize(maxEval, new ObjectiveFunction(function), new InitialGuess(function.guess()), new SimpleBounds(lB, function.getUpperBounds(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY)), new CustomPowellOptimizer.BasisStep(function.step()), GoalType.MAXIMIZE);
        evaluations = powellOptimizer.getEvaluations();
        logger.debug("Powell optimiser fit (N=%d) : MLE = %f (%d evaluations)", n, constrainedSolution.getValue(), powellOptimizer.getEvaluations());
    } catch (TooManyEvaluationsException e) {
        logger.info("Powell optimiser failed to fit (N=%d) : Too many evaluation (%d)", n, powellOptimizer.getEvaluations());
    } catch (TooManyIterationsException e) {
        logger.info("Powell optimiser failed to fit (N=%d) : Too many iterations (%d)", n, powellOptimizer.getIterations());
    } catch (ConvergenceException e) {
        logger.info("Powell optimiser failed to fit (N=%d) : %s", n, e.getMessage());
    }
    if (constrainedSolution == null) {
        logger.info("Trying CMAES optimiser with restarts ...");
        double[] uB = function.getUpperBounds();
        SimpleBounds bounds = new SimpleBounds(lB, uB);
        // Try a bounded CMAES optimiser since the Powell optimiser appears to be 
        // sensitive to the order of the parameters. It is not good when the fast particle
        // is the minority fraction. Could this be due to too low an upper bound?
        // The sigma determines the search range for the variables. It should be 1/3 of the initial search region.
        double[] s = new double[lB.length];
        for (int i = 0; i < s.length; i++) s[i] = (uB[i] - lB[i]) / 3;
        OptimizationData sigma = new CMAESOptimizer.Sigma(s);
        OptimizationData popSize = new CMAESOptimizer.PopulationSize((int) (4 + Math.floor(3 * Math.log(function.x.length))));
        // Iterate this for stability in the initial guess
        CMAESOptimizer cmaesOptimizer = createCMAESOptimizer();
        for (int i = 0; i <= fitRestarts; i++) {
            // Try from the initial guess
            try {
                PointValuePair solution = cmaesOptimizer.optimize(new InitialGuess(function.guess()), new ObjectiveFunction(function), GoalType.MAXIMIZE, bounds, sigma, popSize, maxEval);
                if (constrainedSolution == null || solution.getValue() > constrainedSolution.getValue()) {
                    evaluations = cmaesOptimizer.getEvaluations();
                    constrainedSolution = solution;
                    logger.debug("CMAES optimiser [%da] fit (N=%d) : MLE = %f (%d evaluations)", i, n, solution.getValue(), evaluations);
                }
            } catch (TooManyEvaluationsException e) {
            }
            if (constrainedSolution == null)
                continue;
            // Try from the current optimum
            try {
                PointValuePair solution = cmaesOptimizer.optimize(new InitialGuess(constrainedSolution.getPointRef()), new ObjectiveFunction(function), GoalType.MAXIMIZE, bounds, sigma, popSize, maxEval);
                if (solution.getValue() > constrainedSolution.getValue()) {
                    evaluations = cmaesOptimizer.getEvaluations();
                    constrainedSolution = solution;
                    logger.debug("CMAES optimiser [%db] fit (N=%d) : MLE = %f (%d evaluations)", i, n, solution.getValue(), evaluations);
                }
            } catch (TooManyEvaluationsException e) {
            }
        }
        if (constrainedSolution != null) {
            try {
                // Re-optimise with Powell?
                PointValuePair solution = powellOptimizer.optimize(maxEval, new ObjectiveFunction(function), new InitialGuess(constrainedSolution.getPointRef()), new SimpleBounds(lB, function.getUpperBounds(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY)), new CustomPowellOptimizer.BasisStep(function.step()), GoalType.MAXIMIZE);
                if (solution.getValue() > constrainedSolution.getValue()) {
                    evaluations = cmaesOptimizer.getEvaluations();
                    constrainedSolution = solution;
                    logger.info("Powell optimiser re-fit (N=%d) : MLE = %f (%d evaluations)", n, constrainedSolution.getValue(), powellOptimizer.getEvaluations());
                }
            } catch (TooManyEvaluationsException e) {
            } catch (TooManyIterationsException e) {
            } catch (ConvergenceException e) {
            }
        }
    }
    if (constrainedSolution == null) {
        logger.info("Failed to fit N=%d", n);
        return null;
    }
    double[] fitParams = constrainedSolution.getPointRef();
    ll = constrainedSolution.getValue();
    // Since the fractions must sum to one we subtract 1 degree of freedom from the number of parameters
    ic = Maths.getAkaikeInformationCriterion(ll, jumpDistances.length, fitParams.length - 1);
    double[] d = new double[n];
    double[] f = new double[n];
    double sum = 0;
    for (int i = 0; i < d.length; i++) {
        f[i] = fitParams[i * 2];
        sum += f[i];
        d[i] = fitParams[i * 2 + 1];
    }
    for (int i = 0; i < f.length; i++) f[i] /= sum;
    // Sort by coefficient size
    sort(d, f);
    double[] coefficients = d;
    double[] fractions = f;
    logger.info("Fit Jump distance (N=%d) : %s (%s), MLE = %s, IC = %s (%d evaluations)", n, formatD(d), format(f), Maths.rounded(ll, 4), Maths.rounded(ic, 4), evaluations);
    if (isValid(d, f)) {
        lastIC = ic;
        return new double[][] { coefficients, fractions };
    }
    return null;
}

Example 3 with PointValuePair

use of org.apache.commons.math3.optim.PointValuePair in project GDSC-SMLM by aherbert.

the class BinomialFitter method fitBinomial.

/**
	 * Fit the binomial distribution (n,p) to the input data. Performs fitting assuming a fixed n value and attempts to
	 * optimise p. All n from minN to maxN are evaluated. If maxN is zero then all possible n from minN are evaluated
	 * until the fit is worse.
	 * 
	 * @param data
	 *            The input data (all value must be positive)
	 * @param minN
	 *            The minimum n to evaluate
	 * @param maxN
	 *            The maximum n to evaluate. Set to zero to evaluate all possible values.
	 * @param zeroTruncated
	 *            True if the model should ignore n=0 (zero-truncated binomial)
	 * @return The best fit (n, p)
	 * @throws IllegalArgumentException
	 *             If any of the input data values are negative
	 */
public double[] fitBinomial(int[] data, int minN, int maxN, boolean zeroTruncated) {
    double[] histogram = getHistogram(data, false);
    final double initialSS = Double.POSITIVE_INFINITY;
    double bestSS = initialSS;
    double[] parameters = null;
    int worse = 0;
    int N = (int) histogram.length - 1;
    if (minN < 1)
        minN = 1;
    if (maxN > 0) {
        if (N > maxN) {
            // Limit the number fitted to maximum
            N = maxN;
        } else if (N < maxN) {
            // Expand the histogram to the maximum
            histogram = Arrays.copyOf(histogram, maxN + 1);
            N = maxN;
        }
    }
    if (minN > N)
        minN = N;
    final double mean = getMean(histogram);
    String name = (zeroTruncated) ? "Zero-truncated Binomial distribution" : "Binomial distribution";
    log("Mean cluster size = %s", Utils.rounded(mean));
    log("Fitting cumulative " + name);
    // score several times in succession)
    for (int n = minN; n <= N; n++) {
        PointValuePair solution = fitBinomial(histogram, mean, n, zeroTruncated);
        if (solution == null)
            continue;
        double p = solution.getPointRef()[0];
        log("Fitted %s : N=%d, p=%s. SS=%g", name, n, Utils.rounded(p), solution.getValue());
        if (bestSS > solution.getValue()) {
            bestSS = solution.getValue();
            parameters = new double[] { n, p };
            worse = 0;
        } else if (bestSS != initialSS) {
            if (++worse >= 3)
                break;
        }
    }
    return parameters;
}

Example 4 with PointValuePair

use of org.apache.commons.math3.optim.PointValuePair in project GDSC-SMLM by aherbert.

the class BinomialFitter method fitBinomial.

/**
	 * Fit the binomial distribution (n,p) to the cumulative histogram. Performs fitting assuming a fixed n value and
	 * attempts to optimise p.
	 * 
	 * @param histogram
	 *            The input histogram
	 * @param mean
	 *            The histogram mean (used to estimate p). Calculated if NaN.
	 * @param n
	 *            The n to evaluate
	 * @param zeroTruncated
	 *            True if the model should ignore n=0 (zero-truncated binomial)
	 * @return The best fit (n, p)
	 * @throws IllegalArgumentException
	 *             If any of the input data values are negative
	 * @throws IllegalArgumentException
	 *             If any fitting a zero truncated binomial and there are no values above zero
	 */
public PointValuePair fitBinomial(double[] histogram, double mean, int n, boolean zeroTruncated) {
    if (Double.isNaN(mean))
        mean = getMean(histogram);
    if (zeroTruncated && histogram[0] > 0) {
        log("Fitting zero-truncated histogram but there are zero values - Renormalising to ignore zero");
        double cumul = 0;
        for (int i = 1; i < histogram.length; i++) cumul += histogram[i];
        if (cumul == 0)
            throw new IllegalArgumentException("Fitting zero-truncated histogram but there are no non-zero values");
        histogram[0] = 0;
        for (int i = 1; i < histogram.length; i++) histogram[i] /= cumul;
    }
    int nFittedPoints = Math.min(histogram.length, n + 1) - ((zeroTruncated) ? 1 : 0);
    if (nFittedPoints < 1) {
        log("No points to fit (%d): Histogram.length = %d, n = %d, zero-truncated = %b", nFittedPoints, histogram.length, n, zeroTruncated);
        return null;
    }
    // The model is only fitting the probability p
    // For a binomial n*p = mean => p = mean/n
    double[] initialSolution = new double[] { FastMath.min(mean / n, 1) };
    // Create the function
    BinomialModelFunction function = new BinomialModelFunction(histogram, n, zeroTruncated);
    double[] lB = new double[1];
    double[] uB = new double[] { 1 };
    SimpleBounds bounds = new SimpleBounds(lB, uB);
    // Fit
    // CMAESOptimizer or BOBYQAOptimizer support bounds
    // CMAESOptimiser based on Matlab code:
    // https://www.lri.fr/~hansen/cmaes.m
    // Take the defaults from the Matlab documentation
    int maxIterations = 2000;
    //Double.NEGATIVE_INFINITY;
    double stopFitness = 0;
    boolean isActiveCMA = true;
    int diagonalOnly = 0;
    int checkFeasableCount = 1;
    RandomGenerator random = new Well19937c();
    boolean generateStatistics = false;
    ConvergenceChecker<PointValuePair> checker = new SimpleValueChecker(1e-6, 1e-10);
    // The sigma determines the search range for the variables. It should be 1/3 of the initial search region.
    OptimizationData sigma = new CMAESOptimizer.Sigma(new double[] { (uB[0] - lB[0]) / 3 });
    OptimizationData popSize = new CMAESOptimizer.PopulationSize((int) (4 + Math.floor(3 * Math.log(2))));
    try {
        PointValuePair solution = null;
        boolean noRefit = maximumLikelihood;
        if (n == 1 && zeroTruncated) {
            // No need to fit
            solution = new PointValuePair(new double[] { 1 }, 0);
            noRefit = true;
        } else {
            GoalType goalType = (maximumLikelihood) ? GoalType.MAXIMIZE : GoalType.MINIMIZE;
            // Iteratively fit
            CMAESOptimizer opt = new CMAESOptimizer(maxIterations, stopFitness, isActiveCMA, diagonalOnly, checkFeasableCount, random, generateStatistics, checker);
            for (int iteration = 0; iteration <= fitRestarts; iteration++) {
                try {
                    // Start from the initial solution
                    PointValuePair result = opt.optimize(new InitialGuess(initialSolution), new ObjectiveFunction(function), goalType, bounds, sigma, popSize, new MaxIter(maxIterations), new MaxEval(maxIterations * 2));
                    //		opt.getEvaluations());
                    if (solution == null || result.getValue() < solution.getValue()) {
                        solution = result;
                    }
                } catch (TooManyEvaluationsException e) {
                } catch (TooManyIterationsException e) {
                }
                if (solution == null)
                    continue;
                try {
                    // Also restart from the current optimum
                    PointValuePair result = opt.optimize(new InitialGuess(solution.getPointRef()), new ObjectiveFunction(function), goalType, bounds, sigma, popSize, new MaxIter(maxIterations), new MaxEval(maxIterations * 2));
                    //		opt.getEvaluations());
                    if (result.getValue() < solution.getValue()) {
                        solution = result;
                    }
                } catch (TooManyEvaluationsException e) {
                } catch (TooManyIterationsException e) {
                }
            }
            if (solution == null)
                return null;
        }
        if (noRefit) {
            // Although we fit the log-likelihood, return the sum-of-squares to allow 
            // comparison across different n
            double p = solution.getPointRef()[0];
            double ss = 0;
            double[] obs = function.p;
            double[] exp = function.getP(p);
            for (int i = 0; i < obs.length; i++) ss += (obs[i] - exp[i]) * (obs[i] - exp[i]);
            return new PointValuePair(solution.getPointRef(), ss);
        } else // We can do a LVM refit if the number of fitted points is more than 1
        if (nFittedPoints > 1) {
            // Improve SS fit with a gradient based LVM optimizer
            LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
            try {
                final BinomialModelFunctionGradient gradientFunction = new BinomialModelFunctionGradient(histogram, n, zeroTruncated);
                //@formatter:off
                LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(3000).start(solution.getPointRef()).target(gradientFunction.p).weight(new DiagonalMatrix(gradientFunction.getWeights())).model(gradientFunction, new MultivariateMatrixFunction() {

                    public double[][] value(double[] point) throws IllegalArgumentException {
                        return gradientFunction.jacobian(point);
                    }
                }).build();
                //@formatter:on
                Optimum lvmSolution = optimizer.optimize(problem);
                // Check the pValue is valid since the LVM is not bounded.
                double p = lvmSolution.getPoint().getEntry(0);
                if (p <= 1 && p >= 0) {
                    // True if the weights are 1
                    double ss = lvmSolution.getResiduals().dotProduct(lvmSolution.getResiduals());
                    //	ss += (obs[i] - exp[i]) * (obs[i] - exp[i]);
                    if (ss < solution.getValue()) {
                        //		Utils.rounded(100 * (solution.getValue() - ss) / solution.getValue(), 4));
                        return new PointValuePair(lvmSolution.getPoint().toArray(), ss);
                    }
                }
            } catch (TooManyIterationsException e) {
                log("Failed to re-fit: Too many iterations: %s", e.getMessage());
            } catch (ConvergenceException e) {
                log("Failed to re-fit: %s", e.getMessage());
            } catch (Exception e) {
            // Ignore this ...
            }
        }
        return solution;
    } catch (Exception e) {
        log("Failed to fit Binomial distribution with N=%d : %s", n, e.getMessage());
    }
    return null;
}

Example 5 with PointValuePair

use of org.apache.commons.math3.optim.PointValuePair in project GDSC-SMLM by aherbert.

the class CustomPowellOptimizer method doOptimize.

/** {@inheritDoc} */
@Override
protected PointValuePair doOptimize() {
    final GoalType goal = getGoalType();
    final double[] guess = getStartPoint();
    final int n = guess.length;
    // Mark when we have modified the basis vectors
    boolean nonBasis = false;
    double[][] direc = createBasisVectors(n);
    final ConvergenceChecker<PointValuePair> checker = getConvergenceChecker();
    //int resets = 0;
    //PointValuePair solution = null;
    //PointValuePair finalSolution = null;
    //int solutionIter = 0, solutionEval = 0;
    //double startValue = 0;
    //try
    //{
    double[] x = guess;
    // Ensure the point is within bounds
    applyBounds(x);
    double fVal = computeObjectiveValue(x);
    //startValue = fVal;
    double[] x1 = x.clone();
    while (true) {
        incrementIterationCount();
        final double fX = fVal;
        double fX2 = 0;
        double delta = 0;
        int bigInd = 0;
        for (int i = 0; i < n; i++) {
            fX2 = fVal;
            final UnivariatePointValuePair optimum = line.search(x, direc[i]);
            fVal = optimum.getValue();
            x = newPoint(x, direc[i], optimum.getPoint());
            if ((fX2 - fVal) > delta) {
                delta = fX2 - fVal;
                bigInd = i;
            }
        }
        boolean stop = false;
        if (positionChecker != null) {
            // Check for convergence on the position
            stop = positionChecker.converged(x1, x);
        }
        if (!stop) {
            // Check if we have improved from an impossible position
            if (Double.isInfinite(fX) || Double.isNaN(fX)) {
                if (Double.isInfinite(fVal) || Double.isNaN(fVal)) {
                    // Nowhere to go 
                    stop = true;
                }
            // else: this is better as we now have a value, so continue
            } else {
                stop = DoubleEquality.almostEqualRelativeOrAbsolute(fX, fVal, relativeThreshold, absoluteThreshold);
            }
        }
        final PointValuePair previous = new PointValuePair(x1, fX);
        final PointValuePair current = new PointValuePair(x, fVal);
        if (!stop && checker != null) {
            // User-defined stopping criteria.
            stop = checker.converged(getIterations(), previous, current);
        }
        boolean reset = false;
        if (stop) {
            // Only allow convergence using the basis vectors, i.e. we cannot move along any dimension
            if (basisConvergence && nonBasis) {
                // Reset to the basis vectors and continue
                reset = true;
            //resets++;
            } else {
                //System.out.printf("Resets = %d\n", resets);
                final PointValuePair answer;
                if (goal == GoalType.MINIMIZE) {
                    answer = (fVal < fX) ? current : previous;
                } else {
                    answer = (fVal > fX) ? current : previous;
                }
                return answer;
            // XXX Debugging
            // Continue the algorithm to see how far it goes
            //if (solution == null)
            //{
            //	solution = answer;
            //	solutionIter = getIterations();
            //	solutionEval = getEvaluations();
            //}
            //finalSolution = answer;
            }
        }
        if (reset) {
            direc = createBasisVectors(n);
            nonBasis = false;
        }
        final double[] d = new double[n];
        final double[] x2 = new double[n];
        for (int i = 0; i < n; i++) {
            d[i] = x[i] - x1[i];
            x2[i] = x[i] + d[i];
        }
        applyBounds(x2);
        x1 = x.clone();
        fX2 = computeObjectiveValue(x2);
        // See if we can continue along the overall search direction to find a better value
        if (fX > fX2) {
            // Check if:
            // 1. The decrease along the average direction was not due to any single direction's decrease
            // 2. There is a substantial second derivative along the average direction and we are close to
            // it minimum
            double t = 2 * (fX + fX2 - 2 * fVal);
            double temp = fX - fVal - delta;
            t *= temp * temp;
            temp = fX - fX2;
            t -= delta * temp * temp;
            if (t < 0.0) {
                final UnivariatePointValuePair optimum = line.search(x, d);
                fVal = optimum.getValue();
                if (reset) {
                    x = newPoint(x, d, optimum.getPoint());
                    continue;
                } else {
                    final double[][] result = newPointAndDirection(x, d, optimum.getPoint());
                    x = result[0];
                    final int lastInd = n - 1;
                    direc[bigInd] = direc[lastInd];
                    direc[lastInd] = result[1];
                    nonBasis = true;
                }
            }
        }
    }
//}
//catch (RuntimeException e)
//{
//	if (solution != null)
//	{
//		System.out.printf("Start %f : Initial %f (%d,%d) : Final %f (%d,%d) : %f\n", startValue,
//				solution.getValue(), solutionIter, solutionEval, finalSolution.getValue(), getIterations(),
//				getEvaluations(), DoubleEquality.relativeError(finalSolution.getValue(), solution.getValue()));
//		return finalSolution;
//	}
//	throw e;
//}
}