Projection/Backprojection Adjointedness

src/aspire/covariance/covar.py

@@ -57,7 +57,7 @@ def compute_kernel(self):
                 weights[:, 0, :] = 0
 
             # TODO: This is where this differs from MeanEstimator
-            pts_rot = np.moveaxis(pts_rot[::-1], 1, 0).reshape(-1, 3, L**2)
+            pts_rot = np.moveaxis(pts_rot, 1, 0).reshape(-1, 3, L**2)
             weights = weights.T.reshape((-1, L**2))
 
             batch_n = weights.shape[0]
@@ -67,7 +67,7 @@ def compute_kernel(self):
                 factors[j] = anufft(weights[j], pts_rot[j], (_2L, _2L, _2L), real=True)
 
             factors = Volume(factors).to_vec()
-            kernel += vecmat_to_volmat(factors.T @ factors) / (n * L**8)
+            kernel += (factors.T @ factors).reshape((_2L,) * 6) / (n * L**8)
 
         # Ensure symmetric kernel
         kernel[0, :, :, :, :, :] = 0
@@ -90,6 +90,8 @@ def estimate(self, mean_vol, noise_variance, tol=1e-5, regularizer=0):
         b_coef = self.src_backward(mean_vol, noise_variance)
         est_coef = self.conj_grad(b_coef, tol=tol, regularizer=regularizer)
         covar_est = self.basis.mat_evaluate(est_coef)
+        # Note, notice these cancel out, but requires a lot of changes elsewhere in this file,
+        # basically totally removing all the `utils/matrix` hacks ... todo.
         covar_est = vecmat_to_volmat(make_symmat(volmat_to_vecmat(covar_est)))
         return covar_est
 
@@ -180,7 +182,9 @@ def src_backward(self, mean_vol, noise_variance, shrink_method=None):
                 im_centered_b[j] = self.src.im_backward(im_centered[j], i + j)
             im_centered_b = Volume(im_centered_b).to_vec()
 
-            covar_b += vecmat_to_volmat(im_centered_b.T @ im_centered_b) / self.src.n
+            covar_b += (im_centered_b.T @ im_centered_b).reshape(
+                (self.src.L,) * 6
+            ) / self.src.n
 
         covar_b_coef = self.basis.mat_evaluate_t(covar_b)
         return self._shrink(covar_b_coef, noise_variance, shrink_method)

src/aspire/image/image.py

@@ -537,7 +537,7 @@ def backproject(self, rot_matrices, symmetry_group=None):
             )
             pts_rot = pts_rot.reshape((3, -1))
 
-            vol += anufft(im_f, pts_rot[::-1], (L, L, L), real=True)
+            vol += anufft(im_f, pts_rot, (L, L, L), real=True)
 
         vol /= L
 

src/aspire/reconstruction/estimator.py

@@ -56,6 +56,13 @@ def __init__(
         self.basis = basis
         self.dtype = self.src.dtype
         self.batch_size = batch_size
+        if not preconditioner or preconditioner.lower() == "none":
+            # Resolve None and string nones to None
+            preconditioner = None
+        elif preconditioner not in ["circulant"]:
+            raise ValueError(
+                f"Supplied preconditioner {preconditioner} is not supported."
+            )
         self.preconditioner = preconditioner
         self.boost = boost
 
@@ -128,12 +135,12 @@ def __getattr__(self, name):
     def compute_kernel(self):
         raise NotImplementedError("Subclasses must implement the compute_kernel method")
 
-    def estimate(self, b_coef=None, tol=1e-5, regularizer=0):
+    def estimate(self, b_coef=None, x0=None, tol=1e-5, regularizer=0):
         """Return an estimate as a Volume instance."""
         if b_coef is None:
             b_coef = self.src_backward()
-        est_coef = self.conj_grad(b_coef, tol=tol, regularizer=regularizer)
-        est = Coef(self.basis, est_coef).evaluate().T
+        est_coef = self.conj_grad(b_coef, x0=x0, tol=tol, regularizer=regularizer)
+        est = Coef(self.basis, est_coef).evaluate()
 
         return est
 

src/aspire/reconstruction/mean.py

@@ -64,7 +64,9 @@ def __getattr__(self, name):
                     1.0 / self.kernel.circularize()
                 )
             else:
-                if self.preconditioner.lower() not in (None, "none"):
+                if self.preconditioner and (
+                    self.preconditioner.lower() not in ("none")
+                ):
                     logger.warning(
                         f"Preconditioner {self.preconditioner} is not implemented, resetting to default of None."
                     )
@@ -126,7 +128,7 @@ def _compute_kernel(self):
                         batch_kernel += (
                             1
                             / (self.r * self.src.L**4)
-                            * anufft(weights, pts_rot[::-1], (_2L, _2L, _2L), real=True)
+                            * anufft(weights, pts_rot, (_2L, _2L, _2L), real=True)
                         )
 
                     kernel[k, j] += batch_kernel
@@ -189,7 +191,7 @@ def src_backward(self):
 
         return res
 
-    def conj_grad(self, b_coef, tol=1e-5, regularizer=0):
+    def conj_grad(self, b_coef, x0=None, tol=1e-5, regularizer=0):
         count = b_coef.shape[-1]  # b_coef should be (r, basis.count)
         kernel = self.kernel
 
@@ -240,12 +242,13 @@ def cb(xk):
                 # Construct checkpoint path
                 path = f"{self.checkpoint_prefix}_iter{self.i:04d}.npy"
                 # Write out the current solution
-                np.save(path, xk.reshape(self.r, self.basis.count))
+                np.save(path, xk)
                 logger.info(f"Checkpoint saved to `{path}`")
 
         x, info = cg(
             operator,
             b_coef.flatten(),
+            x0=x0,
             M=M,
             callback=cb,
             tol=tol,