5.5 Applying a function on subexpressions

Sometimes you may want to perform an operation on specific parts of an expression while leaving the general structure of it intact. An example of this would be a matrix trace operation: the trace of a sum is the sum of the traces of the individual terms. That is, the trace should map on the sum, by applying itself to each of the sum's operands. It is possible to do this manually which usually results in code like this:

ex calc_trace(ex e)
{
    if (is_a<matrix>(e))
        return ex_to<matrix>(e).trace();
    else if (is_a<add>(e)) {
        ex sum = 0;
        for (size_t i=0; i<e.nops(); i++)
            sum += calc_trace(e.op(i));
        return sum;
    } else if (is_a<mul>)(e)) {
        ...
    } else {
        ...
    }
}

This is, however, slightly inefficient (if the sum is very large it can take a long time to add the terms one-by-one), and its applicability is limited to a rather small class of expressions. If calc_trace() is called with a relation or a list as its argument, you will probably want the trace to be taken on both sides of the relation or of all elements of the list.

GiNaC offers the map() method to aid in the implementation of such operations:

ex ex::map(map_function & f) const;
ex ex::map(ex (*f)(const ex & e)) const;

In the first (preferred) form, map() takes a function object that is subclassed from the map_function class. In the second form, it takes a pointer to a function that accepts and returns an expression. map() constructs a new expression of the same type, applying the specified function on all subexpressions (in the sense of op()), non-recursively.

The use of a function object makes it possible to supply more arguments to the function that is being mapped, or to keep local state information. The map_function class declares a virtual function call operator that you can overload. Here is a sample implementation of calc_trace() that uses map() in a recursive fashion:

struct calc_trace : public map_function {
    ex operator()(const ex &e)
    {
        if (is_a<matrix>(e))
            return ex_to<matrix>(e).trace();
        else if (is_a<mul>(e)) {
            ...
        } else
            return e.map(*this);
    }
};

This function object could then be used like this:

{
    ex M = ... // expression with matrices
    calc_trace do_trace;
    ex tr = do_trace(M);
}

Here is another example for you to meditate over. It removes quadratic terms in a variable from an expanded polynomial:

struct map_rem_quad : public map_function {
    ex var;
    map_rem_quad(const ex & var_) : var(var_) {}

    ex operator()(const ex & e)
    {
        if (is_a<add>(e) || is_a<mul>(e))
     	    return e.map(*this);
        else if (is_a<power>(e) && 
                 e.op(0).is_equal(var) && e.op(1).info(info_flags::even))
            return 0;
        else
            return e;
    }
};

...

{
    symbol x("x"), y("y");

    ex e;
    for (int i=0; i<8; i++)
        e += pow(x, i) * pow(y, 8-i) * (i+1);
    cout << e << endl;
     // -> 4*y^5*x^3+5*y^4*x^4+8*y*x^7+7*y^2*x^6+2*y^7*x+6*y^3*x^5+3*y^6*x^2+y^8

    map_rem_quad rem_quad(x);
    cout << rem_quad(e) << endl;
     // -> 4*y^5*x^3+8*y*x^7+2*y^7*x+6*y^3*x^5+y^8
}

ginsh offers a slightly different implementation of map() that allows applying algebraic functions to operands. The second argument to map() is an expression containing the wildcard ‘$0’ which acts as the placeholder for the operands:

> map(a*b,sin($0));
sin(a)*sin(b)
> map(a+2*b,sin($0));
sin(a)+sin(2*b)
> map({a,b,c},$0^2+$0);
{a^2+a,b^2+b,c^2+c}

Note that it is only possible to use algebraic functions in the second argument. You can not use functions like ‘diff()’, ‘op()’, ‘subs()’ etc. because these are evaluated immediately:

> map({a,b,c},diff($0,a));
{0,0,0}
  This is because "diff($0,a)" evaluates to "0", so the command is equivalent
  to "map({a,b,c},0)".