klmn
unread,Sep 3, 2010, 12:08:14 PM9/3/10Sign in to reply to author
Sign in to forward
You do not have permission to delete messages in this group
Either email addresses are anonymous for this group or you need the view member email addresses permission to view the original message
to sympy
Dear all,
I am trying to implement various schemes for storing Gaunt
coefficients (integrals of 3 spherical harmonics) calculated
symbolically. Ondrej, I remember you offered me help for this, now I
have to ask for it.
It appears that at moderately big indices, already at l_max~30, the
memory is all eaten up. I tried to force garbage collection, without
success. Interestingly, if I repeatedly calculate just one Gaunt
coefficient for a fixed set of indices, there is no memory problem
(commented out in the script below)…
There is no memory leak per se. When I exit from Python, the memory is
fully recovered. It seems that the garbage is not collected during the
calculations, even when asked, untill it occupies the whole memory.
Below I attach a test module. Sorry for the long post: I don’t see a
way how I could attach it as a file (I am using a web form).
Can you recommend me a memory saving approach?
Best regards,
Konstantin
****************************************
import numpy as np
import sympy
import tables
#import gc
import time
factorials=[sympy.Integer(1)]
def get_factorials(nn):
""" Precalculated factorials """
if nn >= len(factorials):
for ii in range(len(factorials), nn + 1):
factorials.append(factorials[ii - 1] * ii)
return factorials[:(nn + 1)]
def zero_gaunt_l(l_1, l_2, l_3):
"""
Checks l indices for giving an even sum and fulfilling the
triangle realtion.
"""
if np.mod(l_1 + l_2 + l_3, 2) != 0: return True
if l_1 + l_2 - l_3 < 0: return True
if l_1 - l_2 + l_3 < 0: return True
if -l_1 + l_2 + l_3 < 0: return True
return False
def zero_gaunt_m(l_1, l_2, l_3, m_1, m_2, m_3):
"""
Checks m indices for giving 0 sum and being limited by l's
"""
if (m_1 + m_2 + m_3) != 0: return True
if (abs(m_1) > l_1) or (abs(m_2) > l_2) or (abs(m_3) > l_3):
return True
return False
def gaunt(l_1, l_2, l_3, m_1, m_2, m_3, prec=None, include4Pi=True):
"""
Check validity of the indices and calculate the Gaunt coefficient.
"""
if type(l_1) != type(l_2) != type(l_3) != type(1):
raise ValueError("l values must be integer")
if type(m_1) != type(m_2) != type(m_3) != type(1):
raise ValueError("m values must be integer")
if zero_gaunt_l(l_1, l_2, l_3): return 0
if zero_gaunt_m(l_1, l_2, l_3, m_1, m_2, m_3): return 0
return calc_gaunt(l_1, l_2, l_3, m_1, m_2, m_3, prec, include4Pi)
def calc_gaunt(l_1, l_2, l_3, m_1, m_2, m_3, prec=None,
include4Pi=True):
"""
Calculate the Gaunt coefficient.
"""
a1 = l_1 + l_2 - l_3
a2 = l_1 - l_2 + l_3
a3 = -l_1 + l_2 + l_3
bigL = (l_1 + l_2 + l_3) / 2
a4 = -l_3 + l_1 + m_2
a5 = -l_3 + l_2 - m_1
imin = max(a4, a5, 0)
a6 = l_2 + m_2
a7 = l_1 - m_1
imax = min(a6, a7, a1)
maxFactorial = max(l_1 + l_2 + l_3, imax) + 1
get_factorials(maxFactorial)
trySeparate = True # test whether the factors under sqrt must be
kept
if trySeparate:
prefac1 = sympy.sqrt(2 * l_1 + 1) *\
sympy.sqrt(2 * l_2 + 1) *\
sympy.sqrt(2 * l_3 + 1) *\
sympy.sqrt(factorials[a7]) *\
sympy.sqrt(factorials[l_1 + m_1]) *\
sympy.sqrt(factorials[l_2 - m_2]) *\
sympy.sqrt(factorials[a6]) *\
sympy.sqrt(factorials[l_3 - m_3]) *\
sympy.sqrt(factorials[l_3 + m_3])
else:
prefac1 = sympy.sqrt((2 * l_1 + 1) * (2 * l_2 + 1) * (2 * l_3
+ 1) *\
factorials[a7] *\
factorials[l_1 + m_1] *\
factorials[l_2 - m_2] *\
factorials[a6] *\
factorials[l_3 - m_3] *\
factorials[l_3 + m_3])
if include4Pi:
prefac1 /= (sympy.sqrt(sympy.pi) * sympy.Integer(2))
prefac2 = factorials[bigL] * factorials[a1] * factorials[a2] *
factorials[a3] /\
(factorials[2 * bigL + 1] *\
factorials[bigL - l_1] * factorials[bigL - l_2] *
factorials[bigL - l_3])
sumres = 0
for ii in range(imin, imax + 1):
den = factorials[ii] * factorials[a1 - ii] *\
factorials[ii - a4] * factorials[ii - a5] *\
factorials[a6 - ii] * factorials[a7 - ii]
if np.mod(ii, 2) == 0:
sumres += sympy.Integer(1) / den
else:
sumres -= sympy.Integer(1) / den
res = prefac1 * prefac2 * sumres
if np.mod((bigL + l_3 + m_1 - m_2), 2) != 0:
res = -res
if prec is not None:
res = res.evalf(prec)
return res
def store_Rasch_Yu(fileName, l_max):
filters = tables.Filters(complevel=1, complib='zlib')
#complib='zlib'|'blosc'
with tables.openFile(fileName, mode = "w") as fileh:
vlarray0 = fileh.createVLArray(fileh.root, 'vlarray',
tables.ObjectAtom(),
"pickled Gaunt coeffs", filters=filters,
expectedsizeinMB=(l_max*0.1)**4)
zeroList = []
indices = 0
for l_1 in xrange(l_max + 1):
for l_2 in xrange(l_1 + 1):
for l_3 in xrange(l_2 + 1):
if zero_gaunt_l(l_1, l_2, l_3):
for m_3 in xrange(l_3 + 1):
# vlarray0.append(zeroList)
indices += 1
continue
print l_1, l_2, l_3
for m_3 in xrange(l_3 + 1):
m3list = []
for m_2 in xrange(-l_2, min(l_1 - m_3, l_2) +
1):
m_1 = -m_2 - m_3
if zero_gaunt_m(l_1, l_2, l_3, m_1, m_2,
m_3):
curGaunt = sympy.Integer(0)
else:
curGaunt = calc_gaunt(l_1, l_2, l_3,
m_1, m_2, m_3,
include4Pi=False)
# if I comment calc_gaunt above and uncomment the one below with a
dummi
# set of indices, there is no memory problem
# curGaunt =
calc_gaunt(29,29,34,10,-5,-5,
# include4Pi=False)
if type(curGaunt) in \
[sympy.core.numbers.Integer,
sympy.core.numbers.Zero,
sympy.core.numbers.One,
sympy.core.numbers.NegativeOne]:
m2list = curGaunt
elif type(curGaunt) ==
sympy.core.numbers.Rational:
m2list = [curGaunt.p, curGaunt.q]
elif type(curGaunt) == sympy.core.mul.Mul:
m2list = [0, 1, 0]
for bas_arg in
curGaunt.iter_basic_args():
if type(bas_arg) ==
sympy.core.numbers.Rational:
m2list[0] = bas_arg.p
m2list[1] = bas_arg.q
elif type(bas_arg) ==
sympy.core.power.Pow:
if bas_arg.exp ==
sympy.Integer(1)/2:
m2list[2] = bas_arg.base
else:
raise ValueError('wrong type
of Gaunt atoms')
if (m2list[0] == 0) or (m2list[2] ==
0):
raise ValueError('wrong value of
Gaunt')
else:
print l_1, l_2, l_3, m_1, m_2, m_3
print curGaunt, type(curGaunt)
raise ValueError('wrong type of Gaunt
result')
m3list.append(m2list)
# vlarray0.append(m3list)
indices += 1
# vlarray0.flush()
# print 'collecting garbage...',
# gc.collect()
# print ' done'
max_index_expected = (l_max + 1) * (l_max + 2) * (l_max + 3) *
(l_max + 4) / 24
print indices, max_index_expected
if __name__ == '__main__':
#!!!!!!!!!!!!!!!! VERY important !!!!!!!!!!!!!!!!
#in sympy/core/numbers.py, line 1013:
## The following is an algorithm where we collect perfect roots
## from the factors of base
#if b > 4294967296:
# # Prevent from factorizing too big integers
# return None
# <------------------- comment it out!!!
#end of !!!!!!!!!!!!!!!! VERY important !!!!!!!!!!!!!!!!
tstart = time.time()
store_Rasch_Yu('test_Rasch_Yu.h5', 30)
tstop = time.time()
print 'The calculation took {0:.1f} s'.format(tstop - tstart)