axis 改成-1 好一点
np.concatenate((xs[..., None], ys[..., None], zs[..., None]), axis=-1)

不太清楚你的问题细节, 写了个 3d 插值的例子:

# -*- coding: utf-8 -*-
import numpy as np
from scipy.interpolate import interpn
from enthought.mayavi import mlab

mlab.options.offscreen = True

original_points = 20j
interp_points = 40j
figsize = (1280, 720)

x, y, z = np.mgrid[-15:14:original_points, -12:12:original_points, -14:15:original_points]
# x.shape: 20x20x20
# y.shape: 20x20x20
# z.shape: 20x20x20
s = np.sin(x*y*z)/(x*y*z) # s.shape: 20x20x20
mlab.contour3d(x, y, z, s)
mlab.savefig("original.png", size=figsize)

xs, ys, zs = np.mgrid[-15:14:interp_points, -12:12:interp_points, -14:15:interp_points]

# Construct a 3d array of 3-dimensional points.
arr_4d = np.concatenate((xs[..., None], ys[..., None], zs[..., None]), axis=3) # shape: 40x40x40x3
ss = interpn((x[:, 0, 0], y[0, :, 0], z[0, 0, :]), s, arr_4d, method="nearest")

# Weirdly though, this will also work, why?（ﾟдﾟ）
sss = interpn((x[:, 0, 0], y[0, :, 0], z[0, 0, :]), s, (xs, ys, zs), method="nearest")
print np.all(ss == sss) # True

mlab.contour3d(xs, ys, zs, ss)
mlab.savefig("interpolate.png", size=figsize)

解决方法:

http://docs.enthought.com/mayavi/mayavi/auto/mlab_helper_functions.html#mayavi.mlab.contour3d

http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interpn.html#scipy.interpolate.interpn

@dexterzzz
有了 python 跟一系列 py 的工具还有必要学 excel 么
曾经打算学好 office 套, 后来用了 python 跟 latex 后就没管过了...

None 是两个线程中 test 的返回值, 两个线程都在抢着打印.
py3 的 print 不是函数调用不报错?

在 settings > editor > inspections > python 里面各种开关可以自己勾选, 有些勾上还是挺有用的, 例如要写兼容 py2,py3 的代码

@practicer
一般是打印出来放在枕头下面, 寂寞的时候翻一翻 2333

@protream
# For Python2
if sys.version < '3':
reload(sys)
sys.setdefaultencoding('utf-8')

看了下源码, py2 下这样做是很危险的, 会打乱整个程序的运行环境, 尤其是当使用了 sys.setdefaultencoding('utf-8')的代码被用到更大的程序中时, 程序何时崩溃都不奇怪, 一开始处理编码问题时也这样用, 直到程序莫名其妙的退出...
http://stackoverflow.com/search?q=sys.setdefaultencoding%28%27utf-8%27%29

也写个简单总结吧

unicode 字符集说白了就是 0 到 0x10ffff 这一串数字(0 到 1114111),也叫码点, 比如汉字"中", python 中的 unicode 表示为 u'\u4e2d', \u 是转义字符, 码点是 0x4e2d(十进制为 20013). 这些数字需要转化成字节序列(0-255)在内存中表示时, 就要有一个规则将 unicode 码点(比如 0x4e2d)转换为字节序列, 这个规则就是编码, 这个过程反过来时的规则叫做解码, 另外只有 UTF 家族(UTF-8, UTF-16,...)的编码器才能将任意 unicode 码点完全正确的编码为字节序列.

python 处理各种编码的一个准则: early decode, always unicode, late encode.

early decode:
在获取外部输入时(文件, 终端, 网络...)尽早的将其解码为 unicode. 如果已知外部输入的编码时, 就用那个编码将输入 decode 就是; 如果不知道的话就只能靠统计的方法检测其编码了, 没有其它完全可靠的算法(但是 utf 编码的可以直接检测 BOM), 比如大多数浏览器自动检测编码就是靠统计, 最后以某个大概率统计出具体编码, python 的一个扩展库 chardet 也是基于统计, 如果是写爬虫的话, chardet 几乎是一个必须的库吧(requests 都直接把 chardet 放在自己的包里面了), 如果是读取本地文件, py2 要依赖标准库的 codecs 解码, py3 的 open 函数自带解码参数

always unicode:
程序内部处理过程中都统一为 unicode

late encode:
程序处理完, 需要将其输出了, 无论输出在终端是还是在文件, 都是给人看的, 需要将 unicode 码点转换为字节序列, 比如说简中 windows 的控制台吧, 默认代码页是 cp932, 编码数大概是 23,940 个(跟 unicode 的 1114112 个比较下就知道为啥好多字符都不能显示在控制台了), 显示在控制台一般能起个简单观测的作用, 可以简单的 ignore 掉不能显示的字符
eg:
---------------------------------------------------------------------------------------------------------------------------------------------
# -*- encoding = utf-8 -*-
# 上面这一行的作用是提示 python 虚拟机, 它将要读取的这个代码文件是 utf-8 编码的.因为这里面会有英语以外的字符
# 这不是强制的, 如果虚拟机发现不是 utf-8 编码的, 会用默认的 ascii 读取该代码文件, 于是有不是英文的字符就会报错
# 当然你也必须保存为 utf-8, 而且是有 BOM(Byte Order Mark)的 utf-8 文件, 试试用其它编辑器改成 No BOM 或者
# 其它编码, 虚拟机检测不到文件的 BOM, 也转为使用 ascii 解码并读取该文件
import locale
local_encoding = locale.getpreferredencoding() # 获取控制台默认编码, 不一定是 cp932 呀
print local_encoding

alice = u"アリス・マーガトロイド"
# 这个日文的点・在 cp936 的控制台是显示不出来的, 如果想在控制台显示这个字符串,
# 可以在编码时选择 ignore 不能编码的字符,或者 replace(默认使用?替换)
try:
print alice
except Exception as e:
print e
print alice.encode(local_encoding, errors="ignore")
---------------------------------------------------------------------------------------------------------------------------------------------

用 windows 终端运行看看吧, 当然在 IDE 中 print alice 可能能完整显示, 因为输出被重定向了
------------------------------------------
在 python 自带的 IDLE 中运行的结果:
cp936
アリス・マーガトロイド
アリスマーガトロイド
IDLE 中输出被重定向, 点正确显示
使用 cp936 编码后点被 ignore 了
----------------------------------------
控制台运行的结果:
cp936
'gbk' codec can't encode character u'\u30fb' in position 3: illegal multibyte sequence
アリスマーガトロイド
编码失败
忽略不能编码的码点

输出到文件跟输出到控制台是一样的, 都是需要某个具体编码的地方, 当然这时可以指定能完全编码 unicode 的 UTF 系列编码器完全无误的记录下来.

文本的编码解码可以类比音视频的编码解码, 比如用一个音乐播放器播放视频, 解码器给错, 当然不能播放,decode error; 压制视频的时候, 给 video 指定一个 audio 的编码器, encode error

一直用 py2 处理各种编码, 看了 5 分钟这个流程图, 最后还是没明白是啥意思

谢谢, 看来背后要懂的东西还不少

感谢作者

请问一下, 用 ipython notebook 写书后转换成其它格式需要用些什么技巧, 比如转换成 md 或者 pdf 文件.
发现直接用 File 菜单的 Download as...或者基本的转换命令生成的格式都很糟糕, 生成的 md 文件用 markdownPad 打开后格式会出现不对地方, 尤其是代码段, 生成的 pdf 也是相当难看. 看了下 nbconvert 的文档, 试着添了些 options, 结果还是很糟

以前一直用 latex 编写文档, 各种繁琐的格式配置感觉越用越麻烦, 最近开始用 ipython notebook, 确实便利, html 格式比较好看, 转换成其他格式就变丑了

@ruoyu0088
原来您说的这个技巧写在第二版书里面了啊
拿到书刚看完了 cython 那一章, 收获不小 :-D

@ruoyu0088
对 ctypes 不是特别熟悉, ipython notebook 也很少用
翻了下文档, 如果没理解错的话:
--------------------------------------------------------------
In [1]:
%load_ext cythonmagic

In [2]:
%%cython -a
from ctypes import CFUNCTYPE, c_double, cast
from libc.math cimport sin, exp

cdef double cintegrand_simp(double x):
return sin(x)

cdef double cintegrand_comp(double x):
cdef double numerator = x * exp(x) + sin(exp(x)) + x**3 + x**5
cdef double denominator = sin(x) + 1.0 + exp(x)
return sin(numerator / denominator)

func_type = CFUNCTYPE(c_double, c_double)

# get function addresses
cdef int simp_func_addr = <int><void*>cintegrand_simp
cdef int comp_func_addr = <int><void*>cintegrand_comp

# cast the addresses to ctypes functions
simp_integrand = cast(simp_func_addr, func_type)
comp_integrand = cast(comp_func_addr, func_type)

In [3]:
from math import exp, sin, pi
from scipy.integrate import quad

def pyfunc_simp(x):
"""A simple python integrand function."""
return sin(x)

def pyfunc_comp(x):
"""A seemed complicated python integrand function."""
numerator = x * exp(x) + sin(exp(x)) + x**3 + x**5
denominator = sin(x) + 1.0 + exp(x)
return sin(numerator / denominator)

In [5]:
print quad(simp_integrand, 0.0, pi)
print quad(pyfunc_simp, 0.0, pi)
print quad(comp_integrand, 0.0, pi)
print quad(pyfunc_comp, 0.0, pi)
(2.0, 2.220446049250313e-14)
(2.0, 2.220446049250313e-14)
(1.0010180983222745, 8.364968411064218e-09)
(1.0010180983222745, 8.364968411064218e-09)
In [6]:
%timeit quad(simp_integrand, 0.0, pi)
100000 loops, best of 3: 3.55 µs per loop
In [7]:
%timeit quad(pyfunc_simp, 0.0, pi)
100000 loops, best of 3: 7.12 µs per loop
In [8]:
%timeit quad(comp_integrand, 0.0, pi)
10000 loops, best of 3: 32.6 µs per loop
In [9]:
%timeit quad(pyfunc_comp, 0.0, pi)
10000 loops, best of 3: 168 µs per loop
--------------------------------------------------------------
这儿编译好的比之前的稍显慢一点呢? 加上编译器优化选项后还是一样的结果.
翻 ctypes 文档时顺便解决了之前的一个问题:
我用 C 实现了一个算法, 算法需要接收一个可调用的函数对象,之前一直都是在 cython 中完成的,
但是有些时候又想能直接传一个 py 的函数对象过去,之前想尝试的办法是利用 Python.h 里面的
函数得到 py 函数对象的函数指针(scipy 似乎是用的这种方式),
现在发现 ctypes 可以很简单的完成这个任务!ヽ(・∀・)ノ

@yanyuechuixue
@ruoyu0088
根据楼上的提示验证了一下(仅一维的情况):
-----------cintegrand.c--------------------------------------------------
#include <math.h>
double cintegrand_simp(double x)
{
// A simple c function
return sin(x);
}

double cintegrand_comp(double x)
{
// A seemed complicated c function
double numerator = x * exp(x) + sin(exp(x)) + pow(x, 3) + pow(x, 5);
double denominator = sin(x) + 1.0 + exp(x);
return sin(numerator / denominator);
}

-------------compile commands--------------------------------------------
gcc -c -O3 --fast-math cintegrand.c -o cintegrand.o
gcc -shared -fPIC -o cintegrand.dll cintegrand.o

--------------py script--------------------------------------------------
import ctypes
from math import exp, sin, pi
from scipy.integrate import quad

def pyfunc_simp(x):
"""A simple python integrand function."""
return sin(x)

def pyfunc_comp(x):
"""A seemed complicated python integrand function."""
numerator = x * exp(x) + sin(exp(x)) + x**3 + x**5
denominator = sin(x) + 1.0 + exp(x)
return sin(numerator / denominator)

def load_1d_cfunc(dll_name, func_name):
"""Load cytpes function."""
dll = ctypes.CDLL(dll_name)
cfunc = dll.__getattr__(func_name)
cfunc.restype = ctypes.c_double
cfunc.argtypes = (ctypes.c_double,)
return cfunc

cfunc_simp = load_1d_cfunc("cintegrand.dll", "cintegrand_simp")
cfunc_comp = load_1d_cfunc("cintegrand.dll", "cintegrand_comp")
width = 50
print "quadrature of a simple python function.".ljust(width), quad(pyfunc_simp, 0.0, pi)
print "quadrature of a simple c function.".ljust(width), quad(cfunc_simp, 0.0, pi)
print "quadrature of a complicated python function.".ljust(width), quad(pyfunc_comp, 0.0, pi)
print "quadrature of a complicated c function.".ljust(width), quad(cfunc_comp, 0.0, pi)

---------------timeit in ipython------------------------------------------
%paste 到 ipython 中得到结果:
quadrature of a simple python function. (2.0, 2.220446049250313e-14)
quadrature of a simple c function. (2.0, 2.220446049250313e-14)
quadrature of a complicated python function. (1.0010180983222745, 8.364968411064218e-09)
quadrature of a complicated c function. (1.0010180983222743, 8.364968390887983e-09)

In [3]: %timeit quad(pyfunc_simp, 0.0, pi)
100000 loops, best of 3: 7 us per loop

In [4]: %timeit quad(cfunc_simp, 0.0, pi)
100000 loops, best of 3: 3.84 us per loop

In [5]: %timeit quad(pyfunc_comp, 0.0, pi)
10000 loops, best of 3: 176 us per loop

In [6]: %timeit quad(cfunc_comp, 0.0, pi)
10000 loops, best of 3: 24.2 us per loop
-----------------------------------------------------------------------------
简单函数有~2x 的速度提升, 这儿这个复杂函数只有~7x 的速度提升.性能区别主要在函数值计算上面
http://docs.scipy.org/doc/scipy-0.15.0/reference/tutorial/integrate.html#faster-integration-using-ctypes

测试环境:gcc (tdm-1) 5.1.0, Python 2.7.9 32bit, scipy 0.17.1

另外 scipy 是 0.15.0 开始才有支持向 scipy.integrate 传入 ctypes function 的:
http://scipy.github.io/devdocs/release.0.15.0.html#scipy-integrate-improvements
开始用 0.14.0 版本的怎么也不对, wrong ctypes function signature...o(╯□╰)o

@yanyuechuixue 也不叫用不到, 里面每行检测的代码可能都影响着最后结果的正确性跟精确性, 如果真有效率问题就不要把计算密集部分交给 py 解释器, 直接由 c/c++/fortran 这些静态编译语言完成计算, 而且对于积分是最容易并行化处理的问题了, 大多数情况单个积分求和部分没有并行化的必要, 如果是很多积分的结果再累加, 就可以考虑并行了. 如果担心自己写的不够好, 也有很多现成的底层数值计算库可以直接使用, 比如 GNU Scientific Library, C 跟 py 几乎可以完美结合(本来 py 就是一个 C 库), 我喜欢在 cython 中将两者混合使用

效率已经很高了, 但无法满足较为苛刻的计算问题, 毕竟是 python 扩展(套上了 gil 的脚镣)
scipy 的 quad 函数底层是 f77 写的, 它会根据给定的参数决定具体调用哪个 fortran 的 subroutine 积分,
def quad(func, a, b, args=(), full_output=0, epsabs=1.49e-8, epsrel=1.49e-8,
limit=50, points=None, weight=None, wvar=None, wopts=None, maxp1=50,
limlst=50):
if not isinstance(args, tuple):
args = (args,)
if (weight is None):
retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)
else:
retval = _quad_weight(func,a,b,args,full_output,epsabs,epsrel,limlst,limit,maxp1,weight,wvar,wopts)
首先分有权重跟无权重积分的情形, 这里只看看无权重的情形.无权重时调用_quad 函数,
def _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points):
......
省略多行处理 Inf 的代码
......
if points is None:
if infbounds == 0:
return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
else:
return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)
else:
if infbounds != 0:
raise ValueError("Infinity inputs cannot be used with break points.")
else:
nl = len(points)
the_points = numpy.zeros((nl+2,), float)
the_points[:nl] = points
return _quadpack._qagpe(func,a,b,the_points,args,full_output,epsabs,epsrel,limit)
现在可以看到会调用 3 个不同的 fortran 函数 _qagse, _qagie, _qagpe 之一
找到 scipy 源码中的 f77 文件, 定位到 dqagse.f, dqagie.f, dqagpe.f 这几个文件,
发现:
1. dqagse.f 中调用的 subroutine dqk21; dqk21.f 中显示使用的积分方法是 gauss-kronrod rules( https://en.wikipedia.org/wiki/Gauss%E2%80%93Kronrod_quadrature_formula)

2. dqagie.f 中调用的 subroutine dqk15i; dqk15i.f 中显示使用的积分方法也是 gauss-kronrod rules

3. dqagpe.f 中也是调用的 subroutine dqk21

由此可知, 这几个都是用的高斯积分方法, 从数值分析可以知道, 高斯方法的确比梯形法, 辛普森, 龙格库塔法有优越性: 对被积函数的拟合精度高, 收敛快. 从 fortran 代码中可以看出, 里面做了大量的错误检查跟精度的处理, 因为要适用于各种不同的平台, 以及各种不同的奇葩函数.
如果具体问题具体分析, 抛开各种检查的步骤以及平台相关的处理, 只是实现一个积分算法比这种现成的库函数更快还是极有可能的, 我用 c 简单写了一个自适应幸普森积分方法:
double
simp(integrand func, /*in*/
double left_end, /*in*/
double right_end, /*in*/
int n /*in*/
)
{
double result = 0.0;
double result0 = func(left_end) + func(right_end);
double h = (right_end - left_end) / n;
double even_sum = 0.0;
double odd_sum = 0.0;
double x;

for (int i = 0; i < n; ++i){
x = left_end + i*h;
if (i%2){
odd_sum += func(x);
}
else{
even_sum += func(x);
}
}
result = h * (result0 + 2.0 * even_sum + 4.0 * odd_sum) / 3.0;
return result;

}


double
adaptive_simp(integrand func, /*in*/
double left_end, /*in*/
double right_end, /*in*/
double tol /*in*/
)
{
int n = 20;
int max_splits = 20;
double result0, result;

result0 = simp(func, left_end, right_end, n);

for (int i = 0; i < max_splits; ++i){
n = n * 2;
result = simp(func, left_end, right_end, n);
if (fabs(result - result0) < tol)
break;
result0 = result;
}

return result;

并不困难, 很小的两个函数.

scipy 用 quadpack.py 对_quadpack.pyd 做了一层包装, 如果从 scipy import 的 quad 的话, 会比直接调用_quadpack.pyd 中相应的函数多一点点 python 的开销, 在 cython 中使用 quad, _quadpack.pyd 中的_qagse, 以及上面 c 实现的 simp(gcc 添加-O3 --fast-math 选项编译)在给定相同的误差容限下对 sinx(使用 c 标准库 sin)在[0.0, pi]区间积分 100W 次
-------------------------------------------
run quad 1000000 times, use 5.559261s
run _qagse 1000000 times, use 4.717242s
run simp 1000000 times, use 3.680429s
result from quad: 2.0
result from _qagse: 2.0
result from simp: 2.0000000001
-------------------------------------------
可以看到, 结果还是可以的.(没用 fortran 试了, 太久没写生疏了, 而且放到 cython 中还是要转一次 c, 比较麻烦)

另外在上面的 simp c 函数中累加部分可以很容易的用 openmp 并行化,曾经我对一个物体分割后积分,最后所有积分累加,这儿的累加也可以在 cython 中释放掉 gil 利用 openmp 并行化, 充分利用多核. 我那个计算问题一个结果大概要 100W 次积分, 还要放到上万次的迭代中(个人 PC 上算). 当现有的库函数无法满足效率要求时, 手写一个是明智的. 曾经还想过利用 opencl 把这个问题放到显卡上去烧, 但发现难以实现, 因为在迭代的时候总是有太多的逻辑判断跟信息交换, 只把可以纯粹的并行的部分放上去呢, 又会大大增加 CPU 跟 GPU 的信息交换, 效率又会极大的下降, 所以还是压榨的 CPU.

尝试用 memory_profiler 做一下内存分析呢?
py 采用的对象缓冲池, 垃圾回收不是立即的, 有些时候会被缓存一段时间.
曾经在多线程中下载文件, 发现如果用 with 语句打开文件写入后, 即使出了 with 的作用域, 被写入的 buffer 也不会被立即释放, 结果内存会飙升, 显式的 open, close 就没问题了