ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Wed, 12 May 2021 01:17:44 +0200Ranging the Z axishttps://ask.sagemath.org/question/55355/ranging-the-z-axis/I am plotting arrangements of hyperplanes. Like this one:
H3.<x, y, z> = HyperplaneArrangements(QQ)
A = H3([(1, 2, 1), 0],
[(-4, -3, -2), 0],
[(3, 3, -2), 0],
[(1, -4, 3), 0],
[(-2, 2, 1), 0])
A.plot(ranges=[(-1, 1), (-1, 1)], aspect_ratio=(1, 1, 0.25))
![Hyperplane arrangement](/upfiles/16110758982120743.png)
The problem is I need the z axis to range from -1 to 1 rather than from -4 to 4.
Is there a way to change that?LarsTue, 19 Jan 2021 15:36:01 +0100https://ask.sagemath.org/question/55355/A certain polyhedron-related computation: any way to get it quicker?https://ask.sagemath.org/question/57082/a-certain-polyhedron-related-computation-any-way-to-get-it-quicker/(This is a white flag sort of question.) I have a numeric conjecture that I was trying to verify. The conjecture is this:
- take a polyhedron P whose vertices are a poset
- extend this to a (non-reflexive) operation on all faces by saying that F is smaller than G if max vertex of F is
smaller than min vertex of G
- take an algebra R of rational functions on generators x_F for all faces of P, plus extra generator t
- define the following endomorphism f of R:
- for a face F, let Ch(F) be the set of all chains in F, where for a chain C = F1 <... < Fn its total codimension is codim(F,C) = dimF - sum dim Fi
- for a chain C = F1 < ... < Fn, define its monomial m(C) by taking the product of either x_Fi, if dim F_i>0, or of x_Fi/(1-x_Fi), if Fi is a vertex, and then times t^{codim(F,C)}
- define f by sending x_F to the sum of m(C) for all chains C in Ch(C)
- additionally define an endomorhism I by sending x_F to -x_F and t to -t
CONJECTURE: *very often* f*I will be an involution.
The conjecture is proved for simplices, cubes, and all polygons. It's also disproved in some cases. A key example where I need to know the answer is the associahedron, with Tamari order on vertices.
I've asked several questions here and, thanks to tmonteil, wrote a code that checks the conjecture for the 3D associahedron -- and... my code can't run in any reasonable amount of time! Does the problem above sound like something computationally approachable, given that the number of chains in the 3-dimensional face of associahedron is equal to 22663?
How would you write this code?
Attached, my (very bad) code:
---
import itertools
A = Polyhedron(vertices = [[1,2,6,1],[1,6,2,1],[2,1,6,1],[1,2,3,4],[1,4,1,4],[1,6,1,2],[2,1,3,4],[4,1,4,1],[4,3,2,1],[4,3,1,2],[3,1,2,4],[3,2,1,4],[4,1,2,3],[4,2,1,3]])
Fempty = list(A.face_generator())
F = []
for i in Fempty:
if i.dim() != -1:
F.append(i)
R = Frac(PolynomialRing(QQ,x,len(F)+1))
R.inject_variables()
face = {x:f for x,f in list(zip(R.gens(),F))}
x = {f:x for x,f in list(zip(R.gens(),F))}
V = VectorSpace(QQ,4)
W = Cone([[1,-1,0,0],[1,0,-1,0],[1,0,0,-1],[0,1,-1,0],[0,1,0,-1],[0,0,1,-1]])
elms = A.vertices()
rels = []
for i in A.faces(1):
Li = i.vertices()
if vector(list(Li[0]))-vector(list(Li[1])) in W:
rels.append([Li[0],Li[1]])
else:
rels.append([Li[1],Li[0]])
P = Poset((elms, rels), cover_relations = True, facade = True)
def flarger(f1,f2):
L1 = f1.vertices()
L2 = f2.vertices()
for v in L1:
for w in L2:
if (P.lt(w,v) or v == w) == False:
return False
return True
bigelms = A.faces(0)+A.faces(1)+A.faces(2)+A.faces(3)
bigrels = []
for i in bigelms:
for j in bigelms:
if (flarger(i,j) and i != j):
bigrels.append([i,j])
Q = Poset((bigelms, bigrels), cover_relations = False, facade = True)
elll = list(x[i] for i in bigelms)
relll = []
for i in bigrels:
relll.append([x[i[0]],x[i[1]]])
QQ = Poset((elll, relll), cover_relations = False, facade = True)
def inside(f1,f2):
L1 = face[f1].vertices()
L2 = face[f2].vertices()
for v in L1:
if v not in L2:
return False
return True
def subfaces(feis):
result = []
for item in elll:
if inside(item,feis):
result.append(item)
return result
def sfposet(feis):
Res = Poset(QQ.subposet(subfaces(feis)),facade = True)
return Res
def chains(myx):
Huj = sfposet(myx).chains()
result = []
for i in Huj:
if i != []:
result.append(i)
return result
def p(feis):
answer = feis
if face[feis].dim() == 0:
answer = feis/(1-feis)
return answer
def monprod(monlist):
answer = 1
for item in monlist:
answer = answer*p(item)
return answer
def chdim(ch):
dim = 0
for item in ch:
dim = dim+face[item].dim()
return dim
def fmap(feis):
answer = 0
dimf = face[feis].dim()
Flist = chains(feis)
for item in Flist:
power = dimf - chdim(item)
answer = answer + monprod(item)*(x45^power)
return answer
minuslist = []
for feis in F:
minuslist.append(-x[feis])
minuslist.append(-x45)
fmaplist = []
for feis in F:
fmaplist.append(fmap(x[feis]))
fmaplist.append(x45)
f = R.hom(fmaplist)
I = R.hom(minuslist)
f*I*f*IPolydaryaWed, 12 May 2021 01:17:44 +0200https://ask.sagemath.org/question/57082/Subfaces of a face in Polyhedra package?https://ask.sagemath.org/question/57069/subfaces-of-a-face-in-polyhedra-package/I'm using polyhedra package and I need the operation that, for a given face, provides the list of all its subfaces (as faces in the bigger polyhedron). If I'm doing something like face.as_polyhedron().faces(n), the faces stop being recognised as belonging to the bigger polyhedron. What is the correct way to do that?PolydaryaTue, 11 May 2021 15:58:50 +0200https://ask.sagemath.org/question/57069/In Sage is there support for non-convex polyhedra?https://ask.sagemath.org/question/50929/in-sage-is-there-support-for-non-convex-polyhedra/ Which reference page gives help or examples? Can they be triangulated? I want to triangulate them in order to help turn them into stl files.dart2163Wed, 22 Apr 2020 22:46:11 +0200https://ask.sagemath.org/question/50929/Oriented vertices of Polyhedronhttps://ask.sagemath.org/question/56186/oriented-vertices-of-polyhedron/ Hi, I wonder what is the most convenient way to extract the ordered/oriented vertices of a Polyhedron/PolyhedronFace.
Specifically, say that we have a cube:
sage: cb = Polyhedron(vertices=[(0,0,0), (0,0,1), (0,1,0), (0,1,1), (1,0,0), (1,0,1), (1,1,0), (1,1,1)])
which consists of 6 facets:
sage: cb.facets()
(A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices,
A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices,
A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices,
A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices,
A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices,
A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices)
I can easily know the indices of vertices (in `cb.Vrepresentation`) that form each facets, e.g., for the first facet:
sage: cb.facets()[0].ambient_V_indices()
(0, 1, 4, 5)
However, the four vertices are not properly ordered (linked) as I inspect the adjacency matrix:
sage: cb.vertex_adjacency_matrix()
[0 1 1 0 1 0 0 0]
[1 0 0 1 0 1 0 0]
[1 0 0 1 0 0 1 0]
[0 1 1 0 0 0 0 1]
[1 0 0 0 0 1 1 0]
[0 1 0 0 1 0 0 1]
[0 0 1 0 1 0 0 1]
[0 0 0 1 0 1 1 0]
(from the adjacent matrix the ordered indices of the four vertices should be 0->1->5->4, or reversed as 4->5->1->0).
The proper orientation of the vertices is needed for B-rep. Though I can extract the proper ordering with `vertex_adjacency_matrix()`, how to make sure they're consistently CCW orientated (seen from outside) still confuses me.
I'm kinda sure these are well managed in Sage, because the `obj_repr()` method in [plot3d](https://doc.sagemath.org/html/en/reference/plot3d/sage/plot/plot3d/base.html?highlight=obj_repr#sage.plot.plot3d.base.Graphics3d.obj_repr) seems working well providing a consistent facet-vertex representation:
sage: cb_render = cb.render_solid()
sage: cb_render.obj_repr(cb_render.default_render_params())
['g obj_1',
'usemtl texture7',
['v 0 0 0',
'v 0 0 1',
'v 0 1 0',
'v 0 1 1',
'v 1 0 0',
'v 1 0 1',
'v 1 1 0',
'v 1 1 1'],
['f 8 4 2 6',
'f 7 3 4 8',
'f 8 6 5 7',
'f 3 1 2 4',
'f 5 1 3 7',
'f 6 2 1 5'],
['f 6 2 4 8',
'f 8 4 3 7',
'f 7 5 6 8',
'f 4 2 1 3',
'f 7 3 1 5',
'f 5 1 2 6']]
(notice that the vertex index in OBJ representation starts from 1 instead of 0; the line `f 6 2 1 5` corresponds to the ordering `4->5->1->0` that I want)
Is there a native method for Polyhedron or PolyhdronFace that I can directly use for this?zhaiyuTue, 16 Mar 2021 18:20:01 +0100https://ask.sagemath.org/question/56186/Union of touching polyhedrahttps://ask.sagemath.org/question/55769/union-of-touching-polyhedra/Hi, I'm trying to do a union operation over multiple touching polyhedra into one but couldn't find how to achieve this with Sage. For example, how can I merge these two cubes and dissolve the interior face (the touching face):
# two cubes touching
cube_a = Polyhedron(vertices=[(0,0,0), (0,0,1), (0,1,0), (0,1,1), (1,0,0), (1,0,1), (1,1,0), (1,1,1)])
cube_b = Polyhedron(vertices=[(1,1/2,0), (1,1/2,1), (1,3/2,0), (1,3/2,1), (2,1/2,0), (2,1/2,1), (2,3/2,0), (2,3/2,1)])
# visualize
cube_a.plot() + cube_b.plot()
The polyhedra can be assumed all convex, and each touching pair share a polygon interface. I found polyhedron.intersection(other), but no polyhedron.union(other) or so. Any workaround I can do? Thanks in advance.zhaiyuFri, 19 Feb 2021 13:55:30 +0100https://ask.sagemath.org/question/55769/Polyhedra, facets and verticeshttps://ask.sagemath.org/question/55597/polyhedra-facets-and-vertices/I have a 6-dimensional, convex, compact polyhedron in $\mathbb{R}^{12}$, that I am calling `P1`.
The polyhedron `P1` was specified by giving a large number of inequalities.
When I enter `P1.faces(5)` into Sage, I get the following
(A 5-dimensional face of a Polyhedron in RDF^12 defined as the convex hull of 6 vertices,
A 5-dimensional face of a Polyhedron in RDF^12 defined as the convex hull of 7 vertices,
A 5-dimensional face of a Polyhedron in RDF^12 defined as the convex hull of 7 vertices,
A 6-dimensional face of a Polyhedron in RDF^12 defined as the convex hull of 7 vertices,
A 5-dimensional face of a Polyhedron in RDF^12 defined as the convex hull of 6 vertices,
A 6-dimensional face of a Polyhedron in RDF^12 defined as the convex hull of 8 vertices,
...)
How can there be two 6-dimensional faces, each with a different number of vertices,
on the boundary of `P1`? I have checked that `P1` really does have dimension 6.
Also, when I try to find a list of vertices on the boundary of any 5-dimensional face of `P1`,
for example
P1.faces(5)[0].vertices_list()
I get an error message
AttributeError: 'tuple' object has no attribute 'vertices_list'
However, `P1.faces(5)[0].vertices()` works.
My second question is, how do I get a list of vertices on the boundary of an $n$-dimensional face?IngridSat, 06 Feb 2021 10:10:56 +0100https://ask.sagemath.org/question/55597/values of variables from points in the polyhedron defined by a linear programhttps://ask.sagemath.org/question/55520/values-of-variables-from-points-in-the-polyhedron-defined-by-a-linear-program/SageMath provides a [function](https://doc.sagemath.org/html/en/reference/numerical/sage/numerical/mip.html#sage.numerical.mip.MixedIntegerLinearProgram.polyhedron) for constructing the polyhedron defined by a linear program. Given a point in such a polyhedron, how can I tell which coordinate of the point corresponds to which variable of the linear program?
**ADDED:** The documentation suggests that "the polyhedron is built from the variables stored by the LP solver" and that "they usually match the ones created explicitly when defining the LP", except possibly for Gurobi solver. So, the question seems to reduce to finding the order number of each variable in LP (or an accurate bookkeeping of their creations). Things are even more fuzzy for Gurobi, which one of the best solvers out there.Max AlekseyevSun, 31 Jan 2021 19:51:48 +0100https://ask.sagemath.org/question/55520/Find complementary cones from a conehttps://ask.sagemath.org/question/54256/find-complementary-cones-from-a-cone/Suppose I have a cone constructed from a set of points. In 2d for instance:
cone1 = Cone([[-1, 4], [4, -1]])
The rays of this `cone1` are exactly those two vectors I used to define the cone.
Now, I take the cone
cone2 = Cone([[1, 0], [0, 1]])
Again its rays are given by the two vectors used to define it.
Now, `cone1` contains `cone2`, and its intersection with `cone2` is exactly `cone2`.
I would like to know if there is a way to extract the complementary cones from `cone1` "subtracting" `cone2`. Basically I would need:
cone3 = Cone([[-1, 4], [0, 1]])
and
cone4 = Cone([[1, 0], [4, -1]])
Is there a way to extract these cones?
I know that in the 2d example I can plot the cones with .plot() and subtract them visualizing them, of course this is an example. What I would like to know if there is an algorithmic way to extract all the "complementary cones" from a bigger cone that contains (completely) another cone. So I would like to find a way to define `cone3` and `cone4` in higher dimensions, and "algorithmically".amininnoTue, 17 Nov 2020 09:29:43 +0100https://ask.sagemath.org/question/54256/Parallelizing Polyhedral volume computationhttps://ask.sagemath.org/question/36859/parallelizing-polyhedral-volume-computation/ For an application, I need to compute volumes of convex polytopes of dimension of about 300. The polytopes of interest have a restricted form:
$$
P:=\{ (x_1,x_2,\ldots,x_n) \} ~\text{with} \\
0\leq x_i\leq x_{i+1} \leq L,\\
l\leq x_i - x_{i+1} \leq u
$$
where $L,u,l$ are known positive constants.
On my laptop, the polytope volume computation routine takes about 3 hours to compute for $n=9$, and freezes for larger $n$, making it useless for most interesting cases of my application.
> Is there an easy way to parallelize the volume routine?
Perhaps by a straightforward use of the [multiprocessing] module, when running on a cluster ?
A more advanced use case, which also arises in my application, lies in computing integrals over the polytope: $\int_{P} f(\mathbf{x}) d\mathbf{x}$, where $f$ is a symmetric function of its arguments. I would greatly be interested in techniques that extend to this case as well.
dracoTheDragonWed, 08 Mar 2017 07:00:25 +0100https://ask.sagemath.org/question/36859/What is the most efficient way to "look up" a face in the face lattice of a polyhedron?https://ask.sagemath.org/question/34485/what-is-the-most-efficient-way-to-look-up-a-face-in-the-face-lattice-of-a-polyhedron/Say I have a polyhedron `p` with face lattice `L = p.face_lattice()`. I want to define `x` as the element of `L` defined as the convex hull of vertices `<0 1 3>` of `p`. What is the most efficient way to define `x`?
For example, consider
sage: p = polytopes.simplex(3)
sage: for v in p.vertices():
print '\tIndex {}:'.format(v.index()), v
....:
Index 0: A vertex at (0, 0, 0, 1)
Index 1: A vertex at (0, 0, 1, 0)
Index 2: A vertex at (0, 1, 0, 0)
Index 3: A vertex at (1, 0, 0, 0)
We see that `p` has four vertices.
The vertices indexed by `0`, `1`, and `3` are the vertices of a face of `p`. This is confirmed:
sage: L = p.face_lattice()
sage: list(L)
[<>,
<0>,
<1>,
<2>,
<3>,
<0,1>,
<0,2>,
<1,2>,
<0,3>,
<1,3>,
<2,3>,
<0,1,2>,
<0,1,3>,
<0,2,3>,
<1,2,3>,
<0,1,2,3>]
I want to define `x` as the face in `p.face_lattice()` given by these vertices. Of course, I could do this by hand with `x = list(L)[12]`, but I want a way to automate this.done_with_fishTue, 16 Aug 2016 11:06:34 +0200https://ask.sagemath.org/question/34485/Projection along affine hullhttps://ask.sagemath.org/question/35487/projection-along-affine-hull/Let `n_1,...,n_r` be integral points in a polyhedron `P`. The paper I am reading refers to "the projection of `P` along `aff(n_1,...,n_r)`". Here `aff(n_1,...,n_r)` is the affine hull of `n_1,...,n_r`. Can I accomplish this projection in sage?done_with_fishTue, 08 Nov 2016 20:09:57 +0100https://ask.sagemath.org/question/35487/VoronoiDiagram returns empty regionshttps://ask.sagemath.org/question/49749/voronoidiagram-returns-empty-regions/Dear all,
Given the set of 2 points
P = [[-2687.19000000000, -2088.53000000000], [-2686.81000000000, -2084.19000000000]]
I would like to compute the voronoi diagram of $P$. When doing so i saw that the regions of the voronoi are:
{P(-2687.19000000000, -2088.53000000000): The empty polyhedron in RDF^0,
P(-2686.81000000000, -2084.19000000000): A 1-dimensional polyhedron in RDF^2 defined as the convex hull of 1 vertex and 1 ray}
In my code, I do need that those regions are not empty! Why is this happening and how can i fix it?
Please advise and thanks in advance.
P.s. Note that when I make the points a little I get that the voronoi is not empty.
Full Code:
P = [[-2687.19, -2088.53], [-2686.81, -2084.19]]
V = VoronoiDiagram(P)
R = V.regions()muradtukSun, 02 Feb 2020 17:43:31 +0100https://ask.sagemath.org/question/49749/Computing projection of a polyhedronhttps://ask.sagemath.org/question/48155/computing-projection-of-a-polyhedron/How to compute the projection of a given polyhedron onto a given smaller-dimensional plane?
E.g., the projection of Polyhedron(vertices=[(-1,-1), (1,3), (5,-2)]) onto the line $y=x$.Max AlekseyevThu, 03 Oct 2019 02:08:37 +0200https://ask.sagemath.org/question/48155/About algorithm for testing whether a point is in a V-polyhedronhttps://ask.sagemath.org/question/48462/about-algorithm-for-testing-whether-a-point-is-in-a-v-polyhedron/ Hi there,
I have two questions and I wonder if any of you may know the answer.
1 What algorithm does sage use to test whether a given point is in a polyhedron, namely in the function polyhedron.contains()?
2. Is there a difference in the algorithm according to whether the V-polyhedron is a V-polytope or is unbounded?
Thank you, and best regards,
GuillermoguillermoTue, 22 Oct 2019 11:42:20 +0200https://ask.sagemath.org/question/48462/How to find the vertices a a polyhedron define by inequalities ,https://ask.sagemath.org/question/48144/how-to-find-the-vertices-a-a-polyhedron-define-by-inequalities/I would like to know if there is a way to obtain the v-representation of a polyhedron from its definition by inequalities like
$\{x \geq 0, y \geq 0, y + x \leq 10\}$
I find the opposite in the documentation. (If possible in higher dimension, because i can always solve for the solution.CyrilleWed, 02 Oct 2019 19:02:40 +0200https://ask.sagemath.org/question/48144/Memory leak in Polyhedron?https://ask.sagemath.org/question/47243/memory-leak-in-polyhedron/If I run the following code:
import gc
while True:
P = Polyhedron([[1,0,0,0], [0,1,0,0], [0,0,1,0], [0,0,0,1]])
faces = P.faces(3)
del P, faces
print "memory usage: " + str(get_memory_usage()) + ", gc: " + str(gc.collect())
the memory usage keeps increasing.
Do you know what is the cause and how to avoid this problem?GioveMon, 22 Jul 2019 16:15:08 +0200https://ask.sagemath.org/question/47243/Looking for stellation of polyhedrahttps://ask.sagemath.org/question/45291/looking-for-stellation-of-polyhedra/ Hi Sagists,
I'm wondering if anyone knows a way to do stellation with the Polyhedron class?
Mathematica has a Stellate function but they admit it is really an augment operation and not a true stellation
thanks
Pat Browne
PatBSun, 03 Feb 2019 18:45:06 +0100https://ask.sagemath.org/question/45291/Global Dependency in Polyhedron Volumehttps://ask.sagemath.org/question/44695/global-dependency-in-polyhedron-volume/ Hi all,
I am new to sage and want to use it in one of my python scripts to calculate the volume of higher dimensional polyhedrons. When executing my script I get a "ZeroDivisionError" for some of the polyhedrons. However when run in an individual script they work perfectly fine as long as there is no other calculation beforehand.
Here is an example, the file is called DebugSageVolume.py and I execute it with "sage -python DebugSageVolume.py" under Ubuntu 18.04, SageMath version 8.1
#!/usr/lib/python2.7/dist-packages/sage -python
from sage.all import *
# this works fine
matrix1 = matrix(RR, [[0, 1, 0, 0, 0], [1, -1, 0, 0, 0],
[0, 0, 1, 0, 0], [1, 0, -1, 0, 0],
[0, 0, 0, 1, 0], [1, 0, 0, -1, 0],
[0, 0, 0, 0, 1], [1, 0, 0, 0, -1],
[0.5, -0.5, -0.5, -0.5, -0.5], [-0.5, 0.5, 0.5, 0.5, 0.5]])
poly1 = Polyhedron(ieqs=matrix1)
volume1 = poly1.volume(measure="induced")
print(volume1)
#this results in an error when executed after the above. But when commenting the above this works fine
matrix2 = matrix(RR, [[0, 1, 0, 0, 0], [1, -1, 0, 0, 0],
[0, 0, 1, 0, 0], [1, 0, -1, 0, 0],
[0, 0, 0, 1, 0], [1, 0, 0, -1, 0],
[0, 0, 0, 0, 1], [1, 0, 0, 0, -1],
[0.75, -0.5, -0.5, -0.5, -0.5], [-0.75, 0.5, 0.5, 0.5, 0.5]])
poly2 = Polyhedron(ieqs=matrix2)
volume2 = poly2.volume(measure="induced")
print(volume2)
Any help would be greatly appreciated, thank you.
MichaelMichaelRMon, 17 Dec 2018 19:10:30 +0100https://ask.sagemath.org/question/44695/Viewer 3D issueshttps://ask.sagemath.org/question/44620/viewer-3d-issues/I'm trying to draw some 3D polyhedra. A search has brought up
https://ask.sagemath.org/question/41910/viewer-3d-problem/
but I still can't get a live plot. I have constructed my polyhedron `V` from the `Polyhedron` class as an intersection of linear inequalities. The command `V.plot()` works by using the JMol engine, except that when I attempt to make it "live" all the edges and faces disappear and I have collection of vertices only.
I can also plot with tachyon; that works fine but doesn't allow rotating.
When I plot with threejs: `V.plot(viewer='threejs',online=True)` all it plots is the wireframe box surrounding the polyhedron; that is, the xyz "box" in which the polyhedron fits, but not the polyhedron itself.
So - how can I create a 3D polyhedron with a "live" plot that allows rotating? AlasdairMon, 10 Dec 2018 05:49:34 +0100https://ask.sagemath.org/question/44620/Computing polyhedron from MIPhttps://ask.sagemath.org/question/33782/computing-polyhedron-from-mip/ I want to run the following code (I'm using SAGE v6.10 release version 2015-12-18):
P = MixedIntegerLinearProgram()
x = P.new_variable()
A = random_matrix(RR, 3, 2);
P.add_constraint(A*x <= [2.1,1.5,0.4])
P.polyhedron()
However, it outputs the error message:
> AttributeError: type object 'float'
> has no attribute 'fraction_field'
I noticed that if instead of using real variables I use only integer variables, then it works fine. For instance:
P = MixedIntegerLinearProgram()
x = P.new_variable()
A = random_matrix(ZZ, 3, 2);
P.add_constraint(A*x <= [2,1,0])
P.polyhedron()
outputs:
> A 2-dimensional polyhedron in QQ^2
> defined as the convex hull of 1 vertex
> and 2 rays
What do you suggest to overcome this problem?
ThanksmforetsTue, 14 Jun 2016 14:09:09 +0200https://ask.sagemath.org/question/33782/why I get huge number in H-representation in sage?https://ask.sagemath.org/question/43613/why-i-get-huge-number-in-h-representation-in-sage/ I have a set of vertices and I want to find a H-representation of them. I used sage to do that but I got weird number in the inequalities!
here is my code:
vert2 = [[1,2,3,4,2,3,6,0,0,7,6,8,0,1,2,3,4,5,6,7,4,3,2,3,4,1,0,3,0,1],[1,1,2,3,4,2,3,6,0,0,7,6,8,0,1,2,3,4,5,6,7,4,3,2,3,4,1,0,3,0],[0,1,1,2,3,4,2,3,6,0,0,7,6,8,0,1,2,3,4,5,6,7,4,3,2,3,4,1,0,3],[3,0,1,1,2,3,4,2,3,6,0,0,7,6,8,0,1,2,3,4,5,6,7,4,3,2,3,4,1,0],[0,3,0,1,1,2,3,4,2,3,6,0,0,7,6,8,0,1,2,3,4,5,6,7,4,3,2,3,4,1],[1,0,3,0,1,1,2,3,4,2,3,6,0,0,7,6,8,0,1,2,3,4,5,6,7,4,3,2,3,4],[4,1,0,3,0,1,1,2,3,4,2,3,6,0,0,7,6,8,0,1,2,3,4,5,6,7,4,3,2,3],[3,4,1,0,3,0,1,1,2,3,4,2,3,6,0,0,7,6,8,0,1,2,3,4,5,6,7,4,3,2],[2,3,4,1,0,3,0,1,1,2,3,4,2,3,6,0,0,7,6,8,0,1,2,3,4,5,6,7,4,3],[3,2,3,4,1,0,3,0,1,1,2,3,4,2,3,6,0,0,7,6,8,0,1,2,3,4,5,6,7,4],
[4,3,2,3,4,1,0,3,0,1,1,2,3,4,2,3,6,0,0,7,6,8,0,1,2,3,4,5,6,7],[7,4,3,2,3,4,1,0,3,0,1,1,2,3,4,2,3,6,0,0,7,6,8,0,1,2,3,4,5,6],[6,7,4,3,2,3,4,1,0,3,0,1,1,2,3,4,2,3,6,0,0,7,6,8,0,1,2,3,4,5],[5,6,7,4,3,2,3,4,1,0,3,0,1,1,2,3,4,2,3,6,0,0,7,6,8,0,1,2,3,4],[4,5,6,7,4,3,2,3,4,1,0,3,0,1,1,2,3,4,2,3,6,0,0,7,6,8,0,1,2,3],[3,4,5,6,7,4,3,2,3,4,1,0,3,0,1,1,2,3,4,2,3,6,0,0,7,6,8,0,1,2],
[2,3,4,5,6,7,4,3,2,3,4,1,0,3,0,1,1,2,3,4,2,3,6,0,0,7,6,8,0,1],[1,2,3,4,5,6,7,4,3,2,3,4,1,0,3,0,1,1,2,3,4,2,3,6,0,0,7,6,8,0],
[0,1,2,3,4,5,6,7,4,3,2,3,4,1,0,3,0,1,1,2,3,4,2,3,6,0,0,7,6,8],[8,0,1,2,3,4,5,6,7,4,3,2,3,4,1,0,3,0,1,1,2,3,4,2,3,6,0,0,7,6],
[6,8,0,1,2,3,4,5,6,7,4,3,2,3,4,1,0,3,0,1,1,2,3,4,2,3,6,0,0,7],[7,6,8,0,1,2,3,4,5,6,7,4,3,2,3,4,1,0,3,0,1,1,2,3,4,2,3,6,0,0],
[0,7,6,8,0,1,2,3,4,5,6,7,4,3,2,3,4,1,0,3,0,1,1,2,3,4,2,3,6,0],[0,0,7,6,8,0,1,2,3,4,5,6,7,4,3,2,3,4,1,0,3,0,1,1,2,3,4,2,3,6],
[6,0,0,7,6,8,0,1,2,3,4,5,6,7,4,3,2,3,4,1,0,3,0,1,1,2,3,4,2,3],[3,6,0,0,7,6,8,0,1,2,3,4,5,6,7,4,3,2,3,4,1,0,3,0,1,1,2,3,4,2],
[2,3,6,0,0,7,6,8,0,1,2,3,4,5,6,7,4,3,2,3,4,1,0,3,0,1,1,2,3,4],[4,2,3,6,0,0,7,6,8,0,1,2,3,4,5,6,7,4,3,2,3,4,1,0,3,0,1,1,2,3],
[3,4,2,3,6,0,0,7,6,8,0,1,2,3,4,5,6,7,4,3,2,3,4,1,0,3,0,1,1,2],[2,3,4,2,3,6,0,0,7,6,8,0,1,2,3,4,5,6,7,4,3,2,3,4,1,0,3,0,1,1]]
p2=Polyhedron(vertices = vert2)
p2
p2.Hrepresentation()
And here is the part of my answer:
A 29-dimensional polyhedron in ZZ^30 defined as the convex hull of 30 vertices
(An equation (1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1) x - 91 == 0, An inequality (0, -6643755617445331037202609358, -1115451987061522518725310169, -4620824246361587837443723580, -3077838211630051976366550554, -2958690652672660811230373721, -974408181934402275593859480, -3250586099217454711737435326, -1999390120369594790625124513, -1415422103472822894852139772, -2534807215524130595973233650, -2780858391260346768723034201, -4072521014492760746645432088, -1331138469264064061646992630, -4722885096499104848736533209, -3180475210961898948969228316, -3851685143879995858104878858, -3228011163515660410829231425, -3072642994521821646580427696, -1374868705058151838371429430, -4889037447034222629220355689, -2510460065219895800789661380, -2297947035840486945050864530, -3657717103544642567050992273, -3458628241192993528820697768, -295584745844116667823309878, -3868958368255490017832885737, -6000440427999093271578675404, 853080186343081392188416198, -2835167870342318913892483281) x + 260907258623794816872847788563 >= 0, An inequality (0, -3688248056685400306080899479, -853080186343081392188416198, -7496835803788412429391025556, -1968532173404603910913726367, -5473904432704669229632139778, -3930918397973133368554966752, -3811770839015742203418789919, -1827488368277483667782275678, -4103666285560536103925851524, -2852470306712676182813540711, -2268502289815904287040555970, -3387887401867211988161649848, -3633938577603428160911450399, -4925601200835842138833848286, -2184218655607145453835408828, -5575965282842186240924949407, -4033555397304980341157644514, -4704765330223077250293295056, -4081091349858741803017647623, -3925723180864903038768843894, -2227948891401233230559845628, -5742117633377304021408771887, -3363540251562977192978077578, -3151027222183568337239280728, -4510797289887723959239408471, -4311708427536074921009113966, -1148664932187198060011726076, -4722038554598571410021301935, -6853520614342174663767091602) x + 338537555581015223561993662581 >= 0, naghmeh_msFri, 07 Sep 2018 17:55:00 +0200https://ask.sagemath.org/question/43613/Getting convex hull of a set of points and plotting the polygonhttps://ask.sagemath.org/question/33548/getting-convex-hull-of-a-set-of-points-and-plotting-the-polygon/I'm using sagemath cloud.
I'm trying to get the convex hull of a finite set of points, then plotting the polygon. If I just plot the polygon with the vertices directly, I don't get the vertices in the right order to make up a convex polygon (plots edges joining the wrong points). Then I found out about cyclic_sort_vertices_2d which I thought would sort the vertices into the right order, but that's giving me errors. I don't understand the meaning of the errors. Appreciate any help.
here's the relevant code:
hull = Polyhedron(vertices=thePoints);
hull.vertices();
cyclic_sort_vertices_2d(hull.vertices());
And I have this import line at the top of my code for cyclic sort:
from sage.geometry.polyhedron.plot import cyclic_sort_vertices_2d;
Here's my output:
(A vertex at (9464.0, 10816.0), A vertex at (8736.0, 7735.0), A vertex at (7904.0, 8736.0), A vertex at (11088.0, 5607.0), A vertex at (15552.0, 6480.0), A vertex at (8160.0, 9248.0), A vertex at (8568.0, 9792.0), A vertex at (11076.0, 10140.0), A vertex at (12852.0, 5355.0))
And errors:
Error in lines 28-30
Traceback (most recent call last):
File "/projects/sage/sage-6.10/local/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 905, in execute
exec compile(block+'\n', '', 'single') in namespace, locals
File "", line 1, in <module>
File "/projects/sage/sage-6.10/local/lib/python2.7/site-packages/sage/geometry/polyhedron/plot.py", line 245, in cyclic_sort_vertices_2d
adjacency_matrix = Vlist[0].polyhedron().vertex_adjacency_matrix()
File "sage/misc/cachefunc.pyx", line 2215, in sage.misc.cachefunc.CachedMethodCallerNoArgs.__call__ (/projects/sage/sage-6.10/src/build/cythonized/sage/misc/cachefunc.c:14346)
self.cache = f(self._instance)
File "/projects/sage/sage-6.10/local/lib/python2.7/site-packages/sage/geometry/polyhedron/base.py", line 1847, in vertex_adjacency_matrix
return self._vertex_adjacency_matrix()
File "/projects/sage/sage-6.10/local/lib/python2.7/site-packages/sage/geometry/polyhedron/base.py", line 312, in _vertex_adjacency_matrix
face_lattice = self.face_lattice()
File "sage/misc/cachefunc.pyx", line 2215, in sage.misc.cachefunc.CachedMethodCallerNoArgs.__call__ (/projects/sage/sage-6.10/src/build/cythonized/sage/misc/cachefunc.c:14346)
self.cache = f(self._instance)
File "/projects/sage/sage-6.10/local/lib/python2.7/site-packages/sage/geometry/polyhedron/base.py", line 3172, in face_lattice
face_constructor=face_constructor, required_atoms=atoms_vertices)
File "/projects/sage/sage-6.10/local/lib/python2.7/site-packages/sage/geometry/hasse_diagram.py", line 182, in Hasse_diagram_from_incidences
for coatom, atoms in enumerate(coatom_to_atoms)]
KeyError: (frozenset([]), frozenset([0]))ameetnsharmaThu, 26 May 2016 05:45:37 +0200https://ask.sagemath.org/question/33548/find one interior point of a polyhedronhttps://ask.sagemath.org/question/10829/find-one-interior-point-of-a-polyhedron/Hi,
I have a bunch of inequations and I would like to know if there is a solution. What is the simplest way to achieve this in Sage ? I tried using MILP with success... but the point returned is not very fancy (often on the boundary). In an ideal world, I would like to optimize some quadratic function in order for the solution to be the nicest possible.
Here is a sample and simple example with only one inequality (other constraints are equalities) where I reproduced the output of MILP:
Constraints:
2.0 <= x_0 + x_7 + x_8 <= 2.0
2.0 <= x_1 + x_2 <= 2.0
2.0 <= x_3 + x_5 + x_9 <= 2.0
2.0 <= x_4 + x_6 <= 2.0
1.0 <= x_0 + x_2 + x_9 <= 1.0
1.0 <= x_5 + x_6 + x_8 <= 1.0
0.0 <= x_2 + x_6 <= 1.0
Variables:
x_0 is a continuous variable (min=0.0, max=2.0)
x_1 is a continuous variable (min=0.0, max=2.0)
x_2 is a continuous variable (min=0.0, max=2.0)
x_3 is a continuous variable (min=0.0, max=2.0)
x_4 is a continuous variable (min=0.0, max=2.0)
x_5 is a continuous variable (min=0.0, max=2.0)
x_6 is a continuous variable (min=0.0, max=2.0)
x_7 is a continuous variable (min=0.0, max=2.0)
x_8 is a continuous variable (min=0.0, max=2.0)
x_9 is a continuous variable (min=0.0, max=2.0)
There are plenty of *nice* solutions, one is
x0 = x2 = x5 = x6 = x8 = 1/3
x1 = x4 = 5/3
x3 = x7 = 4/3
But with the solver I got:
sage: p.solve()
sage: p.get_values(p[1],p[2],p[3],p[4],p[5],p[6],p[7],p[8],p[9])
[1, 2, 0, 2, 2, 0, 0, 0, 1]
ThanksvdelecroixFri, 13 Dec 2013 08:39:30 +0100https://ask.sagemath.org/question/10829/Speed of enumerating integer points in polytopehttps://ask.sagemath.org/question/24458/speed-of-enumerating-integer-points-in-polytope/ I'm trying to enumerate all integer points in some polytopes, but sage takes way too long.
This is a pretty minimal example of one of these polytopes.
<pre><code>q = MixedIntegerLinearProgram(maximization=False)
x = q.new_variable(integer=True, nonnegative=True)
q.add_constraint( x[0] == 1)
q.add_constraint(88 * x[2] == -306 * x[0] + 1 * x[1])
q.add_constraint( 1 * x[3] == 42636 * x[0] + 42 * x[1])
q.add_constraint(11 * x[4] == 30804 * x[0] + -42 * x[1])
P = q.polyhedron()
P
# A 1-dimensional polyhedron in QQ^5 defined as the convex hull of 2 vertices
P.Vrepresentation()
# (A vertex at (1, 34/7, 5134/7, 73440, 0), A vertex at (1, 0, 306, 55488, 1632))
</code></pre>
Now, I would like to enumerate the integral points in this (1-dimensional) polytope. Looking at the second coordinate of the two boundary points, it is clear there can be at most five such points (in fact, there are exactly five).
However, `P.integral_points()` seems to take forever.
Is there anything I can do to speed this up? I would like to do similar problems that are higher-dimensional, and have a lot more constraints (bounding hyperplanes).
Earlier, I used a MILP library (COIN) to solve the system, giving one integer point, then further constrain the system to exclude this solution, solve again, etc. This is fast enough for this particular example, but becomes unfeasible when the dimension of the polytope and/or the number of integer points increases.ArieFri, 10 Oct 2014 18:53:00 +0200https://ask.sagemath.org/question/24458/cartesian product polyhedrahttps://ask.sagemath.org/question/39952/cartesian-product-polyhedra/Are there any book or paper citations on how the Cartesian product of polyhedra is computed within the polyhedron base class? Consider the following example input:
P1 = Polyhedron(vertices = [[1,1], [2,3], [3,0]], rays = [[3,2], [5,2]])
P2 = Polyhedron(vertices = [[6,0]], rays = [[1,9], [1,3]])
ExampleP = P1.product(P2)
ExampleP.Vrepresentation()
The output is:
A ray in the direction (0, 0, 1, 3),
A ray in the direction (0, 0, 1, 9),
A vertex at (1, 1, 6, 0),
A vertex at (2, 3, 6, 0),
A vertex at (3, 0, 6, 0),
A ray in the direction (3, 2, 0, 0),
A ray in the direction (5, 2, 0, 0)GDRANICTue, 05 Dec 2017 15:03:48 +0100https://ask.sagemath.org/question/39952/force 0/1 in V-representation of a set of inequalitieshttps://ask.sagemath.org/question/36691/force-01-in-v-representation-of-a-set-of-inequalities/Hi,
If I have an array of inequalities (ieqs) that I define as a polyhedron through the command p = Polyhedron(ieqs = ieqs) and then I get the list of vertices of that polyhedron using [p.Vrepresentation( )], however the vertices are fractional so is there any way to force the vertices to be 0-1?
If not, is there a way so that I convert the vertices myself?
Thanks and sorry if my question is so basicmoatiWed, 22 Feb 2017 01:07:12 +0100https://ask.sagemath.org/question/36691/combinatorial equivalence for Polyhedra / isomorphism for latticeshttps://ask.sagemath.org/question/23180/combinatorial-equivalence-for-polyhedra-isomorphism-for-lattices/Is there an easy way to check whether two polyhedra are combinatorial equivalent, i.e. have have isomorphic face lattices.
this does not work:
Poly1=Polyhedron(vertices=[[0,1],[1,0],[1,1]], base_ring=QQ)
Poly2=Polyhedron(vertices=[[2,0],[2,2],[0,2]], base_ring=QQ)
Poly1.face_lattice()==Poly2.face_lattice()
In this case this would work:
str(Poly1.faces(1))==str(Poly2.faces(1))
but what would be a good way to check this in general?
mfThu, 03 Jul 2014 10:20:12 +0200https://ask.sagemath.org/question/23180/Visualizing Polyhedron with specified rangehttps://ask.sagemath.org/question/30561/visualizing-polyhedron-with-specified-range/I am working on visualizing an unbounded polyhedron with plot command. I have a trouble with the way the default visualizer plots the polyhedron. I would like to define range of (x,y,z) domain it visualizes as well as remove big blue arrow that indicates the direction on which the polyhedron extends to because unboundedness is obvious in my application.
The code I have used is following.
`
dd = Polyhedron(ieqs=[[0,1,-sqrt(2),0],[0,1,sqrt(2),0],[0,0,-sqrt(2),1],[0,0,sqrt(2),1]])
dd.plot(ymax=4)
`
yuyuSun, 08 Nov 2015 14:44:04 +0100https://ask.sagemath.org/question/30561/Can the base ring of a polyhedron be restricted?https://ask.sagemath.org/question/34507/can-the-base-ring-of-a-polyhedron-be-restricted/I have a polyhedron `P` with rational vertices
sage: A = matrix([[1,0,0],[1,0,2],[1,1,1],[1,3/2,0]])
sage: A
[ 1 0 0]
[ 1 0 2]
[ 1 1 1]
[ 1 3/2 0]
sage: P = Polyhedron(A)
sage: P
A 2-dimensional polyhedron in QQ^3 defined as the convex hull of 4 vertices
If we scale `P` by a factor of two, then we get a lattice polytope.
sage: (2*P).is_lattice_polytope()
True
However, the base ring of `P` is still `QQ`
sage: (2*P).parent()
Polyhedra in QQ^3
Since `2*P` is a lattice polytope it seems like it should be possible to restrict the base ring of `2*P` to `ZZ`. Is this possible?done_with_fishThu, 18 Aug 2016 03:03:20 +0200https://ask.sagemath.org/question/34507/