Opened 6 years ago
Last modified 18 months ago
#11590 new defect
Integrating the sgn() function can produce incorrect results
Reported by: | mjo | Owned by: | burcin |
---|---|---|---|
Priority: | major | Milestone: | |
Component: | calculus | Keywords: | |
Cc: | nbruin | Merged in: | |
Authors: | Reviewers: | ||
Report Upstream: | Fixed upstream, in a later stable release. | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: | #12731 |
Description
Actual result:
sage: integrate(x * sgn(x^2 - 1/4), x, -1, 0) 1/2
Since the argument to sgn() has only one root, -1/2, on (-1, 0), there are only two subintervals on which sgn() can have different values. In particular,
sage: sgn(x^2 - 1/4)(x = -0.75) 1 sage: sgn(x^2 - 1/4)(x = -0.25) -1
Now, the original, actual result should be equivalent to the sum of the following:
sage: integrate(x, x, -1, -0.5) -0.375 sage: integrate(-x, x, -0.5, 0) 0.125
So, something went wrong during the initial integration.
Change History (19)
comment:1 Changed 6 years ago by
- Type changed from PLEASE CHANGE to defect
comment:2 Changed 6 years ago by
- Report Upstream changed from N/A to Reported upstream. Little or no feedback.
comment:3 follow-up: ↓ 13 Changed 6 years ago by
- Report Upstream changed from Reported upstream. Little or no feedback. to Reported upstream. Developers acknowledge bug.
Even with #11483, this isn't working right. It should, though - see Barton's post at the Maxima ticket:
By the way: (%i4) load(abs_integrate)$ Correct antiderivative: (%i5) 'integrate(x*signum(x^2-1/4),x); (%o5) abs(x^2-1/4)/2 Correct definite integral (%i6) 'integrate(x*signum(x^2-1/4),x,-1,0); (%o6) -1/4
so I'm not sure why this is still returning the "wrong" thing. Probably something about the integration code...
comment:4 follow-up: ↓ 8 Changed 6 years ago by
I'm stumped on this one. We get the correct antiderivative:
sage: integrate(x*sgn(x^2-1/4),x) 1/2*abs(x^2 - 1/4)
And ECL gives us the right answer, so it's not an environment setting:
sage: from sage.interfaces.maxima_lib import ecl_eval sage: ecl_eval("#$'integrate(x*signum(x^2-1/4),x,-1,0);$") <ECL: ((RAT SIMP) -1 4)>
But going through maxima_eval
is still trouble:
sage: integrate(x*sgn(x^2-1/4),x,-1,0) 1/2 sage: from sage.interfaces.maxima_lib import maxima_eval sage: a = '($INTEGRATE ((MTIMES SIMP) $X ((%SIGNUM SIMP) ((MPLUS SIMP) ((RAT SIMP) (- 1) 4) ((MEXPT SIMP) $X 2))) ) $X -1 0)' sage: maxima_eval(a) <ECL: ((RAT SIMP) 1 2)>
comment:5 follow-up: ↓ 6 Changed 6 years ago by
- Cc nbruin added
Well, I'm glad it wasn't just me!
Could it be that because we get an answer without abs_integrate
it just returns the 'regular' answer? But I don't recall the integrate code doing that, I think once we turn on abs_integrate
it should just 'be on'...
comment:6 in reply to: ↑ 5 Changed 6 years ago by
Replying to kcrisman:
Well, I'm glad it wasn't just me!
Could it be that because we get an answer without
abs_integrate
it just returns the 'regular' answer? But I don't recall the integrate code doing that, I think once we turn onabs_integrate
it should just 'be on'...
It's possible. The abs_integrate code looks like it defines extra definite integration methods, and then loops through them until 'integrate
is gone from the expression. Within maxima, the extra methods obviously get tried first, because we get the right answer. But I suppose something in the ECL integration could be trying the default definite integration first.
comment:7 Changed 6 years ago by
Something else is weird here. In standalone maxima-5.26, we don't get back the wrong answer. But through sage, without abs_integrate, we get 1/2
:
sage: integrate(x * sgn(x^2 - 1/4), x, -1, 0) 1/2 sage: maxima.console() ;;; Loading #P"/home/mjo/src/sage-5.0.beta1/local/lib/ecl/sb-bsd-sockets.fas" ;;; Loading #P"/home/mjo/src/sage-5.0.beta1/local/lib/ecl/sockets.fas" ;;; Loading #P"/home/mjo/src/sage-5.0.beta1/local/lib/ecl/defsystem.fas" ;;; Loading #P"/home/mjo/src/sage-5.0.beta1/local/lib/ecl/cmp.fas" Maxima 5.26.0 http://maxima.sourceforge.net using Lisp ECL 11.1.1 Distributed under the GNU Public License. See the file COPYING. Dedicated to the memory of William Schelter. The function bug_report() provides bug reporting information. (%i1) display2d:false; (%o1) false (%i2) 'integrate(x*signum(x^2-1/4),x); (%o2) 'integrate(x*signum(x^2-1/4),x) (%i3) 'integrate(x*signum(x^2-1/4),x,-1,0); (%o3) 'integrate(x*signum(x^2-1/4),x,-1,0)
Intriguing!
comment:8 in reply to: ↑ 4 Changed 6 years ago by
Some random observations that may or may not be relevant. First, let's separate the "Reader" (parser) and the "Evaluator".
sage: from sage.libs.ecl import * sage: I1=EclObject("#$integrate(x*signum(x^2-1/4),x,-1,0);$").cadr().cdr().car() sage: I2=EclObject("#$'integrate(x*signum(x^2-1/4),x,-1,0);$").cadr().cdr().car() sage: I1 <ECL: (($INTEGRATE) ((MTIMES) $X ((%SIGNUM) ((MPLUS) ((MEXPT) $X 2) ((MMINUS) ((MQUOTIENT) 1 4))))) $X ((MMINUS) 1) 0)> sage: I2 <ECL: ((%INTEGRATE) ((MTIMES) $X ((%SIGNUM) ((MPLUS) ((MEXPT) $X 2) ((MMINUS) ((MQUOTIENT) 1 4))))) $X ((MMINUS) 1) 0)>
Note the subtle difference between $INTEGRATE
(function) %INTEGRATE
(inert integral form). Sage's sr_integral
produces the $INTEGRATE
expression, so the two alternatives tested are *not* equivalent. I suspect those trigger different code-paths. You could try changing
interfaces/maxima_lib.py:205:max_integrate=EclObject("$INTEGRATE")
to max_integrate=EclObject("%INTEGRATE")
instead and see if that solves more than it messes up.
Slightly less likely to make a difference: If you look at the full expression resulting for EclObject("#$1+2$")
, you'll notice an MEVAL*
, whereas maxima_eval is defined with MEVAL
. I've never worked out which one is the proper one to use. So, you should test all combinations (these skip setting up the proper maxima-error-catching environment):
sage: EclObject(["meval*",["quote", I1]]).eval() sage: EclObject(["meval*",["quote", I2]]).eval() sage: EclObject(["meval",["quote", I1]]).eval() sage: EclObject(["meval",["quote", I2]]).eval()
comment:9 follow-up: ↓ 10 Changed 6 years ago by
I took a much stupider approach, but I did essentially try replacing $INTEGRATE
with %INTEGRATE
. Once I discovered the (#$$)
syntax, I tried replacing the call to maxima_eval
in sr_integral
with a simple string substitution. It causes more problems than it solves, I think.
If I replace meval
with meval*
in a few places, I'm able to break some additional things but not fix any =)
I guess the Right Thing to do would be to go read the maxima source...
comment:10 in reply to: ↑ 9 Changed 6 years ago by
Replying to mjo:
I tried replacing the call to
maxima_eval
insr_integral
with a simple string substitution. It causes more problems than it solves, I think.
Yes, that would make things much worse (not to mention a major step back in other ways). The sr_integral
interface tries to avoid string-based conversion, but maxima_lib in general is still capable of it in a much better way than crafting your own string-mangling based on #$...$
. For instance:
sage.calculus.calculus.maxima(sage.calculus.calculus.dummy_integrate(x * sgn(x^2 - 1/4), x, -1, 0))
should work and compares nicely to the pexpect-based
maxima(sage.calculus.calculus.dummy_integrate(x * sgn(x^2 - 1/4), x, -1, 0))
You should really try the commands I suggested earlier to see where the problem is, though (no changes to the source required!). I don't have ready access to a sage build with the right patches, but you obviously do. Even on a clearly insufficiently patched sage I observe different behaviour:
sage: maxima_eval(I1) <ECL: ((RAT SIMP) 1 2)> sage: maxima_eval(I2) RuntimeError: ECL says: Error executing code in Maxima: first: empty argument.
comment:11 follow-up: ↓ 12 Changed 6 years ago by
Sorry, I should have been more clear: I did try everything you suggested. Replacing $INTEGRATE
with %INTEGRATE
causes the same test failures in maxima_lib.py that I got when I swapped out the whole thing for string substitution. It does fix this particular answer, though.
Here's a session with just the abs_integrate patch applied, maxima-5.26. I left the error at the top intact, for whatever reason it doesn't work at all until after I've integrated something.
sage: from sage.libs.ecl import * sage: I1=EclObject("#$integrate(x*signum(x^2-1/4),x,-1,0);$").cadr().cdr().car() --------------------------------------------------------------------------- TypeError Traceback (most recent call last) /home/mjo/src/sage-5.0.beta1/devel/sage-devel/<ipython console> in <module>() /home/mjo/src/sage-5.0.beta1/local/lib/python2.7/site-packages/sage/libs/ecl.so in sage.libs.ecl.EclObject.cadr (sage/libs/ecl.c:5528)() TypeError: cadr can only be applied to a cons sage: integrate(x*sgn(x^2-1/4),x,-1,0) 1/2 sage: I1=EclObject("#$integrate(x*signum(x^2-1/4),x,-1,0);$").cadr().cdr().car() sage: I2=EclObject("#$'integrate(x*signum(x^2-1/4),x,-1,0);$").cadr().cdr().car() sage: I1 <ECL: (($INTEGRATE) ((MTIMES) $X ((%SIGNUM) ((MPLUS) ((MEXPT) $X 2) ((MMINUS) ((MQUOTIENT) 1 4))))) $X ((MMINUS) 1) 0)> sage: I2 <ECL: ((%INTEGRATE) ((MTIMES) $X ((%SIGNUM) ((MPLUS) ((MEXPT) $X 2) ((MMINUS) ((MQUOTIENT) 1 4))))) $X ((MMINUS) 1) 0)> sage: EclObject(["meval*",["quote", I1]]).eval() <ECL: ((RAT SIMP) 1 2)> sage: EclObject(["meval*",["quote", I2]]).eval() <ECL: ((RAT SIMP) -1 4)> sage: EclObject(["meval",["quote", I1]]).eval() <ECL: ((RAT SIMP) 1 2)> sage: EclObject(["meval",["quote", I2]]).eval() <ECL: ((RAT SIMP) -1 4)>
And here's a session with $INTEGRATE
switched to %INTEGRATE
:
sage: from sage.libs.ecl import * sage: I1=EclObject("#$integrate(x*signum(x^2-1/4),x,-1,0);$").cadr().cdr().car() --------------------------------------------------------------------------- TypeError Traceback (most recent call last) /home/mjo/src/sage-5.0.beta1/devel/sage-devel/<ipython console> in <module>() /home/mjo/src/sage-5.0.beta1/local/lib/python2.7/site-packages/sage/libs/ecl.so in sage.libs.ecl.EclObject.cadr (sage/libs/ecl.c:5528)() TypeError: cadr can only be applied to a cons sage: integrate(x*sgn(x^2-1/4),x,-1,0) -1/4 sage: I1=EclObject("#$integrate(x*signum(x^2-1/4),x,-1,0);$").cadr().cdr().car() sage: I2=EclObject("#$'integrate(x*signum(x^2-1/4),x,-1,0);$").cadr().cdr().car() sage: I1 <ECL: (($INTEGRATE) ((MTIMES) $X ((%SIGNUM) ((MPLUS) ((MEXPT) $X 2) ((MMINUS) ((MQUOTIENT) 1 4))))) $X ((MMINUS) 1) 0)> sage: I2 <ECL: ((%INTEGRATE) ((MTIMES) $X ((%SIGNUM) ((MPLUS) ((MEXPT) $X 2) ((MMINUS) ((MQUOTIENT) 1 4))))) $X ((MMINUS) 1) 0)> sage: EclObject(["meval*",["quote", I1]]).eval() <ECL: ((RAT SIMP) 1 2)> sage: EclObject(["meval*",["quote", I2]]).eval() <ECL: ((RAT SIMP) -1 4)> sage: EclObject(["meval",["quote", I1]]).eval() <ECL: ((RAT SIMP) 1 2)> sage: EclObject(["meval",["quote", I2]]).eval() <ECL: ((RAT SIMP) -1 4)>
The same thing, except the original integration actually works. Problem is, other tests start to fail with %INTEGRATE
. Here's an example:
$ sage -t sage/interfaces/maxima_lib.py sage -t "devel/sage-devel/sage/interfaces/maxima_lib.py" ********************************************************************** File "/home/mjo/src/sage-5.0.beta1/devel/sage-devel/sage/interfaces/maxima_lib.py", line 661: sage: integral(x^n,x) Expected: Traceback (most recent call last): ... ValueError: Computation failed since Maxima requested additional constraints; using the 'assume' command before integral evaluation *may* help (example of legal syntax is 'assume(n+1>0)', see `assume?` for more details) Is n+1 zero or nonzero? Got: x^(n + 1)/(n + 1) ********************************************************************** File "/home/mjo/src/sage-5.0.beta1/devel/sage-devel/sage/interfaces/maxima_lib.py", line 685: sage: integrate(sgn(x) - sgn(1-x), x) Expected: abs(x - 1) + abs(x) Got: (x - 1)*sgn(x - 1) + x*sgn(x) ********************************************************************** File "/home/mjo/src/sage-5.0.beta1/devel/sage-devel/sage/interfaces/maxima_lib.py", line 700: sage: integrate(cos(x + abs(x)), x) Expected: 1/4*(sgn(x) + 1)*sin(2*x) - 1/2*x*sgn(x) + 1/2*x Got: -1/4*(2*x - sin(2*x))*sgn(x) + 1/2*x + 1/4*sin(2*x) ********************************************************************** File "/home/mjo/src/sage-5.0.beta1/devel/sage-devel/sage/interfaces/maxima_lib.py", line 711: sage: integral(abs(cos(x))*sin(x),(x,pi/2,pi)) Expected: 1/2 Got: -1/2
comment:12 in reply to: ↑ 11 Changed 6 years ago by
Replying to mjo:
I left the error at the top intact, for whatever reason it doesn't work at all until after I've integrated something.
sage: from sage.libs.ecl import *
This doesn't load maxima inside ecl yet, but once you evaluate an integral, maxima is loaded and the #$...$
macro works. Sorry about that.
comment:13 in reply to: ↑ 3 ; follow-up: ↓ 14 Changed 6 years ago by
Replying to kcrisman:
Even with #11483, this isn't working right. It should, though - see Barton's post at the Maxima ticket:
By the way: (%i4) load(abs_integrate)$ Correct antiderivative: (%i5) 'integrate(x*signum(x^2-1/4),x); (%o5) abs(x^2-1/4)/2 Correct definite integral (%i6) 'integrate(x*signum(x^2-1/4),x,-1,0); (%o6) -1/4so I'm not sure why this is still returning the "wrong" thing.
Did you try it with integrate
rather than 'integrate
? Given what we've seen elsewhere in this ticket, I suspect it still gives the wrong answer. In sage, we are interfacing with integrate
. If that is still broken, then the bug is not fixed as far as sage is concerned.
If maxima's position is that we should use a different integration routine (i.e., 'integrate
) then we need a heuristic on when to use what routine ... isn't it maxima's job to figure this out?
comment:14 in reply to: ↑ 13 Changed 6 years ago by
Did you try it with
integrate
rather than'integrate
?
(%i5) display2d:false; (%o5) false (%i6) integrate(x*signum(x^2-1/4),x,-1,0); (%o6) 1/2 (%i7) 'integrate(x*signum(x^2-1/4),x,-1,0); (%o7) -1/4
Note that this was just Barton Willis (Maxima dev) example.
Given what we've seen elsewhere in this ticket, I suspect it still gives the wrong answer. In sage, we are interfacing with
integrate
. If that is still broken, then the bug is not fixed as far as sage is concerned.
I don't think they claimed it was fixed.
If maxima's position is that we should use a different integration routine (i.e.,
'integrate
) then we need a heuristic on when to use what routine ... isn't it maxima's job to figure this out?
Yes, I am a little stumped on this. I thought that the apostrophe just meant "don't evaluate nounform", but apparently I was mistaken. I just reported this example on the Maxima ticket, but I wouldn't hold my breath waiting for it, since there is a way to get it to work correctly there.
comment:15 Changed 6 years ago by
- Stopgaps set to #12731
comment:17 Changed 3 years ago by
- Stopgaps set to #12731
comment:18 Changed 18 months ago by
FYI fixed in Maxima by commit 5a300aa, which I have cherry-picked to branch-5_38.
Related bug reports: https://sourceforge.net/p/maxima/bugs/2242 https://sourceforge.net/p/maxima/bugs/3123
comment:19 Changed 18 months ago by
- Report Upstream changed from Reported upstream. Developers acknowledge bug. to Fixed upstream, in a later stable release.
This appears to be in Maxima, and is reported at their bug tracker.