refactor arithmetic power

[mmatera]

This is another chuck of #766. This part is related to the way in which powers are evaluated.

Part of this code is dealing with the fact that the output from Sympy is not compatible with the way in which WL handles expressions involving powers. Tests were fixed here according to the behavior in WMA.

[mmatera]

@rocky, notice that this part was one of the most polemical in the discussion of #766. The reasons I have for insisting in this are:

• Provides a closer behavior to the WMA, so I can build tests that are compatible with its outputs.
• Despite the implementation of eval_Power_can be improved, I think that splitting it from the symbol definition would help to improve it in the future.
• In particular, this could help to think about how to improve the sympy back and forward conversions, if we still want to use it for these arithmetic evaluations.

[rocky]

@rocky, notice that this part was one of the most polemical in the discussion of #766. The reasons I have for insisting in this are:

• Provides a closer behavior to the WMA, so I can build tests that are compatible with its outputs.
• Despite the implementation of eval_Power_can be improved, I think that splitting it from the symbol definition would help to improve it in the future.
• In particular, this could help to think about how to improve the sympy back and forward conversions, if we still want to use it for these arithmetic evaluations.

Ok - noted. I note this and will try to take a look at over the weekend.

[rocky]

In master, this rule is not needed.: Power[I, -1] returns -I.

In this branch though I am seeing (-I) ^ (-1 / 2). Which one matches WMA. And why the discrepancy?

[mmatera]

The point is that certain rules implemented implicitly by Sympy are not available now because I wrote a different implementation that does not rely too much on Sympy. At this moment I do not remember exactly which test fails if we remove the new rule, but it is easy to check.

[rocky]

I do not consider your different implementation as a good thing but a bad thing. Given the choice of relying on Sympy code or mmatera code we should go with Sympy. It is more used, more maintained ,and generally more well thought out and more complete.

Our job should be to figure out how to get WMA things done in Sympy when it can.

[rocky]

I now see that SymPy has a "Sign" function that works over reals and complexes.

I missed this fact before. We should be using that.

[mmatera]

sympy.Sign works like WL Sign and, even if it is true to force sympy.Sign as WL RealSign, it would mean doing an extra work in conversion, by checking that the expression is a Real expression. This can be done by hacking to_sympy to handle RealSign in a special way, but that work is more or less equivalent to performing the actual computation in pure Python, with the overhead of the back-and-forth conversion.

For now, please you should not worry about the overhead of back-and-forth conversion this way.

I think we can reduce conversions without custom code and actually down the line we will get a better result if we can just have sympy work on larger chunks of sympy code.

[mmatera]

BTW, Sign[3 - 4 I] can not be evaluated by using only sympy. The problems of using sympy.sign are shown in the tests by replacing the implementation of mathics.eval.arithmetic.eval_Sign by

def eval_Sign(expr: BaseElement) -> Optional[BaseElement]:
    """
    if expr is a number, return its sign.
    """
    sympy_expr = expr.to_sympy()
    if sympy_expr is not None:
        print("   sympy_expr", sympy_expr)
        result = from_sympy(sympy.sign(sympy_expr))
        if result is not None:
            return result
    return None

In a similar way, if we have expressions with products of negative rational powers, the work of convert back the product to Mathics is more or less the same than evaluating the product without sympy.

[rocky]

Here is what I see when I evaluate the Python/Sympy equivalent of Sign[3 - 4I]:

>>> from sympy import *
>>> sign(3 - 4j)
0.6 - 0.8*I
>>> 

It looks like Complex.to_sympy() is not doing the right thing.

[rocky]

In this branch I am not seeing parenthesis around the result. I assume WMA adds parenthesis?

I wonder if this is a formatting problem or a precedence problem

[mmatera]

In WMA

In[1]:= a b DirectedInfinity[-I]                                                

Out[1]= a b (-I Infinity)

In[2]:= %//FullForm                                                             

Out[2]//FullForm= Times[a, b, DirectedInfinity[Complex[0, -1]]]

On the other hand, in master,

In[1]:= a b DirectedInfinity[-I]
Out[1]= a b -I Infinity

In[2]:= %//FullForm
Out[2]//FullForm= Times[a, b, DirectedInfinity[Complex[0, -1]]]

Notice that the first output is wrong, and indeed, is a formatting problem.

[rocky]

Thanks for the information. It looks like this is an operator precedence problem.

[mmatera]

OK, I see the problem now: in this branch, the test is marked as "skip". I solved this precedence issue in another branch, and for some reason, the change was lost. I will try to put this fix and tests in another PR.

[rocky]

Same thing here.

[mmatera]

In any case, the test seems to be right. The problem is that if we activate it, it does not pass it. I can revert the change, but it will not have any effect until the underlying problem gets fixed, and the tests get activated.

[rocky]

Hmm. master returns -1 right now, so what we have is wrong.

Do you know what part of the code fixes this?

[mmatera]

I guess this is related to the new rules you mention before...

[mmatera]

Nop. The rule at the end seems to be useless. The computation is done by mathics.eval.arithmetic.eval_Power_number, which now is called from Power.eval_check ( mathics.builtin.arithfns.basic, line 591).

[rocky]

Why the change in test?

[mmatera]

This is the effect of implementing the associativity of powers. In master,

Sqrt[1 / (1 + 1 / (1 + 1 / a))]//FullForm

produces

Power[Power[Plus[1, Power[Plus[1, Power[a, -1]], -1]], -1], Rational[1, 2]]

while here,

Power[Plus[1, Power[Plus[1, Power[a, -1]], -1]], Rational[-1, 2]]

On the other hand, looking at WMA, it seems that it does not implement the associativity in this case either.

[rocky]

Something is weird here in terminology.

In normal terminology, a rule is applied if it matches. period. It has nothing to do with whether things change. The WMA evaluation process has decided to look for fixed points for its own reasons.

[mmatera]

OK, but this is not how Mathics builtin rules work. I guess that the reason to do this is to allow to use simpler patterns in BuiltinRules. For example, suppose you want to define a rule of the form

F[x___/;cond]:=(evaluation...)

Sometimes, to evaluate cond is expensive, and can be done as a part of the evaluation process. Then, to avoid duplicating an expensive computation, in a BuiltinRule we can define
F[x___]:=(evaluation...)
and eventually, if in the evaluation we determine that x1,x2,... does not satisfy the condition, we can return None, so the evaluator understands that the rule didn't match.

A different situation would be to have a rule like F[x___]->(eval) that for a certain input F[x1,x2,...] must return the same function (i.e. is invariant). The evaluator indeed can rebuild the expression, but as it is a new expression, it is very expensive to determine (from outside) that both expressions are the same, and we reach a fixed point. Returning None avoids the need of this comparison.

[mmatera]

To simplify the discussion of this PR, I removed the comment, so we can go back to this in another round.

[rocky]

I don't understand this comment. Based on past experience, I have a sinking feeling this is going off into some misguided direction. (I think of the days when text substitutions on MathML boxes was proposed as a way of addressing formatting?)

[mmatera]

eval_ methods are used as the function attribute of a BuiltinRule. When the rule is applied, the eval_ method is filled with the parameters from the pattern, and produce an output. If the output is None, the rule is assumed to "fail" and the original expression is kept. Otherwise, if the eval_ method returns an Element, then this element is used as the result of the replacement, and the rule is considered "applied".

Here we have certain eval_ functions that are not eval_ methods but can be used as a part of them. If you are up, please help me to rephrase the comment.

[mmatera]

An alternative way to solve the same problem would be to add a parameter expression into the eval_ methods. Then, we can just check if the result of the replacement is the same as the original expression without any extra cost. However, it will take to change the code in many many places...

[mmatera]

To simplify this PR, I am going to remove these comments, that still are available in #766.

[rocky]

Are we sure this kind of thing can't be done in Sympy. Have we tried this?

For another thing brought up and in a previous PR, with Sign, I see in that a problem there mentioned in StackOverflow was just about not tagging a variable as being in the right type (complex vs real).

[mmatera]

The point is not if we can find a way to use Sympy, but if we can do it "natively", without the cost of back-and-forth conversions. Even if we can make these conversions more efficient, avoiding them at all would always be a win, mostly if the native implementation is not too complicated.

[mmatera]

See comment above.

[rocky]

It feels to me there is more here than just

@rocky, notice that this part was one of the most polemical in the discussion of #766. The reasons I have for insisting in this are:

• Provides a closer behavior to the WMA, so I can build tests that are compatible with its outputs.

I would like more detail here. Be specific about what behavior is not close. Elsewhere I see specific differences in output even with this PR. I think we need to understand in better detail what it is that is making the behavior different.

• Despite the implementation of eval_Power_can be improved, I think that splitting it from the symbol definition would help to improve it in the future.

I am not opposed in theory to having something that can let us control how power works. But again, I'd appreciate understanding exactly what. Here I see the associativity of turning repeated powers into multiplications. But I am not convinced Sympy can't do this as well. Or that we couldn't do this via a rule.

• In particular, this could help to think about how to improve the sympy back and forward conversions, if we still want to use it for these arithmetic evaluations.

Until a good cases are made for not using sympy for arithmetic evaluation let's use sympy for arithmetic evaluation. In this PR I see that for I^(2/3) where we get the wrong result, but it is not clear to me that this is the fault of Sympy versus how we feed stuff to Sympy.

So let's take some specific examples like this see what's up here for specific examples.

[rocky]

This code still seems to have bits of code that aren't associated with arithmetic power. Please let's split off anything that is not related to the association stuff and put that in a different PR so we can disentangle whatever issues those have from the association part.