Ideas to improve the performance of Mathics

[mmatera]

During the last month, I was doing some exploration about how to improve one of the biggest issues in the project: the slowness of the evaluation process, compared against WMA and other implementations of WL.

Let's consider a "simple" example:

Do[ll=Join@@Table[{a,b,c},{i,n}];Timing[Plus@@ll][[1]]//Print[n,"->",#1]&,{n,{1,10,100,1000,10000}}]

This is, check how much it takes to evaluate the sum
a+b+c+a+b+c+....
with 3n terms.

In WMA, we get (in my Intel i7, 16GB laptop)

n	WMA	Mathics
1	->4.*^-6	0.000286495
10	->8.*^-6	0.0012344
100	->0.000022	0.007596
1000	->0.000231	0.0675576
10000	->0.003397	0.974039

in a way that a simple addition takes ~ 300 times more time in Mathics than in WMA.

Or, in a more sophisticated example
'Table[Sin[i+Sqrt[2.]],{i,100.}]//Timing//#1[[1]]&'
WMA: 0.00007
mathics: 0.085582

with a difference of ~1200 times or, for a matrix:

Table[Sin[i-j+Sqrt[2.]],{i,100.},{j,100.}]//Timing//#1[[1]]&

WMA: 0.001593
mathics: 10.1881

the factor is around 6800 times.

From these "simple" examples, it follows that the difference between both implementations is not linear, and hence, we can not just blame to the slowness of Python against C, but to something hidden in the evaluation algorithm.

We can expect that using tools like Cython, or JIT could reduce an overall factor (maybe one order of magnitud) in the evaluation time, but not to correct the (polynomial?) scaling. Then, to solve the slowness of Mathics, we need to review more carefully different aspects of the algorithm:

• How does an expression is "walked" inside its leaves be evaluated.
• How the definition of the symbols involved in a function is accessed (to read properties, and rules)
• How, once a definition is found, are the associated rules recovered.
• How a pattern is evaluated to match against a (sub) expression, and how the corresponding rule is applied.
• How well are defined the rules in order to converge to a fixed point expression after a certain number of iterations.
• How new expressions are built using the leaves of the previous one.
• How are implemented all the caching strategies in the middle.

and probably many other important details, ranging from memory management to specific algorithms for certain frequently used symbols (like arithmetic, logic, etc) and the interface to libraries like sympy, mpmath, or numpy.

Ideally, we would like to pick one microscopic aspect of this process at a time, and try to improve it. However, many of these parts are actually quite entangled, and hence, modify one aspect forces to implement a lot of other changes in order to make it works (just at the level of passing the CI tests).
Here I would like to share the ideas I had, how and where were implemented, and what I have learned doing it.

[mmatera]

Symbols, definitions and rules

In WL, each symbol has a definition, that includes different kind of rules (ownvalues, downvalues, upvalues, formats) and values (attributes, options, etc) linked to a "fully qualified name" (FQN) in the form contextname`symbolname. For WL, two symbols with the same FQN are equivalent, so it seems that would be natural to avoid to instantiate two symbols with the same name. However, until the last release, it was perfectly possible to have different instances of the same symbol. For this reason, until that release, to compare two symbols required to access the name of both of them, and compare the strings. Now, we have "singletonized" the Symbol class. This allows to decide if two symbols are the same symbol just by comparing they memory positions (through the is keyword).

In the current implementation, the definition of a symbol is stored in a Definitions object. So, the same instance of the kernel can hold different instances of Definitions. This always have been in this way since mathics 1.0.0 (at least). This architecture implies that obtaining the definition of a symbol requires to look into a dictionary of definitions for a definition that matches with the symbol. Actually, in the current implementation, the Definitions object stores several dictionaries:

• Definitions.builtin, Definitions.user and Definitions.pymathics dictionaries, to store the builtin rules, the user (defined at runtime) rules and the pymathics rules.
• A Definitions.definitions_cache
• A Definitions.lookup_cache

All these dictionaries have the name of a symbol as the key. A definition is then accessed by calling Definitions.get_definition. This method first tries to find the (not necessarily FQN, but possibly the bare name of the symbol) in the definitions_cache dictionary. If it is successful, it returns it, that is what happens in most of the cases. If not, it looks for candidates in builtin, user, and pymathics dictionaries. If it finds more than one candidate, it merges the three definitions into a new definition, which is stored in the definitions_cache and returned. If it does not find one rule, it creates it and stores in the definitions_cache, under the (possibly not ) FQN required, and in the user dictionary.
I guess that much of this overhead could be avoided if each symbol would store its own definition, but this would require to change quite a lot the Mathics API.

By now, the best I can imagine to improve this was implemented in definition3 (PR #59). In this branch, the call to get_definition was avoided always possible in favor of accessing directly to the definitions_cache.

[mmatera]

Following the same line, #64 tries to accelerate the access to the attributes of a symbol, again avoiding calling get_definition. I included also there a shortcut to access the name of the Context and ContextPath (own)values: these values are used each time a non-FQN is passed to the get_definition method. In the current master, each time a non-FQN is asked, both the current Context and the contexts in ContextPath are examined to look for a symbol with the required name. To do this, get_definition is called to get the definition of Context (and ContextPath), then their ownvalues are picked up, and from them, the rule, and finally, the value. In #64, variables in Definitions are stored in a way that always mirrors the values in the corresponding definitions.

[rocky]

Good observations and a good start. Let's start to drill down from a high-level observation into some concrete actions or next steps that will help us address this.

Let me start with:

Table[Sin[i+Sqrt[2.]],{i,100.}]

Actually, the first question I have is: is this the simplest example that shows a general problem worth solving? Would we see dramatic discrepancies still if we removed the outer Sin[]? If I reduced the bound from 100 to 50 would the time be also cut in half? If I am writing tests, having an equivalent test that takes half the time is preferable, especially when the test is in say the range of seconds. And when I am debugging and checking symbolic expansion, I'd rather test using a small number of items if it captures the same idea.

The first step in fixing this is to isolate and pinpoint problems. The simpler and narrower the problem, the better - if Table[i + Sqrt[2.]], {1, 25}] can be improved on its own, we can then extend to look into including Sin[], or maybe what additionally needs to be done going from say 25 to 100 if that is a separate problem needing a separate solution.

Let us assume that we can get dramatic improvements just in the simpler Table[i+Sqrt[2.]],{1,25}] if we did better.

Now let me break down one aspect of this particular expression. In Python or a programming language like that you would write:

[i + sqrt(2.0) for i in range(1, 100)]

And one thing that you notice is that sqrt(2.0) is a constant which needs to be computed only once. So this is the same thing:

j = sqrt(2.0)
[i + j for i in range(1, 100)]

In Python the improvement might not be so great. But in Mathics, since things like evaluating Sqrt takes a while, the linear speedup in evaluating sqrt(2.0) can be large.

In compiler optimization technology this kind of improving code transformation is called "constant propagation" and "moving invariant code out of loops". (here the invariant code is an expression) And in sequences (e.g.
[1 + sqrt(2.0), 2 + sqrt(2.0), 3 + sqrt(2.0), ...]) this kind of thing happens all of the time .

One would hope that a symbolic manipulation system like SymPy already understands that this is is a problem and looks for invariant code and constant propagation.

So a question I have for someone is: Does Sympy understand this issue and handle it in some way? If so, we should try to hook into that.

If not, we would suggest that Sympy look into it and possibly write our own routines for constant propagation and moving invariant expressions out of loops. Especially when functions are "compiled".

[mmatera]

Actually, the first question I have is: is this the simplest example that shows a general problem worth solving? Would we see dramatic discrepancies still if we removed the outer Sin[]? If I reduced the bound from 100 to 50 would the time be also cut in half? If I am

Regarding this question, the choice of the tests maybe was not optimized, but these are the reason to choose them instead of something simpler:

• WMA and Mathics use different cache strategies. For this reason, I want to avoid both running the same expression more than once. Choosing functions that produce irrational numbers is a way to avoid repetitions.
• I am trying to compare two systems, one which is two or three orders of magnitude in runtime, so something that takes 1 second in Mathics, could take 1ms in WMA. For a matter or resolution, I need something that I can resolve in WMA, and does not take forever in Mathics.
• Again regarding the runtime of the tests and cache. Mathics-benchmark runs the same test many times, to measure times of the order of the second. The tests that I propose are designed to take that time (in Mathics) but avoiding to be cacheable. Another possibility would be to run the same test in many instances (let say, call many times mathics -e ....with...` a short test) but then I would expend more time initializing Mathics than carrying out the test.
• A drawback of my approach of working with large expressions is that the low-level memory management starts to entangle with the higher-level aspects of the algorithm. I do not have clear how to disentangle these two effects.

Finally, just to be clear, usually, I show in my comments the simpler and more representative test that I was able to find to communicate an issue/idea. This does not mean that I haven't run several other simpler and more complex tests before picking the one I show.

[mmatera]

Regarding the second part, the problem with symbolic evaluation is that the same expression can be evaluated to different values in different places. For example:

a=1; {F[a], a=a+1; F[a]} returns {Sin[1], Sin[2]}, not {Sin[1],Sin[1]}. Still, more diabolic things could happened. Consider this:

F[x_]=a*x; a=1; {F[2], a=2*a;F[2]}
resulting in {2, 4}
you have here the same expression with a numeric argument. However, it depends implicitly on a variable. These kinds of things are perfectly valid in WL, and make the constant propagation optimization not feasible. I think that the best it can be done is to use some caching method, and indeed it seems to be implemented in Mathics (look at mathics.core.expression.ExpressionCache and related functions). The problem is that I do not still fully understand how is it working.

Regarding sympy, I think that we can do is to try a similar test over pure sympy to see how does it behave and compare it against Mathics.

Finally, something that is shown to help in expressions like Table[N[i+Sqrt[2]],{i,10}] is to numberify it by evaluating N[i+Sqrt[2]] keeping i unevaluated, and then build the table, instead of evaluate directly {N[1+Sqrt[2]], N[2+Sqrt[2]],...}

[rocky]

Responding to this really deserves its own discussion. So see #76 for information on this.

As for the other kinds of things here, over the next few days, I will go over all of this and maybe we can start to drill down into those.

The end goal here might be to split out what each of us is doing so we don't all work in the same area and in cross purposes.

[rocky]

Let's consider a "simple" example:

Do[ll=Join@@Table[{a,b,c},{i,n}];Timing[Plus@@ll][[1]]//Print[n,"->",#1]&,{n,{1,10,100,1000,10000}}]

This is, check how much it takes to evaluate the sum a+b+c+a+b+c+.... with 3n terms.

If this is simple, then would you please write the corresponding thing in sympy, and add to the table what this speed is? This would narrow the slowness between that may exist between Mathics implementation, and Sympy Python and Mathematica.

And/or do likewise for other examples.

I thought I would try this, but when I come down to it, this is too complicated and/or we haven't narrowed things enough.

Yes, writing examples delineating precise boundaries is a bit effort, possibly hard, and maybe requires a bit more skill. However not doing so means flailing around.

I see that there was work in mathics-benchmark to try this, but those examples don't look like they have any kind of sympy interaction at all. Possibly they only show how bad the Mathics Python implementation is for Rule replacement and pattern matching.

And the problem here, is that in say compilation it doesn't matter how bad that is if it is done only one time and subsequently we have a representation that doesn't need to go through those steps anymore. (There would be other representations and other steps that would happen instead.

[rocky]

Symbols, definitions and rules

In WL, each symbol has a definition, that includes different kind of rules (ownvalues, downvalues, upvalues, formats) and values (attributes, options, etc) linked to a "fully qualified name" (FQN) in the form contextname`symbolname. For WL, two symbols with the same FQN are equivalent, so it seems that would be natural to avoid to instantiate two symbols with the same name. However, until the last release, it was perfectly possible to have different instances of the same symbol. For this reason, until that release, to compare two symbols required to access the name of both of them, and compare the strings. Now, we have "singletonized" the Symbol class. This allows to decide if two symbols are the same symbol just by comparing they memory positions (through the is keyword).

One little thing that has been nagging me has been the choice of the name InstanceableBuiltin; "Instanceable" is not a word, and it is confusing.

If I understand correctly, the class InstanceableBuiltin is supposed to implement that Singleton Design pattern, which is a very common design pattern. See for example https://www.tutorialspoint.com/python_design_patterns/python_design_patterns_singleton.htm where you can read descriptions. Google for "instanceable design pattern" and you'll find the same nothing that appears in the docstring for this class.

If that's the case then SingletonBuiltin in my opinion be clearer along with, more importantly, a description of what the intent is.

[mmatera]

OK, but InstanteableBuiltin is just the opposite to a singleton class: Buitlin classes are supposed to be instantiated just once, to store the information about builtin symbols, and objects of the class InstanceableBuiltin can be instantiated many times, like Pattern objects, or GraphicBox objects. Notice also that (except for a typo - I guess-) this class was there before I start getting involved in the project.

[rocky]

Is this a the Factory Design Pattern then?

[mmatera]

Update: If I didn't understand it wrong, "Factory Design Pattern" is similar to the implementation of mathics.core.Pattern. You have an abstract class mathics.core.Pattern, and then a function mathics.core.Pattern_create that creates an instance of a subclass depending on the argument.
On the other hand, InstanceableBuiltin are "Builtins" (classes that represent the definition of a given symbol) that can be instantiated to represent an expression, like "GraphicsBox". In some sense, this is similar to how arithmetic operations are implemented in Sympy: In Sympy, "a+b" is nor represented as an Expression (sympy.core.expr.Expr) but as an object of an specific class sympy.core.add.Add. I do not know if this obeys any standard design pattern.

[TiagoCavalcante]

@mmatera do you think Instanceable builtins could be made singleton?

[TiagoCavalcante]

@rocky @mmatera I think we need to move to a better importing system. Now it is done through a file search and checking every var from the module to see if that is a builtin.

My idea is to use a decorator to do that (otherwise we would need to have a file with hundreds of imports), something like the following:

def contribute(cls):
    cls.contribute()

@contribute
class N(Builtin):
    ...

[TiagoCavalcante]

@mmatera why it shouldn't be called by the decorator?

My guess is that the "contribution" should be done to the symbol instead of to the Definitions collection. But this would change the current implementation in a not backward compatible way

I don't understand this too, it contributes to the definitions so we can assess it from a dict, how it would contribute to the symbol?

[mmatera]

Because the "builtin" definition is the same for a given symbol in any instance of the Definitions object. So, when Definitions needs to lookup a builtin definition, it could look into the Symbol's dictionary instead of its own dictionary.

[TiagoCavalcante]

@mmatera yes, that makes sense. The attributes change is also breaking, so we could accumulate changes to 5.0.0, by now maybe release a new Docker image, @rocky?

[TiagoCavalcante]

Because the "builtin" definition is the same for a given symbol in any instance of the Definitions object. So, when Definitions needs to lookup a builtin definition, it could look into the Symbol's dictionary instead of its own dictionary.

@mmatera if you do this I can make the faster loading system . Currently I don't understand very well where all builtins need to be, seems MakeBoxes and the operators are different from the other builtins.

[mmatera]

Builtins are classes that could perfectly never be instantiated: Are just a way to put together a symbol name, and some rules in the form of WL rules and Python code. InstanceableBuiltins are classes that can be instantiated during evaluations, and behave like Expressions, but also contribute to Definitions, associated to the head of the expression. For example, a GraphicsBox and PatternObjects are InstanceableBuiltins

[mmatera]

The attribute stuff is easy to make backward compatible : just requieres an extra 'if ' in the contribute method.

[TiagoCavalcante]

@mmatera it is not that easy, the routines of getting, setting and clearing attributes are very different.

[mmatera]

No, I mean, in the definition of a Builtin, you can specify the attributes as a list of strings. Then, when the contribute method process its class attributes to build a Definition, you can do the translation.

[TiagoCavalcante]

@mmatera there are 5 Python Mathics libraries as I know, so I could update them all instead of the translation layer.