Fourier–Motzkin elimination

[tvn]

hi,  just wondering if the Fourier-Motzkin algorithm for eliminating variable from a system of linear inequalities is implemented somewhere in Sage ? 

[john_perry_usm]

sage: search_src('motzkin')
>>>>
sage:

Looks to me like it isn't, but that may not be the best approach. :-)
I could use this at some point myself, so if knowledgeable people can
confirm that it isn't, then perhaps a ticket should be created.

regards
john perry

[Robert Dodier]

Well, Maxima has the function fourier_elim. I don't know anything
about it. It is undocumented, but you can ask about on the mailing
list.
Sorry I can't be more helpful.

best

Robert Dodier

[kcrisman]

On Jan 25, 12:16 am, Robert Dodier <robert.dod...@gmail.com> wrote:
> On Jan 23, 11:47 pm, tvn <nguyenthanh...@gmail.com> wrote:
>
> > hi,  just wondering if the Fourier-Motzkin algorithm for eliminating
> > variable from a system of linear inequalities is implemented somewhere in
> > Sage ?
>
> Well, Maxima has the function fourier_elim. I don't know anything
> about it. It is undocumented, but you can ask about on the mailing
> list.

I believe that the Maxima fourier_elim might do this. Barton Willis
is the person to ask about whether the 'Motzkin' is part of it, on the
Maxima list (though I think he reads this as well). It's a nice
package.

- kcrisman

[kcrisman]

tvn wrote:

> I am not familiar with Maxima at all so even it has the FM method,
> how likely will it be available in Sage ? There are several

All of Maxima is in Sage. So you could use Maxima directly, through

sage: maxima_console()

or you could send Maxima commands in Sage but have them come back to
Sage (some documentation at http://www.sagemath.org/doc/reference/sage/interfaces/maxima.html).

Or you could create a Sage object that has Maxima methods and use
that. I don't know if this would work for a system of equations.
However, I think that some of this is even already exposed in Sage -
see, for instance, the patch (now included) at
http://trac.sagemath.org/sage_trac/attachment/ticket/7325/trac_7325-final-for-real.patch
and in particular the "documentation" for fourier_elim at
http://maxima.cvs.sourceforge.net/viewvc/maxima/maxima/share/contrib/fourier_elim/rtest_fourier_elim.mac
for more info. It would be great if there was more to wrap that was
robust which would help something.

> implementations of FM method in languages such as C though.

In which case you could presumably write a Cython extension/
interface. I don't know how that all works, but there are lots of
resources to help you in the developer guide for Sage and Cython.

Good luck!

- kcrisman

[kcrisman]

And here's a reply from Barton on this subject. tvn, I hope this
helps! Keep in mind that while Sage wraps the thing called %solve
below (use the keyword to_poly_solve=True in a single-variable
context, I think it might be automatic in the multivariable one), we
do not have the ability (yet) to parse all the things it outputs,
particularly some of the "nounforms" below. Using Maxima directly is
often a good choice for very fine-grained things regarding symbolic
manipulation.

- kcrisman

The to_poly_solver automatically dispatches Fourier (aka Fourier-
Motzkin) elimination on systems
of linear inequations (<, <=, >, >=, =, #); examples:

(%i7) %solve([x+y < 5, x + 8 * y < 42, x - y > 1],[x,y]);
(%o7) %union([y+1<x,x<5-y,y<2])

The Fourier elimination code has a pre-processor that converts some
nonlinear inequations to
equivalent linear inequations:

(%i8) %solve([min(x,y) > abs(x-y)],[x,y]);
(%o8) %union([x=y,0<y],[y/2<x,x<y,0<y],[y<x,x<2*y,0<y])

A few things that I can think of:

(1) The output can be overly complex

(%i9) %solve(abs(1-abs(1-abs(x-5)))<5,x);
(%o9) %union([-2<x,x<4],[4<x,x<5],[5<x,x<6],[6<x,x<12],[x=4],
[x=5],[x=6])

(2) Depending on your needs, the simplex method might be much(!)
faster than Fourier
elimination (but Fourier elimination gives more information).

(3) I would not expect the Maxima Fourier elimination code to work all
that well
with binary64 numbers for coefficients--the code has no pivoting
strategy. (I've
not seen an article about pivoting strategies for Fourier
elimination, by the way.)

(4) When the to_poly_solver is unable to solve, it should return a
%solve nounform.
But when Fourier elimiation code is dispatched, I see this is
broken :(

(%i10) %solve(x^2+x+1 < 0,x);
(%o10) %union([-x^2-x-1>0])

Also, the pre-processor doesn't try all that hard (doesn't introduce
radicals) to
convert to linear inequations.

--Barton

[Marshall Hampton]

Sage includes cddlib, which I believe uses the Fourier-Motzkin
algorithm (keeping track of both vertices and inequalities. The
easiest way to access this functionality is probably the Polyhedron
class, for example the following defines an unbouded region in R^3
through its inequalities (a list such as [2,0,1,1] means 2 + 0*x_0 +
1*x_1 + 1*x_2 >= 0):

sage: p = Polyhedron(ieqs = [[1, 0, 1,0], [1, 1, 0,0], [1, 1, -1,0],
[0,0,0,1]])
sage: p
A 3-dimensional polyhedron in QQ^3 defined as the convex hull of 2
vertices and 3 rays.

You can get the vertices and rays defining this region:
sage: print p.vertices()
sage: print p.rays()
[[-1, -1, 0], [-1, 0, 0]]
[[1, 1, 0], [1, 0, 0], [0, 0, 1]]

and from that you could eliminate the third variable by restricting
the rays and vertices to the dimensions you want; for instance
eliminating the last variable:
sage: p2 = Polyhedron(vertices = [x[0:2] for x in p.vertices()], rays
= [x[0:2] for x in p.rays()])
sage: p2
A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 2
vertices and 2 rays.

-Marshall Hampton

On Jan 24, 12:47 am, tvn <nguyenthanh...@gmail.com> wrote:

[tvn]

thank you both of you --- I'll play around with these commands

[Philipp Schneider]

Am 24.01.11 07:47, schrieb tvn:

> hi, just wondering if the Fourier-Motzkin algorithm for eliminating
> variable from a system of linear inequalities is implemented somewhere in
> Sage ?
>

Hi,

it shouldn't be too hard to implement Fourier-Motzkin elemimination
yourself. Here is the definition of a function with projects the
Polyhedron P: (A*x<=b) along the vector c. The Fouier-Motzkin
elimination is just the special case with c = e_j (the unit vector with
Komponent j = 1).

Greetings,
Philipp

def proj_Poly(P,c):
A,b = P
m = A.nrows(); M = range(0,m)
n = A.ncols()
N = [i for i in M if A[i,:]*c < 0]
Z = [i for i in M if A[i,:]*c == 0]
P = [i for i in M if A[i,:]*c > 0]
p = Z + [(i,j) for i in N for j in P]
r = len(p)
D = Matrix(r,n); d = Matrix(r,1)
for i in range(0,r):
if not isinstance(p[i],tuple):
D[i,:] = A[p[i],:]
d[i] = b[p[i]]
else:
(s,t) = p[i]
D[i,:] = (A[t,:]*c)*A[s,:] - (A[s,:]*c)*A[t,:]
d[i] = (A[t,:]*c)*b[s] - (A[s,:]*c)*b[t]
return (D,d)