Ask Your Question

jipilab's profile - activity

2019-02-09 02:53:32 -0600 received badge  Necromancer (source)
2019-02-08 08:14:26 -0600 answered a question No Output Given : []

Since your system is linear, I would use matrices to model this. Perhaps I made a typing mistake while translating from what you wrote, but then I get the following. It seems that your system does not have a solution, the matrix is not full-rank over the polynomials ring over the parameters.

sage: PR1 = PolynomialRing(QQ,var('B,V,W,E,Q,H,S,P'))
sage: M = matrix(PR1,[[0,0,0,S+P,B*P,-B*P],[0,0,0,-S,(W*P - S - P),(-W*P - S)],[0,0,0,0,S+P,S],[0,0,S+P,0,0,0],[0,0,-S,0,0,0],[S+P,S,0,0,0,0]])
sage: M
[          0           0           0       S + P         B*P        -B*P]
[          0           0           0          -S W*P - S - P    -W*P - S]
[          0           0           0           0       S + P           S]
[          0           0       S + P           0           0           0]
[          0           0          -S           0           0           0]
[      S + P           S           0           0           0           0]
sage: v = vector(PR1,[-S,2*V*S + V*P,-S,- E*P - S,2*Q*S + Q*P - H*P - S,-S])
sage: M.solve_right(v)
ValueError                                Traceback (most recent call last)
<ipython-input-30-52655db42454> in <module>()
----> 1 M.solve_right(v)

/home/jplabbe/sage/local/lib/python2.7/site-packages/sage/matrix/matrix2.pyx in sage.matrix.matrix2.Matrix.solve_right (build/cythonized/sage/matrix/matrix2.c:7907)()
    447         if self.rank() != self.nrows():
--> 448             X = self._solve_right_general(C, check=check)
    449         else:
    450             X = self._solve_right_nonsingular_square(C, check_rank=False)

/home/jplabbe/sage/local/lib/python2.7/site-packages/sage/matrix/matrix2.pyx in sage.matrix.matrix2.Matrix._solve_right_general (build/cythonized/sage/matrix/matrix2.c:9035)()
    562             # Have to check that we actually solved the equation.
    563             if self*X != B:
--> 564                 raise ValueError("matrix equation has no solutions")
    565         return X

ValueError: matrix equation has no solutions
sage: M.rank()
2019-02-05 15:09:47 -0600 commented answer cartesian product polyhedra

For a reference on the product of polyhedra, see p.390 in the Handbook of Discrete and Computational Geometry, Chapter 15:

Which is the reference which Sage tries to follow for nomenclature.

2019-02-05 15:03:52 -0600 answered a question Looking for stellation of polyhedra


Could give a precise definition of what you mean by stellation?

There is a quick reference page about things you can do with a polyhedron:

Perhaps what you mean is related to the following methods. Say P is a polyhedron object then

  • P.barycentric_subdivision(): constructs a polyhedron out of P whose facets correspond to complete flags of faces of P.
  • P.stack(face): puts a new vertex a little bit outside of the face face.

In general, Sage tries to follow the nomenclature presented in the Handbook of Discrete and Computational Geometry (Chapter 15):

Is the operation described in there?

Otherwise, please let us know what is the operation and we will add it to the desired methods.

2018-12-24 12:20:40 -0600 received badge  Enlightened (source)
2018-12-24 12:20:38 -0600 received badge  Good Answer (source)
2018-11-11 11:00:10 -0600 received badge  Necromancer (source)
2018-11-11 07:06:37 -0600 answered a question I.variety() : missing solution values

According to the documentation of

sage: I.variety?

we read:

Note that with "ring=RR" or "CC", computation is done numerically and potentially inaccurately; in particular, the number of points in the real variety may be miscomputed. With "ring=AA" or "QQbar", computation is done exactly (which may be much slower, of course).

Being miscomputed might also include "missing solutions"... This triangular decomposition also highly depends on the ideal and its groebner basis computation.

Currently, it seems indeed that getting a 0-dimensional variety out of an ideal is rather hazardous.

I am using Bertini to do such computations, eventually is could be interfaced in Sage.

Otherwise, maybe the interface to PHC could do the thing:

2018-07-05 06:47:46 -0600 received badge  Commentator
2018-07-05 06:47:46 -0600 commented answer Can I intersect the boundaries of two polyhedra and display it?

You can also do the intersection operation by using "&":

sage: P1 & P2
2018-07-05 06:44:33 -0600 answered a question Getting convex hull of a set of points and plotting the polygon

Since I see that you have decimal entries, this suggests that you have used "RDF" as a base ring for your polygon.

You may check if it is the case by calling P.base_ring()

There is a known bug in RDF that pops up when plotting:

Otherwise, if you just want to visualize your polygon your should be able to visualize it with the command ".plot()".

More details will be explained in the extended tutorial which will appear in version 8.3:

Under "Geometry", "A longer introduction... " (once version 8.3 is out):

2018-05-02 17:34:21 -0600 received badge  Nice Answer (source)
2018-05-02 17:34:03 -0600 received badge  Nice Answer (source)
2018-05-02 17:33:31 -0600 received badge  Necromancer (source)
2018-05-02 09:19:17 -0600 answered a question find one interior point of a polyhedron

If all you want is some point that would be in the interior of the solution set, you can get one by doing:

sage: q = p.polyhedron()
sage: q.representative_point()

which returns the barycenter of the vertices of the polyhedron and adds the potential generating rays. So it would be in the relative interior of the polyhedron.

2018-05-02 09:09:50 -0600 received badge  Necromancer (source)
2018-05-02 09:08:01 -0600 commented answer Why is this polar polytope incorrect?

Indeed, polarity of polyhedron requires compactness (to get a good duality property, for example that taking twice the dual gives back the original object) and that the origin is in the interior of the (then compact) polyhedron.

Hence, in order to allow polarity of polytope that do not contain the origin, it is automatically shifted so that its barycenter is the origin. This allows for a somehow "canonical" polar object. Taking another point in the interior as the origin would also work, but not implemented.

2018-05-02 08:56:23 -0600 answered a question Speed of enumerating integer points in polytope

In Sage 8.1, I can reproduce your code as follows (notice that I added "ppl" as a solver since the default behavior changed between your post and this answer).

sage: q = MixedIntegerLinearProgram(maximization=False, solver="ppl")
....: 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])

In order to speed things up, I use normaliz.

(You can install it by typing the commands in a terminal:

sage -i normaliz
sage -i pynormaliz


and then use the backend normaliz for the polyhedron:

....: P = q.polyhedron(backend="normaliz")
....: P.integral_points()
((1, 0, 306, 55488, 1632),
 (1, 1, 394, 59184, 1296),
 (1, 2, 482, 62880, 960),
 (1, 3, 570, 66576, 624),
 (1, 4, 658, 70272, 288))

which should appear instantaneously.

2018-05-02 08:35:48 -0600 commented question random polytopes

Note that generating random d-polytopes is a sensitive thing to do: you end up with probability 1-epsilon with a simplicial polytope, which may or may not be what you want to have.

2018-05-02 08:32:31 -0600 received badge  Necromancer (source)
2018-05-02 08:24:44 -0600 answered a question Can sage compute the h* vector of a polytope?

You can use the normaliz backend (requires Normaliz 3.5.4) and its python interface pynormaliz (requires PyNormaliz 1.16).

You can install them by typing:

sage -i normaliz
sage -i pynormaliz

in a terminal once this ticket has been merged. Then you can type in sage:

sage: C = polytopes.hypercube(3, backend="normaliz")
sage: C.hilbert_series().numerator().coefficients()
[1, 3, 6, 7, 6, 3, 1]

Note that this requires the latest features of this ticket.

2018-05-02 08:09:19 -0600 commented question How to solve a large list of simultaneous equations

Which question is it a duplicate of?

2018-01-10 02:58:18 -0600 commented question Obtaining a certain list of posets via SAGE

Since what you are looking for are exactly distributive lattices, perhaps it would be good to change the title of the question to indicate that this is what is looked for. Something like "Obtaining distributive lattices from a list of posets via SAGE"

2017-10-19 01:10:05 -0600 received badge  Necromancer (source)
2017-10-11 10:23:50 -0600 received badge  Nice Answer (source)
2017-10-11 02:06:00 -0600 answered a question Hi, there: I am computing the eigenvalues and eigenvectors of the following complex matrix. It looks that Sage does not give the right answer?

The ring of complex floats CDF is not an exact ring:

sage: CDF.is_exact()

So while doing computation in floating point entries (i.e. in RDF or CDF) some errors may occur (see You can change the base ring to the algebraic field QQbar:

sage: B = matrix(QQbar, [ [-1, 0, 1, 0, 0, 0, 1, I, -1 ], [0, 0, 0, -I, 1, I, 1, I, 0], [-1, -I, 1, -I, 1, 0, 1, 0, I], [0, 0, 0, 0, 1, I, 1, I, -1], [0, -I, 1, -I, 1, I, 1, I, 0], [-1, -I, 1, -I, 1, 0, 0, 0, 0], [-1, 0, 1, 0, 1, I, 1, I, -1], [0, -I, 1, -I, 1, I, 0, 0, 0], [-1, -I, 1, 0, 0, 0, 1, 0, -1]])
sage: B.eigenvalues()
[3.917412364954442? + 0.05429188890228639?*I, 0.6516886324409428? + 0.2353750083092179?*I, 0.3580569417083658? - 0.08398467888210074?*I, -0.1167390153623098? + 0.6697900534213050?*I, -0.1742956154446435? - 0.5989951966469979?*I, -0.7276189211835305? + 0.3726506970413855?*I, -0.8274431770819343? + 1.205961277873422?*I, -0.9961206931117568? - 1.191900422440107?*I, -1.084940516919576? - 0.6631886275784114?*I]

An alternative is to use an algebraic extension of QQusing NumberField:

sage: K.<a> = NumberField(x^2 + 1)
sage: K.is_exact()
sage: B = matrix(K, [ [-1, 0, 1, 0, 0, 0, 1, a, -1 ], [0, 0, 0, -a, 1, a, 1, a, 0], [-1, -a, 1, -a, 1, 0, 1, 0, a], [0, 0, 0, 0, 1, a, 1, a, -1], [0, -a, 1, -a, 1, a, 1, a, 0], [-1,-a, 1, -a, 1, 0, 0, 0, 0], [-1, 0, 1, 0, 1, a, 1, a, -1], [0, -a, 1, -a, 1, a, 0, 0, 0], [-1, -a, 1, 0, 0, 0, 1, 0, -1]])
sage: B.eigenvalues()
[-1.084940516919576? - 0.6631886275784114?*I, -0.9961206931117568? - 1.191900422440107?*I, -0.8274431770819343? + 1.205961277873422?*I, -0.7276189211835305? + 0.3726506970413855?*I, -0.1742956154446435? - 0.5989951966469979?*I, -0.1167390153623098? + 0.6697900534213050?*I, 0.3580569417083658? - 0.08398467888210074?*I, 0.6516886324409428? + 0.2353750083092179?*I, 3.917412364954442? + 0.05429188890228639?*I]
2017-10-11 01:47:59 -0600 commented question Hi, there: I am computing the eigenvalues and eigenvectors of the following complex matrix. It looks that Sage does not give the right answer?

Could you add the values that you obtain and perhaps the version of Sage that you are using? What are the expected eigenvalues? In Sage version 8.1beta7 I get

sage: B.eigenvalues() [3.917412364954441 + 0.054291888902286756I, -0.9961206931117547 - 1.1919004224401037I, -0.8274431770819329 + 1.205961277873422I, -1.0849405169195765 - 0.6631886275784113I, -0.7276189211835324 + 0.3726506970413847I, -0.1167390153623094 + 0.6697900534213056I, 0.6516886324409439 + 0.23537500830921726I, 0.3580569417083653 - 0.08398467888210029I, -0.174295615444644 - 0.598995196646998*I]

Is that correct?

2017-10-10 15:12:28 -0600 commented answer Obtaining lattices quickly in SAGE

This is true. I appended the answer.

2017-10-10 13:51:09 -0600 answered a question Projection along affine hull

Here is an example. It requires to load the Projection function:

sage: from sage.geometry.polyhedron.plot import Projection

Let's say we have the following polytope:

sage: P = Polyhedron(vertices=[[-2,0,3],[2,0,3],[0,2,1],[0,-2,1]])

and we are interested in the affine space spanned by the points

sage: affine_basis = [vector([1,1,2]),vector([1,-1,2]),vector([-1,-1,2])]

we first get a linear subspace and create the projection matrix:

sage: linear_subspace = [ap - affine_points[0] for ap in affine_basis[1:]]
sage: VS = VectorSpace(QQ,3)
sage: proj_matrix = matrix([sum([v.inner_product(lv)*lv/(lv.inner_product(lv)) for lv in linear_subspace]) for v in VS.basis()]).transpose()

Then, we create a function that will emulate the affine projection (notice the addition of the first element of the affine basis).

sage: def my_proj(x): return proj_matrix*x + affine_basis[0]

Then, we project the polytope and use the transformed coordinates of the vertices to create a new polytope.

sage: proj_p = Projection(P,proj=my_proj)
sage: projected_P = Polyhedron(vertices=proj_p.transformed_coords)

Now, we can see the difference between

sage: P
A 3-dimensional polyhedron in ZZ^3 defined as the convex hull of 4 vertices
sage: P.vertices()
(A vertex at (-2, 0, 3),
 A vertex at (0, -2, 1),
 A vertex at (0, 2, 1),
 A vertex at (2, 0, 3))

and the resulting polytope:

sage: projected_P
A 2-dimensional polyhedron in QQ^3 defined as the convex hull of 4 vertices
sage: projected_P.vertices()
(A vertex at (0, -2, 2),
 A vertex at (0, 0, 2),
 A vertex at (2, 2, 2),
 A vertex at (2, 4, 2))
2017-10-10 12:55:27 -0600 commented question Projection along affine hull

Could we have a precise definition of what is meant by "projection of P along the affine hull"? This would make the question more precise.

2017-10-09 10:32:59 -0600 answered a question Obtaining lattices quickly in SAGE

It is not usual to bundle a database of combinatorial objects coming with Sage by default.

It would be nice to have a (fast) enumerated set of all lattices on a given number of elements. Although my impression is that the set of lattices on a given number of vertices is not easy to generate recursively. Indeed, the lattice property of a poset does not behave well with respect to inclusion, so a recursive enumeration would not work as smoothly as it does for posets.

That said, depending on what you want to do, you could use an iterator:

n_elmt = 7
connected_lattice_iterator = (p for p in Posets(n_elmt) if p.is_connected() and p.is_lattice())

to generate all connected lattices one by one, but it does not increase the speed of generation...

2017-10-09 09:14:57 -0600 answered a question force 0/1 in V-representation of a set of inequalities

From what I understand of your question, you could scale the polyhedron by the lcm of the denominators of the coordinates of the vertices. This will give you a polyhedron with integer vertices.

Nevertheless, in general, it is not true that every (combinatorial type of) polyhedron object has a 0/1 V-representation. So your example should have a special structure that allows this to be possible.

Unfortunately, we can not help further since we don't have an explicit example...

2017-10-09 08:59:14 -0600 answered a question combinatorial equivalence for Polyhedra / isomorphism for lattices


Since version 7.6 of Sagemath, there is exactly the required method for Polyhedron objects:

sage: Poly1=Polyhedron(vertices=[[0,1],[1,0],[1,1]], base_ring=QQ)
sage: Poly2=Polyhedron(vertices=[[2,0],[2,2],[0,2]], base_ring=QQ)
sage: Poly1.is_combinatorially_isomorphic(Poly2)

It does not requires the computation of the full face lattice (checks vertex/facet adjacency graph isomorphism), hence should be faster than asking for isomorphic face lattices. The algorithm option can still be set to face_lattice to get another certificate option.

2017-10-09 08:59:09 -0600 received badge  Necromancer (source)
2017-10-09 08:59:09 -0600 received badge  Teacher (source)
2017-10-09 08:48:47 -0600 answered a question Visualizing Polyhedron with specified range

First, to remove the arrows, one possible way is to use render_solid. The following code will not draw the rays (which are by default drawn as arrows).

sage: QF.<a> = QuadraticField(2)
sage: dd = Polyhedron(ieqs=[[0,1,-a,0],[0,1,a,0],[0,0,-a,1],[0,0,a,1]])
sage: dd.render_solid()

Second, to get a bigger picture, there is a work around by multiplying the rays by a certain scalar.

sage: dd2 = Polyhedron(rays = [100*r.vector() for r in dd.rays()])
sage: dd2.render_solid()

This is not as clean as specifying the ranges on each axis, but is a start.

Currently, to plot the polyhedron object, it uses the data of the polyhedron and it is not (yet) possible to set a certain range for the axis.

2017-10-09 07:34:58 -0600 received badge  Editor (source)
2017-10-09 07:34:08 -0600 answered a question Can the base ring of a polyhedron be restricted?

The current difference between the base ring QQ and ZZ for polyhedron objects in Sage is not quite big nor quite well defined...

A polytope (a compact polyhedron) may have all vertices with integer coordinates, so that it could be called a "Lattice polytope", but its H-representation may still contain rational numbers. Of course, it is possible to multiply all coefficients in the H-representation by a certain integer factor to get an integral H-representation, but this is not automatically done in Sage.

The tickets Add .change_ring() method for polyhedra and Add .change_backend() method for polyhedra will then make it possible to apply the desired changes on the base ring and backend, if possible.

2017-10-09 07:19:15 -0600 answered a question Is there a function to map faces of a cone its dual?

As far as I know, there is no such function in Sage yet. But the following ticket contains a discussion about this issue: Ticket 17215.

It may eventually be implemented...

2017-10-09 06:45:57 -0600 commented answer Generating Tikz codes from Sage for drawing scattering points in 3D

You can change the vertex, edge and facet color directly in the tikz method using the parameters "vertex_color", "edge_color" and "facet_color". They can take any string argument which tikz can interpret as a color. One slight issue with the above loop on integral points is that it may not iterate the interior lattice points according to the projection, i.e. from the furthest to closest to the screen, and maybe some point nodes intersect in a weird way. You could sort them as follows:

proj_vector = vector([775,386,500])
linear_proj = [(proj_vector*p,p) for p in P.integral_points()]
for v in linear_proj:
   tex += ( '\n\\node[inner sep=1pt,circle,draw=red!25!black,fill=red!75!black,thick,anchor=base] at (%s, %s, %s) {};' %tuple(v[1])
2017-03-30 10:26:52 -0600 received badge  Supporter (source)