Ondrej Certik wrote:
On Tue, Jun 8, 2010 at 6:28 AM, Andy Ray Terrel <[email protected]> wrote:
Hello Alan,
Do you have anymore details of how the tensor module will work? Way
back when I was coding more with Sympy we tried to do something based
on GiNaC's Indexed system. It would be good to have a more uniform
system such that vectors, matrices, and tensors were implemented in a
very similar manner.
Øyvind is working on code generation for quantum chemistry. This
requires tensors so we are definitely interested in your work.
Yes, I am interested as well. Here is ho ginac does it:
http://www.ginac.de/tutorial/Indexed-objects.html#Indexed-objects
http://www.ginac.de/tutorial/Non_002dcommutative-objects.html#Non_002dcommutative-objects
Ondrej
Ondrej
Try attached not ready for prime time code. Once I clean it up and
generate enough documentation (internal and external) to make it useful
I will submit it.
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To post to this group, send email to [email protected].
To unsubscribe from this group, send email to
[email protected].
For more options, visit this group at
http://groups.google.com/group/sympy?hl=en.
#!/bin/env python
import math,copy,sys,pyx,os,subprocess
from sympy import latex,Symbol,simplify,Function,diff,Rational
import numpy
global MAIN_PROGRAM
MAIN_PROGRAM = ''
ZERO = Rational(0,1)
ONE = Rational(1,1)
HALF = Rational(1,2)
even = True
odd = False
old = True
new = False
def set_main(main_program):
global MAIN_PROGRAM
MAIN_PROGRAM = main_program
return
set_main(sys.modules[__name__])
def perm_lst(str):
if len(str) <=1:
yield str
else:
for perm in perm_lst(str[1:]):
for i in range(len(perm)+1):
yield perm[:i] + str[0:1] + perm[i:]
def even_perm(perm):
"""Determine whether a permutation is even.
An even permutation is a permutation that can be produced by an even number
of exchanges.
Positional arguments :
perm -- permutation (tuple)
"""
dim = len(perm)
even = True
for i in xrange(dim):
for j in xrange(1 + i, dim):
if perm[j] < perm[i]:
even = not even
return even
def find_repeats(lst):
rlst = []
for i in range(len(lst)-1):
lst_i = lst[i]
if lst_i in lst[i+1:]:
if lst_i not in rlst:
rlst.append(lst_i)
return(rlst)
def separate_free_and_repeats(lst):
free_lst = []
repeat_lst = []
for i in range(len(lst)-1):
lst_i = lst[i]
if lst_i in lst[i+1:]:
if lst_i not in repeat_lst:
repeat_lst.append(lst_i)
else:
if lst_i not in repeat_lst:
free_lst.append(lst_i)
if lst[-1] not in repeat_lst:
free_lst.append(lst[-1])
return(free_lst,repeat_lst)
def cycle_list(lst,shift):
n = len(lst)
new_lst = n*[0]
for i in range(n):
new_lst[i] = lst[(i-shift)%n]
return(new_lst)
def compress_free_list(lst):
sort_lst = copy.copy(lst)
sort_lst.sort()
n_lst = len(lst)
comp_lst = n_lst*[0]
for i in range(n_lst):
comp_lst[lst.index(sort_lst[i])] = i+1
return(comp_lst)
def canonical_index(lst):
(free_lst,repeat_lst) = separate_free_and_repeats(lst)
reduced_lst = free_lst+repeat_lst
reduced_lst.sort()
can_lst = len(lst)*[0]
i = 1
for index in reduced_lst:
new_index = lst.index(index)
can_lst[new_index] = i
if index in repeat_lst:
new_index = lst.index(index,new_index+1)
can_lst[new_index] = i
i += 1
return(can_lst)
def minus(lst):
mlst = []
for x in lst:
mlst.append(-x)
return(mlst)
def make_tensor(tensorlst):
"""
make_tensor takes a string of tensor names separated by
blanks and converts them tensors separately
accessible by the main program and addition returns a list
of the symbols.
"""
global MAIN_PROGRAM
tensorlst = tensorlst.split()
tensor_lst = []
isym = 0
for name in tensorlst:
tmp = Tensor(name)
tensor_lst.append(tmp)
setattr(MAIN_PROGRAM,name,tmp)
isym += 1
return(tensor_lst)
class PDF:
mode = 'Pyx'
preamble = '\\documentclass[10pt,letter,fleqn]{report}\n'+\
'\\pagestyle{empty}\n'+\
'\\usepackage[latin1]{inputenc}\n'+\
'\\usepackage[pdftex,landscape,top=1cm,nohead,nofoot]{geometry}\n'+\
'\\usepackage{amsmath}\n'+\
'\\usepackage{bm}\n'+\
'\\usepackage{amsfonts}\n'+\
'\\usepackage{amssymb}\n'+\
'\\usepackage{tensor}\n'+\
'\\setlength{\\parindent}{0pt}\n'+\
'\\newcommand{\\bfrac}[2]{\\displaystyle\\frac{#1}{#2}}\n'+\
'\\newcommand{\\lp}{\\left (}\n'+\
'\\newcommand{\\rp}{\\right )}\n'+\
'\\newcommand{\\half}{\\frac{1}{2}}\n'+\
'\\newcommand{\\llt}{\\left <}\n'+\
'\\newcommand{\\rgt}{\\right >}\n'+\
'\\newcommand{\\abs}[1]{\\left |{#1}\\right | }\n'+\
'\\newcommand{\\pdiff}[2]{\\bfrac{\\partial {#1}}{\\partial {#2}}}\n'+\
'\\newcommand{\\lbrc}{\\left \\{}\n'+\
'\\newcommand{\\rbrc}{\\right \\}}\n'+\
'\\newcommand{\\W}{\\wedge}\n'+\
"\\newcommand{\\prm}[1]{{#1}'}\n"+\
'\\newcommand{\\ddt}[1]{\\bfrac{d{#1}}{dt}}\n'+\
'\\newcommand{\\R}{\\dagger}\n'+\
'\\newcommand{\\be}{\\begin{equation*}}\n'+\
'\\newcommand{\\ee}{\\end{equation*}}\n'+\
'\\begin{document}\n'
postscript = '\\end{document}\n'
@staticmethod
def setup(mode=''):
PDF.mode = mode
if PDF.mode == '':
return
if PDF.mode == 'Pyx':
PDF.filename = 'trash'
PDF.old_filename = 'old_trash'
PDF.file_flg = True
pyx.text.set(mode='latex')
pyx.text.preamble(r'\usepackage{tensor}')
pyx.text.preamble(r'\usepackage{amsmath}')
pyx.text.preamble(r'\usepackage{amsfonts}')
pyx.text.preamble(r'\usepackage{amssymb}')
PDF.canvas = pyx.canvas.canvas()
PDF.vpos = 0.0
PDF.dvpos = 0.5
PDF.file = ''
if PDF.mode == 'latex':
PDF.latex = []
return
@staticmethod
def write(latex_str,inline=True):
if inline:
latex_str = '$'+latex_str+'$'
if PDF.mode == 'Pyx':
PDF.canvas.text(0,PDF.vpos,latex_str)
PDF.vpos -= PDF.dvpos
filename = ''
remove_filename = ''
if PDF.file_flg:
filename = PDF.filename
remove_filename = PDF.old_filename
else:
filename = PDF.old_filename
remove_filename = PDF.filename
PDF.file_flg = not PDF.file_flg
PDF.canvas.writePDFfile(filename)
subprocess.Popen('xpdf -g 500x750 -cont -z 150 -remote myserver '+filename+'.pdf',shell=True)
os.system('rm '+remove_filename+'.pdf')
if PDF.mode == 'latex':
PDF.latex.append(latex_str)
return
@staticmethod
def display(filename='pdflatex'):
if PDF.mode == 'latex':
latex_str = PDF.preamble
for line in PDF.latex:
line = line.replace('$','')
latex_str += '\\be\n\\hspace{-1.5in}'+line+'\n\\ee\n'
latex_str += PDF.postscript
latex_file = open(filename+'.tex','w')
latex_file.write(latex_str)
latex_file.close()
os.system('pdflatex '+filename+'.tex')
subprocess.Popen('xpdf -g 500x750 -cont -z 150 -remote myserver '+filename+'.pdf',shell=True)
return
class BasicTensor:
@staticmethod
def cat_shift(lst1,lst2):
lst = copy.copy(lst1)
n = len(lst1)
for i in lst2:
lst.append(i+n)
return(lst)
@staticmethod
def shift(n,lst1):
lst = []
for i in lst1:
lst.append(n+i)
return(lst)
@staticmethod
def minus(lst):
mlst = []
for x in lst:
mlst.append(-x)
return(mlst)
@staticmethod
def scale(sf,lst):
sflst = []
for x in lst:
sflst.append(sf*x)
return(sflst)
@staticmethod
def setup(pdfmode='latex',dim=2,offset=0):
BasicTensor.index_pool = 'abcdefhijklmnopquvwxyz'
BasicTensor.dim = dim
BasicTensor.offset = offset
BasicTensor.Pd = BasicTensor('Pd_a')
BasicTensor.Pd.op_flg = True
BasicTensor.Cd = BasicTensor('Cd_a')
BasicTensor.Cd.op_flg = True
BasicTensor.op_symbol = {'Pd':r'\partial','Cd':r'\nabla'}
PDF.setup(pdfmode)
return
@staticmethod
def free_indices(free_index,cntrct_index):
f = []
for i in free_index:
if i not in cntrct_index:
f.append(i)
return(f)
@staticmethod
def product_of_list(lst):
plst = lst[0]
for x in lst[1:]:
plst = plst*x
return(plst)
def __init__(self,Tstr,syms=None,sgns=None):
self.Tstr = Tstr #Basic Tensor indexed name
Tstr_lst = Tstr.split('_')
self.base = Tstr_lst[0] #Stem name
self.index = [] #Index list
self.sym = [] #Symmetry list
self.sgn = [] #Signs of symmetries
self.ccflg = [] #Covariant (True)/Contravariant (False) flag list
self.rank = [0,0] # [Contravariant rank,Contravariant rank]
self.cntrct = [] #List of contracted index pairs
self.op_flg = False #Is tensor a derivative operator
i = 1
sgn = -1
for chstr in Tstr_lst[1:]:
if chstr == '':
sgn = -sgn
else:
for ch in chstr:
self.index.append(i)
i += 1
if sgn == -1:
self.rank[0] += 1
self.ccflg.append(True)
else:
self.rank[1] += 1
self.ccflg.append(False)
sgn = -1
if syms != None:
self.add_symmetries(syms,sgns)
def rank(self):
return(self.rank)
def symmetries(self):
return(self.syms,self.sgns)
def add_symmetries(self,sym_lst,sgn_lst=None):
self.sym += sym_lst
if sgn_lst == None:
sgn_lst = len(sym_lst)*[1]
self.sgn += sgn_lst
return
def contract(self,a,b):
if self.ccflg[a-1] != self.ccflg[b-1]:
self.cntrct += [a,b]
self.rank[1] -= 1
self.rank[0] -= 1
return
else:
print 'Cannot contract indices',a,'and',b,'of',self
sys.exit(1)
def pderv(self):
return(BasicTensor.Pd*ProductTensor(self))
def cderv(self):
return(BasicTensor.Cd*ProductTensor(self))
def symbol(self):
Tstr = str(self)
Tsym = sympy.Symbol(Tstr)
return(Tsym)
class ProductTensor:
@staticmethod
def setup(pdfmode='latex',dim=2,offset=0):
BasicTensor.setup(pdfmode,dim,offset)
return
@staticmethod
def mul(T1,T2):
n = len(T1.index)
T1T2 = ProductTensor()
T1T2.plst = T1.plst+T2.plst
T1T2.index = BasicTensor.cat_shift(T1.index,T2.index)
T1T2.ccflg = T1.ccflg+T2.ccflg
T1T2.cntrct = BasicTensor.cat_shift(T1.cntrct,T2.cntrct)
T1T2.sym = copy.copy(T1.sym)
for sym in T2.sym:
T1T2.sym.append(BasicTensor.shift(n,sym))
T1T2.rank[0] = T1.rank[0]+T2.rank[0]
T1T2.rank[1] = T1.rank[1]+T2.rank[1]
return(T1T2)
@staticmethod
def derv(T,op):
dT = Tensor()
Tshft = copy.copy(T)
Tshft.index = BasicTensor.shift(1,Tshft.index)
Tshft.cntrct = BasicTensor.shift(1,Tshft.cntrct)
for i in range(len(T.sym)):
Tshft.sym[i] = BasicTensor.shift(1,Tshft.sym[i])
i = 0
for i_p in range(len(T.plst)):
s = ProductTensor()
s.plst = Tshft.plst[:i_p]+[op]+Tshft.plst[i_p:]
s.index = Tshft.index[:i]+[1]+Tshft.index[i:]
s.ccflg = Tshft.ccflg[:i]+[True]+Tshft.ccflg[i:]
s.cntrct = copy.copy(Tshft.cntrct)
s.sym = copy.copy(Tshft.cntrct)
s.rank[0] = Tshft.rank[0]+1
s.rank[1] = Tshft.rank[1]
dT = dT+Tensor(s)
i += len(T.plst[i_p].index)
dT.rank = dT.slst[0].rank
return(dT)
def __init__(self,Tstr=None):
if Tstr == None:
self.plst = []
self.cntrct = []
self.sym = []
self.ccflg = []
self.index = []
self.rank = [0,0]
else:
if isinstance(Tstr,BasicTensor):
self.plst = [Tstr]
else:
self.plst = [BasicTensor(Tstr)]
self.cntrct = self.plst[0].cntrct
self.sym = [self.plst[0].sym]
self.ccflg = self.plst[0].ccflg
self.index = self.plst[0].index
self.rank = self.plst[0].rank
def __mul__(self,T):
return(ProductTensor.mul(self,T))
def contract(self,a,b):
T = copy.copy(self)
if (T.ccflg[a-1] != T.ccflg[b-1]) and (a not in T.cntrct) and (b not in T.cntrct):
T.cntrct += [a,b]
T.rank[0] -= 1
T.rank[1] -= 1
return(T)
else:
print self,'contract('+str(a)+','+str(b)+') not allowed\n'
sys.exit(1)
def pderv(self):
return(ProductTensor.derv(self,BasicTensor.Pd))
def cderv(self):
return(ProductTensor.derv(self,BasicTensor.Cd))
def permute(self,perm_lst):
T = copy.copy(self)
tmp = copy.copy(T.index)
for i in range(len(perm_lst)):
j = perm_lst[i]-1
tmp[j] = T.index[i]
T.index = tmp
return(T)
def __str__(self):
Tstr = ''
i = 0
for p in self.plst:
j = i+len(p.index)
Tstr += p.base
old_ccflg = not self.ccflg[i]
while i < j:
new_ccflg = self.ccflg[i]
if old_ccflg != new_ccflg:
if new_ccflg:
Tstr += '_'
else:
Tstr += '__'
old_ccflg = new_ccflg
i1 = i+1
if i1 in self.cntrct:
cntrct_pos = self.cntrct.index(i1)
dummy_pos = -((cntrct_pos/2)+1)
Tstr += BasicTensor.index_pool[dummy_pos]
else:
Tstr += BasicTensor.index_pool[self.index[i]-1]
i += 1
Tstr += '*'
return(Tstr[:-1])
def latex(self):
Tstr = ''
i = 0
frst_flg = True
for p in self.plst:
j = i+len(p.index)
if p.op_flg:
Tstr += BasicTensor.op_symbol[p.base]
else:
Tstr += p.base
Tstr += r'\indices{'
old_ccflg = not self.ccflg[i]
while i < j:
new_ccflg = self.ccflg[i]
if old_ccflg != new_ccflg:
if new_ccflg:
Tstr += '_{'
else:
Tstr += '^{'
old_ccflg = new_ccflg
i1 = i+1
if i1 in self.cntrct:
cntrct_pos = self.cntrct.index(i1)
dummy_pos = -((cntrct_pos/2)+1)
Tstr += BasicTensor.index_pool[dummy_pos]
else:
Tstr += BasicTensor.index_pool[self.index[i]-1]
i += 1
Tstr += '}}'
return(Tstr)
def debug(self,i=-1):
if i < 0:
print 'Debug:'
else:
print 'Term:',i
print 'index =',self.index
print 'ccflg =',self.ccflg
print 'rank =',self.rank
print 'cntrct =',self.cntrct
print 'sym =',self.sym
return
class Tensor:
@staticmethod
def setup(pdfmode='latex',dim=2,offset=0):
ProductTensor.setup(pdfmode,dim,offset)
return
def __init__(self,Tstr=None):
self.coef = []
self.slst = []
self.rank = [0,0]
if Tstr != None:
if isinstance(Tstr,ProductTensor):
self.slst = [Tstr]
else:
self.slst = [ProductTensor(Tstr)]
self.coef = [ONE]
self.rank = self.slst[0].rank
@staticmethod
def add(T1,T2):
if T1.rank != T2.rank and T1.rank != [0,0] and T2.rank != [0,0]:
print 'In add '+str(T1)+' and '+str(T2)+' not equal in rank'
sys.exit(1)
T1pT2 = Tensor()
T1pT2.slst = T1.slst+T2.slst
T1pT2.coef = T1.coef+T2.coef
T1pT2.rank[0] = max(T1.rank[0],T2.rank[0])
T1pT2.rank[1] = max(T1.rank[1],T2.rank[1])
return(T1pT2)
@staticmethod
def sub(T1,T2):
if T1.rank != T2.rank and T1.rank != [0,0] and T2.rank != [0,0]:
print 'In sub '+str(T1)+' and '+str(T2)+' not equal in rank'
sys.exit(1)
T1sT2 = Tensor()
T1sT2.slst = T1.slst+T2.slst
T1sT2.coef = T1.coef+minus(T2.coef)
T1sT2.rank[0] = max(T1.rank[0],T2.rank[0])
T1sT2.rank[1] = max(T1.rank[1],T2.rank[1])
return(T1sT2)
@staticmethod
def mul(T1,T2):
T1T2 = Tensor()
for c1,s1 in zip(T1.coef,T1.slst):
for c2,s2 in zip(T2.coef,T2.slst):
T1T2.coef.append(c1*c2)
T1T2.slst.append(s1*s2)
T1T2.rank = T1T2.slst[0].rank
return(T1T2)
@staticmethod
def derv(T,op):
dT = Tensor()
for i in range(len(T.slst)):
dT = dT+ProductTensor.derv(T.slst[i],op)
dT.rank = dT.slst[0].rank
return(dT)
def pderv(self):
return(Tensor.derv(self,BasicTensor.Pd))
def cderv(self):
return(Tensor.derv(self,BasicTensor.Cd))
def __add__(self,T):
return(Tensor.add(self,T))
def __sub__(self,T):
return(Tensor.sub(self,T))
def __mul__(self,T):
if isinstance(T,Tensor):
return(Tensor.mul(self,T))
else:
selfT = Tensor()
selfT.slst = copy.copy(self.slst)
selfT.coef = BasicTensor.scale(T,self.coef)
selfT.rank = copy.cppy(self.rank)
return(selfT)
def __rmul__(self,T):
if isinstance(T,Tensor):
return(Tensor.mul(T,self))
else:
selfT = Tensor()
selfT.slst = copy.copy(self.slst)
selfT.coef = BasicTensor.scale(T,self.coef)
selfT.rank = copy.copy(self.rank)
return(selfT)
def __neg__(self):
T = copy.copy(self)
for i in range(len(T.coef)):
T.coef[i] = -self.coef[i]
return(T)
def contract(self,a,b):
T = copy.copy(self)
for i in range(len(T.slst)):
T.slst[i] = self.slst[i].contract(a,b)
T.rank[0] -= 1
T.rank[1] -= 1
return(T)
def permute(self,perm_lst):
permT = Tensor()
permT.coef = copy.copy(self.coef)
for pt in self.slst:
permT.slst.append(pt.permute(perm_lst))
permT.rank = copy.copy(self.rank)
return(permT)
def __str__(self):
Tsum = ZERO
for i in range(len(self.coef)):
dummy = Symbol('dummy'+str(i))
Tsum += self.coef[i]*dummy
Tstr = str(Tsum)
for i in range(len(self.coef)):
tstr = str(self.slst[i])
Tstr = Tstr.replace('dummy'+str(i),tstr)
return(Tstr)
def latex(self):
Tsum = ZERO
for i in range(len(self.coef)):
dummy = Symbol('dummy'+str(i))
Tsum += self.coef[i]*dummy
Tstr = latex(Tsum)
for i in range(len(self.coef)):
Tstr = Tstr.replace('dummy_{'+str(i)+'}',self.slst[i].latex())
return(Tstr)
def debug(self):
print 'Debug:'
for i in range(len(self.slst)):
self.slst[i].debug(i+1)
return
if __name__ == '__main__':
Tensor.setup()
g_dn = Tensor('g_ab')
g_up = Tensor('g__ab')
dg = g_dn.pderv()
f = HALF*((g_up*(dg+dg.permute([2,1]))).contract(2,5)-((g_up*dg.permute([2,3,1])).contract(2,3)))
PDF.write(f.latex())
PDF.display()
