Repository: sktime/sktime
Pull request: [BUG] Numba/np.median interaction raises error for large data sets run with n_jobs>1 issue #3582
#### Reference Issues/PRs

Fixes #3582

#### What does this implement/fix? Explain your changes.

#3582 identified a bug that occurs with DrCIF classifier when n_jobs > 1 and the dataset is large. This originated in the numbaversion of the median function in utils/numba/stats.py. The problem is that when the data is beyond some size threshold, Parallel uses mmap to write the data to local storage to avoid unnecessary data cloning. When this happens, a bug occurs in the interaction between numba and np.median caused by write permissions. See 
https://github.com/numba/numba/issues/3070

The np.median implementation in Numba doesn't support read-only arrays, and once the mmap condition in Parallel is triggered, the cached arrays are all read-only (and need to be to avoid thread interference). This triggers the error seen in #3582. It was not detected in testing since it only happens with data beyond a certain size. 

This solution simply clones the array before calling median if it is not writable. Alternatives considered are described in #3582. Profiling shows that the cloning does not introduce much of an overhead, and we think this is the safest solution until numba/numpy problem is fixed. 

   
#### Does your contribution introduce a new dependency? If yes, which one?

No

#### What should a reviewer concentrate their feedback on?


#### Any other comments?

Alternative solutions welcome. We could, for example, implement a local median function. However, this bug has been identified by two users, and even if temporary, this solution will move things forward.

#### PR checklist
<!--
Please go through the checklist below. Please feel free to remove points if they are not applicable.
-->

##### For all contributions
- [ ] I've added unit tests and made sure they pass locally.
- [x ] The PR title starts with either [ENH], [MNT], [DOC], or [BUG] indicating whether the PR topic is related to enhancement, maintenance, documentation, or bug.

<!--
Thanks for contributing!
-->


Problem statement:
[BUG] Numba error when setting n_jobs
**Describe the bug**
<!--
-->
Numba error when setting n_jobs

**To Reproduce**
<!--
-->
How can i upload my dataset? It has to be related to my dataset because it works fine with sktime datasets.

```python
estimators = [("oscillator", DrCIF(n_jobs=6), oscillator_columns)]
ensemble_model = ColumnEnsembleClassifier(estimators=estimators, remainder="passthrough")
ensemble_model.fit(exog, endog)
```

**Additional context**
<!--
-->
I get the following error when I fit the classifier and set n_jobs (to 6 for exemple). I have double check my dataset, everything seems fine (ive tried applying np.median to all entries and no problem with that).

```
numba.core.errors.TypingError: Failed in nopython mode pipeline (step: nopython frontend)
No implementation of function Function(<function median at 0x7fecff6e0160>) found for signature:
 
 >>> median(readonly array(float64, 1d, C))
 
There are 2 candidate implementations:
  - Of which 2 did not match due to:
  Overload in function 'np_median': File: numba/np/arraymath.py: Line 1447.
    With argument(s): '(readonly array(float64, 1d, C))':
   Rejected as the implementation raised a specific error:
     TypingError: Failed in nopython mode pipeline (step: nopython frontend)
   No implementation of function Function(<function _median_inner at 0x7fecdee57820>) found for signature:
    
    >>> _median_inner(readonly array(float64, 1d, C), int64)
    
   There are 2 candidate implementations:
         - Of which 2 did not match due to:
         Overload in function 'register_jitable.<locals>.wrap.<locals>.ov_wrap': File: numba/core/extending.py: Line 159.
           With argument(s): '(readonly array(float64, 1d, C), int64)':
          Rejected as the implementation raised a specific error:
            TypingError: Failed in nopython mode pipeline (step: nopython frontend)
          No implementation of function Function(<function _select_two at 0x7fecdee57700>) found for signature:
           
           >>> _select_two(readonly array(float64, 1d, C), int64, Literal[int](0), int64)
           
          There are 2 candidate implementations:
                - Of which 1 did not match due to:
                Overload in function 'register_jitable.<locals>.wrap.<locals>.ov_wrap': File: numba/core/extending.py: Line 159.
                  With argument(s): '(readonly array(float64, 1d, C), int64, int64, int64)':
                 Rejected as the implementation raised a specific error:
                   TypingError: Failed in nopython mode pipeline (step: nopython frontend)
                 No implementation of function Function(<function _partition_factory.<locals>._partition at 0x7fecdee57280>) found for signature:
                  
                  >>> _partition(readonly array(float64, 1d, C), int64, int64)
                  
                 There are 2 candidate implementations:
                       - Of which 2 did not match due to:
                       Overload in function 'register_jitable.<locals>.wrap.<locals>.ov_wrap': File: numba/core/extending.py: Line 159.
                         With argument(s): '(readonly array(float64, 1d, C), int64, int64)':
                        Rejected as the implementation raised a specific error:
                          TypingError: Failed in nopython mode pipeline (step: nopython frontend)
                        No implementation of function Function(<built-in function setitem>) found for signature:
                         
                         >>> setitem(readonly array(float64, 1d, C), int64, float64)
                         
                        There are 16 candidate implementations:
                              - Of which 14 did not match due to:
                              Overload of function 'setitem': File: <numerous>: Line N/A.
                                With argument(s): '(readonly array(float64, 1d, C), int64, float64)':
                               No match.
                              - Of which 2 did not match due to:
                              Overload in function 'SetItemBuffer.generic': File: numba/core/typing/arraydecl.py: Line 176.
                                With argument(s): '(readonly array(float64, 1d, C), int64, float64)':
                               Rejected as the implementation raised a specific error:
                                 NumbaTypeError: Cannot modify readonly array of type: readonly array(float64, 1d, C)
                          raised from /home/henry/anaconda3/lib/python3.9/site-packages/numba/core/typing/arraydecl.py:183
                        
                        During: typing of setitem at /home/henry/anaconda3/lib/python3.9/site-packages/numba/np/arraymath.py (1354)
                        
                        File "../../anaconda3/lib/python3.9/site-packages/numba/np/arraymath.py", line 1354:
                            def _partition(A, low, high):
                                <source elided>
                                if pivotimpl(A[mid], A[low]):
                                    A[low], A[mid] = A[mid], A[low]
                                    ^
                 
                   raised from /home/henry/anaconda3/lib/python3.9/site-packages/numba/core/typeinfer.py:1086
                 
                 During: resolving callee type: Function(<function _partition_factory.<locals>._partition at 0x7fecdee57280>)
                 During: typing of call at /home/henry/anaconda3/lib/python3.9/site-packages/numba/np/arraymath.py (1416)
                 
                 
                 File "../../anaconda3/lib/python3.9/site-packages/numba/np/arraymath.py", line 1416:
                 def _select_two(arry, k, low, high):
                     <source elided>
                         assert high > low  # by construction
                         i = _partition(arry, low, high)
                         ^
          
            raised from /home/henry/anaconda3/lib/python3.9/site-packages/numba/core/typeinfer.py:1086
                - Of which 1 did not match due to:
                Overload in function 'register_jitable.<locals>.wrap.<locals>.ov_wrap': File: numba/core/extending.py: Line 159.
                  With argument(s): '(readonly array(float64, 1d, C), int64, Literal[int](0), int64)':
                 Rejected as the implementation raised a specific error:
                   TypingError: Failed in nopython mode pipeline (step: nopython frontend)
                 No implementation of function Function(<function _partition_factory.<locals>._partition at 0x7fecdee57280>) found for signature:
                  
                  >>> _partition(readonly array(float64, 1d, C), Literal[int](0), int64)
                  
                 There are 2 candidate implementations:
                       - Of which 1 did not match due to:
                       Overload in function 'register_jitable.<locals>.wrap.<locals>.ov_wrap': File: numba/core/extending.py: Line 159.
                         With argument(s): '(readonly array(float64, 1d, C), int64, int64)':
                        Rejected as the implementation raised a specific error:
                          TypingError: Failed in nopython mode pipeline (step: nopython frontend)
                        No implementation of function Function(<built-in function setitem>) found for signature:
                         
                         >>> setitem(readonly array(float64, 1d, C), int64, float64)
                         
                        There are 16 candidate implementations:
                              - Of which 14 did not match due to:
                              Overload of function 'setitem': File: <numerous>: Line N/A.
                                With argument(s): '(readonly array(float64, 1d, C), int64, float64)':
                               No match.
                              - Of which 2 did not match due to:
                              Overload in function 'SetItemBuffer.generic': File: numba/core/typing/arraydecl.py: Line 176.
                                With argument(s): '(readonly array(float64, 1d, C), int64, float64)':
                               Rejected as the implementation raised a specific error:
                                 NumbaTypeError: Cannot modify readonly array of type: readonly array(float64, 1d, C)
                          raised from /home/henry/anaconda3/lib/python3.9/site-packages/numba/core/typing/arraydecl.py:183
                        
                        During: typing of setitem at /home/henry/anaconda3/lib/python3.9/site-packages/numba/np/arraymath.py (1354)
                        
                        File "../../anaconda3/lib/python3.9/site-packages/numba/np/arraymath.py", line 1354:
                            def _partition(A, low, high):
                                <source elided>
                                if pivotimpl(A[mid], A[low]):
                                    A[low], A[mid] = A[mid], A[low]
                                    ^
                 
                   raised from /home/henry/anaconda3/lib/python3.9/site-packages/numba/core/typeinfer.py:1086
                       - Of which 1 did not match due to:
                       Overload in function 'register_jitable.<locals>.wrap.<locals>.ov_wrap': File: numba/core/extending.py: Line 159.
                         With argument(s): '(readonly array(float64, 1d, C), Literal[int](0), int64)':
                        Rejected as the implementation raised a specific error:
                          TypingError: Failed in nopython mode pipeline (step: nopython frontend)
                        No implementation of function Function(<built-in function setitem>) found for signature:
                         
                         >>> setitem(readonly array(float64, 1d, C), Literal[int](0), float64)
                         
                        There are 16 candidate implementations:
                              - Of which 14 did not match due to:
                              Overload of function 'setitem': File: <numerous>: Line N/A.
                                With argument(s): '(readonly array(float64, 1d, C), int64, float64)':
                               No match.
                              - Of which 2 did not match due to:
                              Overload in function 'SetItemBuffer.generic': File: numba/core/typing/arraydecl.py: Line 176.
                                With argument(s): '(readonly array(float64, 1d, C), int64, float64)':
                               Rejected as the implementation raised a specific error:
                                 NumbaTypeError: Cannot modify readonly array of type: readonly array(float64, 1d, C)
                          raised from /home/henry/anaconda3/lib/python3.9/site-packages/numba/core/typing/arraydecl.py:183
                        
                        During: typing of setitem at /home/henry/anaconda3/lib/python3.9/site-packages/numba/np/arraymath.py (1354)
                        
                        File "../../anaconda3/lib/python3.9/site-packages/numba/np/arraymath.py", line 1354:
                            def _partition(A, low, high):
                                <source elided>
                                if pivotimpl(A[mid], A[low]):
                                    A[low], A[mid] = A[mid], A[low]
                                    ^
                 
                   raised from /home/henry/anaconda3/lib/python3.9/site-packages/numba/core/typeinfer.py:1086
                 
                 During: resolving callee type: Function(<function _partition_factory.<locals>._partition at 0x7fecdee57280>)
                 During: typing of call at /home/henry/anaconda3/lib/python3.9/site-packages/numba/np/arraymath.py (1416)
                 
                 
                 File "../../anaconda3/lib/python3.9/site-packages/numba/np/arraymath.py", line 1416:
                 def _select_two(arry, k, low, high):
                     <source elided>
                         assert high > low  # by construction
                         i = _partition(arry, low, high)
                         ^
          
            raised from /home/henry/anaconda3/lib/python3.9/site-packages/numba/core/typeinfer.py:1086
          
          During: resolving callee type: Function(<function _select_two at 0x7fecdee57700>)
          During: typing of call at /home/henry/anaconda3/lib/python3.9/site-packages/numba/np/arraymath.py (1441)
          
          
          File "../../anaconda3/lib/python3.9/site-packages/numba/np/arraymath.py", line 1441:
          def _median_inner(temp_arry, n):
              <source elided>
              if n & 1 == 0:
                  a, b = _select_two(temp_arry, half - 1, low, high)
                  ^
   
     raised from /home/henry/anaconda3/lib/python3.9/site-packages/numba/core/typeinfer.py:1086
   
   During: resolving callee type: Function(<function _median_inner at 0x7fecdee57820>)
   During: typing of call at /home/henry/anaconda3/lib/python3.9/site-packages/numba/np/arraymath.py (1457)
   
   
   File "../../anaconda3/lib/python3.9/site-packages/numba/np/arraymath.py", line 1457:
       def median_impl(a):
           <source elided>
           n = temp_arry.shape[0]
           return _median_inner(temp_arry, n)
           ^

  raised from /home/henry/anaconda3/lib/python3.9/site-packages/numba/core/typeinfer.py:1086

During: resolving callee type: Function(<function median at 0x7fecff6e0160>)
During: typing of call at /home/henry/anaconda3/lib/python3.9/site-packages/sktime/utils/numba/stats.py (17)


File "../../anaconda3/lib/python3.9/site-packages/sktime/utils/numba/stats.py", line 17:
def median(X):
    <source elided>
    """Numba median function for a single time series."""
    return np.median(X)
```

<!--
Please run the following code snippet and paste the output here:
 
from sktime import show_versions; show_versions()
-->

System:
    python: 3.9.12 (main, Apr  5 2022, 06:56:58)  [GCC 7.5.0]
executable: /home/henry/anaconda3/bin/python
   machine: Linux-5.15.0-50-generic-x86_64-with-glibc2.31

Python dependencies:
          pip: 21.2.4
   setuptools: 61.2.0
      sklearn: 1.0.2
       sktime: 0.13.4
  statsmodels: 0.13.2
        numpy: 1.21.5
        scipy: 1.7.3
       pandas: 1.4.2
   matplotlib: 3.5.1
       joblib: 1.1.0
        numba: 0.55.1
     pmdarima: None
      tsfresh: None

<!-- Thanks for contributing! -->

Hints from issue discussion:
Hi, I don't mind taking a look at this. If you are allowed to share your data (which you indicate towards in your issue) you can send it to me at m.middlehurst@uea.ac.uk, else some data characteristics such as the size would be helpful. How are you loading/formatting the data prior to fitting?
Could it be because the dataset is too big?

By the way, I'd like to ask will the model be trained exactly the same way if we fit it with successive chunks of the dataset vs. fitting directly the the whole dataset in one call function call?
(My question applies to these especially: TimeSeriesForestClassifier, SupervisedTimeSeriesForest, CanonicalIntervalForest, DrCIF, RandomIntervalSpectralEnsemble)
Calling fit will overwrite any previous fit calls. Probably some quirk with threading/n_jobs>1 and numba, but not taken a look yet.
Okay I have sent you the dataset and notebook to reproduce the issue by email. thanks!
Seems to be a numba issue related to readonly arrays: https://github.com/numba/numba/issues/3070#issue-337198698

The solution might be:

```python
x = np.array(...)
x.setflags(write=True)
```
thanks @patrickzib, we were just looking at this and were wondering if it were just np.median/numba related or a wider issue with interference with the disk writing used in mmap.  One solution (from @dguijo) was to set 
fit = Parallel(n_jobs=self._threads_to_use, mmap_mode='w+')

I'm just concerned that by making the mmap writable by all threads, they will interfere. This would happen if, for example, first process had not finished writing to the temp file before second looks to see if the mmap is there. It also requires changing a lot of files for a bug localised on median (it now seems, thanks to that thread). 

Another solution is to just not make median a numba operation. 
some possible solutions:
1. set the memory map to be writable everywhere.  mmap_mode='w+'. Dangerous and requires changes to a large number of files 
2. Dont use numba for finding the median until this is (hopefully) fixed. Some performance hit for what is not perhaps a common  use case: parallel run on large data (although I am doing this myself atm!).
3. Always set arrays to be writeable in the median function, as @patrickzib suggests x.setflags(write=True). Seems a bit dangerous and would not really meet the user expectations? 
4. Clone data in median if passed a non writable array (@MatthewMiddlehurst). e.g.
```
def median(X):
    if X.flags.writeable:
        return _median(X)
    else:
        X = X.copy()
        X.setflags(write=True)
        return _median(X)


@njit(fastmath=True, cache=True)
def _median(X):
    """Numba median function for a single time series."""
    return np.median(X)
``` 
Introduces a time and memory overhead, but safer than (3) 

I dont like (1) as a widespread solution. Perhaps we should profile (2), (3) and (4) and try quantify the hit, (4) seems safest
 
Just to confirm, replacing the `median` method in `sktime\utils\numba\stats.py` with the above code allows me to run `DrCIF` with `n_jobs=4` on the data provided by @henrydeclety. With the previous method, this crashed with the same error as the original issue.
(2) may just shift the problem elsewhere, (3) doesnt work, still crashes. I roughly profiled (4) for run time: its not actually faster than (2) but safer, so favoured solution is (4). I'm just running it on FaceDetection
Would this work (avoids cloning)?

```python
def median(X):
    if X.flags.writeable:
        return _median(X)
    else:
        X.setflags(write=True)
        median = _median(X)
        X.setflags(write=False)
        return median
```

Commit message: update comment in median function

Patch to review:
diff --git a/sktime/utils/numba/stats.py b/sktime/utils/numba/stats.py
index 0cab58511c1..6667939a42c 100644
--- a/sktime/utils/numba/stats.py
+++ b/sktime/utils/numba/stats.py
@@ -11,9 +11,27 @@ def mean(X):
     return np.mean(X)
 
 
-@njit(fastmath=True, cache=True)
 def median(X):
-    """Numba median function for a single time series."""
+    """Numba median function for a single time series.
+
+    This function clones the data if the array passed is not writable before finding
+    the median. This is because of a known bug in the interaction between numba and
+    numpy np.median function that causes an error with threaded jobs with large
+    datasets. See
+    https://github.com/numba/numba/issues/3070
+    and
+    https://github.com/sktime/sktime/issues/3582
+    """
+    if X.flags.writeable:
+        return _median(X)
+    else:
+        X = X.copy()
+        X.setflags(write=True)
+        return _median(X)
+
+
+@njit(fastmath=True, cache=True)
+def _median(X):
     return np.median(X)
 
 


Code review comments:
diff_hunk: @@ -11,9 +11,27 @@ def mean(X):
     return np.mean(X)
 
 
-@njit(fastmath=True, cache=True)
 def median(X):
-    """Numba median function for a single time series."""
+    """Numba median function for a single time series.
+
+    This function clones the data if the array passed is not writable before finding
+    the median. This is because of a known bug in the interaction between numba and
+    numpy np.median function that causes an error with threaded jobs with large
+    datasets. See
+    https://github.com/numba/numba/issues/3070
+    and
+    https://github.com/sktime/sktime/issues/3582
+    """
+    if X.flags.writeable:
+        return _median(X)
+    else:
+        X = X.copy()
+        X.setflags(write=True)
+        return _median(X)
+
+
+@njit(fastmath=True, cache=True)
+def _median(X):
     return np.median(X)
path: sktime/utils/numba/stats.py
text: @user2:
Is calling np.median() from within @user1 faster than calling np.median directly?

@author:
> Sorry for asking the obvious? Is calling numpy methods from within numba `@user1` faster than calling numpy directly?

I profiled both options, and the time is about the same. However, if we do not njit this function, then the problem will just occur elsewhere when np.median is called in a njit'ed method, thus shifting the problem elsewhere. @user3 and @user4 may want to comment

@user2:
Ok, but in this case there should be an njit tag in line 14? 

@author:
see above :) I am very much black boxing this, its all about the interaction with the Parallel threading structure of caching. There may be a neater solution here, I don't like the cloning, but this has the benefit of fixing the problem 

@user2:
I am a bit confused :)

As I understand it: To use median() from an njit'ed method, all sub-methods must be njit'ed, too.

As  `median(X)` is not njit'ed, we must not call this from an njit'ed method. But we cannot call `_median(X)` directly from an njit'ed method either, as it has the readonly problem?

@user4:
I dont think you can actually call the new function in other njit functions, which may be a bit of an oversight now that i think about it 🙂 

@user2:
> see above :) I am very much black boxing this, its all about the interaction with the Parallel threading structure of caching. There may be a neater solution here, I don't like the cloning, but this has the benefit of fixing the problem

Which method is calling median? Is this method njit'ed at all? E.g. does it need an njit'ed median at all?

@user4:
This is less about whether it is faster, and more about requiring a numba median function for other numba functions. Just removing the njit is kicking the can down the road to the next numba function that requires the median.

@user2:
>  a numba median function for other numba functions

Sure, but then all methods should be njit'ed? Or am I missing on something?

I do understand that it is not possible to add njit in line 14. But we need a different workaround in this case.

@author:
hmmm, I'm confused too need some input from @user3 maybe? The basic problem is we cannot call np.median directly from any njit'ed method without allowing this bug, so this non-njit wrapper is meant to mitigate that. It definitely runs, tried it on the two cases which previously crashed (mine and the data from =3582). Perhaps we are simply un-njit-ing the whole thing by doing it like this?  

I think the best solution may just be to implement quickselect (or whatever the best median finder is) in numba 

@user4:
yeah

@user2:
Would "@user5" work at median() ? Not @user1

@author:
> I dont think you can actually call the new function in other njit functions, which may be a bit of an oversight now that i think about it 🙂

ok, so it seems its only used here 
![image](https://user-images.githubusercontent.com/9594042/196400928-82a33afa-646a-4dc1-88c5-b3d2f71e4425.png)

@user4 maybe just refactor it out and abandon this method?


@author:
and use np.median directly? 

@user2:
Sounds ok to me.

`numpy.partition` could be a solution to not using median, too. It is based on the quicksort-/quickselect-idea, and can be used to return the median, too by passing the index of the median to it.

@user4:
Can just delete the functions and replace the usage for now as you say

diff_hunk: @@ -11,9 +11,27 @@ def mean(X):
     return np.mean(X)
 
 
-@njit(fastmath=True, cache=True)
 def median(X):
path: sktime/utils/numba/stats.py
text: @user2:
Why is there no @user1 here?

@author:
you cannot test X.flags.writeable within njit, 

    @user1(fastmath=True, cache=True)
    def median4(X):
        if X.flags.writeable:
            return _median(X)
        else:
            X = X.copy()
            X.setflags(write=True)
            return _median(X)

File "debug.py", line 145:
    def median4(X):
        if X.flags.writeable:
        ^

During: typing of get attribute at C:\code\sktime-estimator-evaluation\sktime_estimator_evaluation\debug.py (145)

File "debug.py", line 145:
    def median4(X):
        if X.flags.writeable:
        ^