Numerical evaluation

Symbolica support efficient numerical evaluations of expressions and code generation. We will go through various alternatives and mention their relative advantages.

Single evaluation

To evaluate an expression once, use evaluate. In Python this evaluates with float and in Rust the type of the evaluation is inferred (one could evaluate with 4x4 simd floats or arbitrary-precision floats for example).

Python

from symbolica import *
x, f, cos = Expression.symbol('x', 'f', 'cos')
e = cos(x)*3 + f(x,2)
print(e.evaluate({x: 1}, {f: lambda args: args[0]+args[1]}))
print(e.evaluate_complex({x: 1+2j}, {f: lambda args: args[0]+args[1]}))

Rust

use symbolica::{atom::AtomCore, parse, symbol};

fn main() {
    let a = parse!("cos(x)*3 + f(x,2)").unwrap();
    let mut h = HashMap::default();
    h.insert(symbol!("x"), 1.);
    let mut fn_map = HashMap::default();
    fn_map.insert(symbol!("f"), Box::new(|x, y| x + y));
    println!("{}", a.evaluate(&h, &fn_map));

    let mut h = HashMap::default();
    h.insert(symbol!("x"), Complex::new(1., 2.));
    let mut fn_map = HashMap::default();
    fn_map.insert(symbol!("f"), Box::new(|x, y| x + y));
    println!("{}", a.evaluate(&h, &fn_map));
}

The first argument of evaluate is a map from an atom to a constant and the second argument is a map of a function symbol to a closure that is called with the evaluated arguments.

The above example yields

4.620906917604419
(9.098169021058997-7.155693397455401j)

Evaluators

Naively evaluating an expression can be prohibitively slow for large expressions, especially when it has to be evaluated many times, such as for Monte Carlo simulations/integrations. Symbolica provides state-of-the-art ways to rewrite an expression in a form that makes it very fast to evaluate, often 10-100 times faster than a naive evaluation. Note that during these substitutions, “fast math” operations are applied.

It tries to find the best multivariate Horner scheme, performs common subexpression elimination and recycles intermediate variables and registers. It resembles what a programming language compiler such as gcc would do, but it is much faster as the input is more specific. See this blog post for more information.

Nested expressions

Such an evaluator supports nested expressions; expression that call other expressions. For example, to evaluate \(e\) in

\[ \begin{align*} e &= x + \pi + \cos(x) + f(g(x+1),2x)\\ f(y,z) &=y^2 + z^2 y^2\\ g(y) &=x+y+5 \end{align*} \]

we use the following code:

Python

from symbolica import *
x, y, z, pi, f, g = Expression.symbol(
    'x', 'y', 'z', 'pi', 'f', 'g')

e1 = Expression.parse("x + pi + cos(x) + f(g(x+1),x*2)")
fd = Expression.parse("y^2 + z^2*y^2")
gd = Expression.parse("x + y + 5")

ev = e1.evaluator({pi: Expression.num(22)/7},
             {(f, "f", (y, z)): fd,
              (g, "g", (y, )): gd},
            [x])
res = ev.evaluate([[5.]])
print(res)

Note that at the moment no floats may appear in the expressions or in the map. Use rationalize_coefficients() to convert them.

Rust

fn main() {
    let e1 = parse!("x + pi + cos(x) + f(g(x+1),x*2)").unwrap();
    let f = parse!("y^2 + z^2*y^2").unwrap();
    let g = parse!("x+y+5").unwrap();

    let mut fn_map = FunctionMap::new();

    fn_map.add_constant(
        Atom::new_var(symbol!("pi")),
        Rational::from((22, 7)).into(),
    );
    fn_map
        .add_function(
            symbol!("f"),
            "f".to_string(),
            vec![symbol!("y"), symbol!("z")],
            f.as_view(),
        )
        .unwrap();
    fn_map
        .add_function(
            symbol!("g"),
            "g".to_string(),
            vec![symbol!("y")],
            g.as_view(),
        )
        .unwrap();

    let params = vec![parse!("x").unwrap()];

    let evaluator = e1
        .evaluator(&fn_map, &params, OptimizationSettings::default())
        .unwrap();

    let mut e_f64 = evaluator.map_coeff(&|x| x.into());
    let r = e_f64.evaluate_single(&[5.]);
    println!("{}", r);
}

Which yields 25864.42651932832.

The evaluator takes a hashmap of functions with the following key:

(function symbol, function name, argument list)

where the function_name should be a C++-compatible name.

Multiple expressions

Multiple expressions can be turned into an evaluator at the same time, which means they will be optimized together. As a result, common subexpressions between the expressions will be computed only once.

Python

from symbolica import *
x = Expression.symbol('x')

e1 = x * x + 5
e2 = x * x + 6

ev = Expression.evaluator_multiple([e1, e2], {}, {}, [x])
res = ev.evaluate([[5.]])
print(res)

Rust

fn main() {
    let x = parse!("x").unwrap();
    let e1 = &x * &x + 5;
    let e2 = &x * &x + 6;

    let evaluator = Atom::evaluator_multiple(
        &[e1.as_view(), e2.as_view()],
        &FunctionMap::new(),
        &[x],
        OptimizationSettings::default(),
    )
    .unwrap();

    let mut e_f64 = evaluator.map_coeff(&|x| x.into());
    let mut out = vec![0., 0.];
    e_f64.evaluate(&[5.], &mut out);
    println!("{:?}", out);
}

Yields 30.0, 31.0.

Branching

Conditional functions can be registered and takes three arguments: the first is an expression which is tested for not being equal to 0. The second argument is the true branch and the third argument is the false branch. Since the arguments are arbitrary expressions, they can contain additional branching. Below we register the user function if:

Python

ev = E("if(y, x*x + z*z + x*z*z, x * x + 3)").evaluator({}, {}, [S("x"), S("y"), S("z")], conditionals=[S("if")])
print(ev.evaluate([[3., 1., 2.]])) # 25.0
print(ev.evaluate([[3., 0., 2.]])) # 12.0

Rust

let mut f = FunctionMap::new();
f.add_conditional(symbol!("if")).unwrap();

let mut eval = parse!("if(y, x*x + z*z + x*z*z, x * x + 3)")
    .evaluator(
        &f,
        &vec![crate::parse!("x"), crate::parse!("y"), crate::parse!("z")],
        Default::default(),
    )
    .unwrap()
    .map_coeff(&|x| x.re.to_f64());

println!("{}", val.evaluate_single(&[3., 1., 2.])); // 25.0
println!("{}", val.evaluate_single(&[3., 0., 2.])); // 12.0        

Code generation

Evaluators can be turned into highly optimized C++ code, optionally with inline assembly (x86_64 and arm64). The latter will drastically cut down on C++ compilation time (and may even yield faster code), effectively making the compilation time independent on the size of the expression.

The compiled library accepts evaluations with double and complex<double> and can be directly loaded into Symbolica.

Python

comp = eval.compile('my_fun', 'test.cpp', 'test.so')
res = comp.evaluate([[5.]])
print(res)

Rust

let mut compiled = e_f64
    .export_cpp("test.cpp", "my_fun", true, InlineASM::X64)
    .unwrap()
    .compile("test.so", CompileOptions::default())
    .unwrap()
    .load()
    .unwrap();

let mut out = vec![0.];
compiled.evaluate(&[5.], &mut out);
println!("{}", out[0]);

There is support for SIMD and CUDA generation as well (see the documentation of compile).