Closed kcrisman closed 10 years ago
Ping just because Trac doesn't send auto-updates when people attach things for some reason... Thanks, David! Just think how long it will be once there are doctests ;-)
Something to get the ball rolling ...
Attachment: hyperplane.sage.gz
here is a patch, build from the proposed file
it need a lot of work, first to put coverage to 100% and make all tests pass..
for the moment, the coverage is
SCORE sage/geometry/hyperplane_arrangement.py: 61.6% (45 of 73)
and already 16 tests are failing
here is better patch, passing all tests.
Attachment: trac_14789_hyperplane_arrgt.patch.gz
new patch, with coverage reaching 66% and tests passing
A screenshot.
trac_14789_hyp_arrangement_7_29_13.patch is the only patch that should to be applied.
Description changed:
---
+++
@@ -4,3 +4,7 @@
* Linear Algebra
I'm putting this in combinatorics, but I'm not sure where it belongs.
+
+Apply:
+
+* [attachment: trac_14789_hyp_arrangement_7_29_13.patch](https://github.com/sagemath/sage-prod/files/10658015/trac_14789_hyp_arrangement_7_29_13.patch.gz)
for the bot:
apply trac_14789_hyp_arrangement_7_29_13.patch
something is wrong with the hg patch, it's almost empty.
I made a mistake when creating trac_14789_hyp_arrangement_7_29_13.patch (by not using hg add correctly). I have fixed the problem and replaced trac_14789_hyp_arrangement_7_29_13.patch, using the same name. I applied the patch to a fresh clone, just to make sure, and it is seems fine now.
Attachment: trac_14789_hyp_arrangement_7_29_13.patch.gz
First complete version candidate; passes all tests, 100% coverage
Very nice! I like the functionality but have a couple of suggestions to integrate it with the rest of Sage.
There is a lot of needless copy/deepcopy-ing. Most objects in Sage are immutable and don't need to be copied ever, you can just use/return the existing object. E.g. polyhedra are immutable. You should figure out whether you want your objects to be immutable, too. I suggest you make everything immutable, it is generally a big source of errors to make caching work with mutable stuff.
Caching can be done with the @cached_method
decorator, this is much easier to understand/maintain than a __getattr__
hack.
Please no version number for individual files, Sage already has a version number.
Instead of __repr__
, you should implement _repr_
, see http://www.sagemath.org/doc/developer/coding_in_python.html#special-sage-functions
I think its confusing to have a Hyperplane
in the global namespace that is not a Polyhedron
. How about you replace it essentially with the arrangement consisting of a single-hyperplane. So instead of a special add_hyperplane
, just use union:
A = HyperplaneArrangement([[1,0,0],[0,1,1],[0,1,-1],[1,-1,0],[1,1,0]])
B = A.union(HyperplaneArrangement([[1,1,1]])
B = A.union([[1,1,1]]) # syntactic sugar for adding list of hyperplanes
B = A.union([1,1,1]) # syntactic sugar for adding single hyperplane
B = A & [1,1,1] # operator notation
and pick an operator to overload (probably either plus or ampersand) as an alias for union.
The HyperplaneArrangement should be a Parent and individual hyperplanes its elements. Roughly:
class Hyperplane(AffineSubspace, Element):
def __init__(self, parent, ...):
Element.__init__(self, parent=parent)
class HyperplaneArrangement(Parent):
Element = Hyperplane
def __init__(self, ...):
Parent.__init__(self, base=base_ring, category=Sets())
and always use hyperplane_arrangement.element_class()
instead of the plain Hyperplane
class. This automatically gives you a base_ring()
method (Sage generally uses base_ring
instead of base_field
even if it is only implemented for fields).
for the bot:
apply only trac_14789_hyp_arrangement_7_29_13.patch
Thanks very much for these thoughtful comments! I will get to work on them (although I am on hiatus for a week+ to attend a conference).
Sure, I'm at a conference right now ;-) Implementing the Parent / Element thing the first time is relatively difficult, don't hesitate to ask if you run into any problems.
This is a wonderful addition to Sage. Just a couple comments (not a review, not even close):
sage: b = HyperplaneArrangement([])
---------------------------------------------------------------------------
IndexError Traceback (most recent call last)
<ipython-input-5-df3005fcbc9c> in <module>()
----> 1 b = HyperplaneArrangement([])
/home/darij/sage-5.11.beta3/local/lib/python2.7/site-packages/sage/geometry/hyperplane_arrangement.pyc in __init__(self, A, K, check_for_duplicates)
1398 self._hyperplanes = [Hyperplane(h,K) for h in eqns]
1399 self._base_field = K
-> 1400 self._dim = len(eqns[0])-1 # we trust the user on this point
1401 self._ambient_space = VectorSpace(K,self._dim)
1402
IndexError: list index out of range
sage:
The reference [RS] is available online ( http://math.mit.edu/~rstan/arrangements/arr.html ), and you should also correct its title (it's "An Introduction to Hyperplane Arrangements").
"The default base field is QQ": when you write something like that in a docstring, there is no reason not to put the QQ in double-backticks (QQ
).
"A hyperplane arrangement is essential is the normals to its hyperplane span " the "is" should be an "if".
" Class for an affine space (a translation of a linear subspace). " You want to say "affine subspace", not "affine space"!
Not sure if this is a good thing:
sage: a = HyperplaneArrangement([[2,4,6],[1,0,1]])
sage: b = HyperplaneArrangement([[1,0,1],[2,4,6]])
sage: a == b
True
sage: a.sign_vector((2, -2))
[-1, 1]
sage: b.sign_vector((2, -2))
[1, -1]
This, of course, is hard to avert unless you either make equality much more restrictive or use hyperplane arrangements that are lists, rather than sets, of hyperplanes. The latter might in fact be more useful in the long run -- e. g., what if I need to keep track of what happens to each hyperplane when I am changing base field, and I cannot assume that the hyperplanes keep their indices in the .hyperplanes()
list because some hyperplanes might collapse on base change? I assume restriction to a hyperplane would lead to the same problem.
Double "is" in "A face is is the".
The [GZ] reference, like any other Transactions paper, is online: http://www.ams.org/journals/tran/1983-280-01/S0002-9947-1983-0712251-1/
In the definitions of the Catalan and braid arrangements (in their respective docstrings), replace i \leq j by i < j.
Replying to @darijgr: [...]
- Not sure if this is a good thing:
sage: a = HyperplaneArrangement([[2,4,6],[1,0,1]]) sage: b = HyperplaneArrangement([[1,0,1],[2,4,6]]) sage: a == b True
[...]
to me, this is an indication of an issue that has to be carefully sorted out: namely, a notion of equality vs isomorphism. Especially as one might think about creating arrangements from other data, e.g. often in my work I see arrangements obtained from a set of points by taking all the possible hyperplanes on them.
Until this is sorted out, I'd rather think of disabling ==
all together.
A hyperplane arrangement is a set of hyperplanes. Two hyperplane arrangements are equal if their set of hyperplanes are equal. So "==" perfectly clear and useful as it is.
The issue of isomorphism is separable and it would be nice to implement this. There can be an isomorphism that arises from an affine change of coordinates of the ambient space, and there can be a combinatorial isomorphism, meaning an isomorphism of intersection posets (e.g., an arrangement and its essentialization).
Could you be more precise about the other type of data that should determine a hyperplane arrangement? The input would be a set of points. From this set, how does one create a hyperplane arrangement?
Thanks very much for you comments.
Replying to @dimpase:
Replying to @darijgr: [...]
- Not sure if this is a good thing:
sage: a = HyperplaneArrangement([[2,4,6],[1,0,1]]) sage: b = HyperplaneArrangement([[1,0,1],[2,4,6]]) sage: a == b True
[...]
to me, this is an indication of an issue that has to be carefully sorted out: namely, a notion of equality vs isomorphism. Especially as one might think about creating arrangements from other data, e.g. often in my work I see arrangements obtained from a set of points by taking all the possible hyperplanes on them.
Until this is sorted out, I'd rather think of disabling
==
all together.
Sorry. I meant to write: "Two hyperplane arrangements are equal if their sets of hyperplanes are equal. So "==" is perfectly clear and useful as it is." I am traveling at the moment and a bit distracted.
Replying to @sagetrac-dperkinson:
A hyperplane arrangement is a set of hyperplanes. Two hyperplane arrangements are equal if their set of hyperplanes are equal. So "==" perfectly clear and useful as it is.
The issue of isomorphism is separable and it would be nice to implement this. There can be an isomorphism that arises from an affine change of coordinates of the ambient space, and there can be a combinatorial isomorphism, meaning an isomorphism of intersection posets (e.g., an arrangement and its essentialization).
Could you be more precise about the other type of data that should determine a hyperplane arrangement? The input would be a set of points. From this set, how does one create a hyperplane arrangement?
Thanks very much for you comments.
Replying to @dimpase:
Replying to @darijgr: [...]
- Not sure if this is a good thing:
sage: a = HyperplaneArrangement([[2,4,6],[1,0,1]]) sage: b = HyperplaneArrangement([[1,0,1],[2,4,6]]) sage: a == b True
[...]
to me, this is an indication of an issue that has to be carefully sorted out: namely, a notion of equality vs isomorphism. Especially as one might think about creating arrangements from other data, e.g. often in my work I see arrangements obtained from a set of points by taking all the possible hyperplanes on them.
Until this is sorted out, I'd rather think of disabling
==
all together.
Replying to @sagetrac-dperkinson:
A hyperplane arrangement is a set of hyperplanes.
This needs to be said in the docs then. Indeed, currently one sees:
String Form:<class 'sage.geometry.hyperplane_arrangement.HyperplaneArrangement'>
File: /usr/local/src/sage/sage-5.11.beta3/local/lib/python2.7/site-packages/sage/geometry/hyperplane_arrangement.py
Docstring:
The argument "A" is normally a list of hyperplanes or a list of
lists representing hyperplanes...
Could you be more precise about the other type of data that should determine a hyperplane arrangement? The input would be a set of points. From this set, how does one create a hyperplane arrangement?
Just take all the hyperplanes one obtains as affine spans of d-subsets of these points in d-dimensional space.
I have also seen papers talking about the arrangement associated with a triangulation (there one takes the facets of the simplices in the triangulation as hyperplanes).
I think it would be admissible for ==
to return True
even if the hypersurface equations differ by an overall factor. But it might still be a good idea to pick some convention where the hypersurface equation is unique, e.g. scaling the first non-zero coordinate (or maybe largest absolute value) to one. This would also speed up comparison.
On the other hand, ==
should not check whether two arrangements are isomorphic by a GL rotation, or combinatorially, or any other version of a not necessarily unique isomorphism.
I think the original comment was that it is a little odd that you can have a==b
return True while a.sign_vector((2, -2))
returns something different from b.sign_vector((2, -2))
. If they are really sets of hyperplanes, then the order in which those hyperplanes are specified should have no impact. So are they sets, or are they lists?
(Could you add the new file to the reference manual ?)
Just as an FYI (not necessarily needed to implement this now!) we now have matroids in Sage at #7477. So an enhancement request could be to have an arrangement return its associated matroid - and, for that matter, vice versa, if I understand correctly that matroids may be interpreted thus under some circumstances (but not related to the matroid of a central arrangement). Anyway, whatever is right, we might want to do it someday.
I think the original comment was that it is a little odd that you can have
a==b
return True whilea.sign_vector((2, -2))
returns something different fromb.sign_vector((2, -2))
. If they are really sets of hyperplanes, then the order in which those hyperplanes are specified should have no impact. So are they sets, or are they lists?
+1 to the question
However, in this case maybe (maybe) it's okay. Because the question is about "which side of the hyperplane the point is on" and that of course requires an orientation to each hyperplane, given by the equation, and naturally this list will depend on the order. However, the documentation probably should include this example to point out that it definitely depends on the order you give the hyperplanes in, just in case some unsuspecting user might think otherwise. A sign set probably wouldn't be of very much use.
By the way, I assume that we don't get
a is b
True
which I think would be a bug.
wrong plot
Attachment: plot4.png
There is something strange happening with plotting the following
HyperplaneArrangement([[1, 2, 5], [1, -2, -3], [3, -1, 1], [0, 1, 1], [2, -3, 3], [2, 1, -1]])
See the attached plot. Some lines are cut short...
Please try:
sage: a = HyperplaneArrangement([[1, 2, 5], [1, -2, -3], [3, -1, 1], [0, 1, 1], [2, -3, 3], [2, 1, -1]]) sage: a.show(ranges=5)
For finer control of the ranges, please see the "ranges" option for a.plot and a.show. Try a.plot? or a.show? . Perhaps we could think about computing the "optimal" ranges for an arbitrary arrangement.
Replying to @dimpase:
There is something strange happening with plotting the following
HyperplaneArrangement([[1, 2, 5], [1, -2, -3], [3, -1, 1], [0, 1, 1], [2, -3, 3], [2, 1, -1]])
See the attached plot. Some lines are cut short...
Replying to @sagetrac-dperkinson:
Please try:
sage: a = HyperplaneArrangement([[1, 2, 5], [1, -2, -3], [3, -1, 1], [0, 1, 1], [2, -3, 3], [2, 1, -1]]) sage: a.show(ranges=5)
the layout is still weird, as I'd expect a box to be selected and each line extending to the boundaries of the box. Currently it seems that each line has its own box, and they don't play together well.
Yes. Each hyperplane is parametrized. The "ranges" option allows the plotting parameters to be adjusted individually, e.g.,
sage: a.show(ranges=[[-10,6],[-3,10],[-3,3],[-6,6],[-6,6],[-5,5]])
Perhaps there is a way to make the parameters all large, then constrict the bounding box.
Replying to @dimpase:
Replying to @sagetrac-dperkinson:
Please try:
sage: a = HyperplaneArrangement([[1, 2, 5], [1, -2, -3], [3, -1, 1], [0, 1, 1], [2, -3, 3], [2, 1, -1]]) sage: a.show(ranges=5)
the layout is still weird, as I'd expect a box to be selected and each line extending to the boundaries of the box. Currently it seems that each line has its own box, and they don't play together well.
Replying to @sagetrac-dperkinson:
Yes. Each hyperplane is parametrized. The "ranges" option allows the plotting parameters to be adjusted individually, e.g.,
sage: a.show(ranges=[[-10,6],[-3,10],[-3,3],[-6,6],[-6,6],[-5,5]])
Perhaps there is a way to make the parameters all large, then constrict the bounding box.
Inspect the current layout, hyperplane by hyperplane, and find the box by computing maxima and minima of coordinates. Then adjust each hyperplane to fit in the box.
Well, I'm not saying this must happen on this ticket, but please open a followup ticket to address this.
This looks like a bug in num_regions()
:
(stumbled upon this example in a large set of 2d arrangements)
sage: a=HyperplaneArrangement([[1, 106, 106266],[83, 101, 157866],[111, 110, 186150], [453, 221, 532686], [407, 237, 516882],
[55, 32, 75620], [221, 114, 289346], [452, 115, 474217], [406, 131, 453521], [28, 9, 32446], [287, 19, 271774], [241, 35, 244022],
[231, 1, 210984], [185, 17, 181508], [23, -8, 16609]])
sage: a.num_regions()
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-3-cf364604935f> in <module>()
----> 1 a.num_regions()
/home/dima/sage/sage-5.12.beta4/local/lib/python2.7/site-packages/sage/geometry/hyperplane_arrangement.pyc in num_regions(self)
2283 if self.base_field().characteristic() != 0:
2284 raise TypeError('base field must have characteristic zero')
-> 2285 return (-1)**self._dim*self.characteristic_polynomial(-1)
2286
2287 def num_hyperplanes(self):
/home/dima/sage/sage-5.12.beta4/local/lib/python2.7/site-packages/sage/geometry/hyperplane_arrangement.pyc in characteristic_polynomial(self, a)
2190 return self._characteristic_polynomial
2191 else:
-> 2192 return self._characteristic_polynomial.subs(x=a)
2193
2194 def poincare_polynomial(self, a=None):
[...]
/home/dima/sage/sage-5.12.beta4/local/lib/python2.7/site-packages/sage/geometry/hyperplane_arrangement.pyc in restriction(self, H)
3003 for T in new_hyperplanes:
3004 basis = T.linear_part().matrix()
-> 3005 p = H._isomorphism_with_Kn(T.point())
3006 if basis.nrows() == 0:
3007 eqn = [1] + list(p)
/home/dima/sage/sage-5.12.beta4/local/lib/python2.7/site-packages/sage/geometry/hyperplane_arrangement.pyc in _isomorphism_with_Kn(self, v)
813 q = q.apply_map(lambda x: QQ(round(RR(x),2)))
814 u = q*G # u is the approximately orthonormalized W
--> 815 return u.solve_left(v-self.point())
816 else: # finite field
817 p = self.point()
/home/dima/sage/sage-5.12.beta4/local/lib/python2.7/site-packages/sage/matrix/matrix2.so in sage.matrix.matrix2.Matrix.solve_left (sage/matrix/matrix2.c:4386)()
/home/dima/sage/sage-5.12.beta4/local/lib/python2.7/site-packages/sage/matrix/matrix2.so in sage.matrix.matrix2.Matrix.solve_right (sage/matrix/matrix2.c:5242)()
/home/dima/sage/sage-5.12.beta4/local/lib/python2.7/site-packages/sage/matrix/matrix2.so in sage.matrix.matrix2.Matrix._solve_right_general (sage/matrix/matrix2.c:6182)()
ValueError: matrix equation has no solutions
sage:
I'm still not particularly happy with the syntax. Especially handing around tuples isn't really that userfriendly. I would propose the following:
sage: H.<x,y,z> = HyperplaneArrangements(QQ)
as shortcut for
sage: H = HyperplaneArrangements(QQ, names='x, y, z')
sage: H.inject_variables()
and then construct hyperplanes as
sage: hyperplane = x + 2*y + 3*z - 4
and hyperplane arrangements as
sage: arrangement = (x + 2*y + 3*z - 4) and (y - 7) and (x + z)
sage: arrangement == x + 2*y + 3*z - 4 and y - 7 and x + z # brackets are optional
True
If nobody objects then I'll work on that since I kind of need hyperplane arrangements now...
The syntax you suggest is very nice for very small examples. Using tuples is convenient for larger or for computer-generated examples. (It's similar to the syntax for polyhedra in Sage, of course.) So I would support your suggestion as long as it does not interfere with using tuples.
Thanks!
Yes, agree. That would be natural using conversion syntax:
sage: H.<x,y,z> = HyperplaneArrangements(QQ)
sage: H([4, 1, 2, 3])
x + 2*y + 3*z + 4
On a related note, I think the constant should be first for consistency with polyhedra:
sage: Polyhedron(ieqs=[(4,1,2,3)]).Hrepresentation()
(An inequality (1, 2, 3) x + 4 >= 0,)
Replying to @vbraun:
sage: arrangement = (x + 2*y + 3*z - 4) and (y - 7) and (x + z) sage: arrangement == x + 2*y + 3*z - 4 and y - 7 and x + z # brackets are optional True
If nobody objects then I'll work on that since I kind of need hyperplane arrangements now...
Don't use "and". It explicitly has shortcutting behaviour specified on it:
sage: 1 and 2
2
sage: a = 1; b = 2
sage: b if a else a #it's specified to do this
2
so departing from this behaviour would be very surprising.
The operator possibly available would be the "bitwise" &
, but I doubt you want to repurpose it for union.? Being smart with overloading operators often comes back to bite you (see %
for string formatting in python).
Also, since these are hyperplanes given by equations describing them, and
could easily be interpreted as intersection. That's the most direct interpretation coming from the boolean meaning of the operator in sage.
Multiplication would come closer if you think about the arrangement as a union of hyperplanes. Do you keep multiplicities?
If you think of the hyperplane arrangement as a collection of hyperplanes, then they should probably be modelled as a set of hyperplanes, and python does not overload any binary operation to take unions of sets.
Good point, forgot about the short-circuiting. Of course one could override __nonzero__
but I'd rather not if it can be avoided.
Multiplication has the wrong precedence, so thats out. Also, I don't think we should keep mulitplicities.
Bit ops are probably best. How about bit-or = |
:
sage: arrangement = x + 2*y + 3*z - 4 | y - 7 | x + z
Replying to @vbraun:
Bit ops are probably best. How about bit-or =
|
:
Out of the binary ops available that's probably the one that makes most sense. I don't think overloading binops has much precedent in sage/python, though, so it looks a little odd to me. Given how non-compact this notation is anyway, do we really need a compact way of combining hyperplanes?
There's another small problem: vector spaces in sage do not have variable names. In principle, given a vector space, there is a well-defined set of all hyperplane arrangements in that space. If we insist on giving names to the standard basis of the dual vector space, there would be multiple (determined by the names) and the user would be forced to come up with names if he/she wants to construct the collection.
I agree that describing a hyperplane in n-space with a tuple of length n+1 feels a little uncomfortable: It should really be a vector in the dual space together with a scalar, i.e., a pair of an n-tuple and a scalar, which is awkward to work with by itself.
Would it be better to let (a,b,c,d) mean ax+by+c*z+d=0 (note the sign on the constant)? That's a little more like standard homogenization.
Precedent in other software dealing with hyperplane arrangements might be instructive?
I'm working on making (d,a,b,c) mean ax+by+c*z+d=0. So the sign is expected. And the order as well as the sign matches the notation used for polyhedra. If you want to avoid ambiguity then list/tuple plus scalar also define a hyperplane::
sage: H.<x,y> = HyperplaneArrangements(QQ)
sage: H(x, y, x+y-1)
Arrangement <y | x | x + y - 1>
sage: H([0,1,0], [0,0,1], [-1,1,1])
Arrangement <y | x | x + y - 1>
sage: H([(1,1),0], [(2,3),-1], [(4,5),3])
Arrangement <x + y | 2*x + 3*y - 1 | 4*x + 5*y + 3>
Branch: u/vbraun/hyperplane_arrangements
New commits:
[changeset:2ba6550] | Ported hyperplanes and arrangements to linear expressions |
[changeset:e2f4d13] | implemented linear expressions |
[changeset:d8cb75e] | split up into manageable source files |
[changeset:b528d68] | Trac 14789: hyperplane arrangement package |
Almost finished... here is Dima's example:
sage: H.<x,y> = HyperplaneArrangements(QQ)
sage: h = H([(1, 106), 106266], [(83, 101), 157866], [(111, 110), 186150], [(453, 221), 532686],
....: [(407, 237), 516882], [(55, 32), 75620], [(221, 114), 289346], [(452, 115), 474217],
....: [(406, 131), 453521], [(28, 9), 32446], [(287, 19), 271774], [(241, 35), 244022],
....: [(231, 1), 210984], [(185, 17), 181508], [(23, -8), 16609])
sage: h.n_regions()
85
Replying to @vbraun:
Almost finished... here is Dima's example:
... sage: h.n_regions() 85
Has it magically got fixed? ;–)
Replying to @dimpase:
Has it magically got fixed? ;–)
All sufficiently hard work is indistinguishable from magic ;-)
In case I haven't made it clear, the branch contains quite extensive modifications of the original patch.
Branch pushed to git repo; I updated commit sha1. New commits:
[changeset:f7c4b81] | added remaining functionality, fixed documentation |
This is an active area of research, closely connected to
CC: @sagetrac-dperkinson @sagetrac-Stefan @nathanncohen @rbeezer @vbraun
Component: combinatorics
Keywords: days54
Author: David Perkinson, Volker Braun
Branch/Commit: public/geometry/14789-hyperplane_arrangements @
ba798c0
Reviewer: Travis Scrimshaw
Issue created by migration from https://trac.sagemath.org/ticket/14789