Opened 8 years ago

Closed 7 years ago

#13608 closed defect (fixed)

exp() does not work on mpmath mpf numbers

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

Description (last modified by ddrake)

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 ddrake 8 years ago.
fix trailing whitespace, make patchbot happy
bug13608a.patch (5.3 KB) - added by gagern 8 years ago.
Ensure member is callable, use functions from mpmath module
trac_13608_mpmath.patch (7.6 KB) - added by gagern 8 years ago.

Download all attachments as: .zip

Change History (33)

comment:1 Changed 8 years ago by ddrake

  • Description modified (diff)

comment:2 Changed 8 years ago by kcrisman

  • Cc kcrisman 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 8 years ago by ddrake

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 8 years ago by ddrake

  • Status changed from new to needs_review

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

comment:5 Changed 8 years ago by kcrisman

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 follow-up: Changed 8 years ago by 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 8 years ago by gagern

  • Cc gagern added

comment:8 in reply to: ↑ 6 ; follow-up: Changed 8 years ago by ddrake

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 8 years ago by ddrake

fix trailing whitespace, make patchbot happy

Changed 8 years ago by gagern

Ensure member is callable, use functions from mpmath module

comment:9 in reply to: ↑ 8 ; follow-up: Changed 8 years ago by 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 ; follow-up: Changed 8 years ago by kcrisman

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 (in that later ticket).

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.

Last edited 8 years ago by kcrisman (previous) (diff)

comment:11 in reply to: ↑ 10 Changed 8 years ago by 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 8 years ago by kcrisman

  • Reviewers set to Dan Drake, Karl-Dieter Crisman
  • Status changed from needs_review to needs_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 8 years ago by 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 follow-up: Changed 8 years ago by burcin

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 8 years ago by 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 8 years ago by gagern

comment:16 Changed 8 years ago by gagern

  • Status changed from needs_work to needs_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 7 years ago by eviatarbach

  • Cc eviatarbach added

comment:18 Changed 7 years ago by 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 7 years ago by kcrisman

Let's hope it works!

comment:20 Changed 7 years ago by jdemeyer

  • Milestone changed from sage-5.11 to sage-5.12

comment:21 Changed 7 years ago by 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 7 years ago by chapoton

  • Status changed from needs_review to needs_work

this needs to be rebased

comment:23 Changed 7 years ago by 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 7 years ago by vbraun_spam

  • Milestone changed from sage-6.1 to sage-6.2

comment:25 Changed 7 years ago by gagern

  • Branch set to u/gagern/ticket/13608
  • Created changed from 10/16/12 19:03:08 to 10/16/12 19:03:08
  • Modified changed from 01/30/14 21:20:52 to 01/30/14 21:20:52

comment:26 Changed 7 years ago by gagern

  • Commit set to 6fbe6b98fc6d135c6b7458af171d75481ab69a55
  • Status changed from needs_work to needs_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 7 years ago by rws

  • Branch changed from u/gagern/ticket/13608 to public/13608
  • Commit changed from 6fbe6b98fc6d135c6b7458af171d75481ab69a55 to 8a484a2e07b56b1d98a0b4513619551089de13a0
  • Reviewers changed from Dan Drake, Karl-Dieter Crisman to Dan Drake, Karl-Dieter Crisman, Ralf Stephan
  • Status changed from needs_review to positive_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 follow-up: Changed 7 years ago by vbraun

author name

comment:29 in reply to: ↑ 28 Changed 7 years ago by gagern

  • Authors set to Martin von Gagern

comment:30 Changed 7 years ago by vbraun

  • Branch changed from public/13608 to 8a484a2e07b56b1d98a0b4513619551089de13a0
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.