Repository navigation
proposal to standardize squeezing output from systems #453
Description
Activity
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!
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=Falsesince True is the 99% of the use cases.Hi @ilayn, yes, the proposal includes that
squeeze=Trueby 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 usingnumpy.squeezeon all systems rather than indexing the zeroth element only for SISO systems?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.
ok cool thanks. I'll tally that as a vote for #2 (please correct me if I'm wrong)
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.
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=Trueas acting when there is a single system output (versus a SISO system). - I also found an inconsistency in the way that
return_xhandled: in all time response functions exceptforced_responsethe default return values areT, Yand if you setreturn_x=Truethen you getT, Y, X. I fixed that so that everything is consistent.
- 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
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
mdimensions, each 'row' of which would correspond to the response of the system to e.g. a step or impulse at them'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 belowOne 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.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 responseto 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_responseto 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
squeezeand there is only a single signal as output, you should get back a squeezed output (1D). Thoughts (from all)?-
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=Falseand expect (1, 1, k) or (1, k) even for SISO systems.- linked a pull request that will close this issueStandardize squeeze processing in frequency response functions #507
on Jan 14, 2021 6 remaining items
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=Falsebehavior 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
squeezein "understandable" across those classes of functions.
Assuming those are agreeable, a couple of possibilities:
-
Full numpy: when
squeezeis true, all single dimensional entries in the shape are removed (with anaxiskeyword if you want to only squeeze certain axes?). -
Modified numpy: when
squeezeis 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. -
Current behavior: when
squeezeis 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). -
SISO-only: if
squeezeis 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. -
More customizable: we can of course use the
squeezekeyword more flexibly, so that we could havesqueeze='all',squeeze='inputs', etc to pick up the different cases. -
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
inputsandoutputsinstead of callingissiso.
Preferences? Other options?
Reacted by Ben Greiner-
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.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.
Reacted by Ben Greiner1 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) -> scalariff s is scalar, butfreqresp(sys, omega) -> Tuple[array, array, array]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=Falseand thennp.squeezewith axis to get a 1D list.A middle ground (but perhaps overkill) would be to let
squeezebe a more general keyword:sqeeze=False: keep all dimensionssqueeze=True(default): same as numpy => remove all singleton dimensionssqueeze='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 dimensionssqueeze='all': same as numpy
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.
OK, here's the current plan:
squeeze=True(default): same as numpy => remove all singleton dimensionssqueeze=False: keep all dimensionssqueeze='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
squeezeexplicitly.- Hi, this sounds good to me. S
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=Trueas the default and we decide thatsqueeze=Trueshould do a numpy squeeze (all axes), then bothsys(s)andsys([s])will return the same thing (a 0D or 2D array). This is different than current behavior, where you get 0D/2D forsys(s)and 1D/3D forsys([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=Trueis 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 particularsys(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)andsys([s])will return the same thing
-
- For a SIMO system or a MISO system
sys([s1, s2, ..., sn])returns a 2D array
- For a SIMO system or a MISO system
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=Trueas the default in the various functions, we setsqueeze=Noneand then you have to setsqueeze=Trueto get further squeezing.For time responses, we could do the same thing. In particular, if you set
squeeze=Nonethen you would only remove output dimensions if the system is SISO. If you want to remove dimensions for a MISO system, you need to usesqueeze=True.- If we set
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=Falsedoes 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 scalarsiso_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=Trueperforms a numpy squeeze on the output - Setting
squeeze=Falsetreats 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.
- 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.
- 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.
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=Trueis 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 thenresponse = sys([s], squeeze=True)does nothing.In the current PR #507 , you get the same default behavior (using
squeeze=Noneor just leaving it off) or you can get a (2,)-shaped return value usingsqueeze=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)
- Avoid unneeded dimensions in SISO systems, but keep in MIMO systems (default case, using
- Ok, sounds good, sorry for causing extra work.
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 applyingnumpy.squeezeto output.Currently, there are two ways that python-control squeezes output:
apply 'squeeze' to the output at the end. This is what functions in
timeresp.pydo, e.g.forced_response(edit: andiosystime 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].use indexing to extract the first element if SISO, otherwise leave a full array. This is what the frequency response functions
evalfr,freqrespdo. 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=Truekeyword 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.