/* This file:
 *   http://anggtwu.net/MAXIMA/mysubst.mac.html
 *   http://anggtwu.net/MAXIMA/mysubst.mac
 *          (find-angg "MAXIMA/mysubst.mac")
 * Author: Eduardo Ochs <eduardoochs@gmail.com>
 *
 * (defun e () (interactive) (find-angg "MAXIMA/mysubst.mac"))
 * See: (find-es "maxima" "mysubst")
 *
 * Idea: this is the chain rule,
 *
 *          [ d                          ]
 *   [RC] = [ -- f(g(x)) = f'(g(x))g'(x) ]
 *          [ dx                         ]
 *
 * and this is I how write - in mathematical notation - the
 * substitution that specializes it to a certain particular case:
 *
 *        [  f(u) := sin(u) ]   [ d                             ]
 *   [RC] [ f'(u) := cos(u) ] = [ -- sin(42*(x)) = cos(42*x)*42 ]
 *        [  g(x) := 42*x   ]   [ dx                            ]
 *        [ g'(x) := 42     ]
 *
 * We can't use items like "f(u)=sin(u)" in substitution lists in
 * Maxima, but everything works if we convert them to lambdas, like
 * this:
 *
 *   RC     : ( 'diff(f(g(x)),x) = fp(g(x))*gp(x) );
 *   substs : [  f = lambda([u], sin(u)),
 *              fp = lambda([u], cos(u)), 
 *               g = lambda([x], 42*x),
 *              gp = lambda([x], 42)     ];
 *   subst(substs, RC);
 *
 * The functions in this file implement a preprocessor for
 * substitution lists that convert substitutions declared with `:='s
 * into their lambda forms, and doesn't change the substituions
 * declared with `='s. For example, in
 *
 *   substs : '[ a = 42,
 *               h(x,y) := 10*x+y ];
 *   mysubst_ify(substs);
 *   mysubst    (substs, h(a,3));
 *
 * the results of the "mysubst_ify" and of the "mysubst" are:
 *
 *   (%i16) mysubst_ify(substs);
 *   (%o16)            [a = 42, h = lambda([x, y], y + 10 x)]
 *   (%i17) mysubst    (substs, h(a,3));
 *   (%o17)                              423
*/

mysubst_f(fxye) := block([fxy,f,xy,e],
    fxy : lhs(fxye),
    f   : op(fxy),
    xy  : args(fxy),
    e   : rhs(fxye),
  buildq([f,xy,e], f=lambda([splice(xy)],e)))$

mysubst_1  (ab)      := if is(op(ab) = ":=") then mysubst_f(ab) else ab$
mysubst_ify(fxyes)   := map('mysubst_1, fxyes)$
mysubst    (fxyes,o) := subst(mysubst_ify(fxyes), o)$


/*
 * Tests:

* (eepitch-maxima)
* (eepitch-kill)
* (eepitch-maxima)
load("~/MAXIMA/mysubst.mac");
gradef(f(u), fp(u));
gradef(g(x), gp(x));
RC     : ('diff(f(g(x)), x) = diff(f(g(x)), x));
RC     : ('diff(f(g(x)), x) = fp(g(x)) * gp(x));
substs : '[  f(u) := sin(u),
            fp(u) := cos(u),
             g(x) := 42*x,
            gp(x) := 42     ];
mysubst_ify(substs);
mysubst    (substs, RC);

* (eepitch-maxima)
* (eepitch-kill)
* (eepitch-maxima)
load("~/MAXIMA/mysubst.mac");
mysubst_f('(f(x,y):=10*x+y));

            '[f(x,y):=10*x+y, g(x):=10*x, a=42];
mysubst_ify('[f(x,y):=10*x+y, g(x):=10*x, a=42]);

gradef(g(x), g_x(x));
gradef(F(u), f(u));

BB : F(g(x));
CC : F(u);
aa : diff(BB, x);
dd : diff(CC, u);
AA : integrate(aa, x);
DD : integrate(dd, u);

MM : align_eqs([AA,BB,CC,DD]);

substs2 : '[g(x):=2*x, g_x(x):=2];
substs4 : '[g(x):=2*x, g_x(x):=2, F(u):=sin(u), f(u):=cos(u)];
mysubst_ify(substs2);
mysubst_ify(substs4);
MM2     : mysubst(substs2, MM);
MM2     : mysubst(substs4, MM);
MM3     : ev(MM2, 'integrate);
MM4     : subst([u=2*x], MM3);

*/