Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Mrvasympt #7529

Closed
wants to merge 15 commits into from
Closed

Mrvasympt #7529

wants to merge 15 commits into from

Conversation

avichaldayal
Copy link
Contributor

MrvAsympt algorithm to find asymptotic expansion.

@@ -2468,6 +2468,9 @@ def series(self, x=None, x0=0, n=6, dir="+", logx=None):
if (s1 + o).removeO() == s1:
o = S.Zero

# Try asymptotic expansion
if s1.removeO() == 0:
return self.subs(x, 1/x).aseries(x, n, logx, False).subs(x, 1/x)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It looks for me, this line wasn't covered by your test suite. Perhaps, you should add a test when coeff.series(x, S.Infinity, n, logx=logx) call aseries?

Earlier as soon as coefficient had an infinite series expansion we terminated.
Now other terms are also added as long as they don't have an infinite series expansion
@@ -2646,6 +2649,101 @@ def _eval_nseries(self, x, n, logx):
nseries calls it.""" % self.func)
)

def aseries(self, x, n=6, logx=None, bound=0, hir=False):
"""
This algorithm is directly induced from the limit computational algorithm
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No, that's not a good docstring. See this.

Docstring shouldn't be about the algorithm in first place.

@skirpichev
Copy link
Contributor

@Krastanov

Use the ``hir`` parameter to produce hierarchical series. It stops the recursion
at an early level and may provide nicer and more useful results.

If the most rapidly varrying subexpression of a given expression f is f itself,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

typo *varying

@avichaldayal
Copy link
Contributor Author

better name?

In the thesis, it mentions the parameter as growth-bound

@@ -2646,6 +2649,107 @@ def _eval_nseries(self, x, n, logx):
nseries calls it.""" % self.func)
)

def aseries(self, x, n=6, logx=None, bound=0, hir=False):
"""
Returns the asymptotic expansion of self.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

you can add a reference here.

@avichaldayal
Copy link
Contributor Author

hey, we have n and logx parameters too...

Yes, but the user can call the series method directly. Using aseries provides him the option to mention hir and bound parameters

you can add a reference here.

I added the wikipedia reference in the references section. I'll mention that there.

@skirpichev
Copy link
Contributor

Again. Please keep context of the discussion. Use "Add a line note" buttons or reply to email notifications.

I added the wikipedia reference in the references section.

I'm not sure if you do this right. Please take look to other docstrings and please actually build the documentation.

Yes, but the user can call the series method directly.

I mean that you should document all parameters.

@avichaldayal
Copy link
Contributor Author

The following line : O(x**3).subs(x, exp(-x**2)) works in SymPy live and older commmits.

I tried using git bisect and it points to the following commit:
012a919
But this commit doesn't change anything in order.py. Anyways I compared order.py from this commit to the latest version and created a patch.

When I applied the patch, the substitution works.
I'm unable to find which lines in order.py cause the error though.

@asmeurer
Copy link
Member

The pull request #7636 removed denesting of certain powers x**y**z -> x**(y*z). So something implicitly relies on that working.

@skirpichev
Copy link
Contributor

So something implicitly relies on that working.

No. Someone just don't get familar with git bisect. For commit above the one from #7636:

$ git log --oneline
..
0700b30 watch for improper separation of base and exp in powsimp
8fc5951 Merge pull request #7634 from sushant-hiray/coth-fix
..
$ git checkout 8fc5951
$ cat a.py
from sympy import *
x = Symbol('x')
print(O(x**3).subs(x, exp(-x**2)))
$ python a.py
...
ValueError: expecting a Symbol but got nan

@asmeurer
Copy link
Member

I haven't seen any work on this, so I'm going to bump the milestone. Please let me know if it should be put back on.

@skirpichev
Copy link
Contributor

closed in favor of #8542

@skirpichev skirpichev closed this Dec 1, 2014
@asmeurer asmeurer reopened this Jan 19, 2015
@mihir-w
Copy link
Contributor

mihir-w commented Feb 9, 2015

Hi! I have read the algorithm and would like to work upon the XFAIL'ed test.
I am still new to git. How do I pick these commits in my local repo ( cherry-pick shows bad revision ) ? Should I clone this repo ? In that case, where do I send the PR ( or commit ) ?
Thanks.

@certik
Copy link
Member

certik commented Feb 9, 2015

Do something like this:

git checkout -b avichaldayal-mrvasympt master
git pull https://github.com/avichaldayal/sympy.git mrvasympt

Then commit more changes and push into your own fork at github and create a new PR.

Note that you also have to resolve conflicts with master (git merge master).

@skirpichev
Copy link
Contributor

On Mon, Feb 09, 2015 at 03:26:10PM -0800, Ondřej Čertík wrote:

Note that you also have to resolve conflicts with master (git merge
master).

Conflicts are resolved in the linked pr #8542.

(BTW, a little example how you, Ondrej, read pr content first ;))

Unsubscribed.

@mihir-w
Copy link
Contributor

mihir-w commented Feb 14, 2015

There seems to be 2 issues:

First:
Substitution of the type O(x**3).subs(exp(-x**2)) fails. This was attempted in #8568 but is yet to be resolved.

Second:
Even if the above error is somehow circumvented the error still remains. Consider the expected output of the XFAIL'ed test:

>>> (exp(2*a*x + 1/x)/2 + exp(2*b*x + 1/x)/2 - exp(a*x + b*x + 1/x))*exp(-2*x**2) + (exp(a*x + 1/x) - exp(b*x + 1/x))*exp(-x**2) + O(exp(-3*x**2), (x, oo))

Just entering the expression in the console raise the error that the answer is dependent on the sign of a and b. The expression is attempted to be simplified but since it does not know about a and b, it raises error. So even if the first error is solved, problem in subs in hierarchical series will remain.

The hierarchical series expansion should not be attempted to be simplified, rather it should be printed as it is. It is just an intermediate expansion.( Explained at Page 95). However when we substitute the final variable:

if hir
    return s.subs(d,exp(logw))

it automatically tries to simplify and fails.

Any help to how to proceed ? Is there a way to just substitute and not try to evaluate it ? Or should I return as a string ?

@asmeurer
Copy link
Member

Don't return strings.

What kind of simplification are you trying to avoid?

@mihir-w
Copy link
Contributor

mihir-w commented Feb 15, 2015

Consider the XFAIL'ed test:

>>> e = exp(1/x + exp(-x**2)*(exp(a*x) - exp(b*x))) - exp(1/x)
>>> e.aseries( x, n=3, hir=True)

Following is the expected output:

(exp(2*a*x + 1/x)/2 + exp(2*b*x + 1/x)/2 - exp(a*x + b*x + 1/x))*exp(-2*x**2) + (exp(a*x + 1/x) - exp(b*x + 1/x))*exp(-x**2) + O(exp(-3*x**2), (x, oo))

In the function, one step before the final answer is returned, following is the intermediate output:

_d*(exp(x*a + 1/x) - exp(x*b + 1/x)) + _d**2*(exp(2*x*a + 1/x)/2 + exp(2*x*b + 1/x)/2 - exp(x*a + x*b + 1/x)) + O(_d**3)

with d = exp(-x**2)

Clearly the answer achieved here is correct. Only substitution of d is left and we have the answer.
Currently the substitution fails which is a different issue. However even if the substitution works, let's say d = exp(-x). Then on substituting:

Traceback (most recent call last):
  File "test_aseries.py", line 66, in <module>
    test_issue_7872()
  File "test_aseries.py", line 58, in test_issue_7872
    print (e.aseries(x, n=3, hir=True) )
  File "/home/nightwing/sympy_final/sympy/sympy/core/expr.py", line 2734, in aseries
    return self.subs(x, xpos).aseries(xpos, n, bound, hir).subs(xpos, x)
  File "/home/nightwing/sympy_final/sympy/sympy/core/expr.py", line 2779, in aseries
    return s.subs(d,exp(logw))
  File "/home/nightwing/sympy_final/sympy/sympy/core/basic.py", line 892, in subs
    rv = rv._subs(old, new, **kwargs)
  File "/home/nightwing/sympy_final/sympy/sympy/core/basic.py", line 1006, in _subs
    rv = fallback(self, old, new)
  File "/home/nightwing/sympy_final/sympy/sympy/core/basic.py", line 983, in fallback
    rv = self.func(*args)
  File "/home/nightwing/sympy_final/sympy/sympy/core/operations.py", line 41, in __new__
    c_part, nc_part, order_symbols = cls.flatten(args)
  File "/home/nightwing/sympy_final/sympy/sympy/core/add.py", line 240, in flatten
    if o.contains(t):
  File "/home/nightwing/sympy_final/sympy/sympy/series/order.py", line 383, in contains
    obj = self.func(expr, *self.args[1:])
  File "/home/nightwing/sympy_final/sympy/sympy/series/order.py", line 216, in __new__
    expr = expr.as_leading_term(*args)
  File "/home/nightwing/sympy_final/sympy/sympy/core/expr.py", line 2857, in as_leading_term
    obj = self._eval_as_leading_term(x)
  File "/home/nightwing/sympy_final/sympy/sympy/core/mul.py", line 1455, in _eval_as_leading_term
    return self.func(*[t.as_leading_term(x) for t in self.args])
  File "/home/nightwing/sympy_final/sympy/sympy/core/expr.py", line 2857, in as_leading_term
    obj = self._eval_as_leading_term(x)
  File "/home/nightwing/sympy_final/sympy/sympy/core/add.py", line 703, in _eval_as_leading_term
    plain = expr.func(*[s for s, _ in expr.extract_leading_order(x)])
  File "/home/nightwing/sympy_final/sympy/sympy/core/add.py", line 646, in extract_leading_order
    if o.contains(of) and o != of:
  File "/home/nightwing/sympy_final/sympy/sympy/series/order.py", line 372, in contains
    l = ratio.limit(s, point)
  File "/home/nightwing/sympy_final/sympy/sympy/core/expr.py", line 2806, in limit
    return limit(self, x, xlim, dir)
  File "/home/nightwing/sympy_final/sympy/sympy/series/limits.py", line 42, in limit
    return Limit(e, z, z0, dir).doit(deep=False)
  File "/home/nightwing/sympy_final/sympy/sympy/series/limits.py", line 157, in doit
    r = gruntz(e, z, z0, dir)
  File "/home/nightwing/sympy_final/sympy/sympy/series/gruntz.py", line 651, in gruntz
    r = limitinf(e0, z)
  File "/home/nightwing/sympy_final/sympy/sympy/series/gruntz.py", line 424, in limitinf
    c0, e0 = mrv_leadterm(e, x)
  File "/home/nightwing/sympy_final/sympy/sympy/series/gruntz.py", line 505, in mrv_leadterm
    f, logw = rewrite(exps, Omega, x, w)
  File "/home/nightwing/sympy_final/sympy/sympy/series/gruntz.py", line 577, in rewrite
    raise NotImplementedError('Result depends on the sign of %s' % sig)
NotImplementedError: Result depends on the sign of sign(2*a - 2*b)

Basically when the value gets substituted, the expression becomes univariate in x. subs by it's nature tries to simplify the expression. While simplifying it encounters a and b whose signs are unknown and hence it raises NotImplementedError. If 'aandb` are some defined integers, everything works out fine.

Other tests did not have any symbols apart from x in them and hence I guess it went unnoticed.

So what is needed that d gets substituted in the given expression without further simplification. I have tried xreplace and replace but for them the substituted part should be a symbol which is not the case here.

@avichaldayal
Copy link
Contributor Author

When finding series expansion of a function, x is substituted with a positive x dummy variable.
Can we do the same for other variables in the expression too?

In this case a and b can be replaced by positive a and positive b which might fix the issue.

@skirpichev
Copy link
Contributor

In this case a and b can be replaced by positive a and positive b which might fix the issue.

Great idea! Lets assume something arbitrary about parameters to get random output in turn :trollface:

@asmeurer
Copy link
Member

The issue here is that Add shouldn't ever raise NotImplementedError. It should catch the error and return an unevaluated expression. The way O() works, the expression isn't going to be wrong if terms aren't absorbed into it. So I would add a try, except NotImplementedError in Add.flatten where it calls o.contains(t).

I also don't think it should call gruntz (that is, no O() simplification should happen automatically), but that's a different issue.

@skirpichev
Copy link
Contributor

The issue here is that Add shouldn't ever raise NotImplementedError.

Yes.

It should catch the error and return an unevaluated expression.

If you have one broken place (contains method) - lets fix one in another, completely unrelated. Sympy way :trollface:

So I would add a try, except NotImplementedError in Add.flatten where it calls o.contains(t).

No. The correct way is enforce the contract "foo.contains(t) should return either boolean or None".

Something like this:

diff --git a/sympy/series/order.py b/sympy/series/order.py
index 33049a0..6208798 100644
--- a/sympy/series/order.py
+++ b/sympy/series/order.py
@@ -369,7 +369,10 @@ def contains(self, expr):
             ratio = self.expr/expr.expr
             ratio = powsimp(ratio, deep=True, combine='exp')
             for s in common_symbols:
-                l = ratio.limit(s, point)
+                try:
+                    l = ratio.limit(s, point)
+                except NotImplementedError:
+                    l = None
                 if not isinstance(l, C.Limit):
                     l = l != 0
                 else:

@asmeurer
Copy link
Member

Yes, returning None does sound better, although that would strictly speaking change the API.

@skirpichev
Copy link
Contributor

Mentioned above "contract" is a part of the API. Raising a not implemented error - is not, in any case.

# 6.6
e = sin(1/x + exp(-x)) - sin(1/x)
assert e.aseries(x) == (1/(24*x**4) - 1/(2*x**2) + 1 + O(x**(-6), (x, oo)))*exp(-x)
assert e.series(x, oo) == (1/(24*x**4) - 1/(2*x**2) + 1 + O(x**(-6), (x, oo)))*exp(-x)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This fails on master.

# 6.15
e3 = lambda x: exp(exp(exp(x)))
e = e3(x)/e3(x-1/e3(x))
assert e.aseries(x, n=3) == 1 + exp(x + exp(x))*exp(-exp(exp(x))) + ((-exp(x)/2 - S.Half)*exp(x + exp(x)) + \
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

+1

def test_issue_7872():
a, b = symbols('a b', integer=True)
e = exp(1/x + exp(-x**2) * (exp(a*x) - exp(b*x))) - exp(1/x)
assert e.aseries(x, n=3, hir=True) == (exp(2*a*x + 1/x)/2 + exp(2*b*x + 1/x)/2 - \
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

+1

@0sidharth
Copy link
Member

Closing this as all the code apart from example 6.15 is already there in master. I will open a separate PR with that example

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

8 participants