MAINT: more 1.9.2 backports prep

.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

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 <https://github.com/scipy/scipy/issues/16545>`__: BUG: 1.9.0rc1: \`OptimizeResult\` not populated when \`optimize.milp\`...
 * `#16569 <https://github.com/scipy/scipy/issues/16569>`__: BUG: \`sparse.hstack\` returns incorrect result when the stack...
 * `#16898 <https://github.com/scipy/scipy/issues/16898>`__: BUG: optimize.minimize backwards compatability in scipy 1.9
 * `#16935 <https://github.com/scipy/scipy/issues/16935>`__: BUG: using msvc + meson to build scipy --> cl cannot be used...
 * `#16952 <https://github.com/scipy/scipy/issues/16952>`__: BUG: error from \`scipy.stats.mode\` with \`NaN\`s, \`axis !=...
 * `#16964 <https://github.com/scipy/scipy/issues/16964>`__: BUG: scipy 1.7.3 wheels on PyPI require numpy<1.23 in contradiction...
 * `#17026 <https://github.com/scipy/scipy/issues/17026>`__: BUG: ncf_gen::ppf(..) causes segfault
+* `#17050 <https://github.com/scipy/scipy/issues/17050>`__: Pearson3 PPF does not function properly with negative skew.
 * `#17124 <https://github.com/scipy/scipy/issues/17124>`__: BUG: OSX-64 Test failure test_ppf_against_tables getting NaN
 
+
 Pull requests for 1.9.2
 -----------------------
 
 * `#16628 <https://github.com/scipy/scipy/pull/16628>`__: FIX: Updated dtype resolution in \`_stack_along_minor_axis\`
+* `#16814 <https://github.com/scipy/scipy/pull/16814>`__: FIX: milp: return feasible solutions if available on time out
 * `#16842 <https://github.com/scipy/scipy/pull/16842>`__: ENH: cibuildwheel infrastructure
 * `#16909 <https://github.com/scipy/scipy/pull/16909>`__: MAINT: minimize, restore squeezed ((1.0)) addresses #16898
 * `#16911 <https://github.com/scipy/scipy/pull/16911>`__: REL: prep for SciPy 1.9.2
@@ -58,7 +63,9 @@ Pull requests for 1.9.2
 * `#17011 <https://github.com/scipy/scipy/pull/17011>`__: Rudimentary test for manylinux_aarch64 with cibuildwheel
 * `#17013 <https://github.com/scipy/scipy/pull/17013>`__: BLD: make MKL detection a little more robust, add notes on TODOs
 * `#17046 <https://github.com/scipy/scipy/pull/17046>`__: CI: Update cibuildwheel to 2.10.1
+* `#17055 <https://github.com/scipy/scipy/pull/17055>`__: MAINT: stats.pearson3: fix ppf for negative skew
 * `#17064 <https://github.com/scipy/scipy/pull/17064>`__: BUG: Fix numerical precision error of \`truncnorm.logcdf\` when...
 * `#17096 <https://github.com/scipy/scipy/pull/17096>`__: FIX: ensure a hold on GIL before raising warnings/errors
 * `#17127 <https://github.com/scipy/scipy/pull/17127>`__: TST: stats.studentized_range: fix incorrect test
 * `#17131 <https://github.com/scipy/scipy/pull/17131>`__: MAINT: pyproject.toml: Update build system requirements
+* `#17132 <https://github.com/scipy/scipy/pull/17132>`__: MAINT: 1.9.2 backports

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

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,21 +671,36 @@ 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
     cdef HighsSolution solution
     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': <int> model_status,
             'message': f'model_status is {highs.modelStatusToString(model_status).decode()}; '
-                       f'primal_status is {utilBasisStatusToString(<HighsBasisStatus> info.primal_solution_status)}',
+                       f'primal_status is {utilBasisStatusToString(<HighsBasisStatus> 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': <int> model_status,
             'message': highs.modelStatusToString(model_status).decode(),
-            'unscaled_status': <int> unscaled_model_status,
 
             # Primal solution
             'x': [solution.col_value[ii] for ii in range(numcol)],

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':

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

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.

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

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,)],

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):

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])