`reconstruct3D` PCG backend policy¶

Contract of the code that runs today: the opt-in CTF- and sigma-weighted preconditioned-conjugate-gradient (PCG) reconstruction path. The PCG backend carries no regularization of its own beyond the ordinary FSC/SSNR P_tau replay; nonuniform filtering is the same assembly-owned competition as on the gridding backend (simple_nu_state_filter, policy 2026-09-06). The record of the removed prior experiments -- the binary-envelope solvent prior (2026-08-27) and the direct NU-evidence replay precision Q_NU with its auto-lambda and auto-target controllers (2026-09-06) -- lives in doc/implementation_notes/pcg_priors_history.md. A present-tense four-page description of the backend is doc/implementation_notes/pcg_backend_overview.md and the dated list of decisions with their evidence is doc/implementation_notes/pcg_decision_log.md; new decisions are recorded there, and this policy is updated when a decision changes a contract.

The production reconstruct3D command accepts the selector rec_backend=gridding|pcg, with gridding unchanged as the default. Its pcg branch runs independent kernel accumulation and solves for the even and odd halfsets, writes the standard halfmaps and FSC, and forms the merged map by averaging those two dense halfmaps. It supports both shared-memory execution and distributed worker accumulation followed by a master solve. The former reconstruct3D_pcg command has been retired; there is one production workflow and one backend selector.

The production backend currently supports full-box fixed-pose reconstruction from raw or denoised particle sources, multiple populated states, CTF-aware cc or sigma-weighted euclid, point-group replication, standard halfmap/FSC/merged-map names, project output registration, and the usual final postprocessing handoff. FSC summaries include both the 0.5 and 0.143 resolutions and the cFAR diagnostic, using the same spherical or density-envelope mask as the corresponding FSC. Each selected particle is read once into its state/half accumulation and is never reread during PCG iterations.

In distributed execution, each worker atomically publishes one versioned raw artifact per (state,half,part). It contains the unfolded complex B, real D, particle count, geometry, and provenance. The master validates and adds these artifacts in ascending part-number order. Only after the complete state/half reduction does it fold the RHS, calculate rho floors, finalize Khat, or solve. Empty partitions publish a valid header-only artifact so the association order and completeness check do not depend on particle balance.

Box cropping is supported under the constant-field-of-view contract (box*smpd == box_crop*smpd_crop, enforced at entry). The shared-memory path deliberately rejects projrec=yes, fractional/trailing reconstruction, and conical_fsc=yes regularization; those cases must not silently fall back to gridding or matrix-free PCG. The distributed master integrates the fractional/trailing algebra — raw (B,D) chains blended as (u/f) current + (1-u) previous at full mass, priors applied only after the blend (test=pcg_frac_update is the equivalence gate). The persisted chain is continuous across constant-FOV crop growth: consecutive padded lattices share their frequency step, so the smaller previous grid is an index-aligned central subset and add_raw_accum_weighted embeds it by zero-extension — the PCG analogue of the gridding trailing chain's autoscale ramp. The chain continuation identity (pcg_chain_provenance) therefore excludes the crop fields; a field-of-view change, a shrinking box, or any identity change discards the chain pair and re-seeds through the trailing bootstrap. The one recorded approximation is the old lattice's wrap rim: KB windows that wrapped around the old period leave aliased mass in the outermost old shells, at or beyond the producing stage's matching band and decaying as (1-u)^k — the same approximation the gridding ramp accepts.

1. Production scope and fixed inputs¶

For each populated state and halfset, the backend solves a CTF- and noise-weighted fixed-pose least-squares problem with PCG. Orientations, in-plane shifts, state assignments, half assignments, CTF parameters, phase shifts and ctfflag are read from the project and are never optimized by reconstruction. The Euclidean objective uses per-particle sigma2; correlation is unweighted. mskdiam supplies the spherical fallback support and pgrp is applied by coordinate replication. automsk=yes and automsk=nu replace it with density support; the NU-evidence envelope is never a solve support.

Both execution modes run one route: workers publish raw accumulators and the master reduces, finalizes and solves them. Shared-memory execution runs it in one process, as the only worker and then the master (2026-09-27; the separate in-process solve route and its refusal of fractional and trailing reconstruction are retired). Standard project volumes, halfmaps, FSC, cFAR, resolution fields and postprocessing handoffs are written by the ordinary reconstruct3D workflow. Per-half diagnostics record stop reason, requested/completed iteration count, convergence flag, initial/final true relative residual, final relative update, per-iteration residuals and timings.

Particle images are streamed through a bounded batch. The backend calls begin_accum, repeats accumulate_batch, closes with end_accum, and solves with solve_accum; it never materializes all observed particle planes at once. The two particle-dependent accumulated quantities are the weighted RHS and the |T_i|^2 Gram/sampling-density precursor. The RHS, preconditioner, and kernel are derived from those quantities, so the PCG iteration itself performs no particle-image I/O. Peak image/plane memory is therefore constant in particle count.

Streaming accumulation fuses those two updates. For each particle it evaluates the CTF, rotated coordinate, and KB interpolation window once, then updates both accumulators in the same coloured OpenMP traversal. The monolithic RHS and density routines remain the monolithic test oracle; the streaming-versus- monolithic gate covers both single-batch and multi-batch fused accumulation. It gates B and the D-derived kernel directly at 1e-6 relative error. The 20-step volume comparison has a separate 5e-4 bound because the guarded reciprocal preconditioner and Krylov recurrence amplify single-precision accumulation roundoff; passing the volume gate cannot excuse failure of either direct accumulator gate.

B is complex and D is real. Finalization computes deterministic per-thread rho shell statistics, then produces the reciprocal preconditioner and packed Khat together in one OpenMP pass limited to the reachable Fourier sphere.

2. Where the code lives¶

File	Contents
`src/main/volume/simple_reconstructor_pcg.f90`	`reconstructor_pcg` operator/solver type
`src/main/strategies/parallelization/simple_rec3D_pcg_strategy.f90`	PCG worker accumulation and master reduction and solve, for both execution modes
`src/defs/simple_refine3D_fnames.f90`	raw `(state,half,part)` artifact names
`src/main/commanders/simple/simple_commanders_rec.f90`	`reconstruct3D` commander and backend default
`src/main/ui/simple/simple_ui_refine3D.f90`	`reconstruct3D` UI and backend selector
`src/main/strategies/parallelization/simple_rec3D_strategy.f90`	backend dispatch
`src/main/params/simple_parameters*.f90`	`pcgop`, `rtol` parameters
`src/fileio/simple_sigma2_files.f90`	shared, builder-free sigma2 discovery/loading
`src/main/commanders/test/simple_commanders_test_highlevel.f90`	operator/solver tests

3. Statistical model¶

For volume x and 3-D transform F: G_i extracts an oriented Cartesian Fourier plane by Kaiser-Bessel (KB) interpolation, S_i is the shift phase, C_i the complex CTF (including the phase-flip convention), N_i the diagonal noise covariance from the particle's sigma2. With K_i = N_i^{-1/2} C_i S_i G_i F, the solver targets

H x = b,    H = sum_i K_i^dagger K_i + Lambda,    b = sum_i K_i^dagger N_i^{-1/2} y_i

Lambda = lambda*I is a fixed positive prior, constant for one solve, with the absolute coefficient 1e-3 (PCG_LAMBDA). The relative-lambda CLI (pcg_lambda_rel) was removed as unused; the internal set_lambda_relative mechanism and the deterministic linear fixed-band data scale s_data(D) remain for tests, diagnostics, and future prior anchoring (pcg_priors_history.md). A solve may receive a frozen real-space support. It is fixed before CG starts; no state changes or nonlinear clipping occur during an iteration, so the operator remains linear.

The normal operator carries only the real weight |T_i|^2 = |C_i|^2/sigma2_i. The shift phase is unit-modulus and cancels between the forward and adjoint passes of K_i^dagger K_i. The full complex T_i = C_i S_i / sqrt(sigma2_i) is needed only for b, whose adjoint uses conjg(T_i) — required for shifted particles and phase shifts even when the CTF itself is real.

Output amplitude convention and scale continuity¶

Solved maps are written at the data-quotient convention — the plain weighted least-squares solution, with no box-size scaling — identical to the gridding backend after the legacy division/multiplication pair was retired (doc/implementation_notes/drop_legacy_box_division.md). Deapodization is applied inside the solver; PCG maps must never receive the gridding correction or a second sampling-density correction. Scale continuity is a contract: the euclid/sigma2 equilibrium survives any stable reference amplitude but not an abrupt scale change between consecutive refinement iterations, so any backend handoff, prior, or convention change must preserve the amplitude scale seen by matching (the historical gridding-to-PCG handoff crash was exactly such a jump).

ML two-map contract and starts¶

Refinement solves are two phases from one particle accumulation. The base solve (H_data + lambda I) produces the _unfil half pair; FSC/cFAR and resolution metadata come from that pair. It starts from zero, every iteration: there are NO cross-iteration warm starts (policy 2026-09-10, PfCRT regression gridding_vs_pcg). Warm-starting from the previous iteration's half maps carried an unconverged transient of a slightly different system (new FSC prior, new poses) into the fixed 2-iteration budget, which cannot pull it back; the replay's relative residual then grew iteration over iteration to 10^3-10^4 in the failed runs and never fell below 1 in any run, the shipped ML halves were CG transients, and the base pair, two warm iterations from a stale start, carried no fine-shell evidence into the NU competition at the take-off stage (0% of the mask at 7.96 A where gridding had 6.4%). The fixed small budget from a state-free start is what the simulated-data calibration validated: two iterations beat gridding, and beyond five the residual moves but nothing interpretable in the map does. Every solve line reports INIT= (the relative residual of the start, 1.0 from zero) next to RESID=. A solve from a nonzero start (the replay) that loses positive-definiteness returns stop_reason=indefinite from the solver and is restarted once from zero (solve_with_cold_restart) before the failure is fatal; a solve from zero that loses positive-definiteness fails immediately. Distributed recovery is compute-only inside the even/odd sections; restart reporting and fatal handling occur afterward at the serial finalization boundary. Solver callers that do not request an outcome retain the historical immediate hard failure. The regularized pair deterministically replays kernel finalization from the persisted raw (B,D) with the FSC/SSNR shell-diagonal P_tau installed, and is the closed-form optimum of the diagonal model on that operator (shrink_by_ml_prior, 2026-09-14): the CURRENT same-half base solution's Fourier coefficients on the padded lattice scaled voxelwise by (rho+floor)/(rho+floor+P_tau), i.e. by the ratio of the regularized and base preconditioners. The closed form is shipped as is (maxits_ml=0, the default; ITS=0, STOP=closed_form, RESID/MRES its residuals against the replay system). maxits_ml>0 runs that many coupled PCG iterations of the regularized system FROM the closed form (rtol=0, no start rejection, an indefinite stop falls back to the closed form) and adds a CF MRES a -> b line with the FSC between the start and the solved map inside the band; the sidecar carries closed_form_rel_resid_* and closed_form_vs_solved_*. Tested 2026-09-16 on bgal at box 256, where the closed form's preconditioned residual reads 0.7-1.3 against ~0.1 in the cropped stages: two iterations brought it to 0.3-0.7, changed the map inside the band by FSC >= 0.965 and left every iteration's FSC0.5/0.143, handoff and final map identical to three decimals, so the diagonal model is adequate for the regularized member's jobs and the default is 0. The former replay solve (two CG iterations from a shell-isotropic FSC-shrunk base start, retired 2026-09-14) was prior-dominated over most of Fourier space, its start was rejected on every anisotropically sampled dataset, and its L2 residual above 1 was mostly noise the prior refuses to fit; pcg_decision_log.md. With NU filtering active the base (_unfil) pair seeds the NU candidate bank and the regularized pair joins the competition as the auxiliary member, exactly as on gridding. (The nu_input=gridding|ml alternatives of 2026-09-08 were retired on 2026-09-09; records in doc/implementation_notes/pcg_priors_history.md.) Neither precision nor lambda is ever accumulated into raw B or D.

Every shipped state volume carries a solve-support provenance sidecar (<vol>_pcg_support.txt, solve_support=density|sphere and solve_kind=base|regularized|mixed). The trailing bootstrap follows the same recipe as the gridding bootstrap, applied to the PCG solver's own maps: the FSC prior comes from the lag-one previous shipped pair, and the NU candidate bank is seeded from the CURRENT base pair, volume blended with the previous pair at the applied update weight together with the regularized pair (the gridding trail_restored_halves_if_needed blend). The lag-one pair is never an NU input (fix 2026-09-06). Same recipe and the same kinds of inputs on both backends; the maps differ because the estimators differ. The bootstrap reads the support field for the lag-one FSC pair so the envelope and phase-randomization FSC preprocessing is skipped exactly when that pair was density-constrained in the estimator; a pair without a sidecar is treated as unconstrained, and a bootstrap blend is constrained only if both contributions were. The kind field is provenance only (its former consumer, the base warm-start selector, went with the warm starts). The NU evidence built from a density-constrained pair designates its null on the density envelope's dilation ring (doc/policies/3D/automasking_policy.md).

PCG runs the one static NU competition (2026-09-18: the static ladder [20,15,12,10,8,6,5,4] A capped at fsc/1.5 of the base pair, with the ML-regularized pair as one more member beside the finest retained rung, competing with it at zero prior cost (ml_reg=yes) -- the ed36eb4c abinitio3D machinery, the only NU mechanism since 2026-09-18). Integer Potts coordinates, unit candidate masses, the four fixed evidence bands. The matching handoff is the finest selected label with at least 1% of the signal voxels at that label or finer (nonuniform_filtering_policy.md sections 8 and 12).

Solve support is an automsk feature (policy 2026-09-06). With automsk=no, the default in abinitio3D, every PCG solve, base and regularized replay, runs on the spherical mskdiam support and no density envelope is built. With automsk=yes the conservative density envelope is the support of BOTH the base solve and the regularized replay (policy 2026-09-09; the former envfsc=no split with a spherical base is retired), and envfsc=yes is implied: the FSC pair is envelope-constrained inside the estimator, so no post-hoc mask and no phase-randomized correction are applied to it, and the >>> FSC MODE line in the log and the resolution text says so. This is a deliberate, reported choice, not a claim that a constrained estimate is free of masked-FSC bias. With automsk=nu the same contract holds: the density envelope is the support of both solves, and the NU-evidence envelope is never installed as a solve support (review 2026-09-17). The evidence null of a constrained base pair lies on the density envelope's dilation ring, outside an evidence envelope, so an evidence-supported solve would empty its own null, invalidate the next envelope and fall back to a density envelope derived from a map that is zero outside the evidence support: an oscillating support with no way back for density the evidence excluded (the micelle in the PfCRT collapse, 2026-09-02). Under nu the evidence envelope multiplies the matching references (assembly, matcher fallback) and masks the gridding FSC post hoc. A first solve without a density source bootstraps on the sphere.

Before any reconstruction-derived density source exists, the base necessarily bootstraps on the sphere and its completed pair supplies the conservative replay support. No PCG map is multiplied by either mask after reconstruction. envfsc=yes with automsk=no affects only the phase-randomized FSC evaluation, never a solve. An explicit pcg_mskfile remains the development override and constrains every solve regardless of automsk; it is reported as the state support, so the FSC mode, the support-provenance sidecar and the NU evidence null regime all see a constrained pair.

The original-sampling final reconstructions launched by abinitio3D and refine3D_auto are cold solves. They use a PCG iteration budget of at least five; a larger user-supplied maxits_pcg remains in force. An explicit positive rtol may still stop a converged solve earlier. Ordinary refinement iterations keep the default budget of two iterations from a state-free start (the calibrated regime).

Automatic final-map sharpening estimates its Guinier B-factor from the unregularized half-pair average (_even_unfil/_odd_unfil, no mask of its own since 2026-09-21) between HPLIM_GUINIER = 10 A (2026-09-26, RELION's autob_lowres default; was 20 A) and the FSC=0.143 cutoff. The shipped regularized map is B-sharpened and Butterworth-filtered at the cutoff with no second FSC weighting, since its ML prior already shrank each shell by about its FSC (support-provenance solve_kind=regularized; refine3D policy, final sharpening). The reconstructed map, _lp, _pproc, and mirrored products are not multiplied by any mask; the prohibition on post-hoc PCG masking remains absolute. Iteration-time postprocessing retains its existing behavior.

Beyond-band diagnostic¶

report_beyond_band_excess (module-level in the strategy, both execution paths) compares the RMS of shells beyond the matching band with the band-edge shell and logs >>> PCG BEYOND-BAND EXCESS at ratio >= 10. It is the regression signal for solver defects that park energy above the matched band, where a later stage transition would expose them to euclid matching. The structural mitigation is the closed-form regularized pair (the voxelwise P_tau optimum of the base solution, 2026-09-14), not spectral smoothing (see the removed-experiment record in pcg_priors_history.md).

Backend regression gate¶

test=rec3D_backends reconstructs one fixed particle set with both backends and hard-fails on gated shell-amplitude, FSC, and radial-flatness criteria (band capped at lp); its ground-truth mode adds map-to-truth FSC and radial LS-profile flatness. Mutations restoring the legacy box factor or omitting deapodization must fail it; it is the standing acceptance harness for convention and deapodization changes.

4. Numerical invariants¶

Load-bearing conventions. Changing any of them requires re-passing §9.

Full-symmetric-disk planes. forward_plane/adjoint_plane_add always use an unpacked, both-sign-h disk, never the packed half. This makes them an exact adjoint pair for any orientation, at 2x plane work. The packed half-plane's Nyquist-bin bookkeeping is an orientation-dependent silent energy-loss trap that a single-orientation adjoint test does not catch; the full disk removes that surface permanently. The KB window is used at its natural odd width (2*iwinsz+1); an even width would not be centred on nint(loc) and would break the mirror consistency the fold depends on.

Oversampling. The unknown lives on the native box, but every Fourier operation runs on an OSMPL_PAD_FAC-times padded lattice, centre-padded in and centre-cropped out, pad_vol/crop_vol exact adjoints — the same arrangement gridding uses. Without it KB interpolation is only percent-accurate, the roll-off envelope swings ~30x across the box, and the Gram kernel cannot match the operator. Since the padded lattice is already 2*box, it is the grid the linear (non-wrapping) Gram convolution needs, so kernel and operator share one lattice.

Two KB envelopes, both exact discrete transforms.

Envelope	Origin	Handling
Gather	reading the volume through the KB window	`build_env`, divided out by `deapod_mul`
Deposition	laying `\|T\|^2` onto the kernel grid	divided out in the kernel build

Both use the separable cosine transform of the operator's normalized discrete KB stencil. This is algebraically identical to scattering a unit spike and running the 3-D FFT, without allocating another padded volume. It is not kbinterpol%instr: that continuous transform disagrees with the three discrete renormalized weights by ~2x at the box edge. The matrix-free operator applies the gather envelope twice, so it is A = A_env^{-1}(E T E)A_env^{-1}; deapod_mul brackets it to recover the true T. Not optional for real particles: synthetic data from the operator carries the same envelope and cancels it (an inverse crime), real particles do not, so an uncorrected solve returns E^{-1} x_true.

Soft spherical support. With mskdiam set, the solve is constrained by writing x = P u and minimising ||A P u - y||^2, giving (P H P) u = P b — P on both sides of the operator and once on the RHS. P H P stays symmetric positive semidefinite and its null space is never entered from x = 0. Besides removing edge artefacts where deapodization amplifies hardest, it shrinks the problem: at mskdiam 180 in a 256 box the sphere is ~18% of the volume. Edge profile is production's cosedge_r2_3d.

Scatter reproducibility. Every scatter into a Fourier accumulator is parallelized by h-strided colouring, whose separation guarantee holds on unwrapped coordinates only. The accumulator is periodic, so windows that reach the wrap boundary can collide after folding even in the same colour. All scatter sites therefore split: the parallel colour sweep skips any window win_wraps flags, and a serial pass handles the rim (pre-filtered by sq_rim, a radius below which wrapping is provably impossible). Serial, not atomic — an atomic would fix the race but leave summation order thread-dependent, and bitwise run-to-run reproducibility is what makes a residual trace usable as a regression signal. The rim is a thin outer shell, so the serial cost is small.

5. Preconditioner and kernelized operator¶

Preconditioner. M(k) = rho(k) + floor(shell), with rho = sum_i G_i^dagger |T_i|^2 G_i the sampling density. The floor is a fixed fraction of the shell-mean rho, not an absolute constant: rho spans many orders of magnitude and is genuinely zero between rotated planes, so an absolute floor amplifies the least-constrained modes by six to nine orders of magnitude and PCG fills the map with noise. Modes beyond padf*Rnat are unconstrained and zeroed, making M singular but positive semidefinite, which keeps the Krylov space out of the null space. The envelope belongs in M too: the operator being solved is E^{-1} T E^{-1}, so M^{-1} brackets its Fourier divide with E, not E^{-1}.

Kernelized (Toeplitz/Gram) operator. pcgop=kernel, the default. Replaces the per-iteration particle loop with one padded FFT, a pointwise multiply by a precomputed real Khat, and an inverse FFT — per-iteration cost independent of particle count (~7x faster per iteration). Khat uses the standard NUFFT Gram construction: scatter |T_i|^2 onto the 2x oversampled grid at doubled coordinates. The oversampling is the mechanism, not a safety margin; it resolves the sub-pixel frequency offsets that distinguish the operator from gridding. (A literal impulse-response kernel does not work: KB weights sum to 1, so it reduces exactly to gridding and PCG converges in one step to the gridding answer.) Overall scale is the analytic constant padsc^2; measure_kernel_scale is a test-only check that it is right.

The kernelized operator is shift-invariant and the true operator is not, so it stays an approximation (~3.4% interior error). The matrix-free path is the exact numerical reference for unit tests and small diagnostic fixtures, not a feasible real-workflow backend and not a production fallback. Its per-iteration particle loop becomes prohibitive at experimental particle counts and under symmetry replication.

The approximation is not reliably positive definite under late over-iteration: on the deterministic fixture the residual bottoms out and then curvature fails after the useful map has already saturated. Production therefore defaults to two iterations and rejects maxits>8. Low-level tests may use more iterations when diagnosing the approximation boundary.

pcgop=kernel is therefore the only candidate for production workflows. It must be validated against matrix-free on deterministic, small enough fixtures where both can run. Agreement of a residual trace alone is insufficient: a uniform operator-scale error cancels from the CG recurrence while changing map amplitude. Kernel validation must compare the operator action and scale directly, then compare fixed-iteration solutions built from the same RHS. Real workflow acceptance additionally requires independently reconstructed halfmaps and FSC; a similar operator is not enough.

6. Symmetry¶

Point-group symmetry is applied by coordinate replication inside the operator: each plane pixel is gathered and scattered at all M orientations R_i . S_g. Replication is applied to H (accumulate_absT2, apply_normal_matrixfree) and to b (scatter_plane, reached from apply_adjoint_all and accumulate_rhs_density), so the system solved is consistently symmetrized.

Composition order is production's. matmul(R_i, S_g) in the row-vector convention loc = matmul([h,k,0], rot) is exactly what sym%apply produces and what reconstructor%insert_fplane scatters at. test=pcg_recon stage 9 asserts this by composing the reference path through sym%apply rather than repeating the operator's own expression.

symmats(:,:,1) is exactly the identity for every supported group, so pgrp=c1 (nsym=1) is bit-identical to the pre-symmetry path.

Two consequences to keep in view:

The g loop sits outside the h-strided colour sweep. Inside it, M replicas would write M footprints per colour and separation would have to hold for every pair (g, g'), which it does not.
Cost scales linearly in nsym. Fine for c2/d2; icosahedral (60) multiplies the whole accumulation, including the serial rim pass, by 60. The lattice-exact permutation path — available only for groups whose operators are signed permutations (c2, c4, d2, d4, t, o) — is the way out and is not implemented.

Symmetry-as-constraint and data replication are the same estimator, so there is no SNR difference between them; symmetry's gain is fewer effective unknowns.

7. PCG solver¶

Standard left-preconditioned CG of H x = b (P H P under a support constraint). Zero initialisation. All real-volume dot products use a deterministic double-precision reduction. The solve fails explicitly on non-finite or non-positive dot(p,Hp), a non-positive-definite preconditioner, or a zero RHS, rather than continuing past a broken assumption.

The headline and stopping test are the true relative residual ||r||_2 / ||b||_2, not the preconditioned M-norm — which is not monotone under a singular preconditioner and reads like a diverging solver while the solve is converging. The M-norm is still logged as the diagnostic for the preconditioner: a large gap between the two says M models H poorly. The recurrence residual is periodically audited against a recomputed b - Hx; they agree to six significant figures, which licenses reporting the cheap recurrence norm.

A second criterion, dx/x <= PCG_XTOL, exits on diminishing returns: on noisy real data ||r||/||b|| plateaus above rtol while the map keeps settling, so dx/x is what actually terminates a real solve. It is suppressed by rtol <= 0, the caller's way of demanding exactly maxits iterations — stage 7 depends on that, since comparing two solves stopped by a data-dependent criterion asserts nothing.

Successful completion returns a typed pcg_solver_outcome identifying rtol, xtol, intentional fixed_iterations, or exhausted maxits, with the requested/completed iteration count and initial/final residual/update values. Broken curvature, a non-finite state, a non-positive preconditioner, and a zero RHS are hard errors. Conversion of those failures into recoverable workflow outcomes remains future production work.

The final fixed iteration does not apply the preconditioner because no next search direction will consume it. The terminal residual and update diagnostics remain unchanged; only the unused M-norm is omitted.

8. Sigma2 handling¶

Gated on objfun, matching reconstruct3D: objfun=cc runs unweighted (sigma2 = 1), objfun=euclid requires sigma2 files and weights the fit by them. A missing sigma2 file under objfun=euclid is a hard error, never a silent unweighted fallback — that would quietly change which objective is minimised.

Sigma2 is per-particle-per-shell, read via euclid_sigma2 from the canonical state file (as flex_analysis does), then upsampled to the operator's shell range. Discovery/carry-over/loading lives once in src/fileio/simple_sigma2_files.f90, called by both flex_analysis and the reconstruct3D PCG strategy. That module is deliberately builder-free: builder depends on euclid_sigma2, so the sigma2 side must not depend back on builder; callers pass pftc, esig and orientations explicitly.

9. Tests¶

test=pcg_recon in simple_commanders_test_highlevel.f90 — one gate, nine fail-fast stages, in memory with no project I/O. Each stage gates the ones after it; a visually plausible map is not evidence that the CTF/sigma adjoint is correct.

#	Stage
1	adjoint dot-product identity, `T_i = 1`
2	same with nonzero shift, astigmatic CTF and sigma2 — isolates `build_transfer`
3	normal-operator symmetry and positive-definiteness
4	heterogeneous phantom recovery
5	kernelized-vs-matrix-free operator, scale/energy, and fixed-iteration solution baseline
6	kernel shift-invariance, CTF-dependence, and the preconditioner
7	streaming batches and serialized fixed-order raw reduction reproduce monolithic accumulation
8	deapodization against envelope-free data — the one stage without an inverse crime
9	symmetry replication equals a c1 build of the symmetry-expanded particle set

Stages 1-7 generate observations with forward_plane, so the gather envelope cancels; they gate operator algebra only. Envelope correctness for real particles is gated solely by stage 8.

Matrix-free is the oracle for kernel development. Any change to kernel construction, scale, masks, symmetry, resolution limiting, weights, preconditioning, or box conversion must add or extend a deterministic kernel-versus-matrix-free gate. Production-sized runs exercise kernel only; they do not establish correctness by comparing kernel output with another kernel output.

Not covered, and worth remembering before trusting a change: the backend's own particle I/O loop, and replicated symmetry through the matrix-free operator (stage 9 compares kernel to kernel).

10. Current backend exclusions¶

Not implemented, and hard-errored or absent rather than silently approximated:

no orientation search or pose optimization inside reconstruction;
no online pose update inside reconstruction; fractional/trailing reconstruction is supported through persisted raw accumulator chains, in shared-memory and distributed execution alike;
no post-hoc masking of PCG maps; NU filtering is the assembly-owned competition shared with gridding and writes derived _nu_filt products without touching the primary maps;
no projection-direction compression, conical-FSC regularization, or GPU/offload path;
no reuse of SPIDER BP-CG real-space code — the design is Fourier central-section and the architectures do not transfer. There is no licensing barrier: SIMPLE is GPL-3.0 and SPIDER GPL-2.0-or-later.

11. Reused vs. reimplemented¶

Established data preparation and conventions are reused; the adjoint pair is implemented and tested here rather than repurposed from production gridding.

CTF/sigma physics evaluated directly via ctf/ctfparams (simple_ctf.f90), the same physics image%apply_ctf uses. The per-pixel evaluation is the flat ft_map_ctf_kernel form with per-particle constants hoisted, inlined rather than called: the library routine's memoized (h,k) maps span the h >= 0 half only, while lims2 is a full both-sign-h disk. image%gen_fplane4rec is not reused — it requires a module-global memoised cache (against this path's no-module-global rule) and a packed, padded-box plane convention that clashes with the full-disk native representation.
Particle I/O uses production's batched pattern (prepimgbatch/discrete_read_imgbatch) with norm_noise + taper_edges_particle + fft per particle before plane extraction — the same steps production fuses into norm_noise_taper_edge_pad_fft. Batches stream into the accumulators and are discarded, so peak memory is constant in nptcls.
The adjoint is written fresh, not derived from reconstructor%compress_exp or insert_plane_oversamp: those are production gridding/storage conversions, not established linear adjoints.

12. Execution-path identity and performance rules¶

Backend comparison protocol (review 2026-09-09)¶

A rec_backend=gridding versus rec_backend=pcg comparison in abinitio3D is meaningful because the two paths share everything but the estimator:

stages 1-2 are gridding on both (PCG_REC_START_STAGE=3); the stage ladder, the FSC=0.5 promotion, the uncapped NU-stage handoff, early stopping, the sigma2 handling, the matching references and the final bootstrap_rec3D sequence (gridding bootstrap map, residual sigma pass) are backend-agnostic; only the final map's solver differs;
the worker-side particle preparation is the same (noise normalization against the same mask, edge taper, FFT), the same sigma2 weights and CTF parameters enter both accumulations, and the KB stencil and its deapodization envelope are shared (kb_stencil_centered_crop_inv_envelope_1d);
the ML regularizer is the same formula on both backends: 1/tau2 = <rho>_shell / (tau * fsc/(1-fsc)), FSC clamped to [0.001, 0.999], no prior below hp, driven by the current iteration's unfiltered pair (add_invtausq2rho and build_ml_prior_from_density);
the NU competition, its ladder bound and the handoff run on the unfiltered pair with the regularized pair as the finest member through the one nonuniform_filter_state on both backends; both ship deapodized halves and merged maps carrying the same soft spherical support at msk_crop (the PCG solve support; the gridding restoration applies the identical mask3D_soft after deapodization, 2026-09-09) and compute the FSC on those halves through evaluate_halfmap_pair, which applies no mask of its own;
the master phase gets the same thread budget on local execution (rec3D_master_nthr: nparts x nthr capped at 32, the PCG rule, now also applied to the gridding volassemble instead of NTHR_SHMEM_MAX), and both log one RECONSTRUCTION MASTER PHASE (<backend>): <s> line per iteration; the worker side is in the bench files (partial reconstruction), memory in the peak-RSS fields.

Differences that are the estimator itself and belong in the comparison: PCG solves (H + lambda) x = b with two CG iterations per half (base, from zero) and takes the regularized pair in closed form from it (2026-09-14), a fresh estimate every iteration like the gridding density quotient (cross-iteration warm starts retired 2026-09-10).

The former measurement asymmetry (gridding FSC on the apodized halves for legacy parity, PCG on the solved halves) was removed on 2026-09-09: the gridding path now computes its FSC on the deapodized, support-masked halves it ships. The one-mask contract that came with it:

every reconstruction product (halves, _unfil halves, merged map) carries the soft spherical support at msk_crop exactly once, installed by the estimator (PCG) or by the restoration after deapodization (gridding), and recorded in the support-provenance sidecar <vol>_pcg_support.txt (solve_kind=gridding for gridding products, gridding_regularized when ml_reg=yes put the FSC-derived prior into the restoration; since 2026-09-26 postprocess reads the kind to skip its FSC weighting on ML-regularized maps);
evaluate_halfmap_pair masks nothing of its own (the envfsc envelope + phase-randomization correction is applied only to an unconstrained pair, and the mode actually used is logged as >>> FSC MODE);
postprocess applies no post-hoc mask to a volume carrying the sidecar (previously PCG-only by backend name; an imported map without the sidecar still gets the classical spherical/envelope mask);
the PCG support is a hard solve domain plus one soft window (review 2026-09-09, finding 4.2): set_mask builds the window with mask3D_soft, the solve runs on the domain window > 0 (P^2 = P, exact projections in operator, RHS and preconditioner) and the shipped map is window * u. That is the same "estimate times one soft window" the gridding restoration ships, so the two backends' band treatment is identical and the only estimator difference is that the PCG estimate is zero outside the domain. The soft P H P formulation was not equivalent: where 0 < P < 1 the solved variable compensates for P, so the band was a solver-state dependent mixture (PCG_HARD_SOLVE_SUPPORT in the solver restores it for experiments). Output-space maps entering the solver (nonzero starts, and the base map the closed-form regularized pair is derived from) are never re-masked: the entry converts them back to u = x / window where window >= PCG_SUPPORT_DIV_MIN (zero below) instead of projecting again; the former entry projection squared the edge on every restarted iteration and compounded over a stage. Regressions: test=pcg_recon stage 14 (band profile of the constrained solve against the windowed unconstrained one; band stability under repeated starts);
the matcher still applies mask3D_soft(msk_crop) to its reprojection reference after Fourier filtering (both backends, mask_matching_reference). That is reference preparation, not an estimate: it restores compact support after the filter's ringing and covers user-supplied start volumes. In the 12 px cosine band the reference therefore carries P^2; it is the one remaining second application and is deliberate.

Comparison contract (review 2026-09-09)¶

Two inputs that could change the result independently of the backend are now shared, and the measurement records are complete enough to separate wall time, compute cost and memory:

one observation preparation. Cropping and edge tapering do not commute. Gridding normalized at the native box, Fourier-cropped, then tapered at the crop box; PCG tapered at the native box and cropped afterwards. Both now go through prep_rec_observation (simple_matcher_ptcl_io): normalize at the native box, Fourier-crop, taper at the cropped box; gridding pads and transforms in its fused routine, PCG takes the native plane of the tapered observation. Without a crop the observation is tapered first and normalized second on both backends. test=pcg_recon stage 13 gates the two routes in real space at 1e-5 relative;
the volume and its support sidecar are one artifact. Stage-boundary renames, final copies, symmetric-map copies and the refine3D start-volume rename move or copy the sidecar with the map (copy/rename_support_provenance); noise starts and imported maps remove a stale sidecar (remove_support_provenance); the PCG masters publish the map first and the sidecar second. Before this a PCG half lost its provenance at every stage boundary and the next stage's first base solve was an unintended cold start; a copied final gridding map could be masked a second time by a standalone postprocess;
cost records. Every partition writes REFINE3D_BENCH_ITERnnn_PARTppp.txt (partition 1 also the legacy file) with nparts, worker threads, box, box_crop, backend, maxits_pcg, rtol, peak and phase RSS, the phase timings and their thread-seconds; the master logs one RECONSTRUCTION MASTER PHASE (<backend>): s on n threads = thread-s; master peak RSS line per iteration. Load imbalance is the spread of the per-part partial-reconstruction times; core-hours are the sum of worker thread-seconds plus master thread-seconds.

Baseline profile for the first quality/cost study: automsk=no, envfsc=no, conical_fsc=no, no pcg_mskfile; identical starting project, maps, seed, even/odd assignment, sampling and sigma state, stage schedule and resolution limits, nparts, threads, queue and hardware between the paired runs; each backend in its own clone of the starting project; maxits_pcg, rtol and the operator mode fixed and recorded. Gates before interpreting results: prepared-observation parity (stage 13), the support regressions (stage 14), no reconstruction losing provenance across a stage boundary, and wall time, thread-seconds and peak RSS recomputable from the emitted records. Judge quality on paired replicas (at least 10 seeds), never on a single run.

Shared-memory and distributed execution are two parallelizations of one algorithm. Output conventions, warm starts, and diagnostics are implemented once at module level and invoked identically from both entry points; nothing may depend on which path produced a map.

Durable performance rules (from the retired production-readiness note):

one logical particle read per accumulation phase (direct source reads for standalone reconstruct3D; the validated downscaled cache for cache-enabled refinement); particle residency bounded by MAXIMGBATCHSZ; no particle plane cache — kernel iterations are data-free after (B,D);
the preconditioner uses the padded Toeplitz lattice; do not trade the pad/crop geometry for a native-grid FFT optimization;
process one state/half at a time; do not overlap the largest accumulation and solve scratch allocations; keep D real; keep the fused reciprocal/Khat packing over the reachable Fourier sphere and the exact separable discrete KB-envelope construction;
kernel failure never retries with matrix-free (too slow, hides defects, impossible after distributed reduction); fail structurally instead;
symmetry cost scales with group order (coordinate replication); validate C1/C2/D2 and treat high-order groups as a measured-budget contract; the lattice-exact permutation path for signed-permutation groups is the eventual optimization and is not implemented;
use the per-phase timing diagnostics, not total wall time, to choose the next optimization.

reconstruct3D PCG backend policy¶