List of usage examples for org.apache.commons.math3.linear RealVector toArray
public double[] toArray()
From source file:eu.crisis_economics.abm.markets.clearing.heterogeneous.BoundedQuadraticEstimationClearingAlgorithm.java
@Override public double applyToNetwork(final MixedClearingNetwork network) { Preconditions.checkNotNull(network); final int dimension = network.getNumberOfEdges(); final ResidualCostFunction aggregateCostFunction = super.getResidualScalarCostFunction(network); final RealVector start = new ArrayRealVector(network.getNumberOfEdges()); start.set(1.0); // Initial rate guess. final BOBYQAOptimizer optimizer = new BOBYQAOptimizer(2 * dimension + 1, 1.2, 1.e-8); final PointValuePair result = optimizer.optimize(new MaxEval(maximumEvaluations), new ObjectiveFunction(aggregateCostFunction), GoalType.MINIMIZE, new SimpleBounds(new double[dimension], ArrayUtil.ones(dimension)), new InitialGuess(start.toArray())); final double residualCost = result.getValue(); System.out.println("Network cleared: residual cost: " + residualCost + "."); return residualCost; }
From source file:com.joptimizer.algebra.CholeskySparseFactorizationTest.java
/** * This test shows that the correct check of the inversion accuracy must be done with * the scaled residual, not with the simple norm ||A.x-b|| *//*from w ww. j a va2s . c o m*/ public void testScaledResidual() throws Exception { log.debug("testScaledResidual"); String matrixId = "1"; double[][] A = Utils .loadDoubleMatrixFromFile("factorization" + File.separator + "matrix" + matrixId + ".csv"); RealMatrix Q = MatrixUtils.createRealMatrix(A); int dim = Q.getRowDimension(); RealVector b = new ArrayRealVector(new double[] { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 }); CholeskySparseFactorization cs = new CholeskySparseFactorization(new SparseDoubleMatrix2D(Q.getData())); cs.factorize(); RealVector x = new ArrayRealVector(cs.solve(F1.make(b.toArray())).toArray()); //scaledResidual = ||Ax-b||_oo/( ||A||_oo . ||x||_oo + ||b||_oo ) // with ||x||_oo = max(x[i]) double residual = Utils.calculateScaledResidual(A, x.toArray(), b.toArray()); log.debug("residual: " + residual); assertTrue(residual < Utils.getDoubleMachineEpsilon()); //b - A.x //checking the simple norm, this will fail double n1 = b.subtract(Q.operate(x)).getNorm(); log.debug("||b - A.x||: " + n1); //assertTrue(n1 < 1.E-8); }
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//from w w w .j av a2s .c o m 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); }
From source file:com.github.thorbenlindhauer.factor.GaussianFactorTest.java
@Test public void testFactorDivision() { GaussianFactor quotient = abcFactor.division(abFactor); // then//from w w w. j av a 2 s . c o m Collection<Variable> newVariables = quotient.getVariables().getVariables(); assertThat(newVariables).hasSize(3); assertThat(newVariables).containsAll(abcFactor.getVariables().getVariables()); assertThat(newVariables).containsAll(abFactor.getVariables().getVariables()); // precision matrix RealMatrix precisionMatrix = quotient.getPrecisionMatrix(); assertThat(precisionMatrix.isSquare()).isTrue(); assertThat(precisionMatrix.getColumnDimension()).isEqualTo(3); double[] row = precisionMatrix.getRowVector(0).toArray(); assertThat(row[0]).isEqualTo(-2.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); assertThat(row[1]).isEqualTo(3.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); assertThat(row[2]).isEqualTo(6.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); row = precisionMatrix.getRowVector(1).toArray(); assertThat(row[0]).isEqualTo(2.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); assertThat(row[1]).isEqualTo(4.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); assertThat(row[2]).isEqualTo(7.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); row = precisionMatrix.getRowVector(2).toArray(); assertThat(row[0]).isEqualTo(10.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); assertThat(row[1]).isEqualTo(3.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); assertThat(row[2]).isEqualTo(5.5d, TestConstants.DOUBLE_VALUE_TOLERANCE); // scaled mean vector RealVector scaledMeanVector = quotient.getScaledMeanVector(); assertThat(scaledMeanVector.getDimension()).isEqualTo(3); double[] meanVectorValues = scaledMeanVector.toArray(); assertThat(meanVectorValues[0]).isEqualTo(0.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); assertThat(meanVectorValues[1]).isEqualTo(0.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); assertThat(meanVectorValues[2]).isEqualTo(1.5d, TestConstants.DOUBLE_VALUE_TOLERANCE); assertThat(quotient.getNormalizationConstant()).isEqualTo(6.3d, TestConstants.DOUBLE_VALUE_TOLERANCE); }
From source file:com.github.thorbenlindhauer.factor.GaussianFactorTest.java
@Test public void testFactorProduct() { // when/* w w w. jav a2 s . c o m*/ GaussianFactor product = acFactor.product(abFactor); // then Collection<Variable> newVariables = product.getVariables().getVariables(); assertThat(newVariables).hasSize(3); assertThat(newVariables).containsAll(acFactor.getVariables().getVariables()); assertThat(newVariables).containsAll(abFactor.getVariables().getVariables()); // precision matrix RealMatrix precisionMatrix = product.getPrecisionMatrix(); assertThat(precisionMatrix.isSquare()).isTrue(); assertThat(precisionMatrix.getColumnDimension()).isEqualTo(3); double[] row = precisionMatrix.getRowVector(0).toArray(); assertThat(row[0]).isEqualTo(6.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); assertThat(row[1]).isEqualTo(1.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); assertThat(row[2]).isEqualTo(2.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); row = precisionMatrix.getRowVector(1).toArray(); assertThat(row[0]).isEqualTo(1.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); assertThat(row[1]).isEqualTo(2.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); assertThat(row[2]).isEqualTo(0.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); row = precisionMatrix.getRowVector(2).toArray(); assertThat(row[0]).isEqualTo(3.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); assertThat(row[1]).isEqualTo(0.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); assertThat(row[2]).isEqualTo(1.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); // scaled mean vector RealVector scaledMeanVector = product.getScaledMeanVector(); assertThat(scaledMeanVector.getDimension()).isEqualTo(3); double[] meanVectorValues = scaledMeanVector.toArray(); assertThat(meanVectorValues[0]).isEqualTo(8.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); assertThat(meanVectorValues[1]).isEqualTo(2.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); assertThat(meanVectorValues[2]).isEqualTo(6.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); assertThat(product.getNormalizationConstant()).isEqualTo(7.7d, TestConstants.DOUBLE_VALUE_TOLERANCE); }
From source file:com.joptimizer.util.MPSParserNetlibTest.java
/** * Tests the parsing of a netlib problem. *//*from ww w . j a v a 2s . 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 w w . j a v a 2 s.c o 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 .jav a 2s .c om*/ 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:com.joptimizer.optimizers.LPPresolverTest.java
private void doPresolving(double[] c, double[][] A, double[] b, double[] lb, double[] ub, double s, double[] expectedSolution, double expectedValue, double expectedTolerance) throws Exception { RealMatrix AMatrix = MatrixUtils.createRealMatrix(A); SingularValueDecomposition dec = new SingularValueDecomposition(AMatrix); int rankA = dec.getRank(); log.debug("p: " + AMatrix.getRowDimension()); log.debug("n: " + AMatrix.getColumnDimension()); log.debug("rank: " + rankA); LPPresolver lpPresolver = new LPPresolver(); lpPresolver.setNOfSlackVariables((short) s); lpPresolver.setExpectedSolution(expectedSolution);//this is just for test! //lpPresolver.setExpectedTolerance(expectedTolerance);//this is just for test! lpPresolver.presolve(c, A, b, lb, ub); int n = lpPresolver.getPresolvedN(); double[] presolvedC = lpPresolver.getPresolvedC().toArray(); double[][] presolvedA = lpPresolver.getPresolvedA().toArray(); double[] presolvedB = lpPresolver.getPresolvedB().toArray(); double[] presolvedLb = lpPresolver.getPresolvedLB().toArray(); double[] presolvedUb = lpPresolver.getPresolvedUB().toArray(); double[] presolvedYlb = lpPresolver.getPresolvedYlb().toArray(); double[] presolvedYub = lpPresolver.getPresolvedYub().toArray(); double[] presolvedZlb = lpPresolver.getPresolvedZlb().toArray(); double[] presolvedZub = lpPresolver.getPresolvedZub().toArray(); log.debug("n : " + n); log.debug("presolvedC : " + ArrayUtils.toString(presolvedC)); log.debug("presolvedA : " + ArrayUtils.toString(presolvedA)); log.debug("presolvedB : " + ArrayUtils.toString(presolvedB)); log.debug("presolvedLb : " + ArrayUtils.toString(presolvedLb)); log.debug("presolvedUb : " + ArrayUtils.toString(presolvedUb)); log.debug("presolvedYlb: " + ArrayUtils.toString(presolvedYlb)); log.debug("presolvedYub: " + ArrayUtils.toString(presolvedYub)); log.debug("presolvedZlb: " + ArrayUtils.toString(presolvedZlb)); log.debug("presolvedZub: " + ArrayUtils.toString(presolvedZub)); //check objective function double delta = expectedTolerance; RealVector presolvedX = MatrixUtils.createRealVector(lpPresolver.presolve(expectedSolution)); log.debug("presolved value: " + MatrixUtils.createRealVector(presolvedC).dotProduct(presolvedX)); RealVector postsolvedX = MatrixUtils.createRealVector(lpPresolver.postsolve(presolvedX.toArray())); double value = MatrixUtils.createRealVector(c).dotProduct(postsolvedX); assertEquals(expectedValue, value, delta); //check postsolved constraints for (int i = 0; i < lb.length; i++) { double di = Double.isNaN(lb[i]) ? -Double.MAX_VALUE : lb[i]; assertTrue(di <= postsolvedX.getEntry(i) + delta); }/*from www . j a v a 2 s . c om*/ for (int i = 0; i < ub.length; i++) { double di = Double.isNaN(ub[i]) ? Double.MAX_VALUE : ub[i]; assertTrue(di + delta >= postsolvedX.getEntry(i)); } RealVector Axmb = AMatrix.operate(postsolvedX).subtract(MatrixUtils.createRealVector(b)); assertEquals(0., Axmb.getNorm(), expectedTolerance); //check presolved constraints assertEquals(presolvedLb.length, presolvedX.getDimension()); assertEquals(presolvedUb.length, presolvedX.getDimension()); AMatrix = MatrixUtils.createRealMatrix(presolvedA); RealVector bvector = MatrixUtils.createRealVector(presolvedB); for (int i = 0; i < presolvedLb.length; i++) { double di = Double.isNaN(presolvedLb[i]) ? -Double.MAX_VALUE : presolvedLb[i]; assertTrue(di <= presolvedX.getEntry(i) + delta); } for (int i = 0; i < presolvedUb.length; i++) { double di = Double.isNaN(presolvedUb[i]) ? Double.MAX_VALUE : presolvedUb[i]; assertTrue(di + delta >= presolvedX.getEntry(i)); } Axmb = AMatrix.operate(presolvedX).subtract(bvector); assertEquals(0., Axmb.getNorm(), expectedTolerance); //check rank(A): must be A pXn with rank(A)=p < n AMatrix = MatrixUtils.createRealMatrix(presolvedA); dec = new SingularValueDecomposition(AMatrix); rankA = dec.getRank(); log.debug("p: " + AMatrix.getRowDimension()); log.debug("n: " + AMatrix.getColumnDimension()); log.debug("rank: " + rankA); assertEquals(AMatrix.getRowDimension(), rankA); assertTrue(rankA < AMatrix.getColumnDimension()); }
From source file:com.joptimizer.solvers.DiagonalKKTSolver.java
/** * Returns two vectors v and w./*from w ww . j av a 2 s. c o 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; }