Test is working for DR PCTSP variance constraint

include/pctsp/robust.hh

@@ -5,14 +5,14 @@
 
 #include "algorithms.hh"
 
-std::vector<std::pair<PCTSPvertex, PCTSPvertex>> DistRobustPrizeCollectingTsp(
+std::vector<std::pair<PCTSPvertex, PCTSPvertex>> solveDistRobustPrizeCollectingTsp(
     SCIP * scip,
     PCTSPgraph& graph,
-    EdgeCostMap& cost_mean_map,
-    EdgeCostMap& cost_var_map,
-    VertexPrizeMap& prize_map,
+    EdgeCostMap& costMeanMap,
+    EdgeCostMap& costSigmaMap,
+    VertexPrizeMap& prizeMap,
     PrizeNumberType& quota,
-    PCTSPvertex& root_vertex,
+    PCTSPvertex& rootVertex,
     std::string& name
 );
 
@@ -26,9 +26,9 @@ SCIP_RETCODE SCIPcreateConsSOC(
     SCIP*                 scip,               /**< SCIP data structure */
     SCIP_CONS**           cons,               /**< pointer to hold the created constraint */
     const char*           name,               /**< name of constraint */
-    int                   nvars,              /**< number of variables on left hand side of constraint (n) */
-    SCIP_VAR**            vars,               /**< array with variables on left hand side (x_i) */
-    SCIP_Real*            coefs,              /**< array with coefficients of left hand side variables (alpha_i), or NULL if all 1.0 */
+    int                   nlhsexprs,          /**< number of variables on left hand side of constraint (n) */
+    SCIP_EXPR**           lhsexprs,           /**< array of expressions on left hand side (x_i) */
+    SCIP_Real*            lhscoefs,           /**< array with coefficients of left hand side variables (alpha_i), or NULL if all 1.0 */
     SCIP_Real*            offsets,            /**< array with offsets of variables (beta_i), or NULL if all 0.0 */
     SCIP_Real             constant,           /**< constant on left hand side (gamma) */
     int                   nrhsvars,           /**< number of variables on right hand side of constraint */

src/robust.cpp

@@ -1,28 +1,114 @@
 /** Distributionally Robust Algorithms */
 
+#include <assert.h>
 #include <scip/cons_nonlinear.h>
 #include <scip/expr_sum.h>
 
 #include "pctsp/robust.hh"
 
-void addDistRobustCons() {
-
+void addDistRobustCons(
+    SCIP * scip,
+    PCTSPgraph& graph,
+    EdgeCostMap& costSigmaMap,
+    PCTSPedgeVariableMap& edgeVarMap
+) {
+    // initialize basic variables
+    int numEdges = boost::num_edges(graph);
+    int numVertices = boost::num_vertices(graph);
+
+    // setup the robust variance constraint
+    SCIP_CONS* robustCons;
+    double alpha = 1.0;
+    int nLhsExprs = numEdges + 1;
+    std::vector<double> lhsCoefs (nLhsExprs);
+    std::vector<SCIP_EXPR*> lhsExprs (nLhsExprs);
+
+    // iterate over each variable and make it into an expression
+    int i = 0;
+    for (auto const& [edge, var] : edgeVarMap) {
+        lhsCoefs[i] = sqrt(costSigmaMap[edge]); // sigma is coef for x variables
+        SCIP_EXPR* expr;
+        SCIP_VAR* myVar = var;
+        SCIPcreateExprVar(scip, &expr, myVar, NULL, NULL);
+        lhsExprs[i] = expr;
+        i++;
+    }
+    assert (i == nLhsExprs - 1);
+    // auxilary variables
+    SCIP_VAR* t = NULL;
+    SCIP_VAR* z = NULL;
+    SCIPcreateVarBasic(scip, &t, "t", 0, SCIPinfinity(scip), 1, SCIP_Vartype::SCIP_VARTYPE_CONTINUOUS);
+    SCIPcreateVarBasic(scip, &z, "z", 0, SCIPinfinity(scip), 1, SCIP_Vartype::SCIP_VARTYPE_CONTINUOUS);
+    SCIPaddVar(scip, t);
+    SCIPaddVar(scip, z);
+    SCIP_EXPR* tExpr;
+    SCIP_EXPR* zExpr;
+    SCIPcreateExprVar(scip, &tExpr, t, NULL, NULL);
+    SCIPcreateExprVar(scip, &zExpr, z, NULL, NULL);
+
+    // create expression for last item in vector of the LHS
+    SCIP_EXPR* lastExpr;
+    SCIP_EXPR* sum[2];
+    sum[0] = tExpr;
+    sum[1] = zExpr;
+    std::vector<double> lastCoefs = {1, 1};
+    SCIPcreateExprSum(scip, &lastExpr, 2, sum, lastCoefs.data(), - 2 * alpha, NULL, NULL);
+    lhsCoefs[nLhsExprs - 1] = 1;
+    lhsExprs[nLhsExprs - 1] = lastExpr;
+
+    // set RHS: z + 2a - t
+    std::vector<double> rhsCoefs = {1, -1};
+    std::vector<SCIP_VAR*> rhsVars = {z, t};
+    double rhsOffset = 2 * alpha;   // constant on RHS
+
+    // create the second-order cone constraint
+    SCIPcreateConsSOC(
+        scip,
+        &robustCons,
+        "robust-variance-soc",
+        nLhsExprs,
+        lhsExprs.data(),
+        lhsCoefs.data(),
+        NULL,
+        0,
+        rhsVars.size(),
+        rhsVars.data(),
+        rhsCoefs.data(),
+        rhsOffset
+    );
+    // add the second-order cone constraint
+    SCIPaddCons(scip, robustCons);
 }
 
-std::vector<std::pair<PCTSPvertex, PCTSPvertex>> DistRobustPrizeCollectingTsp(
+std::vector<std::pair<PCTSPvertex, PCTSPvertex>> solveDistRobustPrizeCollectingTsp(
     SCIP * scip,
     PCTSPgraph& graph,
-    EdgeCostMap& cost_mean_map,
-    EdgeCostMap& cost_var_map,
-    VertexPrizeMap& prize_map,
+    EdgeCostMap& costMeanMap,
+    EdgeCostMap& costSigmaMap,
+    VertexPrizeMap& prizeMap,
     PrizeNumberType& quota,
-    PCTSPvertex& root_vertex,
+    PCTSPvertex& rootVertex,
     std::string& name
 ) {
     // add variables, constraints and the SEC cutting plane
-    std::vector<PCTSPedge> heuristic_edges = {};
-    auto edge_var_map = modelPrizeCollectingTSP(
-        scip, graph, heuristic_edges, cost_mean_map, prize_map, quota, root_vertex, name);
+    std::vector<PCTSPedge> heuristicEdges = {};
+    PCTSPedgeVariableMap edgeVarMap = modelPrizeCollectingTSP(
+        scip, graph, heuristicEdges, costMeanMap, prizeMap, quota, rootVertex, name);
+
+    // add distributionally robust constraints
+    addDistRobustCons(scip, graph, costSigmaMap, edgeVarMap);
+
+    // solve the model
+    SCIPsolve(scip);
+
+    // get the solution
+    std::vector<PCTSPedge> solutionEdges = std::vector<PCTSPedge>();
+    if (SCIPgetNSols(scip) > 0) {
+        SCIP_SOL* sol = SCIPgetBestSol(scip);
+        solutionEdges = getSolutionEdges(scip, graph, sol, edgeVarMap);
+    }
+    // return solution
+    return getVertexPairVectorFromEdgeSubset(graph, solutionEdges);
 }
 
 /** creates and captures a nonlinear constraint that is a second-order cone constraint with all its constraint flags set to their default values
@@ -35,9 +121,9 @@ SCIP_RETCODE SCIPcreateConsSOC(
     SCIP*                 scip,               /**< SCIP data structure */
     SCIP_CONS**           cons,               /**< pointer to hold the created constraint */
     const char*           name,               /**< name of constraint */
-    int                   nvars,              /**< number of variables on left hand side of constraint (n) */
-    SCIP_VAR**            vars,               /**< array with variables on left hand side (x_i) */
-    SCIP_Real*            coefs,              /**< array with coefficients of left hand side variables (alpha_i), or NULL if all 1.0 */
+    int                   nlhsexprs,          /**< number of variables on left hand side of constraint (n) */
+    SCIP_EXPR**           lhsexprs,           /**< array of expressions on left hand side (x_i) */
+    SCIP_Real*            lhscoefs,           /**< array with coefficients of left hand side variables (alpha_i), or NULL if all 1.0 */
     SCIP_Real*            offsets,            /**< array with offsets of variables (beta_i), or NULL if all 0.0 */
     SCIP_Real             constant,           /**< constant on left hand side (gamma) */
     int                   nrhsvars,           /**< number of variables on right hand side of constraint */
@@ -52,29 +138,29 @@ SCIP_RETCODE SCIPcreateConsSOC(
     SCIP_Real termcoefs[2];
     int i;
 
-    assert(vars != NULL || nvars == 0);
+    assert(lhsexprs != NULL || nlhsexprs == 0);
 
     SCIP_CALL( SCIPcreateExprSum(scip, &lhssum, 0, NULL, NULL, constant, NULL, NULL) );  /* gamma */
-    for( i = 0; i < nvars; ++i )
+    for( i = 0; i < nlhsexprs; i++ )
     {
-        SCIP_EXPR* varexpr;
         SCIP_EXPR* powexpr;
+        SCIP_EXPR* term;
 
-        SCIP_CALL( SCIPcreateExprVar(scip, &varexpr, vars[i], NULL, NULL) );   /* x_i */
         if( offsets != NULL && offsets[i] != 0.0 )
         {
             SCIP_EXPR* sum;
-            SCIP_CALL( SCIPcreateExprSum(scip, &sum, 1, &varexpr, NULL, offsets[i], NULL, NULL) );  /* x_i + beta_i */
+            SCIP_CALL( SCIPcreateExprSum(scip, &sum, 1, &lhsexprs[i], NULL, offsets[i], NULL, NULL) );  /* x_i + beta_i */
             SCIP_CALL( SCIPcreateExprPow(scip, &powexpr, sum, 2.0, NULL, NULL) );   /* (x_i + beta_i)^2 */
             SCIP_CALL( SCIPreleaseExpr(scip, &sum) );
         }
         else
         {
-            SCIP_CALL( SCIPcreateExprPow(scip, &powexpr, varexpr, 2.0, NULL, NULL) );  /* x_i^2 */
+            // FIXME looks like too many edges, errors at i=20
+            term = lhsexprs[i];
+            SCIP_CALL( SCIPcreateExprPow(scip, &powexpr, term, 2.0, NULL, NULL) );  /* x_i^2 */
         }
-
-        SCIP_CALL( SCIPappendExprSumExpr(scip, lhssum, powexpr, coefs != NULL ? coefs[i]*coefs[i] : 1.0) );  /* + alpha_i^2 (x_i + beta_i)^2 */
-        SCIP_CALL( SCIPreleaseExpr(scip, &varexpr) );
+        SCIP_CALL(SCIPreleaseExpr(scip, &term));
+        SCIP_CALL( SCIPappendExprSumExpr(scip, lhssum, powexpr, lhscoefs != NULL ? lhscoefs[i]*lhscoefs[i] : 1.0) );  /* + alpha_i^2 (x_i + beta_i)^2 */
         SCIP_CALL( SCIPreleaseExpr(scip, &powexpr) );
     }
 

tests/fixtures.cpp

@@ -72,6 +72,19 @@ int GraphFixture::getNumVertices() {
     }
 }
 
+EdgeCostMap GraphFixture::getCostSigmaMap(PCTSPgraph& graph) {
+    EdgeCostMap cost_map = boost::get(edge_weight, graph);
+    for (auto edge : boost::make_iterator_range(boost::edges(graph))) {
+        auto source = boost::source(edge, graph);
+        auto target = boost::target(edge, graph);
+        if (((source == 1) & (target == 4)) | ((source == 3) & (target == 5)))
+            cost_map[edge] = 5;     // make two heavy edges in the middle
+        else
+            cost_map[edge] = 1;
+    }
+    return cost_map;
+}
+
 EdgeCostMap GraphFixture::getCostMap(PCTSPgraph& graph) {
     EdgeCostMap cost_map = boost::get(edge_weight, graph);
     switch (GetParam()) {

tests/fixtures.hh

@@ -19,6 +19,7 @@ class GraphFixture : public::testing::TestWithParam<GraphType> {
 	public:
 		PCTSPgraph getGraph();
 		EdgeCostMap getCostMap(PCTSPgraph& graph);
+		EdgeCostMap getCostSigmaMap(PCTSPgraph& graph);
 		VertexPrizeMap getPrizeMap(PCTSPgraph& graph);
 		VertexPrizeMap getGenOnePrizeMap(PCTSPgraph& graph);
 		std::vector<std::pair<int, int>> getEdgeVector();

tests/test_robust.cpp

@@ -0,0 +1,48 @@
+#include <gtest/gtest.h>
+#include <objscip/objscipdefplugins.h>
+#include "pctsp/graph.hh"
+#include "pctsp/robust.hh"
+#include "fixtures.hh"
+
+
+typedef GraphFixture RobustFixture;
+
+TEST_P(RobustFixture, testSolveRobustDistRobustPrizeCollectingTsp) {
+    auto graph = getGraph();
+    auto quota = getQuota();
+    auto rootVertex = getRootVertex();
+    auto costMeanMap = getCostMap(graph);
+    auto costSigmaMap = getCostSigmaMap(graph);
+    auto prizeMap = getPrizeMap(graph);
+    // addSelfLoopsToGraph(graph);
+    // assignZeroCostToSelfLoops(graph, costMeanMap);
+    // assignZeroCostToSelfLoops(graph, costSigmaMap);
+
+    std::string testName = "testSolveRobustDistRobustPrizeCollectingTsp-";
+    testName.append(getParamName());
+
+    std::vector<PCTSPvertex> tour;
+    auto first = tour.begin();
+    auto last = tour.end();
+    std::vector<PCTSPedge> heuristicEdges = getEdgesInWalk(graph, first, last);
+
+    SCIP* scip = NULL;
+    SCIPcreate(&scip);
+    SCIPincludeDefaultPlugins(scip);    // FIXME!!!!
+    auto solution = solveDistRobustPrizeCollectingTsp(
+        scip,
+        graph,
+        costMeanMap,
+        costSigmaMap,
+        prizeMap,
+        quota,
+        rootVertex,
+        testName
+    );
+}
+
+INSTANTIATE_TEST_SUITE_P(
+    TestAlgorithms,
+    RobustFixture,
+    ::testing::Values(GraphType::SUURBALLE, GraphType::GRID8, GraphType::COMPLETE4, GraphType::COMPLETE5)
+);

tests/test_scipsdp.cpp

@@ -14,9 +14,11 @@ TEST(TestSoc, testSoc) {
 
     int nVars = 4;
     VarVector vars (nVars);
-    for (auto& var : vars) {
-        SCIPcreateVarBasic(scip, &var, NULL, 0, 1, 0, SCIP_Vartype::SCIP_VARTYPE_BINARY);
-        SCIPaddVar(scip, var);
+    SCIP_EXPR* lhsExprs[nVars];
+    for (int i = 0; i < nVars; i++) {
+        SCIPcreateVarBasic(scip, &vars[i], NULL, 0, 1, 0, SCIP_Vartype::SCIP_VARTYPE_BINARY);
+        SCIPaddVar(scip, vars[i]);
+        SCIPcreateExprVar(scip, &lhsExprs[i], vars[i], NULL, NULL);
     }
 
     SCIP_VAR* t = NULL;
@@ -32,7 +34,7 @@ TEST(TestSoc, testSoc) {
     // SCIPcreateConsBasicSOCNonlinear(scip, &socCons, consName.c_str(), nVars, vars.data(), NULL, NULL, 0, t, 1, 0);
     std::vector<SCIP_VAR*> rhsVars = {t, z};
     std::vector<double> rhsCoefs = {1, 1};
-    SCIPcreateConsSOC(scip, &socCons, consName.c_str(), nVars, vars.data(), NULL, NULL, 0, rhsVars.size(), rhsVars.data(), rhsCoefs.data(), 0);
+    SCIPcreateConsSOC(scip, &socCons, consName.c_str(), nVars, lhsExprs, NULL, NULL, 0, rhsVars.size(), rhsVars.data(), rhsCoefs.data(), 0);
     SCIPaddCons(scip, socCons);
 
     SCIP_CONS* knapsack = NULL;