From 78494bc2ea411f6d094b88b975a5c780e6249db2 Mon Sep 17 00:00:00 2001 From: rsemenoff Date: Tue, 22 Nov 2022 16:18:25 -0800 Subject: [PATCH 1/7] neos --- cvxpy/problems/problem.py | 3078 ++++++++++++++++++------------------- 1 file changed, 1515 insertions(+), 1563 deletions(-) diff --git a/cvxpy/problems/problem.py b/cvxpy/problems/problem.py index 380e47ef06..841fadb2ce 100644 --- a/cvxpy/problems/problem.py +++ b/cvxpy/problems/problem.py @@ -1,1563 +1,1515 @@ -""" -Copyright 2013 Steven Diamond, 2017 Akshay Agrawal - -Licensed under the Apache License, Version 2.0 (the "License"); -you may not use this file except in compliance with the License. -You may obtain a copy of the License at - - http://www.apache.org/licenses/LICENSE-2.0 - -Unless required by applicable law or agreed to in writing, software -distributed under the License is distributed on an "AS IS" BASIS, -WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -See the License for the specific language governing permissions and -limitations under the License. -""" -from __future__ import annotations - -import time -import warnings -from collections import namedtuple -from typing import Dict, List, Optional, Union - -import numpy as np - -import cvxpy.utilities as u -import cvxpy.utilities.performance_utils as perf -from cvxpy import Constant, error -from cvxpy import settings as s -from cvxpy.atoms.atom import Atom -from cvxpy.constraints import Equality, Inequality, NonNeg, NonPos, Zero -from cvxpy.constraints.constraint import Constraint -from cvxpy.error import DPPError -from cvxpy.expressions import cvxtypes -from cvxpy.expressions.variable import Variable -from cvxpy.interface.matrix_utilities import scalar_value -from cvxpy.problems.objective import Maximize, Minimize -from cvxpy.reductions import InverseData -from cvxpy.reductions.chain import Chain -from cvxpy.reductions.dgp2dcp.dgp2dcp import Dgp2Dcp -from cvxpy.reductions.dqcp2dcp import dqcp2dcp -from cvxpy.reductions.eval_params import EvalParams -from cvxpy.reductions.flip_objective import FlipObjective -from cvxpy.reductions.solution import INF_OR_UNB_MESSAGE -from cvxpy.reductions.solvers import bisection -from cvxpy.reductions.solvers import defines as slv_def -from cvxpy.reductions.solvers.conic_solvers.conic_solver import ConicSolver -from cvxpy.reductions.solvers.defines import SOLVER_MAP_CONIC, SOLVER_MAP_QP -from cvxpy.reductions.solvers.qp_solvers.qp_solver import QpSolver -from cvxpy.reductions.solvers.solver import Solver -from cvxpy.reductions.solvers.solving_chain import (SolvingChain, - construct_solving_chain,) -from cvxpy.settings import SOLVERS -from cvxpy.utilities import debug_tools -from cvxpy.utilities.deterministic import unique_list - -SolveResult = namedtuple( - 'SolveResult', - ['opt_value', 'status', 'primal_values', 'dual_values']) - - -_COL_WIDTH = 79 -_HEADER = ( - '='*_COL_WIDTH + - '\n' + - ('CVXPY').center(_COL_WIDTH) + - '\n' + - ('v' + cvxtypes.version()).center(_COL_WIDTH) + - '\n' + - '='*_COL_WIDTH -) -_COMPILATION_STR = ( - '-'*_COL_WIDTH + - '\n' + - ('Compilation').center(_COL_WIDTH) + - '\n' + - '-'*_COL_WIDTH -) -_NUM_SOLVER_STR = ( - '-'*_COL_WIDTH + - '\n' + - ('Numerical solver').center(_COL_WIDTH) + - '\n' + - '-'*_COL_WIDTH -) -_FOOTER = ( - '-'*_COL_WIDTH + - '\n' + - ('Summary').center(_COL_WIDTH) + - '\n' + - '-'*_COL_WIDTH -) - - -class Cache: - def __init__(self) -> None: - self.key = None - self.solving_chain: Optional[SolvingChain] = None - self.param_prog = None - self.inverse_data: Optional[InverseData] = None - - def invalidate(self) -> None: - self.key = None - self.solving_chain = None - self.param_prog = None - self.inverse_data = None - - def make_key(self, solver, gp, ignore_dpp): - return (solver, gp, ignore_dpp) - - def gp(self): - return self.key is not None and self.key[1] - - -def _validate_constraint(constraint): - if isinstance(constraint, Constraint): - return constraint - elif isinstance(constraint, bool): - # replace `True` or `False` values with equivalent Expressions. - return (Constant(0) <= Constant(1) if constraint else - Constant(1) <= Constant(0)) - else: - raise ValueError("Problem has an invalid constraint of type %s" % - type(constraint)) - - -class Problem(u.Canonical): - """A convex optimization problem. - - Problems are immutable, save for modification through the specification - of :class:`~cvxpy.expressions.constants.parameters.Parameter` - - Arguments - --------- - objective : Minimize or Maximize - The problem's objective. - constraints : list - The constraints on the problem variables. - """ - - # The solve methods available. - REGISTERED_SOLVE_METHODS = {} - - def __init__( - self, objective: Union[Minimize, Maximize], constraints: Optional[List[Constraint]] = None - ) -> None: - if constraints is None: - constraints = [] - # Check that objective is Minimize or Maximize. - if not isinstance(objective, (Minimize, Maximize)): - raise error.DCPError("Problem objective must be Minimize or Maximize.") - # Constraints and objective are immutable. - self._objective = objective - # Raise warning if objective has too many subexpressions. - if debug_tools.node_count(self._objective) >= debug_tools.MAX_NODES: - warnings.warn("Objective contains too many subexpressions. " - "Consider vectorizing your CVXPY code to speed up compilation.") - self._constraints = [_validate_constraint(c) for c in constraints] - # Raise warning if constraint has too many subexpressions. - for i, constraint in enumerate(self._constraints): - if debug_tools.node_count(constraint) >= debug_tools.MAX_NODES: - warnings.warn(f"Constraint #{i} contains too many subexpressions. " - "Consider vectorizing your CVXPY code to speed up compilation.") - - self._value = None - self._status: Optional[str] = None - self._solution = None - self._cache = Cache() - self._solver_cache = {} - # Information about the shape of the problem and its constituent parts - self._size_metrics: Optional["SizeMetrics"] = None - # Benchmarks reported by the solver: - self._solver_stats: Optional["SolverStats"] = None - self._compilation_time: Optional[float] = None - self._solve_time: Optional[float] = None - self.args = [self._objective, self._constraints] - - @property - def value(self): - """float : The value from the last time the problem was solved - (or None if not solved). - """ - if self._value is None: - return None - else: - return scalar_value(self._value) - - @property - def status(self) -> str: - """str : The status from the last time the problem was solved; one - of optimal, infeasible, or unbounded (with or without - suffix inaccurate). - """ - return self._status - - @property - def solution(self): - """Solution : The solution from the last time the problem was solved. - """ - return self._solution - - @property - def objective(self) -> Union[Minimize, Maximize]: - """Minimize or Maximize : The problem's objective. - - Note that the objective cannot be reassigned after creation, - and modifying the objective after creation will result in - undefined behavior. - """ - return self._objective - - @property - def constraints(self) -> List[Constraint]: - """A shallow copy of the problem's constraints. - - Note that constraints cannot be reassigned, appended to, or otherwise - modified after creation, except through parameters. - """ - return self._constraints[:] - - @property - def param_dict(self): - """ - Expose all parameters as a dictionary - """ - return {parameters.name(): parameters for parameters in self.parameters()} - - @property - def var_dict(self) -> Dict[str, Variable]: - """ - Expose all variables as a dictionary - """ - return {variable.name(): variable for variable in self.variables()} - - @perf.compute_once - def is_dcp(self, dpp: bool = False) -> bool: - """Does the problem satisfy DCP rules? - - Arguments - --------- - dpp : bool, optional - If True, enforce the disciplined parametrized programming (DPP) - ruleset; only relevant when the problem involves Parameters. - DPP is a mild restriction of DCP. When a problem involving - Parameters is DPP, subsequent solves can be much faster than - the first one. For more information, consult the documentation at - - https://www.cvxpy.org/tutorial/advanced/index.html#disciplined-parametrized-programming - - Returns - ------- - bool - True if the Expression is DCP, False otherwise. - """ - return all( - expr.is_dcp(dpp) for expr in self.constraints + [self.objective]) - - @perf.compute_once - def is_dgp(self, dpp: bool = False) -> bool: - """Does the problem satisfy DGP rules? - - Arguments - --------- - dpp : bool, optional - If True, enforce the disciplined parametrized programming (DPP) - ruleset; only relevant when the problem involves Parameters. - DPP is a mild restriction of DGP. When a problem involving - Parameters is DPP, subsequent solves can be much faster than - the first one. For more information, consult the documentation at - - https://www.cvxpy.org/tutorial/advanced/index.html#disciplined-parametrized-programming - - Returns - ------- - bool - True if the Expression is DGP, False otherwise. - """ - return all( - expr.is_dgp(dpp) for expr in self.constraints + [self.objective]) - - @perf.compute_once - def is_dqcp(self) -> bool: - """Does the problem satisfy the DQCP rules? - """ - return all( - expr.is_dqcp() for expr in self.constraints + [self.objective]) - - @perf.compute_once - def is_dpp(self, context: str = 'dcp') -> bool: - """Does the problem satisfy DPP rules? - - DPP is a mild restriction of DGP. When a problem involving - Parameters is DPP, subsequent solves can be much faster than - the first one. For more information, consult the documentation at - - https://www.cvxpy.org/tutorial/advanced/index.html#disciplined-parametrized-programming - - Arguments - --------- - context : str - Whether to check DPP-compliance for DCP or DGP; ``context`` should - be either ``'dcp'`` or ``'dgp'``. Calling ``problem.is_dpp('dcp')`` - is equivalent to ``problem.is_dcp(dpp=True)``, and - `problem.is_dpp('dgp')`` is equivalent to - `problem.is_dgp(dpp=True)`. - - Returns - ------- - bool - Whether the problem satisfies the DPP rules. - """ - if context.lower() == 'dcp': - return self.is_dcp(dpp=True) - elif context.lower() == 'dgp': - return self.is_dgp(dpp=True) - else: - raise ValueError("Unsupported context ", context) - - @perf.compute_once - def is_qp(self) -> bool: - """Is problem a quadratic program? - """ - for c in self.constraints: - if not (isinstance(c, (Equality, Zero)) or c.args[0].is_pwl()): - return False - for var in self.variables(): - if var.is_psd() or var.is_nsd(): - return False - return (self.is_dcp() and self.objective.args[0].is_qpwa()) - - @perf.compute_once - def is_mixed_integer(self) -> bool: - return any(v.attributes['boolean'] or v.attributes['integer'] - for v in self.variables()) - - @perf.compute_once - def variables(self) -> List[Variable]: - """Accessor method for variables. - - Returns - ------- - list of :class:`~cvxpy.expressions.variable.Variable` - A list of the variables in the problem. - """ - vars_ = self.objective.variables() - for constr in self.constraints: - vars_ += constr.variables() - return unique_list(vars_) - - @perf.compute_once - def parameters(self): - """Accessor method for parameters. - - Returns - ------- - list of :class:`~cvxpy.expressions.constants.parameter.Parameter` - A list of the parameters in the problem. - """ - params = self.objective.parameters() - for constr in self.constraints: - params += constr.parameters() - return unique_list(params) - - @perf.compute_once - def constants(self) -> List[Constant]: - """Accessor method for constants. - - Returns - ------- - list of :class:`~cvxpy.expressions.constants.constant.Constant` - A list of the constants in the problem. - """ - const_dict = {} - constants_ = self.objective.constants() - for constr in self.constraints: - constants_ += constr.constants() - # Note that numpy matrices are not hashable, so we use the built-in - # function "id" - const_dict = {id(constant): constant for constant in constants_} - return list(const_dict.values()) - - def atoms(self) -> List[Atom]: - """Accessor method for atoms. - - Returns - ------- - list of :class:`~cvxpy.atoms.Atom` - A list of the atom types in the problem; note that this list - contains classes, not instances. - """ - atoms = self.objective.atoms() - for constr in self.constraints: - atoms += constr.atoms() - return unique_list(atoms) - - @property - def size_metrics(self) -> "SizeMetrics": - """:class:`~cvxpy.problems.problem.SizeMetrics` : Information about the problem's size. - """ - if self._size_metrics is None: - self._size_metrics = SizeMetrics(self) - return self._size_metrics - - @property - def solver_stats(self) -> "SolverStats": - """:class:`~cvxpy.problems.problem.SolverStats` : Information returned by the solver. - """ - return self._solver_stats - - def solve(self, *args, **kwargs): - """Compiles and solves the problem using the specified method. - - Populates the :code:`status` and :code:`value` attributes on the - problem object as a side-effect. - - Arguments - --------- - solver : str, optional - The solver to use. For example, 'ECOS', 'SCS', or 'OSQP'. - verbose : bool, optional - Overrides the default of hiding solver output, and prints - logging information describing CVXPY's compilation process. - gp : bool, optional - If True, parses the problem as a disciplined geometric program - instead of a disciplined convex program. - qcp : bool, optional - If True, parses the problem as a disciplined quasiconvex program - instead of a disciplined convex program. - requires_grad : bool, optional - Makes it possible to compute gradients of a solution with respect to - Parameters by calling ``problem.backward()`` after solving, or to - compute perturbations to the variables given perturbations to Parameters by - calling ``problem.derivative()``. - - Gradients are only supported for DCP and DGP problems, not - quasiconvex problems. When computing gradients (i.e., when - this argument is True), the problem must satisfy the DPP rules. - enforce_dpp : bool, optional - When True, a DPPError will be thrown when trying to solve a non-DPP - problem (instead of just a warning). Only relevant for problems - involving Parameters. Defaults to False. - ignore_dpp : bool, optional - When True, DPP problems will be treated as non-DPP, - which may speed up compilation. Defaults to False. - method : function, optional - A custom solve method to use. - kwargs : keywords, optional - Additional solver specific arguments. See Notes below. - - Notes - ------ - CVXPY interfaces with a wide range of solvers; the algorithms used by these solvers - have arguments relating to stopping criteria, and strategies to improve solution quality. - - There is no one choice of arguments which is perfect for every problem. If you are not - getting satisfactory results from a solver, you can try changing its arguments. The - exact way this is done depends on the specific solver. Here are some examples: - - :: - - prob.solve(solver='ECOS', abstol=1e-6) - prob.solve(solver='OSQP', max_iter=10000). - mydict = {"MSK_DPAR_INTPNT_CO_TOL_NEAR_REL": 10} - prob.solve(solver='MOSEK', mosek_params=mydict). - - You should refer to CVXPY's web documentation for details on how to pass solver - solver arguments, available at - - https://www.cvxpy.org/tutorial/advanced/index.html#setting-solver-options - - Returns - ------- - float - The optimal value for the problem, or a string indicating - why the problem could not be solved. - - Raises - ------ - cvxpy.error.DCPError - Raised if the problem is not DCP and `gp` is False. - cvxpy.error.DGPError - Raised if the problem is not DGP and `gp` is True. - cvxpy.error.DPPError - Raised if DPP settings are invalid. - cvxpy.error.SolverError - Raised if no suitable solver exists among the installed solvers, - or if an unanticipated error is encountered. - """ - func_name = kwargs.pop("method", None) - if func_name is not None: - solve_func = Problem.REGISTERED_SOLVE_METHODS[func_name] - else: - solve_func = Problem._solve - return solve_func(self, *args, **kwargs) - - @classmethod - def register_solve(cls, name: str, func) -> None: - """Adds a solve method to the Problem class. - - Arguments - --------- - name : str - The keyword for the method. - func : function - The function that executes the solve method. This function must - take as its first argument the problem instance to solve. - """ - cls.REGISTERED_SOLVE_METHODS[name] = func - - def get_problem_data( - self, solver, - gp: bool = False, - enforce_dpp: bool = False, - ignore_dpp: bool = False, - verbose: bool = False, - canon_backend: str | None = None, - solver_opts: Optional[dict] = None - ): - """Returns the problem data used in the call to the solver. - - When a problem is solved, CVXPY creates a chain of reductions enclosed - in a :class:`~cvxpy.reductions.solvers.solving_chain.SolvingChain`, - and compiles it to some low-level representation that is - compatible with the targeted solver. This method returns that low-level - representation. - - For some solving chains, this low-level representation is a dictionary - that contains exactly those arguments that were supplied to the solver; - however, for other solving chains, the data is an intermediate - representation that is compiled even further by the solver interfaces. - - A solution to the equivalent low-level problem can be obtained via the - data by invoking the `solve_via_data` method of the returned solving - chain, a thin wrapper around the code external to CVXPY that further - processes and solves the problem. Invoke the unpack_results method - to recover a solution to the original problem. - - For example: - - :: - - objective = ... - constraints = ... - problem = cp.Problem(objective, constraints) - data, chain, inverse_data = problem.get_problem_data(cp.SCS) - # calls SCS using `data` - soln = chain.solve_via_data(problem, data) - # unpacks the solution returned by SCS into `problem` - problem.unpack_results(soln, chain, inverse_data) - - Alternatively, the `data` dictionary returned by this method - contains enough information to bypass CVXPY and call the solver - directly. - - For example: - - :: - - problem = cp.Problem(objective, constraints) - data, _, _ = problem.get_problem_data(cp.SCS) - - import scs - probdata = { - 'A': data['A'], - 'b': data['b'], - 'c': data['c'], - } - cone_dims = data['dims'] - cones = { - "f": cone_dims.zero, - "l": cone_dims.nonneg, - "q": cone_dims.soc, - "ep": cone_dims.exp, - "s": cone_dims.psd, - } - soln = scs.solve(data, cones) - - The structure of the data dict that CVXPY returns depends on the - solver. For details, consult the solver interfaces in - `cvxpy/reductions/solvers`. - - Arguments - --------- - solver : str - The solver the problem data is for. - gp : bool, optional - If True, then parses the problem as a disciplined geometric program - instead of a disciplined convex program. - enforce_dpp : bool, optional - When True, a DPPError will be thrown when trying to parse a non-DPP - problem (instead of just a warning). Defaults to False. - ignore_dpp : bool, optional - When True, DPP problems will be treated as non-DPP, - which may speed up compilation. Defaults to False. - canon_backend : str, optional - 'CPP' (default) | 'SCIPY' - Specifies which backend to use for canonicalization, which can affect - compilation time. Defaults to None, i.e., selecting the default - backend. - verbose : bool, optional - If True, print verbose output related to problem compilation. - solver_opts : dict, optional - A dict of options that will be passed to the specific solver. - In general, these options will override any default settings - imposed by cvxpy. - - Returns - ------- - dict or object - lowest level representation of problem - SolvingChain - The solving chain that created the data. - list - The inverse data generated by the chain. - - Raises - ------ - cvxpy.error.DPPError - Raised if DPP settings are invalid. - """ - - # Invalid DPP setting. - # Must be checked here to avoid cache issues. - if enforce_dpp and ignore_dpp: - raise DPPError("Cannot set enforce_dpp = True and ignore_dpp = True.") - - start = time.time() - # Cache includes ignore_dpp because it alters compilation. - key = self._cache.make_key(solver, gp, ignore_dpp) - if key != self._cache.key: - self._cache.invalidate() - solving_chain = self._construct_chain( - solver=solver, gp=gp, - enforce_dpp=enforce_dpp, - ignore_dpp=ignore_dpp, - canon_backend=canon_backend, - solver_opts=solver_opts) - self._cache.key = key - self._cache.solving_chain = solving_chain - self._solver_cache = {} - else: - solving_chain = self._cache.solving_chain - - if verbose: - print(_COMPILATION_STR) - - if self._cache.param_prog is not None: - # fast path, bypasses application of reductions - if verbose: - s.LOGGER.info( - 'Using cached ASA map, for faster compilation ' - '(bypassing reduction chain).') - if gp: - dgp2dcp = self._cache.solving_chain.get(Dgp2Dcp) - # Parameters in the param cone prog are the logs - # of parameters in the original problem (with one exception: - # parameters appearing as exponents (in power and gmatmul - # atoms) are unchanged. - old_params_to_new_params = dgp2dcp.canon_methods._parameters - for param in self.parameters(): - - if param in old_params_to_new_params: - old_params_to_new_params[param].value = np.log( - param.value) - - data, solver_inverse_data = solving_chain.solver.apply( - self._cache.param_prog) - inverse_data = self._cache.inverse_data + [solver_inverse_data] - self._compilation_time = time.time() - start - if verbose: - s.LOGGER.info( - 'Finished problem compilation ' - '(took %.3e seconds).', self._compilation_time) - else: - if verbose: - solver_name = solving_chain.reductions[-1].name() - reduction_chain_str = ' -> '.join( - type(r).__name__ for r in solving_chain.reductions) - s.LOGGER.info( - 'Compiling problem (target solver=%s).', solver_name) - s.LOGGER.info('Reduction chain: %s', reduction_chain_str) - data, inverse_data = solving_chain.apply(self, verbose) - safe_to_cache = ( - isinstance(data, dict) - and s.PARAM_PROB in data - and not any(isinstance(reduction, EvalParams) - for reduction in solving_chain.reductions) - ) - self._compilation_time = time.time() - start - if verbose: - s.LOGGER.info( - 'Finished problem compilation ' - '(took %.3e seconds).', self._compilation_time) - if safe_to_cache: - if verbose and self.parameters(): - s.LOGGER.info( - '(Subsequent compilations of this problem, using the ' - 'same arguments, should ' 'take less time.)') - self._cache.param_prog = data[s.PARAM_PROB] - # the last datum in inverse_data corresponds to the solver, - # so we shouldn't cache it - self._cache.inverse_data = inverse_data[:-1] - return data, solving_chain, inverse_data - - def _find_candidate_solvers(self, - solver=None, - gp: bool = False): - """ - Find candidate solvers for the current problem. If solver - is not None, it checks if the specified solver is compatible - with the problem passed. - - Arguments - --------- - solver : Union[string, Solver, None] - The name of the solver with which to solve the problem or an - instance of a custom solver. If no solver is supplied - (i.e., if solver is None), then the targeted solver may be any - of those that are installed. If the problem is variable-free, - then this parameter is ignored. - gp : bool - If True, the problem is parsed as a Disciplined Geometric Program - instead of as a Disciplined Convex Program. - - Returns - ------- - dict - A dictionary of compatible solvers divided in `qp_solvers` - and `conic_solvers`. - - Raises - ------ - cvxpy.error.SolverError - Raised if the problem is not DCP and `gp` is False. - cvxpy.error.DGPError - Raised if the problem is not DGP and `gp` is True. - """ - candidates = {'qp_solvers': [], - 'conic_solvers': []} - if isinstance(solver, Solver): - return self._add_custom_solver_candidates(solver) - - if solver is not None: - if solver not in slv_def.INSTALLED_SOLVERS: - raise error.SolverError("The solver %s is not installed." % solver) - if solver in slv_def.CONIC_SOLVERS: - candidates['conic_solvers'] += [solver] - if solver in slv_def.QP_SOLVERS: - candidates['qp_solvers'] += [solver] - else: - candidates['qp_solvers'] = [s for s in slv_def.INSTALLED_SOLVERS - if s in slv_def.QP_SOLVERS] - candidates['conic_solvers'] = [] - # ECOS_BB can only be called explicitly. - for slv in slv_def.INSTALLED_SOLVERS: - if slv in slv_def.CONIC_SOLVERS and slv != s.ECOS_BB: - candidates['conic_solvers'].append(slv) - - # If gp we must have only conic solvers - if gp: - if solver is not None and solver not in slv_def.CONIC_SOLVERS: - raise error.SolverError( - "When `gp=True`, `solver` must be a conic solver " - "(received '%s'); try calling " % solver + - " `solve()` with `solver=cvxpy.ECOS`." - ) - elif solver is None: - candidates['qp_solvers'] = [] # No QP solvers allowed - - if self.is_mixed_integer(): - # ECOS_BB must be called explicitly. - if slv_def.INSTALLED_MI_SOLVERS == [s.ECOS_BB] and solver != s.ECOS_BB: - msg = """ - - You need a mixed-integer solver for this model. Refer to the documentation - https://www.cvxpy.org/tutorial/advanced/index.html#mixed-integer-programs - for discussion on this topic. - - Quick fix 1: if you install the python package CVXOPT (pip install cvxopt), - then CVXPY can use the open-source mixed-integer linear programming - solver `GLPK`. If your problem is nonlinear then you can install SCIP - (pip install pyscipopt). - - Quick fix 2: you can explicitly specify solver='ECOS_BB'. This may result - in incorrect solutions and is not recommended. - """ - raise error.SolverError(msg) - # TODO: provide a useful error message when the problem is nonlinear but - # the only installed mixed-integer solvers are MILP solvers (e.g., GLPK_MI). - candidates['qp_solvers'] = [ - s for s in candidates['qp_solvers'] - if slv_def.SOLVER_MAP_QP[s].MIP_CAPABLE] - candidates['conic_solvers'] = [ - s for s in candidates['conic_solvers'] - if slv_def.SOLVER_MAP_CONIC[s].MIP_CAPABLE] - if not candidates['conic_solvers'] and \ - not candidates['qp_solvers']: - raise error.SolverError( - "Problem is mixed-integer, but candidate " - "QP/Conic solvers (%s) are not MIP-capable." % - (candidates['qp_solvers'] + - candidates['conic_solvers'])) - - return candidates - - def _add_custom_solver_candidates(self, custom_solver: Solver): - """ - Returns a list of candidate solvers where custom_solver is the only potential option. - - Arguments - --------- - custom_solver : Solver - - Returns - ------- - dict - A dictionary of compatible solvers divided in `qp_solvers` - and `conic_solvers`. - - Raises - ------ - cvxpy.error.SolverError - Raised if the name of the custom solver conflicts with the name of some officially - supported solver - """ - if custom_solver.name() in SOLVERS: - message = "Custom solvers must have a different name than the officially supported ones" - raise error.SolverError(message) - - candidates = {'qp_solvers': [], 'conic_solvers': []} - if not self.is_mixed_integer() or custom_solver.MIP_CAPABLE: - if isinstance(custom_solver, QpSolver): - SOLVER_MAP_QP[custom_solver.name()] = custom_solver - candidates['qp_solvers'] = [custom_solver.name()] - elif isinstance(custom_solver, ConicSolver): - SOLVER_MAP_CONIC[custom_solver.name()] = custom_solver - candidates['conic_solvers'] = [custom_solver.name()] - return candidates - - def _construct_chain( - self, - solver: Optional[str] = None, - gp: bool = False, - enforce_dpp: bool = False, - ignore_dpp: bool = False, - canon_backend: str | None = None, - solver_opts: Optional[dict] = None - ) -> SolvingChain: - """ - Construct the chains required to reformulate and solve the problem. - - In particular, this function - - # finds the candidate solvers - # constructs the solving chain that performs the - numeric reductions and solves the problem. - - Arguments - --------- - solver : str, optional - The solver to use. Defaults to ECOS. - gp : bool, optional - If True, the problem is parsed as a Disciplined Geometric Program - instead of as a Disciplined Convex Program. - enforce_dpp : bool, optional - Whether to error on DPP violations. - ignore_dpp : bool, optional - When True, DPP problems will be treated as non-DPP, - which may speed up compilation. Defaults to False. - canon_backend : str, optional - 'CPP' (default) | 'SCIPY' - Specifies which backend to use for canonicalization, which can affect - compilation time. Defaults to None, i.e., selecting the default - backend. - solver_opts: dict, optional - Additional arguments to pass to the solver. - - Returns - ------- - A solving chain - """ - candidate_solvers = self._find_candidate_solvers(solver=solver, gp=gp) - self._sort_candidate_solvers(candidate_solvers) - return construct_solving_chain(self, candidate_solvers, gp=gp, - enforce_dpp=enforce_dpp, - ignore_dpp=ignore_dpp, - canon_backend=canon_backend, - solver_opts=solver_opts) - - @staticmethod - def _sort_candidate_solvers(solvers) -> None: - """Sorts candidate solvers lists according to slv_def.CONIC_SOLVERS/QP_SOLVERS - - Arguments - --------- - candidates : dict - Dictionary of candidate solvers divided in qp_solvers - and conic_solvers - Returns - ------- - None - """ - if len(solvers['conic_solvers']) > 1: - solvers['conic_solvers'] = sorted( - solvers['conic_solvers'], key=lambda s: slv_def.CONIC_SOLVERS.index(s) - ) - if len(solvers['qp_solvers']) > 1: - solvers['qp_solvers'] = sorted( - solvers['qp_solvers'], key=lambda s: slv_def.QP_SOLVERS.index(s) - ) - - def _invalidate_cache(self) -> None: - self._cache_key = None - self._solving_chain = None - self._param_prog = None - self._inverse_data = None - - def _solve(self, - solver: str = None, - warm_start: bool = True, - verbose: bool = False, - gp: bool = False, - qcp: bool = False, - requires_grad: bool = False, - enforce_dpp: bool = False, - ignore_dpp: bool = False, - canon_backend: str | None = None, - **kwargs): - """Solves a DCP compliant optimization problem. - - Saves the values of primal and dual variables in the variable - and constraint objects, respectively. - - Arguments - --------- - solver : str, optional - The solver to use. Defaults to ECOS. - warm_start : bool, optional - Should the previous solver result be used to warm start? - verbose : bool, optional - Overrides the default of hiding solver output. - gp : bool, optional - If True, parses the problem as a disciplined geometric program. - qcp : bool, optional - If True, parses the problem as a disciplined quasiconvex program. - requires_grad : bool, optional - Makes it possible to compute gradients with respect to - parameters by calling `backward()` after solving, or to compute - perturbations to the variables by calling `derivative()`. When - True, the solver must be SCS, and dqcp must be False. - A DPPError is thrown when problem is not DPP. - enforce_dpp : bool, optional - When True, a DPPError will be thrown when trying to solve a non-DPP - problem (instead of just a warning). Defaults to False. - ignore_dpp : bool, optional - When True, DPP problems will be treated as non-DPP, - which may speed up compilation. Defaults to False. - canon_backend : str, optional - 'CPP' (default) | 'SCIPY' - Specifies which backend to use for canonicalization, which can affect - compilation time. Defaults to None, i.e., selecting the default - backend. - kwargs : dict, optional - A dict of options that will be passed to the specific solver. - In general, these options will override any default settings - imposed by cvxpy. - - Returns - ------- - float - The optimal value for the problem, or a string indicating - why the problem could not be solved. - """ - - if verbose: - print(_HEADER) - - for parameter in self.parameters(): - if parameter.value is None: - raise error.ParameterError( - "A Parameter (whose name is '%s') does not have a value " - "associated with it; all Parameter objects must have " - "values before solving a problem." % parameter.name()) - - if verbose: - n_variables = sum(np.prod(v.shape) for v in self.variables()) - n_parameters = sum(np.prod(p.shape) for p in self.parameters()) - s.LOGGER.info( - 'Your problem has %d variables, ' - '%d constraints, and ' '%d parameters.', - n_variables, len(self.constraints), n_parameters) - curvatures = [] - if self.is_dcp(): - curvatures.append('DCP') - if self.is_dgp(): - curvatures.append('DGP') - if self.is_dqcp(): - curvatures.append('DQCP') - s.LOGGER.info( - 'It is compliant with the following grammars: %s', - ', '.join(curvatures)) - if n_parameters == 0: - s.LOGGER.info( - '(If you need to solve this problem multiple times, ' - 'but with different data, consider using parameters.)') - s.LOGGER.info( - 'CVXPY will first compile your problem; then, it will ' - 'invoke a numerical solver to obtain a solution.') - - if requires_grad: - dpp_context = 'dgp' if gp else 'dcp' - if qcp: - raise ValueError("Cannot compute gradients of DQCP problems.") - elif not self.is_dpp(dpp_context): - raise error.DPPError("Problem is not DPP (when requires_grad " - "is True, problem must be DPP).") - elif solver is not None and solver not in [s.SCS, s.DIFFCP]: - raise ValueError("When requires_grad is True, the only " - "supported solver is SCS " - "(received %s)." % solver) - elif s.DIFFCP not in slv_def.INSTALLED_SOLVERS: - raise ModuleNotFoundError( - "The Python package diffcp must be installed to " - "differentiate through problems. Please follow the " - "installation instructions at " - "https://github.com/cvxgrp/diffcp") - else: - solver = s.DIFFCP - else: - if gp and qcp: - raise ValueError("At most one of `gp` and `qcp` can be True.") - if qcp and not self.is_dcp(): - if not self.is_dqcp(): - raise error.DQCPError("The problem is not DQCP.") - if verbose: - s.LOGGER.info( - 'Reducing DQCP problem to a one-parameter ' - 'family of DCP problems, for bisection.') - reductions = [dqcp2dcp.Dqcp2Dcp()] - start = time.time() - if type(self.objective) == Maximize: - reductions = [FlipObjective()] + reductions - # flip objective flips the sign of the objective - low, high = kwargs.get("low"), kwargs.get("high") - if high is not None: - kwargs["low"] = high * -1 - if low is not None: - kwargs["high"] = low * -1 - chain = Chain(problem=self, reductions=reductions) - soln = bisection.bisect( - chain.reduce(), solver=solver, verbose=verbose, **kwargs) - self.unpack(chain.retrieve(soln)) - return self.value - - data, solving_chain, inverse_data = self.get_problem_data( - solver, gp, enforce_dpp, ignore_dpp, verbose, canon_backend, kwargs - ) - - if verbose: - print(_NUM_SOLVER_STR) - s.LOGGER.info( - 'Invoking solver %s to obtain a solution.', - solving_chain.reductions[-1].name()) - start = time.time() - solution = solving_chain.solve_via_data( - self, data, warm_start, verbose, kwargs) - end = time.time() - self._solve_time = end - start - self.unpack_results(solution, solving_chain, inverse_data) - if verbose: - print(_FOOTER) - s.LOGGER.info('Problem status: %s', self.status) - val = self.value if self.value is not None else np.NaN - s.LOGGER.info('Optimal value: %.3e', val) - s.LOGGER.info('Compilation took %.3e seconds', self._compilation_time) - s.LOGGER.info( - 'Solver (including time spent in interface) took ' - '%.3e seconds', self._solve_time) - return self.value - - def backward(self) -> None: - """Compute the gradient of a solution with respect to Parameters. - - This method differentiates through the solution map of the problem, - obtaining the gradient of a solution with respect to the Parameters. - In other words, it calculates the sensitivities of the Parameters - with respect to perturbations in the optimal Variable values. This - can be useful for integrating CVXPY into automatic differentiation - toolkits. - - ``backward()`` populates the ``gradient`` attribute of each Parameter - in the problem as a side-effect. It can only be called after calling - ``solve()`` with ``requires_grad=True``. - - Below is a simple example: - - :: - - import cvxpy as cp - import numpy as np - - p = cp.Parameter() - x = cp.Variable() - quadratic = cp.square(x - 2 * p) - problem = cp.Problem(cp.Minimize(quadratic), [x >= 0]) - p.value = 3.0 - problem.solve(requires_grad=True, eps=1e-10) - # backward() populates the gradient attribute of the parameters - problem.backward() - # Because x* = 2 * p, dx*/dp = 2 - np.testing.assert_allclose(p.gradient, 2.0) - - In the above example, the gradient could easily be computed by hand. - The ``backward()`` is useful because for almost all problems, the - gradient cannot be computed analytically. - - This method can be used to differentiate through any DCP or DGP - problem, as long as the problem is DPP compliant (i.e., - ``problem.is_dcp(dpp=True)`` or ``problem.is_dgp(dpp=True)`` evaluates to - ``True``). - - This method uses the chain rule to evaluate the gradients of a - scalar-valued function of the Variables with respect to the Parameters. - For example, let x be a variable and p a Parameter; x and p might be - scalars, vectors, or matrices. Let f be a scalar-valued function, with - z = f(x). Then this method computes dz/dp = (dz/dx) (dx/p). dz/dx - is chosen as the all-ones vector by default, corresponding to - choosing f to be the sum function. You can specify a custom value for - dz/dx by setting the ``gradient`` attribute on your variables. For example, - - :: - - import cvxpy as cp - import numpy as np - - - b = cp.Parameter() - x = cp.Variable() - quadratic = cp.square(x - 2 * b) - problem = cp.Problem(cp.Minimize(quadratic), [x >= 0]) - b.value = 3. - problem.solve(requires_grad=True, eps=1e-10) - x.gradient = 4. - problem.backward() - # dz/dp = dz/dx dx/dp = 4. * 2. == 8. - np.testing.assert_allclose(b.gradient, 8.) - - The ``gradient`` attribute on a variable can also be interpreted as a - perturbation to its optimal value. - - Raises - ------ - ValueError - if solve was not called with ``requires_grad=True`` - SolverError - if the problem is infeasible or unbounded - """ - if s.DIFFCP not in self._solver_cache: - raise ValueError("backward can only be called after calling " - "solve with `requires_grad=True`") - elif self.status not in s.SOLUTION_PRESENT: - raise error.SolverError("Backpropagating through " - "infeasible/unbounded problems is not " - "yet supported. Please file an issue on " - "Github if you need this feature.") - - # TODO(akshayka): Backpropagate through dual variables as well. - backward_cache = self._solver_cache[s.DIFFCP] - DT = backward_cache["DT"] - zeros = np.zeros(backward_cache["s"].shape) - del_vars = {} - - gp = self._cache.gp() - for variable in self.variables(): - if variable.gradient is None: - del_vars[variable.id] = np.ones(variable.shape) - else: - del_vars[variable.id] = np.asarray(variable.gradient, - dtype=np.float64) - if gp: - # x_gp = exp(x_cone_program), - # dx_gp/d x_cone_program = exp(x_cone_program) = x_gp - del_vars[variable.id] *= variable.value - - dx = self._cache.param_prog.split_adjoint(del_vars) - start = time.time() - dA, db, dc = DT(dx, zeros, zeros) - end = time.time() - backward_cache['DT_TIME'] = end - start - dparams = self._cache.param_prog.apply_param_jac(dc, -dA, db) - - if not gp: - for param in self.parameters(): - param.gradient = dparams[param.id] - else: - dgp2dcp = self._cache.solving_chain.get(Dgp2Dcp) - old_params_to_new_params = dgp2dcp.canon_methods._parameters - for param in self.parameters(): - # Note: if param is an exponent in a power or gmatmul atom, - # then the parameter passes through unchanged to the DCP - # program; if the param is also used elsewhere (not as an - # exponent), then param will also be in - # old_params_to_new_params. Therefore, param.gradient = - # dparams[param.id] (or 0) + 1/param*dparams[new_param.id] - # - # Note that param.id is in dparams if and only if - # param was used as an exponent (because this means that - # the parameter entered the DCP problem unchanged.) - grad = 0.0 if param.id not in dparams else dparams[param.id] - if param in old_params_to_new_params: - new_param = old_params_to_new_params[param] - # new_param.value == log(param), apply chain rule - grad += (1.0 / param.value) * dparams[new_param.id] - param.gradient = grad - - def derivative(self) -> None: - """Apply the derivative of the solution map to perturbations in the Parameters - - This method applies the derivative of the solution map to perturbations - in the Parameters to obtain perturbations in the optimal values of the - Variables. In other words, it tells you how the optimal values of the - Variables would be changed by small changes to the Parameters. - - You can specify perturbations in a Parameter by setting its ``delta`` - attribute (if unspecified, the perturbation defaults to 0). - - This method populates the ``delta`` attribute of the Variables as a - side-effect. - - This method can only be called after calling ``solve()`` with - ``requires_grad=True``. It is compatible with both DCP and DGP - problems (that are also DPP-compliant). - - Below is a simple example: - - :: - - import cvxpy as cp - import numpy as np - - p = cp.Parameter() - x = cp.Variable() - quadratic = cp.square(x - 2 * p) - problem = cp.Problem(cp.Minimize(quadratic), [x >= 0]) - p.value = 3.0 - problem.solve(requires_grad=True, eps=1e-10) - # derivative() populates the delta attribute of the variables - p.delta = 1e-3 - problem.derivative() - # Because x* = 2 * p, dx*/dp = 2, so (dx*/dp)(p.delta) == 2e-3 - np.testing.assert_allclose(x.delta, 2e-3) - - Raises - ------ - ValueError - if solve was not called with ``requires_grad=True`` - SolverError - if the problem is infeasible or unbounded - """ - if s.DIFFCP not in self._solver_cache: - raise ValueError("derivative can only be called after calling " - "solve with `requires_grad=True`") - elif self.status not in s.SOLUTION_PRESENT: - raise ValueError("Differentiating through infeasible/unbounded " - "problems is not yet supported. Please file an " - "issue on Github if you need this feature.") - # TODO(akshayka): Forward differentiate dual variables as well - backward_cache = self._solver_cache[s.DIFFCP] - param_prog = self._cache.param_prog - D = backward_cache["D"] - param_deltas = {} - - gp = self._cache.gp() - if gp: - dgp2dcp = self._cache.solving_chain.get(Dgp2Dcp) - - if not self.parameters(): - for variable in self.variables(): - variable.delta = np.zeros(variable.shape) - return - - for param in self.parameters(): - delta = param.delta if param.delta is not None else np.zeros(param.shape) - if gp: - if param in dgp2dcp.canon_methods._parameters: - new_param_id = dgp2dcp.canon_methods._parameters[param].id - else: - new_param_id = param.id - param_deltas[new_param_id] = ( - 1.0/param.value * np.asarray(delta, dtype=np.float64)) - if param.id in param_prog.param_id_to_col: - # here, param generated a new parameter and also - # passed through to the param cone prog unchanged - # (because it was an exponent of a power) - param_deltas[param.id] = np.asarray(delta, - dtype=np.float64) - else: - param_deltas[param.id] = np.asarray(delta, dtype=np.float64) - dc, _, dA, db = param_prog.apply_parameters(param_deltas, - zero_offset=True) - start = time.time() - dx, _, _ = D(-dA, db, dc) - end = time.time() - backward_cache['D_TIME'] = end - start - dvars = param_prog.split_solution( - dx, [v.id for v in self.variables()]) - for variable in self.variables(): - variable.delta = dvars[variable.id] - if gp: - # x_gp = exp(x_cone_program), - # dx_gp/d x_cone_program = exp(x_cone_program) = x_gp - variable.delta *= variable.value - - def _clear_solution(self) -> None: - for v in self.variables(): - v.save_value(None) - for c in self.constraints: - for dv in c.dual_variables: - dv.save_value(None) - self._value = None - self._status = None - self._solution = None - - def unpack(self, solution) -> None: - """Updates the problem state given a Solution. - - Updates problem.status, problem.value and value of primal and dual - variables. If solution.status is in cvxpy.settins.ERROR, this method - is a no-op. - - Arguments - _________ - solution : cvxpy.Solution - A Solution object. - - Raises - ------ - ValueError - If the solution object has an invalid status - """ - if solution.status in s.SOLUTION_PRESENT: - for v in self.variables(): - v.save_value(solution.primal_vars[v.id]) - for c in self.constraints: - if c.id in solution.dual_vars: - c.save_dual_value(solution.dual_vars[c.id]) - # Eliminate confusion of problem.value versus objective.value. - self._value = self.objective.value - elif solution.status in s.INF_OR_UNB: - for v in self.variables(): - v.save_value(None) - for constr in self.constraints: - for dv in constr.dual_variables: - dv.save_value(None) - self._value = solution.opt_val - else: - raise ValueError("Cannot unpack invalid solution: %s" % solution) - - self._status = solution.status - self._solution = solution - - def unpack_results(self, solution, chain: SolvingChain, inverse_data) -> None: - """Updates the problem state given the solver results. - - Updates problem.status, problem.value and value of - primal and dual variables. - - Arguments - _________ - solution : object - The solution returned by applying the chain to the problem - and invoking the solver on the resulting data. - chain : SolvingChain - A solving chain that was used to solve the problem. - inverse_data : list - The inverse data returned by applying the chain to the problem. - - Raises - ------ - cvxpy.error.SolverError - If the solver failed - """ - - solution = chain.invert(solution, inverse_data) - if solution.status in s.INACCURATE: - warnings.warn( - "Solution may be inaccurate. Try another solver, " - "adjusting the solver settings, or solve with " - "verbose=True for more information." - ) - if solution.status == s.INFEASIBLE_OR_UNBOUNDED: - warnings.warn(INF_OR_UNB_MESSAGE) - if solution.status in s.ERROR: - raise error.SolverError( - "Solver '%s' failed. " % chain.solver.name() + - "Try another solver, or solve with verbose=True for more " - "information.") - - self.unpack(solution) - self._solver_stats = SolverStats(self._solution.attr, - chain.solver.name()) - - def __str__(self) -> str: - if len(self.constraints) == 0: - return str(self.objective) - else: - subject_to = "subject to " - lines = [str(self.objective), - subject_to + str(self.constraints[0])] - for constr in self.constraints[1:]: - lines += [len(subject_to) * " " + str(constr)] - return '\n'.join(lines) - - def __repr__(self) -> str: - return "Problem(%s, %s)" % (repr(self.objective), - repr(self.constraints)) - - def __neg__(self) -> "Problem": - return Problem(-self.objective, self.constraints) - - def __add__(self, other) -> "Problem": - if other == 0: - return self - elif not isinstance(other, Problem): - raise NotImplementedError() - return Problem(self.objective + other.objective, - unique_list(self.constraints + other.constraints)) - - def __radd__(self, other) -> "Problem": - if other == 0: - return self - else: - raise NotImplementedError() - - def __sub__(self, other) -> "Problem": - if not isinstance(other, Problem): - raise NotImplementedError() - return Problem(self.objective - other.objective, - unique_list(self.constraints + other.constraints)) - - def __rsub__(self, other) -> "Problem": - if other == 0: - return -self - else: - raise NotImplementedError() - - def __mul__(self, other) -> "Problem": - if not isinstance(other, (int, float)): - raise NotImplementedError() - return Problem(self.objective * other, self.constraints) - - __rmul__ = __mul__ - - def __div__(self, other) -> "Problem": - if not isinstance(other, (int, float)): - raise NotImplementedError() - return Problem(self.objective * (1.0 / other), self.constraints) - - def is_constant(self) -> bool: - return False - - __truediv__ = __div__ - - -class SolverStats: - """Reports some of the miscellaneous information that is returned - by the solver after solving but that is not captured directly by - the Problem instance. - - Attributes - ---------- - solver_name : str - The name of the solver. - solve_time : double - The time (in seconds) it took for the solver to solve the problem. - setup_time : double - The time (in seconds) it took for the solver to setup the problem. - num_iters : int - The number of iterations the solver had to go through to find a solution. - extra_stats : object - Extra statistics specific to the solver; these statistics are typically - returned directly from the solver, without modification by CVXPY. - This object may be a dict, or a custom Python object. - """ - def __init__(self, results_dict, solver_name: str) -> None: - self.solver_name = solver_name - self.solve_time = None - self.setup_time = None - self.num_iters = None - self.extra_stats = None - - if s.SOLVE_TIME in results_dict: - self.solve_time = results_dict[s.SOLVE_TIME] - if s.SETUP_TIME in results_dict: - self.setup_time = results_dict[s.SETUP_TIME] - if s.NUM_ITERS in results_dict: - self.num_iters = results_dict[s.NUM_ITERS] - if s.EXTRA_STATS in results_dict: - self.extra_stats = results_dict[s.EXTRA_STATS] - - -class SizeMetrics: - """Reports various metrics regarding the problem. - - Attributes - ---------- - - num_scalar_variables : integer - The number of scalar variables in the problem. - num_scalar_data : integer - The number of scalar constants and parameters in the problem. The number of - constants used across all matrices, vectors, in the problem. - Some constants are not apparent when the problem is constructed: for example, - The sum_squares expression is a wrapper for a quad_over_lin expression with a - constant 1 in the denominator. - num_scalar_eq_constr : integer - The number of scalar equality constraints in the problem. - num_scalar_leq_constr : integer - The number of scalar inequality constraints in the problem. - - max_data_dimension : integer - The longest dimension of any data block constraint or parameter. - max_big_small_squared : integer - The maximum value of (big)(small)^2 over all data blocks of the problem, where - (big) is the larger dimension and (small) is the smaller dimension - for each data block. - """ - - def __init__(self, problem: Problem) -> None: - # num_scalar_variables - self.num_scalar_variables = 0 - for var in problem.variables(): - self.num_scalar_variables += var.size - - # num_scalar_data, max_data_dimension, and max_big_small_squared - self.max_data_dimension = 0 - self.num_scalar_data = 0 - self.max_big_small_squared = 0 - for const in problem.constants()+problem.parameters(): - big = 0 - # Compute number of data - self.num_scalar_data += const.size - big = 1 if len(const.shape) == 0 else max(const.shape) - small = 1 if len(const.shape) == 0 else min(const.shape) - - # Get max data dimension: - if self.max_data_dimension < big: - self.max_data_dimension = big - - max_big_small_squared = float(big)*(float(small)**2) - if self.max_big_small_squared < max_big_small_squared: - self.max_big_small_squared = max_big_small_squared - - # num_scalar_eq_constr - self.num_scalar_eq_constr = 0 - for constraint in problem.constraints: - if isinstance(constraint, (Equality, Zero)): - self.num_scalar_eq_constr += constraint.expr.size - - # num_scalar_leq_constr - self.num_scalar_leq_constr = 0 - for constraint in problem.constraints: - if isinstance(constraint, (Inequality, NonPos, NonNeg)): - self.num_scalar_leq_constr += constraint.expr.size +""" +Copyright 2013 Steven Diamond, 2017 Akshay Agrawal + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +""" + +import time +import warnings +from collections import namedtuple +from typing import Dict, List, Optional, Union + +import numpy as np + +import cvxpy.utilities as u +import cvxpy.utilities.performance_utils as perf +from cvxpy import Constant, error +from cvxpy import settings as s +from cvxpy.atoms.atom import Atom +from cvxpy.constraints import Equality, Inequality, NonNeg, NonPos, Zero +from cvxpy.constraints.constraint import Constraint +from cvxpy.error import DPPError +from cvxpy.expressions import cvxtypes +from cvxpy.expressions.variable import Variable +from cvxpy.interface.matrix_utilities import scalar_value +from cvxpy.problems.objective import Maximize, Minimize +from cvxpy.reductions import InverseData +from cvxpy.reductions.chain import Chain +from cvxpy.reductions.dgp2dcp.dgp2dcp import Dgp2Dcp +from cvxpy.reductions.dqcp2dcp import dqcp2dcp +from cvxpy.reductions.eval_params import EvalParams +from cvxpy.reductions.flip_objective import FlipObjective +from cvxpy.reductions.solution import INF_OR_UNB_MESSAGE +from cvxpy.reductions.solvers import bisection +from cvxpy.reductions.solvers import defines as slv_def +from cvxpy.reductions.solvers.conic_solvers.conic_solver import ConicSolver +from cvxpy.reductions.solvers.defines import SOLVER_MAP_CONIC, SOLVER_MAP_QP +from cvxpy.reductions.solvers.qp_solvers.qp_solver import QpSolver +from cvxpy.reductions.solvers.solver import Solver +from cvxpy.reductions.solvers.solving_chain import (SolvingChain, + construct_solving_chain,) +from cvxpy.settings import SOLVERS +from cvxpy.utilities.deterministic import unique_list + +SolveResult = namedtuple( + 'SolveResult', + ['opt_value', 'status', 'primal_values', 'dual_values']) + + +_COL_WIDTH = 79 +_HEADER = ( + '='*_COL_WIDTH + + '\n' + + ('CVXPY').center(_COL_WIDTH) + + '\n' + + ('v' + cvxtypes.version()).center(_COL_WIDTH) + + '\n' + + '='*_COL_WIDTH +) +_COMPILATION_STR = ( + '-'*_COL_WIDTH + + '\n' + + ('Compilation').center(_COL_WIDTH) + + '\n' + + '-'*_COL_WIDTH +) +_NUM_SOLVER_STR = ( + '-'*_COL_WIDTH + + '\n' + + ('Numerical solver').center(_COL_WIDTH) + + '\n' + + '-'*_COL_WIDTH +) +_FOOTER = ( + '-'*_COL_WIDTH + + '\n' + + ('Summary').center(_COL_WIDTH) + + '\n' + + '-'*_COL_WIDTH +) + + +class Cache: + def __init__(self) -> None: + self.key = None + self.solving_chain: Optional[SolvingChain] = None + self.param_prog = None + self.inverse_data: Optional[InverseData] = None + + def invalidate(self) -> None: + self.key = None + self.solving_chain = None + self.param_prog = None + self.inverse_data = None + + def make_key(self, solver, gp, ignore_dpp): + return (solver, gp, ignore_dpp) + + def gp(self): + return self.key is not None and self.key[1] + + +def _validate_constraint(constraint): + if isinstance(constraint, Constraint): + return constraint + elif isinstance(constraint, bool): + # replace `True` or `False` values with equivalent Expressions. + return (Constant(0) <= Constant(1) if constraint else + Constant(1) <= Constant(0)) + else: + raise ValueError("Problem has an invalid constraint of type %s" % + type(constraint)) + + +class Problem(u.Canonical): + """A convex optimization problem. + + Problems are immutable, save for modification through the specification + of :class:`~cvxpy.expressions.constants.parameters.Parameter` + + Arguments + --------- + objective : Minimize or Maximize + The problem's objective. + constraints : list + The constraints on the problem variables. + """ + + # The solve methods available. + REGISTERED_SOLVE_METHODS = {} + + def __init__( + self, objective: Union[Minimize, Maximize], constraints: Optional[List[Constraint]] = None + ) -> None: + if constraints is None: + constraints = [] + # Check that objective is Minimize or Maximize. + if not isinstance(objective, (Minimize, Maximize)): + raise error.DCPError("Problem objective must be Minimize or Maximize.") + # Constraints and objective are immutable. + self._objective = objective + self._constraints = [_validate_constraint(c) for c in constraints] + self._value = None + self._status: Optional[str] = None + self._solution = None + self._cache = Cache() + self._solver_cache = {} + # Information about the shape of the problem and its constituent parts + self._size_metrics: Optional["SizeMetrics"] = None + # Benchmarks reported by the solver: + self._solver_stats: Optional["SolverStats"] = None + self._compilation_time: Optional[float] = None + self._solve_time: Optional[float] = None + self.args = [self._objective, self._constraints] + + @property + def value(self): + """float : The value from the last time the problem was solved + (or None if not solved). + """ + if self._value is None: + return None + else: + return scalar_value(self._value) + + @property + def status(self) -> str: + """str : The status from the last time the problem was solved; one + of optimal, infeasible, or unbounded (with or without + suffix inaccurate). + """ + return self._status + + @property + def solution(self): + """Solution : The solution from the last time the problem was solved. + """ + return self._solution + + @property + def objective(self) -> Union[Minimize, Maximize]: + """Minimize or Maximize : The problem's objective. + + Note that the objective cannot be reassigned after creation, + and modifying the objective after creation will result in + undefined behavior. + """ + return self._objective + + @property + def constraints(self) -> List[Constraint]: + """A shallow copy of the problem's constraints. + + Note that constraints cannot be reassigned, appended to, or otherwise + modified after creation, except through parameters. + """ + return self._constraints[:] + + @property + def param_dict(self): + """ + Expose all parameters as a dictionary + """ + return {parameters.name(): parameters for parameters in self.parameters()} + + @property + def var_dict(self) -> Dict[str, Variable]: + """ + Expose all variables as a dictionary + """ + return {variable.name(): variable for variable in self.variables()} + + @perf.compute_once + def is_dcp(self, dpp: bool = False) -> bool: + """Does the problem satisfy DCP rules? + + Arguments + --------- + dpp : bool, optional + If True, enforce the disciplined parametrized programming (DPP) + ruleset; only relevant when the problem involves Parameters. + DPP is a mild restriction of DCP. When a problem involving + Parameters is DPP, subsequent solves can be much faster than + the first one. For more information, consult the documentation at + + https://www.cvxpy.org/tutorial/advanced/index.html#disciplined-parametrized-programming + + Returns + ------- + bool + True if the Expression is DCP, False otherwise. + """ + return all( + expr.is_dcp(dpp) for expr in self.constraints + [self.objective]) + + @perf.compute_once + def is_dgp(self, dpp: bool = False) -> bool: + """Does the problem satisfy DGP rules? + + Arguments + --------- + dpp : bool, optional + If True, enforce the disciplined parametrized programming (DPP) + ruleset; only relevant when the problem involves Parameters. + DPP is a mild restriction of DGP. When a problem involving + Parameters is DPP, subsequent solves can be much faster than + the first one. For more information, consult the documentation at + + https://www.cvxpy.org/tutorial/advanced/index.html#disciplined-parametrized-programming + + Returns + ------- + bool + True if the Expression is DGP, False otherwise. + """ + return all( + expr.is_dgp(dpp) for expr in self.constraints + [self.objective]) + + @perf.compute_once + def is_dqcp(self) -> bool: + """Does the problem satisfy the DQCP rules? + """ + return all( + expr.is_dqcp() for expr in self.constraints + [self.objective]) + + @perf.compute_once + def is_dpp(self, context: str = 'dcp') -> bool: + """Does the problem satisfy DPP rules? + + DPP is a mild restriction of DGP. When a problem involving + Parameters is DPP, subsequent solves can be much faster than + the first one. For more information, consult the documentation at + + https://www.cvxpy.org/tutorial/advanced/index.html#disciplined-parametrized-programming + + Arguments + --------- + context : str + Whether to check DPP-compliance for DCP or DGP; ``context`` should + be either ``'dcp'`` or ``'dgp'``. Calling ``problem.is_dpp('dcp')`` + is equivalent to ``problem.is_dcp(dpp=True)``, and + `problem.is_dpp('dgp')`` is equivalent to + `problem.is_dgp(dpp=True)`. + + Returns + ------- + bool + Whether the problem satisfies the DPP rules. + """ + if context.lower() == 'dcp': + return self.is_dcp(dpp=True) + elif context.lower() == 'dgp': + return self.is_dgp(dpp=True) + else: + raise ValueError("Unsupported context ", context) + + @perf.compute_once + def is_qp(self) -> bool: + """Is problem a quadratic program? + """ + for c in self.constraints: + if not (isinstance(c, (Equality, Zero)) or c.args[0].is_pwl()): + return False + for var in self.variables(): + if var.is_psd() or var.is_nsd(): + return False + return (self.is_dcp() and self.objective.args[0].is_qpwa()) + + @perf.compute_once + def is_mixed_integer(self) -> bool: + return any(v.attributes['boolean'] or v.attributes['integer'] + for v in self.variables()) + + @perf.compute_once + def variables(self) -> List[Variable]: + """Accessor method for variables. + + Returns + ------- + list of :class:`~cvxpy.expressions.variable.Variable` + A list of the variables in the problem. + """ + vars_ = self.objective.variables() + for constr in self.constraints: + vars_ += constr.variables() + return unique_list(vars_) + + @perf.compute_once + def parameters(self): + """Accessor method for parameters. + + Returns + ------- + list of :class:`~cvxpy.expressions.constants.parameter.Parameter` + A list of the parameters in the problem. + """ + params = self.objective.parameters() + for constr in self.constraints: + params += constr.parameters() + return unique_list(params) + + @perf.compute_once + def constants(self) -> List[Constant]: + """Accessor method for constants. + + Returns + ------- + list of :class:`~cvxpy.expressions.constants.constant.Constant` + A list of the constants in the problem. + """ + const_dict = {} + constants_ = self.objective.constants() + for constr in self.constraints: + constants_ += constr.constants() + # Note that numpy matrices are not hashable, so we use the built-in + # function "id" + const_dict = {id(constant): constant for constant in constants_} + return list(const_dict.values()) + + def atoms(self) -> List[Atom]: + """Accessor method for atoms. + + Returns + ------- + list of :class:`~cvxpy.atoms.Atom` + A list of the atom types in the problem; note that this list + contains classes, not instances. + """ + atoms = self.objective.atoms() + for constr in self.constraints: + atoms += constr.atoms() + return unique_list(atoms) + + @property + def size_metrics(self) -> "SizeMetrics": + """:class:`~cvxpy.problems.problem.SizeMetrics` : Information about the problem's size. + """ + if self._size_metrics is None: + self._size_metrics = SizeMetrics(self) + return self._size_metrics + + @property + def solver_stats(self) -> "SolverStats": + """:class:`~cvxpy.problems.problem.SolverStats` : Information returned by the solver. + """ + return self._solver_stats + + def solve(self, *args, **kwargs): + """Compiles and solves the problem using the specified method. + + Populates the :code:`status` and :code:`value` attributes on the + problem object as a side-effect. + + Arguments + --------- + solver : str, optional + The solver to use. For example, 'ECOS', 'SCS', or 'OSQP'. + verbose : bool, optional + Overrides the default of hiding solver output, and prints + logging information describing CVXPY's compilation process. + gp : bool, optional + If True, parses the problem as a disciplined geometric program + instead of a disciplined convex program. + qcp : bool, optional + If True, parses the problem as a disciplined quasiconvex program + instead of a disciplined convex program. + requires_grad : bool, optional + Makes it possible to compute gradients of a solution with respect to + Parameters by calling ``problem.backward()`` after solving, or to + compute perturbations to the variables given perturbations to Parameters by + calling ``problem.derivative()``. + + Gradients are only supported for DCP and DGP problems, not + quasiconvex problems. When computing gradients (i.e., when + this argument is True), the problem must satisfy the DPP rules. + enforce_dpp : bool, optional + When True, a DPPError will be thrown when trying to solve a non-DPP + problem (instead of just a warning). Only relevant for problems + involving Parameters. Defaults to False. + ignore_dpp : bool, optional + When True, DPP problems will be treated as non-DPP, + which may speed up compilation. Defaults to False. + method : function, optional + A custom solve method to use. + kwargs : keywords, optional + Additional solver specific arguments. See Notes below. + + Notes + ------ + CVXPY interfaces with a wide range of solvers; the algorithms used by these solvers + have arguments relating to stopping criteria, and strategies to improve solution quality. + + There is no one choice of arguments which is perfect for every problem. If you are not + getting satisfactory results from a solver, you can try changing its arguments. The + exact way this is done depends on the specific solver. Here are some examples: + + :: + + prob.solve(solver='ECOS', abstol=1e-6) + prob.solve(solver='OSQP', max_iter=10000). + mydict = {"MSK_DPAR_INTPNT_CO_TOL_NEAR_REL": 10} + prob.solve(solver='MOSEK', mosek_params=mydict). + + You should refer to CVXPY's web documentation for details on how to pass solver + solver arguments, available at + + https://www.cvxpy.org/tutorial/advanced/index.html#setting-solver-options + + Returns + ------- + float + The optimal value for the problem, or a string indicating + why the problem could not be solved. + + Raises + ------ + cvxpy.error.DCPError + Raised if the problem is not DCP and `gp` is False. + cvxpy.error.DGPError + Raised if the problem is not DGP and `gp` is True. + cvxpy.error.DPPError + Raised if DPP settings are invalid. + cvxpy.error.SolverError + Raised if no suitable solver exists among the installed solvers, + or if an unanticipated error is encountered. + """ + func_name = kwargs.pop("method", None) + if func_name is not None: + solve_func = Problem.REGISTERED_SOLVE_METHODS[func_name] + else: + solve_func = Problem._solve + return solve_func(self, *args, **kwargs) + + @classmethod + def register_solve(cls, name: str, func) -> None: + """Adds a solve method to the Problem class. + + Arguments + --------- + name : str + The keyword for the method. + func : function + The function that executes the solve method. This function must + take as its first argument the problem instance to solve. + """ + cls.REGISTERED_SOLVE_METHODS[name] = func + + def get_problem_data( + self, solver, + gp: bool = False, + enforce_dpp: bool = False, + ignore_dpp: bool = False, + verbose: bool = False + ): + """Returns the problem data used in the call to the solver. + + When a problem is solved, CVXPY creates a chain of reductions enclosed + in a :class:`~cvxpy.reductions.solvers.solving_chain.SolvingChain`, + and compiles it to some low-level representation that is + compatible with the targeted solver. This method returns that low-level + representation. + + For some solving chains, this low-level representation is a dictionary + that contains exactly those arguments that were supplied to the solver; + however, for other solving chains, the data is an intermediate + representation that is compiled even further by the solver interfaces. + + A solution to the equivalent low-level problem can be obtained via the + data by invoking the `solve_via_data` method of the returned solving + chain, a thin wrapper around the code external to CVXPY that further + processes and solves the problem. Invoke the unpack_results method + to recover a solution to the original problem. + + For example: + + :: + + objective = ... + constraints = ... + problem = cp.Problem(objective, constraints) + data, chain, inverse_data = problem.get_problem_data(cp.SCS) + # calls SCS using `data` + soln = chain.solve_via_data(problem, data) + # unpacks the solution returned by SCS into `problem` + problem.unpack_results(soln, chain, inverse_data) + + Alternatively, the `data` dictionary returned by this method + contains enough information to bypass CVXPY and call the solver + directly. + + For example: + + :: + + problem = cp.Problem(objective, constraints) + data, _, _ = problem.get_problem_data(cp.SCS) + + import scs + probdata = { + 'A': data['A'], + 'b': data['b'], + 'c': data['c'], + } + cone_dims = data['dims'] + cones = { + "f": cone_dims.zero, + "l": cone_dims.nonneg, + "q": cone_dims.soc, + "ep": cone_dims.exp, + "s": cone_dims.psd, + } + soln = scs.solve(data, cones) + + The structure of the data dict that CVXPY returns depends on the + solver. For details, consult the solver interfaces in + `cvxpy/reductions/solvers`. + + Arguments + --------- + solver : str + The solver the problem data is for. + gp : bool, optional + If True, then parses the problem as a disciplined geometric program + instead of a disciplined convex program. + enforce_dpp : bool, optional + When True, a DPPError will be thrown when trying to parse a non-DPP + problem (instead of just a warning). Defaults to False. + ignore_dpp : bool, optional + When True, DPP problems will be treated as non-DPP, + which may speed up compilation. Defaults to False. + verbose : bool, optional + If True, print verbose output related to problem compilation. + + Returns + ------- + dict or object + lowest level representation of problem + SolvingChain + The solving chain that created the data. + list + The inverse data generated by the chain. + + Raises + ------ + cvxpy.error.DPPError + Raised if DPP settings are invalid. + """ + # Invalid DPP setting. + # Must be checked here to avoid cache issues. + if enforce_dpp and ignore_dpp: + raise DPPError("Cannot set enforce_dpp = True and ignore_dpp = True.") + + start = time.time() + # Cache includes ignore_dpp because it alters compilation. + key = self._cache.make_key(solver, gp, ignore_dpp) + if key != self._cache.key: + self._cache.invalidate() + solving_chain = self._construct_chain( + solver=solver, gp=gp, + enforce_dpp=enforce_dpp, + ignore_dpp=ignore_dpp) + self._cache.key = key + self._cache.solving_chain = solving_chain + self._solver_cache = {} + else: + solving_chain = self._cache.solving_chain + + if verbose: + print(_COMPILATION_STR) + + if self._cache.param_prog is not None: + # fast path, bypasses application of reductions + if verbose: + s.LOGGER.info( + 'Using cached ASA map, for faster compilation ' + '(bypassing reduction chain).') + if gp: + dgp2dcp = self._cache.solving_chain.get(Dgp2Dcp) + # Parameters in the param cone prog are the logs + # of parameters in the original problem (with one exception: + # parameters appearing as exponents (in power and gmatmul + # atoms) are unchanged. + old_params_to_new_params = dgp2dcp.canon_methods._parameters + for param in self.parameters(): + + if param in old_params_to_new_params: + old_params_to_new_params[param].value = np.log( + param.value) + + data, solver_inverse_data = solving_chain.solver.apply( + self._cache.param_prog) + inverse_data = self._cache.inverse_data + [solver_inverse_data] + self._compilation_time = time.time() - start + if verbose: + s.LOGGER.info( + 'Finished problem compilation ' + '(took %.3e seconds).', self._compilation_time) + else: + if verbose: + solver_name = solving_chain.reductions[-1].name() + reduction_chain_str = ' -> '.join( + type(r).__name__ for r in solving_chain.reductions) + s.LOGGER.info( + 'Compiling problem (target solver=%s).', solver_name) + s.LOGGER.info('Reduction chain: %s', reduction_chain_str) + data, inverse_data = solving_chain.apply(self, verbose) + safe_to_cache = ( + isinstance(data, dict) + and s.PARAM_PROB in data + and not any(isinstance(reduction, EvalParams) + for reduction in solving_chain.reductions) + ) + self._compilation_time = time.time() - start + if verbose: + s.LOGGER.info( + 'Finished problem compilation ' + '(took %.3e seconds).', self._compilation_time) + if safe_to_cache: + if verbose and self.parameters(): + s.LOGGER.info( + '(Subsequent compilations of this problem, using the ' + 'same arguments, should ' 'take less time.)') + self._cache.param_prog = data[s.PARAM_PROB] + # the last datum in inverse_data corresponds to the solver, + # so we shouldn't cache it + self._cache.inverse_data = inverse_data[:-1] + return data, solving_chain, inverse_data + + def _find_candidate_solvers(self, + solver=None, + gp: bool = False): + """ + Find candidate solvers for the current problem. If solver + is not None, it checks if the specified solver is compatible + with the problem passed. + + Arguments + --------- + solver : Union[string, Solver, None] + The name of the solver with which to solve the problem or an + instance of a custom solver. If no solver is supplied + (i.e., if solver is None), then the targeted solver may be any + of those that are installed. If the problem is variable-free, + then this parameter is ignored. + gp : bool + If True, the problem is parsed as a Disciplined Geometric Program + instead of as a Disciplined Convex Program. + + Returns + ------- + dict + A dictionary of compatible solvers divided in `qp_solvers` + and `conic_solvers`. + + Raises + ------ + cvxpy.error.SolverError + Raised if the problem is not DCP and `gp` is False. + cvxpy.error.DGPError + Raised if the problem is not DGP and `gp` is True. + """ + candidates = {'qp_solvers': [], + 'conic_solvers': []} + if isinstance(solver, Solver): + return self._add_custom_solver_candidates(solver) + + if solver is not None: + if solver not in slv_def.INSTALLED_SOLVERS: + raise error.SolverError("The solver %s is not installed." % solver) + if solver in slv_def.CONIC_SOLVERS: + candidates['conic_solvers'] += [solver] + if solver in slv_def.QP_SOLVERS: + candidates['qp_solvers'] += [solver] + else: + candidates['qp_solvers'] = [s for s in slv_def.INSTALLED_SOLVERS + if s in slv_def.QP_SOLVERS] + candidates['conic_solvers'] = [] + # ECOS_BB can only be called explicitly. + for slv in slv_def.INSTALLED_SOLVERS: + if slv in slv_def.CONIC_SOLVERS and slv != s.ECOS_BB: + candidates['conic_solvers'].append(slv) + + # If gp we must have only conic solvers + if gp: + if solver is not None and solver not in slv_def.CONIC_SOLVERS: + raise error.SolverError( + "When `gp=True`, `solver` must be a conic solver " + "(received '%s'); try calling " % solver + + " `solve()` with `solver=cvxpy.ECOS`." + ) + elif solver is None: + candidates['qp_solvers'] = [] # No QP solvers allowed + + if self.is_mixed_integer(): + # ECOS_BB must be called explicitly. + if slv_def.INSTALLED_MI_SOLVERS == [s.ECOS_BB] and solver != s.ECOS_BB: + msg = """ + + You need a mixed-integer solver for this model. Refer to the documentation + https://www.cvxpy.org/tutorial/advanced/index.html#mixed-integer-programs + for discussion on this topic. + + Quick fix 1: if you install the python package CVXOPT (pip install cvxopt), + then CVXPY can use the open-source mixed-integer linear programming + solver `GLPK`. If your problem is nonlinear then you can install SCIP + (pip install pyscipopt). + + Quick fix 2: you can explicitly specify solver='ECOS_BB'. This may result + in incorrect solutions and is not recommended. + """ + raise error.SolverError(msg) + # TODO: provide a useful error message when the problem is nonlinear but + # the only installed mixed-integer solvers are MILP solvers (e.g., GLPK_MI). + candidates['qp_solvers'] = [ + s for s in candidates['qp_solvers'] + if slv_def.SOLVER_MAP_QP[s].MIP_CAPABLE] + candidates['conic_solvers'] = [ + s for s in candidates['conic_solvers'] + if slv_def.SOLVER_MAP_CONIC[s].MIP_CAPABLE] + if not candidates['conic_solvers'] and \ + not candidates['qp_solvers']: + raise error.SolverError( + "Problem is mixed-integer, but candidate " + "QP/Conic solvers (%s) are not MIP-capable." % + (candidates['qp_solvers'] + + candidates['conic_solvers'])) + + return candidates + + def _add_custom_solver_candidates(self, custom_solver: Solver): + """ + Returns a list of candidate solvers where custom_solver is the only potential option. + + Arguments + --------- + custom_solver : Solver + + Returns + ------- + dict + A dictionary of compatible solvers divided in `qp_solvers` + and `conic_solvers`. + + Raises + ------ + cvxpy.error.SolverError + Raised if the name of the custom solver conflicts with the name of some officially + supported solver + """ + if custom_solver.name() in SOLVERS: + message = "Custom solvers must have a different name than the officially supported ones" + raise(error.SolverError(message)) + + candidates = {'qp_solvers': [], 'conic_solvers': []} + if not self.is_mixed_integer() or custom_solver.MIP_CAPABLE: + if isinstance(custom_solver, QpSolver): + SOLVER_MAP_QP[custom_solver.name()] = custom_solver + candidates['qp_solvers'] = [custom_solver.name()] + elif isinstance(custom_solver, ConicSolver): + SOLVER_MAP_CONIC[custom_solver.name()] = custom_solver + candidates['conic_solvers'] = [custom_solver.name()] + return candidates + + def _construct_chain( + self, solver: Optional[str] = None, gp: bool = False, + enforce_dpp: bool = False, ignore_dpp: bool = False + ) -> SolvingChain: + """ + Construct the chains required to reformulate and solve the problem. + + In particular, this function + + # finds the candidate solvers + # constructs the solving chain that performs the + numeric reductions and solves the problem. + + Arguments + --------- + solver : str, optional + The solver to use. Defaults to ECOS. + gp : bool, optional + If True, the problem is parsed as a Disciplined Geometric Program + instead of as a Disciplined Convex Program. + enforce_dpp : bool, optional + Whether to error on DPP violations. + ignore_dpp : bool, optional + When True, DPP problems will be treated as non-DPP, + which may speed up compilation. Defaults to False. + + Returns + ------- + A solving chain + """ + candidate_solvers = self._find_candidate_solvers(solver=solver, gp=gp) + self._sort_candidate_solvers(candidate_solvers) + return construct_solving_chain(self, candidate_solvers, gp=gp, + enforce_dpp=enforce_dpp, + ignore_dpp=ignore_dpp) + + @staticmethod + def _sort_candidate_solvers(solvers) -> None: + """Sorts candidate solvers lists according to slv_def.CONIC_SOLVERS/QP_SOLVERS + + Arguments + --------- + candidates : dict + Dictionary of candidate solvers divided in qp_solvers + and conic_solvers + Returns + ------- + None + """ + if len(solvers['conic_solvers']) > 1: + solvers['conic_solvers'] = sorted( + solvers['conic_solvers'], key=lambda s: slv_def.CONIC_SOLVERS.index(s) + ) + if len(solvers['qp_solvers']) > 1: + solvers['qp_solvers'] = sorted( + solvers['qp_solvers'], key=lambda s: slv_def.QP_SOLVERS.index(s) + ) + + def _invalidate_cache(self) -> None: + self._cache_key = None + self._solving_chain = None + self._param_prog = None + self._inverse_data = None + + def _solve(self, + solver: str = None, + warm_start: bool = True, + verbose: bool = False, + gp: bool = False, + qcp: bool = False, + requires_grad: bool = False, + enforce_dpp: bool = False, + ignore_dpp: bool = False, + **kwargs): + """Solves a DCP compliant optimization problem. + + Saves the values of primal and dual variables in the variable + and constraint objects, respectively. + + Arguments + --------- + solver : str, optional + The solver to use. Defaults to ECOS. + warm_start : bool, optional + Should the previous solver result be used to warm start? + verbose : bool, optional + Overrides the default of hiding solver output. + gp : bool, optional + If True, parses the problem as a disciplined geometric program. + qcp : bool, optional + If True, parses the problem as a disciplined quasiconvex program. + requires_grad : bool, optional + Makes it possible to compute gradients with respect to + parameters by calling `backward()` after solving, or to compute + perturbations to the variables by calling `derivative()`. When + True, the solver must be SCS, and dqcp must be False. + A DPPError is thrown when problem is not DPP. + enforce_dpp : bool, optional + When True, a DPPError will be thrown when trying to solve a non-DPP + problem (instead of just a warning). Defaults to False. + ignore_dpp : bool, optional + When True, DPP problems will be treated as non-DPP, + which may speed up compilation. Defaults to False. + kwargs : dict, optional + A dict of options that will be passed to the specific solver. + In general, these options will override any default settings + imposed by cvxpy. + + Returns + ------- + float + The optimal value for the problem, or a string indicating + why the problem could not be solved. + """ + if verbose: + print(_HEADER) + + for parameter in self.parameters(): + if parameter.value is None: + raise error.ParameterError( + "A Parameter (whose name is '%s') does not have a value " + "associated with it; all Parameter objects must have " + "values before solving a problem." % parameter.name()) + + if verbose: + n_variables = sum(np.prod(v.shape) for v in self.variables()) + n_parameters = sum(np.prod(p.shape) for p in self.parameters()) + s.LOGGER.info( + 'Your problem has %d variables, ' + '%d constraints, and ' '%d parameters.', + n_variables, len(self.constraints), n_parameters) + curvatures = [] + if self.is_dcp(): + curvatures.append('DCP') + if self.is_dgp(): + curvatures.append('DGP') + if self.is_dqcp(): + curvatures.append('DQCP') + s.LOGGER.info( + 'It is compliant with the following grammars: %s', + ', '.join(curvatures)) + if n_parameters == 0: + s.LOGGER.info( + '(If you need to solve this problem multiple times, ' + 'but with different data, consider using parameters.)') + s.LOGGER.info( + 'CVXPY will first compile your problem; then, it will ' + 'invoke a numerical solver to obtain a solution.') + + if requires_grad: + dpp_context = 'dgp' if gp else 'dcp' + if qcp: + raise ValueError("Cannot compute gradients of DQCP problems.") + elif not self.is_dpp(dpp_context): + raise error.DPPError("Problem is not DPP (when requires_grad " + "is True, problem must be DPP).") + elif solver is not None and solver not in [s.SCS, s.DIFFCP]: + raise ValueError("When requires_grad is True, the only " + "supported solver is SCS " + "(received %s)." % solver) + elif s.DIFFCP not in slv_def.INSTALLED_SOLVERS: + raise ModuleNotFoundError( + "The Python package diffcp must be installed to " + "differentiate through problems. Please follow the " + "installation instructions at " + "https://github.com/cvxgrp/diffcp") + else: + solver = s.DIFFCP + else: + if gp and qcp: + raise ValueError("At most one of `gp` and `qcp` can be True.") + if qcp and not self.is_dcp(): + if not self.is_dqcp(): + raise error.DQCPError("The problem is not DQCP.") + if verbose: + s.LOGGER.info( + 'Reducing DQCP problem to a one-parameter ' + 'family of DCP problems, for bisection.') + reductions = [dqcp2dcp.Dqcp2Dcp()] + start = time.time() + if type(self.objective) == Maximize: + reductions = [FlipObjective()] + reductions + # flip objective flips the sign of the objective + low, high = kwargs.get("low"), kwargs.get("high") + if high is not None: + kwargs["low"] = high * -1 + if low is not None: + kwargs["high"] = low * -1 + chain = Chain(problem=self, reductions=reductions) + soln = bisection.bisect( + chain.reduce(), solver=solver, verbose=verbose, **kwargs) + self.unpack(chain.retrieve(soln)) + return self.value + + data, solving_chain, inverse_data = self.get_problem_data( + solver, gp, enforce_dpp, ignore_dpp, verbose) + + if verbose: + print(_NUM_SOLVER_STR) + s.LOGGER.info( + 'Invoking solver %s to obtain a solution.', + solving_chain.reductions[-1].name()) + start = time.time() + solution = solving_chain.solve_via_data( + self, data, warm_start, verbose, kwargs) + end = time.time() + self._solve_time = end - start + self.unpack_results(solution, solving_chain, inverse_data) + if verbose: + print(_FOOTER) + s.LOGGER.info('Problem status: %s', self.status) + val = self.value if self.value is not None else np.NaN + s.LOGGER.info('Optimal value: %.3e', val) + s.LOGGER.info('Compilation took %.3e seconds', self._compilation_time) + s.LOGGER.info( + 'Solver (including time spent in interface) took ' + '%.3e seconds', self._solve_time) + return self.value + + def backward(self) -> None: + """Compute the gradient of a solution with respect to Parameters. + + This method differentiates through the solution map of the problem, + obtaining the gradient of a solution with respect to the Parameters. + In other words, it calculates the sensitivities of the Parameters + with respect to perturbations in the optimal Variable values. This + can be useful for integrating CVXPY into automatic differentiation + toolkits. + + ``backward()`` populates the ``gradient`` attribute of each Parameter + in the problem as a side-effect. It can only be called after calling + ``solve()`` with ``requires_grad=True``. + + Below is a simple example: + + :: + + import cvxpy as cp + import numpy as np + + p = cp.Parameter() + x = cp.Variable() + quadratic = cp.square(x - 2 * p) + problem = cp.Problem(cp.Minimize(quadratic), [x >= 0]) + p.value = 3.0 + problem.solve(requires_grad=True, eps=1e-10) + # backward() populates the gradient attribute of the parameters + problem.backward() + # Because x* = 2 * p, dx*/dp = 2 + np.testing.assert_allclose(p.gradient, 2.0) + + In the above example, the gradient could easily be computed by hand. + The ``backward()`` is useful because for almost all problems, the + gradient cannot be computed analytically. + + This method can be used to differentiate through any DCP or DGP + problem, as long as the problem is DPP compliant (i.e., + ``problem.is_dcp(dpp=True)`` or ``problem.is_dgp(dpp=True)`` evaluates to + ``True``). + + This method uses the chain rule to evaluate the gradients of a + scalar-valued function of the Variables with respect to the Parameters. + For example, let x be a variable and p a Parameter; x and p might be + scalars, vectors, or matrices. Let f be a scalar-valued function, with + z = f(x). Then this method computes dz/dp = (dz/dx) (dx/p). dz/dx + is chosen as the all-ones vector by default, corresponding to + choosing f to be the sum function. You can specify a custom value for + dz/dx by setting the ``gradient`` attribute on your variables. For example, + + :: + + import cvxpy as cp + import numpy as np + + + b = cp.Parameter() + x = cp.Variable() + quadratic = cp.square(x - 2 * b) + problem = cp.Problem(cp.Minimize(quadratic), [x >= 0]) + b.value = 3. + problem.solve(requires_grad=True, eps=1e-10) + x.gradient = 4. + problem.backward() + # dz/dp = dz/dx dx/dp = 4. * 2. == 8. + np.testing.assert_allclose(b.gradient, 8.) + + The ``gradient`` attribute on a variable can also be interpreted as a + perturbation to its optimal value. + + Raises + ------ + ValueError + if solve was not called with ``requires_grad=True`` + SolverError + if the problem is infeasible or unbounded + """ + if s.DIFFCP not in self._solver_cache: + raise ValueError("backward can only be called after calling " + "solve with `requires_grad=True`") + elif self.status not in s.SOLUTION_PRESENT: + raise error.SolverError("Backpropagating through " + "infeasible/unbounded problems is not " + "yet supported. Please file an issue on " + "Github if you need this feature.") + + # TODO(akshayka): Backpropagate through dual variables as well. + backward_cache = self._solver_cache[s.DIFFCP] + DT = backward_cache["DT"] + zeros = np.zeros(backward_cache["s"].shape) + del_vars = {} + + gp = self._cache.gp() + for variable in self.variables(): + if variable.gradient is None: + del_vars[variable.id] = np.ones(variable.shape) + else: + del_vars[variable.id] = np.asarray(variable.gradient, + dtype=np.float64) + if gp: + # x_gp = exp(x_cone_program), + # dx_gp/d x_cone_program = exp(x_cone_program) = x_gp + del_vars[variable.id] *= variable.value + + dx = self._cache.param_prog.split_adjoint(del_vars) + start = time.time() + dA, db, dc = DT(dx, zeros, zeros) + end = time.time() + backward_cache['DT_TIME'] = end - start + dparams = self._cache.param_prog.apply_param_jac(dc, -dA, db) + + if not gp: + for param in self.parameters(): + param.gradient = dparams[param.id] + else: + dgp2dcp = self._cache.solving_chain.get(Dgp2Dcp) + old_params_to_new_params = dgp2dcp.canon_methods._parameters + for param in self.parameters(): + # Note: if param is an exponent in a power or gmatmul atom, + # then the parameter passes through unchanged to the DCP + # program; if the param is also used elsewhere (not as an + # exponent), then param will also be in + # old_params_to_new_params. Therefore, param.gradient = + # dparams[param.id] (or 0) + 1/param*dparams[new_param.id] + # + # Note that param.id is in dparams if and only if + # param was used as an exponent (because this means that + # the parameter entered the DCP problem unchanged.) + grad = 0.0 if param.id not in dparams else dparams[param.id] + if param in old_params_to_new_params: + new_param = old_params_to_new_params[param] + # new_param.value == log(param), apply chain rule + grad += (1.0 / param.value) * dparams[new_param.id] + param.gradient = grad + + def derivative(self) -> None: + """Apply the derivative of the solution map to perturbations in the Parameters + + This method applies the derivative of the solution map to perturbations + in the Parameters to obtain perturbations in the optimal values of the + Variables. In other words, it tells you how the optimal values of the + Variables would be changed by small changes to the Parameters. + + You can specify perturbations in a Parameter by setting its ``delta`` + attribute (if unspecified, the perturbation defaults to 0). + + This method populates the ``delta`` attribute of the Variables as a + side-effect. + + This method can only be called after calling ``solve()`` with + ``requires_grad=True``. It is compatible with both DCP and DGP + problems (that are also DPP-compliant). + + Below is a simple example: + + :: + + import cvxpy as cp + import numpy as np + + p = cp.Parameter() + x = cp.Variable() + quadratic = cp.square(x - 2 * p) + problem = cp.Problem(cp.Minimize(quadratic), [x >= 0]) + p.value = 3.0 + problem.solve(requires_grad=True, eps=1e-10) + # derivative() populates the delta attribute of the variables + p.delta = 1e-3 + problem.derivative() + # Because x* = 2 * p, dx*/dp = 2, so (dx*/dp)(p.delta) == 2e-3 + np.testing.assert_allclose(x.delta, 2e-3) + + Raises + ------ + ValueError + if solve was not called with ``requires_grad=True`` + SolverError + if the problem is infeasible or unbounded + """ + if s.DIFFCP not in self._solver_cache: + raise ValueError("derivative can only be called after calling " + "solve with `requires_grad=True`") + elif self.status not in s.SOLUTION_PRESENT: + raise ValueError("Differentiating through infeasible/unbounded " + "problems is not yet supported. Please file an " + "issue on Github if you need this feature.") + # TODO(akshayka): Forward differentiate dual variables as well + backward_cache = self._solver_cache[s.DIFFCP] + param_prog = self._cache.param_prog + D = backward_cache["D"] + param_deltas = {} + + gp = self._cache.gp() + if gp: + dgp2dcp = self._cache.solving_chain.get(Dgp2Dcp) + + if not self.parameters(): + for variable in self.variables(): + variable.delta = np.zeros(variable.shape) + return + + for param in self.parameters(): + delta = param.delta if param.delta is not None else np.zeros(param.shape) + if gp: + if param in dgp2dcp.canon_methods._parameters: + new_param_id = dgp2dcp.canon_methods._parameters[param].id + else: + new_param_id = param.id + param_deltas[new_param_id] = ( + 1.0/param.value * np.asarray(delta, dtype=np.float64)) + if param.id in param_prog.param_id_to_col: + # here, param generated a new parameter and also + # passed through to the param cone prog unchanged + # (because it was an exponent of a power) + param_deltas[param.id] = np.asarray(delta, + dtype=np.float64) + else: + param_deltas[param.id] = np.asarray(delta, dtype=np.float64) + dc, _, dA, db = param_prog.apply_parameters(param_deltas, + zero_offset=True) + start = time.time() + dx, _, _ = D(-dA, db, dc) + end = time.time() + backward_cache['D_TIME'] = end - start + dvars = param_prog.split_solution( + dx, [v.id for v in self.variables()]) + for variable in self.variables(): + variable.delta = dvars[variable.id] + if gp: + # x_gp = exp(x_cone_program), + # dx_gp/d x_cone_program = exp(x_cone_program) = x_gp + variable.delta *= variable.value + + def _clear_solution(self) -> None: + for v in self.variables(): + v.save_value(None) + for c in self.constraints: + for dv in c.dual_variables: + dv.save_value(None) + self._value = None + self._status = None + self._solution = None + + def unpack(self, solution) -> None: + """Updates the problem state given a Solution. + + Updates problem.status, problem.value and value of primal and dual + variables. If solution.status is in cvxpy.settins.ERROR, this method + is a no-op. + + Arguments + _________ + solution : cvxpy.Solution + A Solution object. + + Raises + ------ + ValueError + If the solution object has an invalid status + """ + if solution.status in s.SOLUTION_PRESENT: + for v in self.variables(): + v.save_value(solution.primal_vars[v.id]) + for c in self.constraints: + if c.id in solution.dual_vars: + c.save_dual_value(solution.dual_vars[c.id]) + # Eliminate confusion of problem.value versus objective.value. + self._value = self.objective.value + elif solution.status in s.INF_OR_UNB: + for v in self.variables(): + v.save_value(None) + for constr in self.constraints: + for dv in constr.dual_variables: + dv.save_value(None) + self._value = solution.opt_val + else: + raise ValueError("Cannot unpack invalid solution: %s" % solution) + + self._status = solution.status + self._solution = solution + + def unpack_results(self, solution, chain: SolvingChain, inverse_data) -> None: + """Updates the problem state given the solver results. + + Updates problem.status, problem.value and value of + primal and dual variables. + + Arguments + _________ + solution : object + The solution returned by applying the chain to the problem + and invoking the solver on the resulting data. + chain : SolvingChain + A solving chain that was used to solve the problem. + inverse_data : list + The inverse data returned by applying the chain to the problem. + + Raises + ------ + cvxpy.error.SolverError + If the solver failed + """ + + solution = chain.invert(solution, inverse_data) + if solution.status in s.INACCURATE: + warnings.warn( + "Solution may be inaccurate. Try another solver, " + "adjusting the solver settings, or solve with " + "verbose=True for more information." + ) + if solution.status == s.INFEASIBLE_OR_UNBOUNDED: + warnings.warn(INF_OR_UNB_MESSAGE) + if solution.status in s.ERROR: + raise error.SolverError( + "Solver '%s' failed. " % chain.solver.name() + + "Try another solver, or solve with verbose=True for more " + "information.") + + self.unpack(solution) + self._solver_stats = SolverStats(self._solution.attr, + chain.solver.name()) + + def __str__(self) -> str: + if len(self.constraints) == 0: + return str(self.objective) + else: + subject_to = "subject to " + lines = [str(self.objective), + subject_to + str(self.constraints[0])] + for constr in self.constraints[1:]: + lines += [len(subject_to) * " " + str(constr)] + return '\n'.join(lines) + + def __repr__(self) -> str: + return "Problem(%s, %s)" % (repr(self.objective), + repr(self.constraints)) + + def __neg__(self) -> "Problem": + return Problem(-self.objective, self.constraints) + + def __add__(self, other) -> "Problem": + if other == 0: + return self + elif not isinstance(other, Problem): + raise NotImplementedError() + return Problem(self.objective + other.objective, + unique_list(self.constraints + other.constraints)) + + def __radd__(self, other) -> "Problem": + if other == 0: + return self + else: + raise NotImplementedError() + + def __sub__(self, other) -> "Problem": + if not isinstance(other, Problem): + raise NotImplementedError() + return Problem(self.objective - other.objective, + unique_list(self.constraints + other.constraints)) + + def __rsub__(self, other) -> "Problem": + if other == 0: + return -self + else: + raise NotImplementedError() + + def __mul__(self, other) -> "Problem": + if not isinstance(other, (int, float)): + raise NotImplementedError() + return Problem(self.objective * other, self.constraints) + + __rmul__ = __mul__ + + def __div__(self, other) -> "Problem": + if not isinstance(other, (int, float)): + raise NotImplementedError() + return Problem(self.objective * (1.0 / other), self.constraints) + + def is_constant(self) -> bool: + return False + + __truediv__ = __div__ + + +class SolverStats: + """Reports some of the miscellaneous information that is returned + by the solver after solving but that is not captured directly by + the Problem instance. + + Attributes + ---------- + solver_name : str + The name of the solver. + solve_time : double + The time (in seconds) it took for the solver to solve the problem. + setup_time : double + The time (in seconds) it took for the solver to setup the problem. + num_iters : int + The number of iterations the solver had to go through to find a solution. + extra_stats : object + Extra statistics specific to the solver; these statistics are typically + returned directly from the solver, without modification by CVXPY. + This object may be a dict, or a custom Python object. + """ + def __init__(self, results_dict, solver_name: str) -> None: + self.solver_name = solver_name + self.solve_time = None + self.setup_time = None + self.num_iters = None + self.extra_stats = None + + if s.SOLVE_TIME in results_dict: + self.solve_time = results_dict[s.SOLVE_TIME] + if s.SETUP_TIME in results_dict: + self.setup_time = results_dict[s.SETUP_TIME] + if s.NUM_ITERS in results_dict: + self.num_iters = results_dict[s.NUM_ITERS] + if s.EXTRA_STATS in results_dict: + self.extra_stats = results_dict[s.EXTRA_STATS] + + +class SizeMetrics: + """Reports various metrics regarding the problem. + + Attributes + ---------- + + num_scalar_variables : integer + The number of scalar variables in the problem. + num_scalar_data : integer + The number of scalar constants and parameters in the problem. The number of + constants used across all matrices, vectors, in the problem. + Some constants are not apparent when the problem is constructed: for example, + The sum_squares expression is a wrapper for a quad_over_lin expression with a + constant 1 in the denominator. + num_scalar_eq_constr : integer + The number of scalar equality constraints in the problem. + num_scalar_leq_constr : integer + The number of scalar inequality constraints in the problem. + + max_data_dimension : integer + The longest dimension of any data block constraint or parameter. + max_big_small_squared : integer + The maximum value of (big)(small)^2 over all data blocks of the problem, where + (big) is the larger dimension and (small) is the smaller dimension + for each data block. + """ + + def __init__(self, problem: Problem) -> None: + # num_scalar_variables + self.num_scalar_variables = 0 + for var in problem.variables(): + self.num_scalar_variables += var.size + + # num_scalar_data, max_data_dimension, and max_big_small_squared + self.max_data_dimension = 0 + self.num_scalar_data = 0 + self.max_big_small_squared = 0 + for const in problem.constants()+problem.parameters(): + big = 0 + # Compute number of data + self.num_scalar_data += const.size + big = 1 if len(const.shape) == 0 else max(const.shape) + small = 1 if len(const.shape) == 0 else min(const.shape) + + # Get max data dimension: + if self.max_data_dimension < big: + self.max_data_dimension = big + + max_big_small_squared = float(big)*(float(small)**2) + if self.max_big_small_squared < max_big_small_squared: + self.max_big_small_squared = max_big_small_squared + + # num_scalar_eq_constr + self.num_scalar_eq_constr = 0 + for constraint in problem.constraints: + if isinstance(constraint, (Equality, Zero)): + self.num_scalar_eq_constr += constraint.expr.size + + # num_scalar_leq_constr + self.num_scalar_leq_constr = 0 + for constraint in problem.constraints: + if isinstance(constraint, (Inequality, NonPos, NonNeg)): + self.num_scalar_leq_constr += constraint.expr.size From 468c5913b2042864486b9291bccfa96b8892d2e7 Mon Sep 17 00:00:00 2001 From: rsemenoff Date: Tue, 22 Nov 2022 16:26:00 -0800 Subject: [PATCH 2/7] neos --- cvxpy/reductions/cplex_conif.py | 624 ++++++++++++++++++++++++++++++++ 1 file changed, 624 insertions(+) create mode 100644 cvxpy/reductions/cplex_conif.py diff --git a/cvxpy/reductions/cplex_conif.py b/cvxpy/reductions/cplex_conif.py new file mode 100644 index 0000000000..b9fd15b09e --- /dev/null +++ b/cvxpy/reductions/cplex_conif.py @@ -0,0 +1,624 @@ +""" +Copyright 2013 Steven Diamond, 2017 Robin Verschueren + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +""" + +from collections import namedtuple +from operator import attrgetter + +import numpy as np +from scipy.sparse import dok_matrix + +import cvxpy.settings as s +from cvxpy.constraints import SOC +from cvxpy.reductions.solution import Solution, failure_solution +from cvxpy.reductions.solvers import utilities +from cvxpy.reductions.solvers.conic_solvers.conic_solver import ( + ConicSolver, dims_to_solver_dict,) + +# Values used to distinguish between linear and quadratic constraints. +_LIN, _QUAD = 0, 1 +# For internal bookkeeping, we have to separate linear indices from +# quadratic indices. The "cpx_constrs" member of the solution will +# contain namedtuples of (constr_type, index) where constr_type is either +# _LIN or _QUAD. +_CpxConstr = namedtuple("_CpxConstr", ["constr_type", "index"]) + + +def set_parameters(model, solver_opts) -> None: + """Sets CPLEX parameters.""" + kwargs = sorted(solver_opts.keys()) + if "cplex_params" in kwargs: + for param, value in solver_opts["cplex_params"].items(): + try: + f = attrgetter(param) + parameter = f(model.parameters) + parameter.set(value) + except AttributeError: + raise ValueError( + "invalid CPLEX parameter, value pair ({0}, {1})".format( + param, value)) + kwargs.remove("cplex_params") + if "cplex_filename" in kwargs: + filename = solver_opts["cplex_filename"] + if filename: + model.write(filename) + kwargs.remove("cplex_filename") + if kwargs: + raise ValueError("invalid keyword-argument '{0}'".format(kwargs[0])) + + +def hide_solver_output(model) -> None: + """Set CPLEX verbosity level (either on or off).""" + # By default the output will be sent to stdout. Setting the output + # streams to None will prevent any output from being shown. + model.set_results_stream(None) + model.set_warning_stream(None) + model.set_error_stream(None) + model.set_log_stream(None) + + +def _handle_solve_status(model, solstat): + """Map CPLEX MIP solution status codes to non-MIP status codes.""" + status = model.solution.status + + # With CPLEX 12.10, some benders status codes were removed. + if model.get_versionnumber() < 12100000: + unbounded_status_codes = ( + status.MIP_unbounded, + status.MIP_benders_master_unbounded, + status.benders_master_unbounded + ) + else: + unbounded_status_codes = ( + status.MIP_unbounded, + ) + + if solstat == status.MIP_optimal: + return status.optimal + elif solstat == status.MIP_infeasible: + return status.infeasible + elif solstat in (status.MIP_time_limit_feasible, + status.MIP_time_limit_infeasible): + return status.abort_time_limit + elif solstat in (status.MIP_dettime_limit_feasible, + status.MIP_dettime_limit_infeasible): + return status.abort_dettime_limit + elif solstat in (status.MIP_abort_feasible, + status.MIP_abort_infeasible): + return status.abort_user + elif solstat == status.MIP_optimal_infeasible: + return status.optimal_infeasible + elif solstat == status.MIP_infeasible_or_unbounded: + return status.infeasible_or_unbounded + elif solstat in unbounded_status_codes: + return status.unbounded + elif solstat in (status.feasible_relaxed_sum, + status.MIP_feasible_relaxed_sum, + status.optimal_relaxed_sum, + status.MIP_optimal_relaxed_sum, + status.feasible_relaxed_inf, + status.MIP_feasible_relaxed_inf, + status.optimal_relaxed_inf, + status.MIP_optimal_relaxed_inf, + status.feasible_relaxed_quad, + status.MIP_feasible_relaxed_quad, + status.optimal_relaxed_quad, + status.MIP_optimal_relaxed_quad): + raise AssertionError( + "feasopt status encountered: {0}".format(solstat)) + elif solstat in (status.conflict_feasible, + status.conflict_minimal, + status.conflict_abort_contradiction, + status.conflict_abort_time_limit, + status.conflict_abort_dettime_limit, + status.conflict_abort_iteration_limit, + status.conflict_abort_node_limit, + status.conflict_abort_obj_limit, + status.conflict_abort_memory_limit, + status.conflict_abort_user): + raise AssertionError( + "conflict refiner status encountered: {0}".format(solstat)) + elif solstat == status.relaxation_unbounded: + return status.relaxation_unbounded + elif solstat in (status.feasible, + status.MIP_feasible): + return status.feasible + elif solstat == status.benders_num_best: + return status.num_best + else: + return solstat + + +def get_status(model): + """Map CPLEX status to CPXPY status.""" + pfeas = model.solution.is_primal_feasible() + # NOTE: dfeas is always false for a MIP. + dfeas = model.solution.is_dual_feasible() + status = model.solution.status + solstat = _handle_solve_status(model, model.solution.get_status()) + if solstat in (status.node_limit_infeasible, + status.fail_infeasible, + status.mem_limit_infeasible, + status.fail_infeasible_no_tree, + status.num_best): + return s.SOLVER_ERROR + elif solstat in (status.abort_user, + status.abort_iteration_limit, + status.abort_time_limit, + status.abort_dettime_limit, + status.abort_obj_limit, + status.abort_primal_obj_limit, + status.abort_dual_obj_limit, + status.abort_relaxed, + status.first_order): + if pfeas: + return s.OPTIMAL_INACCURATE + else: + return s.SOLVER_ERROR + elif solstat in (status.node_limit_feasible, + status.solution_limit, + status.populate_solution_limit, + status.fail_feasible, + status.mem_limit_feasible, + status.fail_feasible_no_tree, + status.feasible): + if dfeas: + return s.OPTIMAL + else: + return s.OPTIMAL_INACCURATE + elif solstat in (status.optimal, + status.optimal_tolerance, + status.optimal_infeasible, + status.optimal_populated, + status.optimal_populated_tolerance): + return s.OPTIMAL + elif solstat in (status.infeasible, + status.optimal_relaxed_sum, + status.optimal_relaxed_inf, + status.optimal_relaxed_quad): + return s.INFEASIBLE + elif solstat in (status.feasible_relaxed_quad, + status.feasible_relaxed_inf, + status.feasible_relaxed_sum): + return s.SOLVER_ERROR + elif solstat == status.unbounded: + return s.UNBOUNDED + elif solstat == status.infeasible_or_unbounded: + return s.INFEASIBLE_OR_UNBOUNDED + else: + return s.SOLVER_ERROR + + +class CPLEX(ConicSolver): + """An interface for the CPLEX solver. + """ + + # Solver capabilities. + MIP_CAPABLE = True + SUPPORTED_CONSTRAINTS = ConicSolver.SUPPORTED_CONSTRAINTS + [SOC] + + def name(self): + """The name of the solver. """ + return s.CPLEX + + def import_solver(self) -> None: + """Imports the solver.""" + import cplex + cplex # For flake8 + + def accepts(self, problem) -> bool: + """Can CPLEX solve the problem? + """ + # TODO check if is matrix stuffed. + if not problem.objective.args[0].is_affine(): + return False + for constr in problem.constraints: + if type(constr) not in CPLEX.SUPPORTED_CONSTRAINTS: + return False + for arg in constr.args: + if not arg.is_affine(): + return False + return True + + def apply(self, problem): + """Returns a new problem and data for inverting the new solution. + + Returns + ------- + tuple + (dict of arguments needed for the solver, inverse data) + """ + data, inv_data = super(CPLEX, self).apply(problem) + variables = problem.x + data[s.BOOL_IDX] = [int(t[0]) for t in variables.boolean_idx] + data[s.INT_IDX] = [int(t[0]) for t in variables.integer_idx] + inv_data['is_mip'] = data[s.BOOL_IDX] or data[s.INT_IDX] + + return data, inv_data + + def invert(self, solution, inverse_data): + """Returns the solution to the original problem given the inverse_data. + """ + model = solution["model"] + attr = {} + if s.SOLVE_TIME in solution: + attr[s.SOLVE_TIME] = solution[s.SOLVE_TIME] + + status = get_status(model) + + primal_vars = None + dual_vars = None + if status in s.SOLUTION_PRESENT: + opt_val = (model.solution.get_objective_value() + + inverse_data[s.OFFSET]) + + #x = np.array(model.solution.get_values())#2022-10-20 RWS read .sol file here.XML. + print("reading soln.sol;") + import xml.etree.ElementTree as ET + vd="c:/Users/Rwsin/myproject/" + tree = ET.parse(vd+"soln.sol")#"/mnt/c/Users/Rwsin/myproject/soln.sol") + root = tree.getroot() + x=np.array([child.attrib['value'] for child in root[3]]) + + + primal_vars = {inverse_data[CPLEX.VAR_ID]: x} + + if not inverse_data['is_mip']: + # The dual values are retrieved in the order that the + # constraints were added in solve_via_data() below. We + # must be careful to map them to inverse_data[EQ_CONSTR] + # followed by inverse_data[NEQ_CONSTR] accordingly. + y = -np.array(model.solution.get_dual_values()) + dual_vars = utilities.get_dual_values( + y, + utilities.extract_dual_value, + inverse_data[CPLEX.EQ_CONSTR] + inverse_data[CPLEX.NEQ_CONSTR]) + + return Solution(status, opt_val, primal_vars, dual_vars, attr) + else: + return failure_solution(status) + + def solve_via_data_1(self, data, warm_start: bool, verbose: bool, solver_opts, solver_cache=None):#2022-09-30 + print("We are here!") + import cplex + print("i can write code") + c = data[s.C] + b = data[s.B] + A = dok_matrix(data[s.A]) + # Save the dok_matrix. + data[s.A] = A + dims = dims_to_solver_dict(data[s.DIMS]) + + n = c.shape[0] + + model = cplex.Cplex() + variables = [] + # cpx_constrs will contain CpxConstr namedtuples (see above). + cpx_constrs = [] + vtype = [] + if data[s.BOOL_IDX] or data[s.INT_IDX]: + for i in range(n): + # Set variable type. + if i in data[s.BOOL_IDX]: + vtype.append('B') + elif i in data[s.INT_IDX]: + vtype.append('I') + else: + vtype.append('C') + else: + # If we specify types (even with 'C'), then the problem will + # be interpreted as a MIP. Leaving vtype as an empty list + # here, will ensure that the problem type remains an LP. + pass + # Add the variables in a batch + variables = list(model.variables.add( + obj=[c[i] for i in range(n)], + lb=[-cplex.infinity]*n, # default LB is 0 + ub=[cplex.infinity]*n, + types="".join(vtype), + names=["x_%d" % i for i in range(n)])) + + # Add equality constraints + cpx_constrs += [_CpxConstr(_LIN, x) + for x in self.add_model_lin_constr( + model, variables, + range(dims[s.EQ_DIM]), + 'E', A, b)] + + # Add inequality (<=) constraints + leq_start = dims[s.EQ_DIM] + leq_end = dims[s.EQ_DIM] + dims[s.LEQ_DIM] + cpx_constrs += [_CpxConstr(_LIN, x) + for x in self.add_model_lin_constr( + model, variables, + range(leq_start, leq_end), + 'L', A, b)] + + # Add SOC constraints + soc_start = leq_end + for constr_len in dims[s.SOC_DIM]: + soc_end = soc_start + constr_len + soc_constr, new_leq, new_vars = self.add_model_soc_constr( + model, variables, range(soc_start, soc_end), A, b) + cpx_constrs.append(_CpxConstr(_QUAD, soc_constr)) + cpx_constrs += [_CpxConstr(_LIN, x) for x in new_leq] + variables += new_vars + soc_start += constr_len + + # Set verbosity + if not verbose: + hide_solver_output(model) + + # For CVXPY, we set the qcpduals parameter here, but the user can + # easily override it via the "cplex_params" solver option (see + # set_parameters function). + model.parameters.preprocessing.qcpduals.set( + model.parameters.preprocessing.qcpduals.values.force) + + # Set parameters + reoptimize = solver_opts.pop('reoptimize', False) + set_parameters(model, solver_opts) + + # Solve problem + solution = {"model": model} + #RWS now write_lp() + model.write("svd.mps") + + def solve_via_data_2(self, data, warm_start: bool, verbose: bool, solver_opts, solver_cache=None): + + #model.solve() + #create an empty model + # load solution file + #solution[s.SOLVE_TIME] = model.get_time() - start_time + + # ambiguous_status = get_status(model) == s.INFEASIBLE_OR_UNBOUNDED + # if ambiguous_status and reoptimize: + # model.parameters.preprocessing.presolve.set(0) + # start_time = model.get_time() + # model.solve() + # solution[s.SOLVE_TIME] += model.get_time() - start_time + + # except Exception: + # pass + + # return solution + return + + + + + + def solve_via_data(self, data, warm_start: bool, verbose: bool, solver_opts, solver_cache=None): + import cplex + print("i can write code") + c = data[s.C] + b = data[s.B] + A = dok_matrix(data[s.A]) + # Save the dok_matrix. + data[s.A] = A + dims = dims_to_solver_dict(data[s.DIMS]) + + n = c.shape[0] + + model = cplex.Cplex() + variables = [] + # cpx_constrs will contain CpxConstr namedtuples (see above). + cpx_constrs = [] + vtype = [] + if data[s.BOOL_IDX] or data[s.INT_IDX]: + for i in range(n): + # Set variable type. + if i in data[s.BOOL_IDX]: + vtype.append('B') + elif i in data[s.INT_IDX]: + vtype.append('I') + else: + vtype.append('C') + else: + # If we specify types (even with 'C'), then the problem will + # be interpreted as a MIP. Leaving vtype as an empty list + # here, will ensure that the problem type remains an LP. + pass + # Add the variables in a batch + variables = list(model.variables.add( + obj=[c[i] for i in range(n)], + lb=[-cplex.infinity]*n, # default LB is 0 + ub=[cplex.infinity]*n, + types="".join(vtype), + names=["x_%d" % i for i in range(n)])) + + # Add equality constraints + cpx_constrs += [_CpxConstr(_LIN, x) + for x in self.add_model_lin_constr( + model, variables, + range(dims[s.EQ_DIM]), + 'E', A, b)] + + # Add inequality (<=) constraints + leq_start = dims[s.EQ_DIM] + leq_end = dims[s.EQ_DIM] + dims[s.LEQ_DIM] + cpx_constrs += [_CpxConstr(_LIN, x) + for x in self.add_model_lin_constr( + model, variables, + range(leq_start, leq_end), + 'L', A, b)] + + # Add SOC constraints + soc_start = leq_end + for constr_len in dims[s.SOC_DIM]: + soc_end = soc_start + constr_len + soc_constr, new_leq, new_vars = self.add_model_soc_constr( + model, variables, range(soc_start, soc_end), A, b) + cpx_constrs.append(_CpxConstr(_QUAD, soc_constr)) + cpx_constrs += [_CpxConstr(_LIN, x) for x in new_leq] + variables += new_vars + soc_start += constr_len + + # Set verbosity + if not verbose: + hide_solver_output(model) + + # For CVXPY, we set the qcpduals parameter here, but the user can + # easily override it via the "cplex_params" solver option (see + # set_parameters function). + model.parameters.preprocessing.qcpduals.set( + model.parameters.preprocessing.qcpduals.values.force) + + # Set parameters + reoptimize = solver_opts.pop('reoptimize', False) + set_parameters(model, solver_opts) + + # Solve problem + solution = {"model": model} + try: + start_time = model.get_time() + model.solve() + solution[s.SOLVE_TIME] = model.get_time() - start_time + + ambiguous_status = get_status(model) == s.INFEASIBLE_OR_UNBOUNDED + if ambiguous_status and reoptimize: + model.parameters.preprocessing.presolve.set(0) + start_time = model.get_time() + model.solve() + solution[s.SOLVE_TIME] += model.get_time() - start_time + + except Exception: + pass + + return solution + + def add_model_lin_constr(self, model, variables, + rows, ctype, mat, vec): + """Adds EQ/LEQ constraints to the model using the data from mat and vec. + + Parameters + ---------- + model : CPLEX model + The problem model. + variables : list + The problem variables. + rows : range + The rows to be constrained. + ctype : CPLEX constraint type + The type of constraint. + mat : SciPy COO matrix + The matrix representing the constraints. + vec : NDArray + The RHS part of the constraints. + + Returns + ------- + list + A list of new linear constraint indices. + """ + constr, lin_expr, rhs = [], [], [] + csr = mat.tocsr() + for i in rows: + ind = [variables[x] for x in csr[i].indices] + val = [x for x in csr[i].data] + lin_expr.append([ind, val]) + rhs.append(vec[i]) + # For better performance, we add the constraints in a batch. + if lin_expr: + assert len(lin_expr) == len(rhs) + constr.extend(list( + model.linear_constraints.add( + lin_expr=lin_expr, + senses=ctype * len(lin_expr), + rhs=rhs))) + return constr + + def add_model_soc_constr(self, model, variables, + rows, mat, vec): + """Adds SOC constraint to the model using the data from mat and vec. + + Parameters + ---------- + model : CPLEX model + The problem model. + variables : list + The problem variables. + rows : range + The rows to be constrained. + mat : SciPy COO matrix + The matrix representing the constraints. + vec : NDArray + The RHS part of the constraints. + + Returns + ------- + tuple + A tuple of (a new quadratic constraint index, a list of new + supporting linear constr indices, and a list of new + supporting variable indices). + """ + import cplex + + # Assume first expression (i.e. t) is nonzero. + lin_expr_list, soc_vars, lin_rhs = [], [], [] + csr = mat.tocsr() + for i in rows: + ind = [variables[x] for x in csr[i].indices] + val = [x for x in csr[i].data] + # Ignore empty constraints. + if ind: + lin_expr_list.append((ind, val)) + else: + lin_expr_list.append(None) + lin_rhs.append(vec[i]) + + # Make a variable and equality constraint for each term. + soc_vars, is_first = [], True + for i in rows: + if is_first: + lb = [0.0] + names = ["soc_t_%d" % i] + is_first = False + else: + lb = [-cplex.infinity] + names = ["soc_x_%d" % i] + soc_vars.extend(list(model.variables.add( + obj=[0], + lb=lb, + ub=[cplex.infinity], + types="", + names=names))) + + new_lin_constrs = [] + for i, expr in enumerate(lin_expr_list): + if expr is None: + ind = [soc_vars[i]] + val = [1.0] + else: + ind, val = expr + ind.append(soc_vars[i]) + val.append(1.0) + new_lin_constrs.extend(list( + model.linear_constraints.add( + lin_expr=[cplex.SparsePair(ind=ind, val=val)], + senses="E", + rhs=[lin_rhs[i]]))) + + assert len(soc_vars) > 0 + qconstr = model.quadratic_constraints.add( + lin_expr=cplex.SparsePair(ind=[], val=[]), + quad_expr=cplex.SparseTriple( + ind1=soc_vars, + ind2=soc_vars, + val=[-1.0] + [1.0] * (len(soc_vars) - 1)), + sense="L", + rhs=0.0, + name="") + return (qconstr, new_lin_constrs, soc_vars) From 6f2e51915b2c0c6a74fbcc8f39e56a571b02df86 Mon Sep 17 00:00:00 2001 From: rsemenoff Date: Tue, 22 Nov 2022 16:28:26 -0800 Subject: [PATCH 3/7] neos --- .../solvers/conic_solvers/cplex_conif.py | 1131 +++++++++-------- 1 file changed, 624 insertions(+), 507 deletions(-) diff --git a/cvxpy/reductions/solvers/conic_solvers/cplex_conif.py b/cvxpy/reductions/solvers/conic_solvers/cplex_conif.py index 6ba28b6a16..b9fd15b09e 100644 --- a/cvxpy/reductions/solvers/conic_solvers/cplex_conif.py +++ b/cvxpy/reductions/solvers/conic_solvers/cplex_conif.py @@ -1,507 +1,624 @@ -""" -Copyright 2013 Steven Diamond, 2017 Robin Verschueren - -Licensed under the Apache License, Version 2.0 (the "License"); -you may not use this file except in compliance with the License. -You may obtain a copy of the License at - - http://www.apache.org/licenses/LICENSE-2.0 - -Unless required by applicable law or agreed to in writing, software -distributed under the License is distributed on an "AS IS" BASIS, -WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -See the License for the specific language governing permissions and -limitations under the License. -""" - -from collections import namedtuple -from operator import attrgetter - -import numpy as np -from scipy.sparse import dok_matrix - -import cvxpy.settings as s -from cvxpy.constraints import SOC -from cvxpy.reductions.solution import Solution, failure_solution -from cvxpy.reductions.solvers import utilities -from cvxpy.reductions.solvers.conic_solvers.conic_solver import ( - ConicSolver, dims_to_solver_dict,) - -# Values used to distinguish between linear and quadratic constraints. -_LIN, _QUAD = 0, 1 -# For internal bookkeeping, we have to separate linear indices from -# quadratic indices. The "cpx_constrs" member of the solution will -# contain namedtuples of (constr_type, index) where constr_type is either -# _LIN or _QUAD. -_CpxConstr = namedtuple("_CpxConstr", ["constr_type", "index"]) - - -def set_parameters(model, solver_opts) -> None: - """Sets CPLEX parameters.""" - kwargs = sorted(solver_opts.keys()) - if "cplex_params" in kwargs: - for param, value in solver_opts["cplex_params"].items(): - try: - f = attrgetter(param) - parameter = f(model.parameters) - parameter.set(value) - except AttributeError: - raise ValueError( - "invalid CPLEX parameter, value pair ({0}, {1})".format( - param, value)) - kwargs.remove("cplex_params") - if "cplex_filename" in kwargs: - filename = solver_opts["cplex_filename"] - if filename: - model.write(filename) - kwargs.remove("cplex_filename") - if kwargs: - raise ValueError("invalid keyword-argument '{0}'".format(kwargs[0])) - - -def hide_solver_output(model) -> None: - """Set CPLEX verbosity level (either on or off).""" - # By default the output will be sent to stdout. Setting the output - # streams to None will prevent any output from being shown. - model.set_results_stream(None) - model.set_warning_stream(None) - model.set_error_stream(None) - model.set_log_stream(None) - - -def _handle_solve_status(model, solstat): - """Map CPLEX MIP solution status codes to non-MIP status codes.""" - status = model.solution.status - - # With CPLEX 12.10, some benders status codes were removed. - if model.get_versionnumber() < 12100000: - unbounded_status_codes = ( - status.MIP_unbounded, - status.MIP_benders_master_unbounded, - status.benders_master_unbounded - ) - else: - unbounded_status_codes = ( - status.MIP_unbounded, - ) - - if solstat == status.MIP_optimal: - return status.optimal - elif solstat == status.MIP_infeasible: - return status.infeasible - elif solstat in (status.MIP_time_limit_feasible, - status.MIP_time_limit_infeasible): - return status.abort_time_limit - elif solstat in (status.MIP_dettime_limit_feasible, - status.MIP_dettime_limit_infeasible): - return status.abort_dettime_limit - elif solstat in (status.MIP_abort_feasible, - status.MIP_abort_infeasible): - return status.abort_user - elif solstat == status.MIP_optimal_infeasible: - return status.optimal_infeasible - elif solstat == status.MIP_infeasible_or_unbounded: - return status.infeasible_or_unbounded - elif solstat in unbounded_status_codes: - return status.unbounded - elif solstat in (status.feasible_relaxed_sum, - status.MIP_feasible_relaxed_sum, - status.optimal_relaxed_sum, - status.MIP_optimal_relaxed_sum, - status.feasible_relaxed_inf, - status.MIP_feasible_relaxed_inf, - status.optimal_relaxed_inf, - status.MIP_optimal_relaxed_inf, - status.feasible_relaxed_quad, - status.MIP_feasible_relaxed_quad, - status.optimal_relaxed_quad, - status.MIP_optimal_relaxed_quad): - raise AssertionError( - "feasopt status encountered: {0}".format(solstat)) - elif solstat in (status.conflict_feasible, - status.conflict_minimal, - status.conflict_abort_contradiction, - status.conflict_abort_time_limit, - status.conflict_abort_dettime_limit, - status.conflict_abort_iteration_limit, - status.conflict_abort_node_limit, - status.conflict_abort_obj_limit, - status.conflict_abort_memory_limit, - status.conflict_abort_user): - raise AssertionError( - "conflict refiner status encountered: {0}".format(solstat)) - elif solstat == status.relaxation_unbounded: - return status.relaxation_unbounded - elif solstat in (status.feasible, - status.MIP_feasible): - return status.feasible - elif solstat == status.benders_num_best: - return status.num_best - else: - return solstat - - -def get_status(model): - """Map CPLEX status to CPXPY status.""" - pfeas = model.solution.is_primal_feasible() - # NOTE: dfeas is always false for a MIP. - dfeas = model.solution.is_dual_feasible() - status = model.solution.status - solstat = _handle_solve_status(model, model.solution.get_status()) - if solstat in (status.node_limit_infeasible, - status.fail_infeasible, - status.mem_limit_infeasible, - status.fail_infeasible_no_tree, - status.num_best): - return s.SOLVER_ERROR - elif solstat in (status.abort_user, - status.abort_iteration_limit, - status.abort_time_limit, - status.abort_dettime_limit, - status.abort_obj_limit, - status.abort_primal_obj_limit, - status.abort_dual_obj_limit, - status.abort_relaxed, - status.first_order): - if pfeas: - return s.OPTIMAL_INACCURATE - else: - return s.SOLVER_ERROR - elif solstat in (status.node_limit_feasible, - status.solution_limit, - status.populate_solution_limit, - status.fail_feasible, - status.mem_limit_feasible, - status.fail_feasible_no_tree, - status.feasible): - if dfeas: - return s.OPTIMAL - else: - return s.OPTIMAL_INACCURATE - elif solstat in (status.optimal, - status.optimal_tolerance, - status.optimal_infeasible, - status.optimal_populated, - status.optimal_populated_tolerance): - return s.OPTIMAL - elif solstat in (status.infeasible, - status.optimal_relaxed_sum, - status.optimal_relaxed_inf, - status.optimal_relaxed_quad): - return s.INFEASIBLE - elif solstat in (status.feasible_relaxed_quad, - status.feasible_relaxed_inf, - status.feasible_relaxed_sum): - return s.SOLVER_ERROR - elif solstat == status.unbounded: - return s.UNBOUNDED - elif solstat == status.infeasible_or_unbounded: - return s.INFEASIBLE_OR_UNBOUNDED - else: - return s.SOLVER_ERROR - - -class CPLEX(ConicSolver): - """An interface for the CPLEX solver. - """ - - # Solver capabilities. - MIP_CAPABLE = True - SUPPORTED_CONSTRAINTS = ConicSolver.SUPPORTED_CONSTRAINTS + [SOC] - MI_SUPPORTED_CONSTRAINTS = SUPPORTED_CONSTRAINTS - - def name(self): - """The name of the solver. """ - return s.CPLEX - - def import_solver(self) -> None: - """Imports the solver.""" - import cplex - cplex # For flake8 - - def accepts(self, problem) -> bool: - """Can CPLEX solve the problem? - """ - # TODO check if is matrix stuffed. - if not problem.objective.args[0].is_affine(): - return False - for constr in problem.constraints: - if type(constr) not in CPLEX.SUPPORTED_CONSTRAINTS: - return False - for arg in constr.args: - if not arg.is_affine(): - return False - return True - - def apply(self, problem): - """Returns a new problem and data for inverting the new solution. - - Returns - ------- - tuple - (dict of arguments needed for the solver, inverse data) - """ - data, inv_data = super(CPLEX, self).apply(problem) - variables = problem.x - data[s.BOOL_IDX] = [int(t[0]) for t in variables.boolean_idx] - data[s.INT_IDX] = [int(t[0]) for t in variables.integer_idx] - inv_data['is_mip'] = data[s.BOOL_IDX] or data[s.INT_IDX] - - return data, inv_data - - def invert(self, solution, inverse_data): - """Returns the solution to the original problem given the inverse_data. - """ - model = solution["model"] - attr = {} - if s.SOLVE_TIME in solution: - attr[s.SOLVE_TIME] = solution[s.SOLVE_TIME] - - status = get_status(model) - - primal_vars = None - dual_vars = None - if status in s.SOLUTION_PRESENT: - opt_val = (model.solution.get_objective_value() + - inverse_data[s.OFFSET]) - - x = np.array(model.solution.get_values()) - primal_vars = {inverse_data[CPLEX.VAR_ID]: x} - - if not inverse_data['is_mip']: - # The dual values are retrieved in the order that the - # constraints were added in solve_via_data() below. We - # must be careful to map them to inverse_data[EQ_CONSTR] - # followed by inverse_data[NEQ_CONSTR] accordingly. - y = -np.array(model.solution.get_dual_values()) - dual_vars = utilities.get_dual_values( - y, - utilities.extract_dual_value, - inverse_data[CPLEX.EQ_CONSTR] + inverse_data[CPLEX.NEQ_CONSTR]) - - return Solution(status, opt_val, primal_vars, dual_vars, attr) - else: - return failure_solution(status) - - def solve_via_data(self, data, warm_start: bool, verbose: bool, solver_opts, solver_cache=None): - import cplex - - c = data[s.C] - b = data[s.B] - A = dok_matrix(data[s.A]) - # Save the dok_matrix. - data[s.A] = A - dims = dims_to_solver_dict(data[s.DIMS]) - - n = c.shape[0] - - model = cplex.Cplex() - variables = [] - # cpx_constrs will contain CpxConstr namedtuples (see above). - cpx_constrs = [] - vtype = [] - if data[s.BOOL_IDX] or data[s.INT_IDX]: - for i in range(n): - # Set variable type. - if i in data[s.BOOL_IDX]: - vtype.append('B') - elif i in data[s.INT_IDX]: - vtype.append('I') - else: - vtype.append('C') - else: - # If we specify types (even with 'C'), then the problem will - # be interpreted as a MIP. Leaving vtype as an empty list - # here, will ensure that the problem type remains an LP. - pass - # Add the variables in a batch - variables = list(model.variables.add( - obj=[c[i] for i in range(n)], - lb=[-cplex.infinity]*n, # default LB is 0 - ub=[cplex.infinity]*n, - types="".join(vtype), - names=["x_%d" % i for i in range(n)])) - - # Add equality constraints - cpx_constrs += [_CpxConstr(_LIN, x) - for x in self.add_model_lin_constr( - model, variables, - range(dims[s.EQ_DIM]), - 'E', A, b)] - - # Add inequality (<=) constraints - leq_start = dims[s.EQ_DIM] - leq_end = dims[s.EQ_DIM] + dims[s.LEQ_DIM] - cpx_constrs += [_CpxConstr(_LIN, x) - for x in self.add_model_lin_constr( - model, variables, - range(leq_start, leq_end), - 'L', A, b)] - - # Add SOC constraints - soc_start = leq_end - for constr_len in dims[s.SOC_DIM]: - soc_end = soc_start + constr_len - soc_constr, new_leq, new_vars = self.add_model_soc_constr( - model, variables, range(soc_start, soc_end), A, b) - cpx_constrs.append(_CpxConstr(_QUAD, soc_constr)) - cpx_constrs += [_CpxConstr(_LIN, x) for x in new_leq] - variables += new_vars - soc_start += constr_len - - # Set verbosity - if not verbose: - hide_solver_output(model) - - # For CVXPY, we set the qcpduals parameter here, but the user can - # easily override it via the "cplex_params" solver option (see - # set_parameters function). - model.parameters.preprocessing.qcpduals.set( - model.parameters.preprocessing.qcpduals.values.force) - - # Set parameters - reoptimize = solver_opts.pop('reoptimize', False) - set_parameters(model, solver_opts) - - # Solve problem - solution = {"model": model} - try: - start_time = model.get_time() - model.solve() - solution[s.SOLVE_TIME] = model.get_time() - start_time - - ambiguous_status = get_status(model) == s.INFEASIBLE_OR_UNBOUNDED - if ambiguous_status and reoptimize: - model.parameters.preprocessing.presolve.set(0) - start_time = model.get_time() - model.solve() - solution[s.SOLVE_TIME] += model.get_time() - start_time - - except Exception: - pass - - return solution - - def add_model_lin_constr(self, model, variables, - rows, ctype, mat, vec): - """Adds EQ/LEQ constraints to the model using the data from mat and vec. - - Parameters - ---------- - model : CPLEX model - The problem model. - variables : list - The problem variables. - rows : range - The rows to be constrained. - ctype : CPLEX constraint type - The type of constraint. - mat : SciPy COO matrix - The matrix representing the constraints. - vec : NDArray - The RHS part of the constraints. - - Returns - ------- - list - A list of new linear constraint indices. - """ - constr, lin_expr, rhs = [], [], [] - csr = mat.tocsr() - for i in rows: - ind = [variables[x] for x in csr[i].indices] - val = [x for x in csr[i].data] - lin_expr.append([ind, val]) - rhs.append(vec[i]) - # For better performance, we add the constraints in a batch. - if lin_expr: - assert len(lin_expr) == len(rhs) - constr.extend(list( - model.linear_constraints.add( - lin_expr=lin_expr, - senses=ctype * len(lin_expr), - rhs=rhs))) - return constr - - def add_model_soc_constr(self, model, variables, - rows, mat, vec): - """Adds SOC constraint to the model using the data from mat and vec. - - Parameters - ---------- - model : CPLEX model - The problem model. - variables : list - The problem variables. - rows : range - The rows to be constrained. - mat : SciPy COO matrix - The matrix representing the constraints. - vec : NDArray - The RHS part of the constraints. - - Returns - ------- - tuple - A tuple of (a new quadratic constraint index, a list of new - supporting linear constr indices, and a list of new - supporting variable indices). - """ - import cplex - - # Assume first expression (i.e. t) is nonzero. - lin_expr_list, soc_vars, lin_rhs = [], [], [] - csr = mat.tocsr() - for i in rows: - ind = [variables[x] for x in csr[i].indices] - val = [x for x in csr[i].data] - # Ignore empty constraints. - if ind: - lin_expr_list.append((ind, val)) - else: - lin_expr_list.append(None) - lin_rhs.append(vec[i]) - - # Make a variable and equality constraint for each term. - soc_vars, is_first = [], True - for i in rows: - if is_first: - lb = [0.0] - names = ["soc_t_%d" % i] - is_first = False - else: - lb = [-cplex.infinity] - names = ["soc_x_%d" % i] - soc_vars.extend(list(model.variables.add( - obj=[0], - lb=lb, - ub=[cplex.infinity], - types="", - names=names))) - - new_lin_constrs = [] - for i, expr in enumerate(lin_expr_list): - if expr is None: - ind = [soc_vars[i]] - val = [1.0] - else: - ind, val = expr - ind.append(soc_vars[i]) - val.append(1.0) - new_lin_constrs.extend(list( - model.linear_constraints.add( - lin_expr=[cplex.SparsePair(ind=ind, val=val)], - senses="E", - rhs=[lin_rhs[i]]))) - - assert len(soc_vars) > 0 - qconstr = model.quadratic_constraints.add( - lin_expr=cplex.SparsePair(ind=[], val=[]), - quad_expr=cplex.SparseTriple( - ind1=soc_vars, - ind2=soc_vars, - val=[-1.0] + [1.0] * (len(soc_vars) - 1)), - sense="L", - rhs=0.0, - name="") - return (qconstr, new_lin_constrs, soc_vars) +""" +Copyright 2013 Steven Diamond, 2017 Robin Verschueren + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +""" + +from collections import namedtuple +from operator import attrgetter + +import numpy as np +from scipy.sparse import dok_matrix + +import cvxpy.settings as s +from cvxpy.constraints import SOC +from cvxpy.reductions.solution import Solution, failure_solution +from cvxpy.reductions.solvers import utilities +from cvxpy.reductions.solvers.conic_solvers.conic_solver import ( + ConicSolver, dims_to_solver_dict,) + +# Values used to distinguish between linear and quadratic constraints. +_LIN, _QUAD = 0, 1 +# For internal bookkeeping, we have to separate linear indices from +# quadratic indices. The "cpx_constrs" member of the solution will +# contain namedtuples of (constr_type, index) where constr_type is either +# _LIN or _QUAD. +_CpxConstr = namedtuple("_CpxConstr", ["constr_type", "index"]) + + +def set_parameters(model, solver_opts) -> None: + """Sets CPLEX parameters.""" + kwargs = sorted(solver_opts.keys()) + if "cplex_params" in kwargs: + for param, value in solver_opts["cplex_params"].items(): + try: + f = attrgetter(param) + parameter = f(model.parameters) + parameter.set(value) + except AttributeError: + raise ValueError( + "invalid CPLEX parameter, value pair ({0}, {1})".format( + param, value)) + kwargs.remove("cplex_params") + if "cplex_filename" in kwargs: + filename = solver_opts["cplex_filename"] + if filename: + model.write(filename) + kwargs.remove("cplex_filename") + if kwargs: + raise ValueError("invalid keyword-argument '{0}'".format(kwargs[0])) + + +def hide_solver_output(model) -> None: + """Set CPLEX verbosity level (either on or off).""" + # By default the output will be sent to stdout. Setting the output + # streams to None will prevent any output from being shown. + model.set_results_stream(None) + model.set_warning_stream(None) + model.set_error_stream(None) + model.set_log_stream(None) + + +def _handle_solve_status(model, solstat): + """Map CPLEX MIP solution status codes to non-MIP status codes.""" + status = model.solution.status + + # With CPLEX 12.10, some benders status codes were removed. + if model.get_versionnumber() < 12100000: + unbounded_status_codes = ( + status.MIP_unbounded, + status.MIP_benders_master_unbounded, + status.benders_master_unbounded + ) + else: + unbounded_status_codes = ( + status.MIP_unbounded, + ) + + if solstat == status.MIP_optimal: + return status.optimal + elif solstat == status.MIP_infeasible: + return status.infeasible + elif solstat in (status.MIP_time_limit_feasible, + status.MIP_time_limit_infeasible): + return status.abort_time_limit + elif solstat in (status.MIP_dettime_limit_feasible, + status.MIP_dettime_limit_infeasible): + return status.abort_dettime_limit + elif solstat in (status.MIP_abort_feasible, + status.MIP_abort_infeasible): + return status.abort_user + elif solstat == status.MIP_optimal_infeasible: + return status.optimal_infeasible + elif solstat == status.MIP_infeasible_or_unbounded: + return status.infeasible_or_unbounded + elif solstat in unbounded_status_codes: + return status.unbounded + elif solstat in (status.feasible_relaxed_sum, + status.MIP_feasible_relaxed_sum, + status.optimal_relaxed_sum, + status.MIP_optimal_relaxed_sum, + status.feasible_relaxed_inf, + status.MIP_feasible_relaxed_inf, + status.optimal_relaxed_inf, + status.MIP_optimal_relaxed_inf, + status.feasible_relaxed_quad, + status.MIP_feasible_relaxed_quad, + status.optimal_relaxed_quad, + status.MIP_optimal_relaxed_quad): + raise AssertionError( + "feasopt status encountered: {0}".format(solstat)) + elif solstat in (status.conflict_feasible, + status.conflict_minimal, + status.conflict_abort_contradiction, + status.conflict_abort_time_limit, + status.conflict_abort_dettime_limit, + status.conflict_abort_iteration_limit, + status.conflict_abort_node_limit, + status.conflict_abort_obj_limit, + status.conflict_abort_memory_limit, + status.conflict_abort_user): + raise AssertionError( + "conflict refiner status encountered: {0}".format(solstat)) + elif solstat == status.relaxation_unbounded: + return status.relaxation_unbounded + elif solstat in (status.feasible, + status.MIP_feasible): + return status.feasible + elif solstat == status.benders_num_best: + return status.num_best + else: + return solstat + + +def get_status(model): + """Map CPLEX status to CPXPY status.""" + pfeas = model.solution.is_primal_feasible() + # NOTE: dfeas is always false for a MIP. + dfeas = model.solution.is_dual_feasible() + status = model.solution.status + solstat = _handle_solve_status(model, model.solution.get_status()) + if solstat in (status.node_limit_infeasible, + status.fail_infeasible, + status.mem_limit_infeasible, + status.fail_infeasible_no_tree, + status.num_best): + return s.SOLVER_ERROR + elif solstat in (status.abort_user, + status.abort_iteration_limit, + status.abort_time_limit, + status.abort_dettime_limit, + status.abort_obj_limit, + status.abort_primal_obj_limit, + status.abort_dual_obj_limit, + status.abort_relaxed, + status.first_order): + if pfeas: + return s.OPTIMAL_INACCURATE + else: + return s.SOLVER_ERROR + elif solstat in (status.node_limit_feasible, + status.solution_limit, + status.populate_solution_limit, + status.fail_feasible, + status.mem_limit_feasible, + status.fail_feasible_no_tree, + status.feasible): + if dfeas: + return s.OPTIMAL + else: + return s.OPTIMAL_INACCURATE + elif solstat in (status.optimal, + status.optimal_tolerance, + status.optimal_infeasible, + status.optimal_populated, + status.optimal_populated_tolerance): + return s.OPTIMAL + elif solstat in (status.infeasible, + status.optimal_relaxed_sum, + status.optimal_relaxed_inf, + status.optimal_relaxed_quad): + return s.INFEASIBLE + elif solstat in (status.feasible_relaxed_quad, + status.feasible_relaxed_inf, + status.feasible_relaxed_sum): + return s.SOLVER_ERROR + elif solstat == status.unbounded: + return s.UNBOUNDED + elif solstat == status.infeasible_or_unbounded: + return s.INFEASIBLE_OR_UNBOUNDED + else: + return s.SOLVER_ERROR + + +class CPLEX(ConicSolver): + """An interface for the CPLEX solver. + """ + + # Solver capabilities. + MIP_CAPABLE = True + SUPPORTED_CONSTRAINTS = ConicSolver.SUPPORTED_CONSTRAINTS + [SOC] + + def name(self): + """The name of the solver. """ + return s.CPLEX + + def import_solver(self) -> None: + """Imports the solver.""" + import cplex + cplex # For flake8 + + def accepts(self, problem) -> bool: + """Can CPLEX solve the problem? + """ + # TODO check if is matrix stuffed. + if not problem.objective.args[0].is_affine(): + return False + for constr in problem.constraints: + if type(constr) not in CPLEX.SUPPORTED_CONSTRAINTS: + return False + for arg in constr.args: + if not arg.is_affine(): + return False + return True + + def apply(self, problem): + """Returns a new problem and data for inverting the new solution. + + Returns + ------- + tuple + (dict of arguments needed for the solver, inverse data) + """ + data, inv_data = super(CPLEX, self).apply(problem) + variables = problem.x + data[s.BOOL_IDX] = [int(t[0]) for t in variables.boolean_idx] + data[s.INT_IDX] = [int(t[0]) for t in variables.integer_idx] + inv_data['is_mip'] = data[s.BOOL_IDX] or data[s.INT_IDX] + + return data, inv_data + + def invert(self, solution, inverse_data): + """Returns the solution to the original problem given the inverse_data. + """ + model = solution["model"] + attr = {} + if s.SOLVE_TIME in solution: + attr[s.SOLVE_TIME] = solution[s.SOLVE_TIME] + + status = get_status(model) + + primal_vars = None + dual_vars = None + if status in s.SOLUTION_PRESENT: + opt_val = (model.solution.get_objective_value() + + inverse_data[s.OFFSET]) + + #x = np.array(model.solution.get_values())#2022-10-20 RWS read .sol file here.XML. + print("reading soln.sol;") + import xml.etree.ElementTree as ET + vd="c:/Users/Rwsin/myproject/" + tree = ET.parse(vd+"soln.sol")#"/mnt/c/Users/Rwsin/myproject/soln.sol") + root = tree.getroot() + x=np.array([child.attrib['value'] for child in root[3]]) + + + primal_vars = {inverse_data[CPLEX.VAR_ID]: x} + + if not inverse_data['is_mip']: + # The dual values are retrieved in the order that the + # constraints were added in solve_via_data() below. We + # must be careful to map them to inverse_data[EQ_CONSTR] + # followed by inverse_data[NEQ_CONSTR] accordingly. + y = -np.array(model.solution.get_dual_values()) + dual_vars = utilities.get_dual_values( + y, + utilities.extract_dual_value, + inverse_data[CPLEX.EQ_CONSTR] + inverse_data[CPLEX.NEQ_CONSTR]) + + return Solution(status, opt_val, primal_vars, dual_vars, attr) + else: + return failure_solution(status) + + def solve_via_data_1(self, data, warm_start: bool, verbose: bool, solver_opts, solver_cache=None):#2022-09-30 + print("We are here!") + import cplex + print("i can write code") + c = data[s.C] + b = data[s.B] + A = dok_matrix(data[s.A]) + # Save the dok_matrix. + data[s.A] = A + dims = dims_to_solver_dict(data[s.DIMS]) + + n = c.shape[0] + + model = cplex.Cplex() + variables = [] + # cpx_constrs will contain CpxConstr namedtuples (see above). + cpx_constrs = [] + vtype = [] + if data[s.BOOL_IDX] or data[s.INT_IDX]: + for i in range(n): + # Set variable type. + if i in data[s.BOOL_IDX]: + vtype.append('B') + elif i in data[s.INT_IDX]: + vtype.append('I') + else: + vtype.append('C') + else: + # If we specify types (even with 'C'), then the problem will + # be interpreted as a MIP. Leaving vtype as an empty list + # here, will ensure that the problem type remains an LP. + pass + # Add the variables in a batch + variables = list(model.variables.add( + obj=[c[i] for i in range(n)], + lb=[-cplex.infinity]*n, # default LB is 0 + ub=[cplex.infinity]*n, + types="".join(vtype), + names=["x_%d" % i for i in range(n)])) + + # Add equality constraints + cpx_constrs += [_CpxConstr(_LIN, x) + for x in self.add_model_lin_constr( + model, variables, + range(dims[s.EQ_DIM]), + 'E', A, b)] + + # Add inequality (<=) constraints + leq_start = dims[s.EQ_DIM] + leq_end = dims[s.EQ_DIM] + dims[s.LEQ_DIM] + cpx_constrs += [_CpxConstr(_LIN, x) + for x in self.add_model_lin_constr( + model, variables, + range(leq_start, leq_end), + 'L', A, b)] + + # Add SOC constraints + soc_start = leq_end + for constr_len in dims[s.SOC_DIM]: + soc_end = soc_start + constr_len + soc_constr, new_leq, new_vars = self.add_model_soc_constr( + model, variables, range(soc_start, soc_end), A, b) + cpx_constrs.append(_CpxConstr(_QUAD, soc_constr)) + cpx_constrs += [_CpxConstr(_LIN, x) for x in new_leq] + variables += new_vars + soc_start += constr_len + + # Set verbosity + if not verbose: + hide_solver_output(model) + + # For CVXPY, we set the qcpduals parameter here, but the user can + # easily override it via the "cplex_params" solver option (see + # set_parameters function). + model.parameters.preprocessing.qcpduals.set( + model.parameters.preprocessing.qcpduals.values.force) + + # Set parameters + reoptimize = solver_opts.pop('reoptimize', False) + set_parameters(model, solver_opts) + + # Solve problem + solution = {"model": model} + #RWS now write_lp() + model.write("svd.mps") + + def solve_via_data_2(self, data, warm_start: bool, verbose: bool, solver_opts, solver_cache=None): + + #model.solve() + #create an empty model + # load solution file + #solution[s.SOLVE_TIME] = model.get_time() - start_time + + # ambiguous_status = get_status(model) == s.INFEASIBLE_OR_UNBOUNDED + # if ambiguous_status and reoptimize: + # model.parameters.preprocessing.presolve.set(0) + # start_time = model.get_time() + # model.solve() + # solution[s.SOLVE_TIME] += model.get_time() - start_time + + # except Exception: + # pass + + # return solution + return + + + + + + def solve_via_data(self, data, warm_start: bool, verbose: bool, solver_opts, solver_cache=None): + import cplex + print("i can write code") + c = data[s.C] + b = data[s.B] + A = dok_matrix(data[s.A]) + # Save the dok_matrix. + data[s.A] = A + dims = dims_to_solver_dict(data[s.DIMS]) + + n = c.shape[0] + + model = cplex.Cplex() + variables = [] + # cpx_constrs will contain CpxConstr namedtuples (see above). + cpx_constrs = [] + vtype = [] + if data[s.BOOL_IDX] or data[s.INT_IDX]: + for i in range(n): + # Set variable type. + if i in data[s.BOOL_IDX]: + vtype.append('B') + elif i in data[s.INT_IDX]: + vtype.append('I') + else: + vtype.append('C') + else: + # If we specify types (even with 'C'), then the problem will + # be interpreted as a MIP. Leaving vtype as an empty list + # here, will ensure that the problem type remains an LP. + pass + # Add the variables in a batch + variables = list(model.variables.add( + obj=[c[i] for i in range(n)], + lb=[-cplex.infinity]*n, # default LB is 0 + ub=[cplex.infinity]*n, + types="".join(vtype), + names=["x_%d" % i for i in range(n)])) + + # Add equality constraints + cpx_constrs += [_CpxConstr(_LIN, x) + for x in self.add_model_lin_constr( + model, variables, + range(dims[s.EQ_DIM]), + 'E', A, b)] + + # Add inequality (<=) constraints + leq_start = dims[s.EQ_DIM] + leq_end = dims[s.EQ_DIM] + dims[s.LEQ_DIM] + cpx_constrs += [_CpxConstr(_LIN, x) + for x in self.add_model_lin_constr( + model, variables, + range(leq_start, leq_end), + 'L', A, b)] + + # Add SOC constraints + soc_start = leq_end + for constr_len in dims[s.SOC_DIM]: + soc_end = soc_start + constr_len + soc_constr, new_leq, new_vars = self.add_model_soc_constr( + model, variables, range(soc_start, soc_end), A, b) + cpx_constrs.append(_CpxConstr(_QUAD, soc_constr)) + cpx_constrs += [_CpxConstr(_LIN, x) for x in new_leq] + variables += new_vars + soc_start += constr_len + + # Set verbosity + if not verbose: + hide_solver_output(model) + + # For CVXPY, we set the qcpduals parameter here, but the user can + # easily override it via the "cplex_params" solver option (see + # set_parameters function). + model.parameters.preprocessing.qcpduals.set( + model.parameters.preprocessing.qcpduals.values.force) + + # Set parameters + reoptimize = solver_opts.pop('reoptimize', False) + set_parameters(model, solver_opts) + + # Solve problem + solution = {"model": model} + try: + start_time = model.get_time() + model.solve() + solution[s.SOLVE_TIME] = model.get_time() - start_time + + ambiguous_status = get_status(model) == s.INFEASIBLE_OR_UNBOUNDED + if ambiguous_status and reoptimize: + model.parameters.preprocessing.presolve.set(0) + start_time = model.get_time() + model.solve() + solution[s.SOLVE_TIME] += model.get_time() - start_time + + except Exception: + pass + + return solution + + def add_model_lin_constr(self, model, variables, + rows, ctype, mat, vec): + """Adds EQ/LEQ constraints to the model using the data from mat and vec. + + Parameters + ---------- + model : CPLEX model + The problem model. + variables : list + The problem variables. + rows : range + The rows to be constrained. + ctype : CPLEX constraint type + The type of constraint. + mat : SciPy COO matrix + The matrix representing the constraints. + vec : NDArray + The RHS part of the constraints. + + Returns + ------- + list + A list of new linear constraint indices. + """ + constr, lin_expr, rhs = [], [], [] + csr = mat.tocsr() + for i in rows: + ind = [variables[x] for x in csr[i].indices] + val = [x for x in csr[i].data] + lin_expr.append([ind, val]) + rhs.append(vec[i]) + # For better performance, we add the constraints in a batch. + if lin_expr: + assert len(lin_expr) == len(rhs) + constr.extend(list( + model.linear_constraints.add( + lin_expr=lin_expr, + senses=ctype * len(lin_expr), + rhs=rhs))) + return constr + + def add_model_soc_constr(self, model, variables, + rows, mat, vec): + """Adds SOC constraint to the model using the data from mat and vec. + + Parameters + ---------- + model : CPLEX model + The problem model. + variables : list + The problem variables. + rows : range + The rows to be constrained. + mat : SciPy COO matrix + The matrix representing the constraints. + vec : NDArray + The RHS part of the constraints. + + Returns + ------- + tuple + A tuple of (a new quadratic constraint index, a list of new + supporting linear constr indices, and a list of new + supporting variable indices). + """ + import cplex + + # Assume first expression (i.e. t) is nonzero. + lin_expr_list, soc_vars, lin_rhs = [], [], [] + csr = mat.tocsr() + for i in rows: + ind = [variables[x] for x in csr[i].indices] + val = [x for x in csr[i].data] + # Ignore empty constraints. + if ind: + lin_expr_list.append((ind, val)) + else: + lin_expr_list.append(None) + lin_rhs.append(vec[i]) + + # Make a variable and equality constraint for each term. + soc_vars, is_first = [], True + for i in rows: + if is_first: + lb = [0.0] + names = ["soc_t_%d" % i] + is_first = False + else: + lb = [-cplex.infinity] + names = ["soc_x_%d" % i] + soc_vars.extend(list(model.variables.add( + obj=[0], + lb=lb, + ub=[cplex.infinity], + types="", + names=names))) + + new_lin_constrs = [] + for i, expr in enumerate(lin_expr_list): + if expr is None: + ind = [soc_vars[i]] + val = [1.0] + else: + ind, val = expr + ind.append(soc_vars[i]) + val.append(1.0) + new_lin_constrs.extend(list( + model.linear_constraints.add( + lin_expr=[cplex.SparsePair(ind=ind, val=val)], + senses="E", + rhs=[lin_rhs[i]]))) + + assert len(soc_vars) > 0 + qconstr = model.quadratic_constraints.add( + lin_expr=cplex.SparsePair(ind=[], val=[]), + quad_expr=cplex.SparseTriple( + ind1=soc_vars, + ind2=soc_vars, + val=[-1.0] + [1.0] * (len(soc_vars) - 1)), + sense="L", + rhs=0.0, + name="") + return (qconstr, new_lin_constrs, soc_vars) From 13f5898e4a8e1665ce44447a89d9151167f417b5 Mon Sep 17 00:00:00 2001 From: rsemenoff Date: Sat, 26 Nov 2022 19:33:31 -0800 Subject: [PATCH 4/7] Add files via upload neos --- .../solvers/qp_solvers/cplex_conif.py | 623 ++++++++++++++++++ 1 file changed, 623 insertions(+) create mode 100644 cvxpy/reductions/solvers/qp_solvers/cplex_conif.py diff --git a/cvxpy/reductions/solvers/qp_solvers/cplex_conif.py b/cvxpy/reductions/solvers/qp_solvers/cplex_conif.py new file mode 100644 index 0000000000..1661753e67 --- /dev/null +++ b/cvxpy/reductions/solvers/qp_solvers/cplex_conif.py @@ -0,0 +1,623 @@ +""" +Copyright 2013 Steven Diamond, 2017 Robin Verschueren + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +""" + +from collections import namedtuple +from operator import attrgetter + +import numpy as np +from scipy.sparse import dok_matrix + +import cvxpy.settings as s +from cvxpy.constraints import SOC +from cvxpy.reductions.solution import Solution, failure_solution +from cvxpy.reductions.solvers import utilities +from cvxpy.reductions.solvers.conic_solvers.conic_solver import ( + ConicSolver, dims_to_solver_dict,) + +# Values used to distinguish between linear and quadratic constraints. +_LIN, _QUAD = 0, 1 +# For internal bookkeeping, we have to separate linear indices from +# quadratic indices. The "cpx_constrs" member of the solution will +# contain namedtuples of (constr_type, index) where constr_type is either +# _LIN or _QUAD. +_CpxConstr = namedtuple("_CpxConstr", ["constr_type", "index"]) + + +def set_parameters(model, solver_opts) -> None: + """Sets CPLEX parameters.""" + kwargs = sorted(solver_opts.keys()) + if "cplex_params" in kwargs: + for param, value in solver_opts["cplex_params"].items(): + try: + f = attrgetter(param) + parameter = f(model.parameters) + parameter.set(value) + except AttributeError: + raise ValueError( + "invalid CPLEX parameter, value pair ({0}, {1})".format( + param, value)) + kwargs.remove("cplex_params") + if "cplex_filename" in kwargs: + filename = solver_opts["cplex_filename"] + if filename: + model.write(filename) + kwargs.remove("cplex_filename") + if kwargs: + raise ValueError("invalid keyword-argument '{0}'".format(kwargs[0])) + + +def hide_solver_output(model) -> None: + """Set CPLEX verbosity level (either on or off).""" + # By default the output will be sent to stdout. Setting the output + # streams to None will prevent any output from being shown. + model.set_results_stream(None) + model.set_warning_stream(None) + model.set_error_stream(None) + model.set_log_stream(None) + + +def _handle_solve_status(model, solstat): + """Map CPLEX MIP solution status codes to non-MIP status codes.""" + status = model.solution.status + + # With CPLEX 12.10, some benders status codes were removed. + if model.get_versionnumber() < 12100000: + unbounded_status_codes = ( + status.MIP_unbounded, + status.MIP_benders_master_unbounded, + status.benders_master_unbounded + ) + else: + unbounded_status_codes = ( + status.MIP_unbounded, + ) + + if solstat == status.MIP_optimal: + return status.optimal + elif solstat == status.MIP_infeasible: + return status.infeasible + elif solstat in (status.MIP_time_limit_feasible, + status.MIP_time_limit_infeasible): + return status.abort_time_limit + elif solstat in (status.MIP_dettime_limit_feasible, + status.MIP_dettime_limit_infeasible): + return status.abort_dettime_limit + elif solstat in (status.MIP_abort_feasible, + status.MIP_abort_infeasible): + return status.abort_user + elif solstat == status.MIP_optimal_infeasible: + return status.optimal_infeasible + elif solstat == status.MIP_infeasible_or_unbounded: + return status.infeasible_or_unbounded + elif solstat in unbounded_status_codes: + return status.unbounded + elif solstat in (status.feasible_relaxed_sum, + status.MIP_feasible_relaxed_sum, + status.optimal_relaxed_sum, + status.MIP_optimal_relaxed_sum, + status.feasible_relaxed_inf, + status.MIP_feasible_relaxed_inf, + status.optimal_relaxed_inf, + status.MIP_optimal_relaxed_inf, + status.feasible_relaxed_quad, + status.MIP_feasible_relaxed_quad, + status.optimal_relaxed_quad, + status.MIP_optimal_relaxed_quad): + raise AssertionError( + "feasopt status encountered: {0}".format(solstat)) + elif solstat in (status.conflict_feasible, + status.conflict_minimal, + status.conflict_abort_contradiction, + status.conflict_abort_time_limit, + status.conflict_abort_dettime_limit, + status.conflict_abort_iteration_limit, + status.conflict_abort_node_limit, + status.conflict_abort_obj_limit, + status.conflict_abort_memory_limit, + status.conflict_abort_user): + raise AssertionError( + "conflict refiner status encountered: {0}".format(solstat)) + elif solstat == status.relaxation_unbounded: + return status.relaxation_unbounded + elif solstat in (status.feasible, + status.MIP_feasible): + return status.feasible + elif solstat == status.benders_num_best: + return status.num_best + else: + return solstat + + +def get_status(model): + """Map CPLEX status to CPXPY status.""" + pfeas = model.solution.is_primal_feasible() + # NOTE: dfeas is always false for a MIP. + dfeas = model.solution.is_dual_feasible() + status = model.solution.status + solstat = _handle_solve_status(model, model.solution.get_status()) + if solstat in (status.node_limit_infeasible, + status.fail_infeasible, + status.mem_limit_infeasible, + status.fail_infeasible_no_tree, + status.num_best): + return s.SOLVER_ERROR + elif solstat in (status.abort_user, + status.abort_iteration_limit, + status.abort_time_limit, + status.abort_dettime_limit, + status.abort_obj_limit, + status.abort_primal_obj_limit, + status.abort_dual_obj_limit, + status.abort_relaxed, + status.first_order): + if pfeas: + return s.OPTIMAL_INACCURATE + else: + return s.SOLVER_ERROR + elif solstat in (status.node_limit_feasible, + status.solution_limit, + status.populate_solution_limit, + status.fail_feasible, + status.mem_limit_feasible, + status.fail_feasible_no_tree, + status.feasible): + if dfeas: + return s.OPTIMAL + else: + return s.OPTIMAL_INACCURATE + elif solstat in (status.optimal, + status.optimal_tolerance, + status.optimal_infeasible, + status.optimal_populated, + status.optimal_populated_tolerance): + return s.OPTIMAL + elif solstat in (status.infeasible, + status.optimal_relaxed_sum, + status.optimal_relaxed_inf, + status.optimal_relaxed_quad): + return s.INFEASIBLE + elif solstat in (status.feasible_relaxed_quad, + status.feasible_relaxed_inf, + status.feasible_relaxed_sum): + return s.SOLVER_ERROR + elif solstat == status.unbounded: + return s.UNBOUNDED + elif solstat == status.infeasible_or_unbounded: + return s.INFEASIBLE_OR_UNBOUNDED + else: + return s.SOLVER_ERROR + + +class CPLEX(ConicSolver): + """An interface for the CPLEX solver. + """ + + # Solver capabilities. + MIP_CAPABLE = True + SUPPORTED_CONSTRAINTS = ConicSolver.SUPPORTED_CONSTRAINTS + [SOC] + + def name(self): + """The name of the solver. """ + return s.CPLEX + + def import_solver(self) -> None: + """Imports the solver.""" + import cplex + cplex # For flake8 + + def accepts(self, problem) -> bool: + """Can CPLEX solve the problem? + """ + # TODO check if is matrix stuffed. + if not problem.objective.args[0].is_affine(): + return False + for constr in problem.constraints: + if type(constr) not in CPLEX.SUPPORTED_CONSTRAINTS: + return False + for arg in constr.args: + if not arg.is_affine(): + return False + return True + + def apply(self, problem): + """Returns a new problem and data for inverting the new solution. + + Returns + ------- + tuple + (dict of arguments needed for the solver, inverse data) + """ + data, inv_data = super(CPLEX, self).apply(problem) + variables = problem.x + data[s.BOOL_IDX] = [int(t[0]) for t in variables.boolean_idx] + data[s.INT_IDX] = [int(t[0]) for t in variables.integer_idx] + inv_data['is_mip'] = data[s.BOOL_IDX] or data[s.INT_IDX] + + return data, inv_data + + def invert(self, solution, inverse_data): + """Returns the solution to the original problem given the inverse_data. + """ + model = solution["model"] + attr = {} + if s.SOLVE_TIME in solution: + attr[s.SOLVE_TIME] = solution[s.SOLVE_TIME] + + status = get_status(model) + + primal_vars = None + dual_vars = None + if status in s.SOLUTION_PRESENT: + opt_val = (model.solution.get_objective_value() + + inverse_data[s.OFFSET]) + + #x = np.array(model.solution.get_values())#2022-10-20 RWS read .sol file here.XML. + print("reading soln.sol;") + import xml.etree.ElementTree as ET + vd="c:/Users/Rwsin/myproject/" + tree = ET.parse(vd+"soln.sol")#"/mnt/c/Users/Rwsin/myproject/soln.sol") + root = tree.getroot() + x=np.array([child.attrib['value'] for child in root[3]]) + + + primal_vars = {inverse_data[CPLEX.VAR_ID]: x} + + if not inverse_data['is_mip']: + # The dual values are retrieved in the order that the + # constraints were added in solve_via_data() below. We + # must be careful to map them to inverse_data[EQ_CONSTR] + # followed by inverse_data[NEQ_CONSTR] accordingly. + y = -np.array(model.solution.get_dual_values()) + dual_vars = utilities.get_dual_values( + y, + utilities.extract_dual_value, + inverse_data[CPLEX.EQ_CONSTR] + inverse_data[CPLEX.NEQ_CONSTR]) + + return Solution(status, opt_val, primal_vars, dual_vars, attr) + else: + return failure_solution(status) + + def solve_via_data_1(self, data, warm_start: bool, verbose: bool, solver_opts, solver_cache=None):#2022-09-30 + print("cplex solve_via_data_1!") + import cplex + c = data[s.C] + b = data[s.B] + A = dok_matrix(data[s.A]) + # Save the dok_matrix. + data[s.A] = A + dims = dims_to_solver_dict(data[s.DIMS]) + + n = c.shape[0] + + model = cplex.Cplex() + variables = [] + # cpx_constrs will contain CpxConstr namedtuples (see above). + cpx_constrs = [] + vtype = [] + if data[s.BOOL_IDX] or data[s.INT_IDX]: + for i in range(n): + # Set variable type. + if i in data[s.BOOL_IDX]: + vtype.append('B') + elif i in data[s.INT_IDX]: + vtype.append('I') + else: + vtype.append('C') + else: + # If we specify types (even with 'C'), then the problem will + # be interpreted as a MIP. Leaving vtype as an empty list + # here, will ensure that the problem type remains an LP. + pass + # Add the variables in a batch + variables = list(model.variables.add( + obj=[c[i] for i in range(n)], + lb=[-cplex.infinity]*n, # default LB is 0 + ub=[cplex.infinity]*n, + types="".join(vtype), + names=["x_%d" % i for i in range(n)])) + + # Add equality constraints + cpx_constrs += [_CpxConstr(_LIN, x) + for x in self.add_model_lin_constr( + model, variables, + range(dims[s.EQ_DIM]), + 'E', A, b)] + + # Add inequality (<=) constraints + leq_start = dims[s.EQ_DIM] + leq_end = dims[s.EQ_DIM] + dims[s.LEQ_DIM] + cpx_constrs += [_CpxConstr(_LIN, x) + for x in self.add_model_lin_constr( + model, variables, + range(leq_start, leq_end), + 'L', A, b)] + + # Add SOC constraints + soc_start = leq_end + for constr_len in dims[s.SOC_DIM]: + soc_end = soc_start + constr_len + soc_constr, new_leq, new_vars = self.add_model_soc_constr( + model, variables, range(soc_start, soc_end), A, b) + cpx_constrs.append(_CpxConstr(_QUAD, soc_constr)) + cpx_constrs += [_CpxConstr(_LIN, x) for x in new_leq] + variables += new_vars + soc_start += constr_len + + # Set verbosity + if not verbose: + hide_solver_output(model) + + # For CVXPY, we set the qcpduals parameter here, but the user can + # easily override it via the "cplex_params" solver option (see + # set_parameters function). + model.parameters.preprocessing.qcpduals.set( + model.parameters.preprocessing.qcpduals.values.force) + + # Set parameters + reoptimize = solver_opts.pop('reoptimize', False) + set_parameters(model, solver_opts) + + # Solve problem + solution = {"model": model} + #RWS now write_lp() + model.write("svd.mps") + + def solve_via_data_2(self, data, warm_start: bool, verbose: bool, solver_opts, solver_cache=None): + + #model.solve() + #create an empty model + # load solution file + #solution[s.SOLVE_TIME] = model.get_time() - start_time + + # ambiguous_status = get_status(model) == s.INFEASIBLE_OR_UNBOUNDED + # if ambiguous_status and reoptimize: + # model.parameters.preprocessing.presolve.set(0) + # start_time = model.get_time() + # model.solve() + # solution[s.SOLVE_TIME] += model.get_time() - start_time + + # except Exception: + # pass + + # return solution + return + + + + + + def solve_via_data(self, data, warm_start: bool, verbose: bool, solver_opts, solver_cache=None): + import cplex + print("i can write code") + c = data[s.C] + b = data[s.B] + A = dok_matrix(data[s.A]) + # Save the dok_matrix. + data[s.A] = A + dims = dims_to_solver_dict(data[s.DIMS]) + + n = c.shape[0] + + model = cplex.Cplex() + variables = [] + # cpx_constrs will contain CpxConstr namedtuples (see above). + cpx_constrs = [] + vtype = [] + if data[s.BOOL_IDX] or data[s.INT_IDX]: + for i in range(n): + # Set variable type. + if i in data[s.BOOL_IDX]: + vtype.append('B') + elif i in data[s.INT_IDX]: + vtype.append('I') + else: + vtype.append('C') + else: + # If we specify types (even with 'C'), then the problem will + # be interpreted as a MIP. Leaving vtype as an empty list + # here, will ensure that the problem type remains an LP. + pass + # Add the variables in a batch + variables = list(model.variables.add( + obj=[c[i] for i in range(n)], + lb=[-cplex.infinity]*n, # default LB is 0 + ub=[cplex.infinity]*n, + types="".join(vtype), + names=["x_%d" % i for i in range(n)])) + + # Add equality constraints + cpx_constrs += [_CpxConstr(_LIN, x) + for x in self.add_model_lin_constr( + model, variables, + range(dims[s.EQ_DIM]), + 'E', A, b)] + + # Add inequality (<=) constraints + leq_start = dims[s.EQ_DIM] + leq_end = dims[s.EQ_DIM] + dims[s.LEQ_DIM] + cpx_constrs += [_CpxConstr(_LIN, x) + for x in self.add_model_lin_constr( + model, variables, + range(leq_start, leq_end), + 'L', A, b)] + + # Add SOC constraints + soc_start = leq_end + for constr_len in dims[s.SOC_DIM]: + soc_end = soc_start + constr_len + soc_constr, new_leq, new_vars = self.add_model_soc_constr( + model, variables, range(soc_start, soc_end), A, b) + cpx_constrs.append(_CpxConstr(_QUAD, soc_constr)) + cpx_constrs += [_CpxConstr(_LIN, x) for x in new_leq] + variables += new_vars + soc_start += constr_len + + # Set verbosity + if not verbose: + hide_solver_output(model) + + # For CVXPY, we set the qcpduals parameter here, but the user can + # easily override it via the "cplex_params" solver option (see + # set_parameters function). + model.parameters.preprocessing.qcpduals.set( + model.parameters.preprocessing.qcpduals.values.force) + + # Set parameters + reoptimize = solver_opts.pop('reoptimize', False) + set_parameters(model, solver_opts) + + # Solve problem + solution = {"model": model} + try: + start_time = model.get_time() + model.solve() + solution[s.SOLVE_TIME] = model.get_time() - start_time + + ambiguous_status = get_status(model) == s.INFEASIBLE_OR_UNBOUNDED + if ambiguous_status and reoptimize: + model.parameters.preprocessing.presolve.set(0) + start_time = model.get_time() + model.solve() + solution[s.SOLVE_TIME] += model.get_time() - start_time + + except Exception: + pass + + return solution + + def add_model_lin_constr(self, model, variables, + rows, ctype, mat, vec): + """Adds EQ/LEQ constraints to the model using the data from mat and vec. + + Parameters + ---------- + model : CPLEX model + The problem model. + variables : list + The problem variables. + rows : range + The rows to be constrained. + ctype : CPLEX constraint type + The type of constraint. + mat : SciPy COO matrix + The matrix representing the constraints. + vec : NDArray + The RHS part of the constraints. + + Returns + ------- + list + A list of new linear constraint indices. + """ + constr, lin_expr, rhs = [], [], [] + csr = mat.tocsr() + for i in rows: + ind = [variables[x] for x in csr[i].indices] + val = [x for x in csr[i].data] + lin_expr.append([ind, val]) + rhs.append(vec[i]) + # For better performance, we add the constraints in a batch. + if lin_expr: + assert len(lin_expr) == len(rhs) + constr.extend(list( + model.linear_constraints.add( + lin_expr=lin_expr, + senses=ctype * len(lin_expr), + rhs=rhs))) + return constr + + def add_model_soc_constr(self, model, variables, + rows, mat, vec): + """Adds SOC constraint to the model using the data from mat and vec. + + Parameters + ---------- + model : CPLEX model + The problem model. + variables : list + The problem variables. + rows : range + The rows to be constrained. + mat : SciPy COO matrix + The matrix representing the constraints. + vec : NDArray + The RHS part of the constraints. + + Returns + ------- + tuple + A tuple of (a new quadratic constraint index, a list of new + supporting linear constr indices, and a list of new + supporting variable indices). + """ + import cplex + + # Assume first expression (i.e. t) is nonzero. + lin_expr_list, soc_vars, lin_rhs = [], [], [] + csr = mat.tocsr() + for i in rows: + ind = [variables[x] for x in csr[i].indices] + val = [x for x in csr[i].data] + # Ignore empty constraints. + if ind: + lin_expr_list.append((ind, val)) + else: + lin_expr_list.append(None) + lin_rhs.append(vec[i]) + + # Make a variable and equality constraint for each term. + soc_vars, is_first = [], True + for i in rows: + if is_first: + lb = [0.0] + names = ["soc_t_%d" % i] + is_first = False + else: + lb = [-cplex.infinity] + names = ["soc_x_%d" % i] + soc_vars.extend(list(model.variables.add( + obj=[0], + lb=lb, + ub=[cplex.infinity], + types="", + names=names))) + + new_lin_constrs = [] + for i, expr in enumerate(lin_expr_list): + if expr is None: + ind = [soc_vars[i]] + val = [1.0] + else: + ind, val = expr + ind.append(soc_vars[i]) + val.append(1.0) + new_lin_constrs.extend(list( + model.linear_constraints.add( + lin_expr=[cplex.SparsePair(ind=ind, val=val)], + senses="E", + rhs=[lin_rhs[i]]))) + + assert len(soc_vars) > 0 + qconstr = model.quadratic_constraints.add( + lin_expr=cplex.SparsePair(ind=[], val=[]), + quad_expr=cplex.SparseTriple( + ind1=soc_vars, + ind2=soc_vars, + val=[-1.0] + [1.0] * (len(soc_vars) - 1)), + sense="L", + rhs=0.0, + name="") + return (qconstr, new_lin_constrs, soc_vars) From 6be2ae42fc0906a0eb07aa5c32d434a6e51574ca Mon Sep 17 00:00:00 2001 From: rsemenoff Date: Wed, 30 Nov 2022 13:58:22 -0800 Subject: [PATCH 5/7] Add files via upload solve_via_data_1 --- cvxpy/reductions/solvers/solving_chain.py | 850 +++++++++++----------- 1 file changed, 443 insertions(+), 407 deletions(-) diff --git a/cvxpy/reductions/solvers/solving_chain.py b/cvxpy/reductions/solvers/solving_chain.py index 66a9551028..c741f7aca1 100644 --- a/cvxpy/reductions/solvers/solving_chain.py +++ b/cvxpy/reductions/solvers/solving_chain.py @@ -1,407 +1,443 @@ -from __future__ import annotations - -import warnings -from typing import Any, List, Optional - -import numpy as np - -from cvxpy.atoms import EXP_ATOMS, NONPOS_ATOMS, PSD_ATOMS, SOC_ATOMS -from cvxpy.constraints import (PSD, SOC, Equality, ExpCone, FiniteSet, - Inequality, NonNeg, NonPos, PowCone3D, Zero,) -from cvxpy.constraints.exponential import OpRelConeQuad, RelEntrQuad -from cvxpy.error import DCPError, DGPError, DPPError, SolverError -from cvxpy.problems.objective import Maximize -from cvxpy.reductions.chain import Chain -from cvxpy.reductions.complex2real import complex2real -from cvxpy.reductions.cone2cone.approximations import APPROX_CONES, QuadApprox -from cvxpy.reductions.cone2cone.exotic2common import (EXOTIC_CONES, - Exotic2Common,) -from cvxpy.reductions.cvx_attr2constr import CvxAttr2Constr -from cvxpy.reductions.dcp2cone.cone_matrix_stuffing import ConeMatrixStuffing -from cvxpy.reductions.dcp2cone.dcp2cone import Dcp2Cone -from cvxpy.reductions.dgp2dcp.dgp2dcp import Dgp2Dcp -from cvxpy.reductions.discrete2mixedint.valinvec2mixedint import ( - Valinvec2mixedint,) -from cvxpy.reductions.eval_params import EvalParams -from cvxpy.reductions.flip_objective import FlipObjective -from cvxpy.reductions.qp2quad_form import qp2symbolic_qp -from cvxpy.reductions.qp2quad_form.qp_matrix_stuffing import QpMatrixStuffing -from cvxpy.reductions.solvers import defines as slv_def -from cvxpy.reductions.solvers.constant_solver import ConstantSolver -from cvxpy.reductions.solvers.solver import Solver -from cvxpy.settings import PARAM_THRESHOLD -from cvxpy.utilities.debug_tools import build_non_disciplined_error_msg - - -def _is_lp(self): - """Is problem a linear program? - """ - for c in self.constraints: - if not (isinstance(c, (Equality, Zero)) or c.args[0].is_pwl()): - return False - for var in self.variables(): - if var.is_psd() or var.is_nsd(): - return False - return (self.is_dcp() and self.objective.args[0].is_pwl()) - - -def _solve_as_qp(problem, candidates): - if _is_lp(problem) and \ - [s for s in candidates['conic_solvers'] if s not in candidates['qp_solvers']]: - # OSQP can take many iterations for LPs; use a conic solver instead - # GUROBI and CPLEX QP/LP interfaces are more efficient - # -> Use them instead of conic if applicable. - return False - return candidates['qp_solvers'] and qp2symbolic_qp.accepts(problem) - - -def _reductions_for_problem_class(problem, candidates, gp: bool = False) -> List[Any]: - """ - Builds a chain that rewrites a problem into an intermediate - representation suitable for numeric reductions. - - Parameters - ---------- - problem : Problem - The problem for which to build a chain. - candidates : dict - Dictionary of candidate solvers divided in qp_solvers - and conic_solvers. - gp : bool - If True, the problem is parsed as a Disciplined Geometric Program - instead of as a Disciplined Convex Program. - Returns - ------- - list of Reduction objects - A list of reductions that can be used to convert the problem to an - intermediate form. - Raises - ------ - DCPError - Raised if the problem is not DCP and `gp` is False. - DGPError - Raised if the problem is not DGP and `gp` is True. - """ - reductions = [] - # TODO Handle boolean constraints. - if complex2real.accepts(problem): - reductions += [complex2real.Complex2Real()] - if gp: - reductions += [Dgp2Dcp()] - - if not gp and not problem.is_dcp(): - append = build_non_disciplined_error_msg(problem, 'DCP') - if problem.is_dgp(): - append += ("\nHowever, the problem does follow DGP rules. " - "Consider calling solve() with `gp=True`.") - elif problem.is_dqcp(): - append += ("\nHowever, the problem does follow DQCP rules. " - "Consider calling solve() with `qcp=True`.") - raise DCPError( - "Problem does not follow DCP rules. Specifically:\n" + append) - elif gp and not problem.is_dgp(): - append = build_non_disciplined_error_msg(problem, 'DGP') - if problem.is_dcp(): - append += ("\nHowever, the problem does follow DCP rules. " - "Consider calling solve() with `gp=False`.") - elif problem.is_dqcp(): - append += ("\nHowever, the problem does follow DQCP rules. " - "Consider calling solve() with `qcp=True`.") - raise DGPError("Problem does not follow DGP rules." + append) - - # Dcp2Cone and Qp2SymbolicQp require problems to minimize their objectives. - if type(problem.objective) == Maximize: - reductions += [FlipObjective()] - - if _solve_as_qp(problem, candidates): - reductions += [CvxAttr2Constr(), qp2symbolic_qp.Qp2SymbolicQp()] - else: - # Canonicalize it to conic problem. - if not candidates['conic_solvers']: - raise SolverError("Problem could not be reduced to a QP, and no " - "conic solvers exist among candidate solvers " - "(%s)." % candidates) - - constr_types = {type(c) for c in problem.constraints} - if FiniteSet in constr_types: - reductions += [Valinvec2mixedint()] - - return reductions - - -def construct_solving_chain(problem, candidates, - gp: bool = False, - enforce_dpp: bool = False, - ignore_dpp: bool = False, - canon_backend: str | None = None, - solver_opts: Optional[dict] = None - ) -> "SolvingChain": - """Build a reduction chain from a problem to an installed solver. - - Note that if the supplied problem has 0 variables, then the solver - parameter will be ignored. - - Parameters - ---------- - problem : Problem - The problem for which to build a chain. - candidates : dict - Dictionary of candidate solvers divided in qp_solvers - and conic_solvers. - gp : bool - If True, the problem is parsed as a Disciplined Geometric Program - instead of as a Disciplined Convex Program. - enforce_dpp : bool, optional - When True, a DPPError will be thrown when trying to parse a non-DPP - problem (instead of just a warning). Defaults to False. - ignore_dpp : bool, optional - When True, DPP problems will be treated as non-DPP, - which may speed up compilation. Defaults to False. - canon_backend : str, optional - 'CPP' (default) | 'SCIPY' - Specifies which backend to use for canonicalization, which can affect - compilation time. Defaults to None, i.e., selecting the default - backend. - solver_opts : dict, optional - Additional arguments to pass to the solver. - - Returns - ------- - SolvingChain - A SolvingChain that can be used to solve the problem. - - Raises - ------ - SolverError - Raised if no suitable solver exists among the installed solvers, or - if the target solver is not installed. - """ - if len(problem.variables()) == 0: - return SolvingChain(reductions=[ConstantSolver()]) - reductions = _reductions_for_problem_class(problem, candidates, gp) - - # Process DPP status of the problem. - dpp_context = 'dcp' if not gp else 'dgp' - dpp_error_msg = ( - "You are solving a parameterized problem that is not DPP. " - "Because the problem is not DPP, subsequent solves will not be " - "faster than the first one. For more information, see the " - "documentation on Discplined Parametrized Programming, at\n" - "\thttps://www.cvxpy.org/tutorial/advanced/index.html#" - "disciplined-parametrized-programming") - if ignore_dpp or not problem.is_dpp(dpp_context): - # No warning for ignore_dpp. - if ignore_dpp: - reductions = [EvalParams()] + reductions - elif not enforce_dpp: - warnings.warn(dpp_error_msg) - reductions = [EvalParams()] + reductions - else: - raise DPPError(dpp_error_msg) - elif any(param.is_complex() for param in problem.parameters()): - reductions = [EvalParams()] + reductions - else: # Compilation with DPP. - n_parameters = sum(np.prod(param.shape) for param in problem.parameters()) - if n_parameters >= PARAM_THRESHOLD: - warnings.warn( - "Your problem has too many parameters for efficient DPP " - "compilation. We suggest setting 'ignore_dpp = True'." - ) - - # Conclude with matrix stuffing; choose one of the following paths: - # (1) QpMatrixStuffing --> [a QpSolver], - # (2) ConeMatrixStuffing --> [a ConicSolver] - if _solve_as_qp(problem, candidates): - # Canonicalize as a QP - solver = candidates['qp_solvers'][0] - solver_instance = slv_def.SOLVER_MAP_QP[solver] - reductions += [QpMatrixStuffing(canon_backend=canon_backend), - solver_instance] - return SolvingChain(reductions=reductions) - - # Canonicalize as a cone program - if not candidates['conic_solvers']: - raise SolverError("Problem could not be reduced to a QP, and no " - "conic solvers exist among candidate solvers " - "(%s)." % candidates) - - constr_types = set() - # ^ We use constr_types to infer an incomplete list of cones that - # the solver will need after canonicalization. - for c in problem.constraints: - constr_types.add(type(c)) - ex_cos = [ct for ct in constr_types if ct in EXOTIC_CONES] - approx_cos = [ct for ct in constr_types if ct in APPROX_CONES] - # ^ The way we populate "ex_cos" will need to change if and when - # we have atoms that require exotic cones. - for co in ex_cos: - sim_cos = EXOTIC_CONES[co] # get the set of required simple cones - constr_types.update(sim_cos) - constr_types.remove(co) - - for co in approx_cos: - app_cos = APPROX_CONES[co] - constr_types.update(app_cos) - constr_types.remove(co) - # We now go over individual elementary cones support by CVXPY ( - # SOC, ExpCone, NonNeg, Zero, PSD, PowCone3D) and check if - # they've appeared in constr_types or if the problem has an atom - # requiring that cone. - cones = [] - atoms = problem.atoms() - if SOC in constr_types or any(atom in SOC_ATOMS for atom in atoms): - cones.append(SOC) - if ExpCone in constr_types or any(atom in EXP_ATOMS for atom in atoms): - cones.append(ExpCone) - if any(t in constr_types for t in [Inequality, NonPos, NonNeg]) \ - or any(atom in NONPOS_ATOMS for atom in atoms): - cones.append(NonNeg) - if Equality in constr_types or Zero in constr_types: - cones.append(Zero) - if PSD in constr_types \ - or any(atom in PSD_ATOMS for atom in atoms) \ - or any(v.is_psd() or v.is_nsd() for v in problem.variables()): - cones.append(PSD) - if PowCone3D in constr_types: - # if we add in atoms that specifically use the 3D power cone - # (rather than the ND power cone), then we'll need to check - # for those atoms here as well. - cones.append(PowCone3D) - - # Here, we make use of the observation that canonicalization only - # increases the number of constraints in our problem. - has_constr = len(cones) > 0 or len(problem.constraints) > 0 - - for solver in candidates['conic_solvers']: - solver_instance = slv_def.SOLVER_MAP_CONIC[solver] - # Cones supported for MI problems may differ from non MI. - if problem.is_mixed_integer(): - supported_constraints = solver_instance.MI_SUPPORTED_CONSTRAINTS - else: - supported_constraints = solver_instance.SUPPORTED_CONSTRAINTS - if (all(c in supported_constraints for c in cones) - and (has_constr or not solver_instance.REQUIRES_CONSTR)): - if ex_cos: - reductions.append(Exotic2Common()) - if RelEntrQuad in approx_cos or OpRelConeQuad in approx_cos: - reductions.append(QuadApprox()) - - # Should the objective be canonicalized to a quadratic? - if solver_opts is None: - use_quad_obj = True - else: - use_quad_obj = solver_opts.get("use_quad_obj", True) - quad_obj = use_quad_obj and solver_instance.supports_quad_obj() and \ - problem.objective.expr.has_quadratic_term() - reductions += [ - Dcp2Cone(quad_obj=quad_obj), - CvxAttr2Constr(), - ConeMatrixStuffing(quad_obj=quad_obj, canon_backend=canon_backend), - solver_instance - ] - return SolvingChain(reductions=reductions) - - raise SolverError("Either candidate conic solvers (%s) do not support the " - "cones output by the problem (%s), or there are not " - "enough constraints in the problem." % ( - candidates['conic_solvers'], - ", ".join([cone.__name__ for cone in cones]))) - - -class SolvingChain(Chain): - """A reduction chain that ends with a solver. - - Parameters - ---------- - reductions : list[Reduction] - A list of reductions. The last reduction in the list must be a solver - instance. - - Attributes - ---------- - reductions : list[Reduction] - A list of reductions. - solver : Solver - The solver, i.e., reductions[-1]. - """ - - def __init__(self, problem=None, reductions=None) -> None: - super(SolvingChain, self).__init__(problem=problem, - reductions=reductions) - if not isinstance(self.reductions[-1], Solver): - raise ValueError("Solving chains must terminate with a Solver.") - self.solver = self.reductions[-1] - - def prepend(self, chain) -> "SolvingChain": - """ - Create and return a new SolvingChain by concatenating - chain with this instance. - """ - return SolvingChain(reductions=chain.reductions + self.reductions) - - def solve(self, problem, warm_start: bool, verbose: bool, solver_opts): - """Solves the problem by applying the chain. - - Applies each reduction in the chain to the problem, solves it, - and then inverts the chain to return a solution of the supplied - problem. - - Parameters - ---------- - problem : Problem - The problem to solve. - warm_start : bool - Whether to warm start the solver. - verbose : bool - Whether to enable solver verbosity. - solver_opts : dict - Solver specific options. - - Returns - ------- - solution : Solution - A solution to the problem. - """ - data, inverse_data = self.apply(problem) - solution = self.solver.solve_via_data(data, warm_start, - verbose, solver_opts) - return self.invert(solution, inverse_data) - - def solve_via_data(self, problem, data, warm_start: bool = False, verbose: bool = False, - solver_opts={}): - """Solves the problem using the data output by the an apply invocation. - - The semantics are: - - .. code :: python - - data, inverse_data = solving_chain.apply(problem) - solution = solving_chain.invert(solver_chain.solve_via_data(data, ...)) - - which is equivalent to writing - - .. code :: python - - solution = solving_chain.solve(problem, ...) - - Parameters - ---------- - problem : Problem - The problem to solve. - data : map - Data for the solver. - warm_start : bool - Whether to warm start the solver. - verbose : bool - Whether to enable solver verbosity. - solver_opts : dict - Solver specific options. - - Returns - ------- - raw solver solution - The information returned by the solver; this is not necessarily - a Solution object. - """ - return self.solver.solve_via_data(data, warm_start, verbose, - solver_opts, problem._solver_cache) +import warnings +from typing import Any, List + +import numpy as np + +from cvxpy.atoms import EXP_ATOMS, NONPOS_ATOMS, PSD_ATOMS, SOC_ATOMS +from cvxpy.constraints import (PSD, SOC, Equality, ExpCone, Inequality, NonNeg, + NonPos, PowCone3D, Zero,) +from cvxpy.error import DCPError, DGPError, DPPError, SolverError +from cvxpy.problems.objective import Maximize +from cvxpy.reductions.chain import Chain +from cvxpy.reductions.complex2real import complex2real +from cvxpy.reductions.cone2cone.exotic2common import (EXOTIC_CONES, + Exotic2Common,) +from cvxpy.reductions.cvx_attr2constr import CvxAttr2Constr +from cvxpy.reductions.dcp2cone.cone_matrix_stuffing import ConeMatrixStuffing +from cvxpy.reductions.dcp2cone.dcp2cone import Dcp2Cone +from cvxpy.reductions.dgp2dcp.dgp2dcp import Dgp2Dcp +from cvxpy.reductions.eval_params import EvalParams +from cvxpy.reductions.flip_objective import FlipObjective +from cvxpy.reductions.qp2quad_form import qp2symbolic_qp +from cvxpy.reductions.qp2quad_form.qp_matrix_stuffing import QpMatrixStuffing +from cvxpy.reductions.solvers import defines as slv_def +from cvxpy.reductions.solvers.constant_solver import ConstantSolver +from cvxpy.reductions.solvers.solver import Solver +from cvxpy.settings import PARAM_THRESHOLD +from cvxpy.utilities.debug_tools import build_non_disciplined_error_msg + + +def _is_lp(self): + """Is problem a linear program? + """ + for c in self.constraints: + if not (isinstance(c, (Equality, Zero)) or c.args[0].is_pwl()): + return False + for var in self.variables(): + if var.is_psd() or var.is_nsd(): + return False + return (self.is_dcp() and self.objective.args[0].is_pwl()) + + +def _solve_as_qp(problem, candidates): + if _is_lp(problem) and \ + [s for s in candidates['conic_solvers'] if s not in candidates['qp_solvers']]: + # OSQP can take many iterations for LPs; use a conic solver instead + # GUROBI and CPLEX QP/LP interfaces are more efficient + # -> Use them instead of conic if applicable. + return False + return candidates['qp_solvers'] and qp2symbolic_qp.accepts(problem) + + +def _reductions_for_problem_class(problem, candidates, gp: bool = False) -> List[Any]: + """ + Builds a chain that rewrites a problem into an intermediate + representation suitable for numeric reductions. + + Parameters + ---------- + problem : Problem + The problem for which to build a chain. + candidates : dict + Dictionary of candidate solvers divided in qp_solvers + and conic_solvers. + gp : bool + If True, the problem is parsed as a Disciplined Geometric Program + instead of as a Disciplined Convex Program. + Returns + ------- + list of Reduction objects + A list of reductions that can be used to convert the problem to an + intermediate form. + Raises + ------ + DCPError + Raised if the problem is not DCP and `gp` is False. + DGPError + Raised if the problem is not DGP and `gp` is True. + """ + reductions = [] + # TODO Handle boolean constraints. + if complex2real.accepts(problem): + reductions += [complex2real.Complex2Real()] + if gp: + reductions += [Dgp2Dcp()] + + if not gp and not problem.is_dcp(): + append = build_non_disciplined_error_msg(problem, 'DCP') + if problem.is_dgp(): + append += ("\nHowever, the problem does follow DGP rules. " + "Consider calling solve() with `gp=True`.") + elif problem.is_dqcp(): + append += ("\nHowever, the problem does follow DQCP rules. " + "Consider calling solve() with `qcp=True`.") + raise DCPError( + "Problem does not follow DCP rules. Specifically:\n" + append) + elif gp and not problem.is_dgp(): + append = build_non_disciplined_error_msg(problem, 'DGP') + if problem.is_dcp(): + append += ("\nHowever, the problem does follow DCP rules. " + "Consider calling solve() with `gp=False`.") + elif problem.is_dqcp(): + append += ("\nHowever, the problem does follow DQCP rules. " + "Consider calling solve() with `qcp=True`.") + raise DGPError("Problem does not follow DGP rules." + append) + + # Dcp2Cone and Qp2SymbolicQp require problems to minimize their objectives. + if type(problem.objective) == Maximize: + reductions += [FlipObjective()] + + if _solve_as_qp(problem, candidates): + reductions += [CvxAttr2Constr(), qp2symbolic_qp.Qp2SymbolicQp()] + else: + # Canonicalize it to conic problem. + if not candidates['conic_solvers']: + raise SolverError("Problem could not be reduced to a QP, and no " + "conic solvers exist among candidate solvers " + "(%s)." % candidates) + else: + reductions += [Dcp2Cone(), CvxAttr2Constr()] + return reductions + + +def construct_solving_chain(problem, candidates, + gp: bool = False, + enforce_dpp: bool = False, + ignore_dpp: bool = False) -> "SolvingChain": + """Build a reduction chain from a problem to an installed solver. + + Note that if the supplied problem has 0 variables, then the solver + parameter will be ignored. + + Parameters + ---------- + problem : Problem + The problem for which to build a chain. + candidates : dict + Dictionary of candidate solvers divided in qp_solvers + and conic_solvers. + gp : bool + If True, the problem is parsed as a Disciplined Geometric Program + instead of as a Disciplined Convex Program. + enforce_dpp : bool, optional + When True, a DPPError will be thrown when trying to parse a non-DPP + problem (instead of just a warning). Defaults to False. + ignore_dpp : bool, optional + When True, DPP problems will be treated as non-DPP, + which may speed up compilation. Defaults to False. + + Returns + ------- + SolvingChain + A SolvingChain that can be used to solve the problem. + + Raises + ------ + SolverError + Raised if no suitable solver exists among the installed solvers, or + if the target solver is not installed. + """ + if len(problem.variables()) == 0: + return SolvingChain(reductions=[ConstantSolver()]) + reductions = _reductions_for_problem_class(problem, candidates, gp) + + # Process DPP status of the problem. + dpp_context = 'dcp' if not gp else 'dgp' + dpp_error_msg = ( + "You are solving a parameterized problem that is not DPP. " + "Because the problem is not DPP, subsequent solves will not be " + "faster than the first one. For more information, see the " + "documentation on Discplined Parametrized Programming, at\n" + "\thttps://www.cvxpy.org/tutorial/advanced/index.html#" + "disciplined-parametrized-programming") + if ignore_dpp or not problem.is_dpp(dpp_context): + # No warning for ignore_dpp. + if ignore_dpp: + reductions = [EvalParams()] + reductions + elif not enforce_dpp: + warnings.warn(dpp_error_msg) + reductions = [EvalParams()] + reductions + else: + raise DPPError(dpp_error_msg) + elif any(param.is_complex() for param in problem.parameters()): + reductions = [EvalParams()] + reductions + else: # Compilation with DPP. + n_parameters = sum(np.prod(param.shape) for param in problem.parameters()) + if n_parameters >= PARAM_THRESHOLD: + warnings.warn( + "Your problem has too many parameters for efficient DPP " + "compilation. We suggest setting 'ignore_dpp = True'." + ) + + # Conclude with matrix stuffing; choose one of the following paths: + # (1) QpMatrixStuffing --> [a QpSolver], + # (2) ConeMatrixStuffing --> [a ConicSolver] + if _solve_as_qp(problem, candidates): + # Canonicalize as a QP + solver = candidates['qp_solvers'][0] + solver_instance = slv_def.SOLVER_MAP_QP[solver] + reductions += [QpMatrixStuffing(), + solver_instance] + return SolvingChain(reductions=reductions) + + # Canonicalize as a cone program + if not candidates['conic_solvers']: + raise SolverError("Problem could not be reduced to a QP, and no " + "conic solvers exist among candidate solvers " + "(%s)." % candidates) + + constr_types = set() + # ^ We use constr_types to infer an incomplete list of cones that + # the solver will need after canonicalization. + for c in problem.constraints: + constr_types.add(type(c)) + ex_cos = [ct for ct in constr_types if ct in EXOTIC_CONES] + # ^ The way we populate "ex_cos" will need to change if and when + # we have atoms that require exotic cones. + for co in ex_cos: + sim_cos = EXOTIC_CONES[co] # get the set of required simple cones + constr_types.update(sim_cos) + constr_types.remove(co) + # We now go over individual elementary cones support by CVXPY ( + # SOC, ExpCone, NonNeg, Zero, PSD, PowCone3D) and check if + # they've appeared in constr_types or if the problem has an atom + # requiring that cone. + cones = [] + atoms = problem.atoms() + if SOC in constr_types or any(atom in SOC_ATOMS for atom in atoms): + cones.append(SOC) + if ExpCone in constr_types or any(atom in EXP_ATOMS for atom in atoms): + cones.append(ExpCone) + if any(t in constr_types for t in [Inequality, NonPos, NonNeg]) \ + or any(atom in NONPOS_ATOMS for atom in atoms): + cones.append(NonNeg) + if Equality in constr_types or Zero in constr_types: + cones.append(Zero) + if PSD in constr_types \ + or any(atom in PSD_ATOMS for atom in atoms) \ + or any(v.is_psd() or v.is_nsd() for v in problem.variables()): + cones.append(PSD) + if PowCone3D in constr_types: + # if we add in atoms that specifically use the 3D power cone + # (rather than the ND power cone), then we'll need to check + # for those atoms here as well. + cones.append(PowCone3D) + + # Here, we make use of the observation that canonicalization only + # increases the number of constraints in our problem. + has_constr = len(cones) > 0 or len(problem.constraints) > 0 + + for solver in candidates['conic_solvers']: + solver_instance = slv_def.SOLVER_MAP_CONIC[solver] + if (all(c in solver_instance.SUPPORTED_CONSTRAINTS for c in cones) + and (has_constr or not solver_instance.REQUIRES_CONSTR)): + if ex_cos: + reductions.append(Exotic2Common()) + reductions += [ConeMatrixStuffing(), solver_instance] + return SolvingChain(reductions=reductions) + + raise SolverError("Either candidate conic solvers (%s) do not support the " + "cones output by the problem (%s), or there are not " + "enough constraints in the problem." % ( + candidates['conic_solvers'], + ", ".join([cone.__name__ for cone in cones]))) + + +class SolvingChain(Chain): + """A reduction chain that ends with a solver. + + Parameters + ---------- + reductions : list[Reduction] + A list of reductions. The last reduction in the list must be a solver + instance. + + Attributes + ---------- + reductions : list[Reduction] + A list of reductions. + solver : Solver + The solver, i.e., reductions[-1]. + """ + + def __init__(self, problem=None, reductions=None) -> None: + super(SolvingChain, self).__init__(problem=problem, + reductions=reductions) + if not isinstance(self.reductions[-1], Solver): + raise ValueError("Solving chains must terminate with a Solver.") + self.solver = self.reductions[-1] + + def prepend(self, chain) -> "SolvingChain": + """ + Create and return a new SolvingChain by concatenating + chain with this instance. + """ + return SolvingChain(reductions=chain.reductions + self.reductions) + + def solve(self, problem, warm_start: bool, verbose: bool, solver_opts): + """Solves the problem by applying the chain. + + Applies each reduction in the chain to the problem, solves it, + and then inverts the chain to return a solution of the supplied + problem. + + Parameters + ---------- + problem : Problem + The problem to solve. + warm_start : bool + Whether to warm start the solver. + verbose : bool + Whether to enable solver verbosity. + solver_opts : dict + Solver specific options. + + Returns + ------- + solution : Solution + A solution to the problem. + """ + data, inverse_data = self.apply(problem) + solution = self.solver.solve_via_data(data, warm_start, + verbose, solver_opts) + return self.invert(solution, inverse_data) + + def solve_via_data(self, problem, data, warm_start: bool = False, verbose: bool = False, + solver_opts={}): + """Solves the problem using the data output by the an apply invocation. + + The semantics are: + + .. code :: python + + data, inverse_data = solving_chain.apply(problem) + solution = solving_chain.invert(solver_chain.solve_via_data(data, ...)) + + which is equivalent to writing + + .. code :: python + + solution = solving_chain.solve(problem, ...) + + Parameters + ---------- + problem : Problem + The problem to solve. + data : map + Data for the solver. + warm_start : bool + Whether to warm start the solver. + verbose : bool + Whether to enable solver verbosity. + solver_opts : dict + Solver specific options. + + Returns + ------- + raw solver solution + The information returned by the solver; this is not necessarily + a Solution object. + """ + return self.solver.solve_via_data(data, warm_start, verbose, + solver_opts, problem._solver_cache) + def solve_via_data_1(self, problem, data, warm_start: bool = False, verbose: bool = False, + solver_opts={}): + """Solves the problem using the data output by the an apply invocation. + + The semantics are: + + .. code :: python + + data, inverse_data = solving_chain.apply(problem) + solution = solving_chain.invert(solver_chain.solve_via_data(data, ...)) + + which is equivalent to writing + + .. code :: python + + solution = solving_chain.solve(problem, ...) + + Parameters + ---------- + problem : Problem + The problem to solve. + data : map + Data for the solver. + warm_start : bool + Whether to warm start the solver. + verbose : bool + Whether to enable solver verbosity. + solver_opts : dict + Solver specific options. + + Returns + ------- + raw solver solution + The information returned by the solver; this is not necessarily + a Solution object. + """ + print("here i am.")#2022-10-10 next line to call via_data_1 + r=self.solver.solve_via_data_1(data, warm_start, verbose, + solver_opts, problem._solver_cache) + return r + + def solve_via_data_2(self, problem, data, warm_start: bool = False, verbose: bool = False, + solver_opts={}): + """Solves the problem using the data output by the an apply invocation. + + The semantics are: + + .. code :: python + + data, inverse_data = solving_chain.apply(problem) + solution = solving_chain.invert(solver_chain.solve_via_data(data, ...)) + + which is equivalent to writing + + .. code :: python + + solution = solving_chain.solve(problem, ...) + + Parameters + ---------- + problem : Problem + The problem to solve. + data : map + Data for the solver. + warm_start : bool + Whether to warm start the solver. + verbose : bool + Whether to enable solver verbosity. + solver_opts : dict + Solver specific options. + + Returns + ------- + raw solver solution + The information returned by the solver; this is not necessarily + a Solution object. + """ + print("here i am.") + return self.solver.solve_via_data_2(data, warm_start, verbose, + solver_opts, problem._solver_cache) + From 4d55e1995cb57592438aa98ee7dd11c85085467b Mon Sep 17 00:00:00 2001 From: rsemenoff Date: Thu, 1 Dec 2022 12:43:01 -0800 Subject: [PATCH 6/7] Add files via upload --- .../solvers/qp_solvers/cplex_qpif.py | 459 +++++++++++------- 1 file changed, 288 insertions(+), 171 deletions(-) diff --git a/cvxpy/reductions/solvers/qp_solvers/cplex_qpif.py b/cvxpy/reductions/solvers/qp_solvers/cplex_qpif.py index 7956c42531..1cafcbdf4c 100644 --- a/cvxpy/reductions/solvers/qp_solvers/cplex_qpif.py +++ b/cvxpy/reductions/solvers/qp_solvers/cplex_qpif.py @@ -1,171 +1,288 @@ -import numpy as np - -import cvxpy.interface as intf -import cvxpy.settings as s -from cvxpy.reductions.solution import Solution, failure_solution -from cvxpy.reductions.solvers.conic_solvers.cplex_conif import ( - get_status, hide_solver_output, set_parameters,) -from cvxpy.reductions.solvers.qp_solvers.qp_solver import QpSolver - - -def constrain_cplex_infty(v) -> None: - ''' - Limit values of vector v between +/- infinity as - defined in the CPLEX package - ''' - import cplex as cpx - n = len(v) - - for i in range(n): - if v[i] >= cpx.infinity: - v[i] = cpx.infinity - if v[i] <= -cpx.infinity: - v[i] = -cpx.infinity - - -class CPLEX(QpSolver): - """QP interface for the CPLEX solver""" - - MIP_CAPABLE = True - - def name(self): - return s.CPLEX - - def import_solver(self) -> None: - import cplex - cplex - - def invert(self, results, inverse_data): - model = results["model"] - attr = {} - if "cputime" in results: - attr[s.SOLVE_TIME] = results["cputime"] - attr[s.NUM_ITERS] = \ - int(model.solution.progress.get_num_barrier_iterations()) \ - if not inverse_data[CPLEX.IS_MIP] \ - else 0 - - status = get_status(model) - - if status in s.SOLUTION_PRESENT: - # Get objective value - opt_val = model.solution.get_objective_value() + \ - inverse_data[s.OFFSET] - - # Get solution - x = np.array(model.solution.get_values()) - primal_vars = { - CPLEX.VAR_ID: - intf.DEFAULT_INTF.const_to_matrix(np.array(x)) - } - - # Only add duals if not a MIP. - dual_vars = None - if not inverse_data[CPLEX.IS_MIP]: - y = -np.array(model.solution.get_dual_values()) - dual_vars = {CPLEX.DUAL_VAR_ID: y} - - sol = Solution(status, opt_val, primal_vars, dual_vars, attr) - else: - sol = failure_solution(status, attr) - return sol - - def solve_via_data(self, data, warm_start: bool, verbose: bool, solver_opts, solver_cache=None): - import cplex as cpx - P = data[s.P].tocsr() # Convert matrix to csr format - q = data[s.Q] - A = data[s.A].tocsr() # Convert A matrix to csr format - b = data[s.B] - F = data[s.F].tocsr() # Convert F matrix to csr format - g = data[s.G] - n_var = data['n_var'] - n_eq = data['n_eq'] - n_ineq = data['n_ineq'] - - # Constrain values between bounds - constrain_cplex_infty(b) - constrain_cplex_infty(g) - - # Define CPLEX problem - model = cpx.Cplex() - - # Minimize problem - model.objective.set_sense(model.objective.sense.minimize) - - # Add variables and linear objective - var_idx = list(model.variables.add(obj=q, - lb=-cpx.infinity*np.ones(n_var), - ub=cpx.infinity*np.ones(n_var))) - - # Constrain binary/integer variables if present - for i in data[s.BOOL_IDX]: - model.variables.set_types(var_idx[i], - model.variables.type.binary) - for i in data[s.INT_IDX]: - model.variables.set_types(var_idx[i], - model.variables.type.integer) - - # Add constraints - lin_expr, rhs = [], [] - for i in range(n_eq): # Add equalities - start = A.indptr[i] - end = A.indptr[i+1] - lin_expr.append([A.indices[start:end].tolist(), - A.data[start:end].tolist()]) - rhs.append(b[i]) - if lin_expr: - model.linear_constraints.add(lin_expr=lin_expr, - senses=["E"] * len(lin_expr), - rhs=rhs) - - lin_expr, rhs = [], [] - for i in range(n_ineq): # Add inequalities - start = F.indptr[i] - end = F.indptr[i+1] - lin_expr.append([F.indices[start:end].tolist(), - F.data[start:end].tolist()]) - rhs.append(g[i]) - if lin_expr: - model.linear_constraints.add(lin_expr=lin_expr, - senses=["L"] * len(lin_expr), - rhs=rhs) - - # Set quadratic Cost - if P.count_nonzero(): # Only if quadratic form is not null - qmat = [] - for i in range(n_var): - start = P.indptr[i] - end = P.indptr[i+1] - qmat.append([P.indices[start:end].tolist(), - P.data[start:end].tolist()]) - model.objective.set_quadratic(qmat) - - # Set verbosity - if not verbose: - hide_solver_output(model) - - # Set parameters - reoptimize = solver_opts.pop('reoptimize', False) - set_parameters(model, solver_opts) - - # Solve problem - results_dict = {} - try: - start = model.get_time() - model.solve() - end = model.get_time() - results_dict["cputime"] = end - start - - ambiguous_status = get_status(model) == s.INFEASIBLE_OR_UNBOUNDED - if ambiguous_status and reoptimize: - model.parameters.preprocessing.presolve.set(0) - start_time = model.get_time() - model.solve() - results_dict["cputime"] += model.get_time() - start_time - - except Exception: # Error in the solution - results_dict["status"] = s.SOLVER_ERROR - - results_dict["model"] = model - - return results_dict +import numpy as np + +import cvxpy.interface as intf +import cvxpy.settings as s +from cvxpy.reductions.solution import Solution, failure_solution +from cvxpy.reductions.solvers.conic_solvers.cplex_conif import ( + get_status, hide_solver_output, set_parameters,) +from cvxpy.reductions.solvers.qp_solvers.qp_solver import QpSolver + + +def constrain_cplex_infty(v) -> None: + ''' + Limit values of vector v between +/- infinity as + defined in the CPLEX package + ''' + import cplex as cpx + n = len(v) + + for i in range(n): + if v[i] >= cpx.infinity: + v[i] = cpx.infinity + if v[i] <= -cpx.infinity: + v[i] = -cpx.infinity + + +class CPLEX(QpSolver): + """QP interface for the CPLEX solver""" + + MIP_CAPABLE = True + + def name(self): + return s.CPLEX + + def import_solver(self) -> None: + import cplex + cplex + + def invert(self, results, inverse_data): + print("in cplex_qpif;") + + model = results["model"] + attr = {} + if "cputime" in results: + attr[s.SOLVE_TIME] = results["cputime"] + #this next lines throws even if problem is a MIP...2022-10-26 + attr[s.NUM_ITERS]=0 + # attr[s.NUM_ITERS] = \ + # int(model.solution.progress.get_num_barrier_iterations()) \ + # if not inverse_data[CPLEX.IS_MIP] \ + # else 0 + + #status = get_status(model)#reports solver error....2022-10-26 + status=s.OPTIMAL#SOLUTION_PRESENT#2022-10-26 + + if status in s.SOLUTION_PRESENT: + print("status ok.") + if True: + # Get objective value + # opt_val = model.solution.get_objective_value() + \ + # inverse_data[s.OFFSET] + + # Get solution + #x = np.array(model.solution.get_values())#2022-10-26 RWS read .sol file here.XML. + vd="./"#c:/Users/Rwsin/myproject/" + #2022-11-01 if reading CPLEX neos output : + print("cplex_qpif reading soln.sol;") + + # import xml.etree.ElementTree as ET #2022-11-27 works fine but gurobi preferred (below) + # tree = ET.parse(vd+"soln.sol")#"/mnt/c/Users/Rwsin/myproject/soln.sol") + # root = tree.getroot() + # x=np.array([child.attrib['value'] for child in root[3]]) + + #2022-11-01 if reading gurobi cplex output: + import pandas as pd + cbcsol=pd.read_table(vd+"model.sol",index_col=None,sep=None,names=['NAME','CVX_xpress_qp'],dtype={'NAME':str,'CVX_xpress_qp':float},skiprows=1,engine='python') + x=cbcsol.CVX_xpress_qp.values + + primal_vars = { + CPLEX.VAR_ID: + intf.DEFAULT_INTF.const_to_matrix(np.array(x)) + } + + # Only add duals if not a MIP. + dual_vars = None + if not inverse_data[CPLEX.IS_MIP]: + y = -np.array(model.solution.get_dual_values())#2022-10-26 throws error + dual_vars = {CPLEX.DUAL_VAR_ID: y} + + #sol = Solution(status, opt_val, primal_vars, dual_vars, attr) + sol = Solution(status, 0.0, primal_vars, dual_vars, attr)#20222-10-26 + + else: + sol = failure_solution(status, attr) + return sol + + def solve_via_data_1(self, data, warm_start: bool, verbose: bool, solver_opts, solver_cache=None): + print("in cplex_qpif") + import cplex as cpx + P = data[s.P].tocsr() # Convert matrix to csr format + q = data[s.Q] + A = data[s.A].tocsr() # Convert A matrix to csr format + b = data[s.B] + F = data[s.F].tocsr() # Convert F matrix to csr format + g = data[s.G] + n_var = data['n_var'] + n_eq = data['n_eq'] + n_ineq = data['n_ineq'] + + # Constrain values between bounds + constrain_cplex_infty(b) + constrain_cplex_infty(g) + + # Define CPLEX problem + model = cpx.Cplex()#2022-10-26 added self + + # Minimize problem + model.objective.set_sense(model.objective.sense.minimize) + + # Add variables and linear objective + var_idx = list(model.variables.add(obj=q, + lb=-cpx.infinity*np.ones(n_var), + ub=cpx.infinity*np.ones(n_var))) + + # Constrain binary/integer variables if present + for i in data[s.BOOL_IDX]: + model.variables.set_types(var_idx[i], + model.variables.type.binary) + for i in data[s.INT_IDX]: + model.variables.set_types(var_idx[i], + model.variables.type.integer) + + # Add constraints + lin_expr, rhs = [], [] + for i in range(n_eq): # Add equalities + start = A.indptr[i] + end = A.indptr[i+1] + lin_expr.append([A.indices[start:end].tolist(), + A.data[start:end].tolist()]) + rhs.append(b[i]) + if lin_expr: + model.linear_constraints.add(lin_expr=lin_expr, + senses=["E"] * len(lin_expr), + rhs=rhs) + + lin_expr, rhs = [], [] + for i in range(n_ineq): # Add inequalities + start = F.indptr[i] + end = F.indptr[i+1] + lin_expr.append([F.indices[start:end].tolist(), + F.data[start:end].tolist()]) + rhs.append(g[i]) + if lin_expr: + model.linear_constraints.add(lin_expr=lin_expr, + senses=["L"] * len(lin_expr), + rhs=rhs) + + # Set quadratic Cost + if P.count_nonzero(): # Only if quadratic form is not null + qmat = [] + for i in range(n_var): + start = P.indptr[i] + end = P.indptr[i+1] + qmat.append([P.indices[start:end].tolist(), + P.data[start:end].tolist()]) + model.objective.set_quadratic(qmat) + + # Set verbosity + if not verbose: + hide_solver_output(model) + + # Set parameters + reoptimize = solver_opts.pop('reoptimize', False) + set_parameters(model, solver_opts) + + # Solve problem + results_dict = {} + + model.write("cplex.mps") + self.mymodel=model + + def solve_via_data_2(self, data, warm_start: bool, verbose: bool, solver_opts, solver_cache=None): + #bunch of other solve status that might possible be extractable from sol file ... + print("cplex_qpif via_data2.") + results_dict={} + results_dict["model"] = self.mymodel + + return results_dict + + + def solve_via_data(self, data, warm_start: bool, verbose: bool, solver_opts, solver_cache=None): + import cplex as cpx + P = data[s.P].tocsr() # Convert matrix to csr format + q = data[s.Q] + A = data[s.A].tocsr() # Convert A matrix to csr format + b = data[s.B] + F = data[s.F].tocsr() # Convert F matrix to csr format + g = data[s.G] + n_var = data['n_var'] + n_eq = data['n_eq'] + n_ineq = data['n_ineq'] + + # Constrain values between bounds + constrain_cplex_infty(b) + constrain_cplex_infty(g) + + # Define CPLEX problem + model = cpx.Cplex() + + # Minimize problem + model.objective.set_sense(model.objective.sense.minimize) + + # Add variables and linear objective + var_idx = list(model.variables.add(obj=q, + lb=-cpx.infinity*np.ones(n_var), + ub=cpx.infinity*np.ones(n_var))) + + # Constrain binary/integer variables if present + for i in data[s.BOOL_IDX]: + model.variables.set_types(var_idx[i], + model.variables.type.binary) + for i in data[s.INT_IDX]: + model.variables.set_types(var_idx[i], + model.variables.type.integer) + + # Add constraints + lin_expr, rhs = [], [] + for i in range(n_eq): # Add equalities + start = A.indptr[i] + end = A.indptr[i+1] + lin_expr.append([A.indices[start:end].tolist(), + A.data[start:end].tolist()]) + rhs.append(b[i]) + if lin_expr: + model.linear_constraints.add(lin_expr=lin_expr, + senses=["E"] * len(lin_expr), + rhs=rhs) + + lin_expr, rhs = [], [] + for i in range(n_ineq): # Add inequalities + start = F.indptr[i] + end = F.indptr[i+1] + lin_expr.append([F.indices[start:end].tolist(), + F.data[start:end].tolist()]) + rhs.append(g[i]) + if lin_expr: + model.linear_constraints.add(lin_expr=lin_expr, + senses=["L"] * len(lin_expr), + rhs=rhs) + + # Set quadratic Cost + if P.count_nonzero(): # Only if quadratic form is not null + qmat = [] + for i in range(n_var): + start = P.indptr[i] + end = P.indptr[i+1] + qmat.append([P.indices[start:end].tolist(), + P.data[start:end].tolist()]) + model.objective.set_quadratic(qmat) + + # Set verbosity + if not verbose: + hide_solver_output(model) + + # Set parameters + reoptimize = solver_opts.pop('reoptimize', False) + set_parameters(model, solver_opts) + + # Solve problem + results_dict = {} + try: + start = model.get_time() + model.solve() + end = model.get_time() + results_dict["cputime"] = end - start + + ambiguous_status = get_status(model) == s.INFEASIBLE_OR_UNBOUNDED + if ambiguous_status and reoptimize: + model.parameters.preprocessing.presolve.set(0) + start_time = model.get_time() + model.solve() + results_dict["cputime"] += model.get_time() - start_time + + except Exception: # Error in the solution + results_dict["status"] = s.SOLVER_ERROR + + results_dict["model"] = model + + return results_dict From 132d23b159c6d27a2fe599bb908b5869a24d2a87 Mon Sep 17 00:00:00 2001 From: rsemenoff Date: Thu, 1 Dec 2022 15:16:28 -0800 Subject: [PATCH 7/7] Update cplex_qpif.py can edit in github --- cvxpy/reductions/solvers/qp_solvers/cplex_qpif.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cvxpy/reductions/solvers/qp_solvers/cplex_qpif.py b/cvxpy/reductions/solvers/qp_solvers/cplex_qpif.py index 1cafcbdf4c..45fbed5e83 100644 --- a/cvxpy/reductions/solvers/qp_solvers/cplex_qpif.py +++ b/cvxpy/reductions/solvers/qp_solvers/cplex_qpif.py @@ -53,7 +53,7 @@ def invert(self, results, inverse_data): status=s.OPTIMAL#SOLUTION_PRESENT#2022-10-26 if status in s.SOLUTION_PRESENT: - print("status ok.") + print("cplex_qpif.invert status ok.") if True: # Get objective value # opt_val = model.solution.get_objective_value() + \