Skip to content

Commit 2a919cc

Browse files
murrayrmmurrayrm
authored
Merge pull request #1112 from lkies/streamplot
Use matplotlibs streamplot function for phase_plane_plot
2 parents a042895 + fc1854e commit 2a919cc

File tree

11 files changed

+402
-85
lines changed

11 files changed

+402
-85
lines changed

‎control/phaseplot.py‎

Lines changed: 235 additions & 32 deletions
Large diffs are not rendered by default.

‎control/tests/ctrlplot_test.py‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -116,6 +116,11 @@ def setup_plot_arguments(resp_fcn, plot_fcn, compute_time_response=True):
116116
args2 = (sys2, )
117117
argsc = ([sys1, sys2], )
118118

119+
case (None, ct.phase_plane_plot):
120+
args1 = (sys1, )
121+
args2 = (sys2, )
122+
plot_kwargs = {'plot_streamlines': True}
123+
119124
case _, _:
120125
args1 = (sys1, )
121126
args2 = (sys2, )

‎control/tests/kwargs_test.py‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -173,6 +173,7 @@ def test_unrecognized_kwargs(function, nsssys, ntfsys, moreargs, kwargs,
173173
(control.phase_plane_plot, 1, ([-1, 1, -1, 1], 1), {}),
174174
(control.phaseplot.streamlines, 1, ([-1, 1, -1, 1], 1), {}),
175175
(control.phaseplot.vectorfield, 1, ([-1, 1, -1, 1], ), {}),
176+
(control.phaseplot.streamplot, 1, ([-1, 1, -1, 1], ), {}),
176177
(control.phaseplot.equilpoints, 1, ([-1, 1, -1, 1], ), {}),
177178
(control.phaseplot.separatrices, 1, ([-1, 1, -1, 1], ), {}),
178179
(control.singular_values_plot, 1, (), {})]
@@ -369,6 +370,7 @@ def test_response_plot_kwargs(data_fcn, plot_fcn, mimo):
369370
optimal_test.test_oep_argument_errors,
370371
'phaseplot.streamlines': test_matplotlib_kwargs,
371372
'phaseplot.vectorfield': test_matplotlib_kwargs,
373+
'phaseplot.streamplot': test_matplotlib_kwargs,
372374
'phaseplot.equilpoints': test_matplotlib_kwargs,
373375
'phaseplot.separatrices': test_matplotlib_kwargs,
374376
}

‎control/tests/phaseplot_test.py‎

Lines changed: 108 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@
1212
import warnings
1313
from math import pi
1414

15+
import matplotlib as mpl
1516
import matplotlib.pyplot as plt
1617
import numpy as np
1718
import pytest
@@ -138,29 +139,39 @@ def invpend_ode(t, x, m=0, l=0, b=0, g=0):
138139

139140
# Use callable form, with parameters (if not correct, will get /0 error)
140141
ct.phase_plane_plot(
141-
invpend_ode, [-5, 5, -2, 2], params={'args': (1, 1, 0.2, 1)})
142+
invpend_ode, [-5, 5, -2, 2], params={'args': (1, 1, 0.2, 1)},
143+
plot_streamlines=True)
142144

143145
# Linear I/O system
144146
ct.phase_plane_plot(
145-
ct.ss([[0, 1], [-1, -1]], [[0], [1]], [[1, 0]], 0))
147+
ct.ss([[0, 1], [-1, -1]], [[0], [1]], [[1, 0]], 0),
148+
plot_streamlines=True)
146149

147150

148151
@pytest.mark.usefixtures('mplcleanup')
149152
def test_phaseplane_errors():
150153
with pytest.raises(ValueError, match="invalid grid specification"):
151-
ct.phase_plane_plot(ct.rss(2, 1, 1), gridspec='bad')
154+
ct.phase_plane_plot(ct.rss(2, 1, 1), gridspec='bad',
155+
plot_streamlines=True)
152156

153157
with pytest.raises(ValueError, match="unknown grid type"):
154-
ct.phase_plane_plot(ct.rss(2, 1, 1), gridtype='bad')
158+
ct.phase_plane_plot(ct.rss(2, 1, 1), gridtype='bad',
159+
plot_streamlines=True)
155160

156161
with pytest.raises(ValueError, match="system must be planar"):
157-
ct.phase_plane_plot(ct.rss(3, 1, 1))
162+
ct.phase_plane_plot(ct.rss(3, 1, 1),
163+
plot_streamlines=True)
158164

159165
with pytest.raises(ValueError, match="params must be dict with key"):
160166
def invpend_ode(t, x, m=0, l=0, b=0, g=0):
161167
return (x[1], -b/m*x[1] + (g*l/m) * np.sin(x[0]))
162168
ct.phase_plane_plot(
163-
invpend_ode, [-5, 5, 2, 2], params={'stuff': (1, 1, 0.2, 1)})
169+
invpend_ode, [-5, 5, 2, 2], params={'stuff': (1, 1, 0.2, 1)},
170+
plot_streamlines=True)
171+
172+
with pytest.raises(ValueError, match="gridtype must be 'meshgrid' when using streamplot"):
173+
ct.phase_plane_plot(ct.rss(2, 1, 1), plot_streamlines=False,
174+
plot_streamplot=True, gridtype='boxgrid')
164175

165176
# Warning messages for invalid solutions: nonlinear spring mass system
166177
sys = ct.nlsys(
@@ -171,14 +182,87 @@ def invpend_ode(t, x, m=0, l=0, b=0, g=0):
171182
UserWarning, match=r"initial_state=\[.*\], solve_ivp failed"):
172183
ct.phase_plane_plot(
173184
sys, [-12, 12, -10, 10], 15, gridspec=[2, 9],
174-
plot_separatrices=False)
185+
plot_separatrices=False, plot_streamlines=True)
175186

176187
# Turn warnings off
177188
with warnings.catch_warnings():
178189
warnings.simplefilter("error")
179190
ct.phase_plane_plot(
180191
sys, [-12, 12, -10, 10], 15, gridspec=[2, 9],
181-
plot_separatrices=False, suppress_warnings=True)
192+
plot_streamlines=True, plot_separatrices=False,
193+
suppress_warnings=True)
194+
195+
@pytest.mark.usefixtures('mplcleanup')
196+
def test_phase_plot_zorder():
197+
# some of these tests are a bit akward since the streamlines and separatrices
198+
# are stored in the same list, so we separate them by color
199+
key_color = "tab:blue" # must not be 'k', 'r', 'b' since they are used by separatrices
200+
201+
def get_zorders(cplt):
202+
max_zorder = lambda items: max([line.get_zorder() for line in items])
203+
assert isinstance(cplt.lines[0], list)
204+
streamline_lines = [line for line in cplt.lines[0] if line.get_color() == key_color]
205+
separatrice_lines = [line for line in cplt.lines[0] if line.get_color() != key_color]
206+
streamlines = max_zorder(streamline_lines) if streamline_lines else None
207+
separatrices = max_zorder(separatrice_lines) if separatrice_lines else None
208+
assert cplt.lines[1] == None or isinstance(cplt.lines[1], mpl.quiver.Quiver)
209+
quiver = cplt.lines[1].get_zorder() if cplt.lines[1] else None
210+
assert cplt.lines[2] == None or isinstance(cplt.lines[2], list)
211+
equilpoints = max_zorder(cplt.lines[2]) if cplt.lines[2] else None
212+
assert cplt.lines[3] == None or isinstance(cplt.lines[3], mpl.streamplot.StreamplotSet)
213+
streamplot = max(cplt.lines[3].lines.get_zorder(), cplt.lines[3].arrows.get_zorder()) if cplt.lines[3] else None
214+
return streamlines, quiver, streamplot, separatrices, equilpoints
215+
216+
def assert_orders(streamlines, quiver, streamplot, separatrices, equilpoints):
217+
print(streamlines, quiver, streamplot, separatrices, equilpoints)
218+
if streamlines is not None:
219+
assert streamlines < separatrices < equilpoints
220+
if quiver is not None:
221+
assert quiver < separatrices < equilpoints
222+
if streamplot is not None:
223+
assert streamplot < separatrices < equilpoints
224+
225+
def sys(t, x):
226+
return np.array([4*x[1], -np.sin(4*x[0])])
227+
228+
# ensure correct zordering for all three flow types
229+
res_streamlines = ct.phase_plane_plot(sys, plot_streamlines=dict(color=key_color))
230+
assert_orders(*get_zorders(res_streamlines))
231+
res_vectorfield = ct.phase_plane_plot(sys, plot_vectorfield=True)
232+
assert_orders(*get_zorders(res_vectorfield))
233+
res_streamplot = ct.phase_plane_plot(sys, plot_streamplot=True)
234+
assert_orders(*get_zorders(res_streamplot))
235+
236+
# ensure that zorder can still be overwritten
237+
res_reversed = ct.phase_plane_plot(sys, plot_streamlines=dict(color=key_color, zorder=50), plot_vectorfield=dict(zorder=40),
238+
plot_streamplot=dict(zorder=30), plot_separatrices=dict(zorder=20), plot_equilpoints=dict(zorder=10))
239+
streamlines, quiver, streamplot, separatrices, equilpoints = get_zorders(res_reversed)
240+
assert streamlines > quiver > streamplot > separatrices > equilpoints
241+
242+
243+
@pytest.mark.usefixtures('mplcleanup')
244+
def test_stream_plot_magnitude():
245+
def sys(t, x):
246+
return np.array([4*x[1], -np.sin(4*x[0])])
247+
248+
# plt context with linewidth
249+
with plt.rc_context({'lines.linewidth': 4}):
250+
res = ct.phase_plane_plot(sys, plot_streamplot=dict(vary_linewidth=True))
251+
linewidths = res.lines[3].lines.get_linewidths()
252+
# linewidths are scaled to be between 0.25 and 2 times default linewidth
253+
# but the extremes may not exist if there is no line at that point
254+
assert min(linewidths) < 2 and max(linewidths) > 7
255+
256+
# make sure changing the colormap works
257+
res = ct.phase_plane_plot(sys, plot_streamplot=dict(vary_color=True, cmap='viridis'))
258+
assert res.lines[3].lines.get_cmap().name == 'viridis'
259+
res = ct.phase_plane_plot(sys, plot_streamplot=dict(vary_color=True, cmap='turbo'))
260+
assert res.lines[3].lines.get_cmap().name == 'turbo'
261+
262+
# make sure changing the norm at least doesn't throw an error
263+
ct.phase_plane_plot(sys, plot_streamplot=dict(vary_color=True, norm=mpl.colors.LogNorm()))
264+
265+
182266

183267

184268
@pytest.mark.usefixtures('mplcleanup')
@@ -190,7 +274,7 @@ def test_basic_phase_plots(savefigs=False):
190274
plt.figure()
191275
axis_limits = [-1, 1, -1, 1]
192276
T = 8
193-
ct.phase_plane_plot(sys, axis_limits, T)
277+
ct.phase_plane_plot(sys, axis_limits, T, plot_streamlines=True)
194278
if savefigs:
195279
plt.savefig('phaseplot-dampedosc-default.png')
196280

@@ -203,7 +287,7 @@ def invpend_update(t, x, u, params):
203287
ct.phase_plane_plot(
204288
invpend, [-2*pi, 2*pi, -2, 2], 5,
205289
gridtype='meshgrid', gridspec=[5, 8], arrows=3,
206-
plot_separatrices={'gridspec': [12, 9]},
290+
plot_separatrices={'gridspec': [12, 9]}, plot_streamlines=True,
207291
params={'m': 1, 'l': 1, 'b': 0.2, 'g': 1})
208292
plt.xlabel(r"$\theta$ [rad]")
209293
plt.ylabel(r"$\dot\theta$ [rad/sec]")
@@ -218,7 +302,8 @@ def oscillator_update(t, x, u, params):
218302
oscillator_update, states=2, inputs=0, name='nonlinear oscillator')
219303

220304
plt.figure()
221-
ct.phase_plane_plot(oscillator, [-1.5, 1.5, -1.5, 1.5], 0.9)
305+
ct.phase_plane_plot(oscillator, [-1.5, 1.5, -1.5, 1.5], 0.9,
306+
plot_streamlines=True)
222307
pp.streamlines(
223308
oscillator, np.array([[0, 0]]), 1.5,
224309
gridtype='circlegrid', gridspec=[0.5, 6], dir='both')
@@ -228,6 +313,18 @@ def oscillator_update(t, x, u, params):
228313
if savefigs:
229314
plt.savefig('phaseplot-oscillator-helpers.png')
230315

316+
plt.figure()
317+
ct.phase_plane_plot(
318+
invpend, [-2*pi, 2*pi, -2, 2],
319+
plot_streamplot=dict(vary_color=True, vary_density=True),
320+
gridspec=[60, 20], params={'m': 1, 'l': 1, 'b': 0.2, 'g': 1}
321+
)
322+
plt.xlabel(r"$\theta$ [rad]")
323+
plt.ylabel(r"$\dot\theta$ [rad/sec]")
324+
325+
if savefigs:
326+
plt.savefig('phaseplot-invpend-streamplot.png')
327+
231328

232329
if __name__ == "__main__":
233330
#

‎doc/figures/phaseplot-dampedosc-default.png‎

40.7 KB
Loading

‎doc/figures/phaseplot-invpend-meshgrid.png‎

-33.6 KB
Loading

‎doc/figures/phaseplot-oscillator-helpers.png‎

2 Bytes
Loading

‎doc/functions.rst‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -103,6 +103,7 @@ Phase plane plots
103103
phaseplot.separatrices
104104
phaseplot.streamlines
105105
phaseplot.vectorfield
106+
phaseplot.streamplot
106107

107108

108109
Frequency Response

‎doc/phaseplot.rst‎

Lines changed: 20 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@ functionality is supported by a set of mapping functions that are part of
1212
the `phaseplot` module.
1313

1414
The default method for generating a phase plane plot is to provide a
15-
2D dynamical system along with a range of coordinates and time limit:
15+
2D dynamical system along with a range of coordinates in phase space:
1616

1717
.. testsetup:: phaseplot
1818

@@ -27,8 +27,7 @@ The default method for generating a phase plane plot is to provide a
2727
sys_update, states=['position', 'velocity'],
2828
inputs=0, name='damped oscillator')
2929
axis_limits = [-1, 1, -1, 1]
30-
T = 8
31-
ct.phase_plane_plot(sys, axis_limits, T)
30+
ct.phase_plane_plot(sys, axis_limits)
3231

3332
.. testcode:: phaseplot
3433
:hide:
@@ -39,12 +38,12 @@ The default method for generating a phase plane plot is to provide a
3938
.. image:: figures/phaseplot-dampedosc-default.png
4039
:align: center
4140

42-
By default, the plot includes streamlines generated from starting
43-
points on limits of the plot, with arrows showing the flow of the
44-
system, as well as any equilibrium points for the system. A variety
41+
By default the plot includes streamlines infered from function values
42+
on a grid, equilibrium points and separatrices if they exist. A variety
4543
of options are available to modify the information that is plotted,
46-
including plotting a grid of vectors instead of streamlines and
47-
turning on and off various features of the plot.
44+
including plotting a grid of vectors instead of streamlines, plotting
45+
streamlines from arbitrary starting points and turning on and off
46+
various features of the plot.
4847

4948
To illustrate some of these possibilities, consider a phase plane plot for
5049
an inverted pendulum system, which is created using a mesh grid:
@@ -62,9 +61,7 @@ an inverted pendulum system, which is created using a mesh grid:
6261
invpend = ct.nlsys(invpend_update, states=2, inputs=1, name='invpend')
6362

6463
ct.phase_plane_plot(
65-
invpend, [-2 * np.pi, 2 * np.pi, -2, 2], 5,
66-
gridtype='meshgrid', gridspec=[5, 8], arrows=3,
67-
plot_equilpoints={'gridspec': [12, 9]},
64+
invpend, [-2 * np.pi, 2 * np.pi, -2, 2],
6865
params={'m': 1, 'l': 1, 'b': 0.2, 'g': 1})
6966
plt.xlabel(r"$\theta$ [rad]")
7067
plt.ylabel(r"$\dot\theta$ [rad/sec]")
@@ -79,16 +76,17 @@ an inverted pendulum system, which is created using a mesh grid:
7976

8077
This figure shows several features of more complex phase plane plots:
8178
multiple equilibrium points are shown, with saddle points showing
82-
separatrices, and streamlines generated along a 5x8 mesh of initial
83-
conditions. At each mesh point, a streamline is created that goes 5 time
84-
units forward and backward in time. A separate grid specification is used
85-
to find equilibrium points and separatrices (since the course grid spacing
86-
of 5x8 does not find all possible equilibrium points). Together, the
87-
multiple features in the phase plane plot give a good global picture of the
79+
separatrices, and streamlines generated generated from a rectangular
80+
25x25 grid (default) of function evaluations. Together, the multiple
81+
features in the phase plane plot give a good global picture of the
8882
topological structure of solutions of the dynamical system.
8983

90-
Phase plots can be built up by hand using a variety of helper functions that
91-
are part of the :mod:`phaseplot` (pp) module:
84+
Phase plots can be built up by hand using a variety of helper
85+
functions that are part of the :mod:`phaseplot` (pp) module. For more
86+
precise control, the streamlines can also generated by integrating the
87+
system forwards or backwards in time from a set of initial
88+
conditions. The initial conditions can be chosen on a rectangular
89+
grid, rectangual boundary, circle or from an arbitrary set of points.
9290

9391
.. testcode:: phaseplot
9492
:hide:
@@ -105,7 +103,8 @@ are part of the :mod:`phaseplot` (pp) module:
105103
oscillator = ct.nlsys(
106104
oscillator_update, states=2, inputs=0, name='nonlinear oscillator')
107105

108-
ct.phase_plane_plot(oscillator, [-1.5, 1.5, -1.5, 1.5], 0.9)
106+
ct.phase_plane_plot(oscillator, [-1.5, 1.5, -1.5, 1.5], 0.9,
107+
plot_streamlines=True)
109108
pp.streamlines(
110109
oscillator, np.array([[0, 0]]), 1.5,
111110
gridtype='circlegrid', gridspec=[0.5, 6], dir='both')
@@ -128,6 +127,7 @@ The following helper functions are available:
128127
phaseplot.equilpoints
129128
phaseplot.separatrices
130129
phaseplot.streamlines
130+
phaseplot.streamplot
131131
phaseplot.vectorfield
132132

133133
The :func:`phase_plane_plot` function calls these helper functions

0 commit comments

Comments
 (0)