<jupyter_start><jupyter_text>Multithreaded Cityblock distance matrix function with SciPy and Dask's delayed executionIn this notebook we implement a function to compute the cityblock distance matrix using `scipy.spatial.distance.cdist`. Althought this function is quite fast, it uses a single thread. In cases like this one, it might be convenient to implement a multithreaded version of the function by parallelicing the execution over chunks of data. We are going to use `dask.delayed` to do that.<jupyter_code>import numpy as np
from scipy.spatial.distance import cdist
from dask import compute, delayed, visualize
nsamples = 12000
nfeat = 50

x = 10. * np.random.random([nsamples, nfeat])<jupyter_output><empty_output><jupyter_text>Let's time the `cdist` function and look the `top` command.<jupyter_code># observe here that the funcion `cdist` used to get the cityblock distance
# is not multithreaded

%timeit cdist(x, x, 'cityblock')<jupyter_output>3.96 s ± 231 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)<jupyter_text>With the `top` command we see that `cdist` runs in a single thread. In such cases it could be quite simple write a distributed version of the function. We can do this very easily with `dask.delayed`!  Dask's async delayed executionA simple distributed version of `cdist` can be done as the following:  * Split the array of vectors into chunks. We can use `np.split(x, num_chunks)`  * Compute partial cityblock distance matrices of the complete array with respect to each of the chunks  * Concatenate the resulting list into a single cityblock distance matrix.Note that concatenation is not a fast operation, so probably we will have to continue improving our function.<jupyter_code># define the list of operations to be performed asynchronously
chunks = 12  # we choose on chunk for physical cpu (gpu partition)
partial_distances = [delayed(cdist)(x, xi, 'cityblock') for xi in np.split(x, chunks)]

# at this point nothing is executed
# visualize the copmutational graph
visualize(partial_distances)
cbdm_dask = delayed(np.concatenate)(partial_distances, axis=1)
# visualize the copmutational graph
visualize(cbdm_dask)<jupyter_output><empty_output><jupyter_text>Let's time the compute step and go to the command `top`. Now you can see that computation is executed in parallel resulting in a shorter execution time.<jupyter_code>%timeit compute(cbdm_dask, scheduler='threads')
# check that the resulting matrices are the same
cbdm = compute(cbdm_dask, scheduler='threads')[0]
np.abs(cbdm - cdist(x, x, 'cityblock')).max()<jupyter_output><empty_output><jupyter_text>A problem with this solution, as mentioned above, is that `np.concatenate` is not  a fast operation.Let's check how much time it takes without the concatenation part:<jupyter_code>%timeit compute(*partial_distances, scheduler='threads')<jupyter_output>1.64 s ± 132 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)<jupyter_text>Let's implement the whole thing as a single function:<jupyter_code>def cityblock_dask_concat(x, y, chunks):
    """Implementation using array concatenation"""
    partial_distances = [delayed(cdist)(x, xi, 'cityblock')
              for xi in np.split(x, chunks)]
    cbdm_dask = delayed(np.concatenate)(partial_distances, axis=1)
    return compute(cbdm_dask, scheduler='threads')
print(np.abs(cityblock_dask_concat(x, x, chunks) - cdist(x, x, 'cityblock')).max())<jupyter_output>0.0
StarcoderdataJupyterStructuredCleanDedup