fr., 03.09.2010 kl. 09.08 -0700, skrev klmn:
> 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)…
What do you mean by "force garbage collection"? Have you tried
>>> sympy.cache.clear_cache()
Øyvind
>
> 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.
