update squeeze processing for freq response + unit tests, doc updates

murrayrm · murrayrm · commit 90da4fbdf102 · 2021-01-17T12:33:13.000-08:00
diff --git a/control/config.py b/control/config.py
@@ -16,7 +16,7 @@
 # Package level default values
 _control_defaults = {
     'control.default_dt': 0,
-    'control.squeeze_frequency_response': True
+    'control.squeeze_frequency_response': None
 }
 defaults = dict(_control_defaults)
 
diff --git a/control/frdata.py b/control/frdata.py
@@ -353,25 +353,33 @@ def eval(self, omega, squeeze=None):
 
         Parameters
         ----------
-        omega : float or array_like
+        omega : float or 1D array_like
             Frequencies in radians per second
-        squeeze : bool, optional (default=True)
-            If True and the system is single-input single-output (SISO),
-            return a 1D array rather than a 3D array.  Default value (True)
-            set by config.defaults['control.squeeze_frequency_response'].
+        squeeze : bool, optional
+            If squeeze=True, remove single-dimensional entries from the shape
+            of the output even if the system is not SISO. If squeeze=False,
+            keep all indices (output, input and, if omega is array_like,
+            frequency) even if the system is SISO. The default value can be
+            set using config.defaults['control.squeeze_frequency_response'].
 
         Returns
         -------
-        fresp : (self.outputs, self.inputs, len(x)) or (len(x), ) complex ndarray
-            The frequency response of the system. Array is ``len(x)``
-            if and only if system is SISO and ``squeeze=True``.
+        fresp : complex ndarray
+            The frequency response of the system.  If the system is SISO and
+            squeeze is not True, the shape of the array matches the shape of
+            omega.  If the system is not SISO or squeeze is False, the first
+            two dimensions of the array are indices for the output and input
+            and the remaining dimensions match omega.  If ``squeeze`` is True
+            then single-dimensional axes are removed.
 
         """
-        # Set value of squeeze argument if not set
-        if squeeze is None:
-            squeeze = config.defaults['control.squeeze_frequency_response']
-
         omega_array = np.array(omega, ndmin=1) # array-like version of omega
+
+        # Make sure that we are operating on a simple list
+        if len(omega_array.shape) > 1:
+            raise ValueError("input list must be 1D")
+
+        # Make sure that frequencies are all real-valued
         if any(omega_array.imag > 0):
             raise ValueError("FRD.eval can only accept real-valued omega")
 
@@ -396,7 +404,7 @@ def eval(self, omega, squeeze=None):
 
     def __call__(self, s, squeeze=None):
         """Evaluate system's transfer function at complex frequencies.
-        
+
         Returns the complex frequency response `sys(s)` of system `sys` with
         `m = sys.inputs` number of inputs and `p = sys.outputs` number of
         outputs.
@@ -406,19 +414,24 @@ def __call__(self, s, squeeze=None):
 
         Parameters
         ----------
-        s : complex scalar or array_like
+        s : complex scalar or 1D array_like
             Complex frequencies
         squeeze : bool, optional (default=True)
-            If True and the system is single-input single-output (SISO),
-            return a 1D array rather than a 3D array.  Default value (True)
-            set by config.defaults['control.squeeze_frequency_response'].
+            If squeeze=True, remove single-dimensional entries from the shape
+            of the output even if the system is not SISO. If squeeze=False,
+            keep all indices (output, input and, if omega is array_like,
+            frequency) even if the system is SISO. The default value can be
+            set using config.defaults['control.squeeze_frequency_response'].
 
         Returns
         -------
-        fresp : (p, m, len(s)) complex ndarray or (len(s),) complex ndarray 
-            The frequency response of the system. Array is ``(len(s), )`` if
-            and only if system is SISO and ``squeeze=True``.
-
+        fresp : complex ndarray
+            The frequency response of the system.  If the system is SISO and
+            squeeze is not True, the shape of the array matches the shape of
+            omega.  If the system is not SISO or squeeze is False, the first
+            two dimensions of the array are indices for the output and input
+            and the remaining dimensions match omega.  If ``squeeze`` is True
+            then single-dimensional axes are removed.
 
         Raises
         ------
@@ -427,13 +440,14 @@ def __call__(self, s, squeeze=None):
             :class:`FrequencyDomainData` systems are only defined at imaginary
             frequency values.
         """
-        # Set value of squeeze argument if not set
-        if squeeze is None:
-            squeeze = config.defaults['control.squeeze_frequency_response']
+        # Make sure that we are operating on a simple list
+        if len(np.array(s, ndmin=1).shape) > 1:
+            raise ValueError("input list must be 1D")
 
         if any(abs(np.array(s, ndmin=1).real) > 0):
             raise ValueError("__call__: FRD systems can only accept "
                             "purely imaginary frequencies")
+
         # need to preserve array or scalar status
         if hasattr(s, '__len__'):
             return self.eval(np.asarray(s).imag, squeeze=squeeze)
diff --git a/control/lti.py b/control/lti.py
@@ -131,21 +131,26 @@ def frequency_response(self, omega, squeeze=None):
 
         Parameters
         ----------
-        omega : float or array_like
+        omega : float or 1D array_like
             A list, tuple, array, or scalar value of frequencies in
             radians/sec at which the system will be evaluated.
         squeeze : bool, optional
-            If True and the system is single-input single-output (SISO),
-            return a 1D array rather than a 3D array.  Default value (True)
-            set by config.defaults['control.squeeze_frequency_response'].
+            If squeeze=True, remove single-dimensional entries from the shape
+            of the output even if the system is not SISO. If squeeze=False,
+            keep all indices (output, input and, if omega is array_like,
+            frequency) even if the system is SISO. The default value can be
+            set using config.defaults['control.squeeze_frequency_response'].
 
         Returns
         -------
-        mag : (p, m, len(omega)) ndarray or (len(omega),) ndarray 
+        mag : ndarray
             The magnitude (absolute value, not dB or log10) of the system
-            frequency response. Array is ``(len(omega), )`` if
-            and only if system is SISO and ``squeeze=True``.
-        phase : (p, m, len(omega)) ndarray or (len(omega),) ndarray
+            frequency response.  If the system is SISO and squeeze is not
+            True, the array is 1D, indexed by frequency.  If the system is not
+            SISO or squeeze is False, the array is 3D, indexed by the output,
+            input, and frequency.  If ``squeeze`` is True then
+            single-dimensional axes are removed.
+        phase : ndarray
             The wrapped phase in radians of the system frequency response.
         omega : ndarray
             The (sorted) frequencies at which the response was evaluated.
@@ -482,18 +487,24 @@ def evalfr(sys, x, squeeze=None):
     ----------
     sys: StateSpace or TransferFunction
         Linear system
-    x : complex scalar or array_like
+    x : complex scalar or 1D array_like
         Complex frequency(s)
     squeeze : bool, optional (default=True)
-        If True and the system is single-input single-output (SISO), return a
-        1D array rather than a 3D array.  Default value (True) set by
+        If squeeze=True, remove single-dimensional entries from the shape of
+        the output even if the system is not SISO. If squeeze=False, keep all
+        indices (output, input and, if omega is array_like, frequency) even if
+        the system is SISO. The default value can be set using
         config.defaults['control.squeeze_frequency_response'].
 
     Returns
     -------
-    fresp : (p, m, len(x)) complex ndarray or (len(x),) complex ndarray 
-        The frequency response of the system. Array is ``(len(x), )`` if
-        and only if system is SISO and ``squeeze=True``.
+    fresp : complex ndarray
+        The frequency response of the system.  If the system is SISO and
+        squeeze is not True, the shape of the array matches the shape of
+        omega.  If the system is not SISO or squeeze is False, the first two
+        dimensions of the array are indices for the output and input and the
+        remaining dimensions match omega.  If ``squeeze`` is True then
+        single-dimensional axes are removed.
 
     See Also
     --------
@@ -519,7 +530,7 @@ def evalfr(sys, x, squeeze=None):
 
 def freqresp(sys, omega, squeeze=None):
     """Frequency response of an LTI system at multiple angular frequencies.
-    
+
     In general the system may be multiple input, multiple output (MIMO), where
     `m = sys.inputs` number of inputs and `p = sys.outputs` number of
     outputs.
@@ -528,23 +539,27 @@ def freqresp(sys, omega, squeeze=None):
     ----------
     sys: StateSpace or TransferFunction
         Linear system
-    omega : float or array_like
+    omega : float or 1D array_like
         A list of frequencies in radians/sec at which the system should be
         evaluated. The list can be either a python list or a numpy array
         and will be sorted before evaluation.
-    squeeze : bool, optional (default=True)
-        If True and the system is single-input single-output (SISO), return a
-        1D array rather than a 3D array.  Default value (True) set by
+    squeeze : bool, optional
+        If squeeze=True, remove single-dimensional entries from the shape of
+        the output even if the system is not SISO. If squeeze=False, keep all
+        indices (output, input and, if omega is array_like, frequency) even if
+        the system is SISO. The default value can be set using
         config.defaults['control.squeeze_frequency_response'].
 
     Returns
     -------
-    mag : (p, m, len(omega)) ndarray or (len(omega),) ndarray 
+    mag : ndarray
         The magnitude (absolute value, not dB or log10) of the system
-        frequency response. Array is ``(len(omega), )`` if and only if system 
-        is SISO and ``squeeze=True``.
-
-    phase : (p, m, len(omega)) ndarray or (len(omega),) ndarray 
+        frequency response.  If the system is SISO and squeeze is not True,
+        the array is 1D, indexed by frequency.  If the system is not SISO or
+        squeeze is False, the array is 3D, indexed by the output, input, and
+        frequency.  If ``squeeze`` is True then single-dimensional axes are
+        removed.
+    phase : ndarray
         The wrapped phase in radians of the system frequency response.
     omega : ndarray
         The list of sorted frequencies at which the response was
@@ -601,14 +616,30 @@ def dcgain(sys):
 
 # Process frequency responses in a uniform way
 def _process_frequency_response(sys, omega, out, squeeze=None):
+    # Set value of squeeze argument if not set
+    if squeeze is None:
+        squeeze = config.defaults['control.squeeze_frequency_response']
+
     if not hasattr(omega, '__len__'):
         # received a scalar x, squeeze down the array along last dim
         out = np.squeeze(out, axis=2)
 
+    #
     # Get rid of unneeded dimensions
-    if squeeze is None:
-        squeeze = config.defaults['control.squeeze_frequency_response']
-    if squeeze and sys.issiso():
+    #
+    # There are three possible values for the squeeze keyword at this point:
+    #
+    #   squeeze=None: squeeze input/output axes iff SISO
+    #   squeeze=True: squeeze all single dimensional axes (ala numpy)
+    #   squeeze-False: don't squeeze any axes
+    #
+    if squeeze is True:
+        # Squeeze everything that we can if that's what the user wants
+        return np.squeeze(out)
+    elif squeeze is None and sys.issiso():
+        # SISO system output squeezed unless explicitly specified otherwise
         return out[0][0]
-    else:
+    elif squeeze is False or squeeze is None:
         return out
+    else:
+        raise ValueError("unknown squeeze value")
diff --git a/control/margins.py b/control/margins.py
@@ -294,25 +294,25 @@ def stability_margins(sysdata, returnall=False, epsw=0.0):
             # frequency for gain margin: phase crosses -180 degrees
             w_180 = _poly_iw_real_crossing(num_iw, den_iw, epsw)
             with np.errstate(all='ignore'):  # den=0 is okay
-                w180_resp = evalfr(sys, 1J * w_180, squeeze=True)
+                w180_resp = evalfr(sys, 1J * w_180)
 
             # frequency for phase margin : gain crosses magnitude 1
             wc = _poly_iw_mag1_crossing(num_iw, den_iw, epsw)
-            wc_resp = evalfr(sys, 1J * wc, squeeze=True)
+            wc_resp = evalfr(sys, 1J * wc)
 
             # stability margin
             wstab = _poly_iw_wstab(num_iw, den_iw, epsw)
-            ws_resp = evalfr(sys, 1J * wstab, squeeze=True)
+            ws_resp = evalfr(sys, 1J * wstab)
 
         else:  # Discrete Time
             zargs = _poly_z_invz(sys)
             # gain margin
             z, w_180 = _poly_z_real_crossing(*zargs, epsw=epsw)
-            w180_resp = evalfr(sys, z, squeeze=True)
+            w180_resp = evalfr(sys, z)
 
             # phase margin
             z, wc = _poly_z_mag1_crossing(*zargs, epsw=epsw)
-            wc_resp = evalfr(sys, z, squeeze=True)
+            wc_resp = evalfr(sys, z)
 
             # stability margin
             z, wstab = _poly_z_wstab(*zargs, epsw=epsw)
@@ -437,11 +437,11 @@ def phase_crossover_frequencies(sys):
         omega = _poly_iw_real_crossing(num_iw, den_iw, 0.)
 
         # using real() to avoid rounding errors and results like 1+0j
-        gain = np.real(evalfr(sys, 1J * omega, squeeze=True))
+        gain = np.real(evalfr(sys, 1J * omega))
     else:
         zargs = _poly_z_invz(sys)
         z, omega = _poly_z_real_crossing(*zargs, epsw=0.)
-        gain = np.real(evalfr(sys, z, squeeze=True))
+        gain = np.real(evalfr(sys, z))
 
     return omega, gain
 
diff --git a/control/statesp.py b/control/statesp.py
@@ -645,7 +645,7 @@ def __call__(self, x, squeeze=None):
 
         Returns the complex frequency response `sys(x)` where `x` is `s` for
         continuous-time systems and `z` for discrete-time systems.
-        
+
         In general the system may be multiple input, multiple output
         (MIMO), where `m = self.inputs` number of inputs and `p =
         self.outputs` number of outputs.
@@ -657,18 +657,24 @@ def __call__(self, x, squeeze=None):
 
         Parameters
         ----------
-        x : complex or complex array_like
+        x : complex or complex 1D array_like
             Complex frequencies
         squeeze : bool, optional
-            If True and the system is single-input single-output (SISO),
-            return a 1D array rather than a 3D array.  Default value (True)
-            set by config.defaults['control.squeeze_frequency_response'].
+            If squeeze=True, remove single-dimensional entries from the shape
+            of the output even if the system is not SISO. If squeeze=False,
+            keep all indices (output, input and, if omega is array_like,
+            frequency) even if the system is SISO. The default value can be
+            set using config.defaults['control.squeeze_frequency_response'].
 
         Returns
         -------
-        fresp : (p, m, len(x)) complex ndarray or (len(x),) complex ndarray
-            The frequency response of the system. Array is ``len(x)`` if and
-            only if system is SISO and ``squeeze=True``.
+        fresp : complex ndarray
+            The frequency response of the system.  If the system is SISO and
+            squeeze is not True, the shape of the array matches the shape of
+            omega.  If the system is not SISO or squeeze is False, the first
+            two dimensions of the array are indices for the output and input
+            and the remaining dimensions match omega.  If ``squeeze`` is True
+            then single-dimensional axes are removed.
 
         """
         # Use Slycot if available
@@ -693,9 +699,13 @@ def slycot_laub(self, x):
             Frequency response
         """
         from slycot import tb05ad
+        x_arr = np.atleast_1d(x) # array-like version of x
+
+        # Make sure that we are operating on a simple list
+        if len(x_arr.shape) > 1:
+            raise ValueError("input list must be 1D")
 
         # preallocate
-        x_arr = np.atleast_1d(x) # array-like version of x
         n = self.states
         m = self.inputs
         p = self.outputs
@@ -755,6 +765,11 @@ def horner(self, x):
             # Fall back because either Slycot unavailable or cannot handle
             # certain cases.
             x_arr = np.atleast_1d(x) # force to be an array
+
+            # Make sure that we are operating on a simple list
+            if len(x_arr.shape) > 1:
+                raise ValueError("input list must be 1D")
+
             # Preallocate
             out = empty((self.outputs, self.inputs, len(x_arr)), dtype=complex)
 
diff --git a/control/tests/lti_test.py b/control/tests/lti_test.py
diff --git a/control/xferfcn.py b/control/xferfcn.py

Original file line number	Diff line number	Diff line change
`@@ -16,7 +16,7 @@`
`16`	`16`	`# Package level default values`
`17`	`17`	`_control_defaults = {`
`18`	`18`	`'control.default_dt': 0,`
`19`		`- 'control.squeeze_frequency_response': True`
	`19`	`+ 'control.squeeze_frequency_response': None`
`20`	`20`	`}`
`21`	`21`	`defaults = dict(_control_defaults)`
`22`	`22`