def perf_comp_data(func_list,data_list,rep=3,number=1): from timeit import repeat res_list={} for name in enumerate(func_list): stmt=name[1]+'('+data_list[name[0]]+')' setup="from __main__ import "+name[1]+', '+data_list[name[0]] results=repeat(stmt=stmt,setup=setup,repeat=rep,number=number) res_list[name[1]]=sum(results)/rep res_sort=sorted(res_list.iteritems(),key=lambda (k,v):(v,k)) for item in res_sort: rel=item[1]/res_sort[0][1] print 'function: ' + item[0] + ', av.time sec:%9.5f, ' % item[1] + 'relative: %6,1f' % rel from math import *def f(x): return abs(cos(x))**0.5+sin(2+3*x) I=500000a_py=range(500000) a_py def f1(a): res=[] for x in a: res.append(f(x)) return res def f2(a): return [f(x) for x in a] def f3(a): ex='abs(cos(x))**0.5+sin(2+3*x)' return [eval(ex) for x in a] import numpy as npa_np=np.arange(I)def f4(a): return (np.abs(np.cos(a))**0.5+np.sin(2+3*a)) import numexpr as ne def f5(a): ex='abs(cos(a))**0.5+sin(2+3*a)' ne.set_num_threads(1) return ne.evaluate(ex) def f6(a): ex='abs(cos(a))**0.5+sin(2+3*a)' ne.set_num_threads(16) return ne.evaluate(ex) %%timer1=f1(a_py)r2=f2(a_py)r3=f3(a_py)r4=f4(a_np)r5=f5(a_np)r6=f6(a_np) np.allclose(r1,r2) np.allclose(r1,r3) np.allclose(r1,r4) np.allclose(r1,r5) np.allclose(r1,r6) func_list=['f1','f2','f3','f4','f5','f6']data_list=['a_py','a_py','a_py','a_np','a_np','a_np'] perf_comp_data(func_list,data_list) def perf_comp_data(func_list,data_list,rep=3,number=1): from timeit import repeat res_list={} for name in enumerate(func_list): stmt=name[1]+'('+data_list[name[0]]+')' setup="from __main__ import "+name[1]+', '+data_list[name[0]] results=repeat(stmt=stmt,setup=setup,repeat=rep,number=number) res_list[name[1]]=sum(results)/rep res_sort=sorted(res_list.iteritems(),key=lambda (k,v):(v,k)) for item in res_sort: rel=item[1]/res_sort[0][1] print 'function: ' + item[0] + ', av.time sec:%9.5f, ' % item[1] + 'relative: %6.1f' % rel perf_comp_data(func_list,data_list) import numpy as np np.zero((3,3),dtype=np.float64,order='C') np.zeros((3,3),dtype=np.float64,order='C') c=np.array([[1.,1.,1.], [2.,2.,2.], [3.,3.,3.]],order='C') f=np.array([[1.,1.,1.], [2.,2.,2.], [3.,3.,3.]],order='F') x=np.random.standard_normal((3,1500000))C=np.array(x,order='C')F=np.array(x,order='F')x=0.0 %timeit C.sum(axis=0) %timeit.C.sum(axis=1) %timeit C.sum(axis=1) %timeit C.std(axis=0) %timeit C.std(axis=1) %timeit F.sum(axis=0) %timeit F.sum(axis=1) %timeit F.std(axis=0) %timeit F.std(axis=1) C=0.0 F=0.0 def bsm_mcs_valuation(strike): import numpy as np S0=100.;T=1.0;r=0.05;vola=0.2 M=50;I=20000 dt=T/M rand=np.random.standard_normal((M+1,I)) S=np.zeros((M+1,I));S[0]=S0 for t in range(1,M+1): S[t]=S[t-1]*np.exp((r-0.5*vola**2)*dt+vola*np.sqrt(dt)*rand[t]) value=(np.exp(-r*T)*np.sum(np.maximum(S[-1]-strike,0))/I) return value def seq_value(n): strikes=np.linspace(80,120,n) option_values=[] for strike in strikes: option_values.append(bsm_mcs_valuation(strike)) return strikes,option_values n=100%time strikes,option_values_seq=seq_value(n) import matplotlib.pyplot as plt%matplotlib inlineplt.figure(figsize=(8,4))plt.plot(strikes,option_values_seq,'b')plt.plot(strikes,option_values_seq,'r.')plt.grid(True)plt.xlabel('strikes')plt.ylabel('European call option values') from IPython.parallel import Clientc=Client(profile="default")view.c.load_balanced_view() from IPython.parallel import Clientc=Client(profile="default")view=c.load_balanced_view() view def par_value(n): strikes=np.linspace(80,120,n) option_values=[] for strike in strikes: value=view.apply_async(bsm_mcs_valuation,strike) option_values.append(value) c.wait(option_values) return strikes,option_values %time strikes,option_values_obj=par_value(n) option_values_obj[0].metadata option_values_obj[0].result option_values_par=[]for res in option_values_obj: option_values_par.append(res.result) plt.figure(figsize=(8,4))plt.plot(strikes,option_values_seq,'b',label='Sequential')plt.plot(strikes,option_values_par,'r.',label='Parallel')plt.grid(True);plt.legend(loc=0)plt.xlabel('strikes')plt.ylabel('European call option values') from IPython.parallel import Client c=Client(profile="default")view=c.load_balanced_view() def par_value(n): strikes=np.linspace(80,120,n) option_values=[] for strike in strikes: value=view.apply_async(bsm_mcs_valuation,strike) option_values.append(value) c.wait(option_values) return strikes,option_values %time strikes,option_values_obj=par_value(n) option_values_obj[0].metadata option_values_obj[0].result option_values_par=[]for res in option_values_obj: option_values_par.append(res.result) plt.figure(figsize=(8,4)) plt.figure(figsize=(8,4))plt.plot(strikes,option_values_seq,'b',label='Sequential')plt.plot(strikes,option_values_par,'r.',label='Parallel')plt.grid(True);plt.legend(loc=0)plt.xlabel('strikes')plt.ylabel('European call option values') n=50func_list=['seq_value','par_value']data_list=2*['n'] perf_comp_data(func_list,data_list) import multiprocessing as mp import mathdef simulate_geometric_brownian_motion(p): M,I=p S0=100.;r=0.05;sigma=0.2;T=1.0 dt=T/M paths=np.zeros((M+1,I)) paths[0]=S0 for t in range(1,M+1): paths[t]=paths[t-1]*np.exp((r-0.5*sigma**2)*dt+sigma*math.sqrt(dt)*np.random.standard_normal(I)) return paths paths=simulate_geometric_brownian_motion((5,2))paths from time import timetimes=[]for w in range(1,17): t0=time() pool=mp.Pool(processes=w) result=pool.map(simulate_geometric_brownian_motion,t*[(M,I),])times.append(time()-t0) I=10000 M=100 t=100 from time import timetimes=[]for w in range(1,17): t0=time() pool=mp.Pool(processes=w) result=pool.map(simulate_geometric_brownian_motion,t*[(M,I),])times.append(time()-t0) ##---(Sat Dec 26 14:27:57 2015)--- import multiprocessing as mp import math def simulate_geometric_brownian_motion(p): M,I=p S0=100.;r=0.05;sigma=0.2;T=1.0 dt=T/M paths=np.zeros((M+1,I)) paths[0]=S0 for t in range(1,M+1): paths[t]=paths[t-1]*np.exp((r-0.5*sigma**2)*dt+sigma*math.sqrt(dt)*np.random.standard_normal(I)) return paths I=10000 M=100 t=100 from time import time times=[] for w in range(1,17): t0=time() pool=mp.Pool(processes=w) result=pool.async(simulate_geometric_brownian_motion,t*[(M,I),]) times.append(time()-t0) from time import time times=[] for w in range(1,17): t0=time() pool=mp.Pool(processes=w) result=pool.apply_async(simulate_geometric_brownian_motion,t*[(M,I),]) times.append(time()-t0) import matplotlib.pyplot as plt %matplotlib inline plt.plot(range(1,17),times) plt.plot(range(1,17),times,'ro') plt.grid(True) plt.xlabel('number of processes') plt.ylabel('time in seconds') plt.title('%d Monte Carlo simulations' % t) times from time import time times=[] for w in range(1,17): t0=time() pool=mp.Pool(processes=w) result=pool.apply_async(simulate_geometric_brownian_motion,t*[(M,I),]) times.append(time()-t0) import matplotlib.pyplot as plt %matplotlib inline plt.plot(range(1,17),times) plt.plot(range(1,17),times,'ro') plt.grid(True) plt.xlabel('number of processes') plt.ylabel('time in seconds') plt.title('%d Monte Carlo simulations' % t) from math import cos,log def f_py(I,J): res=0 for i in range(I): for j in range(J): res+=int(cos(log(1))) return res I,J=5000,5000 %time f_py(I,J) def f_np(I,J): a=np.one((I,J),dtype=np.float64) return int(np.sum(np.cos(np.log(a)))),a %time res,a=f_np(I£¬J) %time res,a=f_np(I,J) import numpy as np %time res,a=f_np(I,J) def f_np(I,J): a=np.ones((I,J),dtype=np.float64) return int(np.sum(np.cos(np.log(a)))),a %time res,a=f_np(I,J) a.nbytes import numba as nb f_nb=nb.jit(f_py) %time f_nb(I,J) func_list=['f_py','f_np','f_nb'] data_list=3*['I,J'] perf_comp_data(func_list,data_list) def perf_comp_data(func_list,data_list,rep=3,number=1): from timeit import repeat res_list={} for name in enumerate(func_list): stmt=name[1]+'('+data_list[name[0]]+')' setup="from __main__ import "+name[1]+', '+data_list[name[0]] results=repeat(stmt=stmt,setup=setup,repeat=rep,number=number) res_list[name[1]]=sum(results)/rep res_sort=sorted(res_list.iteritems(),key=lambda (k,v):(v,k)) for item in res_sort: rel=item[1]/res_sort[0][1] print 'function: ' + item[0] + ', av.time sec:%9.5f, ' % item[1] + 'relative: %6,1f' % rel def perf_comp_data(func_list,data_list,rep=3,number=1): from timeit import repeat res_list={} for name in enumerate(func_list): stmt=name[1]+'('+data_list[name[0]]+')' setup="from __main__ import "+name[i]+', '+data_list[name[0]] results=repeat(stmt=stmt,setup=setup,repeat=rep,number=number) res_list[name[1]]=sum(results)/rep res_sort=sorted(res_list.iteritems(),key=lambda(k,v):(v,k)) for item in res_sort: rel=item[1]/res_sort[0][1] print 'function: ' + item[0] + ', av.time sec:%9.5f, ' % item[1] + 'relative: %6,1f' % rel perf_comp_data(func_list,data_list) def perf_comp_data(func_list,data_list,rep=3,number=1): from timeit import repeat res_list={} for name in enumerate(func_list): stmt=name[1]+'('+data_list[name[0]]+')' setup="from __main__ import "+name[1]+', '+data_list[name[0]] results=repeat(stmt=stmt,setup=setup,repeat=rep,number=number) res_list[name[1]]=sum(results)/rep res_sort=sorted(res_list.iteritems(),key=lambda(k,v):(v,k)) for item in res_sort: rel=item[1]/res_sort[0][1] print 'function: ' + item[0] + ', av.time sec:%9.5f, ' % item[1] + 'relative: %6,1f' % rel perf_comp_data(func_list,data_list) def perf_comp_data(func_list,data_list,rep=3,number=1): from timeit import repeat res_list={} for name in enumerate(func_list): stmt=name[1]+'('+data_list[name[0]]+')' setup="from __main__ import "+name[1]+', '+data_list[name[0]] results=repeat(stmt=stmt,setup=setup,repeat=rep,number=number) res_list[name[1]]=sum(results)/rep res_sort=sorted(res_list.iteritems(),key=lambda(k,v):(v,k)) for item in res_sort: rel=item[1]/res_sort[0][1] print 'function: ' + item[0] + ', av.time sec:%9.5f, ' % item[1] + 'relative: %6.1f' % rel perf_comp_data(func_list,data_list) S0=100. T=1. r=0.05 vola=0.20 M=1000 dt=T/M df=exp(-r*dt) from math import * M=1000 dt=T/M df=exp(-r*dt) u=exp(vola*sqrt(dt)) d=1/u q=(exp(r*dt)-d)/(u-d) import numpy as np def binomial_py(strike): S=np.zeros((M+1,M+1),dtype=np.float64) S[0,0]=S0 z1=0 for j in xrange(1,M+1,1): z1=z1+1 for i in xrange(z1+1): S[i,j]=S[0,0]*(u**j)*(d**(i*2)) def binomial_py(strike): S=np.zeros((M+1,M+1),dtype=np.float64) S[0,0]=S0 z1=0 for j in xrange(1,M+1,1): z1=z1+1 for i in xrange(z1+1): S[i,j]=S[0,0]*(u**j)*(d**(i*2)) iv=np.zeros((M+1,M+1),dtype=np.float64) z2=0 for j in xrange(0,M+1,1): for i in xrange(z2+1): iv[i,j]=max(S[i,j]-strike,0) z2=z2+1 pv=np.zeros((M+1,M+1),dtype=np.float64) pv[:,M]=iv[:,M] z3=M+1 for j in xrange(M-1,-1,-1): z3=z3-1 for i in xrange(z3): pv[i,j]=(q*pv[i,j+1]+(1-q)*pv[i+1,j+1])*df return pv[0,0] %time round(binomial_py(100),3) bsm_mcs_valuation(100) def bsm_mcs_valuation(strike): import numpy as np S0=100.;T=1.0;r=0.05;vola=0.2 M=50;I=20000 dt=T/M rand=np.random.standard_normal((M+1,I)) S=np.zeros((M+1,I));S[0]=S0 for t in range(1,M+1): S[t]=S[t-1]*np.exp((r-0.5*vola**2)*dt+vola*np.sqrt(dt)*rand[t]) value=(np.exp(-r*T)*np.sum(np.maximum(S[-1]-strike,0))/I) return value %time round(bsm_mcs_valuation(100),3) def binomial_np(strike): mu=np.arange(M+1) mu=np.resize(mu,(M+1,M+1)) md=np.transpose(mu) mu=u**(mu-md) md=d**md S=S0*mu*md pv=np.maximum(S-strike,0) z=0 for t in range(M-1,-1,-1): pv[0:M-z,t]=(q*pv[0:M-z,t+1]+(1-q)*pv[1:M-z+1,t+1])*df z+=1 return pv[0,0] M=4 mu=np.arange(M+1) mu mu=np.resize(mu,(M+1,M+1)) mu md=np.transpose(mu) md mu=u**(mu-md) mu.round(3) md=d**md md.round(3) S=S0*mu*md S M=1000 %time round(binomial_np(100),3) binomial_nb=nb.jit(binomial_py) %time round(binomial_nb(100),3) func_list=['binomial_py','binomial_np','binomial_nb'] k=100. data_list=3*['k'] perf_comp_data(func_list,data_list) def f_py(I,J): res=0 for i in range(I): for j in range(J*I): res+=1 return res I,J=500,500 %time f_py(I,J) def f_cy(int I,int J): cdef double res=0 for i in range(I): for j in range(J*I): res+=1 return res import pyximport pyximport.install() import sys sys.path.append('C:/Users/zhuto/.spyder2/') from nested_loop import f_cy %time res=f_cy(I,J) res %load_ext cythonmagic %load_ext Cython %%cython def f_cy(int I,int J): cdef double res=0 for i in range(I): for j in range(J*I): res+=1 return res %time res=f_cy(I,J) res import numba as nb f_nb=nb.jit(f_py) %time res=f_nb(I,J) res func_list=['f_py','f_cy','f_nb'] I,J=500,500 data_list=3*['I,J'] perf_comp_data(func_list,data_list) from numbapro.cudalib import curand def get_randoms(x,y): rand=np.random.standard_normal((x,y)) return rand get_randoms(2,2) def get_cuda_randoms(x,y): rand=np.empty((x*y),np.float64) prng=curand.PRNG(rndtype=curand.PRNG.XORWOW) prng.normal(rand,0,1) rand=rand.reshpe((x,y)) return rand get_cuda_randoms(2,2) def get_cuda_randoms(x,y): rand=np.empty((x*y),np.float64) prng=curand.PRNG(rndtype=curand.PRNG.XORWOW) prng.normal(rand,0,1) rand=rand.reshape((x,y)) return rand get_cuda_randoms(2,2) %timeit a=get_randoms(1000,1000) %timeit a=get_cuda_randoms(1000,1000) import time as t step=1000 def time_comparision(factor): cuda_times=list() cpu_times=list() for j in range(1,10002,step): i=j*factor t0=t.time() a=get_randoms(i,1) t1=t.time() cpu_times.append(t1-t0) t2=t.time() a.get_cuda_randoms(i,1) t3=t.time() cuda_times.append(t3-t2) print "Bytes of largest array %i" %a.nbytes return cuda_times,cpu_times def plot_results(cpu_times,cuda_times,factor): plt.plot(x*factor,cpu_times,'b',label='NUMPY') plt.plot(x*factor,cuda_times,'r',label='CUDA') plt.legend(loc=0) plt.grid(True) plt.xlabel('size of random number array') plt.ylabel('time') plt.axis('tight') factor=100 cuda_times,cpu_times=time_comparision(factor) import time as t step=1000 def time_comparision(factor): cuda_times=list() cpu_times=list() for j in range(1,10002,step): i=j*factor t0=t.time() a=get_randoms(i,1) t1=t.time() cpu_times.append(t1-t0) t2=t.time() a=get_cuda_randoms(i,1) t3=t.time() cuda_times.append(t3-t2) print "Bytes of largest array %i" %a.nbytes return cuda_times,cpu_times factor=100 cuda_times,cpu_times=time_comparision(factor) x=np.arange(1,10002,step) plot_results(cpu_times,cuda_times,factor) factor=10 cuda_times,cpu_times=time_comparision(factor) plot_results(cpu_times,cuda_times,factor) factor=5000 cuda_times,cpu_times=time_comparision(factor) plot_results(cpu_times,cuda_times,factor)