Validity of using parallel in core._parallel_rolling_func

[earthgecko]

Hi

This is not a bug report per se rather asking question or just FYI.

This is a difficult issue to report because it is quite difficult to describe.

Note this is has only been tested/found with the stumpy.stump function,
however it is relevant in many other places too via core.preprocess_non_normalized

It appears that the use of parallel=True in core._parallel_rolling_func does
not have the desired effect in all circumstances.

Although there is no facility to invoke cache=True on the njit decorators in
stumpy, there are use cases where cache=True is an absolute requirement (I
shall update #699). So if modifications are made to allow for numba caching it
reveals that on the first run when stumpy.stump is called in a process/thread
it will compile and run fine (FYI with len(T_A) == 1004), the first run with the
njit compilation taking around 13 seconds or so. The numba cache files are
created and each subsequent call of the stumpy.stump function in the same
process/thread with the same T_A will run super fast with no issues in like
0.009974628977943212 seconds. All good.

The problem arises when another process is invoked and imports stumpy and loads
the numba cached jit files. At this point stumpy imports fine but when the
stumpy.stump function is called in jupyter the kernel dies and in a Python
terminal a Segmentation fault (core dumped) is encountered.

As we are all aware debugging with numba can be quite difficult and laborious :)

However, if parallel=True is simply commented out of the core._parallel_rolling_func
njit decorator, it works fine. I went through them all one by one (STRIPPED
down version), disabling both parallel and fastmath all individually and
tested, which eventually revealed core._parallel_rolling_func to be the
culprit.

@njit(
    # parallel=True,
    fastmath={"nsz", "arcp", "contract", "afn", "reassoc"},
    cache=config.STUMPY_NUMBA_CACHE,
)
def _parallel_rolling_func(a, w, func):
    """
    Compute the (embarrassingly parallel) rolling metric by applying a user defined
    function on a 1-D array

I am not certain of why this is the case, perhaps when caching is in play the
func parameter is the problem as what func is numba compiling it, how does it
know what func will be passed to it? However, my understanding of what can and
cannot be achieved with numba is not that verbose.

That said I did try to modify that to actually use the func name and that had
the same kernel died/Segmentation fault result so perhaps not.

For example:

def _parallel_rolling_func(a, w, func_name):
    ...
    ...

    l = a.shape[0] - w + 1
    out = np.empty(l)
    for i in prange(l):
        # out[i] = func(a[i : i + w])
        if func_name == 'np.ptp':
          out[i] = np.ptp(a[i : i + w])

# AND

def _rolling_isconstant(a, w):
    ...
    ...
    out = _parallel_rolling_func(a, w, 'np.ptp')

As I said a difficult one to describe. However although it may run during first
compilation, that the use of parallel in core._parallel_rolling_func is
breaking when compiled and cached does raise the question, is it valid to use
in this function?

That is one of the advantages of cache=True I find it really does validate
njit functions. Anything that is not correct will definitely always break when
the cached version is attempted to be loaded. It is quite handy for development
testing in that way I find (for me at least) :)

I must further caveat this with the fact that I am tested on a STRIPPED down
version, which I stripped down to trace this bug down. The version only has:

__init__.py - ONLY imports required in core.py, stump.py and aamp.py
core.py - ONLY functions required in stump.py and aamp.py
aamp.py
stump.py

Mainly because stump is all I currently need at the moment, but I must have
cache=True otherwise I would have to consider stepping back to
matrixprofile-foundation/matrixprofile but I would have maintain it myself seeing
as they have binned it now :) But it did load and run in under 1 second and
having to wait between 13 and 17 seconds to run stumpy.stump is not an option,
but 0.009974628977943212 seconds is awesome!

So I am not sure what to do with this? Other than advise you of the findings.

I have applied the change to a full v1.11.1 version but is modified to allow for
caching which allows for cache=True to be set on all core, stump and aamp
njit decorators and that it works too, as long as parallel is not passed. Once
again only tested with stumpy.stump.

[seanlaw]

@earthgecko Thank you for writing this up. Please allow me some time to review. In the meantime, considering the short length of your time series, I wonder if using the (nearly) pure numpy version would be "fast enough" for your use case.

from stumpy.stomp import _stomp
import numpy as np

T = np.random.rand(1004)
m = 50
mp = _stomp(T, m)

Note that we really only support direct usage of the functions provided in our user-facing API and _stomp is not one of them (it is left for historical reference). All non-user-facing API functions are not "supported" and are subject to change/removal without notice. In other words, use at your own discretion!

Note: _stomp does not provide a non-normalized version (just use scipy cdist in that case).

I will be back with more comments after I've had time to review your points above.

[seanlaw]

@earthgecko Before I am able to help, I have the following questions/comments:

1. Considering that the issue is resolved when you set parallel=False in core._parallel_rolling_func, this should have little to nothing to do with the stumpy.stump matrix profile computation except that core._parallel_rolling_func is being called within core.preprocess to compute core.rolling_isconstant. No other function calls core._parallel_rolling_func aside from core.rolling_isconstant and so it would be good to isolate the problem
2. Given the point above, can you please provide a minimum-reproducible-example that triggers the same error/segfault when you call core.rolling_isconstant
3. I'm still confused by what exactly do you mean when you say It appears that the use of parallel=True in core._parallel_rolling_func does not have the desired effect in all circumstances. Can you please elaborate? Again, I think a minimum-reproducible-example would really help.
4. I may be misunderstanding how you are doing things but if you are computing matrix profiles using STUMPY, which is already parallelized and leverages all of the threads on your machine, it makes very little sense to further use any Python's multiprocessing/multithreading as you would essentially be competing for the same resources that are being 100% consumed by STUMPY. Can you describe what it is that you are doing and how you are setting things up with regards to your compute environment and using cache=True?
5. Unfortunately, STUMPY is unlikely to provide cache=True support any time soon (if ever) as it would be extremely hard to maintain given our lack of expertise on numba caching. While I certainly appreciate your need for this functionality, this may be the (albeit, unfortunate) tradeoff in STUMPY that you'll have to decide on whether it meets your needs. However, your comment/feedback is welcome and we will do our best to take all suggestions into consideration when prioritizing our work. Thought, given our very limited time/resources (this is essentially 100% volunteer work), it will take time before we can even look at identifying the best path forward and that does not add additional strain to our currently limited time/resources.

[earthgecko]

Hi @seanlaw

Thanks for to reply, I shall try to answer your question as best I can.

1. Yes it is solely related to that function being called via the core.preprocess_non_normalized which is called
   in aamp.py and aamp is included in stump.py and referenced as

@core.non_normalized(aamp)
def stump(T_A, m, T_B=None, ignore_trivial=True, normalize=True, p=2.0, k=1):

The purpose or function of aamp there is beyond my understanding at this point.

2. minimum-reproducible-example

Create the problem.

In master in core.py core._parallel_rolling_func add , cache=True to the njit decorator

@njit(parallel=True, fastmath={"nsz", "arcp", "contract", "afn", "reassoc"}, cache=True)
def _parallel_rolling_func(a, w, func):

Start a jupyter notebook or python terminal:

from timeit import default_timer as timer
import numpy as np
import stumpy

x = np.linspace(start=-np.pi, stop=np.pi, num=1004)
ts = np.sin(x)

start = timer()
profile = stumpy.stump(ts, m=5)
print('took', (timer() - start), 'seconds')

# Runs fine, timer added to display compile time and for comparison when cache does work.

# Run again, super fast OK.
start = timer()
profile = stumpy.stump(ts, m=5)
print('took', (timer() - start), 'seconds')

# Exit the Python terminal or restart the notebook kernel

Now start a new Python terminal or in the notebook with the restarted kernel do the same again.
This time when stumpy.stump is called and the cached compiled jit functions are loaded the
kernel will die or the Python terminal will Segmentation fault.

To fix the problem, in master in core.py core._parallel_rolling_func leave , cache=True but
remove parallel=True from the njit decorator

@njit(fastmath={"nsz", "arcp", "contract", "afn", "reassoc"}, cache=True)
def _parallel_rolling_func(a, w, func):

Purge lib/python3.8/site_packages/stumpy/__pycache__ (whatever your local equivalent is) AND wherever numba created the cache files in my case
with Terminal - /home/$USER/.cache/numba/stumpy_f2c599dea3aece1bf68b6c639cb1190d48aa9219
and with Jupyter notebook - /home/$USER/.cache/ipython/stumpy_f2c599dea3aece1bf68b6c639cb1190d48aa9219
The UUID part _f2c599dea3aece1bf68b6c639cb1190d48aa9219 will/may be different on your system but it will
be stumpy_xxxxxx

Now run the same test again, this time the jit cache files will be created but they will work with the
second run when loaded and no Segmentation fault/kernel die will occur and the numba compilation
overhead will not be incurred and the first execution of stumpy.stump will take seconds.

3. As I said it is not per se a bug with the current stumpy implementation, as the function is
   not declared or tested with cache, but perhaps something is amiss with the use of parallel in
   this specific case? Is it a valid numba function, if it cannot be cached and reloaded?
   So the desired effect would be that there is an expectation that a valid numba function is
   able to be compiled, cached and work, whether it is intended to be cached or not, one should be
   able to compile and cache a numba function without encountering a segmentation fault,
   that suggests it is not valid. One of the advantages of trying to cache and reload is
   that it adds an additional check/test to the function, which is not a bad thing as it
   is quite difficult to test and debug with numba as we know.

4. Using with multiple processes. Consider if it were used on an adhoc basis, say we have
   some process or service that wants to analyse some timeseries on an adhoc basis. Let us
   say that it spawns a new process and assigns it work, check this timeseries with
   stumpy.stump and return the discords. That could be achieved via spawning a new process
   using multiprocessing or even requested from another remote service which spawns up a process
   to do the work. The point being that stumpy.stump could be executed in a single isolated
   process every time. It is here where jit caching comes into play. Perhaps we could spawn
   270 checks at once, without jit cache if the compilation overhead on the function has to be
   incurred by every process it makes it unfeasible to use. Isolating the analysis to a
   process when you are doing 1000s upon 1000s of analysis runs ensures isolation and
   reliability, allowing the parent or master process to assign, kill/terminate if exceeds
   max_execution time if there is a problem, error or unforseen circumstance, etc.
   So it is not about stumpy being parallelized, it is more about system design and stumpy
   not running in a long running process but on-demand.

5. Having added cache=True to all stumpy njit functions in core, aamp and stumpy in v1.11.1
   and it working, without parallel on _parallel_rolling_func of course :) It suggests that
   it may not be that difficult to do and in fact is quite trivial to achieve. Not that
   testing and maintainance thereof will necessarily be trivial, but we will find out as
   this is a must have feature that I need so I will have to maintain my own fork with
   cache support, so I will let you know how it goes :)

   UPDATE - however as you say in Add support for Numba jit/njit cache=True to reduce compilation overhead #699 perhaps I should just get a new cat :) It may be easier than maintaining a fork.

[seanlaw]

Thanks for to reply, I shall try to answer your question as best I can.

Thank you @earthgecko! This is helpful. I will find some time to see if I can reproduce the issue

The purpose or function of aamp there is beyond my understanding at this point.

So, aamp is simply the non-normalized equivalent of stump (the naming comes from the published AAMP paper), which computes the non-normalized (i.e., without applying z-normalization to each subsequence) Euclidean distance. Instead of requiring the user to know when to call stump vs aamp, we created a simple decorator function that simply detects the presence of the normalize argument and calls the stump function when normalize=True OR calls the aamp function when normalize=False. There are some other fancy things like cleaning up the other arguments but it allows us to keep the (very different) normalized vs non-normalized implementations separate while making that purposely opaque to the user.

[seanlaw]

@earthgecko I was able to reproduce the issue with a minimum-reproducible-example that does not rely on STUMPY:

# Inside of a Jupyter Notebook
import numpy as np
from numba import njit, prange
import time
import pathlib


# @njit(parallel=True, fastmath={"nsz", "arcp", "contract", "afn", "reassoc"})
@njit(parallel=True, fastmath={"nsz", "arcp", "contract", "afn", "reassoc"})
def _parallel_rolling_func(a, w, func):
    l = a.shape[0] - w + 1
    out = np.empty(l)
    for i in prange(l):
        out[i] = func(a[i : i + w])

    return out


@njit(cache=True, fastmath={"nsz", "arcp", "contract", "afn", "reassoc"})
def _rolling_isconstant(a, w):
    out = _parallel_rolling_func(a, w, np.ptp)

    return np.where(out == 0.0, True, False)


def rolling_isconstant(a, w):
    axis = a.ndim - 1
    return np.apply_along_axis(
        lambda a_row, w: _rolling_isconstant(a_row, w), axis=axis, arr=a, w=w
    )


def clear_numba_cache():
    [f.unlink() for f in pathlib.Path(".ipython/numba_cache").glob("*") if f.is_file()]

    
def print_numba_cache():
    [print(f) for f in pathlib.Path(".ipython/numba_cache").iterdir()]


T = np.random.rand(500_000_000)
m = 50

start = time.time()
out = rolling_isconstant(T, m)
print(time.time() - start)

print_numba_cache()
# clear_numba_cache()
# print_numba_cache()

I learned that while functions are first class citizens in Python and can be passed around to other functions, sadly, this currently causes issues with numba. As a result, I have fixed this in 45eb3f6 by removing _parallel_rolling_func altogether, which eliminates the passing of a function into a njit function.

Would you mind pulling the latest commit and checking if this solves your problem? Does the segfault/kernel restart go away?

[seanlaw]

In fact, you don't need to get/spawn a new cat. If you want to cache either stump (really, _stump, the z-normalized matrix profile) and/or aamp (really, _aamp, the non-normalized matrix profile) then all you need to do is:

import numpy as np
import stumpy
from stumpy.aamp import _aamp
from stumpy.stump import _stump
import time

_aamp.enable_caching()
_stump.enable_caching()

T = np.random.rand(1004)
m = 50

start = time.time()
out = stumpy.stump(T, m)
print(time.time() - start)

start = time.time()
out = stumpy.stump(T, m, normalize=False)
print(time.time() - start)

The first time you do this, they'll be slow but they will be cached in your ../site-packages/stumpy/__pycache__ and then the next time you call either stump or aamp, they won't be recompiled and therefore are much faster. Note that other njit functions that are being called within stump/aamp like core.preprocess_diagonal can also be cached in a similar way as well.

Finally, to force that cache to be cleared from site-packages, you can use these helper functions:

def clear_numba_cache():
    site_pkg_dir = site.getsitepackages()[0]
    numba_cache_dir = site_pkg_dir + "/stumpy/__pycache__"
    [f.unlink() for f in pathlib.Path(numba_cache_dir).glob("*nb*") if f.is_file()]

    
def print_numba_cache():
    site_pkg_dir = site.getsitepackages()[0]
    numba_cache_dir = site_pkg_dir + "/stumpy/__pycache__"
    [print(f) for f in pathlib.Path(numba_cache_dir).glob("*nb*") if f.is_file()]

where calling clear_numba_cache() will remove all cached functions from the stumpy directory within site-packages and stumpy and aamp (really, _stump and _aamp functions) will be slow again when called the first time within a given process.

Let me know if that makes sense and if I can help explain anything conceptually.

[earthgecko]

Hi @seanlaw

Thanks very much for the efforts, very much appreciated, especially on a Saturday! Apologies for the distraction to your weekend.

The only part of this that you need to know today is that I can confirm that 45eb3f6 has the desired effects and works 👍

Please do not spend any more time on it at this stage, you have done more than enough helping to clarify lots of things on this issue, which I appreciate is well outside the current scope of stumpy. Below is some feedback but can wait as it is just results of some tests and some thoughts thereabouts, so come back it at a later date (not a weekend day) ... if you can ;)

You are still reading? (I hope it is not Sunday 15 Jan 2023, if so close this tab)

Firstly I am not suggesting that implementing caching requires your attention in any way at this point in terms of requesting the ability be added to stumpy. These are just some experiments and notes that may be useful in some way in the future, even if just for my own reference as I move forward.

I have tested it directly and I also took the time to port those changes to the stripped down version I was/am testing with that has cache=True on all the njt decorators just to be certain and it works in that context too, so thanks a stack for working through this with me.

I was going to ask if _parallel_rolling_func was even needed to be declared as it's own function seeing as it is only used in the one place with np.ptp but I figured that perhaps you may have had an idea to use it somewhere else in the future. What really does surprise me is that the test I did yesterday using the func_name and a conditional did not have the same effect as the func is not being passed but is being declared directly in the njit function. Strange ...

Just some feedback regarding your suggestion relating to using the _stump.enable_caching() method, it does indeed create cache resources for stump._stump, e.g.:

ls -1 ~/.cache/ipython/stumpy_f2c599dea3aece1bf68b6c639cb1190d48aa9219/
stump._stump-249.py38.1.nbc
stump._stump-249.py38.nbi

However it does not create cache files for all the other njit functions that are used with stumpy.stump. When cache=True or say cache=config.STUMPY_NUMBA_CACHE_ENABLED numba will compile and cache all njit functions that are called in stumpy.stump, e.g.:

ls -1 ~/.cache/ipython/stumpy_f2c599dea3aece1bf68b6c639cb1190d48aa9219/
core._count_diagonal_ndist-864.py38.1.nbc
core._count_diagonal_ndist-864.py38.nbi
core._get_array_ranges-904.py38.1.nbc
core._get_array_ranges-904.py38.nbi
core._merge_topk_ρI-1075.py38.1.nbc
core._merge_topk_ρI-1075.py38.nbi
core._rolling_isconstant-1013.py38.1.nbc
core._rolling_isconstant-1013.py38.nbi
core._rolling_nanstd_1d-558.py38.1.nbc
core._rolling_nanstd_1d-558.py38.nbi
core._shift_insert_at_index-1151.py38.1.nbc
core._shift_insert_at_index-1151.py38.2.nbc
core._shift_insert_at_index-1151.py38.nbi
stump._compute_diagonal-13.py38.1.nbc
stump._compute_diagonal-13.py38.2.nbc
stump._compute_diagonal-13.py38.nbi
stump._stump-251.py38.1.nbc
stump._stump-251.py38.nbi

The first execution of stumpy.stump for both methods on the first run when jit compiles and caches is around 13 seconds.

The second execution of stumpy.stump with the cache files generated by the _stump.enable_caching() method takes around 1.8207929519703612 seconds
Whereas the second execution of stumpy.stump with the cache files generated by the njit cache=True method takes around 0.28026175301056355 seconds due to the cache files being present for all the other njit functions that are used.

I appreciate that this is somewhat semantics, however I do not think that the enable_caching() method is the best way to achieve performance gains of numba caching. It would be possible to achieve the same as nit cache=True method with the enable_caching() method but that would require declaring:

from stumpy.stump import _compute_diagonal, _stump
from core import (
  _count_diagonal_ndist, _get_array_ranges, _merge_topk_ρI, _rolling_isconstant,
  _rolling_nanstd_1d, _shift_insert_at_index
)

_compute_diagonal.enable_caching()
_stump.enable_caching()
_count_diagonal_ndist.enable_caching()
_get_array_ranges.enable_caching()
_merge_topk_ρI.enable_caching()
_rolling_isconstant.enable_caching()
_rolling_nanstd_1d.enable_caching()
_shift_insert_at_index.enable_caching()

This indeed does produce all the same numba cache files for all the functions called by stumpy.stump.

Granted, it is a solution without any changes to stumpy, however it is quite opaque and obfuscate for the user.
Further to this, to gain the numba cache performance gains the user would have to have in-depth knowledge of what functions where being called by a function to declare enable_caching() on them all and as an added bonus have to assess whether any stumpy changes or additions were made to any core or other related functions that they have to consider in their code when they update the stumpy library.

Just pointing out that although this enable_caching() method works it is probably not something to suggest as method of achieving numba caching as it is quite involved.

Additional to that the path in the helper functions can be variable depending on how both numba and stumpy were run. The NUMBA_CACHE_DIR is determined by numba and the environment in which things are being executed, on my Linux systems all numba cache files are written to ~/.cache/numba and I am not defining the NUMBA_CACHE_DIR envrionment variable specifically (Ubuntu 20.04 and Jupyter and on CentOS Stream 8 via a daemon Python process), so helper functions may not be that simple. I think generally if someone wanted to run numba caching the onus should mostly be on them to manage it, but informing the user of relevant points of consideration does help.

Although at the moment stumpy is not dealing with numba cache, a few questions have been raised and these are not totally invalid queries relating to the ability to use numba caching, there are some upsides to it and in some use cases these can be very important aspects.

All that said, I would suggest if you do consider adding any ability to use numba caching, just go down a cache=config.STUMPY_NUMBA_CACHE_ENABLED type of route and let numba do it all.

I am going to continue down that route myself for the time being as I am also going to port over from using matrixprofile/mass_ts to using the stumpy implementation (at some point after I determine how to convert those functions) so I may want numba caching in that too. So I will fork and try a cache=config.STUMPY_NUMBA_CACHE_ENABLED version of 45eb3f6109cc8588a48ff3ef9e44ae219b6c6681and I will set up a stumpy development environment and run all the tests, etc to see how it pans out and report back here on how that works out at a later date.

Enjoy your day (Sunday?) @seanlaw

[seanlaw]

Thanks @earthgecko, this is a labor of love. Now that I understand caching a little bit better, I've added a (untested, purely experimental) cache (see stumpy/cache.py) module that will automatically search for and cache ALL njit functions (please pull the latest commit). It can be called via:

import numpy as np
from stumpy import cache
import stumpy
import time

cache._enable()

T = np.random.rand(1000)
m = 50

start = time.time()
stumpy.stump(T, m)
print(time.time() - start)

start = time.time()
stumpy.stump(T, m)
print(time.time() - start)

You can also get a list of the cached functions via:

from stumpy import cache

for fname in cache._get_cache():
    print(fname)

and the cached functions (stored in site-packages/stumpy/__pycache__) can be cleared via:

from stumpy import cache

cache._clear()

Note that:

1. njit functions are only cached after being called once along a given execution path (i.e., all sub-functions are also cached)
2. even when all of the functions are cached, calling those functions for the first time still appears to incur a numba start-up/warm-up time (around 0.09 seconds on my M1 Macbook Air) and the second call to the function will still be much, much faster

Please feel free to ask any questions and provide any feedback. Again, this is purely experimental and should not be relied upon but, if it works, it's currently straightforward enough that I likely won't remove it (no guarantees though!). Use it at your own risk!

[earthgecko]

@seanlaw thanks I shall do 👍

[earthgecko]

To move to discussion 👍
Closing as resolved by #779 and 45eb3f6