Search in sources :

Example 11 with MaxEval

use of org.apache.commons.math3.optim.MaxEval in project gatk-protected by broadinstitute.

the class CNLOHCaller method optimizeIt.

private double optimizeIt(final Function<Double, Double> objectiveFxn, final SearchInterval searchInterval) {
    final MaxEval BRENT_MAX_EVAL = new MaxEval(1000);
    final double RELATIVE_TOLERANCE = 0.001;
    final double ABSOLUTE_TOLERANCE = 0.001;
    final BrentOptimizer OPTIMIZER = new BrentOptimizer(RELATIVE_TOLERANCE, ABSOLUTE_TOLERANCE);
    final UnivariateObjectiveFunction objective = new UnivariateObjectiveFunction(x -> objectiveFxn.apply(x));
    return OPTIMIZER.optimize(objective, GoalType.MAXIMIZE, searchInterval, BRENT_MAX_EVAL).getPoint();
}
Also used : MaxEval(org.apache.commons.math3.optim.MaxEval) UnivariateObjectiveFunction(org.apache.commons.math3.optim.univariate.UnivariateObjectiveFunction) BrentOptimizer(org.apache.commons.math3.optim.univariate.BrentOptimizer)

Example 12 with MaxEval

use of org.apache.commons.math3.optim.MaxEval in project vcell by virtualcell.

the class FitTimeSeries method fitToGaussian.

static GaussianFitResults fitToGaussian(double init_center_i, double init_center_j, double init_radius2, FloatImage image) {
    // 
    // do some optimization on the image (fitting to a Gaussian)
    // set initial guesses from ROI operation.
    // 
    ISize imageSize = image.getISize();
    final int num_i = imageSize.getX();
    final int num_j = imageSize.getY();
    final float[] floatPixels = image.getFloatPixels();
    // 
    // initial guess based on previous fit of ROI
    // do gaussian fit in index space for center and standard deviation (later to translate it back to world coordinates)
    // 
    final int window_size = (int) Math.sqrt(init_radius2) * 4;
    // final int window_min_i = 0;       // (int) Math.max(0, Math.floor(init_center_i - window_size/2));
    // final int window_max_i = num_i-1; // (int) Math.min(num_i-1, Math.ceil(init_center_i + window_size/2));
    // final int window_min_j = 0;       // (int) Math.max(0, Math.floor(init_center_j - window_size/2));
    // final int window_max_j = num_j-1; // (int) Math.min(num_j-1, Math.ceil(init_center_j + window_size/2));
    final int window_min_i = (int) Math.max(0, Math.floor(init_center_i - window_size / 2));
    final int window_max_i = (int) Math.min(num_i - 1, Math.ceil(init_center_i + window_size / 2));
    final int window_min_j = (int) Math.max(0, Math.floor(init_center_j - window_size / 2));
    final int window_max_j = (int) Math.min(num_j - 1, Math.ceil(init_center_j + window_size / 2));
    final int PARAM_INDEX_CENTER_I = 0;
    final int PARAM_INDEX_CENTER_J = 1;
    final int PARAM_INDEX_K = 2;
    final int PARAM_INDEX_HIGH = 3;
    final int PARAM_INDEX_RADIUS_SQUARED = 4;
    final int NUM_PARAMETERS = 5;
    double[] initParameters = new double[NUM_PARAMETERS];
    initParameters[PARAM_INDEX_CENTER_I] = init_center_i;
    initParameters[PARAM_INDEX_CENTER_J] = init_center_j;
    initParameters[PARAM_INDEX_HIGH] = 1.0;
    initParameters[PARAM_INDEX_K] = 10;
    initParameters[PARAM_INDEX_RADIUS_SQUARED] = init_radius2;
    PowellOptimizer optimizer = new PowellOptimizer(1e-4, 1e-1);
    MultivariateFunction func = new MultivariateFunction() {

        @Override
        public double value(double[] point) {
            double center_i = point[PARAM_INDEX_CENTER_I];
            double center_j = point[PARAM_INDEX_CENTER_J];
            double high = point[PARAM_INDEX_HIGH];
            double K = point[PARAM_INDEX_K];
            double radius2 = point[PARAM_INDEX_RADIUS_SQUARED];
            double error2 = 0;
            for (int j = window_min_j; j <= window_max_j; j++) {
                // double y = j - center_j;
                double y = j;
                for (int i = window_min_i; i <= window_max_i; i++) {
                    // double x = i - center_i;
                    double x = i;
                    double modelValue = high - FastMath.exp(-K * FastMath.exp(-2 * (x * x + y * y) / radius2));
                    double imageValue = floatPixels[j * num_i + i];
                    double error = modelValue - imageValue;
                    error2 += error * error;
                }
            }
            System.out.println(new GaussianFitResults(center_i, center_j, radius2, K, high, error2));
            return error2;
        }
    };
    PointValuePair pvp = optimizer.optimize(new ObjectiveFunction(func), new InitialGuess(initParameters), new MaxEval(100000), GoalType.MINIMIZE);
    double[] fittedParamValues = pvp.getPoint();
    double fitted_center_i = fittedParamValues[PARAM_INDEX_CENTER_I];
    double fitted_center_j = fittedParamValues[PARAM_INDEX_CENTER_J];
    double fitted_radius2 = fittedParamValues[PARAM_INDEX_RADIUS_SQUARED];
    double fitted_K = fittedParamValues[PARAM_INDEX_K];
    double fitted_high = fittedParamValues[PARAM_INDEX_HIGH];
    double objectiveFunctionValue = pvp.getValue();
    return new GaussianFitResults(fitted_center_i, fitted_center_j, fitted_radius2, fitted_K, fitted_high, objectiveFunctionValue);
}
Also used : MultivariateFunction(org.apache.commons.math3.analysis.MultivariateFunction) InitialGuess(org.apache.commons.math3.optim.InitialGuess) MaxEval(org.apache.commons.math3.optim.MaxEval) ISize(org.vcell.util.ISize) ObjectiveFunction(org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunction) PowellOptimizer(org.apache.commons.math3.optim.nonlinear.scalar.noderiv.PowellOptimizer) PointValuePair(org.apache.commons.math3.optim.PointValuePair)

Example 13 with MaxEval

use of org.apache.commons.math3.optim.MaxEval in project tetrad by cmu-phil.

the class GeneralizedSemIm method simulateDataMinimizeSurface.

public DataSet simulateDataMinimizeSurface(int sampleSize, boolean latentDataSaved) {
    final Map<String, Double> variableValues = new HashMap<>();
    List<Node> continuousVariables = new LinkedList<>();
    final List<Node> variableNodes = pm.getVariableNodes();
    // Work with a copy of the variables, because their type can be set externally.
    for (Node node : variableNodes) {
        ContinuousVariable var = new ContinuousVariable(node.getName());
        var.setNodeType(node.getNodeType());
        if (var.getNodeType() != NodeType.ERROR) {
            continuousVariables.add(var);
        }
    }
    DataSet fullDataSet = new ColtDataSet(sampleSize, continuousVariables);
    final Context context = new Context() {

        public Double getValue(String term) {
            Double value = parameterValues.get(term);
            if (value != null) {
                return value;
            }
            value = variableValues.get(term);
            if (value != null) {
                return value;
            }
            throw new IllegalArgumentException("No value recorded for '" + term + "'");
        }
    };
    final double[] _metric = new double[1];
    MultivariateFunction function = new MultivariateFunction() {

        double metric;

        public double value(double[] doubles) {
            for (int i = 0; i < variableNodes.size(); i++) {
                variableValues.put(variableNodes.get(i).getName(), doubles[i]);
            }
            double[] image = new double[doubles.length];
            for (int i = 0; i < variableNodes.size(); i++) {
                Node node = variableNodes.get(i);
                Expression expression = pm.getNodeExpression(node);
                image[i] = expression.evaluate(context);
                if (Double.isNaN(image[i])) {
                    throw new IllegalArgumentException("Undefined value for expression " + expression);
                }
            }
            metric = 0.0;
            for (int i = 0; i < variableNodes.size(); i++) {
                double diff = doubles[i] - image[i];
                metric += diff * diff;
            }
            for (int i = 0; i < variableNodes.size(); i++) {
                variableValues.put(variableNodes.get(i).getName(), image[i]);
            }
            _metric[0] = metric;
            return metric;
        }
    };
    MultivariateOptimizer search = new PowellOptimizer(1e-7, 1e-7);
    // Do the simulation.
    ROW: for (int row = 0; row < sampleSize; row++) {
        // Take random draws from error distributions.
        for (Node variable : variableNodes) {
            Node error = pm.getErrorNode(variable);
            if (error == null) {
                throw new NullPointerException();
            }
            Expression expression = pm.getNodeExpression(error);
            double value = expression.evaluate(context);
            if (Double.isNaN(value)) {
                throw new IllegalArgumentException("Undefined value for expression: " + expression);
            }
            variableValues.put(error.getName(), value);
        }
        for (Node variable : variableNodes) {
            // RandomUtil.getInstance().nextUniform(-5, 5));
            variableValues.put(variable.getName(), 0.0);
        }
        while (true) {
            double[] values = new double[variableNodes.size()];
            for (int i = 0; i < values.length; i++) {
                values[i] = variableValues.get(variableNodes.get(i).getName());
            }
            PointValuePair pair = search.optimize(new InitialGuess(values), new ObjectiveFunction(function), GoalType.MINIMIZE, new MaxEval(100000));
            values = pair.getPoint();
            for (int i = 0; i < variableNodes.size(); i++) {
                if (isSimulatePositiveDataOnly() && values[i] < 0) {
                    row--;
                    continue ROW;
                }
                if (!Double.isNaN(selfLoopCoef) && row > 0) {
                    values[i] += selfLoopCoef * fullDataSet.getDouble(row - 1, i);
                }
                variableValues.put(variableNodes.get(i).getName(), values[i]);
                fullDataSet.setDouble(row, i, values[i]);
            }
            if (_metric[0] < 0.01) {
                // while
                break;
            }
        }
    }
    if (latentDataSaved) {
        return fullDataSet;
    } else {
        return DataUtils.restrictToMeasured(fullDataSet);
    }
}
Also used : Context(edu.cmu.tetrad.calculator.expression.Context) MultivariateOptimizer(org.apache.commons.math3.optim.nonlinear.scalar.MultivariateOptimizer) InitialGuess(org.apache.commons.math3.optim.InitialGuess) MaxEval(org.apache.commons.math3.optim.MaxEval) ObjectiveFunction(org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunction) PowellOptimizer(org.apache.commons.math3.optim.nonlinear.scalar.noderiv.PowellOptimizer) PointValuePair(org.apache.commons.math3.optim.PointValuePair) MultivariateFunction(org.apache.commons.math3.analysis.MultivariateFunction) Expression(edu.cmu.tetrad.calculator.expression.Expression)

Example 14 with MaxEval

use of org.apache.commons.math3.optim.MaxEval in project tetrad by cmu-phil.

the class SemOptimizerPowell method optimize.

// =========================PUBLIC METHODS==========================//
public void optimize(SemIm semIm) {
    double min = Double.POSITIVE_INFINITY;
    double[] point = null;
    for (int count = 0; count < numRestarts + 1; count++) {
        // System.out.println("Trial " + (count + 1));
        SemIm _sem2 = new SemIm(semIm);
        List<Parameter> freeParameters = _sem2.getFreeParameters();
        double[] p = new double[freeParameters.size()];
        for (int i = 0; i < freeParameters.size(); i++) {
            if (freeParameters.get(i).getType() == ParamType.VAR) {
                p[i] = RandomUtil.getInstance().nextUniform(0, 1);
            } else {
                p[i] = RandomUtil.getInstance().nextUniform(-1, 1);
            }
        }
        _sem2.setFreeParamValues(p);
        MultivariateOptimizer search = new PowellOptimizer(1e-7, 1e-7);
        PointValuePair pair = search.optimize(new InitialGuess(_sem2.getFreeParamValues()), new ObjectiveFunction(fittingFunction(semIm)), GoalType.MINIMIZE, new MaxEval(100000));
        double chisq = _sem2.getChiSquare();
        if (chisq < min) {
            min = chisq;
            point = pair.getPoint();
        }
    }
    if (point == null) {
        throw new NullPointerException("Point could not be found.");
    }
    System.arraycopy(point, 0, semIm.getFreeParamValues(), 0, point.length);
}
Also used : MultivariateOptimizer(org.apache.commons.math3.optim.nonlinear.scalar.MultivariateOptimizer) InitialGuess(org.apache.commons.math3.optim.InitialGuess) MaxEval(org.apache.commons.math3.optim.MaxEval) ObjectiveFunction(org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunction) PowellOptimizer(org.apache.commons.math3.optim.nonlinear.scalar.noderiv.PowellOptimizer) PointValuePair(org.apache.commons.math3.optim.PointValuePair)

Example 15 with MaxEval

use of org.apache.commons.math3.optim.MaxEval in project imagingbook-common by imagingbook.

the class FourierDescriptor method getStartPointPhase.

/**
 * Calculates the 'canonical' start point. This version uses
 * (a) a coarse search for a global maximum of fp() and subsequently
 * (b) a numerical optimization using Brent's method
 * (implemented with Apache Commons Math).
 *
 * @param Mp number of Fourier coefficient pairs
 * @return start point phase
 */
public double getStartPointPhase(int Mp) {
    Mp = Math.min(Mp, (G.length - 1) / 2);
    UnivariateFunction fp = new TargetFunction(Mp);
    // search for the global maximum in coarse steps
    double cmax = Double.NEGATIVE_INFINITY;
    int kmax = -1;
    // number of steps over 180 degrees
    int K = 25;
    for (int k = 0; k < K; k++) {
        // phase to evaluate
        final double phi = Math.PI * k / K;
        final double c = fp.value(phi);
        if (c > cmax) {
            cmax = c;
            kmax = k;
        }
    }
    // optimize using previous and next point as the bracket.
    double minPhi = Math.PI * (kmax - 1) / K;
    double maxPhi = Math.PI * (kmax + 1) / K;
    UnivariateOptimizer optimizer = new BrentOptimizer(1E-4, 1E-6);
    int maxIter = 20;
    UnivariatePointValuePair result = optimizer.optimize(new MaxEval(maxIter), new UnivariateObjectiveFunction(fp), GoalType.MAXIMIZE, new SearchInterval(minPhi, maxPhi));
    double phi0 = result.getPoint();
    // the canonical start point phase
    return phi0;
}
Also used : MaxEval(org.apache.commons.math3.optim.MaxEval) SearchInterval(org.apache.commons.math3.optim.univariate.SearchInterval) UnivariateFunction(org.apache.commons.math3.analysis.UnivariateFunction) UnivariateObjectiveFunction(org.apache.commons.math3.optim.univariate.UnivariateObjectiveFunction) BrentOptimizer(org.apache.commons.math3.optim.univariate.BrentOptimizer) UnivariateOptimizer(org.apache.commons.math3.optim.univariate.UnivariateOptimizer) UnivariatePointValuePair(org.apache.commons.math3.optim.univariate.UnivariatePointValuePair)

Aggregations

MaxEval (org.apache.commons.math3.optim.MaxEval)47 InitialGuess (org.apache.commons.math3.optim.InitialGuess)39 PointValuePair (org.apache.commons.math3.optim.PointValuePair)39 ObjectiveFunction (org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunction)39 TooManyEvaluationsException (org.apache.commons.math3.exception.TooManyEvaluationsException)19 MultivariateOptimizer (org.apache.commons.math3.optim.nonlinear.scalar.MultivariateOptimizer)16 PowellOptimizer (org.apache.commons.math3.optim.nonlinear.scalar.noderiv.PowellOptimizer)15 SimpleBounds (org.apache.commons.math3.optim.SimpleBounds)14 UnivariateObjectiveFunction (org.apache.commons.math3.optim.univariate.UnivariateObjectiveFunction)12 MultivariateFunction (org.apache.commons.math3.analysis.MultivariateFunction)11 TooManyIterationsException (org.apache.commons.math3.exception.TooManyIterationsException)10 OptimizationData (org.apache.commons.math3.optim.OptimizationData)10 SimpleValueChecker (org.apache.commons.math3.optim.SimpleValueChecker)10 CMAESOptimizer (org.apache.commons.math3.optim.nonlinear.scalar.noderiv.CMAESOptimizer)10 UnivariatePointValuePair (org.apache.commons.math3.optim.univariate.UnivariatePointValuePair)10 ConvergenceException (org.apache.commons.math3.exception.ConvergenceException)8 ObjectiveFunctionGradient (org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunctionGradient)6 NelderMeadSimplex (org.apache.commons.math3.optim.nonlinear.scalar.noderiv.NelderMeadSimplex)6 BrentOptimizer (org.apache.commons.math3.optim.univariate.BrentOptimizer)6 SearchInterval (org.apache.commons.math3.optim.univariate.SearchInterval)6