List of usage examples for org.apache.commons.math3.linear RealVector getDimension
public abstract int getDimension();
From source file:edu.stanford.cfuller.imageanalysistools.fitting.BisquareLinearFit.java
/** * Calculates the weight for the next weighted least squares iteration using the bisquare weighting function. * @param stdAdjR the standardized adjusted residuals, as computed by {@link #calculateStandardizedAdjustedResiduals(RealVector, RealVector, RealVector, RealVector)} * @return a RealVector containing weights for each data point suitable for weighted least squares fitting. *//*from w w w . jav a 2s .c o m*/ protected RealVector calculateBisquareWeights(RealVector stdAdjR) { RealVector bisquareWeights = new ArrayRealVector(stdAdjR.getDimension(), 0.0); for (int i = 0; i < bisquareWeights.getDimension(); i++) { if (Math.abs(stdAdjR.getEntry(i)) < 1) { bisquareWeights.setEntry(i, Math.pow(1 - Math.pow(stdAdjR.getEntry(i), 2), 2)); } } return bisquareWeights; }
From source file:com.joptimizer.util.MPSParserNetlibTest.java
/** * Tests the parsing of a netlib problem. *//*from ww w.ja v a2 s .c o m*/ public void xxxtestSingleNetlib() throws Exception { log.debug("testSingleNetlib"); //String problemName = "afiro"; //String problemName = "afiroPresolved"; //String problemName = "adlittle"; //String problemName = "kb2"; //String problemName = "sc50a"; //String problemName = "sc50b"; //String problemName = "blend"; //String problemName = "scorpion"; //String problemName = "recipe"; //String problemName = "recipePresolved"; //String problemName = "sctap1"; //String problemName = "fit1d"; //String problemName = "israel"; //String problemName = "grow15"; //String problemName = "etamacro"; //String problemName = "pilot"; //String problemName = "pilot4"; //String problemName = "osa-14"; //String problemName = "brandyPresolved"; String problemName = "maros"; File f = Utils.getClasspathResourceAsFile("lp" + File.separator + "netlib" + File.separator + problemName + File.separator + problemName + ".mps"); MPSParser mpsParser = new MPSParser(); mpsParser.parse(f); Properties expectedSolProps = null; try { //this is the solution of the mps problem given by Mathematica expectedSolProps = load(Utils.getClasspathResourceAsFile( "lp" + File.separator + "netlib" + File.separator + problemName + File.separator + "sol.txt")); } catch (Exception e) { } log.debug("name: " + mpsParser.getName()); log.debug("n : " + mpsParser.getN()); log.debug("meq : " + mpsParser.getMeq()); log.debug("mieq: " + mpsParser.getMieq()); log.debug("meq+mieq: " + (mpsParser.getMeq() + mpsParser.getMieq())); List<String> variablesNames = mpsParser.getVariablesNames(); log.debug("x: " + ArrayUtils.toString(variablesNames)); // log.debug("c: " + ArrayUtils.toString(p.getC())); // log.debug("G: " + ArrayUtils.toString(p.getG())); // log.debug("h: " + ArrayUtils.toString(p.getH())); // log.debug("A: " + ArrayUtils.toString(p.getA())); // log.debug("b: " + ArrayUtils.toString(p.getB())); // log.debug("lb:" + ArrayUtils.toString(p.getLb())); // log.debug("ub:" + ArrayUtils.toString(p.getUb())); //check consistency: if the problem was correctly parsed, the expectedSol must be its solution double delta = 1.e-7; if (expectedSolProps != null) { //key = variable name //value = sol value assertEquals(expectedSolProps.size(), variablesNames.size()); RealVector expectedSol = new ArrayRealVector(variablesNames.size()); for (int i = 0; i < variablesNames.size(); i++) { expectedSol.setEntry(i, Double.parseDouble(expectedSolProps.getProperty(variablesNames.get(i)))); } log.debug("expectedSol: " + ArrayUtils.toString(expectedSol.toArray())); //check objective function value Map<String, LPNetlibProblem> problemsMap = LPNetlibProblem.loadAllProblems(); LPNetlibProblem problem = problemsMap.get(problemName); RealVector c = new ArrayRealVector(mpsParser.getC().toArray()); double value = c.dotProduct(expectedSol); log.debug("optimalValue: " + problem.optimalValue); log.debug("value : " + value); assertEquals(problem.optimalValue, value, delta); //check G.x < h if (mpsParser.getG() != null) { RealMatrix G = new Array2DRowRealMatrix(mpsParser.getG().toArray()); RealVector h = new ArrayRealVector(mpsParser.getH().toArray()); RealVector Gxh = G.operate(expectedSol).subtract(h); double maxGxh = -Double.MAX_VALUE; for (int i = 0; i < Gxh.getDimension(); i++) { //log.debug(i); maxGxh = Math.max(maxGxh, Gxh.getEntry(i)); assertTrue(Gxh.getEntry(i) <= 0); } log.debug("max(G.x - h): " + maxGxh); } //check A.x = b if (mpsParser.getA() != null) { RealMatrix A = new Array2DRowRealMatrix(mpsParser.getA().toArray()); RealVector b = new ArrayRealVector(mpsParser.getB().toArray()); RealVector Axb = A.operate(expectedSol).subtract(b); double norm = Axb.getNorm(); log.debug("||A.x -b||: " + norm); assertEquals(0., norm, delta * mpsParser.getN());//some more tolerance } //check upper and lower bounds for (int i = 0; i < mpsParser.getLb().size(); i++) { double di = Double.isNaN(mpsParser.getLb().getQuick(i)) ? -Double.MAX_VALUE : mpsParser.getLb().getQuick(i); assertTrue(di <= expectedSol.getEntry(i)); } for (int i = 0; i < mpsParser.getUb().size(); i++) { double di = Double.isNaN(mpsParser.getUb().getQuick(i)) ? Double.MAX_VALUE : mpsParser.getUb().getQuick(i); assertTrue(di >= expectedSol.getEntry(i)); } } Utils.writeDoubleArrayToFile(mpsParser.getC().toArray(), "target" + File.separator + "c.txt"); Utils.writeDoubleMatrixToFile(mpsParser.getG().toArray(), "target" + File.separator + "G.csv"); Utils.writeDoubleArrayToFile(mpsParser.getH().toArray(), "target" + File.separator + "h.txt"); Utils.writeDoubleMatrixToFile(mpsParser.getA().toArray(), "target" + File.separator + "A.csv"); Utils.writeDoubleArrayToFile(mpsParser.getB().toArray(), "target" + File.separator + "b.txt"); Utils.writeDoubleArrayToFile(mpsParser.getLb().toArray(), "target" + File.separator + "lb.txt"); Utils.writeDoubleArrayToFile(mpsParser.getUb().toArray(), "target" + File.separator + "ub.txt"); }
From source file:edu.stanford.cfuller.imageanalysistools.fitting.GaussianFitter3D.java
/** * Fits a 3D Gaussian to a supplied object, starting from an initial guess of the parameters of that Gaussian. * * The Gaussian is contrained to be symmetric in the x and y dimensions (that is, it will have equal variance in both dimensions). * * @param toFit The {@link ImageObject} to be fit to a Gaussian. * @param initialGuess The initial guess at the parameters of the Gaussian. These must be supplied in the order: amplitude, x-y stddev, z stddev, x position, y position, z position, background. Positions should be supplied in absolute coordinates from the original image, not relative to the box around the object being fit. * @param ppg The number of photons corresponding to one greylevel in the original image. * @return The best fit Gaussian parameters, in the same order as the initial guess had to be supplied. *///from w ww . j av a 2 s. co m public RealVector fit(ImageObject toFit, RealVector initialGuess, double ppg) { //parameter ordering: amplitude, stddev x-y, stddev z, x/y/z coords, background double tol = 1.0e-6; SimplexOptimizer nmm = new SimplexOptimizer(tol, tol); NelderMeadSimplex nms = new NelderMeadSimplex(initialGuess.getDimension()); nmm.setSimplex(nms); PointValuePair pvp = nmm.optimize(10000000, new MLObjectiveFunction(toFit, ppg), org.apache.commons.math3.optimization.GoalType.MINIMIZE, initialGuess.toArray()); RealVector result = new ArrayRealVector(pvp.getPoint()); return result; }
From source file:com.joptimizer.functions.SOCPLogarithmicBarrier.java
public double[][] hessian(double[] X) { RealVector x = new ArrayRealVector(X); RealMatrix ret = new Array2DRowRealMatrix(dim, dim); for (int i = 0; i < socpConstraintParametersList.size(); i++) { SOCPConstraintParameters param = socpConstraintParametersList.get(i); double t = this.buildT(param, x); RealVector u = this.buildU(param, x); double t2uu = t * t - u.dotProduct(u); RealVector t2u = u.mapMultiply(-2 * t); RealMatrix Jacob = this.buildJ(param, x); int k = u.getDimension(); RealMatrix H = new Array2DRowRealMatrix(k + 1, k + 1); RealMatrix ID = MatrixUtils.createRealIdentityMatrix(k); H.setSubMatrix(ID.scalarMultiply(t2uu).add(u.outerProduct(u).scalarMultiply(2)).getData(), 0, 0); H.setSubMatrix(new double[][] { t2u.toArray() }, k, 0); for (int j = 0; j < k; j++) { H.setEntry(j, k, t2u.getEntry(j)); }//from ww w .ja va2 s . c o m H.setEntry(k, k, t * t + u.dotProduct(u)); RealMatrix ret_i = Jacob.multiply(H).multiply(Jacob.transpose()).scalarMultiply(2. / Math.pow(t2uu, 2)); ret = ret.add(ret_i); } return ret.getData(); }
From source file:edu.stanford.cfuller.colocalization3d.fitting.P3DFitter.java
/** * Fits the distances between the two channels of a set of objects to a p3d distribution. * /*ww w .j av a 2 s .c o m*/ * @param objects the ImageObjects whose distances will be fit * @param diffs a RealVector containing the scalar distances between the channels of the ImageObjects, in the same order. * * @return a RealVector containing the parameters for the distribution fit: first the mean parameter, second the standard deviation parameter */ public RealVector fit(List<ImageObject> objects, RealVector diffs) { P3dObjectiveFunction of = new P3dObjectiveFunction(); of.setR(diffs); final double tol = 1e-12; NelderMeadMinimizer nmm = new NelderMeadMinimizer(tol); double initialMean = diffs.getL1Norm() / diffs.getDimension(); double initialWidth = diffs.mapSubtract(initialMean) .map(new org.apache.commons.math3.analysis.function.Power(2)).getL1Norm() / diffs.getDimension(); initialWidth = Math.sqrt(initialWidth); RealVector startingPoint = new ArrayRealVector(2, 0.0); startingPoint.setEntry(0, initialMean); startingPoint.setEntry(1, initialWidth); RealVector parametersMin = null; if (this.parameters.hasKey(ROBUST_P3D_FIT_PARAM)) { double cutoff = this.parameters.getDoubleValueForKey(ROBUST_P3D_FIT_PARAM); of.setMinProb(cutoff); } return nmm.optimize(of, startingPoint); }
From source file:com.joptimizer.solvers.DiagonalKKTSolver.java
/** * Returns two vectors v and w./*from w w w . j ava 2s . co m*/ */ @Override public double[][] solve() throws Exception { RealVector v = null; RealVector w = null; if (Log.isLoggable(MainActivity.JOPTIMIZER_LOGTAG, Log.DEBUG)) { Log.d(MainActivity.JOPTIMIZER_LOGTAG, "H: " + ArrayUtils.toString(H.getData())); Log.d(MainActivity.JOPTIMIZER_LOGTAG, "g: " + ArrayUtils.toString(g.toArray())); } v = new ArrayRealVector(g.getDimension()); for (int i = 0; i < v.getDimension(); i++) { v.setEntry(i, -g.getEntry(i) / H.getEntry(i, i)); } // solution checking if (this.checkKKTSolutionAccuracy && !this.checkKKTSolutionAccuracy(v, w)) { Log.e(MainActivity.JOPTIMIZER_LOGTAG, "KKT solution failed"); throw new Exception("KKT solution failed"); } double[][] ret = new double[2][]; ret[0] = v.toArray(); ret[1] = (w != null) ? w.toArray() : null; return ret; }
From source file:edu.stanford.cfuller.imageanalysistools.fitting.BisquareLinearFit.java
/** * Performs a robust least squares fit with bisquare weights to the supplied data. * * @param indVarValues A RealVector containing the values of the independent variable. * @param depVarValues A RealVector containing the values of the dependent variable. * @return a RealVector containing two elements: the slope of the fit and the y-intercept of the fit. *///from ww w . j av a2 s. c om public RealVector fit(RealVector indVarValues, RealVector depVarValues) { RealVector uniformWeights = new ArrayRealVector(indVarValues.getDimension(), 1.0); RealVector lastParams = new ArrayRealVector(2, Double.MAX_VALUE); RealVector currParams = wlsFit(indVarValues, depVarValues, uniformWeights); RealVector weights = uniformWeights; RealVector leverages = this.calculateLeverages(indVarValues); int c = 0; double norm_mult = 1.0; if (!this.noIntercept) { norm_mult = 2.0; } int maxiter = 10000; while (lastParams.subtract(currParams).getNorm() > CONV_NORM * norm_mult && c++ < maxiter) { lastParams = currParams; RealVector stdAdjR = this.calculateStandardizedAdjustedResiduals(indVarValues, depVarValues, leverages, currParams); weights = calculateBisquareWeights(stdAdjR); currParams = wlsFit(indVarValues, depVarValues, weights); } return currParams; }
From source file:com.github.thorbenlindhauer.factor.GaussianFactorTest.java
@Test public void testFactorMarginalCase2() { Scope variables = newScope(new ContinuousVariable("A")); GaussianFactor aMarginal = abcFactor.marginal(variables); // then/* w w w . j a v a 2 s . c o m*/ Collection<Variable> newVariables = aMarginal.getVariables().getVariables(); assertThat(newVariables).hasSize(1); assertThat(newVariables).contains(new ContinuousVariable("A")); // precision matrix: K_xx - K_xy * K_yy^(-1) * K_yx RealMatrix precisionMatrix = aMarginal.getPrecisionMatrix(); assertThat(precisionMatrix.isSquare()).isTrue(); assertThat(precisionMatrix.getColumnDimension()).isEqualTo(1); double precision = precisionMatrix.getRowVector(0).toArray()[0]; assertThat(precision).isEqualTo(-(14.0d / 3.0d), TestConstants.DOUBLE_VALUE_TOLERANCE); // scaled mean vector: h_x - K_xy * K_yy^(-1) * h_y RealVector scaledMeanVector = aMarginal.getScaledMeanVector(); assertThat(scaledMeanVector.getDimension()).isEqualTo(1); double meanValue = scaledMeanVector.toArray()[0]; assertThat(meanValue).isEqualTo(4.0d / 3.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); // normalization constant: g + 0.5 * (log( det( 2 * PI * K_yy^(-1))) + h_y * K_yy^(-1) * h_y) assertThat(aMarginal.getNormalizationConstant()).isEqualTo(9.324590408d, TestConstants.DOUBLE_VALUE_TOLERANCE); }
From source file:com.bolatu.gezkoncsvlogger.GyroOrientation.RotationKalmanFilter.java
/** * Predict the internal state estimation one time step ahead. * * @param u/*from ww w . j ava2s. c o m*/ * the control vector * @throws DimensionMismatchException * if the dimension of the control vector does not match */ public void predict(final RealVector u) throws DimensionMismatchException { // sanity checks if (u != null && u.getDimension() != controlMatrix.getColumnDimension()) { throw new DimensionMismatchException(u.getDimension(), controlMatrix.getColumnDimension()); } // project the state estimation ahead (a priori state) // xHat(k)- = A * xHat(k-1) + B * u(k-1) // stateEstimation = transitionMatrix.operate(stateEstimation); // add control input if it is available // if (u != null) // { // stateEstimation = stateEstimation.add(controlMatrix.operate(u)); // } // We don't need to use the transition matrix or control matrix, since // we have already done all the work... we can just set the state // estimation to u. if (u != null) { stateEstimation = u; } // project the error covariance ahead // P(k)- = A * P(k-1) * A' + Q errorCovariance = transitionMatrix.multiply(errorCovariance).multiply(transitionMatrixT) .add(processModel.getProcessNoise()); }
From source file:com.github.thorbenlindhauer.factor.GaussianFactorTest.java
@Test public void testFactorMarginalCase1() { Scope variables = newScope(new ContinuousVariable("A"), new ContinuousVariable("C")); GaussianFactor acMarginal = abcFactor.marginal(variables); // then// w w w . j ava 2 s .c om Collection<Variable> newVariables = acMarginal.getVariables().getVariables(); assertThat(newVariables).hasSize(2); assertThat(newVariables).containsAll(variables.getVariables()); // precision matrix: K_xx - K_xy * K_yy^(-1) * K_yx RealMatrix precisionMatrix = acMarginal.getPrecisionMatrix(); assertThat(precisionMatrix.isSquare()).isTrue(); assertThat(precisionMatrix.getColumnDimension()).isEqualTo(2); double[] row = precisionMatrix.getRowVector(0).toArray(); assertThat(row[0]).isEqualTo(1, TestConstants.DOUBLE_VALUE_TOLERANCE); assertThat(row[1]).isEqualTo(4.0d / 3.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); row = precisionMatrix.getRowVector(1).toArray(); assertThat(row[0]).isEqualTo(8.5d, TestConstants.DOUBLE_VALUE_TOLERANCE); assertThat(row[1]).isEqualTo(2.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); // scaled mean vector: h_x - K_xy * K_yy^(-1) * h_y RealVector scaledMeanVector = acMarginal.getScaledMeanVector(); assertThat(scaledMeanVector.getDimension()).isEqualTo(2); double[] meanValues = scaledMeanVector.toArray(); assertThat(meanValues[0]).isEqualTo(5.0d / 3.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); assertThat(meanValues[1]).isEqualTo(0.5d, TestConstants.DOUBLE_VALUE_TOLERANCE); // normalization constant: g + 0.5 * (log( det( 2 * PI * K_yy^(-1))) + h_y * K_yy^(-1) * h_y) assertThat(acMarginal.getNormalizationConstant()).isEqualTo(8.856392131, TestConstants.DOUBLE_VALUE_TOLERANCE); }