diff --git a/include/mp/flat/nl_expr/constr_nl.h b/include/mp/flat/nl_expr/constr_nl.h index 29ed21674..0d3938f16 100644 --- a/include/mp/flat/nl_expr/constr_nl.h +++ b/include/mp/flat/nl_expr/constr_nl.h @@ -122,7 +122,7 @@ class NLBaseAssign /// Imitate double coef(int i) const { assert(!i); return -1.0; } /// Imitate - double var(int i) const { assert(!i); return GetVar(); } + int var(int i) const { assert(!i); return GetVar(); } /// Throw - should not be used VarArray1 GetArguments() const { MP_RAISE("No marking for NL items"); } diff --git a/include/mp/nl-reader.h b/include/mp/nl-reader.h index a27e04904..c53c02a69 100644 --- a/include/mp/nl-reader.h +++ b/include/mp/nl-reader.h @@ -2131,7 +2131,8 @@ class NLProblemBuilder { /// Add variables void AddVariables(const NLHeader& h) { // Distinguish NL variable order - // See D.M.Gay, Hooking Your Solver to AMPL; and Writing .NL Files, + // See D.M.Gay, Hooking Your Solver to AMPL; + // D.M.Gay, Writing .NL Files; // and, e.g., // github.com/jump-dev/MathOptInterface.jl/blob/master/src/FileFormats/NL/README.md int k=0; // current block position diff --git a/solvers/mp2nl/mp2nlmodelapi.cc b/solvers/mp2nl/mp2nlmodelapi.cc index e11550fc1..55e8030c5 100644 --- a/solvers/mp2nl/mp2nlmodelapi.cc +++ b/solvers/mp2nl/mp2nlmodelapi.cc @@ -27,6 +27,7 @@ void MP2NLModelAPI::InitProblemModificationPhase(const FlatModelInfo* flat_model flat_model_info->GetNumberOfConstraintsOfGroup(CG_Logical)); sos_info_.reserve( flat_model_info->GetNumberOfConstraintsOfGroup(CG_SOS)); + hdr_is_current_ = false; } void MP2NLModelAPI::AddVariables(const VarArrayDef& vad) { @@ -36,6 +37,8 @@ void MP2NLModelAPI::AddVariables(const VarArrayDef& vad) { if (vad.pnames()) var_names_ = {vad.pnames(), (size_t)vad.size()}; mark_data_.col_sizes_orig_.resize(var_lbs_.size()); + is_var_nlo_.resize(vad.size()); + is_var_nlc_.resize(vad.size()); } @@ -131,19 +134,23 @@ void MP2NLModelAPI::AddConstraint(const SOS2Constraint& sos) template MP2NL_Expr MP2NLModelAPI::AddExpression( const Expr &expr, ExpressionTypeID eid) { + SparsityTmp spstmp; // Only should visit such arguments // which are true expressions, i.e., - // not the pure-variable parts. - // Thus, need to specialize for NLDefVar. - VisitArguments(expr, [this](MP2NL_Expr mp2nle){ + // not the pure-variable parts? + // Thus, need to specialize for NLDefVar? + // But currently seems to be fine as is, + // because RegisterExpression() only takes true expressions. + VisitArguments(expr, [this,&spstmp](MP2NL_Expr mp2nle){ RegisterExpression(mp2nle); + MergeSparsityTmp(spstmp, mp2nle); }); - return StoreMP2NLExprID(expr, eid); + return StoreMP2NLExprID(expr, eid, spstmp); } template MP2NL_Expr MP2NLModelAPI::StoreMP2NLExprID( - const Expr &expr, ExpressionTypeID eid) { + const Expr &expr, ExpressionTypeID eid, const SparsityTmp& spars) { // std::printf(" Storing MP2NL_Expr[%ld]: type %s, logical = %d, @%p\n", // expr_info_.size(), // expr.GetFlatConstraint().GetTypeName(), @@ -151,6 +158,7 @@ MP2NL_Expr MP2NLModelAPI::StoreMP2NLExprID( expr_info_.push_back( MakeItemInfo(expr, eid, IsLogical(expr))); expr_counter_.push_back(0); + expr_sparsity_.push_back( {spars.begin(), spars.end()} ); return MakeExprID( int(expr_info_.size()-1) ); } @@ -158,14 +166,24 @@ void MP2NLModelAPI::RegisterExpression(MP2NL_Expr expr) { CountExpression(expr); } +void MP2NLModelAPI::MergeSparsityTmp( + MP2NLModelAPI::SparsityTmp& spstmp, MP2NL_Expr expr) { + if (expr.IsVariable()) { + spstmp.insert(expr.GetVarIndex()); + } else if (expr.IsExpression()) { + const auto& sps4expr = expr_sparsity_.at(expr.GetExprIndex()); + spstmp.insert(sps4expr.begin(), sps4expr.end()); + } // That's it +} + /// Count expression depending on its kind. /// This duplicates the counters in FlatConverter /// -- could check equality. void MP2NLModelAPI::CountExpression(MP2NL_Expr expr) { if (expr.IsExpression()) { auto index = expr.GetExprIndex(); - assert(index >=0 && index < expr_counter_.size() - && index < expr_info_.size()); + assert(index >=0 && index < (int)expr_counter_.size() + && index < (int)expr_info_.size()); ++expr_counter_[index]; } // @todo here also variable usage dep. on top-level item? } @@ -243,36 +261,223 @@ void MP2NLModelAPI::PrepareModel() { MapExpressions(); MarkVars(); SortVars(); - MarkItems(); } void MP2NLModelAPI::MapExpressions() { - for (const auto& info: obj_info_) { - MapExprTreeFromItemInfo(info); + ResetObjMetaInfo(); + is_alg_con_nl_.resize(alg_con_info_.size()); // grow flags array + for (int i=0; i<(int)obj_info_.size(); ++i) { + MapExprTreeFromItemInfo(i, obj_info_[i], 2); } - for (const auto& info: alg_con_info_) { - MapExprTreeFromItemInfo(info); + for (int i=0; i<(int)alg_con_info_.size(); ++i) { + if (!alg_con_info_[i].IsExprTreeMarked()) { + MapExprTreeFromItemInfo(i, alg_con_info_[i], 0); + alg_con_info_[i].SetExprTreeMarked(); + } } - for (const auto& info: log_con_info_) { - MapExprTreeFromItemInfo(info); + for (int i=0; i<(int)log_con_info_.size(); ++i) { + if (!log_con_info_[i].IsExprTreeMarked()) { + MapExprTreeFromItemInfo(i, log_con_info_[i], 1); + log_con_info_[i].SetExprTreeMarked(); + } } + FinishMapExpressions(); ResetIniExprRetrievedFlags(); } +void MP2NLModelAPI::MapExprTreeFromItemInfo( + int i_item, const ItemInfo& info, int kind) { + auto mp2nlexpr // so that AddExpression() is called + = info.GetDispatcher().GetExpression(info.GetPItem()); + RegisterExpression(mp2nlexpr); // tree root + if (1 != kind) { // alg con / obj + MergeItemSparsity(info, kind, mp2nlexpr); + UpdateAlgebraicMetaInfo(info, kind); // nnz, colsizes, n ranges/eqns + } + MarkNLVars(i_item, mp2nlexpr, kind); // also for logical cons +} + +void MP2NLModelAPI::ResetObjMetaInfo() { + mark_data_.n_obj_nz_ = 0; + mark_data_.nnlo_ = 0; + is_var_nlo_.clear(); + is_var_nlo_.resize(var_lbs_.size()); +} + +void MP2NLModelAPI::MergeItemSparsity( + const ItemInfo &info, int , MP2NL_Expr expr) { + if (expr.IsEmptyExpr()) { // Just use existing linear part + auto lp = GetLinPart(info); + info.SetExtLinPart(std::move(lp)); + } else { // proper expr or variable + std::unordered_map linp_ext; + auto lp = GetLinPart(info); + for (auto i=lp.size(); i--; ) + linp_ext[lp.vars()[i]] += lp.coefs()[i]; + if (expr.IsVariable()) { // @todo should have been inlined + linp_ext[expr.GetVarIndex()] += 1.0; // Merge 1 var + } else { + const auto& spars = expr_sparsity_.at(expr.GetExprIndex()); + for (auto v: spars) + linp_ext[v]; // Merge expression sparsity + } // https://stackoverflow.com/questions/33325470/ + // c-is-there-a-way-to-initialize-unordered-map-using- + // contents-of-vector-as-keys + std::vector vars; // @todo SmallVec + vars.reserve(linp_ext.size()); + std::vector coefs; + coefs.reserve(linp_ext.size()); + for (const auto& val: linp_ext) { + coefs.push_back(val.second); + vars.push_back(val.first); + } + info.SetExtLinPart( // Put joint result into \a info + {std::move(coefs), std::move(vars)} ); + } +} + +inline + MP2NLModelAPI::LinPartRefOrStorage + GetLPRoS(const LinTerms& lt) { + return {lt.coefs(), lt.vars()}; +} + +template +inline + MP2NLModelAPI::LinPartRefOrStorage + GetLPRoS(const NLBaseAssign& nla) { + std::vector c{1, nla.coef(0)}; // @todo: SmallVec + std::vector v{1, nla.var(0)}; + return {std::move(c), std::move(v)}; // std::move to pass storage +} + +MP2NLModelAPI::LinPartRefOrStorage +MP2NLModelAPI::GetLinPart(const ItemInfo &info) { + const auto* pitem = info.GetPItem(); + switch (info.GetStaticTypeID()) { + case StaticItemTypeID::ID_LinearObjective: + case StaticItemTypeID::ID_NLObjective: { + return GetLPRoS( ((LinearObjective*)(pitem))->GetLinTerms() ); + } + case StaticItemTypeID::ID_LinConRange: { + return GetLPRoS( ((LinConRange*)(pitem))->GetLinTerms() ); + } + case StaticItemTypeID::ID_LinConLE: { + return GetLPRoS( ((LinConLE*)(pitem))->GetLinTerms() ); + } + case StaticItemTypeID::ID_LinConEQ: { + return GetLPRoS( ((LinConEQ*)(pitem))->GetLinTerms() ); + } + case StaticItemTypeID::ID_LinConGE: { + return GetLPRoS( ((LinConGE*)(pitem))->GetLinTerms() ); + } + case StaticItemTypeID::ID_NLConstraint: { + return GetLPRoS( + ((const NLConstraint*)(pitem))->GetMainCon().GetLinTerms() ); + } + case StaticItemTypeID::ID_NLAssignLE: { + return GetLPRoS( *((NLAssignLE*)(pitem)) ); + } + case StaticItemTypeID::ID_NLAssignEQ: { + return GetLPRoS( *((NLAssignEQ*)(pitem)) ); + } + case StaticItemTypeID::ID_NLAssignGE: { + return GetLPRoS( *((NLAssignGE*)(pitem)) ); + } + case StaticItemTypeID::ID_NLComplementarity: { + return GetLPRoS( // Do we need the \a compl_var? + ((const NLComplementarity*)(pitem))->GetExpression().GetBody() ); + } + default: + MP_RAISE("Unknown objective or algebraic constraint type"); + } +} + +void MP2NLModelAPI::UpdateAlgebraicMetaInfo( + const ItemInfo &info, int kind) { + const auto& lp_ext = info.GetExtLinPart(); + if (kind) { + mark_data_.n_obj_nz_ += lp_ext.size(); + } else { + mark_data_.n_con_nz_ += lp_ext.size(); + Add2ColSizes(lp_ext.vars()); + info.GetDispatcher().MarkRangeOrEqn(info.GetPItem()); + } +} + +void MP2NLModelAPI::MarkNLVars( + int i_item, MP2NL_Expr expr, int kind) { + if (expr.IsEmptyExpr()) { + // nothing + } else { + std::vector& var_nl_flags + = kind ? is_var_nlo_ : is_var_nlc_; + if (expr.IsVariable()) { // we still count it as non-linear + var_nl_flags.at(expr.GetVarIndex()) = true; + } else { + for (auto v: expr_sparsity_.at(expr.GetExprIndex())) { + assert(v>=0 && v<(int)var_lbs_.size()); + var_nl_flags[v] = true; + } + } + if (0==kind) { + is_alg_con_nl_.at(i_item) = true; // nonlinear cons + ++mark_data_.nnlc_; + } else { + ++mark_data_.nnlo_; + } + } +} + +void MP2NLModelAPI::FinishMapExpressions() { } + void MP2NLModelAPI::MarkVars() { + mark_data_.var_prior_.clear(); mark_data_.var_prior_.resize(var_lbs_.size()); - mark_data_.n_var_lin_bin_ = 0; - mark_data_.n_var_lin_int_ = 0; + mark_data_.n_bin_var_lin_ = 0; + mark_data_.n_int_var_lin_ = 0; + mark_data_.n_cont_var_nl_con_ = 0; + mark_data_.n_cont_var_nl_obj_ = 0; + mark_data_.n_cont_var_nl_both_ = 0; + mark_data_.n_int_var_nl_con_ = 0; + mark_data_.n_int_var_nl_obj_ = 0; + mark_data_.n_int_var_nl_both_ = 0; for (auto i=var_lbs_.size(); i--; ) { mark_data_.var_prior_[i].second = i; + bool f_int_non_bin{}, f_bin{}, f_int_or_bin{}; if (var::Type::INTEGER == var_types_[i]) { - if (!var_lbs_[i] && 1.0==var_ubs_[i]) { + f_int_or_bin = true; + if (!var_lbs_[i] && 1.0==var_ubs_[i]) { // binary var mark_data_.var_prior_[i].first = 1; - ++mark_data_.n_var_lin_bin_; + f_bin = true; } else { mark_data_.var_prior_[i].first = 2; - ++mark_data_.n_var_lin_int_; + f_int_non_bin = true; + } + } + if (is_var_nlc_[i]) { + if (is_var_nlo_[i]) { + mark_data_.var_prior_[i].first -= 9; + f_int_or_bin + ? ++mark_data_.n_int_var_nl_both_ + : ++mark_data_.n_cont_var_nl_both_; + } else { + mark_data_.var_prior_[i].first -= 6; + f_int_or_bin + ? ++mark_data_.n_int_var_nl_con_ + : ++mark_data_.n_cont_var_nl_con_; } + } else if (is_var_nlo_[i]) { + mark_data_.var_prior_[i].first -= 3; + f_int_or_bin + ? ++mark_data_.n_int_var_nl_obj_ + : ++mark_data_.n_cont_var_nl_obj_; + } else { // 'linear' var + if (f_bin) + ++mark_data_.n_bin_var_lin_; + else + if (f_int_non_bin) + ++mark_data_.n_int_var_lin_; } } } @@ -287,13 +492,7 @@ void MP2NLModelAPI::SortVars() { } } -void MP2NLModelAPI::MarkItems() { - ItemMarkingData prm; - for (auto& info: alg_con_info_) - info.GetDispatcher().MarkItem(info.GetPItem(), prm); -} - -void MP2NLModelAPI::Add2ColSizes(const std::vector& vars) { +void MP2NLModelAPI::Add2ColSizes(ArrayRef vars) { for (auto v: vars) ++mark_data_.col_sizes_orig_[v]; } @@ -310,8 +509,8 @@ NLHeader MP2NLModelAPI::DoMakeHeader() { hdr.num_logical_cons = log_con_info_.size(); /** Total number of nonlinear constraints. */ - hdr.num_nl_cons = 0; - hdr.num_nl_objs = 0; + hdr.num_nl_cons = mark_data_.nnlc_; + hdr.num_nl_objs = mark_data_.nnlo_; hdr.num_compl_conds = 0; hdr.num_nl_compl_conds = 0; hdr.num_compl_dbl_ineqs = 0; @@ -321,20 +520,25 @@ NLHeader MP2NLModelAPI::DoMakeHeader() { hdr.num_nl_net_cons = 0; hdr.num_linear_net_cons = 0; + /** Number of nonlinear variables in both constraints and objectives. */ + hdr.num_nl_vars_in_both + = mark_data_.n_cont_var_nl_both_ + mark_data_.n_int_var_nl_both_; + /** Number of nonlinear variables in constraints including nonlinear variables in both constraints and objectives. */ - hdr.num_nl_vars_in_cons = 0; + hdr.num_nl_vars_in_cons + = hdr.num_nl_vars_in_both + + mark_data_.n_cont_var_nl_con_ + mark_data_.n_int_var_nl_con_; /** Number of nonlinear variables in objectives including nonlinear variables in both constraints and objectives. */ - hdr.num_nl_vars_in_objs = 0; - - /** Number of nonlinear variables in both constraints and objectives. */ - hdr.num_nl_vars_in_both = 0; + hdr.num_nl_vars_in_objs + = hdr.num_nl_vars_in_cons + + mark_data_.n_cont_var_nl_obj_ + mark_data_.n_int_var_nl_obj_; // Miscellaneous // ------------- @@ -349,21 +553,21 @@ NLHeader MP2NLModelAPI::DoMakeHeader() { // ------------------------------------ /** Number of linear binary variables. */ - hdr.num_linear_binary_vars = mark_data_.n_var_lin_bin_; + hdr.num_linear_binary_vars = mark_data_.n_bin_var_lin_; /** Number of linear non-binary integer variables. */ - hdr.num_linear_integer_vars = mark_data_.n_var_lin_int_; + hdr.num_linear_integer_vars = mark_data_.n_int_var_lin_; /** Number of integer nonlinear variables in both constraints and objectives. */ - hdr.num_nl_integer_vars_in_both = 0; + hdr.num_nl_integer_vars_in_both = mark_data_.n_int_var_nl_both_; /** Number of integer nonlinear variables just in constraints. */ - hdr.num_nl_integer_vars_in_cons = 0; + hdr.num_nl_integer_vars_in_cons = mark_data_.n_int_var_nl_con_; /** Number of integer nonlinear variables just in objectives. */ - hdr.num_nl_integer_vars_in_objs = 0; + hdr.num_nl_integer_vars_in_objs = mark_data_.n_int_var_nl_obj_; // Information about nonzeros // -------------------------- @@ -437,20 +641,7 @@ int MP2NLModelAPI::ObjType(int i) { template void MP2NLModelAPI::FeedObjGradient(int i, ObjGradWriterFactory& svwf) { - switch (obj_info_.at(i).GetStaticTypeID()) { - case StaticItemTypeID::ID_LinearObjective: - case StaticItemTypeID::ID_NLObjective: { - const auto& obj = *((LinearObjective*)(obj_info_[i].GetPItem())); - if (obj.num_terms()) { - auto svw = svwf.MakeVectorWriter(obj.num_terms()); - for (int j=0; j @@ -548,51 +739,19 @@ void MP2NLModelAPI::FeedConBounds(ConBoundsWriter& cbw) { template void MP2NLModelAPI::FeedLinearConExpr(int i, ConLinearExprWriterFactory& svwf) { - const auto* pitem = alg_con_info_[i].GetPItem(); - switch (alg_con_info_[i].GetStaticTypeID()) { - case StaticItemTypeID::ID_LinConRange: { - FeedLinearConBody( *((LinConRange*)(pitem)), svwf); - } break; - case StaticItemTypeID::ID_LinConLE: { - FeedLinearConBody( *((LinConLE*)(pitem)), svwf); - } break; - case StaticItemTypeID::ID_LinConEQ: { - FeedLinearConBody( *((LinConEQ*)(pitem)), svwf); - } break; - case StaticItemTypeID::ID_LinConGE: { - FeedLinearConBody( *((LinConGE*)(pitem)), svwf); - } break; - case StaticItemTypeID::ID_NLConstraint: { - FeedLinearConBody( - ((const NLConstraint*)(pitem))->GetMainCon(), svwf); - } break; - case StaticItemTypeID::ID_NLAssignLE: { - FeedLinearConBody( *((NLAssignLE*)(pitem)), svwf); - } break; - case StaticItemTypeID::ID_NLAssignEQ: { - FeedLinearConBody( *((NLAssignEQ*)(pitem)), svwf); - } break; - case StaticItemTypeID::ID_NLAssignGE: { - FeedLinearConBody( *((NLAssignGE*)(pitem)), svwf); - } break; - case StaticItemTypeID::ID_NLComplementarity: { - FeedLinearConBody( - ((const NLComplementarity*)(pitem))->GetExpression(), svwf); - } break; - default: - MP_RAISE("Unknown algebraic constraint type"); - } + FeedExtLinPart(alg_con_info_[i], svwf); } -template -void MP2NLModelAPI::FeedLinearConBody( - const ExprBody& algcon, +template +void MP2NLModelAPI::FeedExtLinPart( + const ItemInfo& item, ConLinearExprWriterFactory& svwf) { - if (algcon.size()) { - auto svw = svwf.MakeVectorWriter(algcon.size()); - for (int j=0; j + #include "mp/env.h" #include "mp2nlcommon.h" #include "mp/flat/nl_expr/model_api_base.h" @@ -17,6 +19,9 @@ namespace mp { /// 3. normal expression node. /// /// @note To create MP2NL_Expr, use MakeMP2NL_... below. +/// +/// @note We have no 'constant non-zero' expression +/// because FlatConverter should not produce them class MP2NL_Expr { public: /// Construct @@ -663,11 +668,6 @@ class MP2NLModelAPI template void FeedLinearConExpr(int i, ConLinearExprWriterFactory& svwf); - template - void FeedLinearConBody( - const ExprBody& algcon, - ConLinearExprWriterFactory& svwf); - /** Feed nonlinear expression of constraint \a i. * Algebraic constraints (num_algebraic_cons) * come before logical (num_logical_cons). @@ -947,6 +947,11 @@ class MP2NLModelAPI void PrepareModel(); /// Map expressions. + /// Fill sparsity patterns and merge them with the linear parts, + /// fill nonzero counters, column sizes, and variable/constraint + /// nonlinearity flags. + /// @note Should do this incrementally when called + /// from the MO Emulator. void MapExpressions(); /// Mark variables. @@ -956,16 +961,22 @@ class MP2NLModelAPI /// Sort variables void SortVars(); - /// Mark and count constraint and variable kinds - /// (non/linear/integer, range/eqns.) - void MarkItems(); - NLHeader DoMakeHeader(); /// Parameters passed when marking variables in an expression tree struct ItemMarkingData { - int n_var_lin_bin_ {0}; - int n_var_lin_int_ {0}; + int n_bin_var_lin_ {0}; + int n_int_var_lin_ {0}; + + /// NL continuous vars in constraints/obj/both + int n_cont_var_nl_con_ {0}; + int n_cont_var_nl_obj_ {0}; + int n_cont_var_nl_both_ {0}; + + /// NL integer vars in constraints/objs/both + int n_int_var_nl_con_ {0}; + int n_int_var_nl_obj_ {0}; + int n_int_var_nl_both_ {0}; int n_ranges_ {0}; int n_eqns_ {0}; @@ -973,6 +984,9 @@ class MP2NLModelAPI std::size_t n_obj_nz_ {0}; std::size_t n_con_nz_ {0}; + int nnlo_ {0}; + int nnlc_ {0}; + std::vector< std::pair< int, int > > var_prior_; // new index -> var weight, orig. index std::vector var_order_12_; // new index -> old index std::vector var_order_21_; // old index -> new index @@ -987,76 +1001,52 @@ class MP2NLModelAPI /// Placeholder for item marker template - void MarkItem(const Item& , ItemMarkingData& ) { + void MarkRangeOrEqn(const Item& ) { MP_UNSUPPORTED(std::string("MP2NLModelAPI::MarkItem() not supported for ") + typeid(Item).name()); } /// Mark linear part of any objective - void MarkItem(const LinearObjective& lo, ItemMarkingData& prm) { - mark_data_.n_obj_nz_ += lo.vars().size(); // count nnz - } + void MarkRangeOrEqn(const LinearObjective& ) { } /// Mark LinConRange - void MarkItem(const LinConRange& lcr, ItemMarkingData& prm) { - mark_data_.n_con_nz_ += lcr.size(); + void MarkRangeOrEqn(const LinConRange& lcr) { Add2ColSizes(lcr.vars()); if (lcr.lb() > MinusInfinity() - && lcr.ub() < Infinity() - && lcr.lb() < lcr.ub()) - ++prm.n_ranges_; + && lcr.ub() < Infinity()) + if (lcr.lb() < lcr.ub()) + ++mark_data_.n_ranges_; + else + ++mark_data_.n_eqns_; } /// Mark LinConEQ - void MarkItem(const LinConLE& lcr, ItemMarkingData& prm) { - mark_data_.n_con_nz_ += lcr.size(); - Add2ColSizes(lcr.vars()); - } + void MarkRangeOrEqn(const LinConLE& ) { } /// Mark LinConEQ - void MarkItem(const LinConEQ& lcr, ItemMarkingData& prm) { - mark_data_.n_con_nz_ += lcr.size(); - Add2ColSizes(lcr.vars()); - ++prm.n_eqns_; + void MarkRangeOrEqn(const LinConEQ& ) { + ++mark_data_.n_eqns_; } /// Mark LinConEQ - void MarkItem(const LinConGE& lcr, ItemMarkingData& prm) { - mark_data_.n_con_nz_ += lcr.size(); - Add2ColSizes(lcr.vars()); - } + void MarkRangeOrEqn(const LinConGE& ) { } - /// Mark NLconstraint - void MarkItem(const NLConstraint& nlc, ItemMarkingData& prm) { - MarkItem(nlc.GetMainCon(), prm); + /// Mark NLConstraint + void MarkRangeOrEqn(const NLConstraint& nlc) { + MarkRangeOrEqn(nlc.GetMainCon()); } /// Mark NLAssign - void MarkItem(const NLAssignLE& nlc, ItemMarkingData& prm) { - mark_data_.n_con_nz_ += 1; - Add2ColSizes({ nlc.GetVar() }); - } + void MarkRangeOrEqn(const NLAssignLE& ) { } /// Mark NLAssign - void MarkItem(const NLAssignEQ& nlc, ItemMarkingData& prm) { - mark_data_.n_con_nz_ += 1; - Add2ColSizes({ nlc.GetVar() }); - ++prm.n_eqns_; + void MarkRangeOrEqn(const NLAssignEQ& nlc) { + ++mark_data_.n_eqns_; } /// Mark NLAssign - void MarkItem(const NLAssignGE& nlc, ItemMarkingData& prm) { - mark_data_.n_con_nz_ += 1; - Add2ColSizes({ nlc.GetVar() }); - } + void MarkRangeOrEqn(const NLAssignGE& lc) { } /// Add to col sizes - void Add2ColSizes(const std::vector& vars); - - /// Placeholder for bounds writer - template - void WriteBounds(const Item& , const NLWParams&) { - MP_UNSUPPORTED(std::string("MP2NLModelAPI::WriteBounds() not supported for ") - + typeid(Item).name()); - } + void Add2ColSizes(ArrayRef vars); /// We need item type ID for manual dispatch. @@ -1141,9 +1131,9 @@ class MP2NLModelAPI /// ItemName virtual const char* GetName(void* pitem) = 0; - virtual void MarkItem(void* pitem, ItemMarkingData& vmp) = 0; - virtual void WriteBounds(void* pitem, const NLWParams& nlwp) = 0; - virtual void MarkExprTree(void* pitem) = 0; + virtual void MarkRangeOrEqn(void* pitem) = 0; + /// Get expression + virtual MP2NL_Expr GetExpression(void* pitem) = 0; /// Placeholder for the item category getter: /// static item vs expression @@ -1173,16 +1163,12 @@ class MP2NLModelAPI { return GetMAPI().GetItemName(*(const Item*)pitem); } /// Var marking - void MarkItem(void* pitem, ItemMarkingData& vmp) override - { GetMAPI().MarkItem(*(const Item*)pitem, vmp); } - - /// Write bounds - void WriteBounds(void* pitem, const NLWParams& nlwp) override - { GetMAPI().WriteBounds(*(const Item*)pitem, nlwp); } + void MarkRangeOrEqn(void* pitem) override + { GetMAPI().MarkRangeOrEqn(*(const Item*)pitem); } - /// Add + mark expressions - void MarkExprTree(void* pitem) override - { GetMAPI().MapExprTreeFromItem(*(const Item*)pitem); } + /// Get expression + MP2NL_Expr GetExpression(void* pitem) override + { return GetMAPI().GetExpression(*(const Item*)pitem); } /// Item category getter: /// static item vs expression @@ -1300,6 +1286,37 @@ class MP2NLModelAPI CREATE_EXPRESSION_DISPATCHER(Div) + +public: + /// Reference/storage for (extended) linear part + /// containing the nonlinear expression sparsity pattern + /// @todo SmallVec's for storage in ArrayRef + /// - possibly as a template parameter + class LinPartRefOrStorage { + public: + /// Construct default + LinPartRefOrStorage() = default; + /// Construct + /// @note careful when passing with stored data - use std::move + /// @todo SmallVec's as data + LinPartRefOrStorage(ArrayRef c, ArrayRef v) + : coefs_(std::move(c)), vars_(std::move(v)) + { assert(Check()); } + /// Check + bool Check() const { return coefs_.size()==vars_.size(); } + /// Size + int size() const { assert(Check()); return coefs_.size(); } + /// Coefs + ArrayRef coefs() const { return coefs_; } + /// Vars + ArrayRef vars() const { return vars_; } + private: + ArrayRef coefs_; + ArrayRef vars_; + }; + + +protected: /// Constraint/objective/expression info. /// We rely on the pointers staying valid. class ItemInfo { @@ -1328,6 +1345,11 @@ class MP2NLModelAPI void* GetPItem() const { return p_item_; } /// Is logical? bool IsLogical() const { return f_logical_; } + /// Has the expression tree been marked? + /// This is for incremental model modification. + bool IsExprTreeMarked() const { return f_marked_; } + /// Set expr tree marked status + void SetExprTreeMarked() const { f_marked_=true; } /// Is a static type? bool IsItemTypeStatic() const { return GetDispatcher().IsItemTypeStatic(); } @@ -1338,12 +1360,21 @@ class MP2NLModelAPI ExpressionTypeID GetExprTypeID() const { return GetDispatcher().GetExpressionTypeID(); } + /// Get extended lin part + const LinPartRefOrStorage& GetExtLinPart() const + { return lin_part_ext_; } + /// Set extended linear part + void SetExtLinPart(LinPartRefOrStorage lp) const + { lin_part_ext_=std::move(lp); } + private: BasicItemDispatcher& disp_; void* p_item_; - bool f_logical_ {}; // StaticItemTypeID itemID_ {StaticItemTypeID::ID_None}; // ExpressionTypeID exprID_ {ExpressionTypeID::ID_None}; + bool f_logical_ {}; + mutable bool f_marked_ {}; // @todo just see if we have the extended linear part? + mutable LinPartRefOrStorage lin_part_ext_ {}; }; /// Fill static item info @@ -1370,20 +1401,33 @@ class MP2NLModelAPI void DispatchStaticItem(const ItemInfo& info, Lambda lambda); /// Map expressions from a single item info + /// @param kind: 0 - alg con, 1 - log con, 2 - obj /// @note needs to be in .h to be instantiated/inlined, /// at least for Clang 16. - void MapExprTreeFromItemInfo(const ItemInfo& info) - { info.GetDispatcher().MarkExprTree(info.GetPItem()); } + void MapExprTreeFromItemInfo( + int i_item, const ItemInfo& info, int kind); - /// Map expressions from a single item (type-safe) - /// @note needs to be in .h to be instantiated/inlined, - /// at least for Clang 16. - template - void MapExprTreeFromItem(const Item& item) { - auto mp2nlexpr - = GetExpression(item); // so that AddExpression() is called - RegisterExpression(mp2nlexpr); // mark the tree root - } + /// Reset objective meta info for resolve + /// (NNZ, nlo(i), nlb(i)) + void ResetObjMetaInfo(); + + /// Merge algcon/obj sparsity pattern with jacobian/gradient + void MergeItemSparsity( + const ItemInfo& info, int item_kind, MP2NL_Expr expr); + + /// Retrieve original linear part of an alg con / obj + LinPartRefOrStorage GetLinPart(const ItemInfo& info); + + /// NNZ, col sizes, N ranges and eqns + void UpdateAlgebraicMetaInfo(const ItemInfo& info, int kind); + + /// Register NL variables, + /// in particular a constraint's nonlinearity + void MarkNLVars(int i_item, MP2NL_Expr expr, int kind); + + /// Compute nlvc(i), nlvo(i), nlvb(i) etc + /// from mapping data + void FinishMapExpressions(); /// Reuse GetExpression() using BaseModelAPI::GetExpression; @@ -1413,11 +1457,20 @@ class MP2NLModelAPI MP2NL_Expr AddExpression( const Expr& expr, ExpressionTypeID eid); + /// Typedef intermediate sparsity pattern + using SparsityTmp = std::unordered_set; + + /// Merge sparsity pattern from an expression + /// into the 1st param + void MergeSparsityTmp(SparsityTmp& , MP2NL_Expr ); + + /// Typedef final sparsity pattern for an expression + using Sparsity4Expr = SmallVec; /// Make, store and return MP2NL_Expr for a given expession template MP2NL_Expr StoreMP2NLExprID( - const Expr& expr, ExpressionTypeID eid); + const Expr& expr, ExpressionTypeID eid, const SparsityTmp& ); /// Variable name const char* GetVarName(int v) const { @@ -1426,6 +1479,12 @@ class MP2NLModelAPI (int)var_names_.size() > v ? var_names_[v] : nullptr; } + /// Feed extended linear part + template + void FeedExtLinPart( + const ItemInfo& item, + ConLinearExprWriterFactory& svwf); + private: /// References to the model data. @@ -1443,6 +1502,12 @@ class MP2NLModelAPI std::vector expr_info_; std::vector expr_counter_; // usage counter + std::vector expr_sparsity_; + + /// if a var is nonlinear in obj/con + std::vector is_var_nlo_, is_var_nlc_; + /// If an alg con is nonlinear + std::vector is_alg_con_nl_; ItemMarkingData mark_data_; bool hdr_is_current_ {};