diff --git a/include/mp/flat/constr_general.h b/include/mp/flat/constr_general.h index 824849ce9..db3ac0ce3 100644 --- a/include/mp/flat/constr_general.h +++ b/include/mp/flat/constr_general.h @@ -191,15 +191,15 @@ using ComplementarityQuadratic = ComplementarityConstraint; /// Quadratic cone DEF_STATIC_CONSTR_WITH_PRM( QuadraticConeConstraint, VarArray, DblParamArray, "Quadratic cone p1*x1 >= sqrt((p2*x2)^2 + ...))," - " with factors p1..pn"); + " with factors p1..pn"); /// Rotated quadratic cone DEF_STATIC_CONSTR_WITH_PRM( RotatedQuadraticConeConstraint, VarArray, DblParamArray, "Rotated quadratic cone p1*x1*p2*x2 >= sqrt((p3*x3)^2 + ...))," - " x1, x2 >= 0, with factors p1..pn"); + " x1, x2 >= 0, with factors p1..pn"); /// Exponential cone DEF_STATIC_CONSTR_WITH_PRM( ExponentialConeConstraint, VarArray3, DblParamArray3, "Exponential cone p1*x1 >= p2*x2*exp(p3*x3 / (p2*x2))," - " x1, x2 >= 0, with factors p1..p3"); + " x1, x2 >= 0, with factors p1..p3"); /// Power cone DEF_STATIC_CONSTR_WITH_PRM( PowerConeConstraint, VarArray, DblParamArray, "Power cone with factors "); diff --git a/include/mp/flat/redef/conic/cones.h b/include/mp/flat/redef/conic/cones.h index 3ac467718..1dc6ce7c4 100644 --- a/include/mp/flat/redef/conic/cones.h +++ b/include/mp/flat/redef/conic/cones.h @@ -43,19 +43,19 @@ class Convert1ExpC; template class MCKeeper { public: - /// Typedef MCType - using MCType = ModelConverter; + /// Typedef MCType + using MCType = ModelConverter; - /// Constructor - MCKeeper(MCType& mc) : mc_(mc) { } + /// Constructor + MCKeeper(MCType& mc) : mc_(mc) { } protected: - /// Retrieve the MC - const MCType& MC() const { return mc_; } - MCType& MC() { return mc_; } + /// Retrieve the MC + const MCType& MC() const { return mc_; } + MCType& MC() { return mc_; } private: - MCType& mc_; + MCType& mc_; }; @@ -67,69 +67,69 @@ class MCKeeper { template class ConicConverter : public MCKeeper { public: - virtual ~ConicConverter() { } - - /// Constructor - ConicConverter(MCType& mc) : - MCKeeper(mc) { } - - /// Run conversions - void Run() { - // We _could_ walk everything just once - // and see which cones are there. - // Or, even walk expression trees from exponents, etc. - RunQCones(); - RunExpCones(); - } + virtual ~ConicConverter() { } + + /// Constructor + ConicConverter(MCType& mc) : + MCKeeper(mc) { } + + /// Run conversions + void Run() { + // We _could_ walk everything just once + // and see which cones are there. + // Or, even walk expression trees from exponents, etc. + RunQCones(); + RunExpCones(); + } protected: - void RunQCones() { - if (MC().IfPassSOCPCones()) { // convert everything to QuadraticCones. - Walk(); - Walk(); - Walk(); - - if (MC().GetNumberOfAddable((PowConstraint*)0)>0 || - MC().GetNumberOfAddable((AbsConstraint*)0)>0) { - Walk(); - Walk(); - Walk(); - } - } else - if (MC().IfPassQuadCon() && - (MC().GetNumberOfAddable((PowConstraint*)0)>0 || - MC().GetNumberOfAddable((AbsConstraint*)0)>0)) { - // Still collect QCones expressed by 2-norms. - // They are to be converted to quadratics. - Walk(); - Walk(); - Walk(); - } - } + void RunQCones() { + if (MC().IfPassSOCPCones()) { // convert everything to QuadraticCones. + Walk(); + Walk(); + Walk(); + + if (MC().GetNumberOfAddable((PowConstraint*)0)>0 || + MC().GetNumberOfAddable((AbsConstraint*)0)>0) { + Walk(); + Walk(); + Walk(); + } + } else + if (MC().IfPassQuadCon() && + (MC().GetNumberOfAddable((PowConstraint*)0)>0 || + MC().GetNumberOfAddable((AbsConstraint*)0)>0)) { + // Still collect QCones expressed by 2-norms. + // They are to be converted to quadratics. + Walk(); + Walk(); + Walk(); + } + } - void RunExpCones() { - if (MC().IfPassExpCones() && // convert everything to ExpCones. - MC().GetNumberOfAddable((ExpConstraint*)0)>0) { - Walk(); - Walk(); // also ExpA ?? - Walk(); + void RunExpCones() { + if (MC().IfPassExpCones() && // convert everything to ExpCones. + MC().GetNumberOfAddable((ExpConstraint*)0)>0) { + Walk(); + Walk(); // also ExpA ?? + Walk(); - Walk(); - Walk(); - Walk(); - } + Walk(); + Walk(); + Walk(); } + } - /// Walk a single constraint type + /// Walk a single constraint type template class CvtBase> - void Walk() { - auto& ck = MC().GetConstraintKeeper((Con*)nullptr); - ck.ForEachActive(Convert1 { MC() }); - } + void Walk() { + auto& ck = MC().GetConstraintKeeper((Con*)nullptr); + ck.ForEachActive(Convert1 { MC() }); + } - /// Retrieve the MC - using MCKeeper::MC; + /// Retrieve the MC + using MCKeeper::MC; }; @@ -140,8 +140,8 @@ class ConicConverter : public MCKeeper { template class Convert1QC : public MCKeeper { public: - /// Constructor - Convert1QC(MCType& mc) : MCKeeper(mc) { } + /// Constructor + Convert1QC(MCType& mc) : MCKeeper(mc) { } protected: @@ -155,25 +155,25 @@ class Convert1QC : public MCKeeper { int iSamePos=-1, iSameNeg=-1; }; - /// DoRun. Body: quadratic. - /// - /// Accept non-(+-1) coefficients. - /// - /// @param body: quadratic constraint body. - /// @param sens: -1 for <=, 1 for >=. - /// @param rhs: constraint's rhs. - /// - /// @return true iff converted. - bool DoRun(const QuadAndLinTerms& body, - int sens, double rhs) { - assert((sens==1 || sens==-1) && "sens 1 or -1 only"); + /// DoRun. Body: quadratic. + /// + /// Accept non-(+-1) coefficients. + /// + /// @param body: quadratic constraint body. + /// @param sens: -1 for <=, 1 for >=. + /// @param rhs: constraint's rhs. + /// + /// @return true iff converted. + bool DoRun(const QuadAndLinTerms& body, + int sens, double rhs) { + assert((sens==1 || sens==-1) && "sens 1 or -1 only"); if (body.GetLinTerms().size()<=1) { QPTermsTraits qptt; if (Characterize(body.GetQPTerms(), qptt)) return ProceedQCWithTraits(body, sens, rhs, qptt); - } - return false; - } + } + return false; + } /// Characterize QP terms /// @return false iff not suitable @@ -270,161 +270,161 @@ class Convert1QC : public MCKeeper { /// Rotated SOC. /// 1. sqrt(k*x1*x2) >= ||z [ | const ]||. /// 2. sqrt(k * x2) >= ... (one var is const in lhs) - /// + /// /// Considering const>=0. - /// Accept non-(+-1) coefficients. - /// - /// @param body: linear constraint body. - /// @param sens: -1 for <=, 1 for >=. - /// @param rhs: constraint's rhs. - /// - /// @return true iff converted. - bool DoRun(const LinTerms& lint, - int sens, double rhs) { - assert((sens==1 || sens==-1) && "sens 1 or -1 only"); - if (2 == lint.size() && // x>=y - 0.0 == std::fabs(rhs) && // rhs=0 - 0.0 > lint.coef(0)*lint.coef(1)) { // different signs - // For sens>0 we have x>=y and iX=(index of pos coef) - int iX = (lint.coef(1)>0.0); - if (sens<0) - iX = 1-iX; - int iY = 1-iX; - // Check if var(iY) := || vector-or-var || - if (auto rhs_args = CheckNorm2(lint.var(iY))) { - if (auto lhs_args = CheckSqrtXnXmNonneg(lint.var(iX))) - return ContinueRotatedSOC(lint, iX, iY, - lhs_args, rhs_args); + /// Accept non-(+-1) coefficients. + /// + /// @param body: linear constraint body. + /// @param sens: -1 for <=, 1 for >=. + /// @param rhs: constraint's rhs. + /// + /// @return true iff converted. + bool DoRun(const LinTerms& lint, + int sens, double rhs) { + assert((sens==1 || sens==-1) && "sens 1 or -1 only"); + if (2 == lint.size() && // x>=y + 0.0 == std::fabs(rhs) && // rhs=0 + 0.0 > lint.coef(0)*lint.coef(1)) { // different signs + // For sens>0 we have x>=y and iX=(index of pos coef) + int iX = (lint.coef(1)>0.0); + if (sens<0) + iX = 1-iX; + int iY = 1-iX; + // Check if var(iY) := || vector-or-var || + if (auto rhs_args = CheckNorm2(lint.var(iY))) { + if (auto lhs_args = CheckSqrtXnXmNonneg(lint.var(iX))) + return ContinueRotatedSOC(lint, iX, iY, + lhs_args, rhs_args); return ContinueStdSOC(lint.coef(iX), lint.var(iX), lint.coef(iY), rhs_args); - } + } } else if (1 == lint.size() && // const>=y. 0.0 >= rhs*sens && // either ... <= rhs - // or ... >= rhs(<0) + // or ... >= rhs(<0) 0.0 >= lint.coef(0)*sens) { // In the other case, - // it's -k*sqrt(...) <= rhs(>0), k>0, - // which is always true TODO + // it's -k*sqrt(...) <= rhs(>0), k>0, + // which is always true TODO if (auto rhs_args = CheckNorm2(lint.var(0))) { return ContinueStdSOC(std::fabs(rhs), MC().MakeFixedVar( 1.0 ), lint.coef(0), rhs_args); } } - return false; - } - - /// Typedef for subexpression checkup result, - /// whether it represents some part of an SOCP cone. - struct ConeArgs { - std::vector coefs_; - std::vector vars_; + return false; + } + + /// Typedef for subexpression checkup result, + /// whether it represents some part of an SOCP cone. + struct ConeArgs { + std::vector coefs_; + std::vector vars_; double const_term = 0.0; double coef_extra = 0.0; - /// Result vars of expressions to un-use. - std::vector res_vars_to_delete_; - /// operator bool - operator bool() const - { assert(check()); return !coefs_.empty(); } - /// size() - size_t size() const - { assert(check()); return coefs_.size(); } - /// check() - bool check() const { return coefs_.size()==vars_.size(); } - }; - - /// Check if the variable is defined by an expression - /// representing ||.||. - /// @return pair of (coef, var) vectors - /// so that the cone is ... >= sqrt(sum{ (coef_i * var*i)^2 }) - ConeArgs CheckNorm2(int res_var) { - if (MC().HasInitExpression(res_var)) { - auto init_expr = MC().GetInitExpression(res_var); - if (MC().template IsConInfoType(init_expr)) - return CheckNorm2_Pow(init_expr, res_var); - if (MC().template IsConInfoType(init_expr)) - return CheckNorm2_Abs(init_expr, res_var); - } - return {}; - } - - /// Check if the variable is defined by an expression + /// Result vars of expressions to un-use. + std::vector res_vars_to_delete_; + /// operator bool + operator bool() const + { assert(check()); return !coefs_.empty(); } + /// size() + size_t size() const + { assert(check()); return coefs_.size(); } + /// check() + bool check() const { return coefs_.size()==vars_.size(); } + }; + + /// Check if the variable is defined by an expression + /// representing ||.||. + /// @return pair of (coef, var) vectors + /// so that the cone is ... >= sqrt(sum{ (coef_i * var*i)^2 }) + ConeArgs CheckNorm2(int res_var) { + if (MC().HasInitExpression(res_var)) { + auto init_expr = MC().GetInitExpression(res_var); + if (MC().template IsConInfoType(init_expr)) + return CheckNorm2_Pow(init_expr, res_var); + if (MC().template IsConInfoType(init_expr)) + return CheckNorm2_Abs(init_expr, res_var); + } + return {}; + } + + /// Check if the variable is defined by an expression /// representing sqrt(c1 * x1^2 + ... [ + const ]). - template - ConeArgs CheckNorm2_Pow(const ConInfo& ci, int res_var) { - const auto& con_pow = MC().template - GetConstraint(ci); - const auto arg_pow = con_pow.GetArguments()[0]; - if (0.5 == con_pow.GetParameters()[0] && // sqrt(arg_pow) - MC().HasInitExpression(arg_pow)) { - auto ie_pow = MC().GetInitExpression(arg_pow); - if (MC().template // arg_pow := QFC(...) - IsConInfoType(ie_pow)) { - const auto& con_qdc = MC().template - GetConstraint(ie_pow); - const auto& args_qdc = con_qdc.GetArguments(); + template + ConeArgs CheckNorm2_Pow(const ConInfo& ci, int res_var) { + const auto& con_pow = MC().template + GetConstraint(ci); + const auto arg_pow = con_pow.GetArguments()[0]; + if (0.5 == con_pow.GetParameters()[0] && // sqrt(arg_pow) + MC().HasInitExpression(arg_pow)) { + auto ie_pow = MC().GetInitExpression(arg_pow); + if (MC().template // arg_pow := QFC(...) + IsConInfoType(ie_pow)) { + const auto& con_qdc = MC().template + GetConstraint(ie_pow); + const auto& args_qdc = con_qdc.GetArguments(); // Allowing a nonneg constant term by adding - // an auxiliary variable^2 + // an auxiliary variable^2 if (0.0 <= args_qdc.constant_term() && - args_qdc.GetBody().GetLinTerms().empty()) { - const auto& qpterms = args_qdc.GetBody().GetQPTerms(); - for (auto i=qpterms.size(); i--; ) { - if (qpterms.coef(i) < 0.0 || - qpterms.var1(i) != qpterms.var2(i)) - return {}; - } - ConeArgs result; - result.coefs_ = qpterms.coefs(); - for (auto& c: result.coefs_) - c = std::sqrt(c); - result.vars_ = qpterms.vars1(); + args_qdc.GetBody().GetLinTerms().empty()) { + const auto& qpterms = args_qdc.GetBody().GetQPTerms(); + for (auto i=qpterms.size(); i--; ) { + if (qpterms.coef(i) < 0.0 || + qpterms.var1(i) != qpterms.var2(i)) + return {}; + } + ConeArgs result; + result.coefs_ = qpterms.coefs(); + for (auto& c: result.coefs_) + c = std::sqrt(c); + result.vars_ = qpterms.vars1(); result.const_term = args_qdc.constant_term(); - result.res_vars_to_delete_ = { res_var, arg_pow }; - return result; - } - } - } - return {}; - } - - /// Check if the variable is defined by an expression - /// representing abs( c1*x1 ). - template - ConeArgs CheckNorm2_Abs(const ConInfo& ci, int res_var) { - ConeArgs result; - const auto& con_abs = MC().template - GetConstraint(ci); - const auto arg_abs = con_abs.GetArguments()[0]; - result.coefs_ = { 1.0 }; - result.vars_ = { arg_abs }; - result.res_vars_to_delete_ = { res_var }; - return result; - } - - /// Check if the variable is defined by sqrt(xN*xM) with - /// xN, xM >= 0. - ConeArgs CheckSqrtXnXmNonneg(int res_var) { - ConeArgs result; - if (auto pConPow = MC().template - GetInitExpressionOfType(res_var)) { - if (0.5 == pConPow->GetParameters()[0]) { // sqrt(arg_pow) - const auto arg_pow = pConPow->GetArguments()[0]; - if (auto pConQfc = MC().template + result.res_vars_to_delete_ = { res_var, arg_pow }; + return result; + } + } + } + return {}; + } + + /// Check if the variable is defined by an expression + /// representing abs( c1*x1 ). + template + ConeArgs CheckNorm2_Abs(const ConInfo& ci, int res_var) { + ConeArgs result; + const auto& con_abs = MC().template + GetConstraint(ci); + const auto arg_abs = con_abs.GetArguments()[0]; + result.coefs_ = { 1.0 }; + result.vars_ = { arg_abs }; + result.res_vars_to_delete_ = { res_var }; + return result; + } + + /// Check if the variable is defined by sqrt(xN*xM) with + /// xN, xM >= 0. + ConeArgs CheckSqrtXnXmNonneg(int res_var) { + ConeArgs result; + if (auto pConPow = MC().template + GetInitExpressionOfType(res_var)) { + if (0.5 == pConPow->GetParameters()[0]) { // sqrt(arg_pow) + const auto arg_pow = pConPow->GetArguments()[0]; + if (auto pConQfc = MC().template GetInitExpressionOfType< // arg_pow := QuadExpr - QuadraticFunctionalConstraint>(arg_pow)) { - const auto& qe = pConQfc->GetArguments(); - if (0.0 == std::fabs(qe.constant_term()) && - qe.GetBody().GetLinTerms().empty() && - 1 == qe.GetBody().GetQPTerms().size()) { - const auto& qpterms = qe.GetBody().GetQPTerms(); + QuadraticFunctionalConstraint>(arg_pow)) { + const auto& qe = pConQfc->GetArguments(); + if (0.0 == std::fabs(qe.constant_term()) && + qe.GetBody().GetLinTerms().empty() && + 1 == qe.GetBody().GetQPTerms().size()) { + const auto& qpterms = qe.GetBody().GetQPTerms(); if (qpterms.coef(0) >= 0.0 && // don't check if diff vars - MC().lb(qpterms.var1(0)) >= 0.0 && - MC().lb(qpterms.var2(0)) >= 0.0) { - result.coefs_ = {qpterms.coef(0), 1.0}; - result.vars_ = {qpterms.var1(0), qpterms.var2(0)}; - result.res_vars_to_delete_ = {res_var, arg_pow}; - return result; - } - } + MC().lb(qpterms.var1(0)) >= 0.0 && + MC().lb(qpterms.var2(0)) >= 0.0) { + result.coefs_ = {qpterms.coef(0), 1.0}; + result.vars_ = {qpterms.var1(0), qpterms.var2(0)}; + result.res_vars_to_delete_ = {res_var, arg_pow}; + return result; + } + } } else { // It's just some variable if (MC().lb(arg_pow) >= 0.0) { result.coefs_ = {1.0}; @@ -434,18 +434,18 @@ class Convert1QC : public MCKeeper { return result; } } - } - } - return result; - } + } + } + return result; + } - /// Continue processing a linear constraint x>=y, + /// Continue processing a linear constraint x>=y, /// if y := ||.|| and x is sqrt(xN*xM) or sqrt(k*x). - /// Rotated Qcone. - /// @return true iff converted. - bool ContinueRotatedSOC( - const LinTerms& lint, int iX, int iY, - const ConeArgs& lhs_args, const ConeArgs& rhs_args) { + /// Rotated Qcone. + /// @return true iff converted. + bool ContinueRotatedSOC( + const LinTerms& lint, int iX, int iY, + const ConeArgs& lhs_args, const ConeArgs& rhs_args) { assert(2>=lhs_args.size()); assert(1<=lhs_args.size()); assert(rhs_args); @@ -453,9 +453,9 @@ class Convert1QC : public MCKeeper { (0.0 != std::fabs(rhs_args.const_term)); std::vector c(2+rhs_args.size()+constNZ); std::vector x(c.size()); - auto coefX_abs = std::fabs(lint.coef(iX)); - c[0] = 0.5 * lhs_args.coefs_[0] * coefX_abs; - x[0] = lhs_args.vars_[0]; + auto coefX_abs = std::fabs(lint.coef(iX)); + c[0] = 0.5 * lhs_args.coefs_[0] * coefX_abs; + x[0] = lhs_args.vars_[0]; if (2==lhs_args.size()) { // x is sqrt(xN*xM) c[1] = lhs_args.coefs_[1] * coefX_abs; x[1] = lhs_args.vars_[1]; @@ -463,58 +463,58 @@ class Convert1QC : public MCKeeper { c[1] = lhs_args.coef_extra * coefX_abs; x[1] = MC().MakeFixedVar(1.0); } - auto coefY_abs = std::fabs(lint.coef(iY)); - for (size_t iPush=0; iPush= cY*y, - /// if y := ||.|| and x is not sqrt(xN*xM). - /// Standard Qcone. - /// @return true iff converted. - bool ContinueStdSOC( + /// if y := ||.|| and x is not sqrt(xN*xM). + /// Standard Qcone. + /// @return true iff converted. + bool ContinueStdSOC( double cX, int vX, double cY, - const ConeArgs& rhs_args) { + const ConeArgs& rhs_args) { const size_t constNZ = (0.0 != std::fabs(rhs_args.const_term)); std::vector x(rhs_args.size()+1+constNZ); std::vector c(rhs_args.size()+1+constNZ); x[0] = vX; c[0] = std::fabs(cX); auto coefY_abs = std::fabs(cY); - for (size_t iPush=0; iPush { const auto rhsNon0 = (std::fabs(rhs)!=0.0); std::vector x(qpterms.size()+1+lint.size()+rhsNon0); std::vector c(qpterms.size()+1+lint.size()+rhsNon0); - size_t iPush=1; + size_t iPush=1; if (lint.size()) { x[0] = lint.var(0); x[1] = MC().MakeFixedVar( 1.0 ); c[0] = std::fabs(lint.coef(0)) / 2.0; c[1] = 1.0; // it's 2xy >= z^2 + ... } - for (int i=0; i<(int)qpterms.size(); ++i) { - if (i==iDiffVars) { - x[0] = qpterms.var1(i); - x[1] = qpterms.var2(i); - c[0] = std::fabs(qpterms.coef(i)); + for (int i=0; i<(int)qpterms.size(); ++i) { + if (i==iDiffVars) { + x[0] = qpterms.var1(i); + x[1] = qpterms.var2(i); + c[0] = std::fabs(qpterms.coef(i)); c[1] = 0.5; // it's 2xy >= z^2 + ... - } else { - ++iPush; - x.at(iPush) = qpterms.var1(i); - c.at(iPush) = std::sqrt(std::fabs(qpterms.coef(i))); - } - } + } else { + ++iPush; + x.at(iPush) = qpterms.var1(i); + c.at(iPush) = std::sqrt(std::fabs(qpterms.coef(i))); + } + } if (rhsNon0) { // Add ... >= ... + |rhs| ++iPush; x.at(iPush) = MC().MakeFixedVar( 1.0 ); c.at(iPush) = std::sqrt( std::fabs(rhs) ); } assert(iPush == x.size()-1); - MC().AddConstraint( - RotatedQuadraticConeConstraint( - std::move(x), std::move(c))); - return true; - } + MC().AddConstraint( + RotatedQuadraticConeConstraint( + std::move(x), std::move(c))); + return true; + } /// Add standard cone from pure-quadratic constraint. /// @param iSame1: if <0, use new fixed variable. @@ -570,16 +570,16 @@ class Convert1QC : public MCKeeper { x[0] = MC().MakeFixedVar( 1.0 ); } size_t iPush=0; - for (int i=0; i<(int)qpterms.size(); ++i) { - if (i==iSame1) { - x[0] = qpterms.var1(i); - c[0] = std::sqrt(std::fabs(qpterms.coef(i))); - } else { - ++iPush; - x.at(iPush) = qpterms.var1(i); - c.at(iPush) = std::sqrt(std::fabs(qpterms.coef(i))); - } - } + for (int i=0; i<(int)qpterms.size(); ++i) { + if (i==iSame1) { + x[0] = qpterms.var1(i); + c[0] = std::sqrt(std::fabs(qpterms.coef(i))); + } else { + ++iPush; + x.at(iPush) = qpterms.var1(i); + c.at(iPush) = std::sqrt(std::fabs(qpterms.coef(i))); + } + } if (fNewVar && fRhs2) { ++iPush; c.at(iPush) = std::sqrt(std::fabs(rhs2)); @@ -587,13 +587,13 @@ class Convert1QC : public MCKeeper { } assert(iPush == c.size()-1); MC().AddConstraint( - QuadraticConeConstraint( - std::move(x), std::move(c))); - return true; - } + QuadraticConeConstraint( + std::move(x), std::move(c))); + return true; + } - /// Retrieve the MC - using MCKeeper::MC; + /// Retrieve the MC + using MCKeeper::MC; }; @@ -604,149 +604,149 @@ class Convert1QC : public MCKeeper { template class Convert1ExpC : public MCKeeper { public: - /// Constructor - Convert1ExpC(MCType& mc) : MCKeeper(mc) { } + /// Constructor + Convert1ExpC(MCType& mc) : MCKeeper(mc) { } protected: - /// Characterize target subexpression. - template - struct SubexprTraits { - /// Construct - SubexprTraits( - std::array c={}, - std::array v={}, - std::vector v2d={}) - : coefs_(c), vars_(v), vars2del_(v2d) { } - /// Valid? - operator bool() const { return valid_; } - /// Variables to use. - /// Coefs - std::array coefs_; - /// -1 means 'use constant 1'. - std::array vars_; - /// Vars to delete - std::vector vars2del_; - /// Valid? - bool valid_ {false}; - }; + /// Characterize target subexpression. + template + struct SubexprTraits { + /// Construct + SubexprTraits( + std::array c={}, + std::array v={}, + std::vector v2d={}) + : coefs_(c), vars_(v), vars2del_(v2d) { } + /// Valid? + operator bool() const { return valid_; } + /// Variables to use. + /// Coefs + std::array coefs_; + /// -1 means 'use constant 1'. + std::array vars_; + /// Vars to delete + std::vector vars2del_; + /// Valid? + bool valid_ {false}; + }; - /// LHS of ax >= by * exp( cz / by ) - using LhsTraits = SubexprTraits<1>; - /// RHS - using RhsTraits = SubexprTraits<2>; - - /// DoRun. Body: quadratic. - /// - /// Accept non-(+-1) coefficients. - /// - /// @param body: quadratic constraint body. - /// @param sens: -1 for <=, 1 for >=. - /// @param rhs: constraint's rhs. - /// - /// @return true iff converted. - bool DoRun(const QuadAndLinTerms& body, - int sens, double rhs) { - assert((sens==1 || sens==-1) && "sens 1 or -1 only"); - if (1 == body.GetQPTerms().size()) { + /// LHS of ax >= by * exp( cz / by ) + using LhsTraits = SubexprTraits<1>; + /// RHS + using RhsTraits = SubexprTraits<2>; + + /// DoRun. Body: quadratic. + /// + /// Accept non-(+-1) coefficients. + /// + /// @param body: quadratic constraint body. + /// @param sens: -1 for <=, 1 for >=. + /// @param rhs: constraint's rhs. + /// + /// @return true iff converted. + bool DoRun(const QuadAndLinTerms& body, + int sens, double rhs) { + assert((sens==1 || sens==-1) && "sens 1 or -1 only"); + if (1 == body.GetQPTerms().size()) { const auto& lt = body.GetLinTerms(); // But we have c0*x + c1*y*zz <=> 0 const auto& qt = body.GetQPTerms(); LhsTraits lhst; if (0.0 == std::fabs(rhs) && // rhs==0: 1 == lt.size()) // ax >= by exp(cz / by) ? - lhst = ClassifyLhs(lt.coef(0)*sens, lt.var(0)); + lhst = ClassifyLhs(lt.coef(0)*sens, lt.var(0)); else if (0 == lt.size()) // const >= by exp(cz / by) ? - lhst = ClassifyLhs(-rhs*sens, -1); + lhst = ClassifyLhs(-rhs*sens, -1); if (lhst) { - if (auto rhst = ClassifyRhsQuadr( - -qt.coef(0)*sens, qt.var1(0), qt.var2(0))) + if (auto rhst = ClassifyRhsQuadr( + -qt.coef(0)*sens, qt.var1(0), qt.var2(0))) return AddExpCone(lhst, rhst); - else if (auto rhst = ClassifyRhsQuadr( - -qt.coef(0)*sens, qt.var2(0), qt.var1(0))) + else if (auto rhst = ClassifyRhsQuadr( + -qt.coef(0)*sens, qt.var2(0), qt.var1(0))) return AddExpCone(lhst, rhst); } - } - return false; } + return false; + } - /// DoRun. - /// Body: linear, i.e., var1 (or const) >= var2 (or const). - /// In particular, the following cases. - /// /// - /// Considering const>=0. - /// Accept non-(+-1) coefficients. - /// - /// @param body: linear constraint body. - /// @param sens: -1 for <=, 1 for >=. - /// @param rhs: constraint's rhs. - /// - /// @return true iff converted. - bool DoRun(const LinTerms& lt, - int sens, double rhs) { - assert((sens==1 || sens==-1) && "sens 1 or -1 only"); - if (0.0 == std::fabs(rhs) && // rhs==0: y==const - 2 == lt.size()) { // ax >= b exp(cz / b) + /// DoRun. + /// Body: linear, i.e., var1 (or const) >= var2 (or const). + /// In particular, the following cases. + /// /// + /// Considering const>=0. + /// Accept non-(+-1) coefficients. + /// + /// @param body: linear constraint body. + /// @param sens: -1 for <=, 1 for >=. + /// @param rhs: constraint's rhs. + /// + /// @return true iff converted. + bool DoRun(const LinTerms& lt, + int sens, double rhs) { + assert((sens==1 || sens==-1) && "sens 1 or -1 only"); + if (0.0 == std::fabs(rhs) && // rhs==0: y==const + 2 == lt.size()) { // ax >= b exp(cz / b) if (auto rhst = ClassifyRhsLin(-lt.coef(0)*sens, lt.var(0))) { - // c2*x >= -c1*exp() - if (auto lhst = ClassifyLhs(lt.coef(1)*sens, lt.var(1))) { + // c2*x >= -c1*exp() + if (auto lhst = ClassifyLhs(lt.coef(1)*sens, lt.var(1))) { return AddExpCone(lhst, rhst); - } + } } else if (auto rhst = ClassifyRhsLin(-lt.coef(1)*sens, lt.var(1))) { - // c1*x >= -c2*exp() - if (auto lhst = ClassifyLhs(lt.coef(0)*sens, lt.var(0))) { + // c1*x >= -c2*exp() + if (auto lhst = ClassifyLhs(lt.coef(0)*sens, lt.var(0))) { return AddExpCone(lhst, rhst); - } + } } - } else if (1 == lt.size()) { // const >= b exp(cz / b) + } else if (1 == lt.size()) { // const >= b exp(cz / b) if (auto lhst = ClassifyLhs(-rhs*sens, -1)) { - if (auto rhst = ClassifyRhsLin(-lt.coef(0)*sens, lt.var(0))) { + if (auto rhst = ClassifyRhsLin(-lt.coef(0)*sens, lt.var(0))) { return AddExpCone(lhst, rhst); - } - } } - return false; + } } + return false; + } protected: - LhsTraits ClassifyLhs(double a, int x) { - LhsTraits result {{a}, {x}}; - if ((x<0 && a>=0.0) || // Lhs is const, >=0 - (x>=0 && a>=0.0 && MC().lb(x)>=0.0) || // c*v >= 0 - (x>=0 && a<=0.0 && MC().ub(x)>=0.0)) - result.valid_ = true; - return result; - } + LhsTraits ClassifyLhs(double a, int x) { + LhsTraits result {{a}, {x}}; + if ((x<0 && a>=0.0) || // Lhs is const, >=0 + (x>=0 && a>=0.0 && MC().lb(x)>=0.0) || // c*v >= 0 + (x>=0 && a<=0.0 && MC().ub(x)>=0.0)) + result.valid_ = true; + return result; + } - /// The RHS is just b * (v=exp(z))? - RhsTraits ClassifyRhsLin(double b, int v) { - assert(v>=0); - RhsTraits result {{b, b}, {-1} /*constant*/}; - if (b>=0.0) { - if (auto pConExp = MC().template - GetInitExpressionOfType(v)) { - result.vars_[1] = pConExp->GetArguments()[0]; // the z - result.vars2del_ = {v}; // delete v - result.valid_ = true; - } - } - return result; + /// The RHS is just b * (v=exp(z))? + RhsTraits ClassifyRhsLin(double b, int v) { + assert(v>=0); + RhsTraits result {{b, b}, {-1} /*constant*/}; + if (b>=0.0) { + if (auto pConExp = MC().template + GetInitExpressionOfType(v)) { + result.vars_[1] = pConExp->GetArguments()[0]; // the z + result.vars2del_ = {v}; // delete v + result.valid_ = true; + } } + return result; + } - /// The RHS is by * (v1 = exp( v2 = cz / by )) - /// or by * (v1 = exp( v2 = c / by ))? - /// In our case, c/by is stored as [constvar=c/b]/y - RhsTraits ClassifyRhsQuadr(double b, int y, int v1) { - assert(y>=0 && v1>=0); - RhsTraits result {{b, b}}; - if ((b>=0.0 && MC().lb(y)>=0.0) || // by >= 0 - (b<=0.0 && MC().ub(y)<=0.0)) { - if (auto pConExp = MC().template - GetInitExpressionOfType(v1)) { - int v2 = pConExp->GetArguments()[0]; // the z - if (auto pConDiv = MC().template - GetInitExpressionOfType(v2)) { + /// The RHS is by * (v1 = exp( v2 = cz / by )) + /// or by * (v1 = exp( v2 = c / by ))? + /// In our case, c/by is stored as [constvar=c/b]/y + RhsTraits ClassifyRhsQuadr(double b, int y, int v1) { + assert(y>=0 && v1>=0); + RhsTraits result {{b, b}}; + if ((b>=0.0 && MC().lb(y)>=0.0) || // by >= 0 + (b<=0.0 && MC().ub(y)<=0.0)) { + if (auto pConExp = MC().template + GetInitExpressionOfType(v1)) { + int v2 = pConExp->GetArguments()[0]; // the z + if (auto pConDiv = MC().template + GetInitExpressionOfType(v2)) { int z = pConDiv->GetArguments()[0]; // the z0 int y0 = pConDiv->GetArguments()[1]; if (y0==y) { // v2 = z/y @@ -754,31 +754,31 @@ class Convert1ExpC : public MCKeeper { result.vars2del_ = {v1, v2}; // delete v1, v2 result.valid_ = true; } else if (auto pConLin = MC().template - GetInitExpressionOfType(y0)) { + GetInitExpressionOfType(y0)) { const auto& ae = pConLin->GetAffineExpr(); if (0.0 == std::fabs(ae.constant_term() && 1==ae.GetBody().size())) { - const auto& body = ae.GetBody(); // can we ever get here? - if (y == body.var(0)) { // v2 = z / (c1*y) - result.coefs_ = {b, b/body.coef(0)}; - result.vars_ = {y, z}; - result.vars2del_ = {v1, v2, y0}; // delete - result.valid_ = true; - } + const auto& body = ae.GetBody(); // can we ever get here? + if (y == body.var(0)) { // v2 = z / (c1*y) + result.coefs_ = {b, b/body.coef(0)}; + result.vars_ = {y, z}; + result.vars2del_ = {v1, v2, y0}; // delete + result.valid_ = true; + } } } - } - } } - return result; + } } + return result; + } - bool AddExpCone(LhsTraits l, RhsTraits r) { - return true; - } + bool AddExpCone(LhsTraits l, RhsTraits r) { + return true; + } - /// Retrieve the MC - using MCKeeper::MC; + /// Retrieve the MC + using MCKeeper::MC; }; @@ -787,37 +787,37 @@ class Convert1ExpC : public MCKeeper { * if they are inequalities. */ template class CvtBase> + template class CvtBase> class Convert1, CvtBase > : - public CvtBase { + AlgebraicConstraint, CvtBase > : + public CvtBase { public: - /// Constructor + /// Constructor Convert1(MCType& mc) : CvtBase(mc) { } - /// Constraint type - using ConType = AlgebraicConstraint; - - /// Run - bool operator()(const ConType& con, int i) { - bool lbf = con.lb() >= MC().PracticallyMinusInf(); - bool ubf = con.ub() <= MC().PracticallyInf(); - if (lbf+ubf == 1) { - // Autolink is to pass back duals... Simplify? - // E.g., automatically in ForEachActive? - pre::AutoLinkScope auto_link_scope{ - MC(), MC().template SelectValueNode(i) - }; - return - CvtBase::DoRun(con.GetBody(), - lbf ? 1 : -1, - lbf ? con.lb() : con.ub()); - } - return false; - } - - /// Retrieve the MC - using MCKeeper::MC; + /// Constraint type + using ConType = AlgebraicConstraint; + + /// Run + bool operator()(const ConType& con, int i) { + bool lbf = con.lb() >= MC().PracticallyMinusInf(); + bool ubf = con.ub() <= MC().PracticallyInf(); + if (lbf+ubf == 1) { + // Autolink is to pass back duals... Simplify? + // E.g., automatically in ForEachActive? + pre::AutoLinkScope auto_link_scope{ + MC(), MC().template SelectValueNode(i) + }; + return + CvtBase::DoRun(con.GetBody(), + lbf ? 1 : -1, + lbf ? con.lb() : con.ub()); + } + return false; + } + + /// Retrieve the MC + using MCKeeper::MC; }; @@ -825,31 +825,31 @@ class Convert1 class CvtBase> + template class CvtBase> class Convert1 >, CvtBase > : - public CvtBase { + AlgebraicConstraint< Body, AlgConRhs >, CvtBase > : + public CvtBase { public: - /// Constructor + /// Constructor Convert1(MCType& mc) : CvtBase(mc) { } - /// The constraint type - using ConType = - AlgebraicConstraint< Body, AlgConRhs >; - - /// Run - bool operator()(const ConType& con, int i) { - static_assert (sens==1 || sens==-1, "sens 1 or -1 only"); - pre::AutoLinkScope auto_link_scope{ - MC(), MC().template SelectValueNode(i) - }; - return - CvtBase::DoRun(con.GetBody(), - sens, con.rhs()); - } - - /// Retrieve the MC - using MCKeeper::MC; + /// The constraint type + using ConType = + AlgebraicConstraint< Body, AlgConRhs >; + + /// Run + bool operator()(const ConType& con, int i) { + static_assert (sens==1 || sens==-1, "sens 1 or -1 only"); + pre::AutoLinkScope auto_link_scope{ + MC(), MC().template SelectValueNode(i) + }; + return + CvtBase::DoRun(con.GetBody(), + sens, con.rhs()); + } + + /// Retrieve the MC + using MCKeeper::MC; }; } // namespace mp