ENH: Addition/multiplication operators for StateSpace systems #7483

befelix · 2017-06-11T18:40:31Z

This PR adds addition and multiplication operators for StateSpace systems, see also #2973. This allows for things like

s1 = StateSpace(1, 1, 1, 0)
s2 = SateSpace(2, 2, 2, 0)
new = s1 + s2
new = s1 - s2
new = s1 * s2
new = 2 * s2

These operations have a frequency-domain interpretation, where adding means the systems operate in parallel and their output is added, while multiplying s1 * s2 means that the output of s2 is applied as the input to system s1.

rgommers · 2017-06-12T10:48:01Z

scipy/signal/ltisys.py

+            return StateSpace(a, b, c, d)
+
+        # A scalar/matrix is really just a static system (A=0, B=0, C=0)
+        return StateSpace(self.A, self.B, self.C, self.D + other)


this may not always do the right thing or give a sensible exception. You probably have to check what other is and raise an exception if it's not of the right type.

Wouldn't numpy complain on its own already if other is something that's not compatible?

I'm thinking about things that work when added to an array, e.g. if other is a masked array, or a pandas or xarray class. Or even just a higher-dimensional ndarray with a shape that D broadcasts to.

Anyway, not 100% sure you can write a sensible check, but worth thinking about.

I see, yes broadcasting could cause some really strange bugs here - nice catch! I think

other = np.atleast_2d(other) if self.D.shape != other.shape: raise Error

should cover almost all use cases anyways

rgommers · 2017-06-12T10:51:03Z

Looks like a good idea. What did you decide on for the minimal realisations you suggested in #2973 (comment)?

ilayn · 2017-06-12T11:37:44Z

I wouldn't automatically minreal it. There are many cases where people won't minreal things until the very end such as model reduction and also FEM based modeling. It should at least have an option to not to minreal it .

befelix · 2017-06-13T00:32:09Z

So I'm not doing anything in terms of minimal realisations. This seemed like the best approach - especially for a function that will typically not be called directly.

befelix · 2017-06-18T19:54:09Z

I've removed some functionality, such as __radd__ because it could lead to extremely weird behavior when dealing with numpy arrays and using __radd__ and __add__ is equivalent.

I've kept __rmul__ around because it is useful, but it does allow for things like this

from scipy import signal
sys = signal.StateSpace(1, 1, 1, 1)
np.array([1, 2]) * sys
Out[6]: 
array([ StateSpaceContinuous(
array([[1]]),
array([[1]]),
array([[1]]),
array([[1]]),
dt: None
),
       StateSpaceContinuous(
array([[1]]),
array([[1]]),
array([[2]]),
array([[2]]),
dt: None
)], dtype=object)

So you get back an object array with a system for each entry in the array. Maybe the best would be to strictly only allow multiplication of two StateSpace systems for now.

Down the road one could introduce a convenient classmethod for static systems. Then one could achieve the desired behavior with something like StateSpace.static(2) * sys.

Thoughts?

ilayn · 2017-06-18T20:16:11Z

Oh you are opening a can of worms 😃 And even then, it is still not watertight thanks to NumPy's broadcasting behavior.

shoyer · 2017-06-24T00:37:07Z

@befelix you need to set __array_priority__ to take precedence over numpy arrays in binary operations.

pv · 2017-07-30T12:47:12Z

Some comments: the class should set __array_priority__ = some_large_value and __array_ufunc__ = None to avoid binary op issues caused by Numpy, as noted by shoyer above.

Also, the binary operations should be written in a way that they (i) check if they can deal with type(other) input type and return NotImplemented if not, and (ii) if the input is not handled for some other reason, raise an error.

Made those changes in 520b5cd (and pushed to contributor branch). Was there something else that's still needed for this PR before merging?

larsoner

@ilayn do you want to look again? LGTM aside from my comments about six.

@stefanv @e-q or @endolith you might also be interested in looking if you have time

larsoner · 2017-07-30T19:58:08Z

scipy/signal/ltisys.py

@@ -26,6 +26,7 @@
 # np.linalg.qr fails on some tests with LinAlgError: zgeqrf returns -7
 # use scipy's qr until this is solved

+import six


This should be from scipy._lib.six, no? (And same for the tests code presumably.)

yep, fixed.

ilayn · 2017-07-30T22:55:21Z

scipy/signal/ltisys.py

+        H(s) = E1(s) + E2(s), means that applying H(s) to U(s)
+        is equivalent to applying E1(s) and E2(s) separately and adding up the
+        result.
+        """


I think this explanation is a tautology which already follows from LTIness of the systems. What is more informative might be mentioning that this is parallel connection behavior.

ilayn · 2017-07-30T22:56:18Z

scipy/signal/ltisys.py

+            # A scalar/matrix is really just a static system (A=0, B=0, C=0)
+            return StateSpace(self.A, self.B, self.C, self.D + other)
+        else:
+            raise ValueError('Other needs to have compatible dimensions.')


I guess users won't know what "other" means if they are not into class definitions in Python.

ilayn · 2017-07-30T23:00:50Z

scipy/signal/ltisys.py

+            return NotImplemented
+
+        return (-self).__add__(other)
+


While we are at it and if you agree; Quickly scaling the model with G/6.5 is also possible with __truediv__() method where you can simply test for float or int as you do for mult and then return self*(1/other). You can also explicitly reject __rtruediv__() with a more meaningful error message.

Disable numpy binops and ufuncs via setting __array_priority__ and __array_ufunc__. Check type(other) in binops, and return NotImplemented if it's unsupported. For any other causes of operation not being supported, raise errors.

…omments

…after binops

pv · 2017-08-18T17:57:29Z

Updated with added __truediv__ (only for scalars -- let's not start guessing if the user meant matrix inversion, because / doesn't have that meaning elsewhere). There's no need for __rtruediv__, as the default error messages seem fine.

pv · 2017-10-15T13:17:55Z

other comments? merge?

larsoner · 2017-10-15T19:04:58Z

Nothing from me, let's merge Tuesday if nobody else wants to take a look.

rgommers added enhancement A new feature or improvement scipy.signal labels Jun 12, 2017

rgommers reviewed Jun 12, 2017

View reviewed changes

befelix force-pushed the ss_operators branch from c537dd1 to 81633e1 Compare June 18, 2017 20:02

larsoner reviewed Jul 30, 2017

View reviewed changes

ilayn reviewed Jul 30, 2017

View reviewed changes

befelix and others added 7 commits August 18, 2017 19:34

ENH: StateSpace: Implemented add, multiply and subtract operators

680e6be

TST: StateSpace: Add tests for +,-,* operators

3d0f618

ENH: Stricter type checking

ade1890

BUG: signal: fix up StateSpace binop behavior

cf8d904

Disable numpy binops and ufuncs via setting __array_priority__ and __array_ufunc__. Check type(other) in binops, and return NotImplemented if it's unsupported. For any other causes of operation not being supported, raise errors.

BUG: signal: import scipy._lib.six, not six

0cba9f7

ENH: signal: implement StateSpace.__truediv__ + adjust messages and c…

ffab369

…omments

BUG: signal/StateSpace: ensure coefficient arrays have the same type …

447099b

…after binops

pv force-pushed the ss_operators branch from ec59d97 to 447099b Compare August 18, 2017 17:56

pv merged commit e9b861b into scipy:master Apr 8, 2018

pv added this to the 1.1.0 milestone Apr 8, 2018

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

ENH: Addition/multiplication operators for StateSpace systems #7483

ENH: Addition/multiplication operators for StateSpace systems #7483

befelix commented Jun 11, 2017

rgommers Jun 12, 2017

befelix Jun 13, 2017

rgommers Jun 14, 2017

befelix Jun 18, 2017

rgommers commented Jun 12, 2017

ilayn commented Jun 12, 2017

befelix commented Jun 13, 2017

befelix commented Jun 18, 2017

ilayn commented Jun 18, 2017 •

edited

Loading

shoyer commented Jun 24, 2017

pv commented Jul 30, 2017 •

edited

Loading

larsoner left a comment

larsoner Jul 30, 2017

pv Jul 30, 2017

ilayn Jul 30, 2017

ilayn Jul 30, 2017

ilayn Jul 30, 2017

pv commented Aug 18, 2017 •

edited

Loading

pv commented Oct 15, 2017

larsoner commented Oct 15, 2017

ENH: Addition/multiplication operators for StateSpace systems #7483

ENH: Addition/multiplication operators for StateSpace systems #7483

Conversation

befelix commented Jun 11, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rgommers commented Jun 12, 2017

ilayn commented Jun 12, 2017

befelix commented Jun 13, 2017

befelix commented Jun 18, 2017

ilayn commented Jun 18, 2017 • edited Loading

shoyer commented Jun 24, 2017

pv commented Jul 30, 2017 • edited Loading

larsoner left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

pv commented Aug 18, 2017 • edited Loading

pv commented Oct 15, 2017

larsoner commented Oct 15, 2017

ilayn commented Jun 18, 2017 •

edited

Loading

pv commented Jul 30, 2017 •

edited

Loading

pv commented Aug 18, 2017 •

edited

Loading