{ "cells": [ { "cell_type": "markdown", "id": "c9cffae3-8d9c-4da7-a70d-350e185215fe", "metadata": {}, "source": [ "# Getting to know the flexibility of ibicus: modifying and adjusting debiasers" ] }, { "cell_type": "markdown", "id": "c0679686-c320-4d2e-843d-45b3e2e59faa", "metadata": {}, "source": [ "The following notebook shows how to adjust parameters to modify debiasers. \n", "\n", "For in-depth information about the bias adjustment methods and their usage refer to the documentation that can be found under - API reference -> debias-module." ] }, { "cell_type": "code", "execution_count": 1, "id": "7a4ad5dd-ed67-48d1-8302-67656e8d39b7", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "id": "bd6c5a1f-a7d9-40cc-b39c-2799252acaec", "metadata": { "jp-MarkdownHeadingCollapsed": true, "tags": [] }, "source": [ "## 0. Read in data" ] }, { "cell_type": "markdown", "id": "d72bdb59-914a-43a9-bed9-f43c42c98d09", "metadata": {}, "source": [ "We start by reading in and preprocessing some data. For an explanation of the steps please refer to the [01 Getting started](https://nbviewer.org/github/ecmwf-projects/ibicus/blob/main/notebooks/01%20Getting%20Started.ipynb) notebook." ] }, { "cell_type": "code", "execution_count": 2, "id": "2830b3d6-44f3-4008-bfbb-76f5a4ea057c", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "def get_data(variable, data_path = \"testing_data/\"):\n", " # Load in the data \n", " data = np.load(f\"{data_path}{variable}.npz\", allow_pickle = True)\n", " # Return arrays\n", " return data[\"obs\"], data[\"cm_hist\"], data[\"cm_future\"], data[\"time_obs\"], data[\"time_cm_hist\"], data[\"time_cm_future\"]" ] }, { "cell_type": "markdown", "id": "527787a6", "metadata": {}, "source": [ "We work with daily mean near-surface temperature ('tas') and precipitation flux ('pr') in this notebook. Lets get the testing data for these two variables:" ] }, { "cell_type": "code", "execution_count": 3, "id": "04fc1aca-8029-4db0-8ab9-162fc87397d2", "metadata": {}, "outputs": [], "source": [ "tas_obs, tas_cm_hist, tas_cm_future, tas_time_obs, tas_time_cm_hist, tas_time_cm_future = get_data(\"tas\")\n", "pr_obs, pr_cm_hist, pr_cm_future, pr_time_obs, pr_time_cm_hist, pr_time_cm_future = get_data(\"tas\")" ] }, { "cell_type": "markdown", "id": "a5ab11b8-f183-48c3-928d-5fa727aa953d", "metadata": { "tags": [] }, "source": [ "## 1. A first example" ] }, { "cell_type": "markdown", "id": "77b1a5e0-5b5b-4ca7-b219-3d528ef50d5f", "metadata": {}, "source": [ "As shown in the [01 Getting Started](https://nbviewer.org/github/ecmwf-projects/ibicus/blob/main/notebooks/01%20Getting%20Started.ipynb) notebook, each debiaser is a subclass of the abstract `Debiaser` class. This provides a unified interface for initialising and applying debiasers. Let's read in some temperature data:" ] }, { "cell_type": "markdown", "id": "759ec3a4-120f-42a7-bfec-1b672865fb85", "metadata": {}, "source": [ "We can initialise a debiaser for `\"tas\"` by using the `from_variable` method. If we want to apply ISIMIP and Equidistant-CDF-Matching (ECDFM) we can write:" ] }, { "cell_type": "code", "execution_count": 4, "id": "c6d33b6a-de6c-48aa-9c74-21cd77b23f5f", "metadata": {}, "outputs": [], "source": [ "from ibicus.debias import ECDFM, ISIMIP, DeltaChange\n", "\n", "tas_debiaser_ECDFM = ECDFM.from_variable(\"tas\")\n", "tas_debiaser_ISIMIP = ISIMIP.from_variable(\"tas\")" ] }, { "cell_type": "markdown", "id": "786a6fea-3ef7-435c-9b92-cd6078355aee", "metadata": {}, "source": [ "Applying the debiasers works using the `apply`-classmethod." ] }, { "cell_type": "code", "execution_count": 5, "id": "4afeddf1-46e9-4d06-95b1-3e558a85c95d", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ " 37%|█████████████████████████████████████████████████████████████▌ | 84/225 [00:46<01:15, 1.87it/s]/home/jakobwes/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:684: RuntimeWarning: invalid value encountered in sqrt\n", " sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)\n", " 39%|███████████████████████████████████████████████████████████████▊ | 87/225 [00:48<01:13, 1.87it/s]/home/jakobwes/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:684: RuntimeWarning: invalid value encountered in sqrt\n", " sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)\n", " 39%|████████████████████████████████████████████████████████████████▌ | 88/225 [00:48<01:12, 1.90it/s]/home/jakobwes/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:684: RuntimeWarning: invalid value encountered in sqrt\n", " sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)\n", " 41%|████████████████████████████████████████████████████████████████████▏ | 93/225 [00:51<01:06, 1.99it/s]/home/jakobwes/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:684: RuntimeWarning: invalid value encountered in sqrt\n", " sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)\n", " 51%|███████████████████████████████████████████████████████████████████████████████████▊ | 115/225 [01:03<01:03, 1.72it/s]/home/jakobwes/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:684: RuntimeWarning: invalid value encountered in sqrt\n", " sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)\n", " 57%|█████████████████████████████████████████████████████████████████████████████████████████████▎ | 128/225 [01:10<00:49, 1.95it/s]/home/jakobwes/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:684: RuntimeWarning: invalid value encountered in sqrt\n", " sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)\n", " 62%|█████████████████████████████████████████████████████████████████████████████████████████████████████▎ | 139/225 [01:17<00:53, 1.59it/s]/home/jakobwes/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:684: RuntimeWarning: invalid value encountered in sqrt\n", " sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)\n", " 66%|███████████████████████████████████████████████████████████████████████████████████████████████████████████▉ | 148/225 [01:23<00:56, 1.36it/s]/home/jakobwes/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:684: RuntimeWarning: invalid value encountered in sqrt\n", " sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)\n", " 68%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 152/225 [01:26<00:48, 1.50it/s]/home/jakobwes/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:684: RuntimeWarning: invalid value encountered in sqrt\n", " sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)\n", " 72%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▎ | 161/225 [01:31<00:35, 1.79it/s]/home/jakobwes/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:684: RuntimeWarning: invalid value encountered in sqrt\n", " sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)\n", " 75%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▍ | 168/225 [01:35<00:31, 1.84it/s]/home/jakobwes/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:684: RuntimeWarning: invalid value encountered in sqrt\n", " sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)\n", " 76%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▉ | 170/225 [01:36<00:28, 1.95it/s]/home/jakobwes/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:684: RuntimeWarning: invalid value encountered in sqrt\n", " sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)\n", " 78%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▎ | 176/225 [01:39<00:27, 1.81it/s]/home/jakobwes/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:684: RuntimeWarning: invalid value encountered in sqrt\n", " sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)\n", " 79%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████ | 177/225 [01:40<00:25, 1.85it/s]/home/jakobwes/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:684: RuntimeWarning: invalid value encountered in sqrt\n", " sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)\n", " 81%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▋ | 182/225 [01:42<00:22, 1.91it/s]/home/jakobwes/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:684: RuntimeWarning: invalid value encountered in sqrt\n", " sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)\n", " 82%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 185/225 [01:44<00:19, 2.10it/s]/home/jakobwes/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:684: RuntimeWarning: invalid value encountered in sqrt\n", " sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)\n", " 84%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████ | 188/225 [01:45<00:19, 1.94it/s]/home/jakobwes/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:684: RuntimeWarning: invalid value encountered in sqrt\n", " sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)\n", " 88%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▌ | 197/225 [01:50<00:13, 2.13it/s]/home/jakobwes/anaconda3/lib/python3.9/site-packages/scipy/optimize/_minpack_py.py:175: RuntimeWarning: The iteration is not making good progress, as measured by the \n", " improvement from the last ten iterations.\n", " warnings.warn(msg, RuntimeWarning)\n", " 88%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████ | 199/225 [01:51<00:11, 2.29it/s]/home/jakobwes/anaconda3/lib/python3.9/site-packages/scipy/optimize/_minpack_py.py:175: RuntimeWarning: The iteration is not making good progress, as measured by the \n", " improvement from the last ten iterations.\n", " warnings.warn(msg, RuntimeWarning)\n", " 89%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 200/225 [01:51<00:11, 2.23it/s]/home/jakobwes/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:684: RuntimeWarning: invalid value encountered in sqrt\n", " sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)\n", " 91%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▋ | 204/225 [01:53<00:10, 2.04it/s]/home/jakobwes/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:684: RuntimeWarning: invalid value encountered in sqrt\n", " sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)\n", " 93%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▎ | 209/225 [01:56<00:07, 2.11it/s]/home/jakobwes/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:684: RuntimeWarning: invalid value encountered in sqrt\n", " sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)\n", " 95%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▉ | 214/225 [01:58<00:04, 2.32it/s]/home/jakobwes/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:684: RuntimeWarning: invalid value encountered in sqrt\n", " sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)\n", " 96%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▍ | 216/225 [01:58<00:03, 2.64it/s]/home/jakobwes/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:684: RuntimeWarning: invalid value encountered in sqrt\n", " sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)\n", " 99%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 222/225 [02:01<00:01, 2.03it/s]/home/jakobwes/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:684: RuntimeWarning: invalid value encountered in sqrt\n", " sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)\n", "100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 225/225 [02:03<00:00, 1.82it/s]\n" ] } ], "source": [ "tas_ECDFM = tas_debiaser_ECDFM.apply(tas_obs, tas_cm_hist, tas_cm_future)" ] }, { "cell_type": "markdown", "id": "edc2f6a9-a396-4f7e-bfe8-1d36ec057bd1", "metadata": {}, "source": [ "Some debiasers like ISIMIP require additional arguments in the apply-method like the dates:" ] }, { "cell_type": "code", "execution_count": 6, "id": "b0c24d55-8a4f-41fc-b430-6e8f6ef7ec28", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 225/225 [10:18<00:00, 2.75s/it]\n" ] } ], "source": [ "tas_ISIMIP = tas_debiaser_ISIMIP.apply(tas_obs, tas_cm_hist, tas_cm_future, time_obs = tas_time_obs, time_cm_hist = tas_time_cm_hist, time_cm_future =tas_time_cm_future)" ] }, { "cell_type": "markdown", "id": "46824882-f503-44f7-ad22-2c0cc3a8c60e", "metadata": {}, "source": [ "We can compare the bias corrected output using the evaluation framework. In this notebook, we limit ourselves to plotting the distribution of the debiased data at location [1,1]:" ] }, { "cell_type": "code", "execution_count": 7, "id": "a7de0bf6-e9b2-4393-8219-d2ae992c75c5", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj0AAAGdCAYAAAD5ZcJyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAwr0lEQVR4nO3de3AUZb7G8WfIPSEZSCIZZgkYNSISRESFoCuwXIIC0eMeowsnK8IiHgSMgGB0XdHaTRbOCiiUrHghHBDRWg2H8oKEXQxLRQSCUYIIXiIXSYy6YUIgTgLp8wdF1w4BJDC5vt9P1VQ53b/ufjtdzTy+/Xa3w7IsSwAAAG1cu+ZuAAAAQFMg9AAAACMQegAAgBEIPQAAwAiEHgAAYARCDwAAMAKhBwAAGIHQAwAAjBDY3A1oLHV1dTp06JAiIyPlcDiauzkAAOA8WJalI0eOyO12q107//bNtNnQc+jQIcXHxzd3MwAAwAU4cOCAunTp4td1ttnQExkZKenkHy0qKqqZWwMAAM5HZWWl4uPj7d9xf2qzoefUJa2oqChCDwAArUxjDE1hIDMAADACoQcAABiB0AMAAIzQZsf0AABwJpZl6fjx4zpx4kRzN8VIAQEBCgwMbJbHyRB6AADGqKmpUWlpqY4dO9bcTTFaeHi4OnfurODg4CbdLqEHAGCEuro6lZSUKCAgQG63W8HBwTy8tolZlqWamhp9//33KikpUWJiot8fQHguhB4AgBFqampUV1en+Ph4hYeHN3dzjBUWFqagoCDt27dPNTU1Cg0NbbJtNzhebdq0SaNHj5bb7ZbD4dCaNWvOWjtp0iQ5HA4tXLjQZ7rX69XUqVMVGxuriIgIpaam6uDBgz41FRUVSk9Pl9PplNPpVHp6ug4fPtzQ5gIA4KMpexZwZs11DBq81aNHj6p3795avHjxOevWrFmjjz76SG63u968jIwM5ebmavXq1dq8ebOqqqo0atQon0FlY8aMUVFRkdatW6d169apqKhI6enpDW0uAACApAu4vHXrrbfq1ltvPWfNt99+qylTpuj999/XyJEjfeZ5PB69/PLLWrFihYYOHSpJWrlypeLj47VhwwalpKRo9+7dWrdunbZs2aJ+/fpJkl588UUlJydrz5496t69e0ObDQAADOf3MT11dXVKT0/XI488op49e9abX1hYqNraWg0fPtye5na7lZSUpIKCAqWkpOjDDz+U0+m0A48k9e/fX06nUwUFBWcMPV6vV16v1/5eWVnp5z0DALRVC/L2Ntm2Hh52ZZNtC778flFt7ty5CgwM1LRp0844v6ysTMHBwerYsaPP9Li4OJWVldk1nTp1qrdsp06d7JrTZWdn2+N/nE4nb1gHALQZ48aNk8PhqPcZMWKEXfPxxx/rrrvuUlxcnEJDQ3XllVdq4sSJ2rv3ZKD75ptvfJaNjIxUz5499eCDD+qLL77w2V5OTs4Zt/fSSy/5zO/Ro0e9tr7xxhtyOBy69NJLG+8PcoH8GnoKCwv17LPP2n+MhrAsy2eZMy1/es2/y8zMlMfjsT8HDhxoWOMBAGjBRowYodLSUp/Pa6+9Jkl6++231b9/f3m9Xr366qvavXu3VqxYIafTqSeeeMJnPRs2bFBpaak++eQTZWVlaffu3erdu7f+/ve/+9RFRUXV297YsWPt+RERESovL9eHH37os9wrr7yirl27NtJf4eL49fLWP//5T5WXl/vs7IkTJzRjxgwtXLhQ33zzjVwul2pqalRRUeHT21NeXq4BAwZIklwul7777rt66//+++8VFxd3xm2HhIQoJCTEn7sDAECLERISIpfLVW/6sWPHdN999+m2225Tbm6uPT0hIUH9+vWrd+dzTEyMvZ7LLrtMo0eP1pAhQzRhwgR99dVXCggIkHSy8+FM2zslMDBQY8aM0SuvvKLk5GRJ0sGDB/XBBx/o4YcftgNZS+LX0JOenm4PTj4lJSVF6enpuu+++yRJffv2VVBQkPLy8pSWliZJKi0tVXFxsebNmydJSk5Olsfj0datW3XjjTdKkj766CN5PB47GJmqqa47c80ZAFqH999/Xz/88INmzZp1xvkdOnQ45/Lt2rXTQw89pP/4j/9QYWGh/bt7PiZMmKBbbrlFzz77rMLDw5WTk6MRI0actYOiuTU49FRVVenLL7+0v5eUlKioqEjR0dHq2rWrYmJifOqDgoLkcrnswcdOp1MTJkzQjBkzFBMTo+joaM2cOVO9evWyA1OPHj00YsQITZw4US+88IIk6f7779eoUaO4cwsAYKS3335b7du395k2e/Zs+yrHVVdddcHrPrXsN998Y4cej8fjs7327dvXG1d77bXX6vLLL9ff/vY3paenKycnR/Pnz9fXX399wW1pTA0OPdu3b9fgwYPt79OnT5ck3XvvvcrJyTmvdSxYsECBgYFKS0tTdXW1hgwZopycHLtLTZJeffVVTZs2zb7LKzU19WefDQQAQFs1ePBgLVmyxGdadHS0XnzxxYtet2VZknzH00ZGRmrHjh3297M9UHD8+PFatmyZunbtqqqqKt12220t9ve6waFn0KBB9h/nfHzzzTf1poWGhmrRokVatGjRWZeLjo7WypUrG9o8AADapIiICF1xxRX1pl955cnhCJ9//rk9tqahdu/eLenkOKBT2rVrd8btnW7s2LGaNWuW5syZo9/+9rcKDGy5b7jiWdwAALRiw4cPV2xsrD0u9nQ/9wqnuro6Pffcc0pISFCfPn0avP3o6GilpqYqPz9f48ePb/DyTanlxjEAAGDzer31xtQEBgYqNjZWL730ku666y6lpqZq2rRpuuKKK/TDDz/ojTfe0P79+7V69Wp7mR9//FFlZWU6duyYiouLtXDhQm3dulXvvPOOzzCThsjJydHzzz9fb1xvS0PoAQAYrzXcsbpu3Tp17tzZZ1r37t31+eef6/bbb1dBQYGys7M1ZswYVVZWKj4+Xr/61a/0xz/+0WeZUzcNhYeHq1u3bho8eLCWLl16XpeyziYsLExhYWEXvHxTcVgNGaDTilRWVsrpdMrj8SgqKqq5m+M33LIOABfmp59+UklJiRISEhQaGtrczTHauY5FY/5+M6YHAAAYgdADAACMQOgBAABGIPQAAAAjEHoAAIARCD0AAMAIhB4AAGAEQg8AADACoQcAABiB11AAALAxu+m2NTizwYuMGzdOhw8f1po1a1ReXq4nnnhC7733nr777jt17NhRvXv31pw5c+y3rF966aXKyMhQRkaG/X3fvn167bXXdM899/isu2fPnvrss8+0bNkyjRs37pzLSydfOXHZZZdp6tSpmjRp0oX9DZoJPT0AALQiv/71r/XJJ59o+fLl2rt3r9auXatBgwbpX//61zmXi4+P17Jly3ymbdmyRWVlZYqIiPjZ7T799NMqLS3Vp59+qjvuuEMPPPCAXn/99Yval6ZG6AEAoJU4fPiwNm/erLlz52rw4MHq1q2bbrzxRmVmZmrkyJHnXHbs2LHKz8/XgQMH7GmvvPKKxo4dq8DAn7/wExkZKZfLpSuuuEJ//OMflZiYqDVr1lzsLjUpQg8AAK1E+/bt1b59e61Zs0Zer7dBy8bFxSklJUXLly+XJB07dkyvv/66xo8ff0FtCQ0NVW1t7QUt21wIPQAAtBKBgYHKycnR8uXL1aFDB91000167LHH9Omnn57X8uPHj1dOTo4sy9Lf/vY3XX755br22msb1Ibjx48rJydHO3fu1JAhQy5gL5oPoQcAgFbk17/+tQ4dOqS1a9cqJSVFH3zwga677jrl5OT87LIjR45UVVWVNm3apFdeeaVBvTyzZ89W+/btFRYWpgcffFCPPPIIA5kBAEDjCg0N1bBhw/SHP/xBBQUFGjdunJ588smfXS4wMFDp6el68skn9dFHH2ns2LHnvc1HHnlERUVF2rdvn6qqqjRv3jy1a9e6YkTrai0AAKjn6quv1tGjR8+rdvz48crPz9ftt9+ujh07nvc2YmNjdcUVV8jtdsvhcFxoU5sVz+kBAKCV+PHHH3XXXXdp/PjxuuaaaxQZGant27dr3rx5uv32289rHT169NAPP/yg8PDwRm5ty0PoAQCglWjfvr369eunBQsW6KuvvlJtba3i4+M1ceJEPfbYY+e9npiYmEZsZcvlsCzLau5GNIbKyko5nU55PB5FRUU1d3P8ZkHe3ibZzsPDrmyS7QBAU/npp59UUlKihIQEhYaGNndzjHauY9GYv9+M6QEAAEYg9AAAACMQegAAgBEIPQAAwAiEHgAAYARCDwDAKG30puVWpbmOAaEHAGCEoKAgSSffLo7mdeoYnDomTYWHEwIAjBAQEKAOHTqovLxckhQeHt5qX6fQWlmWpWPHjqm8vFwdOnRQQEBAk26f0AMAMIbL5ZIkO/igeXTo0ME+Fk2J0AMAMIbD4VDnzp3VqVMn1dbWNndzjBQUFNTkPTynEHoAAMYJCAhoth9eNB8GMgMAACMQegAAgBEIPQAAwAiEHgAAYIQGh55NmzZp9OjRcrvdcjgcWrNmjT2vtrZWs2fPVq9evRQRESG3263f/va3OnTokM86vF6vpk6dqtjYWEVERCg1NVUHDx70qamoqFB6erqcTqecTqfS09N1+PDhC9pJAACABoeeo0ePqnfv3lq8eHG9eceOHdOOHTv0xBNPaMeOHXrrrbe0d+9epaam+tRlZGQoNzdXq1ev1ubNm1VVVaVRo0bpxIkTds2YMWNUVFSkdevWad26dSoqKlJ6evoF7CIAAIDksC7iBRgOh0O5ubm64447zlqzbds23Xjjjdq3b5+6du0qj8ejSy65RCtWrNDdd98tSTp06JDi4+P17rvvKiUlRbt379bVV1+tLVu2qF+/fpKkLVu2KDk5WZ9//rm6d+/+s22rrKyU0+mUx+NRVFTUhe5ii7Mgb2+TbOfhYVc2yXYAAPh3jfn73ehjejwejxwOhzp06CBJKiwsVG1trYYPH27XuN1uJSUlqaCgQJL04Ycfyul02oFHkvr37y+n02nXAAAANESjPpzwp59+0qOPPqoxY8bYaa2srEzBwcHq2LGjT21cXJzKysrsmk6dOtVbX6dOneya03m9Xnm9Xvt7ZWWlv3YDAAC0AY3W01NbW6t77rlHdXV1ev7553+23rIsnxe/neklcKfX/Lvs7Gx70LPT6VR8fPyFNx4AALQ5jRJ6amtrlZaWppKSEuXl5flck3O5XKqpqVFFRYXPMuXl5YqLi7Nrvvvuu3rr/f777+2a02VmZsrj8difAwcO+HGPAABAa+f30HMq8HzxxRfasGGDYmJifOb37dtXQUFBysvLs6eVlpaquLhYAwYMkCQlJyfL4/Fo69atds1HH30kj8dj15wuJCREUVFRPh8AAIBTGjymp6qqSl9++aX9vaSkREVFRYqOjpbb7dZ//ud/aseOHXr77bd14sQJewxOdHS0goOD5XQ6NWHCBM2YMUMxMTGKjo7WzJkz1atXLw0dOlSS1KNHD40YMUITJ07UCy+8IEm6//77NWrUqPO6cwsAAOB0DQ4927dv1+DBg+3v06dPlyTde++9mjNnjtauXStJuvbaa32W27hxowYNGiRJWrBggQIDA5WWlqbq6moNGTJEOTk5Pm+8ffXVVzVt2jT7Lq/U1NQzPhsIAADgfFzUc3paMp7Tc3F4Tg8AoDm06uf0AAAAtASEHgAAYARCDwAAMAKhBwAAGIHQAwAAjEDoAQAARiD0AAAAIxB6AACAEQg9AADACIQeAABgBEIPAAAwAqEHAAAYgdADAACMQOgBAABGIPQAAAAjEHoAAIARCD0AAMAIhB4AAGAEQg8AADACoQcAABiB0AMAAIxA6AEAAEYg9AAAACMQegAAgBEIPQAAwAiEHgAAYARCDwAAMAKhBwAAGCGwuRuAlmlB3t4m2c7Dw65sku0AAEBPDwAAMAKhBwAAGIHQAwAAjEDoAQAARiD0AAAAIxB6AACAEQg9AADACIQeAABgBEIPAAAwAqEHAAAYgdADAACM0ODQs2nTJo0ePVput1sOh0Nr1qzxmW9ZlubMmSO3262wsDANGjRIu3bt8qnxer2aOnWqYmNjFRERodTUVB08eNCnpqKiQunp6XI6nXI6nUpPT9fhw4cbvIMAAADSBYSeo0ePqnfv3lq8ePEZ58+bN0/z58/X4sWLtW3bNrlcLg0bNkxHjhyxazIyMpSbm6vVq1dr8+bNqqqq0qhRo3TixAm7ZsyYMSoqKtK6deu0bt06FRUVKT09/QJ2EQAAQHJYlmVd8MIOh3Jzc3XHHXdIOtnL43a7lZGRodmzZ0s62asTFxenuXPnatKkSfJ4PLrkkku0YsUK3X333ZKkQ4cOKT4+Xu+++65SUlK0e/duXX311dqyZYv69esnSdqyZYuSk5P1+eefq3v37j/btsrKSjmdTnk8HkVFRV3oLrY4TfX286bCW9YBAP+uMX+//Tqmp6SkRGVlZRo+fLg9LSQkRAMHDlRBQYEkqbCwULW1tT41brdbSUlJds2HH34op9NpBx5J6t+/v5xOp11zOq/Xq8rKSp8PAADAKX4NPWVlZZKkuLg4n+lxcXH2vLKyMgUHB6tjx47nrOnUqVO99Xfq1MmuOV12drY9/sfpdCo+Pv6i9wcAALQdjXL3lsPh8PluWVa9aac7veZM9edaT2Zmpjwej/05cODABbQcAAC0VX4NPS6XS5Lq9caUl5fbvT8ul0s1NTWqqKg4Z813331Xb/3ff/99vV6kU0JCQhQVFeXzAQAAOMWvoSchIUEul0t5eXn2tJqaGuXn52vAgAGSpL59+yooKMinprS0VMXFxXZNcnKyPB6Ptm7datd89NFH8ng8dg0AAEBDBDZ0gaqqKn355Zf295KSEhUVFSk6Olpdu3ZVRkaGsrKylJiYqMTERGVlZSk8PFxjxoyRJDmdTk2YMEEzZsxQTEyMoqOjNXPmTPXq1UtDhw6VJPXo0UMjRozQxIkT9cILL0iS7r//fo0aNeq87twCAAA4XYNDz/bt2zV48GD7+/Tp0yVJ9957r3JycjRr1ixVV1dr8uTJqqioUL9+/bR+/XpFRkbayyxYsECBgYFKS0tTdXW1hgwZopycHAUEBNg1r776qqZNm2bf5ZWamnrWZwMBAAD8nIt6Tk9LxnN6Wgee0wMA+Het5jk9AAAALRWhBwAAGIHQAwAAjEDoAQAARiD0AAAAIxB6AACAEQg9AADACIQeAABgBEIPAAAwAqEHAAAYgdADAACMQOgBAABGIPQAAAAjEHoAAIARCD0AAMAIhB4AAGAEQg8AADACoQcAABiB0AMAAIxA6AEAAEYg9AAAACMQegAAgBEIPQAAwAiEHgAAYARCDwAAMAKhBwAAGIHQAwAAjEDoAQAARiD0AAAAIxB6AACAEQg9AADACIQeAABgBEIPAAAwAqEHAAAYgdADAACMQOgBAABGIPQAAAAjEHoAAIARCD0AAMAIhB4AAGAEv4ee48eP6/e//70SEhIUFhamyy67TE8//bTq6ursGsuyNGfOHLndboWFhWnQoEHatWuXz3q8Xq+mTp2q2NhYRUREKDU1VQcPHvR3cwEAgCH8Hnrmzp2rv/71r1q8eLF2796tefPm6X/+53+0aNEiu2bevHmaP3++Fi9erG3btsnlcmnYsGE6cuSIXZORkaHc3FytXr1amzdvVlVVlUaNGqUTJ074u8kAAMAAgf5e4Ycffqjbb79dI0eOlCRdeumleu2117R9+3ZJJ3t5Fi5cqMcff1x33nmnJGn58uWKi4vTqlWrNGnSJHk8Hr388stasWKFhg4dKklauXKl4uPjtWHDBqWkpPi72QAAoI3ze0/PzTffrL///e/au3evJOmTTz7R5s2bddttt0mSSkpKVFZWpuHDh9vLhISEaODAgSooKJAkFRYWqra21qfG7XYrKSnJrjmd1+tVZWWlzwcAAOAUv/f0zJ49Wx6PR1dddZUCAgJ04sQJ/elPf9JvfvMbSVJZWZkkKS4uzme5uLg47du3z64JDg5Wx44d69WcWv502dnZeuqpp/y9OwAAoI3we0/P66+/rpUrV2rVqlXasWOHli9frr/85S9avny5T53D4fD5bllWvWmnO1dNZmamPB6P/Tlw4MDF7QgAAGhT/N7T88gjj+jRRx/VPffcI0nq1auX9u3bp+zsbN17771yuVySTvbmdO7c2V6uvLzc7v1xuVyqqalRRUWFT29PeXm5BgwYcMbthoSEKCQkxN+7AwAA2gi/9/QcO3ZM7dr5rjYgIMC+ZT0hIUEul0t5eXn2/JqaGuXn59uBpm/fvgoKCvKpKS0tVXFx8VlDDwAAwLn4vadn9OjR+tOf/qSuXbuqZ8+e+vjjjzV//nyNHz9e0snLWhkZGcrKylJiYqISExOVlZWl8PBwjRkzRpLkdDo1YcIEzZgxQzExMYqOjtbMmTPVq1cv+24uAACAhvB76Fm0aJGeeOIJTZ48WeXl5XK73Zo0aZL+8Ic/2DWzZs1SdXW1Jk+erIqKCvXr10/r169XZGSkXbNgwQIFBgYqLS1N1dXVGjJkiHJychQQEODvJgMAAAM4LMuymrsRjaGyslJOp1Mej0dRUVHN3Ry/WZC3t7mb4FcPD7uyuZsAAGhBGvP3m3dvAQAAIxB6AACAEQg9AADACIQeAABgBEIPAAAwAqEHAAAYgdADAACMQOgBAABGIPQAAAAjEHoAAIARCD0AAMAIhB4AAGAEQg8AADACoQcAABiB0AMAAIwQ2NwNAFqcjdkXv47BmRe/DgCAX9HTAwAAjEDoAQAARiD0AAAAIxB6AACAEQg9AADACIQeAABgBG5ZR7NakLe3ybb18LArm2xbAICWh54eAABgBEIPAAAwAqEHAAAYgdADAACMQOgBAABGIPQAAAAjEHoAAIARCD0AAMAIPJwQaAwbsy9+HYMzL34dAAAbPT0AAMAIhB4AAGAEQg8AADACoQcAABiBgcxoU/rvX3r2mRtjmq4hAIAWh54eAABgBEIPAAAwAqEHAAAYoVFCz7fffqv/+q//UkxMjMLDw3XttdeqsLDQnm9ZlubMmSO3262wsDANGjRIu3bt8lmH1+vV1KlTFRsbq4iICKWmpurgwYON0VwAAGAAv4eeiooK3XTTTQoKCtJ7772nzz77TM8884w6dOhg18ybN0/z58/X4sWLtW3bNrlcLg0bNkxHjhyxazIyMpSbm6vVq1dr8+bNqqqq0qhRo3TixAl/NxkAABjA73dvzZ07V/Hx8Vq2bJk97dJLL7X/27IsLVy4UI8//rjuvPNOSdLy5csVFxenVatWadKkSfJ4PHr55Ze1YsUKDR06VJK0cuVKxcfHa8OGDUpJSfF3swEAQBvn956etWvX6vrrr9ddd92lTp06qU+fPnrxxRft+SUlJSorK9Pw4cPtaSEhIRo4cKAKCgokSYWFhaqtrfWpcbvdSkpKsmtO5/V6VVlZ6fMBAAA4xe+h5+uvv9aSJUuUmJio999/Xw888ICmTZum//3f/5UklZWVSZLi4uJ8louLi7PnlZWVKTg4WB07djxrzemys7PldDrtT3x8vL93DQAAtGJ+Dz11dXW67rrrlJWVpT59+mjSpEmaOHGilixZ4lPncDh8vluWVW/a6c5Vk5mZKY/HY38OHDhwcTsCAADaFL+Hns6dO+vqq6/2mdajRw/t379fkuRyuSSpXo9NeXm53fvjcrlUU1OjioqKs9acLiQkRFFRUT4fAACAU/w+kPmmm27Snj17fKbt3btX3bp1kyQlJCTI5XIpLy9Pffr0kSTV1NQoPz9fc+fOlST17dtXQUFBysvLU1pamiSptLRUxcXFmjdvnr+bDABAo1uQt7dJtvPwsCubZDutkd9Dz8MPP6wBAwYoKytLaWlp2rp1q5YuXaqlS0++E8nhcCgjI0NZWVlKTExUYmKisrKyFB4erjFjxkiSnE6nJkyYoBkzZigmJkbR0dGaOXOmevXqZd/NBQAA0BB+Dz033HCDcnNzlZmZqaeffloJCQlauHChxo4da9fMmjVL1dXVmjx5sioqKtSvXz+tX79ekZGRds2CBQsUGBiotLQ0VVdXa8iQIcrJyVFAQIC/mwwAAAzgsCzLau5GNIbKyko5nU55PJ42Nb6nqbpHW6tzvWU9+bJW9pb1wZnN3QIAfsTlrfPTmL/fvHsLAAAYgdADAACMQOgBAABGIPQAAAAjEHoAAIARCD0AAMAIhB4AAGAEQg8AADACoQcAABiB0AMAAIxA6AEAAEYg9AAAACMQegAAgBEIPQAAwAiEHgAAYARCDwAAMAKhBwAAGIHQAwAAjEDoAQAARiD0AAAAIxB6AACAEQg9AADACIQeAABgBEIPAAAwAqEHAAAYIbC5GwDgLDZmX/w6Bmde/DoAoI2gpwcAABiB0AMAAIzA5S0AANqQBXl7m2xbDw+7ssm25Q/09AAAACPQ0wMAMFpT9oygeRF6AADNj7sV0QS4vAUAAIxA6AEAAEYg9AAAACMQegAAgBEIPQAAwAjcvQVjfPj1j02yneTLYppkOwCAhqGnBwAAGIHQAwAAjNDooSc7O1sOh0MZGRn2NMuyNGfOHLndboWFhWnQoEHatWuXz3Jer1dTp05VbGysIiIilJqaqoMHDzZ2cwEAQBvVqKFn27ZtWrp0qa655hqf6fPmzdP8+fO1ePFibdu2TS6XS8OGDdORI0fsmoyMDOXm5mr16tXavHmzqqqqNGrUKJ04caIxmwwAANqoRgs9VVVVGjt2rF588UV17NjRnm5ZlhYuXKjHH39cd955p5KSkrR8+XIdO3ZMq1atkiR5PB69/PLLeuaZZzR06FD16dNHK1eu1M6dO7Vhw4bGajIAAGjDGi30PPjggxo5cqSGDh3qM72kpERlZWUaPny4PS0kJEQDBw5UQUGBJKmwsFC1tbU+NW63W0lJSXbN6bxeryorK30+AAAApzTKLeurV6/Wjh07tG3btnrzysrKJElxcXE+0+Pi4rRv3z67Jjg42KeH6FTNqeVPl52draeeesofzQcAAG2Q33t6Dhw4oIceekgrV65UaGjoWescDofPd8uy6k073blqMjMz5fF47M+BAwca3ngAANBm+T30FBYWqry8XH379lVgYKACAwOVn5+v5557ToGBgXYPz+k9NuXl5fY8l8ulmpoaVVRUnLXmdCEhIYqKivL5AAAAnOL3y1tDhgzRzp07fabdd999uuqqqzR79mxddtllcrlcysvLU58+fSRJNTU1ys/P19y5cyVJffv2VVBQkPLy8pSWliZJKi0tVXFxsebNm+fvJqOF6L9/aXM3AQDQhvk99ERGRiopKclnWkREhGJiYuzpGRkZysrKUmJiohITE5WVlaXw8HCNGTNGkuR0OjVhwgTNmDFDMTExio6O1syZM9WrV696A6MBAG3Phbw2ZsvxvY3QErQlzfLurVmzZqm6ulqTJ09WRUWF+vXrp/Xr1ysyMtKuWbBggQIDA5WWlqbq6moNGTJEOTk5CggIaI4mA63TxuyLX8fgzItfBwC0AE0Sej744AOf7w6HQ3PmzNGcOXPOukxoaKgWLVqkRYsWNW7jAACAEXjLOgCYjN5AGITQAwC4OP4ITn7gj5shtnS93w8tQUvFW9YBAIARCD0AAMAIhB4AAGAEQg8AADACA5kBnBt39wBoI+jpAQAARiD0AAAAIxB6AACAERjTA/jZhbwo8UIkXxbTJNsBgLaC0AOg8TEYGkALQOgB0DoQnABcJMb0AAAAI9DTAwA4b001Zg1oDPT0AAAAIxB6AACAEbi8BbRSTXmZgdvjWyh/DO4GDELoAWAOf4UE7gIDWiUubwEAACMQegAAgBG4vAUAzYHxOECTI/QA+Fm8TwxAW8DlLQAAYARCDwAAMAKhBwAAGIExPQDQUAxCBlolenoAAIARCD0AAMAIXN4C0GJwa/yFa8p3sQGtFT09AADACIQeAABgBEIPAAAwAqEHAAAYgYHMANCIGGDcuvTfv/Si17Gl6/1+aAkaA6EHgHEIIoCZuLwFAACMQOgBAABGIPQAAAAj+D30ZGdn64YbblBkZKQ6deqkO+64Q3v27PGpsSxLc+bMkdvtVlhYmAYNGqRdu3b51Hi9Xk2dOlWxsbGKiIhQamqqDh486O/mAgAAQ/g99OTn5+vBBx/Uli1blJeXp+PHj2v48OE6evSoXTNv3jzNnz9fixcv1rZt2+RyuTRs2DAdOXLErsnIyFBubq5Wr16tzZs3q6qqSqNGjdKJEyf83WQAAGAAv9+9tW7dOp/vy5YtU6dOnVRYWKhbbrlFlmVp4cKFevzxx3XnnXdKkpYvX664uDitWrVKkyZNksfj0csvv6wVK1Zo6NChkqSVK1cqPj5eGzZsUEpKir+bDQAA2rhGH9Pj8XgkSdHR0ZKkkpISlZWVafjw4XZNSEiIBg4cqIKCAklSYWGhamtrfWrcbreSkpLsmtN5vV5VVlb6fAAAAE5p1NBjWZamT5+um2++WUlJSZKksrIySVJcXJxPbVxcnD2vrKxMwcHB6tix41lrTpednS2n02l/4uPj/b07AACgFWvU0DNlyhR9+umneu211+rNczgcPt8ty6o37XTnqsnMzJTH47E/Bw4cuPCGAwCANqfRnsg8depUrV27Vps2bVKXLl3s6S6XS9LJ3pzOnTvb08vLy+3eH5fLpZqaGlVUVPj09pSXl2vAgAFn3F5ISIhCQkIaY1dwHvzx6HYAABqT33t6LMvSlClT9NZbb+kf//iHEhISfOYnJCTI5XIpLy/PnlZTU6P8/Hw70PTt21dBQUE+NaWlpSouLj5r6AEAADgXv/f0PPjgg1q1apX+7//+T5GRkfYYHKfTqbCwMDkcDmVkZCgrK0uJiYlKTExUVlaWwsPDNWbMGLt2woQJmjFjhmJiYhQdHa2ZM2eqV69e9t1cAAAADeH30LNkyRJJ0qBBg3ymL1u2TOPGjZMkzZo1S9XV1Zo8ebIqKirUr18/rV+/XpGRkXb9ggULFBgYqLS0NFVXV2vIkCHKyclRQECAv5sMAAAM4LAsy2ruRjSGyspKOZ1OeTweRUVFNXdz/GZB3t7mbsIZMaYHAE7a0vX+5m5Ck3l42JV+X2dj/n7z7i0AAGAEQg8AADACoQcAABiB0AMAAIxA6AEAAEYg9AAAACMQegAAgBEIPQAAwAiEHgAAYARCDwAAMAKhBwAAGIHQAwAAjEDoAQAARiD0AAAAIxB6AACAEQKbuwEAALQl/fcvveh1bOl6vx9agtPR0wMAAIxA6AEAAEYg9AAAACMQegAAgBEIPQAAwAiEHgAAYARCDwAAMALP6YFfnikBAEBLR08PAAAwAqEHAAAYgdADAACMQOgBAABGIPQAAAAjEHoAAIARCD0AAMAIhB4AAGAEQg8AADACoQcAABiB11C0crxCAgCA80NPDwAAMAKhBwAAGIHQAwAAjEDoAQAARmAgc3PamN3gRfrv/7ERGgIAaEn8cZPKlq73+6ElbUuL7+l5/vnnlZCQoNDQUPXt21f//Oc/m7tJAACgFWrRoef1119XRkaGHn/8cX388cf65S9/qVtvvVX79+9v7qYBAIBWpkWHnvnz52vChAn63e9+px49emjhwoWKj4/XkiVLmrtpAACglWmxY3pqampUWFioRx991Gf68OHDVVBQUK/e6/XK6/Xa3z0ejySpsrKycRq46ZnGWe/POFrt/fkiAIDxfjpa1ejbaIzf2FPrtCzL7+tusaHnhx9+0IkTJxQXF+czPS4uTmVlZfXqs7Oz9dRTT9WbHh8f32htBACg5Vrc6Ft4rBHXfeTIETmdTr+us8WGnlMcDofPd8uy6k2TpMzMTE2fPt3+XldXp3/961+KiYk5Y71pKisrFR8frwMHDigqKqq5m4Nz4Fi1Hhyr1oNj1TqcOk6fffaZ3G6339ffYkNPbGysAgIC6vXqlJeX1+v9kaSQkBCFhIT4TOvQoUNjNrFVioqK4oRvJThWrQfHqvXgWLUOv/jFL9Sunf+HHbfYgczBwcHq27ev8vLyfKbn5eVpwIABzdQqAADQWrXYnh5Jmj59utLT03X99dcrOTlZS5cu1f79+/XAAw80d9MAAEAr06JDz913360ff/xRTz/9tEpLS5WUlKR3331X3bp1a+6mtTohISF68skn610CRMvDsWo9OFatB8eqdWjs4+SwGuOeMAAAgBamxY7pAQAA8CdCDwAAMAKhBwAAGIHQAwAAjEDoacWys7N1ww03KDIyUp06ddIdd9yhPXv2+NRUVVVpypQp6tKli8LCwtSjR496L2wdNGiQHA6Hz+eee+5pyl1p087nOH333XcaN26c3G63wsPDNWLECH3xxRc+NV6vV1OnTlVsbKwiIiKUmpqqgwcPNuWutHn+OlacU41vyZIluuaaa+yHDSYnJ+u9996z51uWpTlz5sjtdissLEyDBg3Srl27fNbBOdU0/HGs/HZOWWi1UlJSrGXLllnFxcVWUVGRNXLkSKtr165WVVWVXfO73/3Ouvzyy62NGzdaJSUl1gsvvGAFBARYa9assWsGDhxoTZw40SotLbU/hw8fbo5dapN+7jjV1dVZ/fv3t375y19aW7dutT7//HPr/vvvr3csH3jgAesXv/iFlZeXZ+3YscMaPHiw1bt3b+v48ePNtWttjr+OFedU41u7dq31zjvvWHv27LH27NljPfbYY1ZQUJBVXFxsWZZl/fnPf7YiIyOtN99809q5c6d19913W507d7YqKyvtdXBONQ1/HCt/nVOEnjakvLzckmTl5+fb03r27Gk9/fTTPnXXXXed9fvf/97+PnDgQOuhhx5qqmYa7/TjtGfPHkuS/Q+AZVnW8ePHrejoaOvFF1+0LMuyDh8+bAUFBVmrV6+2a7799lurXbt21rp165p2BwxyIcfKsjinmkvHjh2tl156yaqrq7NcLpf15z//2Z73008/WU6n0/rrX/9qWRbnVHNryLGyLP+dU1zeakM8Ho8kKTo62p528803a+3atfr2229lWZY2btyovXv3KiUlxWfZV199VbGxserZs6dmzpypI0eONGnbTXL6cfJ6vZKk0NBQuyYgIEDBwcHavHmzJKmwsFC1tbUaPny4XeN2u5WUlKSCgoKmarpxLuRYncI51XROnDih1atX6+jRo0pOTlZJSYnKysp8zpeQkBANHDjQPl84p5rHhRyrU/xxTrXoJzLj/FmWpenTp+vmm29WUlKSPf25557TxIkT1aVLFwUGBqpdu3Z66aWXdPPNN9s1Y8eOVUJCglwul4qLi5WZmalPPvmk3nvPcPHOdJyuuuoqdevWTZmZmXrhhRcUERGh+fPnq6ysTKWlpZKksrIyBQcHq2PHjj7ri4uLq/dSXvjHhR4riXOqqezcuVPJycn66aef1L59e+Xm5urqq6+2fyxPfzl1XFyc9u3bJ4lzqqldzLGS/HdOEXraiClTpujTTz+t93+bzz33nLZs2aK1a9eqW7du2rRpkyZPnqzOnTtr6NChkqSJEyfa9UlJSUpMTNT111+vHTt26LrrrmvS/WjrznScgoKC9Oabb2rChAmKjo5WQECAhg4dqltvvfVn12dZlhwOR2M22VgXc6w4p5pG9+7dVVRUpMOHD+vNN9/Uvffeq/z8fHv+6efG+ZwvnFON42KPlb/OKS5vtQFTp07V2rVrtXHjRnXp0sWeXl1drccee0zz58/X6NGjdc0112jKlCm6++679Ze//OWs67vuuusUFBRU744UXJyzHSdJ6tu3r/0PQmlpqdatW6cff/xRCQkJkiSXy6WamhpVVFT4LFdeXl7v/5Bw8S7mWJ0J51TjCA4O1hVXXKHrr79e2dnZ6t27t5599lm5XC5Jqtdj8+/nC+dU07qYY3UmF3pOEXpaMcuyNGXKFL311lv6xz/+Ue8f3draWtXW1qpdO9/DHBAQoLq6urOud9euXaqtrVXnzp0bpd2m+bnj9O+cTqcuueQSffHFF9q+fbtuv/12SSd/aIOCgny6cktLS1VcXKwBAwY0+j6Ywh/H6kw4p5qGZVnyer32ZZB/P19qamqUn59vny+cU82rIcfqTC74nLroodBoNv/93/9tOZ1O64MPPvC5je/YsWN2zcCBA62ePXtaGzdutL7++mtr2bJlVmhoqPX8889blmVZX375pfXUU09Z27Zts0pKSqx33nnHuuqqq6w+ffpw26afnM9xeuONN6yNGzdaX331lbVmzRqrW7du1p133umzngceeMDq0qWLtWHDBmvHjh3Wr371K26v9TN/HCvOqaaRmZlpbdq0ySopKbE+/fRT67HHHrPatWtnrV+/3rKsk7dBO51O66233rJ27txp/eY3vznjLeucU43vYo+VP88pQk8rJumMn2XLltk1paWl1rhx4yy3222FhoZa3bt3t5555hmrrq7OsizL2r9/v3XLLbdY0dHRVnBwsHX55Zdb06ZNs3788cdm2qu253yO07PPPmt16dLFCgoKsrp27Wr9/ve/t7xer896qqurrSlTpljR0dFWWFiYNWrUKGv//v1NvDdtmz+OFedU0xg/frzVrVs3Kzg42LrkkkusIUOG2D+ilnXymUpPPvmk5XK5rJCQEOuWW26xdu7c6bMOzqmmcbHHyp/nlMOyLKthfUMAAACtD2N6AACAEQg9AADACIQeAABgBEIPAAAwAqEHAAAYgdADAACMQOgBAABGIPQAAAAjEHoAAIARCD0AAMAIhB4AAGAEQg8AADDC/wNs9QtuE42LNAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.hist(tas_ECDFM[:, 1, 1], bins=\"auto\", alpha = 0.5, label = \"ECDFM\")\n", "plt.hist(tas_ISIMIP[:, 1, 1], bins=\"auto\", alpha = 0.5, label = \"ISIMIP\")\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "b7f98cdd-d056-4cd4-b756-a29d1002d32a", "metadata": {}, "source": [ "We see that the data distributions produced by ECDFM and ISIMIP do differ. For info on how to evaluate this have a look at our notebook on evaluation" ] }, { "cell_type": "markdown", "id": "9f8c657f-1c05-40a6-95cb-2a4497f70bd9", "metadata": {}, "source": [ "## 2. Going deeper: initialising and applying debiasers" ] }, { "cell_type": "markdown", "id": "51825b29", "metadata": {}, "source": [ "Every debiaser that is part of the package consists of a set of parameters controlling the behavior as well as a `from_variable` and `apply`-method." ] }, { "cell_type": "markdown", "id": "e68dd34f-bb4c-4819-87d2-db966b8510b7", "metadata": {}, "source": [ "### 2.1. Initialisation" ] }, { "cell_type": "markdown", "id": "1d88c145-1b49-485f-8661-dfdb561fed56", "metadata": {}, "source": [ "Let's have a look at the tas-ECDFM debiaser defined above using the `from_variable` method." ] }, { "cell_type": "code", "execution_count": 8, "id": "8ee38a89-4459-47a6-ad31-f797be1b98db", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "ECDFM(variable='Daily mean near-surface air temperature', reasonable_physical_range=[100, 400], distribution=, cdf_threshold=1e-10)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "tas_debiaser_ECDFM" ] }, { "cell_type": "markdown", "id": "0d08991e-0107-4c33-8580-6d4ec34fffaf", "metadata": {}, "source": [ "We see that two parameters are set: a variable referring to the variable we are debiasing, as well as a distribution. These parameters fully determine the setting of this particular debiaser, equidistant-CDF-matching.\n", "\n", "Using the `from_variable` method, a debiaser is initialized using the default parameters associated with this variable for this particular debiaser. The same outcome can be achieved by directly setting the required parameters, without using the `from_variable` method:" ] }, { "cell_type": "code", "execution_count": 9, "id": "c7df1d2d-3af7-4297-838b-df53939e586d", "metadata": {}, "outputs": [], "source": [ "import scipy.stats\n", "tas_debiaser_ECDFM_v2 = ECDFM(distribution=scipy.stats.beta)" ] }, { "cell_type": "markdown", "id": "55464614", "metadata": {}, "source": [ "And we can see that these debiasers are absolutely identical:" ] }, { "cell_type": "code", "execution_count": 10, "id": "41a6b923", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "tas_debiaser_ECDFM == tas_debiaser_ECDFM_v2" ] }, { "cell_type": "markdown", "id": "6a01b50b", "metadata": {}, "source": [ "**Learning**: if we set all the parameters needed manually to the default parameters of a specific variable, the debiaser will be equivalent to a debiaser initialized using the `from_variable` method.\n", "\n", "Have a look at the documentation of each bias adjustment method to find out which parameters need to be set for that specific debiaser.\n", "\n", "The table below gives an overview of which variables currently have default setting for which debiasers - 'experimental default settings' are marked with brackets around the x." ] }, { "cell_type": "markdown", "id": "46c0d20a-d3a2-496b-a13d-f1816c8bada9", "metadata": {}, "source": [ "| Variable | `LinearScaling` | `DeltaChange` | `QuantileMapping` | `ScaledDistributionMapping` | `CDFt` | `ECDFM` | `QuantileDeltaMapping` | `ISIMIP` |\n", "|-----------|:-------------------------:|:-----------------------:|:---------------------------:|:-------------------------------------:|:-----------------:|:-----------------:|:--------------------------------:|:------------------:|\n", "| hurs | (x) | (x) | (x) | | (x) | (x) | (x) | x |\n", "| pr | x | x | x | x | x | x | x | x |\n", "| prsnratio | | | | | | | | x |\n", "| psl | (x) | (x) | (x) | | (x) | (x) | (x) | x |\n", "| rlds | (x) | (x) | (x) | | (x) | (x) | (x) | x |\n", "| rsds | (x) | (x) | | | (x) | | | x |\n", "| sfcWind | (x) | (x) | (x) | | (x) | (x) | (x) | x |\n", "| tas | x | x | x | x | x | x | x | x |\n", "| tasmin | x | x | (x) | (x) | x | (x) | (x) | |\n", "| tasmax | x | x | (x) | (x) | x | (x) | (x) | |\n", "| tasrange | | | | | (x) | | | x |\n", "| tasskew | | | | | (x) | | | x |\n" ] }, { "cell_type": "markdown", "id": "c1e9127a-f408-4b3f-891f-113591d73fe6", "metadata": {}, "source": [ "If we try to initialize a debiaser that does not exist for a certain variable using the `from_variable` method, we are basically asking for error message:" ] }, { "cell_type": "code", "execution_count": 11, "id": "40c0f5b1-e52c-4eaa-9398-9d6db556b705", "metadata": {}, "outputs": [ { "ename": "ValueError", "evalue": "Unfortunately currently no default settings exist for the variable tasmin in the debiaser ISIMIP. You can set the required class parameters manually by using the class constructor.", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m/tmp/ipykernel_49752/2091483589.py\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mdebiaser3\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mISIMIP\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfrom_variable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"tasmin\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[0;32m~/anaconda3/lib/python3.9/site-packages/ibicus/debias/_isimip.py\u001b[0m in \u001b[0;36mfrom_variable\u001b[0;34m(cls, variable, **kwargs)\u001b[0m\n\u001b[1;32m 341\u001b[0m \u001b[0;34m@\u001b[0m\u001b[0mclassmethod\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 342\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mfrom_variable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcls\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvariable\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mUnion\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mstr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mVariable\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 343\u001b[0;31m return super()._from_variable(\n\u001b[0m\u001b[1;32m 344\u001b[0m \u001b[0mcls\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 345\u001b[0m \u001b[0mvariable\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mvariable\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m~/anaconda3/lib/python3.9/site-packages/ibicus/debias/_debiaser.py\u001b[0m in \u001b[0;36m_from_variable\u001b[0;34m(child_class, variable, default_settings_variable, experimental_default_setting_variable, default_settings_general, **kwargs)\u001b[0m\n\u001b[1;32m 137\u001b[0m ]\n\u001b[1;32m 138\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 139\u001b[0;31m raise ValueError(\n\u001b[0m\u001b[1;32m 140\u001b[0m \u001b[0;34mf\"Unfortunately currently no default settings exist for the variable {variable} in the debiaser {child_class.__name__}. You can set the required class parameters manually by using the class constructor.\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 141\u001b[0m )\n", "\u001b[0;31mValueError\u001b[0m: Unfortunately currently no default settings exist for the variable tasmin in the debiaser ISIMIP. You can set the required class parameters manually by using the class constructor." ] } ], "source": [ "debiaser3 = ISIMIP.from_variable(\"tasmin\")" ] }, { "cell_type": "markdown", "id": "affead10-d3fd-4447-8a00-1b6a82e1bdfb", "metadata": {}, "source": [ "ISIMIP instead offers the option to debias `tasrange` and `tasskew` and calculate `tasmin` from those." ] }, { "cell_type": "markdown", "id": "5fcf6d02-cbad-4f99-b796-632821482c3b", "metadata": {}, "source": [ "Some bias adjustment methods offer and additional `for_precipitation` method to initialise it for precipitation (`pr`). Precipitation methods can be a bit more complicated and sometimes require the specification of a threshold under which precipitation is assumed to be zero. The `for_precipitation`-method is there to facilitate the choice of method.\n", "\n", "For example we can initialise a `QuantileMapping` debiaser with a precipitation gamma hurdle model. A hurdle model is a two step model where precipitation occurrence is modelled binomially and then a gamma distribution is assumed for the amounts. An alternative model is the censored model where all precipitation amounts under a threshold (so also all dry days) are assumed censored, so labeled 'not known' to the model.\n", "\n", "Let's initialize both:" ] }, { "cell_type": "code", "execution_count": 12, "id": "2062b508-acb6-4b3b-97db-31175cd04c8a", "metadata": {}, "outputs": [], "source": [ "from ibicus.debias import QuantileMapping\n", "\n", "# Initialise debiaser\n", "pr_debiaser_QM_hurdle = QuantileMapping.for_precipitation(model_type = \"hurdle\")\n", "pr_debiaser_QM_censored = QuantileMapping.for_precipitation(model_type = \"censored\", censoring_threshold = 0.1/86400)" ] }, { "cell_type": "markdown", "id": "a687e70a-2a00-4902-9120-ff62ba3f624d", "metadata": {}, "source": [ "### 2.2. Application" ] }, { "cell_type": "markdown", "id": "3bb595a4-da0d-4271-be83-d3639b8778ac", "metadata": {}, "source": [ "The bias adjustment method is applied through the apply-function. Lets initialise and apply a `DeltaChange` debiaser for `\"tas\"`:" ] }, { "cell_type": "code", "execution_count": 13, "id": "2ff9194b-afae-476e-859c-5aa49a795beb", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 225/225 [00:00<00:00, 4552.27it/s]\n" ] } ], "source": [ "tas_debiaser_DC = DeltaChange.from_variable(\"tas\")\n", "tas_debiased_DC = tas_debiaser_DC.apply(tas_obs, tas_cm_hist, tas_cm_future)" ] }, { "cell_type": "markdown", "id": "dfcbf921-e1c4-46ba-bdc8-e9eadcea9c67", "metadata": {}, "source": [ "As you can see, the apply function needs three numpy arrays of data in the form: \n", "\n", "- obs representing observations of a climatological variable\n", "- cm_hist representing climate model values for a climatological variable during a reference/the observational period\n", "- cm_future representing climate model values for a climatological variable during a future or application period that is to be debiased\n", "\n", "All three are assumed to be 3d-numpy arrays where the first dimension corresponds to time and the other two to spatial locations. The locations in obs, cm_hist and cm_future need to be the same and observational data has to be regridded to match the scale of the climate model grid prior to applying the bias adjustment.\n", "\n", "Besides obs, cm_hist and cm_future some debiasers might also require additional information like the dates to which observations and climate model values correspond, for example to apply the debiaser in a running window-mode. We have already seen an example in the first chapter of this notebook, when we initialized and applied ISIMIP. Dates datasets are arrays of dates in one of several formats:" ] }, { "cell_type": "code", "execution_count": 14, "id": "fae174d1-0db7-4668-8781-15a375c40422", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([cftime.DatetimeGregorian(1979, 1, 1, 0, 0, 0, 0, has_year_zero=False),\n", " cftime.DatetimeGregorian(1979, 1, 2, 0, 0, 0, 0, has_year_zero=False),\n", " cftime.DatetimeGregorian(1979, 1, 3, 0, 0, 0, 0, has_year_zero=False),\n", " ...,\n", " cftime.DatetimeGregorian(2005, 12, 29, 0, 0, 0, 0, has_year_zero=False),\n", " cftime.DatetimeGregorian(2005, 12, 30, 0, 0, 0, 0, has_year_zero=False),\n", " cftime.DatetimeGregorian(2005, 12, 31, 0, 0, 0, 0, has_year_zero=False)],\n", " dtype=object)" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "tas_time_obs" ] }, { "cell_type": "markdown", "id": "7b2dc041-1895-42a2-9e88-be8c58e524da", "metadata": {}, "source": [ "ISIMIP also runs without having been passed date arguments, by inference:" ] }, { "cell_type": "code", "execution_count": 15, "id": "58807008-b4f1-4a6c-8a9b-3b94e686c05f", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ " 0%| | 0/225 [00:00