Repository: numpy/numpy
Pull request: BUG: Fix np.quantile([0, 1], 0, method='weibull')
Issue https://github.com/numpy/numpy/issues/24592 was opened because the `np.quantile` method gives incorrect answers for the `weibull` method. The problem:

```   
import numpy as np
x=np.arange(10)
np.quantile(x, 0., method='weibull') # output 0.0, correct
np.quantile(x, 0, method='weibull') # output 9, incorrect!
np.quantile(x, 1, method='weibull') # raises an error
```
The value of `np.quantile(x, q=0, method='weibull')` is incorrect for integer inputs. Two options to address this:
    
1. Do not allow integer input for `np.quantile`. The documentation for the probability argument states that "q : array_like of float", so this seems reasonable. But this would be a behaviour change that might have a large impact and there are unit tests that test for integer input as well.

2. Automatically cast the input argument `q` to float. Also this would be a behaviour change and break existing tests.

We rejected these two options because of the issues mentioned. Applying either one of the options only for `weibull` would make the input handling less uniform.

The problem with `weibull` (and some other methods) is that inside `np.lib._function_base_impl._quantile` there are two paths: one for integer virtual indices and one for other types of virtual indices (in the code the statement `if np.issubdtype(virtual_indexes.dtype, np.integer):`. The integer path is only valid for methods returning integer indices in the range `[0, size(arr)>`. For the `weibull` method something like `q = [0, .5, 1]` is cast to float, which results in float `virtual_indices` (which contains floats with integer values, but that is okay).

We can either:
    
i) Modify the code in the integer path to handle integers < 0 and > size(arr) -1
ii) Check whether the methods supports integer output, and if not use the alternative path

We opted for the second option.

Fixes #24592                                                          
                                                          

Problem statement:
BUG: quantile 1 with 'weibull' gives confusing error
### Describe the issue:

Asking for the 1 quantile with `weibull` gives a confusing error - see below.

### Reproduce the code example:

```python
import numpy as np

np.quantile(np.arange(10), 1, method='weibull')
```


### Error message:

```shell
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In [57], line 1
----> 1 np.quantile(np.arange(10), 1, method='weibull')

File <__array_function__ internals>:180, in quantile(*args, **kwargs)

File ~/Library/Python/3.10/lib/python/site-packages/numpy/lib/function_base.py:4412, in quantile(a, q, axis, out, overwrite_input, method, keepdims, interpolation)
   4410 if not _quantile_is_valid(q):
   4411     raise ValueError("Quantiles must be in the range [0, 1]")
-> 4412 return _quantile_unchecked(
   4413     a, q, axis, out, overwrite_input, method, keepdims)

File ~/Library/Python/3.10/lib/python/site-packages/numpy/lib/function_base.py:4424, in _quantile_unchecked(a, q, axis, out, overwrite_input, method, keepdims)
   4416 def _quantile_unchecked(a,
   4417                         q,
   4418                         axis=None,
   (...)
   4421                         method="linear",
   4422                         keepdims=False):
   4423     """Assumes that q is in [0, 1], and is an ndarray"""
-> 4424     r, k = _ureduce(a,
   4425                     func=_quantile_ureduce_func,
   4426                     q=q,
   4427                     axis=axis,
   4428                     out=out,
   4429                     overwrite_input=overwrite_input,
   4430                     method=method)
   4431     if keepdims:
   4432         return r.reshape(q.shape + k)

File ~/Library/Python/3.10/lib/python/site-packages/numpy/lib/function_base.py:3725, in _ureduce(a, func, **kwargs)
   3722 else:
   3723     keepdim = (1,) * a.ndim
-> 3725 r = func(a, **kwargs)
   3726 return r, keepdim

File ~/Library/Python/3.10/lib/python/site-packages/numpy/lib/function_base.py:4593, in _quantile_ureduce_func(a, q, axis, out, overwrite_input, method)
   4591     else:
   4592         arr = a.copy()
-> 4593 result = _quantile(arr,
   4594                    quantiles=q,
   4595                    axis=axis,
   4596                    method=method,
   4597                    out=out)
   4598 return result

File ~/Library/Python/3.10/lib/python/site-packages/numpy/lib/function_base.py:4683, in _quantile(arr, quantiles, axis, method, out)
   4680     slices_having_nans = np.isnan(arr[-1])
   4681 else:
   4682     # cannot contain nan
-> 4683     arr.partition(virtual_indexes.ravel(), axis=0)
   4684     slices_having_nans = np.array(False, dtype=bool)
   4685 result = take(arr, virtual_indexes, axis=0, out=out)

ValueError: kth(=10) out of bounds (10)
```
```


### Runtime information:

```
[1] import sys, numpy; print(numpy.__version__); print(sys.version)
1.26.0b1
3.10.11 (main, Apr  7 2023, 07:24:53) [Clang 14.0.0 (clang-1400.0.29.202)]
In [2]: print(numpy.show_runtime())
WARNING: `threadpoolctl` not found in system! Install it by `pip install threadpoolctl`. Once installed, try `np.show_runtime` again for more detailed build information
[{'numpy_version': '1.26.0b1',
  'python': '3.10.11 (main, Apr  7 2023, 07:24:53) [Clang 14.0.0 '
            '(clang-1400.0.29.202)]',
  'uname': uname_result(system='Darwin', node='nipraxis.dynevor.org', release='22.6.0', version='Darwin Kernel Version 22.6.0: Wed Jul  5 22:17:35 PDT 2023; root:xnu-8796.141.3~6/RELEASE_ARM64_T8112', machine='arm64')},
 {'simd_extensions': {'baseline': ['NEON', 'NEON_FP16', 'NEON_VFPV4', 'ASIMD'],
                      'found': ['ASIMDHP'],
                      'not_found': ['ASIMDFHM']}}]
```

### Context for the issue:

Writing a tutorial on quantiles, and failing to understand them!

Hints from issue discussion:
The problem is you are passing 1 as an integer argument, while the documentation states `q : array_like of float`
The output of `quantile(x, 1., method='weibull')` is correct: `9`.

For input zero however, the results are:
```
In[83]: quantile(x, 0., method='weibull')
Out[83]: 0.0

In[84]: quantile(x, 0, method='weibull')
Out[84]: 9
```
The fact that we get an incorrect answer when passing `0`, I consider to be a bug in the implementation. I will look if this can be easily addressed.


There are a couple more cases affected:
```
import numpy as np
x=np.arange(10)

for m in ['inverted_cdf', 'averaged_inverted_cdf', 'closest_observation', 
          'interpolated_inverted_cdf', 'hazen', 'weibull', 'linear', 'median_unbiased', 'normal_unbiased']:
    for q in [0,1]:
        try:
            print(f'{m}: q {q}: {np.quantile(x, q, method=m)}')
        except:
            print(f'{m}: could not calculate q={q}')     
```
Output:
```
inverted_cdf: q 0: 0
inverted_cdf: q 1: 9
averaged_inverted_cdf: q 0: 9
averaged_inverted_cdf: q 1: 9
closest_observation: q 0: 0
closest_observation: q 1: 9
interpolated_inverted_cdf: q 0: 9
interpolated_inverted_cdf: q 1: 9
hazen: q 0: 0.0
hazen: q 1: 9.0
weibull: q 0: 9
weibull: could not calculate q=1
linear: q 0: 0
linear: q 1: 9
median_unbiased: q 0: 0.0
median_unbiased: q 1: 9.0
normal_unbiased: q 0: 0.0
normal_unbiased: q 1: 9.0
```

Two simple solutions:

* If the argument `q` is of integer type, cast to float
* If the argument `q` is of integer type, raise a deprecation warning or exception.

The problems arise in this code path: https://github.com/numpy/numpy/blob/main/numpy/lib/_function_base_impl.py#L4731. I suspect the special code path for integers is there for performance reasons.

Is there any reason not to convert to float?  I mean - what is the advantage of the error?
Such a change would be backwards incompatible. For example:
```
from fractions import Fraction
import numpy as np

x=[Fraction(i) for i in range(10)]
print(np.quantile(x, 0)) # works 
print(np.quantile(x, 0.)) # fails
```
Not sure whether the `Fraction` is a use case that needs to be supported with float argument `q`. But there are tests for `Fraction`, that suggest we should support it in some way: https://github.com/numpy/numpy/blob/da8cdcf2f127eb609448944ae33c6cb8eb660138/numpy/lib/tests/test_function_base.py#L3582

Adding @eric-wieser as the author of the tests in the CC.

Sorry @eendebakpt - just checking - you mean automatic conversion to float would be backwards incompatible?

Surely the current behavior can't be intentional, of failing with a floating point 0 (but not integer 0) for a fraction value?  What would the purpose of such a failure be?
Yes, the automatic conversion would break one of the tests I linked to, so it would be backwards incompatible. 

I do agree that the current behavior is strange, and perhaps something we can change. But that is up to a numpy core developer.
Just checking - but the point here is that e.g. `np.quantile(Fraction(5), 0)` gives `Fraction(5, 1)`, but gives an error for `np.quantile(Fraction(5), 0.)`.  I'm not sure at all what anyone would use `np.quantile(Fraction(5), 0)` for.   I notice that you get the same answer for `np.quantile(Fraction(5), 1)`.   Why would `np.quantile(Fraction(5), 0.)` be expected to raise an error?  Isn't it surprising that you can only use integers for the second argument in that case? 
> ```python
> import numpy as np
> 
> np.quantile(np.arange(10), 1, method='weibull')
> ```

I am interested tosolve this can you please assign me this issue
> > ```
> > import numpy as np
> > 
> > np.quantile(np.arange(10), 1, method='weibull')
> > ```
> 
> I am interested tosolve this can you please assign me this issue

@Anianonymous There already is a PRs to fix this issue: https://github.com/numpy/numpy/pull/24710. You can help to review the PR.

Commit message: BUG: Fix np.quantile([0, 1], 0, method='weibull')

Patch to review:
diff --git a/numpy/lib/_function_base_impl.py b/numpy/lib/_function_base_impl.py
index 7013ada34bef..6ad17fa904de 100644
--- a/numpy/lib/_function_base_impl.py
+++ b/numpy/lib/_function_base_impl.py
@@ -50,7 +50,7 @@
 # _QuantileMethods is a dictionary listing all the supported methods to
 # compute quantile/percentile.
 #
-# Below virtual_index refer to the index of the element where the percentile
+# Below virtual_index refers to the index of the element where the percentile
 # would be found in the sorted sample.
 # When the sample contains exactly the percentile wanted, the virtual_index is
 # an integer to the index of this element.
@@ -68,7 +68,7 @@
     # Discrete methods
     inverted_cdf=dict(
         get_virtual_index=lambda n, quantiles: _inverted_cdf(n, quantiles),
-        fix_gamma=lambda gamma, _: gamma,  # should never be called
+        fix_gamma=None,  # should never be called
     ),
     averaged_inverted_cdf=dict(
         get_virtual_index=lambda n, quantiles: (n * quantiles) - 1,
@@ -81,7 +81,7 @@
     closest_observation=dict(
         get_virtual_index=lambda n, quantiles: _closest_observation(n,
                                                                     quantiles),
-        fix_gamma=lambda gamma, _: gamma,  # should never be called
+        fix_gamma=None,  # should never be called
     ),
     # Continuous methods
     interpolated_inverted_cdf=dict(
@@ -121,14 +121,12 @@
     lower=dict(
         get_virtual_index=lambda n, quantiles: np.floor(
             (n - 1) * quantiles).astype(np.intp),
-        fix_gamma=lambda gamma, _: gamma,
-        # should never be called, index dtype is int
+        fix_gamma=None,  # should never be called, index dtype is int
     ),
     higher=dict(
         get_virtual_index=lambda n, quantiles: np.ceil(
             (n - 1) * quantiles).astype(np.intp),
-        fix_gamma=lambda gamma, _: gamma,
-        # should never be called, index dtype is int
+        fix_gamma=None,  # should never be called, index dtype is int
     ),
     midpoint=dict(
         get_virtual_index=lambda n, quantiles: 0.5 * (
@@ -143,7 +141,7 @@
     nearest=dict(
         get_virtual_index=lambda n, quantiles: np.around(
             (n - 1) * quantiles).astype(np.intp),
-        fix_gamma=lambda gamma, _: gamma,
+        fix_gamma=None,
         # should never be called, index dtype is int
     ))
 
@@ -4887,15 +4885,19 @@ def _quantile(
         # Virtual because it is a floating point value, not an valid index.
         # The nearest neighbours are used for interpolation
         try:
-            method = _QuantileMethods[method]
+            method_props = _QuantileMethods[method]
         except KeyError:
             raise ValueError(
                 f"{method!r} is not a valid method. Use one of: "
                 f"{_QuantileMethods.keys()}") from None
-        virtual_indexes = method["get_virtual_index"](values_count, quantiles)
+        virtual_indexes = method_props["get_virtual_index"](values_count, quantiles)
         virtual_indexes = np.asanyarray(virtual_indexes)
 
-        if np.issubdtype(virtual_indexes.dtype, np.integer):
+        int_virtual_indices = np.issubdtype(virtual_indexes.dtype, np.integer)
+        supports_integers = method_props["fix_gamma"] is None
+        supports_integers = supports_integers or (method == 'linear' and int_virtual_indices)
+
+        if supports_integers:
             # No interpolation needed, take the points along axis
             if supports_nans:
                 # may contain nan, which would sort to the end
@@ -4927,7 +4929,7 @@ def _quantile(
             previous = arr[previous_indexes]
             next = arr[next_indexes]
             # --- Linear interpolation
-            gamma = _get_gamma(virtual_indexes, previous_indexes, method)
+            gamma = _get_gamma(virtual_indexes, previous_indexes, method_props)
             result_shape = virtual_indexes.shape + (1,) * (arr.ndim - 1)
             gamma = gamma.reshape(result_shape)
             result = _lerp(previous,
diff --git a/numpy/lib/tests/test_function_base.py b/numpy/lib/tests/test_function_base.py
index e3805b540f0d..e95426cd9390 100644
--- a/numpy/lib/tests/test_function_base.py
+++ b/numpy/lib/tests/test_function_base.py
@@ -3732,6 +3732,13 @@ def test_quantile_preserve_int_type(self, dtype):
                           method="nearest")
         assert res.dtype == dtype
 
+    @pytest.mark.parametrize("method", quantile_methods)
+    def test_q_zero_one(self, method):
+        # gh-24710
+        arr = [10, 11, 12]
+        quantile = np.quantile(arr, q = [0, 1], method=method)
+        assert_equal(quantile,  np.array([10, 12]))
+
     @pytest.mark.parametrize("method", quantile_methods)
     def test_quantile_monotonic(self, method):
         # GH 14685
@@ -3967,6 +3974,13 @@ def test_quantile_weights_raises_unsupported_methods(self, method):
         with pytest.raises(ValueError, match=msg):
             np.quantile(y, 0.5, weights=w, method=method)
 
+    def test_weibull_fraction(self):
+        arr = [Fraction(0, 1), Fraction(1, 10)]
+        quantile = np.quantile(arr, [0, ], method='weibull')
+        assert_equal(quantile, np.array(Fraction(0, 1)))
+        quantile = np.quantile(arr, [Fraction(1, 2)], method='weibull')
+        assert_equal(quantile, np.array(Fraction(1, 20)))
+
 
 class TestLerp:
     @hypothesis.given(t0=st.floats(allow_nan=False, allow_infinity=False,


Code review comments:
diff_hunk: @@ -4887,15 +4885,19 @@ def _quantile(
         # Virtual because it is a floating point value, not an valid index.
         # The nearest neighbours are used for interpolation
         try:
-            method = _QuantileMethods[method]
+            method_props = _QuantileMethods[method]
         except KeyError:
             raise ValueError(
                 f"{method!r} is not a valid method. Use one of: "
                 f"{_QuantileMethods.keys()}") from None
-        virtual_indexes = method["get_virtual_index"](values_count, quantiles)
+        virtual_indexes = method_props["get_virtual_index"](values_count, quantiles)
         virtual_indexes = np.asanyarray(virtual_indexes)
 
-        if np.issubdtype(virtual_indexes.dtype, np.integer):
+        int_virtual_indices = np.issubdtype(virtual_indexes.dtype, np.integer)
path: numpy/lib/_function_base_impl.py
text: @user1:
```suggestion
```

@user1:
Sorry, I thought you didn't use it at all?