diff --git a/solvers/cplex/mip/mip.py b/solvers/cplex/mip/mip.py new file mode 100644 index 00000000..acf30d15 --- /dev/null +++ b/solvers/cplex/mip/mip.py @@ -0,0 +1,34 @@ +import argparse +import cplex +import sys + +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) + sys.exit(1) + + +try: + cpx.solve() +except cplex.exceptions.CplexError as exc: + print("Error during solve:", exc) + sys.exit(1) + +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) diff --git a/solvers/cplex/mip/requirements.txt b/solvers/cplex/mip/requirements.txt new file mode 100644 index 00000000..4cfbd949 --- /dev/null +++ b/solvers/cplex/mip/requirements.txt @@ -0,0 +1 @@ +cplex 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/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_gms_reformulation.gms b/solvers/gams/mip/mip_to_gms_reformulation.gms new file mode 100644 index 00000000..6019d1a7 --- /dev/null +++ b/solvers/gams/mip/mip_to_gms_reformulation.gms @@ -0,0 +1,80 @@ +$setEnv GAMSINP "%INPUT%" +$setEnv GAMSPEN "%PENALTY%" + +$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 or .mps 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 +) + +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' + +# 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}' + +# 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 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 diff --git a/src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java b/src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java new file mode 100644 index 00000000..71bc7172 --- /dev/null +++ b/src/main/java/edu/kit/provideq/toolbox/mip/MipConfiguration.java @@ -0,0 +1,60 @@ +package edu.kit.provideq.toolbox.mip; + +import edu.kit.provideq.toolbox.ResourceProvider; +import edu.kit.provideq.toolbox.exception.MissingExampleException; +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.CplexMipSolver; +import edu.kit.provideq.toolbox.mip.solvers.QuboMipSolver; +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 "Mixed Integer" problem. + */ +@Configuration +public class MipConfiguration { + /** + * MIP (Mixed integer Problem) + * A mixed integer optimization problem. + * For a given integer term, + * find the minimal / maximal variable assignment of the term. + */ + public static final ProblemType MIP = new ProblemType<>( + "mip", + String.class, + String.class + ); + + @Bean + ProblemManager getMipManager( + CplexMipSolver cplexMip, + QuboMipSolver quboMipSolver, + ResourceProvider resourceProvider + ) { + return new ProblemManager<>( + MIP, + Set.of(cplexMip, quboMipSolver), + loadExampleProblems(resourceProvider) + ); + } + + private Set> loadExampleProblems( + ResourceProvider resourceProvider) { + try { + var problemInputStream = Objects.requireNonNull( + getClass().getResourceAsStream("simple.lp"), + "mixed-integer-problem example is unavailable!" + ); + var problem = new Problem<>(MIP); + problem.setInput(resourceProvider.readStream(problemInputStream)); + return Set.of(problem); + } catch (IOException e) { + throw new MissingExampleException(MIP, e); + } + } +} diff --git a/src/main/java/edu/kit/provideq/toolbox/mip/solvers/CplexMipSolver.java b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/CplexMipSolver.java new file mode 100644 index 00000000..eca7a3f9 --- /dev/null +++ b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/CplexMipSolver.java @@ -0,0 +1,67 @@ +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.mip.MipConfiguration; +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; + +/** + * {@link MipConfiguration#MIP} solver using IBM cplex. + */ +@Component +public class CplexMipSolver extends MipSolver { + private final String scriptPath; + private final String venv; + + private final ApplicationContext context; + + @Autowired + public CplexMipSolver( + @Value("${path.cplex.mip}") String scriptPath, + @Value("${venv.cplex.mip}") String venv, + ApplicationContext context) { + this.scriptPath = scriptPath; + this.venv = venv; + this.context = context; + } + + @Override + public String getName() { + return "Cplex"; + } + + @Override + public String getDescription() { + return "This solver uses IBM Cplex Solver to solve MIP problems. " + + "It is called via Python wrapper without any additional parameters"; + } + + @Override + public Mono> solve( + String input, + SubRoutineResolver subRoutineResolver, + SolvingProperties properties + ) { + var solution = new Solution<>(this); + + var processResult = context + .getBean(PythonProcessRunner.class, scriptPath, venv) + .withArguments( + ProcessRunner.INPUT_FILE_PATH, + ProcessRunner.OUTPUT_FILE_PATH + ) + .writeInputFile(input) + .readOutputFile("output.txt") + .run(getProblemType(), solution.getId()); + + // Return if process failed + return Mono.just(processResult.applyTo(solution)); + } +} 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..aedd52f7 --- /dev/null +++ b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/MipSolver.java @@ -0,0 +1,15 @@ +package edu.kit.provideq.toolbox.mip.solvers; + +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() { + return MipConfiguration.MIP; + } +} 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..e48c2c87 --- /dev/null +++ b/src/main/java/edu/kit/provideq/toolbox/mip/solvers/QuboMipSolver.java @@ -0,0 +1,79 @@ +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.mip.MipConfiguration; +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; +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 mipToGamsReformulatorPath; + private final String solverScriptPath; + private final ApplicationContext context; + private static final int PENALTY = 100; + + @Autowired + public QuboMipSolver( + @Value("${path.gams.mip-to-gams-reformulation}") String mipToGamsReformulatorPath, + @Value("${path.gams.gams-to-qubo-reformulation}") String solverScriptPath, + ApplicationContext context) { + this.mipToGamsReformulatorPath = mipToGamsReformulatorPath; + this.solverScriptPath = solverScriptPath; + this.context = context; + } + + @Override + public String getName() { + return "Gams Solver for MIP"; + } + + @Override + public String getDescription() { + return "This solver uses LP / MIP reformulation to QUBO which " + + "is then subsequently solved by gams"; + } + + @Override + public Mono> solve( + String input, + SubRoutineResolver subRoutineResolver, + SolvingProperties properties + ) { + var solution = new Solution<>(this); + + // Run GAMS via console + context + .getBean(GamsProcessRunner.class, mipToGamsReformulatorPath) + .withArguments( + "--INPUT=" + ProcessRunner.INPUT_FILE_PATH, + "--PENALTY=" + PENALTY + ) + .writeInputFile(input, "problem.lp") + .readOutputFile("problem.lp") + .run(getProblemType(), solution.getId()); + + String relativePath = "jobs/mip/" + solution.getId() + "/problem.gms"; + + ProcessResult processResult = context + .getBean(GamsProcessRunner.class, relativePath) + .withArguments( + "--IDIR=" + new File(solverScriptPath).getAbsolutePath(), + "--SOLOUTPUT=" + ProcessRunner.OUTPUT_FILE_PATH + ) + .readOutputFile("output.txt") + .run(getProblemType(), solution.getId()); + return Mono.just(processResult.applyTo(solution)); + } +} 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 86d54fb5..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.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.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 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..533b6950 --- /dev/null +++ b/src/main/resources/edu/kit/provideq/toolbox/mip/simple.lp @@ -0,0 +1,16 @@ +Maximize + obj: 3 x + 4 y + +Subject To + c1: x + 2 y <= 14 + c2: -3 x + y <= 0 + c3: x - y <= 2 + +Bounds + 0 <= x <= 6 + 0 <= y <= 7 + +Generals + x + y +End 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..882d77cf --- /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 testMipSolver(ProblemSolver solver, String input) { + var problem = ApiTestHelper.createProblem(client, solver, input, MIP); + ApiTestHelper.testSolution(problem); + } +}