|
| 1 | +#include <common.h> |
| 2 | +#include <npe.h> |
| 3 | +#include <typedefs.h> |
| 4 | +#include <igl/active_set.h> |
| 5 | + |
| 6 | +const char *ds_active_set = R"igl_Qu8mg5v7( |
| 7 | +ACTIVE_SET Minimize quadratic energy |
| 8 | +
|
| 9 | +0.5*Z'*A*Z + Z'*B + C with constraints |
| 10 | +that Z(known) = Y, optionally also subject to the constraints Aeq*Z = Beq, |
| 11 | +and further optionally subject to the linear inequality constraints that |
| 12 | +Aieq*Z <= Bieq and constant inequality constraints lx <= x <= ux |
| 13 | +
|
| 14 | +Parameters |
| 15 | +---------- |
| 16 | +A n by n matrix of quadratic coefficients |
| 17 | +B n by 1 column of linear coefficients |
| 18 | +known list of indices to known rows in Z |
| 19 | +Y list of fixed values corresponding to known rows in Z |
| 20 | +Aeq meq by n list of linear equality constraint coefficients |
| 21 | +Beq meq by 1 list of linear equality constraint constant values |
| 22 | +Aieq mieq by n list of linear inequality constraint coefficients |
| 23 | +Bieq mieq by 1 list of linear inequality constraint constant values |
| 24 | +lx n by 1 list of lower bounds [] implies -Inf |
| 25 | +ux n by 1 list of upper bounds [] implies Inf |
| 26 | +params struct of additional parameters (see below) |
| 27 | +Z if not empty, is taken to be an n by 1 list of initial guess values (see output) |
| 28 | +Returns |
| 29 | +------- |
| 30 | +Z n by 1 list of solution values |
| 31 | +Returns SOLVER_STATUS_CONVERGED = 0, SOLVER_STATUS_MAX_ITER = 1, SOLVER_STATUS_ERROR = 2, |
| 32 | +
|
| 33 | +See also |
| 34 | +-------- |
| 35 | +
|
| 36 | +Notes |
| 37 | +----- |
| 38 | +For a harmonic solve on a mesh with 325K facets, matlab 2.2 secs, igl / min_quad_with_fixed.h 7.1 secs |
| 39 | +
|
| 40 | +Known Bugs : rows of[Aeq; Aieq] **must **be linearly independent.Should be using QR decomposition otherwise : http : //www.okstate.edu/sas/v8/sashtml/ormp/chap5/sect32.htm |
| 41 | +
|
| 42 | +
|
| 43 | +Examples |
| 44 | +-------- |
| 45 | +
|
| 46 | +)igl_Qu8mg5v7"; |
| 47 | + |
| 48 | +//NONe and dense_like |
| 49 | +npe_function(active_set) |
| 50 | +npe_doc(ds_active_set) |
| 51 | +npe_arg(A, sparse_float, sparse_double) |
| 52 | +npe_arg(B, dense_float, dense_double) |
| 53 | +npe_arg(known, dense_int, dense_long, dense_longlong) |
| 54 | +npe_arg(Y, npe_matches(B)) |
| 55 | +npe_arg(Aeq, npe_matches(A)) |
| 56 | +npe_arg(Beq, npe_matches(B)) |
| 57 | +npe_arg(Aieq, npe_matches(A)) |
| 58 | +npe_arg(Bieq, npe_matches(B)) |
| 59 | +npe_arg(lx, npe_matches(B)) |
| 60 | +npe_arg(ux, npe_matches(B)) |
| 61 | + |
| 62 | +npe_default_arg(Auu_pd, bool, false) |
| 63 | +npe_default_arg(max_iter, int, 100) |
| 64 | +npe_default_arg(inactive_threshold, double, igl::DOUBLE_EPS) |
| 65 | +npe_default_arg(constraint_threshold, double, igl::DOUBLE_EPS) |
| 66 | +npe_default_arg(solution_diff_threshold, double, igl::DOUBLE_EPS) |
| 67 | + |
| 68 | +npe_begin_code() |
| 69 | + assert_cols_equals(A, A.rows(), "A"); |
| 70 | + assert_rows_match(A, B, "A", "B"); |
| 71 | + assert_cols_equals(B, 1, "B"); |
| 72 | + |
| 73 | + assert_cols_match(A, Aeq, "A", "Aeq"); |
| 74 | + assert_rows_match(Aeq, Beq, "Aeq", "Beq"); |
| 75 | + assert_cols_equals(Beq, 1, "Beq"); |
| 76 | + |
| 77 | + assert_rows_match(Aieq, Bieq, "Aieq", "Bieq"); |
| 78 | + |
| 79 | + if(Aieq.size() > 0) |
| 80 | + { |
| 81 | + assert_cols_match(A, Aieq, "A", "Aieq"); |
| 82 | + assert_cols_equals(Bieq, 1, "Bieq"); |
| 83 | + } |
| 84 | + |
| 85 | + if(lx.size() > 0) |
| 86 | + { |
| 87 | + assert_rows_match(A, lx, "A", "lx"); |
| 88 | + assert_cols_equals(lx, 1, "lx"); |
| 89 | + } |
| 90 | + |
| 91 | + if (ux.size() > 0) |
| 92 | + { |
| 93 | + assert_rows_match(A, ux, "A", "ux"); |
| 94 | + assert_cols_equals(ux, 1, "ux"); |
| 95 | + } |
| 96 | + |
| 97 | + Eigen::SparseMatrix<double> A_copy = A.template cast<double>(); |
| 98 | + Eigen::VectorXd B_copy = B.template cast<double>(); |
| 99 | + Eigen::VectorXi known_copy = known.template cast<int>(); |
| 100 | + Eigen::VectorXd Y_copy = Y.template cast<double>(); |
| 101 | + Eigen::SparseMatrix<double> Aeq_copy = Aeq.template cast<double>(); |
| 102 | + Eigen::VectorXd Beq_copy = Beq.template cast<double>(); |
| 103 | + Eigen::SparseMatrix<double> Aieq_copy; |
| 104 | + Eigen::VectorXd Bieq_copy; |
| 105 | + |
| 106 | + if(Aieq.size() > 0) |
| 107 | + { |
| 108 | + Aieq_copy = Aieq.template cast<double>(); |
| 109 | + Bieq_copy = Bieq.template cast<double>(); |
| 110 | + } |
| 111 | + |
| 112 | + Eigen::VectorXd lx_copy; |
| 113 | + if(lx.size() > 0) |
| 114 | + lx_copy = lx.template cast<double>(); |
| 115 | + Eigen::VectorXd ux_copy; |
| 116 | + if (ux.size() > 0) |
| 117 | + ux_copy = ux.template cast<double>(); |
| 118 | + |
| 119 | + igl::active_set_params params; |
| 120 | + params.Auu_pd = Auu_pd; |
| 121 | + params.max_iter = max_iter; |
| 122 | + params.inactive_threshold = inactive_threshold; |
| 123 | + params.constraint_threshold = constraint_threshold; |
| 124 | + params.solution_diff_threshold = solution_diff_threshold; |
| 125 | + |
| 126 | + Eigen::VectorXd Z; |
| 127 | + |
| 128 | + auto res = igl::active_set(A_copy, B_copy, known_copy, Y_copy, Aeq_copy, Beq_copy, Aieq_copy, Bieq_copy, lx_copy, ux_copy, params, Z); |
| 129 | + EigenDenseLike<npe_Matrix_B> Z_row_major = Z.template cast<typename npe_Matrix_B::Scalar>(); |
| 130 | + |
| 131 | + return std::make_tuple(int(res), npe::move(Z_row_major)); |
| 132 | + |
| 133 | + npe_end_code() |
0 commit comments