pyfracval.utils#

Utility functions for vector operations and array manipulation.

Deprecated since version This: module is being split into domain-specific sub-modules. Import from the specific module instead:

All symbols remain importable from this module for backward compatibility.

Module Contents#

pyfracval.utils.batch_check_overlaps_cca(coords1, radii1, coords2_batch, radii2, tolerance)[source]#

Check overlap for batch of cluster2 configurations (CCA).

Uses Numba parallel loops to evaluate multiple cluster configurations simultaneously.

Parameters:

  • coords1 (np.ndarray) – Cluster 1 coordinates (n1, 3)

  • radii1 (np.ndarray) – Cluster 1 radii (n1,)

  • coords2_batch (np.ndarray) – Batch of cluster 2 configurations (n_batch, n2, 3)

  • radii2 (np.ndarray) – Cluster 2 radii (n2,) - same for all configurations

  • tolerance (float) – Overlap tolerance

Returns:

(n_batch,) array of max overlap values for each configuration

Return type:

np.ndarray

pyfracval.utils.batch_rotate_cluster_cca(coords2_in, cm2, cand2_idx, vec_0, i_vec, j_vec, angles)[source]#

Batch rotate cluster2 for multiple angles (CCA).

For each angle, calculates the target position on the intersection circle, then rotates the entire cluster to align the candidate particle with that target.

Parameters:

  • coords2_in (np.ndarray) – Cluster 2 coordinates (n2, 3)

  • cm2 (np.ndarray) – Center of mass of cluster 2 (3,)

  • cand2_idx (int) – Index of candidate particle in cluster 2

  • vec_0 (np.ndarray) – [x0, y0, z0, r0] - center and radius of intersection circle

  • i_vec (np.ndarray) – First basis vector (3,)

  • j_vec (np.ndarray) – Second basis vector (3,)

  • angles (np.ndarray) – Array of rotation angles (n_angles,)

Returns:

(n_angles, n2, 3) array of rotated cluster configurations

Return type:

np.ndarray

pyfracval.utils.calculate_cluster_properties(coords, radii, df, kf, densities=None)[source]#

Calculate aggregate properties: total mass, Rg, center of mass, Rmax.

Parameters:

  • coords (np.ndarray) – Nx3 array of particle coordinates.

  • radii (np.ndarray) – N array of particle radii.

  • df (float) – Fractal dimension used for Rg calculation.

  • kf (float) – Fractal prefactor used for Rg calculation.

  • densities (np.ndarray, optional) – Per-particle densities; None means uniform. Affects the total mass and the center of mass (and hence r_max), which in turn feed the Gamma equation.

Returns:

A tuple containing:

  • total_mass (float): Sum of individual particle masses.

  • rg (float | None): Radius of gyration calculated via calculate_rg, or None if calculation failed.

  • cm (np.ndarray | None): 3D center of mass coordinates, or None if calculation failed.

  • r_max (float): Maximum distance from the center of mass to any particle center in the aggregate.

Returns (0.0, 0.0, np.zeros(3), 0.0) for empty inputs (N=0).

Return type:

tuple[float, float | None, np.ndarray | None, float]

pyfracval.utils.calculate_mass(radii, densities=None)[source]#

Calculate particle mass from radii and (optionally) per-particle density.

m_i = (4/3) pi r_i^3 rho_i.

Parameters:

  • radii (np.ndarray) – Array of particle radii.

  • densities (np.ndarray, optional) – Per-particle densities. None (the default) means uniform density, in which case mass is proportional to r^3 and every density-aware quantity reduces exactly to its single-material form. Supplying densities is what makes heterogeneous aggregates - different materials, not just different sizes - physically meaningful, since the center of mass, radius of gyration and the Gamma equation are all mass-weighted.

Returns:

Array of corresponding particle masses.

Return type:

np.ndarray

pyfracval.utils.calculate_rg(radii, npp, df, kf)[source]#

Calculate the radius of gyration using the fractal scaling law.

Implements the formula Rg = a * (N / kf)^(1/Df), where ‘a’ is the geometric mean radius calculated from the input radii array. See [Morán et al., 2019] and morphology context [Filippov et al., 2000].

Parameters:

  • radii (np.ndarray) – Array of radii of particles in the cluster/aggregate.

  • npp (int) – Number of primary particles (N) in the cluster.

  • df (float) – Fractal dimension (Df).

  • kf (float) – Fractal prefactor (kf).

Returns:

The calculated radius of gyration (Rg). Returns 0.0 if npp is 0, kf or df is zero, or if calculation fails (e.g., log error).

Return type:

float

pyfracval.utils.compute_empirical_rg(coords, radii)[source]#

Compute empirical Rg directly from particle coordinates (mass-weighted).

Unlike calculate_rg which uses the fractal scaling law Rg = a*(N/kf)^(1/Df), this function measures Rg from the actual spatial distribution of particles.

Treats each particle as a point mass. See compute_empirical_rg_polydisperse() for the form that also carries each particle’s own gyration radius (paper Eq. 4) - required wherever the result has to be consistent with the Gamma equation.

Parameters:

  • coords (np.ndarray) – Nx3 array of particle center coordinates.

  • radii (np.ndarray) – N array of particle radii. Mass is proportional to r^3.

Returns:

Empirical (mass-weighted) radius of gyration.

Return type:

float

pyfracval.utils.compute_pair_correlation_dimensions(coords, radii, n_bins=50)[source]#

Estimate fractal dimension from pair-correlation (mass-radius) scaling.

Computes the cumulative mass M(r) as a function of radial distance from the centre of mass. For a fractal aggregate, M(r) ~ r^Df, so a log-log fit gives the empirical Df.

The fit uses raw (un-normalised) cumulative mass because normalised mass fractions are in [0,1] whose logs are non-positive, breaking the log-linear regression.

Parameters:

  • coords (np.ndarray) – Nx3 array of particle center coordinates.

  • radii (np.ndarray) – N array of particle radii.

  • n_bins (int) – Number of radial bins (default 50).

Returns:

r_bins : np.ndarray — bin edge radii (n_bins+1,) r_centers : np.ndarray — bin centre radii (n_bins,) M_r : np.ndarray — cumulative normalised mass fraction within each radius (n_bins+1,) empirical_Df : float — slope of log(M) vs log(r) fit fit_r_squared : float — R^2 of the linear fit empirical_kf : float — estimated kf from the fit

Return type:

dict with keys

pyfracval.utils.gamma_calculation(m1, rg1, radii1, m2, rg2, radii2, df, kf, use_mass=False, all_radii=None)[source]#

Calculates Gamma_pc for adding the next monomer (aggregate 2).

The Gamma_pc relation follows the FracVAL CCA/PCA formulation [Morán et al., 2019].

Parameters:

  • use_mass (bool, default False) – Which form of the Gamma equation to solve. False substitutes particle counts for the masses, giving Filippov et al. (2000) Eq. 7 - what the Fortran PCA does (PCA_cca.f90’s Gamma_calculation takes n1, n2, n3) and what this port has historically done everywhere. True uses the true masses passed in as m1/m2, giving Moran et al. (2019) Eq. 6 - the paper’s central polydisperse contribution, and what the Fortran CCA actually does (CCA_module.f90:301). Identical for monodisperse primary particles; they diverge as polydispersity grows. See NOTE.md 1.2.

  • all_radii (np.ndarray, optional) – If provided, the geometric mean radius for rg3 is computed from this full set of radii (matching Fortran behaviour where R contains all N particles). When None the geometric mean is taken from the local combined set (radii1 + radii2).

  • rg3_override (float, optional) – Use this radius of gyration for the combined aggregate instead of deriving it from the scaling law. Only meaningful together with measured (rather than scaling-law) rg1/rg2 - see cca/pairing.py’s measured-Rg feedback.

pyfracval.utils.validate_fractal_structure(coords, radii, target_df, target_kf, rg_rtol=0.05)[source]#

Validate that generated aggregate matches target fractal parameters.

Compares theoretical Rg (from scaling law) vs empirical Rg (from coordinates), and estimates the actual fractal dimension from mass-radius scaling.

Parameters:

  • coords (np.ndarray) – Nx3 array of particle center coordinates.

  • radii (np.ndarray) – N array of particle radii.

  • target_df (float) – Target fractal dimension.

  • target_kf (float) – Target fractal prefactor.

  • rg_rtol (float) – Relative tolerance for Rg agreement (default 5%).

Returns:

N : int — number of particles theoretical_rg : float — Rg from scaling law Rg = a*(N/kf)^(1/Df) empirical_rg : float — Rg measured from coordinates rg_error_pct : float — (empirical - theoretical)/theoretical * 100 rg_ok : bool — rg_error_pct < rg_rtol * 100 empirical_Df : float — Df estimated from mass-radius scaling target_Df : float — target fractal dimension df_error : float — empirical_Df - target_Df fit_r_squared : float — goodness of fit for Df estimation empirical_kf : float — estimated fractal prefactor target_kf : float — target fractal prefactor

Return type:

dict with keys

pyfracval.utils.random_point_sc(case, spheres_1_ext, spheres_2_ext, rng=None)[source]#

Sample a random point on the appropriate spherical cap for CCA’s ext_case=1 contact-point search (Fortran Random_point_SC, see docs/FracVAL/CCA_module.f90).

Used when the “shell” spheres (defined by [Dmin, Dmax] distance from each cluster’s center of mass) overlap in a way that a single Dmax/Dmax intersection circle (the ext_case=0 path used by two_sphere_intersection()) doesn’t capture - see the case 1/2/3 branch in cca/sticking.py::_cca_sticking_v1.

Parameters:

  • case (int) – Which shell-overlap case applies (1, 2, or 3 - see caller).

  • spheres_1_ext (np.ndarray) – [x, y, z, d_min, d_max] for each cluster’s shell.

  • spheres_2_ext (np.ndarray) – [x, y, z, d_min, d_max] for each cluster’s shell.

  • rng (np.random.Generator, optional) – Random generator to use; a fresh default one if not given.

Returns:

(x, y, z, valid) - a random point on the selected sphere’s surface, within the appropriate cap angle. valid is False for an unrecognized case or a degenerate (coincident) center pair.

Return type:

tuple[float, float, float, bool]

pyfracval.utils.rodrigues_rotation(vectors, axis, angle)[source]#

Rotate vector(s) around an axis using Rodrigues’ rotation formula.

Parameters:

  • vectors (np.ndarray) – A single 3D vector or an Nx3 array of vectors to rotate.

  • axis (np.ndarray) – The 3D rotation axis (does not need to be normalized).

  • angle (float) – The rotation angle in radians.

Returns:

The rotated vector or Nx3 array of rotated vectors. Returns the original vectors if the axis norm is near zero.

Return type:

np.ndarray

Raises:

ValueError – If input vectors is not 1D (3,) or 2D (N, 3).

pyfracval.utils.spherical_cap_angle(sphere_1, sphere_2)[source]#

Critical polar angle of the spherical cap of sphere_1 that lies within sphere_2 (Fortran Spherical_cap_angle, see docs/FracVAL/CCA_module.f90).

Parameters:

  • sphere_1 (np.ndarray) – [x, y, z, r] for each sphere.

  • sphere_2 (np.ndarray) – [x, y, z, r] for each sphere.

Returns:

The cap half-angle in radians, measured from the axis connecting the two sphere centers.

Return type:

float

pyfracval.utils.two_sphere_intersection(sphere_1, sphere_2, rng=None)[source]#

Find the intersection circle of two spheres and pick a random point.

Calculates the center (x0, y0, z0) and radius (r0) of the intersection circle, defines basis vectors (i_vec, j_vec) for the circle’s plane, and returns a random point (x, y, z) on that circle based on a random angle (theta).

Handles edge cases: spheres too far, one contained, coincidence, touching.

Parameters:

  • sphere_1 (np.ndarray) – Definition of the first sphere: [x1, y1, z1, r1].

  • sphere_2 (np.ndarray) – Definition of the second sphere: [x2, y2, z2, r2].

Returns:

A tuple containing:

  • x, y, z (float): Coordinates of a random point on the intersection.

  • theta (float): Random angle (radians) used to generate the point.

  • vec_0 (np.ndarray): [x0, y0, z0, r0] - center and radius of the intersection circle (r0=0 if spheres touch at a point).

  • i_vec (np.ndarray): First basis vector of the intersection plane.

  • j_vec (np.ndarray): Second basis vector of the intersection plane.

  • valid (bool): True if a valid intersection (circle or point) exists, False otherwise (e.g., separate, contained, coincident).

Return type:

tuple[float, float, float, float, np.ndarray, np.ndarray, np.ndarray, bool]

Note

https://mathworld.wolfram.com/Sphere-SphereIntersection.html

pyfracval.utils.calculate_max_overlap_cca(coords1, radii1, coords2, radii2)[source]#

Calculate max overlap between two particle clusters (Numba optimized).

Overlap is defined as 1 - distance / (radius1 + radius2) for overlapping pairs, max(0).

Parameters:

  • coords1 (np.ndarray) – Nx3 coordinates of cluster 1.

  • radii1 (np.ndarray) – N radii of cluster 1.

  • coords2 (np.ndarray) – Mx3 coordinates of cluster 2.

  • radii2 (np.ndarray) – M radii of cluster 2.

Returns:

Maximum overlap fraction found between any particle in cluster 1 and any particle in cluster 2. Returns 0.0 if no overlap.

Return type:

float

pyfracval.utils.calculate_max_overlap_cca_auto(coords1, radii1, coords2, radii2, tolerance=1e-06)[source]#

Check max overlap between two clusters during CCA sticking.

FIX (PyFracVAL-xwx): Always use the sequential early-termination path. The previous parallel dispatch was counterproductive for CCA: in sticking, clusters are placed touching (high overlap probability), so early termination fires almost immediately. The parallel path computes ALL n1*n2 pairs even when the first pair already overlaps, making it 78x slower for large clusters.

Benchmark (n1=n2=256, 65536 pairs): parallel=109µs, fast=1.4µs.

Parameters:

  • coords1 (np.ndarray) – Cluster 1 coordinates (n1, 3)

  • radii1 (np.ndarray) – Cluster 1 radii (n1,)

  • coords2 (np.ndarray) – Cluster 2 coordinates (n2, 3)

  • radii2 (np.ndarray) – Cluster 2 radii (n2,)

  • tolerance (float, optional) – Overlap tolerance for early termination (default: 1e-6)

Returns:

Maximum overlap fraction

Return type:

float

pyfracval.utils.calculate_max_overlap_cca_fast(coords1, radii1, coords2, radii2, tolerance=1e-06)[source]#

Calculate max overlap between clusters with early termination (optimized).

This optimized version includes: 1. Early termination: Returns immediately when overlap exceeds tolerance 2. Bounding sphere pre-check: Avoids sqrt for particles far apart 3. Sequential execution: Trades parallelization for early exit

Overlap is defined as 1 - distance / (radius1 + radius2).

Performance: ~2-3x faster than parallel version when overlap is found early.

Parameters:

  • coords1 (np.ndarray) – Nx3 coordinates of cluster 1.

  • radii1 (np.ndarray) – N radii of cluster 1.

  • coords2 (np.ndarray) – Mx3 coordinates of cluster 2.

  • radii2 (np.ndarray) – M radii of cluster 2.

  • tolerance (float, optional) – Overlap tolerance threshold for early termination (default: 1e-6).

Returns:

Maximum overlap fraction found. Returns immediately if overlap > tolerance.

Return type:

float

pyfracval.utils.calculate_max_overlap_cca_parallel(coords1, radii1, coords2, radii2)[source]#

Calculate max overlap for CCA with parallel execution (no early termination).

This version uses Numba prange to parallelize overlap checks. Computes all pair overlaps in parallel.

Use for n1 * n2 > PARALLEL_OVERLAP_THRESHOLD (~200 pairs).

Parameters:

  • coords1 (np.ndarray) – Cluster 1 coordinates (n1, 3)

  • radii1 (np.ndarray) – Cluster 1 radii (n1,)

  • coords2 (np.ndarray) – Cluster 2 coordinates (n2, 3)

  • radii2 (np.ndarray) – Cluster 2 radii (n2,)

Returns:

Maximum overlap fraction found across all particle pairs

Return type:

float

pyfracval.utils.calculate_max_overlap_pca(coords_agg, radii_agg, coord_new, radius_new)[source]#

Calculate max overlap between a new particle and an aggregate (Numba).

Overlap is defined as 1 - distance / (radius_new + radius_agg) for overlapping pairs, max(0).

Parameters:

  • coords_agg (np.ndarray) – Nx3 coordinates of the existing aggregate.

  • radii_agg (np.ndarray) – N radii of the aggregate particles.

  • coord_new (np.ndarray) – 3D coordinates of the new particle.

  • radius_new (float) – Radius of the new particle.

Returns:

Maximum overlap fraction found between the new particle and any particle in the aggregate. Returns 0.0 if no overlap.

Return type:

float

pyfracval.utils.calculate_max_overlap_pca_auto(coords_agg, radii_agg, coord_new, radius_new, tolerance=1e-06)[source]#

Auto-dispatch to parallel or sequential overlap check based on size.

For large aggregates (n > PARALLEL_OVERLAP_THRESHOLD), uses parallel version without early termination. For small aggregates, uses sequential with early exit.

Parameters:

  • coords_agg (np.ndarray) – Current aggregate coordinates (n_agg, 3)

  • radii_agg (np.ndarray) – Current aggregate radii (n_agg,)

  • coord_new (np.ndarray) – New particle coordinates (3,)

  • radius_new (float) – New particle radius

  • tolerance (float, optional) – Overlap tolerance for early termination (default: 1e-6)

Returns:

Maximum overlap fraction

Return type:

float

pyfracval.utils.calculate_max_overlap_pca_fast(coords_agg, radii_agg, coord_new, radius_new, tolerance=1e-06)[source]#

Calculate max overlap with early termination (optimized for speed).

This optimized version includes: 1. Early termination: Returns immediately when overlap exceeds tolerance 2. Bounding sphere pre-check: Avoids sqrt for particles far apart 3. Sequential execution: Trades parallelization for early exit

Overlap is defined as 1 - distance / (radius_new + radius_agg).

Performance: ~2-3x faster than parallel version when overlap is found early.

Parameters:

  • coords_agg (np.ndarray) – Nx3 coordinates of the existing aggregate.

  • radii_agg (np.ndarray) – N radii of the aggregate particles.

  • coord_new (np.ndarray) – 3D coordinates of the new particle.

  • radius_new (float) – Radius of the new particle.

  • tolerance (float, optional) – Overlap tolerance threshold for early termination (default: 1e-6).

Returns:

Maximum overlap fraction found. Returns immediately if overlap > tolerance.

Return type:

float

pyfracval.utils.calculate_max_overlap_pca_parallel(coords_agg, radii_agg, coord_new, radius_new)[source]#

Calculate max overlap for PCA with parallel execution (no early termination).

This version uses Numba prange to parallelize overlap checks across all aggregate particles. Trade-off: No early termination, but faster for large N.

Use for n_agg > PARALLEL_OVERLAP_THRESHOLD (~200 particles).

Parameters:

  • coords_agg (np.ndarray) – Current aggregate coordinates (n_agg, 3)

  • radii_agg (np.ndarray) – Current aggregate radii (n_agg,)

  • coord_new (np.ndarray) – New particle coordinates (3,)

  • radius_new (float) – New particle radius

Returns:

Maximum overlap fraction found across all particles

Return type:

float

pyfracval.utils.batch_calculate_positions_pca(vec_0, i_vec, j_vec, angles)[source]#

Calculate batch of positions on intersection circle for PCA.

Uses Numba parallel loops to compute multiple rotation positions simultaneously.

Parameters:

  • vec_0 (np.ndarray) – [x0, y0, z0, r0] - center and radius of intersection circle

  • i_vec (np.ndarray) – First basis vector (3D)

  • j_vec (np.ndarray) – Second basis vector (3D)

  • angles (np.ndarray) – Array of rotation angles (1D)

Returns:

(N, 3) array of positions, one per angle

Return type:

np.ndarray

pyfracval.utils.batch_check_overlaps_pca(coords_agg, radii_agg, candidate_positions, radius_new, tolerance)[source]#

Check overlap for batch of candidate positions (PCA).

Uses Numba parallel loops to evaluate multiple positions simultaneously.

Parameters:

  • coords_agg (np.ndarray) – Current aggregate coordinates (n_agg, 3)

  • radii_agg (np.ndarray) – Current aggregate radii (n_agg,)

  • candidate_positions (np.ndarray) – Batch of candidate positions to test (n_candidates, 3)

  • radius_new (float) – Radius of new particle

  • tolerance (float) – Overlap tolerance

Returns:

(n_candidates,) array of max overlap values for each position

Return type:

np.ndarray

pyfracval.utils.shuffle_array(arr, rng=None)[source]#

Randomly shuffle elements of a 1D array in-place (Fisher-Yates).

Modifies the input array directly. Mimics Fortran randsample behavior.

Parameters:

  • arr – The 1D NumPy array to shuffle.

  • rng (np.random.Generator | None, optional) – A NumPy Generator instance for reproducible randomness. If None, a fresh Generator is created.

Returns:

The input arr, modified in-place.

pyfracval.utils.sort_clusters(i_orden)[source]#

Sort cluster information array i_orden by cluster size (count).

Parameters:

i_orden (np.ndarray) – The Mx3 NumPy array [start_idx, end_idx, count].

Returns:

A new Mx3 array sorted by the ‘count’ column (column index 2).

Return type:

np.ndarray

Raises:

ValueError – If i_orden is not an Mx3 array.