Python-数据分析教程-四-

52 阅读24分钟

Python 数据分析教程(四)

原文:Python Data Analysis

协议:CC BY-NC-SA 4.0

十一、Python 生态系统和云计算之外的环境

在 Python 生态系统之外,诸如 R、C、Java 和 Fortran 等编程语言相当流行。在本章中,我们将深入研究与这些环境交换信息的细节。

云计算旨在通过互联网提供计算能力。这意味着我们不需要在本地拥有很多强大的硬件。取而代之的是,我们根据当前的需求,按需付费。我们还将讨论如何在云中获取 Python 代码。在快节奏的世界里,这是一个快速发展的行业。我们在这里有很多选择。亚马逊网络服务 ( AWS )故意不在本书中讨论,因为其他书籍,如用 Python 构建机器学习系统威利·里歇特和路易斯·佩德罗·科埃略帕克特出版都非常详细地涵盖了这个话题。

datasciencetoolbox.org/提供的数据科学工具箱,是一个基于 Linux 的数据分析虚拟环境,可以在本地运行,也可以在 AWS 上运行。数据科学工具箱网站上给出的说明非常清楚,应该可以帮助您建立一个预装了大量 Python 包的环境。

本章将涉及的主题如下:

  • 用 Matlab/Octave 交换信息
  • 安装rpy2包装
  • 与 R 接口
  • 将 NumPy 数组发送到 Java
  • 集成 SWIG 和 NumPy
  • 集成 Boost 和 Python
  • 通过f2py使用 Fortran 代码
  • PythonAnywhere 云

用 Matlab/Octave 交换信息

Matlab 及其开源替代软件 Octave 是流行的数值程序和编程语言。Octave 和 Matlab 的语法与 Python 非常相似。事实上,你可以找到比较它们语法的网站(例如,参见wiki.scipy.org/NumPy_for_M…)。

www.gnu.org/software/oc…下载 Octave。

撰写本文时使用的 Octave 版本是 4.2.0。scipy.io.savemat()函数将数组保存在符合 Octave 和 Matlab 格式的文件中。该函数接受文件名和带有数组名称的字典,并将值作为参数。参考本书代码包中的ch-11.ipynb文件:

import statsmodels.api as sm 
from scipy.io import savemat 

data_loader = sm.datasets.sunspots.load_pandas() 
df = data_loader.data 
savemat("sunspots", {"sunspots": df.values}) 

前面的代码将太阳黑子数据存储在一个名为sunspots.mat的文件中。扩展名会自动添加。启动八度图形用户界面 ( 图形用户界面)或命令行界面。加载我们创建的文件并查看数据,如下所示:

    octave:1> load sunspots.mat
    octave:2> sunspots
    sunspots =
       1.7000e+03   5.0000e+00
       1.7010e+03   1.1000e+01
       1.7020e+03   1.6000e+01

安装 rpy2 包

R 编程语言在统计学家中很受欢迎。它是用 C 和 Fortran 编写的,可以在 GNU 通用公共许可证 ( GPL 下获得。r 支持建模、统计测试、时间序列分析、分类、可视化和聚类。综合 R 档案网 ( CRAN )等资源库网站为各种任务提供数千个 R 包。

www.r-project.org/下载 R。

rpy2包方便了与 Python 中的 R 接口。用pip按照以下步骤安装rpy2:

$ pip3 install rpy2

如果您已经安装了rpy2,请按照rpy.sourceforge.net/rpy2/doc-de…上的说明进行操作,因为升级不是一个简单的过程。

与 R 接口

r 提供了一个包含样本数据集的datasets包。morley数据集有 1879 年光速测量的数据。光速是一个基本的物理常数,它的值目前已知非常精确。数据描述见。光速值可以在scipy.constants模块中找到。R 数据存储在具有三列的 R 数据帧中:

  • 实验编号,从一到五
  • 运行次数,每个实验运行 20 次,使测量总数达到 100
  • 以千米每秒为单位的光速,减去 299,000

rpy2.robjects.r()函数在 Python 环境中执行 R 代码。按照以下方式加载数据:

pandas2ri.activate() 
r.data('morley') 

Pandas 库通过pandas.rpy.common模块的 R 接口被否决,因此建议读者使用rpy2对象模块。将数据加载到 Pandas 数据框中,如下所示:

df = r['morley'] 

让我们用以下代码通过实验对数据进行分组,这将创建一个 5 乘 2 的 NumPy 数组:

samples = dict(list(df.groupby('Expt'))) 
samples = np.array([samples[i]['Speed'].values for i in samples.keys()]) 

当我们有不同实验的数据时,知道这些实验的数据点是否来自同一个分布是很有趣的。 Kruskal-Wallis 单向方差分析(参考http://en . Wikipedia . org/wiki/Kruskal % E2 % 80% 93 Wallis _ 单向方差分析 _of_variance )是一种统计方法,分析样本时不需要对其分布做假设。这个测试的零假设是所有样本的中位数都相等。测试在scipy.stats.kruskal()功能中实现。按照以下步骤进行测试:

print("Kruskal", kruskal(samples[0], samples[1], samples[2], samples[3], samples[4])) 

测试statisticpvalue打印在下面一行:

Kruskal KruskalResult(statistic=15.022124661246552,    pvalue=0.0046555484175328015)

我们可以拒绝零假设,但这并不能告诉我们哪个或哪些实验有偏离的中位数。进一步的分析留给读者去做。如果我们绘制每个实验的最小值、最大值和平均值,我们会得到下图:

Interfacing with R

查看本书代码包中的ch-11.ipynb文件:

import rpy2.robjects as ro 
from rpy2.robjects import pandas2ri 
from rpy2.robjects import r 

from scipy.stats import kruskal 
import matplotlib.pyplot as plt 
import numpy as np 
from scipy.constants import c 

pandas2ri.activate() 
r.data('morley') 

df = r['morley'] 

df['Speed'] = df['Speed'] + 299000 

samples = dict(list(df.groupby('Expt'))) 
samples = np.array([samples[i]['Speed'].values for i in  
samples.keys()]) 
print("Kruskal", kruskal(samples[0], samples[1], samples[2], samples[3], samples[4])) 

plt.title('Speed of light') 
plt.plot(samples.min(axis=1), 'x', label='min') 
plt.plot(samples.mean(axis=1), 'o', label='mean') 
plt.plot(np.ones(5) * samples.mean(), '--', label='All mean') 
plt.plot(np.ones(5) * c/1000, lw=2, label='Actual') 
plt.plot(samples.max(axis=1), 'v', label='max') 
plt.grid(True) 
plt.legend() 
plt.show() 

向 Java 发送 NumPy 数组

和 Python 一样,Java 也是一种非常流行的编程语言。我们在第 8 章中安装了 Java,作为使用 Cassandra 的先决条件。要运行 Java 代码,我们需要 Java 运行时环境 ( JRE )。开发需要 Java 开发工具包 ( JDK )。

Jython 是用 Java 编写的 Python 的实现。Jython 代码可以使用任何 Java 类。但是,用 C 语言编写的 Python 模块不能在 Jython 中导入。这是一个问题,因为许多数值和数据分析 Python 库都有用 c 语言编写的模块。JPype1包提供了一个解决方案,可以从pypi.python.org/pypi/JPype1github.com/originell/j…下载。您可以使用以下命令安装JPype1:

$ pip3 install JPype1 

用以下代码启动 Java 虚拟机 ( JVM ):

jpype.startJVM(jpype.getDefaultJVMPath()) 

用一些随机值创建一个 JPype 数组JArray:

values = np.random.randn(7) 
java_array = jpype.JArray(jpype.JDouble, 1)(values.tolist()) 

按如下方式打印每个数组元素:

for item in java_array: 
   jpype.java.lang.System.out.println(item) 

最后,我们应该关闭 JVM,如下所示:

jpype.shutdownJVM() 

以下是本书代码包中ch-11.ipynb文件的代码列表:

import jpype 
import numpy as np 
from numpy import random 
jpype.startJVM(jpype.getDefaultJVMPath()) 

random.seed(44) 
values = np.random.randn(7) 
java_array = jpype.JArray(jpype.JDouble, 1)(values.tolist()) 

for item in java_array: 
   jpype.java.lang.System.out.println(item) 

jpype.shutdownJVM() 

当您在 Jupyter 笔记本中执行上述代码时,您会在 Jupyter 控制台中获得以下输出:

[W 18:23:45.918 NotebookApp] 404 GET /nbextensions/widgets/notebook/js/extension.js?v=20170305114358 (::1) 4.30ms referer=http://localhost:8888/notebooks/ch-12.ipynb
-0.7506147172558728
1.3163573247118194
1.2461400286434303
-1.6049157412585944
-1.468143678979905
-1.7150704579733684
1.8587836915125544
JVM activity report     :
classes loaded       : 32
JVM has been shutdown

集成 SWIG 和 NumPy

c 是 1970 年左右开发的一种广泛使用的编程语言。C 语言存在多种方言,C 语言影响了其他编程语言。c 不是面向对象的。这导致了 C++的产生,这是一种具有 C 特性的面向对象语言,因为 C 是 C++的子集。C 和 C++是编译语言。我们需要编译源代码来创建所谓的目标文件。之后,我们必须链接对象文件来创建动态共享库。

集成 C 和 Python 的好处是我们有很多选择。第一个选项是简化包装器和接口生成器 ( SWIG )。SWIG 在开发过程中增加了一个额外的步骤,那就是在 Python 和 C(或 C++)之间生成胶水代码。从www.swig.org/download.ht…下载 SWIG。在撰写本文时,最新的 SWIG 版本是 3.0.12。安装 SWIG 的前提是安装 Perl 兼容正则表达式 ( PCRE )。PCRE 是一个 C 正则表达式库。从www.pcre.org/下载 PCRE。撰写本文时的 PCRE 版本为 8.39。打开 PCRE 后,运行以下命令:

$ ./configure
$ make
$ make install

前面片段中的最后一个命令需要rootsudo访问。我们可以用同样的命令安装 SWIG。我们从编写包含函数定义的头文件开始。编写一个头文件,定义以下函数:

double sum_rain(int* rain, int len); 

我们将使用前面的函数对我们在上一章中分析的rain金额值求和。请参考本书代码包中的sum_rain.h文件。该功能在本书的代码包中的sum_rain.cpp文件中实现:

double sum_rain(int* rain, int len) { 

  double sum = 0.; 

  for (int i = 0; i < len; i++){ 
    if(rain[i] == -1) { 
       sum += 0.025; 
    } else { 
      sum += 0.1 * rain[i]; 
    } 
  } 

  return sum; 
} 

定义以下 SWIG 接口文件(参考本书代码包中的sum_rain.i文件):

%module sum_rain 

%{ 
  #define SWIG_FILE_WITH_INIT 
  #include "sum_rain.h" 
%} 

%include "/tmp/numpy.i" 

%init %{ 
  import_array(); 
%} 

%apply (int* IN_ARRAY1, int DIM1) {(int* rain, int len)}; 

%include "sum_rain.h" 

前面的代码依赖于numpy.i界面文件,可以在https://github . com/numpy/numpy/blob/master/tools/swig/numpy . I找到。在这个例子中,文件被放在/tmp目录中,但是我们几乎可以把这个文件放在任何地方。使用以下命令生成 SWIG 粘合代码:

$ swig -c++ -python sum_rain.i

上一步创建了一个sum_rain_wrap.cxx文件。如下编译sum_rain.cpp文件:

$ g++ -O2 -fPIC -c sum_rain.cpp -I<Python headers dir>

在前面的命令中,我们需要指定实际的 Python C 头文件目录。我们可以通过以下命令找到它:

$ python3-config --includes

我们也可以使用以下命令进行编译:

$ g++ -O2 -fPIC -c sum_rain.cpp $(python3-config --includes)

该目录的位置将根据 Python 版本和操作系统而有所不同(类似于/usr/include/python3.6)。如下编译生成的 SWIG 包装文件:

$ g++ -O2 -fPIC -c sum_rain_wrap.cxx $(python3-config --includes)  -I<numpy-dir>/core/include/

前面的命令取决于安装的 NumPy 的位置。从 Python 外壳中找到它,如下所示:

$ python3
>>> import numpy as np
>>> np.__file__

屏幕上打印的字符串应该包含 Python 版本,site-packages,以__init__.pyc结尾。如果去掉最后一部分,我们应该有 NumPy 目录。或者,我们可以使用以下代码:

>>> from imp import find_module
>>> find_module('numpy')

在我的计算机上,发出以下命令:

$ g++ -O2 -fPIC -c sum_rain_wrap.cxx $(python3-config --includes) -I/usr/local/lib/python3.6/site-packages/numpy/core/include/

最后一步是链接通过编译创建的目标文件:

$ g++ -lpython3.6 -dynamiclib sum_rain.o sum_rain_wrap.o -o _sum_rain.so -L /usr/local/Cellar/python3/3.6.0/Frameworks/Python.framework/Versions/3.6/lib

除非我们使用 Cygwin,否则前面的步骤在其他操作系统(如窗口)上的工作方式会有所不同。如果需要,建议您在 SWIG 用户邮件列表(www.swig.org/mail.html)或 StackOverflow 上寻求帮助。

用本书代码包中的swig_demo.py文件测试创建的库:

from _sum_rain import * 
import numpy as np 

rain = np.load('rain.npy') 
print("Swig", sum_rain(rain)) 
rain = .1 * rain 
rain[rain < 0] = .025 
print("Numpy", rain.sum()) 

使用以下命令执行此文件:

$ python3 swig_demo.py

如果一切顺利,并且我们没有混淆 Python 安装,将打印以下行:

Swig 85291.554999999328
Numpy 85291.55

集成 Boost 和 Python

Boost 是一个可以和 Python 接口的 C++库。从www.boost.org/users/downl…下载。撰写本文时的 Boost 版本是 1.63.0。最简单但也是最慢的安装方法包括以下命令:

$ ./bootstrap.sh --prefix=/path/to/boost
$ ./b2 install

prefix参数指定安装目录。在本例中,我们将假设 Boost 安装在名为 Boost 的目录(如~/Boost)中用户的home目录下。在该目录中,将创建一个libinclude目录。对于 Unix 和 Linux,您应该运行以下命令:

export LD_LIBRARY_PATH=$HOME/Boost/lib:${LD_LIBRARY_PATH}

在 Mac OS X 上,设置以下环境变量:

export DYLD_LIBRARY_PATH=$HOME/Boost/lib 

在我们的例子中,我们将这个变量设置如下:

export DYLD_LIBRARY_PATH=/usr/local/Cellar/boost/1.63.0/lib

重新定义本书代码包中boost_rain.cpp文件中给出的 rain 求和函数:

#include <boost/python.hpp> 

double sum_rain(boost::python::list rain, int len) { 

  double sum = 0.; 

  for (int i = 0; i < len; i++){ 
    int val = boost::python::extract<int>(rain[i]); 
    if(val == -1) { 
       sum += 0.025; 
    } else { 
      sum += 0.1 * val; 
    } 
  } 

  return sum; 
} 

BOOST_PYTHON_MODULE(librain) { 
    using namespace boost::python; 

    def("sum_rain", sum_rain); 
} 

该函数接受 Python 列表和列表的大小。从 Python 中调用该函数,如本书代码包中的rain_demo.py文件所示:

import numpy as np 
import librain 

rain_data = np.load('rain.npy') 
print("Boost", librain.sum_rain(rain_data.astype(int).tolist(), len(rain_data))) 
rain_data = .1 * rain_data 
rain_data[rain_data < 0] = .025 
print("Numpy", rain_data.sum()) 

我们将使用本书代码包中的Makefile文件来自动化开发过程:

CC = g++ 
PYLIBPATH = $(shell python3-config --exec-prefix)/lib 
LIB = -L$(PYLIBPATH) $(shell python3-config --libs) -L/usr/local/Cellar/boost/1.63.0/lib -L/usr/local/Cellar/boost-python/1.63.0/lib -lboost_python3 
OPTS = $(shell python3-config --include) -O2 -I/usr/local/Cellar/boost/1.63.0/include 

default: librain.so 
  @python3 ./rain_demo.py 

librain.so: rain.o 
  $(CC) $(LIB)  -Wl,-rpath,$(PYLIBPATH) -shared $< -o $@ 

rain.o: boost_rain.cpp Makefile 
  $(CC) $(OPTS) -c $< -o $@ 

clean: 
  rm -rf *.so *.o 

.PHONY: default clean 

从命令行运行以下命令:

$ make clean;make

结果如预期的一样:

Boost 85291.54999999328
Numpy 85291.55 

通过 f2py 使用 Fortran 代码

Fortran (源自公式翻译)是一种成熟的编程语言,主要用于科学计算。它是在 20 世纪 50 年代开发的,出现了更新的版本,如 Fortran 77、Fortran 90、Fortran 95、Fortran 2003 和 Fortran 2008(参考en.wikipedia.org/wiki/Fortra…)。每个版本都增加了特性和新的编程范例。这个例子我们需要一个 Fortran 编译器。gfortran编译器是 GNU Fortran 编译器,可以从gcc.gnu.org/wiki/GFortr…下载。

NumPy f2py模块作为 Fortran 和 Python 之间的接口。如果有一个 Fortran 编译器,我们可以使用这个模块从 Fortran 代码创建一个共享库。我们将编写一个 Fortran 子程序,用于对前面例子中给出的降雨量值进行求和。定义子例程,并将其存储在 Python 字符串中。之后,我们可以调用f2py.compile()函数,从 Fortran 代码中产生一个共享库。最终产品在本书代码包的fort_src.py文件中:

from numpy import f2py 
fsource = ''' 
       subroutine sumarray(A, N) 
       REAL, DIMENSION(N) :: A 
       INTEGER :: N 
       RES = 0.1 * SUM(A, MASK = A .GT. 0) 
       RES2 = -0.025 * SUM(A, MASK = A .LT. 0) 
       print*, RES + RES2 
       end  
 ''' 
f2py.compile(fsource,modulename='fort_sum',verbose=0) 

使用以下命令执行文件以生成模块:

$ python3 fort_src.py

调用本书代码包中fort_demo.py文件给出的子程序:

import fort_sum 
import numpy as np 
rain = np.load('rain.npy') 
fort_sum.sumarray(rain, len(rain)) 
rain = .1 * rain 
rain[rain < 0] = .025 
print("Numpy", rain.sum()) 

使用以下命令执行文件以生成输出:

$ python3 fort_demo.py

Fortran 和 NumPy 的结果符合预期(我们可以忽略 Fortran 子程序打印的最后两位数字):

85291.5547
Numpy 85291.55

PythonAnywhere 云

Python 这里是 Python 开发的云服务。该界面完全基于网络,模拟了 Bash、Python 和 IPython 控制台。Python 网络环境中预先安装的 Python 库列在www.pythonanywhere.com/batteries_i…上。

软件版本可能会稍微落后于可用的最新稳定版本。在撰写本文时,从 PythonAnywhere Bash 控制台安装 Python 软件似乎有点问题,不建议这样做。

当您第一次访问网址www.pythonanywhere.com/login/时,您将看到以下屏幕以登录到 PythonAnywhere 环境:

PythonAnywhere Cloud

提交登录名和密码后,您将看到以下 web 应用程序屏幕:

PythonAnywhere Cloud

建议您上传 Python 源文件,而不是使用 PythonAnywhere 环境,因为它的响应速度不如我们的本地环境。通过点击网络应用程序中的文件标签上传文件。从本章上传bn_demo.py文件:

PythonAnywhere Cloud

要执行程序,点击bn_demp.py文件,然后点击屏幕右上角的运行按钮:

PythonAnywhere Cloud

您将在控制台的代码列表下方看到输出:

PythonAnywhere Cloud

总结

在这一章中,我们考察了 Python 的边界。在 Python 生态系统之外,诸如 R、C、Java 和 Fortran 等编程语言相当流行。我们查看了提供连接 Python 和外部代码的粘合剂的库,R 的rpy2,C 的【SWIG 和 Boost】,Java 的JPype,以及 Fortran 的f2py。云计算旨在通过互联网提供计算能力。本文还简要介绍了 Python 的一个专门的云计算服务。

下一章,第 12 章性能调优、性能分析和并发,给出了提高性能的提示。通常,我们可以通过使用并行化或用 c 语言重写部分代码来优化代码,从而加快 Python 代码的速度。我们还将讨论几种分析工具和并发 API。

十二、性能调整、性能分析和并发性

|   | “过早优化是万恶之源” |   | |   | - 唐纳德·克努特,著名的计算机科学家和数学家 |

对于现实世界的应用程序,性能与特性、健壮性、可维护性、可测试性和可用性一样重要。性能与应用程序的可伸缩性成正比。结束这本书而不考虑性能增强从来都不是一个选择。事实上,为了避免过早优化,我们将性能这个话题推迟到了本书的最后一章讨论。在本章中,我们将给出使用概要分析作为关键技术来提高性能的提示。我们还将讨论多核分布式系统的相关框架。我们将在本章中讨论以下主题:

  • 分析代码
  • 安装 Cython
  • 调用 C 代码
  • 使用多处理创建池进程
  • 加速与 Joblib 令人尴尬的平行for循环
  • 瓶颈函数与 NumPy 函数的比较
  • 使用 Jug 执行 MapReduce
  • 为 Python 安装 MPI
  • IPython 并行

分析代码

分析包括识别需要性能调整的代码部分,因为它们要么太慢,要么使用大量资源,如处理器功率或内存。我们将从第 9 章分析文本数据和社交媒体中概述情绪分析代码的修改版本。代码经过重构,符合多处理编程准则(您将在本章后面学习多处理)。我们还简化了停止词过滤。第三个变化是在代码中有更少的单词特征,这样减少不会影响准确性。最后一个变化影响最大。原始代码运行了大约 20 秒。新代码运行速度比这更快,并将作为本章的基准。一些更改与性能分析有关,将在本节后面解释。请参考本书代码包中的prof_demo.py文件:

import random 
from nltk.corpus import movie_reviews 
from nltk.corpus import stopwords 
from nltk import FreqDist 
from nltk import NaiveBayesClassifier 
from nltk.classify import accuracy 

import builtins 

# Define profile function so the profile decorator  
#   does not return error when code is not profiled 
try: 
    profile = builtins.profile 
except AttributeError: 
    def profile(func):  
        return func 

@profile 
def label_docs(): 
    docs = [(list(movie_reviews.words(fid)), cat) 
            for cat in movie_reviews.categories() 
            for fid in movie_reviews.fileids(cat)] 
    random.seed(42) 
    random.shuffle(docs) 

    return docs 

@profile 
def isStopWord(word): 
    return word in sw or len(word) == 1 

@profile 
def filter_corpus(): 
    review_words = movie_reviews.words() 
    print("# Review Words", len(review_words)) 
    res = [w.lower() for w in review_words if not isStopWord(w.lower())] 
    print("# After filter", len(res)) 

    return res 

@profile 
def select_word_features(corpus): 
    words = FreqDist(corpus) 
    N = int(.02 * len(words.keys())) 
    return list(words.keys())[:N] 

@profile 
def doc_features(doc): 
    doc_words = FreqDist(w for w in doc if not isStopWord(w)) 
    features = {} 
    for word in word_features: 
        features['count (%s)' % word] = (doc_words.get(word, 0)) 
    return features 

@profile 
def make_features(docs): 
    return [(doc_features(d), c) for (d,c) in docs] 

@profile 
def split_data(sets): 
    return sets[200:], sets[:200] 

if __name__ == "__main__": 
    labeled_docs = label_docs() 

    sw = set(stopwords.words('english')) 
    filtered = filter_corpus() 
    word_features = select_word_features(filtered) 
    featuresets = make_features(labeled_docs) 
    train_set, test_set = split_data(featuresets) 
    classifier = NaiveBayesClassifier.train(train_set) 
    print("Accuracy", accuracy(classifier, test_set)) 
    print(classifier.show_most_informative_features())  
    print(classifier.show_most_informative_features()) 

当我们测量时间时,让尽可能少的进程运行是有帮助的。然而,我们不能确定后台没有运行,所以我们将使用time命令从三次测量中获取最低时间。这是一个可以在各种操作系统和 Cygwin 上使用的命令。按如下方式运行命令:

$ time python3 prof_demo.py

我们得到一个real时间,这是我们用时钟测量的时间。usersys时间测量程序使用的中央处理器时间。sys时间是在内核中度过的时间。在我的机器上,获得了以下以秒为单位的时间(括号内为最低值):

| **时间类型** | **运行 1** | **运行 2** | **运行 3** | | `real` | Eleven point five two one | Ten point eight zero eight | (10.416) | | `user` | Nine point seven five eight | Nine point eight two six | (9.444) | | `sys` | Zero point nine six five | Zero point six four three | (0.620) |

使用标准 Python 分析器对代码进行分析,如下所示:

$ python3 -m cProfile -o /tmp/stat.prof prof_demo.py

-o开关指定输出文件。我们可以用gprof2dot PyPi 包可视化分析器输出。按照以下步骤安装:

$ pip3 install gprof2dot 

使用以下命令创建巴布亚新几内亚可视化:

$ gprof2dot -f pstats /tmp/stat.prof |dot -Tpng -o /tmp/cprof.png

如果出现dot: command not found错误,说明没有安装 Graphviz。你可以从www.graphviz.org/Download.ph…下载 Graphviz。

完整图像太大,无法显示。这是它的一小部分:

Profiling the code

按如下方式查询探查器输出:

$ python3 -m pstats /tmp/stat.prof

使用此命令,我们进入配置文件统计浏览器。从输出中去除文件名,按时间排序,并显示前 10 次:

/tmp/stat.prof% strip
/tmp/stat.prof% sort time
/tmp/stat.prof% stats 10

有关最终结果,请参考以下屏幕截图:

Profiling the code

以下是标题的描述:

| **表头** | **描述** | | `ncalls` | 这是电话号码。 | | `tottime` | 这是在给定函数中花费的总时间(不包括调用子函数的时间)。 | | `percall` | 这是`tottime`除以`ncalls`的商。 | | `cumtime` | 这是在这个和所有子功能中花费的总时间(从调用到退出)。即使对于递归函数,这个数字也是准确的。 | | `percall`(秒) | 这是`cumtime`除以原始调用的商。 |

使用以下命令退出配置文件浏览器:

/tmp/stat.prof% quit

line_profiler是我们可以使用的另一个剖析器。这个分析器可以显示已经用@profile装饰器装饰的函数中每一行的统计数据。使用以下命令安装并运行此探查器:

$ pip3 install line_profiler
$ kernprof -l -v prof_demo.py

完整的报告太长,无法在此复制;相反,下面是每个功能的摘要(有一些重叠):

Function: label_docs at line 9 Total time: 6.19904 s
Function: isStopWord at line 19 Total time: 2.16542 s
File: prof_demo.py Function: filter_corpus at line 23
Function: select_word_features at line 32 Total time: 4.05266 s
Function: doc_features at line 38 Total time: 12.5919 s
Function: make_features at line 46 Total time: 14.566 s
Function: split_data at line 50 Total time: 3.6e-05 s

安装 Cython

Cython 编程语言充当了 Python 和 C/C++之间的粘合剂。有了 Cython 工具,我们可以从普通的 Python 代码生成 C 代码,然后可以编译成二进制,这更接近机器级别。cytoolz包包含由 Cythonizing 方便的 Python toolz包创建的实用程序。以下命令将安装cythoncytoolz:

$ pip3 install cython cytoolz

就像在烹饪节目中一样,我们将在完成相关过程(推迟到下一部分)之前展示 Cythonizing 的结果。timeit Python 模块测量时间。我们将使用这个模块来测量不同的功能。定义以下函数,该函数接受短代码片段、函数调用以及代码作为参数运行的次数:

def time(code, n): 
    times = min(timeit.Timer(code, setup=setup).repeat(3, n)) 

    return round(1000* np.array(times)/n, 3) 

接下来,我们预定义一个包含所有代码的大设置字符串。代码在本书代码包的timeits.py文件中(代码使用cython_module,在您的机器上构建):

import timeit 
import numpy as np 

setup = ''' 
import nltk 
import cython_module as cm 
import collections 
from nltk.corpus import stopwords 
from nltk.corpus import movie_reviews 
from nltk.corpus import names 
import string 
import pandas as pd 
import cytoolz 

sw = set(stopwords.words('english')) 
punctuation = set(string.punctuation) 
all_names = set([name.lower() for name in names.words()]) 
txt = movie_reviews.words(movie_reviews.fileids()[0]) 

def isStopWord(w): 
    return w in sw or w in punctuation 

def isStopWord2(w): 
    return w in sw or w in punctuation or not w.isalpha() 

def isStopWord3(w): 
    return w in sw or len(w) == 1 or not w.isalpha() or w in all_names 

def isStopWord4(w): 
    return w in sw or len(w) == 1 

def freq_dict(words): 
    dd = collections.defaultdict(int) 

    for word in words: 
        dd[word] += 1 

    return dd 

def zero_init(): 
    features = {} 

    for word in set(txt): 
        features['count (%s)' % word] = (0) 

def zero_init2(): 
    features = {} 
    for word in set(txt): 
        features[word] = (0) 

keys = list(set(txt)) 

def zero_init3(): 
    features = dict.fromkeys(keys, 0) 

zero_dict = dict.fromkeys(keys, 0) 

def dict_copy(): 
    features = zero_dict.copy() 
''' 

def time(code, n): 
    times = min(timeit.Timer(code, setup=setup).repeat(3, n)) 

    return round(1000* np.array(times)/n, 3) 

if __name__ == '__main__': 
    print("Best of 3 times per loop in milliseconds") 
    n = 10 
    print("zero_init ", time("zero_init()", n)) 
    print("zero_init2", time("zero_init2()", n)) 
    print("zero_init3", time("zero_init3()", n)) 
    print("dict_copy ", time("dict_copy()", n)) 
    print("\n") 

    n = 10**2 
    print("isStopWord ", time('[w.lower() for w in txt if not isStopWord(w.lower())]', n)) 
    print("isStopWord2", time('[w.lower() for w in txt if not isStopWord2(w.lower())]', n)) 
    print("isStopWord3", time('[w.lower() for w in txt if not isStopWord3(w.lower())]', n)) 
    print("isStopWord4", time('[w.lower() for w in txt if not isStopWord4(w.lower())]', n)) 
    print("Cythonized isStopWord", time('[w.lower() for w in txt if not cm.isStopWord(w.lower())]', n)) 
    print("Cythonized filter_sw()", time('cm.filter_sw(txt)', n)) 
    print("\n") 
    print("FreqDist", time("nltk.FreqDist(txt)", n)) 
    print("Default dict", time('freq_dict(txt)', n)) 
    print("Counter", time('collections.Counter(txt)', n)) 
    print("Series", time('pd.Series(txt).value_counts()', n)) 
    print("Cytoolz", time('cytoolz.frequencies(txt)', n)) 
    print("Cythonized freq_dict", time('cm.freq_dict(txt)', n)) 

因此,我们有几个isStopword()函数版本,运行时间以毫秒为单位如下:

isStopWord  0.843
isStopWord2 0.902
isStopWord3 0.963
isStopWord4 0.869
Cythonized isStopWord 0.924
Cythonized filter_sw() 0.887

作为比较,我们还有一个普通pass语句的运行时间。cytnonizedisStopWord()基于isStopWord3()功能(最精细的滤镜)。如果我们看看prof_demo.py中的doc_features()功能,很明显我们不应该仔细检查每个单词的特征。相反,我们应该只是将文档中的一组单词与作为特征选择的单词相交。所有其他字数都可以安全地设置为零。事实上,最好我们将所有的值初始化为零一次,然后复制这个字典。对于相应的函数,我们得到以下执行时间:

zero_init  0.61
zero_init2 0.555
zero_init3 0.017
dict_copy  0.011

另一个改进是使用 Python defaultdict类,而不是 NLTK FreqDist类。相关例程具有以下运行时间:

FreqDist 2.206
Default dict 0.674
Counter 0.79
Series 7.006
Cytoolz 0.542
Cythonized freq_dict 0.616

正如我们所看到的,Cythonized 版本的速度一直在加快,尽管并没有加快多少。

调用 C 代码

我们可以从 Cython 调用 C 函数。C 字符串strlen()函数相当于 Python len()函数。通过导入以下文件,从 Cython .pyx文件调用该函数:

from libc.string cimport strlen 

然后我们可以从.pyx文件中的其他地方调用strlen().pyx文件可以包含任何 Python 代码。请看一下本书代码包中的cython_module.pyx文件:

from collections import defaultdict 
from nltk.corpus import stopwords 
from nltk.corpus import names 
from libc.string cimport strlen 

sw = set(stopwords.words('english')) 
all_names = set([name.lower() for name in names.words()]) 

def isStopWord(w): 
     py_byte_string = w.encode('UTF-8') 
     cdef char* c_string = py_byte_string 
     truth = (w in sw) or (w in all_names) or (not w.isalpha()) or (strlen(c_string) == 1) 
     return truth 

def filter_sw(words): 
    return [w.lower() for w in words if not isStopWord(w.lower())] 

def freq_dict(words): 
    dd = defaultdict(int) 

    for word in words: 
        dd[word] += 1 

    return dd 

要编译这段代码,我们需要一个包含以下内容的setup.py文件:

from distutils.core import setup 
from Cython.Build import cythonize 

setup( 
    ext_modules = cythonize("cython_module.pyx") 
) 

使用以下命令编译代码:

$ python3 setup.py build_ext --inplace

以下输出显示了cython_module扩展的构建:

$python3 setup.py build_ext --inplace
running build_ext
building 'cython_module' extension
creating build
creating build/temp.macosx-10.12-x86_64-3.6
clang -Wno-unused-result -Wsign-compare -Wunreachable-code -fno-common -dynamic -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -I/usr/local/include -I/usr/local/opt/openssl/include -I/usr/local/opt/sqlite/include -I/usr/local/Cellar/python3/3.6.0/Frameworks/Python.framework/Versions/3.6/include/python3.6m -c cython_module.c -o build/temp.macosx-10.12-x86_64-3.6/cython_module.o
clang -bundle -undefined dynamic_lookup build/temp.macosx-10.12-x86_64-3.6/cython_module.o -L/usr/local/lib -L/usr/local/opt/openssl/lib -L/usr/local/opt/sqlite/lib -o /Users/armando/gdrive/projects/bitbucket/pubs/2016-pda-e2-packt/chapters/ch-12/cython_module.cpython-36m-darwin.so

我们现在可以修改情绪分析程序来调用 Cython 函数。我们还将添加上一节中提到的改进。由于我们要反复使用一些函数,这些函数被提取到本书代码包的core.py文件中。查看本书代码包中的cython_demo.py文件(代码使用cython_module,在您的机器上构建):

from nltk.corpus import movie_reviews 
from nltk import NaiveBayesClassifier 
from nltk.classify import accuracy  
import cython_module as cm 
import cytoolz 
from core import label_docs 
from core import filter_corpus 
from core import split_data 

def select_word_features(corpus): 
    words = cytoolz.frequencies(filtered) 
    sorted_words = sorted(words, key=words.get) 
    N = int(.02 * len(sorted_words)) 

    return sorted_words[-N:] 

def match(a, b): 
    return set(a.keys()).intersection(b) 

def doc_features(doc): 
    doc_words = cytoolz.frequencies(cm.filter_sw(doc)) 

    # initialize to 0 
    features = zero_features.copy() 

    word_matches = match(doc_words, word_features) 

    for word in word_matches: 
        features[word] = (doc_words[word]) 

    return features 

def make_features(docs): 
    return [(doc_features(d), c) for (d,c) in docs] 

if __name__ == "__main__": 
    labeled_docs = label_docs() 
    filtered = filter_corpus() 
    word_features = select_word_features(filtered) 
    zero_features = dict.fromkeys(word_features, 0) 
    featuresets = make_features(labeled_docs) 
    train_set, test_set = split_data(featuresets) 
    classifier = NaiveBayesClassifier.train(train_set) 
    print("Accuracy", accuracy(classifier, test_set)) 
    print(classifier.show_most_informative_features()) 

使用 time 命令执行代码,如下所示:

$ time python3 cython_demo.py

下表总结了time命令的结果(括号内为最低值):

| **时间类型** | **运行 1** | **运行 2** | **运行 3** | | `real` | (9.639) | Nine point eight one seven | Nine point nine one two | | `user` | (9.604) | Nine point six six one | Nine point six eight three | | `sys` | (0.404) | Zero point four two four | Zero point four five one |

与之前的代码执行相比,我们可以看到性能的提升。为了便于比较,下面的时序表是从上一节复制过来的:

| **时间类型** | **运行 1** | **运行 2** | **运行 3** | | `real` | Eleven point five two one | Ten point eight zero eight | (10.416) | | `user` | Nine point seven five eight | Nine point eight two six | (9.444) | | `sys` | Zero point nine six five | Zero point six four three | (0.620) |

创建具有多处理的进程池

多处理是一个标准的 Python 模块,目标是具有多个处理器的机器。多处理通过创建多个进程围绕全局解释器锁 ( GIL )工作。

GIL 锁 Python 字节码,因此只有一个线程可以访问它。

多处理支持进程池、队列和管道。进程池是可以并行执行功能的系统进程池。队列是通常用于存储任务的数据结构。管道连接不同的进程,使得一个进程的输出成为另一个进程的输入。

Windows 没有os.fork()功能,所以我们需要确保只有导入和def块被定义在if __name__ == "__main__"块之外。

创建一个池并注册一个函数,如下所示:

   p = mp.Pool(nprocs) 

该池有一个map()方法,它是 Python map()函数的并行等价物:

p.map(simulate, [i for i in xrange(10, 50)]) 

我们将模拟一维粒子的运动。粒子执行随机游走,我们感兴趣的是计算粒子的平均末端位置。我们对不同的行走长度重复这个模拟。计算本身并不重要。重要的部分是比较多进程和单个进程的加速比。我们将使用 matplotlib 绘制加速图。完整的代码在本书代码包的multiprocessing_sim.py文件中:

from numpy.random import random_integers 
from numpy.random import randn 
import numpy as np 
import timeit 
import argparse 
import multiprocessing as mp 
import matplotlib.pyplot as plt 

def simulate(size): 
    n = 0 
    mean = 0 
    M2 = 0 

    speed = randn(10000) 

    for i in range(1000):  
        n = n + 1 
        indices = random_integers(0, len(speed)-1, size=size) 
        x = (1 + speed[indices]).prod() 
        delta = x - mean 
        mean = mean + delta/n 
        M2 = M2 + delta*(x - mean) 

    return mean 

def serial(): 
    start = timeit.default_timer() 

    for i in range(10, 50): 
        simulate(i) 

    end = timeit.default_timer() - start 
    print("Serial time", end) 

    return end 
def parallel(nprocs): 
    start = timeit.default_timer() 
    p = mp.Pool(nprocs) 
    print(nprocs, "Pool creation time", timeit.default_timer() - start) 

    p.map(simulate, [i for i in range(10, 50)]) 
    p.close() 
    p.join() 

    end = timeit.default_timer() - start 
    print(nprocs, "Parallel time", end) 
    return end 

if __name__ == "__main__": 
    ratios = [] 
    baseline = serial() 

    for i in range(1, mp.cpu_count()): 
        ratios.append(baseline/parallel(i)) 

    plt.xlabel('# processes') 
    plt.ylabel('Serial/Parallel') 
    plt.plot(np.arange(1, mp.cpu_count()), ratios) 
    plt.grid(True) 
    plt.show() 

如果我们取范围从 1 到 8 的进程池大小的加速值(处理器的数量取决于硬件),我们得到下图:

Creating a process pool with multiprocessing

阿姆达尔定律(见en.wikipedia.org/wiki/Amdahl…)最能描述并行化带来的加速。这个定律预测了最大可能的加速。进程的数量限制了绝对最大加速。然而,正如我们在前面的图中看到的,我们没有用两个进程获得双倍的速度,也没有用三个进程获得三倍的速度,但是我们接近了。任何给定 Python 代码的某些部分都可能无法并行化。例如,我们可能需要等待资源变得可用,或者我们可能正在执行一个必须顺序执行的计算。我们还必须考虑来自并行化设置和相关进程间通信的开销。阿姆达尔定律指出,加速的倒数、进程数的倒数和代码中不能并行的部分之间存在线性关系。

加快了与 Joblib 循环的并行速度,令人尴尬

Joblib 是由 scikit-learn 的开发人员创建的 Python 库。它的主要任务是提高长时间运行的 Python 函数的性能。Joblib 通过使用多处理或幕后线程的缓存和并行化实现了这些改进。按照以下步骤安装 Joblib:

$ pip3 install joblib

我们将重复使用前面例子中的代码,只改变parallel()函数。参考本书代码包中的joblib_demo.py文件:

def parallel(nprocs): 
    start = timeit.default_timer() 
    Parallel(nprocs)(delayed(simulate)(i) for i in xrange(10, 50)) 

    end = timeit.default_timer() - start 
    print(nprocs, "Parallel time", end) 
    return end 

最终结果见下图(处理器数量取决于硬件):

Speeding up embarrassingly parallel for loops with Joblib

比较瓶颈和 NumPy 函数

瓶颈是一组受 NumPy 和 SciPy 启发的函数,但在 Cython 中编写时考虑到了高性能。瓶颈为数组维度、轴和数据类型的每个组合提供单独的 Cython 函数。这不会显示给最终用户,瓶颈的限制因素是要确定执行哪个 Cython 功能。安装瓶颈,如下所示:

$ pip3 install bottleneck

我们将比较numpy.median()scipy.stats.rankdata()函数相对于它们的瓶颈对应物的执行时间。在紧密循环或经常调用的函数中使用 Cython 函数之前,手动确定它可能是有用的。

该程序在本书代码包的bn_demo.py文件中给出:

import bottleneck as bn 
import numpy as np 
import timeit 

setup = ''' 
import numpy as np 
import bottleneck as bn 
from scipy.stats import rankdata 

np.random.seed(42) 
a = np.random.randn(30) 
''' 
def time(code, setup, n): 
    return timeit.Timer(code, setup=setup).repeat(3, n) 

if __name__ == '__main__': 
    n = 10**3 
    print(n, "pass", max(time("pass", "", n))) 
    print(n, "min np.median", min(time('np.median(a)', setup, n))) 
    print(n, "min bn.median", min(time('bn.median(a)', setup, n))) 
    a = np.arange(7) 
    print)"Median diff", np.median(a) - bn.median(a)) 

    print(n, "min scipy.stats.rankdata", min(time('rankdata(a)', setup, n))) 
    print(n, "min bn.rankdata", min(time('bn.rankdata(a)', setup, n))) 

以下是带有运行时间和函数名的输出:

$ python3 bn_demo.py 
1000 pass 7.228925824165344e-06
1000 min np.median 0.019842895912006497
1000 min bn.median 0.0003261091187596321
Median diff 0.0
1000 min scipy.stats.rankdata 0.04070987505838275
1000 min bn.rankdata 0.0011222949251532555 

显然,瓶颈很快;不幸的是,由于它的设置,瓶颈还没有那么多功能。下表列出了从pypi.python.org/pypi/Bottle…实现的功能:

| **类别** | **功能** | | NumPy/SciPy | `median`、`nanmedian`、`rankdata`、`ss`、`nansum`、`nanmin`、`nanmax`、`nanmean`、`nanstd`、`nanargmin`和`nanargmax` | | 功能 | `nanrankdata`、`nanvar`、`partsort`、`argpartsort`、`replace`、`nn`、`anynan`和`allnan` | | 移动窗口 | `move_sum`、`move_nansum`、`move_mean`、`move_nanmean`、`move_median`、`move_std`、`move_nanstd`、`move_min`、`move_nanmin`、`move_max`和`move_nanmax` |

用 Jug 执行 MapReduce

Jug 是一个分布式计算框架,使用任务作为中央并行化单元。Jug 使用文件系统或 Redis 服务器作为后端。Redis 服务器在第 8 章使用数据库中进行了讨论。使用以下命令安装 Jug:

$ pip3 install jug

MapReduce (参见【http://en.wikipedia.org/wiki/MapReduce】)是一种分布式算法,用于用一组计算机处理大型数据集。该算法包括一个映射和一个减少阶段。在地图阶段,数据以并行方式处理。数据被分成多个部分,在每个部分上,执行过滤或其他操作。在缩减阶段,来自映射阶段的结果被聚合,例如,创建一个统计报告。

如果我们有一个文本文件列表,我们可以计算每个文件的字数。这可以在映射阶段完成。最后,我们可以将单个的字数合并成一个语料库词频词典。Jug 具有 MapReduce 功能,这在本书的代码包中的jug_demo.py文件中进行了演示(代码取决于cython_module工件):

import jug.mapreduce 
from jug.compound import CompoundTask 
import cython_module as cm 
import cytoolz 
import pickle 

def get_txts(): 
    return [(1, 'Lorem ipsum dolor sit amet, consectetur adipiscing elit.'), (2, 'Donec a elit pharetra, malesuada massa vitae, elementum dolor.'), (3, 'Integer a tortor ac mi vehicula tempor at a nunc.')] 

def freq_dict(file_words): 
    filtered = cm.filter_sw(file_words[1].split()) 

    fd = cytoolz.frequencies(filtered) 

    return fd 

def merge(left, right): 
    return cytoolz.merge_with(sum, left, right) 

merged_counts = CompoundTask(jug.mapreduce.mapreduce, merge, freq_dict, get_txts(), map_step=1) 

在前面的代码中,merge()函数在减少阶段被调用,freq_dict()函数在映射阶段被调用。我们定义了一个由多个子任务组成的 Jug CompoundTask。在运行这段代码之前,我们需要启动一个 Redis 服务器。通过发出以下命令来执行 MapReduce:

$ jug execute jug_demo.py --jugdir=redis://127.0.0.1/&

结尾的&符号(&)表示该命令在后台运行。如果 Redis 服务器可以在网络中访问,我们可以通过这种方式从多台计算机发出命令。在这个例子中,Redis 只在本地机器上运行(127.0.0.1是本地主机的 IP 地址)。但是,我们仍然可以在本地多次运行该命令。我们可以如下检查jug命令的状态:

$ jug status jug_demo.py

默认情况下,如果我们不指定jugdir选项,Jug 会将数据存储在当前工作目录中。使用以下命令清理jug目录:

$ jug cleanup jug_demo.py

为了查询 Redis 并执行剩下的分析,我们将使用另一个程序。在这个程序中,初始化 Jug 如下:

jug.init('jug_demo.py', 'redis://127.0.0.1/') 
import jug_demo 

下面一行是减少阶段的结果:

words = jug.task.value(jug_demo.merged_counts) 

其余代码在本书代码包的jug_redis.py文件中给出:

import jug 

def main(): 
    jug.init('jug_demo.py', 'redis://127.0.0.1/') 
    import jug_demo 
    print("Merged counts", jug.task.value(jug_demo.merged_counts)) 

if __name__ == "__main__": 
    main() 

为 Python 安装 MPI

消息传递接口 ( MPI )(参见en.wikipedia.org/wiki/Messag…)是由专家开发的标准协议,用于广泛的分布式机器。最初,在九十年代,MPI 被用来编写 Fortran 和 c 语言的程序,MPI 独立于硬件和编程语言。MPI 功能包括发送和接收操作、MapReduce 功能和同步。MPI 具有涉及两个处理器的点对点功能,以及涉及所有处理器的操作。MPI 有几种编程语言的绑定,包括 Python。从www.open-mpi.org/software/om…下载 MPI。MPI 2.0.2 在编写本报告时已经安装并使用;我们可以在网站上查看是否有更新的版本。安装 MPI 可能需要一段时间(近 30 分钟)。以下是所涉及的命令,假设我们将其安装在/usr/local目录中:

$ ./configure --prefix=/usr/local
$ make all
$ sudo make install

为 MPI 安装 Python 绑定,如下所示:

$ pip3 install mpi4py

IPython 并行

IPython Parallel 是用于并行计算的 IPython API。我们将设置它使用 MPI 进行消息传递。我们可能必须按如下方式设置环境变量:

$ export LC_ALL=en_US.UTF-8
$ export LANG=en_US.UTF-8

在命令行中发出以下命令:

$ ipython3 profile create --parallel --profile=mpi

前面的命令将在位于home目录的.ipython/profile_mpi文件夹中创建几个文件。

按如下方式启动使用 MPI 配置文件的集群:

$ ipcluster start --profile=mpi --engines=MPI --debug

前面的命令指定我们使用的是带有调试级日志记录的mpi概要文件和 MPI 引擎。我们现在可以通过 IPython 笔记本与集群进行交互。启动笔记本,启用绘图,并自动导入 NumPy、SciPy 和 matplotlib,如下所示:

$ jupyter-notebook --profile=mpi --log-level=DEBUG 

上述命令使用调试日志级别的mpi配置文件。本例的笔记本存储在本书代码包的IPythonParallel.ipynb文件中。导入 IPython 并行Client类和statsmodels.api模块,如下所示:

  In [1]:from ipyparallel import Client 
  import statsmodels.api as sm 

加载太阳黑子数据并计算平均值:

 In [2]: data_loader = sm.datasets.sunspots.load_pandas() 
 vals = data_loader.data['SUNACTIVITY'].values 
 glob_mean = vals.mean() 
 glob_mean 

输出如下:

Out [2]: 49.752103559870541

按照以下方式创建客户端:

In [3]: c = Client(profile='mpi') 

使用以下行创建客户端视图:

In [4]: view=c[:] 

IPython 使用魔法的概念。这些是 IPython 笔记本特有的特殊命令。按如下方式启用 magics:

In [5]: view.activate() 

在本书的代码包中加载mpi_ipython.py文件:

from mpi4py import MPI 
from numpy.random import random_integers 
from numpy.random import randn 
import numpy as np 
import statsmodels.api as sm 
import bottleneck as bn 
import logging 

def jackknife(a, parallel=True): 
    data_loader = sm.datasets.sunspots.load_pandas() 
    vals = data_loader.data['SUNACTIVITY'].values 

    results = [] 

    for i in a: 
        tmp = np.array(vals.tolist()) 
        tmp[i] = np.nan 
        results.append(bn.nanmean(tmp)) 
    results = np.array(results) 

    if parallel: 
        comm = MPI.COMM_WORLD 
        rcvBuf = np.zeros(0.0, 'd') 
        comm.gather([results, MPI.DOUBLE], [rcvBuf, MPI.DOUBLE]) 

   return results 

if __name__ == "__main__": 
    skiplist = np.arange(39, dtype='int') 
    print(jackknife(skiplist, False)) 

前面的程序包含一个执行刀切重采样的功能。刀切重采样是重采样的一种,我们省略样本中的一个观测值,然后计算我们感兴趣的统计估计量。在这种情况下,我们对平均值感兴趣。我们将一个观察值设置为 NumPy NaN,从而将其忽略。然后,我们在新样本上调用瓶颈nanmean()函数。以下是加载命令:

In [6]: view.run('mpi_ipython.py') 

接下来,我们分割并展开一个包含太阳黑子数组所有指数的数组:

In [7]: view.scatter('a',np.arange(len(vals),dtype='int')) 

a数组可以在笔记本中显示如下:

In [8]: view['a'] 

以下是前面命令的输出:

Out[8]:[array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,  17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,  34, 35, 36, 37, 38]), ... TRUNCATED ...]

在所有客户端上调用jackknife()功能:

In [9]: %px means = jackknife(a) 

完成所有工作进程后,我们可以查看结果:

In [10]: view['means'] 

结果列表显示了与我们开始时一样多的进程。每个进程返回一个 NumPy 数组,其中包含通过刀切重采样计算的平均值。这个结构不是很有用,所以把它转换成一个平面列表:

In [11]: all_means = [] 

for v in view['means']: 
    all_means.extend(v) 

mean(all_means) 

您将获得以下输出:

Out [11]: 49.752103559870577

我们也可以计算标准差,但这很容易,所以我们将跳过它。画出折刀的直方图更有趣:

In [13]: hist(all_means, bins=sqrt(len(all_means))) 

最终结果见下图:

对于故障排除,我们可以使用下面一行显示工作进程的错误消息:

In [14]: [(k, c.metadata[k]['started'], c.metadata[k]['pyout'], c.metadata[k]['pyerr']) for k in c.metadata.keys()] 

总结

在这一章中,我们从第 9 章分析文本数据和社交媒体调整了情绪分析脚本的表现。使用概要分析、Cython 和各种改进,我们将该示例的执行速度提高了一倍。我们还通过 IPython Parallel 使用了多处理、Joblib、Jug 和 MPI 来利用并行化。

这是这本书的最后一章。当然,学习的过程不应该停止。更改代码以满足您的需求。有一个私人的数据分析项目总是很好的,即使只是为了练习。如果想不出项目,可以参加www.kaggle.com/上的比赛。他们有好几个比赛,奖品很好。

十三、附录 a:关键概念

本附录给出了本书使用的技术概念的简要概述和词汇表。

阿姆达尔定律预测了由于并行化可能带来的最大加速。进程的数量限制了绝对最大加速。任何给定 Python 代码的某些部分都可能无法并行化。我们还必须考虑并行化设置和相关进程间通信的开销。阿姆达尔定律指出,加速的倒数、进程数的倒数和不能并行化的代码部分之间存在线性关系。

ARMA 模型结合了自回归和移动平均模型。它们用于预测时间序列的未来值。

人工神经网络 ( ANN )是受动物大脑启发的模型。神经网络是由神经元组成的网络,神经元是具有输入和输出的单元。一个神经元的输出可以传递给一个神经元等等,从而形成一个多层网络。神经网络包含自适应元素,使其适合处理非线性模型和模式识别问题。

增广 Dickey Fuller ( ADF )检验是一个与协整相关的统计检验。ADF 检验用于检验时间序列的平稳性。

自相关是数据集内的相关性,可以指示趋势。例如,如果我们有一个周期的滞后,我们可以检查前一个值是否影响当前值。要做到这一点,自相关值必须相当高。

自相关绘制不同滞后时间序列数据的自相关图。自相关是时间序列与相同滞后时间序列的相关性。

自回归模型是一种使用(通常是线性)回归来使用以前的值预测时间序列的未来值的模型。自回归模型是 ARMA 模型的一个特例。它们相当于零移动平均成分的 ARMA 模型。

单词包模型是一个简化的文本模型,其中文本由一个单词包来表示。在此表示中,单词的顺序被忽略。通常,字数或某些单词的存在被用作该模型中的特征。

气泡图是散点图的延伸。在气泡图中,第三个变量的值由数据点周围气泡的大小表示。

卡珊德拉查询语言 ( CQL )是一种针对 Apache 卡珊德拉的查询语言,语法类似于 SQL。

协整类似于相关性,是时间序列数据的一个统计特征。协整是衡量两个时间序列同步程度的指标。

聚类旨在将数据划分为称为聚类的组。在训练数据没有标记的意义上,聚类通常是无监督的。一些聚类算法需要猜测聚类的数量,而其他算法则不需要。

层叠样式表 ( CSS )是一种用于设计网页元素样式的语言。CSS 由万维网联盟维护和开发。

CSS 选择器是用于选择网页内容的规则。

字符代码包含在 NumPy 中,以向后兼容数字。NumPy 的前身是 NumPy。

数据类型对象是 numpy.dtype 类的实例。它们为 NumPy 数据类型的操作提供了一个面向对象的接口。

特征值是方程 Ax = ax 的标量解,其中 A 是二维矩阵,x 是一维向量。

特征向量是对应于特征值的向量。

指数移动平均线是一种权重随时间呈指数递减的移动平均线。

快速傅里叶变换 ( FFT )是一种计算傅里叶变换的快速算法。FFT 是 O(N log N),这是对旧算法的巨大改进。

滤波是一种信号处理技术,包括去除或抑制部分信号。存在许多滤波器类型,包括中值滤波器和维纳滤波器。

傅立叶分析是基于以数学家约瑟夫·傅立叶命名的傅立叶级数。傅立叶级数是一种将函数表示为正弦和余弦项的无穷级数的数学方法。所讨论的函数可以是实数或复数。

遗传算法基于生物进化论。这种类型的算法对于搜索和优化很有用。

图形处理器单元 ( 图形处理器)是用于高效显示图形的专用电路。最近,图形处理器被用来执行大规模并行计算(例如,训练神经网络)。

分层数据格式 ( HDF )是存储大数值数据的规范和技术。HDF 集团维护着一个相关的软件库。

希尔伯特-黄变换是一种分解信号的数学算法。该方法可用于检测时间序列数据中的周期性周期。它被成功地用于确定太阳黑子周期。

超文本标记语言 ( HTML )是用来创建网页的基础技术。它为媒体、文本和超链接定义标签。

互联网工程任务组 ( IETF )是一个致力于维护和发展互联网的开放团体。IETF 是开放的,原则上任何人都可以加入。

JavaScript 对象符号 ( JSON )是一种数据格式。在这种格式中,数据是用 JavaScript 符号写下来的。JSON 比 XML 等其他数据格式更简洁。

k-fold 交叉验证是一种涉及 k(一个小整数)个随机数据分区的交叉验证形式,称为 fold。在 k 次迭代中,每个折叠使用一次进行验证,其余的数据用于训练。迭代的结果可以在最后合并。

Kruskal-Wallis 单向方差分析是一种分析样本方差的统计方法,不需要对样本的分布做假设。

滞后图是一个时间序列和同一时间序列滞后的散点图。滞后图显示了特定滞后时间序列数据中的自相关。

学习曲线是一种可视化学习算法行为的方式。它是一系列训练数据大小的训练和测试分数图。

对数图(或对数图)是使用对数刻度的图。当数据变化很大时,这种类型的图很有用,因为它们显示数量级。

逻辑回归是一种分类算法。这种算法可以用来预测一个类或事件发生的概率。逻辑回归基于逻辑函数,其值在 0 到 1 之间,就像概率一样。因此,逻辑函数可用于将任意值转换为概率。

MapReduce 是一种分布式算法,用于处理具有一组计算机的大型数据集。该算法由映射阶段和缩减阶段组成。在地图阶段,数据以并行方式处理。数据被分成多个部分,在每个部分上,执行过滤或其他操作。在缩减阶段,将汇总地图阶段的结果。

摩尔定律是指现代计算机芯片中的晶体管数量每两年翻一番的观察结果。自 1970 年左右摩尔定律形成以来,这一趋势一直在持续。还有第二个摩尔定律,也叫洛克定律。该定律指出,集成电路的研发和制造成本呈指数级增长。

移动平均线指定了一个以前看到的数据窗口,每次该窗口向前滑动一个周期时,该窗口将被平均。不同类型的移动平均在用于平均的权重上有本质的不同。

朴素贝叶斯分类是一种基于概率论和统计学中贝叶斯定理的概率分类算法。它被称为幼稚,因为它有很强的独立性假设。

对象关系映射 ( ORM )是一种用于在数据库模式和面向对象编程语言之间进行转换的软件架构模式。

观点挖掘情感分析是一个以高效发现和评价文本中的观点和情感为目标的研究领域。

词性 ( POS )标签是句子中每个单词的标签。这些标签有语法意义,如动词或名词。

表示状态转移 ( REST )是一种用于 web 服务的架构风格。

真正简单的联合 ( RSS )是博客等网络订阅源的发布和检索标准。

散点图是显示笛卡尔坐标系中两个变量之间关系的二维图。一个变量的值用一个轴表示,另一个变量的值用另一个轴表示。我们可以通过这种方式快速可视化相关性。

信号处理是工程和应用数学的一个领域,处理模拟和数字信号的分析,对应随时间变化的变量。

SQL 是关系数据库查询和操作的专用语言。这包括创建表、在表中插入行和删除表。

Stopwords 是低信息值的常用词。通常在分析文本之前会删除停用词。虽然过滤停用词是一种常见的做法,但停用词没有标准的定义。

监督学习是一种需要标注训练数据的机器学习类型。

支持向量机 ( SVM )可用于回归(SVR)和分类(SVC)。SVM 将数据点映射到多维空间中的点。映射由所谓的内核函数执行。核函数可以是线性的或非线性的。

词频-逆文档频率 ( tf-idf )是一个度量一个单词在语料库中重要性的指标。它由术语频率号和逆文档频率号组成。术语“频率”统计单词在文档中出现的次数。逆文档频率计算单词出现的文档数,取该数的倒数。

一个时间序列是一个有序的数据点列表,从最早的测量值开始。通常,每个数据点都有一个相关的时间戳。时间序列可以是平稳的,也可以是非平稳的。

十四、附录 b:实用功能

本附录列出了 matplotlib、NumPy、Pandas、scikit-learn 和 SciPy 的软件包组织的有用功能。

Matplotlib

以下是有用的 matplotlib 函数:

matplotlib.pyplot.axis(*v, **kwargs):这是获取或设置轴属性的方法。例如,轴(“关闭”)关闭轴线和标签。

matplotlib.pyplot.figure(num=None, figsize=None, dpi=None, facecolor=None, edgecolor=None, frameon=True, FigureClass=<class 'matplotlib.figure.Figure'>, **kwargs):此功能创建一个新的图形。

matplotlib.pyplot.grid(b=None, which='major', axis='both', **kwargs):此功能打开或关闭绘图网格。

matplotlib.pyplot.hist(x, bins=10, range=None, normed=False, weights=None, cumulative=False, bottom=None, histtype='bar', align='mid', orientation='vertical', rwidth=None, log=False, color=None, label=None, stacked=False, hold=None, **kwargs):该函数绘制直方图。

matplotlib.pyplot.imshow(X, cmap=None, norm=None, aspect=None, interpolation=None, alpha=None, vmin=None, vmax=None, origin=None, extent=None, shape=None, filternorm=1, filterrad=4.0, imlim=None, resample=None, url=None, hold=None, **kwargs):此功能显示类似数组数据的图像。

matplotlib.pyplot.legend(*args, **kwargs):该功能在可选的指定位置(例如plt.legend(loc='best'))显示一个图例。

matplotlib.pyplot.plot(*args, **kwargs):这个函数创建一个二维图,其中有一个或多个(x,y)对和一个相应的可选格式字符串。

matplotlib.pyplot.scatter(x, y, s=20, c='b', marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, verts=None, hold=None, **kwargs):该函数创建两个数组的散点图。

matplotlib.pyplot.show(*args, **kw):此功能显示一个图。

matplotlib.pyplot.subplot(*args, **kwargs):如果给定了图的行号、列号和索引号,这个函数会创建子图。所有这些数字都是从一开始的。例如,plt.subplot(221)在一个二乘二的网格中创建第一个子场景。

matplotlib.pyplot.title(s, *args, **kwargs):这个功能给剧情加上一个标题。

NumPy

以下是有用的 NumPy 函数:

numpy.arange([start,] stop[, step,], dtype=None):此函数创建一个 NumPy 数组,其值在指定范围内均匀分布。

numpy.argsort(a, axis=-1, kind='quicksort', order=None):该函数返回对输入数组进行排序的索引。

numpy.array(object, dtype=None, copy=True, order=None, subok=False, ndmin=0):这个函数从类似数组的序列比如 Python 列表中创建一个 NumPy 数组。

numpy.dot(a, b, out=None):这个函数计算两个数组的点积。

numpy.eye(N, M=None, k=0, dtype=<type 'float'>):此函数返回单位矩阵。

numpy.load(file, mmap_mode=None):该函数从加载 NumPy 数组或腌制对象。npy,。npz,或者泡菜。内存映射数组存储在文件系统中,不必完全加载到内存中。这对于大型数组尤其有用。

numpy.loadtxt(fname, dtype=<type 'float'>, comments='#', delimiter=None, converters=None, skiprows=0, usecols=None, unpack=False, ndmin=0):该函数将文本文件中的数据加载到 NumPy 数组中。

numpy.mean(a, axis=None, dtype=None, out=None, keepdims=False):此函数计算沿给定轴的算术平均值。

numpy.median(a, axis=None, out=None, overwrite_input=False):此函数计算沿给定轴的中位数。

numpy.ones(shape, dtype=None, order='C'):这个函数创建一个指定形状和数据类型的 NumPy 数组,包含一个。

numpy.polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False):该函数执行最小二乘多项式拟合。

numpy.reshape(a, newshape, order='C'):这个函数改变 NumPy 数组的形状。

numpy.save(file, arr):该功能将 NumPy 数组保存为 NumPy .npy格式的文件。

numpy.savetxt(fname, X, fmt='%.18e', delimiter=' ', newline='\n', header='', footer='', comments='# '):该函数将 NumPy 数组保存到文本文件中。

numpy.std(a, axis=None, dtype=None, out=None, ddof=0, keepdims=False):此函数返回沿给定轴的标准偏差。

numpy.where(condition, [x, y]):该函数根据布尔条件从输入数组中选择数组元素。

numpy.zeros(shape, dtype=float, order='C'):这个函数创建一个指定形状和数据类型的 NumPy 数组,包含零。

Pandas

以下是有用的 Pandas 功能:

pandas.date_range(start=None, end=None, periods=None, freq='D', tz=None, normalize=False, name=None, closed=None):这个函数创建一个固定频率的日期时间索引

pandas.isnull(obj):此函数查找 NaN 和 None 值

pandas.merge(left, right, how='inner', on=None, left_on=None, right_on=None, left_index=False, right_index=False, sort=False, suffixes=('_x', '_y'), copy=True):此函数将数据框对象与列或索引上类似数据库的连接合并

pandas.pivot_table(data, values=None, rows=None, cols=None, aggfunc='mean', fill_value=None, margins=False, dropna=True):这个函数创建一个类似电子表格的数据透视表作为 Pandas 数据框

pandas.read_csv(filepath_or_buffer, sep=',', dialect=None, compression=None, doublequote=True, escapechar=None, quotechar='"', quoting=0, skipinitialspace=False, lineterminator=None, header='infer', index_col=None, names=None, prefix=None, skiprows=None, skipfooter=None, skip_footer=0, na_values=None, na_fvalues=None, true_values=None, false_values=None, delimiter=None, converters=None, dtype=None, usecols=None, engine='c', delim_whitespace=False, as_recarray=False, na_filter=True, compact_ints=False, use_unsigned=False, low_memory=True, buffer_lines=None, warn_bad_lines=True, error_bad_lines=True, keep_default_na=True, thousands=Nment=None, decimal='.', parse_dates=False, keep_date_col=False, dayfirst=False, date_parser=None, memory_map=False, nrows=None, iterator=False, chunksize=None, verbose=False, encoding=None, squeeze=False, mangle_dupe_cols=True, tupleize_cols=False, infer_datetime_format=False):这个函数从 CSV 文件创建一个数据帧

pandas.read_excel(io, sheetname, **kwds):该函数将 Excel 工作表读入数据框

pandas.read_hdf(path_or_buf, key, **kwargs):这个函数从 HDF 商店返回一个 Pandas 对象

pandas.read_json(path_or_buf=None, orient=None, typ='frame', dtype=True, convert_axes=True, convert_dates=True, keep_default_dates=True, numpy=False, precise_float=False, date_unit=None):这个函数从一个 JSON 字符串创建一个 pandas 对象

pandas.to_datetime(arg, errors='ignore', dayfirst=False, utc=None, box=True, format=None, coerce=False, unit='ns', infer_datetime_format=False):此函数将字符串或字符串列表转换为日期时间

科学学习

以下是有用的 scikit 学习功能:

sklearn.cross_validation.train_test_split(*arrays, **options):该功能将数组拆分为随机训练集和测试集

sklearn.metrics.accuracy_score(y_true, y_pred, normalize=True, sample_weight=None):此功能返回准确率分类得分

sklearn.metrics.euclidean_distances (X, Y=None, Y_norm_squared=None, squared=False):该函数计算输入数据的距离矩阵

黑桃

本节显示了有用的 SciPy 函数:

scipy.fftpack

fftshift(x, axes=None):该功能将零频率分量移至频谱中心

rfft(x, n=None, axis=-1, overwrite_x=0):该函数对包含实值的数组进行离散傅里叶变换

scipy .信号

detrend(data, axis=-1, type='linear', bp=0):此函数从数据中移除线性趋势或常数

medfilt(volume, kernel_size=None):该函数对数组应用中值滤波

wiener(im, mysize=None, noise=None):该函数对数组应用维纳滤波器

scipy . state

anderson(x, dist='norm'):该函数对来自指定分布的数据执行安德森-达林测试

kruskal(*args):该功能对数据进行克鲁斯卡尔-沃利斯 H 测试

normaltest(a, axis=0):该功能测试数据是否符合正态分布

scoreatpercentile(a, per, limit=(), interpolation_method='fraction'):该函数计算输入数组中指定百分比的分数

shapiro(x, a=None, reta=False):该函数应用夏皮罗-维尔克检验进行正态性

十五、附录 c:在线资源

以下是文档、论坛、文章和其他信息的链接列表:

阿帕奇卡珊德拉数据库:cassandra.apache.org

靓汤:www.crummy.com/software/Be…

HDF 集团网站:www.hdfgroup.org/

有趣的 ipython 笔记本图库:https://github . com/IPython/IPython/wiki/A-有趣的 IPython 笔记本图库

Graphviz 开源图形可视化软件:graphviz.org/

IPython 网站:ipython.org/

Jupyter 网站:jupyter.org/

Matplotlib(一个 Python 绘图库):matplotlib.org/

MongoDB(开源文档数据库):www.mongodb.org

mpi4 py 文档:mpi4py.scipy.org/docs/usrman…

自然语言工具包(NLTK):www.nltk.org/

NumPy 和 SciPy 文档:docs.scipy.org/doc/

NumPy 和 SciPy 邮件列表:www.scipy.org/Mailing_Lis…

打开 MPI(高性能消息传递库):www.open-mpi.org

Packt 发布帮助和支持:www.packtpub.com/support

Pandas 主页:pandas.pydata.org

Python 性能提示:wiki.python.org/moin/Python…

Redis(开源键值存储):redis.io/

Scikit-learn(Python 机器学习):scikit-learn.org/stable/

Scikit-学习性能提示:scikit-learn.org/stable/deve…

SciPy 性能提示:wiki.scipy.org/Performance…

SQLAlchemy(Python SQL 工具包和对象关系映射器):www.sqlalchemy.org

Toolz 实用程序功能文档:toolz.readthedocs.org/en/latest/

plotplotlib 图形转换器:plot.ly/matplotlib/…

离线使用 Python Plotly:plot.ly/python/offl…

保存静态图片(PNG、PDF 等):plot.ly/python/stat…

用 Python 创建 HTML 或 PDF 报告:plot.ly/python/#rep…

贵公司的安全和 Plotly 的服务器:plot.ly/products/on…

用 Plotly 和 Python 创建仪表盘:plot.ly/python/dash…

连接数据库:plot.ly/python/#dat…

Plotly 和 IPython / Jupyter 笔记本:plot.ly/ipython-not…