"""
Stacked extinction coefficient computation across multiple cameras/LEDSA runs.
Usage example::
from ledsa.analysis.ConfigDataStacked import ConfigDataStacked
from ledsa.analysis.StackedExtinctionCoefficients import StackedExtinctionCoefficients
cfg = ConfigDataStacked() # reads config_stacked.ini
for reference_property in cfg.reference_properties:
for arr_id in cfg.get_led_array_ids():
for channel in cfg.camera_channels:
sec = StackedExtinctionCoefficients(cfg, arr_id, channel, reference_property)
sec.calc_and_set_coefficients()
sec.save()
"""
from __future__ import annotations
import datetime
from pathlib import Path
from typing import Dict, List, Optional
import numpy as np
import pandas as pd
from ledsa.analysis.ConfigDataStacked import ConfigDataStacked
from ledsa.analysis.Experiment import Camera, Experiment, Layers
from ledsa.analysis.ExtinctionCoefficientsLinear import ExtinctionCoefficientsLinear
from ledsa.analysis.ExtinctionCoefficientsNonLinear import ExtinctionCoefficientsNonLinear
from ledsa.analysis.temperature_correction import TemperatureCorrection
from ledsa.core.file_handling import get_extinction_coefficients_file_path, read_experiment_times
from importlib.metadata import version
[docs]
class StackedExtinctionCoefficients:
"""
Joint inversion of extinction coefficients using ray data stacked across
multiple cameras (LEDSA simulation directories).
For a given virtual LED-array ID and colour channel the class:
1. Builds one :class:`~ledsa.analysis.Experiment.Experiment` per
(simulation, ``led_array_id``) pair.
2. Constructs a stacked system matrix by vertically concatenating the
per-experiment distance-per-layer matrices and writes it once per
virtual array (channel-independent) to ``<output_path>/``.
3. Computes the common time axis across all simulations, writes a
time-alignment report to ``<output_path>/``, and stores it for reuse.
4. Solves the augmented Beer–Lambert system with the configured solver.
5. Writes extinction coefficient CSVs to
``<output_path>/analysis/extinction_coefficients/<solver>/<reference_property>/``.
:ivar config: Stacked-analysis configuration.
:vartype config: ConfigDataStacked
:ivar led_array_id: Virtual LED-array ID (integer key in ``[led_arrays]``).
:vartype led_array_id: int
:ivar channel: Colour channel index.
:vartype channel: int
:ivar reference_property: Quantity the coefficients are computed from.
:vartype reference_property: str
:ivar layers: Shared layer discretisation.
:vartype layers: Layers
:ivar led_array_mapping: ``{sim_idx: led_array_id}`` for this virtual array.
:vartype led_array_mapping: dict[int, int]
:ivar sub_experiments: Per-simulation :class:`Experiment` objects.
:vartype sub_experiments: list[Experiment]
:ivar distances_per_led_and_layer: Stacked distance matrix (n_leds_total × n_layers).
:vartype distances_per_led_and_layer: np.ndarray
:ivar ref_intensities: Stacked reference intensities (n_leds_total,).
:vartype ref_intensities: np.ndarray
:ivar coefficients_per_image_and_layer: Computed extinction coefficients.
:vartype coefficients_per_image_and_layer: list[np.ndarray]
:ivar common_times: Experiment times [s] shared by all simulations.
:vartype common_times: np.ndarray or None
"""
def __init__(self, config: ConfigDataStacked, led_array_id: int, channel: int,
reference_property: Optional[str] = None) -> None:
"""
:param config: Parsed ``config_stacked.ini``.
:type config: ConfigDataStacked
:param led_array_id: Integer key defined in the ``[led_arrays]`` section.
:type led_array_id: int
:param channel: Colour channel (0 = R, 1 = G, 2 = B).
:type channel: int
:param reference_property: Quantity the coefficients are computed from. Defaults to the
first quantity the configuration names.
:type reference_property: str or None
"""
self.config = config
self.led_array_id = led_array_id
self.channel = channel
self.reference_property = (reference_property if reference_property is not None
else config.reference_property)
self.layers = Layers(config.num_layers, *config.domain_bounds)
self.led_array_mapping: Dict[int, int] = config.get_led_array_mapping(led_array_id)
# The factors belong to the LEDs, so every camera of this inversion is corrected with them
reduction_file = config.get_reduction_file(channel)
self.temperature_correction: Optional[TemperatureCorrection] = (
TemperatureCorrection.from_file(reduction_file) if reduction_file is not None else None
)
self.sub_experiments: List[Experiment] = []
for sim_idx, arr_id in sorted(self.led_array_mapping.items()):
cam_pos = config.get_camera_position(sim_idx)
exp = Experiment(
layers=self.layers,
led_array=arr_id,
camera=Camera(*cam_pos, pane_distance=config.get_pane_distance(sim_idx)),
path=config.get_simulation_path(sim_idx),
channel=channel,
merge_led_arrays='None',
)
self.sub_experiments.append(exp)
# sub_experiments follows sorted(self.led_array_mapping), so the simulation index of every
# experiment is the one its reference image settings have to be read from
self._sub_solvers: List = [
self._build_sub_solver(exp, sim_idx)
for sim_idx, exp in zip(sorted(self.led_array_mapping), self.sub_experiments)
]
self.distances_per_led_and_layer: np.ndarray = np.array([])
self.ref_intensities: np.ndarray = np.array([])
self.coefficients_per_image_and_layer: List[np.ndarray] = []
self.common_times: Optional[np.ndarray] = None
self._time_maps: Optional[List[Dict[int, float]]] = None # set by _compute_time_alignment
self.solver: str = config.solver
# ------------------------------------------------------------------
# Sub-solver factory
# ------------------------------------------------------------------
def _build_sub_solver(self, experiment: Experiment, sim_idx: int):
"""Create a single-experiment solver instance for *experiment*.
:param experiment: Experiment of the simulation the solver belongs to.
:type experiment: Experiment
:param sim_idx: Index of that simulation, used to read its own reference images.
:type sim_idx: int
"""
cfg = self.config
if cfg.solver == 'linear':
return ExtinctionCoefficientsLinear(
experiment=experiment,
reference_property=self.reference_property,
num_ref_imgs=cfg.get_num_ref_images(sim_idx),
ref_img_indices=cfg.get_ref_img_indices(sim_idx),
average_images=cfg.average_images,
n_summarize=cfg.n_summarize,
lambda_reg=cfg.lambda_reg,
temperature_correction=self.temperature_correction,
)
elif cfg.solver == 'nonlinear':
return ExtinctionCoefficientsNonLinear(
experiment=experiment,
reference_property=self.reference_property,
num_ref_imgs=cfg.get_num_ref_images(sim_idx),
ref_img_indices=cfg.get_ref_img_indices(sim_idx),
average_images=cfg.average_images,
n_summarize=cfg.n_summarize,
weighting_preference=cfg.weighting_preference,
weighting_curvature=cfg.weighting_curvature,
num_iterations=cfg.num_iterations,
temperature_correction=self.temperature_correction,
)
else:
raise ValueError(
f'Unknown solver {cfg.solver!r}. Choose "linear" or "nonlinear".'
)
# ------------------------------------------------------------------
# Main entry points
# ------------------------------------------------------------------
[docs]
def calc_and_set_coefficients(self) -> None:
"""
Compute extinction coefficients for all time-aligned images.
Populates :attr:`coefficients_per_image_and_layer` and
:attr:`common_times`.
"""
self._ensure_stacked_variables()
ref_prop = self.reference_property
for _, aligned_intensities in self._iter_aligned_images(ref_prop):
rel_intensities = aligned_intensities / self.ref_intensities
if self.config.solver == 'nonlinear' and self.coefficients_per_image_and_layer:
self._sub_solvers[0].coefficients_per_image_and_layer = (
self.coefficients_per_image_and_layer
)
sigmas = self._sub_solvers[0].calc_coefficients_of_img(rel_intensities)
self.coefficients_per_image_and_layer.append(sigmas)
[docs]
def calc_and_set_coefficients_mp(self, cores: Optional[int] = None) -> None:
"""Parallel computation using multiprocessing.
:param cores: Number of worker processes. Defaults to ``config.num_cores``.
:type cores: int or None
"""
from multiprocessing import Pool
if cores is None:
cores = self.config.num_cores
self._ensure_stacked_variables()
ref_prop = self.reference_property
all_rel_intensities = [
aligned / self.ref_intensities
for _, aligned in self._iter_aligned_images(ref_prop)
]
with Pool(processes=cores) as pool:
self.coefficients_per_image_and_layer = pool.map(
self._sub_solvers[0].calc_coefficients_of_img, all_rel_intensities
)
# ------------------------------------------------------------------
# Lazy initialisation
# ------------------------------------------------------------------
def _ensure_stacked_variables(self) -> None:
"""Load/compute distance matrix, time alignment, and reference intensities."""
if self.distances_per_led_and_layer.size == 0:
self._build_stacked_distance_matrix()
if self._time_maps is None:
self._compute_time_alignment()
if self.ref_intensities.size == 0:
self._build_stacked_ref_intensities()
def _build_stacked_distance_matrix(self) -> None:
"""
Populate each sub-solver's distance matrix, then stack them.
The stacked distance matrix is written once per virtual array
(channel-independent) to::
<output_path>/led_array_<id>_distances_per_layer.txt
A guard prevents redundant writes when multiple channels are processed.
"""
for sim_idx, (sub, arr_id) in enumerate(
zip(self._sub_solvers,
[self.led_array_mapping[i] for i in sorted(self.led_array_mapping)])
):
exp = sub.experiment
if exp.num_leds is None:
raise FileNotFoundError(
f"Could not load LED data for simulation {sim_idx} "
f"(led_array_id={arr_id}, path={exp.path}).\n"
f"Expected files:\n"
f" {exp.path / 'analysis' / f'led_array_indices_{arr_id:03d}.csv'}\n"
f" {exp.path / 'analysis' / 'led_search_areas_with_coordinates.csv'}\n"
"Make sure LEDSA steps 1-3 have been completed in that directory "
"and that the led_array_id in [led_arrays] is correct."
)
sub.set_all_member_variables(save_distances=False)
self.distances_per_led_and_layer = np.vstack(
[sub.distances_per_led_and_layer for sub in self._sub_solvers]
)
dist_file = self.config.output_path / f'led_array_{self.led_array_id}_distances_per_layer.txt'
if not dist_file.exists():
self.config.output_path.mkdir(parents=True, exist_ok=True)
np.savetxt(dist_file, self.distances_per_led_and_layer)
self._sub_solvers[0].distances_per_led_and_layer = self.distances_per_led_and_layer
def _build_stacked_ref_intensities(self) -> None:
"""Concatenate per-sub-solver reference intensities."""
self.ref_intensities = np.concatenate(
[sub.ref_intensities for sub in self._sub_solvers]
)
def _compute_time_alignment(self) -> None:
"""
Compute the common time axis across all simulations and write a
time-alignment report to::
<output_path>/time_alignment_led_array_<id>.txt
The report is written once per virtual array (channel-independent).
Sets :attr:`_time_maps` and :attr:`common_times`.
Images that were skipped during the LED intensity extraction have no data to invert, so
they take no part in the alignment and are left out of the report as well.
"""
tol = self.config.time_sync_tolerance
self._time_maps = [
self._load_time_map(sub.experiment, sub.average_images, sub.calculated_img_data)
for sub in self._sub_solvers
]
all_time_arrays = [np.array(list(tm.values())) for tm in self._time_maps]
common_times = self._intersect_time_axes(all_time_arrays, tol)
if len(common_times) == 0:
raise RuntimeError(
'No overlapping experiment times found across simulations for '
f'led_array_id={self.led_array_id}. '
'Check time_sync_tolerance in config_stacked.ini.'
)
self.common_times = common_times
print(
f' [array {self.led_array_id}] '
f'{len(common_times)} time-aligned images from '
f'{self.config.get_num_simulations()} simulations.'
)
report_file = self.config.output_path / f'time_alignment_led_array_{self.led_array_id}.txt'
if not report_file.exists():
self._write_time_alignment_report(report_file, self._time_maps, all_time_arrays, common_times, tol)
def _write_time_alignment_report(
self,
path: Path,
time_maps: List[Dict[int, float]],
all_time_arrays: List[np.ndarray],
common_times: np.ndarray,
tol: float,
) -> None:
"""
Write a plain-text time-alignment report.
Includes a per-image table showing, for every common time step, which
image ID was selected from each simulation and the corresponding time
difference.
:param path: Destination file path.
:param time_maps: ``{img_id: time}`` per simulation.
:param all_time_arrays: Flat time arrays per simulation.
:param common_times: Intersected common time axis.
:param tol: Matching tolerance used [s].
"""
self.config.output_path.mkdir(parents=True, exist_ok=True)
sim_indices = sorted(self.led_array_mapping)
lines = [
f'Time-alignment report — led_array_id {self.led_array_id}',
f'Generated : {datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S")}',
f'Tolerance : {tol} s',
'',
'Simulations',
'-----------',
]
for sim_idx, arr_id, times in zip(
sim_indices, self.led_array_mapping.values(), all_time_arrays
):
sim_path = self.config.get_simulation_path(sim_idx)
t_min = float(times.min()) if len(times) else float('nan')
t_max = float(times.max()) if len(times) else float('nan')
lines += [
f' sim {sim_idx} (led_array_id={arr_id})',
f' path : {sim_path}',
f' images : {len(times)}',
f' range : {t_min:.1f} s – {t_max:.1f} s',
]
lines += [
'',
'Alignment result',
'----------------',
f' Common images : {len(common_times)}',
]
if len(common_times):
lines += [
f' Time range : {float(common_times[0]):.1f} s – {float(common_times[-1]):.1f} s',
]
lines += ['', 'Per-simulation match rate', '-------------------------']
for sim_idx, times in zip(sim_indices, all_time_arrays):
matched = sum(1 for t in times if np.any(np.abs(common_times - t) <= tol))
pct = 100.0 * matched / len(times) if len(times) else 0.0
lines.append(f' sim {sim_idx} : {matched}/{len(times)} images matched ({pct:.1f} %)')
# ------------------------------------------------------------------
# Per-image matching table
# ------------------------------------------------------------------
# Build reverse map: time → img_id for each simulation
rev_maps: List[Dict[float, int]] = [
{t: img_id for img_id, t in tm.items()} for tm in time_maps
]
# Column widths
col_t_common = 14 # "common t / s"
col_per_sim = 26 # "img_id t_sim / s Δt / s" per simulation
# Header
sim_headers = ''.join(
f' {"sim " + str(si):^{col_per_sim}}'
for si in sim_indices
)
sub_headers = ''.join(
f' {"img_id":>8} {"t_sim / s":>10} {"Δt / s":>6}'
for _ in sim_indices
)
separator = '-' * (col_t_common + len(sim_indices) * (col_per_sim + 2))
lines += [
'',
'Per-image match table',
'---------------------',
f' {"common t / s":>{col_t_common}}' + sim_headers,
f' {" " * col_t_common}' + sub_headers,
' ' + separator,
]
for t_common in common_times:
row = f' {t_common:>{col_t_common}.2f}'
for tm in time_maps:
best_id = self._nearest_img_id(tm, t_common, tol)
if best_id is not None:
t_sim = tm[best_id]
dt = t_sim - t_common
row += f' {best_id:>8d} {t_sim:>10.2f} {dt:>+6.2f}'
else:
row += f' {"—":>8} {"—":>10} {"—":>6}'
lines.append(row)
path.write_text('\n'.join(lines) + '\n', encoding='utf-8')
print(f' Time-alignment report: {path}')
# ------------------------------------------------------------------
# Time-aligned image iteration
# ------------------------------------------------------------------
def _iter_aligned_images(self, reference_property: str):
"""
Yield ``(virtual_img_id, stacked_intensity_vector)`` for every
time-aligned image set, using the pre-computed :attr:`common_times`
and :attr:`_time_maps`.
"""
tol = self.config.time_sync_tolerance
for v_img_id, t_common in enumerate(self.common_times, start=1):
blocks: List[np.ndarray] = []
skip = False
for sub, time_map in zip(self._sub_solvers, self._time_maps):
orig_id = self._nearest_img_id(time_map, t_common, tol)
if orig_id is None:
print(
f' Warning: no image within {tol} s of t={t_common:.2f} s '
f'in simulation at {sub.experiment.path} — skipping time step.'
)
skip = True
break
try:
img_slice = sub.calculated_img_data.loc[orig_id]
except KeyError:
print(
f' Warning: img_id={orig_id} missing in calculated_img_data '
f'for {sub.experiment.path} — skipping time step.'
)
skip = True
break
blocks.append(img_slice[reference_property].to_numpy())
if skip:
continue
yield v_img_id, np.concatenate(blocks)
# ------------------------------------------------------------------
# Save
# ------------------------------------------------------------------
[docs]
def save(self) -> None:
"""
Write computed extinction coefficients to CSV.
Output path::
<output_path>/analysis/extinction_coefficients/<solver>/<ref_prop>/
extinction_coefficients_<solver>_channel_<ch>_<ref_prop>_led_array_<id>.csv
"""
if not self.coefficients_per_image_and_layer:
raise RuntimeError('No coefficients computed yet. Call calc_and_set_coefficients first.')
out_file = get_extinction_coefficients_file_path(self.solver, self.channel,
self.reference_property, self.led_array_id,
root=self.config.output_path)
out_file.parent.mkdir(parents=True, exist_ok=True)
layers = self.layers
sim_info = ' | '.join(
f'sim{idx}:arr{arr}' for idx, arr in sorted(self.led_array_mapping.items())
)
header = (
f'Stacked inversion | {sim_info} | {self._pane_info()}solver: {self.solver} | '
f'ref_prop: {self.reference_property} | '
f'ref_imgs: {self._ref_img_info()} | '
f'temperature_correction: {self._sub_solvers[0].temperature_correction_info()} | '
f'LEDSA {version("ledsa")}\n'
)
header += 'Layer_bottom_height[m],' + ','.join(map(str, layers.borders[:-1])) + '\n'
header += 'Layer_top_height[m],' + ','.join(map(str, layers.borders[1:])) + '\n'
header += 'Layer_center_height[m],' + ','.join(
f'{c:.4f}' for c in (layers.borders[:-1] + layers.borders[1:]) / 2
) + '\n'
header += 'Experiment_Time[s],layer0'
header += ''.join(f',layer{i + 1}' for i in range(layers.amount - 1))
times = (
self.common_times
if self.common_times is not None
else np.arange(len(self.coefficients_per_image_and_layer), dtype=float)
)
if len(times) != len(self.coefficients_per_image_and_layer):
raise ValueError(
f'Got {len(times)} common experiment times for '
f'{len(self.coefficients_per_image_and_layer)} computed images. The time axis and '
f'the analysis results of the simulations are out of sync.'
)
data = np.column_stack((times, self.coefficients_per_image_and_layer))
np.savetxt(out_file, data, delimiter=',', header=header)
print(f'Saved: {out_file}')
def __str__(self) -> str:
sim_info = ', '.join(
f'sim{idx}(arr{arr})' for idx, arr in sorted(self.led_array_mapping.items())
)
return (
f'StackedExtinctionCoefficients | led_array_id={self.led_array_id} | '
f'channel={self.channel} | solver={self.solver} | '
f'ref_prop={self.reference_property} | '
f'simulations=[{sim_info}] | '
f'n_leds={sum(e.num_leds for e in self.sub_experiments if e.num_leds is not None)} | '
f'n_layers={self.layers.amount}\n'
)
# ------------------------------------------------------------------
# Internal helpers
# ------------------------------------------------------------------
def _pane_info(self) -> str:
"""Header entry of the cameras behind a window pane, e.g. ``pane_distance: sim0 0.5 | ``.
Only a camera behind a pane is named, so the header of a run without panes stays as it was.
:return: The entry, ending in a separator, or an empty string if no camera is behind a pane.
:rtype: str
"""
distances = [
f'sim{idx} {exp.camera.pane_distance}'
for idx, exp in zip(sorted(self.led_array_mapping), self.sub_experiments)
if exp.camera.pane_distance is not None
]
return f'pane_distance: {", ".join(distances)} | ' if distances else ''
def _ref_img_info(self) -> str:
"""Describe the reference images of every simulation, e.g. ``sim0 first 10, sim1 [385, 386]``.
Every camera can use reference images of its own, so the result only states which ones
went into I₀ if it names them per simulation.
:return: One entry per simulation, in the order of the simulation indices.
:rtype: str
"""
descriptions = []
for sim_idx, sub in zip(sorted(self.led_array_mapping), self._sub_solvers):
if sub.ref_img_indices is None:
descriptions.append(f'sim{sim_idx} first {sub.num_ref_imgs}')
else:
descriptions.append(f'sim{sim_idx} {list(sub.ref_img_indices)}')
return ', '.join(descriptions)
@staticmethod
def _load_time_map(experiment: Experiment, average_images: bool,
img_data: Optional[pd.DataFrame] = None) -> Dict[int, float]:
"""
Read the experiment times of the images of one simulation.
:param experiment: Experiment whose image infos file is read.
:type experiment: Experiment
:param average_images: Whether the averaged image infos file is used.
:type average_images: bool
:param img_data: Intensity data of the simulation. If given, only the images it contains
are kept.
:type img_data: pd.DataFrame or None
:return: ``{img_id: experiment_time}`` of the images available for matching.
:rtype: dict[int, float]
"""
time_map = read_experiment_times(experiment.path, average_images)
if img_data is not None and not img_data.empty:
# Images without extraction results have no LED intensities and must not be matched
available_img_ids = set(img_data.index.get_level_values('img_id'))
time_map = {i: t for i, t in time_map.items() if i in available_img_ids}
return time_map
@staticmethod
def _intersect_time_axes(time_arrays: List[np.ndarray], tol: float) -> np.ndarray:
if not time_arrays:
return np.array([])
common = time_arrays[0]
for times in time_arrays[1:]:
mask = np.array([np.any(np.abs(times - t) <= tol) for t in common])
common = common[mask]
return common
@staticmethod
def _nearest_img_id(time_map: Dict[int, float], target: float, tol: float) -> Optional[int]:
best_id, best_diff = None, float('inf')
for img_id, t in time_map.items():
diff = abs(t - target)
if diff < best_diff:
best_diff = diff
best_id = img_id
return best_id if best_diff <= tol else None