List of usage examples for org.apache.commons.math.linear RealVector add
RealVector add(double[] v);
From source file:gda.images.camera.SampleMovementServiceBase.java
@Override public void moveSampleByMicronsAlongBeamlineAxes(double x, double y, double z) throws DeviceException { logger.debug(String.format("move in microns (beamline axes): (x=%.2f, y=%.2f, z=%.2f)", x, y, z)); double[] currentPos = getPosition(); RealVector currentPosVector = MatrixUtils.createRealVector(currentPos); logger.debug(String.format("current position is %s", currentPosVector)); RealVector move = MatrixUtils.createRealVector(new double[] { x, y, z }); RealVector newPosition = currentPosVector.add(move); logger.debug(String.format("new position is %s", newPosition)); setPosition(newPosition.getData()); }
From source file:gda.scan.TrajectoryScanLine.java
List<Map<Scannable, double[]>> generateTrajectoryForDetectorsThatIntegrateBetweenTriggers() { List<Map<Scannable, double[]>> triggers = new LinkedList<Map<Scannable, double[]>>(); List<Map<Scannable, RealVector>> binCentres = scanPositionRecorder.getPoints(); HashMap<Scannable, double[]> pointToAdd; Set<Scannable> scannables = binCentres.get(0).keySet(); // Add first trigger: xtrig[0] = bincentre[0] - (bincentre[1] - bincentre[0]) / 2 pointToAdd = new HashMap<Scannable, double[]>(); for (Scannable scannable : scannables) { RealVector first = binCentres.get(0).get(scannable); RealVector second = binCentres.get(1).get(scannable); double[] firstTrigger = first.subtract(second.subtract(first).mapDivide(2.)).toArray(); pointToAdd.put(scannable, firstTrigger); }/* w ww. j ava 2 s.c o m*/ triggers.add(pointToAdd); // Add middle triggers: xtrig[i] = (bincentre[i] + bincentre[i+1]) / 2 for (int i = 0; i < binCentres.size() - 1; i++) { // not the last one pointToAdd = new HashMap<Scannable, double[]>(); for (Scannable scannable : scannables) { RealVector current = binCentres.get(i).get(scannable); RealVector next = binCentres.get(i + 1).get(scannable); double[] trigger = current.add(next).mapDivide(2.).toArray(); pointToAdd.put(scannable, trigger); } triggers.add(pointToAdd); } // Add last trigger: xtrig[n+1] = bincentre[n] + (bincentre[n] - bincentre[n-1]) / 2 pointToAdd = new HashMap<Scannable, double[]>(); for (Scannable scannable : scannables) { int lastIndex = binCentres.size() - 1; RealVector last = binCentres.get(lastIndex).get(scannable); RealVector secondLast = binCentres.get(lastIndex - 1).get(scannable); double[] lastTrigger = last.add(last.subtract(secondLast).mapDivide(2.)).toArray(); pointToAdd.put(scannable, lastTrigger); } triggers.add(pointToAdd); return triggers; }
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 w w .j av a 2s .co 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]); }
From source file:rb.app.RBSystem.java
/** * Perform online solve with the N RB basis functions, for the set of * parameters in current_params, where 1 <= N <= RB_size. *//*from ww w . j a v a 2 s.co m*/ public double RB_solve(int N) { current_N = N; 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"); } // Assemble the RB system RealMatrix RB_system_matrix_N = new Array2DRowRealMatrix(N, N); for (int q_a = 0; q_a < get_Q_a(); q_a++) { RB_system_matrix_N = RB_system_matrix_N .add(RB_A_q_vector[q_a].getSubMatrix(0, N - 1, 0, N - 1).scalarMultiply(eval_theta_q_a(q_a))); } // Assemble the RB rhs RealVector RB_rhs_N = new ArrayRealVector(N); for (int q_f = 0; q_f < get_Q_f(); q_f++) { // Note getSubVector takes an initial index and the number of entries // i.e. the interface is a bit different to getSubMatrix RB_rhs_N = RB_rhs_N.add(RB_F_q_vector[q_f].getSubVector(0, N).mapMultiply(eval_theta_q_f(q_f))); } // Solve the linear system DecompositionSolver solver = new LUDecompositionImpl(RB_system_matrix_N).getSolver(); RB_solution = solver.solve(RB_rhs_N); // Evaluate the dual norm of the residual for RB_solution_vector double epsilon_N = compute_residual_dual_norm(N); // Get lower bound for coercivity constant double alpha_LB = get_SCM_lower_bound(); // If SCM lower bound is negative if (alpha_LB <= 0) { // Get an upper bound instead alpha_LB = get_SCM_upper_bound(); } // Store (absolute) error bound double abs_error_bound = epsilon_N / residual_scaling_denom(alpha_LB); // Compute the norm of RB_solution double RB_solution_norm = RB_solution.getNorm(); // Now compute the outputs and associated errors RealVector RB_output_vector_N = new ArrayRealVector(N); for (int i = 0; i < get_n_outputs(); i++) { RB_outputs[i] = 0.; for (int q_l = 0; q_l < get_Q_l(i); q_l++) { RB_output_vector_N = RB_output_vectors[i][q_l].getSubVector(0, N); RB_outputs[i] += eval_theta_q_l(i, q_l) * RB_solution.dotProduct(RB_output_vector_N); } RB_output_error_bounds[i] = compute_output_dual_norm(i) * abs_error_bound; } return (return_rel_error_bound ? abs_error_bound / RB_solution_norm : abs_error_bound); }
From source file:rb.app.TransientRBSystem.java
/** * Perform online solve with the N RB basis functions, for the set of * parameters in current_params, where 1 <= N <= RB_size. *///ww w . ja v a 2s . co m @Override public double RB_solve(int N) { current_N = N; 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"); } // 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_LHS_matrix = new Array2DRowRealMatrix(N, N); RealMatrix RB_RHS_matrix = new Array2DRowRealMatrix(N, N); RB_LHS_matrix = RB_LHS_matrix.add(RB_mass_matrix_N.scalarMultiply(1. / get_dt())); RB_RHS_matrix = RB_RHS_matrix.add(RB_mass_matrix_N.scalarMultiply(1. / get_dt())); for (int q_a = 0; q_a < get_Q_a(); q_a++) { RB_LHS_matrix = RB_LHS_matrix.add(RB_A_q_vector[q_a].getSubMatrix(0, N - 1, 0, N - 1) .scalarMultiply(get_euler_theta() * eval_theta_q_a(q_a))); RB_RHS_matrix = RB_RHS_matrix.add(RB_A_q_vector[q_a].getSubMatrix(0, N - 1, 0, N - 1) .scalarMultiply(-(1. - get_euler_theta()) * eval_theta_q_a(q_a))); } set_time_level(0); // Sets the member variable _k to zero // Zero initial condition RealVector RB_solution_N = new ArrayRealVector(N); RealVector old_RB_solution_N = new ArrayRealVector(N); // Initialize rhs RealVector RB_rhs_N = new ArrayRealVector(N); // Load the initial condition into RB_temporal_solution_data RB_temporal_solution_data[_k] = RB_solution; 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_N)); } RB_output_error_bounds_all_k[i][_k] = compute_output_dual_norm(i) * error_bound_all_k[_k]; } double alpha_LB = get_SCM_lower_bound(); // Precompute time-invariant parts of the dual norm of the residual. cache_online_residual_terms(N); for (int time_level = 1; time_level <= _K; time_level++) { set_time_level(time_level); // This updates the member variable _k old_RB_solution_N = RB_solution_N; // Compute RB_rhs, as RB_LHS_matrix x old_RB_solution RB_rhs_N = RB_RHS_matrix.operate(old_RB_solution_N); // Add forcing terms for (int q_f = 0; q_f < get_Q_f(); q_f++) { RB_rhs_N = RB_rhs_N.add(RB_F_q_vector[q_f].getSubVector(0, N).mapMultiply(eval_theta_q_f(q_f))); } // Solve the linear system DecompositionSolver solver = new LUDecompositionImpl(RB_LHS_matrix).getSolver(); RB_solution_N = solver.solve(RB_rhs_N); // Save RB_solution for current time level RB_temporal_solution_data[_k] = RB_solution_N; // Evaluate the dual norm of the residual for RB_solution_vector RB_solution = RB_solution_N; old_RB_solution = old_RB_solution_N; double epsilon_N = compute_residual_dual_norm(N); error_bound_sum += residual_scaling_numer(alpha_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(alpha_LB)); // Now compute the outputs and associated errors 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_N)); } RB_output_error_bounds_all_k[i][_k] = compute_output_dual_norm(i) * error_bound_all_k[_k]; } } // Now compute the L2 norm of the RB solution at time-level _K // to normalize the error bound // We reuse RB_rhs here RealMatrix RB_L2_matrix_N = RB_L2_matrix.getSubMatrix(0, N - 1, 0, N - 1); double final_RB_L2_norm = Math.sqrt(RB_L2_matrix_N.operate(RB_solution_N).dotProduct(RB_solution_N)); if (apply_temporal_filter_flag) { apply_temporal_filter(); } return (return_rel_error_bound ? error_bound_all_k[_K] / final_RB_L2_norm : error_bound_all_k[_K]); }