Ctf3d ci

[bHimes]

putting claude code review to the test

[claude]

🔴 The CUDA grid/block calculation on line 102 is inverted: blocks = (threads / *n_elements + 1024 - 1) / threads computes 0 blocks for any realistic image size (n_elements > 1024) due to integer division truncation. Additionally, the kernel launches at lines 124 and 168 swap grid and block dimensions (<<<threads, blocks>>> instead of <<<blocks, threads>>>). This means both convert_fp32_to_fp16 and convert_fp16_to_fp32 GPU kernels silently produce no output, causing the emc_halfcast GPU path to return all-zeros and breaking template matching when test_half is enabled.

Extended reasoning...

The Bug

The block count formula on line 102 of mexFP16.cu is:

const size_t threads = 1024;
const size_t blocks = (threads / *n_elements + 1024 - 1) / threads;

This divides threads by *n_elements instead of the other way around. The standard CUDA ceiling-division idiom should be (*n_elements + threads - 1) / threads. Because threads = 1024 and any realistic cryo-EM image has millions of elements, the integer division 1024 / n_elements truncates to 0, making blocks = (0 + 1023) / 1024 = 0.

Swapped Launch Parameters

Beyond the formula bug, the kernel launches at lines 124 and 168 pass the variables in the wrong order:

convert_fp32_to_fp16<<<threads, blocks, 0, cudaStreamPerThread>>>(...);  // line 124
convert_fp16_to_fp32<<<threads, blocks, 0, cudaStreamPerThread>>>(...);  // line 168

CUDA syntax is <<<gridDim, blockDim>>>. Here threads (1024) is used as gridDim and blocks (0) as blockDim. Even if the formula were correct, the dimensions would be swapped.

Step-by-Step Proof

Consider a typical cryo-EM tomogram slice of 4000×4000 pixels, so n_elements = 16,000,000:

1. threads = 1024
2. blocks = (1024 / 16000000 + 1024 - 1) / 1024
3. 1024 / 16000000 = 0 (integer division)
4. (0 + 1023) / 1024 = 0
5. Kernel launch: <<<1024, 0>>> — blockDim of 0 means zero threads execute
6. The output buffer is never written, so it retains its uninitialized/zero contents

This affects both the FP32→FP16 path (line 124) and the FP16→FP32 path (line 168).

Impact

In BH_templateSearch3d_2.m, emc_halfcast is called at line 453 to store tomogram data in half-precision on the GPU, and at line 530 to convert it back to single-precision for template matching. When both arrays are on the GPU, the buggy kernel path is taken. The conversion silently produces all-zeros, meaning tomoFou at line 530 is a zero array, making the entire template search produce meaningless cross-correlation results. The CPU fallback path using half_float::half (lines 147-149 and 187-189) works correctly, so this bug only manifests when both arrays reside on the GPU.

Fix

Replace line 102 with the correct ceiling-division formula and swap the launch parameters:

const size_t threads = 1024;
const size_t blocks = (*n_elements + threads - 1) / threads;

And change both kernel launches to <<<blocks, threads, 0, cudaStreamPerThread>>>.

[claude]

🔴 PYTHON_STYLE_GUIDE.md references Black, isort, flake8, mypy, and autopep8 throughout, but the project actually uses Ruff and Pyright as correctly documented in CLAUDE.md, pyproject.toml, .pre-commit-config.yaml, and the CI workflows. Developers following this guide will install the wrong tools and get different behavior than CI.

Extended reasoning...

Wrong Tooling Referenced in PYTHON_STYLE_GUIDE.md

The newly added PYTHON_STYLE_GUIDE.md consistently references Black, isort, flake8, mypy, and autopep8 as the project's Python tooling. However, the actual project tooling—as correctly configured elsewhere in this same PR—is Ruff (for linting and formatting, replacing black, isort, and flake8) and Pyright (for type checking, replacing mypy).

Evidence of the correct tooling

• CLAUDE.md (lines ~219-224) states: "Linting/Formatting: ruff (replaces black, isort, flake8)" and "Type checking: pyright"
• .pre-commit-config.yaml uses ruff-pre-commit (v0.12.11) for linting/formatting and pyright-python (v1.1.405) for type checking
• .github/workflows/code-style.yml installs and runs ruff check and ruff format
• .github/workflows/type-checking.yml installs and runs pyright
• pyproject.toml lists ruff>=0.1.0 and pyright>=1.1.330 as dev dependencies

Scope of incorrect references in PYTHON_STYLE_GUIDE.md

The wrong tools are referenced pervasively throughout the file:

• Overview (lines 7-11): Lists Black, isort, flake8, mypy
• Install instructions (line 17): pip install black isort flake8 mypy autopep8
• Pre-commit hook (lines 22-37): Uses black ., isort ., flake8 .
• Auto-fix commands (lines 232-239): Uses isort ., black ., autopep8
• IDE configuration (lines 249-265): Configures VS Code for black/isort/flake8/mypy
• CI Integration (lines 276-285): Shows black --check, isort --check-only, flake8, mypy
• Quick Reference (lines 289-293): Uses black . && isort . && flake8 .
• Resources (lines 309-315): Links to black/isort/flake8/mypy documentation

Step-by-step proof of impact

1. A new developer reads PYTHON_STYLE_GUIDE.md and runs pip install black isort flake8 mypy autopep8 (line 17)
2. They format their code with black . and isort . as instructed (lines 289-293)
3. They push their changes, and CI runs ruff check and ruff format --check instead
4. Ruff and Black have slightly different formatting rules (e.g., magic trailing commas, quote styles), so CI may fail even though the developer followed the style guide exactly
5. Similarly, mypy and pyright have different type checking behaviors, so type errors flagged locally may differ from CI

Recommended fix

Rewrite PYTHON_STYLE_GUIDE.md to reference Ruff and Pyright throughout, matching the actual tooling configured in pyproject.toml, .pre-commit-config.yaml, and the CI workflows. Alternatively, since CLAUDE.md already documents the correct tooling, consider whether this separate style guide is needed at all—having two documents increases the risk of them diverging.

[claude]

🔴 CUM_e_DOSE is hardcoded to 0 at line 99, replacing the original pBH.("CUM_e_DOSE") lookup without wiring it into the new emc struct. This silently disables dose weighting in the fallback dose path (lines 252-259), where exposure = CUM_e_DOSE./nPrjs evaluates to zero regardless of the value users set in their parameter files (e.g., CUM_e_DOSE=60). The fix is to add CUM_e_DOSE to the emc struct in BH_parseParameterFile.m and read it back here, or if the old dose path is being intentionally deprecated, remove the dead code.

Extended reasoning...

Bug Analysis: CUM_e_DOSE hardcoded to 0 disables dose weighting

During the migration from the pBH parameter struct to the new emc struct, the line CUM_e_DOSE = pBH.("CUM_e_DOSE"); was replaced with CUM_e_DOSE = 0; at line 99 of ctf/BH_ctf_Estimate.m. However, CUM_e_DOSE was never added to the new emc struct in BH_parseParameterFile.m, so there is no way for the user-specified value to reach this code.

How the bug manifests

The fallback dose calculation path is activated when flgOldDose && flgStandardOrdeDoCalc is true (lines 251-259). In this path, CUM_e_DOSE determines the per-tilt exposure:

if CUM_e_DOSE < 0
  flgCosineDose = 1;
  exposure = abs(CUM_e_DOSE);
else
  flgCosineDose = 0;
  exposure = CUM_e_DOSE./nPrjs;
end

With CUM_e_DOSE = 0, the condition CUM_e_DOSE < 0 is false, so exposure = 0 / nPrjs = 0. Zero exposure means no dose weighting is applied to the CTF estimation.

Step-by-step proof

1. User sets CUM_e_DOSE = 60 in their parameter file (as shown in example .param files).
2. BH_parseParameterFile.m parses the file but does not store CUM_e_DOSE in the emc struct.
3. BH_ctf_Estimate.m line 99 sets CUM_e_DOSE = 0 (hardcoded).
4. The fallback dose path is entered (when flgOldDose=1 and flgStandardOrdeDoCalc=1).
5. exposure = 0 / nPrjs = 0 — dose weighting is silently disabled.
6. The user gets CTF estimates without dose weighting, with no warning or error.

Impact

This is a silent correctness issue. Users who rely on the fallback dose path (which is the active path when newer dose parameters are at their defaults) will get CTF estimates computed without dose weighting. The user-configured CUM_e_DOSE value is simply ignored. Since there is no error or warning, users have no indication that dose weighting is not being applied, which can degrade reconstruction quality.

Suggested fix

Add CUM_e_DOSE to the parsed parameters in BH_parseParameterFile.m (following the same pattern as other migrated parameters), and replace CUM_e_DOSE = 0; with a read from the emc struct. Alternatively, if the old dose calculation path is intentionally being deprecated in favor of the new dose mechanism, the dead code in lines 251-259 should be removed and the hardcoded 0 should be documented as intentional.

[claude]

🟡 firstLoopOverAngles = false at line 689 is inside the if (measure_noise_variance) block, but measure_noise_variance defaults to false (line 128), so firstLoopOverAngles is never set to false during normal operation. This causes ANGLE_LIST to be redundantly rebuilt on every tomogram chunk iteration. Move firstLoopOverAngles = false; after the if (measure_noise_variance)...end block (i.e., between the current lines 690 and 691).

Extended reasoning...

What the bug is

In BH_templateSearch3d_2.m, the variable firstLoopOverAngles is initialized to true at line 493 and is used as a guard at line 569 to control whether ANGLE_LIST is constructed during the angle search loops. The intent is for ANGLE_LIST to be built only once (on the first tomogram chunk iteration) and then reused for all subsequent chunks. However, the statement firstLoopOverAngles = false; at line 689 is scoped inside the if (measure_noise_variance) block (lines 661-690).

Why existing code does not prevent it

measure_noise_variance defaults to false at line 128 (measure_noise_variance = false). Since the if (measure_noise_variance) block is never entered during normal operation, firstLoopOverAngles remains true for the entire execution. There is no other location in the code that sets firstLoopOverAngles = false.

Step-by-step proof

1. Line 128: measure_noise_variance = false; (the default, unless overridden by the parameter file)
2. Line 493: firstLoopOverAngles = true; (initialization before the main loop)
3. Line 569: if (firstLoopOverAngles) gates the construction of ANGLE_LIST(currentGlobalAngle,:) = [phi, theta, psi - phi];
4. Lines 661-690: The if (measure_noise_variance) block contains firstLoopOverAngles = false; at line 689
5. Since measure_noise_variance is false, the block is skipped, and firstLoopOverAngles stays true
6. On the next tomogram chunk iteration, the if (firstLoopOverAngles) guard at line 569 is still true, so ANGLE_LIST is rebuilt identically

Comparison with old code

In the old version of BH_templateSearch3d_2.m, firstLoopOverAngle = false (note the slightly different variable name) was unconditionally executed at the end of the chunk loop body (visible in the diff context around line 1222 of the old code). The refactoring in this PR moved the loop structure so that the tomo chunks are now the outer loop and angles are the inner loop, but in doing so, the firstLoopOverAngles = false statement was accidentally placed inside the if (measure_noise_variance) conditional instead of after it.

Impact

The practical impact is a performance issue, not a correctness issue. ANGLE_LIST receives identical values each time it is rebuilt because the angle search parameters are deterministic. The redundant work includes re-evaluating the angle loop conditionals and re-writing the same values into ANGLE_LIST on every tomogram chunk iteration. For large datasets with many chunks, this adds unnecessary computation.

Fix

Move firstLoopOverAngles = false; from line 689 (inside the if (measure_noise_variance) block) to after line 690 (after the end that closes the if block but before the end that closes the for iTomo loop). This ensures it is unconditionally set to false after the first complete pass over all angles, regardless of the measure_noise_variance setting.

[claude]

🟡 Bounds checks at lines 66-77 still use old-format Y-position logic (yMin/yEnd) but the coordinate format changed to (NX, NY, NZ, dX, dY, dZ). Now read_in_Coords(2) is NY (a size, always positive) so the < -75 check is dead code, and read_in_Coords(3) is NZ (Z thickness) which gets incorrectly clamped against yMax (the Y image dimension) — a semantically wrong comparison that could silently truncate Z thickness for thick reconstructions. The stale comment at lines 9-14 also still describes the old format. The fix is to remove these old bounds checks (correct positive-size validation already exists at lines 96-104) and update the comment.

Extended reasoning...

What the bug is

The function BH_multi_recGeom.m reads reconstruction coordinate files and builds geometry structs for each tomogram. The coordinate format was changed from the old layout (NX, yStart, yEnd, NZ, xShift, zShift) to the new layout (NX, NY, NZ, dX, dY, dZ). The struct assignment at lines 80-87 correctly uses the new format, but the bounds-checking block at lines 63-77 was never updated and still applies old-format Y-position logic.

How it manifests

There are two specific issues in the bounds-checking code:

1. Line 66: read_in_Coords(2) < -75 — In the new format, read_in_Coords(2) is NY, a positive dimension size. A size value can never be negative, so this check is dead code. The follow-up clamp at line 70 (read_in_Coords(2) = 1 if it is less than 1) is also unreachable under normal conditions since NY should always be a positive integer.

2. Lines 72-77: read_in_Coords(3) > yMax + 75 and the clamp read_in_Coords(3) = yMax — In the new format, read_in_Coords(3) is NZ (the Z thickness of the reconstruction), but it is being compared against yMax = tilt_geometry(1,21) (the Y dimension of the tilt image). These are semantically unrelated quantities. If a user requests a Z thickness larger than the tilt image Y dimension, the code would silently clamp NZ to yMax, producing a reconstruction with incorrect Z extent.

Step-by-step proof

Consider a tilt series where tilt_geometry(1,21) = 4000 (Y dimension) and a reconstruction coordinate file specifying NX=3000, NY=2000, NZ=5000, dX=0, dY=0, dZ=0. Walking through the code:

• read_in_Coords = [3000, 2000, 5000, 0, 0, 0]
• Line 66: read_in_Coords(2) < -75 → 2000 < -75 → false (harmless)
• Line 72: yMax = 4000
• Line 73: read_in_Coords(3) > yMax + 75 → 5000 > 4075 → true → error thrown, citing "ymin" and "ymax" in the error message even though the value is NZ
• If the value were, say, 4050 instead: the error would not trigger, but line 75-76 would clamp it: read_in_Coords(3) = 4000, silently reducing NZ from 4050 to 4000 with no warning

This means thick reconstructions (NZ > tilt Y dimension) would either error out with a misleading message or be silently truncated.

Why existing code does not prevent it

The correct positive-size validation at lines 96-104 checks that NX, NY, NZ are all positive, which is the appropriate validation for dimension sizes. However, the old bounds-checking block at lines 63-77 runs first and can modify read_in_Coords(3) (NZ) before it reaches the struct assignment, corrupting the value.

How to fix

Remove the old bounds-checking block at lines 63-77 entirely (the correct validation at lines 96-104 already covers the meaningful case). Update the stale comment at lines 9-14 to reflect the new coordinate format (NX, NY, NZ, dX, dY, dZ). If Y-range bounds checking is still desired, it should operate on the computed y_i and y_f values at lines 90-93, which are the actual Y-range coordinates derived from NY and dY.

[claude]

🟡 Line 838: fprint(...) should be fprintf(...). The old code had print(...) which was also wrong; this PR changed one typo to another. If gather() fails in this catch block, MATLAB will throw an "Undefined function fprint" error instead of printing the diagnostic message, and the break will never execute.

Extended reasoning...

What the bug is

In alignment/BH_templateSearch3d_2.m at line 838, there is a catch block that attempts to print a diagnostic message when gather([i,j,k]) fails (e.g., due to GPU memory issues). The code uses fprint(...), which is not a valid MATLAB built-in function. The correct function is fprintf(...).

How it manifests

The old code (visible in the diff for BH_templateSearch3d.m, the predecessor file) used print(...), which was also not a valid MATLAB function. This PR reorganized the code into BH_templateSearch3d_2.m and changed print to fprint, but the correct call is fprintf. When gather() throws an exception and execution enters the catch block, MATLAB will immediately throw a new error: "Undefined function or variable 'fprint'". This means:

1. The intended diagnostic message ("Ran into some trouble gathering the i,j,k. Breaking out") is never displayed
2. The break statement on line 839 never executes
3. Instead of a graceful exit from the loop, an unhandled error crashes the template search

Step-by-step proof

1. During template matching, the while loop at ~line 803 iterates over peaks
2. At line 833, [i,j,k] = ind2sub(sizeTomo, coord) computes indices
3. At line 835, c = gather([i,j,k]) attempts to transfer GPU data to CPU
4. If gather() fails (e.g., GPU out of memory), execution enters the catch block at line 837
5. Line 838 calls fprint("Ran into some trouble...") — but fprint is undefined
6. MATLAB throws: Undefined function 'fprint' for input arguments of type 'char'
7. The break on line 839 never executes, and the error propagates up the call stack

Impact

This is in a rare error-handling path (GPU memory failure during gather), so the practical impact is low. However, when it does trigger, the user gets an unhelpful "Undefined function fprint" error instead of the intended diagnostic message and graceful loop exit.

Fix

Change fprint to fprintf on line 838:

fprintf("Ran into some trouble gathering the i,j,k. Breaking out\n");

This is a pre-existing issue (the old code also had the wrong function name), but since this PR directly modified the line, it had the opportunity to fix it and instead introduced a slightly different typo.

[claude]

🟣 Bug: save.m line 57 uses typecast(mRCImage.volume, 'uint16') to convert the volume before writing mode-12 (half-precision) files, but after loadVolume the volume is already single-precision (4 bytes/element). typecast reinterprets raw bytes rather than converting numerically, so each 4-byte single becomes two uint16 values — doubling the element count and writing completely corrupt data. The fix is to replace typecast(...) with emc_halfcast(mRCImage.volume), which properly converts single→FP16 stored as uint16, matching what loadVolume does in reverse.

Extended reasoning...

What the bug is

In @MRCImage/save.m, when saving mode-12 (half-precision / FP16) MRC files, line 57 uses typecast(mRCImage.volume, 'uint16') to prepare the volume data for writing via fwrite. However, typecast performs raw byte reinterpretation — it does not do numeric conversion. If the volume contains single-precision floats (4 bytes each), typecasting to uint16 (2 bytes each) reinterprets each single as two uint16 values, doubling the element count and producing completely meaningless data.

The specific code path that triggers it

1. Load: loadVolume.m reads mode-12 data as uint16 via fread (lines 36-38, since fread cannot handle half), then converts to single-precision via emc_halfcast(mRCImage.volume) at line 72. The header mode remains 12 (there is even a FIXME comment at line 70 acknowledging this).
2. Save: save.m calls getModeString() which returns 'half' for mode 12 (line 50). Lines 51-57 detect the 'half' mode string and substitute 'uint16', then call typecast(mRCImage.volume, 'uint16'). Since the volume is now single-precision (from step 1), this reinterprets 4-byte floats as pairs of 2-byte uint16 values.
3. Write: fwrite at line 101 writes the corrupted, doubled-length array to disk.

Step-by-step proof with a concrete example

Consider a mode-12 MRC file with a single voxel whose FP16 value is 1.0 (stored on disk as uint16 value 0x3C00, i.e. 2 bytes):

1. loadVolume reads it as uint16 0x3C00 via fread.
2. emc_halfcast converts uint16 0x3C00 → single 1.0 (stored as 4 bytes: 0x3F800000).
3. Header mode stays 12. Volume is now single 1.0.
4. save.m sees mode 12, calls typecast(single(1.0), 'uint16'). This reinterprets the 4 bytes 0x3F800000 as two uint16 values: [0x0000, 0x3F80] (on little-endian). That is 2 elements instead of 1, and neither represents 1.0 in FP16.
5. fwrite writes 4 bytes (2 uint16s) for what should be a single 2-byte FP16 value. The file is now corrupt — wrong size and wrong values.

Why existing code does not prevent it

The count ~= nElements check at line 102 compares against numel(mRCImage.volume), but after typecast the volume array has already been doubled in size. So numel returns 2N and fwrite also writes 2N elements — the check passes despite the corruption. The FIXME comment in loadVolume.m line 70 (% FIXME: should the header mode be changed?) acknowledges the mode/data mismatch but was never addressed.

Impact

Any round-trip of mode-12 (half-precision) MRC data through loadVolume followed by save will produce a corrupt output file with doubled size and scrambled voxel values. This is silent data corruption — no error is raised.

How to fix

Replace line 57 in save.m:

mRCImage.volume = typecast(mRCImage.volume, 'uint16');

with:

mRCImage.volume = emc_halfcast(mRCImage.volume);

emc_halfcast properly converts single→FP16 (stored as uint16), which is the exact inverse of what loadVolume does at line 72. This is already used at line 114 of save.m for reconverting the in-memory volume after writing, confirming it is the correct function for this purpose.

[claude]

🟣 Pre-existing bugs in D-symmetry case of SetSymmetryLimits: (1) min() on line 162 is called with a single scalar argument, making it a no-op that never clamps psi_max, unlike the correct two-argument min(obj.psi_max, ...) pattern used for C-symmetry on line 153. (2) On line 165, obj.symmetry_symbol(2:end)*2 multiplies the ASCII codes of the character array by 2 before passing to EMC_str2double, so for D4 symmetry it computes EMC_str2double(104) = 104 instead of the correct value 8. These bugs predate this PR, which only reformatted indentation.

Extended reasoning...

Bug 1: min() with single scalar argument (line 162)

The C-symmetry case at line 153 correctly uses min(obj.psi_max, 360.0 / EMC_str2double(...)) to clamp psi_max to the symmetry-derived limit. The D-symmetry case at line 162 was intended to follow the same pattern but is missing the first argument:

% C-symmetry (correct):
obj.psi_max = min(obj.psi_max, 360.0 / EMC_str2double(obj.symmetry_symbol(2:end)));
% D-symmetry (bug):
obj.psi_max = min(360.0 / EMC_str2double(obj.symmetry_symbol(2:end)));

In MATLAB, min(scalar) returns the scalar unchanged — it is a no-op. This means psi_max is always overwritten with the symmetry-derived value rather than being clamped to the smaller of the two. If the user specifies a smaller search range than the symmetry allows, their preference is silently ignored.

Bug 2: ASCII multiplication instead of numeric multiplication (line 165)

The intent is to compute number_of_asymmetric_units as twice the symmetry order for dihedral groups (e.g., D4 has 8 asymmetric units). However, the *2 is applied inside the EMC_str2double call:

obj.number_of_asymmetric_units = EMC_str2double(obj.symmetry_symbol(2:end)*2);

In MATLAB, obj.symmetry_symbol(2:end) returns a character array. Multiplying a char array by 2 multiplies the ASCII code values, not the parsed number. For example, with D4 symmetry: obj.symmetry_symbol(2:end) is the char 4 (ASCII 52), so 4 * 2 = 104, and EMC_str2double(104) = 104. The correct result should be 8.

Step-by-step proof for D4:

1. obj.symmetry_symbol = D4
2. obj.symmetry_symbol(2:end) = 4 (char, ASCII value 52)
3. 4 * 2 = 52 * 2 = 104 (MATLAB promotes char to double via ASCII)
4. EMC_str2double(104) = 104 (already numeric)
5. obj.number_of_asymmetric_units = 104 (should be 8)

The fix should be: obj.number_of_asymmetric_units = EMC_str2double(obj.symmetry_symbol(2:end)) * 2;

Pre-existing nature

The % FIXME is this right? comment at line 158 already acknowledges uncertainty about the D-symmetry code. The PR diff shows only whitespace/indentation changes to these lines — the logic was not modified. These bugs predate this PR entirely.

[claude]

🟣 Pre-existing bug: sprintf('%aliStacks/%s_ali%d.fixed', ...) has a stray % before aliStacks, causing MATLAB to interpret %a as a hexadecimal float format specifier. This consumes tiltList{iTilt} as a numeric argument, producing a garbage file path that will fail at BH_multi_loadOrBin. The fix is to remove the leading % to match the correct else branch at line 408: aliStacks/%s_ali%d.fixed.

Extended reasoning...

Bug Description

At line 388 of ctf/BH_ctf_Correct3d.m, the sprintf format string begins with '%aliStacks/...' instead of 'aliStacks/...'. In MATLAB's sprintf, %a is a valid format specifier that formats a number as a hexadecimal floating-point value. This means the format string is parsed as: %a (hex float specifier) + liStacks/%s_ali%d.fixed (literal and remaining specifiers).

Code Path and Trigger Condition

This code path is reached when samplingRate > 1 (line 381). The buggy line is:

fullStack = sprintf('%aliStacks/%s_ali%d.fixed', tiltList{iTilt}, mapBackIter + 1);

The %a specifier consumes tiltList{iTilt} (a string/char array) as its argument. MATLAB will attempt to convert the string to a hex float representation, producing unexpected output like 0x0p+0liStacks/ or similar garbage. The %s specifier then has no remaining string argument and consumes mapBackIter + 1 (a number), and %d has no argument left at all.

Why Existing Code Doesn't Prevent This

The else branch at line 408 has the correct format string:

inputStack = sprintf('aliStacks/%s_ali%d.fixed', tiltList{iTilt}, mapBackIter + 1);

This shows the intended behavior — no leading %. The buggy line was likely introduced by an accidental keystroke that added % before aliStacks, and it has persisted because the samplingRate > 1 code path may not be exercised in routine testing.

Impact

When samplingRate > 1, the fullStack variable will contain a malformed path. This path is then passed to BH_multi_loadOrBin(fullStack, samplingRate, 2, false) at line 395, which will fail to locate the file, causing a runtime error during CTF correction. This effectively breaks the downsampled-stack processing path.

Step-by-Step Proof

1. Assume tiltList{iTilt} = 'TS_001' and mapBackIter = 3.
2. The format string '%aliStacks/%s_ali%d.fixed' is parsed by sprintf as: format specifier %a, literal liStacks/, format specifier %s, literal _ali, format specifier %d, literal .fixed.
3. %a consumes 'TS_001' — MATLAB converts char codes to hex floats, producing something like 0x1.18p+6 0x1.4cp+6 ...liStacks/.
4. %s then consumes mapBackIter + 1 = 4, converting it to a string '4'.
5. %d has no argument, causing an error or empty output.
6. Result: a garbage path instead of 'aliStacks/TS_001_ali4.fixed'.

Fix

Remove the stray % from the format string on line 388:

fullStack = sprintf('aliStacks/%s_ali%d.fixed', tiltList{iTilt}, mapBackIter + 1);