Hi Andrew,
Il 04/12/23 02:02, Andrew Simper ha scritto:
I chose to base my symbolic system on how I can see Mathematica
represents things: an n-ary expression tree, which is working out well
for me. I'm using a sorted n-ary tree structure where each node has a
"head" saying what the node is and a set of links to other nodes. I
use these fundamental mathematical operations: Add, Times, Power, but
derive subtraction and division from them: a - b = add(a, times(-1
,b)), a / b = times(a, power(b, -1)), which really helps when doing
simplifications, and is how Mathematica does things. Another powerful
concept in Mathematica are "rules" for pattern matching and
expressions manipulation, for example a_*b_ + a_*c_ -> a_*(b_ + c_),
where the a_ means "any expressions" and something like a_Integer, or
a_Symbol just means a single integer or variable respectively. Finder
matching sub expressions, including negative ones, is still something
I want to optimise further.
I've got something similar: n-ary trees where each node is an operation
or a constant/variable, the basic operations are sum (+ and -),
multiplication (* and /), power, and "function" which represents any
other kind of operation specialized via JavaScript's prototype
inheritance. Sums and multiplications can have an arbitrary number of
operands. There's a bunch of very specific details on why I chose this
specific representation but, TBH, I don't remember everything (I
developed this stuff around end of 2020 and didn't touch it since then)
- what I do remember is that this allows to efficiently solve linear
equations and to check for equality of expressions at least.
I remember I also initially tried implementing a pattern matching
mechanism to specify substitutions, but IIRC I ditched the idea as it
would have been really complex and it would be hard anyway to understand
which substitutions would need to applied when in the first place. And
indeed...
I think it's an interesting optimization problem to solve these
systems efficiently, that from my point of view looks NP-Hard, much
like the traveling salesman problem. The problem statement goes
something like: "Solve A.x = z for x (x and z being vectors of length
n, and A being an n by n matrix) in the least "cost"
in fundamental operations =, +, -, *, / where the cost of operations
are zero for =, one for + - *, and some higher number for eg 5 for /.
The resultant solution will usually can be written as a series of
assignments of temporary variables to (possibly common)
sub-expressions of the final solution, ending with assignments to the
n elements of x being the solution to the system of linear equations.
... that is what happens. :-)
In practice, I noticed that for the application at hand (solving the
linear part of circuit models):
1. matrix triangulation via Gauss elimination seems to reduce the number
of operations in total compared to e.g. substitution (in any case, you
get one division per row, which is always the worst offender and in most
cases unavoidable);
2. choosing the order of rows is extremely important w.r.t. the number
of operations you actually end up performing, especially in case you
have several parameters - that is why the software I built let's the
user play with such order freely.
Nailing these two aspects leads to good enough results IMO, without
getting into automatic expression optimization (you can always do that
by hand afterwards if really needed).
I just coupled all of this with some rudimentary a SPICE netlist parser
with the possibility to have expressions as component values (including
unknown parameters), some code to automatically derive circuit equations
and "categorize" them in different sets, a rudimentary HTML GUI, a
flexible discretization approach [1], and Ciaramella
(https://ciaramella.dev/) code generation from the triangulated matrix,
and voila! :-)
I'm surprised, given the usefulness, that there aren't more papers
published in this area, perhaps I just don't know what to look for!
I also couldn't find much back then. I know there were a couple of books
on the matter, one of which too basic to be useful, so I just did it
from scratch my way. Maybe one could go and study Maxima's source code
or something.
Best,
Stefano
[1] F. Germain, K. Werner, "Design Principles for Lumped Model
Discretisation Using Mobius Transforms", DAFx-15.
