List of usage examples for org.apache.commons.math3.linear Array2DRowRealMatrix operate
@Override public double[] operate(final double[] v) throws DimensionMismatchException
From source file:fp.overlapr.algorithmen.StressMajorization.java
@Deprecated private static ArrayRealVector conjugateGradientsMethod(Array2DRowRealMatrix A, ArrayRealVector b, ArrayRealVector werte) {/*from w w w.ja va 2s .co m*/ Array2DRowRealMatrix preJacobi = new Array2DRowRealMatrix(A.getRowDimension(), A.getColumnDimension()); // Predconditioner berechnen preJacobi.walkInRowOrder(new DefaultRealMatrixChangingVisitor() { @Override public double visit(int row, int column, double value) { if (row == column) { return 1 / A.getEntry(row, column); } else { return 0.0; } } }); // x_k beliebig whlen ArrayRealVector x_k = new ArrayRealVector(werte); // r_k berechnen ArrayRealVector r_k = b.subtract(A.operate(x_k)); // h_k berechnen ArrayRealVector h_k = (ArrayRealVector) preJacobi.operate(r_k); // d_k = r_k ArrayRealVector d_k = h_k; // x_k+1 und r_k+1 und d_k+1, sowie alpha und beta und z erzeugen ArrayRealVector x_k1; ArrayRealVector r_k1; ArrayRealVector d_k1; ArrayRealVector h_k1; double alpha; double beta; ArrayRealVector z; do { // Speichere Matrix-Vektor-Produkt, um es nur einmal auszurechnen z = (ArrayRealVector) A.operate(d_k); // Finde von x_k in Richtung d_k den Ort x_k1 des Minimums der // Funktion E // und aktualisere den Gradienten bzw. das Residuum alpha = r_k.dotProduct(h_k) / d_k.dotProduct(z); x_k1 = x_k.add(d_k.mapMultiply(alpha)); r_k1 = r_k.subtract(z.mapMultiply(alpha)); h_k1 = (ArrayRealVector) preJacobi.operate(r_k1); // Korrigiere die Suchrichtung d_k1 beta = r_k1.dotProduct(h_k1) / r_k.dotProduct(h_k); d_k1 = h_k1.add(d_k.mapMultiply(beta)); // Werte "eins" weitersetzen x_k = x_k1; r_k = r_k1; d_k = d_k1; h_k = h_k1; } while (r_k1.getNorm() > TOL); return x_k1; }
From source file:fi.smaa.jsmaa.model.MultivariateGaussianCriterionMeasurement.java
@Override public void deleteAlternative(Alternative alt) { final int n = alternatives.size(); final int index = alternatives.indexOf(alt); alternatives.remove(alt);//from w ww . ja v a 2 s .c om double select[][] = new double[n - 1][n]; for (int i = 0; i < n; ++i) { if (i < index) { select[i][i] = 1; } else if (i > index) { select[i - 1][i] = 1; } } Array2DRowRealMatrix m = new Array2DRowRealMatrix(select); final RealVector newMeanVector = m.operate(meanVector); final RealMatrix newCovarianceMatrix = covarianceMatrix.multiply(m.transpose()).preMultiply(m); fireEvents(setMeanVectorInternal(newMeanVector), setCovarianceMatrixInternal(newCovarianceMatrix)); }
From source file:fi.smaa.jsmaa.model.MultivariateGaussianCriterionMeasurement.java
@Override public void reorderAlternatives(List<Alternative> alts) { final int n = alternatives.size(); double permute[][] = new double[n][n]; for (int i = 0; i < n; ++i) { int j = alternatives.indexOf(alts.get(i)); permute[i][j] = 1;// w ww . jav a 2s . c o m } Array2DRowRealMatrix m = new Array2DRowRealMatrix(permute); alternatives.clear(); alternatives.addAll(alts); final RealVector newMeanVector = m.operate(meanVector); final RealMatrix newCovarianceMatrix = covarianceMatrix.multiply(m.transpose()).preMultiply(m); fireEvents(setMeanVectorInternal(newMeanVector), setCovarianceMatrixInternal(newCovarianceMatrix)); }
From source file:fp.overlapr.algorithmen.StressMajorization.java
/** * Fhrt die Stress-Majorization durch. siehe: Gansner, Koren, North: Graph * Drawing by Stress Majorization, 2004/*from w w w . ja v a 2 s . c o m*/ * * @param graph * Graph, dessen Knoten-Positionen neu berechnet werden sollen * @param d * Matrix, die die idealen Distanzen (d_ij) zwischen den Knoten * enthlt * @return Matrix, die die neuen x- und y-Werte der einzelnen Knoten enthlt */ public static double[][] doStressMajorization(Graph graph, double[][] d) { int iter = 0; // X holen Array2DRowRealMatrix X = new Array2DRowRealMatrix(graph.getKnotenAnzahl(), 2); for (int i = 0; i < graph.getKnotenAnzahl(); i++) { X.setEntry(i, 0, graph.getKnoten().get(i).getX()); X.setEntry(i, 1, graph.getKnoten().get(i).getY()); } // D holen Array2DRowRealMatrix D = new Array2DRowRealMatrix(d); // W berechnen Array2DRowRealMatrix W = new Array2DRowRealMatrix(D.getRowDimension(), D.getColumnDimension()); W.walkInRowOrder(new DefaultRealMatrixChangingVisitor() { @Override public double visit(int row, int column, double value) { if (D.getEntry(row, column) == 0) { return 0.0; } else { return 1.0 / (D.getEntry(row, column) * D.getEntry(row, column)); } } }); // LW berechnen Array2DRowRealMatrix LW = new Array2DRowRealMatrix(D.getRowDimension(), D.getColumnDimension()); LW.walkInRowOrder(new DefaultRealMatrixChangingVisitor() { @Override public double visit(int row, int column, double value) { if (row != column) { return (-1) * W.getEntry(row, column); } else { return value; } } }); LW.walkInRowOrder(new DefaultRealMatrixChangingVisitor() { @Override public double visit(int row, int column, double value) { if (row == column) { double sum = 0; for (int k = 0; k < LW.getColumnDimension(); k++) { if (k != row) { sum = sum + W.getEntry(row, k); } } return sum; } else { return value; } } }); double[][] x = null; while (iter < ITER) { iter++; // LX berechnen Array2DRowRealMatrix LX = new Array2DRowRealMatrix(D.getRowDimension(), D.getColumnDimension()); LX.walkInRowOrder(new DefaultRealMatrixChangingVisitor() { @Override public double visit(int row, int column, double value) { if (row != column) { // norm 2 double term1 = FastMath.pow((X.getEntry(row, 0) - X.getEntry(column, 0)), 2); double term2 = FastMath.pow((X.getEntry(row, 1) - X.getEntry(column, 1)), 2); double abst = Math.sqrt(term1 + term2); return (-1) * W.getEntry(row, column) * D.getEntry(row, column) * inv(abst); } else { return value; } } }); LX.walkInRowOrder(new DefaultRealMatrixChangingVisitor() { @Override public double visit(int row, int column, double value) { if (row == column) { double sum = 0; for (int k = 0; k < LX.getColumnDimension(); k++) { if (k != row) { sum = sum + LX.getEntry(row, k); } } return (-1) * sum; } else { return value; } } }); /* * Lineare Gleichungssysteme lsen */ // x-Werte holen ArrayRealVector xWerte = new ArrayRealVector(X.getColumn(0)); // y-Werte holen ArrayRealVector yWerte = new ArrayRealVector(X.getColumn(1)); // b_x berechnen ArrayRealVector b_x = (ArrayRealVector) LX.operate(xWerte); // b_y berechnen ArrayRealVector b_y = (ArrayRealVector) LX.operate(yWerte); /* * CG-Verfahren anwenden */ // neue x-Werte berechnen mittels PCG // xWerte = conjugateGradientsMethod(LW, b_x, xWerte); // neue y-Werte berechnen mittels PCG // yWerte = conjugateGradientsMethod(LW, b_y, yWerte); ConjugateGradient cg = new ConjugateGradient(Integer.MAX_VALUE, TOL, false); xWerte = (ArrayRealVector) cg.solve(LW, JacobiPreconditioner.create(LW), b_x); yWerte = (ArrayRealVector) cg.solve(LW, JacobiPreconditioner.create(LW), b_y); /* * neue Positiones-Werte zurckgeben */ x = new double[X.getRowDimension()][2]; for (int i = 0; i < x.length; i++) { x[i][0] = xWerte.getEntry(i); x[i][1] = yWerte.getEntry(i); X.setEntry(i, 0, xWerte.getEntry(i)); X.setEntry(i, 1, yWerte.getEntry(i)); } } return x; }
From source file:org.apache.hadoop.hive.ql.abm.simulation.MCSimNode.java
public List<SimulationResult> simulate(IntArrayList[] requests) { init(requests);//from w ww. java2s. co m if (dimension == 0) { return null; } if (parent == null) { return defaultSimulate(); } List<SimulationResult> parRet = parent.simulate(targets); if (parRet == null) { return defaultSimulate(); } List<SimulationResult> ret = new ArrayList<SimulationResult>(); for (SimulationResult res : parRet) { ArrayList<IntArrayList> condIds = res.lastCondId; boolean[] fake = new boolean[dimension]; double[] mu = new double[dimension]; double[][] A = new double[dimension][dimension]; double[][] B = new double[res.invSigma.getColumnDimension()][dimension]; double[] zero = new double[dimension]; for (int i = 0; i < within1.length; ++i) { IntArrayList cIds = condIds.get(i); within1[i].fill(cIds, fake, mu, A); InterDistOracle[] w2s = within2[i]; for (int j = i + 1; j < within2.length; ++j) { w2s[j].fillSym(cIds, condIds.get(j), fake, mu, A); } } int dif = between.size() - res.condIds.size(); double[] pmu = res.getMean(); for (int i = 0; i < targets.length; ++i) { int cum = 0; IntArrayList cIds = condIds.get(i); for (int k = 0; k < res.condIds.size(); ++k) { ArrayList<IntArrayList> condIds2 = res.condIds.get(k); InterDistOracle[] os = between.get(k + dif)[i]; for (int j = 0; j < os.length; ++j) { cum += os[j].fillAsym(cIds, condIds2.get(j), fake, mu, pmu, B, cum); } } } for (int k = 0, cum = 0; k < res.condIds.size(); ++k) { ArrayList<IntArrayList> pCIds = res.condIds.get(k); InterDistOracle[][] oss = between.get(k + dif); for (int j = 0; j < oss.length; ++j) { IntArrayList cIds = pCIds.get(j); InterDistOracle[] os = oss[j]; int off = 0; for (int i = 0; i < os.length; ++i) { off = os[i].fillAsym(cIds, condIds.get(i), fake, pmu, mu, B, cum); } cum += off; } } Array2DRowRealMatrix a = new Array2DRowRealMatrix(A); Array2DRowRealMatrix b = new Array2DRowRealMatrix(B); Array2DRowRealMatrix c = (Array2DRowRealMatrix) b.transpose(); Array2DRowRealMatrix tmp = c.multiply(res.invSigma); Array2DRowRealMatrix sigma = a.subtract(tmp.multiply(b)); double[] scale = correct(sigma.getDataRef()); MultivariateNormalDistribution dist = new MultivariateNormalDistribution(zero, sigma.getDataRef()); double[][] smpls = dist.sample(res.samples.size()); int pos = 0; for (double[] smpl : smpls) { if (scale != null) { restore(smpl, scale); } fixFake(fake, mu, smpl); double[] s = res.getSample(pos); subtract(s, pmu); add(smpl, tmp.operate(s)); res.samples.get(pos)[level] = smpl; ++pos; } res.means.add(mu); res.invSigma = new Array2DRowRealMatrix( new LUDecomposition(concat(a, c, b, res.invSigma)).getSolver().getInverse().getData()); dispatch(res, ret); } return ret; }
From source file:org.eclipse.dawnsci.analysis.dataset.roi.fitting.EllipseCoordinatesFunction.java
/** * least-squares using algebraic cost function * <p>/*from ww w .ja v a 2 s . c o m*/ * This uses the method in "Numerically stable direct least squares fitting of ellipses" * by R. Halir and J. Flusser, Proceedings of the 6th International Conference in Central Europe * on Computer Graphics and Visualization. WSCG '98. CZ, Pilsen 1998, pp 125-132. * <p> * @param ix * @param iy * @return geometric parameters */ private static double[] quickfit(IDataset ix, IDataset iy) { Dataset x = DatasetUtils.convertToDataset(ix); Dataset y = DatasetUtils.convertToDataset(iy); final Dataset xx = Maths.square(x); final Dataset yy = Maths.square(y); final Dataset xxx = Maths.multiply(xx, x); final Dataset yyy = Maths.multiply(yy, y); final Dataset xy = Maths.multiply(x, y); Matrix S1 = new Matrix(3, 3); S1.set(0, 0, LinearAlgebra.dotProduct(xx, xx).getDouble()); S1.set(0, 1, LinearAlgebra.dotProduct(xxx, y).getDouble()); S1.set(0, 2, LinearAlgebra.dotProduct(xx, yy).getDouble()); S1.set(1, 0, S1.get(0, 1)); S1.set(1, 1, S1.get(0, 2)); S1.set(1, 2, LinearAlgebra.dotProduct(x, yyy).getDouble()); S1.set(2, 0, S1.get(0, 2)); S1.set(2, 1, S1.get(1, 2)); S1.set(2, 2, LinearAlgebra.dotProduct(yy, yy).getDouble()); Matrix S2 = new Matrix(3, 3); S2.set(0, 0, ((Number) xxx.sum()).doubleValue()); S2.set(0, 1, LinearAlgebra.dotProduct(xx, y).getDouble()); S2.set(0, 2, ((Number) xx.sum()).doubleValue()); S2.set(1, 0, S2.get(0, 1)); S2.set(1, 1, LinearAlgebra.dotProduct(x, yy).getDouble()); S2.set(1, 2, ((Number) xy.sum()).doubleValue()); S2.set(2, 0, S2.get(1, 1)); S2.set(2, 1, ((Number) yyy.sum()).doubleValue()); S2.set(2, 2, ((Number) yy.sum()).doubleValue()); Matrix S3 = new Matrix(3, 3); S3.set(0, 0, S2.get(0, 2)); S3.set(0, 1, S2.get(1, 2)); S3.set(0, 2, ((Number) x.sum()).doubleValue()); S3.set(1, 0, S3.get(0, 1)); S3.set(1, 1, S2.get(2, 2)); S3.set(1, 2, ((Number) y.sum()).doubleValue()); S3.set(2, 0, S3.get(0, 2)); S3.set(2, 1, S3.get(1, 2)); S3.set(2, 2, x.getSize()); Matrix T = S3.solve(S2.transpose()).uminus(); Matrix M = S1.plus(S2.times(T)); Matrix Cinv = new Matrix(new double[] { 0, 0, 0.5, 0, -1.0, 0, 0.5, 0, 0 }, 3); Matrix Mp = Cinv.times(M); // System.err.println("M " + Arrays.toString(Mp.getRowPackedCopy())); Matrix V = Mp.eig().getV(); // System.err.println("V " + Arrays.toString(V.getRowPackedCopy())); double[][] mv = V.getArray(); ArrayRealVector v1 = new ArrayRealVector(mv[0]); ArrayRealVector v2 = new ArrayRealVector(mv[1]); ArrayRealVector v3 = new ArrayRealVector(mv[2]); v1.mapMultiplyToSelf(4); ArrayRealVector v = v1.ebeMultiply(v3).subtract(v2.ebeMultiply(v2)); double[] varray = v.getDataRef(); int i = 0; for (; i < 3; i++) { if (varray[i] > 0) break; } if (i == 3) { throw new IllegalArgumentException("Could not find solution that satifies constraint"); } v = new ArrayRealVector(new double[] { mv[0][i], mv[1][i], mv[2][i] }); varray = v.getDataRef(); final double ca = varray[0]; final double cb = varray[1]; final double cc = varray[2]; Array2DRowRealMatrix mt = new Array2DRowRealMatrix(T.getArray(), false); varray = mt.operate(varray); final double cd = varray[0]; final double ce = varray[1]; final double cf = varray[2]; // System.err.println(String.format("Algebraic: %g, %g, %g, %g, %g, %g", ca, cb, cc, cd, ce, cf)); final double disc = cb * cb - 4. * ca * cc; if (disc >= 0) { throw new IllegalArgumentException("Solution is not an ellipse"); } if (cb == 0) { throw new IllegalArgumentException("Solution is a circle"); } double[] qparameters = new double[PARAMETERS]; qparameters[3] = (2. * cc * cd - cb * ce) / disc; qparameters[4] = (2. * ca * ce - cb * cd) / disc; final double sqrt = Math.sqrt((ca - cc) * (ca - cc) + cb * cb); qparameters[0] = -2. * (ca * ce * ce + cc * cd * cd + cf * cb * cb - cb * cd * ce - 4. * ca * cc * cf) / disc; qparameters[1] = qparameters[0] / (ca + cc + sqrt); qparameters[0] /= (ca + cc - sqrt); qparameters[0] = Math.sqrt(qparameters[0]); qparameters[1] = Math.sqrt(qparameters[1]); if (cb == 0) { qparameters[2] = 0.; } else { qparameters[2] = 0.5 * Math.atan2(cb, ca - cc); } if (qparameters[0] < qparameters[1]) { final double t = qparameters[0]; qparameters[0] = qparameters[1]; qparameters[1] = t; } else { qparameters[2] += Math.PI * 0.5; } return qparameters; }