List of usage examples for org.apache.commons.math.linear RealMatrix copy
RealMatrix copy();
From source file:juicebox.tools.utils.juicer.apa.APAUtils.java
public static RealMatrix standardNormalization(RealMatrix matrix) { return matrix.copy().scalarMultiply(1. / Math.max(1., MatrixTools.mean(matrix))); }
From source file:juicebox.tools.utils.juicer.apa.APAUtils.java
public static RealMatrix centerNormalization(RealMatrix matrix) { int center = matrix.getRowDimension() / 2; double centerVal = matrix.getEntry(center, center); if (centerVal == 0) { centerVal = MatrixTools.minimumPositive(matrix); if (centerVal == 0) centerVal = 1;/* w w w .ja v a 2 s .co m*/ } return matrix.copy().scalarMultiply(1. / centerVal); }
From source file:math2605.qr_fact_givens.java
public qr_fact_givens(RealMatrix A) { int m = A.getRowDimension(); int n = A.getColumnDimension(); R = A.copy(); RealMatrix a = A.copy();//from ww w .j a v a 2 s . c o m Q = null; boolean firstQ = true; for (int j = 0; j < n; j++) { int row = j; for (int i = j; i < m - 1; i++) { if (a.getEntry(row, j) == 0) { //System.out.println("in break"); break; } if (upperTriangular(a)) { //System.out.println("is ut"); continue; } double[] cs = new double[2]; cs = givensRoatation(a.getEntry(i, j), a.getEntry(i + 1, j)); RealMatrix G = createG(i, i + 1, cs, n); R = G.multiply(R); if (firstQ) { Q = G.transpose(); firstQ = false; } else { Q = Q.multiply(G.transpose()); } a = G.multiply(a); row++; } } }
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 ////from w ww. j a va 2s . 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./* w w w . j a va 2 s . com*/ */ @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]); }
From source file:tools.help.java
/** * initialize matrix with zero values based on size of another matrix * @param source_mtx matrix defining the size of the initiating matrix * @return matrix with zero values with size of source_mtx *///from w w w .java2 s .c o m public static RealMatrix initMatrix(RealMatrix source_mtx) { RealMatrix final_mtx = source_mtx.copy(); return clearMatrix(final_mtx); }