Opened 10 years ago

Closed 8 years ago

#13608 closed defect (fixed)

exp() does not work on mpmath mpf numbers

Reported by: Dan Drake Owned by: Burcin Erocal
Priority: major Milestone: sage-6.2
Component: symbolics Keywords: exp mpmath mpf
Cc: Karl-Dieter Crisman, Martin von Gagern, Eviatar Bach Merged in:
Authors: Martin von Gagern Reviewers: Dan Drake, Karl-Dieter Crisman, Ralf Stephan
Report Upstream: N/A Work issues:
Branch: 8a484a2 (Commits, GitHub, GitLab) Commit: 8a484a2e07b56b1d98a0b4513619551089de13a0
Dependencies: Stopgaps:

Status badges

Description (last modified by Dan Drake)

As reported in https://groups.google.com/forum/?fromgroups=#!topic/sage-support/_SGY3ohrmk8 :

sage: import mpmath
sage: exp(mpmath.mpf('-0.0712959029907420240935')) 

fails with "TypeError: 'int' object is not callable". As Karl-Dieter Crisman discovered, this is because exp(a) by default calls a.exp(), and for mpmath floats, a.exp (no parentheses) returns a Python int representing its binary exponent -- which isn't callable.

Our top-level exp() function should properly evaluate mpf numbers!

Attachments (3)

trac13608-workaround.patch (1.8 KB) - added by Dan Drake 10 years ago.
fix trailing whitespace, make patchbot happy
bug13608a.patch (5.3 KB) - added by Martin von Gagern 10 years ago.
Ensure member is callable, use functions from mpmath module
trac_13608_mpmath.patch (7.6 KB) - added by Martin von Gagern 9 years ago.

Download all attachments as: .zip

Change History (33)

comment:1 Changed 10 years ago by Dan Drake

Description: modified (diff)

comment:2 Changed 10 years ago by Karl-Dieter Crisman

Cc: Karl-Dieter Crisman added

I feel like we should have a more solid workaround for mpmath numbers. Are there other things which don't work? Plus, what about those people who DO want this exp attribute?

comment:3 Changed 10 years ago by Dan Drake

I agree, we need better mpmath integration. The "workaround" patch does fix the problem, although it seems like not a very good solution.

I would say anyone who *does* want the mpmath exp attribute can call it directly.

As for other things that don't work:

  • log and trig functions (sin, arccos, as well as hyperbolics) all fail because they can't convert the mpf to the symbolic ring. That's a different problem than we have for exp.
  • gamma works but throws a strange deprecation error.
  • sqrt works.

So the problem here is unique and due to the "exp" collision.

comment:4 Changed 10 years ago by Dan Drake

Status: newneeds_review

Marking as "needs review" in case that attracts attention of someone who has a better solution...

comment:5 Changed 10 years ago by Karl-Dieter Crisman

I'm wondering whether it's even appropriate to convert mpf things to non-mpf things here. Presumably mpmath does support exp and trigs!

sage: import mpmath
sage: a = mpmath.mpf('-0.0712959029907420240935')
sage: mpmath.gamma(a)
mpf('-14.678778643873606')
sage: mpmath.acos(a)
mpf('1.6421527691351614')

I'd hate to have a check in every single function for mpmath objects, very inelegant. But this seems like the "right" way to handle it (i.e., use mpmath versions of functions for mpf objects).

comment:6 Changed 10 years ago by Martin von Gagern

kcrisman pointed out this issue here in response to my question on Ask Sage. I'd like to see this fixed eventually. I agree with ddrake that the exp issue is special, and might therefore need special attention. If this ticket here should focus on exp alone, please point out or create a new ticket for other functions.

The problem with exp being an attribute but not callable are due to source:sage/symbolic/function.pyx@16735#L721. Perhaps this could be fixed by changing hasattr(args[0], self._name) to callable(getattr(args[0], self._name, None)) to ensure that the attribute not only exists but actually is a function. We might cache the result from getattr to avoid calling that again. With this modification, exp should be treated like most builtin functions, e.g. sin.

kcrisman suggested using functions from the mpmath module where appropriate. According to source:sage/symbolic/function.pyx@16735#L348 there is some customized handling of function evaluation for numpy objects. I guess the same approach should be possible for mpmath as well: detect whether one of the arguments (or the first argument) is an mpf object. Check the namespace first to avoid unneccessary imports. If an object is mpf, then either try a _eval_mpmath_ method of the function object, or a function with the same name from the mpmath module. The _eval_mpmath_ function might default to an implementation that tries to convert arguments to RR and results back to mpf.

As an alternative, we might try to support conversion from mpf to SR. Again there is a precedence with numpy, this time found in source:sage/symbolic/ring.pyx@16681#L140. Providing another morphism from mpf to a RR of the current working precision should be easy. One benefit of this approach is that it allows mixing of argument types. On the other hand, a drawback of this approach might be that the result of the function call will remain a Sage type, not an mpf. To change that, an additional clause in the back conversion code from source:sage/symbolic/function.pyx@16735#L733 might be used.

Should I write a patch for any of this?

comment:7 Changed 10 years ago by Martin von Gagern

Cc: Martin von Gagern added

comment:8 in reply to:  6 ; Changed 10 years ago by Dan Drake

Replying to gagern:

Should I write a patch for any of this?

Sure! If you think you have a good solution, just write a patch. A little improvement is (almost always) better than none...

Changed 10 years ago by Dan Drake

Attachment: trac13608-workaround.patch added

fix trailing whitespace, make patchbot happy

Changed 10 years ago by Martin von Gagern

Attachment: bug13608a.patch added

Ensure member is callable, use functions from mpmath module

comment:9 in reply to:  8 ; Changed 10 years ago by Martin von Gagern

OK, this bug13608a.patch does the callable check described in comment:6, and also the special-casing to use mpmath functions where available, and fall back to using implicit conversions otherwise. There still is no automatic coercion of mpf to SR. If that is still desired, I believe it should be in another commit.

As a drive-by fix, this patch does limit the scope of a number of tryexcept in such a way that an exception encountered while searching for a function does not get confused with an exception caused by executing that function. This could lead to more useful error messages in some cases, and perhaps even affect correctness of behaviour in some other cases. I don't have concrete examples, though, as I don't know a sure way to trigger one of these exceptions in a builtin function.

comment:10 in reply to:  9 ; Changed 10 years ago by Karl-Dieter Crisman

Replying to gagern:

OK, this bug13608a.patch does the callable check described in comment:6, and also the special-casing to use mpmath functions where available, and fall back to using implicit conversions otherwise. There still is no automatic coercion of mpf to SR. If that is still desired, I believe it should be in another commit.

That seems very reasonable. Also, perhaps we should first coerce mpmath numbers to RR(prec=foo) as you suggest above.

Might want to add a test that the exp of mpmath works.

sage: mpmath.mpf('0.5').exp
-1

or something. I think more tests in general would be good - maybe for a complex mpc and some interesting function, for some interesting input. Similarly, an example where the except AttributeError: branch is reached would be nice, since in Sage 5.2 I get

sage: import mpmath
sage: getattr(mpmath,'exp')
<built-in method _sage_exp of MPContext object at 0x117f21600>
sage: modulefn = getattr(mpmath,'exp')
sage: modulefn(mpmath.mpf('0.5'))
mpf('1.6487212707001282')

so it seems the first branch is the one that fixes the problem in question.

As a drive-by fix, this patch does limit the scope of a number of tryexcept in such a way that an exception encountered while searching for a function does not get confused with an exception caused by executing that function. This could lead to more useful error messages in some cases, and perhaps even affect correctness of behaviour in some other cases. I don't have concrete examples, though, as I don't know a sure way to trigger one of these exceptions in a builtin function.

That's a great idea; I always forget about the else clause but it is a good idea to do this here, I think. These are some really nice idiomatic fixes.

There seems to be something missing, by the way, immediately above if len(args) == 1: - maybe hg blame can help here. I think it's just

                #     exp(M)

and I think I may have even been involved on that ticket, but I'm not sure how that didn't get in properly.

Version 0, edited 10 years ago by Karl-Dieter Crisman (next)

comment:11 in reply to:  10 Changed 10 years ago by Martin von Gagern

Replying to kcrisman:

That seems very reasonable. Also, perhaps we should first coerce mpmath numbers to RR(prec=foo) as you suggest above (in that later ticket).

Agreed: coercion to SR should go via the RealField of appropriate precision.

Might want to add a test that the exp of mpmath works.

Good idea, will add that later.

or something. I think more tests in general would be good - maybe for a complex mpc and some interesting function, for some interesting input.

OK, will try one or two later on.

Similarly, an example where the except AttributeError: branch is reached would be nice,

Currently, my arcsin example should cater for that, as there is only mpmath.asin, not mpmath.arcsin. But in the long run, there should be some aliasing mechanism that invokes arin in this case.

One problem here is telling the difference: in an ideal world, both versions of a function would yield identical results, so it's hard to tell which function actually got called. And if there is a case where results differ, then that difference is likely a bug, and nothing we'd want to rely on in tests. So how can we see which implementation got called, short of monkey-patching the module to use a different function?

There seems to be something missing, by the way, immediately above if len(args) == 1: - maybe hg blame can help here. I think it's just

                #     exp(M)

and I think I may have even been involved on that ticket, but I'm not sure how that didn't get in properly.

Can't follow you: why do you believe something is missing? My patch was the result of a hg diff based on 17335:d06cf4b2215d.

comment:12 Changed 10 years ago by Karl-Dieter Crisman

Reviewers: Dan Drake, Karl-Dieter Crisman
Status: needs_reviewneeds_work

All sounds good. You could create a custom function that had two different versions, I suppose, and then delete it, in the doctest... the point is just to make sure that all branches are indeed functioning, even if you comment in the explanation.

As for the exp(M), that wasn't something you made missing, but rather something that at some point went missing, so it might as well be added back. If you look at the syntax in that comment you'll see what I mean; it is nonsensical without that extra line.

comment:13 Changed 10 years ago by Martin von Gagern

Sorry, appears I lost track of this one here. As I'm pretty busy just now, is there any chance someone will pick up my patch and polish it ready for inclusion? If not, it might still be some time before I get around to that.

comment:14 Changed 10 years ago by Burcin Erocal

I don't have time to investigate right now, so just a quick comment.

The code path for evaluating symbolic functions is very time sensitive. We should check for speed regressions before accepting anything that modifies it. I'm afraid adding a callable check, or a new Python function call for mpmath function might be too expensive.

comment:15 in reply to:  14 Changed 9 years ago by Martin von Gagern

Replying to burcin:

We should check for speed regressions before accepting anything that modifies it. I'm afraid adding a callable check […] might be too expensive.

I did the following check:

a = [RDF(i*0.01) for i in range(100)]
timeit('[sin(i) for i in a]', repeat=100, number=10000)

and studied the effect of the following code change:

-        if len(args) == 1 and not hold and not dont_call_method_on_arg and \
-                hasattr(args[0], self._name):
-            return getattr(args[0], self._name)()
+        if len(args) == 1 and not hold and not dont_call_method_on_arg:
+            memberfn = getattr(args[0], self._name, None)
+            if callable(memberfn): 
+                return memberfn()

The time per loop went down from 47 µs to 42 µs. So the cost of that callable seems to be far less than the overhead due to the hasattr in the original code. And in any case the change is very small. Avoiding hasattr but doing a memberfn is not None instead of the callable(memberfn) I got those 42 µs again, so the cost of callable gets lost in the noise. For reference: the direct method invocation i.sin() amounts to 29 µs per loop.

or a new Python function call for mpmath function

To judge the cost of the mpmath check, I used the two-argument function atan2 like this:

a = [(RDF(i), RDF(i + 2)) for i in range(100)]
timeit('[atan2(y, x) for y, x in a]', repeat=10, number=1000)

Subject to the whole patch (basically as submitted before, although I'll attach an updated version shortly) the running time increased from 5.7 ms to 5.8 ms per loop. Not much, perhaps not even significant. And this won't affect GiNaC functions that translate to member invocations, which should cover quite a large number of use cases. Even more after @a526457085e3f07440dc39cb77f93ef3e32ae2af which isn't included in my 5.10 code base yet.

Talking of turning function calls to member invocations: in source:src/sage/symbolic/function.pyx@8118b2b39e3a129a5a6186fdc9917940f92b87a6#L387 the Function.__call__ implementation claims to contain work which enables exp(M) for matrices M. But my experiments showed the GinacFunction.__call__ implementation at source:src/sage/symbolic/function.pyx@8118b2b39e3a129a5a6186fdc9917940f92b87a6#L738 is actually responsible for this behavior. After @a526457085e3f07440dc39cb77f93ef3e32ae2af that will be BuiltinFunction.__call__. I'm a bit surprised about this duplicate code, haven't really figured out whether or not this makes sense.

Changed 9 years ago by Martin von Gagern

Attachment: trac_13608_mpmath.patch added

comment:16 Changed 9 years ago by Martin von Gagern

Status: needs_workneeds_review

The file I just attached has improved doctests compared to my previous attempt. I hope this is ready for inclusion now, in the light of my performance checks from comment:15.

comment:17 Changed 9 years ago by Eviatar Bach

Cc: Eviatar Bach added

comment:18 Changed 9 years ago by Martin von Gagern

Patchbot: Apply trac_13608_mpmath.patch​

(I hope this is the right syntax, and I hope that I may issue such directives as an unprivileged contributor.)

comment:19 Changed 9 years ago by Karl-Dieter Crisman

Let's hope it works!

comment:20 Changed 9 years ago by Jeroen Demeyer

Milestone: sage-5.11sage-5.12

comment:21 Changed 9 years ago by Frédéric Chapoton

yes, in principle anybody can tell the bot which patch to apply. But the bot may choose not to obey its orders... :)

apply trac_13608_mpmath.patch

comment:22 Changed 9 years ago by Frédéric Chapoton

Status: needs_reviewneeds_work

this needs to be rebased

comment:23 Changed 9 years ago by Martin von Gagern

I need some assistance here. As far as I understand things, Sage development is currently moving from mercurial to git, right? So the way to grab the latest sources is not via hg (where I never found out how to update to the latest sources), but instead via git?

So I got myself a checkout of ssh://git@trac.sagemath.org:2222/sage.git, as outlined here. I rebased the patch to that, and then wanted to run the tests locally before committing and pushing them. But a bare source tree won't build sage due to lack of a spkg subtree. And using the sage 5.11 tarball as the basis for git-based development won't work either, I assume, since that tarball is still based on mercurial and doesn't provide the git metadata directories, afaics.

How do I obtain a working snapshot of the current sage development tree? Do I manually copy the spkg tree from the 5.11 tarball? Can I check out a repository of spkgs from somewhere, and if so, where and how? Or should I rebase to the 5.11 release instead of a current checkout? Or upgrade my 5.10-based hg dev tree to 5.11 sources somehow?

comment:24 Changed 9 years ago by For batch modifications

Milestone: sage-6.1sage-6.2

comment:25 Changed 9 years ago by Martin von Gagern

Branch: u/gagern/ticket/13608
Created: Oct 16, 2012, 7:03:08 PMOct 16, 2012, 7:03:08 PM
Modified: Jan 30, 2014, 9:20:52 PMJan 30, 2014, 9:20:52 PM

comment:26 Changed 9 years ago by Martin von Gagern

Commit: 6fbe6b98fc6d135c6b7458af171d75481ab69a55
Status: needs_workneeds_review

OK, I finally got my sage git set up properly.


New commits:

fa26d97Trac #14753: Revert ATLAS package
96bf331mpmath integration for builtin functions.
d40615aMerge tag '6.0' into ticket/13608
343bf74Merge branch 'develop' into ticket/13608
661fc7cAdjust doctest to actual output.
6fbe6b9Renamed memberfn to method.

comment:27 Changed 8 years ago by Ralf Stephan

Branch: u/gagern/ticket/13608public/13608
Commit: 6fbe6b98fc6d135c6b7458af171d75481ab69a558a484a2e07b56b1d98a0b4513619551089de13a0
Reviewers: Dan Drake, Karl-Dieter CrismanDan Drake, Karl-Dieter Crisman, Ralf Stephan
Status: needs_reviewpositive_review

Rebased on 6.2.beta6. Long tests pass in symbolic/ and functions/ (and tutorials).


New commits:

8a484a2Merge branch 'u/gagern/ticket/13608' of ssh://trac.sagemath.org/sage into tmp

comment:28 Changed 8 years ago by Volker Braun

author name

comment:29 in reply to:  28 Changed 8 years ago by Martin von Gagern

Authors: Martin von Gagern

comment:30 Changed 8 years ago by Volker Braun

Branch: public/136088a484a2e07b56b1d98a0b4513619551089de13a0
Resolution: fixed
Status: positive_reviewclosed
Note: See TracTickets for help on using tickets.