# lec/manager.py
import logging
import shutil
import numpy as np
import numba as nb
from contextlib import ExitStack
from pathlib import Path
from oasislmf.pytools.summary.manager import SUMMARY_META_SIZE
import pyarrow as pa
import pyarrow.parquet as pq
from oasislmf.pytools.common.data import (DEFAULT_BUFFER_SIZE,
summary_stream_index_dtype, def_to_type_and_size)
from oasislmf.pytools.common.event_stream import MAX_LOSS_IDX, MEAN_IDX, NUMBER_OF_AFFECTED_RISK_IDX, SUMMARY_STREAM_ID, init_streams_in, mv_read
from oasislmf.pytools.common.input_files import PERIODS_FILE, occ_get, read_occurrence, read_periods, read_returnperiods
from oasislmf.pytools.lec.data import (AEP, AEPTVAR, AGG_FULL_UNCERTAINTY, AGG_SAMPLE_MEAN, AGG_WHEATSHEAF, AGG_WHEATSHEAF_MEAN,
OCC_FULL_UNCERTAINTY, OCC_SAMPLE_MEAN, OCC_WHEATSHEAF, OCC_WHEATSHEAF_MEAN, OEP, OEPTVAR,
OUTLOSS_DTYPE, EPT_dtype, EPT_fmt, EPT_headers, PSEPT_dtype, PSEPT_fmt, PSEPT_headers)
from oasislmf.pytools.lec.aggreports import AggReports, LecConfig, make_output_fn, output_for_summary_idx
from oasislmf.pytools.lec.utils import get_outloss_mean_idx, get_outloss_sample_idx
from oasislmf.pytools.utils import redirect_logging
[docs]
logger = logging.getLogger(__name__)
event_id_dtype, event_id_dtype_size = def_to_type_and_size('event_id')
item_id_dtype, item_id_dtype_size = def_to_type_and_size('item_id')
summary_id_dtype, summary_id_dtype_size = def_to_type_and_size('summary_id')
sidx_dtype, sidx_size = def_to_type_and_size('sidx')
loss_dtype, loss_dtype_size = def_to_type_and_size('loss')
@nb.njit(cache=True, error_model="numpy")
[docs]
def get_max_summary_id(file_handles):
"""Get max summary_id from all summary files
Args:
file_handles (List[np.memmap]): List of memmaps for summary files data
Returns:
max_summary_id (int): Max summary ID
"""
max_summary_id = 0
for fin in file_handles:
cursor = SUMMARY_META_SIZE
valid_buff = len(fin)
while cursor < valid_buff:
# header (event_id, summary_id, exposure_value)
_, cursor = mv_read(fin, cursor, event_id_dtype, event_id_dtype_size)
summary_id, cursor = mv_read(fin, cursor, summary_id_dtype, summary_id_dtype_size)
_, cursor = mv_read(fin, cursor, loss_dtype, loss_dtype_size)
max_summary_id = max(max_summary_id, summary_id)
while cursor < valid_buff:
sidx, cursor = mv_read(fin, cursor, sidx_dtype, sidx_size)
_, cursor = mv_read(fin, cursor, loss_dtype, loss_dtype_size)
if sidx == 0:
break
return max_summary_id
@nb.njit(cache=True, fastmath=True, error_model="numpy")
[docs]
def do_lec_output_agg_summary(
summary_id,
sidx,
loss,
filtered_occ_map,
outloss_mean,
row_used_mean,
outloss_sample,
row_used_sample,
num_sidxs,
max_summary_id,
):
"""Populate outloss_mean and outloss_sample with aggregate and max losses
Args:
summary_id (summary_id_dtype): summary_id
sidx (sidx_dtype): Sample ID
loss (loss_dtype): Loss value
filtered_occ_map (ndarray[occ_map_dtype]): Filtered numpy map of event_id, period_no, occ_date_id from the occurrence file_
outloss_mean (ndarray[OUTLOSS_DTYPE]): ndarray indexed by summary_id, period_no containing aggregate and max losses
row_used_mean (ndarray[bool]): bool mask for outloss_mean
outloss_sample (ndarray[OUTLOSS_DTYPE]): ndarray indexed by summary_id, sidx, period_no containing aggregate and max losses
row_used_sample (ndarray[bool]): bool mask for outloss_sample
num_sidxs (int): Number of sidxs to consider for outloss_sample
max_summary_id (int): Max summary ID
"""
for row in filtered_occ_map:
period_no = row["period_no"]
if sidx == MEAN_IDX:
idx = get_outloss_mean_idx(period_no, summary_id, max_summary_id)
outloss_mean[idx]["agg_out_loss"] += loss
if loss > outloss_mean[idx]["max_out_loss"]:
outloss_mean[idx]["max_out_loss"] = loss
row_used_mean[idx] = True
else:
idx = get_outloss_sample_idx(period_no, sidx, summary_id, num_sidxs, max_summary_id)
outloss_sample[idx]["agg_out_loss"] += loss
if loss > outloss_sample[idx]["max_out_loss"]:
outloss_sample[idx]["max_out_loss"] = loss
row_used_sample[idx] = True
@nb.njit(cache=True, error_model="numpy")
[docs]
def process_summary_entries(
fin,
offsets,
occ_csr,
use_return_period,
outloss_mean_s,
row_used_mean_s,
outloss_sample_s,
row_used_sample_s,
num_sidxs,
):
"""Process all indexed event blocks for one (summary_id, file) pair in a single call.
Eliminates per-event Python→numba overhead by looping over all offsets inside numba.
offsets should be sorted ascending for best OS page-cache utilisation.
"""
valid_buff = len(fin)
for offset in offsets:
cursor = offset
event_id, cursor = mv_read(fin, cursor, event_id_dtype, event_id_dtype_size)
_, cursor = mv_read(fin, cursor, summary_id_dtype, summary_id_dtype_size) # summary_id known from idx
_, cursor = mv_read(fin, cursor, loss_dtype, loss_dtype_size) # expval
filtered_occ_map = occ_get(occ_csr, event_id)
if len(filtered_occ_map) == 0:
continue
while cursor < valid_buff:
sidx, cursor = mv_read(fin, cursor, sidx_dtype, sidx_size)
loss, cursor = mv_read(fin, cursor, loss_dtype, loss_dtype_size)
if sidx == 0:
break
if sidx == NUMBER_OF_AFFECTED_RISK_IDX or sidx == MAX_LOSS_IDX:
continue
if loss > 0 or use_return_period:
do_lec_output_agg_summary(
1, sidx, loss, filtered_occ_map,
outloss_mean_s, row_used_mean_s,
outloss_sample_s, row_used_sample_s,
num_sidxs, 1,
)
_MERGED_IDX_DTYPE = np.dtype([
("summary_id", np.int32),
("file_index", np.int32),
("offset", np.int64),
])
[docs]
def build_merged_idx(idx_handles):
"""Merge per-file .idx memmaps into one array sorted by summary_id."""
parts = []
for file_index, idx in enumerate(idx_handles):
if idx is None or len(idx) == 0:
continue
chunk = np.empty(len(idx), dtype=_MERGED_IDX_DTYPE)
chunk["summary_id"] = idx["summary_id"]
chunk["file_index"] = file_index
chunk["offset"] = idx["offset"]
parts.append(chunk)
if not parts:
return np.empty(0, dtype=_MERGED_IDX_DTYPE)
merged = np.concatenate(parts)
merged.sort(order="summary_id")
return merged
@nb.njit(cache=True, error_model="numpy")
@nb.njit(cache=True, error_model="numpy")
[docs]
def run_lec(
file_handles,
outloss_mean,
row_used_mean,
outloss_sample,
row_used_sample,
occ_csr,
use_return_period,
num_sidxs,
max_summary_id,
):
"""Process each summary file and populate outloss_mean and outloss_sample
Args:
file_handles (List[np.memmap]): List of memmaps for summary files data
outloss_mean (ndarray[OUTLOSS_DTYPE]): ndarray indexed by summary_id, period_no containing aggregate and max losses
row_used_mean (ndarray[bool]): bool mask for outloss_mean
outloss_sample (ndarray[OUTLOSS_DTYPE]): ndarray indexed by summary_id, sidx, period_no containing aggregate and max losses
row_used_sample (ndarray[bool]): bool mask for outloss_sample
occ_csr (OccurrenceCSR): id_index-backed CSR occurrence map
use_return_period (bool): Use Return Period file.
num_sidxs (int): Number of sidxs to consider for outloss_sample
max_summary_id (int): Max summary ID
"""
for fin in file_handles:
process_input_file(
fin,
outloss_mean,
row_used_mean,
outloss_sample,
row_used_sample,
occ_csr,
use_return_period,
num_sidxs,
max_summary_id,
)
def _open_output_files(outmap, stack, output_binary, output_parquet, noheader):
if output_binary:
for out_type in outmap:
if not outmap[out_type]["compute"]:
continue
outmap[out_type]["file"] = stack.enter_context(open(outmap[out_type]["file_path"], 'wb'))
elif output_parquet:
for out_type in outmap:
if not outmap[out_type]["compute"]:
continue
dtype = outmap[out_type]["dtype"]
schema = pa.schema([(name, pa.from_numpy_dtype(dtype[name])) for name in dtype.names])
outmap[out_type]["schema"] = schema
outmap[out_type]["file"] = stack.enter_context(pq.ParquetWriter(outmap[out_type]["file_path"], schema))
else:
for out_type in outmap:
if not outmap[out_type]["compute"]:
continue
out_file = stack.enter_context(open(outmap[out_type]["file_path"], 'w'))
if not noheader:
out_file.write(','.join(outmap[out_type]["headers"]) + '\n')
outmap[out_type]["file"] = out_file
[docs]
def run(
run_dir,
subfolder,
ept_output_file=None,
psept_output_file=None,
agg_full_uncertainty=False,
agg_wheatsheaf=False,
agg_sample_mean=False,
agg_wheatsheaf_mean=False,
occ_full_uncertainty=False,
occ_wheatsheaf=False,
occ_sample_mean=False,
occ_wheatsheaf_mean=False,
use_return_period=False,
noheader=False,
output_format="csv",
):
"""Runs LEC calculations
Args:
run_dir (str | os.PathLike): Path to directory containing required files structure
subfolder (str): Workspace subfolder inside <run_dir>/work/<subfolder>
ept_output_file (str, optional): Path to EPT output file. Defaults to None
psept_output_file (str, optional): Path to PSEPT output file. Defaults to None
agg_full_uncertainty (bool, optional): Aggregate Full Uncertainty. Defaults to False.
agg_wheatsheaf (bool, optional): Aggregate Wheatsheaf. Defaults to False.
agg_sample_mean (bool, optional): Aggregate Sample Mean. Defaults to False.
agg_wheatsheaf_mean (bool, optional): Aggregate Wheatsheaf Mean. Defaults to False.
occ_full_uncertainty (bool, optional): Occurrence Full Uncertainty. Defaults to False.
occ_wheatsheaf (bool, optional): Occurrence Wheatsheaf. Defaults to False.
occ_sample_mean (bool, optional): Occurrence Sample Mean. Defaults to False.
occ_wheatsheaf_mean (bool, optional): Occurrence Wheatsheaf Mean. Defaults to False.
use_return_period (bool, optional): Use Return Period file. Defaults to False.
noheader (bool): Boolean value to skip header in output file
output_format (str): Output format extension. Defaults to "csv".
"""
outmap = {
"ept": {
"compute": ept_output_file is not None,
"file_path": ept_output_file,
"fmt": EPT_fmt,
"headers": EPT_headers,
"file": None,
"dtype": EPT_dtype,
},
"psept": {
"compute": psept_output_file is not None,
"file_path": psept_output_file,
"fmt": PSEPT_fmt,
"headers": PSEPT_headers,
"file": None,
"dtype": PSEPT_dtype,
},
}
output_format = "." + output_format
output_binary = output_format == ".bin"
output_parquet = output_format == ".parquet"
# Check for correct suffix
for path in [v["file_path"] for v in outmap.values()]:
if path is None:
continue
if Path(path).suffix == "": # Ignore suffix for pipes
continue
if (Path(path).suffix != output_format):
raise ValueError(f"Invalid file extension for {output_format}, got {path},")
if not all([v["compute"] for v in outmap.values()]):
logger.warning("No output files specified")
with ExitStack() as stack:
workspace_folder = Path(run_dir, "work", subfolder)
if not workspace_folder.is_dir():
raise RuntimeError(f"Error: Unable to open directory {workspace_folder}")
# work folder for lec files
lec_files_folder = Path(workspace_folder, "lec_files")
lec_files_folder.mkdir(parents=False, exist_ok=True)
# Find summary binary files (sorted for stable pairing with .idx files)
files = sorted(workspace_folder.glob("*.bin"))
file_handles = [np.memmap(file, mode="r", dtype="u1") for file in files]
streams_in, (stream_source_type, stream_agg_type, sample_size) = init_streams_in(files, stack)
if stream_source_type != SUMMARY_STREAM_ID:
raise RuntimeError(f"Error: Not a summary stream type {stream_source_type}")
# Detect .idx files for per-summary streaming path
idx_handles_raw = []
for f in files:
idx_path = f.with_suffix(".idx")
if idx_path.exists():
if idx_path.stat().st_size > 0:
idx_handles_raw.append(np.memmap(str(idx_path), mode="r", dtype=summary_stream_index_dtype))
else:
idx_handles_raw.append(np.empty(0, dtype=summary_stream_index_dtype))
else:
idx_handles_raw.append(None)
n_idx = sum(h is not None for h in idx_handles_raw)
# All files must have .idx: the merged index covers all files simultaneously,
# so a gap would silently drop events from the missing file.
use_idx_path = n_idx == len(files) and n_idx > 0
if use_idx_path:
logger.info("Found %d/%d .idx file(s) — using per-summary-id indexing (reduced disk)", n_idx, len(files))
merged = build_merged_idx(idx_handles_raw)
unique_summary_ids = np.unique(merged["summary_id"])
max_summary_id = int(unique_summary_ids[-1]) if len(unique_summary_ids) > 0 else 0
else:
if 0 < n_idx < len(files):
logger.warning(
"%d/%d .idx files found — all .bin files must have .idx to use index path; "
"falling back to sequential scan", n_idx, len(files)
)
max_summary_id = get_max_summary_id(file_handles)
if max_summary_id == 0:
return
file_data, use_return_period, agg_wheatsheaf_mean, occ_wheatsheaf_mean = read_input_files(
run_dir,
use_return_period,
agg_wheatsheaf_mean,
occ_wheatsheaf_mean,
sample_size,
)
output_flags = [
agg_full_uncertainty,
agg_wheatsheaf,
agg_sample_mean,
agg_wheatsheaf_mean,
occ_full_uncertainty,
occ_wheatsheaf,
occ_sample_mean,
occ_wheatsheaf_mean,
]
# Check output_flags against output files
handles_agg = [AGG_FULL_UNCERTAINTY, AGG_SAMPLE_MEAN, AGG_WHEATSHEAF_MEAN]
handles_occ = [OCC_FULL_UNCERTAINTY, OCC_SAMPLE_MEAN, OCC_WHEATSHEAF_MEAN]
handles_psept = [AGG_WHEATSHEAF, OCC_WHEATSHEAF]
hasAGG = any([output_flags[idx] for idx in handles_agg])
hasOCC = any([output_flags[idx] for idx in handles_occ])
hasEPT = hasAGG or hasOCC
hasPSEPT = any([output_flags[idx] for idx in handles_psept])
outmap["ept"]["compute"] = outmap["ept"]["compute"] and hasEPT
outmap["psept"]["compute"] = outmap["psept"]["compute"] and hasPSEPT
if not outmap["ept"]["compute"]:
logger.warning("WARNING: no valid output stream to fill EPT file")
if not outmap["psept"]["compute"]:
logger.warning("WARNING: no valid output stream to fill PSEPT file")
if not (outmap["ept"]["compute"] or outmap["psept"]["compute"]):
return
# ── Per-summary streaming path (idx files present) ───────────────────────
if use_idx_path:
num_sidxs = int(sample_size) + 2
no_of_periods = int(file_data["no_of_periods"])
# Open output files before the loop so headers are written once
_open_output_files(outmap, stack, output_binary, output_parquet, noheader)
# Output buffers allocated once and reused across every per-summary generator
# call (the write_* generators overwrite each row before yielding, so no zeroing).
ept_buffer = np.empty(DEFAULT_BUFFER_SIZE, dtype=EPT_dtype)
psept_buffer = np.empty(DEFAULT_BUFFER_SIZE, dtype=PSEPT_dtype)
idx_config = LecConfig(
period_weights=file_data["period_weights"],
max_summary_id=1,
sample_size=int(sample_size),
no_of_periods=no_of_periods,
num_sidxs=num_sidxs,
use_return_period=use_return_period,
returnperiods=file_data["returnperiods"],
ept_buffer=ept_buffer,
psept_buffer=psept_buffer,
)
output_fn = make_output_fn(outmap, output_binary, output_parquet)
# Per-summary arrays — no × max_summary_id factor
outloss_mean_s = np.zeros(no_of_periods, dtype=OUTLOSS_DTYPE)
outloss_sample_s = np.zeros(no_of_periods * num_sidxs, dtype=OUTLOSS_DTYPE)
row_used_mean_s = np.zeros(no_of_periods, dtype=np.bool_)
row_used_sample_s = np.zeros(no_of_periods * num_sidxs, dtype=np.bool_)
# Pre-compute uint8 views once — numpy's struct[:]=0 is ~4x slower than raw byte zeroing
_mean_bytes = outloss_mean_s.view(np.uint8)
_sample_bytes = outloss_sample_s.view(np.uint8)
for summary_id in unique_summary_ids:
_mean_bytes[:] = 0
_sample_bytes[:] = 0
row_used_mean_s[:] = False
row_used_sample_s[:] = False
lo = int(np.searchsorted(merged["summary_id"], summary_id, side="left"))
hi = int(np.searchsorted(merged["summary_id"], summary_id, side="right"))
entries = merged[lo:hi]
for fi in np.unique(entries["file_index"]):
offsets = np.sort(entries[entries["file_index"] == fi]["offset"])
process_summary_entries(
file_handles[int(fi)],
offsets,
file_data["occ_csr"],
use_return_period,
outloss_mean_s,
row_used_mean_s,
outloss_sample_s,
row_used_sample_s,
num_sidxs,
)
output_for_summary_idx(
int(summary_id),
outloss_mean_s, row_used_mean_s,
outloss_sample_s, row_used_sample_s,
output_flags, hasOCC, hasAGG,
outmap=outmap, config=idx_config, output_fn=output_fn,
)
return # skip sequential path below
# ── Sequential path (no .idx files) ──────────────────────────────────────
# Check required disk space for work bdat files
num_sidxs = int(sample_size) + 2
no_of_periods = int(file_data["no_of_periods"])
_mean_elems = no_of_periods * max_summary_id
_sample_elems = no_of_periods * num_sidxs * max_summary_id
_required_bytes = (
_mean_elems * OUTLOSS_DTYPE.itemsize
+ _sample_elems * OUTLOSS_DTYPE.itemsize
+ _mean_elems # row_used_mean (bool_)
+ _sample_elems # row_used_sample (bool_)
)
_free_bytes = shutil.disk_usage(lec_files_folder).free
if _required_bytes > _free_bytes:
raise RuntimeError(
f"Insufficient disk space for lec .bdat files: "
f"{_required_bytes / 2**30:.2f} GiB required "
f"({no_of_periods} periods × {num_sidxs} sidxs × {max_summary_id} summary_ids), "
f"but only {_free_bytes / 2**30:.2f} GiB free at {lec_files_folder}. "
f"Run summarypy with -m to generate .idx files and avoid this pre-allocation."
)
# outloss_mean has only -1 SIDX
outloss_mean_file = Path(lec_files_folder, "lec_outloss_mean.bdat")
outloss_mean = np.memmap(
outloss_mean_file,
dtype=OUTLOSS_DTYPE,
mode="w+",
shape=(_mean_elems),
)
# outloss_sample has all SIDXs plus -2 and -3
outloss_sample_file = Path(lec_files_folder, "lec_outloss_sample.bdat")
outloss_sample = np.memmap(
outloss_sample_file,
dtype=OUTLOSS_DTYPE,
mode="w+",
shape=(_sample_elems),
)
row_used_mean_file = Path(lec_files_folder, "lec_row_used_mean.bdat")
row_used_mean = np.memmap(row_used_mean_file, dtype=np.bool_, mode="w+", shape=(len(outloss_mean),))
row_used_sample_file = Path(lec_files_folder, "lec_row_used_sample.bdat")
row_used_sample = np.memmap(row_used_sample_file, dtype=np.bool_, mode="w+", shape=(len(outloss_sample),))
# Run LEC calculations to populate outloss arrays
run_lec(
file_handles,
outloss_mean,
row_used_mean,
outloss_sample,
row_used_sample,
file_data["occ_csr"],
use_return_period,
num_sidxs,
max_summary_id,
)
# Initialise output files LEC
_open_output_files(outmap, stack, output_binary, output_parquet, noheader)
# Output aggregate reports to CSVs
ept_buffer = np.empty(DEFAULT_BUFFER_SIZE, dtype=EPT_dtype)
psept_buffer = np.empty(DEFAULT_BUFFER_SIZE, dtype=PSEPT_dtype)
seq_config = LecConfig(
period_weights=file_data["period_weights"],
max_summary_id=max_summary_id,
sample_size=int(sample_size),
no_of_periods=int(file_data["no_of_periods"]),
num_sidxs=num_sidxs,
use_return_period=use_return_period,
returnperiods=file_data["returnperiods"],
ept_buffer=ept_buffer,
psept_buffer=psept_buffer,
)
agg = AggReports(
outmap,
outloss_mean, row_used_mean,
outloss_sample, row_used_sample,
seq_config,
lec_files_folder,
output_binary,
output_parquet,
)
# Output Mean Damage Ratio
if outmap["ept"]["compute"]:
if hasOCC:
agg.output_mean_damage_ratio(OEP, OEPTVAR, "max_out_loss")
if hasAGG:
agg.output_mean_damage_ratio(AEP, AEPTVAR, "agg_out_loss")
# Output Full Uncertainty
if output_flags[OCC_FULL_UNCERTAINTY]:
agg.output_full_uncertainty(OEP, OEPTVAR, "max_out_loss")
if output_flags[AGG_FULL_UNCERTAINTY]:
agg.output_full_uncertainty(AEP, AEPTVAR, "agg_out_loss")
# Output Wheatsheaf and Wheatsheaf Mean
if output_flags[OCC_WHEATSHEAF] or output_flags[OCC_WHEATSHEAF_MEAN]:
agg.output_wheatsheaf_and_wheatsheafmean(
OEP, OEPTVAR, "max_out_loss",
output_flags[OCC_WHEATSHEAF], output_flags[OCC_WHEATSHEAF_MEAN]
)
if output_flags[AGG_WHEATSHEAF] or output_flags[AGG_WHEATSHEAF_MEAN]:
agg.output_wheatsheaf_and_wheatsheafmean(
AEP, AEPTVAR, "agg_out_loss",
output_flags[AGG_WHEATSHEAF], output_flags[AGG_WHEATSHEAF_MEAN]
)
# Output Sample Mean
if output_flags[OCC_SAMPLE_MEAN]:
agg.output_sample_mean(OEP, OEPTVAR, "max_out_loss")
if output_flags[AGG_SAMPLE_MEAN]:
agg.output_sample_mean(AEP, AEPTVAR, "agg_out_loss")
@redirect_logging(exec_name='lecpy')
[docs]
def main(
run_dir='.',
subfolder=None,
ept=None,
psept=None,
agg_full_uncertainty=False,
agg_wheatsheaf=False,
agg_sample_mean=False,
agg_wheatsheaf_mean=False,
occ_full_uncertainty=False,
occ_wheatsheaf=False,
occ_sample_mean=False,
occ_wheatsheaf_mean=False,
use_return_period=False,
noheader=False,
ext="csv",
**kwargs
):
run(
run_dir,
subfolder,
ept_output_file=ept,
psept_output_file=psept,
agg_full_uncertainty=agg_full_uncertainty,
agg_wheatsheaf=agg_wheatsheaf,
agg_sample_mean=agg_sample_mean,
agg_wheatsheaf_mean=agg_wheatsheaf_mean,
occ_full_uncertainty=occ_full_uncertainty,
occ_wheatsheaf=occ_wheatsheaf,
occ_sample_mean=occ_sample_mean,
occ_wheatsheaf_mean=occ_wheatsheaf_mean,
use_return_period=use_return_period,
noheader=noheader,
output_format=ext,
)