Examples with Evaluation org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation used on opensource projects

Example 1 with Evaluation

use of org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation 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 2 with Evaluation

use of org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation in project GDSC-SMLM by aherbert.

the class BenchmarkFilterAnalysis method filterAnalysis.

private int filterAnalysis(FilterSet filterSet, int setNumber, DirectFilter currentOptimum, double rangeReduction) {
    // Check if the filters are the same so allowing optimisation
    final boolean allSameType = filterSet.allSameType();
    this.ga_resultsList = resultsList;
    Chromosome<FilterScore> best = null;
    String algorithm = "";
    // All the search algorithms use search dimensions.
    // Create search dimensions if needed (these are used for testing if the optimum is at the limit).
    ss_filter = null;
    ss_lower = null;
    ss_upper = null;
    FixedDimension[] originalDimensions = null;
    boolean rangeInput = false;
    boolean[] disabled = null;
    double[][] seed = null;
    boolean nonInteractive = false;
    if (allSameType) {
        // There should always be 1 filter
        ss_filter = (DirectFilter) filterSet.getFilters().get(0);
        int n = ss_filter.getNumberOfParameters();
        // Option to configure a range
        rangeInput = filterSet.getName().contains("Range");
        double[] range = new double[n];
        if (rangeInput && filterSet.size() == 4) {
            originalDimensions = new FixedDimension[n];
            // This is used as min/lower/upper/max
            final Filter minF = ss_filter;
            final Filter lowerF = filterSet.getFilters().get(1);
            final Filter upperF = filterSet.getFilters().get(2);
            final Filter maxF = filterSet.getFilters().get(3);
            for (int i = 0; i < n; i++) {
                double min = minF.getParameterValue(i);
                double lower = lowerF.getParameterValue(i);
                double upper = upperF.getParameterValue(i);
                range[i] = upper - lower;
                double max = maxF.getParameterValue(i);
                double minIncrement = ss_filter.getParameterIncrement(i);
                try {
                    originalDimensions[i] = new FixedDimension(min, max, minIncrement, lower, upper);
                } catch (IllegalArgumentException e) {
                    Utils.log(TITLE + " : Unable to configure dimension [%d] %s: " + e.getMessage(), i, ss_filter.getParameterName(i));
                    originalDimensions = null;
                    rangeInput = false;
                    break;
                }
            }
        }
        if (rangeInput && (filterSet.size() == 3 || filterSet.size() == 2)) {
            originalDimensions = new FixedDimension[n];
            // This is used as lower/upper[/increment]
            final Filter lowerF = ss_filter;
            final Filter upperF = filterSet.getFilters().get(1);
            for (int i = 0; i < n; i++) {
                // Do not disable if the increment is not set. This is left to the user to decide.
                //					if (incF.getParameterValue(i) == incF.getDisabledParameterValue(i) ||
                //							Double.isInfinite(incF.getParameterValue(i)))
                //					{
                //						// Not enabled
                //						dimensions[i] = new SearchDimension(incF.getDisabledParameterValue(i));
                //						continue;
                //					}
                double lower = lowerF.getParameterValue(i);
                double upper = upperF.getParameterValue(i);
                range[i] = upper - lower;
                ParameterType type = ss_filter.getParameterType(i);
                double min = BenchmarkSpotFit.getMin(type);
                double max = BenchmarkSpotFit.getMax(type);
                double minIncrement = ss_filter.getParameterIncrement(i);
                try {
                    originalDimensions[i] = new FixedDimension(min, max, minIncrement, lower, upper);
                } catch (IllegalArgumentException e) {
                    Utils.log(TITLE + " : Unable to configure dimension [%d] %s: " + e.getMessage(), i, ss_filter.getParameterName(i));
                    originalDimensions = null;
                    rangeInput = false;
                    break;
                }
            }
        }
        // Get the dimensions from the filters
        if (originalDimensions == null) {
            originalDimensions = new FixedDimension[n];
            // Allow inputing a filter set (e.g. saved from previous optimisation)
            // Find the limits in the current scores
            final double[] lower = ss_filter.getParameters().clone();
            final double[] upper = lower.clone();
            // Allow the SearchSpace algorithms to be seeded with an initial population 
            // for the first evaluation of the optimum. This is done when the input filter 
            // set is not a range.
            seed = new double[filterSet.size()][];
            int c = 0;
            for (Filter f : filterSet.getFilters()) {
                final double[] point = f.getParameters();
                seed[c++] = point;
                for (int j = 0; j < lower.length; j++) {
                    if (lower[j] > point[j])
                        lower[j] = point[j];
                    if (upper[j] < point[j])
                        upper[j] = point[j];
                }
            }
            // Min/max must be set using values from BenchmarkSpotFit.
            for (int i = 0; i < n; i++) {
                if (lower[i] == upper[i]) {
                    // Not enabled
                    originalDimensions[i] = new FixedDimension(lower[i]);
                    continue;
                }
                ParameterType type = ss_filter.getParameterType(i);
                double min = BenchmarkSpotFit.getMin(type);
                double max = BenchmarkSpotFit.getMax(type);
                double minIncrement = ss_filter.getParameterIncrement(i);
                if (min > lower[i])
                    min = lower[i];
                if (max < upper[i])
                    max = upper[i];
                try {
                    originalDimensions[i] = new FixedDimension(min, max, minIncrement, lower[i], upper[i]);
                } catch (IllegalArgumentException e) {
                    Utils.log(TITLE + " : Unable to configure dimension [%d] %s: " + e.getMessage(), i, ss_filter.getParameterName(i));
                    originalDimensions = null;
                    break;
                }
            }
            if (originalDimensions == null) {
                // Failed to work out the dimensions. No optimisation will be possible.
                // Sort so that the filters are in a nice order for reporting
                filterSet.sort();
                // This will not be used when the dimensions are null
                seed = null;
            }
        }
        if (originalDimensions != null) {
            // Use the current optimum if we are doing a range optimisation
            if (currentOptimum != null && rangeInput && currentOptimum.getType().equals(ss_filter.getType()) && evolve != 0) {
                // Suppress dialogs and use the current settings
                nonInteractive = true;
                double[] p = currentOptimum.getParameters();
                // Range search uses SearchDimension and we must centre on the optimum after creation.							
                for (int i = 0; i < originalDimensions.length; i++) {
                    double centre = p[i];
                    double r = 0;
                    if (originalDimensions[i].isActive()) {
                        // Set the range around the centre.
                        // This uses the range for each param when we read the filters.
                        r = range[i];
                        // Optionally reduce the width of the dimensions. 
                        if (rangeReduction > 0 && rangeReduction < 1)
                            r *= rangeReduction;
                    }
                    double lower = centre - r * 0.5;
                    double upper = centre + r * 0.5;
                    originalDimensions[i] = originalDimensions[i].create(lower, upper);
                }
            }
            // Store the dimensions so we can do an 'at limit' check
            disabled = new boolean[originalDimensions.length];
            ss_lower = new double[originalDimensions.length];
            ss_upper = new double[originalDimensions.length];
            for (int i = 0; i < disabled.length; i++) {
                disabled[i] = !originalDimensions[i].isActive();
                ss_lower[i] = originalDimensions[i].lower;
                ss_upper[i] = originalDimensions[i].upper;
            }
        }
    } else {
        // Sort so that the filters are in a nice order for reporting
        filterSet.sort();
    }
    analysisStopWatch = StopWatch.createStarted();
    if (evolve == 1 && originalDimensions != null) {
        // Collect parameters for the genetic algorithm
        pauseFilterTimer();
        // Remember the step size settings
        double[] stepSize = stepSizeMap.get(setNumber);
        if (stepSize == null || stepSize.length != ss_filter.length()) {
            stepSize = ss_filter.mutationStepRange().clone();
            for (int j = 0; j < stepSize.length; j++) stepSize[j] *= delta;
            // See if the same number of parameters have been optimised in other algorithms
            boolean[] enabled = searchRangeMap.get(setNumber);
            if (enabled != null && enabled.length == stepSize.length) {
                for (int j = 0; j < stepSize.length; j++) if (!enabled[j])
                    stepSize[j] *= -1;
            }
        }
        GenericDialog gd = null;
        int[] indices = ss_filter.getChromosomeParameters();
        boolean runAlgorithm = nonInteractive;
        if (!nonInteractive) {
            // Ask the user for the mutation step parameters.
            gd = new GenericDialog(TITLE);
            String prefix = setNumber + "_";
            gd.addMessage("Configure the genetic algorithm for [" + setNumber + "] " + filterSet.getName());
            gd.addNumericField(prefix + "Population_size", populationSize, 0);
            gd.addNumericField(prefix + "Failure_limit", failureLimit, 0);
            gd.addNumericField(prefix + "Tolerance", tolerance, -1);
            gd.addNumericField(prefix + "Converged_count", convergedCount, 0);
            gd.addSlider(prefix + "Mutation_rate", 0.05, 1, mutationRate);
            gd.addSlider(prefix + "Crossover_rate", 0.05, 1, crossoverRate);
            gd.addSlider(prefix + "Mean_children", 0.05, 3, meanChildren);
            gd.addSlider(prefix + "Selection_fraction", 0.05, 0.5, selectionFraction);
            gd.addCheckbox(prefix + "Ramped_selection", rampedSelection);
            gd.addCheckbox(prefix + "Save_option", saveOption);
            gd.addMessage("Configure the step size for each parameter");
            for (int j = 0; j < indices.length; j++) {
                // Do not mutate parameters that were not expanded, i.e. the input did not vary them.
                final double step = (originalDimensions[indices[j]].isActive()) ? stepSize[j] : 0;
                gd.addNumericField(getDialogName(prefix, ss_filter, indices[j]), step, 2);
            }
            gd.showDialog();
            runAlgorithm = !gd.wasCanceled();
        }
        if (runAlgorithm) {
            // Used to create random sample
            FixedDimension[] dimensions = Arrays.copyOf(originalDimensions, originalDimensions.length);
            if (!nonInteractive) {
                populationSize = (int) Math.abs(gd.getNextNumber());
                if (populationSize < 10)
                    populationSize = 10;
                failureLimit = (int) Math.abs(gd.getNextNumber());
                tolerance = gd.getNextNumber();
                // Allow negatives
                convergedCount = (int) gd.getNextNumber();
                mutationRate = Math.abs(gd.getNextNumber());
                crossoverRate = Math.abs(gd.getNextNumber());
                meanChildren = Math.abs(gd.getNextNumber());
                selectionFraction = Math.abs(gd.getNextNumber());
                rampedSelection = gd.getNextBoolean();
                saveOption = gd.getNextBoolean();
                for (int j = 0; j < indices.length; j++) {
                    stepSize[j] = gd.getNextNumber();
                }
                // Store for repeat analysis
                stepSizeMap.put(setNumber, stepSize);
            }
            for (int j = 0; j < indices.length; j++) {
                // A zero step size will keep the parameter but prevent range mutation.
                if (stepSize[j] < 0) {
                    dimensions[indices[j]] = new FixedDimension(ss_filter.getDisabledParameterValue(indices[j]));
                    disabled[indices[j]] = true;
                }
            }
            //				// Reset negatives to zero
            //				stepSize = stepSize.clone();
            //				for (int j = 0; j < stepSize.length; j++)
            //					if (stepSize[j] < 0)
            //						stepSize[j] = 0;
            // Create the genetic algorithm
            RandomDataGenerator random = new RandomDataGenerator(new Well44497b());
            SimpleMutator<FilterScore> mutator = new SimpleMutator<FilterScore>(random, mutationRate);
            // Override the settings with the step length, a min of zero and the configured upper
            double[] upper = ss_filter.upperLimit();
            mutator.overrideChromosomeSettings(stepSize, new double[stepSize.length], upper);
            Recombiner<FilterScore> recombiner = new SimpleRecombiner<FilterScore>(random, crossoverRate, meanChildren);
            SelectionStrategy<FilterScore> selectionStrategy;
            // If the initial population is huge ensure that the first selection culls to the correct size
            final int selectionMax = (int) (selectionFraction * populationSize);
            if (rampedSelection)
                selectionStrategy = new RampedSelectionStrategy<FilterScore>(random, selectionFraction, selectionMax);
            else
                selectionStrategy = new SimpleSelectionStrategy<FilterScore>(random, selectionFraction, selectionMax);
            ToleranceChecker<FilterScore> ga_checker = new InterruptChecker(tolerance, tolerance * 1e-3, convergedCount);
            // Create new random filters if the population is initially below the population size
            List<Filter> filters = filterSet.getFilters();
            if (filterSet.size() < populationSize) {
                filters = new ArrayList<Filter>(populationSize);
                // Add the existing filters if they are not a range input file
                if (!rangeInput)
                    filters.addAll(filterSet.getFilters());
                // Add current optimum to seed
                if (nonInteractive)
                    filters.add(currentOptimum);
                // The GA does not use the min interval grid so sample without rounding
                double[][] sample = SearchSpace.sampleWithoutRounding(dimensions, populationSize - filters.size(), null);
                filters.addAll(searchSpaceToFilters(sample));
            }
            ga_population = new Population<FilterScore>(filters);
            ga_population.setPopulationSize(populationSize);
            ga_population.setFailureLimit(failureLimit);
            selectionStrategy.setTracker(this);
            // Evolve
            algorithm = EVOLVE[evolve];
            ga_statusPrefix = algorithm + " [" + setNumber + "] " + filterSet.getName() + " ... ";
            ga_iteration = 0;
            ga_population.setTracker(this);
            createGAWindow();
            resumeFilterTimer();
            best = ga_population.evolve(mutator, recombiner, this, selectionStrategy, ga_checker);
            if (best != null) {
                // In case optimisation was stopped
                IJ.resetEscape();
                // The GA may produce coordinates off the min interval grid
                best = enumerateMinInterval(best, stepSize, indices);
                // Now update the filter set for final assessment
                filterSet = new FilterSet(filterSet.getName(), populationToFilters(ga_population.getIndividuals()));
                // Option to save the filters
                if (saveOption)
                    saveFilterSet(filterSet, setNumber, !nonInteractive);
            }
        } else
            resumeFilterTimer();
    }
    if ((evolve == 2 || evolve == 4) && originalDimensions != null) {
        // Collect parameters for the range search algorithm
        pauseFilterTimer();
        boolean isStepSearch = evolve == 4;
        // The step search should use a multi-dimension refinement and no range reduction
        SearchSpace.RefinementMode myRefinementMode = SearchSpace.RefinementMode.MULTI_DIMENSION;
        // Remember the enabled settings
        boolean[] enabled = searchRangeMap.get(setNumber);
        int n = ss_filter.getNumberOfParameters();
        if (enabled == null || enabled.length != n) {
            enabled = new boolean[n];
            Arrays.fill(enabled, true);
            // See if the same number of parameters have been optimised in other algorithms
            double[] stepSize = stepSizeMap.get(setNumber);
            if (stepSize != null && enabled.length == stepSize.length) {
                for (int j = 0; j < stepSize.length; j++) if (stepSize[j] < 0)
                    enabled[j] = false;
            }
        }
        GenericDialog gd = null;
        boolean runAlgorithm = nonInteractive;
        if (!nonInteractive) {
            // Ask the user for the search parameters.
            gd = new GenericDialog(TITLE);
            String prefix = setNumber + "_";
            gd.addMessage("Configure the " + EVOLVE[evolve] + " algorithm for [" + setNumber + "] " + filterSet.getName());
            gd.addSlider(prefix + "Width", 1, 5, rangeSearchWidth);
            if (!isStepSearch) {
                gd.addCheckbox(prefix + "Save_option", saveOption);
                gd.addNumericField(prefix + "Max_iterations", maxIterations, 0);
                String[] modes = SettingsManager.getNames((Object[]) SearchSpace.RefinementMode.values());
                gd.addSlider(prefix + "Reduce", 0.01, 0.99, rangeSearchReduce);
                gd.addChoice("Refinement", modes, modes[refinementMode]);
            }
            gd.addNumericField(prefix + "Seed_size", seedSize, 0);
            // Add choice of fields to optimise
            for (int i = 0; i < n; i++) gd.addCheckbox(getDialogName(prefix, ss_filter, i), enabled[i]);
            gd.showDialog();
            runAlgorithm = !gd.wasCanceled();
        }
        if (runAlgorithm) {
            SearchDimension[] dimensions = new SearchDimension[n];
            if (!nonInteractive) {
                rangeSearchWidth = (int) gd.getNextNumber();
                if (!isStepSearch) {
                    saveOption = gd.getNextBoolean();
                    maxIterations = (int) gd.getNextNumber();
                    refinementMode = gd.getNextChoiceIndex();
                    rangeSearchReduce = gd.getNextNumber();
                }
                seedSize = (int) gd.getNextNumber();
                for (int i = 0; i < n; i++) enabled[i] = gd.getNextBoolean();
                searchRangeMap.put(setNumber, enabled);
            }
            if (!isStepSearch)
                myRefinementMode = SearchSpace.RefinementMode.values()[refinementMode];
            for (int i = 0; i < n; i++) {
                if (enabled[i]) {
                    try {
                        dimensions[i] = originalDimensions[i].create(rangeSearchWidth);
                        dimensions[i].setPad(true);
                        // Prevent range reduction so that the step search just does a single refinement step
                        dimensions[i].setReduceFactor((isStepSearch) ? 1 : rangeSearchReduce);
                        // Centre on current optimum
                        if (nonInteractive)
                            dimensions[i].setCentre(currentOptimum.getParameterValue(i));
                    } catch (IllegalArgumentException e) {
                        IJ.error(TITLE, String.format("Unable to configure dimension [%d] %s: " + e.getMessage(), i, ss_filter.getParameterName(i)));
                        return -1;
                    }
                } else {
                    dimensions[i] = new SearchDimension(ss_filter.getDisabledParameterValue(i));
                }
            }
            for (int i = 0; i < disabled.length; i++) disabled[i] = !dimensions[i].active;
            // Check the number of combinations is OK
            long combinations = SearchSpace.countCombinations(dimensions);
            if (!nonInteractive && combinations > 10000) {
                gd = new GenericDialog(TITLE);
                gd.addMessage(String.format("%d combinations for the configured dimensions.\n \nClick 'Yes' to optimise.", combinations));
                gd.enableYesNoCancel();
                gd.hideCancelButton();
                gd.showDialog();
                if (!gd.wasOKed()) {
                    combinations = 0;
                }
            }
            if (combinations == 0) {
                resumeFilterTimer();
            } else {
                algorithm = EVOLVE[evolve] + " " + rangeSearchWidth;
                ga_statusPrefix = algorithm + " [" + setNumber + "] " + filterSet.getName() + " ... ";
                ga_iteration = 0;
                es_optimum = null;
                SearchSpace ss = new SearchSpace();
                ss.setTracker(this);
                if (seedSize > 0) {
                    double[][] sample;
                    // Add current optimum to seed
                    if (nonInteractive) {
                        sample = new double[1][];
                        sample[0] = currentOptimum.getParameters();
                        seed = merge(seed, sample);
                    }
                    int size = (seed == null) ? 0 : seed.length;
                    // Sample without rounding as the seed will be rounded
                    sample = SearchSpace.sampleWithoutRounding(dimensions, seedSize - size, null);
                    seed = merge(seed, sample);
                }
                // Note: If we have an optimum and we are not seeding this should not matter as the dimensions 
                // have been centred on the current optimum					
                ss.seed(seed);
                ConvergenceChecker<FilterScore> checker = new InterruptConvergenceChecker(0, 0, maxIterations);
                createGAWindow();
                resumeFilterTimer();
                SearchResult<FilterScore> optimum = ss.search(dimensions, this, checker, myRefinementMode);
                if (optimum != null) {
                    // In case optimisation was stopped
                    IJ.resetEscape();
                    best = ((SimpleFilterScore) optimum.score).r.filter;
                    if (seedSize > 0) {
                    // Not required as the search now respects the min interval
                    // The optimum may be off grid if it was from the seed
                    //best = enumerateMinInterval(best, enabled);
                    }
                    // Now update the filter set for final assessment
                    filterSet = new FilterSet(filterSet.getName(), searchSpaceToFilters((DirectFilter) best, ss.getSearchSpace()));
                    // Option to save the filters
                    if (saveOption)
                        saveFilterSet(filterSet, setNumber, !nonInteractive);
                }
            }
        } else
            resumeFilterTimer();
    }
    if (evolve == 3 && originalDimensions != null) {
        // Collect parameters for the enrichment search algorithm
        pauseFilterTimer();
        boolean[] enabled = searchRangeMap.get(setNumber);
        int n = ss_filter.getNumberOfParameters();
        if (enabled == null || enabled.length != n) {
            enabled = new boolean[n];
            Arrays.fill(enabled, true);
            // See if the same number of parameters have been optimised in other algorithms
            double[] stepSize = stepSizeMap.get(setNumber);
            if (stepSize != null && enabled.length == stepSize.length) {
                for (int j = 0; j < stepSize.length; j++) if (stepSize[j] < 0)
                    enabled[j] = false;
            }
        }
        GenericDialog gd = null;
        boolean runAlgorithm = nonInteractive;
        if (!nonInteractive) {
            // Ask the user for the search parameters.
            gd = new GenericDialog(TITLE);
            String prefix = setNumber + "_";
            gd.addMessage("Configure the enrichment search algorithm for [" + setNumber + "] " + filterSet.getName());
            gd.addCheckbox(prefix + "Save_option", saveOption);
            gd.addNumericField(prefix + "Max_iterations", maxIterations, 0);
            gd.addNumericField(prefix + "Converged_count", convergedCount, 0);
            gd.addNumericField(prefix + "Samples", enrichmentSamples, 0);
            gd.addSlider(prefix + "Fraction", 0.01, 0.99, enrichmentFraction);
            gd.addSlider(prefix + "Padding", 0, 0.99, enrichmentPadding);
            // Add choice of fields to optimise
            for (int i = 0; i < n; i++) gd.addCheckbox(getDialogName(prefix, ss_filter, i), enabled[i]);
            gd.showDialog();
            runAlgorithm = !gd.wasCanceled();
        }
        if (runAlgorithm) {
            FixedDimension[] dimensions = Arrays.copyOf(originalDimensions, originalDimensions.length);
            if (!nonInteractive) {
                saveOption = gd.getNextBoolean();
                maxIterations = (int) gd.getNextNumber();
                convergedCount = (int) gd.getNextNumber();
                enrichmentSamples = (int) gd.getNextNumber();
                enrichmentFraction = gd.getNextNumber();
                enrichmentPadding = gd.getNextNumber();
                for (int i = 0; i < n; i++) enabled[i] = gd.getNextBoolean();
                searchRangeMap.put(setNumber, enabled);
            }
            for (int i = 0; i < n; i++) {
                if (!enabled[i])
                    dimensions[i] = new FixedDimension(ss_filter.getDisabledParameterValue(i));
            }
            for (int i = 0; i < disabled.length; i++) disabled[i] = !dimensions[i].active;
            algorithm = EVOLVE[evolve];
            ga_statusPrefix = algorithm + " [" + setNumber + "] " + filterSet.getName() + " ... ";
            ga_iteration = 0;
            es_optimum = null;
            SearchSpace ss = new SearchSpace();
            ss.setTracker(this);
            // Add current optimum to seed
            if (nonInteractive) {
                double[][] sample = new double[1][];
                sample[0] = currentOptimum.getParameters();
                seed = merge(seed, sample);
            }
            ss.seed(seed);
            ConvergenceChecker<FilterScore> checker = new InterruptConvergenceChecker(0, 0, maxIterations, convergedCount);
            createGAWindow();
            resumeFilterTimer();
            SearchResult<FilterScore> optimum = ss.enrichmentSearch(dimensions, this, checker, enrichmentSamples, enrichmentFraction, enrichmentPadding);
            if (optimum != null) {
                // In case optimisation was stopped
                IJ.resetEscape();
                best = ((SimpleFilterScore) optimum.score).r.filter;
                // Not required as the search now respects the min interval
                // Enumerate on the min interval to produce the final filter
                //best = enumerateMinInterval(best, enabled);
                // Now update the filter set for final assessment
                filterSet = new FilterSet(filterSet.getName(), searchSpaceToFilters((DirectFilter) best, ss.getSearchSpace()));
                // Option to save the filters
                if (saveOption)
                    saveFilterSet(filterSet, setNumber, !nonInteractive);
            }
        } else
            resumeFilterTimer();
    }
    IJ.showStatus("Analysing [" + setNumber + "] " + filterSet.getName() + " ...");
    // Do not support plotting if we used optimisation
    double[] xValues = (best != null || isHeadless || (plotTopN == 0)) ? null : new double[filterSet.size()];
    double[] yValues = (xValues == null) ? null : new double[xValues.length];
    SimpleFilterScore max = null;
    // It can just assess the top 1 required for the summary.
    if (best != null) {
        // Only assess the top 1 filter for the summary
        List<Filter> list = new ArrayList<Filter>();
        list.add((DirectFilter) best);
        filterSet = new FilterSet(filterSet.getName(), list);
    }
    // Score the filters and report the results if configured.
    FilterScoreResult[] scoreResults = scoreFilters(setUncomputedStrength(filterSet), showResultsTable);
    if (scoreResults == null)
        return -1;
    analysisStopWatch.stop();
    for (int index = 0; index < scoreResults.length; index++) {
        final FilterScoreResult scoreResult = scoreResults[index];
        if (xValues != null) {
            xValues[index] = scoreResult.filter.getNumericalValue();
            yValues[index] = scoreResult.score;
        }
        final SimpleFilterScore result = new SimpleFilterScore(scoreResult, allSameType, scoreResult.criteria >= minCriteria);
        if (result.compareTo(max) < 0) {
            max = result;
        }
    }
    if (showResultsTable) {
        BufferedTextWindow tw = null;
        if (resultsWindow != null) {
            tw = new BufferedTextWindow(resultsWindow);
            tw.setIncrement(Integer.MAX_VALUE);
        }
        for (int index = 0; index < scoreResults.length; index++) addToResultsWindow(tw, scoreResults[index].text);
        if (resultsWindow != null)
            resultsWindow.getTextPanel().updateDisplay();
    }
    // Check the top filter against the limits of the original dimensions
    char[] atLimit = null;
    if (allSameType && originalDimensions != null) {
        DirectFilter filter = max.r.filter;
        int[] indices = filter.getChromosomeParameters();
        atLimit = new char[indices.length];
        StringBuilder sb = new StringBuilder(200);
        for (int j = 0; j < indices.length; j++) {
            atLimit[j] = ComplexFilterScore.WITHIN;
            final int p = indices[j];
            if (disabled[p])
                continue;
            final double value = filter.getParameterValue(p);
            double lowerLimit = originalDimensions[p].getLower();
            double upperLimit = originalDimensions[p].getUpper();
            int c1 = Double.compare(value, lowerLimit);
            if (c1 <= 0) {
                atLimit[j] = ComplexFilterScore.FLOOR;
                sb.append(" : ").append(filter.getParameterName(p)).append(' ').append(atLimit[j]).append('[').append(Utils.rounded(value));
                if (c1 == -1) {
                    atLimit[j] = ComplexFilterScore.BELOW;
                    sb.append("<").append(Utils.rounded(lowerLimit));
                }
                sb.append("]");
            } else {
                int c2 = Double.compare(value, upperLimit);
                if (c2 >= 0) {
                    atLimit[j] = ComplexFilterScore.CEIL;
                    sb.append(" : ").append(filter.getParameterName(p)).append(' ').append(atLimit[j]).append('[').append(Utils.rounded(value));
                    if (c2 == 1) {
                        atLimit[j] = ComplexFilterScore.ABOVE;
                        sb.append(">").append(Utils.rounded(upperLimit));
                    }
                    sb.append("]");
                }
            }
        }
        if (sb.length() > 0) {
            if (max.criteriaPassed) {
                Utils.log("Warning: Top filter (%s @ %s|%s) [%s] at the limit of the expanded range%s", filter.getName(), Utils.rounded((invertScore) ? -max.score : max.score), Utils.rounded((invertCriteria) ? -minCriteria : minCriteria), limitFailCount + limitRange, sb.toString());
            } else {
                Utils.log("Warning: Top filter (%s @ -|%s) [%s] at the limit of the expanded range%s", filter.getName(), Utils.rounded((invertCriteria) ? -max.criteria : max.criteria), limitFailCount + limitRange, sb.toString());
            }
        }
    }
    // Note that max should never be null since this method is not run with an empty filter set
    // We may have no filters that pass the criteria
    String type = max.r.filter.getType();
    if (!max.criteriaPassed) {
        Utils.log("Warning: Filter does not pass the criteria: %s : Best = %s using %s", type, Utils.rounded((invertCriteria) ? -max.criteria : max.criteria), max.r.filter.getName());
        return 0;
    }
    // This could be an option?
    boolean allowDuplicates = true;
    // XXX - Commented out the requirement to be the same type to store for later analysis. 
    // This may break the code, however I think that all filter sets should be able to have a best filter
    // irrespective of whether they were the same type or not.
    //if (allSameType)
    //{
    ComplexFilterScore newFilterScore = new ComplexFilterScore(max.r, atLimit, algorithm, analysisStopWatch.getTime(), "", 0);
    addBestFilter(type, allowDuplicates, newFilterScore);
    // Add spacer at end of each result set
    if (isHeadless) {
        if (showResultsTable && filterSet.size() > 1)
            IJ.log("");
    } else {
        if (showResultsTable && filterSet.size() > 1)
            resultsWindow.append("");
        if (plotTopN > 0 && xValues != null) {
            // Check the xValues are unique. Since the filters have been sorted by their
            // numeric value we only need to compare adjacent entries.
            boolean unique = true;
            for (int ii = 0; ii < xValues.length - 1; ii++) {
                if (xValues[ii] == xValues[ii + 1]) {
                    unique = false;
                    break;
                }
            }
            String xAxisName = filterSet.getValueName();
            if (unique) {
                // Check the values all refer to the same property
                for (Filter filter : filterSet.getFilters()) {
                    if (!xAxisName.equals(filter.getNumericalValueName())) {
                        unique = false;
                        break;
                    }
                }
            }
            if (!unique) {
                // If not unique then renumber them and use an arbitrary label
                xAxisName = "Filter";
                for (int ii = 0; ii < xValues.length; ii++) xValues[ii] = ii + 1;
            }
            String title = filterSet.getName();
            // Check if a previous filter set had the same name, update if necessary
            NamedPlot p = getNamedPlot(title);
            if (p == null)
                plots.add(new NamedPlot(title, xAxisName, xValues, yValues));
            else
                p.updateValues(xAxisName, xValues, yValues);
            if (plots.size() > plotTopN) {
                Collections.sort(plots);
                p = plots.remove(plots.size() - 1);
            }
        }
    }
    return 0;
}

Example 3 with Evaluation

use of org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation in project GDSC-SMLM by aherbert.

the class PoissonGammaGaussianFunction method likelihood.

/**
	 * Compute the likelihood
	 * 
	 * @param o
	 *            The observed count
	 * @param e
	 *            The expected count
	 * @return The likelihood
	 */
public double likelihood(final double o, final double e) {
    // Use the same variables as the Python code
    final double cij = o;
    // convert to photons
    final double eta = alpha * e;
    if (sigma == 0) {
        // No convolution with a Gaussian. Simply evaluate for a Poisson-Gamma distribution.
        final double p;
        // Any observed count above zero
        if (cij > 0.0) {
            // The observed count converted to photons
            final double nij = alpha * cij;
            // The limit on eta * nij is therefore (709/2)^2 = 125670.25
            if (eta * nij > 10000) {
                // Approximate Bessel function i1(x) when using large x:
                // i1(x) ~ exp(x)/sqrt(2*pi*x)
                // However the entire equation is logged (creating transform),
                // evaluated then raised to e to prevent overflow error on 
                // large exp(x)
                final double transform = 0.5 * Math.log(alpha * eta / cij) - nij - eta + 2 * Math.sqrt(eta * nij) - Math.log(twoSqrtPi * Math.pow(eta * nij, 0.25));
                p = FastMath.exp(transform);
            } else {
                // Second part of equation 135
                p = Math.sqrt(alpha * eta / cij) * FastMath.exp(-nij - eta) * Bessel.I1(2 * Math.sqrt(eta * nij));
            }
        } else if (cij == 0.0) {
            p = FastMath.exp(-eta);
        } else {
            p = 0;
        }
        return (p > minimumProbability) ? p : minimumProbability;
    } else if (useApproximation) {
        return mortensenApproximation(cij, eta);
    } else {
        // This code is the full evaluation of equation 7 from the supplementary information  
        // of the paper Chao, et al (2013) Nature Methods 10, 335-338.
        // It is the full evaluation of a Poisson-Gamma-Gaussian convolution PMF. 
        // Read noise
        final double sk = sigma;
        // Gain
        final double g = 1.0 / alpha;
        // Observed pixel value
        final double z = o;
        // Expected number of photons
        final double vk = eta;
        // Compute the integral to infinity of:
        // exp( -((z-u)/(sqrt(2)*s)).^2 - u/g ) * sqrt(vk*u/g) .* besseli(1, 2 * sqrt(vk*u/g)) ./ u;
        // vk / g
        final double vk_g = vk * alpha;
        final double sqrt2sigma = Math.sqrt(2) * sk;
        // Specify the function to integrate
        UnivariateFunction f = new UnivariateFunction() {

            public double value(double u) {
                return eval(sqrt2sigma, z, vk_g, g, u);
            }
        };
        // Integrate to infinity is not necessary. The convolution of the function with the 
        // Gaussian should be adequately sampled using a nxSD around the maximum.
        // Find a bracket containing the maximum
        double lower, upper;
        double maxU = Math.max(1, o);
        double rLower = maxU;
        double rUpper = maxU + 1;
        double f1 = f.value(rLower);
        double f2 = f.value(rUpper);
        // Calculate the simple integral and the range
        double sum = f1 + f2;
        boolean searchUp = f2 > f1;
        if (searchUp) {
            while (f2 > f1) {
                f1 = f2;
                rUpper += 1;
                f2 = f.value(rUpper);
                sum += f2;
            }
            maxU = rUpper - 1;
        } else {
            // Ensure that u stays above zero
            while (f1 > f2 && rLower > 1) {
                f2 = f1;
                rLower -= 1;
                f1 = f.value(rLower);
                sum += f1;
            }
            maxU = (rLower > 1) ? rLower + 1 : rLower;
        }
        lower = Math.max(0, maxU - 5 * sk);
        upper = maxU + 5 * sk;
        if (useSimpleIntegration && lower > 0) {
            // remaining points in the range
            for (double u = rLower - 1; u >= lower; u -= 1) {
                sum += f.value(u);
            }
            for (double u = rUpper + 1; u <= upper; u += 1) {
                sum += f.value(u);
            }
        } else {
            // Use Legendre-Gauss integrator
            try {
                final double relativeAccuracy = 1e-4;
                final double absoluteAccuracy = 1e-8;
                final int minimalIterationCount = 3;
                final int maximalIterationCount = 32;
                final int integrationPoints = 16;
                // Use an integrator that does not use the boundary since u=0 is undefined (divide by zero)
                UnivariateIntegrator i = new IterativeLegendreGaussIntegrator(integrationPoints, relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
                sum = i.integrate(2000, f, lower, upper);
            } catch (TooManyEvaluationsException ex) {
                return mortensenApproximation(cij, eta);
            }
        }
        // Compute the final probability
        //final double 
        f1 = z / sqrt2sigma;
        final double p = (FastMath.exp(-vk) / (sqrt2pi * sk)) * (FastMath.exp(-(f1 * f1)) + sum);
        return (p > minimumProbability) ? p : minimumProbability;
    }
}

Example 4 with Evaluation

use of org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation 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.
 */
@Nullable
private double[][] doFitJumpDistancesMle(double[] jumpDistances, double estimatedD, int n) {
    final MaxEval maxEval = new MaxEval(20000);
    final 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
            final SimpleBounds bounds = new SimpleBounds(function.getLowerBounds(), function.getUpperBounds(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY));
            final PointValuePair solution = powellOptimizer.optimize(maxEval, new ObjectiveFunction(function), new InitialGuess(function.guess()), bounds, new CustomPowellOptimizer.BasisStep(function.step()), GoalType.MAXIMIZE);
            final double[] fitParams = solution.getPointRef();
            ll = solution.getValue();
            lastFitValue = fitValue = MathUtils.getAkaikeInformationCriterion(ll, 1);
            final double[] coefficients = fitParams;
            final double[] fractions = new double[] { 1 };
            LoggerUtils.log(logger, Level.INFO, "Fit Jump distance (N=1) : %s, MLE = %s, Akaike IC = %s (%d evaluations)", formatD(fitParams[0]), MathUtils.rounded(ll, 4), MathUtils.rounded(fitValue, 4), powellOptimizer.getEvaluations());
            return new double[][] { coefficients, fractions };
        } catch (final TooManyEvaluationsException ex) {
            LoggerUtils.log(logger, Level.INFO, "Powell optimiser failed to fit (N=1) : Too many evaluation (%d)", powellOptimizer.getEvaluations());
        } catch (final TooManyIterationsException ex) {
            LoggerUtils.log(logger, Level.INFO, "Powell optimiser failed to fit (N=1) : Too many iterations (%d)", powellOptimizer.getIterations());
        } catch (final ConvergenceException ex) {
            LoggerUtils.log(logger, Level.INFO, "Powell optimiser failed to fit (N=1) : %s", ex.getMessage());
        }
        return null;
    }
    final MixedJumpDistanceFunction function = new MixedJumpDistanceFunction(jumpDistances, estimatedD, n);
    final 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();
        LoggerUtils.log(logger, Level.FINE, "Powell optimiser fit (N=%d) : MLE = %f (%d evaluations)", n, constrainedSolution.getValue(), powellOptimizer.getEvaluations());
    } catch (final TooManyEvaluationsException ex) {
        LoggerUtils.log(logger, Level.INFO, "Powell optimiser failed to fit (N=%d) : Too many evaluation (%d)", n, powellOptimizer.getEvaluations());
    } catch (final TooManyIterationsException ex) {
        LoggerUtils.log(logger, Level.INFO, "Powell optimiser failed to fit (N=%d) : Too many iterations (%d)", n, powellOptimizer.getIterations());
    } catch (final ConvergenceException ex) {
        LoggerUtils.log(logger, Level.INFO, "Powell optimiser failed to fit (N=%d) : %s", n, ex.getMessage());
    }
    if (constrainedSolution == null) {
        LoggerUtils.log(logger, Level.INFO, "Trying CMAES optimiser with restarts ...");
        final double[] uB = function.getUpperBounds();
        final 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.
        final double[] s = new double[lB.length];
        for (int i = 0; i < s.length; i++) {
            s[i] = (uB[i] - lB[i]) / 3;
        }
        final OptimizationData sigma = new CMAESOptimizer.Sigma(s);
        final OptimizationData popSize = new CMAESOptimizer.PopulationSize((int) (4 + Math.floor(3 * Math.log(function.x.length))));
        // Iterate this for stability in the initial guess
        final CMAESOptimizer cmaesOptimizer = createCmaesOptimizer();
        for (int i = 0; i <= fitRestarts; i++) {
            // Try from the initial guess
            try {
                final 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;
                    LoggerUtils.log(logger, Level.FINE, "CMAES optimiser [%da] fit (N=%d) : MLE = %f (%d evaluations)", i, n, solution.getValue(), evaluations);
                }
            } catch (final TooManyEvaluationsException ex) {
            // No solution
            }
            if (constrainedSolution == null) {
                continue;
            }
            // Try from the current optimum
            try {
                final 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;
                    LoggerUtils.log(logger, Level.FINE, "CMAES optimiser [%db] fit (N=%d) : MLE = %f (%d evaluations)", i, n, solution.getValue(), evaluations);
                }
            } catch (final TooManyEvaluationsException ex) {
            // No solution
            }
        }
        if (constrainedSolution != null) {
            try {
                // Re-optimise with Powell?
                final 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;
                    LoggerUtils.log(logger, Level.INFO, "Powell optimiser re-fit (N=%d) : MLE = %f (%d evaluations)", n, constrainedSolution.getValue(), powellOptimizer.getEvaluations());
                }
            } catch (final TooManyEvaluationsException ex) {
            // No solution
            } catch (final TooManyIterationsException ex) {
            // No solution
            } catch (final ConvergenceException ex) {
            // No solution
            }
        }
    }
    if (constrainedSolution == null) {
        LoggerUtils.log(logger, Level.INFO, "Failed to fit N=%d", n);
        return null;
    }
    final 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
    fitValue = MathUtils.getAkaikeInformationCriterion(ll, fitParams.length - 1);
    final double[] d = new double[n];
    final 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);
    final double[] coefficients = d;
    final double[] fractions = f;
    LoggerUtils.log(logger, Level.INFO, "Fit Jump distance (N=%d) : %s (%s), MLE = %s, Akaike IC = %s (%d evaluations)", n, formatD(d), format(f), MathUtils.rounded(ll, 4), MathUtils.rounded(fitValue, 4), evaluations);
    if (isValid(d, f)) {
        lastFitValue = fitValue;
        return new double[][] { coefficients, fractions };
    }
    return null;
}

Example 5 with Evaluation

use of org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation in project GDSC-SMLM by aherbert.

the class CameraModelAnalysis method execute.

/**
 * Execute the analysis.
 */
private boolean execute() {
    dirty = false;
    final CameraModelAnalysisSettings settings = this.settings.build();
    if (!(getGain(settings) > 0)) {
        ImageJUtils.log(TITLE + "Error: No total gain");
        return false;
    }
    if (!(settings.getPhotons() > 0)) {
        ImageJUtils.log(TITLE + "Error: No photons");
        return false;
    }
    // Avoid repeating the same analysis
    if (settings.equals(lastSettings)) {
        return true;
    }
    lastSettings = settings;
    final IntHistogram h = getHistogram(settings);
    // Build cumulative distribution
    final double[][] cdf1 = cumulativeHistogram(h);
    final double[] x1 = cdf1[0];
    final double[] y1 = cdf1[1];
    // Interpolate to 300 steps faster evaluation?
    // Get likelihood function
    final LikelihoodFunction f = getLikelihoodFunction(settings);
    // Create likelihood cumulative distribution
    final double[][] cdf2 = cumulativeDistribution(settings, cdf1, f);
    // Compute Komolgorov distance
    final double[] distanceAndValue = getDistance(cdf1, cdf2);
    final double distance = distanceAndValue[0];
    final double value = distanceAndValue[1];
    final double area = distanceAndValue[2];
    final double[] x2 = cdf2[0];
    final double[] y2 = cdf2[1];
    // Fill y1
    int offset = 0;
    while (x2[offset] < x1[0]) {
        offset++;
    }
    final double[] y1b = new double[y2.length];
    System.arraycopy(y1, 0, y1b, offset, y1.length);
    Arrays.fill(y1b, offset + y1.length, y2.length, y1[y1.length - 1]);
    // KolmogorovSmirnovTest
    // n is the number of samples used to build the probability distribution.
    final int n = (int) MathUtils.sum(h.histogramCounts);
    // From KolmogorovSmirnovTest.kolmogorovSmirnovTest(RealDistribution distribution, double[]
    // data, boolean exact):
    // Returns the p-value associated with the null hypothesis that data is a sample from
    // distribution.
    // E.g. If p<0.05 then the null hypothesis is rejected and the data do not match the
    // distribution.
    double pvalue = Double.NaN;
    try {
        pvalue = 1d - kolmogorovSmirnovTest.cdf(distance, n);
    } catch (final MathArithmeticException ex) {
    // Cannot be computed to leave at NaN
    }
    // Plot
    final WindowOrganiser wo = new WindowOrganiser();
    String title = TITLE + " CDF";
    Plot plot = new Plot(title, "Count", "CDF");
    final double max = 1.05 * MathUtils.maxDefault(1, y2);
    plot.setLimits(x2[0], x2[x2.length - 1], 0, max);
    plot.setColor(Color.blue);
    plot.addPoints(x2, y1b, Plot.BAR);
    plot.setColor(Color.red);
    plot.addPoints(x2, y2, Plot.BAR);
    plot.setColor(Color.magenta);
    plot.drawLine(value, 0, value, max);
    plot.setColor(Color.black);
    plot.addLegend("CDF\nModel");
    plot.addLabel(0, 0, String.format("Distance=%s @ %.0f (Mean=%s) : p=%s", MathUtils.rounded(distance), value, MathUtils.rounded(area), MathUtils.rounded(pvalue)));
    ImageJUtils.display(title, plot, ImageJUtils.NO_TO_FRONT, wo);
    // Show the histogram
    title = TITLE + " Histogram";
    plot = new Plot(title, "Count", "Frequency");
    plot.setLimits(x1[0] - 0.5, x1[x1.length - 1] + 1.5, 0, MathUtils.max(h.histogramCounts) * 1.05);
    plot.setColor(Color.blue);
    plot.addPoints(x1, SimpleArrayUtils.toDouble(h.histogramCounts), Plot.BAR);
    plot.setColor(Color.red);
    final double[] x = floatHistogram[0].clone();
    final double[] y = floatHistogram[1].clone();
    final double scale = n / (MathUtils.sum(y) * (x[1] - x[0]));
    for (int i = 0; i < y.length; i++) {
        y[i] *= scale;
    }
    plot.addPoints(x, y, Plot.LINE);
    plot.setColor(Color.black);
    plot.addLegend("Sample\nExpected");
    ImageJUtils.display(title, plot, ImageJUtils.NO_TO_FRONT, wo);
    wo.tile();
    return true;
}