SOCP: move QP terms into aux SOC #229 #192

When any cones found, move QP objective terms into an aux SOC, in the easy case

Specially for Mosek

include/mp/flat/expr_quadratic.h

@@ -107,6 +107,13 @@ class QuadTerms {
   /// Sort and eliminate duplicates
   void sort_terms();
 
+  /// Clear
+  void clear() {
+    coefs_.clear();
+    vars1_.clear();
+    vars2_.clear();
+  }
+
   /// Test equality
   bool equals(const QuadTerms& qt) const {
     return *this == qt;

include/mp/flat/obj_std.h

@@ -26,8 +26,10 @@ class LinearObjective {
     name_(std::move(nm)){ }
   /// Get sense
   obj::Type obj_sense() const { return sense_; }
-  /// Get lin terms
+  /// Get lin terms, const
   const LinTerms& GetLinTerms() const { return lt_; }
+  /// Get lin terms
+  LinTerms& GetLinTerms() { return lt_; }
   /// Get N terms
   int num_terms() const { assert(check()); return (int)lt_.size(); }
   /// Validate
@@ -55,8 +57,10 @@ class QuadraticObjective : public LinearObjective {
   QuadraticObjective(LinearObjective&& lc, QuadTerms&& qt) :
     LinearObjective(std::move(lc)), qt_(std::move(qt)) { sort_qp_terms(); }
 
-  /// Get QP terms
+  /// Get QP terms, const
   const QuadTerms& GetQPTerms() const { return qt_; }
+  /// Get QP terms
+  QuadTerms& GetQPTerms() { return qt_; }
 
   /// Sort QP terms
   void sort_qp_terms() {

include/mp/flat/redef/conic/cones.h

@@ -85,6 +85,7 @@ class ConicConverter : public MCKeeper<MCType> {
       RunQConesFromNonQC();
     if (MC().IfPassExpCones())
       RunExpCones();
+    TryReformulateQPObjective();
     if (IfNeedSOCP2QC())
       SetupSOCP2QC();
     WarnOnMix();
@@ -119,11 +120,62 @@ class ConicConverter : public MCKeeper<MCType> {
     }
   }
 
+  /// Ideally also for exp/log terms
+  void TryReformulateQPObjective() {
+    if (HasAnyCones()
+        && MC().HasQPObjective() && 1==MC().num_objs()) {
+      // See if the QP objective is SOCP-easy
+      auto& obj = MC().get_objectives().at(0);
+      auto& qp = obj.GetQPTerms();
+      auto objsns = (obj::MAX==obj.obj_sense()) ? 1.0 : -1.0;
+      for (auto i=qp.size(); i--; ) {
+        if (qp.coef(i) * objsns > 0.0
+            || qp.var1(i) != qp.var2(i))       // not x1[i]==x2[i]
+          return;
+      }
+      // Set up AutoLink
+      auto obj_src =              // source value node for this obj
+          MC().GetValuePresolver().GetSourceNodes().GetObjValues()().Select(0);
+      pre::AutoLinkScope<MCType> auto_link_scope{ MC(), obj_src };
+      // Aux vars
+      int z = (int)MC().AddVar(0.0, MC().Infty());
+      int z1 = (int)MC().MakeFixedVar(1.0);
+      // Add cone
+      std::vector<double> c = {{1.0, 0.5}};
+      c.insert(c.end(), qp.coefs().begin(), qp.coefs().end());
+      if (objsns > 0)       // negate coefs if maximizing
+        for (size_t i=2; i<c.size(); ++i)
+          c[i] = -c[i];
+      std::vector<int> x = {{z, z1}};
+      x.insert(x.end(), qp.vars1().begin(), qp.vars1().end());
+      MC().AddConstraint(
+            RotatedQuadraticConeConstraint(
+              std::move(x), std::move(c)));
+      // Add linear term, remove QP terms
+      qp.clear();
+      obj.GetLinTerms().add_term(-objsns, z);
+    }
+  }
+
+  /// Any cones at all
+  bool HasAnyCones() const {
+    return MC().GetNumberOfAddable((QuadraticConeConstraint*)0)>0 ||
+        MC().GetNumberOfAddable((RotatedQuadraticConeConstraint*)0)>0 ||
+        HasAnyNonSOCPCones();
+  }
+
+  /// any non-SOCP cones?
+  bool HasAnyNonSOCPCones() const {
+    return MC().GetNumberOfAddable((ExponentialConeConstraint*)0)>0 ||
+        MC().GetNumberOfAddable((PowerConeConstraint*)0)>0 ||
+        MC().GetNumberOfAddable((GeometricConeConstraint*)0)>0;
+  }
+
   bool IfNeedSOCP2QC() {
     return
         MC().SOCP2QCMode() >= 2     // compulsory
         || (1 == MC().SOCP2QCMode()
-            && 0 == MC().NumExpConesRecognized()
+            && !HasAnyNonSOCPCones()
                // Might also have SOCP from sqrt(), abs().
                // Some QC -> SOCP but not all, even if 0 succeeded.
             && (MC().NumQC2SOCPAttempted() > MC().NumQC2SOCPSucceeded()
@@ -137,14 +189,10 @@ class ConicConverter : public MCKeeper<MCType> {
 
   void WarnOnMix() {
     if ( !MC().ModelAPICanMixConicQCAndQC()) {  // cannot mix
-      if ((MC().NumExpConesRecognized()         // exp cones
+      if ((HasAnyCones()                        // exp cones
            && (MC().NumQC2SOCPAttempted() > MC().NumQC2SOCPSucceeded()
                || MC().HasQPObjective()))       // and quadratics left in
-          ||                // Some QC -> SOCP but not all
-          ((MC().NumQC2SOCPAttempted() > MC().NumQC2SOCPSucceeded()
-            || MC().HasQPObjective())           // or a quadratic obj
-           && MC().NumQC2SOCPSucceeded()        // Warn only if some succeeded
-           && !MC().IfConvertSOCP2QC())         // and not decided to convert
+           && !MC().IfConvertSOCP2QC()          // and not decided to convert
           ) {
         MC().AddWarning("Mix QC+cones",
                         "Not all quadratic constraints could "
@@ -153,7 +201,7 @@ class ConicConverter : public MCKeeper<MCType> {
                         "additionally, further convertion back to QC\n"
                         "not desired (option cvt:socp2qc) or other cone types present;\n"
                         "solver might not accept the model.\n"
-                        "Try to express all SOCP cones in standard forms,\n"
+                        "Try to express all cones in standard forms,\n"
                         "not in the objective.\n"
                         "See mp.ampl.com/model-guide.html.");
       }

test/end2end/cases/categorized/fast/conic/modellist.json

@@ -28,13 +28,13 @@
     "name" : "TotalVariation2D",
     "objective" : 9.569949819774543,
     "files": ["totalvariation2d.mod", "totalvar2d_4e_4.dat"],
-    "tags" : ["socp", "socp_hard_to_recognize"]
+    "tags" : ["socp"]
   },
   {
     "name" : "TotalVariation2D_Quadr",
     "objective" : 9.569949819774543,
     "files": ["totalvariation2d_quad.mod", "totalvar2d_4e_4.dat"],
-    "tags" : ["socp", "socp_hard_to_recognize"]
+    "tags" : ["socp"]
   },
   {
     "name" : "socp_06_rsoc_constants",
@@ -77,6 +77,14 @@
       "solve_result_num": 0
     }
   },
+  {
+    "name" : "socp_12__issue229_2",
+    "tags" : ["socp"],
+    "objective" : -1607.1428,
+    "values": {
+      "solve_result_num": 0
+    }
+  },
   {
     "name" : "expcones_01__plain",
     "objective" : 0.7821882953,

test/end2end/cases/categorized/fast/conic/socp_12__issue229_2.mod

@@ -0,0 +1,7 @@
+var x {1..7} >= 0;
+var y >= 0;
+
+maximize obj: sum {j in 1..7} -x[j]^2 - y^2;
+
+subj to conQ: sum {j in 1..7} x[j]^2 <= y^2;
+subj to conL: sum {j in 1..7} x[j] = 75;

test/end2end/scripts/python/Solver.py

@@ -729,6 +729,7 @@ def __init__(self, exeName, timeout=None, nthreads=None,
                  ModelTags.continuous, ModelTags.integer, ModelTags.binary,
                  ModelTags.sos,
                  ModelTags.quadratic, ModelTags.quadratic_obj,
+                 ModelTags.socp,
 
                  ModelTags.return_mipgap, ModelTags.multiobj,