{ "cells": [ { "attachments": {}, "cell_type": "markdown", "id": "18517e3a-e898-4a19-826b-c417fcbc5f47", "metadata": { "tags": [] }, "source": [ "# Continuous Ranked Probability Score (CRPS) for Ensembles" ] }, { "cell_type": "markdown", "id": "54e7f1a6", "metadata": {}, "source": [ "General information about the CRPS, and how to use `crps_cdf` with a forecast expressed as a cumulative distribution function (CDF), is available in the tutorial \"Continuous Ranked Probability Score\". \n", "We recommend you work through that tutorial before looking at the `crps_for_ensembles` and this tutorial." ] }, { "cell_type": "markdown", "id": "595cd38f", "metadata": {}, "source": [ "In weather forecasting it is common to run an ensemble of deterministic models. Each ensemble member provides a deterministic forecast. The idea is that the combined ensemble forecasts can be used to estimate a probability density function (PDF), and hence a CDF. \n", "\n", "The `crps_for_ensemble` function interprets ensemble forecast input as a CDF in one of two ways and calculates the CRPS as demonstrated below.\n", "\n", "A CDF is a function of the probability of non-exceedence. It is relevant to real-valued parameters such as temperature, rainfall amount, wind speed etc." ] }, { "cell_type": "code", "execution_count": 1, "id": "9859bbe2", "metadata": {}, "outputs": [], "source": [ "from scores.probability import crps_for_ensemble\n", "import matplotlib.pyplot as plt\n", "\n", "import numpy\n", "import xarray" ] }, { "cell_type": "markdown", "id": "c34ff15a", "metadata": {}, "source": [ "## Using the `crps_for_ensemble` function" ] }, { "cell_type": "code", "execution_count": 2, "id": "6e5251e4", "metadata": {}, "outputs": [], "source": [ "# Uncomment the following line and run to read the documentation regarding the crps_for_ensemble function in scores.\n", "# crps_for_ensemble?" ] }, { "cell_type": "markdown", "id": "149a6ba4", "metadata": {}, "source": [ "For the purpose of this tutorial, we will create a set of 10 equally likely forecasts of temperature which we will compare to a single verifying observation." ] }, { "cell_type": "code", "execution_count": 3, "id": "244c5cf2", "metadata": {}, "outputs": [], "source": [ "# In this example, each ensemble member is given a name 0 to 9\n", "ensemble_forecast = xarray.DataArray(\n", " [1.2, 2.0, 2.7, 2.9, 3.0, 3.0, 3.1, 3.2, 3.8, 5.0],\n", "\tcoords=[[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]], dims=[\"ensemble_member\"])\n", "\n", "# The observation is assumed to be 4.5\n", "obs_array = xarray.DataArray(4.5)" ] }, { "cell_type": "code", "execution_count": 4, "id": "e2078a9f", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAE8CAYAAAAsfWGYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABPdElEQVR4nO3dd1gU19cH8O+ywC4dpYMIqKAYEBAUkShGUbAQu0ZRAcsvRrARX1sU7CRWNKDGhjHBoFHTLBBE0USxIkZjBUWMSrOAgLCwO+8fhAkrIAsCuzDn8zwk7uydmTPDzuHunTv38hiGYUAIIYQTlOQdACGEkKZDSZ8QQjiEkj4hhHAIJX1CCOEQSvqEEMIhlPQJIYRDKOkTQgiHUNInhBAOoaRPCCEcQkm/mUtMTASPx0NiYqJM5fv06YM+ffo0WjyWlpbw9/dvtO03BR6Ph6CgIHmHUSd79+4Fj8dDenq6vEPhBH9/f2hqaspUlsfjYdmyZY0bUB1Q0q+Digurpp8LFy7IO8Rmpbi4GJs2bYKrqyt0dHQgFAphY2ODoKAg3Lt3jy23bNkyqfOsrq6Otm3bwsfHB1FRUSgpKamybX9//xp/T7GxsU15mJz0rvMvFArlHR6nKcs7gOZoxYoVsLKyqrK8Q4cOTR5L79698ebNG6iqqspU/vfff2/kiGSTm5sLb29vXL16FUOGDMH48eOhqamJu3fvIiYmBjt27IBIJJJaZ9u2bdDU1ERJSQmePHmCuLg4TJ48GeHh4Th69CjMzc2lygsEAuzatavKvh0cHBr12Ei5ms4/n8+XQzSkAiX9ehg4cCBcXFzkHQYAQElJSaaaU1FREdTV1WX+49DY/P39ce3aNRw6dAgjR46Uem/lypX44osvqqwzatQo6Ovrs69DQkIQHR2NSZMmYfTo0VW+aSkrK2PChAmNcwCkVnT+FRM17zSC9PR08Hg8rF+/HpGRkWjXrh3U1dUxYMAAPH78GAzDYOXKlWjTpg3U1NQwdOhQvHjxQmoblpaWGDJkCH7//Xc4OjpCKBSic+fOOHLkiFS56tr0+/TpAzs7O1y9ehW9e/eGuro6Fi9ezL73dpt+cXExli1bBhsbGwiFQpiYmGDEiBFIS0tjy6xfvx49e/aEnp4e1NTU4OzsjEOHDtXr/Fy8eBHHjh3DlClTqiR8oLyGuH79epm25evri6lTp+LixYuIj4+vVzw1iY6ORseOHSEUCuHs7IyzZ89Kvf/o0SPMmDEDHTt2hJqaGvT09DB69Ogq7eqlpaVYvnw5rK2tIRQKoaenhw8//LBKvHfu3MGoUaPQunVrCIVCuLi44Ndff60S199//42+fftCTU0Nbdq0wapVqyCRSGQ+rlOnTqFXr17Q0NCArq4uhg4ditu3b0uVqWhSS01Nhb+/P3R1daGjo4OAgAAUFRXJvK/aVDSZnjt3DsHBwTAwMICGhgaGDx+OnJwcqbJXrlyBl5cX9PX1oaamBisrK0yePFmqjEQiQXh4OD744AMIhUIYGRnh008/xcuXL6XKVVxfiYmJcHFxgZqaGuzt7dnr6MiRI7C3t2d/99euXas2/gcPHsDLywsaGhowNTXFihUrIMvAxU+ePMHkyZNhZGQEgUCADz74AHv27KnDmas/qunXQ15eHnJzc6WW8Xg86OnpSS2Ljo6GSCTCzJkz8eLFC6xduxZjxoxB3759kZiYiAULFiA1NRVff/015s2bV+WXfv/+fYwdOxbTp0+Hn58foqKiMHr0aMTGxqJ///7vjPH58+cYOHAgPvnkE0yYMAFGRkbVlhOLxRgyZAgSEhLwySefYPbs2Xj9+jXi4+Nx8+ZNtG/fHgCwefNmfPzxx/D19YVIJEJMTAxGjx6No0ePYvDgwXU6fxWJbOLEiXVaryYTJ07Ejh078Pvvv1c5L2//nlRUVKCjo1PrNs+cOYMDBw5g1qxZEAgE2Lp1K7y9vXHp0iXY2dkBAC5fvozz58/jk08+QZs2bZCeno5t27ahT58+uHXrFtTV1QGUJ9CwsDBMnToV3bt3R35+Pq5cuYLk5GQ23r///hvu7u4wMzPDwoULoaGhgYMHD2LYsGE4fPgwhg8fDgDIzMzERx99hLKyMrbcjh07oKamJtO5OnnyJAYOHIh27dph2bJlePPmDb7++mu4u7sjOTkZlpaWUuXHjBkDKysrhIWFITk5Gbt27YKhoSG++uormfb39vkHAFVVVWhra0stmzlzJlq1aoXQ0FCkp6cjPDwcQUFBOHDgAAAgOzsbAwYMgIGBARYuXAhdXV2kp6dXqQR9+umn2Lt3LwICAjBr1iw8fPgQERERuHbtGs6dOwcVFRW2bGpqKsaPH49PP/0UEyZMwPr16+Hj44Pt27dj8eLFmDFjBgAgLCwMY8aMwd27d6Gk9F89WSwWw9vbGz169MDatWsRGxuL0NBQlJWVYcWKFTWek6ysLPTo0YPtMGBgYIATJ05gypQpyM/Px5w5c2Q6t/XGEJlFRUUxAKr9EQgEbLmHDx8yABgDAwPm1atX7PJFixYxABgHBwemtLSUXT5u3DhGVVWVKS4uZpdZWFgwAJjDhw+zy/Ly8hgTExPGycmJXXb69GkGAHP69Gl2mYeHBwOA2b59e5Vj8PDwYDw8PNjXe/bsYQAwGzdurFJWIpGw/y4qKpJ6TyQSMXZ2dkzfvn2llltYWDB+fn5VtlXZ8OHDGQDMy5cv31muQmhoKAOAycnJqfb9ly9fMgCY4cOHs8v8/Pyq/T1VPvaaVJS9cuUKu+zRo0eMUCiU2sfb54RhGCYpKYkBwOzbt49d5uDgwAwePPid++zXrx9jb28v9RmQSCRMz549GWtra3bZnDlzGADMxYsX2WXZ2dmMjo4OA4B5+PDhO/fj6OjIGBoaMs+fP2eXXb9+nVFSUmImTZrELqs455MnT5Zaf/jw4Yyent4798EwNZ9/AIyXlxdbruKa8vT0lPq8zZ07l+Hz+ez189NPPzEAmMuXL9e4zz/++IMBwERHR0stj42NrbK84vo6f/48uywuLo4BwKipqTGPHj1il3/zzTdVrrGK45s5cya7TCKRMIMHD2ZUVVWlPqsAmNDQUPb1lClTGBMTEyY3N1cqzk8++YTR0dGp9nPVkKh5px4iIyMRHx8v9XPixIkq5UaPHi1Vq3R1dQUATJgwAcrKylLLRSIRnjx5IrW+qakpW8MDAG1tbUyaNAnXrl1DZmbmO2MUCAQICAio9VgOHz4MfX19zJw5s8p7PB6P/XflmuTLly+Rl5eHXr16ITk5udZ9vC0/Px8AoKWlVed1q1PRde7169dSy4VCYZXf04YNG2TappubG5ydndnXbdu2xdChQxEXFwexWAxA+pyUlpbi+fPn6NChA3R1daXOi66uLv7++2/cv3+/2n29ePECp06dwpgxY/D69Wvk5uYiNzcXz58/h5eXF+7fv89+No4fP44ePXqge/fu7PoGBgbw9fWt9ZiePXuGlJQU+Pv7o3Xr1uzyLl26oH///jh+/HiVdaZPny71ulevXnj+/Dn7O3yX6s5/fHw8vvzyyypl//e//0l93nr16gWxWIxHjx4BKD+HAHD06FGUlpZWu78ff/wROjo66N+/P3sOc3Nz4ezsDE1NTZw+fVqqfOfOneHm5sa+rrg++/bti7Zt21ZZ/uDBgyr7rNy1t6LmLhKJcPLkyWpjZBgGhw8fho+PDxiGkYrTy8sLeXl59bqm6oKad+qhe/fuMt3IrfzBAcD+AXi7l0nF8rfbHTt06CB1IQCAjY0NgPL7BsbGxjXu28zMTKabtmlpaejYsaPUH6HqHD16FKtWrUJKSopUF8m345NFxVf7169fsxfz+ygoKABQ9Y8In8+Hp6dnvbZpbW1dZZmNjQ2KioqQk5MDY2NjvHnzBmFhYYiKisKTJ0+k2nLz8vLYf69YsQJDhw6FjY0N7Ozs4O3tjYkTJ6JLly4AypsZGIbB0qVLsXTp0mrjyc7OhpmZGR49esQmoco6duxY6zFVJNDqytra2iIuLg6FhYXQ0NBgl7/9GW7VqhWA8s/q2000b6vL+X/XfgDAw8MDI0eOxPLly7Fp0yb06dMHw4YNw/jx4yEQCACUN4fm5eXB0NCw2n1kZ2e/c591vT6VlJTQrl07qWWVr8/q5OTk4NWrV9ixYwd27NghU5wNjZJ+I6qpa1pNy5kGnLlS1jZeWfzxxx/4+OOP0bt3b2zduhUmJiZQUVFBVFQU9u/fX+ftderUCQBw48YN9OrV673ju3nzJoCm7zI7c+ZMREVFYc6cOXBzc4OOjg54PB4++eQTqRurvXv3RlpaGn755Rf8/vvv2LVrFzZt2oTt27dj6tSpbNl58+bBy8ur2n3Jozsw0DSfVVn2w+PxcOjQIVy4cAG//fYb2113w4YNuHDhAjQ1NSGRSGBoaIjo6Ohqt2VgYCDTPhvzmCt+1xMmTICfn1+1ZSoqA42Fkr4Cq6gBVq5NVzy09PYNt/pq3749Ll68iNLSUqmbXJUdPnwYQqEQcXFxbK0KAKKiouq1Tx8fH4SFheH7779vkKT/3XffAUCNCbM+qmuKuXfvHtTV1dnkcejQIfj5+Uk1GRUXF+PVq1dV1m3dujUCAgIQEBCAgoIC9O7dG8uWLcPUqVPZ2qKKikqtNWMLC4tqY7t7926tx2RhYVFj2Tt37kBfX1+qlq+IevTogR49emD16tXYv38/fH19ERMTg6lTp6J9+/Y4efIk3N3dG7TSUxOJRIIHDx6wtXug9uvTwMAAWlpaEIvF9f4W+r6oTV+BPX36FD/99BP7Oj8/H/v27YOjo+M7m3bqYuTIkcjNzUVERESV9ypqNnw+Hzwej23LBsq/vv7888/12qebmxu8vb2xa9euarchEokwb948mba1f/9+7Nq1C25ubujXr1+94qlOUlKSVNvq48eP8csvv2DAgAFsTZDP51ep/X399ddS5wko70lVmaamJjp06MA2kxkaGqJPnz745ptv8OzZsyqxVO66OGjQIFy4cAGXLl2Ser+m2m1lJiYmcHR0xLfffiv1h+nmzZv4/fffMWjQoFq3IS8vX76scq4dHR0BgD2PY8aMgVgsxsqVK6usX1ZWVu0f4/dV+bphGAYRERFQUVGp8bPI5/MxcuRIHD58mP2GWtnb3VQbA9X06+HEiRO4c+dOleU9e/as0sb3PmxsbDBlyhRcvnwZRkZG2LNnD7Kysupdw67OpEmTsG/fPgQHB+PSpUvo1asXCgsLcfLkScyYMQNDhw7F4MGDsXHjRnh7e2P8+PHIzs5GZGQkOnTogL/++qte+923bx8GDBiAESNGwMfHB/369YOGhgbu37+PmJgYPHv2rEpf/UOHDkFTU5O96R0XF4dz587BwcEBP/74Y0OcDpadnR28vLykumwCwPLly9kyQ4YMwXfffQcdHR107twZSUlJOHnyZJWuu507d0afPn3g7OyM1q1b48qVKzh06JDUTcDIyEh8+OGHsLe3x7Rp09CuXTtkZWUhKSkJ//zzD65fvw4AmD9/Pr777jt4e3tj9uzZbJdNCwsLmX4X69atw8CBA+Hm5oYpU6awXTZ1dHQafHyYsrIyfP/999W+N3z48Dp9q/j222+xdetWDB8+HO3bt8fr16+xc+dOaGtrs3+sPDw88OmnnyIsLAwpKSkYMGAAVFRUcP/+ffz444/YvHkzRo0a1SDHBpTfqI6NjYWfnx9cXV1x4sQJHDt2DIsXL67SlFTZl19+idOnT8PV1RXTpk1D586d8eLFCyQnJ+PkyZNVntlpcI3aN6iFeVeXTQBMVFQUwzD/ddlct26d1PoV3St//PHHardbuTuahYUFM3jwYCYuLo7p0qULIxAImE6dOlVZt6Yumx988EG1x/B2l02GKe96+MUXXzBWVlaMiooKY2xszIwaNYpJS0tjy+zevZuxtrZm44iKimK79VUmS5fNyvtdv349061bN0ZTU5NRVVVlrK2tmZkzZzKpqalsuYr9VPwIhUKmTZs2zJAhQ5g9e/ZIdXOs4Ofnx2hoaMgUx9sAMIGBgcz333/PHrOTk5PUOWaY8q6iAQEBjL6+PqOpqcl4eXkxd+7cqXIOVq1axXTv3p3R1dVl1NTUmE6dOjGrV69mRCKR1PbS0tKYSZMmMcbGxoyKigpjZmbGDBkyhDl06JBUub/++ovx8PBghEIhY2ZmxqxcuZLZvXu3TF02GYZhTp48ybi7uzNqamqMtrY24+Pjw9y6dUuqTE3dZCs+q7Xt511dNiuvX91nn2Gqfq6Tk5OZcePGMW3btmUEAgFjaGjIDBkyRKpbbYUdO3Ywzs7OjJqaGqOlpcXY29sz8+fPZ54+fcqWqbi+3lbxu6+suuu54vOVlpbGDBgwgFFXV2eMjIyY0NBQRiwWV9lm5S6bDMMwWVlZTGBgIGNubs5ec/369WN27NjxzvPaEHj/BkUUjKWlJezs7HD06FF5h0IIaUGoTZ8QQjiEkj4hhHAIJX1CCOEQatMnhBAOoZo+IYRwCCV9QgjhEM49nCWRSPD06VNoaWnVa7AwQghRNAzD4PXr1zA1NZUa8786nEv6T58+rTKKHiGEtASPHz9GmzZt3lmGc0m/Yvjdx48f1zo0LCGENAf5+fkwNzeXaY4KziX9iiYdbW1tSvqEkBZFliZrupFLCCEcQkmfEEI4hJI+IYRwiFzb9M+ePYt169bh6tWrePbsGX766ScMGzbsneskJiYiODgYf//9N8zNzbFkyRL4+/s3aFwMw6CsrKzKZBiE1IbP50NZWZm6AxOFJdekX1hYCAcHB0yePBkjRoyotfzDhw8xePBgTJ8+HdHR0UhISMDUqVNhYmLSYFPliUQiPHv2DEVFRQ2yPcI96urqMDExkWliekKamsKMvcPj8Wqt6S9YsADHjh2Tmmbsk08+watXrxAbGyvTfvLz86Gjo4O8vLwqvXckEgnu378PPp8PAwMDqKqqUo2NyIxhGIhEIuTk5EAsFsPa2rrWB2UIaQjvymtva1ZdNpOSkqpMJuzl5YU5c+bUuE5JSQk7hyZQfnJqIhKJIJFIYG5uDnV19feOl3CPmpoaVFRU8OjRI4hEIgiFQnmHVKPiUjFm/nANj1803LfawKKtsCu73WDb4yrx0K3o4ODeKNtuVkk/MzMTRkZGUsuMjIyQn5+PN2/eQE1Nrco6YWFhUvOayoJqZ+R9NJfPz1//5CH+VlaDbU8PefARnmiw7XFZ7N2/KenX16JFixAcHMy+rnhyjRCuE0vKW3Y1eCX4UCX9vbfXCnnl22WUEMGbWG0ljNSOr8zHIPeGuUdZnWaV9I2NjZGVJV0zycrKgra2do0fMIFAAIFA0BThEdKsMChP+iqQwJT/+r23p8uUN52KeXz4z10JHR2d994maXjN43vov9zc3JCQkCC1LD4+Hm5ubnKKqPmxtLREeHi4vMNoMImJieDxeHj16pW8Q+E8PiQAADGUwOfz5RwNqYlck35BQQFSUlKQkpICoLxLZkpKCjIyMgCUN81MmjSJLT99+nQ8ePAA8+fPx507d7B161YcPHgQc+fOlUf4Cufx48eYPHkyTE1NoaqqCgsLC8yePRvPnz+Xd2gNok+fPlVu2vfs2RPPnj2jWmV9NHC/PaV/k74ESs3mvgYXyfU3c+XKFTg5OcHJyQkAEBwcDCcnJ4SEhAAAnj17xv4BAAArKyscO3YM8fHxcHBwwIYNG7Br164G66PfnD148AAuLi64f/8+fvjhB6SmpmL79u1ISEiAm5sbXrx4IZe4xGIxJBJJo21fVVUVxsbG1LVWAVSu6VPSV1xy/c306dMHDMNU+dm7dy8AYO/evUhMTKyyzrVr11BSUoK0tLQGfxr3bQzDoEhUJpefujxCERgYCFVVVfz+++/w8PBA27ZtMXDgQJw8eRJPnjzBF198wZZ9/fo1xo0bBw0NDZiZmSEyMlLqeJctW4a2bdtCIBDA1NQUs2bNYt8vKSnBvHnzYGZmBg0NDbi6ukr9jvbu3QtdXV38+uuv6Ny5MwQCAXbt2gWhUFilCWb27Nno27cvAOD58+cYN24czMzMoK6uDnt7e/zwww9sWX9/f5w5cwabN28Gj8cDj8dDenp6tc07hw8fxgcffACBQABLS0ts2LBBar+WlpZYs2YNJk+eDC0tLbRt2xY7duyQ+Vy3FEw1/3ofSih/gl0CPjXvKLBmdSNXHt6UitE5JE4u+761wgvqqrX/il68eIG4uDisXr26yg1tY2Nj+Pr64sCBA9i6dSsAYN26dVi8eDGWL1+OuLg4zJ49GzY2Nujfvz8OHz6MTZs2ISYmBh988AEyMzNx/fp1dntBQUG4desWYmJiYGpqip9++gne3t64ceMGrK2tAQBFRUX46quvsGvXLujp6aFNmzYICQnB4cOHMWXKFADl3wAOHDiA1atXAwCKi4vh7OyMBQsWQFtbG8eOHcPEiRPRvn17dO/eHZs3b8a9e/dgZ2eHFStWAAAMDAyQnp4udbxXr17FmDFjsGzZMowdOxbnz5/HjBkzoKenJ1VB2LBhA1auXInFixfj0KFD+Oyzz+Dh4YGOHTvW7ZdEWFTTbx4o6bcA9+/fB8MwsLW1rfZ9W1tbvHz5Ejk5OQAAd3d3LFy4EABgY2ODc+fOYdOmTejfvz8yMjJgbGwMT09PqKiooG3btujevTsAICMjA1FRUcjIyICpqSkAYN68eYiNjUVUVBTWrFkDACgtLcXWrVvh4ODAxvDJJ59g//79bNJPSEjAq1evMHLkSACAmZkZ5s2bx5afOXMm4uLicPDgQXTv3h06OjpQVVWFuro6jI2NazwXGzduRL9+/bB06VL2+G7duoV169ZJJf1BgwZhxowZAMqf9N60aRNOnz7NqaRf8UWyoRrGqE2/eaCkXws1FT5urZDPPQM1lbp9RZa1Oejt3k5ubm5sj57Ro0cjPDwc7dq1g7e3NwYNGgQfHx8oKyvjxo0bEIvFsLGxkVq/pKQEenp67GtVVVV06dJFqoyvry969OiBp0+fwtTUFNHR0Rg8eDB0dXUBlNf816xZg4MHD+LJkycQiUQoKSmp85PRt2/fxtChQ6WWubu7Izw8HGKxmG12qBwfj8eDsbExsrOz67QvIo1fKenTPRbFRUm/FjweT6YmFnnq0KEDeDwebt++jeHDh1d5//bt22jVqhUMDAxq3Za5uTnu3r2LkydPIj4+HjNmzMC6detw5swZFBQUgM/n4+rVq1XabDU1Ndl/q6mpVbnou3Xrhvbt2yMmJgafffYZfvrpJ/beDVDe5LR582aEh4fD3t4eGhoamDNnDkQiUR3PhmxUVFSkXvN4vEa94ayImAbuvvNfTZ/a8xWZYmczIhM9PT30798fW7duxdy5c6Xa9TMzMxEdHY1JkyaxifjChQtS61+4cEGqaUhNTQ0+Pj7w8fFBYGAgOnXqhBs3bsDJyQlisRjZ2dno1atXneP09fVFdHQ02rRpAyUlJQwePJh979y5cxg6dCgmTJgAoHzwu3v37qFz585sGVVV1VqHu7a1tcW5c+eklp07dw42NjZ0c7GRsUmfR007iox+Oy1EREQESkpK4OXlhbNnz+Lx48eIjY1F//79YWZmxt4wBcqT4Nq1a3Hv3j1ERkbixx9/xOzZswGU977ZvXs3bt68iQcPHuD777+HmpoaLCwsYGNjA19fX0yaNAlHjhzBw4cPcenSJYSFheHYsWO1xujr64vk5GSsXr0ao0aNknpS2traGvHx8Th//jxu376NTz/9tMrT15aWlrh48SLS09ORm5tbbc38888/R0JCAlauXIl79+7h22+/RUREhNT9AlKuodv0+ZV67xDFRUm/hbC2tsaVK1fQrl07jBkzBu3bt8f//vc/fPTRR0hKSkLr1q3Zsp9//jn7jMSqVauwceNG9lkHXV1d7Ny5E+7u7ujSpQtOnjyJ3377jW2zj4qKwqRJk/D555+jY8eOGDZsGC5fvoy2bdvWGmOHDh3QvXt3/PXXX/D19ZV6b8mSJejatSu8vLzQp08fGBsbVxlme968eeDz+ejcuTMMDAyknuGo0LVrVxw8eBAxMTGws7NDSEgIVqxY0ehde0nlmj4lfUWmMOPpN5V3jTtdXFyMhw8fwsrKSqGHxCWKrbl8js7cy4HfnktozSvCUOGt995eF+YWhiMW6fx2sFx6rQEiJLKqy3j6VNMnhDQIquk3D5T0CeGoii/5DdemX5H0qX+IIqOkTwhpEBXDMDBU01dolPQJ4aiGHnunoqZPSV+xUdInhDSIijZ9RomadxQZ/XYIURS5qcBvs4A3r5pkd91KynBCtQjKkECbKXnv7Wng3wnWqaav0CjpE6Iobv8KPDpXe7kGognAthG+6xeoGjb8RkmDoaRPiKIQl48zdBftcBFOjb67HIkGksvMoINizPqw5pFL60LCF8CsC01qpMgo6ROiKMSlAIA8aOMhz6LRd/cYOjgnsYa+UhE6DRrd6PsjioFu5BKiKCTlSV/cxJclDYLMLZT0Wwh/f392GsHKP6mpqfIOrV4qpl3kFHEZgPLx6AlpLNS804J4e3sjKipKapksY+i/TSQSQVVVtaHCIrJia/pN0/ulofvpk+aBqhS1YRhAVCifnzqOhScQCGBsbCz1w+fzcebMGXTv3h0CgQAmJiZYuHAhysrK2PX69OmDoKAgzJkzB/r6+uyImzdv3sTAgQOhqakJIyMjTJw4Ebm5uex6EokEa9euRYcOHSAQCNC2bVupIZwXLFgAGxsbqKuro127dli6dClKS0vZ969fv46PPvoIWlpa0NbWhrOzM65cuYLExEQEBAQgLy+P/caybNmyev4Cm5F/2/Sppk8aE9X0a1NaBKwxlc++Fz8FVDXeaxNPnjzBoEGD4O/vj3379uHOnTuYNm0ahEKhVCL99ttv8dlnn7ETkLx69Qp9+/bF1KlTsWnTJrx58wYLFizAmDFjcOrUKQDAokWLsHPnTmzatAkffvghnj17hjt37rDb1NLSwt69e2FqaoobN25g2rRp0NLSwvz58wGUj6/v5OSEbdu2gc/nIyUlBSoqKujZsyfCw8MREhKCu3fvApCemavFklebPjXqcwol/Rbk6NGjUslx4MCBsLGxgbm5OSIiIsDj8dCpUyc8ffoUCxYsQEhICDuBtbW1NdauXcuuu2rVKjg5ObGTnQPAnj17YG5ujnv37sHExASbN29GREQE/Pz8AADt27fHhx9+yJZfsmQJ+29LS0vMmzcPMTExbNLPyMjA//3f/6FTp05sDBV0dHTYuWs5g9r0SROgpF8bFfXyGre89l0HH330EbZt28a+1tDQQGBgINzc3KTmrHV3d0dBQQH++ecfdvITZ2dnqW1dv34dp0+frraGnZaWhlevXqGkpAT9+vWrMZ4DBw5gy5YtSEtLQ0FBAcrKyqTG+g4ODsbUqVPx3XffwdPTE6NHj0b79u3rdMwtSpO36VMVn4so6deGx3vvJpamoqGhgQ4dOtR73coKCgrg4+ODr776qkpZExMTPHjw4J3bS0pKgq+vL5YvXw4vLy/o6OggJiYGGzZsYMssW7YM48ePx7Fjx3DixAmEhoYiJiam2sndOUFCNX3S+Cjpt3C2trY4fPgwGIZha/vnzp2DlpYW2rRpU+N6Xbt2xeHDh2FpaQll5aofE2tra6ipqSEhIQFTp06t8v758+dhYWGBL774gl326NGjKuVsbGxgY2ODuXPnYty4cYiKisLw4cNlmgS9xfm3eYf66ZPGRFWKFm7GjBl4/PgxZs6ciTt37uCXX35BaGgogoOD2fb86gQGBuLFixcYN24cLl++jLS0NMTFxSEgIABisRhCoRALFizA/PnzsW/fPqSlpeHChQvYvXs3gPI/ChkZGYiJiUFaWhq2bNmCn376id3+mzdvEBQUhMTERDx69Ajnzp3D5cuXYWtrC6D8HkBBQQESEhKQm5uLoqKixj1RikBCvXdI46NPVwtnZmaG48eP49KlS3BwcMD06dMxZcoUqZus1TE1NcW5c+cgFosxYMAA2NvbY86cOdDV1WX/WCxduhSff/45QkJCYGtri7FjxyI7OxsA8PHHH2Pu3LkICgqCo6Mjzp8/j6VLl7Lb5/P5eP78OSZNmgQbGxuMGTMGAwcOxPLlywEAPXv2xPTp0zF27FgYGBhI3WRuscTy6b1DuIUmRq+kuUxoTRRbvT9He4cA6X/gEAbjb17HxgvwX4/Eujgl6gAjfiEurh7T6PsjjYcmRiekOZJTTZ/a9LmFkj4hiqKJ2/Q59RWfsOSe9CMjI2FpaQmhUAhXV1dcunTpneXDw8PRsWNHqKmpwdzcHHPnzkVxcXETRUtII5LXMAxU1ecUuSb9AwcOIDg4GKGhoUhOToaDgwO8vLzYm4Fv279/PxYuXIjQ0FDcvn0bu3fvxoEDB7B48eImjpyQRiCRT5dNwi1y/XRt3LgR06ZNQ0BAADp37ozt27dDXV0de/bsqbb8+fPn4e7ujvHjx8PS0hIDBgzAuHHjav12QEizIKeaPlX0uUVuD2eJRCJcvXoVixYtYpcpKSnB09MTSUlJ1a7Ts2dPfP/997h06RK6d++OBw8e4Pjx45g4cWKN+ykpKUFJyX+TPufn5zfcQRBuKc4HXj8DGMm7y5VKgNfZwP7FQHGWzJtnXqaDh6YbhoHSPTfJLenn5uZCLBbDyMhIarmRkZHUSI2VjR8/Hrm5ufjwww/BMAzKysowffr0dzbvhIWFsX2/CXkvRc/LR12tjZgpn+/2xX2g4LHMm+cBkICHPGjVP0ZCatGshmFITEzEmjVrsHXrVri6uiI1NRWzZ8/GypUrpR78qWzRokUIDg5mX+fn58Pc3LypQiYtSnl/lwKooxiCGkuVQowCFOFnDIAIr+q0h5fQwWte0yR96r3DTXJL+vr6+uDz+cjKkv76m5WVVeNwukuXLsXEiRPZsV7s7e1RWFiI//3vf/jiiy+qHVZAIBBAIKj5ApVVXl5ekw4FoK6uDh0dnSbbH5HBv88xlkEZItQ8s1gZylAGPv6BKQp5iv87pEYebpFb0ldVVYWzszMSEhIwbNgwAOUzMSUkJCAoKKjadYqKiqokdj6/vP2zMR8szsvLQ0REhNRsU41NWVkZQUFBzTbx37lzB/7+/khJSUGnTp2QkpIi75AaAFPpvy1AizkQUhdy7b0THByMnTt34ttvv8Xt27fx2WefobCwEAEBAQCASZMmSd3o9fHxwbZt2xATE4OHDx8iPj4eS5cuhY+PD5v8G0NRUVGTJnwAKCsra9aDjIWGhkJDQwN3795FQkJCg29fJBJh7dq1cHBwgLq6OvT19eHu7o6oqCh2SsbKk8WrqKjAysoK8+fPr/JcR+WJ5HV0dODu7s7ODgYAOTk5+Oyzz9DWoTcEVq7o5NgD48ePx+XLlxv8uAhpbHJt0x87dixycnIQEhKCzMxMODo6IjY2lr25m5GRIVWzX7JkCXg8HpYsWYInT57AwMAAPj4+UvOyEtmVlpZCRUWlUbadlpaGwYMHw8LCot7bqGmCdpFIBC8vL1y/fh0rV66Eu7s7tLW1ceHCBaxfvx5OTk5wdHQE8N9k8aWlpbh69Sr8/PzA4/GqzBMQFRUFb29v5Obm4osvvsCQIUNw8+ZNtGvXDiNHjoRIJMK3EV+inake7uWIcOrPS3j58mW9j00RUEWfm+T+FEhQUBAePXqEkpISXLx4Ea6urux7iYmJ2Lt3L/taWVkZoaGhSE1NxZs3b5CRkYHIyEjo6uo2feAKJjY2Fh9++CF0dXWhp6eHIUOGIC0tjX0/PT0dPB4PBw4cgIeHB4RCIaKjowEAu3btgq2tLYRCITp16oStW7dKbbu2Cc7fxuPxcPXqVaxYsUJqUvMbN26gb9++UFNTg56eHv73v/+hoKCAXc/f3x/Dhg3D6tWrYWpqio4dqx90LDw8HGfPnkVCQgICAwPh6OiIdu3aYfz48bh48aLUtIsVk8Wbm5tj2LBh8PT0RHx8fJVt6urqwtjYGHZ2dti2bRvevHmD+Ph4vHr1Cn/88Qe++uorfOTuCos2pujq5IiZM2diwIABtf9imgFq0+eWZtV7h9SssLAQwcHB6NKlCwoKChASEoLhw4cjJSVF6tvSwoULsWHDBjg5ObGJPyQkBBEREXBycsK1a9cwbdo0aGhosHPf1jbB+duePXsGT09PeHt7Y968edDU1ERhYSG8vLzg5uaGy5cvIzs7G1OnTkVQUJDUH/aEhARoa2tXm5grREdHw9PTE05OTlXeU1FRqfHby82bN9nJXd5FTU0NQPk3Ck1NTWhqauLnn39Gj/YBEMi9mkTI+6Gk30KMHDlS6vWePXtgYGCAW7duwc7Ojl0+Z84cjBgxgn0dGhqKDRs2sMusrKxw69YtfPPNN2zSr22C87cZGxtDWVkZmpqabE+snTt3ori4GPv27WOnZoyIiGCnZKxo0tPQ0MCuXbuqbdapcP/+ffTp00em81IxWXxZWRlKSkqgpKSEiIiIGssXFRVhyZIl4PP58PDwgLKyMvbu3Ytp06Zh+/Zt6GrXCd16uGPI0OHo3LmzTDEQokgo6bcQ9+/fR0hICC5evIjc3FxIJOVPjWZkZEglfRcXF/bfhYWFSEtLw5QpUzBt2jR2eVlZmVSvodomOJfF7du34eDgIDUXr7u7OyQSCe7evcsmfXt7+3cmfKBuPbUqJosvLCzEpk2boKysXOUPJACMGzcOfD4fb968gYGBAXbv3o0uXboAKP+DOnjwYPzx2/e4cCkZv52+gIht32DdunUYO3aszLEoHmrY4SJK+i2Ej48PLCwssHPnTpiamkIikcDOzg4ikUiqXOWkW9GevnPnTql7KcB/XWFlmeC8Ib09QXt1bGxsanxqu7rtVUwWv2fPHjg4OGD37t2YMmWKVLlNmzbB09MTOjo6MDAwqLIdoVAIz9490d/dGZ/NnY9Z8xZjw4YNzTzpl+NR7ucUaqFsAZ4/f467d+9iyZIl6NevH2xtbWXqWWJkZARTU1M8ePAAHTp0kPqxsrICID3BuYuLC6ytraud4Lw2tra2uH79OgoLC9ll586dg5KSUo03bGsyfvx4nDx5EteuXavyXmlpqdQ+KlNSUsLixYuxZMkSvHnzRuo9Y2NjdOjQodqE/5+Kfvo8WFtbN+sutQD13uEqSvotQKtWraCnp4cdO3YgNTUVp06dkhp64l2WL1+OsLAwbNmyBffu3cONGzcQFRWFjRs3Aqh9gnNZ+fr6QigUws/PDzdv3sTp06cxc+ZMTJw4scr4S7WZM2cO3N3d0a9fP0RGRuL69et48OABDh48iB49euD+/fs1rjt69Gjw+XxERkbKtK/nz5+jb9+++P777/HX33fwMOMJfv3tGLZt2wYvL686xU2IIqCkLwN1dXUoKzdtS5iysjLU1dVlKqukpISYmBhcvXoVdnZ2mDt3LtatWyfTulOnTsWuXbsQFRUFe3t7eHh4YO/evWxNv7YJzmWlrq6OuLg4vHjxAt26dcOoUaPQr1+/d95UrYlAIEB8fDzmz5+Pb775Bj169EC3bt2wZcsWzJo1S+oextsqnnReu3Ztjd8IKtPU1ISrqys2bdoEj+EBsOs7Gl+u24jx48dj1apVdY5dEfGozs8pNDF6Je+a0JrG3iHMs7/AY8TIhh7K3nE7rKysDE+ePMG5c+dk+sMiL2llrXG2tB3MVQrxx0qaGL05q8vE6HQjV0Y6OjqUhDnvvzZ9Qporat4hRFbc+lJMWihK+oTIrGXV9Cv+hFGbPrdQ0idEFkxLSfWE6yjpV4Nj97aJLCp9Jmr7dDSfzw/9GeMiupFbScVAXUVFReygW6QFe/MSeJ2Fuj+m9O5kWVpaCrFYXGXcfoVFj+RyCiX9Svh8PnR1dZGdnQ2gvNskjy6IlutVlmwTnVcihhJKIa72PYZhUFpaihcvXuDRo0cQi6svpyioTZ+bKOm/pWJUyIrET1qwgiygrATFELyz331lYiiBwZua3xeL8ejRI6SmpjZUlIQ0KEr6b+HxeDAxMYGhoeE7JwohLcChVUDmdRzHR3iA+s/wVVlxcbHC1/AJt1HSrwGfz2/UeXeJAijOBgoeQ4RXKOTpyzsaOaImTC6h3juEw6gtG6A2fa6hpE+4i2lZD1sRIgtK+oS7mPLZxbia9Kl+z02U9AmHUdoj3ENJn3AX27zDVeXfcOhRFG6pV9I/ffp0Q8dBiBz893gSIVxRr6Tv7e2N9u3bY9WqVXj8+HFDx0RI02Db9LmN/uRxS72S/pMnTxAUFIRDhw6hXbt28PLywsGDByESiRo6PkIaD1PxP0p7hDvqlfT19fUxd+5cpKSk4OLFi7CxscGMGTNgamqKWbNm4fr16w0dJyGNgJp3CPe8943crl27YtGiRQgKCkJBQQH27NkDZ2dn9OrVC3///XdDxEhI4+D4jVyuHjfX1Tvpl5aW4tChQxg0aBAsLCwQFxeHiIgIZGVlITU1FRYWFhg9enRDxkpIA6OaPsD1o+eeeo29M3PmTPzwww9gGAYTJ07E2rVrYWdnx76voaGB9evXw9TUtMECJaTB0Y1cwkH1Svq3bt3C119/jREjRkAgEFRbRl9fn7p2EsVGwzCU4/jhc029mndCQ0MxevToKgm/rKwMZ8+eBQAoKyvDw8Oj1m1FRkbC0tISQqEQrq6uuHTp0jvLv3r1CoGBgTAxMYFAIICNjQ2OHz9en8MgnMftOj63j5676pX0P/roI7x48aLK8ry8PHz00Ucyb+fAgQMIDg5GaGgokpOT4eDgAC8vrxonMBGJROjfvz/S09Nx6NAh3L17Fzt37oSZmVl9DoNwHdX0AVBFn2vq1bzDMEy10wg+f/4cGhoaMm9n48aNmDZtGgICAgAA27dvx7Fjx7Bnzx4sXLiwSvk9e/bgxYsXOH/+PDufraWlZX0OgRDOD7hG6Z6b6pT0R4wYAaB8dil/f3+p5h2xWIy//voLPXv2lGlbIpEIV69exaJFi9hlSkpK8PT0RFJSUrXr/Prrr3Bzc0NgYCB++eUXGBgYYPz48ViwYEGNE56UlJSgpKSEfZ2fny9TfKTlYzic7ivj0VnglDolfR0dHQDlNX0tLS2oqamx76mqqqJHjx6YNm2aTNvKzc2FWCyGkZGR1HIjIyPcuXOn2nUePHiAU6dOwdfXF8ePH0dqaipmzJiB0tJShIaGVrtOWFgYli9fLlNMhGOoeYdwUJ2SflRUFIDyJpV58+bVqSmnIUgkEhgaGmLHjh3g8/lwdnbGkydPsG7duhqT/qJFixAcHMy+zs/Ph7m5eVOFTBQa3cosR+eBS+rVpl9Tgq0LfX198Pl8ZGVlSS3PysqCsbFxteuYmJhARUVFqinH1tYWmZmZEIlEUFVVrbKOQCCosVsp4TiO99Pn6nFzncxJv2vXrkhISECrVq3g5ORU7Y3cCsnJybVuT1VVFc7OzkhISMCwYcMAlNfkExISEBQUVO067u7u2L9/PyQSCZSUyjse3bt3DyYmJtUmfELeiZp3AOCd1zJpeWRO+kOHDmVrzBVJ+n0FBwfDz88PLi4u6N69O8LDw1FYWMj25pk0aRLMzMwQFhYGAPjss88QERGB2bNnY+bMmbh//z7WrFmDWbNmNUg8hGMYGoaBcI/MSb9yk05DNO8AwNixY5GTk4OQkBBkZmbC0dERsbGx7M3djIwMtkYPAObm5oiLi8PcuXPRpUsXmJmZYfbs2ViwYEGDxEO4htsDrv2HzgCX1KtNvyEFBQXV2JyTmJhYZZmbmxsuXLjQyFERTuB4TZ9SPTfJnPRbtWolc9tfdU/rEqJwOP9wVjnqp88tMif98PDwRgyDEHngevMOJXsukjnp+/n5NWYchDQ9jjfvVOD20XOPzEk/Pz8f2tra7L/fpaIcIYqN6zV9wkV1atN/9uwZDA0NoaurW237fsVAbGKxuEGDJKRRUJt+OY4fPtfInPRPnTqF1q1bAwBNjkJaBobbdXxuHz13yZz0K0+IIsvkKIQ0F1yv6XP76Lmn3v30X758id27d+P27dsAgM6dOyMgIID9NkCIwuN4TZ9wU71mzjp79iwsLS2xZcsWvHz5Ei9fvsSWLVtgZWXFTpdIiOKjsXcAqulzTb1q+oGBgRg7diy2bdvGjngpFosxY8YMBAYG4saNGw0aJCGNgvM3crl63NxWr5p+amoqPv/8c6khjvl8PoKDg5GamtpgwRHSqKh5BwBAg2xyS72SfteuXdm2/Mpu374NBweH9w6KkKbB7X769DePm2Ru3vnrr7/Yf8+aNQuzZ89GamoqevToAQC4cOECIiMj8eWXXzZ8lIQ0Bnoil3CQzEnf0dERPB4PTKXqwfz586uUGz9+PMaOHdsw0RHSmDjfpk+4SOak//Dhw8aMgxA54HbzTgVq0+cWmZO+hYVFY8ZBSNPjePMO1//YcdV7TaJy69YtZGRkQCQSSS3/+OOP3ysoQpoSJT/CJfVK+g8ePMDw4cNx48YNqXb+ikHYaMA10ixQmz7hoHp12Zw9ezasrKyQnZ0NdXV1/P333zh79ixcXFyqneKQEMVEdXyAq41b3FWvmn5SUhJOnToFfX19KCkpQUlJCR9++CHCwsIwa9YsXLt2raHjJKThMVwfhoGrx81t9arpi8ViaGlpAQD09fXx9OlTAOU3e+/evdtw0RHSqKimD1Dq55p61fTt7Oxw/fp1WFlZwdXVFWvXroWqqip27NiBdu3aNXSMhDQOjtf06U8eN9Ur6S9ZsgSFhYUAgBUrVmDIkCHo1asX9PT0cODAgQYNkJDGwqNRNstRR31OqVfS9/LyYv/doUMH3LlzBy9evECrVq2qnUaREIVDA88QjnqvfvoA8PjxYwCAubn5ewdDSJOplPS5XtPnUUMPp9TrRm5ZWRmWLl0KHR0dWFpawtLSEjo6OliyZAlKS0sbOkZCGgElOsJN9arpz5w5E0eOHMHatWvh5uYGoLwb57Jly/D8+XNs27atQYMkpMH9+2AWQDV9apLllnol/f379yMmJgYDBw5kl3Xp0gXm5uYYN24cJX2i+KSadwjhjno17wgEAlhaWlZZbmVlBVVV1feNiZAmUDnVc7Omy+3h5rirXkk/KCgIK1euRElJCbuspKQEq1evRlBQUIMFR0ijoZo+4SiZm3dGjBgh9frkyZNo06YNOz3i9evXIRKJ0K9fv4aNkJDGUKlNn7t1Xa4eN7fJXNPX0dGR+hk5ciSGDBkCc3NzmJubY8iQIRgxYgR0dHTqHERkZCQsLS0hFArh6uqKS5cuybReTEwMeDwehg0bVud9Eq6jLpuEm2Su6UdFRTVKAAcOHEBwcDC2b98OV1dXhIeHw8vLC3fv3oWhoWGN66Wnp2PevHno1atXo8RFWjhq3qE2fY6qV5t+hZycHPz555/4888/kZOTU69tbNy4EdOmTUNAQAA6d+6M7du3Q11dHXv27KlxHbFYDF9fXyxfvpzG+iH1xNVUXw3K+pxSr6RfWFiIyZMnw8TEBL1790bv3r1hamqKKVOmoKioSObtiEQiXL16FZ6env8FpKQET09PJCUl1bjeihUrYGhoiClTptS6j5KSEuTn50v9EEJP5BKuqlfSDw4OxpkzZ/Dbb7/h1atXePXqFX755RecOXMGn3/+uczbyc3NhVgshpGRkdRyIyMjZGZmVrvOn3/+id27d2Pnzp0y7SMsLEzqXgQNF0EA0MNZlXD76LmnXkn/8OHD2L17NwYOHAhtbW1oa2tj0KBB2LlzJw4dOtTQMbJev36NiRMnYufOndDX15dpnUWLFiEvL4/9qRgriHAdNe8QbqrXE7lFRUVVaucAYGhoWKfmHX19ffD5fGRlZUktz8rKgrGxcZXyaWlpSE9Ph4+PD7tMIimvsSkrK+Pu3bto37691DoCgQACgUDmmAhHUPMO4ah61fTd3NwQGhqK4uJidtmbN2+wfPlydiweWaiqqsLZ2RkJCQnsMolEgoSEhGq306lTJ9y4cQMpKSnsz8cff4yPPvoIKSkp1HRDSB3Qdx1uqldNPzw8HN7e3lUezhIKhYiLi6vTtoKDg+Hn5wcXFxd0794d4eHhKCwsREBAAABg0qRJMDMzQ1hYGIRCIezs7KTW19XVBYAqywl5J2rTZ9GAa9xSr6Rvb2+P+/fvIzo6Gnfu3AEAjBs3Dr6+vlBTU6vTtsaOHYucnByEhIQgMzMTjo6OiI2NZZuPMjIyoKT0Xj1LCamq8iQqnE16XD1ubuMxTN2mECotLUWnTp1w9OhR2NraNlZcjSY/Px86OjrIy8uDtra2vMMh8lKQDay3BgAs5wXLORj5uFZqipQyU3TTfYMfF46SdzjkPdQlr9W5Cq2ioiLVlk9Is8ROik4It9Sr3SQwMBBfffUVysrKGjoeQprGv236XG/PBzjcusVR9WrTv3z5MhISEvD777/D3t4eGhoaUu8fOXKkQYIjpPFU1PS5m/HoWw431Svp6+rqYuTIkQ0dCyFNp263sghpMeqU9CUSCdatW4d79+5BJBKhb9++WLZsWZ177BAif1TTJ9xUpzb91atXY/HixdDU1ISZmRm2bNmCwMDAxoqNkMZDNX0WtelzS52S/r59+7B161bExcXh559/xm+//Ybo6Gh2KARCmg26kUtt+hxVp6SfkZGBQYMGsa89PT3B4/Hw9OnTBg+MkMZFXTYr8Dj8h4+L6pT0y8rKIBQKpZapqKigtLS0QYMipNExNG8Ut4+du+p0I5dhGPj7+0uNWllcXIzp06dLddukLptE8VFNvwKlfm6pU9L38/OrsmzChAkNFgwhTYah3juEm+qU9BtrcnRCmhw177Co9w630PCVhKOoeYfLx85l9Xoil5Bm4Vo0cGGb1Nj5rDIaNLACVfS5hZI+abkubAWybr6zyEvoNFEwhCgGSvqk5SorAQDEwQNZ0K+2yFNUnYuZc6iqzymU9EmLxUhKwQPwD0zwD89U3uEoHGrT5ya6kUtaLnH5Q4Ni+pi/Ez2Ryy10NZCW69+kLwFfzoEoKkr2XERJn7RckvKZ3aim/26U+rmFrgbScrE1ffqYE1KBrgbSckmoTV8W9EQut9DVQFquf5t3qKZfPeq9w010NZCWiWHAo6QvE6rpcwtdDaRlkojZf1LzDiH/oauBtEyS/yb2oS6bhPyHkj5pmcT/JX2q6RPyH7oaSMv0b3s+QG36Nfr3Ti6PGvU5ha4G0jL9W9NnADA8+pgTUoGuBtIysX30qT2/JhVTRVI9n1so6ZOWiZ7GJaRaCnFFREZGwtLSEkKhEK6urrh06VKNZXfu3IlevXqhVatWaNWqFTw9Pd9ZnnAUjbsjM6rpc4vcr4gDBw4gODgYoaGhSE5OhoODA7y8vJCdnV1t+cTERIwbNw6nT59GUlISzM3NMWDAADx58qSJIycKjWr6hFRL7lfExo0bMW3aNAQEBKBz587Yvn071NXVsWfPnmrLR0dHY8aMGXB0dESnTp2wa9cuSCQSJCQkNHHkRKHRuDu1YodhoKo+p8h15iyRSISrV69i0aJF7DIlJSV4enoiKSlJpm0UFRWhtLQUrVu3rvb9kpISlJSUsK/z8/PfL2hSf3FfAGmnm2ZfpUUAqKZPyNvkmvRzc3MhFothZGQktdzIyAh37tyRaRsLFiyAqakpPD09q30/LCwMy5cvf+9YyXsSFQJJEU2+25fQbfJ9Njc0cxa3NOs5cr/88kvExMQgMTERQqGw2jKLFi1CcHAw+zo/Px/m5uZNFSKpIBax//weIyBpkkTDwz808TkhUuSa9PX19cHn85GVlSW1PCsrC8bG775Y169fjy+//BInT55Ely5daiwnEAggEAgaJF7yHsT/PSGbBgsa2lGB0K+CW+Ta4KmqqgpnZ2epm7AVN2Xd3NxqXG/t2rVYuXIlYmNj4eLi0hShkvdV+cYqZRlC5EbuzTvBwcHw8/ODi4sLunfvjvDwcBQWFiIgIAAAMGnSJJiZmSEsLAwA8NVXXyEkJAT79++HpaUlMjMzAQCamprQ1NSU23GQWlAXSoVDT+Ryk9yT/tixY5GTk4OQkBBkZmbC0dERsbGx7M3djIwMKCn9lyi2bdsGkUiEUaNGSW0nNDQUy5Yta8rQSV3Qw1KEKAS5J30ACAoKQlBQULXvJSYmSr1OT09v/IBIw6OaPiEKga5A0jRoADRCFAIlfdI0qKavsOi+OrfQFUiaBk1STohCoCuQNA26katwKsbeoYo+t9AVSJoGNe8QohDoCiRNg0a9VFjUps8tdAWSpiGuaNOn3juKh7I+l1DSJ01DQs07ioaeyOUmugJJ0xBT8w4hioCuQNI0JNS8o7Coqs8plPRJ06CavsKinM8tdAWSpkFt+oQoBLoCSdOgfvqEKASFGGWTNIGS18ChyUDeE/nsv+g5AEDCo6SvKOiJXG6ipM8Vj84D93+XdxR4wehSllE09HQWp1DS54qyEgBANvQQiz7yCQHK+Acmctk3qRmlfG6hpM8V/95ILYI6HvIs5BwMIUReqIGVK8Q0yiV5G9XxuYgyAFdQl0lSA2rS5xbKAFzBdpmkJ2JJOeq9w02U9LmCJjEhhICSPndU1PSpnzwhnEYZgCtoEhNCCCjpc0dFTZ+hXzkpx9RehLRAlAG4QiIGQDV9UhWPuu9wCmUArqAum4QQUNLnDhrPntSA6vncQhmAK9iZq+hXTgiXUQbgCramTw9nkQr/ToxOVX1OoaTPFdSmTwgBJX3uoDZ98hbqsslNCpEBIiMjYWlpCaFQCFdXV1y6dOmd5X/88Ud06tQJQqEQ9vb2OH78eBNF2oxRmz4hBAqQ9A8cOIDg4GCEhoYiOTkZDg4O8PLyQnZ2drXlz58/j3HjxmHKlCm4du0ahg0bhmHDhuHmzZtNHHkzQ236pAZK1KjPKXJP+hs3bsS0adMQEBCAzp07Y/v27VBXV8eePXuqLb9582Z4e3vj//7v/2Bra4uVK1eia9euiIiIaOLImxlq0yeEQM4zZ4lEIly9ehWLFi1ilykpKcHT0xNJSUnVrpOUlITg4GCpZV5eXvj555+rLV9SUoKSkhL2dX5+fr1ivbDtUxjlXKjXuorAUJIFDQDJpWb4WdJZ3uEQBVDIqMo7BCIHck36ubm5EIvFMDIyklpuZGSEO3fuVLtOZmZmteUzMzOrLR8WFobly5e/d6yqBU9hJUl/7+3I2y3GHC8ZdXmHQRSIvho173BJi58jd9GiRVLfDPLz82Fubl7n7egOXILYvy4j4/HjhgyvSRUoacGtsxP6qlPSJ+W0VHkY0I2++XGJXJO+vr4++Hw+srKypJZnZWXB2Ni42nWMjY3rVF4gEEAgELx3rO3sXNHOzvW9t0MIIfIk17t6qqqqcHZ2RkJCArtMIpEgISEBbm5u1a7j5uYmVR4A4uPjayxPCCHkP3Jv3gkODoafnx9cXFzQvXt3hIeHo7CwEAEBAQCASZMmwczMDGFhYQCA2bNnw8PDAxs2bMDgwYMRExODK1euYMeOHfI8DEIIaRbknvTHjh2LnJwchISEIDMzE46OjoiNjWVv1mZkZEBJ6b8vJD179sT+/fuxZMkSLF68GNbW1vj5559hZ2cnr0MghJBmg8cwDKeexs7Pz4eOjg7y8vKgra0t73AIIeS91SWv0ZM6hBDCIZT0CSGEQyjpE0IIh8j9Rm5Tq7iFUd/hGAghRNFU5DNZbtFyLum/fv0aAOr1VC4hhCiy169fQ0dH551lONd7RyKR4OnTp9DS0gKvDkPKVgzf8PjxY+r1U090Dt8Pnb/311LPIcMweP36NUxNTaW6uFeHczV9JSUltGnTpt7ra2trt6gPizzQOXw/dP7eX0s8h7XV8CvQjVxCCOEQSvqEEMIhlPRlJBAIEBoa2iAjdnIVncP3Q+fv/dE55OCNXEII4TKq6RNCCIdQ0ieEEA6hpE8IIRxCSZ8QQjiEkr6MIiMjYWlpCaFQCFdXV1y6dEneITUbYWFh6NatG7S0tGBoaIhhw4bh7t278g6r2fryyy/B4/EwZ84ceYfSbDx58gQTJkyAnp4e1NTUYG9vjytXrsg7LLmgpC+DAwcOIDg4GKGhoUhOToaDgwO8vLyQnZ0t79CahTNnziAwMBAXLlxAfHw8SktLMWDAABQWFso7tGbn8uXL+Oabb9ClSxd5h9JsvHz5Eu7u7lBRUcGJEydw69YtbNiwAa1atZJ3aHJBXTZl4Orqim7duiEiIgJA+fg95ubmmDlzJhYuXCjn6JqfnJwcGBoa4syZM+jdu7e8w2k2CgoK0LVrV2zduhWrVq2Co6MjwsPD5R2Wwlu4cCHOnTuHP/74Q96hKASq6ddCJBLh6tWr8PT0ZJcpKSnB09MTSUlJcoys+crLywMAtG7dWs6RNC+BgYEYPHiw1GeR1O7XX3+Fi4sLRo8eDUNDQzg5OWHnzp3yDktuKOnXIjc3F2KxmJ2ovYKRkREyMzPlFFXzJZFIMGfOHLi7u9Nk9nUQExOD5ORkhIWFyTuUZufBgwfYtm0brK2tERcXh88++wyzZs3Ct99+K+/Q5IJzo2wS+QoMDMTNmzfx559/yjuUZuPx48eYPXs24uPjIRQK5R1OsyORSODi4oI1a9YAAJycnHDz5k1s374dfn5+co6u6VFNvxb6+vrg8/nIysqSWp6VlQVjY2M5RdU8BQUF4ejRozh9+vR7DW/NNVevXkV2dja6du0KZWVlKCsr48yZM9iyZQuUlZUhFovlHaJCMzExQefOnaWW2draIiMjQ04RyRcl/VqoqqrC2dkZCQkJ7DKJRIKEhAS4ubnJMbLmg2EYBAUF4aeffsKpU6dgZWUl75CalX79+uHGjRtISUlhf1xcXODr64uUlBTw+Xx5h6jQ3N3dq3QRvnfvHiwsLOQUkXxR844MgoOD4efnBxcXF3Tv3h3h4eEoLCxEQECAvENrFgIDA7F//3788ssv0NLSYu+F6OjoQE1NTc7RKT4tLa0q9z80NDSgp6dH90VkMHfuXPTs2RNr1qzBmDFjcOnSJezYsQM7duyQd2jywRCZfP3110zbtm0ZVVVVpnv37syFCxfkHVKzAaDan6ioKHmH1mx5eHgws2fPlncYzcZvv/3G2NnZMQKBgOnUqROzY8cOeYckN9RPnxBCOITa9AkhhEMo6RNCCIdQ0ieEEA6hpE8IIRxCSZ8QQjiEkj4hhHAIJX1CCOEQSvqEEMIhlPQJIYRDKOmTZoXH473zZ9myZfIOscFZWlrSDFmkwdCAa6RZefbsGfvvAwcOICQkRGoERU1NTXmEVWcMw0AsFkNZuekuQZFIBFVV1SbbH1FMVNMnzYqxsTH7o6OjAx6PJ7UsJiYGtra2EAqF6NSpE7Zu3cqum56eDh6Ph4MHD6JXr15QU1NDt27dcO/ePVy+fBkuLi7Q1NTEwIEDkZOTw67n7++PYcOGYfny5TAwMIC2tjamT58OkUjElpFIJAgLC4OVlRXU1NTg4OCAQ4cOse8nJiaCx+PhxIkTcHZ2hkAgwJ9//om0tDQMHToURkZG0NTURLdu3XDy5El2vT59+uDRo0eYO3cu+20GAJYtWwZHR0epcxMeHg5LS8sqca9evRqmpqbo2LEjgPJJWcaMGQNdXV20bt0aQ4cORXp6ekP8ekgzQEmftBjR0dEICQnB6tWrcfv2baxZswZLly6tMi1eaGgolixZguTkZCgrK2P8+PGYP38+Nm/ejD/++AOpqakICQmRWichIQG3b99GYmIifvjhBxw5cgTLly9n3w8LC8O+ffuwfft2/P3335g7dy4mTJiAM2fOSG1n4cKF+PLLL3H79m106dIFBQUFGDRoEBISEnDt2jV4e3vDx8eHneDjyJEjaNOmDVasWIFnz55JfdORRUJCAu7evYv4+HgcPXoUpaWl8PLygpaWFv744w+cO3cOmpqa8Pb2lvojRlowOY/ySUi9RUVFMTo6Ouzr9u3bM/v375cqs3LlSsbNzY1hGIZ5+PAhA4DZtWsX+/4PP/zAAGASEhLYZWFhYUzHjh3Z135+fkzr1q2ZwsJCdtm2bdsYTU1NRiwWM8XFxYy6ujpz/vx5qX1PmTKFGTduHMMwDHP69GkGAPPzzz/XelwffPAB8/XXX7OvLSwsmE2bNkmVCQ0NZRwcHKSWbdq0ibGwsJCK28jIiCkpKWGXfffdd0zHjh0ZiUTCLispKWHU1NSYuLi4WmMjzR+16ZMWobCwEGlpaZgyZQqmTZvGLi8rK4OOjo5U2S5durD/rpjw3t7eXmpZdna21DoODg5QV1dnX7u5uaGgoACPHz9GQUEBioqK0L9/f6l1RCIRnJycpJa5uLhIvS4oKMCyZctw7NgxPHv2DGVlZXjz5k2DTeVnb28v1Y5//fp1pKamQktLS6pccXEx0tLSGmSfRLFR0ictQkFBAQBg586dcHV1lXrv7ekEVVRU2H9XtJG/vUwikdR538eOHYOZmZnUewKBQOq1hoaG1Ot58+YhPj4e69evR4cOHaCmpoZRo0bV2tSipKQE5q2pMEpLS6uUe3t/BQUFcHZ2RnR0dJWyBgYG79wnaRko6ZMWwcjICKampnjw4AF8fX0bfPvXr1/Hmzdv2OkdL1y4AE1NTZibm6N169YQCATIyMiAh4dHnbZ77tw5+Pv7Y/jw4QDKk/LbN1VVVVWrTH5uYGCAzMxMMAzD/uFKSUmpdX9du3bFgQMHYGhoCG1t7TrFSloGupFLWozly5cjLCwMW7Zswb1793Djxg1ERUVh48aN771tkUiEKVOm4NatWzh+/DhCQ0MRFBQEJSUlaGlpYd68eZg7dy6+/fZbpKWlITk5GV9//XWVm8hvs7a2xpEjR5CSkoLr169j/PjxVb5lWFpa4uzZs3jy5Alyc3MBlPfqycnJwdq1a5GWlobIyEicOHGi1uPw9fWFvr4+hg4dij/++AMPHz5EYmIiZs2ahX/++af+J4g0G5T0SYsxdepU7Nq1C1FRUbC3t4eHhwf27t0LKyur9952v379YG1tjd69e2Ps2LH4+OOPpR4EW7lyJZYuXYqwsDDY2trC29sbx44dq3XfGzduRKtWrdCzZ0/4+PjAy8sLXbt2lSqzYsUKpKeno3379mwTjK2tLbZu3YrIyEg4ODjg0qVLmDdvXq3Hoa6ujrNnz6Jt27YYMWIEbG1tMWXKFBQXF1PNnyNojlxCauHv749Xr17h559/lncohLw3qukTQgiHUNInhBAOoeYdQgjhEKrpE0IIh1DSJ4QQDqGkTwghHEJJnxBCOISSPiGEcAglfUII4RBK+oQQwiGU9AkhhEP+H9sLUKmRRmtNAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# The ensemble forecast can be converted to a CDF\n", "# When converted in a naive way we get the empirical CDF illustrated here. The larger the ensemble the smoother and more sensible this will be.\n", "# The plot below also shows the CDF corresponding to the observation, and the area corrsponding to the CRPS.\n", "fcst_thresholds = numpy.linspace(0, 7, 700) \n", "empirical_cdf = xarray.DataArray(coords={'temperature': fcst_thresholds}, data=[0]*120 + [0.1]*80 + [0.2]*70 + [0.3]*20 +[0.7]*30 + [0.8]*60 + [0.9]*120 +[1]*200)\n", "observed_cdf = numpy.heaviside(fcst_thresholds-4.5, 1)\n", "plt.figure(figsize=(4, 3))\n", "plt.plot(fcst_thresholds, observed_cdf, label='Observation')\n", "plt.plot(fcst_thresholds, empirical_cdf, label='Forecast')\n", "plt.fill_between(fcst_thresholds, empirical_cdf, observed_cdf, color='gray', label='area for CRPS')\n", "plt.title(\"Empirical CDF based on Ensemble\")\n", "plt.xlabel(\"Temperature\")\n", "plt.ylabel(\"Probability\")\n", "plt.legend(loc=\"upper left\")" ] }, { "cell_type": "code", "execution_count": 5, "id": "c11d30f5", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray ()>\n",
       "array(1.111)
" ], "text/plain": [ "\n", "array(1.111)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Specifying method='ecdf' assumes the empirical CDF\n", "crps_for_ensemble(ensemble_forecast, obs_array, ensemble_member_dim='ensemble_member', method='ecdf').round(3)" ] }, { "cell_type": "code", "execution_count": 6, "id": "1d296010", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray ()>\n",
       "array(1.056)
" ], "text/plain": [ "\n", "array(1.056)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Specifying method='fair' gives an approximated CRPS assuming that the ensemble values can be interpreted as a random sample from an underlying predictive distribution\n", "crps_for_ensemble(ensemble_forecast, obs_array, ensemble_member_dim='ensemble_member', method='fair').round(3)" ] }, { "cell_type": "markdown", "id": "2517d17e", "metadata": {}, "source": [ "### Things to try next\n", "\n", "* Use crps_for_ensemble with some real forecasts with extra dimensions in space and/or time. \n", "* Read about the 'fair' method in C. Ferro (2014), \"Fair scores for ensemble forecasts\", Q J R Meteorol Soc\n", " 140(683):1917-1923. https://doi.org/10.1002/qj.2270" ] }, { "cell_type": "code", "execution_count": null, "id": "2b4027aa", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "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.11.6" } }, "nbformat": 4, "nbformat_minor": 5 }