#!/usr/bin/python

#galgebra.py

import os
import sys
import math
import copy
import re
import types

from string import *

import Numeric
import LinearAlgebra

float_type = type(1.0)
int_type   = type(1)
list_type  = type([])
str_type   = type('abc')

"""
Calculates n!
"""

def factorial(n):
    if n == 0 or n == 1: return(1)
    nfac = 1
    for i in range(1,n+1):
        nfac = nfac*i
    return(nfac)

"""
Calculates (m,n) = m!/(n!(m-n)!).
"""

def binomial(m,r):
    bcoef = factorial(m)/(factorial(r)*factorial(m-r))
    return(bcoef)

"""
Generates a list of the binomial coefficients (m,n),
where 0 <= n <= m.
"""

def binomial_array(m):
    mfac = factorial(m)
    barray = []
    for r in range(m+1):
        barray.append(mfac/(factorial(r)*factorial(m-r)))
    return(barray)

"""
Make a list of N entrys, each entry being the null string, ''.
"""

def make_empty_list(N):
    lst = []
    for i in range(N):
        lst.append('')
    return(lst)

"""
Takes a list of integers and removes duplicate values. Assumes
initial list is sorted.
"""

def remove_duplicates(lst):
    lst.sort()
    new_lst = []
    xold = ''
    for x in lst:
        if x != xold:
            new_lst.append(x)
        xold = x
    return(new_lst)

def union_int_list(lst1,lst2):
    lst = lst1+lst2
    return(remove_duplicates(lst))

def member_index(lst,value):
    index = lst.index(value)
    try:
        return(index)
    except ValueError:
        return(-1)

def insert_in_order(lst,x):
    if len(lst) == 0:
        lst.append(x)
        return((0,lst))
    i = 0
    for j in lst:
        if x < j:
            lst.insert(i,x)
            return((i,lst))
        i += 1
    lst.append(x)
    return((len(lst)-1,lst))

def inv_int_lst(lst,N):
    inv_lst = [-1 for i in range(N)]
    if len(lst) > 0:
        i = 0
        for j in lst:
            inv_lst[j] = i
            i += 1
    return(inv_lst)

"""
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 combinations(N,P):
    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)

"""
Multivectors are stored as a list of numeric arrays (default type is
float). If a grade of the multivector is empty (zero) the entry in the
list is a null string, ''. This is done to speed up operations on
multivectors (so that one would not waste time adding or multiplying
empty grades) and in realization that for very many calculations the
multivector is of a pure grade (vector,bivector,pseudo scalar,etc.).
The Geometric Algebra is initialized with the constructor:

     GeometicAlgebra(p,q)

Where p is the dimension of the positive signature basis and q is the
dimension of the negative signature basis.

In this representation of the Geometric Algebra a basis blade is denoted
by a list of integers. For example [1,2,3] represents the basis blade
which is the product of the basis vectors e1e2e3. Note that [] represents
the scalar 1.
"""

class GeometricAlgebra:
    def __init__(self,p,q=0,tol=1.0e-12):
        self.p = p
        self.q = q
        #Relative tolerance for setting component blade to zero
        self.tol = tol
        #Maximum grade in multivector space
        self.N = p+q 
        self.N1 = self.N+1
        #Dimension of multivector space
        self.dim = 2**self.N
        #List of dimensions of each grade subspace
        self.grade_dim_array = binomial_array(self.N)
        #Metric for base vector space
        self.metric = [1 for i in range(self.N)]
        for i in range(p,self.N):
            self.metric[i] = -1
        #Generate double indexed list (grade,blade) of all basis blades
        self.basis_blades = self.basis(self.N)            
        #Table to determine what grades are produced by multiplication of pure grades
        self.grade_mul_table()
        #Multiplication table for basis blades
        self.mul_index_table()

    """
    basis returns a list of lists containing all the basis blades of the
    geometric algebra. The return value base[r][s] returns a list of unit
    vector indices of the s'th component of the r'th grade basis blade.
    """

    def basis(self,N):
        base = []
        for n in range(N+1):
            base.append(combinations(N,n))
        return base
        
    """
    Multiply two multivector basis blades (blade1,blade2) together.
    The blades are represented as a list of integers that are the
    indices of each unit vector in the respective blade. The result
    is the basis blade that is the result of the product and the
    appropriate sign (+1 or -1). The indices of all basis blades are
    in ascending order.
    """

    def mul_blades(self,blade1,blade2):
        (sign,blade) = self.basis_product(blade1,blade2)
        if len(blade) == 1: return((sign,blade))
        i = 0
        n1 = len(blade)-1
        reduce_blade = []
        while i < n1:
            i1 = i+1
            if blade[i] == blade[i1]:
                if self.metric[blade[i]] == -1: sign = -sign
                i = i+2
            else:
                reduce_blade.append(blade[i])
                i = i+1
            if i == n1: reduce_blade.append(blade[-1])
        return((sign,reduce_blade))

    """
    grade_mul_table calculates what grade multivectors are generated when two pure
    grade multivectors are multiplied. The output is a double indexed list with
    integer list elements so that grade_mul[r][s] is a list of the grades in the
    product of a grade r and grade s multivector. This information is used when
    multivectors are multiplied to minimize the number of operatons.
    """

    def grade_mul_table(self):
        self.grade_mul = []
        for r in range(self.N1):
            self.grade_mul.append([])
            for s in range(self.N1):
                self.grade_mul[r].append([])
                for i in range(abs(r-s),min(r+s+1,2*self.N-(r+s)+1),2):
                    self.grade_mul[r][s].append(i)
        return

    """
    mul_index_table the multiplication table for all basis blades. The output 'mul_index' is
    a quadruply indexs list that contains the basis blade and correct sign for the product.
    The list is of form mul_index[g1][b1][g2][b2] = (g,b,sign) where (g1,b1) is the grade and
    blade index of the first basis blade,  (g2,b2) is the grade and blade index of the second
    basis blade, and (g,b,sign) is the grade and blade index and the sign of the product of the
    two basis blades.
    """

    def mul_index_table(self):
        self.mul_index = []
        for g1 in range(self.N1):
            self.mul_index.append([])
            for i1 in range(self.grade_dim_array[g1]):
                self.mul_index[g1].append([])
                for g2 in range(self.N1):
                    self.mul_index[g1][i1].append([])
                    for i2 in range(self.grade_dim_array[g2]):
                        blade1 = self.basis_blades[g1][i1]
                        blade2 = self.basis_blades[g2][i2]
                        (sign,blade) = self.mul_blades(blade1,blade2)
                        blade_grade = len(blade)
                        blade_index = self.basis_blades[blade_grade].index(blade)
                        if sign == 1:
                            sign_flg = 1
                        else:
                            sign_flg = 0
                        self.mul_index[g1][i1][g2].append((blade_grade,blade_index,sign_flg))
        return

    """
    basis_product calculates the product of two basis blades and
    returns the product as another basis blade with the correct
    sign for the product.
    """

    def basis_product(self,blade1,blade2):
        p = blade1+blade2
        sign = 1
        if not blade1: return((sign,blade2)) #Is blade1 a scalar
        if not blade2: return((sign,blade1)) #Is blade2 a scalar
        n1 = len(blade1)
        n2 = len(blade2)
        np = n1+n2
        for uv in range(n1,np):
            uv1 = uv
            while uv1 >= 1 and p[uv1] < p[uv1-1]:
                sign =  -sign
                tmp = p[uv1]
                p[uv1] = p[uv1-1]
                p[uv1-1] = tmp
                uv1 = uv1-1
        return((sign,p))


    """
    pretty_print prints the multiplication table for the Geometric
    Algebra. 'uvsymb' is the symbol used for the basis vector
    (default is 'e'). If 'filename' is given a TeX file of the
    multiplication table is generated so that 'uvsymb' can be
    any TeX symbol such as '\\gamma' for the greek gamma. The
    double backslash is required because of the way text strings
    are processed.
    """

    ############################Output Functions#################################

    def pretty_print(self,uvsymb='e',filename='',index_shift = 0):
        #Generate single indexed list of all basis blades
        self.basis_lst = []
        for blades in self.basis_blades:
            for blade in blades:
                self.basis_lst.append(blade)
        #Generate doulbe indexed list of basis blade products (mutiplication table)
        self.mul_table = []
        for blade1 in self.basis_lst:
            row = []
            for blade2 in self.basis_lst:
                row.append(self.mul_blades(blade1,blade2))
            self.mul_table.append(row)
        if filename:
            #Write TeX file for basis blades multiplication table
            TeX_file = open(filename+'.tex','w')
            TeX_str = '\\documentclass[8pt]{report}\n'+\
                      '\\usepackage[dvips]{graphicx,epsfig,rotating}\n'+\
                      '\\usepackage{amsmath}\n'+\
                      '\\begin{document}\n'+\
                      'Multiplication Table for $G('+str(self.p)+','+str(self.q)+')$\n'
            TeX_str = TeX_str+self.TeX_mul_table(uvsymb,index_shift)+'\\end{document}\n'
            TeX_file.write(TeX_str)
            TeX_file.close()
            return
        #Print pretty basis blades multiplication table with 'uvsymb' for basis vector symbol
        self.print_mul_table(uvsymb)
        return

    def TeX_mul_table(self,uvsymb,index_shift,max_cols=8):
        TeX_str = ''
        icol = 1
        jcol = min(max_cols,self.dim)
        while 1:
            TeX_str = TeX_str+self.TeX_table_block_str(uvsymb,icol,jcol,index_shift)
            if jcol == self.dim: break
            icol = jcol+1
            jcol = min(icol+max_cols-1,self.dim)
        return(TeX_str)

    def TeX_table_block_str(self,uvsymb,icol,jcol,index_shift):
        TeX_str = ''
        ncol = jcol-icol+1
        TeX_str = TeX_str+'\\begin{equation} \\hspace{-1in} \\nonumber '+\
        '\n\\begin{array}{c|*{'+str(ncol)+'}{c}} \n'
        for blade in self.basis_lst[icol-1:jcol]:
            TeX_str = TeX_str+' & '+self.TeX_table_str(blade,uvsymb,index_shift,0)
        TeX_str = TeX_str+' \\\\ \\hline \n'
        irow = 0
        for row in self.mul_table:
            TeX_str = TeX_str+self.TeX_table_str(self.basis_lst[irow],uvsymb,index_shift,0)+' '
            for blade in row[icol-1:jcol]:
                TeX_str = TeX_str+' & '+self.TeX_table_str(blade,uvsymb,index_shift)
            TeX_str = TeX_str+' \\\\ \n'
            irow += 1
        TeX_str = TeX_str+'\\end{array}\n\\end{equation}\n'
        return(TeX_str)

    def TeX_table_str(self,blade,uvsymb,index_shift,sgnflg=1):
        blade_str = ''
        if sgnflg:
            if blade[0] < 0: blade_str = blade_str+'-'
            if not blade[1]:
                blade_str = blade_str+'1'
            else:
                for uv in blade[1]:
                    if type(uvsymb) == types.StringType:
	                    blade_str = blade_str+uvsymb+'_{'+str(uv+index_shift)+'}'
                    else:
                        blade_str = blade_str+uvsymb[uv]
        else:
            if not blade:
                blade_str = blade_str+'1'
            else:
                for uv in blade:
                    if type(uvsymb) == types.StringType:
                        blade_str = blade_str+uvsymb+'_{'+str(uv+index_shift)+'}'
                    else:
                        blade_str = blade_str+uvsymb[uv]
        return(blade_str)

    def print_mul_table(self,uvsymb='e'):
        print 'Multiplication Table for Ga('+str(self.p)+','+str(self.q)+'):'
        print '\t\t'+self.blade_str(self.basis_lst,0)
        irow = 0
        for row in self.mul_table:
            print self.uv_product_str(self.basis_lst[irow],0)+'\t\t'+self.blade_str(row)
            irow += 1
        return

    """
    uv_product_str converts the list representation and a sign to a
    string representation for understandable printed output. The
    character uvshmb (default 'e') is used to represent unit vectors.
    """
    
    def uv_product_str(self,blade,sgnflg=1,uvsymb='e'):
        bstr = ''
        if sgnflg:
            if blade[0] < 0:
                bstr = bstr+'-'
            if not blade[1]:
                bstr = bstr+'1'
            else:
                for uv in blade[1]:
                    bstr = bstr+uvsymb+str(uv)
        else:
            if not blade:
                bstr = bstr+'1'
            else:
                for uv in blade:
                    bstr = bstr+uvsymb+str(uv)
        return(bstr)


    def blade_str(self,blades,sgnflg=1):
        bstr = ''
        if sgnflg:
            for blade in blades:
                if blade[0] < 0: bstr = bstr+'-'
                if not blade[1]:
                    bstr = bstr+'1'
                else: 
                    for uv in blade[1]:
                        bstr = bstr+'e'+str(uv)
                    bstr = bstr+'\t'
        else:
            for blade in blades:
                if not blade:
                    bstr = bstr+'1'
                else:
                    for uv in blade:
                        bstr = bstr+'e'+str(uv)
                    bstr = bstr+'\t'
        return(bstr)


    """
    Input parameters to constuct multivector:

        mvtype       grade             value
       
       'scalar'       NA               scalar
       'vector'       NA             scalar list
      'bivector'      NA             scalar list
      'pscalar'       NA               scalar
         ''           NA          scalar double list(1)
      'grades'    integer list(2) scalar double list
    'pure grade'    integer          scalar list

    (1) Use '' as outer list elements to denote empty grades.
    (2) List contains indices of non-empty grades
    """

class MultiVector:
    def __init__(self,value='',mvtype='',grades=-1):
        if mvtype == 'init' and value:
            p = value[0]
            if len(value) > 1:
                q = value[1]
            else:
                q = 0
            if len(value) > 2:
                tol = value[2]
            else:
                tol = 1.0e-12
            MultiVector.GA   = GeometricAlgebra(p,q,tol)
            MultiVector.N    = MultiVector.GA.N
            MultiVector.N1   = MultiVector.GA.N1
            MultiVector.tol  = MultiVector.GA.tol
            MultiVector.tol2 = MultiVector.tol**2
            MultiVector.I    = MultiVector(1.0,'pscalar')
            return
        self.mv  = []
        self.grades = []
        if mvtype == 'scalar':
            self.grades = [0]
            self.mv.append(Numeric.array([value],Numeric.Float))
        if mvtype == 'vector':
            self.grades = [1]
            self.mv.append(Numeric.array(value,Numeric.Float))
        if mvtype == 'zerovector':
            self.grades = [1]
            self.mv.append(Numeric.zeros(MultiVector.N,Numeric.Float))                    
        if mvtype == 'bivector':
            self.grades = [2]
            self.mv.append(Numeric.array(value,Numeric.Float))
        if mvtype == 'pscalar':
            self.grades = [MultiVector.N]
            self.mv.append(Numeric.array([value],Numeric.Float))
        if mvtype == 'grades':
            self.grades = grades
            i = 0
            for v in value:
                self.mv.append(Numeric.array(v,Numeric.Float))
        if mvtype == 'pure grade':
            if  grades >= 0:
                self.grades = [grades]
                self.mv.append(Numeric.array(value,Numeric.Float))
        if mvtype == 'zero':
            self.grades = grades
            for grade in self.grades:
                self.mv.append(Numeric.zeros(MultiVector.GA.grade_dim_array[grade],Numeric.Float))
        if mvtype == 'unitvector':
            self.grades = [1]
            self.mv.append(Numeric.zeros(MultiVector.N,Numeric.Float))
            self.mv[0][value] = 1.0
        self.inv_grades = inv_int_lst(self.grades,MultiVector.N1)
        return

    def scalar(self):
        if self.grades[0] != 0:
            return(0.0)
        else:
            return(self.mv[0][0])

    def dual(self):
        return(MultiVector.I*self)

    def get_grade(self,g):
        MV = MultiVector('','pure grade')
        gi = self.inv_grades[g]
        if gi >= 0:
            MV.mv.append(copy.deepcopy(self.mv[gi]))
            self.grades = [g]
            self.inv_grades = inv_int_lst(self.grades,MultiVector.N1)
        return(MV)

    def grade(self,g):
        gi = self.inv_grades[g]
        return(copy.deepcopy(self.mv[gi]))

    def grade_magnitudes2(self):
        if self.is_scalar():
            return(Numeric.array([self.mv[0][0]**2],Numeric.Float))
        if self.is_vector():
            grade_mag2 = 0.0
            for i in range(MultiVector.N):
                grade_mag2 += self.mv[0][i]**2
            return(Numeric.array([grade_mag2],Numeric.Float))
        grade_mag2 = Numeric.zeros(len(self.grades),Numeric.Float)
        for grade in self.grades:
            index = self.inv_grades[grade]
            for blade in range(MultiVector.GA.grade_dim_array[grade]):
                grade_mag2[index] += +self.mv[index][blade]**2
        return(grade_mag2)

    def sum_array(self,array):
        sum = 0.0
        for x in array:
            sum += x
        return(sum)
    
    def zero_residules(self,tol=1.0e-12):
        grade_mag2 = self.grade_magnitudes2()
        norm = self.sum_array(grade_mag2)
        if norm == 0.0:
            self.grades = [0]            
            self.inv_grades = inv_int_lst(self.grades,MultiVector.N1)
            self.mv = [Numeric.zeros([1],Numeric.Float)]
            return
        norm = tol*norm
        grades_rev = copy.deepcopy(self.grades)
        grades_rev.reverse()
        for grade in grades_rev:
            index = self.inv_grades[grade]
            if grade_mag2[index] < norm:
                self.mv.pop(index)
                self.grades.remove(grade)
        self.inv_grades = [-1 for i in range(MultiVector.N1)]
        i = 0
        for grade in self.grades:
            self.inv_grades[grade] = i
            i += 1
        return

    def is_scalar(self):
        if len(self.grades) > 1:
            return(0)
        if self.grades[0] == 0:
            return(1)
        return(0)

    def is_vector(self):
        if len(self.grades) > 1:
            return(0)
        if self.grades[0] == 1:
            return(1)
        return(0)

    def is_bivector(self):
        if len(self.grades) > 1:
            return(0)
        if self.grades[0] == 2:
            return(1)
        return(0)

    def is_pscalar(self):
        if len(self.grades) > 1:
            return(0)
        if self.grades[0] == MultiVector.N:
            return(1)
        return(0)    

    def __call__(self,grade=0,blade=0):
        index = self.inv_grades[grade]
        if index >= 0:
            return(self.mv[index][blade])
        else:
            return('')

    def __neg__(self):
        MV = copy.deepcopy(self)
        for grade in MV.grades:
            index = MV.inv_grades[grade]
            MV.mv[index] = -MV.mv[index]
        return(MV)

    def __pos__(self):
        return(self)

    def inverse(self):
        selfsq = self*self
        if selfsq.mv[0][0] != 0.0:
            return((1.0/selfsq.mv[0][0])*self)
        else:
            return(0)

    def __invert__(self):
        return(self.reverse())
    
    def __add__(self,x):
        if float_type == type(x) or int_type == type(x):
            MV = copy.deepcopy(self)
            if self.grades[0] != 0:
                MV.grades.insert(0,0)
                MV.inv_grades = inv_int_lst(MV.grades,MultiVector.N1)
                MV.mv.insert(0,Numeric.array(x,Numeric.Float))
            else:
                MV.mv[0] = MV.mv[0]+Numeric.array([x],Numeric.Float)
            return(MV)
        MV = MultiVector('','pure grade')
        MV.grades = union_int_list(self.grades,x.grades)
        MV.inv_grades = inv_int_lst(MV.grades,MultiVector.N1)
        for g in MV.grades:
            g1i = self.inv_grades[g] 
            g2i = x.inv_grades[g]
            if g1i >= 0 and g2i >= 0:
                MV.mv.append(self.mv[g1i]+x.mv[g2i])
            if g1i >= 0 and g2i < 0:
                MV.mv.append(copy.deepcopy(self.mv[g1i]))
            if g1i < 0 and g2i >= 0:
                MV.mv.append(copy.deepcopy(x.mv[g2i]))                
        return(MV)

    def __radd__(self,x):
        return(self+x)
 
    def __sub__(self,x):
        if float_type == type(x) or int_type == type(x):
            MV = copy.deepcopy(self)
            if self.grades[0] != 0:
                MV.grades.insert(0,0)
                MV.inv_grades = inv_int_lst(MV.grades,MultiVector.N1)
                MV.mv.insert(0,Numeric.array([x],Numeric.Float))
            else:
                MV.mv[0] = MV.mv[0]+Numeric.array([x],Numeric.Float)
            return(MV)
        MV = MultiVector('','pure grade')
        MV.grades = union_int_list(self.grades,x.grades)
        MV.inv_grades = inv_int_lst(MV.grades,MultiVector.N1)
        for g in MV.grades:
            g1i = self.inv_grades[g] 
            g2i = x.inv_grades[g]
            if g1i >= 0 and g2i >= 0:
                MV.mv.append(self.mv[g1i]-x.mv[g2i])
            if g1i >= 0 and g2i < 0:
                MV.mv.append(-copy.deepcopy(self.mv[g1i]))
            if g1i < 0 and g2i >= 0:
                MV.mv.append(-copy.deepcopy(x.mv[g2i]))                
        return(MV)

    def __rsub__(self,x):
        if float_type == type(x) or int_type == type(x):
            MV = copy.deepcopy(self)
            if self.grades[0] != 0:
                MV.grades.insert(0,0)
                MV.inv_grades = inv_int_lst(MV.grades,MultiVector.N1)
                MV.mv.insert(0,Numeric.array([x],Numeric.Float))
            else:
                MV.mv[0] = MV.mv[0]-Numeric.array([x],Numeric.Float)
            for i in range(1,len(MV.grades)):
                MV.mv[i] = -MV.mv[i]
            return(MV)
        MV = MultiVector('','pure grade')
        MV.grades = union_int_list(self.grades,x.grades)
        MV.inv_grades = inv_int_lst(MV.grades,MultiVector.N1)
        for g in MV.grades:
            g1i = self.inv_grades[g] 
            g2i = x.inv_grades[g]
            if g1i >= 0 and g2i >= 0:
                MV.mv.append(x.mv[g2i]-self.mv[g1i])
            if g1i >= 0 and g2i < 0:
                MV.mv.append(copy.deepcopy(x.mv[g1i]))
            if g1i < 0 and g2i >= 0:
                MV.mv.append(-copy.deepcopy(self.mv[g2i]))                
        return(MV)        

    def __mul__(self,x):

        if type(x) == list_type:
            lst = []
            for mv in x:
                lst.append(self*mv)
            return(lst)
        
        #scalar multiplication

        if float_type == type(x) or int_type == type(x):
            MV = MultiVector('','pure grade')
            MV.grades = self.grades
            MV.inv_grades = self.inv_grades
            for array in self.mv:
                    MV.mv.append(x*array)
            return(MV)

        #geometric product

        """
        Determine what grades are generated when grade1 multiplies grade2 and
        generate the properly dimensioned index and zero filled numerical arrays to
        accumulate the products of the grade components (blades).
        """
        grades = []
        for grade1 in self.grades:
            for grade2 in x.grades:
                new_grades = MultiVector.GA.grade_mul[grade1][grade2]
                grades = grades+new_grades
        grades.sort()
        grades = remove_duplicates(grades)
        MV = MultiVector('','zero',grades)
 
        for grade1 in self.grades:
            index1 = self.inv_grades[grade1]
            for grade2 in x.grades:
                index2 = x.inv_grades[grade2]
                """
                Determine grade and blade index and sign for blade products. Multiply blades and
                accumulate product according to proper sign (add or subtract from accumulated
                component).
                """
                for blade1 in range(MultiVector.GA.grade_dim_array[grade1]):
                    for blade2 in range(MultiVector.GA.grade_dim_array[grade2]):
                        b1b2 = self.mv[index1][blade1]*x.mv[index2][blade2]
                        (grade,blade,sign_flg) = MultiVector.GA.mul_index[grade1][blade1][grade2][blade2]
                        index = MV.inv_grades[grade]
                        if sign_flg:
                            MV.mv[index][blade] = MV.mv[index][blade]+b1b2
                        else:
                            MV.mv[index][blade] = MV.mv[index][blade]-b1b2
        MV.zero_residules()
        return(MV)    

    def __rmul__(self,x):
        if float_type == type(x) or int_type == type(x):
            return(self*x)
        return('')


    def __xor__(self,MV):
        return(wedge([self,MV]))

    def __or__(self,MV):
        return(dot(self,MV))
    
##################Printing Functions for Multivectors##############################

    def printmv(self,name='',uvsymb='e',mode=1):
        mv_str = ''
        if name:
            mv_str = mv_str+name+' = '
        if self.is_scalar():
            x = self.mv[0][0]
            if x != 0.0:
                mv_str = mv_str+str(self.mv[0][0])
            else:
                mv_str = mv_str+'0.0'
            if mode:
                print mv_str
                return
            else:
                print(mv_str)
            
        first_flg = 1
        for grade in self.grades:
            gi = self.inv_grades[grade]
            for blade in range(MultiVector.GA.grade_dim_array[grade]):
                xgb = self.mv[gi][blade]
                base = MultiVector.GA.basis_blades[grade][blade]
                if xgb > 0.0:
                    if first_flg:
                        first_flg = 0
                        mv_str = mv_str+str(xgb)
                    else:
                        mv_str = mv_str+'+'+str(xgb)
                if xgb < 0.0:
                    if first_flg:
                        first_flg = 0
                    mv_str = mv_str+str(xgb)
                if xgb != 0.0:
                    mv_str = mv_str+self.blade_str(base,uvsymb)
        if mode:
            print mv_str
            return
        return(mv_str)

    def blade_str(self,blade,uvsymb):
        bstr = ''
        if len(blade) == 0:
            return(bstr)
        if blade:
            bstr = bstr+'*'
            for uv in blade:
                if type(uvsymb) == list_type:
                    bstr = bstr+uvsymb[uv]
                else:
                    bstr = bstr+uvsymb+str(uv)
        return(bstr)

def is_pure_grade(mv_lst):
    for mv in mv_lst:
        if len(mv.grades) > 1:
            return(0)
    return(1)

def reverse(MV):
    MV_rev = copy.deepcopy(MV)
    i = 0
    for grade in MV.grades:
        if ((grade*(grade-1))/2)%2 == 1:
            MV_rev.mv[i] = -MV_rev.mv[i]
        i += 1
    return(MV_rev)

def dot(MV1,MV2):
    if MV1.is_scalar() or MV2.is_scalar():
        return(0.0)
    if MV1.is_vector() and MV2.is_vector():
        MV1sq = (MV1.grade_magnitudes2())[0]
        MV2sq = (MV2.grade_magnitudes2())[0]
        v1 = MV1.mv[0]
        v2 = MV2.mv[0]
        dot_prod = 0.0
        for i in range(MultiVector.N):
            if MultiVector.GA.metric[i] > 0:
                dot_prod += v1[i]*v2[i]
            else:
                dot_prod -= v1[i]*v2[i]
        if dot_prod**2 < MV1sq*MV2sq*MultiVector.tol2:
            dot_prod = 0.0
        return(dot_prod)
        #return(MultiVector(dot_prod,'scalar'))
    dot_grade = abs(MV1.grades[0]-MV2.grades[0])
    """
    Determine what grades are generated when grade1 multiplies grade2 and
    generate the properly dimensioned index and zero filled numerical arrays to
    accumulate the products of the grade components (blades).
    """
    MV = MultiVector('','zero',[dot_grade])
    for grade1 in MV1.grades:
        index1 = MV1.inv_grades[grade1]
        for grade2 in MV2.grades:
            index2 = MV2.inv_grades[grade2]
            """
            Determine grade and blade index and sign for blade products. Multiply blades and
            accumulate product according to proper sign (add or subtract from accumulated
            component).
            """
            for blade1 in range(MultiVector.GA.grade_dim_array[grade1]):
                for blade2 in range(MultiVector.GA.grade_dim_array[grade2]):
                    (grade,blade,sign_flg) = MultiVector.GA.mul_index[grade1][blade1][grade2][blade2]
                    if grade == dot_grade:
                        b1b2 = MV1.mv[index1][blade1]*MV2.mv[index2][blade2]
                        if sign_flg:
                            MV.mv[0][blade] = MV.mv[0][blade]+b1b2
                        else:
                            MV.mv[0][blade] = MV.mv[0][blade]-b1b2
    MV.zero_residules()
    if dot_grade == 0:
        return(MV.mv[0][0])
    else:
        return(MV)

"""
wedge: calculate the outer product of a list of pure grade multivectors.
"""

def wedge(mv_lst):
    if is_pure_grade(mv_lst):
        if len(mv_lst) == 2:
            MV1 = mv_lst[0]
            MV2 = mv_lst[1]
            wedge_grade = MV1.grades[0]+MV2.grades[0]
            if wedge_grade > MultiVector.N:
                return(MultiVector(0,'scalar'))
  
            #outer product
            """
            Determine what grades are generated when grade1 multiplies grade2 and
            generate the properly dimensioned index and zero filled numerical arrays to
            accumulate the products of the grade components (blades).
            """
            MV = MultiVector('','zero',[wedge_grade])
            for grade1 in MV1.grades:
                index1 = MV1.inv_grades[grade1]
                for grade2 in MV2.grades:
                     index2 = MV2.inv_grades[grade2]
                     """
                     Determine grade and blade index and sign for blade products. Multiply blades and
                     accumulate product according to proper sign (add or subtract from accumulated
                     component).
                     """
                     for blade1 in range(MultiVector.GA.grade_dim_array[grade1]):
                         for blade2 in range(MultiVector.GA.grade_dim_array[grade2]):
                             b1b2 = MV1.mv[index1][blade1]*MV2.mv[index2][blade2]
                             (grade,blade,sign_flg) = MultiVector.GA.mul_index[grade1][blade1][grade2][blade2]
                             if grade == wedge_grade:
                                 if sign_flg:
                                     MV.mv[0][blade] = MV.mv[0][blade]+b1b2
                                 else:
                                     MV.mv[0][blade] = MV.mv[0][blade]-b1b2
        else:
            MV = mv_lst[0]
            for mv in mv_lst[1:]:
                MV = wedge([MV,mv])
        MV.zero_residules()
        return(MV)
    else:
        return(0.0)

def commutator(A,B):
    return(0.5*(A*B-B*A))

def rotor(v1,v2,theta):
    B = wedge([v1,v2])
    Bsq = (B*B)(0)
    B = (1.0/Numeric.sqrt(abs(Bsq)))*B
    if Bsq < 0.0:
        U = Numeric.cos(theta/2)+Numeric.sin(theta/2)*B
    else:
        U = Numeric.cosh(theta/2)+Numeric.sinh(theta/2)*B
    return((U,reverse(U)))

def bexp(B):
    if not B.is_bivector():
        return(0.0)
    Bsq = B*B()
    Bmag = sqrt(abs(Bsq))
    if Bsq > 0.0:
       Bnorm = (1.0/Bmag)*B
       return(cosh(Bmag)+Bnorm*sinh(Bmag))        
    if Bsq < 0.0:
        Bnorm = (1.0/Bmag)*B
        return(cos(Bmag)+Bnorm*sin(Bmag))
    return(1.0+B)

def reciprocal_frame(frame):
    E = wedge(frame)
    Einv = E.inverse()
    rec_frame = []
    for i in range(len(frame)):
        fi = []
        for j in range(len(frame)):
            if j != i:
                fi.append(frame[j])
        if i%2 == 0:
            rec_frame.append(wedge(fi)*Einv)
        else:
            rec_frame.append(-wedge(fi)*Einv)
    return(rec_frame)

def dot_vector_array(varr1,varr2):
    dp = []
    i = 0
    for v1 in varr1:
        dp.append([])
        for v2 in varr2:
            dp[i].append((dot(v1,v2))())
        i += 1
    return(dp)

def printmv(mv_lst,Name='',uvsymb='e',mode=1):
    mv_str = ''
    if Name:
        mv_str += Name+' = ['
    else:
        mv_str += '['
    for MV in mv_lst[:-1]:
        mv_str += MV.printmv('',uvsymb,0)+',\n'
    mv_str+= mv_lst[-1].printmv('',uvsymb,0)+']'
    if mode:
        print mv_str
        return
    return(mv_str)

class LinearTransformation:
    GA = ''
    bases = ''
    def __init__(self,F='',lttype=''):
        if lttype == 'init':
            LinearTransformation.GA = F
            LinearTransformation.bases = self.generate_bases()
            return
        self.lttype = lttype
        self.F = F
        if lttype=='matrix':
            self.F = ''
            self.Fmv_matrix = [F]
            self.col_bases = []
            for icol in range(LinearTransformation.GA.N):
                MV = MultiVector('','zero',[1])
                MV.mv[0] = self.Fmv_matrix[0][:,icol]
                self.col_bases.append(MV)
            self.Build_Fmv_matrix()
            return
        if lttype=='function':
            self.F = F
            self.col_bases = []
            col_arrays = []
            N = LinearTransformation.GA.N
            for icol in range(LinearTransformation.GA.N):
                MV = self.F(LinearTransformation.bases[0][icol])
                col_arrays.append(MV.mv[0])
                self.col_bases.append(MV)
            Fmatrix = Numeric.concatenate(col_arrays)
            Fmatrix = Numeric.resize(self.Fmatrix,(N,N))
            Fmatrix = Numeric.swapaxes(self.Fmatrix,0,1)
            self.Fmv_matrix = [Fmatrix]
            self.Build_Fmv_matrix()
            return
        return

    def Build_Fmv_matrix(self):
        basis = LinearTransformation.GA.basis_blades
        N1 = LinearTransformation.GA.N1
        for grade in range(2,N1):
            blades = basis[grade]
            Nb = len(blades)
            array_lst = []
            for blade in blades:
                v_lst = []
                for index in blade:
                    v_lst.append(self.col_bases[index])
                array_lst.append(wedge(v_lst).mv[0])
            array = Numeric.concatenate(array_lst)
            array = Numeric.resize(array,(Nb,Nb))
            array = Numeric.swapaxes(array,0,1)
            self.Fmv_matrix.append(array)
        return

    def generate_bases(self):
        bases = []
        N1 = LinearTransformation.GA.N1
        print 
        for grade in range(1,N1):
            bases.append([])
            nblades = LinearTransformation.GA.grade_dim_array[grade]
            for blade in range(nblades):
                array = Numeric.zeros(nblades,Numeric.Float)
                array[blade] = 1.0
                MV = MultiVector()
                MV.mv = []
                MV.mv.append(array)
                MV.grades = [grade]
                MV.inv_grades = inv_int_lst(MV.grades,N1)
                bases[grade-1].append(MV)
        return(bases)

    def __add__(self,x):
        matrix = self.matrix()+x.matrix()
        return(LinearTransformation(matrix,'matrix'))

    def __radd__(self,x):
        return(self+x)

    def __mul__(self,x):
        if type(x) == float_type or type(x) == int_type:
            matrix = (1.0*x)*self.matrix()
            return(LinearTransformation(matrix,'matrix'))
        matrix = LinearAlgebra.matrixmultiply(self.matrix(),x.matrix())
        return(LinearTransformation(matrix,'matrix'))

    def __rmul__(self,x):
        if type(x) == float_type or type(x) == int_type:
            return(self*x)
        return        

    def print_array_mv(self,mv_array,uvsymb='e',mode=1):
        array_mv_str = '['
        for grade in mv_array:
            if type(grade) == list_type:
                array_mv_str += '['
                for mv in grade:
                    if mv == grade[-1]:
                        array_mv_str += mv.printmv('',uvsymb,0)
                    else:
                        array_mv_str += mv.printmv('',uvsymb,0)+','
                if grade == mv_array[-1]:
                    array_mv_str += ']'
                else:
                    array_mv_str += ']\n'
            else:
                if grade == mv_array[-1]:
                    array_mv_str += grade.printmv('',uvsymb,0)
                else:
                    array_mv_str += grade.printmv('',uvsymb,0)+',\n'                
        array_mv_str += ']'
        if mode:
            print array_mv_str
            return
        return(array_mv_str)

    def adjoint(self):
        Fadj = Numeric.swapaxes(self.Fmv_matrix[0],0,1)
        return(LinearTransformation(Fadj,'matrix'))

    def matrix(self):
        return(self.Fmv_matrix[0])

    def det(self):
        return(self.Fmv_matrix[MultiVector.N][0][0])

    def __call__(self,MV):
        if type(MV) == list_type:
            lst = []
            for mv in MV:
                lst.append(self(mv))
            return(lst)
        if MV.is_vector():
            LTMV = MultiVector('','zero',[1])
            LTMV.mv[0] = Numeric.matrixmultiply(self.Fmv_matrix[0],MV.mv[0])
            return(LTMV)
        if MV.is_scalar():
            return(MV)
        LTMV = MultiVector()
        for grade in MV.grades:
            LTMV.grades.append(grade)
            matrix = self.Fmv_matrix[grade-1]
            v      = MV.mv[MV.inv_grades[grade]]
            LTMV.mv.append(Numeric.matrixmultiply(matrix,v))
        LTMV.inv_grades = inv_int_lst(LTMV.grades,MV.N1)
        return(LTMV)

def unit_vectors(GA):
    N = GA.N
    uv = []
    index = [ 0 for i in range(N)]
    for i in range(N):
        if i > 0:
            index[i-1] = 0
        index[i] = 1
        uv.append(MultiVector(index,'vector'))
    return(uv)