Revision: 7976 http://matplotlib.svn.sourceforge.net/matplotlib/?rev=7976&view=rev Author: fer_perez Date: 2009年11月21日 00:26:20 +0000 (2009年11月21日) Log Message: ----------- Updated to modern numpy/mpl code, cleaned up interface Modified Paths: -------------- trunk/py4science/examples/logistic/exercise01.py trunk/py4science/examples/logistic/exercise02.py trunk/py4science/examples/logistic/maplib.py trunk/py4science/examples/logistic/maplib.pyc Modified: trunk/py4science/examples/logistic/exercise01.py =================================================================== --- trunk/py4science/examples/logistic/exercise01.py 2009年11月20日 19:08:18 UTC (rev 7975) +++ trunk/py4science/examples/logistic/exercise01.py 2009年11月21日 00:26:20 UTC (rev 7976) @@ -1,29 +1,28 @@ +import numpy as np +import matplotlib.pyplot as plt + from maplib import Logistic -from matplotlib.numerix import arange, absolute, log, Float -from matplotlib.mlab import polyfit, polyval -from pylab import subplot, plot, show epsilon = 1e-10 x0 = 0.4 y0 = x0 + epsilon logmap = Logistic(0.9) -x = logmap.iterate(x0, 100) -y = logmap.iterate(y0, 100) -ind = arange(len(x), typecode=Float) +x = logmap.trajectory(x0, 100) +y = logmap.trajectory(y0, 100) +ind = np.arange(len(x), dtype=float) # x-y \sim epsilon exp(lambda * t) # log(|x-y|) = log(epsilon) + lambda*t (b=log(epsilon) and m=lambda) -d = log(absolute(x-y)) -coeffs = polyfit(ind, d, 1) -lambda_, b = coeffs -print 'lyapunov exponent= %1.3f'%lambda_ -print 'log(epsilon)=%1.3f, b = %1.3f' %(log(epsilon), b) +d = np.log(abs(x-y)) +coeffs = np.polyfit(ind, d, 1) +lyap_exp, b = coeffs +print 'lyapunov exponent= %1.3f' % lyap_exp +print 'log(epsilon)=%1.3f, b = %1.3f' % (np.log(epsilon), b) # now plot the result -subplot(211) -plot(ind, x, '-r', ind, x, '--g', ) -subplot(212) -plot(ind, d, ind, polyval(coeffs, ind)) -show() - +plt.subplot(211) +plt.plot(ind, x, '-r', ind, x, '--g', ) +plt.subplot(212) +plt.plot(ind, d, ind, np.polyval(coeffs, ind)) +plt.show() Modified: trunk/py4science/examples/logistic/exercise02.py =================================================================== --- trunk/py4science/examples/logistic/exercise02.py 2009年11月20日 19:08:18 UTC (rev 7975) +++ trunk/py4science/examples/logistic/exercise02.py 2009年11月21日 00:26:20 UTC (rev 7976) @@ -1,39 +1,46 @@ -from maplib import Logistic,Sine -from matplotlib.numerix import arange, pi, sqrt, sort, zeros,ones, Float -from matplotlib.numerix.mlab import rand -from pylab import figure, draw, show, ion, ioff,frange,gcf,rcParams,rc +from maplib import Logistic, Sine -def bifurcation_diagram(map_type,param0=0,param1=1,nparam=300, - ntransients=100,ncycles=200, dotcolor="0.5", +import numpy as np +import matplotlib +import matplotlib.pyplot as plt + +def bifurcation_diagram(map_type, params, + ntransients=100, ncycles=200, dotcolor="0.5", fig=None, - nboundaries=0): + nboundaries=2): + """Plot a bifurcation diagram for an iterated map object. - if nboundaries: - boundaries = zeros((nparam,nboundaries),Float) - params = frange(param0,param1,npts=nparam) + Parameters + ---------- + + map_type : functor + An iterated map constructor. + """ + nparam = len(params) + if nboundaries>0: + boundaries = np.zeros((nparam, nboundaries)) bound_rng = range(nboundaries) xs = [] ys = [] for idx,param in enumerate(params): m = map_type(param) - y = m.iterate(m.iterate(0.5,ntransients,lastonly=True), - ncycles, lastonly=False) - xs.extend(param*ones(len(y))) + y = m.trajectory(m.iterate_from(0.5, ntransients), ncycles) + xs.extend(param*np.ones_like(y)) ys.extend(y) if nboundaries: # the boundaries are the iterates of the map's maximum, we assume # here that it's located at 0.5 - boundaries[idx] = m.iterate(0.5,nboundaries,lastonly=False)[1:] + boundaries[idx] = m.trajectory(0.5, nboundaries+1)[1:] + + # Make final figure if fig is None: - fig = figure() + fig = plt.figure() ax = fig.add_subplot(111) - # save state (John, is there a cleaner way to do this?) - ax.plot(xs, ys, '.', mfc=dotcolor, mec=dotcolor, ms=1, mew=0) + ax.set_xlim(params[0], params[-1]) - bound_lines = [] - for i in bound_rng: - bound_lines.extend(ax.plot(params,boundaries[:,i])) + if nboundaries>0: + bound_lines = ax.plot(params, boundaries) def toggle(event): if event.key!='t': return @@ -46,25 +53,27 @@ fig.canvas.mpl_connect('key_press_event', toggle) return fig + def cobweb(mu, walkers=10, steps=7): - f = figure() + f = plt.figure() ax = f.add_subplot(111) - interval = frange(0.0, 1.0, npts=100) + interval = np.linspace(0.0, 1.0, 100) logmap = Logistic(mu) logmap.plot(ax, interval, lw=2) - for x0 in rand(walkers): + for x0 in np.random.rand(walkers): logmap.plot_cobweb(ax, x0, steps, lw=2) ax.set_title('Ex 2A: Random init. cond. for mu=%1.3f'%mu) return f + def invariant_density(mu, x0,cycles=1000000,ret_all=False): transients = 1000 bins = 500 - f = figure() + f = plt.figure() ax = f.add_subplot(111) logmap = Logistic(mu) - y = logmap.iterate(x0, transients, lastonly=True) - y = logmap.iterate(y, cycles, lastonly=False) + y0 = logmap.iterate_from(x0, transients) + y = logmap.trajectory(y0, cycles) n, bins, patches = ax.hist(y, bins, normed=1) ax.set_xlim(0,1) if ret_all: @@ -82,14 +91,14 @@ def ex2B(): def rho(x): - return 1./(pi * sqrt( x*(1.-x))) + return 1./(np.pi * np.sqrt( x*(1.-x))) # Don't start from 0.5, which is a fixed point! f = invariant_density(1.0,0.567) ax = f.gca() # avoid the edges: rho(x) is singular at 0 and 1! - x0 = frange(0.001, 0.999, npts=1000) - l, = ax.plot(x0, rho(x0), 'r-', lw=3, alpha=0.5) + x0 = np.linspace(0.001, 0.999, 1000) + l, = ax.plot(x0, rho(x0), 'r-', lw=3, alpha=0.7) ax.set_title('Ex 2B: invariant density') ax.legend((ax.patches[0], l), ('empirical', 'analytic'), loc='upper center') ax.set_xlim(0,1) @@ -98,27 +107,26 @@ def ex2CD(mu=0.9,x0=0.64): - f,logmap,n = invariant_density(mu,x0,ret_all=True) - ax = f.gca() - ax.set_xticks(frange(0,1,npts=10)) + fig, logmap, n = invariant_density(mu, x0, ret_all=True) + ax = fig.gca() + ax.set_xticks(np.linspace(0, 1, 10)) ax.grid(True) # Now, iterate x=1/2 a few times and plot this 'orbit', which corresponds # to peaks in the invariant density. x0 = 0.5 - pts = logmap.iterate(x0,10) - pts_y = 0.5*frange(1,max(n),npts=len(pts)) + pts = logmap.trajectory(x0,10) + pts_y = 0.5*np.linspace(1, max(n), len(pts)) ax.plot(pts,pts_y,'ro-') - ax.set_title('Ex 2C/D: Analytics of cusps at mu=%0.2g' % mu) + ax.set_title('**Ex 2C/D: Analytics of cusps at mu=%0.2g' % mu) def ex2E(): - par0 = 0.8 - par1 = 1.0 - fig = bifurcation_diagram(Logistic,par0,par1,nparam=1000, - nboundaries=8) + # Parameter grid to sample each map on + params = np.linspace(0.5, 1, 500) + fig = bifurcation_diagram(Logistic, params, nboundaries=8) fig.gca().set_title('Ex 2E: Bifurcation diag. with boundaries (press t)') - fig = bifurcation_diagram(Logistic,par0,par1,dotcolor='blue') - fig = bifurcation_diagram(Sine,par0,par1,dotcolor='red',fig = fig) + fig = bifurcation_diagram(Logistic, params, dotcolor='blue') + fig = bifurcation_diagram(Sine, params, dotcolor='red', fig = fig) fig.gca().set_title('Ex 2E: Bifurcation diag. logistic/sine maps') @@ -127,4 +135,4 @@ ex2B() ex2CD() ex2E() - show() + plt.show() Modified: trunk/py4science/examples/logistic/maplib.py =================================================================== --- trunk/py4science/examples/logistic/maplib.py 2009年11月20日 19:08:18 UTC (rev 7975) +++ trunk/py4science/examples/logistic/maplib.py 2009年11月21日 00:26:20 UTC (rev 7976) @@ -1,8 +1,8 @@ -import matplotlib.numerix as nx -from matplotlib.mlab import polyfit +import numpy as np + from matplotlib.cbook import iterable -class SomeMap: +class IteratedMap(object): """ Define an interface for a map """ @@ -43,7 +43,7 @@ kwargs are passed onto mpl plot """ - iterates = self.iterate(x0, numsteps) + iterates = self.trajectory(x0, numsteps) vertices = [] lasty = 0 for this, next in zip(iterates[:-1], iterates[1:]): @@ -55,65 +55,66 @@ x, y = zip(*vertices) ax.plot(x, y, **kwargs) - def iterate(self, x0, numsteps, lastonly=False): - """ - iterate self starting at x0 for numsteps + def iterator_from(self, x0): + while 1: + x0 = self(x0) + yield x0 + + + def iterate_from(self, x0, numsteps): + for i in xrange(numsteps): + x0 = self(x0) + return x0 + + + def trajectory(self, x0, numsteps): + """iterate self starting at x0 for numsteps, returning whole trajectory + Return value is an array of the time-series. If x0 is a scalar, a - numsteps+1 length 1D vector is returned with x0 as the first - value. If x0 is a vector, an numsteps+1 x len(x0) 2D array is - returned with x0 as the first row - - if lastonly is True, only return the last iterate + numsteps length 1D vector is returned with x0 as the first value. If x0 + is a vector, a 2D array with shape (numsteps, len(x0)) is returned, with + x0 as the first row. """ - if not lastonly: - if iterable(x0): # return a 2D array - ret = nx.zeros( (numsteps+1,len(x0)), typecode=nx.Float ) - else: # return a 1D array - ret = nx.zeros( (numsteps+1,), typecode=nx.Float ) + if iterable(x0): # return a 2D array + ret = np.zeros( (numsteps,len(x0)) ) + else: # return a 1D array + ret = np.zeros(numsteps) # assign the initial condtion to the 0-th element - if not lastonly: ret[0] = x0 - + ret[0] = x0 # iterate the map for numsteps - last = x0 - for i in range(1,numsteps+1): - this = self(last) - if not lastonly: ret[i] = this - last = this - - if lastonly: - return last - else: - return ret + for i in range(1,numsteps): + ret[i] = self(ret[i-1]) + return ret -class Logistic(SomeMap): +class Logistic(IteratedMap): def __init__(self, mu): - self.R = 4.*mu + self.R = 4.0*mu def __call__(self, x): 'iterate self one step starting at x. x can be a scalar or array' - return self.R*x*(1.-x) + return self.R*x*(1.0-x) -class Sine(SomeMap): +class Sine(IteratedMap): def __init__(self, B): self.B = B def __call__(self, x): 'iterate self one step starting at x. x can be a scalar or array' - return self.B*nx.sin(nx.pi*x) + return self.B*np.sin(np.pi*x) def test(): m = Logistic(0.9) - x0 = nx.mlab.rand(100) + x0 = np.random.rand(100) ret = m.iterate(x0, 3) - assert( ret.shape == 4,100) + assert ret.shape == 4,100 x0 = 0.2 ret = m.iterate(x0, 3) - assert( ret.shape == 4,) + assert ret.shape == 4 print 'all tests passed' Modified: trunk/py4science/examples/logistic/maplib.pyc =================================================================== (Binary files differ) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.