Made package installable, added README, docs

.gitignore

+*.ini
+pytorchbridge.egg-info
+__pycache__/
\ No newline at end of file

README.md

+# `pyStateSpace`
+
+`pyStateSpace` is a package for simulating  a model using state space variables.
+
+## System modeling
+
+A system is described by its state-space variables `x`. The dynamics are described by the rate of change of those variables `dx/dt`.
+
+A system can respond to inputs `u`. The observable measurements of the system are its outputs `y`.
+
+In the general case, with no time invariance and linearity constraints:
+
+```
+dx/dt   = function(time, x, u)  --  (1)
+
+y       = function(time, x, u)  --  (2)
+```
+
+The value of `x` at each time is found by integrating `(1)` over the time step.
+
+## API
+
+All classes adhere to the `scikit-learn Estimator` API with a `predict(t: np.ndarray, x0:np.ndarray, u: np.ndarray) -> Tuple[np.ndarray, np.ndarray]` method. 
+
+All classes have these methods:
+
+* `dxdt(time: float, x: np.ndarray, u: np.ndarray) -> np.ndarray`
+* `y(time: float, x: np.ndarray, u: np.ndarray) -> np.ndarray`
+* `predict(t: np.ndarray, x0:np.ndarray, u: np.ndarray) -> Tuple[np.ndarray, np.ndarray]` The return value is a tuple of an array of state variables and an array of output variables `(x, y)`.
+
+Where `x`, `x0`, and `u` are vectors. `x0` is the initial state.
+
+Three classes are defined:
+
+### `Trapezoid`
+
+Uses trapezoid integration rule to integrate `(1)`. The `dxdt` and `y` methods need to be overridden.
+
+### `LTISystem`
+
+Same as `Trapezoid` but explicitly accepts linear time-invarnant system matrices `A, B, C, D` such that:
+
+```
+dx/dt   =   Ax + Bu
+y       =   Cx + Du
+```
+
+No methods need to be overridden.
+
+### `SS_ODE`
+
+An alternative to `Trapezoid` that uses `scipy`'s `odeint` function to integrate `dx/dt`. `dxdt` and `y` methods can be overridden or provided as arguments.
\ No newline at end of file

pystatespace/state_space.py

+"""
+Defines state space class for rough linear system modelling.
+"""
+from numbers import Number
+from typing import Callable, Tuple, Any
+
+import numpy as np
+from scipy.integrate import odeint
+
+
+
+class Trapezoid:
+    """
+    A state-space model that uses trapezoid integration to predict state variables.
+    `dxdt` and `y` methods need to be overridden.
+    """
+
+
+    def __init__(self, dims: int, outputs: int, tstep: float=1e-4):
+        """ 
+        Arguments:
+            dims {int} -- Number of state-space variables.
+            outputs {int} -- Number of output variables.
+        
+        Keyword Arguments:
+            tstep {float} -- Minimum size of simulation step (default: {1e-4})
+        """
+        self.dims = dims
+        self.outputs = outputs
+        self.tstep = tstep
+
+
+    def dxdt(self, t: float, x: np.ndarray, u: np.ndarray) -> np.ndarray:
+        raise NotImplementedError
+
+
+    def y(self, t: float, x: np.ndarray, u: np.ndarray) -> np.ndarray:
+        raise NotImplementedError
+
+
+    def predict(self, t: np.ndarray, x0: np.ndarray=None, u: np.ndarray=None) \
+        -> Tuple[np.ndarray, np.ndarray]:
+        """
+        Generate system response to inputs. The system response is calculated through
+        trapezoid rule.
+
+        Arguments:
+            t {np.ndarray} -- The end time (float) or an array of time stamps for
+                the corresponding element in `u`.
+            x0 {np.ndarray} -- The initial state of the system.
+            u {np.ndarray} -- An `N x Inputs` array of inputs,
+
+        Returns:
+            np.ndarray -- array of shape `N x Outputs`,
+        """
+        xdim = self.dims
+        if isinstance(t, Number):
+            t = np.linspace(0, t, len(u), endpoint=False)
+
+        xt = np.zeros(xdim) if x0 is None else x0
+        X = np.zeros((len(t), self.dims))
+        X[0] = xt
+        Y = np.zeros((len(t), self.outputs))
+        Y[0] = self.y(t[0], X[0], u[0])
+
+        # for t_, ut in zip(t[1:], u[1:]):
+        for i in range(1, len(t)):
+            Dt = t[i] - t[i - 1]        # Time step in supplied data
+            dt = min(Dt, self.tstep)    # Time step to take per calculation
+            while True:
+                dx = self.dxdt(t[i] - Dt, xt, u[i])
+                # xt += 0.5 * dt * (dx + dxprev)
+                xt += dt * dx
+
+                Dt -= dt
+                dt = min(Dt, self.tstep)
+                if dt == 0: break
+            X[i] = xt
+            Y[i] = self.y(t[i], xt, u[i])
+        return X, Y
+
+
+
+class LTISystem(Trapezoid):
+    """
+    A simple state-space simulator. Takes state-space matrices A, B, C, D and
+    simulates a linear time-invariant system response to provided inputs.
+    """
+
+    def __init__(self, A: np.ndarray, B: np.ndarray, C: np.ndarray, D: np.ndarray,
+                 tstep: float=1e-4):
+        """
+        Arguments:
+            A {np.ndarray} -- State variables x State variables,
+            B {np.ndarray} -- State variables x Inputs,
+            C {np.ndarray} -- Output Variables x State Variables,
+            D {np.ndarray} -- Output Variables x Inputs
+
+        Keyword Arguments:
+            tstep {float} -- Resolution of simulation (default: {1e-4})
+        """
+        self.A = A
+        self.B = B
+        self.C = C
+        self.D = D
+        self.tstep = tstep
+        self.dims = self.A.shape[0]
+        self.outputs = self.C.shape[0]
+
+
+    def dxdt(self, t: float, x: np.ndarray, u: np.ndarray) -> np.ndarray:
+        return self.A @ x + self.B @ u
+
+
+    def y(self, t: float, x: np.ndarray, u: np.ndarray) -> np.ndarray:
+        return self.C @ x + self.D @ u
+
+
+
+class SS_ODE:
+
+
+    def __init__(self, dims: int, outputs: int, dx: Callable, out: Callable,
+        tstep: float=1e-4):
+        """
+        Arguments:
+            dims {int} -- Number of state-space variables.
+            outputs {int} -- Number of output variables.
+            dx {Callable} -- A function with a signature (time, x, u) -> np.ndarray
+                that returns dx/dt. `x, u` are vectors. `time` is a float.
+            out {Callable} -- A function with a signature (time, x, u) -> np.ndarray
+                that returns y. `x, u` are vectors. `time` is a float.
+        
+        Keyword Arguments:
+            tstep {float} -- Maximum size of integration step (default: {1e-4})
+        """
+        self.dx = dx
+        self.out = out
+        self.tstep = tstep
+        self.dims = dims
+        self.outputs = outputs
+
+
+    def dxdt(self, t: float, x: np.ndarray, u: np.ndarray) -> np.ndarray:
+        return self.dx(t, x, u)
+
+
+    def y(self, t: float, x: np.ndarray, u: np.ndarray) -> np.ndarray:
+        return self.out(t, x, u)
+
+
+    def predict(self, t: Any[float, np.ndarray], x0: np.ndarray, u: np.ndarray) \
+        -> Tuple[np.ndarray, np.ndarray]:
+        if isinstance(t, Number):
+            t = np.linspace(0, t, len(u), endpoint=False)
+        t = np.concatenate(([t[0]], t), axis=0)
+        u = np.concatenate(([u[0]-u[0]], u), axis=0)
+        X = np.zeros((len(t), self.dims))
+        Y = np.zeros((len(t), self.outputs))
+        X[0] = np.zeros(self.dims) if x0 is None else x0
+        Y[0] = self.y(t[0], X[0], u[0])
+        for idx in range(1, len(t)):
+            X[idx] = odeint(func=self.dxdt,
+                            y0=np.squeeze(X[idx-1]),
+                            t=t[idx-1:idx+1],
+                            args=(u[idx],),
+                            tfirst=True,
+                            hmax=self.tstep)[-1]
+            Y[idx] = self.y(t[idx], X[idx], u[idx])
+        return np.asarray(X[1:]), np.asarray(Y[1:])

pystatespace/test.py

+"""
+Unit tests for code.
+"""
+
+
+
+from unittest import TestCase, main
+
+import numpy as np
+
+from .state_space import Trapezoid, LTISystem, SS_ODE
+
+
+
+class TestTrapezoid(TestCase):
+
+
+    def setUp(self):
+        class _Trapezoid(Trapezoid):
+            def dxdt(self, t, x, u):
+                return np.asarray([1., 1.])
+            def y(self, t, x, u):
+                return 1*x
+        self.solver = _Trapezoid(dims=2, outputs=2, tstep=0.5)
+
+
+    def test_integration(self):
+        u = np.zeros(10)
+        t = np.arange(10) * 2
+        x, y = self.solver.predict(t, np.asarray([0., 0.]), u)
+        self.assertEqual(len(y), len(t), 'Input/output length mismatch.')
+
+
+
+class TestSS_ODE(TestTrapezoid):
+
+
+    def setUp(self):
+        self.solver = SS_ODE(dims=2, outputs=2,
+                             dx=lambda t, x, u: np.asarray([1., 1.]),
+                             out=lambda t, x, u: 1*x,
+                             tstep=1.0)
+
+
+
+if __name__ == '__main__':
+    main()
\ No newline at end of file

setup.py

+from distutils.core import setup
+
+
+
+setup(name='pystatespace',
+      version='0.1.0',
+      packages=['pystatespace'],
+      install_requires=['numpy', 'scipy'],
+      author='Ibrahim Ahmed',
+      author_email='ibrahim.ahmed@vanderbilt.edu',
+      description='Scikit-learn Estimator API for state-space models',
+      url='https://git.isis.vanderbilt.edu/ahmedi/pyStateSpace',
+      classifiers=[
+          'Development Status :: 3 - Alpha',
+          'Programming Language :: Python :: 3.6',
+          'Programming Language :: Python :: 3.7'
+      ]
+)
\ No newline at end of file