Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
""" Diamond cutting implementation
AUTHORS:
- Jan Poeschko (2012-07-02): initial version """
#***************************************************************************** # Copyright (C) 2012 Jan Poeschko <jan@poeschko.com> # # Distributed under the terms of the GNU General Public License (GPL) # as published by the Free Software Foundation; either version 2 of # the License, or (at your option) any later version. # http://www.gnu.org/licenses/ #*****************************************************************************
""" Return the inequality for points on the same side as the origin with respect to the plane through ``v`` normal to ``v``.
EXAMPLES::
sage: from sage.modules.diamond_cutting import plane_inequality sage: ieq = plane_inequality([1, -1]); ieq [2, -1, 1] sage: ieq[0] + vector(ieq[1:]) * vector([1, -1]) 0 """
r""" Compute the upper-triangular part of the Cholesky/Jacobi decomposition of the symmetric matrix ``M``.
Let `M` be a symmetric `n \times n`-matrix over a field `F`. Let `m_{i,j}` denote the `(i,j)`-th entry of `M` for any `1 \leq i \leq n` and `1 \leq j \leq n`. Then, the upper-triangular part computed by this method is the upper-triangular `n \times n`-matrix `Q` whose `(i,j)`-th entry `q_{i,j}` satisfies
.. MATH::
q_{i,j} = \begin{cases} \frac{1}{q_{i,i}} \left( m_{i,j} - \sum_{r<i} q_{r,r} q_{r,i} q_{r,j} \right) & i < j, \\ a_{i,j} - \sum_{r<i} q_{r,r} q_{r,i}^2 & i = j, \\ 0 & i > j, \end{cases}
for all `1 \leq i \leq n` and `1 \leq j \leq n`. (These equalities determine the entries of `Q` uniquely by recursion.) This matrix `Q` is defined for all `M` in a certain Zariski-dense open subset of the set of all `n \times n`-matrices.
EXAMPLES::
sage: from sage.modules.diamond_cutting import jacobi sage: jacobi(identity_matrix(3) * 4) [4 0 0] [0 4 0] [0 0 4]
sage: def testall(M): ....: Q = jacobi(M) ....: for j in range(3): ....: for i in range(j): ....: if Q[i,j] * Q[i,i] != M[i,j] - sum(Q[r,i] * Q[r,j] * Q[r,r] for r in range(i)): ....: return False ....: for i in range(3): ....: if Q[i,i] != M[i,i] - sum(Q[r,i] ** 2 * Q[r,r] for r in range(i)): ....: return False ....: for j in range(i): ....: if Q[i,j] != 0: ....: return False ....: return True
sage: M = Matrix(QQ, [[8,1,5], [1,6,0], [5,0,3]]) sage: Q = jacobi(M); Q [ 8 1/8 5/8] [ 0 47/8 -5/47] [ 0 0 -9/47] sage: testall(M) True
sage: M = Matrix(QQ, [[3,6,-1,7],[6,9,8,5],[-1,8,2,4],[7,5,4,0]]) sage: testall(M) True """ raise ValueError("the matrix must be square")
r""" Perform diamond cutting on polyhedron ``V`` with basis matrix ``GM`` and radius ``C``.
INPUT:
- ``V`` -- polyhedron to cut from
- ``GM`` -- half of the basis matrix of the lattice
- ``C`` -- radius to use in cutting algorithm
- ``verbose`` -- (default: ``False``) whether to print debug information
OUTPUT:
A :class:`Polyhedron` instance.
EXAMPLES::
sage: from sage.modules.diamond_cutting import diamond_cut sage: V = Polyhedron([[0], [2]]) sage: GM = matrix([2]) sage: V = diamond_cut(V, GM, 4) sage: V.vertices() (A vertex at (2), A vertex at (0)) """ # coerce to floats print("Cut\n{}\nwith radius {}".format(GM, C))
raise ValueError("the matrix must be square")
# calculate the Gram matrix print( "q:\n{}".format(q.n()) ) # apply Cholesky/Jacobi decomposition print( "q:\n{}".format(q.n()) )
print("Dimension: {}".format(i)) print("Z: {}".format(Z)) print("L: {}".format(L))
print("x: {}".format(x)) else:
print("\n%d) Cut using normal vector %s" % (cut_count, hv)) #cut = Polyhedron(ieqs=[plane_inequality(hv)]) #V = V.intersection(cut)
print("Final cut")
print("End")
""" Calculate the Voronoi cell of the lattice defined by basis
INPUT:
- ``basis`` -- embedded basis matrix of the lattice
- ``radius`` -- radius of basis vectors to consider
- ``verbose`` -- whether to print debug information
OUTPUT:
A :class:`Polyhedron` instance.
EXAMPLES::
sage: from sage.modules.diamond_cutting import calculate_voronoi_cell sage: V = calculate_voronoi_cell(matrix([[1, 0], [0, 1]])) sage: V.volume() 1 """
# introduce "artificial" basis points (representing infinity) # LLL-reduce to get quadratic matrix raise ValueError("invalid matrix")
# twice the length of longest vertex in Q is a safe choice
# remove inequalities introduced by artificial basis points not V._is_zero(v.A() * (-w) / 2 - v.b())) for w in additional_vectors)]
|