series: support Piecewise expressions in mrv()

[skirpichev]

No description provided.

[skirpichev]

@anutosh491 do you still think, that sympy/sympy#18492 is solved in the SymPy? ;-)

[skirpichev]

I sense I could share the link with you here , so that you could maybe give it a quick read.

Thanks, I did. But I don't know how much I can help you: my knowledge of the current state of the SymPy's series module (and related things) is very limited. In my opinion, the limit() function and series() helpers tend to be more and more complex in the SymPy. (BTW, as a first point, I would suggest you to change pictures to the code snippets, if it possible. In this way I can better understand which example really poses a problem for the Diofant.)

But maybe I'm wrong. For instance, you can try to compare tests for limit() for the SymPy vs the Diofant. It's not too hard to extract them with the ast module, here I can help you. That you can compare with the code complexity. E.g.:

$ wc -l diofant/series/gruntz.py diofant/series/limits.py 
  352 diofant/series/gruntz.py
  214 diofant/series/limits.py
  566 total

vs

$ wc -l ../sympy/sympy/series/gruntz.py ../sympy/sympy/series/limits.py 
  719 ../sympy/sympy/series/gruntz.py
  369 ../sympy/sympy/series/limits.py
 1088 total

I can guess, that your chances to be accepted will slightly increase if you find the SymPy - more successful.

see if some of these things can be improved in diofant

Feel free to open a bug if you have found any issue. There are definitely a lot of things to be improved (and a related project).

[anutosh491]

Thanks @skirpichev , whenever I see any bug in sympy, the first thing I usually do is to try it out in diofant and I have been shocked multiple times when a bug (related to the series) module has worked perfectly fine in diofant and has performed poorly in sympy . Hence I'll keep this in mind( to open any issues if spotted) for the future .

Btw, I have a pr trying to implement the leading term method for Piecewise functions, which addresses all these cases perfectly and hence there's no real need to do it through gruntz if that goes in . Though there is some minor confusion regarding handling bi-directional cases , I am somewhat happy with how it performs .
Pr - sympy/sympy#22555
Comment explaining logic behind implementation - sympy/sympy#22555 (comment)

The tests are quite extensive so that it caters to any form of piecewise expression that can be formed . Maybe you could have a look at it and let me know if I can improve something . Could be added to diofant too if you like , currently I find no failing piecewise limit through the code as I ended up trying multiple cases !!

[skirpichev]

I have a pr trying to implement the leading term method for Piecewise functions, which addresses all these cases perfectly and hence there's no real need to do it through gruntz if that goes in .

A good example why the implementation of the limit() in SymPy tends to be more and more complex...

Maybe you could have a look at it and let me know if I can improve something.

As I said before, I don't know how limit() capabilities are implemented in the SymPy now.

Your code seems to be too complex. Does it handle more Piecewise limits than the Diofant?

Could be added to diofant too if you like

Yeah, I should implement support for arbitrary boolean expressions in limit(), i.e. And(x > 0, x < 1). Also, the implementation of limits for relational expressions is wrong in the Diofant. E.g. limit(x > 0, x, 0) should be True.

[anutosh491]

I have a pr trying to implement the leading term method for Piecewise functions, which addresses all these cases perfectly and hence there's no real need to do it through gruntz if that goes in .

A good example why the implementation of the limit() in SymPy tends to be more and more complex...

Correct correct , I've tried to explain a bit regarding why it ended up becoming complex .

Your code seems to be too complex. Does it handle more Piecewise limits than the Diofant?

I assume yes , cause I see that code answering any sort of Piecewise expression , the test go as follows(shown below) - and it took me good amount of time to answer all sort of limits for all these functions. Like i've shown below there are about 6 ways to represent the same Piecewise function(or rather the meaning behind it) in sympy currently , hence to address all these limits, it took some careful observation from my side and eventually the code became complex.

Hence either we standardize how we declare peicewise functions ( maybe moving from -oo to oo, like first the left most piece and moving to the right most piece) or end with complicated code.

def test_piecewise_basic2():
    from sympy.logic.boolalg import (And, Or)
    p0 = Piecewise((0, Or(x < -2, x > 2)), (1, True))
    p1 = Piecewise((0, x < -2), (0, x > 2), (1, True))
    p2 = Piecewise((0, x > 2), (0, x < -2), (1, True))
    p3 = Piecewise((0, x < -2), (1, x < 2), (0, True))
    p4 = Piecewise((0, x > 2), (1, x > -2), (0, True))
    p5 = Piecewise((1, And(-2< x, x < 2)), (0, True))
    Functions = [p0, p1, p2, p3, p4, p5]
    for function in Functions:
        assert limit(function, x, -oo) == 0
        assert limit(function, x, -2, '-') == 0
        assert limit(function, x, -2, '+') == 1
        raises(ValueError, lambda: limit(function, x, -2, dir='+-'))
        assert limit(function, x, 0, '+-') == 1
        assert limit(function, x, 2, '-') == 1
        assert limit(function, x, 2, '+') == 0
        raises(ValueError, lambda: limit(function, x, 2, dir='+-'))
        assert limit(function, x, oo) == 0


def test_piecewise_basic3():
    p0 = Piecewise((0, x < -2), (2, x > 2), (1, True))
    p1 = Piecewise((2, x > 2), (0, x < -2), (1, True))
    p2 = Piecewise((0, x < -2), (1, x < 2), (2, True))
    p3 = Piecewise((2, x > 2), (1, x > -2), (0, True))
    Functions = [p0, p1, p2, p3]
    for function in Functions:
        assert limit(function, x, -oo) == 0
        assert limit(function, x, -2, '-') == 0
        assert limit(function, x, -2, '+') == 1
        raises(ValueError, lambda: limit(function, x, -2, dir='+-'))
        assert limit(function, x, 0, '+-') == 1
        assert limit(function, x, 2, '-') == 1
        assert limit(function, x, 2, '+') == 2
        raises(ValueError, lambda: limit(function, x, 2, dir='+-'))
        assert limit(function, x, oo) == 2

[skirpichev]

I assume yes

I'll look then. But I would appreciate, if you point to a specific example (except for above notes about BooleanFunction's and #1217).

[skirpichev]

BTW, now I've an implementation for Boolean/Relational limits (#1218), which supports all mentioned above your tests. Few lines less.

[anutosh491]

BTW, now I've an implementation for Boolean/Relational limits (#1218), which supports all mentioned above your tests. Few lines less.

That's nice to hear, I'll give it a look sometime soon . Would like to know better ways to tackle these sort of limits !

[anutosh491]

@skirpichev I had written to you on gitter few days back , continuing our discussion on the performance of the series module in diofant and sympy . I tried to reach out to you there so that I do not to disturb the purpose of the pr which our discussion would end up doing eventually .
Would be glad to know whether you're active on gitter or I could write to you here itself ?
I wanted to talk about the implementation of order term for multivariate expressions in sympy and diofant ! and some other things too . Would be glad to hear your views on the issues i would like to present !

[skirpichev]

Thank you for the report, I've deleted my gitter account.

Feel free to reach me in any public space of the Diofant project. If you see an issue - report, please.

[anutosh491]

Thanks @skirpichev for the reply , I had some doubts regarding implementation of order based on multivariate expressions in sympy and diofant . There were some comments pointed on a pr of mine by jksoum (didn't use @ to prevent unnecessary reference)
He writes
1)
I'm not sure that O(x + y) should be x + y. It seems to imply that x + y could be regarded as a single variable while O(1/(x + y)) seems to give O(1/x, x, y).
2)
We also have

In [51]: O(x - y)
Out[51]: O(x + y; (x, y) → (0, 0))

That makes no sense, x - y and x + y can take any values independently of each other.

3)
It is not clear how Order should be defined in case of several variables. I would expect that examples like this

>>> O(x + 2*y)
O(x + y, x, y)

would mean that |x + 2y| < M|x + y| for some constant M, at least when x and y are close to 0, but that is obviously not true: x + 2*y need not be 0 when x + y is 0 (except when both variables are 0). This type of definition might be meaningful but probably impossible to implement in general.

I think that the only feasible definition in multivariable case is the one where leading terms are taken in a definite order. First, the leading term of f(x, y) wrt. x would be something like x**n*g(y) and then the leading term wrt. y would be constant times x**n*y**m. I don't know how useful that would be but at least it could be implemented once the leading term methods are in place.

The points raised here are quite legit and would I be great if I could look into these and address these during Gsoc . But I am not clear regarding what should be the probable answers then
What should be the output for the following O(x + 2*y) , If the outputs are clearer obviously I could address the coding part better .
I had also asked as a reply whether returning self here would be good enough to which he writes

I don't think so. The result should rather be the one coming from taking leading terms in given variable order.

Hence I would like to know your views on this and how could we maybe incorporate what he says ,cause the things being pointed out are quite valid and I feel should be addressed !

[skirpichev]

I'm not sure that O(x + y) should be x + y.

I see here

In [6]: O(x + y)
Out[6]: O(x + y; (x, y) → (0, 0))

What else do you expect?

O(x - y) ... O(1/(x + y)) ... O(x + 2*y)

All this seems to be wrong, indeed. See also sympy/sympy#6822

It is not clear how Order should be defined in case of several variables.

Well, there is a reference (wikipedia article). f(x)=O(g(x)) as x->p iff there exists M: |f(x)|<=M |g(x)| in a some neighbourhood of x=p. Most code for multivariate O probably doesn't use this definition. For most expressions, I doubt there are too much possible correct (and simple!) transformations for the O constructor: you can drop a constant for the Mul expression, or take a leading term if one factor is an univariate expression... Probably, it's fine to use original expressions in all examples above.

BTW, I'm not sure, that the multivariate O notation is useful at all.

[skirpichev]

FYI: #1228

[anutosh491]

Thanks @skirpichev , I feel maybe we could go with the sequential leading term protocol that we have been following .
Which says
O(x + 2*y, x, y) == (x + 2*y).as_leading_term(x).as_leading_term(y) == (2*y).as_leading_term(y) == 2*y . This is a concept which felt more interesting than other many other issues and would surely like to solve this if there is some good direction to proceed here .

[anutosh491]

Another issue I had but can't comprehend of a solution is something which diofant currently lacks implementation for . It is based on the order term and the methods(based on order) used in the is_convergent block .
Consider the following code (This should be answered by dirichlet's test which is also broken so I would expect a NotImplementedError if this first step is fixed in sympy)

        order = O(sequence_term, (sym, S.Infinity))

        ### --------- p-series test (1/n**p) ---------- ###
        p_series_test = order.expr.match(sym**p)
        if p_series_test is not None:
            if p_series_test[p] < -1:
                return S.true
            if p_series_test[p] >= -1:
                return S.false

        ### ------------- comparison test ------------- ###
        # 1/(n**p*log(n)**q*log(log(n))**r) comparison
        n_log_test = order.expr.match(1/(sym**p*log(sym)**q*log(log(sym))**r))
        if n_log_test is not None:
            if (n_log_test[p] > 1 or
                (n_log_test[p] == 1 and n_log_test[q] > 1) or
                (n_log_test[p] == n_log_test[q] == 1 and n_log_test[r] > 1)):
                    return S.true
            return S.false

Now Sum(sin(n)/n, (n, 1, oo)) is convergent but when we calculate order here , we have O(sin(n)/n , (n, oo)) which is nothing but O(1/n, (n, oo)) , hence this will never pass the p-series or the comparison or limit comparison tests in the first place ( because what simply happens in p series tests in based on power of order of expression, every power below -1 would be convergent and every power including -1 and above will not be)

I would be glad if I you could guide me tackle this as there is no clear way for me to go here . This also seems to be a very interesting problem to me as this makes me question both the standard p-series test and also the result from O(sin(n)/n), (n, oo)) which I still feel should be O(1/n, (n, oo))

[skirpichev]

I feel maybe we could go with the sequential leading term protocol that we have been following.

Any meaningful definition of the multivariate Order should a restriction of the quoted above. Is this happens in your case?

Second, in which case your definition could be useful practically? Currently, the univariate order is used only for series (i.e. to get a form Add of terms + Order)... Does you definition will work for multivariate series? (I doubt.)

if p_series_test[p] >= -1

This branch seems to be wrong.

the result from O(sin(n)/n), (n, oo)) which I still feel should be O(1/n, (n, oo))

That's a correct simplification.

[anutosh491]

Second, in which case your definition could be useful practically? Currently, the univariate order is used only for series (i.e. to get a form Add of terms + Order)... Does you definition will work for multivariate series? (I doubt.)

Okay , I guess I need to give more thought about these things . I haven't taken into account the usage and other thing it could affect as of now .I'll give more thought on this and get back to you . Till then I would like to hear your views on the same . Did you stumble upon any approach or idea since when I told you about the issue here ?? If you think there's a somewhat correct direction to proceed here , I would like to build upon that.

if p_series_test[p] >= -1

Could you maybe tell me how though ? I mean p_series test is standard right ! whenever the power is greater than or equal to -1 the series is divergent !

[skirpichev]

whenever the power is greater than or equal to -1 the series is divergent !

Your above sum with the term sin(n)/n is a counterexample.

[anutosh491]

Another issue I have been trying to solve is based on a log rewrite for series expansion at a complex branch cut . Diofant is close but fails due to different reasons , sympy maybe needs more work here .
Consider this example

>>> asinh(x).rewrite(log)
log(x + sqrt(x**2 + 1))
>>> asinh(x).rewrite(log).series(x, I*oo)
nan
>>> asinh(x).rewrite(log).series(x, oo)
-3/(32*x**4) + 1/(4*x**2) + log(2) + log(x) + O(x**(-6), (x, oo))

The results aren't too different as wolfram shows

The case of -oo and -I*oo are also pretty similar , I will be trying to approach this in the coming days and will keep you updated if I find a solution for the same . Sympy gives a huge traceback of errors hence I feel diofant would be closer here !
@skirpichev

[skirpichev]

First result is odd. Second is correct, isn't? -oo case seems incorrect too (by -1 factor?).

[anutosh491]

Yupp the first result needs improvements, diofant and sympy both haven't focused too much on complex infinities (not zoo but something like I*oo which are important branch cuts in quite some functions) though diofant looks closer to me cause sympy introduced a system of cdir for handling oo and -oo separately ( -1 for oo and + 1 for -oo), now it gets weird considering that theses cdir constants cannot be used for something like I*oo so we need to maybe introduce something new or some related solution .

[skirpichev]

#1230
#1231

[anutosh491]

Another issue I had seen in diofant and thought of letting you know was the following . Series for the inverse error functions fails at two prominent points which are +1 and -1

>>> erfinv(x).series(x, 1)
Traceback (most recent call last):
sympy.core.function.PoleError: Cannot expand erfinv(_x + 1) around 0

[skirpichev]

Series for the inverse error functions fails at two prominent points which are +1 and -1

Thanks, PoleError is a wrong error here. See #1233.

I would appreciate a bugreport next time, however.