From: <mur...@us...> - 2011-07-26 02:41:10
|
Revision: 170 http://python-control.svn.sourceforge.net/python-control/?rev=170&view=rev Author: murrayrm Date: 2011-07-26 02:41:03 +0000 (Tue, 26 Jul 2011) Log Message: ----------- * Updated phaseplot module with new (more pythonic) call signature * Added boxplot function to generate initial conditions around a box * New genetic switch example (from FBS), demonstrating phase plots Modified Paths: -------------- trunk/ChangeLog trunk/examples/phaseplots.py trunk/src/__init__.py trunk/src/phaseplot.py Added Paths: ----------- trunk/examples/genswitch.py Modified: trunk/ChangeLog =================================================================== --- trunk/ChangeLog 2011-07-25 05:00:41 UTC (rev 169) +++ trunk/ChangeLog 2011-07-26 02:41:03 UTC (rev 170) @@ -1,3 +1,17 @@ +2011-07-25 Richard Murray <murray@malabar.local> + + * examples/phaseplots.py: updated calls to PhasePlot to use new + argument structure + + * src/phaseplot.py (PhasePlot): Updated call signature to be + more pythonic and fixed up documentation. + + * examples/genswitch.py (genswitch): added new example showing + PhasePlot functionality + + * src/phaseplot.py (boxgrid): added function to compute initial + conditions around the edges of a box + 2011-07-24 Richard Murray <murray@malabar.local> * tests/margin_test.py: added simple unit tests for margin functions Added: trunk/examples/genswitch.py =================================================================== --- trunk/examples/genswitch.py (rev 0) +++ trunk/examples/genswitch.py 2011-07-26 02:41:03 UTC (rev 170) @@ -0,0 +1,81 @@ +# genswitch_plot.m - run the Collins genetic switch model +# RMM, 24 Jan 07 +# +# This file contains an example from FBS of a simple dynamical model +# of a genetic switch. Plots time traces and a phase portrait using +# the python-control library. + +import numpy as np +import matplotlib.pyplot as mpl +from scipy.integrate import odeint +from matplotlib.mlab import frange +from control import PhasePlot, boxgrid + +# Simple model of a genetic switch +# +# This function implements the basic model of the genetic switch +# Parameters taken from Gardner, Cantor and Collins, Nature, 2000 +def genswitch(y, t, mu=4, n=2): + return (mu / (1 + y[1]**n) - y[0], mu / (1 + y[0]**n) - y[1]) + +# Run a simulation from an initial condition +tim1 = np.linspace(0, 10, 100) +sol1 = odeint(genswitch, [1, 5], tim1) + +# Extract the equlibirum points +mu = 4; n = 2; # switch parameters +eqpt = np.empty(3); +eqpt[0] = sol1[0,-1] +eqpt[1] = sol1[1,-1] +eqpt[2] = 0; # fzero(@(x) mu/(1+x^2) - x, 2); + +# Run another simulation showing switching behavior +tim2 = np.linspace(11, 25, 100); +sol2 = odeint(genswitch, sol1[-1,:] + [2, -2], tim2) + +# First plot out the curves that define the equilibria +u = frange(0, 4.5, 0.1) +f = np.divide(mu, (1 + u**n)) # mu / (1 + u^n), elementwise + +mpl.figure(1); mpl.clf(); +mpl.axis([0, 5, 0, 5]); # box on; +mpl.plot(u, f, '-', f, u, '--') # 'LineWidth', AM_data_linewidth); +mpl.legend('z1, f(z1)', 'z2, f(z2)') # legend(lgh, 'boxoff'); +mpl.plot([0, 3], [0, 3], 'k-') # 'LineWidth', AM_ref_linewidth); +mpl.plot(eqpt[0], eqpt[1], 'k.', eqpt[1], eqpt[0], 'k.', + eqpt[2], eqpt[2], 'k.') # 'MarkerSize', AM_data_markersize*3); +mpl.xlabel('z1, f(z2)'); +mpl.ylabel('z2, f(z1)'); + +# Time traces +mpl.figure(3); mpl.clf(); # subplot(221); +mpl.plot(tim1, sol1[:,0], 'b-', tim1, sol1[:,1], 'g--'); +# set(pl, 'LineWidth', AM_data_linewidth); +mpl.hold(True); +mpl.plot([tim1[-1], tim1[-1]+1], + [sol1[-1,0], sol2[0,1]], 'ko:', + [tim1[-1], tim1[-1]+1], [sol1[-1,1], sol2[0,0]], 'ko:'); +# set(pl, 'LineWidth', AM_data_linewidth, 'MarkerSize', AM_data_markersize); +mpl.plot(tim2, sol2[:,0], 'b-', tim2, sol2[:,1], 'g--'); +# set(pl, 'LineWidth', AM_data_linewidth); +mpl.axis([0, 25, 0, 5]); + +mpl.xlabel('Time {\itt} [scaled]'); +mpl.ylabel('Protein concentrations [scaled]'); +mpl.legend('z1 (A)', 'z2 (B)') # 'Orientation', 'horizontal'); +# legend(legh, 'boxoff'); + +# Phase portrait +mpl.figure(2); mpl.clf(); # subplot(221); +mpl.axis([0, 5, 0, 5]); # set(gca, 'DataAspectRatio', [1, 1, 1]); +PhasePlot(genswitch, X0 = boxgrid([0, 5, 6], [0, 5, 6]), T = 10, + timepts = [0.2, 0.6, 1.2]) + +# Add the stable equilibrium points +mpl.hold(True); +mpl.plot(eqpt[0], eqpt[1], 'k.', eqpt[1], eqpt[0], 'k.', + eqpt[2], eqpt[2], 'k.') # 'MarkerSize', AM_data_markersize*3); + +mpl.xlabel('Protein A [scaled]'); +mpl.ylabel('Protein B [scaled]'); # 'Rotation', 90); + Modified: trunk/examples/phaseplots.py =================================================================== --- trunk/examples/phaseplots.py 2011-07-25 05:00:41 UTC (rev 169) +++ trunk/examples/phaseplots.py 2011-07-26 02:41:03 UTC (rev 170) @@ -27,16 +27,16 @@ # Outer trajectories PhasePlot(invpend_ode, - 'logtime', (3, 0.7), None, - [ [-2*pi, 1.6], [-2*pi, 0.5], [-1.8, 2.1], - [-1, 2.1], [4.2, 2.1], [5, 2.1], - [2*pi, -1.6], [2*pi, -0.5], [1.8, -2.1], - [1, -2.1], [-4.2, -2.1], [-5, -2.1] ], - np.linspace(0, 40, 200)) + X0 = [ [-2*pi, 1.6], [-2*pi, 0.5], [-1.8, 2.1], + [-1, 2.1], [4.2, 2.1], [5, 2.1], + [2*pi, -1.6], [2*pi, -0.5], [1.8, -2.1], + [1, -2.1], [-4.2, -2.1], [-5, -2.1] ], + T = np.linspace(0, 40, 200), + logtime = (3, 0.7) ) # Separatrices mpl.hold(True); -PhasePlot(invpend_ode, 'auto', 0, None, [[-2.3056, 2.1], [2.3056, -2.1]], 6) +PhasePlot(invpend_ode, X0 = [[-2.3056, 2.1], [2.3056, -2.1]], T=6, lingrid=0) mpl.show(); # @@ -57,10 +57,10 @@ mpl.figure(); mpl.clf(); mpl.axis([-1, 1, -1, 1]); # set(gca, 'DataAspectRatio', [1, 1, 1]); PhasePlot(oscillator_ode, - 'timepts', [0.25, 0.8, 2, 3], None, [ - [-1, 1], [-0.3, 1], [0, 1], [0.25, 1], [0.5, 1], [0.75, 1], [1, 1], - [1, -1], [0.3, -1], [0, -1], [-0.25, -1], [-0.5, -1], [-0.75, -1], [-1, -1] - ], np.linspace(0, 8, 80)) + X0 = [ + [-1, 1], [-0.3, 1], [0, 1], [0.25, 1], [0.5, 1], [0.75, 1], [1, 1], + [1, -1], [0.3, -1], [0, -1], [-0.25, -1], [-0.5, -1], [-0.75, -1], [-1, -1] + ], T = np.linspace(0, 8, 80), timepts = [0.25, 0.8, 2, 3]) mpl.hold(True); mpl.plot([0], [0], 'k.'); # 'MarkerSize', AM_data_markersize*3); # set(gca,'DataAspectRatio',[1,1,1]); mpl.xlabel('x1'); mpl.ylabel('x2'); @@ -81,11 +81,13 @@ m = 1; b = 1; k = 1; # default values mpl.figure(); mpl.clf(); mpl.axis([-1, 1, -1, 1]); # set(gca, 'DataAspectRatio', [1 1 1]); -PhasePlot(oscillator_ode, 'timepts', [0.3, 1, 2, 3], None, - [[-1,1], [-0.3,1], [0,1], [0.25,1], [0.5,1], [0.7,1], [1,1], [1.3,1], - [1,-1], [0.3,-1], [0,-1], [-0.25,-1], [-0.5,-1], [-0.7,-1], [-1,-1], - [-1.3,-1]], - np.linspace(0, 10, 100), parms = (m, b, k)); +PhasePlot(oscillator_ode, + X0 = [ + [-1,1], [-0.3,1], [0,1], [0.25,1], [0.5,1], [0.7,1], [1,1], [1.3,1], + [1,-1], [0.3,-1], [0,-1], [-0.25,-1], [-0.5,-1], [-0.7,-1], [-1,-1], + [-1.3,-1] + ], T = np.linspace(0, 10, 100), + timepts = [0.3, 1, 2, 3], parms = (m, b, k)); mpl.hold(True); mpl.plot([0], [0], 'k.'); # 'MarkerSize', AM_data_markersize*3); # set(gca,'FontSize', 16); mpl.xlabel('{\itx}_1'); mpl.ylabel('{\itx}_2'); @@ -93,14 +95,14 @@ # Saddle mpl.figure(); mpl.clf(); mpl.axis([-1, 1, -1, 1]); # set(gca, 'DataAspectRatio', [1 1 1]); -PhasePlot(saddle_ode, 'timepts', [0.2, 0.5, 0.8], None, +PhasePlot(saddle_ode, scale = 2, timepts = [0.2, 0.5, 0.8], X0 = [ [-1, -1], [1, 1], [-1, -0.95], [-1, -0.9], [-1, -0.8], [-1, -0.6], [-1, -0.4], [-1, -0.2], [-0.95, -1], [-0.9, -1], [-0.8, -1], [-0.6, -1], [-0.4, -1], [-0.2, -1], [1, 0.95], [1, 0.9], [1, 0.8], [1, 0.6], [1, 0.4], [1, 0.2], [0.95, 1], [0.9, 1], [0.8, 1], [0.6, 1], [0.4, 1], [0.2, 1], [-0.5, -0.45], [-0.45, -0.5], [0.5, 0.45], [0.45, 0.5], - [-0.04, 0.04], [0.04, -0.04] ], np.linspace(0, 2, 20)); + [-0.04, 0.04], [0.04, -0.04] ], T = np.linspace(0, 2, 20)); mpl.hold(True); mpl.plot([0], [0], 'k.'); # 'MarkerSize', AM_data_markersize*3); # set(gca,'FontSize', 16); mpl.xlabel('{\itx}_1'); mpl.ylabel('{\itx}_2'); @@ -109,10 +111,10 @@ m = 1; b = 0; k = 1; # zero damping mpl.figure(); mpl.clf(); mpl.axis([-1, 1, -1, 1]); # set(gca, 'DataAspectRatio', [1 1 1]); -PhasePlot(oscillator_ode, 'timepts', - [pi/6, pi/3, pi/2, 2*pi/3, 5*pi/6, pi, 7*pi/6, 4*pi/3, 9*pi/6, 5*pi/3, 11*pi/6, 2*pi], None, - [ [0.2,0], [0.4,0], [0.6,0], [0.8,0], [1,0], [1.2,0], [1.4,0] ], - np.linspace(0, 20, 200), parms = (m, b, k)); +PhasePlot(oscillator_ode, timepts = + [pi/6, pi/3, pi/2, 2*pi/3, 5*pi/6, pi, 7*pi/6, 4*pi/3, 9*pi/6, 5*pi/3, 11*pi/6, 2*pi], + X0 = [ [0.2,0], [0.4,0], [0.6,0], [0.8,0], [1,0], [1.2,0], [1.4,0] ], + T = np.linspace(0, 20, 200), parms = (m, b, k)); mpl.hold(True); mpl.plot([0], [0], 'k.') # 'MarkerSize', AM_data_markersize*3); # set(gca,'FontSize', 16); mpl.xlabel('{\itx}_1'); mpl.ylabel('{\itx}_2'); Modified: trunk/src/__init__.py =================================================================== --- trunk/src/__init__.py 2011-07-25 05:00:41 UTC (rev 169) +++ trunk/src/__init__.py 2011-07-26 02:41:03 UTC (rev 170) @@ -64,6 +64,7 @@ from mateqn import * from modelsimp import * from nichols import * +from phaseplot import PhasePlot, boxgrid from rlocus import * from statefbk import * from statesp import * Modified: trunk/src/phaseplot.py =================================================================== --- trunk/src/phaseplot.py 2011-07-25 05:00:41 UTC (rev 169) +++ trunk/src/phaseplot.py 2011-07-26 02:41:03 UTC (rev 170) @@ -1,7 +1,8 @@ # phaseplot.py - generate 2D phase portraits # # Author: Richard M. Murray -# Date: Fall 2002 (MATLAB version), based on a version by Kristi Morgansen +# Date: 24 July 2011, converted from MATLAB version (2002); based on +# a version by Kristi Morgansen # # Copyright (c) 2011 by California Institute of Technology # All rights reserved. @@ -38,99 +39,107 @@ from exception import ControlNotImplemented from scipy.integrate import odeint -def PhasePlot(odefun, xlimv, ylimv, scale=1, xinit=None, T=None, parms=(), - verbose=True): +def PhasePlot(odefun, X=None, Y=None, scale=1, X0=None, T=None, + lingrid=None, lintime=None, logtime=None, timepts=None, + parms=(), verbose=True): """ - PhasePlot Phase plot for 2D dynamical systems + Phase plot for 2D dynamical systems - PhasePlot(F, X1range, X2range, scale) produces a quiver plot for - the function F. X1range and X2range have the form [min, max, num] - and specify the axis limits for the plot, along with the number of - subdivisions in each axis, used to produce an quiver plot for the - vector field. The vector field is scaled by a factor 'scale' - (default = 1). + Produces a vector field or stream line plot for a planar system. - The function F should be the same for as used for scipy.integrate. - Namely, it should be a function of the form dxdt = F(x, t) that - accepts a state x of dimension 2 and returns a derivative dxdt of - dimension 2. + Call signatures: + PhasePlot(func, X, Y, ...) - display vector field on meshgrid + PhasePlot(func, X, Y, scale, ...) - scale arrows + PhasePlot(func. X0=(...), T=Tmax, ...) - display stream lines + PhasePlot(func, X, Y, X0=[...], T=Tmax, ...) - plot both + PhasePlot(func, X0=[...], T=Tmax, lingrid=N, ...) - plot both + PhasePlot(func, X0=[...], lintime=N, ...) - stream lines with arrows - PhasePlot(F, X1range, X2range, scale, Xinit) produces a phase - plot for the function F, consisting of the quiver plot plus stream - lines. The streamlines start at the initial conditions listed in - Xinit, which should be a matrix whose rows give the desired inital - conditions for x1 and x2. X1range or X2range is 'auto', the arrows - are produced based on the stream lines. If 'scale' is negative, - dots are placed at the base of the arrows. If 'scale' is zero, no - dots are produced. - - PhasePlot(F, X1range, X2range, scale, boxgrid(X1range2, X2range2)) - produces a phase plot with stream lines generated at the edges of - the rectangle defined by X1range2, X2range2. These ranges are in - the same form as X1range, X2range. + Parameters + ---------- + func : callable(x, t, ...) + Computes the time derivative of y (compatible with odeint). + The function should be the same for as used for + scipy.integrate. Namely, it should be a function of the form + dxdt = F(x, t) that accepts a state x of dimension 2 and + returns a derivative dx/dt of dimension 2. - PhasePlot(F, X1range, X2range, scale, Xinit, T) produces a phase - plot where the streamlines are simluted for time T (default = 50). + X, Y: ndarray, optional + Two 1-D arrays representing x and y coordinates of a grid. + These arguments are passed to meshgrid and generate the lists + of points at which the vector field is plotted. If absent (or + None), the vector field is not plotted. - PhasePlot(F, X1range, X2range, scale, Xinit, T, P1, P2, ...) - passes additional parameters to the function F, in the same way as - ODE45. + scale: float, optional + Scale size of arrows; default = 1 - Instead of drawing arrows on a grid, arrows can also be drawn on - streamlines by usinge the X1range and X2range arguments as follows: + X0: ndarray of initial conditions, optional + List of initial conditions from which streamlines are plotted. + Each initial condition should be a pair of numbers. - X1range X2range - ------- ------- - 'auto' N Draw N arrows using equally space time points - 'logtime' [N, lam] Draw N arrows using exponential time constant lam - 'timepts' [t1, t2, ...] Draw arrows at the list times - """ + T: array-like or number, optional + Length of time to run simulations that generate streamlines. + If a single number, the same simulation time is used for all + initial conditions. Otherwise, should be a list of length + len(X0) that gives the simulation time for each initial + condition. Default value = 50. - # - # Parameters defining the plot - # - # The constants below define the parameters that control how the - # plot looks. These can be modified to customize the look of the - # phase plot. - # + lingrid = N or (N, M): integer or 2-tuple of integers, optional + If X0 is given and X, Y are missing, a grid of arrows is + produced using the limits of the initial conditions, with N + grid points in each dimension or N grid points in x and M grid + points in y. - #! TODO: convert to keywords - #! TODO: eliminate old parameters that aren't used - # PP color = ['m', 'c', 'r', 'g', 'b', 'k', 'y']; - PP_stream_color = ('b'); # color array for streamlines - PP_stream_linewidth = 1; # line width for stream lines + lintime = N: integer, optional + Draw N arrows using equally space time points - PP_arrow_linewidth = 1; # line width for arrows (quiver) - PP_arrow_markersize = 10; # size of arrow base marker + logtime = (N, lambda): (integer, float), optional + Draw N arrows using exponential time constant lambda + timepts = [t1, t2, ...]: array-like, optional + Draw arrows at the given list times + + parms: tuple, optional + List of parameters to pass to vector field: func(x, t, *parms) + + See also + -------- + boxgrid(X, Y): construct box-shaped grid of initial conditions + + Examples + -------- + """ + # # Figure out ranges for phase plot (argument processing) # - auto = 0; logtime = 0; timepts = 0; Narrows = 0; - if (isinstance(xlimv, str) and xlimv == 'auto'): - auto = 1; - Narrows = ylimv; + #! TODO: need to add error checking to arguments + autoFlag = False; logtimeFlag = False; timeptsFlag = False; Narrows = 0; + if (lingrid != None): + autoFlag = True; + Narrows = lingrid; if (verbose): print 'Using auto arrows\n'; - elif (isinstance(xlimv, str) and xlimv == 'logtime'): - logtime = 1; - Narrows = ylimv[0]; - timefactor = ylimv[1]; + elif (logtime != None): + logtimeFlag = True; + Narrows = logtime[0]; + timefactor = logtime[1]; if (verbose): print 'Using logtime arrows\n'; - elif (isinstance(xlimv, str) and xlimv == 'timepts'): - timepts = 1; - Narrows = len(ylimv); + elif (timepts != None): + timeptsFlag = True; + Narrows = len(timepts); else: # Figure out the set of points for the quiver plot + #! TODO: Add sanity checks (x1, x2) = np.meshgrid( - frange(xlimv[0], xlimv[1], float(xlimv[1]-xlimv[0])/xlimv[2]), - frange(ylimv[0], ylimv[1], float(ylimv[1]-ylimv[0])/ylimv[2])); + frange(X[0], X[1], float(X[1]-X[0])/X[2]), + frange(Y[0], Y[1], float(Y[1]-Y[0])/Y[2])); - if ((not auto) and (not logtime) and (not timepts)): + if ((not autoFlag) and (not logtimeFlag) and (not timeptsFlag)): # Now calculate the vector field at those points (nr,nc) = x1.shape; dx = np.empty((nr, nc, 2)) @@ -151,16 +160,16 @@ #! TODO: Tweak the shape of the plot # a=gca; set(a,'DataAspectRatio',[1,1,1]); - # set(a,'XLim',xlimv(1:2)); set(a,'YLim',ylimv(1:2)); + # set(a,'XLim',X(1:2)); set(a,'YLim',Y(1:2)); mpl.xlabel('x1'); mpl.ylabel('x2'); # See if we should also generate the streamlines - if (xinit == None or len(xinit) == 0): + if (X0 == None or len(X0) == 0): return # Convert initial conditions to a numpy array - xinit = np.array(xinit); - (nr, nc) = np.shape(xinit); + X0 = np.array(X0); + (nr, nc) = np.shape(X0); # Generate some empty matrices to keep arrow information x1 = np.empty((nr, Narrows)); x2 = np.empty((nr, Narrows)); @@ -183,12 +192,12 @@ ymin = alim[2]; ymax = alim[3]; else: # Use the maximum extent of all trajectories - xmin = np.min(xinit[:,0]); xmax = np.max(xinit[:,0]); - ymin = np.min(xinit[:,1]); ymax = np.max(xinit[:,1]); + xmin = np.min(X0[:,0]); xmax = np.max(X0[:,0]); + ymin = np.min(X0[:,1]); ymax = np.max(X0[:,1]); # Generate the streamlines for each initial condition for i in range(nr): - state = odeint(odefun, xinit[i], TSPAN, args=parms); + state = odeint(odefun, X0[i], TSPAN, args=parms); time = TSPAN mpl.hold(True); mpl.plot(state[:,0], state[:,1]) @@ -197,7 +206,7 @@ # set(h[i], 'LineWidth', PP_stream_linewidth); # Plot arrows if quiver parameters were 'auto' - if (auto or logtime or timepts): + if (autoFlag or logtimeFlag or timeptsFlag): # Compute the locations of the arrows #! TODO: check this logic to make sure it works in python for j in range(Narrows): @@ -206,18 +215,18 @@ k = -1 if scale == None else 0; # Figure out what time index to use for the next point - if (auto): + if (autoFlag): # Use a linear scaling based on ODE time vector tind = np.floor((len(time)/Narrows) * (j-k)) + k; - elif (logtime): + elif (logtimeFlag): # Use an exponential time vector # MATLAB: tind = find(time < (j-k) / lambda, 1, 'last'); tarr = find(time < (j-k) / timefactor); tind = tarr[-1] if len(tarr) else 0; - elif (timepts): + elif (timeptsFlag): # Use specified time points - # MATLAB: tind = find(time < ylimv[j], 1, 'last'); - tarr = find(time < ylimv[j]); + # MATLAB: tind = find(time < Y[j], 1, 'last'); + tarr = find(time < timepts[j]); tind = tarr[-1] if len(tarr) else 0; # For tailless arrows, skip the first point @@ -263,3 +272,20 @@ # set(bp, 'MarkerSize', PP_arrow_markersize); return; + +# Utility function for generating initial conditions around a box +def boxgrid(xlimp, ylimp): + """BOXGRID generate list of points on edge of box + + list = BOXGRID([xmin xmax xnum], [ymin ymax ynum]) generates a + list of points that correspond to a uniform grid at the end of the + box defined by the corners [xmin ymin] and [xmax ymax]. + """ + + sx10 = frange(xlimp[0], xlimp[1], float(xlimp[1]-xlimp[0])/xlimp[2]) + sy10 = frange(ylimp[0], ylimp[1], float(ylimp[1]-ylimp[0])/ylimp[2]) + + sx1 = np.hstack((0, sx10, 0*sy10+sx10[0], sx10, 0*sy10+sx10[-1])) + sx2 = np.hstack((0, 0*sx10+sy10[0], sy10, 0*sx10+sy10[-1], sy10)) + + return np.transpose( np.vstack((sx1, sx2)) ) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |