Skip to content

Commit

Permalink
Compatibility with cvxpy 1.4.1
Browse files Browse the repository at this point in the history
  • Loading branch information
@maxschaller
maxschaller committed Oct 13, 2023
1 parent d2a61ae commit 34826c1
Show file tree
Hide file tree
Showing 5 changed files with 87 additions and 76 deletions.
44 changes: 31 additions & 13 deletions cvxpygen/cpg.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@
from cvxpy.expressions.variable import upper_tri_to_full


def generate_code(problem, code_dir='CPG_code', solver=None, enable_settings=[], unroll=False, prefix='', wrapper=True):
def generate_code(problem, code_dir='CPG_code', solver=None, solver_opts=None, enable_settings=[], unroll=False, prefix='', wrapper=True):
"""
Generate C code for CVXPY problem and (optionally) python wrapper
"""
Expand All @@ -40,10 +40,6 @@ def generate_code(problem, code_dir='CPG_code', solver=None, enable_settings=[],
create_folder_structure(code_dir)

# problem data
solver_opts = {}
# TODO support quadratic objective for SCS.
if solver == 'SCS':
solver_opts['use_quad_obj'] = False
data, solving_chain, inverse_data = problem.get_problem_data(
solver=solver,
gp=False,
Expand Down Expand Up @@ -81,16 +77,16 @@ def generate_code(problem, code_dir='CPG_code', solver=None, enable_settings=[],

constraint_info = get_constraint_info(solver_interface)

adjacency, parameter_canon = process_canonical_parameters(constraint_info, param_prob,
adjacency, parameter_canon, canon_p_ids = process_canonical_parameters(constraint_info, param_prob,
parameter_info, solver_interface,
solver_name, problem)
solver_opts, problem, cvxpy_interface_class)

cvxpygen_directory = os.path.dirname(os.path.realpath(__file__))
solver_code_dir = os.path.join(code_dir, 'c', 'solver_code')
solver_interface.generate_code(code_dir, solver_code_dir, cvxpygen_directory, parameter_canon)

parameter_canon.user_p_name_to_canon_outdated = {
user_p_name: [solver_interface.canon_p_ids[j] for j in np.nonzero(adjacency[:, i])[0]]
user_p_name: [canon_p_ids[j] for j in np.nonzero(adjacency[:, i])[0]]
for i, user_p_name in enumerate(parameter_info.names)}

write_c_code(problem, configuration, variable_info, dual_variable_info, parameter_info,
Expand All @@ -101,11 +97,30 @@ def generate_code(problem, code_dir='CPG_code', solver=None, enable_settings=[],
if wrapper:
compile_python_module(code_dir)

def process_canonical_parameters(constraint_info, param_prob, parameter_info, solver_interface, solver_name, problem):
adjacency = np.zeros(shape=(len(solver_interface.canon_p_ids), parameter_info.num), dtype=bool)

def get_quad_obj(problem, solver_type, solver_opts, solver_class):
if solver_type == 'quadratic':
return True
if solver_opts is None:
use_quad_obj = True
else:
use_quad_obj = solver_opts.get('use_quad_obj', True)
return use_quad_obj and solver_class().supports_quad_obj() and \
problem.objective.expr.has_quadratic_term()


def process_canonical_parameters(constraint_info, param_prob, parameter_info, solver_interface, solver_opts, problem, cvxpy_interface_class):
parameter_canon = ParameterCanon()
parameter_canon.quad_obj = get_quad_obj(problem, solver_interface.solver_type, solver_opts, cvxpy_interface_class)

if not parameter_canon.quad_obj:
canon_p_ids = [p_id for p_id in solver_interface.canon_p_ids if p_id != 'P']
else:
canon_p_ids = solver_interface.canon_p_ids

adjacency = np.zeros(shape=(len(canon_p_ids), parameter_info.num), dtype=bool)
# compute affine mapping for each canonical parameter
for i, (p_id, p_sign) in enumerate(zip(solver_interface.canon_p_ids, solver_interface.canon_p_signs)):
for i, p_id in enumerate(canon_p_ids):

affine_map = solver_interface.get_affine_map(p_id, param_prob, constraint_info)

Expand All @@ -119,13 +134,16 @@ def process_canonical_parameters(constraint_info, param_prob, parameter_info, so

adjacency = update_adjacency_matrix(adjacency, i, parameter_info, affine_map.mapping)

# take sign into account
affine_map.mapping = sparse.csc_matrix(affine_map.mapping.toarray() * affine_map.sign) # be able to use broadcasting

# take sparsity into account
affine_map.mapping = affine_map.mapping[:, parameter_info.sparsity_mask]

# compute default values of canonical parameters
affine_map, parameter_canon = set_default_values(affine_map, p_id, parameter_canon, parameter_info, solver_interface)

parameter_canon.p_id_to_mapping[p_id] = p_sign * affine_map.mapping.tocsr()
parameter_canon.p_id_to_mapping[p_id] = affine_map.mapping.tocsr()
parameter_canon.p_id_to_changes[p_id] = affine_map.mapping[:, :-1].nnz > 0
parameter_canon.p_id_to_size[p_id] = affine_map.mapping.shape[0]

Expand All @@ -136,7 +154,7 @@ def process_canonical_parameters(constraint_info, param_prob, parameter_info, so
parameter_canon.p_id_to_size[p_id] = 0

parameter_canon.is_maximization = type(problem.objective) == Maximize
return adjacency, parameter_canon
return adjacency, parameter_canon, canon_p_ids


def update_to_dense_mapping(affine_map, param_prob):
Expand Down
2 changes: 2 additions & 0 deletions cvxpygen/mappings.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@ class Configuration:
class AffineMap:
mapping_rows: list = field(default_factory=list)
mapping: list = field(default_factory=list)
sign: int = 1
indices: list = field(default_factory=list)
indptr: list = field(default_factory=list)
shape = ()
Expand All @@ -31,6 +32,7 @@ class ParameterCanon:
nonzero_d: bool = True
is_maximization: bool = False
user_p_name_to_canon_outdated: dict[str, list[str]] = field(default_factory=dict)
quad_obj: bool = True


@dataclass
Expand Down
101 changes: 45 additions & 56 deletions cvxpygen/solvers.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
from platform import system

import numpy as np
import scipy as sp

from cvxpygen.utils import replace_in_file, write_struct_prot, write_struct_def, write_vec_prot, write_vec_def
from cvxpygen.mappings import PrimalVariableInfo, DualVariableInfo, ConstraintInfo, AffineMap, \
Expand Down Expand Up @@ -94,7 +95,13 @@ def configure_settings(self) -> None:
def get_affine_map(self, p_id, param_prob, constraint_info: ConstraintInfo) -> AffineMap:
affine_map = AffineMap()

if p_id == 'c':
if p_id == 'P':
if self.indices_obj is None: # problem is an LP
return None
affine_map.mapping = param_prob.reduced_P.reduced_mat
affine_map.indices = self.indices_obj
affine_map.shape = (self.n_var, self.n_var)
elif p_id in ['q', 'c']:
affine_map.mapping = param_prob.c[:-1]
elif p_id == 'd':
affine_map.mapping = param_prob.c[-1]
Expand Down Expand Up @@ -122,8 +129,9 @@ def get_affine_map(self, p_id, param_prob, constraint_info: ConstraintInfo) -> A
elif p_id in ['G', 'h']:
affine_map.indices = self.indices_constr[affine_map.mapping_rows] - self.n_eq

if p_id.isupper():
affine_map.mapping = -param_prob.reduced_A.reduced_mat[affine_map.mapping_rows]
if p_id in ['A', 'G']:
affine_map.mapping = param_prob.reduced_A.reduced_mat[affine_map.mapping_rows]
affine_map.sign = -1

return affine_map

Expand Down Expand Up @@ -160,18 +168,17 @@ def generate_code(self, code_dir, solver_code_dir, cvxpygen_directory,
parameter_canon: ParameterCanon) -> None:
pass

def declare_workspace(self, f, prefix) -> None:
def declare_workspace(self, f, prefix, parameter_canon) -> None:
pass

def define_workspace(self, f, prefix) -> None:
def define_workspace(self, f, prefix, parameter_canon) -> None:
pass


class OSQPInterface(SolverInterface):
solver_name = 'OSQP'
solver_type = 'quadratic'
canon_p_ids = ['P', 'q', 'd', 'A', 'l', 'u']
canon_p_signs = [1, 1, 1, 1, -1, -1]
canon_p_ids_constr_vec = ['l', 'u']
parameter_update_structure = {
'PA': ParameterUpdateLogic(
Expand Down Expand Up @@ -316,6 +323,8 @@ def get_affine_map(self, p_id, param_prob, constraint_info: ConstraintInfo) -> A
affine_map = AffineMap()

if p_id == 'P':
if self.indices_obj is None: # problem is an LP
return None
affine_map.mapping = param_prob.reduced_P.reduced_mat
affine_map.indices = self.indices_obj
affine_map.shape = (self.n_var, self.n_var)
Expand All @@ -327,16 +336,19 @@ def get_affine_map(self, p_id, param_prob, constraint_info: ConstraintInfo) -> A
affine_map.mapping = param_prob.reduced_A.reduced_mat[
:constraint_info.n_data_constr_mat]
affine_map.indices = self.indices_constr[:constraint_info.n_data_constr_mat]
affine_map.mapping[affine_map.indices >= self.n_eq] *= -1
affine_map.shape = (self.n_eq + self.n_ineq, self.n_var)
elif p_id == 'l':
mapping_rows_eq = np.nonzero(self.indices_constr < self.n_eq)[0]
affine_map.mapping_rows = mapping_rows_eq[
mapping_rows_eq >= constraint_info.n_data_constr_mat] # mapping to the finite part of l
affine_map.sign = -1
affine_map.indices = self.indices_constr[affine_map.mapping_rows]
affine_map.shape = (self.n_eq, 1)
elif p_id == 'u':
affine_map.mapping_rows = np.arange(constraint_info.n_data_constr_mat,
constraint_info.n_data_constr)
affine_map.sign = np.vstack((-np.ones((self.n_eq, 1)), np.ones((self.n_ineq, 1))))
affine_map.indices = self.indices_constr[affine_map.mapping_rows]
affine_map.shape = (self.n_eq + self.n_ineq, 1)
else:
Expand All @@ -354,14 +366,15 @@ def augment_vector_parameter(self, p_id, vector_parameter):
class SCSInterface(SolverInterface):
solver_name = 'SCS'
solver_type = 'conic'
canon_p_ids = ['c', 'd', 'A', 'b']
canon_p_signs = [1, 1, 1, 1]
canon_p_ids = ['P', 'c', 'd', 'A', 'b']
canon_p_ids_constr_vec = ['b']
parameter_update_structure = {
'init': ParameterUpdateLogic(
update_pending_logic = UpdatePendingLogic(
['A'], extra_condition = '!{prefix}Scs_Work', extra_condition_operator = '||', functions_if_false = ['bc']
),
update_pending_logic = UpdatePendingLogic([], extra_condition = '!{prefix}Scs_Work', functions_if_false = ['PA']),
function_call = '{prefix}Scs_Work = scs_init(&{prefix}Scs_D, &{prefix}Scs_K, &{prefix}Canon_Settings)',
),
'PA': ParameterUpdateLogic(
update_pending_logic = UpdatePendingLogic(['P', 'A'], '||', functions_if_false = ['bc']),
function_call = '{prefix}Scs_Work = scs_init(&{prefix}Scs_D, &{prefix}Scs_K, &{prefix}Canon_Settings)',
),
'bc': ParameterUpdateLogic(
Expand Down Expand Up @@ -521,9 +534,11 @@ def generate_code(self, code_dir, solver_code_dir, cvxpygen_directory,
with open(os.path.join(code_dir, 'setup.py'), 'w') as f:
f.write(setup_text)

def declare_workspace(self, f, prefix) -> None:
f.write(f'\n// SCS matrix A\n')
write_struct_prot(f, f'{prefix}scs_A', 'ScsMatrix')
def declare_workspace(self, f, prefix, parameter_canon) -> None:
matrices = ['P', 'A'] if parameter_canon.quad_obj else ['A']
for m in matrices:
f.write(f'\n// SCS matrix {m}\n')
write_struct_prot(f, f'{prefix}scs_{m}', 'ScsMatrix')
f.write(f'\n// Struct containing SCS data\n')
write_struct_prot(f, f'{prefix}Scs_D', 'ScsData')
if self.canon_constants['qsize'] > 0:
Expand All @@ -545,20 +560,23 @@ def declare_workspace(self, f, prefix) -> None:
write_struct_prot(f, f'{prefix}Scs_Work', 'ScsWork*')


def define_workspace(self, f, prefix) -> None:
f.write(f'\n// SCS matrix A\n')
scs_A_fields = ['x', 'i', 'p', 'm', 'n']
scs_A_casts = ['(cpg_float *) ', '(cpg_int *) ', '(cpg_int *) ', '', '']
scs_A_values = [f'&{prefix}canon_A_x', f'&{prefix}canon_A_i',
f'&{prefix}canon_A_p', str(self.canon_constants['m']),
str(self.canon_constants['n'])]
write_struct_def(f, scs_A_fields, scs_A_casts, scs_A_values, f'{prefix}Scs_A', 'ScsMatrix')
def define_workspace(self, f, prefix, parameter_canon) -> None:
matrices = ['P', 'A'] if parameter_canon.quad_obj else ['A']
scs_PA_fields = ['x', 'i', 'p', 'm', 'n']
scs_PA_casts = ['(cpg_float *) ', '(cpg_int *) ', '(cpg_int *) ', '', '']
for m in matrices:
f.write(f'\n// SCS matrix {m}\n')
scs_PA_values = [f'&{prefix}canon_{m}_x', f'&{prefix}canon_{m}_i',
f'&{prefix}canon_{m}_p', str(self.canon_constants[('n' if m == 'P' else 'm')]),
str(self.canon_constants['n'])]
write_struct_def(f, scs_PA_fields, scs_PA_casts, scs_PA_values, f'{prefix}Scs_{m}', 'ScsMatrix')

f.write(f'\n// Struct containing SCS data\n')
scs_d_fields = ['m', 'n', 'A', 'P', 'b', 'c']
scs_d_casts = ['', '', '', '', '(cpg_float *) ', '(cpg_float *) ']
scs_d_values = [str(self.canon_constants['m']), str(self.canon_constants['n']),
f'&{prefix}Scs_A', 'SCS_NULL', f'&{prefix}canon_b', f'&{prefix}canon_c']
f'&{prefix}Scs_A', (f'&{prefix}Scs_P' if parameter_canon.quad_obj else 'SCS_NULL'),
f'&{prefix}canon_b', f'&{prefix}canon_c']
write_struct_def(f, scs_d_fields, scs_d_casts, scs_d_values, f'{prefix}Scs_D', 'ScsData')

if self.canon_constants['qsize'] > 0:
Expand Down Expand Up @@ -611,7 +629,6 @@ class ECOSInterface(SolverInterface):
solver_name = 'ECOS'
solver_type = 'conic'
canon_p_ids = ['c', 'd', 'A', 'b', 'G', 'h']
canon_p_signs = [1, 1, 1, 1, 1, 1]
canon_p_ids_constr_vec = ['b', 'h']
solve_function_call = '{prefix}ecos_flag = ECOS_solve({prefix}ecos_workspace)'

Expand Down Expand Up @@ -781,7 +798,7 @@ def generate_code(self, code_dir, solver_code_dir, cvxpygen_directory,
with open(os.path.join(code_dir, 'setup.py'), 'w') as f:
f.write(setup_text)

def declare_workspace(self, f, prefix) -> None:
def declare_workspace(self, f, prefix, parameter_canon) -> None:
if self.canon_constants['n_cones'] > 0:
f.write('\n// ECOS array of SOC dimensions\n')
write_vec_prot(f, self.canon_constants['q'], f'{prefix}ecos_q', 'cpg_int')
Expand All @@ -790,7 +807,7 @@ def declare_workspace(self, f, prefix) -> None:
f.write('\n// ECOS exit flag\n')
f.write(f'extern cpg_int {prefix}ecos_flag;\n')

def define_workspace(self, f, prefix) -> None:
def define_workspace(self, f, prefix, parameter_canon) -> None:
if self.canon_constants['n_cones'] > 0:
f.write('\n// ECOS array of SOC dimensions\n')
write_vec_def(f, self.canon_constants['q'], f'{prefix}ecos_q', 'cpg_int')
Expand All @@ -804,7 +821,6 @@ class ClarabelInterface(SolverInterface):
solver_name = 'Clarabel'
solver_type = 'conic'
canon_p_ids = ['P', 'q', 'd', 'A', 'b']
canon_p_signs = [1, 1, 1, -1, 1, 1]
canon_p_ids_constr_vec = ['b']

# header and source files
Expand Down Expand Up @@ -976,35 +992,8 @@ def generate_code(self, code_dir, solver_code_dir, cvxpygen_directory,
"extra_objects=[cpg_lib, os.path.join(cpg_dir, 'solver_code', 'rust_wrapper', 'target', 'debug', 'libclarabel_c.a')])")
with open(os.path.join(code_dir, 'setup.py'), 'w') as f:
f.write(setup_text)

def get_affine_map(self, p_id, param_prob, constraint_info: ConstraintInfo) -> AffineMap:
affine_map = AffineMap()

if p_id == 'P':
if self.indices_obj is None: # problem is an LP
return None
affine_map.mapping = param_prob.reduced_P.reduced_mat
affine_map.indices = self.indices_obj
affine_map.shape = (self.n_var, self.n_var)
elif p_id == 'q':
affine_map.mapping = param_prob.c[:-1]
elif p_id == 'd':
affine_map.mapping = param_prob.c[-1]
elif p_id == 'A':
affine_map.mapping = param_prob.reduced_A.reduced_mat[:constraint_info.n_data_constr_mat]
affine_map.indices = self.indices_constr[:constraint_info.n_data_constr_mat]
affine_map.shape = (self.n_eq, self.n_var) # only equality constraints
elif p_id == 'b': # provide 'mapping_rows' instead of 'mapping' for vector parameters since mapping will be decompressed
affine_map.mapping_rows = constraint_info.mapping_rows_eq[
constraint_info.mapping_rows_eq >= constraint_info.n_data_constr_mat]
affine_map.indices = self.indices_constr[affine_map.mapping_rows]
affine_map.shape = (self.n_eq, 1) # only equality constraints
else:
raise ValueError(f'Unknown parameter name {p_id}.')

return affine_map

def declare_workspace(self, f, prefix) -> None:
def declare_workspace(self, f, prefix, parameter_canon) -> None:
f.write('\n// Clarabel workspace\n')
f.write(f'extern ClarabelCscMatrix {prefix}P;\n')
f.write(f'extern ClarabelCscMatrix {prefix}A;\n')
Expand All @@ -1013,7 +1002,7 @@ def declare_workspace(self, f, prefix) -> None:
f.write(f'extern ClarabelDefaultSolver *{prefix}solver;\n')
f.write(f'extern ClarabelDefaultSolution {prefix}solution;\n')

def define_workspace(self, f, prefix) -> None:
def define_workspace(self, f, prefix, parameter_canon) -> None:
f.write('\n// Clarabel workspace\n')
f.write(f'ClarabelCscMatrix {prefix}P;\n')
f.write(f'ClarabelCscMatrix {prefix}A;\n')
Expand Down
Loading

0 comments on commit 34826c1

Please sign in to comment.