Solving rational inequality should give simplified result

[Robert Pollak]

Hello list!

I have already posted the following question to the Maxima mailing list [1],
but I am not sure whether corresponding improvements should go into Sage or Maxima.

In Sage 5.4, when I enter

solve(abs((x-1)/(x-5)) <= 1/3, x)

I get the following DNF as output

#0: solve_rat_ineq(ineq=abs(x-1)/abs(x-5) <= 1/3)
[[x == -1, -6 != 0, -6 != 0], [x == -1, -6 != 0, -6 != 0, -6 != 0], [x
== -1, -6 != 0, -6 != 0], [x == -1, -6 != 0, -6 != 0, -6 != 0], [x == 2,
-3 != 0, -3 != 0], [x == 2, -3 != 0, -3 != 0, -3 != 0], [x == 2, -3 !=
0, -3 != 0], [x == 2, -3 != 0, -3 != 0, -3 != 0], [x == 1], [1 < x, x
< 2], [-1 < x, x < 1]]

This has two issues:

The minor one is that as a Sage newbie, I was confused by the line containing 'solve_rat_ineq'. I the meantime I learned that this just seems to be a byproduct of a failed try to use this function (because solve_rat_ineq cannot deal with 'abs'). Maybe Sage should give something more explicit like "Failed to use solve_rat_ineq" instead, or not print the line at all?

The main issue is that the result should be

[[-1 <= x, x <= 2]]

, where terms like '-6 != 0' are evaluated and the intervals are merged.

I am a Sage newbie, just having browsed its code the first time, and I see that solve_rat_ineq is called by solve_ineq_univar in sage/symbolic/relation.py, which itself is called by solve_ineq. (Therefore I am CCing the relevant author Robert Marik.)

If we had a function that eliminates the constant inequalities and (for the univariate case) merges the intervals, this function could be applied to the result of solve_ineq before returning. Would this be the correct way to fix this issue?

By the way, I do not really understand which Maxima function generates the current result. According to the Maxima list thread, both 'to_poly_solve' and 'fourier_elim' do not generate the '-6 != 0' terms:

Stavros Macrakis proposed:

load(to_poly_solve)
solve(abs((x-1)/(x-5))<=1/3,x)
=>
union([−1<x,x<1],[1<x,x<2],[x=−1],[x=1],[x=2])

and Volker van Nek gave:

load(fourier_elim);
fourier_elim([abs((x-1)/(x-5))<=1/3],[x]);
=>
[x = - 1] or [x = 2] or [x = 1] or [1 < x, x < 2] or [- 1 < x, x < 1]

Best regards,
Robert


[1] http://thread.gmane.org/gmane.comp.mathematics.maxima.general/42727/

[kcrisman]

On Monday, July 1, 2013 5:37:02 AM UTC-4, Robert Pollak wrote:

Hello list!

I have already posted the following question to the Maxima mailing list [1],
but I am not sure whether corresponding improvements should go into Sage or Maxima.

In Sage 5.4, when I enter

solve(abs((x-1)/(x-5)) <= 1/3, x)

I get the following DNF as output

#0: solve_rat_ineq(ineq=abs(x-1)/abs(x-5) <= 1/3)
[[x == -1, -6 != 0, -6 != 0], [x == -1, -6 != 0, -6 != 0, -6 != 0], [x
== -1, -6 != 0, -6 != 0], [x == -1, -6 != 0, -6 != 0, -6 != 0], [x == 2,
-3 != 0, -3 != 0], [x == 2, -3 != 0, -3 != 0, -3 != 0], [x == 2, -3 !=
0, -3 != 0], [x == 2, -3 != 0, -3 != 0, -3 != 0], [x == 1], [1 < x, x
< 2], [-1 < x, x < 1]]

 

The main issue is that the result should be

[[-1 <= x, x <= 2]]

, where terms like '-6 != 0' are evaluated and the intervals are merged.

As you see in your examples below, to_poly_solve and fourier_elim in Maxima simply don't return their information like that, they return them as a union of open intervals and points.  It would be possible to try to merge those intervals, but likely at a large computational cost (since a generic solution probably won't piece together so nicely).  So I don't know whether it would be worth it to try to do that.

 

If we had a function that eliminates the constant inequalities and (for the univariate case) merges the intervals, this function could be applied to the result of solve_ineq before returning. Would this be the correct way to fix this issue?

Presumably, yes.

 

By the way, I do not really understand which Maxima function generates the current result. According to the Maxima list thread, both 'to_poly_solve' and 'fourier_elim' do not generate the '-6 != 0' terms:

sage: F = abs((2*x-2)/(x-5)) <= 2/3

sage: F0 = [F._maxima_()]

sage: F0

[2*abs(x-1)/abs(x-5)<=2/3]

sage: F0[0]

2*abs(x-1)/abs(x-5)<=2/3

sage: F0[0].parent()

Maxima_lib

sage: F0[0].parent().fourier_elim(F0,[x])

[x=-1,-6#0,-6#0]or[x=-1,-6#0,-6#0,-6#0]or[x=-1,-6#0,-6#0]or[x=-1,-6#0,-6#0,-6#0]or[x=2,-3#0,-3#0]or[x=2,-3#0,-3#0,-3#0]or[x=2,-3#0,-3#0]or[x=2,-3#0,-3#0,-3#0]or[x=1]or[1<x,x<2]or[-1<x,x<1]

which is a Maxima object, and that is what generates it.  But I have trouble getting this to happen in the Maxima console. 

(%i5) g:2*abs(x-1)/abs(x-5)<=2/3;

                               2 abs(x - 1)    2

(%o5)                          ------------ <= -

                                abs(x - 5)     3

(%i6) fourier_elim(g,[x]);

(%o6) [x = - 1] or [x = 2] or [x = 1] or [1 < x, x < 2] or [- 1 < x, x < 1]

(%i7) fourier_elim([g],[x]);

(%o7) [x = - 1] or [x = 2] or [x = 1] or [1 < x, x < 2] or [- 1 < x, x < 1]

I feel like this might have something to do with the fact we're using the binary ECL interface, maybe?  

[Nils Bruin]

On Monday, July 1, 2013 8:33:47 AM UTC-7, kcrisman wrote:

I feel like this might have something to do with the fact we're using the binary ECL interface, maybe?  

 Nope. It may have to do with our interface in general, or with the options we set/packages we load. All of the following give the same result:

 F = abs((2*x-2)/(x-5)) <= 2/3

maxima.eval("load(fourier_elim)")
maxima.fourier_elim(F,[x])               #pexpect interface

maxima_calculus.fourier_elim(F,[x]) #string-based interface to maxima_lib

from sage.interfaces.maxima_lib import *
maxima_eval(EclObject([["$fourier_elim"],sr_to_max([F]),sr_to_max([x])])) #assemble appropriate maxima construct and evaluate
Even

maxima.eval("load(fourier_elim);")
maxima.eval("F: abs((2*x-2)/(x-5)) <= 2/3;")
maxima.eval("fourier_elim(F,[x]);")

produces the same output, so I'm pretty sure the state of maxima under the interface is the cause, not something in the way we communicate with maxima. So next thing is probably to try the usual: see if the maxima console behaves the same if you paste in all of sage's maxima init code. If not, we have a real riddle. If so, you could try to pare the init code down to the line(s) that cause the change in behaviour.

[Nils Bruin]

On Thursday, July 10, 2014 11:46:58 AM UTC-7, Nils Bruin wrote:

So next thing is probably to try the usual: see if the maxima console behaves the same if you paste in all of sage's maxima init code. If not, we have a real riddle. If so, you could try to pare the init code down to the line(s) that cause the change in behaviour.

And we have a winner:

(%i1) display2d : false;

(%o1) false
(%i2) domain : complex;

(%o2) complex
(%i3) load(fourier_elim);

(%o3) "/usr/local/sage/sage-git/local/share/maxima/5.33.0/share/fourier_elim/fourier_elim.lisp"
(%i4) F: abs((2*x-2)/(x-5)) <= 2/3;

(%o4) abs(2*x-2)/abs(x-5) <= 2/3
(%i5) fourier_elim(F,[x]);

(%o5) [x = -1,-6 # 0,-6 # 0] or [x = -1,-6 # 0,-6 # 0,-6 # 0]
                             or [x = -1,-6 # 0,-6 # 0]
                             or [x = -1,-6 # 0,-6 # 0,-6 # 0]
                             or [x = 2,-3 # 0,-3 # 0]
                             or [x = 2,-3 # 0,-3 # 0,-3 # 0]
                             or [x = 2,-3 # 0,-3 # 0]
                             or [x = 2,-3 # 0,-3 # 0,-3 # 0] or [x = 1]
                             or [1 < x,x < 2] or [-1 < x,x < 1]

It should really be no surprise that inequalities don't play nice with "domain: complex"

[Robert Pollak]

Hello Nils, thank you for your analysis!

You wrote:

It should really be no surprise that inequalities don't play nice with "domain: complex"

Hm. "domain : complex" is set both in maxima_lib.py and maxima.py. Which one is responsible here?
It has been in there "forever" (I followed it back to commit 0daf6b in 2009.) and I assume it is necessary for many Maxima calls.

Now I have two questions:

* How can Maxima be called without "domain : complex" as a workaround?
The following does not work, it still gives the "!= 0" terms:

maxima.eval("domain : real;")
solve(abs((x-1)/(x-5)) <= 1/3, x)

* What would happen if Maxima were called with the default "domain : real"?

[Nils Bruin]

On Friday, July 11, 2014 2:01:33 AM UTC-7, Robert Pollak wrote:

The following does not work, it still gives the "!= 0" terms:

maxima.eval("domain : real;")
solve(abs((x-1)/(x-5)) <= 1/3, x)

sage: maxima_calculus("domain: real")
real
sage: solve(abs((x-1)/(x-5)) <= 1/3, x)
#0: solve_rat_ineq(ineq=abs(x-1)/abs(x-5) <= 1/3)
[[x == -1], [x == 2], [x == 1], [1 < x, x < 2], [-1 < x, x < 1]]

 

* What would happen if Maxima were called with the default "domain : real"?

If you mean: "what would happen if sage would initialize maxima_calculus (which is maxima_lib)  with domain: real?" -- a lot of doctests would break. Compare

sage: maxima_calculus("domain: real")
real
sage: maxima_calculus("sqrt(x^2)")
abs(x)
sage: maxima_calculus("domain: complex")
complex
sage: maxima_calculus("sqrt(x^2)")
sqrt(x^2)

[Robert Pollak]

If you mean: "what would happen if sage would initialize maxima_calculus (which is maxima_lib)  with domain: real?" -- a lot of doctests would break.

Yes, that's what I meant, and you are confirming my assumption.

But shouldn't Maxima even give an error message when you ask it to solve an inequality in the complex field? I have asked this in the Maxima tracker: https://sourceforge.net/p/maxima/bugs/2783/
Why does "sage: solve(x^2 <= -1, x)" give []? Ether 'solve' or Maxima seems to have switched to "domain: real" to make sense of the input. Maybe that's what 'solve' should always automatically do when it is called on an inequality?

[parisse]

sage: %giac

  --> Switching to Giac <--

giac: solve(abs((x-1)/(x-5)) <= 1/3, x)
list[((x>=-1) and (x<=2))]

[Robert Pollak]

I wrote

> * What would happen if Maxima were called with the default "domain : real"?

Of course I mean: What would happen if Maxima were *always* called like
that.

I have created the Maxima issue 'fourier_elim and "domain : complex"':
https://sourceforge.net/p/maxima/bugs/2783/

[Robert Pollak]

This looks great, and I could confirm this in Giac/Xcas.

However, there is no Sage instance easily available, where this works:
Both the current sage-6.2.ova and Sage online at http://www.sagenb.org,
https://cloud.sagemath.com and http://sagecell.sagemath.org/ don't have
giac support enabled.

[parisse]

That's unfortunate for sage, since giac has other features that sage could take benefit of, like much faster Groebner basis on Q, rational univariate representation, and other multivariate polynomial operations (*, gcd, factorization).

[Dima Pasechnik]

there apparently is an spkg at
http://www.math.jussieu.fr/~han/xcas/sage/giac-1.1.1.spkg

Why don't you giac guys propose it as an experimental/optional Sage package?

[Dima Pasechnik]

oops, there is in fact http://trac.sagemath.org/ticket/12375
(needing review)

Was there a vote carried out for inclusion of such an spkg?
I can't seem to find one.

Dima

>
>

[Robert Pollak]

On Monday, July 1, 2013 11:37:02 AM UTC+2, I wrote:

solve(abs((x-1)/(x-5)) <= 1/3, x)

[...]

The main issue is that the result should be

[[-1 <= x, x <= 2]]

 
I have found a way to get this! A little term massaging gives me the equivalent:
sage: qepcad((x-1)^2 <= (1/3)^2 * (x-5)^2, vars='(x)')
x + 1 >= 0 /\ x - 2 <= 0

Also, qepcad can be used for cleaning up some results of 'solve':
sage: qepcad(qepcad_formula.or_(x == 1, qepcad_formula.and_(1 < x, x < 2)), vars='(x)')
x - 1 >= 0 /\ x - 2 < 0
but it cannot deal with the '-6 != 0' terms.



[Dima Pasechnik]

perhaps a newer qepcad can deal with more things (Sage's qepcad 1.50
is quite outdated by now). Care to make an update of the corresponding
spkg?

[Robert Pollak]

Am 22.07.2014 10:29, schrieb Dima Pasechnik:
> perhaps a newer qepcad can deal with more things (Sage's qepcad 1.50
> is quite outdated by now). Care to make an update of the corresponding
> spkg?

I have been thinking about it.

As a first step I have built qepcad-B-1.69 and its dependency, the
current saclib-2.2.6 (as static lib), under a clean install of
Linuxmint-17-32bit (a derivative of Ubuntu 14.04 LTS). I could
successfully compute one of the examples from the qepcad website. Then
I have added qepcad to the path.

Now for some reason my sage instance (from sage-6.2-i686-Linux-
Ubuntu_13.04_i686.tar.lzma) still says "RuntimeError: Unable to start
QEPCAD". How can I debug this?

Here is my log:

robert@rpmint32 ~/bin/sage-6.2-i686-Linux $ PATH=$PATH:/home/robert/bin/qepcad-B.1.69/bin
robert@rpmint32 ~/bin/sage-6.2-i686-Linux $ export PATH
robert@rpmint32 ~/bin/sage-6.2-i686-Linux $ which qepcad
/home/robert/bin/qepcad-B.1.69/bin/qepcad
robert@rpmint32 ~/bin/sage-6.2-i686-Linux $ ./sage
┌────────────────────────────────────────────────────────────────────┐
│ Sage Version 6.2, Release Date: 2014-05-06 │
│ Type "notebook()" for the browser-based notebook interface. │
│ Type "help()" for help. │
└────────────────────────────────────────────────────────────────────┘
sage: qepcad(x^2<4,vars='(x)')
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
<ipython-input-1-40446c40694e> in <module>()
----> 1 qepcad(x**Integer(2)<Integer(4),vars='(x)')

/home/robert/bin/sage-6.2-i686-Linux/local/lib/python2.7/site-packages/sage/interfaces/qepcad.pyc in qepcad(formula, assume, interact, solution, vars, **kwargs)
1406 use_witness = True
1407
-> 1408 qe = Qepcad(formula, vars=vars, **kwargs)
1409 if assume is not None:
1410 qe.assume(assume)

/home/robert/bin/sage-6.2-i686-Linux/local/lib/python2.7/site-packages/sage/interfaces/qepcad.pyc in __init__(self, formula, vars, logfile, verbose, memcells, server)
764
765 qex = Qepcad_expect(logfile=logfile)
--> 766 qex._send('[ input from Sage ]')
767 qex._send('(' + ','.join(varlist) + ')')
768 qex._send(str(free_vars))

/home/robert/bin/sage-6.2-i686-Linux/local/lib/python2.7/site-packages/sage/interfaces/expect.pyc in _send(self, cmd)
206 def _send(self, cmd):
207 if self._expect is None:
--> 208 self._start()
209 E = self._expect
210 self.__so_far = ''

/home/robert/bin/sage-6.2-i686-Linux/local/lib/python2.7/site-packages/sage/interfaces/expect.pyc in _start(self, alt_message, block_during_init)
446 self._session_number = BAD_SESSION
447 failed_to_start.append(self.name())
--> 448 raise RuntimeError("Unable to start %s"%self.name())
449 self._expect.timeout = None
450 with gc_disabled():

RuntimeError: Unable to start QEPCAD

[Dima Pasechnik]

On 2014-07-22, Robert Pollak <robert...@jku.at> wrote:
> Am 22.07.2014 10:29, schrieb Dima Pasechnik:
>> perhaps a newer qepcad can deal with more things (Sage's qepcad 1.50
>> is quite outdated by now). Care to make an update of the corresponding
>> spkg?
>
> I have been thinking about it.
>
> As a first step I have built qepcad-B-1.69 and its dependency, the
> current saclib-2.2.6 (as static lib), under a clean install of
> Linuxmint-17-32bit (a derivative of Ubuntu 14.04 LTS). I could
> successfully compute one of the examples from the qepcad website. Then
> I have added qepcad to the path.
>
> Now for some reason my sage instance (from sage-6.2-i686-Linux-
> Ubuntu_13.04_i686.tar.lzma) still says "RuntimeError: Unable to start
> QEPCAD". How can I debug this?

perhaps the qepad executables (or links to them) should be in
SAGE_LOCAL/bin ?
Then, there is a script sage-qepad installed by qepad 1.50 spkg,
which does stuff, and might interfere here.
(it also sets an env. var. QE...)

[Robert Pollak]

I wrote that qepcad cannot deal with terms like '-6 != 0'. This shows as follows:

<start of log>
sage: qepcad(-6 != 0, vars='(x)')
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-1-5a8fdb21f19e> in <module>()
----> 1 qepcad(-Integer(6) != Integer(0), vars='(x)')

/home/robert/bin/sage-6.2-i686-Linux/local/lib/python2.7/site-packages/sage/interfaces/qepcad.pyc in qepcad(formula, assume, interact, solution, vars, **kwargs)
1406 use_witness = True
1407
-> 1408 qe = Qepcad(formula, vars=vars, **kwargs)
1409 if assume is not None:
1410 qe.assume(assume)

/home/robert/bin/sage-6.2-i686-Linux/local/lib/python2.7/site-packages/sage/interfaces/qepcad.pyc in __init__(self, formula, vars, logfile, verbose, memcells, server)

749 # and ensure they match up with the variables in the formula.
750 if frozenset(varlist) != (fvars | frozenset(fqvars)):
--> 751 raise ValueError("specified vars don't match vars in formula")
752 if len(fqvars) and varlist[-len(fqvars):] != fqvars:
753 raise ValueError("specified vars don't match quantified vars")

ValueError: specified vars don't match vars in formula
<end of log>

(I see the same on sagecell.sagemath.org, which probably has qepcad-1.50.)

I have also tried the following version:
sage: qepcad(-6 != 0)
Here qepcad.py comes as far as to to try starting the qepcad binary, which unfortunately fails on my machine.
On sagecell.sagemath.org, this version causes the evaluation to hang forever.


However, standalone qepcad can solve this, as shown in the following log. (Probably qepcad-1.50 could also solve this.)
So the first version certainly could be improved in qepcad.py.


<start of log>
robert@rpmint32 ~/bin $ qepcad
=======================================================
Quantifier Elimination
in
Elementary Algebra and Geometry
by
Partial Cylindrical Algebraic Decomposition

Version B 1.69, 16 Mar 2012

by
Hoon Hong
(hh...@math.ncsu.edu)

With contributions by: Christopher W. Brown, George E.
Collins, Mark J. Encarnacion, Jeremy R. Johnson
Werner Krandick, Richard Liska, Scott McCallum,
Nicolas Robidoux, and Stanly Steinberg
=======================================================
Enter an informal description between '[' and ']':
[]
Enter a variable list:
(x)
Enter the number of free variables:
1
Enter a prenex formula:
[-6 /= 0]
.


=======================================================

Before Normalization >
go

At the end of projection phase >
go

Before Choice >
go

Before Solution >
go

An equivalent quantifier-free formula:

TRUE


===================== The End =======================

-----------------------------------------------------------------------------
0 Garbage collections, 0 Cells and 0 Arrays reclaimed, in 0 milliseconds.
493089 Cells in AVAIL, 500000 Cells in SPACE.

System time: 12 milliseconds.
System time after the initialization: 8 milliseconds.
-----------------------------------------------------------------------------
<end of log>