<p>根据一些评论,使用cython似乎是最好的方法。我愚蠢地从来没有考虑过使用赛顿。事实证明,生成工作代码相对容易。在</p>
<p>经过一番搜索,我把下面的cython代码放在一起。这不是最通用的代码,可能不是最好的编写方法,而且可能会变得更高效。即使如此,它只比原始问题中的<code>einsum()</code>代码慢25%,所以它不是太糟糕!它的目的是显式地处理在原始问题中创建的数组(因此是输入数组的假定模式)。<br/>
尽管有警告,但它确实为原始问题提供了一个合理有效的解决方案,并且可以作为类似情况下的起点。在</p>
<pre><code>import numpy as np
cimport numpy as np
import cython
DTYPE = np.float64
ctypedef np.float64_t DTYPE_t
cdef inline double d_abs (double a) : return a if a >= 0 else -a
@cython.boundscheck(False)
@cython.wraparound(False)
def process_vectors (np.ndarray[DTYPE_t, ndim=2, mode="fortran"] nhat not None,
np.ndarray[DTYPE_t, ndim=2, mode="c"] m not None) :
if nhat.shape[1] != m.shape[1] :
raise ValueError ("Arrays must contain vectors of the same dimension")
cdef Py_ssize_t imax = nhat.shape[0]
cdef Py_ssize_t jmax = m.shape[0]
cdef Py_ssize_t kmax = nhat.shape[1] # same as m.shape[1]
cdef np.ndarray[DTYPE_t, ndim=1] S = np.zeros(imax, dtype=DTYPE)
cdef Py_ssize_t i, j, k
cdef DTYPE_t val, tmp
for i in range(imax) :
val = 0
for j in range(jmax) :
tmp = 0
for k in range(kmax) :
tmp += nhat[i,k] * m[j,k]
val += d_abs(tmp)
S[i] = val / jmax
return S
</code></pre>