List of usage examples for org.apache.commons.math.optimization.direct NelderMead optimize
public RealPointValuePair optimize(final MultivariateRealFunction function, final GoalType goalType, final double[] startPoint) throws FunctionEvaluationException, OptimizationException, IllegalArgumentException
From source file:ch.algotrader.option.SABR.java
/** * Perfors a SABR calibartion based on specified volatilities. * * @return SABRSmileVO The SABR smile// ww w . ja v a 2 s . c om */ public static SABRSmileVO calibrate(final Double[] strikes, final Double[] volatilities, final double atmVol, final double forward, final double years) throws SABRException { MultivariateRealFunction estimateRhoAndVol = x -> { double r = x[0]; double v = x[1]; double alpha = findAlpha(forward, forward, atmVol, years, beta, x[0], x[1]); double sumErrors = 0; for (int i = 0; i < volatilities.length; i++) { double modelVol = vol(forward, strikes[i], years, alpha, beta, r, v); sumErrors += Math.pow(modelVol - volatilities[i], 2); } if (Math.abs(r) > 1) { sumErrors = 1e100; } return sumErrors; }; NelderMead nelderMead = new NelderMead(); RealPointValuePair result; try { result = nelderMead.optimize(estimateRhoAndVol, GoalType.MINIMIZE, new double[] { -0.5, 2.6 }); } catch (MathException ex) { throw new SABRException(ex.getMessage(), ex); } double rho = result.getPoint()[0]; double volVol = result.getPoint()[1]; SABRSmileVO params = new SABRSmileVO(); params.setYears(years); params.setRho(rho); params.setVolVol(volVol); params.setAlpha(findAlpha(forward, forward, atmVol, years, beta, rho, volVol)); params.setAtmVol(atmVol); return params; }
From source file:endrov.nucAutoJH.FitGaussian.java
private static double[] fitGaussian2D_(EvPixels p, double sigmaInit, final double midxInit, final double midyInit) { //sigma00, sigma01, sigma11, mu_x, mu_y, c p = p.getReadOnly(EvPixelsType.DOUBLE); final double[] arrPixels = p.getArrayDouble(); final int w = p.getWidth(); final int h = p.getHeight(); int extent = (int) Math.round(3 * sigmaInit); extent = Math.max(extent, 2); final int sx = Math.max(0, (int) (midxInit - extent)); final int ex = Math.min(w, (int) (midxInit + extent + 1)); //+1 to the right? final int sy = Math.max(0, (int) (midyInit - extent)); final int ey = Math.min(h, (int) (midyInit + extent + 1)); double minIntensity = Double.MAX_VALUE; double maxIntensity = Double.MIN_VALUE; for (int y = sy; y < ey; y++) { int base = y * w; double dy2 = y - midyInit; dy2 = dy2 * dy2;/*from w w w .j a va 2s. c om*/ for (int x = sx; x < ex; x++) { double dx2 = x - midxInit; dx2 = dx2 * dx2; double t = arrPixels[base + x]; //if(dx2+dy2<=extent*extent) { if (t < minIntensity) minIntensity = t; if (t > maxIntensity) maxIntensity = t; } } } //double[] weights=new double[]{1}; double[] startPoint = new double[] { sigmaInit, 0, sigmaInit, midxInit, midyInit, minIntensity, maxIntensity - minIntensity }; //double[] output=new double[startPoint.length]; try { MultivariateRealFunction func = new MultivariateRealFunction() { // opt.optimize( public double value(double[] arg) throws FunctionEvaluationException, IllegalArgumentException { double sigma00 = arg[0]; double sigma01 = arg[1]; double sigma11 = arg[2]; double mu0 = arg[3]; double mu1 = arg[4]; double C = arg[5]; double D = arg[6]; double sumError = 0; Matrix2d sigma = new Matrix2d(sigma00, sigma01, sigma01, sigma11); Matrix2d sigmaInv = new Matrix2d(); sigma.invert(sigmaInv); double sigmaDet = sigma.determinant(); double front = 1.0 / (2 * Math.PI * Math.sqrt(sigmaDet)); //System.out.println("front: "+front); //System.out.println("sigma inv "+sigmaInv); if (mu0 < sx || mu0 > ex) sumError += 1000000; if (mu1 < sy || mu1 > ey) sumError += 1000000; if (sigma00 < 1) sumError += 1000000; //if(sigma01<0) sumError+=1000000; if (sigma11 < 1) sumError += 1000000; if (D <= 0) sumError += 1000000; for (int y = sy; y < ey; y++) { int base = y * w; double dy2 = y - midyInit; dy2 = dy2 * dy2; for (int x = sx; x < ex; x++) { double dx2 = x - midxInit; dx2 = dx2 * dx2; double thisReal = arrPixels[base + x]; // if(dx2+dy2<=extent*extent) { // DoubleMatrix2D sigma=new DenseDoubleMatrix2D(new double[][]{{sigma00,sigma01},{sigma01,sigma11}}); //double sigmaDet=sigma00*sigma11-sigma01*sigma01; double dx0 = x - mu0; double dx1 = y - mu1; //http://en.wikipedia.org/wiki/Multivariate_normal_distribution Vector2d vX = new Vector2d(dx0, dx1); Vector2d op = new Vector2d(vX); sigmaInv.transform(op); double upper = -0.5 * op.dot(vX); double exp = Math.exp(upper); //System.out.println("front "+front+" "+exp+" C "+C+" thisreal"+thisReal+" upper "+upper); if (upper > -0.4) exp = 1; else exp = Math.max(0, 1 + upper + 0.4); /* if(exp<0.7) exp=0; else exp=1; */ double thisExpected = D * front * exp + C; double diff = thisExpected - thisReal; sumError += diff * diff; } } } //System.out.println(sigma00+"\t"+sigma01+"\t"+sigma11+"\tC"+C+"\tmu "+mu0+","+mu1+"\terr "+sumError); return sumError; // return new double[]{sumError}; } }; NelderMead opt = new NelderMead(); //LevenbergMarquardtOptimizer opt=new LevenbergMarquardtOptimizer(); opt.setMaxIterations(10000); RealPointValuePair pair = opt.optimize(func, GoalType.MINIMIZE, startPoint); int numit = opt.getIterations(); System.out.println("#it " + numit); System.out.println("err " + func.value(pair.getPointRef())); return pair.getPointRef(); // for(int i=0;i<startPoint.length;i++) // System.out.println("i: "+i+" "+output[i]); //, output, weights, startPoint); } /* catch (MaxIterationsExceededException e) { System.out.println("max it reached"); }*/ catch (Exception e) { e.printStackTrace(); } //Maybe this is a bad point? System.out.println("max it reached"); return startPoint; // return output; }
From source file:endrov.typeLineageAutoNucJH.FitGaussian.java
private static double[] fitGaussian2D_(EvPixels p, double sigmaInit, final double midxInit, final double midyInit) { //sigma00, sigma01, sigma11, mu_x, mu_y, c p = p.getReadOnly(EvPixelsType.DOUBLE); final double[] arrPixels = p.getArrayDouble(); final int w = p.getWidth(); final int h = p.getHeight(); int extent = (int) Math.round(3 * sigmaInit); extent = Math.max(extent, 2); final int sx = Math.max(0, (int) (midxInit - extent)); final int ex = Math.min(w, (int) (midxInit + extent + 1)); //+1 to the right? final int sy = Math.max(0, (int) (midyInit - extent)); final int ey = Math.min(h, (int) (midyInit + extent + 1)); double minIntensity = Double.MAX_VALUE; double maxIntensity = -Double.MAX_VALUE; for (int y = sy; y < ey; y++) { int base = y * w; double dy2 = y - midyInit; dy2 = dy2 * dy2;// w w w .j a v a2 s . c om for (int x = sx; x < ex; x++) { double dx2 = x - midxInit; dx2 = dx2 * dx2; double t = arrPixels[base + x]; //if(dx2+dy2<=extent*extent) { if (t < minIntensity) minIntensity = t; if (t > maxIntensity) maxIntensity = t; } } } //double[] weights=new double[]{1}; double[] startPoint = new double[] { sigmaInit, 0, sigmaInit, midxInit, midyInit, minIntensity, maxIntensity - minIntensity }; //double[] output=new double[startPoint.length]; try { MultivariateRealFunction func = new MultivariateRealFunction() { // opt.optimize( public double value(double[] arg) throws FunctionEvaluationException, IllegalArgumentException { double sigma00 = arg[0]; double sigma01 = arg[1]; double sigma11 = arg[2]; double mu0 = arg[3]; double mu1 = arg[4]; double C = arg[5]; double D = arg[6]; double sumError = 0; Matrix2d sigma = new Matrix2d(sigma00, sigma01, sigma01, sigma11); Matrix2d sigmaInv = new Matrix2d(); sigma.invert(sigmaInv); double sigmaDet = sigma.determinant(); double front = 1.0 / (2 * Math.PI * Math.sqrt(sigmaDet)); //System.out.println("front: "+front); //System.out.println("sigma inv "+sigmaInv); if (mu0 < sx || mu0 > ex) sumError += 1000000; if (mu1 < sy || mu1 > ey) sumError += 1000000; if (sigma00 < 1) sumError += 1000000; //if(sigma01<0) sumError+=1000000; if (sigma11 < 1) sumError += 1000000; if (D <= 0) sumError += 1000000; for (int y = sy; y < ey; y++) { int base = y * w; double dy2 = y - midyInit; dy2 = dy2 * dy2; for (int x = sx; x < ex; x++) { double dx2 = x - midxInit; dx2 = dx2 * dx2; double thisReal = arrPixels[base + x]; // if(dx2+dy2<=extent*extent) { // DoubleMatrix2D sigma=new DenseDoubleMatrix2D(new double[][]{{sigma00,sigma01},{sigma01,sigma11}}); //double sigmaDet=sigma00*sigma11-sigma01*sigma01; double dx0 = x - mu0; double dx1 = y - mu1; //http://en.wikipedia.org/wiki/Multivariate_normal_distribution Vector2d vX = new Vector2d(dx0, dx1); Vector2d op = new Vector2d(vX); sigmaInv.transform(op); double upper = -0.5 * op.dot(vX); double exp = Math.exp(upper); //System.out.println("front "+front+" "+exp+" C "+C+" thisreal"+thisReal+" upper "+upper); if (upper > -0.4) exp = 1; else exp = Math.max(0, 1 + upper + 0.4); /* if(exp<0.7) exp=0; else exp=1; */ double thisExpected = D * front * exp + C; double diff = thisExpected - thisReal; sumError += diff * diff; } } } //System.out.println(sigma00+"\t"+sigma01+"\t"+sigma11+"\tC"+C+"\tmu "+mu0+","+mu1+"\terr "+sumError); return sumError; // return new double[]{sumError}; } }; NelderMead opt = new NelderMead(); //LevenbergMarquardtOptimizer opt=new LevenbergMarquardtOptimizer(); opt.setMaxIterations(10000); RealPointValuePair pair = opt.optimize(func, GoalType.MINIMIZE, startPoint); int numit = opt.getIterations(); System.out.println("#it " + numit); System.out.println("err " + func.value(pair.getPointRef())); return pair.getPointRef(); // for(int i=0;i<startPoint.length;i++) // System.out.println("i: "+i+" "+output[i]); //, output, weights, startPoint); } /* catch (MaxIterationsExceededException e) { System.out.println("max it reached"); }*/ catch (Exception e) { e.printStackTrace(); } //Maybe this is a bad point? System.out.println("max it reached"); return startPoint; // return output; }
From source file:edu.valelab.GaussianFit.ZCalibrator.java
/** * Use the fitfunction to estimate the z position given width in x and y * //from w w w .j a v a 2 s. co m * minimize the distance D in sqrt wx and sqrt wy space * D = sqrt ( square (sqrt wx - sqrt wx, calib) + sqr(sqrt wy - sqrt w, calib) ) * * */ public double getZ(double wx, double wy) { if (!hasFitFunctions()) return 0.0; NelderMead nmx = new NelderMead(); SimpleScalarValueChecker convergedChecker_ = new SimpleScalarValueChecker(1e-6, -1); MultiVariateZFunction mz = new MultiVariateZFunction(fitFunctionWx_, fitFunctionWy_, wx, wy); double[] params0_ = new double[1]; // initial estimates: params0_[0] = 15; // TODO: Need the middle z value of the stack here!!! nmx.setStartConfiguration(params0_); nmx.setConvergenceChecker(convergedChecker_); nmx.setMaxIterations(maxIterations_); double[] paramsOut = { 0.0 }; try { RealPointValuePair result = nmx.optimize(mz, GoalType.MINIMIZE, params0_); paramsOut = result.getPoint(); } catch (java.lang.OutOfMemoryError e) { throw (e); } catch (Exception e) { ij.IJ.log(" " + e.toString()); } return paramsOut[0]; }
From source file:edu.valelab.gaussianfit.fitting.ZCalibrator.java
/** * Use the fitfunction to estimate the z position given width in x and y * /*from ww w. jav a 2 s. c o m*/ * minimize the distance D in sqrt wx and sqrt wy space * D = sqrt ( square (sqrt wx - sqrt wx, calib) + sqr(sqrt wy - sqrt w, calib) ) * * * @param wx - width in x * @param wy - width in y * @return - calculated z position */ public double getZ(double wx, double wy) { if (!hasFitFunctions()) return 0.0; NelderMead nmx = new NelderMead(); SimpleScalarValueChecker convergedChecker_ = new SimpleScalarValueChecker(1e-6, -1); MultiVariateZFunction mz = new MultiVariateZFunction(fitFunctionWx_, fitFunctionWy_, wx, wy); double[] params0_ = new double[1]; // initial estimates: params0_[0] = 15; // TODO: Need the middle z value of the stack here!!! nmx.setStartConfiguration(params0_); nmx.setConvergenceChecker(convergedChecker_); nmx.setMaxIterations(maxIterations_); double[] paramsOut = { 0.0 }; try { RealPointValuePair result = nmx.optimize(mz, GoalType.MINIMIZE, params0_); paramsOut = result.getPoint(); } catch (java.lang.OutOfMemoryError e) { throw (e); } catch (FunctionEvaluationException e) { ij.IJ.log(" " + e.toString()); } catch (OptimizationException e) { ij.IJ.log(" " + e.toString()); } catch (IllegalArgumentException e) { ij.IJ.log(" " + e.toString()); } return paramsOut[0]; }
From source file:edu.valelab.gaussianfit.fitting.ZCalibrator.java
/** * Creates fitFunctionWx_ and fitFunctionWy_ based on data in data_ * //from w w w.j a v a2 s . co m * * @throws org.apache.commons.math.FunctionEvaluationException * @throws org.apache.commons.math.optimization.OptimizationException */ public void fitFunction() throws FunctionEvaluationException, OptimizationException { NelderMead nmx = new NelderMead(); SimpleScalarValueChecker convergedChecker_ = new SimpleScalarValueChecker(1e-6, -1); double[][] wxData = getDataAsArray(0); MultiVariateZCalibrationFunction mvcx = new MultiVariateZCalibrationFunction(wxData); double[] params0_ = new double[5]; // initial estimates: params0_[0] = 37; // TODO: better estimate for c params0_[1] = 200; // Estimate for w0 params0_[2] = 10; // TODO: better estimate for d params0_[3] = 1; // TODO: better estimate for A params0_[4] = 1; // TODO: better estimate for B nmx.setStartConfiguration(params0_); nmx.setConvergenceChecker(convergedChecker_); nmx.setMaxIterations(maxIterations_); double[] paramsOut; RealPointValuePair result = nmx.optimize(mvcx, GoalType.MINIMIZE, params0_); paramsOut = result.getPoint(); // write fit result to Results Table: ResultsTable res = new ResultsTable(); res.incrementCounter(); res.addValue("c", paramsOut[0]); res.addValue("w0", paramsOut[1]); res.addValue("d", paramsOut[2]); res.addValue("A", paramsOut[3]); res.addValue("B", paramsOut[4]); fitFunctionWx_ = paramsOut; double[][] yxData = getDataAsArray(1); MultiVariateZCalibrationFunction yvcx = new MultiVariateZCalibrationFunction(yxData); nmx.setStartConfiguration(params0_); result = nmx.optimize(yvcx, GoalType.MINIMIZE, params0_); paramsOut = result.getPoint(); res.incrementCounter(); res.addValue("c", paramsOut[0]); res.addValue("w0", paramsOut[1]); res.addValue("d", paramsOut[2]); res.addValue("A", paramsOut[3]); res.addValue("B", paramsOut[4]); res.show("Fit Parameters"); fitFunctionWy_ = paramsOut; plotFitFunctions(); }
From source file:edu.valelab.GaussianFit.ZCalibrator.java
/** * Creates fitFunctionWx_ and fitFunctionWy_ based on data in data_ * /*w w w. jav a 2s.c o m*/ * */ public void fitFunction() throws FunctionEvaluationException, OptimizationException { NelderMead nmx = new NelderMead(); SimpleScalarValueChecker convergedChecker_ = new SimpleScalarValueChecker(1e-6, -1); double[][] wxData = getDataAsArray(0); MultiVariateZCalibrationFunction mvcx = new MultiVariateZCalibrationFunction(wxData); double[] params0_ = new double[5]; // initial estimates: params0_[0] = 37; // TODO: better estimate for c params0_[1] = 200; // Estimate for w0 params0_[2] = 10; // TODO: better estimate for d params0_[3] = 1; // TODO: better estimate for A params0_[4] = 1; // TODO: better estimate for B nmx.setStartConfiguration(params0_); nmx.setConvergenceChecker(convergedChecker_); nmx.setMaxIterations(maxIterations_); double[] paramsOut = { 0.0 }; RealPointValuePair result = nmx.optimize(mvcx, GoalType.MINIMIZE, params0_); paramsOut = result.getPoint(); //for (int i = 0; i < paramsOut.length; i++) { // System.out.println("Result " + i + " value: " + (int) paramsOut[i]); //} // write fit result to Results Table: ResultsTable res = new ResultsTable(); res.incrementCounter(); res.addValue("c", paramsOut[0]); res.addValue("w0", paramsOut[1]); res.addValue("d", paramsOut[2]); res.addValue("A", paramsOut[3]); res.addValue("B", paramsOut[4]); fitFunctionWx_ = paramsOut; double[][] yxData = getDataAsArray(1); MultiVariateZCalibrationFunction yvcx = new MultiVariateZCalibrationFunction(yxData); nmx.setStartConfiguration(params0_); result = nmx.optimize(yvcx, GoalType.MINIMIZE, params0_); paramsOut = result.getPoint(); System.out.println("Y:"); //for (int i = 0; i < paramsOut.length; i++) { // System.out.println("Result " + i + " value: " + (int) paramsOut[i]); //} res.incrementCounter(); res.addValue("c", paramsOut[0]); res.addValue("w0", paramsOut[1]); res.addValue("d", paramsOut[2]); res.addValue("A", paramsOut[3]); res.addValue("B", paramsOut[4]); res.show("Fit Parameters"); fitFunctionWy_ = paramsOut; plotFitFunctions(); }
From source file:org.renjin.primitives.optimize.Optimizations.java
/** * General-purpose optimization based on NelderMead, quasi-Newton and conjugate-gradient algorithms. * It includes an option for box-constrained * optimization and simulated annealing. * * @param par initial parameters/*from ww w. j a va 2s .co m*/ * @param fn * @param gradientFunction * @param method * @param controlParameters * @param lower * @param upper * @return */ @Primitive public static ListVector optim(@Current Context context, @Current Environment rho, DoubleVector par, Function fn, SEXP gradientFunction, String method, ListVector controlParameters, DoubleVector lower, DoubleVector upper) { MultivariateRealClosure g = new MultivariateRealClosure(context, rho, fn); if (method.equals("Nelder-Mead")) { NelderMead optimizer = new NelderMead(); try { RealPointValuePair res = optimizer.optimize(g, GoalType.MINIMIZE, par.toDoubleArray()); ListVector.Builder result = new ListVector.Builder(); result.add(new DoubleArrayVector(res.getPoint())); result.add(new DoubleArrayVector(res.getValue())); result.add(new IntArrayVector(IntVector.NA, IntVector.NA)); result.add(new IntArrayVector(0)); result.add(Null.INSTANCE); return result.build(); } catch (FunctionEvaluationException e) { throw new EvalException(e); } catch (OptimizationException e) { throw new EvalException(e); } } else { throw new EvalException("method '" + method + "' not implemented."); } }
From source file:org.renjin.stats.internals.optimize.Optimizations.java
/** * General-purpose optimization based on NelderMead, quasi-Newton and conjugate-gradient algorithms. * It includes an option for box-constrained * optimization and simulated annealing. * * @param par initial parameters/*from ww w . ja v a 2 s.c om*/ * @param fn * @param gradientFunction * @param method * @param controlParameters * @param lower * @param upper * @return */ @Internal public static ListVector optim(@Current Context context, @Current Environment rho, DoubleVector par, Function fn, SEXP gradientFunction, String method, ListVector controlParameters, DoubleVector lower, DoubleVector upper) { MultivariateRealClosure g = new MultivariateRealClosure(context, rho, fn); if (method.equals("Nelder-Mead")) { NelderMead optimizer = new NelderMead(); try { RealPointValuePair res = optimizer.optimize(g, GoalType.MINIMIZE, par.toDoubleArray()); ListVector.Builder result = new ListVector.Builder(); result.add(new DoubleArrayVector(res.getPoint())); result.add(new DoubleArrayVector(res.getValue())); result.add(new IntArrayVector(IntVector.NA, IntVector.NA)); result.add(new IntArrayVector(0)); result.add(Null.INSTANCE); return result.build(); } catch (FunctionEvaluationException e) { throw new EvalException(e); } catch (OptimizationException e) { throw new EvalException(e); } } else { throw new EvalException("method '" + method + "' not implemented."); } }
From source file:r.base.optimize.Optimizations.java
/** * General-purpose optimization based on NelderMead, quasi-Newton and conjugate-gradient algorithms. * It includes an option for box-constrained * optimization and simulated annealing. * * @param par initial parameters// w ww . j a v a2s. c o m * @param fn * @param gradientFunction * @param method * @param controlParameters * @param lower * @param upper * @return */ @Primitive public static ListVector optim(@Current Context context, @Current Environment rho, DoubleVector par, Function fn, SEXP gradientFunction, String method, ListVector controlParameters, DoubleVector lower, DoubleVector upper) { MultivariateRealClosure g = new MultivariateRealClosure(context, rho, fn); if (method.equals("Nelder-Mead")) { NelderMead optimizer = new NelderMead(); try { RealPointValuePair res = optimizer.optimize(g, GoalType.MINIMIZE, par.toDoubleArray()); ListVector.Builder result = new ListVector.Builder(); result.add(new DoubleVector(res.getPoint())); result.add(new DoubleVector(res.getValue())); result.add(new IntVector(IntVector.NA, IntVector.NA)); result.add(new IntVector(0)); result.add(Null.INSTANCE); return result.build(); } catch (FunctionEvaluationException e) { throw new EvalException(e); } catch (OptimizationException e) { throw new EvalException(e); } } else { throw new EvalException("method '" + method + "' not implemented."); } }