Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/optimize/_linprog_highs.py: 15%
73 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""HiGHS Linear Optimization Methods
3Interface to HiGHS linear optimization software.
4https://highs.dev/
6.. versionadded:: 1.5.0
8References
9----------
10.. [1] Q. Huangfu and J.A.J. Hall. "Parallelizing the dual revised simplex
11 method." Mathematical Programming Computation, 10 (1), 119-142,
12 2018. DOI: 10.1007/s12532-017-0130-5
14"""
16import inspect
17import numpy as np
18from ._optimize import _check_unknown_options, OptimizeWarning, OptimizeResult
19from warnings import warn
20from ._highs._highs_wrapper import _highs_wrapper
21from ._highs._highs_constants import (
22 CONST_I_INF,
23 CONST_INF,
24 MESSAGE_LEVEL_NONE,
25 HIGHS_OBJECTIVE_SENSE_MINIMIZE,
27 MODEL_STATUS_NOTSET,
28 MODEL_STATUS_LOAD_ERROR,
29 MODEL_STATUS_MODEL_ERROR,
30 MODEL_STATUS_PRESOLVE_ERROR,
31 MODEL_STATUS_SOLVE_ERROR,
32 MODEL_STATUS_POSTSOLVE_ERROR,
33 MODEL_STATUS_MODEL_EMPTY,
34 MODEL_STATUS_OPTIMAL,
35 MODEL_STATUS_INFEASIBLE,
36 MODEL_STATUS_UNBOUNDED_OR_INFEASIBLE,
37 MODEL_STATUS_UNBOUNDED,
38 MODEL_STATUS_REACHED_DUAL_OBJECTIVE_VALUE_UPPER_BOUND
39 as MODEL_STATUS_RDOVUB,
40 MODEL_STATUS_REACHED_OBJECTIVE_TARGET,
41 MODEL_STATUS_REACHED_TIME_LIMIT,
42 MODEL_STATUS_REACHED_ITERATION_LIMIT,
44 HIGHS_SIMPLEX_STRATEGY_CHOOSE,
45 HIGHS_SIMPLEX_STRATEGY_DUAL,
47 HIGHS_SIMPLEX_CRASH_STRATEGY_OFF,
49 HIGHS_SIMPLEX_EDGE_WEIGHT_STRATEGY_CHOOSE,
50 HIGHS_SIMPLEX_EDGE_WEIGHT_STRATEGY_DANTZIG,
51 HIGHS_SIMPLEX_EDGE_WEIGHT_STRATEGY_DEVEX,
52 HIGHS_SIMPLEX_EDGE_WEIGHT_STRATEGY_STEEPEST_EDGE,
54 HIGHS_VAR_TYPE_CONTINUOUS,
55)
56from scipy.sparse import csc_matrix, vstack, issparse
59def _highs_to_scipy_status_message(highs_status, highs_message):
60 """Converts HiGHS status number/message to SciPy status number/message"""
62 scipy_statuses_messages = {
63 None: (4, "HiGHS did not provide a status code. "),
64 MODEL_STATUS_NOTSET: (4, ""),
65 MODEL_STATUS_LOAD_ERROR: (4, ""),
66 MODEL_STATUS_MODEL_ERROR: (2, ""),
67 MODEL_STATUS_PRESOLVE_ERROR: (4, ""),
68 MODEL_STATUS_SOLVE_ERROR: (4, ""),
69 MODEL_STATUS_POSTSOLVE_ERROR: (4, ""),
70 MODEL_STATUS_MODEL_EMPTY: (4, ""),
71 MODEL_STATUS_RDOVUB: (4, ""),
72 MODEL_STATUS_REACHED_OBJECTIVE_TARGET: (4, ""),
73 MODEL_STATUS_OPTIMAL: (0, "Optimization terminated successfully. "),
74 MODEL_STATUS_REACHED_TIME_LIMIT: (1, "Time limit reached. "),
75 MODEL_STATUS_REACHED_ITERATION_LIMIT: (1, "Iteration limit reached. "),
76 MODEL_STATUS_INFEASIBLE: (2, "The problem is infeasible. "),
77 MODEL_STATUS_UNBOUNDED: (3, "The problem is unbounded. "),
78 MODEL_STATUS_UNBOUNDED_OR_INFEASIBLE: (4, "The problem is unbounded "
79 "or infeasible. ")}
80 unrecognized = (4, "The HiGHS status code was not recognized. ")
81 scipy_status, scipy_message = (
82 scipy_statuses_messages.get(highs_status, unrecognized))
83 scipy_message = (f"{scipy_message}"
84 f"(HiGHS Status {highs_status}: {highs_message})")
85 return scipy_status, scipy_message
88def _replace_inf(x):
89 # Replace `np.inf` with CONST_INF
90 infs = np.isinf(x)
91 x[infs] = np.sign(x[infs])*CONST_INF
92 return x
95def _convert_to_highs_enum(option, option_str, choices):
96 # If option is in the choices we can look it up, if not use
97 # the default value taken from function signature and warn:
98 try:
99 return choices[option.lower()]
100 except AttributeError:
101 return choices[option]
102 except KeyError:
103 sig = inspect.signature(_linprog_highs)
104 default_str = sig.parameters[option_str].default
105 warn(f"Option {option_str} is {option}, but only values in "
106 f"{set(choices.keys())} are allowed. Using default: "
107 f"{default_str}.",
108 OptimizeWarning, stacklevel=3)
109 return choices[default_str]
112def _linprog_highs(lp, solver, time_limit=None, presolve=True,
113 disp=False, maxiter=None,
114 dual_feasibility_tolerance=None,
115 primal_feasibility_tolerance=None,
116 ipm_optimality_tolerance=None,
117 simplex_dual_edge_weight_strategy=None,
118 mip_rel_gap=None,
119 mip_max_nodes=None,
120 **unknown_options):
121 r"""
122 Solve the following linear programming problem using one of the HiGHS
123 solvers:
125 User-facing documentation is in _linprog_doc.py.
127 Parameters
128 ----------
129 lp : _LPProblem
130 A ``scipy.optimize._linprog_util._LPProblem`` ``namedtuple``.
131 solver : "ipm" or "simplex" or None
132 Which HiGHS solver to use. If ``None``, "simplex" will be used.
134 Options
135 -------
136 maxiter : int
137 The maximum number of iterations to perform in either phase. For
138 ``solver='ipm'``, this does not include the number of crossover
139 iterations. Default is the largest possible value for an ``int``
140 on the platform.
141 disp : bool
142 Set to ``True`` if indicators of optimization status are to be printed
143 to the console each iteration; default ``False``.
144 time_limit : float
145 The maximum time in seconds allotted to solve the problem; default is
146 the largest possible value for a ``double`` on the platform.
147 presolve : bool
148 Presolve attempts to identify trivial infeasibilities,
149 identify trivial unboundedness, and simplify the problem before
150 sending it to the main solver. It is generally recommended
151 to keep the default setting ``True``; set to ``False`` if presolve is
152 to be disabled.
153 dual_feasibility_tolerance : double
154 Dual feasibility tolerance. Default is 1e-07.
155 The minimum of this and ``primal_feasibility_tolerance``
156 is used for the feasibility tolerance when ``solver='ipm'``.
157 primal_feasibility_tolerance : double
158 Primal feasibility tolerance. Default is 1e-07.
159 The minimum of this and ``dual_feasibility_tolerance``
160 is used for the feasibility tolerance when ``solver='ipm'``.
161 ipm_optimality_tolerance : double
162 Optimality tolerance for ``solver='ipm'``. Default is 1e-08.
163 Minimum possible value is 1e-12 and must be smaller than the largest
164 possible value for a ``double`` on the platform.
165 simplex_dual_edge_weight_strategy : str (default: None)
166 Strategy for simplex dual edge weights. The default, ``None``,
167 automatically selects one of the following.
169 ``'dantzig'`` uses Dantzig's original strategy of choosing the most
170 negative reduced cost.
172 ``'devex'`` uses the strategy described in [15]_.
174 ``steepest`` uses the exact steepest edge strategy as described in
175 [16]_.
177 ``'steepest-devex'`` begins with the exact steepest edge strategy
178 until the computation is too costly or inexact and then switches to
179 the devex method.
181 Curently, using ``None`` always selects ``'steepest-devex'``, but this
182 may change as new options become available.
184 mip_max_nodes : int
185 The maximum number of nodes allotted to solve the problem; default is
186 the largest possible value for a ``HighsInt`` on the platform.
187 Ignored if not using the MIP solver.
188 unknown_options : dict
189 Optional arguments not used by this particular solver. If
190 ``unknown_options`` is non-empty, a warning is issued listing all
191 unused options.
193 Returns
194 -------
195 sol : dict
196 A dictionary consisting of the fields:
198 x : 1D array
199 The values of the decision variables that minimizes the
200 objective function while satisfying the constraints.
201 fun : float
202 The optimal value of the objective function ``c @ x``.
203 slack : 1D array
204 The (nominally positive) values of the slack,
205 ``b_ub - A_ub @ x``.
206 con : 1D array
207 The (nominally zero) residuals of the equality constraints,
208 ``b_eq - A_eq @ x``.
209 success : bool
210 ``True`` when the algorithm succeeds in finding an optimal
211 solution.
212 status : int
213 An integer representing the exit status of the algorithm.
215 ``0`` : Optimization terminated successfully.
217 ``1`` : Iteration or time limit reached.
219 ``2`` : Problem appears to be infeasible.
221 ``3`` : Problem appears to be unbounded.
223 ``4`` : The HiGHS solver ran into a problem.
225 message : str
226 A string descriptor of the exit status of the algorithm.
227 nit : int
228 The total number of iterations performed.
229 For ``solver='simplex'``, this includes iterations in all
230 phases. For ``solver='ipm'``, this does not include
231 crossover iterations.
232 crossover_nit : int
233 The number of primal/dual pushes performed during the
234 crossover routine for ``solver='ipm'``. This is ``0``
235 for ``solver='simplex'``.
236 ineqlin : OptimizeResult
237 Solution and sensitivity information corresponding to the
238 inequality constraints, `b_ub`. A dictionary consisting of the
239 fields:
241 residual : np.ndnarray
242 The (nominally positive) values of the slack variables,
243 ``b_ub - A_ub @ x``. This quantity is also commonly
244 referred to as "slack".
246 marginals : np.ndarray
247 The sensitivity (partial derivative) of the objective
248 function with respect to the right-hand side of the
249 inequality constraints, `b_ub`.
251 eqlin : OptimizeResult
252 Solution and sensitivity information corresponding to the
253 equality constraints, `b_eq`. A dictionary consisting of the
254 fields:
256 residual : np.ndarray
257 The (nominally zero) residuals of the equality constraints,
258 ``b_eq - A_eq @ x``.
260 marginals : np.ndarray
261 The sensitivity (partial derivative) of the objective
262 function with respect to the right-hand side of the
263 equality constraints, `b_eq`.
265 lower, upper : OptimizeResult
266 Solution and sensitivity information corresponding to the
267 lower and upper bounds on decision variables, `bounds`.
269 residual : np.ndarray
270 The (nominally positive) values of the quantity
271 ``x - lb`` (lower) or ``ub - x`` (upper).
273 marginals : np.ndarray
274 The sensitivity (partial derivative) of the objective
275 function with respect to the lower and upper
276 `bounds`.
278 mip_node_count : int
279 The number of subproblems or "nodes" solved by the MILP
280 solver. Only present when `integrality` is not `None`.
282 mip_dual_bound : float
283 The MILP solver's final estimate of the lower bound on the
284 optimal solution. Only present when `integrality` is not
285 `None`.
287 mip_gap : float
288 The difference between the final objective function value
289 and the final dual bound, scaled by the final objective
290 function value. Only present when `integrality` is not
291 `None`.
293 Notes
294 -----
295 The result fields `ineqlin`, `eqlin`, `lower`, and `upper` all contain
296 `marginals`, or partial derivatives of the objective function with respect
297 to the right-hand side of each constraint. These partial derivatives are
298 also referred to as "Lagrange multipliers", "dual values", and
299 "shadow prices". The sign convention of `marginals` is opposite that
300 of Lagrange multipliers produced by many nonlinear solvers.
302 References
303 ----------
304 .. [15] Harris, Paula MJ. "Pivot selection methods of the Devex LP code."
305 Mathematical programming 5.1 (1973): 1-28.
306 .. [16] Goldfarb, Donald, and John Ker Reid. "A practicable steepest-edge
307 simplex algorithm." Mathematical Programming 12.1 (1977): 361-371.
308 """
310 _check_unknown_options(unknown_options)
312 # Map options to HiGHS enum values
313 simplex_dual_edge_weight_strategy_enum = _convert_to_highs_enum(
314 simplex_dual_edge_weight_strategy,
315 'simplex_dual_edge_weight_strategy',
316 choices={'dantzig': HIGHS_SIMPLEX_EDGE_WEIGHT_STRATEGY_DANTZIG,
317 'devex': HIGHS_SIMPLEX_EDGE_WEIGHT_STRATEGY_DEVEX,
318 'steepest-devex': HIGHS_SIMPLEX_EDGE_WEIGHT_STRATEGY_CHOOSE,
319 'steepest':
320 HIGHS_SIMPLEX_EDGE_WEIGHT_STRATEGY_STEEPEST_EDGE,
321 None: None})
323 c, A_ub, b_ub, A_eq, b_eq, bounds, x0, integrality = lp
325 lb, ub = bounds.T.copy() # separate bounds, copy->C-cntgs
326 # highs_wrapper solves LHS <= A*x <= RHS, not equality constraints
327 lhs_ub = -np.ones_like(b_ub)*np.inf # LHS of UB constraints is -inf
328 rhs_ub = b_ub # RHS of UB constraints is b_ub
329 lhs_eq = b_eq # Equality constaint is inequality
330 rhs_eq = b_eq # constraint with LHS=RHS
331 lhs = np.concatenate((lhs_ub, lhs_eq))
332 rhs = np.concatenate((rhs_ub, rhs_eq))
334 if issparse(A_ub) or issparse(A_eq):
335 A = vstack((A_ub, A_eq))
336 else:
337 A = np.vstack((A_ub, A_eq))
338 A = csc_matrix(A)
340 options = {
341 'presolve': presolve,
342 'sense': HIGHS_OBJECTIVE_SENSE_MINIMIZE,
343 'solver': solver,
344 'time_limit': time_limit,
345 'highs_debug_level': MESSAGE_LEVEL_NONE,
346 'dual_feasibility_tolerance': dual_feasibility_tolerance,
347 'ipm_optimality_tolerance': ipm_optimality_tolerance,
348 'log_to_console': disp,
349 'mip_max_nodes': mip_max_nodes,
350 'output_flag': disp,
351 'primal_feasibility_tolerance': primal_feasibility_tolerance,
352 'simplex_dual_edge_weight_strategy':
353 simplex_dual_edge_weight_strategy_enum,
354 'simplex_strategy': HIGHS_SIMPLEX_STRATEGY_DUAL,
355 'simplex_crash_strategy': HIGHS_SIMPLEX_CRASH_STRATEGY_OFF,
356 'ipm_iteration_limit': maxiter,
357 'simplex_iteration_limit': maxiter,
358 'mip_rel_gap': mip_rel_gap,
359 }
361 # np.inf doesn't work; use very large constant
362 rhs = _replace_inf(rhs)
363 lhs = _replace_inf(lhs)
364 lb = _replace_inf(lb)
365 ub = _replace_inf(ub)
367 if integrality is None or np.sum(integrality) == 0:
368 integrality = np.empty(0)
369 else:
370 integrality = np.array(integrality)
372 res = _highs_wrapper(c, A.indptr, A.indices, A.data, lhs, rhs,
373 lb, ub, integrality.astype(np.uint8), options)
375 # HiGHS represents constraints as lhs/rhs, so
376 # Ax + s = b => Ax = b - s
377 # and we need to split up s by A_ub and A_eq
378 if 'slack' in res:
379 slack = res['slack']
380 con = np.array(slack[len(b_ub):])
381 slack = np.array(slack[:len(b_ub)])
382 else:
383 slack, con = None, None
385 # lagrange multipliers for equalities/inequalities and upper/lower bounds
386 if 'lambda' in res:
387 lamda = res['lambda']
388 marg_ineqlin = np.array(lamda[:len(b_ub)])
389 marg_eqlin = np.array(lamda[len(b_ub):])
390 marg_upper = np.array(res['marg_bnds'][1, :])
391 marg_lower = np.array(res['marg_bnds'][0, :])
392 else:
393 marg_ineqlin, marg_eqlin = None, None
394 marg_upper, marg_lower = None, None
396 # this needs to be updated if we start choosing the solver intelligently
397 solvers = {"ipm": "highs-ipm", "simplex": "highs-ds", None: "highs-ds"}
399 # Convert to scipy-style status and message
400 highs_status = res.get('status', None)
401 highs_message = res.get('message', None)
402 status, message = _highs_to_scipy_status_message(highs_status,
403 highs_message)
405 x = np.array(res['x']) if 'x' in res else None
406 sol = {'x': x,
407 'slack': slack,
408 'con': con,
409 'ineqlin': OptimizeResult({
410 'residual': slack,
411 'marginals': marg_ineqlin,
412 }),
413 'eqlin': OptimizeResult({
414 'residual': con,
415 'marginals': marg_eqlin,
416 }),
417 'lower': OptimizeResult({
418 'residual': None if x is None else x - lb,
419 'marginals': marg_lower,
420 }),
421 'upper': OptimizeResult({
422 'residual': None if x is None else ub - x,
423 'marginals': marg_upper
424 }),
425 'fun': res.get('fun'),
426 'status': status,
427 'success': res['status'] == MODEL_STATUS_OPTIMAL,
428 'message': message,
429 'nit': res.get('simplex_nit', 0) or res.get('ipm_nit', 0),
430 'crossover_nit': res.get('crossover_nit'),
431 }
433 if np.any(x) and integrality is not None:
434 sol.update({
435 'mip_node_count': res.get('mip_node_count', 0),
436 'mip_dual_bound': res.get('mip_dual_bound', 0.0),
437 'mip_gap': res.get('mip_gap', 0.0),
438 })
440 return sol