diff --git a/algorithm-engineering-common/src/main/java/berlin/tu/algorithmengineering/common/Graph.java b/algorithm-engineering-common/src/main/java/berlin/tu/algorithmengineering/common/Graph.java index 399e886e0c8b2abf5b13ad98c2ec1bd84c3d6f2a..76a672c77469554f9d2636496524e157bb4f91b5 100644 --- a/algorithm-engineering-common/src/main/java/berlin/tu/algorithmengineering/common/Graph.java +++ b/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 */ diff --git a/algorithm-engineering-heuristics/src/main/java/berlin/tu/algorithmengineering/heuristics/SimulatedAnnealing.java b/algorithm-engineering-heuristics/src/main/java/berlin/tu/algorithmengineering/heuristics/SimulatedAnnealing.java index 2b7e2410ebf5213eeded4072bfb6694db4087954..cca4b670b91aae54779f66f143e33a96facc27ac 100644 --- a/algorithm-engineering-heuristics/src/main/java/berlin/tu/algorithmengineering/heuristics/SimulatedAnnealing.java +++ b/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++) { diff --git a/algorithm-engineering-search-tree/src/main/java/berlin/tu/algorithmengineering/searchtree/Main.java b/algorithm-engineering-search-tree/src/main/java/berlin/tu/algorithmengineering/searchtree/Main.java index aeae47f514eda4c388bc02700b4c0665b65849ce..13ddf2b9795967494344d6e674693ccf308bc51c 100644 --- a/algorithm-engineering-search-tree/src/main/java/berlin/tu/algorithmengineering/searchtree/Main.java +++ b/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); diff --git a/algorithm-engineering-search-tree/src/main/java/berlin/tu/algorithmengineering/searchtree/lp/LpLowerBound.java b/algorithm-engineering-search-tree/src/main/java/berlin/tu/algorithmengineering/searchtree/lp/LpLowerBound.java index 437bf1779f7b0e22f158f691940414f7085412d6..ca28ae4ed6da32e6290c21fc16451a3b209e4bee 100644 --- a/algorithm-engineering-search-tree/src/main/java/berlin/tu/algorithmengineering/searchtree/lp/LpLowerBound.java +++ b/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; } } diff --git a/scripts/compute-averages.py b/scripts/compute-averages.py new file mode 100644 index 0000000000000000000000000000000000000000..b5f242e7808c9107fcac5fd9d459050b2816899c --- /dev/null +++ b/scripts/compute-averages.py @@ -0,0 +1,54 @@ +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}") diff --git a/scripts/compute-lower-bound-solution.py b/scripts/compute-lower-bound-solution.py new file mode 100644 index 0000000000000000000000000000000000000000..e90aa8d30b9c29016e01c5d686141803dff54d0e --- /dev/null +++ b/scripts/compute-lower-bound-solution.py @@ -0,0 +1,33 @@ +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) diff --git a/scripts/compute-relative-solution.py b/scripts/compute-relative-solution.py new file mode 100644 index 0000000000000000000000000000000000000000..f85b9a512e033b04c39c1b7007009c9856f59144 --- /dev/null +++ b/scripts/compute-relative-solution.py @@ -0,0 +1,26 @@ +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) diff --git a/scripts/compute_heuristic_score.py b/scripts/compute_heuristic_score.py new file mode 100644 index 0000000000000000000000000000000000000000..c6119729e5954d5cf7a930c64011fdd768aac559 --- /dev/null +++ b/scripts/compute_heuristic_score.py @@ -0,0 +1,42 @@ +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}') diff --git a/scripts/compute_heuristic_score_sa_values.py b/scripts/compute_heuristic_score_sa_values.py new file mode 100644 index 0000000000000000000000000000000000000000..4c0a4a05e70ae52aad1acf8ca54338b9a19d2ede --- /dev/null +++ b/scripts/compute_heuristic_score_sa_values.py @@ -0,0 +1,48 @@ +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}') diff --git a/scripts/create-cactus-plot-file.py b/scripts/create-cactus-plot-file.py new file mode 100644 index 0000000000000000000000000000000000000000..f48883106838ef96b36521fa3c502d1320c3b658 --- /dev/null +++ b/scripts/create-cactus-plot-file.py @@ -0,0 +1,23 @@ +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=';') + + diff --git a/scripts/create-compare.py b/scripts/create-compare.py new file mode 100644 index 0000000000000000000000000000000000000000..e879a3713804290d3a2931bab513fc0040526902 --- /dev/null +++ b/scripts/create-compare.py @@ -0,0 +1,38 @@ +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=';') + diff --git a/scripts/create-sa-csv.py b/scripts/create-sa-csv.py new file mode 100644 index 0000000000000000000000000000000000000000..abfc0eff7b201abc1e1bface021ea825997404bf --- /dev/null +++ b/scripts/create-sa-csv.py @@ -0,0 +1,26 @@ +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() +