diff --git a/docs/source/_static/custom.css b/docs/source/_static/custom.css new file mode 100644 index 0000000..61ba797 --- /dev/null +++ b/docs/source/_static/custom.css @@ -0,0 +1,3 @@ +img.transparent-logo { + background-color: transparent !important; +} \ No newline at end of file diff --git a/docs/source/_static/profsea-logo-dark.png b/docs/source/_static/profsea-logo-dark.png new file mode 100644 index 0000000..a355518 Binary files /dev/null and b/docs/source/_static/profsea-logo-dark.png differ diff --git a/docs/source/_static/profsea-logo.png b/docs/source/_static/profsea-logo-light.png similarity index 100% rename from docs/source/_static/profsea-logo.png rename to docs/source/_static/profsea-logo-light.png diff --git a/docs/source/conf.py b/docs/source/conf.py index 3109f1b..b776318 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -79,10 +79,17 @@ "type": "fontawesome", } ], + "logo": { + "image_light": "_static/profsea-logo-light.png", + "image_dark": "_static/profsea-logo-dark.png", + }, } html_sidebars = { + "index": [], # This removes the left sidebar on the landing page "**": ["search-field", "sidebar-nav-bs", "page-toc"], } html_static_path = ["_static"] -html_logo = "_static/profsea-logo.png" html_favicon = "_static/favicon.png" +html_css_files = [ + "custom.css", +] diff --git a/docs/source/index.rst b/docs/source/index.rst index 07ca9ff..6afa9fa 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -1,8 +1,64 @@ -ProFSea -======= +.. title:: ProFSea -**ProFSea is sea-level rise simulator based on statistical emulations of physical modelling experiments and lines of evidence from the IPCC and across the literature. It's modular, easy to setup and self-contained - all you need is global mean surface temperature and ocean heat content forcing anomalies, using any baseline period, and you're good to go.** +.. grid:: 1 1 2 2 + :margin: 4 4 0 0 + :padding: 0 + + .. grid-item:: + :columns: 12 12 8 8 + + .. rst-class:: display-2 font-weight-bold + + ProFSea + + .. rst-class:: lead + + **🌊 A modular, easy-to-setup, and self-contained sea-level rise simulator based on statistical emulations of physical modelling experiments and lines of evidence from the IPCC.** + + .. container:: d-flex gap-3 pt-3 + + .. button-ref:: user_guide + :ref-type: doc + :color: primary + :shadow: + :class: font-weight-bold + + Get Started + + .. button-link:: [https://github.com/MetOffice/ProFSea-tool](https://github.com/MetOffice/ProFSea-tool) + :color: secondary + :shadow: + :outline: + :class: font-weight-bold + + View on GitHub + + .. grid-item:: + :columns: 12 12 4 4 + + .. image:: /_static/logo.png + :alt: ProFSea Logo + :class: align-center transparent-logo + :width: 100% + + +---- +ProFSea makes complex sea-level rise simulations accessible and fast. All you need is global mean surface temperature and ocean heat content forcing anomalies, using any baseline period, and you're good to go. ProFSea development is supported by the MetOffice. + +Quick Install +------------- + +ProFSea is available as a Python package. To install it: + +.. code-block:: bash + + pip install profsea + +.. note:: + ProFSea relies on standard scientific libraries including ``numpy``, ``xarray``, and ``dask``. Check the :doc:`user_guide` for detailed dependency requirements. + +---- .. grid:: 1 1 2 2 :gutter: 2 @@ -35,6 +91,18 @@ ProFSea Academic and software references for the project. +Getting in Touch +---------------- + +Whether you need help getting started with ProFSea, found a bug, want the emulator to be more awesome, or just want to chat about sea-level modelling, you have a few options: + +* **Open an issue** on the `GitHub repository `_ if you find a bug, need a missing feature, or spot a typo in this documentation 👀. +* **Start a discussion** on GitHub for general questions about emulations, statistical methods, or setting up a new experiment. + +Citing ProFSea +-------------- + +If you use ProFSea for your research, teaching, or analysis, please credit the project by citing the relevant academic references. See the :doc:`references` page for the full bibliography, methodology papers, and software DOIs. .. toctree:: :caption: User Guide @@ -47,6 +115,8 @@ ProFSea :caption: Development :maxdepth: 1 :hidden: + + documentation .. toctree:: :caption: Reference diff --git a/profsea/components/core/global_model.py b/profsea/components/core/global_model.py index aacd519..1c44d29 100644 --- a/profsea/components/core/global_model.py +++ b/profsea/components/core/global_model.py @@ -2,18 +2,19 @@ import concurrent.futures import logging -import os -from pathlib import Path import numpy as np import xarray as xr +from rich.console import Console +from rich.progress import track -from profsea.utils import check_shapes, sample_members_2D +from profsea.utils import check_shapes, sample_members_2D, save_components from profsea.utils.ui import print_global_preflight from .base import Component from .state import ClimateState +console = Console() logger = logging.getLogger(__name__) @@ -68,6 +69,7 @@ def __init__( output_percentiles: list | np.ndarray = None, palmer_method: bool = True, random_sample: bool = False, + dtype: np.dtype | str = np.float32, ): self.components = components self.end_yr = end_yr @@ -78,11 +80,15 @@ def __init__( self.output_percentiles = output_percentiles self.palmer_method = palmer_method self.random_sample = random_sample + self.dtype = np.dtype(dtype) self.endofhistory = 2006 self.endofAR5 = 2100 self.nyr = self.end_yr - self.endofhistory + # Inject method! + save_components = save_components + def _arr_to_xr(self, arr_dict: dict[str, np.ndarray]) -> dict[str, xr.DataArray]: """Convert a dictionary of numpy/dask arrays to xarray DataArrays. @@ -152,10 +158,11 @@ def run( print_global_preflight(self, scenario) + T_change = T_change.astype(self.dtype) T_ens, T_int_ens, T_int_med = self._calculate_drivers(T_change) # Shared physical correlation state - fraction = run_rng.random(self.num_members * self.nt) + fraction = run_rng.random(self.num_members * self.nt).astype(self.dtype) state = ClimateState( scenario=scenario, @@ -170,6 +177,7 @@ def run( nyr=self.nyr, nt=self.nt, num_members=self.num_members, + dtype=self.dtype, ) # Child RNGs for each component @@ -193,7 +201,9 @@ def run( except Exception as e: raise RuntimeError(f"Component '{comp_name}' failed.") from e else: - for name, comp in self.components.items(): + for name, comp in track( + self.components.items(), description="Projecting components..." + ): results[name] = comp.project(state, comp_rngs[name]) # Random Sampling @@ -209,42 +219,138 @@ def run( f"Sampling {len(self.output_percentiles)} members per component..." ) for comp_name, data in results.items(): - results[comp_name] = sample_members_2D(data, self.output_percentiles) + results[comp_name] = sample_members_2D( + data, self.output_percentiles, dtype=self.dtype + ) self.results = self._arr_to_xr(results) return self.results - def sum_components(self, components: dict[str, xr.DataArray]) -> xr.DataArray: - """Sum the components to get total GMSLR.""" - gmslr = xr.concat( - [components[name] for name in components.keys()], dim="component" - ).sum(dim="component") - gmslr.attrs["units"] = "m" - gmslr.attrs["description"] = "Total global mean sea level rise" - components["total_gmslr"] = gmslr - return gmslr + # def save_components( + # self, + # components: dict[str, xr.DataArray], + # scenario_name: str, + # output_prefix: str = "global", + # output_dir: str = ".", + # output_format: str = "netcdf", + # ) -> None: + # """ + # Stream all global sea level projections to disk in a single file/store. + + # Parameters + # ---------- + # components: dict[str, xr.DataArray] + # Dictionary of component names and their corresponding Xarray DataArrays. + # scenario_name: str + # Name of the scenario you've run the emulator for. + # output_prefix: str + # Prefix for the output file name (default: 'global'). + # output_dir: str + # Directory to save components to. + # output_format: str + # Format to save the output in. Must be either 'netcdf' or 'zarr'. + + # Returns + # ------- + # None + # """ + # ds = xr.Dataset(components) + + # # Add ProFSea version and scenario metadata + # ds.attrs["source"] = "ProFSea v3.0" + # ds.attrs["scenario"] = scenario_name + # ds.attrs["description"] = "Global sea level rise projections" + + # output_format = output_format.lower() + # if output_format not in ["netcdf", "zarr"]: + # raise ValueError("output_format must be either 'netcdf' or 'zarr'.") + + # Path(output_dir).mkdir(parents=True, exist_ok=True) + # encoding = {} + + # # Sort out Zarr encoding + # if output_format == "zarr": + # import numcodecs + # from numcodecs.zarr3 import Blosc + + # compressor = Blosc( + # cname="zstd", clevel=5, shuffle=numcodecs.Blosc.BITSHUFFLE + # ) + + # # Set the encoding/compression dynamically based on the component's actual dtype + # for name, component in components.items(): + # comp_dtype = ( + # component.dtype.name + # ) # Captures 'float32' or 'float64' dynamically + + # if output_format == "netcdf": + # encoding[name] = {"zlib": True, "complevel": 1, "dtype": comp_dtype} + # elif output_format == "zarr": + # encoding[name] = {"compressor": compressor, "dtype": comp_dtype} + + # file_name = f"{scenario_name}_{output_prefix}" + + # # Stream the computation and write to disk + # if output_format == "netcdf": + # out_path = os.path.join(output_dir, f"{file_name}.nc") + # with console.status( + # "[bold cyan]Computing and saving Global NetCDF...[/bold cyan]", + # spinner="dots", + # ): + # # If using Dask, .compute() is required before .to_netcdf() + # # If arrays are already eager NumPy arrays, .compute() is a harmless no-op + # if hasattr(ds, "compute"): + # ds.compute().to_netcdf(out_path, encoding=encoding) + # else: + # ds.to_netcdf(out_path, encoding=encoding) + + # logger.info( + # f"[bold green]✓ Successfully saved NetCDF:[/bold green] {out_path}" + # ) + + # elif output_format == "zarr": + # out_path = os.path.join(output_dir, f"{file_name}.zarr") + # with console.status( + # "[bold cyan]Streaming computation and saving Global Zarr...[/bold cyan]", + # spinner="dots", + # ): + # ds.to_zarr(out_path, encoding=encoding, mode="w", compute=True) + # logger.info( + # f"[bold green]✓ Successfully saved Zarr:[/bold green] {out_path}" + # ) + + # # Log the shape of the total_gmslr (or the first available component) + # sample_name = ( + # "total_gmslr" if "total_gmslr" in ds else list(ds.data_vars.keys())[0] + # ) + # dims_str = ", ".join(ds[sample_name].dims) + # logger.info(f"Global output shape was {ds[sample_name].shape} ({dims_str})") - def save_components( - self, components: dict[str, xr.DataArray], output_dir: str, scenario_name: str - ) -> None: - """Save SLR components as nc files to a directory. + def sum_components(self, components: dict[str, xr.DataArray]) -> xr.DataArray: + """ + Sum the components in-place to get total GMSLR. Parameters ---------- - components: Dict[str, xr.DataArray] - Dictionary of component names and their corresponding xarray DataArrays. - output_dir: str - Directory to save components to. - scenario_name: str - Name of the scenario you've run the emulator for. + components: dict[str, xr.DataArray] + Dictionary of component names and their corresponding Xarray DataArrays. Returns ------- - None + xr.DataArray + DataArray of the summed global projections. """ - Path(output_dir).mkdir(parents=True, exist_ok=True) - ds = xr.Dataset(components) - ds.to_netcdf(os.path.join(output_dir, f"{scenario_name}_global.nc")) + + iterator = iter(components.values()) + gmslr = next(iterator).copy() + + for comp in iterator: + gmslr += comp + + gmslr.attrs["units"] = "m" + gmslr.attrs["description"] = "Total global mean sea level rise" + components["total_gmslr"] = gmslr + return gmslr def _calculate_drivers(self, T_change: np.ndarray) -> tuple: """Calculate the drivers of GMSLR: temperature change and @@ -264,4 +370,4 @@ def _calculate_drivers(self, T_change: np.ndarray) -> tuple: # Time-integral of temperature anomaly T_int_ens = np.cumsum(T_ens, axis=1) T_int_med = np.cumsum(np.median(T_ens, axis=0)) - return T_ens, T_int_ens, T_int_med + return (T_ens, T_int_ens, T_int_med) diff --git a/profsea/components/core/local_model.py b/profsea/components/core/local_model.py index 1efb777..d3ce230 100644 --- a/profsea/components/core/local_model.py +++ b/profsea/components/core/local_model.py @@ -238,10 +238,6 @@ def sum_components(self, components: dict[str, xr.DataArray]) -> xr.DataArray: total_rsl = xr.concat(components.values(), dim="component").sum( dim="component", skipna=False ) - total_rsl.attrs = { - "units": "m", - "long_name": "Local total sea-level projections", - "source": "ProFSea-Climate v0.1", - } + total_rsl.attrs["units"] = "m" components["total_rsl"] = total_rsl return total_rsl diff --git a/profsea/components/core/spatial_model.py b/profsea/components/core/spatial_model.py index 6c22662..1303754 100644 --- a/profsea/components/core/spatial_model.py +++ b/profsea/components/core/spatial_model.py @@ -37,6 +37,7 @@ def __init__( end_year: int = 2301, baseline_yrs: tuple = (1995, 2014), output_percentiles: list | np.ndarray = [5, 17, 50, 83, 95], + dtype: np.dtype = np.float32, ) -> None: """ Parameters @@ -54,6 +55,8 @@ def __init__( Tuple defining the start and end years of the baseline period for calculating anomalies. Default is (1995, 2014). output_percentiles: list or np.ndarray, optional List or array of percentiles to sample from the ensemble for output. If None, outputs all members. Default is [5, 17, 50, 83, 95]. + dtype: np.dtype, optional + Data type for the output arrays. Default is np.float32. """ # Define the path where the data should live @@ -69,6 +72,7 @@ def __init__( self.output_percentiles = output_percentiles self.start_year = 2006 self.n_years = self.end_year - self.start_year + self.dtype = dtype if self.output_percentiles is not None and len(self.output_percentiles) > 0: self.num_members = len(self.output_percentiles) @@ -196,6 +200,7 @@ def run(self, member_seed: int = 42) -> None: grid_interpolation="linear", output_percentiles=self.output_percentiles, baseline_yrs=self.baseline_yrs, + dtype=self.dtype, ) child_seeds = seed_seq.spawn(len(self.components)) @@ -240,13 +245,6 @@ def sum_components(self, components: dict[str, xr.DataArray]) -> xr.DataArray: total_rsl = xr.concat(components.values(), dim="component").sum( dim="component", skipna=False ) - - # Optionally, apply attributes so it matches the other DataArrays - total_rsl.attrs = { - "units": "m", - "long_name": "Regional total sea-level projections", - "source": "ProFSea-Climate v0.1", - } - + total_rsl.attrs["units"] = "m" components["total_rsl"] = total_rsl return total_rsl diff --git a/profsea/components/core/state.py b/profsea/components/core/state.py index 7fc7296..442c731 100644 --- a/profsea/components/core/state.py +++ b/profsea/components/core/state.py @@ -21,6 +21,7 @@ class ClimateState: nyr: int nt: int num_members: int + dtype: np.dtype = np.float32 @dataclass @@ -34,6 +35,7 @@ class SpatialState: grid_interpolation: str output_percentiles: list[int] | np.ndarray baseline_yrs: tuple[int, int] + dtype: np.dtype = np.float32 @dataclass @@ -47,3 +49,4 @@ class LocalState: interpolation_method: str output_percentiles: list[int] | np.ndarray baseline_yrs: tuple[int, int] + dtype: np.dtype = np.float32 diff --git a/profsea/components/core/time_projection.py b/profsea/components/core/time_projection.py index bb38d1c..b1372e8 100644 --- a/profsea/components/core/time_projection.py +++ b/profsea/components/core/time_projection.py @@ -38,10 +38,10 @@ def time_projection( """ # Create a field of elapsed time since start in years timeendofAR5 = state.endofAR5 - state.endofhistory + 1 - time = np.arange(state.end_yr - state.endofhistory) + 1 + time = np.arange(state.end_yr - state.endofhistory, dtype=state.dtype) + 1 if fraction is None: - fraction = rng.random((state.num_members, state.nt)) + fraction = rng.random((state.num_members, state.nt), dtype=state.dtype) elif fraction.size != state.num_members * state.nt: raise ValueError("fraction is the wrong size") @@ -49,9 +49,9 @@ def time_projection( # Convert inputs to startrate (m yr-1) and afinal (m), where both are # arrays with the size of fraction - startrate = ( - startratemean + startratepm * np.array([-1, 1], dtype=float) - ) * 1e-3 # convert mm yr-1 to m yr-1 + startrate = ((startratemean + startratepm * np.array([-1, 1])) * 1e-3).astype( + state.dtype + ) # convert mm yr-1 to m yr-1 finalisrange = isinstance(final, Sequence) if finalisrange: @@ -71,7 +71,9 @@ def time_projection( # a = S/t**2-b/t = (S-b*t)/t**2 # If nfinal=1, the following two lines are equivalent to # halfacc=(final-startyr*nyr)/nyr**2 - finalyr = np.arange(nfinal) - nfinal + 94 + 1 # last element ==nyr + finalyr = ( + np.arange(nfinal, dtype=state.dtype) - nfinal + 94 + 1 + ) # last element ==nyr halfacc = (afinal - startrate * finalyr.mean()) / (finalyr**2).mean() quadratic = halfacc[:, :, np.newaxis] * (time**2) linear = startrate[:, :, np.newaxis] * time diff --git a/profsea/components/global_/antarctica.py b/profsea/components/global_/antarctica.py index 453e224..4b77270 100644 --- a/profsea/components/global_/antarctica.py +++ b/profsea/components/global_/antarctica.py @@ -60,6 +60,7 @@ def _impulse_response_term( gamma: float, params_1d: np.ndarray, dt: float, + state: ClimateState, ) -> np.ndarray: """Computes the slow response term for a single trajectory.""" n_time = tas_1d.shape[0] @@ -68,7 +69,7 @@ def _impulse_response_term( forcing_base = np.sign(tas_1d) * (np.abs(tas_1d) ** gamma) # decay_factors1 = np.exp(-np.arange(n_time) * dt / tau1) * (dt / tau1) - t_arr = np.arange(n_time) + t_arr = np.arange(n_time, dtype=state.dtype) decay_factors1 = (t_arr * dt / tau1**2) * np.exp(-t_arr * dt / tau1) * dt rate_delayed1 = fftconvolve(forcing_base, decay_factors1, mode="full")[:n_time] term_slow1 = alpha1 * (np.cumsum(rate_delayed1) * dt) @@ -81,21 +82,21 @@ def _impulse_response_term( return term_slow1 + term_slow2 def _precompute_delayed_rates( - self, tas: np.ndarray, dt: float + self, tas: np.ndarray, dt: float, state: ClimateState ) -> tuple[np.ndarray, np.ndarray]: """ Precomputes the cumulative delayed rates for all models across all trajectories. Returns two arrays of shape (n_models, n_traj, n_time). """ n_traj, n_time = tas.shape - cum_rate1 = np.zeros((self.n_models, n_traj, n_time)) - cum_rate2 = np.zeros((self.n_models, n_traj, n_time)) - t_arr = np.arange(n_time) + cum_rate1 = np.zeros((self.n_models, n_traj, n_time), dtype=state.dtype) + cum_rate2 = np.zeros((self.n_models, n_traj, n_time), dtype=state.dtype) + t_arr = np.arange(n_time, dtype=state.dtype) for m_idx in range(self.n_models): - tau1 = float(self.param_ds.tau1[m_idx].values) - tau2 = float(self.param_ds.tau2[m_idx].values) - gamma = float(self.param_ds.gamma[m_idx].values) + tau1 = np.array(self.param_ds.tau1[m_idx].values, dtype=state.dtype) + tau2 = np.array(self.param_ds.tau2[m_idx].values, dtype=state.dtype) + gamma = np.array(self.param_ds.gamma[m_idx].values, dtype=state.dtype) # Forcing base is model-specific due to gamma. Shape: (n_traj, n_time) forcing_base = np.sign(tas) * (np.abs(tas) ** gamma) @@ -132,25 +133,29 @@ def project(self, state: ClimateState, rng: np.random.Generator) -> np.ndarray: nm = state.nt * state.num_members # 1. Precompute expensive convolutions for unique (Model x Trajectory) combos - cum_rate1, cum_rate2 = self._precompute_delayed_rates(tas, dt) + cum_rate1, cum_rate2 = self._precompute_delayed_rates(tas, dt, state) tas_int = np.cumsum(tas, axis=1) * dt # 2. Generate random model/residual assignments for all members model_indices = rng.integers(0, self.n_models, size=nm) - all_residuals = self.param_ds.param_residuals.values + all_residuals = np.asarray( + self.param_ds.param_residuals.values, dtype=state.dtype + ) n_train_scenarios = all_residuals.shape[1] residual_indices = rng.integers(0, n_train_scenarios, size=nm) # 3. Create flat mapping array to match the (Trajectory x Member) layout # This groups by trajectory: [Traj0_Mem0...Traj0_Mem999, Traj1_Mem0...] - t_indices = np.repeat(np.arange(n_traj), state.num_members) + t_indices = np.repeat(np.arange(n_traj, dtype=np.int32), state.num_members) # NOTE: If your required grouping is interleaved [Traj0_Mem0, Traj1_Mem0...] # uncomment the line below instead: - # t_indices = np.tile(np.arange(n_traj), state.num_members) + # t_indices = np.tile(np.arange(n_traj, dtype=state.dtype), state.num_members) # 4. Vectorized Parameter Extraction - general_p = self.param_ds.general_params.values[model_indices] + general_p = np.asarray( + self.param_ds.general_params.values[model_indices], dtype=state.dtype + ) sampled_residuals = all_residuals[model_indices, residual_indices, :] total_params = general_p + sampled_residuals @@ -168,58 +173,6 @@ def project(self, state: ClimateState, rng: np.random.Generator) -> np.ndarray: return term_slow + term_fast - # def project(self, state: ClimateState, rng: np.random.Generator) -> np.ndarray: - # """ - # Projects AIS response using empirical additive bootstrapping aligned - # to the ClimateState ensemble size. - # """ - # tas = state.T_ens - # if tas.ndim > 2: - # tas = np.squeeze(tas) - # if tas.ndim == 1: - # tas = np.expand_dims(tas, axis=0) - - # dt = 1.0 - # nm = state.nt * state.num_members - # n_time = tas.shape[1] - - # preds = np.zeros((nm, n_time)) - # tas_int = np.cumsum(tas, axis=1) * dt - - # # Randomly assign an ISMIP6 model and residual draw to each ensemble member - # model_indices = rng.integers(0, self.n_models, size=nm) - # all_residuals = self.param_ds.param_residuals.values - # n_train_scenarios = all_residuals.shape[1] - # residual_indices = rng.integers(0, n_train_scenarios, size=nm) - - # for i in range(nm): - # t_idx = i // state.num_members - # m_idx = model_indices[i] - # r_idx = residual_indices[i] - - # # Extract assigned model parameters - # tau1 = float(self.param_ds.tau1[m_idx].values) - # tau2 = float(self.param_ds.tau2[m_idx].values) - # gamma = float(self.param_ds.gamma[m_idx].values) - - # general_p = self.param_ds.general_params[m_idx].values - # sampled_residuals = all_residuals[m_idx, r_idx, :] - # total_params = general_p + sampled_residuals - - # # Slow response - # term_slow = self._impulse_response_term( - # tas[t_idx], tau1, tau2, gamma, total_params, dt - # ) - - # # Fast response - # beta = total_params[2] - # term_fast = beta * tas_int[t_idx] - - # # Combine - # preds[i, :] = term_fast + term_slow - - # return preds - class AntarcticaDynAR5(Component): """ @@ -266,11 +219,18 @@ def project(self, state: ClimateState, rng: np.random.Generator) -> np.ndarray: ) lcoeff = lcoeff[state.scenario] - ascale = norm.ppf(state.fraction) + ascale = norm.ppf(state.fraction).astype(state.dtype) final = np.exp(lcoeff[2] * ascale**2 + lcoeff[1] * ascale + lcoeff[0]) final = final.reshape(state.num_members, state.nt) return ( - time_projection(state, 0.41, 0.20, final, rng, fraction=state.fraction) + time_projection( + state, + 0.41, + 0.20, + final, + rng, + fraction=state.fraction, + ) + self.d_ant ) @@ -285,8 +245,7 @@ class AntarcticaSMBAR5(Component): """ def __init__(self): - # Conversion factor for Gt to m SLE - self.mSLEoGt = 1e12 / 3.61e14 * 1e-3 + pass def project( self, @@ -307,21 +266,29 @@ def project( antsmb: np.ndarray Antarctic SMB contribution to GMSLR. """ + # Conversion factor for Gt to m SLE + mSLEoGt = 1e12 / 3.61e14 * 1e-3 # The following are [mean,SD] pcoK = [5.1, 1.5] # % change in Ant SMB per K of warming from G&H06 KoKg = [1.1, 0.2] # ratio of Antarctic warming to global warming from G&H06 # Generate a distribution of products of the above two factors pcoKg = ( - pcoK[0] + rng.standard_normal([state.num_members, state.nt]) * pcoK[1] - ) * (KoKg[0] + rng.standard_normal([state.num_members, state.nt]) * KoKg[1]) + pcoK[0] + + rng.standard_normal([state.num_members, state.nt], dtype=state.dtype) + * pcoK[1] + ) * ( + KoKg[0] + + rng.standard_normal([state.num_members, state.nt], dtype=state.dtype) + * KoKg[1] + ) meansmb = 1923 # model-mean time-mean 1979-2010 Gt yr-1 from 13.3.3.2 moaoKg = ( - -pcoKg * 1e-2 * meansmb * self.mSLEoGt + -pcoKg * 1e-2 * meansmb * mSLEoGt ) # m yr-1 of SLE per K of global warming if state.fraction is None: - fraction = rng.random((state.num_members, state.nt)) + fraction = rng.random((state.num_members, state.nt), dtype=state.dtype) elif state.fraction.size != state.num_members * state.nt: raise ValueError("fraction is the wrong size") else: diff --git a/profsea/components/global_/expansion.py b/profsea/components/global_/expansion.py index 1bf5b05..b6945ff 100644 --- a/profsea/components/global_/expansion.py +++ b/profsea/components/global_/expansion.py @@ -22,6 +22,8 @@ def __init__(self, OHC_change: np.ndarray, distribution_scaler: float = 1.0): self.distribution_scaler = distribution_scaler def project(self, state: ClimateState, rng: np.random.Generator) -> np.ndarray: + self.OHC_change = np.asarray(self.OHC_change, dtype=state.dtype) + # check the shape here check_shapes(self.OHC_change, state.nyr) @@ -39,7 +41,7 @@ def project(self, state: ClimateState, rng: np.random.Generator) -> np.ndarray: exp_efficiency = ( rng.normal(loc=mean_eff, scale=std_eff, size=(state.nt, state.num_members)) * 1e-24 - ) # m/YJ + ).astype(state.dtype) # m/YJ ohc_3d = self.OHC_change[:, None, :] # Efficiency shape: (nt, num_members, 1) diff --git a/profsea/components/global_/glacier.py b/profsea/components/global_/glacier.py index e525044..4f159f9 100644 --- a/profsea/components/global_/glacier.py +++ b/profsea/components/global_/glacier.py @@ -66,24 +66,26 @@ def project(self, state: ClimateState, rng: np.random.Generator) -> np.ndarray: ngl = len(glparm) model_indices = rng.integers(0, ngl, size=(state.nt, state.num_members)) - base_factors = np.array([p["factor"] for p in glparm]) - base_exponents = np.array([p["exponent"] for p in glparm]) - base_cvgls = np.array([p["cvgl"] for p in glparm]) + base_factors = np.array([p["factor"] for p in glparm], dtype=state.dtype) + base_exponents = np.array([p["exponent"] for p in glparm], dtype=state.dtype) + base_cvgls = np.array([p["cvgl"] for p in glparm], dtype=state.dtype) factors = base_factors[model_indices][:, :, None] exponents = base_exponents[model_indices][:, :, None] cvgls = base_cvgls[model_indices][:, :, None] - r = rng.standard_normal((state.nt, state.num_members))[:, :, None] + r = rng.standard_normal((state.nt, state.num_members), dtype=state.dtype)[ + :, :, None + ] T_int_ens_3d = state.T_int_ens[:, None, :] # Median shape: (1, 1, nyr) T_int_med_3d = state.T_int_med[None, None, :] - zgl = self._project_glacier1(T_int_ens_3d, factors, exponents) + zgl = self._project_glacier1(T_int_ens_3d, factors, exponents, state) # Passes (1, 1, nyr) + (nt, nm, 1) -> Returns (nt, nm, nyr) - mgl = self._project_glacier1(T_int_med_3d, factors, exponents) + mgl = self._project_glacier1(T_int_med_3d, factors, exponents, state) # 6. Apply variance and clip using 3D matrix math glacier = zgl + (mgl * r * cvgls) @@ -94,7 +96,11 @@ def project(self, state: ClimateState, rng: np.random.Generator) -> np.ndarray: return glacier.reshape(state.nt * state.num_members, state.nyr) def _project_glacier1( - self, T_int: np.ndarray, factor: np.ndarray, exponent: np.ndarray + self, + T_int: np.ndarray, + factor: np.ndarray, + exponent: np.ndarray, + state: ClimateState, ) -> np.ndarray: """Project glacier contribution by one glacier method.""" scale = 1e-3 # mm to m diff --git a/profsea/components/global_/greenland.py b/profsea/components/global_/greenland.py index 8910dac..0ecc215 100644 --- a/profsea/components/global_/greenland.py +++ b/profsea/components/global_/greenland.py @@ -40,7 +40,7 @@ def project(self, state: ClimateState, rng: np.random.Generator) -> np.ndarray: tas = np.expand_dims(tas, axis=0) nt = tas.shape[0] - time_delta = np.arange(state.nyr) + time_delta = np.arange(state.nyr, dtype=state.dtype) df = self.df n_models = len(df) @@ -48,12 +48,12 @@ def project(self, state: ClimateState, rng: np.random.Generator) -> np.ndarray: model_indices = rng.integers(0, n_models, size=(nt, state.num_members)) # Extract parameters and reshape to 3D: (nt, nm, 1) - b0 = df["b0"].values[model_indices][:, :, None] - b1 = df["b1"].values[model_indices][:, :, None] - b2 = df["b2"].values[model_indices][:, :, None] - b3 = df["b3"].values[model_indices][:, :, None] - b4 = df["b4"].values[model_indices][:, :, None] - b5 = df["b5"].values[model_indices][:, :, None] + b0 = df["b0"].values[model_indices][:, :, None].astype(state.dtype) + b1 = df["b1"].values[model_indices][:, :, None].astype(state.dtype) + b2 = df["b2"].values[model_indices][:, :, None].astype(state.dtype) + b3 = df["b3"].values[model_indices][:, :, None].astype(state.dtype) + b4 = df["b4"].values[model_indices][:, :, None].astype(state.dtype) + b5 = df["b5"].values[model_indices][:, :, None].astype(state.dtype) # GIS trend values taken from FACTS GitHub repo trend_mean = 0.19 @@ -68,7 +68,8 @@ def project(self, state: ClimateState, rng: np.random.Generator) -> np.ndarray: b=b_bound, loc=trend_mean, scale=trend_std, - ) + ).astype(state.dtype) + trend_sle = (trend[:, :, None] * time_delta[None, None, :]) * 1e-3 tas_3d = tas[:, None, :] @@ -131,9 +132,12 @@ def project(self, state: ClimateState, rng: np.random.Generator) -> np.ndarray: febound = [1, 1.15] # bounds of uniform pdf of SMB elevation feedback factor # random log-normal factor - fn = np.exp(rng.standard_normal(state.num_members) * fnlogsd) + fn = np.exp(rng.standard_normal(state.num_members, dtype=state.dtype) * fnlogsd) # elevation feedback factor - fe = rng.random(state.num_members) * (febound[1] - febound[0]) + febound[0] + fe = ( + rng.random(state.num_members, dtype=state.dtype) * (febound[1] - febound[0]) + + febound[0] + ) ff = fn * fe ztgreen = state.T_ens - dtgreen diff --git a/profsea/components/global_/landwater.py b/profsea/components/global_/landwater.py index d1f262b..cd434c2 100644 --- a/profsea/components/global_/landwater.py +++ b/profsea/components/global_/landwater.py @@ -36,14 +36,14 @@ def project(self, state: ClimateState, rng: np.random.Generator) -> np.ndarray: np.ndarray Land water storage contribution to GMSLR. """ - - lw_ds = self.lw_ds + self.lw_ds = self.lw_ds.astype(state.dtype) # Interpolate to annual projections - interp_ds = lw_ds.interp( - years=np.arange(2005, 2301, 1), method="linear" - ).squeeze() - interp_ds["sea_level_change"].values * 1e-3 # mm to m + interp_ds = ( + self.lw_ds.interp(years=np.arange(2005, 2301, 1), method="linear") + .squeeze() + .astype(state.dtype) + ) # Base array: Shape (n_samples_in_nc, 296) lw_base = interp_ds["sea_level_change"].values * 1e-3 diff --git a/profsea/scm_profsea_experiments.py b/profsea/scm_profsea_experiments.py index e62b592..5f1645e 100644 --- a/profsea/scm_profsea_experiments.py +++ b/profsea/scm_profsea_experiments.py @@ -382,7 +382,13 @@ def main(args): "glacier": Glacier(), } - global_model = Global(components=slr_components, end_yr=2301, num_members=1) + global_model = Global( + components=slr_components, + end_yr=2301, + num_members=1000, + dtype="float32", + parallel=False, + ) projections = global_model.run( scenario=scenario, T_change=tas_scen, diff --git a/profsea/tutorial_notebooks/worked_example.ipynb b/profsea/tutorial_notebooks/worked_example.ipynb index 1fd20ae..3194d14 100644 --- a/profsea/tutorial_notebooks/worked_example.ipynb +++ b/profsea/tutorial_notebooks/worked_example.ipynb @@ -79,7 +79,7 @@ "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAADZCAYAAADVJMz+AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAANf9JREFUeJztnQecFOX5x59tt9c5jt6rSFekozSp0gRJbCSI+tFgKCoxJug/QYMRTTSxx4CKIRARAwZsCChSBGmCIAKKVBGkHtdvb2/n//k9s7O3u7e7t7u3c7fl+foZ73Z25r2Xudv5zfM+zaAoikKCIAiC4AOjr52CIAiCICIhCIIgBEQsCUEQBMEvIhKCIAiCX0QkBEEQBL+ISAiCIAh+EZEQBEEQ/CIiIQiCIPjF7P+txMXhcNCPP/5IGRkZZDAYano6giAIVQZ503l5edS4cWMyGoO3D0QkfACBaNasWdV/K4IgCFHGyZMnqWnTpkEfLyLhA1gQ2sXMzMyM3G9HEAShhsjNzeWHX+3+FiwiEj7QlpggECISgiDEE6EuoYvjWhAEQfCLiIQgCIIQnSKxceNGGjt2LHvbYQL973//q+CNf+yxx/j9lJQUGjRoEO3fv7/ScZcvX04dO3Ykq9XKX999910d/xWCIAjxS436JAoKCuiqq66iO++8kyZOnFjh/b/85S/0t7/9jd58801q164dPfHEEzRs2DA6dOiQX+fL1q1b6ZZbbqG5c+fShAkTWCBuvvlm2rx5M/Xu3bsa/lWCEDrnc+x0Mbcs7EuXk6eem5Vhivjl13NsvceP17GzM01UN6t6bt+GaGk6BEsCN/Tx48fza0wLFsQDDzxAv/vd73hfSUkJNWjQgJ5++mn61a9+5XMcCAS8+B999JFr38iRI6l27dr01ltvBTUXnF+rVi26fPmyOK4F3Tl7sZQuOm8KDkfo5zsUhXLz1RNrpRsjmtuj59h6jx+vYxsNBqdImMhoDP7nhntfi1qfxNGjR+nMmTM0fPhw1z4sHw0cOJC2bNni9zxYEu7ngBEjRgQ8RxBiVSAEQW+iNgQWAgFgObiD18ePHw94nq9ztPF8AQsFm7viCoLey0vnL9tdr0UghGglai0JDW9TC8tQlZl2oZ4zb948NsO0TbKtBb2tB00gIA4iEEI0E7Ui0bBhQ/7qbQGcPXu2gqXgfV6o58yePZvX6bQNmdaCoAeyvCTEGlErEq1ateIb/tq1a137bDYbbdiwgfr16+f3vL59+3qcA9asWRPwHPg6tOxqybIW9EIEQohFatQnkZ+fT4cPH/ZwVu/Zs4eys7OpefPmHNn05JNP0hVXXMEbvk9NTaXbb7/ddc7kyZOpSZMmvGQE7r//fhowYABHQN144420cuVKWrduHYfACkJN4XAo7KBGaGNmmmdYo+LmsA4n1LCsrHzJym4nMkUwIlPPsfUeP17HNpjUvxOHUj1P+TUqEjt37qTBgwe7Xs+aNYu/3nHHHZwb8fDDD1NRURH9+te/pkuXLnGeA6wC9xyJEydOeJS9hcWwdOlS+r//+z/6wx/+QG3atKG3335bciSEGrcitNh3dx8ZBMLuTI/AB99eFrpM2EoV3kBhsYMslsjdOvQcW+/x43Vsq4Wo1K7wZjYZoi9PAlFA27dvp2PHjlFhYSHVq1ePunXrxstD8YLkSQgRd1Tn2jn2HU+IWux7JAQC4IZSXKI+eqZYDRG/Yek1tt7jx+vYVouBameaqE4tE6VYjbrf14K2JJBn8OKLL3LpDPgGsrKyuFTGxYsXWThat25N9957L02dOjXkUrSCkGh+iEgJhCDoTVAigbX9HTt2sC/g448/ph49erBvQOPIkSO0adMmzmhGGY1FixZx+QxBSGR8CQTsdghEmUMVB4UUsjuXFvi4MLSi1O4gu/NEm53IaIyc4Og5tt7jx+vYDrOhWkOngxIJZDC/8847lJSU5PN9WBHY4EtAAT50dhOERMafQEAZ2IIwqAJhsynqfjgs+f3wfBJlTkvEhg+1KXI3LD3H1nv8eB3b4VB4KwvniUIvkZg2bVrQA3bq1Ik3QUhUfAoEvipOCwKfbUO5QOClut/5IkRwrvZzyoy4gUTu36Ln2HqPH69jK0r5Vh0E7ZNAdNHixYvZWvB2esARgiUmX+8JQiLhzweBsEZ8+EFpmUKKXf2Q44ERT40QCH4dxs2m2FYuLqV2IksEn2r1HFvv8eN17GT4s9x8WlEjEi+99BLt3buXZsyYUeE9eMzhk4D3/NFHH430HAUh5p3Upc4PdJmihi6ajEYWhDLn0gEEAudgCSpUvAMUI1nYWc+x9R4/XsdW+IFC3aqDoOOn0MgHkUv+QOnu//73v5GalyDEVRSTajUoVMrrzeo+OChhRWhiEY5ACELUWBLff/89Zz37A+/hGEFINIIJc4W1AIGADMCqMCmOKlkPglBdBG1JmEymgFFLeM8981kQEoGgBIIUKrHh/+oxyIkQgRBihaDv6siq9u5B7Q66yuEYQUiUfhAHjxcHJRCIYsJrHMMO7DI1YkUsCCGulpumT59Ot956KzVt2pTuu+8+tixAWVkZvfLKK/T3v/+d/vOf/+g5V0GIbuvBLQ9CS5TTwlyxDwIBZyO+T7LIEpMQZyIxceJELrg3c+ZMjmBC8hzqz8APgWquv/3tb+lnP/uZvrMVhCheXgKaQKBPMfsgnGGtiGjiiBSOg+dEiRr6FwiCjlVg//znP3OJjiVLlnCJb/zRoyw3ynX06tUrxB8tCPEjELjvq0tI6nISiwLHuashrhANWWISYpGQS4VDDEQQhEQjkEBoWdRcqA91d5zWhCYQ4qQWYpmQwpFKS0u5PLggJBLB5EDAQV1aKgIhJLBIvP/++9zAp3PnzvTss8/qOytBiBKCDXGFg1qzKMSCEBJSJOCYXrhwITccguO6oKBA35kJQgwJhFZwzea2xCRZ1EJCiYTD4eBkOWz4HpsgxCvhCkR5oT4JcRUSzHH917/+lau8Ij8C/aOl+5yQqL0gtBBXFghxUAtxTtAiMW7cOBo5ciS3Lk1PT9d3VoIQbb0gSF+BwDjnLpWFdXZhsadVn1sQufI4eo6t9/jxOnZOnoMu5Tm4/3XUhcCiM52/7nSCEOvgZg+ByMkro8w0U9A5EFUViB/P2WnFZwWUWyBLuELw/OlXdahZg6ToEIkTJ05Q8+bNgx701KlT1KRJk6rMSxBqxIqAQAC1S5zBw/+AXtSuCKayyPSB2HWwhNbtKOIxkpMMZE0K/enQu68AKiFECj3H1nv8eB3baEDfazQjiiJLomfPnrzcdM899/hNpEN3umXLltHzzz/PvSV8NScShKheZsovX2bSkuT0clCjj/GHWwrom6Ol/PrKFhYafW0aC0Wo5BV6WiAZqZFb+tBzbL3Hj9exM9NMlJVhouxM1dqNCpE4cOAAPfnkk+yTsFgs1KNHD2rcuDElJydzW9NvvvmG9u/fz/vh4L7hhhv0n7kg6OSH0CKY8ACnh0BczC2j1VsL6cJlB/+M63ukUK+O1og/pQtCJAhK/rKzs+mZZ57hnhH/+Mc/qF27dnT+/Hn67rvv+P1JkybRrl276PPPP4+4QLRs2ZI/PN7btGnTfB7/2Wef+Tz+4MGDEZ2XEJ/lvt0tCHSSg0A4nOKg5UCo/arDE4hvT9ho2bp8Foj0FAP9YmQG9e6ULAIhRC0hOa5hOdx00028VRc7duzgcuQaX3/9NQ0bNox+/vOfBzzv0KFDlJmZ6Xpdr149XecpxEcUk0sg0EnOHrkaTHB4b/6qmPYetvHrFg3NdOPANEpPkUZdQpwV+KtuvG/uTz31FJcHGThwYMDz6tevT1lZWTrPToinRDmIAu9TFCouhZltiIiDGlFLy9cX0E8X1Z/Zvb2VhvVKISM8kIIQ5cTUYwxyNBYvXkx33XVXpeY5uuQ1atSIhgwZQuvXrw94bElJCeXm5npsQmJZEOyHsCu8xKT1f3AXiHBLbBw5VUqvr8plgbBaiMZcl0r9uiaLQAgxQ0yJBNqn5uTk0JQpU/weA2GYP38+LV++nFasWEFXXnklC8XGjRv9njNv3jyqVauWa2vWrJlO/wIhGgUC7gU4qu1OfwPEAg2D3AUiVODc3rSniJauzaeiEoXq1TbSLcMyqFVjS6T/WYKgKwbFOyg3ihkxYgQn87333nshnTd27Fi2PFatWuXXksCmAUsCQoGwXne/hhB/AuGeB5GTW0bFzkgmLQ49HOsBGbOrNhbQkR/V9aur2yVR3y7JZHaLa4+XcMxoHj8RQmDrZgXvMcB9DQ/Bod7Xot4noXH8+HFat24dWweh0qdPH16m8ofVauVNSFyBQBSTDY5qp+/BXqaQ2RRu9nQ+5Rao54/sm0pd21orfPAFIVYIS/7+/e9/07XXXsu5Erh5g+eee45WrlxJeoEy5XBGjx49OuRzd+/ezctQQuJRmUDAB1FSooa5wg+BKCSt1WgowCDfeaCYFn2UxwKRnWmkKaMzWCAEIaFEAnkSs2bNolGjRrF/QAtPRSQRhEIPUJYcIoEqtGazp/Eze/Zsmjx5sus15gDfBXI4kOCH9+GfmD59ui5zE2JbILREOa7HxM7p0KOYkD29cmMBrdmmltdA9vSUMZlUPztmDHVB8EvIf8UvvvgiLViwgMaPH8/hqBrItn7ooYdID7DMhPpRiGry5vTp0/yeewQU5oH6USkpKdSpUyf64IMPWNSExCEkgdCimJz5EaFwPqeMlq/Pr0L2dPlx4SZc889ym3YkM7f1HFvv8eN1bINbknBUisTRo0c5vNQbrOnr1a1u+PDhFYpeabz55pserx9++GHehMQlkEBACjiD2lsgwohi2n/ExvWXkF+B7OkJg9KpWYNQPlLqhxyfde0DH87n3mQy8L+hvPhb5G4eeo6t9/jxOraB/17Cf6jQXSRatWpFe/bsoRYtWnjs/+ijj6hjx46RnJsgVF0gnDWXIBBwUGMfwlvhmFYruqphr6EsL8Fv8cnOItp5oCSs7GmnHLjEwej2gQ/nXmMxEWk1CUxmInMEA5D0HFvv8eN1bKOBCIFyJmOUigR6XaNuUnFxMX/A0PP6rbfe4lyD1157TZ9ZCkKYhfq00t5aqW879jkFAk9rofofkD29Yn0+/Xhe/Tn9uiTTgG7BJ8dpR0EUjLxkQAQ3G/YbnE+RoeJQDFSiFpMlq8VASZbIPWLqObbe48fr2CaTgTf3kOqoEok777yT7HY7L+kUFhbS7bffzr0jUCL81ltv1WeWghBmJVd+TWoEk7akBCsinDIbyJ6GgxrJcSjpPbZ/Kl3RLLimL66PM6wGfNCNBqc44MNOZDQZOIokHEuCz3S2z0tJMpLFHMH1cR3H1nv8eB3bbCT+mwknRDscwgq/QF8JbKgEi8gjhKYKQjQKhE8HdYj+B1jMKM63aU8xv25Yx0QTBqVR7YxQPqVQB4UFAhYEBAJCwWvOVRAIQdCbKsXo1a1bN3IzEYQwSn2fv+ysyqeTQPjKnh7eK5XMQT85Oj0QRoULBkIYTCY1m1sTB/e15XAiVnC+wakwPGYknag6jq33+PE6tiHao5sQ2eRrctiHUuJt27bl2kqDBw+O1BwFoUYimPxlTwdPecSS2YiFJXz4VWdkpARCEKJOJNCdDgl1Xbp04VamnGm6cyft3buXxQFd6oYOHcrlM2688UZ9Zi0kNJUJRFkVI5hw/L7DNtr0VTGPj+zpmwalhZAc5xnaigfB5GR1Tvie15md4qBGNmnWBoWN0VYeGQNLJZLoObbe48fj2AajWlsMW1SKBPwQv/nNb+gPf/iDx/4nnniCS3SsWbOG5syZQ3PnzhWRECIOLAIIRE5eGRc68yUQyIAON4IJ567fVUTfngi393S5QLDPweh0UiMU16E6G5NgRWjiwD4J77NDhGPp1W/ZGR5Jg0TPsfUeP07HNji36vJhhaxFy5Yto9tuu63CfkQ24T2A99EZThD0sCIgEIDLaDgFgkNcy9QlJggErIhQBQK9p5d9ks8CgQ/mkJ4pbEEEJxDqR5dDW+F3cEYv8Vdn9BJ/1RKhsPxkLBcI7YMvCNFGyJYE/A5btmxh34M72If3ACKepKqqoMsyU75zmanM00GNHAheZgrT/+CePZ2WbGD/w5UtksKyHrA8wBFMbqGtJvSqcKqAwQyntffZghAnIjFjxgyaOnUq7dq1i3r27MkmMxLqkEj3yCOP8DEff/yxz9IdghAJP0SZ3TNJrioOalgd63YU0a6DavZ00/omGtEnlVKTjWH5HmA9QCDgkGbrwWk58GxgPRjUdebyswUhDpsOLVmyhF566SXXkhK6v0E8kFgHioqKXNFOsUi4zTkE/QUCzt9Luaq3Oi3FwE/+vLTEyXGhOai9s6fRe7pPZ6sr5DBwIxn1hk8BEuO0yCUMpxjUZjI4JyvN6ApxVN93+z5MZ6S2BIdmNJFGz7H1Hj9ex84OseFQVe5rMdWZrroQkYjeSKbiEoVvtviztVgMrvajoQqEr+zphnU8P3S+RcLgs6wGSiSwM9FkqBC5BN+DYlAoN9/Bx9RKN7oimjSBwIdeSosLelLtnelQkvvs2bPsf3CnefPm4Q4pCJVnUtvVJSWOXFIUfooPxUHtK3sazmk8sQXuHud7acmX9aAZCbAWNN8DHsXcl5fcrYe6WaE/FQpCdRHyXyaa+aCvAxzV3h8+fHi0JkSCoEepDfgatDBXi7MwXrD+h/Czp4NwTLtbD07fg5b34B65hPM0gRBxEOJSJJAwh+5w77//PrcElSxRoboEApZCiV0VCGCz48EkuDFPnbPTuyFnTwfvmHZfWnIXB/IhDnVqmWVpSYhfkUAvCUQ2tW/fXp8ZCQmNh0CUlfeCQJgrBKOkxME+CK3VKFkCjwcLF5FLiGAKJXvavaR3IOvB19KS+/kaOBdLWuJ7EOJeJNBYCFnXgqC3QLD1oJTnQcByKMU+55KTKhP+TQlYHMh9+OZo8NnT7iW9EeWkWQ/uFVstlSwtuSO+ByHhROLpp5/mXhJPPvkk12+yWDwf5SRkVIhEHoR3oT4tzNVTICLbe9pdINDPBT4Pd+sh2KUlIOIgJKxIoHgfGDJkiMd+cVwLkRYIdlCjzIYWzRSkQITXe1rrIaqQyWDgRi8Ia2VLAp3AglxaEnEQKNFFYv369frMRKBE7wdRuUAoAQXCO3s6+N7T6m0eQmAyIIfB2VvY7GU9BFhaEnEQ4pWQRWLgwIH6zERIKAJFMOErMqvdW41WFuYaXu9pz+illGQj52GAVKtazls7HctNlS0tiVNaiEfCzuBBf+sTJ05wUp07Xbt2jcS8hDgm0gIRXu/pimU1HA4DOcoUtiCsKMLntB5kaUlIZEKuFHPu3DkaM2YMZWRkUKdOnbiQn/sWSR577DGPVn3YGjZsGPCcDRs2UPfu3bluVOvWrenVV1+N6JwEfQQC2dP4Wlzi4EJ93DDImTznD/jBNu0poqVr81kgkD1919iMSgSivKQ3xCHJovof0E4UZT74Nce3OkNdTZUnxLVvaZWMaSFuCVkkHnjgAbp06RJ98cUXlJKSQqtXr6Z//etfdMUVV9CqVasiPkEI0enTp13bvn37/B579OhRGjVqFPXv3592797NVWlnzpxJy5cvj/i8hAjkQGhZ1Ogcx3WZHEHXYSoqcdCqTYWu8hrInp58Q0YlxdacPgWnQGg1l7TopST2Sag+B67FhOglEQchwQl5uenTTz+llStXcplwo9FILVq0oGHDhnHo67x582j06NGRnaDZXKn1oAGrAbWjnnvuOX7doUMHbq36zDPP0MSJEyM6L6HqORCwHrxDXNWOcoGXl85csNNHWwspvzDY7OnAZTVc5bxRXwkqAUNCbRrnQvwOQqISsiVRUFBA9evX5++zs7N5+QkgZ+LLL7+M+ARRK6px48bUqlUr7n535MgRv8du3bqVhg8f7rFvxIgRLBSlpWpClS9KSkq4QqL7JugnEAhNdTith5KS8l7UlQkElpe++q6Elq8vYIHISjfSlNEZlQqEZj3AOtBCW5EDkYRlJm4v6hQPc/mxWslu76UlqdQqJBohiwR6R2h9JK6++mr65z//SadOneKneNRyiiS9e/emRYsWcROjBQsW0JkzZ6hfv3504cIFn8fj/QYNGnjsw2u73R4wSxwWEEroaluzZs0i+u9IZLxzICAQmoNaC3HVWo0GEghkT8M5vXF3Mfsy2jQx081D0wPctCv6HrSlJbPTB6HlPrDvAcUC3ZaXxO8gCGEuN8EnAd8AmDNnDj+powlRUlISvfnmmxRJbrjhBtf3sFT69u1Lbdq0YR/IrFmzfJ7jnVGrtcsIlGk7e/Zsj/FgSYhQ6Jckh68l3E2u3IoI1NbEO3v62q7J7IPw/zstj1yCOPAN3xhc3oPkOwhCFUVi0qRJru8RzXTs2DE6ePAg+wLq1q1LepKWlsZigSUoX8B3AWvCHfS8gF+jTp06fsdFP27pyV19AgEHNdAimAI5qL2zp9FatHE9c9C+B19Z0/5Kaqiio7Dzu13zpEpyLAQhMahyp5PU1FS65pprqDqA7+DAgQMcveQLWBrvvfeex741a9ZQjx49KtSYEvQDS0cQCLRgzEgxuZzUoUQw+cue9m9wqFaBWpDPUMEx7er34CcpTrMgIBDwPYhACEKYIoGmQlhW+uSTT3x2pkP0U6R46KGHaOzYsWyl4Gc98cQTvBR0xx13uJaJ4A+B3wJMnTqVe29j6eiee+5hR/brr79Ob731VsTmJARnRUAgEOJqs+Pvw8DVW202VSAqc1Bfzndw7wdf2dMVu8d5Wg9sMXD0kuqYrqykhiwvCUKEReL+++9nkUCoa+fOnXVtOvTDDz/Qbbfdxk7nevXqUZ8+fTg/A2G3AL4RZH1rIALqww8/pAcffJBefvlljop64YUXJPy1BpaZVIGATwB5EKpAuNdgikT2NDuZw7QeRBwEITgMSiCPoQ/gd8CTO5LW4pVwG4YnOppAIGrpUp6DBcGahB7UVKlABOo97Q4sCfeS3tmZZrVTnInI5Oz1IOIgCJG7r4VsSSCKqW3btqGeJiSKo1qB3wH+BrVBkENx8DJQIP9DKL2neY9mPRjVvAdvxzQsBz7MGdKqIQlxglANeRK/+c1v6Pnnnw8YsigkVrnvg8eLXQJRYiMqdagCwb4H/upfINB7+o33clkgsFQ05rpUGtUvzadAsHMawgCLwWRUQ1rNBjUpzumH0JaW3AVCEuIEIXyCsiRuuummCs7pjz76iOsqeUcNrVixogrTEWI9zBVRTJolgZUl5EOkJEei97TqeE62GNQIJwNRSrKBrM7VKPE7CEINigTWsdyZMGGCTtMRYl0gYEFAGFxhruyDMFSh97RnYhyWrkpKHbzMlGw2+lxaEqe0IFSzSCxcuDCCP1KIu1Ibzr4P7mGu/pzUwfee9hfaqvCxWCc1WvyHtEoDIEGIDCE7rlGOG7WQUBrcHWRBY+mpZcuWEZqaECsCgeUipENUFuYafO9pz8Q4zp7WynfD0V2Ccq2qb0KsB0GIMsf1lClTaMuWLRX2b9u2jd8TEiOCSRMILtRnV1ggUGYD7T+9BQL+iY+/KOT8BwgEsqfvGpfpQyDKi/JxIyCzUS3K5+acRu9pHKBaGFKITxCizpJAM59rr722wn4kuk2fPj1S8xKiWCBspeWd5LTlJX91mPIKHLT808Igek97Wg/uNZcgDgD+B4VrLylk8irhXTeryhVmBEHwQcifLKwH5+XlVdiPBA2U7BDiWyAQ4uqrUB82b46fKaU1XxRRsS1Q9nS57wGOafwH3wMsB/e8B833gDlAIDS/Q4M6UpNLEPQk5OUmFNdD/wV3QcD32HfddddFen5ClAkELAiU29AK9ZXaKwoEwlu37S+mVRsLWSD8954uXzayWoyuZkAIc+WlJWeJDfc+D2ZnNyBkYotACEIUWhJ/+ctfaMCAAdx8SKvGumnTJk75jmRxPyHKsqhtFSOYfBXq4+zpTYV04oyaPd2pdRKN7uedPe1pPWj9HmA9IKHO4KPekntYa9tmSbK8JAjRakl07NiR9u7dSzfffDNXZsXS0+TJk7mnBAr+CfFXqA8CgQgmzUGNCCZYEN4CoWVPQyBwsx/aK4VDXL0FooL1YCJKSdKsBwOZLOXWg7fvoX1LqwiEIERzgb9EINEL/LkLREmp6pDmSCbn0pK3g9o7exq9p2/ol8o3dZCRagwY1sqVW41uzYCc40rOgyDEYIE/IfEEAhYElpZ8Oah9ZU8P7JbC1V/L0cJaceM3BnRMSwtRQYguRCQET4HILWPfAwQCEUylAQTCX/Z0fpGzr7jzfxAD5DxwOW9nDgQLhFdJDSmnIQjRh4iE4CEQEAaIASKY7KXlEUze/geP7OlUA00Y6J49rbgEAlnRSRajh2Pa1SnOaT1oIa1Ach4EIUZFAh3gmjVrpmsnOqFmBUJdVgocweSv93R6SnkMBNdWMijOHg+qQ1qNYvJcWhJxEIQ4Egm0BkW70Pr16+s7I6Fae0Gcv2xXs6jhlLZrXyEMxF/dHdSBek+XgxIaangrEuPgm8Dm7ZgWp7QgxJlISBBUfIEifBCInNwySrYaXSGu/gSi8t7T5dFLKVYj2ezoSEeUngxfhFgPghCriE8igZeYLl1WfRBmU3knOW8HNcQEvaex+e89rQoEfA6qVaHwV1gP5iRVIGRpSRASQCRee+01Sk9PD3jMzJkzqzonoRoE4oLTSQ2nNKKOuNy3l0AgexrWw1G/vafLQ1tNxnLnNL4WIMLJYKAkLtQnhfgEIe6T6YxGIzVt2pRMuKP4G8xgoCNHjlCsE8/JdFoeRFGxg0NXYSkkJRk4L8LdQY3safgfcgvggCa6oW8qdWlrdRupPBkOy0lmo9Y5jrhSa16hg8UBiXV1apkDtCUVBCFukul27twpjus4yYPILyyPWiouVshsDrb3tKf1oCXGqVnTFX0PIhCCkCC1m2oi9BWVZXv27EkZGRksTuPHj6dDhw4FPOezzz5To2i8NtSWSmTcw1xRi6kINZh4uUmtzaRlT2N5ac02VSCQPT1lTKZLIHAdYSlwAyCzkUNbuWKruVwgIA7Y4LNAIT6xIAQhtonq6KYNGzbQtGnTWCjQMvXRRx+l4cOH0zfffENpaWkBz4WYuJtU9erVo0TFXSAQ5srZ1HbVD1GKUCYOh3UE6D1dufWApSXtMUIS4gQhAUVizpw5lTqtI83q1as9Xi9cuJAtil27dnG58kDguKysLEp03AWitLS8DhNHMjmXm749YaNPdxb5zJ5WE+OcpTScvgfkPXC4q7OshnsTILEcBCFBl5tmzJhBFy9e9Ni3f/9+uvPOO7ls+H/+8x/SGzhcQHZ2dqXHduvWjRo1akRDhgyh9evXU6ILhM1NILhRkEP9fsOXRfTxF6pAtGhkprvGar2n1WU6WAwoq4GIJpTztlrVMFcIhLa0pFkPIhCCkMCWBJZ9cNP929/+xq/RSwJNhxo3bkxt2rShKVOmcIe6X/7yl7pMFMtds2bN4u53gfpWYI7z58+n7t27U0lJCf373/9moYCvwp/1geOwuUcBxJNAoKOcWuJbFQg4rNF7+qOthfTTRWf2dNdkGnC1mj0NcVArthpYHGRpSRASl6BF4osvvuDlHo1FixbxE/2ePXvIbDbTM888Qy+//LJuIjF9+nRudrR58+aAx6FjHjaNvn370smTJ3l+/kQCDvLHH3+c4ioP4nIZ94DwJRDInv7fhgJ2WFstRMN6p1JXDm9VcxxQlA9LS66aS86KreJ3EITEI+jlpjNnznD9Jg20Kp0wYQILBBg3bhx99913ukwSS12rVq3iZSPkaoRKnz59As5t9uzZvJSlbRCVWBcIWBBFxWpyHMShxIavDtq4u4iWrs1ngahX20i3DMug1o0t6tKSRe0WZ7EYuPQGlpbYD4Ey307HtHSHE4TEImhLApFCOTk51KJFC369fft2uvvuu13v4ybjvmQTqSUmCMS7777Ly0XuIhUKu3fv5mUof1itVt5iHU0gEL2EGktaWCvEwTt7Gr2nUZwPVgMsBavF0zHtHbUkTmlBSEyCFolevXrRCy+8QAsWLKAVK1Zwb+vrr7/e9f63337LpcQjCfwgcIivXLmScyVgzQBkDaakpLisgFOnTvHyF3juueeoZcuW1KlTJ7LZbLR48WJavnw5b/EMBOLMRTuX2oDV4L685Ct7ulVjC5dj5TIazuglWA+c6yBLS4IghCoSc+fOpaFDh/JNFzkLjzzyCNWuXdv1/tKlS2ngwIEUSf7xj3/w10GDBnnsh28EjnKA8uXodaEBYXjooYdYOCAkEIsPPviARo0aRXFb7jvHruZAeAmEvcxRMXt6cDrVr21iSwOCwfkPFiOlJKs+B67Y6hxb8h0EQQi6dhM4d+4cbdmyhRo2bEi9e/f2eA834o4dO4a9JBRNxErtpkARTEUljgq9p0dfm0bJSUZnrwf1eHxfO8NEaalqTS7N71A3S2otCUI8Ee59LSSRSBRiQSTcBQIOaq2LHCyJczl2n9nTKNKohbaivEYJKsAaYWGYuYOc+B0EIX7J1bvAn7bmXxmTJ08O+ocLVReIwiK1NRAEotjm8NN72uIq4605py1JREpReaMgsR4EQahyqXCU5UDIq79TEOHknZUdi0SzJaEJBEcwFasCwctLxWWevacbmenGAWmUkWpyWQ9wTnM5bxMc1UbKySvjpab2LWM/sksQhBq2JDp06EA//fQT/eIXv6C77rqLunbtGvQPESIb4lpkU8hmUwWiuMTBouHRe5qzp1Nc2dKuxLgkWBFGV59ptB8V34MgCBFJpkOdJjini4qKOHO5R48eHH0UDyUsYk0gStwE4ruTNnrjvVwWCISw/nxIGg3unsp9pvEa1kOq1cB9rGE9QCDqZKrWgwiEIAi6OK4hFO+88w6HoiKpDn0e3njjjbhISIvG5SYtBwJ+BiTHqQ5qB23YXeTVezqd6tU2sRMaHeOwtISEeIgDkMglQUhccmsiumnjxo1cQhxfz58/75E3EctEk0hoAgHrAXkQEIhLuWX0v435Hr2nR/ROo7QUtc80BCLZ6rm0JI5pQUhscqujfSlAktq//vUvtiIKCgrYR4Flp3gRiGhCE4jiEmROqwLx/SmbV/Z0GnXvYGXrARtKaqD+ktZGtJ7kPAiCUAWCFolly5axMKBb3IgRI+jZZ5+l0aNHkwmhMkLEQTkNCMTZi3ZKSzGR3e6gz/cWeWRPTxyshrci1wHO6dRkdI1T86VFHARBqPYQ2ObNm9OkSZOoQYMGfo+bOXMmxTrRsNx05FQJHTll4+/Ri/q9zfke2dPwP0AUtKUli6ncKS3NfwRBqHafBIrmqf2OAwxmMNCRI0co1qlpkcAy0/6jJby8hKZAWF66mKtmTw/pmUL9r0rhqCVVIMR6EAQhCnwSx44dC3pQITyg1z9dsNMP5+wsEPu+t9HqreXZ07cMVXs/oIw3nNSytCQIgt5IFbcoweFw0E8Xyujk2VLOoF69tZC+PKRmT6Pn9O3DM6hWupGSk02UhEqtRvE7CIIQRSKB3IhPPvmExowZ4+rj4N5kCA5slBNPTk7WZ6ZxLhCnztrp9AU7h7cuXp3HPSBAn85W7j2N8hmpTusBhfgaZIu+C4KgPyEV+Hv//fddIvHSSy9xrwat+c/BgwepcePG9OCDD+o32zjk3KVSOn3ezu1EDx230dtr86iwROFQ1jHXpdIVTS2UkmykjHQ1ikyilgRBiEqRWLJkSQUBQNe41q1b8/doRvTyyy+LSITooIZAoPfDuh2F9OnOIlf29K1D09lyQAugrHSziIMgCNEtEmhP2q5dO9drLCshLNa9vSnajQrB8dOFUl5SupxfRm+tzaPDJ9Xw1h7trTS2fxpZk4xcgqNOlpk6SJVWQRCiXSQQNoUy4e5d6rzX1d19FEJggTh+ppROnCmlJR/n0eV8B2dP3zggnbq3t3I70VSrkbJrie9BEIQYEYmmTZvS119/TVdeeaXP9/fu3cvHCJULxLHTNtqyr5g+/LyAy23UqWWkSSMyqHE9C1dvRYlv8T0IghBTIjFq1Cj64x//yKU4vCOYEPn0+OOP83tC4BwIlPZe8VkeffWdmk3dqVUS/WxIGtXKsJDVYqC6tUxUr7ZELgmCEB0EnXGNhkNXX301JSUl0fTp09k/gQxrRDUh0slut9Pu3bsDluxI1IxrLMWdPl9Guw8VcXjr2Utl3DJ0ZN9U6n91CmWmmdmagDhUltUuCIIQlRnXuPlv2bKF7rvvPvr973/vamGKm9qwYcPolVdeiQuBiDRlZQ46+ZOd1mwroBXr88hmJ8pMM9DtwzOpQysrNchWLQcRB0EQopGQ1jVatWpFq1ev5j7Whw8f5n1t27al7OxsveYX06ByK5aXFq/Opa371OZAbZqY6dbhmdSptZXq1zZ5RIgJgiBEG2EtfkMUEPIq+Ke01EHbDxTT/BU5dPKsmj09uHsKTRiYQfXFehAEIUaIicdYLGXBioHDvHv37rRp06aAx6PnBY7D8Uj2e/XVV6k6sZU6aNXmfJr35gUWCJTynjI6k+4dn0Wd2lipfrZFlpcEQYgJol4k3n77bXrggQfo0UcfZcd4//796YYbbqATJ074PP7o0aMciYXjcPwjjzzCPS6WL19eLfO12R304rJL9Mo7OVRYrFCTeiaa9+u6NGlkJjWoI+IgCEJsUaUe19VB79696ZprruEWqRodOnSg8ePH07x58yoc/7vf/Y5WrVpFBw4ccO2bOnUqffXVV7R161ZdowCQPT33jQv05UHV/zD4mlR64NYsSk8zieUgCEKNEu59LaotCZvNRrt27aLhw4d77MdrRFr5AkLgfTzare7cuZNKS9XSF94gUxwX0H0Lp93og38/ywKRZDHQfTfVokfurEMZ6RK5JAhC7BLVInH+/HkqKyurEFqL12fOnPF5Dvb7Oh55HBjPF7BIoLDa1qxZs5DnajIa6I7RtahpfTO98nAD+vnQWq6mQIIgCLFKVIuEhncOAVbIAuUV+Dre134N9MaACaZtJ0+eDGueA69JpdcebUStmySFdb4gCEK0EdX1H+rWrcvNjLythrNnz/pN3GvYsKHP41GcsE6dOj7PsVqtvEUCLDUJgiDEC1FtSaAECEJZ165d67Efr/v16+fznL59+1Y4fs2aNdSjRw+yWCy6zlcQBCHeiGqRALNmzaLXXnuN3njjDY5YQuMjhL8iYklbKpo8ebLreOw/fvw4n4fjcd7rr79ODz30UA3+KwRBEGKTqF5uArfccgtduHCB/vSnP9Hp06epc+fO9OGHH1KLFi34fexzz5lA0h3eh5igUx5aqr7wwgs0ceLEGvxXCIIgxCZRnydRE8B5nZWVxQ7sSFSBFQRBqGkQ2o/IzZycHI7ijBtLoibIy8vjr+GEwgqCIET7/S0UkRBLwk//hx9//JEyMjL8hs1qqhyr1kYszz+W5w5k/nL9a+LvB4tGEAgswYdSfVosCR/gAgbbihW/pFi8UcXD/GN57kDmL9e/uv9+QrEgYia6SRAEQag5RCQEQRAEv4hIhAkytOfMmROxTO3qJpbnH8tzBzJ/uf6x9PcjjmtBEATBL2JJCIIgCH4RkRAEQRD8IiIhCIIg+EVEQhAEQfBLwooEutH17NmTs6rr16/PPbMPHTpUIUPxscce4wzFlJQUGjRoEO3fv79C69MZM2Zw74u0tDQaN24c/fDDDx7HXLp0iX75y1+6Ot/he9RPiYb5Yx+yyt23W2+9NSrmv2LFCm49i2uLee3Zs6fCONF8/YOZf7Ref7T6Rb/4Ll268HXF3xCqLaMSQSxc/2DnXxPXf14Qfzv43LZv357nXrt2bRo6dCht27atZq69kqCMGDFCWbhwofL1118re/bsUUaPHq00b95cyc/Pdx3z1FNPKRkZGcry5cuVffv2KbfccovSqFEjJTc313XM1KlTlSZNmihr165VvvzyS2Xw4MHKVVddpdjtdtcxI0eOVDp37qxs2bKFN3w/ZsyYqJj/wIEDlXvuuUc5ffq0a8vJyfH4WTU1/0WLFimPP/64smDBAhShVHbv3l1hnGi+/sHMP1qvP+YwdOhQ5e2331YOHjyobN26Vendu7fSvXv3mLj+wc6/Jq7/iCD+dpYsWcLX9Pvvv+fj7r77biUzM1M5e/ZstV/7hBUJb3Dx8UHesGEDv3Y4HErDhg35RqtRXFys1KpVS3n11Vf5Nf6YLBaLsnTpUtcxp06dUoxGo7J69Wp+/c033/C4X3zxhesY/MFiH/54a3L+2ofk/vvv9ztuTc3fnaNHj/q8yUbz9Q9m/rFy/TW2b9/Oxxw/fjymrr+/+UfL9T8bxNwvX77Mx6xbt67ar33CLjf5Kg8OsrOz+evRo0e5Derw4cNdxyB5ZeDAgbRlyxZ+vWvXLjZr3Y+BWYueF9oxW7duZTOvd+/ermP69OnD+7Rjamr+GkuWLGGTtVOnTtycSauCW5PzD4Zovv6hECvXH8dgOQZl9GPx+nvPP1qu/+VK5m6z2Wj+/Pn8M6+66qpqv/ZS4M+5do9Odtdddx1fZKD1yfbupY3X6HynHYMWq1gz9D5GOx9fse7oDfZ59+Ku7vmDSZMmcaMm9Ab/+uuvudPfV1995WoBW1PzD4Zovv7BEivXv7i4mH7/+9/T7bff7ioqF0vX39f8o+H6KwHm/v7777N/pLCwkBo1asRzgphV97UXkSCi6dOn0969e2nz5s0VLpB3qXD8Uv2VD/d3jK/jgxmnOuZ/zz33uL7HH+kVV1zB/cC//PJLuuaaa2p8/uEQTde/MmLh+uOJFTcrlNB/5ZVXYu76B5p/TV//6QHmPnjwYA52OH/+PC1YsIBuvvlmdl77uvH7m1ck5p7wy02IDli1ahWtX7/eozw4niyAt+KePXvW9XSOY2AKIoIg0DE//fRThQt/7ty5Ck/51T1/X+CDYbFY6LvvvqvR+QdDNF//cIm2648bLG5OWL7Ek6z7U3gsXP9A86/p6z+jkrkjYqlt27a8RPT666+T2Wzmr9V+7ZUEBY7dadOmKY0bN1a+/fZbn+/D8fv000+79pWUlPh0XCOCQuPHH3/06Tzatm2b6xg4kqrq+IrE/H2BKCh3J1pNzT8Ux3U0Xv9g5h/t199msynjx49XOnXq5BFVEyvXv7L519T1d4Twt+NOmzZtlDlz5lT7tU9Ykbjvvvv4hvnZZ595hL8VFha6jkFkEI5ZsWIF//HcdtttPkNgmzZtylEHCEO7/vrrfYahde3alSMLsHXp0qXKIYCRmP/hw4c5RHPHjh18I/vggw+U9u3bK926dYuK+V+4cIFvrJgX/rARyYHXOC4Wrn9l84/m619aWqqMGzeOry3CNN2PwcNGtF//YOZfU9f/vkrmjlDY2bNn8886duyYsmvXLg6BtVqtHA5b3dc+YUUCH1pfG+KX3RUfyo0ncvyCBgwYwDdbd4qKipTp06cr2dnZSkpKCv8CTpw44XEMbhaTJk3inAVs+P7SpUs1Pn/ME/sw96SkJH5SmTlzJs83GuaP730doz1NRfv1r2z+0Xz9NevH17Z+/fqov/7BzL+mrj9VMndc0wkTJrClgXnhwQ6ChxBed6rr2kupcEEQBMEvCe+4FgRBEPwjIiEIgiD4RURCEARB8IuIhCAIguAXEQlBEATBLyISgiAIgl9EJARBEAS/iEgIgiAIfhGREIQIgWRatJlEy1JvUH0UdfxPnDgh11uIKUQkBCFCoPzywoULuZzzP//5T9d+VCBFv+Xnn3+emjdvHtHrjSqngqAnIhKCEEGaNWvGYoAOZxAHWBd33303DRkyhHr16kWjRo2i9PR0LtWMpvToFaCxevVqbj6Dzml16tShMWPG0Pfff+96/9ixYyxEy5Yto0GDBlFycjItXrxYfn+CrkjtJkHQgfHjx1NOTg5NnDiR5s6dSzt27OBmNmhyM3nyZCoqKmLrwm6306effsrnLF++nEWgS5cuVFBQQH/84x9ZGNB4xmg08vfootayZUt69tlnqVu3btySFm0rBUEvRCQEQQfQ/AWdzi5cuED//e9/affu3bwM9fHHH7uO+eGHH9jyOHToELVr185ncxh0Idu3bx+PpYnEc889R/fff7/83oRqQZabBEEHcHO/9957qUOHDjRhwgRuXI8OZFhq0rb27dvzsdqSEr6iB3Pr1q25gxoEAXg7u2GRCEJ1IT2uBUGvD5fZzBtAf+WxY8fS008/XeE4NLkHeB+WBfoZYwkJ58CCQJtK77aWglBdiEgIQjWA3snwOcCfoAmHO1iWOnDgAEdF9e/fn/dt3rxZfjdCjSPLTYJQDUybNo0uXrxIt912G23fvp2OHDlCa9asobvuuovKysqodu3aHNE0f/58Onz4MDuzZ82aJb8bocYRkRCEagDLR59//jkLApLtsIwE5zMS7BC5hG3p0qXsu8B7Dz74IP31r3+V341Q40h0kyAIguAXsSQEQRAEv4hICIIgCH4RkRAEQRD8IiIhCIIg+EVEQhAEQfCLiIQgCILgFxEJQRAEwS8iEoIgCIJfRCQEQRAEv4hICIIgCH4RkRAEQRD8IiIhCIIgkD/+H4XC2g5z0l8BAAAAAElFTkSuQmCC", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAADcCAYAAACF6V1NAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAAOWNJREFUeJztnQeYU2XWx0/aZCoMQx16k16kSVHK0AVBEFcBV2yPuyrYWLegu4t87qfiqh+6FtYCroprA1YWpAhSBaSrICB96H0YpiaZ3O/5n5ubSTLJTDIkmZTz271OcnPzzsud5P7veU/TKYqikCAIgiB4Qe9tpyAIgiCISAiCIAjlIpaEIAiC4BMRCUEQBMEnIhKCIAiCT0QkBEEQBJ+ISAiCIAg+EZEQBEEQfCIiIQiCIPhEREIQBEGITJFYt24djRo1iurXr086nY7+85//uL2OiiHPPvssv56UlEQDBgygPXv2VDju/PnzqV27dmQ2m/nnwoULQ/ivEARBiF2MVfnL8/PzqXPnznTffffRuHHjyrz+0ksv0auvvkoffPABtWrViv72t7/RkCFDaP/+/ZSWluZ1zE2bNtGdd95Jzz33HI0dO5YF4o477qANGzZQz549/ZqX3W6nU6dO8e+AeAmCIEQ7uOm+evUq33Tr9QHYB0qEgKksXLjQ+dxutyv16tVTXnzxRee+oqIipXr16srs2bN9jnPHHXcow4cPd9s3bNgwZfz48X7P5fjx4zwf2eQcyGdAPgMUY+cA17dAqFJLojyOHDlCZ86coaFDhzr3Yfmof//+tHHjRvrtb3/r05J48skn3fYNGzaMZs2a5fN3FRcX86ahFcY9fvw4VatWLQj/GkEQhKolNzeXGjVq5HMVxhcRKxIQCFC3bl23/Xh+7Nixct/n7T3aeN544YUXaMaMGWX2QyBEJARBiCUCXULXR9s/CHf5Ff0jA33PtGnT6MqVK84NFoQgCIIQwZZEvXr1+CcsgMzMTOf+c+fOlbEUPN/naTVU9B4sY2ETBEEQosSSaNasGV/wv/nmG+c+i8VCa9eupT59+vh8X+/evd3eA1asWFHuewRBEIQItCTy8vLo4MGDbs7qXbt2UUZGBjVu3JieeOIJev755+m6667jDY+Tk5Np4sSJzvdMmjSJGjRowH4F8Pjjj1O/fv1o5syZdOutt9JXX31FK1eu5BBYQYhULuTY6FJuSaXfn3NVfW96miGIswr92KEeP1bHzqhmoFrpxtgXiW3btlFWVpbz+dSpU/nnPffcw7kRf/jDH6iwsJAeeeQRunz5Muc5wCpw9c5nZ2e7xfzCYvj000/pz3/+M/3lL3+hFi1a0GeffeZ3joQghJtzl6x0yXFRsNsDf79dUchWokbkWW32oOb2hHLsUI8fq2PrdTr+nNjtCun1oc/j0jlyFASPULHq1auzE1uim4RIFgh+n6JQbp765uqp+qBfsEI1dqjHj9Wx9Tqdw5IwBCQSlb2uRazjWhBifXnpwhWb+kQpFYjK3LEpdkTwOR4H+ZYvlGOHevxYHVvRlWbGhQMRCUGoSuuhhKjERSBwgQgUm610jCILkSmIS+ShHDvU48fq2CajelNRUkJk0MdxdJMgxCLBFghBCDViSQhCFQoElhPsimOJgRR+HCgWm0I2m8PRabUH9d4vlGOHevxYHVuvx2dH4Y0o9I5rsSQEIQwgEgUCgdDGYAqEIIQasSQEIUxWBMe+K+QIb9Q5BaLEbnc6KCsjFAiTtDneaLHhTjN4ahPKsUM9fqyObTdqIbAUmSKBaqlbtmyho0ePUkFBAdWuXZu6dOnCGdKCIPhYZsrTBEK1IDQfBAQC+wC+80olVMJiVajEEVdvwZfaELwLVijHDvX4sTq23a7wpi43RZBIoDz3P/7xD+4eh/IY6enp3C3u0qVLLBzNmzen3/zmN/TQQw8FXIpWEGLeD+EqEI5NEwh813mzK5UKp7S6RMMgxj6h8onbYR071OPH6th2fHb480NhwS+fBMpb3H777Vz+Yvny5dzd6OLFi3TixAm2Jg4cOMAZzqtWreIOcp61kwQh3h3VVqu6PAARsNkhGHZeStAEAneOktYqRCJ+WRJo/PPFF19QQkKC19dhRWBDOQ30oEbrT0GIZzwjmayOu0EYClaLnRRHFi0EQivBwK9XwpLA250rD1irDuIqRCjHDvX4sTq24mKNRoxITJ482e8B27dvz5sgxCveQl3Z0ehYT7YoRCaTzikQWmWcyn7xef3a8b4Sx+8IFqEcO9Tjx+rYSphFwu8QWBTYg08C9T88QS0QX68JQrz6IDSBUMNdVWcknJL4bmNfMARCECJGJN544w1at26d18JQKBq1fv16FgpBiFdcBQLORzilWSDsRMUWO1m1BCl2TKoCEe67QkEImUjMnz+fI5d88dvf/pa+/PLLgCcgCLEoEK6JcohisljV47DsZCmxuwsEx76qm7q/Epvr/yo7RlWMHc1zp6ofO6JCYA8dOsSNf3yB13CMIMQbFQkEh7k6BIJj23Uu1oPL91xdiBKEKLUkDAZDuVFLeM21+Y8gxAMVWhBamKsdEU5q6Q0RCCGa8PuqjqxqJNL5YuHChXyMIMRLP4h9x4oqFAiAfRAIDefykvZcLAghgvF7uWnKlCk0fvx4atiwIT388MNsWYCSkhJ666236P/+7//ok08+CeVcBSEKqrmq5RJcBcIzikkEQohJkRg3bhz3nH7sscfomWee4eQ5tNWDHyIvL49+//vfc1a2IMRjDoRaiwniAJFQj0VTmBKnk9FdIMR6EKKFgAr8/e///i+X6Jg3bx4dPHiQP/z9+vWjiRMn0g033BC6WQpCBFCRQBRbS8t9ox8ABAJ480EIQrQQcBVYiIEIghBvlCcQSJJDDoSzUF9J+QKhWhGhbxYjCMEgoHAkq9XKBf0EIZ4IvkAIQgyKxOLFi6lFixbUoUMHeuWVV0I7K0GIwBDXEtcsal5SsjuXmDQHtQiEEGvoFD/T9tq2bculOTp37swRTigVnpKSQrEIalCh1AhqUnkrQyLEB/7kQACvhfrs7tbD1QL34v9pye73Z/uPWejgCUdadoBo5T40TMbgLWWFcuxQjx+rY5uMOkow6ejOwdWoY8vEkF/X/PZJ2O12TpbDhsfYBIHiXSBQqM+lDhPwFIjyQLjsqq2FtG1vcSj/OUIMMrx3eG7S/RaJv//979wvAvkRaDAUru5zTZs2pWPHjpXZ/8gjj9Cbb75ZZv+aNWsoKyurzP69e/dSmzZtQjZPIR4EwiMHgrOo7ZX2P+QV2mnhmnw6flYdsEurBKqequYfBUKx1f2GzWwKXuWDUI4d6vFjdezEBB0lmvVUOz3guKNK4fdvGT16NA0fPpxbl6amplK42Lp1KyfsaezevZuGDBlCv/rVr8p93/79+91MKvTiFoRrEQgsKWnNg/CRtNkrLxAnz9toweo8ulqgUIKJaNRNKdS6ifemXhVR0VLWtRDKsUM9fqyOXS3FQOlpBsqoFvgNRWUISIrQmc5Xd7pQ4Xlxf/HFF9mB3r9//3LfV6dOHe7DLQj+guYuEIic3BJKSTJUOgeiIoHYfchC63YWsgO8ZnU9jctKpVrp4fnCC0Kg+CV/2dnZAQ168uRJCgWwYj7++GO6//77Odu7PFBHKjMzkwYNGkSrV68u99ji4mJ26rhuQnxaERAIZ89pR4irWwSTp0C41GHSSjv7AmN9u62AVm9XBaJVYxPdO7KaCIQQ/SLRo0cPevDBB2nLli0+j4HH/N133+UQ2QULFlAoQIHBnJwcuvfee30eA2F45513uP8F5tG6dWsWCjRM8sULL7zAXn9ta9SoUUjmL0T+MpMqBOq+MgLhGeIagIM6r8BO81fn057DagRT/y6JNC4rhcwJklQnxEAI7KVLl+j555+nOXPmkMlkou7du1P9+vUpMTGR25r+/PPPtGfPHt4Pp/bNN98ckskOGzaMl7v++9//BvS+UaNGseWxaNEin5YENg1YEhAKCYGNL4GwWhW6nKte+ZPMOqf/4VojmLLPWmnB6nwqKFLIbCIa2iuZOrYwU7yvvYd6/HjwSdQKwHkd0hDYjIwMevnll+lvf/sbff3119yq9OjRo1RYWEi1atWiu+66iy/gsCJCBSKcVq5cWSkrpVevXrxM5Quz2cybEH/lvi9ccZgNCqJK1IewFgqLiYxG3TVFMEFUENqKEFdYIvA/jLwxuVIRTIJQVQTkuIblcNttt/EWbubOncvO6JEjRwb83p07d/IylCB4i2JiSwGZ1I7lJFgUOhP26yotEEiIWrqpgJ3UAP6Hgd2Tgp6MJgihJjyBttcIEvcgEsjTMBrdpzxt2jR2lH/44Yf8fNasWZxb0b59e6ejG/4JbIJQXpgr12Fy+CAgGgo3HQ3c/5BztYT9D2cvlRDiKwZ1T6I2TU0VBlsIQiQSFSKBZSZEWCGqyZPTp0+7RV9BGJ566ikWjqSkJBaLJUuW0IgRI8I8ayHa8iCKih1hrlh6stkp2aAPWCCOnLLSf9bmU2Gxwn6NsQNSqGmmqcw6syuu4qHXBy4kuG9C+C6/X4dWw8FbHw/l2KEeP1bH1ut1zi2iajfFE1K7Kf4EAstDOfl2ruLKrUZ1RCmJer8FAl+jzbuLac2OQh43s6aBbstKpeqp+nKdkRAI7auOL31ljI38QjtngqtjEKUmBc/nEcqxQz1+rI6dmqx3Oq5rpBkip3aTIMSSQHgt9V2i9qLGY10AFoTFqtDi7/Jp31HV892pZQIN75XMju+KrAf812BwPNZVrssELiTarR66CgfzBjOUY4d6/FgdW+fIXQiTISEiIcS3QCAPAmCfpcTudFC72tflCcSl3BL68ts8upCDAphEQ25Ipq6tE8r1P2hioDgEQrso4EtfGUvC4HpR0anjBYtQjh3q8WN1bJ3jcxIuF1elLImPPvqIZs+eTUeOHKFNmzZRkyZN2GHcrFkzbm8qCNEmEK5RTP72oT5w3EKL1uVz6Gxqko6XlxrW8f2V0r7UmvWgfdFZHBwXgtLFJ/8x6NVQXdJEJ4hXj1COHerxY3VsHZYoIR5hUomA9e/tt9+mqVOnsiMY2c9a8T3USYJQCEIsCER5wP+wflchfbFKFYiGdQx0/6hqfgkEcLUYXB9XRiAEIeIsiX/84x9cfmPMmDFcbE8D2daIKhKESBUIO4JaS8inQJT1QZS9aBdZFPrv+nw6cFz1P3RrY6bBPZKcfgXvqHd9dkUhI6JSDKowaEEruHPUqEyQDPI5NG0zG3VkCqKnMZRjh3r8WB3boC/dwkHA/zQsMaF4nifIWM7Pzw/WvAQh6AKB9qMWR9evMqW+/XBSn7+M/Ic8upRr5y/ozb2TqdN15WfqOx3UOmKBcIY16vHTseSk15GmMZVxRuIigrpS6mOdm+hcK6EcO9Tjx+rYeixR4vMTqSIBv8OuXbvYD+HK0qVLqV27dsGcmyBck0BYrI4wV1LIalHIhuUmpTSz2jOLujwfxN6jFlq8IZ9Fp1qKjst7Z9Yq/+ujhbcqji87Yt/x+5xffIdTUguCUn0UlfFJqEID9OwID+b6eOjGDvX4sTq2zuGTCFdyZsAi8fvf/54mT55MRUVFvDaLyrD//ve/uZLqe++9F5pZCoK/ApGr+sggEJwD4SEQthK7ut9P6wEX9rU7CmnTbrUAZJNMI43pn8I5FL5w/fJqDmrcGWKv0SEQmvXg5pOQjGwhAglYJO677z6y2Wz0hz/8gQoKCmjixInUoEEDeu2112j8+PGhmaUg+CkQcDFgKclZZkPLgajQ/1CWgiI7fbUun46cUosA9mxvpqxuSeXeNXpaD9p1n5cQ4JrWqSKhdaSsrPUgCOGiUu4W9JbAduHCBa6rhMJ7ghBLAgH/Awr0Xcmz8wV+5I0p1K6Z766M2oVeXS9WxUGzErApej3Z9Qo/TjCKOAjRwzX55FEmXBCqvNQ3RAGVXG0Vh7j6IxD7jlno222FLDjpaXq6fWAq1alh8Mt60JaXAPwNnB2LuHdkztoVh3CEbz1ZEMIuEohs8vYBxz6UEm/ZsiV3jsvKyrrmyQlCsEJcOZCpAgc16uV890MR/XBALe/dooGRRvdLoSSzvlLWA0cvOfabSWER42M9wmVdV6+0pLpA0atTJpNBLeUQTEI5dqjHj8WxdfrwRjcF/GuGDx9Ohw8fppSUFBaCAQMGUGpqKh06dIjbnKIq6+DBg+mrr74KzYyFuKe8EFcIhBrBpAoERzf5IRB5hXau3qoJRI+2ZvrVoNRyBUK7tmulNVgYOLRVzX1QI1OI/Q+agHjWUtCyrfmLbyz98geyGVzeZ6zkGFUxdjTP3VCFY/PHKJJrN8EP8bvf/Y7+8pe/uO1H1zp0j1uxYgVNnz6dnnvuOSnRIQQdRBtBIHJySyglyeAUCJtVLc5XmQimk+dsNH9NHuUVKOx/GNIzmVo0MPlwUDtyG5wlm11Ka3hYD2Wc064lwbWIJodAqCMLQuQRsEh8/vnntH379jL7EdnUrVs3zsaeMGECvfrqq8GaoyC4WREQCM05rSC89BoEYucvxbRicwH7NGqk6WnEjclcgtkbrsusWoVXLXNaEw1NIIxlQlvV368ZE07rwa1UeGX/0DoyGrTkq2CvQYRy7FCPH7tj69mqiNA8CfgdNm7cyL4HV7APrwFEPEnPaCFUy0xaFJNiVqi42JEPgeQ5m/8OauQtrPi+gHb9UtpeFOGtCSZductLLAoO57Rf1oPrGC6WAzZU+AQQpToZUrVfiEwC/mQ++uij9NBDD7E1AR8EvghIqEMi3dNPP83HLF++3GvpDkG4VoGA81cNc1WooFgtiheogzo3304LVufRqQtq4l3/ronUp2Mi5RUqQbYePLNo1ccJBiTSqa/XSjdQrXQRCCFyqVRnunnz5tEbb7xB+/fv5+etW7dm8UBiHSgsLHRGO0Uj0pkucjOprxYolFeoikWCCRdffUACkX3GSgvW5FNBkUKJCTq6tV8KtWho8to9rlqK4ZqtB21FQNEplO8QofRUPdWuYRRxEKLiuibtS4N4MoXQJ8ohCzoXF3O7QkaTepH2pwYT7oW27S2mVVsLeXkKeQ/jBqa4tX/UREK7zqclGzhr2pf14FpWw1tYuFufCD3R1Tw7t51s1TghbOvJglBl7UstFgudO3eO/Q+uNG7cuLJDCoJPgYDlgCimYgtyIVQZsFsUSkjQVeh/QOb10o0FtPuw6n9o18xEI/qkePU/aJFLbCEYr8168HROQyCwvCQCIUQTAYvEgQMH6P7772dHteedGr4sWhMiQQi2QDjzIDjxTfVLmOy6cgUi5yrKe+fT2UslfEEf1D2JerQz+0wIxSiq1VAqENpjCERlxAGI70GIG5FANrXRaKTFixdTZmamlBcQwiIQWpgrDFdEJmHJqCL/w5FTVk6QKyxWKDlRR2P7p1CTTNX/4IlaeE8pfYxMV31g1oO2tMT+C8d+EQch7kQCvSQQ2dSmTZvQzEiIa8oTCM6DsJfmQUAbfPkfNu8upjU7Cvm4zJrwP6RStZTyy2uoDeh1zvh0raR3ZawHCWsV4lYk0FgIWdeCEA6BsBSTs0EQ8iC4T3U5vaixHLVkQz7tO6a2F+18XQIN65nsDF91xfWCD/+DrURHNhuyrpHMVLq85C2sVZaWhHghYJGYOXMm95J4/vnnqWPHjmQyuZvvEg0kBEMgIAxWCywHVSCcdZjKLDOVXrjx/i+/zaMLOXbOSB16QzJ1aZ3g0/+gXeixPAQrgh3Vjuglo6FUIGRpSYhnAs4nR/G+zZs306BBg7iPRI0aNXhLT0/nn8Hk2WefdWvVh61evXrlvmft2rVcHgQ5Gs2bN6fZs2cHdU5C6ATCVlIqEBbuJqdwZVYIBFsQjjwIiIPnMtOB4xaa+99cFojUJB39engadW1T1kHtWqabxcGR/8A+B0erSO2xp0BwBVeXkhqa7wFLS22amiXvQYhJArYkVq9eTeGkffv2tHLlSudzQzk1eY8cOUIjRozghkgff/wxfffdd/TII49Q7dq1ady4cWGasVCZfhAsEA6HNASC8yIcIlFekhz8D+t3FdGGH4r4ecM6BrptQCqlJusrLK2hVdTUopfwxMimikMIPARCHNNCPBKwSPTv35/CCSKpKrIeNGA1IE9j1qxZ/Lxt27a0bds2evnll0UkIrjct9qPWhUGq9V/gUDOBPIfDp5Q/Q/d2phpcI8kFgB/SmuoloGjpDeEwK6QzaEirr0iNAHRxsV/JWpJiBcqnUyH/tbZ2dmcVOdKp06dKJggL6N+/fpcMLBnz57sC8Eykjc2bdpEQ4cOdds3bNgwev/998lqtZbxn2gUFxfz5pqZKISnzEaxxZEDUQKxcK/kWp5AXLxSQl9/V0A5eXb2J9zcJ5k6tTQHbD24Ri+VsIvDXUQkakmIdwIWifPnz9N9991HS5cu9fp6MJPpIAoffvghtWrVis6ePcs9K/r06UN79uyhmjVrljn+zJkzVLduXbd9eG6z2TgiC3kd3njhhRdoxowZQZu3cA0hrn4IxMHjVlq5tYAbDiGsdVxWCmXW8vwol7YRdfU7wBhgp7WhbFkNtXiA4latVXIehHgnYMf1E088QZcvX2bndVJSEi1btoz+9a9/0XXXXUeLFi0K6uRuvvlmXiZCFBUc5kuWLOH9+H2+8HRUavULy+spPG3aNK5nom3Hjx8P2r9BqDjE1WIrDXGFb8Luw0GNhkPf/VhESzepAgH/w32j0soIhOqcVgUAEUvc7QvRS7AkDHoWCFgPnnWXXK0HI2ozORzTWFoSx7QQrwRsSXz77bfcmhRlwvV6PTVp0oSGDBnCoa+4Ix85cmRoZkrELVMhGFiC8gZ8F7AmXEF9Kfg1vFkeGljKkv4XVSMQaogrIpjK9z+gsB+yp4+etvFzhLaivHdKot5nYpwWqYS+Lb6sBw0OhzW4v6e2lPEWhMAtifz8fA59BRkZGbz8BHDx3rFjR0hPKfwGe/fu9bls1Lt3b/rmm2/c9qGdavfu3X36I4SqyYHQBEINcS1fIM5ctNHc/15lgUCS27BeSXRT5yS3QnlaaKtmPThzHyqwHoBmPWjLUqgMK5aDIFRSJNA7Qusjcf3119M///lPOnnyJEcW+bp4V5annnqK8x4Q2vr999/T7bffzk7le+65x7lMNGnSJOfxaIaEPttTp05lMZkzZw47rTGOEHlJcv5EMP10qJg+/PoqXcm3c3vRXw1K5VLbvqwH7UIPSwDOaYgDHmuJcaqlUBq5pDUD0nImIBBSqVUQrmG5CT6J06dP8+Pp06dz9BCaECUkJNAHH3xAweTEiRPcLxtOZ+Q69OrVi30hWOICmAcirDSaNWtGX3/9NT355JP05ptvclTU66+/LuGvESAQKNcNh7S/AoHX0PsBPSBAi4ZGGt03hX0WGppFoPdoCqSJA+6AWAS8LS05nNMS1ioIIW46hFDYffv2cX5CrVq1KBaQpkNBzqIuKZsDgX2+ivTlFdpp4Zp8On5W9T/c2DmR+l2fyBd6z8ZA6amGgH0PnklxUoxPiAdyw910SCM5OZm6du16rcMIMQSikCAQyGVISdKzUEAg7I4QV6c14EUgTp6z0fw1eZRXoHB7UlgPrstL2vVeK5HBuQ9+Wg+SFCcIgROwSCAPAstKq1at8tqZDtFPQnwDKwICgY+Ga4mNihzUO/cX04rvC/j4mtX1dPvAVKpZ3b0MC8RAGwDZ0xVZD7K0JAhhFonHH3+cRQKhrh06dJCmQ0IZgUAXOIgBymZgPQeX6/IEAktPKzYX0K4DavZ+6yYmuuWmFDJ7tBfVGgNhHDzmng9GffnWgywtCUJ4ReLTTz+lzz//nAvpCYI3CwLlNYogEOiFblU4mc2XQOTm22nB6jw6dUFdgxrQNZF6d1T9D76jl9TwVoTDemsIJEtLglCFIoEoppYtWwZxCkKsCMT5HNWCKLCoTmksMZENnxmdV/9D9hkrLViTTwVFCiUm6GhM/xRq3sA9n8VTIAx2iIDi7PkgS0uCEGF5Er/73e/otddec5a7EOIblPved7TIKRBYVkKYqxriqj73FAh8drb+XESfLM9jgahTQy2v4SoQWnKc1gwIAoHEODOX2ShNlHOzHjzqLdWUPg+CEB5L4rbbbivjnEaBP/R68MxkXrBgwbXPSojaPAhuFuQIecVP1c1cKhA4DuW9dx9W/Q/tmyfQiD7JLAS+rAegZk7rCJ82CBBHNzkimyRqSRCqWCQQW+vK2LFjQzUfIUoFQu0qpwoEh7k6i/SVknO1hOavzmfHNi7sg7onUY92pd3j3DKhnZnTsCLUvg8QDRgmWjkOtdS3+hpbEdLnQRCqRiTmzp0b/N8sxKhAoO2oZw6Ejo6csnKBvsJihZITdTS2fwo1yXRfXvJlPXhGL+lQztshJNqxUoxPECLEcY06SujPgNLgrqAyK5aemjZtGsz5CRGeSY2uciW8vKQ4S367CgT8D5t+KqI1Owr5tcxaBhqXlcp9IPy1Hlyd00aD4rQeTAbV71An45pzQgVB8EHAjut7772XNm7cWGY/CvDhNSF2CVQgEP66bFMBrd6uCkTn6xLo7uFpZQSCw1mNqkDAejDAOe1FIFhIHK9hV+10owiEIISYgG/Bdu7cSTfeeGOZ/Si+N2XKlGDNS4hQgUCzHwgEfmKpyVcdpstX1fail3LtnCcxtGcydWmV4OZ/0JzOsCDKsx6047hGk56oRjWDLC8JQqSKBBdZu3q1zH4UjQpm61IhcgUCVgOX/UapDZszQ84pEPA/oLwGLI2UJB0vLzWsYwzY9+AtckmWlwQhwpeb+vbtyx3oXAUBj7HvpptuCvb8hAgRCAgDch74J1dzdQgE/q84Go0qCq3bWUiLN6gCAf/DnYPdBULr4aBZDyi9gaQ4kyODWuv3oOU9cE6EI+8BzmnxPwhChFsSL730EvXr14+bD0EwwPr167kMrRT3i80salUQ1EqurjkQrtZDUbGdFq0voIMnrPy8U8sEuqlzorOMRtmmQKUZ05zzoPe+tAQkckkQosiSaNeuHf344490xx13cBVYLD2hOxx6SqDgnxB7AgGrQBMIbXMViPOXS+iDJVdZIGANDL4hifp3TfJoBFRqPSSgequjvajWi9qtSxyc2Hp1aaltUzPVSpfoJUGI2qZDsUi8Nx1yFYiiYrU/hOaL8HRQ7z1qocUb8tnSQNTS7QNTuIeEVp8PF/70NGPF1oMkxQlCbDYdEuJTILAfuQ+bd6vtRZtmGrlAX3KinjvLOZeN8CHzEbkkS0uCEPmISAhuAnHucgnZFYUzo+0lagRTsbU0MQ4UFNk5e/roabW9aK8OZhrQNclhGajWAhRFr9NTghl9H1TrQRUIiVoShGhCREJwEwi2INBNDmGujrIbrstLZy7aaP63+XQl304mI9HIG1OoXTO1vagmAOgWp5CeRSHB6CWsVZaWBCH2RCI7O5saNWoknehiuZucogqE1o/aUyB+OlTMFVyRcV0jTU/jBqZymW+glfWGGJjNes7Exr4Ek97pmJaoJUGIYZFo1qwZnT59murUqRPaGQlh7QWBTfM3+BIIWBerthbStr2q/6FFQyPd2jeFEs36MuGtaskMhQrtjnIbhtIIJglpFYQYFgkJgoot4HiGQCB8Fc5mG3IhXPIgNIHIK7DTwrX5dPys6n9A7kPf69X2oq7Wg9Gg5+UniISCPAoU4sPyEoTDIRCSLS0I0Yf4JOI5iukySm2owoCfmkWh3RCcPGej+WvyKK9AIbOJaFTfFGrVOMGr9QDfgyYaCsQBS0wIeXUIhCTECUIciMR7771Hqamp5R7z2GOPXeuchDAIxJlL6jITWo4WWe1ciwloArFzfzEt/76AX69ZHfkPqVSzusGn9aDtM3BZDYyp/qybYeRNEIQYT6bT6/XUsGFDMiDg3ddgOh0dPnyYop1YTqbT8iDQJS6vUK23hKUlhKlieQlWxYrNBbTrgNpetHUTE91yUwolJujLtR6czmmXyCVZXhKEOEum27ZtW1gd1ygaiJ7ZKPmRlJREffr0oZkzZ3LdKF+sWbOGsrKyyuzfu3cvtWnThuIZ10S5gmJVIGApYKnJaCTKzbfTgtV5dOqCWrwRuQ+9O5r5BsEf60Gc04IQe/gtEto6dDhZu3YtTZ48mXr06MHd8J555hkaOnQo/fzzz5SSklLue/fv3++mlrVr16Z4RhMIWAqcB8E1mNRCfRCLY2dstHBNPhUUKZSYoOPs6eYNTGWWlzytB+0nEN+DIMQeER3dtGzZsjK9tmHJbN++nSvRlgeOS09PD/EMozhRziEQdrudfjhgoQ0/FLHTGnkPqL9Uo5oxIOsho5pBfA+CEM8iMX369Aqd1qEGa2kgIyOjwmO7dOlCRUVFXLX2z3/+s9clKI3i4mLeXNfuYk0gtOglG/pBOCKaLFY7rd5WSPuz1fLe7Zsn0Ig+yWRO0DuT3xKdj1FuQ5aWBCHe8NtxfenSJSooKGDntcaePXvo5Zdfpvz8fBozZgxNnDgxZBPFNG+99Va6fPky968ob5lp3bp11K1bN77wf/TRRzR79mz2VfiyPp599lmaMWNGmf3R7rh2RjFxkyC13LdWh+n4WSt9vTGfLuTYWQwG9UiiG9olOpeQ2HpwNAGSpSVBiF/Htd8iMWHCBMrMzKRXX32Vn6OXBBzB9evXpxYtWtDSpUvp/fffp7vvvptCAXwTS5YsoQ0bNrgJlT+MGjWK74QXLVrktyWBEiTRLBKuAsFLSzaHQChEh05auEAflp6SzDq6uXcytW1mdi4loZQG6i+5Wg9a+1DxOwhCdBLy6KbNmzezT0Djww8/5GWfXbt2kdFoZIvizTffDIlIPProo3yBh4UQqECAXr160ccff+zzdbPZzFus4CoQ6PMAJzUEAv6HTT8V09qdhex/qJthoBE3JnMfCIiBp/Xg7P0gfgdBiFv87kx35swZrt+kgValY8eOZYEAo0ePpgMHDgR1cjBypkyZwmGw+H2uvz8Qdu7cyVZQPOBLIIosdlqwJp97QEAg2jUz0biBKZSeqieDw3pA5BJnSRvUvAnOh3DJmJakOEGIP/y2JGCe5OTkUJMmTfj5li1b6IEHHnC+jqUJ1yWbYC0xffLJJ/TVV19RWloaCxWAyYS8CTBt2jQ6efIkWzZg1qxZ1LRpU2rfvj1ZLBa2IObPn89bvAoEajTNX53H/ge9nmhYz2Rq3TSB9CjordeRyaQjk6OdKLcVlaUlQRACFYkbbriBXn/9dXr33Xf5zh69rQcOHOh8/ZdffuF1/GDy9ttv888BAwa47cey17333suPUZkWZcw1IAxPPfUUCweEBGIBX8aIESMoljlzEeW+bRy5BIFQfyr0S3YxLVqXT8VWotRkHY0bkEoN65rIarWTQqqlYDap2dNavSUJaRUEIWDHNXwPgwcPZnFAYtvTTz9Nzz33nPN1+CKQ4IZIomgnmspywEqAODiL9FnJGd66bmch5z+ARnWNNLZ/ClVLNbDFgExrHAdRSE/RU1Kioy+EFOMThJgkN9SO6+uvv55LW2zcuJHq1atHPXv2dHt9/PjxnJMghD+LGstLiFSC3EMc8gpKaNH6Ajp4Qs1/6N7WTINvSFb9DHr4G/SkoGkEGgRhn0kv4iAIwrVZEvFENFgSWpIcHNJYXkJ/CAgFrIr53+bR5at2dkAjvLXTdYn8GBnT8D+gfLdep1CRqiHUtL6J6tcyVfU/SRCEaLYkNMdwRUyaNMnvXy4E10H985FiWrwhn/dVT1Hbizaoo5bX4DaiqOJqUH0QKNpXZC2hGtUMIhCCIASnVDjKciDk1ddbEOGEzOxoJ5ItCW8CUVhsp9XbC2jzbjW6rGmmkcYOSHXkP6jRSrAe4Jw2IktOfA+CEHfkhtqSaNu2LZ09e5Z+/etf0/3330+dOnWq7FyFIArE5asltHBNHh09rXYN6t0hkQb2SKIEo56XliAJCVhicsl5kGJ8giAEPZkOdZoQSlpYWMg1kLp3784hqrFUDC8aBKLYUioQR09bac6iXBYIVGkdOyCFhvZK5qJ8mu8BZTd4qUmS4gRBCJfjGkLxxRdfcL4CkupQ3G/OnDkxU9oi0pabIBCnL6oCwYX6bApt31tISzcVcF2mGml6umNwKvsWEMGEqq2wHiAOQMJaBUHIDXWBP2+glhJKiOPnhQsXqEaNGjHxl4gkkYBAnDyvJslBIAqKSmj55gLatlf1P7RsaKJxWSmUmmxg3wMyp9E0CP4hEQdBEMLavhQgk/lf//oXWxEoEQ4fBZadYkUgIgkIRPZZKye+QSDOXbLRgjV5dPys6n/o1yWRsrqp+Q9YWjIn6DgHAtSUWkuCIAQBv0Xi888/Z2FAS9Fhw4bRK6+8QiNHjiSDQc3UFYJLscXOPofc/BJKTjTQwWwLfbnmKuUVKBzCCv9D++ZmsR4EQYicENjGjRvTXXfdRXXr1vV53GOPPUbRTlUvN6Gk994jFjp90cpZ1Dv2F9M3WwrYoqiVbqAJQ9Kobk2D6phO1HO/B1Anw0C10gM2DgVBiANyQ+2TQGVVrHOXO5hOR4cPH6ZopypFAgJx8pyNC/Oh3eiyTQX00yELv9a2qYluH5hGSWY9O6YRxQRkaUkQhCr3SRw9etTvQYVrEwhEMuXk2bm895mLJWp70e7J1L9rIuc/JJpV60Ec04IghBpZm4hAgTh0wkrzludSQZHCkUqj+iZTl1ZJZDbruSAfkKUlQRAiSiSQG7Fq1Sq65ZZbnM1+XJsMwYGN0uGJiYmhmWmMCwQiltAT4rsfimjJxnz2RdSpYaAx/ZOpdg0jpSbruUCfLC0JghBOAirwt3jxYqdIvPHGG9zQR+sQt2/fPqpfvz49+eSToZttDHL+spWjmJAot2B1Hu06oApvh+YJNLRXEvsdqqcaOYoJLUTFMS0IQkSKxLx588oIAFqLNm/enB+jTeibb74pIhFgHgQE4lJuCX20NJdOX1D9Dyjv3btjIneTQ7hrg9omyqwlK4OCIIQfv688aE/aqlUr53MsKyEs1rW9KXpSC/5x9qKVjp2x0i/ZFvr3iqtUWKxQSqKOxg9NoxYNTGQ2G6hGmiwvCYIQJSKBsCmUCdc4f/58mXV1Vx+FUL5AHD1tobU7Cmn59wXsf2hYx0gTh6ZSrXQTJSdKpzhBEKJMJBo2bEi7d++m1q1be339xx9/5GME30BIz14qoYPHi+mLVXm0+7Ca/9CtjZn7P6DPNKKXJHJJEISoE4kRI0bQX//6Vy7F4RnBhMinGTNm8GtCOSGu523048Fi+ujrXG49io5xo/ulUO8OyWw91K5h4EgmQRCESMHvjGs0HLr++uspISGBpkyZwv4JZFgjqgmRTjabjXbu3FluyY54zbiGQGSfsdGaHfn02co8jmRC17hfD0+jlo3MVK+mgepmGCvMaBcEQYjYjGtc/Ddu3EgPP/ww/elPf3K2MMWFbciQIfTWW2/FhEAEGwjEwRMW+mLlVVq1rdDZXvTum6uzMIj1IAhCJBPQ2kazZs1o2bJl3Mf64MGDvK9ly5aUkZERqvlFNVabnZeX3v8qh/Yds/K+Pp0SaVxWGjWoDYEQ60EQhMimUgvgEAWEvAq+KbLYaeOPhfT2/Mt08YqdjAai27JSaXTfNF5ecg0fFgRBiFSi4kqFpSxYMXCYd+vWjdavX1/u8eh5geNwPJL9Zs+eTeGk2GqnL7+9Si99dJEFIj1NT0+Mr0Fj+6dxUpwIhCAI0ULEi8Rnn31GTzzxBD3zzDPsGO/bty/dfPPNlJ2d7fX4I0eOcCQWjsPxTz/9NPe4mD9/fljma7Ha6dVPLtGcRVfIYlXbi86cXJuG9UqhujVN4pwWBCGquKYe1+GgZ8+e1LVrV26RqtG2bVsaM2YMvfDCC2WO/+Mf/0iLFi2ivXv3Ovc99NBD9MMPP9CmTZtCGgVwJa+EZrx3gXb9oiYV3nJjCj18e3VKTDCIOAiCUKVU9roW0ZaExWKh7du309ChQ9324zkirbwBIfA8Hu1Wt23bRlar6jz2BJniOIGuW6CU2BWaOuscCwR6TT8xPp0en5BBSWZxTguCEL1EtEhcuHCBSkpKyoTW4vmZM2e8vgf7vR2PPA6M5w1YJFBYbWvUqFHAc0UToHtGVufyGm/9oS6N7lfN2VZUEAQhWolokdDwTDLDCll5iWfejve2XwO9MWCCadvx48crNc9+XZLpvWcyqVn9hEq9XxAEIdKI6BoQtWrV4mZGnlbDuXPnfCbu1atXz+vxKE5Ys2ZNr+8xm828BQP0nhYEQYgVItqSQAkQhLJ+8803bvvxvE+fPl7f07t37zLHr1ixgrp3704mkymk8xUEQYg1IlokwNSpU+m9996jOXPmcMQSGh8h/BURS9pS0aRJk5zHY/+xY8f4fTge73v//ffpqaeeqsJ/hSAIQnQS0ctN4M4776SLFy/S//zP/9Dp06epQ4cO9PXXX1OTJk34dexzzZlA0h1eh5igUx5aqr7++us0bty4KvxXCIIgRCcRnydRFcB5nZ6ezg7sYFSBFQRBqGoQ2o/IzZycHI7ijBlLoiq4evUq/6xMKKwgCEKkX98CEQmxJHyU9z516hSlpaX5DJvVVDlarY1onn80zx3I/OX8V8XnB4tGEAgswQdSP04sCS/gBPrbihV/pGi8UMXC/KN57kDmL+c/3J+fQCyIqIluEgRBEKoOEQlBEATBJyISlQQZ2tOnTw9apna4ieb5R/Pcgcxfzn80fX7EcS0IgiD4RCwJQRAEwSciEoIgCIJPRCQEQRAEn4hICIIgCD6JW5FAN7oePXpwVnWdOnW4Z/b+/fvLZCg+++yznKGYlJREAwYMoD179pRpffroo49y74uUlBQaPXo0nThxwu2Yy5cv09133+3sfIfHqJ8SCfPHPmSVu27jx4+PiPkvWLCAW8/i3GJeu3btKjNOJJ9/f+YfqecfrX7RL75jx458XvEZQrVlVCKIhvPv7/yr4vy/4MdnB9/bNm3a8Nxr1KhBgwcPpu+//75qzr0SpwwbNkyZO3eusnv3bmXXrl3KyJEjlcaNGyt5eXnOY1588UUlLS1NmT9/vvLTTz8pd955p5KZmank5uY6j3nooYeUBg0aKN98842yY8cOJSsrS+ncubNis9mcxwwfPlzp0KGDsnHjRt7w+JZbbomI+ffv31958MEHldOnTzu3nJwct99VVfP/8MMPlRkzZijvvvsuilAqO3fuLDNOJJ9/f+Yfqecfcxg8eLDy2WefKfv27VM2bdqk9OzZU+nWrVtUnH9/518V53+YH5+defPm8Tk9dOgQH/fAAw8o1apVU86dOxf2cx+3IuEJTj6+yGvXruXndrtdqVevHl9oNYqKipTq1asrs2fP5uf4MJlMJuXTTz91HnPy5ElFr9cry5Yt4+c///wzj7t582bnMfjAYh8+vFU5f+1L8vjjj/sct6rm78qRI0e8XmQj+fz7M/9oOf8aW7Zs4WOOHTsWVeff1/wj5fyf82PuV65c4WNWrlwZ9nMft8tN3sqDg4yMDP555MgRboM6dOhQ5zFIXunfvz9t3LiRn2/fvp3NWtdjYNai54V2zKZNm9jM69mzp/OYXr168T7tmKqav8a8efPYZG3fvj03Z9Kq4Fbl/P0hks9/IETL+ccxWI5BGf1oPP+e84+U83+lgrlbLBZ65513+Hd27tw57OdeCvw51u7Rye6mm27ikwy0PtmevbTxHJ3vtGPQYhVrhp7HaO/HT6w7eoJ9nr24wz1/cNddd3GjJvQG3717N3f6++GHH5wtYKtq/v4QyeffX6Ll/BcVFdGf/vQnmjhxorOoXDSdf2/zj4Tzr5Qz98WLF7N/pKCggDIzM3lOELNwn3sRCSKaMmUK/fjjj7Rhw4YyJ8izVDj+qL7Kh/s6xtvx/owTjvk/+OCDzsf4kF533XXcD3zHjh3UtWvXKp9/ZYik818R0XD+cceKixVK6L/11ltRd/7Lm39Vn/8p5cw9KyuLgx0uXLhA7777Lt1xxx3svPZ24fc1r2DMPe6XmxAdsGjRIlq9erVbeXDcWQBPxT137pzz7hzHwBREBEF5x5w9e7bMiT9//nyZu/xwz98b+GKYTCY6cOBAlc7fHyL5/FeWSDv/uMDi4oTlS9zJut6FR8P5L2/+VX3+H61g7ohYatmyJS8Rvf/++2Q0Gvln2M+9EqfAsTt58mSlfv36yi+//OL1dTh+Z86c6dxXXFzs1XGNCAqNU6dOeXUeff/9985j4Ei6VsdXMObvDURBuTrRqmr+gTiuI/H8+zP/SD//FotFGTNmjNK+fXu3qJpoOf8Vzb+qzr89gM+OKy1atFCmT58e9nMftyLx8MMP8wVzzZo1buFvBQUFzmMQGYRjFixYwB+eCRMmeA2BbdiwIUcdIAxt4MCBXsPQOnXqxJEF2Dp27HjNIYDBmP/Bgwc5RHPr1q18IVuyZInSpk0bpUuXLhEx/4sXL/KFFfPCBxuRHHiO46Lh/Fc0/0g+/1arVRk9ejSfW4Rpuh6Dm41IP//+zL+qzv/DFcwdobDTpk3j33X06FFl+/btHAJrNps5HDbc5z5uRQJfWm8b4pddFR/KjTty/IH69evHF1tXCgsLlSlTpigZGRlKUlIS/wGys7PdjsHF4q677uKcBWx4fPny5SqfP+aJfZh7QkIC36k89thjPN9ImD8eeztGu5uK9PNf0fwj+fxr1o+3bfXq1RF//v2Zf1Wdf6pg7jinY8eOZUsD88KNHQQPIbyuhOvcS6lwQRAEwSdx77gWBEEQfCMiIQiCIPhEREIQBEHwiYiEIAiC4BMRCUEQBMEnIhKCIAiCT0QkBEEQBJ+ISAiCIAg+EZEQBEEQfCIiIQhBAhUX0IsYfa09QYlqNHvJzs6W8y1EFSISghAkUKN/7ty5XPP/n//8p3M/ylT/8Y9/pNdee40aN24c1PONUtiCEEpEJAQhiDRq1IjFAG0wIQ6wLh544AEaNGgQ3XDDDTRixAhKTU3lev533303N5TRWLZsGXcoQ3vNmjVr0i233EKHDh1yvn706FEWos8//5wGDBhAiYmJ9PHHH8vfTwgpUuBPEELAmDFjKCcnh8aNG0fPPfccbd26lTueoRPapEmTqLCwkK0Lm81G3377Lb9n/vz5LAIdO3ak/Px8+utf/8rCgO5ker2eH6PVZtOmTemVV16hLl26cN9y9DYWhFAhIiEIIQAdwtAO8+LFi/Tll1/Szp07eRlq+fLlzmNOnDjBlsf+/fupVatWXjuIoVXlTz/9xGNpIjFr1ix6/PHH5e8mhAVZbhKEEICL+29+8xtq27YtjR07lrZv385tKrHUpG1t2rThY7UlJfycOHEiNW/enNtsQhCAp7MbFokghAtj2H6TIMQZ6EmMDdjtdho1ahTNnDmzzHGZmZn8E6/DskDTeywh4T2wINDL2LP3sSCECxEJQQgDXbt2ZZ8D/AmacLiCZam9e/dyVFTfvn1534YNG+RvI1Q5stwkCGFg8uTJdOnSJZowYQJt2bKFDh8+TCtWrKD777+fSkpKqEaNGhzR9M4779DBgwfZmT116lT52whVjoiEIIQBLB999913LAhItsMyEpzPSLBD5BK2Tz/9lH0XeO3JJ5+kv//97/K3EaociW4SBEEQfCKWhCAIguATEQlBEATBJyISgiAIgk9EJARBEASfiEgIgiAIPhGREARBEHwiIiEIgiD4RERCEARB8ImIhCAIguATEQlBEATBJyISgiAIAvni/wG7AXnvwTZiYQAAAABJRU5ErkJggg==", "text/plain": [ "
" ] @@ -123,10 +123,22 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": null, "id": "a1c9b208", "metadata": {}, - "outputs": [], + "outputs": [ + { + "ename": "TypeError", + "evalue": "AntarcticaISMIP6.__init__() got an unexpected keyword argument 'ais_region'", + "output_type": "error", + "traceback": [ + "\u001b[31m---------------------------------------------------------------------------\u001b[39m", + "\u001b[31mTypeError\u001b[39m Traceback (most recent call last)", + "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[3]\u001b[39m\u001b[32m, line 7\u001b[39m\n\u001b[32m 3\u001b[39m \u001b[33m\"greenland\"\u001b[39m: GreenlandAR6(),\n\u001b[32m 4\u001b[39m \"expansion\": ThermalExpansion(\n\u001b[32m 5\u001b[39m ohc_change\n\u001b[32m 6\u001b[39m ), \u001b[38;5;66;03m# only the thermal expansion needs OHC change, so we'll pass it in here\u001b[39;00m\n\u001b[32m----> \u001b[39m\u001b[32m7\u001b[39m \u001b[33m\"wais\"\u001b[39m: AntarcticaISMIP6(ais_region=\u001b[33m\"wais\"\u001b[39m),\n\u001b[32m 8\u001b[39m \u001b[33m\"eais\"\u001b[39m: AntarcticaISMIP6(ais_region=\u001b[33m\"eais\"\u001b[39m),\n\u001b[32m 9\u001b[39m \u001b[33m\"peninsula\"\u001b[39m: AntarcticaISMIP6(ais_region=\u001b[33m\"peninsula\"\u001b[39m),\n\u001b[32m 10\u001b[39m \u001b[33m\"glacier\"\u001b[39m: Glacier(),\n", + "\u001b[31mTypeError\u001b[39m: AntarcticaISMIP6.__init__() got an unexpected keyword argument 'ais_region'" + ] + } + ], "source": [ "global_components = {\n", " \"landwater\": LandwaterAR6(),\n", @@ -134,9 +146,9 @@ " \"expansion\": ThermalExpansion(\n", " ohc_change\n", " ), # only the thermal expansion needs OHC change, so we'll pass it in here\n", - " \"wais\": AntarcticaISMIP6(ais_region=\"wais\"),\n", - " \"eais\": AntarcticaISMIP6(ais_region=\"eais\"),\n", - " \"peninsula\": AntarcticaISMIP6(ais_region=\"peninsula\"),\n", + " \"wais\": AntarcticaISMIP6(region=\"wais\"),\n", + " \"eais\": AntarcticaISMIP6(region=\"eais\"),\n", + " \"peninsula\": AntarcticaISMIP6(region=\"peninsula\"),\n", " \"glacier\": Glacier(),\n", "}\n", "global_model = Global(components=global_components, end_yr=2301, num_members=1000)" @@ -152,7 +164,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": null, "id": "84894d79", "metadata": {}, "outputs": [ @@ -215,7 +227,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": null, "id": "5b63dfb5", "metadata": {}, "outputs": [ @@ -256,7 +268,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": null, "id": "47694def", "metadata": {}, "outputs": [], @@ -279,7 +291,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": null, "id": "09f1fff4", "metadata": {}, "outputs": [ @@ -445,7 +457,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": null, "id": "60d81a5b", "metadata": {}, "outputs": [ @@ -502,7 +514,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": null, "id": "e62d28d9", "metadata": {}, "outputs": [ @@ -559,7 +571,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": null, "id": "5512d606", "metadata": {}, "outputs": [ @@ -720,7 +732,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": null, "id": "061f92d5", "metadata": {}, "outputs": [ diff --git a/profsea/utils/utils.py b/profsea/utils/utils.py index c4c680e..848fd0f 100644 --- a/profsea/utils/utils.py +++ b/profsea/utils/utils.py @@ -1,9 +1,11 @@ from __future__ import annotations +import logging import os import zipfile from pathlib import Path +import dask import dask.array as da import numpy as np import requests @@ -20,29 +22,49 @@ from scipy.spatial.distance import cdist console = Console() +logger = logging.getLogger(__name__) -def sample_members_2D(array: np.ndarray, percentiles: list | np.ndarray) -> np.ndarray: +def sample_members_2D( + array: np.ndarray | da.Array, + percentiles: list | np.ndarray, + dtype: np.dtype = np.float32, +) -> np.ndarray | da.Array: """ - Sample real ensemble members from a 2D numpy array. + Sample real ensemble members from a 2D numpy or dask array lazily. Parameters ---------- - array: np.ndarray + array: np.ndarray | da.Array Input 2D array of shape (realisation, time). percentiles: list | np.ndarray List of percentiles to sample from the input array. Returns ------- - np.ndarray - Sampled array of shape (len(percentiles), time) corresponding to the closest real ensemble members to the specified percentiles. + np.ndarray | da.Array + Sampled array of shape (len(percentiles), time). Maintains lazy + evaluation if a dask array is provided. """ - # Caculate statistical timeseries, then match with closest real timeseries - array_percentiles = np.nanpercentile(array, percentiles, axis=0) - distances = cdist(array_percentiles, array) - mem_indices = np.argmin(distances, axis=1) - return array[mem_indices] + + def _eager_sample(arr, percs): + # Calculate statistical timeseries, then match with closest real timeseries + array_percentiles = np.nanpercentile(arr, percs, axis=0) + distances = cdist(array_percentiles, arr) + mem_indices = np.argmin(distances, axis=1) + return arr[mem_indices].astype(dtype) + + if isinstance(array, da.Array): + # Tell Dask to delay this operation until the graph is computed + lazy_result = dask.delayed(_eager_sample)(array, percentiles) + + # Reconstruct into a Dask array so downstream Xarray operations continue to work lazily + shape = (len(percentiles), array.shape[1]) + return da.from_delayed(lazy_result, shape=shape, dtype=array.dtype) + + else: + # Fallback for standard numpy arrays + return _eager_sample(array, percentiles) def interpolate(data: da.array, lats: int, lons: int) -> da.array: @@ -185,14 +207,14 @@ def fetch_zenodo_fingerprints( # 1. Check if data already exists if target_dir.exists() and any(target_dir.iterdir()): - console.log("[bold green]✓ ProFSea assets found locally![/bold green]") + logger.info("[bold green]✓ ProFSea assets found locally![/bold green]") return # Create the base directory if it doesn't exist data_dir.mkdir(parents=True, exist_ok=True) zip_path = data_dir / "temp_fingerprints.zip" - console.log(f"Initiating download from {zenodo_url}...") + logger.info(f"Initiating download from {zenodo_url}...") # 2. Stream the download with a rich progress bar try: @@ -220,12 +242,12 @@ def fetch_zenodo_fingerprints( progress.update(download_task, advance=len(chunk)) except requests.exceptions.RequestException as e: - console.log(f"[bold red]Failed to download data: {e}[/bold red]") + logger.error(f"[bold red]Failed to download data: {e}[/bold red]") if zip_path.exists(): zip_path.unlink() # Clean up partial downloads raise - console.log("Extracting data...") + logger.info("Extracting data...") try: with zipfile.ZipFile(zip_path, "r") as zip_ref: # Filter out the __MACOSX directory and its contents @@ -236,11 +258,11 @@ def fetch_zenodo_fingerprints( ] zip_ref.extractall(data_dir, members=valid_members) - console.log( + logger.info( f"[bold green]✓ Successfully extracted data to {data_dir}[/bold green]" ) except zipfile.BadZipFile: - console.log( + logger.error( "[bold red]Error: Downloaded file is not a valid zip archive.[/bold red]" ) raise @@ -257,29 +279,45 @@ def save_components( output_prefix: str = "projection", output_dir: str = ".", output_format: str = "zarr", + output_dtype: str | None = None, + description: str = "Spatial sea level rise projections", ) -> None: """ - Stream all regional sea level projections to disk in a single file/store. + Stream sea level projections to disk in a single file/store. Parameters ---------- components: dict[str, xr.DataArray] Dictionary of component names and their corresponding Xarray DataArrays. - output_format: str - Format to save the output in. Must be either 'netcdf' or 'zarr'. - output_dir: str - Directory to save components to. scenario_name: str Name of the scenario you've run the emulator for. output_prefix: str - Prefix for the output file name (e.g., 'projection' will result in 'ssp + Prefix for the output file name (e.g., 'projection' or 'global'). + output_dir: str + Directory to save components to. + output_format: str + Format to save the output in. Must be either 'netcdf' or 'zarr'. + output_dtype: str | None + Data type to force the output to (e.g., 'float32'). If None, the dtype + is dynamically inferred from each individual component. + description: str + Metadata description attached to the dataset attributes. Returns ------- None """ + if not components: + logger.warning("No components provided to save.") + return + ds = xr.Dataset(components) + # Add ProFSea version and scenario metadata + ds.attrs["source"] = "ProFSea v3.0" + ds.attrs["scenario"] = scenario_name + ds.attrs["description"] = description + output_format = output_format.lower() if output_format not in ["netcdf", "zarr"]: raise ValueError("output_format must be either 'netcdf' or 'zarr'.") @@ -294,32 +332,45 @@ def save_components( compressor = Blosc(cname="zstd", clevel=5, shuffle=numcodecs.Blosc.BITSHUFFLE) - # Set the encoding/compression for each variable based on the output format + # Set the encoding dynamically, overriding with output_dtype if provided for name, component in components.items(): + comp_dtype = output_dtype if output_dtype else component.dtype.name + if output_format == "netcdf": - encoding[name] = {"zlib": True, "complevel": 1, "dtype": "float32"} + encoding[name] = {"zlib": True, "complevel": 1, "dtype": comp_dtype} elif output_format == "zarr": - encoding[name] = {"compressor": compressor, "dtype": "float32"} + encoding[name] = {"compressor": compressor, "dtype": comp_dtype} file_name = f"{scenario_name}_{output_prefix}" + # Dynamically adjust console message string based on the prefix + log_prefix = "Global " if "global" in output_prefix.lower() else "" + # Stream the computation and write to disk if output_format == "netcdf": out_path = os.path.join(output_dir, f"{file_name}.nc") with console.status( - "[bold cyan]Computing and saving NetCDF...[/bold cyan]", spinner="dots" + f"[bold cyan]Computing and saving {log_prefix}NetCDF...[/bold cyan]", + spinner="dots", ): + # Let xarray handle streaming directly to avoid RAM overload ds.compute().to_netcdf(out_path, encoding=encoding) - console.log(f"[bold green]✓ Successfully saved NetCDF:[/bold green] {out_path}") + + logger.info(f"[bold green]✓ Successfully saved NetCDF:[/bold green] {out_path}") elif output_format == "zarr": out_path = os.path.join(output_dir, f"{file_name}.zarr") with console.status( - "[bold cyan]Streaming computation and saving Zarr...[/bold cyan]", + f"[bold cyan]Streaming computation and saving {log_prefix}Zarr...[/bold cyan]", spinner="dots", ): ds.to_zarr(out_path, encoding=encoding, mode="w", compute=True) - console.log(f"[bold green]✓ Successfully saved Zarr:[/bold green] {out_path}") - dims_str = ", ".join(ds[name].dims) - console.log(f"Output shape was {ds[name].shape} ({dims_str})") + logger.info(f"[bold green]✓ Successfully saved Zarr:[/bold green] {out_path}") + + # Safely get a sample name to log the shape (fixes the scope-leak bug) + sample_name = "total_gmslr" if "total_gmslr" in ds else list(ds.data_vars.keys())[0] + dims_str = ", ".join(ds[sample_name].dims) + logger.info( + f"{log_prefix}Output shape for '{sample_name}' was {ds[sample_name].shape} ({dims_str})" + ) diff --git a/pyproject.toml b/pyproject.toml index 5bc38a1..894724f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,7 +5,17 @@ build-backend = "setuptools.build_meta" [project] name = "profsea" version = "3.0" -dependencies = ["numpy", "xarray", "scipy", "dask", "rich", "requests", "regionmask"] +dependencies = [ + "numpy", + "xarray", + "scipy", + "dask", + "rich", + "requests", + "regionmask", + "netcdf4", + "h5netcdf", +] [tool.setuptools.packages.find] include = ["profsea*"] @@ -14,8 +24,8 @@ exclude = ["docs*", "tests*"] [project.optional-dependencies] test = [ "pytest>=7.0", - "pytest-cov", # Useful for checking test coverage later - "pytest-github-actions-annotate-failures" + "pytest-cov", # Useful for checking test coverage later + "pytest-github-actions-annotate-failures", ] # Settings for Ruff formatting diff --git a/tests/core/test_global_model.py b/tests/core/test_global_model.py index 6309c09..cf3744b 100644 --- a/tests/core/test_global_model.py +++ b/tests/core/test_global_model.py @@ -9,7 +9,7 @@ class MockGlobalComponent(Component): def project(self, state: ClimateState, rng: np.random.Generator) -> np.ndarray: # Just return an array of 1s with the correct shape - return np.ones((state.nt * state.num_members, state.nyr)) + return np.ones((state.nt * state.num_members, state.nyr), dtype=state.dtype) def test_calculate_drivers_math(): @@ -53,7 +53,11 @@ def test_save_components(tmp_path): # Save the output to the temporary directory global_model.save_components( - components, output_dir=str(tmp_path), scenario_name="test" + components, + output_dir=str(tmp_path), + scenario_name="test", + output_prefix="global", + output_format="netcdf", ) expected_file = tmp_path / "test_global.nc" @@ -64,3 +68,53 @@ def test_save_components(tmp_path): assert "mock_comp" in ds.data_vars assert ds["mock_comp"].shape == (5, 4) assert list(ds.dims) == ["member", "time"] + + +def test_return_types(): + # Import a number of real components to test their return types + from profsea.components.global_ import ( + AntarcticaDynAR5, + AntarcticaISMIP6, + AntarcticaSMBAR5, + Glacier, + GreenlandAR6, + GreenlandDynAR5, + GreenlandSMBAR5, + LandwaterAR5, + LandwaterAR6, + ThermalExpansion, + ) + + tas = np.zeros((2, 95)) # 2 trajectories, 4 years + ohc = np.zeros((2, 95)) + + components = { + "AntarcticaISMIP6": AntarcticaISMIP6(region="wais"), + "AntarcticaSMBAR5": AntarcticaSMBAR5(), + "AntarcticaDynAR5": AntarcticaDynAR5(), + "GreenlandAR6": GreenlandAR6(), + "GreenlandDynAR5": GreenlandDynAR5(), + "GreenlandSMBAR5": GreenlandSMBAR5(), + "Glacier": Glacier(), + "LandwaterAR5": LandwaterAR5(), + "LandwaterAR6": LandwaterAR6(), + "ThermalExpansion": ThermalExpansion(OHC_change=ohc), + } + + global_model_float64 = Global( + components=components, end_yr=2101, nt=2, num_members=3, dtype=np.float64 + ) + results = global_model_float64.run(scenario="rcp26", T_change=tas) + + for comp_name, data in results.items(): + assert isinstance(data, xr.DataArray) + assert data.dtype == np.float64 + + global_model_float32 = Global( + components=components, end_yr=2101, nt=2, num_members=3, dtype=np.float32 + ) + results = global_model_float32.run(scenario="rcp26", T_change=tas) + + for comp_name, data in results.items(): + assert isinstance(data, xr.DataArray) + assert data.dtype == np.float32 diff --git a/ukcp18_example_script.py b/ukcp18_example_script.py index a0d05f3..de29cda 100644 --- a/ukcp18_example_script.py +++ b/ukcp18_example_script.py @@ -37,6 +37,12 @@ ) global_model.sum_components(projections) gmslr = global_model.results["total_gmslr"] +global_model.save_components( + global_model.results, + scenario_name="rcp85", + output_format="zarr", + output_prefix="global_", +) ### Now spatial projections ### spatial_components = {