def __init__(self, projection_direction, height = 1.1):
        """
        Initializes the projection.

        EXAMPLES::

            sage: from sage.geometry.polyhedron.plot import ProjectionFuncSchlegel
            sage: proj = ProjectionFuncSchlegel([2,2,2])
            sage: proj.__init__([2,2,2])
            sage: proj(vector([1.1,1.1,1.11]))[0]
            0.0302...
            sage: TestSuite(proj).run(skip='_test_pickling')
        """
        self.projection_dir = vector(RDF, projection_direction)
        if norm(self.projection_dir).is_zero():
            raise ValueError, "projection direction must be a non-zero vector."
        self.dim = self.projection_dir.degree()
        spcenter = height * self.projection_dir/norm(self.projection_dir)
        self.height = height
        v = vector(RDF, [0.0]*(self.dim-1) + [self.height]) - spcenter
        polediff = matrix(RDF,v).transpose()
        denom = (polediff.transpose()*polediff)[0][0]
        if denom.is_zero():
            self.house = identity_matrix(RDF,self.dim)
        else:
            self.house = identity_matrix(RDF,self.dim) \
            - 2*polediff*polediff.transpose()/denom #Householder reflector

    def mutate(self, k, mutating_F=True):
        r"""
        mutate seed
        bla bla ba
        """
        n = self.parent().rk

        if k not in xrange(n):
            raise ValueError('Cannot mutate in direction ' + str(k) + '.')

        # store mutation path
        if self._path != [] and self._path[-1] == k:
            self._path.pop()
        else:
            self._path.append(k)

        # find sign of k-th c-vector
        # Will this be used enough that it should be a built-in function?
        if any(x > 0 for x in self._C.column(k)):
            eps = +1
        else:
            eps = -1

        # store the g-vector to be mutated in case we are mutating F-polynomials also
        old_g_vector = self.g_vector(k)

        # G-matrix
        J = identity_matrix(n)
        for j in xrange(n):
            J[j,k] += max(0, -eps*self._B[j,k])
        J[k,k] = -1
        self._G = self._G*J

        # g-vector path list
        g_vector = self.g_vector(k)
        if g_vector not in self.parent().g_vectors_so_far():
            self.parent()._path_dict[g_vector] = copy(self._path)
            # F-polynomials
            if mutating_F:
                self.parent()._F_poly_dict[g_vector] = self._mutated_F(k, old_g_vector)

        # C-matrix
        J = identity_matrix(n)
        for j in xrange(n):
            J[k,j] += max(0, eps*self._B[k,j])
        J[k,k] = -1
        self._C = self._C*J

        # B-matrix
        self._B.mutate(k)

        # exchange relation
        if self.parent()._store_exchange_relations:
            ex_pair = frozenset([g_vector,old_g_vector])
            if ex_pair not in self.parent()._exchange_relations:
                coefficient = self.coefficient(k).lift()
                variables = zip(self.g_vectors(), self.b_matrix().column(k))
                Mp = [ (g,p) for (g,p) in variables if p > 0 ]
                Mm = [ (g,-p) for (g,p) in variables if p < 0 ]
                self.parent()._exchange_relations[ex_pair] = ( (Mp,coefficient.numerator()), (Mm,coefficient.denominator()) )

    def mutate(self, sequence, inplace=True, mutate_F=True):
        if not isinstance(mutate_F, bool):
            raise ValueError('mutate_F must be boolean.')

        if not isinstance(inplace, bool):
            raise ValueError('inplace must be boolean.')
        if inplace:
            seed = self.current_seed
        else:
            seed = copy(self.current_seed)
            
        if sequence in xrange(seed._n):
            seq = iter([sequence])
        else:
            seq = iter(sequence)
            
        for k in seq:
            if k not in xrange(seed._n):
                raise ValueError('Cannot mutate in direction' + str(k) + '.')

            # G-matrix
            J = identity_matrix(seed._n)
            if any(x > 0 for x in seed._C.column(k)):
                eps = +1
            else:
                eps = -1
            for j in xrange(seed._n):
                J[j,k] += max(0, -eps*seed._B[j,k])
            J[k,k] = -1
            seed._G = seed._G*J

            # F-polynomials
            if mutating_F:
                gvect = tuple(seed._G.column(k))
                if not self._F_dict.has_key(gvect):
                    self._F_dict.setdefault(gvect,self._mutated_F(k))
                seed._F[k] = self._F_dict[gvect]
            
            # C-matrix
            J = identity_matrix(seed._n)
            if any(x > 0 for x in seed._C.column(k)):
                eps = +1
            else:
                eps = -1
            for j in xrange(seed._n):
                J[k,j] += max(0, eps*seed._B[k,j])
            J[k,k] = -1
            seed._C = seed._C*J
            
            # B-matrix
            seed._B.mutate(k)
            
        if not inplace:
            return seed

    def _extend1(self):
        "Increases the matrix size by 1"

        #save current self
        self.iv0 = copy(self)
        self.iv0.abel_raw0 = lambda z: self.iv0.abel0poly(z-self.iv0.x0)

        N = self.N

        A = self.A
        #Carleman matrix without 0-th row:
        Ct = A + identity_matrix(self.R,N).submatrix(1,0,N-1,N-1)
        AI = self.AI
        #assert AI*self.A == identity_matrix(N-1)

        if isinstance(self.f,FormalPowerSeries):
            coeffs = [ self.f[n] for n in xrange(0,N) ]
        else:
            x = self.f.args()[0]
            coeffs = taylor(self.f.substitute({x:x+self.x0sym})-self.x0sym,x,0,N).polynomial(self.R)

        self.Ct = matrix(self.R,N,N)
        self.Ct.set_block(0,0,Ct)
        self.Ct[0,N-1] = coeffs[N-1]
        for m in range(1,N-1):
            self.Ct[m,N-1] = psmul_at(self.Ct[0],self.Ct[m-1],N-1)
        self.Ct[N-1] = psmul(self.Ct[0],self.Ct[N-2])

        print "C extended"
        self.A = self.Ct - identity_matrix(self.R,N+1).submatrix(1,0,N,N)

        av=self.A.column(N-1)[:N-1]
        ah=self.A.row(N-1)[:N-1]
        a_n=self.A[N-1,N-1]

        # print 'A:';print A
        # print 'A\''; print self.A
        # print 'ah:',ah
        # print 'av:',av
        # print 'a_n:',a_n

        AI0 = matrix(self.R,N,N)
        AI0.set_block(0,0,self.AI)

        horiz = matrix(self.R,1,N)
        horiz.set_block(0,0,(ah*AI).transpose().transpose())
        horiz[0,N-1] = -1

        vert = matrix(self.R,N,1)
        vert.set_block(0,0,(AI*av).transpose())
        vert[N-1,0] = -1

        self.N += 1
        self.AI = AI0 +  vert*horiz/(a_n-ah*AI*av)

    def mutate_initial(self, k):
        r"""
        Mutate ``self`` in direction `k` at the initial cluster.

        INPUT:
        - ``k`` -- integer in between 0 and ``self.rk``
        """
        n = self.rk

        if k not in xrange(n):
            raise ValueError('Cannot mutate in direction %s, please try a value between 0 and %s.'%(str(k),str(n-1)))

        #modify self._path_dict using Nakanishi-Zelevinsky (4.1) and self._F_poly_dict using CA-IV (6.21)
        new_path_dict = dict()
        new_F_dict = dict()
        new_path_dict[tuple(identity_matrix(n).column(k))] = []
        new_F_dict[tuple(identity_matrix(n).column(k))] = self._U(1)

        poly_ring = PolynomialRing(ZZ,'u')
        h_subs_tuple = tuple([poly_ring.gen(0)**(-1) if j==k else poly_ring.gen(0)**max(-self._B0[k][j],0) for j in xrange(n)])
        F_subs_tuple = tuple([self._U.gen(k)**(-1) if j==k else self._U.gen(j)*self._U.gen(k)**max(-self._B0[k][j],0)*(1+self._U.gen(k))**(self._B0[k][j]) for j in xrange(n)])

        for g_vect in self._path_dict:
            #compute new path
            path = self._path_dict[g_vect]
            if g_vect == tuple(identity_matrix(n).column(k)):
                new_path = [k]
            elif path != []:
                if path[0] != k:
                    new_path = [k] + path
                else:
                    new_path = path[1:]
            else:
                new_path = []

            #compute new g-vector
            new_g_vect = vector(g_vect) - 2*g_vect[k]*identity_matrix(n).column(k)
            for i in xrange(n):
                new_g_vect += max(sign(g_vect[k])*self._B0[i,k],0)*g_vect[k]*identity_matrix(n).column(i)
            new_path_dict[tuple(new_g_vect)] = new_path

            #compute new F-polynomial
            h = 0
            trop = tropical_evaluation(self._F_poly_dict[g_vect](h_subs_tuple))
            if trop != 1:
                h = trop.denominator().exponents()[0]-trop.numerator().exponents()[0]
            new_F_dict[tuple(new_g_vect)] = self._F_poly_dict[g_vect](F_subs_tuple)*self._U.gen(k)**h*(self._U.gen(k)+1)**g_vect[k]

        self._path_dict = new_path_dict
        self._F_poly_dict = new_F_dict

        self._B0.mutate(k)

    def burau_matrix(self, var='t'):
        """
        Return the Burau matrix of the braid.

        INPUT:

        - ``var`` -- string (default: ``'t'``). The name of the
          variable in the entries of the matrix.

        OUTPUT:

        The Burau matrix of the braid. It is a matrix whose entries
        are Laurent polynomials in the variable ``var``.

        EXAMPLES::

            sage: B = BraidGroup(4)
            sage: B.inject_variables()
            Defining s0, s1, s2
            sage: b=s0*s1/s2/s1
            sage: b.burau_matrix()
            [     -t + 1           0    -t^2 + t         t^2]
            [          1           0           0           0]
            [          0           0           1           0]
            [          0        t^-2 t^-1 - t^-2    1 - t^-1]
            sage: s2.burau_matrix('x')
            [     1      0      0      0]
            [     0      1      0      0]
            [     0      0 -x + 1      x]
            [     0      0      1      0]

        REFERENCES:

            http://en.wikipedia.org/wiki/Burau_representation
        """
        R = LaurentPolynomialRing(IntegerRing(), var)
        t = R.gen()
        M = identity_matrix(R, self.strands())
        for i in self.Tietze():
            A = identity_matrix(R, self.strands())
            if i>0:
                A[i-1, i-1] = 1-t
                A[i, i] = 0
                A[i, i-1] = 1
                A[i-1, i] = t
            if i<0:
                A[-1-i, -1-i] = 0
                A[-i, -i] = 1-t**(-1)
                A[-1-i, -i] = 1
                A[-i, -1-i] = t**(-1)
            M=M*A
        return M

    def invariant_form(self):
        """
        Return the hermitian form preserved by the unitary
        group.

        OUTPUT:

        A square matrix describing the bilinear form

        EXAMPLES::

            sage: SU4 = SU(4,QQ)
            sage: SU4.invariant_form()
            [1 0 0 0]
            [0 1 0 0]
            [0 0 1 0]
            [0 0 0 1]
        """
        if self._invariant_form is not None:
            return self._invariant_form

        from sage.matrix.constructor import identity_matrix
        m = identity_matrix(self.base_ring(), self.degree())
        m.set_immutable()
        return m

    def calculate_generators(self):
        """
        If generators haven't already been computed, calculate generators
        for this homspace. If they have been computed, do nothing.

        EXAMPLES::

            sage: E = End(J0(11))
            sage: E.calculate_generators()
        """

        if self._gens is not None:
            return

        if (self.domain() == self.codomain()) and (self.domain().dimension() == 1):
            self._gens = tuple([ identity_matrix(ZZ,2) ])
            return

        phi = self.domain()._isogeny_to_product_of_powers()
        psi = self.codomain()._isogeny_to_product_of_powers()

        H_simple = phi.codomain().Hom(psi.codomain())
        im_gens = H_simple._calculate_product_gens()

        M = phi.matrix()
        Mt = psi.complementary_isogeny().matrix()

        R = ZZ**(4*self.domain().dimension()*self.codomain().dimension())
        gens = R.submodule([ (M*self._get_matrix(g)*Mt).list()
                             for g in im_gens ]).saturation().basis()
        self._gens = tuple([ self._get_matrix(g) for g in gens ])

    def ARP(self):
        from sage.matrix.constructor import identity_matrix
        A1 = matrix(3, [1,1,1, 0,1,0, 0,0,1])
        A2 = matrix(3, [1,0,0, 1,1,1, 0,0,1])
        A3 = matrix(3, [1,0,0, 0,1,0, 1,1,1])
        P12 = matrix(3, [1,0,1, 1,1,1, 0,0,1])
        P13 = matrix(3, [1,1,0, 0,1,0, 1,1,1])
        P23 = matrix(3, [1,0,0, 1,1,0, 1,1,1])
        P21 = matrix(3, [1,1,1, 0,1,1, 0,0,1])
        P31 = matrix(3, [1,1,1, 0,1,0, 0,1,1])
        P32 = matrix(3, [1,0,0, 1,1,1, 1,0,1])
        gens = (A1, A2, A3, P23, P32, P13, P31, P12, P21)
        alphabet = [1, 2, 3, 123, 132, 213, 231, 312, 321]
        gens = dict(zip(alphabet, gens))

        cone = {}
        cone[123] = H23 = matrix(3, [1,0,1, 0,1,0, 0,0,1])
        cone[132] = H32 = matrix(3, [1,1,0, 0,1,0, 0,0,1])
        cone[213] = H13 = matrix(3, [1,0,0, 0,1,1, 0,0,1])
        cone[231] = H31 = matrix(3, [1,0,0, 1,1,0, 0,0,1])
        cone[312] = H12 = matrix(3, [1,0,0, 0,1,0, 0,1,1])
        cone[321] = H21 = matrix(3, [1,0,0, 0,1,0, 1,0,1])
        cone[1] = cone[2] = cone[3] = identity_matrix(3) 

        from .language import languages
        return MatrixCocycle(gens, cone, language=languages.ARP())

 def initial_seed(self):
     r"""
     Return the initial seed
     """
     n = self.rk
     I = identity_matrix(n)
     return ClusterAlgebraSeed(self._B0, I, I, self)

    def initial_pair(self):
        """
        Return an initial double description pair.

        Picks an initial set of rays by selecting a basis. This is
        probably the most efficient way to select the initial set.

        INPUT:

        - ``pair_class`` -- subclass of
          :class:`DoubleDescriptionPair`.

        OUTPUT:

        A pair consisting of a :class:`DoubleDescriptionPair` instance
        and the tuple of remaining unused inequalities.

        EXAMPLES::

            sage: A = matrix([(-1, 1), (-1, 2), (1/2, -1/2), (1/2, 2)])
            sage: from sage.geometry.polyhedron.double_description import Problem
            sage: DD, remaining = Problem(A).initial_pair()
            sage: DD.verify()
            sage: remaining
            [(1/2, -1/2), (1/2, 2)]
        """
        pivot_rows = self.A_matrix().pivot_rows()
        A0 = [self.A()[pivot] for pivot in pivot_rows]
        Ac = [self.A()[i] for i in range(len(self.A())) if i not in pivot_rows]
        from sage.matrix.constructor import identity_matrix, matrix
        I = identity_matrix(self.base_ring(), self.dim())
        R = matrix(self.base_ring(), A0).solve_right(I)
        return self.pair_class(self, A0, R.columns()), list(Ac)

def calculate_voronoi_cell(basis, radius=None, verbose=False):
    """
    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
    """

    dim = basis.dimensions()
    artificial_length = None
    if dim[0] < dim[1]:
        # introduce "artificial" basis points (representing infinity)
        artificial_length = max(abs(v) for v in basis).ceil() * 2
        additional_vectors = identity_matrix(dim[1]) * artificial_length
        basis = basis.stack(additional_vectors)
        # LLL-reduce to get quadratic matrix
        basis = basis.LLL()
        basis = matrix([v for v in basis if v])
        dim = basis.dimensions()
    if dim[0] != dim[1]:
        raise ValueError("invalid matrix")
    basis = basis / 2

    ieqs = []
    for v in basis:
        ieqs.append(plane_inequality(v))
        ieqs.append(plane_inequality(-v))
    Q = Polyhedron(ieqs=ieqs)

    # twice the length of longest vertex in Q is a safe choice
    if radius is None:
        radius = 2 * max(abs(v) ** 2 for v in basis)

    V = diamond_cut(Q, basis, radius, verbose=verbose)

    if artificial_length is not None:
        # remove inequalities introduced by artificial basis points
        H = V.Hrepresentation()
        H = [v for v in H if all(not V._is_zero(v.A() * w / 2 - v.b() and
             not V._is_zero(v.A() * (-w) / 2 - v.b())) for w in additional_vectors)]
        V = Polyhedron(ieqs=H)

    return V

    def matrix_power(self,t):
        eigenvalues = self.eigenvalues
        iprec = self.iprec
        CM = self.CM
        
        ev = [ e.n(iprec) for e in eigenvalues]
        n = len(ev)
    
        Char = [CM - ev[k] * identity_matrix(n) for k in range(n)]
    
        #product till k-1
        prodwo = n * [0]
        prod = identity_matrix(n)

	#if we were to start here with prod = IdentityMatrix
	#prodwo[k]/sprodwo[k] would be the component projector of ev[k]
        #component projector of ev[k] is a matrix Z such that
        #CM * Z = ev[k] * Z and Z*Z=Z
        #then f(CM)=sum_k f(ev[k])*Z[k]
        #as we are only interested in the first line we can start 
        #left with (0,1,...) instead of the identity matrix
        
        for k in range(n):
            prodwo[k] = prod
            for i in range(k+1,n):
                prodwo[k] = prodwo[k] * Char[i]
    
            if k == n:
                break
            prod = prod * Char[k]
    
        sprodwo = n * [0]
        for k in range(n):
            if k == 0:
                sprodwo[k] = ev[k] - ev[1]
                start = 2
            else:
                sprodwo[k] = ev[k] - ev[0]
                start = 1
            
            for i in range(start,n):
                if not i == k:
                    sprodwo[k] = sprodwo[k] * (ev[k]-ev[i])
    
        return sum([ev[k]**t/sprodwo[k]*prodwo[k] for k in range(n)])

 def reduce(self, t) :
     ## We compute the rational diagonal form of t. Whenever a zero entry occures we
     ## find a primitive isotropic vector and apply a base change, such that t finally
     ## has the form diag(0,0,...,P) where P is positive definite. P will then be a
     ## LLL reduced.
     
     ot = t
     t = matrix(QQ, t)
     n = t.nrows()
     u = identity_matrix(QQ, n)
     
     for i in xrange(n) :
         if t[i,i] < 0 :
             return None
         elif t[i,i] == 0 :
             ## get the isotropic vector
             v = u.column(i)
             v = v.denominator() * v
             v = v / gcd(list(v))
             u = self._gln_lift(v, 0)
             
             t = u.transpose() * ot * u
             ts = self.reduce(t.submatrix(1,1,n-1,n-1))
             if ts is None :
                 return None
             t.set_block(1, 1, ts[0])
             
             t.set_immutable()
             
             return (t,1)
         else :
             for j in xrange(i + 1, n) :
                 us = identity_matrix(QQ, n, n)
                 us[i,j] = -t[i,j]/t[i,i]
                 u = u * us
                  
                 t.add_multiple_of_row(j, i, -t[j,i]/t[i,i])
                 t.add_multiple_of_column(j, i, -t[i,j]/t[i,i])
     
     u = ot.LLL_gram()
     t = u.transpose() * ot * u
     t.set_immutable()
     
     return (t, 1)

    def _rank(self, K) :
        
        if K is QQ or K in NumberFields() :
            return len(_jacobi_forms_by_taylor_expansion_coords(self.__index, self.__weight, 0))

            ## This is the formula used by Poor and Yuen in Paramodular cusp forms
            if self.__weight == 2 :
                delta = len(self.__index.divisors()) // 2 - 1
            else :
                delta = 0
                
            return sum( ModularForms(1, self.__weight + 2 * j).dimension() + j**2 // (4 * self.__index)
                        for j in xrange(self.__index + 1) ) \
                   + delta
            

            ## This is the formula given by Skoruppa in 
            ## Jacobi forms of critical weight and Weil representations
            ##FIXME: There is some mistake here
            if self.__weight % 2 != 0 :
                ## Otherwise the space X(i**(n - 2 k)) is different
                ## See: Skoruppa, Jacobi forms of critical weight and Weil representations
                raise NotImplementedError
            
            m = self.__index
            K = CyclotomicField(24 * m, 'zeta')
            zeta = K.gen(0)
            
            quadform = lambda x : 6 * x**2
            bilinform = lambda x,y : quadform(x + y) - quadform(x) - quadform(y)
            
            T = diagonal_matrix([zeta**quadform(i) for i in xrange(2*m)])
            S =   sum(zeta**(-quadform(x)) for x in xrange(2 * m)) / (2 * m) \
                * matrix([[zeta**(-bilinform(j,i)) for j in xrange(2*m)] for i in xrange(2*m)])
            subspace_matrix_1 = matrix( [ [1 if j == i or j == 2*m - i else 0 for j in xrange(m + 1) ]
                                        for i in xrange(2*m)] )
            subspace_matrix_2 = zero_matrix(ZZ, m + 1, 2*m)
            subspace_matrix_2.set_block(0,0,identity_matrix(m+1))
            
            T = subspace_matrix_2 * T * subspace_matrix_1
            S = subspace_matrix_2 * S * subspace_matrix_1
            
            sqrt3 = (zeta**(4*m) - zeta**(-4*m)) * zeta**(-6*m) 
            rank =   (self.__weight - 1/2 - 1) / 2 * (m + 1) \
                   + 1/8 * (   zeta**(3*m * (2*self.__weight - 1)) * S.trace()
                             + zeta**(3*m * (1 - 2*self.__weight)) * S.trace().conjugate() ) \
                   + 2/(3*sqrt3) * (   zeta**(4 * m * self.__weight) * (S*T).trace()
                                     + zeta**(-4 * m * self.__weight) * (S*T).trace().conjugate() ) \
                   - sum((j**2 % (m+1))/(m+1) -1/2 for j in range(0,m+1))
            
            if self.__weight > 5 / 2 :
                return rank
            else :
                raise NotImplementedError
            
        raise NotImplementedError

def homogeneous_components(self):
    deg_matrix = block_matrix([[identity_matrix(self.parent().rk),-self.parent()._B0]])
    components = dict()
    x = self.lift()
    monomials = x.monomials()
    for m in monomials:
        g_vect = tuple(deg_matrix*vector(m.exponents()[0]))
        if g_vect in components:
            components[g_vect] += self.parent().retract(x.monomial_coefficient(m)*m)
        else:
            components[g_vect] = self.parent().retract(x.monomial_coefficient(m)*m)
    return components

def lift(A, N):
    r"""
    Lift a matrix A from SL_m(Z/NZ) to SL_m(Z).

    Follows Shimura, Lemma 1.38, p21.

    """
    assert A.is_square()
    assert A.determinant() != 0
    m = A.nrows()
    if m == 1:
        return identity_matrix(1)

    D, U, V = A.smith_form()
    if U.determinant() == -1 :
        U = matrix(2,2,[-1,0,0,1])* U
    if V.determinant() == -1 :
        V = V *matrix(2,2,[-1,0,0,1])
    D = U*A*V
    assert U.determinant() == 1
    assert V.determinant() == 1
    a = [ D[i, i] for i in range(m) ]
    b = prod(a[1:])
    W = identity_matrix(m)
    W[0, 0] = b
    W[1, 0] = b-1
    W[0, 1] = 1
    X = identity_matrix(m)
    X[0, 1] = -a[1]
    Ap = D.parent()(D)
    Ap[0, 0] = 1
    Ap[1, 0] = 1-a[0]
    Ap[1, 1] *= a[0]
    assert (W*U*A*V*X).change_ring(Zmod(N)) == Ap.change_ring(Zmod(N))
    Cp = diagonal_matrix(a[1:])
    Cp[0, 0] *= a[0]
    C = lift(Cp, N)
    Cpp = block_diagonal_matrix(identity_matrix(1), C)
    Cpp[1, 0] = 1-a[0]
    return (~U * ~W * Cpp * ~X * ~V).change_ring(ZZ)

    def cone_points_iter(self):
        """
        Iterate over the open torus orbits and yield distinct points.
        
        OUTPUT:

        For each open torus orbit (cone): A triple consisting of the
        cone, the nonzero homogeneous coordinates in that orbit (list
        of integers), and the nonzero log coordinates of distinct
        points as a cokernel.

        EXAMPLES::

            sage: fan = NormalFan(ReflexivePolytope(2, 0))
            sage: X = ToricVariety(fan, base_ring=GF(7))
            sage: point_set = X.point_set()
            sage: ffe = point_set._finite_field_enumerator()
            sage: cpi = ffe.cone_points_iter()
            sage: cone, nonzero_points, cokernel = list(cpi)[5]
            sage: cone
            1-d cone of Rational polyhedral fan in 2-d lattice N
            sage: cone.ambient_ray_indices()
            (2,)
            sage: nonzero_points
            [0, 1]
            sage: cokernel
            Finitely generated module V/W over Integer Ring with invariants (2)
            sage: list(cokernel)
            [(0), (1)]
            sage: [p.lift() for p in cokernel]
            [(0, 0), (0, 1)]
        """
        from sage.matrix.constructor import matrix, block_matrix, identity_matrix
        from sage.rings.all import ZZ
        nrays = len(self.rays())
        N = self.multiplicative_group_order()
        # Want cokernel of the log rescalings in (ZZ/N)^(#rays). But
        # ZZ/N is not a integral domain. Instead: work over ZZ
        log_generators = self.rescaling_log_generators()
        log_relations = block_matrix(2, 1, [
            matrix(ZZ, len(log_generators), nrays, log_generators), 
            N * identity_matrix(ZZ, nrays)])
        for cone in self.cone_iter():
            nrays = self.fan().nrays() + len(self.fan().virtual_rays())
            nonzero_coordinates = [i for i in range(nrays)
                                   if i not in cone.ambient_ray_indices()]
            log_relations_nonzero = log_relations.matrix_from_columns(nonzero_coordinates)
            image = log_relations_nonzero.image()
            cokernel = image.ambient_module().quotient(image)
            yield cone, nonzero_coordinates, cokernel

def permutation_simplicial_action(r,u,n,w):
    r"""
    From a word in 'l','r' return the simplicial action on homology as
    well as the obtained origami.

    INPUT:

    - ``r`` and ``u`` - permutations

    - ``n`` - the degree of the permutations

    - ``w`` - a string in 'l' and 'r' or a list of 2-tuples which are made of
      the letter 'l' or 'r' and a positive integer

    """
    if w is None:
        w = []
    elif isinstance(w,list):
        w = flatten_word(w)

    res = identity_matrix(2*n)

    for letter in reversed(w):
        if letter == 'l':
            u = u*~r
            m = identity_matrix(2*n)
            m[:n,n:] = u.matrix()
            res = m * res
        elif letter == 'r':
            r = r*~u
            m = identity_matrix(2*n)
            m[n:,:n] = r.matrix()
            res = m * res
        else:
            raise ValueError, "does not understand the letter %s" %str(letter)

    return r,u,res

def construct_from_dim_degree(dim, max_degree, base_ring, check):
    """
    Construct a filtered vector space.

    INPUT:

    - ``dim`` -- integer. The dimension.

    - ``max_degree`` -- integer or infinity. The maximal degree where
      the vector subspace of the filtration is still the entire space.

    EXAMPLES::

        sage: V = FilteredVectorSpace(2, 5);  V
        QQ^2 >= 0
        sage: V.get_degree(5)
        Vector space of degree 2 and dimension 2 over Rational Field
        Basis matrix:
        [1 0]
        [0 1]
        sage: V.get_degree(6)
        Vector space of degree 2 and dimension 0 over Rational Field
        Basis matrix:
        []

        sage: FilteredVectorSpace(2, oo)
        QQ^2
        sage: FilteredVectorSpace(2, -oo)
        0 in QQ^2

    TESTS::

        sage: from sage.modules.filtered_vector_space import construct_from_dim_degree
        sage: V = construct_from_dim_degree(2, 5, QQ, True);  V
        QQ^2 >= 0
    """
    if dim not in ZZ:
        raise ValueError('dimension must be an integer')
    dim = ZZ(dim)
    from sage.matrix.constructor import identity_matrix
    generators = identity_matrix(base_ring, dim).columns()
    filtration = dict()
    if max_degree is None:
        max_degree = infinity
    filtration[normalize_degree(max_degree)] = range(dim)
    return construct_from_generators_indices(generators, filtration, base_ring, check)