From 9f20c7fd55ef6786fe09242c0bbf72af6d491369 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Wed, 2 Jul 2025 20:06:01 -0500 Subject: [PATCH 1/7] Draft --- .../solvers/GravityLinearOpSolver.py | 189 +++++++++++++ .../pygeos_tools/solvers/GravitySolver.py | 169 ++++++++++++ .../pygeos_tools/solvers/InversionUtils.py | 255 ++++++++++++++++++ .../src/geos/pygeos_tools/solvers/Solver.py | 6 +- .../src/geos/pygeos_tools/solvers/__init__.py | 2 + 5 files changed, 620 insertions(+), 1 deletion(-) create mode 100644 pygeos-tools/src/geos/pygeos_tools/solvers/GravityLinearOpSolver.py create mode 100644 pygeos-tools/src/geos/pygeos_tools/solvers/GravitySolver.py create mode 100644 pygeos-tools/src/geos/pygeos_tools/solvers/InversionUtils.py diff --git a/pygeos-tools/src/geos/pygeos_tools/solvers/GravityLinearOpSolver.py b/pygeos-tools/src/geos/pygeos_tools/solvers/GravityLinearOpSolver.py new file mode 100644 index 00000000..6be54e3c --- /dev/null +++ b/pygeos-tools/src/geos/pygeos_tools/solvers/GravityLinearOpSolver.py @@ -0,0 +1,189 @@ +# ------------------------------------------------------------------------------------------------------------ +# SPDX-License-Identifier: LGPL-2.1-only +# +# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC +# Copyright (c) 2018-2024 TotalEnergies +# Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University +# Copyright (c) 2023-2024 Chevron +# Copyright (c) 2019- GEOS/GEOSX Contributors +# All rights reserved +# +# See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. +# ------------------------------------------------------------------------------------------------------------ + +import numpy as np +from typing import Optional, Tuple +from scipy.sparse.linalg import LinearOperator + +from geos.pygeos_tools.solvers.InversionUtils import computeResidual, computeL2Loss +from geos.pygeos_tools.solvers.GravitySolver import GravitySolver + + +class GravityLinearOpSolver( LinearOperator ): + """ + A SciPy-compatible linear operator that wraps GEOS-based forward and adjoint gravity modeling. + + This class provides a bridge between GEOS's GravitySolver and SciPy's optimization and inversion + tools by exposing the forward and adjoint operations as a LinearOperator. It supports: + + - Forward modeling via `_matvec` + - Adjoint modeling via `_rmatvec` + - Loss and gradient computation for inversion workflows + + Parameters + ---------- + rank : int, optional + MPI rank of the current process (default is 0). + scaleData : float, optional + Scaling factor applied to the forward model output (default is 1.0). + xml : object, optional + Parsed XML configuration for GEOS. + nm : int + Number of model parameters (must be > 0). + nd : int, optional + Number of data points. If not provided, it will be inferred on first forward call. + dtype : np.dtype, optional + Data type for the operator (default is np.float64). + """ + + def __init__( self, + rank: int = 0, + scaleData: float = 1.0, + xml: Optional[ object ] = None, + nm: int = 0, + nd: Optional[ int ] = None, + dtype: np.dtype = np.float64 ) -> None: + + if nm <= 0: + raise ValueError( f"Invalid dimensions: nm={nm}, must be a positive integer." ) + + self.rank = rank + self.scaleData = scaleData + self.xml = xml + self.nm = nm + self.nd = nd + self.dtype = dtype + self._shape_initialized = False + self.solver = GravitySolver( solverType="GravityFE" ) + + if self.nd is not None: + super().__init__( dtype=self.dtype, shape=( self.nd, self.nm ) ) + self._shape_initialized = True + + def _initialize_shape( self, nd: int ) -> None: + if not self._shape_initialized: + super().__init__( dtype=self.dtype, shape=( nd, self.nm ) ) + self._shape_initialized = True + + def _resetSolver( self ) -> None: + self.solver.initialize( rank=self.rank, xml=self.xml ) + + def _matvec( self, x: np.ndarray ) -> np.ndarray: + self._resetSolver() + y = self.solver.modeling( x, scale_data=self.scaleData ) + self._initialize_shape( nd=y.size ) + return y + + def _rmatvec( self, y: np.ndarray ) -> np.ndarray: + self._resetSolver() + self._initialize_shape( nd=y.size ) + return self.solver.adjoint( self.nm, y ) + + def getData( self, model: np.ndarray ) -> np.ndarray: + """ + Run forward modeling on the given model. + + Parameters + ---------- + model : np.ndarray + Model parameters. + + Returns + ------- + np.ndarray + Predicted data. + """ + return self._matvec( model ) + + def getLoss( self, model: np.ndarray, dObs: np.ndarray, scale: float = 1.0 ) -> float: + """ + Compute the L2 loss between predicted and observed data. + + Parameters + ---------- + model : np.ndarray + Model parameters. + dObs : np.ndarray + Observed data. + scale : float, optional + Scaling factor for the loss (default is 1.0). + + Returns + ------- + float + Scaled L2 loss. + """ + dPred = self._matvec( model ) + return computeL2Loss( dPred, dObs, scale ) + + def getGradient( self, model: np.ndarray, dObs: np.ndarray, scale: float = 1.0 ) -> np.ndarray: + """ + Compute the gradient of the loss with respect to the model. + + Parameters + ---------- + model : np.ndarray + Model parameters. + dObs : np.ndarray + Observed data. + scale : float, optional + Scaling factor for the gradient (default is 1.0). + + Returns + ------- + np.ndarray + Gradient vector. + """ + dPred = self._matvec( model ) + residue = computeResidual( dPred, dObs ) + return self._rmatvec( residue ) * scale + + def getLossAndGradient( self, + model: np.ndarray, + dObs: np.ndarray, + scale: float = 1.0 ) -> Tuple[ float, np.ndarray ]: + """ + Compute both the loss and its gradient. + + Parameters + ---------- + model : np.ndarray + Model parameters. + dObs : np.ndarray + Observed data. + scale : float, optional + Scaling factor (default is 1.0). + + Returns + ------- + Tuple[float, np.ndarray] + Tuple containing the loss and the gradient. + """ + dPred = self._matvec( model ) + residue = computeResidual( dPred, dObs ) + loss = computeL2Loss( dPred, dObs, scale ) + gradient = self._rmatvec( residue ) * scale + return loss, gradient + + def finalize( self ) -> None: + """ + Finalize the GEOS solver if it was initialized. + """ + if self.solver.alreadyInitialized: + self.solver.finalize() + + def __enter__( self ) -> "GravityLinearOpSolver": + return self + + def __exit__( self, excType, excValue, traceback ) -> None: + self.finalize() diff --git a/pygeos-tools/src/geos/pygeos_tools/solvers/GravitySolver.py b/pygeos-tools/src/geos/pygeos_tools/solvers/GravitySolver.py new file mode 100644 index 00000000..1fb6bfdc --- /dev/null +++ b/pygeos-tools/src/geos/pygeos_tools/solvers/GravitySolver.py @@ -0,0 +1,169 @@ +# ------------------------------------------------------------------------------------------------------------ +# SPDX-License-Identifier: LGPL-2.1-only +# +# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC +# Copyright (c) 2018-2024 TotalEnergies +# Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University +# Copyright (c) 2023-2024 Chevron +# Copyright (c) 2019- GEOS/GEOSX Contributors +# All rights reserved +# +# See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. +# ------------------------------------------------------------------------------------------------------------ + +import numpy as np +import numpy.typing as npt +from typing import Dict, List +from typing_extensions import Self +from geos.pygeos_tools.solvers.Solver import Solver +from geos.pygeos_tools.wrapper import get_matching_wrapper_path, set_wrapper_to_value, allgather_wrapper +from pygeosx import run, COMPLETED +import copy + +__doc__ = """ +GravitySolver class inherits from Solver class. + +This adds accessor methods for gravity modeling and adjoint. +""" + + +class GravitySolver( Solver ): + """ + Gravity solver object containing methods to run gravity simulations in GEOS + + Attributes + ----------- + The ones inherited from Solver class + """ + + def __init__( self: Self, solverType: str = "GravityFE", **kwargs ): + """ + Parameters + ----------- + solverType : str + The solver used in the XML file + """ + super().__init__( solverType, **kwargs ) + + def getGzAtStations( self: Self ) -> npt.NDArray: + """ + Get gz values at station coordinates + + Returns + ------ + numpy Array : Array containing gz values at all stations coordinates + """ + return self.getGeosWrapperByName( "gzAtStations", [ "Solvers" ] ) + + def getDensityModel( self: Self, filterGhost: bool = False, **kwargs ) -> npt.NDArray: + """ + Get the density values + WARNING: this function aims to work in the specific case of having only 1 CellElementRegion in your XML file + and that this CellElementRegion contains only one cellBlock. + + Parameters + ----------- + densityName : str + Name of density array in GEOS + filterGhost : bool + Filter the ghost ranks + + Returns + ------- + Numpy Array : Array containing the density values + """ + density = self.getSolverFieldWithPrefix( "mediumDensity", **kwargs ) + + if density is not None: + if filterGhost: + density_filtered = self.filterGhostRankFor1RegionWith1CellBlock( density, **kwargs ) + if density_filtered is not None: + return density_filtered + else: + print( "getDensityModel->filterGhostRank: No ghostRank was found.", flush=True ) + else: + return density + else: + print( "getDensityModel: No velocity was found.", flush=True ) + return None + + def updateDensityModel( self: Self, density: npt.NDArray ) -> None: + """ + Update density values in GEOS + + Parameters + ----------- + density : array + New values for the density + """ + self.setGeosWrapperValueByName( "mediumDensity", value=density, filters=[ self.discretization ] ) + + def setMode( self: Self, mode: str ) -> None: + mode_key = get_matching_wrapper_path( self.geosx, [ 'mode' ] ) + self.geosx.get_wrapper( mode_key ).set_value( mode ) + + def modeling( self: Self, model: npt.NDArray, scale_data: float = 1.0 ) -> npt.NDArray: + + # Make sure we are in modeling mode. + self.setMode( "modeling" ) + self.applyInitialConditions() + + # Send model. + self.updateDensityModel( model ) + + # Run. + while run() != COMPLETED: + pass + + # Get vertical component gz at all stations. + gz = copy.deepcopy( self.getGzAtStations() ) * scale_data + + return gz + + def adjoint( self: Self, nm: int, residue: npt.NDArray ) -> npt.NDArray: + + # Make sure we are in adjoint mode. + self.setMode( "adjoint" ) + self.applyInitialConditions() + + # Send the residues to Geos. + self.setGeosWrapperValueByName( "residue", value=residue ) + + # Run + while run() != COMPLETED: + pass + + # Retrieve adjoint: + # * If any, gather all subdomains. + # * Make sure to remove ghosts. + adjoint = np.zeros( nm ) + + try: + localToGlobal_key = get_matching_wrapper_path( + self.geosx, [ self.discretization, 'elementSubRegions', 'localToGlobalMap' ] ) + ghost_key = get_matching_wrapper_path( self.geosx, + [ self.discretization, 'elementSubRegions', 'ghostRank' ] ) + adjoint_key = get_matching_wrapper_path( self.geosx, + [ self.discretization, 'elementSubRegions', 'adjoint' ] ) + except Exception as e: + raise RuntimeError( f"Failed to resolve wrapper paths: {e}" ) + + adjoint_field = copy.deepcopy( allgather_wrapper( self.geosx, adjoint_key, ghost_key=ghost_key ) ) + + localToGlobal_noGhost = allgather_wrapper( self.geosx, localToGlobal_key, ghost_key=ghost_key ).astype( int ) + adjoint[ localToGlobal_noGhost ] = adjoint_field + + #print(f'size adjoint {adjoint.shape}', flush=True) + #localToGlobal_noGhost = self.getLocalToGlobalMapFor1RegionWith1CellBlock(filterGhost=True) + #print(f'size localToGlobal_noGhost {localToGlobal_noGhost.shape} {np.min(localToGlobal_noGhost)} {np.max(localToGlobal_noGhost)}') + + # adjoint_field = self.getAdjoint() + #adjoint_field = self.getSolverFieldWithPrefix("adjoint") + #adjoint_key = get_matching_wrapper_path(self.geosx, [self.discretization, 'elementSubRegions', 'adjoint']) + #adjoint_field = allgather_wrapper(self.geosx, adjoint_key) + + #print(f'size adjoint_field {adjoint_field.shape} {np.min(adjoint_field)} {np.max(adjoint_field)}', flush=True) + + #adjoint[localToGlobal_noGhost] = adjoint_field + + return adjoint diff --git a/pygeos-tools/src/geos/pygeos_tools/solvers/InversionUtils.py b/pygeos-tools/src/geos/pygeos_tools/solvers/InversionUtils.py new file mode 100644 index 00000000..6a0c795e --- /dev/null +++ b/pygeos-tools/src/geos/pygeos_tools/solvers/InversionUtils.py @@ -0,0 +1,255 @@ +# ------------------------------------------------------------------------------------------------------------ +# SPDX-License-Identifier: LGPL-2.1-only +# +# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC +# Copyright (c) 2018-2024 TotalEnergies +# Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University +# Copyright (c) 2023-2024 Chevron +# Copyright (c) 2019- GEOS/GEOSX Contributors +# All rights reserved +# +# See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. +# ------------------------------------------------------------------------------------------------------------ + +import numpy as np +import importlib.util +import warnings +from math import isclose +from typing import Callable, Optional, Tuple + + +def computeResidual( predicted: np.ndarray, observed: np.ndarray ) -> np.ndarray: + """ + Compute the residual vector between predicted and observed data. + + Parameters + ---------- + predicted : np.ndarray + Predicted data vector. + observed : np.ndarray + Observed data vector. + + Returns + ------- + np.ndarray + Residual vector (predicted - observed). + """ + return predicted - observed + + +def computeL2Loss( predicted: np.ndarray, observed: np.ndarray, scale: float = 1.0 ) -> float: + """ + Compute the L2 loss (squared L2 norm) between predicted and observed data. + + Parameters + ---------- + predicted : np.ndarray + Predicted data vector. + observed : np.ndarray + Observed data vector. + scale : float, optional + Scaling factor applied to the loss (default is 1.0). + + Returns + ------- + float + Scaled L2 loss value. + """ + residual = computeResidual( predicted, observed ) + loss = 0.5 * np.linalg.norm( residual )**2 + + return scale * loss + + +def gradientTest( f: Callable[ [ np.ndarray ], float ], + g: Callable[ [ np.ndarray ], np.ndarray ], + m0: np.ndarray, + dm: np.ndarray, + flag_verbose: bool = True, + nk: int = 25, + dk: float = 1.8, + history_filename: Optional[ str ] = None ) -> Tuple[ bool, np.ndarray, np.ndarray ]: + r""" + Tests that the gradient g of an objective function f satisfies the Taylor expansion: + f(m0 + h·Δm) ≈ f(m0) + h·g(m0)ᵀ·Δm + O(h²) + + Parameters + ---------- + f : function + Objective function of the form m->f(m). + g : function + Gradient function of the form m->g(m). + m0 : :obj:`np.ndarray` + Reference model. + dm : :obj:`np.ndarray` + Model perturbation. + flag_verbose : :obj:`logical`, optional + Print convergence report. + nk : :obj:`int`, optional + Number of steps (values of h). + dk : :obj:`float`, optional + Ratio between two consecutive steps. + + Returns + ------- + passed : :obj:`logical` + True if the gradient test passed, False otherwise. + slope : :obj:`np.ndarray` + Slope of the error curve at every step (should be 2. if the jacobian is correct). + history: :obj:`np.ndarray` + Full convergence history (steps and errors), useful for further display. + + """ + f0 = f( m0 ) + g0_dm = np.dot( g( m0 ), dm ) + + e0 = np.zeros( nk ) + e1 = np.zeros( nk ) + h = np.zeros( nk ) + + for k in range( nk ): + if flag_verbose: + print( f"Gradient test: Refinement {k} / {nk}", flush=True ) + + h[ k ] = pow( dk, -( k + 1 ) ) + f1 = f( m0 + h[ k ] * dm ) + + e0[ k ] = abs( f1 - f0 ) / abs( f0 ) + e1[ k ] = abs( f1 - f0 - h[ k ] * g0_dm ) / abs( f0 ) + + if flag_verbose: + print( " h h2 e0 e1", flush=True ) + for k in range( nk ): + print( f'{h[k]:.4e} {h[k]*h[k]:.4e} {e0[k]:.4e} {e1[k]:.4e}', flush=True ) + + # Keep meaningful part of the curve, and estimate its slope... + # Should be 2 since quadratic behaviour is expected. + ind = np.where( ( 1e-15 < e1 ) & ( e1 < 1e-2 ) )[ 0 ] + if ind.size > 1: + slope = np.diff( np.log( e1[ ind ] ) ) / np.diff( np.log( h[ ind ] ) ) + mean_slope = np.mean( slope ) + passed = ( isclose( mean_slope, 2.0, abs_tol=0.1 ) or np.count_nonzero( slope > 1.9 ) > 2 + or np.all( e1 < 1e-15 ) ) + else: + slope = np.array( [] ) + passed = False + + history = np.vstack( ( h, h * h, e0, e1 ) ) + + if history_filename is not None: + np.savetxt( history_filename, history ) + + return passed, slope, history + + +def plotGradientTest( h: np.ndarray, + e1: np.ndarray, + save_path: Optional[ str ] = None, + style: Optional[ dict ] = None ) -> Optional[ "matplotlib.figure.Figure" ]: + """ + Plot the gradient test error curve and reference trends. + + This function visualizes the convergence behavior of a gradient test by plotting + the first-order Taylor error (e1) against the step size (h) on a log-log scale. + It also overlays reference lines for linear and quadratic convergence trends. + + Parameters + ---------- + h : np.ndarray + Array of step sizes used in the gradient test. + e1 : np.ndarray + Array of first-order Taylor error values. + save_path : str, optional + If provided, the plot will be saved to this file path (e.g., "plot.png"). + style : dict, optional + Dictionary of matplotlib style arguments (e.g., color, marker, linestyle) to customize the error curve. + + Returns + ------- + matplotlib.figure.Figure or None + The matplotlib Figure object if plotting is successful, otherwise None if matplotlib is not available or data is invalid. + """ + + if importlib.util.find_spec( "matplotlib" ) is None: + warnings.warn( "matplotlib is not installed." ) + return None + + import matplotlib + matplotlib.use( "Agg" ) # Safe for headless environments + import matplotlib.pyplot as plt + + # Filter out invalid values + h = np.asarray( h ) + e1 = np.asarray( e1 ) + valid = np.isfinite( h ) & np.isfinite( e1 ) & ( h > 0 ) & ( e1 > 0 ) + + if not np.any( valid ): + warnings.warn( "No valid data points for plotting." ) + return None + + h_valid = h[ valid ] + e1_valid = e1[ valid ] + + try: + fig, ax = plt.subplots() + except Exception as e: + print( f"Subplot creation failed: {e}" ) + fig = plt.figure() + ax = fig.add_subplot( 1, 1, 1 ) + + style = style or {} + try: + ax.loglog( h_valid, e1_valid, label="Error", **style ) + ax.loglog( h_valid, h_valid, label="Linear trend", linestyle='--' ) + ax.loglog( h_valid, np.power( h_valid, 2 ), label="Quadratic trend", linestyle=':' ) + except Exception as e: + print( f"Plotting failed: {e}" ) + return None + + ax.grid( True ) + ax.legend() + ax.set_title( "Gradient test" ) + ax.set_xlabel( "h" ) + ax.set_ylabel( "Error" ) + + if save_path: + try: + fig.savefig( save_path, bbox_inches='tight' ) + print( f"Plot saved to {save_path}" ) + except Exception as e: + print( f"Saving plot failed: {e}" ) + return None + + return fig + + +def plotGradientTestFromFile( filename: str, + column_h: int = 0, + column_error: int = 3, + **kwargs ) -> Optional[ "matplotlib.figure.Figure" ]: + """ + Load gradient test history from a file and plot the error curve. + + This function reads a saved gradient test history file (typically generated by `gradientTest`) + and plots the first-order Taylor error against the step size using `plotGradientTest`. + + Parameters + ---------- + filename : str + Path to the file containing the gradient test history (as a 2D array). + column_h : int, optional + Index of the column containing step sizes (default is 0). + column_error : int, optional + Index of the column containing first-order Taylor errors (default is 3). + **kwargs : dict + Additional keyword arguments passed to `plotGradientTest`, such as `style` or `save_path`. + + Returns + ------- + matplotlib.figure.Figure or None + The matplotlib Figure object if plotting is successful, otherwise None if matplotlib is not available. + """ + history = np.loadtxt( filename ) + h = history[ column_h, : ] + e1 = history[ column_error, : ] + return plotGradientTest( h, e1, **kwargs ) diff --git a/pygeos-tools/src/geos/pygeos_tools/solvers/Solver.py b/pygeos-tools/src/geos/pygeos_tools/solvers/Solver.py index e562836d..c28ff3a8 100644 --- a/pygeos-tools/src/geos/pygeos_tools/solvers/Solver.py +++ b/pygeos-tools/src/geos/pygeos_tools/solvers/Solver.py @@ -143,8 +143,8 @@ def initialize( self: Self, rank: int = 0, xml: XML = None ) -> None: raise ValueError( f"The solver type '{self.type}' does not exist in your XML '{self.xml.filename}'." ) + geosState: int = self._getGEOSState() if not self.alreadyInitialized: - geosState: int = self._getGEOSState() if geosState == GEOS_STATE.UNINITIALIZED.value: self.geosx = pygeosx.initialize( rank, self.geosxArgs.getCommandLine() ) self.alreadyInitialized = True @@ -168,6 +168,10 @@ def initialize( self: Self, rank: int = 0, xml: XML = None ) -> None: self.updateOutputs() self.updateTimeVariables() + else: + if geosState == GEOS_STATE.COMPLETED.value: + self.geosx = pygeosx.reinit( self.geosxArgs.getCommandLine() ) + """ Accessors from pygeosx and xml """ diff --git a/pygeos-tools/src/geos/pygeos_tools/solvers/__init__.py b/pygeos-tools/src/geos/pygeos_tools/solvers/__init__.py index 6918bc50..b28890ce 100644 --- a/pygeos-tools/src/geos/pygeos_tools/solvers/__init__.py +++ b/pygeos-tools/src/geos/pygeos_tools/solvers/__init__.py @@ -18,3 +18,5 @@ from geos.pygeos_tools.solvers.ReservoirSolver import ReservoirSolver from geos.pygeos_tools.solvers.Solver import Solver from geos.pygeos_tools.solvers.WaveSolver import WaveSolver +from geos.pygeos_tools.solvers.GravitySolver import GravitySolver +from geos.pygeos_tools.solvers.GravityLinearOpSolver import GravityLinearOpSolver From 12a4adb90534d5c95529ff16405048adcc2eeaae Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Wed, 2 Jul 2025 20:09:29 -0500 Subject: [PATCH 2/7] added example --- .../adjoint-test/gravityAdjointTest.py | 98 +++++++++++++++++ .../gradient-test/gravityGradientTest.py | 101 ++++++++++++++++++ .../gravity/modeling/gravityModeling.py | 62 +++++++++++ .../modeling/gravityModelingLinearOp.py | 63 +++++++++++ 4 files changed, 324 insertions(+) create mode 100644 pygeos-tools/examples/solvers/gravity/adjoint-test/gravityAdjointTest.py create mode 100644 pygeos-tools/examples/solvers/gravity/gradient-test/gravityGradientTest.py create mode 100644 pygeos-tools/examples/solvers/gravity/modeling/gravityModeling.py create mode 100644 pygeos-tools/examples/solvers/gravity/modeling/gravityModelingLinearOp.py diff --git a/pygeos-tools/examples/solvers/gravity/adjoint-test/gravityAdjointTest.py b/pygeos-tools/examples/solvers/gravity/adjoint-test/gravityAdjointTest.py new file mode 100644 index 00000000..302b12a3 --- /dev/null +++ b/pygeos-tools/examples/solvers/gravity/adjoint-test/gravityAdjointTest.py @@ -0,0 +1,98 @@ +# ------------------------------------------------------------------------------------------------------------ +# SPDX-License-Identifier: LGPL-2.1-only +# +# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC +# Copyright (c) 2018-2024 TotalEnergies +# Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University +# Copyright (c) 2023-2024 Chevron +# Copyright (c) 2019- GEOS/GEOSX Contributors +# Copyright (c) 2019- INRIA project-team Makutu +# All rights reserved +# +# See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. +# ------------------------------------------------------------------------------------------------------------ + +import argparse +import numpy as np +import sys +from mpi4py import MPI + +from geos.pygeos_tools.input import XML +from geos.pygeos_tools.solvers import GravityLinearOpSolver + +__doc__ = """ +This is an example of how to set up and run an adjoint test for GEOS simulation using the GravitySolver. +""" + + +def parse_args(): + """ + Parse command-line arguments. + + Returns + ------- + argparse.Namespace + Parsed arguments including XML input, dimensions and scaling factor. + """ + parser = argparse.ArgumentParser( description="Run a GEOS gravity adjoint test." ) + parser.add_argument( '--xml', type=str, required=True, help="Input XML file for GEOSX." ) + parser.add_argument( '--n_model', type=int, required=True, help="Number of model samples." ) + parser.add_argument( '--n_data', type=int, required=True, help="Number of data samples." ) + parser.add_argument( '--scale', type=float, default=1.0, help="Scaling factor." ) + return parser.parse_args() + + +def main(): + rank = MPI.COMM_WORLD.Get_rank() + verbose = ( rank == 0 ) + + args = parse_args() + xml_file = args.xml + scale = args.scale + nm = args.n_model + nd = args.n_data + + # Load XML configuration + xml = XML( xml_file ) + + if rank == 0: + x = np.random.rand( nm ) + y = np.random.rand( nd ) + else: + x = None + y = None + + x = MPI.COMM_WORLD.bcast( x, root=0 ) + y = MPI.COMM_WORLD.bcast( y, root=0 ) + + # Initialize solver as a linear operator + solver = GravityLinearOpSolver( rank=rank, xml=xml, nm=args.n_model, nd=args.n_data, scaleData=scale ) + + Ax = solver._matvec( x ) + ATy = solver._rmatvec( y ) + + IP1 = np.dot( Ax.T, y ) + IP2 = np.dot( x.T, ATy ) + + nAx = np.linalg.norm( Ax, 2 ) + nx = np.linalg.norm( x, 2 ) + nATy = np.linalg.norm( ATy, 2 ) + ny = np.linalg.norm( y, 2 ) + + e = abs( IP1 - IP2 ) / np.max( [ nAx * ny, nATy * nx ] ) + + if verbose: + print( "\n=== Adjoint Test Summary ===" ) + print( f"IP1 = {IP1}" ) + print( f"IP2 = {IP2}" ) + print( f"Passed: {e < 1e-13}" ) + print( f"Error: {e}" ) + + solver.finalize() + + if e >= 1e-13: + sys.exit( 1 ) + + +if __name__ == "__main__": + main() diff --git a/pygeos-tools/examples/solvers/gravity/gradient-test/gravityGradientTest.py b/pygeos-tools/examples/solvers/gravity/gradient-test/gravityGradientTest.py new file mode 100644 index 00000000..f14a1029 --- /dev/null +++ b/pygeos-tools/examples/solvers/gravity/gradient-test/gravityGradientTest.py @@ -0,0 +1,101 @@ +# ------------------------------------------------------------------------------------------------------------ +# SPDX-License-Identifier: LGPL-2.1-only +# +# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC +# Copyright (c) 2018-2024 TotalEnergies +# Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University +# Copyright (c) 2023-2024 Chevron +# Copyright (c) 2019- GEOS/GEOSX Contributors +# Copyright (c) 2019- INRIA project-team Makutu +# All rights reserved +# +# See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. +# ------------------------------------------------------------------------------------------------------------ + +import argparse +import numpy as np +import sys +from mpi4py import MPI + +from geos.pygeos_tools.input import XML +from geos.pygeos_tools.solvers import GravityLinearOpSolver +from geos.pygeos_tools.solvers.InversionUtils import gradientTest, plotGradientTestFromFile + +__doc__ = """ +This is an example of how to set up and run a gradient test for GEOS simulation using the GravitySolver. +""" + + +def parse_args(): + """ + Parse command-line arguments. + + Returns + ------- + argparse.Namespace + Parsed arguments including XML input, models, and scaling factor. + """ + parser = argparse.ArgumentParser( description="Run a GEOS gravity simulation and gradient test." ) + parser.add_argument( '--xml', type=str, required=True, help="Input XML file for GEOSX." ) + parser.add_argument( '--m_true', type=str, required=True, help="True model file (.npy)." ) + parser.add_argument( '--m0', type=str, required=True, help="Initial model file (.npy)." ) + parser.add_argument( '--dm', type=str, required=True, help="Perturbation model file (.npy)." ) + parser.add_argument( '--scale', type=float, default=1.0, help="Scaling factor for the loss." ) + parser.add_argument( '--history', type=str, default='grad_test.txt', help="Output file for gradient test history." ) + parser.add_argument( '--plot', action='store_true', help="Plot the gradient test results." ) + return parser.parse_args() + + +def main(): + + rank = MPI.COMM_WORLD.Get_rank() + verbose = ( rank == 0 ) + + args = parse_args() + xml_file = args.xml + scale = args.scale + + # Load XML configuration + xml = XML( xml_file ) + + # Load models from .npy files + try: + m_true = np.load( args.m_true ) + m0 = np.load( args.m0 ) + dm = np.load( args.dm ) + except FileNotFoundError as e: + raise RuntimeError( f"Could not load input model file: {e}" ) + + # Initialize solver as a linear operator + solver = GravityLinearOpSolver( rank=rank, xml=xml, nm=m_true.size, scaleData=scale ) + + # Generate observed data + uobs = solver.getData( m_true ) + + # Run gradient test + passed, slope, history = gradientTest( lambda m: solver.getLoss( m, uobs, scale=scale ), + lambda m: solver.getGradient( m, uobs, scale=scale ), + m0, + dm, + history_filename=args.history, + flag_verbose=verbose ) + + if verbose: + print( "\n=== Gradient Test Summary ===" ) + print( f"Passed: {passed}" ) + print( f"Slope: {slope if slope.size > 0 else 'N/A'}" ) + print( f"History saved to: {args.history}" ) + + if args.plot: + plotGradientTestFromFile( args.history, save_path="grad_test.png", style={ "color": "red", "marker": "o" } ) + + # Finalize GEOS + solver.finalize() + + # Exit with non-zero code if test failed + if not passed: + sys.exit( 1 ) + + +if __name__ == "__main__": + main() diff --git a/pygeos-tools/examples/solvers/gravity/modeling/gravityModeling.py b/pygeos-tools/examples/solvers/gravity/modeling/gravityModeling.py new file mode 100644 index 00000000..169e30b1 --- /dev/null +++ b/pygeos-tools/examples/solvers/gravity/modeling/gravityModeling.py @@ -0,0 +1,62 @@ +# ------------------------------------------------------------------------------------------------------------ +# SPDX-License-Identifier: LGPL-2.1-only +# +# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC +# Copyright (c) 2018-2024 TotalEnergies +# Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University +# Copyright (c) 2023-2024 Chevron +# Copyright (c) 2019- GEOS/GEOSX Contributors +# All rights reserved +# +# See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. +# ------------------------------------------------------------------------------------------------------------ +import argparse +import copy + +from mpi4py import MPI + +import numpy as np +from geos.pygeos_tools.input import XML +from geos.pygeos_tools.solvers import GravitySolver +from pygeosx import run, COMPLETED + +__doc__ = """ +“This is an example of how to set up and run your GEOS simulation using the GravitySolver. +""" + + +def parse_args(): + """Get arguments + + Returns: + argument '--xml': Input xml file for GEOSX + """ + parser = argparse.ArgumentParser( description="Gravity simulation example" ) + parser.add_argument( '--xml', type=str, required=True, help="Input xml file for GEOS" ) + parser.add_argument( '--m_true', type=str, required=True, help="True model (.npy)" ) + + args, _ = parser.parse_known_args() + return args + + +def main(): + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + + args = parse_args() + xmlfile = args.xml + xml = XML( xmlfile ) + + solver = GravitySolver() + solver.initialize( rank=rank, xml=xml ) + + density = np.load( args.m_true ) + + gz = solver.modeling( density ) + print( f"gz min={np.min(gz)}, max={np.max(gz)}" ) + + solver.finalize() + + +if __name__ == "__main__": + main() diff --git a/pygeos-tools/examples/solvers/gravity/modeling/gravityModelingLinearOp.py b/pygeos-tools/examples/solvers/gravity/modeling/gravityModelingLinearOp.py new file mode 100644 index 00000000..304de08e --- /dev/null +++ b/pygeos-tools/examples/solvers/gravity/modeling/gravityModelingLinearOp.py @@ -0,0 +1,63 @@ +# ------------------------------------------------------------------------------------------------------------ +# SPDX-License-Identifier: LGPL-2.1-only +# +# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC +# Copyright (c) 2018-2024 TotalEnergies +# Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University +# Copyright (c) 2023-2024 Chevron +# Copyright (c) 2019- GEOS/GEOSX Contributors +# All rights reserved +# +# See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. +# ------------------------------------------------------------------------------------------------------------ +import argparse +import copy + +from mpi4py import MPI + +import numpy as np +from geos.pygeos_tools.input import XML +from geos.pygeos_tools.solvers import GravityLinearOperatorSolver +from pygeosx import run, COMPLETED + +__doc__ = """ +“This is an example of how to set up and run your GEOS simulation using the GravitySolver. +""" + + +def parse_args(): + """Get arguments + + Returns: + argument '--xml': Input xml file for GEOSX + """ + parser = argparse.ArgumentParser( description="Gravity simulation example" ) + parser.add_argument( '--xml', type=str, required=True, help="Input xml file for GEOS" ) + parser.add_argument( '--m_true', type=str, required=True, help="True model (.npy)" ) + + args, _ = parser.parse_known_args() + return args + + +def main(): + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + + args = parse_args() + xmlfile = args.xml + xml = XML( xmlfile ) + + m_true = np.load( args.m_true ) + + solver = GravityLinearOperatorSolver( rank=rank, xml=xml, nm=m_true.size ) + + print( f"Density min={np.min(m_true)}, max={np.max(m_true)}", flush=True ) + + gz = solver.getData( m_true ) + print( f"gz min={np.min(gz)}, max={np.max(gz)}" ) + + solver.finalize() + + +if __name__ == "__main__": + main() From 211cd77d5fb2329328b02e5ef18a3229e5fc3a14 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Thu, 3 Jul 2025 11:44:09 -0500 Subject: [PATCH 3/7] Clean up. --- .../adjoint-test/gravityAdjointTest.py | 34 ++++-------- .../gradient-test/gravityGradientTest.py | 1 - .../pygeos_tools/solvers/GravitySolver.py | 18 +------ .../pygeos_tools/solvers/InversionUtils.py | 54 +++++++++++++++++++ 4 files changed, 66 insertions(+), 41 deletions(-) diff --git a/pygeos-tools/examples/solvers/gravity/adjoint-test/gravityAdjointTest.py b/pygeos-tools/examples/solvers/gravity/adjoint-test/gravityAdjointTest.py index 302b12a3..ed8314e8 100644 --- a/pygeos-tools/examples/solvers/gravity/adjoint-test/gravityAdjointTest.py +++ b/pygeos-tools/examples/solvers/gravity/adjoint-test/gravityAdjointTest.py @@ -6,7 +6,6 @@ # Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University # Copyright (c) 2023-2024 Chevron # Copyright (c) 2019- GEOS/GEOSX Contributors -# Copyright (c) 2019- INRIA project-team Makutu # All rights reserved # # See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. @@ -19,6 +18,7 @@ from geos.pygeos_tools.input import XML from geos.pygeos_tools.solvers import GravityLinearOpSolver +from geos.pygeos_tools.solvers.InversionUtils import adjointTest __doc__ = """ This is an example of how to set up and run an adjoint test for GEOS simulation using the GravitySolver. @@ -43,7 +43,8 @@ def parse_args(): def main(): - rank = MPI.COMM_WORLD.Get_rank() + comm = MPI.COMM_WORLD + rank = comm.Get_rank() verbose = ( rank == 0 ) args = parse_args() @@ -55,6 +56,8 @@ def main(): # Load XML configuration xml = XML( xml_file ) + # Generate random fields + np.random.seed( 42 ) # Set random seed for reproducibility if rank == 0: x = np.random.rand( nm ) y = np.random.rand( nd ) @@ -62,35 +65,18 @@ def main(): x = None y = None - x = MPI.COMM_WORLD.bcast( x, root=0 ) - y = MPI.COMM_WORLD.bcast( y, root=0 ) + x = comm.bcast( x, root=0 ) + y = comm.bcast( y, root=0 ) # Initialize solver as a linear operator solver = GravityLinearOpSolver( rank=rank, xml=xml, nm=args.n_model, nd=args.n_data, scaleData=scale ) - Ax = solver._matvec( x ) - ATy = solver._rmatvec( y ) - - IP1 = np.dot( Ax.T, y ) - IP2 = np.dot( x.T, ATy ) - - nAx = np.linalg.norm( Ax, 2 ) - nx = np.linalg.norm( x, 2 ) - nATy = np.linalg.norm( ATy, 2 ) - ny = np.linalg.norm( y, 2 ) - - e = abs( IP1 - IP2 ) / np.max( [ nAx * ny, nATy * nx ] ) - - if verbose: - print( "\n=== Adjoint Test Summary ===" ) - print( f"IP1 = {IP1}" ) - print( f"IP2 = {IP2}" ) - print( f"Passed: {e < 1e-13}" ) - print( f"Error: {e}" ) + # Run the test + passed, error = adjointTest( solver._matvec, solver._rmatvec, x, y, flag_verbose=verbose ) solver.finalize() - if e >= 1e-13: + if not passed: sys.exit( 1 ) diff --git a/pygeos-tools/examples/solvers/gravity/gradient-test/gravityGradientTest.py b/pygeos-tools/examples/solvers/gravity/gradient-test/gravityGradientTest.py index f14a1029..ca68499b 100644 --- a/pygeos-tools/examples/solvers/gravity/gradient-test/gravityGradientTest.py +++ b/pygeos-tools/examples/solvers/gravity/gradient-test/gravityGradientTest.py @@ -6,7 +6,6 @@ # Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University # Copyright (c) 2023-2024 Chevron # Copyright (c) 2019- GEOS/GEOSX Contributors -# Copyright (c) 2019- INRIA project-team Makutu # All rights reserved # # See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. diff --git a/pygeos-tools/src/geos/pygeos_tools/solvers/GravitySolver.py b/pygeos-tools/src/geos/pygeos_tools/solvers/GravitySolver.py index 1fb6bfdc..e86ed6c6 100644 --- a/pygeos-tools/src/geos/pygeos_tools/solvers/GravitySolver.py +++ b/pygeos-tools/src/geos/pygeos_tools/solvers/GravitySolver.py @@ -84,7 +84,7 @@ def getDensityModel( self: Self, filterGhost: bool = False, **kwargs ) -> npt.ND else: return density else: - print( "getDensityModel: No velocity was found.", flush=True ) + print( "getDensityModel: No density was found.", flush=True ) return None def updateDensityModel( self: Self, density: npt.NDArray ) -> None: @@ -136,8 +136,6 @@ def adjoint( self: Self, nm: int, residue: npt.NDArray ) -> npt.NDArray: # Retrieve adjoint: # * If any, gather all subdomains. # * Make sure to remove ghosts. - adjoint = np.zeros( nm ) - try: localToGlobal_key = get_matching_wrapper_path( self.geosx, [ self.discretization, 'elementSubRegions', 'localToGlobalMap' ] ) @@ -151,19 +149,7 @@ def adjoint( self: Self, nm: int, residue: npt.NDArray ) -> npt.NDArray: adjoint_field = copy.deepcopy( allgather_wrapper( self.geosx, adjoint_key, ghost_key=ghost_key ) ) localToGlobal_noGhost = allgather_wrapper( self.geosx, localToGlobal_key, ghost_key=ghost_key ).astype( int ) + adjoint = np.zeros( nm ) adjoint[ localToGlobal_noGhost ] = adjoint_field - #print(f'size adjoint {adjoint.shape}', flush=True) - #localToGlobal_noGhost = self.getLocalToGlobalMapFor1RegionWith1CellBlock(filterGhost=True) - #print(f'size localToGlobal_noGhost {localToGlobal_noGhost.shape} {np.min(localToGlobal_noGhost)} {np.max(localToGlobal_noGhost)}') - - # adjoint_field = self.getAdjoint() - #adjoint_field = self.getSolverFieldWithPrefix("adjoint") - #adjoint_key = get_matching_wrapper_path(self.geosx, [self.discretization, 'elementSubRegions', 'adjoint']) - #adjoint_field = allgather_wrapper(self.geosx, adjoint_key) - - #print(f'size adjoint_field {adjoint_field.shape} {np.min(adjoint_field)} {np.max(adjoint_field)}', flush=True) - - #adjoint[localToGlobal_noGhost] = adjoint_field - return adjoint diff --git a/pygeos-tools/src/geos/pygeos_tools/solvers/InversionUtils.py b/pygeos-tools/src/geos/pygeos_tools/solvers/InversionUtils.py index 6a0c795e..fff84642 100644 --- a/pygeos-tools/src/geos/pygeos_tools/solvers/InversionUtils.py +++ b/pygeos-tools/src/geos/pygeos_tools/solvers/InversionUtils.py @@ -253,3 +253,57 @@ def plotGradientTestFromFile( filename: str, h = history[ column_h, : ] e1 = history[ column_error, : ] return plotGradientTest( h, e1, **kwargs ) + + +def adjointTest( A: callable, + AT: callable, + x: np.ndarray, + y: np.ndarray, + tol: float = 1e-13, + flag_verbose: bool = True ) -> Tuple[ bool, float ]: + """ + Perform an adjoint test to verify that ≈ . + + Parameters + ---------- + A : callable + Function that computes the forward operation A(x). + AT : callable + Function that computes the adjoint operation A^T(y). + x : np.ndarray + Input vector for the forward operation. + y : np.ndarray + Input vector for the adjoint operation. + tol : float, optional + Tolerance for the error to consider the test passed. + flag_verbose : bool, optional + If True, prints the test summary. + + Returns + ------- + passed : bool + True if the adjoint test passed, False otherwise. + error : float + Relative error between the inner products. + """ + Ax = A( x ) + ATy = AT( y ) + + IP1 = np.dot( Ax.T, y ) + IP2 = np.dot( x.T, ATy ) + + nAx = np.linalg.norm( Ax, 2 ) + nx = np.linalg.norm( x, 2 ) + nATy = np.linalg.norm( ATy, 2 ) + ny = np.linalg.norm( y, 2 ) + + error = abs( IP1 - IP2 ) / max( nAx * ny, nATy * nx ) + passed = error < tol + if flag_verbose: + print( "\n=== Adjoint Test Summary ===" ) + print( f" = {IP1}" ) + print( f" = {IP2}" ) + print( f"Passed = {passed}" ) + print( f"Error = {error}" ) + + return passed, error From f55961c8119ace54705cdc2a068eb4d002dacaa3 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Thu, 3 Jul 2025 18:47:00 -0500 Subject: [PATCH 4/7] Added modeling data. --- .../gravity/data/rectangularPrism_C3D4/dm.npy | Bin 0 -> 83072 bytes .../gravity/data/rectangularPrism_C3D4/m0.npy | Bin 0 -> 83072 bytes .../data/rectangularPrism_C3D4/m_true.npy | Bin 0 -> 83072 bytes .../rectangularPrism_C3D4.xml | 2625 +++++++++++++++++ .../gravity/data/rectangularPrism_C3D8/dm.npy | Bin 0 -> 13952 bytes .../gravity/data/rectangularPrism_C3D8/m0.npy | Bin 0 -> 13952 bytes .../data/rectangularPrism_C3D8/m_true.npy | Bin 0 -> 13952 bytes .../rectangularPrism_C3D8.xml | 2625 +++++++++++++++++ .../solvers/gravity/modeling/README.md | 25 + .../modeling/gravityModelingLinearOp.py | 4 +- 10 files changed, 5277 insertions(+), 2 deletions(-) create mode 100644 pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D4/dm.npy create mode 100644 pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D4/m0.npy create mode 100644 pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D4/m_true.npy create mode 100644 pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D4/rectangularPrism_C3D4.xml create mode 100644 pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D8/dm.npy create mode 100644 pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D8/m0.npy create mode 100644 pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D8/m_true.npy create mode 100644 pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D8/rectangularPrism_C3D8.xml create mode 100644 pygeos-tools/examples/solvers/gravity/modeling/README.md diff --git a/pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D4/dm.npy b/pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D4/dm.npy new file mode 100644 index 0000000000000000000000000000000000000000..4e1695ee93e05503dd183bbdc95443fea20c312b GIT binary patch literal 83072 zcmbT7`Cp9h_s5I0)1I_W`(9JiGOcrFA1x$pDq4i3NJ3IXsSt%E3h_>gkcy)16p!<}zY|>QTr_M&L*D}yQNc?`uaQor&s?fG|3=hV)D5Y%` z5WozqG7V2q12@XVS5luIHlK8RNw;9YmgQyC<{TdKC5-tl{vZfWvzxp9b;!6OO7>4j zJ_EK!M@s12x4}{GOCs}R2FSS1jY$hyVd3z_*Mjb55Neyty)ecAPKC8@{9`V9oVsx2 z{bzvN#aBVembj4__2e!kTcEA9@7SEnMXm7t-%c*F19L~EZw_PDC}_QJqw6g?knZIQ zl77t5qJ1?eT?IU7;`L8N!K!u~Bb5P9yBz_!)0(AIRMwvw%YU+YsZh92i( zSSrHqQvCkw$V2 zbT)#3gOT%12YAEKW-SV{1=6do!9yksG<>x)u5<|pWofjS;B!<^j#wAYou*)%yz=h7 z_vnDmLD5>yil7|UE!5jXf&R%0Uq{vG7UAn)1xe^*YyJcH>mDo3dFF)k3wgtV{KYLFcl!4Iy&aF5;BoA(1 zQ-cd{Ek@SJNxpj<>`>Y;Yqseh9X33+ur2MP;vQl5%K^WLbpq-p))$RI{>Za8H#X4G zLc~Kjx=$apcm4iHdPPTfyC?Y+WhRiA8jJ^_Jj5@GboUiYLE0;!>1GrQJ!YGpQavd! zeep-N3y%pQqHkA@ez$<~h08f{V(Lf^=}6n!V~Q~@GE;W0#u(bhOY0A?2Ej^7$G}D& z9P1;mDtSo90n?yMyY{h=clwQ|i5wjs9ABT8+r&iohKj~VYFy~29D4uU&jb}RH(0Mq zvje`*^>XiO1mM%d3nSNS?T{~gAzVgA9TMVCtLbd9MSR+Nq^?H_N{eTH6(q30CBUHb z$8k;Qo|TVI*u+E9iW2cQU0uAlHU3ttry~x<|LmFJmx2rx#-u^FJSxvUKB>BySl@kn z3i?)<0bl*UXSKZ=FwThUyD>_KR^G44dM_5fmC@c<{})i9;KYE=ZU*$*FMNG<*8~f0 z7ufny*$~ON@%MI*65N$LN@0~*qT7k+p$%jm+)eIHx$uRH&sF?1YPRukd#rr7wHpVg z;>2Hizmx&{a~_(04Gx@-&yLi4v&1cU$rucVirV`m3Pwj5!>;E*4<0 zm<0{HlXZT_+F+ooxJJgOsk*!)BOzI`O5mfj^I7&MCb0F^sV-l^#eKsr^31D%eEx?i z9{D4*-c8fC{9?A~mFnT2D$9b4zl-JirLD2}*y^7n)waOb^YYt@r+;ZboW%@ESX>;W zX+>5`kYTS`@tQa_28KuN)I1o?1_vo;h->CxUJg5`E&gX+!!hTw-Etggy0z_QH$xxo zrRtS~9F8Cq?Z2)~0EUatN0tuLFyWYFNc}G> zsC@5!>a3Iv=F6pjTvwxo!xg{tf?EXZCr zJ?xgg0;%uRi!C>jky`Tlj+v7kRF24Kw0m2_m!k~@M<$uzofElp@lOgyiVRH0`moTX zk@?wv3kUZXh=fSJVqoBnW9^B-mYA7NSC4vcgb9}m*2|AkaQu`)L2-pP=w(Cr7E_(OM@F$0Z zu{|Q3q01}~I`V$3}ze|6qb>u*T<+oK!iZ~#y9~F7P+8j!@4F1tr$AMgHqt|h| zhM>ABmjnu>wsiu9q|AV+* zvRKlxOBZ^!US~F&YC|!}+sJd22a}VlIzzaQ$lQ0R?95g+N?y+FJ{M*RYLj{!*M0yr zW}c*JQydY?9M*>h>H%N-NbKI9KWWi7O0xuy*#obRc6mn8439H}&&+4hu{_@Yx`V$7 zUQ%D^G>PDX%j$)k^=<}O$Z_;PmmL_{SI?%vr9| zJdI3@Ss%OXRX+giRK4%J#Sj7oLlQ=p{i=($h~Mp=N5;1GN0U_a)Pd~t?dz5Rdk{Uj z@L1{D0bqR8<=GQ=~ zPfgYHy=Kre_FzxYBpH-9dp^54z{Zm^(m&(IF6)5Vr-^FJHxjWEBLHoZ6B7T#LSOT0%KhFblleiCv3*H$&1lD~8m=@VS} zFv1dR&VllWWFF?NNt0Hd(ZwDoy3)lTW{~hROyP^78b(D z+3RdEbolAD4-__BZTczSq00mR4#lDKLwY!VpLx+?TmYJ1HAShH0Q$W-Yth9B`R&Gg|S=D3$kz88Vw@W%U9<<=;tFVi=()g1G7#AJxRs_Lgw+i$L1c`fXm7cS2BKE zKxL3t_q!T3v@iczH29JR1I(hyg$E3j6@C(?+@S^iU({RpIeXsP_g>gV{M|!685Dhvn zkyloWGl4gHuKeK%X^cF7SVm)q1>Wy`YymrX5Sp+)KW`5WG>j}IInOu{-1KVxr=lZ< zv$pQb>WT(C3 z78umpYu!U*x^Bc)A{mMVlX~-Xi-@k8xv(WIH>KQW4P&=_HaEW`giA@S=#Lg7bJRH zZk_|6Zmqbtua$*HuT3wdO;e!L$N%wmaXf{i5J zo}?;%CEVgiTb(pQd_E<@db?$8V7ADb-;?Mg-9p>msfO^dB`NoMy^bmR&M20sM_FP* z%zLe`QB=6PnX2{}#9zGIX`F?@qQdJb)eLO^Zj?nSpu<>nhKv@ki-Ty`JPKjj9ucAfNLKwR);R)y52r+}C)@^2!p>5>CEn+h7ZAcQ^0J-Dw4~5rtMe-R)4k z;@9rw@9ltPlrpI>OZ3Up?klC%8ez%VjK11KR0w@%{-h|^3}PbYE(OQ15yM)x3_Mwe zG5@vn8{V}B?9zFW>Z6BLlOvABXL)G7`)u64LrlDz5zr$!ZH^|o#Zp&}vmwLuOiH^e z531Wz#Sa1#DrNNRF8eWIssFC(UO@(yzrDgqJEQ}=E_T6f_vJ{sH#8aa%>j8azOp-i zabS;y5AZX1ActCiC`$io~=AIpjUZ(?nmtN;@+CoKH(alva zZ2>L*+U7pJq=SWm@u~AfAD*gFQ&#hGIaXWWkH19pl}@)-$o7OXpr!luzt*n?IIQ^m zxOadp4!w9Cb~)1qv4hIzfXbEC7rKag;kiQP*0eblKczJGB|1CeW~~77+cpdQQrC2kzKh^C z<%d@J@hyj*nl~?=%Nl?MdH1+svOS8X(ogY39B}@t^0o&rI3SFMlIJ!nVP$z@3+Jm2 zi05k@J$ix*ddaD`eu%T-dtd$y9V&qCtl{;;7EH`~WU@p28V|_SFC}|FSz|+|QF5RN z8`_3Lg-SGafO@j3^WqC5$XF&PWp`2z<9_b@&+dsW#Ob_OzzO192_*Dm0S8Fv+H3J<8v{;8S|(PW zU_hH%Q!j0eF_Q0g*=nv4!&LX;KBH@7`1i zms&%lN}%u3m-aAm>gni68v`ZFYYs~d8e+y8gTu$pn?MV z9_6v~BA!-&R7WRrT9fgj)d5=ZQs9w>|Kk zSsr~Uiie@$V%~Rm7y-Fd(6;2S0X7})8QXN!4)5GJkh)lLWTSnFy8`>Eq~MTBh4# zLI*5nNMtSJpm5d3_s?FaLHC-GKSv0C-CuC=-K}*rC~e+pLr()R@qhhqK#>Lb^=Wq1 zV`ESrolEu9V>C4-%S*R!K z3}==yp)Ks`Er+*c*sSPx=s+wF((IvfXqOepS3Hk5JPH^(w~&>y(h_h>;2Dv3f|&b3 z`dIn{M=-LSjwF*wFz*`H+?PP`#GN(U}NE0=-&KU99ejd5Cr#L@AGNchUoOWj%4EU*O=rD}30)P+cvGLAr3a5_X6qs??XV-Rtvz$UJ^fm*&LdTi`o@NZYgMH*{Q9zDi|gd2ia?{qOh z*fuG)q0ScbHdObUfDlS`^A;W;2YVi;e>3@~3G;&LcHuf4D0$W&*yW=SB*UHiPt+|3 zeDW@I%K!~Bp6ZSEwsY|uO~`-rAQRzG7?(U@1bWZOC&Fa)AzbCp&zu%5Ac@=D4o#7Q zSt`A6*BLuV-@CJv`GbjWaTi18))?aOqs#H~S_J>z@X@oo@pqlkfNO;j$q0j~Zck+r z0QU6+Df=heVRp8yrGb(<4n6Cty12*+9!@^HRzs{!w>n8T|d{#f1Kx zn4bf+@*p;)|MGvjjdv@5`j3(gpqt(cB^QDNhwi$ik!V5anEk@Pp7BAy-Ib)wmm-)g zd+As69!W?zpDorN=v=>y@OJr0{4G-_`Y67_%mCl)^2?Pko_!a zCgzO=x)+pPEAL~1xiopXYd;$j_H8=9b>RyQ73Br);tEU+TT(F6Xo8VzPfG}OE`n@V zeMi|`z>q0l_R&Wqr26(ONqHrQd`WXDL3hj`b=MuCI01rdo!_)oDaHoEMPiSHtFd8R z;m8Ih1$}r$+4SjnlO^^qjaN9MBZeNQZ=M?|qC%?tx6SfafYP~gdJmld97m)YPVaQU z;1E7LNfjpK9XZ?)a(=Q7tNzv3@1Cb!kN>4sjEb{ z+XfRx+hTUlh?1*4kR>oVxLam_HTmE!#K|rrMmU2c9?ZZ}H zDY3W0{+5t8zgq?H>BGXfIf8ph>vYwoEaKq&*~Kqp-V^)EK4(^l=vVZ@mIWtW*Tupe z@4s&xH-hgGj}L4U;h^NlBa;=)M7=&<6yx<)7Yh4?QigqKI5?+G-9qTwmWs+>cMVvW zHdtzUJlp{KKc(s)?ofu7&zmxrWzj(SYXte4yA8bg5+0fr0#Lih?$>k;1HKGhd{)=U z0+JG}FElpBo-ygqv`THTX|ZK=m2%KRVDF>)TV{|^Hr$i*jtX64$7Rp{$HeJa9g%=( z3#?f8DtB$H6;3)moZdm)r;<-%W-m zVRN)6LsA|t4qS|>b6P#sWe>TECS^Malc_iB)SE5J%zl zMd@hZs9*e3z6mCJX|9*)oaSIt&Om_hXBtrNwF&HhOPv3|H=GOFZ;FA{1H1AaY_P$z zHNVE4&=F#p77_nu>q1m_p4UEY2Nh-IQ!(j)9^!}G!rqa=M`lpWr%V?j{W-p0jk$1h zY2BpeCObU(So+w&O*#&pc9d;>q>LE6Uq_MfF&>H@d3k?jo=u}WbTQTMn_=P|1&r=YwblC(PID!hx|X=V(MI-P z*`*|OS2aISPN(1?dGG3a5gwKees`J*w84qsGT)iC1gBh{IkfJw0Sezc2v=UJVY__v zH`jeczk3GWXZ9+g1zUn&>y!o#_H^!@bkV|Ek%D?sIs=BwDj$yrvmrCLf@6N#9^@wf zc{2zfamq6?Bs_=$sa(-@XC5%&LiWY)_g-+YvTa^hzQF(|#?{}gHZ_Nh{p-R@tvRS_ zT_QgHVgpJ+!57_DJ(9g@YKm>Pvij5p8ECuex_RR%W#oIXdWC`fc%7=itE_J~8F=-% zW}u%Pv3~-?1$}POQ6Vz0_kNEIg3%+DHu7RL`=KA;^4u1Mi&P>nPUs+3)E{}BM}?*~ z=X3cl3EorqY)L=CCqqg)Q^MzT(B*^BPZe`bq|TkJNES0fpVNY*3Q1`UNr`*aciIsW zL@E=r%V%gWCa(*G|c92jC%f^;K)sl_XbC(FuVRLm5~J?=_q|HyvPQOJDQh{ z4Ny=Zj87r|I~|I(HZh%-&;UoiSC;PL!o;iar{2kC;P#{`KBLA2-rJYA9Fudz&@gYM z`BELIetyfjhTvA>^D1wn+j)?vRP(*Y6woC=IBW2qF}gVn4@LU1Q8>R*B<+nZq}|%L zXrG@2G_Z@+7$IgD*S>Q8wjL9fs8SRrPB5^J@@2{S4njXI{j$PIObmj;qt-VYTcLla zLA$Rb8Dyjr>DMKB7=9*iFjk3y`H9y843zazS?n0OlgEbsf=0oF?*a(N>dxDzlEF=F znV82?Gu(ZsVDYnVJLI+AJXW*D7(}Bd_{#_`H@bQ&H-bv!9kjw)w5RPbckrD)|1wJ~ zIa1-~5T^)Qo~hb>{|LT1p|nsM$pP8B|8BL*&~d~5M|ry}DVV+beO!LKBa(%yYxlP- zCw$(U65@oe+F0p}z~HwLGyAP9Xz9PC_w=p}Y=}7F>-&SqpA5J^ z@};t{q)751zabY2k5At1HI~New0p)+$o42^>NDm2ii?7-nRT9*9B^^++x5G?GH~-H zVZBR&Or%b&6jQHJhVrbmfsYsTK|cA&vR!^Gg5W2loF=$$jDOMBp_k^MocQ|y%heK* zn>pPo)K8mEa_j!u!Ud0(x#@gYDp0*vC>!l1c%|Zgn;UPI;gb^`PyKQUj_;^*ST)%i z_fORFweXo^(%6f>uG`koqIc*mUpnFIuxo3E@(n>CqD7=>{rKJ8j^SQWLp(-(b!|- z6V1JL2wN3>pPl(x*QA|yBJ_bGrroz%i0Gk%a=8z5vdm$bdvE0VbS~rx2;^Paxdb2A z)HaqS5qkUTT*KuEYqN6cj@|_jj_i!>*N@HO!646Sg69W8Fh0w}ML@|LB+ zoqNTyiWhmvOPQ8atCq#oZ-wVqts->mZ=P?ihZb5bsF;3#Pluh`-fR;kbMbbHUwiFp z9%i3e{HAaX54D5;?Vafkthlo+OrFgB_gf^AI{| zVuOWZ#rc6v#vtIU(w#%RJ}EUNBH4Z!4j-^R`d*F;_anJ~4bH&*cvvtP~O6mljHJ07M)?at;G?a6 z8250xj%qUWa;}8%$M`QDU%o~aAD-TFH2NG-w{BJBjUUp*p<_ZTPCQjc^7cKAe=;3W z_{`w2#1_JjE{jSdy)#5=^&|S<6FTs)gC|sQk^%YO4r&b%{bBl!Ki7A^qhaRZJ&{v9 zO*BdT6tI!#!%|zE8mrpmusCI9&ZBuEZ{&ro4<^U!1R^9%V`a5bcI!qFu^n1O9xV9s zIU+}u`B1lg@UkH!i1J)-DHA%gut>^`!^9pNkp=!74w77kENzGE0Ut|DY7qPExA|<+ z!EP(u-C4Oupqz%a<%jhGi>z_ftLpVcn*~hW_O*I8FNJ)6n|DrKTLK;m(&}|?x;UsU zf4DWw8hka~s`!lT(c*!KWTrpy|G$oX7A{<%^}gU9QMhe^BvPuE=)g2h_WGN|%Y}B( zbMHRs^=duLDzd!i{MZ6J)dOU!8`wCym-X`b3kN7%x9wj0GaCq8yISSTei3|QHAOD} z%LN|)g9#G`6&}+xk^^RVkp5F@L0{VrU!8id_BDi!$zS>x&nVI%?8!y{6-73<(u!)O zRKi22ilyJvJs99$;dA)OZDY(U6H`9eEsmEn7kaNm*kM!hG0Xe|)+ql-XeZU43vS1w zr3|yRz|ChOZp$t`;Jd8i6&_3YHOZfEhJ7bE*=a-n0HzFdPj7$zeGB1FRW(w@+_WKc zxykL)Id!!7BdvV!4i8!aM5l`jG(mP}PPlUr8~C=6c1h@bq-BoHC)}y#+cKP@ZJ4hHTu<|^z0!Q1p(sj+6@a5a| zRL66oJ{pINtq_-mJbJg!_glp4^c&>vnxH_dV%CNp2O__5xBbuXMh=8j`^-w)asXm! zzdr47zy#;2h_RodNIIUSncGR!4fSg`Yt`9ECTZ*)BDiYm(L4&joe}nUd4DI5F~H$z zrduP)0mE~%?Ef2QV%`bP4p$-J^KKqKHA47B|yHb!f1TsW^y+(zinK&iq$X?U+4@yclr9o)pzwFPxW*zSy&8J^YFtE%A3C+}vPYJ%K*ZKBt@m)<^ z{xhD;_|F0@uH~j^WF?4%(>D!AV zT97+fGJ5e17iDDKk`84vkyZWIifcuJ;`6m0I>*U4WnlMCxYHV%M@)9buA-oxk&ll0 z853AsFYqg{iiLcZ`^r8hX`|QTy1m27_K1BZB{T8kbvLHJExx##=zn92bY%}Rfs}KH z{$xxAY?mEq@t!93QRqjtLz|aEVaKwi?%qU>XGhmMTXAg+v|d%Rj873Vg4uNx!pFAn znps_OTL;Bs`BDT-c-a3{G5tQlZxkM;S~d+E6FJ2-)T8^vue|K2T0Q744$otdN;bnL%%=XHd-6f~s$8pFMW z-^*~r)XmcMI)c(5I$QKsS&@zYPX)`J z3mh@Y;p9=-iyYuJH7h%^r3igTBW;=k6dxHbELY*6TfhGPVhJ$}AJAWxmal{1;++0; zK|8Q0PhXUks|~*oQUoMB9kAG4+xVZF9)?R=RHaQ?;(;`y_RbtSJbsk+$$Oa{M9XK2 ztUU;LM{}rtKGza0UK?bj#sL7c=%D*;Hlbfk)t+>*v9137ePKt!N7=hP&E3!n`oCYf zma^IeMkOw#+#>oCo~K4D=QSI&`0ZV)mzkr7u;Ar4U$vp()uLlp!p*Sq=zw(Qml7ifO@YMKlgIb!i@OW33>z%ce8kPX2)S$II+^o)1-%m z3B|XTXB68)yP-$!$twVnl{w3p4@IG0Z2zigqJM3&II+Dq-wKzS?(jK79&E~DbWQX7b;jxP?tKo8)BZu zk5xUt2;uy4AGPhpLY?rFGwL6W?Vu&1`}buT4knzE^!Y`xhH^pQWJv`w{MyIw8|Mi) z{-3z#CZb-G9_V2bxgGMA+skCv zF9)B(ZSA>ZvOsbOI}}Y<1CNfin>M{7`U}^Q%e&pF_;cg#D?L~2F#5Pby+<+;zg0N8 zF!Pv;YTu*|H8--cIOnv7OCn%vKxgmm5*AE_&$%5GHiF)@$$wSaDLBm^vaaN-1I8b; z_cjZ%!b@RuORmNc{ygJMQ9Hrmo4!B$(mKJyMMJBP))M*WnEX_!vjTLO4anq6Tt)}s z1HT>#JF8)}`qRZvDQrNM+99z>6X?{=dgwgLB-Z75*DNAuCG<(`QrkWUd~v@>qCc^o>+bIJ_uwaZ^pTett4WZc)Y#}G_mh^?a~W5cvM^WX_KuHTWK5>( zizxK~ZaDhuNalZBnBV@_R7-&iPLGmyJKC6Ge|);3GruWV@Na$iAX^74f+Z-=g!SQu z=mkf0whfN&=M)Slkg!>L@$S4@qTZ&cEc(8A5rzxRJuIV;p`?E4;GL&3aN~x+_BuJj zpZdM~L?^*HeVkMuc%|tB-@91DqwgB)5`Nt?6h{&aRWK>x-9)ZnB4saX zF!9ty;V$=S;`MY*Pg_3F#?T8(Q#y80LBah+T%MabTC8j3buBc1|`psG1{hAUv-oY zoliVR> zdk$QFbNTPWHft>Jknz_4rVegKBKw+ZZLw|T_`XG_HQ}}JRoB1YY>?_~xhP8zpk$=1 zY9qLzri29~+z4D=wya z+rosT>*$ydaZXRKx%`}_3Pu?>G{ZL&K1Khmsjr6-Mq3xGzwp@sNH>(9KR#%V;dOsC zS4-RA_p~FQtTK-248oGl@kr~d%=Nv5pS**1TI{qExJ%o_k>3O4|DIbpP{=~rb2hb%AI3N|=P5k4 zT?371Zc0bW%plLQHS{Qx4*oHGheZrXP$=HxpioJMqdzU5bi84qt&RPm_0@pA|1Fkn zu3%!|XF+=s!58%wryTsc+6dh~3zfbywE#y-JKeC53vy?x5~GRSM0nu#Ztus0Use&b zg>i%p+YOHs8p9I9!^ppe?KMHTdN0%z$zpqK<5#CVLWkT`h`Svt2`4@^A6>h_5i{;F z3LK^Mp|EMqzod0$&~{=?(S)EGM!FO|*M4FLk|K8e_7gH#yrKIPY;=I0M*FzD=YSb6 zs`K7&G9l&==AF)^)A8`S#Pm;kWYlt3+bR^Tk39O`AFLb|RD5;5sICZLET%|4ON8LY zeP>G}cr@UozJ9x)PxxIWJ33-6+W~K~!q$C`1toedw1^}ooW3*i`e!JSV}V%jxVs`? zptI{i5{r%baZ&~9s&vq|n1G76#GHk+=-mgC95C6Ij#tX5cMsY2{j?pSvy? zV3VT_l@r#tlMblEZ1vj>mc)9_Q$O&LqG}7%m#c&y68lp8%}S9dZ*v$PIeI0*_;O2+7+z#8V%ZfzC_ILm*NAsQT`qF+b0NSFh! zq(4xz)E2JZ?0O++Z;PoqMeI&dYqjF+2QvOZ6Lw z99~4yh6v5?wAsYxyT#jx=M?aS#5gDb)lTjm4=kD3w?ia#yxta@M7qnx9 zYf~bhdFv5*IY07!Nbo(IK2k?sBX6@noa?ez<>^a_KHuY>&$<{@!XMF$eDHjm96X%+ zE_Ln@Q7?+8CXO@*m9%g3gDw25Z=BjJaO^a~>b@7Rw-LOPqaS#aMC_-~)bKdcW*u}Z9LfwL zINDIY&xsSL1;`9m{Joe97u+Wwp)M2jEM`5ExF#rBH2W)NK^H@t{2Nq?2(Fk&5hq^+ zOyDS|$Z>mV?`!zo4eTgHcts+RcY&x^uT#te$%KEf!zu8Tm>t+>{8O=*=V0+Y-E`kJ zg6l6F3?Rj+px4s?Jsm$Y%;Gr+A2?0Sn`HgeV}CM$+O&h64#YW^y)S9;d=3NBuPHaH zMibmF_qV0XO;e<9txGw{P(Y8y(1!180QkDk8QUqB()gbCn@rqzMHBvBao1E?3QO~S zUz2`O5a)fROQ|XlynU;kXB4rY&U_1)9@c@E#|#aksdi}lZGLHL9SfpeVv-(DsiSCj z-21R99^n_c$J|-G82YQ1HB-9eFt_1cDwp5^hejF%wW4j&!ZV;_E`Z>7W#2~+)(|c$L;O^hsbFNmW2W6DRl+~*yj$Y2l*k?KcgZqR z#;MXkJHb{ml7s8^9O>hrOTVl53AzQ4j$AC;WN!;i^MPkuI&8twq*K@V2?I`V>N!ML z=OEdr*1IrG2R^ozcsw9-aI$r3BDT>w80wem>UYZ!eptRF8|7G|o_$;4>c4u>e@$ny zd&VB57f8#Ee(->FjVJ3{Nyn09RCfC^W6%n7ru#+HaJ^yyDXrH6XNPLtkG!;j&WR&c zUN^}w<$2(?QVtua2XAj71yDfMg5CHshlOpI(qHrj>0{PZ%G6s;3TiCw%6Mr>?CHO$Qao1 z=T5{gIS|-9$6Gb3ji`VBeZJpFoqot4wQS;C@0ToneElpNE>!+^%%hj!d6aucN9MIL zPx!F*xyI=_fm&Fv5y?SCnT%KJ@^o0dD7fI-Zw_w1U6t4F!h_aQn=wHr4vfEkb;`xm z7+K?{#}rKI=xX}b*qPuEE|<2?%0AM@{)fBz-X5bMUz!r{8OI2%zMXMgkNU!#CNJn1K}&(FO2t{Hn^Bi^c68@A|Ia} zbn_S!+_+(0j%y6z>Zi(3kx&biRQl6ZTW*EVFQk2zKgLCuF$P6xqcOhenhB_ zMKXcEA&T&ik54DtJ!23$C8pJ)ZpvaK``u%cNo=+s! zy?4m1HEShtw)Y2a*{DgZ+o~+rP-D1lIlUp$mH50HMnm%MI>P2h2DRQw0GW5b9$u-X z3o$9yRz)|gG1d0k=ZnaN+3!lF?V&`@Rdni^O@tl%UfZ~x8eoCHmQ{J47~o<}#vXFa zIVvW%D_^?v>z#PmuyB!c(36Svq-(J~M1FE;Q+apQVi{b={Bua6n2RN) zp@yMv>>(pCy3BsF8D8~}bGopKi;;|BzStr<$UC3eQ9y8)V(plFjDvb;JL$+>@FCU_ z_sp)S6HJt(YyN8vq2Mp>CxbO(#9WB2)%~(h!bq*r|2}Kv2zRCW`wcp6(CfF*%N$}K zFgHYf8q{TB|CSy9A>o4;_|NDpts>?Ko*e(tHzor<0h|c-(6mt2jRy#RBZWB+f zWgk)3r4!2jC|kngqHQ;~|8>Mv%Ts%N!yF1H*GBdAOM=VUspyw=1g~0k z()7Nx4dGmGvvl5U29q>J4RfX)n%(PFWxwKJTZF00X_h?TGOvxxHkg2-g&bL3h{##5 zA(=$`Sp(Ivfuy$J2(vj?=TeBg;v%7M3*k?R9Nq7@r4Mqb(2^y6U*?4lPDupZAlF)f zTlP1}eq&RjA8ZNn+$)7{I}_*gJ`=w8!l&jh5wakz5;On7%npXS-c|oOAPd1>bdOyV zJm|f0e)?h#6K?)48fA=IBRL?$CDmIE`#;^IyzNm1@`9mpIadTRD}i}?=ifSuMt-`= z8Usi^H_$21HOEH7lFdqgh}_%24 zI$c5J)FVB!1h*KVtjw8&Iy-H6ktE=9%+C^K=ftjiZ|8tdQ*MZ%j}b~=`I>)IkPa~| z|4v3Sh&)Y93B}BKw>--@Q&0`@kTmTmLn{zi2;`9p~QIGcsSd|J+1f9Rlp=H-XW zREYCET}~qOlqgisFU#m9=0+Kb#SN@Y49x2pu}{9WK(o}8kET51Vb2xnr|Rcgz$bNo zsp*L(+D+s8OHLo3Q>kk@x%ecX`C`{uV(wpCfgLPhhtnz=)Mx!%%-g^AdO#uJhxsek zEqrOIOIUO4Tk-q^E%dF>KgrJwjJ8d3lvSdj&r`{mU2Dv+c$>%6Z~zxc_ebDPiX@zl zI$I}e$3h35ceVO<27V`B6>F}tLGmcnju7()foq~A%|0!|qzRMvvre{fVTn~pr2$}I zuSHkrzrQqr(&m~LFCwoS)U)$b3c#29XBHbr+avV^HF#E(0hQHmK1*g8$m_WGVT&-oF1s= z2SM`1Uf+h0ad~jgF1h;iA(t)^cpL%DJBnPC9MtzIdu)w^i_&wV~ zdEjx?xxXe~6Z(r@fB7d&g=+U1(=Z~pkn!zGN$w6i^tin)IGe+TrciP+^COY>t#4S7 zPxPJR#$VQ^za`%1hv^ziQ&l)Qe_EzCn~MUS#gg$7_7F-jEHLHc0N*(ec0r9Ih*s!I z$P=96vBk~nH7^KG&wq+}Th$sXtLtOMw#wnWNzl1_;zZ2$&c03OW?7(XLH$39&ODx~ zw~OLPNXb0U^E}Vv-FI*aMWG}~C8Plnk|{}&gL+C?bhU6d7;ZTs0MA(F00?OlkY>nCdzQ50EgUaYw^tMCO66RXVC*p@Vi z{&~7~(j{muve=FN5!UIWhz+(hx875OfFCDMz0J`9Z+-8i+N*#ahdcgDWRXV?5w4)h zI6hPpbQKo3V1-Np%AfiwOANhiUHR~=6_zIn|CC83KCRlZ@5io#a)kE%TY9}Uq3R~++j=dKgP))r2#mM z7QOv-8Nf6SDf?7MgLHaq^sBv!_*j#5Q~y&YJkMO^Q_D(prr%~8N1{|fDI!Af@Ezh; zf4fDv(POyg`Y(@y)2B`G@;2Qf;RNCrD8AZPeL@r*grwYHJIS#xPCj7vsKUojieI`? z$@3E^Qge&s>{n4$dCH}&u+3+9%cpSz@Cm%Lx91PJ{!ugjM~+aTma8`;eb|niZ?&3Q z>i~RjUHje7Nqh-u7TVfx{NN!uwQGF?FT+k*AxHVu zg2yq27#rC}IV26JrTVS1^9K{}aj00tFWO>D^xD(>0+t|$ZfWp z4~u8=s!`Fe*!B6UZYDIzy){>PO>+6L9&fua@N7-v6itwU^S>^6dVQxsYM#O!ZU;5= z8+Y#w_B6l^v@0ULoD7US>s@*^Q5roL9?jK6*}=r!^{vMSbWx0>U9Xd>3LY!N?Itc6 zW68h)vF(TLpz3npKB}KN`jlF~Yn(F&H`etrx`zO*vzqCJR(AMp@2jgVN@O0vlJs`X zJxwg71Sbwn>Jr^WXghZu;Vul63o)h1{8q^+ZOdjo2zCk$t5i0{iAUq%Ezz8)<(p=c zpk)Op?8Ua)zyz!hU5(_5G(vWnv))$Q)gbt;lp(c25WCB`K9r_PL&LQ9{Wtp!AX|At zB`rk-+zzbaaCDVNMTdPwx4jIJHrc+g#h3xDdR)udYwfXh=~j(-ydC7txep4Y$RLyi zNL>CS28s16gacgFL9hDxwr?AFA@FT>*Jni|G~f<;xpe~*j_PX0a3nA>(Yx&Q+S9f; zShYz`YCXW`RcCxy#L1ZM3Y+4&WEDtm`yhSkBbh6t4P2;5GzWM7S0e4C@1HiB7~WfF z4cdxqD=Y4j`<4;G)4AIaCN^zqX}rh^W;?&HT#T@Q0=CY1lTS2c5u!aEzg%8Z^_NQ} zc#k=F*0=}C`5U6#;-TN?onXQv4R)=e$C zdB_t#l=H1b$0PyR#OC|&fDC|(&p8jJFcDb)K^ z!^E~STFwKO#tA|z59_-ZwGoY1N2pYRDpVsIaYnb zhRoe5yb?_-Ffvr4Gn(*ver6oodH4{?SDsk&WVXc?jl^XR_XsNEtF`qf#Oe&O;YRM| zUC))E`^^1|4&UhDz3**9m>C`Y6ItIX)!Jc)Td`y)@s~{0yx|RV=0%q#`NmFTI~;6S z%wQ%c0#oqcZ=D^)Z?(h;>W8#|>uX)Oh^02f@8WuYYK1X2ICws0J3{>HZg}9y?Zuk( zq|D|_VUp{8@1weONDfnG@4UWD`X(0vtNI6$?2zUB<%C%t@l!9%|ND>3Q-Q^INVi59 zLB!(sc#jzRq^{51ugnJNRri1RsIa5i51#)9q6APN7|LV?0qRQKA5)K;!KmY&wT;*5 z;M4D-RZVpB;0wZ6e=m_dU(pYmBMRcdxcEd3Ydt)=ZNO?))eKqQdf#$>H(s+eI1{0x zpg?kE%R#nb)|fQBe9I)0a1*+d15f6OV$&OZwd;*7$)$WR+YwCW-6Kme+YE$&xQ6rWs$#oNdeN9g&&@pHpjK5%=Ytyzfos* zi?=n8j-y3(Gx{&&;r`Y0D`b0Z(YN_)-aytT2I}#`g8~2`0eC_R#ODd)%g;B9hUy#qp#anvOgRm^}SEqhEmg{nvUP=bf*}am}Vt&^n;QeOx@Ku*WGE3 z`F5YpEVUR&27Kxj4huCC;kmmP>Lnn?DC*q7I1}9a-y~DAn~KeSV{&68FV>u2Hd`=2 zat&8ox2(NHd>ZGQI;Dm!$=tVKlZ!Af&Nj~ORH>n3*`cSF{^{zN*K_hwyeX+`+bVDU zC#M5KwNJja9oB$Ettb2D%P25pH~xBeB&m-(n>x|j7Bhlt4^JPkg*@&r1y2+h5L_>0 zqTwKhrOp*X$ER(e_3z&Q8jp}XP_B((W*rZT^E5mN`(%SZT~^P;p&o>L(YZ`TZSkIZ z*%b{6sf(`)l$Es!!)W^I^`aD0*iV^FDJ6NK9DZm2Lz)UO7FjN|xt8QwhcB-#I$@9f zdy-l=^AYaJPqt4qE+(|?f8kq{E{#tHcwa{TA^w2UzL9csD_~Sej(^#wgGG#uS!^nH zc+p}^Z1N5RJy>;DCdgTXecAYrwXP%anN13 zYzP{s+l4tQ3^8)C(lC}weUW3 zEsX{~-4Y8UA>{lhc+K0q-VolrU-?)#$qrb!RR22|Mu9S|*P***tzlPjs-xrtL;Psul)}kGn6v(gCAl?t_su16ZQB|1KQThq9N(dmgIM0IL`8uA?iVar{Tu z>T3*a^PE3hAS@3o0{mvd!7VlQb``!>kErOhTid^5&J@jMn3iBJ6EF>fa}<6I2c2G zL@M0UrfZ0=ZedgRqnAWi+M+v|NH`Jk3Nv=O19VtwUb;WypafQ40j%M}`Y_SxBInC% z4Qich?~5d8KwH;>&_uoY zHvKgNJ4(ZR{~e?v*Q8Vz?@c_xQ^|LbhLUfGZ7V!1imQ~oGkn`e~tPUK^ zcqFl*$A;*r9<1|98jz>ESSCh`bAo?}|)S77|HY2VzDQb}Ms`hDVfjwT*j$MlTN&3p~oIi}00XZ5>ber%> z+h|V4lILW}b8)qds#0lquqhaR|@i|W% zZTRWm$cI{QOKP^vX<&U04{yb98pIc#H{NDCLV3+sVZQ$s6V)sO;6T0r~TJyjr#Gy9OrUKP|Z$>}NjMRZ=8dry8S zC+d0hZ27D#iA#?Jd}K?gV1FXo&REP4T7I#vxRGrRa&I=RkCGKame$r)V(hCiRlZAl zdP*PN{5~@`_mTd=r7PDthY9mJ3BK_!?IDqavby{M@kbRK=&vYeVu|w7;s1yZlue6# zqitszE4IsvQC)8`C(p$hqpC03ljd0mJA4`tgTTAR=}G-lE!z}QL!p5 zu4?9>21pq-KseU<_) ztuZB%2BK)<8OYCnnu>5GWhm>04bC^*<5!5#LvF{Md>1p4U%P8syx}d$H(%A>bcArR z=Kr#K8(!0d;DMptVN^0NSYN~T>O85}vgH1@wlE>G)o1=+xGFNG@)yqP)6r_1&lRBn zKyF>bz0>;`SoLLTu#xa;tfF52%z9-8pXFYLba)!0k6_5=)L0St_=R!d?j#j^M3Mq3 z3uz$d#5DXoY5?=p4c~hV^zflg(7AC54JM6eHmz2(1(*9uX-SJ1k(=2JmR;L>jf#w+xmXO6Lf zhp&%eT$2US-5XcT6F!HO4`X0bM**XkQwJs+<#4Gik#}gjF-*(5-Bw*^1A9CF^O_s7 zL(h(f!dHo&)@=7%BlI)L�JW6>ic3wpFVd)mjLDYQnhXMYuZ3`JU!=hUuERC)%n$ zn^e*0!dn3zXG1uWByc}z#^(QZfon+#FCNr$Kd@+Gf)I{PZ+W8*SVR?MCD!qN&mWHlFG5P@p2Ng+>|_+*Us z<;b)tYeHUt&C8&%|^R@aFzYFD17yA%zh>zsZ>9IR*~MLe67 zA>*6c?S>2nmPyAr4^688%fp?SMf?)*=(BZh8_AQ*Z#?|wQ;;d(Qu30*EpxER|6$g` zr;qqcd9N*fgz|v-p{-wz@Y;$)#d}7HUKl9v5@cx)W~-l;=5mS6CwXwcq;=P9bujzQf~tNS$z{wB z-kdLx2ai3<&-~3P(6FlZ!0%is#6HIrmtrNM++F?J2~t-_b_pF>Tqg};*P|3w-WkEN z4)ffNWF{6SUQ`ewJVUOP_na1(+Sp%vBJl~~AqK;I)YkPYkx%@nx~r%PmK@zIC9g>G zK=-GOXO%5r?3R0yJw53yBu)PXDTvUKhy3b>`K*(nN`LALSYbZNg0bPxIZU%Na(nGjAgteB9*HHGP& zD=a8@?N7oSH|aAL5-;!0K1TnJCB^mmMl!$nA!sAjPWY0h3;lv^fcEG2 z-M4Shfu0?Ys*HD$dY5`;*RMOPv0;R3%72dNE|rpN-Z$yPp#tGOo2RHavAF66_YAqN z9~0$9HGbCUo(vn}^dj|hlkA>|Y)gMxq&P`FP!$~Hn*Q~O!JyQyO=68_cG6pz zS3o|eOT5%QRB#%(dw=tH4QP42d7j_g4q9Hls+E>C0G5EIne|aE6b9#q#h0w)`#%5p zoIltG#uU!09F4UHtc`s9yTuZIF6=0EB!2(yEW4C4aU;bWl;QpA8TEZkWbf@wdt0N3Sku01HOZy5#@Tuft|h*! z-q^YTDxKtX4BzXmRfv3D z;vXNIlZLf$e^3q%lPDiJFz0GGWt)<#Q=I zB##i-Ht?LxSqp8$g9M0w$9HeC!FMxLU>nq4QAYk>!3+0tL!(TwF;+?XEy-1;?#gVg zWl})&R=+&kAn{8@{n%>Acu~`FbL4(!BL&@r1=3|=l;QSuj|=v*1~@<2q2_yE74r_R zqVM>R9l5=RgW`2*IBm1aJ-ARGUd)?sUQ6_9zHHWt>=8*638AcfK%Sp?zZQ{Ol^l>^ z94FRwj}HIZPM$DKBpgq?LHm+P&b_!blfBmn@5jdRuVTJAj9S@@{GKoY7M3Xoo|{EA z^sL6uikxOR%er4?d5<1`yfXRzGx6s+9rj)jQlX;RuC8q>Yz6Vbi3>~oH;FGaEXetc zhBOXdIa+ZciUCg}9$)U+V}|U-i)$&-mN3;Y++6e39E&|E{DC`3KS8_qulSE4giW3N z!F*^1Z0D-p#``f~jkceBW(KK;UW6Rm+f4Eqw_~r^+H0YTPWDyJ)pR@(6FDBoYmXb8 zln>~%6K-AV#rD1LbkNp1$ik`L5_#_TXK@>n9Gsex@3a+}yM|WDy;i5drUU-cr+$%q zdkVkeiCQ{5Ka;^QOk$wb!zr$VLWF-5x4!4zuq}w}EtTEBn*q(-$G<<(5QhJbTL}fO zreOVpcQewgbO^TVzjFpeuOMi z^x|gtkG7L@Vp&mHNR{{>8-2{~KW3r=>$y4Zk5;H4aKv1W_zTi_nd6BebO_KrQ(u2x z7WD>>mbj|}tUWjB?4oN+&JCZ4q7;7o^vdNzMyM6a#SC?A9|ZYq*;|UvrbRHpMaQJ> z&V6mL>U{E@cc(GLvx%AKo-l^!PGiH?jRufD9Ae~?D-PSJ>WZ-&7|1?6xKWh&G{m;H zdUUlZLw&#Dio7&yblLp)znE+$B;R*^!Y@Pmk9Fs)9&vIL9b@C`Z(Hm@YlU`p9PyoI zi)gNQk{7|e+xn>~*A(#4%GtrPMJB442!tpFlKhs+`|`=a3Yt4^X<9(F#>p+UHcJcUB5+DZT08C#jz zuLQoQIA0GfD??lPO+~jDI>u%Ra{iLAM3!GJ_c%C)C?5>I>`Kou#~+UT|AZ-Y{N_II zXOn4(saGO)OzGG_*GBbqpUe!gL^R9#9`Pl)w06+LiN1z0PrimJau9CYu7VTsL{A^m z@w$A>8mg)^9zNWq4^;s;y_(Ef8dRG%Mc^n3EEBJRpXKJzpaijlK|n zOI?1)(oZVd+_X1ce5^ZcQ!SI2or7q7ePCL4RC3k`_ZHgU?ukV zklJN^aLXadbxRxE*U)@&DAgA2isvgebIEm7?-V=y+5o&l*<4mD5PoQ6MgFh#YG8Cr z*pVtto=^9+eL+RUUl?_^m3CJFVhov%B-M1yV6MPu0BzjkxJuZxNDNt{+@(CtDJgwWTN$-A z_r3t6Cmp#^`B4Z)ADcO9EK$JfNZ6_Ga{3TxTk$5Skq$Dp{|;1>xj~ibr1nc5buvKm z*k6)4K^eRAVa8RH2{J$bp82e1fmwE85rY9Ff7rmGwEn(4EIfRmd;I~?xqM@HCp$@i z%ggEyx);PD*~?Dq9GUaSfA?~i-${NSbG|G3qe;KF?aD~vJ~>85kq?v8w?z=)6xeqj=fXfE2Xg8cWcx3lWLfLs$6weRS`Cv z?1BOE{>hu z58X>rcNC?u+t^LZT$&qVboK}&c@l21u;T5`SzBDgj^n8js+ce?b6z~k78j;wINd}^ zp8UJM1e1*e3Zm79H(Hwmi^f}#9%pfk)|y`T;Ho?}r1wtW{7L5S2BHT9h|k;8hovH$ zJO^@t75aLINDgVPmbK2E06*d-|sjX}pVA3O&P?(4Fo3(Db|ja|vB zznB1wPJEucBa7BU!w1D)lX{%@V)m}%Cg4F4Sgx~Vf>J9lSIivoA>>%hglg$x&euEB zKf-OWX~AISY`htse6i{6m#+-0arJr-l0eRZ@Q^>jKE~L5Io@jTpe^QXXWLtONgl49 zddGSGCh>*Fjqzp3DwABHZ`HSNW-xVr^9QlpgtJNu6RLPg_@6VC0>>qYe@=V#lRJuB zFwat{5w$FU)el0ybBEJl>B;pMc^9pLWFuXR)g-ZPZa3GOPjb-E(o?g4CGmBqYv`Lu z>LZ)lYq3T)8v5`~e&aZ;17orSyEKOgNAIw<*RMuJ3^5%!6}4=K)qHS1WRc|A1U^vO z$Blta#Mbb|aVnNC%XcbnFu<_e-yMvyn80}NK4Gy;eCi8)i&`s0ktJfc2fMu(Mha~h zd9%j^#C~PDnLnVxo57C!%1+Yuwntun+N=hde8UyG-HMRNdAU6FjSVC&IB#&$pu^4F z?oQ)K2DD`vwRhS9u9Me? za_vM>vqozOOr&eB9wFyXQh3bSID2$D&6o?=tPXrTM6FJwD?>^4zy4mrS6@0^vo~@* zxz1aD_^hJIqVK0n-_LPGR~)-z!n#HexKnB?oJrmK?)*{OJh^Y=>M|rQ)hGaq@TYSZ z_e|GJ=p91QI9^!!QNmzogBkifJG?VCoE;`y-$ya)OrWf&BKVEnz2@S`0 zNM3zU{6;?&YpymCUW}Ku<{uk7GbVN!HpQpM$F_ zsf+`y#6Oz#wSMRC&2|uM(cyIHkpO^Oyu6$;6{OcY<)^M9*Yy>J%g})M!KA(yhgWC; z`wnmma-QT_B3q{O4OW{C_3aa5t%T`lal*?gN{_ijGi_?!}%MED2k=e)} zA?`F{BMZnDY$3cP(DGPiuMPWd>)+~pK*ii!IqeVYnOH78F>)ZD@c0vzxOSFMaazwV zX-eM`^?D~sRC(q6o5x5W zvqGw6Rf;7Xx|93rt{MX~FO2tE6P=Tdk~kCDpof9yMCY1BNj~9@+h6I=*6`q|xor;N zgy)^0KI_VxptyMmIe?!kuH7bmSIM7>(xw+5Z6(~jZ2WzH*47?_XxoB~ikL*X;`QFV zK=Q=0PZw-6S3qdZ%jv&6jInB3A|>Ik9zqUhpBBl_vz-0QHP-TlqNK3cIuS+%pYrKE zrz#%wLDlSzZ&G+x-1g5vgFU)^ijY{vW&!ctr$VF`*P?4plIBS=ABmhix@Wt-0^%Lx zuhoUb7dVwyaQ3kn7A~YeZztf^ydM@^XPv09{owkTXCDEhqcb0|ts?%kSL}5uE7kFN zd~nE*ixlKb=;D`DmIb~q&egdaRzURe%~yK<*n#r#S+b#o)OBoc*sIr&zIv;%TGy;8 zgkN_by1~Z8=-8JY&veb<;>#~dR|r@6{ehkzq2&GGct%uuve*b|+}|{R{Kvp#-`ehZ zGM7+Fk$ApBoZNqP4IJ;qt${`1fK$LsOHEd5nB#%RbO?#A}{zlhFu$7Ehful>*>gn$K#V<*g_$= zX~&>pX*cYhwDB&{(NuUsU# zgW<@)K*I56N!`gCwW5rYF4Jtp@yCEGx9nD&kxPHqa;t!Uuph!FZsEyVgI#JQWwJPV5V-7 z7L>12{Pm$ipUgeS12;EnVE6BcBtB9P2J`%R=AWSmk0nmn*Uk}MaIe6tgec;Bf6sVp zeBToE79^*kUmW*0&u%Zdt^`57Vh?KHkUBG&mF zd>x;h8BSQ-86Hs3hUv3wFBt#dKc2YP$e^C^-qe*YYQEvaZd0B|-i}K(7(ODhjog<$ zMzVX;&8%Q4BYtaX4w;iT4<6_FNOTPU;uOAdX{ge=BxZPy==|Ti_$b9PSmOL}LHDpR za`oLiOGRONGadUN=YmrGFermQ3Vo z`C0M%fI7q{FPp2c{zLH*jo9Q}KVMV)?q}CYVk4?5w3lW#pyFYLg(QrbTaGxCJxs?>a?5RoQK2tIWtXnAZ`o%1t>gr_yhCjCjSR#e5c#x>#D}Q3#4GTb zhI}}D{ogMZXn1g7#Wo!#{Pez)xDJb>&pKA%o*eZcW)KBu|86gdx~_+5yqfjb z!}KuYw2kyZJrl$q1GhA46=CAGFaP%<(kDDtY|UCwg*B@Z-0k+!5al*q8&MYqkLYun z`R`S6i~ieq?oPr_kk4X~3N39{INoVJm;#nemz-*E)!t-byAw-jznt$pV08H+x-5N?H|cO;U!dGL2cX9x$%VW;aSTXEiJd>e4;J?Ei|A-i{JO$s%m@RN zt(M@7-$q5x>XxkPdRxqUwyw-9+YIWBcBteulDz;s7AJ&>u45^5><8Z-CUU~7Bl1pG zsJvx^!juCQEsQT@m;NLBgfv^XF^FE8ZN=FY?P&#Vq8cZk6R!O(Cxe8J2LO)~+%IJj zPMm?pm|W>bCb)KpFMg^a^Jm4z*LPPlP@y!gk8_I5GtNn7eC*eTd&y3hjn^@8+X6q% zsnhYHgkX(G0TWrcx)`JDr)shb`rEI66u{p~90QtKcE}Q=ne*C96&r?Xe@(_)1Ao9} z-oFP(uFY;kX5~L?+)(Q~FLH>Ew}*Y-1q~TvfJk;je~k>9clmd{{7!gqLZ52l9_YcK z*luc`A(@ljk7xA>vO%lEoA|emY2y0+CtF%gDfsa5nQ_4)Tioxa|Hmj3@YLbCZARBE z&|&%@o5wXKaG4w)k-23GVSnOFLf&dXLEDx)dvyqZkIJUOO?(yVYcBTe&@jgN;I}fg zCX$y>%3q|JOA^C(Yv`tyc}nwb=1%@gRG6<*(WR2}TyDUqt6Y@=6Q3Pz-(OOI{SAVe ztHz9>JgDHPWSk|i*t~5Jz5A~wc)cad_D^)UoO!fyM9&a`Z8_E7))WICsy%#ul z_0vHI&A_AD$8+bTK3;M&QDCJJ-J$aN`M+eY_)f)*s##_W#S5*I+x$(SYG8wjsG1S5 zNW5ybd#?r}9~Y7uD8%1>!d0q`gK!{z{vySgJ;_&Fm&eNLL-2HF(xnVZv|&r*vK!Dx zt6H6J8%Qq4IW(m$`++6q4Ow*VQPqL@l&7a1NKSausrR$lZvpD>bdH37q@!?p=*6uY zX*l?OtUB9T0j%!gpx`zeXt|_&`TlXTC&6Un$D5y7u<-Nq-p(VI(7f=X^P)5l1fN`0 zqLnrLRX%CusoJt zd>#HKRT#Yv1+Q?wt%vIB8#iqGYmCWmhT_`00gIWLGb(q<9B9vF%PsjNw?0SD9nQ0Y zXeq<1cLz;@#d&K>k0BjyX_iZzil<@MDeCHjeE>0@$@2f%Q8410<<8JhI?k_O+LAwK z1PP%%+8%Lq==Yi{3wUjh1N zME}LU`7k5@e-sVLHkMrCLqEs)c`|kiFp*C?O6Ms?3qaN)!z4&c zBv*9Nz~O1;m1ANwoGkhm*~p~h#&VvWN5xEW#;5IRauEZ4>(*smnIrzy2;P|w4@tjx zevChd+?Pw+Lra|!2oJ^e@GUJ5DxQqrpQsZ<`0pWYzgYE1j!o{s-antTF!C>FaP1*7 z*R~H%))>))yr|ssr2>Xf#RUzb@#LKRrgedb%$Fj6Y-RSIBJ?WO9el1D{3v$s%dcCG zRE%Fd;c}JmQeD`lUv59cjllu;H*b%#A9;OIfC`M1l2Cu~n)tSMcMd9p|H zyJGp#+o}K0;~#VsOR{=$l<>-pSz_pyRA~6z@_J?5c`6j183?X)B7Un}{vGcR5;bU*2~4jXIQ&TspeBfxq?q@&pt2KSoJ?RrO-*{XFB%1{nv#sy+^b7+8pt< z_y_i1`%dQ4zwNm9{nCZZYu}H`tfNAuP?O~FMhk3HEb#8qvmp9mnp6nMad=eS=Lt_D z{w!X;`lx9g$Y#mq>uu48)b`6OJhc_z#h}Wm_i{EE74iDW=`k8ieaW$PjWg ziEnP6)kXh;1<7&gk7_ryvtbIyp0ZI^!be#nv&B(d6Ij^UVw!f%QQXP{V_x3X!u_FN zLV2Ps(d>=aeOT1Q=1*NmE7|42q4AuX;Su6TQ1i=8;7k{4)wxVHj=n0d2QENDE=5K&0glDO=O{S28AB8KjFN*#k`?Mr~igW}} z5y!T6Q`Nc9`gRUi-5lxbRT6xT+5^}dt=_F3rVD8)>2*=$JzW0I@3n`j6h?macD>0@ z`lQRcIwRk;A-<Dwv9YQm$DeY z*slO$oxbToM2`uMx%TlpLmgz9qcY9GcIaXEPD9K=5yU=fFve%pNG?vL|2FxvOLH3= zH|eX=L)h0PoAcu&$Co+*Ya0x=jI5O+pLeOa?}R|4B}n|N*kG<`fWAFW2A1437|LUN zskMS!Kk?VJTcYy7_dSd<>S6?=SA$H4-5B^m>V;auxE=Ck^vwp6=kebF2X9$~A&w`% z`KL#?Pd**-aYbHA&~bNg&yy>fz^%>wW?h0Merz4z8%xIK-nuH$dU90sc`I=DDes?} z3{RC`q3ot87ux(`{ZR?r@t@+c2w8Ju~%MdP~(8TO3(Vuhb^?ald2!krw3p4M?d~FtO|U25uG^GM5@(;!ZbW;Eff# zD$WxfGdAyA|Itp;-BIS+iS@apXaQ5KwDdK3UwG7Q+u(oA2!1}E5j-18Ly;}jlDk@H=*;}`^eyp$ z=WQuaSrZ|KS&795MUPXVnVjaw*h!t0Ua9y!*c>|4{li8d&@q_bvTphaH+Bo@oGmCM zzG}sT%l8SO7a-%!U21*d1KolmrjM#75{eE>?FL0=DzOL z&kUeWZI&xaM*)1!h>V{6vpY zIdenl8Q~`6$h}cqp+WZG?Bee15+goacxn9~A4<)>$@w5aa>a39(%p%k%beekbn6H} zkMH~Qe&bZ=D$r}aK4S-MCfN!huZX|6^*)}sClBep`F}dTksNck_GdEIL_S&eT`XI) zN#EXWbn?w21unR{Kl|2Ck#o&EneAx|Zk(nzX+JHYQ%v+G-!>aGGAXayZ(;#lT|0$& zS?CBX2b8xZ{I0Rfdf#0!z5lI^rAS^BU=VX*AWc6;% zJljZkz8S+?wa;2Ws={)7))`F{obrfw_-lzpTn}4&+6_?cKp(qa3^`}&f9mSK6M~z( zD_-#sU%>tdGfL%Ya;~u0Z)|WNIga_vDK0z$SR%(c?Y~{B_}pC||BArrtEwy6s2yU1gs? z^w3-k;#az^7!dwZW7fZ_47F|W+EvZRUh%}QEBSe~jVzfvjv7~QCwuQygm%iGBf7E6 zg{HWImm*MDRQKlbh#0QviRF88T7l#~(`D_0W+?F-TOuMW)N!zW@d0-^gUqAYc>I!h zQR42mLvPxs5I(b8Yhb4}u5~NfV3}`$^PGmFjbuLfZeVxF^}7`OrM^$(p}Z|--b3qM z{z~|TC!_H51se>WxK2r-%V60*50NSd@;-HY#`Wca67t)$1hNK_d7N)^e~zRznB6*> z8gteVrj+xxKDuv-iZR*TNvnvj=39_B)5{J%)I0WGUN*w{j7Qu2tW-d+Pxm^((?hAU z^vR!Fl|f9PRZ&+{8-7nelwypT;F;F{qzcl_FiiB$i<4XBp!S|k)(UY8bo1JEmTwc` zPQU+}Qy{B{!DYvu29W3H{lH<+j+V| zV;9l=T{gbxe4r0j`!;yrpHjx+&Z`d#txaKOh@yTWiS+yR7vAQy+JjFqtyT8yZ%QQf zWz3|I5b`-+&wgyK2Ulaq+9z!oC}$ST`?qnHvPPbF+b$_(jCY9JR3JxqI%(PVj3MF! zkhsxWvSy%W;*l`7o|FZ8OtyM0=m|pf1?CocSqosRYkJl}d}FCRiVm-ce&Q0Odg6Mt zEX2D`uc+8DRFisa!o#|s>{awA`fqMA-zn|IGN1#iAQOtlfh zZiBs*Xf6TC&Mnz`gEnA$-6`nj7dm-gD*g3lqky-G<;=eC)~J_uwAlBh8Q8m=>TBAo z3+zgjFL*avf@5y^`+j%_98_OAVX9Ii0(S&!_y=7eAd97=MO; zS8{_qvUqIf$}m~3samhLv+|5Fw8j2;%}*6aulmIwI_`9^f0CD_yGVz_GGmG7ju2jn zn!;bBfo*f(0g47ZpRkL&C`?Y@k9Yu7K zuK)VGvqOxqA-OzdkFq9ux>bHUba%`l$)bz(*@<_QEJ=U1)AK1_wE(%# z8sCY93g@3xt=jR?0(~z1e5Y)yj`7ESMXQ+fQVy}|S`6lqd0eCH*hk{;;M=V~64gPt z3^^{<`__^ip`L!pFQT{o+pzB6CBg%mprmNnTp)h?YVCw>7d=qADkSSLLwtSxe-0n6 zw1?c`zh{D->A(_JR=V}GDh8jQlG%No=&ZVBH$L+JqL^_$?|Jjg0B@Ih@Q1yyMwULo zDub~48V$8N-ZY|Dig8)&7Ge_pUL^a5jFJ|yFR`@l%r?VdzN-g0`1P?3Si7D(OJc$K zZEyQ*EYT%{?M=cCcJP_vRr9=|hRz?qxy~iopq5OdZQmC%XS&2^!kI<#!y-oX#54BL z^GtKwnZPlwq$m^YE5GRyZRZ{tHo?fmztavqGN7j1U-jOB3LDIn-W?t%81+oi zp)hmu{ySH6%U#JBy5Eg)yt|;xTR zJ3CmPk^QV(I#0(k@|l1#W^;XszbFevT8(04j?&FL8P_Vsi5{z)cG|{MKuct*Ecl=y z=GbdKVmwlSa`sy~8cu4&r*~($dk`=cdI~y>R3I^;?r5s4Iovd>Xs-BZfKoys37=lu z0YggCaB4g8$-K6eK4D36h+g*Ah7@CvQ@`aHMD&!oEmE&SZ3y@3sdLrcDHXiSm9AtZ zWQeKt=0{hPzO-9k;GW9`vPVLCHD}&&J1l=@&1>#$12@Zm)Bw@b%4P2)B{DUUMbrOj z)gA^8ZU3e;FHCYqde7Q)=JYY|U31>X)x^hfQRwFIlP@*I%=0b>OgM(@e?@YMVjr87e3OYk;N|n9=VWxyYjAUjodDsnB%h3R zcb7R9GW8s5a(8=-r^id0@nGDlAvc6T!11GOe8 z)5C^jztE#Z_aRR}#jWkT4&SlG;R^p#1!FcS!}Yj@N_Z@4fk{=t5=Kx(b+Rt)Bj-}s zH{T=4Q#3E{x`p(mc?Kt zjtXpY93(o0hS0SOQ~!~@1~HX;@mMSv^%N*qakEZ`?%4~Ckl$tm$tJ$iqr|d z?YHdb+tr9a*f3l5o)YS%Cmq%&{l2sk3;L8<;J;7zGlNf)xxM)O>DA$8WUu7T`_H}^ zfaQ@jcWZWAV|mJF_Vw9k? zJCuU&wV$G;2L&f)KKm6vQi4laF6x(iEwNzj{q_VJnRii7JZ3eKMHYU){oAaDDYvX8 z%)Yq-uB+V>yg+oi?o*E2XFq8{)jG5Ka|yZ_Uuh$Ksmcua{(Ba#k!%h=m8{XFTS*@g zdXraw^&};=X>zXN0235s^JbDh(;#xLtTwG$0s0L>rk!8fLABvgm)<-2=-qgyQ>UBw z6PxsCJZDHh7iQwUSgeTCI}Gbu(=0G-<%epe7NXabELs{(3&5gw=$NOO9scwM z1xj*wF3)%Y7>PIPrY4d1Wuxej^H)3Ud7sqt=Cd_^Jo~~dWj&E7DLOT*NeYNk+b#FZ z|Ej^*?QeZNdnqjT#eJ^w+7RIN@b?_aH)j@{^jT*?_^w~JH+~S-#dV<-sJPCa&YMZ^zMT-|r2N|plkiBAe(ryYjw?pRs*_WyP1tydraBV3(c zwHav%IV|BAgJn>(EVyMDaSar>9Gpj)`eOZ%&cxe$>*YFOvgGkjJ_e(yoy2ctM85-*2pINlI7Xmn^edi*QXgz z6SOrjp5^O_Q@eU=#%jJ9VBvVKkFfW*7;7(56@QW!Oyotza_Kfk&)&-u_f9F z;_DnQaPQ~Cf8hsPo`4COxBZWz^A4oCeZx4ivUj#)Z;pM;^El3NG(;H*Sy3t#r9xH` zBBP8Xt0>Y!RuYd=Nmii{vhs>Z5<>aizrXtHE${I?-{*d=`?@|?`@qZd2&bO8@p3H- zbnv0&^7Tv+B$_(ygyE){P z(q-_j$Y94A@NFTL*M-{iLSQ5ym`??=L+?(ELn%og+y9C7Z2Ce%UNviMcg{#6S&uKa zGi)Y!jlDB^iI^e3DBGC&mh85DGK75uz7I@R1Pp$cTL-f1 z+3FCMSX4Z|UC9!QPCRPZf5HZ@9^Ttx*P{kL5n<=j_tIEA{+a3oHw{-6>7EIyBcaBY z-}>#o;CGGX|L!+yidxtr4lRNHi5-vBCrVJT;qlFeWmPJ+PkR0OdokowVy3AvZZ>FC zaIt5f8JtrGt{2bfY(orFSywE8Tg^^4Y0IipMn{EuNqZzH=#}=|!&1l@6N4`09m{}w zp`4_%%_f|CQlq9(rqHAG?ef(S$amFu>mv1Gm*x;Z^>-%7PA zvab(s5vv8iU!ufkvJV^5O`3V;jH&33)q|uTkUy^Dxt*yEesg`-`Xwx^f|V%3f_^8h zaH*(NyhR#9(UNwg8@H#)W{1b>b=yu zlh;iLpLIMKHk%C`|A%MZhMpI~kz&KT0$ea(Glg}KVia+%XTkcKLt@tv09?^1pdGYrp()u2zZ zI4jXXH%Qz6?cVWr_&r7*-%@jq1yh|YXU@3mVBf{qM}i~@@cE6nl$yX_ko-LRFU)5t zf8s8UtF37JOxgSsZu8eSU-%4z(9?IG?dLs2M>@6Pi-Dhlvg;X@tg;(HD zVPVa;dSPIRTVC<-WJH>wkVCON=Fe2|!f~&$Z4Io5P-DZreYYfT4O_eQ{h<|F?Fc^+ zj({T{{-0cNzcQjGzxH(xQA4A)4kfD!aL)|d8}+wW70qN>zdsoczS2vyLn#U-NT~5M zxBaj=`Y5CQyd~5KC)eqZ``**Xq*$!@tK@vxl7F>@Tl;P6HRpd#B#_+JdiKAO(7)DBuY<^Rwg@IxaTl{Z?K} z1x{pvK=Yyj@<<3WI2J;~%Zn4ivO%U;RAkZe5`?G4!%B9@ZU>)Q3ui?wJM8iM&n09i zLmy+e9>>WfHAI}LQu;V2j#_8rRU+%9aCmajv7MnNs8>(@cit|@TNdl+va!=JVS7iA zN}eKW2$C|`@>2_k7bW6}W({0spZV1p_KFJv9C{-rJQy`Pvv@Wxk*_sSm7C!F^smSY zEP{QQK|xeNU#=b|y{jJRpM!n^m%O_xfr{um+4IQFDB$Vnb`Y9GNl5&cVy&(s1qF}u zC42?LkpJ@$#XB0{W4WFDy}F~XDdEh%bKQ0PC`|b01vhy9tBPI7dxPXK_s|``-I-?i zuE2Cm8;cE6*&XFyOcKTJMrXJzw*aqw8(Wak8|Z~T_;_5D*#z<0lE2$!T0kHC(lpNk zoFf*;$2mph@o0}W*V9_y&Ivp_ar`g=9X&uhk@Ce96EakUlN=Oq9o4|z`5JI;l&_Lr zOR0g2YI~djaFXmNCayHWeMj=0`1?nXp;y8C6;l;670<}*yE3I@gGvZymv$rq|C%q# zzxoE~ng*?``awk0|D|)U;Thn+8VfPd=)hHYmYH_L5_IU{C?#+3F_fb_^TTKKkgxcY zeSF5&cuT`$=k!tvz9LkXwf{Qw3GtX;SS-{=UJuI|uK_2Jp~39H3(4*#nOp1zjPZ2b zo@F4D2=`PqAN6@kJ{9*AzgiD?4?HZdm-qFKTcI?a!~Y#BH%5h9ZPd#zknm{SW+lhe zT9dVuR_=Hxa3EgCcOu9YFJ;DgE&Y%O_3>(~axYM!y0j( zzbwQ49`4mFogx+E##l5qMNlMN2?Y?B#6A$M(KSiO-M_5hzI*<1+UB!OGJ8{Wov|nQ zI=Z8~|EOAGUJr4VVMN1~{P#qzR2yOEAEqJK!8fqTHNvi4iV*uDb$+LH0UU8_*s<)f z9_EQK7?k{Bh1lnND%_jlee|Gl_I8ma8k*fS!7*}(ufF5;ayA$YcbS>ze*M*yFeCG9 z|Kw6bGjO9c z|9y0Ua~QYI5#MyrlKuOv-bmT%LSL{4UyGOq@aP?6wvCvhQ|x_9$Hw4!QFY6?J8z9h zC$BUnSjwZL@;X0)U~hx=hK~nb;K!0V-v!0(El}KCn$`{A6oj1r*ny*DQ29jfS!yjA zHPc=zu6pX>h=%l!&Ro{m$gAD9Z;Ku#Xxr>Ka@Yz#IAYUxI8+}O%!Dga^9}HWIJ;R7 zE*HrSff zRaWc({W{zDs-L`IMD53#PM?L(li`9`NWoR;X z>CBOt;CG+oNg*FfQb50*(D2w63tX2Vboh3KKKlIah>Z}3IcAVhm`jS_!{4sGzS$h4 zi^YGfxEdQ1afs(@rJnEVxUh~7J=;b`F@E0zQfwf{o_m=`f(QI95wrI*M=kN>j;aGa z?nwG+qS0r?f{-lQ*D8 zgkQ!qXdmdy|5$$Xaw6#0*<5w((ZM%4*gl6hL!QCB`q+968I8sXo*n99K+(wWOQ)(f z9#yX^AP>^8NW9njmQUtbL+-yC84K9cP7dch%QHpMiQn(26aSLkw~y+~{xU^JhxU8K z(yY;yD=kWK9-7#?Q7g@4iiXeJm|8T;0YA7?#S^xC;5XVd4Y9^TkG;h^3BoOY6zqQF z&1-*kl!bT)DhCm!xCKxxuEAU~%~7=VjSV`HetO>7hY073ZCWrDIEPoKCx!{|ept3K zbAJQ=HgDc@_Q4BHj$iDOAfSp1Z7&`PA+n;B?U(O3Dq7;h-1$YFmPVK_GN^teNgnr7 z_g31HjX_^d9vcSDEPE|y<88bawsIMr{Y|4`zB_duaD76B5pnSs+XeBXfuCvV&EP|3 z`6vBH;w+(9#K@PmOG+D#S1`R#`mtp@<87VLk4{eC6V0-_cd*8DtVV9LXr>}DA zqFa}i>c4*k{o>-QwW)Bp-!Yw7x08^?KNAoC`{oPt2x&pgpy^MOy@T1AW)5BC(iI)L zVJC_T+50c*N(m!T;pV+Tl{AF%;taiYIZ?m0?5j)KirD3$R#wWADex^Lm1eg8Y+^_n z^vV6u-4tuYMxQ@qgks97g+hdM@XS>5@)ct`(qD}`BRgt|2z%O_&#A%ZzhC2V`3^B; zzq3sBW{E5g^a!ftP=LM;A8Y^eR5-6YuMVppQb2{T?1b~%mdJbg_|Ba(w88NoeQs3l zC*h5!1ywKP=qN#WJ-`Y2s6q~IX_74$Mfy5NPMxYzLB7YHZGF6{i*{<{8U!CQ!?ptP zl$BK2qZGOs`?wKNuiEGHD}O;R+8P&KTVRY^8lrb;{G)?`Wy_$!pI>Bm_7;b;kAZJD zb;={<9}P`qcvnzpR>+gTH0jAs8?0XRhQ;`XHS|jL+kV|03Fz0y!PBmt3JPADNjk+6g z6_5*AeO+W4%3_QezH7ZaQB&RIeIP)h^NJbXxhMHzR4(XWZ*=x4lNB+$rt}wM4|w0* zHY`Z3GC^a<-XoD2eSDK)Qj}W~bo;yleuWEp^e+LEa zxwJQKV}gp0-DN5JnyQc4^9k&B-!u{XagDMC@Vob~3i>`itdHLEK55?0u7_VmLApH* zd_V+|iUjb#F?h}f)xMZ(a<5ywvusU9i9dCR6m2b0?#Pz+H;$V?U-zZLeLuiHa4b7p zM99O3kdy4W){Q>sS;DuDsaA(<*gfB>}zsm%uzcBzJtQ1{5Jc2Td>af zRrlU(@XNJKOD7YwQO#9p(zbUrq%oml-)jp#DCuDtIdy)tdrMr><7n`Um}k9ZeXWNB zm)zLj`I{m`#nt(paG%?)cfN$42Y!FsvDhSaWlYks_;+7`j#7`=)&8bhVv*+T!Dn}M z@K5RT0@d@N4|xjl7M|3^!f%g7{aK^oOh2>f3~2-GoqcP+X#jA>-wF6~>nos9xmR!U z^;z-V7oCF#c2KZa{Auptt73Rp>8+$MUV2zRMY)Q02Y#LnwtZ}&3P@5o_G4`U*o11m znD*YYLZTXFdzKQ#aO}kkls8P!uXl|rJ6%Nvy;xC5biQMW`PwJBOx|i^m*T>+b9U-D zBzopv?@Lo$?HF1=kZyppzW<5h`G-(bhw~NQD@O312)yl^;ki~YfNWUrICEU2{A6CT&eP#4l$F6DNJ7Kqc--_trLguKcpx^Ak2kZQI zS>6)u751OzH-?^pY%Av%%o@0nCAV!4yyx=|oUXa#OvG!q|F#o=qcwDmKhUTfa_y4q z-^Kcpi$ynMpuVMm6XB~p$0^R4N@ zyTuH#5JOzntSK{}M8Rq2jLP2amBAtE!;B5lrugyhv=8O{z%@*WdCwfBi-_!2$0BS+ zk%rgVSQlg1YtH@3uy~13v_)X|I(eXpy5HZUHy?}vg$l5cH$gt-AwNC4bl_f_Pp{;W z29u`7fK-whI$HlcS}%@{lFsL~t?`)P?CW!qT@^4uuQ0JN_bPzC`l_>8nICQJsp==_ znBf>@zJF*6;Q-E_?en5AcPw!lQEyoz+lTkF=JwF>n&AB2ZOT@6_!d&@|z9!(^+gSIUmySr9JV7fBRVA;sd2JXN3(=wvzkjA#KR> zg#791cLx1uyVTuCdm)r|>XP_}9jvJEY7ij;-Yb|}-22)M@Dlb~H*l&>~LIl62-lu&t&^yI{Bapp|f)Bc<&OVR*vZH+T# zh5vI0u2iW_s@vgLf=HG#@TQZ^YOsZ_8~>2 z;A~y+*GUB}Y(IF(9XLK{=GZ*}KT%Zaq#o%BK74z%8)LT+9cx9V%estP;~*Z)GsQQUeSXK-4UgU zhFNm7RpgT+*pB?6AH-z2<|E2ZC{Q^?BzrE(?gHKc(cu?%}vk7%vr$1|`{GTR1 z7x#XP?zR7O|912qg&w%l?p0~QI!)-$+&H~fZi%}ctOM z4zgIGY5ShC!2%M3jZf`C7pQgc^LMv~IiX(t9_WIG!b2(@3YM67N@w5V4MX&x!s$fVp9-9OR0!%7-$>D6r?~ zPwsWL!76$`XR_J|sK;VqpkNgAi|j&&TpGp&64bAud2c^@n4;nUiDT6t7irTu9 zrvT=&P}8ZuZ)q4C&UY6`$b!Ez?$5th%;-|i%XT?NW9XTlGVX_dPZUo}Vg&s!*&@p{ z73S;iw_hi$vNh1qdDRC@CGb5sfAux{L+}f{FWJ#xpn|*kK9-$`2fyUr&1Am-;9-mt z-YQB+V?sBn$8kUxsg(CS6}xI8;&X;7<_G2|FjdF)>Jx;TgAXT5xtrmCq-}D8a#XCl zIV5np&HxwBjlZ=N10JF4E$)=Zz}FrbO*;-gW3QU>%C<4!#2I!7x)1AvURxvf`!F9O z@+LF)3_2hxL%2`24iSZ{2cS2}h$Q2e2f65*mb}?YL{zn}wgbf*;_MJ_!9A?- zJ$#XF?(|6?W$(PL(6vVsW82AntLJ|=iSB9LANEfLMv*vZsm`^ z8ltO?!n^Z*^dToYxsN_rC9jX1sU~Xr`erzeWwhg9st#UwbB58eDB`%6_;sa%52=X@IWfsG=F0_I%|rN z?BOj{8_?@UmH61uMn!FQza#=-zim_HV-_f2g-@lf=tw;Tu7A>zUm%teiZ zZ(?<;ReF7z84fse;A@Y80dSelapt=MxBmSYK`jY5c@gRM6a8k$*R?(*z{?nYc>eM< z^`{kDJIwAB&j$T6|J`GcX;i>g6un=o4P>0ZTbrxol_h2l4fw^y1Nm~H@$^Y^Ito^$ ze;WI$h!UEGzJ{;pVxnN3Y=W5(&Zu@Ln{cS3>B*ip5&3QCJEf+}fCu`dLOuGpGEGqe zWx2mt34EM)7iGmqETS$}t2qK$>q_=IL+vLc2U7d%0IKy-n3^*PwlUGkm5L;+%03Nd zSYj@{ms4YNvWP*oVD`d38=Q9c+tmbDPT{CV}`+!$LBlB$@R|6YPOO(I3;c# z>7O;j`M=tCOu*k`uQD{k|Hu?ows`m*i-vyDO}Urxn^KsaLpL!!4}3U`+6NU&fxqMV z$LOIC%tI%%{5=;eP%Dz=e+_(<;6}e7CvA-THPaskmujQdmlss4??WFG zq0rrXN+18^XOR4!PQ!(ZZ~n%XnIZ9km}h!gbd+}W$Mq5}c4Shha=Yi`|9!ta3sQ|6 zWUreZtmE)!-9i4W{4MI}et6WaYH!%Td}Kf7=?4#rkx;;g?a*J#WLY?pMMSD&1G4ZT zU=+RO{6gytIn4FAkSVh@j*ba>FvFydpGJqb8GohX-C17Cb`B)WUOLDfeF^pqla~XO z|3er@-Z*M!E`=4Y?a>H?eMPsWO!{_LRV2xGGvhhrFRv+blPq~4N88GLKq?vf4;V82 zcW}z%1O??C>#%2!_T_%&o-<4?FFp60t4I`O`wqM+x{pw{s55dh(#BU0`?-}u59?h? zQwz-uD&p-wUhjE{f?J<5mc8c(p4d-{G6o+`*jQOm{&m>1UCQ1ip$j>&e3chpBlMAf zrtWjYKQLF{j+4JXZ;7`R?^^ueL&rhKqW-3E8sTuxH=b&t6uh7+df`)}GA_*fTc;99 z#_S>MH6h$|to!kZw8~j){NaSK+v*44$Xp%`jI*%9oH9mV;!i++&GRb9eRpki%IwIE zyNQ{G(tX^-D$%$VYY_K;}raCJJuR5GX;tlOXp-&xEQ47vL$laWe(cCmy`Tx@p8f1lw!FS|4ypfEx9X9?ul{#_R5YDdKN-!wwI z9Mj8Y_rjyRD!|}0HrYf(GRC$NCTRVli1%;EM+hr8(VGceXX_jd(e1J+3q3R!>8PWyc;d*XkqkmlJz5U{7!Vn>$$#{1>htj-&`$ zBFx%QwEuD}6*bltUn_%i|H@2X?lpG;A}UqNPE-ryj~-Rx=J#Rl_C4@9>AVs~fkI>x zJ#I{zla_WCfLw6NwjrrUS~y2a!#b?o26Oy1y&krchVdlKtBicJLqZ^6lj5%C_d8KU2s3f}}aOH2A`}60M zQJBbBFpDk~HH@BC|9@T|>GTmE{E|lgZ4VVMO{n43he?^iQxx1Feay)4lNs{4@S&hu z%nCCPRwW9&B%>^cYgP1fDC36-sQEn6vbe>epteva9F zV-G@JFRvwayB_D0}EB;apQA+8ZIhlhKsbM+o80 zJ?b&;a;S#Hd;0+Rq!@HEIm7ko=z+-NXO$Ypxa8Tl&5})1oE8-xMGDhJN3$u1=Nzrk z`_rY4=I}i{EW9BTJ%Vt^`#k*+b?7B5FZtmByavib{vHE3pL}E;6aziz6g(t%gI!~U483Q{%G)7#Bw8u(kkO5d zNf#$lRgHCkhxT4;b_F=@;$4I?7k#u5*3*6%_Kytp83oEt(9?EuNHwV%_&Iar)2E08 zRQIa+5T`NxUd{LWd~(#nA!XaxA}_<|Bj%m4JU!RsT74MQ&7J#NvE1kSTanipwzeQ+6`c2Ww#$1v8T*?k9*auS!b8WnHxf?RAkor5 zuF!)*Xvq4zMG~U{YMoW3z4}Z?kxxCo+JJu;$9Lu*Y?elicay(9J*JN{So{mKdZ8zL z^GfQcJn+$9=dIM625zY5O_P&WdU&LMzVO}{%;#fz_gNiHp#MqR+N4t%wZ&vnMan@p ze)Pz%qK}5U{1=qm=BRj;p4F6aK@0WL*%Od`s8i7NT;<*C5t>M~#eSs= z@-F_W&L4}f0yw79cPpiHD-vZ;eJpy~5^3EWJw^yZxNbd(-1tHUbdI*j4rGMV12rQb z_*!APY>_W%H0bYo{jjaqQy(YXe^cmdD~J*n<_g~@!(L}rXY1bQ(_{ud#e;r_%9@b; zT|s3miu7vIL)eReuiRLCKcSY6M(giNM9wK9-}$w3@t0u!OYIDUPc$(*-#C;eI(4uV$}njY@R{x4W^D z@ZgKXw=O?{yj;TP-!0{e$V)TJXzLX^vSmKKjdMW{Utk&ezHiY8m#bZVC||;bTDbf} zgMIZ;g5<;DZNWO|(Ie)HJ<4>%M_|icfSl8~zr~JE4B&(KI}ua8R|_{C;jN%KQE`q; z_kYCQpf77!Nfr&leY5O_Voo|gu8+u^9jr0K3~mewB;YcX@?T-mXriHvy*bV?&K5{% zvxsmU5T<^SBQpgqIC940xE$|!rlPkw}c1oaju~rfh z&dqDx2keWx;Wc!3NWd0+mZiFRr(De8_muC>sCiBfPtG~bwd^M%yAkCA^;s%z=b#4u zfjnK|t7cl3qY$dn>bcmVOhE!{d*1AthFsH)cl_(#pvPEq&!nU(WAnlv;Nzv@V!0bU zoZ5znZD&N+<>%&DA*t-{UC>*^4z{xE2Ljiacd%meKRsOG>yVQ8z#O*>e{w$MZi!#- zHL%(?1AF9XLaRF?2Vy>RE+R)B=4LV>+Oiz{P*LvwV%Cdf#fVl5(*_$HO?V>a)-p>j z-b%4^X8}F0ibJ@@%No698%_3q1b&C&Qv>%)7RlIWp)>sSWK+xevWZ;CUB2YLGQTUD zj;9v-#a~;(^L)-(++LA@BbRzwF1y;`5gqEMH{WcK?)$0evVXv_t}?MeWBO>g?Z{zg zXA+Jst@?Mj{4ZJM(6q?@X>k-@qZ74St%W7a1f5Nl;rpIiy1(V1F~UA|xz!I1p)agN zQEo5HJ9-wXm>YZ%=A(Hdvm_ig_|DZt7Pyv2vX#boWU*46ck%up5>5`W18Gtowl4Rxu8ImTE& z{G;hAFXWvf%DMYYfZtqDT3jjtKEU>1QcSo3>N}&K?seZBqfj#|m8#XIjc?v^raj=# zr0qN6a2I$LH8fNEZQ!fmn8`SNVt_pPI^^pvIRjiuj{NbZO%tt7uK6jCKyQ<%W?zXT z_{dJzl{|5nBxAomfsR~3@L_*%AG<=uYPaXgBg!msX|>_^FHX?cG?97T(pe3ci;0CD zt~NsqTZ^t7P;-17Yg7YVPVCvQI}ZiQD8 z-rwpe(Zjetl~>kQ3lVH)4!!X>VuAz}k4MBYhqTiX z)5dL@>j@en9DjCixqd9acx*#Wg29rOeP-?J=rG3-pU{WM92t$cL6*a2`G5 z4*7^HEkQ%q)bU1)#v@4plB?dQj_&aTZp~)4HP0m6|Gyox+gn1%TvXO?`kcVkGpSM? z6c9o9h)oD)q`)5c?(NITM8uVIo7$^J#U{TE?~FkjS%xU`Y^;EtBtiohMVCx!!N<@MP$!5U+B6 zZ4^4ge)tQ`46-|VR(4S=_ynoNVEmjF(%O|~=HUW!NJWy<4cK?C=3I&I`bbA-zYx<5 z1|escd_XBe8T0U9{?lG^7t@b8Y6zT3*o!;5npZ@WCct%Mdxu7U@7q4&r^sw5ib@<1tW{q)~- z{Gnvw*-LvmN`KaQbY{j9;|mv>qQ7e46>4UaE98P0G@^=!!G9PM+Z(NG0y+1ZomS7C z2nesjkEX$o^hGj#H*>)D@c6>6wgLlu+a&h8aV8bTRw~7Yh=30Agm|<|P6;Pq>cR#| z3a2^0bNS1sjT#E1*VI5-UrgW{ZoH?1`EtMRQsA~g-{<}{uL2M4!o33mAHu1aM{?tp z(;Lw5TF?LG7_dMd>!fmj@Gl8Xoo?=YYKZ<(+9M>N=;5WGUj@q<%@E(rkt6t>ESgp$ z&b%AbMQh^a2GTB8*ndm>Za#f+Ox$2N=n8#vL{E-4SL=aa*GS=vinK&GyDYwvZkGuQ@`UTC_LYLpe_xQsffu@h( z?-!+54Ri})XD4b*++8YCi4O^@Di;BcyeU!OGy$tTkKt>q(ZP~dp)zZVz{z-8MqaQJ z1m4-{gfouXD7fNM=;mL9K7^Kgo%N$5Qd^LL{OJF=0Z)#n)zi`R+aAVUE0UPt#cgF? zut=sPe8_wQ@k@%l_Me60mZ>pkGVW!Sian3d)wYQq1@wjbd+VJ|79PLQIT>Lywc8-;45Teetvn-pXRdGa?56 zz063(D+fiaf8n*5^FL!W9W=4@^Vc?%z4)X1w7(GM>gr1NW;H|ENqsh+U8q|{_XPZ9l5VS^ti^ubtX033z(~`-6CK2)UyG`=MLTZ7~_U0ZRvMaH1x>Al#34U zud@>^eISp(9{nZh$Za;@5o`beOHA zBfro9Guu&Qok0h1ZyxCWd=>WY?H5RybW7-)3iVw7@_)TNTSHT*mb}z(XHJ3D1TEez zG$R2ANaz}0UPm|Zm*a0mCZEv9m6^TbdRxS@geDy2Nb;8)^nsdu*opBP7@ z8bkd)1DxGe>ozU}92<@K3aSAeiSo!FC;*;p=fydm!$j~yZ*TpN*NYeUO1{+eT}IfK z(Stqdp9wCHcotE=0{qJX+nJJSHC%T^wqcHBfzssb_8q+mzA*ddG-44M&4fQOBE#On z!%(#{n+^1}^Sz-OuE0r46YgWE0Dsj!8AXx;=zeMq{~T6L(fWf@5n1q`%eOHb%nayZ zmcRQV;~l|g-19kSU0DZv(SIIJ)cV`R{B(a|$Yn*8f9>Y=luavqZI^Dy_Wcw*x+|8C ziNOpNRwfNSG5OWh&N!BT!$b~6inmFgbEn|*moC1I`bWbSq1AEBa9<5?zObl52Oo-O zMf#IV^3WG3%3Hou0Bf5^>_qFpDO0Zg@3n{lisVUY&FHj7qpM-j)6dt*`kra8sU~vx z#Fnw?uiLFqIEOZmX&>~0r`k4rzO0C`gW^E&E&_hQUhE#n4!l*o_ihgb_P*A`^AtJY za<;VO+4_ngGP$(g8s{tdm{SD zG<4J=VdW|0NF)meFT1lTqQYCA4|8%&5l=`xpR1EK(iQr<=Y>BFEqi-kzIfdVxW>%l zapGoZ^g&0ejk*Au72j?1@ErwfJkrnn8wL46ZjX5J#|Y^LIkW6(Lnx_*nI87f3bS39 z=}=j>#2>|L7aNTYF+*ttTeIgRnQzHmtoaEmc7DOF>g5Cb^oIjHCCZFgE#(clokGM5 zQszF(k{_teGWclLE;~Lr z2ztVed1kfKVlcO>#nE6cXF9?>^5zocUz&Ea(%Fdk5zl6^5D|JH*{hG{CdlH!liXbb zHB{^4+f4|J5#av&lcxlkXdig4tUZro6vW(HmE; z@CQ^HZW2-VU+aLs&rDqf1X|2D_JA|=1Q}(1TiLs5X^nGF9k@u6Az)H_Yhcb}Nlg02 z)>5&ci}?Dxoz(QV;t&ga#!=v7@;!Q3rz{4X6aD>CF0p|6_T>Ir_SzJdvqda<35!=sRC9~-P9oSzrT02Fq7Cv`Z7xGZMCYpb|Vg-)A*y=#{O=V1!@pcx0 z-n#B1o7riZEU2?$=ljdBPqUG#Y4|1z=WS@Wrv>DSvGdm7PKjHwmBCwutH8xW;}@OI z__B9F01NepfCNZih4pc#uLLI zZt9wtVZzF*jUxou`?40tUKoexsG?Pq>yHSM3^P%1C2Gt1Ha8_OK6hX|HX1NJ83oZcjFF+u-tgo4VjQ z$cxUsQE&u5>1ci;qvLHF4$@sgjXzCrz+ZRGwU0VD;-N`yTRZd_oQsV(>7|AJduT`G zmJkX%Yq(YOtpyf3Ua~1009?>Evih}eHdxs8eXE6^31$&A+EuF${1iQ--Vb-o(8J>0 z@?W&6_$nnRs?CyuRd{N5Du&Xq(CZR&V?TXl!4cG$TCIa$?55Wx%~~R_&?p1>>q5Ze zBRnIAkj$JqKl%Dr)hEGV?P&|YJX zeX8tMbE>UT^5Z^UGpZg|c_K6@Fhs|v_;#yy0GHs=NTPWg_x9;S_zDE$`Jg{fI=p3v3*=%)U#tDR z{Y?OusbxOPN=C>-+P2nzRv$NPd|SE7tcSOGM1_ge1l3QiXo-FLs6j?FxTTV9nxe{_b@!|~7Nm|?fxlC7C8e7~Q&Z*;RG&7Z>&j<>9^ zmrIw`4Brx&AtOw@&7OjMqyOWvT2aN6A&s~D55pczzHY@}w+M3F?Q4#0jj&|O?s4`{ z>WH10^Tpr}ZFGLZCrA9bA$DJDlX1UpgfhCm*mbQbp$F&xZvGWCK;_TR4D#)lf}SN+ zzp|s&h*+_}U)%|O@or(?8gp(PdG;OXnK$n4#V>*;kGI9rR2m-`3gw72qalk-wV8LpSa3*BG`gQ6TgUiZ!5D##%Ik33+}qpP$7CpnMz%4 z3cv3==9Q*|?iPfPRTa}v=i z`g}qWD^$5xWDA%{5U zX%EHIpig3uyP`;m3k7tuHHpN6Ui)PCr3nKfPP`raZ0Bz(o_9<;_{baXQFRZhzPSKT z>GK@wy>EreE_2&(4O39z-bz>Ya&7#}T}N8uJM`XVYQ=m}m&HUax3&a-;NSTsjz*hm zU>3#n*S;%s{Ql+-`W!R(uANoTzb}$VX-f{ZJHZ-F-*o-E`IZX3Swtjo!~nhc`0I!2 zKHx09D_?x^*BUjj>2bTm+~lhKW-A80vz=WSr9J`Q&w%r~>`rbxxMqJW{4oX3#dXph zKEeB|79Y8?^Dj9(F+O+;_!>ulD%R15RdM;Fie<^0!l?C2&?4Jj5~`7U>Rosobj`e{ zL`MkLCy-SCcr(jkf59-jMkO*Tu3Wrh@x~aqL0j`VTP0A;3F6L|-oUS&Sgv6D|G7Y2 zd^@}f`GJ^I9v97?o3inS<#HUmGCH$lk$~SI#YOvS_*%CC9@P&^y&%JkvahJ$8*?^7MB$s|YpEpEjlZ=Im`Y)x7cR^V&>zgu7gVJ) z`j%Wg89aLOvG+bAMHsyxJQw4M#WW1vvc^#QOiX>~>q%fnM=#yJ0mqx>w}i3eZP+sb0Rv0Qmd~@|x=1(A!iXeNI724>9ygu_Y`nl2y(SB0A|N zIL(mmngJX$v&W97gVe#-*H@ge^DGxCeJd9#c@KKc#q`>^-z#Ge0{3(z=+lJ(aY1*w zw2^&Ej*)=7GURs8ZqmDH=qY(ucFQSlAv%UgJAL`NhBg_Ozm7dMF9rQd zV@zMw`hPZ^y+Uc)x@ z$k*}beR~>ejm!=qMgHdosMv@3r-38n)CdH)0Fx1epUb4qS6#G3Q7$!y-Vxo1_Vb#3 zM4Wp$*;$Lj1V=ON$i69!uyyLIB$Z`roboR=rZN(G1D@4IgjHE0hVu<2H!gyHF!p}e zqe5fc9Ti5n#sfaWmX0)m1Hi#dKQZ!t7x<-D_vOw+Dq}SCAG^+1Nu<7{@ix4`0w;sQ zFzIWEm+or@cC?X@{n_QjH8&<4Yx(s1b;!XGS{POTQvx18kO5yMD4`*yu;aA^=r_G0 z```gD%g6G|xi=02|&0&N`Ldnx>KXfs(PR_a+%q^2sjXC9z|6^EB zwdM4a!9v$q9E5e%vEvi3pJnj~_NIomKNOdmJ`SIrWJo^4|Ie}Di%RuhS-IWEBxM77WG={Y0Ua- z<0a+Q-?>?2)IG%HWuVD}W4{!4*1j-=y?H%%j;$Hm!SRmy`Ac)$@Bo){Z&ydJt@5r_ zgWgOvufG=2V9Bi@~zXX&p@!S3%Kd$PD0qhGh&CTy7Co?TaJ zRyTw2Lf-Ps$FJZk>}lt_9t>QbZD*Iu*>PwKH-y(;8}leuocI&O@~ z$`U#DpbwoIN_<~=n~sNwXfq<&3KLig)cuM?Q1P<8S}fQZ)Fhv0_5wFv?nQTa$CNd0 z?~-q>-h}x@EUL_H%T^?7lWez-!5ovk2*3I~RB)Zz!h*_OJv7qr)Weg-9QQw2{xA!i zvK^oN-u8W5dxG6Smt|&QL9lX0Dd#yv9)VqyJxbyEKNXmd6zi%DOHh8&tbmv^06&l^Uyr5F>hdlYQeD*a1_=J~pyN;wV;Xvu6;pwg5?=i|6 zx&K`rb>CoPKDtQ2r3(hiv4ip`v}tJXTB9D$Hvh?RM+EX}><;U1ZRJtbkFnM&e;O_q zUU~DS9r8Y3y2v7Ky2!AN68Ct*4D(U$K3y$<`xIx>rbX8lq!K%D?d}CLbjUs`s1|Yz z_^ks{57}&L)noZbZ#GBYR72j4zo+6Of!ntaSQ%scm+}s02qI`Qn=f+TX?c9&?Qo9Q zLkjlt^v-%rf}fkSK&@9q02(V3s+5kW&F_%_rlhxfamH+4^2M8}$=3Nb`-*APL(8g;<5A z;45Iiytbr7f3zvz5}$wZleJ;%(Kkl#-CDN+--EEC z2i4kkBsH7%Dr6&`hkE8C%N_0D<$#N-|2fomj1T;EG{U^l^8@eaDh)`>&+i-_zEt45 zWWQYCrxFN%lMa4*K%czkZQP`eoV2)?KC?-VqwLbIE0@9xZNX;!*9KWxv=|c+f2`BVR|?{9XKmdsq+WKYA*YAWb@Tk}}}18PO1&oBvcn2RWr^$vFuIdE}FIzi6{5 zIZ~M$tz|>y1}=;7t4HOOOqa#ZsGw>a`fPjcSsK#Q~1X?(z122xTsrMbl55C7e zx&MAXYKtlH%WT`%7Q-TE-0p1Rb~Ym3`J2vv&9n!B6Dvf2{k12o;DDfnRP=ETM-Cd~ z{-tSYtAB7r9$|6bfj0&%*bfKy>h4&HzMoal{|VOrQP9}bIqg!iG4BDl5}FiB0%Z1t5fMa+XdX%ngWgBC*H zm$~SKdc+>L^b>RovZ*-d=iM|LFnF9Wd~_U0%i|*(6Kr7b z3X2F7?|tiEhmEsn|){WTlPUf+v%>)PIsF54I%=p55WwU(B;fiMg6Q?rpRS*z*Rf z4at_!g}f|Y1M~;b#hzVQ2?B*llmiAV7;-l{p@s9`h-zbn%W70YE`60s4_RQ(-RqwH zyqpOVQ7-S`rzwb8Cw=v4cR~Nig!5PoKPzdl?*1UL$Y8{ZIPM{KmzAkWpg1&aepb}->C{8 zHQItWb}@P4wbDVF%Wy9_U@|W{|nX!$U#_7I1E^cC)<#a_^tjwRYMG z68DQ8mz>|x;A>!dRCt^@`l22md6_LqYVV$S>ep!t-l-Zri@2B49oeDUv>Q2WpUk6} zJa&ZjY9M?%k3#acjvf6O%Y?n1a;%sA2Befn=gr#3qGaG}L5ueh6G(3!rOj=@yx?xH z`)QV_XVFULc&Ba(z5@d;8qMQuqd<4>3asCZ)){_(1C}6O_`*eB0m$}~k8_q?Gz6*O zJqK?(TM#2k-npH@R4BFxTqiQ5L{fZThI~JvLuQ;tU4lal;PBPXM~D^j{OU&9D9Em-z&s`&aW>9Z4aHGGBrHAjqXo=_b z5h3z+dGC%dA?P37xo*06OdW={%wN0o2w-Hs&(JQK3Eo*7FADRc5A2Nm!f=)d=6zkV z%FAPc-J1e|=6*WNrtdFKj>kPzbk+N7E2V+uFx|LpQIX92(knmRX$3*IY-)G68N)6S z>4&$!twEimm!{E-G}ejwr>O@lKyTk1=f6=S5-_WNku>36I9+aZ3g=|314Zh&FCL{9v^ebF~$pl=nyha-(WodJPw7 znqJ)NKl)Hfnc$)X$29hlp9Pj&?V!V^u>Wgh@wdHiEU{#A`q1K zpd|%$e6_mT=C_SusUdN%uUW-$}>se6`gxg$<&xvSt2j_;x!oDf9bU&K4#pyl*mZ_TdL^ z`c84?FzQyGKKYlmkA-@S#Z0NDNt)JM&eI{kEl8vLF?RNLMdDtoe8TaS2Fx5OdwUUe zOF_;@)CS`vFfY>7(sP_a_H*#dEcZqIgyqop3}wt8cbwKUzODnF8NqEWCvCu8#IpS| zFU}hR#ct^k3Dle1{+9or4P{ySAerqtWnw!mIQ6JRKQxR*(AGKjl{D z>E&CIjA-+wGa-nlxk8b3Icf$O7FivB=s%lFjj6g{jyV%s+mCsgW1og4JO~WS&$WwLOOZh^z5o^~)^a#`BL{TSpDZNZiCL!>gD>PzVRZ`pqHZ*#O6(y@rI& z;LPAZhJLx4fNf5}vgF5_vpx^7{^yWnLv(%nQOqu(vTgD*$U774pU+f`GqDCtl zIS(lfGE0hZUtW`XbEkf*JSdEF>Ms2^#@@cAT}9zAe&^-i)zc4srn!5{5;f3c-TKfXt2xU)o`Be&>5 zDZ^G+mAtwt()KGDbp$l&;5A7ymyK74qP2uRgaq&7Vo?L-QPDrp(|*_pUcnyG@|G$BFx@oFegWvPqf4 zzel=NpGR9F<^~sTSa$fL4p42=hc@(L{j!ZOTzoT~IBij)Y@Wva>_YdPMlCu5Irgbele>Amchs~ze1@ZQ%Gs{mR@O!_rjuO?pMsojhQ|hzdw>n3f8-AIaX>9!5?mj zdQ2IRz3=!$D6*Or@;k1~O!f~9uE$b0{p5hV&apZ3e#QXKikV+TTZo<#m! z$lBLN?=?s*#~F_0xz->jNAYz0W=1?-(2m*PV}Y*^t-k5xQt0c-T*eILB{_~Jmjk#= z0P|Bj)^C?3Bk>2e>s_Q1@R6A5AJHHj!m4_apHy8#$g}IYpOXbdZOJ2E z-*M!r7xV2MDzhMFE9%Y+hS8vO?z?f=4;wPQIUvp%{ZX1>{YyW5H6;u84$Nu3L7kaO zj@ZRGJF@-NzRUMEVP0MfM@Od|=EA+KNxL0v54pchbXI<1K--wK3iYEt>{87QH~DNy zDwat8ESf^U-G2?G$D-{>(eo|(s@SL0e-!S@2{M6SBeGEMfOVI+rEh!G3_B-$SY~_* z|Gsi7+@kycKQrk% zfPMrOgYwX4W|*J+g8UqJMjwhe=ls4i8e~rJ)JKgZ+=E3Yt#+!Gh2)TL=brCDZdQlT zvY1PTB=IM|rcev|z&>s9e6tlfNxCr#eNSj4G0l&NB9~HFowvX_07#d0cjznhv->{! zcr^m8iJO}i-`vE1A2<3yDkn<@ge)aOpS)+mW)EX^vE!KUTobxsgP{VU+P&1V?Y1M- zM8~A((fZKCy?FWAZWGeI+WAc8Bon?EY*_Dy`o@O5p_|9#S(xvd_A~i{J{hTBqjmn6 z1$ge_4O!uZ94}9hXit}fUdKT^=xiYK-LPrB3+~%ZoqdivV4uAGV)jEFdBVTzoV24N z?&p2`O}O5xVxH|0+1bN-;JJ3EU9tn#_Y0;rTs8W{^MD11je{j=PEfB-s#Zk(`0!Gr zC5n)vC2?}$ya`z3mt`n#u!WXs@iLDe04>Ak{kgR;?`}pe<@?W7Fm<*;Po;xNN-vHL zuhZctvJIb>9WXfa^lUzGy+|H9xB^rLO6jm)@9(mHT`CzX<2;qhP$YbQ zhr(nV(9ubJ{zufFP;NEaY^~?O zc|kDWD&%)?M?0UVR52jt8BJ-|T~nBtb}M$nIfj9WMd=k{L6ji%^S#x!0BduUf|t*- zh1Gg@Oo`Z%WYtE$*E;B1i8vi4xxoTf&2vesUN$2p^WQf6IvHT@-0YE|<(%M9HyayJ zZAUx~9ekrVMImt=u@P#{f`ns#0CW6}A~ryIX!g*RvsUf7fbOY!U>%Fn_AQ zjWWck73^74w+cj>KmLm6vc_CCTHLdTa(JGOCAA{YA}%Lm9<_q>m-^738yi@fOeYmQb}-9pOZI5_IN6{+sD84fMGJZPxS%WTOZ@-& zCV%qDBO0mP8{^f=Z%f9H=xeN=HX`n}k8BA4~8>pXDeH zNW?67yXQk`Egxnc_aWR-$zhd9{sNibq!Tn39fA*7U^o zH6W(T8Ki*vkC-dWW{X68k6Vc%C+;nHoh7RR(C6RKTHx?o-kuZ~OP>(XHi1m}0Oy`4 zZP=wmONer^Bs;zeM0@h1|AiHKGv5P!zp)q7!&^lmXL(ymG3tt<+Lle-^I?(j(IXQS z2I`9|y;5TzVg68dVEF-@Uxj~7nO=884@&8*?Xf}r)3bKgqoTjACn+*QUI=cay|}F(jJb%PQ;&Rqi*=LqDf^DVdG^RXE1$a8sQW$p zY|c=9oObKs>s~J8n#Ov$Rc5?1CT*n;ub)ZNBZ=R7HaO?XLT+OH)gttf-wcWWNIMKb zQBlhFs4|55r^^M3Vyqyy>P*&f6qTrioxLZlZVeo3@9*_M{(Z_tDLVxP^w0DC{2hS% z-_K8@%9tV6r12B~_JqG8Fd&${^*fganJCM#e)GeCOw*EPZxciE!R%y}VYvZu)o0rK z8>vG0!;d#dRc&EuW-cc!#|DxOBGp{A(SI?!VQc9bPNFD3uQ-77&!Fms!5y+j@R#>! z$u$#u5|jRR+u0#a!acquV|2G7Nm2A_*Qw@)0_8bP_p4MQZx=0fi;247r!!}4Z|jp5 zgOz;kN0i`(#?piPCNv@CHk0!Y>aMkpe|!~r7w5zkDcihW>Vlc$xmyW!3M3{c(eWI^ zh9FSTZ{!>&Wt5wlaqtoO26wpOYz>fTphPI zqJb87b7SMugq6@QQQabc7wf}U=eGRUO9%I!j(b| z=6eI!y~V;_r3v+E*WZE8loGhi3eb1_ML${Vl2=Oa%t^*eLFus`oZDpoa<9lXBBgs< zmWNh|64tlE;<0v1^yCm-qo7uC)eWc$6cBg!PV5x0OoPY1z>F~%N z`t}d=ANz!Jhr`SFi_{oGMdr&5R~S}AO<_uVZzYhwjum}F|1bwpo;FzCjQ8(fwIh@p zYoIxBm!%l$bETh`cdWwuLrhK0>h3-p(AY8hpr?iga$eOvahW)O{P?d~{g)}6I$T#0LR8Os^f$zC=DE6=WdBQ-ExyQFVvs%nc0o&@JYgTb-$ehu|V#%FE0rw zH6U8QCQn-?p?*+)aVb4X5$B;D8WFX&pc0YsuY_$+f>ya(E8z3nzl1JC$2>VF-Ub5) z0SywKy6QqC=0@2X=?WG;w!sdo)MZfgH?=ZxKkJtza%s&tf(PCq;*AZWURN44K49(oYzh7^%@CEi;n}8=Q~_7uU?S2D$be4r|pX>c6p1?57sF9>QF+ zx8EvOSLs9It)^I=4T@m=)_8GH5(qo|egABq9f@0&y(r|%114hAg^Rl=B*5#$_z4|b zQsApY*DVks^(I?PPKBAl`h39}ebfgIe3)OyMZ^5h)XBTLsO#>_NdKwmfq8izn+EBHdCP+{OcL+6Ya8J(A%WEg0|v3~6p8O$m--ZSfNB3d+Qfl#wLzgOB|B>f8f6|% zIgfK0HTH1o2@|mErv|j1K|XGk<~5lGGsvxYaoNI64{kYVJ-8xaOKSX&Y2Mt4_3>k+ zv;W2ELyG_Dc!4h_1bQYp=sGH-{<-a=fLLWvIiiv4hg>P~toxbO9+t%S7JuW$gQ!!e zeh&7J?1*o@s%LHw6@AXP0-D_nFkjFqKTyU9Cn?4Fk!~usJ2&DC&`;ci=+E- zZ!Br7^4iJ=E#3Ljo*Fwtf6L z=r#*RgsYnyq{X1`U#<7MugI%;t3hb{=JpHPaeup>~~J=59&&7ZK_|vg&<_zLoa9u&t^VU9msHTAVN-A1v60J1K0goN~ zD|m@-_V-A(j5VYO#l0A?R7H+b$U53M&H<`;?I!GnKxVF^S`N>B7s}P*P+uycG^}3v z7df%H4~uNqRml@6o1AcqRzN6nN5VdvLPw`KlF19#Ig=HGIC&>Ec6 z;{C(`kBRLYs)d$K;M#8Hz3U4MI9PLo0|e7s_IB+H(}KNN7~)(PjS2PpcRBw)b?Cln@lg}~xEJl_Mb**2)A)h^ z1J@s{muiBq9ei&F;g#{5uNT?EB#G4K=4HXB)N|Vs(ZAKw;1{?uOo?z%rEdP_{6O>U zRWE&FZUQ|wr9=0kUZU?x#kTMH!l1x?I$``aA8hG}JxUS7-_InXf4vZc#NP|MgY{6Y(Y(s~WPGl?ACJfO3K8Oc=z%v| z0TRrP&FDqYh(@GHTf#>=A?aIVRDJ)@RDOr3Uf+UsT=wkc@o6c@xpeyVcM&xrEX2D$ z=%g_rOC3>*D+Gy`w3t76rAeVr!L`L4yuYj?Be(yjL6&~L^yA7P3-cVf3k3OCuprnO z?`CWdjm#&JZpjSj(n&qcTZuZgm+_$V0H6PSzmJ=&JY?DzQF{(aEhJ>JyF+#ala?f6y|A_Gz1$Qh@Ay!S)N`A|K8uyj!sm zIS9OE3!88}Y15p&%Ac5!Aj6)!lk3#UuD=!=)E{e;0*gS(NH21MS1JxR&MA-*=W?IA zMBE2{vud(O9r5|&{|z6~rjb>5!n+eXScJpy*Rl9S0rL9pshf;E>>uKFf88zNBNZG0 zny0n`v~{t!U$aqy1mmX~B5W2Ed9HdQgq$_cn;~-Fwj1NX{p;F>O*nsRc=gi7@c;hp z<2mK&4B+Tdl8{cP!VO=Q0w*P$C-1zoE+Sl+#3^?!P_E08g{)q=qf41^Mf*qVOFtHL zx}IMqERKG*S}(i0JVQvhK6mt663*`>Cru9gP=R8P-q{}w05gI$Sj)&DzoiNdfK;#x z<_LU>dwLV+1q$;39TWnec?^nZke3lG(;_Y8_#O-nlaf}2j@^1^BT>iD-I71=VWdo= z6JNWT@5LOhOFb`}+tDu)lD08g*$QTMNq8n-5QT5XeWx8zN4GV?HHG~H`B|QP{EVZ* zBqeOoS;%xX=+^69`)rK-{P>!toqfnNH#1HP8|5XrNrCCBR_h{HKYcnV-h%KSoD;af zLHj??@Rp=9a;ix3oey+L(n<4Q%_+(tv{ei9d%iFTDc$mFY!ZE}$sIG#DE8!Luk{s= zCgidlo?l*$oTu&$T_ZIy=5QnZq)#7*8sXF4U1jq?1_}*ct-FjmdXa~N=e2v8kSQ&& zGdmdf*L^w{awXLXN79}1DQD)`&mPu2`TQGsWh|bZak>nsv7vbK9sLpLwq;hcaFG?=h<r7eBMjEm)LL_ZIDp@yC8+`IaAJ z7A8dia)QI1ZN>0Sm~clS|TEK@V`Iep4`t%QUo0)`EDDckH>q~uG|pH!WXYU zJOS{|$XitdeO>#NJ^tIyX-9Gc($92UHHOL3%U`O@HAyj}EH6pot@xlzPt z@4kb$$A0SNO6Q-awTTGq`u52L!t+dWvdTC>rCZRaB+LjrD+fHnC`!chtl0HQ%y+Z< z#2ZtJe5l8lE`L~seoyU;y6>Zv0FlR)hQ&>pWO{CAP0(AMOYeC1K-~dzNNkH@c}#Uk z=6k{Z$$v(WY_^U=1^ej~uGqK8=Yfv2>}3rn86;URu&8mRH8H4P+*=%GOKzH{MSh94 zhYy96mVuXMM1M$gYi6h&xoov)^APtkBb}{RxN$#wgWCq5hbG}KKbAdn(*kBqgI=vf zy(NY*8zv+FMQD(IpPOn=NZ$9F$G9qCy)z~{fgA(wj@VwrCE(L zatOE<^me|q0hPlo#|3%mWG40GeVZM0Ql?D3FgI&NDCt&8cl(gj(7+Y4XHE{r22YN6 zxUD317NcX$>9!>6rtRRpZYJT_mXM<+)J$ti%jNdCg!SEA)JNH3E+Td0RH$8s8F_tH zWKl+yMTC4$Uq01?=g^5)WA=+yP~jNh7PyiE5~&4ZJrSbB{#wcsv24^qX~_Qhg*n8& zKhIp{FvtAEi9@O)u^Ny|&C?WlrVk3u^3Q4mEa5os1@4{&79nP>e0JNl$;+saZ5DW!vD$##y6|V;0{V8S z?tb+py(UCf!MbSO4ohORadv83tO6t_UI?wh`Zb0AOn8$Q`uwYnm;E#{$9wAqBYW}> zyLwmZvYS_s$Gk`x()KbTDd~(ccSBy{+#=eziC+f$fCmy0iIz}UTCwC=hB+DAx2b(% zYAIy=A%|@E7(`cXk>YwA^C4I5XQ&hbae2nNbE=+7=GHpDNxi{>IF6=wviFe-SGFyO z5Bs<{ZS5_IKV-<;4u{jbG+88vXYU)!5^b_8Y579@j1{<7oZHb^uYx?0EnjlV(a+l# z`<7P4Ngj)Lb*FhE&)-e>#?$)>L`IUf(6(I_dgtHWpY%gN{@Zg)GJ*jVBtC`&mZM($ z>?MKA7WmwokL2V#@h%?!-ie`|p`F;@E)G5tDWSl4N!*#i z%SI$8b~}6YgFa!IAMkAZ&Vcn3eQzF|Mh+aY`s(*enNYd}7ex2k!+h72jKoqyA~6(t zJU)~H=^KMN79EkRc|qCn(g7BDNVKi(*ophO%|+IgsO!4$J{4TAoBX9hc3){&Ix(#5}o*9IT;1x5vF&L@3gAS+$)yi3=^E z9&ONvGTlVeZyHQ6I{Co+dpG*?CkEpVifF>Ex;vdJZy6-TJ9PS4+DcOM4O>PrWukF+ z{@lf9`2D3Xm7`?HV}1-bKQl~}a5$|!;qaOtSi%`j*-vopQP=yLeMlBsjFT%b)yWZ& zBZ;yWVWv=-qP0w2$dL5ye)i}y&JXU6f8J=C1JKmZJ9uye=fnT?2p+GNA*^VD?-Li0 zdkvj>G7e0#(m*_gTMzqXSJ^PzOmlGBsL?coGU6G*>GMW%6i`_5N<&vh0{h8pN3PwF$SET^LxZV>U{CYic+)0o*nBOK?pIp2LUPV?yT{5r=cZ$^Yd3x1cs zQ$tQCU+v(mIqpfRTw_ssr2a=%v@IL&`S%f_^Jf2OW}RunzZ$hs|GqQ31p5UFN9vJv^uN&;k~@ zE*CC(0r8N&uD^dhm84X}+sN=yz@zX(bM-|-Qg=HgSzP#a;0k z8~>*oJ9=>*FcR%verykN9pU)QrdS&K8GJtS-#|UT{n`&Z2Q-K=-w&SI_sDBL>B{Yk zzJ}^+Pl}zq1`If3>vHK7Ynyz-E+Hvseu24=nKD0?YYU|SAXUfcRG=M!`E-- zYe6QhCuW7h=){1s=4Wg;oqWI1bdawY^@+zXUM)OEC(Y*JFH4XsrC;|{LevF$>;2Ev zzJfiZ`wpZHk0T!}-rMY!nmyF4dD#|=^Ky|_x_vnZjlir`Z>ej&6u6H%8uskBgNn$# z9JQ_6= z|9I;ZZMp;wmUQ$bs@RdG;a%pM*NsWsraj@cil1m4r#_mwR86o;U-n7DE=iDEmhE-g zA2}0i`}E@MP>1t*?&$@rKU1n!WIsy;lE0gnXCF5u5|jH14QQC3?exw>s|WjAZWj~z zJ>sNJu_AL5K3`w9%K?$R3HIC-!5$up5=^XLxwXU!^@Mj)zCYT;AaQLgzU;3XU`qt0 z78awwFVeM*dSu)HIJ+y4({O+AP*%3WFE|bzEBnKq%Ua*;x`_gPg8U+uZo&}f z*dUV1=w{0#8c|GMVTz8}HvP~pCd9;ejw<~5DRYOs<#A5-E>hjH z#~L)h{Z8XQX-+sl+Afn3G$z;d61119qJMeuqP{@)|Gzs`@e3n{&~^F2=DYr=dj+zO zS!@V%|1I^_-;8^%#vPJ}?O~doy6$khEksWx z-ahh}MsBGdwKA2(@8HBR?F$3X{d_+EaLngxxZ1%pLubGP>3~%i29Oiu{WbQ#VSeK5 z<}#SvE=me0m-g@9ZU7&jPeopNVn%@Hibp=~bqfpy0}Fa+AmpRIddFV;IseWL^&aFS zT|u5Luh*K9$HrkO!a$ydx4qDPXPlFkBz`IVi1~0IQ<`sPsURQDRWVgq43fSqCY@CR zFlryZc0WoL1f*5Mq#98-yM4n~?R7Y>rJeiNs%lRbAMRMW-;_$sIE$xOQjyzqM0Bhy z){=B64HUOuLS3n#{=gXT5^%F^tSHJgBr)|b4^4*4klGw>=5E}xjWM*pN37C;0o9X_ zZC#PScXkV>ANtXMn7{w%ET9Ejdsvx&l8}%3H}RhG2CM_5c>=TCR7s3cnVzegJS1c& z7BG>k-nU}c4$cW7d~aqO5{Ri z%kAdeW(!B#4p;IeAupeQbo*xfebP@B&iuu#Wk!KT=xNb|;q%M5Bqrv9Rz&+#5OZ@B{?3c~BZw4gV zlMPqfP%et+N@B+8zwh~p*6V!5d}S7Z;@&If|D?zni$QCLA`1c{I|Fy57y-!|Tv*AN zp>4W(V_*#DJ~>WQr-L@um>c`jP;C?LKaNh*hLB(O^H%Dz)4yoM?$(d(rRXmR<*9gN zhrXUMjz`Q5OO}wt=y(a^A$3Rzxmk3eo(g4EHZRmx=#v}eKRaul;B!mq!UU-3zms$Qf_u6a>9gl8Be7}Sntc3xD+R*$HH03%i1}mXqvy4Y%%HjS(vkQ% zYvQWMXPy;|`>msqJ7cx@h?&f@_p8RyKcK3=$p5Q_-SM}X8DH=Y_?Bk^c1b*LoI`I3iDtO%xqPYQzc2F zIeXmyp?+I3PbRAxIk-+~l^NPL5Z;lmR-G#JShKc`jpRPg*Xydp|jC0IY3q_0}~T zf#Q?#TaTa9N%$YXWnP#|FqeN^;dGud;T9hqxPKn!h@UU(R<1Q6jf#hkS#(&z2IFJn zp5An_8|W#oZ=0ikWAOHBU2$?c>U=IYX_#xX(5jbiO!&V< z-TKgF3<9Fwt2`drfzVo4p-tAv@%!A~d#1z?gxvWBG0zNQY$Q(}RaAmnj-+ST0#WbY z_17?#&zeX*-(kVE5#Wd1`~K_}^h1S;x9+dgB;lLY*Zk4}R6+B$TAk3h7YH zUCFpinGj$gaIr_B`3`chiitye6bsaw4#@l3;=J{B`Ba01F5!sYRj#}9HM=qm3wuDBQ4HVID|w@)NsF z3N%7&WJaV9t%N>h@1Uc(3v3Q&QBCdkQCdq4xwVqm!fdKOQ*aRLNUmpB_ni2?J;`n`Wx%`w@lZTXX|09*x6~W2r#ZH%H$f*xYE)d$oPt-yJenmPl zA-l;fE@&Keiyq4RIuvz@r}q`lf$zvQFAyob*Q*DUTaPa`1ga6Cz$KA}2Pvd>&xOFV z7E0uV3LNjGGLRc7p;fK*jUB$}d~nBCoX0w-i}Jc)&W>*D-MUFVGQ*DRG;7m`zPT3R zGs{UTP@U~y^}65hzwuiW1)Zs1wFhKDqnlD1f;l=_Qm;~1<2kkc&J@G>sU?Ih z-I&F8vB$YHYd6$0rMwm;SKXPxCGEu1Pi6R}*nPG)d6kabNd$Uz2D-eu4>nwgRC%s(oBn+(ZzD@g@~^>}~9 zy!X*E=O*3LDc8O6{7=cy@G1xpC%zo*E_ZE_2S_QC$h@UQ8qF^JbC)tEE7~^EOYhi{ z#raErh5sS1CA;nMBF}1g_4M)YcXsG^t=cV

s6ccMs`xb}X{vPwitTanwKWA*5sn zav8D==8ub5K^*hS^55z|X+boD&fzp6$SYT*aH-?D7stLivyusCigg{H*xHks=q1Cv z{w(5o^?QF=9R=o^F1ypWngI2-rmYZ{3Uuf1*&2#^$scm-@16OlKK1&R8*s{H?-R+Jd^|T;sI@!uY(J)J`0DY6m8a zD($RVbFzoO`A)619h9H4WWT_Dea%(xcB-^4+z6lZiW_0z;LjyWTxpOUt3khFlwgE; z<5%dr^-PKCmZM6hxbKE1-Ond-1aSWEVzVq(ANVwRgln$hy)iZypM-v240tt}S;K;X zR6m(L{uHw4soKMvZa^~4{QkZiGy{fvLeXw+eAzirN7&@e- zH-kFl>tn8;O<6#gx21>Q;CxB0bX##a&gVYr*Pq~&A-}WQdv;MVcW}0c$B0)I3i#&j zeY{s1xr7zD5nY{a$7f1r-rRj=VP}pyfrvd@Fz1ZUwTA8Zn-0I?*4q>-8vAH_Rz$< z>$?|q1EdTvFLAT`fF*LKQ;#W^M45ngL=SHU?&Wkng!I~fVgApg4d;@>LS&bjy-zCg z5+d_NZT|IIz|DBCGQ|y8|9MnqcRQLu-|mE~XF6?QB5D826IPhRR;>7_aKIFnDvK>} zJO+};^FercGv@fTY`ft@@O?Jv1&L71`MY)PQO+(EDCORAJW3Dw^#u`kr&j~1cpP{u zz!rH1H+D)!&MT9&fxU%!@pz8%$7zeXqMj+Z!;);MDe!`efvqG9?d0und6D zK`QEV8;gt$7Ur*+K_3{W$6>MU5~Nb_lCSC(e6Q4hq&(_0B^=4-ZO@xV*+xdGoQsaC za925_;qX<=u@+ZP^?rkTYu6=$@gFsyG}3d0`ddxn{y14NS`BkcdxlOXm7>{E&r5U# zj}{o|OSUqt)xpfKpm_}Uk@U9wqJxt*V71iAY6$DiOTOtc<2L3bU0}S^rOz0y$kz`Y zAHm;0qU$pS^(gE23uV8lu|-}STWQtZf9#y5-EV`xS-=~!88ruAN@-~0|)K+Rs7Ek z$z0Er;g_Y@fBRVm^jKhh)zoPY3o>9Vw%11r*-_40coq0KBD=X z1wyaxmnVraF&83HZQmnxVlr6&V*Lolu07_=jtkNxJ-p_B&V92YPC_om3zl}oziq3- zT3y_ue+LAi4&zJa%quv{T{Zs2VU_@Kmv8&TpBtTy3zk|y1}JLas?r)l+HNdeAp~XXk$>N7HL~` zelzb71){GdRD1L@VL{ZU%Ewrrn7mWB*Sp1lTw{S@olD{n>9O4A#58if^@shV(Kqqc zpk(cCDIj%+c31wIw;&_@%cQRSMBRSPgw5KOie#wz>haII_K=fwxk_<^DKxK|n6%-? zoIQC5A)8zP4z7-&ci!J$A@|X%t$!>j&}vIB4>dF=LVJH6 zblJs8?z!$^K7BXA$e_sRJvi<1qB2j<3tD!g7LlOqE9nU_D(jety``7wfG02&oi~H~WVUpaLu5@dC zK^S_J67c1wE%Do$z`eZ?^`}!S7+tffBy-K#Ht84>cwZ?g>A%P#l$xJ=Hmt61uF|*+!Rg6dw%oReq|}tC9W90y~IcYoVjDOf7|hsmM108%jax~ z;Q9si@p(JK_o8nHGe?9>gf)~WpJR}h7Md3qhUjE&O!RfNrxje=V`36|+MJvl^8eb- zV@#@;GPN-*TgVg$Zacn{1ycGP9G?UM#!no-5QO{Nm6BEa*F@oY8f&>Huuqf(-4iiL zZZv~}#M}JuEOEb8D_bI#hx?U_df(PC>A-cW;zs(M0nu9(70b*;zxZ-lIXf2W7urbF}Wm+u9lXHOYe~u2TNvROBia-{nWY67}KYBM7HL%%H*j5@lN;0)pbHJCwml z;g4gWwjs21TGc*8eq2p(L@2XLnTRaR(vvW!t8r~@f{E{HIH-11uB6KZI&4$LvbiWw zx}JSk@w6z>W_w*3Qb0YQea4y~C9Ln4#EwP2rNZ&czpvfu#yyRPAL}>rp<#Y!>&q+h zFddtCGWLfJbn-sPC~Cmxm}WAWUC3rn{4s7mltL%mM)E7@Xmbni+-srgAWU2zCg<>C zeq&!?x$w2`E8yUnlbt!pTjTgAZkmt2vbpJ7j*Xiv$cfq9f(_?v;c39>s4E&&!t>H2 z_z|}qkv*n;p-zKIG9M4EYso=f(AJ9s#jkDAMzYk(kYY*{QZM<@2bUH_HmN`u$3ilM@B#v~?7VqNAZOZajBw`o zla5D4Zr46BNyk4WnocO{|4yr~q--&P3B`cXWb}8epXi`^t+ya9KgxNRZ7?8wp0m!@ z4iqRn@nCOBB#ktZdAU3LM2POhfq9pH6Y>);&u{48CW*dxDkhe}kI}~2%K+33$y?B005lm<+r%Lhbo&tJT8vZEmx1zjEAgr8d5lOM-q zoaeTxk(!oDv&Z7v5VUt(B*4Z3gbw_0=nm&1W6|q$|zbu^-_sU!EtZIvTZb+`~eQ>M#A%nDeDaJ>fGy-Z=aBCRT zmblC4Y?5TCgJjyjTEi-w^Xu=Eu3T+~`-a;ATD&}P(`!koVguFLv2m*`rn}} zP*ASgb*uvQp!x$V0y;$?H;P(NVSs+7?dR@gZlJ*CRQK8-^cj@We~LU^W{7$&k>zK+ zlt}ZK!PoM`I;8ykfI?Coa>V%OHIwd9;Lg23r@=oga=C-6Xe~bjxtiM>%wNev?d9%^ zCX=dA(=U?0Ax4Je{5hFg`pp#PUIZss@ECxAchH}6#&oF4qFvtizlzR0Am+D=;!{eg zv~OCcecGpeuiMO2Lq!`R5=sdPNeB@s6{#deLO-o4N|I_G38@g3LMcfS60Nj&@B5d( zFg4$q=RW72^Es%ygm39A-O{2-XrEnm0{R_6w=82}{3VmTwglG0CKt#{qy2iKY!8*h zFlrz~8`QVer2T4OL+wD|m=|)%m6>*Tc|Nd7$e)iQZ?i3k%lN;9g~z!6lznZTB!xc9 zXv1YM=z|X7y{jyO&!KMe%Aiysoe8bQe|1<-k=yv^nli`Cfy8)yR9>;=+3_NgJTpFDf@lk&OOMFpO3; zCkm$g_D4)DLH)@w73Duv!gJS0u)xI(!iT83Z>AK<+{?RyF?8gmS1xxv{0Wc;mtwMr z-yGkMw(zn@)O>q?6m4y@C+C7qcY2yx63X8*Z`OVJ#W^4Cu~qmY&JPFWPaXfPOJdw5 zc}jxxNN()d`tf5^obWmO(;9SD;%Q#b`L1D&dEMckGf_??#d)b!wE!DlbH!!I3<9Lb zHT<{Jkp<<$Jv|$|mXOAVQ(Z;KQQj+kAx+VoMPknFcD{I6j^rF#dHV8VOG58c6V$}q zvFN**AEyEBGx}_Y{lN%3HHji zFXE(l?v>71f)#lXJ05GrV?nUL(54d47ul*_Z7#gu^PF!!)Vb0aeY8<6ef|y*e?nN0 z^e6zS*&DQdiyf4E@z^bn6osbDob_oJkrVXt3T=qZfFzB~8zwIuNbP%J{uc|#%kYa7 z7UFXN&e`7@TkP$K=q2WSy^S1%?um#Tv}ePp!~xmgO-{t;xN=X{f;P!r_hr(zX$j2n z{##tggFS9)`^Ddst0Vt4KXxWG|K6VJ#asZfi`Q(O zE;S`7`m-{W5^kC&*N!JI@;#lJ(Z)rcIYCqGgE{WEzvF{7^95Ywfk<)?{Uz3WK>0?3^`A*ub4 zqsb)8tNPlF9=MRP_D$Pfa?6pd*6w54+W~vaSC`D5TShd`f2zQmZt|W=|5W<}^Pc__ z9Pe!cWN=NU>o~qn_O1>ZIbjT~a+?Jo-eo{4`4YHz>7^psN2-*6$tacAdQm9V51IxmVaR-^2n|TrKDqxaS0{ zvZ=K;MOxrHvguduX(O=Nn{>Lf0DVL^9~mXA(E|F6$RWko22f+9tsRUy!BcI$ec?$= z^4)&$*x^!jpvF~=_vnb*KIjIoRZ+5hv z?gRny8=Y@!A~#L)p5dUC6wvNi`vi^9;lO!~r;l%%5V@>!F%KICSrpT2z8JYh*X?_M zJ|1Ji4i}}{h1ic7{ogWWp&mLZYAv zmk`9z(X}KaJS#ODmF&pM+4DXt8}PnDy{07_V+Y^A?vt$d(}YG?^<;Mfm9To#vqe|X zut)q|B5eWlb2q4An`ZEwsD3YbPI?(pkc=C%t+$7$O;6-4OvOQVQ>1#Wqcu6Y>BpqB z7wXv}lX`1ql}Xc_#K#}8wuFC$x8!+qLz3P(vil)c3%zP z&7d1X%h)X+?hVNOHl=wEZI>rfZgp|DKmdVl@^gD>1D}txSJ{*O)EXzRnE3bO-yTY}c2S ztXo2UckCQCmt>KNq*v003OKLW^(qKl$herQa2RdTmqlN4K5H$@bsf=ad21yl1&7$p4@mm)p9@VF}qX zv_f75_Xn)i9yZE0YUJ=XA)7$_eCP#kzkhiX72IY+S$B2Rpn9|=Gnpb!&S)Mxcx(=H zhJ5SW{LfOM|4IFU0qpfES+hy-cen{ya>{m39i;)!KKf5 zpg$^j&?CXwnpEE$zx<0|1={XWUfaKS0BTQwU!sa3nNv|0t=WM4#Br8^)4&2J_IJiPqpMNhv&=9SRxJN_Z3Z~u|C-k<~By0h;J2TVwot6p-t zHFC#eRfggs=>+b5W}0kdk{tqPTvflAkcGgL|3!7NNMc`6Ml+tL0eOScV;a~uzWc53 z7b$0;{MBCUzxEAh=&j&WZ8bVkpAT_u%F%%m!Ti%>1-j(y4)whAvnbn!eO)2{341oq zP3Rm_!abI-Qhn>JCKT+v*|_e83;1SeX{r4+Af+DZ{1`DYBjvYlr1TtSlbBMiQz}}hn;sjh+)028Tzw<(wE%g8fgc@1P-loi9-P=? zL!8$>we57ZfbdQH{qceOm$dIukM$JJ)+bL!C60UL4yJqem*A2fQuDi=t-p&a~!rOhV^Y zt2%N<0YW-l&6`Ix$)HuJ!t!+*gqEW_w-WvJp)Z)cQU>_C#y*5zc{ImSIxQ!s{6~WKL;NObf>#=VnS9SmEuFqy5Yfr0l9Hf$o2YbiGUKzt{m%g93 zaZg>++-o%}Z2=!NgpQjw;rsdbqiE^_d)V;#Y)4R|F%0LLDBIeif8GDktFU%!ka?H* zX|E^xIv;MOtf7N#<%sz7@l-T3*-7fY}UhVl+pn4IMTg~^l8Kf z=Snqbmpem1Sr{eO-UV#OM=f+3R3Lyl|L2EGJQjmRV?ec8RNvd7lYow3H8qAgC~KVDtI zI2W=lNW1>Z3ny|?x7bvG26)Wmh)l@4v@C{)gNfs$B8zjB^r%w#u$2Kf{cGE0g}_74C}yL8d@&Qk3)NO@cjTYaf&l5o{rMj? zfXBuceBkkyfk{i%|*A@T%>t2mxIOk<<3$8BVQ33^*=nLDtbcm1S#}TuS zW~6e(kI~atXk^X(4^7^vH&q+(eQCVmLi*o*>1me5yyW7>Hq9i;bCu`;1w7d`FGHJj9HA zNNykDgP6L|Q$3dyNzunW({9**>EWFjW)``ayyDkZn98srmIl|qJq_U1owwcX5`VHijyX&%QOG@)5NqK ziP(T7A0V%A%SOF7p%K_4UvXJPaKIACwW>pwc5@ExFQiUBIbuw}B69EA3PX~4kUi>r zObSv8^sJ+<2$QZu@mX%j^PM{(@#of%Bb>kOkd^zufwUjktaRARk_>o9tdm3k7W>dP zE#YW!lF*p>zChiAe4prjdejiPT6~lDD16ANTkU;%V%iY;<@0t9QRg|#Q`FN=YxuQ?()E)>#;QlCly;pFhGh~(28>=G6 zf_bO--jjCJ{an9=j6A^p!&RH9fcKd1hTr08&W0qZCc8;1!5MNN|C{+rlLP+Yq77f| z4MxbPnveKGfBVyEQQV4{GxaQ!vp{giMqWWA8@n*K#2 zC|iwGn;2II6i8yPS=kR&P3(`lU;2{Aod&tw8B9KJZgA7Cez2B?JnwBxx0c)LB)rT| zk1Z)oc=DHDHh7}}QAbWbzaOaz9iI&n-ILi6QnBo6KKjepu zjLuv6N04K@-aY@3lOu^+McCni&5TA7 zX%`4sdq(i11r^*M4)hB8qkmE0*8XWtQwX>6*In8#hQ0MyA4oXsgLHHqx2zkToODnq z$<{_5!k-eaL}hb`wqGsE&Y=O_AZt1y!5Y%rcYDl)qi^!j!u^xj1CZ0wz|33BAi4MU zUi$oi4!^$$rww8b!6!91OwbtpHNP28v&NA-sQ%@P@DU-BX=vnH?fbkOmXLE5cUj3{%Jcd$O3Al@mH>2MkJbSRntFf2OY>+PV=V`=5&np;}`h< z{TS6g{T;FCzA6I7A~wW}UbD&u=i`3mV$Jg9`2GYPC1WKN05hRme~m{u9dcgxX40u7 zW)JZuhYz@nszIHHb?$)xLT; z?WYrr&i=ch7KA?OhiWJGPa2SF={TzuMOs9UP-vt9ImGd+RmdD*NfRL~) z=T0IItL+}PDI^*|lhZm%xd98d*)6Wvp^e-$XQzuFyD_)4UNeMl>;f|%f_-g9)FJSx zUi9Zl{CzL`k6VhuO9Fal-zJ}T1WU6o^q=|Y*ZyqRvHGttNm1;wSB{Z}czS-&j(ymd z{L+&vmr4b4c*n0hny62-3w-8S8AA-E(>){424X%RPub?91@R)=%~s`_1I?;nn?Z*Z zX=*h*xJjX8unkfsiw!Na%JJ@$@8%BS*5&JgX zX`v=GM~L4z)**RW26HCicRTWo$@u!xR_hY%qp)c>=vz#~`Jk@Hb5aB(cQwz_S82jv zs+D`Wz8R>`b0w28dC@2;h*dh<}U~0Q?XRmUKSZH64wqsVg{BE51vrMxwPM(M{)NwJ{NwYW;1=F zBt`nVKfR1W>YuyGZO2}Wl#%Ka59(~e<7i2jT(AZBs;SnxOly+7OhMCE24-Mu1f4Zg z$fAmNdQ9(T5a(?M-9^Wx_cJMVAU*ks?##b8UY3$`K3ryoi?f9{2bChuX&3 zXD2vM?@DzR^EeTy+P56L9t+YoR@e4OoR{R8R4rdFYXc_^TApv-YfOaZ56MsFFo3S) zu-eH6`H2_5B=;=hB9yROp)%N)GiBW{bI8h>9L+eh^@G0!+29%!nDP)nw!?Lm{0;>O zc;!g7n6?8C(7#rki(f1|KS(ztN3o#>oaxbf`r3PkK`F z3uzeB|4}_XunZV2=c+q&P`9C;-C@2~k3_SnlZoQCpmS^TlfzR6QD5_EVdYvHG3t8o zHhBnh5HgLCvgn5-d+y3}*I=L7zyq40t2w!I`a0M4?JUA$@2B##4d?Te71H0+S~wG` z3EFF)8Ipi&yKH%pUoDrinVm6!`{j!9O>^N?BL3@`j?g+FYH!2p>$pv^kLg)hMWhKV zvP>}JPO>Nc)tP5CPvRcHu)wNBrS8$KhbnFbrx!#m#^Eupq1e;iYGvEQ+e=n-^ zVsybF$Fuk2DfAWkg`U~?Tn|)_9~Lzx#mn{hPSr%iJ@+}b8oVx9QwRoHhc)dJ%GX?!*xbX zb?jFSFaP3UK%%)l)CP=@(;yb@8T{HD{&ID%P=Oji{Tj*buEX z<>2vm<-4UtmN0y0*gnP|{fEL4^8fID+2(C7Y8os`UVO;gKNg8|j;CCozJM}$k<-X8 z=e7AocA&8+GFGhSS`lA_z@CS!^RGI`E5TwO^K#f{28w zS-B}D(Bd1rZ~eu4fl6BD+;K*#jKEBU5wC3Uc=H zt!7OC^fTvg)6XE^=G~cV37b)eQUTWlxN* zN?xE=V(ww%&zPAaG0?sl;w1bH=cFc)SKksbPbEfqV(|+5$7RgUEJ?!WVbw3$;4EFz z?mZjU!Havnk)?lD8=_yeQ{uwMb!<4pt1;3k!y*bRN=2TpQU^EdnYbJwZOkRb?NY76 zx&7$ouTmS8A>?i3i0PaKv^`uE<(bb-_H)sc)K{~~%m$b07ekJuE%!{}?M8JHaB=%m zb&ef)N#F9UmShs?#9o1OPfQ4Xo7bKIV?IKA-^b2SP$NxZ+$~2wX~AAW*==JbPC(oI zcl0Df7NoSZjHmGXX~ww0H3)mqCy?al>_mn5)~L>RYh)1MAA0dk;s7UVqkptUqykCX z&&PWjeW<=>H;ej`En)BU9hrK0fPm~PwpNpBptz#8*BSK&q1EM&|Ls#J>6#P2I>W7j zwWi%PVYMX)3S@-6vbG>TA08E|;5i(~xZP1xZAA`KMNEG38j~#_uJ+m;VM3-5)9-Y- z3;8?tt)jigtOU?W;O(|Nm4*@SjB=P3A}O0Uqx0Od}}~+>tk66NuSqQ zxGsme<%-IwrSgRQ>pPn7EejLl*PE(U)X2f%k&%u~Y;r34MX28}&ZSgdMrNrksp%;C z@54?uv=p&a4R^4KcdXd6S_NY=cd{Tp`Wl;1G$fzz{qfjg{Mo>9UuX2GJzKPLG`^ZSyk zcCf{I+li&$8PI)EWtXgsBrqc4s6h|WcU;7x*JPvawc$W(*IpM$y0qoOt{4^aLqtd9 z=}HsQeIUQoajzo?)jvz+R>nQ`-(MNud})Lt@Y6@N|2-#VBt?10oD9hnJ00Y`j{5S5E_}+j;Qb0M` zc!&xC6yK}|=$rPxox&SCVMcPN($7`tD3L*_LY;m0bO>wIpg+9M5h$DYDAa{Hkn`&u z%N39d^+va6c^Kwdg|0+&Y0qPCY^#)~rIiVIgmU!M4|2onzlw=DEw~>QUy+|Q%1;7j z@9&siq6Oiv&ipuMi<}ORDi=O}W2olScxIyn{k&4fYChb%~kT?#H{(|F(2-YDd?))x1S+wS_`aU?t4 z`QOTlIglRrf;knAGw>Ag4mSD#2+Qr7x7&a`1%*>d4`!T*=!KE+sdhE;uUX>9DOvQn zw)-=@)GS~$!$Vjz4$tYxGNX5;R1$qPcgr0aI_W$1@ag0GCIq%=oHFs&!(6y*u)%&O z60j=qp4EjhPE>2*`65MiSZvHSbBf;)tc9u@mOf)a?OP@IzArQVMO|}_E~x#MnAsS zepUlGz{Ho~`IWvVK;ccFS6GeQb&)_7mZ~!8cc*)Ql$_ywj%(dAzZJ+KM)#ioWH3iD zHu+qasFRbgn?FBwu_P%WCBB-d|M<-7fBL=O34ZWyxtA9QBq8>PJnv0Y62Y(KTNr8% zDIN~`3_m{N!*wK2{|Fs=qJp(5B9Pb3v_6s8$OXxRO0BcI&~Lwl@7eteCS>(gteGI& zicqgQzzRz~)N0h}-`=k-*|ES0)GaZ{T!}rePs$IK zG++++Y?QT5g(<8&%2Q+#hjR_@$H<+(%|K!1setcO6)3nAFJZ8VN>sNy>3w(hYarf;%7`D@aM`o8p`vUnC*wC!Nkb>t0g zc=6uqWsD<4t+_CD|BMpm!!OTee8hh6y;Ie#sJA7B3VOf#5BV6n!bS9HMM&X{tPu*6 zCndF&eG-@xPLF+CoFD5#u3DzA^D{IhO#<5L{%0+LYMENz@ZXd}O5CpK{VQ8rM zt7@2&sG?&Bqz#~MmcBBrG!5iU?X(z7V2hk<;# zzR2a%OwyF0BNlkk9Bwa}{Wn;SzN$~UGY2lxAgVU{s*`h_-^`%o;I(|FzZ&MaC}%8>f9C{a=uKTve(wAK9+t<3GIfq5zDp8@8&ET*=KrNe4urbx6Dr{_h8mS+8?R)ChYmAFb}aCQo{?J3r^I*+-vIKRBmXNYJRpR0K@9SqtW0+Ym_WI}PNMgQg-2U|qeac-vo9`mWYGl`U z*yr)c4QHrV`S6;1x>!bO>Wwq8Jo}~N#<&t0kfgs?Qp7nyBJD?i)PO@Q zchZ*aZ|w>2Z{BOv!+?$@2W*3zaXy;rDxbpp6P%-59oqPlleEr5>2oLEV@@W`tFO}r zM%Cp>(-m>QtL5r$J6i{@m~8Gm;U(U z%%@xfB6L7!qozFWA)hqcyY*t8A+X}m{=>*qZr{b?t26*pj|~NV=pSL17XOnnmLx=i z_B-&14RT7ql^oS!k&e1!e&)z`i8-f{Fpc`%@t@u$nT4jL;Jvctfsa(OkiR-Tz|N8M z+h3U!KwiN_)p8@coiULepfqtcp+76p^q6HC=GV7=-}`(c@+|rvbMthdZjf}oEQUdc z7`47Fe-{mM1a)!-O%oVo=TfElh#}-LH`=Cp&l|%?dpvj6QznV2Zxs1@Mw(m+%l!9S z(u@dSUvBhK$_d`Sn`ZsQ^DUxY_;*JIix4-3^D&-QFkC5dgDZ>)6WJfC$8qkhU8bE; za0&fvu0dm3Y1$-yMPO=GvlGN@zuC^-_Rm4KP%7rfQCVOwUXjT(FhHHO$jDut1^iNC zJRj3g@2GQ=k?XR8vATnY&WA0=d1LnRQZXm;a4PPNt(7^cQmAdBB)Sms!k{0EGLYMU zF_2mc*rzMCC~}dQ6lAIHs9cM^rI`|in)){xIH%hs=dGX;x%{!Gin|=h)+}0=sUDL| zjE!vYxord0ZLR|oHjdE0_(Z^y5jhxqIruqJl!hG4sy{BD%;4mXw_lDo>f)_|zvQR{ z_LPdpPu1Hi5TRtLS1#x~=jX5T4>7|*h3n6?^+G~$*7oZ}qBiRIH9y}fvYFueF^|{& z1o{t|NPnhuUcfJ2x%#2AD$^X|BZ^M6m@a7cZb z6I9O^n8id&BllHzX{4JGwDEbbF13~-hRf&eALBVPx-q6JJrsNEUn+k6hWQqo9K~A} zk4?#rvw|`L{5IrSq4m#|3&<~WU;Wa;6}hqd@6MfBWCoQkjJ?@AbV$J8nsh%DD)3Ee z;>YJMZ?TDLcZ)4Cme>Z7- z&iJtQC(gM;UfD{5yd=}=Rg2FrYs3gVq6_F+5UJzQ{8Q&FNX&{U&wF74m|Oh%hPI`R zBPx3|(^Je0=vOj^3nxWLLC3#}s!>%a2oYYVdD)q0n`JC_%(f)rmzME-&%pDe@n`Gv zEy(YY*#G?iFXs8BrPqI`WD&1ro)#U`>^7Q$r|<9ao@fUf1z{3pT;rl8|M zzdqg@JJ~7wu-Bn+d?!0fSD!?^SoiQdFY<7x^2UcRSirZU$(rEL=!>v>vSa%40*441 z?LA91ge1?>&HAb~5N-abwG;mxa8G_G=W2`dk>>}syAnjwwW~YvHgYDf99#VM4epO* ze@itAXp^ADJp@ z-0_~>qI5m+H|9ikRDJ7OBu-*pycJbAJk6;e+*#GHg!#N0fU@Z_2mn##9pHy zxgDk6bkdsFV-d0)`x#D9w+?-`1C`QgpN{i(B>JF-AGaFrC60o^Fh3BzVYm9I>*|nI9&><# zeS)9-@%j0FBYx+H9K4?<&qrQEf40XxD0xY^z+0fp_V_gx$+fuHzAPdC7?SU^ z*{!$xLGnnRv8>u-6Bi*aPgzVuJk`y;DXxdQDXTtT6_>W^^7dTZs%{=lC$o8$j??f8 zKh`kmD}U#@J-7pZzyS_$fCC)h00%h00S<70103K02ROh14sd`29N+*4IKTl8aDW3G z-~b0WzyS_$fCC)h00%h00S<70103K02ROh14sd`29N+*4IKTl8aDW3G-~b0WzyS_$ zfCC)h00%h00S<70103K02ROh14sd`29N+*4IKTl8aDW3G-~b0WzyS_$fCC)h00%h0 z0S<70103K02ROh14sd`29N+*4IKTl8aDW3G-~b0WzyS_$fCC)h00%h00S<70103K0 z2ROh14sd`29N+*4IKTl8aDW3G-~b0WzyS_$fCC)h00%h00S<70103K02ROh14sd`2 z9N+*4IKTl8aDW3G-~b0WzyS_$fCC)h00%h00S<70103K02ROh14sd`29N+*4IKTl8 zaDW3G-~b0WzyS_$fCC)h00%h00S<70103K02ROh14sd`29N+*4IKTl8aDW3G-~b0W zzyS_$fCC)h00%h00S<70103K02ROh14sd`29N+*4IKTl8aDW3G-~b0WzyS_$fCC)h n00%h00S<70103K02ROh14sd`29N+*4IKTl8aDW3G_>Ti$xg&Lj literal 0 HcmV?d00001 diff --git a/pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D4/m_true.npy b/pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D4/m_true.npy new file mode 100644 index 0000000000000000000000000000000000000000..29a954babdd544e2370c164cbb2bce9e5df51d19 GIT binary patch literal 83072 zcmeIyp$fua5C`yCeTuJr3JZ!LLkuRHVlkLxxDA6i*~B8gfz20odj*T-cjdT$Jnnv< z?S8X643bCkj78a=nz#sYS!FQ|v8tPUlV6W@Q + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D8/dm.npy b/pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D8/dm.npy new file mode 100644 index 0000000000000000000000000000000000000000..d374ba6470b755704f8d148a715c854c697ae7b5 GIT binary patch literal 13952 zcmbVT`8$-~`yabVi0sSQ_kA$dxo4lGQY0cFN+qdODtib~Dnu!ikSIk`sksYTqLd|- zN|Lm^qf!y^eLnxfcYd3>uDPx`=Q-zot#fk4n&lq80;Ij9g9ag?yEpDKa3mWzZ=xEQ zk_|TP*cG)aC}P8oT_K_W`xBU#rUVyog z7JA@+ayrXmvkR!*oZ1(&*BCV(Ev;9x zquaLFK-y}H(!efj7?c_)d>To|A>ZJ1g?%%8ODFEZ1q!Oz%GxO>U~hxuS3OVeDQH08 zs*8_;ef5y^sOOY@xhteHv?BTbOz5;Q$f(bvLX_+68r`pC^fHb}2zYA)FEtlxIRXz# z7ktsH5oN)`0Pgg^S{{--)Ze_87SGr2h42MGT*Rn{cT@ z3mpWRx2yP0Sk&Ndm+NPPDIlxeCd0rT#ZSJn85H#4mY*;XxwNEC`?eM0q4= zj9Z&;lIDuIkSE~0`mrh(diR-K5T&~UJo(Oly=MXXuFk#l&4Z3Eqdx;{gY40Kjl++@ zdP`iX9Xg714i4=okU&owTDlg}Z`)Ynif?rtnLnMes(iYvXEznYY;Nj!{#3!@yp^8W z>GlBAo@yWTb-_K8f39lK6o<5%#-%N%`8m3W+B^!_FfSg;_<59$HFInBOzWGX=}yV) z_a#PP{dyuKHGvD%I^2SX@0^u1GiZ#@Sh&TQUZ;|S7Kt=ts(7Tr^8))U6A=Hnk`Qk08^^`iEKk_D1Ptl_8N83z+ zIE>DBdS1TDL8S)1MN63qhUP_WoGxeM#PaUA{cBAjQu$+@{2pt(U9WOz=`t#=WTdsd z(sstdS{a3VBx^YLBHa77hb2xATd%gTv%&gAsmJ5VO!#{3%%R9!7OIr=U*d(^U_$VL zF|!#FP_k0(u{q}qtgmt0kB@Y)>dyDh*%=n*c$YljZfBsg`O1mc+sGKQOyr?hgD$Qm z?RxqC9R)?Z%!6acIA|Ff*Hs=2FeP)p-8hYdS6(gI=@CSHKTy~E>KG6E_0wOJv~ZyK zWtn31s4=!}GS#~nVg(^Gt2|ePP;u88nKfsknMheamvv;DEhJBV=k%*^5ySQ@YzVXl z%C?m!uL@A0Zlr%pRG|}S+%ph5;Xp(AmmcPu_M4vs#BVWd)Hqyhc*TzZP$I#Ht&K7=laT%M=T*HY;^8( zF%^CHnks)+Vq%4t+5EESI6UzKCV!>T>c!d#C~V3+RX$eS=0M%C_lp4Fs<%iPRG|FCVar)#@Lg}Mzg z9G;Bk9bjT-le_uq%{1u0iJNLq1JwD4Yq`C00wtfT$H%(VA+WwZFk5aR$UJRpcm8RE zUT4L{qh=`3o#tA@)m{LxEA=FW>>ZIb_lk4)s1zow3_0$qsEkcKz0pOIY#bTVQAm5H zhvWY#Z)_{4;_|AEn|^+@#ys=Tsx9Z}&=^-L^qlUD3E8$$QV)4RT4te`UMq-s*G~FH zc3L8tbjkOp_&-|aYRY-fs{~G5?Axajzyk`)isBn#11j--4sP|Xm^{$;{)P$-5@Zy{ z4m`Gj1)nI=#xX3UH|@46yU)T2Kl0+mx8$I5|C{5dwpwBFXODoS^U6>-Ao${dHx&y% zRZ@k2IzYnxv73=?FZp;`ci~KeDIVEuv#9PU1+zUq=s2FX2ajJhL%+LJAmg3i`fn?l zAgSgm*1N$Ki&{^;PbcPX*dA?xtTU=Okt|x0d6EW>O69p5DvA5(DV*F~;|vWe6b=Mi zvZ34GFS4Rn2%arV9N+E2!qGca&YD;n$^^&yfAb^oP5A7On=Trd@KdWC4mzNN{s}ty z3=f}`{!m%t%LZQahD0HATfkp4ws{WUXtQ%BzL`7<$b2v^9}Q#R4T_KEk10AH{ov<3 z|HlcFJS3Ycw(($S+&3-wqBRDlK2(2w#RSY2%u2|)0os|=$e4R`&?HHx=g(?;>}+jP zH1P!t6HGB3vXTWE&D%jI6Uo?SVQqL~uQmAm+L&_bARV&0E^mJ8#lY-FyVAWEEx?`r zRPT?O6O2|HHFkgIpr`2BYcEWlkaXb5&v#EcX{*+yiI;@hVR7l?-chP4bgDbie;j+feM%(XPl^sxtZLaKZiY%;Qc*z( zn%GuVv%X`uD(aqU-MvLv8>EF6%h+Au;!uL1Y+ivX$Q-KD_#;Ds5QnKR=FSG#upvQx zbt?-?+ipEFt#ARK;gUetPF=V~A2f^&wZW*Z|GG+ID8TC$ar8Q<3a^tC!*zGqp%VMm zzv2ltaMs_^xg5>F4A-eg5rI66xCQ;;*$nuaN`EbT$sS0T=TY7NA%7;e_vVsq9L)GS zusiDxv5w%?1!u15L-94^=Vw2u;jzOReWgY;%-NaUDSN>Xz0*c?zH4*9R3Wl2E7t{N z^eQf|&R4`LCiTtpQab1s1$tV$>A}?YZvn9+J@g9r(X;)p6U_!#@ zpvT`AOF*gKxZ4tvJw6}WF`C4gBdj{2GJpRQ z1;WC|ev5uk0@Csq5`v$(U{ZNimyzR&Wy6c+TIHNk{<^Wjg&-Da^G+V_38J8|*i5zC zBWG+}{coe3vKgwkD4adLgo}R_w#j*~WJAWKq$B@)S#b5?8`EtI*f=YA@kve#4c>*l zmi~Gi&~cmTrU3&$gTmRD-i>rr`ugy0aX1x5qn}*$HUTf5hm2uJd+EyD>!zSmt~q^p7Ypjc=5O-faKX~M+bv6x2hSZX~|;Xy1H>jy_ff$&)KCyJGasMkm`<4DcDb+f?Glglr>qTTD>_HK(y%vsRYyx2z~! zVTKOkYMV-cuxHC3j{CDAr~cG| zlqWRkZEK3#w9y6I>{X1rwy_~mGEpEn+zv^T74>%AvwX9X3%nGLBc^_eJmGzW0o&Tk z-OGr1o|*W!`#p<|%`m%D`GYY=>$FccY;pwhM8m=5hbL(hN$sDNW7NT8QQHbHdkf57 zE#*?wBLeHnlQNA5Iq{t>*l{%`H6i!n*}H0L$jw3S)ojh-!-vy zvVe+N`g@P6z@`_vCA4TRWPkS;Nm<4QyH&Lt3bmXtCV#<(rhG?GT0#*H+Nlm8c;vrw zKUw6}%6+}y=LkcKC)|E{E25p3u8u2(4p*Mv{j>cX!Q-fMF8Rj|Kse4L%~r%3R{5kj z93}Q^-GeVFqf{=Qdg;lTNOi=q%iQ@3W-f4K*8JXHZ7TA@gFlT(TY=hXiv>po=|Iwt z{5-vTkj5A6_oM!`MpIkKrifxDnx5S!tY2>l9((plt$Hnm{zWb$MO0gucwl}q=c?HcD z)h~IdTuQRX%zu=WYqJgz_^9heaQr-NLc&gO`x*_TynWd0chDBry%qfMX^xI4`04ki z21gL5eR%cZG6R|~$Dd78(#N<3DckRVqC+vKRrbg>Ggvg~dPHvx2V@SsI~6u(1Q}zO zu0(3FaBy3qT+B{D=x_Y~Vps_$QiZSY9Mgs{osa|OXU#B7Xakdck`8r=e;)}3IsnVw z$@1STXUH%9m>|4X6{`f~Et}+=VPO6AdA@=XMD2+f^78;NZ9e0n)M*x&z zxTmydT&XsJ<=$s?+qzsZkVzJLlOlwHTE*O#1qyI-YZ$T)gOT9h0Q5tJTAGmNlnMgh&8oYg_3wVf$|Cv}Tg#Q1POr(ku z_wLIBhJK_mOx=3;@Vh!4pf+vBY`q}f_Kupz0xrZ^YkfaR=nclqrQgn@mY6Rs5MJ4* ziM*NaH>!)Ype^WJ{`@9olx^*-rjOcTg8j)|j~qoXV(9dGcM~3F`^=Qxd?*DqhvmZ0 zUlfDhkRPvht~N&2{BZs?>QTC#|7eC1=s^>Jc!y6WI1JHj$5m{ z3%*_A0Q6+n8zoAD<3#K3G-n1zl~IejF0x?aVEpPcDil1&aJJw0jE@E%$d7a3w%S5a9JVK%a@qHN=6ykq-*050v5FA+|Yp&Jy=toQcm2 zqC>?|<6HaWSs0e{Rr;e16&8QVw|waYNa~}%OMmSIb#cC*rfw0>b>^&zhPDXYUWBSe zW(=&NZFg;crUNxb&&#SJ<$*UVEgn;0f@9JW%Y5I{pmF1>$a4X<7F5@4WDM^vhF2@T#Z`E)BQkd0$imvbIzn?r`{kdxqdI)40D&n3OIL65i`PJW&g zh+CZ&x%im{L-i$;3zbv+!{!dqOBKKSj>g+FX8h%o99ljS1i)z%M_5_UkMdR%i6 z{zjwJ4a*2KdHA*Ir9*op4}{n4uosD?K=5o$^g|D8|RM zIav`HWdTp|t=!5=3pjVmp#SDRg0Gyk zsLAcnf--xFpLqiuG&H`id&$xq(qfMuy#Lw`1wWeh6jsoI>>25>OcP+&Sqsfx`NeFt4rZpmEcX7!f%BTejsnjKAMAX`scjD0;Jr4mmBqHg z&Qi&fN5?I2_G%LXHsG`0i9wsbZ|I6Tk23Qo0zESwV z#;Jsl5&wj>fWmnex^}A;NXG2G`Z1Y`Ulk62%w@2#zkAt-%u+kFOb$N$r`7>W{v1^l z)aPP1#gjI8gbk#3QF6n^TDbQ0#q-=d1kTZ^?HSb+@KVeg-rVQ1Q&Po>yHh*;wD zA{ATcW~~cI(Ut&}snYPA0UoGXizTi-#{|utnb;BKj5W`^^(^izfrb;|)v+0N5b2`) zF?+i+;>U)wYxilPVb=GCMj?BMx^!1e;tvmH`btxiE~vo6E|Tt@<1BDHZaA<~mVu0-pzio4dDC^gG(mgaR&DwPPb(|4q$WtX>w9s*7jfUqq6OdEy|9vi- z4W3FpK8g2vSiCd*#i=0?zzE+JW^`3tTfBDxd$kLa`wZv)O{=Xn0w`QPR+ZnO?VF@flo<(MwH= z`m2Qv-&Dr+jm;or<#wjaPhIFfKKo)0M9_n|L^$NM7IGs4jP8En;Iq@cZ!5K!7#%F2 zxzv^lPMT+>Q=JIk;C}1qhH822KC1HKmHZL}-F=(H9Y*=1OIJ}M{?M?|DOOH-DI4{> z_*)No5T zqYDGyRy%=wT*g`{cP$Jo7dvdCFAmA?k61`;A!DS@`A|JCEeP!W7y4lFB)>4*MkZau z9Q_Z7pQ|2bz=P?>E4bYRj+amR-piu^Dd>4+=BKy(d^?FZd{Z7$Wmg)wJ+;D+Y>mwe zdMV(M>C|NZR~rQVO&7eIqJwa_k=Tc|Mwpgy+)`-eVp#Y1rN--OE}oeEe0`A6<5Wk% z8~Yy0V41RV_AO&B+&eFK((^PE#PdtCjrK9nF6FOJ*%TY)4=otgSw`TFeWiw1v<6Bo zdd2F@c0_TNT`e75)=0hU^~`mWz`-%q@c$S~fYiP~W@y_Oe{@)noOjO#jYCtv=ssg( z)Ya7XXqFCAKgKypmT02k37xJIVR;EH6_Ecnv)V%_9F zCg%Nhi}+!{1ExNo)~LcoLEkTb7SRBJ@8h{|)EUzbllx+|tWbAl@6lj4Cn&{Lg?0yR zAR#GDI&1m|AFhQPN4-!7mSp9T9}HKF(PqBh80`vFh59vHu4qB_ke-tHFB)9G|Lx*= z5*q_;&m^`O+C#sV_b|PMhLIvo`WxkK;Nf>EaXUGJKeX6;m~R~BlagA-J|9{Dkzb2(}9d z)1v-*+}1|5hejp~ZRkiA?^-57_{=W4KhxO<37=kQp5pma8xqVSYNqRc@E_i7K}SMY z8r%-4uGHnB663P^RpNYR-_IO5QmF@|hL@Tv>L1f))gOIbD#FHvOIEyFK;%B7FJ8Z` zs3d&Qt##L%9s|^z)?GuFvVqb%z3tJ>6x8*kT(6&|_Fib~ge~f|?cxYU` z=sU!M{It{g>E0TUka3i=d`<(dE7<1Ea|wN&wP6)2kq41W{%y^CW)8mt(eu473q+zi zPSOpWP={tWE~yD9bL_vo+*}2;47lGd%#kPdx1}!4%Nl;Sx5rsswMG36)waWZbjZWY z?TXP>=$QE1;?O+7C#*lYKbW%t}Nv9 z?&kA{EFduYx_9rqF#b%`J9VuMU&b3l+iRGi&XjptbqcuaoL<9G%fukVQ`_ zo&WtWl7gsEXwn|131dg7?-{*1kWH2A`yxrfkbZ0aE=5P!fBa9LL9aEQo0ZAc?6br` zi8)mrVM3=HpIExd)DpVG{)^wCWd`zRPMuy|#DVQkBv)qJFmdZUhb;$!X(%Q2!Lr(v ziKjKGj~KB;-lpyGXVgRtObb;89{E~AVE@3)6-ByW^J_}2W{Qa!>->s4MTxwn-S1GV zb1$Ei-=h%v<~_~5MF9$uc%U>@^MicF5~(BqR%dpb!itZxb)JOZ(#>)UuO3u}O`h@o zfwwu(Oo=huc1suKGj}nW6f21N{(V7LI~nyAg*#8R(=n~5x6|T@HcaFktMnNa2GNh= z))^CYSZ>5Ec{fA&jnS{un3Gm$DR5KJ`7IMF-j@ke?5%*b!N1UdA03nRf;t90wXihu zdQ-G24YTtCrDbGY!SB(#^7#|a*nQP+G+csuT9+K$?IX=C6mWK=$3!mlB<0~1H6 znC)_O@|OUSV_%T+!Zcz&q>i)QS)za{qnUx;A8A#rrDlqA78o5C*C0XU+Fri5X?@8o z97{xhs*^R$70@Nr3Hgno7Z+toF@UP9k*@K91>JuOSDm=5fgm0(v{`__uSvV%)+|D= z1Vm3AIzoo*!JK+iA_p{FRCy+5j{`(ClTOIrvxZ3hNAvXiiU`YJ1a5Y)gRjEpUB(5i zp-%n%=N4l!%6vR2yevr$^56fV7pqucm4%C~^CSg0&&w8-ujZn#ZwEam+zJXNL$2th zxT0aYRngHTE4Vc2u{wd^jp(x3Vx_|`T4eU7^}KHukg$D4`OG$HSo?{?GB@F1f9aFf zuy`82%$Kh%NM#`3ck8okDMw5)P6|o#vBCxR|4pWJGB7fAZ%^N}Cd$Z5TYCR8#_qVN z;n{;OP?*=FGq{9@|7jbX4t`{f`G-b3vIfo2NA7Y$M6oT*-gVQa{jdU$yQdo(4_M8H}8=A^dQ{RzmAI7l2Qt(M_>ZB@kJe{9WCw!%-7fB#H z#|Cm|m}>>)(0gpouDxS8m-q$wb=8`}+8ujsDe~pFOCed+Am?TQHbEjovRy@h6p1*(|j-4<2zsc$K1u-S{qFJNe!+)rH&!@ z#4jeD(1YxiQOEvG=m6&2c~32n0_DDm;en@Sp4ag7B%olDV zV?vAzCoqKI>V5EAK;0pOGuJ9gle7`%WH;Dr*N? z30>bN`*f#d=wf_aZ5{0Xk_9!1ANO7CA``jOQp(pn8j^=g+jpjN@I+VVstPWV=S5Wt zxNPNNT5`i{?;#cF%@OV{U1ETyGsno2xh$AZ&KNl=W{2I6l4|;OoN!PyZ0qMo#CrNi z*%(k8&Mi{)Ng3)o2A>)TrU#1NI&Ei@K5V}p*T4B2X5*|2bP zS#dq#pR`VzEd1KeMvu+0p}v*M5aQo?PNa+=c&g;_#&lcEJRh?oWupmt?YA7T8YXmg zOrlZbQYJo>J-qybzBN*AzkU4v9g#mCz1osi!$pT5uW9YX+{^spnf;%;BP>uj@uafF z8Y902}|BbeU+hLTT>kk1d&U78o zbR>9<@6neFMhvjgW>dAx0u$_dy*6v%7#)3PD%~#XaB#$fF|$p{7L}IUY!Njf@({<+ z@{ndG<|mZ7&Gk6L1Df0C=_C%^-syS%ethN=S%{YCl3Dr~s2_@9aCYiC^j>F#mzB?6a=tWLL< z&hWcKUBvRT6=9J@mB2u?Gr~);2Nz?+U^2MCt8u*}GHP_4T9bJgMY0YwRACR?b<&fNqc4Uxg83K$UvZ_@Gt^%1Uyd z7OiDL|BN7di4s0+%9V#Yt{kM!mJg}>6F%Pznn(vLLz_TsNsW;)p{pzhW>Rf%LP$1n z_N4}rSG-M5GiJa;&5oUdgl*8B+rt`+7lZ5-6WxyIwXoxYLZV_g9lE#JM%NIz4XH0O zBju(J;xmTK_{53y5HqGNQIAk z>s*SC08Tbv`Z?XH4_8Y*Z&D%hpvfJNVk1}BgZ-Vw-hw}AIF&ORx-Ny_F_&8dKGqPq z{5+p#EoKFgJ8l|8?6iRG$1L0RuMj?IW>J$_zaw@JQB5~ay1-zwu^2be6zjj=@;f#H zm|4O+>sN1pm6hI?D8;UX{kwGb(qRWo5YK9mSY!uNZ-WkIYT3Z@6&ZWl37tI@c&E7T zj|ddrk*&Y>hw$I4TrEmc0B6S@dYo}JMUshnZQ(2plC4@!mLwWu)x(HItwt35xz>sA z?&XM9H!>fqMHBe@M#;8Wc%1h9M%%`=E_>|mA31770zY)QLfzwPpk5zyrUr13LKFFa1ElCkD|kGk-BsZUL0ql>Kh* zWaNEG5_b>hK*ERff1WBDpxnk&tSLgjc0$NQf);|H@&(-mY77Vq%XUAoNgoy#s}J0; zF(!1QccTiG3TcUUWk%yf-jV+KqcNF@nWpq3Vop38kIA#H|w4*!KD3p4P9bU?DF z95EMK)Vq38ox!km&MC3U2|5?_9Cs-ubSss;$nCrZ26pru@cZuvf1|Us-!Kc^ZDb;;C4X$ARdH zwJis_IZ*R&Z?*q=WiZn)IdwCO=zGM!Zo2Wu9*^cFBs@85ho$!2O$-$tHj3|ey*vOQ zU(FT{?R5ase7}WE#myGCUtl$^`ZCM8Y<@NRi3#rYK0A>E^ilDn!>pAhIW@s zDz1Jxd-(~KgK2a2yLW7`hJ0T5!-iCV3Dc|%7E_wAB4Y2^{O5qqSEL`u6*HkqVfo>F zBG>p*N{MX$$v{Kj!GQExGP<)E8f*BjP`9g0Ft-J8-4c@no%sxy*>yHbq>qB@?!MYv zp8;^~()h1_a7BxCMR%H9>~QdYNytRe|8P8CS2=c?KY7+phGS2K#)c355|PgE@=eb~ zvnUH@i#jE55c$wt-+%odHxhYEhrptpht*+l=doBnH%B0CzOh{L(I1-O+~ionDs7w? zkJHp|U_xfzlJE__CScxZCR}X9L95n@ZRt;mUdD=r5h|gM2@URHV1Dq-vO7JN zD0C~x4vBS|zf~>_11_+-|M0$l05~9dpmtj(9fqSuqXv%>dhwY1M9)httSG$eL-@Tw zU1rDXeCi6_4xdw3Jk$cESMi@;#+gI6R>Wt8H5NFXe|M9c0v&~q8N@F-q=zdXDBhV& zvPS3T=VVhNhgA}MmOxu!g*}cJmrZ&T`r*pH>EyGfU|G#rbz6g&iz zTn^C#?3>&GK)3hno_sf2p)4{{SDF2b{0_5p@)swhx z0l^cBxa$B9@|Dvk$JP`6R0(rw(HRyF^e)`x`ON}l{^3%GcSg_|ps>Fq&I+?FrTz(> zqJ#W&>2@_DCmS>T9OwSm5@ugzl)OB|0Ni5Nm3mzUYc$i+_rHFd!-Il6R1<1I#a<8^56GZjB(_yh79(xx+h2_w+r;7(C^nOZibtH5u4ELw36NQlHCYMzoIO3zGuB>PQ8_W#I z+^n5#2)k0&t7o%Wpi!SUw*8wW!X#w1YK-tnUpS;Uxdb0f4Y?#f;RIC6X1$$%WzcX- zjW6wj4S{nvU#TB)MAJPq^Z90H%)DN-@5w$BSgiWbj1j_u;?~L|$KTP3d@R+bG~W`2 zihgxJx=ezpPI1w-igq~GBp-G#hz4Ol8(2HltRT&{C%R6-4m2j#*v)Tu#WSq+M})3& z(Y0cu=(ZP5SV7N{_kmFRz1#TICb z5OW~u`(CjWTPTxp*L#*_h0AoVG>eqbabdV&onp2&nWl6Q32mX2+?Nhz61nxHhV7<@*K@T&vv`C8`?eu{ON zc3>0H$L5bj|04J;?0Io7pXgUL7fgxBtuV%BVM=NGFCKIz54SqKw*{Z(n2w22Jxr6f zFOT1^3_0WbSAHrY`Xhoa!?iQ!Fj2$3qeeG}ug%K_Ms%DZoWE1Xk8K0Ja4q*(EfiKDFtd3^ww!!SD1(L=YzjWA*7p7pv0MvzpSTM^v zBP3y>&zH1H^RCWTymW=a4qe<8xQr9gMC2^d^L2Z_5AjLternH~ zPc+KmP?p+39jtk;df;`o9l)ylDu0WXz--p{m=glJkT81gp@B3Lsa8gpZa8uwka^HC zar19}TGJus)FvWFeI`_wnQw`fx~I;ox41y$wUoGIu?0x6;wUr*+k@rpGrgnU3MkGl zzUA`E6yDqKFfti)!u)#|7e|mZ!1OMaKS8XU=la^l`ae6cv^YmyhGy{jt%$6MxIK>D ztL}ME^xYgdazXmuj@aoUIrbq>2i!dq!W4Hf!Lss+mgA%tlI-3s+bj8rPude;O;`TO zC!gNrzE|fjA9EeH>9}wqS=3W|yPY*0xL4zs_lr1>_1iy;#frcc6a62fM9m84jTfCA zWI^_{9w&xl0Z-2tkJb7xPI>x>;3)9^Nmexp{D(WBoio0-|tzwiicj8vZ`&P_29=R?*(EV zOwgj*Nj7d|Vq@>M$CRB!|1j0R;B+IATc}BQq^z<+s`TqSOS_3aaqGUP>at^ejiJLs zLqsn)dvyEuhOO#w-hNT)ij5pxD78L)`ZXDoPrS+q&S65`#goi?EDF)HrS6Iu=Rs_A zN&M=OI~f;-XWnEleH1;GQaGiMw9JzLy;4VE5cQ-L^$c zV2%qkt@y@6ja`o;t9Wc60>3xkW~7kuHde3xkUdf&(>=p3vB3V({H{GJju_S{;UrLD z4w(-+&xzZxA)9qLF?_!$q_u@-U4Lu|V+*-UuJhILL}WzBVSff>?t6AaHqH?OpYIQ^ z-mw^tnnevO^)R6=?D7`(U`1Y@++Y|f11p`{*LK*V_318A8 zo~iw%W&G-C`~z zylk7;x#bQ2?y7~xR6>W7zO`(wB-hbA>V6b?pJV}rNg*2@Vqu7;{P!l3J$#p1vrc&( z2O}LKYabImzJd6z6S7%E?tNwL}tzz40*ki|zrFOEV_V`tB zV@And6Y#kiMk+L9!H}ePzUL)1Ovuj`jyYn2q)W;!I|d1V|QC7lKiF=iw$2etZ0CihpUIBp*QRD~P%b2#% zIsaSK<^TmI@9Hl1Cwy7f-w(g^uLFwTIChZoTphDt&D*4fDuTH0m+yPjoZ!u9jh*Lx zCrl0xO0k>KA^vY=x;*hN8F&Y&hY~*8!EdQKI8W%&g_*`6m&ih=SJJ&g`za`MQ>MbM z-v%bz$TY3hR=8xzs+w)#9K67eJ1_Z=iyeRDBF*KAK4|;zpu2~3U{Fp}N9_gCgU+98 g&QepuD#NGlkrM`Js9BcEyCP5YYn)eXYWzj}KkRNIcK`qY literal 0 HcmV?d00001 diff --git a/pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D8/m0.npy b/pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D8/m0.npy new file mode 100644 index 0000000000000000000000000000000000000000..31f9f7bacc4ac9b9fef56277143b5f32e9b01f89 GIT binary patch literal 13952 zcmeIyAr8V&3LFj3m#*&f?p)sLexA%0^E{ho=_~!U z#>r5{yU^{y9ryzdaDW3G-~b0WzyS_$fCC)h00%h00S<70103K02ROh14sd`29N+*4 hIKTl8aDW3G-~b0WzyS_$fCC)h00%h00S^3^17FC()i(eD literal 0 HcmV?d00001 diff --git a/pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D8/m_true.npy b/pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D8/m_true.npy new file mode 100644 index 0000000000000000000000000000000000000000..c0262740d6e085468b5bfe565966e6574679358e GIT binary patch literal 13952 zcmeIyAr8Vo5Cza;ImNn#qy`KPL7;FD3<5 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/pygeos-tools/examples/solvers/gravity/modeling/README.md b/pygeos-tools/examples/solvers/gravity/modeling/README.md new file mode 100644 index 00000000..24c0dc8f --- /dev/null +++ b/pygeos-tools/examples/solvers/gravity/modeling/README.md @@ -0,0 +1,25 @@ +# Gravity Modeling + + +## Typical Usage + +### Rectangular Prism Anomaly Meshed with Tetrahedra + +```bash +python gravityModeling.py --xml ../data/rectangularPrism_C3D4/rectangularPrism_C3D4.xml --m_true ../data/rectangularPrism_C3D4/m_true.npy +``` +or equivalently: +```bash +python gravityModelingLinearOp.py --xml ../data/rectangularPrism_C3D4/rectangularPrism_C3D4.xml --m_true ../data/rectangularPrism_C3D4/m_true.npy +``` + +### Rectangular Prism Anomaly Meshed with Hexahedra: + +```bash +python gravityModeling.py --xml ../data/rectangularPrism_C3D8/rectangularPrism_C3D8.xml --m_true ../data/rectangularPrism_C3D8/m_true.npy +``` +or equivalently: +```bash +python gravityModelingLinearOp.py --xml ../data/rectangularPrism_C3D8/rectangularPrism_C3D8.xml --m_true ../data/rectangularPrism_C3D8/m_true.npy +``` + diff --git a/pygeos-tools/examples/solvers/gravity/modeling/gravityModelingLinearOp.py b/pygeos-tools/examples/solvers/gravity/modeling/gravityModelingLinearOp.py index 304de08e..9995a05a 100644 --- a/pygeos-tools/examples/solvers/gravity/modeling/gravityModelingLinearOp.py +++ b/pygeos-tools/examples/solvers/gravity/modeling/gravityModelingLinearOp.py @@ -17,7 +17,7 @@ import numpy as np from geos.pygeos_tools.input import XML -from geos.pygeos_tools.solvers import GravityLinearOperatorSolver +from geos.pygeos_tools.solvers import GravityLinearOpSolver from pygeosx import run, COMPLETED __doc__ = """ @@ -49,7 +49,7 @@ def main(): m_true = np.load( args.m_true ) - solver = GravityLinearOperatorSolver( rank=rank, xml=xml, nm=m_true.size ) + solver = GravityLinearOpSolver( rank=rank, xml=xml, nm=m_true.size ) print( f"Density min={np.min(m_true)}, max={np.max(m_true)}", flush=True ) From 46148441f45e5e59093994633229bed87bccb497 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Thu, 3 Jul 2025 19:17:23 -0500 Subject: [PATCH 5/7] Added data for gradient and adjoint test --- .../solvers/gravity/adjoint-test/README.md | 27 +++++++++ .../solvers/gravity/gradient-test/README.md | 55 +++++++++++++++++++ 2 files changed, 82 insertions(+) create mode 100644 pygeos-tools/examples/solvers/gravity/adjoint-test/README.md create mode 100644 pygeos-tools/examples/solvers/gravity/gradient-test/README.md diff --git a/pygeos-tools/examples/solvers/gravity/adjoint-test/README.md b/pygeos-tools/examples/solvers/gravity/adjoint-test/README.md new file mode 100644 index 00000000..469f101f --- /dev/null +++ b/pygeos-tools/examples/solvers/gravity/adjoint-test/README.md @@ -0,0 +1,27 @@ +# Adjoint Test for GEOS Gravity modeling + + +## Typical Usage + +### Rectangular Prism Anomaly Meshed with Tetrahedra + +```bash +python gravityAdjointTest.py --xml ../data/rectangularPrism_C3D4/rectangularPrism_C3D4.xml --n_model 10368 --n_data 1271 +``` + +### Rectangular Prism Anomaly Meshed with Hexahedra: + +```bash +python gravityAdjointTest.py --xml ../data/rectangularPrism_C3D8/rectangularPrism_C3D8.xml --n_model 1728 --n_data 1271 +``` + +## Typical Output + +```text +=== Adjoint Test Summary === + = 7.585858080970467e-07 + = 7.585858080970469e-07 +Passed = True +Error = 1.2332569442944734e-16 +``` + diff --git a/pygeos-tools/examples/solvers/gravity/gradient-test/README.md b/pygeos-tools/examples/solvers/gravity/gradient-test/README.md new file mode 100644 index 00000000..4795b34b --- /dev/null +++ b/pygeos-tools/examples/solvers/gravity/gradient-test/README.md @@ -0,0 +1,55 @@ +# Gradient Test for GEOS Gravity modeling + + +## Typical Usage + +### Rectangular Prism Anomaly Meshed with Tetrahedra + +```bash +python gravityGradientTest.py --xml ../data/rectangularPrism_C3D4/rectangularPrism_C3D4.xml --m_true ../data/rectangularPrism_C3D4/m_true.npy --m0 ../data/rectangularPrism_C3D4/m0.npy --dm ../data/rectangularPrism_C3D4/dm.npy +``` + +### Rectangular Prism Anomaly Meshed with Hexahedra: + +```bash +python gravityGradientTest.py --xml ../data/rectangularPrism_C3D8/rectangularPrism_C3D8.xml --m_true ../data/rectangularPrism_C3D8/m_true.npy --m0 ../data/rectangularPrism_C3D8/m0.npy --dm ../data/rectangularPrism_C3D8/dm.npy +``` + + +## Typical Output showing the expected quadratic error convergence up to machine precision +```text + h h2 e0 e1 +5.5556e-01 3.0864e-01 6.8482e-04 1.4697e-07 +3.0864e-01 9.5260e-02 3.8049e-04 4.5360e-08 +1.7147e-01 2.9401e-02 2.1140e-04 1.4000e-08 +9.5260e-02 9.0744e-03 1.1745e-04 4.3210e-09 +5.2922e-02 2.8008e-03 6.5249e-05 1.3336e-09 +2.9401e-02 8.6443e-04 3.6250e-05 4.1162e-10 +1.6334e-02 2.6680e-04 2.0139e-05 1.2704e-10 +9.0744e-03 8.2346e-05 1.1188e-05 3.9214e-11 +5.0414e-03 2.5415e-05 6.2157e-06 1.2104e-11 +2.8008e-03 7.8442e-06 3.4532e-06 3.7367e-12 +1.5560e-03 2.4211e-06 1.9184e-06 1.1550e-12 +8.6443e-04 7.4724e-07 1.0658e-06 3.5841e-13 +4.8024e-04 2.3063e-07 5.9211e-07 1.1246e-13 +2.6680e-04 7.1182e-08 3.2895e-07 3.5991e-14 +1.4822e-04 2.1970e-08 1.8275e-07 1.2474e-14 +8.2346e-05 6.7808e-09 1.0153e-07 6.0212e-15 +4.5748e-05 2.0928e-09 5.6404e-08 3.2005e-15 +2.5415e-05 6.4594e-10 3.1336e-08 2.5012e-15 +1.4120e-05 1.9936e-10 1.7409e-08 2.8565e-15 +7.8442e-06 6.1532e-11 9.6715e-09 2.3721e-15 +4.3579e-06 1.8991e-11 5.3731e-09 2.8261e-15 +2.4211e-06 5.8615e-12 2.9850e-09 2.9543e-15 +1.3450e-06 1.8091e-12 1.6583e-09 2.4471e-15 +7.4724e-07 5.5837e-13 9.2130e-10 2.8471e-15 +4.1513e-07 1.7234e-13 5.1183e-10 2.1396e-15 + +=== Gradient Test Summary === +Passed: True +Slope: [ 1.99999993 1.99999983 1.99999957 1.9999963 1.99999456 1.99997363 + 1.99988094 1.99984551 1.99960117 1.99755535 1.99079356 1.97189443 + 1.93835375 1.80267853 1.2392288 1.07519441 0.41943078 -0.22598208 + 0.31615789 -0.29794799 -0.07551366 0.32048653 -0.25757919 0.48605047] +History saved to: grad_test.txt +``` From 1213ea83ec0a1bf8c83bef4baf0a4a04c12f5203 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Fri, 4 Jul 2025 16:09:53 -0500 Subject: [PATCH 6/7] Fixed mypy. --- .../pygeos_tools/solvers/GravityLinearOpSolver.py | 6 ++++-- .../src/geos/pygeos_tools/solvers/GravitySolver.py | 5 +++-- .../geos/pygeos_tools/solvers/InversionUtils.py | 14 ++++++-------- 3 files changed, 13 insertions(+), 12 deletions(-) diff --git a/pygeos-tools/src/geos/pygeos_tools/solvers/GravityLinearOpSolver.py b/pygeos-tools/src/geos/pygeos_tools/solvers/GravityLinearOpSolver.py index 6be54e3c..7ff0e5c2 100644 --- a/pygeos-tools/src/geos/pygeos_tools/solvers/GravityLinearOpSolver.py +++ b/pygeos-tools/src/geos/pygeos_tools/solvers/GravityLinearOpSolver.py @@ -12,9 +12,11 @@ # ------------------------------------------------------------------------------------------------------------ import numpy as np +from numpy.typing import DTypeLike from typing import Optional, Tuple from scipy.sparse.linalg import LinearOperator +from geos.pygeos_tools.input.Xml import XML from geos.pygeos_tools.solvers.InversionUtils import computeResidual, computeL2Loss from geos.pygeos_tools.solvers.GravitySolver import GravitySolver @@ -49,10 +51,10 @@ class GravityLinearOpSolver( LinearOperator ): def __init__( self, rank: int = 0, scaleData: float = 1.0, - xml: Optional[ object ] = None, + xml: Optional[ XML ] = None, nm: int = 0, nd: Optional[ int ] = None, - dtype: np.dtype = np.float64 ) -> None: + dtype: DTypeLike = np.float64 ) -> None: if nm <= 0: raise ValueError( f"Invalid dimensions: nm={nm}, must be a positive integer." ) diff --git a/pygeos-tools/src/geos/pygeos_tools/solvers/GravitySolver.py b/pygeos-tools/src/geos/pygeos_tools/solvers/GravitySolver.py index e86ed6c6..f4b59aee 100644 --- a/pygeos-tools/src/geos/pygeos_tools/solvers/GravitySolver.py +++ b/pygeos-tools/src/geos/pygeos_tools/solvers/GravitySolver.py @@ -13,7 +13,7 @@ import numpy as np import numpy.typing as npt -from typing import Dict, List +from typing import Dict, List, Optional from typing_extensions import Self from geos.pygeos_tools.solvers.Solver import Solver from geos.pygeos_tools.wrapper import get_matching_wrapper_path, set_wrapper_to_value, allgather_wrapper @@ -55,7 +55,7 @@ def getGzAtStations( self: Self ) -> npt.NDArray: """ return self.getGeosWrapperByName( "gzAtStations", [ "Solvers" ] ) - def getDensityModel( self: Self, filterGhost: bool = False, **kwargs ) -> npt.NDArray: + def getDensityModel( self: Self, filterGhost: bool = False, **kwargs ) -> Optional[np.ndarray]: """ Get the density values WARNING: this function aims to work in the specific case of having only 1 CellElementRegion in your XML file @@ -81,6 +81,7 @@ def getDensityModel( self: Self, filterGhost: bool = False, **kwargs ) -> npt.ND return density_filtered else: print( "getDensityModel->filterGhostRank: No ghostRank was found.", flush=True ) + return None else: return density else: diff --git a/pygeos-tools/src/geos/pygeos_tools/solvers/InversionUtils.py b/pygeos-tools/src/geos/pygeos_tools/solvers/InversionUtils.py index fff84642..9b034b02 100644 --- a/pygeos-tools/src/geos/pygeos_tools/solvers/InversionUtils.py +++ b/pygeos-tools/src/geos/pygeos_tools/solvers/InversionUtils.py @@ -58,7 +58,7 @@ def computeL2Loss( predicted: np.ndarray, observed: np.ndarray, scale: float = 1 residual = computeResidual( predicted, observed ) loss = 0.5 * np.linalg.norm( residual )**2 - return scale * loss + return float(scale * loss) def gradientTest( f: Callable[ [ np.ndarray ], float ], @@ -128,7 +128,7 @@ def gradientTest( f: Callable[ [ np.ndarray ], float ], if ind.size > 1: slope = np.diff( np.log( e1[ ind ] ) ) / np.diff( np.log( h[ ind ] ) ) mean_slope = np.mean( slope ) - passed = ( isclose( mean_slope, 2.0, abs_tol=0.1 ) or np.count_nonzero( slope > 1.9 ) > 2 + passed = bool( isclose( mean_slope, 2.0, abs_tol=0.1 ) or np.count_nonzero( slope > 1.9 ) > 2 or np.all( e1 < 1e-15 ) ) else: slope = np.array( [] ) @@ -174,8 +174,6 @@ def plotGradientTest( h: np.ndarray, warnings.warn( "matplotlib is not installed." ) return None - import matplotlib - matplotlib.use( "Agg" ) # Safe for headless environments import matplotlib.pyplot as plt # Filter out invalid values @@ -255,8 +253,8 @@ def plotGradientTestFromFile( filename: str, return plotGradientTest( h, e1, **kwargs ) -def adjointTest( A: callable, - AT: callable, +def adjointTest( A: Callable, + AT: Callable, x: np.ndarray, y: np.ndarray, tol: float = 1e-13, @@ -266,9 +264,9 @@ def adjointTest( A: callable, Parameters ---------- - A : callable + A : function Function that computes the forward operation A(x). - AT : callable + AT : function Function that computes the adjoint operation A^T(y). x : np.ndarray Input vector for the forward operation. From 2526cb9e004e9539534560823fb65b6bb975e104 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Tue, 8 Jul 2025 13:21:46 -0500 Subject: [PATCH 7/7] Update. --- .../adjoint-test/gravityAdjointTest.py | 2 +- .../rectangularPrism_C3D4.xml | 3 +- .../rectangularPrism_C3D8.xml | 3 +- .../gradient-test/gravityGradientTest.py | 2 +- .../solvers/gravity/modeling/README.md | 14 +++-- .../gravity/modeling/gravityModeling.py | 36 ++++++++---- .../modeling/gravityModelingLinearOp.py | 55 ++++++++++++++----- .../pygeos_tools/solvers/GravitySolver.py | 7 ++- .../pygeos_tools/solvers/InversionUtils.py | 7 +-- 9 files changed, 88 insertions(+), 41 deletions(-) diff --git a/pygeos-tools/examples/solvers/gravity/adjoint-test/gravityAdjointTest.py b/pygeos-tools/examples/solvers/gravity/adjoint-test/gravityAdjointTest.py index ed8314e8..7c117a34 100644 --- a/pygeos-tools/examples/solvers/gravity/adjoint-test/gravityAdjointTest.py +++ b/pygeos-tools/examples/solvers/gravity/adjoint-test/gravityAdjointTest.py @@ -21,7 +21,7 @@ from geos.pygeos_tools.solvers.InversionUtils import adjointTest __doc__ = """ -This is an example of how to set up and run an adjoint test for GEOS simulation using the GravitySolver. +Example script demonstrating how to configure and run an adjoint test for GEOS gravity simulation using the GravityLinearOpSolver. """ diff --git a/pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D4/rectangularPrism_C3D4.xml b/pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D4/rectangularPrism_C3D4.xml index fc39ee68..5b7d9213 100644 --- a/pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D4/rectangularPrism_C3D4.xml +++ b/pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D4/rectangularPrism_C3D4.xml @@ -72,8 +72,7 @@ targetRegions="{ region }" mode="adjoint" logLevel="1" - outputGz="1" - outputGzBasename="gz_C3D4" + outputGz="0" stationCoordinates="{ { 0, 0, 10 }, { 500, 0, 10 }, diff --git a/pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D8/rectangularPrism_C3D8.xml b/pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D8/rectangularPrism_C3D8.xml index 5a40c18e..fe64ac73 100644 --- a/pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D8/rectangularPrism_C3D8.xml +++ b/pygeos-tools/examples/solvers/gravity/data/rectangularPrism_C3D8/rectangularPrism_C3D8.xml @@ -72,8 +72,7 @@ targetRegions="{ region }" mode="adjoint" logLevel="1" - outputGz="1" - outputGzBasename="gz_C3D8" + outputGz="0" stationCoordinates="{ { 0, 0, 10 }, { 500, 0, 10 }, diff --git a/pygeos-tools/examples/solvers/gravity/gradient-test/gravityGradientTest.py b/pygeos-tools/examples/solvers/gravity/gradient-test/gravityGradientTest.py index ca68499b..efd24a3b 100644 --- a/pygeos-tools/examples/solvers/gravity/gradient-test/gravityGradientTest.py +++ b/pygeos-tools/examples/solvers/gravity/gradient-test/gravityGradientTest.py @@ -21,7 +21,7 @@ from geos.pygeos_tools.solvers.InversionUtils import gradientTest, plotGradientTestFromFile __doc__ = """ -This is an example of how to set up and run a gradient test for GEOS simulation using the GravitySolver. +Example script demonstrating how to configure and run a gradient test for GEOS gravity simulation using the GravityLinearOpSolver. """ diff --git a/pygeos-tools/examples/solvers/gravity/modeling/README.md b/pygeos-tools/examples/solvers/gravity/modeling/README.md index 24c0dc8f..2606dc03 100644 --- a/pygeos-tools/examples/solvers/gravity/modeling/README.md +++ b/pygeos-tools/examples/solvers/gravity/modeling/README.md @@ -5,21 +5,27 @@ ### Rectangular Prism Anomaly Meshed with Tetrahedra +To use the density model defined in the XML file: ```bash -python gravityModeling.py --xml ../data/rectangularPrism_C3D4/rectangularPrism_C3D4.xml --m_true ../data/rectangularPrism_C3D4/m_true.npy +python gravityModeling.py --xml ../data/rectangularPrism_C3D4/rectangularPrism_C3D4.xml --save_gz gz_C3D4.npy +``` + +Alternatively, override the XML-defined model by providing a custom model as a NumPy array: +```bash +python gravityModeling.py --xml ../data/rectangularPrism_C3D4/rectangularPrism_C3D4.xml --model ../data/rectangularPrism_C3D4/m_true.npy --save_gz gz_C3D4.npy ``` or equivalently: ```bash -python gravityModelingLinearOp.py --xml ../data/rectangularPrism_C3D4/rectangularPrism_C3D4.xml --m_true ../data/rectangularPrism_C3D4/m_true.npy +python gravityModelingLinearOp.py --xml ../data/rectangularPrism_C3D4/rectangularPrism_C3D4.xml --model ../data/rectangularPrism_C3D4/m_true.npy --save_gz gz_C3D4.npy ``` ### Rectangular Prism Anomaly Meshed with Hexahedra: ```bash -python gravityModeling.py --xml ../data/rectangularPrism_C3D8/rectangularPrism_C3D8.xml --m_true ../data/rectangularPrism_C3D8/m_true.npy +python gravityModeling.py --xml ../data/rectangularPrism_C3D8/rectangularPrism_C3D8.xml --model ../data/rectangularPrism_C3D8/m_true.npy --save_gz gz_C3D8.npy ``` or equivalently: ```bash -python gravityModelingLinearOp.py --xml ../data/rectangularPrism_C3D8/rectangularPrism_C3D8.xml --m_true ../data/rectangularPrism_C3D8/m_true.npy +python gravityModelingLinearOp.py --xml ../data/rectangularPrism_C3D8/rectangularPrism_C3D8.xml --model ../data/rectangularPrism_C3D8/m_true.npy --save_gz gz_C3D8.npy ``` diff --git a/pygeos-tools/examples/solvers/gravity/modeling/gravityModeling.py b/pygeos-tools/examples/solvers/gravity/modeling/gravityModeling.py index 169e30b1..91ee91a8 100644 --- a/pygeos-tools/examples/solvers/gravity/modeling/gravityModeling.py +++ b/pygeos-tools/examples/solvers/gravity/modeling/gravityModeling.py @@ -21,20 +21,27 @@ from pygeosx import run, COMPLETED __doc__ = """ -“This is an example of how to set up and run your GEOS simulation using the GravitySolver. +Example script for running a gravity simulation using GEOS and the GravitySolver interface. """ def parse_args(): - """Get arguments + """ + Parse command-line arguments for the simulation. Returns: - argument '--xml': Input xml file for GEOSX + argparse.Namespace: Parsed arguments including: + --xml (str): Path to the GEOS XML input file (required). + --model (str, optional): Path to a .npy file containing the density model. + If not provided, the model defined in the XML will be used. + --save_gz (str, optional): Path to save the computed gz output as a .npy file. """ - parser = argparse.ArgumentParser( description="Gravity simulation example" ) - parser.add_argument( '--xml', type=str, required=True, help="Input xml file for GEOS" ) - parser.add_argument( '--m_true', type=str, required=True, help="True model (.npy)" ) - + parser = argparse.ArgumentParser( description="Gravity modeling example" ) + parser.add_argument( "--xml", type=str, required=True, help="Input xml file for GEOS" ) + parser.add_argument( "--model", type=str, default=None, + help="True model file (.npy). If not provided, the density model from the xml will be used") + parser.add_argument( "--save_gz", type=str, default=None, + help="Optional output file to save gz (.npy)") args, _ = parser.parse_known_args() return args @@ -50,10 +57,19 @@ def main(): solver = GravitySolver() solver.initialize( rank=rank, xml=xml ) - density = np.load( args.m_true ) + if args.model is not None: + model = np.load( args.model ) + print( f"Density min={np.min(model)}, max={np.max(model)}", flush=True ) + else: + model = None + + gz = solver.modeling( model ) - gz = solver.modeling( density ) - print( f"gz min={np.min(gz)}, max={np.max(gz)}" ) + if rank == 0: + print(f"gz min={np.min(gz)}, max={np.max(gz)}") + if args.save_gz: + np.save(args.save_gz, gz) + print(f"gz saved to {args.save_gz}") solver.finalize() diff --git a/pygeos-tools/examples/solvers/gravity/modeling/gravityModelingLinearOp.py b/pygeos-tools/examples/solvers/gravity/modeling/gravityModelingLinearOp.py index 9995a05a..cb47f7da 100644 --- a/pygeos-tools/examples/solvers/gravity/modeling/gravityModelingLinearOp.py +++ b/pygeos-tools/examples/solvers/gravity/modeling/gravityModelingLinearOp.py @@ -21,20 +21,31 @@ from pygeosx import run, COMPLETED __doc__ = """ -“This is an example of how to set up and run your GEOS simulation using the GravitySolver. +Example script demonstrating how to configure and run a GEOS gravity simulation using the GravityLinearOpSolver. """ def parse_args(): - """Get arguments + """ + Parse command-line arguments for the simulation. Returns: - argument '--xml': Input xml file for GEOSX + argparse.Namespace: Parsed arguments including: + --xml (str): Path to the GEOS input XML file (required). + --model (str, optional): Path to a .npy file containing the density model. + --nm (int, optional): Number of model elements. Required if --model is not provided. + --save_gz (str, optional): Path to save the computed gz output as a .npy file. """ - parser = argparse.ArgumentParser( description="Gravity simulation example" ) - parser.add_argument( '--xml', type=str, required=True, help="Input xml file for GEOS" ) - parser.add_argument( '--m_true', type=str, required=True, help="True model (.npy)" ) + parser = argparse.ArgumentParser(description="Gravity modeling example") + parser.add_argument('--xml', type=str, required=True, help="Input xml file for GEOS") + parser.add_argument('--model', type=str, default=None, + help="Density model file (.npy). If not provided, --nm must be specified") + parser.add_argument('--nm', type=int, default=None, + help="Model size. Required if --model is not provided") + parser.add_argument("--save_gz", type=str, default=None, + help="Optional output file to save gz (.npy)") + args, _ = parser.parse_known_args() return args @@ -47,14 +58,30 @@ def main(): xmlfile = args.xml xml = XML( xmlfile ) - m_true = np.load( args.m_true ) - - solver = GravityLinearOpSolver( rank=rank, xml=xml, nm=m_true.size ) - - print( f"Density min={np.min(m_true)}, max={np.max(m_true)}", flush=True ) - - gz = solver.getData( m_true ) - print( f"gz min={np.min(gz)}, max={np.max(gz)}" ) + if args.model is not None: + model = np.load(args.model) + print( f"Density min={np.min(model)}, max={np.max(model)}", flush=True ) + if args.nm is not None and model.size != args.nm: + raise ValueError(f"Mismatch: model has size {model.size}, but nm={args.nm}") + nm = model.size + else: + if args.nm is None: + raise ValueError("Either --model or --nm must be provided") + model = None + nm = args.nm + + # Initialize solver + solver = GravityLinearOpSolver( rank=rank, xml=xml, nm=nm ) + + # Perform modeling + gz = solver.getData( model ) + + # Save output + if rank == 0: + print(f"gz min={np.min(gz)}, max={np.max(gz)}") + if args.save_gz: + np.save(args.save_gz, gz) + print(f"gz saved to {args.save_gz}") solver.finalize() diff --git a/pygeos-tools/src/geos/pygeos_tools/solvers/GravitySolver.py b/pygeos-tools/src/geos/pygeos_tools/solvers/GravitySolver.py index f4b59aee..48095176 100644 --- a/pygeos-tools/src/geos/pygeos_tools/solvers/GravitySolver.py +++ b/pygeos-tools/src/geos/pygeos_tools/solvers/GravitySolver.py @@ -55,7 +55,7 @@ def getGzAtStations( self: Self ) -> npt.NDArray: """ return self.getGeosWrapperByName( "gzAtStations", [ "Solvers" ] ) - def getDensityModel( self: Self, filterGhost: bool = False, **kwargs ) -> Optional[np.ndarray]: + def getDensityModel( self: Self, filterGhost: bool = False, **kwargs ) -> Optional[ np.ndarray ]: """ Get the density values WARNING: this function aims to work in the specific case of having only 1 CellElementRegion in your XML file @@ -103,14 +103,15 @@ def setMode( self: Self, mode: str ) -> None: mode_key = get_matching_wrapper_path( self.geosx, [ 'mode' ] ) self.geosx.get_wrapper( mode_key ).set_value( mode ) - def modeling( self: Self, model: npt.NDArray, scale_data: float = 1.0 ) -> npt.NDArray: + def modeling( self: Self, model: Optional[ npt.NDArray ], scale_data: float = 1.0 ) -> npt.NDArray: # Make sure we are in modeling mode. self.setMode( "modeling" ) self.applyInitialConditions() # Send model. - self.updateDensityModel( model ) + if model is not None: + self.updateDensityModel( model ) # Run. while run() != COMPLETED: diff --git a/pygeos-tools/src/geos/pygeos_tools/solvers/InversionUtils.py b/pygeos-tools/src/geos/pygeos_tools/solvers/InversionUtils.py index 9b034b02..61898c4c 100644 --- a/pygeos-tools/src/geos/pygeos_tools/solvers/InversionUtils.py +++ b/pygeos-tools/src/geos/pygeos_tools/solvers/InversionUtils.py @@ -58,7 +58,7 @@ def computeL2Loss( predicted: np.ndarray, observed: np.ndarray, scale: float = 1 residual = computeResidual( predicted, observed ) loss = 0.5 * np.linalg.norm( residual )**2 - return float(scale * loss) + return float( scale * loss ) def gradientTest( f: Callable[ [ np.ndarray ], float ], @@ -98,7 +98,6 @@ def gradientTest( f: Callable[ [ np.ndarray ], float ], Slope of the error curve at every step (should be 2. if the jacobian is correct). history: :obj:`np.ndarray` Full convergence history (steps and errors), useful for further display. - """ f0 = f( m0 ) g0_dm = np.dot( g( m0 ), dm ) @@ -128,8 +127,8 @@ def gradientTest( f: Callable[ [ np.ndarray ], float ], if ind.size > 1: slope = np.diff( np.log( e1[ ind ] ) ) / np.diff( np.log( h[ ind ] ) ) mean_slope = np.mean( slope ) - passed = bool( isclose( mean_slope, 2.0, abs_tol=0.1 ) or np.count_nonzero( slope > 1.9 ) > 2 - or np.all( e1 < 1e-15 ) ) + passed = bool( + isclose( mean_slope, 2.0, abs_tol=0.1 ) or np.count_nonzero( slope > 1.9 ) > 2 or np.all( e1 < 1e-15 ) ) else: slope = np.array( [] ) passed = False