Make CbcPolyhedronSolver more robust

5 files changed, 134 insertions, 84 deletions

diff --git a/Solvers/ILP-Solver/hu.bme.mit.inf.dslreasoner.ilp.cbc/cpp/viatracbc.cpp b/Solvers/ILP-Solver/hu.bme.mit.inf.dslreasoner.ilp.cbc/cpp/viatracbc.cpp
index 49994244..ffd35759 100644
--- a/Solvers/ILP-Solver/hu.bme.mit.inf.dslreasoner.ilp.cbc/cpp/viatracbc.cpp
+++ b/Solvers/ILP-Solver/hu.bme.mit.inf.dslreasoner.ilp.cbc/cpp/viatracbc.cpp
@@ -36,9 +36,10 @@ static const jint kCbcError = 5;
 static CoinModel CreateModel(JNIEnv *env, jdoubleArray columnLowerBoundsArray,
     jdoubleArray columnUpperBoundsArray, jintArray rowStartsArray, jintArray columnIndicesArray,
     jdoubleArray entriedArray, jdoubleArray rowLowerBoundsArray, jdoubleArray rowUpperBoundsArray,
-    jdoubleArray objectiveArray);
+    jdoubleArray objectiveArray, jboolean lpRelaxation);
 static void CreateModelColumns(JNIEnv *env, jdoubleArray columnLowerBoundsArray,
-    jdoubleArray columnUpperBoundsArray, jdoubleArray objectiveArray, CoinModel &build);
+    jdoubleArray columnUpperBoundsArray, jdoubleArray objectiveArray, jboolean lpRelaxation,
+    CoinModel &build);
 static void CreateModelRows(JNIEnv *env, jintArray rowStartsArray, jintArray columnIndicesArray,
     jdoubleArray entriesArray, jdoubleArray rowLowerBoundsArray, jdoubleArray rowUpperBoundsArray,
     CoinModel &build);
@@ -83,11 +84,11 @@ jint Java_hu_bme_mit_inf_dslreasoner_ilp_cbc_CbcSolver_solveIlpProblem(
     JNIEnv *env, jclass klazz, jdoubleArray columnLowerBoundsArray, jdoubleArray columnUpperBoundsArray,
     jintArray rowStartsArray, jintArray columnIndicesArray, jdoubleArray entriesArray,
     jdoubleArray rowLowerBoundsArray, jdoubleArray rowUpperBoundsArray, jdoubleArray objectiveArray,
-    jdoubleArray outputArray, jdouble timeoutSeconds, jboolean silent) {
+    jdoubleArray outputArray, jboolean lpRelaxation, jdouble timeoutSeconds, jboolean silent) {
     try {
         auto build = CreateModel(env, columnLowerBoundsArray, columnUpperBoundsArray,
             rowStartsArray, columnIndicesArray, entriesArray, rowLowerBoundsArray, rowUpperBoundsArray,
-            objectiveArray);
+            objectiveArray, lpRelaxation);
         double value;
         jint result = SolveModel(build, timeoutSeconds, silent, value);
         if (result == kCbcSolutionBounded) {
@@ -106,16 +107,18 @@ jint Java_hu_bme_mit_inf_dslreasoner_ilp_cbc_CbcSolver_solveIlpProblem(
 CoinModel CreateModel(JNIEnv *env, jdoubleArray columnLowerBoundsArray,
     jdoubleArray columnUpperBoundsArray, jintArray rowStartsArray, jintArray columnIndicesArray,
     jdoubleArray entriesArray, jdoubleArray rowLowerBoundsArray, jdoubleArray rowUpperBoundsArray,
-    jdoubleArray objectiveArray) {
+    jdoubleArray objectiveArray, jboolean lpRelaxation) {
     CoinModel build;
-    CreateModelColumns(env, columnLowerBoundsArray, columnUpperBoundsArray, objectiveArray, build);
+    CreateModelColumns(env, columnLowerBoundsArray, columnUpperBoundsArray, objectiveArray,
+        lpRelaxation, build);
     CreateModelRows(env, rowStartsArray, columnIndicesArray, entriesArray, rowLowerBoundsArray,
         rowUpperBoundsArray, build);
     return build;
 }
 
 void CreateModelColumns(JNIEnv *env, jdoubleArray columnLowerBoundsArray,
-    jdoubleArray columnUpperBoundsArray, jdoubleArray objectiveArray, CoinModel &build) {
+    jdoubleArray columnUpperBoundsArray, jdoubleArray objectiveArray, jboolean lpRelaxation,
+    CoinModel &build) {
     int numColumns = env->GetArrayLength(columnLowerBoundsArray);
     PinnedDoubleArray columnLowerBounds{env, columnLowerBoundsArray};
     PinnedDoubleArray columnUpperBounds{env, columnUpperBoundsArray};
@@ -123,7 +126,9 @@ void CreateModelColumns(JNIEnv *env, jdoubleArray columnLowerBoundsArray,
     for (int i = 0; i < numColumns; i++) {
         build.setColumnBounds(i, columnLowerBounds[i], columnUpperBounds[i]);
         build.setObjective(i, objective[i]);
-        build.setInteger(i);
+        if (!lpRelaxation) {
+            build.setInteger(i);
+        }
     }
 }
 
@@ -215,6 +220,9 @@ jint SolveModel(CoinModel &build, jdouble timeoutSeconds, jboolean silent, jdoub
     if (model.isInitialSolveProvenPrimalInfeasible()) {
         return kCbcUnsat;
     }
+    if (model.isInitialSolveProvenDualInfeasible()) {
+        return kCbcSolutionUnbounded;
+    }
     if (model.isInitialSolveAbandoned()) {
         return kCbcTimeout;
     }
@@ -226,20 +234,26 @@ jint SolveModel(CoinModel &build, jdouble timeoutSeconds, jboolean silent, jdoub
 
     model.branchAndBound();
 
-    if (model.isProvenInfeasible()) {
-        return kCbcUnsat;
-    }
-    if (model.isProvenDualInfeasible()) {
-        return kCbcSolutionUnbounded;
-    }
-    if (model.isProvenOptimal()) {
-        value = model.getMinimizationObjValue();
-        return kCbcSolutionBounded;
-    }
-    if (model.maximumSecondsReached()) {
+    switch (model.status()) {
+    case 0:
+        if (model.isProvenInfeasible()) {
+            return kCbcUnsat;
+        }
+        if (model.isProvenDualInfeasible()) {
+            return kCbcSolutionUnbounded;
+        }
+        if (model.isProvenOptimal()) {
+            value = model.getMinimizationObjValue();
+            return kCbcSolutionBounded;
+        }
+        throw std::runtime_error("CBC status is 0, but no solution is found");
+    case 1:
         return kCbcTimeout;
+    case 2:
+        return kCbcAbandoned;
+    default:
+        throw std::runtime_error("Unknown CBC status");
     }
-    return kCbcAbandoned;
 }
 
 void ThrowException(JNIEnv *env, const char *message) {
diff --git a/Solvers/ILP-Solver/hu.bme.mit.inf.dslreasoner.ilp.cbc/cpp/viatracbc.hpp b/Solvers/ILP-Solver/hu.bme.mit.inf.dslreasoner.ilp.cbc/cpp/viatracbc.hpp
index c65f71e3..12198c8b 100644
--- a/Solvers/ILP-Solver/hu.bme.mit.inf.dslreasoner.ilp.cbc/cpp/viatracbc.hpp
+++ b/Solvers/ILP-Solver/hu.bme.mit.inf.dslreasoner.ilp.cbc/cpp/viatracbc.hpp
@@ -9,7 +9,7 @@ JNIEXPORT jint JNICALL Java_hu_bme_mit_inf_dslreasoner_ilp_cbc_CbcSolver_solveIl
     JNIEnv *env, jclass klazz, jdoubleArray columnLowerBoundsArray, jdoubleArray columnUpperBoundsArray,
     jintArray rowStartsArray, jintArray columnIndicesArray, jdoubleArray entriesArray,
     jdoubleArray rowLowerBoundsArray, jdoubleArray rowUpperBoundsArray, jdoubleArray objectiveArray,
-    jdoubleArray outputArray, jdouble timeoutSeconds, jboolean silent);
+    jdoubleArray outputArray, jboolean lpRelaxation, jdouble timeoutSeconds, jboolean silent);
 
 }
 
diff --git a/Solvers/ILP-Solver/hu.bme.mit.inf.dslreasoner.ilp.cbc/lib/libviatracbc.so b/Solvers/ILP-Solver/hu.bme.mit.inf.dslreasoner.ilp.cbc/lib/libviatracbc.so
index 21fd2ff2..4eae7de6 100755
--- a/Solvers/ILP-Solver/hu.bme.mit.inf.dslreasoner.ilp.cbc/lib/libviatracbc.so
+++ b/Solvers/ILP-Solver/hu.bme.mit.inf.dslreasoner.ilp.cbc/lib/libviatracbc.so
diff --git a/Solvers/ILP-Solver/hu.bme.mit.inf.dslreasoner.ilp.cbc/src/hu/bme/mit/inf/dslreasoner/ilp/cbc/CbcSolver.java b/Solvers/ILP-Solver/hu.bme.mit.inf.dslreasoner.ilp.cbc/src/hu/bme/mit/inf/dslreasoner/ilp/cbc/CbcSolver.java
index 39b9d537..085d4448 100644
--- a/Solvers/ILP-Solver/hu.bme.mit.inf.dslreasoner.ilp.cbc/src/hu/bme/mit/inf/dslreasoner/ilp/cbc/CbcSolver.java
+++ b/Solvers/ILP-Solver/hu.bme.mit.inf.dslreasoner.ilp.cbc/src/hu/bme/mit/inf/dslreasoner/ilp/cbc/CbcSolver.java
@@ -15,14 +15,14 @@ public class CbcSolver {
         }
 
         public static CbcResult solve(double[] columnLowerBounds, double[] columnUpperBounds, int[] rowStarts,
-                        int[] columnIndices, double[] entries, double[] rowLowerBounds, double[] rowUpperBounds,
-                        double[] objective, double timeoutSeconds, boolean silent) {
+                        int[] columnIndices, double[] entries, double[] rowLowerBounds, double[] rowUpperBounds, double[] objective,
+                        boolean lpRelaxation, double timeoutSeconds, boolean silent) {
                 loadNatives();
                 validate(columnLowerBounds, columnUpperBounds, rowStarts, columnIndices, entries, rowLowerBounds,
                                 rowUpperBounds, objective);
                 double[] output = new double[1];
                 int result = solveIlpProblem(columnLowerBounds, columnUpperBounds, rowStarts, columnIndices, entries,
-                                rowLowerBounds, rowUpperBounds, objective, output, timeoutSeconds, silent);
+                                rowLowerBounds, rowUpperBounds, objective, output, lpRelaxation, timeoutSeconds, silent);
                 if (result == CBC_SOLUTION_BOUNDED) {
                         return new CbcResult.SolutionBounded(output[0]);
                 } else if (result == CBC_SOLUTION_UNBOUNDED) {
@@ -67,5 +67,5 @@ public class CbcSolver {
 
         private static native int solveIlpProblem(double[] columnLowerBounds, double[] columnUpperBounds, int[] rowStarts,
                         int[] columnIndices, double[] entries, double[] rowLowerBounds, double[] rowUpperBounds, double[] objective,
-                        double[] output, double timeoutSeconds, boolean silent);
+                        double[] output, boolean lpRelaxation, double timeoutSeconds, boolean silent);
 }
diff --git a/Solvers/VIATRA-Solver/hu.bme.mit.inf.dslreasoner.viatrasolver.logic2viatra/src/hu/bme/mit/inf/dslreasoner/viatrasolver/logic2viatra/cardinality/CbcPolyhedronSolver.xtend b/Solvers/VIATRA-Solver/hu.bme.mit.inf.dslreasoner.viatrasolver.logic2viatra/src/hu/bme/mit/inf/dslreasoner/viatrasolver/logic2viatra/cardinality/CbcPolyhedronSolver.xtend
index 7753e68e..4bd46fbf 100644
--- a/Solvers/VIATRA-Solver/hu.bme.mit.inf.dslreasoner.viatrasolver.logic2viatra/src/hu/bme/mit/inf/dslreasoner/viatrasolver/logic2viatra/cardinality/CbcPolyhedronSolver.xtend
+++ b/Solvers/VIATRA-Solver/hu.bme.mit.inf.dslreasoner.viatrasolver.logic2viatra/src/hu/bme/mit/inf/dslreasoner/viatrasolver/logic2viatra/cardinality/CbcPolyhedronSolver.xtend
@@ -1,27 +1,32 @@
 package hu.bme.mit.inf.dslreasoner.viatrasolver.logic2viatra.cardinality
 
+import com.google.common.collect.ImmutableList
 import com.google.common.collect.ImmutableMap
 import hu.bme.mit.inf.dslreasoner.ilp.cbc.CbcResult
 import hu.bme.mit.inf.dslreasoner.ilp.cbc.CbcSolver
+import java.util.HashSet
 import java.util.List
 import java.util.Map
+import java.util.Set
 import org.eclipse.xtend.lib.annotations.FinalFieldsConstructor
 
 @FinalFieldsConstructor
 class CbcPolyhedronSolver implements PolyhedronSolver {
+        val boolean lpRelaxation
         val double timeoutSeconds
         val boolean silent
 
         new() {
-                this(10, true)
+                this(false, -1, true)
         }
 
         override createSaturationOperator(Polyhedron polyhedron) {
-                new CbcSaturationOperator(polyhedron, timeoutSeconds, silent)
+                new CbcSaturationOperator(polyhedron, lpRelaxation, timeoutSeconds, silent)
         }
 }
 
 class CbcSaturationOperator extends AbstractPolyhedronSaturationOperator {
+        val boolean lpRelaxation
         val double timeoutSeconds
         val boolean silent
         val double[] columnLowerBounds
@@ -29,8 +34,9 @@ class CbcSaturationOperator extends AbstractPolyhedronSaturationOperator {
         val double[] objective
         val Map<Dimension, Integer> dimensionsToIndicesMap
 
-        new(Polyhedron polyhedron, double timeoutSeconds, boolean silent) {
+        new(Polyhedron polyhedron, boolean lpRelaxation, double timeoutSeconds, boolean silent) {
                 super(polyhedron)
+                this.lpRelaxation = lpRelaxation
                 this.timeoutSeconds = timeoutSeconds
                 this.silent = silent
                 val numDimensions = polyhedron.dimensions.size
@@ -56,6 +62,12 @@ class CbcSaturationOperator extends AbstractPolyhedronSaturationOperator {
                 rowStarts.set(numConstraints, numEntries)
                 val columnIndices = newIntArrayOfSize(numEntries)
                 val entries = newDoubleArrayOfSize(numEntries)
+                val unconstrainedDimensions = new HashSet
+                for (dimension : polyhedron.dimensions) {
+                        if (dimension.lowerBound === null && dimension.upperBound === null) {
+                                unconstrainedDimensions += dimension
+                        }
+                }
                 var int index = 0
                 for (var int i = 0; i < numConstraints; i++) {
                         rowStarts.set(i, index)
@@ -69,6 +81,7 @@ class CbcSaturationOperator extends AbstractPolyhedronSaturationOperator {
                                 val dimension = polyhedron.dimensions.get(j)
                                 val coefficient = constraint.coefficients.get(dimension)
                                 if (coefficient !== null && coefficient != 0) {
+                                        unconstrainedDimensions -= dimension
                                         columnIndices.set(index, j)
                                         entries.set(index, coefficient)
                                         index++
@@ -79,71 +92,94 @@ class CbcSaturationOperator extends AbstractPolyhedronSaturationOperator {
                         throw new AssertionError("Last entry does not equal the number of entries in the constraint matrix")
                 }
                 for (expressionToSaturate : polyhedron.expressionsToSaturate) {
-                        for (var int j = 0; j < numDimensions; j++) {
-                                objective.set(j, 0)
+                        val result = saturate(expressionToSaturate, rowStarts, columnIndices, entries, rowLowerBounds,
+                                rowUpperBounds, unconstrainedDimensions, constraints)
+                        if (result != PolyhedronSaturationResult.SATURATED) {
+                                return result
                         }
-                        switch (expressionToSaturate) {
-                                Dimension: {
-                                        val j = getIndex(expressionToSaturate)
-                                        objective.set(j, 1)
+                }
+                PolyhedronSaturationResult.SATURATED
+        }
+
+        protected def saturate(LinearBoundedExpression expressionToSaturate, int[] rowStarts, int[] columnIndices,
+                double[] entries, double[] rowLowerBounds, double[] rowUpperBounds, Set<Dimension> unconstrainedDimensions,
+                ImmutableList<LinearConstraint> constraints) {
+                val numDimensions = objective.size
+                for (var int j = 0; j < numDimensions; j++) {
+                        objective.set(j, 0)
+                }
+                switch (expressionToSaturate) {
+                        Dimension: {
+                                // CBC will return nonsensical results or call free() with invalid arguments if
+                                // it is passed a fully unconstrained (-Inf lower and +Int upper bound, no inequalities) variable
+                                // in the objective function.
+                                if (unconstrainedDimensions.contains(expressionToSaturate)) {
+                                        return PolyhedronSaturationResult.SATURATED
                                 }
-                                LinearConstraint: {
-                                        for (pair : expressionToSaturate.coefficients.entrySet) {
-                                                val j = getIndex(pair.key)
-                                                objective.set(j, pair.value)
+                                val j = getIndex(expressionToSaturate)
+                                objective.set(j, 1)
+                        }
+                        LinearConstraint: {
+                                for (pair : expressionToSaturate.coefficients.entrySet) {
+                                        val dimension = pair.key
+                                        // We also have to check for unconstrained dimensions here to avoid crashing.
+                                        if (unconstrainedDimensions.contains(dimension)) {
+                                                expressionToSaturate.lowerBound = null
+                                                expressionToSaturate.upperBound = null
+                                                return PolyhedronSaturationResult.SATURATED
                                         }
+                                        val j = getIndex(dimension)
+                                        objective.set(j, pair.value)
                                 }
-                                default:
-                                        throw new IllegalArgumentException("Unknown expression: " + expressionToSaturate)
                         }
-                        val minimizationResult = CbcSolver.solve(columnLowerBounds, columnUpperBounds, rowStarts, columnIndices,
-                                entries, rowLowerBounds, rowUpperBounds, objective, timeoutSeconds, silent)
-                        switch (minimizationResult) {
-                                CbcResult.SolutionBounded: {
-                                        val value = Math.floor(minimizationResult.value)
-                                        expressionToSaturate.lowerBound = value as int
-                                        setBound(expressionToSaturate, constraints, value, columnLowerBounds, rowLowerBounds)
-                                }
-                                case CbcResult.SOLUTION_UNBOUNDED: {
-                                        expressionToSaturate.lowerBound = null
-                                        setBound(expressionToSaturate, constraints, Double.NEGATIVE_INFINITY, columnLowerBounds,
-                                                rowLowerBounds)
-                                }
-                                case CbcResult.UNSAT:
-                                        return PolyhedronSaturationResult.EMPTY
-                                case CbcResult.ABANDONED,
-                                case CbcResult.TIMEOUT:
-                                        return PolyhedronSaturationResult.UNKNOWN
-                                default:
-                                        throw new RuntimeException("Unknown CbcResult: " + minimizationResult)
+                        default:
+                                throw new IllegalArgumentException("Unknown expression: " + expressionToSaturate)
+                }
+                val minimizationResult = CbcSolver.solve(columnLowerBounds, columnUpperBounds, rowStarts, columnIndices,
+                        entries, rowLowerBounds, rowUpperBounds, objective, lpRelaxation, timeoutSeconds, silent)
+                switch (minimizationResult) {
+                        CbcResult.SolutionBounded: {
+                                val value = Math.floor(minimizationResult.value)
+                                expressionToSaturate.lowerBound = value as int
+                                setBound(expressionToSaturate, constraints, value, columnLowerBounds, rowLowerBounds)
                         }
-                        for (var int j = 0; j < numDimensions; j++) {
-                                val objectiveCoefficient = objective.get(j)
-                                objective.set(j, -objectiveCoefficient)
+                        case CbcResult.SOLUTION_UNBOUNDED: {
+                                expressionToSaturate.lowerBound = null
+                                setBound(expressionToSaturate, constraints, Double.NEGATIVE_INFINITY, columnLowerBounds, rowLowerBounds)
                         }
-                        val maximizationResult = CbcSolver.solve(columnLowerBounds, columnUpperBounds, rowStarts, columnIndices,
-                                entries, rowLowerBounds, rowUpperBounds, objective, timeoutSeconds, silent)
-                        switch (maximizationResult) {
-                                CbcResult.SolutionBounded: {
-                                        val value = Math.ceil(-maximizationResult.value)
-                                        expressionToSaturate.upperBound = value as int
-                                        setBound(expressionToSaturate, constraints, value, columnUpperBounds, rowUpperBounds)
-                                }
-                                case CbcResult.SOLUTION_UNBOUNDED: {
-                                        expressionToSaturate.upperBound = null
-                                        setBound(expressionToSaturate, constraints, Double.POSITIVE_INFINITY, columnUpperBounds,
-                                                rowUpperBounds)
-                                }
-                                case CbcResult.UNSAT:
-                                        throw new RuntimeException("Minimization was SAT, but maximization is UNSAT")
-                                case CbcResult.ABANDONED,
-                                case CbcResult.TIMEOUT:
-                                        return PolyhedronSaturationResult.UNKNOWN
-                                default:
-                                        throw new RuntimeException("Unknown CbcResult: " + maximizationResult)
+                        case CbcResult.UNSAT:
+                                return PolyhedronSaturationResult.EMPTY
+                        case CbcResult.ABANDONED,
+                        case CbcResult.TIMEOUT:
+                                return PolyhedronSaturationResult.UNKNOWN
+                        default:
+                                throw new RuntimeException("Unknown CbcResult: " + minimizationResult)
+                }
+                for (var int j = 0; j < numDimensions; j++) {
+                        val objectiveCoefficient = objective.get(j)
+                        objective.set(j, -objectiveCoefficient)
+                }
+                val maximizationResult = CbcSolver.solve(columnLowerBounds, columnUpperBounds, rowStarts, columnIndices,
+                        entries, rowLowerBounds, rowUpperBounds, objective, lpRelaxation, timeoutSeconds, silent)
+                switch (maximizationResult) {
+                        CbcResult.SolutionBounded: {
+                                val value = Math.ceil(-maximizationResult.value)
+                                expressionToSaturate.upperBound = value as int
+                                setBound(expressionToSaturate, constraints, value, columnUpperBounds, rowUpperBounds)
                         }
+                        case CbcResult.SOLUTION_UNBOUNDED: {
+                                expressionToSaturate.upperBound = null
+                                setBound(expressionToSaturate, constraints, Double.POSITIVE_INFINITY, columnUpperBounds, rowUpperBounds)
+                        }
+                        case CbcResult.UNSAT:
+                                throw new RuntimeException("Minimization was SAT, but maximization is UNSAT")
+                        case CbcResult.ABANDONED,
+                        case CbcResult.TIMEOUT:
+                                return PolyhedronSaturationResult.UNKNOWN
+                        default:
+                                throw new RuntimeException("Unknown CbcResult: " + maximizationResult)
                 }
-                PolyhedronSaturationResult.SATURATED
+                return PolyhedronSaturationResult.SATURATED
         }
 
         private def toDouble(Integer nullableInt, double defaultValue) {