List of usage examples for org.apache.commons.math.linear RealVector getEntry
double getEntry(int index);
From source file:fi.smaa.libror.MaximalVectorComputation.java
/** * Checks whether v1 dominates v2.//from w w w .j a v a 2 s .c o m * * @param v1 * @param v2 * @return true, if v1 \succ v2, false otherwise */ public static boolean dominates(RealVector v1, RealVector v2) { assert (v1.getDimension() == v2.getDimension()); boolean largerFound = false; for (int i = 0; i < v1.getDimension(); i++) { if (v1.getEntry(i) > v2.getEntry(i)) { largerFound = true; } else if (v1.getEntry(i) < v2.getEntry(i)) { return false; } } return largerFound; }
From source file:gda.images.camera.SampleMovementServiceBase.java
@Override public void moveSampleByMicronsAlongStandardAxes(double h, double v, double b) throws DeviceException { logger.debug(String.format("move in microns (standard axes): (h=%.2f, v=%.2f, b=%.2f)", h, v, b)); RealVector originalMovement = MatrixUtils.createRealVector(new double[] { h, v, b }); RealVector beamlineMovement = axisOrientationMatrix.operate(originalMovement); moveSampleByMicronsAlongBeamlineAxes(beamlineMovement.getEntry(0), beamlineMovement.getEntry(1), beamlineMovement.getEntry(2)); }
From source file:datafu.pig.hash.lsh.p_stable.AbstractStableDistributionFunction.java
/** * Compute the LSH for a given vector./*from ww w .j av a 2s . c o m*/ */ public long apply(RealVector vector) { /* * The hash is just floor(<v, a>/w) */ double ret = b; for (int i = 0; i < dim; ++i) { ret += vector.getEntry(i) * a[i]; } return (long) Math.floor(ret / w); }
From source file:com.opengamma.analytics.math.linearalgebra.SVDecompositionCommonsResultTest.java
private void assertRealVectorEquals(final RealVector v1, final DoubleMatrix1D v2) { final int n = v1.getDimension(); assertEquals(n, v2.getNumberOfElements()); for (int i = 0; i < n; i++) { assertEquals(v1.getEntry(i), v2.getEntry(i), EPS); }//from w ww .j a v a 2 s . c om }
From source file:geogebra.common.kernel.statistics.AlgoFitImplicit.java
private final void makeFunction() { int order = (int) orderGeo.evaluateDouble(); double[][] coeffs = new double[order + 1][order + 1]; // App.debug("row/cols = "+V.getRowDimension() + " "+ // V.getColumnDimension()+" "+(order * (order + 3) / 2 -1)); // App.debug(V.toString()); RealVector coeffsRV = V.getColumnVector(V.getColumnDimension() - 1); int c = 0;//from ww w . j a v a 2 s .c o m // create powers eg x^2y^0, x^1y^1, x^0*y^2, x, y, 1 for (int i = 0; i <= order; i++) { for (int j = 0; j <= i; j++) { coeffs[j][i - j] = coeffsRV.getEntry(c++); } } fitfunction.setCoeff(coeffs); fitfunction.setDefined(); }
From source file:eu.amidst.core.exponentialfamily.EF_Normal_NormalParents.java
/** * {@inheritDoc}/*from w w w . ja v a 2 s . c o m*/ */ @Override public void updateNaturalFromMomentParameters() { /* * First step: means and convariances */ CompoundVector globalMomentParam = (CompoundVector) this.momentParameters; double mean_X = globalMomentParam.getXYbaseMatrix().getEntry(0); RealVector mean_Y = globalMomentParam.getTheta_beta0BetaRV(); double cov_XX = globalMomentParam.getcovbaseMatrix().getEntry(0, 0) - mean_X * mean_X; RealMatrix cov_YY = globalMomentParam.getcovbaseMatrix().getSubMatrix(1, nOfParents, 1, nOfParents) .subtract(mean_Y.outerProduct(mean_Y)); RealVector cov_XY = globalMomentParam.getcovbaseMatrix().getSubMatrix(0, 0, 1, nOfParents).getRowVector(0) .subtract(mean_Y.mapMultiply(mean_X)); //RealVector cov_YX = cov_XY; //outerProduct transposes the vector automatically /* * Second step: betas and variance */ RealMatrix cov_YYInverse = new LUDecompositionImpl(cov_YY).getSolver().getInverse(); RealVector beta = cov_YYInverse.preMultiply(cov_XY); this.betas = new double[beta.getDimension()]; for (int i = 0; i < beta.getDimension(); i++) { this.betas[i] = beta.getEntry(i); } this.beta0 = mean_X - beta.dotProduct(mean_Y); this.variance = cov_XX - beta.dotProduct(cov_XY); }
From source file:net.ostis.scpdev.scg.geometry.GraphLayout.java
private void calculateForces() { int n = nodes.size(); RealVector[] forces = new RealVector[n]; for (int i = 0; i < n; ++i) forces[i] = new ArrayRealVector(new double[] { 0, 0 }); Map<IFigure, Integer> obj_f = new HashMap<IFigure, Integer>(); Point[] o_pos = new Point[n]; ////from w w w. j a v a2 s. c om // Calculation repulsion forces. // for (int idx = 0; idx < n; ++idx) { obj_f.put(nodes.get(idx), idx); Point p1 = nodes.get(idx).getBounds().getLocation(); RealVector p1v = new ArrayRealVector(2); p1v.setEntry(0, p1.x); p1v.setEntry(0, p1.y); double l = nullVector.getDistance(p1v); RealVector f = p1v.mapMultiply(gravity * (l - 3.0)); forces[idx].subtract(f); for (int jdx = idx + 1; jdx < n; ++jdx) { Point p2 = nodes.get(idx).getBounds().getLocation(); RealVector p2v = new ArrayRealVector(2); p2v.setEntry(0, p2.x); p2v.setEntry(0, p2.y); l = p1v.getDistance(p2v); if (l > max_rep_length) continue; if (l > 0.5) { f = p1v.subtract(p2v).mapMultiply(repulsion / l / l); } else { f = new ArrayRealVector(new double[] { Math.cos(0.17 * idx) * length * 7, Math.sin(0.17 * (idx + 1)) * length * 7 }); } forces[idx].add(f); forces[jdx].subtract(f); } } for (int idx = 0; idx < n; ++idx) { RealVector f = forces[idx]; f.mapMultiply(stepSize); nodes.get(idx).setLocation(new Point(f.getEntry(0), f.getEntry(1))); } }
From source file:net.ostis.scpdev.scg.geometry.LayoutRunnable.java
private void calculateForces() { int n = nodes.size(); RealVector[] forces = new RealVector[n]; for (int i = 0; i < n; ++i) forces[i] = new ArrayRealVector(nullVector); Map<IFigure, Integer> obj_f = new HashMap<IFigure, Integer>(); ////from w w w .j a v a 2s . c o m // Calculation repulsion forces. // for (int idx = 0; idx < n; ++idx) { SCgNodeShape obj = nodes.get(idx); obj_f.put(obj, idx); Point p1 = translateToCenter(obj.getBounds().getLocation()); RealVector p1v = new ArrayRealVector(2); p1v.setEntry(0, p1.x); p1v.setEntry(1, p1.y); double l = nullVector.getDistance(p1v); RealVector f = p1v.mapMultiply(gravity * (l - 3.0)); forces[idx] = forces[idx].subtract(f); for (int jdx = idx + 1; jdx < n; ++jdx) { Point p2 = translateToCenter(nodes.get(jdx).getBounds().getLocation()); RealVector p2v = new ArrayRealVector(2); p2v.setEntry(0, p2.x); p2v.setEntry(1, p2.y); l = p1v.getDistance(p2v) / 50; if (l > max_rep_length) continue; if (l > 0.5) { f = p1v.subtract(p2v).mapMultiply(repulsion / l / l); } else { f = new ArrayRealVector(new double[] { Math.cos(0.17 * idx) * length * 7, Math.sin(0.17 * (idx + 1)) * length * 7 }); } forces[idx] = forces[idx].add(f); forces[jdx] = forces[jdx].subtract(f); } } // // Calculation springs. // for (SCgPairConnection line : lines) { SCgPair pair = line.getModel(); SCgObject begin = pair.getBegin(); SCgObject end = pair.getEnd(); if (begin == null || end == null) continue; IFigure beginFigure = obj2figure.get(begin); IFigure endFigure = obj2figure.get(end); Point p1 = translateToCenter(beginFigure.getBounds().getLocation()); Point p2 = translateToCenter(endFigure.getBounds().getLocation()); RealVector p1v = new ArrayRealVector(2); p1v.setEntry(0, p1.x); p1v.setEntry(1, p1.y); RealVector p2v = new ArrayRealVector(2); p2v.setEntry(0, p2.x); p2v.setEntry(1, p2.y); double l = p1v.getDistance(p2v) / 50; RealVector f = null; if (l > 0) { RealVector pv = p2v.subtract(p1v); pv.unitize(); int cnt = begin.getInputCount() + end.getOutputCount(); f = pv.mapMultiply(rigidity * (l - Math.max(length, cnt / 3.0)) / l); if (nullVector.getDistance(f) > 10) f = pv.mapMultiply(rigidity / l); } else { f = new ArrayRealVector(nullVector); } if (obj_f.containsKey(beginFigure)) { int index = obj_f.get(beginFigure); forces[index] = forces[index].add(f); } if (obj_f.containsKey(endFigure)) { int index = obj_f.get(endFigure); forces[index] = forces[index].subtract(f); } } double maxf = 0.0; for (int idx = 0; idx < n; ++idx) { RealVector f = forces[idx]; f.mapMultiplyToSelf(stepSize); maxf = Math.max(maxf, nullVector.getDistance(f)); IFigure node = nodes.get(idx); Point location = translateToCenter(node.getBounds().getLocation()); location.x += f.getEntry(0); location.y += f.getEntry(1); node.setLocation(translateFromCenter(location)); } if (maxf > maxForce) { stepSize = stepMaxSize; } else { stepSize *= 0.97; } needLayout = stepSize > stepMinSize; }
From source file:org.caleydo.view.bicluster.elem.layout.MDSLayout.java
@Override public boolean doLayout(List<? extends IGLLayoutElement> children, float w_p, float h_p, IGLLayoutElement parent, int deltaTimeMs) { final int n = children.size(); RealMatrix a = MatrixUtils.createRealMatrix(n, n); // Sind als Startwerte anstatt Distanzen hnlichkeitsmae c_{ij} zwischen Objekten gegeben, so lassen sich diese // durch die Transformation ///* w w w . j av a 2 s . c o m*/ // d_{ij} = \sqrt{c_{ii}+c_{jj}-2c_{ij}} // // in Distanzen berfhren. for (int i = 0; i < n; ++i) { ClusterElement ci = (ClusterElement) children.get(i).asElement(); int c_ii = ci.getDimSize() + ci.getRecSize(); for (int j = i + 1; j < n; ++j) { ClusterElement cj = (ClusterElement) children.get(j).asElement(); int recOverlap = ci.getRecOverlap(cj); int dimOverlap = ci.getDimOverlap(cj); int c_jj = cj.getDimSize() + cj.getRecSize(); int c_ij = recOverlap + dimOverlap; double d_ij = Math.sqrt(c_ii + c_jj - 2 * c_ij); a.setEntry(i, j, d_ij); a.setEntry(j, i, d_ij); } } //#1. negative squared dissimilarity matrix Q //q = as.matrix( -0.5 * d ** 2 ) RealMatrix q = a.copy(); for (int i = 0; i < n; ++i) { q.getRowVector(i).mapPowToSelf(2); } q = q.scalarMultiply(-0.5); //#2. centering matrix H //h = diag(n) - 1/n * 1 RealMatrix h = MatrixUtils.createRealMatrix(n, n); for (int i = 0; i < n; ++i) h.setEntry(i, i, n - 1. / n * 1); //#3. double-center matrix B //b = h %*% q %*% h RealMatrix b = h.copy().multiply(q).multiply(h); // #4. eigen decomposition of B // eig = eigen(b) EigenDecompositionImpl eig = new EigenDecompositionImpl(b, 0); // #5. use top k values/vectors to compute projected points // points = eig$vectors[,1:k] %*% diag(sqrt(eig$values[1:k])) for (int i = 0; i < n; ++i) { RealVector v = eig.getEigenvector(i).getSubVector(0, 2); double x = v.getEntry(0) * eig.getRealEigenvalue(0); double y = v.getEntry(1) * eig.getRealEigenvalue(1); IGLLayoutElement child = children.get(i); child.setLocation((float) x, (float) y); } return false; }
From source file:rb.app.QNTransientRBSystem.java
/** * Perform online solve for current_params * with the N basis functions. Overload this * to solve the nonlinear RB system using * Newton's method./*from w ww. j a va 2 s . c o m*/ */ @Override public double RB_solve(int N) { current_N = N; _N_in_RB_solve = N; // Initialize tau_prod_k tau_prod_k = 1.; if (N > get_n_basis_functions()) { throw new RuntimeException( "ERROR: N cannot be larger than the number " + "of basis functions in RB_solve"); } if (N == 0) { throw new RuntimeException("ERROR: N must be greater than 0 in RB_solve"); } double dt = get_dt(); // First assemble the mass matrix RealMatrix RB_mass_matrix_N = new Array2DRowRealMatrix(N, N); for (int q_m = 0; q_m < get_Q_m(); q_m++) { RB_mass_matrix_N = RB_mass_matrix_N .add(RB_M_q_vector[q_m].getSubMatrix(0, N - 1, 0, N - 1).scalarMultiply(eval_theta_q_m(q_m))); } RealMatrix RB_RHS_Aq_matrix = new Array2DRowRealMatrix(N, N); RealMatrix Base_RB_LHS_matrix = RB_mass_matrix_N.scalarMultiply(1. / dt); for (int q_a = 0; q_a < get_Q_a(); q_a++) { double cached_theta_q_a = eval_theta_q_a(q_a); Base_RB_LHS_matrix = Base_RB_LHS_matrix.add(RB_A_q_vector[q_a].getSubMatrix(0, N - 1, 0, N - 1) .scalarMultiply(get_euler_theta() * cached_theta_q_a)); RB_RHS_Aq_matrix = RB_RHS_Aq_matrix .add(RB_A_q_vector[q_a].getSubMatrix(0, N - 1, 0, N - 1).scalarMultiply(-cached_theta_q_a)); } // Set system time level to 0 set_time_level(0); // Initialize a vector to store our current Newton iterate RealVector RB_u_bar = new ArrayRealVector(N); // and load the _k=0 data RB_solution = RB_u_bar; RB_temporal_solution_data[_k] = RB_u_bar; // Should use .copy() here! double error_bound_sum = 0.; // Set error bound at _k=0 error_bound_all_k[_k] = Math.sqrt(error_bound_sum); // Compute the outputs and associated error bounds at _k=0 for (int i = 0; i < get_n_outputs(); i++) { RB_outputs_all_k[i][_k] = 0.; RB_output_error_bounds_all_k[i][_k] = 0.; for (int q_l = 0; q_l < get_Q_l(i); q_l++) { RB_outputs_all_k[i][_k] += eval_theta_q_l(i, q_l) * (RB_output_vectors[i][q_l].getSubVector(0, N).dotProduct(RB_solution)); } RB_output_error_bounds_all_k[i][_k] = compute_output_dual_norm(i) * error_bound_all_k[_k]; } // Initialize a vector to store the solution from the old time-step RealVector RB_u_old = new ArrayRealVector(N); // Initialize a vector to store the Newton increment, RB_delta_u RealVector RB_delta_u = new ArrayRealVector(N); // Pre-compute eval_theta_c() double cached_theta_c = eval_theta_c(); for (int time_level = 1; time_level <= _K; time_level++) { set_time_level(time_level); // update the member variable _k // Set RB_u_old to be the result of the previous Newton loop RB_u_old = RB_u_bar.copy(); // Now we begin the nonlinear loop for (int l = 0; l < n_newton_steps; ++l) { // Get u_euler_theta = euler_theta*RB_u_bar + (1-euler_theta)*RB_u_old RealVector RB_u_euler_theta = RB_u_bar.mapMultiply(get_euler_theta()) .add(RB_u_old.mapMultiply(1. - get_euler_theta())); // Assemble the left-hand side for the RB linear system RealMatrix RB_LHS_matrix = Base_RB_LHS_matrix.copy(); // Add the trilinear term for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) { double new_entry = RB_LHS_matrix.getEntry(i, j); for (int n = 0; n < N; n++) { new_entry += cached_theta_c * get_euler_theta() * RB_u_euler_theta.getEntry(n) * (RB_trilinear_form[n][i][j] + RB_trilinear_form[j][i][n]); } RB_LHS_matrix.setEntry(i, j, new_entry); } } // Assemble the right-hand side for the RB linear system (the residual) // First add forcing terms RealVector RB_rhs = new ArrayRealVector(N); for (int q_f = 0; q_f < get_Q_f(); q_f++) { RB_rhs = RB_rhs.add(RB_F_q_vector[q_f].getSubVector(0, N).mapMultiply(eval_theta_q_f(q_f))); } // Now add -1./dt * M * (RB_u_bar - RB_u_old) RB_rhs = RB_rhs.add(RB_mass_matrix_N.operate(RB_u_bar).mapMultiply(-1. / dt)); RB_rhs = RB_rhs.add(RB_mass_matrix_N.operate(RB_u_old).mapMultiply(1. / dt)); // Now add the Aq stuff RB_rhs = RB_rhs.add(RB_RHS_Aq_matrix.operate(RB_u_euler_theta)); // Finally add the trilinear term for (int i = 0; i < N; i++) { double new_entry = RB_rhs.getEntry(i); for (int j = 0; j < N; j++) { double RB_u_euler_theta_j = RB_u_euler_theta.getEntry(j); for (int n = 0; n < N; n++) { new_entry -= cached_theta_c * RB_u_euler_theta.getEntry(n) * RB_u_euler_theta_j * RB_trilinear_form[n][i][j]; } } RB_rhs.setEntry(i, new_entry); } DecompositionSolver solver = new LUDecompositionImpl(RB_LHS_matrix).getSolver(); RB_delta_u = solver.solve(RB_rhs); // update the Newton iterate RB_u_bar = RB_u_bar.add(RB_delta_u); // Compute the l2 norm of RB_delta_u double RB_delta_u_norm = RB_delta_u.getNorm(); if (RB_delta_u_norm < nonlinear_tolerance) { break; } if ((l == (n_newton_steps - 1)) && (RB_delta_u_norm > nonlinear_tolerance)) { throw new RuntimeException("RB Newton loop did not converge"); } } // Load RB_solution into RB_solution_vector for residual computation RB_solution = RB_u_bar; old_RB_solution = RB_u_old; RB_temporal_solution_data[_k] = RB_u_bar; //should use copy here! double rho_LB = (mRbScmSystem == null) ? get_nominal_rho_LB() : get_SCM_lower_bound(); // Evaluate the dual norm of the residual for RB_solution_vector double epsilon_N = compute_residual_dual_norm(N); error_bound_sum += residual_scaling_numer(rho_LB) * Math.pow(epsilon_N, 2.); // store error bound at time-level _k error_bound_all_k[_k] = Math.sqrt(error_bound_sum / residual_scaling_denom(rho_LB)); // Now compute the outputs and associated error bounds for (int i = 0; i < get_n_outputs(); i++) { RB_outputs_all_k[i][_k] = 0.; RB_output_error_bounds_all_k[i][_k] = 0.; for (int q_l = 0; q_l < get_Q_l(i); q_l++) { RB_outputs_all_k[i][_k] += eval_theta_q_l(i, q_l) * (RB_output_vectors[i][q_l].getSubVector(0, N).dotProduct(RB_solution)); } RB_output_error_bounds_all_k[i][_k] = compute_output_dual_norm(i) * error_bound_all_k[_k]; } Log.d(DEBUG_TAG, "output = " + RB_outputs_all_k[0][_k] + ", bound=" + RB_output_error_bounds_all_k[0][_k]); } // Now compute the L2 norm of the RB solution at time-level _K // to normalize the error bound // We reuse RB_rhs here double final_RB_L2_norm = Math.sqrt(RB_mass_matrix_N.operate(RB_solution).dotProduct(RB_solution)); return (return_rel_error_bound ? error_bound_all_k[_K] / final_RB_L2_norm : error_bound_all_k[_K]); }