-
Notifications
You must be signed in to change notification settings - Fork 60
Symbolic algebra
Sam Yates edited this page Nov 24, 2016
·
6 revisions
Related to work on issue #64, how can we well represent symbolic descriptions of ODEs as described by NMODL DERIVATIVE or KINETIC blocks, or by our as yet hypothetical alternative mechanism description language?
Symbolic manipulation tasks:
- identification of a system of ODEs as being linear or diagonal;
- determination of null-space of linear ODE system (or more generally, null-space with respect to monomials in a polynomial ODE system; relevant to kinetic schemes);
- reduction of a consistent over-determined linear (or polynomial) ODE system, or at least perform symbolic substitution;
- polynomial expansion and simplification;
- constant folding.
We need to be able to translate our AST descriptions of an NMODL-specified state evolution into this symbolic representation, apply the above tests and manipulations, and then translate this into code for evaluation within the context of some choice of integrator.
These seem to constitute the majority: GiNaC, CoCoA, Giac, Maxima, SymbolicC++. Due to licensing, we have to exclude them.
YaCaS is licensed under LGPL 2.1. From 1.6.0, Yacas builds a C++ library libyacas.a that provides the back end for its GUI and console drivers, but API documentation is scarce.
-
Axiom has a modified BSD license (refer to the
readmeon GitHub). It's very large, and has a long history. Embedding/linking would have to done via Embedded Common Lisp (ECL), or via communication with an external process. - FriCAS is an actively maintained fork of Axiom, and also appears to be under a modified BSD license. Same caveats apply.
- Reduce also has a modified BSD license; it can run under Codemist Standard Lisp (included in the distribution) which apparently is suitable for linking with C.
- Symengine (note the use a reference-counted pointer implementation from Trilinos which is under the BSD license). It is native C++, but API documentation is very scant. Fairly large distribution at circa 120k lines, 75k in library proper.
This is definitely a possibility: our demands of a symbolic algebra system are fairly simple.
It looks like symengine will indeed do everything we need. Exploratory/demo code below.
#include <iostream>
#include <symengine/add.h>
#include <symengine/basic.h>
#include <symengine/derivative.h>
#include <symengine/expression.h>
#include <symengine/functions.h>
#include <symengine/integer.h>
#include <symengine/matrix.h>
#include <symengine/mul.h>
#include <symengine/subs.h>
#include <symengine/symbol.h>
namespace se = ::SymEngine;
inline se::RCP<const se::Symbol> operator "" _S(const char* s, std::size_t) {
return se::symbol(s);
}
template <typename T>
inline se::RCP<const se::Basic> minus(T x) { return se::sub(se::integer(0), x); }
inline se::RCP<const se::Basic> sum() {
return se::integer(0);
}
inline se::RCP<const se::Basic> sum(const se::RCP<const se::Basic>& x) {
return x;
}
template <typename... Args>
inline se::RCP<const se::Basic> sum(const se::RCP<const se::Basic>& x, const Args&... rest) {
return add(x, sum(rest...));
}
template <typename X, std::size_t n>
inline constexpr std::size_t size(X (&)[n]) { return n; }
void run_lintest() {
using expr = se::Expression;
using std::cout;
cout << "\nrun_lintest\n===\n";
// 'state' variables
auto a = "a"_S;
auto b = "b"_S;
auto c = "c"_S;
// ode system from kinetic scheme
// a <-> b (f1, r1)
// b <-> c (f2, r2)
auto adot = sum(mul(minus("f1"_S), a), mul("r1"_S, b));
auto bdot = sum(mul("f1"_S, a), mul(minus( "r1"_S), b),
mul(minus("f2"_S), b), mul("r2"_S, c));
auto cdot = sum(mul("f2"_S, b), mul(minus( "r2"_S), c));
// rates depend on v; opaque function
auto f1_rate = se::function_symbol("f1", "v"_S);
auto r1_rate = se::function_symbol("r1", "v"_S);
auto f2_rate = se::function_symbol("f2", "v"_S);
auto r2_rate = se::function_symbol("r2", "v"_S);
// what do we get when these are subbed into adot etc?
auto rates_lintest = [&](const se::map_basic_basic& rates) {
auto adot_rate = subs(adot, rates);
auto bdot_rate = subs(bdot, rates);
auto cdot_rate = subs(cdot, rates);
cout << "adot: " << expr(adot) << "\n" << "adot_rate: " << expr(adot_rate) << "\n";
cout << "bdot: " << expr(bdot) << "\n" << "bdot_rate: " << expr(bdot_rate) << "\n";
cout << "cdot: " << expr(cdot) << "\n" << "cdot_rate: " << expr(cdot_rate) << "\n";
// is bdot linear in a, b, and c?
cout << "\n";
se::RCP<const se::Symbol> state_vars[] = {a, b, c};
for (unsigned i = 0; i < size(state_vars); ++i) {
for (unsigned j = i; j < size(state_vars); ++j) {
auto x = state_vars[i];
auto y = state_vars[j];
auto dbdot_dxdy = diff(diff(bdot_rate, x), y);
cout << "∂²bdot/∂" << expr(x) << "∂" << expr(y) << ": "
<< expr(dbdot_dxdy) << "\n";
}
}
};
cout << "\n";
se::map_basic_basic rates = {
{"f1"_S, f1_rate},
{"r1"_S, r1_rate},
{"f2"_S, f2_rate},
{"r2"_S, r2_rate}
};
rates_lintest(rates);
// tricksy not a nonlinear state dependency:
rates["f2"_S] = sum(mul(mul(sin(a),sin(a)),div(exp(sum("v"_S,b)),exp(b))),
mul(mul(cos(a),cos(a)),exp("v"_S)));
cout << "\n";
rates_lintest(rates);
// actual dependency:
rates["f2"_S] = sum(mul(mul(sin(a),sin(a)),div(exp(sum("v"_S,b)),exp(b))),
mul(mul(cos(a),sin(a)),exp("v"_S)));
cout << "\n";
rates_lintest(rates);
}
void run_jacobian() {
using std::cout;
cout << "\nrun_jacobian\n===\n";
// 'state' variables
std::size_t n = 3;
se::vec_basic x = { "a"_S, "b"_S, "c"_S };
// ODE system
auto adot = sum(mul(minus("f1"_S), "a"_S), mul("r1"_S, "b"_S));
auto bdot = sum(mul("f1"_S, "a"_S), mul(minus( "r1"_S), "b"_S),
mul(minus("f2"_S), "b"_S), mul("r2"_S, "c"_S));
auto cdot = sum(mul("f2"_S, "b"_S), mul(minus( "r2"_S), "c"_S));
se::vec_basic xdot = { adot, bdot, cdot };
se::DenseMatrix J(n, n);
jacobian(se::DenseMatrix(n, 1, xdot), se::DenseMatrix(n, 1, x), J);
std::cout << J << "\n";
// Let's solve for implicit Euler!
se::DenseMatrix I(n, n);
eye(I);
se::DenseMatrix A(n, n);
J.mul_scalar(minus("dt"_S), A);
J.add_matrix(I, A);
se::DenseMatrix xnew(n, 1);
fraction_free_LU_solve(A, se::DenseMatrix(n, 1, x), xnew);
cout << "lu solve for implicit Euler:\n" << xnew;
}
int main() {
run_lintest();
run_jacobian();
}
Output from above code:
run_lintest
===
adot: -a*f1 + b*r1
adot_rate: -a*f1(v) + b*r1(v)
bdot: a*f1 - b*f2 - b*r1 + c*r2
bdot_rate: a*f1(v) - b*f2(v) - b*r1(v) + c*r2(v)
cdot: b*f2 - c*r2
cdot_rate: b*f2(v) - c*r2(v)
∂²bdot/∂a∂a: 0
∂²bdot/∂a∂b: 0
∂²bdot/∂a∂c: 0
∂²bdot/∂b∂b: 0
∂²bdot/∂b∂c: 0
∂²bdot/∂c∂c: 0
adot: -a*f1 + b*r1
adot_rate: -a*f1(v) + b*r1(v)
bdot: a*f1 - b*f2 - b*r1 + c*r2
bdot_rate: a*f1(v) - b*r1(v) + c*r2(v) - (E**v*sin(a)**2 + E**v*cos(a)**2)*b
cdot: b*f2 - c*r2
cdot_rate: -c*r2(v) + (E**v*sin(a)**2 + E**v*cos(a)**2)*b
∂²bdot/∂a∂a: 0
∂²bdot/∂a∂b: 0
∂²bdot/∂a∂c: 0
∂²bdot/∂b∂b: 0
∂²bdot/∂b∂c: 0
∂²bdot/∂c∂c: 0
adot: -a*f1 + b*r1
adot_rate: -a*f1(v) + b*r1(v)
bdot: a*f1 - b*f2 - b*r1 + c*r2
bdot_rate: a*f1(v) - b*(E**v*sin(a)**2 + E**v*sin(a)*cos(a)) - b*r1(v) + c*r2(v)
cdot: b*f2 - c*r2
cdot_rate: b*(E**v*sin(a)**2 + E**v*sin(a)*cos(a)) - c*r2(v)
∂²bdot/∂a∂a: -b*(-2*E**v*sin(a)**2 + 2*E**v*cos(a)**2 - 4*E**v*sin(a)*cos(a))
∂²bdot/∂a∂b: E**v*sin(a)**2 - E**v*cos(a)**2 - 2*E**v*sin(a)*cos(a)
∂²bdot/∂a∂c: 0
∂²bdot/∂b∂b: 0
∂²bdot/∂b∂c: 0
∂²bdot/∂c∂c: 0
run_jacobian
===
[-f1, r1, 0]
[f1, -f2 - r1, r2]
[0, f2, -r2]
lu solve for implicit Euler:
[(a - r1*(-a*f1 + b*(1 - f1) - r2*(1 - f1)*(c*(1 - f1)*(-r1*f1 + (1 - f2 - r1)*(1 - f1)) - f2*(1 - f1)*(-a*f1 + b*(1 - f1)))/(-r2*f2*(1 - f1)**2 + (1 - f1)*(1 - r2)*(-r1*f1 + (1 - f2 - r1)*(1 - f1))))/(-r1*f1 + (1 - f2 - r1)*(1 - f1)))/(1 - f1)]
[(-a*f1 + b*(1 - f1) - r2*(1 - f1)*(c*(1 - f1)*(-r1*f1 + (1 - f2 - r1)*(1 - f1)) - f2*(1 - f1)*(-a*f1 + b*(1 - f1)))/(-r2*f2*(1 - f1)**2 + (1 - f1)*(1 - r2)*(-r1*f1 + (1 - f2 - r1)*(1 - f1))))/(-r1*f1 + (1 - f2 - r1)*(1 - f1))]
[(c*(1 - f1)*(-r1*f1 + (1 - f2 - r1)*(1 - f1)) - f2*(1 - f1)*(-a*f1 + b*(1 - f1)))/(-r2*f2*(1 - f1)**2 + (1 - f1)*(1 - r2)*(-r1*f1 + (1 - f2 - r1)*(1 - f1)))]