List of usage examples for org.apache.commons.math3.analysis.solvers UnivariateSolver solve
double solve(int maxEval, FUNC f, double min, double max);
From source file:jat.examples.CommonsMath.RootFinder.java
/** * @param args/*from ww w . j a va2s.c o m*/ */ public static void main(String[] args) { final double relativeAccuracy = 1.0e-12; final double absoluteAccuracy = 1.0e-8; UnivariateFunction function = new MyFunction(); UnivariateSolver nonBracketing = new BrentSolver(relativeAccuracy, absoluteAccuracy); // double baseRoot = bs.solve(100, function, -2.0, 0); System.out.println("Roots of f(x)=x^2-2: "); double baseRoot; baseRoot = nonBracketing.solve(100, function, -2.0, 0); System.out.println("root1: " + baseRoot); baseRoot = nonBracketing.solve(100, function, 0, 2.0); System.out.println("root2: " + baseRoot); // System.out.println("rs: " + Math.sqrt(2.)); }
From source file:gamlss.smoothing.PB.java
/** * The main fitting method, initiate appropriate smoothing * function according to incoming parameters. * @param additiveFit -object of AdditiveFoit class * @return reiduals// ww w .java2 s.c o m */ //gamlss.pb <- function(x, y, w, xeval = NULL, ...) public ArrayRealVector solve(final AdditiveFit additiveFit) { Double lambda = Controls.LAMBDAS_USER_DEFINED; colNum = additiveFit.getColNum(); y = additiveFit.getZ(); w = additiveFit.getW(); whichDistParameter = additiveFit.getWhichDistParameter(); basisM = (BlockRealMatrix) getBasisM().get(colNum); penaltyM = (BlockRealMatrix) getPenaltyM().get(colNum); df = (Integer) getDfValues().get(colNum); //n <- nrow(X) # the no of observations n = basisM.getRowDimension(); //lambdaS <- get(startLambdaName, envir=gamlss.env) //geting the starting value lambdaS = getLambdas().get(whichDistParameter); //if (lambdaS>=1e+07) lambda <- 1e+07 if (lambdaS >= 1e+07) { lambda = 1e+07; } //if (lambdaS<=1e-07) lambda <- 1e-07 if (lambdaS <= 1e-07) { lambda = 1e-07; } //if (is.null(df)&&!is.null(lambda)||!is.null(df)&&!is.null(lambda)) if (lambda != null) { //fit <- regpen(y, X, w, lambda, D) beta = regpen(lambda); //fv <- X %*% fit$beta fv = (ArrayRealVector) basisM.operate(beta); //else if (is.null(df)&&is.null(lambda)) } else if (df == null) { //lambda <- lambdaS lambda = lambdaS; //switch(control$c," ML"={ switch (Controls.SMOOTH_METHOD) { case DistributionSettings.GAIC: fv = functionGAIC(lambda); break; case DistributionSettings.ML: fv = functionML(lambda); break; case DistributionSettings.ML1: fv = functionML1(lambda); break; case DistributionSettings.EM: System.err.println("EM has not been implemented yet"); break; case DistributionSettings.GCV: fv = functionGCV(lambda); break; default: System.err.println( " Cannot recognise the " + "smoothing method or it has" + "not been implemented yet"); } } else { //QR <- qr(sqrt(w)*X) //Rinv <- solve(qr.R(QR)) rM = (BlockRealMatrix) new QRDecomposition( MatrixFunctions.multVectorMatrix(MatrixFunctions.sqrtVec(w), basisM)).getR(); rM = rM.getSubMatrix(0, rM.getColumnDimension() - 1, 0, rM.getColumnDimension() - 1); rMinv = (BlockRealMatrix) new QRDecomposition(rM).getSolver().getInverse(); //S <- t(D)%*%D sM = penaltyM.transpose().multiply(penaltyM); //UDU <- eigen(t(Rinv)%*%S%*%Rinv) uduV = new EigenDecomposition(rMinv.transpose().multiply(sM).multiply(rMinv)).getRealEigenvalues(); //lambda <- if (sign(edf1_df(0))==sign(edf1_df(100000))) 100000 //in case they have the some sign edfTemp1 = edf1_df(0, df); edfTemp2 = edf1_df(100000.0, df); if (FastMath.signum(edfTemp1) == FastMath.signum(edfTemp2)) { lambda = 100000.0; } else { //else uniroot(edf1_df, c(0,100000))$root uniRootObj.setDf(df); final double relativeAccuracy = 1.0e-12; final double absoluteAccuracy = 1.0e-8; UnivariateSolver uniRootSolver = new BrentSolver(relativeAccuracy, absoluteAccuracy); lambda = uniRootSolver.solve(1000, uniRootObj, 0.0, 100000.0); } //fit <- regpen(y, X, w, lambda, D) beta = regpen(lambda); fv = (ArrayRealVector) basisM.operate(beta); } if (Controls.IF_NEED_THIS) { //waug <- as.vector(c(w, rep(1,nrow(D)))) waug = w.append(MatrixFunctions.repV(1.0, penaltyM.getRowDimension())); //xaug <- as.matrix(rbind(X,sqrt(lambda)*D)) xaug = MatrixFunctions.appendMatricesRows(basisM, (BlockRealMatrix) penaltyM.scalarMultiply(FastMath.sqrt(lambda))); //lev <- hat(sqrt(waug)*xaug,intercept=FALSE)[1:n] //get the hat matrix lev = (ArrayRealVector) MatrixFunctions.getMainDiagonal(new BlockRealMatrix(MatrixFunctions .calculateHat(MatrixFunctions.multVectorMatrix(MatrixFunctions.sqrtVec(waug), xaug)).getData())) .getSubVector(0, n); //lev <- (lev-.hat.WX(w,x)) //subtract the linear since is already fitted lev = lev.subtract( hatWX((ArrayRealVector) getSmootherMatrices().get(whichDistParameter).getColumnVector(colNum))); // var <- lev/w //the variance of the smootherz var = lev.ebeDivide(w); } //if (is.null(xeval)) # if no prediction if (Controls.XEVAL_USER_DEFINED == null) { // Residuals return y.subtract(fv); //else # for prediction } else { //ll <- dim(as.matrix(attr(x,"X")))[1] //nx <- as.matrix(attr(x,"X"))[seq(length(y)+1,ll),] tempArr = ArithmeticSeries.getSeries(y.getDimension() + 1, basisM.getRowDimension(), 1); BlockRealMatrix nx = new BlockRealMatrix(tempArr.length, basisM.getColumnDimension()); for (int i = 0; i < tempArr.length; i++) { nx.setRowVector(i, basisM.getRowVector((int) tempArr[i])); } //pred <- drop(nx %*% fit$beta) ArrayRealVector pred = (ArrayRealVector) nx.operate(beta); // Residuals return y.subtract(fv); } }
From source file:org.micromanager.asidispim.fit.Fitter.java
/** * Finds the x value corresponding to the maximum function value within the * range of the provided data set.// w w w . j a v a 2s . co m * * @param type one of the Fitter.FunctionType predefined functions * @param parms parameters describing the function. These need to match the * selected function or an IllegalArgumentEception will be thrown * @param data JFreeChart series, used to bracket the range in which the * maximum will be found * * @return x value corresponding to the maximum function value */ public static double getXofMaxY(XYSeries data, FunctionType type, double[] parms) { double xAtMax = 0.0; double minX = data.getMinX(); double maxX = data.getMaxX(); switch (type) { case NoFit: // find the position in data with the highest y value double highestScore = data.getY(0).doubleValue(); int highestIndex = 0; for (int i = 1; i < data.getItemCount(); i++) { double newVal = data.getY(i).doubleValue(); if (newVal > highestScore) { highestScore = newVal; highestIndex = i; } } return data.getX(highestIndex).doubleValue(); case Pol1: case Pol2: case Pol3: checkParms(type, parms); PolynomialFunction derivativePolFunction = (new PolynomialFunction(parms)).polynomialDerivative(); final double relativeAccuracy = 1.0e-12; final double absoluteAccuracy = 1.0e-8; final int maxOrder = 5; UnivariateSolver solver = new BracketingNthOrderBrentSolver(relativeAccuracy, absoluteAccuracy, maxOrder); xAtMax = solver.solve(100, derivativePolFunction, minX, maxX); break; case Gaussian: // for a Gaussian we can take the mean and be sure it is the maximum // note that this may be outside our range of X values, but // this will be caught by our sanity checks below xAtMax = parms[1]; } // sanity checks if (xAtMax > maxX) xAtMax = maxX; if (xAtMax < minX) xAtMax = minX; return xAtMax; }
From source file:org.orekit.frames.TopocentricFrame.java
/** * Compute the limit visibility point for a satellite in a given direction. * <p>//from w ww . ja va 2 s. co m * This method can be used to compute visibility circles around ground stations * for example, using a simple loop on azimuth, with either a fixed elevation * or an elevation that depends on azimuth to take ground masks into account. * </p> * @param radius satellite distance to Earth center * @param azimuth pointing azimuth from station * @param elevation pointing elevation from station * @return limit visibility point for the satellite * @throws OrekitException if point cannot be found */ public GeodeticPoint computeLimitVisibilityPoint(final double radius, final double azimuth, final double elevation) throws OrekitException { try { // convergence threshold on point position: 1mm final double deltaP = 0.001; final UnivariateSolver solver = new BracketingNthOrderBrentSolver( deltaP / Constants.WGS84_EARTH_EQUATORIAL_RADIUS, deltaP, deltaP, 5); // find the distance such that a point in the specified direction and at the solved-for // distance is exactly at the specified radius final double distance = solver.solve(1000, new UnivariateFunction() { /** {@inheritDoc} */ public double value(final double x) { try { final GeodeticPoint gp = pointAtDistance(azimuth, elevation, x); return parentShape.transform(gp).getNorm() - radius; } catch (OrekitException oe) { throw new OrekitExceptionWrapper(oe); } } }, 0, 2 * radius); // return the limit point return pointAtDistance(azimuth, elevation, distance); } catch (TooManyEvaluationsException tmee) { throw new OrekitException(tmee); } catch (OrekitExceptionWrapper lwe) { throw lwe.getException(); } }
From source file:org.orekit.models.earth.Geoid.java
/** * {@inheritDoc}/* w ww.ja v a 2s . c o m*/ * * <p> The intersection point is computed using a line search along the * specified line. This is accurate when the geoid is slowly varying. */ @Override public GeodeticPoint getIntersectionPoint(final Line lineInFrame, final Vector3D closeInFrame, final Frame frame, final AbsoluteDate date) throws OrekitException { /* * It is assumed that the geoid is slowly varying over it's entire * surface. Therefore there will one local intersection. */ // transform to body frame final Frame bodyFrame = this.getBodyFrame(); final Transform frameToBody = frame.getTransformTo(bodyFrame, date); final Vector3D close = frameToBody.transformPosition(closeInFrame); final Line lineInBodyFrame = frameToBody.transformLine(lineInFrame); // set the line's direction so the solved for value is always positive final Line line; if (lineInBodyFrame.getAbscissa(close) < 0) { line = lineInBodyFrame.revert(); } else { line = lineInBodyFrame; } final ReferenceEllipsoid ellipsoid = this.getEllipsoid(); // calculate end points // distance from line to center of earth, squared final double d2 = line.pointAt(0.0).getNormSq(); // the minimum abscissa, squared final double minAbscissa2 = FastMath.pow(ellipsoid.getPolarRadius() + MIN_UNDULATION, 2) - d2; // smaller end point of the interval = 0.0 or intersection with // min_undulation sphere final double lowPoint = FastMath.sqrt(FastMath.max(minAbscissa2, 0.0)); // the maximum abscissa, squared final double maxAbscissa2 = FastMath.pow(ellipsoid.getEquatorialRadius() + MAX_UNDULATION, 2) - d2; // larger end point of the interval final double highPoint = FastMath.sqrt(maxAbscissa2); // line search function final UnivariateFunction heightFunction = new UnivariateFunction() { @Override public double value(final double x) { try { final GeodeticPoint geodetic = transform(line.pointAt(x), bodyFrame, date); return geodetic.getAltitude(); } catch (OrekitException e) { // due to frame transform -> re-throw throw new RuntimeException(e); } } }; // compute answer if (maxAbscissa2 < 0) { // ray does not pierce bounding sphere -> no possible intersection return null; } // solve line search problem to find the intersection final UnivariateSolver solver = new BracketingNthOrderBrentSolver(); try { final double abscissa = solver.solve(MAX_EVALUATIONS, heightFunction, lowPoint, highPoint); // return intersection point return this.transform(line.pointAt(abscissa), bodyFrame, date); } catch (TooManyEvaluationsException e) { // no intersection return null; } catch (MathIllegalArgumentException e) { // no intersection return null; } }