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)
--
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.
