Menu
▾
▴
SF.net SVN: matplotlib:[7976] trunk/py4science/examples/logistic
|
From: <fer...@us...> - 2009-11-21 00:26:27
|
Revision: 7976
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=7976&view=rev
Author: fer_perez
Date: 2009-11-21 00:26:20 +0000 (Sat, 21 Nov 2009)
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.
|
