Integrating the sgn() function can produce incorrect results

[mjo]

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.

[kcrisman]

This appears to be in Maxima, and is reported at  their bug tracker.

(%i11) integrate(x*signum(x^2-1/4),x,-1,0);
                                       1
(%o11)                                 -
                                       2
(%i15) integrate(x*signum(x^2-1/4),x);
                                   ! 2   1!
                                   !x  - -!
                                   !     4!
(%o15)                             --------
                                      2

[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/4

so I'm not sure why this is still returning the "wrong" thing. Probably something about the integration code...

[mjo]

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)>

[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 on abs_integrate it should just 'be on'...

[mjo]

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 on abs_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.

[mjo]

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!

[nbruin]

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()

[mjo]

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...

[nbruin]

Replying to mjo:

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.

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.

[mjo]

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

[nbruin]

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.

[nbruin]

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/4

so 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?

[kcrisman]

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.