Spherical functions

In [1]:

    

from __future__ import print_function, division, absolute_import

import numpy as np
import random
import quaternion
import numbapro as nb
from numbapro import *
import spherical_functions as sp

ru = lambda : random.uniform(-1,1)


    

In [2]:

    

j1,j2,j3,m1,m2,m3 = 1, 0, 1, 0, 0, (-0 - 0)
sp.Wigner3j(j1,j2,j3,m1,m2,m3), sp.Wigner3j(j1,(j2+j3-m1)/2,(j2+j3+m1)/2,j3-j2,(j2-j3-m1)/2+m1+m2,(j2-j3+m1)/2-m1-m2)


    

    
Out[2]:





(-0.5773502691896257, 0.0)

In [3]:

    

j1,(j2+j3-m1)/2,(j2+j3+m1)/2,j3-j2,(j2-j3-m1)/2+m1+m2,(j2-j3+m1)/2-m1-m2


    

    
Out[3]:





(1, 0.5, 0.5, 1, -0.5, -0.5)

In [5]:

    

j2,j3,m1,(j2-j3-m1)/2


    

    
Out[5]:





(0, 1, 0, -0.5)

In [2]:

    

%%timeit
for j1 in range(8):
    for j2 in range(8):
        for j3 in range(8):
            for m1 in range(-j1,j1+1):
                for m2 in range(-j2,j2+1):
                    sp.wigner3j(j1,j2,j3,m1,m2,-m1-m2)


    

    




10 loops, best of 3: 22.5 ms per loop

In [3]:

    

i=0
for j1 in range(8):
    for j2 in range(8):
        for j3 in range(8):
            for m1 in range(-j1,j1+1):
                for m2 in range(-j2,j2+1):
                    i+=1
i


    

    
Out[3]:





32768

In [4]:

    

22.5e-3/_


    

    
Out[4]:





6.866455078125e-07

In [2]:

    

q = np.quaternion(ru(), ru(), ru(), ru())
q = q/q.abs()
q


    

    
Out[2]:





quaternion(0.542904926009884, -0.513641541910117, 0.265393949689424, -0.609091667326559)

In [9]:

    

# Using slow wrapper

ells = range(16+1)
evals = np.empty_like(ells, dtype=int)
nanoseconds = np.empty_like(ells, dtype=float)
for i,ell_max in enumerate(ells):
    indices = np.array([[ell, mp, m] for ell in range(ell_max+1) for mp in range(-ell,ell+1) for m in range(-ell,ell+1)], dtype=int)
    elements = np.empty((indices.shape[0],), dtype=complex)
    result = %timeit -o sp.WignerD(q, indices)
    evals[i] = len(indices)
    nanoseconds[i] = 1e9*result.best/len(indices)
    print("With ell_max={0}, and {1} evaluations, each D component averages {2:.0f} ns".format(ell_max,evals[i],nanoseconds[i]))
    sys.stdout.flush()


    

    




100000 loops, best of 3: 5.55 µs per loop
With ell_max=0, and 1 evaluations, each D component averages 5555 ns
10000 loops, best of 3: 21.9 µs per loop
With ell_max=1, and 10 evaluations, each D component averages 2194 ns
10000 loops, best of 3: 67.2 µs per loop
With ell_max=2, and 35 evaluations, each D component averages 1920 ns
10000 loops, best of 3: 154 µs per loop
With ell_max=3, and 84 evaluations, each D component averages 1833 ns
1000 loops, best of 3: 295 µs per loop
With ell_max=4, and 165 evaluations, each D component averages 1786 ns
1000 loops, best of 3: 522 µs per loop
With ell_max=5, and 286 evaluations, each D component averages 1825 ns
1000 loops, best of 3: 825 µs per loop
With ell_max=6, and 455 evaluations, each D component averages 1813 ns
1000 loops, best of 3: 1.28 ms per loop
With ell_max=7, and 680 evaluations, each D component averages 1877 ns
1000 loops, best of 3: 1.77 ms per loop
With ell_max=8, and 969 evaluations, each D component averages 1830 ns
100 loops, best of 3: 2.32 ms per loop
With ell_max=9, and 1330 evaluations, each D component averages 1746 ns
100 loops, best of 3: 3.07 ms per loop
With ell_max=10, and 1771 evaluations, each D component averages 1731 ns
100 loops, best of 3: 4.03 ms per loop
With ell_max=11, and 2300 evaluations, each D component averages 1750 ns
100 loops, best of 3: 4.99 ms per loop
With ell_max=12, and 2925 evaluations, each D component averages 1705 ns
100 loops, best of 3: 6.16 ms per loop
With ell_max=13, and 3654 evaluations, each D component averages 1686 ns
100 loops, best of 3: 8.07 ms per loop
With ell_max=14, and 4495 evaluations, each D component averages 1796 ns
100 loops, best of 3: 9.39 ms per loop
With ell_max=15, and 5456 evaluations, each D component averages 1722 ns
100 loops, best of 3: 11.6 ms per loop
With ell_max=16, and 6545 evaluations, each D component averages 1772 ns

In [8]:

    

# Directly using jitted inner function, also creating array

ells = range(16+1)
evals = np.empty_like(ells, dtype=int)
nanoseconds = np.empty_like(ells, dtype=float)
for i,ell_max in enumerate(ells):
    indices = np.array([[ell, mp, m] for ell in range(ell_max+1) for mp in range(-ell,ell+1) for m in range(-ell,ell+1)], dtype=int)
    result = %timeit -o sp._WignerD(q.a, q.b, indices, np.empty((indices.shape[0],), dtype=complex))
    evals[i] = len(indices)
    nanoseconds[i] = 1e9*result.best/len(indices)
    print("With ell_max={0}, and {1} evaluations, each D component averages {2:.0f} ns".format(ell_max,evals[i],nanoseconds[i]))
    sys.stdout.flush()


    

    




1000000 loops, best of 3: 1.28 µs per loop
With ell_max=0, and 1 evaluations, each D component averages 1284 ns
100000 loops, best of 3: 3.06 µs per loop
With ell_max=1, and 10 evaluations, each D component averages 306 ns
100000 loops, best of 3: 8.69 µs per loop
With ell_max=2, and 35 evaluations, each D component averages 248 ns
10000 loops, best of 3: 20.3 µs per loop
With ell_max=3, and 84 evaluations, each D component averages 242 ns
10000 loops, best of 3: 41.6 µs per loop
With ell_max=4, and 165 evaluations, each D component averages 252 ns
10000 loops, best of 3: 74.3 µs per loop
With ell_max=5, and 286 evaluations, each D component averages 260 ns
10000 loops, best of 3: 122 µs per loop
With ell_max=6, and 455 evaluations, each D component averages 268 ns
10000 loops, best of 3: 186 µs per loop
With ell_max=7, and 680 evaluations, each D component averages 274 ns
1000 loops, best of 3: 271 µs per loop
With ell_max=8, and 969 evaluations, each D component averages 280 ns
1000 loops, best of 3: 381 µs per loop
With ell_max=9, and 1330 evaluations, each D component averages 286 ns
1000 loops, best of 3: 515 µs per loop
With ell_max=10, and 1771 evaluations, each D component averages 291 ns
1000 loops, best of 3: 680 µs per loop
With ell_max=11, and 2300 evaluations, each D component averages 295 ns
1000 loops, best of 3: 877 µs per loop
With ell_max=12, and 2925 evaluations, each D component averages 300 ns
1000 loops, best of 3: 1.12 ms per loop
With ell_max=13, and 3654 evaluations, each D component averages 307 ns
1000 loops, best of 3: 1.37 ms per loop
With ell_max=14, and 4495 evaluations, each D component averages 305 ns
1000 loops, best of 3: 1.65 ms per loop
With ell_max=15, and 5456 evaluations, each D component averages 303 ns
100 loops, best of 3: 1.99 ms per loop
With ell_max=16, and 6545 evaluations, each D component averages 304 ns

In [3]:

    

# After re-configuring the calculation

ells = range(16+1)
evals = np.empty_like(ells, dtype=int)
nanoseconds = np.empty_like(ells, dtype=float)
for i,ell_max in enumerate(ells):
    indices = np.array([[ell, mp, m] for ell in range(ell_max+1) for mp in range(-ell,ell+1) for m in range(-ell,ell+1)], dtype=int)
    elements = np.empty((indices.shape[0],), dtype=complex)
    result = %timeit -o sp._WignerD(q.a, q.b, indices, elements)
    evals[i] = len(indices)
    nanoseconds[i] = 1e9*result.best/len(indices)
    print("With ell_max={0}, and {1} evaluations, each D component averages {2:.0f} ns".format(ell_max,evals[i],nanoseconds[i]))
    sys.stdout.flush()


    

    




1000000 loops, best of 3: 558 ns per loop
With ell_max=0, and 1 evaluations, each D component averages 558 ns
100000 loops, best of 3: 2.04 µs per loop
With ell_max=1, and 10 evaluations, each D component averages 204 ns
100000 loops, best of 3: 7.18 µs per loop
With ell_max=2, and 35 evaluations, each D component averages 205 ns
100000 loops, best of 3: 17.7 µs per loop
With ell_max=3, and 84 evaluations, each D component averages 211 ns
10000 loops, best of 3: 35.2 µs per loop
With ell_max=4, and 165 evaluations, each D component averages 213 ns
10000 loops, best of 3: 66.7 µs per loop
With ell_max=5, and 286 evaluations, each D component averages 233 ns
10000 loops, best of 3: 109 µs per loop
With ell_max=6, and 455 evaluations, each D component averages 239 ns
10000 loops, best of 3: 165 µs per loop
With ell_max=7, and 680 evaluations, each D component averages 243 ns
1000 loops, best of 3: 243 µs per loop
With ell_max=8, and 969 evaluations, each D component averages 250 ns
1000 loops, best of 3: 336 µs per loop
With ell_max=9, and 1330 evaluations, each D component averages 252 ns
1000 loops, best of 3: 456 µs per loop
With ell_max=10, and 1771 evaluations, each D component averages 258 ns
1000 loops, best of 3: 620 µs per loop
With ell_max=11, and 2300 evaluations, each D component averages 270 ns
1000 loops, best of 3: 792 µs per loop
With ell_max=12, and 2925 evaluations, each D component averages 271 ns
1000 loops, best of 3: 1.01 ms per loop
With ell_max=13, and 3654 evaluations, each D component averages 277 ns
1000 loops, best of 3: 1.25 ms per loop
With ell_max=14, and 4495 evaluations, each D component averages 279 ns
1000 loops, best of 3: 1.57 ms per loop
With ell_max=15, and 5456 evaluations, each D component averages 287 ns
1000 loops, best of 3: 1.87 ms per loop
With ell_max=16, and 6545 evaluations, each D component averages 285 ns

In [3]:

    

# Directly using jitted inner function

ells = range(16+1)
evals = np.empty_like(ells, dtype=int)
nanoseconds = np.empty_like(ells, dtype=float)
for i,ell_max in enumerate(ells):
    indices = np.array([[ell, mp, m] for ell in range(ell_max+1) for mp in range(-ell,ell+1) for m in range(-ell,ell+1)], dtype=int)
    elements = np.empty((indices.shape[0],), dtype=complex)
    result = %timeit -o sp._WignerD(q.a, q.b, indices, elements)
    evals[i] = len(indices)
    nanoseconds[i] = 1e9*result.best/len(indices)
    print("With ell_max={0}, and {1} evaluations, each D component averages {2:.0f} ns".format(ell_max,evals[i],nanoseconds[i]))
    sys.stdout.flush()


    

    




1000000 loops, best of 3: 615 ns per loop
With ell_max=0, and 1 evaluations, each D component averages 615 ns
100000 loops, best of 3: 2.26 µs per loop
With ell_max=1, and 10 evaluations, each D component averages 226 ns
100000 loops, best of 3: 7.87 µs per loop
With ell_max=2, and 35 evaluations, each D component averages 225 ns
100000 loops, best of 3: 19.5 µs per loop
With ell_max=3, and 84 evaluations, each D component averages 233 ns
10000 loops, best of 3: 41.5 µs per loop
With ell_max=4, and 165 evaluations, each D component averages 251 ns
10000 loops, best of 3: 72.8 µs per loop
With ell_max=5, and 286 evaluations, each D component averages 255 ns
10000 loops, best of 3: 120 µs per loop
With ell_max=6, and 455 evaluations, each D component averages 264 ns
10000 loops, best of 3: 186 µs per loop
With ell_max=7, and 680 evaluations, each D component averages 273 ns
1000 loops, best of 3: 269 µs per loop
With ell_max=8, and 969 evaluations, each D component averages 278 ns
1000 loops, best of 3: 381 µs per loop
With ell_max=9, and 1330 evaluations, each D component averages 287 ns
1000 loops, best of 3: 520 µs per loop
With ell_max=10, and 1771 evaluations, each D component averages 294 ns
1000 loops, best of 3: 686 µs per loop
With ell_max=11, and 2300 evaluations, each D component averages 298 ns
1000 loops, best of 3: 876 µs per loop
With ell_max=12, and 2925 evaluations, each D component averages 300 ns
1000 loops, best of 3: 1.12 ms per loop
With ell_max=13, and 3654 evaluations, each D component averages 307 ns
1000 loops, best of 3: 1.37 ms per loop
With ell_max=14, and 4495 evaluations, each D component averages 304 ns
1000 loops, best of 3: 1.65 ms per loop
With ell_max=15, and 5456 evaluations, each D component averages 302 ns
100 loops, best of 3: 1.98 ms per loop
With ell_max=16, and 6545 evaluations, each D component averages 303 ns

In [20]:

    

# Before wrapping key math parts in jit

ells = range(1,16)
evals = np.empty_like(ells, dtype=int)
microseconds = np.empty_like(ells, dtype=float)
for i,ell_max in enumerate(ells):
    indices = np.array([[ell, mp, m] for ell in range(ell_max+1) for mp in range(-ell,ell+1) for m in range(-ell,ell+1)])
    result = %timeit -o sp.wignerD(q.a, q.b, indices)
    evals[i] = len(indices)
    microseconds[i] = 1e6*result.best/len(indices)
    print("With ell_max={0}, and {1} evaluations, each D component averages {2:.3} µs".format(ell_max,evals[i],microseconds[i]))
    sys.stdout.flush()


    

    




10000 loops, best of 3: 48.5 µs per loop
With ell_max=1, and 10 evaluations, each D component averages 4.85 µs
10000 loops, best of 3: 151 µs per loop
With ell_max=2, and 35 evaluations, each D component averages 4.32 µs
1000 loops, best of 3: 351 µs per loop
With ell_max=3, and 84 evaluations, each D component averages 4.18 µs
1000 loops, best of 3: 675 µs per loop
With ell_max=4, and 165 evaluations, each D component averages 4.09 µs
1000 loops, best of 3: 1.18 ms per loop
With ell_max=5, and 286 evaluations, each D component averages 4.13 µs
1000 loops, best of 3: 1.85 ms per loop
With ell_max=6, and 455 evaluations, each D component averages 4.07 µs
100 loops, best of 3: 2.6 ms per loop
With ell_max=7, and 680 evaluations, each D component averages 3.82 µs
100 loops, best of 3: 3.73 ms per loop
With ell_max=8, and 969 evaluations, each D component averages 3.85 µs
100 loops, best of 3: 5.12 ms per loop
With ell_max=9, and 1330 evaluations, each D component averages 3.85 µs
100 loops, best of 3: 6.78 ms per loop
With ell_max=10, and 1771 evaluations, each D component averages 3.83 µs
100 loops, best of 3: 8.8 ms per loop
With ell_max=11, and 2300 evaluations, each D component averages 3.83 µs
100 loops, best of 3: 11.1 ms per loop
With ell_max=12, and 2925 evaluations, each D component averages 3.79 µs
100 loops, best of 3: 14.4 ms per loop
With ell_max=13, and 3654 evaluations, each D component averages 3.95 µs
100 loops, best of 3: 17.4 ms per loop
With ell_max=14, and 4495 evaluations, each D component averages 3.86 µs
10 loops, best of 3: 20.3 ms per loop
With ell_max=15, and 5456 evaluations, each D component averages 3.73 µs

In [22]:

    

plot(evals, microseconds)


    

    
Out[22]:





[<matplotlib.lines.Line2D at 0x102fb5690>]






    

In [77]:

    

@njit #('void(complex128, complex128, int[:,:], complex128[:])')
def _wignerD(Ra, Rb, indices, elements):
    """Main work function for computing Wigner D matrix elements

    This is the core function that does all the work in the
    computation, but it is strict about its input, and does not check
    them for validity.

    Input arguments
    ===============
    _wignerD(Ra, Rb, indices, elements)

      * Ra, Rb are the complex components of the rotor
      * indices is an array of integers [ell,mp,m]
      * elements is an array of complex with length equal to the first
        dimension of indices

    The `elements` variable is needed because numba cannot create
    arrays at the moment, but this is modified in place and returned.

    """
    N = indices.shape[0]

    # These constants are the recurring quantities in the computation
    # of the matrix elements, so we calculate them here just once
    absRa = abs(Ra)
    absRb = abs(Rb)
    absRRatioSquared = (absRb*absRb/(absRa*absRa) if absRa>epsilon else 0.0)
    absRa_exp = (int(np.log(absRa)/_log2) if absRa>epsilon else min_exp)
    absRb_exp = (int(np.log(absRb)/_log2) if absRb>epsilon else min_exp)

    for i in range(N):
        ell = indices[i,0]
        mp = indices[i,1]
        m = indices[i,2]

        if(abs(mp)>ell or abs(m)>ell):
            elements[i] = 0.0+0.0j

        elif(absRa<=epsilon or 2*absRa_exp*(mp-m)<min_exp+mant_dig):
            if mp!=m:
                elements[i] = 0.0j
            else:
                if (ell+mp)%2==0:
                    elements[i] = Rb**(2*m)
                else:
                    elements[i] = -Rb**(2*m)

        elif(absRb<=epsilon or 2*absRb_exp*(mp-m)<min_exp+mant_dig):
            if mp!=m:
                elements[i] = 0.0j
            else:
                elements[i] = Ra**(2*m)

        else:
            rhoMin = max(0,mp-m)
            rhoMax = min(ell+mp,ell-m)
            if(absRa < 1.e-3):
                # In this branch, we deal with NANs in certain cases when
                # absRa is small by separating it out from the main sum,
                Prefactor = coeff(ell, mp, m) * Ra**(m+mp) * Rb**(m-mp)
                absRaSquared = absRa*absRa
                absRbSquared = absRb*absRb
                Sum = 0.0
                for rho in range(rhoMax, rhoMin-1, -1):
                    aTerm = absRaSquared**(ell-m-rho);
                    if(aTerm != aTerm or aTerm<1.e-100):
                        Sum *= absRbSquared
                    else:
                        if rho%2==0:
                            Sum = ( binomial_coefficient(ell+mp,rho) * binomial_coefficient(ell-mp, ell-rho-m) * aTerm
                                    + Sum * absRbSquared )
                        else:
                            Sum = ( -binomial_coefficient(ell+mp,rho) * binomial_coefficient(ell-mp, ell-rho-m) * aTerm
                                    + Sum * absRbSquared )
                elements[i] = Prefactor * Sum * absRbSquared**rhoMin
            else:
                Prefactor = coeff(ell, mp, m) * absRa**(2*ell-2*m) * Ra**(m+mp) * Rb**(m-mp)
                Sum = 0.0
                for rho in range(rhoMax, rhoMin-1, -1):
                    if rho%2==0:
                        Sum = (  binomial_coefficient(ell+mp,rho) * binomial_coefficient(ell-mp, ell-rho-m)
                                 + Sum * absRRatioSquared )
                    else:
                        Sum = ( -binomial_coefficient(ell+mp,rho) * binomial_coefficient(ell-mp, ell-rho-m)
                                 + Sum * absRRatioSquared )
                elements[i] = Prefactor * Sum * absRRatioSquared**rhoMin

#     return elements


    

In [78]:

    

ell_max=8
indices = np.array([[ell, mp, m] for ell in range(ell_max+1) for mp in range(-ell,ell+1) for m in range(-ell,ell+1)], dtype=int)
elements = np.empty((indices.shape[0],), dtype=complex)
elements_orig = np.copy(elements)
_wignerD(q.a, q.b, indices, elements)
#print(elements)
print(np.allclose(elements, elements_orig))


    

    




False

In [54]:

    

indices


    

    
Out[54]:





array([[0, 0, 0]])

In [55]:

    

ell_max


    

    
Out[55]:





0

In [43]:

    

elements, indices.shape[0]


    

    
Out[43]:





(array([  3.10503618e+231 -2.31584193e+77j]), 1)

In [46]:

    

elements.shape


    

    
Out[46]:





(1,)

In [ ]:

    




    

In [12]:

    

_wignerD.inspect_types()


    

    




_wignerD (complex128, complex128, array(int32, 1d, A, nonconst), array(complex128, 1d, A, nonconst))
--------------------------------------------------------------------------------
# File: <ipython-input-11-9925f89de8f1>
# --- LINE 5 --- 

@jit('void(c16, c16, i4[:], c16[:])')

# --- LINE 6 --- 

def _wignerD(Ra, Rb, indices, elements):

    # --- LINE 7 --- 

    """Main work function for computing Wigner D matrix elements

# --- LINE 8 --- 



    # --- LINE 9 --- 

    This is the core function that does all the work in the

    # --- LINE 10 --- 

    computation, but it is strict about its input, and does not check

    # --- LINE 11 --- 

    them for validity.

# --- LINE 12 --- 



    # --- LINE 13 --- 

    Input arguments

    # --- LINE 14 --- 

    ===============

    # --- LINE 15 --- 

    _wignerD(Ra, Rb, indices, elements)

# --- LINE 16 --- 



      # --- LINE 17 --- 

      * Ra, Rb are the complex components of the rotor

      # --- LINE 18 --- 

      * indices is an array of integers [ell,mp,m]

      # --- LINE 19 --- 

      * elements is an array of complex with length equal to the first

        # --- LINE 20 --- 

        dimension of indices

# --- LINE 21 --- 



    # --- LINE 22 --- 

    The `elements` variable is needed because numba cannot create

    # --- LINE 23 --- 

    arrays at the moment, but this is modified in place and returned.

# --- LINE 24 --- 



    # --- LINE 25 --- 

    """

# --- LINE 26 --- 



    # --- LINE 27 --- 

    # These constants are the recurring quantities in the computation

    # --- LINE 28 --- 

    # of the matrix elements, so we calculate them here just once

    # --- LINE 29 --- 
    # label 0
    #   Ra.1 = Ra  :: pyobject
    #   del Ra
    #   Rb.1 = Rb  :: pyobject
    #   del Rb
    #   indices.1 = indices  :: pyobject
    #   del indices
    #   elements.1 = elements  :: pyobject
    #   del elements
    #   $0.1 = global(abs: <built-in function abs>)  :: pyobject
    #   $0.3 = call $0.1(Ra.1, )  :: pyobject
    #   del Ra.1
    #   del $0.1
    #   absRa = $0.3  :: pyobject
    #   del $0.3

    absRa = abs(Ra)

    # --- LINE 30 --- 
    #   $0.4 = global(abs: <built-in function abs>)  :: pyobject
    #   $0.6 = call $0.4(Rb.1, )  :: pyobject
    #   del Rb.1
    #   del $0.4
    #   absRb = $0.6  :: pyobject
    #   del $0.6

    absRb = abs(Rb)

    # --- LINE 31 --- 
    #   $0.8 = global(epsilon: 1e-14)  :: pyobject
    #   $0.9 = absRa > $0.8  :: pyobject
    #   del $0.8
    #   branch $0.9, 36, 54
    # label 36
    #   del $0.9
    #   $absRb36.1 = absRb  :: pyobject
    #   $36.3 = $absRb36.1 * absRb  :: pyobject
    #   del $absRb36.1
    #   $absRa36.4 = absRa  :: pyobject
    #   $36.6 = $absRa36.4 * absRa  :: pyobject
    #   del $absRa36.4
    #   $36.7 = $36.3 /? $36.6  :: pyobject
    #   del $36.6
    #   del $36.3
    #   $phi57.1 = $36.7  :: pyobject
    #   del $36.7
    #   jump 57
    # label 54
    #   del $0.9
    #   $const54.1 = const(float, 0.0)  :: pyobject
    #   $phi57.1 = $const54.1  :: pyobject
    #   del $const54.1
    #   jump 57
    # label 57
    #   absRRatioSquared = $phi57.1  :: pyobject
    #   del absRRatioSquared
    #   del $phi57.1

    absRRatioSquared = (absRb*absRb/(absRa*absRa) if absRa>epsilon else 0.0)

    # --- LINE 32 --- 
    # label 97
    #   del absRa
    #   del $57.4
    #   $97.1 = global(min_exp: -1021)  :: pyobject
    #   $phi100.1 = $97.1  :: pyobject
    #   del $97.1
    #   jump 100
    # label 100
    #   absRa_exp = $phi100.1  :: pyobject
    #   del absRa_exp
    #   del $phi100.1
    # label 72
    #   del $57.4
    #   $72.1 = global(int: <type 'int'>)  :: pyobject
    #   $72.2 = global(np: <module 'numpy' from '/Users/boyle/.continuum/anaconda/lib/python2.7/site-packages/numpy/__init__.pyc'>)  :: pyobject
    #   $72.3 = getattr(attr=log, value=$72.2)  :: pyobject
    #   del $72.2
    #   $72.5 = call $72.3(absRa, )  :: pyobject
    #   del absRa
    #   del $72.3
    #   $72.6 = global(_log2: 0.69314718056)  :: pyobject
    #   $72.7 = $72.5 /? $72.6  :: pyobject
    #   del $72.6
    #   del $72.5
    #   $72.8 = call $72.1($72.7, )  :: pyobject
    #   del $72.7
    #   del $72.1
    #   $phi100.1 = $72.8  :: pyobject
    #   del $72.8
    #   jump 100
    #   $57.3 = global(epsilon: 1e-14)  :: pyobject
    #   $57.4 = absRa > $57.3  :: pyobject
    #   del $57.3
    #   branch $57.4, 72, 97

    absRa_exp = (int(np.log(absRa)/_log2) if absRa>epsilon else min_exp)

    # --- LINE 33 --- 
    #   $100.3 = global(epsilon: 1e-14)  :: pyobject
    #   $100.4 = absRb > $100.3  :: pyobject
    #   del $100.3
    #   branch $100.4, 115, 140
    # label 140
    #   del absRb
    #   del $100.4
    #   $140.1 = global(min_exp: -1021)  :: pyobject
    #   $phi143.1 = $140.1  :: pyobject
    #   del $140.1
    #   jump 143
    # label 143
    #   absRb_exp = $phi143.1  :: pyobject
    #   del absRb_exp
    #   del $phi143.1
    # label 115
    #   del $100.4
    #   $115.1 = global(int: <type 'int'>)  :: pyobject
    #   $115.2 = global(np: <module 'numpy' from '/Users/boyle/.continuum/anaconda/lib/python2.7/site-packages/numpy/__init__.pyc'>)  :: pyobject
    #   $115.3 = getattr(attr=log, value=$115.2)  :: pyobject
    #   del $115.2
    #   $115.5 = call $115.3(absRb, )  :: pyobject
    #   del absRb
    #   del $115.3
    #   $115.6 = global(_log2: 0.69314718056)  :: pyobject
    #   $115.7 = $115.5 /? $115.6  :: pyobject
    #   del $115.6
    #   del $115.5
    #   $115.8 = call $115.1($115.7, )  :: pyobject
    #   del $115.7
    #   del $115.1
    #   $phi143.1 = $115.8  :: pyobject
    #   del $115.8
    #   jump 143

    absRb_exp = (int(np.log(absRb)/_log2) if absRb>epsilon else min_exp)

# --- LINE 34 --- 



    # --- LINE 35 --- 
    # label 149.1
    #   $const149.1.1 = const(LiftedLoop, LiftedLoop(<function _wignerD at 0x111f70aa0>))  :: pyobject
    #   $149.1.4 = call $const149.1.1(indices.1, elements.1, )  :: pyobject
    #   del indices.1
    #   del elements.1
    #   del $const149.1.1
    #   del $149.1.4
    #   jump 233
    #   jump 149.1

    for i,(ell,mp,m) in enumerate(indices):

# --- LINE 36 --- 



        # --- LINE 37 --- 

        if(abs(mp)>ell or abs(m)>ell):

            # --- LINE 38 --- 
            # label 233
            #   $const233.1 = const(NoneType, None)  :: pyobject
            #   $233.2 = cast(value=$const233.1)  :: pyobject
            #   del $const233.1
            #   return $233.2

            elements[i] = 0.0+0.0j

# The function contains lifted loops
# Loop at line 35
# Has 0 overloads

================================================================================

In [ ]: