2.2.2.2 练习：从0开始构建一个ufunc

Mandelbrot分形由如下迭代定义：

2.2.2.2 练习：从0开始构建一个ufunc - 图1

C=X+iy是一个复数，只要Z仍然是有限的，无论迭代要跑多久，迭代都会重复。C属于Mandelbrot集。

  • ufunc调用mandel(z0, c)计算：

In [ ]:

  1. z = z0
  2. for k in range(iterations):
  3.     z = z*z + c

比如，一百次迭代或者直到z.real2 + z.imag2 > 1000。用它来决定哪个C是在Mandelbrot集中。

  • 我们的函数是很简单的，因此，请利用PyUFunc_*助手。
  • 用Cython来完成

也可以看一下mandel.pyx，mandelplot.py

提醒：一些预设Ufunc循环：

PyUfunc_f_ffloat elementwise_func(float input_1)
PyUfunc_ff_ffloat elementwise_func(float input_1, float input_2)
PyUfunc_d_ddouble elementwise_func(double input_1)
PyUfunc_dd_ddouble elementwise_func(double input_1, double input_2)
PyUfunc_D_Delementwise_func(complex_double input, complex_double output)
PyUfunc_DD_Delementwise_func(complex_double in1, complex_double in2, complex_double* out)

打印代码：

NPY_BOOL, NPY_BYTE, NPY_UBYTE, NPY_SHORT, NPY_USHORT, NPY_INT, NPY_UINT, NPY_LONG, NPY_ULONG, NPY_LONGLONG, NPY_ULONGLONG, NPY_FLOAT, NPY_DOUBLE, NPY_LONGDOUBLE, NPY_CFLOAT, NPY_CDOUBLE, NPY_CLONGDOUBLE, NPY_DATETIME, NPY_TIMEDELTA, NPY_OBJECT, NPY_STRING, NPY_UNICODE, NPY_VOID

2.2.2.3 答案：从0开始创建一个ufunc

In [ ]:

  1. # The elementwise function
  2. #  ------------------------
  4. cdef void mandel_single_point(double complex *z_in, 
  5.                               double complex *c_in,
  6.                               double complex *z_out) nogil:
  7.     #
  8.     # The Mandelbrot iteration
  9.     #
  11.     #
  12.     # Some points of note:
  13.     #
  14.     # - It's *NOT* allowed to call any Python functions here.
  15.     #
  16.     #   The Ufunc loop runs with the Python Global Interpreter Lock released.
  17.     #   Hence, the ``nogil``.
  18.     #
  19.     # - And so all local variables must be declared with ``cdef``
  20.     #
  21.     # - Note also that this function receives *pointers* to the data;
  22.     #   the "traditional" solution to passing complex variables around
  23.     #
  25.     cdef double complex z = z_in[0]
  26.     cdef double complex c = c_in[0]
  27.     cdef int k  # the integer we use in the for loop
  29.     # Straightforward iteration
  31.     for k in range(100):
  32.         z = z*z + c
  33.         if z.real**2 + z.imag**2 > 1000:
  34.             break
  36.     # Return the answer for this point
  37.     z_out[0] = z
  39. # Boilerplate Cython definitions
  40. #
  41. # You don't really need to read this part, it just pulls in
  42. # stuff from the Numpy C headers.
  43. #  ------------- ------------- ------------- -------------------
  45. cdef extern from "numpy/arrayobject.h":
  46.     void import_array()
  47.     ctypedef int npy_intp
  48.     cdef enum NPY_TYPES:
  49.         NPY_CDOUBLE
  51. cdef extern from "numpy/ufuncobject.h":
  52.     void import_ufunc()
  53.     ctypedef void (*PyUFuncGenericFunction)(char**, npy_intp*, npy_intp*, void*)
  54.     object PyUFunc_FromFuncAndData(PyUFuncGenericFunction* func, void** data,
  55.         char* types, int ntypes, int nin, int nout,
  56.         int identity, char* name, char* doc, int c)
  58.     void PyUFunc_DD_D(char**, npy_intp*, npy_intp*, void*)
  60. # Required module initialization
  61. #  ------------- -----------------
  63. import_array()
  64. import_ufunc()
  66. # The actual ufunc declaration
  67. #  ------------- ---------------
  69. cdef PyUFuncGenericFunction loop_func[1]
  70. cdef char input_output_types[3]
  71. cdef void *elementwise_funcs[1]
  73. loop_func[0] = PyUFunc_DD_D
  75. input_output_types[0] = NPY_CDOUBLE
  76. input_output_types[1] = NPY_CDOUBLE
  77. input_output_types[2] = NPY_CDOUBLE
  79. elementwise_funcs[0] = <void*>mandel_single_point
  81. mandel = PyUFunc_FromFuncAndData(
  82.     loop_func,
  83.     elementwise_funcs,
  84.     input_output_types,
  85.     1, # number of supported input types
  86.     2, # number of input args
  87.     1, # number of output args
  88.     0, # `identity` element, never mind this
  89.     "mandel", # function name
  90.     "mandel(z, c) -> computes iterated z*z + c", # docstring
  91.     0 # unused
  92.     )

In [ ]:

  1. import numpy as np
  2. import mandel
  3. x = np.linspace(-1.7, 0.6, 1000)
  4. y = np.linspace(-1.4, 1.4, 1000)
  5. c = x[None,:] + 1j*y[:,None]
  6. z = mandel.mandel(c, c)
  8. import matplotlib.pyplot as plt
  9. plt.imshow(abs(z)**2 < 1000, extent=[-1.7, 0.6, -1.4, 1.4])
  10. plt.gray()
  11. plt.show()

2.2.2.2 练习：从0开始构建一个ufunc - 图2

笔记：大多数模板可以由下列Cython模块来自动完成： http://wiki.cython.org/MarkLodato/CreatingUfuncs

一些可接受的输入类型

例如：支持小数点后一位及两位两个准确度版本

In [ ]:

  1. cdef void mandel_single_point(double complex *z_in,
  2.                               double complex *c_in,
  3.                               double complex *z_out) nogil:
  4.    ...
  6. cdef void mandel_single_point_singleprec(float complex *z_in,
  7.                                          float complex *c_in,
  8.                                          float complex *z_out) nogil:
  9.    ...
  11. cdef PyUFuncGenericFunction loop_funcs[2]
  12. cdef char input_output_types[3*2]
  13. cdef void *elementwise_funcs[1*2]
  15. loop_funcs[0] = PyUFunc_DD_D
  16. input_output_types[0] = NPY_CDOUBLE
  17. input_output_types[1] = NPY_CDOUBLE
  18. input_output_types[2] = NPY_CDOUBLE
  19. elementwise_funcs[0] = <void*>mandel_single_point
  21. loop_funcs[1] = PyUFunc_FF_F
  22. input_output_types[3] = NPY_CFLOAT
  23. input_output_types[4] = NPY_CFLOAT
  24. input_output_types[5] = NPY_CFLOAT
  25. elementwise_funcs[1] = <void*>mandel_single_point_singleprec
  27. mandel = PyUFunc_FromFuncAndData(
  28.     loop_func,
  29.     elementwise_funcs,
  30.     input_output_types,
  31.     2, # number of supported input types   < ----------------
  32.     2, # number of input args
  33.     1, # number of output args
  34.     0, # `identity` element, never mind this
  35.     "mandel", # function name
  36.     "mandel(z, c) -> computes iterated z*z + c", # docstring
  37.     0 # unused
  38.     )