Skip to content

python-api.md

Latest commit

 

History

History
1155 lines (896 loc) · 32.5 KB

python-api.md

File metadata and controls

1155 lines (896 loc) · 32.5 KB
Raw

Python API guide

This guide documents the Python interface of structure-factor 1.0.0. It is intended both as an API reference and as an explanation of the numerical choices behind the returned arrays.

The library analyzes particle coordinates in a periodic, orthorhombic, two-dimensional box. It provides:

  • a direct microscopic particle structure factor;
  • FFT spectra of occupancy and number-density grids;
  • an FFT spectrum of a thresholded binary phase field;
  • FFT spectra of fields smoothed by Gaussian, square top-hat, or disk kernels;
  • circular radial averages in physical wave-number coordinates;
  • publication-style Matplotlib figures.

The implementation uses float64 real arrays and complex128 Fourier work arrays. It does not perform file output, plotting, or display as a side effect of the numerical calculation.

Imports and API organization

The package root exports the five symbols needed by most programs:

from structure_factor import (
    StructureFactorConfig,
    StructureFactorResult,
    compute_structure_factor,
    load_positions,
    occupancy_grid,
)

These names are also available from their defining modules:

from structure_factor.core import (
    StructureFactorConfig,
    StructureFactorResult,
    compute_structure_factor,
    occupancy_grid,
)
from structure_factor.io import load_positions

Two additional helpers are intentionally imported from their modules:

from structure_factor.core import radial_average
from structure_factor.plotting import create_figures

radial_average is useful when rebinning a computed spectrum or averaging an independently produced two-dimensional spectrum. create_figures is kept in the plotting module so numerical users can import the package root without importing Matplotlib.

The installed version is available as:

import structure_factor

print(structure_factor.__version__)

Names beginning with an underscore are implementation details and should not be used as application APIs.

Minimal end-to-end calculation

The following example loads a coordinate file, computes the microscopic particle structure factor, and inspects the radial result:

import numpy as np

from structure_factor import (
    StructureFactorConfig,
    compute_structure_factor,
    load_positions,
)

positions = load_positions("positions.dat")

config = StructureFactorConfig(
    box_size=(128.0, 128.0),
    method="particle",
    max_mode=16,
)
result = compute_structure_factor(positions, config)

print(f"N = {result.particle_count}")
print(f"2D spectrum shape = {result.spectrum.shape}")

assert result.radial_wave_numbers is not None
assert result.radial_spectrum is not None
peak = int(np.argmax(result.radial_spectrum))
print(f"strongest radial bin: k = {result.radial_wave_numbers[peak]:.6g}")

The assert statements are useful to type checkers because radial attributes are optional. They are present by default, but are None if compute_radial_average=False was requested.

Coordinate, indexing, and unit conventions

Understanding these conventions prevents the most common interpretation errors.

Periodic box and coordinates

Coordinates are supplied as an array with shape (N, 2):

positions[:, 0]  -> x coordinates
positions[:, 1]  -> y coordinates

For box_size=(Lx, Ly), every point must already satisfy:

$$ 0 \leq x < L_x, \qquad 0 \leq y < L_y. $$

The upper boundary is excluded because it is periodically identical to zero. The library validates these bounds but deliberately does not wrap coordinates. If an upstream simulation can emit unwrapped positions, wrap them explicitly:

positions = positions % np.asarray((Lx, Ly))

Do this only when periodic wrapping is scientifically appropriate for the input data.

X-first two-dimensional arrays

Every numerical two-dimensional array is stored with x first:

value = result.spectrum[ix, iy]
k_vector = (result.kx[ix], result.ky[iy])

Likewise, a real-space field uses:

value = result.source_field[ix, iy]

This is different from Matplotlib's image convention, in which the first array index is the vertical row. The package plotting functions transpose an array exactly once before displaying it. If making a custom image, normally use field.T or spectrum.T.

Units

The library does not prescribe a particular length unit. It only requires consistent inputs.

Quantity Unit
positions, box_size, grid_size caller's length unit
particle_diameter, kernel_length caller's length unit
kx, ky, radial_wave_numbers inverse length
occupancy source field particles per cell
density source field particles per length squared
kernel source field smoothed particle mass per cell
binary source field dimensionless, exactly -1 or +1
computed spectra dimensionless
radial mode counts number of Fourier modes

Wave numbers are angular wave numbers: reciprocal-lattice spacings contain 2*pi, not just 1/L.

Connected spectra and the zero mode

Every estimator returns a connected spectrum. Grid methods subtract the spatial mean before the FFT, and the direct method discards the particle amplitude at $\mathbf{k}=0$. The DC spectrum element is then assigned exactly zero.

The unconnected microscopic definition would instead have $S(\mathbf{0})=N$. Neither the forced connected value zero nor the unconnected fixed-$N$ value $N$ is, by itself, an estimate of compressibility from one fixed-particle-number snapshot.

Loading coordinate files

load_positions

load_positions(path: str | pathlib.Path) -> numpy.ndarray

load_positions reads a whitespace-delimited numeric text file with numpy.loadtxt. Each nonempty row must contain exactly two values:

0.25 1.50
2.75 3.10
4.00 0.20

The returned array:

  • has shape (N, 2), including the single-particle case;
  • has float64 dtype;
  • is C-contiguous;
  • contains at least one row;
  • contains only finite values.

The loader does not know the box dimensions, so it does not check periodic bounds. Bound checking occurs in occupancy_grid or compute_structure_factor.

Example:

from pathlib import Path

from structure_factor import load_positions

positions = load_positions(Path("examples") / "xy.dat")
print(positions.shape)

Errors:

  • OSError is raised when the file cannot be opened or read;
  • ValueError is raised for an empty file, nonnumeric content, a column count other than two, or NaN/infinite coordinates.

Configuring an estimator

StructureFactorConfig

StructureFactorConfig(
    box_size,
    method="particle",
    grid_size=2.0,
    density_threshold=0.7,
    particle_diameter=1.0,
    max_mode=8,
    kernel="gaussian",
    kernel_length=None,
    direct_chunk_size=4096,
)

The configuration is a frozen dataclass. Construction validates and normalizes every field, including fields unused by the selected method. This gives one consistent configuration contract and catches invalid values early.

Field Default Meaning
box_size required Positive periodic lengths (Lx, Ly)
method "particle" "particle", "occupancy", "density", "binary", or "kernel"
grid_size 2.0 Requested maximum grid-cell spacing
density_threshold 0.7 Binary local packing-fraction threshold
particle_diameter 1.0 Diameter used to calculate local packing fraction
max_mode 8 Direct reciprocal cutoff, scalar or (mx, my)
kernel "gaussian" "gaussian", "top-hat", or "disk"
kernel_length None Gaussian standard deviation, square half-width, or disk radius
direct_chunk_size 4096 Number of particles in each direct-sum work chunk

All lengths must be finite and positive. density_threshold may be zero but not negative. Mode cutoffs and direct_chunk_size must be positive integers; booleans are rejected rather than being interpreted as zero or one.

If kernel_length=None, the constructor resolves it to the current grid_size and stores that positive float:

config = StructureFactorConfig((10.0, 10.0), grid_size=1.25)
assert config.kernel_length == 1.25

This is a construction-time default. The configuration is immutable after it has been created.

Derived grid properties

grid_size is a target maximum spacing, not necessarily the exact cell width. For direction $\alpha$:

$$ M_\alpha = \left\lceil \frac{L_\alpha}{\Delta} \right\rceil, \qquad h_\alpha = \frac{L_\alpha}{M_\alpha}, $$

where $\Delta$ is grid_size.

The configuration exposes both derived quantities:

config = StructureFactorConfig(
    box_size=(5.0, 7.0),
    method="occupancy",
    grid_size=2.0,
)

assert config.grid_shape == (3, 4)
assert config.effective_cell_size == (5.0 / 3.0, 7.0 / 4.0)
  • config.grid_shape returns (nx, ny);
  • config.effective_cell_size returns (hx, hy) in length units;
  • config.mode_cutoffs returns an explicit (mx, my) pair.

For example:

isotropic = StructureFactorConfig((8.0, 8.0), max_mode=12)
anisotropic = StructureFactorConfig((8.0, 12.0), max_mode=(12, 18))

assert isotropic.mode_cutoffs == (12, 12)
assert anisotropic.mode_cutoffs == (12, 18)

The grid properties are available even for the particle method, although that method does not use a real-space grid.

Constructing the occupancy grid directly

occupancy_grid

occupancy_grid(
    positions,
    config: StructureFactorConfig,
) -> numpy.ndarray

This helper assigns each particle to a regular cell:

$$ i_x = \left\lfloor \frac{x}{h_x} \right\rfloor, \qquad i_y = \left\lfloor \frac{y}{h_y} \right\rfloor. $$

It returns an x-first float64 array with shape config.grid_shape. Values are numerically integer counts and their sum equals N.

import numpy as np

from structure_factor import StructureFactorConfig, occupancy_grid

positions = np.array(
    [
        [0.1, 0.1],
        [0.8, 0.3],
        [2.9, 1.9],
    ]
)
config = StructureFactorConfig(
    box_size=(3.0, 2.0),
    method="occupancy",
    grid_size=1.0,
)

counts = occupancy_grid(positions, config)
assert counts.shape == (3, 2)
assert counts.sum() == positions.shape[0]

The configured method does not change this helper; only box and grid geometry matter. Assignment is cell membership without interpolation. Cloud-in-cell, triangular-shaped-cloud, and assignment-window corrections are not implemented in this release.

occupancy_grid accepts a single coordinate pair with shape (2,) and promotes it to (1, 2). It rejects empty, malformed, nonfinite, and out-of-box coordinates.

Computing a structure factor

compute_structure_factor

compute_structure_factor(
    positions,
    config: StructureFactorConfig,
    *,
    compute_radial_average: bool = True,
) -> StructureFactorResult

positions may be any numeric array-like object with shape (N, 2). A single pair with shape (2,) is accepted. The function validates the coordinates, dispatches to the configured estimator, and optionally computes a radial average.

Every method produces a two-dimensional spectrum with matching kx and ky axes. Grid methods additionally retain both the raw occupancy and the real-space source_field that was Fourier transformed.

Estimator definitions and normalization

Method Fourier amplitude Power normalization
particle $\sum_j e^{-i\mathbf{k}\cdot\mathbf{r}_j}$ at nonzero modes $N$
occupancy FFT of $n_{ij}-\langle n\rangle$ $N$
density cell area times FFT of $\rho_{ij}-\langle\rho\rangle$ $N$
binary FFT of $\psi_{ij}-\langle\psi\rangle$ $M_xM_y$
kernel FFT of centered, smoothed cell mass $N$

The different normalizations are part of the scientific definitions. In particular, the binary spectrum is an order-parameter spectrum and its absolute amplitude should not be compared directly to a particle-normalized spectrum.

Direct particle method

config = StructureFactorConfig(
    box_size=(64.0, 48.0),
    method="particle",
    max_mode=(16, 12),
    direct_chunk_size=4096,
)
result = compute_structure_factor(positions, config)

The sampled reciprocal vectors are:

$$ \mathbf{k}_{mn}

\left( \frac{2\pi m}{L_x}, \frac{2\pi n}{L_y} \right), $$

with $-m_x \leq m \leq m_x$ and $-m_y \leq n \leq m_y$. At nonzero modes:

$$ S(\mathbf{k})

\frac{1}{N} \left| \sum_{j=1}^{N} e^{-i\mathbf{k}\cdot\mathbf{r}_j} \right|^2. $$

The spectrum has shape (2*mx + 1, 2*my + 1). It is computed directly from the coordinates without grid quantization. direct_chunk_size changes temporary phase-matrix memory, but not modes, normalization, or results apart from possible floating-point summation-order differences.

Occupancy method

config = StructureFactorConfig(
    box_size=(64.0, 48.0),
    method="occupancy",
    grid_size=1.0,
)
result = compute_structure_factor(positions, config)

Particles are binned into counts $n_{ij}$. The method transforms the centered counts and divides squared amplitudes by N:

$$ S_n(\mathbf{k})

\frac{ \left|\operatorname{FFT}(n-\langle n\rangle)\right|^2 }{N}. $$

The returned source_field and occupancy both contain cell counts.

Number-density method

config = StructureFactorConfig(
    box_size=(64.0, 48.0),
    method="density",
    grid_size=1.0,
)
result = compute_structure_factor(positions, config)

The real-space source field is:

$$ \rho_{ij} = \frac{n_{ij}}{h_xh_y}. $$

The density Fourier amplitude is multiplied by the cell area before its power is divided by N. Therefore:

$$ (h_xh_y)\widehat{\delta\rho}

\widehat{\delta n}, $$

and density and occupancy return the same spectrum up to floating-point round-off. They retain different source_field arrays: one is number density and the other is particles per cell.

occupancy_result = compute_structure_factor(
    positions,
    StructureFactorConfig((64.0, 48.0), method="occupancy", grid_size=1.0),
)
density_result = compute_structure_factor(
    positions,
    StructureFactorConfig((64.0, 48.0), method="density", grid_size=1.0),
)

np.testing.assert_allclose(
    density_result.spectrum,
    occupancy_result.spectrum,
    rtol=1e-12,
    atol=1e-12,
)

Binary phase-field method

config = StructureFactorConfig(
    box_size=(64.0, 48.0),
    method="binary",
    grid_size=2.0,
    particle_diameter=1.0,
    density_threshold=0.7,
)
result = compute_structure_factor(positions, config)

For particle diameter $\sigma$, the local packing fraction is:

$$ \phi_{ij}

\frac{n_{ij}\pi\sigma^2}{4h_xh_y}. $$

It is mapped to:

$$ \psi_{ij}

\begin{cases} -1, & \phi_{ij}<\phi_c,\ +1, & \phi_{ij}\geq\phi_c. \end{cases} $$

A cell exactly equal to the threshold is dense (+1). The connected phase field is transformed and squared amplitudes are divided by the total number of cells. This estimator is useful for domain morphology and coarsening. It is not the microscopic particle structure factor.

Kernel-smoothed method

config = StructureFactorConfig(
    box_size=(64.0, 48.0),
    method="kernel",
    grid_size=1.0,
    kernel="gaussian",
    kernel_length=2.0,
)
result = compute_structure_factor(positions, config)

The method samples a nonnegative kernel on the periodic grid, normalizes its discrete weights to sum to one, and circularly convolves it with occupancy. Consequently, the smoothed source field conserves particle mass:

np.testing.assert_allclose(
    result.source_field.sum(),
    positions.shape[0],
)

kernel_length has a shape-dependent meaning:

  • Gaussian: standard deviation $\ell$ in $K(r)\propto\exp[-r^2/(2\ell^2)]$;
  • square top-hat: half-width $\ell$, with support $|x|\leq\ell$ and $|y|\leq\ell$;
  • disk: radius $\ell$, with support $r\leq\ell$.

The output is the filtered occupancy spectrum:

$$ S_K(\mathbf{k})

|W_{\mathrm{grid}}(\mathbf{k})|^2 S_{\mathrm{occupancy}}(\mathbf{k}). $$

No kernel deconvolution is applied. In particular, correction would be unstable around transform zeros of top-hat and disk windows and near the Nyquist limit.

Computing several methods

Each calculation receives a complete configuration. This makes method comparisons explicit and reproducible:

from structure_factor import (
    StructureFactorConfig,
    compute_structure_factor,
    load_positions,
)

box = (1284.0, 1284.0)
grid_size = 2.0
positions = load_positions("examples/xy.dat")

configs = {
    "particle": StructureFactorConfig(
        box,
        method="particle",
        max_mode=8,
    ),
    "occupancy": StructureFactorConfig(
        box,
        method="occupancy",
        grid_size=grid_size,
    ),
    "density": StructureFactorConfig(
        box,
        method="density",
        grid_size=grid_size,
    ),
    "binary": StructureFactorConfig(
        box,
        method="binary",
        grid_size=grid_size,
        density_threshold=0.7,
        particle_diameter=1.0,
    ),
    "gaussian": StructureFactorConfig(
        box,
        method="kernel",
        grid_size=grid_size,
        kernel="gaussian",
        kernel_length=2.0,
    ),
    "top-hat": StructureFactorConfig(
        box,
        method="kernel",
        grid_size=grid_size,
        kernel="top-hat",
        kernel_length=2.0,
    ),
    "disk": StructureFactorConfig(
        box,
        method="kernel",
        grid_size=grid_size,
        kernel="disk",
        kernel_length=2.0,
    ),
}

results = {
    name: compute_structure_factor(positions, config)
    for name, config in configs.items()
}

The complete executable comparison is available in examples/xy_all_methods.ipynb.

Understanding StructureFactorResult

compute_structure_factor returns a frozen StructureFactorResult dataclass. Freezing prevents attribute rebinding, but it does not make the contained NumPy arrays read-only. Treat result arrays as calculation outputs; copy an array before intentionally modifying it.

Its complete constructor is:

StructureFactorResult(
    config,
    particle_count,
    spectrum,
    kx,
    ky,
    radial_wave_numbers=None,
    radial_spectrum=None,
    radial_mode_counts=None,
    source_field=None,
    occupancy=None,
)

Applications normally receive this object from compute_structure_factor rather than constructing it directly. The dataclass constructor does not revalidate array shapes, axis alignment, or method-dependent fields.

Result fields

Attribute Shape Meaning and units
config scalar object Validated configuration used by the calculation
particle_count integer Number of input particles, N
spectrum (len(kx), len(ky)) Connected dimensionless 2D spectrum
kx (len(kx),) First-dimension angular wave numbers, inverse length
ky (len(ky),) Second-dimension angular wave numbers, inverse length
radial_wave_numbers (B,) or None Mean actual radius in each radial bin, inverse length
radial_spectrum (B,) or None Mean spectrum in each radial bin
radial_mode_counts (B,) or None Number of 2D modes in each radial bin
source_field grid shape or None Real-space field transformed by a grid estimator
occupancy grid shape or None Raw particles-per-cell grid

Method-dependent shapes are:

Method Spectrum shape source_field occupancy
particle (2*mx + 1, 2*my + 1) None None
any grid method config.grid_shape present present

For a grid method, shifted FFT axes are constructed as:

kx = 2 * np.pi * np.fft.fftshift(np.fft.fftfreq(nx, d=hx))
ky = 2 * np.pi * np.fft.fftshift(np.fft.fftfreq(ny, d=hy))

Thus spectrum[i, j] is always aligned with (kx[i], ky[j]). On an even FFT grid there is one negative Nyquist mode and no separate positive partner. The direct particle grid is instead explicitly symmetric because it contains integer modes from -max_mode to +max_mode, inclusive.

The DC location can be found robustly from the axes:

dc = (
    int(np.argmin(np.abs(result.kx))),
    int(np.argmin(np.abs(result.ky))),
)
assert result.spectrum[dc] == 0.0

Optional radial arrays

Use the convenience property when code accepts results computed with or without radial averaging:

if result.has_radial_average:
    assert result.radial_wave_numbers is not None
    assert result.radial_spectrum is not None
    assert result.radial_mode_counts is not None
    print(result.radial_wave_numbers)

The property tests that all three attributes are non-None. A valid radial calculation can still return three empty arrays if the sampled spectrum contains no retained nonzero modes.

To avoid radial work:

result = compute_structure_factor(
    positions,
    config,
    compute_radial_average=False,
)

assert result.radial_wave_numbers is None
assert result.radial_spectrum is None
assert result.radial_mode_counts is None

Radial averaging

Automatic radial average

By default, compute_structure_factor calls radial_average with:

$$ \Delta k_\mathrm{bin}

\min\left(\frac{2\pi}{L_x},\frac{2\pi}{L_y}\right). $$

For every Fourier sample:

$$ k = \sqrt{k_x^2+k_y^2}. $$

The algorithm:

  1. retains the largest origin-centered circle sampled symmetrically by both axes;
  2. thereby excludes unmatched even-grid Nyquist rows or columns and reciprocal-space corners with anisotropic angular sampling;
  3. excludes the zero mode;
  4. assigns each retained mode to the nearest integer-width bin using floor(k / bin_width + 0.5);
  5. omits bin zero and all empty bins;
  6. returns the mean actual wave number, mean power, and contributing mode count for every retained bin.

Using the actual mean radius matters for rectangular boxes, where a bin can contain modes at several nearby but unequal radii.

radial_average

radial_average(
    spectrum,
    kx,
    ky,
    *,
    bin_width,
) -> tuple[radial_wave_numbers, radial_spectrum, radial_mode_counts]

The function accepts any finite x-first two-dimensional spectrum whose shape is (len(kx), len(ky)), not only a StructureFactorResult.

To rebin a result more coarsely:

from structure_factor.core import radial_average

fundamental = min(
    2.0 * np.pi / result.config.box_size[0],
    2.0 * np.pi / result.config.box_size[1],
)

radial_k, radial_s, counts = radial_average(
    result.spectrum,
    result.kx,
    result.ky,
    bin_width=2.0 * fundamental,
)

bin_width must be finite and positive. A very large custom width can place small positive radii into bin zero; because bin zero is omitted, those modes will not appear in the result. If no valid modes remain, the function returns three empty one-dimensional arrays.

No characteristic length, integral moment, uncertainty estimate, or error bar is imposed by the library. Those choices depend on the physical application. For a visibly isolated positive peak, one possible diagnostic is:

positive = (
    (radial_k > 0.0)
    & np.isfinite(radial_s)
    & (radial_s > 0.0)
)
if np.any(positive):
    shown_k = radial_k[positive]
    shown_s = radial_s[positive]
    peak_k = shown_k[np.argmax(shown_s)]
    wavelength = 2.0 * np.pi / peak_k

This reports a wavelength associated with the largest sampled radial bin; it does not establish the physical origin or statistical significance of that peak.

Plotting and saving figures

create_figures

from structure_factor.plotting import create_figures

figures = create_figures(result)

The return value is a mapping from semantic names to newly created Matplotlib Figure objects:

Key When present Content
"source_field" result.source_field is not None Real-space grid field
"spectrum" always Two-dimensional reciprocal-space spectrum
"radial_spectrum" radial wave numbers and powers are not None Log-log radial curve

No window is opened and no file is written by create_figures. Figure ownership passes to the caller.

The current plotting behavior is:

  • binary source fields use two categorical colors for dilute -1 and dense +1 cells;
  • occupancy, density, and kernel fields use a continuous color scale;
  • 2D spectra use logarithmic normalization, mask nonpositive values and the connected DC cell, and can add an adaptive low-$k$ inset;
  • the spectrum color floor uses the 0.5th percentile of positive, non-DC power so isolated numerical speckles do not flatten the visible range;
  • radial figures use logarithmic axes and annotate the largest displayed positive bin with both $k_\mathrm{peak}$ and $2\pi/k_\mathrm{peak}$;
  • plotting style is applied in a local Matplotlib context and does not change the application's global style settings.

Because Matplotlib uses row-first images, the plotting functions transpose x-first arrays internally. Custom plotting code should follow the same rule.

Saving and closing figures

from pathlib import Path

import matplotlib.pyplot as plt

from structure_factor.plotting import create_figures

output_dir = Path("results") / "figures"
output_dir.mkdir(parents=True, exist_ok=True)

figures = create_figures(result)
try:
    for name, figure in figures.items():
        figure.savefig(
            output_dir / f"{name}.png",
            dpi=200,
            bbox_inches="tight",
        )
finally:
    for figure in figures.values():
        plt.close(figure)

Closing figures is important in loops, notebooks that recompute many results, and long-running services.

For a noninteractive machine, select a headless backend before importing pyplot or structure_factor.plotting:

import matplotlib

matplotlib.use("Agg")

import matplotlib.pyplot as plt

from structure_factor.plotting import create_figures

Alternatively, set the backend in the shell:

MPLBACKEND=Agg python3 analysis.py

The full-data plotting workflow is implemented in examples/xy_plotting.py.

Saving numerical arrays

Version 1.0.0 does not expose a separate Python result-serialization function. NumPy can save the arrays and essential metadata directly:

payload = {
    "spectrum": result.spectrum,
    "kx": result.kx,
    "ky": result.ky,
    "particle_count": np.asarray(result.particle_count),
    "box_size": np.asarray(result.config.box_size),
    "method": np.asarray(result.config.method),
}

if result.radial_wave_numbers is not None:
    payload["radial_wave_numbers"] = result.radial_wave_numbers
if result.radial_spectrum is not None:
    payload["radial_spectrum"] = result.radial_spectrum
if result.radial_mode_counts is not None:
    payload["radial_mode_counts"] = result.radial_mode_counts
if result.source_field is not None:
    payload["source_field"] = result.source_field
if result.occupancy is not None:
    payload["occupancy"] = result.occupancy

np.savez_compressed("structure_factor_result.npz", **payload)

For reproducibility, also save all configuration fields used by the analysis, especially grid size, threshold, particle diameter, mode cutoffs, and kernel parameters.

Warnings and errors

Configuration errors

StructureFactorConfig rejects:

  • a box_size that does not contain exactly two finite positive lengths;
  • a method outside particle, occupancy, density, binary, and kernel;
  • nonfinite or nonpositive grid_size, particle_diameter, or kernel_length;
  • a nonfinite or negative density_threshold;
  • a max_mode that is not a positive integer or a two-positive-integer tuple;
  • an unknown kernel name;
  • a non-positive or non-integral direct_chunk_size;
  • lengths so large that squaring them would overflow.

Every field is validated even if the selected method does not use it.

Position errors

occupancy_grid and compute_structure_factor raise:

  • TypeError when config is not a StructureFactorConfig;
  • ValueError for an empty position set, an incompatible shape, nonfinite coordinates, or a coordinate outside the half-open box.

Nonnumeric objects can also fail during conversion to float64.

Radial-average errors

radial_average raises TypeError when an input cannot be converted to a numeric NumPy array. It raises ValueError when:

  • bin_width is not finite and positive;
  • kx or ky is not one-dimensional;
  • spectrum.shape != (len(kx), len(ky));
  • spectrum or axes contain NaN or infinite values.

Plotting errors

For normal results returned by compute_structure_factor, plotting inputs are already consistent. If a StructureFactorResult is constructed or modified manually, create_figures can raise ValueError for:

  • a source field that is not a finite, nonempty two-dimensional array;
  • a spectrum shape not aligned with one-dimensional kx and ky;
  • nonfinite wave-number axes;
  • radial wave-number and power arrays with incompatible dimensions.

Nonfinite or nonpositive spectrum values are masked in the logarithmic display. Nonfinite or nonpositive radial values are omitted from the log-log curve.

Runtime warnings

Grid estimators emit RuntimeWarning in two scientifically uninformative or under-resolved cases:

  1. If config.grid_shape == (1, 1), subtracting the only cell's mean leaves an identically zero connected spectrum.
  2. For the kernel method, if kernel_length is smaller than the larger effective cell spacing, the sampled kernel may be under-resolved and can collapse to the occupancy estimator.

Warnings are not errors: a result is still returned. They indicate that the chosen discretization should be reconsidered.

Performance and convergence guidance

Direct particle method

For cutoffs (mx, my), the computational cost scales as:

$$ O\left( N(2m_x+1)(2m_y+1) \right). $$

The default max_mode=8 is deliberately a low-wave-number calculation. Doubling an isotropic mode cutoff approximately quadruples the number of sampled mode pairs.

Particles are processed in chunks. Temporary phase matrices scale roughly with:

$$ \texttt{direct_chunk_size} \left[(2m_x+1)+(2m_y+1)\right]. $$

Decrease direct_chunk_size when temporary memory is constrained. Increase it only after measuring performance on the target machine; it does not change reciprocal resolution.

Grid and kernel methods

Occupancy construction is linear in particle count plus grid storage. The FFT work scales approximately as:

$$ O(M_xM_y\log(M_xM_y)). $$

Grid memory scales with $M_xM_y$. In two dimensions, halving grid_size in both directions produces roughly four times as many cells when box lengths are large compared with the spacing.

The kernel method performs additional FFT convolution work and constructs a full periodic kernel array. compute_structure_factor also retains both occupancy and source_field, so account for both arrays in memory planning.

Radial averaging constructs a full array of wave-number radii comparable in shape to the spectrum. Set compute_radial_average=False when only the anisotropic 2D spectrum is needed.

Scientific convergence

Defaults make the package immediately runnable; they are not universal convergence parameters.

  • Increase max_mode until the direct spectrum covers the scientifically relevant wave-number range.
  • Refine grid_size and confirm that conclusions are stable. Grid assignment aliases high-wave-number power, and version 1.0.0 applies no assignment-window correction.
  • Vary kernel_length and ensure the kernel is resolved by several grid cells when a smooth continuum-like field is intended.
  • Treat Gaussian, top-hat, and disk kernels as different filters, not interchangeable numerical implementations.
  • Use density and occupancy as an internal normalization cross-check: their spectra should agree at the same grid resolution.
  • Do not compare the absolute binary phase-spectrum amplitude directly with particle-normalized estimators.
  • Inspect the full 2D spectrum before relying only on a radial average; radial averaging intentionally removes directional information.

For the million-particle examples/xy.dat dataset, start with the provided defaults and examples before increasing direct mode cutoffs or refining the grid:

For the estimator definitions and first-release scope, see Scientific methods. The command-line guide documents batch execution and the CLI archive schema, while Examples collects the complete xy.dat workflows.