39

In Python or NumPy, what is the best way to find out the first occurrence of a subarray?

For example, I have

a = [1, 2, 3, 4, 5, 6]
b = [2, 3, 4]

What is the fastest way (run-time-wise) to find out where b occurs in a? I understand for strings this is extremely easy, but what about for a list or numpy ndarray?

Thanks a lot!

[EDITED] I prefer the numpy solution, since from my experience numpy vectorization is much faster than Python list comprehension. Meanwhile, the big array is huge, so I don't want to convert it into a string; that will be (too) long.

1
  • Could you just convert the list to a string to make the comparison? x=''.join(str(x) for x in a) Then use the find method with the resulting strings? Or do they have to remain lists?
    – danem
    Commented Aug 17, 2011 at 22:57

10 Answers 10

26

I'm assuming you're looking for a numpy-specific solution, rather than a simple list comprehension or for loop. One straightforward approach is to use the rolling window technique to search for windows of the appropriate size.

This approach is simple, works correctly, and is much faster than any pure Python solution. It should be sufficient for many use cases. However, it is not the most efficient approach possible, for a number of reasons. For an approach that is more complicated, but asymptotically optimal in the expected case, see the numba-based rolling hash implementation in norok2's answer.

Here's the rolling_window function:

>>> def rolling_window(a, size):
...     shape = a.shape[:-1] + (a.shape[-1] - size + 1, size)
...     strides = a.strides + (a. strides[-1],)
...     return numpy.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)
... 

Then you could do something like

>>> a = numpy.arange(10)
>>> numpy.random.shuffle(a)
>>> a
array([7, 3, 6, 8, 4, 0, 9, 2, 1, 5])
>>> rolling_window(a, 3) == [8, 4, 0]
array([[False, False, False],
       [False, False, False],
       [False, False, False],
       [ True,  True,  True],
       [False, False, False],
       [False, False, False],
       [False, False, False],
       [False, False, False]], dtype=bool)

To make this really useful, you'd have to reduce it along axis 1 using all:

>>> numpy.all(rolling_window(a, 3) == [8, 4, 0], axis=1)
array([False, False, False,  True, False, False, False, False], dtype=bool)

Then you could use that however you'd use a boolean array. A simple way to get the index out:

>>> bool_indices = numpy.all(rolling_window(a, 3) == [8, 4, 0], axis=1)
>>> numpy.mgrid[0:len(bool_indices)][bool_indices]
array([3])

For lists you could adapt one of these rolling window iterators to use a similar approach.

For very large arrays and subarrays, you could save memory like this:

>>> windows = rolling_window(a, 3)
>>> sub = [8, 4, 0]
>>> hits = numpy.ones((len(a) - len(sub) + 1,), dtype=bool)
>>> for i, x in enumerate(sub):
...     hits &= numpy.in1d(windows[:,i], [x])
... 
>>> hits
array([False, False, False,  True, False, False, False, False], dtype=bool)
>>> hits.nonzero()
(array([3]),)

On the other hand, this will probably be somewhat slower.

2
  • 1
    The problem with this approach is that ,while the return of rolling_window doesn't require any new memory, and reuses that of the original array, when doing the == operation you instantiate a new boolean array that is size times the full size of your original array. If the array is big enough, this can kill performance big time.
    – Jaime
    Commented Dec 19, 2013 at 17:09
  • 1
    You're right, this is not an asymptotically optimal solution. However, it strikes a nice balance between simplicity and efficiency — it's straightforward, correct, and faster than any pure-Python approach by orders of magnitude. For those who need a provably optimal solution, norok2's very detailed answer has several candidates, including a numba-based rolling hash approach that is asymptotically optimal in the expected case.
    – senderle
    Commented Mar 17, 2021 at 13:56
26

(EDITED to include a deeper discussion, better code and more benchmarks)


Summary

For raw speed and efficiency, one can use a Cython or Numba accelerated version (when the input is a Python sequence or a NumPy array, respectively) of one of the classical algorithms.

The recommended approaches are:

  • find_kmp_cy() for Python sequences (list, tuple, etc.)
  • find_kmp_nb() for NumPy arrays

Other efficient approaches, are find_rk_cy() and find_rk_nb() which, are more memory efficient but are not guaranteed to run in linear time.

If Cython / Numba are not available, again both find_kmp() and find_rk() are a good all-around solution for most use cases, although in the average case and for Python sequences, the naïve approach, in some form, notably find_pivot(), may be faster. For NumPy arrays, find_conv() (from @Jaime answer) outperforms any non-accelerated naïve approach.

(Full code is below, and here and there.)


Theory

This is a classical problem in computer science that goes by the name of string-searching or string matching problem. The naive approach, based on two nested loops, has a computational complexity of O(n + m) on average, but worst case is O(n m). Over the years, a number of alternative approaches have been developed which guarantee a better worst case performances.

Of the classical algorithms, the ones that can be best suited to generic sequences (since they do not rely on an alphabet) are:

This last algorithm relies on the computation of a rolling hash for its efficiency and therefore may require some additional knowledge of the input for optimal performance. Eventually, it is best suited for homogeneous data, like for example numeric arrays. A notable example of numeric arrays in Python is, of course, NumPy arrays.

Remarks

  • The naïve algorithm, by being so simple, lends itself to different implementations with various degrees of run-time speed in Python.
  • The other algorithms are less flexible in what can be optimized via language tricks.
  • Explicit looping in Python may be a speed bottleneck and several tricks can be used to perform the looping outside of the interpreter.
  • Cython is especially good at speeding up explicit loops for generic Python code.
  • Numba is especially good at speeding up explicit loops on NumPy arrays.
  • This is an excellent use-case for generators, so all the code will be using those instead of regular functions.

Python Sequences (list, tuple, etc.)

Based on the Naïve Algorithm

  • find_loop(), find_loop_cy() and find_loop_nb() which are the explicit-loop only implementation in pure Python, Cython and with Numba JITing respectively. Note the forceobj=True in the Numba version, which is required because we are using Python object inputs.
def find_loop(seq, subseq):
    n = len(seq)
    m = len(subseq)
    for i in range(n - m + 1):
        found = True
        for j in range(m):
            if seq[i + j] != subseq[j]:
                found = False
                break
        if found:
            yield i
%%cython -c-O3 -c-march=native -a
#cython: language_level=3, boundscheck=False, wraparound=False, initializedcheck=False, cdivision=True, infer_types=True


def find_loop_cy(seq, subseq):
    cdef Py_ssize_t n = len(seq)
    cdef Py_ssize_t m = len(subseq)
    for i in range(n - m + 1):
        found = True
        for j in range(m):
            if seq[i + j] != subseq[j]:
                found = False
                break
        if found:
            yield i
find_loop_nb = nb.jit(find_loop, forceobj=True)
find_loop_nb.__name__ = 'find_loop_nb'
  • find_all() replaces the inner loop with all() on a comprehension generator
def find_all(seq, subseq):
    n = len(seq)
    m = len(subseq)
    for i in range(n - m + 1):
        if all(seq[i + j] == subseq[j] for j in range(m)):
            yield i
  • find_slice() replaces the inner loop with direct comparison == after slicing []
def find_slice(seq, subseq):
    n = len(seq)
    m = len(subseq)
    for i in range(n - m + 1):
        if seq[i:i + m] == subseq:
            yield i
  • find_mix() and find_mix2() replaces the inner loop with direct comparison == after slicing [] but includes one or two additional short-circuiting on the first (and last) character which may be faster because slicing with an int is much faster than slicing with a slice().
def find_mix(seq, subseq):
    n = len(seq)
    m = len(subseq)
    for i in range(n - m + 1):
        if seq[i] == subseq[0] and seq[i:i + m] == subseq:
            yield i
def find_mix2(seq, subseq):
    n = len(seq)
    m = len(subseq)
    for i in range(n - m + 1):
        if seq[i] == subseq[0] and seq[i + m - 1] == subseq[m - 1] \
                and seq[i:i + m] == subseq:
            yield i
  • find_pivot() and find_pivot2() replace the outer loop with multiple .index() call using the first item of the sub-sequence, while using slicing for the inner loop, eventually with additional short-circuiting on the last item (the first matches by construction). The multiple .index() calls are wrapped in a index_all() generator (which may be useful on its own).
def index_all(seq, item, start=0, stop=-1):
    try:
        n = len(seq)
        if n > 0:
            start %= n
            stop %= n
            i = start
            while True:
                i = seq.index(item, i)
                if i <= stop:
                    yield i
                    i += 1
                else:
                    return
        else:
            return
    except ValueError:
        pass


def find_pivot(seq, subseq):
    n = len(seq)
    m = len(subseq)
    if m > n:
        return
    for i in index_all(seq, subseq[0], 0, n - m):
        if seq[i:i + m] == subseq:
            yield i
def find_pivot2(seq, subseq):
    n = len(seq)
    m = len(subseq)
    if m > n:
        return
    for i in index_all(seq, subseq[0], 0, n - m):
        if seq[i + m - 1] == subseq[m - 1] and seq[i:i + m] == subseq:
            yield i

Based on Knuth–Morris–Pratt (KMP) Algorithm

  • find_kmp() is a plain Python implementation of the algorithm. Since there is no simple looping or places where one could use slicing with a slice(), there is not much to be done for optimization, except using Cython (Numba would require again forceobj=True which would lead to slow code).
def find_kmp(seq, subseq):
    n = len(seq)
    m = len(subseq)
    # : compute offsets
    offsets = [0] * m
    j = 1
    k = 0
    while j < m: 
        if subseq[j] == subseq[k]: 
            k += 1
            offsets[j] = k
            j += 1
        else: 
            if k != 0: 
                k = offsets[k - 1] 
            else: 
                offsets[j] = 0
                j += 1
    # : find matches
    i = j = 0
    while i < n: 
        if seq[i] == subseq[j]: 
            i += 1
            j += 1
        if j == m:
            yield i - j
            j = offsets[j - 1] 
        elif i < n and seq[i] != subseq[j]: 
            if j != 0: 
                j = offsets[j - 1] 
            else: 
                i += 1
  • find_kmp_cy() is Cython implementation of the algorithm where the indices use C int data type, which result in much faster code.
%%cython -c-O3 -c-march=native -a
#cython: language_level=3, boundscheck=False, wraparound=False, initializedcheck=False, cdivision=True, infer_types=True


def find_kmp_cy(seq, subseq):
    cdef Py_ssize_t n = len(seq)
    cdef Py_ssize_t m = len(subseq)
    # : compute offsets
    offsets = [0] * m
    cdef Py_ssize_t j = 1
    cdef Py_ssize_t k = 0
    while j < m: 
        if subseq[j] == subseq[k]: 
            k += 1
            offsets[j] = k
            j += 1
        else: 
            if k != 0: 
                k = offsets[k - 1] 
            else: 
                offsets[j] = 0
                j += 1
    # : find matches
    cdef Py_ssize_t i = 0
    j = 0
    while i < n: 
        if seq[i] == subseq[j]: 
            i += 1
            j += 1
        if j == m:
            yield i - j
            j = offsets[j - 1] 
        elif i < n and seq[i] != subseq[j]: 
            if j != 0: 
                j = offsets[j - 1] 
            else: 
                i += 1

Based on Rabin-Karp (RK) Algorithm

  • find_rk() is a pure Python implementation, which relies on Python's hash() for the computation (and comparison) of the hash. Such hash is made rolling by mean of a simple sum(). The roll-over is then computed from the previous hash by subtracting the result of hash() on the just visited item seq[i - 1] and adding up the result of hash() on the newly considered item seq[i + m - 1].
def find_rk(seq, subseq):
    n = len(seq)
    m = len(subseq)
    if seq[:m] == subseq:
        yield 0
    hash_subseq = sum(hash(x) for x in subseq)  # compute hash
    curr_hash = sum(hash(x) for x in seq[:m])  # compute hash
    for i in range(1, n - m + 1):
        curr_hash += hash(seq[i + m - 1]) - hash(seq[i - 1])   # update hash
        if hash_subseq == curr_hash and seq[i:i + m] == subseq:
            yield i
  • find_rk_cy() is Cython implementation of the algorithm where the indices use the appropriate C data type, which results in much faster code. Note that hash() truncates "the return value based on the bit width of the host machine."
%%cython -c-O3 -c-march=native -a
#cython: language_level=3, boundscheck=False, wraparound=False, initializedcheck=False, cdivision=True, infer_types=True


def find_rk_cy(seq, subseq):
    cdef Py_ssize_t n = len(seq)
    cdef Py_ssize_t m = len(subseq)
    if seq[:m] == subseq:
        yield 0
    cdef Py_ssize_t hash_subseq = sum(hash(x) for x in subseq)  # compute hash
    cdef Py_ssize_t curr_hash = sum(hash(x) for x in seq[:m])  # compute hash
    cdef Py_ssize_t old_item, new_item
    for i in range(1, n - m + 1):
        old_item = hash(seq[i - 1])
        new_item = hash(seq[i + m - 1])
        curr_hash += new_item - old_item  # update hash
        if hash_subseq == curr_hash and seq[i:i + m] == subseq:
            yield i

Benchmarks

The above functions are evaluated on two inputs:

  • random inputs
def gen_input(n, k=2):
    return tuple(random.randint(0, k - 1) for _ in range(n))
  • (almost) worst inputs for the naïve algorithm
def gen_input_worst(n, k=-2):
    result = [0] * n
    result[k] = 1
    return tuple(result)

The subseq has fixed size (32). Since there are so many alternatives, two separate grouping have been done and some solutions with very small variations and almost identical timings have been omitted (i.e. find_mix2() and find_pivot2()). For each group both inputs are tested. For each benchmark the full plot and a zoom on the fastest approach is provided.

Naïve on Random

bm_full_naive_random bm_zoom_naive_random

Naïve on Worst

bm_full_naive_worst bm_zoom_naive_worst

Other on Random

bm_full_other_random bm_zoom_other_random

Other on Worst

bm_full_other_worst bm_zoom_other_worst

(Full code is available here.)


NumPy Arrays

Based on the Naïve Algorithm

  • find_loop(), find_loop_cy() and find_loop_nb() which are the explicit-loop only implementation in pure Python, Cython and with Numba JITing respectively. The code for the first two are the same as above and hence omitted. find_loop_nb() now enjoys fast JIT compilation. The inner loop has been written in a separate function because it can then be reused for find_rk_nb() (calling Numba functions inside Numba functions does not incur in the function call penalty typical of Python).
@nb.jit
def _is_equal_nb(seq, subseq, m, i):
    for j in range(m):
        if seq[i + j] != subseq[j]:
            return False
    return True


@nb.jit
def find_loop_nb(seq, subseq):
    n = len(seq)
    m = len(subseq)
    for i in range(n - m + 1):
        if _is_equal_nb(seq, subseq, m, i):
            yield i
  • find_all() is the same as above, while find_slice(), find_mix() and find_mix2() are almost identical to the above, the only difference is that seq[i:i + m] == subseq is now the argument of np.all(): np.all(seq[i:i + m] == subseq).

  • find_pivot() and find_pivot2() share the same ideas as above, except that now uses np.where() instead of index_all() and the need for enclosing the array equality inside an np.all() call.

def find_pivot(seq, subseq):
    n = len(seq)
    m = len(subseq)
    if m > n:
        return
    max_i = n - m
    for i in np.where(seq == subseq[0])[0]:
        if i > max_i:
            return
        elif np.all(seq[i:i + m] == subseq):
            yield i


def find_pivot2(seq, subseq):
    n = len(seq)
    m = len(subseq)
    if m > n:
        return
    max_i = n - m
    for i in np.where(seq == subseq[0])[0]:
        if i > max_i:
            return
        elif seq[i + m - 1] == subseq[m - 1] \
                and np.all(seq[i:i + m] == subseq):
            yield i
  • find_rolling() express the looping via a rolling window and the matching is checked with np.all(). This vectorizes all the looping at the expenses of creating large temporary objects, while still substantially appling the naïve algorithm. (The approach is from @senderle answer).
def rolling_window(arr, size):
    shape = arr.shape[:-1] + (arr.shape[-1] - size + 1, size)
    strides = arr.strides + (arr.strides[-1],)
    return np.lib.stride_tricks.as_strided(arr, shape=shape, strides=strides)


def find_rolling(seq, subseq):
    bool_indices = np.all(rolling_window(seq, len(subseq)) == subseq, axis=1)
    yield from np.mgrid[0:len(bool_indices)][bool_indices]
  • find_rolling2() is a slightly more memory efficient variation of the above, where the vectorization is only partial and one explicit looping (along the expected shortest dimension -- the length of subseq) is kept. (The approach is also from @senderle answer).
def find_rolling2(seq, subseq):
    windows = rolling_window(seq, len(subseq))
    hits = np.ones((len(seq) - len(subseq) + 1,), dtype=bool)
    for i, x in enumerate(subseq):
        hits &= np.in1d(windows[:, i], [x])
    yield from hits.nonzero()[0]

Based on Knuth–Morris–Pratt (KMP) Algorithm

  • find_kmp() is the same as above, while find_kmp_nb() is a straightforward JIT-compilation of that.
find_kmp_nb = nb.jit(find_kmp)
find_kmp_nb.__name__ = 'find_kmp_nb'

Based on Rabin-Karp (RK) Algorithm

  • find_rk() is the same as the above, except that again seq[i:i + m] == subseq is enclosed in an np.all() call.

  • find_rk_nb() is the Numba accelerated version of the above. Uses _is_equal_nb() defined earlier to definitively determine a match, while for the hashing, it uses a Numba accelerated sum_hash_nb() function whose definition is pretty straightforward.

@nb.jit
def sum_hash_nb(arr):
    result = 0
    for x in arr:
        result += hash(x)
    return result


@nb.jit
def find_rk_nb(seq, subseq):
    n = len(seq)
    m = len(subseq)
    if _is_equal_nb(seq, subseq, m, 0):
        yield 0
    hash_subseq = sum_hash_nb(subseq)  # compute hash
    curr_hash = sum_hash_nb(seq[:m])  # compute hash
    for i in range(1, n - m + 1):
        curr_hash += hash(seq[i + m - 1]) - hash(seq[i - 1])  # update hash
        if hash_subseq == curr_hash and _is_equal_nb(seq, subseq, m, i):
            yield i
  • find_conv() uses a pseudo Rabin-Karp method, where initial candidates are hashed using the np.dot() product and located on the convolution between seq and subseq with np.where(). The approach is pseudo because, while it still uses hashing to identify probable candidates, it is may not be regarded as a rolling hash (it depends on the actual implementation of np.correlate()). Also, it needs to create a temporary array the size of the input. (The approach is from @Jaime answer).
def find_conv(seq, subseq):
    target = np.dot(subseq, subseq)
    candidates = np.where(np.correlate(seq, subseq, mode='valid') == target)[0]
    check = candidates[:, np.newaxis] + np.arange(len(subseq))
    mask = np.all((np.take(seq, check) == subseq), axis=-1)
    yield from candidates[mask]

Benchmarks

Like before, the above functions are evaluated on two inputs:

  • random inputs
def gen_input(n, k=2):
    return np.random.randint(0, k, n)
  • (almost) worst inputs for the naïve algorithm
def gen_input_worst(n, k=-2):
    result = np.zeros(n, dtype=int)
    result[k] = 1
    return result

The subseq has fixed size (32). This plots follow the same scheme as before, summarized below for convenience.

Since there are so many alternatives, two separate grouping have been done and some solutions with very small variations and almost identical timings have been omitted (i.e. find_mix2() and find_pivot2()). For each group both inputs are tested. For each benchmark the full plot and a zoom on the fastest approach is provided.

Naïve on Random

bm_full_naive_random bm_zoom_naive_random

Naïve on Worst

bm_full_naive_worst bm_zoom_naive_worst

Other on Random

bm_full_other_random bm_zoom_other_random

Other on Worst

bm_full_other_worst bm_zoom_other_worst

(Full code is available here.)

1
  • Thanks for running all these tests! I linked to this from my answer. I still like mine because it's reasonably fast, easy to reason about, and doesn't add any dependencies. But for people who need a truly optimal solution, this is great. I disagree that the KMP approach is the best, though. For the vast majority of practical use cases, RK is faster, and very few people will truly need the worst-case guarantees that KMP provides.
    – senderle
    Commented Mar 17, 2021 at 15:05
25

The following code should work:

[x for x in xrange(len(a)) if a[x:x+len(b)] == b]

Returns the index at which the pattern starts.

3
  • 2
    This might not be the fastest solution, but +1 for the simplest answer. This might fit the needs of many users, especially if numpy is not available.
    – David
    Commented Jun 18, 2014 at 16:31
  • 4
    In Python 3 use range instead of xrange.
    – Samoth
    Commented Mar 2, 2018 at 23:07
  • 2
    For improved performance, you could replace len(a) with len(a) - len(b) + 1
    – norok2
    Commented Mar 23, 2020 at 23:49
20

A convolution based approach, that should be more memory efficient than the stride_tricks based approach:

def find_subsequence(seq, subseq):
    target = np.dot(subseq, subseq)
    candidates = np.where(np.correlate(seq,
                                       subseq, mode='valid') == target)[0]
    # some of the candidates entries may be false positives, double check
    check = candidates[:, np.newaxis] + np.arange(len(subseq))
    mask = np.all((np.take(seq, check) == subseq), axis=-1)
    return candidates[mask]

With really big arrays it may not be possible to use a stride_tricks approach, but this one still works:

haystack = np.random.randint(1000, size=(1e6))
needle = np.random.randint(1000, size=(100,))
# Hide 10 needles in the haystack
place = np.random.randint(1e6 - 100 + 1, size=10)
for idx in place:
    haystack[idx:idx+100] = needle

In [3]: find_subsequence(haystack, needle)
Out[3]: 
array([253824, 321497, 414169, 456777, 635055, 879149, 884282, 954848,
       961100, 973481], dtype=int64)

In [4]: np.all(np.sort(place) == find_subsequence(haystack, needle))
Out[4]: True

In [5]: %timeit find_subsequence(haystack, needle)
10 loops, best of 3: 79.2 ms per loop
2
  • While I really like this approach, I should note that in general finding candidates by l2 norm is not better than finding a particular symbol from needle. But after a small modification by computing dot product with randomized pattern of the same length as needle, this method will be just awesome.
    – Alleo
    Commented Apr 27, 2015 at 11:56
  • @Alleo I do not see what the problem with using the actual subseq is per se. The problem I see is that you may have more collisions if the subseq has duplicates or zeros, but a random sequence may have the same problem.
    – norok2
    Commented Jul 27, 2022 at 17:35
15

you can call tostring() method to convert an array to string, and then you can use fast string search. this method maybe faster when you have many subarray to check.

import numpy as np

a = np.array([1,2,3,4,5,6])
b = np.array([2,3,4])
print a.tostring().index(b.tostring())//a.itemsize
2
  • 1
    this solution is very quick and elegant, thanks a lot! Slightly related, I had a project grabbing np arrays of about 1e8 elements from C++ using a SWIG wrapper and the array creation was very slow. Working with them as strings boosted performance into real-time
    – fr_andres
    Commented Jun 7, 2019 at 19:52
  • Method is incorrect though. See np.array([0, 1]).tostring().index(np.array([256]).tostring()) Commented May 29, 2021 at 17:29
2

Another try, but I'm sure there is more pythonic & efficent way to do that ...

def array_match(a, b):
    for i in xrange(0, len(a)-len(b)+1):
        if a[i:i+len(b)] == b:
            return i
    return None
a = [1, 2, 3, 4, 5, 6]
b = [2, 3, 4]

print array_match(a,b)
1

(This first answer was not in scope of the question, as cdhowie mentionned)

set(a) & set(b) == set(b)
3
  • Two problems: This would also match [1, 3, 2, 4, 5, 6] (sets are not ordered; arrays are), and it doesn't report the location of the match (which should be index 1).
    – cdhowie
    Commented Aug 17, 2011 at 22:28
  • Yeah my bad, answered too quickly :-/
    – Stéphane
    Commented Aug 17, 2011 at 22:37
  • You can simplify your code a bit by replacing first_occurence=i with return i, and return first_occurence with return None.
    – Nayuki
    Commented Aug 17, 2011 at 23:06
1

Here is a rather straight-forward option:

def first_subarray(full_array, sub_array):
    n = len(full_array)
    k = len(sub_array)
    matches = np.argwhere([np.all(full_array[start_ix:start_ix+k] == sub_array) 
                   for start_ix in range(0, n-k+1)])
    return matches[0]

Then using the original a, b vectors we get:

a = [1, 2, 3, 4, 5, 6]
b = [2, 3, 4]
first_subarray(a, b)
Out[44]: 
array([1], dtype=int64)
1
  • You would probably add in some logic to take care of cases where there are no matches... Commented Apr 2, 2017 at 9:56
1
a = [1, 2, 3, 4, 5, 6]
b = [2, 3, 4]

np.concatenate((np.all(np.array([a[i:len(a)-(len(b)-1-i)] for i in range(len(b))]).T == b, axis = 1), np.full((len(b)-1), False)))

array([False,  True, False, False, False, False])

UPD: for relatively small subarrays this algorithm is faster than with rolling window, and the output is an array of the same size as a

the same idea:

np.all(np.lib.stride_tricks.sliding_window_view(a, len(b)) == b, axis = 1)
2
  • 1
    Your answer could be improved with additional supporting information. Please edit to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center. Commented Sep 28, 2023 at 14:39
  • Thank you for your interest in contributing to the Stack Overflow community. This question already has quite a few answers—including one that has been extensively validated by the community. Are you certain your approach hasn’t been given previously? If so, it would be useful to explain how your approach is different, under what circumstances your approach might be preferred, and/or why you think the previous answers aren’t sufficient. Can you kindly edit your answer to offer an explanation? Commented Sep 28, 2023 at 18:16
0

Quick comparison of three of the proposed solutions (average time of 100 iteration for randomly created vectors.):

import time
import collections
import numpy as np


def function_1(seq, sub):
    # direct comparison
    seq = list(seq)
    sub = list(sub)
    return [i for i in range(len(seq) - len(sub)) if seq[i:i+len(sub)] == sub]

def function_2(seq, sub):
    # Jamie's solution
    target = np.dot(sub, sub)
    candidates = np.where(np.correlate(seq, sub, mode='valid') == target)[0]
    check = candidates[:, np.newaxis] + np.arange(len(sub))
    mask = np.all((np.take(seq, check) == sub), axis=-1)
    return candidates[mask]

def function_3(seq, sub):
    # HYRY solution
    return seq.tostring().index(sub.tostring())//seq.itemsize


# --- assessment time performance
N = 100

seq = np.random.choice([0, 1, 2, 3, 4, 5, 6], 3000)
sub = np.array([1, 2, 3])

tim = collections.OrderedDict()
tim.update({function_1: 0.})
tim.update({function_2: 0.})
tim.update({function_3: 0.})

for function in tim.keys():
    for _ in range(N):
        seq = np.random.choice([0, 1, 2, 3, 4], 3000)
        sub = np.array([1, 2, 3])
        start = time.time()
        function(seq, sub)
        end = time.time()
        tim[function] += end - start

timer_dict = collections.OrderedDict()
for key, val in tim.items():
    timer_dict.update({key.__name__: val / N})

print(timer_dict)

Which would result (on my old machine) in:

OrderedDict([
('function_1', 0.0008518099784851074), 
('function_2', 8.157730102539063e-05), 
('function_3', 6.124973297119141e-06)
])
-1

First, convert the list to string.

a = ''.join(str(i) for i in a)
b = ''.join(str(i) for i in b)

After converting to string, you can easily find the index of substring with the following string function.

a.index(b)

Cheers!!

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.