Fixed Qiskit QAOA solver (reimplemtented for Qiskit-Optimization-0.6)…

…, refactored Qrisp QAOA (was in the wrong directory)

solvers/qiskit/qubo/qp_converter/QpConverter.py

@@ -0,0 +1,125 @@
+import re
+
+from qiskit_optimization import QuadraticProgram
+
+def parse_lp_file(file_path):
+    # Read the LP file
+    with open(file_path, 'r') as file:
+        lines = file.readlines()
+
+    # Initialize components for parsing
+    objective_sense = 'minimize'  # default
+    objective = ''
+    constraints = []
+    variable_types = {}
+    in_objective_section = False
+    in_constraints_section = False
+    in_binary_section = False
+
+    # Parsing logic
+    for line in lines:
+        line = line.strip()
+
+        # Objective function
+        if line.lower().startswith('minimize') or line.lower().startswith('maximize'):
+            in_objective_section = True
+            objective_sense = line.split()[0].lower()
+            continue
+
+        if in_objective_section and (line.lower().startswith('subject to') or line.lower().startswith('binary')):
+            in_objective_section = False
+            continue
+
+        if in_objective_section:
+            objective += line
+            continue
+
+        # Constraints (no constraints in this example, but you may extend it later)
+        if line.lower().startswith('subject to'):
+            in_constraints_section = True
+            continue
+        elif in_constraints_section and line.lower().startswith('binary'):
+            in_constraints_section = False
+
+        # Binary variables
+        if line.lower().startswith('binary'):
+            in_binary_section = True
+            continue
+
+        if in_binary_section:
+            variables = line.split()
+            for var in variables:
+                variable_types[var] = 'binary'
+
+    return objective_sense, objective, constraints, variable_types
+
+
+def parse_objective(objective_str):
+    """ Parse the objective function to extract linear and quadratic coefficients """
+    # Remove the 'obj:' prefix
+    objective_str = re.sub(r'^obj:\s*', '', objective_str)
+
+    # Extract and remove the division factor (e.g., "/ 2", "/ 3", "/ 10")
+    division_factor_match = re.search(r'/\s*([\d\.]+)', objective_str)
+    division_factor = float(division_factor_match.group(1)) if division_factor_match else 1.0
+    objective_str = re.sub(r'/\s*[\d\.]+', '', objective_str)
+
+    # Remove square brackets
+    objective_str = objective_str.replace('[', '').replace(']', '')
+
+    # Extract terms using regular expression
+    term_pattern = re.compile(r'([+-]?\s*\d*\.?\d+(?:e[+-]?\d+)?)\s*\*?\s*([\w\d]+)\s*(?:\*\s*([\w\d]+))?')
+    terms = term_pattern.findall(objective_str)
+    linear = {}
+    quadratic = {}
+
+    # Parse each term
+    for coeff_str, var1, var2 in terms:
+        coeff = float(coeff_str.replace(' ', ''))  # Remove any spaces in coefficient
+
+        if var2:  # Quadratic term
+            if (var1, var2) in quadratic:
+                quadratic[(var1, var2)] += coeff
+            else:
+                quadratic[(var1, var2)] = coeff
+        else:  # Linear term
+            if var1 in linear:
+                linear[var1] += coeff
+            else:
+                linear[var1] = coeff
+
+    # Apply division factor
+    if division_factor != 1.0 and division_factor != 0:
+        for key in linear:
+            linear[key] /= division_factor
+        for key in quadratic:
+            quadratic[key] /= division_factor
+
+    return linear, quadratic
+
+
+def convert_to_quadratic_program(lp_file_path):
+    # Parse the LP file
+    objective_sense, objective, constraints, variable_types = parse_lp_file(lp_file_path)
+
+    # Parse objective components
+    linear, quadratic = parse_objective(objective)
+
+    # Create an instance of QuadraticProgram
+    qp = QuadraticProgram()
+
+    # Add variables
+    for var in variable_types:
+        if variable_types[var] == 'binary':
+            qp.binary_var(var)
+
+    # Set the objective sense
+    if objective_sense == 'minimize':
+        qp.minimize(linear=linear, quadratic=quadratic)
+    else:
+        qp.maximize(linear=linear, quadratic=quadratic)
+
+    # Add constraints (empty in this example)
+    # You can extend this logic to handle any constraints
+
+    return qp

solvers/qiskit/qubo/qubo_qiskit.py

@@ -1,9 +1,9 @@
 import sys
+from qp_converter import QpConverter
 
-from qiskit.algorithms.minimum_eigensolvers import QAOA
-from qiskit.algorithms.optimizers import COBYLA
 from qiskit.primitives import Sampler
-from qiskit_optimization import QuadraticProgram
+from qiskit_algorithms import QAOA
+from qiskit_algorithms.optimizers import COBYLA
 from qiskit_optimization.algorithms import MinimumEigenOptimizer
 
 if len(sys.argv) != 3:
@@ -12,16 +12,20 @@
 input_path = sys.argv[1]
 output_path = sys.argv[2]
 
-qp = QuadraticProgram()
-qp.read_from_lp_file(input_path)
-
-qaoa_mes = QAOA(Sampler(), optimizer=COBYLA(), initial_point=[0.0, 1.0])
+qubo = QpConverter.convert_to_quadratic_program(input_path)
+# parse LP file:
+print(qubo.prettyprint())
 
+# TODO: Sampler() has to be replaces with StatevectorSampler() in newer versions. 
+# (currently not yet supported by qiskit-optimization)
+# TODO: add a dedicated mixer
+qaoa_mes = QAOA(sampler=Sampler(), optimizer=COBYLA())
 qaoa = MinimumEigenOptimizer(qaoa_mes)
+qaoa.solve(qubo)
 
-qaoa_result = qaoa.solve(qp)
+qaoa_result = qaoa.solve(qubo)
 print(qaoa_result.prettyprint())
 
 f = open(output_path, 'w')
 f.write(qaoa_result.prettyprint())
-f.close()
+f.close()

solvers/qrisp/vrp/qaoa.py → solvers/qrisp/qubo/qaoa.py

@@ -12,7 +12,7 @@
     create_QUBO_cost_operator,
     def_backend,
 )
-from qrisp_solver.permutation import create_perm_specifiers, eval_perm_old
+from qrisp_qubo_util.permutation import create_perm_specifiers, eval_perm_old
 
 from qrisp import QuantumArray, QuantumFloat, QuantumVariable, cyclic_shift, h, x
 

solvers/qrisp/vrp/grover.py

@@ -4,9 +4,9 @@
 
 import numpy as np
 from qrisp.grover import grovers_alg
-from qrisp_solver.oracle import eval_distance_threshold
-from qrisp_solver.permutation import create_perm_specifiers, eval_perm
-from qrisp_solver.vrp import calc_paths, normalize_vrp
+from qrisp_vrp_util.oracle import eval_distance_threshold
+from qrisp_vrp_util.permutation import create_perm_specifiers, eval_perm
+from qrisp_vrp_util.vrp import calc_paths, normalize_vrp
 from tsplib95 import load
 from tsplib95.models import StandardProblem
 

solvers/qrisp/vrp/qrisp_vrp_util/permutation.py

@@ -0,0 +1,93 @@
+import numpy as np
+from qrisp import QuantumArray, QuantumFloat
+from qrisp.core import demux
+from qrisp.environments import invert
+
+
+# Create a function that generates a state of superposition of all permutations
+def swap_to_front(qa, index):
+    with invert():
+        # The keyword ctrl_method = "gray_pt" allows the controlled swaps to be synthesized
+        # using Margolus gates. These gates perform the same operation as a regular Toffoli
+        # but add a different phase for each input. This phase will not matter though,
+        # since it will be reverted once the ancilla values of the oracle are uncomputed.
+        demux(qa[0], index, qa, permit_mismatching_size=True)
+
+
+def eval_perm(perm_specifiers, city_amount, qa=None):
+    N = len(perm_specifiers)
+
+    # To filter out the cyclic permutations, we impose that the first city is always city 0
+    # We will have to consider this assumption later when calculating the route distance
+    # by manually adding the trip distance of the first trip (from city 0) and the
+    # last trip (to city 0)
+    if qa is None:
+        qa = QuantumArray(
+            QuantumFloat(int(np.ceil(np.log2(city_amount)))), city_amount - 1
+        )
+
+    for i in range(city_amount - 1):
+        qa[i] += i + 1
+
+    for i in range(N):
+        swap_to_front(qa[i:], perm_specifiers[i])
+
+    return qa
+
+
+def eval_perm_backward(perm_specifiers, city_amount, qa=None):
+    N = len(perm_specifiers)
+
+    # To filter out the cyclic permutations, we impose that the first city is always city 0
+    # We will have to consider this assumption later when calculating the route distance
+    # by manually adding the trip distance of the first trip (from city 0) and the
+    # last trip (to city 0)
+
+    for i in reversed(range(N)):
+        demux(qa[i], perm_specifiers[i], qa[i:], permit_mismatching_size=True)
+
+    for i in range(city_amount - 1):
+        qa[i] -= i + 1
+
+    return qa
+
+
+def eval_perm_old(perm_specifiers, city_amount, qa=None):
+    N = len(perm_specifiers)
+
+    # To filter out the cyclic permutations, we impose that the first city is always city 0
+    # We will have to consider this assumption later when calculating the route distance
+    # by manually adding the trip distance of the first trip (from city 0) and the
+    # last trip (to city 0)
+    if qa is None:
+        qa = QuantumArray(QuantumFloat(int(np.ceil(np.log2(city_amount)))), city_amount)
+
+    add = np.arange(0, city_amount)
+
+    for i in range(city_amount):
+        qa[i] += int(add[i])
+
+    for i in range(N):
+        swap_to_front(qa[i:], perm_specifiers[i])
+
+    return qa
+
+
+# Create function that returns QuantumFloats specifying the permutations (these will be in uniform superposition)
+def create_perm_specifiers(city_amount, init_seq=None) -> list[QuantumFloat]:
+    perm_specifiers = []
+
+    for i in range(city_amount - 1):
+        qf_size = int(np.ceil(np.log2(city_amount - i)))
+
+        if i == 0:
+            continue
+
+        temp_qf = QuantumFloat(qf_size)
+
+        if not init_seq is None:
+            temp_qf[:] = init_seq[i - 1]
+
+        perm_specifiers.append(temp_qf)
+
+    return perm_specifiers

src/main/java/edu/kit/provideq/toolbox/qubo/QuboConfiguration.java

@@ -48,8 +48,8 @@ ProblemManager<String, String> getQuboManager(
   private Set<Problem<String, String>> loadExampleProblems(ResourceProvider resourceProvider) {
     try {
       var problemInputStream = Objects.requireNonNull(
-          getClass().getResourceAsStream("linear-problem.txt"),
-          "linear-problem example for QUBO is unavailable!"
+          getClass().getResourceAsStream("quadratic-problem.txt"),
+          "quadratic-problem example for QUBO is unavailable!"
       );
       var problem = new Problem<>(QUBO);
       problem.setInput(resourceProvider.readStream(problemInputStream));

src/main/java/edu/kit/provideq/toolbox/qubo/solvers/QrispQuboSolver.java

@@ -20,7 +20,7 @@ public class QrispQuboSolver extends QuboSolver {
 
   @Autowired
   public QrispQuboSolver(
-      @Value("${qrisp.directory.vrp}") String vrpPath,
+      @Value("${qrisp.directory.qubo}") String vrpPath,
       ApplicationContext context) {
     this.vrpPath = vrpPath;
     this.context = context;
@@ -42,14 +42,14 @@ public Mono<Solution<String>> solve(
     // it is used to prevent denial of service issues for large simulations.
     // Default value is 4, higher values are possible but might take much longer to simulate.
     // TODO: allow user to pass a custom gate size as a solver setting
-    int maxNumberOfCities = 4;
+    int maxNumberOfVariables = 4;
 
     var processResult = context.getBean(
             PythonProcessRunner.class,
             vrpPath,
             "qaoa.py",
             new String[] {"%1$s", "--output-file", "%2$s",
-                "--size-gate", String.valueOf(maxNumberOfCities)}
+                "--size-gate", String.valueOf(maxNumberOfVariables)}
         )
         .problemFileName("problem.lp")
         .solutionFileName("problem.bin")

src/main/resources/application.properties

@@ -1,6 +1,6 @@
 # default spring profile, correct one will be set during runtime (see ToolboxServerApplication.java)
 # options: mac, windows, linux
-spring.profiles.active=linux
+spring.profiles.active=mac
 
 working.directory=jobs
 examples.directory=examples
@@ -22,6 +22,7 @@ cirq.directory.max-cut=${cirq.directory}/max-cut
 
 qrisp.directory=${solvers.directory}/qrisp
 qrisp.directory.vrp=${qrisp.directory}/vrp
+qrisp.directory.qubo=${qrisp.directory}/qubo
 
 dwave.directory=${solvers.directory}/dwave
 dwave.directory.qubo=${dwave.directory}/qubo

src/main/resources/edu/kit/provideq/toolbox/qubo/quadratic-problem.txt

@@ -0,0 +1,7 @@
+Minimize
+  obj: 0 + [ -32.984845004941285 x1 * x1 + 32.984845004941285 x1 * x2 + 32.984845004941285 x1 * x3 + 8.246211251235321 x1 * x4 - 32.984845004941285 x2 * x2 + 8.246211251235321 x2 * x3 + 32.984845004941285 x2 * x4 - 32.984845004941285 x3 * x3 + 32.984845004941285 x3 * x4 - 32.984845004941285 x4 * x4 ] / 2
+
+Subject To
+
+Binary
+x1 x2 x3 x4