Skip to content

Commit 53dfeb0

Browse files
murrayrmmurrayrm
committed
* 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
1 parent 19edd7b commit 53dfeb0

File tree

5 files changed

+236
-112
lines changed

5 files changed

+236
-112
lines changed

‎ChangeLog

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,17 @@
1+
2011-07-25 Richard Murray <murray@malabar.local>
2+
3+
* examples/phaseplots.py: updated calls to PhasePlot to use new
4+
argument structure
5+
6+
* src/phaseplot.py (PhasePlot): Updated call signature to be
7+
more pythonic and fixed up documentation.
8+
9+
* examples/genswitch.py (genswitch): added new example showing
10+
PhasePlot functionality
11+
12+
* src/phaseplot.py (boxgrid): added function to compute initial
13+
conditions around the edges of a box
14+
115
2011-07-24 Richard Murray <murray@malabar.local>
216

317
* tests/margin_test.py: added simple unit tests for margin functions

‎examples/genswitch.py

Lines changed: 81 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,81 @@
1+
# genswitch_plot.m - run the Collins genetic switch model
2+
# RMM, 24 Jan 07
3+
#
4+
# This file contains an example from FBS of a simple dynamical model
5+
# of a genetic switch. Plots time traces and a phase portrait using
6+
# the python-control library.
7+
8+
import numpy as np
9+
import matplotlib.pyplot as mpl
10+
from scipy.integrate import odeint
11+
from matplotlib.mlab import frange
12+
from control import PhasePlot, boxgrid
13+
14+
# Simple model of a genetic switch
15+
#
16+
# This function implements the basic model of the genetic switch
17+
# Parameters taken from Gardner, Cantor and Collins, Nature, 2000
18+
def genswitch(y, t, mu=4, n=2):
19+
return (mu / (1 + y[1]**n) - y[0], mu / (1 + y[0]**n) - y[1])
20+
21+
# Run a simulation from an initial condition
22+
tim1 = np.linspace(0, 10, 100)
23+
sol1 = odeint(genswitch, [1, 5], tim1)
24+
25+
# Extract the equlibirum points
26+
mu = 4; n = 2; # switch parameters
27+
eqpt = np.empty(3);
28+
eqpt[0] = sol1[0,-1]
29+
eqpt[1] = sol1[1,-1]
30+
eqpt[2] = 0; # fzero(@(x) mu/(1+x^2) - x, 2);
31+
32+
# Run another simulation showing switching behavior
33+
tim2 = np.linspace(11, 25, 100);
34+
sol2 = odeint(genswitch, sol1[-1,:] + [2, -2], tim2)
35+
36+
# First plot out the curves that define the equilibria
37+
u = frange(0, 4.5, 0.1)
38+
f = np.divide(mu, (1 + u**n)) # mu / (1 + u^n), elementwise
39+
40+
mpl.figure(1); mpl.clf();
41+
mpl.axis([0, 5, 0, 5]); # box on;
42+
mpl.plot(u, f, '-', f, u, '--') # 'LineWidth', AM_data_linewidth);
43+
mpl.legend('z1, f(z1)', 'z2, f(z2)') # legend(lgh, 'boxoff');
44+
mpl.plot([0, 3], [0, 3], 'k-') # 'LineWidth', AM_ref_linewidth);
45+
mpl.plot(eqpt[0], eqpt[1], 'k.', eqpt[1], eqpt[0], 'k.',
46+
eqpt[2], eqpt[2], 'k.') # 'MarkerSize', AM_data_markersize*3);
47+
mpl.xlabel('z1, f(z2)');
48+
mpl.ylabel('z2, f(z1)');
49+
50+
# Time traces
51+
mpl.figure(3); mpl.clf(); # subplot(221);
52+
mpl.plot(tim1, sol1[:,0], 'b-', tim1, sol1[:,1], 'g--');
53+
# set(pl, 'LineWidth', AM_data_linewidth);
54+
mpl.hold(True);
55+
mpl.plot([tim1[-1], tim1[-1]+1],
56+
[sol1[-1,0], sol2[0,1]], 'ko:',
57+
[tim1[-1], tim1[-1]+1], [sol1[-1,1], sol2[0,0]], 'ko:');
58+
# set(pl, 'LineWidth', AM_data_linewidth, 'MarkerSize', AM_data_markersize);
59+
mpl.plot(tim2, sol2[:,0], 'b-', tim2, sol2[:,1], 'g--');
60+
# set(pl, 'LineWidth', AM_data_linewidth);
61+
mpl.axis([0, 25, 0, 5]);
62+
63+
mpl.xlabel('Time {\itt} [scaled]');
64+
mpl.ylabel('Protein concentrations [scaled]');
65+
mpl.legend('z1 (A)', 'z2 (B)') # 'Orientation', 'horizontal');
66+
# legend(legh, 'boxoff');
67+
68+
# Phase portrait
69+
mpl.figure(2); mpl.clf(); # subplot(221);
70+
mpl.axis([0, 5, 0, 5]); # set(gca, 'DataAspectRatio', [1, 1, 1]);
71+
PhasePlot(genswitch, X0 = boxgrid([0, 5, 6], [0, 5, 6]), T = 10,
72+
timepts = [0.2, 0.6, 1.2])
73+
74+
# Add the stable equilibrium points
75+
mpl.hold(True);
76+
mpl.plot(eqpt[0], eqpt[1], 'k.', eqpt[1], eqpt[0], 'k.',
77+
eqpt[2], eqpt[2], 'k.') # 'MarkerSize', AM_data_markersize*3);
78+
79+
mpl.xlabel('Protein A [scaled]');
80+
mpl.ylabel('Protein B [scaled]'); # 'Rotation', 90);
81+

‎examples/phaseplots.py

Lines changed: 24 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -27,16 +27,16 @@ def invpend_ode(x, t, m=1., l=1., b=0.2, g=1):
2727

2828
# Outer trajectories
2929
PhasePlot(invpend_ode,
30-
'logtime', (3, 0.7), None,
31-
[ [-2*pi, 1.6], [-2*pi, 0.5], [-1.8, 2.1],
32-
[-1, 2.1], [4.2, 2.1], [5, 2.1],
33-
[2*pi, -1.6], [2*pi, -0.5], [1.8, -2.1],
34-
[1, -2.1], [-4.2, -2.1], [-5, -2.1] ],
35-
np.linspace(0, 40, 200))
30+
X0 = [ [-2*pi, 1.6], [-2*pi, 0.5], [-1.8, 2.1],
31+
[-1, 2.1], [4.2, 2.1], [5, 2.1],
32+
[2*pi, -1.6], [2*pi, -0.5], [1.8, -2.1],
33+
[1, -2.1], [-4.2, -2.1], [-5, -2.1] ],
34+
T = np.linspace(0, 40, 200),
35+
logtime = (3, 0.7) )
3636

3737
# Separatrices
3838
mpl.hold(True);
39-
PhasePlot(invpend_ode, 'auto', 0, None, [[-2.3056, 2.1], [2.3056, -2.1]], 6)
39+
PhasePlot(invpend_ode, X0 = [[-2.3056, 2.1], [2.3056, -2.1]], T=6, lingrid=0)
4040
mpl.show();
4141

4242
#
@@ -57,10 +57,10 @@ def oscillator_ode(x, t, m=1., b=1, k=1):
5757
mpl.figure(); mpl.clf();
5858
mpl.axis([-1, 1, -1, 1]); # set(gca, 'DataAspectRatio', [1, 1, 1]);
5959
PhasePlot(oscillator_ode,
60-
'timepts', [0.25, 0.8, 2, 3], None, [
61-
[-1, 1], [-0.3, 1], [0, 1], [0.25, 1], [0.5, 1], [0.75, 1], [1, 1],
62-
[1, -1], [0.3, -1], [0, -1], [-0.25, -1], [-0.5, -1], [-0.75, -1], [-1, -1]
63-
], np.linspace(0, 8, 80))
60+
X0 = [
61+
[-1, 1], [-0.3, 1], [0, 1], [0.25, 1], [0.5, 1], [0.75, 1], [1, 1],
62+
[1, -1], [0.3, -1], [0, -1], [-0.25, -1], [-0.5, -1], [-0.75, -1], [-1, -1]
63+
], T = np.linspace(0, 8, 80), timepts = [0.25, 0.8, 2, 3])
6464
mpl.hold(True); mpl.plot([0], [0], 'k.'); # 'MarkerSize', AM_data_markersize*3);
6565
# set(gca,'DataAspectRatio',[1,1,1]);
6666
mpl.xlabel('x1'); mpl.ylabel('x2');
@@ -81,26 +81,28 @@ def saddle_ode(x, t):
8181
m = 1; b = 1; k = 1; # default values
8282
mpl.figure(); mpl.clf();
8383
mpl.axis([-1, 1, -1, 1]); # set(gca, 'DataAspectRatio', [1 1 1]);
84-
PhasePlot(oscillator_ode, 'timepts', [0.3, 1, 2, 3], None,
85-
[[-1,1], [-0.3,1], [0,1], [0.25,1], [0.5,1], [0.7,1], [1,1], [1.3,1],
86-
[1,-1], [0.3,-1], [0,-1], [-0.25,-1], [-0.5,-1], [-0.7,-1], [-1,-1],
87-
[-1.3,-1]],
88-
np.linspace(0, 10, 100), parms = (m, b, k));
84+
PhasePlot(oscillator_ode,
85+
X0 = [
86+
[-1,1], [-0.3,1], [0,1], [0.25,1], [0.5,1], [0.7,1], [1,1], [1.3,1],
87+
[1,-1], [0.3,-1], [0,-1], [-0.25,-1], [-0.5,-1], [-0.7,-1], [-1,-1],
88+
[-1.3,-1]
89+
], T = np.linspace(0, 10, 100),
90+
timepts = [0.3, 1, 2, 3], parms = (m, b, k));
8991
mpl.hold(True); mpl.plot([0], [0], 'k.'); # 'MarkerSize', AM_data_markersize*3);
9092
# set(gca,'FontSize', 16);
9193
mpl.xlabel('{\itx}_1'); mpl.ylabel('{\itx}_2');
9294

9395
# Saddle
9496
mpl.figure(); mpl.clf();
9597
mpl.axis([-1, 1, -1, 1]); # set(gca, 'DataAspectRatio', [1 1 1]);
96-
PhasePlot(saddle_ode, 'timepts', [0.2, 0.5, 0.8], None,
98+
PhasePlot(saddle_ode, scale = 2, timepts = [0.2, 0.5, 0.8], X0 =
9799
[ [-1, -1], [1, 1],
98100
[-1, -0.95], [-1, -0.9], [-1, -0.8], [-1, -0.6], [-1, -0.4], [-1, -0.2],
99101
[-0.95, -1], [-0.9, -1], [-0.8, -1], [-0.6, -1], [-0.4, -1], [-0.2, -1],
100102
[1, 0.95], [1, 0.9], [1, 0.8], [1, 0.6], [1, 0.4], [1, 0.2],
101103
[0.95, 1], [0.9, 1], [0.8, 1], [0.6, 1], [0.4, 1], [0.2, 1],
102104
[-0.5, -0.45], [-0.45, -0.5], [0.5, 0.45], [0.45, 0.5],
103-
[-0.04, 0.04], [0.04, -0.04] ], np.linspace(0, 2, 20));
105+
[-0.04, 0.04], [0.04, -0.04] ], T = np.linspace(0, 2, 20));
104106
mpl.hold(True); mpl.plot([0], [0], 'k.'); # 'MarkerSize', AM_data_markersize*3);
105107
# set(gca,'FontSize', 16);
106108
mpl.xlabel('{\itx}_1'); mpl.ylabel('{\itx}_2');
@@ -109,10 +111,10 @@ def saddle_ode(x, t):
109111
m = 1; b = 0; k = 1; # zero damping
110112
mpl.figure(); mpl.clf();
111113
mpl.axis([-1, 1, -1, 1]); # set(gca, 'DataAspectRatio', [1 1 1]);
112-
PhasePlot(oscillator_ode, 'timepts',
113-
[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,
114-
[ [0.2,0], [0.4,0], [0.6,0], [0.8,0], [1,0], [1.2,0], [1.4,0] ],
115-
np.linspace(0, 20, 200), parms = (m, b, k));
114+
PhasePlot(oscillator_ode, timepts =
115+
[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],
116+
X0 = [ [0.2,0], [0.4,0], [0.6,0], [0.8,0], [1,0], [1.2,0], [1.4,0] ],
117+
T = np.linspace(0, 20, 200), parms = (m, b, k));
116118
mpl.hold(True); mpl.plot([0], [0], 'k.') # 'MarkerSize', AM_data_markersize*3);
117119
# set(gca,'FontSize', 16);
118120
mpl.xlabel('{\itx}_1'); mpl.ylabel('{\itx}_2');

‎src/__init__.py

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -64,6 +64,7 @@
6464
from mateqn import *
6565
from modelsimp import *
6666
from nichols import *
67+
from phaseplot import PhasePlot, boxgrid
6768
from rlocus import *
6869
from statefbk import *
6970
from statesp import *

0 commit comments

Comments
 (0)