月度归档:2014年02月

试验性NumPy教程 学习笔记 1

numpy里最主要是提供了数组的支持,让一些循环操作转换为数组操作,刚开始接触的时候,习惯了常规循环的朋友,可能脑子还转不过弯来。我本人也是这样的。

所以我去找来这个教程,重点先看的数组操作这块!
根据上面的例子输入:
>>> from numpy import *
>>> a=array([20,30,40,50])
>>> b=arange(4)
>>> b
array([0, 1, 2, 3])
>>> c=a-b
>>> c
array([20, 29, 38, 47])
>>> b**2
array([0, 1, 4, 9])
>>> 10*sin(a)
array([ 9.12945251, -9.88031624, 7.4511316 , -2.62374854])
>>> a<35
array([ True, True, False, False], dtype=bool)
>>> f=a<35

对数组的运算有了一个大概的了解。

文中提到当数组非常大的时候,print会略去中间部分:
>>> print (arange(10000))
[ 0 1 2 …, 9997 9998 9999]
>>> print (arange(10000).reshape (100,100))
[[ 0 1 2 …, 97 98 99]
[ 100 101 102 …, 197 198 199]
[ 200 201 202 …, 297 298 299]
…,
[9700 9701 9702 …, 9797 9798 9799]
[9800 9801 9802 …, 9897 9898 9899]
[9900 9901 9902 …, 9997 9998 9999]]

在测试上面例子的时候,突然想到这个100×100的数组,是否可以用matplotlib来画图呢?

于是:
>>> import matplotlib.pyplot
>>> matplotlib.pyplot.subplot(221)

>>> a= (arange(10000).reshape (100,100))
>>> matplotlib.pyplot.imshow(a)

>>> matplotlib.pyplot.show

大家猜我看到了什么?
我看到了一个从上面是蓝色,到下面是红色彩虹渐变色!
加上import语句,也只用了6句话! 而这如果要在ps里,需要用渐进笔去拉一下才行。而且关键我们这个是数组,后面python程序对它有完全的掌控!

用Weave嵌入C语言 zt

用Weave嵌入C语言
Python作为动态语言其功能虽然强大,但是在数值计算方面有一个最大的缺点:速度不够快。在Python级别的循环和计算的速度只有C语言程序的百分之一。因此才有了NumPy, SciPy这样的函数库,将高度优化的C、Fortran的函数库进行包装,以供Python程序调用。如果这些高度优化的函数库无法实现我们的算法,必须从头开始写循环、计算的话,那么用Python来做显然是不合适的。因此SciPy提供了快速调用C++语言程序的方法– Weave。下面是对NumPy的数组求和的例子:

# -*- coding: utf-8 -*-
import scipy.weave as weave
import numpy as np
import time

def my_sum(a):
    n=int(len(a))
    code="""
    int i;

    double counter;
    counter =0;
    for(i=0;i<n;i++){
        counter=counter+a(i);
    }
    return_val=counter;
    """

    err=weave.inline(
        code,['a','n'],
        type_converters=weave.converters.blitz,
        compiler="gcc"
    )
    return err

a = np.arange(0, 10000000, 1.0)
# 先调用一次my_sum,weave会自动对C语言进行编译,此后直接运行编译之后的代码
my_sum(a)

start = time.clock()
for i in xrange(100):
    my_sum(a)  # 直接运行编译之后的代码
print "my_sum:", (time.clock() - start) / 100.0

start = time.clock()
for i in xrange(100):
    np.sum( a ) # numpy中的sum,其实现也是C语言级别
print "np.sum:", (time.clock() - start) / 100.0

start = time.clock()
print sum(a) # Python内部函数sum通过数组a的迭代接口访问其每个元素,因此速度很慢
print "sum:", time.clock() - start

 

此例子在我的电脑上的运行结果为:

my_sum: 0.0294527349146
np.sum: 0.0527649547638
sum: 9.11022322669
可以看到用Weave编译的C语言程序比numpy自带的sum函数还要快。而Python的内部函数sum使用数组的迭代器接口进行运算,因此速度是Python语言级别的,只有Weave版本的1/300。

weave.inline函数的第一个参数为需要执行的C++语言代码,第二个参数是一个列表,它告诉weave要把Python中的两个变量a和n传递给C++程序,注意我们用字符串表示变量名。converters.blitz是一个类型转换器,将numpy的数组类型转换为C++的blitz类。C++程序中的变量a不是一个数组,而是blitz类的实例,因此它使用a(i)获得其各个元素的值,而不是用a[i]。最后我们通过compiler参数告诉weave要采用gcc为C++编译器。如果你安装的是python(x,y)的话,gcc(mingw32)也一起安装好了,否则你可能需要手工安装gcc编译器或者微软的Visual C++。

Note 在我的电脑上,虽然安装了Visual C++ 2008 Express,但仍然提示找不到合适的Visual C++编译器。似乎必须使用编译Python的编译器版本。因此还是用gcc来的方便。
本书的进阶部分还会对weave进行详细介绍。这段程序先给了我们一个定心丸:你再也不必担心Python的计算速度不够快了。

2014.2.21 日补充:
在wakari.io网站测试,结果如下:
creating /tmp/scipy-w_skywalk-yrP2fa/python27_intermediate/compiler_ee052248824614d4d0d9f34acc183870a3803cfcd7f5ea6b10d6861381bff3cc
my_sum: 0.017
np.sum: 0.0163
4.9999995e+13
sum: 0.02

显然三者相差不多

我在自己机器上测试,暂时没有weave,结果如下:
np.sum: 0.03719650715901786
4.9999995e+13
sum: 14.177532618725465

可见numpy确实快了很多!

看看你的机器能运行吗

import numpy

LIM = 10 ** 6
N = 10 ** 9
P = 10000
primes = []
p = 2
dch=0
#By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13.
tmp=0
#What is the 10 001st prime number?

def check_primes(a, p):
   #2. Sieve out multiples of p
   a = a[a % p != 0]

   return a

for i in xrange(3, N, LIM):
   #1. Create a list of consecutive integers
   a = numpy.arange(i, i + LIM, 2)

   while len(primes) < P:
      a = check_primes(a, p)
      primes.append(p)
      p = a[0]

print len(primes), primes[P-1]
print primes[500:600]

for i in xrange (10**10+1,10**10+10**8,LIM):
    a =numpy.arange(i,i+LIM,2)

    while (tmp<10**10 * 5 ):
        a= check_primes(a, p)
        tmp=tmp+p

print tmp

我的两台机器都出错,大约是内存不够,看来4G的内存是不行的。在https://www.wakari.io 网站是可以运行的。 该程序的由来是有人问到这个题目:如何计算100亿到101亿之间的素数之和。

我的思路是,先计算101亿的开平方,也就是10万多一点,然后列出从2-10万多所有的素数,然后再从100亿开始,每个尝试去除那些素数,是素数的就累加…..现在还不太会操作….

答案是:
10000 104729
[3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727, 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821, 3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989, 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057, 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139, 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231, 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297, 4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409]
50000013737

测试

import numpy

LIM = 10 ** 6
N = 10 ** 9
P = 10000
primes = []
p = 2
dch=0
#By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13.
tmp=0
#What is the 10 001st prime number?

def check_primes(a, p):
#2. Sieve out multiples of p
a = a[a % p != 0]

return a

for i in xrange(3, N, LIM):
#1. Create a list of consecutive integers
a = numpy.arange(i, i + LIM, 2)

while len(primes) < P:
a = check_primes(a, p)
primes.append(p)
p = a[0]

print len(primes), primes[P-1]
print primes[500:600]

for i in xrange (10**10+1,10**10+10**8,LIM):
a =numpy.arange(i,i+LIM,2)

while (tmp<10**10 * 5 ):
a= check_primes(a, p)
tmp=tmp+p

print tmp

pyzo之python学习

首先pyzo里所有的包的链接:

http://www.pyzo.org/packages.html#packages

几个练习网站:

Practice

Once you know the basics, it’s good to practice your skills in order to improve them. Here are some sites where you can do that:

 

其中code academy 中可以一步一步的学习python,很推荐:http://www.codecademy.com/learn

而checkio是很专业,注册的时候需要先完成编程的小测验才行!

国内的话,可以上这个网站学习:

http://www.pythoner.cn/

普通的编程测试,可以去粉笔网站的闪电编程:

http://www.fenby.com/index

而真正用起来,当然是推荐pyzo!

python科学讲义:

http://scipy-lectures.github.io/

google+的pyzo论坛:https://groups.google.com/forum/#!forum/pyzo

pyzo国外网站下载慢的朋友,可以到360云盘这里下:

http://yunpan.cn/QpIRFbqx8TNpy

ipython 在win python3下的安装

ipython 在win python3下的安装

 

先安装python3

然后下载  distribute_setup.py

运行 python distribute_setup.py 安装setuptools

然后下载最新的ipython版本:

https://github.com/ipython/ipython/releases

解开后,到目录里,运行

python setup.py install

 

Installed d:program filespy3libsite-packagespyreadline-2.0-py3.3-win32.egg
Finished processing dependencies for ipython==1.2.0

运行下看看:

D:downloadipython-1.2.0>python ipython -pylab
WARNING: `-pylab` flag has been deprecated.
Use `–matplotlib <backend>` and import pylab manually.
Python 3.3.0 (v3.3.0:bd8afb90ebf2, Sep 29 2012, 10:55:48) [MSC v.1600 32 bit (In
tel)]
Type “copyright”, “credits” or “license” for more information.

IPython 1.2.0 — An enhanced Interactive Python.
? -> Introduction and overview of IPython’s features.
%quickref -> Quick reference.
help -> Python’s own help system.
object? -> Details about ‘object’, use ‘object??’ for extra details.
[TerminalIPythonApp] WARNING | Eventloop or matplotlib integration failed. Is ma
tplotlib installed?

奥,原来matplotlib还没有安装,这个简单

下载地址在这里:

http://matplotlib.org/downloads.html

有提示:

ImportError: matplotlib requires dateutil

如果一个一个找包太麻烦了,除了用easy_install来安装外,还可以 到这里找所有常用的已经编译好的包:

http://www.lfd.uci.edu/~gohlke/pythonlibs/

然后下载并安装dateutil

python-dateutil-2.2.win32-py3.3.exe

然后还要安装pyparsing

在matplotlib包处,就告诉了要安装以下包:

Matplotlib is a 2D plotting library
Requires numpydateutilpytzpyparsingsix, and optionally pillowpycairotornadowxpythonpysidepyqt, ghostscript, miktex, ffmpeg, mencoder, avconv, or imagemagick.

然后还要安装numpy…..

怒了,不能这样!

我们来个一揽子解决方案吧!

Scientific Python distributions

For most users, especially on Windows and Mac, the easiest way to install the packages of the SciPy stack is to download one of these Python distributions, which includes all the key packages:

  • Anaconda: A free distribution for the SciPy stack. Supports Linux, Windows and Mac.
  • Enthought Canopy: The free and commercial versions include the core SciPy stack packages. Supports Linux, Windows and Mac.
  • Python(x,y): A free distribution including the SciPy stack, based around the Spyder IDE. Windows only.
  • WinPython: A free distribution including the SciPy stack. Windows only.
  • Pyzo: A free distribution based on Python 3 (see Note on Python 3) with the IEP editor. Supports Linux and Windows.
  • Algorete Loopy: A free, community oriented distribution for the SciPy stack maintained by researches at Dartmouth College. Loopy supports both Python 2 and 3 on Linux, Windows and Mac OSX. The distribution is derived from Anaconda with additional packages (e.g. Space Physics, BioInformatics).

我选择使用pyzo,原因是这个词最短,最简单。enthought这个我没找到下载整个包的地方.

 

这是它的网站:http://www.pyzo.org/

这是下载地址:http://www.pyzo.org/downloads.html

http://bitbucket.org/pyzo/pyzo/downloads/pyzo_distro-2013c.win32.exe

但是的联通宽带下载非常慢,qq旋风无法加速,用360云盘离线也出现异常

稍后使用360浏览器下载,网速又到了较可接受的地步。

终于安装好了pyzo!

运行了一下,界面尚可,看着很简单!

IEP编辑器的特点:

IEP
shell 和编辑器
编辑器里,每个文件一个标签,可以右键选择:运行 存盘 关闭等
shell可以同时运行多个,并且可以是不同的python版本

shell配置,5种GUI 工具箱,比如选Matplotlib还是Visvis渲染

Shell >Edit shell configurations 可以修改或添加shell配置

运行代码
run selection可以运行选中的代码
run cell 可以运行两行由‘##’开始的语句中间的所有程序
run file
run progect man file 运行当前工程里的main文件

交互模式和运行脚本模式
交互模式 sys.path[0] 是一个空字符串,sys.argv设置为[”]
脚本模式__file__ 和sys.argv[0] 设置为脚本文件名,sys.path[0]和工作目录设置为包含这个脚本的目录。

便利的工具
在Tools menu 可以选择使用哪些工具
工具系统是被用设计成容易创建你自己的工具的。具体可以看在线wiki或者就看实际的工具例子。

建议的工具
Source structure tool 源文件构建工具
File browser tool 文件浏览工具

ps python3中range就是xrange,所以xrange废弃了。

 

ps,我在360云盘里放置了pyzo的exe安装文件和zip压缩文件,该压缩文件解压后就可以直接使用了

http://yunpan.cn/QpIRFbqx8TNpy

学习递归想到的:有些问题不能用递归解决!

正在学习计算概论里的递归,突然灵光一闪,我以前思考的一个问题,终于找到答案了!

我们知道,汉诺塔可以用递归解决

切蛋糕可以用递归解决

而我本人,一直认为注明的微软面试题:海盗分金子,不能用递归解决! 以前,我的理由是:因为海盗都极度聪明的前提,导致第一个海盗根本就不会做出错误决定,这样也就不会被扔下海,也就不会出现最后只剩2个海盗的情况。竟然不会出现这样的情况,那用不能用倒推(也就是递归)。

上面理由已经很充分了,但没有归纳!

现在归纳下:

所有从1开始,到N的情况(数值,时间/重量/个数等),是可以用递归解决的!如汉诺塔/黄金数列

而只存在N的情况,可能不存在n-1情况的,不可以用递归解决,如海盗分金子,考试日期不可测等问题!

未来可测的那部分,可以用递归,不可测的那部分,不能用递归(倒推)!

利用Python中的matplotlib模块抓取yahoo finance里的历史数据

参见这篇文章:

http://blog.csdn.net/yinyao1992/article/details/8208629

 

from pylab import figure, show
from matplotlib.finance import quotes_historical_yahoo
from matplotlib.dates import YearLocator, MonthLocator, DateFormatter
import datetime
date1 = datetime.date( 2012, 1, 1 )
date2 = datetime.date( 2012, 11, 11 )

daysFmt = DateFormatter(‘%m-%d-%Y’)

quotes = quotes_historical_yahoo(‘MSFT’, date1, date2)

 

也可以一下子抓多个股票的:

symbols = [“CSCO”, “IBM”, “MSFT”]

quotes = [finance.quotes_historical_yahoo(symbol, date1, date2, asobject=True)
for symbol in symbols]

要注意的是:如果多个股票里,有不存在的,可能就会报错,导致程序无法运行(也就是没有raise error语句)

 

用scikits-learn对股票进行聚类分析:

import datetime
import numpy
import sklearn.cluster
from matplotlib import finance

#1. Download price data

# 2011 to 2012
start = datetime.datetime(2011, 01, 01)
end = datetime.datetime(2012, 01, 01)

#Dow Jones symbols
symbols = [“AA”, “CSCO”, “HPQ”, “IBM”, “MSFT”]

quotes = [finance.quotes_historical_yahoo(symbol, start, end, asobject=True)
for symbol in symbols]

close = numpy.array([q.close for q in quotes]).astype(numpy.float)
print close.shape

#2. Calculate affinity matrix
logreturns = numpy.diff(numpy.log(close))
print logreturns.shape

logreturns_norms = numpy.sum(logreturns ** 2, axis=1)
S = – logreturns_norms[:, numpy.newaxis] – logreturns_norms[numpy.newaxis, :] + 2 * numpy.dot(logreturns, logreturns.T)

#3. Cluster using affinity propagation
aff_pro = sklearn.cluster.AffinityPropagation().fit(S)
labels = aff_pro.labels_

for i in xrange(len(labels)):

print ‘%s in Cluster %d’ % (symbols[i], labels[i])

输出:

(5, 252)
(5, 251)
AA in Cluster 0
CSCO in Cluster 2
HPQ in Cluster 1
IBM in Cluster 2
MSFT in Cluster 2
/opt/anaconda/envs/np17py27-1.8/lib/python2.7/site-packages/sklearn/cluster/affinity_propagation_.py:264: UserWarning: The API of AffinityPropagation has changed.Now ``fit`` constructs an affinity matrix from the data. To use a custom affinity matrix, set ``affinity=precomputed``.
  warnings.warn("The API of AffinityPropagation has changed."

Numpy 第10001个素数

import numpy

LIM = 10 ** 6
N = 10 ** 9
P = 10001
primes = []
p = 2

#By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13.

#What is the 10 001st prime number?

def check_primes(a, p):
#2. Sieve out multiples of p
a = a[a % p != 0]

return a

for i in xrange(3, N, LIM):
#1. Create a list of consecutive integers
a = numpy.arange(i, i + LIM, 2)

while len(primes) < P:
a = check_primes(a, p)
primes.append(p)

p = a[0]

print len(primes), primes[P-1]

 

10001 104743