add wigner9j symbols

[tgorordo]

Hi!

I've taken a rough pass at minimally implementing the 9j symbols based on L. Wei (1998) with binomials optimized according to the methods in Johansson & Forssen (2015), per the README's TODO. I'm new to package-quality Julia, but tried to use mostly similar patterns to the older code - for readability and minimal impact on that existing code.

Based on a bit of profiling, the computation runs reasonably quickly, is exact, and tests well against the decomposition in terms of sums of products of 6j symbols. There's likely still some performance to be gained.

I'd be happy to take any feedback/make any edits! :) Sorry for any trouble due to inexperience packaging.

[Jutho]

Thanks for this and my apologies for the late response. I was on holidays in the second half of August. At a first glance, this looks of great quality. I have to admit that I will have to quickly read up on 9j symbols before doing a proper review; I have no experience with them (or not knowingly at least; I decompose everything in 6js, maybe without realising that what I am doing are 9j calculations).

[Jutho]

Was there a specific reason for calling the last one β₀ and not β₄?

[Jutho]

Also, the construction of the α and β is the same for b₁, b₂ and b₃. Was there a specific reason not to include this as the start of the wei9jbracket function, such that here you could simply have

b₁ = wei9jbracket(a, b, c, f, j, k)
b₂ = wei9jbracket(f, d, e, h, b, k)
b₃ = wei9jbracket(h, j, g, a, d, k)

[tgorordo]

Was there a specific reason for calling the last one β₀ and not β₄?

No real reason - loosely the last one stood out to me as a break in a pattern since it doesn't have a matching αᵢ (wrt. the mᵢ used to form it), so I picked the 0 subscript to distinguish it a bit at a glance. I'd be happy to relabel it to β₄ if you prefer!

[tgorordo]

Also, the construction of the α and β is the same for b₁, b₂ and b₃...

wrt. the b₁, b₂, b₃ I was trying to maintain readability of the various bracketed (a,b,c,d,e,f,g,h,j, k) combinations against the Wei paper - but was also generally following the pattern from the 3j and 6j code where conversions tended to be performed before passing into a more specialized computational function. In this case, the mᵢ are generally BigHalfInts but the relevant linear combinations always come out convertible to BigInts, which are what the primebinomial calls in wei9jbracket were going to need. This could be refactored pretty easily if you prefer, since that's usually just a good type-stability pattern but things should be type stable either way here.

[Jutho]

Not sure if this needs to be exported. I assume it is an internal function and users would only call wigner9j?

[tgorordo]

Oops - that's a leftover from some spot-testing while I was writing my first pass. I don't think it needs exporting either. Will remove.

[Jutho]

if you simply do

divexact!(num, primefactorial(k))
divexact!(num, primefactorial(n-k))

there is no need to take a copy of primefactorial(k). divexact! should be as quick as mul!, as it is just subtracting rather than adding prime powers, so there is no need to avoid doing it twice by first doing a mul!

[tgorordo]

Thanks! Duh, will improve that!

[Jutho]

I would prefer an explicit return in front of this sum, or first assigning it to a variable and then returning that. It took me a few seconds to realise this was the return value.

[tgorordo]

Can do! I was split on whether the implicit return value was bad style myself.

[Jutho]

So if I understand correctly, the result of compute9jseries here is a BigInt; there are no denominators involved. Maybe then it could be slightly faster still to first take out primefactors from sden and only then construct the rational, i.e. something like

snum2 = compute9jseries(j₁, j₂, j₃, j₄, j₅, j₆, j₇, j₈, j₉)
for n = 1:length(sden.powers)
    p = bigprime(n)
    while sden.powers[n] > 0
        q, r = divrem(snum2, p)
        if iszero(r)
            snum2 = q
            sden.powers[n] -= 1
        else
            break
        end
    end
end
Base.unsafe_rational(convert(BigInt, snum)*snum2, convert(BigInt, sden))

[tgorordo]

Yeah - the series factor in the 9j formula given by Wei is always just an integer so compute9jseries returns a BigInt (actually I think it returns a BigHalfInt, but it should actually always be integer so I'll add that conversion to the end of compute9jseries for posterity). So that's a good catch - I'll add that optimization.

[Jutho]

I have to admit that I am confused by how the permutations are encoded here and so how this code functions exactly. Would it be possible to follow the 3j and 6j philosophy and try to impose a canonical order. For example, given all the symmetries of the 9j symbol, it seems like one could always impose:

j₁ <= j₂ <= j₄ <= min(j₃, j₅, j₆, j₇, j₈, j₉)

When the inequalities are strict, I believe this fixes the order completely, i.e. any permuation or reflection would destroy this order. However, when there are equalities, further conditions could be imposed.

[tgorordo]

I was following the approach mentioned by the wigxjpf authors in this slide-deck (3rd slide from the end), where the symmetry checks for the 9j are done explicitly for the 72 equivalent permutations+transposition (up to sign) against the cache, rather than using a canonical ordering like the 3j or 6j.

The rationale there being (I believe) that unlike what happens in the 3j or 6j cases, the semimagic square symmetries give an under-constraining 6 conditions on 9 linearly independent entries rather than over-constraining the system to give a single canonical ordering.

My implementation could probably be cleaner or made more readable, but I'm using that Julia arrays can be indexed with vectors to get a permutation, e.g.

A = ['a', 'b', 'c']

A[[2, 3, 1]] # = ['b', 'c', 'a']

B = [1 2 3; 4 5 6; 7 8 9]

# swap the first two rows, and last two columns
B[[2, 1, 3], [1, 3, 2]] # = [4 6 5; 1 3 2; 7 9 8]

and then hard-coded all the column and row permutations accordingly (in lexicographic order) as well as the parity of each permutation. The transpositions get checked too but aren't encoded in the _perms9j table of indices, hence the kk and kkT cases.

I think you're right though, and like your suggested condition as a starting point; it might be plausible to make a choice of ordering that amounts to a choice of the last few constraints, and handle a few cases to do better than checking all permutations. I'm out of time to spare today, so I'll sit down with this one later in the week to work out the cases (for when there are equalities) and how to cast an arbitrary symbol into the resulting form to see exactly how this can be improved!

[Jutho]

Great. Looking at the slides, I indeed remember that for 3j and 6j you can go further and really compute a unique linear index for each of the symbols (such that those that are related up to a sign by symmetries have indeed the same linear index). That's probably too much to ask for the 9j symbols, but I don't use that anyway. I simply use the symmetries to reduce them to canonical order, but then still rely on a hash table.

[tgorordo]

I've been busier than planned, so I've only been able to take a brief of a stab at a canonical form. The condition you suggest as a starting point unfortunately doesn't work: if a symbol can be canonicalized to meet it, then it's a good condition in that any row or col swap or transposition would break it - however, not all symbols can be canonicalized to satisfy that condition in the first place.

The issue is that anything that sets a particular entry over-constrains the permutations that can follow. e.g. in the suggested canonicalization j₁ must be made the smallest number in the symbol but doing so disallows the 1st row and 1st column from being involved in future row or column swaps to set the other entries, but in order to move the 2nd smallest number into j₂ it must have landed in the 1st row or the 1st column during the permutation that set j₁ - which there's no guarantee it will have. Because of patterns like this I'm finding most canonicalization schemes I can come up with in the same vein require specifying many cases/conditions and get a bit ugly (and aren't looking much better than just checking the permutations).

I'll see if I can find some better way to approach canonicalization over the week here - apologies for the delay.

[Jutho]

Maybe replace this one with

let ...

to make Julia 1.6 happy. However, I am happy to bump the minimum Julia version.

[tgorordo]

Oops - can do! I should have paid attention to that compat check; my bad - that made failing test messages a bit more readable, but I see no need to regress for that.

[tgorordo]

Awesome, thanks for the feedback! I've gone through and tweaked functionality/performance things as recommended, and replied to comments wrt. more stylistic things (if you confirm preferences for them I'm also happy to make those changes definitive too!). I'll need to loop back later in the week to take a look at improving the caching scheme/checks though - I'm out of bulk time to spare today :)

[Jutho]

Ok I see, none of the symmetry transformations can actually change what is living on a given row or column, except for interchanging rows and columns. I guess you can still have that:

j1 <= everything (bring smallest element to the upper left corner)
row 1 is sorted such that j1 <= j2 <= j3 (by permuting column 2 and 3 if necessary)
column 1 is sorted such that j1 <= j4 <= j7 (by permuting row 2 and 3 if necessary)
j2 <= j4 (by transposing if necessary)

I think that fixes all the symmetry transformations when all the j's are different, but still leaves the degenerate cases which require further attention.

I might try to take a stab at this.