{ "cells": [ { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "%matplotlib inline\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", ">作者: Pauli Virtanen \n", ">原文: http://www.scipy-lectures.org/advanced/advanced_numpy/index.html \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Numpy是Python科学工具栈的基础。它的目的很简单：在一个内存块上实现针对多个条目（items）的高效操作。了解它的工作细节有助于有效的使用它的灵活性，使用有用的快捷方式，基于它构建新的工作。\n", "\n", "这个指南的目的包括：\n", "- 剖析Numpy数组，以及它的重要性。提示与技巧。\n", "- 通用函数：是什么、为什么以及如果你需要一个全新的该做什么。\n", "- 与其他工具整合：Numpy提供了一些方式将任意数据封装为ndarray，而不需要不必要的复制。\n", "- 新近增加的功能，对我来说他们包含什么：PEP 3118 buffers、广义ufuncs, ...\n", "\n", "\n", "---\n", "**先决条件**\n", "- Numpy (>= 1.2; 越新越好...)\n", "- Cython (>= 0.12, 对于Ufunc例子)\n", "- Pillow (Python 处理图形的库，在一些例子中使用)\n", "\n", "---\n", "\n", "在这个部分，numpy将被如下引入：" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "章节内容\n", "- ndarry的一生\n", " - 它是...\n", " - 内存块\n", " - 数据类型\n", " - 索引体系：strides\n", " - 剖析中的发现\n", "- 通用函数\n", " - 他们是什么？\n", " - 练习：从零开始构建一个ufunc\n", " - 答案：从零开始构建一个ufunc\n", " - 广义ufuncs\n", "- 协同工作功能\n", " - 共享多维度，类型数据\n", " - 旧的buffer协议\n", " - 旧的buffer协议\n", " - 数组接口协议\n", "- 数组切片：`chararray`、`maskedarray`、`matrix`\n", " - `chararray`：向量化字符操作\n", " - `masked_array` 缺失值\n", " - recarray：纯便利\n", " - `matrix`：便利？\n", "- 总结\n", "- 为Numpy/Scipy做贡献\n", " - 为什么\n", " - 报告bugs\n", " - 贡献文档\n", " - 贡献功能\n", " - 如何帮忙，总的来说" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2.2.1 ndarray的一生\n", "\n", "### 2.2.1.1 它是...\n", "\n", "**ndarray** =\n", "\n", " 内存块 + 索引体系 + 数据类型描述符\n", " \n", "- 原始数据\n", "- 如何定位(locate)一个元素\n", "- 如何解释(interpret)一个元素\n", "![ndarray](http://www.scipy-lectures.org/_images/threefundamental.png)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "```C\n", "typedef struct PyArrayObject {\n", " PyObject_HEAD\n", "\n", " /* 内存块 */\n", " char *data;\n", "\n", " /* 数据类型描述符 */\n", " PyArray_Descr *descr;\n", "\n", " /* 索引体系 */\n", " int nd;\n", " npy_intp *dimensions;\n", " npy_intp *strides;\n", "\n", " /* 其他属性 */\n", " PyObject *base;\n", " int flags;\n", " PyObject *weakreflist;\n", "} PyArrayObject;\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.2.1.2 内存块" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = np.array([1, 2, 3, 4], dtype=np.int32)\n", "x.data" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "''" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "str(x.data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "数据的内存地址：" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "140264211656336" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x.__array_interface__['data'][0] " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "完整的`__array_interface__`：" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'data': (140264211656336, False),\n", " 'descr': [('', '" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = b'1234' # 'b' 指代 \"bytes\", Python 3中要这么使用\n", "y = np.frombuffer(x, dtype=np.int8)\n", "y.data " ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y.base is x" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ " C_CONTIGUOUS : True\n", " F_CONTIGUOUS : True\n", " OWNDATA : False\n", " WRITEABLE : False\n", " ALIGNED : True\n", " UPDATEIFCOPY : False" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y.flags" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`owndata`和`writeable`标记表明了内存块的状态。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.2.1.3 数据类型\n", "#### 2.2.1.3.1 描述符\n", "`dtype`描述了数组里的一个项目：\n", "\n", "| | |\n", "|---|---|\n", "|type|数据的**标量类型**，int8、int16、float64等之一（固定大小），str、unicode、void（可变大小）|\n", "|itemsize|数据块的**大小**|\n", "|byteorder|**字节序**: big-endian `>` / little-endian `<` / 不可用 |\n", "|fields|子-dtypes，如果是一个**结构化的数据类型**|\n", "|shape|数组的形状，如果是一个**子数组**|" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "numpy.int64" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.dtype(int).type" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "8" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.dtype(int).itemsize" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'='" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.dtype(int).byteorder" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "#### 2.2.1.3.2 例子：读取.wav文件\n", "\n", "The`.wav` file header:\n", "\n", "| | |\n", "|---|---|\n", "|chunk_id|\"RIFF\"|\n", "|chunk_size|4字节无符号little-endian整型|\n", "|format|\"WAVE\"|\n", "|fmt_id|\"fmt \"|\n", "|fmt_size|4字节无符号little-endian整型|\n", "|audio_fmt|2字节无符号little-endian整型|\n", "|num_channels|2字节无符号little-endian整型|\n", "|sample_rate|4字节无符号little-endian整型|\n", "|byte_rate|4字节无符号little-endian整型|\n", "|block_align|2字节无符号little-endian整型|\n", "|bits_per_sample|2字节无符号little-endian整型|\n", "|data_id|\"data\"|\n", "|data_size|4字节无符号little-endian整型|\n", "\n", "- 44字节块的原始数据（在文件的开头）\n", "- ...接下来是`data_size` 实际声音数据的字节。\n", "\n", "`.wav`文件头是Numpy结构化数据类型：" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": true }, "outputs": [], "source": [ "wav_header_dtype = np.dtype([\n", " (\"chunk_id\", (str, 4)), # flexible-sized scalar type, item size 4\n", " (\"chunk_size\", \"'" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = np.array([[1, 2, 3],\n", " [4, 5, 6],\n", " [7, 8, 9]], dtype=np.int8)\n", "str(x.data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "item x[1,2]开始在`x.data`中的哪个字节？" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**答案**（在Numpy）\n", "- **步幅**：寻找一下个元素跳跃的字节数\n", "- 每个维度一个步幅" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(3, 1)" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x.strides" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "6" ] }, "execution_count": 44, "metadata": {}, "output_type": "execute_result" } ], "source": [ "byte_offset = 3*1 + 1*2 # 查找x[1,2]\n", "x.flat[byte_offset]" ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "6" ] }, "execution_count": 45, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x[1, 2]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- 简单、**灵活**\n", "\n", "##### 2.2.1.4.1.1 C 和 Fortran 顺序" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(6, 2)" ] }, "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = np.array([[1, 2, 3],\n", " [4, 5, 6],\n", " [7, 8, 9]], dtype=np.int16, order='C')\n", "x.strides" ] }, { "cell_type": "code", "execution_count": 48, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "str(x.data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- 需要跳跃6个字节寻找下一行\n", "- 需要跳跃2个字节寻找下一列" ] }, { "cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(2, 6)" ] }, "execution_count": 49, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y = np.array(x, order='F')\n", "y.strides" ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "''" ] }, "execution_count": 50, "metadata": {}, "output_type": "execute_result" } ], "source": [ "str(y.data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- 需要跳跃2个字节寻找下一行\n", "- 需要跳跃6个字节寻找下一列\n", "\n", "更高维度也类似：\n", " - C：最后的维度变化最快（=最小的步幅）\n", " - F：最早的维度变化最快\n", " \n", "![png](http://www.scipy-lectures.org/_images/math/f9d27734438964d22cb2d6449ef13c33c2b76f82.png)\n", "\n", "---\n", "\n", "**注意**：现在我们可以理解`.view()`的行为：\n" ] }, { "cell_type": "code", "execution_count": 51, "metadata": { "collapsed": true }, "outputs": [], "source": [ "y = np.array([[1, 3], [2, 4]], dtype=np.uint8).transpose()\n", "x = y.copy()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "变换顺序不影响数据的内部布局，只是步幅" ] }, { "cell_type": "code", "execution_count": 52, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(2, 1)" ] }, "execution_count": 52, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x.strides" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(1, 2)" ] }, "execution_count": 53, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y.strides" ] }, { "cell_type": "code", "execution_count": 54, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "''" ] }, "execution_count": 54, "metadata": {}, "output_type": "execute_result" } ], "source": [ "str(x.data)" ] }, { "cell_type": "code", "execution_count": 55, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "''" ] }, "execution_count": 55, "metadata": {}, "output_type": "execute_result" } ], "source": [ "str(y.data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- 当解释为int16时结果会不同\n", "- `.copy()`以C顺序（默认）创建新的数组\n", "\n", "---\n", "\n", "##### 2.2.1.4.1.2 用整数切片\n", "\n", "- 通过仅改变形状、步幅和可能调整数据指针可以代表任何东西！\n", "- 不用制造数据的副本" ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([6, 5, 4, 3, 2, 1], dtype=int32)" ] }, "execution_count": 56, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = np.array([1, 2, 3, 4, 5, 6], dtype=np.int32)\n", "y = x[::-1]\n", "y" ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(-4,)" ] }, "execution_count": 57, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y.strides" ] }, { "cell_type": "code", "execution_count": 58, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "8" ] }, "execution_count": 58, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y = x[2:]\n", "y.__array_interface__['data'][0] - x.__array_interface__['data'][0]" ] }, { "cell_type": "code", "execution_count": 59, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(800, 80, 8)" ] }, "execution_count": 59, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = np.zeros((10, 10, 10), dtype=np.float)\n", "x.strides" ] }, { "cell_type": "code", "execution_count": 60, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(1600, 240, 32)" ] }, "execution_count": 60, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x[::2,::3,::4].strides" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- 类似的，变换顺序绝不会创建副本（只是交换的步幅）" ] }, { "cell_type": "code", "execution_count": 61, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(800, 80, 8)" ] }, "execution_count": 61, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = np.zeros((10, 10, 10), dtype=np.float)\n", "x.strides" ] }, { "cell_type": "code", "execution_count": 62, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(8, 80, 800)" ] }, "execution_count": 62, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x.T.strides" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "但是：并不是所有的重排操作都可以通过操纵步幅来完成。" ] }, { "cell_type": "code", "execution_count": 63, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(1, 2)" ] }, "execution_count": 63, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a = np.arange(6, dtype=np.int8).reshape(3, 2)\n", "b = a.T\n", "b.strides" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "到目前为止，都很不错，但是：" ] }, { "cell_type": "code", "execution_count": 64, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "''" ] }, "execution_count": 64, "metadata": {}, "output_type": "execute_result" } ], "source": [ "str(a.data)" ] }, { "cell_type": "code", "execution_count": 65, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0, 2, 4],\n", " [1, 3, 5]], dtype=int8)" ] }, "execution_count": 65, "metadata": {}, "output_type": "execute_result" } ], "source": [ "b" ] }, { "cell_type": "code", "execution_count": 66, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0, 2, 4, 1, 3, 5], dtype=int8)" ] }, "execution_count": 66, "metadata": {}, "output_type": "execute_result" } ], "source": [ "c = b.reshape(3*2)\n", "c" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "这里，没办法用一个给定的步长和`a`的内存块来表示数组`c`。因此，重排操作在这里需要制作一个副本。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 2.2.1.4.2 例子：用步长伪造维度\n", "\n", "**步长操作**" ] }, { "cell_type": "code", "execution_count": 67, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on function as_strided in module numpy.lib.stride_tricks:\n", "\n", "as_strided(x, shape=None, strides=None, subok=False, writeable=True)\n", " Create a view into the array with the given shape and strides.\n", " \n", " .. warning:: This function has to be used with extreme care, see notes.\n", " \n", " Parameters\n", " ----------\n", " x : ndarray\n", " Array to create a new.\n", " shape : sequence of int, optional\n", " The shape of the new array. Defaults to ``x.shape``.\n", " strides : sequence of int, optional\n", " The strides of the new array. Defaults to ``x.strides``.\n", " subok : bool, optional\n", " .. versionadded:: 1.10\n", " \n", " If True, subclasses are preserved.\n", " writeable : bool, optional\n", " .. versionadded:: 1.12\n", " \n", " If set to False, the returned array will always be readonly.\n", " Otherwise it will be writable if the original array was. It\n", " is advisable to set this to False if possible (see Notes).\n", " \n", " Returns\n", " -------\n", " view : ndarray\n", " \n", " See also\n", " --------\n", " broadcast_to: broadcast an array to a given shape.\n", " reshape : reshape an array.\n", " \n", " Notes\n", " -----\n", " ``as_strided`` creates a view into the array given the exact strides\n", " and shape. This means it manipulates the internal data structure of\n", " ndarray and, if done incorrectly, the array elements can point to\n", " invalid memory and can corrupt results or crash your program.\n", " It is advisable to always use the original ``x.strides`` when\n", " calculating new strides to avoid reliance on a contiguous memory\n", " layout.\n", " \n", " Furthermore, arrays created with this function often contain self\n", " overlapping memory, so that two elements are identical.\n", " Vectorized write operations on such arrays will typically be\n", " unpredictable. They may even give different results for small, large,\n", " or transposed arrays.\n", " Since writing to these arrays has to be tested and done with great\n", " care, you may want to use ``writeable=False`` to avoid accidental write\n", " operations.\n", " \n", " For these reasons it is advisable to avoid ``as_strided`` when\n", " possible.\n", "\n" ] } ], "source": [ "from numpy.lib.stride_tricks import as_strided\n", "help(as_strided) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "**警告**：`as_strided`并不检查你是否还待在内存块边界里..\n", "\n", "---" ] }, { "cell_type": "code", "execution_count": 68, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([1, 3], dtype=int16)" ] }, "execution_count": 68, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = np.array([1, 2, 3, 4], dtype=np.int16)\n", "as_strided(x, strides=(2*2, ), shape=(2, ))" ] }, { "cell_type": "code", "execution_count": 69, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([1, 3], dtype=int16)" ] }, "execution_count": 69, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x[::2]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**练习**" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "array([1, 2, 3, 4], dtype=np.int8)\n", "\n", "-> array([[1, 2, 3, 4],\n", " [1, 2, 3, 4],\n", " [1, 2, 3, 4]], dtype=np.int8)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "仅使用`as_strided`.:\n", "\n", "提示：`byte_offset = stride[0]\\*index[0] + stride[1]\\*index[1] + ...`\n", "\n", "解密：\n", "\n", "步长可以设置为0：" ] }, { "cell_type": "code", "execution_count": 70, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[1, 2, 3, 4],\n", " [1, 2, 3, 4],\n", " [1, 2, 3, 4]], dtype=int8)" ] }, "execution_count": 70, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = np.array([1, 2, 3, 4], dtype=np.int8)\n", "y = as_strided(x, strides=(0, 1), shape=(3, 4))\n", "y" ] }, { "cell_type": "code", "execution_count": 71, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 71, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y.base.base is x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 2.2.1.4.3 广播\n", "- 用它来做一些有用的事情：[1, 2, 3, 4]和[5, 6, 7]的外积" ] }, { "cell_type": "code", "execution_count": 72, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[1, 2, 3, 4],\n", " [1, 2, 3, 4],\n", " [1, 2, 3, 4]], dtype=int16)" ] }, "execution_count": 72, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = np.array([1, 2, 3, 4], dtype=np.int16)\n", "x2 = as_strided(x, strides=(0, 1*2), shape=(3, 4))\n", "x2" ] }, { "cell_type": "code", "execution_count": 74, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[5, 5, 5, 5],\n", " [6, 6, 6, 6],\n", " [7, 7, 7, 7]], dtype=int16)" ] }, "execution_count": 74, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y = np.array([5, 6, 7], dtype=np.int16)\n", "y2 = as_strided(y, strides=(1*2, 0), shape=(3, 4))\n", "y2" ] }, { "cell_type": "code", "execution_count": 75, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 5, 10, 15, 20],\n", " [ 6, 12, 18, 24],\n", " [ 7, 14, 21, 28]], dtype=int16)" ] }, "execution_count": 75, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x2 * y2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**...看起来有一些熟悉...**" ] }, { "cell_type": "code", "execution_count": 76, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 5, 10, 15, 20],\n", " [ 6, 12, 18, 24],\n", " [ 7, 14, 21, 28]], dtype=int16)" ] }, "execution_count": 76, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = np.array([1, 2, 3, 4], dtype=np.int16)\n", "y = np.array([5, 6, 7], dtype=np.int16)\n", "x[np.newaxis,:] * y[:,np.newaxis]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* 在内部，数组**广播**的确使用0步长来实现的。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 2.2.1.4.4 更多技巧：对角线\n", "\n", "**挑战**\n", "\n", "- 提取矩阵对角线的起点：（假定是C内存顺序）：" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "```python\n", "x = np.array([[1, 2, 3],\n", " [4, 5, 6],\n", " [7, 8, 9]], dtype=np.int32)\n", "\n", "x_diag = as_strided(x, shape=(3,), strides=(???,)) \n", "```" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "- 提取第一个超级-对角线的起点[2,6]。\n", "- 那么子对角线呢？\n", "\n", "（后两个问题的提示：切片首先移动步长起点的点。）\n", "\n", "答案\n", "\n", "提取对角线：" ] }, { "cell_type": "code", "execution_count": 77, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 1, -29117, 2], dtype=int16)" ] }, "execution_count": 77, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x_diag = as_strided(x, shape=(3, ), strides=((3+1)*x.itemsize, ))\n", "x_diag" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "首先切片，调整数据指针：" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "as_strided(x[0, 1:], shape=(2, ), strides=((3+1)*x.itemsize, ))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "as_strided(x[1:, 0], shape=(2, ), strides=((3+1)*x.itemsize, ))" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "**备注： 使用 np.diag**" ] }, { "cell_type": "code", "execution_count": 80, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0, 1, 0, 0, 0],\n", " [0, 0, 2, 0, 0],\n", " [0, 0, 0, 3, 0],\n", " [0, 0, 0, 0, 4],\n", " [0, 0, 0, 0, 0]], dtype=int16)" ] }, "execution_count": 80, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y = np.diag(x, k=1)\n", "y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "但是" ] }, { "cell_type": "code", "execution_count": 81, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 81, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y.flags.owndata" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*** 注意 ***：Numpy 1.9之前，`np.diag` 会做一个拷贝。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**挑战**\n", "\n", "计算张量的迹：" ] }, { "cell_type": "code", "execution_count": 83, "metadata": { "collapsed": true }, "outputs": [], "source": [ "x = np.arange(5*5*5*5).reshape(5,5,5,5)\n", "s = 0\n", "for i in range(5):\n", " for j in range(5):\n", " s += x[j,i,j,i]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "通过跨越并且在结果上使用`sum()`。" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "```python\n", "y = as_strided(x, shape=(5, 5), strides=(TODO, TODO)) \n", "s2 = ... \n", "assert s == s2\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "答案" ] }, { "cell_type": "code", "execution_count": 84, "metadata": { "collapsed": true }, "outputs": [], "source": [ "y = as_strided(x, shape=(5, 5), strides=((5*5*5 + 5)*x.itemsize,\n", " (5*5 + 1)*x.itemsize))\n", "s2 = y.sum()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 2.2.1.4.5 CPU缓存效果\n", "\n", "内存布局可以影响性能：" ] }, { "cell_type": "code", "execution_count": 85, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "((20000,), (20000,))" ] }, "execution_count": 85, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = np.zeros((20000,))\n", "y = np.zeros((20000*67,))[::67]\n", "x.shape, y.shape" ] }, { "cell_type": "code", "execution_count": 88, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "10 µs ± 20.3 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "%timeit x.sum()" ] }, { "cell_type": "code", "execution_count": 89, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "37.2 µs ± 1.08 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n" ] } ], "source": [ "%timeit y.sum()" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "((8,), (536,))" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x.strides, y.strides" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**小步长更快？**\n", "\n", "![](http://www.scipy-lectures.org/_images/cpu-cacheline.png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- CPU从主内存中拉取数据到缓存块(pulls data from main memory to its cache in blocks)\n", "- 如果需要数据项连续操作适应于一个内存块（小步长）：\n", " - 需要更少的迁移\n", " - 更快\n", " \n", "**也可以看一下**：[numexpr](https://github.com/pydata/numexpr)设计用来减轻数组计算时的缓存效果。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 2.2.1.4.6 例子：原地操作（买者当心）\n", "\n", "有时，" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "`a -= b`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "并不等同于" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "`a -= b.copy() `" ] }, { "cell_type": "code", "execution_count": 90, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 0, -1],\n", " [ 1, 0]])" ] }, "execution_count": 90, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = np.array([[1, 2], [3, 4]])\n", "x -= x.transpose()\n", "x" ] }, { "cell_type": "code", "execution_count": 91, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 0, -1],\n", " [ 1, 0]])" ] }, "execution_count": 91, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y = np.array([[1, 2], [3, 4]])\n", "y -= y.T.copy()\n", "y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- `x`和`x.transpose()`共享数据\n", "- `x -= x.transpose()`逐个元素修改数据...\n", "- 因为`x`和`x.transpose()`步长不同，修改后的数据重新出现在RHS" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.2.1.5 剖析上的发现\n", "\n", "![](http://www.scipy-lectures.org/_images/threefundamental.png)\n", "\n", "- *内存块*：可以共享，`.base`、`.data`\n", "- *数据类型描述符*：结构化数据、子数组、字节顺序、投射、视图、`.astype()`、`.view()`\n", "- *步长索引*：跨越、C/F-order、w/ 整数切片、`as_strided`、广播、跨越技巧、`diag`、CPU缓存一致性" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2.2.2 通用函数\n", "\n", "### 2.2.2.1 他们是什么？\n", "\n", "- Ufunc在数组的所有元素上进行元素级操作。\n", "\n", "例子：\n", "\n", "`np.add`、`np.subtract`、`scipy.special`.*, ...\n", "\n", "- 自动话支持：广播、投射...\n", "- ufunc的作者只提供了元素级操作，Numpy负责剩下的。\n", "- 元素级操作需要在C中实现（或者比如Cython）\n", "\n", "#### 2.2.2.1.1 Ufunc的部分\n", "\n", "- 由用户提供" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "```C\n", "void ufunc_loop(void **args, int *dimensions, int *steps, void *data)\n", "{\n", " /*\n", " * int8 output = elementwise_function(int8 input_1, int8 input_2)\n", " *\n", " * This function must compute the ufunc for many values at once,\n", " * in the way shown below.\n", " */\n", " char *input_1 = (char*)args[0];\n", " char *input_2 = (char*)args[1];\n", " char *output = (char*)args[2];\n", " int i;\n", "\n", " for (i = 0; i < dimensions[0]; ++i) {\n", " *output = elementwise_function(*input_1, *input_2);\n", " input_1 += steps[0];\n", " input_2 += steps[1];\n", " output += steps[2];\n", " }\n", "}\n", "```" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "- Numpy部分，由下面的代码创建" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "```C\n", "char types[3]\n", "\n", "types[0] = NPY_BYTE /* type of first input arg */\n", "types[1] = NPY_BYTE /* type of second input arg */\n", "types[2] = NPY_BYTE /* type of third input arg */\n", "\n", "PyObject *python_ufunc = PyUFunc_FromFuncAndData(\n", " ufunc_loop,\n", " NULL,\n", " types,\n", " 1, /* ntypes */\n", " 2, /* num_inputs */\n", " 1, /* num_outputs */\n", " identity_element,\n", " name,\n", " docstring,\n", " unused)\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " - ufunc也可以支持多种不同输入输出类型组合。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 2.2.2.1.2 简化一下\n", "\n", "`ufunc_loop`是非常通用的模式，Numpy提供了预制\n", "\n", "| | |\n", "|---|---|\n", "|`PyUfunc_f_f`|`float elementwise_func(float input_1)`|\n", "|`PyUfunc_ff_f`|`float elementwise_func(float input_1, float input_2)`|\n", "|`PyUfunc_d_d`|`double elementwise_func(double input_1)`|\n", "|`PyUfunc_dd_d`|`double elementwise_func(double input_1, double input_2)`|\n", "|`PyUfunc_D_D`|`elementwise_func(npy_cdouble \\*input, npy_cdouble\\* output)`|\n", "|`PyUfunc_DD_D`|`elementwise_func(npy_cdouble \\*in1, npy_cdouble \\*in2, npy_cdouble\\* out)`|\n", "\n", " - 只有需要提供`elementwise_func`\n", " - ... 除非当你的元素级函数不是上面的形式中的任何一种" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " \n", "### 2.2.2.2 练习：从0开始构建一个ufunc\n", "\n", "Mandelbrot分形由如下迭代定义：\n", "\n", "![](http://www.scipy-lectures.org/_images/math/a8ced34b782ab71770c847a2aaa149b2dd41c7bc.png)\n", "\n", "$C=X+iy$ 是一个复数，只要$Z$仍然是有限的，无论迭代要跑多久，迭代都会重复。$C$属于Mandelbrot集。\n", "\n", "- ufunc调用`mandel(z0, c)`计算：" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "```python\n", "z = z0\n", "for k in range(iterations):\n", " z = z*z + c\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "比如，一百次迭代或者直到`z.real**2 + z.imag**2 > 1000`。用它来决定哪个`C`是在Mandelbrot集中。\n", "\n", "- 我们的函数是很简单的，因此，请利用`PyUFunc_*`助手。\n", "- 用Cython来完成\n", "\n", "\n", "提醒：一些预设Ufunc循环：\n", "\n", "| | |\n", "|---|---|\n", "|`PyUfunc_f_f`|`float elementwise_func(float input_1)`|\n", "|`PyUfunc_ff_f`|\t`float elementwise_func(float input_1, float input_2)`|\n", "|`PyUfunc_d_d`|\t`double elementwise_func(double input_1)`|\n", "|`PyUfunc_dd_d`|\t`double elementwise_func(double input_1, double input_2)`|\n", "|`PyUfunc_D_D`|\t`elementwise_func(complex_double *input, complex_double* output)`|\n", "|`PyUfunc_DD_D`|\t`elementwise_func(complex_double *in1, complex_double *in2, complex_double* out)`|\n", "\n", "类型代码：\n", "\n", "---\n", "\n", "NPY_BOOL, NPY_BYTE, NPY_UBYTE, NPY_SHORT, NPY_USHORT, NPY_INT, NPY_UINT,\n", "NPY_LONG, NPY_ULONG, NPY_LONGLONG, NPY_ULONGLONG, NPY_FLOAT, NPY_DOUBLE,\n", "NPY_LONGDOUBLE, NPY_CFLOAT, NPY_CDOUBLE, NPY_CLONGDOUBLE, NPY_DATETIME,\n", "NPY_TIMEDELTA, NPY_OBJECT, NPY_STRING, NPY_UNICODE, NPY_VOID" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "#### 2.2.2.3 答案：从0开始创建一个ufunc" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "```python\n", "# The elementwise function\n", "# ------------------------\n", "\n", "cdef void mandel_single_point(double complex *z_in, \n", " double complex *c_in,\n", " double complex *z_out) nogil:\n", " #\n", " # The Mandelbrot iteration\n", " #\n", "\n", " #\n", " # Some points of note:\n", " #\n", " # - It's *NOT* allowed to call any Python functions here.\n", " #\n", " # The Ufunc loop runs with the Python Global Interpreter Lock released.\n", " # Hence, the ``nogil``.\n", " #\n", " # - And so all local variables must be declared with ``cdef``\n", " #\n", " # - Note also that this function receives *pointers* to the data;\n", " # the \"traditional\" solution to passing complex variables around\n", " #\n", "\n", " cdef double complex z = z_in[0]\n", " cdef double complex c = c_in[0]\n", " cdef int k # the integer we use in the for loop\n", "\n", " # Straightforward iteration\n", "\n", " for k in range(100):\n", " z = z*z + c\n", " if z.real**2 + z.imag**2 > 1000:\n", " break\n", "\n", " # Return the answer for this point\n", " z_out[0] = z\n", "\n", "\n", "# Boilerplate Cython definitions\n", "#\n", "# You don't really need to read this part, it just pulls in\n", "# stuff from the Numpy C headers.\n", "# ----------------------------------------------------------\n", "\n", "cdef extern from \"numpy/arrayobject.h\":\n", " void import_array()\n", " ctypedef int npy_intp\n", " cdef enum NPY_TYPES:\n", " NPY_CDOUBLE\n", "\n", "cdef extern from \"numpy/ufuncobject.h\":\n", " void import_ufunc()\n", " ctypedef void (*PyUFuncGenericFunction)(char**, npy_intp*, npy_intp*, void*)\n", " object PyUFunc_FromFuncAndData(PyUFuncGenericFunction* func, void** data,\n", " char* types, int ntypes, int nin, int nout,\n", " int identity, char* name, char* doc, int c)\n", "\n", " void PyUFunc_DD_D(char**, npy_intp*, npy_intp*, void*)\n", "\n", "\n", "# Required module initialization\n", "# ------------------------------\n", "\n", "import_array()\n", "import_ufunc()\n", "\n", "\n", "# The actual ufunc declaration\n", "# ----------------------------\n", "\n", "cdef PyUFuncGenericFunction loop_func[1]\n", "cdef char input_output_types[3]\n", "cdef void *elementwise_funcs[1]\n", "\n", "loop_func[0] = PyUFunc_DD_D\n", "\n", "input_output_types[0] = NPY_CDOUBLE\n", "input_output_types[1] = NPY_CDOUBLE\n", "input_output_types[2] = NPY_CDOUBLE\n", "\n", "elementwise_funcs[0] = mandel_single_point\n", "\n", "mandel = PyUFunc_FromFuncAndData(\n", " loop_func,\n", " elementwise_funcs,\n", " input_output_types,\n", " 1, # number of supported input types\n", " 2, # number of input args\n", " 1, # number of output args\n", " 0, # `identity` element, never mind this\n", " \"mandel\", # function name\n", " \"mandel(z, c) -> computes iterated z*z + c\", # docstring\n", " 0 # unused\n", " )\n", "\n", "```" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "```python\n", "import numpy as np\n", "import mandel\n", "\n", "\n", "x = np.linspace(-1.7, 0.6, 1000)\n", "y = np.linspace(-1.4, 1.4, 1000)\n", "c = x[None,:] + 1j*y[:,None]\n", "z = mandel.mandel(c, c)\n", "\n", "import matplotlib.pyplot as plt\n", "plt.imshow(abs(z)**2 < 1000, extent=[-1.7, 0.6, -1.4, 1.4])\n", "plt.gray()\n", "plt.show()\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![](http://scipy-lectures.github.io/_images/mandelbrot.png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**笔记**：大多数模板可以由下列Cython模块来自动完成：\n", "http://wiki.cython.org/MarkLodato/CreatingUfuncs\n", "\n", "**一些可接受的输入类型**\n", "\n", "例如：支持小数点后一位及两位两个准确度版本" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "```python\n", "cdef void mandel_single_point(double complex *z_in,\n", " double complex *c_in,\n", " double complex *z_out) nogil:\n", " ...\n", "\n", "cdef void mandel_single_point_singleprec(float complex *z_in,\n", " float complex *c_in,\n", " float complex *z_out) nogil:\n", " ...\n", "\n", "cdef PyUFuncGenericFunction loop_funcs[2]\n", "cdef char input_output_types[3*2]\n", "cdef void *elementwise_funcs[1*2]\n", "\n", "loop_funcs[0] = PyUFunc_DD_D\n", "input_output_types[0] = NPY_CDOUBLE\n", "input_output_types[1] = NPY_CDOUBLE\n", "input_output_types[2] = NPY_CDOUBLE\n", "elementwise_funcs[0] = mandel_single_point\n", "\n", "loop_funcs[1] = PyUFunc_FF_F\n", "input_output_types[3] = NPY_CFLOAT\n", "input_output_types[4] = NPY_CFLOAT\n", "input_output_types[5] = NPY_CFLOAT\n", "elementwise_funcs[1] = mandel_single_point_singleprec\n", "\n", "mandel = PyUFunc_FromFuncAndData(\n", " loop_func,\n", " elementwise_funcs,\n", " input_output_types,\n", " 2, # number of supported input types <----------------\n", " 2, # number of input args\n", " 1, # number of output args\n", " 0, # `identity` element, never mind this\n", " \"mandel\", # function name\n", " \"mandel(z, c) -> computes iterated z*z + c\", # docstring\n", " 0 # unused\n", " )\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.2.2.4 广义ufuncs\n", "\n", "**ufunc**\n", "\n", "`output = elementwise_function(input)`\n", "\n", "`output`和`input`都可以只是一个数组元素。\n", "\n", "**广义ufunc**\n", "\n", "`output`和`input`可以是有固定维度数的数组\n", "\n", "\n", "例如，矩阵迹（对象线元素的sum）：" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "```\n", "input shape = (n, n)\n", "output shape = () i.e. scalar\n", "\n", "(n, n) -> ()\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "矩阵乘积：" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "```\n", "input_1 shape = (m, n)\n", "input_2 shape = (n, p)\n", "output shape = (m, p)\n", "\n", "(m, n), (n, p) -> (m, p)\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- 这是广义ufunc的”签名“\n", "- g-ufunc发挥作用的维度是“核心维度”\n", "\n", "**Numpy中的状态**\n", "\n", "- g-ufuncs已经在Numpy中...\n", "- 新的可以用`PyUFunc_FromFuncAndDataAndSignature`来创建\n", "- 大部分的线性代数函数用 g-ufuncts 来实现:" ] }, { "cell_type": "code", "execution_count": 119, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([-0.11206352, -0.02984284, -0.08754103])" ] }, "execution_count": 119, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import numpy as np\n", "np.linalg.det(np.random.rand(3, 5, 5))" ] }, { "cell_type": "code", "execution_count": 120, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'(m,m)->()'" ] }, "execution_count": 120, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.linalg._umath_linalg.det.signature" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- 也有一些 g-ufuncs 用来测试，ATM:" ] }, { "cell_type": "code", "execution_count": 121, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'(m,n),(n,p)->(m,p)'" ] }, "execution_count": 121, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import numpy.core.umath_tests as ut\n", "ut.matrix_multiply.signature" ] }, { "cell_type": "code", "execution_count": 122, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(10, 2, 5)" ] }, "execution_count": 122, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = np.ones((10, 2, 4))\n", "y = np.ones((10, 4, 5))\n", "ut.matrix_multiply(x, y).shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- 后两个维度成为了核心维度，并且根据每个*签名*去修改\n", "- 否则，g-ufunc“按元素级”运行\n", "- 这种方式的矩阵乘法对一次在许多小矩阵是非常有用\n", "\n", "**广义ufunc循环**\n", "\n", "矩阵相乘 `(m,n),(n,p) -> (m,p)`" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "```C\n", "void gufunc_loop(void **args, int *dimensions, int *steps, void *data)\n", "{\n", " char *input_1 = (char*)args[0]; /* these are as previously */\n", " char *input_2 = (char*)args[1];\n", " char *output = (char*)args[2];\n", "\n", " int input_1_stride_m = steps[3]; /* strides for the core dimensions */\n", " int input_1_stride_n = steps[4]; /* are added after the non-core */\n", " int input_2_strides_n = steps[5]; /* steps */\n", " int input_2_strides_p = steps[6];\n", " int output_strides_n = steps[7];\n", " int output_strides_p = steps[8];\n", "\n", " int m = dimension[1]; /* core dimensions are added after */\n", " int n = dimension[2]; /* the main dimension; order as in */\n", " int p = dimension[3]; /* signature */\n", "\n", " int i;\n", "\n", " for (i = 0; i < dimensions[0]; ++i) {\n", " matmul_for_strided_matrices(input_1, input_2, output,\n", " strides for each array...);\n", "\n", " input_1 += steps[0];\n", " input_2 += steps[1];\n", " output += steps[2];\n", " }\n", "}\n", "```" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## 2.2.3 互操性功能\n", "\n", "### 2.2.3.1 多维度类型数据贡献\n", "\n", "假设你\n", "\n", "1. 写一个库处理（多维度）二进制数据，\n", "2. 想要它可以用Numpy或者其他库来简单的操作数据，\n", "3. ... 但是并不像依赖Numpy。\n", "\n", "目前，三个解决方案：\n", "\n", "- “旧的” buffer接口\n", "- 数组接口\n", "- “新的” buffer接口([PEP 3118](http://www.python.org/dev/peps/pep-3118))\n", "\n", "### 2.2.3.2 旧的buffer协议\n", "\n", "- 只有1-D buffers\n", "- 没有数据类型信息\n", "- C-级接口；`PyBufferProcs` `tp_as_buffer`在类型对象中\n", "- 但是它被整合在Python中（比如，字符支持这个协议）\n", "\n", "使用PIL(Python Imaging Library)的小练习：\n" ] }, { "cell_type": "code", "execution_count": 98, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/linyong/anaconda/envs/pydata/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: the frombuffer defaults may change in a future release; for portability, change the call to read:\n", " frombuffer(mode, size, data, 'raw', mode, 0, 1)\n", " \n" ] } ], "source": [ "from PIL import Image\n", "data = np.zeros((200, 200, 4), dtype=np.int8)\n", "data[:, :] = [255, 0, 0, 255] # Red\n", "# In PIL, RGBA images consist of 32-bit integers whose bytes are [RR,GG,BB,AA]\n", "data = data.view(np.int32).squeeze()\n", "img = Image.frombuffer(\"RGBA\", (200, 200), data)\n", "img.save('test.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Q**：\n", "检查一下如果`data`修改的话，再保存一下`img`看一下会发生什么。\n", "\n", "### 2.2.3.3 旧的buffer协议" ] }, { "cell_type": "code", "execution_count": 100, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/linyong/anaconda/envs/pydata/lib/python3.6/site-packages/ipykernel_launcher.py:14: RuntimeWarning: the frombuffer defaults may change in a future release; for portability, change the call to read:\n", " frombuffer(mode, size, data, 'raw', mode, 0, 1)\n", " \n" ] } ], "source": [ "import numpy as np\n", "# import Image \n", "from PIL import Image\n", "\n", "# Let's make a sample image, RGBA format\n", "\n", "x = np.zeros((200, 200, 4), dtype=np.int8)\n", "\n", "x[:,:,0] = 254 # red\n", "x[:,:,3] = 255 # opaque\n", "\n", "data = x.view(np.int32) # Check that you understand why this is OK!\n", "\n", "img = Image.frombuffer(\"RGBA\", (200, 200), data)\n", "img.save('test.png')\n", "\n", "#\n", "# Modify the original data, and save again.\n", "#\n", "# It turns out that PIL, which knows next to nothing about Numpy,\n", "# happily shares the same data.\n", "#\n", "\n", "x[:,:,1] = 254\n", "img.save('test2.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![](http://scipy-lectures.github.io/_images/test.png)![](http://scipy-lectures.github.io/_images/test2.png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.2.3.4 数组接口协议\n", "\n", "- 多维度buffers\n", "- 存在数据类型信息\n", "- Numpy-特定方法；慢慢的废弃（不过不会消失）\n", "- 然而，没有整合在Python中\n", "\n", "**也可以看一下**：文档：http://docs.scipy.org/doc/numpy/reference/arrays.interface.html" ] }, { "cell_type": "code", "execution_count": 101, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'data': (140264175850144, False),\n", " 'descr': [('', '= 1903) & (year <= 1910))\n", " | ((year >= 1917) & (year <= 1918)))\n", "# '&' means 'and' and '|' means 'or'\n", "populations[bad_years, 0] = np.ma.masked\n", "populations[bad_years, 1] = np.ma.masked\n", "\n", "populations.mean(axis=0)" ] }, { "cell_type": "code", "execution_count": 117, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "masked_array(data = [21087.656489006717 15625.799814240254 3322.5062255844787],\n", " mask = [False False False],\n", " fill_value = 1e+20)" ] }, "execution_count": 117, "metadata": {}, "output_type": "execute_result" } ], "source": [ "populations.std(axis=0)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "注意，Matplotlib了解masked数组：" ] }, { "cell_type": "code", "execution_count": 124, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD8CAYAAACcjGjIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8VNX5+PHPk40kkJBASIBAWCRB2ZeIUFQUlMWVulXb\nX6XWlmr1p9VvsdpFrXu1X6201f5sxWLrhqiIICIFcakIhsWELYBsSQgkEBIghGxzfn/cO2FIZpJJ\nMpOZZJ7365XXzJw5996TyWSeueec+xwxxqCUUkq5Cgt0A5RSSgUfDQ5KKaUa0OCglFKqAQ0OSiml\nGtDgoJRSqgENDkoppRrwKjiIyD0iskVENovIGyISLSIDRGStiOwUkbdEJMqu28l+vMt+vr/Lfh6w\ny3NFZJpL+XS7bJeI3O/rX1IppVTzNBkcRCQVuAvINMYMA8KBG4E/AM8ZY9KBo8Ct9ia3AkeNMYOA\n5+x6iMgQe7uhwHTgBREJF5Fw4K/ADGAIcJNdVymlVIB4260UAcSISAQQCxQCk4GF9vPzgZn2/avt\nx9jPTxERscvfNMZUGmP2ALuAcfbPLmPMbmNMFfCmXVcppVSARDRVwRhTICJ/BPYDFcDHwHqg1BhT\nY1fLB1Lt+6lAnr1tjYiUAd3t8q9cdu26TV698vOaaldSUpLp379/U9WUUkrZ1q9ff9gY08Obuk0G\nBxFJxPomPwAoBd7G6gKqz5mHQzw856nc3dmL25weIjIbmA2QlpZGVlZWo21XSil1mojs87auN91K\nlwB7jDHFxphq4F3gO0CC3c0E0Ac4YN/PB/raDYkAugIlruX1tvFU3oAx5iVjTKYxJrNHD6+Cn1JK\nqRbwJjjsB8aLSKw9djAF2Ap8Alxn15kFvG/fX2w/xn5+lbGy+y0GbrRnMw0A0oF1wNdAuj37KQpr\n0Hpx6381pZRSLeXNmMNaEVkIbABqgI3AS8BS4E0Recwue9ne5GXgXyKyC+uM4UZ7P1tEZAFWYKkB\n7jDG1AKIyJ3AcqyZUPOMMVt89ysqpZRqLmmvKbszMzONjjkopZT3RGS9MSbTm7p6hbRSSqkGmuxW\nUkp1bIs2FvDM8lwOlFbQOyGGOdMGM3N0atMbqg5Ng4NSIWzRxgIeeDeHiupaAApKK3jg3RwADRAh\nTruVlAphzyzPrQsMThXVtTyzPDdALVLBQoODUiHsQGlFs8pV6NDgoFQI650Q06xyFTo0OCgVwuZM\nG4zUS2wTExnOnGmDA9MgFTR0QFqpEDZxUBLGQHiYUOswJMd14teXnaOD0UrPHJQKZZ/tKAbgf68f\nCcAvdRqrsmlwUCqErd5RTFKXTlw5sjddYyLZsO9ooJukgoQGB6VCVK3D8PnOYiZl9CA8TBiTlkCW\nBgdl0+CgVIjalFdK6clqLhpspb/P7N+NXUUnKD1ZFeCWqWCgwUGpEPVpbhFhAhekJwEwtl8iAOv1\n7EGhwUGpkLV6RzGj0xJJiI0CYGSfBCLCRLuWFKDBQamQdPhEJdn5ZVyUcXpFxZiocIamdmX9Xg0O\nSoODUiHJOYX1osHJZ5Rn9kvkm/xSqmocgWiWCiIaHJQKQatzi0nqEsXQ3vFnlGf2S6SyxsHmA2UB\napkKFhoclAoxtQ7DZzuLuTCjB2FhZ+bOGNvfGpTW6x2UBgelQsw3+c4prMkNnkuOiyatWyxZOu4Q\n8poMDiIyWEQ2ufwcE5FfiEg3EVkhIjvt20S7vojIXBHZJSLZIjLGZV+z7Po7RWSWS/lYEcmxt5kr\nUj8VWPu2aGMBE59axYD7lzLxqVUs2lgQ6CapELY6t5gwgQvtKaz1je2XSNa+o7TX9eWVbzQZHIwx\nucaYUcaYUcBY4CTwHnA/sNIYkw6stB8DzADS7Z/ZwIsAItINeAg4DxgHPOQMKHad2S7bTffJbxcE\nnCttFZRWYDi90pYGCBUon+YWMapvQt0U1vrG9kvk8IlK9pecbOOWqWDS3G6lKcC3xph9wNXAfLt8\nPjDTvn818KqxfAUkiEgvYBqwwhhTYow5CqwAptvPxRtj1hjrq8qrLvtq93SlLRVMjpyoJLugzG2X\nklOmPe6gXUuhrbnB4UbgDft+ijGmEMC+db7bUoE8l23y7bLGyvPdlDcgIrNFJEtEsoqLi5vZ9MDQ\nlbZUMPlsZzHGUJcyw52M5DjioiP0YrgQ53VwEJEo4Crg7aaquikzLShvWGjMS8aYTGNMZo8ent/c\nwURX2lLBxDmFdVjvrh7rhIUJY9ISWb+vpA1bpoJNc84cZgAbjDGH7MeH7C4h7Nsiuzwf6OuyXR/g\nQBPlfdyUdwhzpg0mKuLMl1lX2lKBUOswfLajmAvTG05hrS+zXyI7Dp2g7GR1G7VOBZvmBIebON2l\nBLAYcM44mgW871J+sz1raTxQZnc7LQemikiiPRA9FVhuP3dcRMbbs5RudtlXuzdzdCqXnmP1uAmQ\nmhDDk9cM1wVVVJvLzi/l6MlqJjXSpeRUd73Dfu1aClVeLRMqIrHApcDPXIqfAhaIyK3AfuB6u/xD\n4DJgF9bMplsAjDElIvIo8LVd7xFjjPO89Xbgn0AMsMz+6TCOV9aSkdKFj++ZFOimqBB2egpr08Fh\nVN8EwsOE9fuOcvHZngevVcflVXAwxpwEutcrO4I1e6l+XQPc4WE/84B5bsqzgGHetKW9qal1sH5v\nCd8do2cKKrBW7yhmZN8EEju7n8LqKjYqgqG948nScYeQpVdI+9nWwmOUV9Vy3oDuTVdWyk+OnKgk\nO7+UizK8PwsYk5bIprxSqms1CV8o0uDgZ+v2WN+8xg3oFuCWqFD2+c7DTU5hrS+zfyKnqh1sPXDM\njy1TwUqDg5+t3VNC/+6xpMRHB7opKoStzi2ie+cohqd6nsJaX2Y/6wuNXu8QmjQ4+JHDYfh6b4me\nNaiAcjgMn+087DYLa2N6do0mNSFGr3cIURoc/GhH0XFKT1breIMKqOyCMkrKq5rVpeSU2T+RrL2a\nhC8UaXDwIx1vUMFgdW4RInCBF1NY68vsl0jR8Uryj2q6l1CjwcGP1u4poXfXaPokaqoMFTirc4sZ\n2SeBbl5MYa1vrD3usF7HHUKOBgc/Mcawdrc13tDBlqdQ7UhJeRXf5Je2qEsJYHDPOOI6Rej1DiFI\ng4Of7DlczuETlYzT8QYVQJ/XZWFt2VXO4WHCqLQETd8dgjQ4+IlzvOG8gTreoAJndW4x3TpHMaIZ\nU1jrG9svkdxDxzl2SpPwhRINDn6ybk8JSV2iGJjUOdBNUSHKUZeFNalZU1jry+zXDWNg4/5SH7ZO\nBTsNDn6ydo+ON6jAyiko40h5VYu7lJxGpSUQJrB+r447hBINDn6Qf/QkBaUVjOuvXUoqcFbnFiMC\nF2a0bmGsLp0iOKdXvF4pHWI0OPjB6esbdDBaBc7qHUWMaOEU1voy+1lJ+Go0CV/I0ODgB+v2lBAf\nHcHZPeMC3RQVoo6WV7Epr5SLWnnW4DS2fzdOVtWyrfC4T/angp8GBz9YZ483tGYQUKnW+KxuCqtv\ngkNmP2tlOM2zFDo0OPhY0bFT7D5crikzVEB9mltMYmwkI/ok+GR/vRNi6N01WscdQogGBx9bt1fH\nG1RgORyGT3cUc2FGD8J9ePY6tn83TaMRQrwKDiKSICILRWS7iGwTkQki0k1EVojITvs20a4rIjJX\nRHaJSLaIjHHZzyy7/k4RmeVSPlZEcuxt5ko7nv+5bk8JsVHhDOsdH+imqBC1+YBzCqtvupScxqYl\nUFh2ioJSTcIXCrw9c3ge+MgYczYwEtgG3A+sNMakAyvtxwAzgHT7ZzbwIoCIdAMeAs4DxgEPOQOK\nXWe2y3bTW/drBc66PSWM7ZdIRLielKnAqJvC2oIsrI3JtKdmZ+n1DiGhyU8wEYkHLgReBjDGVBlj\nSoGrgfl2tfnATPv+1cCrxvIVkCAivYBpwApjTIkx5iiwAphuPxdvjFljrKTxr7rsq10pPVnF9oPH\nOU/HG1QArc4tYkRqV7p36eTT/Z7dM47YqHDtWgoR3ny9HQgUA6+IyEYR+YeIdAZSjDGFAPat8zLM\nVCDPZft8u6yx8nw35Q2IyGwRyRKRrOLiYi+a3rb0+gYVaKUnrSmsk1p5VbQ7EeFhjNYkfCHDm+AQ\nAYwBXjTGjAbKOd2F5I678QLTgvKGhca8ZIzJNMZk9ujh21NmX1i3p4SoiDBG9Gl5kjOlWuOznYdx\n+HAKa31j+3Vj+8FjnKis8cv+VfDwJjjkA/nGmLX244VYweKQ3SWEfVvkUr+vy/Z9gANNlPdxU97u\nrNtbwui+CURHhge6KSpErc4tIjE2kpE+msJaX2a/RBwGNmkSvg6vyeBgjDkI5InIYLtoCrAVWAw4\nZxzNAt637y8GbrZnLY0Hyuxup+XAVBFJtAeipwLL7eeOi8h4e5bSzS77ajdOVNawuaBMxxtUwDiz\nsF6Q7tsprK5G20n4dPGfji/Cy3r/F3hNRKKA3cAtWIFlgYjcCuwHrrfrfghcBuwCTtp1McaUiMij\nwNd2vUeMMc532O3AP4EYYJn9066s33cUh9HxBhU4Ww4c4/AJ309hdRUXHcngnvE6KB0CvAoOxphN\nQKabp6a4qWuAOzzsZx4wz015FjDMm7YEq7W7jxARJozp55/TeaWasjrX6tltbRbWpmT2S+S9jQXU\nOozfzlBU4OlkfB9Zt6eEYaldiY3y9mRMKd9YtLGAiU+t4n9X7CAyXPhi52G/Hm9sv0ROVNaw/eAx\nvx5HBZYGBx84VV3LN/mluiSoanOLNhbwwLs5dVctV9caHng3h0UbC/x2zLF1Sfi0a6kj0+DgAxv3\nl1Jda3QwWrW5Z5bnUlFde0ZZRXUtzyzP9dsx+yTGkBLfSa936OA0OPjA2j1HELHmgCvVlg54yHPk\nqdwXRITMfpqEr6PT4OAD6/aUcE7PeLrGRAa6KSrE9E6IaVa5r4ztl0hBaQWFZZqEr6PS4NBKVTUO\nNuw/quMNKiDmTBtMZPiZM4ZiIsOZM22why18I7O/jjt0dBocWimnoIxT1Q4db1ABMXN0KoNT4ggT\nKw9NakIMT14znJmj3aYn85lzesUTExmu4w4dmM67bCVnsr1z+2twUG3v+KlqdhSd4OYJ/Xn4qqFt\ndtzI8DBG9U3QM4cOTM8cWmntniMMSu7i8/TISnlj1fYiqmocXDGiV5sfO7N/IlsLj1GuSfg6JA0O\nrVDrMGTtParrRauAWZJdSM/4aMakJTZd2cfG9Euk1mH4Jk+T8HVEGhxaYVuhlbpYxxtUIBw/Vc2n\nucXMGN6TsACksRiTlogIZGnXUoekwaEV1tYt7qPBQbW9/2w7RFVtYLqUALrGRJKRHKfBoYPS4NAK\n6/YcIa1bLL26+ndOuVLuLM0+SK+u0Yzu2/ZdSk7dukTy+Y5iBty/lIlPrfJr2g7VtjQ4tJDDYVi3\np0TPGlRAHDtVzWc7ipkxrFdAupTAyuuUtfcoBmvpxoLSCr/ndVJtR4NDC+0qPsHRk9UaHFRA/Ger\n1aV0eYC6lMDK61Rde+aKvv7O66TajgaHFnKON4zXxX1UAHyYU0jvrtGM7hu49UMCkddJtR0NDi20\nbk8JPeOj6dtNxxtU2yqrqOazHYeZMTxwXUrgg7xO2QvguWHwcIJ1m73Ah61TraXBoQWMMazdfYRx\nA7phLXutVNsJhi4lsPI6xUSGn1HmdV6n7AXwwV1QlgcY6/aDuzRABBGvgoOI7BWRHBHZJCJZdlk3\nEVkhIjvt20S7XERkrojsEpFsERnjsp9Zdv2dIjLLpXysvf9d9rZB/Ym778hJio5X6niDCogPcwpJ\nTYgJaJcSWHmdnrxmOAl2NuKU+E7e53Va+QhU1+t+qq6wylVQaM6Zw8XGmFHGGOda0vcDK40x6cBK\n+zHADCDd/pkNvAhWMAEeAs4DxgEPOQOKXWe2y3bTW/wbtQFnPqXxmolVtbGyimo+21nMjGE9g+Ks\ndeboVObdci4Aj81sRsK/svzmlas215pupauB+fb9+cBMl/JXjeUrIEFEegHTgBXGmBJjzFFgBTDd\nfi7eGLPGGGOAV132FZTW7imhW+cozurRJdBNUSFmxdZDVNeagHcpuUpPtv4Pdhw67t0GJ4ogLNz9\nc137+KhVqrW8DQ4G+FhE1ovIbLssxRhTCGDfJtvlqUCey7b5dllj5fluyoPWur1HGNdfxxtU23N2\nKY0KcJeSq7joSFITYrwLDscOwCuXgREIr5esMjIGpjzon0aqZvM2OEw0xozB6jK6Q0QubKSuu09M\n04LyhjsWmS0iWSKSVVxc3FSb/eJAaQV5JRU63qDaXNnJaj7fWcxlw4OjS8lVRkoXcg82ERxK98Mr\nM+B4IfxoMVz9F+ja13ouLAKunAsjbvB/Y5VXvAoOxpgD9m0R8B7WmMEhu0sI+7bIrp4P9HXZvA9w\noInyPm7K3bXjJWNMpjEms0ePHt40/QxLdy9l6sKpjJg/gqkLp7J099Jm72Od5lNSAfLx1oN2l1Lv\nQDelgYyecewuLqem1uG+wpFvrTOGk0fhh4ug33esQHDPZrj0UXDUQP/z27bRqlFNBgcR6Swicc77\nwFRgM7AYcM44mgW8b99fDNxsz1oaD5TZ3U7LgakikmgPRE8FltvPHReR8fYspZtd9uUzS3cv5eEv\nH6awvBCDobC8kIe/fNjrAOEMLL/Lnk7coKfYU/G5r5uoVKOcXUoj+3QNdFMaGJwSR1Wtg71HTjZ8\nsniHFRiqymHWYuh77pnPD5xk3e75zP8NVV7zZiW4FOA9+zQ2AnjdGPORiHwNLBCRW4H9wPV2/Q+B\ny4BdwEngFgBjTImIPAp8bdd7xBhTYt+/HfgnEAMss3986vkNz3Oq9tQZZadqT/HkuiepcdQQERZB\neFg4kRJJeFg4EWERVpmEs65wHfM2z6PKUWVtGFnKI1/9nrAw4fKBl/u6qUo1YHUpHebW8wcEXZcS\nQEZKHAA7Dx1nULLLRI1DW+DVq637P1oKKUMabpwyHGK6we5PYeSNbdBa5Y0mg4MxZjcw0k35EWCK\nm3ID3OFhX/OAeW7Ks4BhXrS3xQ6WH3RbXlZZxm//+9tm7+9U7Sme3/C8BgfVJpZvPUiNw3DZ8OCZ\npeRqUHIXRCD30HFmONt4YCP867sQEQ03L4YeGe43DguDARfAnk/BGAjC4BcMlu5eyvMbnudg+UF6\ndu7J3WPu9uvnT8isId2zc08KywsblCfHJvPP6f+k1lFLjaOGWmPd1pga67Gjlls/vtXtPgs9BByl\nmqupf/yl2YX0SYxhRBB2KQFER4bTr1vs6RlLeV/Dv6+F6HirK6nbwMZ3MOBC2Pq+NTaRNMj/DW5n\nnN3izt4PZ7c44LcAETLB4e4xd5/x4gJEh0dz79h76RvXt5EtQWoSMRENFzSRmuCZTqjar6b+8UtP\nVvHfXYe59YLg7FJyykiJY8ehE7D3v/D6DdC5B8z6ABIa//8CYMBF1u2e1Roc3PDULe7P3ouQya10\n+cDLefg7D9Orcy8EoVfnXjz8nYe9emErDk3FOCIblFed6O+HlqpQ09g/PsDHWw5R4zBcHqRdSk6D\ne8bR68hXmH9fC/G94ZZl3gUGgO5nQXyqDkp74Klb3FO5L4TMmQNYAaIlUTY57DscKoROPZYjkaWY\n6gQctbFEdt3EFwVfcH6qTsFTjXcN1ThqKDxRyP7j+9l3bB95x/PYd2wf+4/vd9vdCaf/8ZfkFNK3\nWwzDU4OwSyl7gZUPqSyfO6MSCYsopTI+g+gfLYEuzZhuLgIDJsGOj8DhsMYhVJ3E6ERKTpU0KO/Z\nuaffjhlSwaGl5kwbzC/fPkX5sdF1ZTFRNfQcOo85n87htcteY2BCE32qqkNz1zX0my9+wys5r1Dp\nqCT/RD41jpq6+jERMaTFpZGRmMGRiiOcqD7RYJ9R4VGsP7CVL3cd5icXDAy+LiVnZlU7gV6nqhIc\nCJv7fo/RzQkMTgMnwTevw6Ec6NVgDkzI2lS0ieOVxxEE43J9cHR4NHePudtvxxVrclH7k5mZabKy\nstrseOOf+A8l5dVU1zronRDDnGmDOS89jJuW3kRsZCyvX/Y6CdE6BhGqpi6c6vYMICIsgov7Xkxa\nXBr94vuRFp9GWlwaSTFJdR/29QMLQIRY06grHVVUl47hxcsf4OJBXqTCbkvPDbNTbp+pLKonXX/d\ngtXgjhXCs2dbF8VNvMsHDWz/cktyuWX5LSR2SuSHQ37IvM3zWjVbSUTWuyRPbZSeOXhhd/EJDh6r\n5MErhvDj8wec8dzzk5/nxx/9mHtW38NLl75EZHjDsQnVsR2vOu6xa6jWUcuzFz3b6PbOf/D6XVIT\ne0/kmjcepbjrKn655vvcdPQmfjL8J8HzJcRDBtX4qkMt2198L0jKsKa0anBg/7H9/GzFz4iJiOHv\nU/9O7y69ufHstrsORDv2vLBss9X3O31Yw/69kT1G8vuJvyfrUBaPr32c9nom1pH4Ik2Kt74++DXX\nLr7W4/Pe9glfPvByPr7uY7JnZfPxdR9z+cDLMbWx7Ns5he/2mMuMATP417Z/MePdGfw9+++crHZz\nJXJb85BB9ZAktXyfAy6EfV9CTVXL99EBHCo/xOwVs6k1tfz9UiswtDUNDl74aPNBRvVN8Lj84RUD\nr+Cnw3/KOzvf4bVtr7Vx65Sr1qZJ8VZVbRXPZj3LrctvJTIskp+P/DnR4dFn1Gltn/DyLQepdRhu\nGDWCx85/jHeufIfMnpnM3TiXK967ggW5C6h2VLf2V2m5KQ9CeNQZRdVh0TxZdT0VVbUt2+eASVB9\nEgrarss42JSeKuVnK35GaWUpf7vkbwEbz9Tg0IS8kpPkFJQxw81Zg6s7R9/JlLQpPJP1DJ/na96l\nQDhVc4qn1j3V6LRQX9hxdAc3Lb2JV7a8wrUZ1/L2lW9z+6jbWzxV2pOlOYX06x7L0N7xAAxKHMSf\nJ/+Z+dPn0yeuD49+9SjXvH8NT659ss3OlM4w4gZIm4CVWFmga1+2Zj7K+7Xns6uo4QC7V/qfb+0r\nRKe0lleXc/t/bifveB5/nvxnhiYNDVhbdMyhCR/ZXUozhjU+xzxMwnji/CeY9dEs7vvsPv592b85\nK+GstmhiyMstyeWdne+wZPcSjle5TxtdWF7I0VNHSYxOdPu8NxzGwb+2/ovnNzxPXFQcf5n8Fyb1\nnVT3fEunSrtTUl7Fl98e4WcXNpylNCZlDPOnz2d13moe++oxXt/+et1zbXHl7BkqjlpdQbMWAxBX\nfAI++5TcQ8cZ3pKruWO7WTOVdn8KF93fdP0A8Fcai8raSu5adRfbSrbxp4v/xLk9z216Iz/SM4cm\nLNtcyJBe8aR1j22ybmxkLH+e/Gc6hXfizpV3cvRUw6uqlW+UV5ezcMdCblpyE9d9cB0Ldyzk/NTz\n6R7d3eM2l7x9CQ/+90G2l2xv9vEKTxTy049/yh+z/sjE1Im8e9W7ZwQGX3N2KXnKpSQiXJx2MeFu\nVlTz9ZmSR9UVVmK91LF1Rf26dyYqIsz7VeHcGTgJ8r+2srgGGX91W9Y4apjz6RzWHVzHoxMf5aK+\nF/mmwa2gwaERB8tOsWF/aZNdSq56du7J3MlzKTpZxD2r76G6NoB9wh2MMYbs4mwe+vIhLl5wMb9f\n83tO1Z7iV+f+ilXXr+LpC59mzrlz3Pb9/2LML5g5aCYf7f2I6z+4nh999CNW7FtxxrUHno65ZPcS\nrl18LTmHc/j9d37P3Ivn0j3GcxDyhaXZhfR36VLyxNMVsoXlheSWtGA6aXMUfgOmFvqcnhkZHiYM\n6tGldcFhwCRwVMO+NT5opO+cqjnF018/7fNuS4dx8NCXD/FJ3ic8MO4BrjzrytY21Se0W6kRy7fY\nXUrNTFswoscIHpn4CPd/fj+PrX2Mhyc8HHwXMAWx+qftPx3+U6ocVbyz8x12Ht1JTEQM0/tP55r0\naxjZY+QZr62naaHO8rvG3MWiXYt4Y/sb3Lv6Xnp27smNg2/k2vRrSYhOOOPYybHJpMSmkH04m1E9\nRvHE+U/QN97LdBCtcOREJWt2H+G2SU1f+OYpoaQgXPfBdUzvP52fj/o5A7oOcLN1KxWst25dzhzA\nSqOxdveRlu83bQKERVp5ltIvafl+WslhHGwr2cZXB75iTeEaNh7aeDptfz2F5YXsOLqDjEQPmWc9\nMMbw9NdPs/jbxdwx6g6+f873fdF0n9CL4Bpx40trOHKiihX3tqz7YO6Gufw95+/MyZzDzUNv9nHr\nOiZ3F4Q5Dek+hGvTr+WyAZfRJaqLm629V+uo5dP8T3l92+usPbiWTuGdGJ40nOzD2VTVnvkBMLXf\nVP5w4R+ICGub71Kvr93Pr9/LYeld5zO0d+P99u5er+jwaO479z4Kywv597Z/U1lbyZUDr+T2UbeT\n2sWHy7O/fQvkrYN7t5xR/MLqXTz9US45D08lLrqF1/28chlUHofb/DO5w9O4QcGJAtYcWMNXhV+x\ntnAtpZWlAGQkZjC+13iW7F7iNo2FU0ZiBpcPvJzLBlzm1TTmFze9yAvfvMAPh/yQOZlz/P4lUi+C\n84HDJypZt6eEOy9ueYbIO0ffyZ6yPTyT9Qz/2PwPSk+Vtkke9vbMXRI6gKSYJN664i2fHSc8LJzJ\naZOZnDaZnUd38sb2N3h7x9tu6+YczmmzwACwNOcAA5I6M6RX411K0PSZ0g/O+QEvb36Zt7a/xdI9\nS7k2/Vpmj5hNcmxy6xtakAV9xjYoHmwv/LPj0AnG9mvhBIABk2D1k3CyxBqk9iFPqU6eXvc0JZXW\nB39yTDIX9rmQCb0nML7XeJJirGs3hnQf4jYY//LcX+IwDpbuXspz65/jT+v/RGbPTK4YeAWX9LuE\n+Kj4umM7/1ZxUXEcqzrG1WddzS8zfxl0vQsaHDz4eMshHAamNzFLqTFhEsYFqRewcv/KusHpNp9N\n0o4YYzxeaXykohXdFE1IT0znwQkPsnDHwjNy1zj5M/NlfYdPVLLm2yP8/KJBXn9YNDZLqntMd+47\n9z5uHnIzL2W/xDs73mHRrkXcdPZN/HjYj0mMTmzZ7JsTxVC6H879SYOnMuqCw/GWB4eBk2D1E7D3\ncxhydcvEw2x7AAAgAElEQVT24YG7LyC1ppbymnLuH3c/E3pNYEBX9+nRmwrGN519E3nH8liyZwkf\n7v6Qh758iMe/epxJfSeREpvCwh0L6459rOoYYRLGuJ7jCJPgG/7V4ODBss3WHPNzesW1aj9/y/5b\ngw+ctlxFrq1Xj2qp/OP5dUHTHX9mn3Q9hrvg1BbHdlq+5SAOg89XfOvZuScPTniQW4bewovfvMj8\nLfNZkLuA8b3G898D/6WythJoxpeXuvGGhj0UqQkxxEaFt25QOnUsRHWxprT6ODh4+gJSVVvFD875\nQZPbNzVluW98X24feTu3jbiNLUe2sHT3Uj7c86Hb7iiHcfCXTX/hqkFXef8LtBGvw5WIhIvIRhFZ\nYj8eICJrRWSniLwlIlF2eSf78S77+f4u+3jALs8VkWku5dPtsl0iEvDJzWUnq1nz7RFmDOvV6lO9\nQORhd2qrq4Vbw2EcvLH9Da5ZfA2bj2zmu4O+6/Mrjb1195i7A3bsRRsLmPjUKn7z3mYiwoTcg8f8\ncpy+8X154oIneO/q95iYOpFVeavqAoOTV7NvCtaDhEHvUQ2eCgsT0lPiWhccwiOh33esPEs+tDpv\ntcfnfP0lQEQYljSMX437FSuvX+mxXluemTZHc85l7ga2uTz+A/CcMSYdOAo419K8FThqjBkEPGfX\nQ0SGADcCQ4HpwAt2wAkH/grMAIYAN9l1A2bFNmtxleZMYfXE0xvOYHg261nKq/03l7upRWSa4u8c\nRXnH8rh1+a08sfYJRieP5r2r3uORiY/4/Epjb7VmQajWWLSxgAfezaGg1Ep9XeMw/Pq9zSzaWOC3\nY56VcBbPXvQsgvsvP4XlhY2/NwuyIHkIRHV2+3RGchdyD7bwKmmnAZPgyC4o883r8N7O9/jFJ7+g\nT+c+dArvdMZz/v4SEBEWQa/O7s8G2/LMtDm86lYSkT7A5cDjwL1ifZ2eDDjnXc0HHgZeBK627wMs\nBP5i178aeNMYUwnsEZFdwDi73i5jzG77WG/adbe26jdrhWU5haQm+Ga9XnfLk3YK78Tw7sN5Zcsr\nLNm9hHvG3sMVA6/w6YDUliNbPJ4+F5YXcsMHN9A3ru8ZaaTT4tPoHt0dEfHrmrXOs4XnNzxPuITz\n++/8nu8O+m7d7+/LK42bKxDHfmZ5LhXVZ+Yiqqiu5Znlucwc7cPZRW546koDOP+N8xnRYwQTek9g\nQu8JDO0+1BqYdzisM4chMz3ud3DPON5en09JeRXdOkd5rNeogfYswT2fwqjWTfGct3kez61/jgm9\nJvCni//EJ3mftHl3q6elitvizLQlvB1z+BNwH+DsgO8OlBpjnFcQ5QPOd3EqkAdgjKkRkTK7firw\nlcs+XbfJq1d+XjN+B586fqqaz3ce5v+M7+eTD+vGBrCyi7N5cu2T/PqLX/P2jrf59Xm/5uxuZ7f4\nWNWOalbuX8nr215nY9HGBouDOMVGxNI9pjvbS7azcv9Kas3pD6bOkZ1Ji0tjT9kev6xZu+/YPh78\n74NsKNrA+ann89CEh4L2m1NbOWCfMXhb7kuePrB+cM4PMBjWHFjDC5te4K+b/kpcZBzn9jyXCfGD\nmFBTTlrvMXzoYUzLdVB6/MAWXjCYPBRiu1vjDi0MDg7j4NmsZ5m/dT7T+0/nifOfIDI8MiBfApoa\nzA42TQYHEbkCKDLGrBeRi5zFbqqaJp7zVO6ua8vtxRciMhuYDZCWltZIq1tu1fYiqmodzBjuuw8s\nT2/EET1G8Nrlr7Fo1yKe3/A831vyPa7PuJ47R93ZrJz9JadKWLhjIW/lvkXRySL6dOnDfefeR0xE\nDH9Y94cG//gPTniwrj3VjmoKTxTWLVm5/9h+9h3fx7aSbW6PVVheyEd7P2JE0gir+8XLAFrrqOW1\nba8xd+NcosKjeGziY1x11lVBN30vEHonxNR1KdUv97emPrDuGXsPR08dZe3BtXx14Cu+KvyKVXmr\noG9vuub+P07UVtR9uXA9uxybMgWAnY0EhyYnS4SFQf8LrDMHY6ylRJuh2lHNw18+zOJvF3PT2Tdx\n/7j7Az4rKJBnxc3V5EVwIvIk8EOgBogG4oH3gGlAT/vsYALwsDFmmogst++vEZEI4CDQA7gfwBjz\npL3f5ZzufnrYGDPNLn/AtZ4n/roI7vZ/r2f9vqN89cAUwsLa7oPrWNUxXtj0Am9uf5MuUV24a/Rd\nXJt+rdvcOU7bjmzjtW2vsWzPMqocVUzoNYEfnPMDzk89v267ls5W8rSymavu0d0ZnjSc4T2GMzxp\nOMOShhEXFdfguEkxSURHRJN3PI9JfSbx4IQHfTPPvoNwjjm4di3FRIbz5DXD/d6t1FzGGPKW3MGa\nb5fxx6Qkt9ekACTHJlNUGk5STCKjUnvTtVPX0z9RXdlZupO3c98+44rj6PDohmM8WfNgyT1wZxYk\npXvdzoqaCn756S/5LP8z7hx1J7NHzNYvIjTvIrhmXSFtnzn80hhzhYi8DbxjjHlTRP4GZBtjXhCR\nO4DhxpjbRORG4BpjzA0iMhR4HWucoTewEkjHOqPYAUwBCoCvge8bY7Y0aIALfwSHiqpaxjy6guvG\n9uHRmcN8um9v7Ti6gyfXPknWoSzO6XYOk/pM4v1v36/7cL9j1B10iujEG9veYEPRBmIiYrjqrKv4\n/tnf92ned09X3v52/G8ZlDiInOIccg7nkF2czd5jewErZcOArgNI6JRAzuGcBmsNXJ9xPb8b/zv9\nJ3Vj0cYCnlmey4HSirplaIMtMNR56WKIjGVEWJ7bbkuA7w76Liu276FWykntbiirLKO0srTJ9Sd6\nde7Fx9d9fLrgyLfw5zFw2R9h3E+9al5ZZRl3rLyDnMM5/Oa833DD4Bu8/tU6ura6QvpXwJsi8hiw\nEXjZLn8Z+Jc94FyCNUMJY8wWEVmANdBcA9xhjHU+KiJ3AsuBcGBeU4HBXz7dUURFda1PZim1VEZi\nBvOmzWP53uU8uuZR/pb9t7rnCssL+e1/fwtAapdU5mTOYWb6zLqrL32pqe6God2HcqP1p6Wssowt\nh7eQfTibnMM5fJH/BQ4cDfb5RcEXGhg8mDk6NXiDgavqU3AwByb8nJ5l1W7PLnt17sUjEx+hpiiH\nD3MKefeWSxERjDFU1FRwrOoYUxdOdRtYCssLeX/X+0zrP43oiGjoNhC69rW6lrwIDgfLD3LbitvY\nf3w/f5z0Ry7td6lPfu1QpLmVXNz1xkY+31nM17+5hIjwwF+xeOnbl3LwZMM50ImdEvnkhk8a7XIK\npBHzR7j9xxeE7FnZAWiR8pn8LPjHFLjhVZZGR7g9u3R2Dc3/ci8PLd7Cul9PITn+zOtHPHVbhks4\ntaaWrp26MvOsmdww+AbSVv0Bti+B+3ZDI+/53WW7uW3FbRyrOsbci+cyrtc4j3VDVXPOHAL/CRgk\nKmtqWbW9iKlDegZFYAA4dNL9Qu2llaVBGxjA87ztUJ+V1CHk21/IUjObvC4kPcVKjpjr5mI4Txcc\nPjbxMV6e+jLjeo7jtW2vcfl7l/MzRz4rwyqpObDRY7NyinOYtWwWlbWVvDLtFQ0MPqDpM2xf7DzM\nicoan85Saq1gSOfQEu1tPrdqhoL1ENcLulpdYI3NvnFNwHdBeo8znmuq23Jcr3EUnyzmnZ3v8Pb2\nt/hFSg9SVt/JdcN/xLXp17Lu4Lq6bROjEzleeZzkzsm8dOlLpMX7ZyZjqNHgYFu2+SBx0RF856yk\nQDelTnv9kG1v87lVMxRkNVi/wZPuXTqR1CWKHQfdp9Foalpnj9ge3DbyNn4y/Cd8+tJ5LIhx8NdN\nf+WFTS8gIjiMNa5VcqoEQZg1ZJYGBh/S4ABU1zpYsfUQl56TQlREcHQpQfv+kG1P87mVl06WQMlu\nGP1DrzdJT45jR1ErcixhpZ6YkjaZKRv+xb6ff8H3lv0fymvOTO1hMLyy5RVuOuemVh1LnabBAVjz\n7RHKKqqZHsBZSp7oh6wKGgUbrNs+Xo1nAnYajaw8jDGtm6k2YBKse4l+xw5xsuak2yrBmsCuvQqe\nr8kBtGzzQWKjwrkwo0fTlZUKVQVZgECvhplYPUlP6UJ5Va3bK8Cbpf/5VhbYPZ/qhIc2EvLBodZh\nWLH1IJPPTiY6MnhnACkVcAXrocfZEO39dTXOQemdh1qZoTUmwQpKuz8NaGr1UBLyweHrvSUcPlHF\njFas+KZUh2eMNY3Vy8Fop3Q7OLibztpsAydBQRaXp04KWFr3UBLyYw4fbT5Ip4gwLhqsXUpKeXR0\nD1SUuF0zujFdYyLpGR/tccZSswy4EL54DvZ9yeUZOhbnbyF95uBwGJZtLmRSRg86dwr5OKmUZ87B\naDfLgjYlo2frZywB0Hc8hEf5fHU45V5IB4eNeaUcOlYZVBe+KRWU8rMgIsZa/a2ZBqd0YeehE9Q6\nWpmqJyoW+p6nwaGNhHRw+GhzIZHhwuSzUwLdFKWCW8F6a73o8OafYaenxFFZ4yCvxP0U1GYZMMlK\n/Fd+pPX7Uo0K2eBgjGHZ5oOcPyiJrjGRgW6OUsGrpgoKv2n2YLTTYF8PSgPs/az1+1KNCtngsOXA\nMfKPVugsJaWacmgz1Fa2ODgMSrYS8PlkULr3GIiKs5YOVX4VssHhw5xCwsOES4dol5JSjSpYb902\n48poV507RdC3Www7ilp5rQNY3Vr9J+q4QxsIyeBgjOGjzQcZP7AbiZ2jAt0cpYJbwXro3MNadKeF\nBqfE+ebMAawprSW7oTTPN/tTboVkcNhx6AS7D5czXbuUlGpawXprCmsrciOlp8Sx+/AJqmsbrhDY\nbAPscQc9e/CrkAwOyzYXIgLThmqXklKNqiiFwzuaffFbfYNT4qiuNew9XN505aYkD4HYJNijg9L+\n1GRwEJFoEVknIt+IyBYR+b1dPkBE1orIThF5S0Si7PJO9uNd9vP9Xfb1gF2eKyLTXMqn22W7ROR+\n3/+aZ/po80HO7deN5LjopisrFcoOOC9+a11wyPDljKWwMKtrafenVloP5RfenDlUApONMSOBUcB0\nERkP/AF4zhiTDhwFbrXr3wocNcYMAp6z6yEiQ4AbgaHAdOAFEQkXkXDgr8AMYAhwk13X5xZtLGDc\n4/9h+8Hj5B46xqKNBf44jFIdh3MwuveYVu1mYI/OhInVpesTAyfBiYPWWY3yiyaDg7E4/6KR9o8B\nJgML7fL5wEz7/tX2Y+znp4iVyP1q4E1jTKUxZg+wCxhn/+wyxuw2xlQBb9p1fWrRxgIeeDeHouOV\nAJRV1PDAuzkaIJRqTP566J5uZUVthejIcPondfbdoHSV3T3113Hw3DDIXuCb/ao6Xo052N/wNwFF\nwArgW6DUGFNjV8kHUu37qUAegP18GdDdtbzeNp7KfeqZ5blUVNeeUVZRXcszy3N9fSilOgZjrDOH\nFk5hrS8jOY4dvuhWyl4Aqx49/bgsDz64SwOEj3kVHIwxtcaYUUAfrG/657irZt+6m9JgWlDegIjM\nFpEsEckqLi5uuuEuDnhYbMRTuVIhrywPyotaPd7glNEzjr1HyjlV70tas618BKrr/d9WV1jlymea\nNVvJGFMKrAbGAwki4ky00gc4YN/PB/oC2M93BUpcy+tt46nc3fFfMsZkGmMye/RoXort3gkxzSpX\nKuTlZ1m3PgoOg1PicBj4triV4w5l+c0rVy3izWylHiKSYN+PAS4BtgGfANfZ1WYB79v3F9uPsZ9f\nZYwxdvmN9mymAUA6sA74Gki3Zz9FYQ1aL/bFL+dqzrTBxNRb6S0mMpw50wb7+lBKdQwF6yG8E6QM\n88nuMlLsNBqt7Vrq2qd55apFvEmx2AuYb88qCgMWGGOWiMhW4E0ReQzYCLxs138Z+JeI7MI6Y7gR\nwBizRUQWAFuBGuAOY0wtgIjcCSwHwoF5xpgtPvsNbTNHW8MYzyzP5UBpBb0TYpgzbXBduVKqnoL1\n0GsERPgmi0D/pM5EhkvrZyxNedAaY3DtWoqMscqVzzQZHIwx2cBoN+W7scYf6pefAq73sK/Hgcfd\nlH8IfOhFe1tl5uhUDQZKeaO2Bg5sgrE/8tkuI8PDOKtHl9bPWBpxg3W78hGrK6lrHyswOMuVT+jy\nZ0qphoq2Qk2Fz2YqOaWnxLEp72jrdzTiBg0GfhaS6TOUUk0ocA5Gt+7it/oGp3Qhr6SC8sqapiur\ngNLgoJRqqGA9xHaHxAE+3W26nUZjpy/Sdyu/0uCglGoof701hbUVmVjdca4K55OL4ZRfaXBQSp2p\n8jgUb/fZ9Q2u+naLJToyzHdpNJTfaHBQSp3pwEbAWGs4+Fh4mDAouYtvsrMqv9LgoJQ6U75/BqOd\nMlLi2Omr7KzKb3Qqq1LqTAXrodtAiO3ml91npMTx7oYCyiqq6RoT6ZdjdESLNha06UW8euaglDqT\nc1lQP3EOSu/UriWvOZccKCitwAAFpRV+X3JAg4NS6rRjB+B4oV8Go50yevpwVbgQEYglBzQ4KKVO\nc443+PjKaFe9u0bTpVOEzlhqhkAsOaDBQSl1WkEWhEVCz+F+O4SIkJ7SxXdLhoaAlHj36937c8kB\nDQ5KqdMKNliBIaKTXw8zOMVHq8KFAGMMSV0aZsb195IDGhyUUhZHrXWNgx+7lJzSU+I4Ul7F4ROV\nfj9We7cku5DNB45x1YhepCbEIEBqQgxPXjPcr7OVdCqrUspSnAtVJ/w6GO3kmkYjqYt/z1LasyMn\nKnlo8RZG9unKs98bRUR4232f1zMHpZSlLhOr/88c6laF00HpRj38wVaOn6rm6etGtmlgAA0OSimn\n/CyI7grdz/L7oXrEdSIhNpIdmp3Vo+VbDvLBNwe4a3I6g+3pv21Jg4NSylKwwS+ZWN0RETJS4vTM\nwYPSk1X8dtFmhvSK57aL/B+s3dHgoJSCqnIo2tImXUpOGSlWAj5jTJsds714dMk2jpZX8fR1I4hs\n4+4kpyaPKiJ9ReQTEdkmIltE5G67vJuIrBCRnfZtol0uIjJXRHaJSLaIjHHZ1yy7/k4RmeVSPlZE\ncuxt5oq0wVcXpdRphd+AcbTJYLTT4JQ4jp+q4dAxnbHk6pPcIt7ZkM/tF53FsNSuAWuHNyGpBvgf\nY8w5wHjgDhEZAtwPrDTGpAMr7ccAM4B0+2c28CJYwQR4CDgPGAc85Awodp3ZLttNb/2vppTyWhtc\nGV1fRoqm0ajv2Klqfv1uDunJXbhz8qCAtqXJ4GCMKTTGbLDvHwe2AanA1cB8u9p8YKZ9/2rgVWP5\nCkgQkV7ANGCFMabEGHMUWAFMt5+LN8asMdb55asu+1JK+Vv2AvjkCev+SxdZj9vALnsweta8dUx8\napVfk8i1F09+uJ1Dx07xzPUj6RQRHtC2NKszS0T6A6OBtUCKMaYQrAACJNvVUoE8l83y7bLGyvPd\nlLs7/mwRyRKRrOLi4uY0XSnlTvYC+OAuqLFz9JTlWY/9HCAWbSzgsaXb6h63RZbRYPffXYd5Y91+\nfnrBQEb1TQh0c7wPDiLSBXgH+IUx5lhjVd2UmRaUNyw05iVjTKYxJrNHjx5NNVkp1ZSVj0B1veRt\n1RVWuR8FIstoMCuvrOFX72QzIKkz91yaEejmAF4GBxGJxAoMrxlj3rWLD9ldQti3RXZ5PtDXZfM+\nwIEmyvu4KVdK+VtZfvPKfSQQWUaD2TPLcykoreDp60YQHRnY7iQnb2YrCfAysM0Y86zLU4sB54yj\nWcD7LuU327OWxgNldrfTcmCqiCTaA9FTgeX2c8dFZLx9rJtd9qWU8heHw3OCva593Jf7iKdsov7M\nMhqs1u0p4Z9f7mXWhP6c298/q++1hDdnDhOBHwKTRWST/XMZ8BRwqYjsBC61HwN8COwGdgF/B34O\nYIwpAR4FvrZ/HrHLAG4H/mFv8y2wzAe/m1KqMZ//L9ScslJ0u4qMgSkP+vXQc6YNJqbeN2QB7rkk\n3a/HDTYVVbXct/Ab+naL4b7p/suw2hJNJt4zxnyB+3EBgClu6hvgDg/7mgfMc1OeBQxrqi1KKR/Z\ntRI+eRxGfA8GXWKNMZTlW2cMUx6EETf49fDObKLONZG7dY7iSHkVe4+c9Otxg81z/9nB3iMnef0n\n5xEbFVx5UIOrNUop/yvNg3d+AsnnwBXPQVRnvwcDd2aOTj0j5fSct7/hxU+/5dIhKYwMgtk6/rZx\n/1H+8fluvn9eGt8ZlBTo5jSg6TOUCiU1lbDgZnDUwPf+bQWGIPG7K4eQHNeJ/3n7G07Vm8nU0VTW\n1DJnYTYp8dE8MOPsQDfHLQ0OSoWSj+6HAxtg5ottkn21OeKjI/nDtSPYVXSC51bsCHRz/GLRxgIm\nPrWKwb/9iF1FJ7hyZC/ioiOb3jAANDgoFSo2vQFZ82Di3XDOFYFujVsXZvTgpnFpvPT5btbvK2l6\ng3Zk0cYCHng3hwKX6br/WrM/aC/80+CgVCg4mANLfgH9L4DJ/p2J1Fq/ufwceneN4ZdvZ1NR1XG6\nl9rbhX8aHJTq6CpK4a0fQkwiXDcPwoN7HkqXThE8c/0I9hwuD9oPzpZobxf+aXBQqiNzOGDR7VbO\npOvnQ5fkprcJAt85K4lZE/rxypd7WLv7SKCb4xPt7cI/DQ5KdWT/fQ5yP4Spj0PaeYFuTbP8asbZ\n9E2MZc7CbMorawLdnFZzd+FfTGQ4c6YF18VvThoclOqodq+GVY/BsOvgvJ8FujXNFhsVwR+vH0ne\n0ZP84aPtgW5Oq80cncqT1wwnNSEGAVITYnjymuFnXOsRTIK781Ep1TJlBbDwVkjKgCufb5N1of1h\n3IBu/HjiAF7+Yg/Th/YMyovFmqP+hX/BTM8clOpoaqrg7VlW3qQb/gWdugS6Ra0yZ9pgBiZ1Zs7C\nbI6fqg50c0KGBgelOprlv4b8r+Hqv0KP4FgboDWiI8P54w0jKSyr4IkP23/3Unuh3UpKtXfZC04n\nzotJhIoSmHAnDO04q+2OSUvkpxcO5P99upvpw3oyKUMX+/I3PXNQqj1zLvNZlgcYKzBIGKQMD3TL\nfO6eSzJIT+7CrxZmU1ah3Uv+psFBqfbM3TKfxgGfPBaY9vhRdGQ4f7x+JMUnKnlsydZAN6fD024l\npdqzAC3zGSgj+yZw+6Sz+Msnu1i1vYiS8ip6J8QwZ9rgdjMLqL3QMwel2itjIMbDugd+XuYzkPp3\nj0WAI+VVGKCgtIIH3s0J2gR27ZUGB6Xao6qTsOjnUHHUGmNw1QbLfAbSc//ZialXFswJ7NqrJoOD\niMwTkSIR2exS1k1EVojITvs20S4XEZkrIrtEJFtExrhsM8uuv1NEZrmUjxWRHHubuSLt9GodpdrK\nkW/h5Uvhmzdg0v3W2gxd+wJi3V45NyAru7WV9pbA7gzZC+C5YfBwgnWbvSDQLfLImzGHfwJ/AV51\nKbsfWGmMeUpE7rcf/wqYAaTbP+cBLwLniUg34CEgEzDAehFZbIw5ateZDXwFfAhMB5a1/lcLIq5T\nDdtojV7VQW1dDO/fAWHh8IOFkH6JVT7yxsC2qw31Tog5Y00E1/Kg5pxZ5pxAUJZnPYag/Dxo8szB\nGPMZUH/VjauB+fb9+cBMl/JXjeUrIEFEegHTgBXGmBI7IKwAptvPxRtj1hhjDFYA6jiTs6HhVEPn\nGyKIvzGoIFRbDct/Awt+CEnp8LPPTweGENPeEtjVcTezrLrCKg9CLR1zSDHGFALYt848wKlAnku9\nfLussfJ8N+UdRzt7Q6ggdKwQ5l8Ja/4C5/4UblkGCX0D3aqAaW8J7Op4nFmWB3u/sNKeBBFfT2V1\nN15gWlDufucis7G6oEhLS2tJ+9peiE01VD625zNY+GOoKodr/gEjrg90i4JCe0pgV6drH7sHwY1/\nXg6RsdBvIpx1MQy8GJLPOTNhYht3T7c0OBwSkV7GmEK7a6jILs8HXL/S9AEO2OUX1StfbZf3cVPf\nLWPMS8BLAJmZmR6DSNDYvtT64xo3Te3AUw2VDzgc1loMqx6D7oNg1hJIPjvQrVKtMeXBM8ccwJpZ\nNv0p6JwMuz+Bbz+xcmMBdEmBgRdZgaLyOPznwTYdr2hpcFgMzAKesm/fdym/U0TexBqQLrMDyHLg\nCeesJmAq8IAxpkREjovIeGAtcDPw5xa2KXiUH4Fl98HmhRDfB04WQ03l6ec7+FRD1QKu3wrje0Ns\ndziYDUOvgavmQqe4QLdQtZbzQ9zTt/+zL7Nuy/KttTi+/QR2rYTst9zvz9k97afgIMbdt1rXCiJv\nYH3rTwIOYc06WgQsANKA/cD19ge9YM1smg6cBG4xxmTZ+/kxYIdEHjfGvGKXZ2LNiIrBmqX0f01T\njcI6c8jKymrO79o2tr4PS//HWrf3wjlwwb2w5T2draQ8qz+LxWnEjfDdv7XbtRiUDzgcULQF/na+\nhwoCD5d6vTsRWW+MyfSqrhefw0Ep6ILDiWL48JewdRH0GglXvwA9hwW6Vao9eG6Y+77orn3hns0N\ny1Xo8dF7pDnBQa+Qbi1jYPO78MJ51lq9k38HP1mpgUF5TyctqKZMedDqjnbl5+5pTbzXGieKYOm9\nsO0D6D0GZr5gzTBQqjk8zWLRSQvKqanxCj/Q4OCt+tPI0qfBlnesHDeX/N5aXCVcX07VAp5mseik\nBeVqxA1tOlapn2becHfZe9Y/IHEA/Hg59AjyKzNVcAvAt0KlmqLBwRvurnIGcFRrYFC+0cbfCpVq\nigaHxtTWWBemeLqqsUzzxyulOiYNDu4U58Km1+Cbt+DEQStfvnE0rKcDhkqpDkqDg1PFUWtK6qbX\noGA9SDhkTINR37cuXV96rw4YKqVCRmgFh/ozjib/FmKTrICwfSnUVkLyEJj6uNX/2yX59LZhETpg\nqJQKGaFzhbSnFAUAMYkw/AbrLKHXSE1XoJTqkJpzhXTonDl4mnEU2x3u3QYRndq+TUopFaRCJ32G\np1QEJ0s0MCilVD2hExw8zSzSGUdKKdVA6ASHACSuUkqp9ip0gsOIG+DKuVaKW8S6vXKuzjhSSik3\nQjF5Z1IAAAbxSURBVGdAGjRFgVJKeSl0zhyUUkp5TYODUkqpBjQ4KKWUakCDg1JKqQY0OCillGqg\n3eZWEpFiYF8LN08CDvuwOb6i7WoebVfzaLuapyO2q58xpoc3FdttcGgNEcnyNvlUW9J2NY+2q3m0\nXc0T6u3SbiWllFINaHBQSinVQKgGh5cC3QAPtF3No+1qHm1X84R0u0JyzEEppVTjQvXMQSmlVGOM\nMe3uB5gHFAGbXcpGAmuAHOADIN7luQeAXUAuMM2lfLpdtgu438OxOgFv2XXWAv391SagL/AJsA3Y\nAtzt4VgXAWXAJvvnwTZ6vfba9TcBWR6OJcBce/tsYIw/2wUMdnkdNgHHgF+05jVrTruA7vbf7ATw\nl3r7GWvX32W/JtJWr5endgGxwFJgu/0ee8rDsfoDFS6v19/a4PVabf9tncdM9nA8t+9PP71ecfXe\nX4eBP7Xh63UpsN4uXw9M9tf7q8G23lYMph/gQmBMvRf3a2CSff/HwKP2/SHAN1gf8gOAb4Fw++db\nYCAQZdcZ4uZYP3f+oYEbgbf82KZezj+e/abc4aFNFwFL2vL1sp/bCyQ1cazLgGX2m3I8sNbf7XLZ\nNhw4iDWXu8WvWTPb1Rk4H7iNhh9264AJ9muxDJjRhq+X23ZhBYeL7ftRwOce2tXf9Tht9HqtBjKb\nOFaT7wNft6vePtcDF7bh6zUa6G3fHwYU+Ov9Vf+nXXYrGWM+A0rqFQ8GPrPvrwCute9fDbxpjKk0\nxuzBiqDj7J9dxpjdxpgq4E27bn1XA/Pt+wuBKSIi/miTMabQGLPB3t9xrDOIVM+vhHd89Hp562rg\nVWP5CkgQkV5t1K4pwLfGmJZeHNnsdhljyo0xXwCnXCvbv3O8MWaNsf5LXwVmujmcX14vT+0yxpw0\nxnxi368CNgCtWg7RF+1qBq/fn75ul4ikA8lYAbXFmtmujcaYA3b5FiBaRDr54/1VX7sMDh5sBq6y\n71+P1UUD1odrnku9fLvMU3l9dfWMMTVYXRPd/dSmOiLSH+tbw1oP+54gIt+IyDIRGeple1rbNgN8\nLCLrRWS2h/16+7r6sl1ONwJvNLLv1rxmntrlSardRqcm319N1PNVu+qISAJwJbDSQ5UBIrJRRD4V\nkQua0abWtOsVEdkkIr9z9yWMAL5ewE1YPQeeZvH4+/W6FthojKmkDd5fHSk4/Bi4Q0TWY3XJVNnl\n7t5gppHy+ryt54s2WU+KdAHeweo7P+am7gasrpORwJ+BRV62p7Vtm2iMGQPMsLe90E3d1rxeLW0X\nIhKF9c/1tof9tvY189QuT9ri/dWSdlkHFYnACqRzjTG73VQpBNKMMaOBe4HXRSTez+36gTFmOHCB\n/fNDd013U+b318vW2JcPv75e9peZPwA/cxa52YdP318dZiU4Y8x2YCqAiGQAl9tP5XNmFO4DOE/T\nPJW7cm6fb/9DdaXhKaHP2iQikViB4TVjzLse9nvM5f6HIvKCiCQZY7zOt9KStjlPb40xRSLyHtbp\n/GecqbHX2y/tss0A/n/75u8aRRTE8c9oBPEQGw0WNglExMZGJPgniKbRIjax0CKFoK3EKlhYiYKC\nCDYWtmLANp0gaALxB2L0xEIJFhYSLYLFWMycbm7vNiF3GyV8P7DsY9+P/d7ssLPvvbl5d//aZdye\nbFahqxufWb1cs5Z/rdWuX7pa3APeu/vNLuOuACtZnjOzJnAQeFGXLnf/kudlM3tI+NeDtmb/xF5m\ndgQYcPe5LuPWZi8zOwA8AibcvZmXa/evLTNzMLPBPG8DrgJ3s2oGGM91uiFghNjIeQ6MmNlQfnWO\nZ9t2ZoBzWT4DzFZMK3vSlNPo+8Bbd79RMe7+1pTbzI4Rz/HbejT1oK1hZruzT4Nw5Ncdhp4BJiwY\nBb67+1Jdugpdz1KxpNSrzSp0dSR/87KZjeZ9J4DHHZrWZa+qPteIj5zLFW32mdn2LA8T9u40w+iL\nLjMbMLO9Wd4BnKS7f1X5QV91FVjLv2qxVy79PQGuuPvTVvtN8a/17lz/TwfxkJaAX0RkPA9cIrJ7\nFoHrFNK6gCkiq+EdhR19Yid/MeumCtengbEs7ySWKj4QTjhclyYiW8KJlLNWStyJrJsEJrN8kdic\nWgCeAcfrtheR1bWQx5s2exW1GXAn+7+iIvukj89xF/Gi39M2/oZstgFdn4jZ5I9sfzivHyVecE3g\ndqvPJtqrpIv4cnQi2aHlYxey/RgwneXTBXvNA6dq1tUgMoFe5n1v8TdL7o+uKj+o6zlm3UfgUNv4\ntduLCBQ/WZ1OO1iHf7Uf+oe0EEKIEltmWUkIIUT/UHAQQghRQsFBCCFECQUHIYQQJRQchBBClFBw\nEEIIUULBQQghRAkFByGEECV+A0m6jdPg1anNAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "plt.plot(year, populations, 'o-') \n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.2.4.3 recarray：仅仅方便" ] }, { "cell_type": "code", "execution_count": 125, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([b'a', b'b'], \n", " dtype='|S1')" ] }, "execution_count": 125, "metadata": {}, "output_type": "execute_result" } ], "source": [ "arr = np.array([('a', 1), ('b', 2)], dtype=[('x', 'S1'), ('y', int)])\n", "arr2 = arr.view(np.recarray)\n", "arr2.x" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([1, 2])" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "arr2.y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.2.4.4 矩阵：方便？\n", "\n", "- 通常是2-D\n", "- \\* 是矩阵的积，不是元素级的积" ] }, { "cell_type": "code", "execution_count": 126, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "matrix([[1, 2],\n", " [3, 4]])" ] }, "execution_count": 126, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.matrix([[1, 0], [0, 1]]) * np.matrix([[1, 2], [3, 4]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2.2.5 总结\n", "\n", "- ndarray的剖析：data、dtype, 步长\n", "- 通用函数：元素级操作，如何常见一个新的通用函数\n", "- Ndarray子类\n", "- 整合其他工具的多种buffer接口\n", "- 最近的补充：PEP 3118，广义ufuncs\n", "\n", "## 2.2.6 为Numpy/Scipy做贡献\n", "\n", "看一下这篇教程：http://www.euroscipy.org/talk/882\n", "\n", "### 2.2.6.1 为什么\n", "\n", "- “这有个bug?”\n", "- “我不理解这个要做什么？”\n", "- “我有这个神器的代码。你要吗？”\n", "- “我需要帮助！我应该怎么办？”\n", "\n", "### 2.2.6.2 报告bugs\n", "\n", "- Bug跟踪（推荐这种方式）\n", " - http://projects.scipy.org/numpy\n", " - http://projects.scipy.org/scipy\n", " - 点击“注册”链接获得一个帐号\n", "\n", "- 邮件列表 ( scipy.org/Mailing_Lists )\n", " - 如果你不确定\n", " - 在一周左右还没有任何回复？去开一个bug ticket吧。\n", " \n", "#### 2.2.6.2.1 好的bug报告\n", "\n", "---\n", "\n", "Title: numpy.random.permutations fails for non-integer arguments\n", "\n", "I'm trying to generate random permutations, using numpy.random.permutations\n", "\n", "When calling numpy.random.permutation with non-integer arguments it fails with a cryptic error message::\n", "\n", " >>> np.random.permutation(12)\n", " array([ 6, 11, 4, 10, 2, 8, 1, 7, 9, 3, 0, 5])\n", " >>> np.random.permutation(12.)\n", " Traceback (most recent call last):\n", " File \"\", line 1, in \n", " File \"mtrand.pyx\", line 3311, in mtrand.RandomState.permutation\n", " File \"mtrand.pyx\", line 3254, in mtrand.RandomState.shuffle\n", " TypeError: len() of unsized object\n", "\n", "This also happens with long arguments, and so np.random.permutation(X.shape[0]) where X is an array fails on 64 bit windows (where shape is a tuple of longs).\n", "\n", "It would be great if it could cast to integer or at least raise a proper error for non-integer types.\n", "\n", "I'm using Numpy 1.4.1, built from the official tarball, on Windows 64 with Visual studio 2008, on Python.org 64-bit Python.\n", "\n", "---\n", "\n", "0. 你要做什么？\n", "1. **重现bug的小代码段**（如果可能的话）\n", " - 实际上发生了什么\n", " - 期望发生什么\n", "2. 平台（Windows / Linux / OSX, 32/64 bits, x86/PPC, ...）\n", "3. Numpy/Scipy的版本" ] }, { "cell_type": "code", "execution_count": 128, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.12.1\n" ] } ], "source": [ "print(np.__version__)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**检查下面的文件是你所期望的**" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "print(np.__file__)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "以免你想要旧/损坏的Numpy安装在哪里\n", "\n", "如果不确定，试着删除现有的Numpy安装文件，并且重新安装...\n", "\n", "### 2.2.6.3 为文档做贡献\n", "\n", "1. 文档编辑器\n", " - http://docs.scipy.org/numpy\n", " - 注册\n", " - 注册一个帐号\n", " - 订阅scipy-dev邮件列表（仅限订阅者）\n", " - 邮件列表有问题：你可以发邮件\n", " - 但是：**你可以关闭邮件送达**\n", " - 在http://mail.scipy.org/mailman/listinfo/scipy-dev 底部“改变你的订阅选项”\n", " - 给@`scipy-dev`邮件列表发一封邮件；要求激活：\n", "\n", " ---\n", " > To: scipy-dev@scipy.org\n", "\n", " > Hi,\n", "\n", " > I'd like to edit Numpy/Scipy docstrings. My account is XXXXX\n", "\n", " > Cheers,\n", " \n", " > N. N.\n", " \n", " ---\n", "\n", " - 检查一下风格指南：\n", " - http://docs.scipy.org/numpy/\n", " - 不要被吓住；要修补一个小事情，就修补它\n", " - 编辑 \n", "2. 编辑源码发送补丁（和bug一样）\n", "3. 向邮件列表抱怨\n", "\n", "\n", "### 2.2.6.4 贡献功能\n", "\n", "0. 在邮件列表上询问，如果不确定应该它应该在哪里\n", "1. 写一个补丁，在bug跟踪器上添加一个增强的ticket\n", "2. 或者，创建一个实现了这个功能的Git分支 + 添加增强ticket。\n", " - 特别是对于大的/扩散性的功能添加\n", " - http://projects.scipy.org/numpy/wiki/GitMirror\n", " - http://www.spheredev.org/wiki/Git_for_the_lazy" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# 克隆numpy仓库\n", "git clone --origin svn http://projects.scipy.org/git/numpy.git numpy\n", "cd numpy\n", "\n", "# 创建功能分支\n", "git checkout -b name-of-my-feature-branch svn/trunk\n", "\n", "\n", "\n", "git commit -a" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- 在http://github.com （或者其他地方）创建一个帐号\n", "- @ Github创建一个新仓库\n", "- 将你的工作推送到github" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "```shell\n", "git remote add github git@github:YOURUSERNAME/YOURREPOSITORYNAME.git\n", "git push github name-of-my-feature-branch\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.2.6.5 如何帮助，总体而言\n", "\n", "- 永远欢迎修补bug！\n", " - 什么让你最恼怒\n", " - 浏览一下跟踪器\n", "- 文档工作\n", " - API文档：改善文档字符串\n", " - 很好的了解了一些Scipy模块\n", " - 用户指南\n", " - 最终需要完成\n", " - 想要想一下？看一下目录\n", " http://scipy.org/Developer_Zone/UG_Toc\n", "- 在沟通渠道上询问：\n", " - `numpy-discussion` 列表\n", " - `scipy-dev` 列表" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.1" } }, "nbformat": 4, "nbformat_minor": 1 }