from functools import reduce


def _orthogonalize(cls, *vecs, **kwargs):
    """Apply the Gram-Schmidt orthogonalization procedure
    to vectors supplied in ``vecs``.

    Parameters
    ==========

    vecs
        vectors to be made orthogonal

    normalize : bool
        If ``True``, return an orthonormal basis.

    rankcheck : bool
        If ``True``, the computation does not stop when encountering
        linearly dependent vectors.

        If ``False``, it will raise ``ValueError`` when any zero
        or linearly dependent vectors are found.

    Returns
    =======

    list
        List of orthogonal (or orthonormal) basis vectors.

    Examples
    ========

    >>> from sympy import I, Matrix
    >>> v = [Matrix([1, I]), Matrix([1, -I])]
    >>> Matrix.orthogonalize(*v)
    [Matrix([
    [1],
    [I]]), Matrix([
    [ 1],
    [-I]])]

    See Also
    ========

    MatrixBase.QRdecomposition

    References
    ==========

    .. [1] https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process
    """

    normalize = kwargs.get('normalize', False)
    rankcheck = kwargs.get('rankcheck', False)

    def project(a, b):
        return b * (a.dot(b, hermitian=True) / b.dot(b, hermitian=True))

    def perp_to_subspace(vec, basis):
        """projects vec onto the subspace given
        by the orthogonal basis ``basis``"""

        components = [project(vec, b) for b in basis]

        if len(basis) == 0:
            return vec

        return vec - reduce(lambda a, b: a + b, components)

    ret  = []
    vecs = list(vecs) # make sure we start with a non-zero vector

    while len(vecs) > 0 and vecs[0].is_zero_matrix:
        if rankcheck is False:
            del vecs[0]
        else:
            raise ValueError("GramSchmidt: vector set not linearly independent")

    for vec in vecs:
        perp = perp_to_subspace(vec, ret)

        if not perp.is_zero_matrix:
            ret.append(cls(perp))
        elif rankcheck is True:
            raise ValueError("GramSchmidt: vector set not linearly independent")

    if normalize:
        ret = [vec / vec.norm() for vec in ret]

    return ret