7

いくつかの型を追加してコンパイルするだけで、Python関数をcythonに変換しました。Python 関数と cython 関数の結果の間に小さな数値の違いがありました。いくつかの作業の後、違いは int の代わりに unsigned int を使用して numpy 配列にアクセスしたことによるものであることがわかりました。

http://docs.cython.org/src/userguide/numpy_tutorial.html#tuning-indexing-furtherに従って、アクセスを高速化するために unsigned int インデックスを使用していました。

とにかく、unsigned int を使用しても無害だと思いました。

このコードを参照してください:

cpdef function(np.ndarray[np.float32_t, ndim=2] response, max_loc):
    cdef unsigned int x, y   
    x, y = int(max_loc[0]), int(max_loc[1])
    x2, y2 = int(max_loc[0]), int(max_loc[1])
    print response[y,x], type(response[y,x]), response.dtype
    print response[y2,x2], type(response[y2,x2]), response.dtype   
    print 2*(response[y,x] - min(response[y,x-1], response[y,x+1]))
    print 2*(response[y2,x2] - min(response[y2,x2-1], response[y2,x2+1]))  

プリント:

0.959878861904 <type 'float'> float32
0.959879 <type 'numpy.float32'> float32
1.04306024313
1.04306030273   

なぜこれが起こるのですか?バグですか?

OK、ここで要求されているのは、元の関数で使用したのと同じ型と値を持つ SSCCE です

cpdef function():
    cdef unsigned int x, y  
    max_loc2 = np.asarray([ 15., 25.], dtype=float) 
    cdef np.ndarray[np.float32_t, ndim=2] response2 = np.zeros((49,49), dtype=np.float32)    
    x, y = int(max_loc2[0]), int(max_loc2[1])
    x2, y2 = int(max_loc2[0]), int(max_loc2[1])

    response2[y,x] = 0.959878861904  
    response2[y,x-1] = 0.438348740339
    response2[y,x+1] = 0.753262758255  


    print response2[y,x], type(response2[y,x]), response2.dtype
    print response2[y2,x2], type(response2[y2,x2]), response2.dtype
    print 2*(response2[y,x] - min(response2[y,x-1], response2[y,x+1]))
    print 2*(response2[y2,x2] - min(response2[y2,x2-1], response2[y2,x2+1]))  

版画

0.959878861904 <type 'float'> float32
0.959879 <type 'numpy.float32'> float32
1.04306024313
1.04306030273

私はpython 2.7.3 cython 0.18とmsvc9 expressを使用しています

4

2 に答える 2

7

モジュール用に生成された C ソースを読みやすくするために、質問の例を修正しました。配列からオブジェクトfloatを取得するのではなく、Python オブジェクトを作成するロジックにのみ興味があります。np.float32response

pyximport拡張モジュールのコンパイルに使用しています。生成された C ファイルを~/.pyxbld(おそらく%userprofile%\.pyxbldWindows の場合) のサブディレクトリに保存します。

import numpy as np
import pyximport
pyximport.install(setup_args={'include_dirs': [np.get_include()]})

open('_tmp.pyx', 'w').write('''
cimport numpy as np
cpdef function(np.ndarray[np.float32_t, ndim=2] response, max_loc):
    cdef unsigned int p_one, q_one
    p_one = int(max_loc[0])
    q_one = int(max_loc[1])
    p_two = int(max_loc[0])
    q_two = int(max_loc[1])
    r_one = response[q_one, p_one]
    r_two = response[q_two, p_two]
''')

import _tmp
assert(hasattr(_tmp, 'function'))

関心のあるセクションの生成された C コードを次に示します (読みやすくするために少し再フォーマットされています)。Cunsigned intインデックス変数を使用すると、生成されたコードが配列バッファーから直接データを取得し、 を呼び出してPyFloat_FromDouble、 に強制することがわかりdoubleます。一方、Pythonintインデックス変数を使用する場合は、一般的なアプローチが取られます。タプルを形成し、 を呼び出しますPyObject_GetItem。このようにして、 がdtypendarrayを正しく尊重できるようになります。np.float32

#define __Pyx_BufPtrStrided2d(type, buf, i0, s0, i1, s1) \
    (type)((char*)buf + i0 * s0 + i1 * s1)

  /* "_tmp.pyx":9
 *     p_two = int(max_loc[0])
 *     q_two = int(max_loc[1])
 *     r_one = response[q_one, p_one]             # <<<<<<<<<<<<<<
 *     r_two = response[q_two, p_two]
 */
  __pyx_t_3 = __pyx_v_q_one;
  __pyx_t_4 = __pyx_v_p_one;
  __pyx_t_5 = -1;

  if (unlikely(__pyx_t_3 >= (size_t)__pyx_bshape_0_response))
    __pyx_t_5 = 0;
  if (unlikely(__pyx_t_4 >= (size_t)__pyx_bshape_1_response))
    __pyx_t_5 = 1;

  if (unlikely(__pyx_t_5 != -1)) {
    __Pyx_RaiseBufferIndexError(__pyx_t_5);
    {
      __pyx_filename = __pyx_f[0];
      __pyx_lineno = 9;
      __pyx_clineno = __LINE__;
      goto __pyx_L1_error;
    }
  }

  __pyx_t_1 = PyFloat_FromDouble((
    *__Pyx_BufPtrStrided2d(
      __pyx_t_5numpy_float32_t *,
      __pyx_bstruct_response.buf,
      __pyx_t_3, __pyx_bstride_0_response,
      __pyx_t_4, __pyx_bstride_1_response)));

  if (unlikely(!__pyx_t_1)) {
    __pyx_filename = __pyx_f[0];
    __pyx_lineno = 9;
    __pyx_clineno = __LINE__;
    goto __pyx_L1_error;
  }

  __Pyx_GOTREF(__pyx_t_1);
  __pyx_v_r_one = __pyx_t_1;
  __pyx_t_1 = 0;

  /* "_tmp.pyx":10
 *     q_two = int(max_loc[1])
 *     r_one = response[q_one, p_one]
 *     r_two = response[q_two, p_two]             # <<<<<<<<<<<<<<
 */
  __pyx_t_1 = PyTuple_New(2);

  if (unlikely(!__pyx_t_1)) {
    __pyx_filename = __pyx_f[0];
    __pyx_lineno = 10;
    __pyx_clineno = __LINE__;
    goto __pyx_L1_error;
  }

  __Pyx_GOTREF(((PyObject *)__pyx_t_1));
  __Pyx_INCREF(__pyx_v_q_two);
  PyTuple_SET_ITEM(__pyx_t_1, 0, __pyx_v_q_two);
  __Pyx_GIVEREF(__pyx_v_q_two);
  __Pyx_INCREF(__pyx_v_p_two);
  PyTuple_SET_ITEM(__pyx_t_1, 1, __pyx_v_p_two);
  __Pyx_GIVEREF(__pyx_v_p_two);

  __pyx_t_2 = PyObject_GetItem(
    ((PyObject *)__pyx_v_response),
    ((PyObject *)__pyx_t_1));

  if (!__pyx_t_2) {
    __pyx_filename = __pyx_f[0];
    __pyx_lineno = 10;
    __pyx_clineno = __LINE__;
    goto __pyx_L1_error;
  }

  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(((PyObject *)__pyx_t_1));
  __pyx_t_1 = 0;
  __pyx_v_r_two = __pyx_t_2;
  __pyx_t_2 = 0;
于 2013-03-10T12:28:00.377 に答える
2

私のマシンでこれをいじっても、違いはわかりません。私は cython マジックで ipython ノートブックを使用しています:

In [1]:

%load_ext cythonmagic

In [12]:

%%cython

import numpy as np
cimport numpy as np

cpdef function(np.ndarray[np.float32_t, ndim=2] response, max_loc):
    cdef unsigned int x, y   
    x, y = int(max_loc[0]), int(max_loc[1])
    x2, y2 = int(max_loc[0]), int(max_loc[1])
    #return 2*(response[y,x] - min(response[y,x-1], response[y,x+1])), 2*(response[y2,x2] - min(response[y2,x2-1], response[y2,x2+1]))
    print response[y,x], type(response[y,x]), response.dtype
    print response[y2,x2], type(response[y2,x2]), response.dtype   
    print 2*(response[y,x] - min(response[y,x-1], response[y,x+1]))
    print 2*(response[y2,x2] - min(response[y2,x2-1], response[y2,x2+1])) 

In [13]:

a = np.random.normal(size=(10,10)).astype(np.float32)
m = [3,2]
function(a,m)

0.586090564728 <type 'float'> float32
0.586091 <type 'numpy.float32'> float32
4.39655685425
4.39655685425

結果の最初のペアです。違いは、print ステートメントの出力精度だけです。Cython のどのバージョンを使用していますか? numpy 配列のデータ属性が格納している固定長のメモリにアクセスしているだけなので、インデックスが答えに影響を与える可能性はほとんどありません。

于 2013-03-10T01:50:20.683 に答える