diff --git a/AGV_quantum/LinearProg.py b/AGV_quantum/LinearProg.py index 4c3e27e..a908c69 100644 --- a/AGV_quantum/LinearProg.py +++ b/AGV_quantum/LinearProg.py @@ -1,5 +1,5 @@ import dimod -from cpp_pyqubo import Binary, Constraint, Placeholder +from cpp_pyqubo import Constraint, Placeholder from pyqubo import LogEncInteger from pyqubo import Binary @@ -48,35 +48,35 @@ def _to_bqm_qubo_ising(self, pdict=None): ) H = 0 ind = 0 - vars = [] - for (lb, ub) in self.bounds: + model_vars = [] + for (lb, ub) in bounds: if lb == 0 and ub == 1: - vars.append(Binary(f"x_{ind}")) + model_vars.append(Binary(f"x_{ind}")) else: - vars.append(LogEncInteger(f"x_{ind}", (lb, ub))) + model_vars.append(LogEncInteger(f"x_{ind}", (lb, ub))) ind += 1 - pyqubo_obj = sum(var * coef for var, coef in zip(vars, c) if coef != 0) + pyqubo_obj = sum(var * coef for var, coef in zip(model_vars, c) if coef != 0) H += Placeholder("obj") * pyqubo_obj num_eq = 0 if A_eq is not None: - for i in range(len(A_eq)): + for i, _ in enumerate(A_eq): expr = sum( - A_eq[i][j] * vars[j] for j in range(self.nvars) if A_eq[i][j] != 0 + A_eq[i][j] * model_vars[j] for j in range(self.nvars) if A_eq[i][j] != 0 ) expr -= b_eq[i] H += Constraint(Placeholder(f"eq_{num_eq}") * expr ** 2, f"eq_{num_eq}") num_eq += 1 if A_ub is not None: - for i in range(len(A_ub)): + for i, _ in enumerate(A_ub): expr = sum( - A_ub[i][j] * vars[j] for j in range(self.nvars) if A_ub[i][j] != 0 + A_ub[i][j] * model_vars[j] for j in range(self.nvars) if A_ub[i][j] != 0 ) expr -= b_ub[i] slack = LogEncInteger( f"eq_{num_eq}_slack", - (0, LinearProg._get_slack_ub(vars, A_ub[i], b_ub[i])), + (0, LinearProg._get_slack_ub(model_vars, A_ub[i], b_ub[i])), ) H += Constraint( @@ -88,9 +88,9 @@ def _to_bqm_qubo_ising(self, pdict=None): self.num_eq = num_eq pyqubo_model = H.compile() - if pdict == None: + if pdict is None: pdict = {f"eq_{i}": 2 for i in range(self.num_eq)} - elif type(pdict) == int or type(pdict) == float: + elif isinstance(pdict, int) or isinstance(pdict, float): pdict = {f"eq_{i}": pdict for i in range(self.num_eq)} pdict["obj"] = 1 self.qubo = pyqubo_model.to_qubo(feed_dict=pdict) @@ -137,7 +137,7 @@ def _get_slack_ub(vars: list, coefs: list, offset: int) -> int: """ ub = 0 for var, coef in zip(vars, coefs): - if type(var) == LogEncInteger: + if isinstance(var, LogEncInteger): ub += coef * (var.value_range[1] if coef < 0 else var.value_range[0]) else: ub += coef * (1 if coef < 0 else 0) @@ -158,23 +158,23 @@ def _to_cqm(self): ) ind = 0 - vars = [] - for (lb, ub) in self.bounds: + model_vars = [] + for (lb, ub) in bounds: if lb == 0 and ub == 1: - vars.append(dimod.Binary(f"x_{ind}")) + model_vars.append(dimod.Binary(f"x_{ind}")) else: - vars.append(dimod.Integer(f"x_{ind}", lower_bound=lb, upper_bound=ub)) + model_vars.append(dimod.Integer(f"x_{ind}", lower_bound=lb, upper_bound=ub)) ind += 1 cqm = dimod.ConstrainedQuadraticModel() - dimod_obj = sum(var * coef for var, coef in zip(vars, c) if coef != 0) + dimod_obj = sum(var * coef for var, coef in zip(model_vars, c) if coef != 0) cqm.set_objective(dimod_obj) num_eq = 0 if A_eq is not None: for i in range(len(A_eq)): expr = sum( - A_eq[i][j] * vars[j] for j in range(self.nvars) if A_eq[i][j] != 0 + A_eq[i][j] * model_vars[j] for j in range(self.nvars) if A_eq[i][j] != 0 ) new_c = expr == b_eq[i] cqm.add_constraint(new_c, label=f"eq_{num_eq}") @@ -182,7 +182,7 @@ def _to_cqm(self): if A_ub is not None: for i in range(len(A_ub)): expr = sum( - A_ub[i][j] * vars[j] for j in range(self.nvars) if A_ub[i][j] != 0 + A_ub[i][j] * model_vars[j] for j in range(self.nvars) if A_ub[i][j] != 0 ) new_c = expr <= b_ub[i] cqm.add_constraint(new_c, label=f"eq_{num_eq}") @@ -191,9 +191,8 @@ def _to_cqm(self): self.cqm = cqm def _count_qubits(self): - vars = self.bqm.variables - return len(vars) - + model_vars = self.bqm.variables + return len(model_vars) def _count_quadratic_couplings(self): """ @@ -205,7 +204,6 @@ def _count_quadratic_couplings(self): count = count + 1 return count - def _count_linear_fields(self): """ return number of local fields hs @@ -215,6 +213,5 @@ def _count_linear_fields(self): if h != 0: count = count + 1 return count - - + diff --git a/AGV_quantum/__init__.py b/AGV_quantum/__init__.py index 9d6c71f..639e7f2 100644 --- a/AGV_quantum/__init__.py +++ b/AGV_quantum/__init__.py @@ -1,13 +1,8 @@ from .LinearProg import LinearProg - from .process_results import ( - get_results, load_results, store_result, print_results -) - -from .quadratic_solver_CPLEX import ( - load_linear_prog_object, process_result, quadratic_solve_qubo + get_results, load_results, store_result, get_objective, analyze_constraints, process_result ) @@ -17,7 +12,6 @@ create_t_iterator, create_y_iterator, create_z_iterator ) - from .linear_solver import ( LinearAGV, print_ILP_size ) @@ -28,6 +22,5 @@ sim_anneal, annealing, constrained_solver, hybrid_anneal ) - from .train_diagram import plot_train_diagram diff --git a/AGV_quantum/linear_solver.py b/AGV_quantum/linear_solver.py index 71b7e06..4566595 100644 --- a/AGV_quantum/linear_solver.py +++ b/AGV_quantum/linear_solver.py @@ -1,13 +1,14 @@ -from curses import A_LEFT -import itertools -import numpy as np -from typing import Callable, Dict, List, Tuple +"""implementation of linear solver on ILP""" -from AGV_quantum import see_non_zero_variables, create_graph, create_iterators, create_agv_list +import itertools from typing import Optional +import numpy as np from docplex.mp.model import Model from docplex.mp.solution import SolveSolution +from AGV_quantum import see_non_zero_variables, create_graph, create_iterators, create_agv_list + + class LinearAGV: """ diff --git a/AGV_quantum/process_results.py b/AGV_quantum/process_results.py index 785520f..9fe87a8 100644 --- a/AGV_quantum/process_results.py +++ b/AGV_quantum/process_results.py @@ -1,7 +1,7 @@ +"""analyse results given quantum/hybrid/simulation approach""" + import os import pickle -from typing import Any -import dimod from dimod import SampleSet from AGV_quantum import LinearProg @@ -106,26 +106,18 @@ def analyze_constraints(lp, sample): result[f"eq_{num_eq}"] = expr <= lp.b_ub[i] num_eq += 1 - return result, sum(x == False for x in result.values()) - - -def print_results(dict_list): - soln = next((l for l in dict_list if l["feasible"]), None) - if soln is not None: - print("obj:", soln["objective"], "x:", list(soln["sample"].values())) - print("First 10 solutions") - for d in dict_list[:10]: - print(d) - else: - print("No feasible solution") - for d in dict_list[:10]: - print( - "Energy:", - d["energy"], - "Objective:", - d["objective"], - "Feasible", - d["feasible"], - "Broken constraints:", - d["feas_constraints"][1], - ) + return result, sum(x is False for x in result.values()) + + +def process_result(sampleset: list): + energy = sampleset[0]["energy"] + objective = sampleset[0]["objective"] + feasible = sampleset[0]["feasible"] + broken_constrains = [] + if not feasible: + for eq, feas in sampleset[0]["feas_constraints"][0].items(): + if not feas: + broken_constrains.append(eq) + num_broken = len(broken_constrains) + + return {"energy": energy, "objective": objective, "feasible": feasible, "num_broken": num_broken} \ No newline at end of file diff --git a/AGV_quantum/quadratic_solver.py b/AGV_quantum/quadratic_solver.py index 403625e..787c1b2 100644 --- a/AGV_quantum/quadratic_solver.py +++ b/AGV_quantum/quadratic_solver.py @@ -1,3 +1,4 @@ +""" implementation of quadratic (or D-Wace cqm) solver for quantum, hybrid and simulations """ import dimod from cpp_pyqubo import Binary, Constraint, Placeholder from pyqubo import LogEncInteger @@ -7,7 +8,7 @@ class QuadraticAGV: - + """class containing quadratic model of the problem i.e. QUBO, bqm or Ising""" def __init__(self, lin_agv: LinearAGV): self.lin_agv = lin_agv self.c = lin_agv.c @@ -41,35 +42,35 @@ def to_bqm_qubo_ising(self, pdict=None): ) H = 0 ind = 0 - vars = [] + model_vars = [] for (lb, ub) in bounds: if lb == 0 and ub == 1: - vars.append(Binary(f"x_{ind}")) + model_vars.append(Binary(f"x_{ind}")) else: - vars.append(LogEncInteger(f"x_{ind}", (lb, ub))) + model_vars.append(LogEncInteger(f"x_{ind}", (lb, ub))) ind += 1 - pyqubo_obj = sum(var * coef for var, coef in zip(vars, c) if coef != 0) + pyqubo_obj = sum(var * coef for var, coef in zip(model_vars, c) if coef != 0) H += Placeholder("obj") * pyqubo_obj num_eq = 0 if A_eq is not None: - for i in range(len(A_eq)): + for i, _ in enumerate(A_eq): expr = sum( - A_eq[i][j] * vars[j] for j in range(self.nvars) if A_eq[i][j] != 0 + A_eq[i][j] * model_vars[j] for j in range(self.nvars) if A_eq[i][j] != 0 ) expr -= b_eq[i] H += Constraint(Placeholder(f"eq_{num_eq}") * expr ** 2, f"eq_{num_eq}") num_eq += 1 if A_ub is not None: - for i in range(len(A_ub)): + for i, _ in enumerate(A_ub): expr = sum( - A_ub[i][j] * vars[j] for j in range(self.nvars) if A_ub[i][j] != 0 + A_ub[i][j] * model_vars[j] for j in range(self.nvars) if A_ub[i][j] != 0 ) expr -= b_ub[i] slack = LogEncInteger( f"eq_{num_eq}_slack", - (0, self._get_slack_ub(vars, A_ub[i], b_ub[i])), + (0, self._get_slack_ub(model_vars, A_ub[i], b_ub[i])), ) H += Constraint( @@ -81,9 +82,9 @@ def to_bqm_qubo_ising(self, pdict=None): self.num_eq = num_eq pyqubo_model = H.compile() - if pdict == None: + if pdict is None: pdict = {f"eq_{i}": 2 for i in range(self.num_eq)} - elif type(pdict) == int or type(pdict) == float: + elif isinstance(pdict, int) or isinstance(pdict, float): pdict = {f"eq_{i}": pdict for i in range(self.num_eq)} pdict["obj"] = 1 self.qubo = pyqubo_model.to_qubo(feed_dict=pdict) @@ -116,10 +117,10 @@ def interpreter(sampleset: dimod.SampleSet): self.interpreter = lambda ss: interpreter(ss) @staticmethod - def _get_slack_ub(vars: list, coefs: list, offset: int) -> int: + def _get_slack_ub(model_vars: list, coefs: list, offset: int) -> int: """Returns upper bound for slack variables - :param vars: List of variables (can be integer or binary) + :param model_vars: List of variables (can be integer or binary) :type vars: list :param coefs: List of coefficients for the inequality :type coefs: list @@ -129,7 +130,7 @@ def _get_slack_ub(vars: list, coefs: list, offset: int) -> int: :rtype: int """ ub = 0 - for var, coef in zip(vars, coefs): + for var, coef in zip(model_vars, coefs): if type(var) == LogEncInteger: ub += coef * (var.value_range[1] if coef < 0 else var.value_range[0]) else: @@ -152,7 +153,7 @@ def to_cqm(self): ind = 0 vars = [] - for (lb, ub) in self.bounds: + for (lb, ub) in bounds: if lb == 0 and ub == 1: vars.append(dimod.Binary(f"x_{ind}")) else: @@ -184,8 +185,8 @@ def to_cqm(self): self.cqm = cqm def _count_qubits(self): - vars = self.bqm.variables - return len(vars) + model_vars = self.bqm.variables + return len(model_vars) def _count_quadratic_couplings(self): diff --git a/AGV_quantum/quadratic_solver_CPLEX.py b/AGV_quantum/quadratic_solver_CPLEX.py deleted file mode 100644 index fa7384b..0000000 --- a/AGV_quantum/quadratic_solver_CPLEX.py +++ /dev/null @@ -1,232 +0,0 @@ -import pickle -from docplex.mp.model import Model -from docplex.mp.solution import SolveSolution -import dimod -from docplex.mp.progress import ProgressListener, ProgressClock, TextProgressListener -from typing import Union -from math import inf - -from AGV_quantum import get_results -from AGV_quantum import LinearProg - - -class AutomaticAborter(ProgressListener): - """ a simple implementation of an automatic search stopper. - WIP - """ - - def __init__(self, lp_instance: LinearProg, out_path:str, name:str, given_time:float = 3600., gap_fmt=None, obj_fmt=None): - super(AutomaticAborter, self).__init__(ProgressClock.All) - self.last_obj = None - self.given_time = given_time - self.lp = lp_instance - self.path = out_path - self.name = name - self.sol_found = False - self._gap_fmt = gap_fmt or "{:.2%}" - self._obj_fmt = obj_fmt or "{:.4f}" - self._count = 0 - - def notify_start(self): - super(AutomaticAborter, self).notify_start() - self.last_obj = inf - - def process_result(self, sampleset: list): - energy = sampleset[0]["energy"] - objective = sampleset[0]["objective"] - feasible = sampleset[0]["feasible"] - broken_constrains = [] - if not feasible: - for eq, feas in sampleset[0]["feas_constraints"][0].items(): - if not feas: - broken_constrains.append(eq) - num_broken = len(broken_constrains) - - return {"energy": energy, "objective": objective, "feasible": feasible, "broken_constrains": broken_constrains, - "num_broken": num_broken} - - def check_solution(self, sol: SolveSolution): - - solved_list = [v.name for v in sol.iter_variables()] - sample = {v: 1 if v in solved_list else 0 for v in self.lp.bqm.variables} - sampleset = dimod.SampleSet.from_samples(dimod.as_samples(sample), 'BINARY', sol.objective_value) - sampleset = self.lp.interpreter(sampleset) - results = get_results(sampleset, prob=self.lp) - results = process_result(results) - - return results["feasible"], results - - def save_results(self, results: dict): - - with open(self.path, "a") as f: - f.write(f"{self.name}: ") - f.write(str(results)) - - def is_improving(self, new_obj, eps=1e-2): - last_obj = self.last_obj - return abs(new_obj - last_obj) >= eps - - def notify_solution(self, sol): - feasible, results = self.check_solution(sol) - print("test") - if feasible: - self.sol_found = True - self.save_results(results) - print("found feasible solution!") - - def notify_progress(self, pdata): - super(AutomaticAborter, self).notify_progress(pdata) - - if pdata.has_incumbent and self.is_improving(pdata.current_objective): - - self.last_obj = pdata.current_objective - - self._count += 1 - pdata_has_incumbent = pdata.has_incumbent - incumbent_symbol = '+' if pdata_has_incumbent else ' ' - # if pdata_has_incumbent: - # self._incumbent_count += 1 - current_obj = pdata.current_objective - if pdata_has_incumbent: - objs = self._obj_fmt.format(current_obj) - else: - objs = "N/A" # pragma: no cover - best_bound = pdata.best_bound - nb_nodes = pdata.current_nb_nodes - remaining_nodes = pdata.remaining_nb_nodes - if pdata_has_incumbent: - gap = self._gap_fmt.format(pdata.mip_gap) - else: - gap = "N/A" # pragma: no cover - raw_time = pdata.time - rounded_time = round(raw_time, 1) - - print("{0:>3}{7}: Node={4} Left={5} Best Integer={1}, Best Bound={2:.4f}, gap={3}, ItCnt={8} [{6}s]" - .format(self._count, objs, best_bound, gap, nb_nodes, remaining_nodes, rounded_time, - incumbent_symbol, pdata.current_nb_iterations)) - - if pdata.time > self.given_time: - if self.sol_found: - self.abort() - print("Feasible solution found, aborting search") - - - - -def add_zero_h_qubo(lp: LinearProg) -> LinearProg: - for var in lp.bqm.variables: - if (var, var) not in lp.qubo[0]: - lp.qubo[0][(var, var)] = 0 - return lp - - -def add_zero_h_ising(lp: LinearProg) -> LinearProg: - for var in lp.bqm.variables: - if var not in lp.ising[0]: - lp.ising[0][var] = 0 - return lp - - -def load_linear_prog_object(lp_location: str) -> LinearProg: - with open(lp_location, "rb") as f: - lp = pickle.load(f) - p = 2.75 - lp._to_bqm_qubo_ising(p) - return lp - - -def count_vertices(qubo: dict) -> int: - s = 0 - for key in qubo.keys(): - if key[0] == key[1]: - s += 1 - return s - - -def count_edges(qubo: dict) -> int: - s = 0 - for key in qubo.keys(): - if key[0] != key[1]: - s += 1 - return s - - -def quadratic_solve_qubo(lp_location: str, num_threads: int = None) -> (SolveSolution, LinearProg): - lp = load_linear_prog_object(lp_location) - p = 2.75 - lp._to_bqm_qubo_ising(p) - lp = add_zero_h_qubo(lp) - qubo = lp.qubo[0] - - m = Model(name='qubo') - if num_threads: - m.context.cplex_parameters.threads = num_threads - variables = m.binary_var_dict(lp.bqm.variables, name="", key_format="%s") - obj_fnc = sum(variables[k1] * variables[k2] * qubo[(k1, k2)] for k1, k2 in qubo.keys()) - m.set_objective("min", obj_fnc) - #m.print_information() - m.add_progress_listener(TextProgressListener(clock='Objective')) - sol = m.solve() - m.print_solution() - return sol, lp - - -def quadratic_model(lp_location: str) -> Model: - """ - :param lp_location: path to lp file - :return: docplex model for given qubo - """ - lp = load_linear_prog_object(lp_location) - lp = add_zero_h_qubo(lp) - qubo = lp.qubo[0] - - model = Model(name='qubo') - variables = model.binary_var_dict(lp.bqm.variables, name="", key_format="%s") - obj_fnc = sum(variables[k1] * variables[k2] * qubo[(k1, k2)] for k1, k2 in qubo.keys()) - model.set_objective("min", obj_fnc) - - return model - -def process_result(sampleset: list): - energy = sampleset[0]["energy"] - objective = sampleset[0]["objective"] - feasible = sampleset[0]["feasible"] - broken_constrains = [] - if not feasible: - for eq, feas in sampleset[0]["feas_constraints"][0].items(): - if not feas: - broken_constrains.append(eq) - num_broken = len(broken_constrains) - - return {"energy": energy, "objective": objective, "feasible": feasible, "num_broken": num_broken} - - -def check_solution(sol: Union[SolveSolution, dict], lp: LinearProg): - if isinstance(sol, dict): - raise NotImplementedError - else: - solved_list = [v.name for v in sol.iter_variables()] - sample = {v: 1 if v in solved_list else 0 for v in lp.bqm.variables} - sampleset = dimod.SampleSet.from_samples(dimod.as_samples(sample), 'BINARY', sol.objective_value) - sampleset = lp.interpreter(sampleset) - results = get_results(sampleset, prob=lp) - results = process_result(results) - - return results["feasible"], results - - -def save_results(results: dict, name:str, output_path: str): - - with open(output_path, "a") as f: - f.write(f"{name}: ") - f.write(str(results)) - - -if __name__ == "__main__": - for name in ["tiny", "smallest", "small", "medium_small"]: - sol, lp = quadratic_solve_qubo(f"lp_{name}.pkl") - sol.export(f"sol_{name}.json") - feasible, results = check_solution(sol, lp) - save_results(results, f"{name}", "results.txt") - - diff --git a/AGV_quantum/qubo_solver.py b/AGV_quantum/qubo_solver.py index 7ea9aab..66e0a0a 100644 --- a/AGV_quantum/qubo_solver.py +++ b/AGV_quantum/qubo_solver.py @@ -1,4 +1,5 @@ -from typing import Callable, Dict, List, Tuple +"""D-Wave implementations""" +from typing import Optional import os import neal @@ -9,11 +10,10 @@ LeapHybridSampler, LeapHybridCQMSampler, ) -#from scipy.optimize import linprog from AGV_quantum import LinearProg from AGV_quantum import get_results, load_results, store_result -from typing import Optional + def sim_anneal( bqm: dimod.BinaryQuadraticModel, @@ -140,7 +140,7 @@ def get_parameters(real_anneal_var_dict): :return: Number of reads, annealing_time and chain strength :rtype: Tuple[int, int, int] """ - if real_anneal_var_dict == None: + if real_anneal_var_dict is None: num_reads = 1000 annealing_time = 250 chain_strength = 4 @@ -243,4 +243,3 @@ def annealing( sampleset = lp.interpreter(sampleset) return get_results(sampleset, prob=lp) - diff --git a/AGV_quantum/qubo_to_matrix.py b/AGV_quantum/qubo_to_matrix.py deleted file mode 100644 index e095fa6..0000000 --- a/AGV_quantum/qubo_to_matrix.py +++ /dev/null @@ -1,13 +0,0 @@ -from AGV_quantum import load_linear_prog_object -from AGV_quantum import qubo_to_matrix -import scipy.sparse as sparse - -if __name__ == "__main__": - size = "smallest" - p = 2.75 - lp = load_linear_prog_object(f"../lp_files/lp_{size}.pkl") - lp._to_bqm_qubo_ising(p) - qubo = lp.qubo[0] - matrix = qubo_to_matrix(qubo, lp) - matrix = sparse.coo_matrix(matrix) - sparse.save_npz(f"../qubo/{size}_qubo_coo.npz", matrix) diff --git a/AGV_quantum/train_diagram.py b/AGV_quantum/train_diagram.py index 9d10dab..d1f0fa0 100644 --- a/AGV_quantum/train_diagram.py +++ b/AGV_quantum/train_diagram.py @@ -1,17 +1,16 @@ +"""plots train diagram for given solutions""" import matplotlib.pyplot as plt import matplotlib as mpl - mpl.rc("text", usetex=True) mpl.rc("font", family="serif") mpl.rc("font", size=10) def zones_location(track_len, n_zones, s_ofset): - + """determines zones locations and borders in space for the plot""" marks = {f"s{k}":0 for k in range(n_zones)} zone_borders = [] - prev_s = -100 x = 0 for s in marks: @@ -19,12 +18,12 @@ def zones_location(track_len, n_zones, s_ofset): x = x + s_ofset + track_len[(f"{prev_s}", f"{s}")] marks[s] = x + s_ofset/2 zone_borders.append(x) - zone_borders.append(x+1) + zone_borders.append(x+s_ofset) prev_s = s return marks, zone_borders - def AGVS_coordinates(sol, agv_routes, marks, s_ofset): + """ determines coordinates of AGVs in space and time to get their paths on the plot """ times = {} spaces = {} for agv in agv_routes: @@ -38,34 +37,26 @@ def AGVS_coordinates(sol, agv_routes, marks, s_ofset): return times, spaces - def plot_train_diagram(sol, agv_routes, track_len, n_zones): + """plots and saves train diagram""" plt.figure(figsize=(4, 3)) - - s_ofset = 1 + s_ofset = 1.75 # the size of the station marks, zone_borders = zones_location(track_len, n_zones, s_ofset) - times, spaces = AGVS_coordinates(sol, agv_routes, marks, s_ofset) - - - colors = {0: "black", 1: "red", 2: "green", 3: "blue", 4: "orange", 5: "brown", 6: "cyan"} for agv in agv_routes: - plt.plot(times[agv], spaces[agv], "o-", label=f"AGV {agv} ", color=colors[agv], linewidth=0.85, markersize=2) + plt.plot(times[agv], spaces[agv], "o-", label=f"$AGV_{agv}$ ", color=colors[agv], linewidth=0.85, markersize=2) plt.legend(loc='upper center', bbox_to_anchor=(0.5, 1.45), ncol = 3) for el in zone_borders: plt.axhline(y = el, color="gray", linewidth=0.5, linestyle=":") - locs = [marks[k] for k in marks] + + our_marks = [f"$s_{i}$" for i, _ in enumerate(marks) ] - plt.yticks(locs, [k for k in marks]) - - + plt.yticks(locs, our_marks) plt.xlabel("time") plt.ylabel("zones") - plt.subplots_adjust(bottom=0.19, top = 0.75) - - - plt.savefig("train_diagram.pdf") \ No newline at end of file + plt.subplots_adjust(bottom=0.19, top = 0.75) + plt.savefig("train_diagram.pdf") diff --git a/AGV_quantum/utils.py b/AGV_quantum/utils.py index 2546ea1..cbbeaae 100644 --- a/AGV_quantum/utils.py +++ b/AGV_quantum/utils.py @@ -24,6 +24,7 @@ def number_gen(): def create_stations_list(tracks): + """return list of zones""" stations = [] for track in tracks: for station in track: @@ -32,7 +33,8 @@ def create_stations_list(tracks): return list(set(stations)) -def agv_routes_as_edges(agv_routes): # TODO write test +def agv_routes_as_edges(agv_routes): + """ """ return_dict = {} for j in agv_routes.keys(): if len(agv_routes[j]) > 1: diff --git a/README.md b/README.md index 329f53b..9825363 100644 --- a/README.md +++ b/README.md @@ -12,13 +12,17 @@ In directory ```examples``` there are AGVs scheduling problems. In order of incr - example_large.py - example_largest.py -There are optional boolean parameters (```1``` yes, ```0``` no): ```--solve_linear``` - solve on CPLEX , ```--train_diagram``` - plot "train diagram" for given problem ```--solve_quadratic``` - solve on hybrid quantum classical (the particular solver and penalty parameter can be set in the script). - ## Usage To run this project in terminal use path/to/project> python -m examples.file +There are optional boolean parameters (```1``` yes, ```0``` no): ```--solve_linear``` - solve on CPLEX , ```--train_diagram``` - plot "train diagram" for given problem ```--solve_quadratic``` - solve on hybrid quantum classical (the particular solver and penalty parameter can be set in the script). + +Example: + +```python -m examples.example_smallest --solve_linear 1 --train_diagram 1 ``` + To run tests use path/to/project> python3 -m unittest ## Computational results diff --git a/examples/example_large.py b/examples/example_large.py index a16cb5b..db9b2f8 100644 --- a/examples/example_large.py +++ b/examples/example_large.py @@ -3,10 +3,10 @@ import time import os -from src import create_stations_list, create_agv_list, create_graph, create_same_way_dict, agv_routes_as_edges -from src import print_ILP_size, LinearAGV -from src import QuadraticAGV -from src import annealing, constrained_solver, hybrid_anneal +from AGV_quantum import create_stations_list, create_agv_list, create_graph, create_same_way_dict, agv_routes_as_edges +from AGV_quantum import print_ILP_size, LinearAGV +from AGV_quantum import QuadraticAGV +from AGV_quantum import annealing, constrained_solver, hybrid_anneal diff --git a/examples/example_largest.py b/examples/example_largest.py index e853fcc..6fe3971 100644 --- a/examples/example_largest.py +++ b/examples/example_largest.py @@ -3,12 +3,11 @@ import time import os -from src import create_stations_list, create_agv_list, create_graph, create_same_way_dict, agv_routes_as_edges -from src import print_ILP_size, LinearAGV -from src import QuadraticAGV -from src import annealing, constrained_solver, hybrid_anneal +from AGV_quantum import create_stations_list, create_agv_list, create_graph, create_same_way_dict, agv_routes_as_edges +from AGV_quantum import print_ILP_size, LinearAGV +from AGV_quantum import QuadraticAGV +from AGV_quantum import annealing, constrained_solver, hybrid_anneal -from src.LinearProg import LinearProg cwd = os.getcwd() diff --git a/examples/example_medium.py b/examples/example_medium.py index 96513f8..443c851 100644 --- a/examples/example_medium.py +++ b/examples/example_medium.py @@ -3,11 +3,11 @@ import time import os -from src import create_stations_list, create_agv_list, create_graph, create_same_way_dict, agv_routes_as_edges -from src import plot_train_diagram -from src import annealing, constrained_solver, hybrid_anneal -from src import print_ILP_size, LinearAGV -from src import QuadraticAGV +from AGV_quantum import create_stations_list, create_agv_list, create_graph, create_same_way_dict, agv_routes_as_edges +from AGV_quantum import plot_train_diagram +from AGV_quantum import annealing, constrained_solver, hybrid_anneal +from AGV_quantum import print_ILP_size, LinearAGV +from AGV_quantum import QuadraticAGV diff --git a/examples/example_small.py b/examples/example_small.py index 8721763..7688b33 100644 --- a/examples/example_small.py +++ b/examples/example_small.py @@ -3,11 +3,11 @@ import time import os -from src import create_stations_list, create_agv_list, create_graph, create_same_way_dict, agv_routes_as_edges -from src import plot_train_diagram -from src import print_ILP_size, LinearAGV -from src import QuadraticAGV -from src import annealing, constrained_solver, hybrid_anneal +from AGV_quantum import create_stations_list, create_agv_list, create_graph, create_same_way_dict, agv_routes_as_edges +from AGV_quantum import plot_train_diagram +from AGV_quantum import print_ILP_size, LinearAGV +from AGV_quantum import QuadraticAGV +from AGV_quantum import annealing, constrained_solver, hybrid_anneal cwd = os.getcwd() diff --git a/examples/example_smallest.py b/examples/example_smallest.py index fe85aee..132f03f 100644 --- a/examples/example_smallest.py +++ b/examples/example_smallest.py @@ -4,11 +4,11 @@ import time import os -from src import create_stations_list, create_agv_list, create_graph, create_same_way_dict, agv_routes_as_edges -from src import plot_train_diagram -from src import print_ILP_size, LinearAGV -from src import QuadraticAGV -from src import annealing, constrained_solver, hybrid_anneal +from AGV_quantum import create_stations_list, create_agv_list, create_graph, create_same_way_dict, agv_routes_as_edges +from AGV_quantum import plot_train_diagram +from AGV_quantum import print_ILP_size, LinearAGV +from AGV_quantum import QuadraticAGV +from AGV_quantum import annealing, constrained_solver, hybrid_anneal M = 20 diff --git a/tests/test_converters.py b/tests/test_converters.py index 15da146..f968030 100644 --- a/tests/test_converters.py +++ b/tests/test_converters.py @@ -1,6 +1,6 @@ import unittest from AGV_quantum import LinearProg -from AGV_quantum import sim_anneal, annealing, process_results +from AGV_quantum import sim_anneal, annealing, get_objective, analyze_constraints import dimod from scipy.optimize import linprog import numpy as np @@ -38,16 +38,20 @@ def test_feasibility_check(self): bqm = self.lp.bqm sampleset = sim_anneal(bqm, beta_range=(0.01, 10), num_sweeps=1000, num_reads=1000) sampleset = self.lp.interpreter(sampleset) - prob=self.lp dict_list = [] + + assert analyze_constraints(prob, {'x_0': 5, 'x_1': 4}) == ({'eq_0': True, 'eq_1': True, 'eq_2': True, 'eq_3': True}, 0) + assert get_objective(prob, {'x_0': 5, 'x_1': 4}) == -13 + + for data in sampleset.data(): rdict = {} sample = data.sample rdict["energy"] = data.energy - rdict["objective"] = round(process_results.get_objective(prob, sample), 2) - rdict["feasible"] = all(process_results.analyze_constraints(prob, sample)[0].values()) - rdict["feas_constraints"] = process_results.analyze_constraints(prob, sample) + rdict["objective"] = round(get_objective(prob, sample), 2) + rdict["feasible"] = all(analyze_constraints(prob, sample)[0].values()) + rdict["feas_constraints"] = analyze_constraints(prob, sample) dict_list.append(rdict) ret = sorted(dict_list, key=lambda d: d["energy"])[0] if make_probabilistic_test: @@ -55,6 +59,7 @@ def test_feasibility_check(self): assert ret["feas_constraints"][1] == 0 assert ret["objective"] < -9. + def test_bqm_soln_sim(self): dict_list = annealing(self.lp, "sim", "test_1", load=False, store=False) soln = (next((l for l in dict_list if l["feasible"]), None)) diff --git a/tests/test_example.py b/tests/test_example.py index 1b84bd4..d638b8c 100644 --- a/tests/test_example.py +++ b/tests/test_example.py @@ -3,7 +3,7 @@ import unittest from AGV_quantum import create_stations_list, create_agv_list -from AGV_quantum import print_ILP_size, LinearAGV +from AGV_quantum import print_ILP_size, LinearAGV, QuadraticAGV @@ -63,6 +63,13 @@ def test_line_fragemnt(self): assert sol["t_in_0_s1"] == 6. assert sol["t_out_0_s1"] == 7. + # quadratic testing + model_q = QuadraticAGV(AGV) + p = 5. + model_q.to_bqm_qubo_ising(p) + assert model_q.qubo[0][('x_0[0]', 'x_0[0]')] == -10.0 + model_q.to_cqm() + print(model_q.cqm.constraints['eq_0']) if __name__ == "__main__": diff --git a/tests/test_utility_functions.py b/tests/test_utility_functions.py index 10a509e..ee54894 100644 --- a/tests/test_utility_functions.py +++ b/tests/test_utility_functions.py @@ -1,7 +1,7 @@ import unittest import networkx as nx from AGV_quantum import create_stations_list, create_agv_list, create_graph, create_t_iterator, create_y_iterator, create_z_iterator -from AGV_quantum import create_iterators, create_v_in_out +from AGV_quantum import create_iterators, create_v_in_out, agv_routes_as_edges class SingleStation(unittest.TestCase): @@ -62,6 +62,9 @@ def setUpClass(cls): cls.stations = create_stations_list(cls.tracks) cls.tau_operation = {(agv, station): 2 for agv in cls.J for station in cls.stations} + def test_AGVs_routes(self): + assert agv_routes_as_edges(self.agv_routes) == {0: [('s0', 's1')], 1: [('s0', 's2'), ('s2', 's3')]} + def test_create_lists_multi(self): stations = create_stations_list(self.tracks) J = create_agv_list(self.agv_routes)