Merge pull request #262 from roryyorke/rory/no-mplib-find

murrayrm · web-flow · commit bcf538287b37 · 2019-01-04T00:09:56.000-07:00
Fix: don't called deprecated Matplotlib functions find and frange
diff --git a/control/phaseplot.py b/control/phaseplot.py
@@ -39,12 +39,20 @@
 
 import numpy as np
 import matplotlib.pyplot as mpl
-from matplotlib.mlab import frange, find
+
 from scipy.integrate import odeint
 from .exception import ControlNotImplemented
 
 __all__ = ['phase_plot', 'box_grid']
 
+
+def _find(condition):
+    """Returns indices where ravel(a) is true.
+    Private implementation of deprecated matplotlib.mlab.find
+    """
+    return np.nonzero(np.ravel(condition))[0]
+
+
 def phase_plot(odefun, X=None, Y=None, scale=1, X0=None, T=None,
               lingrid=None, lintime=None, logtime=None, timepts=None,
               parms=(), verbose=True):
@@ -70,11 +78,11 @@ def phase_plot(odefun, X=None, Y=None, scale=1, X0=None, T=None,
         dxdt = F(x, t) that accepts a state x of dimension 2 and
         returns a derivative dx/dt of dimension 2.
 
-    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.
+    X, Y: 3-element sequences, optional, as [start, stop, npts]
+        Two 3-element sequences specifying x and y coordinates of a
+        grid.  These arguments are passed to linspace and meshgrid to
+        generate the points at which the vector field is plotted.  If
+        absent (or None), the vector field is not plotted.
 
     scale: float, optional
         Scale size of arrows; default = 1
@@ -145,8 +153,10 @@ def phase_plot(odefun, X=None, Y=None, scale=1, X0=None, T=None,
     #! TODO: Add sanity checks
     elif (X is not None and Y is not None):
         (x1, x2) = np.meshgrid(
-            frange(X[0], X[1], float(X[1]-X[0])/X[2]),
-            frange(Y[0], Y[1], float(Y[1]-Y[0])/Y[2]));
+            np.linspace(X[0], X[1], X[2]),
+            np.linspace(Y[0], Y[1], Y[2]))
+        Narrows = len(x1)
+
     else:
         # If we weren't given any grid points, don't plot arrows
         Narrows = 0;
@@ -234,12 +244,12 @@ def phase_plot(odefun, X=None, Y=None, scale=1, X0=None, T=None,
                 elif (logtimeFlag):
                     # Use an exponential time vector
                     # MATLAB: tind = find(time < (j-k) / lambda, 1, 'last');
-                    tarr = find(time < (j-k) / timefactor);
+                    tarr = _find(time < (j-k) / timefactor);
                     tind = tarr[-1] if len(tarr) else 0;
                 elif (timeptsFlag):
                     # Use specified time points
                     # MATLAB: tind = find(time < Y[j], 1, 'last');
-                    tarr = find(time < timepts[j]);
+                    tarr = _find(time < timepts[j]);
                     tind = tarr[-1] if len(tarr) else 0;
 
                 # For tailless arrows, skip the first point
@@ -295,8 +305,8 @@ def box_grid(xlimp, ylimp):
     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])
+    sx10 = np.linspace(xlimp[0], xlimp[1], xlimp[2])
+    sy10 = np.linspace(ylimp[0], ylimp[1], 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))
diff --git a/examples/genswitch.py b/examples/genswitch.py
@@ -8,7 +8,6 @@
 import numpy as np
 import matplotlib.pyplot as mpl
 from scipy.integrate import odeint
-from matplotlib.mlab import frange
 from control import phase_plot, box_grid
 
 # Simple model of a genetic switch
@@ -34,7 +33,7 @@ def genswitch(y, t, mu=4, n=2):
 sol2 = odeint(genswitch, sol1[-1,:] + [2, -2], tim2)
 
 # First plot out the curves that define the equilibria
-u = frange(0, 4.5, 0.1)
+u = np.linspace(0, 4.5, 46)
 f = np.divide(mu, (1 + u**n))   # mu / (1 + u^n), elementwise
 
 mpl.figure(1); mpl.clf();
diff --git a/examples/phaseplots.py b/examples/phaseplots.py
@@ -23,7 +23,7 @@ def invpend_ode(x, t, m=1., l=1., b=0.2, g=1):
 # Set up the figure the way we want it to look
 mpl.figure(); mpl.clf(); 
 mpl.axis([-2*pi, 2*pi, -2.1, 2.1]);
-mpl.title('Inverted pendlum')
+mpl.title('Inverted pendulum')
 
 # Outer trajectories
 phase_plot(invpend_ode,
@@ -35,9 +35,7 @@ def invpend_ode(x, t, m=1., l=1., b=0.2, g=1):
     logtime = (3, 0.7) )
 
 # Separatrices
-mpl.hold(True);
 phase_plot(invpend_ode, X0 = [[-2.3056, 2.1], [2.3056, -2.1]], T=6, lingrid=0)
-mpl.show();
 
 #
 # Systems of ODEs: damped oscillator example (simulation + phase portrait)
@@ -49,9 +47,10 @@ def oscillator_ode(x, t, m=1., b=1, k=1):
 # Generate a vector plot for the damped oscillator
 mpl.figure(); mpl.clf();
 phase_plot(oscillator_ode, [-1, 1, 10], [-1, 1, 10], 0.15);
-mpl.hold(True); mpl.plot([0], [0], '.');
+#mpl.plot([0], [0], '.');
 # a=gca; set(a,'FontSize',20); set(a,'DataAspectRatio',[1,1,1]);
-mpl.xlabel('x1'); mpl.ylabel('x2');
+mpl.xlabel('$x_1$'); mpl.ylabel('$x_2$');
+mpl.title('Damped oscillator, vector field')
 
 # Generate a phase plot for the damped oscillator
 mpl.figure(); mpl.clf(); 
@@ -61,11 +60,10 @@ def oscillator_ode(x, t, m=1., b=1, k=1):
     [-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);
+mpl.plot([0], [0], 'k.'); # 'MarkerSize', AM_data_markersize*3);
 # set(gca,'DataAspectRatio',[1,1,1]);
-mpl.xlabel('x1'); mpl.ylabel('x2');
-
-mpl.show()
+mpl.xlabel('$x_1$'); mpl.ylabel('$x_2$');
+mpl.title('Damped oscillator, vector field and stream lines')
 
 #
 # Stability definitions
@@ -88,9 +86,10 @@ def saddle_ode(x, t):
     [-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);
+mpl.plot([0], [0], 'k.'); # 'MarkerSize', AM_data_markersize*3);
 # set(gca,'FontSize', 16); 
-mpl.xlabel('{\itx}_1'); mpl.ylabel('{\itx}_2');
+mpl.xlabel('$x_1$'); mpl.ylabel('$x_2$');
+mpl.title('Asymptotically stable point')
 
 # Saddle
 mpl.figure(); mpl.clf();
@@ -103,9 +102,10 @@ def saddle_ode(x, t):
     [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] ], T = np.linspace(0, 2, 20));
-mpl.hold(True); mpl.plot([0], [0], 'k.'); # 'MarkerSize', AM_data_markersize*3);
+mpl.plot([0], [0], 'k.'); # 'MarkerSize', AM_data_markersize*3);
 # set(gca,'FontSize', 16); 
-mpl.xlabel('{\itx}_1'); mpl.ylabel('{\itx}_2');
+mpl.xlabel('$x_1$'); mpl.ylabel('$x_2$');
+mpl.title('Saddle point')
 
 # Stable isL
 m = 1; b = 0; k = 1;			# zero damping
@@ -115,8 +115,9 @@ def saddle_ode(x, t):
   [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);
+mpl.plot([0], [0], 'k.') # 'MarkerSize', AM_data_markersize*3);
 # set(gca,'FontSize', 16); 
-mpl.xlabel('{\itx}_1'); mpl.ylabel('{\itx}_2');
+mpl.xlabel('$x_1$'); mpl.ylabel('$x_2$');
+mpl.title('Undamped system\nLyapunov stable, not asympt. stable')
 
 mpl.show()
