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",
(),