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();
}
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);
}
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);
}
}
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);
}
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;
}
Aggregations