4
$\begingroup$

I am trying to collect elements with Reap and Sow, where I am interested only in non-zero elements, that is I am trying something like that

Flatten@Last@Reap@Do[
        Sow[Cases[{preThreeJSymbols[l1,m1,l2,m2,k]},Except[0]]],
        {l1,0,3},{m1,-3,3},{l2,0,3},{m2,-3,3},{k,0,6}]

where preThreeJSymbols is given by

preThreeJSymbols[l1_, m1_, l2_, m2_, k_] := 
 If[Abs[l1 - l2] <= k <= l1 + l2 && 
    Abs[m1 - m2] <= k && Abs[m1] <= l1 && 
    Abs[m2] <= l2,
    ThreeJSymbol[{l1, 0.}, {k, 0.}, {l2, 0.}]*
    ThreeJSymbol[{l1, -m1}, {k, m1 - m2}, {l2, m2}], 0]

I have two issues here:

First, is this the most efficient/natural way of doing it in Mathematica, I have the feeling that there will be a solution with pure functions (some condition/pattern matching?) Secondly, in my case preThreeJSymbols[l1,m1,l2,m2,k] seems to return both 0 and 0.. (I don't know why - I am confused here), which are interpreted differently by Except[], hence not all zeroes are discarded.

Probably, I have to add that I do not want to use Table/ParallelTable, as in my application the ranges of l1, m1, l2, m2 and k are much larger and the resulting table has more than 100 million elements, which makes it too big for my RAM to store. However, I know that most of the elements will be zero (around 5-6% are non-zero), so I want to go over all the elements and pick only the non-vanishing ones without actually generating a massive table.

Late Update

Created a new question as per a suggestion

I have implemented Mr. Wizard's solution as

precomputeThreeJTable[Lmax_] := (Clear[precomputedSymbols]; 
precomputedSymbols = 
Reap[Do[preThreeJSymbols[l1, m1, l2, m2, k] // 
   If[# != 0, 
     Sow[{pairFunction[pairFunction[l1 + 1, m1 + Lmax + 1], 
        pairFunction[l2 + 1, pairFunction[m2 + Lmax + 1, k + 1]]],
        N[#]}]] &,
   {l1, 0, Lmax}, {m1, -Lmax, Lmax}, {l2, 0, Lmax}, {m2, Lmax, Lmax},
   {k, 0, 2 Lmax}]][[2, 1]]; 
threeJTable = 
AssociationThread[precomputedSymbols[[All, 1]], precomputedSymbols[[All, 2]]]; 
Export[StringInsert["somewhere", ToString[Lmax], -6], threeJTable, "WDX"];)

where preThreeJSymbols is as above and pairFunction is given by

pairFunction[x_, y_] := ((x + y) (x + y + 1))/2 + y

And I want to do this for Lmax around 45. At the moment I can manage only 30 - the reason being that Mathematica eats all my RAM. I have noticed that when I run precomputeThreeJTable, at the end of the calculation (when the Association is being created I assume) there is a massive spike in RAM usage, but after the command is evaluated the spike is gone and I can go on using threeJTable in other calculations with far less RAM being taken. This spike in RAM usage is also seen when I Import the exported file (the file is around 200mb). For Lmax=35 my laptop can't take the memory spike. My question is whether it is possible to avoid this somehow - maybe create (and later Import) the Association in pieces or something else? I also tried using a SparseArray instead, replacing the AssociationThread code by

SparseArray[Table[precomputedSymbols[[i, 1]] -> precomputedSymbols[[i, 2]],
           {i, 1, Length[precomputedSymbols]}]]

The spkie is still there. Not only that but another problem comes up, due to pairFunction's nature - it creates very huge integers pretty easily, which means that the dimension of the SparseArray will be enormous, so when I try using it in a calculation, even for Lmax=3, my Mathematica just resets (that is - clears all definitions - like when it runs out of memory). I can show the code for the way it is used, though it is a bit long. Is this expected behaviour or I am probably doing something wrong.

Thank you

$\endgroup$
7
  • 1
    $\begingroup$ Perhaps add the code for f? $\endgroup$
    – ciao
    May 25, 2015 at 8:55
  • $\begingroup$ Also, the code doesn't have correct opening/closing {} and []. I guess this is correct: Flatten@Last@Reap@Do[Sow[Cases[{f[x, y, z, w]}, Except[0]]], {x, 0, 5}, {y, 0, 5}, {z, 0, 5}, {w, 0, 5}] $\endgroup$
    – Mahdi
    May 25, 2015 at 8:58
  • $\begingroup$ It looks to me like all zeroes disappear from your result? Can you explain further what you mean by "both 0 and 0"? $\endgroup$
    – C. E.
    May 25, 2015 at 10:35
  • 1
    $\begingroup$ Perhaps a simple With[{el = preThre..}, If[el != 0, Sow[el]] would be better. $\endgroup$
    – C. E.
    May 25, 2015 at 11:03
  • $\begingroup$ Yep, that works fine! The zeroes were due to me putting in 0. as an argument. $\endgroup$ May 25, 2015 at 12:41

1 Answer 1

5
$\begingroup$

Inexact zeros

preThreeJSymbols[l1,m1,l2,m2,k] seems to return both 0 and 0.

You introduced inexact arithmatic in the expression:

ThreeJSymbol[{l1, 0.}, {k, 0.}, {l2, 0.}]

You can either add Chop to your result or make these zeros exact:

ThreeJSymbol[{l1, 0}, {k, 0}, {l2, 0}]

However if you use expr != 0 as Pickett recommended you do not need either.

Cases

You are using Cases on an expression with only one element when it is ostensibly intended to be used on lists of elements therefore other methods are likely to be faster. Your use also results in sowing an empty list {} for every zero which actually takes more RAM:

ByteCount[0]
ByteCount[{}]
16

40

Try this instead:

Reap[Do[
  preThreeJSymbols[l1, m1, l2, m2, k] // If[# != 0, Sow[#]] &,
  {l1, 0, 3}, {m1, -3, 3}, {l2, 0, 3}, {m2, -3, 3}, {k, 0, 6}
]][[2, 1]]

Possible optimization

It seems that ThreeJSymbol is defined in top-level Mathematica code. You may read its definition in version 10.1 with:

Needs["GeneralUtilities`"]

PrintDefinitions @ ThreeJSymbol

(See What is the most convenient way to read definitions of in-memory symbols when we don't have the source files? (Spelunking tools) for other versions.)

If you can vectorize this code you may be able produce a significant performance improvement.

$\endgroup$
3
  • $\begingroup$ Your answer helped me quite a lot and I have been playing with various things since then and now I have a different issue, but I didn't know whether to create a new question or not, so I just edited my post $\endgroup$ Jun 17, 2015 at 19:13
  • $\begingroup$ @user1482714 I think it is better to post a late update like this as a new question. The update is about performance optimization of a specific piece of code rather than a general method, IMO. $\endgroup$
    – Mr.Wizard
    Jun 18, 2015 at 5:58
  • $\begingroup$ I created a new question and updated this one as well with the link $\endgroup$ Jun 18, 2015 at 8:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.