GA.py

#!/usr/bin/python

"""
The module GA implements symbolic Geometric Algebra in python.
The relevant references for this module are:

    1. "Geometric Algebra for Physicists" by C. Doran and A. Lazenby,
       Cambridge University Press, 2003.

    2. "Geometric Algebra for Computer Science" by Leo Dorst,
       Daniel Fontijne, and Stephen Mann, Morgan Kaufmann Publishers, 2007.
"""
import sys
import numpy
import sympy
import re as regrep
import sympy.galgebra.latex_ex
from sympy.core.decorators import deprecated

NUMPAT = regrep.compile( '([\-0-9])|([\-0-9]/[0-9])')
"""Re pattern for rational number"""

ZERO = sympy.Rational(0)
ONE = sympy.Rational(1)
TWO = sympy.Rational(2)
HALF = sympy.Rational(1, 2)

from sympy.core import Pow as pow_type
from sympy import Abs as abs_type
from sympy.core import Mul as mul_type
from sympy.core import Add as add_type


@sympy.vectorize(0)
def substitute_array(array, *args):
    return(array.subs(*args))


def is_quasi_unit_numpy_array(array):
    """
    Determine if a array is square and diagonal with
    entries of +1 or -1.
    """
    shape = numpy.shape(array)
    if len(shape) == 2 and (shape[0] == shape[1]):
        n = shape[0]
        ix = 0
        while ix < n:
            iy = 0
            while iy < ix:
                if array[ix][iy] != ZERO:
                    return(False)
                iy += 1
            if sympy.Abs(array[ix][ix]) != ONE:
                return(False)
            ix += 1
        return(True)
    else:
        return(False)


@deprecated(value="This function is no longer needed.", issue=3379,
            deprecated_since_version="0.7.2")
def set_main(main_program):
    pass


def plist(lst):
    if type(lst) == list:
        for x in lst:
            plist(x)
    else:
        sys.stderr.write(lst + '\n')
    return


def numeric(num_str):
    """
    Returns rational numbers compatible with symbols.
    Input is a string representing a fraction or integer
    or a simple integer.
    """
    if type(num_str) == int:
        a = num_str
        b = 1
    else:
        tmp = num_str.split('/')
        if len(tmp) == 1:
            a = int(tmp[0])
            b = 1
        else:
            a = int(tmp[0])
            b = int(tmp[1])
    return(sympy.Rational(a, b))


def collect(expr, lst):
    """
    Wrapper for sympy.collect.
    """
    lst = MV.scalar_to_symbol(lst)
    return(sympy.collect(expr, lst))


def sqrt(expr):
    return(sympy.sqrt(expr))


def isint(a):
    """
    Test for integer.
    """
    return(type(a) == int)


def make_null_array(n):
    """
    Return list of n empty lists.
    """
    a = []
    for i in range(n):
        a.append([])
    return(a)


def test_int_flgs(lst):
    """
    Test if all elements in list are 0.
    """
    for i in lst:
        if i:
            return(1)
    return(0)


def comb(N, P):
    """
    Calculates the combinations of the integers [0,N-1] taken P at a time.
    The output is a list of lists of integers where the inner lists are
    the different combinations. Each combination is sorted in ascending
    order.
    """
    def rloop(n, p, combs, comb):
        if p:
            for i in range(n):
                newcomb = comb + [i]
                np = p - 1
                rloop(i, np, combs, newcomb)
        else:
            combs.append(comb)
    combs = []
    rloop(N, P, combs, [])
    for comb in combs:
        comb.sort()
    return(combs)


def diagpq(p, q=0):
    """
    Returns string equivalent metric tensor for signature (p, q).
    """
    n = p + q
    D = []
    for i in xrange(p):
        D.append((i*'0 ' +'1 '+ (n-i-1)*'0 ')[:-1])
    for i in xrange(p,n):
        D.append((i*'0 ' +'-1 '+ (n-i-1)*'0 ')[:-1])
    return ','.join(D)

def arbitrary_metric(n):
    """
    Returns string equivalent metric tensor for arbitrary signature.
    """
    return ','.join(n*[(n*'# ')[:-1]])

def arbitrary_metric_conformal(n):
    """
    Returns string equivalent metric tensor for arbitrary signature (n+1,1).
    """
    str1 = ','.join(n*[n*'# '+'0 0'])
    return ','.join([str1, n*'0 '+'1 0', n*'0 '+'0 -1'])


def make_scalars(symnamelst):
    """
    make_scalars takes a string of symbol names separated by
    blanks and converts them to MV scalars and returns a list
    of the symbols.
    """
    symlst = sympy.symbols(symnamelst)
    scalar_lst = []
    for s in symlst:
        tmp = MV(s, 'scalar')
        scalar_lst.append(tmp)
    return(scalar_lst)


@deprecated(useinstead="sympy.symbols()", issue=3379,
            deprecated_since_version="0.7.2")
def make_symbols(symnamelst):
    return sympy.symbols(symnamelst)


def israt(numstr):
    """
    Test if string represents a rational number.
    """
    global NUMPAT
    if NUMPAT.match(numstr):
        return(1)
    return(0)


def dualsort(lst1, lst2):
    """
    Inplace dual sort of lst1 and lst2 keyed on sorted lst1.
    """
    _indices = range(len(lst1))
    _indices.sort(key=lst2.__getitem__)
    lst1[:] = map(lst1.__getitem__, _indices)
    lst2[:] = map(lst2.__getitem__, _indices)
    return


def cp(A, B):
    """
    Calculates the commutator product (A*B-B*A)/2 for
    the objects A and B.
    """
    return(HALF*(A*B - B*A))


def reduce_base(k, base):
    """
    If base is a list of sorted integers [i_1,...,i_R] then reduce_base
    sorts the list [k,i_1,...,i_R] and calculates whether an odd or even
    number of permutations is required to sort the list. The sorted list
    is returned and +1 for even permutations or -1 for odd permutations.
    """
    if k in base:
        return(0, base)
    grade = len(base)
    if grade == 1:
        if k < base[0]:
            return(1, [k, base[0]])
        else:
            return(-1, [base[0], k])
    ilo = 0
    ihi = grade - 1
    if k < base[0]:
        return(1, [k] + base)
    if k > base[ihi]:
        if grade % 2 == 0:
            return(1, base + [k])
        else:
            return(-1, base + [k])
    imid = ihi + ilo
    if grade == 2:
        return(-1, [base[0], k, base[1]])
    while True:
        if ihi - ilo == 1:
            break
        if base[imid] > k:
            ihi = imid
        else:
            ilo = imid
        imid = (ilo + ihi)/2
    if ilo % 2 == 1:
        return(1, base[:ihi] + [k] + base[ihi:])
    else:
        return(-1, base[:ihi] + [k] + base[ihi:])


def sub_base(k, base):
    """
    If base is a list of sorted integers [i_1,...,i_R] then sub_base returns
    a list with the k^th element removed. Note that k=0 removes the first
    element.  There is no test to see if k is in the range of the list.
    """
    n = len(base)
    if n == 1:
        return([])
    if n == 2:
        if k == base[0]:
            return([base[1]])
        else:
            return([base[0]])
    return(base[:k] + base[k + 1:])


def magnitude(vector):
    """
    Calculate magnitude of vector containing trig expressions
    and simplify.  This is a hack because of way he sign of
    magsq is determined and because of the way that absolute
    values are removed.
    """
    magsq = sympy.expand((vector | vector)())
    magsq = sympy.trigsimp(magsq, deep=True, recursive=True)
    #print magsq
    magsq_str = sympy.galgebra.latex_ex.LatexPrinter()._print(magsq)
    if magsq_str[0] == '-':
        magsq = -magsq
    mag = unabs(sqrt(magsq))
    #print mag
    return(mag)


def LaTeX_lst(lst, title=''):
    """
    Output a list in LaTeX format.
    """
    if title != '':
        sympy.galgebra.latex_ex.LaTeX(title)
    for x in lst:
        sympy.galgebra.latex_ex.LaTeX(x)
    return


def unabs(x):
    """
    Remove absolute values from expressions so a = sqrt(a**2).
    This is a hack.
    """
    if type(x) == mul_type:
        y = unabs(x.args[0])
        for yi in x.args[1:]:
            y *= unabs(yi)
        return(y)
    if type(x) == pow_type:
        if x.args[1] == HALF and type(x.args[0]) == add_type:
            return(x)
        y = 1/unabs(x.args[0])
        return(y)
    if len(x.args) == 0:
        return(x)
    if type(x) == abs_type:
        return(x.args[0])
    return(x)


def function_lst(fstr, xtuple):
    sys.stderr.write(fstr + '\n')
    fct_lst = []
    for xstr in fstr.split():
        f = sympy.Function(xstr)(*xtuple)
        fct_lst.append(f)
    return(fct_lst)


def vector_fct(Fstr, x):
    """
    Create a list of functions of arguments x.  One function is
    created for each variable in x.  Fstr is a string that is
    the base name of each function while each function in the
    list is given the name Fstr+'__'+str(x[ix]) so that if
    Fstr = 'f' and str(x[1]) = 'theta' then the LaTeX output
    of the second element in the output list would be 'f^{\\theta}'.
    """
    nx = len(x)
    Fvec = []
    for ix in range(nx):
        ftmp = sympy.Function(Fstr + '__' + sympy.galgebra.latex_ex.LatexPrinter.str_basic(x[ix]))(*tuple(x))
        Fvec.append(ftmp)
    return(Fvec)


def print_lst(lst):
    for x in lst:
        print x
    return


def normalize(elst, nname_lst):
    """
    Normalize a list of vectors and rename the normalized vectors. 'elist' is the list
    (or array) of vectors to be normalized and nname_lst is a list of the names for the
    normalized vectors.  The function returns the numpy arrays enlst and mags containing
    the normalized vectors (enlst) and the magnitudes of the original vectors (mags).
    """
    i = 0
    mags = numpy.array(MV.n*[ZERO], dtype=numpy.object)
    enlst = numpy.array(MV.n*[ZERO], dtype=numpy.object)
    for (e, nname) in zip(elst, nname_lst):
        emag = magnitude(e)
        emaginv = 1/emag
        mags[i] = emag
        enorm = emaginv*e
        enorm.name = nname
        enlst[i] = enorm
        i += 1
    return(enlst, mags)


def build_base(base_index, base_vectors, reverse=False):
    base = base_vectors[base_index[0]]
    if len(base_index) > 1:
        for i in base_index[1:]:
            base = base ^ base_vectors[i]
    if reverse:
        base = base.rev()
    return(base)


class MV(object):

    is_setup = False
    basislabel_lst = 0
    curvilinear_flg = False
    coords = None

    @staticmethod
    def pad_zeros(value, n):
        """
        Pad list with zeros to length n. If length is > n
        truncate list.  Return padded list.
        """
        nvalue = len(value)
        if nvalue < n:
            value = value + (n - nvalue)*[ZERO]
        if nvalue > n:
            value = value[:n]
        return(value)

    @staticmethod
    def define_basis(basis):
        """
        Calculates all the MV static variables needed for
        basis operations.  See reference 5 section 2.
        """
        MV.vbasis = basis
        MV.vsyms = sympy.symbols(MV.vbasis)
        MV.n = len(MV.vbasis)
        MV.nrg = range(MV.n)
        MV.n1 = MV.n + 1
        MV.n1rg = range(MV.n1)
        MV.npow = 2**MV.n
        MV.index = range(MV.n)
        MV.gabasis = [[]]
        MV.basis = (MV.n + 1)*[0]
        MV.basislabel = (MV.n + 1)*[0]
        MV.basis[0] = []
        MV.basislabel[0] = '1'
        MV.basislabel_lst = [['1']]
        MV.nbasis = numpy.array((MV.n + 1)*[1], dtype=numpy.object)
        for igrade in range(1, MV.n + 1):
            tmp = comb(MV.n, igrade)
            MV.gabasis += [tmp]
            ntmp = len(tmp)
            MV.nbasis[igrade] = ntmp
            MV.basis[igrade] = tmp
            gradelabels = []
            gradelabel_lst = []
            for i in range(ntmp):
                tmp_lst = []
                bstr = ''
                for j in tmp[i]:
                    bstr += MV.vbasis[j]
                    tmp_lst.append(MV.vbasis[j])
                gradelabel_lst.append(tmp_lst)
                gradelabels.append(bstr)
            MV.basislabel_lst.append(gradelabel_lst)
            MV.basislabel[igrade] = gradelabels
        MV.basis_map = [{'':0}]
        igrade = 1
        while igrade <= MV.n:
            tmpdict = {}
            bases = MV.gabasis[igrade]
            ibases = 0
            for base in bases:
                tmpdict[str(base)] = ibases
                ibases += 1
            MV.basis_map.append(tmpdict)
            igrade += 1

        if MV.debug:
            print 'basis strings =', MV.vbasis
            print 'basis symbols =', MV.vsyms
            print 'basis labels  =', MV.basislabel
            print 'basis         =', MV.basis
            print 'grades        =', MV.nbasis
            print 'index         =', MV.index
        return

    @staticmethod
    def define_metric(metric):
        """
        Calculates all the MV static variables needed for
        metric operations.  See reference 5 section 2.
        """
        if MV.metric_str:
            MV.g = []
            MV.metric = numpy.array(MV.n*[MV.n*[ZERO]], dtype=numpy.object)
            if metric == '':
                metric = numpy.array(MV.n*[MV.n*['#']], dtype=numpy.object)
            for i in MV.index:
                for j in MV.index:
                    gij = metric[i][j]
                    if israt(gij):
                        MV.metric[i][j] = numeric(gij)
                    else:
                        if gij == '#':
                            if i == j:
                                gij = '(' + MV.vbasis[j] + '**2)'
                            else:
                                gij = '(' + MV.vbasis[min(
                                    i, j)] + '.' + MV.vbasis[max(i, j)] + ')'
                        tmp = sympy.Symbol(gij)
                        MV.metric[i][j] = tmp
                        if i <= j:
                            MV.g.append(tmp)
        else:
            MV.metric = metric
            MV.g = []
            for row in metric:
                g_row = []
                for col in metric:
                    g_row.append(col)
                    MV.g.append(g_row)
        if MV.debug:
            print 'metric =', MV.metric
        return

    @staticmethod
    def define_reciprocal_frame():
        """
        Calculates unscaled reciprocal vectors (MV.brecp) and scale
        factor (MV.Esq). The ith scaled reciprocal vector is
        (1/MV.Esq)*MV.brecp[i].  The pseudoscalar for the set of
        basis vectors is MV.E.
        """
        if MV.tables_flg:
            MV.E = MV.bvec[0]
            MV.brecp = []
            for i in range(1, MV.n):
                MV.E = MV.E ^ MV.bvec[i]
            for i in range(MV.n):
                tmp = ONE
                if i % 2 != 0:
                    tmp = -ONE
                for j in range(MV.n):
                    if i != j:
                        tmp = tmp ^ MV.bvec[j]
                tmp = tmp*MV.E
                MV.brecp.append(tmp)
            MV.Esq = (MV.E*MV.E)()
            MV.Esq_inv = ONE/MV.Esq
            for i in range(MV.n):
                MV.brecp[i] = MV.brecp[i]*MV.Esq_inv
        return

    @staticmethod
    def reduce_basis_loop(blst):
        """
        Makes one pass through basis product representation for
        reduction of representation to normal form.
        See reference 5 section 3.
        """
        nblst = len(blst)
        if nblst <= 1:
            return(1)
        jstep = 1
        while jstep < nblst:
            istep = jstep - 1
            if blst[istep] == blst[jstep]:
                i = blst[istep]
                if len(blst) > 2:
                    blst = blst[:istep] + blst[jstep + 1:]
                else:
                    blst = []
                if len(blst) <= 1 or jstep == nblst - 1:
                    blst_flg = 0
                else:
                    blst_flg = 1
                return(MV.metric[i][i], blst, blst_flg)
            if blst[istep] > blst[jstep]:
                blst1 = blst[:istep] + blst[jstep + 1:]
                a1 = TWO*MV.metric[blst[jstep]][blst[istep]]
                blst = blst[:istep] + [blst[jstep]] + [blst[istep]] + blst[jstep + 1:]
                if len(blst1) <= 1:
                    blst1_flg = 0
                else:
                    blst1_flg = 1
                return(a1, blst1, blst1_flg, blst)
            jstep += 1
        return(1)

    @staticmethod
    def reduce_basis(blst):
        """
        Repetitively applies reduce_basis_loop to basis
        product representation until normal form is
        realized.  See reference 5 section 3.
        """
        if blst == []:
            blst_coef = []
            blst_expand = []
            for i in MV.n1rg:
                blst_coef.append([])
                blst_expand.append([])
            blst_expand[0].append([])
            blst_coef[0].append(ONE)
            return(blst_coef, blst_expand)
        blst_expand = [blst]
        blst_coef = [ONE]
        blst_flg = [1]
        while test_int_flgs(blst_flg):
            for i in range(len(blst_flg)):
                if blst_flg[i]:
                    tmp = MV.reduce_basis_loop(blst_expand[i])
                    if tmp == 1:
                        blst_flg[i] = 0
                    else:
                        if len(tmp) == 3:
                            blst_coef[i] = tmp[0]*blst_coef[i]
                            blst_expand[i] = tmp[1]
                            blst_flg[i] = tmp[2]
                        else:
                            blst_coef[i] = -blst_coef[i]
                            blst_flg[i] = 1
                            blst_expand[i] = tmp[3]
                            blst_coef.append(-blst_coef[i]*tmp[0])
                            blst_expand.append(tmp[1])
                            blst_flg.append(tmp[2])
        (blst_coef, blst_expand) = MV.combine_common_factors(
            blst_coef, blst_expand)
        return(blst_coef, blst_expand)

    @staticmethod
    def combine_common_factors(blst_coef, blst_expand):
        new_blst_coef = []
        new_blst_expand = []
        for i in range(MV.n1):
            new_blst_coef.append([])
            new_blst_expand.append([])
        nfac = len(blst_coef)
        for ifac in range(nfac):
            blen = len(blst_expand[ifac])
            new_blst_coef[blen].append(blst_coef[ifac])
            new_blst_expand[blen].append(blst_expand[ifac])
        for i in range(MV.n1):
            if len(new_blst_coef[i]) > 1:
                MV.contract(new_blst_coef[i], new_blst_expand[i])
        return(new_blst_coef, new_blst_expand)

    @staticmethod
    def contract(coefs, bases):
        dualsort(coefs, bases)
        n = len(bases) - 1
        i = 0
        while i < n:
            j = i + 1
            if bases[i] == bases[j]:
                coefs[i] += coefs[j]
                bases.pop(j)
                coefs.pop(j)
                n -= 1
            else:
                i += 1
        n = len(coefs)
        i = 0
        while i < n:
            if coefs[i] == ZERO:
                coefs.pop(i)
                bases.pop(i)
                n -= 1
            else:
                i += 1
        return

    @staticmethod
    def convert(coefs, bases):
        mv = MV()
        mv.bladeflg = 0
        for igrade in MV.n1rg:
            coef = coefs[igrade]
            base = bases[igrade]
            if len(coef) > 0:
                nbases = MV.nbasis[igrade]
                mv.mv[igrade] = numpy.array(nbases*[ZERO], dtype=numpy.object)
                nbaserg = range(len(base))
                for ibase in nbaserg:
                    if igrade > 0:
                        k = MV.basis[igrade].index(base[ibase])
                        mv.mv[igrade][k] = coef[ibase]
                    else:
                        mv.mv[igrade] = numpy.array(
                            [coef[0]], dtype=numpy.object)
        return(mv)

    @staticmethod
    def set_str_format(str_mode=0):
        MV.str_mode = str_mode
        return

    @staticmethod
    def str_rep(mv):
        """
         Converts internal representation of a multivector to a string
         for outputting.  If lst_mode = 1, str_rep outputs a list of
         strings where each string contains one multivector coefficient
         concatenated with the corresponding base or blade symbol.

            MV.str_mode     Effect
                0           Print entire multivector on one line (default)
                1           Print each grade on a single line
                2           Print each base on a single line
         """

        if MV.bladeprint:
            mv.convert_to_blades()
            labels = MV.bladelabel
        else:
            if not mv.bladeflg:
                labels = MV.basislabel
            else:
                labels = MV.bladelabel
        mv.compact()
        if isinstance(mv.mv[0], int):
            value = ''
        else:
            value = (mv.mv[0][0]).__str__()
            value = value.replace(' ', '')
        dummy = sympy.Dummy('dummy')
        for igrade in MV.n1rg[1:]:
            if isinstance(mv.mv[igrade], numpy.ndarray):
                j = 0
                for x in mv.mv[igrade]:
                    if x != ZERO:
                        xstr = (x*dummy).__str__()
                        xstr = xstr.replace(' ', '')
                        if xstr[0] != '-' and len(value) > 0:
                            xstr = '+' + xstr
                        if xstr.find('_dummy') < 2 and xstr[-5:] != '_dummy':
                            xstr = xstr.replace(
                                '_dummy*', '') + '*' + labels[igrade][j]
                        else:
                            xstr = xstr.replace('_dummy', labels[igrade][j])
                        if MV.str_mode == 2:
                            xstr += '\n'
                        value += xstr
                    j += 1
                if MV.str_mode == 1:
                    value += '\n'
        if value == '':
            value = '0'
        #value = value.replace(' ','')
        value = value.replace('dummy', '1')
        return(value)

    @staticmethod
    def xstr_rep(mv, lst_mode=0):
        """
         Converts internal representation of a multivector to a string
         for outputing.  If lst_mode = 1, str_rep outputs a list of
         strings where each string contains one multivector coefficient
         concatenated with the corresponding base or blade symbol.
         """
        if lst_mode:
            outlst = []
        if MV.bladeprint:
            mv.convert_to_blades()
            labels = MV.bladelabel
        else:
            if not mv.bladeflg:
                labels = MV.basislabel
            else:
                labels = MV.bladelabel
        value = ''
        for igrade in MV.n1rg:
            tmp = []
            if isinstance(mv.mv[igrade], numpy.ndarray):
                j = 0
                for x in mv.mv[igrade]:
                    if x != ZERO:
                        xstr = x.__str__()
                        if xstr == '+1' or xstr == '1' or xstr == '-1':
                            if xstr == '+1' or xstr == '1':
                                xstr = '+'
                            else:
                                xstr = '-'
                        else:
                            if xstr[0] != '+':
                                xstr = '+(' + xstr + ')'
                            else:
                                xstr = '+(' + xstr[1:] + ')'
                        value += xstr + labels[igrade][j]
                        if MV.str_mode and not lst_mode:
                            value += value + '\n'
                        if lst_mode:
                            tmp.append(value)
                    j += 1
            if lst_mode:
                if len(tmp) > 0:
                    outlst.append(tmp)
                value = ''
        if not lst_mode:
            if len(value) > 1 and value[0] == '+':
                value = value[1:]
            if len(value) == 0:
                value = '0'
        else:
            value = outlst
        return(value)

    @staticmethod
    def setup(basis, metric='', rframe=False, coords=None, debug=False, offset=0):
        """
        MV.setup initializes the MV class by calculating the static
        multivector tables required for geometric algebra operations
        on multivectors.  See reference 5 section 2 for details on
        basis and metric arguments.
        """
        MV.is_setup = True
        MV.metric_str = False
        MV.debug = debug
        MV.bladeprint = 0
        MV.tables_flg = 0
        MV.str_mode = 0
        MV.basisroot = ''
        MV.index_offset = offset
        if coords is None:
            MV.coords = None
        else:
            MV.coords = tuple(coords)
            rframe = True
        if type(basis) == str:
            basislst = basis.split()
            if len(basislst) == 1:
                MV.basisroot = basislst[0]
                basislst = []
                for coord in coords:
                    basislst.append(MV.basisroot + '_' + str(coord))
            MV.define_basis(basislst)
        if type(metric) == str:
            MV.metric_str = True
            if len(metric) > 0:
                if metric[0] == '[' and metric[-1] == ']':
                    tmps = metric[1:-1].split(',')
                    N = len(tmps)
                    metric = []
                    itmp = 0
                    for tmp in tmps:
                        xtmp = N*['0']
                        xtmp[itmp] = tmp
                        itmp += 1
                        metric.append(xtmp)
                else:
                    tmps = metric.split(',')
                    metric = []
                    for tmp in tmps:
                        xlst = tmp.split()
                        xtmp = []
                        for x in xlst:
                            xtmp.append(x)
                        metric.append(xtmp)

        MV.define_metric(metric)
        MV.multiplication_table()
        MV.blade_table()
        MV.inverse_blade_table()
        MV.tables_flg = 1
        isym = 0
        MV.bvec = []
        for name in MV.vbasis:
            bvar = MV(value=isym, mvtype='basisvector', mvname=name)
            bvar.bladeflg = 1
            MV.bvec.append(bvar)
            isym += 1
        if rframe:
            MV.define_reciprocal_frame()
        MV.I = MV(ONE, 'pseudo', 'I')
        MV.ZERO = MV()
        Isq = (MV.I*MV.I)()
        MV.Iinv = (1/Isq)*MV.I
        return MV.bvec

    @staticmethod
    def set_coords(coords):
        MV.coords = coords
        return

    @staticmethod
    def scalar_fct(fct_name):
        """
        Create multivector scalar function with name fct_name (string) and
        independent variables coords (list of variable).  Default variables are
        those associated with each dimension of vector space.
        """
        phi = sympy.Function(fct_name)(*MV.coords)
        Phi = MV(phi, 'scalar')
        Phi.name = fct_name
        return(Phi)

    @staticmethod
    def vector_fct(fct_name, vars=''):
        """
        Create multivector vector function with name fct_name (string) and
        independent variables coords (list of variable).  Default variables are
        those associated with each dimension of vector space.
        """
        if isinstance(vars, str):
            Acoefs = vector_fct(fct_name, MV.coords)
        else:
            Acoefs = numpy.array(MV.n*[ZERO], dtype=numpy.object)
            x = MV.coords
            if isinstance(vars, sympy.core.symbol.Symbol):
                for icoef in MV.nrg:
                    Acoefs[icoef] = sympy.Function(fct_name + '__' +
                                   sympy.galgebra.latex_ex.LatexPrinter.str_basic(x[icoef]))(vars)
            else:
                for icoef in MV.nrg:
                    Acoefs[icoef] = sympy.Function(fct_name + '__' +
                                   sympy.galgebra.latex_ex.LatexPrinter.str_basic(x[icoef]))(*tuple(vars))
        A = MV(Acoefs, 'vector', fct_name)
        return(A)

    @staticmethod
    def rebase(x, coords, base_name='', debug=False, debug_level=0):
        """
        Define curvilinear coordinates for previously defined vector (multivector) space (MV.setup has been run)
        with position vector, x, that is a vector function of the independent coordinates, coords (list of
        sympy variables equal in length to dimension of vector space), and calculate:

            1. Frame (basis) vectors
            2. Normalized frame (basis) vectors.
            3. Metric tensor
            4. Reciprocal frame vectors
            5. Reciprocal metric tensor
            6. Connection multivectors

        The basis vectors are named with the base_name (string) and a subscript derived from the name of each
        coordinate.  So that if the base name is 'e' and the coordinated are [r,theta,z] the variable names
        of the frame vectors would be e_r, e_theta, and e_z.  For LaTeX output the names of the frame vectors
        would be e_{r}, e_{\theta}, and e_{z}.  Everything needed to compute the geometric, outer, and inner
        derivatives of multivector functions in curvilinear coordinates is calculated.

        If debug is True all the quantities in the above list are output in LaTeX format.

        Currently rebase works with cylindrical and spherical coordinates in any dimension.  The limitation is the
        ability to automatically simplify complex sympy expressions generated while calculating the quantities in
        the above list.  This is why the debug option is included.  The debug_level can equal 0,1,2, or 3 and
        determines how far in the list to calculate (input 0 to do the entire list) while debugging.
        """
        #Form root names for basis, reciprocal basis, normalized basis, and normalized reciprocal basis

        if base_name == '':
            base_name = MV.basisroot + 'prm'

        LaTeX_base = sympy.galgebra.latex_ex.LatexPrinter.extended_symbol(
            base_name)
        bm = '\\bm{' + LaTeX_base + '}'
        bmhat = '\\hat{' + bm + '}'
        bstr = bmhat + '_{[i_{1},\dots, i_{R}]}'
        base_name += 'bm'
        base_name_hat = base_name + 'hat'

        base_name_lst = []
        nbase_name_lst = []
        rbase_name_lst = []
        rnbase_name_lst = []
        coords_lst = []

        for coord in coords:
            coord_str = sympy.galgebra.latex_ex.LatexPrinter.str_basic(coord)
            coords_lst.append(coord_str)
            base_name_lst.append(base_name + '_' + coord_str)
            rbase_name_lst.append(base_name + '__' + coord_str)
            nbase_name_lst.append(base_name_hat + '_' + coord_str)
            rnbase_name_lst.append(base_name_hat + '__' + coord_str)

        if not (MV.n == len(coords) == len(base_name_lst)):
            print 'rebaseMV inputs not congruent:'
            print 'MV.n =', MV.n
            print 'coords =', coords
            print 'bases =', base_name
            sys.exit(1)

        if isinstance(x, MV):

            #Calculate basis vectors from derivatives of position vector x

            bases = numpy.array(MV.n*[ZERO], dtype=numpy.object)
            i = 0
            for coord in coords:
                ei = x.diff(coords[i])
                ei.set_name(base_name_lst[i])
                bases[i] = ei
                i += 1

            #Calculate normalizee basis vectors and basis vector magnitudes

            if debug:
                print 'Coordinate Generating Vector'
                print x
                print 'Basis Vectors'
                for base in bases:
                    print base

        else:

            #Input basis vectors as N vector fields

            bases = x

            for (base, name) in zip(bases, base_name_lst):
                base.set_name(name)

            if debug:
                print 'Basis Vectors'
                for base in bases:
                    print base
                if debug_level == 1:
                    return

        if debug_level == 1:
            return

        #Calculate normalized basis vectors and magnitudes of
        #unormalized basis vectors

        (nbases, mags) = normalize(bases, nbase_name_lst)

        if debug:
            print 'Magnitudes'
            print '\\abs{' + LaTeX_base + '_{i}} = ', mags
            print 'Normalized Basis Vectors'
            for nbase in nbases:
                print nbase

        g = numpy.array(MV.n*[MV.n*[ZERO]], dtype=numpy.object)

        for irow in MV.nrg:
            for icol in MV.nrg:
                magsq = sympy.expand((nbases[irow] | nbases[icol])())
                g[irow][icol] = sympy.simplify(
                    sympy.trigsimp(magsq, deep=True, recursive=True))

        if debug:
            print 'Metric $\\hat{g}_{ij} = \\hat{' + LaTeX_base + \
                '}_{i}\\cdot \\hat{' + LaTeX_base + '}_{j}$'
            print r'\hat{g}_{ij} =', sympy.galgebra.latex_ex.LaTeX(g)

        if debug_level == 2:
            return

        #Calculate reciprocal normalized basis vectors

        rnbases = []

        if is_quasi_unit_numpy_array(g):
            ibasis = 0
            while ibasis < MV.n:
                base = g[ibasis][ibasis]*nbases[ibasis]
                base.set_name(rnbase_name_lst[ibasis])
                base.simplify()
                base.trigsimp()
                rnbases.append(base)
                ibasis += 1
        else:
            rnbases = reciprocal_frame(nbases, rnbase_name_lst)
            ibase = 0
            for base in rnbases:
                base.simplify()
                base.trigsimp()
                rnbases[ibase] = base
                ibase += 1

        if debug:
            if debug_level != 0:
                sympy.galgebra.latex_ex.MV_format(1)
            print 'Reciprocal Normalized Basis Vectors'
            for rnbase in rnbases:
                print rnbase

        if debug_level == 3:
            return

        #Calculate components of inverse vectors

        Acoef = []

        for ibasis in MV.nrg:
            evec = numpy.array(MV.n*[ZERO], dtype=numpy.object)
            for jbasis in MV.nrg:
                evec[jbasis] = (MV.bvec[ibasis] | rnbases[jbasis])()
            Acoef.append(evec)

        #Calculat metric tensors

        gr = numpy.array(MV.n*[MV.n*[ZERO]], dtype=numpy.object)

        for irow in MV.nrg:
            for icol in MV.nrg:
                magsq = sympy.expand((rnbases[irow] | rnbases[icol])())
                gr[irow][icol] = sympy.simplify(
                    sympy.trigsimp(magsq, deep=True, recursive=True))

        if debug:
            print 'Metric $\\hat{g}^{ij} = \\hat{' + LaTeX_base + \
                '}^{i}\\cdot \\hat{' + LaTeX_base + '}^{j}$'
            print r'\hat{g}^{ij} =', sympy.galgebra.latex_ex.LaTeX(gr)

        if debug_level == 4:
            return

        #Calculate bases and reciprocal bases for curvilinear mulitvectors

        MV_bases = [[ONE]]
        MV_rbases = [[ONE]]
        igrade = 1
        while igrade <= MV.n:
            base_index = MV.gabasis[igrade]
            grade_bases = []
            rgrade_bases = []
            for index in base_index:
                base = build_base(index, nbases)
                base.simplify()
                base.trigsimp()
                rbase = build_base(index, rnbases, True)
                rbase.simplify()
                rbase.trigsimp()
                grade_bases.append(base)
                rgrade_bases.append(rbase)
            igrade += 1
            MV_bases.append(grade_bases)
            MV_rbases.append(rgrade_bases)

        #Calculate connection multivectors for geometric derivative

        MV_connect = [[ZERO]]
        igrade = 1
        while igrade <= MV.n:
            grade_connect = []
            ibase = 0
            for base in MV_bases[igrade]:
                sum = MV()
                itheta = 0
                for (theta, etheta) in zip(coords, rnbases):
                    psum = (1/mags[itheta])*etheta*base.diff(theta)
                    psum.trigsimp()
                    sum += psum
                    itheta += 1
                sum.simplify()
                sum.trigsimp()
                grade_connect.append(sum)
                ibase += 1
            MV_connect.append(grade_connect)
            igrade += 1

        if debug:
            print 'Curvilinear Bases: $' + bstr + ' = ' + bmhat + \
                '_{i_{1}}\\W\\dots\\W' + bmhat + '_{i_{R}}$'
            igrade = 1
            for grade in MV_bases[1:]:
                ibase = 0
                for base in grade:
                    index = MV.gabasis[igrade][ibase]
                    sub_str = ''
                    for i in index:
                        sub_str += sympy.galgebra.latex_ex.LatexPrinter.extended_symbol(coords_lst[i])
                    base_str = bmhat + '_{[' + sub_str + ']} = '
                    print base_str, base
                    ibase += 1
                igrade += 1

        if debug_level == 5:
            return

        #Calculate representation of connection multivectors in curvilinear system

        MV_Connect = [[ZERO]]
        igrade = 1
        while igrade <= MV.n:
            grade_connect = []
            ibase = 0
            nbase = len(MV_bases[igrade])

            if igrade < MV.n:
                ibase = 0
                p1base = len(MV_bases[igrade + 1])
                m1base = len(MV_bases[igrade - 1])
                while ibase < nbase:
                    Cm1 = numpy.array(m1base*[ZERO], dtype=numpy.object)
                    Cp1 = numpy.array(p1base*[ZERO], dtype=numpy.object)
                    C = MV_connect[igrade][ibase]
                    if igrade == 1:
                        X = C(0)
                    else:
                        X = MV()
                    jbase = 0
                    while jbase < m1base:
                        Cm1[jbase] = sympy.trigsimp((MV.inner_product(MV_rbases[igrade - 1][jbase], C))(), deep=True, recursive=True)
                        jbase += 1
                    jbase = 0
                    while jbase < p1base:
                        Cp1[jbase] = sympy.trigsimp((MV.inner_product(MV_rbases[igrade + 1][jbase], C))(), deep=True, recursive=True)
                        jbase += 1
                    X += MV(
                        (igrade - 1, Cm1), 'grade') + MV((igrade + 1, Cp1), 'grade')
                    X.simplify()
                    X.trigsimp()
                    grade_connect.append(X)
                    ibase += 1
            else:
                ibase = 0
                m1base = len(MV_bases[igrade - 1])
                while ibase < nbase:
                    Cm1 = numpy.array(m1base*[ZERO], dtype=numpy.object)
                    C = MV_connect[igrade][ibase]
                    jbase = 0
                    while jbase < m1base:
                        Cm1[jbase] = sympy.trigsimp((MV.inner_product(MV_rbases[igrade - 1][jbase], C))(), deep=True, recursive=True)
                        jbase += 1
                    X = MV()
                    X.mv[MV.n - 1] = Cm1
                    X.simplify()
                    X.trigsimp()
                    grade_connect.append(X)
                    ibase += 1

            MV_Connect.append(grade_connect)
            igrade += 1

        base_str = ''
        for coord in coords:
            base_str += base_name + '_' + \
                sympy.galgebra.latex_ex.LatexPrinter.str_basic(coord) + ' '
        base_str = base_str[:-1]

        old_names = MV.vbasis

        MV.setup(base_str, g, True, coords)

        MV.curvilinear_flg = True
        MV.Connect = MV_Connect
        sympy.galgebra.latex_ex.LatexPrinter.latex_bases()
        MV.Rframe = numpy.array(MV.n*[ZERO], dtype=numpy.object)
        ibasis = 0
        while ibasis < MV.n:
            base = MV()
            jbasis = 0
            while jbasis < MV.n:
                base.add_in_place(gr[ibasis][jbasis]*MV.bvec[jbasis])
                jbasis += 1
            base.scalar_mul_inplace(1/mags[ibasis])
            MV.Rframe[ibasis] = base
            ibasis += 1

        MV.org_basis = []
        for ibasis in MV.nrg:
            evec = MV(Acoef[ibasis], 'vector', old_names[ibasis])
            MV.org_basis.append(evec)

        if MV.coords[0] == sympy.Symbol('t'):
            MV.dedt = []
            # I can't fix this no name error because I have no clue what the
            # original writer was trying to do.
            for coef in dedt_coef:
                MV.dedt.append(MV(coef, 'vector'))
        else:
            MV.dedt = None

        if debug:
            print 'Representation of Original Basis Vectors'
            for evec in MV.org_basis:
                print evec

            print 'Renormalized Reciprocal Vectors ' + '$\\bfrac{' + bmhat + \
                '^{k}}{\\abs{\\bm{' + LaTeX_base + '}_{k}}}$'

            ibasis = 0
            while ibasis < MV.n:
                c_str = sympy.galgebra.latex_ex.LatexPrinter.extended_symbol(
                    coords_lst[ibasis])
                print '\\bfrac{\\bm{\\hat{' + LaTeX_base + '}}^{' + c_str + \
                    '}}{\\abs{\\bm{' + LaTeX_base + '}_{' + c_str + '}}} =', \
                    MV.Rframe[ibasis]
                ibasis += 1

            title_str = 'Connection Multivectors: $C\\lbrc' + bstr + \
                '\\rbrc = ' + '\\bfrac{' + bmhat + '^{k}}{\\abs{' + bmhat + \
                '_{k}}}\\pdiff{' + bstr + '}{\\theta^{k}}$'

            print title_str
            igrade = 1
            for grade in MV.Connect[1:]:
                ibase = 0
                for base in grade:
                    index = MV.gabasis[igrade][ibase]
                    sub_str = ''
                    for i in index:
                        sub_str += sympy.galgebra.latex_ex.LatexPrinter.extended_symbol(coords_lst[i])

                    base_str = 'C\\lbrc\\hat{' + \
                        LaTeX_base + '}_{[' + sub_str + ']}\\rbrc = '
                    print base_str, base
                    ibase += 1
                igrade += 1
        return

    @staticmethod
    def print_blades():
        """
        Set multivector output to blade representation.
        """
        MV.bladeprint = 1
        return

    @staticmethod
    def print_bases():
        """
        Set multivector output to base representation.
        """
        MV.bladeprint = 0
        return

    @staticmethod
    def multiplication_table():
        """
        Calculate geometric product base multiplication table.
        See reference 5 section 3 for details.
        """
        MV.mtable = []
        for igrade in MV.n1rg:
            MV.mtable.append([])
            for ibase in range(MV.nbasis[igrade]):
                MV.mtable[igrade].append([])
                if igrade == 0:
                    base1 = []
                else:
                    base1 = MV.basis[igrade][ibase]
                for jgrade in MV.n1rg:
                    MV.mtable[igrade][ibase].append([])
                    for jbase in range(MV.nbasis[jgrade]):
                        if jgrade == 0:
                            base2 = []
                        else:
                            base2 = MV.basis[jgrade][jbase]
                        base = base1 + base2
                        (coefs, bases) = MV.reduce_basis(base)
                        product = MV.convert(coefs, bases)
                        product.name = '(' + MV.basislabel[igrade][ibase] + \
                            ')(' + MV.basislabel[jgrade][jbase] + ')'
                        MV.mtable[igrade][ibase][jgrade].append(product)
        if MV.debug:
            print 'Multiplication Table:'
            for level1 in MV.mtable:
                for level2 in level1:
                    for level3 in level2:
                        for mv in level3:
                            mv.printmv()
        return

    @staticmethod
    def geometric_product(mv1, mv2):
        """
        MV.geometric_product(mv1,mv2) calculates the geometric
        product the multivectors mv1 and mv2 (mv1*mv2). See
        reference 5 section 3.
        """
        product = MV()
        if isinstance(mv1, MV) and isinstance(mv2, MV):
            bladeflg1 = mv1.bladeflg
            bladeflg2 = mv2.bladeflg
            if bladeflg1:
                mv1.convert_from_blades()
            if bladeflg2:
                mv2.convert_from_blades()
            for igrade in MV.n1rg:
                gradei = mv1.mv[igrade]
                if isinstance(gradei, numpy.ndarray):
                    for ibase in range(MV.nbasis[igrade]):
                        xi = gradei[ibase]
                        if xi != ZERO:
                            for jgrade in MV.n1rg:
                                gradej = mv2.mv[jgrade]
                                if isinstance(gradej, numpy.ndarray):
                                    for jbase in range(MV.nbasis[jgrade]):
                                        xj = gradej[jbase]
                                        if xj != ZERO:
                                            xixj = MV.mtable[igrade][ibase][jgrade][jbase].scalar_mul(xi*xj)
                                            product.add_in_place(xixj)
            product.bladeflg = 0
            if bladeflg1:
                mv1.convert_to_blades()
            if bladeflg2:
                mv1.convert_to_blades()
            if bladeflg1 and bladeflg2:
                product.convert_to_blades()
        else:
            if isinstance(mv1, MV):
                product = mv1.scalar_mul(mv2)
            else:
                product = mv2.scalar_mul(mv1)
        return(product)

    @staticmethod
    def wedge(igrade1, blade1, vector2, name=''):
        """
        Calculate the outer product of a multivector blade
        and a vector.  See reference 5 section 5 for details.
        """
        w12 = blade1*vector2
        w21 = vector2*blade1
        if igrade1 % 2 != 0:
            w = w12 - w21
        else:
            w = w12 + w21
        w.name = name
        return(w*HALF)

    @staticmethod
    def blade_table():
        """
        Calculate basis blades in terms of bases. See reference 5
        section 5 for details. Used to convert from blade to base
        representation.
        """
        MV.btable = []
        MV.bladelabel = []
        basis_str = MV.basislabel[1]
        for igrade in MV.n1rg:
            MV.bladelabel.append([])
            if igrade == 0:
                MV.bladelabel[0].append('1')
                tmp = [MV(value=ONE, mvtype='scalar', mvname='1')]
            if igrade == 1:
                tmp = []
                for ibase in range(MV.nbasis[1]):
                    MV.bladelabel[1].append(basis_str[ibase])
                    tmp.append(MV(value=ibase,
                               mvtype='basisvector', mvname=basis_str[ibase]))
            if igrade >= 2:
                tmp = []
                basis = MV.basis[igrade]
                for blade in basis:
                    name = ''
                    for i in blade:
                        name += basis_str[i] + '^'
                    name = name[:-1]
                    MV.bladelabel[igrade].append(name)
                    lblade = MV.basis[igrade - 1].index(blade[:-1])
                    rblade = blade[-1]
                    igrade1 = igrade - 1
                    blade1 = MV.btable[igrade1][lblade]
                    vector2 = MV.btable[1][rblade]
                    b1Wv2 = MV.wedge(igrade1, blade1, vector2, name)
                    tmp.append(b1Wv2)
            MV.btable.append(tmp)
        if MV.debug:
            print 'Blade Tabel:'
            for grade in MV.btable:
                for mv in grade:
                    print mv
            print 'Blade Labels:'
            print MV.bladelabel
        return

    @staticmethod
    def inverse_blade_table():
        """
        Calculate bases in terms of basis blades. See reference 5
        section 5 for details. Used to convert from base to blade
        representation.
        """
        MV.ibtable = []
        for igrade in MV.n1rg:
            if igrade == 0:
                tmp = [MV(value=ONE, mvtype='scalar', mvname='1')]
            if igrade == 1:
                tmp = []
                for ibase in range(MV.nbasis[1]):
                    tmp.append(MV(value=ibase, mvtype='basisvector'))
            if igrade >= 2:
                tmp = []
                iblade = 0
                for blade in MV.btable[igrade]:
                    invblade = MV()
                    for i in range(igrade - 1):
                        invblade.mv[i] = -blade.mv[i]
                        invblade.mv[igrade] = +blade.mv[igrade]
                        invblade.bladeflg = 1
                    if igrade >= 4:
                        jgrade = igrade - 2
                        while jgrade > 1:
                            for ibase in range(MV.nbasis[jgrade]):
                                invblade.substitute_base(
                                    jgrade, ibase, MV.ibtable[jgrade][ibase])
                            jgrade -= 2
                    invblade.name = MV.basislabel[igrade][iblade]
                    iblade += 1
                    tmp.append(invblade)
            MV.ibtable.append(tmp)
        if MV.debug:
            print 'Inverse Blade Tabel:'
            for grade in MV.ibtable:
                for mv in grade:
                    mv.printmv()
        return

    @staticmethod
    def outer_product(mv1, mv2):
        """
        MV.outer_product(mv1,mv2) calculates the outer (exterior,wedge)
        product of the multivectors mv1 and mv2 (mv1^mv2). See
        reference 5 section 6.
        """
        if type(mv1) == type(MV) and type(mv2) == type(MV):
            if mv1.is_scalar() and mv2.is_scalar():
                return(mv1*mv2)

        if isinstance(mv1, MV) and isinstance(mv2, MV):
            product = MV()
            product.bladeflg = 1
            mv1.convert_to_blades()
            mv2.convert_to_blades()
            for igrade1 in MV.n1rg:
                if isinstance(mv1.mv[igrade1], numpy.ndarray):
                    pg1 = mv1.project(igrade1)
                    for igrade2 in MV.n1rg:
                        igrade = igrade1 + igrade2
                        if igrade <= MV.n:
                            if isinstance(mv2.mv[igrade2], numpy.ndarray):
                                pg2 = mv2.project(igrade2)
                                pg1pg2 = pg1*pg2
                                product.add_in_place(pg1pg2.project(igrade))
        else:
            if isinstance(mv1, MV):
                product = mv1.scalar_mul(mv2)
            if isinstance(mv2, MV):
                product = mv2.scalar_mul(mv1)
        return(product)

    @staticmethod
    def inner_product(mv1, mv2, mode='s'):
        """
        MV.inner_product(mv1,mv2) calculates the inner

        mode = 's' - symmetric (Doran & Lasenby)
        mode = 'l' - left contraction (Dorst)
        mode = 'r' - right contraction (Dorst)
        """
        if type(mv1) == type(MV) and type(mv2) == type(MV):
            if mv1.is_scalar() and mv2.is_scalar():
                return(mv1*mv2)

        if isinstance(mv1, MV) and isinstance(mv2, MV):
            product = MV()
            product.bladeflg = 1
            mv1.convert_to_blades()
            mv2.convert_to_blades()
            for igrade1 in range(MV.n1):
                if isinstance(mv1.mv[igrade1], numpy.ndarray):
                    pg1 = mv1.project(igrade1)
                    for igrade2 in range(MV.n1):
                        igrade = igrade1 - igrade2
                        if mode == 's':
                            igrade = igrade.__abs__()
                        else:
                            if mode == 'l':
                                igrade = -igrade
                        if igrade >= 0:
                            if isinstance(mv2.mv[igrade2], numpy.ndarray):
                                pg2 = mv2.project(igrade2)
                                pg1pg2 = pg1*pg2
                                product.add_in_place(pg1pg2.project(igrade))
            return(product)
        else:
            if mode == 's':
                if isinstance(mv1, MV):
                    product = mv1.scalar_mul(mv2)
                if isinstance(mv2, MV):
                    product = mv2.scalar_mul(mv1)
            else:
                product = None
        return(product)

    @staticmethod
    def addition(mv1, mv2):
        """
        MV.addition(mv1,mv2) calculates the sum
        of the multivectors mv1 and mv2 (mv1+mv2).
        """
        sum = MV()
        if isinstance(mv1, MV) and isinstance(mv2, MV):
            if mv1.bladeflg or mv2.bladeflg:
                mv1.convert_to_blades()
                mv2.convert_to_blades()
                sum.bladeflg = 1
            for i in MV.n1rg:
                if isinstance(mv1.mv[i], numpy.ndarray) and isinstance(mv2.mv[i], numpy.ndarray):
                    sum.mv[i] = mv1.mv[i] + mv2.mv[i]
                else:
                    if isinstance(mv1.mv[i], numpy.ndarray) and not isinstance(mv2.mv[i], numpy.ndarray):
                        sum.mv[i] = +mv1.mv[i]
                    else:
                        if isinstance(mv2.mv[i], numpy.ndarray) and not isinstance(mv1.mv[i], numpy.ndarray):
                            sum.mv[i] = +mv2.mv[i]
            return(sum)
        else:
            if isinstance(mv1, MV):
                return(mv1 + MV(mv2, 'scalar'))
            else:
                return(MV(mv1, 'scalar') + mv2)

    @staticmethod
    def subtraction(mv1, mv2):
        """
        MV.subtraction(mv1,mv2) calculates the difference
        of the multivectors mv1 and mv2 (mv1-mv2).
        """
        diff = MV()
        if isinstance(mv1, MV) and isinstance(mv2, MV):
            if mv1.bladeflg or mv2.bladeflg:
                mv1.convert_to_blades()
                mv2.convert_to_blades()
                diff.bladeflg = 1
            for i in MV.n1rg:
                if isinstance(mv1.mv[i], numpy.ndarray) and isinstance(mv2.mv[i], numpy.ndarray):
                    diff.mv[i] = mv1.mv[i] - mv2.mv[i]
                else:
                    if isinstance(mv1.mv[i], numpy.ndarray) and not isinstance(mv2.mv[i], numpy.ndarray):
                        diff.mv[i] = +mv1.mv[i]
                    else:
                        if not isinstance(mv1.mv[i], numpy.ndarray) and isinstance(mv2.mv[i], numpy.ndarray):
                            diff.mv[i] = -mv2.mv[i]
            return(diff)
        else:
            if isinstance(mv1, MV):
                return(mv1 - MV(mv2, 'scalar'))
            else:
                return(MV(mv1, 'scalar') - mv2)

    @staticmethod
    def vdiff(vec, x):
        dvec = numpy.array(len(vec)*[ZERO])
        ivec = 0
        for veci in vec:
            dvec[ivec] = sympy.diff(veci, x)
            ivec += 1
        return(dvec)

    @staticmethod
    def scalar_to_symbol(scalar):
        if isinstance(scalar, MV):
            return(scalar.mv[0][0])
        if type(scalar) == list:
            sym = []
            for x in scalar:
                sym.append(MV.scalar_to_symbol(x))
            return(sym)
        return(scalar)

    def __init__(self, value='', mvtype='', mvname='', fct=False, vars=None):
        """
        Initialization of multivector X. Inputs are as follows

        mvtype           value                       result

        default         default                  Zero multivector
        'basisvector'   int i                    ith basis vector
        'basisbivector' int i                    ith basis bivector
        'scalar'        symbol x                 scalar of value x
                        string s
        'grade'         [int i, symbol array A]  X.grade(i) = A
                        [int i, string s]
        'vector'        symbol array A           X.grade(1) = A
                        string s
        'grade2'        symbol array A           X.grade(2) = A
                        string s
        'pseudo'        symbol x                 X.grade(n) = x
                        string s
        'spinor'        string s                 spinor with coefficients
                                                 s__indices and name sbm

        mvname is name of multivector.
        If fct is 'True' and MV.coords is defined in MV.setup then a
        multivector field of MV.coords is instantiated.
        """

        self.name = mvname
        self.mv = MV.n1*[0]
        self.bladeflg = 0  # 1 for blade expansion
        self.puregrade = 1
        if mvtype == 'basisvector':
            self.mv[1] = numpy.array(MV.nbasis[1]*[ZERO], dtype=numpy.object)
            self.mv[1][value] = ONE
        if mvtype == 'basisbivector':
            self.mv[2] = numpy.array(MV.nbasis[2]*[ZERO], dtype=numpy.object)
            self.mv[2][value] = ONE
        if mvtype == 'scalar':
            if isinstance(value, str):
                value = sympy.Symbol(value)
            if isinstance(value, int):
                value = sympy.Rational(value)
            self.mv[0] = numpy.array([value], dtype=numpy.object)
        if mvtype == 'pseudo':
            if isinstance(value, str):
                value = sympy.Symbol(value)
            self.mv[MV.n] = numpy.array([value], dtype=numpy.object)
        if mvtype == 'vector':
            if isinstance(value, str):  # Most general vector
                symbol_str = ''
                for ibase in MV.nrg:
                    if MV.coords is None:
                        symbol = value + '__' + str(ibase + MV.index_offset)
                        symbol_str += symbol + ' '
                    else:
                        symbol = value + '__' + (MV.coords[ibase]).name
                        symbol_str += symbol + ' '
                symbol_lst = sympy.symbols(symbol_str)
                self.mv[1] = numpy.array(symbol_lst, dtype=numpy.object)
                self.name = value
            else:
                value = MV.pad_zeros(value, MV.nbasis[1])
                self.mv[1] = numpy.array(value, dtype=numpy.object)
        if mvtype == 'grade2':
            if isinstance(value, str):  # Most general grade-2 multivector
                if value != '':
                    symbol_str = ''
                    for base in MV.basis[2]:
                        symbol = value + MV.construct_index(base)
                        symbol_str += symbol + ' '
                    symbol_lst = sympy.symbols(symbol_str)
                    self.mv[2] = numpy.array(symbol_lst, dtype=numpy.object)
            else:
                value = MV.pad_zeros(value, MV.nbasis[2])
                self.mv[2] = numpy.array(value, dtype=numpy.object)
        if mvtype == 'grade':
            igrade = value[0]
            coefs = value[1]
            if isinstance(coefs, str):  # Most general pure grade multivector
                base_symbol = coefs
                coefs = []
                bases = MV.basis[igrade]
                if igrade == 0:
                    self.mv[0] = numpy.array(
                        [sympy.Symbol(base_symbol)], dtype=numpy.object)
                else:
                    for base in bases:
                        coef = base_symbol + MV.construct_index(base)
                        coef = sympy.Symbol(coef)
                        coefs.append(coef)
                    self.mv[igrade] = numpy.array(coefs, dtype=numpy.object)
            else:
                self.mv[igrade] = coefs
        if mvtype == 'base':
            self.mv[value[0]] = numpy.array(
                MV.nbasis[value[0]]*[ZERO], dtype=numpy.object)
            self.mv[value[0]][value[1]] = ONE
        if mvtype == 'spinor':
            if isinstance(value, str):  # Most general spinor
                for grade in MV.n1rg:
                    if grade % 2 == 0:
                        symbol_str = ''
                        if grade != 0:
                            symbol_lst = []
                            for base in MV.basis[grade]:
                                symbol = sympy.Symbol(
                                    value + MV.construct_index(base))
                                symbol_lst.append(symbol)
                            self.mv[grade] = numpy.array(
                                symbol_lst, dtype=numpy.object)
                        else:
                            self.mv[0] = numpy.array(
                                [sympy.Symbol(value)], dtype=numpy.object)
                self.name = value + 'bm'
        if isinstance(value, str) and mvtype == '':  # Most general multivector
            if value != '':
                for grade in MV.n1rg:
                    symbol_str = ''
                    if grade != 0:
                        symbol_lst = []
                        for base in MV.basis[grade]:
                            symbol = sympy.Symbol(
                                value + MV.construct_index(base))
                            symbol_lst.append(symbol)
                        self.mv[grade] = numpy.array(
                            symbol_lst, dtype=numpy.object)
                    else:
                        self.mv[0] = numpy.array(
                            [sympy.Symbol(value)], dtype=numpy.object)
                self.name = value + 'bm'
        if fct:
            if vars is not None:
                vars = tuple(vars)
            for grade in MV.n1rg:
                if not isinstance(self.mv[grade], int):
                    if grade == 0:
                        coef = sympy.galgebra.latex_ex.LatexPrinter.str_basic(
                            self.mv[0][0])
                        if vars is None and MV.coords is not None:
                            self.mv[0] = numpy.array([sympy.Function(
                                coef)(*MV.coords)], dtype=numpy.object)
                        else:
                            self.mv[0] = numpy.array([sympy.Function(
                                coef)(*vars)], dtype=numpy.object)
                    else:
                        for base in range(MV.nbasis[grade]):
                            coef = sympy.galgebra.latex_ex.LatexPrinter.str_basic(self.mv[grade][base])
                            if vars is None and MV.coords is not None:
                                self.mv[grade][
                                    base] = sympy.Function(coef)(*MV.coords)
                            else:
                                self.mv[
                                    grade][base] = sympy.Function(coef)(*vars)

    @staticmethod
    def construct_index(base):
        index_str = ''
        if len(base) == 0:
            return('')
        if MV.coords is None:
            for ix in base:
                index_str += str(ix + MV.index_offset)
        else:
            for ix in base:
                index_str += (MV.coords[ix]).name
        return('__' + index_str)

    def set_name(self, namestr):
        self.name = namestr
        return

    def max_grade(self):
        """
        X.max_grade() is maximum grade of non-zero grades of X.
        """
        for i in range(MV.n, -1, -1):
            if isinstance(self.mv[i], numpy.ndarray):
                return(i)
        return(-1)

    @staticmethod
    def coord(xname, offset=0):
        xi_str = ''
        for i in MV.nrg:
            xi_str += xname + str(i + offset) + ' '
        xi = sympy.symbols(xi_str)
        x = MV(xi, 'vector')
        return(x)

    def x(self, i):
        if isint(self.mv[1]):
            return(ZERO)
        return(self.mv[1][i])

    def set_coef(self, grade, base, value):
        if isinstance(self.mv[grade], int):
            self.mv[grade] = numpy.array(
                MV.nbasis[grade]*[ZERO], dtype=numpy.object)
        self.mv[grade][base] = value
        return

    @staticmethod
    def named(mvname, value='', mvtype=''):
        name = mvname
        tmp = MV(value=value, mvtype=mvtype, mvname=name)
        setattr(sys.modules[__name__], name, tmp)
        return

    @staticmethod
    def printnm(tpl):
        for a in tpl:
            print a.name, ' =', a.mv
        return

    def __str__(self):
        return(MV.str_rep(self))

    def printmv(self, name=''):
        title = ''
        if name:
            title += name + ' = '
        else:
            if self.name:
                title += self.name + ' = '
        print title + MV.str_rep(self)
        return

    def set_value(self, igrade, ibase, value):
        if isinstance(self.mv[igrade], numpy.ndarray):
            self.mv[igrade][ibase] = value
        else:
            self.mv[igrade] = numpy.array(
                MV.nbasis[igrade]*[ZERO], dtype=numpy.object)
            self.mv[igrade][ibase] = value
        return

    def add_in_place(self, mv):
        """
        X.add_in_place(mv) increments multivector X by multivector
        mv.
        """
        if self.bladeflg or mv.bladeflg:
            self.convert_to_blades()
            mv.convert_to_blades()
        for i in MV.n1rg:
            if isinstance(self.mv[i], numpy.ndarray) and isinstance(mv.mv[i], numpy.ndarray):
                self.mv[i] += mv.mv[i]
            else:
                if not isinstance(self.mv[i], numpy.ndarray) and isinstance(mv.mv[i], numpy.ndarray):
                    self.mv[i] = +mv.mv[i]
        return

    def sub_in_place(self, mv):
        """
        X.sub_in_place(mv) decrements multivector X by multivector
        mv.
        """
        if self.bladeflg or mv.bladeflg:
            self.convert_to_blades()
            mv.convert_to_blades()
        for i in MV.n1rg:
            if isinstance(self.mv[i], numpy.ndarray) and isinstance(mv.mv[i], numpy.ndarray):
                self.mv[i] -= mv.mv[i]
            else:
                if not isinstance(self.mv[i], numpy.ndarray) and isinstance(mv.mv[i], numpy.ndarray):
                    self.mv[i] = +mv.mv[i]
        return

    def __pos__(self):
        p = MV()
        p.bladeflg = self.bladeflg
        for i in MV.n1rg:
            if isinstance(self.mv[i], numpy.ndarray):
                p.mv[i] = +self.mv[i]
        return(p)

    def __neg__(self):
        n = MV()
        n.bladeflg = self.bladeflg
        for i in MV.n1rg:
            if isinstance(self.mv[i], numpy.ndarray):
                n.mv[i] = -self.mv[i]
        return(n)

    def __add_ab__(self, mv):
        self.add_in_place(mv)
        return

    def __sub_ab__(self, mv):
        self.sub_in_place(mv)
        return

    def __add__(self, mv):
        """See MV.addition(self,mv)"""
        return(MV.addition(self, mv))

    def __radd__(self, mv):
        """See MV.addition(mv,self)"""
        return(MV.addition(mv, self))

    def __sub__(self, mv):
        """See MV.subtraction(self,mv)"""
        return(MV.subtraction(self, mv))

    def __rsub__(self, mv):
        """See MV.subtraction(mv,self)"""
        return(MV.subtraction(mv, self))

    def __xor__(self, mv):
        """See MV.outer_product(self,mv)"""
        return(MV.outer_product(self, mv))

    def __pow__(self, mv):
        """See MV.outer_product(self,mv)"""
        return(MV.outer_product(self, mv))

    def __rxor__(self, mv):
        """See MV.outer_product(mv,self)"""
        return(MV.outer_product(mv, self))

    def __or__(self, mv):
        """See MV.inner_product(self,mv)"""
        return(MV.inner_product(self, mv))

    def __ror__(self, mv):
        """See MV.inner_product(mv,self)"""
        return(MV.inner_product(mv, self))

    def __lt__(self, mv):
        """See MV.inner_product(self,mv)"""
        return(MV.inner_product(self, mv, 'l'))

    def __lshift__(self, mv):
        """See MV.inner_product(self,mv)"""
        return(MV.inner_product(self, mv, 'l'))

    def __rlshift__(self, mv):
        """See MV.inner_product(self,mv)"""
        return(MV.inner_product(mv, self, 'l'))

    def lc(self, mv):
        return(MV.inner_product(self, mv, 'l'))

    def __gt__(self, mv):
        """See MV.inner_product(self,mv)"""
        return(MV.inner_product(self, mv, 'r'))

    def __rshift__(self, mv):
        """See MV.inner_product(self,mv)"""
        return(MV.inner_product(self, mv, 'r'))

    def __rrshift__(self, mv):
        """See MV.inner_product(self,mv)"""
        return(MV.inner_product(mv, self, 'r'))

    def rc(self, mv):
        return(MV.inner_product(self, mv, 'r'))

    def scalar_mul(self, c):
        """
        Y = X.scalar_mul(c), multiply multivector X by scalar c and return
        result.
        """
        mv = MV()
        mv.bladeflg = self.bladeflg
        mv.puregrade = self.puregrade
        for i in MV.n1rg:
            if isinstance(self.mv[i], numpy.ndarray):
                #print self.mv[i]
                #print c,type(c)
                mv.mv[i] = self.mv[i]*c
        return(mv)

    def scalar_mul_inplace(self, c):
        """
        X.scalar_mul_inplace(c), multiply multivector X by scalar c and save
        result in X.
        """
        for i in MV.n1rg:
            if isinstance(self.mv[i], numpy.ndarray):
                self.mv[i] = self.mv[i]*c
        return

    def __mul__(self, mv):
        """See MV.geometric_product(self,mv)"""
        return(MV.geometric_product(self, mv))

    def __rmul__(self, mv):
        """See MV.geometric_product(mv,self)"""
        return(MV.geometric_product(mv, self))

    def __div__(self, scalar):
        div = MV()
        div.bladeflg = self.bladeflg
        for i in MV.n1rg:
            if isinstance(self.mv[i], numpy.ndarray):
                div.mv[i] = self.mv[i]/scalar
        return(div)

    def __div_ab__(self, scalar):
        for i in MV.n1rg:
            if isinstance(self.mv[i], numpy.ndarray):
                self.mv[i] /= scalar
        return

    def __call__(self, igrade=0, ibase=0):
        """
        X(i,j) returns symbol in ith grade and jth base or blade of
        multivector X.
        """
        if not isinstance(self.mv[igrade], numpy.ndarray):
            return(ZERO)
        return(self.mv[igrade][ibase])

    @staticmethod
    def equal(mv1, mv2):
        mv1.compact()
        if isinstance(mv2, MV):
            mv2.compact()
        pure_grade = mv1.is_pure()
        if not isinstance(mv2, MV) and pure_grade != 0:
            return(False)
        if not isinstance(mv2, MV) and pure_grade == 0:
            if isinstance(mv1.mv[0], int):
                return(mv2 == 0)
            else:
                return(mv1.mv[0][0] == mv2)
        for (mvi, mvj) in zip(mv1.mv, mv2.mv):
            if isint(mvi) ^ isint(mvj):
                return(False)
            if isinstance(mvi, numpy.ndarray) and isinstance(mvj, numpy.ndarray):
                for (x, y) in zip(mvi, mvj):
                    if x != y:
                        return(False)
        return(True)

    def __eq__(self, mv):
        return(MV.equal(self, mv))

    def copy(self, sub=0):
        """
        Y = X.copy(), make a deep copy of multivector X in multivector
        Y so that Y can be modified without affecting X.
        """
        cpy = MV()
        cpy.name = self.name
        cpy.bladeflg = self.bladeflg
        cpy.puregrade = self.puregrade
        for i in MV.n1rg:
            if sub:
                if isinstance(self.mv[i], numpy.ndarray):
                    cpy.mv[i] = -self.mv[i]
            else:
                if isinstance(self.mv[i], numpy.ndarray):
                    cpy.mv[i] = +self.mv[i]
        return(cpy)

    def substitute_base(self, igrade, base, mv):
        if not isinstance(self.mv[igrade], numpy.ndarray):
            return
        if isinstance(base, numpy.ndarray):
            ibase = MV.basis[igrade].index(base)
        else:
            ibase = base
        coef = self.mv[igrade][ibase]
        if coef == ZERO:
            return
        self.mv[igrade][ibase] = ZERO
        self.add_in_place(mv*coef)
        return

    def convert_to_blades(self):
        """
        X.convert_to_blades(), inplace convert base representation
        to blade representation.  See reference 5 section 5.
        """
        if self.bladeflg:
            return
        self.bladeflg = 1
        for igrade in range(2, MV.n1):
            if isinstance(self.mv[igrade], numpy.ndarray):
                for ibase in range(MV.nbasis[igrade]):
                    coef = self.mv[igrade][ibase]
                    if (not coef == ZERO):
                        self.mv[igrade][ibase] = ZERO
                        self.add_in_place(MV.ibtable[igrade][ibase]*coef)
        return

    def convert_from_blades(self):
        """
        X.convert_from_blades(), inplace convert blade representation
        to base representation.  See reference 5 section 5.
        """
        if not self.bladeflg:
            return
        self.bladeflg = 0
        for igrade in range(2, MV.n1):
            if isinstance(self.mv[igrade], numpy.ndarray):
                for ibase in range(MV.nbasis[igrade]):
                    coef = self.mv[igrade][ibase]
                    if (not coef == ZERO):
                        self.mv[igrade][ibase] = ZERO
                        self.add_in_place(MV.btable[igrade][ibase]*coef)
        return

    def project(self, r):
        """
        Grade projection operator. For multivector X, X.project(r)
        returns multivector of grade r components of X if r is an
        integer. If r is a multivector X.project(r) returns a
        multivector consisting of the grade of X for which r has non-
        zero grades. For example if X is a general multivector and
        r is a general spinor then X.project(r) will return the even
        grades of X.
        """
        if isinstance(r, int):
            grade_r = MV()
            if r > MV.n:
                return(grade_r)
            self.convert_to_blades()
            if not isinstance(self.mv[r], numpy.ndarray):
                return(grade_r)
            grade_r.bladeflg = 1
            grade_r.puregrade = 1
            grade_r.mv[r] = +self.mv[r]
            return(grade_r)
        if isinstance(r, MV):
            self.convert_to_blades()
            r.convert_to_blades()
            proj = MV()
            for i in MV.n1rg:
                if not isinstance(r.mv[i], int):
                    proj.mv[i] = self.mv[i]
            proj.bladeflg = self.bladeflg
            return(proj)
        return(None)

    def even(self):
        """
        Even grade projection operator. For multivector X, X.even()
        returns multivector of even grade components of X.
        """
        egrades = MV()
        self.convert_to_blades()
        egrades.bladeflg = self.bladeflg
        egrades.puregrade = self.puregrade
        for igrade in range(0, MV.n1, 2):
            egrades.mv[igrade] = +self.mv[igrade]
        return(egrades)

    def odd(self):
        """
        Odd grade projection operator. For multivector X, X.odd()
        returns multivector of odd grade components of X.
        """
        ogrades = MV()
        self.convert_to_blades()
        ogrades.bladeflg = self.bladeflg
        ogrades.puregrade = self.puregrade
        for igrade in range(1, MV.n1, 2):
            ogrades.mv[igrade] = +self.mv[igrade]
        return(ogrades)

    def rev(self):
        """
        Revision operator. For multivector X, X.rev()
        returns reversed multivector of X.
        """
        revmv = MV()
        self.convert_to_blades()
        revmv.bladeflg = self.bladeflg
        revmv.puregrade = self.puregrade
        for igrade in MV.n1rg:
            if isinstance(self.mv[igrade], numpy.ndarray):
                if igrade < 2 or (not (((igrade*(igrade - 1))/2) % 2)):
                    revmv.mv[igrade] = +self.mv[igrade]
                else:
                    revmv.mv[igrade] = -self.mv[igrade]
        return(revmv)

    def cse(self, grade):
        cse_lst = []
        if isinstance(self.mv[grade], numpy.ndarray):
            for ibase in range(MV.nbasis[grade]):
                if self.mv[grade][ibase] != ZERO:
                    cse_lst.append(sympy.cse(self.mv[grade][ibase]))
        return(cse_lst)

    def div(self, grade, divisor):
        div_lst = []
        if isinstance(self.mv[grade], numpy.ndarray):
            for ibase in range(MV.nbasis[grade]):
                if self.mv[grade][ibase] != ZERO:
                    div_lst.append(
                        self.mv[grade][ibase].as_coefficient(divisor))
        return(div_lst)

    # don't know which one is correctly named

    def div(self):
        return self.grad_int()

    def collect(self, lst):
        """
        Applies sympy collect function to each component
        of multivector.
        """
        for igrade in MV.n1rg:
            if isinstance(self.mv[igrade], numpy.ndarray):
                for ibase in range(MV.nbasis[igrade]):
                    if self.mv[igrade][ibase] != ZERO:
                        self.mv[igrade][ibase] = \
                            sympy.collect(self.mv[igrade][ibase], lst)
        return

    def sqrfree(self, lst):
        """
        Applies sympy sqrfree function to each component
        of multivector.
        """
        for igrade in MV.n1rg:
            if isinstance(self.mv[igrade], numpy.ndarray):
                for ibase in range(MV.nbasis[igrade]):
                    if self.mv[igrade][ibase] != ZERO:
                        self.mv[igrade][ibase] = \
                            sympy.sqrfree(self.mv[igrade][ibase], lst)
        return

    def flatten(self):
        flst = []
        for igrade in MV.n1rg:
            if isinstance(self.mv[igrade], int):
                flst += MV.nbasis[igrade]*[ZERO]
            else:
                for coef in self.mv[igrade]:
                    flst.append(coef)
        return(flst)

    def subs(self, *args):
        X = MV()
        X.bladeflg = self.bladeflg
        X.puregrade = self.puregrade
        for igrade in MV.n1rg:
            if isinstance(self.mv[igrade], numpy.ndarray):
                X.mv[igrade] = numpy.array(
                    substitute_array(self.mv[igrade], *args))
        return(X)

    def sub_mv(self, mv1, mv2):
        mv1_flat = mv1.flatten()
        mv2_flat = mv2.flatten()
        self.sub_scalar(mv1_flat, mv2_flat)
        return

    def sub_scalar(self, expr1, expr2):
        if (isinstance(expr1, list) and isinstance(expr2, list)) or \
           (isinstance(expr1, tuple) and isinstance(expr2, tuple)):
            for (var1, var2) in zip(expr1, expr2):
                self.sub_scalar(var1, var2)
        else:
            for igrade in MV.n1rg:
                if not isinstance(self.mv[igrade], int):
                    for ibase in range(MV.nbasis[igrade]):
                        if expr1 != ZERO:
                            self.mv[igrade][ibase] = self.mv[
                                igrade][ibase].subs(expr1, expr2)
        return

    def simplify(self):
        """
        Applies sympy simplify function
        to each component of multivector.
        """
        for igrade in MV.n1rg:
            if isinstance(self.mv[igrade], numpy.ndarray):
                for ibase in range(MV.nbasis[igrade]):
                    if self.mv[igrade][ibase] != ZERO:
                        self.mv[igrade][ibase] = \
                            sympy.simplify(self.mv[igrade][ibase])
        return

    def trigsimp(self):
        """
        Applies sympy trigsimp function
        to each component of multivector
        using the recursive option.
        """
        for igrade in MV.n1rg:
            if isinstance(self.mv[igrade], numpy.ndarray):
                for ibase in range(MV.nbasis[igrade]):
                    if self.mv[igrade][ibase] != ZERO:
                        self.mv[igrade][ibase] = \
                            sympy.trigsimp(self.mv[igrade]
                                           [ibase], deep=True, recursive=True)
        return

    def cancel(self):
        """
        Applies sympy cancel function
        to each component of multivector.
        """
        for igrade in MV.n1rg:
            if isinstance(self.mv[igrade], numpy.ndarray):
                for ibase in range(MV.nbasis[igrade]):
                    if self.mv[igrade][ibase] != ZERO:
                        self.mv[igrade][ibase] = \
                            sympy.cancel(self.mv[igrade][ibase])
        return

    def expand(self):
        """
        Applies sympy expand function to each component
        of multivector.
        """
        for igrade in MV.n1rg:
            if isinstance(self.mv[igrade], numpy.ndarray):
                for ibase in range(MV.nbasis[igrade]):
                    if self.mv[igrade][ibase] != ZERO:
                        self.mv[igrade][ibase] = \
                            self.mv[igrade][ibase].expand()
        return

    def is_pure(self):
        igrade = -1
        ngrade = 0
        self.compact()
        for i in MV.n1rg:
            if isinstance(self.mv[i], numpy.ndarray):
                for base in self.mv[i]:
                    if base != 0:
                        igrade = i
                        ngrade += 1
                        break
                if ngrade > 1:
                    return(-1)
        if igrade == -1:
            return(0)
        return(igrade)

    def is_scalar(self):
        if self.is_pure() == 0:
            return(True)
        return(False)

    def compact(self):
        """
        Convert zero numpy arrays to single integer zero place holder
        in grade list for instantiated multivector. For example if
        numpy array of grade one components is a zero array then replace
        with single integer equal to zero.
        """
        for i in MV.n1rg:
            if not isint(self.mv[i]):
                zero_flg = True
                for x in self.mv[i]:
                    if x != 0:
                        zero_flg = False
                        break
                if zero_flg:
                    self.mv[i] = 0
        return

    def diff(self, x):
        """
        Calculate partial derivative of multivector with respect to
        argument x.
        """
        D = MV()
        igrade = 0
        for grade in self.mv:
            if not isinstance(grade, int):
                D.mv[igrade] = MV.vdiff(grade, x)
            igrade += 1
        return(D)

    def ddt(self):
        if MV.coords[0] != sympy.Symbol('t'):
            return(MV())
        dxdt = self.diff(MV.coords[0])
        for ibase in MV.nrg:
            dxdt += self.mv[1][ibase]*MV.dedt[ibase]
        #dxdt.simplify()
        #dxdt.trigsimp()
        return(dxdt)

    def grad(self):
        """
        Calculate geometric (grad) derivative of multivector function
        """
        D = []
        dD = MV()
        for theta in MV.coords:
            D.append(self.diff(theta))
        if MV.curvilinear_flg:
            recp = MV.Rframe
        else:
            recp = MV.brecp
        for (rbase, iD) in zip(recp, D):
            dD.add_in_place(rbase*iD)
        if MV.curvilinear_flg:  # Add Connection
            igrade = 1
            while igrade <= MV.n:
                coefs = self.mv[igrade]
                if type(coefs) != int:
                    for (coef, connect) in zip(coefs, MV.Connect[igrade]):
                        dD.add_in_place(coef*connect)
                igrade += 1
        return(dD)

    def grad_ext(self):
        """
        Calculate outer (exterior,curl) derivative of multivector function.
        """
        D = []
        dD = MV()
        for ix in MV.coords:
            D.append(self.diff(ix))
        if MV.curvilinear_flg:
            recp = MV.Rframe
        else:
            recp = MV.brecp
        for (irbase, iD) in zip(recp, D):
            dD.add_in_place(irbase ^ iD)
        if MV.curvilinear_flg:  # Add Connection
            igrade = 1
            while igrade <= MV.n:
                coefs = self.mv[igrade]
                if type(coefs) != int:
                    for (coef, connect) in zip(coefs, MV.Connect[igrade]):
                        if igrade < MV.n:
                            dD.add_in_place(coef*connect.project(igrade + 1))
                igrade += 1
        return(dD)

    def curl(self):
        return(self.grad_ext())

    def grad_int(self):
        """
        Calculate inner (interior,div) derivative of multivector function.
        """
        D = []
        dD = MV()
        for ix in MV.coords:
            D.append(self.diff(ix))
        if MV.curvilinear_flg:
            recp = MV.Rframe
        else:
            recp = MV.brecp
        for (irbase, iD) in zip(recp, D):
            dD.add_in_place(irbase | iD)
        if MV.curvilinear_flg:  # Add Connection
            igrade = 1
            while igrade <= MV.n:
                coefs = self.mv[igrade]
                if type(coefs) != int:
                    for (coef, connect) in zip(coefs, MV.Connect[igrade]):
                        dD.add_in_place(coef*connect.project(igrade - 1))
                igrade += 1
        return(dD)

    def mag2(self):
        """
        Calculate scalar component of square of multivector.
        """
        return((self | self)())

    def Name(self):
        """
        Get LaTeX name of multivector.
        """
        return(sympy.galgebra.latex_ex.LatexPrinter.extended_symbol(self.name))


def set_names(var_lst, var_str):
    """
    Set the names of a list of multivectors (var_lst) for a space delimited
    string (var_str) containing the names.
    """
    var_str_lst = var_str.split()
    if len(var_lst) == len(var_str_lst):
        for (var, var_name) in zip(var_lst, var_str_lst):
            var.set_name(var_name)
        return
    sys.stderr.write('Error in set_names. Lists incongruent!\n')
    sys.exit()
    return


def reciprocal_frame(vlst, names=''):
    """
    Calculate reciprocal frame of list (vlst) of vectors.  If desired name each
    vector in list of reciprocal vectors with names in space delimited string
    (names).
    """
    E = vlst[0]
    recp = []
    if type(names) != str:
        name_lst = names
    else:
        if names != '':
            name_lst = names.split()
    for i in range(1, MV.n):
        E = E ^ vlst[i]
    for i in range(MV.n):
        tmp = ONE
        if i % 2 != 0:
            tmp = -ONE
        for j in range(MV.n):
            if i != j:
                tmp = tmp ^ vlst[j]
        tmp = tmp*E
        recp.append(tmp)
    Esq = sympy.trigsimp(E.mag2(), deep=True, recursive=True)
    print Esq
    print sympy.simplify(Esq)
    Esq_inv = ONE/Esq
    i = 0
    for i in range(MV.n):
        recp[i].trigsimp()
        recp[i] = recp[i]*Esq_inv
        if names != '':
            recp[i].set_name(name_lst[i])
            i += 1
    return(recp)


def S(value):
    return(MV(value, 'scalar'))