Reconstruction of clean Sim with True rotations in Single Precision

[garrettwrong]

Making an example for a user, I found that perfect reconstructions are hitting some degenerate condition causing them to crash (example below). This only appears so far as an issue in the most simple case of a perfect reconstruction. Clean images and actual projection rotations as stored in Simulation source.

Using fewer projections (say 200) seems to work. Also adding noise and/or CTF corruptions also seems to avoid.

I checked that this is not a regression from the weighted volume PR. Problem existed well before that. Using doubles in the code below seems to resolve.

>>> import aspire
>>> src = aspire.source.Simulation(L=32, n=1024, offsets=0, amplitudes=1)
2023-07-18 07:53:18,210 INFO [aspire.source.image] Creating Simulation with 1024 images.
>>> bas = aspire.basis.FFBBasis3D(src.L,dtype=src.dtype)
2023-07-18 07:53:46,291 INFO [aspire.basis.ffb_3d] Expanding 3D map in a frequency-domain Fourier–Bessel basis using the fast method.
>>> recon = aspire.reconstruction.MeanEstimator(src,bas).estimate()
2023-07-18 07:54:18,056 INFO [aspire.reconstruction.mean] Determined weighted adjoint mappings. Shape = (1, 8254)
2023-07-18 07:54:18,056 INFO [aspire.reconstruction.estimator] Computing kernel
/Users/garrett/work/ASPIRE-wts/src/aspire/reconstruction/mean.py:270: RuntimeWarning: overflow encountered in square
  sq_filters_f = evaluate_src_filters_on_grid(self.src, _range) ** 2
2023-07-18 07:54:18,154 INFO [aspire.reconstruction.mean] Computing non-centered Fourier Transform
2023-07-18 07:54:18,814 INFO [aspire.reconstruction.estimator] Computing Preconditioner kernel
2023-07-18 07:54:18,814 INFO [aspire.reconstruction.kernel] Circularizing kernel
2023-07-18 07:54:19,339 INFO [aspire.reconstruction.kernel] Circularizing dimension 0
2023-07-18 07:54:19,339 INFO [aspire.reconstruction.kernel] Circularizing dimension 1
2023-07-18 07:54:19,339 INFO [aspire.reconstruction.kernel] Circularizing dimension 2
2023-07-18 07:54:21,168 CRITICAL [root] array must not contain infs or NaNs
Traceback:
  File "<stdin>", line 1, in <module>
  File "/Users/garrett/work/ASPIRE-wts/src/aspire/reconstruction/estimator.py", line 65, in estimate
    est_coeff = self.conj_grad(b_coeff, tol=tol, regularizer=regularizer)
  File "/Users/garrett/work/ASPIRE-wts/src/aspire/reconstruction/mean.py", line 201, in conj_grad
    x, info = cg(operator, b_coeff.flatten(), M=M, callback=cb, tol=tol, atol=0)
  File "<decorator-gen-3>", line 2, in cg
  File "/Users/garrett/miniconda3/envs/test_fftw_src/lib/python3.8/site-packages/scipy/_lib/_threadsafety.py", line 44, in caller
    return func(*a, **kw)
  File "/Users/garrett/miniconda3/envs/test_fftw_src/lib/python3.8/site-packages/scipy/sparse/linalg/_isolve/iterative.py", line 357, in cg
    callback(x)
  File "/Users/garrett/work/ASPIRE-wts/src/aspire/reconstruction/mean.py", line 198, in cb
    f"Delta {norm(b_coeff - self.apply_kernel(xk))} (target {target_residual})"
  File "/Users/garrett/miniconda3/envs/test_fftw_src/lib/python3.8/site-packages/scipy/linalg/_misc.py", line 146, in norm
    a = np.asarray_chkfinite(a)
  File "/Users/garrett/miniconda3/envs/test_fftw_src/lib/python3.8/site-packages/numpy/lib/function_base.py", line 628, in asarray_chkfinite
    raise ValueError(

Error in sys.excepthook:
Traceback (most recent call last):
  File "/Users/garrett/work/ASPIRE-wts/src/aspire/exceptions.py", line 87, in handle_exception
    raise exc_value
  File "<stdin>", line 1, in <module>
  File "/Users/garrett/work/ASPIRE-wts/src/aspire/reconstruction/estimator.py", line 65, in estimate
    est_coeff = self.conj_grad(b_coeff, tol=tol, regularizer=regularizer)
  File "/Users/garrett/work/ASPIRE-wts/src/aspire/reconstruction/mean.py", line 201, in conj_grad
    x, info = cg(operator, b_coeff.flatten(), M=M, callback=cb, tol=tol, atol=0)
  File "<decorator-gen-3>", line 2, in cg
  File "/Users/garrett/miniconda3/envs/test_fftw_src/lib/python3.8/site-packages/scipy/_lib/_threadsafety.py", line 44, in caller
    return func(*a, **kw)
  File "/Users/garrett/miniconda3/envs/test_fftw_src/lib/python3.8/site-packages/scipy/sparse/linalg/_isolve/iterative.py", line 357, in cg
    callback(x)
  File "/Users/garrett/work/ASPIRE-wts/src/aspire/reconstruction/mean.py", line 198, in cb
    f"Delta {norm(b_coeff - self.apply_kernel(xk))} (target {target_residual})"
  File "/Users/garrett/miniconda3/envs/test_fftw_src/lib/python3.8/site-packages/scipy/linalg/_misc.py", line 146, in norm
    a = np.asarray_chkfinite(a)
  File "/Users/garrett/miniconda3/envs/test_fftw_src/lib/python3.8/site-packages/numpy/lib/function_base.py", line 628, in asarray_chkfinite
    raise ValueError(
ValueError: array must not contain infs or NaNs

Original exception was:
Traceback (most recent call last):
  File "/Users/garrett/work/ASPIRE-wts/src/aspire/exceptions.py", line 87, in handle_exception
    raise exc_value
  File "<stdin>", line 1, in <module>
  File "/Users/garrett/work/ASPIRE-wts/src/aspire/reconstruction/estimator.py", line 65, in estimate
    est_coeff = self.conj_grad(b_coeff, tol=tol, regularizer=regularizer)
  File "/Users/garrett/work/ASPIRE-wts/src/aspire/reconstruction/mean.py", line 201, in conj_grad
    x, info = cg(operator, b_coeff.flatten(), M=M, callback=cb, tol=tol, atol=0)
  File "<decorator-gen-3>", line 2, in cg
  File "/Users/garrett/miniconda3/envs/test_fftw_src/lib/python3.8/site-packages/scipy/_lib/_threadsafety.py", line 44, in caller
    return func(*a, **kw)
  File "/Users/garrett/miniconda3/envs/test_fftw_src/lib/python3.8/site-packages/scipy/sparse/linalg/_isolve/iterative.py", line 357, in cg
    callback(x)
  File "/Users/garrett/work/ASPIRE-wts/src/aspire/reconstruction/mean.py", line 198, in cb
    f"Delta {norm(b_coeff - self.apply_kernel(xk))} (target {target_residual})"
  File "/Users/garrett/miniconda3/envs/test_fftw_src/lib/python3.8/site-packages/scipy/linalg/_misc.py", line 146, in norm
    a = np.asarray_chkfinite(a)
  File "/Users/garrett/miniconda3/envs/test_fftw_src/lib/python3.8/site-packages/numpy/lib/function_base.py", line 628, in asarray_chkfinite
    raise ValueError(
ValueError: array must not contain infs or NaNs

I quickly looked at the first warning above, and indeed evaluate_src_filters_on_grid is returning some extremal values that would cause trouble squared in singles. Likely this is just a small bug.

[janden]

That's odd. Just to clarify, when you say “perfect reconstruction”, you are referring to the case of zero noise, known rotations, and a very large number of images?

[garrettwrong]

Yes, taking a Simulation with no noise, no ctf, no offsets, unity amplitudes (and so using the simulated projection rotation). This seems to fail mainly in single precision. The traceback is long, but the sample code above is pretty light to reproduce.

I think its just a small issue in evaluate_src_filters_on_grid. We get back some extremal values, and then in our use, we square them. Some of the values then become inf under singles ~~> crash. It's concerning, but once I noticed where the values were coming from I stopped panicking 😅 .

[janden]

Got it.

[j-c-c]

Resolved by #1076.