From 98fa32c24bdbefa0d31f6deb576e509a94f601ed Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Sierro?= Date: Tue, 8 Aug 2023 10:59:08 +0200 Subject: [PATCH] rudimentary new io data structure --- src/scgenerator/io.py | 212 +++++++++++- src/scgenerator/parameter.py | 64 ++-- src/scgenerator/physics/fiber.py | 14 +- src/scgenerator/physics/pulse.py | 7 +- src/scgenerator/spectra.py | 400 ++++++---------------- tests/data/PM2000D_2 extrapolated 4 0.npz | Bin 0 -> 32558 bytes tests/data/PM2000D_A_eff_marcuse.npz | Bin 0 -> 2916 bytes tests/test_io_handlers.py | 124 +++++++ 8 files changed, 475 insertions(+), 346 deletions(-) create mode 100644 tests/data/PM2000D_2 extrapolated 4 0.npz create mode 100644 tests/data/PM2000D_A_eff_marcuse.npz create mode 100644 tests/test_io_handlers.py diff --git a/src/scgenerator/io.py b/src/scgenerator/io.py index 53ef5df..406062c 100644 --- a/src/scgenerator/io.py +++ b/src/scgenerator/io.py @@ -1,8 +1,14 @@ +from __future__ import annotations + import datetime import json +import os +from dataclasses import dataclass from pathlib import Path -from typing import Sequence +from typing import BinaryIO, Protocol, Sequence +from zipfile import ZipFile +import numpy as np import pkg_resources @@ -11,27 +17,207 @@ def data_file(path: str) -> Path: return Path(pkg_resources.resource_filename("scgenerator", path)) -class DatetimeEncoder(json.JSONEncoder): +class CustomEncoder(json.JSONEncoder): def default(self, obj): if isinstance(obj, (datetime.date, datetime.datetime)): return obj.isoformat() + elif isinstance(obj, np.ndarray): + return tuple(obj) + elif isinstance(obj, DataFile): + return obj.__json__() -def decode_datetime_hook(obj): +def _decode_datetime(s: str) -> datetime.date | datetime.datetime | str: + try: + return datetime.datetime.fromisoformat(s) + except Exception: + pass + try: + return datetime.date.fromisoformat(s) + except Exception: + pass + return s + + +def custom_decode_hook(obj): + """ + Some simple, non json-compatible objects are encoded as str, this function reconstructs the + original object and sets it in the decoded json structure + """ for k, v in obj.items(): - if not isinstance(v, str): - continue - try: - dt = datetime.datetime.fromisoformat(v) - except Exception: - try: - dt = datetime.date.fromisoformat(v) - except Exception: - continue - obj[k] = dt + if isinstance(v, str): + obj[k] = _decode_datetime(v) + elif isinstance(v, list): + obj[k] = tuple(v) return obj +class PropagationIOHandler(Protocol): + def __len__(self) -> int: + ... + + def keys(self) -> list[str]: + ... + + def save_spectrum(self, index: int, spectrum: np.ndarray): + ... + + def load_spectrum(self, index: int) -> np.ndarray: + ... + + def save_data(self, name: str, data: bytes): + ... + + def load_data(self, name: str) -> bytes: + ... + + +class MemoryIOHandler: + spectra: dict[int, np.ndarray] + data: dict[str, bytes] + + def __init__(self): + self.spectra = {} + self.data = {} + + def __len__(self) -> int: + return len(self.spectra) + + def keys(self) -> list[str]: + return list(self.data.keys()) + + def save_spectrum(self, index: int, spectrum: np.ndarray): + self.spectra[index] = spectrum + + def load_spectrum(self, index: int) -> np.ndarray: + return self.spectra[index] + + def save_data(self, name: str, data: bytes): + self.data[name] = data + + def load_data(self, name: str) -> bytes: + return self.data[name] + + +class ZipFileIOHandler: + file: BinaryIO + SPECTRUM_FN = "spectra/spectrum_{}.npy" + + def __init__(self, file: os.PathLike): + """ + Create a IO handler to be used in Propagation. + This handler saves spectra as numpy `.npy` files. + + Parameters + ---------- + file : os.PathLike + path to the desired zip file. Will be created if it doesn't exist. + Passing in an already opened file-like object is not supported at the moment + """ + self.file = Path(file) + if not self.file.exists(): + ZipFile(self.file, "w").close() + + def __len__(self) -> int: + with ZipFile(self.file, "r") as zip_file: + i = 0 + while True: + try: + zip_file.open(self.SPECTRUM_FN.format(i)) + except KeyError: + return i + i += 1 + + def keys(self) -> list[str]: + with ZipFile(self.file, "r") as zip_file: + return zip_file.namelist() + + def save_spectrum(self, index: int, spectrum: np.ndarray): + with ZipFile(self.file, "a") as zip_file, zip_file.open( + self.SPECTRUM_FN.format(index), "w" + ) as file: + np.lib.format.write_array(file, spectrum, allow_pickle=False) + + def load_spectrum(self, index: int) -> np.ndarray: + with ZipFile(self.file, "r") as zip_file, zip_file.open( + self.SPECTRUM_FN.format(index), "r" + ) as file: + return np.load(file) + + def save_data(self, name: str, data: bytes): + with ZipFile(self.file, "a") as zip_file, zip_file.open(name, "w") as file: + file.write(data) + + def load_data(self, name: str) -> bytes: + with ZipFile(self.file, "r") as zip_file, zip_file.open(name, "r") as file: + return file.read() + + +@dataclass +class DataFile: + """ + Holds information about external files necessary for a simulation. + In the current implementation, only reading data from the file is supported + """ + + prefix: str | None + path: str + io: PropagationIOHandler + + PREFIX_SEP = "::" + + @classmethod + def from_str(cls, s: str, io: PropagationIOHandler | None = None) -> DataFile: + if cls.PREFIX_SEP in s: + prefix, path = s.split(cls.PREFIX_SEP, 1) + else: + prefix = None + path = s + return cls(prefix, path, io) + + @classmethod + def validate(cls, name: str, data: str | DataFile) -> DataFile: + """To be used to automatically construct a DataFile when creating a Parameters obj""" + if isinstance(data, cls): + return data + elif isinstance(data, str): + return cls.from_str(data) + elif isinstance(data, Path): + return cls.from_str(os.fspath(data)) + else: + raise TypeError( + f"{name!r} must be a path or a bundled file specifier, not a {type(data)!r}" + ) + + def __json__(self) -> str: + if self.prefix is None: + return os.fspath(Path(self.path)) + return self.prefix + self.PREFIX_SEP + self.path + + def load_data(self) -> bytes: + if self.prefix is not None and self.io is None: + raise ValueError( + f"a bundled file prefixed with {self.prefix} " + "must have a PropagationIOHandler attached" + ) + if self.io is not None: + return self.io.load_data(self.path) + else: + return Path(self.path).read_bytes() + + +def unique_name(base_name: str, existing: set[str]) -> str: + name = base_name + p = Path(base_name) + base = p.stem + ext = p.suffix + i = 0 + while name in existing: + name = f"{base}_{i}{ext}" + i += 1 + return name + + def format_graph(left_elements: Sequence[str], middle: str, right_elements: Sequence[str]): if len(left_elements) == 0: left_elements = [""] diff --git a/src/scgenerator/parameter.py b/src/scgenerator/parameter.py index 71318ba..7d36d04 100644 --- a/src/scgenerator/parameter.py +++ b/src/scgenerator/parameter.py @@ -4,7 +4,7 @@ import datetime as datetime_module import json import os from copy import copy, deepcopy -from dataclasses import dataclass, field, fields +from dataclasses import dataclass, field from functools import lru_cache, wraps from math import isnan from pathlib import Path @@ -12,12 +12,10 @@ from typing import Any, Callable, ClassVar, Iterable, Iterator, Set, Type, TypeV import numpy as np -from scgenerator import utils from scgenerator.const import MANDATORY_PARAMETERS, __version__ from scgenerator.evaluator import Evaluator, EvaluatorError -from scgenerator.io import DatetimeEncoder, decode_datetime_hook +from scgenerator.io import CustomEncoder, DataFile, custom_decode_hook from scgenerator.operators import Qualifier, SpecOperator -from scgenerator.utils import update_path_name T = TypeVar("T") DISPLAY_INFO = {} @@ -77,6 +75,11 @@ def string(name, n): raise ValueError(f"{name!r} must not be empty") +def low_string(name, n): + string(name, n) + return n.lower() + + def in_range_excl(_min, _max): @type_checker(float, int) def _in_range(name, n): @@ -255,9 +258,6 @@ class Parameter: del instance._param_dico[self.name] def __set__(self, instance: Parameters, value): - if instance._frozen: - raise AttributeError("Parameters instance is frozen and can no longer be modified") - if isinstance(value, Parameter): if self.default is not None: instance._param_dico[self.name] = copy(self.default) @@ -288,9 +288,9 @@ class Parameter: except TypeError: is_value = True if is_value: - if self.converter is not None: - v = self.converter(v) - self._validator(self.name, v) + ret_val = self._validator(self.name, v) + if ret_val is not None: + v = ret_val return is_value, v @@ -307,7 +307,6 @@ class Parameters: # root name: str = Parameter(string, default="no name") - output_path: Path = Parameter(type_checker(Path), converter=Path) # fiber input_transmission: float = Parameter(in_range_incl(0, 1), default=1.0) @@ -315,10 +314,10 @@ class Parameters: n2: float = Parameter(non_negative(float, int)) chi3: float = Parameter(non_negative(float, int)) loss: str = Parameter(literal("capillary")) - loss_file: str = Parameter(string) + loss_file: DataFile = Parameter(DataFile.validate) effective_mode_diameter: float = Parameter(positive(float, int)) effective_area: float = Parameter(non_negative(float, int)) - effective_area_file: str = Parameter(string) + effective_area_file: DataFile = Parameter(DataFile.validate) numerical_aperture: float = Parameter(in_range_excl(0, 1)) pcf_pitch: float = Parameter(in_range_excl(0, 1e-3), display_info=(1e6, "μm")) pcf_pitch_ratio: float = Parameter(in_range_excl(0, 1)) @@ -326,7 +325,7 @@ class Parameters: he_mode: tuple[int, int] = Parameter(int_pair, default=(1, 1)) fit_parameters: tuple[int, int] = Parameter(float_pair, default=(0.08, 200e-9)) beta2_coefficients: Iterable[float] = Parameter(num_list) - dispersion_file: str = Parameter(string) + dispersion_file: DataFile = Parameter(DataFile.validate) model: str = Parameter( literal("pcf", "marcatili", "marcatili_adjusted", "hasan", "custom"), ) @@ -345,7 +344,7 @@ class Parameters: capillary_nested: int = Parameter(non_negative(int), default=0) # gas - gas_name: str = Parameter(string, converter=str.lower, default="vacuum") + gas_name: str = Parameter(low_string, default="vacuum") pressure: float = Parameter(non_negative(float, int), display_info=(1e-5, "bar")) pressure_in: float = Parameter(non_negative(float, int), display_info=(1e-5, "bar")) pressure_out: float = Parameter(non_negative(float, int), display_info=(1e-5, "bar")) @@ -353,7 +352,7 @@ class Parameters: plasma_density: float = Parameter(non_negative(float, int), default=0) # pulse - field_file: str = Parameter(string) + field_file: DataFile = Parameter(DataFile.validate) input_time: np.ndarray = Parameter(type_checker(np.ndarray)) input_field: np.ndarray = Parameter(type_checker(np.ndarray)) repetition_rate: float = Parameter( @@ -380,11 +379,9 @@ class Parameters: # simulation full_field: bool = Parameter(boolean, default=False) integration_scheme: str = Parameter( - literal("erk43", "erk54", "cqe", "sd", "constant"), - converter=str.lower, - default="erk43", + literal("erk43", "erk54", "cqe", "sd", "constant"), default="erk43" ) - raman_type: str = Parameter(literal("measured", "agrawal", "stolen"), converter=str.lower) + raman_type: str = Parameter(literal("measured", "agrawal", "stolen")) raman_fraction: float = Parameter(non_negative(float, int)) spm: bool = Parameter(boolean, default=True) repeat: int = Parameter(positive(int), default=1) @@ -437,7 +434,8 @@ class Parameters: @classmethod def from_json(cls, s: str) -> Parameters: - return cls(**json.loads(s, object_hook=decode_datetime_hook)) + decoded = json.loads(s, object_hook=custom_decode_hook) + return cls(**decoded) @classmethod def load(cls, path: os.PathLike) -> Parameters: @@ -462,18 +460,32 @@ class Parameters: self.__post_init__() def __setattr__(self, k, v): - if self._frozen: + if self._frozen and not k.endswith("_file"): raise AttributeError( f"cannot set attribute to frozen {self.__class__.__name__} instance" ) object.__setattr__(self, k, v) - def copy(self) -> Parameters: - return Parameters(**deepcopy(self.strip_params_dict())) + def items(self) -> Iterator[tuple[str, Any]]: + for k, v in self._param_dico.items(): + if v is None: + continue + yield k, v + + def copy(self, freeze: bool = False) -> Parameters: + """create a deep copy of self. if freeze is True, the returned copy is read-only""" + params = Parameters(**deepcopy(self.strip_params_dict())) + if freeze: + params.freeze() + return params + + def freeze(self): + """render the current instance read-only. This is not reversible""" + self._frozen = True def to_json(self) -> str: d = self.dump_dict() - return json.dumps(d, cls=DatetimeEncoder, default=list) + return json.dumps(d, cls=CustomEncoder, indent=4) def get_evaluator(self): evaluator = Evaluator.default(self.full_field) @@ -611,7 +623,7 @@ class Parameters: "linear_op", "c_to_a_factor", } - types = (np.ndarray, float, int, str, list, tuple, Path) + types = (np.ndarray, float, int, str, list, tuple, Path, DataFile) c = deepcopy if copy else lambda x: x out = {} for key, value in self._param_dico.items(): diff --git a/src/scgenerator/physics/fiber.py b/src/scgenerator/physics/fiber.py index 45a3772..e7aeb8d 100644 --- a/src/scgenerator/physics/fiber.py +++ b/src/scgenerator/physics/fiber.py @@ -1,3 +1,4 @@ +from io import BytesIO from typing import Iterable, TypeVar import numpy as np @@ -7,6 +8,7 @@ from numpy.polynomial.chebyshev import Chebyshev, cheb2poly from scipy.interpolate import interp1d from scgenerator import io +from scgenerator.io import DataFile from scgenerator.math import argclosest, u_nm from scgenerator.physics import materials as mat from scgenerator.physics import units @@ -653,7 +655,7 @@ def saitoh_paramters(pcf_pitch_ratio: float) -> tuple[float, float]: return A, B -def load_custom_effective_area(effective_area_file: str, l: np.ndarray) -> np.ndarray: +def load_custom_effective_area(effective_area_file: DataFile, l: np.ndarray) -> np.ndarray: """ loads custom effective area file @@ -669,14 +671,14 @@ def load_custom_effective_area(effective_area_file: str, l: np.ndarray) -> np.nd np.ndarray, shape (n,) wl-dependent effective mode field area """ - data = np.load(effective_area_file) + data = np.load(BytesIO(effective_area_file.load_data())) effective_area = data.get("A_eff", data.get("effective_area")) wl = data["wavelength"] return interp1d(wl, effective_area, fill_value=1, bounds_error=False)(l) -def load_custom_dispersion(dispersion_file: str) -> tuple[np.ndarray, np.ndarray]: - disp_file = np.load(dispersion_file) +def load_custom_dispersion(dispersion_file: DataFile) -> tuple[np.ndarray, np.ndarray]: + disp_file = np.load(BytesIO(dispersion_file.load_data())) wl_for_disp = disp_file["wavelength"] interp_range = (np.min(wl_for_disp), np.max(wl_for_disp)) D = disp_file["dispersion"] @@ -684,7 +686,7 @@ def load_custom_dispersion(dispersion_file: str) -> tuple[np.ndarray, np.ndarray return wl_for_disp, beta2, interp_range -def load_custom_loss(l: np.ndarray, loss_file: str) -> np.ndarray: +def load_custom_loss(l: np.ndarray, loss_file: DataFile) -> np.ndarray: """ loads a npz loss file that contains a wavelength and a loss entry @@ -700,7 +702,7 @@ def load_custom_loss(l: np.ndarray, loss_file: str) -> np.ndarray: np.ndarray, shape (n,) loss in 1/m units """ - loss_data = np.load(loss_file) + loss_data = np.load(BytesIO(loss_file.load_data())) wl = loss_data["wavelength"] loss = loss_data["loss"] return interp1d(wl, loss, fill_value=0, bounds_error=False)(l) diff --git a/src/scgenerator/physics/pulse.py b/src/scgenerator/physics/pulse.py index 07c785e..4ece788 100644 --- a/src/scgenerator/physics/pulse.py +++ b/src/scgenerator/physics/pulse.py @@ -12,6 +12,7 @@ n is the number of spectra at the same z position and nt is the size of the time import itertools import os from dataclasses import astuple, dataclass +from io import BytesIO from pathlib import Path from typing import Literal, Tuple, TypeVar @@ -25,6 +26,7 @@ from scipy.optimize._optimize import OptimizeResult from scgenerator import math from scgenerator.defaults import default_plotting +from scgenerator.io import DataFile from scgenerator.physics import units c = 299792458.0 @@ -410,8 +412,9 @@ def interp_custom_field( return field_0 -def load_custom_field(field_file: str) -> tuple[np.ndarray, np.ndarray]: - field_data = np.load(field_file) +def load_custom_field(field_file: DataFile) -> tuple[np.ndarray, np.ndarray]: + data = field_file.load_data() + field_data = np.load(BytesIO(data)) return field_data["time"], field_data["field"] diff --git a/src/scgenerator/spectra.py b/src/scgenerator/spectra.py index acd9e51..d6219d6 100644 --- a/src/scgenerator/spectra.py +++ b/src/scgenerator/spectra.py @@ -1,37 +1,32 @@ from __future__ import annotations import os +import warnings from pathlib import Path -from typing import Callable, Iterator, Optional, Union -import matplotlib.pyplot as plt import numpy as np from scgenerator import math -from scgenerator.const import PARAM_FN, SPEC1_FN, SPEC1_FN_N -from scgenerator.logger import get_logger +from scgenerator.io import DataFile, PropagationIOHandler, ZipFileIOHandler, unique_name from scgenerator.parameter import Parameters from scgenerator.physics import pulse, units -from scgenerator.physics.units import PlotRange -from scgenerator.plotting import ( - mean_values_plot, - propagation_plot, - single_position_plot, - transform_1D_values, - transform_2D_propagation, -) -from scgenerator.utils import load_spectrum, load_toml, simulations_list class Spectrum(np.ndarray): - params: Parameters + w: np.ndarray + t: np.ndarray + l: np.ndarray + ifft: Callable[[np.ndarray], np.ndarray] def __new__(cls, input_array, params: Parameters): # Input array is an already formed ndarray instance # We first cast to be our class type obj = np.asarray(input_array).view(cls) # add the new attribute to the created instance - obj.params = params + obj.w = params.compute("w") + obj.t = params.compute("t") + obj.t = params.compute("t") + obj.ifft = params.compute("ifft") # Finally, we must return the newly created object: return obj @@ -40,14 +35,17 @@ class Spectrum(np.ndarray): # see InfoArray.__array_finalize__ for comments if obj is None: return - self.params = getattr(obj, "params", None) + self.w = getattr(obj, "w", None) + self.t = getattr(obj, "t", None) + self.l = getattr(obj, "l", None) + self.ifft = getattr(obj, "ifft", None) def __getitem__(self, key) -> "Spectrum": return super().__getitem__(key) @property def wl_int(self): - return units.to_WL(math.abs2(self), self.params.l) + return units.to_WL(math.abs2(self), self.l) @property def freq_int(self): @@ -59,13 +57,13 @@ class Spectrum(np.ndarray): @property def time_int(self): - return math.abs2(self.params.ifft(self)) + return math.abs2(self.ifft(self)) def amplitude(self, unit): if unit.type in ["WL", "FREQ", "AFREQ"]: - x_axis = unit.inv(self.params.w) + x_axis = unit.inv(self.w) else: - x_axis = unit.inv(self.params.t) + x_axis = unit.inv(self.t) order = np.argsort(x_axis) func = dict( @@ -84,7 +82,7 @@ class Spectrum(np.ndarray): np.sqrt( units.to_WL( math.abs2(self), - self.params.l, + self.l, ) ) * self @@ -106,311 +104,115 @@ class Spectrum(np.ndarray): @property def wl_max(self): if self.ndim == 1: - return self.params.l[np.argmax(self.wl_int, axis=-1)] + return self.l[np.argmax(self.wl_int, axis=-1)] return np.array([s.wl_max for s in self]) def mask_wl(self, pos: float, width: float) -> Spectrum: - return self * np.exp( - -(((self.params.l - pos) / (pulse.fwhm_to_T0_fac["gaussian"] * width)) ** 2) - ) + return self * np.exp(-(((self.l - pos) / (pulse.fwhm_to_T0_fac["gaussian"] * width)) ** 2)) def measure(self) -> tuple[float, float, float]: - return pulse.measure_field(self.params.t, self.time_amp) + return pulse.measure_field(self.t, self.time_amp) -class SimulationSeries: - """ - SimulationsSeries are the interface the user should use to load and - interact with simulation data. The object loads each fiber of the simulation - into a separate object and exposes convenience methods to make the series behave - as a single fiber. - - It should be noted that the last spectrum of a fiber and the first one of the next - fibers are identical. Therefore, SimulationSeries object will return fewer datapoints - than when manually mergin the corresponding data. - - """ - - path: Path - fibers: list[SimulatedFiber] +class Propagation: + io: PropagationIOHandler params: Parameters - z_indices: list[tuple[int, int]] - fiber_positions: list[tuple[str, float]] + _current_index: int - def __init__(self, path: os.PathLike): + PARAMS_FN = "params.json" + + def __init__( + self, + io_handler: PropagationIOHandler, + params: Parameters | None = None, + bundle_data: bool = False, + ): """ - Create a SimulationSeries + A propagation is the object that manages IO for one single propagation. + It is recommended to use the "memory_propagation" and "zip_propagation" convenience + factories instead of creating this object directly. Parameters ---------- - path : os.PathLike - path to the last fiber of the series - - Raises - ------ - FileNotFoundError - No simulation found in specified directory + io : PropagationIOHandler + object that implements the PropagationIOHandler Protocol. + params : Parameters + simulations parameters. Those will be saved via the """ - self.logger = get_logger() - for self.path in simulations_list(path): - break - else: - raise FileNotFoundError(f"No simulation in {path}") - self.fibers = [SimulatedFiber(self.path)] - while (p := self.fibers[-1].params.prev_data_dir) is not None: - p = Path(p) - if not p.is_absolute(): - p = Path(self.fibers[-1].params.output_path) / p - self.fibers.append(SimulatedFiber(p)) - self.fibers = self.fibers[::-1] + self.io = io_handler + self._current_index = len(self.io) - self.fiber_positions = [(self.fibers[0].params.name, 0.0)] - self.params = Parameters(**self.fibers[0].params.dump_dict(False, False)) - z_targets = list(self.params.z_targets) - self.z_indices = [(0, j) for j in range(self.params.z_num)] - for i, fiber in enumerate(self.fibers[1:]): - self.fiber_positions.append((fiber.params.name, z_targets[-1])) - z_targets += list(fiber.params.z_targets[1:] + z_targets[-1]) - self.z_indices += [(i + 1, j) for j in range(1, fiber.params.z_num)] - self.params.z_targets = np.array(z_targets) - self.params.length = self.params.z_targets[-1] - self.params.z_num = len(self.params.z_targets) + new_params = params is not None + if not new_params: + if bundle_data: + raise ValueError( + "cannot bundle data to existing Propagation. Create a new one instead" + ) + try: + params_data = self.io.load_data(self.PARAMS_FN) + params = Parameters.from_json(params_data.decode()) + except KeyError: + raise ValueError(f"Missing Parameters in {self.__class__.__name__}.") from None - def spectra( - self, z_descr: Union[float, int, None] = None, sim_ind: Optional[int] = 0 - ) -> Spectrum: - ... - if z_descr is None: - out = [self.fibers[i].spectra(j, sim_ind) for i, j in self.z_indices] - else: - if isinstance(z_descr, (float, np.floating)): - fib_ind, z_ind = self.z_ind(z_descr) - else: - fib_ind, z_ind = self.z_indices[z_descr] - out = self.fibers[fib_ind].spectra(z_ind, sim_ind) - return Spectrum(out, self.params) + self.params = params - def fields( - self, z_descr: Union[float, int, None] = None, sim_ind: Optional[int] = 0 - ) -> Spectrum: - return self.params.ifft(self.spectra(z_descr, sim_ind)) + if bundle_data: + if not params._frozen: + warnings.warn( + f"Parameters instance {params.name!r} is not frozen but will be saved " + "in an unmodifiable state anyway." + ) + _bundle_external_files(self.params, self.io) - def z_ind(self, pos: float) -> tuple[int, int]: - if 0 <= pos <= self.params.length: - ind = np.argmin(np.abs(self.params.z_targets - pos)) - return self.z_indices[ind] - else: - raise ValueError(f"cannot match z={pos} with max length of {self.params.length}") + if new_params: + self.io.save_data(self.PARAMS_FN, self.params.to_json().encode()) - def plot_2D( - self, - left: float, - right: float, - unit: Union[Callable[[float], float], str], - ax: plt.Axes, - sim_ind: int = 0, - **kwargs, - ): - plot_range = PlotRange(left, right, unit) - vals = self.retrieve_plot_values(plot_range, None, sim_ind) - return propagation_plot(vals, plot_range, self.params, ax, **kwargs) + def __len__(self) -> int: + return self._current_index - def plot_values_2D( - self, - left: float, - right: float, - unit: Union[Callable[[float], float], str], - sim_ind: int = 0, - **kwargs, - ): - plot_range = PlotRange(left, right, unit) - vals = self.retrieve_plot_values(plot_range, None, sim_ind) - return transform_2D_propagation(vals, plot_range, self.params, **kwargs) + def __getitem__(self, key: int | slice) -> Spectrum: + if isinstance(key, slice): + return self._load_slice(key) + if isinstance(key, (float, np.floating)): + key = math.argclosest(self.params.compute("z_targets"), key) + array = self.io.load_spectrum(key) + return Spectrum(array, self.params) - def plot_1D( - self, - left: float, - right: float, - unit: Union[Callable[[float], float], str], - ax: plt.Axes, - z_pos: int, - sim_ind: int = 0, - **kwargs, - ): - plot_range = PlotRange(left, right, unit) - vals = self.retrieve_plot_values(plot_range, z_pos, sim_ind) - return single_position_plot(vals, plot_range, self.params, ax, **kwargs) + def __setitem__(self, key: int, value: np.ndarray): + if not isinstance(key, int): + raise TypeError(f"Cannot save a spectrum at index {key!r} of type {type(key)!r}") + self.io.save_spectrum(key, np.asarray(value)) - def plot_values_1D( - self, - left: float, - right: float, - unit: Union[Callable[[float], float], str], - z_pos: int, - sim_ind: int = 0, - **kwargs, - ) -> tuple[np.ndarray, np.ndarray]: - """ - gives the desired values already tranformes according to the give range + def _load_slice(self, key: slice) -> Spectrum: + _iter = range(len(self))[key] + out = Spectrum(np.zeros((len(_iter), self.params.t_num), dtype=complex), self.params) + for i in _iter: + out[i] = self.io.load_spectrum(i) + return out - Parameters - ---------- - left : float - leftmost limit in unit - right : float - rightmost limit in unit - unit : Union[Callable[[float], float], str] - unit - z_pos : Union[int, float] - position either as an index (int) or a real position (float) - sim_ind : Optional[int] - which simulation to take when more than one are present + def append(self, spectrum: np.ndarray): + self.io.save_spectrum(self._current_index, spectrum.asarray()) + self._current_index += 1 - Returns - ------- - np.ndarray - x axis - np.ndarray - y values - """ - plot_range = PlotRange(left, right, unit) - vals = self.retrieve_plot_values(plot_range, z_pos, sim_ind) - return transform_1D_values(vals, plot_range, self.params, **kwargs) - - def plot_mean( - self, - left: float, - right: float, - unit: Union[Callable[[float], float], str], - ax: plt.Axes, - z_pos: int, - **kwargs, - ): - plot_range = PlotRange(left, right, unit) - vals = self.retrieve_plot_values(plot_range, z_pos, None) - return mean_values_plot(vals, plot_range, self.params, ax, **kwargs) - - def retrieve_plot_values( - self, plot_range: PlotRange, z_pos: Optional[Union[int, float]], sim_ind: Optional[int] - ): - if plot_range.unit.type == "TIME": - return self.fields(z_pos, sim_ind) - else: - return self.spectra(z_pos, sim_ind) - - def rin_propagation( - self, left: float, right: float, unit: str - ) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - """ - returns the RIN as function of unit and z - - Parameters - ---------- - left : float - left limit in unit - right : float - right limit in unit - unit : str - unit descriptor - - Returns - ------- - x : np.ndarray, shape (nt,) - x axis - y : np.ndarray, shape (z_num, ) - y axis - rin_prop : np.ndarray, shape (z_num, nt) - RIN - """ - spectra = [] - for spec in np.moveaxis(self.spectra(None, None), 1, 0): - x, z, tmp = transform_2D_propagation(spec, (left, right, unit), self.params, False) - spectra.append(tmp) - return x, z, pulse.rin_curve(np.moveaxis(spectra, 0, 1)) - - # Magic methods - - def __iter__(self) -> Iterator[Spectrum]: - for i, j in self.z_indices: - yield self.fibers[i].spectra(j, None) - - def __repr__(self) -> str: - return f"{self.__class__.__name__}(path={self.path})" - - def __eq__(self, other: SimulationSeries) -> bool: - return self.path == other.path and self.params == other.params - - def __contains__(self, fiber: SimulatedFiber) -> bool: - return fiber in self.fibers - - def __getitem__(self, key) -> Spectrum: - if isinstance(key, tuple): - return self.spectra(*key) - else: - return self.spectra(key, None) + def load_all(self) -> Spectrum: + return self._load_slice(slice(None)) -class SimulatedFiber: - params: Parameters - t: np.ndarray - w: np.ndarray +def _bundle_external_files(params: Parameters, io: PropagationIOHandler): + """copies every external file specified in the parameters and saves it""" + existing_files = set(io.keys()) + for _, value in params.items(): + if isinstance(value, DataFile): + data = value.load_data() - def __init__(self, path: os.PathLike): - self.path = Path(path) - self.params = Parameters(**load_toml(self.path / PARAM_FN)) - self.params.output_path = str(self.path.resolve()) - self.t = self.params.t - self.w = self.params.w - self.z = self.params.z_targets + value.io = io + value.prefix = "zip" + value.path = unique_name(Path(value.path).name, existing_files) + existing_files.add(value.path) - def spectra( - self, z_descr: Union[float, int, None] = None, sim_ind: Optional[int] = 0 - ) -> np.ndarray: - if z_descr is None: - out = [self.spectra(i, sim_ind) for i in range(self.params.z_num)] - else: - if isinstance(z_descr, (float, np.floating)): - return self.spectra(self.z_ind(z_descr), sim_ind) - else: - z_ind = z_descr + io.save_data(value.path, data) - if z_ind < 0: - z_ind = self.params.z_num + z_ind - if sim_ind is None: - out = [self._load_1(z_ind, i) for i in range(self.params.repeat)] - else: - out = self._load_1(z_ind) - return Spectrum(out, self.params) - - def z_ind(self, pos: float) -> int: - if 0 <= pos <= self.z[-1]: - return np.argmin(np.abs(self.z - pos)) - else: - raise ValueError(f"cannot match z={pos} with max length of {self.params.length}") - - def _load_1(self, z_ind: int, sim_ind=0) -> np.ndarray: - """ - loads a spectrum file - - Parameters - ---------- - z_ind : int - z_index relative to the entire simulation - sim_ind : int, optional - simulation index, used when repeated simulations with same parameters are ran, - by default 0 - - Returns - ------- - np.ndarray - loaded spectrum file - """ - if sim_ind > 0: - return load_spectrum(self.path / SPEC1_FN_N.format(z_ind, sim_ind)) - else: - return load_spectrum(self.path / SPEC1_FN.format(z_ind)) - psd = np.fft.rfft(signal) / np.sqrt(0.5 * len(time) / dt) - - def __repr__(self) -> str: - return f"{self.__class__.__name__}(path={self.path})" +def load_all(path: os.PathLike) -> Spectrum: + io = ZipFileIOHandler(path) + return Propagation(io).load_all() diff --git a/tests/data/PM2000D_2 extrapolated 4 0.npz b/tests/data/PM2000D_2 extrapolated 4 0.npz new file mode 100644 index 0000000000000000000000000000000000000000..1a9da5e0e5c8c10d6ae4aa5b924310765d818a75 GIT binary patch literal 32558 zcmd4YdsNKdA3yqvh?FAIjqXWNB1$}olv_j>Qly9|QKBRwMWkGcC=rn^Q%!Se?(?4e zOc4CvQu^ORQ%OfDHuz9zL}=u$n6O!q`wuHCY5n)i|DyjF@xPY2*4NX^Us>sp(h;j2 zp$E1{TP-lQTDa56%GTIw=f3Ee=-|kpebGBY|2OX*9C0A@-~55F;QgWh#?u|_?Hz2b zZH`ecowB9hgy!v}m*ErgMuSQku8u$)aNPXRIhV;ZY2Gqku89*B4_~>E3-U zn~QP7`S`Q*JBsm4e#vEVWHHA1-M!TnTa0Dz7^82TC`NpA&0f#6Vo;P`*AHeD<738G zXVQ&g7#peN2NV@!W2w$4bxJWZoUJ0*oMK4j&QEqri!tcC#!2r{G3LD6dPDJ|7<;0l zM#i-k9#v2vSJ2g}c_fpPv$1x>1>ikfpeL@Ln@|LhH!xD7* zet)>kvIM5@R2{$DmB2rGO73mv5?oNSQQoku1d5DByOcdkFl6XeOZF?lywZ?4+qRWp zuk*oc|`R0U19xcIycclA9r%SLhnwVXCp#;a3o@6KFlz@@(?vKr# z5_B2%ZGTl(0*lfyHD?(m2y`~IUmz?&rhL|wk2NKz@pT=%`m_Yf@4Q1+G?&0BdZ(iA zLkadN9iDZsrv!x=X<2^1OVDh1V?b546q8EH+vu95SmP{~?U+)Elk#V_I;N#y`*vIs z+mxcGZQxh5V=1ObYXq4tE=7=%vGnnZQe4S!usN}z6t#xSFWLu|Vr1#YzBi$zaB<$V z<@~-DMm*JoQwTZ3MfT> zv+OIy!HnaZ8k9&N8D{cR#*)xha@X8ZhXj3R=Gk{fB&?U$bZ1(Tkm~!|Z`o`Th_oN|Vm4UrWM{46Cn7n@GqtT;NL%CgDlR8ct9+31ggtOvWE1VTnBI zH2)Y0alXl4_NJ0h(RRh#=n@I6>IzBb8M-;i*lZP6!>k0d;c^788YM#4Cy?NxMfT4zn^C%L&@VnPu?veWEX(SB8wXcgLfbltB{J z_ugbx83u>Ocs}+n!<-C5>hUdQ*kfp?XTPfqw@O?S-|R2L3n!oUjQBE)m+xGSX}WpYWuQgfIJ)Lm89Il^ZNE#(V45LbSxPH|zhOhgW_}qilypoQt0;rq z>Gu)V<1+k}X|#sDDubi1@ro(!Wr%2VD3^DY;ZD?Yo!DPxcsaB&-cqF;6EgO+JX0^n z3d5uBDU-@^yyQG-?zD0+oNiBeH@zHPGFn`wLpdyb<;_bMmLssO(e2yva%4t*F1fX? z95q9KwbuKWLpeh`R%v@VoD9uhS45O!U&*{>L5Ipw=(Ms}^LRO4%L2yp&y+*QH}c@# z%yO)5JO0WruN+BHnXV5C%fTMHd+!Lj9NpY33D3amsGkKC1#;S@FV&P8Im>!@VtXt-y@7hicKQD-ay@Hqykm0#}E= zKdswZfx7fjF2}4FvNyS65(W$&^%dCR)INgMRDoPs z-|n#Y6?oz^wobRZ0%O_?=gIpkkk8az%TXqSJ7lb!If{(u?lWk6)XDgsK4;fbEiyEp zE!5YaM8*t5_u3vkGL~;&cf5Kk8Nnr+=G-zPBcUhwecW_1t~!NZ^O{KpH8Ez*bO$mX z$&L;U%_HNJ#%aoP7cxfqUI?vNLWb$}Yg0}yCu33D?fY9-k+Eq8>Dat=WW+?#XKVV9 zae>Ky*S?92(jj>^Zwnc6_s6R)hmg^l{_;=6E;9Z;Z6muykfCSzWxHV%8IIe3PVPNO zhF8f@_5FA=TnQ#xQHZT|tKV_PnyGwPY+QDctt+F&P0p6%!vekP+*|lHGkp z#-&8j;iD~Nl*{hh`nHo%t?}%Q?I$wYeVQ+g{7Qz>^$*LNzLR0t*8Pk6iwvh31Epzy z$ncID5wu;U5|PYt6I@1D;_T2Q>4b5WD0H6|_erY~!u07jqDhr#dOExLntmmI8Ms`G znpO#&?QU)>%`0J7;@M|xU5V8_ekH$VR$`~qwk>rIl}Jk5tyMI?5_z&{apJ;Curv}5 zZFH@~Q=b%TJNHU-=VmmGT3reCw(RrG>ndS2$!TKAr0R3crbx49&<5=9zj8cA0xA@Q*l_~%vPO|IjC*|#h4 zyA_sW?p0#)j1`S7rIm1qT7Qm7uEbjA=EWE2mDoMB{o5{fC6e9u-d#*oB0v4mCS6G- zxK9(+zbY#6+~5pfR$Ga0A(x_VJgJ0cNv`>!=arbzbLY8dQze!=m7leETZ!OA#-cy% zl}L~Yx}Sck#8r)&0@Bw?P<`tCl6x!hD7R@$V1Fe(wYGET{H?^O8DICQkD$OL>X&KT zXbKiFl^fXODA+VK_RJ-13S!)K7Veo$!G(0AudW6Zls>h*t!GSu%)s8a$DD%J5NEY& zYYP4rFXiOVqCl@_b!6Ne3g$Zbns_-;;FY-bsnsG1B4oQx54losS|e(~b9V{~eB!^9 zucm;Xdour&7X>d{&v|e4rQnD4l~MCHQ=lDnlcfg#bV*7?GzB?NYn|L;DX21d-f56P!IO}d8{dyp@TK_U`uiyqsP%jsS#XvDbEiQ@ z{CNtNB&vpQ$f6)Xra5)SH3|-COnIz)gMv#wCdrMrDJajKF~9O21=X!{KBbZ$ZkSHL1cQkVRkhIXP?GAihM{xp~0~v z_j(G1A*UUuzM!C~_(I2zCJKIaU(2n1OM#A4!P>j;DX>cCmVH*Uu$KPrk^ zPv|yprb1$!RznY>;?4f7qv<=S_|3?(58X}0-BGF$M?RWiO~WACl0Y_nL~* z;uO-sRw`uO89|;OsA!#kb%MobD*hh3B^~UhLQhr__w)x9b2Vr-qyZ|te0a?#hp33i zm0b)NNyF*ZM{aY*&~V$jv2Uy<4gCGBC9M-^c**G8!q%nX$KPJ9ONKOPyZ;sMHle{b zeRQm=<^SEcUQr%6H-m<&`M(oh%%Z`wg`7#Z@nL9n4 za%dQvx2@*n4H^W6n>`=hrr{9Tr%wMK4L0nx-aAWacq?AfP+37kcI{ICKU5mlJYU$f zf=PqQ+j&7*TpAdk>|5UmXxQ63Bh*GpgV|s6jziTnycj(#TvSU#hPGap#uFOc3?@eT zH_*^;uF-qvB@Gp`M#pwH(-7*UJmAtsgZ|Rr38z2MP`CQW(34IYPWgUKn%Yf+OVCHv zJ-syagtev8f6;I+`n9_99~uG^o@cC4rengX$J$p_>5!klpVc;oj-%J4x^|j$%qbvp z;wRA2L1G!nCev}Fikh!&K!=Z@%xudvI>yx8D=adjgQ&k{-D^e1!PmK@#kO?J_;A@S z&7O{y?(>v}xpZ6^IOAyQOvmbx$;`cr=up-?=E8ENgP|M8AL&jZ*h2C=XccJ!uCWJ3?19Yzw19Z{JdX zQUo23;ubdb?WZFpXgLEu7XWyD}n2zqtGeTb;rK9ksd51+39RbBt!=p~op-r9A z#XU=hj5{Gp?L7V8_sh6m?@T%z>PE%h%%2$m?dYFl-wJf89kLWPjM9mk})6o!IW;X5x9p}RD6>e;zW9h+L)&(te^c~A3 zb-kmbEcLS8f{%1;zjU6G@`a8m*UvaU?xEx1-IGkipL8Tw9CZmBpu?FJ%cuUOj;h0*5og1*hH0&k{$>mr9L9GE{N^vt{}Om{W!c)Fkp z&8wz{8!xUxmd})~@MTq4xphJm!=nmAyTM$YE+aT&Ti$jeCWo zE>}S^`IhyDYgLFh&L!pLRl$D7WxJ2Js?a{?JY{ZS6>=AyaXem9h4mgMnfJ@9Flzl# z7hP%a_3LG2j(blhr+k0o7v?-GxpJa8C1cPA_C2>LSyqehCAXxs?2=%Nc08 zLo(a5l7UO*g@tra20WPg*2)_g_#@0At?^@kdOyqV>ShMQo@G$lwlXmFO{$~ab_VJ@ zl9=(k7&!AS!9^C%z~aF}eC_=V^r}X=Zi!)_WPG@!D2{<`dZ8Y@2@Fg!3#wUsf`R+C z{+?+m44jzfUDt4ifq6?jy-hP1_`K4iVecgdZhO1>v#v0(af?gS$Xo_Ab~y%l-DE(# z->x;cfPuJp>(CGP7?_o8*5N>6pe@}f{1};mYuUP8)iegy-O`Sl#A5urUs3Pf#$$j@ zQ;j7F8Q9Ms>hF^>V4>(wSXRTptHb~8)As3L;K`(#tec-1NS!9p?d)a%*8H4#KN$G#z%)wgXP|f? zCBOC$13~U2Gd&e1bi4`+cZ_19#y{VhJeG;$+jB_2HJNbSn`P%Qk%`Vj8I()9Ox!x2 z>ey_+gx}dDrnNB>1JX5xaH zW6;n>CYDdIYh4w{#BT%Z&@0=RsI)Nac)Np%opwgyGj}s#=&ainx0i`0OSPjUQB0g( zqux9I5EG01RAU3;nfMko)L(dviK5;83Ez^K*m~gm(4x~!Og#D}=}a0E)u%qFK0D7u z;)S=V##v0v&3UCBag~Wr1a4YQnNSysb#scD5Z3T=+RKX`nNiAI|w=8G3hoOeucF@Md3+u}p~eQ%iq#e~7mtu=Yyn0UNzljp~uOq@Erp>FPPCKjGt<9&RHiJo)r4fjW|aPP_ze_b^e zw&X8p3K_@3gp#>I7RGjs?d8v6LC`lccI+G$4k{1z`_5;s z{y{Z)1q;d!Z&M$wVS(ZDO5I>R3wxJ8%h>75g6XaO3k)=Xd_rSKZb|2PXS z(WGt(bgshOoM8+mi@7Dm~#u>m)&d2`uN zTXBt~=FA3veWsnyA~p^Lq*HFXvSA&1+Ocyv8*d^{Fz2meBP;HROVV04RwW(e*KS~A zC~cps-bOa)mv>8c1h5f)bGrvQhz;Z7Ej7Qlv(Z4==;;y0#yRf#x=RskxXM<0H}7Yo z@6qxG>w|2RH7@o)c$kflHs>b6Q8uQ0nG-ZFiH(OpXSQxU#l}gc>7fN@*>E0f+R>H4 zM%P5c@CBJ{+%cNmm6FZIW-G0z$GL2bpFOtMFrN+Sf|0Rdci1@M_P3u}#D@LaUkO7b zHr{XiHnggejqAamldjO&SRejgwUx!j=z}e(Gx==rj=fYrEMg-%^=XDg#>VtZ54E-K zv(cQZ$O?GGM&@0S?!9_8R#tFxzCCB-FSE)BO>EFa-5trGc zIjA)|$WKt`AZg}4SH*Y^oaXPAOw{4v%aZLLK~p#=Shc05#E=7jpN*bBjXBWVy1vfU zf&ALvoCO?sU7guFYB2|+Zd-+J zSjGXj)TATNgM%o#LHNhj99R*Py5@Ru@LDlG>bMUFm!6F2y}yZr6-^^zb+>Tvr~OZV zNH7P~uYC#Sp&ab)?-}af&B4?WpOcnHa`1HAd)1549GsoplG=2LgT=-#)u$if{JXDs znsFeJgVH$*VQo;zfkxU@(#{eN#95c@DCHa+zHyH7hr+?EqEn75863Q; zJkHGGaFD}^cWEIwSSO9)+ekQ2eH7^$tKfk1B1|HBz`_335D(4A99Vu1tnq)wLDLUE z&%2ErTpaSM>we9F$C%aL3tKrDoUpv%^m`5{hKv2{KXDLd>D)B+D+fmQbAtAK=b+ws zW-GmqgEPylLX`(OfTu~v8YM1z{S3mdj^v_b+oY~GH7rQ7RjXP7#1G4olicH}}X+M4BAtR-Ar>k#UyF6UzHH%^Y%N-jnY(v7Zra=}(5 z=YLqw#lG<+W^??wu+Y0xcx*EluS{=P*KFnDg6&n(r0ra|&%0z7yo-z9i_cL=;apU% zJmvUn9~Zm4k29CWaACM5-sOB87f*J^@Lwiyk-9(9)$#-vi{is1Q7K$}KN;e|JHthB zdSH!O1{Yhi{XBgxaiMd|tM29%E^10vd3WY;kw|lEn0J#4M?U~!qnf2$T%9zg z*Sdxa&uJrK4?g5#g!P|(!4oc64t)vZo^uhouxDuFD=y64JCh3DaPe|oyXx0AF3$VE zNnOyv{dW)EsDA46zkNs4XFTraf_$h}+vo=uJCDn=!uq)|I4jhp{o&$qCMRb|nTOMP zbfeXycvyI^BLB)*9=?%_%~~~iC}I~B&YZ}@7D=A_U?9S zO4!DO*#U2?&cxmlx1kpULM>o7$E28SLMe#)He#p9!4vJoLQo9vYR!!@Uomk~UoBA+YD2 z>W%9>O!(cL`Z1pe#mGi=$Gbcn)vV7rUd+RsDYe=U%6Rx_BG1yL@NmOMsJp$22OmdH zPC1*0F^lO&1AHEc6&3l*#XKC`P;7Ql&clqrg2JZzJhX)7Sx>Ly;mW>j(t)QutUi3v zj(EX?^2s#H*w;KT&ZRi|z2za|N+R=CI}fJ$hh4sWjG+a4)L(#rLQMt1Rr|u*46zP&BvoJE4^2a<0GYSSwof<9}AS> z-=f1u_gJSU8$CV>bsU0Xjra(dYTGI@;X`}6WvHelAM)A89sV=;NN~{)zdMT$hvhn5 z-41+otksNKIG>NaO=`WTUHI@0QHiZz%!gXUVEFEU$loR6|Uw+qvf`Pe=>&-!^PA5*lm zNv7$1Jk-Bv7kQD7Wb-r%>oOnCvr-&YbNKi=KauHmgO9tehh47U=412f1N;w#d}#UZ zb)8enhivOE$*~GPj)Vn!)KK}D9UV|JiOI)@Bfg%&Tt2R+tg9ml_}FlMrS~t%zkOd@ z*08La58mxX{^x7?I6#`;^ztzu)2n6&SvK&|OxUzWz2qaS+9H(K%*V?5X&q{9d<-?| zh5LNqgZ_SE*Ue5o_H=7Rb$;c;cwltzyk0&UMySLl{o*4{W3az=kPlbgp9y-(1p17- zhjyqEAkFxcL>@yRWbQlFK@9@BXioK*K%my+g?i>>0w>o$$!OLm;1uvcd&V>ZT{~o1 z2h9lFi4^DxtqAzXv2(`R5*VLEGumWNK$=#bUoe+I{N*CEuTBK)Z{99ku!z9>qU+YD zTnXe-u83#1g zfwK#iHGDfkV2Rrzf1D!l({p~)nX?2+H_i@vmO&sm*rwGulfdL~i_nN{0uN%Qbue-X zBpuTWACXUBe(J=owRZ@7xu_A9Q$(O3cXV$%iNL12%CWO6320XQ?oXf-kT8EFC|Cpz zi@pv`(-6~1g832Aszoq;L|LH?oMR^q;r-xri~EbrjsbAQdNNXMT-NUj~2jpDc9)F7y;}( zoZpy^6F_}URsIT10p5Gf35nDaK;TO@%bFlS?&g{8EFA$31(g)OnJmEi9adqgdIH$) zzGH1;Ai(HIlTI%q0p3R6ARU?}05{ek^16ut*$G$eL}mg+pP1D1!9oDfQ|f z8y_=gh5#xT&N=$q3h+8>Oy99t0x++gV&1VAAoIqE_!sw!r{pcblF#)QEqw(r z=<)Jivr&MbKWcTN{ROBSSk;gnAQ=ANR)A9@m;1K`2@pI+q(3G^fQ6cin`Z40VDdz6 zp3hDJdUTzG;==@ZU|40Gzej+3CUaV)5dtJxk_$Q`1=wOUGgLcDfcf?%mh%q?Fu~EP zBk-UAU!3n0C&db&fJu1KVF3!3-LQRlM1W%}47$D_72tor&h(B8FlWQ0sKrSFXl}gZ z5SlE&$AIy@XHN+r*>;XjP8Hxr=$P1svjQC6bINJ(oB%%iMhuvo7r-v&ILG6n0AmmT zN!XhyK>N{nWL_44Na`D6W($yWI!4%hO@M>vdXh$77r^Ucq}z-e0@z&XRQ1XiU{r3H z{NQZ?T5h(d=H3y2b0=h#uuy<2#c$L<6bTSj9#}iBRDdG^Zr85?~@4-_RS&0s9Jz!eHCWQ?h7z=(6&9SR)Bt`lEREf z0z6f<3Zp$1pnU8d>zDNcoEdM@IrL0`9XdBi7B2)?tY;9p`lSH+Q?J@ZH3{&;bW+ci zW&s{qU7~PX1Sqi`AJf_@K+5cMj%w`!Y@0WxZ{`O9TwG2uy*~;tX~~HA!<_;5fO%VbR04kqpdAt^IWZu>=bfD? z_XY)!hlhF8{uSW%{&uzRNGc<%X$X;b$;7lYjsQ}2;qHwRfD^Z5VLN{vm+)8F{aSX|I!p8-jxXT znfgNTE70`XP>5?ZPTq)VLd39~f~-x1SjVRuuQL;3hS;Gs#zF{HMMXi5l@M=e$Da(WyR6pfZx?G6oS|bMP zJcOW3I?nmIQiyc@KM96wgxEDL9!u8>vDCb8XqT4|M%FRH^bJD%n%R>?^%0`pA=0hU zPl&SlovK5dgh*W&CN~cdV!La5>Z&b5z&&Ku{vaXrR=-id5-ddTy1-h_b|D`6G-kXF z6{6VR&ues;5XoEXwP)@TB51o;{e}o37KGJi#YGA+Y40k(8~cU$8YS12L<>=K$gS~X zj1YH^2y@2A3UT~60v!(v5s<<$3OFK!<5{OS$Bzo3b)KGIm?%VNrbEbs6GF(YR+xP| zDa5Tj+jiYkLL}TSE?krY^c6S3BM-9%!isW7juPBt50)eBs*GhyK`Qv9IK@y_5|A1?Ixe%z|v|m+Omg zs3)LykD&m)*? zqeIA%1tP3mP+@j!p$KM+ZQErK;qS8I!p}=YXk1|#Hes0vG|vKSr{yA?-(cLi#Y2QW z8}mpfSBl^ks2^FpT7+rA*>;aSMHmRx>FHf3!m~XWDf$~kAn((RS>hu?T1=YbPCpSs z53BX1Z4zP0u@q+IW)Tcds>DAJ6yfLTM3+BXMR#K=$) zf^!eJX6+JTVScv?d$$O>clS!(gp2U4_*05%qzDhncX`C{Bcl@)w#N@gjV=?^~yFR0KuclZ$f_MJRZ-&fEWl2*+MN&`CTg zg8!S94R=yR|Jcb^gA;}?Pc2G z7GcfuG&+wWg4HUu*fyF7D(h05#xO)^@>Us`#S($BDUsvD5#iF7!Gw682zx^gBcBjq z#jc-2QsKYzjyNFf6pQd@e|M6$OoSH)_qxqjh(L}1q`IX>gp9;p@}vhM>`s1{TJ%ta zWoLp{J**R9YDTkq?-LRFF9p=xnUI&S8F^ff!?)ekK}C731B)1Fp-A#UNa}Rd$<-k>kEsl3^}J%<4}mG)pnodF}Fe zIbDnyzVFnOY{XFA99&~DQ;e3a&1tLc#9(g^@Qj)*#^taVn%Q&2*uU4ej_WAK>gXpI zTjz^m8N1F~ZGjld2@iB;EfnL`@s$nU5ThzZmVJ1M7#Ghj^S`-NjEM6BeW{xm9$AZ; zKDmowdX1f@wNi}1y!k=%R*UicHtpXhxfqmtvs+KB6C<5eUT|-{7-5w*p|##(EUhZC z{O&7;5yzrq$|f;>5x0x6S&VwgwD28)Vw6{3w>`U6j581Qx+=DbvAzCE#j_AG@M2=r z?;T?3y}saJvP+EKw;H|fyTy3){wzHrT#Vw*(Xp5Iijmx%?8MwBM$ixCf!9%DEEqV( z84)AKq@mvl)`!IC9&s4!;>4&?`#BU7FUDPs1HzntWc&Dp?j*r6F#>e=y1hRx#yrDM zs^d48DK08H>PK$S`8&8WNn-RR~=ov9?*)^*dq=}I*H=tILF2+XZ7a3jW#h49W zFP%$bj9d0ZdqI{MA3WC82VD_EwB|up%2hG)*01y{%@HHcPp13$x)>V*mNovoA;!#Y zf*iwJVyJ~K3S3$sM%!+-QP^EE_>uG9q~8FSE|Iwx;(9u!xUp*&UMmTwiv5!>P3#`iD7Z) zirq{?45i|UJsX5#yezvwi4%)Kr)b37kcx4Ean?~H7bBcIy6Q+NQ!&mp{z!b#AjYodXxDEq#BgoVaq{@zyev|Wr4rRKDO4`M7*4e&JnB*v7n zFErghi}794w=TR(jE57ST)f!*@4Tn1^JaV#BP z)|MCejg{cj2Aj}?aS}*37FpiblpsIQqC>7F!I9wG#a|{!;J0&H_{2#P*oR-Yb=H+Y zeZOAURy_&c$6TpMHjqGgcw$tEkp$O|U2v$ICPC~;jou$7|94M}U3q1eg#@E7B|Gh! zE`gVx;N?vl$*{+gpv7wbG;Ie7rdLD@S2{`%HLrYol9L3i|6G;?%F8W^zC(hQ54Nqc zT`IvP@9Up;ESKQre?ChB(~ya?8&*rO=i}OX&2Xk*l^!xTcQMATSoU?O_E?<`_gNVQzXNVO9J_K z*6XF`Bp83>N8*w55^PkNZA!Z&8TMQfbZ8arw7B*!zgOGTwqBPYZrZfcD>o&>u1kV3 zv#xv{T`0kN`s0gBiY2(VP-E_4l4RI-NnpJ?Iqnlhg6JoNzFCz7Z2#Zf0G4Ffc}cJ; zELypikRbDGS%Iopf>(z=E{9Aq?7buiPuX_1;(-LTao24-9!W5m`TU`Ay=2&ZNpSA= z+KHDMCI7tSS$a*A1V)w8jU!ql!~RQxGGgKNnD-L=aF0XTM+xwd+S2}6GVH)4c-S<< zV&iuSChyoOKmSvLpw3&n75$Q74<^Ca{wd|oN>a=};o}uOLW<*R7rvE@l0xyH3zI@i zcl3Ni4JkI=SlZ~LB}Kl)u_NgdrTF-t50hf{+#keWJt+=z>_g@oNg-SsF>=4L6l4E$ zVp42a-?rkdr4%_YrgrMhkfJqcbNYsv(qS(q#evAjkHvGO;QSaDGdNF*5l2@{c3dDG zc4Jaxr4gG8mq^hx;rC$kGAYb+=H1lwkPiDXDdDgRNK5>oqSq4 z?98O_U*L1$N4gZZh36e+U6kU}iu+G@W=V&=nH2FGm+HO9kwVmbl%sV+3iX|=Eh}$H zhuxVJxiR*)+4rPq`yKA~y+jJzllR`+luL*GnH1cMQ*H8SQjFBwT>Fe6#hU!f`!qSy zVTUHg>+*55M4=QG^H%sVB~t9=o&41;mkxV0DgNFcSQzq93U}{$Z*m?>k?}I>MEz4K zp8x04q?q=xL+bWQie3B7!;UsfQPCGPw(6}Ezy9-SQY;zue63Z76shSGdV)HosMA@S zd9_PA?9`+PHj_Rd^FxZ_5~cX1zoh7yz0mN;Z|SgClOkyeb=zlU8EWKRDi*3TOjtML z?iMxKuv?R%V9TwK4>e`zeAkpVdV&mdBBsn*GD$Y<*JO~ypO3CIkYSw4eVtEJW$-yO zieYLh8+L3myvsiNJ9D}WGpDig9@xlmu&{Trs-0}uv&o=J4^KYiD8t%?_be-%WVj-3 zReW%f4ZAiOET8ym_%D?q(*JVB#pNzaSa<^~<)O50MRfHyO@Yn-?kVmf4Ej#F zpQHB4hTWSCCGHbuw;hz>Tdilk!C@H|_?|iFlOP-RZ!+B9uEdrm$uP0gCGc;G3<1%x zf9IW%4Ldj)J|FK|T6A89xoR`oTQ13va6a_3-j#oO`_IG4pmC$={@Lp?_*zVf6yKB~ zkL0a2cw2_||G79BW^rnK_7=-Z`df9r^?W&6l3W{u7sxSv z?9n4Ri{!)ZP7cDaoW%S#SZo>8Zzj~sj4 z{L_Or%0cz+4RxJ6dunw6*J z=wfPDempB5c6xFglc)aJd{GYhT1j?hmK@_7hMXT{%ZI(595>$`I=1MB939~l<3qRP zu>bx=RDMT3?DpgkD2I+VF8!CcGx^1v%H>$E)wKFzrF_`$$?5&gRKRHg3`xox}CC9`0j&JV$mSZxp|3u4Q`LOp>;9l(q>3UTKzWSJjomEp{ ze$&>mVs*u^`%^&i>Dk)36BN*j*6xXzq`;>BHJNv(D3Jf3|5M=OXvy=*(-fF}e(3Nz zQw0uBb}>w~Pz*aj1;$!Xw)NX6uz~bN#lcR2oH^Ea_c$nqJ)iQ`wx z3l-pa>CW<8tQdBI3akmrh~_U>AnSd#PTvXzn)a$P?A9oTeV_sn2?@XRHz+_K!OUy) zRp8Ir?~5n+D~6q*0_oTGB%cgYpuxD%k{hDH)S|bF-cZG`7gV6UYLmv!2nBvF%Br}r zPl3hKN9&(QD~8>m0*~t5udFz%z?6U_rxOVZY-{DcWF;zw{h$Kf-E*heoKnCkZohE* znSXg3D&2nloMPA!Dxj@tR-}=sfdA#KtK6?BaN98F^Ra7+VNa-lgN^p=?wbn4Q`Xc^ zFZh=?=hTD2g^FQUsDSs1A=VR;0=bVC1geu2X!AStcNtYN>%{~Y4Ky!B;%H5pJK_}HV1oBt?Ktg1dYQ>l8`BUWR9PEy>+ z(bY)G<>@aPTa6mie(phy>S32yjm@+7D}R_+joZx9f@!+d=v@3^`6m79VV_uyBWt&w ztv0TPY_74vG0+w$z;0;}p_r&tZI|5rUdXI1T!|7Sh@j@<|LhejXRy)W|rul4jD zO1Cd&ovn~ z-|d|Ru}$^5vh=A1(X>`~p{LY>xG<_r?KRzkh_l_YE2hwbxGJ#PG%?$P_$(Utigo(m z{G4&!yW%VerwNwXM!PKt`=HHrFa0fuK>gy}o2x7c!LkYc;Vu@$rdO9c7Ta17)8Ee7 zV_;-K*!q4x`d7n(2%Arsef(=qyf{q?c=gqsIPOlfc-CT0Y&+0!qT#VQv0roI{x>3X zf~&iy`zzI)a611=MeVLR(H*nvp4F9q>l7NzYfqaKzs^1^JQ`Lqkt6=a z^A_xg?KLBsYDcv8w3`uY)_mwW@Z5~}<6|6K2GwM$(=j;Of%!sg0vAxq*m=V)0 zj;;Rd{I7p_W9}=PfAQ^&6@-x)!QErWD;;k}1a7?@ddp0S3zIUbhWVz%i};4NHJ40@ZENQi zCMTN`Q!bDO#D`1?wIz$+j@@laOyBz1-#@^VNd83KQQ>Jy=(Jx})>vXnGok>!j>V@wIH#piBq7&IX|9*!%0-DN^_Ke*?<|E&pO@4O{u z>{An>XyVMdB$)}}6R3O2mt{gYO*Gsyiey4~Ik~j(ZZ zTer#hZKi<z+XbGY43Lqs0B&!-pJlN1WVr(=6ztu%i33gPFMPd~v0i?w^g?+K^pcylmrRSFFod053yI`BRbt7&^xoTZJls|COvLu`3FoV#2ea(p&+wMSSDq zbQK^kS>4QQtpda>T&jLSQUN9|?9ezjFAr%E;pUE=^3WB}Ne=!Y4_*Nc64_bEy+4x% zqvc`n@5Bi;e|adgI;0+ACl68|IyR4*qCD{O==p6}KL0GXUkU5?{^GwbA`e$GwIY{J zFNCI})dtnA3n8sOVMTb!Lg>;+T4J5K5JH`rcI=8;2uXqa{Z(%+1Vyib!F6^ELCEmS zHS)+pu$%e0=IHi?pmIZmeq-fAu;{)yONlK6+N$-k#GD)udxVc$bjpF3&azKm%jCd9 z_Tql>y&R-{@0tD-kMfDA@1{X=;9qp(^Nh0`NQuTZy|R!4A+=!zhdpv&r|#;rLsbq8 zB~|&B%F99M&m~!+gdC`y)w?M&BnyVWW8FCQvQVR}a&OBQS;#v)Ty`!^7C7cT>k=bm zVQ`03)_^a{jH?#AFJXC{(3;mrWT8y1UO`z$7W6LHZ%kh!3o1{A?R2GO!7J14doPbH zsKv?{#rDg9#ku}*t6CWt6d3=wxd0ig`da!G9`AknnKwcPw3qA8^4^dEFD>OGGIlbc zvR1QJ?U)Sc&3Rorvr7iblBajZuakkobuZQr$jN}>gODA22yCZ1=}N}1G^mV>{aoEB z4P||<25*a`q3P4QfBWA`gN3tv)<~>0=rOcop9e|9L%RpBF1h0I>_sg*PD_J-eY4sM zBWd9LseU891^J$KKvqc_LNA+3tP+=ov+Uv7yyrYDqHOhdwNl_e+8xCF zA_YRTM$A*MrGRMei~AlW1$MF2sol4wK$|Skt8kG5udF#ETPrEx2we0OGm-+ets*~* zw@5+S{+?Z7%2E)cHrZq?DFrl|R`fo)6y(Lf@z(w)34=>4{C6}+La5iaH;0Om?@YK| zvL%7`BhKkrq9km`KyH@wTSibL%a=p3) z&|Ewo>`|0}v}*^K7fN9J7Up#a8Q2ad`t{(jIM{s_$l&x-`*=r>R zUC&3pS{}snO2hWw)DZ)_lGEcks#yN)!sq#=Vi05gJj+;I44Q^zzh^M8p4CTU`Isog zq(8xxktj(0PU+tGOBD5H=uYKVQ84^ud;EH~C~&UJ?%$e>^1^uwz6enmyeIjn;}){Z zYP8Hv6zrrsor-Kl(O#ltDvpan=z?3V-TEk3Jueg1#&Q>0i|#tCKYB|tXt5~JjP}?! zi;9A&KRKYo5CxS;>pNLvB5*eFvx$0_2yl*=ZZE180TphJ(fMK#$i6A(E1fF>F@N{} z{*@{MgDWx)C&q|C-gbjA?=TUVP*-?=!B+%A1(wF0bryk#xyMr-&WeE4o|1;!N3s6v zOB*uwiGcQm?#muh1ZXmccCTM20*X@My#7l>AZP+~tmNo`p)hD?hWhN%L3zWi zvX`5MK_w`A+v1hNQ1(go*Hd|6D6D=Hsv&~)zIa{y!w?3cue}+1I)NWtC&QD(FJM`TgTt!tesn7Cc~( z{|7m>b>Ym{1we~_u;I|_1rRfHIIJXc0dPdU)lc4B0A4L(N)iru{Je-|m)QbHd*Rw# zql?EA_dM@hkLR1-$PtrY0C}I3#4Xq;cWuu4HX;ZXpI^lqH3@>(u{)LH#V9|oH~5ex z2->gjh(^Z?LY_mKLQDwCy=so&hV|rE3eQ^OarVtI(}RK#ld^AVg*G1Nb&s%FEeITr z_>sjjf{@lPX)(?e1e&w&tL|X|P!&Wx{H1 D8Oj}IQprn(ut;e(CUcPnC^@ygZ58mBuG&b_)0|UKSKLsZ~7;f@A#($0vW~DtINFL{dM@r2z+D3d( z@Z0I2|4wZGgYBt4bv|J4uo%3wh7aOs2R|-Y!Uu`p-^7$i@_`-gz=5aySpL8!E|QB6 zQtX}_el^Amashp64L!VYZ^dU3rDk3TTKF#|q?#A_xMC}XOL)Ql+!>=+`MhAbMs>hC zix*BkJMOBP!V97O_j8oud13F80q^yXc)|Hq?ZZPMyuc1#QW|!H7ygEg_xyF^g}RLT zk4GJN!NW9Dq5Ui`{D`E={LOj6COITw=V4yhG$(alc0VtyD2N|q>+pgQ=O;^8lNWsa z9(k&)=LMGVKE<<2yf9fARG!7*gPD4Q%~Ze;O-!>d=VEpu$R z@p88A*%%vMyp5xL`q?m>b=%UpgAF6s#vXV!u)#pietmQ`8yt*kSykWJ5Z808hw~Nd zt#wSl4r~bC-m_*Ziw!G7eyQEW2e92#BNLR8*|7Oimu_<`8&0#bCQ>5UU`p)|3c1S$ zX-cd9Mj#t<*e5=O-eAM_fl7Ls2OBnB6$ol|VMA|c&Z?F6_#n4jeU9IGHt0&Fi;SIO z!%q7Gw(D^=SbuAGU0}inqut6{?+n>+OQu)Jc@G; z!}FJHdboT&8xC%mTe@u(8#KB4rY2XF? z3)tYlfa7qCmkn+$+;DGeCg5d(NZ3 z-2_+{2=;Zf6Y#}eg&Juhz&6^YN2DIlvvNML?H2)VXC^B>%L({ zCZMb#e`Col0y6h4HoZw*Rl=~pmE^P+Hzw8?jAcaGGRc#iK2|xEB6zy{(+z2 z`8@>a(AQ*S=n`<@@VWxtodjgKZnw7GPQY*cNrT$0c;9AgxQ(|EQ0S&pSfzpE6WEYz zrG`;#nQzhDM!arnM^*d=0zQf~e?72{04+xsf5Ei`JZ-F=uUJLESFS@63CaWn*I&x; zP$YnUt@Pumz2?DmX)sbrA1Y8v2lzNF0u%W)BvRs&eAQ5BjH9`c`PUuX9 z2@+tzZ{9}}AYfzOHTo4k0$McVlqcAzM>7BRcoPIXRC;YF$|7K!vf+wU1_4>2rtgjD z1jOug8=d9GakIWv@R|$tM&PQb9gTn}rp5P-^DI~$ce$Bojs-0j1k~$hSYSBwG(UTq z1s+Lt*3nZe_;U1=RPZDVy7$ZF`A)DPyOOrgXPgC-3agR>##o@ps4WN|Wx?*xmVqfF zELc>zwXbNH1)4E@RRcpTm^tdcP7bM`Z#`)6j|F$z^5<_2vOsqva)0Rn3;J0~g%U{a zzi*CP_p^YvA!&x{!|OYR-d%{)F%lB?>1Dywq3%-yJuJB8H-GgcQsVXQo%P);Xg6<8 zH$y50O}%RU%YtIb_BHlM>p*^a26Aj<+%~$41%A(V^JydJ#i6UClLg98R7US2Hy^sg z(MA4h?mtRHuH0dN4?0*7|MY0CFVa_y^~C_0JK7t?K?-|3RvT|;fpJHUPX*HBwvc@$ z@@~81zeuFTu|J~$$U$+FOYTT^70bX8dFIn8-T8CJkng`2C&wZ0TgH_1Ax{a@k8Z(n@To7;xr6+)a?JBD^5DxWb31Xo zvXvd)JVmBm5Dw$Rag7s>2y#SLoU(eXW64PQ@E!+gy9<3DmQtHL3;tsCQxx6D$h(G=%q( z1L~uG^WlCa)K6wd-PJba*W%hvc_dcD z0`04H17qKPQ^~c|~{jvn4j6|kzFC<__F!`bs`mv{v z{GP2`L=oTdc{3XM$TB1MQ6ke3=w~~N065;|W=x3L{i@#L4ngFS$ zS2Z(h2v{NPsl8lVmqYcq7@9od-b6bz`G_`w468h_x&KESQR0;4k`>4dTi2(UU z&fV)a6A-Sj$jK7@d#9T9$2;l-+{qg+`iOqt;9Stn5lsT7*YyuB$9QntxM%(dNr0S^ zWl9jn3GEe|hw`-vI1T1s$F|{p8P_aZwFCXLTr%^d4gsD^+fP2+iSwxVe8Z1j1a#cy zI571Hu;1k-N$w^f?B}?gGsYtg$p*L7y%?t~yjFMa!}>ckS90_T*mc+G;t`Br4f|sT z!wd;+uHah@->_Fl!fXMc{;%RG$jKHqFVtL1UyTYX#HeK zKnd-Z+}ugbFLKweQ?o+-?3a~wIF0#0icRH%VOc6Uf06m=xvW`Nj@-A;I&*AA42GXXM<7FXW55O9CcJ*D#s0poqPrm|NFh#YIUtnWs^ zfNH>$7v>cUe=2ppzJ~cktDkC{2in8Ai+{wtP~PyiUdJ2bdZXrQ$Lj=q`rP&{4)c;r z{ZCTLd@(LF%qnJYV7vQeH!0sFK=oc@*&%-diWX-?c?Dp8koh9w`7Hw1sZCIE8|VMa zW%f)Uj=OKL!;)Yehu0QkyF+l?0vCl|3MD|+b%Ex+I~d2OFAZgd5m5g;B)#q~j?bcx zmuKz~P(Eg>=unK0YJ+hydOjtBUVOpnq7sc5C`$0xW;% zKPZkw{f-qcX?ue9eMPK(CJM*(Q%k;R49?FLc~_-av}+OD*)4Gd@OV4#F^b1LMPcCB zsRRQ2o(HXVOvHR(lTeP&QvzZaYtq7=5nv)xI}w$Hd06B7nCHm^4BgEd&Uj8h+|ieP zxi4^@uDaj){1W|;LlB6oZO41pG@~={D-fM>j(6A*JPBeb1|iqEF!?$g;D&k83^>aPU37^Zh3RwsyNYJuX21JMm@5)z1VN6Bha=U(lbK zez~*$EBZ$b4_^L4Y=Ixi(%PJDxNRz$u(utL*XDSskC#wFtKRkI>?!O3lB$kxF zx(4jP@h|2=1Z-B3rs*_t)l^++4BRv!O(vdp$-nt3%&eII{DP$o}8Q1M* zyw4vUdeSWfbkAMUPejhay9H}o2~b+!6z~?QVX@-9W*g3@#`h-_(m7ai{nmB@(&>8~ zbCFLZ(k#_G@Hp?o@z==G-UAzyJJG-2QWuFv=AVzdBix1KUD%WAja(e!Y(0RiF?g1J z^e+M1wL79eA)l>}mQn6TzZ1G`c`#D$pmyE>@-bbcL%#>@I#uQV3*?^$>0(AN=EYC+ zZyZN5&mG8rjSOF%=FaFNz;wr8#v$aunK0&&ijBXJGu2%wUy;Yw9etL83~j&qEeffp zx=+LpInFBgJdf-h6cO5otmH`LEB_~-=?i5dBUi>=YsCAzY^!r46DgD3o#>C;s{H(q zDe{)}W98+@gh#cv265cvZLZH`A>Z)4z2=T2>EZk&^5gYA&t`D^3z!povyn5t-zMyk zE6jRcEJxlp3_Mkf^Y+#M&%yxY-2jo=jmXpME$=qsd|!RXwbu`MJZa6G60+~t(icUj zA2Jt~7@kGWR@(G%BhB--xkaEJtzFhLy9QaVl4hTW`c|y}yHX!Hpzv;49qOaWm4_!! zA}<^(3g|?=ot%h$bPlPN>=@dGdj0S3BU>wE)KQyNe^Bpj-|2ieMkY1-7=A{($fx|% zRgpnvnFBFsHwI-F*u2P#Tcu_%pwGyvkD_&<`A4^P_bh{X%GE@x84` z>zrslck~;(eM)v^q92ieS*6p1{v|t1n=X%hec_(Ze)Kz5dyeM2qCYy$yYNFC`lm$} z=Nt;rUp*};ZRtiowjr{Pj~AKNx~O;APt=ENo3E2)xNaxXOAh=%J!x51box61#PYSS zj-{A)J-4p$EWvr3v+1~5OhC0&U2I?xj$4&6JNO&wF(U?KVoneJ-8uLJ~aH`4I< zf_b0V0Ud|W|HYT^;gbb;Ujg5C8Ga&QLCRe%4fONMm*Xoq==YTqpYhT$9(+tNY-+_g zQF6}kKpw`8Y0sl4A7LE1Kwrv#3FAs@+$n=??+JK*_ouQD#+{*uE1vwuc;q<4NsPz% zv?J*1wzC+oMk0solrVn1pVeF6iE%A;p7f1PNBz7Qe%dq*&kK>L;(LYhLgD()oD`f# zgV-ZyUlK5%uVpLz0`th*@d^cYQRr%O90^CJOVZ6V{Erk!Zj3Ne?bP#&)-B74<(t{VAWm za{eKXQ_82FJIMTyynIQz6Dr-M- zE(rH|j-Gy3cpLMF@#%o2w=mvo-YIhTN4uC7xcbwNfET;X!WZ8_dv&s2^xrf3gB;C@ho_vP#ExbI-LPJp_K0jrkmC~!N8`nOc+k*@{jK^4os1e`#>m8G@kwi)_y8;7L< zN70`>DLLR{ih9gv{LA$)+M^e5s_j7>?_YD($BZ#Q7*;dVHN^N>@JV*9KK4iZ#YWM6 z1e~#&c{02k^REq`50vQ=;FKljh_T61t z*l+#152+exCwo6T+pD2IjpR71Y{EQ!T~$uc1_GA+QN+(X)ZUt>Q>Y3}K*4QO@-COqSwl2bP$g=uSBagDiU7w3`n2)ZP z_oK<+`Y&$9r(j9+r%vVS(&A{xj@$o6i=aL1i!56#g!X-tF%T<&elIv-$3i~zPsuT* z_XzZ-lZyfgCINe{8hKpdK|eTe8sA34JZRIYseQPPeHHWgdOEImz3;8j6v6dyfY$CM zmvG(8dtX|k5ZBiSZN7Uf#PzqdS@_^ZT(4V}Cj@`Mb^aQrqA3&i1+=D(_wB=dhGI|C z^TDmSkFwG>y`+f+`DPzEtOi`)YlLQP`OSjZG(IQ08WtS8Xk!{(h2EN1p}+~-$ZBQK1!c0Jt30?H}-2dB)n#URlwME%q!ea2yKaZ^pXW!3~h_xWELcA z7^!!3j>Ef9~%*#)99DtB+hY#rcwYR7@Pg{h;ZiJMS9dJjH+cEUl0G zV2UdHBKNT1#pd2$GP=0mC1QB?{tg!O1g6jMY2$tcuk}3-E$qMN*EM}=EJ(3eE;iYS z?GAUIrPkqo)-zMqqBYpBcU6g=$}HH_mc65HIi7dHtVv}F?%yntb@Nfc{UUQi(Q;W9 zlyI9S$w}e9lXI~8DKQqzCK-N76vFWnp4s2VkNY>t{XgW{EI2XfxAzc(1>d`tWHp9$w*X??NmW`cxTz=TT& z6ZS{7$=_;WLgtCrY7ZNjAZ&klXIvc<8g)17Bvms(|Hm8E=jBXDvE3|@^qmRkW|wuv z7BL~Tx;;7k3!eAs&(VO7Oo(GTO*vCc&^j#YZ}FZ9U;l7;cD-eSi5Aml#cL*Td{+TI zl?i9nXXpMTGeOe0Md?i<6RK2qZS{?1f@}6&(&z~jzm&6Sh0G%+cql!ZXt>7&-+sl{ z(VC7RZEV3-8}z{!G|>dFh^SK1_JWpB&}t$%IQb*J_rzF=5YUq482@yl?Na zNqYw-z=V>PH(;RHAv0#F+h`;e!GbY?N4~!c&VZ!gO z^75+(nD}+d(w13$Cj5~Z{OP_2@7wNW^!QFD*c^4S|)Vd6o25Y%!Ikux#t^~F`=Y7QCfWwj>{pPxL`RZQ2h(&O_JE( z51!;|Q6|t;`HY+enQ(f!vStpO2{&>N&P>ypu=pC!JvADR@3vQ5_R|a~-uhD~af|_N zJgsI`{}^!G{={ysJ_fA$n9Ek`Vt|UD+mn5*418`zgl^lwfC{5It(&zB=v?$(H>#2W zvZMJ+(tj`@nif_DMGP=6dAh&oGXu=lUyCTqV}Q(wu(#zo43I846GX76I+SxwSch$+L#@1Q=|xRC;JSxNpkx}h|9);o zXfh8B9Qroh|=%W(g0R z{NU;9BgX^NVu5v!#d+Z6U{LWZK_0jm_g2`Z>uBmraK6RSa>% zF8Of3pFP}g!}z^WaT_<}*o>^t`@;=AJb!xAs=2}87a#fL2bLE|Ex%F74NHn2C!Eda zh64W8?{wdD!*sFLsl{)&!Nc}<+{6oR@OMyhDN5jm(d(U!;ZL~X`N`CJiwE2gJtlEs zSqL|{?|2{DeiN^MVV!QQ7dL#%5)Uo_2&nbzv$OoHi7^Z2FW7zPscl21jwht1h;A-UBWWCT~Xd z1#`ikl~t!CZgRn0x6o7#Pc9G=IqY%Ng$qpoek^yo#0ANsjxPeOxnRqCCJ}AU1@;T2 zx?Y+f@9o_C)_@BZ2d!F{qss+vZ=ZjZvy}_}C_nMZ+ROzCvL^>p*Kz@SiAjCz3NEOW zm)R4lfY+n#R=X<61y{d&KQ&)~?I)hRr_JI5sjL95h4VD1(YxO_GDZWB%q6CU{WMr% zTcI7^K?C){?awS5XrNRpbbCcL4R)0O?&&I}K|#d!Z;792u!+&O_(U!ZevF@(5YD7Q zk&W^0{FgNFy>MLIHh~7latdZ5ku*?zzG-2`T^c+(w!_ruHVtmyY~$~~P6PSycN!j7 zX^_P@*(~5d1AU8Yv`6P?@O$Tra0N>md^>adX{;#?=6!QlDHzfq<#MyqLtPrsCBvWa zZ>7OzE03G*n`mHQxw*4*H4XZfp8TM{j0W4PqnPP(G>CW^^jT1h2JZ*HPgwKOKx3Rc z<^wkkoDSV{k{+3->bP&+No$#>uF)0ml0WCEcE4AypJ1LkS!JoSHFciSJmNK!9yw19 zJlf|YdwZVxrA3Unxz1B&^wHfF*7MZrvV$Ir4$o8UC+EV>=*&|}zD_sd*3VOg@82J; z;>=S|!rU11g7ehtD^sBgvvbtMsO!Y0o;hkhMN4yM?Htu*e#3On*Ex#T{=i^Y<{V|b zE<-^*evZm9<`gczGe_~;%j@xY&QX6$|3%f=;`yKc6ek^>qc~reeZQbLM~PYqBr0v1 zqlC0;4C z@qoRVv(%{QjqcpIS?bxL3F#G~v(&}XHxf6m%~B_cnZ`foXQ}vPai^80v((dbJF713 zoTc_ll+8b1KTB1HO%2y8%u-Un*y#*`S?bNe_D##CW~k!X?(3SJGt}Kpx)1a!XQ;ES zk(c)6&rmxo*9q#T&QR|-d+gO7&rlYLs~&RvXQ*3#^4I4bXQ=24r`J?j&QN-ix{=rzsi|whH?$mr-}DZQ{C%=x{~XrDd%j% zpp9RrsjhR4eOYg&sp~PG>~%5Il=)Ye#+blq3XTX~V!2LJ+^qpwR;Q<_usN0NY@=x^ z^#8KH*Z!4_>NFLz*g#KOa+)fMF1(&UIYl{IHJq@jpQ7&Bg`yXjqLvmkOh1a9qU6lK zeOTp-$F~_OygxHV`P*DxvVHdy)mFd$TfxdHD(!x-(00KoYT~}%Y=4~->IvU)Y7khb9rV+N5yoht&`Nfvu7AKizX?@1nuc8 zo=K`?MON5+`ves^^sRVX;RIDrcYW&de1Zx*QLy-J-~^RD|3bU#;sj;dwU#DsG(nw} zGPKm%FhRA2@LQb_pP-ubPZMtAux3m%iOGPMMl^Z%A1-PIVbCYmZ`$Q-k+rHsAg`Mmfsct#T+HqihE} zKORUKqY7E8rq+gxQ3=+L6-@gvYDDmeevR=M^}B7^Pkwy3!bDF-q}!{U&s_k5X;_>M9zS zj8aQqIBPpIMyakBsrszW5lU+0`dnDy2$knEY#@zWJvB#ZJ0fn6P-d>lzr-&f8~X(U z^+%}6z_Rw4wIkFAZ+-K#!XuOrb@uA_p<$|No5JR`mBZ930f+myGKZ;p1JBrimCuPJ4D6UmR4VG7@~YHb8k<8 zA*w0ALjOn15S3)IZFa;PW!1g=#ZL`U_m2D{Rdj}^51%)b?_4%S9ds389cB(uW{yFB zEIR*DO>YwpohkfB)o7buI+Og5(iTZlw!HO^;;g)P^RUf7%J#OL->!ZCs8t)bEnl}9 z<;C1*B?SLbQb)}9j`a^xD!)w2%6<$|d6E5NiD`pWNiqkV?hI1<=K2-29R?}Yh$E_u z1B2AVYZlD~8wM%cBpww{u|ev}MnCPd-2+rt(iW?)F9s;X7pse9FAPvIT8EuZuNk0f z%(6}<4D?eId`5R`(vZ8Cu(%!isV0dv6lX&}mF-}ktvb?2T{+%2x+SxZB7Vg>YdQB( zF^$0*Dx3Q#A(xD8(i6Ru-4l_(@$6pep>eWwv1>0?_~1(SU5#F}SFbsv=^m=9IPVbS zLl2eqL47dAy@&GB>M7m3rH4}4Th-q%+fCJUjT>82-ISjAfvGl+ZYrs)qs?S%H${_u zBa%1&mr@zY@>!VomohaQ-E8OemohwHvo>4%FQr{8eSVbomy+tIb*;$nqMGh-k`27O zDF1>hDYo0XD8i3D>cQ1Tc~zHt`sH_0O_`PZeY`s<#ia|(owjvSdW){VG2`l_EK*kQ zRnPCBEDD8}^Lckr+74@U%C$SFiRV6v!L$w>2jM;1dF_|K)M_)It?iV` z62W4z`8LXbxrh5>s*NhlUT!bx(MEOEvXcF`v{4fQi!>pIjE`T`>!_IeT^BUo)>0}@ z_C!5gUrX7!!J+2V8Y<5v@vhj)8Y=D1jLG_Ezo?jR1#ucneo-n4bC*{~S5u*3NyWUf z)l}1xg5>fCRa9BIq2LXXDk?8pe{xZ9C1rRk_k1$Dk}3;(a#Hq21;u&S-{eB8pbTfz z+9|hkO5426ZgA`;Ro38to!{XnW#=9xBidI+m2F<8PCr{lH7(Nk0H1$QUJLj|_vL)2 zIO=)MDQTsYMbzQW*2EIZV(Sm9frw(NNmhB_M@SLnwZYo&hVM6ubJI^*$hDC2lJNh1 z>B3jatJ&sJg5?*gsdJrG>Y>jRO}S1pXm^UR8~ZObRh%krmik;+Gkc9+ZR|nJg2+{Ew)}$drHxUS@E5-u@vpoCC7a?9#b^0 zX10UCU5d7;h~G=`21QG%GdV+G7*Pj$Xju!4)}kH>EI-jJeR)Q&X2u z<^_dqUfcJPNY`Uphx?$ND(&C+A`$n~|B<J z{`V)okxiW@O?0bzw>=zDLw#*?9hUQoF^n+OYZu z8G{ezoOdcC4SPHNKhBqtv|Q=$Yj6D|2QSqHz7{JdISl2`d!oxpe8$0LLaBljs&3UQ|gAu9cj-aj1$^5y@K+-Bd+}R*y%N+EkMk*QjXW@oLhL zDTW`h{UQzRC&cRnYRELo=H`hT=lcGS~7H4KGRmKj*R(M%CB2k zN0xc>HEcTgn-r3(QPgVsP3Aq^(0ACjo(vsdxbo&iJ?W(u-csWGhYURyEvhKoKnhuI zq8>hJAY=5*TI7`*u|L9I@833(LgvqlGZ!_HQcp$2q|Y{ydM9+GZ>Bck@7rtro*!u< zYZlEM@zHE1Emr&77rx$1b|wD|dS2K}o(=EZuD_s#R2eJCm^su!PJBMsni<(b4l=91 zySBEFUD<3GJ>^!CsHfL2vTr5T?7q%2vRg@iC!_k2xmI$pC(fmRXB%l+k@I6Hu#K!) z8QI2F*+xDb-SJjdwjJMF^hQVfWIL%+|2O1PayzLw*%_WZ*iLe8rkM3_?jSAvG#a*e zc98yK!Ez73bdXKIo=eUObmDtMgMZi_>LhFW)O0!{J4w66FPv=KI!VRCwo@}JyU3=7 zcERw=U1ZmpRA;T6F0#g)eB8_Zmz?-^-~HL1zob-RhRlVqzof;#2&XM|f62mBy&$Q@ z-K5%|T+6w0-DHjX)KGtVH+dy!PP==uo77%Vw0Ce@4@qzz7v;LuLr$E|3YRNKS+llg zyG$=B<^8ex^2uJ(LhbUMSI>LNY>OUZY^awMno{1PtJz0p)7P(w_vs_`ylD08l0I_K z_Q5VUv3@eUX!ydoSw9*2NXpnFv7bz;r4KLYN8ZXde7R|W#OLB|4tNfbQrgbj8HEF+ z&SLEz3u2I@|4F)AuxgOhF?5PpZ8S)VxhKv)b{ZsYE!00Ph!`ZzW>wn#sX?;dIPf&9 zb&#yI*sF7!^^de$anW0B)ju-u>4;psA?6$n4HBE3{*hMSS6wQ5^pBKtX;^lNLb+yG zL$viDsWW7o`hhh>CKc@VvRO4m-ZQwBreHKg%H3VFwbp5fE$;BtnaE*$f2(_) zVg4{FR?m>u?ieQT2@Ged@{N!Q<2%yVsEm+m5fi*?4vmnuqe7q7UmYRE*g>^g(Ie!+ zDqH=11ta9rqiKsQ|Bm4A%$gi<6C5S2f=;$ZY#1d=G=KgC(^0Y@E?&O-+9;{>vG$vE z{3v-g;&O$4;V3zyY_A#6H%clxmUDd*86!uW!j}uGjghJ&WwFO(W5_ocU7_1`j1+PZuWn>aa6@{jScFZz#@ zTC3mw=}H+VN8}2G&Q+lNd1GAv1eP!73v-p9AlG_iFiNL%6N!02-mWc}mjaE*`& zvbSL;*Y6B0uM_y@RX0IK2CEOOq)p=Y>AASYc>aHHsQ;gL(Es<>i|fC?|NDLQ|L?i~ m``h^c)&HM=UR;LwHv9km-F(Bn4Ca5O-1z?rGREUS>Hh%!t&sHq literal 0 HcmV?d00001 diff --git a/tests/data/PM2000D_A_eff_marcuse.npz b/tests/data/PM2000D_A_eff_marcuse.npz new file mode 100644 index 0000000000000000000000000000000000000000..41b572e390483b701c5c7384ac04c342b9d37f92 GIT binary patch literal 2916 zcmd6pYcv$<8plVXi=8fpWS^0?wyDN(X(!tCJg!Nlve{&GBM~ZQkfPlMxl_ocXdCxX z$aM_U7&BvxWT=Qrc4Wj3yNgN5)XCZFoVE5@>wGz%&i`HS^RDN4{_oe{`rGZ4pGzeC zsZ|8oj-A;C3IxLJRVEM#r`&=(j(hkX4fN9Ur3K3omi%)w>z;MY>dJ1%U3OGCLJ;8$ z`G`k=yFd9yE%Fvm6SA%r+0)NI(BIAXke~k%kJ-7k+wlO8uk!#eH=4)SevP4#f$mye zt+QJHvlb?d*JMu#voY1KUbcL14*pujaWQ$BgN15*rt_&>EDLTS2G7gG+P*gWwe&pX zn0AU;q!|No*}E4 z`;?>4Azhl=bUAKvyE<&P^)}u~)=>Q^;Wk`;#i6gz8o zB$Za6eQbd&M_Pf88pyhpP>C+%#V)*vO6=Fvjd{>siHGgVn1iZSI2>v){PsW<`czdC z{jCwoTbToU6u=;0Cd!uo^#?fXqR;2E+RrDMy@YaGtV- zOJG(F#+eGDFAmq>V!AamMXwgG(In?fAq#Lb;wTYueILBLNmUr z#OX2%1(Sm+M|)UkqvK(5YGpm%+xwChb+jJsV!RpY<@MNA|Em4QR6RP6otR;oHees= z^*T{h1G?_y46|dR&-dbk9*5nR` zoIX$6bm0zuE|z0#7u~_Ir{V3?B^*>-q&Va4%)wZUUKf(hL1o%PYU~IHlXK#uG7K6~ z?V)mQ$>~PSn42_I&uzrTnK3aQbGdlG>TdRlHC()MPpR&tEf*anwjF0_TwMM+VK^d% zi;vYhWC?X#%+z0^n%2j~gEl|uW)YfD%PrZmaCH+720wHuvuVQoWTJoV@h0@V=@`RF zY{E5yv}|E@6aL=U!xZ;6;Vr3p$D^4h^i_5l?$_d>0Xb7PY{kR3rhUXeym(mcOwyf- z=V1{2pyij_JT!@AxXkP3;dJhE`r?l~WH+peSx(~Pxvs<6E6w?6HlD{Md+>3NqC~78 z!$*Op&aiPQA0v&A%1rO`(Z;TbXd~sLlJ_XZZiN8vg|4%71_55mIOgK&CP2HYTXe4o z0WP~Y9!)D2VE2$-_UU#3rpx@87v2bvs#YP6R&B;__1{VoH#Ot4Z9mAa9%x2}YalTv zv>DxkYbeD9&8VIH-lC$l86_pgF0Aopyum+BZ(1ru?><&^tC0{X(rLzhsu0H)ZDBq> zE5vf0v*LjqA<|7blF?=%{^0y&=FO-Or8G0*^r9AI#)MMltZ%{KT%N^3rxx7YAn&p) zqy=RUt>_wAEyx{M_Mu{CQ{zJz!lc+Wm#M?gLH-{>9q` zUwLbO$ivg~{}r{d6rq-5csP(AWSXCT5e}-AsXZI5FT#P&lHiP+7a``I!NpRQ2w-h+ z(>il50!;jt?OHq*0b=*Ok)ZvNprf>u-PjQcC!QNeeb^cW1CwUThK*6Mk_Kx+Hb=wL zyOUx`eKfRZz0!|@7}#>KgyYYTfoJ8mHXPen2)X1wKlE`dq~9`XeCQnq!(V(_J7sa; z#53HKk{S<6H^*LYHcWuX{48g-K3IC_3X~7-9ldk^3Y=ZFe&U^V3XG)f(r;``f#R|x{`M`Yu!Tc#;<8ggQ`Vbq3g#k%7VNj$p+E=S@7M? zFX2Ka1IF6o`A==L!FRPM`0%pfU8Kg8cuNi}%bHofo1Ftk6{#=u&2vGlP5;1c%Y~x~ zr_${A=K)#k?$D3#@}S79&a^u*A4ok5V#ACIpr|}Qym+7hf~``%XQmbcQM*lkrTq;! z>MX6^x4Z~kBkIilkQ70MCO^fUQw;vJrE)uSN?@*SSA*H*n=pOP=fr6AEnv6}3*KES zg@@*8sw;BKz=%cqOyZYAUvfvhckpeHX$zk%Tu}j5JnFs?ZUr1VN4$EDUJ0ciog06m zRKZfu+e#A?RS@{iH&!F;YIxO7HS15UfrP-0lcz#zVUdwzoX0UH;Hh8UAMvgO`ToEh zILU&tl%Nd$<$AcXTC=@})d0B{HVeHb+3+YSGc3T215x{^+vhSGA$(+YYi2YTUUen> zLb7dw>>!Ee4GkV7Zi`wwILU+k#kCOjln*aRJ@s2V1h8k>?fvgNn;}t58uT6z!uBsI z>@YbI6gUow4A!@TwC~DN|G>KtHhwnKp|uU1l^$$)YuFCw<1NmgDHlVC&138IRh>}s z^QSbusZLnTOLyP+_&%h@71kdUJ%AGga%x)hL)fqLyKn6ME-3tak3?s@8yLnjuFflZ zAgz$wsJ^!sx{Q;S(r!G4?452-CzYOpNP5U-&&59YwaG%?NTnb4hs53SWBmpmL+UfS z?C0<#@)EJ4XaFp?`}r#54T8_d#qXszU%)+nmWV1C0@Yns%_`#(sLSqMtwMeYaRguW zp94m~=dRTE`GZkdJ+RfH)A|*>-q~_YFjEjaX3YQI&Qgrf1OCu{Vi(EN+NPq*89*pSGWBK`d=VhQTq*FpDUqCIO^wHQU0utJf9c=TI!wsE(n^th2 zgT3qxarP(VzoJ%d4pIKUh1);l>ukLsH(UQ2b^r74UqL