~kdsch/ode2dsp

ce5e508b6fb05322f476d3d2095e1226f9bbdd41 — Karl Schultheisz a month ago 121fcab master
Implement stability report; fix some bugs
4 files changed, 177 insertions(+), 39 deletions(-)

M ode2dsp/example.py
M ode2dsp/faust.py
M ode2dsp/math.py
M ode2dsp/test.py
M ode2dsp/example.py => ode2dsp/example.py +14 -3
@@ 1,18 1,18 @@
from os.path import join
from subprocess import run
from sympy import symbols, Function, Matrix, Tuple, sinh
from .math import symbol_matrix, function_matrix, ODE, backward_euler
from .math import symbol_matrix, function_matrix, ODE, backward_euler, Report
from .faust import File, Printer

x, s = symbols("x s", cls=Function, real=True)
t = symbols("t", real=True)
k = symbols("k", integer=True)
Omega, beta = symbols("Omega beta", positive=True)
Omega, beta, f_s = symbols("Omega beta f_s", positive=True)
y = function_matrix("y", 2, t)
n = symbol_matrix("n", 2)

nonlinear_svf = ODE(
    t,
    Tuple(x(t)),
    y,
    Matrix(
        [


@@ 22,6 22,17 @@ nonlinear_svf = ODE(
    ).subs(s, sinh),
)

r = Report.from_spec(
    nonlinear_svf,
    backward_euler,
    k,
    f_s,
    Matrix([1, 1]),
    symbol_matrix("p", 2),
)

print(r)

faust_file = File.with_function_from_ode(
    "nonlinear_svf",
    k,

M ode2dsp/faust.py => ode2dsp/faust.py +6 -1
@@ 99,6 99,7 @@ def symbolize(f: Expr):

    raise BaseException(f"can't symbolize {f}")


class Delay(Expr):
    """
    Represents a time-delay expression, where expr is the


@@ 185,7 186,11 @@ def _advance_var(func, state, time) -> Expr:


def _is_signal(expr, k):
    return isinstance(expr, AppliedUndef) and len(expr.args) == 1 and Tuple(expr.args).has(k)
    return (
        isinstance(expr, AppliedUndef)
        and len(expr.args) == 1
        and Tuple(expr.args).has(k)
    )


def _get_shift(signal, k):

M ode2dsp/math.py => ode2dsp/math.py +99 -23
@@ 6,12 6,15 @@ from sympy import (
    Eq,
    Function,
    Matrix,
    pprint,
    pretty,
    re,
    solve,
    Subs,
    Symbol,
    symbols,
    Tuple,
    ZeroMatrix,
)
from typing import Any
from sympy.core.function import AppliedUndef


@@ 23,21 26,49 @@ def symbol_matrix(name, length):
    return Matrix(symbols(f"{name}_:{length}", real=True))


def function_matrix(name, length, arg):
def function_matrix(name, length, *args):
    return Matrix(
        [func(arg) for func in symbols(f"{name}_:{length}", cls=Function, real=True)]
        [func(*args) for func in symbols(f"{name}_:{length}", cls=Function, real=True)]
    )


def _everything(ode, discretize, time, sample_freq, n, point):
    fde = FDE.from_ode(ode, discretize, time, sample_freq)
    solved_fde = SolvedFDE.from_fde(fde, n)
    return (
        ode.stable(point),
        fde.stable(point),
        solved_fde.stable(point),
        solved_fde,
    )
class Report:
    def __init__(self, ode, fde, sfde, ode_stable, fde_stable, sfde_stable):
        self.ode = ode
        self.fde = fde
        self.sfde = sfde
        self.ode_stable = ode_stable
        self.fde_stable = fde_stable
        self.sfde_stable = sfde_stable

    @classmethod
    def from_spec(cls, ode, discretize, time, sample_freq, n, point):
        fde = FDE.from_ode(ode, discretize, time, sample_freq)
        solved_fde = SolvedFDE.from_fde(fde, n)
        return Report(
            ode,
            fde,
            solved_fde,
            ode.stable(point),
            fde.stable(point),
            solved_fde.stable(point),
        )

    def __str__(self):
        items = (
            ("ODE", self.ode),
            ("ODE stability condition", self.ode_stable),
            ("FDE", self.fde),
            ("FDE stability condition", self.fde_stable),
            ("Solved FDE", self.sfde),
            ("Solved FDE stability condition", self.sfde_stable),
        )

        def format(arg):
            label, value = arg
            return f"{label}:\n\n{pretty(value)}"

        return "\n\n\n".join(map(format, items))


class ODE(Basic):


@@ 49,24 80,37 @@ class ODE(Basic):
    #
    # See https://ccrma.stanford.edu/~jos/pasp/General_Nonlinear_ODE.html
    #
    def __init__(self, time: Symbol, state: Matrix, expr: Matrix):
    def __init__(self, time: Symbol, inputs, state: Matrix, expr: Matrix):
        self.time = time
        self.inputs = inputs
        self.state = state
        self.expr = expr

    def stable(self, operating_point: Matrix):
        return _ode_stability_condition(
            self.expr, self.state, self.time, operating_point
    def _autonomize(self):
        return ODE(
            self.time,
            (),
            self.state,
            self.expr.subs([(input, 0) for input in self.inputs]),
        )

    def stable(self, operating_point: Matrix):
        ode = self._autonomize()
        return _ode_stability_condition(ode.expr, ode.state, ode.time, operating_point)

    def __str__(self):
        rep = Eq(self.state.diff(self.time), self.expr)
        return pretty(rep)


class FDE(Basic):
    # Assumes the form
    #
    #    0 = expr(state, time)
    #
    def __init__(self, time: Symbol, state: Matrix, expr: Matrix):
    def __init__(self, time: Symbol, inputs, state: Matrix, expr: Matrix):
        self.time = time
        self.inputs = inputs
        self.state = state
        self.expr = expr



@@ 74,15 118,27 @@ class FDE(Basic):
    def from_ode(cls, ode: ODE, discretize, time: Symbol, sample_freq: Symbol):
        return FDE(
            time,
            ode.inputs.subs(ode.time, time),
            ode.state.subs(ode.time, time),
            discretize(ode.state, ode.expr, ode.time, time, sample_freq),
        )

    def stable(self, operating_point: Matrix):
        return _fde_stability_condition(
            self.expr, self.state, self.time, operating_point
    def _autonomize(self):
        return FDE(
            self.time,
            (),
            self.state,
            self.expr.subs([(input, 0) for input in self.inputs]),
        )

    def stable(self, operating_point: Matrix):
        fde = self._autonomize()
        return _fde_stability_condition(fde.expr - fde.state, fde.state, fde.time, operating_point)

    def __str__(self):
        rep = Eq(ZeroMatrix(*self.state.shape), self.expr)
        return pretty(rep)


class SolvedFDE(Basic):
    # Assumes the form


@@ 91,8 147,9 @@ class SolvedFDE(Basic):
    #
    # which must be an explicit equation.
    #
    def __init__(self, time: Symbol, state: Matrix, expr: Matrix):
    def __init__(self, time: Symbol, inputs, state: Matrix, expr: Matrix):
        self.time = time
        self.inputs = inputs
        self.state = state
        self.expr = expr



@@ 100,15 157,32 @@ class SolvedFDE(Basic):
    def from_fde(cls, fde: FDE, n):
        return SolvedFDE(
            fde.time,
            fde.inputs,
            fde.state,
            _solve_fde(fde.state, fde.expr, fde.time, n),
        )

    def _autonomize(self):
        return SolvedFDE(
            self.time,
            (),
            self.state,
            self.expr.subs([(input, 0) for input in self.inputs]),
        )

    def stable(self, operating_point: Matrix):
        sfde = self._autonomize()
        return _fde_stability_condition(
            self.expr, self.state, self.time, operating_point
            sfde.expr,
            sfde.state,
            sfde.time,
            operating_point,
        )

    def __str__(self):
        rep = Eq(self.state, self.expr)
        return pretty(rep)


def backward_euler(y, f, t, k, f_s):
    return Matrix(


@@ 186,7 260,7 @@ def _linearize(func, state, point):
    f[m](x) ≅  Σ  f[m](p[n]) +  ----- f[m] |   (x[n] - p[n])
              n=0               ∂x[n]      |x=p
    """
    substs = [(state[n], point[n]) for n in range(func.shape[0])]
    substs = [(state[n], point[n]) for n in range(state.shape[0])]
    return Matrix(
        [
            Add(


@@ 305,4 379,6 @@ def _fde_stable(y_k, f_k_, k):


def _fde_stability_condition(f_k, y_k, k, point):
    return _fde_stable(y_k, _block_dc(_linearize(f_k, y_k, point) - y_k, k), k)
    state = Matrix(list(y_k) + list(y_k.subs(k, k - 1)))
    point_ = Matrix(list(point) + list(point))
    return _fde_stable(y_k, _block_dc(_linearize(f_k, state, point_), k), k)

M ode2dsp/test.py => ode2dsp/test.py +58 -12
@@ 57,9 57,13 @@ class Test:

def test_linearize(ts):
    t = symbols("t", real=True)
    k = symbols("k", integer=True)
    Omega, beta, f_s = symbols("Omega beta f_s", positive=True)
    s = symbols("s", cls=Function, real=True)
    y = function_matrix("y", 2, t)
    p = symbol_matrix("p", 2)
    p_0, p_1 = tuple(p)
    y_0, y_1 = map(lambda f: f.func, y)

    cases = (
        {


@@ 69,6 73,8 @@ def test_linearize(ts):
                    s(y[0]) + s(y[1]),
                ]
            ).subs(s, sinh),
            "state": y,
            "point": p,
            "want": Matrix(
                [
                    (y[0] - p[0]) * cosh(p[0] + p[1])


@@ 81,9 87,43 @@ def test_linearize(ts):
                ]
            ),
        },
        {
            "func": Matrix(
                [
                    y_0(k - 1)
                    - (Omega * y_1(k) - sinh(y_0(k - 1))) / (-f_s - cosh(y_0(k - 1))),
                    (-2 * Omega * sinh(y_0(k)) + f_s * y_1(k - 1)) / (beta + f_s),
                ]
            ),
            "state": Matrix(list(y.subs(t, k)) + list(y.subs(t, k - 1))),
            "point": Matrix(list(p) + list(p)),
            # This is wrong because it's nonlinear with respect to delayed states.
            "want": Matrix(
                [
                    [
                        -Omega * (-p_1 + y_1(k)) / (-f_s - cosh(p_0))
                        + p_0
                        + (-p_0 + y_0(k - 1))
                        * (
                            1
                            + cosh(p_0) / (-f_s - cosh(p_0))
                            - (Omega * p_1 - sinh(p_0))
                            * sinh(p_0)
                            / (-f_s - cosh(p_0)) ** 2
                        )
                        - (Omega * p_1 - sinh(p_0)) / (-f_s - cosh(p_0))
                    ],
                    [
                        -2 * Omega * (-p_0 + y_0(k)) * cosh(p_0) / (beta + f_s)
                        + f_s * (-p_1 + y_1(k - 1)) / (beta + f_s)
                        + (-2 * Omega * sinh(p_0) + f_s * p_1) / (beta + f_s)
                    ],
                ]
            ),
        },
    )
    for i, c in enumerate(cases):
        got = _linearize(c["func"], y, p)
        got = _linearize(c["func"], c["state"], c["point"])
        if got != c["want"]:
            ts.fail(i, c["want"], got)



@@ 99,6 139,7 @@ def test_ode_stable(ts):
        {
            "ode": ODE(
                t,
                (),
                y,
                Matrix(
                    [


@@ 112,6 153,7 @@ def test_ode_stable(ts):
        {
            "ode": ODE(
                t,
                (),
                y,
                Matrix(
                    [


@@ 128,6 170,7 @@ def test_ode_stable(ts):
        {
            "ode": ODE(
                t,
                (),
                y,
                Matrix(
                    [


@@ 141,6 184,7 @@ def test_ode_stable(ts):
        {
            "ode": ODE(
                t,
                (),
                y,
                Matrix(
                    [


@@ 213,14 257,14 @@ def test_fde_from_ode(ts):

    cases = (
        {
            "ode": ODE(t, y, y),
            "ode": ODE(t, Tuple(), y, y),
            "discretize": backward_euler,
            "want": FDE(k, yy, -f_s * (yy - yy_) + yy),
            "want": FDE(k, Tuple(), yy, -f_s * (yy - yy_) + yy),
        },
        {
            "ode": ODE(t, y, y),
            "ode": ODE(t, Tuple(), y, y),
            "discretize": trapezoidal,
            "want": FDE(k, yy, -2 * f_s * (yy - yy_) + yy + yy_),
            "want": FDE(k, Tuple(), yy, -2 * f_s * (yy - yy_) + yy + yy_),
        },
    )
    for i, c in enumerate(cases):


@@ 236,11 280,11 @@ def test_fde_stable(ts):

    cases = (
        {
            "fde": FDE(k, y, y.subs(k, k - 1)),
            "fde": FDE(k, (), y, y.subs(k, k - 1)),
            "want": False,
        },
        {
            "fde": FDE(k, y, Matrix([0, 0])),
            "fde": FDE(k, (), y, Matrix([0, 0])),
            "want": True,
        },
    )


@@ 258,12 302,13 @@ def test_solved_fde_from_fde(ts):
    y_ = y.subs(k, k - 1)
    cases = (
        {
            "fde": FDE(k, y, -y + y_),
            "want": SolvedFDE(k, y, y_),
            "fde": FDE(k, (), y, -y + y_),
            "want": SolvedFDE(k, Tuple(), y, y_),
        },
        {
            "fde": FDE(
                k,
                (),
                y,
                Matrix(
                    [


@@ 274,6 319,7 @@ def test_solved_fde_from_fde(ts):
            ),
            "want": SolvedFDE(
                k,
                Tuple(),
                y,
                Matrix(
                    [


@@ 302,11 348,11 @@ def test_solved_fde_stable(ts):
    y_ = y.subs(k, k - 1)
    cases = (
        {
            "fde": SolvedFDE(k, y, y_),
            "fde": SolvedFDE(k, Tuple(), y, y_),
            "want": False,
        },
        {
            "fde": SolvedFDE(k, y, Matrix([0, 0])),
            "fde": SolvedFDE(k, Tuple(), y, Matrix([0, 0])),
            "want": True,
        },
    )


@@ 401,7 447,7 @@ def test_function_from_spec(ts):
        {
            "params": Tuple(),
            "outputs": Tuple(*y),
            "fde": SolvedFDE(k, y, y.subs(k, k - 1)),
            "fde": SolvedFDE(k, Tuple(), y, y.subs(k, k - 1)),
            "want": FunctionDefn(
                "func",
                (),