List of usage examples for org.apache.commons.math3.optim.univariate UnivariatePointValuePair getPoint
public double getPoint()
From source file:com.google.cloud.genomics.dataflow.functions.verifybamid.Solver.java
/** * Maximizes a univariate function using a grid search followed by Brent's algorithm. * * @param fn the likelihood function to minimize * @param gridStart the lower bound for the grid search * @param gridEnd the upper bound for the grid search * @param gridStep step size for the grid search * @param relErr relative error tolerance for Brent's algorithm * @param absErr absolute error tolerance for Brent's algorithm * @param maxIter maximum # of iterations to perform in Brent's algorithm * @param maxEval maximum # of Likelihood function evaluations in Brent's algorithm * * @return the value of the parameter that maximizes the function *//*from w w w . jav a2 s . co m*/ public static double maximize(UnivariateFunction fn, double gridStart, double gridEnd, double gridStep, double relErr, double absErr, int maxIter, int maxEval) { Interval interval = gridSearch(fn, gridStart, gridEnd, gridStep); BrentOptimizer bo = new BrentOptimizer(relErr, absErr); UnivariatePointValuePair max = bo.optimize(new MaxIter(maxIter), new MaxEval(maxEval), new SearchInterval(interval.getInf(), interval.getSup()), new UnivariateObjectiveFunction(fn), GoalType.MAXIMIZE); return max.getPoint(); }
From source file:com.wwidesigner.optimization.ObjectiveFunctionOptimizer.java
protected static PointValuePair doSingleStart(BaseObjectiveFunction objective, double[] startPoint, int maxEvaluations, double[] nextStart) { singleRunEvaluations = 0;//from w ww .j a va 2 s . c o m PointValuePair result = null; try { int numVariables = objective.getNrDimensions(); if (numVariables > 1) // Use BOBYQA { double trustRegion = objective.getInitialTrustRegionRadius(nextStart); double stoppingTrustRegion = objective.getStoppingTrustRegionRadius(); BOBYQAOptimizer optimizer = new BOBYQAOptimizer(objective.getNrInterpolations(), trustRegion, stoppingTrustRegion); result = runBobyqa(optimizer, objective, nextStart, maxEvaluations); singleRunEvaluations = optimizer.getEvaluations(); } else // Use Brent { BrentOptimizer optimizer = new BrentOptimizer(1.e-6, 1.e-14); UnivariatePointValuePair outcome = runBrent(optimizer, objective, startPoint); result = new PointValuePair(new double[] { outcome.getPoint() }, outcome.getValue()); singleRunEvaluations = optimizer.getEvaluations(); } double value = result.getValue(); if (value == Double.POSITIVE_INFINITY) { System.out.print("no valid solution found"); } else { System.out.print("optimum " + result.getValue()); } } catch (TooManyEvaluationsException e) { System.out.println("Exception: " + e.getMessage()); } // Thrown by BOBYQA for no apparent reason: a bug? catch (NoSuchElementException e) { System.out.println("no valid solution found"); } catch (OperationCancelledException e) { // Restore starting point. objective.setGeometryPoint(startPoint); // Re-throw the exception to give up the whole multi-start // optimization. throw new OperationCancelledException(e.getMessage()); } catch (Exception e) { System.out.println("Exception: " + e.getMessage()); // e.printStackTrace(); } finally { System.out.println(" at start point " + Arrays.toString(nextStart)); } return result; }
From source file:eu.crisis_economics.abm.algorithms.optimization.BrentLineSearch.java
/** * Perform a line search minimization. This function accepts as input: * (a) a starting point (a vector),/*from w w w .jav a 2 s. co m*/ * (b) a direction in which to travel (a vector), * (c) limits on the total distance to travel along (b). * * With these inputs the function attempts to find the minimum of a * scalar-valued multivariate function along the line starting at * (a) and pointing in the direction of (b). * * @param function * A scalar-valued multivariate function to minimize, * @param startingPoint * A vector starting point from which to begin the minimization (P), * @param vectorDirection * A vector direction along which to travel from P, (V) * @param maximumDistanceToTravel * The maximum distance to travel in the direction of V, * @param maximumEvaluations * The maximum number of function evaluations to identify the minimum, * @param relativeErrorGoal * The relative error target of the minimization, * @param absoluteErrorGoal * The absolute error target of the minimization. * @return * A lightweight immutable struct containing the vector solution and * the evaluation of the function at this point. */ static public LineSearchResult doLineSearch(final MultivariateFunction function, final double[] startingPoint, final double[] vectorDirection, final double maximumDistanceToTravel, final int maximumEvaluations, final double relativeErrorGoal, final double absoluteErrorGoal) { Preconditions.checkArgument(maximumEvaluations > 0); Preconditions.checkArgument(relativeErrorGoal > 0. || absoluteErrorGoal > 0.); Preconditions.checkArgument(maximumDistanceToTravel > 0.); final LineSearchObjectiveFunction lineSearcher = new LineSearchObjectiveFunction(function, startingPoint, vectorDirection); final UnivariateOptimizer optimizer = new BrentOptimizer(relativeErrorGoal, absoluteErrorGoal); UnivariatePointValuePair result = optimizer.optimize(new MaxEval(maximumEvaluations), new UnivariateObjectiveFunction(lineSearcher), GoalType.MINIMIZE, new SearchInterval(0, maximumDistanceToTravel, 0)); final double[] vectorSolution = new double[startingPoint.length]; for (int i = 0; i < vectorDirection.length; ++i) vectorSolution[i] = lineSearcher.startingPoint[i] + lineSearcher.normalizedDirection[i] * result.getPoint(); final LineSearchResult solution = new LineSearchResult(vectorSolution, result.getValue()); return solution; }
From source file:edu.uchc.octane.GaussianFitAstigmatism.java
/** * calculate the Z coordinate based on astigmatism *//*from w w w . j a va 2s. c o m*/ void calculateZ() { if (calibration_ == null) { z_ = 0; return; } UnivariateFunction func = new UnivariateFunction() { @Override public double value(double z) { double sigmax = FastMath.sqrt(pvp_.getPoint()[3] / 2); double sigmay = FastMath.sqrt(pvp_.getPoint()[4] / 2); double vx = calibration_[0] + (z - calibration_[1]) * (z - calibration_[1]) * calibration_[2] - sigmax; double vy = calibration_[3] + (z - calibration_[4]) * (z - calibration_[4]) * calibration_[5] - sigmay; return vx * vx + vy * vy; } }; BrentOptimizer optimizer = new BrentOptimizer(1e-4, 1e-4); UnivariatePointValuePair upvp = optimizer.optimize(new UnivariateObjectiveFunction(func), GoalType.MINIMIZE, MaxEval.unlimited(), new SearchInterval(2 * z0min_ - z0max_, 2 * z0max_ - z0min_)); if (upvp.getValue() > errTol_) { throw (new ConvergenceException()); } z_ = upvp.getPoint(); }
From source file:com.itemanalysis.psychometrics.irt.equating.IrtScaleLinking.java
public void computeCoefficients() { raschFamily = checkRaschModel();// ww w. j a v a2 s.c o m ms.setPrecision(precision); mm.setPrecision(precision); hb.setPrecision(precision); sl.setPrecision(precision); if (raschFamily) { double[] sv = { mm.getIntercept() }; UnivariateOptimizer underlying = new BrentOptimizer(1e-10, 1e-14); JDKRandomGenerator g = new JDKRandomGenerator(); //Haebara method MultiStartUnivariateOptimizer optimizer = new MultiStartUnivariateOptimizer(underlying, 5, g);//Five random starts to Brent optimizer. UnivariatePointValuePair hbPair = optimizer.optimize(new MaxEval(500), new UnivariateObjectiveFunction(hb), GoalType.MINIMIZE, new SearchInterval(-4, 4), new InitialGuess(sv)); hb.setIntercept(hbPair.getPoint()); hb.setScale(1.0); fHB = hbPair.getValue(); //Stocking-Lord method UnivariatePointValuePair slPair = optimizer.optimize(new MaxEval(500), new UnivariateObjectiveFunction(sl), GoalType.MINIMIZE, new SearchInterval(-4, 4), new InitialGuess(sv)); sl.setIntercept(slPair.getPoint()); sl.setScale(1.0); fSL = slPair.getValue(); } else { double[] hbStartValues = { mm.getIntercept(), mm.getScale() }; double[] slStartValues = { mm.getIntercept(), mm.getScale() }; if (useUncmin) { DefaultUncminOptimizer optimizer = new DefaultUncminOptimizer(); try { optimizer.minimize(hb, hbStartValues); double[] param = optimizer.getParameters(); fHB = optimizer.getFunctionValue(); hb.setIntercept(param[0]); if (param.length > 1) { hb.setScale(param[1]); } else { hb.setScale(1.0);//Rasch family of models } optimizer.minimize(sl, slStartValues); param = optimizer.getParameters(); fSL = optimizer.getFunctionValue(); sl.setIntercept(param[0]); if (param.length > 1) { sl.setScale(param[1]); } else { sl.setScale(1.0);//Rasch family of models } } catch (UncminException ex) { ex.printStackTrace(); } } else { int numIterpolationPoints = 2 * 2;//two dimensions A and B BOBYQAOptimizer underlying = new BOBYQAOptimizer(numIterpolationPoints); RandomGenerator g = new JDKRandomGenerator(); RandomVectorGenerator generator = new UncorrelatedRandomVectorGenerator(2, new GaussianRandomGenerator(g)); MultiStartMultivariateOptimizer optimizer = new MultiStartMultivariateOptimizer(underlying, 10, generator); PointValuePair hbOptimum = optimizer.optimize(new MaxEval(1000), new ObjectiveFunction(hb), GoalType.MINIMIZE, SimpleBounds.unbounded(2), new InitialGuess(hbStartValues)); double[] hbCoefficients = hbOptimum.getPoint(); hb.setIntercept(hbCoefficients[0]); hb.setScale(hbCoefficients[1]); fHB = hbOptimum.getValue(); PointValuePair slOptimum = optimizer.optimize(new MaxEval(1000), new ObjectiveFunction(sl), GoalType.MINIMIZE, SimpleBounds.unbounded(2), new InitialGuess(slStartValues)); double[] slCoefficients = slOptimum.getPoint(); sl.setIntercept(slCoefficients[0]); sl.setScale(slCoefficients[1]); fSL = slOptimum.getValue(); } } }
From source file:com.itemanalysis.psychometrics.irt.estimation.IrtExaminee.java
/** * Maximum likelihood estimate (MLE) of examinee ability. * * @param thetaMin smallest possible ability estimate (lower bound on BrentOptimizer) * @param thetaMax largest possible ability estimate (upper bound on BrentOptimizer) * @return MLE of examinee ability// w ww .jav a 2s .c o m */ public double maximumLikelihoodEstimate(double thetaMin, double thetaMax) { method = EstimationMethod.ML; UnivariateOptimizer optimizer = new BrentOptimizer(1e-10, 1e-14); UnivariatePointValuePair pair = optimizer.optimize(new MaxEval(100), new UnivariateObjectiveFunction(this), GoalType.MAXIMIZE, new SearchInterval(thetaMin, thetaMax)); estimatedTheta = pair.getPoint(); return estimatedTheta; }
From source file:com.itemanalysis.psychometrics.irt.estimation.IrtExaminee.java
/** * Maximum a Posteriori (MAP) estimate of examinee ability using a normal prior * distribution./*from ww w. j a v a 2s . c o m*/ * * @param mean mean of normal prior distribution * @param sd standard deviation of prior distribution * @param thetaMin smallest possible ability estimate (lower bound on BrentOptimizer) * @param thetaMax largest possible ability estimate (upper bound on BrentOptimizer) * @return MAP estimate of examinee ability */ public double mapEstimate(double mean, double sd, double thetaMin, double thetaMax) { mapPrior = new NormalDistribution(mean, sd); method = EstimationMethod.MAP; UnivariateOptimizer optimizer = new BrentOptimizer(1e-10, 1e-14); UnivariatePointValuePair pair = optimizer.optimize(new MaxEval(100), new UnivariateObjectiveFunction(this), GoalType.MAXIMIZE, new SearchInterval(thetaMin, thetaMax)); estimatedTheta = pair.getPoint(); return estimatedTheta; }
From source file:com.wwidesigner.modelling.PlayingRange.java
/** * Find fmin for a playing range, given fmax. * fmin is the highest frequency <= fmax that satisfies * either gain(fmin) == MinimumGain//from w w w . ja v a 2s .c om * or fmin is a local minimum of Im(Z)/Re(Z). * @param fmax - maximum frequency, as returned by findFmax(). */ public double findFmin(double fmax) { final double stepSize = fmax * Granularity; // Step size for search. // Upper bound on fmin is fmax. // findFmax ensures Im(Z(fmax)) == 0.0. double lowerFreq = fmax; Complex z_lo = calculator.calcZ(fmax, fingering); double g_lo = calculator.calcGain(lowerFreq, z_lo); double ratio = z_lo.getImaginary() / z_lo.getReal(); double minRatio = ratio + 1.0; if (g_lo < MinimumGain) { // Loop gain is too small, even at fmax. // There is no playing range here. throw new NoPlayingRange(fmax); } // Lower bound on fmin either has gain < MinimumGain // or is past a local minimum of Im(Z)/Re(Z). while (g_lo >= MinimumGain && ratio < minRatio) { minRatio = ratio; lowerFreq -= stepSize; if (lowerFreq < fmax / SearchBoundRatio) { throw new NoPlayingRange(fmax); } z_lo = calculator.calcZ(lowerFreq, fingering); g_lo = calculator.calcGain(lowerFreq, z_lo); ratio = z_lo.getImaginary() / z_lo.getReal(); } double freqGain; // Frequency at which gain == MinimumGain. double freqRatio; // Frequency of local minimum of Im(Z)/Re(Z). if (g_lo < MinimumGain) { // Find the point at which gain == MinimumGain. try { freqGain = solver.solve(50, gainOne, lowerFreq, fmax); } catch (Exception e) { System.out.println("Exception solving for fmin (gain): " + e.getMessage()); // e.printStackTrace(); throw new NoPlayingRange(fmax); } } else { freqGain = lowerFreq; } // Find the local minimum of Im(Z)/Re(Z). try { UnivariatePointValuePair minimum; minimum = optimizer.optimize(GoalType.MINIMIZE, new UnivariateObjectiveFunction(zRatio), new MaxEval(50), MaxIter.unlimited(), new SearchInterval(lowerFreq, fmax, 0.5 * (lowerFreq + fmax))); freqRatio = minimum.getPoint(); } catch (Exception e) { System.out.println("Exception solving for fmin (ratio): " + e.getMessage()); // e.printStackTrace(); throw new NoPlayingRange(fmax); } if (freqRatio > freqGain) { return freqRatio; } return freqGain; }
From source file:imagingbook.pub.fd.FourierDescriptor.java
/** * 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). *//*from www . j a v a2s . c o m*/ 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; int K = 25; // number of steps over 180 degrees for (int k = 0; k < K; k++) { final double phi = Math.PI * k / K; // phase to evaluate 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(); return phi0; // the canonical start point phase }
From source file:com.itemanalysis.psychometrics.irt.equating.HaebaraMethodTest.java
/** * This test uses a combination of 3PL and PCM items. Results compared to plink * * plink code:/* w ww .jav a 2 s . com*/ * * library(plink) * * fX<-matrix(c( * 1, -3.188047976, 0,NA,NA, * 1, 1.031760328, 0,NA,NA, * 1, 0.819040914, 0,NA,NA, * 1, -2.706947360, 0,NA,NA, * 1, -0.094527077, 0,NA,NA, * 1, 0.689697135, 0,NA,NA, * 1, -0.551837153, 0,NA,NA, * 1, -0.359559276, 0,NA,NA, * 1, -1.451470831, -0.146619694, -0.636399040, 0.783018734), * nrow=9, byrow=TRUE) * fX<-as.data.frame(fX) * names(fX)<-c("aparam", "bparam","cparam","s1","s2") * * fY<-matrix(c( * 1,-3.074599226,0,NA,NA, * 1,1.01282435,0,NA,NA, * 1,0.868538408,0,NA,NA, * 1,-2.404483603,0,NA,NA, * 1,0.037402866,0,NA,NA, * 1,0.70074742,0,NA,NA, * 1,-0.602555046,0,NA,NA, * 1,-0.350426446,0,NA,NA, * 1,-1.267744832,-0.185885988,-0.61535623,0.801242218), * nrow=9, byrow=TRUE) * fY<-as.data.frame(fY) * names(fY)<-c("aparam", "bparam","cparam","s1","s2") * * common<-cbind(1:9, 1:9) * cat<-c(rep(2,8),4) * * pmX <- as.poly.mod(9,c("drm","gpcm"),list(1:8,9)) * pmY <- as.poly.mod(9,c("drm","gpcm"),list(1:8,9)) * * pars <- as.irt.pars(list(fx=fX,fy=fY), common, cat=list(fx=cat,fy=cat), * poly.mod=list(pmX,pmY), location=c(TRUE,TRUE)) * * plink(pars, startvals=c(1,0), rescale="SL", base.grp=2, D=1.0, symmetric=TRUE) * * */ @Test public void mixedFormtRaschPCMTest() { System.out.println("Mixed format Haebara test: Rasch and PCM"); LinkedHashMap<String, ItemResponseModel> itemFormX = new LinkedHashMap<String, ItemResponseModel>(); LinkedHashMap<String, ItemResponseModel> itemFormY = new LinkedHashMap<String, ItemResponseModel>(); itemFormX.put("Item1", new Irm3PL(-3.188047976, 1.0)); itemFormX.put("Item2", new Irm3PL(1.031760328, 1.0)); itemFormX.put("Item3", new Irm3PL(0.819040914, 1.0)); itemFormX.put("Item4", new Irm3PL(-2.706947360, 1.0)); itemFormX.put("Item5", new Irm3PL(-0.094527077, 1.0)); itemFormX.put("Item6", new Irm3PL(0.689697135, 1.0)); itemFormX.put("Item7", new Irm3PL(-0.551837153, 1.0)); itemFormX.put("Item8", new Irm3PL(-0.359559276, 1.0)); double[] step1x = { -0.146619694, -0.636399040, 0.783018734 }; itemFormX.put("Item9", new IrmPCM(-1.451470831, step1x, 1.0)); itemFormY.put("Item1", new Irm3PL(-3.074599226, 1.0)); itemFormY.put("Item2", new Irm3PL(1.012824350, 1.0)); itemFormY.put("Item3", new Irm3PL(0.868538408, 1.0)); itemFormY.put("Item4", new Irm3PL(-2.404483603, 1.0)); itemFormY.put("Item5", new Irm3PL(0.037402866, 1.0)); itemFormY.put("Item6", new Irm3PL(0.700747420, 1.0)); itemFormY.put("Item7", new Irm3PL(-0.602555046, 1.0)); itemFormY.put("Item8", new Irm3PL(-0.350426446, 1.0)); double[] step1y = { -0.185885988, -0.61535623, 0.801242218 }; itemFormY.put("Item9", new IrmPCM(-1.267744832, step1y, 1.0)); double f = 0; double[] param1 = new double[2]; double[] param2 = new double[2]; UniformDistributionApproximation distX = new UniformDistributionApproximation(-4.0, 4.0, 161);//plink default UniformDistributionApproximation distY = new UniformDistributionApproximation(-4.0, 4.0, 161);//plink default HaebaraMethod hb = new HaebaraMethod(itemFormX, itemFormY, distX, distY, EquatingCriterionType.Q1Q2); hb.setPrecision(6); double[] startValues = { 0.0 };//Using a start value for the intercept only for the Rasch family of models. DefaultUncminOptimizer optimizer = new DefaultUncminOptimizer(); try { optimizer.minimize(hb, startValues); double[] r = optimizer.getParameters(); param1[0] = r[0]; param1[1] = 1.0; f = optimizer.getFunctionValue(); hb.setIntercept(param1[0]); hb.setScale(param1[1]); } catch (UncminException ex) { ex.printStackTrace(); } //Check UNCMIN values against results from plink. System.out.println(" UNCMIN Optimization"); System.out.println(" Iterations: "); System.out.println(" fmin: " + f); System.out.println(" B = " + hb.getIntercept() + " A = " + hb.getScale()); //TODO there's something wrong with plink. I'm getting different results, but sirt agree with mine (at least for the Rasch model). // assertEquals(" Intercept test", 0.105526, hb.getIntercept(), 1e-4);//True results from plink package in R // assertEquals(" Scale test", 1.0, hb.getScale(), 1e-4); //Check Brent optimizer values against plink BrentOptimizer brentOptimizer = new BrentOptimizer(1e-8, 1e-8); UnivariatePointValuePair pair = brentOptimizer.optimize(new MaxEval(200), new UnivariateObjectiveFunction(hb), org.apache.commons.math3.optim.nonlinear.scalar.GoalType.MINIMIZE, new SearchInterval(-3, 3)); param2[0] = pair.getPoint(); param2[1] = 1.0; hb.setIntercept(param2[0]); hb.setScale(param2[1]); f = pair.getValue(); System.out.println(); System.out.println(" Brent Optimization"); System.out.println(" Iterations: "); System.out.println(" fmin: " + f); System.out.println(" B = " + hb.getIntercept() + " A = " + hb.getScale()); //Check UNCMIN values against results from plink. // assertEquals(" Intercept test", 0.105526, hb.getIntercept(), 1e-4); // assertEquals(" Scale test", 1.0, hb.getScale(), 1e-4); // // // //Compare UMNCMIN and Brent optimizer results. // assertEquals(" UNCMIN/Brent intercept test", param1[0], param2[0], 1e-6); // assertEquals(" UNCMIN/Brent slope test", param1[1], param2[1], 1e-6); }