SOCP #192: const >= ||.||

Recognize std SOCs with constant lhs, by substituting a fixed variable

include/mp/flat/converter.h

@@ -182,6 +182,8 @@ class FlatConverter :
 	}
 
 	/// Unuse result variable.
+  /// Actually this is to 'unuse' the init expression
+  /// - might change naming.
 	/// Throw if already not used.
 	void DecrementVarUsage(int v) {
 		assert(VarUsageRef(v)>0);

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

@@ -170,7 +170,7 @@ class Convert1QC : public MCKeeper<MCType> {
 
 	/// DoRun. Body: linear.
 	///
-	/// Currently considering rhs=0.
+  /// Currently considering rhs=0 or rhs>0.
 	/// Accept non-(+-1) coefficients.
 	///
 	/// @param body: linear constraint body.
@@ -194,9 +194,21 @@ class Convert1QC : public MCKeeper<MCType> {
 				if (auto lhs_args = CheckSqrtXnXmNonneg(lint.var(iX)))
 					return ContinueRotatedSOC(lint, iX, iY,
 																		lhs_args, rhs_args);
-				return ContinueStdSOC(lint, iX, iY, 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)
+               0.0 >= lint.coef(0)*sens) {   // In the other case,
+                                             // it's -k*sqrt(...) <= rhs(>0), k>0,
+                                             // which is always true TODO
+      if (auto rhs_args = CheckNorm2(lint.var(0))) {
+        return ContinueStdSOC(1.0,
+                              MC().MakeFixedVar(std::fabs(rhs)),
+                              lint.coef(0), rhs_args);
+      }
+    }
 		return false;
 	}
 
@@ -347,18 +359,18 @@ class Convert1QC : public MCKeeper<MCType> {
 		return true;
 	}
 
-	/// Continue processing a linear constraint x>=y,
+  /// Continue processing a linear constraint cX*x >= cY*y,
 	/// if y := ||.|| and x is not sqrt(xN*xM).
 	/// Standard Qcone.
 	/// @return true iff converted.
 	bool ContinueStdSOC(
-			const LinTerms& lint, int iX, int iY,
+      double cX, int vX, double cY,
 			const ConeArgs& rhs_args) {
 		std::vector<int> x(rhs_args.size()+1);
 		std::vector<double> c(rhs_args.size()+1);
-		x[0] = lint.var(iX);
-		c[0] = std::fabs(lint.coef(iX));
-		auto coefY_abs = std::fabs(lint.coef(iY));
+    x[0] = vX;
+    c[0] = std::fabs(cX);
+    auto coefY_abs = std::fabs(cY);
 		for (size_t iPush=0; iPush<rhs_args.size(); ++iPush) {
 			x[iPush+1] = rhs_args.vars_[iPush];
 			c[iPush+1] = coefY_abs * rhs_args.coefs_[iPush];

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

@@ -23,5 +23,11 @@
     "name" : "socp_05_rqcSqrt_abs",
     "objective" : 3.623966621,
     "tags" : ["socp"]
+  },
+  {
+    "name" : "TotalVariation2D",
+    "objective" : 9.569949819774543,
+    "files": ["totalvariation2d.mod", "totalvar2d_4e_4.dat"],
+    "tags" : ["socp", "socp_hard_to_recognize"]
   }
 ]

test/end2end/cases/categorized/fast/conic/totalvar2d_4e_4.dat

@@ -0,0 +1,6 @@
+data;
+
+param n := 10;
+param m := 20;
+
+param sigmaVal := 0.004;

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

@@ -0,0 +1,35 @@
+## 2D Total Variation
+## This case study is based mainly on the paper by Goldfarb and Yin [GY05]
+## https://docs.mosek.com/latest/pythonfusion/case-studies-total-variation.html
+
+param n>1;
+param m>1;
+
+param f {i in 1..n, j in 1..m} default        ## signal + noise
+    1.0*(i+j-2)/(n+m) + Normal(0.0, 0.08);
+param sigmaVal default 0.0004;
+param sigma >=0.0 default sigmaVal*n*m;
+
+var u {1..n, 1..m} >=0.0, <=1.0;
+var delta_uf {1..n, 1..m};
+var t {1..n-1, 1..m-1};
+var delta_plus {1..n-1, 1..m-1, 1..2};
+
+minimize TotalVar:
+    sum {i in 1..n-1, j in 1..m-1} t[i, j];
+
+s.t. DefDeltaPlusUp {i in 1..n-1, j in 1..m-1}:
+    delta_plus[i, j, 1] == u[i+1, j] - u[i, j];
+
+s.t. DefDeltaPlusRight {i in 1..n-1, j in 1..m-1}:
+    delta_plus[i, j, 2] == u[i, j+1] - u[i, j];
+
+s.t. DefDeltaUF {i in 1..n, j in 1..m}:
+    delta_uf[i, j] == u[i, j] - f[i, j];
+
+s.t. VariationBound {i in 1..n-1, j in 1..m-1}:
+    -t[i, j] <= -sqrt(delta_plus[i, j, 1]^2 + delta_plus[i, j, 2]^2);
+
+s.t. Total2Norm:
+    sigma >= sqrt( sum {i in 1..n, j in 1..m} delta_uf[i, j]^2 );
+

test/end2end/scripts/python/Model.py

@@ -12,6 +12,7 @@ class ModelTags(enum.Enum):
     quadratic_obj = 2.5
     quadraticnonconvex = 3
     socp = 4
+    socp_hard_to_recognize = 4.1  ## For solvers recognizing from quadratics
     nonlinear = 5
     complementarity = 6
     arc = 7

test/end2end/scripts/python/Solver.py

@@ -757,7 +757,7 @@ def __init__(self, exeName, timeout=None, nthreads=None, otherOptions=None):
                  ModelTags.binary,
 
                  ModelTags.quadratic, ModelTags.quadratic_obj,
-                 ModelTags.socp,
+                 ModelTags.socp, ModelTags.socp_hard_to_recognize,
 
                  ModelTags.warmstart, ModelTags.mipstart,
                  ModelTags.return_mipgap, ModelTags.sens, ModelTags.sstatus}