{ "cells": [ { "cell_type": "markdown", "id": "e6597bb3-15b6-4639-82ed-841386591567", "metadata": {}, "source": [ "# Quickstart\n", "\n", "**Acknowledgments:** This notebook adapts the content in [this NCAR tutorial](https://github.com/NCAR/GPU_workshop/blob/workshop/12_CuPyAndLegate/12_CuPyAndLegate.ipynb) to Xarray, and uses it to illustrate `cupy-xarray` and working with cupy arrays and Xarray objects in general." ] }, { "cell_type": "markdown", "id": "71235ea5-c0fb-4afb-a30e-72275ad0b7f0", "metadata": {}, "source": [ "## Setup" ] }, { "cell_type": "code", "execution_count": 1, "id": "860567c3-18b2-444c-93a9-96ce620ea08f", "metadata": {}, "outputs": [], "source": [ "import cupy as cp\n", "import cupy_xarray # Adds .cupy to Xarray objects\n", "import numpy as np\n", "import xarray as xr" ] }, { "cell_type": "markdown", "id": "be9ca22a-8e41-46b1-a43d-6cf93e4fd677", "metadata": {}, "source": [ "## Creating Arrays\n", "\n", "First we create arrays on the CPU and GPU" ] }, { "cell_type": "code", "execution_count": 2, "id": "8aba4551-3b83-4f2b-9fc6-541b0ae071e6", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "On the CPU: [0. 0.5 1. 1.5 2. ]\n", "On the GPU: [2. 2.5 3. 3.5 4. ]\n" ] } ], "source": [ "# NumPy data (host / cpu)\n", "x_cpu = np.linspace(0, 2, 5)\n", "print(\"On the CPU: \", x_cpu)\n", "\n", "# CuPy data\n", "x_gpu = cp.linspace(2, 4, 5)\n", "print(\"On the GPU: \", x_gpu)" ] }, { "cell_type": "markdown", "id": "3976d9b8-d028-4d87-b5c7-70943d807e85", "metadata": {}, "source": [ "And now wrap those in a Xarray DataArray" ] }, { "cell_type": "code", "execution_count": 3, "id": "0bacb2e0-3367-4efe-98af-2176937f84ab", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (x: 5)>\n",
       "array([2. , 2.5, 3. , 3.5, 4. ])\n",
       "Dimensions without coordinates: x
" ], "text/plain": [ "\n", "array([2. , 2.5, 3. , 3.5, 4. ])\n", "Dimensions without coordinates: x" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "da_gpu = xr.DataArray(x_gpu, dims=\"x\")\n", "da_gpu" ] }, { "cell_type": "markdown", "id": "0e536b2b-e540-4a68-b7c8-627fd015832c", "metadata": {}, "source": [ "That was easy! Xarray seamlessly wraps numpy array-like objects that [support specific protocols](https://docs.xarray.dev/en/stable/internals/duck-arrays-integration.html)." ] }, { "cell_type": "markdown", "id": "6fb23b66-4cd9-4e1f-82ca-77b6235dc1d0", "metadata": {}, "source": [ "For array-specific functionality Xarray recommends adding new packages that provide [\"accessors\"](https://docs.xarray.dev/en/stable/internals/extending-xarray.html) on Xarray objects. \n", "\n", "For example, the [pint-xarray](https://pint-xarray.readthedocs.io/en/latest/) package that wraps unit-aware pint arrays and provides a `.pint` for unit-specific functionality.\n", "\n", "In this tutorial, we demonstrate `cupy-xarray` which provides a `cupy` accessor that in turn provides access to cupy-specific functionality." ] }, { "cell_type": "markdown", "id": "8fed5c08-a729-4d55-a167-4fc557c67d74", "metadata": {}, "source": [ "## Checking for cupy arrays\n", "\n", "Unfortunately the text representation of CuPy arrays isn't [very informative](https://github.com/cupy/cupy/issues/6926) so it isn't obvious that this DataArray wraps a CuPy array on the GPU." ] }, { "cell_type": "code", "execution_count": 4, "id": "c582c151-9399-4dd0-a067-2125589bb605", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (x: 5)>\n",
       "array([2. , 2.5, 3. , 3.5, 4. ])\n",
       "Dimensions without coordinates: x
" ], "text/plain": [ "\n", "array([2. , 2.5, 3. , 3.5, 4. ])\n", "Dimensions without coordinates: x" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "da_gpu" ] }, { "cell_type": "markdown", "id": "f089cf99-9d60-49da-a74b-bee9be52fd27", "metadata": {}, "source": [ "Instead we'll use the `is_cupy` property provided by the `cupy` accessor" ] }, { "cell_type": "code", "execution_count": 5, "id": "951776ff-dad4-45ba-8adf-a2a83c68de85", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "da_gpu.cupy.is_cupy" ] }, { "cell_type": "markdown", "id": "78044088-9853-4fbc-ae23-05b8853254d8", "metadata": {}, "source": [ "## Accessing the underlying array\n", "\n", "Use the `DataArray.data` property to access the underlying CuPy Array" ] }, { "cell_type": "code", "execution_count": 6, "id": "8b419a4c-59de-45bd-84f7-b4c4e9268dd2", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([2. , 2.5, 3. , 3.5, 4. ])" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "da_gpu.data" ] }, { "cell_type": "markdown", "id": "0beda1aa-2989-49a1-948d-ee42f96cffb3", "metadata": {}, "source": [ "This means we now have access to CuPy-specific properties" ] }, { "cell_type": "code", "execution_count": 7, "id": "63522d11-fe3b-44dd-b51c-6ba4ba0c98cb", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "da_gpu.data.device" ] }, { "cell_type": "markdown", "id": "10fa58ed-0911-4dea-80fb-d4061b95e2d3", "metadata": {}, "source": [ "## Moving data between CPU and GPU (or host and device)\n", "\n", "Xarray provides [DataArray.as_numpy](https://docs.xarray.dev/en/stable/generated/xarray.Dataset.as_numpy.html#xarray.Dataset.as_numpy) to convert all kinds of arrays to numpy arrays" ] }, { "cell_type": "code", "execution_count": 8, "id": "11b7f877-2087-4d6e-a134-e0fcd084abd7", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (x: 5)>\n",
       "array([2. , 2.5, 3. , 3.5, 4. ])\n",
       "Dimensions without coordinates: x
" ], "text/plain": [ "\n", "array([2. , 2.5, 3. , 3.5, 4. ])\n", "Dimensions without coordinates: x" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Move data to host\n", "da_cpu = da_gpu.as_numpy()\n", "da_cpu" ] }, { "cell_type": "markdown", "id": "86ac2adf-da1b-4bf6-82fb-3472292f2c12", "metadata": {}, "source": [ "Let's make sure this isn't a cupy array anymore" ] }, { "cell_type": "code", "execution_count": 9, "id": "1669e378-1757-4051-9073-d78c5b0b9ff1", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "da_cpu.cupy.is_cupy" ] }, { "cell_type": "markdown", "id": "1c504ad1-7527-4757-8042-bc1641d55506", "metadata": {}, "source": [ "To convert a numpy array to a CuPy array (move data to GPU) use `cupy.as_cupy()`" ] }, { "cell_type": "code", "execution_count": 10, "id": "8b23a051-d936-49e8-b351-aae75e7a88bf", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (x: 5)>\n",
       "array([2. , 2.5, 3. , 3.5, 4. ])\n",
       "Dimensions without coordinates: x
" ], "text/plain": [ "\n", "array([2. , 2.5, 3. , 3.5, 4. ])\n", "Dimensions without coordinates: x" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Move data to GPU\n", "da_cpu.cupy.as_cupy()" ] }, { "cell_type": "code", "execution_count": 11, "id": "3f1e19f6-7084-4e03-90e1-48e418c6d6a0", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "da_cpu.as_cupy().cupy.is_cupy" ] }, { "cell_type": "markdown", "id": "b46e6a52-5de1-4aac-94d0-627ef22f47de", "metadata": {}, "source": [ "## Most Xarray operations preserve array type\n" ] }, { "cell_type": "code", "execution_count": 12, "id": "91b51f09-524f-4797-b7f2-0dafe6185622", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "expanded = da_gpu.expand_dims(y=3)\n", "expanded.cupy.is_cupy" ] }, { "cell_type": "markdown", "id": "f0a735f9-cca7-4cbd-b635-ad57811468a1", "metadata": {}, "source": [ "### Alignment\n", "\n", "Alignment is a fundamental Xarray operation. It preserves array types" ] }, { "cell_type": "code", "execution_count": 13, "id": "0ce4dd47-267a-469a-b539-8aa557254a29", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[True, True]" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "aligned = xr.align(da_gpu, expanded)\n", "[a.cupy.is_cupy for a in aligned]" ] }, { "cell_type": "markdown", "id": "fad69278-82da-40b1-aba2-2b2baa31d859", "metadata": {}, "source": [ "### Broadcasting" ] }, { "cell_type": "code", "execution_count": 14, "id": "9cc6e7a8-9c1d-4429-a962-8827c78cc363", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[True, True]" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "da_gpu2 = da_gpu.rename({\"x\": \"y\"})\n", "broadcasted = xr.broadcast(da_gpu, da_gpu2)\n", "[a.cupy.is_cupy for a in broadcasted]" ] }, { "cell_type": "markdown", "id": "b6829643-c7ea-462d-ac7a-38daa4ab0f24", "metadata": {}, "source": [ "### Basic Arithmetic" ] }, { "cell_type": "code", "execution_count": 15, "id": "a4ad6ac6-c0d1-4e82-a2c9-22d0f1adb2b8", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "is_gpu: True\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (x: 5)>\n",
       "array([3. , 3.5, 4. , 4.5, 5. ])\n",
       "Dimensions without coordinates: x
" ], "text/plain": [ "\n", "array([3. , 3.5, 4. , 4.5, 5. ])\n", "Dimensions without coordinates: x" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# works on both CPU and GPU\n", "print(\"is_gpu: \", (da_gpu + 1).cupy.is_cupy)\n", "da_gpu + 1" ] }, { "cell_type": "markdown", "id": "0cc3b197-5732-489c-829f-43a168ad15a5", "metadata": {}, "source": [ "### Numpy universal functions" ] }, { "cell_type": "code", "execution_count": 16, "id": "9c571c2f-5d61-443e-897c-9b7a098bb18f", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.min(da_gpu.mean()).cupy.is_cupy" ] }, { "cell_type": "markdown", "id": "1c075cc9-e2c3-464f-a504-32d9d3a879d5", "metadata": {}, "source": [ "We can use `np.round` which dispatches" ] }, { "cell_type": "code", "execution_count": 17, "id": "bd0713dc-d86d-4afb-bd10-563a9eb60883", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.round(da_gpu.mean(), 2).cupy.is_cupy" ] }, { "cell_type": "markdown", "id": "19899ada-0871-4cfb-baa4-feb46c9a6981", "metadata": {}, "source": [ "## High-level Xarray functions" ] }, { "cell_type": "markdown", "id": "42394855-a20f-4ec4-b816-3e1cb1a5be54", "metadata": { "tags": [] }, "source": [ "### Groupby works\n", "\n", "Though this is a slow for loop over groups. We could add an explicit parallel algorithm to [flox](https://github.com/xarray-contrib/flox)" ] }, { "cell_type": "code", "execution_count": 18, "id": "fe994a8a-1c45-4c14-add5-6ed38fe95585", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (x: 5)>\n",
       "array([2. , 2.5, 3. , 3.5, 4. ])\n",
       "Dimensions without coordinates: x
" ], "text/plain": [ "\n", "array([2. , 2.5, 3. , 3.5, 4. ])\n", "Dimensions without coordinates: x" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "da_gpu.groupby(\"x\").mean(...)" ] }, { "cell_type": "markdown", "id": "141771cc-0594-411f-9fe2-8a18f289ef32", "metadata": {}, "source": [ "Since groupby works; groupby_bins and resample *should* also work" ] }, { "cell_type": "markdown", "id": "7c7286e4-b651-4778-bbb8-c3f8cc24c78a", "metadata": { "tags": [] }, "source": [ "### Rolling windows do not work\n", "\n", "cupy needs to add support for [sliding_window_view](https://numpy.org/devdocs/reference/generated/numpy.lib.stride_tricks.sliding_window_view.html)" ] }, { "cell_type": "code", "execution_count": 19, "id": "fd5da1e2-26ad-40b3-9fac-2559d33c345a", "metadata": { "tags": [ "raises-exception" ] }, "outputs": [ { "ename": "TypeError", "evalue": "no implementation found for 'numpy.lib.stride_tricks.sliding_window_view' on types that implement __array_function__: []", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", "Input \u001b[0;32mIn [19]\u001b[0m, in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[43mda_gpu\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrolling\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m3\u001b[39;49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mmean\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241m.\u001b[39mcupy\u001b[38;5;241m.\u001b[39mis_cupy\n", "File \u001b[0;32m~/miniconda3/envs/gpu/lib/python3.10/site-packages/xarray/core/rolling.py:155\u001b[0m, in \u001b[0;36mRolling._reduce_method..method\u001b[0;34m(self, keep_attrs, **kwargs)\u001b[0m\n\u001b[1;32m 151\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mmethod\u001b[39m(\u001b[38;5;28mself\u001b[39m, keep_attrs\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs):\n\u001b[1;32m 153\u001b[0m keep_attrs \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_get_keep_attrs(keep_attrs)\n\u001b[0;32m--> 155\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_numpy_or_bottleneck_reduce\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m 156\u001b[0m \u001b[43m \u001b[49m\u001b[43marray_agg_func\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 157\u001b[0m \u001b[43m \u001b[49m\u001b[43mbottleneck_move_func\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 158\u001b[0m \u001b[43m \u001b[49m\u001b[43mrolling_agg_func\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 159\u001b[0m \u001b[43m \u001b[49m\u001b[43mkeep_attrs\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mkeep_attrs\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 160\u001b[0m \u001b[43m \u001b[49m\u001b[43mfillna\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mfillna\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 161\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 162\u001b[0m \u001b[43m \u001b[49m\u001b[43m)\u001b[49m\n", "File \u001b[0;32m~/miniconda3/envs/gpu/lib/python3.10/site-packages/xarray/core/rolling.py:580\u001b[0m, in \u001b[0;36mDataArrayRolling._numpy_or_bottleneck_reduce\u001b[0;34m(self, array_agg_func, bottleneck_move_func, rolling_agg_func, keep_attrs, fillna, **kwargs)\u001b[0m\n\u001b[1;32m 576\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_bottleneck_reduce(\n\u001b[1;32m 577\u001b[0m bottleneck_move_func, keep_attrs\u001b[38;5;241m=\u001b[39mkeep_attrs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs\n\u001b[1;32m 578\u001b[0m )\n\u001b[1;32m 579\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m rolling_agg_func:\n\u001b[0;32m--> 580\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mrolling_agg_func\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mkeep_attrs\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_get_keep_attrs\u001b[49m\u001b[43m(\u001b[49m\u001b[43mkeep_attrs\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 581\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m fillna \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m 582\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m fillna \u001b[38;5;129;01mis\u001b[39;00m dtypes\u001b[38;5;241m.\u001b[39mINF:\n", "File \u001b[0;32m~/miniconda3/envs/gpu/lib/python3.10/site-packages/xarray/core/rolling.py:169\u001b[0m, in \u001b[0;36mRolling._mean\u001b[0;34m(self, keep_attrs, **kwargs)\u001b[0m\n\u001b[1;32m 168\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_mean\u001b[39m(\u001b[38;5;28mself\u001b[39m, keep_attrs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs):\n\u001b[0;32m--> 169\u001b[0m result \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msum\u001b[49m\u001b[43m(\u001b[49m\u001b[43mkeep_attrs\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mFalse\u001b[39;49;00m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m \u001b[38;5;241m/\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcount(keep_attrs\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mFalse\u001b[39;00m)\n\u001b[1;32m 170\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m keep_attrs:\n\u001b[1;32m 171\u001b[0m result\u001b[38;5;241m.\u001b[39mattrs \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mobj\u001b[38;5;241m.\u001b[39mattrs\n", "File \u001b[0;32m~/miniconda3/envs/gpu/lib/python3.10/site-packages/xarray/core/rolling.py:155\u001b[0m, in \u001b[0;36mRolling._reduce_method..method\u001b[0;34m(self, keep_attrs, **kwargs)\u001b[0m\n\u001b[1;32m 151\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mmethod\u001b[39m(\u001b[38;5;28mself\u001b[39m, keep_attrs\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs):\n\u001b[1;32m 153\u001b[0m keep_attrs \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_get_keep_attrs(keep_attrs)\n\u001b[0;32m--> 155\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_numpy_or_bottleneck_reduce\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m 156\u001b[0m \u001b[43m \u001b[49m\u001b[43marray_agg_func\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 157\u001b[0m \u001b[43m \u001b[49m\u001b[43mbottleneck_move_func\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 158\u001b[0m \u001b[43m \u001b[49m\u001b[43mrolling_agg_func\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 159\u001b[0m \u001b[43m \u001b[49m\u001b[43mkeep_attrs\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mkeep_attrs\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 160\u001b[0m \u001b[43m \u001b[49m\u001b[43mfillna\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mfillna\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 161\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 162\u001b[0m \u001b[43m \u001b[49m\u001b[43m)\u001b[49m\n", "File \u001b[0;32m~/miniconda3/envs/gpu/lib/python3.10/site-packages/xarray/core/rolling.py:589\u001b[0m, in \u001b[0;36mDataArrayRolling._numpy_or_bottleneck_reduce\u001b[0;34m(self, array_agg_func, bottleneck_move_func, rolling_agg_func, keep_attrs, fillna, **kwargs)\u001b[0m\n\u001b[1;32m 586\u001b[0m kwargs\u001b[38;5;241m.\u001b[39msetdefault(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mskipna\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;28;01mFalse\u001b[39;00m)\n\u001b[1;32m 587\u001b[0m kwargs\u001b[38;5;241m.\u001b[39msetdefault(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mfillna\u001b[39m\u001b[38;5;124m\"\u001b[39m, fillna)\n\u001b[0;32m--> 589\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mreduce\u001b[49m\u001b[43m(\u001b[49m\u001b[43marray_agg_func\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mkeep_attrs\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mkeep_attrs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n", "File \u001b[0;32m~/miniconda3/envs/gpu/lib/python3.10/site-packages/xarray/core/rolling.py:472\u001b[0m, in \u001b[0;36mDataArrayRolling.reduce\u001b[0;34m(self, func, keep_attrs, **kwargs)\u001b[0m\n\u001b[1;32m 470\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m 471\u001b[0m obj \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mobj\n\u001b[0;32m--> 472\u001b[0m windows \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_construct\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m 473\u001b[0m \u001b[43m \u001b[49m\u001b[43mobj\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mrolling_dim\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mkeep_attrs\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mkeep_attrs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mfill_value\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mfillna\u001b[49m\n\u001b[1;32m 474\u001b[0m \u001b[43m\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 476\u001b[0m result \u001b[38;5;241m=\u001b[39m windows\u001b[38;5;241m.\u001b[39mreduce(\n\u001b[1;32m 477\u001b[0m func, dim\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mlist\u001b[39m(rolling_dim\u001b[38;5;241m.\u001b[39mvalues()), keep_attrs\u001b[38;5;241m=\u001b[39mkeep_attrs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs\n\u001b[1;32m 478\u001b[0m )\n\u001b[1;32m 480\u001b[0m \u001b[38;5;66;03m# Find valid windows based on count.\u001b[39;00m\n", "File \u001b[0;32m~/miniconda3/envs/gpu/lib/python3.10/site-packages/xarray/core/rolling.py:389\u001b[0m, in \u001b[0;36mDataArrayRolling._construct\u001b[0;34m(self, obj, window_dim, stride, fill_value, keep_attrs, **window_dim_kwargs)\u001b[0m\n\u001b[1;32m 384\u001b[0m window_dims \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_mapping_to_list(\n\u001b[1;32m 385\u001b[0m window_dim, allow_default\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mFalse\u001b[39;00m, allow_allsame\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mFalse\u001b[39;00m \u001b[38;5;66;03m# type: ignore[arg-type] # https://github.com/python/mypy/issues/12506\u001b[39;00m\n\u001b[1;32m 386\u001b[0m )\n\u001b[1;32m 387\u001b[0m strides \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_mapping_to_list(stride, default\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m1\u001b[39m)\n\u001b[0;32m--> 389\u001b[0m window \u001b[38;5;241m=\u001b[39m \u001b[43mobj\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mvariable\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrolling_window\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m 390\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdim\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mwindow\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mwindow_dims\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcenter\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mfill_value\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mfill_value\u001b[49m\n\u001b[1;32m 391\u001b[0m \u001b[43m\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 393\u001b[0m attrs \u001b[38;5;241m=\u001b[39m obj\u001b[38;5;241m.\u001b[39mattrs \u001b[38;5;28;01mif\u001b[39;00m keep_attrs \u001b[38;5;28;01melse\u001b[39;00m {}\n\u001b[1;32m 395\u001b[0m result \u001b[38;5;241m=\u001b[39m DataArray(\n\u001b[1;32m 396\u001b[0m window,\n\u001b[1;32m 397\u001b[0m dims\u001b[38;5;241m=\u001b[39mobj\u001b[38;5;241m.\u001b[39mdims \u001b[38;5;241m+\u001b[39m \u001b[38;5;28mtuple\u001b[39m(window_dims),\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 400\u001b[0m name\u001b[38;5;241m=\u001b[39mobj\u001b[38;5;241m.\u001b[39mname,\n\u001b[1;32m 401\u001b[0m )\n", "File \u001b[0;32m~/miniconda3/envs/gpu/lib/python3.10/site-packages/xarray/core/variable.py:2319\u001b[0m, in \u001b[0;36mVariable.rolling_window\u001b[0;34m(self, dim, window, window_dim, center, fill_value)\u001b[0m\n\u001b[1;32m 2315\u001b[0m axis \u001b[38;5;241m=\u001b[39m [\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mget_axis_num(d) \u001b[38;5;28;01mfor\u001b[39;00m d \u001b[38;5;129;01min\u001b[39;00m dim]\n\u001b[1;32m 2316\u001b[0m new_dims \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdims \u001b[38;5;241m+\u001b[39m \u001b[38;5;28mtuple\u001b[39m(window_dim)\n\u001b[1;32m 2317\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m Variable(\n\u001b[1;32m 2318\u001b[0m new_dims,\n\u001b[0;32m-> 2319\u001b[0m \u001b[43mduck_array_ops\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msliding_window_view\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m 2320\u001b[0m \u001b[43m \u001b[49m\u001b[43mpadded\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdata\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mwindow_shape\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mwindow\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43maxis\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43maxis\u001b[49m\n\u001b[1;32m 2321\u001b[0m \u001b[43m \u001b[49m\u001b[43m)\u001b[49m,\n\u001b[1;32m 2322\u001b[0m )\n", "File \u001b[0;32m~/miniconda3/envs/gpu/lib/python3.10/site-packages/xarray/core/duck_array_ops.py:640\u001b[0m, in \u001b[0;36msliding_window_view\u001b[0;34m(array, window_shape, axis)\u001b[0m\n\u001b[1;32m 638\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m dask_array_compat\u001b[38;5;241m.\u001b[39msliding_window_view(array, window_shape, axis)\n\u001b[1;32m 639\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m--> 640\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mnpcompat\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msliding_window_view\u001b[49m\u001b[43m(\u001b[49m\u001b[43marray\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mwindow_shape\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43maxis\u001b[49m\u001b[43m)\u001b[49m\n", "File \u001b[0;32m<__array_function__ internals>:180\u001b[0m, in \u001b[0;36msliding_window_view\u001b[0;34m(*args, **kwargs)\u001b[0m\n", "\u001b[0;31mTypeError\u001b[0m: no implementation found for 'numpy.lib.stride_tricks.sliding_window_view' on types that implement __array_function__: []" ] } ], "source": [ "da_gpu.rolling(x=3).mean().cupy.is_cupy" ] }, { "cell_type": "markdown", "id": "5296f2b6-d8c7-49c1-8d39-c7231f7ad70b", "metadata": { "tags": [] }, "source": [ "### Weighted operations work" ] }, { "cell_type": "code", "execution_count": 20, "id": "7a8ab1a5-48c0-448f-a062-dd27b3956b3f", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "weights = xr.DataArray(cp.asarray([0, 0.5, 1, 0.5, 0]), dims=\"y\")\n", "da_gpu.weighted(weights).sum().cupy.is_cupy" ] }, { "cell_type": "markdown", "id": "5aaa21ee-781a-4f76-9a8c-61291f3ecf57", "metadata": {}, "source": [ "## Plotting works\n", "\n", "Automatically moves data to the CPU before passing on to matplotlib" ] }, { "cell_type": "code", "execution_count": 21, "id": "153b4e92-55a8-4a0a-9ed3-eb825483936c", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvQAAAILCAYAAACHGAjaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAABYlAAAWJQFJUiTwAABe5UlEQVR4nO3dd3hU55n38e8jISE6pttgm96RHfeSuMe9UbybfdN34zTvOs7GMbgQ47gB6T2b5myS3ThZcO+9O3HiJjqmY0wxYDpC7Xn/mGGQCAKVEUcjfT/Xpeuge8555p7jY/jp0SkhxogkSZKk3JSXdAOSJEmSGs5AL0mSJOUwA70kSZKUwwz0kiRJUg4z0EuSJEk5zEAvSZIk5TADvSRJkpTDDPSSJElSDjPQS5IkSTnMQC9JkiTlMAO9JEmSlMMM9JIkSVIOa5N0A81dCGEp0BlYlnArkiRJatn6A1tijAPqs5GB/sA6t2vXrtuIESO6Jd2IJEmSWq558+axc+fOem9noD+wZSNGjOj2+uuvJ92HJEmSWrBjjz2WN954Y1l9t/McekmSJCmHGeglSZKkHGaglyRJknKYgV6SJEnKYQZ6SZIkKYcZ6CVJkqQcZqCXJEmScpiBXpIkScphBnpJkiQphxnoJUmSpBxmoJckSZJymIFekiRJymFNEuhDCJ8MIcT01+fquW2/EMJvQgjvhRB2hRCWhRC+H0I4ZD/bnBJCeCSEsDGEsCOEUBJCuCaEkN/4TyNJkiQ1X1kP9CGEw4EfAdsasO0g4HXgs8BrwPeAJcBXgFdDCN33sc1lwAvAacC9wE+AwvS2dzfsU0iSJEm5IauBPoQQgLuADcDPGzDET4FewNUxxstjjJNijGeRCufDgNv3er/OwC+BSuCMGOO/xRi/DhwNvApMCCF8rKGfR5IkSWrusj1DfzVwFqkZ9u312TCEMBA4F1hGapa9upvT430yhNChWn0C0BO4O8b4993FGGMpcFP62y/Vpw9JkiS1XltLy3lk1uqk26iXNtkaKIQwApgK/CDG+EII4ax6DrF7/SdijFXVX4gxbg0hvEwq8J8EPL3XNo/tY7wXgB3AKSGEtjHGXQfo//VaXhpel+YlSZKU256dv44b7p3Fmi2lzPzSKRxzRK2XcDYrWZmhDyG0AX4PrABuaOAww9LLhbW8/k56ObQu28QYK4ClpH5oGdjAniRJktTCbdxexjV3v8lnf/s3Vm8uJUaYOKOEsoqqA2/cDGRrhv4bwIeAD8cYdzZwjC7p5eZaXt9d79rIbfYpxnjsvurpmftjDrS9JEmSckuMkYdKVjPlgTls2F6WqXfvUMjVZw+hID8k2F3dNTrQhxBOIDUr/50Y46uNb6n2t0ovYxNvI0mSpBZu7ZZSbrx3Nk/NW1ujfvnRh/GNS0bRrUNhQp3VX6MCfbVTbRYCkxvZy+7Z9C61vN55r/Uauo0kSZJaqRgjf/rbSm5/ZB5bSysy9UO7FHH72NGcNbx3gt01TGNn6Duy55z20tRdK//BL0MIvyR1sew1+xlrQXo5tJbXh6SX1c+XXwAcl96mxkWt6R82BgAVpO5lL0mSpFZsxYYdTLqnhFcWb6hR/8RJRzDx/OF0KipIqLPGaWyg3wX8upbXjiF1Xv1LpIL3gU7HeTa9PDeEkFf9TjchhE7AqcBO4C/VtnkG+DhwPvDHvcY7DWgPvHCgO9xIkiSp5aqsitz18lK+/cQCSsv3XOjav3t7po4v5qSB//Ds0pzSqECfvgD2c/t6LYQwhVSg/+8Y46+q1QuAQUB5jHFxtbEWhxCeIHVryqtIPW12t1uADsB/xRir399+BjAN+FgI4Ue770UfQigCbkuv87PGfEZJkiTlrgVrtnLdzBLeXrkpU8sLcOVpA/nqOUMpKshPrrksydp96OuhLzAPWA703+u1LwOvAD8MIZydXu9E4ExSp9rcWH3lGOOWEMKVpIL9cyGEu4GNwKWkbmk5A/hTk30SSZIkNUtlFVX89LlF/OTZRZRX7rk/yvA+nZg+oZjifl2Tay7Lkgj0tUrP0h8HfJPUaTQXAquBHwK3xBg37mOb+0IIp5MK++OBImAR8J/AD2OM3uFGkiSpFXl75Saum1HCgrVbM7XC/Dz+46zBfOH0QRS2ycqjmJqNJgv0McYpwJR91Jex53aS+9puJfDZer7Xy6TCvyRJklqpnWWVfPfJBfz6paVUVZvS/dARXZk+vpghvTsl11wTalYz9JIkSVJDvLJ4PdffM4vlG3Zkau0K8vn6ecP49Cn9yc/LjYdENYSBXpIkSTlrS2k5dz4ynz++tqJG/cODe3DnuDEc3q19Qp0dPAZ6SZIk5aSn5q7lxvtmsXbLnjuUdypqw+SLRnLFcf2o5RlJLY6BXpIkSTllw7Zd3PLgXB54+70a9XNH9ubWy0fTu3NRQp0lw0AvSZKknBBj5IG332PKA3P4YEd5pt6jYyHfvGw0F4zu02pm5asz0EuSJKnZW715JzfdO5un56+rUR93TF8mXzSSQzoUJtRZ8gz0kiRJaraqqiJ//NsK7nxkPtt2VWTqfbu24/axozljWK8Eu2seDPSSJElqlpau386kmSX8dWnNZ4t+6uQjue784XRsa5QFA70kSZKamYrKKn7z8lK+88RCdlVUZeoDe3Rg6vhiThjQLcHumh8DvSRJkpqNeau3MHFmCSXvbs7U8vMCnz9tIF85ewhFBfkJdtc8GeglSZKUuF0Vlfzk2cX89NlFVFTFTH3koZ2ZPqGY0X27JNhd82aglyRJUqLeWPEBE2eU8M66bZlaYZs8vnL2ED5/2kAK8vMS7K75M9BLkiQpETvKKvjOEwv5zctLiXsm5Tn2yEOYNr6Ywb06JtdcDjHQS5Ik6aB7edF6Jt1TwsqNOzO19oX5TDx/OJ886Ujy8lrfA6IaykAvSZKkg2bzznLueHgef/r7yhr1jwzpwR1jx3B4t/YJdZa7DPSSJEk6KJ6Ys4ab7pvNuq27MrUu7QqYfPFIxh/TlxCclW8IA70kSZKa1PtbdzHlwTk8XLK6Rv3CMX2YcukoenUqSqizlsFAL0mSpCYRY+S+t1Zxy4Nz2bSjPFPv0bEtt10+ivNHH5pgdy2HgV6SJElZt2rTTm68dxbPLXi/Rv2KY/tx00Uj6dK+IKHOWh4DvSRJkrKmqiryP39dztRH57O9rDJT79u1HXeOG8NpQ3sm2F3LZKCXJElSVix5fxuTZs7itWUbM7UQ4NMn9+fr5w2jQ1ujZ1Nwr0qSJKlRKiqr+OWLS/neUwspq6jK1Af17MC08cUc179bgt21fAZ6SZIkNdjc97Zw3cy3mb1qS6bWJi/wxdMH8e9nDaaoID/B7loHA70kSZLqrbS8kh8/s4ifP7+YiqqYqY/u25lp44sZdViXBLtrXQz0kiRJqpfXl2/kuhklLH5/e6ZW2CaPr54zlCs/MoA2+XkJdtf6GOglSZJUJ9t3VfCtxxfw368uI+6ZlOeE/t2YOn4MA3t2TK65VsxAL0mSpAN6YeH7XH/PLFZt2pmpdSjMZ9IFw/n4iUeSlxcS7K51M9BLkiSpVpt3lHPrw3OZ8fq7NepnDOvJ7WPH0Ldru4Q6024GekmSJO3TY7NXM/n+Oby/dVem1rV9ATdfMpLLj+5LCM7KNwcGekmSJNWwbmspN98/h0dnr6lRv6j4UG65dBQ9OrZNqDPti4FekiRJAMQYmfnGKm59aC6bd5Zn6r06teXWy0dz3qg+CXan2hjoJUmSxMqNO7jh3lm8+M76GvV/Pu5wbrhoBF3aFSTUmQ7EQC9JktSKVVVFfvfqMqY/voAdZZWZ+uHd2jF1XDGnDu6RYHeqCwO9JElSK7Vo3TYmzSzh78s/yNRCgH89dQBfO3co7QuNirkgK/+VQgjTgOOAoUAPYCewHLgP+HGMcUMdxvgMcNcBVquKMeZX26Y/sHQ/6/8pxvixA723JElSa1JeWcUvXljCD556h7LKqkx9SK+OTJtQzDFHHJJgd6qvbP3Y9VXgDeBJYB3QATgJmAJ8PoRwUoxx5QHGeAu4pZbXPgKcBTxay+tvk/rhYW+zD/CekiRJrcrsVZu5bkYJc1dvydTa5AW+fOZgrjpzEG3b5O9nazVH2Qr0nWOMpXsXQwi3AzcA1wNf3t8AMca3SIX6fxBCeDX9x1/UsvlbMcYpdexVkiSp1Sktr+QHT7/DL15YQmVVzNSL+3Vh2vhiRhzaOcHu1BhZCfT7CvNpfyYV6Ic0dOwQwmhSs/2rgIcbOo4kSVJr9bdlG5k4o4Ql67dnam3b5PG1c4fyr6cOoE1+XoLdqbGa+kqHS9LLkkaM8YX08tcxxspa1jkshPAFoDuwAXg1xtiY95QkScp523ZVMP2x+fzu1eU16icO6Ma08cX079Ehoc6UTVkN9CGEa4GOQBdSF8l+mFSYn9rA8doBnwCqgF/tZ9WPpr+qb/sc8OkY44o6vtfrtbw0vC7bS5IkNSfPLVjHjffOZtWmnZlax7ZtuP7C4fzL8UeQlxcS7E7ZlO0Z+muB3tW+fwz4TIzx/QaO909AV+DhWi6q3QHcSuqC2CXpWjGpi3HPBJ4OIRwdY9y+j20lSZJanA+2l3Hrw3O5541VNepnDe/F7WNHc2iXdgl1pqaS1UAfY+wDEELoDZxCamb+zRDCxTHGNxow5OfTy/+q5f3WAd/Yq/xCCOFc4CXgROBzwA/q0Pux+6qnZ+6PqWvDkiRJSYgx8ujsNXzj/tms31aWqXfrUMjNl4zk0qMOIwRn5VuiJjmHPsa4Frg3hPAGsBD4HTC6PmOEEEaS+qHgXeCRer5/RQjhV6QC/WnUIdBLkiTlqnVbSpl8/2wen7O2Rv3Sow7j5ktG0r1j24Q608HQpBfFxhiXhxDmAkeHEHrEGNfXY/O6XAy7P7tP8/FqD0mS1CLFGPm/v7/LrQ/PZWtpRabep3MRt10+mnNG9t7P1mopDsbzfA9LL+scykMIRcAnSV0M++sGvu9J6eWS/a4lSZKUg1Zu3MH198zipUU150v/5YQjuP7C4XQuKkioMx1sjQ70IYThwKYY45q96nmkLljtBbwSY/wgXS8ABgHlMcbFtQx7BXAI8ND+njAbQjgReDPGWLZX/SxST68F+EP9P5UkSVLzVFkV+e9XlvGtxxews3zPfOmR3dtz57gxnDKoR4LdKQnZmKE/H/hWCOEFYDGp+8D3Bk4HBgJrgCurrd8XmAcsB/rXMubui2FrezLsbtOAUelbVL6brhUDZ6X/PDnG+EpdP4gkSVJz9s7arUycWcIbKzZlankBPveRgXz1nKG0K8xPrjklJhuB/ilSwftU4ChSt5ncTupi2N8DP4wxbqzrYCGEEaTuX1+Xi2F/D4wFjgcuAAqAtaSeUPvjGOOL9fkgkiRJzVFZRRX/9fxifvTMIsoqqzL1Yb07MW1CMUcf3jW55pS4Rgf6GONs4Kp6rL8MqPWeSTHGeft7fa91f03Dz7GXJElq9kre3cR1M0qYv2ZrplaQH/j3M4fwpTMGUdgmL8Hu1BwcjItiJUmSVE+l5ZV878mF/PLFJVTFPfWjDu/K9PHFDOvTKbnm1KwY6CVJkpqZvyzZwKSZJSzbsCNTKyrI49pzh/HZUweQn+cDorSHgV6SJKmZ2FpaztRH5/M/f11Ro37KoO7cOW4MR3b38Tr6RwZ6SZKkZuDZ+eu44d5ZrN5cmql1atuGGy8awT8ffzghOCuvfTPQS5IkJWjj9jK++eAc7nvrvRr1c0b05rbLR9OnS1FCnSlXGOglSZISEGPkwZLVTHlgDhu373lGZvcOhUy5dBQXFx/qrLzqxEAvSZJ0kK3ZXMpN983mqXlra9THfqgvky8eSbcOhQl1plxkoJckSTpIYozc/beV3PHwPLbuqsjUD+1SxO1jR3PW8N4JdqdcZaCXJEk6CJZv2M6kmbN4dcmGGvVPnHQEE88fTqeigoQ6U64z0EuSJDWhyqrIXS8v5dtPLKC0vCpTH9CjA1PHjeHEgd0T7E4tgYFekiSpiSxYs5XrZpbw9spNmVpegCtPG8hXzxlKUUF+cs2pxTDQS5IkZVlZRRU/fW4RP3l2EeWVMVMf3qcT0ycUU9yva3LNqcUx0EuSJGXRWys3MXFGCQvWbs3UCvPz+I+zBvPFMwZRkJ+XYHdqiQz0kiRJWbCzrJLvPrmAX7+0lKo9k/J86IiuTB9fzJDenZJrTi2agV6SJKmRXlm8nkkzZ7Fi445MrV1BPtedP4xPndyf/DwfEKWmY6CXJElqoC2l5dz5yHz++NqKGvUPD+7BnePGcHi39gl1ptbEQC9JktQAT81dy433zWLtll2ZWueiNtx08UiuOLYfITgrr4PDQC9JklQPG7btYsqDc3nw7fdq1M8b1ZtbLxtNr85FCXWm1spAL0mSVAcxRh54+z2mPDCHD3aUZ+o9OhbyzctGc8HoPs7KKxEGekmSpAN4b9NObrpvNs/MX1ejPv6Yftx00QgO6VCYUGeSgV6SJKlWVVWRP/5tBXc+Mp9tuyoy9b5d23HHuDGcPrRngt1JKQZ6SZKkfVi6fjuTZpbw16Uba9Q/ffKRfP384XRsa4xS8+CRKEmSVE1FZRW/fmkp331yIbsqqjL1gT07MG18Mcf375Zgd9I/MtBLkiSlzVu9hYkzSyh5d3Omlp8X+MJpA7n67CEUFeQn2J20bwZ6SZLU6u2qqOQnzyzip88tpqIqZuojD+3M9AnFjO7bJcHupP0z0EuSpFbtjRUfMHFGCe+s25apFbbJ4ytnD+Hzpw2kID8vwe6kAzPQS5KkVmlHWQXffnwhd72ylLhnUp7jjjyEqeOLGdyrY3LNSfVgoJckSa3Oy4vWM+meElZu3JmptS/MZ+L5w/nkSUeSl+cDopQ7DPSSJKnV2LyznDsensef/r6yRv20oT25Y+xo+h3SPqHOpIYz0EuSpFbh8TlrmHzfbNZt3ZWpdWlXwOSLRzL+mL6E4Ky8cpOBXpIktWjvb93FlAfm8PCs1TXqF47pw5RLR9GrU1FCnUnZYaCXJEktUoyRe99cxTcfmsumHeWZes9Obbn1slGcP/rQBLuTssdAL0mSWpxVm3Zy472zeG7B+zXqVxzbj5suGkmX9gUJdSZln4FekiS1GFVVkf/563KmPjqf7WWVmXq/Q9px57gxfGRIzwS7k5qGgV6SJLUIi9/fxvUzZ/Haso2ZWgjw6ZP78/XzhtGhrbFHLVNWjuwQwjTgOGAo0APYCSwH7gN+HGPcUMdxlgFH1vLy2hhjn1q2OwW4CTgJKAIWAb8BfhRjrNzXNpIkqWWoqKziFy8u4ftPvUNZRVWmPqhnB6ZPKObYI7sl2J3U9LL1o+pXgTeAJ4F1QAdS4XoK8PkQwkkxxpW1b17DZuD7+6hv20eNEMJlwEygFPgTsBG4BPgecCpwRV0/hCRJyi1z3tvMxJklzF61JVNrkxf44umD+PezBlNUkJ9gd9LBka1A3znGWLp3MYRwO3ADcD3w5TqOtSnGOKUuK4YQOgO/BCqBM2KMf0/XJwPPABNCCB+LMd5dx/eWJEk5oLS8kh898w4/f34JlVUxUx/dtzPTxhcz6rAuCXYnHVx52RhkX2E+7c/p5ZBsvM8+TAB6AnfvDvPV+rkp/e2Xmui9JUlSAl5fvpGLfvgiP3l2cSbMt22Tx6QLhnPfl081zKvVaeqrQy5JL0vqsU3bEMIngCOA7eltX6jlXPiz0svH9vHaC8AO4JQQQtsY4659rJMRQni9lpeG161tSZLUlLbvquBbjy/gv19dRtwzKc8J/bsxdfwYBvbsmFxzUoKyGuhDCNcCHYEupC6S/TCpQD61HsP0AX6/V21pCOGzMcbn96oPSy8X7j1IjLEihLAUGAUMBObVowdJktSMvLDwfa6/ZxarNu3M1DoU5jPpwhF8/IQjyMsLCXYnJSvbM/TXAr2rff8Y8JkY4/u1rL+3u4AXgTnAVlJB/N+BzwOPhhBOjjG+XW393b9T21zLeLvrXQ/0xjHGY/dVT8/cH3PAziVJUtZt2lHGbQ/PY8br79aonzGsJ7ePHUPfru0S6kxqPrIa6HffVjKE0Bs4hdTM/JshhItjjG/UYftb9irNBr4YQtgGfI3UXXPG1qOl3T+ux/2uJUmSmp1HZ61m8v1zWL9tz1mzXdsXcPMlI7n86L6E4Ky8BE10Dn2McS1wbwjhDVKnw/wOGN2IIX9OKtCftld99wx8bVe/dN5rPUmS1Myt21rKzffP4dHZa2rULy4+lCmXjqJHx7YJdSY1T016UWyMcXkIYS5wdAihR4xxfQOHWpdedtirvoA9D7SqcVFrCKENMACoAJY08H0lSdJBEmNkxuvvctvD89i8szxT79WpLbddPppzR+3z+ZJSq3cwnoF8WHrZmCe2npxe7h3MnwE+DpwP/HGv104D2pO6Q85+73AjSZKStXLjDm64dxYvvlNz7u9jxx/O9ReOoEu7goQ6k5q/Rgf6EMJwUg+DWrNXPQ+4FegFvBJj/CBdLwAGAeUxxsXV1h8FrI4xbtxrnCOBH6e//cNebz8DmAZ8LITwo2oPlioCbkuv87PGfkZJktQ0qqoiv3t1GdMfX8COsj1zf4d3a8fUccWcOrhHgt1JuSEbM/TnA98KIbwALAY2kLrTzemk7lKzBriy2vp9Sd1CcjnQv1r9CmBSCOFZYCmpu9wMAi4CioBHgG9Xf+MY45YQwpWkgv1zIYS7gY3ApaRuaTkD+FMWPqMkScqyReu2MnHmLF5f/kGmFgL866kD+Nq5Q2lfeDBOJJByXzb+T3kK+AVwKnAUqVtEbid1MezvgR/uPetei2dJhfAPkTrFpgOwCXgpPc7vY4z/cLeaGON9IYTTgRuB8aTC/yLgP9Pv7R1uJElqRsorq/jFC0v4wVPvUFZZlakP6dWRaROKOeaIQxLsTso9jQ70McbZwFX1WH8Ze24nWb3+PLD3g6PqOubLwIUN2VaSJB08s1dt5roZJcxdvSVTa5MXuOrMwXz5zEG0bZOfYHdSbvJ3WZIkqcmVllfyg6ff4RcvLKGyas8vz4v7dWH6hGKG9+m8n60l7Y+BXpIkNanXlm5k0swSlqzfnqm1bZPHtecO47On9qdNfl6C3Um5z0AvSZKaxLZdFUx7dD6//8vyGvUTB3Rj2vhi+vfY+/EykhrCQC9JkrLu2QXruPGeWby3uTRT69i2DTdcOIKPHX84eXn/cDmdpAYy0EuSpKz5YHsZtz40l3veXFWjfvbwXtw2djSHdmmXUGdSy2WglyRJjRZj5JFZa7j5gdms31aWqXfrUMjNl4zk0qMOIwRn5aWmYKCXJEmNsnZLKZPvm80Tc9fWqF929GF84+KRdO/YNqHOpNbBQC9Jkhokxsif/76S2x6ex9bSiky9T+cibrt8NOeM7J1gd1LrYaCXJEn1tmLDDq6/t4SXF22oUf9/Jx7BpAuG07moIKHOpNbHQC9Jkuqssiry21eW8e3HF7CzvDJTP7J7e+4cN4ZTBvVIsDupdTLQS5KkOnln7Vaum1nCmys2ZWp5AT73kYF89ZyhtCvMT645qRUz0EuSpP0qq6ji588v5sfPLKKssipTH96nE9PGF3PU4V2Ta06SgV6SJNWu5N1NXDejhPlrtmZqBfmBfz9zCF86YxCFbfIS7E4SGOglSdI+7Cyr5PtPLeSXLy6hKu6pH314V6ZPKGZo707JNSepBgO9JEmq4S9LNjBpZgnLNuzI1IoK8rj23GF89tQB5Of5gCipOTHQS5IkALaWljP10fn8z19X1KifMqg7U8cVc0T39gl1Jml/DPSSJIln5q/lxntns3pzaabWqagNN100gn867nBCcFZeaq4M9JIktWIbtu3imw/N5f633qtR/+jI3tx2+Wh6dy5KqDNJdWWglySpFYox8mDJaqY8MIeN28sy9e4dCrnlslFcNOZQZ+WlHGGglySplVmzuZSb7pvFU/PW1aiP/VBfvnHxSA7pUJhQZ5IawkAvSVIrEWPk7r+t5I6H57F1V0WmfmiXIu4YO4Yzh/dKsDtJDWWglySpFVi+YTuTZs7i1SUbatQ/edKRXHf+MDoVFSTUmaTGMtBLktSCVVZF7np5Kd9+YgGl5VWZ+oAeHZg6bgwnDuyeYHeSssFAL0lSC7VgzVaum1nC2ys3ZWr5eYErPzKQa84ZQlFBfnLNScoaA70kSS1MWUUVP3l2ET99bhHllTFTH3FoZ6aPL2ZMvy4Jdicp2wz0kiS1IG+t3MR1M95m4dptmVphfh5Xnz2YL5w+iIL8vAS7k9QUDPSSJLUAO8sq+c4TC/jNy0up2jMpzzFHdGX6hGIG9+qUXHOSmpSBXpKkHPfK4vVMmjmLFRt3ZGrtC/O57rxhfPLk/uTn+YAoqSUz0EuSlKO2lJZz5yPz+ONrK2vUPzKkB3eMHcPh3don1Jmkg8lAL0lSDnpy7lpuum8Wa7fsytQ6F7Vh8sUjmXBsP0JwVl5qLQz0kiTlkPXbdjHlgTk8VLK6Rv28Ub259bLR9OpclFBnkpJioJckKQfEGLn/rfe45cE5fLCjPFPv0bEtt142igvGHJpgd5KSZKCXJKmZe2/TTm66bzbPzF9Xoz7+mH5MvngEXdsXJtSZpObAQC9JUjNVVRX539dWMPXR+WzbVZGp9+3ajjvGjeH0oT0T7E5Sc5GVQB9CmAYcBwwFegA7geXAfcCPY4wb6jBGd2AscBEwBugLlAGzgLuAu2KMVXtt0x9Yup9h/xRj/Fg9P44kSYlbun47E2eW8NrSjZlaCPCpk47k6+cPp2Nb5+QkpWTrb4OvAm8ATwLrgA7AScAU4PMhhJNijCtr3xyAK4CfAauBZ4EVQG9gHPAr4IIQwhUxxriPbd8m9cPD3mbX+5NIkpSgisoqfvXSUr735EJ2VeyZxxrYswPTxhdzfP9uCXYnqTnKVqDvHGMs3bsYQrgduAG4HvjyAcZYCFwKPFx9Jj6EcAPwGjCeVLifuY9t34oxTmlY65IkNQ9z39vCxJklzFq1OVPLzwt84bSBXH32EIoK8hPsTlJzlZVAv68wn/ZnUoF+SB3GeKaW+poQws+B24Ez2HeglyQpZ+2qqOTHzyziZ88tpqJqzy+iRx7amekTihndt0uC3Ulq7pr6BLxL0suSRo6z+/5cFbW8flgI4QtAd2AD8GqMsbHvKUlSk3t9+QdMnFnConXbMrXCNnlcc84QrvzIQAry8xLsTlIuyGqgDyFcC3QEupC6SPbDpML81EaM2Qb4VPrbx2pZ7aPpr+rbPQd8Osa4oo7v83otLw2vy/aSJNXHjrIKvvX4An77yjKqXx12fP9DmDq+mEE9OybXnKScku0Z+mtJXci622PAZ2KM7zdizKnAaOCRGOPje722A7iV1AWxS9K1YlIX454JPB1CODrGuL0R7y9JUla99M56Jt1Twrsf7MzUOhTmM/GC4XzixCPJywsJdicp12Q10McY+wCEEHoDp5AK42+GEC6OMb5R3/FCCFcDXwPmA5/cx/utA76xV/mFEMK5wEvAicDngB/Uofdja+nhdeCY+nUuSdI/2ryjnNsfmcuf//5ujfppQ3tyx9jR9DukfUKdScplTXIOfYxxLXBvCOENUnev+R2pWfY6CyFcRSqIzwXOjjFuPMAm1d+/IoTwK1KB/jTqEOglSWpKj81ew+T7Z/P+1l2ZWpd2BXzj4pGMO6YvITgrL6lhmvSi2Bjj8hDCXODoEEKPGOP6umwXQrgG+B6p+8ifnZ6Jr6/dp/l0aMC2kiRlxftbdzHlgTk8PGt1jfpFYw5lyqWj6NmpbUKdSWopDsZj5g5LLyvrsnIIYSKpU3XeAj5a1x8C9uGk9HLJfteSJKkJxBi5541VfPOhuWzeWZ6p9+zUllsvG835o/sk2J2klqTRgT6EMBzYFGNcs1c9j9QFq72AV2KMH6TrBcAgoDzGuHivbSYD3wReB8490Gk2IYQTgTdjjGV71c8i9fRagD809LNJktQQqzbt5IZ7ZvH8wpr3hPin4/px44Uj6dK+IKHOJLVE2ZihPx/4VgjhBWAxqfvA9wZOBwYCa4Arq63fF5gHLAf67y6GED5NKsxXAi8CV+/jfMJlMcbfVvt+GjAqfYvK3VcYFQNnpf88Ocb4SqM+nSRJdVRVFfnDX5cz7dH5bC/b84vpfoe0Y+q4Yj48pEeC3UlqqbIR6J8CfgGcChwFdAW2k7oY9vfAD+t4QeuA9DIfuKaWdZ4Hflvt+98DY4HjgQuAAmAtqSfU/jjG+GLdP4YkSQ23+P1tTJpZwt+WfZCphQCfOaU/1547jA5tD8ZZrpJao0b/7RJjnA1cVY/1lwH/MPUeY5xC6v7x9XnvXwO/rs82kiRlU3llFb98cQnff+odyiqqMvXBvToybXwxxx55SILdSWoNnC6QJKmBZq/azMSZJcx5b0um1iYv8OUzBnHVWYNp2yY/we4ktRYGekmS6qm0vJIfPfMOP39+CZVVMVMf07cL08YXM/Kwzgl2J6m1MdBLklQPf1+2ketmlrDk/e2ZWts2efznR4fybx8eQJv8vAS7k9QaGeglSaqDbbsq+NZj8/ndX5YT90zKc8KAbkwdN4aBPTsm15ykVs1AL0nSATy/8H1uuGcWqzbtzNQ6tm3DpAuG8/9OOIK8vH+414MkHTQGekmSarFpRxm3PjSPmW+8W6N+xrCe3DF2DId1bZdQZ5K0h4FekqR9eHTWaibfP4f123Zlaoe0L+DmS0Zx2dGHsY+HH0pSIgz0kiRVs25LKd+4fw6PzVlTo37JUYdx8yUj6dGxbUKdSdK+GeglSQJijMx4/V1ufWguW0orMvXendty2+Vj+OjI3gl2J0m1M9BLklq9lRt3cMO9s3jxnfU16v9ywuFMumAEXdoVJNSZJB2YgV6S1GpVVkV+9+oyvvX4AnaUVWbqR3Rrz9RxYzhlcI8Eu5OkujHQS5JapUXrtjJx5ixeX/5BppYX4F9PHcB/njuU9oX+EykpN/i3lSSpVSmvrOK/nl/MD59eRFllVaY+tHdHpo0v5kNHHJJgd5JUfwZ6SVKrMevdzVw3s4R5q7dkagX5gS+fMZirzhxMYZu8BLuTpIYx0EuSWrzS8kq+/9Q7/PLFJVRWxUz9qH5dmDahmOF9OifYnSQ1joFektSi/XXJBibdM4ul67dnakUFeXzto8P41w8PID/PB0RJym0GeklSi7S1tJzpjy3g939ZXqN+0sBuTB1XTP8eHRLqTJKyy0AvSWpxnl2wjhvvmcV7m0sztU5t23DDRSP45+MOJ89ZeUktiIFektRifLC9jFsfmss9b66qUT9nRC9uu3wMfboUJdSZJDUdA70kKefFGHl41mpuvn8OG7aXZerdOhQy5dJRXFJ8KCE4Ky+pZTLQS5Jy2totpdx032yenLu2Rv2yow/j5ktG0a1DYUKdSdLBYaCXJOWkGCN//vtKbnt4HltLKzL1Pp2LuH3saM4e0TvB7iTp4DHQS5JyzooNO5h0TwmvLN5Qo/7xE49g4gXD6VxUkFBnknTwGeglSTmjsiry21eW8e3HF7CzvDJT79+9PVPHF3PSwO4JdidJyTDQS5JywsK1W7luRglvrdyUqeUFuPIjA7nmnKG0K8xPrjlJSpCBXpLUrJVVVPGz5xbz42ffobwyZurD+3Ri+oRiivt1Ta45SWoGDPSSpGbr7ZWbmDizhPlrtmZqBfmB/zhrCF88fRCFbfIS7E6SmgcDvSSp2dlZVsn3nlrIr15cQtWeSXmOPrwr0ycUM7R3p+Sak6RmxkAvSWpWXl28gevvKWHZhh2ZWruCfK49bxifOaU/+Xk+IEqSqjPQS5KahS2l5Ux9dD7/+9cVNeqnDu7OnWOLOaJ7+4Q6k6TmzUAvSUrc0/PWcuO9s1mzpTRT61TUhskXjeSK4/oRgrPyklQbA70kKTEbtu3ilgfn8sDb79Wof3Rkb267fDS9Oxcl1Jkk5Q4DvSTpoIsx8sDb73HLg3PZuL0sU+/RsZBbLh3NhWP6OCsvSXVkoJckHVSrN+/kpntn8/T8dTXq4z7Ul8kXj+SQDoUJdSZJuSkrN/ANIUwLITwdQlgZQtgZQtgYQngzhHBzCKFez+EOIfQLIfwmhPBeCGFXCGFZCOH7IYRD9rPNKSGER9LvuyOEUBJCuCaE4GMDJamZqKqK/O9fV3Dud1+oEeYP61LEXZ89nu/+89GGeUlqgGzN0H8VeAN4ElgHdABOAqYAnw8hnBRjXHmgQUIIg4BXgF7A/cB84ATgK8D5IYRTY4wb9trmMmAmUAr8CdgIXAJ8DzgVuCILn0+S1AjL1m9n0j0l/GXJxhr1T518JNedP5yObf2FsSQ1VLb+Bu0cYyzduxhCuB24Abge+HIdxvkpqTB/dYzxR9XG+S6pHxpuB75Yrd4Z+CVQCZwRY/x7uj4ZeAaYEEL4WIzx7oZ+MElSw1VUVnHXy8v4zpMLKC2vytQH9OjAtPHFnDCgW4LdSVLLkJVTbvYV5tP+nF4OOdAYIYSBwLnAMuAne718M7Ad+GQIoUO1+gSgJ3D37jBfrZ+b0t9+6UDvLUnKvvlrtjD+Z69w+yPzMmE+Py/wpTMG8ehXPmKYl6QsaerfcV6SXpbUYd2z0ssnYoxV1V+IMW4NIbxMKvCfBDy91zaP7WO8F4AdwCkhhLYxxl316lyS1CC7Kir5ybOL+emzi6ioipn6iEM7M318MWP6dUmwO0lqebIa6EMI1wIdgS7AccCHSYX5qXXYfFh6ubCW198hFeiHsifQ17pNjLEihLAUGAUMBOYdoPfXa3lp+P62kyTt8eaKD5g4s4SFa7dlaoX5eXzlnCF8/rSBFORn5RfDkqRqsj1Dfy3Qu9r3jwGfiTG+X4dtd0/ZbK7l9d31ro3cRpKUZTvKKvjOEwv5zctLiXsm5Tn2yEOYNr6Ywb06JtecJLVwWQ30McY+ACGE3sAppGbm3wwhXBxjfKORw+9+wkjc71oN3CbGeOw+B0jN3B9Tj/eUpFbllUXrmXTPLFZs3JGptS/M5+vnDeNTJ/cnP88HRElSU2qSc+hjjGuBe0MIb5A6HeZ3wOgDbLZ7Nr22kys777VeQ7eRJGXB5p3l3PnIPO7+W827En9kSA/uGDuGw7u1T6gzSWpdmvSi2Bjj8hDCXODoEEKPGOP6/ay+IL0cWsvru++UU/18+QWkztUfCtQ4Bz6E0AYYAFQAS+rbuySpdk/MWcNN981m3dY99xvoXNSGyRePZMKx/QjBWXlJOlgOxpM8DksvKw+w3rPp5bkhhLzqd7oJIXQi9ZConcBfqm3zDPBx4Hzgj3uNdxrQHnjBO9xIUnas37aLKQ/M4aGS1TXq54/qwzcvH0WvTkUJdSZJrVejbzcQQhgeQuizj3pe+sFSvYBXYowfpOsF6W0GVV8/xrgYeALoD1y113C3kHr67O9ijNur1WcA64GPhRCOq/beRcBt6W9/1pjPJ0mCGCP3vvku53z3+RphvkfHtvzs48fw808ea5iXpIRkY4b+fOBbIYQXgMXABlJ3ujmd1O0i1wBXVlu/L6lbSC4nFd6r+zLwCvDDEMLZ6fVOBM4kdarNjdVXjjFuCSFcSSrYPxdCuBvYCFxK6paWM4A/ZeEzSlKr9d6mndx47yyeXVDzhmUTju3HTReNoGv7woQ6kyRBdgL9U8AvSJ0ScxSpW0RuJxXAfw/8MMa4sS4DxRgXp2fav0nqB4ULgdXAD4Fb9jVOjPG+EMLppML+eKAIWAT8Z/q963NXHElSWlVV5H9eW8G0R+ezbVdFpt63azvuHDeG04b2TLA7SdJujQ70McbZ/OMpMvtbfxl7bie5r9dXAp+tZw8vkwr/kqQsWPL+NibNnMVry/bMo4QAnz65P18/bxgd2h6MS7AkSXXh38iSpIyKyip+9dJSvvfkQnZVZO5NwMCeHZg+vpjj+ndLsDtJ0r4Y6CVJAMx9bwvXzXyb2au2ZGr5eYEvnj6Q/zhrCEUF+Ql2J0mqjYFeklq5XRWV/PiZRfzsucVUVO257GjUYZ2ZPqGYUYfV9uw+SVJzYKCXpFbs9eUfMHFmCYvWbcvUCtvk8dVzhnLlRwbQJr/RdzeWJDUxA70ktULbd1Xw7ScW8NtXllH9XmDH9z+EqeOLGdSzY3LNSZLqxUAvSa3Mi++8z/X3zOLdD3Zmah0K85l0wXA+fuKR5OXVeiMySVIzZKCXpFZi845ybnt4Lv/3+rs16qcP7cntY0fT75D2CXUmSWoMA70ktQKPzV7D5Ptn8/7WXZla1/YFfOPikYz9UF9CcFZeknKVgV6SWrB1W0uZ8sAcHpm1pkb9ouJDmXLJKHp2aptQZ5KkbDHQS1ILFGNk5huruPWhuWzeWZ6p9+zUltsuH815o/ok2J0kKZsM9JLUwrz7wQ5uuHc2Lyx8v0b9n487nBsuHEGX9gUJdSZJagoGeklqIaqqIr//y3KmPTafHWWVmXq/Q9oxdVwxHx7SI8HuJElNxUAvSS3A4ve3MXFGCX9f/kGmFgJ89pQBXHveUNoX+te9JLVU/g0vSTmsvLKKX7ywhB88/Q5lFVWZ+uBeHZk2vphjjzwkwe4kSQeDgV6SctTsVZuZOLOEOe9tydTa5AW+fMYgrjprMG3b5CfYnSTpYDHQS1KOKS2v5IdPv8N/vbCEyqqYqY/p24XpE4oZcWjnBLuTJB1sBnpJyiF/W7aRiTNKWLJ+e6bWtk0e//nRofzbhwfQJj8vwe4kSUkw0EtSDti2q4Lpj83nd68ur1E/YUA3po0vZkCPDgl1JklKmoFekpq55xe+zw33zGLVpp2ZWse2bZh0wXD+3wlHkJcXEuxOkpQ0A70kNVObdpTxzYfmcs8bq2rUzxzWk9vHjuGwru0S6kyS1JwY6CWpGXpk1mq+cf9s1m8ry9QOaV/AlEtHcelRhxGCs/KSpBQDvSQ1I+u2lDL5/tk8PmdtjfolRx3GlEtG0r1j24Q6kyQ1VwZ6SWoGYoz83+vvcttDc9lSWpGp9+7cltsuH8NHR/ZOsDtJUnNmoJekhK3cuIPr75nFS4vW16j/ywmHc/2FI+hcVJBQZ5KkXGCgl6SEVFZFfvfqMqY/toCd5ZWZ+hHd2jN13BhOGdwjwe4kSbnCQC9JCVi0bivXzSjhjRWbMrW8AP966gC+du4w2hXmJ9ecJCmnGOgl6SAqr6zi588t5kfPLKKssipTH9q7I9PGF/OhIw5JsDtJUi4y0EvSQTLr3c18fcbbzF+zNVMryA9cdeZgvnzGYArb5CXYnSQpVxnoJamJlZZX8r2nFvLLF5ZQFffUjzq8K9PHFzOsT6fkmpMk5TwDvSQ1ob8u2cCke2axdP32TK2oII9rzx3GZ08dQH6eD4iSJDWOgV6SmsDW0nKmPTafP/xlRY36yQO7M3X8GI7s3iGhziRJLY2BXpKy7Nn567jx3lm8t7k0U+vUtg03XDSCjx1/OCE4Ky9Jyh4DvSRlycbtZdz60FzufXNVjfo5I3px2+Vj6NOlKKHOJEktmYFekhopxshDJauZ8sAcNmwvy9S7dyhkyqWjuLj4UGflJUlNxkAvSY2wdkspN947m6fmra1Rv/zow/jGJaPo1qEwoc4kSa1FowN9CKE7MBa4CBgD9AXKgFnAXcBdMcaq2kfIjPOZ9Pr7UxVjzDw+MYTQH1i6n/X/FGP82IHeW5LqK8bIn/62ktsfmcfW0opM/dAuRdw+djRnDe+dYHeSpNYkGzP0VwA/A1YDzwIrgN7AOOBXwAUhhCtijLH2IQB4C7illtc+ApwFPFrL628D9+2jPvsA7ylJ9bZiww4m3VPCK4s31Kh//MQjmHTBcDoVFSTUmSSpNcpGoF8IXAo8XH0mPoRwA/AaMJ5UuJ+5v0FijG+RCvX/IITwavqPv6hl87dijFPq07Qk1VdlVeSul5fy7ScWUFq+5xeP/bu3Z+r4Yk4a2D3B7iRJrVWjA32M8Zla6mtCCD8HbgfO4ACBvjYhhNHAScAq4OEGtilJjbJgzVaum1nC2ys3ZWp5Aa78yECuOWco7Qrza99YkqQm1NQXxZanlxX7XWv/vpBe/jrGWFnLOoeFEL4AdAc2AK/GGEsa8Z6SBEBZRRU/fW4RP3l2EeWVe84cHN6nE9MnFFPcr2tyzUmSRBMG+hBCG+BT6W8fa+AY7YBPAFWkzsevzUfTX9W3fQ74dIxxxT63+Mf3er2Wl4bXZXtJLc/bKzdx3YwSFqzdmqkV5ufxH2cN5gunD6KwTV6C3UmSlNKUM/RTgdHAIzHGxxs4xj8BXUmdn79yH6/vAG4ldUHsknStGJgCnAk8HUI4Osa4vYHvL6kV2llWyXefXMCvX1pKVbXL+T90RFemjy9mSO9OyTUnSdJemiTQhxCuBr4GzAc+2YihPp9e/te+XowxrgO+sVf5hRDCucBLwInA54AfHOiNYozH7quenrk/pq4NS8ptry7ewKR7Sli+YUem1q4gn6+fN4xPn9Kf/DwfECVJal6yHuhDCFeRCtBzgbNjjBsbOM5I4BTgXeCR+mwbY6wIIfyKVKA/jToEekmt25bScu58ZD5/fK3mWXqnDu7OnWOLOaJ7+4Q6kyRp/7Ia6EMI1wDfI3X/97PTM+gNVZeLYffn/fSyQyN6kNQKPDV3LTfeN4u1W3Zlap2K2jD5opFccVw/QnBWXpLUfGUt0IcQJpI6b/4t4KMxxvWNGKuI1Kk6VcCvGzjMSenlkv2uJanV2rBtF7c8OJcH3n6vRv3ckb259fLR9O5clFBnkiTVXVYCfQhhMvBN4HXg3P2dZhNCKAAGAeUxxsW1rHYFcAjwUC0Xw+4e60TgzRhj2V71s4Cvpr/9Q50/iKRWIcbIA2+/x5QH5vDBjvJMvUfHQm65dDQXjunjrLwkKWc0OtCHED5NKsxXAi8CV+/jH8JlMcbfpv/cF5gHLAf61zLs7otha3sy7G7TgFHpW1S+m64VA2el/zw5xvjKAT+EpFZj9ead3HTvbJ6eX/OMwHHH9GXyRSM5pENhQp1JktQw2ZihH5Be5gPX1LLO88Bv6zJYCGEE8GHqdjHs74GxwPHABUABsBb4M/DjGOOLdXlPSS1fVVXkj39bwZ2PzGfbrj3PujusSxG3jxvDmcN6JdidJEkN1+hAH2OcQuq+73VdfxlQ6++yY4zz9vf6Xuv+moafYy+plVi2fjuT7inhL0tqng34qZOP5Lrzh9OxbVM/NFuSpKbjv2KSWqyKyip+8/JSvvPEQnZVVGXqA3t0YOr4Yk4Y0C3B7iRJyg4DvaQWad7qLUycWULJu5sztfy8wOdPG8hXzh5CUUF+gt1JkpQ9BnpJLcquikp+8uxifvrsIiqqYqY+8tDOTJ9QzOi+XRLsTpKk7DPQS2ox3ljxARNnlPDOum2ZWmF+Hl85ZwifP20gBfl5CXYnSVLTMNBLynk7yir4zhML+c3LS4l7JuU59shDmDa+mMG9OibXnCRJTcxALymnvbxoPZPuKWHlxp2ZWvvCfK47bxifOrk/eXk+IEqS1LIZ6CXlpM07y7nj4Xn86e81Hyb9kSE9uGPsGA7v1j6hziRJOrgM9JJyzhNz1nDTfbNZt3VXptalXQGTLx7J+GP6so+nVUuS1GIZ6CXljPe37mLKg3N4uGR1jfoFo/twy2Wj6NWpKKHOJElKjoFeUrMXY+S+t1Zxy4Nz2bSjPFPv0bEtt142igvGHJpgd5IkJctAL6lZW7VpJzfeO4vnFrxfoz7h2H7cdNEIurYvTKgzSZKaBwO9pGapqiryP39dztRH57O9rDJT79u1HXeOG8NpQ3sm2J0kSc2HgV5Ss7Pk/W1MmjmL15ZtzNRCgE+f3J+vnzeMDm39q0uSpN38V1FSs1FRWcUvX1zK955aSFlFVaY+qGcHpo0v5rj+3RLsTpKk5slAL6lZmPveFq6b+TazV23J1PLzAl86fRD/ftZgigryE+xOkqTmy0AvKVGl5ZX8+JlF/Pz5xVRUxUx91GGdmT6hmFGHdUmwO0mSmj8DvaTEvL58I9fNKGHx+9sztcI2eXz1nKFc+ZEBtMnPS7A7SZJyg4Fe0kG3fVcF33p8Af/96jLinkl5ju9/CFPHFzOoZ8fkmpMkKccY6CUdVC8sfJ/r75nFqk07M7UOhflMumA4Hz/xSPLyQoLdSZKUewz0kg6KzTvKufXhucx4/d0a9dOH9uSOcWPo27VdQp1JkpTbDPSSmtxjs1cz+f45vL91V6bWtX0B37h4JGM/1JcQnJWXJKmhDPSSmsy6raXcfP8cHp29pkb9ouJDmXLJKHp2aptQZ5IktRwGeklZF2Nk5huruPWhuWzeWZ6p9+zUltsuH815o/ok2J0kSS2LgV5SVq3cuIMb7p3Fi++sr1H/5+MO54YLR9ClfUFCnUmS1DIZ6CVlRVVV5HevLmP64wvYUVaZqR/erR13ji3mw0N6JNidJEktl4FeUqMtWreNSTNL+PvyDzK1EOCzpwzg2vOG0r7Qv2okSWoq/isrqcHKK6v4xQtL+MFT71BWWZWpD+nVkWkTijnmiEMS7E6SpNbBQC+pQWav2sx1M0qYu3pLptYmL/DlMwZx1VmDadsmP8HuJElqPQz0kuqltLySHzz9Dr94YQmVVTFTH9O3C9MnFDPi0M4JdidJUutjoJdUZ39btpGJM0pYsn57pta2TR7/+dGh/NuHB9AmPy/B7iRJap0M9JIOaNuuCqY/Np/fvbq8Rv3EAd2YOr6YAT06JNSZJEky0Evar+cWrOPGe2ezatPOTK1j2zZcf+Fw/uX4I8jLCwl2J0mSDPSS9umD7WXc+vBc7nljVY36WcN7cfvY0RzapV1CnUmSpOoM9JJqiDHyyKw13PzAbNZvK8vUD2lfwJRLR3HpUYcRgrPykiQ1F40O9CGE7sBY4CJgDNAXKANmAXcBd8UYq2ofocZYy4Aja3l5bYyxTy3bnQLcBJwEFAGLgN8AP4oxVu5rG0n/aN2WUm66bzZPzF1bo37JUYcx5ZKRdO/YNqHOJElSbbIxQ38F8DNgNfAssALoDYwDfgVcEEK4IsYYax+ihs3A9/dR37avlUMIlwEzgVLgT8BG4BLge8Cp6f4k7UeMkf/7+7vc+vBctpZWZOq9O7fltsvH8NGRvRPsTpIk7U82Av1C4FLg4eoz8SGEG4DXgPGkwv3MOo63KcY4pS4rhhA6A78EKoEzYox/T9cnA88AE0IIH4sx3l3H95ZanZUbd3D9PbN4adH6GvV/OeEIrr9wOJ2LChLqTJIk1UWjbxodY3wmxvjg3qfVxBjXAD9Pf3tGY9+nFhOAnsDdu8N8+r1LSZ2CA/ClJnpvKadVVkV+89JSzv3eCzXC/JHd2/O/V57InePGGOYlScoBTX1RbHl6WbHftWpqG0L4BHAEsB0oAV6o5Vz4s9LLx/bx2gvADuCUEELbGOOuevQgtWjvrN3KxJklvLFiU6aWF+DfPjyA//zoMNoV5ifXnCRJqpcmC/QhhDbAp9Lf7itw16YP8Pu9aktDCJ+NMT6/V31Yerlw70FijBUhhKXAKGAgMO8A/b5ey0vDD9yylBvKKqr4r+cX86NnFlFWueeXasN6d2LahGKOPrxrcs1JkqQGacoZ+qnAaOCRGOPjddzmLuBFYA6wlVQQ/3fg88CjIYSTY4xvV1u/S3q5uZbxdte71qNvqUUqeXcT180oYf6arZlaQX7gqjMH8+UzBlPYptFn4EmSpAQ0SaAPIVwNfA2YD3yyrtvFGG/ZqzQb+GIIYVt6vCmkbpFZ51Z2D12H9z52nwOkZu6Pqcd7Ss1KaXkl33tyIb98cQlV1f5POOrwrkwfX8ywPp2Sa06SJDVa1gN9COEq4AfAXODsGOPGLAz7c1KB/rS96rtn4Luwb533Wk9qVf6yZAOTZpawbMOOTK2oII9rzx3GZ08dQH6eD4iSJCnXZTXQhxCuIXX/99mkwvy6LA29e5wOe9UXAMcBQ4Ea58Cnz+EfQOqC3CVZ6kPKCVtLy5n66Hz+568ratRPHtidqePHcGT3vf9XkiRJuSprgT6EMJHUefNvAR+NMa7f/xb1cnJ6uXcwfwb4OHA+8Me9XjsNaE/qDjne4UatxrPz13HDvbNYvbk0U+vUtg03XjSCfz7+cEJwVl6SpJYkK4E+/SCnb5KaJT93f6fZhBAKgEFAeYxxcbX6KGD13tuGEI4Efpz+9g97DTcDmAZ8LITwo2oPlioCbkuv87MGfzAph2zcXsY3H5zDfW+9V6N+zohe3Hb5GPp0KUqoM0mS1JQaHehDCJ8mFeYrSd2h5up9zAAuizH+Nv3nvqRuIbkc6F9tnSuASSGEZ4GlpO5yMwi4CCgCHgG+XX3QGOOWEMKVpIL9cyGEu4GNpJ5cOyxd/1NjP6PUnMUYebBkNVMemMPG7WWZevcOhUy5dBQXFx/qrLwkSS1YNmboB6SX+cA1tazzPPDbA4zzLKkQ/iFSp9h0ADYBL5G6L/3vY4z/cLeaGON9IYTTgRuB8aTC/yLgP4Ef7msbqaVYs7mUm+6bzVPz1taoX370YXzjklF061CYUGeSJOlgaXSgjzFOIXU7ybquv4w9t5OsXn+eVPBvSA8vAxc2ZFspF8UYuftvK7nj4Xls3bXnQcyHdini9rGjOWt47wS7kyRJB1NTPlhKUhNYvmE7k2bO4tUlG2rUP3HSEUw8fzidigoS6kySJCXBQC/liMqqyF0vL+XbTyygtLwqU+/fvT1Txxdz0sDuCXYnSZKSYqCXcsCCNVu5bmYJb6/clKnlBbjytIF89ZyhFBXkJ9ecJElKlIFeasbKKqr46XOL+Mmziyiv3HN99/A+nZg+oZjifl2Ta06SJDULBnqpmXpr5SYmzihhwdqtmVphfh7/cdZgvnD6IArb5CXYnSRJai4M9FIzs7Osku8+uYBfv7SUqmo3Xf3QEV2ZPr6YIb07JdecJElqdgz0UjPyyuL1TJo5ixUbd2Rq7Qry+fp5w/j0Kf3Jz/MBUZIkqSYDvdQMbCkt585H5vPH11bUqH94cA/uHDeGw7u1T6gzSZLU3BnopYQ9NXctN943i7VbdmVqnYraMPmikVxxXD9CcFZekiTVzkAvJWTDtl1MeXAuD779Xo36uSN7c+vlo+nduSihziRJUi4x0EsHWYyRB95+jykPzOGDHeWZeo+Ohdxy6WguHNPHWXlJklRnBnrpIHpv005uum82z8xfV6M+7pi+TL5oJId0KEyoM0mSlKsM9NJBUFUV+ePfVnDnI/PZtqsiU+/btR23jx3NGcN6JdidJEnKZQZ6qYktXb+dSTNL+OvSjTXqnzr5SK47fzgd2/q/oSRJajiThNREKiqr+PVLS/nukwvZVVGVqQ/s0YGp44s5YUC3BLuTJEkthYFeagLzVm9h4swSSt7dnKnl5wU+f9pAvnL2EIoK8hPsTpIktSQGeimLdlVU8pNnFvHT5xZTURUz9ZGHdmb6hGJG9+2SYHeSJKklMtBLWfLGig+YOKOEd9Zty9QK2+TxlbOH8PnTBlKQn5dgd5IkqaUy0EuNtKOsgm8/vpC7XllK3DMpz7FHHsK08cUM7tUxueYkSVKLZ6CXGuGld9Zz/b0lrNy4M1NrX5jPxPOH88mTjiQvzwdESZKkpmWglxpg885ybn94Ln/++7s16h8Z0oM7xo7h8G7tE+pMkiS1NgZ6qZ4en7OGyffNZt3WXZlal3YFTL54JOOP6UsIzspLkqSDx0Av1dH7W3cx5YE5PDxrdY36BaP7cMtlo+jVqSihziRJUmtmoJcOIMbIvW+u4psPzWXTjvJMvUfHttx62SguGHNogt1JkqTWzkAv7ceqTTu58d5ZPLfg/Rr1K47tx00XjaRL+4KEOpMkSUox0Ev7UFUV+Z+/Lmfqo/PZXlaZqfft2o47x43htKE9E+xOkiRpDwO9tJfF729j0swS/rbsg0wtBPj0yf35+nnD6NDW/20kSVLzYTKR0ioqq/jFi0v4/lPvUFZRlakP6tmBaeOLOa5/twS7kyRJ2jcDvQTMeW8zE2eWMHvVlkwtPy/wpdMH8e9nDaaoID/B7iRJkmpnoFerVlpeyY+eeYefP7+EyqqYqY/u25lp44sZdViXBLuTJEk6MAO9Wq3Xl2/kuhklLH5/e6ZW2CaPr54zlCs/MoA2+XkJdidJklQ3Bnq1Ott3VfCtxxfw368uI+6ZlOeE/t2YOn4MA3t2TK45SZKkejLQq1V5YeH7XH/PLFZt2pmpdSjMZ9IFw/n4iUeSlxcS7E6SJKn+DPRqFTbtKOO2h+cx4/V3a9RPH9qTO8aNoW/Xdgl1JkmS1DgGerV4j85azeT757B+265MrWv7Ar5x8UjGfqgvITgrL0mSclejA30IoTswFrgIGAP0BcqAWcBdwF0xxqraR2j4OCGE/sDS/Qz7pxjjx+r5kdRCrNtays33z+HR2Wtq1C8qPpQpl4yiZ6e2CXUmSZKUPdmYob8C+BmwGngWWAH0BsYBvwIuCCFcEWP1yw+zPs7bwH37qM+u96dRzosxMuP1d7nt4Xls3lmeqffq1JZbLx/NeaP6JNidJElSdmUj0C8ELgUerj6DHkK4AXgNGE8qlM9swnHeijFOacRnUAuxcuMObrh3Fi++s75G/Z+PO5wbLhpBl3YFCXUmSZLUNBod6GOMz9RSXxNC+DlwO3AGBwj02RpHrVNVVeR3ry5j+uML2FFWmakf3q0dU8cVc+rgHgl2J0mS1HSa+qLY3ec7VDTxOIeFEL4AdAc2AK/GGEvq8wYhhNdreWl4fcbRwbdo3VYmzpzF68s/yNRCgM+eMoBrzxtK+0Kv/ZYkSS1XkyWdEEIb4FPpbx9r4nE+mv6qvt1zwKdjjCsa+t5q3sorq/jFC0v4wVPvUFa553rpIb06Mm1CMccccUiC3UmSJB0cTTl1ORUYDTwSY3y8icbZAdxK6oLYJelaMTAFOBN4OoRwdIxx+4HeJMZ47L7q6Zn7YxrUuZrM7FWbuW5GCXNXb8nU2uQFvnzmYK46cxBt2+Qn2J0kSdLB0ySBPoRwNfA1YD7wyaYaJ8a4DvjGXuUXQgjnAi8BJwKfA37Q0B7UvJSWV/KDp9/hFy8sobJqzw2Pivt1Ydr4YkYc2jnB7iRJkg6+rAf6EMJVpAL0XODsGOPGgz1OjLEihPArUoH+NAz0LcJrSzcyaWYJS9bv+YVL2zZ5fO3cofzrqQNok5+XYHeSJEnJyGqgDyFcA3yP1P3fz07PoCc1zvvpZYeG9KDmY9uuCqY9Op/f/2V5jfqJA7oxbXwx/Xv4n1iSJLVeWQv0IYSJpM53fwv4aIxx/f63aNpxgJPSyyX7XUvN2rML1nHjPbN4b3NpptaxbRuuv3A4/3L8EeTlhQS7kyRJSl5WAn0IYTLwTeB14Nz9nR4TQigABgHlMcbFDR0nvf6JwJsxxrK96mcBX01/+4d6fhw1Ax9sL+PWh+Zyz5uratTPGt6L28eO5tAu7RLqTJIkqXlpdKAPIXyaVAivBF4Erg7hH2ZNl8UYf5v+c19gHrAc6N+IcQCmAaPSt6h8N10rBs5K/3lyjPGVhn0yJSHGyCOz1nDzA7NZv23Pz2ndOhRy8yUjufSow9jHcSFJktRqZWOGfkB6mQ9cU8s6zwO/bYJxfg+MBY4HLgAKgLXAn4EfxxhfPMB7qhlZu6WUyffN5om5a2vULz3qMG6+ZCTdO7ZNqDNJkqTmq9GBPsY4hdR93+u6/jLgH6ZY6ztOeptfA7+uzzZqfmKM/PnvK7nt4XlsLd3zMOA+nYu47fLRnDOyd4LdSZIkNW9N+WAp6YBWbNjB9feW8PKiDTXq/+/EI5h0wXA6FxUk1JkkSVJuMNArEZVVkd++soxvP76AneWVmfqR3dtz57gxnDKoR4LdSZIk5Q4DvQ66d9Zu5bqZJby5YlOmlhfgcx8ZyFfPGUq7wvzkmpMkScoxBnodNGUVVfz8+cX8+JlFlFVWZerDendi2oRijj68a3LNSZIk5SgDvQ6Kknc3cd2MEuav2ZqpFeQH/v3MIXzpjEEUtslLsDtJkqTcZaBXk9pZVsn3n1rIL19cQlXcUz/q8K5MH1/MsD6dkmtOkiSpBTDQq8n8ZckGJs0sYdmGHZlaUUEe1547jM+eOoD8PB8QJUmS1FgGemXd1tJypj46n//564oa9VMGdWfquGKO6N4+oc4kSZJaHgO9suqZ+Wu58d7ZrN5cmql1atuGGy8awT8ffzghOCsvSZKUTQZ6ZcWGbbv45kNzuf+t92rUzxnRm9suH02fLkUJdSZJktSyGejVKDFGHixZzZQH5rBxe1mm3r1DIVMuHcXFxYc6Ky9JktSEDPRqsDWbS7npvlk8NW9djfrYD/Vl8sUj6dahMKHOJEmSWg8Dveotxsjdf1vJHQ/PY+uuikz90C5F3DF2DGcO75Vgd5IkSa2LgV71snzDdibNnMWrSzbUqH/ipCOYeP5wOhUVJNSZJElS62SgV51UVkXuenkp335iAaXlVZn6gB4dmDpuDCcO7J5gd5IkSa2XgV4HtGDNVq6bWcLbKzdlankBrjxtIF89ZyhFBfnJNSdJktTKGehVq7KKKn7y7CJ++twiyitjpj68TyemTyimuF/X5JqTJEkSYKBXLd5auYnrZrzNwrXbMrXC/DyuPnswXzh9EAX5eQl2J0mSpN0M9KphZ1kl33liAb95eSlVeyblOeaIrkyfUMzgXp2Sa06SJEn/wECvjFcWr2fSzFms2LgjU2tXkM915w/jUyf3Jz/PB0RJkiQ1NwZ6saW0nDsfmccfX1tZo/7hwT24c9wYDu/WPqHOJEmSdCAG+lbuyblruem+WazdsitT61zUhpsuHskVx/YjBGflJUmSmjMDfSu1ftsupjwwh4dKVteonzeqN7deNppenYsS6kySJEn1YaBvZWKM3P/We9zy4Bw+2FGeqffo2JZbLxvFBWMOTbA7SZIk1ZeBvhV5b9NObrpvNs/MX1ejPv6Yfky+eARd2xcm1JkkSZIaykDfClRVRf73tRVMfXQ+23ZVZOp9u7bjjnFjOH1ozwS7kyRJUmMY6Fu4peu3M2lmCX9dujFTCwE+ddKRfP384XRs6yEgSZKUy0xzLVRFZRW/fmkp331yIbsqqjL1gT07MG18Mcf375Zgd5IkScoWA30LNPe9LUycWcKsVZsztfy8wBdOG8jVZw+hqCA/we4kSZKUTQb6FmRXRSU/fmYRP3tuMRVVMVMfeWhnpk8oZnTfLgl2J0mSpKZgoG8hXl/+ARNnlrBo3bZMrbBNHl85ewifP20gBfl5CXYnSZKkpmKgz3E7yir41uML+O0ry4h7JuU57shDmDq+mMG9OibXnCRJkpqcgT6HvfTOeibdU8K7H+zM1DoU5jPxguF84sQjycsLCXYnSZKkg8FAn4M27yjn9kfm8ue/v1ujftrQntwxdjT9DmmfUGeSJEk62Bp9YnUIoXsI4XMhhHtDCItCCDtDCJtDCC+FEP4thFCv9wgh9Ash/CaE8F4IYVcIYVkI4fshhEP2s80pIYRHQggbQwg7QgglIYRrQggt7nYuj81ewznfe75GmO/SroDvXHEU//3Z4w3zkiRJrUw2ZuivAH4GrAaeBVYAvYFxwK+AC0IIV8RY/QzvfQshDAJeAXoB9wPzgROArwDnhxBOjTFu2Guby4CZQCnwJ2AjcAnwPeDUdH857/2tu5jywBwenrW6Rv3CMX245dLR9OzUNqHOJEmSlKRsBPqFwKXAwzHGzBOMQgg3AK8B40mF+5l1GOunpML81THGH1Ub67vAV4HbgS9Wq3cGfglUAmfEGP+erk8GngEmhBA+FmO8u1GfMEExRu55YxXffGgum3eWZ+o9O7Xl1stGcf7oQxPsTpIkSUlr9Ck3McZnYowPVg/z6foa4Ofpb8840DghhIHAucAy4Cd7vXwzsB34ZAihQ7X6BKAncPfuMJ9+71LgpvS3X6rzh2lmVm3ayWfu+htf+7+3a4T5fzquH0999XTDvCRJkpr8otjdKbSiDuuelV4+sY8fDraGEF4mFfhPAp7ea5vH9jHeC8AO4JQQQtsY4656dZ6gqqrIH/66nGmPzmd7WWWm3u+QdkwdV8yHh/RIsDtJkiQ1J00W6EMIbYBPpb/dV+De27D0cmEtr79DKtAPZU+gr3WbGGNFCGEpMAoYCMw7QL+v1/LS8P1t1xS+//Q7/PDpdzLfhwCfOaU/1547jA5tvTGRJEmS9mjKx4dOBUYDj8QYH6/D+l3Sy821vL673rWR2zR7nzjpCLq0KwBgcK+OzPjiKdx8ySjDvCRJkv5BkyTEEMLVwNdI3aXmk9kaNr084N1yGrJNjPHYfQ6Qmrk/ph7v2Wi9OhUx5dKRLHl/O/9+1mDatmlxd9+UJElSlmQ90IcQrgJ+AMwFzo4xbqzjprtn07vU8nrnvdZr6DY5YeyH+iXdgiRJknJAVk+5CSFcA/wYmA2cmb7TTV0tSC+H1vL6kPSy+vnytW6TPod/AKkLcpfUow9JkiQpZ2Qt0IcQJpJ6mNNbpML8unoO8Wx6ee7eT5cNIXQi9ZConcBfqr30THp5/j7GOw1oD7ySS3e4kSRJkuojK4E+/SCnqcDrpE6zWb+fdQtCCMPTT4XNiDEuBp4A+gNX7bXZLUAH4Hcxxu3V6jOA9cDHQgjHVXuPIuC29Lc/a9CHkiRJknJAo8+hDyF8Gvgmqae1vghcHULYe7VlMcbfpv/cl9QtJJeTCu/VfRl4BfhhCOHs9HonAmeSOtXmxuorxxi3hBCuJBXsnwsh3A1sJPXk2mHp+p8a+xklSZKk5iobF8UOSC/zgWtqWed54LcHGijGuDg90/5NUqfRXAisBn4I3LKvC2xjjPeFEE4nFfbHA0XAIuA/gR/GGOtzVxxJkiQppzQ60McYpwBT6rH+MvbcTnJfr68EPlvPHl4mFf4lSZKkVqUpHywlSZIkqYkZ6CVJkqQcZqCXJEmScpiBXpIkScphBnpJkiQphxnoJUmSpBxmoJckSZJymIFekiRJymEGekmSJCmHGeglSZKkHBZijEn30KyFEDa0a9eu24gRI5JuRZIkSS3YvHnz2Llz58YYY/f6bGegP4AQwlKgM7DsIL/18PRy/kF+31zl/qo/91n9uL/qx/1VP+6v+nF/1Y/7q36S3F/9gS0xxgH12chA30yFEF4HiDEem3QvucD9VX/us/pxf9WP+6t+3F/14/6qH/dX/eTi/vIcekmSJCmHGeglSZKkHGaglyRJknKYgV6SJEnKYQZ6SZIkKYd5lxtJkiQphzlDL0mSJOUwA70kSZKUwwz0kiRJUg4z0EuSJEk5zEAvSZIk5TADvSRJkpTDDPSSJElSDjPQH0QhhH4hhN+EEN4LIewKISwLIXw/hHBIEuM0d9n4nOltYi1fa5qy/4MphDAhhPCjEMKLIYQt6c/3hwaO1eKPr2ztr9ZwfIUQuocQPhdCuDeEsCiEsDOEsDmE8FII4d9CCPX6d6SlH1/Z3F+t4fjaLYQwLYTwdAhhZXqfbQwhvBlCuDmE0L2eY7XoYwyyt79a0zFWXQjhk9U+5+fquW2zPL58sNRBEkIYBLwC9ALuB+YDJwBnAguAU2OMGw7WOM1dFvfXMqAr8P19vLwtxvjt7HScrBDCW8BRwDbgXWA48D8xxk/Uc5zWcny9RXb21zJa+PEVQvgi8DNgNfAssALoDYwDugAzgStiHf4xaQ3HV5b31zJa+PG1WwihDHgDmAusAzoAJwHHAe8BJ8UYV9ZhnBZ/jEFW99cyWskxtlsI4XBgFpAPdASujDH+qo7bNt/jK8bo10H4Ah4HIvAfe9W/m67//GCO09y/sri/lgHLkv48B2F/nQkMAQJwRnof/SGp/d7cv7K4v1r88QWcBVwC5O1V70MqrEZgfB3HavHHV5b3V4s/vqp91qJa6ren99lP6zhOiz/Gsry/Ws0xlv68AXgKWAx8K72vPleP7Zvt8ZX4zm0NX8DA9H/opfv4S74TqVnC7UCHgzFOc//K5udsbX9ZpT9zgwJqazm+srW/0tu2uuNrr89/Q3rf/agO67bK46uh+yu9fqs+vtL74Kj0PnuyDut6jNVjf6XXb1XHGPAVoAo4DZhSn0Df3I8vz6E/OM5KL5+IMVZVfyHGuBV4GWhP6tdlB2Oc5i7bn7NtCOETIYQbQghfCSGcGULIz2K/LUVrOb6yrTUfX+XpZUUd1vX4qt/+2q01H1+Q+m0HQEkd1vUYq9/+2q1VHGMhhBHAVOAHMcYXGjBEsz6+2iTxpq3QsPRyYS2vvwOcCwwFnj4I4zR32f6cfYDf71VbGkL4bIzx+Ya12CK1luMr21rl8RVCaAN8Kv3tY3XYpFUfXw3YX7u1quMrhHAtqfOau5A6H/zDpMLp1Dps3uqOsUbur91a/DGW/v/v96ROe7uhgcM06+PLGfqDo0t6ubmW13fXux6kcZq7bH7Ou4CzSf2F1QEYA/wX0B94NIRwVIO7bHlay/GVTa35+JoKjAYeiTE+Xof1W/vxVd/9Ba3z+LoWuBm4hlQ4fQw4N8b4fh22bY3HWGP2F7SeY+wbwIeAz8QYdzZwjGZ9fBnom4eQXjb2lkPZGqe5q/PnjDHeEmN8Jsa4Nsa4I8Y4O8b4RVIXsLQjdQ6d6qa1HF911lqPrxDC1cDXSN3h4ZPZGja9bHHHV0P3V2s8vmKMfWKMgVTAHEfqvOU3QwjHZGH4FneMNXZ/tYZjLIRwAqlZ+e/EGF9tyrdKLxM5vgz0B8fun9q61PJ6573Wa+pxmruD8Tl/nl6e1ogxWprWcnwdDC32+AohXAX8gNTt8s6MMW6s46at8vhqxP7anxZ7fO2WDpj3kjqFoTvwuzps1iqPMWjw/tqfFnGMVTvVZiEwuZHDNevjy0B/cCxIL4fW8vqQ9LK287KyPU5zdzA+57r0skMjxmhpWsvxdTC0yOMrhHAN8GNgNqlwWp8Hz7S646uR+2t/WuTxtS8xxuWkfhgaFULocYDVW90xtrd67q/9aSnHWEdSx8MIoLT6g7NInaoE8Mt07fsHGKtZH19eFHtwPJtenhtCyKt+dXQIoRNwKrAT+MtBGqe5Oxif8+T0ckkjxmhpWsvxdTC0uOMrhDCR1HngbwEfjTGur+cQrer4ysL+2p8Wd3wdwGHpZeUB1mtVx9h+1HV/7U9LOcZ2Ab+u5bVjSJ1X/xKpsH6g03Ga9fHlDP1BEGNcDDxB6iKTq/Z6+RZSPwH/Lsa4HSCEUBBCGJ5+IlmDx8lV2dpfIYRRIYRue48fQjiS1KwZwB+y3H6z19qPr/ry+IIQwmRS4fR14Oz9hVOPr+zsr1Z2fA0PIfTZRz0vhHA7qadyvhJj/CBdb9XHWLb2V2s4xmKMO2OMn9vXF/BAerX/Ttf+BLl7fIUYW8y1Ic3aPh4XPA84kdQTKxcCp8T044JDCP1JPbhgeYyxf0PHyWXZ2F8hhCnAJFI/VS8FtgKDgIuAIuARYGyMsexgfKamFEK4HLg8/W0f4DxSMysvpmvrY4zXptftj8fX5TRyf7WW4yuE8Gngt6Rm+37Evs8PXRZj/G16/f604uMrW/urtRxfkDk16VvAC6Se4LkB6A2cTuoizzWkfjCam16/P637GLuGLOyv1nSM7Uv6898MXBlj/FW1en9y8fiKzeDJXa3lCzic1C2iVgNlwHJSF0t122u9/qSukl7WmHFy/aux+4vUX25/JHVniU2kHuryPvAkqftBh6Q/Yxb31ZT0Pqjta1m1dVv98ZWN/dVajq867KsIPOfxld391VqOr/RnHQ38hNTpSetJPXhrM/C39P7038gm2F+t6RirZT/u/n/1c3vVc/L4coZekiRJymGeQy9JkiTlMAO9JEmSlMMM9JIkSVIOM9BLkiRJOcxAL0mSJOUwA70kSZKUwwz0kiRJUg4z0EuSJEk5zEAvSZIk5TADvSRJkpTDDPSSJElSDjPQS5IkSTnMQC9JkiTlMAO9JEmSlMMM9JIkSVIOM9BLkhothHBfCCGGEP5jH6/dmn7tV0n0JkktXYgxJt2DJCnHhRC6AW8CvYGTY4xvputnA08A84HjY4w7kutSklomA70kKStCCKcAzwNLgWOA9sDbQBdSYX5Ogu1JUovlKTeSpKyIMb4CTAaGAP8F/AHoA1xtmJekpuMMvSQpa0IIAXgUOC9d+mOM8f8l2JIktXjO0EuSsiamZonurVb6fkKtSFKr4Qy9JClrQghDgDeAclLnzs8BTogxlibamCS1YM7QS5KyIoTQFvgT0AH4GHAnMAZn6SWpSRnoJUnZ8m3gQ8D0GOMTwM3Ay8AXQgj/lGhnktSCecqNJKnRQgiXkzp3/q/Ah2OMFen64cBbQBvgQzHGJUn1KEktlYFektQoIYQjSIX2PFKhfeler18G3Af8jVTYLzvYPUpSS2aglyRJknKY59BLkiRJOcxAL0mSJOUwA70kSZKUwwz0kiRJUg4z0EuSJEk5zEAvSZIk5TADvSRJkpTDDPSSJElSDjPQS5IkSTnMQC9JkiTlMAO9JEmSlMMM9JIkSVIOM9BLkiRJOcxAL0mSJOUwA70kSZKUwwz0kiRJUg4z0EuSJEk57P8Dm/W+2Pz4J6wAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "image/png": { "height": 261, "width": 378 }, "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "da_gpu.plot()" ] }, { "cell_type": "markdown", "id": "cf09471d-6df5-4bd4-b78d-3ab5ffefc371", "metadata": {}, "source": [ "## Apply custom kernels with apply_ufunc\n", "\n", "(This kernel was copied from [this NCAR tutorial](https://github.com/NCAR/GPU_workshop/blob/workshop/12_CuPyAndLegate/12_CuPyAndLegate.ipynb))" ] }, { "cell_type": "code", "execution_count": 22, "id": "50389d84-49e6-44a9-9bd8-679457752280", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 0., 2., 4.],\n", " [ 0., 4., 10.]], dtype=float32)" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = cp.arange(6, dtype=\"f\").reshape(2, 3)\n", "y = cp.arange(3, dtype=\"f\")\n", "\n", "kernel = cp.ElementwiseKernel(\n", " \"float32 x, float32 y\",\n", " \"float32 z\",\n", " \"\"\"\n", " if (x - 2 > y) {\n", " z = x * y;\n", " } else {\n", " z = x + y;\n", " }\n", " \"\"\",\n", " \"my_kernel\",\n", ")\n", "\n", "kernel(x, y)" ] }, { "cell_type": "markdown", "id": "da96291a-28b3-468d-abdb-98392c050e00", "metadata": {}, "source": [ "We can apply these and other custom kernels using `xarray.apply_ufunc`" ] }, { "cell_type": "code", "execution_count": 23, "id": "4fdc7331-0e1a-415c-a2c0-ca172e9e2573", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "is_gpu: True\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (a: 2, b: 3)>\n",
       "array([[ 0.,  2.,  4.],\n",
       "       [ 0.,  4., 10.]], dtype=float32)\n",
       "Dimensions without coordinates: a, b
" ], "text/plain": [ "\n", "array([[ 0., 2., 4.],\n", " [ 0., 4., 10.]], dtype=float32)\n", "Dimensions without coordinates: a, b" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "xda = xr.DataArray(x, dims=(\"a\", \"b\"))\n", "yda = xr.DataArray(y, dims=(\"b\"))\n", "result = xr.apply_ufunc(\n", " kernel,\n", " xda,\n", " yda,\n", ")\n", "print(\"is_gpu:\", result.cupy.is_cupy)\n", "result" ] } ], "metadata": { "kernelspec": { "display_name": "Python [conda env:miniconda3-gpu]", "language": "python", "name": "conda-env-miniconda3-gpu-py" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.5" }, "widgets": { "application/vnd.jupyter.widget-state+json": { "state": {}, "version_major": 2, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 5 }