diff --git a/include/pyoptinterface/knitro_model.hpp b/include/pyoptinterface/knitro_model.hpp index 18d82df..46b81ae 100644 --- a/include/pyoptinterface/knitro_model.hpp +++ b/include/pyoptinterface/knitro_model.hpp @@ -128,64 +128,95 @@ enum ConstraintSenseFlags CON_UPBND = 1 << 1, // 0x02 }; +template struct CallbackPattern { - std::vector indexCons; - std::vector objGradIndexVars; - std::vector jacIndexCons; - std::vector jacIndexVars; - std::vector hessIndexVars1; - std::vector hessIndexVars2; + std::vector indexCons; + std::vector objGradIndexVars; + std::vector jacIndexCons; + std::vector jacIndexVars; + std::vector hessIndexVars1; + std::vector hessIndexVars2; }; -template +using namespace CppAD; + +template struct CallbackEvaluator { - static inline const std::string jac_coloring_ = "cppad"; - static inline const std::string hess_coloring_ = "cppad.symmetric"; - std::vector indexVars; - std::vector indexCons; - - CppAD::ADFun fun; - CppAD::sparse_rc> jac_pattern_; - CppAD::sparse_rcv, std::vector> jac_; - CppAD::sparse_jac_work jac_work_; - CppAD::sparse_rc> hess_pattern_; - CppAD::sparse_rc> hess_pattern_symm_; - CppAD::sparse_rcv, std::vector> hess_; - CppAD::sparse_hes_work hess_work_; - - std::vector x; - std::vector w; + + static inline constexpr const char *CLRNG = "cppad"; + + std::vector indexVars; + std::vector indexCons; + + ADFun fun; /// < CppAD tape. + ADFun jfun; /// < CppAD tape for Aggregated Jacobian + + /// Sparsity patterns + sparse_rc> jp; + sparse_rc> hp; + + /// Workspaces for Jacobian and Hessian calculations + sparse_jac_work jw; + sparse_jac_work hw; + + /// Temporary vectors for evaluations + vector x; + vector xw; + sparse_rcv, vector> jac; + sparse_rcv, vector> hes; void setup() { fun.optimize(); - auto nx = fun.Domain(); - auto ny = fun.Range(); - CppAD::sparse_rc> jac_pattern_in(nx, nx, nx); + size_t nx = fun.Domain(); + size_t ny = fun.Range(); + + vector dom(nx, true); + vector rng(ny, true); + fun.subgraph_sparsity(dom, rng, false, jp); + + ADFun, V> af = fun.base2ad(); + vector> jaxw(nx + ny); + Independent(jaxw); + vector> jax(nx); + vector> jaw(ny); + vector> jaz(nx); for (size_t i = 0; i < nx; i++) { - jac_pattern_in.set(i, i, i); + jax[i] = jaxw[i]; } - fun.for_jac_sparsity(jac_pattern_in, false, false, false, jac_pattern_); - std::vector select_rows(ny, true); - fun.rev_hes_sparsity(select_rows, false, false, hess_pattern_); - auto &hess_rows = hess_pattern_.row(); - auto &hess_cols = hess_pattern_.col(); - for (size_t k = 0; k < hess_pattern_.nnz(); k++) + for (size_t i = 0; i < ny; i++) { - size_t row = hess_rows[k]; - size_t col = hess_cols[k]; - if (row <= col) + jaw[i] = jaxw[nx + i]; + } + af.Forward(0, jax); + jaz = af.Reverse(1, jaw); + jfun.Dependent(jaxw, jaz); + jfun.optimize(); + vector jdom(nx + ny, false); + for (size_t i = 0; i < nx; i++) + { + jdom[i] = true; + } + vector jrng(nx, true); + sparse_rc> hsp; + jfun.subgraph_sparsity(jdom, jrng, false, hsp); + + auto &hrow = hsp.row(); + auto &hcol = hsp.col(); + for (size_t k = 0; k < hsp.nnz(); k++) + { + if (hrow[k] <= hcol[k]) { - hess_pattern_symm_.push_back(row, col); + hp.push_back(hrow[k], hcol[k]); } } - x.resize(nx, 0.0); - w.resize(ny, 0.0); - jac_ = CppAD::sparse_rcv, std::vector>(jac_pattern_); - hess_ = CppAD::sparse_rcv, std::vector>(hess_pattern_symm_); + x.resize(nx); + xw.resize(nx + ny); + jac = sparse_rcv, vector>(jp); + hes = sparse_rcv, vector>(hp); } bool is_objective() const @@ -195,107 +226,108 @@ struct CallbackEvaluator void eval_fun(const V *req_x, V *res_y) { - copy_ptr(req_x, indexVars.data(), x); + copy(fun.Domain(), req_x, indexVars.data(), x.data()); auto y = fun.Forward(0, x); - copy_vec(y, res_y, is_objective()); + int mode = is_objective() ? 2 : 0; + copy(fun.Range(), y.data(), (const I *)nullptr, res_y, mode); } void eval_jac(const V *req_x, V *res_jac) { - copy_ptr(req_x, indexVars.data(), x); - fun.sparse_jac_rev(x, jac_, jac_pattern_, jac_coloring_, jac_work_); - auto &jac = jac_.val(); - copy_vec(jac, res_jac); + copy(fun.Domain(), req_x, indexVars.data(), x.data()); + fun.sparse_jac_rev(x, jac, jp, CLRNG, jw); + copy(jac.nnz(), jac.val().data(), (const I *)nullptr, res_jac); } void eval_hess(const V *req_x, const V *req_w, V *res_hess) { - copy_ptr(req_x, indexVars.data(), x); - copy_ptr(req_w, indexCons.data(), w, is_objective()); - fun.sparse_hes(x, w, hess_, hess_pattern_, hess_coloring_, hess_work_); - auto &hess = hess_.val(); - copy_vec(hess, res_hess); + copy(fun.Domain(), req_x, indexVars.data(), xw.data()); + int mode = is_objective() ? 1 : 0; + copy(fun.Range(), req_w, indexCons.data(), xw.data() + fun.Domain(), mode); + jfun.sparse_jac_rev(xw, hes, hp, CLRNG, hw); + copy(hes.nnz(), hes.val().data(), (const I *)nullptr, res_hess); } - CallbackPattern get_callback_pattern() const + CallbackPattern get_callback_pattern() const { - CallbackPattern pattern; - pattern.indexCons = indexCons; + CallbackPattern p; + p.indexCons = indexCons; - auto &jac_rows = jac_pattern_.row(); - auto &jac_cols = jac_pattern_.col(); + auto &jrow = jp.row(); + auto &jcol = jp.col(); if (indexCons.empty()) { - for (size_t k = 0; k < jac_pattern_.nnz(); k++) + for (size_t k = 0; k < jp.nnz(); k++) { - pattern.objGradIndexVars.push_back(indexVars[jac_cols[k]]); + p.objGradIndexVars.push_back(indexVars[jcol[k]]); } } else { - for (size_t k = 0; k < jac_pattern_.nnz(); k++) + for (size_t k = 0; k < jp.nnz(); k++) { - pattern.jacIndexCons.push_back(indexCons[jac_rows[k]]); - pattern.jacIndexVars.push_back(indexVars[jac_cols[k]]); + p.jacIndexCons.push_back(indexCons[jrow[k]]); + p.jacIndexVars.push_back(indexVars[jcol[k]]); } } - auto &hess_rows = hess_pattern_symm_.row(); - auto &hess_cols = hess_pattern_symm_.col(); - for (size_t k = 0; k < hess_pattern_symm_.nnz(); k++) + auto &hrow = hp.row(); + auto &hcol = hp.col(); + for (size_t k = 0; k < hp.nnz(); k++) { - pattern.hessIndexVars1.push_back(indexVars[hess_rows[k]]); - pattern.hessIndexVars2.push_back(indexVars[hess_cols[k]]); + p.hessIndexVars1.push_back(indexVars[hrow[k]]); + p.hessIndexVars2.push_back(indexVars[hcol[k]]); } - return pattern; + return p; } private: - template - static void copy_ptr(const T *src, const I *idx, std::vector &dst, bool duplicate = false) + // Copy mode: + // - 0: normal copy + // - 1: duplicate (copy first element of src to all elements of dst) + // - 2: aggregate (sum all elements of src and copy to all elements of dst) + static void copy(const size_t n, const V *src, const I *idx, V *dst, int mode = 0) { - for (size_t i = 0; i < dst.size(); i++) + if (mode == 1) { - if (duplicate) + for (size_t i = 0; i < n; i++) { dst[i] = src[0]; } - else - { - dst[i] = src[idx[i]]; - } } - } - - template - static void copy_vec(const std::vector &src, T *dst, bool aggregate = false) - { - if (aggregate) + else if (mode == 2) { - dst[0] = 0.0; + if (n == 0) + { + return; + } + dst[0] = src[0]; + for (size_t i = 1; i < n; i++) + { + dst[0] += src[i]; + } } - for (size_t i = 0; i < src.size(); i++) + else { - if (aggregate) + if (idx == nullptr) { - dst[0] += src[i]; + for (size_t i = 0; i < n; i++) + { + dst[i] = src[i]; + } } else { - dst[i] = src[i]; + for (size_t i = 0; i < n; i++) + { + dst[i] = src[idx[i]]; + } } } } }; -struct Outputs -{ - std::vector objective_outputs; - std::vector constraint_outputs; - std::vector constraints; -}; - inline bool is_name_empty(const char *name) { return name == nullptr || name[0] == '\0'; @@ -577,8 +609,17 @@ class KNITROModel : public OnesideLinearConstraintMixin, std::unordered_map m_con_sense_flags; uint8_t m_obj_flag = 0; + struct Outputs + { + std::vector objective_outputs; + std::vector constraint_outputs; + std::vector constraints; + }; + + using Evaluator = CallbackEvaluator; + std::unordered_map m_pending_outputs; - std::vector>> m_evaluators; + std::vector> m_evaluators; bool m_has_pending_callbacks = false; int m_solve_status = 0; bool m_is_dirty = true; @@ -604,7 +645,7 @@ class KNITROModel : public OnesideLinearConstraintMixin, void _add_callbacks(const ExpressionGraph &graph, const Outputs &outputs); void _add_callback(const ExpressionGraph &graph, const std::vector &outputs, const std::vector &constraints); - void _register_callback(CallbackEvaluator *evaluator); + void _register_callback(Evaluator *evaluator); void _update(); void _pre_solve(); void _solve(); diff --git a/lib/knitro_model.cpp b/lib/knitro_model.cpp index 1ffc098..cde9672 100644 --- a/lib/knitro_model.cpp +++ b/lib/knitro_model.cpp @@ -836,11 +836,11 @@ double KNITROModel::get_obj_value() const return _get_value(knitro::KN_get_obj_value); } -void KNITROModel::_register_callback(CallbackEvaluator *evaluator) +void KNITROModel::_register_callback(Evaluator *evaluator) { auto f = [](KN_context *, CB_context *cb, KN_eval_request *req, KN_eval_result *res, void *data) -> int { - auto evaluator = static_cast *>(data); + auto evaluator = static_cast(data); if (evaluator->is_objective()) { evaluator->eval_fun(req->x, res->obj); @@ -854,7 +854,7 @@ void KNITROModel::_register_callback(CallbackEvaluator *evaluator) auto g = [](KN_context *, CB_context *cb, KN_eval_request *req, KN_eval_result *res, void *data) -> int { - auto evaluator = static_cast *>(data); + auto evaluator = static_cast(data); if (evaluator->is_objective()) { evaluator->eval_jac(req->x, res->objGrad); @@ -868,7 +868,7 @@ void KNITROModel::_register_callback(CallbackEvaluator *evaluator) auto h = [](KN_context *, CB_context *cb, KN_eval_request *req, KN_eval_result *res, void *data) -> int { - auto evaluator = static_cast *>(data); + auto evaluator = static_cast(data); if (evaluator->is_objective()) { evaluator->eval_hess(req->x, req->sigma, res->hess); @@ -900,7 +900,7 @@ void KNITROModel::_register_callback(CallbackEvaluator *evaluator) void KNITROModel::_add_callback(const ExpressionGraph &graph, const std::vector &outputs, const std::vector &constraints) { - auto evaluator_ptr = std::make_unique>(); + auto evaluator_ptr = std::make_unique(); auto *evaluator = evaluator_ptr.get(); evaluator->indexVars.resize(graph.n_variables()); for (size_t i = 0; i < graph.n_variables(); i++) diff --git a/tests/test_lukvle10.py b/tests/test_lukvle10.py new file mode 100644 index 0000000..3c128dd --- /dev/null +++ b/tests/test_lukvle10.py @@ -0,0 +1,54 @@ +import pytest + +import pyoptinterface as poi +from pyoptinterface import ipopt, nl, copt + + +def test_nlp_lukvle10(nlp_model_ctor): + # LUKSAN-VLCEK Problem 10 + # + # min sum[i=1..n] (x[2i]^2)^(x[2i+1]^2 + 1) + (x[2i+1]^2)^(x[2i]^2 + 1) + # s.t. (3 - 2*x[i+1])*x[i+1] + 1 - x[i] - 2*x[i+2] = 0, i = 1,...,2n-2 + # + # Starting point: x[2i] = 1, x[2i+1] = -1 + model = nlp_model_ctor() + if isinstance(model, ipopt.Model): + # LUKVLE10 is too large and IPOPT raises a bad_alloc error. + pytest.skip("lukvle10 is too large to be supported with IPOPT") + if isinstance(model, copt.Model): + # LUKVLE10 is too large the current license of COpt supports up + # to 2000 variables. + pytest.skip("lukvle10 is too large to be supported with COpt") + + n = 100000 + x = model.add_m_variables(2 * n, name="x") + + for i in range(2 * n - 2): + model.add_quadratic_constraint( + (3 - 2 * x[i + 1]) * x[i + 1] + 1 - x[i] - 2 * x[i + 2], + poi.ConstraintSense.Equal, + 0.0, + name=f"c{i}", + ) + + with nl.graph(): + for i in range(n): + x2i_sq = x[2 * i] * x[2 * i] + x2ip1_sq = x[2 * i + 1] * x[2 * i + 1] + model.add_nl_objective( + nl.pow(x2i_sq, x2ip1_sq + 1) + nl.pow(x2ip1_sq, x2i_sq + 1) + ) + + for i in range(n): + model.set_variable_attribute(x[2 * i], poi.VariableAttribute.PrimalStart, 1.0) + model.set_variable_attribute( + x[2 * i + 1], poi.VariableAttribute.PrimalStart, -1.0 + ) + + model.optimize() + + termination_status = model.get_model_attribute(poi.ModelAttribute.TerminationStatus) + assert ( + termination_status == poi.TerminationStatusCode.LOCALLY_SOLVED + or termination_status == poi.TerminationStatusCode.OPTIMAL + )