from sage.matrix.constructor import matrix
from sage.misc.cachefunc import cached_method
from sage.rings.all import QQ

from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra


class FiniteDimensionalEJAElementSubalgebra(FiniteDimensionalEJASubalgebra):
    def __init__(self, elt, orthonormalize=True, **kwargs):
        superalgebra = elt.parent()

        powers = tuple( elt**k for k in range(superalgebra.dimension()) )
        power_vectors = ( p.to_vector() for p in powers )
        P = matrix(superalgebra.base_ring(), power_vectors)

        if orthonormalize:
            basis = powers # let god sort 'em out
        else:
            # Echelonize the matrix ourselves, because otherwise the
            # call to P.pivot_rows() below can choose a non-optimal
            # row-reduction algorithm. In particular, scaling can
            # help over AA because it avoids the RecursionError that
            # gets thrown when we have to look too hard for a root.
            #
            # Beware: QQ supports an entirely different set of "algorithm"
            # keywords than do AA and RR.
            algo = None
            if superalgebra.base_ring() is not QQ:
                algo = "scaled_partial_pivoting"
                P.echelonize(algorithm=algo)

                # In this case, we just need to figure out which elements
                # of the "powers" list are redundant... First compute the
                # vector subspace spanned by the powers of the given
                # element.

                # Figure out which powers form a linearly-independent set.
                ind_rows = P.pivot_rows()

                # Pick those out of the list of all powers.
                basis = tuple(map(powers.__getitem__, ind_rows))


        super().__init__(superalgebra,
                         basis,
                         associative=True,
                         **kwargs)

        # The rank is the highest possible degree of a minimal
        # polynomial, and is bounded above by the dimension. We know
        # in this case that there's an element whose minimal
        # polynomial has the same degree as the space's dimension
        # (remember how we constructed the space?), so that must be
        # its rank too.
        self.rank.set_cache(self.dimension())


    @cached_method
    def one(self):
        """
        Return the multiplicative identity element of this algebra.

        The superclass method computes the identity element, which is
        beyond overkill in this case: the superalgebra identity
        restricted to this algebra is its identity. Note that we can't
        count on the first basis element being the identity -- it might
        have been scaled if we orthonormalized the basis.

        SETUP::

            sage: from mjo.eja.eja_algebra import (HadamardEJA,
            ....:                                  random_eja)

        EXAMPLES::

            sage: J = HadamardEJA(5)
            sage: J.one()
            e0 + e1 + e2 + e3 + e4
            sage: x = sum(J.gens())
            sage: A = x.subalgebra_generated_by()
            sage: A.one()
            f0
            sage: A.one().superalgebra_element()
            e0 + e1 + e2 + e3 + e4

        TESTS:

        The identity element acts like the identity over the rationals::

            sage: set_random_seed()
            sage: x = random_eja(field=QQ,orthonormalize=False).random_element()
            sage: A = x.subalgebra_generated_by()
            sage: x = A.random_element()
            sage: A.one()*x == x and x*A.one() == x
            True

        The identity element acts like the identity over the algebraic
        reals with an orthonormal basis::

            sage: set_random_seed()
            sage: x = random_eja().random_element()
            sage: A = x.subalgebra_generated_by(orthonormalize_basis=True)
            sage: x = A.random_element()
            sage: A.one()*x == x and x*A.one() == x
            True

        The matrix of the unit element's operator is the identity over
        the rationals::

            sage: set_random_seed()
            sage: x = random_eja(field=QQ,orthonormalize=False).random_element()
            sage: A = x.subalgebra_generated_by()
            sage: actual = A.one().operator().matrix()
            sage: expected = matrix.identity(A.base_ring(), A.dimension())
            sage: actual == expected
            True

        The matrix of the unit element's operator is the identity over
        the algebraic reals with an orthonormal basis::

            sage: set_random_seed()
            sage: x = random_eja().random_element()
            sage: A = x.subalgebra_generated_by(orthonormalize_basis=True)
            sage: actual = A.one().operator().matrix()
            sage: expected = matrix.identity(A.base_ring(), A.dimension())
            sage: actual == expected
            True

        """
        if self.dimension() == 0:
            return self.zero()

        return self(self.superalgebra().one())