NaNIsaac docs

Appearance

Python API ​

Signatures and docstrings of the public functions and classes in qoraphase at docs/site (b77d3f1), read from the source without importing it. The guides on this site say how the pieces fit together; this page is for looking up an argument.

qoraphase.algorithms.er ​

Error Reduction.

ER ​

python
@dataclass(frozen=True)
class ER(Iterative):
    ...
``x <- P_S P_M x``.  Monotone in the Fourier error, and slow.

ER.step ​

python
def step(self, state, problem)

qoraphase.algorithms.hio ​

Hybrid Input-Output.

HIO ​

python
@dataclass(frozen=True)
class HIO(Iterative):
    beta: float = 0.9
``x <- P_M x`` inside the support, ``x - beta P_M x`` outside.

HIO.step ​

python
def step(self, state, problem)

qoraphase.algorithms.raar ​

Relaxed Averaged Alternating Reflections.

RAAR ​

python
@dataclass(frozen=True)
class RAAR(Iterative):
    beta: float = 0.9
``x <- (beta/2)(R_S R_M x + x) + (1 - beta) P_M x``

with ``R_M = 2 P_M - I`` and ``R_S = 2 P_S - I``.  ``beta -> 1`` is pure
reflection (fast, unstable), ``beta -> 0`` degenerates to modulus projection.

RAAR.step ​

python
def step(self, state, problem)

qoraphase.api ​

Public entry point.

default_recipe ​

python
def default_recipe() -> Recipe
RAAR, a detwin pass, then shrinkwrapped RAAR and ER.

The detwin runs before shrinkwrap starts updating the support.  At the
shrinkwrap settings shipped here the order makes little difference either
way; it matters a great deal at looser settings, where a support that is
still bloated gets cut in half and shrinkwrap -- which can only shrink --
locks the fragment in.  Doing the cut first is simply the order with fewer
ways to go wrong.

run ​

python
def run(problem: CDI, recipe=None, *, starts=1, batch=None, seed=None)
Execute *recipe* on *problem*.  See :func:`reconstruct`.

reconstruct ​

python
def reconstruct(intensity, *, starts: Literal[1]=1, **kwargs) -> Reconstruction
def reconstruct(intensity, *, starts: int, **kwargs) -> Ensemble
Array in, reconstruction out.

Args:
    intensity: measured intensity, center-origin, 2D or 3D.
    support: boolean, ``True`` inside the object.  ``None`` derives an
        autocorrelation support.
    mask: boolean, ``True`` marks *unusable* pixels.
    recipe: ``None`` uses :func:`default_recipe`.
    starts: number of random starts; >1 returns an :class:`Ensemble`.
    batch: starts run simultaneously under ``vmap``.  ``None`` runs them all
        at once, which is fastest but can exhaust GPU memory on large 3D
        data -- lower it if you hit OOM.
    seed: PRNG seed; ``None`` means 0, not "random".

qoraphase.core.ensemble ​

Multi-start results.

Ensemble ​

python
@dataclass(frozen=True)
class Ensemble:
    runs: tuple
    problem: object = None
    history: tuple = ()
A sequence of reconstructions that behaves like a Python collection.

Runs are *not* mutually aligned as they come out of the solver: each one is
the object up to its own global phase, translation and twin.
Compare them by error or by agreement (:attr:`consensus`), or call
:meth:`align` before doing anything that combines them.

``problem`` is carried so that a combined result can be scored against the
original measurement; it is absent for hand-built ensembles.  ``history``
holds the earlier generations of a :func:`~qoraphase.evolve` run, oldest first.

Ensemble.ranking ​

python
@property
def ranking(self) -> tuple
Run indices ordered best (lowest error) first.

Ensemble.best ​

python
@property
def best(self) -> Reconstruction
The run with the lowest Fourier error.

Runs stuck in different basins can fit the data equally well: on BCDI
benchmark data the phase error of single runs split into two groups
(about 0.08 and 0.2 rad) whose Fourier errors differed by 1e-3, so the
lowest error picked either.  :attr:`consensus` does not have that problem.

Ensemble.spread ​

python
@cached_property
def spread(self) -> tuple
Mean distance of each run from the others, after aligning them onto it.

``||align_to(other, run) - run|| / ||run||``, averaged over the other runs.
A run in the basin most starts reach is close to the rest.

Ensemble.consensus ​

python
@property
def consensus(self) -> Reconstruction
The run closest to all the others (smallest :attr:`spread`).

Needs no ground truth, only agreement between starts.  With fewer than
three runs there is no majority, and this is :attr:`best`.

Ensemble.align ​

python
def align(self, upsample: int=1, reference: Reconstruction | None=None) -> 'Ensemble'
Bring every run onto *reference* (default :attr:`best`), in phase, position and branch.

The reference supplies the convention, nothing more: if it happens to sit
in the conjugate-flipped branch, the whole aligned ensemble does too.
That is fine for combining runs, which is what this is for, and it is not
a statement about which branch is real.

``upsample`` is passed to :func:`~qoraphase.postprocess.align_to`; the
default aligns to whole pixels.

Ensemble.average ​

python
def average(self, upsample: int=1) -> Reconstruction
Coherent mean of the runs aligned onto :attr:`consensus`.

Aligning first is not optional: averaging raw runs cancels them against
each other, because they differ by exactly the symmetries the solver
cannot see.  The reference is the consensus run rather than the best one:
a lowest-error run from a minority basin would set the branch and
position, and the majority would then be averaged onto it.

qoraphase.core.problem ​

The CDI problem: fixed reconstruction inputs.

CDI ​

python
@dataclass(frozen=True, init=False)
class CDI:
    amplitude: jnp.ndarray
    mask: jnp.ndarray
    support: jnp.ndarray
    model: object
    def __init__(
        self,
        intensity=None,
        *,
        support=None,
        mask=None,
        model=None,
        support_threshold: float = 0.1,
        amplitude=None,
    )
Immutable reconstruction inputs.  Holds no iteration state.

Constructed from *intensity*; the stored field is the amplitude
``sqrt(clip(intensity, 0))`` so the square root is not recomputed inside the
iteration loop.

Registered as a pytree so it can cross ``jit`` / ``scan`` /
``vmap`` boundaries.  ``model`` is a meta field: it is static and hashable.

CDI.intensity ​

python
@property
def intensity(self) -> jnp.ndarray

CDI.shape ​

python
@property
def shape(self) -> tuple

CDI.reconstruct ​

python
def reconstruct(self, recipe=None, *, starts=1, batch=None, seed=None)

qoraphase.core.result ​

Single reconstruction result.

Reconstruction ​

python
@dataclass(frozen=True)
class Reconstruction:
    object: jnp.ndarray
    support: jnp.ndarray
    metrics: jnp.ndarray
    coherence: jnp.ndarray | None = None
One reconstructed object plus its error history.

The solution is only defined up to a global phase, a translation, and a
conjugate flip (twin image).  ``phase`` is therefore *not*
an absolute phase; remove the ramp before interpreting it physically.

Reconstruction.amplitude ​

python
@property
def amplitude(self) -> jnp.ndarray

Reconstruction.phase ​

python
@property
def phase(self) -> jnp.ndarray

Reconstruction.error ​

python
@property
def error(self) -> float
Final Fourier error; the ranking key.

Reconstruction.plot ​

python
def plot(self, axes=None)
Amplitude / phase panels.  3D shows the central slice of each axis.

qoraphase.core.state ​

Iteration state and initialisation.

ReconState ​

python
class ReconState(NamedTuple):
    obj: jnp.ndarray
    support: jnp.ndarray
    iteration: jnp.ndarray
    key: jax.Array
    coherence: jnp.ndarray | None = None
    lowpass: tuple | None = None
Everything that changes during iteration.

Metrics are deliberately absent: the carry must be fixed-size, and metrics
grow per iteration.  They come out of ``lax.scan`` as stacked outputs
instead.

ReconState.sigma ​

python
@property
def sigma(self) -> jnp.ndarray | None
Blur width of the data at the current iteration, or None.

random_init ​

python
def random_init(key: jax.Array, problem, support: jnp.ndarray) -> jnp.ndarray
Measured amplitude with uniformly random phase, back-transformed.

Masked pixels start at zero magnitude: what is stored under the mask is not a
measurement and must not steer the start.

initial_state ​

python
def initial_state(key: jax.Array, problem, obj: jnp.ndarray | None=None) -> ReconState
Start from *obj* (cut to the support) if given, else from random phases.

qoraphase.core.validate ​

Host-side validation.

Runs on concrete numpy arrays before anything is traced, so failures surface as
plain Python errors with useful messages instead of XLA shape errors.

check_inputs ​

python
def check_inputs(intensity, support, mask) -> None

check_run ​

python
def check_run(recipe, starts: int, batch: int | None) -> None

qoraphase.diagnostics ​

Pre-flight checks on a measurement.

Everything here takes arrays and returns numbers -- no file handling, no
reconstruction.  The point is to catch the mistakes that otherwise surface as a
plausible-looking but wrong result: a pattern that is not centred, an object
that is not oversampled enough to be unique, a background subtraction that went
negative.

None of these are fatal on their own.  They are the questions worth answering
before believing a reconstruction.

looks_inverted ​

python
def looks_inverted(intensity: np.ndarray, mask: np.ndarray) -> bool
True when the mask covers most pixels *and* most of the counts.

A real mask covers zeros (module gaps, padded frames) or a few saturated pixels; a mask with
its polarity inverted covers the pattern itself.

robust_peak_offset ​

python
def robust_peak_offset(intensity: np.ndarray) -> tuple
Like :func:`peak_offset`, after capping every pixel at its brightest face neighbour.

An isolated hot pixel or zinger then drops to the level around it, so it cannot decide where
the pattern sits; a real peak, whose neighbours are bright too, barely moves. Two adjacent hot
pixels still count as a peak.

inverted_mask_message ​

python
def inverted_mask_message(masked_fraction: float) -> str

flip_about_origin ​

python
def flip_about_origin(a: np.ndarray) -> np.ndarray
Reflect ``r -> -r`` for a center-origin array.

peak_offset ​

python
def peak_offset(intensity: np.ndarray) -> tuple
Offset of the brightest pixel from the array centre, in pixels.

friedel_offset ​

python
def friedel_offset(intensity: np.ndarray) -> tuple
Centring error inferred from Friedel symmetry, in pixels.

A correctly centred pattern satisfies ``I(q) == I(-q)``.  If instead the
pattern sits at ``c``, correlating it with its own reflection peaks at
``2c`` -- so this resolves a half-pixel error that the brightest pixel
cannot, and it works on broad patterns with no single dominant peak.

Returns half the correlation peak position, rounded, so an odd result means
the true offset is a half pixel.

autocorrelation_fill ​

python
def autocorrelation_fill(intensity: np.ndarray, threshold: float=0.1) -> tuple
Autocorrelation extent per axis, as a fraction of the array.

This is the quantity actually measured, and the honest one to read: a
well-oversampled measurement puts the autocorrelation well inside the array,
and a value approaching 1 means it no longer fits.

An estimate, not a measurement -- it inherits the arbitrariness of
*threshold*, and a noisy autocorrelation floor inflates the extent.

oversampling ​

python
def oversampling(intensity: np.ndarray, threshold: float=0.1) -> tuple
Linear oversampling ratio per axis, from :func:`autocorrelation_fill`.

The autocorrelation spans twice the object, so ``ratio = 2 / fill``.

Read this together with the fill, not on its own.  The autocorrelation
cannot extend past the array, so ``fill <= 1`` and this ratio is **always at
least 2** -- including on data that is badly undersampled, where the
autocorrelation has simply been clipped.  A ratio of 2.0 is not a pass; it
means the estimate ran out of room.

Diagnosis ​

python
@dataclass(frozen=True)
class Diagnosis:
    shape: tuple
    dtype: str
    peak_offset: tuple
    friedel_offset: tuple
    autocorrelation_fill: tuple
    oversampling: tuple
    total_oversampling: float
    negative_fraction: float
    masked_fraction: float
    support_fraction: float
    dynamic_range: float
    warnings: tuple = field(default_factory=tuple)
What the checks found.  ``print()`` it.

diagnose ​

python
def diagnose(intensity, mask=None, support=None, threshold: float=0.1) -> Diagnosis
Check a measurement against the assumptions QoraPhase makes about it.

qoraphase.genetic ​

Genetic algorithm over generations of reconstructions.

Several random starts are run, the best one is kept as the reference, every other
run is aligned to it and bred with it, and the next generation restarts from the
children instead of from random phases.  Breeding formulas and the idea follow
cohere's ``ga_utils`` / ``dvc_utils.breed``.

breed ​

python
def breed(mode: str, alpha: jnp.ndarray, beta: jnp.ndarray) -> jnp.ndarray
Child of *alpha* (the best run) and *beta*, which must already be aligned.

Alignment (position, twin branch, global phase) is what
:meth:`Ensemble.align` does; breeding unaligned runs mixes unrelated phases.
Both are first rotated by the same global phase so that the mean phase of
*alpha* is zero, which keeps the phase averages away from the +-pi wrap.

``_pa`` modes keep the phase of *alpha* and mix only the amplitudes.

evolve ​

python
def evolve(
    problem: CDI,
    recipe=None,
    *,
    starts: int = 8,
    generations: int = 3,
    breed: str = 'sqrt_ab',
    cull: int | Sequence[int] = 0,
    batch: int | None = None,
    seed: int | None = None,
    upsample: int = 10,
) -> Ensemble
Run *recipe* for several generations, breeding with the best run in between.

Generation 0 is ``starts`` random starts.  After each generation the
``cull[g]`` worst runs (by final error) are dropped, the rest are aligned to
the best (:meth:`Ensemble.align`, to ``1 / upsample`` pixel) and bred with it,
and the next generation runs *recipe* again from the children.  The best run
is carried over unchanged.  The returned :class:`Ensemble` is the last
generation; earlier ones are in its ``history``.

Args:
    recipe: one recipe used for every generation, or a sequence with one per
        generation; ``None`` uses :func:`default_recipe`.
    generations: total generations, including the first.
    breed: a key of :data:`BREEDERS`.
    cull: an int (dropped after every generation) or a sequence indexed by
        generation; missing entries are 0.  The population must stay at one
        run or more, and ``batch`` must divide every generation's size.
    batch, seed: as in :func:`~qoraphase.reconstruct`.  The same seed gives
        the same result.

Children restart on the problem's support, not on the shrinkwrapped one of
the previous generation, and a partial-coherence kernel starts afresh.

qoraphase.metrics.errors ​

Error metrics.

fourier_error ​

python
def fourier_error(
    computed_magnitude: jnp.ndarray,
    amplitude: jnp.ndarray,
    mask: jnp.ndarray,
) -> jnp.ndarray
Normalised Fourier-space error over unmasked pixels only.

    E_F = sqrt( sum_valid (|F(obj)| - amplitude)^2 / sum_valid amplitude^2 )

Takes the already-computed magnitude rather than the object, so that
ModulusProjection can hand over the FFT it has already paid for
(computing this separately would double the FFT cost).

qoraphase.models.fft ​

Forward models: object space <-> detector space.

FFT ​

python
@dataclass(frozen=True)
class FFT:
    ...
Unitary centered FFT.

All public arrays are center-origin: index ``n // 2`` is the
origin.  The ``ifftshift`` / ``fftshift`` sandwich is entirely internal, so
users never call ``fftshift`` themselves.

``adjoint(forward(x)) == x`` because ``norm="ortho"`` makes the transform
unitary.

FFT.forward ​

python
def forward(self, obj: jnp.ndarray) -> jnp.ndarray

FFT.adjoint ​

python
def adjoint(self, wave: jnp.ndarray) -> jnp.ndarray

qoraphase.operators.coherence ​

Partial coherence: the measured intensity is the coherent one blurred by a kernel.

With a partially coherent beam the detector records ``|F(obj)|^2 * K``, the
coherent intensity convolved with a coherence kernel ``K`` (Clark et al.,
Nat. Commun. 3, 993 (2012)).  Fitting a coherent model to such data fills the
missing contrast with spurious structure in the object.  Here ``K`` is estimated
during the run by Richardson-Lucy deconvolution, and the modulus constraint
compares the measurement with the blurred model instead.

convolve ​

python
def convolve(intensity: jnp.ndarray, kernel: jnp.ndarray) -> jnp.ndarray
Circular convolution of center-origin arrays; *kernel* may be smaller.

lucy_update ​

python
def lucy_update(kernel, coherent, measured, mask, iterations: int)
Richardson-Lucy steps for ``measured = coherent * kernel``.

The usual normalising sum is a constant here and drops out when the kernel
is rescaled to unit sum after each step.  Masked pixels do not vote.

blurred_magnitude ​

python
def blurred_magnitude(wave: jnp.ndarray, coherence) -> jnp.ndarray
``sqrt(|wave|^2 * K)``, or ``|wave|`` when there is no kernel.

UpdateCoherence ​

python
@dataclass(frozen=True)
class UpdateCoherence(Stage):
    iterations: int = 20
Re-estimate the coherence kernel from the current object.

The model intensity is scaled to the measured total first, so the kernel
carries only the shape of the blur.  Produces no error metric.

UpdateCoherence.run ​

python
def run(self, state, problem)

PartialCoherence ​

python
@dataclass(frozen=True)
class PartialCoherence(Stage):
    body: Recipe
    size: int | tuple
Run *body* with a coherence kernel, starting from a flat one.

The kernel stays in the state afterwards, so later stages keep using it and
the result reports it as ``Reconstruction.coherence``.

PartialCoherence.run ​

python
def run(self, state, problem)

partial_coherence ​

python
def partial_coherence(*stages, kernel=7, lucy: int=20) -> PartialCoherence
``qp.partial_coherence(qp.repeat(8, qp.RAAR(50), qp.Shrinkwrap()), kernel=7)``

Every iterative algorithm inside *stages* is preceded by a kernel update
(*lucy* Richardson-Lucy steps), so the update interval is the length of the
algorithm stages.  *kernel* is the kernel size in pixels, an odd int or one
per axis; it must not exceed the array.  Start it after a few hundred
coherent iterations: the first estimate comes from the current object.

qoraphase.operators.lowpass ​

Low-pass resolution progression: fit blurred data first, the data itself last.

Early iterations see a Gaussian-smoothed version of the measured amplitude, and
the blur narrows to the raw data over the run.  This is cohere's
``lowpass_filter`` (``features.LowPassFilter``): the filter is applied to the
*data*, not to the computed wave, and the sigma ramps linearly with the
iteration number.

Why the data and not the wave: the modulus constraint is a projection onto
"amplitude equals the target".  Smoothing the target gives a different, smoother
constraint set that every algorithm projects onto unchanged, and the iterate is
never filtered, so nothing but the target differs between filtered and
unfiltered iterations.  Smoothing the wave instead would put a lossy operator in
the middle of the iteration, and it would not be a projection.

A Gaussian blur of the amplitude in the detector is a Gaussian envelope on the
object in real space.  Speckle contrast drops and the large-scale structure of
the pattern survives, so the early iterations are asked to explain less.  The
blur is in detector pixels, so a good sigma depends on the fringe spacing,
roughly the array size divided by the object size.

filtered_amplitude ​

python
def filtered_amplitude(problem, sigma) -> jnp.ndarray
The measured amplitude blurred to *sigma* pixels, ignoring masked pixels.

Normalised convolution: masked pixels (beamstop, hot pixels) are left out of
the average instead of dragging their stored value into their neighbours.
Where no unmasked pixel is in reach the measured value is kept.

LowPass ​

python
@dataclass(frozen=True)
class LowPass(Stage):
    body: Recipe
    sigma: tuple
Run *body* with the data blurred by a sigma that ramps from ``sigma[0]`` to
``sigma[1]`` over the iterations of *body*, then go back to the data.

The schedule lives in the iteration state, so the body can be any mix of
algorithms, shrinkwrap, detwin, phase limits and partial coherence.  Error
metrics inside the body are against the blurred data.

LowPass.run ​

python
def run(self, state, problem)

low_pass ​

python
def low_pass(*stages, sigma=(2.0, 0.5)) -> LowPass
``qp.low_pass(qp.RAAR(100), qp.Shrinkwrap(), sigma=(1.5, 0.0))``

*sigma* is the blur width in detector pixels at the first and at the last
iteration of *stages*, interpolated linearly in between; a single number is a
constant blur.  A last sigma of 0 ends on the unblurred data.  Metrics from
these iterations are against the blurred data.  To use it in one generation
of :func:`~qoraphase.evolve` only, give that generation its own recipe.

qoraphase.operators.modulus ​

Modulus (Fourier magnitude) projection.

modulus_project ​

python
def modulus_project(obj: jnp.ndarray, problem, coherence=None, sigma=None)
Project onto the measured-magnitude set.

Returns ``(projected_obj, fourier_error)``.  The error comes for free from
the forward transform this projection already computed; a separate metric
operator would double the FFT cost.

Masked pixels keep their *computed* magnitude -- they are unconstrained, not
zero.  With a *coherence* kernel the measurement is compared with the
blurred model intensity, and the wave is rescaled by the same ratio.  With a
*sigma* the measured amplitude is first low-passed to that width (see
:func:`~qoraphase.low_pass`), and the error is against the filtered data.

ModulusProjection ​

python
@dataclass(frozen=True)
class ModulusProjection:
    ...
Public handle for :func:`modulus_project`.

qoraphase.operators.phase ​

Phase constraint: confine the object's phase to a range during a run.

PhaseLimited ​

python
@dataclass(frozen=True)
class PhaseLimited(Stage):
    algorithm: Iterative
    low: float
    high: float
Run *algorithm* with its support narrowed, at every iteration, to the
pixels whose current phase lies in ``[low, high]`` (radians).

The narrowed support is used for that iteration only; the stored support
is left alone, so shrinkwrap and later stages see the full one.  Phase is
measured from zero, so a range around zero also pins the global phase.
This is cohere's ``phc`` and the "phase constraint" of the legacy MATLAB
code; a weakly strained crystal starts close to real, and keeping it there
early on discourages solutions with wild phase.

PhaseLimited.run ​

python
def run(self, state, problem)

phase_limit ​

python
def phase_limit(low: float, high: float, *stages) -> Recipe
``qp.phase_limit(-0.5, 0.5, qp.repeat(10, qp.HIO(50), qp.Shrinkwrap()))``

Every iterative algorithm inside *stages*, however deeply nested, runs
under :class:`PhaseLimited`; everything else runs unchanged.

qoraphase.operators.shrinkwrap ​

Shrinkwrap support update.

gaussian_blur ​

python
def gaussian_blur(array: jnp.ndarray, sigma: float) -> jnp.ndarray
Separable Gaussian blur applied in Fourier space.

Works for any ndim.  ``sigma`` is in pixels; the transfer function is
``exp(-2 (pi sigma f)^2)`` with ``f`` in cycles/pixel, which is 1 at f=0 and
therefore preserves the sum.

Shrinkwrap ​

python
@dataclass(frozen=True)
class Shrinkwrap(Stage):
    threshold: float = 0.3
    sigma: float = 1.0
    method: str = 'max'
Re-derive the support from the current object magnitude.

The blurred amplitude is thresholded at ``threshold`` times a reference
level, chosen by *method*:

``"max"``
    the largest blurred amplitude.  The classic Marchesini formulation, and
    the reference drifts down as the reconstruction sharpens.
``"rms"``
    the root-mean-square blurred amplitude inside the current support.  A
    bright outlier lifts this too -- it is inside the support after all --
    but diluted across every pixel rather than setting the reference on its
    own, so the support update is harder to destabilise.

The two thresholds are on different scales and do not transfer: PyNX uses
0.28 for ``"rms"``, and the classic Marchesini ``"max"`` value is 0.2.

The defaults here are 0.3 / 1.0, not 0.2 / 2.0.  On a loose autocorrelation
support -- the case that actually needs shrinkwrap -- 0.2 / 2.0 never
converged in testing while 0.3 / 1.0 did, across every object tried.  That
measurement is on synthetic data, so treat it as a starting point and not a
law; these two numbers are the first thing to re-tune on a real dataset.

Produces no error metric, so it contributes nothing to the metric history.

Shrinkwrap.run ​

python
def run(self, state, problem)

qoraphase.operators.support ​

Support projection and autocorrelation support estimate.

support_project ​

python
def support_project(obj: jnp.ndarray, support: jnp.ndarray) -> jnp.ndarray
Zero everything outside the support.

SupportProjection ​

python
@dataclass(frozen=True)
class SupportProjection:
    ...
Public handle for :func:`support_project`.

v0.1 applies no real/positivity constraint -- the object is complex.  A
``constraint=`` argument belongs here when that changes.

half_support ​

python
def half_support(support: jnp.ndarray, axis: int=0) -> jnp.ndarray
Keep the half of *support* below its centre of mass along *axis*.

Used for detwinning: a support cut in half is no longer symmetric under
``r -> -r``, so it fits only one of the two conjugate solutions.  Cutting at
the centre of mass rather than the array centre keeps the cut on the object
even when the object sits off-centre.

autocorrelation_support ​

python
def autocorrelation_support(intensity: jnp.ndarray, threshold: float=0.1) -> jnp.ndarray
Default support when the user supplies none.

The autocorrelation is about twice the object size, so this is a loose
starting support and relies on shrinkwrap to tighten it.  A measured support
converges better; pass one when you have it.

qoraphase.postprocess ​

Resolving the CDI ambiguities after a reconstruction.

A reconstruction is determined only up to a global phase, a translation, and a
conjugate flip.  Nothing during iteration can fix that -- all three are exact
symmetries of the constraints -- so it is dealt with here, once, on concrete
arrays.  These run on the host and are not part of any traced loop.

A linear phase across the object means its Fourier transform is centred off the
array centre.  :func:`remove_ramp` measures that offset from the data and cancels
it; it does not know the diffraction geometry, so it removes whatever ramp is
there, including one that carries physical information (strain, a displaced
Bragg peak).

flip_about_origin ​

python
def flip_about_origin(x: jnp.ndarray) -> jnp.ndarray
Reflect ``r -> -r`` for a center-origin array.

Reversing the axes is not enough: with the origin at index ``n // 2`` the
reflection maps ``i`` to ``(n - i) % n``, i.e. the reversal rolled by one.

twin ​

python
def twin(x: jnp.ndarray) -> jnp.ndarray
The conjugate-flipped solution, indistinguishable from ``x`` in the data.

register ​

python
def register(obj, reference, upsample: int=1)
Shift and peak value that best match *obj* onto *reference*.

Returns ``(shift, value)``: ``obj`` is *reference* displaced by ``+shift``
(in pixels), and ``value`` is the complex correlation at that shift, whose
angle is the global phase between the two.

The whole-pixel peak comes from the FFT cross-correlation.  With
``upsample > 1`` it is then refined in decades around the peak, evaluating
the correlation only on a small grid by matrix-multiply DFT (Guizar-Sicairos
et al., Opt. Lett. 33, 156, 2008), down to a step of ``1 / upsample`` pixel
or finer.

align_to ​

python
def align_to(obj: jnp.ndarray, reference: jnp.ndarray, upsample: int=1) -> jnp.ndarray
Transform *obj* into the equivalent solution closest to *reference*.

The translation and global phase come from the cross-correlation peak, and
the branch is chosen by taking whichever of ``obj`` and ``twin(obj)``
correlates more strongly.

This is a *relative* operation, and the twin part of it especially so: it
adopts whatever branch *reference* is in.  It does not, and cannot, tell you
which branch is physically correct -- both fit the measurement exactly.  What
it does give you is a common representation, which is all that averaging an
ensemble needs.

With the default ``upsample=1`` the translation is to whole pixels.  A larger
value locates the peak to about ``1 / upsample`` pixel (see :func:`register`)
and applies the remainder as a Fourier shift, which is exact only for an
object that does not reach the array edge.

remove_ramp ​

python
def remove_ramp(
    obj,
    support=None,
    *,
    pad: int = 2,
    iterations: int = 3,
    zero_phase: bool = False,
)
Remove the linear phase ramp by centring the object's Fourier peak.

The ramp is read off as the centre of mass of ``|F|^2`` of the zero-padded
object (``pad`` times larger per axis, which sets how finely the offset is
resolved), and cancelled with a phase ramp that is zero at the array centre.
The offset is a real number, so this is a sub-pixel shift of the spectrum.
The centre of mass is unbiased when the ramp-free spectrum is centrosymmetric
(a real object); for a general complex object it is an estimate.

A spectrum displaced from the centre is cut off asymmetrically by the array
edge, which pulls the centre of mass towards zero.  Each pass removes most of
the ramp and therefore most of that bias, so ``iterations`` passes are made.

*support*, if given, restricts the measurement to the object (noise outside
it would otherwise pull the peak); the returned array is not masked.  With
``zero_phase`` the global phase is then set so that the mean phase inside the
support, weighted by amplitude, is zero.  Without a support the voxels above
a tenth of the peak amplitude stand in for it.

Traceable: usable under ``jit`` with ``pad``, ``iterations`` and
``zero_phase`` static.

qoraphase.prtf ​

Phase retrieval transfer function.

Independent reconstructions agree on the Fourier components the data pins down
and disagree elsewhere.  Averaging their transforms, after removing the
ambiguities, keeps the first and washes out the second, so the ratio of the
averaged amplitude to the measured one falls from 1 towards 0 as the runs stop
agreeing (Chapman et al., J. Opt. Soc. Am. A 23, 1179, 2006).

PRTF ​

python
@dataclass(frozen=True)
class PRTF:
    frequency: jnp.ndarray
    value: jnp.ndarray
    crossing: float
Radially averaged PRTF.

``frequency`` holds the shell centres in cycles per object pixel and
``value`` the PRTF in each shell (nan for a shell with no usable pixel).
``crossing`` is the frequency at which ``value`` first falls below 1/e,
linearly interpolated, or nan if it never does; ``resolution`` is its
reciprocal, a full period in object pixels.

PRTF.resolution ​

python
@property
def resolution(self) -> float

prtf_curve ​

python
def prtf_curve(objects, intensity, mask=None)
PRTF of already-aligned *objects* as ``(frequency, value, crossing)``.

``value`` is the ratio, per shell, of the summed amplitude of the mean
transform to the summed measured amplitude, over unmasked pixels.  Summing
before dividing keeps near-empty pixels from dominating a shell.  Pure
``jnp``: traceable under ``jit``.

prtf ​

python
def prtf(objects, intensity, mask=None, *, align: bool=True, upsample: int=100) -> PRTF
PRTF of an :class:`~qoraphase.Ensemble` or a sequence of objects.

With ``align=True`` every run is first brought onto the reference (the best
run of an ensemble, the first of a sequence) in phase, position and twin
branch, to ``1 / upsample`` pixel.  Sub-pixel alignment matters here: a
residual shift of a fraction of a pixel is a phase error that grows with
frequency and would read as lost resolution.  Pass ``align=False`` for
objects that are already aligned.

*intensity* and *mask* are the measurement (``mask=True`` marks an unusable
pixel, which is left out).  The curve is meaningful only for two or more
runs; a single run has a PRTF of 1 by construction.

qoraphase.recipes.recipe ​

Stage protocol and declarative Recipe.

Code order is execution order.  No operator-overloading DSL.

Stage ​

python
class Stage:
    ...
Anything a Recipe can run.

``run(state, problem) -> (new_state, metrics)`` where ``metrics`` is a 1-D
float array with one entry per iteration that produced one.  Stages that
produce no error metric (Shrinkwrap) return an empty array.

Stage.run ​

python
def run(self, state, problem)

Iterative ​

python
@dataclass(frozen=True)
class Iterative(Stage):
    n: int
A stage that repeats ``step`` ``n`` times under ``lax.scan``.

``n`` is the scan length and therefore static: changing it recompiles.

Iterative.step ​

python
def step(self, state, problem)

Iterative.run ​

python
def run(self, state, problem)

Recipe ​

python
@dataclass(frozen=True, init=False)
class Recipe(Stage):
    stages: tuple
    def __init__(self, *stages, **kwargs)
Ordered sequence of stages.

Recipe.run ​

python
def run(self, state, problem)

Repeat ​

python
@dataclass(frozen=True)
class Repeat(Stage):
    n: int
    body: Recipe
Run ``body`` ``n`` times.

Repeat.run ​

python
def run(self, state, problem)

Detwin ​

python
@dataclass(frozen=True)
class Detwin(Stage):
    body: Recipe
    axis: int = 0
Run ``body`` against a support cut in half, then restore the support.

The twin ``conj(obj(-r))`` fits the measurement exactly, so nothing in the
data prefers one branch; a centrosymmetric support does not separate them
either and the iteration stalls in a mixture.  Halving the support removes
that symmetry for a few cycles and biases the solution into one branch, after
which the full support is restored.

Which branch it lands in follows from the cut, not from the data.

Detwin.run ​

python
def run(self, state, problem)

repeat ​

python
def repeat(n: int, *stages) -> Repeat
``qp.repeat(10, qp.RAAR(50), qp.Shrinkwrap())``

detwin ​

python
def detwin(*stages, axis: int=0) -> Detwin
``qp.detwin(qp.HIO(20))`` -- break the twin symmetry, then carry on.

Each project section is built from that project’s docs/ directory.

Isaac Yong