Order unity errors when reconstructing real-space function about spinsfast HOT 4 CLOSED

This is expected behavior; it is a mathematical result, rather than just something about the code. There are two things you are doing "wrong" with your test function which makes it impossible for a finite expansion in s=1 spin-weighted spherical harmonics (SWSHs) to reproduce your input.

First, you are demanding that the function have spin weight 1, but a nonzero average value. Fields of nonzero spin weight must have average value 0, since modes with ell < abs(s) are necessarily zero by the nature of spin-weighted fields.

Second, you can't expect to resolve a discontinuous function very well with just lmax=8.
What you're seeing is related to the Gibbs phenomenon. Just as a truncated Fourier expansion of a step function has wiggles, the truncated SWSH expansion of your two-dimensional step function has two-dimensional wiggles. To see this more clearly, you can increase the resolution, and you'll find that the errors are mostly at the edges of your nonzero pixel.

There's a little subtlety here: as you add pixels, if you keep defining your function with f_in[2, 4] = 1.0 + 1.0j, you're changing the function. So I'll redefine your function as being 1.0 + 1.0j anywhere in the space occupied by the [2, 4] pixel when lmax=8:

theta_phi_8 = np.array([[[theta, phi]
                         for phi in np.linspace(0.0, 2*np.pi, num=2*8+1, endpoint=False)]
                        for theta in np.linspace(0.0, np.pi, num=2*8+1, endpoint=True)])
theta_min = (theta_phi_8[1, 4, 0] + theta_phi_8[2, 4, 0])/2
theta_max = (theta_phi_8[2, 4, 0] + theta_phi_8[3, 4, 0])/2
phi_min = (theta_phi_8[2, 3, 1] + theta_phi_8[2, 4, 1])/2
phi_max = (theta_phi_8[2, 4, 1] + theta_phi_8[2, 5, 1])/2

def f(theta, phi):
    if (theta > theta_min and theta < theta_max and
        phi > phi_min and phi < phi_max):
        return 1.0 + 1.0j
    else:
        return 0.0j

Then, we can remake your images above with a higher lmax:

lmax = 256
f_in = np.array([[f(theta, phi)
                  for phi in np.linspace(0.0, 2*np.pi, num=2*lmax+1, endpoint=False)]
                 for theta in np.linspace(0.0, np.pi, num=2*lmax+1, endpoint=True)])
alm2 = spinsfast.map2salm(f_in, s, lmax)
f_out = spinsfast.salm2map(alm2, s, lmax, 2*lmax+1, 2*lmax+1)
print("Maximum deviation: ", np.max(abs(f_in - f_out))/abs(np.max(f_in)))
plt.close('all')
fig, axes = plt.subplots(figsize=[18, 4], ncols=3)
draw(axes[0], f_in.real, 'input')
draw(axes[1], f_out.real, 'output')
draw(axes[2], abs(f_in - f_out), 'difference')
plt.show()

Maximum deviation: 0.514373151752

As expected, the errors are localized around that original pixel.

A much better test is to see if the difference converges to smaller values as you increase the lmax value. Just as with Gibbs phenomena, we don't expect the maximum value of the difference to keep getting smaller with increasing resolution; we only expect the L2 norm of the difference to get smaller. As with most spectral methods, I would expect roughly exponential convergence. So, with those caveats, the following code is a better test for your function:

pixel_area = - (phi_max - phi_min) * (np.cos(theta_max) - np.cos(theta_min))
f_average = (1.0 + 1.0j) * pixel_area / (4*np.pi)

ells = np.arange(8, 256+1, 8)
L2norms_of_differences = np.zeros_like(ells, dtype=float)

for i, lmax in enumerate(ells):
    f_in = np.array([[f(theta, phi) - f_average
                      for phi in np.linspace(0.0, 2*np.pi, num=2*lmax+1, endpoint=False)]
                     for theta in np.linspace(0.0, np.pi, num=2*lmax+1, endpoint=True)])
    alm2 = spinsfast.map2salm(f_in, s, lmax)
    f_out = spinsfast.salm2map(alm2, s, lmax, 2*lmax+1, 2*lmax+1)
    
    difference_abs_squared_map = np.array(np.abs(f_in - f_out)**2, dtype=complex)
    difference_abs_squared_alm = spinsfast.map2salm(difference_abs_squared_map, 0, lmax)
    difference_abs_squared_alm[1:] = 0.0j
    L2norms_of_differences[i] = spinsfast.salm2map(difference_abs_squared_alm, 0, lmax, 2*lmax+1, 2*lmax+1)[0, 0].real

    print("lmax:", lmax,
          "\tMaximum deviation:", np.max(abs(f_in - f_out))/abs(np.max(f_in)),
          "\tL2 norm of deviation:", L2norms_of_differences[i])

plt.semilogy(ells, L2norms_of_differences)
plt.xlabel(r'$\ell_{\mathrm{max}}$')
plt.ylabel(r'$L^2$ norm of difference between input and output');

The convergence is very slow because the function is discontinuous, but it is converging. It might be a little slower than exponential, but it still looks close to me. So I would say that this looks right.

from spinsfast.

Thanks for the extremely thorough response @moble, this makes a lot of sense!

from spinsfast.

As an unrelated note: at one point I was (stupidly) passing in floats instead of complexes, which led to both incorrect and occasionally undefined behavior. I'm happy to put in a PR, but I'm not sure how to do such checks without explicitly writing python wrapper functions.

from spinsfast.

Yeah, wrappers to check inputs are pretty reasonable. I wanted to minimize the changes I made to this code, but over the years I've done enough that there's no reason to avoid this. So the newest version has those two wrappers. Thanks.

from spinsfast.