numpy中的广义通用函数

Question

我正在尝试使用 numpy API 制作一个通用的 ufunc。输入是一个(n x m) 矩阵和一个标量，输出是两个矩阵（(n x p) 和(p x m)）。但我不知道该怎么做。有人可以帮助我吗？ 在初始化函数中，我使用带有签名的PyUFunc_FromFuncAndDataAndSignature 函数：

"(n,m),()->(n,p),(p,m)"

我可以读取输入（矩阵和标量），但我想使用标量输入，例如签名中的维度 p。有可能吗？

这是一个仅打印输入的示例代码：

#include "Python.h"
#include "math.h"
#include "numpy/ndarraytypes.h"
#include "numpy/ufuncobject.h"

static PyMethodDef nmfMethods[] = {
        {NULL, NULL, 0, NULL}
};


static void double_nmf(char **args, npy_intp *dimensions,
                            npy_intp* steps, void* data)
{
    npy_intp i, j, 
             n = dimensions[1], //dimensions of input matrix
             m = dimensions[2]; //

    printf("scalar: %d\n",*(int*)args[1]); // input scalar

    // just print input matrix
    printf("Input matrix:\n");
    for(i=0;i<n;i++){
        for(j=0;j<m;j++){
            printf("%.1f ",*(double*)(args[0]+8*(i*m+j)));
        }
    printf("\n");
    }
    return;

}

static PyUFuncGenericFunction nmf_functions[] = { double_nmf };
static void * nmf_data[] = { (void *)NULL };
static char nmf_signatures[] = { PyArray_DOUBLE, PyArray_INT, PyArray_DOUBLE, PyArray_DOUBLE };
char *nmf_signature = "(n,m),()->(n,p),(p,m)";

PyMODINIT_FUNC initnmf(void)
{
    PyObject *m, *d, *version, *nmf;

    m = Py_InitModule("nmf", nmfMethods);
    if (m == NULL) {
        return;
    }

    import_array();
    import_umath();
    d = PyModule_GetDict(m);
    version = PyString_FromString("0.1");
    PyDict_SetItemString(d, "__version__", version);
    Py_DECREF(version);

    nmf = PyUFunc_FromFuncAndDataAndSignature(nmf_functions, nmf_data, nmf_signatures, 1,
                                    2, 2, PyUFunc_None, "nmf",
                                    "", 0, nmf_signature);
    PyDict_SetItemString(d, "nmf", nmf);
    Py_DECREF(nmf);
}

此代码可以编译并运行。 python脚本在这里：

#/usr/bin/python

import numpy as np
import nmf

x = np.array([[1,2,3,4,5],[6,7,8,9,10],[11,12,13,14,15],[16,17,18,19,20]])
y,z = nmf.nmf(x,2)
print "Shapes of outputs: ", y.shape, z.shape

终端输出为：

scalar: 2
Input matrix:
1.0 2.0 3.0 4.0 5.0 
6.0 7.0 8.0 9.0 10.0 
11.0 12.0 13.0 14.0 15.0 
16.0 17.0 18.0 19.0 20.0 
Shapes of outputs:  (4, 1) (1, 5)

我的疑问是如何使用标量输入（在这种情况下为 2），如输出矩阵的维度 p。在示例中 p = 1，我没有设置它。

Answer 1

除了提供一定大小的数组之外，没有办法在 gufunc 中设置维度。 1 是所有维度在内部初始化的值，你不应该依赖它不会改变。我个人对此的看法是，未定义的维度应该会引发错误。

设置p 的唯一方法是创建一个具有正确形状的空数组并将其作为输出数组传入。要实现它，您需要重新定义您的 nmf 以具有签名 "(m,n)->(m,p),(p,n)" 并用一些类似于以下的 Python 包装它：

def nmf_wrap(x, p):
    x = np.asarray(x)
    assert x.ndim >= 2
    shape = x.shape[:-2]
    m, n = x.shape[-2:]
    out1 = np.empty(shape+(m, p), dtype=x.dtype)
    out2 = np.empty(shape+(p, n), dtype=x.dtype)
    return nmf.nmf(x, out1, out2)

关于扩展 gufuncs 签名 on the numpy-dev mailing list recently 提供的功能的讨论正在进行中。你所描述的与我所说的“计算维度”有一些相似之处。如果您想看到在 numpy 1.10 中实现的某些东西，如果您能在该列表中更详细地解释您的用例，那就太好了：我们不知道有很多（任何？） gufunc 编码器在野外！

Answer 2

感谢@jaime 的回答，对我帮助很大。我进行了您建议的更改，并且有效。 这里是 C 示例代码。它只是将一些输入矩阵元素复制到输出。

#include "Python.h"
#include "numpy/ndarraytypes.h"
#include "numpy/ufuncobject.h"

static PyMethodDef nmfMethods[] = {
        {NULL, NULL, 0, NULL}
};


static void double_nmf(char **args, npy_intp *dimensions,
                            npy_intp* steps, void* data)
{
    npy_intp i, j,
            n = dimensions[1],
            m = dimensions[2],
            p = dimensions[3];
    char *in = args[0], *out1 = args[1], *out2 = args[2];

    for(i=0; i<n; i++){
        for(j=0; j<p; j++){
            *(double*)(out1 + 8*(j + p*i)) = *(double*)(in + 8*(j + m*i));
        }
    }
    for(i=0; i<p; i++){
        for(j=0; j<m; j++){
            *(double*)(out2 + 8*(j + m*i)) = *(double*)(in + 8*(j + m*i));
        }
    }
    return;

}

static PyUFuncGenericFunction nmf_functions[] = { double_nmf };
static void * nmf_data[] = { (void *)NULL };
static char nmf_signatures[] = { PyArray_DOUBLE, PyArray_DOUBLE, PyArray_DOUBLE };
char *nmf_signature = "(n,m)->(n,p),(p,m)";

PyMODINIT_FUNC initnmf(void)
{
    PyObject *m, *d, *version, *nmf;

    m = Py_InitModule("nmf", nmfMethods);
    if (m == NULL) {
        return;
    }

    import_array();
    import_umath();
    d = PyModule_GetDict(m);
    version = PyString_FromString("0.1");
    PyDict_SetItemString(d, "__version__", version);
    Py_DECREF(version);

    nmf = PyUFunc_FromFuncAndDataAndSignature(nmf_functions, nmf_data, nmf_signatures, 1,
                                    1, 2, PyUFunc_None, "nmf",
                                    "", 0, nmf_signature);
    PyDict_SetItemString(d, "nmf", nmf);
    Py_DECREF(nmf);
}

这里是 Python 示例代码和终端输出。

#/usr/bin/python

import numpy as np
import nmf

def nmf_wrap(x,p):
    x = np.asarray(x)
    assert x.ndim >=2
    shape = x.shape[-2:]
    n,m = shape[-2:]
    out1 = np.empty((n, p), dtype=x.dtype)
    out2 = np.empty((p, m), dtype=x.dtype)
    return nmf.nmf(x, out1, out2)

x = np.array([[1,2,3,4,5],[6,7,8,9,10],[11,12,13,14,15],[16,17,18,19,20]])
y,z = nmf_wrap(x,2)
print 'Input:\n', x
print 'Output 1:\n', y
print 'Output 2:\n', z

---

Input:
[[ 1  2  3  4  5]
 [ 6  7  8  9 10]
 [11 12 13 14 15]
 [16 17 18 19 20]]
Output 1:
[[ 1  2]
 [ 6  7]
 [11 12]
 [16 17]]
Output 2:
[[ 1  2  3  4  5]
 [ 6  7  8  9 10]]

现在我可以继续编程了。

我也在编写非 ufunc，使用 Python/C API（没有 numpy），正如@nathaniel-j-smith 建议的那样（这就是我得到的。）。我希望能有一些结果，并简要比较两种方法。