Mr.CAS— A Minimalistic (pure) Ruby CAS for Fast Prototyping and Code Generation

Matteo Ragni

Abstract

There are Computer Algebra System (CAS) systems on the market with complete solutions for manipulation of analytical models. But exporting a model that implements specific algorithms on specific platforms, for target languages or for particular numerical library, is often a rigid procedure that requires manual post-processing.

This work presents a Ruby library that exposes core CAS capabilities, i.e. simplification, substitution, evaluation, etc. The library aims at programmers that need to rapidly prototype and generate numerical code for different target languages, while keeping separated mathematical expression from the code generation rules, where best practices for numerical conditioning are implemented.

The library is written in pure Ruby language and is compatible with most Ruby interpreters.

Motivation and significance

Ruby (Flanagan and Matsumoto 2008) is a purely object-oriented scripting language designed in the mid-1990s by Yukihiro Matsumoto, internationally standardized since 2012 as ISO/IEC 30170.

With the advent of the Internet of Things, a compact version of the Ruby interpreter called mRuby(eMbedded Ruby) (Tanaka, Nagumanthri, and Matsumoto 2015) was published on GitHub by Matsumoto, in 2014. The new interpreter is a lightweight implementation, aimed at both low power devices and personal computers, and complies with the standard (“ISO/IEC 30170 – Information technology – Programming languages – Ruby” 2000). mRubyhas a completely new API, and it is designed to be embedded in complex projects as a front-end interface—for example, a numerical optimization suite may use mRubyfor problem definition.

The Ruby code-base exposes a large set of utilities in core and standard libraries, that can be furthermore expanded through third party libraries, or gems. Among the large number of available gems, Ruby still lacks an Automatic and Symbolic Differentiation (ASD) (Tolsma and Barton 1998) engine that handles basic computer algebra routines, compatible with all different Ruby interpreters.

Nowadays Ruby is mainly known thanks to the web-oriented Rails framework. Its expressiveness and elegance make it interesting for use in the scientific and technical field. An ASD-capable gem would prove a fundamental step in this direction, including the support for flexible code generation for high-level software—for example, IPOPT (Wächter and Laird 2009; Wächter and Biegler 2006).

Mr.CAS¹ is a gem implemented in pure Ruby that supports symbolic differentiation (SD) and fundamentals computer algebra operations (Von Zur Gathen and Gerhard 2013). The library aims at supporting programmers in rapid prototyping of numerical algorithms and in code generation, for different target languages. It permits to implement mathematical models with a clean separation between actual mathematical formulations and conditioning rules for numerical instabilities, in order to support generation of code that is more robust with respect to issues that can be introduced by specific applications. As a long-term effort, it will become a complete open-source CAS system for the standard Ruby language.

Other CAS libraries for Ruby are available:

Rucas (Lees-Miller 2010), Symbolic (Bayramgalin 2012)

: milestone gems, yet at an early stage and with discontinued development status. Both offer basic simplification routines, although they lack differentiation.

Symengine (Certik et al. 2016)

: is a wrapper of the symengine C++ library. The back-end library is very complete, but it is compatible only with the vanilla C Ruby interpreter and has several dependencies. At best of Author knowledge, the gem does not work with Ruby 2.x interpreter.

In Section [sec:description], Mr.CAScontainers and tree structure are explained in detail and applied to basic CAS tasks. In Section [sec:examples], examples on how to use the library as code generator or as interface are described. Finally, the reasons behind the implementation and the long term desired impact are depicted in Section [sec:impact]. All code listings are available at http://bit.ly/Mr_CAS_examples.

Software description

Software Architecture

Mr.CAS is an object oriented ASD gem that supports computer algebra routines such as simplifications and substitutions. When gem is required, it overloads methods of Fixnum and Float classes, making them compatible with fundamental symbolic classes.

Each symbolic expression (or operation) is the instance of an object, that inherits from a common virtual ancestor: CAS::Op. An operation encapsulates sub-operations recursively, building a tree, that is the mathematical equivalent of function composition: \[\left( f \, \circ \, g \right)\]

When a new operation is created, it is appended to the tree. The number of branches are determined by the parent container class of the current symbolic function. There are three possible containers:

CAS::Op

single sub-tree operation—e.g. \(\sin(\cdot)\).

CAS::BinaryOp

dual sub-tree operation—e.g. exponent \(x^y\)—that inherits from CAS::Op.

CAS::NaryOp

operation with arbitrary number of sub-tree—e.g. sum \(x_1 + \cdots + x_N\)—that inherits from CAS::Op.

Fig. [fig:graph] contains a graphical representation of a expression tree. The different kind of containers allows to introduce some properties—i.e. associativity and commutativity for sums and multiplications (Cohen 2003). Each container exposes the sub-tree as instance properties. Basic containers interfaces and inheritances are shown in Fig. [fig:uml-container]. For a complete overview of all classes and inheritance, please refer to software documentation.

The terminal leaves of the graph are the classes CAS::Constant, CAS::Variable and CAS::Function. The first models a simple numerical value, while the second represents an independent variable, that can be used to perform derivatives and evaluations, and the latter is a prototype of implicit functions. Those leaves exemplify only real scalar expressions, with definition of complex, vectorial, and matricial extensions as milestones for the next major release.

The symbolic differentiation (CAS::Op#diff) explores the graph until it reaches ending nodes. A terminal node is the starting point for derivatives accumulation, the mathematical equivalent of the chain rule: \[\left( f \, \circ \, g \right)' \: = \: \left( f' \, \circ \, g \right) \: g'\] The recursiveness is used also for simplifications (CAS::Op#simplify), substitutions (CAS::Op#subs), evaluations (CAS::Op#call) and code generation.

Software Functionalities

Basic Functionalities

No additional dependencies are required. The gem can be installed through the rubygems.org provider². Gem functionalities are required using the Kernel method: require ’Mr.CAS’. All methods and classes are encapsulated in the module CAS.

Symbolic Differentiation (SD) is performed with respect to independent variables (CAS::Variable) through forward accumulation, even for implicit functions. The differentiation is done by the method CAS::Op#diff, having a CAS::Variable as argument, as shown in Listing [code:example-diff].

Automatic differentiation (AD) is included as a plugin and exploits the properties of dual numbers to efficiently perform differentiation, see (Bartholomew-Biggs et al. 2000) for further details. The AD strategy is useful in case of complex expressions, where explicit derivative’s tree may exceed the call stack depth.

Simplifications are not executed automatically, after differentiation. Each node of the tree knows rules for simplify itself, and rules are called recursively, exactly like ASD. Simplifications that require a heuristic expansion of the sub-graph—i.e. some trigonometric identities—are not defined for now, but can be easily achieved through substitutions, as shown in Listing [code:example-simp].

The tree is numerically evaluated when the independent variables values are provided in a feed dictionary. The graph is reduced recursively to a single numeric value, as shown in Listing [code:example-call].

Symbolic expressions can be used to create comparative expressions, that are stored in special container classes, modeled by the ancestor CAS::Condition—for example, \(f(\cdot) \geq g(\cdot)\). This allow the definition of piecewise functions, in CAS::Piecewise. Internally, \(\max(\cdot)\) and \(\min(\cdot)\) functions are declared as operations that inherits from CAS::Piecewise—for example, \(\max(f(\cdot), g(\cdot))\). Usage is shown in Listing [code:example-expr].

Meta-programming and Code-Generation

Mr.CAS is developed explicitly for metaprogramming and code generation. Expressions can be exported as source code or used as prototypes for callable closures (the Proc object in Listing [code:example-proc]):

Compiling a closure of a tree is like making its snapshot, thus any further manipulation of the expression does not update the callable object. This drawback is balanced by the faster execution time of a Proc: when a graph needs only to be evaluated, transforming it in a closure reduces the execution time—for example, in a iterative algorithm, where a closure is called at each iteration.

Code generation should be flexible enough to export expression trees in a user’s target language. Generation methods for common languages are included in specific plugins. Users can furthermore expand exporting capabilities by writing specific exportation rules, overriding method for existing plugin, or designing their own exporter, like the one shown in Listing [code:example-exporting]:

Illustrative Examples

Code Generation as C Library

This example shows how a user of Mr.CAS can export a mathematical model as a C library. The c-opt plugin implements advanced features such as code optimization and generation of libraries.

The library example implements the model: \[\label{eq:ex1model} f(x, y) = x^y + g(x)\, \log(\sin(x^y))\] where the expression \(g(x)\) belongs to a external object, declared as g_impl, which interface is described in g_impl.h header. What should be noted is the corpus of the exported code: the intermediate operation \(x^y\) is evaluated once, even if appears twice in eq. [eq:ex1model]. The C function that implements \(f(x,y)\) is declared with the token f_impl. The exporter uses as default type double for variables and function returned values. Library created by CLib contains the code shown in Listing [code:c.ex.1].

The function \(g(x)\) models the following operation: \[g(x) = (\sqrt{x + a} - \sqrt{x}) + \sqrt{\pi + x}\] and may suffer from catastrophic numerical cancellation (Higham 2002) when the \(x\) value is considerably greater than \(a\). The user may decide to specialize code generation rules for this particular expression, stabilizing it through rationalization. Without modifying the actual model, \(g(x)\) the rationalization for differences of square roots³ is inserted into the exportation rules, as in Listing [code:example-exporting-C-2]. The rules are valid only for the current user script. For more insight about __to_c and __to_c_impl, refer to the software manual.

It should be noted the separation between the model, which does not contain stabilization, and the code generation rule. For this particular case, the code generation rule in Listing [code:example-exporting-C-2] overloads the predefined one, in order to obtain the conditioned code. Obviously, the user can decide to apply directly the conditioning on the model itself, but this may change the calculus behavior in further manipulation.

Using the module as interface

As example, an implementation of an algorithm that estimates the order of convergence for trapezoidal integration scheme (Weideman 2002) is provided, using the symbolic differentiation as interface.

Given a function \(f(x)\), the trapezoidal rule for primitive estimation for the interval \([a,\,b]\) is: \[I_{n}(a, b) = h\, \left( \dfrac{f(a) + f(b)}{2} + \sum\limits_{k = 1}^{n - 1}{f \left( a + k \,h \right)} \right)\] with \(h = (b - a) / n \), where \(n\) mediates the step size of the integration. When exact primitive \(F(x)\) is known, approximation error is: \[E[n] = F(b) - F(a) - I_{n}(a, b)\] The error has an asymptotic expansion of the form: \[E[n] \propto C\,{n}^{-p}\] where \(p\) is the convergence order. Using a different value for \(n\), for example \(2\,n\), the ratio [eq:error.ratio] takes the approximate vale: \[\label{eq:error.ratio} \dfrac{E[n]}{E[2\,n]} \approx 2^{p} \quad \rightarrow \quad p \approx \log_2 \left( \dfrac{E[n]}{E[2\,n]} \right)\] The Listings [code:example-integrate-ruby] and [code:example-integrate-python] contain the implementation of the described procedure using the proposed gem and the well known Python (Van Rossum and Drake 2011) library SymPy (Smith et al. 2016).

ODE Solver with Taylor’s series

In this example, a solving step for a specific ordinary differential equation (ODE) using Taylor’s series method (Butcher 2008) is derived. Given an ODE in the form: \[\label{eq:taylor.ode} y'(x) = f\left(x, y(x)\right)\] the integration step with order \(n\) has the form: \[y(x + h) = y(x) + h\,y'(x) + \dots + \dfrac{h^{n}}{n!} \, y^{(n)}(x) + E_{n}(x)\] where it is possible to substitute equation [eq:taylor.ode]: \[{y}^{(i)}(x) = \dfrac{\partial {y}^{(i-1)}(x)}{\partial x} + \dfrac{\partial {y}^{(i-1)}(x)}{\partial y} {y'}(x)\] For this algorithm, three methods are defined. The first evaluates the factorial, the second evaluates the list of required derivatives, and the third returns the integration step in a symbolic form. The result of the third method is transformed in a C function. In this particular case, the ODE is \(y' = x y\). For the resulting C code of Listing [code:example-ode], refer to the online version of the examples.

Other examples are available online⁴:

adding a user defined CAS::Op that implements the \(\mathrm{sign}(\cdot)\) function with the appropriate optimized C generation rule;

exporting the operation as a continuous function through overloading or substitutions;

performing a symbolic Taylor’s series;

writing an exporter for the LaTeX language;

a Newton-Raphson algorithm using automatic differentiation plugin.

Impact

Mr.CAS is a midpoint between a CAS and an ASD library. It allows one to manipulate expressions while maintaining the complete control on how the code is exported. Each rule is overloaded and applied run-time, without the need of compilation. Each user’s model may include the mathematical description, code generation rules and high level logic that should be intrinsic to such a rule—for example, exporting a Hessian as pattern instead of matrix.

Our research group is including Mr.CAS in a solver for optimal control problem with indirect methods, as interface for problems description (Biral, Bertolazzi, and Bosetti 2016).

As a long term ambitious impact, this library will become a complete CAS for Ruby language, filling the empty space reported by SciRuby for symbolic math engines.

Conclusions

This work presents a pure Ruby library that implements a minimalistics CAS with automatic and symbolic differentiation that is aimed at code generation and meta-programming. Although at an early developing stage, Mr.CAS has promising feature, some of them shown in Section [sec:examples]. Also, this is the only gem that implements symbolic manipulation for this language.

Language features and lack of dependencies simplify the use of the module as interface, extending model definition capabilities for numerical algorithms. All core functionalities and basic mathematics are defined, with the plan to include more features in next releases. Reopening a class guarantees a liquid behaviour, in which users are free to modify core methods at their needs.

Library is published in rubygems.org repository and versioned on github.com, under MIT license. It can be included easily in projects and in inline interpreter, or installed as a standalone gem.

Current code version

Code metadata (mandatory)

Nr.	Code metadata description	Please fill in this column
C1	Current code version	0.2.7
C2	Permanent link to code/repository used for this code version	github.com/ ElsevierSoftwareX/SOFTX-D-17-00013
C3	Legal Code License	MIT
C4	Code versioning system used	git (GitHub)
C5	Software code languages, tools, and services used	Ruby language
C6	Compilation requirements, operating environments	Ruby \(\geq 2.x\)
C7	If available Link to developer documentation/manual	rubydoc.info/gems/Mr.CAS
C8	Support email for questions	[email protected]