44884


ある配列 `a` のサブセットの配列を `b` としたときに、 `b` の各要素が `a` の何番目にあるか、を調べたいです。  
具体的には以、下のコードを `numpy` を使って高速に実現したいのですが、やり方がわかりません。  
何卒ご教示お願いいたします。

[code]

    a = [1,4,5,2,3,6]
    b = [1,4,3]
    
    I = [-1] * len(b)
    last = 0
    
    for i,n in enumerate(b):
        I[i] = a.index(n, last)
        last = I[i] + 1
    print(I)
    
[/code]

2018/06/20  
言葉足らずでしたので補足させていただきます。  
任意の `i` について `b[i] == a[j]` を満たす `j` が必ず存在します。  
また、a と b の関係は、順番も保存されます。  
つまり `b[i] == a[j] and b[i+1] == a[k]` のとき、常に `j<k` です。  
何卒よろしくお願いします。

2018/06/23  
言葉足らずでしたので補足させていただきます。  
bの各要素が、「aの中での順番」通りにならんでいることは保証されますが、a自体はソートされているわけではありません。  
合わせてサンプルコードを修正しました。  
度々申し訳ありません。


CC BY-SA 4.0
2018-06-19T07:38:34.417
44880
2018-06-24T03:01:12.277
2018-06-23T13:58:44.090
28998
28998
question
python アルゴリズム numpy
numpyで2つの配列をシーケンシャルに比較したい


NumPy だけでやる場合、配列の要素に対して処理する際に `last`
をグローバル変数として保持しておきながら高速に処理する方法を思い付けませんでした。`numpy.vectorize` や
`numpy.ufunc.reduce` を使ってできないかと思ったのですが、上手く行きませんでした。また、NumPy 1.14.0
には[「配列から、ある要素が含まれる最初のインデックスだけを出力する」という関数が未だに無い](https://stackoverflow.com/q/7632963/5989200)ので、`a`
がソートされていない場合は NumPy のみでは素早い実行ができなさそうです。

代わりに Cython を使って高速化することはできたので、参考までに IPython での実装を共有いたします。こちらの環境では生 Python
での実装に比べ 4 倍ほど速く実行されました。

[code]

    In [1]: import numpy as np
    
            # データの準備
            len_a = 10000
            len_b = 5000
            a_arr = np.random.randint(low=0, high=10000, size=len_a)
            b_mask = np.concatenate((np.ones((len_b,), dtype=bool), np.zeros((len_a - len_b,), dtype=bool)))
            np.random.shuffle(b_mask)
            b_arr = a_arr[b_mask]
    
            a = list(a_arr)
            b = list(b_arr)
    In [2]: def search_plain(a, b):
                result = [-1] * len(b)
                last = 0
                for i, n in enumerate(b):
                    result[i] = a.index(n, last)
                    last = result[i] + 1
                return result
    
            %timeit search_plain(a, b)
            1.93 ms ± 150 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
    In [3]: %reload_ext Cython
    In [4]: %%cython
            # cython: boundscheck=False, wraparound=False
    
            def search_cython(a, b):
                cdef int i, j
                cdef int last = 0
                cdef int len_a = len(a)
                cdef int len_b = len(b)
                result = [-1] * len_b
                for i in range(len_b):
                    for j in range(last, len_a):
                        if a[j] == b[i]:
                            result[i] = j
                            last = j + 1
                            break
                return result
    In [5]: %timeit search_cython(a, b)
            462 µs ± 14.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
    In [6]: # asserting...
    
            result_plain = search_plain(a, b)
            result_cython = search_cython(a, b)
            print("search_cython is valid:", set(result_plain) == set(result_cython))
            search_cython is valid: True
    
[/code]

* * *

## `b` の順序が保証されない場合

**注意: 以下は私がご質問の内容を誤解して回答した際の、古いものです。**

NumPy array
から特定の値を持つインデックスを得るには、[`numpy.where`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.where.html)
を使います。その値が複数ヶ所にある場合、`numpy.where` はそれら全てを返します。

[code]

    >>> import numpy as np
    >>> a = np.array([1, 1, 2, 3, 5, 8])
    >>> np.where(a == 3)
    (array([3]),)
    >>> np.where(a == 1)
    (array([0, 1]),)
    
[/code]

これを単純に `b` の各要素について実行すれば良いのであれば、次のように書けます。

[code]

    >>> a = np.array([1, 2, 3, 4, 5, 6])
    >>> b = np.array([1, 4, 6])
    >>> np.array([np.where(a == x) for x in b])
    array([[[0]],
    
           [[3]],
    
           [[5]]])
    
[/code]

`a` が一次元配列のときに限るのであれば、次のように書いた方が冗長性を省けます。

[code]

    >>> np.array([np.where(a == x)[0] for x in b])
    array([[0],
           [3],
           [5]])
    
[/code]

リスト内包表記では遅いということであれば、[`numpy.frompyfunc`](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.frompyfunc.html)
や
[`numpy.fromiter`](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.fromiter.html)
を使う方法もあります
([参考1](https://stackoverflow.com/q/14372613/5989200)、[2](https://qiita.com/ysk24ok/items/6736c8d8cc6eb06c6aa1))。

[code]

    >>> where = lambda x: np.where(a == x)[0]
    >>> np.frompyfunc(where, 1, 1)(b)
    array([array([0]), array([3]), array([5])], dtype=object)
    
[/code]


CC BY-SA 4.0
2018-06-19T08:33:47.207
44884
2018-06-22T07:07:38.023
2018-06-22T07:07:38.023
19110
19110
44880
answer


生Pythonではなく、Numpyを使うようにしてから、NumbaかCythonを使うようにすると30倍以上高速化します。途中までは、nekketsuuu♦さんの回答にある
IPython での実装とほぼ同じだったので省略して、追加部分のみ書いておきます。

Numbaの方は事前コンパイルが必要でないので手軽に使えます。

[code]

    import numba
    
    @numba.jit
    def search_numba(a, b):
        last = 0
        result = np.full(len_b, -1)
        for i in range(len_b):
            for j in range(last, len_a):
                if a[j] == b[i]:
                    result[i] = j
                    last = j + 1
                    break
        return result
    
    %timeit search_numba(a_arr, b_arr)
    10000 loops, best of 3: 57.5 µs per loop
    
    %%cython -a
    # cython: boundscheck=False, wraparound=False
    cimport numpy as np
    import numpy as np
    
    cpdef np.ndarray[np.int_t] search_cython2(np.ndarray[np.int_t] a, np.ndarray[np.int_t] b):
        cdef int i, j
        cdef int last = 0
        cdef int len_a = len(a)
        cdef int len_b = len(b)
        cdef np.ndarray[np.int_t] result = np.full(len_b, -1)
        for i in range(len_b):
            for j in range(last, len_a):
                if a[j] == b[i]:
                    result[i] = j
                    last = j + 1
                    break
        return result
    
    %timeit search_cython2(a_arr, b_arr)
    10000 loops, best of 3: 57.7 µs per loop
    
[/code]

NumbaもCythonもほぼ同じ処理速度になりましたが、Cythonの方の最適化が不十分なためだと思われます。

### 以下は、サンプルコードを修正する前の回答です。

この問題には、リストに値を挿入する際のインデックスを取得する`numpy.searchsorted`が使えます。aが既にソートされていること、aに同じ値が存在しない場合でも`np.nan`にはならないのですが、任意の
i について b[i] == a[j] を満たす j が必ず存在するということなので、単純に以下で計算できます。

[code]

    >>> import numpy as np
    >>> a = np.array([1, 2, 3, 4, 5, 6])
    >>> b = np.array([1, 4, 6])
    >>> numpy.searchsorted(a, b)
    array([0, 3, 5])
    
[/code]

重複がある場合は、左側と右側の位置を計算すれば求められます。また、値が存在しない場合は左側と右側の位置が同じになるのでそれもわかります。

[code]

    >>> a = np.array([1, 2, 3, 4, 4, 5, 6, 6, 6])
    >>> b = np.array([1, 4, 6])
    >>> numpy.searchsorted(a, b)
    array([0, 3, 6])
    >>> numpy.searchsorted(a, b, side='right')
    array([1, 5, 9])
    
[/code]


CC BY-SA 4.0
2018-06-22T05:07:25.170
44983
2018-06-24T03:01:12.277
2018-06-24T03:01:12.277
15171
15171
44880
answer