proposal to standardize squeezing output from systems

[sawyerbfuller]

Proposal: When computing frequency and time response output of systems, standardize on indexing first element so output is a 1D array if system is SISO and squeeze=True. This is instead of applying numpy.squeeze to output.

Currently, there are two ways that python-control squeezes output:

1. apply 'squeeze' to the output at the end. This is what functions in timeresp.py do, e.g. forced_response (edit: and iosys time responses). An m x n MIMO system (m-output, n-input, m,n>1) produces a [m x n x k] array, where k is the number of data points in time. But a 1 x n system's output gets squeezed to [n x k] and a m x 1 is squeezed to [m x k].

2. use indexing to extract the first element if SISO, otherwise leave a full array. This is what the frequency response functions evalfr, freqresp do. If SISO, give [k] length array; MIMO systems always produce [m x n x k] arrays.

I propose to standardize on the latter: a MIMO system should always output [m x n x k] arrays. Rationale: this facilitates interconnection and keeps indexing the system and its outputs consistent. Along with the proposal: retain the squeeze=True keyword argument because it's already in use, but do the 'squeeze' by indexing the first element, and only do it if SISO.

Documentation would read something like:

"If squeeze is True (default) and sys is single input single output (SISO), returns a 1D array equal to the length k of the input, otherwise returns an (n_outputs, n_inputs, k) array for multiple input multiple output (MIMO) systems."

This is relevant for #449 and popped up in #442.

[sawyerbfuller]

By the way, as far as I can tell this last came up ... almost ten years ago.

https://sourceforge.net/p/python-control/mailman/message/27671663/

Looks like somehow the conclusion was to do something like the above options, but it was not explicitly spelled out. Somehow the frequency response functions ended up doing something different than the time response functions.

Just realized, one of the participants, who I didn’t know at the time, is now a colleague of mine. Hi @eigensteve!

[ilayn]

This has been always an annoying thing. I would actually expect the squeeze to be the default behavior so it has to be turned off with squeeze=False since True is the 99% of the use cases.

[sawyerbfuller]

Hi @ilayn, yes, the proposal includes that squeeze=True by default, so that SISO systems should always return 1d arrays unless you specifically want them to behave like MIMO systems -- because SISO systems are something like 99% of use cases.

Just to confirm, would that be your preference?

Can you see any reason to you would rather, when squeeze=True, to keep using numpy.squeeze on all systems rather than indexing the zeroth element only for SISO systems?

[ilayn]

As fas as I experienced, you only need keepdims or squeeze True if you are doing a batch loop and you don't want to care about if the system in that spin of the loop is SISO or MIMO. For the rest, I can't think of any other usage.

[sawyerbfuller]

ok cool thanks. I'll tally that as a vote for #2 (please correct me if I'm wrong)

[murrayrm]

I think this is something we should address in 0.9.0, since it is somewhat related to return of array vs matrix objects. Tagging it with that milestone.

[murrayrm]

I'm taking a cut at this. PR coming soon, but some comments on things I found:

• For the time response systems, the input is given and so even for a MIMO system with m inputs and p outputs, we get a 2D output of shape (p, k). This means that you might have a multi-input system and still get a system output. So what I did was to interpret squeeze=True as acting when there is a single system output (versus a SISO system).
• I also found an inconsistency in the way that return_x handled: in all time response functions except forced_response the default return values are T, Y and if you set return_x=True then you get T, Y, X. I fixed that so that everything is consistent.

[sawyerbfuller]

I see, i didnt realize that you have to select an input in the forced response functions, so you always get a (p, k) sized output for a MIMO system unless squeeze.

I think i was under the mistaken impression that the output would have m dimensions, each 'row' of which would correspond to the response of the system to e.g. a step or impulse at the m'th input.

Assuming we keep things as they are, that is, that an input must be selected, then I think this issue is actually a non-issue. edit: see below

One minor difference is that in freq response functions, squeezing is achieved by selectig the zeroth element iff SISO whereas in time response it is currently accomplished by squeeze and no SISO test.

I think best to use freq response behavior in time response functions, ie test for SISO and only squeeze if SISO and squeeze. reason: to get consistent behavior eg in batch jobs with multiple systems with varying numbers of outputs, that @ilayn refers to above.

[murrayrm]

If we require that the system be SISO, then you could have the following scenarios:

• Forced response for a 2-input, 1-output system: you provide the two inputs (explicitly) and call forced response to get the single output. Because the system is not SISO, you would not be able to get a squeezed (1D) output. So you would have to index the (single) output in order to get back a 1D vector.

• Forced response for a 1-input, 2-output system: you select the second output as the output that you care about for step_response to get the single output. Because the system is not SISO, you would not be able to get a squeezed (1D) output. So again you would have to index the (single) output in order to get back a 1D vector.

These are counter-intuitive to me: it seems like if you call the function with squeeze and there is only a single signal as output, you should get back a squeezed output (1D). Thoughts (from all)?

[bnavigator]

I agree with @murrayrm. If you set squeeze=True (the default) and expect only one series you should get it as (k,)

If you batch jobs with multiple systems, the most natural way is to explicitly set squeeze=False and expect (1, 1, k) or (1, k) even for SISO systems.

[sawyerbfuller]

Let's see ..

• dominant use case: code that is always going to use siso systems: using defalts works
• rarer use case: code that generally always works regardless of mimo or siso: can use squeeze=false and be assured output arrays are same num dims

So the corner case is code with interspersed SISO systems, and MIMO systems that sometimes have one output and using teh defaults. Output as currently coded will have different num dims depending on num outputs. Thsi is a pretty insignificant use case.

I believe that SISO and MIMO are sort of like two different classes of system and that mimo systems should always produce mimo-like outputs.

But if we are keeping the squeeze term then that suggests the numpy squeeze function and Im ok with functionality as it is currently in #507 in time resp functions. But then the freq response functions should have the same behavior and not test for SISO before squeezing along first two dimensions.

[murrayrm]

Good point on frequency response. To make sure I understand it: if you compute the frequency response for a 2 input, 1 output system, right now you get a 3D array (roughly Y[output][input][freq]), independent of whether squeeze is true or false. What you are proposing, I think, is that if you specify squeeze=True, this would become Y[input][freq].

Similarly, if you have a 1 input, 2 output system and you squeeze, you get Y[output][freq] as the response.

Right?

[bnavigator]

I see a difference in time series vs frequency responses: A time series should always be at least (k, ) even if k=1. For frequency responses, I can see value in the current behavior, that if I give one frequency as input, then I want a scalar value returned.

[sawyerbfuller]

@murrayrm Yes that sounds right.

Not sure exactly how the code would look to squeeze along the first two dimensions in the freqresp functions thoigh.

[sawyerbfuller]

@bnavigator right now the docs and behavior return an array even for a single freuqency.

[sawyerbfuller]

i am a little worried that the addition of squeeze as a defaukt may break code that assumes a SISO system returns a 1D output when the default value is false. (might be hard to catch all instances of sys.__call__ in code that makes this assumption)

I think the "right" way to do this might be to be object-oriebted: to eliminate the squeeze keyword altogether and add an attribute to systems that allow a siso system to behave like a MIMO system (that is, sys.issiso() returns false even if it actually siso). don't know how soon I could get to that though, but i could put it on my list if it siunds interesting

[murrayrm]

I think the following are all agreeable:

• For SISO systems, the default should be that you don't need to subscript inputs, outputs, etc. I think all variants have this property.

• It should be possible to treat a SISO system like a MIMO systems (the current squeeze=False behavior currently satisfies that).

• H(s) should return a scalar for SISO and 2D array for non-SISO (including SIMO and MISO).

• Time and frequency responses should be handled in a consistent fashion, so that that the behavior of squeeze in "understandable" across those classes of functions.

Assuming those are agreeable, a couple of possibilities:

1. Full numpy: when squeeze is true, all single dimensional entries in the shape are removed (with an axis keyword if you want to only squeeze certain axes?).

2. Modified numpy: when squeeze is true, all single dimensional input/output entries in the shape of the output variables are removed, but single dimensional entries in time/frequency are maintained.

3. Current behavior: when squeeze is true, if all input/output entries are single dimensional, then single dimensional input/output entries are removed in the shape. This one pushes a bit on "consistency" because for a MISO system, the frequency response would return a 3D system ([input][output][freq]) but a time response would return a 1D system (indexed by time alone).

4. SISO-only: if squeeze is true and the system is SISO, then single dimensional input/output entries are removed in the shape. This one pushes a bit on "understandable" because for a MISO system, the "squeezed" forced response would return a 2D system ([output][time]) even if though there is only one output.

5. More customizable: we can of course use the squeeze keyword more flexibly, so that we could have squeeze='all', squeeze='inputs', etc to pick up the different cases.

6. Additional variant: allow a SISO system to be "locked" as a MIMO system (this would probably require a lot of care since I suspect we have various places where we check the value of inputs and outputs instead of calling issiso.

Preferences? Other options?

[murrayrm]

FYI, I have moved all of the processing of time and frequency responses into separate functions (lti._process_frequency_response, timeresp._process_time_response), so it should be easy to make everything consistent once we decide what we want things to look like.

[sawyerbfuller]

Nice synopsis. I agree on the four main points. and in particular agree that a frequency response of a siso system to a scalar should probably be a scalar.

I think 5 is interesting but too complicated.

1 is probably good combo of simplicity & least change.

My favorite is still 6 because it feels most object oriented, is opinionated, and would require the shortest docstring.

But advocates of a feature have a responsibility to implement it, and my doc says it will be a few more weeks till I can type with both hands again (broken bone). And as @murrayrm suggested, the devil may be in the details. Changes I foresee include an update to the ‘conventions’ section to describe the ‘as_mimo’ feature and adding links to conventions on the time and freq response functions.

My suggestion: do 1 and if a few weeks have gone by and 0.9 still isn’t released then maybe a PR will appear & you can take a look. But this isn’t important enough to hold up a a release waiting on my PR.

[bnavigator]

1 sounds good to me.

Except that I would prefer to keep it 1-d for every function when there is a vector involved in the output. In time series or in frequency responses. So H(s) -> scalar iff s is scalar, but freqresp(sys, omega) -> Tuple[array, array, array]

[murrayrm]

I feel like #2 or #3 is most useful, since that lets you squeeze out dimensions that you don't need when you only have a "single" response. However, in all cases except when you happen to evaluate at a single time or frequency (using a list), then I think it is equivalent to #1.

A related question is what the default should be. If we go with #1 and we leave squeeze=True as the default, then asking for a frequency response on a singleton list (which could come up) will not generate a list of frequency responses, but just a scalar. So you would have to use squeeze=False and then np.squeeze with axis to get a 1D list.

A middle ground (but perhaps overkill) would be to let squeeze be a more general keyword:

• sqeeze=False: keep all dimensions
• squeeze=True (default): same as numpy => remove all singleton dimensions
• squeeze='inputs': squeeze input dimensions, if possible (only matters for frequency response)
• squeeze='outputs': squeeze output dimensions, if possible (time and frequency response)
• squeeze='signals': squeeze input and output dimensions
• squeeze='all': same as numpy

[sawyerbfuller]

I like 2 better than 3 in that it keeps the output looking like the input, wheter it is a unit length array vs scalar. (I dont think its necessary to return a list if the user suppled frequencies as a list. am ok with returning all non-scalar outputs as numpy arrays.)

i dont folow @bnavigator 's comment: should freqresp always return arays regardless of whether given omega is a scalar, or should it return a scalar for a scalar input? both freqresp and call are now vectorized.

ultimately, I think @murrayrm has most experience with the lobrary and as cheif architect he should get final say.

[murrayrm]

OK, here's the current plan:

• squeeze=True (default): same as numpy => remove all singleton dimensions
• squeeze=False: keep all dimensions
• squeeze='siso': squeeze input and output dimensions iff the system is SISO (leave time/freq alone)
• squeeze='io': squeeze input/output dimensions only (leave time/freq alone)

This has the advantage that it captures all current behaviors in both time and frequency response, so if someone was relying on one of those then they can recover it by specifying squeeze explicitly.

[sawyerbfuller]

Hi, this sounds good to me. S

[murrayrm]

I started to implement the changes required to implement the scheme above, but it was pretty messy and it lead to some counter-intuitive behavior (to me):

• If we set squeeze=True as the default and we decide that squeeze=True should do a numpy squeeze (all axes), then both sys(s) and sys([s]) will return the same thing (a 0D or 2D array). This is different than current behavior, where you get 0D/2D for sys(s) and 1D/3D for sys([s]).

Because of this, I decided to pull out the time response portion of PR #507 into a separate PR (#511), focused on tightening the code and fixing an inconsistency with return_x. I've converted #507 into a draft PR for now while I/we sort out how to best proceed. It should be possible to make things consistent (the original issue that @sawyerbfuller pointed out at the top of this) once we sort out what we really want.

A couple of paths forward:

• Make the proposed change that squeeze=True is the same as numpy, and so it squeezes everything it can. Set this to the default case, changing current behavior.
• Keep something closer to what we currently have, which squeezes inputs and outputs, but not frequencies (so H(s) returns a 0D or 2D array but H([s]) returns a 1D or 3D array).
• Something more complicated (which was what I was starting to do and then thought better of it).

One thought is that we might think about having the following three behaviors:

• squeeze=None: squeeze inputs and outputs, but leave frequencies alone. This would be the default and is consistent with current behavior. In particular
  • sys(s) returns a scalar or a 2D array (SISO vs MIMO)
  • sys([s]) returns a 1D/3D array (with the frequency axis having dimension 1)
• squeeze=True: numpy squeeze => all singleton axes are removed. In particular
• • sys(s) and sys([s]) will return the same thing
• • For a SIMO system or a MISO system sys([s1, s2, ..., sn]) returns a 2D array
• squeeze=False: full MIMO, even for SISO. You always get a m x p x k array.

The main change in the code is that instead of having squeeze=True as the default in the various functions, we set squeeze=None and then you have to set squeeze=True to get further squeezing.

For time responses, we could do the same thing. In particular, if you set squeeze=None then you would only remove output dimensions if the system is SISO. If you want to remove dimensions for a MISO system, you need to use squeeze=True.

[murrayrm]

PR's are now in place: #507 for frequency response, #511 for time response. The behavior described above is implemented, with one modification: the result from squeeze=False does not modify the default behavior for the frequency axis. In particular:

• siso_sys(s, squeeze=False) returns an m x p array if s is a scalar
• siso_sys([s], squeeze=False) returns an m x p x 1 array

I think this addresses the original problem that @sawyerbfuller raised. In particular, both time and frequency responses have the following behavior:

• Default case squeezes the inputs and outputs only for SISO systems
• Setting squeeze=True performs a numpy squeeze on the output
• Setting squeeze=False treats SISO like MIMO

I have set up separate issue (#512) to track the creation of arrays of step and impulse responses. When that is implemented, time and frequency responses should have completely consistent behavior.

[sawyerbfuller]

Richard’s proposal seems good to me. In #507, docstrings don’t yet indicate behavior of squeeze =none but I assume that is to come assuming we go forward with this.

[sawyerbfuller]

A further thought: am I right that the only difference now in Richard’s proposal is that squeeze=true squeezes all axes, changing current behavior? Just want to make sure that is an important enough use case to bother with changing anything. I think we are all in agreement that the output should match input (scalar-> scalar for siso, array-like -> array for siso), and squeeze=true as proposed would change that, which nobody is asking for. Feel free to let me know if I’m missing something. Is a bit complicated.

[murrayrm]

For frequency responses, you are basically correct. The "main" cases (SISO or full MIMO) are essentially the same. The main differences are what happens when you have a SIMO systems, a MISO system.

The current behavior for a frequency response is that squeeze=True is the default and it squeezes the input and output axes if the system is SISO, leaving the frequency axis alone. So if you have a 1-output, 2-input system then response = sys([s], squeeze=True) does nothing.

In the current PR #507 , you get the same default behavior (using squeeze=None or just leaving it off) or you can get a (2,)-shaped return value using squeeze=True.

For time responses, things are currently different, but PR #511 will make things consistent with whatever we decide (it is currently set up to match PR #507).

What I think I like about the proposed system is that you can basically get three basic behaviors that are (relatively) easy to articulate and understand:

• Avoid unneeded dimensions in SISO systems, but keep in MIMO systems (default case, using squeeze=None)
• Expand SISO systems to return the same as MIMO systems (using squeeze=False)
• Squeeze everything consistent with numpy.squeeze (using squeeze=True)

[sawyerbfuller]

Ok, sounds good, sorry for causing extra work.