Glyph: Symbolic Regression Tools

We present Glyph – a Python package for genetic programming based symbolic regression. Glyph is designed for usage in numerical simulations as well as real world experiments. For experimentalists, glyph-remote provides a separation of tasks: a ZeroMQ interface splits the genetic programming optimization task from the evaluation of an experimental (or numerical) run. Glyph can be accessed at https://github. com/Ambrosys/glyph. Domain experts are able to employ symbolic regression in their experiments with ease, even if they are not expert programmers. The reuse potential is kept high by a generic interface design. Glyph is available on PyPI and Github.


Introduction
Symbolic regression [1] is an optimization method to find an optimal representation of a function. The method is "symbolic", because building blocks of the functions, i.e. variables, primitive functions, and operators, are represented symbolically on the computer. Genetic programming (GP) [2] can be implemented to find such a function for system identification [3,4] or fluid dynamical control [5,6]. Glyph is an effort to separate optimization method and optimization task allowing domainexperts without special programming skills to employ symbolic regression in their experiments. We adopt this separation of concerns implementing a client-server architecture; a minimal communication protocol eases its use. Throughout this paper "experiment" is meant as a synonym for any symbolic regression task including a labexperiment, a numerical simulation or data fitting. Previous work on system identification and reverse engineering of conservation laws was reported in [1,7]. As a substantial extension, we put the algorithms in the context of control problems and in contrast to these works, we publish the implementation. Modern algorithms also include multi-objective optimization [4] and advances like age fitness based genetic programming [8] or epigenetic local search [9]. There exist various approaches to the representation of multi IO problems, including stack-or graph-based representations and pointers [9,10]. As a special feature we extended our implementation such that it can be easily used for finding an optimal control law for a given, measured and actuated system. This is, however just one application of our software.

Implementation and architecture
Glyph is intended as a lightweight framework to build an application finding an optimal system representation given measurement data. The main application is intended as system control, consequently a control law is determined and returned. Glyph is built on the idea of loose coupling such that dependencies can be released if wanted.
A typical control application consists of a system and its controller, possibly separated, cf. Figure 1. Glyph has three main module to build such an application: i) the assessment, which holds all methods and data structures belonging to the experiment, ii) the GP which is responsible for the system identification, and iii) the application components, which constitute an application.

Abstraction Layers
Glyph is conceptualized around different layers of abstraction: 1 1. Representation layer: This layer consists of a datastructure representing a GP solution, i.e. an individual, and methods to manipulate such data-structures. This abstraction layer also allows for an interchangeable representation. 2. Algorithm layer: This layers encapsulates the selection and tweaking of individuals.
3. AssessmentRunner layer: Abstracts calculation of cost functions. 4. GPRunner layer: This is an integration layer which lets the user iterate the evolutionary algorithm calculating statistics after each iteration. 5. Application layer: This layer sets the context and control flow.
The first three layers are sufficient to formulate and optimize a symbolic regression problem. The GPRunner and Application let you gap the bridge between prototypes and real-world-problem-solving applications.

Building an Application
An application consists of a GP callable, the gp_runner, an assessment callable for input, the assessment_ runner, and the application which uses both of these classes and holds all application-relevant details. A command-line application is built by assessment_runner = AssessmentRunner (assess_args) gp_runner = glyph.application.GPRunner (gp_args) app = command_line_application (app_args) The assessment_runner has one argument, the parallel_factory which implements a map() method, possibly parallel. For an application one needs to implement setup, assign_fitness, and measure: setup is self-explaining, measure is a key method which takes as input a set of measurement functions and combines them into a tuple of callable measures for multiobjective optimization. The measures are used eventually in assign_fitness where the return values are used to assign a fitness to an individual from GP. The interface is freely extensible. A gp_runner forwards the evolutionary iteration. It takes as arguments, gp_ args, an individual class, a gp_algorithm, and an assessment_runner. The individual class contains the representation of a function, the individual; it is based on DEAP's tree-based implementation. The gp_ algorithm takes care for the breeding and selection steps, its principles are described in [2].
The application is run in the main function with app. run().
In the application and gp_runner, the user has freedom to add functionality using the list of callbacks in the arguments, say, to implement other logging or streaming options. This allows for very flexible programming. We constructed the components that way to allow users to specialize for their particular experiments and possibly increase performance or extend the symbolic regression, e.g. by replacing the DEAP tree-based representation of an individual.

Remote Control
One main objective of glyph is its use in a real experiment. In this case, the GP loop is separated from the experimental loop in a client-server setup using ZeroMQ [11], cf. Figure 2.
Consequently, one should implement the interface to the experiment using the protocol described in the following subsection. Having the implementation of the experiment, the server, one needs to implement the client, i.e. the interface to the gp_runner. In essence this means connecting the correct sockets with ZeroMQ and ensuring that the gp_runner and the assessment_runner use the corresponding sockets. Then, the main application is assembled as before, now using a RemoteApp for the main application, which in turn uses a gp_runner, which then uses now a RemoteAssessmentRunner. That is it, we can run remotely our GP evaluation from some client and the experiment in place of the experiment.

Communication Protocol
The communication is encoded in json [12]. A message is a json object with two members: { "action": "value", "payload": "value" } The possible values are listed in Table 1. The config action is performed prior to the evolutionary loop. Entering The corresponding control law a = u(s, t) is determined by symbolic regression. Right: gp-based symbolic regression finds different candidate control laws. Each candidate solution is given a fitness score Γ which is used to compare different solutions and to advance the search in function space. Figure adapted from [5] with permission.
s a Γ the loop, discovered solutions will be batched and a experiment action will be requested. You can configure optional caching for re-discovered solutions. This includes persistent caching between different runs. The shutdown action will let the experiment program know that the GP loop is finished and you can safely stop the hardware.
Configuration settings are sent as a json object in key:value form, where the keys contain the option to be set, there is only one mandatory option: the primitive set. To configure the primitive set, the primitive names are passed as content of the key config, whose values specify the corresponding arities, both fields described again as json object.
The experiment action sends a list of expressions, encoded as string in prefix (also: polish) notation [13]. For each expression sent, the experiment returns a fitness tuple.
Additionally, one can define the type of algorithm, error metric, representation, hyper-parameters, etc. A comprehensive up to date list can be found at http:// glyph.readthedocs.io/en/latest/usr/glyph_remote/.

Application example: control of the chaotic Lorenz System
In the following, we demonstrate the application and use of Glyph by the determination of an unknown optimal control law for a chaotic system. As an example, we study the control of the potentially chaotic Lorenz system. Chaotic systems are very hard to predict and control in practice due to their sensitivity towards small changes in the initial state which may lead to exponential divergence of trajectories. The Lorenz model [14] consists of a system of three ordinary differential equations: with two nonlinearities, xy and xz. Here x, y, and z make up the system state and s, r, b are parameters: s is the Prandtl number, r is the Rayleigh number, and b is related to the aspect ratio of the air rolls. For a certain choice of parameters and initial conditions chaotic behavior emerges. We present two examples where the target is to learn control of bring a chaotic Lorenz system to a complete stop, that is, (x, y, z) = 0 (t ∈ R). In the first example, "control in y", the actuator term is applied to ẏ. This allows for a more direct control of the system, since y appears in every equation of (1) and, thus, influence all three state components, x, y, and z. In the second example, "control in z", the actuator term is applied to ż, which leads to a more indirect control, since the flow of information from z to x is only through y.
The system setup is summarized in Tables 2 and 3. When r = 28, s =10, and b = 8/3, the Lorenz system produces chaotic solutions (not all solutions are chaotic). Almost all initial points will tend to an invariant setthe Lorenz attractor -a strange attractor and a fractal. When plotted the chaotic trajectory of the Lorenz system resembles a butterfly (blue graph in Figure 3). The target of control is, again, formulated as RMSE 2 of the system state with respect to zero (separately for each component). The control function u can make use of ideal measurements of the state components. Additionally, we allow for a single symbolic constant k to be used as argument for  SHUTDOWN --the control law. Given an expression, this constant is optimized separately, see also [4]. The respective GP runs for control in y and control in z are conducted with the corresponding random seeds labeled "in y" and "in z".  Table 4. The wide spread of the RMSE values is a sign of conflicting objectives that are hard to satisfy in conjunction. Interestingly, almost all solutions, u, commonly introduce a negative growth rate into ẏ. This effectively drives y to zero and suppresses the growth terms, sy and xy, in the equations for ẋ and ż respectively, in turn, driving x and z to zero as well. As would be    Table 3: Control of the Lorenz system: system setup. expected, minimal expressions, of length 1 or 2, cannot compete in terms of the RMSE. For example, the simple solution, u(x, y, z) = -ky (fourth row), is almost as good as the lengthier one, u(x, y, z) = -exp(x) + ky (first row), and even better in RMSE y . Table 4 shows the results from the GP run. One solution immediately stands out: u = k · x + z, with k = -27.84 (second row). It is exactly what one might expect as a control term for the chaotic Lorenz system with control in y. This control law effectively reduces the Rayleigh number r to a value close to zero (k ≈ r), pushing the Lorenz system past the first pitchfork bifurcation, at r = 1, back into the stable-origin regime. If r < 1 then there is only one equilibrium point, which is at the origin. This point corresponds to no convection. All orbits converge to the origin, which is a global attractor, when r < 1.

Control in
The phase portrait of the solution from the first and second row of Table 4 are illustrated in Figure 3. After a short excursion in negative y direction (t ≈ 5), the green trajectory quickly converges to zero. The red trajectory seems to take a shorter path in phase space, but, it is actually slower to converge to the origin. This is verified by a plot of the trajectories for the separate dimensions x, y and z over time Figure 4. Control in z: For control in z the actuator term u is added to the left side of the equation for ż in the uncontrolled system (1)

( )
, , x y bz u x y z = -+ ż Selected Pareto-front individuals from the GP run are displayed in Table 5. As mentioned at the beginning of this section, effective control is hindered by the indirect influence of z on the other state variables, hence, it is not surprising that the control laws here are more involved than in the previous case. Also, they generally do not perform well in the control of z, which is expressed by the relatively high values in RMSE z . This is confirmed by the phase portrait of the solution u(x, y, z) = -(k · (-y) + x · z + y + z) shown in figure Figure 5: While going straight to the origin in the xy-plane there are strong oscillations of the trajectory along the z-axis. The dynamics caused by the actuation, e.g. for the best control law found, can be explained qualitatively: there  is a strong damping in all variables but y. This reflects the tendency to suppress z-oscillations and, at the same time, to add damping in y through the xz term: if y grows, the z contribution to damping on the right hand side of the Lorenz equations (1) grows and, in turn, damps y. This is, however, only possible to some extent, hence, the oscillations observed in figure Figure 5.
We conclude the demonstration with a short summary: Using Glyph we can find complex control laws, even for unknown systems. This cannot be easily achieved with other frameworks. The control laws found can be studied analytically in contrast to several other methods which have black-box character. The usage is straightforward, as we have described above. The above example can be found online as an example.

Other symbolic regression libraries
Due to its popularity, symbolic regression is implemented by most genetic programming libraries. A semi-curated list can be found at http://geneticprogramming.com/ software/. In contrast to other implementations, Glyph implements higher concepts, such as symbolic constant optimization, and also offers parallel execution for complex examples (control simulation, system identification). This is very important for practical applications. Glyph is well tested, cf. Table 6 and currently applied in two  experiments and several numerical problems. For control, apart from our implementation, there exists only one more alternative as a dedicated matlab toolbox (with python interface), openMLC [15], which contains much of the material treated in [6]. Apart from the fact that matlab itself is no free software, one difference is that a toolbox is not as powerful as a library made for developers to extend into various directions. Our first application is control, however the library is built in an extensible way and now more and different applications are ongoing. Further, we offer multi-objectivity, multi-output and constant optimization which make symbolic regression applicable and useful for practical tasks.

(3) Reuse potential
The potential to use Glyph is twofold: on one hand applications can be easily written and the elegant core functionality can be extended; on the other hand, researchers can use the code as core for symbolic regression and extend its functionality in a very generic way. With respect to applications, currently two main directions are targeted: modeling using genetic programming-based symbolic regression and the control of complex system, where a control law can be found generically, using genetic programming.
The detailed examples and tutorial allow usage from beginner to experienced level, i.e. undergraduate research projects to faculty research. The design of Glyph is such that generic interfaces are provided allowing for very flexible extension.