From 6fb67007dd7105755057f3379fb7ef423eae524e Mon Sep 17 00:00:00 2001 From: Nicholas McKibben Date: Sun, 2 Oct 2022 17:45:48 -0700 Subject: [PATCH 1/4] FIX: optimize.milp: return feasible solution if available on timeout/node limit (#16814) * FIX: optimize.milp: return feasible solution if available on timeout/node limit Co-authored-by: Matt Haberland Co-authored-by: laenNoCode <36074576+laenNoCode@users.noreply.github.com> --- scipy/optimize/_highs/cython/src/Highs.pxd | 1 - .../_highs/cython/src/_highs_wrapper.pyx | 28 ++++++++++--- scipy/optimize/_linprog_highs.py | 6 +++ scipy/optimize/_milp.py | 8 +++- scipy/optimize/tests/test_milp.py | 42 +++++++++++++++++++ 5 files changed, 76 insertions(+), 9 deletions(-) diff --git a/scipy/optimize/_highs/cython/src/Highs.pxd b/scipy/optimize/_highs/cython/src/Highs.pxd index 7312b799004e..b5615907f8e4 100644 --- a/scipy/optimize/_highs/cython/src/Highs.pxd +++ b/scipy/optimize/_highs/cython/src/Highs.pxd @@ -30,7 +30,6 @@ cdef extern from "Highs.h": # split up for cython below #const HighsModelStatus& getModelStatus(const bool scaled_model = False) const const HighsModelStatus & getModelStatus() const - const HighsModelStatus & getModelStatus(const bool scaled_model) const const HighsInfo& getHighsInfo "getInfo" () const string modelStatusToString(const HighsModelStatus model_status) const diff --git a/scipy/optimize/_highs/cython/src/_highs_wrapper.pyx b/scipy/optimize/_highs/cython/src/_highs_wrapper.pyx index 422a8e3aa638..2069763119ea 100644 --- a/scipy/optimize/_highs/cython/src/_highs_wrapper.pyx +++ b/scipy/optimize/_highs/cython/src/_highs_wrapper.pyx @@ -17,6 +17,7 @@ from .HighsIO cimport ( kWarning, ) from .HConst cimport ( + HIGHS_CONST_INF, HighsModelStatus, HighsModelStatusNOTSET, HighsModelStatusLOAD_ERROR, @@ -139,6 +140,7 @@ cdef apply_options(dict options, Highs & highs): 'ipm_iteration_limit', 'keep_n_rows', 'max_threads', + 'mip_max_nodes', 'highs_debug_level', 'min_threads', 'simplex_crash_strategy', @@ -669,8 +671,6 @@ def _highs_wrapper( # Extract what we need from the solution cdef HighsModelStatus model_status = highs.getModelStatus() - cdef HighsModelStatus scaled_model_status = highs.getModelStatus(True) - cdef HighsModelStatus unscaled_model_status = model_status # We might need an info object if we can look up the solution and a place to put solution cdef HighsInfo info = highs.getHighsInfo() # it should always be safe to get the info object @@ -678,12 +678,29 @@ def _highs_wrapper( cdef HighsBasis basis cdef double[:, ::1] marg_bnds = np.zeros((2, numcol)) # marg_bnds[0, :]: lower - # If the status is bad, don't look up the solution - if model_status != HighsModelStatusOPTIMAL: + # Failure modes: + # LP: if we have anything other than an Optimal status, it + # is unsafe (and unhelpful) to read any results + # MIP: has a non-Optimal status or has timed out/reached max iterations + # 1) If not Optimal/TimedOut/MaxIter status, there is no solution + # 2) If TimedOut/MaxIter status, there may be a feasible solution. + # if the objective function value is not Infinity, then the + # current solution is feasible and can be returned. Else, there + # is no solution. + mipFailCondition = model_status not in { + HighsModelStatusOPTIMAL, + HighsModelStatusREACHED_TIME_LIMIT, + HighsModelStatusREACHED_ITERATION_LIMIT, + } or (model_status in { + HighsModelStatusREACHED_TIME_LIMIT, + HighsModelStatusREACHED_ITERATION_LIMIT, + } and (info.objective_function_value == HIGHS_CONST_INF)) + lpFailCondition = model_status != HighsModelStatusOPTIMAL + if (highs.getLp().isMip() and mipFailCondition) or (not highs.getLp().isMip() and lpFailCondition): return { 'status': model_status, 'message': f'model_status is {highs.modelStatusToString(model_status).decode()}; ' - f'primal_status is {utilBasisStatusToString( info.primal_solution_status)}', + f'primal_status is {utilBasisStatusToString( info.primal_solution_status).decode()}', 'simplex_nit': info.simplex_iteration_count, 'ipm_nit': info.ipm_iteration_count, 'fun': None, @@ -705,7 +722,6 @@ def _highs_wrapper( res = { 'status': model_status, 'message': highs.modelStatusToString(model_status).decode(), - 'unscaled_status': unscaled_model_status, # Primal solution 'x': [solution.col_value[ii] for ii in range(numcol)], diff --git a/scipy/optimize/_linprog_highs.py b/scipy/optimize/_linprog_highs.py index c4015e559321..1f8b1878ace8 100644 --- a/scipy/optimize/_linprog_highs.py +++ b/scipy/optimize/_linprog_highs.py @@ -115,6 +115,7 @@ def _linprog_highs(lp, solver, time_limit=None, presolve=True, primal_feasibility_tolerance=None, ipm_optimality_tolerance=None, simplex_dual_edge_weight_strategy=None, + mip_max_nodes=None, **unknown_options): r""" Solve the following linear programming problem using one of the HiGHS @@ -179,6 +180,10 @@ def _linprog_highs(lp, solver, time_limit=None, presolve=True, Curently, using ``None`` always selects ``'steepest-devex'``, but this may change as new options become available. + mip_max_nodes : int + The maximum number of nodes allotted to solve the problem; default is + the largest possible value for a ``HighsInt`` on the platform. + Ignored if not using the MIP solver. unknown_options : dict Optional arguments not used by this particular solver. If ``unknown_options`` is non-empty, a warning is issued listing all @@ -339,6 +344,7 @@ def _linprog_highs(lp, solver, time_limit=None, presolve=True, 'dual_feasibility_tolerance': dual_feasibility_tolerance, 'ipm_optimality_tolerance': ipm_optimality_tolerance, 'log_to_console': disp, + 'mip_max_nodes': mip_max_nodes, 'output_flag': disp, 'primal_feasibility_tolerance': primal_feasibility_tolerance, 'simplex_dual_edge_weight_strategy': diff --git a/scipy/optimize/_milp.py b/scipy/optimize/_milp.py index 99836bd9697f..adf77fe89cba 100644 --- a/scipy/optimize/_milp.py +++ b/scipy/optimize/_milp.py @@ -130,13 +130,14 @@ def _milp_iv(c, integrality, bounds, constraints, options): # options IV options = options or {} - supported_options = {'disp', 'presolve', 'time_limit'} + supported_options = {'disp', 'presolve', 'time_limit', 'node_limit'} unsupported_options = set(options).difference(supported_options) if unsupported_options: message = (f"Unrecognized options detected: {unsupported_options}. " "These will be passed to HiGHS verbatim.") warnings.warn(message, RuntimeWarning, stacklevel=3) - options_iv = {'log_to_console': options.get("disp", False)} + options_iv = {'log_to_console': options.pop("disp", False), + 'mip_max_nodes': options.pop("node_limit", None)} options_iv.update(options) return c, integrality, lb, ub, indptr, indices, data, b_l, b_u, options_iv @@ -225,6 +226,9 @@ def milp(c, *, integrality=None, bounds=None, constraints=None, options=None): disp : bool (default: ``False``) Set to ``True`` if indicators of optimization status are to be printed to the console during optimization. + node_limit : int, optional + The maximum number of nodes (linear program relaxations) to solve + before stopping. Default is no maximum number of nodes. presolve : bool (default: ``True``) Presolve attempts to identify trivial infeasibilities, identify trivial unboundedness, and simplify the problem before diff --git a/scipy/optimize/tests/test_milp.py b/scipy/optimize/tests/test_milp.py index f179029e0b60..694401b220a7 100644 --- a/scipy/optimize/tests/test_milp.py +++ b/scipy/optimize/tests/test_milp.py @@ -272,6 +272,48 @@ def test_infeasible_prob_16609(): np.testing.assert_equal(res.status, 2) +_msg_time = "Time limit reached. (HiGHS Status 13:" +_msg_iter = "Iteration limit reached. (HiGHS Status 14:" + + +@pytest.mark.skipif(np.intp(0).itemsize < 8, + reason="Unhandled 32-bit GCC FP bug") +@pytest.mark.slow +@pytest.mark.timeout(360) +@pytest.mark.parametrize(["options", "msg"], [({"time_limit": 10}, _msg_time), + ({"node_limit": 1}, _msg_iter)]) +def test_milp_timeout_16545(options, msg): + # Ensure solution is not thrown away if MILP solver times out + # -- see gh-16545 + rng = np.random.default_rng(5123833489170494244) + A = rng.integers(0, 5, size=(100, 100)) + b_lb = np.full(100, fill_value=-np.inf) + b_ub = np.full(100, fill_value=25) + constraints = LinearConstraint(A, b_lb, b_ub) + variable_lb = np.zeros(100) + variable_ub = np.ones(100) + variable_bounds = Bounds(variable_lb, variable_ub) + integrality = np.ones(100) + c_vector = -np.ones(100) + res = milp( + c_vector, + integrality=integrality, + bounds=variable_bounds, + constraints=constraints, + options=options, + ) + + assert res.message.startswith(msg) + assert res["x"] is not None + + # ensure solution is feasible + x = res["x"] + tol = 1e-8 # sometimes needed due to finite numerical precision + assert np.all(b_lb - tol <= A @ x) and np.all(A @ x <= b_ub + tol) + assert np.all(variable_lb - tol <= x) and np.all(x <= variable_ub + tol) + assert np.allclose(x, np.round(x)) + + def test_three_constraints_16878(): # `milp` failed when exactly three constraints were passed # Ensure that this is no longer the case. From ed9760e60a28b8f13e5644494033e2dab9aafbcd Mon Sep 17 00:00:00 2001 From: Matt Haberland Date: Mon, 3 Oct 2022 12:39:02 -0700 Subject: [PATCH 2/4] MAINT: stats.pearson3: fix ppf for negative skew (#17055) --- scipy/stats/_continuous_distns.py | 4 +++- scipy/stats/_distr_params.py | 1 + scipy/stats/tests/test_distributions.py | 19 +++++++++++++++++++ scipy/stats/tests/test_fit.py | 4 ++-- 4 files changed, 25 insertions(+), 3 deletions(-) diff --git a/scipy/stats/_continuous_distns.py b/scipy/stats/_continuous_distns.py index 48645ba11dbc..5bc098aded66 100644 --- a/scipy/stats/_continuous_distns.py +++ b/scipy/stats/_continuous_distns.py @@ -6910,7 +6910,9 @@ def _ppf(self, q, skew): ans, q, _, mask, invmask, beta, alpha, zeta = ( self._preprocess(q, skew)) ans[mask] = _norm_ppf(q[mask]) - ans[invmask] = sc.gammaincinv(alpha, q[invmask])/beta + zeta + q = q[invmask] + q[beta < 0] = 1 - q[beta < 0] # for negative skew; see gh-17050 + ans[invmask] = sc.gammaincinv(alpha, q)/beta + zeta return ans @_call_super_mom diff --git a/scipy/stats/_distr_params.py b/scipy/stats/_distr_params.py index 280b83f24ad3..c2d7f05e8360 100644 --- a/scipy/stats/_distr_params.py +++ b/scipy/stats/_distr_params.py @@ -87,6 +87,7 @@ ['norminvgauss', (1.25, 0.5)], ['pareto', (2.621716532144454,)], ['pearson3', (0.1,)], + ['pearson3', (-2,)], ['powerlaw', (1.6591133289905851,)], ['powerlognorm', (2.1413923530064087, 0.44639540782048337)], ['powernorm', (4.4453652254590779,)], diff --git a/scipy/stats/tests/test_distributions.py b/scipy/stats/tests/test_distributions.py index 4d49034009bb..e8f899832ab0 100644 --- a/scipy/stats/tests/test_distributions.py +++ b/scipy/stats/tests/test_distributions.py @@ -1788,6 +1788,25 @@ def test_return_array_bug_11746(self): assert_equal(moment, 0) assert_equal(type(moment), float) + def test_ppf_bug_17050(self): + # incorrect PPF for negative skews were reported in gh-17050 + # Check that this is fixed (even in the array case) + skews = [-3, -1, 0, 0.5] + x_eval = 0.5 + res = stats.pearson3.ppf(stats.pearson3.cdf(x_eval, skews), skews) + assert_allclose(res, x_eval) + + # Negation of the skew flips the distribution about the origin, so + # the following should hold + skew = np.array([[-0.5], [1.5]]) + x = np.linspace(-2, 2) + assert_allclose(stats.pearson3.pdf(x, skew), + stats.pearson3.pdf(-x, -skew)) + assert_allclose(stats.pearson3.cdf(x, skew), + stats.pearson3.sf(-x, -skew)) + assert_allclose(stats.pearson3.ppf(x, skew), + -stats.pearson3.isf(x, -skew)) + class TestKappa4: def test_cdf_genpareto(self): diff --git a/scipy/stats/tests/test_fit.py b/scipy/stats/tests/test_fit.py index 143cc4d23daa..c210b477a293 100644 --- a/scipy/stats/tests/test_fit.py +++ b/scipy/stats/tests/test_fit.py @@ -377,8 +377,8 @@ def test_basic_fit(self, dist_name): dist = getattr(stats, dist_name) shapes = np.array(dist_data[dist_name]) bounds = np.empty((len(shapes) + 2, 2), dtype=np.float64) - bounds[:-2, 0] = shapes/10**np.sign(shapes) - bounds[:-2, 1] = shapes*10**np.sign(shapes) + bounds[:-2, 0] = shapes/10.**np.sign(shapes) + bounds[:-2, 1] = shapes*10.**np.sign(shapes) bounds[-2] = (0, 10) bounds[-1] = (0, 10) loc = rng.uniform(*bounds[-2]) From a6ba7cad3b54c35d2ccb55c595691689004742c1 Mon Sep 17 00:00:00 2001 From: Tyler Reddy Date: Tue, 4 Oct 2022 21:01:23 -0600 Subject: [PATCH 3/4] MAINT: misc 1.9.2 updates * update 1.9.2 relnotes * uncomment wheel upload code --- .github/workflows/wheels.yml | 35 ++++++++++++++++------------------- doc/release/1.9.2-notes.rst | 17 ++++++++++++----- 2 files changed, 28 insertions(+), 24 deletions(-) diff --git a/.github/workflows/wheels.yml b/.github/workflows/wheels.yml index 7cce67b87990..7a9d9ba5283b 100644 --- a/.github/workflows/wheels.yml +++ b/.github/workflows/wheels.yml @@ -199,22 +199,19 @@ jobs: path: ./wheelhouse/*.whl name: ${{ matrix.python[0] }}-${{ matrix.buildplat[1] }} - # TODO uncomment when those responsible for uploading - # nightly/release wheels want to make this script live. - - # - name: Upload wheels - # if: success() - # shell: bash - # env: - # SCIPY_STAGING_UPLOAD_TOKEN: ${{ secrets.SCIPY_STAGING_UPLOAD_TOKEN }} - # SCIPY_NIGHTLY_UPLOAD_TOKEN: ${{ secrets.SCIPY_NIGHTLY_UPLOAD_TOKEN }} - # run: | - # source tools/wheels/upload_wheels.sh - # set_upload_vars - # # trigger an upload to - # # https://anaconda.org/scipy-wheels-nightly/scipy - # # for cron jobs or "Run workflow" (restricted to main branch). - # # Tags will upload to - # # https://anaconda.org/multibuild-wheels-staging/scipy - # # The tokens were originally generated at anaconda.org - # upload_wheels + - name: Upload wheels + if: success() + shell: bash + env: + SCIPY_STAGING_UPLOAD_TOKEN: ${{ secrets.SCIPY_STAGING_UPLOAD_TOKEN }} + SCIPY_NIGHTLY_UPLOAD_TOKEN: ${{ secrets.SCIPY_NIGHTLY_UPLOAD_TOKEN }} + run: | + source tools/wheels/upload_wheels.sh + set_upload_vars + # trigger an upload to + # https://anaconda.org/scipy-wheels-nightly/scipy + # for cron jobs or "Run workflow" (restricted to main branch). + # Tags will upload to + # https://anaconda.org/multibuild-wheels-staging/scipy + # The tokens were originally generated at anaconda.org + upload_wheels diff --git a/doc/release/1.9.2-notes.rst b/doc/release/1.9.2-notes.rst index 2cdbde8a54dc..f5a0cb4f6844 100644 --- a/doc/release/1.9.2-notes.rst +++ b/doc/release/1.9.2-notes.rst @@ -13,37 +13,42 @@ Authors * Hood Chatham (1) * Thomas J. Fan (1) -* Ralf Gommers (7) -* Matt Haberland (3) +* Ralf Gommers (22) +* Matt Haberland (5) * Julien Jerphanion (1) * Loïc Estève (1) -* Nicholas McKibben (1) +* Nicholas McKibben (2) * Naoto Mizuno (1) * Andrew Nelson (3) -* Tyler Reddy (21) +* Tyler Reddy (28) * Pamphile Roy (1) * Ewout ter Hoeven (2) +* Warren Weckesser (1) * Meekail Zain (1) + -A total of 13 people contributed to this release. +A total of 14 people contributed to this release. People with a "+" by their names contributed a patch for the first time. This list of names is automatically generated, and may not be fully complete. Issues closed for 1.9.2 ----------------------- +* `#16545 `__: BUG: 1.9.0rc1: \`OptimizeResult\` not populated when \`optimize.milp\`... * `#16569 `__: BUG: \`sparse.hstack\` returns incorrect result when the stack... * `#16898 `__: BUG: optimize.minimize backwards compatability in scipy 1.9 * `#16935 `__: BUG: using msvc + meson to build scipy --> cl cannot be used... * `#16952 `__: BUG: error from \`scipy.stats.mode\` with \`NaN\`s, \`axis !=... * `#16964 `__: BUG: scipy 1.7.3 wheels on PyPI require numpy<1.23 in contradiction... * `#17026 `__: BUG: ncf_gen::ppf(..) causes segfault +* `#17050 `__: Pearson3 PPF does not function properly with negative skew. * `#17124 `__: BUG: OSX-64 Test failure test_ppf_against_tables getting NaN + Pull requests for 1.9.2 ----------------------- * `#16628 `__: FIX: Updated dtype resolution in \`_stack_along_minor_axis\` +* `#16814 `__: FIX: milp: return feasible solutions if available on time out * `#16842 `__: ENH: cibuildwheel infrastructure * `#16909 `__: MAINT: minimize, restore squeezed ((1.0)) addresses #16898 * `#16911 `__: REL: prep for SciPy 1.9.2 @@ -58,7 +63,9 @@ Pull requests for 1.9.2 * `#17011 `__: Rudimentary test for manylinux_aarch64 with cibuildwheel * `#17013 `__: BLD: make MKL detection a little more robust, add notes on TODOs * `#17046 `__: CI: Update cibuildwheel to 2.10.1 +* `#17055 `__: MAINT: stats.pearson3: fix ppf for negative skew * `#17064 `__: BUG: Fix numerical precision error of \`truncnorm.logcdf\` when... * `#17096 `__: FIX: ensure a hold on GIL before raising warnings/errors * `#17127 `__: TST: stats.studentized_range: fix incorrect test * `#17131 `__: MAINT: pyproject.toml: Update build system requirements +* `#17132 `__: MAINT: 1.9.2 backports From 6b098c25223e224ff44101f86bbc86efecffe1d9 Mon Sep 17 00:00:00 2001 From: Matt Haberland Date: Thu, 6 Oct 2022 11:00:27 -0700 Subject: [PATCH 4/4] TST: optimize.milp: remove problematic timeout/iteration test --- scipy/optimize/tests/test_milp.py | 42 ------------------------------- 1 file changed, 42 deletions(-) diff --git a/scipy/optimize/tests/test_milp.py b/scipy/optimize/tests/test_milp.py index 694401b220a7..f179029e0b60 100644 --- a/scipy/optimize/tests/test_milp.py +++ b/scipy/optimize/tests/test_milp.py @@ -272,48 +272,6 @@ def test_infeasible_prob_16609(): np.testing.assert_equal(res.status, 2) -_msg_time = "Time limit reached. (HiGHS Status 13:" -_msg_iter = "Iteration limit reached. (HiGHS Status 14:" - - -@pytest.mark.skipif(np.intp(0).itemsize < 8, - reason="Unhandled 32-bit GCC FP bug") -@pytest.mark.slow -@pytest.mark.timeout(360) -@pytest.mark.parametrize(["options", "msg"], [({"time_limit": 10}, _msg_time), - ({"node_limit": 1}, _msg_iter)]) -def test_milp_timeout_16545(options, msg): - # Ensure solution is not thrown away if MILP solver times out - # -- see gh-16545 - rng = np.random.default_rng(5123833489170494244) - A = rng.integers(0, 5, size=(100, 100)) - b_lb = np.full(100, fill_value=-np.inf) - b_ub = np.full(100, fill_value=25) - constraints = LinearConstraint(A, b_lb, b_ub) - variable_lb = np.zeros(100) - variable_ub = np.ones(100) - variable_bounds = Bounds(variable_lb, variable_ub) - integrality = np.ones(100) - c_vector = -np.ones(100) - res = milp( - c_vector, - integrality=integrality, - bounds=variable_bounds, - constraints=constraints, - options=options, - ) - - assert res.message.startswith(msg) - assert res["x"] is not None - - # ensure solution is feasible - x = res["x"] - tol = 1e-8 # sometimes needed due to finite numerical precision - assert np.all(b_lb - tol <= A @ x) and np.all(A @ x <= b_ub + tol) - assert np.all(variable_lb - tol <= x) and np.all(x <= variable_ub + tol) - assert np.allclose(x, np.round(x)) - - def test_three_constraints_16878(): # `milp` failed when exactly three constraints were passed # Ensure that this is no longer the case.