Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Solve quadprog #1896

Closed
wants to merge 14 commits into from
Prev Previous commit
Next Next commit
any arg can be a parameter
  • Loading branch information
rok-cesnovar committed Apr 18, 2020
commit 4c3a855dd6262e24eda88c0a4e18d0f6361d2a2b
26 changes: 14 additions & 12 deletions stan/math/rev/fun/solve_quadprog.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -277,14 +277,14 @@ solve_quadprog(const Eigen::Matrix<T0, -1, -1> &G,
int n = g0.size();
int p = ce0.size();
int m = ci0.size();
stan::math::matrix_v R(G.rows(), G.cols());

stan::math::vector_v np(n), u(m + p), d(n), z(n), r(m + p), s(m + p);
stan::math::vector_v x_old(n), u_old(m + p);

using T_return = return_type_t<T0, T1, T2, T3, T4, T5>;
Eigen::Matrix<T_return, -1, 1> np(n), u(m + p), d(n), z(n), r(m + p), s(m + p), u_old(m + p);
Eigen::Matrix<T_return, -1, 1> x_old(n);
const double inf = std::numeric_limits<double>::infinity();
var c2, sum, psi, ss, t1;
/* t is the step length, which is the minimum of the partial
* step length t1 and the full step length t2 */

T_return sum, psi, ss, t1;

Eigen::VectorXi A(m + p), A_old(m + p), iai(m + p);
int q;
int iq, iter = 0;
Expand All @@ -307,14 +307,16 @@ solve_quadprog(const Eigen::Matrix<T0, -1, -1> &G,

/* initialize the matrix R */
d.setZero();

Eigen::Matrix<T0, -1, -1> R(G.rows(), G.cols());
R.setZero();
var R_norm = 1.0; /* this variable will hold the norm of the matrix R */
T_return R_norm = 1.0; /* this variable will hold the norm of the matrix R */

/* compute the inverse of the factorized matrix G^-1, this is the initial
* value for H */
//J = L^-T
auto J = stan::math::transpose(stan::math::mdivide_left_tri_low(chol));
c2 = J.trace();
auto c2 = J.trace();
#ifdef TRACE_SOLVER
print_matrix("J", J, n);
#endif
Expand All @@ -337,8 +339,8 @@ solve_quadprog(const Eigen::Matrix<T0, -1, -1> &G,
#endif
/* Add equality constraints to the working set A */
iq = 0;
var t2 = 0.0;
var t = 0.0;
T_return t2 = 0.0;
T_return t = 0.0;
for (i = 0; i < me; i++) {
np = CE.col(i);
compute_d(d, J, np);
Expand Down Expand Up @@ -467,7 +469,7 @@ solve_quadprog(const Eigen::Matrix<T0, -1, -1> &G,
t1 = inf; /* +inf */
/* find the index l s.t. it reaches the minimum of u+(x) / r */
for (k = me; k < iq; k++) {
var tmp = u(k) / r(k);
T_return tmp = u(k) / r(k);
if (r(k) > 0.0 && (tmp < t1)) {
t1 = tmp;
l = A(k);
Expand Down
172 changes: 170 additions & 2 deletions test/unit/math/rev/fun/solve_quadprog_test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,35 @@
#include <test/unit/math/rev/fun/util.hpp>
#include <test/unit/math/rev/util.hpp>

TEST(AgradRevMatrix, multiply_val_vv_cl) {
TEST(AgradRevMatrix, dddddd) {
stan::math::matrix_d G(3,3);
stan::math::vector_d g0(3);
Eigen::MatrixXd CE(3,1);
Eigen::VectorXd ce0(1);
Eigen::MatrixXd CI(3,4);
Eigen::VectorXd ci0(4);

G << 2.1, 1.5, 1.2,
1.5, 2.2, 1.3,
1.2, 1.3, 3.1;

g0 << 6, 1, 1;

CE << 1, 2, -1;

ce0(0)=-4;

CI << 1, 0, 0, -1,
0, 1, 0, -1,
0, 0, 1, 0;

ci0 << 0, 0, 0, 10;


stan::math::vector_d x = stan::math::solve_quadprog(G, g0, CE, ce0, CI, ci0);
}

TEST(AgradRevMatrix, vddddd) {
stan::math::matrix_v G(3,3);
stan::math::vector_d g0(3);
Eigen::MatrixXd CE(3,1);
Expand All @@ -28,5 +56,145 @@ TEST(AgradRevMatrix, multiply_val_vv_cl) {
ci0 << 0, 0, 0, 10;


stan::math::vector_v x = solve_quadprog(G, g0, CE, ce0, CI, ci0);
stan::math::vector_v x = stan::math::solve_quadprog(G, g0, CE, ce0, CI, ci0);
}

TEST(AgradRevMatrix, vvdddd) {
stan::math::matrix_v G(3,3);
stan::math::vector_v g0(3);
Eigen::MatrixXd CE(3,1);
Eigen::VectorXd ce0(1);
Eigen::MatrixXd CI(3,4);
Eigen::VectorXd ci0(4);

G << 2.1, 1.5, 1.2,
1.5, 2.2, 1.3,
1.2, 1.3, 3.1;

g0 << 6, 1, 1;

CE << 1, 2, -1;

ce0(0)=-4;

CI << 1, 0, 0, -1,
0, 1, 0, -1,
0, 0, 1, 0;

ci0 << 0, 0, 0, 10;


stan::math::vector_v x = stan::math::solve_quadprog(G, g0, CE, ce0, CI, ci0);
}

TEST(AgradRevMatrix, vvvddd) {
stan::math::matrix_v G(3,3);
stan::math::vector_v g0(3);
stan::math::matrix_v CE(3,1);
Eigen::VectorXd ce0(1);
Eigen::MatrixXd CI(3,4);
Eigen::VectorXd ci0(4);

G << 2.1, 1.5, 1.2,
1.5, 2.2, 1.3,
1.2, 1.3, 3.1;

g0 << 6, 1, 1;

CE << 1, 2, -1;

ce0(0)=-4;

CI << 1, 0, 0, -1,
0, 1, 0, -1,
0, 0, 1, 0;

ci0 << 0, 0, 0, 10;


stan::math::vector_v x = stan::math::solve_quadprog(G, g0, CE, ce0, CI, ci0);
}

TEST(AgradRevMatrix, vvvvdd) {
stan::math::matrix_v G(3,3);
stan::math::vector_v g0(3);
stan::math::matrix_v CE(3,1);
stan::math::vector_v ce0(1);
Eigen::MatrixXd CI(3,4);
Eigen::VectorXd ci0(4);

G << 2.1, 1.5, 1.2,
1.5, 2.2, 1.3,
1.2, 1.3, 3.1;

g0 << 6, 1, 1;

CE << 1, 2, -1;

ce0(0)=-4;

CI << 1, 0, 0, -1,
0, 1, 0, -1,
0, 0, 1, 0;

ci0 << 0, 0, 0, 10;


stan::math::vector_v x = stan::math::solve_quadprog(G, g0, CE, ce0, CI, ci0);
}

TEST(AgradRevMatrix, vvvvvd) {
stan::math::matrix_v G(3,3);
stan::math::vector_v g0(3);
stan::math::matrix_v CE(3,1);
stan::math::vector_v ce0(1);
stan::math::matrix_v CI(3,4);
Eigen::VectorXd ci0(4);

G << 2.1, 1.5, 1.2,
1.5, 2.2, 1.3,
1.2, 1.3, 3.1;

g0 << 6, 1, 1;

CE << 1, 2, -1;

ce0(0)=-4;

CI << 1, 0, 0, -1,
0, 1, 0, -1,
0, 0, 1, 0;

ci0 << 0, 0, 0, 10;


stan::math::vector_v x = stan::math::solve_quadprog(G, g0, CE, ce0, CI, ci0);
}

TEST(AgradRevMatrix, vvvvvv) {
stan::math::matrix_v G(3,3);
stan::math::vector_v g0(3);
stan::math::matrix_v CE(3,1);
stan::math::vector_v ce0(1);
stan::math::matrix_v CI(3,4);
stan::math::vector_v ci0(4);

G << 2.1, 1.5, 1.2,
1.5, 2.2, 1.3,
1.2, 1.3, 3.1;

g0 << 6, 1, 1;

CE << 1, 2, -1;

ce0(0)=-4;

CI << 1, 0, 0, -1,
0, 1, 0, -1,
0, 0, 1, 0;

ci0 << 0, 0, 0, 10;


stan::math::vector_v x = stan::math::solve_quadprog(G, g0, CE, ce0, CI, ci0);
}