Merge branch 'master' into improve-connected-components

algorithm-engineering-common/src/main/java/berlin/tu/algorithmengineering/common/Graph.java

@@ -9,6 +9,8 @@ import java.util.*;
 
 public class Graph {
 
+    public static final int FORBIDDEN_VALUE = (int) -Math.pow(2, 16);
+
     private int numberOfVertices;
     private int[][] edgeWeights;
     private boolean[][] edgeExists;
@@ -128,6 +130,25 @@ public class Graph {
             return false;
     }
 
+    public int flipEdgeAndSetForbidden(int i, int j) {
+        flipEdge(i, j);
+        int costToFlip = edgeWeights[i][j];
+        edgeWeights[i][j] = FORBIDDEN_VALUE;
+        edgeWeights[j][i] = FORBIDDEN_VALUE;
+        absoluteNeighborhoodWeights[i] += -FORBIDDEN_VALUE + costToFlip;
+        absoluteNeighborhoodWeights[j] += -FORBIDDEN_VALUE + costToFlip;
+
+        return costToFlip;
+    }
+
+    public void flipBackForbiddenEdge(int i, int j, int cost) {
+        edgeWeights[i][j] = cost;
+        edgeWeights[j][i] = cost;
+        absoluteNeighborhoodWeights[i] -= -FORBIDDEN_VALUE + cost;
+        absoluteNeighborhoodWeights[j] -= -FORBIDDEN_VALUE + cost;
+        flipEdge(i, j);
+    }
+
     public MergeVerticesInfo mergeVertices(int a, int b) {
         if (a > b) {
             int temp = a;
@@ -539,6 +560,175 @@ public class Graph {
         return lowerBound;
     }
 
+    public int getLowerBound2BasedOnMaximumWeightIndependentSetMoreEfficient(List<P3> p3List) {
+        int[] vertexWeights = new int[p3List.size()];
+        boolean[][] edges = new boolean[vertexWeights.length][vertexWeights.length];
+
+        for (int i = 0; i < vertexWeights.length; i++) {
+            vertexWeights[i] = getSmallestAbsoluteWeight(p3List.get(i));
+        }
+
+        for (int i = 0; i < vertexWeights.length; i++) {
+            for (int j = i+1; j < vertexWeights.length; j++) {
+                boolean edgeExists = getNumberOfSharedVertices(p3List.get(i), p3List.get(j)) > 1;
+                edges[i][j] = edgeExists;
+                edges[j][i] = edgeExists;
+            }
+        }
+
+        boolean[] independentSet = new boolean[vertexWeights.length];
+        boolean[] freeVertices = new boolean[vertexWeights.length];
+        int numberOfFreeVertices = vertexWeights.length;
+        Arrays.fill(freeVertices, true);
+
+        createInitialSelection(edges, independentSet, freeVertices);
+
+        return performSimulatedAnnealingOnMaximumWeightIndependent(vertexWeights, edges, independentSet, freeVertices);
+    }
+
+    private static int performSimulatedAnnealingOnMaximumWeightIndependent(
+            int[] vertexWeights, boolean[][] edges, boolean[] independentSet, boolean[] freeVertices
+    ) {
+        final double T_START = 2;
+        double t = T_START;
+        final int ITERATIONS = 10_000;
+
+        int currentSolution = getIndependentSetWeight(vertexWeights, independentSet);
+        int highestSolution = currentSolution;
+
+        for (int i = 0; i < ITERATIONS; i++) {
+            boolean[] perturbedIndependentSet = Arrays.copyOf(independentSet, independentSet.length);
+            boolean[] perturbedFreeVertices = Arrays.copyOf(freeVertices, freeVertices.length);
+            int deltaSolution = perturbIndependentSet(vertexWeights, edges, perturbedIndependentSet, perturbedFreeVertices);
+            if (deltaSolution >= 0 || Utils.RANDOM.nextDouble() < Math.exp(deltaSolution / t)) {
+                currentSolution += deltaSolution;
+                independentSet = perturbedIndependentSet;
+                freeVertices = perturbedFreeVertices;
+                if (currentSolution > highestSolution) {
+                    highestSolution = currentSolution;
+                }
+            }
+            t -= T_START / ITERATIONS;
+        }
+
+        return highestSolution;
+    }
+
+    private static int perturbIndependentSet(int[] vertexWeights, boolean[][] edges, boolean[] independentSet, boolean[] freeVertices) {
+        int deltaSolution = 0;
+
+        int sizeOfIndependentSet = getNumberOfTrueValues(independentSet);
+        int index = Utils.randInt(0, sizeOfIndependentSet);
+
+        for (int i = 0; i < independentSet.length; i++) {
+            if (independentSet[i]) {
+                if (index == 0) {
+                    independentSet[i] = false;
+                    freeVertices[i] = true;
+                    deltaSolution -= vertexWeights[i];
+                    for (int j = 0; j < edges.length; j++) {
+                        if (i != j && edges[i][j] && isFreeVertex(j, edges, independentSet)) {
+                            freeVertices[j] = true;
+                        }
+                    }
+                    break;
+                } else {
+                    index--;
+                }
+            }
+        }
+
+        int numberOfFreeVertices = getNumberOfTrueValues(freeVertices);
+
+        while (numberOfFreeVertices > 0) {
+            index = Utils.randInt(0, numberOfFreeVertices);
+
+            for (int i = 0; i < freeVertices.length; i++) {
+                if (freeVertices[i]) {
+                    if (index == 0) {
+                        freeVertices[i] = false;
+                        numberOfFreeVertices--;
+                        independentSet[i] = true;
+                        deltaSolution += vertexWeights[i];
+                        for (int j = 0; j < independentSet.length; j++) {
+                            if (i != j && freeVertices[j] && edges[i][j]) {
+                                freeVertices[j] = false;
+                                numberOfFreeVertices--;
+                            }
+                        }
+                        break;
+                    }
+                    index--;
+                }
+            }
+        }
+
+        return deltaSolution;
+    }
+
+    private static boolean isFreeVertex(int vertex, boolean[][] edges, boolean[] independentSet) {
+        for (int i = 0; i < independentSet.length; i++) {
+            if (i != vertex && edges[vertex][i] && independentSet[i]) {
+                return false;
+            }
+        }
+        return true;
+    }
+
+    private static int getNumberOfTrueValues(boolean[] independentSet) {
+        int size = 0;
+        for (int i = 0; i < independentSet.length; i++) {
+            if (independentSet[i]) {
+                size++;
+            }
+        }
+        return size;
+    }
+
+    private static int getIndependentSetWeight(int[] vertexWeights, boolean[] independentSet) {
+        int weight = 0;
+        for (int i = 0; i < independentSet.length; i++) {
+            if (independentSet[i]) {
+                weight += vertexWeights[i];
+            }
+        }
+        return weight;
+    }
+
+    private static void createInitialSelection(boolean[][] edges, boolean[] independentSet, boolean[] freeVertices) {
+        int numberOfFreeVertices = freeVertices.length;
+
+        while (numberOfFreeVertices > 0) {
+            for (int i = 0; i < freeVertices.length; i++) {
+                if (freeVertices[i]) {
+                    independentSet[i] = true;
+                    freeVertices[i] = false;
+                    numberOfFreeVertices--;
+                    for (int j = 0; j < edges.length; j++) {
+                        if (edges[i][j]) {
+                            freeVertices[j] = false;
+                            numberOfFreeVertices--;
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    public int getNumberOfSharedVertices(P3 a, P3 b) {
+        int sharedVertices = 0;
+        if (a.getU() == b.getU() || a.getU() == b.getV() || a.getU() == b.getW()) {
+            sharedVertices++;
+        }
+        if (a.getV() == b.getU() || a.getV() == b.getV() || a.getV() == b.getW()) {
+            sharedVertices++;
+        }
+        if (a.getW() == b.getU() || a.getW() == b.getV() || a.getW() == b.getW()) {
+            sharedVertices++;
+        }
+        return sharedVertices;
+    }
+
     public int getLowerBound2BasedOnMaximumWeightIndependentSet(List<P3> p3List) {
         Map<P3, P3WithSharedEdgeP3s> p3WithSharedEdgeP3sMap = P3WithSharedEdgeP3s.fromP3List(this, p3List);
         Set<P3> freeP3s = new HashSet<>(p3List);
@@ -600,7 +790,6 @@ public class Graph {
     /**
      *
      * @param p3WithSharedEdgeP3sMap
-     * @param freeP3s
      * @param selectedP3s
      * @return deltaSolution
      */

algorithm-engineering-heuristics/src/main/java/berlin/tu/algorithmengineering/heuristics/SimulatedAnnealing.java

@@ -108,6 +108,10 @@ public class SimulatedAnnealing {
     }
 
     private static int getDeltaCost(Graph graph, boolean[][] resultEdgeExists, int vertex, int moveToVertex) {
+        if (resultEdgeExists[vertex][moveToVertex]) {
+            return 0;
+        }
+
         int deltaCost = 0;
 
         for (int i = 0; i < graph.getNumberOfVertices(); i++) {

algorithm-engineering-search-tree/src/main/java/berlin/tu/algorithmengineering/searchtree/Main.java

@@ -11,6 +11,7 @@ import berlin.tu.algorithmengineering.common.model.ResultEdgeExistsWithSolutionS
 import berlin.tu.algorithmengineering.heuristics.Heuristics;
 import berlin.tu.algorithmengineering.heuristics.SimulatedAnnealing;
 import berlin.tu.algorithmengineering.searchtree.lp.LpLowerBound;
+import gurobi.*;
 
 import java.util.*;
 
@@ -37,7 +38,7 @@ public class Main {
         System.out.printf("#output cost: %d\n", Utils.getCostToChange(graph, resultEdgeExistsWithSolutionSize.getResultEdgeExists()));
 
         System.out.printf("#recursive steps: %d\n", recursiveSteps);
-//        System.out.printf("#recursive steps: %d\n", graph.getLowerBound2(graph.findAllP3()));
+//        System.out.printf("#recursive steps: %d\n", getLowerBound2BasedOnMaximumWeightIndependentSetIlp(graph, graph.findAllP3()));
     }
 
     private static boolean[][] weightedClusterEditingBranch(Graph graph, int k) {
@@ -113,24 +114,67 @@ public class Main {
             return resultEdgeExists;
         }
 
-        graph.flipEdge(p3.getU(), p3.getV());
-        int cost = graph.getEdgeWeights()[p3.getU()][p3.getV()];
-        graph.getEdgeWeights()[p3.getU()][p3.getV()] = FORBIDDEN_VALUE;
-        graph.getEdgeWeights()[p3.getV()][p3.getU()] = FORBIDDEN_VALUE;
-        graph.getAbsoluteNeighborhoodWeights()[p3.getU()] += -FORBIDDEN_VALUE + cost;
-        graph.getAbsoluteNeighborhoodWeights()[p3.getV()] += -FORBIDDEN_VALUE + cost;
+        int cost = graph.flipEdgeAndSetForbidden(p3.getU(), p3.getV());
         resultEdgeExists = weightedClusterEditingBranch(graph, k + cost);
-        graph.getEdgeWeights()[p3.getU()][p3.getV()] = cost;
-        graph.getEdgeWeights()[p3.getV()][p3.getU()] = cost;
-        graph.getAbsoluteNeighborhoodWeights()[p3.getU()] -= -FORBIDDEN_VALUE + cost;
-        graph.getAbsoluteNeighborhoodWeights()[p3.getV()] -= -FORBIDDEN_VALUE + cost;
-        graph.flipEdge(p3.getU(), p3.getV());
+        graph.flipBackForbiddenEdge(p3.getU(), p3.getV(), cost);
 
         DataReduction.revertHeavyNonEdgeReduction(graph, originalWeightsBeforeHeavyNonEdgeReduction);
 
         return resultEdgeExists;
     }
 
+    public static int getLowerBound2BasedOnMaximumWeightIndependentSetIlp(Graph graph, List<P3> p3List) {
+        int[] vertexWeights = new int[p3List.size()];
+        boolean[][] edges = new boolean[vertexWeights.length][vertexWeights.length];
+
+        for (int i = 0; i < vertexWeights.length; i++) {
+            vertexWeights[i] = graph.getSmallestAbsoluteWeight(p3List.get(i));
+        }
+
+        for (int i = 0; i < vertexWeights.length; i++) {
+            for (int j = i+1; j < vertexWeights.length; j++) {
+                boolean edgeExists = graph.getNumberOfSharedVertices(p3List.get(i), p3List.get(j)) > 1;
+                edges[i][j] = edgeExists;
+                edges[j][i] = edgeExists;
+            }
+        }
+
+        try {
+            GRBEnv env = new GRBEnv(true);
+            env.set(GRB.IntParam.LogToConsole, 0);
+            env.start();
+
+            GRBModel model = new GRBModel(env);
+
+            GRBVar[] x = new GRBVar[vertexWeights.length];
+
+            for (int i = 0; i < vertexWeights.length; i++) {
+                x[i] = model.addVar(0, 1, vertexWeights[i], GRB.INTEGER, String.valueOf(i));
+            }
+
+            model.set(GRB.IntAttr.ModelSense, GRB.MAXIMIZE);
+
+            for (int i = 0; i < vertexWeights.length; i++) {
+                for (int j = 0; j < vertexWeights.length; j++) {
+                    if (i != j && edges[i][j]) {
+                        GRBLinExpr constraint = new GRBLinExpr();
+                        constraint.addTerm(1., x[i]);
+                        constraint.addTerm(1., x[j]);
+                        model.addConstr(1., GRB.GREATER_EQUAL, constraint, String.format("%d %d", i, j));
+                    }
+                }
+            }
+
+            model.optimize();
+
+            return (int) Math.round(model.get(GRB.DoubleAttr.ObjVal));
+        } catch (GRBException e) {
+            e.printStackTrace();
+        }
+
+        return 0;
+    }
+
     private static P3 getBiggestWeightP3(Graph graph, List<P3> p3List) {
         int biggestWeight = Integer.MIN_VALUE;
         P3 biggestWeightP3 = null;
@@ -268,10 +312,6 @@ public class Main {
             return new ResultEdgeExistsWithSolutionSize(Utils.copy(graph.getEdgeExists(), graph.getNumberOfVertices()), costToEdit);
         }
 
-        if (costToEdit + graph.getLowerBound2(p3List) >= upperBound) {
-            return new ResultEdgeExistsWithSolutionSize(upperBoundSolutionEdgeExists, upperBound);
-        }
-
         // call recursively for all connected components
         ArrayList<ArrayList<Integer>> connectedComponents = graph.getConnectedComponents();
         if (connectedComponents.size() > 1) {
@@ -290,9 +330,9 @@ public class Main {
 
                 boolean[][] parentUpperBoundEdgeExists = getEdgeExistsOfSubGraph(upperBoundSolutionEdgeExists, subGraphIndices);
                 //to test: quick:
-                int parentUpperBoundCost = upperBound - costToEdit; //reduce because the costToEdit would be added for each component otherwise
+//                int parentUpperBoundCost = upperBound - costToEdit; //reduce because the costToEdit would be added for each component otherwise
                 //to test: middle slow/quick:
-//                int parentUpperBoundCost = Math.min( upperBound - costToEdit, Utils.getCostToChange(subGraph, parentUpperBoundEdgeExists) );
+                int parentUpperBoundCost = Math.min( upperBound - costToEdit, Utils.getCostToChange(subGraph, parentUpperBoundEdgeExists) );
                 //to test: slow:
                 //TODO use parameter or other heuristic: don't recalculate upper bound ALWAYS again
 //                ResultEdgeExistsWithSolutionSize resultEdgeExistsWithSolutionSizeOfUpperBound = getUpperBound(subGraph, 4, false);
@@ -324,39 +364,78 @@ public class Main {
             return new ResultEdgeExistsWithSolutionSize(resultEdgeExists, costToEdit);
         }
 
+        //LP lower bound
+//        if (recursiveSteps % 2 == 1){ //TODO use level or probablity instead?
+        double lowerBound = LpLowerBound.getLowerBoundOrTools(graph);
+        if (costToEdit + lowerBound >= upperBound) {
+            return new ResultEdgeExistsWithSolutionSize(upperBoundSolutionEdgeExists, upperBound);
+        }
+//        }
+        //Lower Bound 2
+//        if (recursiveSteps % 2 == 0){ //TODO use level or probablity instead?
+//            if (costToEdit + graph.getLowerBound2(p3List) >= upperBound) {
+//                return new ResultEdgeExistsWithSolutionSize(upperBoundSolutionEdgeExists, upperBound);
+//            }
+//        }
+
         P3 p3 = getBiggestWeightP3(graph, p3List);
 
         //merge or delete, TODO heuristic which to do first
         MergeVerticesInfo mergeVerticesInfo = graph.mergeVertices(p3.getU(), p3.getV());
-        ResultEdgeExistsWithSolutionSize solutionAfterMerge = weightedClusterEditingOptim(
-                graph, costToEdit + mergeVerticesInfo.getCost(), upperBoundSolutionEdgeExists, upperBound
-        );
-        if (solutionAfterMerge.getSolutionSize() < upperBound) {
-            upperBoundSolutionEdgeExists = Utils.reconstructMergeForResultEdgeExists(
-                    solutionAfterMerge.getResultEdgeExists(), graph, mergeVerticesInfo
-            );
-            upperBound = solutionAfterMerge.getSolutionSize();
-        }
+        ResultEdgeExistsWithSolutionSize resultEdgeExistsWithSolutionSizeUpperBoundAfterMerge = getUpperBound(graph, 16, false);
         graph.revertMergeVertices(mergeVerticesInfo);
 
-        graph.flipEdge(p3.getU(), p3.getV());
-        int costToFlip = graph.getEdgeWeights()[p3.getU()][p3.getV()];
-        graph.getEdgeWeights()[p3.getU()][p3.getV()] = FORBIDDEN_VALUE;
-        graph.getEdgeWeights()[p3.getV()][p3.getU()] = FORBIDDEN_VALUE;
-        graph.getAbsoluteNeighborhoodWeights()[p3.getU()] += -FORBIDDEN_VALUE + costToFlip;
-        graph.getAbsoluteNeighborhoodWeights()[p3.getV()] += -FORBIDDEN_VALUE + costToFlip;
-        ResultEdgeExistsWithSolutionSize solutionAfterDeletion = weightedClusterEditingOptim(
-                graph, costToEdit - costToFlip, upperBoundSolutionEdgeExists, upperBound
-        );
-        graph.getEdgeWeights()[p3.getU()][p3.getV()] = costToFlip;
-        graph.getEdgeWeights()[p3.getV()][p3.getU()] = costToFlip;
-        graph.getAbsoluteNeighborhoodWeights()[p3.getU()] -= -FORBIDDEN_VALUE + costToFlip;
-        graph.getAbsoluteNeighborhoodWeights()[p3.getV()] -= -FORBIDDEN_VALUE + costToFlip;
-        graph.flipEdge(p3.getU(), p3.getV());
-
-        if (solutionAfterDeletion.getSolutionSize() < upperBound) {
-            upperBoundSolutionEdgeExists = solutionAfterDeletion.getResultEdgeExists();
-            upperBound = solutionAfterDeletion.getSolutionSize();
+        int costToFlip = graph.flipEdgeAndSetForbidden(p3.getU(), p3.getV());
+        ResultEdgeExistsWithSolutionSize resultEdgeExistsWithSolutionsSizeUpperBoundAfterDeletion = getUpperBound(graph, 16, false);
+        graph.flipBackForbiddenEdge(p3.getU(), p3.getV(), costToFlip);
+
+        if (resultEdgeExistsWithSolutionsSizeUpperBoundAfterDeletion.getSolutionSize() < resultEdgeExistsWithSolutionSizeUpperBoundAfterMerge.getSolutionSize()) {
+            costToFlip = graph.flipEdgeAndSetForbidden(p3.getU(), p3.getV());
+            ResultEdgeExistsWithSolutionSize solutionAfterDeletion = weightedClusterEditingOptim(
+                    graph, costToEdit - costToFlip, upperBoundSolutionEdgeExists, upperBound
+            );
+            if (solutionAfterDeletion.getSolutionSize() < upperBound) {
+                upperBoundSolutionEdgeExists = solutionAfterDeletion.getResultEdgeExists();
+                upperBound = solutionAfterDeletion.getSolutionSize();
+            }
+            graph.flipBackForbiddenEdge(p3.getU(), p3.getV(), costToFlip);
+
+            mergeVerticesInfo = graph.mergeVertices(p3.getU(), p3.getV());
+            ResultEdgeExistsWithSolutionSize solutionAfterMerge = weightedClusterEditingOptim(
+                    graph, costToEdit + mergeVerticesInfo.getCost(), upperBoundSolutionEdgeExists, upperBound
+            );
+
+            if (solutionAfterMerge.getSolutionSize() < upperBound) {
+                upperBoundSolutionEdgeExists = Utils.reconstructMergeForResultEdgeExists(
+                        solutionAfterMerge.getResultEdgeExists(), graph, mergeVerticesInfo
+                );
+                upperBound = solutionAfterMerge.getSolutionSize();
+            }
+
+            graph.revertMergeVertices(mergeVerticesInfo);
+        } else {
+            mergeVerticesInfo = graph.mergeVertices(p3.getU(), p3.getV());
+            ResultEdgeExistsWithSolutionSize solutionAfterMerge = weightedClusterEditingOptim(
+                    graph, costToEdit + mergeVerticesInfo.getCost(), upperBoundSolutionEdgeExists, upperBound
+            );
+            if (solutionAfterMerge.getSolutionSize() < upperBound) {
+                upperBoundSolutionEdgeExists = Utils.reconstructMergeForResultEdgeExists(
+                        solutionAfterMerge.getResultEdgeExists(), graph, mergeVerticesInfo
+                );
+                upperBound = solutionAfterMerge.getSolutionSize();
+            }
+            graph.revertMergeVertices(mergeVerticesInfo);
+
+            costToFlip = graph.flipEdgeAndSetForbidden(p3.getU(), p3.getV());
+            ResultEdgeExistsWithSolutionSize solutionAfterDeletion = weightedClusterEditingOptim(
+                    graph, costToEdit - costToFlip, upperBoundSolutionEdgeExists, upperBound
+            );
+            graph.flipBackForbiddenEdge(p3.getU(), p3.getV(), costToFlip);
+
+            if (solutionAfterDeletion.getSolutionSize() < upperBound) {
+                upperBoundSolutionEdgeExists = solutionAfterDeletion.getResultEdgeExists();
+                upperBound = solutionAfterDeletion.getSolutionSize();
+            }
         }
 
         return new ResultEdgeExistsWithSolutionSize(upperBoundSolutionEdgeExists, upperBound);

algorithm-engineering-search-tree/src/main/java/berlin/tu/algorithmengineering/searchtree/lp/LpLowerBound.java

@@ -6,13 +6,13 @@ import com.google.ortools.linearsolver.MPConstraint;
 import com.google.ortools.linearsolver.MPObjective;
 import com.google.ortools.linearsolver.MPSolver;
 import com.google.ortools.linearsolver.MPVariable;
-import gurobi.*;
+//import gurobi.*;
 
 import java.util.Set;
 
 
 public class LpLowerBound {
-
+/*
     public static double getLowerBound(Graph graph) {
         return getLowerBound(graph, false);
     }
@@ -72,6 +72,7 @@ public class LpLowerBound {
         }
         return Double.MAX_VALUE;
     }
+    */
 
 //    public static double getLowerBoundSubGraphs(Graph graph) {
 //        final int SUB_GRAPH_SIZE = 50;
@@ -111,7 +112,7 @@ public class LpLowerBound {
 
         for (int i = 0; i < graph.getNumberOfVertices(); i++) {
             for (int j = i+1; j < graph.getNumberOfVertices(); j++) {
-                x[i][j] = mpSolver.makeNumVar(0.0, 1.0, String.format("%d %d", graph.getNumberOfVertices(), graph.getNumberOfVertices()));
+                x[i][j] = mpSolver.makeNumVar(0, 1, String.format("%d %d", graph.getNumberOfVertices(), graph.getNumberOfVertices()));
             }
         }
 
@@ -119,10 +120,10 @@ public class LpLowerBound {
             for (int j = 0; j < graph.getNumberOfVertices(); j++) {
                 for (int k = 0; k < graph.getNumberOfVertices(); k++) {
                     if (i != j && j != k && i != k) {
-                        MPConstraint constraint = mpSolver.makeConstraint(Double.NEGATIVE_INFINITY, 1.0);
-                        constraint.setCoefficient(x[Math.min(i, j)][Math.max(i, j)], 1.0);
-                        constraint.setCoefficient(x[Math.min(j, k)][Math.max(j, k)], 1.0);
-                        constraint.setCoefficient(x[Math.min(i, k)][Math.max(i, k)], -1.0);
+                        MPConstraint constraint = mpSolver.makeConstraint(Double.NEGATIVE_INFINITY, 1);
+                        constraint.setCoefficient(x[Math.min(i, j)][Math.max(i, j)], 1);
+                        constraint.setCoefficient(x[Math.min(j, k)][Math.max(j, k)], 1);
+                        constraint.setCoefficient(x[Math.min(i, k)][Math.max(i, k)], -1);
                     }
                 }
             }
@@ -133,17 +134,21 @@ public class LpLowerBound {
 
         for (int i = 0; i < graph.getNumberOfVertices(); i++) {
             for (int j = i + 1; j < graph.getNumberOfVertices(); j++) {
-                objective.setCoefficient(x[i][j], -graph.getEdgeWeights()[i][j]);
-                constant += Math.max(graph.getEdgeWeights()[i][j], 0.0);
+                if (graph.getEdgeWeights()[i][j] != 0) {
+                    objective.setCoefficient(x[i][j], -graph.getEdgeWeights()[i][j]);
+                    if (graph.getEdgeExists()[i][j]) {
+                        constant += graph.getEdgeWeights()[i][j];
+                    }
+                }
             }
         }
 
-        objective.setOffset(constant);
+        //objective.setOffset(constant);
 
         objective.setMinimization();
 
         MPSolver.ResultStatus resultStatus = mpSolver.solve();
-
-        return objective.value();
+//        System.out.printf("#walltime %d: %s\n",graph.getNumberOfVertices(),mpSolver.wallTime());
+        return objective.value() + constant;
     }
 }

scripts/compute-averages.py

+import pandas as pd
+import sys
+
+
+INSTANCE_TYPE = 'action-seq'
+
+
+def type_of(filename):
+    if 'random' in filename:
+        return 'random'
+    if 'real-world' in filename:
+        return 'real-world'
+    if 'actionseq' in filename:
+        return 'action-seq'
+
+
+if len(sys.argv) < 2:
+    raise Exception('Number of files needs to be provided')
+
+if len(sys.argv) < int(sys.argv[1]) + 1:
+    raise Exception(f"Input filename{'' if int(sys.argv[1]) == 1 else 's'} of {int(sys.argv[1])} file{'' if int(sys.argv[1]) == 1 else 's'} need to be provided as console arguments")
+
+data = []
+
+for i in range(int(sys.argv[1])):
+    data.append(pd.read_csv(sys.argv[i+2], delimiter=';', index_col=False))
+
+total_time = 0
+total_recsteps = 0
+instances = 0
+
+
+def correct_in_all_files(filename, data):
+    for i in range(len(data)):
+        row = data[i][data[i]['file'] == filename]
+        if row.empty or row['verified'].item() != 'correct':
+            return False
+    return True
+
+
+for index, row in data[0].iterrows():
+    if INSTANCE_TYPE is None or type_of(row['file']) == INSTANCE_TYPE:
+        if correct_in_all_files(row['file'], data):
+            total_time += row['time']
+            total_recsteps += row['recsteps']
+            instances += 1
+
+
+time_per_instance = total_time / instances
+recsteps_per_instance = total_recsteps / instances
+
+print(f"time per instance: {time_per_instance}")
+print(f"recursive steps per instance: {recsteps_per_instance}")
+print(f"time per recursive step: {1000 * time_per_instance / recsteps_per_instance}")

scripts/compute-lower-bound-solution.py

+import pandas as pd
+import sys
+
+
+if len(sys.argv) < 2:
+    raise Exception('Input filename needs to be provided.')
+
+data = []
+
+for i in range(1, len(sys.argv)):
+    data.append(pd.read_csv(sys.argv[i], delimiter=';'))
+
+
+instances = []
+
+for index, row in data[0].iterrows():
+    solved_in_all_files = True
+    for dataframe in data:
+        other_row = dataframe[dataframe['file'] == row['file']]
+        if other_row.empty or (other_row['verified'].item() != 'correct' and other_row['verified'].item() != '>>BAD COST<< '):
+            solved_in_all_files = False
+            break
+    if solved_in_all_files:
+        exact_solution = row['solsize'] - (0 if row['verified'] == 'correct' else int(row['costdiff']))
+        solution_ratios = []
+        for dataframe in data:
+            other_row = dataframe[dataframe['file'] == row['file']]
+            solution_ratios.append(other_row['recsteps'].item() / exact_solution if exact_solution > 0 else 1)
+        instances.append(solution_ratios)
+
+result = pd.DataFrame(instances)
+
+result.to_csv('lb_solutions.csv', sep=';', index=True)

scripts/compute-relative-solution.py

+import pandas as pd
+import sys
+
+
+if len(sys.argv) < 3:
+    raise Exception('Input filename needs to be provided.')
+
+upper_bound_data = pd.read_csv(sys.argv[1], delimiter=';')
+lower_bound_data = pd.read_csv(sys.argv[2], delimiter=';')
+
+instances = []
+
+for index, upper_bound_row in upper_bound_data.iterrows():
+    if upper_bound_row['finished'] == 1 and (upper_bound_row['verified'] == 'correct' or upper_bound_row['verified'] == '>>BAD COST<< '):
+        file = upper_bound_row['file']
+        lower_bound_row = lower_bound_data[lower_bound_data['file'] == file]
+        if not lower_bound_row.empty and (lower_bound_row['verified'].item() == 'correct' or lower_bound_row['verified'].item() == '>>BAD COST<< '):
+            heuristic_solution = upper_bound_row['solsize']
+            exact_solution = heuristic_solution - (0 if upper_bound_row['verified'] == 'correct' else upper_bound_row['costdiff'])
+            lower_bound_solution = lower_bound_row['recsteps'].item()
+            if exact_solution > 0:
+                instances.append([file, heuristic_solution / exact_solution, lower_bound_solution / exact_solution])
+
+result = pd.DataFrame(instances, columns=['file', 'upperbound', 'lowerbound'])
+
+result.to_csv('lb_ub.csv', sep=';', index=True)

scripts/compute_heuristic_score.py

+import pandas as pd
+import sys
+
+
+if len(sys.argv) < 2:
+    raise Exception('Input filename needs to be provided.')
+
+data = []
+
+for i in range(1, len(sys.argv)):
+    data.append(pd.read_csv(sys.argv[i], delimiter=';'))
+
+
+avg_score = 0
+avg_running_time = 0
+n = 0
+
+
+def finished_in_all_files(filename, data):
+    for i in range(len(data)):
+        row = data[i][data[i]['file'] == filename]
+        if row.empty or not (row['finished'].item() == 1 and (row['verified'].item() == 'correct' or row['verified'].item() == '>>BAD COST<< ')):
+            return False
+    return True
+
+
+indices = []
+
+for index, row in data[0].iterrows():
+    verified = row['verified']
+    if finished_in_all_files(row['file'], data):
+        n += 1
+        indices.append(index)
+        solsize = row['solsize']
+        score = 1 if verified == 'correct' else (solsize - float(row['costdiff'])) / solsize
+        avg_score = (n-1) / n * avg_score + 1/n * score
+        avg_running_time = (n-1) / n * avg_running_time + 1/n * row['time']
+
+
+print(n)
+print(f'avg score = {avg_score}')
+print(f'avg running time = {avg_running_time}')

scripts/compute_heuristic_score_sa_values.py

+import pandas as pd
+import sys
+
+
+if len(sys.argv) < 2:
+    raise Exception('Input filename needs to be provided.')
+
+data = []
+
+for i in range(1, len(sys.argv)):
+    data.append(pd.read_csv(sys.argv[i], delimiter=';'))
+
+
+avg_score = 0
+avg_running_time = 0
+n = 0
+
+
+def finished_in_all_files(filename, data):
+    for i in range(len(data)):
+        row = data[i][data[i]['file'] == filename]
+        sa_values_str = row['recsteps'].item()
+        if not isinstance(sa_values_str, str):
+            return False
+        sa_values = sa_values_str.split(',')
+        if row.empty or not (row['finished'].item() == 1 and (row['verified'].item() == 'correct' or row['verified'].item() == '>>BAD COST<< ') and len(sa_values) == 30_000):
+            return False
+    return True
+
+
+indices = []
+
+for index, row in data[0].iterrows():
+    verified = row['verified']
+    if index == 149:
+        print()
+    if finished_in_all_files(row['file'], data):
+        n += 1
+        indices.append(index)
+        solsize = row['solsize']
+        score = 1 if verified == 'correct' else (solsize - row['costdiff']) / solsize
+        avg_score = (n-1) / n * avg_score + 1/n * score
+        avg_running_time = (n-1) / n * avg_running_time + 1/n * row['time']
+
+
+print(n)
+print(f'avg score = {avg_score}')
+print(f'avg running time = {avg_running_time}')

scripts/create-cactus-plot-file.py

+import pandas as pd
+import sys
+
+data = []
+
+for i in range(1, len(sys.argv)):
+    df = pd.read_csv(sys.argv[i], delimiter=';', index_col=False)
+    values = []
+
+    for index, row in df.iterrows():
+        if row['verified'] == 'correct':
+            values.append(row['time'])
+
+    values.sort()
+
+    data.append(values)
+
+compare_df = pd.DataFrame(data).transpose()
+compare_df.set_axis([chr(65 + i) + "-sorted" for i in range(len(data))], axis=1, inplace=True)
+
+compare_df.to_csv('compare-solvers.csv', index=False, sep=';')
+
+

scripts/create-compare.py

+import pandas as pd
+import sys
+
+
+TIMEOUT = 30
+
+
+def type_of(filename):
+    if 'random' in filename:
+        return 'random'
+    if 'real-world' in filename:
+        return 'real-world'
+    if 'actionseq' in filename:
+        return 'action-seq'
+
+
+if len(sys.argv) < 3:
+    raise Exception('Input filenames of both files need to be provided as console arguments')
+
+data1 = pd.read_csv(sys.argv[1], delimiter=';', index_col=False)
+data2 = pd.read_csv(sys.argv[2], delimiter=';', index_col=False)
+
+compare = []
+
+for index, row in data1.iterrows():
+    row2 = data2[data2['file'] == row['file']]
+    if row['verified'] == 'correct' and (not row2.empty and row2['verified'].item() == 'correct'):
+        compare.append([
+            row['file'],
+            type_of(row['file']),
+            row['time'] if row['verified'] == 'correct' else TIMEOUT,
+            row2['time'].item() if not row2.empty and row2['verified'].item() == 'correct' else TIMEOUT,
+            row['recsteps'] if row['verified'] == 'correct' else None,
+            row2['recsteps'].item() if not row2.empty and row2['verified'].item() == 'correct' else None
+        ])
+
+pd.DataFrame(compare, columns=['file', 'Type', 'Atime', 'Btime', 'Arecsteps', 'Brecsteps']).to_csv('compare.csv', index=False, sep=';')
+

scripts/create-sa-csv.py

+import math
+import matplotlib.pyplot as plt
+import pandas as pd
+import numpy as np
+import sys
+
+
+if len(sys.argv) < 2:
+    raise Exception('Input filename needs to be provided.')
+
+df = pd.read_csv(sys.argv[1], delimiter=';')
+
+data = np.zeros((0, 30_000))
+
+for index, row in df.iterrows():
+    if row['verified'] == 'correct' or row['verified'] == '>>BAD COST<< ':
+        sa_values_str = row['recsteps']
+        if isinstance(sa_values_str, str):
+            sa_values = sa_values_str.split(',')
+            costdiff = row['costdiff']
+            exact_solution = row['solsize'] - (0 if math.isnan(costdiff) else float(costdiff))
+            data = np.vstack((data, np.asarray([[exact_solution / int(cost) if exact_solution != 0 else int(cost) for cost in sa_values]])))
+
+plt.plot(np.arange(0, 30_000, 1), data.mean(axis=0))
+plt.show()
+