attachment:hyperplane.sage of education5

   1 r"""
   2 AUTHOR:
   3     -- David Perkinson (2013-6-21)
   4 
   5 VERSION: -1
   6 
   7 EXAMPLES: (see HArrangement, below)
   8 
   9 #*****************************************************************************
  10 #       Copyright (C) 2013 David Perkinson <davidp@reed.edu>
  11 #
  12 #  Distributed under the terms of the GNU General Public License (GPL)
  13 #                  http://www.gnu.org/licenses/
  14 #*****************************************************************************
  15 """
  16 
  17 class AffSpace():
  18     r"""
  19     Class for an affine space (a translation of a linear subspace).
  20 
  21     INPUT:
  22 
  23     - ``W`` -- vector space
  24     - ``p`` -- list representing a point in W
  25 
  26     OUTPUT:
  27 
  28     - AffSpace
  29 
  30     EXAMPLES:
  31 
  32         sage: a = AffSpace([1,0,0,0],VectorSpace(QQ,4))
  33         sage: a.dim()
  34         4
  35         sage: a.point()
  36         (1, 0, 0, 0)
  37         sage: a.linear_part()
  38         Vector space of dimension 4 over Rational Field
  39         sage: a
  40         Affine space p + W where:
  41            p = [1, 0, 0, 0]
  42            W = Vector space of dimension 4 over Rational Field.
  43         sage: b = AffSpace((1,0,0,0),matrix(QQ,[[1,2,3,4]]).right_kernel())
  44         sage: c = AffSpace((0,2,0,0),matrix(QQ,[[0,0,1,2]]).right_kernel())
  45         sage: b.intersection(c)
  46         Affine space p + W where:
  47            p = (-3, 2, 0, 0)
  48            W = Vector space of degree 4 and dimension 2 over Rational Field
  49         Basis matrix:
  50         [  1   0  -1 1/2]
  51         [  0   1  -2   1].
  52         sage: b < a
  53         True
  54         sage: c < b
  55         False
  56     """
  57 
  58     def __init__(self, p, V):
  59         self._linear_part = V
  60         self._point = vector(V.base_field(),p)
  61 
  62     def __repr__(self):
  63         return "Affine space p + W where:\n   p = "+str(self._point)+"\n   W = "+str(self._linear_part)+"."
  64 
  65     def __eq__(self, other):
  66         V = self._linear_part
  67         W = other._linear_part
  68         if V==W and self._point-other._point in V:
  69             return True
  70         else:
  71             return False
  72 
  73     def __ne__(self, other):
  74         return not self==other
  75 
  76     def __le__(self, other):
  77         V = self._linear_part
  78         W = other._linear_part
  79         if V.is_subspace(W) and self._point-other._point in W:
  80             return True
  81         else:
  82             return False
  83 
  84     def __lt__(self, other):
  85         if self<=other and not self==other:
  86             return True
  87         else:
  88             return False
  89 
  90     def __ge__(self, other):
  91         return other<=self
  92 
  93     def __gt__(self, other):
  94         return other<self
  95 
  96     def __contains__(self,q):
  97         return self._point - q in self._linear_part
  98 
  99     def linear_part(self):
 100         r"""
 101         The linear part of the affine space.
 102 
 103         INPUT:
 104 
 105         - None
 106 
 107         OUTPUT:
 108 
 109         - vector space
 110 
 111         EXAMPLES::
 112 
 113             sage: A = AffSpace([2,3,1], matrix(QQ,[[1,2,3]]).right_kernel()) 
 114             sage: A.linear_part()
 115             Vector space of degree 3 and dimension 2 over Rational Field
 116             Basis matrix:
 117             [   1    0 -1/3]
 118             [   0    1 -2/3]
 119         """
 120         return copy(self._linear_part)
 121 
 122     def point(self):
 123         r"""
 124         A point ``p`` in the affine space.
 125 
 126         INPUT:
 127 
 128         - None
 129 
 130         OUTPUT:
 131 
 132         - vector
 133 
 134         EXAMPLES::
 135 
 136             sage: A = AffSpace([2,3,1], VectorSpace(QQ,3))
 137             sage: A.point()
 138             (2, 3, 1)
 139         """
 140         return copy(self._point)
 141 
 142     def dim(self):
 143         r"""
 144         The dimension of the affine space.
 145 
 146         INPUT:
 147 
 148         - None
 149 
 150         OUTPUT:
 151 
 152         - Integer
 153 
 154         EXAMPLES::
 155 
 156             sage: a = AffSpace([1,0,0,0],VectorSpace(QQ,4))
 157             sage: a.dim()
 158             4
 159         """
 160         return self._linear_part.dimension()
 161 
 162 
 163     def intersection(self,other):
 164         r"""
 165 
 166 
 167         INPUT:
 168 
 169         -
 170 
 171         OUTPUT:
 172 
 173         -
 174 
 175         EXAMPLES::
 176 
 177 
 178         """
 179         m = self._linear_part.matrix()
 180         n = other._linear_part.matrix()
 181         p = self._point
 182         q = other._point
 183         A = m.stack(n)
 184         v = q-p
 185         if A.rank()!=(A.stack(v)).rank():
 186             return -1
 187         else:
 188             t = A.solve_left(v)
 189             new_p = p + t[:m.nrows()]*m
 190             new_V = self._linear_part.intersection(other._linear_part)
 191             return AffSpace(new_p,new_V)
 192 
 193     def isomorphism_with_Kn(self,v):
 194         p = self._point
 195         W = self._linear_part
 196         V = VectorSpace(W.base(),W.dimension())
 197         f = W.hom(V.gens())
 198         return f(v-p)
 199 
 200 
 201 class Hyperplane(AffSpace):
 202     def __init__(self,H,K=QQ):
 203         r"""
 204         The argument `H` is a hyperplane, given as a list ``[a1,...an,a]``
 205         representing the equation ``a1x1 + ... + anxn = a``. An optional field,
 206         `K`, may also be provided.
 207 
 208         INPUT:
 209 
 210         - `H` -- list of integers representing a hyperplane
 211         - `K` -- field (default = `QQ`)
 212 
 213         OUTPUT:
 214 
 215         - Hyperplane
 216 
 217         EXAMPLES::
 218 
 219             sage: h = Hyperplane([1,3,2,2,8])
 220             sage: h
 221             Hyperplane over Rational Field
 222             (1, 3, 2, 2)*X = 8
 223         """
 224         m = matrix(K,H[:-1])
 225         p = m.solve_right(vector([H[-1]]))
 226         AffSpace.__init__(self,p,m.right_kernel())
 227         self._raw_data = H
 228         self._normal = vector(K,H[:-1])
 229         self._field = K
 230 
 231     def __repr__(self):
 232         return "Hyperplane over "+str(self._field)+"\n"+str(self._normal)+"*X = "+str(self._raw_data[-1])
 233 
 234     def equation(self,suppress_printing=False):
 235         if not suppress_printing:
 236             print str(self._field)+"\n"+str(self._normal)+"*X = "+str(self._raw_data[-1])
 237         return deepcopy(self._raw_data)
 238 
 239     def normal(self):
 240         return self._normal
 241 
 242     # orthogonalize target basis to get less distorted picture?
 243     # random color option?
 244     def show(self,r=None,**kwds):
 245         if r==None:
 246             r = 1
 247         if self.dim()==1: # so the dimension of the ambient space is 2
 248             p = self.point()
 249             w = self.linear_part().matrix()
 250             var('t')
 251             parametric_plot(p+t*w[0],(t,-r,r)).show(**kwds)
 252         elif self.dim()==2: # so the dimension of the ambient space is 3
 253             p = self.point()
 254             w = self.linear_part().matrix()
 255             var('s t')
 256             parametric_plot3d(p+s*w[0]+t*w[1],(s,-r,r),(t,-r,r)).show(**kwds)
 257 
 258     def plot(self,r=None,**kwds):
 259         if r==None:
 260             r = 1
 261         if self.dim()==1: # so the dimension of the ambient space is 2
 262             p = self.point()
 263             w = self.linear_part().matrix()
 264             var('t')
 265             return parametric_plot(p+t*w[0],(t,-r,r),**kwds)
 266         elif self.dim()==2: # so the dimension of the ambient space is 3
 267             p = self.point()
 268             w = self.linear_part().matrix()
 269             var('s t')
 270             return parametric_plot3d(p+s*w[0]+t*w[1],(s,-r,r),(t,-r,r),**kwds)
 271 
 272 class HArrangement():
 273 
 274     def __init__(self,A,K=QQ,check_for_duplicates=True):
 275         r"""
 276         The argument `A` is a list of hyperplanes. Each hyperplane is given as a
 277         list ``[a1,...an,a]`` representing the equation ``a1x1 + ... + anxn =
 278         a``. An optional field, `K`, may also be provided.  Multiple copies of a
 279         hyperplane are not supported, and a check is performed by default to
 280         remove duplicates.
 281 
 282         INPUT:
 283 
 284         - `A` -- list
 285         - `K` -- field (default = QQ)
 286         - check_for_duplicates -- boolean
 287 
 288         OUTPUT:
 289 
 290         - HArrangement
 291 
 292         EXAMPLE::
 293 
 294             sage: A = HArrangement(([1,0,0],[2,1,5]))
 295             sage: A
 296             Hyperplane arrangement of 2 hyperplanes over Rational Field of dimension 2, rank 2.
 297             sage: C = B.essentialization()
 298             sage: B.characteristic_polynomial()
 299             x^4 - 6*x^3 + 11*x^2 - 6*x
 300             sage: C.characteristic_polynomial()
 301             x^3 - 6*x^2 + 11*x - 6
 302             sage: G = graphs.CycleGraph(4)
 303             sage: B = hyperplane_arrangements.graphical(G)
 304             sage: B.is_essential()
 305             False
 306             sage: C = B.essentialization()
 307             sage: B.characteristic_polynomial()
 308             x^4 - 4*x^3 + 6*x^2 - 3*x
 309             sage: G.chromatic_polynomial()
 310             x^4 - 4*x^3 + 6*x^2 - 3*x
 311             sage: C.characteristic_polynomial()
 312             x^3 - 4*x^2 + 6*x - 3
 313             sage: B.num_regions()
 314             14
 315             sage: C.num_regions()
 316             14
 317             sage: D = hyperplane_arrangements.semiorder(3)
 318             sage: D.is_essential()
 319             False
 320             sage: D.num_regions()
 321             19
 322             sage: D.num_bounded_regions()
 323             7
 324         """
 325         if hasattr(A,'Hrepresentation'):
 326             self._from_polytope = True  # test to see if there is a polytope for the arrangement?
 327             self._polytope = deepcopy(A)
 328             self._hrepresentation = A.Hrepresentation()
 329             self._raw_data = []
 330             for h in A.Hrepresentation():
 331                 hh = list(h)
 332                 self._raw_data.append(hh[1:]+[-hh[0]])
 333             A = self._raw_data
 334         else:
 335             self._from_polytope = False
 336             self._raw_data = A
 337         if check_for_duplicates:
 338             self._hyperplanes = []
 339             for h in A:
 340                 H = Hyperplane(h,K)
 341                 if any (H==k for k in self._hyperplanes)==False:
 342                     self._hyperplanes.append(H)
 343         else:
 344             self._hyperplanes = [Hyperplane(H,K) for H in A]
 345         self._base_field = K
 346         self._dim = len(A[0])-1
 347         self._ambient_space = VectorSpace(K,self._dim)
 348 
 349     def __getattr__(self, name):
 350 
 351         if not self.__dict__.has_key(name):
 352             if name == '_characteristic_polynomial':
 353                 self._set_characteristic_polynomial()
 354                 return deepcopy(self.__dict__[name])
 355             elif name == '_essentialization':
 356                 self._set_essentialization()
 357                 return self.__dict__[name]
 358             elif name == '_intersection_poset':
 359                 self._set_intersection_poset()
 360                 return deepcopy(self.__dict__[name])
 361             elif name == '_poincare_polynomial':
 362                 self._set_poincare_polynomial()
 363                 return deepcopy(self.__dict__[name])
 364             elif name == '_rank':
 365                 self._set_rank()
 366                 return self.__dict__[name]
 367             elif name == '_regions':
 368                 self._set_regions()
 369                 return self.__dict__[name]
 370             else:
 371                 raise AttributeError, name
 372 
 373     def __repr__(self):
 374         return "Hyperplane arrangement of "+str(len(self._hyperplanes))+ " hyperplanes over "+str(self._base_field)+" of dimension "+str(self._dim)+", rank "+str(self.rank())+"."
 375 
 376     def hyperplanes(self):
 377         r"""
 378         Returns the hyperplanes of the arrangement.
 379 
 380         INPUT:
 381 
 382         - None
 383 
 384         OUTPUT:
 385 
 386         - list of Hyperplanes
 387 
 388         EXAMPLES::
 389 
 390             sage: H = HArrangement([[1,0,0],[0,1,1],[0,1,-1],[1,-1,0],[1,1,0]])
 391             sage: H.hyperplanes()
 392             [Hyperplane over Rational Field
 393             (1, 0)*X = 0,
 394              Hyperplane over Rational Field
 395             (0, 1)*X = 1,
 396              Hyperplane over Rational Field
 397             (0, 1)*X = -1,
 398              Hyperplane over Rational Field
 399             (1, -1)*X = 0,
 400              Hyperplane over Rational Field
 401             (1, 1)*X = 0]
 402         """
 403         return copy(self._hyperplanes)
 404 
 405     def cone(self):
 406         r"""
 407         NOTE: INCOMPLETE
 408 
 409         INPUT:
 410 
 411         - None
 412 
 413         OUTPUT:
 414 
 415         - hyperplane arrangement
 416 
 417         EXAMPLES::
 418 
 419             sage: HA.cone()
 420             Hyperplane arrangement of 3 hyperplanes over Rational Field of dimension 3, rank 3.
 421         """
 422         hyps = []
 423         for h in self.hyperplanes():
 424             new = h.equation(suppress_printing=True)
 425             new[-1] = -new[-1]
 426             new.append(0)
 427             hyps.append(new)
 428         hyps.append([0]*self.dim()+[1,0])
 429         return HArrangement(hyps,self.base_field())
 430 
 431     def dim(self):
 432         return self._dim
 433 
 434     def _set_rank(self):
 435         self._rank= matrix([list(H.normal()) for H in self._hyperplanes]).rank()
 436 
 437     def rank(self):
 438         return self._rank
 439 
 440     def has_good_reduction(self,p):
 441         raise UserWarning, 'not implemented'
 442 
 443     def is_essential(self):
 444         r"""
 445         Indicates whether the hyperplane arrangement is essential. A hyperplane
 446         arrangement is essential if the span of the normals of its hyperplanes
 447         spans the ambient space.
 448 
 449         INPUT:
 450 
 451         - None
 452 
 453         OUTPUT:
 454 
 455         - boolean
 456 
 457         EXAMPLES::
 458 
 459             sage: HArrangement([[1,0,0],[1,0,1]]).is_essential()
 460             False
 461         """
 462         r = matrix([list(H.normal()) for H in self._hyperplanes]).rank()
 463         return r==self._dim
 464 
 465     def _set_intersection_poset(self):
 466         ## initiate with codim = 0, K^n
 467         K = self.base_field()
 468         Kn = AffSpace(vector(K,[0 for i in
 469             range(self._dim)]),self._ambient_space)
 470         L = [[Kn]]
 471         active = True
 472         codim = 0
 473         while active:
 474             active = False
 475             new_level = []
 476             for T in L[codim]:
 477                 for H in self._hyperplanes:
 478                     I = H.intersection(T)
 479                     if type(I)!=type(-1) and I!=T and I not in new_level:
 480                         new_level.append(I)
 481                         active = True
 482             if active:
 483                 L.append(new_level)
 484             codim += 1
 485         L = flatten(L)
 486         t = {}
 487         for i in range(len(L)):
 488             t[i] = L[i]
 489         cmp_fn = lambda p,q: t[q]<t[p]
 490         self._intersection_poset = Poset((t,cmp_fn))
 491 
 492     def intersection_poset(self):
 493         return deepcopy(self._intersection_poset)
 494 
 495     def _set_characteristic_polynomial(self):
 496         #P = self._intersection_poset
 497         #n = self._dim
 498         #self._characteristic_polynomial = sum([P.mobius_function(0,p)*x^(n-P.rank(p)) for p in P])
 499         if self.rank()==1:
 500             self._characteristic_polynomial = x^(self.dim()-1)*(x-self.num_hyperplanes())
 501         else:
 502             H = self.hyperplanes()[0]
 503             R = self.restriction(H)
 504             r = R.characteristic_polynomial()
 505             D = self.deletion(H)
 506             d = D.characteristic_polynomial()
 507             self._characteristic_polynomial = expand(d - r)
 508             return
 509 
 510     def _set_poincare_polynomial(self):
 511         p = (-x)^self.dim()*self.characteristic_polynomial(-1/x)
 512         self._poincare_polynomial = expand(p)
 513 
 514     def change_base_field(self):
 515         raise UserWarning, 'not implemented'
 516 
 517     def characteristic_polynomial(self,a=None):
 518         if a==None:
 519             return self._characteristic_polynomial
 520         else:
 521             return self._characteristic_polynomial.subs(x=a)
 522 
 523     def poincare_polynomial(self,a=None):
 524         if a==None:
 525             return self._poincare_polynomial
 526         else:
 527             return self._poincare_polynomial.subs(x=a)
 528 
 529     def num_bounded_regions(self):
 530         r"""
 531         Returns the number of (relatively) bounded regions of the hyperplane
 532         arrangement as a subset of ``\mathbb{RR}``.
 533 
 534         INPUT:
 535 
 536         - None
 537 
 538         OUTPUT:
 539 
 540         - integer
 541 
 542         EXAMPLES::
 543 
 544             sage: A = hyperplane_arrangements.semiorder(3)
 545             sage: A.num_bounded_regions()
 546             7
 547         """
 548         if self.base_field().characteristic()!=0:
 549             raise UserWarning, 'base field must have characteristic zero'
 550         return (-1)^self.rank()*self.characteristic_polynomial(1)
 551 
 552     def num_regions(self):
 553         r"""
 554         Returns the number of regions of the hyperplane arrangement.
 555 
 556         INPUT:
 557 
 558         - None
 559 
 560         OUTPUT:
 561 
 562         - integer
 563 
 564         EXAMPLES::
 565 
 566             sage: A = hyperplane_arrangements.semiorder(3)
 567             sage: A.num_regions()
 568             19
 569         """
 570         if self.base_field().characteristic()!=0:
 571             raise UserWarning, 'base field must have characteristic zero'
 572         return (-1)^self._dim*self.characteristic_polynomial(-1)
 573 
 574     def num_hyperplanes(self):
 575         return len(self.hyperplanes())
 576 
 577     def ambient_space(self):
 578         r"""
 579         The ambient space of the hyperplane arrangement.
 580 
 581         INPUT:
 582 
 583         - None
 584 
 585         OUTPUT:
 586 
 587         - vector space
 588 
 589         EXAMPLES::
 590 
 591             sage: A = HArrangement([[1,0,0],[0,1,1],[0,1,-1],[1,-1,0],[1,1,0]])
 592             sage: A.ambient_space()
 593             Vector space of dimension 2 over Rational Field
 594         """
 595         return deepcopy(self._ambient_space)
 596 
 597     def base_field(self):
 598         r"""
 599         The base field of the hyperplane arrangement.
 600 
 601         INPUT:
 602 
 603         - None
 604 
 605         OUTPUT:
 606 
 607         - field
 608 
 609         EXAMPLES::
 610 
 611             sage: A = HArrangement([[1,0,0],[0,1,1],[0,1,-1],[1,-1,0],[1,1,0]])
 612             sage: A.base_space()
 613             Rational Field
 614 
 615 	.. SEEALSO::
 616 
 617 		:meth:`.change_base_field`
 618         """
 619         return self._base_field
 620 
 621     def deletion(self,h):
 622         r"""
 623         The hyperplane arrangement obtained by deleting `h` from the hyperplane arrangement.
 624 
 625         INPUT:
 626 
 627         - `h` -- Hyperplane
 628 
 629         OUTPUT:
 630 
 631         - HArrangement
 632 
 633         EXAMPLES::
 634 
 635 	    sage: A = HArrangement([[1,0,0],[0,1,1],[0,1,-1],[1,-1,0],[1,1,0]])
 636 	    sage: h = Hyperplane([1,0,0])
 637 	    sage: A
 638 	    Hyperplane arrangement of 5 hyperplanes over Rational Field of dimension 2, rank 2.
 639 	    sage: A.deletion(h)
 640 	    Hyperplane arrangement of 4 hyperplanes over Rational Field of dimension 2, rank 2.
 641 
 642 
 643 	.. SEEALSO::
 644 
 645 		:meth:`.restriction`
 646         """
 647         if h in self._hyperplanes:
 648             A = [k.equation(True) for k in self.hyperplanes()]
 649             A.remove(h._raw_data)
 650             return HArrangement(A,self.base_field())
 651         else:
 652             raise UserWarning, 'hyperplane is not in the arrangement'
 653 
 654     def _set_essentialization(self):
 655         r"""
 656         The essentialization of the hyperplane arrangement.  See the
 657         documentation for ``essentialization``.
 658         """
 659         if self.rank()==self.dim():
 660             raise UserWarning, 'already essential'
 661         else:
 662             n = [h.normal() for h in self._hyperplanes]
 663             N = self.ambient_space().subspace(n)
 664             # we need to be careful finding complementary spaces when the
 665             # characteristic is not 0
 666             if N.base_field().characteristic()!=0:
 667                 m = N.echelonized_basis_matrix()
 668                 c = [m.nonzero_positions_in_row(i)[0] for i in range(m.nrows())]
 669                 new_basis = [self.ambient_space().basis()[i] for i in c]
 670                 N = self.ambient_space().subspace(new_basis)
 671             V = VectorSpace(self.base_field(),N.dimension())
 672             N = AffSpace(vector([0]*self.dim()),N)
 673             new_hyperplanes = []
 674             for h in self.hyperplanes():
 675                 I = h.intersection(N)
 676                 basis = I.linear_part().matrix()
 677                 p = N.isomorphism_with_Kn(I.point())
 678                 new_basis = matrix(self.base_field(),[N.isomorphism_with_Kn(b) for b in basis])
 679                 if new_basis.nrows()==0:
 680                     eqn = [1,p[0]]
 681                 else:
 682                     eqn = new_basis.right_kernel_matrix()[0]
 683                     q = [eqn*p]
 684                     eqn = list(eqn)+q
 685                 new_hyperplanes.append(eqn)
 686             self._essentialization = HArrangement(new_hyperplanes,self.base_field())
 687 
 688     def essentialization(self):
 689         r"""
 690         The essentialization of a hyperplane arrangement.  If the characteristic
 691         of the base field is 0, this returns the hyperplane arrangement obtained
 692         by intersecting the hyperplanes by the space spanned by their normal
 693         vectors.
 694 
 695         INPUT:
 696 
 697         - None
 698 
 699         OUTPUT:
 700 
 701         - HArrangement
 702 
 703         EXAMPLES::
 704 
 705 	    sage: A = hyperplane_arrangements.braid(3)
 706 	    sage: A.is_essential()
 707 	    False
 708 	    sage: A.essentialization()
 709 	    Hyperplane arrangement of 3 hyperplanes over Rational Field of dimension 2, rank 2.
 710 	    sage: B = HArrangement([[1,0,1],[1,0,-1]])
 711 	    sage: B.is_essential()
 712 	    False
 713 	    sage: B.essentialization()
 714 	    Hyperplane arrangement of 2 hyperplanes over Rational Field of dimension 1, rank 1.
 715 	    sage: B.essentialization().hyperplanes()
 716 	    [Hyperplane over Rational Field
 717 	    (1)*X = 1,
 718 	     Hyperplane over Rational Field
 719 	    (1)*X = -1]
 720             sage: C = HArrangement([[1,1,1],[1,1,0]],GF(2))
 721             sage: C.essentialization().hyperplanes()
 722             [Hyperplane over Finite Field of size 2
 723             (1)*X = 1,
 724              Hyperplane over Finite Field of size 2
 725             (1)*X = 0]
 726         """
 727         return self._essentialization
 728 
 729     def polytope(self):
 730         if self._from_polytope:
 731             return deepcopy(self._polytope)
 732 
 733     def _set_regions(self):
 734         r"""
 735         Computes the regions of the hyperplane arrangement.  See the
 736         documentation for ``regions``.
 737         """
 738         if self.base_field().characteristic()!=0:
 739             raise UserWarning, 'base field must have characteristic zero'
 740         result = []
 741         num = self.num_hyperplanes()
 742         h = [i.equation(suppress_printing=True) for i in self.hyperplanes()]
 743         for pos in powerset(range(num)):
 744             q = []
 745             for i in range(num):
 746                 if i in pos:
 747                     new = [-h[i][-1]]+h[i][:-1]
 748                     q.append(new)
 749                 else:
 750                     new = [h[i][-1]]+[-t for t in h[i][:-1]]
 751                     q.append(new)
 752             P = Polyhedron(ieqs=q)
 753             if P.dim()==self.dim():
 754                 result.append(P)
 755         self._regions = result
 756 
 757     def regions(self):
 758         r"""
 759         The regions of the hyperplane arrangement, i.e., the connected components of the
 760         complement of the union of the hyperplanes as a subset of `\mathbb{R}^n`.
 761 
 762         INPUT:
 763 
 764         - None
 765 
 766         OUTPUT:
 767 
 768         - a list of polyhedra
 769 
 770         EXAMPLES::
 771 
 772             sage: A = HArrangement([[1,0,0],[0,1,1]])
 773             sage: A.regions()
 774             [[A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 1 vertex and 2 rays],
 775             [A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 1 vertex and 2 rays],
 776             [A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 1 vertex and 2 rays],
 777             [A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 1 vertex and 2 rays]]
 778         """
 779         return self._regions
 780 
 781     def region_containing_point(self,p):
 782         r"""
 783         The region in the hyperplane arrangement containing a given point.
 784 
 785         INPUT:
 786 
 787         - a point
 788 
 789         OUTPUT:
 790 
 791         - a polyhedron
 792 
 793         EXAMPLES::
 794 
 795             sage: A = HArrangement([[1,0,0],[0,1,1],[0,1,-1],[1,-1,0],[1,1,0]])
 796             sage: A.region_containing_point((1,2))
 797             A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 2 vertices and 2 rays
 798         """
 799         if self.base_field().characteristic()!=0:
 800             raise UserWarning, 'base field must have characteristic zero'
 801         h = [i.equation(suppress_printing=True) for i in self.hyperplanes()]
 802         q = []
 803         sign_vector = self.sign_vector(p)
 804         if 0 in sign_vector:
 805             raise UserWarning, 'point sits on a hyperplane'
 806             return
 807         for i in range(self.num_hyperplanes()):
 808             if sign_vector[i]==1:
 809                 new = [-h[i][-1]]+h[i][:-1]
 810                 q.append(new)
 811             else:
 812                 new = [h[i][-1]]+[-t for t in h[i][:-1]]
 813                 q.append(new)
 814         return Polyhedron(ieqs=q)
 815 
 816     def sign_vector(self,p):
 817         r"""
 818         Sign_vector indicates on which side of each hyperplane the given point `p` lies.
 819 
 820         INPUT:
 821 
 822         - point
 823 
 824         OUTPUT:
 825 
 826         - vector
 827 
 828         EXAMPLES::
 829 
 830             sage: A = HArrangement([[1,0,0],[0,1,1]])
 831             sage: A.sign_vector((2,-2))
 832             [1, -1]
 833             A.sign_vector((1,1))
 834             [1, 0]
 835         """
 836         if self.base_field().characteristic()!=0:
 837             raise UserWarning, 'characteristic must be zero'
 838         else:
 839             p = vector(p)
 840             eqns = [i.equation(suppress_printing=True) for i in self.hyperplanes()]
 841             return [sign(p*vector(e[:-1])-e[-1]) for e in eqns]
 842 
 843     def restriction(self,H):
 844         r"""
 845         The restriction of the hyperplane arrangement to the given hyperplane, `H`.
 846         This is obtained by intersecting  `H` with all of the other distinct
 847         hyperplanes.  The hyperplane `H` must be in the arrangement.
 848 
 849 
 850         INPUT:
 851 
 852         - a Hyperplane
 853 
 854         OUTPUT:
 855 
 856         - an HArrangement
 857 
 858         EXAMPLES::
 859 
 860 	    sage: A = hyperplane_arrangements.braid(4)
 861 	    sage: H = A.hyperplanes()[0]
 862 	    sage: H
 863 	    Hyperplane over Rational Field
 864 	    (1, -1, 0, 0)*X = 0
 865 	    sage: R = A.restriction(H)
 866 	    sage: D = A.deletion(H)
 867 	    sage: ca = A.characteristic_polynomial()
 868 	    sage: cr = R.characteristic_polynomial()
 869 	    sage: cd = D.characteristic_polynomial()
 870 	    sage: ca
 871 	    x^4 - 6*x^3 + 11*x^2 - 6*x
 872 	    sage: cd - cr
 873 	    x^4 - 6*x^3 + 11*x^2 - 6*x
 874 
 875 	.. SEEALSO::
 876 
 877 		:meth:`.deletion`
 878         """
 879         if H not in self.hyperplanes():
 880             raise UserWarning, 'hyperplane not in arrangement'
 881         else:
 882             new_hyperplanes = []
 883             for T in self.hyperplanes():
 884                 if T!=H:
 885                     I = H.intersection(T)
 886                     if type(I)!=Integer and I not in new_hyperplanes:
 887                         new_hyperplanes.append(I)
 888         final_hyperplanes = []
 889         for T in new_hyperplanes:
 890             basis = T.linear_part().matrix()
 891             p = H.isomorphism_with_Kn(T.point())
 892             if basis.nrows()==0:
 893                 eqn = [1] + list(p)
 894             else:
 895                 new_basis = matrix(self.base_field(),[H.isomorphism_with_Kn(b+H.point())
 896                     for b in basis])
 897                 eqn = new_basis.right_kernel_matrix()[0]
 898                 q = list(eqn*p)
 899                 eqn = list(eqn)+q
 900             final_hyperplanes.append(eqn)
 901         return HArrangement(final_hyperplanes,self.base_field())
 902 
 903     def show(self,r=None,**kwds):
 904         if r==None:
 905             r = 1
 906         if self.dim()==2 or self.dim()==3: # so the dimension of the ambient space is 3
 907             sum([h.plot(r) for h in self.hyperplanes()]).show(**kwds)
 908 
 909     def union(self,A):
 910         r"""
 911         The union of one HArrangement with another.
 912 
 913         INPUT:
 914 
 915         - `A` -- HArrangement
 916 
 917         OUTPUT:
 918 
 919         - HArrangement
 920 
 921         EXAMPLES::
 922 
 923             sage: A = HArrangement([[1,0,0],[0,1,1],[0,1,-1],[1,-1,0],[1,1,0]])
 924             sage: B = HArrangement([[1,1,1],[1,-1,1],[1,0,-1]])
 925             sage: A.union(B)
 926 	    Hyperplane arrangement of 8 hyperplanes over Rational Field of dimension 2, rank 2.
 927         """
 928         hyps = [i.equation(true) for i in self.hyperplanes()] + [i.equation(true)
 929                 for i in A.hyperplanes()]
 930         return HArrangement(hyps,self.base_field())
 931 
 932 # Idea: calculate characteristic polynomials and set
 933 # self._characteristic_polynomial while constructing the special hyperplane
 934 # arrangement.  Also set poset?
 935 class HArrangementGenerators():
 936 
 937     def semiorder(self,n,K=QQ):
 938         r"""
 939         The semiorder arrangement is the set of `n(n-1)` hyperplanes: `\{ x_i - x_j = -1,1 : 1\leq i \leq j\leq n\}.`
 940 
 941         INPUT:
 942 
 943         - `n` -- integer
 944         - `K` (default `QQ`)
 945 
 946         OUTPUT:
 947 
 948         - HArrangement
 949 
 950         EXAMPLES::
 951 
 952             sage: hyperplane_arrangements.semiorder(4)
 953             Hyperplane arrangement of 12 hyperplanes over Rational Field of dimension 4, rank 3.
 954         """
 955         hyperplanes = []
 956         for i in range(n):
 957             for j in range(i+1,n):
 958                 new = [0]*n
 959                 new[i] = 1
 960                 new[j] = -1
 961                 for k in [-1,1]:
 962                     h = deepcopy(new)
 963                     h.append(k)
 964                     hyperplanes.append(h)
 965         return HArrangement(hyperplanes,K)
 966 
 967     def Shi(self,n,K=QQ):
 968         r"""
 969         The Shi arrangement is the set of `n(n-1)` hyperplanes: `\{ x_i - x_j = 0,1 : 1\leq i \leq j\leq n\}.`
 970 
 971         INPUT:
 972 
 973         - `n` -- integer
 974 
 975         OUTPUT:
 976 
 977         - HArrangement
 978 
 979         EXAMPLES::
 980 
 981             sage: hyperplane_arrangements.Shi(4)
 982             Hyperplane arrangement of 12 hyperplanes over Rational Field of dimension 4, rank 3.
 983         """
 984         hyperplanes = []
 985         for i in range(n):
 986             for j in range(i+1,n):
 987                 new = [0]*n
 988                 new[i] = 1
 989                 new[j] = -1
 990                 for k in [0,1]:
 991                     h = deepcopy(new)
 992                     h.append(k)
 993                     hyperplanes.append(h)
 994         return HArrangement(hyperplanes,K)
 995 
 996     def Catalan(self,n,K=QQ):
 997         r"""
 998         The Catalan arrangement is the set of `3n(n-1)/2` hyperplanes: `\{ x_i - x_j = -1,0,1 : 1\leq i \leq j\leq n\\}.`
 999 
1000         INPUT:
1001 
1002         - `n` -- integer
1003 
1004         OUTPUT:
1005 
1006         - HArrangement
1007 
1008         EXAMPLES::
1009 
1010             sage: hyperplane_arrangements.Catalan(5)
1011             Hyperplane arrangement of 30 hyperplanes over Rational Field of dimension 5, rank 4.
1012         """
1013         hyperplanes = []
1014         for i in range(n):
1015             for j in range(i+1,n):
1016                 new = [0]*n
1017                 new[i] = 1
1018                 new[j] = -1
1019                 for k in [-1,0,1]:
1020                     h = deepcopy(new)
1021                     h.append(k)
1022                     hyperplanes.append(h)
1023         Cn = HArrangement(hyperplanes,K)
1024         Cn._characteristic_polynomial = expand(x*prod([x-n-i for i in range(1,n)]))
1025         return Cn
1026 
1027     def coordinate(self,n,K=QQ):
1028         r"""
1029         The coordinate hyperplane arrangement is the central hyperplane arrangement consisting of the coordinate hyperplanes `x_i=0`.
1030 
1031         INPUT:
1032 
1033         - `n` -- integer
1034 
1035         OUTPUT:
1036 
1037         - HArrangement
1038 
1039         EXAMPLES::
1040 
1041             sage: hyperplane_arrangements.coordinate(5)
1042             Hyperplane arrangement of 5 hyperplanes over Rational Field of dimension 5, rank 5.
1043         """
1044         hyperplanes = []
1045         for i in range(n):
1046             new = [0]*(n+1)
1047             new[i]=1
1048             hyperplanes.append(new)
1049         return HArrangement(hyperplanes,K)
1050 
1051     def graphical(self,G,K=QQ):
1052         hyperplanes = []
1053         n = G.num_verts()
1054         for e in G.edges():
1055             i = G.vertices().index(e[0])
1056             j = G.vertices().index(e[1])
1057             new = [0]*(n+1)
1058             new[i] = 1
1059             new[j] = -1
1060             hyperplanes.append(new)
1061         A = HArrangement(hyperplanes,K)
1062         A._characteristic_polynomial = G.chromatic_polynomial()
1063         return A
1064 
1065     def braid(self,n,K=QQ):
1066         r"""
1067         The braid arrangement is the set of `n(n-1)/2` hyperplanes: `\{ x_i - x_j = 0 : 1\leq i \leq j\leq n\\}.`
1068 
1069         INPUT:
1070 
1071         - `n` integer
1072 
1073         OUTPUT:
1074 
1075         - HArrangement
1076 
1077         EXAMPLES::
1078 
1079             sage: hyperplane_arrangements.braid(4)
1080             Hyperplane arrangement of 6 hyperplanes over Rational Field of dimension 4, rank 3.
1081         """
1082         A = self.graphical(graphs.CompleteGraph(n),K)
1083         A._characteristic_polynomial = expand(prod(x-i for i in range(n)))
1084         return A
1085 
1086     def G_Shi(self,G,K=QQ):
1087         hyperplanes = []
1088         n = G.num_verts()
1089         for e in G.edges():
1090             i = G.vertices().index(e[0])
1091             j = G.vertices().index(e[1])
1092             new = [0]*(n+1)
1093             new[i] = 1
1094             new[j] = -1
1095             hyperplanes.append(new)
1096             new = deepcopy(new)
1097             new[-1]=1
1098             hyperplanes.append(new)
1099         return HArrangement(hyperplanes,K)
1100 
1101     def G_semiorder(self,G,K=QQ):
1102         hyperplanes = []
1103         n = G.num_verts()
1104         for e in G.edges():
1105             i = G.vertices().index(e[0])
1106             j = G.vertices().index(e[1])
1107             new = [0]*(n+1)
1108             new[i] = 1
1109             new[j] = -1
1110             new[-1] = -1
1111             hyperplanes.append(new)
1112             new = deepcopy(new)
1113             new[-1]=1
1114             hyperplanes.append(new)
1115         return HArrangement(hyperplanes,K)
1116 
1117     def G_arrangement(self,G,A=None,K=QQ):
1118         r"""
1119         A - None (semiorder), 'generic'
1120         """
1121         n = G.num_verts()
1122         if A==None:  #default to G-semiorder arrangement
1123             A = matrix(K,n,lambda i,j: 1)
1124         elif A=='generic':
1125             A = random_matrix(ZZ,n,x=10000)
1126             A = matrix(K,A)
1127         hyperplanes = []
1128         for e in G.edges():
1129             i = G.vertices().index(e[0])
1130             j = G.vertices().index(e[1])
1131             new = [0]*(n+1)
1132             new[i] = 1
1133             new[j] = -1
1134             new[-1] = A[i,j]
1135             hyperplanes.append(new)
1136             new = [0]*(n+1)
1137             new[j] = 1
1138             new[i] = -1
1139             new[-1] = A[j,i]
1140             hyperplanes.append(new)
1141         return HArrangement(hyperplanes,K)
1142 
1143 hyperplane_arrangements = HArrangementGenerators()
1144 harrangements = HArrangementGenerators()
1145 
1146 # idea: store information about some of these hyperplanes
1147 
1148 def repr_point(P,scale=1):
1149     c = P.center()
1150     if P.n_rays()==0:
1151         return c
1152     else:
1153         r = sum([vector(QQ,i) for i in P.rays()])
1154         return c + scale*r/r.norm()