From 139e51c73b3746247c7ac3e0aed0829d3f1ae746 Mon Sep 17 00:00:00 2001 From: Piotr Malkowski Date: Mon, 10 Mar 2025 11:44:14 +0100 Subject: [PATCH 01/37] feat: add LP / MIP Problems to backend, add google OR-Tools MIP Problem wrapper --- .../provideq/toolbox/lp/LpConfiguration.java | 63 +++++++++++++++++++ .../provideq/toolbox/lp/solvers/LpSolver.java | 12 ++++ .../toolbox/lp/solvers/OrToolsMipCbc.java | 58 +++++++++++++++++ 3 files changed, 133 insertions(+) create mode 100644 src/main/java/edu/kit/provideq/toolbox/lp/LpConfiguration.java create mode 100644 src/main/java/edu/kit/provideq/toolbox/lp/solvers/LpSolver.java create mode 100644 src/main/java/edu/kit/provideq/toolbox/lp/solvers/OrToolsMipCbc.java diff --git a/src/main/java/edu/kit/provideq/toolbox/lp/LpConfiguration.java b/src/main/java/edu/kit/provideq/toolbox/lp/LpConfiguration.java new file mode 100644 index 00000000..1347bfe0 --- /dev/null +++ b/src/main/java/edu/kit/provideq/toolbox/lp/LpConfiguration.java @@ -0,0 +1,63 @@ +package edu.kit.provideq.toolbox.lp; + +import edu.kit.provideq.toolbox.ResourceProvider; +import edu.kit.provideq.toolbox.exception.MissingExampleException; +import edu.kit.provideq.toolbox.lp.solvers.OrToolsMipCbc; +import edu.kit.provideq.toolbox.meta.Problem; +import edu.kit.provideq.toolbox.meta.ProblemManager; +import edu.kit.provideq.toolbox.meta.ProblemType; +import edu.kit.provideq.toolbox.qubo.solvers.DwaveQuboSolver; +import edu.kit.provideq.toolbox.qubo.solvers.KipuQuboSolver; +import edu.kit.provideq.toolbox.qubo.solvers.QiskitQuboSolver; +import edu.kit.provideq.toolbox.qubo.solvers.QrispQuboSolver; +import edu.kit.provideq.toolbox.qubo.solvers.QuantagoniaQuboSolver; +import java.io.IOException; +import java.util.Objects; +import java.util.Set; +import org.springframework.context.annotation.Bean; +import org.springframework.context.annotation.Configuration; + +/** + * Definition and registration of the "Quadratic Unconstrained Binary Optimization" problem. + */ +@Configuration +public class LpConfiguration { + /** + * QUBO (Quadratic Unconstrained Binary Optimization) + * A combinatorial optimization problem. + * For a given quadratic term with binary decision variables, + * find the minimal variable assignment of the term. + */ + public static final ProblemType LP = new ProblemType<>( + "lp", + String.class, + String.class + ); + + @Bean + ProblemManager getLpManager( + OrToolsMipCbc orToolsMipCbc, + ResourceProvider resourceProvider + ) { + return new ProblemManager<>( + LP, + Set.of(orToolsMipCbc), + loadExampleProblems(resourceProvider) + ); + } + + private Set> loadExampleProblems( + ResourceProvider resourceProvider) { + try { + var problemInputStream = Objects.requireNonNull( + getClass().getResourceAsStream("qubo.lp"), + "quadratic-problem example for QUBO is unavailable!" + ); + var problem = new Problem<>(LP); + problem.setInput(resourceProvider.readStream(problemInputStream)); + return Set.of(problem); + } catch (IOException e) { + throw new MissingExampleException(LP, e); + } + } +} diff --git a/src/main/java/edu/kit/provideq/toolbox/lp/solvers/LpSolver.java b/src/main/java/edu/kit/provideq/toolbox/lp/solvers/LpSolver.java new file mode 100644 index 00000000..228936f5 --- /dev/null +++ b/src/main/java/edu/kit/provideq/toolbox/lp/solvers/LpSolver.java @@ -0,0 +1,12 @@ +package edu.kit.provideq.toolbox.lp.solvers; + +import edu.kit.provideq.toolbox.lp.LpConfiguration; +import edu.kit.provideq.toolbox.meta.ProblemSolver; +import edu.kit.provideq.toolbox.meta.ProblemType; + +public abstract class LpSolver implements ProblemSolver { + @Override + public ProblemType getProblemType() { + return LpConfiguration.LP; + } +} diff --git a/src/main/java/edu/kit/provideq/toolbox/lp/solvers/OrToolsMipCbc.java b/src/main/java/edu/kit/provideq/toolbox/lp/solvers/OrToolsMipCbc.java new file mode 100644 index 00000000..79b09a56 --- /dev/null +++ b/src/main/java/edu/kit/provideq/toolbox/lp/solvers/OrToolsMipCbc.java @@ -0,0 +1,58 @@ +package edu.kit.provideq.toolbox.lp.solvers; + +import edu.kit.provideq.toolbox.Solution; +import edu.kit.provideq.toolbox.meta.SolvingProperties; +import edu.kit.provideq.toolbox.meta.SubRoutineResolver; +import edu.kit.provideq.toolbox.process.ProcessRunner; +import edu.kit.provideq.toolbox.process.PythonProcessRunner; +import org.springframework.beans.factory.annotation.Autowired; +import org.springframework.beans.factory.annotation.Value; +import org.springframework.context.ApplicationContext; +import org.springframework.stereotype.Component; +import reactor.core.publisher.Mono; + +@Component +public class OrToolsMipCbc extends LpSolver { + private final String scriptPath; + private final ApplicationContext context; + + @Autowired + public OrToolsMipCbc( + @Value("${ortools.script.ormip}") String scriptPath, + ApplicationContext context) { + this.scriptPath = scriptPath; + this.context = context; + } + + @Override + public String getName() { + return "(OR-Tools) Cbc Solver for MIP"; + } + + @Override + public String getDescription() { + return "This solver uses OR-Tools CBC Solver to solve MIP problems"; + } + + @Override + public Mono> solve( + String input, + SubRoutineResolver subRoutineResolver, + SolvingProperties properties + ) { + var solution = new Solution<>(this); + + var processResult = context + .getBean(PythonProcessRunner.class, scriptPath) + .withArguments( + ProcessRunner.INPUT_FILE_PATH, + ProcessRunner.OUTPUT_FILE_PATH + ) + .writeInputFile(input, "problem.lp") + .readOutputFile() + .run(getProblemType(), solution.getId()); + + // Return if process failed + return Mono.just(processResult.applyTo(solution)); + } +} From 7e2e249c730682f7752faf289252a776f0b6ab50 Mon Sep 17 00:00:00 2001 From: Piotr Malkowski Date: Mon, 10 Mar 2025 11:45:15 +0100 Subject: [PATCH 02/37] feat: add OR-Tools MIP wrapper using build-in cbc solver --- solvers/ortools/mips/cbc.py | 33 ++++++++++++++++++++++++++++++++ solvers/ortools/requirements.txt | 2 ++ 2 files changed, 35 insertions(+) create mode 100644 solvers/ortools/mips/cbc.py create mode 100644 solvers/ortools/requirements.txt diff --git a/solvers/ortools/mips/cbc.py b/solvers/ortools/mips/cbc.py new file mode 100644 index 00000000..ca717a5d --- /dev/null +++ b/solvers/ortools/mips/cbc.py @@ -0,0 +1,33 @@ +import argparse +import os +from ortools.linear_solver import pywraplp +from lp_mips_converter import convert_model + + +parser = argparse.ArgumentParser( + prog="OR-Tools CBC MIP Solver", + description="A CLI Program to solve MIP problems in lp or mips format with Or-Tools CBC" +) + +parser.add_argument("input_file") +parser.add_argument("--output-file") +args = parser.parse_args() + +input_ext = os.path.splitext(args.input_file)[1].lower() +if input_ext == ".lp": + convert_model(args.input_file, args.output_file) + mps_file = args.output_file +else: + mps_file = args.input_file + + +solver = pywraplp.Solver.CreateSolver("CBC") +model = solver.loadModelFromFile(mps_file) +status = solver.Solve() + +if status == pywraplp.Solver.OPTIMAL: + print("Solution:") + print("Objective value =", solver.Objective().Value()) +else: + print("The problem does not have an optimal solution.") + diff --git a/solvers/ortools/requirements.txt b/solvers/ortools/requirements.txt new file mode 100644 index 00000000..5367027c --- /dev/null +++ b/solvers/ortools/requirements.txt @@ -0,0 +1,2 @@ +ortools +cplex From 93d038cc922db6e578f9fe0993e0853c5e18ce36 Mon Sep 17 00:00:00 2001 From: Piotr Malkowski Date: Mon, 10 Mar 2025 11:45:46 +0100 Subject: [PATCH 03/37] feat: add CPLEX converter from and to .lp and .mps files --- solvers/ortools/mips/lp_mips_converter.py | 11 +++++++++++ 1 file changed, 11 insertions(+) create mode 100644 solvers/ortools/mips/lp_mips_converter.py diff --git a/solvers/ortools/mips/lp_mips_converter.py b/solvers/ortools/mips/lp_mips_converter.py new file mode 100644 index 00000000..2896c21a --- /dev/null +++ b/solvers/ortools/mips/lp_mips_converter.py @@ -0,0 +1,11 @@ +import cplex + + +def convert_model(input_file, output_file): + """ + Loads a model (LP or MPS) into CPLEX and writes it back out + in the requested format (.lp or .mps). + """ + # Create a CPLEX object and read in the model + model = cplex.Cplex(input_file) + model.write(output_file) From f7899713ba22c53b8504c0fe5b966478c38ef14b Mon Sep 17 00:00:00 2001 From: Piotr Malkowski Date: Mon, 10 Mar 2025 11:46:04 +0100 Subject: [PATCH 04/37] chore: add paths for OR-Tools solvers --- src/main/resources/application.properties | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/main/resources/application.properties b/src/main/resources/application.properties index f3659bb8..a717df3d 100644 --- a/src/main/resources/application.properties +++ b/src/main/resources/application.properties @@ -21,6 +21,9 @@ qiskit.script.max-cut=${qiskit.directory}/max-cut/maxCut_qiskit.py qiskit.script.qubo=${qiskit.directory}/qubo/qubo_qiskit.py qiskit.script.knapsack=${qiskit.directory}/knapsack/knapsack_qiskit.py +ortools.directory=${solvers.directory}/ortools +ortools.script.ormip=${ortools.directory}/mips/cbc.py + cirq.directory=${solvers.directory}/cirq cirq.script.max-cut=${cirq.directory}/max-cut/max_cut_cirq.py From c50ef0f98796a7b4197174701fdc96cfe981206d Mon Sep 17 00:00:00 2001 From: Piotr Malkowski Date: Mon, 10 Mar 2025 11:47:05 +0100 Subject: [PATCH 05/37] refactor: remove unused statements --- .../java/edu/kit/provideq/toolbox/lp/LpConfiguration.java | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/main/java/edu/kit/provideq/toolbox/lp/LpConfiguration.java b/src/main/java/edu/kit/provideq/toolbox/lp/LpConfiguration.java index 1347bfe0..a600a797 100644 --- a/src/main/java/edu/kit/provideq/toolbox/lp/LpConfiguration.java +++ b/src/main/java/edu/kit/provideq/toolbox/lp/LpConfiguration.java @@ -6,11 +6,6 @@ import edu.kit.provideq.toolbox.meta.Problem; import edu.kit.provideq.toolbox.meta.ProblemManager; import edu.kit.provideq.toolbox.meta.ProblemType; -import edu.kit.provideq.toolbox.qubo.solvers.DwaveQuboSolver; -import edu.kit.provideq.toolbox.qubo.solvers.KipuQuboSolver; -import edu.kit.provideq.toolbox.qubo.solvers.QiskitQuboSolver; -import edu.kit.provideq.toolbox.qubo.solvers.QrispQuboSolver; -import edu.kit.provideq.toolbox.qubo.solvers.QuantagoniaQuboSolver; import java.io.IOException; import java.util.Objects; import java.util.Set; From fc914861588923848ef25c028b0ddb7b6432dae9 Mon Sep 17 00:00:00 2001 From: Piotr Malkowski Date: Mon, 10 Mar 2025 11:50:49 +0100 Subject: [PATCH 06/37] feat: add simple resource problem --- .../provideq/toolbox/lp/LpConfiguration.java | 4 ++-- .../edu/kit/provideq/toolbox/lp/simple.mps | 18 ++++++++++++++++++ 2 files changed, 20 insertions(+), 2 deletions(-) create mode 100644 src/main/resources/edu/kit/provideq/toolbox/lp/simple.mps diff --git a/src/main/java/edu/kit/provideq/toolbox/lp/LpConfiguration.java b/src/main/java/edu/kit/provideq/toolbox/lp/LpConfiguration.java index a600a797..f74e1e9c 100644 --- a/src/main/java/edu/kit/provideq/toolbox/lp/LpConfiguration.java +++ b/src/main/java/edu/kit/provideq/toolbox/lp/LpConfiguration.java @@ -45,8 +45,8 @@ private Set> loadExampleProblems( ResourceProvider resourceProvider) { try { var problemInputStream = Objects.requireNonNull( - getClass().getResourceAsStream("qubo.lp"), - "quadratic-problem example for QUBO is unavailable!" + getClass().getResourceAsStream("simple.mps"), + "mixed-integer-problem example for LP is unavailable!" ); var problem = new Problem<>(LP); problem.setInput(resourceProvider.readStream(problemInputStream)); diff --git a/src/main/resources/edu/kit/provideq/toolbox/lp/simple.mps b/src/main/resources/edu/kit/provideq/toolbox/lp/simple.mps new file mode 100644 index 00000000..24b8d833 --- /dev/null +++ b/src/main/resources/edu/kit/provideq/toolbox/lp/simple.mps @@ -0,0 +1,18 @@ +NAME EXAMPLE +ROWS + N COST <- (N) indicates the objective row + G C1 <- (G) 'Greater than or equal' constraint + L C2 <- (L) 'Less than or equal' constraint +COLUMNS + X COST 1 <- Objective coefficient for X is 1 + X C1 1 <- X's coefficient in constraint C1 + Y COST 3 <- Objective coefficient for Y is 3 + Y C1 1 <- Y's coefficient in constraint C1 + Y C2 1 <- Y's coefficient in constraint C2 +RHS + RHS C1 4 <- Right-hand side of constraint C1 + RHS C2 3 <- Right-hand side of constraint C2 +BOUNDS + LO BND X 0 <- Lower bound of X is 0 + LO BND Y 0 <- Lower bound of Y is 0 +ENDATA From bab5da1798413588c89f7ba7db1d2487e88feaba Mon Sep 17 00:00:00 2001 From: Piotr Malkowski Date: Mon, 10 Mar 2025 11:52:14 +0100 Subject: [PATCH 07/37] refactor: adapt problem description --- .../java/edu/kit/provideq/toolbox/lp/LpConfiguration.java | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/main/java/edu/kit/provideq/toolbox/lp/LpConfiguration.java b/src/main/java/edu/kit/provideq/toolbox/lp/LpConfiguration.java index f74e1e9c..465861a2 100644 --- a/src/main/java/edu/kit/provideq/toolbox/lp/LpConfiguration.java +++ b/src/main/java/edu/kit/provideq/toolbox/lp/LpConfiguration.java @@ -13,14 +13,14 @@ import org.springframework.context.annotation.Configuration; /** - * Definition and registration of the "Quadratic Unconstrained Binary Optimization" problem. + * Definition and registration of the "Linear and Mixed Integer" problems. */ @Configuration public class LpConfiguration { /** - * QUBO (Quadratic Unconstrained Binary Optimization) - * A combinatorial optimization problem. - * For a given quadratic term with binary decision variables, + * LP (Linear Problem) + * A linear optimization problem. + * For a given linear term, * find the minimal variable assignment of the term. */ public static final ProblemType LP = new ProblemType<>( From 011be8ddcc7b0dee566b277c9e2ddec53f2de9c2 Mon Sep 17 00:00:00 2001 From: Piotr Malkowski Date: Mon, 10 Mar 2025 12:08:31 +0100 Subject: [PATCH 08/37] refactor: rename lp to mip. --- solvers/ortools/mips/cbc.py | 2 +- ..._mips_converter.py => lp_mps_converter.py} | 0 .../provideq/toolbox/lp/solvers/LpSolver.java | 12 ------- .../MipConfiguration.java} | 32 +++++++++---------- .../toolbox/mip/solvers/MipSolver.java | 12 +++++++ .../solvers/OrToolsCbc.java} | 7 ++-- .../provideq/toolbox/{lp => mip}/simple.mps | 0 7 files changed, 33 insertions(+), 32 deletions(-) rename solvers/ortools/mips/{lp_mips_converter.py => lp_mps_converter.py} (100%) delete mode 100644 src/main/java/edu/kit/provideq/toolbox/lp/solvers/LpSolver.java rename src/main/java/edu/kit/provideq/toolbox/{lp/LpConfiguration.java => mip/MipConfiguration.java} (62%) create mode 100644 src/main/java/edu/kit/provideq/toolbox/mip/solvers/MipSolver.java rename src/main/java/edu/kit/provideq/toolbox/{lp/solvers/OrToolsMipCbc.java => mip/solvers/OrToolsCbc.java} (90%) rename src/main/resources/edu/kit/provideq/toolbox/{lp => mip}/simple.mps (100%) diff --git a/solvers/ortools/mips/cbc.py b/solvers/ortools/mips/cbc.py index ca717a5d..f7e9bdbe 100644 --- a/solvers/ortools/mips/cbc.py +++ b/solvers/ortools/mips/cbc.py @@ -1,7 +1,7 @@ import argparse import os from ortools.linear_solver import pywraplp -from lp_mips_converter import convert_model +from lp_mps_converter import convert_model parser = argparse.ArgumentParser( diff --git a/solvers/ortools/mips/lp_mips_converter.py b/solvers/ortools/mips/lp_mps_converter.py similarity index 100% rename from solvers/ortools/mips/lp_mips_converter.py rename to solvers/ortools/mips/lp_mps_converter.py diff --git a/src/main/java/edu/kit/provideq/toolbox/lp/solvers/LpSolver.java b/src/main/java/edu/kit/provideq/toolbox/lp/solvers/LpSolver.java deleted file mode 100644 index 228936f5..00000000 --- a/src/main/java/edu/kit/provideq/toolbox/lp/solvers/LpSolver.java +++ /dev/null @@ -1,12 +0,0 @@ -package edu.kit.provideq.toolbox.lp.solvers; - -import edu.kit.provideq.toolbox.lp.LpConfiguration; -import edu.kit.provideq.toolbox.meta.ProblemSolver; -import edu.kit.provideq.toolbox.meta.ProblemType; - -public abstract class LpSolver implements ProblemSolver { - @Override - public ProblemType getProblemType() { - return LpConfiguration.LP; - } -} diff --git a/src/main/java/edu/kit/provideq/toolbox/lp/LpConfiguration.java b/src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java similarity index 62% rename from src/main/java/edu/kit/provideq/toolbox/lp/LpConfiguration.java rename to src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java index 465861a2..fcfa966e 100644 --- a/src/main/java/edu/kit/provideq/toolbox/lp/LpConfiguration.java +++ b/src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java @@ -1,8 +1,8 @@ -package edu.kit.provideq.toolbox.lp; +package edu.kit.provideq.toolbox.mip; import edu.kit.provideq.toolbox.ResourceProvider; import edu.kit.provideq.toolbox.exception.MissingExampleException; -import edu.kit.provideq.toolbox.lp.solvers.OrToolsMipCbc; +import edu.kit.provideq.toolbox.mip.solvers.OrToolsCbc; import edu.kit.provideq.toolbox.meta.Problem; import edu.kit.provideq.toolbox.meta.ProblemManager; import edu.kit.provideq.toolbox.meta.ProblemType; @@ -13,30 +13,30 @@ import org.springframework.context.annotation.Configuration; /** - * Definition and registration of the "Linear and Mixed Integer" problems. + * Definition and registration of the "Mixed Integer" problem. */ @Configuration -public class LpConfiguration { +public class MipConfiguration { /** - * LP (Linear Problem) - * A linear optimization problem. - * For a given linear term, + * MIP (Mixed integer Problem) + * A mixed integer optimization problem. + * For a given integer term, * find the minimal variable assignment of the term. */ - public static final ProblemType LP = new ProblemType<>( - "lp", + public static final ProblemType MIP = new ProblemType<>( + "mip", String.class, String.class ); @Bean - ProblemManager getLpManager( - OrToolsMipCbc orToolsMipCbc, + ProblemManager getMipManager( + OrToolsCbc orToolsCbc, ResourceProvider resourceProvider ) { return new ProblemManager<>( - LP, - Set.of(orToolsMipCbc), + MIP, + Set.of(orToolsCbc), loadExampleProblems(resourceProvider) ); } @@ -46,13 +46,13 @@ private Set> loadExampleProblems( try { var problemInputStream = Objects.requireNonNull( getClass().getResourceAsStream("simple.mps"), - "mixed-integer-problem example for LP is unavailable!" + "mixed-integer-problem example is unavailable!" ); - var problem = new Problem<>(LP); + var problem = new Problem<>(MIP); problem.setInput(resourceProvider.readStream(problemInputStream)); return Set.of(problem); } catch (IOException e) { - throw new MissingExampleException(LP, e); + throw new MissingExampleException(MIP, e); } } } diff --git a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/MipSolver.java b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/MipSolver.java new file mode 100644 index 00000000..cd345e70 --- /dev/null +++ b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/MipSolver.java @@ -0,0 +1,12 @@ +package edu.kit.provideq.toolbox.mip.solvers; + +import edu.kit.provideq.toolbox.mip.MipConfiguration; +import edu.kit.provideq.toolbox.meta.ProblemSolver; +import edu.kit.provideq.toolbox.meta.ProblemType; + +public abstract class MipSolver implements ProblemSolver { + @Override + public ProblemType getProblemType() { + return MipConfiguration.MIP; + } +} diff --git a/src/main/java/edu/kit/provideq/toolbox/lp/solvers/OrToolsMipCbc.java b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/OrToolsCbc.java similarity index 90% rename from src/main/java/edu/kit/provideq/toolbox/lp/solvers/OrToolsMipCbc.java rename to src/main/java/edu/kit/provideq/toolbox/mip/solvers/OrToolsCbc.java index 79b09a56..148dff6d 100644 --- a/src/main/java/edu/kit/provideq/toolbox/lp/solvers/OrToolsMipCbc.java +++ b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/OrToolsCbc.java @@ -1,4 +1,4 @@ -package edu.kit.provideq.toolbox.lp.solvers; +package edu.kit.provideq.toolbox.mip.solvers; import edu.kit.provideq.toolbox.Solution; import edu.kit.provideq.toolbox.meta.SolvingProperties; @@ -12,12 +12,12 @@ import reactor.core.publisher.Mono; @Component -public class OrToolsMipCbc extends LpSolver { +public class OrToolsCbc extends MipSolver { private final String scriptPath; private final ApplicationContext context; @Autowired - public OrToolsMipCbc( + public OrToolsCbc( @Value("${ortools.script.ormip}") String scriptPath, ApplicationContext context) { this.scriptPath = scriptPath; @@ -42,6 +42,7 @@ public Mono> solve( ) { var solution = new Solution<>(this); + // TODO: enable both .lp and .mps Problems from the start var processResult = context .getBean(PythonProcessRunner.class, scriptPath) .withArguments( diff --git a/src/main/resources/edu/kit/provideq/toolbox/lp/simple.mps b/src/main/resources/edu/kit/provideq/toolbox/mip/simple.mps similarity index 100% rename from src/main/resources/edu/kit/provideq/toolbox/lp/simple.mps rename to src/main/resources/edu/kit/provideq/toolbox/mip/simple.mps From 9e8dbb3096f21972d4ffcfd4fd661e7cc53c3dbc Mon Sep 17 00:00:00 2001 From: Piotr Malkowski Date: Mon, 10 Mar 2025 12:53:16 +0100 Subject: [PATCH 09/37] feat: add converters --- scripts/mps-lp-to-gms.py | 32 ++++++++++++++++++++++++ solvers/ortools/mips/lp_mps_converter.py | 1 - 2 files changed, 32 insertions(+), 1 deletion(-) create mode 100644 scripts/mps-lp-to-gms.py diff --git a/scripts/mps-lp-to-gms.py b/scripts/mps-lp-to-gms.py new file mode 100644 index 00000000..c8fe8bea --- /dev/null +++ b/scripts/mps-lp-to-gms.py @@ -0,0 +1,32 @@ +import os +import subprocess +import argparse + +parser = argparse.ArgumentParser( + prog="OR-Tools CBC MIP Solver", + description="A CLI Program to solve MIP problems in lp or mips format with Or-Tools CBC" +) + +parser.add_argument("input_file") +parser.add_argument("--output-file") +args = parser.parse_args() + + +input_file_abs = os.path.abspath(args.input_file) +output_file_abs = os.path.abspath(args.output_file) +# Directory where we'll run GAMS so it creates .gms in that same folder +work_dir = os.path.dirname(input_file_abs) +# GAMS by default names the new .gms file after the base of the input file +base_name = os.path.splitext(os.path.basename(input_file_abs))[0] +intermediate_gms = os.path.join(work_dir, base_name + ".gms") +# 1) Run GAMS to convert the LP/MPS -> GMS +# "task=ConvertIn" instructs GAMS to read the file and produce a .gms. +subprocess.run(["gams", input_file_abs, "task=ConvertIn"], check=True, cwd=work_dir) +# 2) Move or rename the newly created .gms to the user-specified output_file +# If you prefer to copy instead, you can use shutil.copyfile(). +if not os.path.exists(intermediate_gms): + raise FileNotFoundError(f"Expected GAMS output file not found: {intermediate_gms}") +os.replace(intermediate_gms, output_file_abs) + + + diff --git a/solvers/ortools/mips/lp_mps_converter.py b/solvers/ortools/mips/lp_mps_converter.py index 2896c21a..6dbea81f 100644 --- a/solvers/ortools/mips/lp_mps_converter.py +++ b/solvers/ortools/mips/lp_mps_converter.py @@ -1,6 +1,5 @@ import cplex - def convert_model(input_file, output_file): """ Loads a model (LP or MPS) into CPLEX and writes it back out From 8cd7955931881ac8bbb9b02bd7729d618483dfc0 Mon Sep 17 00:00:00 2001 From: Piotr Malkowski Date: Wed, 2 Apr 2025 12:11:21 +0200 Subject: [PATCH 10/37] feat: add qubo gams translation --- solvers/gams/mip/LICENSE.txt | 22 + solvers/gams/mip/qubo_solve.gms | 865 ++++++++++++++++++++++++++++++++ 2 files changed, 887 insertions(+) create mode 100644 solvers/gams/mip/LICENSE.txt create mode 100644 solvers/gams/mip/qubo_solve.gms diff --git a/solvers/gams/mip/LICENSE.txt b/solvers/gams/mip/LICENSE.txt new file mode 100644 index 00000000..6308781c --- /dev/null +++ b/solvers/gams/mip/LICENSE.txt @@ -0,0 +1,22 @@ +MIT License + +Copyright (c) 2024 GAMS Software GmbH +Copyright (c) 2024 GAMS Development Corp. + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/solvers/gams/mip/qubo_solve.gms b/solvers/gams/mip/qubo_solve.gms new file mode 100644 index 00000000..cd0ced7f --- /dev/null +++ b/solvers/gams/mip/qubo_solve.gms @@ -0,0 +1,865 @@ +$offEolCom +$offListing +$eolcom # + +$set modelName %1 +$set modelType %2 +$set direction %3 +$set obj %4 +$set penalty %5 +$shift shift shift shift shift + +$set method classic +$set solver cplex +$set maxIter 1 +$set timeLimit 10 +$eval num_threads min(8,numcores) +$set log_on 0 +$set examinerOn 0 +$set get_Q "n" +$set stop_process 0 + +$label ProcessNamedArguments +$ splitOption "%1" key val +$ if x%key%==x $goto ProcessNamedArgumentsDone +$ ifThenI.qubo_solve__arguments %key%==method +$ set method %val% +$ ifE (not(sameas('%method%','classic'))and(not(sameas('%method%','qpu')))) $abort 'Not a valid method. Vaild values are: [classic, qpu].' +$ elseIfI.qubo_solve__arguments %key%==solver +$ set solver %val% +$ elseIfI.qubo_solve__arguments %key%==maxIter +$ set maxIter %val% +$ elseIfI.qubo_solve__arguments %key%==timeLimit +$ set timeLimit %val% +$ elseIfI.qubo_solve__arguments %key%==numThreads +$ set num_threads %val% +$ elseIfI.qubo_solve__arguments %key%==logOn +$ set log_on %val% +$ ifE ((%log_on%<>0)and(%log_on%<>1)and(%log_on%<>2)) $abort 'Not a valid log number. Valid values are [0,1,2].' +$ elseIfI.qubo_solve__arguments %key%==examinerOn +$ set examinerOn %val% +$ IfThenE.examiner_file ((%examinerOn%<>1)and(%examinerOn%<>0)) $abort 'Not a valid examinerOn option. Valid values are [0,1].' +$ elseIfE.examiner_file %examinerOn%==1 $echo examineGamsPoint 1 > examiner.opt +$ endIf.examiner_file +$ elseIfI.qubo_solve__arguments %key%==getQ +$ set get_Q %val% +$ ifE (not(sameas('%get_Q%','y'))and(not(sameas('%get_Q%','n')))) $abort 'Wrong flag for >getQ<. Valid options are [y,n]' +$ eval stop_process sameas('%get_Q%','y') +$ else.qubo_solve__arguments +$ abort Unknown option `%key%`. +$ endif.qubo_solve__arguments +$ shift +$ goTo ProcessNamedArguments +$label ProcessNamedArgumentsDone + +$log *** Options (required):modelName=%modelName%, modelType=%modelType%, objective=%direction%, objectiveVariable=%obj%, penalty=%penalty% +$log *** Options (optional, all): method=%method%, solver=%solver%, maxIter=%maxIter%, timeLimit=%timeLimit%, numThreads=%num_threads%, logOn=%log_on%, examinerOn=%examinerOn%, getQ=%get_Q% + +$onEcho > convert.opt +dumpgdx %modelName%.gdx +GDXQuadratic 1 +$offEcho + +option %modelType%=convert; + +%modelName%.optfile = 1; + +Solve %modelName% use %modelType% %direction% %obj%; # dumping the problem data in a gdx file +put_utility$(%log_on% > 1) 'log' / 'Starting QUBO reformulation.'; + +* QUBO Reformulations +EmbeddedCode Python: +import logging +import re +import warnings +from typing import Tuple, Optional + +import numpy as np +import pandas as pd +from gams import transfer as gt + +log_level_dict = {"0": logging.WARN, "1": logging.INFO, "2": logging.DEBUG} +log_level = log_level_dict.get('%log_on%', logging.WARN) + +if log_level < logging.WARN: + logging.basicConfig(filename='%modelName%_reformulation.log', filemode='w', format='%(message)s', level=log_level, force=True) + +warnings.simplefilter(action='ignore', category=pd.errors.PerformanceWarning) + +is_max = "x" in "%direction%".lower() +gdx_file = r"%modelName%.gdx" +container = gt.Container(gdx_file) +obj_eq_name = container['iobj'].records + +if obj_eq_name is None: + raise Exception("The objective is not defined using a scalar equation. `iobj` in gdx is empty. Quitting.") + +obj_var = container['jobj'].records # fetch the objective variable name +all_vars = container['j'].records # fetches all variable names +raw_a = container['A'].records # A coefficients +eq_data = container['e'].records # fetches equation data + +if raw_a[-raw_a['i'].isin(obj_eq_name['i'].tolist())]['value'].mod(1).sum(axis=0) > 0: # floating point coeffs in objective function is accepted + raise Exception("Reformulation with Non-Integer Coefficients not possible. Quitting.") + +raw_a = raw_a.pivot(index="i", columns="j", values="value").fillna(0) # arranging in a matrix + +if eq_data['lower'].mod(1).sum(axis=0) > 0 or eq_data['upper'].mod(1).sum(axis=0) > 0: + raise Exception("Reformulation with Non-Integer RHS not possible. Quitting.") + +bin_vars = container['jb'].records # fetches binary variable names, if any +int_vars = container['ji'].records # fetches integer variable names, if any +bin_vars = [] if bin_vars is None else bin_vars['j'].to_list() # check if any bin_vars are present +int_vars = [] if int_vars is None else int_vars['j'].to_list() # check if any int_vars are present +obj_var = obj_var['j'].to_list() +all_var_vals = container['x'].records # get all variable values, viz., [level, marginal, lower, upper, scale] + +if len(all_vars) - len(bin_vars) - len(int_vars) != 1: # Continuous variables are not allowed + raise Exception("There are continuous variables. Quitting.") + +obj_eq_name = obj_eq_name['i'].to_list() + +check_quad = container['ANL'].records + +""" +Check if there are any fixed variables in the gdx, i.e., lb=ub=level of any variable. +If such variables exist, separate them from the list of non-fixed vairables and treat them as constanst in the objective function. + +We also need to check if the level of variables are set and handle them separately +""" + +vars_with_lower_bounds = {var.j: var.lower for _, var in all_var_vals.iterrows() if (var.lower > 0) and (var.lower != var.upper)} # would only contain, integer varaibles with lower bound defined +fixed_vars = {var.j: var.level for _, var in all_var_vals.iterrows() if (var.level == var.lower) and (var.level == var.upper)} # check for fixed variables +fixed_and_lower_bounds = {**vars_with_lower_bounds, **fixed_vars} +sum_fixed_obj_var_coeffs = 0 + +if check_quad is not None: + rawquad = container['Q'].records # fetch quadratic terms from the original problem, if any. + + if len(int_vars) != 0: + raise Exception("Quadratic Program with integer variables are not supported.") + + if any(check_quad['j'].isin(fixed_and_lower_bounds.keys())): + raise Exception("Quadratic terms with non-zero variable levels are not supported at the moment.") + +logging.debug("Coefficient matrix: raw_a\n"+raw_a.to_string()) +logging.debug("\nEquation Data: eq_data\n"+eq_data.to_string()) +logging.debug("\nVariable Data: all_var_vals\n"+all_var_vals.to_string()) + +def var_contribution(A: pd.DataFrame, vars: dict, cons: Optional[list] = None) -> np.array: + """ + helper function to calculate the contribution of given variables + in a constraint or set of constraints + + Args: + A: df of coefficients + vars: contributing variables + cons: participating constraints + + Returns: + np.array of Total contribution of all variables for that constraint + """ + cons = slice(None) if cons is None else cons + coeffs_of_vars_in_constraint = A.loc[cons, vars.keys()].to_numpy() + lb_var_levels = np.array(list(vars.values())).reshape((len(vars), 1)) + if coeffs_of_vars_in_constraint.size > 0: + return coeffs_of_vars_in_constraint@lb_var_levels + + return np.array([0]) + +if fixed_and_lower_bounds: # adjust the rhs of equations when level of variables > 0 + logging.info(f"\nList of variables with lower bounds:\n{vars_with_lower_bounds}") + contribution = var_contribution(raw_a, fixed_and_lower_bounds) + eq_data.loc[:,['lower', 'upper']] -= contribution + if fixed_vars: + logging.info(f"\nList of Fixed Variables:\n{fixed_vars}") + # remove the fixed variables from computation + bin_vars = [var for var in bin_vars if var not in fixed_vars] + int_vars = [var for var in int_vars if var not in fixed_vars] + sum_fixed_obj_var_coeffs += np.ndarray.item(var_contribution(raw_a, fixed_vars, cons=obj_eq_name)) + raw_a.drop(fixed_vars, axis=1, inplace=True) # dropping columns from the coefficient matrix + fixed_var_vals = all_var_vals[all_var_vals['j'].isin(fixed_vars)].copy(deep=True) + + logging.debug("\nAfter removing fixed variables and adjusting for non-zero levels: raw_a\n"+raw_a.to_string()) + logging.debug("\nAfter removing fixed variables and adjusting for non-zero levels: eq_data\n"+eq_data.to_string()) + +""" +Check if there exist a row in coefficient matrix with all zero values. This can happen if all vars in a constraint are fixed. +Such row is irrelevant for QUBO and can be dropped out of the matrix and set of constraints. +""" +redundant_cons = list(raw_a[raw_a.apply(abs).sum(axis=1)==0].index) +if len(redundant_cons) > 0: + logging.info(f"\nDropping these redundant constraint: \n{redundant_cons}") + raw_a.drop(redundant_cons, axis=0, inplace=True) + eq_data.drop(eq_data[eq_data['i'].isin(redundant_cons)].index, axis=0, inplace=True) + + +def gen_slacks(var_range: float) -> np.array: + """ + helper function to generate slacks depending on the range of variables or rhs + + Args: + var_range: upper bound of variable + + Returns: + Numpy array containing slack co-efficients + + example: + if var_range=5, then gen_slacks(5) returns [1, 2, 2] + """ + if var_range >= 1e+4: + raise Exception("The Upper bound is greater than or equal to 1e+4, Quitting!") + + power = int(np.log2(var_range)) if var_range > 0 else 0 + bounded_coef = var_range - (2**power - 1) + D_val = [2**i for i in range(power)] + [bounded_coef] + return np.array(D_val) + + +""" +If integer variables exist, convert all integers to binary with '@' as a delimiter of variable names +If Integer variables with lower bound exist, i.e., lb >=1 and lb!=ub, then convert binary variable for that range +If these variables contribute to the objective function, their lower bounds are added as a constant +""" +sum_lower_bound_of_int_vars = 0 +if len(int_vars) != 0: + if vars_with_lower_bounds: + sum_lower_bound_of_int_vars += np.ndarray.item(var_contribution(raw_a, vars_with_lower_bounds, obj_eq_name)) + + int_var_vals = all_var_vals[all_var_vals['j'].isin(int_vars)] + int_to_bin_bounds = {row['j']: gen_slacks(row['upper']-row['lower']) for _,row in int_var_vals.iterrows()} # generate coeffs for converted binary vars + int_bin_vals = pd.DataFrame(columns=['intName', 'binName', 'value']) + for var, bin_bounds in int_to_bin_bounds.items(): + for i in range(len(bin_bounds)): # naming the converted binary vars + new_row = pd.DataFrame({'intName': var, 'binName': f"{var}@_bin{i}", 'value': bin_bounds[i]}, index=[0]) + int_bin_vals = pd.concat([int_bin_vals, new_row], ignore_index=True) if not int_bin_vals.empty else new_row.copy() + + binName_list = int_bin_vals['binName'].to_list() # list of all converted binary variable names + int_bin_name_map = list(int_bin_vals[['intName', 'binName']].itertuples(index=None, name=None)) + logging.info("\nInteger to Binary Mapping: int_bin_vals\n"+int_bin_vals.to_string()) + int_bin_vals = int_bin_vals.pivot(index='intName', columns='binName', values='value').fillna(0) # mapping each binary var to its integer var component + int_bin_vals = int_bin_vals.reindex(labels=int_vars, axis='index') + int_bin_vals = int_bin_vals.reindex(labels=binName_list, axis='columns') + # int_bin_vals.columns = pd.MultiIndex.from_tuples(int_bin_name_map) + + raw_a_int = raw_a[int_vars] + raw_a_int = raw_a_int.dot(int_bin_vals) # updating the "A" coeff matrix with the new coeffs for converted binary vars + raw_a_rest = raw_a[obj_var+bin_vars] + raw_a = pd.concat([raw_a_rest, raw_a_int], axis='columns') # new "A" coeff matrix + logging.info("\nInteger to Binary Mapping: raw_a\n"+raw_a.to_string()) + bin_vars += binName_list # append the list of original binary variables with the list of converted binary variables + +cons = eq_data[-eq_data['i'].isin(obj_eq_name)].reset_index(drop=True) # fetch only the constrainsts and not the objective equation +nvars = len(bin_vars) +nslacks = 0 +obj_var_direction = raw_a[obj_var].loc[obj_eq_name].to_numpy() +obj_var_coeff = raw_a[bin_vars].loc[obj_eq_name].to_numpy() +if obj_var_direction > 0: + obj_var_coeff = -1*obj_var_coeff +obj = np.zeros((nvars, nvars)) +np.fill_diagonal(obj, obj_var_coeff) + +""" +Pre-processing in case special constraints exist +Remove the constraint from the "A" matrix and include the special penalty directly in the objective +Doing so, reduces the number of slack variables used in the final reformulation +special penalty case 1: sum(x_i | 1 <= i <= n) <= 1 => P*sum(x_i*x_j | i < j) +special penalty case 2: x_i + x_j >= 1 => P*(1 - x_i - x_j + x_i*x_j) +""" + +def check_row_entries(df: pd.DataFrame) -> pd.DataFrame: + """ + helper function to filter DataFrame having either 0 or 1 entries in each row. + Args: + df: A Pandas DataFrame. + + Returns: + Filtered DataFrame with rows having either 0 or 1 + """ + row_contains_only_0s_or_1s = df.isin([0, 1]).all(axis=1) + return df[row_contains_only_0s_or_1s] + +# Case 1 implementation +special_cons_case_1_lable = [ele.i for _, ele in cons.iterrows() if ele.upper == 1 and ele.lower != 1] +if special_cons_case_1_lable: + case1_cons = raw_a[bin_vars].loc[special_cons_case_1_lable] + case1_cons = check_row_entries(case1_cons.copy()) + case1_cons_index_lable = list(case1_cons.index) + case1_penalty = case1_cons.to_numpy() + if case1_penalty.size > 0: + case1_penalty = (case1_penalty.T@case1_penalty)/2 + np.fill_diagonal(case1_penalty, np.zeros((1, len(bin_vars)))) + else: # if there are no rows with only 0/1 entries + case1_penalty = np.zeros((nvars, nvars)) + logging.debug(f"\nSpecial constraint case 1:\n{special_cons_case_1_lable}") +else: + case1_cons_index_lable = [] + case1_penalty = np.zeros((nvars, nvars)) + +# Case 2 implementation +special_cons_case_2_lable = [ele.i for _, ele in cons.iterrows() if ele.lower == 1 and ele.upper != 1] +if special_cons_case_2_lable: + case2_cons = raw_a[bin_vars].loc[special_cons_case_2_lable] + case2_cons = check_row_entries(case2_cons.copy()) + case2_cons = case2_cons[case2_cons.sum(axis=1)==2] + case2_cons_index_lable = list(case2_cons.index) + case2_penalty = case2_cons.to_numpy() + if case2_penalty.size > 0: + case2_penalty = (case2_penalty.T@case2_penalty)/2 + case2_diag = np.diag_indices_from(case2_penalty) + case2_penalty[case2_diag] *= -2 + else: # if there are no rows with two 1s in them + case2_penalty = np.zeros((nvars, nvars)) + logging.debug("\nSpecial constraint case 2:\n{special_cons_case_2_lable}") + +else: + case2_cons_index_lable = [] + case2_penalty = np.zeros((nvars, nvars)) + +final_special_cons = case1_cons_index_lable + case2_cons_index_lable +final_special_penalty = case1_penalty + case2_penalty +case2_penalty_offset_factor = len(case2_cons_index_lable) + +P = -1 * %penalty% if is_max else %penalty% # penalty term for classic solvers +P_qpu = %penalty% # penalty term for qpu, since we always going to minimize the qubo using qpu no treatment is required + +if '%method%' == 'classic': + obj += P*final_special_penalty +else: + obj = P_qpu*final_special_penalty - obj if is_max else P_qpu*final_special_penalty + obj # Here, we update the obj and fix the direction to always Minimize for 'qpu' + +cons.drop(cons[cons["i"].isin(final_special_cons)].index, axis=0, inplace=True) +raw_a.drop(final_special_cons, axis=0, inplace=True) + +A_coeff = raw_a.loc[cons['i'], bin_vars] + +quad = None + +def fetch_quadratic_coeff(raw_df: pd.DataFrame) -> np.array: + """ + helper function to convert the original Q matrix of the problem to a symmetric matrix + + Args: + raw_df: Original problem Q data in a pd.DataFrame + + Returns: + Numpy Q matrix + """ + raw_df['value'] /= 2 + mask = raw_df['j_1'].astype(str) == raw_df['j_2'].astype(str) + filtered_quad = raw_df.loc[mask, :].copy() + raw_df = raw_df.loc[~mask].copy() + diag_quad = filtered_quad.reset_index(drop=True) + quad = raw_df.copy(deep=True) + quad['j_1'], quad['j_2'] = raw_df['j_2'], raw_df['j_1'] + quad = pd.concat([raw_df, quad, diag_quad], axis=0) + quad = quad.pivot(index="j_1", columns="j_2", values="value").fillna(0) + quad = quad.reindex(labels=bin_vars, axis='index') + quad = quad.reindex(labels=bin_vars, axis='columns') + return quad.to_numpy() + + +if check_quad is not None: # check if quadratic terms are present in the original problem + logging.debug("\nRaw Q data from GDX: Q\n"+rawquad.to_string()) + rawquad_obj = rawquad[rawquad['i_0'].isin(obj_eq_name)].copy(deep=True) + if len(rawquad_obj.index) != 0: # check if quadratic terms exist in the objective function + rawquad_obj.drop(['i_0'],axis=1,inplace=True) + quad = fetch_quadratic_coeff(raw_df=rawquad_obj) + sum_fixed_obj_var_coeffs /= 2 + + rawquad_cons = rawquad[-rawquad['i_0'].isin(obj_eq_name)] # non-linear constraints without objective equation + if len(rawquad_cons.index) != 0: # non-linear constraints exists + raise Exception("There are non-linear constraints. Quitting.") + + ### Removed the support for quadratic constraints. + # mask = rawquad_cons['j_1'].astype(str) == rawquad_cons['j_2'].astype(str) + # filtered_quad = rawquad_cons.loc[mask, :].reset_index(drop=True) + # if len(filtered_quad.index) == len(rawquad_cons.index): # if pair-wise quadratic terms are not present, i.e., only terms like x1*x1 and not x1*x2 + # filtered_quad.drop(['j_2'], axis=1, inplace=True) + # filtered_quad['value'] /= 2 + # filtered_quad = filtered_quad.pivot(index='i_0', columns='j_1', values='value').fillna(0) + # if len(filtered_quad.columns) != len(bin_vars): + # remain_cols = [var for var in bin_vars if var not in set(filtered_quad.columns)] + # filtered_quad[remain_cols] = 0 + # filtered_quad = filtered_quad[bin_vars] + # A_coeff.loc[filtered_quad.index] += filtered_quad.loc[filtered_quad.index].values + # else: # if pair-wise quadratic terms present, i.e., x1*x2. Quit + # raise Exception("There are non-linear constraints. Quitting.") + +if quad is not None: # add the old quadratic terms/matrix to the new objective + logging.debug("\nUpdate Objective by adding Q: \n"+np.array2string(quad)) + obj += -1*quad if obj_var_direction > 0 else quad + logging.debug("\nNew Q: \n"+np.array2string(obj)) + + +def modify_matrix( + b_vec: np.array, + rhs: float, + slacks: np.array, + A_coeff: pd.DataFrame, + ele: pd.Series, + nslacks: int + ) -> Tuple[np.array, pd.DataFrame, int]: + """ + helper function to update the original "A" matrix of coeffs + + Args: + b_vec: The n*1 vector + rhs : The Right hand side of a constraint + slacks: result of gen_slacks() + A_coeff: "A" matrix + ele: constraint + nslacks: number of slacks + + Returns: + updated b_vec, A_coeff and number of slacks + """ + con_index = ele.i + slack_names = [f'slack_{con_index}_{i}' for i in range(nslacks+1, nslacks + len(slacks) + 1)] + A_coeff[slack_names] = 0 + A_coeff.loc[con_index, slack_names] = slacks + nslacks += len(slacks) + return np.append(b_vec, [rhs]), A_coeff, nslacks + +def get_lhs_bounds(ele: pd.DataFrame) -> Tuple[float, float]: + """ + helper function to find the bounds of a constraint + + Args: + ele: The coefficents of the constraint + + Returns: + lower_bound, upper_bound + """ + return ele[ele < 0].sum(), ele[ele > 0].sum() + +b_vec = np.array([]) +logging.info("\nFinal Cons: \n"+cons.to_string()) +for _, ele in cons.iterrows(): + + if ele.upper == ele.lower: # equal-to type constraint + rhs = ele.lower + lhs_min_lb, lhs_max_ub = get_lhs_bounds(A_coeff.loc[ele.i]) + if (rhs - lhs_min_lb) < 0 or (lhs_max_ub - rhs) < 0: + raise Exception(f"Constraint is infeasible: {ele.i}") + else: + b_vec = np.append(b_vec, [rhs]) + slacks = [] # do not introduce slacks for equality type constraints + + elif ele.upper == np.inf: # greater than type constraint + rhs = ele.lower + _, lhs_max_ub = get_lhs_bounds(A_coeff.loc[ele.i]) + slacks_range = lhs_max_ub - rhs + if slacks_range > 0: + slacks = -1*gen_slacks(slacks_range) + b_vec, A_coeff, nslacks = modify_matrix(b_vec, rhs, slacks, A_coeff, ele, nslacks) + elif slacks_range == 0: + b_vec = np.append(b_vec, [rhs]) + slacks = [] + else: + raise Exception(f"Constraint is infeasible: {ele.i}") + + else: # less-than type constraint + rhs = ele.upper + lhs_min_lb, _ = get_lhs_bounds(A_coeff.loc[ele.i]) + slacks_range = rhs - lhs_min_lb + if slacks_range > 0: + slacks = gen_slacks(slacks_range) + b_vec, A_coeff, nslacks = modify_matrix(b_vec, rhs, slacks, A_coeff, ele, nslacks) + elif slacks_range == 0: + b_vec = np.append(b_vec, [rhs]) + slacks = [] + else: + raise Exception(f"Constraint is infeasible: {ele.i}") + + +logging_a_mat = A_coeff.unstack().reset_index() +logging_a_mat = logging_a_mat[logging_a_mat[0]!= 0] +logging.debug("\nFinal coefficient matrix: \n"+logging_a_mat.to_string()) +logging.debug(f"\nFinal RHS: \n{b_vec}") +logging.debug(f"Constant RHS term: {b_vec.T@b_vec}") +logging.debug(f"Case 2 Offset Penalty Factor: {case2_penalty_offset_factor}") +logging.debug(f"Fixed Variable contribution to Objective Function: {sum_fixed_obj_var_coeffs}") +logging.debug(f"Integer variable lower bound contribution: {sum_lower_bound_of_int_vars}") + +# A matrix and b_vec are available. Now, for penalization: $(A.X - B)^{2}$ = $(A.X - B)^{T} * (A.X - B)$ +a_mat = A_coeff.to_numpy() +nvars += nslacks # increment the total number of variables by total number of slack variable used +X_diag = np.zeros((nvars, nvars)) +np.fill_diagonal(X_diag,-b_vec.T@a_mat) +new_x = a_mat.T@a_mat + 2*X_diag + +newobj = np.zeros((nvars,nvars)) +newobj[:len(bin_vars), :len(bin_vars)] = obj # define the new objective: Q for the qubo + + +def qubo_to_ising(Q: dict, offset: float = 0.0) -> Tuple[dict, dict, float]: + """ + This is the Qubo to Ising Reformulation. Here, the variable X in {-1,1} + + Args: + Q: in a form of dict, {(i,j): val} + offset: offset from the Qubo reformulation + + Returns: + h: the bias vector + J: the coupling matrix + offset: adjusted offset for the Ising model + """ + + h = {} + J = {} + linear_offset = 0.0 + quadratic_offset = 0.0 + + for (u, v), bias in Q.items(): + if u == v: + if u in h: + h[u] += .5 * bias + else: + h[u] = .5 * bias + linear_offset += bias + else: + if bias != 0.0: + J[(u, v)] = .25 * bias + + if u in h: + h[u] += .25 * bias + else: + h[u] = .25 * bias + + if v in h: + h[v] += .25 * bias + else: + h[v] = .25 * bias + + quadratic_offset += bias + + offset += .5 * linear_offset + .25 * quadratic_offset + + return h, J, offset + + +def qubo_to_maxcut(Q: np.array) -> np.array: + """ + This is the Qubo to Maxcut Reformulation. This can be used for SDP procedures. + + Args: + Q: a n x n symmetric numpy matrix + + Returns: + n x 1 vector associated with the extra variable required in max cut transformation + """ + + return -1*np.sum(Q, axis=1) + +### Section to get the qubo for submitting it to the dwave-hybrid method +if '%method%' == 'qpu': + if '%get_Q%'.lower() == 'y': + raise Exception(f"Specify >method=classic< to export Q matrix in CSV files.") + Q = P_qpu*new_x + newobj # Note: We are always minimizing in the case of qpu and we have already fixed the direction of obj on line 184 + offset = P_qpu*b_vec.T@b_vec + P_qpu*case2_penalty_offset_factor + logging.debug(f"\nPenaly: {P_qpu} | Total Offset: {offset}") + if Q.size > 0: + Qdf = pd.DataFrame(Q, columns=list(A_coeff.columns), index=list(A_coeff.columns)) + Qdf = Qdf.unstack() + Qdf = Qdf.reset_index() + Qdf_dict = {(row['level_0'], row['level_1']): row[0] for _,row in Qdf.iterrows()} + else: + Qdf_dict = {} + +### Section to solve the qubo using miqcp +elif '%method%' in ['classic', 'sdp']: + const = P*b_vec.T@b_vec + P*case2_penalty_offset_factor + sum_fixed_obj_var_coeffs + sum_lower_bound_of_int_vars + logging.debug(f"\nPenalty: {P} | Total Offset: {const}\n") + Q = newobj + P*new_x + continue_forward = True + #sparsity = 1.0 - (np.count_nonzero(Q) / float(Q.size)) + if '%get_Q%'.lower() == 'y': + np.savetxt(fr'%modelName%_p{P:.0f}_c{const:.0f}.csv', Q, delimiter=",") + logging.debug(f"\nExported Q matrix as >%modelName%_p{P:.0f}_c{const:.0f}.csv<\n") + continue_forward = False + ### For producing the qs format for QUBOWL + # non_zero_indices = np.tril_indices_from(Q) + # non_zero_values = Q[non_zero_indices] + # with open('QMat_%modelName%.qs', 'w') as f: + # f.write(f"{Q.shape[0]} {len(non_zero_values)} {const}\n") + # for i, j, value in zip(*non_zero_indices, non_zero_values): + # if value != 0: + # f.write(f"{i+1} {j+1} {value}\n") + if continue_forward: + Qdf = pd.DataFrame(Q, columns=list(A_coeff.columns), index=list(A_coeff.columns)) + Qdf = Qdf.unstack() + Qdf = Qdf.reset_index() + Qdf = list(Qdf.itertuples(index=None, name=None)) + + m = gt.Container() + qconst = gt.Parameter(m, "qconst", None, const, description="Constant term to offset the objective value to original") + if Qdf: + qi = gt.Set(m, "qi", records=A_coeff.columns, description="QUBO variables") + qd = gt.Parameter(m, "qd", [qi, qi], records=Qdf, description="Q matrix") + if '%method%' == 'sdp': + max_cut_var = qubo_to_maxcut(Q) + gt.Parameter(m, "c_sun", [qi], records=list(zip(A_coeff.columns, max_cut_var)), description="extra variable for max cut transformation, ") + else: + max_cut_var = np.zeros((len(A_coeff.columns),1)) + gt.Parameter(m, "c_sun", [qi], records=list(zip(A_coeff.columns, max_cut_var)), description="extra variable for max cut transformation, all set to zero since SDP method is not chosen.") + else: + qi = gt.Set(m, "qi", records=['all_variables_fixed'], description="QUBO variables") + qd = gt.Parameter(m, "qd", [qi, qi], records=[('all_variables_fixed', 'all_variables_fixed', 0)], description="Q matrix") + if '%method%' == 'sdp': + raise Exception("All variables are fixed. Q matrix is Empty. Quitting SDP SOLVE.") + m.write('qout_%modelName%.gdx') +endEmbeddedCode + +abort$execError 'Error occured in Reformulation!'; +put_utility$(%log_on% > 1) 'log' / 'Reformulation complete. The Q matrix should be available in qout_%modelName%.gdx'; +abort.noError$[%stop_process%] '--getQ=Y. Skipping rest of the code.'; + +* Solve the QUBO using a miqcp +$ifThenE.method sameas('%method%','classic') +embeddedCode GAMS: lo=%gams.lo% +set qi; + +Alias (qi,i,j); + +Parameter qd(i,j); + +Scalar qconst; + +Parameter c_sun(i); + +$gdxload 'qout_%modelName%.gdx' qconst qi qd c_sun + +Equation defQubo; + +Binary variable x(i); + +Free variable z; + +defQubo.. z =E= sum((i,j), x(i)*qd(i,j)*x(j)) + qconst; + +Model qubo /all/; + +option miqcp=%solver%, threads=%num_threads%; + +qubo.reslim = %timeLimit%; + +put_utility$(%log_on% > 1) 'log' / 'Solving QUBO classically using -solver=%solver%'; +Solve qubo %direction% z using miqcp; + +execute_unload 'qout_%modelName%.gdx'; +endembeddedCode + +* Map the QUBO's solution to the original problem +EmbeddedCode Python: +q_container = gt.Container('qout_%modelName%.gdx') +obj_var_coeff = q_container['z'].records +obj_var = container['jobj'].records['j'].values[0] + +all_vars = container['j'].records +original_obj_sym = all_vars[all_vars['uni'] == obj_var]['element_text'].values[0] +rem_syms = all_vars[all_vars['uni'] != obj_var]['uni'].to_list() +optimized_vals = q_container['x'].records + +if fixed_vars: + fixed_var_vals.rename({'j': 'i'}, axis=1, inplace=True) + optimized_vals = pd.concat([optimized_vals, fixed_var_vals], ignore_index=True) + +if len(int_vars) != 0: # check if integer variable exist. If yes, combine and merge the solution of converted binary variables to their integer representation + int_bin_vals_unstack = int_bin_vals.unstack().reset_index() + int_bin_vals_unstack.drop(int_bin_vals_unstack[int_bin_vals_unstack[0] == 0].index, inplace=True) + bin_to_int_vals = pd.merge(int_bin_vals_unstack, optimized_vals, how="left", left_on="binName", right_on="i") + bin_to_int_vals['final_level'] = bin_to_int_vals[0]*bin_to_int_vals['level'] + bin_to_int_vals = bin_to_int_vals.groupby('intName')['final_level'].sum().reset_index() + + original_int_vals = container['x'].records + original_int_vals = original_int_vals[original_int_vals['j'].isin(bin_to_int_vals['intName'])].copy(deep=True) + original_int_vals = pd.merge(original_int_vals, bin_to_int_vals, left_on="j", right_on="intName", how="left") + original_int_vals.drop(['level', 'intName'], axis=1, inplace=True) + original_int_vals.rename({'j':'i', 'final_level': 'level'}, axis=1, inplace=True) + original_int_vals = original_int_vals[["i", "level", "marginal", "lower", "upper", "scale"]] + + optimized_vals.drop(optimized_vals[optimized_vals.i.isin(binName_list)].index, inplace=True) + optimized_vals = pd.concat([optimized_vals, original_int_vals], ignore_index=True) + if vars_with_lower_bounds: + optimized_vals.loc[optimized_vals["i"].isin(vars_with_lower_bounds.keys()), 'level'] += list(vars_with_lower_bounds.values()) + + +mapper = {} +for _, ele in all_vars.iterrows(): # extract the domain and symbols from the gdx + domain = re.findall(r'\((.*?)\)', ele['element_text']) + var_sym = re.findall(r'(.+)(?=\()', ele['element_text']) + if len(domain) > 0: + domain = domain[0].strip(r"\'") + if var_sym[0] not in mapper: + mapper[var_sym[0]] = {ele['uni'] : domain} + else: + mapper[var_sym[0]][ele['uni']] = domain + +for vars, ele in mapper.items(): + newsol = optimized_vals.copy(deep=True) + newsol = newsol[newsol['i'].isin(ele.keys())] + newsol.index = newsol['i'] + newsol.drop(['i'], axis=1, inplace=True) + newsol['oldlabel'] = None + + for key, val in ele.items(): + newsol.loc[key, 'oldlabel'] = val + newsol = newsol[['oldlabel', 'level', 'marginal', 'lower','upper','scale']] + split_labels = newsol['oldlabel'].str.split(',', expand=True) + new_columns = [f'newlabel{i+1}' for i in range(split_labels.shape[1])] + split_labels.columns = new_columns + newsol = pd.concat([split_labels, newsol], axis=1) + newsol.drop(['oldlabel'], axis=1, inplace=True) + gams.set(vars, list(newsol.itertuples(index=None,name=None))) + +if len(mapper) == 0: # support for flat gms file + opt_list = optimized_vals.copy(deep=True) + opt_list.drop(['i'], axis=1, inplace=True) + flat_var_mapping = all_vars.set_index('uni')['element_text'].to_dict() + for key, val in flat_var_mapping.items(): + if val != original_obj_sym: + gams.set(val, list(opt_list.iloc[optimized_vals[optimized_vals['i']==key].index].itertuples(index=None,name=None))) + +""" +The code below is required to calculate the contribution of variables towards the objective using the new levels obtained from the QUBO solve. +Since the QUBO solve returns a different level for the objective variable when the optimal solution is not returned, for example, it includes the penalty for every constraint not satisfied. +""" +x_l = optimized_vals[optimized_vals['i'].isin(rem_syms)]['level'].to_numpy() +orig_syms_w_new_levels = optimized_vals[optimized_vals['i'].isin(rem_syms)].set_index('i')['level'].to_dict() +quad_contribution = x_l.T@quad@x_l if check_quad is not None else 0 #at the moment `quad` only contains contribution from the objective row. +orig_jacobian = container['A'].records # A coefficients +orig_jacobian = orig_jacobian.pivot(index="i", columns="j", values="value").fillna(0) # arranging in a matrix +linear_contribution = var_contribution(orig_jacobian, orig_syms_w_new_levels, cons=obj_eq_name).flatten()[0] +total_objective_contribution = linear_contribution + quad_contribution +total_objective_contribution = -1*total_objective_contribution if obj_var_direction > 0 else total_objective_contribution +obj_var_coeff.loc[0, 'level'] = total_objective_contribution + +gams.set(original_obj_sym, list(obj_var_coeff.itertuples(index=None, name=None))) +endEmbeddedCode + +* Solve the QUBO on the QPU +$elseIfE.method sameas('%method%','qpu') +embeddedCode Python: +try: + from dwave.system import LeapHybridSampler +except: + raise Exception("Package 'Dwave.system' is not available") + +try: + from dimod import BinaryQuadraticModel +except: + raise Exception("Package 'Dimod' is not available") + +gams.printLog("Submitting problem to QPU") +solver = LeapHybridSampler.default_solver +sampler = LeapHybridSampler(solver=solver) +from datetime import datetime + +cur_date = datetime.now().strftime("%d.%m.%Y") + +if Q.size > 0:# check if QMatrix exist, if not skip model generation + bqm = BinaryQuadraticModel.from_qubo(Qdf_dict, offset=offset) + max_iter = %maxIter% + all_sol = pd.DataFrame() + sampler_label_base = f'QUBO-%modelName%-{cur_date}' + + for counter in range(max_iter): # For repeated sampling, a default option is not avaiable for the chosen sampler + gams.printLog(f"Solve number: {counter+1}") + answer = sampler.sample(bqm, label=f'{sampler_label_base}-{counter+1}', time_limit=%timeLimit%) + new_df = answer.to_pandas_dataframe() + all_sol = pd.concat([all_sol, new_df], ignore_index=True) + + best_val = all_sol['energy'].min() + best_sol = all_sol[all_sol['energy']==best_val].reset_index(drop=True) + best_sol.drop(['energy', 'num_occurrences'], axis=1, inplace=True) + best_sol = best_sol.loc[0].to_dict() + + gams.printLog(f"Best Objective Value: {best_val}") + gams.printLog(f"All solutions:\n{all_sol}") + + sample = pd.DataFrame(best_sol.items(), columns=['j', 'level']) +else: + warnings.warn("QMatrix size is zero. Skipping model creation.") + sample = pd.DataFrame(columns=['j', 'level']) + best_val = 0 + +if fixed_vars: + sample = pd.concat([sample, fixed_var_vals[['j', 'level']]], ignore_index=True) + +if len(int_vars) != 0: # check if integer variable exist. If yes, combine and merge the solution of converted binary variables to their integer representation + int_bin_vals_unstack = int_bin_vals.unstack().reset_index() + int_bin_vals_unstack.drop(int_bin_vals_unstack[int_bin_vals_unstack[0] == 0].index, inplace=True) + + bin_to_int_vals = pd.merge(int_bin_vals_unstack, sample, how="left", left_on="binName", right_on="j") + bin_to_int_vals['final_level'] = bin_to_int_vals[0]*bin_to_int_vals['level'] + bin_to_int_vals = bin_to_int_vals.groupby('intName')['final_level'].sum().reset_index() + + original_int_vals = container['x'].records + original_int_vals = original_int_vals[original_int_vals['j'].isin(bin_to_int_vals['intName'])].copy(deep=True) + original_int_vals = pd.merge(original_int_vals, bin_to_int_vals, left_on="j", right_on="intName", how="left") + original_int_vals.drop(['level', 'intName', 'marginal','lower','upper','scale'], axis=1, inplace=True) + original_int_vals.rename({'final_level': 'level'}, axis=1, inplace=True) + + if vars_with_lower_bounds: + original_int_vals.loc[original_int_vals["j"].isin(vars_with_lower_bounds.keys()), 'level'] += list(vars_with_lower_bounds.values()) + + sample.drop(sample[sample.j.isin(binName_list)].index, inplace=True) + sample = pd.concat([sample, original_int_vals], ignore_index=True) + +def map_vars(container: gt.Container, solution: pd.DataFrame, obj_val: float) -> None: + + """ + helper function to map the solution returned from the qpu to the original problem and set the respective gams symbols + """ + + oldvars = container['x'].records + oldvars.drop(['level'], inplace=True, axis=1) + + res = solution.merge(oldvars, how='right', on='j') + vardict = container['j'].records + separate_sym_domain = vardict['element_text'].str.split('(', expand=True) + + if len(separate_sym_domain.columns) == 1: # check if all variables are flat, i.e., no domain + vardict['symbol'] = separate_sym_domain[0] + vardict['domain'] = None + else: + vardict[['symbol','domain']] = separate_sym_domain[[0,1]] + vardict['domain'] = vardict['domain'].str.rstrip(')') + vardict['domain'] = vardict['domain'].str.strip(r"\'") + + vardict.drop(columns=['element_text'], inplace=True) + vardict.rename({"uni": "j"}, axis=1, inplace=True) + + final = res.merge(vardict, how='right', on='j') + final.loc[0, ['level', 'marginal']].to_list() + final.loc[final['symbol']=='%obj%', 'level'] = obj_val + + for symbol in final['symbol'].unique(): + new_vals = [] + temp = final[final['symbol']==symbol].reset_index(drop=True) + temp = temp[['domain','level', 'marginal', 'lower', 'upper', 'scale']] + for i in range(len(temp)): + domain = temp.loc[i, 'domain'].split(",") if type(temp.loc[i, 'domain']) == str else [] + new_vals += [(*domain, *temp.loc[i, ['level', 'marginal', 'lower', 'upper', 'scale']].to_list())] + gams.set(symbol, new_vals, mergeType=MergeType.REPLACE) + +if is_max: + best_val = (sum_fixed_obj_var_coeffs + sum_lower_bound_of_int_vars) - best_val +else: + best_val += (sum_fixed_obj_var_coeffs + sum_lower_bound_of_int_vars) + +map_vars(container, sample, best_val) +endEmbeddedCode +$endIf.method + +$ifThenE.run_examiner %examinerOn% + +option solver=examiner; +%modelName%.optfile = 1; +solve %modelName% %direction% %obj% use %modelType%; + +$endIf.run_examiner + +$onListing \ No newline at end of file From e0b244268004cd44280681d1d79fa1ac37f8bfb2 Mon Sep 17 00:00:00 2001 From: Piotr Malkowski Date: Wed, 2 Apr 2025 12:11:37 +0200 Subject: [PATCH 11/37] feat: add lp testcase --- .../edu/kit/provideq/toolbox/mip/simple.lp | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) create mode 100644 src/main/resources/edu/kit/provideq/toolbox/mip/simple.lp diff --git a/src/main/resources/edu/kit/provideq/toolbox/mip/simple.lp b/src/main/resources/edu/kit/provideq/toolbox/mip/simple.lp new file mode 100644 index 00000000..89e74067 --- /dev/null +++ b/src/main/resources/edu/kit/provideq/toolbox/mip/simple.lp @@ -0,0 +1,17 @@ +Maximize + obj: 3 x + 4 y + +Subject To + c1: x + 2 y <= 14 + c2: 3 x - y >= 0 + c3: x - y <= 2 + +Bounds + x >= 0 + y >= 0 + +Generals + x + y + +End From 757166bc3ace32fce1ac22eef2dc2baecf44eebc Mon Sep 17 00:00:00 2001 From: Piotr Malkowski Date: Wed, 2 Apr 2025 12:12:03 +0200 Subject: [PATCH 12/37] feat: add qubo mip solver --- solvers/gams/mip/mip.gms | 96 +++++++++++++++++++ .../toolbox/mip/solvers/QuboMipSolver.java | 58 +++++++++++ 2 files changed, 154 insertions(+) create mode 100644 solvers/gams/mip/mip.gms create mode 100644 src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java diff --git a/solvers/gams/mip/mip.gms b/solvers/gams/mip/mip.gms new file mode 100644 index 00000000..cb811c0a --- /dev/null +++ b/solvers/gams/mip/mip.gms @@ -0,0 +1,96 @@ +*----------------------------------------------------------- +* Cross-Platform MPS -> GAMS -> QUBO in One Script +* Usage: +* gams convert_qb.gms --INPUT=myModel.mps +* or +* gams convert_qb.gms --INPUT=myModel.lp (if mps2gms supports LP) +*----------------------------------------------------------- + +$if not set INPUT $set INPUT myModel.mps + +* Turn the GAMS macro %INPUT% into a real OS environment variable GAMSINP: +$setEnv GAMSINP "%INPUT%" + + +$onEmbeddedCode Python: +import os +import sys +import subprocess +import re + +# 1) Get the input file from GAMS macro or default +inp_file = os.environ.get('GAMSINP', 'myModel.mps') + +# 2) Remove old "gams.gms" if it exists +if os.path.exists('gams.gms'): + os.remove('gams.gms') + +# 3) Invoke mps2gms on the input file => "gams.gms" +if not os.path.exists(inp_file): + sys.exit(f"*** ERROR: input file not found: {inp_file}") + +cmd = ['mps2gms', inp_file, 'gams.gms'] +print(f"--- Running: {' '.join(cmd)}") +res = subprocess.run( + cmd, + stdout=subprocess.PIPE, + stderr=subprocess.PIPE +) +out = res.stdout.decode('utf-8', errors='ignore') +err = res.stderr.decode('utf-8', errors='ignore') + + +# 4) Rename gams.gms => convertedModel.gms +if os.path.exists('convertedModel.gms'): + os.remove('convertedModel.gms') + +if not os.path.exists('gams.gms'): + sys.exit('*** ERROR: "gams.gms" not found - did mps2gms succeed?') + +os.rename('gams.gms', 'convertedModel.gms') + +# 5) Parse "convertedModel.gms" for Solve statement +new_name = 'convertedModel.gms' +with open(new_name, 'r') as f: + content = f.read() + +pattern = re.compile( + r'Solve\s+(\S+)\s+using\s+(\S+)\s+(maximizing|minimizing)\s+(\S+)\s*;', + re.IGNORECASE +) +match = pattern.search(content) +if not match: + sys.exit('*** ERROR: No "Solve ... using ..." statement found in convertedModel.gms') + +model_name = match.group(1) # e.g. "m" +model_type = match.group(2) # e.g. "LP" or "MIP" +sense_text = match.group(3) # "maximizing" or "minimizing" +obj_var = match.group(4) # e.g. "x7" + +if sense_text.lower().startswith('max'): + gams_sense = 'max' +else: + gams_sense = 'min' + +# 6) Show info to GAMS listing: +print(f'--- Detected from "{new_name}":') +print(f' Model name = {model_name}') +print(f' Model type = {model_type}') +print(f' Sense = {gams_sense}') +print(f' Objective var = {obj_var}') + +# 7) Emit GAMS macros so the main code can pick them up +print(f'$set PARSED_MODELNAME {model_name}') +print(f'$set PARSED_MODELTYPE {model_type}') +print(f'$set PARSED_SENSE {gams_sense}') +print(f'$set PARSED_OBJVAR {obj_var}') + +# Also define a hardcoded penalty (example: 10) +print('$set HARDCODED_PENALTY 10') +$offEmbeddedCode + +*----------------------------------------------------------- +* 8) Finally, call qubo_solve.gms with the macros we just set +*----------------------------------------------------------- +$batinclude qubo_solve.gms %PARSED_MODELNAME% %PARSED_MODELTYPE% \ + %PARSED_SENSE% %PARSED_OBJVAR% %HARDCODED_PENALTY% -logOn=2 diff --git a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java new file mode 100644 index 00000000..6d820c28 --- /dev/null +++ b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java @@ -0,0 +1,58 @@ +package edu.kit.provideq.toolbox.mip.solvers; + +import edu.kit.provideq.toolbox.Solution; +import edu.kit.provideq.toolbox.meta.SolvingProperties; +import edu.kit.provideq.toolbox.meta.SubRoutineResolver; +import edu.kit.provideq.toolbox.process.GamsProcessRunner; +import edu.kit.provideq.toolbox.process.ProcessResult; +import edu.kit.provideq.toolbox.process.ProcessRunner; +import org.springframework.beans.factory.annotation.Autowired; +import org.springframework.beans.factory.annotation.Value; +import org.springframework.context.ApplicationContext; +import org.springframework.stereotype.Component; +import reactor.core.publisher.Mono; + +@Component +public class QuboMipSolver extends MipSolver { + private final String translationScriptPath; + private final ApplicationContext context; + + @Autowired + public QuboMipSolver( + @Value("${ortools.script.ormip}") String translationScriptPath, + ApplicationContext context) { + this.translationScriptPath = translationScriptPath; + this.context = context; + } + + @Override + public String getName() { + return "TODO"; + } + + @Override + public String getDescription() { + return "TODO"; + } + + @Override + public Mono> solve( + String input, + SubRoutineResolver subRoutineResolver, + SolvingProperties properties + ) { + var solution = new Solution<>(this); + + // Run GAMS via console + ProcessResult processResult = context + .getBean(GamsProcessRunner.class, translationScriptPath) + .withArguments( + "--INPUT=" + ProcessRunner.INPUT_FILE_PATH, + "--SOLOUTPUT=" + ProcessRunner.OUTPUT_FILE_PATH + ) + .writeInputFile(input, "problem.lp") + .readOutputFile() + .run(getProblemType(), solution.getId()); + return Mono.just(processResult.applyTo(solution)); + } +} From d29c94715a13149ee153c59d017e8cb2f0fc41a4 Mon Sep 17 00:00:00 2001 From: "piotr.malkowski" Date: Wed, 9 Apr 2025 09:14:23 +0200 Subject: [PATCH 13/37] feat: adapt gams parsing logic --- solvers/gams/mip/lp-to-gms.gms | 64 ++++++++++++++++ solvers/gams/mip/mip.gms | 74 +++++-------------- .../toolbox/mip/solvers/QuboMipSolver.java | 13 +++- .../toolbox/process/GamsProcessRunner.java | 2 +- .../toolbox/process/PythonProcessRunner.java | 3 +- src/main/resources/application.properties | 3 + .../edu/kit/provideq/toolbox/mip/simple.lp | 16 ++-- 7 files changed, 106 insertions(+), 69 deletions(-) create mode 100644 solvers/gams/mip/lp-to-gms.gms diff --git a/solvers/gams/mip/lp-to-gms.gms b/solvers/gams/mip/lp-to-gms.gms new file mode 100644 index 00000000..de4f8e15 --- /dev/null +++ b/solvers/gams/mip/lp-to-gms.gms @@ -0,0 +1,64 @@ +$setEnv GAMSINP "%INPUT%" + +$onEmbeddedCode Python: +import os +import sys +import subprocess +import re + +inp_file = os.environ.get('GAMSINP') +print(f"{inp_file}") + +if not os.path.exists(inp_file): + sys.exit(f"*** ERROR: input file not found: {inp_file}") + +# Replace .lp with .gms +base, ext = os.path.splitext(inp_file) +new_name = base + '.gms' + +cmd = ['mps2gms', inp_file] +print(f"--- Running: {' '.join(cmd)}") +res = subprocess.run( + cmd, + stdout=subprocess.PIPE, + stderr=subprocess.PIPE +) +#out = res.stdout.decode('utf-8', errors='ignore') +#err = res.stderr.decode('utf-8', errors='ignore') + +with open(new_name, 'r', encoding='utf-8', errors='ignore') as f: + content = f.read() + +pattern = re.compile( + r'Solve\s+(\S+)\s+using\s+(\S+)\s+(maximizing|minimizing)\s+(\S+)\s*;', + re.IGNORECASE +) +match = pattern.search(content) +if not match: + sys.exit('*** ERROR: No "Solve ... using ..." statement found in generated .gms file') + +model_name = match.group(1) # e.g. "m" +model_type = match.group(2) # e.g. "LP" or "MIP" +sense_text = match.group(3) # "maximizing" or "minimizing" +obj_var = match.group(4) # e.g. "x7" + +if sense_text.lower().startswith('max'): + gams_sense = 'max' +else: + gams_sense = 'min' + +# 6) Show info to GAMS listing: +print(f'--- Detected from "{new_name}":') +print(f' Model name = {model_name}') +print(f' Model type = {model_type}') +print(f' Sense = {gams_sense}') +print(f' Objective var = {obj_var}') + +with open("parsed_macros.gms", "w") as macrofile: + macrofile.write(f"$set PARSED_MODELNAME {model_name}\n") + macrofile.write(f"$set PARSED_MODELTYPE {model_type}\n") + macrofile.write(f"$set PARSED_SENSE {gams_sense}\n") + macrofile.write(f"$set PARSED_OBJVAR {obj_var}\n") + macrofile.write(f"$set HARDCODED_PENALTY 10\n") + +$offEmbeddedCode diff --git a/solvers/gams/mip/mip.gms b/solvers/gams/mip/mip.gms index cb811c0a..7dfa68a6 100644 --- a/solvers/gams/mip/mip.gms +++ b/solvers/gams/mip/mip.gms @@ -1,16 +1,5 @@ -*----------------------------------------------------------- -* Cross-Platform MPS -> GAMS -> QUBO in One Script -* Usage: -* gams convert_qb.gms --INPUT=myModel.mps -* or -* gams convert_qb.gms --INPUT=myModel.lp (if mps2gms supports LP) -*----------------------------------------------------------- - -$if not set INPUT $set INPUT myModel.mps - -* Turn the GAMS macro %INPUT% into a real OS environment variable GAMSINP: $setEnv GAMSINP "%INPUT%" - +$setEnv GAMSPEN "%PENALTY%" $onEmbeddedCode Python: import os @@ -18,40 +7,25 @@ import sys import subprocess import re -# 1) Get the input file from GAMS macro or default -inp_file = os.environ.get('GAMSINP', 'myModel.mps') +inp_file = os.environ.get('GAMSINP') +print(f"{inp_file}") -# 2) Remove old "gams.gms" if it exists -if os.path.exists('gams.gms'): - os.remove('gams.gms') - -# 3) Invoke mps2gms on the input file => "gams.gms" if not os.path.exists(inp_file): sys.exit(f"*** ERROR: input file not found: {inp_file}") -cmd = ['mps2gms', inp_file, 'gams.gms'] +# Replace .lp with .gms +base, ext = os.path.splitext(inp_file) +new_name = base + '.gms' + +cmd = ['mps2gms', inp_file] print(f"--- Running: {' '.join(cmd)}") res = subprocess.run( cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE ) -out = res.stdout.decode('utf-8', errors='ignore') -err = res.stderr.decode('utf-8', errors='ignore') - -# 4) Rename gams.gms => convertedModel.gms -if os.path.exists('convertedModel.gms'): - os.remove('convertedModel.gms') - -if not os.path.exists('gams.gms'): - sys.exit('*** ERROR: "gams.gms" not found - did mps2gms succeed?') - -os.rename('gams.gms', 'convertedModel.gms') - -# 5) Parse "convertedModel.gms" for Solve statement -new_name = 'convertedModel.gms' -with open(new_name, 'r') as f: +with open(new_name, 'r', encoding='utf-8', errors='ignore') as f: content = f.read() pattern = re.compile( @@ -60,7 +34,7 @@ pattern = re.compile( ) match = pattern.search(content) if not match: - sys.exit('*** ERROR: No "Solve ... using ..." statement found in convertedModel.gms') + sys.exit('*** ERROR: No "Solve ... using ..." statement found in generated .gms file') model_name = match.group(1) # e.g. "m" model_type = match.group(2) # e.g. "LP" or "MIP" @@ -72,25 +46,15 @@ if sense_text.lower().startswith('max'): else: gams_sense = 'min' -# 6) Show info to GAMS listing: -print(f'--- Detected from "{new_name}":') -print(f' Model name = {model_name}') -print(f' Model type = {model_type}') -print(f' Sense = {gams_sense}') -print(f' Objective var = {obj_var}') +# Create the replacement line using the matched groups +penalty = os.environ.get("GAMSPEN", "100") # Default to 100 if not set +new_line = f'$batinclude %IDIR% {model_name} {model_type} {gams_sense} {obj_var} {penalty} -logOn=2;' -# 7) Emit GAMS macros so the main code can pick them up -print(f'$set PARSED_MODELNAME {model_name}') -print(f'$set PARSED_MODELTYPE {model_type}') -print(f'$set PARSED_SENSE {gams_sense}') -print(f'$set PARSED_OBJVAR {obj_var}') +# Replace the matching "Solve ..." line with the new line (only the first occurrence) +new_content = pattern.sub(new_line, content, count=1) -# Also define a hardcoded penalty (example: 10) -print('$set HARDCODED_PENALTY 10') -$offEmbeddedCode +# Write the modified content back to the .gms file +with open(new_name, 'w', encoding='utf-8', errors='ignore') as f: + f.write(new_content) -*----------------------------------------------------------- -* 8) Finally, call qubo_solve.gms with the macros we just set -*----------------------------------------------------------- -$batinclude qubo_solve.gms %PARSED_MODELNAME% %PARSED_MODELTYPE% \ - %PARSED_SENSE% %PARSED_OBJVAR% %HARDCODED_PENALTY% -logOn=2 +$offEmbeddedCode diff --git a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java index 6d820c28..b300f40d 100644 --- a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java +++ b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java @@ -14,14 +14,18 @@ @Component public class QuboMipSolver extends MipSolver { - private final String translationScriptPath; + private final String scriptPath; + private final String dependencyScriptPath; private final ApplicationContext context; + public static final int PENALTY = 2; @Autowired public QuboMipSolver( - @Value("${ortools.script.ormip}") String translationScriptPath, + @Value("${gams.script.mip}") String scriptPath, + @Value("${gams.script.qubo_solve}") String dependencyScriptPath, ApplicationContext context) { - this.translationScriptPath = translationScriptPath; + this.scriptPath = scriptPath; + this.dependencyScriptPath = dependencyScriptPath; this.context = context; } @@ -45,9 +49,10 @@ public Mono> solve( // Run GAMS via console ProcessResult processResult = context - .getBean(GamsProcessRunner.class, translationScriptPath) + .getBean(GamsProcessRunner.class, scriptPath) .withArguments( "--INPUT=" + ProcessRunner.INPUT_FILE_PATH, + "--IDIR=" + dependencyScriptPath, "--SOLOUTPUT=" + ProcessRunner.OUTPUT_FILE_PATH ) .writeInputFile(input, "problem.lp") diff --git a/src/main/java/edu/kit/provideq/toolbox/process/GamsProcessRunner.java b/src/main/java/edu/kit/provideq/toolbox/process/GamsProcessRunner.java index a00921b3..52f33b96 100644 --- a/src/main/java/edu/kit/provideq/toolbox/process/GamsProcessRunner.java +++ b/src/main/java/edu/kit/provideq/toolbox/process/GamsProcessRunner.java @@ -26,7 +26,7 @@ public class GamsProcessRunner extends ProcessRunner { /** * The name of the file to call for running a GAMS script. */ - private static final String GAMS_EXECUTABLE_NAME = "gams"; + private static final String GAMS_EXECUTABLE_NAME = "/home/piotr/gams49.3_linux_x64_64_sfx/gams"; /** * Creates a process runner for a GAMS task. diff --git a/src/main/java/edu/kit/provideq/toolbox/process/PythonProcessRunner.java b/src/main/java/edu/kit/provideq/toolbox/process/PythonProcessRunner.java index b4b7c277..f01f699c 100644 --- a/src/main/java/edu/kit/provideq/toolbox/process/PythonProcessRunner.java +++ b/src/main/java/edu/kit/provideq/toolbox/process/PythonProcessRunner.java @@ -13,7 +13,8 @@ public class PythonProcessRunner extends ProcessRunner { /** * The name of the file to call for running a GAMS script. */ - private static final String PYTHON_EXECUTABLE_NAME = "python"; + private static final String PYTHON_EXECUTABLE_NAME + = "/home/piotr/anaconda3/envs/provideq-toolbox-server/bin/python"; /** * Creates a process runner for a Python script. diff --git a/src/main/resources/application.properties b/src/main/resources/application.properties index a717df3d..ef3c4daf 100644 --- a/src/main/resources/application.properties +++ b/src/main/resources/application.properties @@ -15,6 +15,9 @@ solvers.directory=solvers gams.directory=${solvers.directory}/gams gams.script.max-cut=${gams.directory}/max-cut/maxcut.gms gams.script.sat=${gams.directory}/sat/sat.gms +gams.script.mip=${gams.directory}/mip/mip.gms +gams.script.qubo_solve=${gams.directory}/mip/qubo_solve.gms + qiskit.directory=${solvers.directory}/qiskit qiskit.script.max-cut=${qiskit.directory}/max-cut/maxCut_qiskit.py diff --git a/src/main/resources/edu/kit/provideq/toolbox/mip/simple.lp b/src/main/resources/edu/kit/provideq/toolbox/mip/simple.lp index 89e74067..2582d3bf 100644 --- a/src/main/resources/edu/kit/provideq/toolbox/mip/simple.lp +++ b/src/main/resources/edu/kit/provideq/toolbox/mip/simple.lp @@ -1,17 +1,17 @@ Maximize - obj: 3 x + 4 y + obj: 3 x + 4 y Subject To - c1: x + 2 y <= 14 - c2: 3 x - y >= 0 - c3: x - y <= 2 + c1: x + 2 y <= 14 + c2: -3 x + y <= 0 + c3: x - y <= 2 Bounds - x >= 0 - y >= 0 + 0 <= x <= 6 + 0 <= y <= 7 Generals - x - y + x + y End From cb0dec4f745d1e96d249c4aa33b033a55810a236 Mon Sep 17 00:00:00 2001 From: Piotr Malkowski Date: Wed, 9 Apr 2025 13:20:57 +0200 Subject: [PATCH 14/37] fix: undo hardcoded gams and python paths --- .../edu/kit/provideq/toolbox/process/GamsProcessRunner.java | 2 +- .../edu/kit/provideq/toolbox/process/PythonProcessRunner.java | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/main/java/edu/kit/provideq/toolbox/process/GamsProcessRunner.java b/src/main/java/edu/kit/provideq/toolbox/process/GamsProcessRunner.java index 52f33b96..a00921b3 100644 --- a/src/main/java/edu/kit/provideq/toolbox/process/GamsProcessRunner.java +++ b/src/main/java/edu/kit/provideq/toolbox/process/GamsProcessRunner.java @@ -26,7 +26,7 @@ public class GamsProcessRunner extends ProcessRunner { /** * The name of the file to call for running a GAMS script. */ - private static final String GAMS_EXECUTABLE_NAME = "/home/piotr/gams49.3_linux_x64_64_sfx/gams"; + private static final String GAMS_EXECUTABLE_NAME = "gams"; /** * Creates a process runner for a GAMS task. diff --git a/src/main/java/edu/kit/provideq/toolbox/process/PythonProcessRunner.java b/src/main/java/edu/kit/provideq/toolbox/process/PythonProcessRunner.java index f01f699c..a06f7d39 100644 --- a/src/main/java/edu/kit/provideq/toolbox/process/PythonProcessRunner.java +++ b/src/main/java/edu/kit/provideq/toolbox/process/PythonProcessRunner.java @@ -14,7 +14,7 @@ public class PythonProcessRunner extends ProcessRunner { * The name of the file to call for running a GAMS script. */ private static final String PYTHON_EXECUTABLE_NAME - = "/home/piotr/anaconda3/envs/provideq-toolbox-server/bin/python"; + = "python"; /** * Creates a process runner for a Python script. From 10fc65c9f5676bd4705785dd673680e3c2fce9f0 Mon Sep 17 00:00:00 2001 From: Piotr Malkowski Date: Wed, 9 Apr 2025 13:27:05 +0200 Subject: [PATCH 15/37] chore: remove unused scripts --- scripts/mps-lp-to-gms.py | 32 ----------------- solvers/gams/mip/lp-to-gms.gms | 64 ---------------------------------- 2 files changed, 96 deletions(-) delete mode 100644 scripts/mps-lp-to-gms.py delete mode 100644 solvers/gams/mip/lp-to-gms.gms diff --git a/scripts/mps-lp-to-gms.py b/scripts/mps-lp-to-gms.py deleted file mode 100644 index c8fe8bea..00000000 --- a/scripts/mps-lp-to-gms.py +++ /dev/null @@ -1,32 +0,0 @@ -import os -import subprocess -import argparse - -parser = argparse.ArgumentParser( - prog="OR-Tools CBC MIP Solver", - description="A CLI Program to solve MIP problems in lp or mips format with Or-Tools CBC" -) - -parser.add_argument("input_file") -parser.add_argument("--output-file") -args = parser.parse_args() - - -input_file_abs = os.path.abspath(args.input_file) -output_file_abs = os.path.abspath(args.output_file) -# Directory where we'll run GAMS so it creates .gms in that same folder -work_dir = os.path.dirname(input_file_abs) -# GAMS by default names the new .gms file after the base of the input file -base_name = os.path.splitext(os.path.basename(input_file_abs))[0] -intermediate_gms = os.path.join(work_dir, base_name + ".gms") -# 1) Run GAMS to convert the LP/MPS -> GMS -# "task=ConvertIn" instructs GAMS to read the file and produce a .gms. -subprocess.run(["gams", input_file_abs, "task=ConvertIn"], check=True, cwd=work_dir) -# 2) Move or rename the newly created .gms to the user-specified output_file -# If you prefer to copy instead, you can use shutil.copyfile(). -if not os.path.exists(intermediate_gms): - raise FileNotFoundError(f"Expected GAMS output file not found: {intermediate_gms}") -os.replace(intermediate_gms, output_file_abs) - - - diff --git a/solvers/gams/mip/lp-to-gms.gms b/solvers/gams/mip/lp-to-gms.gms deleted file mode 100644 index de4f8e15..00000000 --- a/solvers/gams/mip/lp-to-gms.gms +++ /dev/null @@ -1,64 +0,0 @@ -$setEnv GAMSINP "%INPUT%" - -$onEmbeddedCode Python: -import os -import sys -import subprocess -import re - -inp_file = os.environ.get('GAMSINP') -print(f"{inp_file}") - -if not os.path.exists(inp_file): - sys.exit(f"*** ERROR: input file not found: {inp_file}") - -# Replace .lp with .gms -base, ext = os.path.splitext(inp_file) -new_name = base + '.gms' - -cmd = ['mps2gms', inp_file] -print(f"--- Running: {' '.join(cmd)}") -res = subprocess.run( - cmd, - stdout=subprocess.PIPE, - stderr=subprocess.PIPE -) -#out = res.stdout.decode('utf-8', errors='ignore') -#err = res.stderr.decode('utf-8', errors='ignore') - -with open(new_name, 'r', encoding='utf-8', errors='ignore') as f: - content = f.read() - -pattern = re.compile( - r'Solve\s+(\S+)\s+using\s+(\S+)\s+(maximizing|minimizing)\s+(\S+)\s*;', - re.IGNORECASE -) -match = pattern.search(content) -if not match: - sys.exit('*** ERROR: No "Solve ... using ..." statement found in generated .gms file') - -model_name = match.group(1) # e.g. "m" -model_type = match.group(2) # e.g. "LP" or "MIP" -sense_text = match.group(3) # "maximizing" or "minimizing" -obj_var = match.group(4) # e.g. "x7" - -if sense_text.lower().startswith('max'): - gams_sense = 'max' -else: - gams_sense = 'min' - -# 6) Show info to GAMS listing: -print(f'--- Detected from "{new_name}":') -print(f' Model name = {model_name}') -print(f' Model type = {model_type}') -print(f' Sense = {gams_sense}') -print(f' Objective var = {obj_var}') - -with open("parsed_macros.gms", "w") as macrofile: - macrofile.write(f"$set PARSED_MODELNAME {model_name}\n") - macrofile.write(f"$set PARSED_MODELTYPE {model_type}\n") - macrofile.write(f"$set PARSED_SENSE {gams_sense}\n") - macrofile.write(f"$set PARSED_OBJVAR {obj_var}\n") - macrofile.write(f"$set HARDCODED_PENALTY 10\n") - -$offEmbeddedCode From 0d1d9ab0923fd1c0de13327bdfe82aad121fcbbb Mon Sep 17 00:00:00 2001 From: Piotr Malkowski Date: Wed, 9 Apr 2025 13:27:25 +0200 Subject: [PATCH 16/37] feat: adapt MIP solvers logic --- .../java/edu/kit/provideq/toolbox/mip/MipConfiguration.java | 6 ++++-- .../edu/kit/provideq/toolbox/mip/solvers/OrToolsCbc.java | 3 +-- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java b/src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java index fcfa966e..b284332d 100644 --- a/src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java +++ b/src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java @@ -6,6 +6,7 @@ import edu.kit.provideq.toolbox.meta.Problem; import edu.kit.provideq.toolbox.meta.ProblemManager; import edu.kit.provideq.toolbox.meta.ProblemType; +import edu.kit.provideq.toolbox.mip.solvers.QuboMipSolver; import java.io.IOException; import java.util.Objects; import java.util.Set; @@ -32,11 +33,12 @@ public class MipConfiguration { @Bean ProblemManager getMipManager( OrToolsCbc orToolsCbc, + QuboMipSolver quboMipSolver, ResourceProvider resourceProvider ) { return new ProblemManager<>( MIP, - Set.of(orToolsCbc), + Set.of(orToolsCbc, quboMipSolver), loadExampleProblems(resourceProvider) ); } @@ -45,7 +47,7 @@ private Set> loadExampleProblems( ResourceProvider resourceProvider) { try { var problemInputStream = Objects.requireNonNull( - getClass().getResourceAsStream("simple.mps"), + getClass().getResourceAsStream("simple.lp"), "mixed-integer-problem example is unavailable!" ); var problem = new Problem<>(MIP); diff --git a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/OrToolsCbc.java b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/OrToolsCbc.java index 148dff6d..58664f35 100644 --- a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/OrToolsCbc.java +++ b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/OrToolsCbc.java @@ -42,14 +42,13 @@ public Mono> solve( ) { var solution = new Solution<>(this); - // TODO: enable both .lp and .mps Problems from the start var processResult = context .getBean(PythonProcessRunner.class, scriptPath) .withArguments( ProcessRunner.INPUT_FILE_PATH, ProcessRunner.OUTPUT_FILE_PATH ) - .writeInputFile(input, "problem.lp") + .writeInputFile(input) .readOutputFile() .run(getProblemType(), solution.getId()); From 654724febaad162a38f9045ba9c9f67cddd0e795 Mon Sep 17 00:00:00 2001 From: Piotr Malkowski Date: Wed, 9 Apr 2025 13:27:43 +0200 Subject: [PATCH 17/37] refactor: remove whitelines --- src/main/resources/edu/kit/provideq/toolbox/mip/simple.lp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/main/resources/edu/kit/provideq/toolbox/mip/simple.lp b/src/main/resources/edu/kit/provideq/toolbox/mip/simple.lp index 2582d3bf..533b6950 100644 --- a/src/main/resources/edu/kit/provideq/toolbox/mip/simple.lp +++ b/src/main/resources/edu/kit/provideq/toolbox/mip/simple.lp @@ -13,5 +13,4 @@ Bounds Generals x y - End From 98a7d2f8d911a9c68f237633f84cfe48e596db5e Mon Sep 17 00:00:00 2001 From: Piotr Malkowski Date: Wed, 9 Apr 2025 13:28:48 +0200 Subject: [PATCH 18/37] fix: adapt input arguments --- .../toolbox/mip/solvers/QuboMipSolver.java | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java index b300f40d..14fc2f33 100644 --- a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java +++ b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java @@ -17,7 +17,7 @@ public class QuboMipSolver extends MipSolver { private final String scriptPath; private final String dependencyScriptPath; private final ApplicationContext context; - public static final int PENALTY = 2; + public static final int PENALTY = 100; @Autowired public QuboMipSolver( @@ -48,14 +48,25 @@ public Mono> solve( var solution = new Solution<>(this); // Run GAMS via console - ProcessResult processResult = context + ProcessResult fileGeneratorResult = context .getBean(GamsProcessRunner.class, scriptPath) .withArguments( "--INPUT=" + ProcessRunner.INPUT_FILE_PATH, + "--PENALTY=" + String.valueOf(PENALTY) + ) + .writeInputFile(input) + .readOutputFile() + .run(getProblemType(), solution.getId()); + + String modifiedInputFilePath = ProcessRunner.INPUT_FILE_PATH.replaceFirst("\\.[^.]*$", ".gms"); + + ProcessResult processResult = context + .getBean(GamsProcessRunner.class, modifiedInputFilePath) + .withArguments( "--IDIR=" + dependencyScriptPath, "--SOLOUTPUT=" + ProcessRunner.OUTPUT_FILE_PATH ) - .writeInputFile(input, "problem.lp") + .writeInputFile(input) .readOutputFile() .run(getProblemType(), solution.getId()); return Mono.just(processResult.applyTo(solution)); From a9ab795c9abac6b0c3f90dc0cf69f056102dc5d3 Mon Sep 17 00:00:00 2001 From: Piotr Malkowski Date: Wed, 9 Apr 2025 13:29:18 +0200 Subject: [PATCH 19/37] feat: add testing for mip / lp solvers --- .../provideq/toolbox/api/MipSolversTest.java | 51 +++++++++++++++++++ 1 file changed, 51 insertions(+) create mode 100644 src/test/java/edu/kit/provideq/toolbox/api/MipSolversTest.java diff --git a/src/test/java/edu/kit/provideq/toolbox/api/MipSolversTest.java b/src/test/java/edu/kit/provideq/toolbox/api/MipSolversTest.java new file mode 100644 index 00000000..f09887b5 --- /dev/null +++ b/src/test/java/edu/kit/provideq/toolbox/api/MipSolversTest.java @@ -0,0 +1,51 @@ +package edu.kit.provideq.toolbox.api; + + +import static edu.kit.provideq.toolbox.mip.MipConfiguration.MIP; + +import edu.kit.provideq.toolbox.meta.ProblemManagerProvider; +import edu.kit.provideq.toolbox.meta.ProblemSolver; +import java.time.Duration; +import java.util.stream.Stream; +import org.junit.jupiter.api.BeforeEach; +import org.junit.jupiter.api.TestInstance; +import org.junit.jupiter.params.ParameterizedTest; +import org.junit.jupiter.params.provider.Arguments; +import org.junit.jupiter.params.provider.MethodSource; +import org.springframework.beans.factory.annotation.Autowired; +import org.springframework.boot.test.autoconfigure.web.servlet.AutoConfigureMockMvc; +import org.springframework.boot.test.context.SpringBootTest; +import org.springframework.test.web.reactive.server.WebTestClient; + +@TestInstance(TestInstance.Lifecycle.PER_CLASS) +@SpringBootTest +@AutoConfigureMockMvc +class MipSolversTest { + @Autowired + private WebTestClient client; + + @Autowired + private ProblemManagerProvider problemManagerProvider; + + @BeforeEach + void beforeEach() { + this.client = this.client.mutate() + .responseTimeout(Duration.ofSeconds(20)) + .build(); + } + + @SuppressWarnings("OptionalGetWithoutIsPresent") + Stream provideArguments() { + var problemManager = problemManagerProvider.findProblemManagerForType(MIP).get(); + + return ApiTestHelper.getAllArgumentCombinations(problemManager) + .map(list -> Arguments.of(list.get(0), list.get(1))); + } + + @ParameterizedTest + @MethodSource("provideArguments") + void testSatSolver(ProblemSolver solver, String input) { + var problem = ApiTestHelper.createProblem(client, solver, input, MIP); + ApiTestHelper.testSolution(problem); + } +} From 6f91ee3c7b1af0f351bd629bfb2131ea7657e563 Mon Sep 17 00:00:00 2001 From: Piotr Malkowski Date: Wed, 9 Apr 2025 13:29:39 +0200 Subject: [PATCH 20/37] fix: adapt output parsing logic in gams mip solver --- solvers/gams/mip/mip.gms | 24 ++++++++++++++++++++++-- 1 file changed, 22 insertions(+), 2 deletions(-) diff --git a/solvers/gams/mip/mip.gms b/solvers/gams/mip/mip.gms index 7dfa68a6..3c06a48a 100644 --- a/solvers/gams/mip/mip.gms +++ b/solvers/gams/mip/mip.gms @@ -48,13 +48,33 @@ else: # Create the replacement line using the matched groups penalty = os.environ.get("GAMSPEN", "100") # Default to 100 if not set -new_line = f'$batinclude %IDIR% {model_name} {model_type} {gams_sense} {obj_var} {penalty} -logOn=2;' +new_line = f'$batinclude %IDIR% {model_name} {model_type} {gams_sense} {obj_var} {penalty} -logOn=2' # Replace the matching "Solve ..." line with the new line (only the first occurrence) new_content = pattern.sub(new_line, content, count=1) +solution_block = r''' +$setNames "%INPUT%" fp fn fe +$if not set SOLOUTPUT $set SOLOUTPUT %fp%solution%fe% +file solFile / "%SOLOUTPUT%" /; +put solFile 'c Solution of %INPUT%'; +put solFile / 's m '; + +if((m.modelstat = %modelStat.Optimal% or m.modelstat = %modelStat.IntegerSolution%) or m.solvestat = %solveStat.NormalCompletion%, + put solFile / 'Optimal solution found.'; + put solFile / 'Objective value: ' obj.l:0:0; + put solFile / 'Variable values:'; + loop(j, + put solFile / j.tl:0 ' = ' xi.l(j):0:0; + ); +else + put solFile / 'Error: no solution was returned.'; +); +''' + +# Append the solution parsing code at the end of the .gms file. +new_content += "\n" + solution_block # Write the modified content back to the .gms file with open(new_name, 'w', encoding='utf-8', errors='ignore') as f: f.write(new_content) - $offEmbeddedCode From 2d14e3464b54f1b41d0161f4043a7ffa4f6c4f33 Mon Sep 17 00:00:00 2001 From: "piotr.malkowski" Date: Wed, 9 Apr 2025 17:31:03 +0200 Subject: [PATCH 21/37] fix: repair input paths for qubo mip solver --- .../toolbox/mip/solvers/QuboMipSolver.java | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java index 14fc2f33..a094998e 100644 --- a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java +++ b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java @@ -6,6 +6,7 @@ import edu.kit.provideq.toolbox.process.GamsProcessRunner; import edu.kit.provideq.toolbox.process.ProcessResult; import edu.kit.provideq.toolbox.process.ProcessRunner; +import java.io.File; import org.springframework.beans.factory.annotation.Autowired; import org.springframework.beans.factory.annotation.Value; import org.springframework.context.ApplicationContext; @@ -54,20 +55,19 @@ public Mono> solve( "--INPUT=" + ProcessRunner.INPUT_FILE_PATH, "--PENALTY=" + String.valueOf(PENALTY) ) - .writeInputFile(input) - .readOutputFile() + .writeInputFile(input, "problem.lp") + .readOutputFile("problem.lp") .run(getProblemType(), solution.getId()); - String modifiedInputFilePath = ProcessRunner.INPUT_FILE_PATH.replaceFirst("\\.[^.]*$", ".gms"); + String relativePath = "jobs/mip/" + solution.getId() + "/problem.gms"; ProcessResult processResult = context - .getBean(GamsProcessRunner.class, modifiedInputFilePath) + .getBean(GamsProcessRunner.class, relativePath) .withArguments( - "--IDIR=" + dependencyScriptPath, + "--IDIR=" + new File(dependencyScriptPath).getAbsolutePath(), "--SOLOUTPUT=" + ProcessRunner.OUTPUT_FILE_PATH ) - .writeInputFile(input) - .readOutputFile() + .readOutputFile("output.txt") .run(getProblemType(), solution.getId()); return Mono.just(processResult.applyTo(solution)); } From 33ae13715555ea71bfa1671616739f4f36761092 Mon Sep 17 00:00:00 2001 From: "piotr.malkowski" Date: Wed, 9 Apr 2025 18:01:38 +0200 Subject: [PATCH 22/37] chore: remove unused files --- solvers/ortools/mips/lp_mps_converter.py | 10 ---------- .../edu/kit/provideq/toolbox/mip/simple.mps | 18 ------------------ 2 files changed, 28 deletions(-) delete mode 100644 solvers/ortools/mips/lp_mps_converter.py delete mode 100644 src/main/resources/edu/kit/provideq/toolbox/mip/simple.mps diff --git a/solvers/ortools/mips/lp_mps_converter.py b/solvers/ortools/mips/lp_mps_converter.py deleted file mode 100644 index 6dbea81f..00000000 --- a/solvers/ortools/mips/lp_mps_converter.py +++ /dev/null @@ -1,10 +0,0 @@ -import cplex - -def convert_model(input_file, output_file): - """ - Loads a model (LP or MPS) into CPLEX and writes it back out - in the requested format (.lp or .mps). - """ - # Create a CPLEX object and read in the model - model = cplex.Cplex(input_file) - model.write(output_file) diff --git a/src/main/resources/edu/kit/provideq/toolbox/mip/simple.mps b/src/main/resources/edu/kit/provideq/toolbox/mip/simple.mps deleted file mode 100644 index 24b8d833..00000000 --- a/src/main/resources/edu/kit/provideq/toolbox/mip/simple.mps +++ /dev/null @@ -1,18 +0,0 @@ -NAME EXAMPLE -ROWS - N COST <- (N) indicates the objective row - G C1 <- (G) 'Greater than or equal' constraint - L C2 <- (L) 'Less than or equal' constraint -COLUMNS - X COST 1 <- Objective coefficient for X is 1 - X C1 1 <- X's coefficient in constraint C1 - Y COST 3 <- Objective coefficient for Y is 3 - Y C1 1 <- Y's coefficient in constraint C1 - Y C2 1 <- Y's coefficient in constraint C2 -RHS - RHS C1 4 <- Right-hand side of constraint C1 - RHS C2 3 <- Right-hand side of constraint C2 -BOUNDS - LO BND X 0 <- Lower bound of X is 0 - LO BND Y 0 <- Lower bound of Y is 0 -ENDATA From 539e8d0862cda0ac300b09ec31c8e04cda486b85 Mon Sep 17 00:00:00 2001 From: "piotr.malkowski" Date: Wed, 9 Apr 2025 18:01:51 +0200 Subject: [PATCH 23/37] chore: remove unused files --- solvers/ortools/mips/cbc.py | 33 --------------------------------- 1 file changed, 33 deletions(-) delete mode 100644 solvers/ortools/mips/cbc.py diff --git a/solvers/ortools/mips/cbc.py b/solvers/ortools/mips/cbc.py deleted file mode 100644 index f7e9bdbe..00000000 --- a/solvers/ortools/mips/cbc.py +++ /dev/null @@ -1,33 +0,0 @@ -import argparse -import os -from ortools.linear_solver import pywraplp -from lp_mps_converter import convert_model - - -parser = argparse.ArgumentParser( - prog="OR-Tools CBC MIP Solver", - description="A CLI Program to solve MIP problems in lp or mips format with Or-Tools CBC" -) - -parser.add_argument("input_file") -parser.add_argument("--output-file") -args = parser.parse_args() - -input_ext = os.path.splitext(args.input_file)[1].lower() -if input_ext == ".lp": - convert_model(args.input_file, args.output_file) - mps_file = args.output_file -else: - mps_file = args.input_file - - -solver = pywraplp.Solver.CreateSolver("CBC") -model = solver.loadModelFromFile(mps_file) -status = solver.Solve() - -if status == pywraplp.Solver.OPTIMAL: - print("Solution:") - print("Objective value =", solver.Objective().Value()) -else: - print("The problem does not have an optimal solution.") - From 5ebd946353d0d81b1fd14a3862697e3ed848c822 Mon Sep 17 00:00:00 2001 From: "piotr.malkowski" Date: Wed, 9 Apr 2025 18:02:20 +0200 Subject: [PATCH 24/37] chore: adapt dependencies --- solvers/cplex/requirements.txt | 1 + 1 file changed, 1 insertion(+) create mode 100644 solvers/cplex/requirements.txt diff --git a/solvers/cplex/requirements.txt b/solvers/cplex/requirements.txt new file mode 100644 index 00000000..4cfbd949 --- /dev/null +++ b/solvers/cplex/requirements.txt @@ -0,0 +1 @@ +cplex From 0a435fd3c2060b04188d480a5179ae8031b9a547 Mon Sep 17 00:00:00 2001 From: "piotr.malkowski" Date: Wed, 9 Apr 2025 18:02:36 +0200 Subject: [PATCH 25/37] feat: add cplex mip solver --- solvers/cplex/mips/mip.py | 40 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) create mode 100644 solvers/cplex/mips/mip.py diff --git a/solvers/cplex/mips/mip.py b/solvers/cplex/mips/mip.py new file mode 100644 index 00000000..363771b4 --- /dev/null +++ b/solvers/cplex/mips/mip.py @@ -0,0 +1,40 @@ +import argparse +import cplex + +parser = argparse.ArgumentParser( + description="Solve a MIP problem from an MPS file using IBM CPLEX." +) +parser.add_argument("input_file", help="Path to the input MPS file.") +parser.add_argument("output_file", help="Path to write the solution.") +args = parser.parse_args() + +try: + cpx = cplex.Cplex() + cpx.read(args.input_file) +except cplex.exceptions.CplexError as exc: + print("Error reading the MPS file:", exc) + exit() + +try: + # Solve the model. + cpx.solve() +except cplex.exceptions.CplexError as exc: + print("Error during solve:", exc) + exit() + +status = cpx.solution.get_status() +# See the full list of status codes in the CPLEX documentation. +with open(args.output_file, 'w') as out_file: + if status == cpx.solution.status.optimal: + out_file.write("Solution is optimal.\n") + out_file.write("Objective value = {}\n".format(cpx.solution.get_objective_value())) + # Get the variable names and values. + var_names = cpx.variables.get_names() + var_values = cpx.solution.get_values() + for name, val in zip(var_names, var_values): + out_file.write("{:<15} = {}\n".format(name, val)) + else: + out_file.write("No optimal solution found. Status: {}\n".format(status)) + # Optionally, you can also write additional status information: + out_file.write("Solution status string: {}\n".format(cpx.solution.get_status_string())) +print("Solution written to", args.output_file) From a3bce39e5aaa44a6e85c6834a276c0eaf7a7556b Mon Sep 17 00:00:00 2001 From: "piotr.malkowski" Date: Wed, 9 Apr 2025 18:02:50 +0200 Subject: [PATCH 26/37] chore: adapt properties for cplex solver --- src/main/resources/application.properties | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/main/resources/application.properties b/src/main/resources/application.properties index ef3c4daf..af02e2df 100644 --- a/src/main/resources/application.properties +++ b/src/main/resources/application.properties @@ -24,8 +24,8 @@ qiskit.script.max-cut=${qiskit.directory}/max-cut/maxCut_qiskit.py qiskit.script.qubo=${qiskit.directory}/qubo/qubo_qiskit.py qiskit.script.knapsack=${qiskit.directory}/knapsack/knapsack_qiskit.py -ortools.directory=${solvers.directory}/ortools -ortools.script.ormip=${ortools.directory}/mips/cbc.py +cplex.solver.directory=${solvers.directory}/cplex +cplex.script.mip=${cplex.solver.directory}/mips/mip.py cirq.directory=${solvers.directory}/cirq cirq.script.max-cut=${cirq.directory}/max-cut/max_cut_cirq.py From 10e54a68f3326c391fc0d2ce961082ad885c5e50 Mon Sep 17 00:00:00 2001 From: "piotr.malkowski" Date: Wed, 9 Apr 2025 18:03:09 +0200 Subject: [PATCH 27/37] refactor: rename directory --- solvers/ortools/requirements.txt | 2 -- 1 file changed, 2 deletions(-) delete mode 100644 solvers/ortools/requirements.txt diff --git a/solvers/ortools/requirements.txt b/solvers/ortools/requirements.txt deleted file mode 100644 index 5367027c..00000000 --- a/solvers/ortools/requirements.txt +++ /dev/null @@ -1,2 +0,0 @@ -ortools -cplex From a736af46937e08822bc342388666cdb4dea5f18e Mon Sep 17 00:00:00 2001 From: "piotr.malkowski" Date: Wed, 9 Apr 2025 18:04:08 +0200 Subject: [PATCH 28/37] refactor: rename ortools occurances with cplex --- .../edu/kit/provideq/toolbox/mip/MipConfiguration.java | 6 +++--- .../mip/solvers/{OrToolsCbc.java => CplexMip.java} | 8 ++++---- 2 files changed, 7 insertions(+), 7 deletions(-) rename src/main/java/edu/kit/provideq/toolbox/mip/solvers/{OrToolsCbc.java => CplexMip.java} (91%) diff --git a/src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java b/src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java index b284332d..2e37f442 100644 --- a/src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java +++ b/src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java @@ -2,7 +2,7 @@ import edu.kit.provideq.toolbox.ResourceProvider; import edu.kit.provideq.toolbox.exception.MissingExampleException; -import edu.kit.provideq.toolbox.mip.solvers.OrToolsCbc; +import edu.kit.provideq.toolbox.mip.solvers.CplexMip; import edu.kit.provideq.toolbox.meta.Problem; import edu.kit.provideq.toolbox.meta.ProblemManager; import edu.kit.provideq.toolbox.meta.ProblemType; @@ -32,13 +32,13 @@ public class MipConfiguration { @Bean ProblemManager getMipManager( - OrToolsCbc orToolsCbc, + CplexMip cplexMip, QuboMipSolver quboMipSolver, ResourceProvider resourceProvider ) { return new ProblemManager<>( MIP, - Set.of(orToolsCbc, quboMipSolver), + Set.of(cplexMip, quboMipSolver), loadExampleProblems(resourceProvider) ); } diff --git a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/OrToolsCbc.java b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/CplexMip.java similarity index 91% rename from src/main/java/edu/kit/provideq/toolbox/mip/solvers/OrToolsCbc.java rename to src/main/java/edu/kit/provideq/toolbox/mip/solvers/CplexMip.java index 58664f35..01496df2 100644 --- a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/OrToolsCbc.java +++ b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/CplexMip.java @@ -12,13 +12,13 @@ import reactor.core.publisher.Mono; @Component -public class OrToolsCbc extends MipSolver { +public class CplexMip extends MipSolver { private final String scriptPath; private final ApplicationContext context; @Autowired - public OrToolsCbc( - @Value("${ortools.script.ormip}") String scriptPath, + public CplexMip( + @Value("${cplex.script.mip}") String scriptPath, ApplicationContext context) { this.scriptPath = scriptPath; this.context = context; @@ -49,7 +49,7 @@ public Mono> solve( ProcessRunner.OUTPUT_FILE_PATH ) .writeInputFile(input) - .readOutputFile() + .readOutputFile("output.txt") .run(getProblemType(), solution.getId()); // Return if process failed From 1957b3c214e7c9c4bde9282694f78bc412f5c55d Mon Sep 17 00:00:00 2001 From: "piotr.malkowski" Date: Wed, 9 Apr 2025 19:40:07 +0200 Subject: [PATCH 29/37] refactor: add comments, expand explanations --- solvers/cplex/mips/mip.py | 4 ---- .../provideq/toolbox/mip/MipConfiguration.java | 2 +- .../provideq/toolbox/mip/solvers/CplexMip.java | 9 +++++++-- .../provideq/toolbox/mip/solvers/MipSolver.java | 5 ++++- .../toolbox/mip/solvers/QuboMipSolver.java | 15 ++++++++++----- 5 files changed, 22 insertions(+), 13 deletions(-) diff --git a/solvers/cplex/mips/mip.py b/solvers/cplex/mips/mip.py index 363771b4..3eb90393 100644 --- a/solvers/cplex/mips/mip.py +++ b/solvers/cplex/mips/mip.py @@ -16,25 +16,21 @@ exit() try: - # Solve the model. cpx.solve() except cplex.exceptions.CplexError as exc: print("Error during solve:", exc) exit() status = cpx.solution.get_status() -# See the full list of status codes in the CPLEX documentation. with open(args.output_file, 'w') as out_file: if status == cpx.solution.status.optimal: out_file.write("Solution is optimal.\n") out_file.write("Objective value = {}\n".format(cpx.solution.get_objective_value())) - # Get the variable names and values. var_names = cpx.variables.get_names() var_values = cpx.solution.get_values() for name, val in zip(var_names, var_values): out_file.write("{:<15} = {}\n".format(name, val)) else: out_file.write("No optimal solution found. Status: {}\n".format(status)) - # Optionally, you can also write additional status information: out_file.write("Solution status string: {}\n".format(cpx.solution.get_status_string())) print("Solution written to", args.output_file) diff --git a/src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java b/src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java index 2e37f442..3babe15c 100644 --- a/src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java +++ b/src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java @@ -22,7 +22,7 @@ public class MipConfiguration { * MIP (Mixed integer Problem) * A mixed integer optimization problem. * For a given integer term, - * find the minimal variable assignment of the term. + * find the minimal / maximal variable assignment of the term. */ public static final ProblemType MIP = new ProblemType<>( "mip", diff --git a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/CplexMip.java b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/CplexMip.java index 01496df2..b424e4d0 100644 --- a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/CplexMip.java +++ b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/CplexMip.java @@ -3,14 +3,19 @@ import edu.kit.provideq.toolbox.Solution; import edu.kit.provideq.toolbox.meta.SolvingProperties; import edu.kit.provideq.toolbox.meta.SubRoutineResolver; +import edu.kit.provideq.toolbox.mip.MipConfiguration; import edu.kit.provideq.toolbox.process.ProcessRunner; import edu.kit.provideq.toolbox.process.PythonProcessRunner; +import edu.kit.provideq.toolbox.sat.SatConfiguration; import org.springframework.beans.factory.annotation.Autowired; import org.springframework.beans.factory.annotation.Value; import org.springframework.context.ApplicationContext; import org.springframework.stereotype.Component; import reactor.core.publisher.Mono; +/** + * {@link MipConfiguration#MIP} solver using IBM cplex. + */ @Component public class CplexMip extends MipSolver { private final String scriptPath; @@ -26,12 +31,12 @@ public CplexMip( @Override public String getName() { - return "(OR-Tools) Cbc Solver for MIP"; + return "IBM Cplex Solver for MIP"; } @Override public String getDescription() { - return "This solver uses OR-Tools CBC Solver to solve MIP problems"; + return "This solver uses IBM Cplex Solver to solve MIP problems"; } @Override diff --git a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/MipSolver.java b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/MipSolver.java index cd345e70..aedd52f7 100644 --- a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/MipSolver.java +++ b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/MipSolver.java @@ -1,9 +1,12 @@ package edu.kit.provideq.toolbox.mip.solvers; -import edu.kit.provideq.toolbox.mip.MipConfiguration; import edu.kit.provideq.toolbox.meta.ProblemSolver; import edu.kit.provideq.toolbox.meta.ProblemType; +import edu.kit.provideq.toolbox.mip.MipConfiguration; +/** + * A solver for Mixed Integer Optimization problems. + */ public abstract class MipSolver implements ProblemSolver { @Override public ProblemType getProblemType() { diff --git a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java index a094998e..e97e1f3b 100644 --- a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java +++ b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java @@ -3,6 +3,7 @@ import edu.kit.provideq.toolbox.Solution; import edu.kit.provideq.toolbox.meta.SolvingProperties; import edu.kit.provideq.toolbox.meta.SubRoutineResolver; +import edu.kit.provideq.toolbox.mip.MipConfiguration; import edu.kit.provideq.toolbox.process.GamsProcessRunner; import edu.kit.provideq.toolbox.process.ProcessResult; import edu.kit.provideq.toolbox.process.ProcessRunner; @@ -13,12 +14,15 @@ import org.springframework.stereotype.Component; import reactor.core.publisher.Mono; +/** + * {@link MipConfiguration#MIP} solver using QUBO Reformulation implementing gams. + */ @Component public class QuboMipSolver extends MipSolver { private final String scriptPath; private final String dependencyScriptPath; private final ApplicationContext context; - public static final int PENALTY = 100; + private static final int PENALTY = 100; @Autowired public QuboMipSolver( @@ -32,12 +36,13 @@ public QuboMipSolver( @Override public String getName() { - return "TODO"; + return "Gams Solver for MIP"; } @Override public String getDescription() { - return "TODO"; + return "This solver uses LP / MIP reformulation to QUBO which " + + "is then subsequently solved by gams"; } @Override @@ -49,11 +54,11 @@ public Mono> solve( var solution = new Solution<>(this); // Run GAMS via console - ProcessResult fileGeneratorResult = context + context .getBean(GamsProcessRunner.class, scriptPath) .withArguments( "--INPUT=" + ProcessRunner.INPUT_FILE_PATH, - "--PENALTY=" + String.valueOf(PENALTY) + "--PENALTY=" + PENALTY ) .writeInputFile(input, "problem.lp") .readOutputFile("problem.lp") From 71f571113beecd770fdf3f8c4fb11bd642849545 Mon Sep 17 00:00:00 2001 From: "piotr.malkowski" Date: Wed, 9 Apr 2025 19:43:20 +0200 Subject: [PATCH 30/37] refactor: import checkstyle --- .../java/edu/kit/provideq/toolbox/mip/MipConfiguration.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java b/src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java index 3babe15c..60dc88dd 100644 --- a/src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java +++ b/src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java @@ -2,10 +2,10 @@ import edu.kit.provideq.toolbox.ResourceProvider; import edu.kit.provideq.toolbox.exception.MissingExampleException; -import edu.kit.provideq.toolbox.mip.solvers.CplexMip; import edu.kit.provideq.toolbox.meta.Problem; import edu.kit.provideq.toolbox.meta.ProblemManager; import edu.kit.provideq.toolbox.meta.ProblemType; +import edu.kit.provideq.toolbox.mip.solvers.CplexMip; import edu.kit.provideq.toolbox.mip.solvers.QuboMipSolver; import java.io.IOException; import java.util.Objects; From 4c53204b2dcb0ed8218ffe1c6b4c2392f67ae5b5 Mon Sep 17 00:00:00 2001 From: "piotr.malkowski" Date: Wed, 9 Apr 2025 21:39:18 +0200 Subject: [PATCH 31/37] fix: adapt cplex logic to always write solution, not only in optimal case --- solvers/cplex/mips/mip.py | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/solvers/cplex/mips/mip.py b/solvers/cplex/mips/mip.py index 3eb90393..4c293448 100644 --- a/solvers/cplex/mips/mip.py +++ b/solvers/cplex/mips/mip.py @@ -15,6 +15,7 @@ print("Error reading the MPS file:", exc) exit() + try: cpx.solve() except cplex.exceptions.CplexError as exc: @@ -23,14 +24,10 @@ status = cpx.solution.get_status() with open(args.output_file, 'w') as out_file: - if status == cpx.solution.status.optimal: - out_file.write("Solution is optimal.\n") - out_file.write("Objective value = {}\n".format(cpx.solution.get_objective_value())) - var_names = cpx.variables.get_names() - var_values = cpx.solution.get_values() - for name, val in zip(var_names, var_values): - out_file.write("{:<15} = {}\n".format(name, val)) - else: - out_file.write("No optimal solution found. Status: {}\n".format(status)) - out_file.write("Solution status string: {}\n".format(cpx.solution.get_status_string())) + out_file.write("Solution status: {}\n".format(cpx.solution.get_status_string())) + out_file.write("Objective value = {}\n".format(cpx.solution.get_objective_value())) + var_names = cpx.variables.get_names() + var_values = cpx.solution.get_values() + for name, val in zip(var_names, var_values): + out_file.write("{:<15} = {}\n".format(name, val)) print("Solution written to", args.output_file) From d805bdd3db9d685bec06ce3051106739947292c1 Mon Sep 17 00:00:00 2001 From: "piotr.malkowski" Date: Thu, 10 Apr 2025 12:05:10 +0200 Subject: [PATCH 32/37] refactor: line brakes --- .../edu/kit/provideq/toolbox/process/PythonProcessRunner.java | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/main/java/edu/kit/provideq/toolbox/process/PythonProcessRunner.java b/src/main/java/edu/kit/provideq/toolbox/process/PythonProcessRunner.java index a06f7d39..b4b7c277 100644 --- a/src/main/java/edu/kit/provideq/toolbox/process/PythonProcessRunner.java +++ b/src/main/java/edu/kit/provideq/toolbox/process/PythonProcessRunner.java @@ -13,8 +13,7 @@ public class PythonProcessRunner extends ProcessRunner { /** * The name of the file to call for running a GAMS script. */ - private static final String PYTHON_EXECUTABLE_NAME - = "python"; + private static final String PYTHON_EXECUTABLE_NAME = "python"; /** * Creates a process runner for a Python script. From 79e7652054c915465901b6db118971e663bc3bb0 Mon Sep 17 00:00:00 2001 From: "piotr.malkowski" Date: Thu, 10 Apr 2025 12:44:05 +0200 Subject: [PATCH 33/37] refactor: add comment in gams mip solver --- solvers/gams/mip/mip.gms | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/solvers/gams/mip/mip.gms b/solvers/gams/mip/mip.gms index 3c06a48a..0153bb29 100644 --- a/solvers/gams/mip/mip.gms +++ b/solvers/gams/mip/mip.gms @@ -13,7 +13,7 @@ print(f"{inp_file}") if not os.path.exists(inp_file): sys.exit(f"*** ERROR: input file not found: {inp_file}") -# Replace .lp with .gms +# Replace .lp or .mps with .gms base, ext = os.path.splitext(inp_file) new_name = base + '.gms' From efb770f8f1118b405f8e5d873f20ba3c481bb258 Mon Sep 17 00:00:00 2001 From: Piotr Malkowski Date: Tue, 22 Apr 2025 13:05:59 +0200 Subject: [PATCH 34/37] fix: undo merge PythonProcessRunner change --- .../toolbox/process/PythonProcessRunner.java | 52 ++++++++++++++++++- 1 file changed, 50 insertions(+), 2 deletions(-) diff --git a/src/main/java/edu/kit/provideq/toolbox/process/PythonProcessRunner.java b/src/main/java/edu/kit/provideq/toolbox/process/PythonProcessRunner.java index b4b7c277..212317f1 100644 --- a/src/main/java/edu/kit/provideq/toolbox/process/PythonProcessRunner.java +++ b/src/main/java/edu/kit/provideq/toolbox/process/PythonProcessRunner.java @@ -1,5 +1,7 @@ package edu.kit.provideq.toolbox.process; +import java.util.List; +import java.util.Optional; import org.springframework.beans.factory.config.ConfigurableBeanFactory; import org.springframework.context.annotation.Scope; import org.springframework.stereotype.Component; @@ -19,10 +21,56 @@ public class PythonProcessRunner extends ProcessRunner { * Creates a process runner for a Python script. * * @param scriptPath the filepath of the Python script to run. + * @param venvName the name of the virtual environment to use. */ - public PythonProcessRunner(String scriptPath) { + public PythonProcessRunner(String scriptPath, String venvName) { super(new ProcessBuilder()); - withArguments(PYTHON_EXECUTABLE_NAME, scriptPath); + String osName = System.getProperty("os.name").toLowerCase(); + if (osName.contains("win")) { + withArguments( + String.format("./venv/%s/Scripts/activate.bat", venvName), + "&&", + PYTHON_EXECUTABLE_NAME, + scriptPath); + } else { + processBuilder.command( + "sh", + "-c", + String.format(". ./venv/%s/bin/activate && python %s", venvName, scriptPath)); + } + } + + @Override + public ProcessRunner withArguments(String... arguments) { + String osName = System.getProperty("os.name").toLowerCase(); + if (osName.contains("win")) { + return super.withArguments(arguments); + } + + preProcessors.add((problemType, solutionId) -> { + for (String argument : arguments) { + for (var transformer : argumentTransformers) { + argument = transformer.apply(argument); + } + + if (processBuilder.command().isEmpty()) { + processBuilder.command(argument); + } else { + List existingCommands = processBuilder.command(); + int lastIndex = existingCommands.size() - 1; + String lastCommand = existingCommands.get(lastIndex); + existingCommands.remove(lastIndex); + + var newArgument = lastCommand + " " + argument; + existingCommands.add(newArgument); + processBuilder.command(existingCommands); + } + } + + return Optional.empty(); + }); + + return this; } } From 9f4b02076a14bcd7a409bbda0e26406e535e1c1d Mon Sep 17 00:00:00 2001 From: Piotr Malkowski Date: Tue, 22 Apr 2025 13:06:58 +0200 Subject: [PATCH 35/37] feat: add mip cplex solver --- solvers/cplex/mip/mip.py | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) create mode 100644 solvers/cplex/mip/mip.py diff --git a/solvers/cplex/mip/mip.py b/solvers/cplex/mip/mip.py new file mode 100644 index 00000000..4c293448 --- /dev/null +++ b/solvers/cplex/mip/mip.py @@ -0,0 +1,33 @@ +import argparse +import cplex + +parser = argparse.ArgumentParser( + description="Solve a MIP problem from an MPS file using IBM CPLEX." +) +parser.add_argument("input_file", help="Path to the input MPS file.") +parser.add_argument("output_file", help="Path to write the solution.") +args = parser.parse_args() + +try: + cpx = cplex.Cplex() + cpx.read(args.input_file) +except cplex.exceptions.CplexError as exc: + print("Error reading the MPS file:", exc) + exit() + + +try: + cpx.solve() +except cplex.exceptions.CplexError as exc: + print("Error during solve:", exc) + exit() + +status = cpx.solution.get_status() +with open(args.output_file, 'w') as out_file: + out_file.write("Solution status: {}\n".format(cpx.solution.get_status_string())) + out_file.write("Objective value = {}\n".format(cpx.solution.get_objective_value())) + var_names = cpx.variables.get_names() + var_values = cpx.solution.get_values() + for name, val in zip(var_names, var_values): + out_file.write("{:<15} = {}\n".format(name, val)) +print("Solution written to", args.output_file) From e126e753d89184bdad60ae3a7e91ea5021b36af0 Mon Sep 17 00:00:00 2001 From: Piotr Malkowski Date: Tue, 6 May 2025 10:51:08 +0200 Subject: [PATCH 36/37] refactor: unused imports, naming schemes --- solvers/cplex/mip/mip.py | 5 +++-- .../gams/mip/{mip.gms => mip_to_qubo_reformulation.gms} | 0 .../edu/kit/provideq/toolbox/mip/solvers/CplexMip.java | 1 - .../kit/provideq/toolbox/mip/solvers/QuboMipSolver.java | 8 ++++---- src/main/resources/application.properties | 2 +- .../java/edu/kit/provideq/toolbox/api/MipSolversTest.java | 2 +- 6 files changed, 9 insertions(+), 9 deletions(-) rename solvers/gams/mip/{mip.gms => mip_to_qubo_reformulation.gms} (100%) diff --git a/solvers/cplex/mip/mip.py b/solvers/cplex/mip/mip.py index 4c293448..acf30d15 100644 --- a/solvers/cplex/mip/mip.py +++ b/solvers/cplex/mip/mip.py @@ -1,5 +1,6 @@ import argparse import cplex +import sys parser = argparse.ArgumentParser( description="Solve a MIP problem from an MPS file using IBM CPLEX." @@ -13,14 +14,14 @@ cpx.read(args.input_file) except cplex.exceptions.CplexError as exc: print("Error reading the MPS file:", exc) - exit() + sys.exit(1) try: cpx.solve() except cplex.exceptions.CplexError as exc: print("Error during solve:", exc) - exit() + sys.exit(1) status = cpx.solution.get_status() with open(args.output_file, 'w') as out_file: diff --git a/solvers/gams/mip/mip.gms b/solvers/gams/mip/mip_to_qubo_reformulation.gms similarity index 100% rename from solvers/gams/mip/mip.gms rename to solvers/gams/mip/mip_to_qubo_reformulation.gms diff --git a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/CplexMip.java b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/CplexMip.java index 8101819c..7569dae6 100644 --- a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/CplexMip.java +++ b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/CplexMip.java @@ -6,7 +6,6 @@ import edu.kit.provideq.toolbox.mip.MipConfiguration; import edu.kit.provideq.toolbox.process.ProcessRunner; import edu.kit.provideq.toolbox.process.PythonProcessRunner; -import edu.kit.provideq.toolbox.sat.SatConfiguration; import org.springframework.beans.factory.annotation.Autowired; import org.springframework.beans.factory.annotation.Value; import org.springframework.context.ApplicationContext; diff --git a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java index dc8bd849..a4f9df76 100644 --- a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java +++ b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java @@ -26,11 +26,11 @@ public class QuboMipSolver extends MipSolver { @Autowired public QuboMipSolver( - @Value("${path.gams.mip}") String scriptPath, - @Value("${path.gams.qubo-solve}") String dependencyScriptPath, + @Value("${path.gams.mip-to-qubo-reformulation}") String mipToQuboReformulatorPath, + @Value("${path.gams.qubo-solve}") String solverScriptPath, ApplicationContext context) { - this.scriptPath = scriptPath; - this.dependencyScriptPath = dependencyScriptPath; + this.scriptPath = mipToQuboReformulatorPath; + this.dependencyScriptPath = solverScriptPath; this.context = context; } diff --git a/src/main/resources/application.properties b/src/main/resources/application.properties index 05d4d027..63bdcae2 100644 --- a/src/main/resources/application.properties +++ b/src/main/resources/application.properties @@ -1 +1 @@ -# default spring profile, correct one will be set during runtime (see ToolboxServerApplication.java) # options: mac, windows, linux spring.profiles.active=linux springdoc.swagger-ui.operationsSorter=alpha springdoc.swagger-ui.tagsSorter=alpha working.directory=jobs examples.directory=examples springdoc.swagger-ui.path=/ # Solvers name.solvers=solvers # Non OS-specific solvers: (typically GAMS and Python) name.gams=gams path.gams=${name.solvers}/${name.gams} name.gams.max-cut=max-cut path.gams.max-cut=${path.gams}/${name.gams.max-cut}/maxcut.gms name.gams.sat=sat path.gams.sat=${path.gams}/${name.gams.sat}/sat.gms name.gams.mip=mip path.gams.mip=${path.gams}/${name.gams.mip}/mip.gms path.gams.qubo-solve=${path.gams}/${name.gams.mip}/qubo_solve.gms name.qiskit=qiskit path.qiskit=${name.solvers}/${name.qiskit} name.qiskit.knapsack=knapsack path.qiskit.knapsack=${path.qiskit}/${name.qiskit.knapsack}/knapsack_qiskit.py venv.qiskit.knapsack=${name.solvers}_${name.qiskit}_${name.qiskit.knapsack} name.qiskit.max-cut=max-cut path.qiskit.max-cut=${path.qiskit}/${name.qiskit.max-cut}/maxCut_qiskit.py venv.qiskit.max-cut=${name.solvers}_${name.qiskit}_${name.qiskit.max-cut} name.qiskit.qubo=qubo path.qiskit.qubo=${path.qiskit}/${name.qiskit.qubo}/qubo_qiskit.py venv.qiskit.qubo=${name.solvers}_${name.qiskit}_${name.qiskit.qubo} name.cirq=cirq path.cirq=${name.solvers}/${name.cirq} name.cirq.max-cut=max-cut path.cirq.max-cut=${path.cirq}/${name.cirq.max-cut}/max_cut_cirq.py venv.cirq.max-cut=${name.solvers}_${name.cirq}_${name.cirq.max-cut} name.cplex=cplex path.cplex=${name.solvers}/${name.cplex} name.cplex.mip=mip path.cplex.mip=${path.cplex}/${name.cplex.mip}/mip.py venv.cplex.mip=${name.solvers}_${name.cplex}_${name.cplex.mip} name.qrisp=qrisp path.qrisp=${name.solvers}/${name.qrisp} name.qrisp.vrp=vrp path.qrisp.vrp=${path.qrisp}/${name.qrisp.vrp}/grover.py venv.qrisp.vrp=${name.solvers}_${name.qrisp}_${name.qrisp.vrp} name.qrisp.qubo=qubo path.qrisp.qubo=${path.qrisp}/${name.qrisp.qubo}/qaoa.py venv.qrisp.qubo=${name.solvers}_${name.qrisp}_${name.qrisp.qubo} name.qrisp.sat=sat path.qrisp.sat.grover=${path.qrisp}/${name.qrisp.sat}/grover.py path.qrisp.sat.exact=${path.qrisp}/${name.qrisp.sat}/exact_grover.py venv.qrisp.sat=${name.solvers}_${name.qrisp}_${name.qrisp.sat} name.dwave=dwave path.dwave=${name.solvers}/${name.dwave} name.dwave.qubo=qubo path.dwave.qubo=${path.dwave}/${name.dwave.qubo}/main.py venv.dwave.qubo=${name.solvers}_${name.dwave}_${name.dwave.qubo} # Non OS-specific custom solvers: (solvers that are not part of a framework) name.custom=custom path.custom=${name.solvers}/${name.custom} name.custom.hs-knapsack=hs-knapsack path.custom.hs-knapsack=${path.custom}/${name.custom.hs-knapsack}/knapsack.py venv.custom.hs-knapsack=${name.solvers}_${name.custom}_${name.custom.hs-knapsack} name.custom.lkh=lkh path.custom.lkh=${path.custom}/${name.custom.lkh}/vrp_lkh.py venv.custom.lkh=${name.solvers}_${name.custom}_${name.custom.lkh} name.custom.berger-vrp=berger-vrp name.custom.sharp-sat-bruteforce=sharp-sat-bruteforce path.custom.sharp-sat-bruteforce=${path.custom}/${name.custom.sharp-sat-bruteforce}/exact-solution-counter.py venv.custom.sharp-sat-bruteforce=${name.solvers}_${name.custom}_${name.custom.sharp-sat-bruteforce} name.custom.sharp-sat-ganak=sharp-sat-ganak venv.custom.sharp-sat-ganak=${name.solvers}_${name.custom}_${name.custom.sharp-sat-ganak} # Demonstrators name.demonstrators=demonstrators name.demonstrators.cplex=cplex path.demonstrators.cplex=${name.demonstrators}/${name.demonstrators.cplex} name.demonstrators.cplex.mip=mip-solver path.demonstrators.cplex.mip=${path.demonstrators.cplex}/${name.demonstrators.cplex.mip}/mip-solver.py venv.demonstrators.cplex.mip=${name.demonstrators}_${name.demonstrators.cplex}_${name.demonstrators.cplex.mip} \ No newline at end of file +# default spring profile, correct one will be set during runtime (see ToolboxServerApplication.java) # options: mac, windows, linux spring.profiles.active=linux springdoc.swagger-ui.operationsSorter=alpha springdoc.swagger-ui.tagsSorter=alpha working.directory=jobs examples.directory=examples springdoc.swagger-ui.path=/ # Solvers name.solvers=solvers # Non OS-specific solvers: (typically GAMS and Python) name.gams=gams path.gams=${name.solvers}/${name.gams} name.gams.max-cut=max-cut path.gams.max-cut=${path.gams}/${name.gams.max-cut}/maxcut.gms name.gams.sat=sat path.gams.sat=${path.gams}/${name.gams.sat}/sat.gms name.gams.mip=mip path.gams.mip-to-qubo-reformulation=${path.gams}/${name.gams.mip}/mip_to_qubo_reformulation.gms path.gams.qubo-solve=${path.gams}/${name.gams.mip}/qubo_solve.gms name.qiskit=qiskit path.qiskit=${name.solvers}/${name.qiskit} name.qiskit.knapsack=knapsack path.qiskit.knapsack=${path.qiskit}/${name.qiskit.knapsack}/knapsack_qiskit.py venv.qiskit.knapsack=${name.solvers}_${name.qiskit}_${name.qiskit.knapsack} name.qiskit.max-cut=max-cut path.qiskit.max-cut=${path.qiskit}/${name.qiskit.max-cut}/maxCut_qiskit.py venv.qiskit.max-cut=${name.solvers}_${name.qiskit}_${name.qiskit.max-cut} name.qiskit.qubo=qubo path.qiskit.qubo=${path.qiskit}/${name.qiskit.qubo}/qubo_qiskit.py venv.qiskit.qubo=${name.solvers}_${name.qiskit}_${name.qiskit.qubo} name.cirq=cirq path.cirq=${name.solvers}/${name.cirq} name.cirq.max-cut=max-cut path.cirq.max-cut=${path.cirq}/${name.cirq.max-cut}/max_cut_cirq.py venv.cirq.max-cut=${name.solvers}_${name.cirq}_${name.cirq.max-cut} name.cplex=cplex path.cplex=${name.solvers}/${name.cplex} name.cplex.mip=mip path.cplex.mip=${path.cplex}/${name.cplex.mip}/mip.py venv.cplex.mip=${name.solvers}_${name.cplex}_${name.cplex.mip} name.qrisp=qrisp path.qrisp=${name.solvers}/${name.qrisp} name.qrisp.vrp=vrp path.qrisp.vrp=${path.qrisp}/${name.qrisp.vrp}/grover.py venv.qrisp.vrp=${name.solvers}_${name.qrisp}_${name.qrisp.vrp} name.qrisp.qubo=qubo path.qrisp.qubo=${path.qrisp}/${name.qrisp.qubo}/qaoa.py venv.qrisp.qubo=${name.solvers}_${name.qrisp}_${name.qrisp.qubo} name.qrisp.sat=sat path.qrisp.sat.grover=${path.qrisp}/${name.qrisp.sat}/grover.py path.qrisp.sat.exact=${path.qrisp}/${name.qrisp.sat}/exact_grover.py venv.qrisp.sat=${name.solvers}_${name.qrisp}_${name.qrisp.sat} name.dwave=dwave path.dwave=${name.solvers}/${name.dwave} name.dwave.qubo=qubo path.dwave.qubo=${path.dwave}/${name.dwave.qubo}/main.py venv.dwave.qubo=${name.solvers}_${name.dwave}_${name.dwave.qubo} # Non OS-specific custom solvers: (solvers that are not part of a framework) name.custom=custom path.custom=${name.solvers}/${name.custom} name.custom.hs-knapsack=hs-knapsack path.custom.hs-knapsack=${path.custom}/${name.custom.hs-knapsack}/knapsack.py venv.custom.hs-knapsack=${name.solvers}_${name.custom}_${name.custom.hs-knapsack} name.custom.lkh=lkh path.custom.lkh=${path.custom}/${name.custom.lkh}/vrp_lkh.py venv.custom.lkh=${name.solvers}_${name.custom}_${name.custom.lkh} name.custom.berger-vrp=berger-vrp name.custom.sharp-sat-bruteforce=sharp-sat-bruteforce path.custom.sharp-sat-bruteforce=${path.custom}/${name.custom.sharp-sat-bruteforce}/exact-solution-counter.py venv.custom.sharp-sat-bruteforce=${name.solvers}_${name.custom}_${name.custom.sharp-sat-bruteforce} name.custom.sharp-sat-ganak=sharp-sat-ganak venv.custom.sharp-sat-ganak=${name.solvers}_${name.custom}_${name.custom.sharp-sat-ganak} # Demonstrators name.demonstrators=demonstrators name.demonstrators.cplex=cplex path.demonstrators.cplex=${name.demonstrators}/${name.demonstrators.cplex} name.demonstrators.cplex.mip=mip-solver path.demonstrators.cplex.mip=${path.demonstrators.cplex}/${name.demonstrators.cplex.mip}/mip-solver.py venv.demonstrators.cplex.mip=${name.demonstrators}_${name.demonstrators.cplex}_${name.demonstrators.cplex.mip} \ No newline at end of file diff --git a/src/test/java/edu/kit/provideq/toolbox/api/MipSolversTest.java b/src/test/java/edu/kit/provideq/toolbox/api/MipSolversTest.java index f09887b5..882d77cf 100644 --- a/src/test/java/edu/kit/provideq/toolbox/api/MipSolversTest.java +++ b/src/test/java/edu/kit/provideq/toolbox/api/MipSolversTest.java @@ -44,7 +44,7 @@ Stream provideArguments() { @ParameterizedTest @MethodSource("provideArguments") - void testSatSolver(ProblemSolver solver, String input) { + void testMipSolver(ProblemSolver solver, String input) { var problem = ApiTestHelper.createProblem(client, solver, input, MIP); ApiTestHelper.testSolution(problem); } From beb374b529692dc4b370b8a775c2fdce7be03964 Mon Sep 17 00:00:00 2001 From: Piotr Malkowski Date: Thu, 12 Jun 2025 16:31:44 +0200 Subject: [PATCH 37/37] feat: split translation functionality from solve functionality --- .../gams/mip/gams_to_qubo_reformulation.gms | 446 ++++++++++++++++++ ...ation.gms => mip_to_gms_reformulation.gms} | 2 +- .../toolbox/mip/MipConfiguration.java | 4 +- .../{CplexMip.java => CplexMipSolver.java} | 9 +- .../toolbox/mip/solvers/QuboMipSolver.java | 16 +- .../toolbox/process/GamsProcessRunner.java | 2 +- src/main/resources/application.properties | 2 +- 7 files changed, 464 insertions(+), 17 deletions(-) create mode 100644 solvers/gams/mip/gams_to_qubo_reformulation.gms rename solvers/gams/mip/{mip_to_qubo_reformulation.gms => mip_to_gms_reformulation.gms} (98%) rename src/main/java/edu/kit/provideq/toolbox/mip/solvers/{CplexMip.java => CplexMipSolver.java} (90%) diff --git a/solvers/gams/mip/gams_to_qubo_reformulation.gms b/solvers/gams/mip/gams_to_qubo_reformulation.gms new file mode 100644 index 00000000..915e672b --- /dev/null +++ b/solvers/gams/mip/gams_to_qubo_reformulation.gms @@ -0,0 +1,446 @@ +$offEolCom +$offListing +$eolcom # + +$set modelName %1 +$set modelType %2 +$set direction %3 +$set obj %4 +$set penalty %5 + + +$onEcho > convert.opt +dumpgdx %modelName%.gdx +GDXQuadratic 1 +$offEcho + +option %modelType%=convert; +%modelName%.optfile = 1; +Solve %modelName% use %modelType% %direction% %obj%; # dumping the problem data in a gdx file + +* QUBO Reformulations +EmbeddedCode Python: +import logging +import warnings +from typing import Tuple, Optional + +import numpy as np +import pandas as pd +from gams import transfer as gt + + + +warnings.simplefilter(action='ignore', category=pd.errors.PerformanceWarning) + +is_max = "x" in "%direction%".lower() +gdx_file = r"%modelName%.gdx" +container = gt.Container(gdx_file) +obj_eq_name = container['iobj'].records + +if obj_eq_name is None: + raise Exception("The objective is not defined using a scalar equation. `iobj` in gdx is empty. Quitting.") + +obj_var = container['jobj'].records # fetch the objective variable name +all_vars = container['j'].records # fetches all variable names +raw_a = container['A'].records # A coefficients +eq_data = container['e'].records # fetches equation data + +if raw_a[-raw_a['i'].isin(obj_eq_name['i'].tolist())]['value'].mod(1).sum(axis=0) > 0: # floating point coeffs in objective function is accepted + raise Exception("Reformulation with Non-Integer Coefficients not possible. Quitting.") + +raw_a = raw_a.pivot(index="i", columns="j", values="value").fillna(0) # arranging in a matrix + +if eq_data['lower'].mod(1).sum(axis=0) > 0 or eq_data['upper'].mod(1).sum(axis=0) > 0: + raise Exception("Reformulation with Non-Integer RHS not possible. Quitting.") + +bin_vars = container['jb'].records # fetches binary variable names, if any +int_vars = container['ji'].records # fetches integer variable names, if any +bin_vars = [] if bin_vars is None else bin_vars['j'].to_list() # check if any bin_vars are present +int_vars = [] if int_vars is None else int_vars['j'].to_list() # check if any int_vars are present +obj_var = obj_var['j'].to_list() +all_var_vals = container['x'].records # get all variable values, viz., [level, marginal, lower, upper, scale] + +if len(all_vars) - len(bin_vars) - len(int_vars) != 1: # Continuous variables are not allowed + raise Exception("There are continuous variables. Quitting.") + +obj_eq_name = obj_eq_name['i'].to_list() + +check_quad = container['ANL'].records + +""" +Check if there are any fixed variables in the gdx, i.e., lb=ub=level of any variable. +If such variables exist, separate them from the list of non-fixed vairables and treat them as constanst in the objective function. + +We also need to check if the level of variables are set and handle them separately +""" + +vars_with_lower_bounds = {var.j: var.lower for _, var in all_var_vals.iterrows() if (var.lower > 0) and (var.lower != var.upper)} # would only contain, integer varaibles with lower bound defined +fixed_vars = {var.j: var.level for _, var in all_var_vals.iterrows() if (var.level == var.lower) and (var.level == var.upper)} # check for fixed variables +fixed_and_lower_bounds = {**vars_with_lower_bounds, **fixed_vars} +sum_fixed_obj_var_coeffs = 0 + +if check_quad is not None: + rawquad = container['Q'].records # fetch quadratic terms from the original problem, if any. + + if len(int_vars) != 0: + raise Exception("Quadratic Program with integer variables are not supported.") + + if any(check_quad['j'].isin(fixed_and_lower_bounds.keys())): + raise Exception("Quadratic terms with non-zero variable levels are not supported at the moment.") + +logging.debug("Coefficient matrix: raw_a\n"+raw_a.to_string()) +logging.debug("\nEquation Data: eq_data\n"+eq_data.to_string()) +logging.debug("\nVariable Data: all_var_vals\n"+all_var_vals.to_string()) + +def var_contribution(A: pd.DataFrame, vars: dict, cons: Optional[list] = None) -> np.array: + """ + helper function to calculate the contribution of given variables + in a constraint or set of constraints + + Args: + A: df of coefficients + vars: contributing variables + cons: participating constraints + + Returns: + np.array of Total contribution of all variables for that constraint + """ + cons = slice(None) if cons is None else cons + coeffs_of_vars_in_constraint = A.loc[cons, vars.keys()].to_numpy() + lb_var_levels = np.array(list(vars.values())).reshape((len(vars), 1)) + if coeffs_of_vars_in_constraint.size > 0: + return coeffs_of_vars_in_constraint@lb_var_levels + + return np.array([0]) + +if fixed_and_lower_bounds: # adjust the rhs of equations when level of variables > 0 + logging.info(f"\nList of variables with lower bounds:\n{vars_with_lower_bounds}") + contribution = var_contribution(raw_a, fixed_and_lower_bounds) + eq_data.loc[:,['lower', 'upper']] -= contribution + if fixed_vars: + logging.info(f"\nList of Fixed Variables:\n{fixed_vars}") + # remove the fixed variables from computation + bin_vars = [var for var in bin_vars if var not in fixed_vars] + int_vars = [var for var in int_vars if var not in fixed_vars] + sum_fixed_obj_var_coeffs += np.ndarray.item(var_contribution(raw_a, fixed_vars, cons=obj_eq_name)) + raw_a.drop(fixed_vars, axis=1, inplace=True) # dropping columns from the coefficient matrix + fixed_var_vals = all_var_vals[all_var_vals['j'].isin(fixed_vars)].copy(deep=True) + + logging.debug("\nAfter removing fixed variables and adjusting for non-zero levels: raw_a\n"+raw_a.to_string()) + logging.debug("\nAfter removing fixed variables and adjusting for non-zero levels: eq_data\n"+eq_data.to_string()) + +""" +Check if there exist a row in coefficient matrix with all zero values. This can happen if all vars in a constraint are fixed. +Such row is irrelevant for QUBO and can be dropped out of the matrix and set of constraints. +""" +redundant_cons = list(raw_a[raw_a.apply(abs).sum(axis=1)==0].index) +if len(redundant_cons) > 0: + logging.info(f"\nDropping these redundant constraint: \n{redundant_cons}") + raw_a.drop(redundant_cons, axis=0, inplace=True) + eq_data.drop(eq_data[eq_data['i'].isin(redundant_cons)].index, axis=0, inplace=True) + + +def gen_slacks(var_range: float) -> np.array: + """ + helper function to generate slacks depending on the range of variables or rhs + + Args: + var_range: upper bound of variable + + Returns: + Numpy array containing slack co-efficients + + example: + if var_range=5, then gen_slacks(5) returns [1, 2, 2] + """ + if var_range >= 1e+4: + raise Exception("The Upper bound is greater than or equal to 1e+4, Quitting!") + + power = int(np.log2(var_range)) if var_range > 0 else 0 + bounded_coef = var_range - (2**power - 1) + D_val = [2**i for i in range(power)] + [bounded_coef] + return np.array(D_val) + + +""" +If integer variables exist, convert all integers to binary with '@' as a delimiter of variable names +If Integer variables with lower bound exist, i.e., lb >=1 and lb!=ub, then convert binary variable for that range +If these variables contribute to the objective function, their lower bounds are added as a constant +""" +sum_lower_bound_of_int_vars = 0 +if len(int_vars) != 0: + if vars_with_lower_bounds: + sum_lower_bound_of_int_vars += np.ndarray.item(var_contribution(raw_a, vars_with_lower_bounds, obj_eq_name)) + + int_var_vals = all_var_vals[all_var_vals['j'].isin(int_vars)] + int_to_bin_bounds = {row['j']: gen_slacks(row['upper']-row['lower']) for _,row in int_var_vals.iterrows()} # generate coeffs for converted binary vars + int_bin_vals = pd.DataFrame(columns=['intName', 'binName', 'value']) + for var, bin_bounds in int_to_bin_bounds.items(): + for i in range(len(bin_bounds)): # naming the converted binary vars + new_row = pd.DataFrame({'intName': var, 'binName': f"{var}@_bin{i}", 'value': bin_bounds[i]}, index=[0]) + int_bin_vals = pd.concat([int_bin_vals, new_row], ignore_index=True) if not int_bin_vals.empty else new_row.copy() + + binName_list = int_bin_vals['binName'].to_list() # list of all converted binary variable names + int_bin_name_map = list(int_bin_vals[['intName', 'binName']].itertuples(index=None, name=None)) + logging.info("\nInteger to Binary Mapping: int_bin_vals\n"+int_bin_vals.to_string()) + int_bin_vals = int_bin_vals.pivot(index='intName', columns='binName', values='value').fillna(0) # mapping each binary var to its integer var component + int_bin_vals = int_bin_vals.reindex(labels=int_vars, axis='index') + int_bin_vals = int_bin_vals.reindex(labels=binName_list, axis='columns') + # int_bin_vals.columns = pd.MultiIndex.from_tuples(int_bin_name_map) + + raw_a_int = raw_a[int_vars] + raw_a_int = raw_a_int.dot(int_bin_vals) # updating the "A" coeff matrix with the new coeffs for converted binary vars + raw_a_rest = raw_a[obj_var+bin_vars] + raw_a = pd.concat([raw_a_rest, raw_a_int], axis='columns') # new "A" coeff matrix + logging.info("\nInteger to Binary Mapping: raw_a\n"+raw_a.to_string()) + bin_vars += binName_list # append the list of original binary variables with the list of converted binary variables + +cons = eq_data[-eq_data['i'].isin(obj_eq_name)].reset_index(drop=True) # fetch only the constrainsts and not the objective equation +nvars = len(bin_vars) +nslacks = 0 +obj_var_direction = raw_a[obj_var].loc[obj_eq_name].to_numpy() +obj_var_coeff = raw_a[bin_vars].loc[obj_eq_name].to_numpy() +if obj_var_direction > 0: + obj_var_coeff = -1*obj_var_coeff +obj = np.zeros((nvars, nvars)) +np.fill_diagonal(obj, obj_var_coeff) + +""" +Pre-processing in case special constraints exist +Remove the constraint from the "A" matrix and include the special penalty directly in the objective +Doing so, reduces the number of slack variables used in the final reformulation +special penalty case 1: sum(x_i | 1 <= i <= n) <= 1 => P*sum(x_i*x_j | i < j) +special penalty case 2: x_i + x_j >= 1 => P*(1 - x_i - x_j + x_i*x_j) +""" + +def check_row_entries(df: pd.DataFrame) -> pd.DataFrame: + """ + helper function to filter DataFrame having either 0 or 1 entries in each row. + Args: + df: A Pandas DataFrame. + + Returns: + Filtered DataFrame with rows having either 0 or 1 + """ + row_contains_only_0s_or_1s = df.isin([0, 1]).all(axis=1) + return df[row_contains_only_0s_or_1s] + +# Case 1 implementation +special_cons_case_1_lable = [ele.i for _, ele in cons.iterrows() if ele.upper == 1 and ele.lower != 1] +if special_cons_case_1_lable: + case1_cons = raw_a[bin_vars].loc[special_cons_case_1_lable] + case1_cons = check_row_entries(case1_cons.copy()) + case1_cons_index_lable = list(case1_cons.index) + case1_penalty = case1_cons.to_numpy() + if case1_penalty.size > 0: + case1_penalty = (case1_penalty.T@case1_penalty)/2 + np.fill_diagonal(case1_penalty, np.zeros((1, len(bin_vars)))) + else: # if there are no rows with only 0/1 entries + case1_penalty = np.zeros((nvars, nvars)) + logging.debug(f"\nSpecial constraint case 1:\n{special_cons_case_1_lable}") +else: + case1_cons_index_lable = [] + case1_penalty = np.zeros((nvars, nvars)) + +# Case 2 implementation +special_cons_case_2_lable = [ele.i for _, ele in cons.iterrows() if ele.lower == 1 and ele.upper != 1] +if special_cons_case_2_lable: + case2_cons = raw_a[bin_vars].loc[special_cons_case_2_lable] + case2_cons = check_row_entries(case2_cons.copy()) + case2_cons = case2_cons[case2_cons.sum(axis=1)==2] + case2_cons_index_lable = list(case2_cons.index) + case2_penalty = case2_cons.to_numpy() + if case2_penalty.size > 0: + case2_penalty = (case2_penalty.T@case2_penalty)/2 + case2_diag = np.diag_indices_from(case2_penalty) + case2_penalty[case2_diag] *= -2 + else: # if there are no rows with two 1s in them + case2_penalty = np.zeros((nvars, nvars)) + logging.debug("\nSpecial constraint case 2:\n{special_cons_case_2_lable}") + +else: + case2_cons_index_lable = [] + case2_penalty = np.zeros((nvars, nvars)) + +final_special_cons = case1_cons_index_lable + case2_cons_index_lable +final_special_penalty = case1_penalty + case2_penalty +case2_penalty_offset_factor = len(case2_cons_index_lable) + +P = -1 * %penalty% if is_max else %penalty% # penalty term for classic solvers +P_qpu = %penalty% # penalty term for qpu, since we always going to minimize the qubo using qpu no treatment is required + +obj += P*final_special_penalty + + +cons.drop(cons[cons["i"].isin(final_special_cons)].index, axis=0, inplace=True) +raw_a.drop(final_special_cons, axis=0, inplace=True) + +A_coeff = raw_a.loc[cons['i'], bin_vars] + +quad = None + +def fetch_quadratic_coeff(raw_df: pd.DataFrame) -> np.array: + """ + helper function to convert the original Q matrix of the problem to a symmetric matrix + + Args: + raw_df: Original problem Q data in a pd.DataFrame + + Returns: + Numpy Q matrix + """ + raw_df['value'] /= 2 + mask = raw_df['j_1'].astype(str) == raw_df['j_2'].astype(str) + filtered_quad = raw_df.loc[mask, :].copy() + raw_df = raw_df.loc[~mask].copy() + diag_quad = filtered_quad.reset_index(drop=True) + quad = raw_df.copy(deep=True) + quad['j_1'], quad['j_2'] = raw_df['j_2'], raw_df['j_1'] + quad = pd.concat([raw_df, quad, diag_quad], axis=0) + quad = quad.pivot(index="j_1", columns="j_2", values="value").fillna(0) + quad = quad.reindex(labels=bin_vars, axis='index') + quad = quad.reindex(labels=bin_vars, axis='columns') + return quad.to_numpy() + + +if check_quad is not None: # check if quadratic terms are present in the original problem + logging.debug("\nRaw Q data from GDX: Q\n"+rawquad.to_string()) + rawquad_obj = rawquad[rawquad['i_0'].isin(obj_eq_name)].copy(deep=True) + if len(rawquad_obj.index) != 0: # check if quadratic terms exist in the objective function + rawquad_obj.drop(['i_0'],axis=1,inplace=True) + quad = fetch_quadratic_coeff(raw_df=rawquad_obj) + sum_fixed_obj_var_coeffs /= 2 + + rawquad_cons = rawquad[-rawquad['i_0'].isin(obj_eq_name)] # non-linear constraints without objective equation + if len(rawquad_cons.index) != 0: # non-linear constraints exists + raise Exception("There are non-linear constraints. Quitting.") + + ### Removed the support for quadratic constraints. + # mask = rawquad_cons['j_1'].astype(str) == rawquad_cons['j_2'].astype(str) + # filtered_quad = rawquad_cons.loc[mask, :].reset_index(drop=True) + # if len(filtered_quad.index) == len(rawquad_cons.index): # if pair-wise quadratic terms are not present, i.e., only terms like x1*x1 and not x1*x2 + # filtered_quad.drop(['j_2'], axis=1, inplace=True) + # filtered_quad['value'] /= 2 + # filtered_quad = filtered_quad.pivot(index='i_0', columns='j_1', values='value').fillna(0) + # if len(filtered_quad.columns) != len(bin_vars): + # remain_cols = [var for var in bin_vars if var not in set(filtered_quad.columns)] + # filtered_quad[remain_cols] = 0 + # filtered_quad = filtered_quad[bin_vars] + # A_coeff.loc[filtered_quad.index] += filtered_quad.loc[filtered_quad.index].values + # else: # if pair-wise quadratic terms present, i.e., x1*x2. Quit + # raise Exception("There are non-linear constraints. Quitting.") + +if quad is not None: # add the old quadratic terms/matrix to the new objective + logging.debug("\nUpdate Objective by adding Q: \n"+np.array2string(quad)) + obj += -1*quad if obj_var_direction > 0 else quad + logging.debug("\nNew Q: \n"+np.array2string(obj)) + + +def modify_matrix( + b_vec: np.array, + rhs: float, + slacks: np.array, + A_coeff: pd.DataFrame, + ele: pd.Series, + nslacks: int + ) -> Tuple[np.array, pd.DataFrame, int]: + """ + helper function to update the original "A" matrix of coeffs + + Args: + b_vec: The n*1 vector + rhs : The Right hand side of a constraint + slacks: result of gen_slacks() + A_coeff: "A" matrix + ele: constraint + nslacks: number of slacks + + Returns: + updated b_vec, A_coeff and number of slacks + """ + con_index = ele.i + slack_names = [f'slack_{con_index}_{i}' for i in range(nslacks+1, nslacks + len(slacks) + 1)] + A_coeff[slack_names] = 0 + A_coeff.loc[con_index, slack_names] = slacks + nslacks += len(slacks) + return np.append(b_vec, [rhs]), A_coeff, nslacks + +def get_lhs_bounds(ele: pd.DataFrame) -> Tuple[float, float]: + """ + helper function to find the bounds of a constraint + + Args: + ele: The coefficents of the constraint + + Returns: + lower_bound, upper_bound + """ + return ele[ele < 0].sum(), ele[ele > 0].sum() + +b_vec = np.array([]) +logging.info("\nFinal Cons: \n"+cons.to_string()) +for _, ele in cons.iterrows(): + + if ele.upper == ele.lower: # equal-to type constraint + rhs = ele.lower + lhs_min_lb, lhs_max_ub = get_lhs_bounds(A_coeff.loc[ele.i]) + if (rhs - lhs_min_lb) < 0 or (lhs_max_ub - rhs) < 0: + raise Exception(f"Constraint is infeasible: {ele.i}") + else: + b_vec = np.append(b_vec, [rhs]) + slacks = [] # do not introduce slacks for equality type constraints + + elif ele.upper == np.inf: # greater than type constraint + rhs = ele.lower + _, lhs_max_ub = get_lhs_bounds(A_coeff.loc[ele.i]) + slacks_range = lhs_max_ub - rhs + if slacks_range > 0: + slacks = -1*gen_slacks(slacks_range) + b_vec, A_coeff, nslacks = modify_matrix(b_vec, rhs, slacks, A_coeff, ele, nslacks) + elif slacks_range == 0: + b_vec = np.append(b_vec, [rhs]) + slacks = [] + else: + raise Exception(f"Constraint is infeasible: {ele.i}") + + else: # less-than type constraint + rhs = ele.upper + lhs_min_lb, _ = get_lhs_bounds(A_coeff.loc[ele.i]) + slacks_range = rhs - lhs_min_lb + if slacks_range > 0: + slacks = gen_slacks(slacks_range) + b_vec, A_coeff, nslacks = modify_matrix(b_vec, rhs, slacks, A_coeff, ele, nslacks) + elif slacks_range == 0: + b_vec = np.append(b_vec, [rhs]) + slacks = [] + else: + raise Exception(f"Constraint is infeasible: {ele.i}") + + +logging_a_mat = A_coeff.unstack().reset_index() +logging_a_mat = logging_a_mat[logging_a_mat[0]!= 0] +logging.debug("\nFinal coefficient matrix: \n"+logging_a_mat.to_string()) +logging.debug(f"\nFinal RHS: \n{b_vec}") +logging.debug(f"Constant RHS term: {b_vec.T@b_vec}") +logging.debug(f"Case 2 Offset Penalty Factor: {case2_penalty_offset_factor}") +logging.debug(f"Fixed Variable contribution to Objective Function: {sum_fixed_obj_var_coeffs}") +logging.debug(f"Integer variable lower bound contribution: {sum_lower_bound_of_int_vars}") + +# A matrix and b_vec are available. Now, for penalization: $(A.X - B)^{2}$ = $(A.X - B)^{T} * (A.X - B)$ +a_mat = A_coeff.to_numpy() +nvars += nslacks # increment the total number of variables by total number of slack variable used +X_diag = np.zeros((nvars, nvars)) +np.fill_diagonal(X_diag,-b_vec.T@a_mat) +new_x = a_mat.T@a_mat + 2*X_diag + +newobj = np.zeros((nvars,nvars)) +newobj[:len(bin_vars), :len(bin_vars)] = obj # define the new objective: Q for the qubo + + +np.savetxt(fr'%modelName%_p{P:.0f}_c{const:.0f}.csv', Q, delimiter=",") +logging.debug(f"\nExported Q matrix as >%modelName%_p{P:.0f}_c{const:.0f}.csv<\n") + + +endEmbeddedCode + +abort$execError 'Error occured in Reformulation!'; \ No newline at end of file diff --git a/solvers/gams/mip/mip_to_qubo_reformulation.gms b/solvers/gams/mip/mip_to_gms_reformulation.gms similarity index 98% rename from solvers/gams/mip/mip_to_qubo_reformulation.gms rename to solvers/gams/mip/mip_to_gms_reformulation.gms index 0153bb29..6019d1a7 100644 --- a/solvers/gams/mip/mip_to_qubo_reformulation.gms +++ b/solvers/gams/mip/mip_to_gms_reformulation.gms @@ -48,7 +48,7 @@ else: # Create the replacement line using the matched groups penalty = os.environ.get("GAMSPEN", "100") # Default to 100 if not set -new_line = f'$batinclude %IDIR% {model_name} {model_type} {gams_sense} {obj_var} {penalty} -logOn=2' +new_line = f'$batinclude %IDIR% {model_name} {model_type} {gams_sense} {obj_var} {penalty}' # Replace the matching "Solve ..." line with the new line (only the first occurrence) new_content = pattern.sub(new_line, content, count=1) diff --git a/src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java b/src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java index 60dc88dd..71bc7172 100644 --- a/src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java +++ b/src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java @@ -5,7 +5,7 @@ import edu.kit.provideq.toolbox.meta.Problem; import edu.kit.provideq.toolbox.meta.ProblemManager; import edu.kit.provideq.toolbox.meta.ProblemType; -import edu.kit.provideq.toolbox.mip.solvers.CplexMip; +import edu.kit.provideq.toolbox.mip.solvers.CplexMipSolver; import edu.kit.provideq.toolbox.mip.solvers.QuboMipSolver; import java.io.IOException; import java.util.Objects; @@ -32,7 +32,7 @@ public class MipConfiguration { @Bean ProblemManager getMipManager( - CplexMip cplexMip, + CplexMipSolver cplexMip, QuboMipSolver quboMipSolver, ResourceProvider resourceProvider ) { diff --git a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/CplexMip.java b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/CplexMipSolver.java similarity index 90% rename from src/main/java/edu/kit/provideq/toolbox/mip/solvers/CplexMip.java rename to src/main/java/edu/kit/provideq/toolbox/mip/solvers/CplexMipSolver.java index 7569dae6..eca7a3f9 100644 --- a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/CplexMip.java +++ b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/CplexMipSolver.java @@ -16,14 +16,14 @@ * {@link MipConfiguration#MIP} solver using IBM cplex. */ @Component -public class CplexMip extends MipSolver { +public class CplexMipSolver extends MipSolver { private final String scriptPath; private final String venv; private final ApplicationContext context; @Autowired - public CplexMip( + public CplexMipSolver( @Value("${path.cplex.mip}") String scriptPath, @Value("${venv.cplex.mip}") String venv, ApplicationContext context) { @@ -34,12 +34,13 @@ public CplexMip( @Override public String getName() { - return "IBM Cplex Solver for MIP"; + return "Cplex"; } @Override public String getDescription() { - return "This solver uses IBM Cplex Solver to solve MIP problems"; + return "This solver uses IBM Cplex Solver to solve MIP problems. " + + "It is called via Python wrapper without any additional parameters"; } @Override diff --git a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java index a4f9df76..e48c2c87 100644 --- a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java +++ b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java @@ -19,18 +19,18 @@ */ @Component public class QuboMipSolver extends MipSolver { - private final String scriptPath; - private final String dependencyScriptPath; + private final String mipToGamsReformulatorPath; + private final String solverScriptPath; private final ApplicationContext context; private static final int PENALTY = 100; @Autowired public QuboMipSolver( - @Value("${path.gams.mip-to-qubo-reformulation}") String mipToQuboReformulatorPath, - @Value("${path.gams.qubo-solve}") String solverScriptPath, + @Value("${path.gams.mip-to-gams-reformulation}") String mipToGamsReformulatorPath, + @Value("${path.gams.gams-to-qubo-reformulation}") String solverScriptPath, ApplicationContext context) { - this.scriptPath = mipToQuboReformulatorPath; - this.dependencyScriptPath = solverScriptPath; + this.mipToGamsReformulatorPath = mipToGamsReformulatorPath; + this.solverScriptPath = solverScriptPath; this.context = context; } @@ -55,7 +55,7 @@ public Mono> solve( // Run GAMS via console context - .getBean(GamsProcessRunner.class, scriptPath) + .getBean(GamsProcessRunner.class, mipToGamsReformulatorPath) .withArguments( "--INPUT=" + ProcessRunner.INPUT_FILE_PATH, "--PENALTY=" + PENALTY @@ -69,7 +69,7 @@ public Mono> solve( ProcessResult processResult = context .getBean(GamsProcessRunner.class, relativePath) .withArguments( - "--IDIR=" + new File(dependencyScriptPath).getAbsolutePath(), + "--IDIR=" + new File(solverScriptPath).getAbsolutePath(), "--SOLOUTPUT=" + ProcessRunner.OUTPUT_FILE_PATH ) .readOutputFile("output.txt") diff --git a/src/main/java/edu/kit/provideq/toolbox/process/GamsProcessRunner.java b/src/main/java/edu/kit/provideq/toolbox/process/GamsProcessRunner.java index a00921b3..30aa0ada 100644 --- a/src/main/java/edu/kit/provideq/toolbox/process/GamsProcessRunner.java +++ b/src/main/java/edu/kit/provideq/toolbox/process/GamsProcessRunner.java @@ -26,7 +26,7 @@ public class GamsProcessRunner extends ProcessRunner { /** * The name of the file to call for running a GAMS script. */ - private static final String GAMS_EXECUTABLE_NAME = "gams"; + private static final String GAMS_EXECUTABLE_NAME = "/home/piotr/gams49.4_linux_x64_64_sfx/gams"; /** * Creates a process runner for a GAMS task. diff --git a/src/main/resources/application.properties b/src/main/resources/application.properties index 63bdcae2..cba6146a 100644 --- a/src/main/resources/application.properties +++ b/src/main/resources/application.properties @@ -1 +1 @@ -# default spring profile, correct one will be set during runtime (see ToolboxServerApplication.java) # options: mac, windows, linux spring.profiles.active=linux springdoc.swagger-ui.operationsSorter=alpha springdoc.swagger-ui.tagsSorter=alpha working.directory=jobs examples.directory=examples springdoc.swagger-ui.path=/ # Solvers name.solvers=solvers # Non OS-specific solvers: (typically GAMS and Python) name.gams=gams path.gams=${name.solvers}/${name.gams} name.gams.max-cut=max-cut path.gams.max-cut=${path.gams}/${name.gams.max-cut}/maxcut.gms name.gams.sat=sat path.gams.sat=${path.gams}/${name.gams.sat}/sat.gms name.gams.mip=mip path.gams.mip-to-qubo-reformulation=${path.gams}/${name.gams.mip}/mip_to_qubo_reformulation.gms path.gams.qubo-solve=${path.gams}/${name.gams.mip}/qubo_solve.gms name.qiskit=qiskit path.qiskit=${name.solvers}/${name.qiskit} name.qiskit.knapsack=knapsack path.qiskit.knapsack=${path.qiskit}/${name.qiskit.knapsack}/knapsack_qiskit.py venv.qiskit.knapsack=${name.solvers}_${name.qiskit}_${name.qiskit.knapsack} name.qiskit.max-cut=max-cut path.qiskit.max-cut=${path.qiskit}/${name.qiskit.max-cut}/maxCut_qiskit.py venv.qiskit.max-cut=${name.solvers}_${name.qiskit}_${name.qiskit.max-cut} name.qiskit.qubo=qubo path.qiskit.qubo=${path.qiskit}/${name.qiskit.qubo}/qubo_qiskit.py venv.qiskit.qubo=${name.solvers}_${name.qiskit}_${name.qiskit.qubo} name.cirq=cirq path.cirq=${name.solvers}/${name.cirq} name.cirq.max-cut=max-cut path.cirq.max-cut=${path.cirq}/${name.cirq.max-cut}/max_cut_cirq.py venv.cirq.max-cut=${name.solvers}_${name.cirq}_${name.cirq.max-cut} name.cplex=cplex path.cplex=${name.solvers}/${name.cplex} name.cplex.mip=mip path.cplex.mip=${path.cplex}/${name.cplex.mip}/mip.py venv.cplex.mip=${name.solvers}_${name.cplex}_${name.cplex.mip} name.qrisp=qrisp path.qrisp=${name.solvers}/${name.qrisp} name.qrisp.vrp=vrp path.qrisp.vrp=${path.qrisp}/${name.qrisp.vrp}/grover.py venv.qrisp.vrp=${name.solvers}_${name.qrisp}_${name.qrisp.vrp} name.qrisp.qubo=qubo path.qrisp.qubo=${path.qrisp}/${name.qrisp.qubo}/qaoa.py venv.qrisp.qubo=${name.solvers}_${name.qrisp}_${name.qrisp.qubo} name.qrisp.sat=sat path.qrisp.sat.grover=${path.qrisp}/${name.qrisp.sat}/grover.py path.qrisp.sat.exact=${path.qrisp}/${name.qrisp.sat}/exact_grover.py venv.qrisp.sat=${name.solvers}_${name.qrisp}_${name.qrisp.sat} name.dwave=dwave path.dwave=${name.solvers}/${name.dwave} name.dwave.qubo=qubo path.dwave.qubo=${path.dwave}/${name.dwave.qubo}/main.py venv.dwave.qubo=${name.solvers}_${name.dwave}_${name.dwave.qubo} # Non OS-specific custom solvers: (solvers that are not part of a framework) name.custom=custom path.custom=${name.solvers}/${name.custom} name.custom.hs-knapsack=hs-knapsack path.custom.hs-knapsack=${path.custom}/${name.custom.hs-knapsack}/knapsack.py venv.custom.hs-knapsack=${name.solvers}_${name.custom}_${name.custom.hs-knapsack} name.custom.lkh=lkh path.custom.lkh=${path.custom}/${name.custom.lkh}/vrp_lkh.py venv.custom.lkh=${name.solvers}_${name.custom}_${name.custom.lkh} name.custom.berger-vrp=berger-vrp name.custom.sharp-sat-bruteforce=sharp-sat-bruteforce path.custom.sharp-sat-bruteforce=${path.custom}/${name.custom.sharp-sat-bruteforce}/exact-solution-counter.py venv.custom.sharp-sat-bruteforce=${name.solvers}_${name.custom}_${name.custom.sharp-sat-bruteforce} name.custom.sharp-sat-ganak=sharp-sat-ganak venv.custom.sharp-sat-ganak=${name.solvers}_${name.custom}_${name.custom.sharp-sat-ganak} # Demonstrators name.demonstrators=demonstrators name.demonstrators.cplex=cplex path.demonstrators.cplex=${name.demonstrators}/${name.demonstrators.cplex} name.demonstrators.cplex.mip=mip-solver path.demonstrators.cplex.mip=${path.demonstrators.cplex}/${name.demonstrators.cplex.mip}/mip-solver.py venv.demonstrators.cplex.mip=${name.demonstrators}_${name.demonstrators.cplex}_${name.demonstrators.cplex.mip} \ No newline at end of file +# default spring profile, correct one will be set during runtime (see ToolboxServerApplication.java) # options: mac, windows, linux spring.profiles.active=linux springdoc.swagger-ui.operationsSorter=alpha springdoc.swagger-ui.tagsSorter=alpha working.directory=jobs examples.directory=examples springdoc.swagger-ui.path=/ # Solvers name.solvers=solvers # Non OS-specific solvers: (typically GAMS and Python) name.gams=gams path.gams=${name.solvers}/${name.gams} name.gams.max-cut=max-cut path.gams.max-cut=${path.gams}/${name.gams.max-cut}/maxcut.gms name.gams.sat=sat path.gams.sat=${path.gams}/${name.gams.sat}/sat.gms name.gams.mip=mip path.gams.mip-to-gams-reformulation=${path.gams}/${name.gams.mip}/mip_to_gms_reformulation.gms path.gams.gams-to-qubo-reformulation=${path.gams}/${name.gams.mip}/gams_to_qubo_reformulation.gms name.qiskit=qiskit path.qiskit=${name.solvers}/${name.qiskit} name.qiskit.knapsack=knapsack path.qiskit.knapsack=${path.qiskit}/${name.qiskit.knapsack}/knapsack_qiskit.py venv.qiskit.knapsack=${name.solvers}_${name.qiskit}_${name.qiskit.knapsack} name.qiskit.max-cut=max-cut path.qiskit.max-cut=${path.qiskit}/${name.qiskit.max-cut}/maxCut_qiskit.py venv.qiskit.max-cut=${name.solvers}_${name.qiskit}_${name.qiskit.max-cut} name.qiskit.qubo=qubo path.qiskit.qubo=${path.qiskit}/${name.qiskit.qubo}/qubo_qiskit.py venv.qiskit.qubo=${name.solvers}_${name.qiskit}_${name.qiskit.qubo} name.cirq=cirq path.cirq=${name.solvers}/${name.cirq} name.cirq.max-cut=max-cut path.cirq.max-cut=${path.cirq}/${name.cirq.max-cut}/max_cut_cirq.py venv.cirq.max-cut=${name.solvers}_${name.cirq}_${name.cirq.max-cut} name.cplex=cplex path.cplex=${name.solvers}/${name.cplex} name.cplex.mip=mip path.cplex.mip=${path.cplex}/${name.cplex.mip}/mip.py venv.cplex.mip=${name.solvers}_${name.cplex}_${name.cplex.mip} name.qrisp=qrisp path.qrisp=${name.solvers}/${name.qrisp} name.qrisp.vrp=vrp path.qrisp.vrp=${path.qrisp}/${name.qrisp.vrp}/grover.py venv.qrisp.vrp=${name.solvers}_${name.qrisp}_${name.qrisp.vrp} name.qrisp.qubo=qubo path.qrisp.qubo=${path.qrisp}/${name.qrisp.qubo}/qaoa.py venv.qrisp.qubo=${name.solvers}_${name.qrisp}_${name.qrisp.qubo} name.qrisp.sat=sat path.qrisp.sat.grover=${path.qrisp}/${name.qrisp.sat}/grover.py path.qrisp.sat.exact=${path.qrisp}/${name.qrisp.sat}/exact_grover.py venv.qrisp.sat=${name.solvers}_${name.qrisp}_${name.qrisp.sat} name.dwave=dwave path.dwave=${name.solvers}/${name.dwave} name.dwave.qubo=qubo path.dwave.qubo=${path.dwave}/${name.dwave.qubo}/main.py venv.dwave.qubo=${name.solvers}_${name.dwave}_${name.dwave.qubo} # Non OS-specific custom solvers: (solvers that are not part of a framework) name.custom=custom path.custom=${name.solvers}/${name.custom} name.custom.hs-knapsack=hs-knapsack path.custom.hs-knapsack=${path.custom}/${name.custom.hs-knapsack}/knapsack.py venv.custom.hs-knapsack=${name.solvers}_${name.custom}_${name.custom.hs-knapsack} name.custom.lkh=lkh path.custom.lkh=${path.custom}/${name.custom.lkh}/vrp_lkh.py venv.custom.lkh=${name.solvers}_${name.custom}_${name.custom.lkh} name.custom.berger-vrp=berger-vrp name.custom.sharp-sat-bruteforce=sharp-sat-bruteforce path.custom.sharp-sat-bruteforce=${path.custom}/${name.custom.sharp-sat-bruteforce}/exact-solution-counter.py venv.custom.sharp-sat-bruteforce=${name.solvers}_${name.custom}_${name.custom.sharp-sat-bruteforce} name.custom.sharp-sat-ganak=sharp-sat-ganak venv.custom.sharp-sat-ganak=${name.solvers}_${name.custom}_${name.custom.sharp-sat-ganak} # Demonstrators name.demonstrators=demonstrators name.demonstrators.cplex=cplex path.demonstrators.cplex=${name.demonstrators}/${name.demonstrators.cplex} name.demonstrators.cplex.mip=mip-solver path.demonstrators.cplex.mip=${path.demonstrators.cplex}/${name.demonstrators.cplex.mip}/mip-solver.py venv.demonstrators.cplex.mip=${name.demonstrators}_${name.demonstrators.cplex}_${name.demonstrators.cplex.mip} \ No newline at end of file