RUBI integrator infrastructure

[certik]

@naveensaigit and I brainstormed steps that could be done to get the infrastructure for rule rewriting and decision tree generation working in SymEngine. Our goal is to start with a few rules, but get them working all the way, and then keep adding more rules.

• Last status update on RUBI: https://github.com/symengine/symengine/issues/1597#issuecomment-521721092

Here are the steps that we could do:

1. Take 3 rules from SymPy from here: https://github.com/sympy/sympy/blob/d7a854fb789591114d7d51c3f23bba2fd4d520cc/sympy/integrals/rubi/rubi.py#L71 (for example pick some rules from this file: https://github.com/sympy/sympy/blob/d7a854fb789591114d7d51c3f23bba2fd4d520cc/sympy/integrals/rubi/rules/piecewise_linear.py), replace ManyToOneReplacer with ManyToOneMatcher and create a simple script. Use the following file: https://github.com/HPAC/matchpy/blob/5cae3f275e3a1f725518516bf3c700dc3be03a56/tests/test_code_generation.py as an example how to use the ManyToOneMatcher class. Update: it seems one can use ManyToOneReplacer also.

2. Generate a Python decision tree with those 3 rules using https://github.com/HPAC/matchpy/blob/a90cc4684c85f6afea0a6386d2338374eda7e1d9/matchpy/matching/code_generation.py or using https://github.com/Upabjojr/rubi_generated

3. Test that it works.

4. Document well into SymPy (and add tests) how to generate a working decision tree for a user defined set of rules. This will be a nice capability to have, well documented with a tutorial.

5. Generate C++ decision tree with those 3 rules using https://github.com/symengine/symengine/blob/ed7479baf59a36636061acc35b676bcf9073a932/symengine/utilities/matchpycpp/cpp_code_generation.py

6. Test that it works, perhaps some utility functions might need to be ported from SymPy: https://github.com/sympy/sympy/blob/d7a854fb789591114d7d51c3f23bba2fd4d520cc/sympy/integrals/rubi/utility_function.py

7. Document the process well, just like in 4., but for C++ / SymEngine.

[certik]

The branch w and this file: https://github.com/certik/rubi_generated/blob/79613660c77db95736b1d4eb6fbfb29081230364/run_example.py contains a test for a few simple rules that seem to work.

The branch rubi_test and this file: https://github.com/certik/sympy/blob/0773a99bc991742790c033fbc108db579787f8b5/sympy/integrals/rubi/rules/integrand_simplification.py#L138 contains our custom 3 rules.

This finishes the steps 1., 2. and 3. The step 4. (documentation) we can do later.

Going forward, one option is to continue with the step 5.

An alternative path forward is to try to simply the generated decision tree, which currently looks like this: https://github.com/certik/rubi_generated/blob/79613660c77db95736b1d4eb6fbfb29081230364/gen/generated_part000001.py, and make it look like this: https://github.com/certik/sympy/blob/3b34a308887e611add9b95acde4b6b5010b859f2/sympy/integrals/rubi/SineIntegrationRules.py. So in the simple example above, matchpy would be used to match the pattern Integral(x_**WC("n", S(1))*WC("b", S(1)), x_) with not constraints, and then the following function would be called:

def replacement1(b, n, x):
    if (ZeroQ(n+S(1))):      
        return b*log(x)
    else:
        return x**(n+1)*b/(n+1)

or more specifically, it should look like this:

def generated_rule(b, n, x):
    if (ZeroQ(n+S(1))):      
        return replacement2(b, n, x)
    else:
        return replacement1(b, n, x)

or actually even like this:

def generated_rule(b, n, x):
    if cons3(b, x) and cons4(n, x):
        if cons_a2(n):      
            return replacement2(b, n, x)
        elif cons_a1(n):
            return replacement1(b, n, x)
        else:
           raise Exception("Unhandled case")
    else:
        raise Exception("Another unhandled case")

And later we can introduce a Not class in our rules, so that instead of using cons_a2() we can just do Not(cons_a1()), and from that the generator can easily directly output:

def generated_rule(b, n, x):
    if cons3(b, x) and cons4(n, x):
        if cons_a1(n):
            return replacement1(b, n, x)
        else:
            return replacement2(b, n, x)
    else:
        raise Exception("Another unhandled case")

[certik]

@Upabjojr, do you know if the decision tree generated could be modified to output something like:

def generated_rule(b, n, x):
    if (ZeroQ(n+S(1))):      
        return replacement2(b, n, x)
    else:
        return replacement1(b, n, x)

or rather something like:

def generated_rule(b, n, x):
    if cons3(b, x) and cons4(n, x):
        if cons_a1(n):
            return replacement1(b, n, x)
        else:
            return replacement2(b, n, x)
    else:
        raise Exception("Another unhandled case")

instead of the current output, such as in https://github.com/certik/rubi_generated/blob/79613660c77db95736b1d4eb6fbfb29081230364/gen/generated_part000001.py, which seems really complex and probably slower?

[Upabjojr]

@Upabjojr, do you know if the decision tree generated could be modified to output something like:

Do you mean the decision tree generated by MatchPy? The code follows the exact way in which MatchPy operates internally, so it would be quite tough to customize the generated code.

[certik]

Yes, the tree generated by the file https://github.com/Upabjojr/rubi_generated/blob/cd35e9e51722b04fb159ada3d5811d62a423e429/decision_tree_rubi_generator.py. I thought you wrote it.

[Upabjojr]

It's subclassing MatchPy's CodeGenerator object in order to do some customizations of the generated code.

def generated_rule(b, n, x):
    if (ZeroQ(n+S(1))):      
        return replacement2(b, n, x)
    else:
        return replacement1(b, n, x)

Where is this code coming from?

[certik]

That's the code that I would like to get generated. If it can be generated, together with the pattern, then it will greatly simplify the decision tree and also probably make it faster.

[Upabjojr]

Well, I'm not sure it can be generated like that because it is missing all the logic to match the elements of the expression tree.

[Upabjojr]

By the way, I feel that the decision tree is being generated more complicated than what it ought be. Incidentally, I believe it's an issue of replication of objects in the rule definition.

The rules may not be optimal when they are added to MatchPy.

[certik]

@Upabjojr would you have time to do a video call with @naveensaigit and I? We have been working on this lately. We think it should be possible to generate, but we have to write the generator, either using MatchPy or even from scratch.

[Upabjojr]

@Upabjojr would you have time to do a video call with @naveensaigit and I? We have been working on this lately. We think it should be possible to generate, but we have to write the generator, either using MatchPy or even from scratch.

Yes, sure!

[Upabjojr]

Please remember that unless there is support for the coroutines introduced in C++20, any attempt at using RUBI in SymEngine will either timeout or cause a memory error.

[rikardn]

There is a repo setup for the upcoming Rubi 5 here: https://github.com/RuleBasedIntegration/Rubi-5

This version of Rubi will require virtually no pattern matching (only 42 rules as compared to 3000 for Rubi-4) and would be both faster to run and easier to implement in symengine. What is currently available is a prototype for a small subset of rules. The author of Rubi has called for help to work on Rubi 5.

Perhaps the 42 patterns could be manually translated into pattern checking functions in symengine and the underlying Intnnn functions that use no pattern matching could be automatically translated into symengine C++.

I know you guys have been working hard on the Rubi-4, but perhaps it would be worthwhile to work on Rubi-5 in parallel? I might have some time to start working on Rubi-5 in that case.

[rikardn]

I implemented the first two patterns of Rubi-5 in https://github.com/symengine/symengine/pull/1801.

[Upabjojr]

Maybe we need a way to automatically translate the code into C++.

[Upabjojr]

(only 42 rules as compared to 3000 for Rubi-4)

Have you tried to use the submodule matchPyCpp on those 42 rules? It's a Python submodule in SymEngine that is able to transform Mathematica pattern matching into a structure of nested if-/for-loops in C++.

[Upabjojr]

Also, as a reminder, this work could be useful to SymPy as well...

[rikardn]

Thanks for your input @Upabjojr! It's a good idea to try matchPyCpp. Since it was "only" 42 rules I was think they could be translated manually, but I might be naïve. The thousands of if/then/else functions must still be machine translated. This translation could perhaps also generate code for SymPy.

(We also have to keep in mind that Rubi-5 is far from ready.)