From 4d424e52dd75d58ea6519606c5ce6e85e0d3c4b5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Sierro?= Date: Tue, 2 May 2023 12:32:24 +0200 Subject: [PATCH] Simplified setup, not working yet --- developement_help.md | 18 +- src/scgenerator/__init__.py | 50 +- src/scgenerator/evaluator.py | 18 +- src/scgenerator/math.py | 10 +- src/scgenerator/operators.py | 424 ++++-------- src/scgenerator/parameter.py | 8 +- src/scgenerator/physics/pulse.py | 6 +- src/scgenerator/solver.py | 624 +++++------------- src/scgenerator/utils.py | 17 +- tests/Optica_PM2000D/Optica_PM2000D.toml | 18 + .../PM2000D_2 extrapolated 4 0.npz | Bin 0 -> 32558 bytes .../Optica_PM2000D/PM2000D_A_eff_marcuse.npz | Bin 0 -> 2916 bytes tests/Optica_PM2000D/Pos30000New.npz | Bin 0 -> 49656 bytes tests/Optica_PM2000D/test_Optica_PM2000D.py | 27 + 14 files changed, 440 insertions(+), 780 deletions(-) mode change 100644 => 100755 developement_help.md mode change 100644 => 100755 src/scgenerator/operators.py create mode 100755 tests/Optica_PM2000D/Optica_PM2000D.toml create mode 100644 tests/Optica_PM2000D/PM2000D_2 extrapolated 4 0.npz create mode 100644 tests/Optica_PM2000D/PM2000D_A_eff_marcuse.npz create mode 100644 tests/Optica_PM2000D/Pos30000New.npz create mode 100644 tests/Optica_PM2000D/test_Optica_PM2000D.py diff --git a/developement_help.md b/developement_help.md old mode 100644 new mode 100755 index 5dbb76c..4e2bbd4 --- a/developement_help.md +++ b/developement_help.md @@ -3,9 +3,23 @@ - add it to README.md - add the necessary Rules in ```utils.parameters``` - optional : add a default value -- optional : add to valid_variable - optional : add to mandatory_parameters ## complicated Rule conditions - add the desired parameters to the init of the operator -- raise OperatorError if the conditions are not met \ No newline at end of file +- raise OperatorError if the conditions are not met + +## operators +There are 3 kinds of operators +- `SpecOperator(spec: np.ndarray, z: float) -> np.ndarray` : operate on the spectrum + - Envelope nonlinear operator used in the solver + - Full field nonlinear operator used in the solver +- `FieldOperator(field: np.ndarray, z: float) -> np.ndarray` : operate on the field + - SPM + - Raman + - Ionization +- `VariableQuantity(z: float) -> float | np.ndarray` : return the value of a certain quantity along the fiber depending on z + - dispersion + - refractive index + - full field nonlinear prefactor + - nonlinear parameter (chi3, n2, gamma) diff --git a/src/scgenerator/__init__.py b/src/scgenerator/__init__.py index 4831082..c4fb4c6 100644 --- a/src/scgenerator/__init__.py +++ b/src/scgenerator/__init__.py @@ -1,27 +1,27 @@ -# flake8: noqa -from scgenerator import math, operators -from scgenerator.evaluator import Evaluator -from scgenerator.helpers import * -from scgenerator.math import abs2, argclosest, normalized, span, tspace, wspace +# # flake8: noqa +# from scgenerator import math, operators +# from scgenerator.evaluator import Evaluator +# from scgenerator.helpers import * +# from scgenerator.math import abs2, argclosest, normalized, span, tspace, wspace from scgenerator.parameter import FileConfiguration, Parameters -from scgenerator.physics import fiber, materials, plasma, pulse, simulate, units -from scgenerator.physics.simulate import RK4IP, parallel_RK4IP, run_simulation -from scgenerator.physics.units import PlotRange -from scgenerator.plotting import ( - get_extent, - mean_values_plot, - plot_spectrogram, - propagation_plot, - single_position_plot, - transform_1D_values, - transform_2D_propagation, - transform_mean_values, -) -from scgenerator.spectra import SimulationSeries, Spectrum +# from scgenerator.physics import fiber, materials, plasma, pulse, simulate, units +# from scgenerator.physics.simulate import RK4IP, parallel_RK4IP, run_simulation +# from scgenerator.physics.units import PlotRange +# from scgenerator.plotting import ( +# get_extent, +# mean_values_plot, +# plot_spectrogram, +# propagation_plot, +# single_position_plot, +# transform_1D_values, +# transform_2D_propagation, +# transform_mean_values, +# ) +# from scgenerator.spectra import SimulationSeries, Spectrum from scgenerator.utils import Paths, _open_config, open_single_config, simulations_list -from scgenerator.variationer import ( - DescriptorDict, - VariationDescriptor, - Variationer, - VariationSpecsError, -) +# from scgenerator.variationer import ( +# DescriptorDict, +# VariationDescriptor, +# Variationer, +# VariationSpecsError, +# ) diff --git a/src/scgenerator/evaluator.py b/src/scgenerator/evaluator.py index 2b3f216..7210285 100644 --- a/src/scgenerator/evaluator.py +++ b/src/scgenerator/evaluator.py @@ -8,7 +8,7 @@ import numpy as np from scgenerator import math, operators, utils from scgenerator.const import MANDATORY_PARAMETERS from scgenerator.errors import EvaluatorError, NoDefaultError -from scgenerator.physics import fiber, materials, pulse, units, plasma +from scgenerator.physics import fiber, materials, plasma, pulse, units from scgenerator.utils import _mock_function, func_rewrite, get_arg_names, get_logger @@ -412,8 +412,12 @@ default_rules: list[Rule] = [ Rule("n_eff_op", operators.vincetti_refractive_index), Rule("gas_op", operators.ConstantGas), Rule("gas_op", operators.PressureGradientGas), - Rule("loss_op", operators.constant_array_operator, ["alpha_arr"]), - Rule("loss_op", operators.no_op_freq, priorities=-1), + Rule("square_index", lambda gas_op: gas_op.square_index), + Rule("number_density", lambda gas_op: gas_op.number_density), + Rule("n2_op", lambda gas_op: gas_op.n2), + Rule("chi3_op", lambda gas_op: gas_op.chi3), + Rule("loss_op", operators.constant_quantity, ["alpha_arr"]), + Rule("loss_op", lambda: operators.constant_quantity(0), priorities=-1), ] envelope_rules = default_rules + [ @@ -435,12 +439,10 @@ envelope_rules = default_rules + [ ), # Operators Rule("gamma_op", operators.variable_gamma, priorities=2), - Rule("gamma_op", operators.constant_array_operator, ["gamma_arr"], priorities=1), - Rule( - "gamma_op", lambda w_num, gamma: operators.constant_array_operator(np.ones(w_num) * gamma) - ), + Rule("gamma_op", operators.constant_quantity, ["gamma_arr"], priorities=1), + Rule("gamma_op", lambda w_num, gamma: operators.constant_quantity(np.ones(w_num) * gamma)), Rule("gamma_op", operators.no_op_freq, priorities=-1), - Rule("ss_op", lambda w_c, w0: operators.constant_array_operator(w_c / w0)), + Rule("ss_op", lambda w_c, w0: operators.constant_quantity(w_c / w0)), Rule("ss_op", operators.no_op_freq, priorities=-1), Rule("spm_op", operators.envelope_spm), Rule("spm_op", operators.no_op_freq, priorities=-1), diff --git a/src/scgenerator/math.py b/src/scgenerator/math.py index bd451ef..7d2a569 100644 --- a/src/scgenerator/math.py +++ b/src/scgenerator/math.py @@ -26,7 +26,7 @@ def span(*vec: np.ndarray) -> tuple[float, float]: raise ValueError(f"did not provide any value to span") for x in vec: x = np.atleast_1d(x) - out = (np.min([np.min(x), out[0]]), np.max([np.max(x), out[1]])) + out = (min(np.min(x), out[0]), max(np.max(x), out[1])) return out @@ -193,7 +193,7 @@ def tspace(time_window=None, t_num=None, dt=None): def build_envelope_w_grid(t: np.ndarray, w0: float) -> tuple[np.ndarray, np.ndarray, np.ndarray]: """ - convenience function to + convenience function to Parameters ---------- @@ -289,8 +289,8 @@ def polynom_extrapolation(x: np.ndarray, y: np.ndarray, degree: int) -> np.ndarr y array with zero values on either side replaces with polynomial extrapolation Example - - + + """ out = y.copy() @@ -348,7 +348,7 @@ def envelope_ind( signal: np.ndarray, dmin: int = 1, dmax: int = 1, split: bool = False ) -> tuple[np.ndarray, np.ndarray]: """ - Attempts to separate the envolope from a periodic signal and return the indices of the + Attempts to separate the envolope from a periodic signal and return the indices of the top/bottom envelope of a signal Parameters diff --git a/src/scgenerator/operators.py b/src/scgenerator/operators.py old mode 100644 new mode 100755 index d05e5be..9ea0e1f --- a/src/scgenerator/operators.py +++ b/src/scgenerator/operators.py @@ -4,189 +4,56 @@ Nothing except the solver should depend on this file """ from __future__ import annotations -from abc import abstractmethod -from copy import deepcopy -from dataclasses import dataclass -from functools import cached_property -from typing import Any, Callable, Protocol +from typing import Callable import numpy as np -from matplotlib.cbook import operator - from scgenerator import math from scgenerator.logger import get_logger from scgenerator.physics import fiber, materials, plasma, pulse, units - -class SimulationState: - length: float - z: float - current_step_size: float - conversion_factor: np.ndarray | float - converter: Callable[[np.ndarray], np.ndarray] - stats: dict[str, Any] - spectrum: np.ndarray - spec2: np.ndarray - field: np.ndarray - field2: np.ndarray - - def __init__( - self, - spectrum: np.ndarray, - length: float = 10.0, - current_step_size: float = 1.0, - z: float = 0.0, - conversion_factor: np.ndarray | float = 1.0, - converter: Callable[[np.ndarray], np.ndarray] = np.fft.ifft, - spectrum2: np.ndarray | None = None, - field: np.ndarray | None = None, - field2: np.ndarray | None = None, - stats: dict[str, Any] | None = None, - ): - self.length = length - self.z = z - self.current_step_size = current_step_size - self.conversion_factor = conversion_factor - self.converter = converter - - if spectrum2 is None and field is None and field2 is None: - self.set_spectrum(spectrum) - elif any(el is None for el in (spectrum2, field, field2)): - raise ValueError( - "You must provide either all three of (spectrum2, field, field2) or none of them" - ) - else: - self.spectrum = spectrum - self.spectrum2 = spectrum2 - self.field = field - self.field2 = field2 - self.stats = stats or {} - - @property - def z_ratio(self) -> float: - return self.z / self.length - - @property - def actual_spectrum(self) -> np.ndarray: - return self.conversion_factor * self.spectrum - - @cached_property - def spectrum2(self) -> np.ndarray: - return math.abs2(self.spectrum) - - @cached_property - def field(self) -> np.ndarray: - return self.converter(self.spectrum) - - @cached_property - def field2(self) -> np.ndarray: - return math.abs2(self.field) - - def set_spectrum(self, new_spectrum: np.ndarray)->SimulationState: - """sets the new spectrum and clears cached properties""" - self.spectrum = new_spectrum - for el in ["spectrum2", "field", "field2"]: - if el in self.__dict__: - delattr(self, el) - return self - - def clear(self): - """clears cached properties and stats dict""" - self.stats = {} - for el in ["spectrum2", "field", "field2"]: - if el in self.__dict__: - delattr(self, el) - - def copy(self) -> SimulationState: - return SimulationState( - self.spectrum.copy(), - self.length, - self.current_step_size, - self.z, - self.conversion_factor, - self.converter, - self.spectrum2.copy(), - self.field.copy(), - self.field2.copy(), - deepcopy(self.stats), - ) +SpecOperator = Callable[[np.ndarray, float], np.ndarray] +FieldOperator = Callable[[np.ndarray, float], np.ndarray] +VariableQuantity = Callable[[float], float | np.ndarray] +Qualifier = Callable[[np.ndarray, float, float], float] -Operator = Callable[[SimulationState], np.ndarray] -Qualifier = Callable[[SimulationState], float] - - -def no_op_time(t_num) -> Operator: +def no_op_time(t_num) -> SpecOperator: arr_const = np.zeros(t_num) - def operate(state: SimulationState) -> np.ndarray: + def operate(spec: np.ndarray, z: float) -> np.ndarray: return arr_const return operate -def no_op_freq(w_num) -> Operator: +def no_op_freq(w_num) -> SpecOperator: arr_const = np.zeros(w_num) - def operate(state: SimulationState) -> np.ndarray: + def operate(spec: np.ndarray, z: float) -> np.ndarray: return arr_const return operate -def constant_array_operator(array: np.ndarray) -> Operator: - def operate(state: SimulationState) -> np.ndarray: +def constant_array_operator(array: np.ndarray) -> SpecOperator: + def operate(spec: np.ndarray, z: float) -> np.ndarray: return array return operate +def constant_quantity(qty: float | np.ndarray) -> VariableQuantity: + def quantity(z: float) -> float | np.ndarray: + return qty + + return quantity + + ################################################## ###################### GAS ####################### ################################################## -class GasOp(Protocol): - def pressure(self, state: SimulationState) -> float: - """returns the pressure at the current - - Parameters - ---------- - state : CurrentState - - Returns - ------- - float - pressure un bar - """ - - def number_density(self, state: SimulationState) -> float: - """returns the number density in 1/m^3 of at the current state - - Parameters - ---------- - state : CurrentState - - Returns - ------- - float - number density in 1/m^3 - """ - - def square_index(self, state: SimulationState) -> np.ndarray: - """returns the square of the material refractive index at the current state - - Parameters - ---------- - state : CurrentState - - Returns - ------- - np.ndarray - n^2 - """ - - class ConstantGas: gas: materials.Gas pressure_const: float @@ -205,16 +72,21 @@ class ConstantGas: self.pressure_const = pressure self.number_density_const = self.gas.number_density(temperature, pressure, ideal_gas) self.n_gas_2_const = self.gas.sellmeier.n_gas_2(wl_for_disp, temperature, pressure) + self.n2_const = self.gas.n2(temperature, pressure) + self.chi3_const = self.gas.chi3(temperature, pressure) - def pressure(self, state: SimulationState = None) -> float: - return self.pressure_const - - def number_density(self, state: SimulationState = None) -> float: + def number_density(self, z: float) -> float: return self.number_density_const - def square_index(self, state: SimulationState = None) -> float: + def square_index(self, z: float) -> np.ndarray: return self.n_gas_2_const + def n2(self, z: float) -> np.ndarray: + return self.n2_const + + def chi3(self, z: float) -> np.ndarray: + return self.chi3_const + class PressureGradientGas: gas: materials.Gas @@ -230,6 +102,7 @@ class PressureGradientGas: temperature: float, ideal_gas: bool, wl_for_disp: np.ndarray, + length: float, ): self.p_in = pressure_in self.p_out = pressure_out @@ -237,15 +110,22 @@ class PressureGradientGas: self.temperature = temperature self.ideal_gas = ideal_gas self.wl_for_disp = wl_for_disp + self.z_max = length - def pressure(self, state: SimulationState) -> float: - return materials.pressure_from_gradient(state.z_ratio, self.p_in, self.p_out) + def _pressure(self, z: float) -> float: + return materials.pressure_from_gradient(z / self.z_max, self.p_in, self.p_out) - def number_density(self, state: SimulationState) -> float: - return self.gas.number_density(self.temperature, self.pressure(state), self.ideal_gas) + def number_density(self, z: float) -> float: + return self.gas.number_density(self.temperature, self._pressure(z), self.ideal_gas) - def square_index(self, state: SimulationState) -> np.ndarray: - return self.gas.sellmeier.n_gas_2(self.wl_for_disp, self.temperature, self.pressure(state)) + def square_index(self, z: float) -> np.ndarray: + return self.gas.sellmeier.n_gas_2(self.wl_for_disp, self.temperature, self._pressure(z)) + + def n2(self, z: float) -> np.ndarray: + return self.gas.n2(self.temperature, self._pressure(z)) + + def chi3(self, z: float) -> np.ndarray: + return self.gas.chi3(self.temperature, self._pressure(z)) ################################################## @@ -253,33 +133,26 @@ class PressureGradientGas: ################################################## -def constant_refractive_index(n_eff: np.ndarray) -> Operator: - def operate(state: SimulationState) -> np.ndarray: - return n_eff - - return operate - - def marcatili_refractive_index( - gas_op: GasOp, core_radius: float, wl_for_disp: np.ndarray -) -> Operator: - def operate(state: SimulationState) -> np.ndarray: - return fiber.n_eff_marcatili(wl_for_disp, gas_op.square_index(state), core_radius) + square_index: VariableQuantity, core_radius: float, wl_for_disp: np.ndarray +) -> VariableQuantity: + def operate(z: float) -> np.ndarray: + return fiber.n_eff_marcatili(wl_for_disp, square_index(z), core_radius) return operate def marcatili_adjusted_refractive_index( - gas_op: GasOp, core_radius: float, wl_for_disp: np.ndarray -) -> Operator: - def operate(state: SimulationState) -> np.ndarray: - return fiber.n_eff_marcatili_adjusted(wl_for_disp, gas_op.square_index(state), core_radius) + square_index: VariableQuantity, core_radius: float, wl_for_disp: np.ndarray +) -> VariableQuantity: + def operate(z: float) -> np.ndarray: + return fiber.n_eff_marcatili_adjusted(wl_for_disp, square_index(z), core_radius) return operate def vincetti_refractive_index( - gas_op: GasOp, + square_index: VariableQuantity, core_radius: float, wl_for_disp: np.ndarray, wavelength: float, @@ -288,12 +161,12 @@ def vincetti_refractive_index( gap: float, n_tubes: int, n_terms: int, -) -> Operator: - def operate(state: SimulationState) -> np.ndarray: +) -> VariableQuantity: + def operate(z: float) -> np.ndarray: return fiber.n_eff_vincetti( wl_for_disp, wavelength, - gas_op.square_index(state), + square_index(z), wall_thickness, tube_radius, gap, @@ -312,7 +185,7 @@ def vincetti_refractive_index( def constant_polynomial_dispersion( beta2_coefficients: np.ndarray, w_c: np.ndarray, -) -> Operator: +) -> VariableQuantity: """ dispersion approximated by fitting a polynome on the dispersion and evaluating on the envelope @@ -322,10 +195,7 @@ def constant_polynomial_dispersion( ) disp_arr = fiber.fast_poly_dispersion_op(w_c, beta2_coefficients, w_power_fact) - def operate(state: SimulationState) -> np.ndarray: - return disp_arr - - return operate + return constant_quantity(disp_arr) def constant_direct_dispersion( @@ -335,7 +205,7 @@ def constant_direct_dispersion( n_eff: np.ndarray, dispersion_ind: np.ndarray, w_order: np.ndarray, -) -> Operator: +) -> VariableQuantity: """ Direct dispersion for when the refractive index is known """ @@ -347,25 +217,22 @@ def constant_direct_dispersion( disp_arr[w_order], left_ind, right_ind, 1 ) - def operate(state: SimulationState) -> np.ndarray: - return disp_arr - - return operate + return constant_quantity(disp_arr) def direct_dispersion( w_for_disp: np.ndarray, w0: np.ndarray, t_num: int, - n_eff_op: Operator, + n_eff_op: SpecOperator, dispersion_ind: np.ndarray, -) -> Operator: +) -> VariableQuantity: disp_arr = np.zeros(t_num, dtype=complex) w0_ind = math.argclosest(w_for_disp, w0) - def operate(state: SimulationState) -> np.ndarray: + def operate(z: float) -> np.ndarray: disp_arr[dispersion_ind] = fiber.fast_direct_dispersion( - w_for_disp, w0, n_eff_op(state), w0_ind + w_for_disp, w0, n_eff_op(z), w0_ind )[2:-2] return disp_arr @@ -391,10 +258,7 @@ def constant_wave_vector( beta_arr[w_order], left_ind, right_ind, 1.0 ) - def operate(state: SimulationState) -> np.ndarray: - return beta_arr - - return operate + return constant_quantity(beta_arr) ################################################## @@ -402,26 +266,23 @@ def constant_wave_vector( ################################################## -def envelope_raman(raman_type: str, raman_fraction: float, t: np.ndarray) -> Operator: +def envelope_raman(raman_type: str, raman_fraction: float, t: np.ndarray) -> FieldOperator: hr_w = fiber.delayed_raman_w(t, raman_type) - def operate(state: SimulationState) -> np.ndarray: - return raman_fraction * np.fft.ifft(hr_w * np.fft.fft(state.field2)) + def operate(field: np.ndarray, z: float) -> np.ndarray: + return raman_fraction * np.fft.ifft(hr_w * np.fft.fft(math.abs2(field))) return operate def full_field_raman( raman_type: str, raman_fraction: float, t: np.ndarray, w: np.ndarray, chi3: float -) -> Operator: +) -> FieldOperator: hr_w = fiber.delayed_raman_w(t, raman_type) factor_in = units.epsilon0 * chi3 * raman_fraction - factor_out = 1j * w**2 / (2.0 * units.c**2 * units.epsilon0) - def operate(state: SimulationState) -> np.ndarray: - return factor_out * np.fft.rfft( - factor_in * state.field * np.fft.irfft(hr_w * np.fft.rfft(state.field2)) - ) + def operate(field: np.ndarray, z: float) -> np.ndarray: + return factor_in * field * np.fft.irfft(hr_w * np.fft.rfft(math.abs2(field))) return operate @@ -431,22 +292,21 @@ def full_field_raman( ################################################## -def envelope_spm(raman_fraction: float) -> Operator: +def envelope_spm(raman_fraction: float) -> FieldOperator: fraction = 1 - raman_fraction - def operate(state: SimulationState) -> np.ndarray: - return fraction * state.field2 + def operate(field: np.ndarray, z: float) -> np.ndarray: + return fraction * math.abs2(field) return operate -def full_field_spm(raman_fraction: float, w: np.ndarray, chi3: float) -> Operator: +def full_field_spm(raman_fraction: float, w: np.ndarray, chi3: float) -> FieldOperator: fraction = 1 - raman_fraction - factor_out = 1j * w**2 / (2.0 * units.c**2 * units.epsilon0) factor_in = fraction * chi3 * units.epsilon0 - def operate(state: SimulationState) -> np.ndarray: - return factor_out * np.fft.rfft(factor_in * state.field2 * state.field) + def operate(field: np.ndarray, z: float) -> np.ndarray: + return factor_in * field**3 return operate @@ -456,12 +316,11 @@ def full_field_spm(raman_fraction: float, w: np.ndarray, chi3: float) -> Operato ################################################## -def variable_gamma(gas_op: GasOp, w0: float, A_eff: float, t_num: int) -> Operator: +def variable_gamma(n2_op: VariableQuantity, w0: float, A_eff: float, t_num: int) -> SpecOperator: arr = np.ones(t_num) - def operate(state: SimulationState) -> np.ndarray: - n2 = gas_op.square_index(state) - return arr * fiber.gamma_parameter(n2, w0, A_eff) + def operate(z: float) -> np.ndarray: + return arr * fiber.gamma_parameter(n2_op(z), w0, A_eff) return operate @@ -471,15 +330,15 @@ def variable_gamma(gas_op: GasOp, w0: float, A_eff: float, t_num: int) -> Operat ################################################## -def ionization(w: np.ndarray, gas_op: GasOp, plasma_obj: plasma.Plasma) -> Operator: - factor_out = -w / (2.0 * units.c**2 * units.epsilon0) - - def operate(state: SimulationState) -> np.ndarray: - N0 = gas_op.number_density(state) - plasma_info = plasma_obj(state.field, N0) - state.stats["ionization_fraction"] = plasma_info.electron_density[-1] / N0 - state.stats["electron_density"] = plasma_info.electron_density[-1] - return factor_out * np.fft.rfft(plasma_info.polarization) +def ionization( + w: np.ndarray, number_density: VariableQuantity, plasma_obj: plasma.Plasma +) -> FieldOperator: + def operate(field: np.ndarray, z: float) -> np.ndarray: + N0 = number_density(z) + plasma_info = plasma_obj(field, N0) + # state.stats["ionization_fraction"] = plasma_info.electron_density[-1] / N0 + # state.stats["electron_density"] = plasma_info.electron_density[-1] + return plasma_info.polarization return operate @@ -489,51 +348,41 @@ def ionization(w: np.ndarray, gas_op: GasOp, plasma_obj: plasma.Plasma) -> Opera ################################################## -def photon_number_with_loss(w: np.ndarray, gamma_op: Operator, loss_op: Operator) -> Qualifier: +def photon_number_with_loss( + w: np.ndarray, gamma: VariableQuantity, loss: VariableQuantity +) -> Qualifier: w = w dw = w[1] - w[0] - def qualify(state: SimulationState) -> float: - return pulse.photon_number_with_loss( - state.spectrum2, - w, - dw, - gamma_op(state), - loss_op(state), - state.current_step_size, - ) + def qualify(spec: np.ndarray, z: float, h: float) -> float: + return pulse.photon_number_with_loss(math.abs2(spec), w, dw, gamma(z), loss(z), h) return qualify -def photon_number_without_loss(w: np.ndarray, gamma_op: Operator) -> Qualifier: +def photon_number_without_loss(w: np.ndarray, gamma: VariableQuantity) -> Qualifier: dw = w[1] - w[0] - def qualify(state: SimulationState) -> float: - return pulse.photon_number(state.spectrum2, w, dw, gamma_op(state)) + def qualify(spec: np.ndarray, z: float, h: float) -> float: + return pulse.photon_number(math.abs2(spec), w, dw, gamma(z)) return qualify -def energy_with_loss(w: np.ndarray, loss_op: Operator) -> Qualifier: +def energy_with_loss(w: np.ndarray, loss: SpecOperator, conversion_factor: float) -> Qualifier: dw = w[1] - w[0] - def qualify(state: SimulationState) -> float: - return pulse.pulse_energy_with_loss( - math.abs2(state.conversion_factor * state.spectrum), - dw, - loss_op(state), - state.current_step_size, - ) + def qualify(spec: np.ndarray, z: float, h: float) -> float: + return pulse.pulse_energy_with_loss(math.abs2(conversion_factor * spec), dw, loss(z), h) return qualify -def energy_without_loss(w: np.ndarray) -> Qualifier: +def energy_without_loss(w: np.ndarray, conversion_factor: float) -> Qualifier: dw = w[1] - w[0] - def qualify(state: SimulationState) -> float: - return pulse.pulse_energy(math.abs2(state.conversion_factor * state.spectrum), dw) + def qualify(spec: np.ndarray, z: float, h: float) -> float: + return pulse.pulse_energy(math.abs2(conversion_factor * spec), dw) return qualify @@ -542,25 +391,26 @@ def conserved_quantity( adapt_step_size: bool, raman: bool, loss: bool, - gamma_op: Operator, - loss_op: Operator, + gamma: VariableQuantity, + loss_op: VariableQuantity, w: np.ndarray, + conversion_factor: float, ) -> Qualifier: if not adapt_step_size: return lambda state: 0.0 logger = get_logger(__name__) if raman and loss: logger.debug("Conserved quantity : photon number with loss") - return photon_number_with_loss(w, gamma_op, loss_op) + return photon_number_with_loss(w, gamma, loss_op) elif raman: logger.debug("Conserved quantity : photon number without loss") - return photon_number_without_loss(w, gamma_op) + return photon_number_without_loss(w, gamma) elif loss: logger.debug("Conserved quantity : energy with loss") - return energy_with_loss(w, loss_op) + return energy_with_loss(w, loss_op, conversion_factor) else: logger.debug("Conserved quantity : energy without loss") - return energy_without_loss(w) + return energy_without_loss(w, conversion_factor) ################################################## @@ -568,23 +418,32 @@ def conserved_quantity( ################################################## -def envelope_linear_operator(dispersion_op: Operator, loss_op: Operator) -> Operator: - def operate(state: SimulationState) -> np.ndarray: - return dispersion_op(state) - loss_op(state) / 2 +def envelope_linear_operator( + dispersion_op: VariableQuantity, loss_op: VariableQuantity +) -> VariableQuantity: + def operate(z: float) -> np.ndarray: + return dispersion_op(z) - loss_op(z) / 2 return operate def full_field_linear_operator( - beta_op: Operator, - loss_op: Operator, + beta_op: VariableQuantity, + loss_op: VariableQuantity, frame_velocity: float, w: np.ndarray, -) -> operator: +) -> VariableQuantity: delay = w / frame_velocity - def operate(state: SimulationState) -> np.ndarray: - return 1j * (beta_op(state) - delay) - loss_op(state) / 2 + def operate(z: float) -> np.ndarray: + return 1j * (beta_op(z) - delay) - loss_op(z) / 2 + + return operate + + +def fullfield_nl_prefactor(w: np.ndarray, n_eff: VariableQuantity) -> VariableQuantity: + def operate(z: float) -> np.ndarray: + return w / (2 * units.c * units.epsilon0 * n_eff(z)) return operate @@ -595,14 +454,18 @@ def full_field_linear_operator( def envelope_nonlinear_operator( - gamma_op: Operator, ss_op: Operator, spm_op: Operator, raman_op: Operator -) -> Operator: - def operate(state: SimulationState) -> np.ndarray: + gamma_op: VariableQuantity, + ss_op: VariableQuantity, + spm_op: FieldOperator, + raman_op: FieldOperator, +) -> SpecOperator: + def operate(spec: np.ndarray, z: float) -> np.ndarray: + field = np.fft.ifft(spec) return ( -1j - * gamma_op(state) - * (1 + ss_op(state)) - * np.fft.fft(state.field * (spm_op(state) + raman_op(state))) + * gamma_op(z) + * (1 + ss_op(z)) + * np.fft.fft(field * (spm_op(field, z) + raman_op(field, z))) ) return operate @@ -610,13 +473,14 @@ def envelope_nonlinear_operator( def full_field_nonlinear_operator( w: np.ndarray, - raman_op: Operator, - spm_op: Operator, - plasma_op: Operator, - beta_op: Operator, -) -> Operator: - def operate(state: SimulationState) -> np.ndarray: - total_nonlinear = spm_op(state) + raman_op(state) + plasma_op(state) - return total_nonlinear / beta_op(state) + raman_op: FieldOperator, + spm_op: FieldOperator, + plasma_op: FieldOperator, + fullfield_nl_prefactor: VariableQuantity, +) -> SpecOperator: + def operate(spec: np.ndarray, z: float) -> np.ndarray: + field = np.fft.irfft(spec) + total_nonlinear = spm_op(field) + raman_op(field) + plasma_op(field) + return 1j * fullfield_nl_prefactor(z) * np.fft.rfft(total_nonlinear) return operate diff --git a/src/scgenerator/parameter.py b/src/scgenerator/parameter.py index 8c0d690..9bfcf3b 100644 --- a/src/scgenerator/parameter.py +++ b/src/scgenerator/parameter.py @@ -18,8 +18,7 @@ from scgenerator.const import MANDATORY_PARAMETERS, PARAM_FN, VALID_VARIABLE, __ from scgenerator.errors import EvaluatorError from scgenerator.evaluator import Evaluator from scgenerator.logger import get_logger -from scgenerator.operators import Qualifier -from scgenerator.solver import Integrator +from scgenerator.operators import Qualifier, SpecOperator from scgenerator.utils import DebugDict, fiber_folder, update_path_name from scgenerator.variationer import VariationDescriptor, Variationer @@ -393,9 +392,8 @@ class Parameters: worker_num: int = Parameter(positive(int)) # computed - linear_operator: Operator = Parameter(is_function) - nonlinear_operator: Operator = Parameter(is_function) - integrator: Integrator = Parameter(type_checker(Integrator)) + linear_operator: SpecOperator = Parameter(is_function) + nonlinear_operator: SpecOperator = Parameter(is_function) conserved_quantity: Qualifier = Parameter(is_function) fft: Callable[[np.ndarray], np.ndarray] = Parameter(is_function) ifft: Callable[[np.ndarray], np.ndarray] = Parameter(is_function) diff --git a/src/scgenerator/physics/pulse.py b/src/scgenerator/physics/pulse.py index 6cfe1dc..a5d3fad 100644 --- a/src/scgenerator/physics/pulse.py +++ b/src/scgenerator/physics/pulse.py @@ -1107,7 +1107,7 @@ def measure_properties(spectra, t, compress=True, return_limits=False, debug="") # Compute mean coherence mean_g12 = avg_g12(spectra) - fwhm_abs = math.length(fwhm_lim) + fwhm_abs = math.span(fwhm_lim) # To compute amplitude and fwhm fluctuations, we need to measure every single peak P0 = [] @@ -1119,7 +1119,7 @@ def measure_properties(spectra, t, compress=True, return_limits=False, debug="") lobe_lim, fwhm_lim, _, big_spline = find_lobe_limits(t, f, debug) all_lims.append((lobe_lim, fwhm_lim)) P0.append(big_spline(lobe_lim[2])) - fwhm.append(math.length(fwhm_lim)) + fwhm.append(math.span(fwhm_lim)) t_offset.append(lobe_lim[2]) energies.append(np.trapz(fields, t)) fwhm_var = np.std(fwhm) / np.mean(fwhm) @@ -1165,7 +1165,7 @@ def measure_field(t: np.ndarray, field: np.ndarray) -> Tuple[float, float, float else: intensity = field _, fwhm_lim, _, _ = find_lobe_limits(t, intensity) - fwhm = math.length(fwhm_lim) + fwhm = math.span(fwhm_lim) peak_power = intensity.max() energy = np.trapz(intensity, t) return fwhm, peak_power, energy diff --git a/src/scgenerator/solver.py b/src/scgenerator/solver.py index f321fab..d04f119 100644 --- a/src/scgenerator/solver.py +++ b/src/scgenerator/solver.py @@ -1,415 +1,34 @@ from __future__ import annotations import warnings -from typing import Any, Callable, Iterator, Protocol, Type +from collections import defaultdict +from typing import Any, Iterator, Sequence import numba import numpy as np -from scgenerator import math -from scgenerator.logger import get_logger -from scgenerator.operators import Operator, Qualifier, SimulationState -from scgenerator.utils import get_arg_names - -Integrator = None +from scgenerator.math import abs2 +from scgenerator.operators import SpecOperator +from scgenerator.utils import TimedMessage -class Stepper(Protocol): - def set_state(self, state: SimulationState): - """set the initial state of the stepper""" +class SimulationResult: + spectra: np.ndarray + size: int + stats: dict[str, list[Any]] + z: np.ndarray - def take_step( - self, state: SimulationState, step_size: float, new_step: bool - ) -> tuple[SimulationState, float]: - """ - Paramters - --------- - state : SimulationState - state of the simulation at the start of the step - step_size : float - step size - new_step : bool - whether it's the first time this particular step is attempted - - Returns - ------- - SimulationState - final state at the end of the step - float - estimated numerical error - """ - - -class StepJudge(Protocol): - def __call__(self, error: float, step_size: float) -> tuple[float, bool]: - """ - Parameters - ---------- - error : float - relative error - step_size : float - step size that lead to `error` - - Returns - ------- - float - new step size - bool - the given `error` was acceptable and the step should be accepted - """ - - -def no_judge(error: float, step_size: float) -> tuple[float, bool]: - """always accept the step and keep the same step size""" - return step_size, True - - -def adaptive_judge( - target_error: float, order: int, step_size_change_bounds: tuple[float, float] = (0.5, 2.0) -) -> Callable[[float, float], tuple[float, bool]]: - """ - smoothly adapt the step size to stay within a range of tolerated error - - Parameters - ---------- - target_error : float - desirable relative local error - order : float - order of the integration method - step_size_change_bounds : tuple[float, float], optional - lower and upper limit determining how fast the step size may change. By default, the new - step size it at least half the previous one and at most double. - """ - exponent = 1 / order - smin, smax = step_size_change_bounds - - def judge_step(error: float, step_size: float) -> tuple[float, bool]: - if error > 0: - next_h_factor = max(smin, min(smax, (target_error / error) ** exponent)) + def __init__(self, spectra: Sequence[np.ndarray], stats: dict[str, list[Any]]): + if "z" in stats: + self.z = np.array(stats["z"]) else: - next_h_factor = 2.0 - h_next_step = step_size * next_h_factor - accepted = error <= 2 * target_error - return h_next_step, accepted + self.z = np.arange(len(spectra), dtype=float) + self.size = len(self.z) + self.spectra = spectra + self.stats = stats - return judge_step - - -def decide_step_alt(self, h: float) -> tuple[float, bool]: - """decides if the current step must be accepted and computes the next step - size regardless - - Parameters - ---------- - h : float - size of the step used to set the current self.current_error - - Returns - ------- - float - next step size - bool - True if the step must be accepted - """ - error = self.current_error / self.target_error - if error > 2: - accepted = False - next_h_factor = 0.5 - elif 1 < error <= 2: - accepted = True - next_h_factor = 2 ** (-1 / self.order) - elif 0.1 < error <= 1: - accepted = True - next_h_factor = 1.0 - else: - accepted = True - next_h_factor = 2 ** (1 / self.order) - h_next_step = h * next_h_factor - if not accepted: - self.steps_rejected += 1 - self.logger.info( - f"step {self.state.step} rejected : {h=}, {self.current_error=}, {h_next_step=}" - ) - return h_next_step, accepted - - -class ConservedQuantityIntegrator: - last_qty: float - conserved_quantity: Qualifier - current_error: float = 0.0 - - def __init__( - self, - init_state: SimulationState, - linear_operator: Operator, - nonlinear_operator: Operator, - tolerated_error: float, - conserved_quantity: Qualifier, - ): - super().__init__(init_state, linear_operator, nonlinear_operator, tolerated_error) - self.conserved_quantity = conserved_quantity - self.last_qty = self.conserved_quantity(self.state) - - def __iter__(self) -> Iterator[SimulationState]: - while True: - h_next_step = self.state.current_step_size - lin = self.linear_operator(self.state) - nonlin = self.nonlinear_operator(self.state) - self.record_tracked_values() - while True: - h = h_next_step - new_spec = rk4ip_step( - self.nonlinear_operator, - self.state, - self.state.spectrum, - h, - lin, - nonlin, - ) - new_state = self.state.replace(new_spec) - - new_qty = self.conserved_quantity(new_state) - self.current_error = np.abs(new_qty - self.last_qty) / self.last_qty - h_next_step, accepted = self.decide_step(h) - if accepted: - break - self.last_qty = new_qty - yield self.accept_step(new_state, h, h_next_step) - - def values(self) -> dict[str, float]: - return super().values() | dict(cons_qty=self.last_qty, relative_error=self.current_error) - - -class RK4IPSD: - """Runge-Kutta 4 in Interaction Picture with step doubling""" - - next_h_factor: float = 1.0 - current_error: float = 0.0 - - def __iter__(self) -> Iterator[SimulationState]: - while True: - h_next_step = self.state.current_step_size - lin = self.linear_operator(self.state) - nonlin = self.nonlinear_operator(self.state) - self.record_tracked_values() - while True: - h = h_next_step - new_fine_inter = self.take_step(h / 2, self.state.spectrum, lin, nonlin) - new_fine_inter_state = self.state.replace(new_fine_inter) - new_fine = self.take_step( - h / 2, - new_fine_inter, - self.linear_operator(new_fine_inter_state), - self.nonlinear_operator(new_fine_inter_state), - ) - new_coarse = self.take_step(h, self.state.spectrum, lin, nonlin) - self.current_error = compute_diff(new_coarse, new_fine) - h_next_step, accepted = self.decide_step(h) - if accepted: - break - - self.state.spectrum = new_fine - yield self.accept_step(self.state, h, h_next_step) - - def take_step( - self, h: float, spec: np.ndarray, lin: np.ndarray, nonlin: np.ndarray - ) -> np.ndarray: - return rk4ip_step(self.nonlinear_operator, self.state, spec, h, lin, nonlin) - - def values(self) -> dict[str, float]: - return super().values() | dict( - z=self.state.z, - local_error=self.current_error, - ) - - -class ERKIP43Stepper: - k5: np.ndarray - - fine: SimulationState - coarse: SimulationState - tmp: SimulationState - - def __init__(self, linear_operator: Operator, nonlinear_operator: Operator): - self.linear_operator = linear_operator - self.nonlinear_operator = nonlinear_operator - - def set_state(self, state: SimulationState): - self.k5 = self.nonlinear_operator(state) - self.fine = state.copy() - self.tmp = state.copy() - self.coarse = state.copy() - - def take_step( - self, state: SimulationState, step_size: float, new_step: bool - ) -> tuple[SimulationState, float]: - if not new_step: - self.k5 = self.nonlinear_operator(state) - lin = self.linear_operator(state) - - t = self.tmp - t.z = state.z - expD = np.exp(step_size * 0.5 * lin) - A_I = expD * state.spectrum - k1 = expD * self.k5 - - t.set_spectrum(A_I + 0.5 * step_size * k1) - t.z += step_size * 0.5 - k2 = self.nonlinear_operator(t) - - t.set_spectrum(A_I + 0.5 * step_size * k2) - k3 = self.nonlinear_operator(t) - - t.set_spectrum(expD * A_I + step_size * k3) - t.z += step_size * 0.5 - k4 = self.nonlinear_operator(t) - - r = expD * (A_I + step_size / 6 * (k1 + 2 * k2 + 2 * k3)) - - self.fine.set_spectrum(r + step_size / 6 * k4) - - self.k5 = self.nonlinear_operator(self.fine) - - self.coarse.set_spectrum(r + step_size / 30 * (2 * k4 + 3 * self.k5)) - - error = compute_diff(self.coarse.spectrum, self.fine.spectrum) - return self.fine, error - - -class ERKIP54Stepper: - """ - Reference - --------- - [1] BALAC, Stéphane. High order embedded Runge-Kutta scheme for adaptive step-size control in - the Interaction Picture method. Journal of the Korean Society for Industrial and Applied - Mathematics, 2013, vol. 17, no 4, p. 238-266. - """ - - k7: np.ndarray - fine: SimulationState - coarse: SimulationState - tmp: SimulationState - - def __init__( - self, - linear_operator: Operator, - nonlinear_operator: Operator, - error: Callable[[np.ndarray, np.ndarray], float] = None, - ): - self.error = error or press_error(1e-6, 1e-6) - self.linear_operator = linear_operator - self.nonlinear_operator = nonlinear_operator - - def set_state(self, state: SimulationState): - self.k7 = self.nonlinear_operator(state) - self.fine = state.copy() - self.tmp = state.copy() - self.coarse = state.copy() - - def take_step( - self, state: SimulationState, step_size: float, new_step: bool - ) -> tuple[SimulationState, float]: - if not new_step: - self.k7 = self.nonlinear_operator(state) - lin = self.linear_operator(state) - - t = self.tmp - expD2p = np.exp(step_size * 0.5 * lin) - expD4p = np.exp(step_size * 0.25 * lin) - expD4m = np.exp(-step_size * 0.25 * lin) - - A_I = expD2p * state.spectrum - k1 = expD2p * self.k7 - - t.set_spectrum(A_I + 0.5 * step_size * k1) - t.z += step_size * 0.5 - k2 = self.nonlinear_operator(t) - - t.set_spectrum(expD4m * (A_I + 0.0625 * step_size * (3 * k1 + k2))) - t.z -= step_size * 0.25 - k3 = expD4p * self.nonlinear_operator(t) - - t.set_spectrum(A_I + 0.25 * step_size * (-k1 - k2 + 4 * k3)) - t.z += step_size * 0.25 - k4 = self.nonlinear_operator(t) - - t.set_spectrum(expD4p * (A_I + 0.1875 * step_size * (k1 + 3 * k4))) - t.z += step_size * 0.25 - k5 = expD4m * self.nonlinear_operator(t) - - t.set_spectrum( - expD2p * (A_I + 1 / 7 * step_size * (-2 * k1 + k2 + 12 * k3 - 12 * k4 + 8 * k5)) - ) - t.z += step_size * 0.25 - k6 = self.nonlinear_operator(t) - - self.fine.set_spectrum( - expD2p * (A_I + step_size / 90 * (7 * k1 + 32 * k3 + 12 * k4 + 32 * k5)) - + step_size / 90 * k6 - ) - - self.k7 = self.nonlinear_operator(self.fine) - - self.coarse.set_spectrum( - expD2p * (A_I + step_size / 42 * (3 * k1 + 16 * k3 + 4 * k4 + 16 * k5)) - + step_size / 14 * self.k7 - ) - - error = compute_diff(self.coarse.spectrum, self.fine.spectrum) - return self.fine, error - - -def press_error(atol: float, rtol: float): - def compute(coarse, fine): - scale = atol + np.maximum(np.abs(coarse), np.abs(fine)) * rtol - return np.sqrt(np.mean(math.abs2((coarse - fine) / scale))) - - return compute - - -def press_judge(error: float, step_size: float) -> tuple[float, bool]: - return 0.99 * step_size * error**-5, error < 1 - - -def rk4ip_step( - nonlinear_operator: Operator, - init_state: SimulationState, - spectrum: np.ndarray, - h: float, - init_linear: np.ndarray, - init_nonlinear: np.ndarray, -) -> np.ndarray: - """Take a normal RK4IP step - - Parameters - ---------- - nonlinear_operator : NonLinearOperator - non linear operator - init_state : CurrentState - state at the start of the step - h : float - step size - spectrum : np.ndarray - spectrum to propagate - init_linear : np.ndarray - linear operator already applied on the initial state - init_nonlinear : np.ndarray - nonlinear operator already applied on the initial state - - Returns - ------- - np.ndarray - resutling spectrum - """ - expD = np.exp(h * 0.5 * init_linear) - A_I = expD * spectrum - - k1 = h * expD * init_nonlinear - k2 = h * nonlinear_operator(init_state.replace(A_I + k1 * 0.5)) - k3 = h * nonlinear_operator(init_state.replace(A_I + k2 * 0.5)) - k4 = h * nonlinear_operator(init_state.replace(expD * (A_I + k3))) - - return expD * (A_I + k1 / 6 + k2 / 3 + k3 / 3) + k4 / 6 + def stat(self, stat_name: str) -> np.ndarray: + return np.array(self.stats[stat_name]) @numba.jit(nopython=True) @@ -419,64 +38,169 @@ def compute_diff(coarse_spec: np.ndarray, fine_spec: np.ndarray) -> float: return np.sqrt(diff2.sum() / (fine_spec.real**2 + fine_spec.imag**2).sum()) -class ConstantEuler: +def weaknorm(fine: np.ndarray, coarse: np.ndarray, rtol: float, atol: float) -> float: + alpha = max(max(np.sqrt(abs2(fine).sum()), np.sqrt(abs2(coarse).sum())), atol) + return 1 / (alpha * rtol) * np.sqrt(abs2(coarse - fine).sum()) + + +def norm_hairer(fine: np.ndarray, coarse: np.ndarray, rtol: float, atol: float) -> float: + alpha = np.maximum(np.abs(fine), np.abs(coarse)) + return np.sqrt(abs2((fine - coarse) / (atol + rtol * alpha)).mean()) + + +def pi_step_factor(error: float, last_error: float, order: int, eps: float = 0.8): """ - Euler method with constant step size. This is for testing purposes, please do not use this - method to carry out actual simulations + computes the next step factor based on the current and last error. + + Parameters + ---------- + error : float + error on which to base the new step size + last_error : float + error of the last accepted step size + order : int + order of the integration method + eps : arbitrary factor + + Returns + ------- + float + factor such that `h_new = factor * h_old`. The factor is smoothly limited to avoid + increasing or decreasing the step size too fast. + + Reference + --------- + [1] SÖDERLIND, Gustaf et WANG, Lina. Adaptive time-stepping and computational stability. + Journal of Computational and Applied Mathematics, 2006, vol. 185, no 2, p. 225-243. """ + b1 = 3 / 5 / order + b2 = -1 / 5 / order + last_error = last_error or error + fac = (eps / error) ** b1 * (eps / last_error) ** b2 + return 1 + np.arctan(fac - 1) - def __init__(self, rhs: Callable[[SimulationState], np.ndarray]): - self.rhs = rhs - def set_state(self, state: SimulationState): - self.state = state.copy() +def solve43( + spec: np.ndarray, + linear: SpecOperator, + nonlinear: SpecOperator, + z_max: float, + atol: float, + rtol: float, + safety: float, + h_const: float | None = None, +) -> Iterator[tuple[np.ndarray, dict[str, Any]]]: + """ + Solve the GNLSE using an embedded Runge-Kutta of order 4(3) in the interaction picture. - def take_step( - self, state: SimulationState, step_size: float, new_step: bool - ) -> tuple[SimulationState, float]: - self.state.spectrum = state.spectrum * (1 + step_size * self.rhs(self.state)) - return self.state, 0.0 + Parameters + ---------- + spec : np.ndarray + initial spectrum + linear : Operator + linear operator + nonlinear : Operator + nonlinear operator + z_max : float + stop propagation when z >= z_max (the last step is not guaranteed to be exactly on z_max) + atol : float + absolute tolerance + rtol : float + relative tolerance + safety : float + safety factor when computing new step size + h_const : float | None, optional + constant step size to use, by default None (automatic step size based on atol and rtol) + + Yields + ------ + np.ndarray + last computed spectrum + dict[str, Any] + stats about the last step, including `z` + """ + if h_const is not None: + h = h_const + const_step_size = True + else: + h = 0.000664237859 # from Luna + const_step_size = False + k5 = nonlinear(spec, 0) + z = 0 + stats = {"z": z} + + yield spec, stats + + step_ind = 1 + msg = TimedMessage(2) + running = True + last_error = 0 + while running: + expD = np.exp(h * 0.5 * linear(z)) + + A_I = expD * spec + k1 = expD * k5 + k2 = nonlinear(A_I + 0.5 * h * k1, z + 0.5 * h) + k3 = nonlinear(A_I + 0.5 * h * k2, z + 0.5 * h) + k4 = nonlinear(expD * (A_I + h * k3), z + h) + + r = expD * (A_I + h / 6 * (k1 + 2 * k2 + 2 * k3)) + + fine = r + h / 6 * k4 + + new_k5 = nonlinear(fine, z + h) + + coarse = r + h / 30 * (2 * k4 + 3 * new_k5) + + error = weaknorm(fine, coarse, rtol, atol) + if 0 < error <= 1: + next_h_factor = safety * pi_step_factor(error, last_error, 4, 0.8) + else: + next_h_factor = max(0.1, safety * error ** (-0.25)) + + if const_step_size or error <= 1: + k5 = new_k5 + spec = fine + z += h + stats["z"] = z + step_ind += 1 + last_error = error + yield fine, stats + if z > z_max: + return + if const_step_size: + continue + else: + print(f"{z = :.3f} rejected step {step_ind} with {h = :.2g}, {error = :.2g}") + + h = h * next_h_factor + if msg.ready(): + print(f"step {step_ind}, {z = :.3f}, {error = :g}, {h = :.3g}") def integrate( - stepper: Stepper, - initial_state: SimulationState, - step_judge: StepJudge = no_judge, - min_step_size: float = 1e-6, - max_step_size: float = float("inf"), -) -> Iterator[SimulationState]: - state = initial_state.copy() - state.stats |= dict(rejected_steps=0, z=state.z) - yield state.copy() - new_step = True - num_rejected = 0 - z = 0 - step_size = state.current_step_size - stepper.set_state(state) + initial_spectrum: np.ndarray, + length: float, + linear: SpecOperator, + nonlinear: SpecOperator, + atol: float = 1e-6, + rtol: float = 1e-6, + safety: float = 0.9, +) -> SimulationResult: + spec0 = initial_spectrum.copy() + all_spectra = [] + stats = defaultdict(list) + msg = TimedMessage(2) with warnings.catch_warnings(): warnings.filterwarnings("error") - while True: - new_state, error = stepper.take_step(state, step_size, new_step) - new_h, step_is_valid = step_judge(error, step_size) - if step_is_valid: - z += step_size - new_state.z = z - new_state.stats |= dict(rejected_steps=num_rejected, z=z) + for i, (spec, new_stat) in enumerate( + solve43(spec0, linear, nonlinear, length, atol, rtol, safety) + ): + print(new_stat) + if msg.ready(): + print(f"step {i}, z = {new_stat['z']*100:.2f}cm") + all_spectra.append(spec) + for k, v in new_stat.items(): + stats[k].append(v) - num_rejected = 0 - - yield new_state.copy() - - state = new_state - state.clear() - new_step = True - else: - if num_rejected > 1 and step_size == min_step_size: - raise RuntimeError("Solution got rejected even with smallest allowed step size") - print(f"rejected with h = {step_size:g}") - num_rejected += 1 - new_step = False - - step_size = min(max_step_size, max(min_step_size, new_h)) - state.current_step_size = step_size - state.z = z + return SimulationResult(all_spectra, stats) diff --git a/src/scgenerator/utils.py b/src/scgenerator/utils.py index f1627a1..7ef3477 100644 --- a/src/scgenerator/utils.py +++ b/src/scgenerator/utils.py @@ -5,6 +5,7 @@ scgenerator module but some function may be used in any python program """ from __future__ import annotations +import datetime import inspect import itertools import os @@ -20,7 +21,6 @@ import numpy as np import pkg_resources as pkg import tomli import tomli_w - from scgenerator.const import PARAM_FN, PARAM_SEPARATOR, ROOT_PARAMETERS, SPEC1_FN, Z_FN from scgenerator.errors import DuplicateParameterError from scgenerator.logger import get_logger @@ -33,6 +33,19 @@ class DebugDict(dict): return super().__setitem__(k, v) +class TimedMessage: + def __init__(self, interval: float = 10.0): + self.interval = datetime.timedelta(seconds=interval) + self.next = datetime.datetime.now() + + def ready(self) -> bool: + t = datetime.datetime.now() + if self.next <= t: + self.next = t + self.interval + return True + return False + + class Paths: _data_files = [ "materials.toml", @@ -148,7 +161,7 @@ def conform_toml_path(path: os.PathLike) -> Path: def open_single_config(path: os.PathLike) -> dict[str, Any]: d = _open_config(path) - f = d.pop("Fiber")[0] + f = d.pop("Fiber", [{}])[0] return d | f diff --git a/tests/Optica_PM2000D/Optica_PM2000D.toml b/tests/Optica_PM2000D/Optica_PM2000D.toml new file mode 100755 index 0000000..2df3172 --- /dev/null +++ b/tests/Optica_PM2000D/Optica_PM2000D.toml @@ -0,0 +1,18 @@ +name = "PM2000D" +mean_power = 0.23 +field_file = "Pos30000New.npz" +repetition_rate = 40e-6 +wavelength = 1546e-9 +dt = 1e-15 +t_num = 8192 +tolerated_error = 1e-6 +quantum_noise = true +# raman_type = "agrawal" +z_num = 128 +length = 0.3 +dispersion_file = "PM2000D_2 extrapolated 4 0.npz" +interpolation_degree = 12 +A_eff_file = "PM2000D_A_eff_marcuse.npz" +n2 = 4.5e-20 + + diff --git a/tests/Optica_PM2000D/PM2000D_2 extrapolated 4 0.npz b/tests/Optica_PM2000D/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/Optica_PM2000D/PM2000D_A_eff_marcuse.npz b/tests/Optica_PM2000D/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^r74UqLloPB-HM2S zDBYo#-~X(+Gk4~Gx*skxjw8a$d-mSziM^hM<*I(NBc!;GNm6oSx8F@PmXiAKkG#|f zsZcLJx0(KdQPNVw{^xW5ecgZm_uq$TjqPg7_0m$|QZbsYZXqtgn)4@WE_652)S0O1 z9uOQF?Bu^CAlTLIzt5XG`G&aleLlp~DbTI&-&%UQ21Yv5bSB14{D1vWDAd`LR+L(R z=Bd}Wbj25-^rz28>!S+rynfH{OCbg5ojs;@q<;ZUs4%)6=2?KR=Gzy)axTCOVOY5S zrUHatA`78?0XCVO(Q33Rz=PYKNn4s1phIKbhZCj+IBbRDIYpBKyt7+7*LOhy#&oQ3 zFEcDaT}#hJ({&55@?e^3=JW#G_HTb<)#L(PV_*8vbV334J6XQ(z~}XC zZXR#vGNb?}U068lmV5!eQQRm$sb2wRyG4~J^h(ee3$J`>mtf0?3nFNe;9=iazB_(P zaLetcW!*m|I6`gsgbn2qydN^FxUBD@#2t! z!5s-&PQPa;x-LPv1UsrO|IHLi@YtevZ>w`8=v3Ht$}C-iBkM?nx zJFOPFhD*?(EkLC-P=Yn9^Xew~NO0$YWA_r>B)GQs;qI?nBsjpfCW_f0!E46`tSz>Y zV36!oot`xk)O1)Xw{e98Kb>{{a(Rga`Gb-!sxlH>=CX?)VI;w>E9mt`Pl6|vJ}ow! zAwd`KuCZcG398;wZEcz)!N;TLKCx1hU~-V{!BZn8IQKzF;?Us|teq%a?>|U_yCYBO zJ(iQ8?XxF?bfhFWNbAp!tS&L$NEpP{w1_eE?X<$>^8&8X|zFBhn%26@yTXn8@bcq<(?Rya(xm%2i-3=CRw~O(X_0VZ^ zBw~y>GDAj^FGg({^N-D$V*Iw=`=+*86rlX3b{>D@?xC0VDg-IGGcr!Ts(Ac zuL!e#IaW*BMYw2j!p)XO5jJlt-oEy?2oE)04LS2ugqv1;FdtqnLZv+&QvyGV@NUP* ze@|YEFvimG-Hc}<)IGHJbk0K&R{jghs=Y13B745e%IhMuIC)~>;mabF9q?F1{+tLe zZ?3I#KOsW@3j^=pIV3{$AzFJT?-Su$H`C~pViD$`>)P*yA~aD+(J>KS)LD?AZvEd?2(3rO*I#7iAG4^Wjy+v64BCK_; zn+S`ii=Py15utVRnS-tCMcDttvqT#k5ni2Lzy9nR5e8-|>MN}f;goXi!9hzz_`zsp z!09^aFlX(YlW)oF#fdLrytbRh7zjtGwzmabT(DMF|Ea`hvVL^#S!sz+W;gpYQP zf9)|+gh_1+Puv|Q!Z~X;rcW6p!e0lX98={)xU08tUd2Bl+Sp!DTHGbXfyZCf?rsqx zE8BegU!4$x9ft4RTq8uybF(6_LWrLRt+7=5B1D0U@3feALR@w=OXkB1A$BVr{%G`A zh^M^ooh!U2M7LX&xotOvIA)BjhwU{XJ_%AcK6gQgDG!WCE1wpk(WEU+K}Uu7Co-<= z*#RN$c~-Q4)@~u%X(i0iwt`onx7Dz5!>bDp6P zZyo8{zIB!mBV|U1woen{j16;F+f5c?`Dxp!=f?}NU|{G!<NS$N%+3#``-lU@Z{!N&5r^cHualM+G_#ciTf*8`BZ?>uhqUTc__eHGZtLhds~1N zsp|#O*9ExsQ-rtUWdW`+C|Gj+tN>+l&Z&(*A;60jFIr;{3D9>zT0T!K%d3|iX04)b>KQYKgfO1z=W8|=RYY3(BS3byajRs zteNiQvHc$(i<1+KJ3IMk{h@gDx)weT(7V=jp^lH&GC!25)bKH=ykmc81s|u*8x{ZZ zGaoj4{UhE$Ie3j_cizUc%uHq)#Eq$=wkL*Jop+PRd?0; z`CQ=Rg#e7td zyRoTXJ|C|;e3|2%$;Xg$JwtD%@KI~<*qRCPeEjS(?^b*yAB9)#i#`YQ(M&l!bfF&~ zd%VS~i#+&v>eiX5T~2&-8}t0%ItMKhFBIQU)ZeT!LqygBRpo)gpf7@jT_qd1w5+F!@p z_>Sk}w>dDQOqGw~yp8>*EAw&1 %ihTSFTQ5}&;Nw}z1%asyAD?|V{CZ_)KF-Vf znJ!@<0lkKfGRTmHL|k1N|%`r-L}JXd0UKKe*L zj2i-ow~6wx@tpjlvl;og(REpTP+UIVR!-KP6`YSVOWE(g-Scs4aJTr# z#(b1kH!}0I%txQ+p`BBf<>Twb(?=`j<>MmVx=nj#=HsES(}!)*$Va7ok6UW0`543O z5B)GWA1nWSo4QpxA1#W?k^fAEy%;bh}fqrFAqOVy_^&q zorg=_wCc|f$iw42y`Swac{s{2VC(sHd6@L$=*p1Qc~~o|?$KSGhc-*5oUEIZhisFR z(=qKl)ZDvM*?VFh3amcf(NxL9E-BT>%7J<4cFJb@zQ4Knbc^7VV{hcp?@y^ z2);k#M0XCZQJ4GX(~yG~pPLA$Rp#KNr1+&(?{hFq_iEezCpp;qwf&HDX%24AH`qAw zQV!l{f)u|T%R#-mIViVUbMmO{9Q5yTdGjhE2j3svogoa(!NnUt z&$sl-!D9+z{wQwFL6s|ZyUMI{Fu_Y?m1&lPHLB10nJvsg>-&Q*$j;6|j9B4$Yf27I z)l3_cI3@@AZ*Dz;p*h&e`xoEcF9%%=7tX%ik&TakL{^0V$;LV2^Adya*|>X2lUdW7 zY*c8T*>$ok8$sUt_azo&SiO;-GB2InIjN=$?(*k8>MTH)f+a_Tk<)7TNe$dw{KYX*PO( zTr$9VUN*kSO1wODW;QOIk3Nspv++Pp=eX=q**J9DoM+1iXJh2@kmUZ-+4#Nvq`~dh zEHpp(`)AUxEWBVhZQG)+Svax3`>NiTS(tf#Z|{}+S=j9Q^<*Swp~DC@=Q(Gx@NVh4 z5iKQIs2d`=d%7qKx2wO14B%&>?294VGk96(m$d5huc$11tDBL3$Uh5>zuh)*cgeyd z`O+;K>#|UJ5mQobo`vys(d%}ZWMOsjMTJdsve0UE%e67uS$MTaH{jibES!AIZ&JaC zEX><@n#T))N6P>SU7(H#u#7AC^zjCTF(QtI}j+LJ>ao7F#YviA2qI~2? z>AQC_F<5i$*_11p`0-7?o5{&cG);Rn>fgRhJZU)a;kB)qIJ$CKY;m}%fx;D<&*ae&%`0)$G_d&KNBM!Z^%^Z&cO26!UZ4dGjNs8 zt3QQ5GVt8TVY_YLXW;nk)z->SGBAC9mhAIe8Q4^F_d?#q4BWV_pO^Wu47{`4c&x&n z44l~!`}m$D1GgT$oWRST;((Nh7+f6etHpF?$`2`vHQ+?O)5Zw&4c=54RcS;6cNm3238=Zj~vu(AG4avaV zZ-NhA{W7pE|4EK!M>;w#8nmeLcRD_-Gi%&io{suEQ};Q(NynY5Z`zJ8OUD7dy#qd7 zPsgBR^RH|>myRDcM)=wvPRFH#&W%^unU2S=G(3ADOvh2)GgA2J>6kp)+i*>6I@aDV zsTvfNj<%5%+a9>3Bda;lJbgntYQ1&%yVN2b1!>z($t+FBZo@auH;mHJy;4afenvVz z6Z-pH1h3f4r>BuW=aDbxkss%g z59g8p=8^B_k>BQ#&*qW8=8>=Fk)P&~kLHno=8BN!&rKzNn@YYmmHccf`PfwQuc_o)Q^~KUl21(~f0|0ZG?n~lD*4b<@}H^X zJ5$MTrjpN0C4ZSpzA}~kWGeZ{RPv9he*I$?qhS z&q*eKlT5xQnfy#L`Isc~FG=KElE|+lkxxk?f09JLB#Hb;68VrM@*hd$JCew6B$3ZZ zB7c!Yz9NbIL=yRkB=QeQvqNFslbM7|)2{6Hf4fJD;&iKP1zN$)3;&QBzL zpGdksk@S2b>G(v_?}?<_6G^Wpl1@(~eV$0VJdyNxBI)o%(%*@syAw%oCz8%iBz@iY zegC7Y6G%@dkd96u{hUC$If3+Y0_o%g(#Hv;ixWr>Cy)+KApILpx;LKmZanGSc+$7= zq-*0z&&HFEjVJvYPr5ap^lCim)Hu?oaimM*NRP&m4vi!I8ArM^j`U_6>C8CNmvN*k z<48}&k&cWd{TNHSF_!dVEa}8p(uc963u8$S#*z+4K7F{I04NRLI64vQxJ6-~M;n)FsQ>8xndSJ9-aqDfCh zla7ie{S-~QDT?$`6zQZW(nnFGi=s#mMUf7QBK;FZx+jYCP88{!DAG5Pq-!Ec&qR`r zi6s3JNxCJH^hzY@lt|Jik)%r^NsmO54v8fF5ka~mg7iiN>5K@{7ZIc@B1lg}kdBBT z{SZOAA%gTm1nGou(g)$B3&Kecgp&>kC;ks7-VZ0f4=0`vCw>nnUJoZe4-UuST2qK;cB7O)W zUI-#S2qGQ`B>WE~+z%wY4WB}+zuqX4kVloAbbuWTn->S z4j>#3ApG?w-1R5C^(UP5Cw%oMT=geB^&=ehBmDFu-1H;7^dp?~BYgBFT=XS8^d%hh zCH(Ux-18y4^C6t`A$;>8T=OA3^C2AbCj9ay-0~*8@+O?}CVcWHT=F73@**7aBK+|p z-0>p3@g$t_Bz*BCT=674@gyAaB>eCo-0&d0@F1Msar}}rNx_75~ zcc(gcqxyEEx^|;_cB49WqxyBDx^<;`b)`CWrTTQGx^$&_bfr3Uq55;7x^tm=bD=tO zq55*Bx^kv^a;7?RruuQFx^bp@aiThLqWW;6x^SX;aH2YJqVw-a=iZUdyCa=*M>^k* zbgmugJa3_MyoJv1W;(Z<>AY^HbGnJn=O#Lro9H|`&^dIV^XEY4ZX=zyjdad7()rpz z=V}9;rww$D*3SGotJfVPS(-+u%~ljPv^m&&Ve26za8zn9qqR*?XxZIuPyDX z4eh55?V}Ct-&)$YwX|Q>v`^NwKi0G_RrA8X zpGx1SL0>vVv-SDr~3`9>$E`hZ%79 zQp@T6tN`wwZY901mB8IwuA=v64%|JiIla%-z};)Fq4&E6xO;vJdf%47^*t=<{aXUp z2eP8)umY|xW<}3q4P2klnx4xVxW4IHdcL*5^^tApIc)Y^L?y3|t>%3$61O;QBI- zwBC-u^@$v5-5r7JTRGDDI|A27bE1840zljK{<#3x$9JK9bOkO~;7a@H3S7>?mG;#YxZH&+?XMef zIS@D6XE)$-DQ>jiZouVa+-Tq3fy?c<)Bd{ycTea}=fEAfdq)pC4<5kXV|viJ@Br>! z)Pv542XObSo^(z;fxEZ$r1Rnl+&!=-of}W!?v=gh{CEL(PwhqL$P2i8Z!bDe-oV|X zd(*k{2JT+oo6eUvaQFP)bk2N$>l^scdGi6T58*@S&Ih=@h7X-TU*P&AzH|4gju?7(nMb5V*c%Af4|(;QE|_bj|~T>zfABc@G4x4;x75J`lLRZXlii zK;Zhsfm8=U!1bMjs2+lV>thE|T?7Hw7Z0NP2m-Fp9z=B#1YF-fi0UN>xEw$*)lD#P zxq@J-pJ3o}3c*xI!NBDnf~lT@fy+?@Q(XlEm&*vI`U(av=Mh4676M#uB!ucM1h^bZ z2-RH(aJiNcs=pB6axx)Qhate_ZbGOYLxIcjgi>9G0+$O4rTPp7E@u=ZYh-N zH59lUR4CPLC~&!|Fsk1$;Bs1FRL5b!<-WqGp2L94k%duRhXI#M3#0lD11{$lMs*$r zTy8F$>OCB|99}rpeK>HrzHqAlaNzm^;e-R>!1WQr2@k@7>pO%KE<^y=r-&eYhybpy z5kWW+0bCy>g76{&xV}jQ;YI{-eU=Eqj|kxUGLeKMk-+tFA_-3-f$RH360Sr7*C&c3 ze2E0EuM|l*6A4@&Dw6Of61cuq6yZ)3aDA>Q!k;MM`eIRpLs7u>(V_^CqJZnWMG-DV z0oSLCB7BMluCEtOI28?CA26EmDjK-HVKm`ZG;n>!Xu_{(;QEr$gk#ab^)X`z&tib< zd&Us1#Q@hQjUjxC0j{qaLpT=$Tpu=u@Gb_pzHJQQUJP)3-dMuFSm64?v4n%M!1a-1 z2@hj|>pRC1F2(}ar;a6jj0LW*9ZNVF2V5UKj_@)LxW0KD;bt6gefBuQ&p6=v@^OTt zalrNQ;|NdVfb09m5w6AqmlKF5e2oV#R}fD)8xLF#A)fFy9=O~>JmGFUa5;x~!ryq{ zauErH!wJCUC=v*d6M)NIBoHnq0GHE9Abd^$F4vJjIGq4=Url(O2wZL?k#IW^xSUBM z;ddf%xs*i0@kHQqEQy5YiNNJv5((E6fy>Dx622z_m#ax6oKFNUhm%Nnp9oxTCy{VJ z5xAUBBH@1`aJisF;(iRTA+>5^%Y!B;t`I;Bs6^#3xC><-U@LSCW9si6s%gBmtKzOCp|01}=w| zOnj3JTy8CycqbXSoLe&SPcm@1xMbp?WZ-gi$;3y=z~%0eiIq{n{ zN(L?mm`r??3|ww7nRqK1xSU}!@mDf%xx^IWu@vBPj48xtDZu3(Q;646fXhjy5Wl4W zm#a)6o=X8PhnYfrmjYaFGlh6B1-P7N3h`eGaJkSF;=vT)a-=E5hbh42PE&{%Q-I5< zrV>A<0+(w|C7w(LE(e=Re3=SdZZ?&8GZna;Z7T6+DsZ{nRN~Q8;BvgF#HXpi<$hC% zS5tw@38xajrUI8MP9>gA1ulo2N_?9NTy8m)c$WuU&Y4I2%L6VK%_AP>0hgoZ5g+q_ z%U$z`mwCYDw0XqOJm7NOJmP5{a5-=u@ih;)+&GVTn+IIZoJai411^`&BOd1gmt*G< zpYwpry{8ecrvaCfPa}R$11?veMm(PeTn;~t_&yD|+3}q# zya4HeG@v{I>4G$%yaDNhG@v{J>4Y?(yaMTkbf7!~>4tQmyaVZnz62zfhaer14lI|K zAU%-|l&2tFkq(r%AbpVzl*b^Qkq(sCAia?el;3zlt(0;m;;noB)yme zlxHN}m;;n|B>k8Jl!qi8nFExUBt4l6l&2(JnG2M+Bz>6+l*c5UnG2NHB)ypnl;kBSlm{gpnhTT{B|Vx8lqV%!nhTUSC4HI;lt(3(eRe4u#(@&y8* zc?0qX0-$*W@(BW0-$*k@+AVGc@y#{0-$*m z@+ktKc@^?20-$*o@+|_Oc^C380-$*q@-YISc^UFELZEpX@-;%Bc^mRKLZEpZ@;O4F zc^&dQLZEpb@;ySJc^~pWLZEpd@iLZEph@<~FVc_s2o zLZEpj@=ZdZc_;EuLZEpl@=-#dc`5Q!LZEpn@>N2hc`Nc)LZEpp@>xQlc`fo=LZEpr z@?Aopc`x!`LZEpt@?j#Nc`@>1BA|IP@?|2Rc{B27BA|IR@@XQVc{TEDBA|IT@@*oZ zc{lQJBA|IV@^K=dc{%cPBA|IX@^vDhc{}oVBA|IZ@_8blc|G!bBA|Ib@_izpc|Y=h zBA|Id@_{0tc|r1nBA|If@`WOxc|-DtBA|Ih@`)m#c}4PzBA|Ij@{MAkc}Mb(VxW0Q z@{wYoc}en;~5}^63RY^XlZ+3xMX?$+s5(&AXF-F94c{Cm&w`G%rtnz5r;Ro_u`)(7ZkQ`vRbO zeDe7PK=b(#(YzrxHik3YZ~l`C<4|Y4**%w6arDx6#U8l<}@vig7n@#Yt_?@$B3g zYd+2Yi{*!fs#7mjf&p(x*rIP8@UPfM>|oY_Dlh7$RxfD6{x2$>4HH{nj*c;0DX0Vo zhg(lZE^Eh_vbPaWb?;%#pS4Zzw|2m=9}riU_!39`(nf{A4pcDXb==mfhAoA5Y|*ne zFhBbvhQFx^-(7!lviN5+-d*6H^JhgncyW0e~kZ$&)SZ0;v$}ggEvFk-u=9G3KT5$6ed`a&TAgr0C+}X%_g8pNa&Ej;R|_r* zHY&3DSpjh>Ui^rlI&fNOzCFaF7IFWyqix5(;dmzvqtkO6z<(EGaAMG5NN_H4dNr*P ze?QdS6Zh*Q97?y7BqudMVWf#g>E7R1&fh)AFyuETUOrK1{-g!i6Yp&Hz4-x+2Gi1& zKXl@QIOTEAPrk-wholOd7xhA`>4NnZ8~5Ni%c6d~bANGbz}fYCK7541yL;E38vPfH zhHZ;H)cyvlw`KniRXk+#~i;{tZ0QWc$eDj}+spJA9IV;yawD@LfFYn-u%e zX|mmZkKgck+-OO(tQ1poBu4%Ip;jFIy65ZK*FE_DY>Mt;`5qXb9;f!pumeuS=j%Hr zNU{3gU1qHhYQpWiW+`T=$S`k5jXliYQw4bv+dn6d%do*K|8{0=dw~jdr&ddu$uJjB z{xYq*a2q{!Wt%-prC8TLi;hIgS3tsGWx25zdthC%MaKQudYs(w_~V7uZ8**N*>b;I zEzl{eqCCLz4;(OTYRcT)fky@R^RvU=wkq0e9k8wo7Zttx z+F6a{weRnMvih}0XR9^h3M13uqWP`(e4>}}3-xvg5Z$)Akk?Osym0Cezh2xF6+NY2$aUDl(RH^hr2PTf#t{qoa3t}N-J8^Y0lr+4s3fBcm-jupvZ4uptP1j-xv&|RuMmu|c;5z70^V*48(Rls13d4(uI@l>+Z9XY zbbcT&sG;=Ux=yHHc@Q5CeF;;Prmt$sYs0kK;q2N;Ct>B}1n-4rEwJ5WoYtIeZ}DbA z@vtw6e=zxoXX~7r3TXWD?1VyWIZRM1((7DZi|3ckvzqzwB2H0D8usGCA2_^UK0ABl zJJkK{H23fOI@D)A?3Mjp1NL&eM?8L73j<#GhaYih#-P*=gRk!^@j~sMiMO@e`zXJ4 zQaXHsb}RjTPV0N{^VP6Hd*b$@YUk}WA+vg+zl>Is+wPY*qE4m6`DPE+Ca#pZHlZ4Z zNR2($&#wz?8kTrFiW_m?;C>_T{cc5_=l|{|HMYYZcjM@H0S#c)zimX}xgMOcEAYD! z|2ul`*Q`BhF2(3c!Y(<8S9_x6rx=8`G!S_$Ej*yJroR8#SN} zTDHFKI5@8#>o>Rk^1-e;e33uv&eGws%-BH5;Hr$z@N3F!=j5mT*hTt9-#7Rkh7)r) z*ZSO-VMgBY|DJZa0xu7|;vkwU#U35q|LKB+2GH0vafZvYPN=>yubcO?9cQ(6g{j?d z!nj*_YfQ3p0{tX11$*D{T0Nt;4qY z9Bc2Z_ev*)Qkw=)wNpKGuCLw=H}6z=A@0PD5f{HqR_#IE4WD0(D(iyzr!2QkI^7OI z{*$k@%Sf?%M+Teqw{Af5+JQZ(knJvY-JHGSTvE?diZQ3;B z;H7_HRhg?dYr#`2kG5c@tm?*$9Zoe%m41Te%p0%$UbaD7!up4BsR7G9J8Wg@8nE9w zLz`f<;x-Knh$1UHJvoiOV;sJT@#5e7o&^oQF!C&+kYd5)Eac%8H zbCV;oE1&!XJEwi)uC%v9|8vhPR@Kzwd&|LZt5RAqvmLfQbZLQy-FFw5jBNnkt5b%a zD((2;PN_6#{=^?Arp~dfZ--mqjUir&FX6|MQWo#2S=Nm6 z3ckI(?DZ6`6-+H2_jcmr7*Seqsgt_Z+MpS6ZV`?1s3 z{ZJaW2>ydxhB*h#l_^e;n_jAX<6I1oR^SJJe^s%r#jnJ5zoqjC59d<2Wp|q-B z7fx-PoV)8tD;lUx7$ZHm7Z{^|Z}$KA1}~n6-M+s}ihVw!!X=-1foPIYr~X!&$#0vL zp}guje!P12?Ua*JY-oDi>o0DVaNyvtyd~?U7?~4~*7a65qow85eX@I+(eC-Q=rJKJ zp#JBwg^7a{^X_%GoBXXdWUKhShU*&9=D~dHDEnSGW-upgnqfV(m!GjbaJm%@FZrCC zp8pWX`(J*l1zj-5w`}ik%Tl!NJUcOETQl}IZ;??PUJm+)<&-mDHbQD>?~e?lKj;?G zw)geIpE$Q{%8*}Mn?W=1t z_|`WoT}{3NZR<9i_*w7?To1*_==(=FJeK#UP30MW>Q~meLA?&1 z=G_ihz*fV>(Be~#KVP7SmUL{!zLRKtYV-Mo;gyj8#_3Uc&`lUAiiv2s-hk>4ABw_? zzCfb-5_~NG3m(*r)X*xc?b8V_^3Sz2U_gqEk=gJ@NHWdqJs(mF@v7#YgGBA9GGdwD zuX7c+_}cgU^9LFMPF?Cq9C;I$n{OGbx<-~Q7tKEVBIXHB_`31m*4KX_alqOltByC= z*774}s*en7XYN1abZY}Bo6TOWcp1(5*g-Xpj*O)KAe!1dNSMNLk&LJKVGM)@9_igJ^EoK)ud>DL~+?K?j9d{fa{RE_Ij$u5Ti8FqU`*`nDtQtSt*oi|Ko+{YC!6y*MV z6~L0Fhl`){y3xkkV%o%h-|&Y0*L_pf|6tpg7QH)ht#EL0tEFjN2dZm^^Q6Chgp^5e zdaU|8xI3|D-Rx%#pzAPTV@bbebRAc>d5%Rbl#R(X-EG(l?(FTH?&2@dyZHUbn5W9j z?6VyfbAqIq;_3G@f)6P&h69v#xNMhZziqfs|Fh^1T-rRSXQ{$x7<+QyF9Xk3Q2OEi z$9B~ywzp*Hj|J)%u

1@uJ}`QIc|f!G=s(Hty~%n^i{@+2UbaR+rv<4qJ>e_CDC( z1+f!shoI80}k= z%68rD&sc0bo1XBj8dP~xM{cvMg;jT>qpiC7@VhEmNB`U)=E`aHVZk@bG3%p?R=e+S z6c|h}I$!V(hiQG7646it#Sd>4c}8~QtT!GWJHGX2r*Dc~uu-QQ25l;i-nO*`GahJV zp0lXO{@JeEU#AUVuZ`L>SU6@FleFs0G8ykq(7!d#y=#aU%v)Gz^Ezq(Gf*>l^5=yL zj9%x2U=8zuj8*ysJr$uG8z^)5`1*!{%(Adw!~aJ8K^S)Ju-@2j2yQ(pr}IXzdMZ=i zJ)HOkizjSdmbc_DZg}tT>Y}9*TV)aPDr-nHnpcH*_n#@to*5dmB+OKU@qMpo`N&#r&bsZb)4D(jp~afjM!%M6dAl z7nGj4{awTKVN9-Cx0YwhRQ5oLg`D^1cF5{=keTyAi=EaGdbj=21L!%jN8NqC0vq~s z>)E?;gIKMq(J}*zC$f1T5|7XMI*3hvEAu+LVi@~H9UFCTj$n&)W4gop>NsPVabRzh zJUj0FL+usKuOPSUt#;MIf$U**6;Dz8$%(#+I`qe?-SmDtyzV_jFjlx7TSf?rhj_G6D4?{Y|g zpu{+Tefr?j-WFW??CPJAF~gZrDg$8t3t8r)&bBGTH5HgoDkYw33#1t5=*tCTbn4;s zw(OpNl0i)OTlwNr?_b#Ks^^p+r@~lQs3sKOl4dT2{2lN0W;k>Gee0qJ%kPGEI+^Z)p$0lDBw()Su^?@ zEn+K2DY82k?U?m>fgF3gv3pacrxJT>O2?A7-v=@Jje^)i*{yi|aK_5->l(07cz1I8 zr`Mo-Zt`sJ8b#*tjl=(@Umwg|*!8&Q&eCzrc$e`JVPj@7BUStgAKA`fMg}i-Q$G0t z<~JR@nb0Z6^0#d`IDf!M=E!^1ofje%n2Qc`pS}<3#Ebh~#rkvq;Mud8>BsF`A!qoq zrQMnGOx-w>FBc+}Si8Bh%cFJVnH~e%6}n6N)_?4y^1}YZ8Nt%cqXH!|Y^$a7->WNR znDNfXcBvc~#on72@_fKvS@wC}G>sX&(d1DTaPn->GD1~YAD=8M%&%dnDJ z*VXG6DX}(F-OSdwPh#E&E%;sIBgdTSNL!+v@EU*2FAz@e*$-vSIgf&)$1|oehYeKk z4PdHw_%^FF^q{qN(XL0&rI^Z_mF;E+lo-d@&Z=3O+Dw2Idod(u5}Q15sEb#F7Snz| zIQR2Pbyn&8gJm;QF53dU)LFlCV^3~5qr&<(&YE&y z{RB4X{_io~X2V$@lU-xddOOf&k*E2wjY_PW+USIDu{C(0$;VZ2duMTBT3OsK(I<-OU+s`^pgJjs#=cbCa6qMMOTfQcYP@bLr z=EMTWFP)g%aK=YuaxJ4dsvGwaUnDVfOp(Dj^l z-|i=BjEXrNraYLD*N<5tS2~qVUA)oi<%Jj6c0}Q4>DUkO!{mToc+e=;fM-$ct)a>) zloWS1w<)s9tIw5cFPzMjynP>izgUfRGdkJY5%dgQtHbUFt(9lk;_^SwohLB&XNFEL zuxmxtZ}+Cwd#W(*pd+4qi=gdp3Z{k&b4*naDAB4i1_AM%<0h%+&QhtB+?? zW-eR%+Fh9)@we>Nrz@k`k+ECn{r;uOKFJCf7i!OA@9Dgom~W-b-0k{X+B96BEwDcB zntgQ~{CM z=op<|`&2ZH{dsJ)=wrr2CVcD{lhMWvcx225r~0QIxbMT&QwJ|7u{j3S?~44?n13b< z0+n;;GT#FhX8Q@JGSwR!Wr8e?nZ`ql6wdKyFm;Z9vScfinXh~EKeX)`%^WRL8pzKw zV*Gx_${7{>f`_ASR2s0l%!H~j%d|U2v1XQ;Z>9cfF_WBnXI;Oc!&>#sn$FuI&*U%M ze!16f7JJg6G;&1#NVb1jME=f0T5ROKf>+DK+Wa@*4d`!$%> zGYR_FMoBXjD_yUqrf4v$dB>z~mZ~$>iuW=T#tmfRcU^On+p5IKPJOvi<^4dW!04q$ z<+~}&-O%F!Q_VHl=WWgFbaQ7hcdP4a6{?oAX@^ap9WIn-{9eDC{tk`V82KLp{XU+A zMd$PSn-5^S8^f&&hiR}Qx4ETPsu;4}MdFCDz|7$p_vEus> zbLT0nzF~Mw^sXMzbXjq0`J-<9?U0U>UXEkCV$Ii`F`mjkH7}dt8#Rr6)Ksi9bfXq? zVfQ7E=d!cegQ&3MY>PgV@{1Wd?7k{{BV}TllJXSB7&c5S)vU%z&u`W^>dj#;G7DZT ztkP!gK6`WT=wn^RoCu#Yg4eHE^AOXPD<+bWaiSz1E1=D_29RzMfE#xXfo1Qtz$m_ z8^K;Y{h;D{lO9t&#^Pno{d#z@YrAL2;OWdev8`hB6jkPWY353gmf1{QgAul@?Z@y} z+WXE{p2(bR7*}}NWh&FPrQBw^!9?b0;oz*|S(<(Qiq+`~&AO~&=gQ)1GxeG84yz2R z1@qa5NB+%=EE&d}w8GU=(1;7OD-(|17|rai@AMsX-H>(8x;fFm`wv)3xAe~xjAiTF zyj7Ry$+IZ)eTiA2HtVot`TH3M4?;|)Vu*iX6O^o2{P%oc{cd*;?Z_R~j-3bn2EQMn z&d$32JMv59U-%;Mc(P%K8N2*?)rXo3lNg(l6a2xAebJCm*}|zIhhT(jw_#|MHYB0afFot+eK za_#8zNlbEXjl%?s`RrKt2-{P#X3VDpa|2bjD>G)M7xT?_7%)u>%0|c*O<^7-SuHAE zvY0t(GSZ;%&K$thNu8tS%ir)|F#5a5XNd!|bD-ME_z$a@1CfzQT8kW+ntM<7EI+o2d9-Jn`?`f1 zY+R$xSzBv0X40A2X1kuQW(+m&4@k9G%QUgqD*I^%F$ypC89qJ5V|dS^KkqdNW&S$b zxP5*T!3;G$anZ!vf=OFfcU;U{#k~K1P;q^0-{(VGycWiJGHd2eD;X2$#*B7ezGmKL zGv;8^Ou3z!KFrU>%in+Nse!+0lkC!aEtpkHui3^QH)W63YE;i2s>xWbJZs$XY8e~7 zCp^i?G@beSVd4t^o=Dcu?PkAGJ$cN*Yo9`7D$`iKyfGub2u&Ers%q;_Z3{LEyFRE_ z1~WaA^~#6U$Ff6!SI{?x%6R)?MwT0m@tpEWPxwDyoAF4w) zSuSBs``R^Y;B7qyUkoPCjDgYG3!ZM?5ng-jiL?u>}SoY^5`s2c2oU=G^h$;);}q(iMcLe z4oB`xE45E$j#$r)%=O#Bd{`0Z3l3WSO=R}WUQ%-XayaX-mHGFU^q;Q&2|OtHBKwp13wpe zoLjb>4GUA4dvaa|Tf3xFrsJp+yGPpPZH==(J9VMTc=zj@SgYtUmwpcoW)J>_4VR|d zu|eO|etSDlVtj6RV8GEitW0<1vnzUL%z%O$_rIuSu+i6!FST>q+J|E&nEM*V?Dc0U zc18O8*wS%2145l_*tScI^sp2YmhqMuyFpFJTD+Yt-8>|nT|9^1U=W|rhJ05r|7sV) z&c0A><##ZZy?-gGEo7)0+it$cXNyG?t83SLbf2i?3E*a z_5!PPcH8>02ED}zY_)*}YuB@f?OgbHW#Pu{teoI@vY~4(t9dXl_Nm<}ws5^_jBgT; zeSO)iFs^eGdug@k!MZpBo4^*nnzt&Lk!>@46IT$$o_jm_Z~Ip_=C;YsI>GZ0)~)(^ zjD1ETWBaPCJgLx`DGGcPW4morAC4@PZvSk_3?)2t};elpM66a ze_+bTNX#y*hO6BLvOvBP;swSCll*&NA=s?MX~td~c-aF(q%c2*)v8%8GojI&4HzG6VIIDh!D9`Rrs$?JAI%k%S&~COb0<3A>+bYKrnYzF_xrE5GOwomdRL+w$80%q zSoYz<66RIQy%j3Uw=$z69&Ziv+|D#^tFG_3!(-oR-g_$foy&yTe*GppJh#t}Hrnv! z@R`kF5byh z^^@rnTK6%p%#0&co{O2%>HEse?;T~Vo@=)LdYQ`9_3qd*W$Y0qecs~WJ*&gm<+|${ z9{ZnW)LhoZBqtqaXHLP$&8Lqr?pIIStz2+{mA|>wTYjpD`KI~w)2M}~Sa08Vv)k1=7e#5EWliqPEb(wV!3MPaikF{%n7un3eysIA%U)8w zv<@_OuqHCaCOzZM{x62kI~)imj^l_BN+M+MjIw3q&d3VcviBa@dt_I#x2(+Umn8cR znN5|2YN{e_+02e&uOuGx9Zdbx-B21;zKC zF$4$aqs#{up)7V!0Ps(JEq^xyX|_8&5M}KF1CbdOW38FMBt^9(ilrN9ew!N%cX|jI z{BF&%D)xh-g$@JO^~XRg=AN9-QY)~F9UduV!oc{?S)p#FW^l7%!J>@46Z913mR>Nm zMA?$+;jm}hqwOB+YcBs8_CmZE57mo@TMgjj92MbdCMQFRNM@~{L8|0(p#g9)* zQP=wX)kc{dV5FMzGp)ZC`7TXe8T^xtyfanQY9dO}ySadd5zQo2Bc?rh{!}fRnZMjT zI9rBT5AKq`e_MuHH|(lIjOvgS{paGv+oeakv_)&|OB2$yQTqJr(=F7`?^pFM<01OD zUAD8d9u4R|`%0%9JVIe2c1d<$ynyMBXu6S48=_OQbMuL`LjF_5DRswdk>%Io8sWo- zfLI}Na`1E&qPg9BO?K%KkU5iI%J?@0eY0b&T_NiL{jo!hEnT@s`Aph=5^)zeA@Gxf zouLUtExwU2Fzf(Do}v_cn%w{fR_K|1dj!J1eU<;Gjezs$RNf(BIWQnL7LxJp0B_VU zE*HsVBUA0l&ZcPuQt-&+G39E+E>f0AM1g^sq}=%9Mkm%gx!;d|e-5Z|ifu+! z@-GyIjeF4M+>qtdjUJR{7WYzcw+Z2_nVfKbjJ!gQxsszCMDLa6;iTSM8S4 z9Epyu{?As#%&kXctyTcs%qT0*do-Z24>zeExaR{B=30u6%q58FQHIN%*;3Fv(_c{J zl#FD{B*$+!WdW`LNy6R7@qqmPm4HL`Oz^;Iet4+R5UKZ)-0n%U0dt$4R1thP5mUW| zX^gBlD!rA}bZajU?9ieL?tlkqyQFaZAy*=Be8MKH^`s8L``a7|R5pP1>dPm+zK;;a z)Rk~dmyg=`h_Xj7bs%ZW8U=oZYIGabCkqU9p@wS@99M!KA?ZgP*x$AX-T0UiT*B0i zhP=iWpFe0v{D%xQmF#_pZ*<)DaY!?Y34UntzV!)elz)12jkO91q+A$Rbr?iM9vo8! z?s4er)`wtIvjKEA(~9GewE)QLQsfao??*MWV=DitG=cT8rkNrkgd8d8wGWye1DXQU zt-pU-(c;&UCUd$jz*P_=N>o~d6b23c)lGE)nt;ETs=5==e|gMzjpRE)%H^dVCi8gU zf6EtKu5AVThHE~AwD~~&1vojU-vo#Yd)2>mmH_37?>ftN72tuNK6yu85n!d2xB_2f z0Y2mM$!yzfz^SL7c`GgmxEVdT*P|N`B7GbVvzBdCd0y#}dHihwrV{Y~sjsoS8O#Fcvrn@c0E?4nRwG*a(1|KEA52DB4j+perEemMr(aLJd6teehx~JCE8NjT z^@EdHrkUv2+1@~cUQ4jrgzAoqWumh`Hq2=kBSH0=P2LPgI&$<8pd(wZv(fZWjlXmFzDkyvaf2f3N$0pz{WjOcbxX zCGrsEi`APHel}>HtN{bcujgU z2C!(*Utx-fMxWm%7#Hlt119Pfj^T%Xs4lfEJM8)`u>bN)Ag_fU;;sC!pqrBb%H+EC zLvh}=D4Zm?%uhDCQG(UtqJqdOTYHC zOvx1#;hLYfz5Br9P(nMqrq^sphI#E$R-|jSfgp1XlmK>YGv2qE|@+6pJ+X zfF0*~=OnUf^s$kcFwCkP*gksd<7!udilvg6Z?o3`QJHauR$Pq!MbgZaEH;3Z&{Gdn znln)HD(`Xq`X*2;b#--xuhl)u0<6=? zKNF60s`K~YW#)_b!FCMu^Nr9Zki~1ySfW=6UU_N1b8>0~bve4d>q9l*!x%9;eNY`3 zb|{MfZ?pj<6kAF>jjRT4>x;+9-!%c#5BnX~g_U4L>`bin+(S@8Q9I+&TL!X5H5I8Y zw1PyzYWjDsC4l{egr3Y)8<-*%d?Le@2j19A(~+!nf~lYXN*HvK!Gh}H^-d!M;^OV} zjgUF$RkdUtTkQqH0h}W5O;eF&b*qVIdOzr0G}fJf<)|g`jS zu4}Fx049Nz2dC0Hkhs`SNygj(KuAVG2Hy6e*37b>VlGcWc>dyai2X1+Hzf8S{?rRv zSmVAux%2{EBenUQ-Pi%X-R0%nhgvHAa;mw*a3 zQhEC1lcnbqr0qVPsC(}VxCX=2O-^cr^`5i{@`$eDWFHVJh%3XZr-!zsm((>->m5XoT*vN`3>m zi9_zBvr{PbgQ3a6*Hut2bAiV^`WezQNMC>Z{RePVXGnRj>WBQlZJVA9+XTWaJ5{x= zFM%%&?Jsk)pI|Rv_*sVVM{r`ZH;KsiC%9w(vnNq&85|qEE=*pu3B=;VgKe}nfn&r8 z0TlfMSa3!!hrihYT=BFu!D-7t{r2dd`~ENB8GP2;O@9HLe9}zp6nzNB0z0Sum$E)C90*p0~TqaTtZ0r^ntHCV-907H59jenC_ZTWRTs2q67pz>6V?pJ@B5Z4>>} zKS2FKQv39y-w5}=)gDwh0HGZw384*y`0l{{xAGOcpx}|boW&OsETH7*s6YQ5BpOU^ z`V5fc!D1p-jh1ns+oWlkOh<*qv$nMBm}byg?EcyI{8RYB*yXkCA6qE-=VIAr0yQ2j z&-{GKj1Y^EQ7vp0QR8MptLt|~$gwx~pcMU&Q#k5}^4dfmHRi1AD*E%A63gj6V7wDg zi{*1$nNzNi;rag}ug+ayz|F@V7k}&~#4Q8S=NN35@auTMyFBp|{g=rW61&2JCwAy2 zgC0BtEv6|uKeO5J&EM^25$E@Tq2t642S#=*R;O~!bchV^UgES@W1z;v9aaGs~M9#15IyRpvEm(&2Jjcb%V2TyR^#tT1qf z3g4H$f7(Bp7mnPZH~l_%fRYFPGkf;+Jp3-ESnfPN3$An5-F%uV0BuGmmz4l11l+=L z{EmW9b8zTcF7H`*CAV5(SX~fKvx*ANg|I^b!n=}E90E}0FCVBOeq} zmEYwU7l0b2xv%cUa6`3J!Ab^QVMw{FV!KSu4yl^U@*c#A!eilD;bls6(5tGKO#X`) z?A?~J6)_-$m#4A?Q|iQ_mUO@JB>xV&&u;jIWJDakR^O0W;-kflE*ww21jOO_I5U>| zRt`KU{b)pDR21$T_mMvY{MfkB|KTNdVR-t?Vvn_`2yQmHE%P&)A6}C)>_6iyj%U23 z<2$T4;VZQ_%!7>jWko$Nd>EVpyy>o2gWY6sPq(L_uqF+@ zVElOH`G^eWn0!wtX~%`18c_^x3(Mdu5;V0R?_R*3bMJUJ(=K7{>V$r&VsT7mIl1hz zB7x^_q|VJ8lg8%9>NKU4MX`r<XOzvbc(Zj-mPp&0e z5-BTSnZNB}E?!&cxwY8v@d5>$UvV{QQk4##Xp@}2rlWvwh@THukvI=iNTv^)bVO1Btsvtg0}(GFY=% zT?F2_p{!&vrUDf+dFx-_kb%owU*Da`zXDI}CD6@tC_;1?jmNx|;Eb!!e!7GT{CcoQ zTHgU++GP3P&~9}&{h(GY>$Md8z{D+}dRG&2%W6(sdL#fuKS{WaKJZ+QZSNaEy!D&Eq3Ef^rvGM>GofD0CwEc+%j z;Fqb6>Akv)A9Q1fhv6eEJ!WyjR{t8b7T~1)P>RR zqTDanjBtd~ovoO=Mo?7w8Mqk$UiznB9Z)(>9Pdt;kyX{&hGtTPj>fY6*Xpzf`CevxGm^UmVCF6PPC! z*C~By0|QO#M6rzlJT32akU4J$k(qJ8w+sL%8FOQdnuIEx}?4A`5&uRa25%FpcZ57GD8XPpgT471^TJrzgnbo~kK^GkHFzK}Xb!mKq(chWlh70CWIVD`Ic@39M zWZp`ja>nWW+RaM8u49TXo;sJ^oUq%MFy)&pHn=@tn5v7<0Y7t(-?e&Wk8_iaHCi}q zahk>fW74V<4p=y*+&pcD6-qRR>-Amn1v&fbx@-D)ZRF5zY26KT4*30Cc&&(q!uAth zJ#fbc|1lpd*Hc5)q8+sfUw8aW^_G!`t11)&U)jlY-SE)OG?hbnBbaOdf|*>=1%nxq zo1)LGAYrU%7oCwK9_jv+mj1*Z3hR#Cmn*f!SFJ34ZPT0~#{>rqj5Whwd+RT&Rk^`` ztq(a{+jOx)ZPefKH4k{1CB169K^lh~D_AJMdIPHVklZ}WDGm#NzsL`5xB&+X7X|av zbs+n9nbNOkZorU-%q)T-W-$C}j{kbN2kaA1ZJ7afP-Sr5Pw}TK+_U&c5!2@cDHXQ@ zGH*FSi&|Kpd(REV-Yc+|G_-^5Bg=Z<>O5hflwf}Tt{H4Re`m3M)(fiqGIJVI)`hpN zWu+EwdBfp;ox2Bd;&8=A2d$KN!|nBtfq^=5cso#s__TyKY*03eyB?*Frvi4qsNJ~% z=@w`FjtyGk9FB`Cybs->SjkdB3AY3O`k;I4V9Oa=OU_Dso_4`$GpRjH?)EVMbg9bh ztOxe$AU=%)EFg8jKnn$_H;xgkq?)7DgVTYt@#}m(*s=X)CeN5SOr)bB7#Q-wC!(%U zqq|u4*|SOecR%uO!PF z|B6eJbc*uD*JC_?TErV*m9dxU)(U=@m})NnVxKIgq)!U(^Yz1ey{>lU^0LrZucLkB zj2{NN_o7TE44^)_+_wijzF6#2ry8xT4XjPI5NfaU#^2&i2~dC&RQx7*-Oq1Dbq;v^T$V&BK!?CcG{vi51Tez(IZ7b4zD*!V&h%3+T_S5sWh_;`0|#t&K% zIoI%CQpb-4wl?2C_lK$eZYiJk&p{R(edsF`0LiJ|73A~kz!Vl%LbH8;7<$zsMws6c zdTK}3%nSQNTZ@V}e58&?@qx~YecBg_*D%Wdx90|50HHcc0v|a2v@~yx$s69a`9nr) z=?Mqv3-2lG`9gDwfx)U(XQ(0Q;>mokcbD@V)S5-5p(#{O(i%L#yz#6sKBItExQcs1#Bra#oU zwf~=$lMPl|O|r9~@PlMSf!?;*35QH*5>Ao$KwB>I&b$N<{EY6@dLExAwP}{Gv$^ulvhrFS-0Vf{Rp<35|84n=I2|_@(cYj!LI=F|8$HXX4Ilz9Sy*S3Mvlf zFXds-a|UNcuRzQmmTP%2L>8Z4|3z{wDG;muI$lhfWr**SzGtN748$H{itJ_NwpfQU zv{zX;08_$B9OvYWe*_jbo$K+#b7VPTT9-X>x%!iim2n?@iuJ*gnwSqJ=#x$$BJ{%N z2Pza7Z2d50zM{)|u`8zS(cRu)2{?M6Bu&;!?eXHGmcNNbAm;5SuU4`$!`y4!>}`xe z_-T^99i!Tc(sG@)tIH#Hpf#?m%Z@()15&4AjxGvxL znhCwA?g)gki6m4J$2f63&&gx=C4*pp{L`l)w(3xbtbg&oU=UmnUz7)aX7I{0SI6`B z0wD|29mLJ-0L@opdKlsZV5X??z2q!6DEY92xo+JLem69c(O2+>864fGLl%5sxk@8f z$VXq;v)Nv4Gj;=-JL@OjKN|q)y4H5`7+j&oWj}$|$v`N1!stnVj~)Co#Co=tF&Jv> zgx`7JdmVn#d$dRUIv7H-P>z};B}jSwP0rWf!7%JZLHTfj6iy(`+*Hd5hOuUKI2K4bcZt91%yAk;Q z=yx^pSGG8V`@k!YJreg_tR$dkbi)^UUyWLgL}2#r=aLmvKGv zX4?1~jl;HsR?0c-@YA0gY9aH{xI&?#q%P4IiyeCXRl6OH<)bZ=Gwc*G-zIl+;X)Ml zBz zC#cfF?nW6Gh(kYjWgMS$hY=a`{B!!gcyc@TM0=4pj6G?d$Y9`(7oJNvK0Moa8Q1mI#8*7^{A8Re^1(4}BMBg5foKifjfCTi7mm z=?uAaC?qlUM=uvV;fTufu{ymlcoU0tO=<_gqmtn1JLkh;)cKaZ;r1{{_x`H*>7HF@v^#W-k5dVA-iLImVq4!bpNbqik7OuY9~CIYVas}L?b zCc?ReJli>{2xvnW_r+f`84A>XJ(?Xp+CM#7YK9~Yw(P|Z)!T$aa!#?4WW98FWl{Rf z!B!Z&_iw>MBPIiu(nJUf+zf-v|5zgxb2A~qAFhD{sxTOQUUY;$H48GLe)id_P`G|| zL0!@@8}iPQOfFsvh0A7kPoFU6K+Gs;hS3a9O<^$JcAOyRRu0Thv)hUC2#2@5 z&aX^uWyATyughg10tN&g+Pi~n7*kJGOvDxmy^=mFjCy9lGpv6s5&2E{&N1^lacm}3 z4|2X{x_uLt|9fV^7MTIbqY`ONT%#cM7E@uJV>+x1Ok1ZNi-NDu+-4t?OM@4G{JCs% zJ{q=Hb^knbJO$SEFZO)%h=zV5Z+|?Fx(gQ;qxYTfMMGM(w3c5hw;_3fxB+clG*rvh z|8+|z0XA+%mCF`I!;+5W=Q3@vP<8GTrJPqZyu_6kpQS08iXCDkd5SxZ%kVir7*885}oqjOx)Qm1~W+a4_=^M8zZa{B$e-ndO z;qWovaV9}0C+J#H@)$l1g+i2|`9!WI{N~pwpqCg7`)|K-a?R6$f-UpNf+-MM+bw)V6+T8q(%o?12L5jH{jcY;EuQr7(fGn)4ApxHyy{x)ux?N(&jM0~ z6GAV8)r0JDaY)~%bY~h^UsknxBhw!12JT4+y->ipggd)m&e&s3nd#~UUOgAgE9X2aTG?3af!Gs*< zR4`SB3G|#-(c!EUx&?-A8gQALx$qu!YXymoAt)D7g>a+7gf86*&*%c#&Y$KmvU4=)s5?^3VU!!caq$8h{E$IgXr6BhwX`p;rqf)4xjW+}zrl6bD^TzD zQB?^tNJ$uJCL2O`cGm+(PIY)m<-VZhJw4daej2Tu)`rVtiR@zQS}?R-D2L{jF0`T* za5;XICoLE61^pn`hknFXw|lc7d?b760wt{h9MC3T(y|qU#~I#pqACMOLL(@k@s$o* zR;v|d~-#sIc8CR{pKB!;ht+yDA+QXkrO*#>4fj?Hy1qZs#OwCm(w-Nt!e!A% z8uibvVSV>cDLVp3c_a-? zcscO$b&5hE(adT0%li0)9k*Dz1P7EVCJ}4k(8ocQ?`}FAPM~LIz2w9ldblsO(J7gL z6JH6DEZ`#3!&*&Zgp)#|cxQ`FF_lgiSNbRsyu;FXMK@F6*$Zv#(&Mx6A`fEqt*Tsd zDot!0P;)u=wGv)8{AuZUbRJeA-U4n%6nj zs-LBWzBq$-X~ywn%2oXOhX=Q|!e5YIIC?W_Km+^byXC*nV}~ErqE8&lRmV$<+O|6e zLNMsI*gheZ8t#Cie7~H;;gY_WV!fmae!|Tbbp7WgIQy>4^fJ3LewvhHqtY!0v2RPA z__P8ZHX?Sug8=N+^*+Y>=`wzo8SbkUssOu>N|~=$FX7iWJ&*kpR)k^OIpG0oVmN1l z&e|+N5k^~P&gT*fVJkArJ6lnT@WzD;b1t_y@q;EKIyyc@c(1NdNbuYl9Chx}F|s%X zxT~ooBhdXB?eg6e;#7q2{l}lt0dnMU|2@H$2DLmqxaazQf1eT7T)1cdCrK81^5h0a zAWryiu34Pw`6Z|rv1srokq?qT=6yr{O#(LTwr*EA2*Ndf!>0$EqA;F7D~(+DB0S0W zNn`%c1^DMw`c3v!`JLiNV?{HxYPsiQT1ao(MJ*wjErF5QZJU8ScvZegxsp1d49> z0_^>hvRp^JiF)@7NemtfLK}``3B$icxG+%9=A=JAOwpNlV|SszA3!=A3pF3?uW!w< zT&Kq7&zUCHin-w)mHa_DJcAX*Du~U6IpJm%D=xycxLM-%8JRLR=pD@jixcQD!HeY8 zMph>H9z=#fSx|k|0k)MHSAKDGnGHI~~EyEc;wo@?r zj^AWr)ETVvin5?@j|}E<=ym;-qQRRT?#DDrkHPtamgs{~YV1TEU|J{m7uez;t-`#M zSaA9Hi@ls(z!+|~`?7!n=WPZt#<*+%&-L9wmR2%s#pAwX{AUqx1jy^Y{zHNt1$76W z4t@Zl8{ZatBaUIx8#{LIeBOem)~_-mlnAkkrWNsC>k!bRP5%{C{uj9pmwJ*Fwt>Bm zdeLJhe^Al6nJF8&3J`K#J4NySKH||Ak3CL!3z*f7UC8s@MRU7qx>vu6f>M?0N~ykW z)bzWdJz7x~sB1J9-H6yismce)vMUtOc&Odd%lV1+N0xH@m{Y(9(mR*tb~lmW9Qn-` z968|Ld7bRPNt?(bJ1+?}bnZu^bB|HI>QM00lJ6+vEp6PUcMr1340pQ3u!b(P+00*b z>PL-bj?$xsD`;(=Kcr!A5E;PwGeZ_jD2zO_>6G>;qU|fF9+g={=l1Q4PN+Rc@8nx> z_0DIMyL6kkd+a4jRt(o^DOy026^RB64rA!d{MQMa@DFHD|5tPu#cMQS`;)+8U>cEA z-|3>Wc#Vvo$jDA!e1}TR254tBUL#Y|z~8stjv?Wk2Ma-O#t^X#zL1zQil$OUo-l5` zLaIA1+rbL`C{bt9U@7Va^4-gzIe6HHn5*Jh_}reN5#MuO??2X{+j_Q~>-|Hhp8Lxq z%SY*`QcL=|M|nSbzjNM6U^@aRLxY03gl?p_lbcBVBo}Nwaxi62X+h)L`2iJrB_O9+ z{J~qDN|ZXG>SZZf0;F$9*ox0YqviC7yfX1T(BDMHtvQ_yZe(|)%HFsIeEEF|7PP9t zd1~QZ-Lrm3{tfpNjLvd)2%@BUyfE+Y848*5@l)dupOuwrD|@u zHz9=)HN^|UogkS?sa8e5743vQJLT)z3AirO#*@}`pvE&M>}z2iK)yY6w_d3SeKtIG zswTDVXwD|{1Rqf^TBx}(YW?OBAWWF`9QxXetjn3+O9VB8PC;Q)Y}AKjOMXXN7#{Vf zjR=bQ4EoUAs|QoW!sS4EBu-s!wijJ_p-nO0mkkyjFaG8t>^+)K<{aob6#yLB;}63N zdyu?P_7^Uxcy!?;^NSzZorvOU!`dawe8d*AG?q!(hRiF!n;*9>M}`ZD;*O_U(6O*@ zjQV^v$iec(>S><_#1^Sc{ARTQ8CmBTr7Be*FTvF{iswxzB-259-u7s&Y=fq3V1Ah{6`T%|T)#$*$v#959{+q4c*?*LWB8lODmO^Rxrc*r%?S930Y}pB$(zFZgSMAYfRy2Zgy*k$)?>a!|c;nEtWebShZa3a& zYXk2jsUNmw#oB{s6OMR;NaJE@coiRqw{Dp zy8KA?|f|!n6p?0!pa(vPil16n@)GsE7TPPck7V<1C^z! z>?ma4slc(IQG=?SrL&eKZzHu*vJ<{dRVeXykw@jJ6g0fIYLy~Yfu`x++5M7EN320H z7b`y<&2?sqdrZE`K;^VgAK1E=p#Rj4moaQ-A{V&04NMA9XnhSUY|lbJ>_!5pqH__~ z{^-o}i`mGvZ+yLHCKI8($$yiU+31GnEgw_MG<3g&|7omBHZt)P(5*O?h*p1;(U;9< zAw^shCDR*=o-UcqRPkmZj^w2A&F5jLYuYi;iXju>S1-jxwSCad%;W2)+tShXXNR1X z**1u9wx4|G7_|49y}T9bfkIraUmKD)1KL9rFP|v6Bl5s*!xU0$u=Ae- z%bP7`# zgYAFb_lsqe0GW(n&+h?s(0}dgY;l}A5Q=7uj9}*lX{M)p%igMkWMKy$KGRDGw@xzg zzEB3)2Kw`_gEde<4-pZStt9wC*>ko$N)H_~tv_Rv_pj&Q$%|dxgoa4^{4^VZsT`{5 zQS()EG(_cn;Y`85G|=w`(}b>H`sg}|q(P~(K3eJbvl!viMb01HcABW*h45$Yjz z)L)?PEZ%I2BJ~B_Zip%%L*k~A7pmrHRg7Ec>35flTwrp3Bq~K@vf%A@#st5KB3c`KZ4@6X2ToEGJYV%t zM0vw^7Yq9UkPeH@{!djAy^cOMGZ083BAjR3>O;s z=&1qM#HW(~93)Zfw+*HnXEeZ+?Fl!jl=Gkx#L?^hv^?66jr<_5Ao!27I|+wOo$nLD?}!I`U2c zMz3xvyd#wYdX)s@L$M0LJ2r2|H{&8;lyoOE%TNG+^SqQv8MwioN>N7fBs^;KY0|0P zWdw0gxlN+h<-z_5_RI4`)L^Ex;LT-*%iv>eUmwF4BA`+mm?UEwxrF=PaBuO zT*~2mRu3lR-GFt1$SmT2{W| zEQ;)c)?RzB@*Me7#^dhji=osr4RpdW^nm2td$F}6-zBveduh7xqI_)G-_Hh_BIw@m z=?jOG)M$HY>cxS85UNR^+#(_8Kp)bMhmYUoK}Lk%bJU5>qY7G<-}8}7XjN)NqJ{AS z8tT;?PwON`F}cF5=Xyq1p#h*iYNK-j%6)S>%7?gPo7&8LpAnLfLUg0C19K(6Y z$^*D22XyWbR3*vVa;Eh&uOIgfj&vWX`pX~XxpnD9dk=@^Rv) zR(S>Jfk4MjG{H$SaO$2yJ3|B=h}+(?zVwL}DEO|3=NM3fJI~iw(wdk+9Mjhffg_%e zCqs?zoERH87t80lb;JohKKZmSo{|j+CbQn_G27@#PA7lzUo{ip8LR8%cR%Q{f12OV zqJ0)@5xi==m~h1D;UL|77ETQ|3SEkf-+u3TzC>5pbn}SILFXN&bA4PsM$%R8vjGvH z3G?1U#((6CEaK_BT!_GRNBUa=XGsBn(rWO~by8qp;_|5d9VHO5XcJNSatb`>r+>5f zg$fi(tF`Bl(E+0U+)gn}4w_p@u8Kcr1P7ZWY~V5>c$V|M<(Ul&*u9*ppLp+^{5!41 zzCuAZ5T~^fuGqKI^ZeRrfj_6&Kz!uHygKnQ6jxSM>0!GutX7!u}{*NcTMjtSsV2)?}=8iOgO!UVl>0M^THMJxA`H0Iy;{U>C zYL*Rg+X<@Cu@M37-K7XNO%Bx1XDLC&eSLA4HD<;h-&M432!` z9HM)~2~41Yt=M60iWV?LR{qz{!3$1{`P)s75dv<75x)KrA;9&~ck6lWubu^R?}LqZ z!oVh#MS{bS4sHLMKWox>5wIw*(Of*ij@+IUof+ET2LUk^-i-Ym=ya?Cv*0!-s1Lr= zwd=r!%;gs97XleTL!9D!%Pu-p()(8YV1W{Zdqx&q?IK1U^W;3#--*GN7`RrJ^sdJV zG)BI>Km-I6PF}wqy(v%G;yO_~bx`Xpgo;& zv|*7OENf2NON=oCcm38oM0@PucUYt}&C$9&r;r!;9}PH9l_BSyf5aa!mjCbbt!?=t z`1iPPDK99Yb>#_u^|QxR2UY8Kae|Sa;iZSXe|yMIPT2H)KL>I|QmzhPTzB>4AY+hD>)d5A zFmaW$WB8g9NRrQAGujaVE{^+wOm(Ut{8yNAeDpbRr>n|-VO|5&=a4niTf4YW(?6wfPCAWIQD(S3e-K*$`ir}UZ^JenGBWUW^Qo4HhspAG=%Y?8=6 zF?aFrG2SRt}>v{n<`5_&m0rpMx!XJK70J9qcZk^F~ApVtM!tqQxaNllW`H_S( zAog|Vq5PtTevD-L(w4XZ6tHvhVUZ=J?#dTCPmJa?KmS@RpQD@LnpAUUCMAX z@&+2v|M**x%>rDG)XNLI?v1=;LN3uc@B;ibgmA~$11-x8ii9udpmXGpB4d>7QLW+} z-LX7dWHGGxvy4F#?WTP{F)8hadS!`U`5c{--nUyW$lD32FxaViS=#~$T1wV0bT)`% zw^N>X*cqTkk__HdWu%-)B=~7b+pZt)h_Kll=ymZT!S3&gl#^k3Bq}e#{X>jG3tn z%zJ?M|K!+DSdxH(|IknV9lZ_a$KJg87d4GrjQIKD2Cm#8WPbE-lt`J76Z_5;V2$XMq<>$(tQ=~Alddn>r3igOO~~)eriLS}Q?_d4*9%Gb;`i8juyyP#K^(iu~;(FPzLc(&yBwaNUwq}er>j5E z^WJA4Ffa$*n)!t>VFBRODx68`R|MUznuAsDeqdQ*{b{4q73AhRJ2i6X3o1Jtgg-eL zAe&K_r=mvQfZFj2k*u&Csz`q6V`Ab36nUqoxNDtJL(TG9CZh|8O1{Wd{pkjx;m)u$ zD6|FN<FXu6tnr~(u@x1($VyixeFN$=OZD@S6wg57M}6%pz4k<4dVqdWqY z%2|6mByg^8VzJW!b*cX;BW<~go~|u#-X?QF@LF)56eAa)Ux@5qn0G)z*N^o^eb5I2 z|6FddvDzTU83k9`9V5`JK{=ZKPZwqXsm?r4Wdhi|md852ilI4aHdfXiEx=o9TYOVq z0)!4}&x2EnV5Bu7t2;*#RCP2<{-R<$6SRRdoO6r=})mUGPkBj{U?D=ODXru(0Qu0WdgzouMrA0-E7T7yey`7yfyB=%Wao%Yk*Kl6YdtP0 zklIlzpu%Q`z7>bR63x{Gk{64+8K3K;+y7bnDa7dlTjq6t`XnuM)PnB2Olb-{h}OGD zn#9pvS+l$<8Uv6fb!^Cnfepx$zg^8pH3E#LuJR>;O5lNHkyLT8DmWJ6y@<24K!%Dy zomTf1pqf)|68bhv#-pGYa*TklSm*Kr18|!I!aLyK`GV>Lk)i9y2?D88hIJz5ZC$y zeW3%l{;d=JIje*Oa`fmDY6Pu3wnA(6ztYk+2=w5AY7Ox`d*kHd`}Gdm2SWZ zbSX@phCiYMmpctkG&4vX`D=+MuL>XW;o=l)n@YHVg@JBY2njv1x?vY&ICB9={4k~T zQ{hBA6MiNFuB0GT`=&2p3Og!U390yaq+5-IKQy?r@F43Qgn{363LpR2#IW! zB3qjzDio4dNh+n3rNwrRQrXw+`@ZkXwYz@z`}}@?{brt-d1g7!Gc(W3XFkiB_sbQh zZw%rUhn<~ek(U<9*q@-+ZIvmDy>l6kb9vJ6rNWmjcb<~q3(Zn7sTl#xTscip$gZMrC$LFUoATn zHLXU3egW1N7vs0#8=ZwNhdxTcxG=F#Pd(SrI*0XH)+r0%p+l)Pm!ddMa0%nK+9Lz= zTAxQ)b5()yODO{mk!v6j9|%v#%tYU~;>xxdY=Zo5XZvsNVCkU>Cmc@bX+(=*p0ggC z`QZiG?X$jZ956jg#5a$J2!-{33w!)tLQKx9{_%NasJC6-Yb|shX*s89$ne*O~W3AA%mShp;eA5+^^Cfemz09f6&YbT>0 zXgX8V!wTv_L{MBnT0aktmQ3|*InL78_YA&~7(sy#-==*^J3kBpI`V$I9%RSO4}-;! z(kfu_8O|SD;)7GSW3Re(cY?&~KS}oI$XKm9PA_`bKM*#cMI5KDgIRJ0?M8~dDAal6 zSd#=9@ArAdp&!7~`J6Q3J~u4}yVXVf=k+-8<#AiD_YaxKTt2&8hu96GimdOy`X~Sg z`Nrs^p;c59EYs+GKmd0vzskS(bPCjuIt9y03&2)UhorngD&E+0uR53}h;NcUpE_Sc zflMtI(}p^3Xk|#f8|uT2zw^~`bLsJ76FR@;?iJQ~mmKvzLgRqueCDft-vqIHde1!V zAQNp}J@(+%c0p)U*C6gGL4tKNr4Ox?1#mMD3gti7h62A9wPi*LK;`mhG1H?&E|D2c^xbYe4r=^_6!AIN;ij zFzMF;G^}SPcr0-n8_N%peZp@e05`U9-Z*GA3(|^mJN0=5@Zrji5XqZ;fO2s5krPYT zjhb?3pK{?6N~I^bOeu+AB1vKYAtOQ9->W&~#>Ebq1ID_+-c;`U5WeqW7MEhrhR<>qmkd7M1dH;xYG$+9q3R)~>%Cj7cy&XV z*yYL*#3r8Jf9!)WRD=m>U7|3%@OU1*oii z!N>Cvd&gBumfK(;)N+i&u~r!3eKjAmh_UN83G)t$~m8`SdU{J526fF|6CM z0%XkB7~DsCQ2yst7ni|S)Hnog_ZSeMU!#@v>!l%ZYVU+!C+84S-kH8W+jRy#{bu=# zN?idVt#M!UEr)>VYD$?Hoo{w6gt znX|xHQ80Qej)J2kzUjmnjGzvi=JXg7AxPuo?JYbq2Ueuwepy*T45icIB zqj|Va?Z*;eBrsYf-T2|&rY^0mc67X_Twqjsh>ZPSs{J*&#R)6k%%qb-p!|ItcMKMGi;naD=ZD=6Bn679Qe@PzuvT`-vH;@Kc3At{E&QP zzG9uoEON$n3qi92_`uG6{N1^0z-s%~_~N&G@S$NTmscYd$9*hOp0ee~!57 zJ|R(`t~7=YvYMz>uw$LTP4AC5PJ*cYf1m8RNQR&6pG|x%n?+__J6wW!$T(2*V%rbD zaTXGAox{r|28g@KAzbvJ8F>ZKE(TmA;A7c^1lt1y2!o*^qJkOFd{meB zVZaEAC+Qb&s;EcnKh`fC*c^n6?)Eo)b(jU;-dz3^xupw9*}ktYG3iIY420eAty*B> zU2R_eVF|QNsJws7;xZV#{5W+)fq`zt+BGSjo&s6=DWT1=bF6y8@AWk6ZlsW`q1w^9 zg8oG-SE*zzfVWR_Cmw9=2kg3kx%2H-k))%Cw%>Xd7eTy8|K7S4(7&yhk+xV3TK}}U zW*{o|7UF*LEPnv_`1UD?-e7~(?0w56cHKzc+l1JApNgeS4$}+$3D9_)d1oPo1DeV0 zefda_j2}7%9#AeJ<51@*K}QcZSgPcboYzf*hVR;*jaKvF^Z;5Uha?%#>ojB~-{gdr zNjFkg`xijFc=i6<>-n(%>3=zSvU4bxd{)-y20vW4-FU7hvJQwA@ZLEe%Z2m5kyBE|>UVg9u?n$lGY zbowJ>@#70S?i{?Hn(=^u-L4A#+*VG9Z)APcj97J`oI>8|At->k-rmu!h-pO+n{15! z9pZ$M6V0(wngf70X=LU>CO3ZVB>bC_+6F#Uhn30)Q{k=0%J%~XCeWqyMQ+;yD)y^g zFS@sE1_)Kva%(!W>Px3}gVZTR{3`Hkz0mL?N;qmFhx|z}S?(NNnKNP$=HiQlwg`~^AMPk!wPBI2N`r+>@3%h3GA;%kghGR*HdN@;o0 zjGml!dYNO#MC;yL{29oK0r!&6Dw)O4gJ(xQ17q^4QTd_nok^D%$WBG(M3-+a*dKf( z`h!b9_#k5<&9JBjpYq|!p7v%GWvTI~FQovP6Bl;9OsfJDZ&Tp$eZxSW!sd6~sFQ{6 z-JB~*A4e_mQ_pX5)Sz>rDf=vs&Vr;g!};joPO#&U2s&Q3fDYZh^ykux6X0qnB5rI(KX2A?^*S}W_u9zfa z>QU|J5P3Sb>klb(h@?WJ;+61~UvwDGm}cAflM6=;r_N0@lCWAp^BM6*F1V4L=@of^ z2tQlhH=HeJ=|OLkzbJcr0RaQ`dR{eND6(|P(Rg$aeC_nMi09+RC$E&++#ksYj$Rix z-+jgM+jbmo2~r@}jgwy6+qazr4xZ;SuY8CRQEROlH|%E~_nOifnjNz`nrp5gQaa zH6~b+$ONEzmvjuIV{>|wZooAL+RMg!-_e5u+8oQ?tzxn!h*GEzhsZZc7HlvN_5n!+}R7}xx8bHcTl0}d{D^;?;!-IhWa+N270z+ z9KPU2Vd+R89cj7Ogudq+*f=?}<9hn`J?io-?#$`641Y2m?vlKsJ! zGhsM*lSr`Clg~ENo*jPq@GkhB9y`X#hyF!Lknj~HL;Yw-hjVBA&Pp9#0|BMRM_;*d z;c=A%M8l~WG{mEnzjYTk96hD7=-ku_DqP~k;aP5co!nmY>|6~fRR|7k;pT#~hJNRw zlc&*2n8|zDZW?YiF>~onV1kkIop06eu;L>&yH24`Y`C$t2=`+Gj#8J8NVTEE>`QAm zv$ionzrnlZav?5UwalAQ|5Q5S%P+J@N#=nG zavnpis2eE>D80^R0{&mMbsW5__&~ff*bd7oo)?W z1;KYi3V*v0VC8|`ZkKcjc<p-N$`yZsY|hAXNNe$O)P&oRZZ}kb^#eSy7ja9JtSeh+$GQ7eoS;) zFG9s}D+TU&>_UC^Xa)#`Jr(1BOT}+wbz>>3?Z~1bB6YN#3R9!)Y(EIqAdd#!@Jrt* zxIM4SBeP=?a93?w+U-V$Z$(#z{?Jy?bL&b^jhZ#o7R!@EZ6!gMX==6~eIC$qRRtoV zsQ7tM$|hpf5YkdJ*Wa~{4$1T>>-#bYxVI}uDyeYbYy5F_eKpC}`yXbxuPu-J&zDU{IFxGdC@`W(l`J@jV zkuFYG`$L0EqK5rG&soo3@+}!N87daJYTjZs%LK_|Rh?|1B&gG?eac;mjKj*8UwLS- zu3Iq`YA`B&~mlT=YS{F%|l)YoERMuMC|Jui1FN zgbn}U_7&$%7y+&vDM=MPBv^aKq~jcY9y!VxaPALXK|`X6jcqwhkm*{jut9eQ?2}0F zHUz9ZuVqO|TSgb+*vTA;S|CH=+SO>E)t|^nNj>eW9~<6!J}B|_8Y^#eDbPYIiwuoj zuBBX2m`7DF-bXzuBjV<8yCg?;0#sR4`ZV)#6$BY<|JnDLjO}+E<-a66htdZ+*i1K2 zVVY@)<;I+0Fm9X1@#F$K-pWfUP_U{(m%imm7Dlnd_F|tcJ#n?@#qZMk&U;i`D!vp) zzc~Wd#Sa|Z?8^pk6i+IzsWOl!C%NX#Eh0{f89Qakw+ynjJkxU7+#7W=aVed+n?$g{Kx9vv(^lmDR- z*-vF)t^128MnO_ckh2$@+`@cc5;z0;pWNV;85~8?xv|=zzJtiqC`@sFdJ?2OE33QE zUj_E>dchDoIfdHNO1iVED?zgLkXhxP3GgQP?vrOP29QZA=<`Y*KpLCv{=Lqa23^`s z`y;4LU^c4dxuwzqQdudwqotGtvO|aI{l&|m!r^{bj2ZQjTNaz$9RU z6d@nS6(HD@AnHEC;uU|dGivlZ9kcnx)C`y&)Ry1>E z!(V4UHY(}Q1J$&+Oa_ev(==xaHXU6<7s8d|-(Oip`o$k7o==h>y!I%n^#TL99c=cT z(WT&r1?NH&$_5em?6u9!j|xvq#O4*%7J%!Em+HL(sMsrnlOU1V4D@3YbysvL@FPFR zt01Qd6hH`hg~v>~aau=tCzp|C|cwBH^cE^3vzW z7{K9(;KRRSMA#HO+jNFPz$Z7A+gW9-pe(J9L>>b+D0t$%v0o?yh-d7RPh{1(!Y_2( z6j2^VBUbU0Hz9P`7FDzL2do7Aj~wS@lsNE7-#p1|LOqD=ikHpvrNdMAavjn(j3H%Z z>R-o8?D$52Rw&!IIUs53=1n=s1|lf|~P`r!t)*sEcxBaq{9YLdJ>eGXxfA^RYy2l}8I z@Tx9$nG}LJFAqY=!zr+VU-xa>N+0^*=2_M|HHHL-cx!bmron8J2yfEPJ}~HInkH(z zfUX4gOJ3MkgF4;~e;d2M0{l&wsV*!&C%60q--n2a=1&<qEizH;)2DwX64pRQuk0^bwnab=NT!85l5l&czVt*R z>-?@%ipngLA+hdO$XsC$(C*R>v1Zv@O0A*k^*@;?XYzvgyGk@1Xw5_fZMUV&AGm}{TOZUtbDoGB5^5O4~eQ&J+ zyYPRVpA;Er{pGt>2WFO$>0u&Ch%7|--$9B%5qnk;s7)aJH_`rY>wk+t;(v2**KnPSl|L?&68v=>h+602LtCPo>|7}L0{