"""Implements Particle-Cluster Aggregation (PCA) for initial subclusters.
This module follows the PCA stage of the FracVAL workflow described in
:cite:p:`Moran2019FracVAL`, with fractal scaling context from
:cite:p:`Filippov2000Tunable`.
.. note::
PCA solves the Gamma equation with particle **counts**, while
:mod:`pyfracval.cca` solves it with true **masses**. That is not an
oversight in either place. The scaling law supplying ``Rg`` is a
function of particle count, and mass-weighting is only consistent with
it when both merging bodies are aggregates the law describes - here
the second body is a single monomer, whose scaling-law ``Rg`` is
meaningless at n=1 while its mass can still be a large fraction of the
growing cluster's. Using masses here builds 1/150 subclusters where
counts build 93/150. See ``PCAggregator._gamma_calculation`` for
the full argument and the measurement.
Per-particle densities still matter in this stage: they feed
``self.mass``, ``self.m1`` and the running center of mass.
"""
import logging
import numpy as np
from scipy.spatial import cKDTree
from . import fractal, geometry, overlap, pca_kernels, utils
from .config import OrchestratorAlgorithmConfig
from .logs import TRACE_LEVEL_NUM
logger = logging.getLogger(__name__)
# Spatial indexing threshold: use k-d tree when aggregate has more than this many particles
KDTREE_THRESHOLD = 50 # Empirically tuned: tree overhead < linear search benefit
[docs]
class PCAggregator:
"""Performs Particle-Cluster Aggregation (PCA).
Builds a single cluster by sequentially adding individual primary
particles (monomers) to a growing aggregate, attempting to match
the target Df and kf via the Gamma_pc calculation at each step.
Includes overlap checking and rotation to find valid placements.
Parameters
----------
initial_radii : np.ndarray
1D array of radii for the primary particles to be aggregated.
df : float
Target fractal dimension for the aggregate.
kf : float
Target fractal prefactor for the aggregate.
tol_ov : float
Maximum allowable overlap fraction between particles.
Attributes
----------
N : int
Total number of particles to aggregate.
initial_mass : np.ndarray
Calculated initial masses corresponding to `initial_radii`.
coords : np.ndarray
Nx3 array storing coordinates of particles as they are placed.
radii : np.ndarray
N array storing radii of particles as they are placed.
mass : np.ndarray
N array storing masses of particles as they are placed.
n1 : int
Number of particles currently in the aggregate.
m1 : float
Mass of the current aggregate.
rg1 : float
Radius of gyration of the current aggregate.
cm : np.ndarray
3D center of mass of the current aggregate.
r_max : float
Maximum distance from CM to any particle center in the aggregate.
not_able_pca : bool
Flag indicating if the aggregation process failed.
"""
def __init__(
self,
initial_radii: np.ndarray,
df: float,
kf: float,
tol_ov: float,
rng: np.random.Generator | None = None,
algorithm_config: OrchestratorAlgorithmConfig | None = None,
densities: np.ndarray | None = None,
):
[docs]
self.N = len(initial_radii)
if self.N < 2:
raise ValueError("PCA requires at least 2 particles.")
self.algorithm_config: OrchestratorAlgorithmConfig = (
algorithm_config
if algorithm_config is not None
else OrchestratorAlgorithmConfig()
)
# Do NOT sort radii — use them in the order passed in (after shuffling in
# main_runner). The Fortran processes subcluster particles in whatever
# global order they appear; sorting (ascending or descending) creates
# pathological cases where the first two particles are either both tiny
# (making rg1 tiny and gamma_pc huge) or both huge (creating immediate
# overlap). Random order from shuffling avoids both extremes.
self.initial_radii = initial_radii.copy()
# Optional per-particle densities. Kept as a parallel array and
# swapped in lockstep with initial_radii/initial_mass below, so a
# particle's density always follows the particle rather than the
# index it happened to start at. None means uniform density.
self.initial_densities = fractal.resolve_densities(
densities, self.N, context="PCAggregator densities"
)
if self.initial_densities is not None:
self.initial_densities = self.initial_densities.copy()
# Calculate initial mass using utils consistently
[docs]
self.initial_mass = fractal.calculate_mass(
self.initial_radii, self.initial_densities
)
self.df = df
self.kf = kf
self.tol_ov = tol_ov
self._rng = rng if rng is not None else np.random.default_rng()
# State variables for the growing cluster
[docs]
self.coords = np.zeros((self.N, 3), dtype=float)
[docs]
self.radii = np.zeros(self.N, dtype=float)
[docs]
self.mass = np.zeros(self.N, dtype=float)
# Placement-ordered densities, filled alongside self.radii as each
# particle is accepted. This is what callers must read to know a
# finished subcluster's densities, since PCA reorders particles.
self.densities = (
np.zeros(self.N, dtype=float)
if self.initial_densities is not None
else None
)
[docs]
self.n1: int = 0 # Number of particles currently in the aggregate
[docs]
self.m1: float = 0.0 # Mass of the current aggregate
[docs]
self.rg1: float = 0.0 # Radius of gyration of the current aggregate
[docs]
self.cm = np.zeros(3) # Center of mass of the current aggregate
[docs]
self.r_max: float = 0.0 # Max distance from CM in the current aggregate
# Incremental Rg calculation: track sum of log(radii) for O(1) updates
self.sum_log_radii: float = 0.0 # Running sum of log(radii) for geo mean
[docs]
self.not_able_pca: bool = False
#: Populated when PCA gives up, describing *which* mechanism
#: failed and where. Read by Subclusterer to emit a structured
#: pca_failure event; None while the run is healthy.
self.failure_info: dict | None = None
def _random_point_sphere(self) -> tuple[float, float]:
"""Generates random angles (theta, phi) for a point on a sphere."""
u, v = self._rng.random(2)
# Use constant from config
theta = 2.0 * np.pi * u
phi = np.arccos(2.0 * v - 1.0)
return theta, phi
def _update_rg_incremental(self, new_radius: float) -> float:
"""Incrementally update radius of gyration when adding a particle.
This is O(1) vs O(n) for recalculating from scratch.
Uses running sum of log(radii) for geometric mean.
Parameters
----------
new_radius : float
Radius of the particle being added.
Returns
-------
float
Updated radius of gyration.
"""
if new_radius <= 1e-12:
# Skip near-zero radii (same filter as calculate_rg)
return self.rg1
# Update running sum of log(radii)
self.sum_log_radii += np.log(new_radius)
# Calculate geometric mean from running sum
# geo_mean_r = exp(sum_log_radii / n)
geo_mean_r = np.exp(self.sum_log_radii / self.n1)
# Calculate Rg using fractal scaling law
# Rg = geo_mean_r * (n / kf)^(1/df)
if self.kf > 0 and self.df > 0:
rg = geo_mean_r * (self.n1 / self.kf) ** (1.0 / self.df)
return max(rg, 0.0)
else:
return 0.0
def _first_two_monomers(self):
"""Places the first two monomers."""
if self.N < 2:
return # Should be caught by __init__ but safe check
self.radii[0] = self.initial_radii[0]
self.radii[1] = self.initial_radii[1]
if self.densities is not None and self.initial_densities is not None:
self.densities[0] = self.initial_densities[0]
self.densities[1] = self.initial_densities[1]
self.mass[0] = self.initial_mass[0]
self.mass[1] = self.initial_mass[1]
# Place first particle at origin
self.coords[0, :] = 0.0
# Place second particle touching the first (deterministically on X for consistency)
# distance = self.radii[0] + self.radii[1]
# self.coords[1, :] = [distance, 0.0, 0.0]
# Alternative: random orientation (original Fortran way)
theta, phi = self._random_point_sphere()
distance = self.radii[0] + self.radii[1]
self.coords[1, 0] = self.coords[0, 0] + distance * np.cos(theta) * np.sin(phi)
self.coords[1, 1] = self.coords[0, 1] + distance * np.sin(theta) * np.sin(phi)
self.coords[1, 2] = self.coords[0, 2] + distance * np.cos(phi)
self.n1 = 2
self.m1 = self.mass[0] + self.mass[1]
# Initialize incremental Rg calculation for first two particles
valid_radii = [r for r in [self.radii[0], self.radii[1]] if r > 1e-12]
if len(valid_radii) > 0:
self.sum_log_radii = sum(np.log(r) for r in valid_radii)
else:
self.sum_log_radii = 0.0
# Use utils for rg calculation (initial setup)
self.rg1 = fractal.calculate_rg(
self.radii[: self.n1], self.n1, self.df, self.kf
)
if self.m1 > geometry.FLOATING_POINT_ERROR: # Use utils tolerance
self.cm = (
self.coords[0] * self.mass[0] + self.coords[1] * self.mass[1]
) / self.m1
else:
self.cm = np.mean(self.coords[: self.n1], axis=0)
# Initial r_max (max distance from CM)
dist_0_cm = np.linalg.norm(self.coords[0] - self.cm)
dist_1_cm = np.linalg.norm(self.coords[1] - self.cm)
self.r_max = float(max(dist_0_cm, dist_1_cm))
def _gamma_calculation(
self,
m2: float,
rg2: float,
use_mass: bool = False,
) -> tuple[bool, float]:
r"""Distance Gamma_pc at which to place the next monomer.
Solves the **count** form of the Gamma equation, unlike CCA which
solves the **mass** form. That asymmetry is deliberate, is what
the original Fortran does, and is load-bearing rather than
cosmetic - see the module docstring for the short version and the
derivation below for why.
Background
----------
Both forms come from the same identity (:cite:p:`Moran2019FracVAL`
Eq. 6, derived in its Appendix A), relating two bodies being
merged to the result:
.. math::
m^2 R_g^2 = m\,(m_1 R_{g1}^2 + m_2 R_{g2}^2) + \Gamma^2 m_1 m_2
Substituting particle counts for the masses gives Filippov et al.
(2000) Eq. 7, which is exact only when every particle has the same
mass. FracVAL's contribution was to use true masses so that
*polydisperse* aggregates preserve Df and kf individually, and its
CCA stage does exactly that.
Why counts are nonetheless correct here
---------------------------------------
The identity is exact in the masses, but it is not solved in
isolation: :math:`R_g` for the *result* is supplied by the
scaling law
.. math:: R_{g} = a\,(n/k_f)^{1/D_f}
which is a function of the particle **count** :math:`n`, not of
mass. The two are mutually consistent only while every body
involved is an aggregate that the scaling law actually describes.
In CCA that holds: both bodies are clusters of many particles, so
their scaling-law :math:`R_g` is meaningful and the mass form is
both more faithful and better behaved.
In PCA it does not. The second body is a *single monomer*, for
which:
- its scaling-law :math:`R_g` is meaningless at :math:`n=1` (the
code substitutes the sphere's own :math:`\sqrt{3/5}\,r`), while
- its **mass** is not small at all. Under a wide size
distribution a single large monomer can carry a sizeable
fraction of the whole growing cluster's mass, since
:math:`m \propto r^3`.
So the mass-weighted left-hand side and the count-derived
:math:`R_g` on the right describe different objects, and
:math:`\Gamma^2` comes out of that mismatch either negative
(reported as "gamma not real") or far too large to admit any
candidate monomer on the existing cluster's surface. The
aggregation then stalls.
This is measurable rather than theoretical. At
:math:`\sigma_{p,geo}=1.9`, N=12 particles per subcluster, 150
seeds:
=================== ==========================
PCA Gamma form Subclusters built
=================== ==========================
counts (this one) 93/150 (62.0%)
masses 1/150 (0.7%)
=================== ==========================
Hence the Fortran's split between its two stages is not the
inconsistency it first appears to be, and this behaviour is fixed
per stage rather than exposed as a configuration flag someone
could set to a value that cannot work.
Scope
-----
This governs **only** which scalars enter the Gamma equation.
Everything mass-weighted in PCA's own bookkeeping - ``self.mass``,
``self.m1``, and the running center of mass - is computed from
real masses and is density-aware, so supplying per-particle
densities still shapes the resulting subcluster's geometry.
Parameters
----------
m2 : float
Mass of the monomer being added.
rg2 : float
Radius of gyration of that monomer, :math:`\sqrt{3/5}\,r`.
use_mass : bool, default False
Escape hatch for experiments only. See above for why the
default is what it is.
Returns
-------
tuple[bool, float]
``(gamma_real, gamma_pc)``; ``gamma_real`` is False when the
equation has no real solution for this pairing.
"""
return fractal.gamma_calculation(
self.m1,
self.rg1,
self.radii[: self.n1],
m2,
rg2,
np.array([self.initial_radii[self.n1]]),
self.df,
self.kf,
use_mass=use_mass,
all_radii=self.initial_radii,
)
# n1 = self.n1
# n2 = 1
# n3 = n1 + n2
# m1 = self.m1
# m3 = m1 + m2
# if heuristic:
# m1 = n1
# m2 = n2
# m3 = n3
# # Ensure index is valid before accessing initial_radii
# if n1 >= self.N:
# logger.error(
# f"Gamma calculation requested for particle index {n1} >= N ({self.N})"
# )
# return False, 0.0
# # Radii of particles already in cluster + the next one to be added
# combined_radii = np.concatenate((self.radii[:n1], [self.initial_radii[n1]]))
# rg3 = fractal.calculate_rg(combined_radii, n3, self.df, self.kf)
# # Heuristic from Fortran: ensure rg3 is not smaller than rg1
# # (avoids issues if rg calculation is noisy for small N)
# if self.rg1 > 0 and rg3 < self.rg1:
# logger.info(
# f"Gamma calc: Adjusted rg3 from {rg3:.2e} to match rg1 {self.rg1:.2e}"
# )
# rg3 = self.rg1
# gamma_pc = 0.0
# gamma_real = False
# term1 = (m3**2) * (rg3**2)
# term2 = m3 * (m1 * self.rg1**2 + m2 * rg2**2) # rg2 is for monomer
# denominator = m1 * m2
# radicand = term1 - term2
# try:
# gamma_pc = np.sqrt(radicand / denominator)
# gamma_real = True
# except (ValueError, ZeroDivisionError, OverflowError) as e:
# logger.warning(f"Gamma calculation internal failed: {e}")
# logger.warning(
# f"Gamma_pc calculation non-real or denominator zero: "
# f"n1={n1}, m1={m1:.2e}, rg1={self.rg1:.2e}, "
# f"m2={m2:.2e}, rg2={rg2:.2e}, "
# f"m3={m3:.2e}, rg3={rg3:.2e} -> "
# f"radicand={radicand:.2e}, denominator={denominator:.2e}"
# )
# gamma_real = False
# return gamma_real, gamma_pc
def _select_candidates(
self, radius_k: float, gamma_pc: float, gamma_real: bool
) -> tuple[np.ndarray, float]:
"""
Generates the list of candidate particles (indices within 0 to n1-1)
that monomer 'k' could stick to, based on Gamma_pc geometry.
Returns the candidate indices and Rmax (max distance from CM).
"""
candidates = []
r_max_current = 0.0
# If gamma is not real, sticking based on this criterion is impossible.
if not gamma_real:
logger.debug(
"Gamma_pc not real in PCA candidate selection. No candidates selected."
)
return np.array([], dtype=int), self.r_max # Return current r_max
# Rmax needs to be tracked based on *all* particles in the cluster
if self.n1 > 0:
distances_sq = np.sum((self.coords[: self.n1] - self.cm) ** 2, axis=1)
self.r_max = np.sqrt(np.max(distances_sq)) if distances_sq.size > 0 else 0.0
logger.debug(
f" _select_candidates: Checking N1={self.n1} particles against Gamma_pc={gamma_pc:.4f}, R_k={radius_k:.4f}"
)
# Use spatial indexing (k-d tree) for large aggregates, vectorized for small
if self.n1 > KDTREE_THRESHOLD:
# SPATIAL INDEXING: O(log n) candidate search with k-d tree
# Build k-d tree from current aggregate coordinates
tree = cKDTree(self.coords[: self.n1])
# Conservative search radius: particles within gamma_pc ± maximum possible radius_sum
# This includes all potentially valid candidates
max_radius_in_aggregate = (
np.max(self.radii[: self.n1]) if self.n1 > 0 else 0.0
)
search_radius = gamma_pc + radius_k + max_radius_in_aggregate
# Query k-d tree for particles within search radius of CM
candidate_indices = tree.query_ball_point(self.cm, search_radius)
# Refine candidates with exact geometric constraints
if len(candidate_indices) > 0:
# Compute distances from CM for candidates only (not all particles)
distances = np.linalg.norm(
self.coords[candidate_indices] - self.cm, axis=1
)
radii_candidates = self.radii[candidate_indices]
# Apply Fortran geometric conditions
radius_sum = radius_k + radii_candidates
lower_dist_bound = gamma_pc - radius_sum
upper_dist_bound = gamma_pc + radius_sum
radius_sum_check = radius_sum <= (
gamma_pc + geometry.FLOATING_POINT_ERROR
)
lower_bound_check = distances > (
lower_dist_bound - geometry.FLOATING_POINT_ERROR
)
upper_bound_check = distances <= (
upper_dist_bound + geometry.FLOATING_POINT_ERROR
)
all_checks = radius_sum_check & lower_bound_check & upper_bound_check
valid_candidate_mask = np.where(all_checks)[0]
# Map back to original indices
candidates = np.array(
[candidate_indices[i] for i in valid_candidate_mask]
)
# Debug logging (only if enabled)
if logger.isEnabledFor(logging.DEBUG):
logger.debug(
f" Spatial indexing: {len(candidate_indices)} in radius {search_radius:.2f}, "
f"{len(candidates)} passed geometric checks"
)
for i in candidates:
logger.debug(f" -> Candidate {i} ADDED (via k-d tree).")
else:
candidates = np.array([], dtype=int)
elif self.n1 > 0:
# VECTORIZED: O(n) search for small aggregates (faster due to lower overhead)
# Compute distances from CM for all particles (vectorized)
distances = np.linalg.norm(self.coords[: self.n1] - self.cm, axis=1)
# Get radii for all particles
radii_all = self.radii[: self.n1]
# Vectorized geometric checks
radius_sum = radius_k + radii_all
lower_dist_bound = gamma_pc - radius_sum
upper_dist_bound = gamma_pc + radius_sum
# Apply all three Fortran conditions (vectorized)
radius_sum_check = radius_sum <= (gamma_pc + geometry.FLOATING_POINT_ERROR)
lower_bound_check = distances > (
lower_dist_bound - geometry.FLOATING_POINT_ERROR
)
upper_bound_check = distances <= (
upper_dist_bound + geometry.FLOATING_POINT_ERROR
)
# Combine all conditions
all_checks = radius_sum_check & lower_bound_check & upper_bound_check
# Get indices of candidates
candidates = np.where(all_checks)[0]
# Debug logging (only if enabled)
if logger.isEnabledFor(logging.DEBUG):
for i in candidates:
logger.debug(
f" Cand i={i}: Dist={distances[i]:.4f}, R_i={radii_all[i]:.4f} | "
f"Cond1 (Rk+Ri <= G): {radius_sum[i]:.4f} <= {gamma_pc:.4f}? -> {radius_sum_check[i]} | "
f"Cond2 (Dist > G-Rk-Ri): {distances[i]:.4f} > {lower_dist_bound[i]:.4f}? -> {lower_bound_check[i]} | "
f"Cond3 (Dist <= G+Rk+Ri): {distances[i]:.4f} <= {upper_dist_bound[i]:.4f}? -> {upper_bound_check[i]}"
)
logger.debug(f" -> Candidate {i} ADDED.")
else:
candidates = np.array([], dtype=int)
logger.debug(
f"PCA selecting candidates for radius {radius_k:.2f} (gamma={gamma_pc:.3f}): Found {len(candidates)} candidates from {self.n1} particles."
)
return candidates, self.r_max # Return updated r_max
def _search_and_select_candidate(
self,
k: int,
considered_indices: list[int],
force_swap: bool = False,
tried_swaps: set[int] | None = None,
) -> tuple[int, float, float, bool, float, np.ndarray]:
"""
Handles the complex logic of selecting a candidate, potentially swapping
monomer 'k' with another if the initial attempt yields no candidates.
Corresponds roughly to the loop calling `Search_list` and `Random_select_list`.
Args:
k: The index of the particle being added to the aggregate.
considered_indices: List of particle indices already successfully placed.
force_swap: If True, forces a particle swap even if candidates exist.
This is needed when all candidates fail the overlap check.
tried_swaps: Persistent set of monomer indices already tried at position k.
When provided (not None) it is updated in-place so the caller
can pass the same set on subsequent calls and avoid retrying
the same swaps. If None a fresh set is created (legacy use).
Returns:
tuple: (selected_idx, m2, rg2, gamma_real, gamma_pc, candidate_list)
Returns -1 for selected_idx if no candidate found after all attempts.
candidate_list is the list of indices (0 to n1-1) found for the *final* particle k.
"""
available_monomers = list(range(k, self.N)) # Indices of unprocessed monomers
if tried_swaps is None:
tried_swaps = {k} # Monomers tried *at position k* (fresh set)
else:
tried_swaps.add(k) # Ensure current k is always in the tried set
# Track whether we have already done the forced swap this call.
# After the first forced swap the remaining iterations should check
# candidates normally (just like in the Fortran Search_list loop).
swap_done = False
while True:
# --- Try with current monomer k ---
current_k_idx = k # The actual index in the original list being processed
current_k_radius = self.initial_radii[current_k_idx]
current_k_mass = self.initial_mass[current_k_idx]
# Rg of a single sphere = sqrt(3/5) * R => sqrt(0.6) * R
current_k_rg = np.sqrt(0.6) * current_k_radius
gamma_real, gamma_pc = self._gamma_calculation(current_k_mass, current_k_rg)
logger.debug(
f"PCA search k={k}: Radius={current_k_radius:.2f}, Gamma_real={gamma_real}, Gamma_pc={gamma_pc:.4f}"
)
candidates = np.array([], dtype=int)
if gamma_real:
# Rmax is updated inside _select_candidates
candidates, self.r_max = self._select_candidates(
current_k_radius, gamma_pc, gamma_real
)
logger.debug(
f"PCA search k={k}: Found {len(candidates)} candidates: {candidates}"
)
# force_swap only applies to the very first iteration (before any swap
# has been done). Once we have performed a swap the candidates should
# be evaluated normally so we don't keep swapping indefinitely.
need_swap = (force_swap and not swap_done) or len(candidates) == 0
if len(candidates) > 0 and not need_swap:
# Select one candidate randomly (will be used as starting point in run loop)
idx_in_candidates = int(self._rng.integers(len(candidates)))
selected_initial_candidate = candidates[idx_in_candidates]
logger.debug(
f"PCA search k={k}: Initial candidate {selected_initial_candidate} selected from {len(candidates)} options."
)
# Return the *full list* of candidates found for this k
return (
selected_initial_candidate,
current_k_mass,
current_k_rg,
gamma_real,
gamma_pc,
candidates, # Return the list of all candidates
)
else:
# Need to swap: either no candidates OR a forced initial swap
if force_swap and not swap_done:
logger.debug(
f"PCA search k={k}: force_swap=True - will swap particle even though {len(candidates)} candidates exist."
)
else:
logger.debug(
f"PCA search k={k}: No candidates found or gamma not real. Looking for swap."
)
swap_done = (
True # mark that we have performed (or are about to perform) a swap
)
# --- No candidates: Try swapping k with an untried, available monomer ---
# Find monomers eligible for swapping (not k itself, not already tried at pos k,
# and not already successfully placed in the aggregate)
eligible_for_swap = [
idx
for idx in available_monomers
if idx not in tried_swaps and idx not in considered_indices
]
if not eligible_for_swap:
# No more monomers to swap with
logger.warning(
f"PCA k={k}: No candidates found and no more available monomers to swap with."
)
return (
-1, # Indicate failure
current_k_mass,
current_k_rg,
gamma_real,
gamma_pc,
np.array([], dtype=int), # Empty candidate list
)
# Select a random monomer to swap with k
swap_idx_in_eligible = int(self._rng.integers(len(eligible_for_swap)))
swap_target_original_idx = eligible_for_swap[swap_idx_in_eligible]
# Store original values for logging
radius_before = self.initial_radii[k]
radius_after = self.initial_radii[swap_target_original_idx]
logger.info(
f" PCA k={k}: SWAP - Particle radius {radius_before:.2f} → {radius_after:.2f} "
f"(swapping with index {swap_target_original_idx})"
)
# Perform the swap in the initial_radii and initial_mass arrays
# Swap particle at index k with particle at swap_target_original_idx
self.initial_radii[k], self.initial_radii[swap_target_original_idx] = (
self.initial_radii[swap_target_original_idx],
self.initial_radii[k],
)
self.initial_mass[k], self.initial_mass[swap_target_original_idx] = (
self.initial_mass[swap_target_original_idx],
self.initial_mass[k],
)
if self.initial_densities is not None:
(
self.initial_densities[k],
self.initial_densities[swap_target_original_idx],
) = (
self.initial_densities[swap_target_original_idx],
self.initial_densities[k],
)
# Note: The particle originally at k is now at swap_target_original_idx
# The particle originally at swap_target_original_idx is now at k
# Mark the monomer *originally* at swap_target_original_idx as having been tried *at position k*
tried_swaps.add(swap_target_original_idx)
# Loop continues, recalculating gamma/candidates with the new monomer now at index k
logger.debug(
f"PCA search k={k}: Swapped monomer from index {swap_target_original_idx} into position {k}. Retrying gamma/candidate search."
)
# The state of self.initial_radii/mass at index k has changed, restart loop
def _sticking_process(
self, k: int, selected_idx: int, gamma_pc: float
) -> tuple[np.ndarray | None, float, np.ndarray, np.ndarray, np.ndarray]:
"""
Calculates the geometric parameters for placing monomer k based on
intersection of two spheres.
Sphere 1: Center = selected particle coords, Radius = R_sel + R_k
Sphere 2: Center = CM of aggregate, Radius = Gamma_pc
Returns:
tuple: (coord_k_initial, theta_a, vec_0, i_vec, j_vec) or (None, ...) if intersection fails.
coord_k_initial is *one* point on the intersection circle.
Other return values define the circle for rotation (_reintento).
"""
if selected_idx < 0 or selected_idx >= self.n1:
logger.error(
f"Sticking process called with invalid selected_idx={selected_idx} (n1={self.n1})"
)
return None, 0.0, np.zeros(4), np.zeros(3), np.zeros(3)
if k < 0 or k >= self.N:
logger.error(f"Sticking process called with invalid k={k} (N={self.N})")
return None, 0.0, np.zeros(4), np.zeros(3), np.zeros(3)
coord_sel = self.coords[selected_idx]
radius_sel = self.radii[selected_idx]
# Use initial radius of particle k before it's placed
radius_k = self.initial_radii[k]
# Define the two spheres for intersection
sphere1_center = coord_sel
sphere1_radius = radius_sel + radius_k
sphere1 = np.concatenate((sphere1_center, [sphere1_radius]))
sphere2_center = self.cm # Use current aggregate CM
sphere2_radius = gamma_pc
sphere2 = np.concatenate((sphere2_center, [sphere2_radius]))
intersection_valid = False
try:
# This utility finds the circle and returns *one* random point on it
x_k, y_k, z_k, theta_a, vec_0, i_vec, j_vec, intersection_valid = (
geometry.two_sphere_intersection(sphere1, sphere2, rng=self._rng)
)
except Exception as e:
logger.error(
f"Error during two_sphere_intersection call: {e}", exc_info=True
)
intersection_valid = False # Ensure it's false on exception
if not intersection_valid:
logger.warning(
f"PCA sticking sphere intersection failed for k={k}, sel={selected_idx}."
)
# Log details to help diagnose intersection failures
dist_centers = np.linalg.norm(sphere1_center - sphere2_center)
radius_sum = sphere1_radius + sphere2_radius
radius_diff = abs(sphere1_radius - sphere2_radius)
logger.warning(
f" Intersection Fail Details: Center1={sphere1_center}, R1={sphere1_radius:.4f} | "
f"Center2={sphere2_center}, R2={sphere2_radius:.4f} | "
f"Dist={dist_centers:.4f}, R1+R2={radius_sum:.4f}, |R1-R2|={radius_diff:.4f}"
)
# Check conditions violated by intersection check in utils
if dist_centers > radius_sum + geometry.FLOATING_POINT_ERROR:
logger.warning(" -> Spheres too far apart.")
if dist_centers < radius_diff - geometry.FLOATING_POINT_ERROR:
logger.warning(" -> Sphere contained.")
if (
dist_centers < geometry.FLOATING_POINT_ERROR
and abs(radius_diff) < geometry.FLOATING_POINT_ERROR
):
logger.warning(" -> Spheres coincide.")
return None, 0.0, np.zeros(4), np.zeros(3), np.zeros(3) # Indicate failure
# Return the initial point found and the circle parameters for rotation
coord_k_initial = np.array([x_k, y_k, z_k])
return coord_k_initial, theta_a, vec_0, i_vec, j_vec
# Remove internal overlap check, use overlap.calculate_max_overlap_pca
# def _overlap_check(self, k: int) -> float: ...
def _reintento(
self,
k: int,
vec_0: np.ndarray,
i_vec: np.ndarray,
j_vec: np.ndarray,
attempt: int = 0,
) -> tuple[np.ndarray, float]:
"""
Calculates a *new* point on the intersection circle defined by
vec_0 (center, radius) and basis vectors i_vec, j_vec.
This is used to rotate monomer k to try and resolve overlaps.
Uses Fibonacci spiral sampling for optimal angular coverage, which avoids
redundant angle sampling and provides better geometric exploration than
pure random sampling.
Args:
k: Particle index
vec_0: [x0, y0, z0, r0] - center and radius of intersection circle
i_vec: First basis vector for the circle plane
j_vec: Second basis vector for the circle plane
attempt: Rotation attempt number (0-indexed) for Fibonacci spiral
Returns:
tuple: (coord_k_new, theta_a_new) - the new coordinates and the angle used.
"""
x0, y0, z0, r0 = vec_0
# If radius of intersection is near zero (touching point case), rotation is meaningless
if r0 < geometry.FLOATING_POINT_ERROR:
logger.debug(
f"Reintento k={k}: Intersection radius near zero, no rotation possible."
)
# Return the center of the 'circle' (the touch point)
return np.array([x0, y0, z0]), 0.0
# Generate angle using Fibonacci spiral for optimal coverage
# Golden ratio provides quasi-uniform distribution without repetition
golden_ratio = (1.0 + np.sqrt(5.0)) / 2.0
theta_a_new = 2.0 * np.pi * attempt / golden_ratio
# Calculate new position using the circle equation
coord_k_new = np.zeros(3)
coord_k_new = (
np.array([x0, y0, z0])
+ r0 * np.cos(theta_a_new) * i_vec
+ r0 * np.sin(theta_a_new) * j_vec
)
# coord_k_new[0] = (
# x0
# + r0 * np.cos(theta_a_new) * i_vec[0]
# + r0 * np.sin(theta_a_new) * j_vec[0]
# )
# coord_k_new[1] = (
# y0
# + r0 * np.cos(theta_a_new) * i_vec[1]
# + r0 * np.sin(theta_a_new) * j_vec[1]
# )
# coord_k_new[2] = (
# z0
# + r0 * np.cos(theta_a_new) * i_vec[2]
# + r0 * np.sin(theta_a_new) * j_vec[2]
# )
return coord_k_new, theta_a_new
def _pca_coarse_scan(
self,
k: int,
vec_0: np.ndarray,
i_vec: np.ndarray,
j_vec: np.ndarray,
n_coarse: int = 20,
) -> tuple:
"""Phase 1 of bisection: evaluate n_coarse Fibonacci-spiral positions.
Scans steps 1..n_coarse on the intersection circle and returns enough
information for the bisection phase to home in on the best bracket.
Parameters
----------
k : int
Particle index being placed.
vec_0 : np.ndarray
[x0, y0, z0, r0] — intersection circle centre and radius.
i_vec, j_vec : np.ndarray
Basis vectors of the intersection circle plane.
n_coarse : int
Number of Fibonacci steps to evaluate (default 20).
Returns
-------
tuple: (found, best_coord, best_overlap, best_angle,
left_angle, right_angle, steps_used)
found : bool — True if a valid position was found early.
best_coord : np.ndarray — coords of the best position seen.
best_overlap : float — overlap at best_coord.
best_angle : float — angle (radians) at best_coord.
left_angle : float — left bracket angle for bisection.
right_angle : float — right bracket angle for bisection.
steps_used : int — number of Fibonacci steps consumed.
"""
golden_ratio = (1.0 + np.sqrt(5.0)) / 2.0
best_overlap = np.inf
best_coord = self.coords[k].copy()
best_angle = 0.0
best_step = 1
# Track per-step overlaps so we can find the bracket neighbours
step_overlaps: list[float] = []
step_angles: list[float] = []
for step in range(1, n_coarse + 1):
coord_new, angle = self._reintento(k, vec_0, i_vec, j_vec, attempt=step)
self.coords[k] = coord_new
ov = overlap.calculate_max_overlap_pca_auto(
self.coords[: self.n1],
self.radii[: self.n1],
self.coords[k],
self.radii[k],
tolerance=self.tol_ov,
)
step_overlaps.append(ov)
step_angles.append(angle)
if ov < best_overlap:
best_overlap = ov
best_coord = coord_new.copy()
best_angle = angle
best_step = step
if ov <= self.tol_ov:
# Valid position found — return early
return (
True,
best_coord,
best_overlap,
best_angle,
best_angle,
best_angle,
step,
)
# Build bracket: find left/right neighbours of best_step with higher overlap
# Left neighbour: the step just before best_step (wrap to n_coarse if step==1)
left_idx = (best_step - 2) % n_coarse # step is 1-indexed, list is 0-indexed
right_idx = best_step % n_coarse # step after best_step (mod wraps)
left_angle = step_angles[left_idx]
right_angle = step_angles[right_idx]
# Restore coords to best position before returning
self.coords[k] = best_coord
return (
False,
best_coord,
best_overlap,
best_angle,
left_angle,
right_angle,
n_coarse,
)
def _true_overlap_at(self, k: int, threshold: float) -> float:
"""Overlap of particle ``k`` re-evaluated with early exit at ``threshold``.
``calculate_max_overlap_pca_auto`` stops at the first pair exceeding
the ``tolerance`` it is handed and returns *that* pair's overlap.
The returned value is therefore only a lower bound on the true
maximum once it exceeds that tolerance - which is exactly the
regime the adaptive-tolerance path operates in, since it triggers
on values already above ``tol_ov``.
Comparing such a lower bound against the 10x-larger ``relaxed_tol``
silently accepted placements whose real worst-case overlap was
enormous: the scan would return, say, 2.6e-6 from the first
offending pair (under relaxed_tol, so accepted) while a later pair
overlapped by 0.43. That produced "successful" subclusters
containing deeply interpenetrating particles, which CCA then
carried into the final aggregate untouched - it only ever checks
cluster-against-cluster, never within a cluster. See
docs/source/catalog_overlap_leak.md.
Re-running the scan with early exit at the threshold actually being
compared against makes the comparison sound: either no early exit
happens and the result is the true maximum, or it exceeds
``threshold`` and the placement is rejected regardless.
"""
return overlap.calculate_max_overlap_pca_auto(
self.coords[: self.n1],
self.radii[: self.n1],
self.coords[k],
self.radii[k],
tolerance=threshold,
)
def _pca_bisection(
self,
k: int,
vec_0: np.ndarray,
i_vec: np.ndarray,
j_vec: np.ndarray,
best_coord: np.ndarray,
best_overlap: float,
left_angle: float,
right_angle: float,
n_bisect: int = 15,
) -> tuple:
"""Phase 2 of bisection: binary search within the bracket from Phase 1.
Evaluates the midpoint of [left_angle, right_angle] and narrows the
bracket toward the half that has lower overlap, homing in on the
minimum-overlap arc.
Parameters
----------
k : int
Particle index being placed.
vec_0, i_vec, j_vec : np.ndarray
Intersection circle geometry (same as passed to coarse scan).
best_coord : np.ndarray
Current best-known coordinates for particle k (from coarse scan).
best_overlap : float
Overlap at best_coord.
left_angle, right_angle : float
Angular bracket around the best coarse-scan position (radians).
n_bisect : int
Maximum bisection iterations (default 15).
Returns
-------
tuple: (found, best_coord, best_overlap, steps_used)
found : bool — True if a valid position was found.
best_coord : np.ndarray — coordinates of the best position seen.
best_overlap : float — overlap at best_coord.
steps_used : int — number of bisection steps consumed.
"""
x0, y0, z0, r0 = vec_0
for step in range(n_bisect):
mid_angle = (left_angle + right_angle) / 2.0
coord_mid = (
np.array([x0, y0, z0])
+ r0 * np.cos(mid_angle) * i_vec
+ r0 * np.sin(mid_angle) * j_vec
)
self.coords[k] = coord_mid
ov = overlap.calculate_max_overlap_pca_auto(
self.coords[: self.n1],
self.radii[: self.n1],
self.coords[k],
self.radii[k],
tolerance=self.tol_ov,
)
if ov < best_overlap:
best_overlap = ov
best_coord = coord_mid.copy()
if ov <= self.tol_ov:
self.coords[k] = best_coord
return (True, best_coord, best_overlap, step + 1)
# Narrow bracket: move the endpoint whose midpoint gave lower overlap
# toward mid. Since we don't know the landscape shape, we simply
# split the larger half (standard bisection on the half-interval
# closest to mid).
mid_left = (left_angle + mid_angle) / 2.0
coord_mid_left = (
np.array([x0, y0, z0])
+ r0 * np.cos(mid_left) * i_vec
+ r0 * np.sin(mid_left) * j_vec
)
self.coords[k] = coord_mid_left
ov_left = overlap.calculate_max_overlap_pca_auto(
self.coords[: self.n1],
self.radii[: self.n1],
self.coords[k],
self.radii[k],
tolerance=self.tol_ov,
)
if ov_left < best_overlap:
best_overlap = ov_left
best_coord = coord_mid_left.copy()
if ov_left <= self.tol_ov:
self.coords[k] = best_coord
return (True, best_coord, best_overlap, step + 1)
if ov_left < ov:
right_angle = mid_angle
else:
left_angle = mid_angle
self.coords[k] = best_coord
return (False, best_coord, best_overlap, n_bisect)
[docs]
def run(self) -> np.ndarray | None:
"""Run the complete PCA process for all N particles.
Sequentially adds particles from index 2 to N-1. For each particle k,
it calculates Gamma_pc, finds potential sticking partners (`candidates`)
in the existing aggregate (0..k-1), potentially swaps particle k with
an unused one if no candidates are found initially.
It then attempts to stick particle k to each candidate partner,
calculating an initial position based on sphere intersections defined
by Gamma_pc. If overlap occurs, it rotates particle k around the
intersection circle (`_reintento`) up to `max_rotations` times.
If a non-overlapping position is found for any candidate, the particle
is successfully added, and aggregate properties are updated. If all
candidates and all rotations fail for a particle k, or if the initial
search/swap fails, the aggregation stops, `not_able_pca` is set True,
and None is returned.
Returns
-------
np.ndarray | None
An Nx4 NumPy array [X, Y, Z, R] of the final aggregate if
successful, otherwise None.
"""
if self.N < 2:
return None
self._first_two_monomers()
considered_indices = list(range(self.n1))
for k in range(self.n1, self.N):
logger.debug(f"--- PCA Step: Aggregating particle k={k} ---")
# --- Outer loop to allow re-searching/swapping if all candidates fail overlap ---
search_attempt = 0
max_search_attempts = self.N
sticking_successful = False
# Recorded for the failure event: how many partners the last
# search actually offered, which separates "nothing to try"
# from "tried several and all overlapped".
last_candidate_count = 0
# Persistent set of monomers already tried at position k.
# Must survive across search attempts so that each retry swaps in a
# DIFFERENT monomer (mirrors Fortran's `considerados` array which
# accumulates tried particles across all Search_list calls for a
# given k step).
tried_swaps_for_k: set[int] = set()
while not sticking_successful and search_attempt < max_search_attempts:
search_attempt += 1
logger.debug(f"PCA k={k}: Search/Swap Attempt #{search_attempt}")
# --- Perform Search/Swap ---
# Force particle swap on retry attempts (when previous candidates failed)
force_swap = search_attempt > 1
search_result = self._search_and_select_candidate(
k,
considered_indices,
force_swap=force_swap,
tried_swaps=tried_swaps_for_k,
)
(
initial_candidate_idx,
m2,
rg2,
gamma_real,
gamma_pc,
candidates_list,
) = search_result
# Check if search failed completely (no candidates even after swaps)
# Don't need len(candidates_list) check here
if initial_candidate_idx < 0 or not gamma_real:
logger.error(
f"PCA failed Search/Swap for k={k} (Attempt {search_attempt}). No valid gamma/candidates found even after swaps."
)
# Structured record of *which* PCA mechanism failed.
# "no candidate at a workable Gamma distance" and
# "candidates existed but all overlapped" have
# different causes and different fixes, and a failure
# taxonomy that cannot separate them is not much use.
self.failure_info = {
"particle_index": int(k),
"reason": (
"gamma_not_real" if not gamma_real else "no_candidates"
),
"search_attempts": int(search_attempt),
"n_candidates": 0,
"gamma_real": bool(gamma_real),
"gamma_pc": float(gamma_pc),
}
self.not_able_pca = True
return None # Cannot continue if search itself fails
# Store radius/mass for the current particle at index k
radius_k_current = self.initial_radii[k]
mass_k_current = self.initial_mass[k]
# --- Try Sticking with Found Candidates ---
if len(candidates_list) == 0:
logger.debug(
f"PCA k={k}, Attempt {search_attempt}: Search yielded Gamma but no candidates list. Retrying search/swap."
)
# Force the outer while loop to continue (effectively re-swaps)
# No need to do anything else, the while loop condition handles it
# This case might happen if _select_candidates fails geometrically
# even if gamma was real after a swap.
continue # Go to next iteration of the outer while loop
candidates_to_try = utils.shuffle_array(
candidates_list.copy(), rng=self._rng
)
last_candidate_count = len(candidates_to_try)
logger.debug(
f"PCA k={k}, Attempt {search_attempt}: Trying {len(candidates_to_try)} candidates: {candidates_to_try}"
)
all_candidates_failed_overlap = True # Assume failure until success
for current_selected_idx in candidates_to_try:
logger.debug(
f"PCA k={k}: Trying candidate partner index {current_selected_idx}"
)
stick_result = self._sticking_process(
k, current_selected_idx, gamma_pc
)
if stick_result is None or stick_result[0] is None:
logger.debug(
f" PCA k={k}, cand={current_selected_idx}: Sticking geometry failed."
)
continue # Try next candidate
coord_k_initial, theta_a, vec_0, i_vec, j_vec = stick_result
self.coords[k] = coord_k_initial
self.radii[k] = radius_k_current
self.mass[k] = mass_k_current
if (
self.densities is not None
and self.initial_densities is not None
):
self.densities[k] = self.initial_densities[k]
cov_max = overlap.calculate_max_overlap_pca_auto(
self.coords[: self.n1],
self.radii[: self.n1],
self.coords[k],
self.radii[k],
tolerance=self.tol_ov,
)
logger.debug(
f" PCA k={k}, cand={current_selected_idx}: Initial overlap = {cov_max:.4e}"
)
intento = 0
max_rotations = 360
adaptive_tol_threshold = (
180 # Relax tolerance after this many attempts
)
relaxed_tol = 1.0e-5 # Relaxed tolerance (10x more lenient)
used_adaptive_tol = False
# Check TRACE logging once (optimization: avoid check on every rotation)
trace_enabled = logger.isEnabledFor(TRACE_LEVEL_NUM)
# Choose rotation strategy based on configuration
if self.algorithm_config.use_batch_rotation:
# Batch rotation evaluation (Phase 3 - experimental, slower for N<1000)
batch_size = self.algorithm_config.rotation_batch_size
golden_ratio = (1.0 + np.sqrt(5.0)) / 2.0
# If intersection radius is near zero, skip rotation attempts
if vec_0[3] < geometry.FLOATING_POINT_ERROR:
logger.debug(
f" PCA k={k}, cand={current_selected_idx}: Intersection radius near zero, no rotation needed."
)
# cov_max already computed above, will be checked later
else:
# Batch rotation loop
while cov_max > self.tol_ov and intento < max_rotations:
# Determine batch range
batch_start = intento
batch_end = min(intento + batch_size, max_rotations)
batch_count = batch_end - batch_start
if batch_count == 0:
break
# Generate batch of angles using Fibonacci spiral
attempts = np.arange(batch_start, batch_end)
angles = 2.0 * np.pi * attempts / golden_ratio
# Calculate all positions in batch (parallel)
candidate_positions = (
pca_kernels.batch_calculate_positions_pca(
vec_0, i_vec, j_vec, angles
)
)
# Check overlaps for all positions in batch (parallel)
overlaps = pca_kernels.batch_check_overlaps_pca(
self.coords[: self.n1],
self.radii[: self.n1],
candidate_positions,
self.radii[k],
self.tol_ov,
)
# Find first valid position (overlap <= tolerance)
valid_indices = np.where(overlaps <= self.tol_ov)[0]
if len(valid_indices) > 0:
# Found valid position
best_idx = valid_indices[0]
intento = batch_start + best_idx + 1
self.coords[k] = candidate_positions[best_idx]
cov_max = overlaps[best_idx]
if trace_enabled:
logger.log(
TRACE_LEVEL_NUM,
f" PCA k={k}, cand={current_selected_idx}, Batch rotation {intento}: Found valid position with overlap={cov_max:.4e}",
)
break # Exit rotation loop
else:
# No valid position in this batch
# Use best (minimum overlap) from batch
best_idx = np.argmin(overlaps)
intento = batch_start + best_idx + 1
self.coords[k] = candidate_positions[best_idx]
cov_max = overlaps[best_idx]
# Check adaptive tolerance
if (
intento >= adaptive_tol_threshold
and cov_max <= relaxed_tol
):
true_cov = self._true_overlap_at(k, relaxed_tol)
if true_cov > relaxed_tol:
intento = batch_end
continue
cov_max = true_cov
logger.info(
f" PCA k={k}, cand={current_selected_idx}: Accepting relaxed tolerance "
f"(overlap={cov_max:.4e} <= {relaxed_tol:.4e}) after {intento} rotations."
)
used_adaptive_tol = True
break
if trace_enabled:
logger.log(
TRACE_LEVEL_NUM,
f" PCA k={k}, cand={current_selected_idx}, Batch {batch_start}-{batch_end}: Best overlap={cov_max:.4e} at attempt {intento}",
)
# Continue to next batch
intento = batch_end
else:
# Sequential rotation with 3-phase bisection (PyFracVAL-cut)
# Phase 1: coarse Fibonacci scan (steps 1..N_COARSE)
N_COARSE = 20
N_BISECT = 15
(
found_coarse,
best_coord,
cov_max,
best_angle,
left_angle,
right_angle,
intento,
) = self._pca_coarse_scan(
k, vec_0, i_vec, j_vec, n_coarse=N_COARSE
)
if found_coarse:
# Valid position found in coarse scan — done
logger.debug(
f" PCA k={k}, cand={current_selected_idx}: "
f"Bisect Phase 1 found valid pos at step {intento}, "
f"overlap={cov_max:.4e}"
)
else:
# Phase 2: bisection within the bracket
(
found_bisect,
best_coord,
cov_max,
bisect_steps,
) = self._pca_bisection(
k,
vec_0,
i_vec,
j_vec,
best_coord,
cov_max,
left_angle,
right_angle,
n_bisect=N_BISECT,
)
intento += bisect_steps
if found_bisect:
logger.debug(
f" PCA k={k}, cand={current_selected_idx}: "
f"Bisect Phase 2 found valid pos after {bisect_steps} bisect steps, "
f"overlap={cov_max:.4e}"
)
else:
# Phase 3: fallback — continue uniform Fibonacci
# from step N_COARSE+1 up to max_rotations
fallback_start = N_COARSE + 1
while cov_max > self.tol_ov and intento < max_rotations:
intento += 1
coord_k_new, _ = self._reintento(
k,
vec_0,
i_vec,
j_vec,
attempt=fallback_start
+ (intento - N_COARSE - N_BISECT - 1),
)
self.coords[k] = coord_k_new
cov_max = overlap.calculate_max_overlap_pca_auto(
self.coords[: self.n1],
self.radii[: self.n1],
self.coords[k],
self.radii[k],
tolerance=self.tol_ov,
)
if trace_enabled:
ov_details = []
for idx_agg in range(self.n1):
ov_agg = 1 - (
np.linalg.norm(
self.coords[k]
- self.coords[idx_agg]
)
/ (self.radii[k] + self.radii[idx_agg])
)
ov_details.append(
f"vs{idx_agg}:{ov_agg:.2e}"
)
logger.log(
TRACE_LEVEL_NUM,
f" PCA k={k}, cand={current_selected_idx}, "
f"Fallback Rot {intento}: Overlap = {cov_max:.4e} "
f"({', '.join(ov_details)})",
)
# Adaptive tolerance: relax constraint after many attempts
if (
intento >= adaptive_tol_threshold
and cov_max <= relaxed_tol
):
true_cov = self._true_overlap_at(k, relaxed_tol)
if true_cov > relaxed_tol:
continue
cov_max = true_cov
logger.info(
f" PCA k={k}, cand={current_selected_idx}: Accepting relaxed tolerance "
f"(overlap={cov_max:.4e} <= {relaxed_tol:.4e}) after {intento} rotations."
)
used_adaptive_tol = True
break
if cov_max <= self.tol_ov or used_adaptive_tol:
logger.debug(
f"PCA k={k}: Sticking successful with cand {current_selected_idx} after {intento} rotations."
)
sticking_successful = True # Set flag for outer loop
all_candidates_failed_overlap = (
False # Mark success for this attempt
)
break # Exit the 'for current_selected_idx' loop
else:
logger.debug(
f" PCA k={k}, cand={current_selected_idx}: Failed overlap after {max_rotations} rotations."
)
# Continue to the next candidate in candidates_to_try
# --- After trying all candidates for this search attempt ---
if all_candidates_failed_overlap:
logger.warning(
f"PCA k={k}, Attempt {search_attempt}: All {len(candidates_to_try)} candidates failed overlap check. Retrying search/swap..."
)
# Reset temporary placement before potentially swapping particle k
self.coords[k] = 0.0
self.radii[k] = 0.0
self.mass[k] = 0.0
if self.densities is not None:
self.densities[k] = 0.0
# The outer `while not sticking_successful` loop will continue
# else: sticking_successful is True, outer while loop will exit
# --- After the outer while loop ---
if not sticking_successful:
# This happens if max_search_attempts was reached
logger.error(
f"PCA failed at k={k}. Could not find non-overlapping position "
f"after {max_search_attempts} search/swap attempts."
)
self.failure_info = {
"particle_index": int(k),
"reason": "all_candidates_overlapped",
"search_attempts": int(max_search_attempts),
"n_candidates": int(last_candidate_count),
"gamma_real": True,
"gamma_pc": float(gamma_pc),
}
self.not_able_pca = True
return None # Critical failure
# --- Update aggregate properties (only if sticking was successful) ---
self.n1 += 1
m_old = self.m1
self.m1 += self.mass[k] # Use mass that was set during successful attempt
if self.m1 > geometry.FLOATING_POINT_ERROR:
self.cm = (self.cm * m_old + self.coords[k] * self.mass[k]) / self.m1
else:
self.cm = np.mean(self.coords[: self.n1], axis=0)
# Update rg1 using the FULL subcluster geomean (Fortran line 145):
# rg1 = geomean(R_all) * (n1/kf)^(1/Df)
# This matches the Fortran exactly: after each step the growing
# aggregate's Rg is estimated from the theoretical fractal law
# using the global subcluster geometric mean, not just the placed
# particles. Using only the placed particles underestimates rg1,
# making gamma_pc too large and breaking condition 2.
all_geomean = np.exp(np.mean(np.log(self.initial_radii)))
self.rg1 = all_geomean * (self.n1 / self.kf) ** (1.0 / self.df)
# Keep sum_log_radii in sync (used by _update_rg_incremental if called elsewhere)
self.sum_log_radii = np.sum(np.log(self.initial_radii))
considered_indices.append(k)
logger.debug(
f"--- PCA Step: Successfully added particle k={k}. Aggregate size n1={self.n1} ---"
)
# --- End of k loop ---
# ... (final checks and return as before) ...
if self.not_able_pca:
return None
final_data = np.hstack(
(self.coords[: self.N], self.radii[: self.N].reshape(-1, 1))
)
if np.any(np.isnan(final_data)):
logger.error("NaN detected in final PCA data.")
self.not_able_pca = True
return None
logger.info(f"PCA run completed successfully for N={self.N} particles.")
return final_data