NumPy 中的 `ix_()` 函数到底是什么意思？怎么用？

[daa233]

The ix_() function in Quickstart tutorial of NumPy.

numpy.ix in NumPy's documentation.

[daa233]

"""
ix_(*args)
    Construct an open mesh from multiple sequences.

    This function takes N 1-D sequences and returns N outputs with N
    dimensions each, such that the shape is 1 in all but one dimension
    and the dimension with the non-unit shape value cycles through all
    N dimensions.

    Using `ix_` one can quickly construct index arrays that will index
    the cross product. ``a[np.ix_([1,3],[2,5])]`` returns the array
    ``[[a[1,2] a[1,5]], [a[3,2] a[3,5]]]``.

    Parameters
    ----------
    args : 1-D sequences
        Each sequence should be of integer or boolean type.
        Boolean sequences will be interpreted as boolean masks for the
        corresponding dimension (equivalent to passing in
        ``np.nonzero(boolean_sequence)``).

    Returns
    -------
    out : tuple of ndarrays
        N arrays with N dimensions each, with N the number of input
        sequences. Together these arrays form an open mesh.

    See Also
    --------
    ogrid, mgrid, meshgrid

    Examples
    --------
    >>> a = np.arange(10).reshape(2, 5)
    >>> a
    array([[0, 1, 2, 3, 4],
           [5, 6, 7, 8, 9]])
    >>> ixgrid = np.ix_([0, 1], [2, 4])
    >>> ixgrid
    (array([[0],
           [1]]), array([[2, 4]]))
    >>> ixgrid[0].shape, ixgrid[1].shape
    ((2, 1), (1, 2))
    >>> a[ixgrid]
    array([[2, 4],
           [7, 9]])

    >>> ixgrid = np.ix_([True, True], [2, 4])
    >>> a[ixgrid]
    array([[2, 4],
           [7, 9]])
    >>> ixgrid = np.ix_([True, True], [False, False, True, False, True])
    >>> a[ixgrid]
    array([[2, 4],
           [7, 9]])
"""

[daa233]

In [20]: a = np.arange(25).reshape(5, 5)

In [21]: a
Out[21]:
array([[ 0,  1,  2,  3,  4],
       [ 5,  6,  7,  8,  9],
       [10, 11, 12, 13, 14],
       [15, 16, 17, 18, 19],
       [20, 21, 22, 23, 24]])

In [23]: a[np.ix_([1,3],[2,4])]
Out[23]:
array([[ 7,  9],
       [17, 19]])

In [24]: np.ix_([1,3],[2,4])    # cross product indices
Out[24]:
(array([[1],
        [3]]), array([[2, 4]]))

[daa233]

引入 True 和 False 的情况，会根据值转换为具体的索引值：

In [25]: ixgrid = np.ix_([True, True], [2, 4])

In [26]: ixgrid
Out[26]:
(array([[0],
        [1]]), array([[2, 4]]))

In [27]: np.ix_([True, True, True], [2, 3, 4])
Out[27]:
(array([[0],
        [1],
        [2]]), array([[2, 3, 4]]))

In [28]: np.ix_([True, True], [False, False, True, False, True])
Out[28]:
(array([[0],
        [1]]), array([[2, 4]]))

In [29]: np.ix_([True, True], [False, False, True, True, True])
Out[29]:
(array([[0],
        [1]]), array([[2, 3, 4]]))

[daa233]

下面是 NumPy 的 Quickstart tutorial 中的例子，当时觉得有些复杂，没有看懂：

>>> a = np.array([2,3,4,5])
>>> b = np.array([8,5,4])
>>> c = np.array([5,4,6,8,3])
>>> ax,bx,cx = np.ix_(a,b,c)
>>> ax
array([[[2]],
       [[3]],
       [[4]],
       [[5]]])
>>> bx
array([[[8],
        [5],
        [4]]])
>>> cx
array([[[5, 4, 6, 8, 3]]])
>>> ax.shape, bx.shape, cx.shape
((4, 1, 1), (1, 3, 1), (1, 1, 5))
>>> result = ax+bx*cx
>>> result
array([[[42, 34, 50, 66, 26],
        [27, 22, 32, 42, 17],
        [22, 18, 26, 34, 14]],

       [[43, 35, 51, 67, 27],
        [28, 23, 33, 43, 18],
        [23, 19, 27, 35, 15]],

       [[44, 36, 52, 68, 28],
        [29, 24, 34, 44, 19],
        [24, 20, 28, 36, 16]],

       [[45, 37, 53, 69, 29],
        [30, 25, 35, 45, 20],
        [25, 21, 29, 37, 17]]])
>>> bx*cx  # add by me
array([[[40, 32, 48, 64, 24],
        [25, 20, 30, 40, 15],
        [20, 16, 24, 32, 12]]])
>>> result[3,2,4]
17
>>> a[3]+b[2]*c[4]
17

有了上面的基础，现在再看，可以看出 np.ix_ 可以推广到更高的维度，并且可以快速地对向量中的元素做同样的运算：

>>> ax,bx,cx = np.ix_(a,b,c)
>>> result = f(ax, bx, cx)
>>> assert result(x, y, z) == f(a[x], b[y], c[z])

[daa233]

np.ufunc.reduce 可以根据指定的运算方式减小向量维度：

In [56]: np.multiply.reduce([2, 3, 5])
Out[56]: 30

In [58]: X = np.arange(8).reshape((2, 2, 2))

In [59]: X
Out[59]:
array([[[0, 1],
        [2, 3]],

       [[4, 5],
        [6, 7]]])

In [60]: np.add.reduce(X, 0)
Out[60]:
array([[ 4,  6],
       [ 8, 10]])

In [61]: np.add.reduce(X, 1)
Out[61]:
array([[ 2,  4],
       [10, 12]])

In [62]: np.add.reduce(X, 2)
Out[62]:
array([[ 1,  5],
       [ 9, 13]])

也可以通过下面的方法来实现 np.ufunc.reduce 的功能：

In [47]: a = np.array([2,3,4,5])
In [48]: b = np.array([8,5,4])
In [49]: c = np.array([5,4,6,8,3])
In [50]: def ufunc_reduce(ufct, *vectors):
    ...:     vs = np.ix_(*vectors)
    ...:     r = ufct.identity
    ...:     for v in vs:
    ...:         r = ufct(r, v)
    ...:     return r
In [51]: result = ufunc_reduce(np.add, a, b, c)
Out[51]: result
array([[[15, 14, 16, 18, 13],
        [12, 11, 13, 15, 10],
        [11, 10, 12, 14,  9]],

       [[16, 15, 17, 19, 14],
        [13, 12, 14, 16, 11],
        [12, 11, 13, 15, 10]],

       [[17, 16, 18, 20, 15],
        [14, 13, 15, 17, 12],
        [13, 12, 14, 16, 11]],

       [[18, 17, 19, 21, 16],
        [15, 14, 16, 18, 13],
        [14, 13, 15, 17, 12]]])
In [52]: result[0, 0, 0] == np.add.reduce([a[0], b[0], c[0]])  # a[0] + b[0] + c[0] = 15
Out[52]: True
In [53]: result[0, 2, 4] == np.add.reduce([a[0], b[2], c[4]])  # a[0] + b[2] + c[4] = 9
Out[53]: True

这样实现的好处是利用了广播规则，就不用为数组大小创建额外参数了。我觉得 np.ix_ 的实质就是创建交叉项的操作元，之后便可以利用广播规则，能够快速构建交叉项索引。

关于 np.ufunc.identity，我不知道该怎么翻译，但有下面的用法：

>>> np.add.identity
0
>>> np.multiply.identity
1
>>> np.power.identity
1
>>> print(np.exp.identity)
None
>>> print(np.subtract.identity)
None
>>> print(np.divide.identity)
None

[knightstream]

谢谢，对官方文档有补充，更容易理解！

[daa233]

@knightstream 谢谢鼓励！哈哈，其实我现在回过头看还是要捋半天 😂