第三章 使用Python进行数字计算

尽管IPython强大的shell和扩展的控制台能被任何Python程序员使用,但是这个包最初是科学奖为科学家设计的。它的主要涉及目标就是为使用Python进行交互式科学计算提供一个方便的方式。

IPython只是为NumPy、Scipy、Pandas、Matplotlib等包提供一个交互式接口。其本身并不提供科学计算的功能。这些工具组合在一起就形成了可以匹敌如Matlab、Mathmatic这样被广泛使用的商业产品的科学计算框架。

NumPy提供了支持响亮操作的多维数组对象。SciPy基于NumPy提供了大量的科学计算算法(信号处理、最优化求解等等)。Pandas为处理来自真是世界数据集的表数据提供了方便使用的数据结构。Matplotlib允许你轻易的绘制图像以交互式可视化任何格式的数据,生成的图像甚至能达到出版水平。IPython为使用这些工具提供了绝佳的流式交互式框架。

在这一章节,我们将:

  • 探索NumPy和Pandas的交互式计算;
  • 理解为什么多维数组能很好的是英语高性能计算;
  • 认识怎样将数组应用于具体的应用;
  • 寻找包含更多高级示例与应用的参考资料。

向量计算基础

在这一章节我们将会介绍向量计算这一概念。这个概念相当重要,因为它是Python这样的高级语言获得高性能的最简单的方法。

Python循环计算的简单示例

现今的科学与工程都是关于数字的。大多数的数据处理和数字仿真就是大量数字数据的元素操作的成功而电脑十分擅长做这方面的事。然而数据必须要适当合理的被组织起来。数字数据的组织结构一般是向量和矩阵,更普遍的是多维数组。

在我们解释关于数字数组的更多细节之前,让我们来看一个简单的示例来对这些对象有一个直观的认识。现在假设我们获得了大量地点带有坐标(经度,维度)的地理数据,我们要在这里面寻找出距离我们感兴趣的地方最近的地点。例如,我们也许想要找到距离一个智能手机用户的GPS定位地点最近的餐馆。

如果这个地点存储在一个Python元祖列表中,我们可以写出如下的代码:

def closest(position, positions):
    x0, y0 = position
    dbest, ibest = None, None
    for i, (x, y) in enumerate(positions):
        # squared Euclidean distance from every position to the position of interest
        d = (x - x0) ** 2 + (y - y0) ** 2
        if dbest is None or d < dbest:
            dbest, ibest = d, i
    return ibest

在这里我们遍历寻找了每一个地点。变量i保存当前位置(x,y)的索引。我们感兴趣的地点是position(x0,y0)。在首次迭代中,当前位置一直保存最佳目标,在接下来的迭代中如果有比当前坐标更近的就将其更新为当前坐标。循环结束后,最近地点的索引是ibest,也即我们寻找到的最近地点是position[ibest]。我们感兴趣的地点与选找到的最近地点之间的平方距离保存在dbest。我们可以使用欧几里得距离公式来计算这个距离:$D=(X-x_0)^2+(Y-y_0)^2$

这是一个基础的标准算法。让我们在一个大数据集上来评估一下这个算法的性能。我们首先像下面这样随机产生1000W的随机地点:

In[1]:import random
In[2]:positions = [(random.random(),random.random()) for _ in xrange(10000000)]

我们定义了一个有1000W个0到1之间的随机数字对组成的数组——position。现在,让我们使用下面的命令来评测一下:

In[3]:%timeit closest((.5,.5),position)
1 loops, best of 3: 16.4 s per loop

这个算法处理1000W个地点花费了16.4秒。让我们来看一下这是否接近一个CPU的理论最大值。 代码运行在2GHz的单核处理器上。每一个时钟周期理论上能够处理四个浮点操作,也就是每秒能进行80亿次操作。在我们的算法中,每一次迭代包含5次数学运算和一次对比操作,所以共有500W次操作(只计算数学运算)。理论上的最佳性能应该是6.25ms。这意味着我们的算法比理论值糟糕2600倍!

当然,这是一次非常幼稚的评测,且理论值在实际中很难达到。但是2600倍的差距还是非常蛋疼糟糕的。我们能做的更好吗?我们将会在下一章节进行讲述。

什么是数组

在先前的例子中,相同的计算两点距离的运算需要使用在大量的数字上。Numpy提供了一个完美适应这种情况的新数据类型:多维数组。那么,数组是什么呢?

数组是将数据组织成若干个维度的数据块。一维数组是一个向量,它是有一个整数进行索引的有序元素序列(元素通常是数字)。一个二维数组是一个矩阵,包含用一对整数索引的元素,整数对即行索引和列索引。更一般地,一个n维数组是一组由n个整数的元组进行索引的具有相同数据类型的元素集合。

【图】

数组中的所有元素必须是相同类型的:这被称作数据类型(dtype)。NumPy中可用的类型可以是:布尔值、有/无符号整数、单/双精度浮点数、负数、字符串等等。也可以自定义数据类型。

一个数组中的元素被存储在内存的一块连续区域。例如大小为10的向量会占据10个连续的内存地址。当数组的维度大于等于2时,组合数组元素的方式并不是唯一的。对于一个矩阵,根据索引值增长的快慢可以按行主序存储,也可以按照列主序进行存储。这一概念在三维数组以至更高维度的数组都是通用的。NumPy中默认的存储方式是行主序,但是在创建数组的时候可以使用key关键字参数对其进行改变。

【图】

这一概念可以扩展到任意维度。stride定义了贯穿每个维度中所有元素的步数。NumPy会自动处理所有这些底层细节,为创建,操作和计算这些数组的提供了十分方便的方法。大部分情况下,我们只需把我们的变量当作多维数组,而不必理会这些底层细节。然而,知道其内部是怎么工作的话, 我们就能比较清楚怎样修复错误,在涉及非常大的数组时明白怎样来优化代码。

相较于原生的Python类型,数组的优势是可以以向量运算替换Python的循环操作来进行非常搞笑的计算。 他们之间的不同是NumPy中的循环是用C代替Python进行实现的,这样就节省了循环中Python代码解释的时间。说实话,Python作为一种可交互的动态类型语言,其每一次的迭代操作包含着大量的低级操作(类型检查等等)。这些操作花费的时间通常可以忽略不计,但是当他们重复几十万次的时候就会对性能造成很大的影响。

除此之外,现代的CPU大部分通过能够包含几个单精度或双精度的大寄存器(128位或256位)实现了矢量指令(SSE、AVX、XOP等等)。如果NumPy采用适当的编译选项,数组操作可以直接使用这些矢量CPU指令从而使运算速度有四倍还多的提升。

这些就是为什么NumPy的矢量运算比Python自身的循环操作高效的多的主要原因。由于NumPy数组操作中相同的计算操作应用到了多个元素上面,从而参照实现了单指令、多数据(SIMD)计算范式。我们将在先前示例的帮助下来证明这一点。

使用数组重新实现示例

让我们使用数组来重写我们的示例。首先,我们需要导入NumPy。在IPython中我们可以使用%pylab魔法命令(或以ipython --pylab形式启动IPython)将NumPy和Matplotlib加载到交互命名空间(numpy被简写作np,Matplotlib.pyplot被简写做plt)。这是在IPython交互式会话中使用BumPy的最方便的方法。也可以使用import numpyfrom numpy import *来导入NumPy(当然也可以import numpy as np做一下输入简化)。前者的语法是脚本编写中的首选,而后者则更适合交互式会话。在这里以及下面所有的章节中,我们将一直假设pylab模式已经像下面这样被激活:

In[1]:%pylab

首先,我们需要生成一些随机数据。NumPy为这样的操作提供了如下的高效方法:

In[2]:position=rand(10000000,2)

position数组是100W行100W列包含0到1之间独立随机数字的二维数组。我们注意到我们这次创建数组没有使用循环。能使用NumPy操作就尽量避免循环操作。让我们看一下这个对象的属性:

In[3]:type(position)
Out[3]:numpy.ndarray
In[4]:position.ndim,position.shape
Out[4]:2,(10000000,2)

shape属性以一种整数元组的形式包含着数组的形状。一个数组包含的其他重要属性:

  • ndim:维度数,等于len(position.shape)
  • size:元素的总数(position.shape中所有值)
  • itemsize:以bytes为单位的每个元素的大小(int32数据类型是4,float64是8,等等)

现在我们可以在两部之内计算出两点之间的距离。首先,我们输入下面的命令:

In[5]:x,y = position[:,0],position[:,1]

这里xy包含着所有点的xy坐标。事实上,变量position[:,0]是指position第一列(Python中索引从0开始)。这是Python/NumPy中一种特殊的索引语法。中括号[]允许从Python容器对象中访问元素。中括号中的冒号:,0是指数组中的第一列的所有元素。因此,在NumPy中第一个维度始终引用行,第二位始终引用列,在这里我们特指第一列。相似的,position[:,1]引用第二列中的所有数据。变量xy都是一维向量。让我们使用下面的命令计算变量distances的值:

In [6]: distances = (x - .5) ** 2 + (y - .5) ** 2

这里我们计算地点(0.5,0.5)距离其他所有地点的距离向量。事实上,x-.5这个表达式是将position的第一列元素全都减去0.5。原因是x是一个有100W个元素的一维向量,0.5只是一个浮点数。NumPy中约定采用矢量微积分的数学约定,也就是说减法会应用于数组中的所有元素。

同样,(x-.5)**2计算出现咋括号中的向量中的所有元素的平方。最后,+运算符对两个100W长的向量进行主调操作。

我们可以看到NumPy允许以一种十分简单的语法来进行响亮操作。使用数组进行计算是一种非常特殊的编程方式,掌握她需要一些时间。在大多数语言中这是完全不同于编程的标准顺序方式,但是,正如下面所示,在Python中这是相当高效的:

In [7]: %timeit exec(In[6])
1 loops, best of 3: 508 ms per loop

但我们再次使用%timeit魔法函数来评测计算变量distances的效率时,我们可以发现计算效率明显的比纯Python版本的高。即使我们添加上最小元素的计算,我们依然发现比纯Python版本的快30多倍:

In [8]: %timeit ibest = distances.argmin()
1 loops, best of 3: 20 ms per loop

故而,在我们进行大数量的数字计算的时候应尽量避免使用Python循环。一个算法进行向量化计算也许有事比较困难,但通常对于性能提升是值得的。

数组的创建和加载

在这一章节,我们将会看到怎样创建和加载数组,以及怎样根据现有数据进行创建和加载。这往往是使用Python进行数据分析的第一步。

创建数组

创建数组有多种方式,在这一章节我们将对其一一查看。

逐个元素

首先,我们可以通过采用手动指定对应系数的方式创建数组。这是创建数组的最直接方式,但是在实际中很少使用这种办法。NumPyarray函数接受一个元素列并返回相应的NumPy数组,就像下面这样(IPython的pylab模式需要被激活):

In [1]: x = array([1, 2, 3])
In [2]: x.shape
Out[2]: (3,)
In [3]: x.dtype
Out[3]: dtype('int32')

这里我们创建了包含3个32位整数的一维数组(也就是向量,32位系统中整数的默认类型是32位)。创建的数组的数据类型根据array函数接受的元素自动指定的。我们也可以使用dtype关键字参数强制指定数据类型,就像下面这样:

In [4]: x = array([1, 2, 3], dtype=float64)
In [5]: x.dtype
Out[5]: dtype('float64')

要创建二维数组(矩阵),我们需要提供一个列表的嵌套列表,里面的每一个列表包含一行数据,如下:

In [6]: array([[1, 2, 3], [4, 5, 6]])
Out[6]:
array([[1, 2, 3],
[4, 5, 6]])

要创建一个n维数组,我们需要n级递归的嵌套列表。让我们使用两个嵌套列表解析创建一个乘法表:

def mul1(n):
    return array([[(i + 1) * (j + 1) for i in xrange(n)] for j in xrange(n)])

这个函数接受表的大小作为参数,根据行列表创建乘法表,如下:

In [7]: mul1(4)
Out[7]:
array([[ 1, 2, 3, 4],
[ 2, 4, 6, 8],
[ 3, 6, 9, 12],
[ 4, 8, 12, 16]])
In [8]: %timeit mul1(100)
100 loops, best of 3: 5.14 ms per loop

我们将会在后面砍到更多创建乘法表的高效方法。

使用预定义的模板 使用给每个元素手动指定系数的方法创建数组是非常不实用的。有一种比较方便的方法是可以使用NumPy中预定义的几个函数,根据给定的数组形状(shape)进行数组的创建。比如,要创建有100个0填充的向量,我们可以使用下面的命令:

In[1]:x = zeros(100)

要创建2维矩阵,我们需要提供一个包含两个整数元素的元组作为参数:

In[2]:x= zeros((10,10))

数据类型缺省状态下是float64.与此相似,one函数用来创建由1填充的数组。identityeyediag函数用以创建对角线矩阵。

也有一些比较方便的函数可创建元素均匀分布的向量,如:

In[5]:arrange(2,10,2)
Out[5]:
array([2, 4, 6, 8])

这里,我们创建了一个元素从2到10步长为2线性分布的向量。请记住第一个数字包含在内(第一个2),但队列中的最后一个数字,10,却不在其中。这是Python的一种通用约定。另外一个与此相关的函数是linspace,此函数与arange相比,除了以输出向量的大小而不是步长作为第三个参数,其它都很相似。这一次,序列中的第一个和最后一个元素都包含于内。

函数声明:包含参数序列和关键字列表的函数声明可以在IPython中使用?help()获得。除此之外,在Qt控制台会notebook中,输入linspace将会自动给出函数提示(可以使用Tab展开提示)

根据随机数创建数组 NumPy提供了能够产生不同独立分布的随机函数。例如,我们可以像下面这样使用rand函数创建由0到1之间的随机浮点数填充的2x5数组:

In [1]: rand(2, 5)
Out[1]:
array([[ 0.925, 0.849, 0.858, 0.269, 0.644],
[ 0.796, 0.001, 0.183, 0.397, 0.788]])

注意一下rand函数的双括号内用以指定数组形状的两个数两边是没有元元括号的,这一点和array函数以一个元组为参数不同。

IPython中的数字格式化:IPython中数字的显示形式可以通过%precision魔法命令进行指定。比如,要显示浮点数小数点后三位,我们可以在IPython输入%precision 3。实际上,任何格式都可以被指定,更多信息请使用%presicion?查阅其文档

其他随机函数还有randn(生成高斯分布随机数),randint(随机整数),exponential(指数分布)等等。与此相关的函数包括shufflepermutation用于对已存在的数组进行随机重排列。

数组加载

数组结构最有趣的一点就是可以从Python或其他外部资源中加载已经存在的数据。MumPy针对从文本、二进制缓冲或文件中加载多维数组提供了方便而高效的方法(文本可以是Python字符串、text/CSV文件)。除此之外,Pandas包在加载包含多种数据类型的表数据时非常有用。

从本地Python对象加载数据

我们时常需要将本地Python对象数据转换为NumPy数组。标准的方法是使用array函数。当我们直接指定他们的值来创建数组的时候,实际上我们是将Python数字列表转换成了数组。

从缓冲区或外部文件加载

另一个会经常遇到的情况是需要从内存缓冲区或一个外部文件加载二进制或字符串元素数据来创建数组。从一个我们知道其确定类型的Python对象加载数据,我们可以使用frombuffer函数获得一个NumPy数组。相似地,fromstring函数接受由界定符或任意数据类型的二进制数据分隔的ASICII文本。示例如下:

In [1]: np.fromstring('1 2 5 10', dtype=int, sep=' ')
Out[1]: array([ 1, 2, 5, 10])

fromfileloadtxtgenfromtxt等函数允许我们从文本或二进制文件加载数据并将其转换成NumPy数组。loadtxtgenfromtxt的简化版本,当文件格式比较简单时比较有用。fromfile函数加载二进制数据时非常高效。例如,我们可以使用下面的命令加载Facebook数据集中文本文件中的数据:

In [1]: cd fbdata
In [2]: loadtxt('0.edges')
Out[2]:
array([[ 236., 186.],
...,
[ 291., 339.]])

最后,将数组保存都文件和加载NumPy数组一样简单。这里基本有两个函数——savesavetxt,用来将数组保存到文本或二进制文件。与此相关,可以使用loadzsavez函数很方便地保存任意类型变量的字典(包括NumPy数组类型)。所有的这些函数都使用平台无关的文件格式。

使用Pandas

Pandas是另外一个比较新的用以方便高效加载和处理多种数据源的数据集的Python包。尤其是当处理与纯数字数据(数字矩阵或数组)相对的表数据集时相当有用。他可以处理缺失值和数据队列问题(例如,时间序列)。NumPy能使用被加载的数据集进行高效的计算。简单来说,Pandas提供了处理表数据的高级方法,NumPy提供了处理原始多维数组的低级方法。(译者注:这里的高级如同编程语言的高级低级,不是优劣之分,只是处理数据的层次不同)

'NumPy的未来发展:'NumPy的创造者Travis Oliphant现在正进行着NumPy的后续开发与推广。这个项目将会把NumPy、Pandas、SciPy、Numba、Theano等包所提供的功能统一到一个单一的框架中。

这里有一个示例来掩饰怎样使用Pandas加载数据。我们将下载和分析一个包含这世界上大量城市及其人口的数据集。该数据集由MaxMind创建,可以从http://www.maxmind.com免费下载。

'在线公开数据集':随着开放数据集的流行,越来越多的数据可以公开得到。使用从本书中学到的工具来分析有趣的数据是获得实践经验的一个良好途径。然而,在线寻找优质的数据资源往往不是意见容易的事。下面给出的链接网站上有不少高质量的数据集,这些数据主要由政府部门、国际组织、大学以及研究机构所维护。

我们先像下面这样下载ZIP文件并将其解压到一个文件夹(ZIP文件大约有40MB,所以下载它可能会需要一些时间):

In [1]: import urllib2, zipfile
In [2]: url = 'http://ipython.rossant.net/'
In [3]: filename = 'cities.zip'
In [4]: downloaded = urllib2.urlopen(url + filename)
In [5]: folder = 'data'
In [6]: mkdir $folder
In [7]: with open(filename, 'wb') as f:
            f.write(downloaded.read())
In [8]: with zipfile.ZipFile(filename) as zip:
            zip.extractall(folder)

为了方便起见,我们使用%bookmark citiesdata data命令为新创建的文件夹创建一个alias。现在我们就使用Pandas加载已经被解压的CSV文件。Pandas的read_csv函数能够打开任意的CSV文件,就像下面这样:

In[9]:import pandas as pd
In[10]:filename = 'data/worldcitiespop.txt'
In[11]:data = pd.read_csv(filename)

现在让我们来探索一下这个新创建的data对象:

In[12]:type[data]
In[13]:pandas.core.frame.DataFrame

data是一个DataFrame对象,该Pandas对象类型包含一个二维标记型数据结构,每行的数据类型可能不同(就像一个Excel工作簿)。就像一个NumPy数组,shape属性返回表的形状。但是和NumPy数组不同的是,DataFrame对象包含更丰富的结构,尤其是keys方法能够返回不同数据列的名字。示例如下:

In[13]:data.shape,data.keys()
Out[13]:((3173958, 7),
Index([Country, City, AccentCity, Region, Population, Latitude, 
Longitude], dtype=object))

我们可以看到data包含30W多行数据,每行有七列,分别表示每个城市所对应的国家、城市、人口和每个城市的地理坐标等数据。headtail方法允许我们大致浏览一下data中的首位部分的数据。请记住,当在IPython notebook中使用Pandas时,要显示的数据可以像下面示例这样格式化为HTML表格以方便阅读:

In[14]:data.tail()

下面是示例在notebook中使用的HTML表格形式: 【图】

我们可以从看到有些城市的人口是一个NAN(不是一个数字)。这是因为数据集中部分城市的人口值是缺失的,而Pandas将这些确实值以NAN进行标记。

我们将会在下一章节看到我们能够对数据集进行怎样的操作与计算以获得有用的信息。

使用数组

一旦NumPy数组被创建或被加载,我们就可以进行下面的这三种基础操作:

  • 选取
  • 处理
  • 计算

选取

选取就是访问和使用数组中的一个或多个元素。这一操作使用NumPy和Pandas都能完成

使用Pandas

让我们使用示例中Pandas打开的数据继续我们的实验。DataFrame对象data的每一列都能通过它的名字进行访问。在IPython中,tab补全能提示数据的不同列。在接下来的示例中,我们获得了所有城市的名字(·AccentCity是城市的全称,包含大写字母和重音符号):

In [15]: data.AccentCity
Out[15]:
0 Aixas
1 Aixirivali
...
3173956 Zuzumba
3173957 Zvishavane
Name: AccentCity, Length: 3173958

这一列是Series类的一个实例。我么可以通过索引访问确定的行。在下面的示例中,我们将获得第30001个城市的名字(牢记索引是从0开始的):

In [16]: data.AccentCity[30000]
Out[16]: 'Howasiyan'

所以我们可以使用索引访问一个元素。但是我们怎样才能根据城市的名字获得该城市的信息呢?例如,我们想要获得纽约的人口和GPS坐标。一个可能是循环遍历所有的城市并检查它们的名字,但是这样是非常慢的,因为Python循环即使忘得元素却没有充分利用他们。Pandas和NumPy提供了一个更为优雅和高效的方式——布尔索引

典型的,我们会在同一行代码中进行两步操作。首先,我们使用布尔值创建一个数组,针对每一个元素,布尔值能够标识它是否满足条件(在这里就是城市名字是否是纽约)。接着,我们将这个布尔数组作为索引传递给院士数组。结果就是整个数组中满足元素匹配条件为True的部分。示例如下:

In [17]: data[data.AccentCity=='New York']
Out[17]:
Country City AccentCity Region Population Latitude 
Longitude
998166 gb new york New York H7 NaN 53.083333 
-0.150000
...
2990572 us new york New York NY 8107916 40.714167 
-74.006389

相同的语法在NumPy和Pandas中都适用。结果我们发现有十几个城市名字叫纽约,但是只有一个在纽约州。为了利用Pandas对单个元素进行访问,我们可以使用.ix属性(ix用于索引),示例如下:

In [18]: ny = 2990572
In [19]: data.ix[ny]
Out[19]:
Country us
City new york
AccentCity New York
Region NY
Population 8107916
Latitude 40.71417
Longitude -74.00639
Name: 2990572

使用NumPy

现在,我们将这个序列对象转换成一个纯净的NumPy数组。我们将要从Pandas的世界转到NumPy的世界(记住Pandas是基于NumPy构建的的)。我们大部分操作只针对所有城市的人口数目,就像下面这样:

In [20]: population = array(data.Population)
In [21]: population.shape
Out[21]: (3173958,)

population数组是表示所有城市人口数的一维向量(如果人口值丢失就是NAN)。NumPy中可以使用基本的索引操作来访问纽约的人口。如下:

In [22]: population[ny]
Out[22]: 8107916.0

我来找出有多少没有丢失人口数据的城市。要做到这一点,我们将从population数组中选取所有人口数与NAN不同的元素。我们可以像下面这样使用NumPy的isnan函数:

In [23]: isnan(population)
Out[23]: array([ True, True, True, ..., True, True, False], 
dtype=bool)
In [24]: x = population[~_]
In [25]: len(x), len(x) / float(len(population))
Out[25]: (47980, 0.015)

注意~_包含isnan(population)的负值。我们发现在数据集中大约占1.5%的48000个城市有确定的人口数据。

更多的索引功能

更普遍的来说,索引允许我们使用数组中的任意数据。我们在先前的章节中已经看到怎样使用布尔条件来过滤一个数组。我们也可以直接指定我们想要使用的索引列表。如,如果x是一个一维NumPy数组,x[i:j:k]则表示已步长为k选取从第i行到第j行的数据(其中包含i不包含j)。如果i缺省,则会将其指派为0.如果j缺省,将会把数组在这一维的长度指派给他。如果使用负值则表示我们从后向前索引。最后,k的默认值是1.这些概念在多维数组中也是适用的。如,M[i:j:1]创建了二维数组M的一个子矩阵。同时我们使用x[::-1]来获得x向量的反转向量(逆序)。

像包含i不包含j这一年过的惯例使得我们处理数组的连续部分的时候非常方便。例如,假设x有2n个元素,他的前一半和后一半就是x[:n]x[n:]。除此之外,x[i:j]的长度直接就是j-i。最后,一般情况下不应让+1或-1这样的值出现在索引中。

关于数组视图非常重要的一点就是他们都指向内存的相同地址。所以对一个大数组进行索引视图不会进行内存的分配。改变视图中的值的同时也会改变原始数组的值。示例如下:


In [1]: x = rand(5)
In [2]: x
Out[2]: array([ 0.5 , 0.633, 0.158, 0.862, 0.35 ])
In [3]: y = x[::2]
In [4]: y
Out[4]: array([ 0.5 , 0.158, 0.35 ])
In [5]: y[0] = 1
In [6]: x
Out[6]: array([ 1. , 0.633, 0.158, 0.862, 0.35 ])

在这个示例中,y包含x中的所有元素甚至索引(在这里,索引是0、2、4)。改变y[0]的值时,同时也改变了x[0]的值。因为y[0]引用了x的第一个元素。如果我们不像只是引用原来的数组,我们可以使用y=x.copy()y=array(x)来强制创建一个数组。在后续部分,我们甚至可以使用dtype关键字参数来改变x的数据类型。

最后,另外一个进行数组选取的办法是直接传递一个明确的整数值数组进行索引。这被称作趣味索引(fancy indexing)。如果x是一维向量,indices是另外一个由正整数组成的一维向量(或列表),那么x[indices]将会返回一个包含x[indices[0]]``x[indicies[1]``···这样的一个向量。因此,x[indices]的长度等于indices的长度而不等于x的长度。就像下面这样:


In [7]: ind = [0, 1, 0, 2]
In [8]: x[ind]
Out[8]: array([ 1. , 0.633, 1. , 0.158])

记住,一个给定的索引可以在索引数组中出现多次。(译者注:比如其中的数组ind中0就出现了两次)

操作处理

数组可以被操作处理和整形,这有时使得在进行向量化运算的时候非常有用。通过操作处理从一个原始数组构建一个新的数组也是可能的。完整的方法列表可以在NumPy的参考文档中查阅: http://docs.scipy.org/doc/numpy/reference/ routines.html

整形

首先,如果数组中的元素是连续的,可以使用reshape方法对该数组的形状进行改变。就像下面这样:

In [1]: rand(6)
Out[1]: array([ 0.872, 0.257, 0.083, 0.788, 0.931, 0.232])
In [2]: x.reshape((2, 3))
array([[ 0.872, 0.257, 0.083],
[ 0.788, 0.931, 0.232]])

对于大多数的一维数组,reshape参数可以用-1指示它的值能自动推测出来。如x.reshape((2,-1))等价于`x.reshap((2,3))。

维度数可以使用ravelsqueeze、和expand_dims来改变,(ravel移除一个数组中的多为结构并返回一个;sueeze从一个数组的形状中移除一维内容;expand_dims在一个数组中插入一个新的维度)

重复和连接

titlerepeat函数可以用来创建一个数组的拷贝。或根据指定的轴来连接他相同的拷贝或多次复制每一系数。


In [1]: x = arange(3)
In [2]: tile(x, (2, 1))
Out[2]:
array([[0, 1, 2],
[0, 1, 2]])
In [3]: repeat(x, 2)
Out[3]:
array([0, 0, 1, 1, 2, 2])

在这里,我们首先使用x的两个完全相同的拷贝组成的垂直栈来创建一个新的数组,接着我们通过重复三次循环使用x的每一个元素创建了一个新数组。repeat的第二个参数也可以是一个表示系数x[i]被重复的reps[i]次的reps列表。

例如,让我们使用reshapetitle来创建一个多维表。思路是首先定义一个行向量和一个列向量,每个向量的元素都是1到n的整数排列。将这两个向量相称,如下:

def mul2(n):
M = arange(1, n + 1).reshape((-1, 1))
M = tile(M, (1, n))
N = arange(1, n + 1).reshape((1, -1))
N = tile(N, (n, 1))
return M * N

让我们用下面的命令来评测一下这个函数的运行时间


In [1]: %timeit mul2(100)
10000 loops, best of 3: 188 us per loop

这个函数的运行速度大队是先前使用Python循环的mul1的27倍。

而且我们可以使用hastackvstackdstack、或concatenate将几个数组逐维连接到一个单一数组。

与此相似,hsplitvsplitdsplit、 或split函数可以逐维将一个数组分解为几个数组。示例如下:

In [1]: x = arange(6)
In [2]: split(x, 2)
Out[2]:
[array([0, 1, 2]), array([3, 4, 5])]
In [3]: split(x, [2,5])
Out[3]:
[array([0, 1]), array([2, 3, 4]), array([5])]

split的第二个参数可以是一个整数,比如n,用来将数组分成n等分的数组,或者一个标识数组的那部分应该被分解的索引列表(也就是除第一个子数组外,其他每个子数组中的第一个元素的索引)。

广播

在先前的乘法表的实例中,我们不得不重复使用行和列的相同副本,这样我们才可以使两个数组具有相同的形状(n,n)。实际上,repeat这一步不是必要的,因为不同形状的数组在一定的条件下也能结合在一起;而这被称为“广播”。一般的规则是两维数组在他们的形状是相同的或一方为1时可以结合在一起。例如,形状分别是(1,n)和(n,1)的M和N两个数组就可以相乘结合在一起。因为在第一维M数组的形状是1,在第二维度N数组的形状为1.形状为1的维度会自动转换以匹配另一个维度,而且这个操作不牵涉内存拷贝。

因此,我们可以在乘法表实例中丢弃tile操作:

def mul3(n):
    M = arange(1, n + 1).reshape((-1, 1))
    N = arange(1, n + 1).reshape((1, -1))
    return M * N

接着来评测一下这个函数的运行时间:

In [1]: timeit mul3(100)
10000 loops, best of 3: 71.8 us per loop

结果就是mul3的速度大约是mul2的2.6倍,mul1的70倍。原因是tile操作牵涉到了数组复制和内存分配,而mul3只是相乘而已。

交换

有很多函数可以用来整轴交换某一数组的数据。例如,transpose函数用来交换一个数组的维度。可以使用axes关键字参数给出表示具体的交换索引。

其它的转换函数,如fliprflipud,可以用来从左到右或从上到下反转一个数组。roll可以根据指定的轴 进行元素的循环交换,rot90可以将一个数组顺时针选钻90度。

计算

创建和处理素组的最终关键是可以使数组进行高效的向量化计算。在具有合适形状的条件下,数组有四种基础运算。除此之外,针对NumPy数组在进行向量形式的处理中有大量的数学函数可以使用。

假如A和B两个NumPy数组具有匹配的形状,A+B,A-B,A X B, A/B是四个基础操作。需要注意的是,当A和B是二维矩阵时,A X B不是矩阵乘积。矩阵乘积运算是由dot函数提供的。是的,dot函数更普遍地应用是计算两个数组的点乘。

一般的单元运算包含-A,A**X(A的X次方),abs(A)(绝对值),sign(A)(符号函数,结果只有-1,0,1,结果取决于每一个元素的符号),floor(A)(每一个元素的基底,小于等于自身的最大整数),sqrt(A)(开方),log(A)(求自然对数),exp(A)(指数),和其他很多的数学函数(三角函数、双曲线函数以及一些算数运算函数等等)。

NumPy还提供了一下儿函数用于计算一个函数中所有或指定维度的元素的加和(sum)和乘积(prod)。axis关键字参数指定对哪一维度的元素进行家和运算。这个函数返回一个一维数组。

maxmin函数返回一个数组或一个给定维度中的最大值和最小值。argminargmax函数返回数组元素的最小值和最大值的索引。例如,继续我们先前的cities的例子,我们可以像下面这样写出locate(x,y)函数:


In [26]: def locate(x, y):
    # locations is a Ncities x 2 array with the cities positions
    locations = data[['Latitude','Longitude']].as_matrix()
    d = locations - array([x, y])
    # squared distances from every city to the position (x,y)
    distances = d[:,0] ** 2 + d[:,1] ** 2
    # closest in the index of the city achieving the minimum 
    distance to the position (x,y)
    closest = distances.argmin()
    # we return the name of that city
    return data.AccentCity[closest]
    In [27]: print(locate(48.861, 2.3358))
Paris

locate函数接受一个地点的经纬坐标值做参数,并返回最近的城市的名字。argmin函数返回距离指定地点最近的城市的索引。

最后,meanmedianstd、和var这样的函数用来计算给定维度或整个数组中的元素的中位数、标准差、方差。而且Pandas 对象的describe方法可以提供一些有用的统计结果(包含平均数、标准差、50%分位数/中位数、25%分位数和75%分位数)。示例如下:


In [28]: population.describe()
count 47980.000000
mean 47719.570634
std 302888.715626
min 7.000000
25% 3732.000000
50% 10779.000000
75% 27990.500000
max 31480498.000000

与数学模型模拟相关的有用函数包括diff(离散区别)、cumsun(累积和)和cumprod(累积积)。diff函数可以用来计算信号的离散导数(达到一个标量数目)。cumsum可以用来计算信号的不定积分。

高级数学处理

NumPy提供了使用Python进行高效数值计算的所有必须的类型与惯例。基于NumPy的SciPy实现了大量的高级数学处理算法。这些算法可以延伸到数值计算的多个领域,例如最优化、线性代数、信号处理、统计学、等等。而且大量的科学工具包(SciKits)包(scikit-learn、scikit-image等)和更高级的包实现了许多专门领域的特有算法(机器学习,图像处理等)。

我们在这里给出一个关于SciPy和一些其他包提供的科学计算功能的概要列表。完整的功能列表可以在官方参考指南中找到: http://docs.scipy.org/doc/scipy/reference/。给出实际实例应用已经超出了本书的范围,感兴趣的读者可以在NumPy Pack出版社出版的Ivan Idris的NumPy和Francico Blanco-Silva的Learning SciPy for Numerical and Scientific Computing中找到大量的实例。

  • 线性代数:
  • 最优化求解
  • 数值积分
  • 信号处理
  • 统计

总结

在这一章中我们讲述了NumPy提供的多维数组,并展示了如何使用它对数值数据集进行高效计算。尤其是它能很好的适用于加载任意类型的数据,Pandas使得这一任务变得跟简单,即使面对非常复杂的数据文件也能胜任。在NumPy、SciPy这样的数据包和SciKit库使得我们在IPython中使用高级算法成为可能。然而,这个主题已经超出了本书的范围,感兴趣的读者可以在NumPy Pack出版社出版的Ivan Idris的NumPy和Francico Blanco-Silva的Learning SciPy for Numerical and Scientific Computing中找到大量的实例。

在下一章我们展示IPython和Matplotlib提供的与可视化相关的功能。我们经常将他们和NumPy结合在一起来进行数据的交互式可视化。

end