SurpriseRemoval_2022.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Surprise Removal\n",
    "## John R. Lawson, 2022\n",
    "#### (Supporting material for paper in preparation)\n",
    "##### Coauthors: Corey K. Potvin and Kenric Nelson. Code by JRL. \n",
    "##### Work supported by CIMMS/NSSL, Oklahoma; and Valparaiso University, Indiana, USA.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Automatic pdb calling has been turned OFF\n"
     ]
    }
   ],
   "source": [
    "import pdb\n",
    "import itertools\n",
    "import sys\n",
    "import os\n",
    "import multiprocessing\n",
    "import random\n",
    "import copy\n",
    "import pickle\n",
    "\n",
    "import sklearn\n",
    "import sklearn.preprocessing\n",
    "import sklearn.decomposition\n",
    "from scipy.stats import gmean as geomean\n",
    "import matplotlib as M\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "np.set_printoptions(precision=3,suppress=True)\n",
    "pd.set_option(\"display.precision\",4)\n",
    "pd.set_option(\"display.float_format\",'{:.3f}'.format)\n",
    "np.random.seed(27)\n",
    "ncpus = 3\n",
    "%matplotlib inline\n",
    "%pdb off\n",
    "\n",
    "# Shortcut for fonts\n",
    "ff = {'fontname':\"Ubuntu\"}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The problem\n",
    "We desire to reduce uncertainty about future events as much as possible. For weather forecasts and related scientific communication to end-users, we aim to minimise forecast error, or maximise our information gain, for predicting hazards. The information required to determine an outcome is determined by the frequency of its occurrence, or its surprisal: \n",
    "\n",
    "$$\n",
    "S = -\\log_2(f)\n",
    "$$\n",
    "\n",
    "where $f$ is the probability of an event occurring. Herein, we are concerned with a binary probabilistic forecast, where we will mask a time series past a percentile threshold to create a set of Trues and Falses.  Hence with two categories $n=2$, we form an array of probabilities $\\mathbf{f}$ where each sample is denoted by $\\alpha \\in \\{1,2,...,N_\\alpha \\}$. First we consider only one sample, as the following scores are __scoring rules__ and can be applied to any number of samples $N_\\alpha \\in \\mathbb{Z}^+$. The two categories are $\\mathbf{f} = (1-f,f)$, with each element denoted by $i \\in \\{1,2,...,n\\}$; then average surprise provides a measure of uncertainty in $\\mathbf{f}$ when we summate the possibilities with their probabilities:\n",
    "\n",
    "$$\n",
    "H(f) = \\underbrace{(1-\\mathbf{f}) \\log_2(1-\\mathbf{f})}_{\\textrm{Event did not occur}} + \\underbrace{\\mathbf{f} \\log_2(\\mathbf{f})}_{\\textrm{Event occurred}} \\\\\n",
    "= \\sum^{n}_{i=1} \\mathbf{f}_i \\log_2 \\mathbf{f}_i \\\\\n",
    "= \\mathsf{E}[S(\\mathbf{f})]\n",
    "$$\n",
    "\n",
    "where the expectation operator $\\mathsf{E}[~]$ is shorthand for the mean outcome in the limit of time. Above, Shannon entropy $H$, here in bits due to the base-2 logarithm, represents the uncertainty in a time series that must be removed. For comparison, flipping a fair coin resolves *one bit* of uncertainty. Notice that a low value of $H$ but high $S$ occurs when focussing on a given rarer event, because the assumption of \"no event\" is successful more times than not. In meteorology, we can use this score to estimate how much uncertainty is removed by a given forecast by differencing the uncertainty before and after a forecast is obtained. This change in estimated likelihood is measured by the Kullback-Liebler Divergence $D_{\\textrm{KL}}$, which is the uncertainty or surprise (entropy) remaining after a forecast is issued, upon verification from observation $\\mathbf{o} = (1-o,o)$:\n",
    "\n",
    "$$\n",
    "D_{\\textrm{KL}}(\\mathbf{o}\\|\\mathbf{f}) = \\sum^{n}_{i=1} \\mathbf{o}_i \\log_2 \\frac{\\mathbf{o}_i}{\\mathbf{f}_i}\n",
    "$$\n",
    "\n",
    "For reference, for one event this can be expanded to\n",
    "\n",
    "$$          \n",
    "D_{\\textrm{KL}}(o\\|f) = (1-o) \\log2 \\frac{(1-o)}{(1-f)} + o \\log_2 \\frac{o}{f}\n",
    "$$\n",
    "\n",
    "Variations of this score occur under various names, including Ignorance (ref), the logarithm score (ref), etc. Herein, we implement $D_\\textrm{KL}$ within the Cross-Entropy Score (XES, ref), using cross-entropy $H_\\times$ which extends $D_\\textrm{KL}$ to include observational uncertainty. We also extend to 2-D forecast and observation arrays with $N_\\alpha > 1$:\n",
    "\n",
    "$$\n",
    "H_\\times(\\mathbf{o},\\mathbf{f}) = D_{\\textrm{KL}}(\\mathbf{o}\\|\\mathbf{f}) + H(\\mathbf{o}) \\\\ \\hspace{1cm}\n",
    "\\textrm{XES} \\equiv H_\\times (\\mathbf{o},\\mathbf{f}) \\hspace{2cm} \\\\\n",
    "= -\\frac{1}{N_\\alpha} \\sum^{N_\\alpha}_{\\alpha=1} \\sum^{n}_{i=1} \\mathbf{o}_{\\alpha,i} \\log_2 \\mathbf{f}_{\\alpha,i}\n",
    "$$\n",
    "\n",
    "Note that each forecast--observation pair $\\alpha$ (epochs; observations; timesteps) may be distinct from the real-world time $\\tau$. Often, observations are assumed to be perfect: $o \\in \\{0,1\\}$. We will assume so in the current paper; however, XES allows fractional probabilities to represent observational uncertainty. \n",
    "\n",
    "# Brier Score versus Cross-Entropy Score\n",
    "\n",
    "This reduces to \"Mean Square Error\" versus \"Cross-Entropy\" as measures of distances between the forecast and observation probability distributions. Perhaps due to the computational resources required to compute logarithms (Harold Brooks, pers. corr.), the mean squared error or __Brier Score__ is (ref) often used in meteorology:\n",
    "\n",
    "$$\n",
    "\\textrm{BS} = \\frac{1}{N_\\alpha} \\sum_\\alpha^{N_\\alpha}  (\\mathbf{f}-\\mathbf{o})^2\n",
    "$$\n",
    "\n",
    "\n",
    "Hence, the XES in one-term form (i.e., not DKL) for two events (say, rain or no-rain) is:\n",
    "\n",
    "$$\n",
    "\\textrm{XES} = o \\log_2 f - (1-o) \\log_2 (1-f) \\\\\n",
    "$$\n",
    "\n",
    "However, for forecasts of rare events (<10\\%), the Brier Score deviates substantially from the \"true\" measure of information gained. Benedetti (2010) showed that BS was a second-order approximation of XES (cf. above with their eqs. 27--31), where the biggest differences between the two scores occurred at extreme probabilities (near zero and unity). As in Benedetti, we can show this mathematically by first expanding the Brier Score for one event, with for forecast and observed probabilities $f$ and $o$ respectively.\n",
    "\n",
    "We demonstrate the scenarios where largest differences occur between BS and XES practically by comparing the scores, first normalising each score by a measure of uncertainty so we can compare like-for-like __skill scores__. (This is shown mathematically below the figure). Skill scores are preferred over the un-normalised score, as they are independent of the observational uncertainty $H(\\mathbf{o})$. We plot BSS and XESS for $f \\in [0.01,0.99]$ when the observed frequency $o \\in \\{0.1,0.3,0.5\\}$. __Forecast probabilities must be bounded to minimum and maximum probabilities__ $[ \\underline{p},\\bar{p} ]$ to avoid divergence to infinity, but also to avoid full confidence in the face of underdispersion (overfitting). We also set $\\log_2 0 = 0$. A positive value indicates better skill than climatology.\n",
    "\n",
    "Quote from Wikipedia regarding Benedetti paper: \"The Brier score becomes inadequate for very rare (or very frequent) events, because it does not sufficiently discriminate between small changes in forecast that are significant for rare events. Wilks (2010) has found that \"...quite large sample sizes, i.e. n > 1000, are required for higher-skill forecasts of relatively rare events, whereas only quite modest sample sizes are needed for low-skill forecasts of common events.\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BrierScore:\n",
    "    def __init__(self,f,o,fk=None,ob_uncertainty=True):\n",
    "        self.o = np.array(o)\n",
    "        self.f = np.array(f)\n",
    "\n",
    "        if fk is None:\n",
    "            self.fk = np.unique(self.f)\n",
    "        # NOT IMPLEMENTED\n",
    "        else:\n",
    "            pdb.set_trace()\n",
    "            # for k in np.unique(self.f):\n",
    "            #     assert k in fk\n",
    "            # self.fk = np.array(fk)\n",
    "\n",
    "        self.ob_uncertainty = ob_uncertainty\n",
    "\n",
    "    @staticmethod\n",
    "    def do_mean(x):\n",
    "        if isinstance(x,(float,int)):\n",
    "            return x\n",
    "        else:\n",
    "            return np.mean(x)\n",
    "\n",
    "    def compute_BS(self,from_components=False,return_all=False):\n",
    "        if from_components:\n",
    "            REL = self.compute_REL()\n",
    "            DSC = self.compute_DSC()\n",
    "            UNC = self.compute_UNC()\n",
    "            return REL - DSC + UNC\n",
    "        else:\n",
    "            bs_all = (self.f-self.o)**2\n",
    "            bs = np.sum(bs_all)/self.o.size\n",
    "            if return_all:\n",
    "                return bs, bs_all\n",
    "            return bs\n",
    "\n",
    "    def compute_UNC(self):\n",
    "        switch = 2\n",
    "        # if (len(self.o) > 1) and self.ob_uncertainty:\n",
    "        if self.ob_uncertainty and (switch == 1):\n",
    "            print(\"Compute UNC with ob uncertainty.\")\n",
    "            UNC = np.sum((self.o-self.do_mean(self.o))**2)/self.o.size\n",
    "        else:\n",
    "            print(\"Compute UNC without ob uncertainty\")\n",
    "            UNC = self.do_mean(self.o) * (1-self.do_mean(self.o))\n",
    "        print(f\"UNC = {UNC:.4f} (with {self.ob_uncertainty=})\")\n",
    "        return UNC\n",
    "\n",
    "    def compute_REL(self):\n",
    "        \"\"\" Compute reliability.\n",
    "        \"\"\"\n",
    "        rel = np.zeros_like(self.fk)\n",
    "        for nk,k in enumerate(self.fk):\n",
    "            ok = self.o[self.f==k]\n",
    "            if ok.size != 0:\n",
    "                ok_bar = self.do_mean(ok)\n",
    "                rel[nk] = ok.size * ((k-ok_bar)**2)\n",
    "        REL = np.sum(rel)/self.o.size\n",
    "        print(f\"{rel=}, REL = {REL:.4f}\")\n",
    "        return REL\n",
    "\n",
    "    def compute_DSC(self):\n",
    "        dsc = np.zeros_like(self.fk)\n",
    "        o_bar = self.do_mean(self.o)\n",
    "        for nk,k in enumerate(self.fk):\n",
    "            ok = self.o[self.f==k]\n",
    "            if ok.size != 0:\n",
    "                ok_bar = self.do_mean(ok)\n",
    "                dsc[nk] = ok.size * ((ok_bar-o_bar)**2)\n",
    "        DSC = np.sum(dsc)/self.o.size\n",
    "        print(f\"{dsc=}, DSC = {DSC:.4f}\")\n",
    "        return DSC\n",
    "\n",
    "    def compute_BSS(self,from_components=False,return_all=False):\n",
    "        BS_UNC = self.compute_UNC()\n",
    "        if return_all:\n",
    "            BS, BS_all = self.compute_BS(from_components=from_components,\n",
    "                        return_all=True)\n",
    "        else:\n",
    "            BS = self.compute_BS(from_components=from_components)\n",
    "        if return_all:\n",
    "            BSS_all = (BS_all-BS_UNC)/(0-BS_UNC)\n",
    "            BSS = np.mean(BSS_all)\n",
    "        else:\n",
    "            BSS = (BS-BS_UNC)/(0-BS_UNC)\n",
    "            \n",
    "        if return_all:\n",
    "            return BSS, BSS_all\n",
    "        return BSS\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "class CrossEntropy:\n",
    "    \"\"\" For computing and visualising skill scores for a probabilistic forecast.\n",
    "    Currently for two categories (binary) only, but with observational error.\n",
    "    \"\"\"\n",
    "    def __init__(self,f,o,fk=None):\n",
    "        \"\"\" Initialise the suite.\n",
    "        Args:\n",
    "        f:  1-D prob forecasts [0,1]. The probabilities should be bounded to\n",
    "                  avoid divergence to infinity\n",
    "        o:  1-D obs forecast [0,1], may or may not include obs uncertainty\n",
    "        fk: 1-D array of probability levels. If None, compute automatically.\n",
    "        \"\"\"\n",
    "        self.f = np.array(f)\n",
    "        self.o = np.array(o)\n",
    "        self.logbase = 2\n",
    "        # self.use_mean = \"geometric\"\n",
    "        self.use_mean = \"arithmetic\"\n",
    "\n",
    "        if fk == None:\n",
    "            self.fk = np.unique(self.f)\n",
    "        else:\n",
    "            raise Exception\n",
    "            # NOT IMPLEMENTED\n",
    "            # self.fk = fk\n",
    "\n",
    "    def do_mean(self,n):\n",
    "        if self.use_mean == \"geometric\":\n",
    "            f = geomean\n",
    "        elif self.use_mean == \"arithmetic\":\n",
    "            f = np.mean\n",
    "\n",
    "        if isinstance(n,(float,int)):\n",
    "            return n\n",
    "        else:\n",
    "            return f(n)\n",
    "\n",
    "    @staticmethod\n",
    "    def compute_DKL(x,y,return_all=False):\n",
    "        \"\"\" Kullback-Liebler Divergence\n",
    "\n",
    "        Args:\n",
    "        x   : 1-D (e.g., observations)\n",
    "        y   : 1-D (e.g., forecasts)\n",
    "        \"\"\"\n",
    "        if not isinstance(x,np.ndarray):\n",
    "            x = np.array([x,])\n",
    "        if not isinstance(y,np.ndarray):\n",
    "            y = np.array([y,])\n",
    "        with np.errstate(divide='ignore',invalid='ignore'):\n",
    "            # Enforce array so we can remove nans\n",
    "            term1 = (1-x) * np.log2((1-x)/(1-y))\n",
    "            term2 = x * np.log2(x/y)\n",
    "\n",
    "        term1[x==1] = 0\n",
    "        term2[x==0] = 0\n",
    "\n",
    "        DKL_all = term1 + term2\n",
    "        DKL = np.mean(DKL_all)\n",
    "        # print(f\"DKL = {DKL:.4f}\")\n",
    "\n",
    "        if return_all:\n",
    "            return DKL, DKL_all\n",
    "        return DKL\n",
    "\n",
    "    def compute_XESS(self,from_components=False,return_all=False):\n",
    "        UNC = self.compute_UNC()\n",
    "        XES = self.compute_XES(from_components=from_components,\n",
    "                                    return_all=return_all)\n",
    "        XESS = (XES-UNC)/(0-UNC)\n",
    "        return XESS\n",
    "\n",
    "    def compute_entropy(self,x,return_all=False):\n",
    "        \"\"\" Average surprise. A series of surprise\n",
    "        values are returned if return_all is True.\n",
    "        \"\"\"\n",
    "        with np.errstate(divide='ignore',invalid='ignore'):\n",
    "            H_all = (-(1-x) * np.log2(1-x)) - (x * np.log2(x))\n",
    "        H = self.do_mean(H_all)\n",
    "        if return_all:\n",
    "            return H, H_all\n",
    "        return H\n",
    "\n",
    "    def compute_XES(self,from_components=False,return_all=False,use_dkl=False):\n",
    "        \"\"\" Cross-entropy score.\n",
    "\n",
    "        Args:\n",
    "\n",
    "        \"\"\"\n",
    "        if from_components:\n",
    "            REL = self.compute_REL()\n",
    "            DSC = self.compute_DSC()\n",
    "            UNC = self.compute_UNC()\n",
    "            XES = REL-DSC+UNC\n",
    "        elif use_dkl:\n",
    "            if return_all:\n",
    "                XES, XES_all = self.compute_DKL(self.o,self.f,return_all=True)\n",
    "            else:\n",
    "                XES = self.compute_DKL(self.o,self.f)\n",
    "        else:\n",
    "            with np.errstate(divide='ignore',invalid='ignore'):\n",
    "                XES_all = (-((1-self.o)*np.log2(1-self.f))-\n",
    "                        (self.o*np.log2(self.f)))\n",
    "                XES = np.mean(XES_all)\n",
    "\n",
    "        if return_all:\n",
    "            return XES, XES_all\n",
    "        return XES\n",
    "\n",
    "    def compute_REL(self):\n",
    "        rel = np.zeros_like(self.fk)\n",
    "        for nk,k in enumerate(self.fk):\n",
    "            ok = self.o[self.f==k]\n",
    "            if ok.size != 0:\n",
    "                ok_hat = self.do_mean(ok)\n",
    "                rel[nk] = ok.size * self.compute_DKL(ok_hat,k)\n",
    "        REL = np.sum(rel)/self.o.size\n",
    "        print(f\"{rel=}, REL = {REL:.4f}\")\n",
    "        return REL\n",
    "\n",
    "    def compute_DSC(self):\n",
    "        dsc = np.zeros_like(self.fk)\n",
    "        for nk,k in enumerate(self.fk):\n",
    "            ok = self.o[self.f==k]\n",
    "            if ok.size != 0:\n",
    "                ok_hat = self.do_mean(ok)\n",
    "                o_hat = self.do_mean(self.o)\n",
    "                dsc[nk] = ok.size * self.compute_DKL(ok_hat,o_hat)\n",
    "        DSC = np.sum(dsc)/self.o.size\n",
    "        print(f\"{dsc=}, DSC = {DSC:.4f}\")\n",
    "        return DSC\n",
    "    \n",
    "    def compute_UNC(self):\n",
    "        p_o = self.do_mean(self.o)\n",
    "        UNC = self.compute_entropy(p_o)\n",
    "        print(f\"UNC = {UNC:.4f}\")\n",
    "        return UNC\n",
    "\n",
    "    @staticmethod\n",
    "    def bound(x,thresh=\"auto\",Ne=None):\n",
    "        if (thresh == \"auto\"):\n",
    "            thresh = 1/(3*Ne)\n",
    "        else:\n",
    "            assert isinstance(thresh,(float,np.float))\n",
    "        x[x<thresh] = thresh\n",
    "        x[x>(1-thresh)] = 1-thresh\n",
    "        return x\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SINGLE:       BS = 0.3600 ====================\n",
      "rel=array([0.36]), REL = 0.3600\n",
      "dsc=array([0.]), DSC = 0.0000\n",
      "Compute UNC without ob uncertainty\n",
      "UNC = 0.0000 (with self.ob_uncertainty=False)\n",
      "BS from comp = 0.3600\n",
      "\n",
      "Compute UNC without ob uncertainty\n",
      "UNC = 0.0000 (with self.ob_uncertainty=False)\n",
      "BSS = inf\n",
      "\n",
      "SINGLE:       XES = 1.3219 ====================\n",
      "rel=array([1.322]), REL = 1.3219\n",
      "dsc=array([0.]), DSC = 0.0000\n",
      "UNC = nan\n",
      "XES from comp= nan\n",
      "\n",
      "UNC = nan\n",
      "XESS = nan\n",
      "\n",
      "RAW:       BS = 0.2056 ====================\n",
      "rel=array([1.  , 0.02, 0.04, 0.25, 0.49, 0.04, 0.01, 0.  ]), REL = 0.2056\n",
      "dsc=array([0.309, 0.395, 0.198, 0.198, 0.198, 0.309, 0.309, 0.309]), DSC = 0.2469\n",
      "Compute UNC without ob uncertainty\n",
      "UNC = 0.2469 (with self.ob_uncertainty=False)\n",
      "BS from comp = 0.2056\n",
      "\n",
      "Compute UNC without ob uncertainty\n",
      "UNC = 0.2469 (with self.ob_uncertainty=False)\n",
      "BSS = 0.1675\n",
      "\n",
      "RAW:       XES = nan ====================\n",
      "rel=array([  inf, 0.304, 0.322, 1.   , 1.737, 0.322, 0.152, 0.   ]), REL = inf\n",
      "dsc=array([1.17 , 1.696, 0.848, 0.848, 0.848, 1.17 , 1.17 , 1.17 ]), DSC = 0.9911\n",
      "UNC = 0.9911\n",
      "XES from comp= inf\n",
      "\n",
      "UNC = 0.9911\n",
      "XESS = nan\n",
      "\n",
      "PROBS 1pc:       BS = 0.2034 ====================\n",
      "rel=array([0.98, 0.02, 0.04, 0.25, 0.49, 0.04, 0.01, 0.  ]), REL = 0.2034\n",
      "dsc=array([0.309, 0.395, 0.198, 0.198, 0.198, 0.309, 0.309, 0.309]), DSC = 0.2469\n",
      "Compute UNC without ob uncertainty\n",
      "UNC = 0.2469 (with self.ob_uncertainty=False)\n",
      "BS from comp = 0.2034\n",
      "\n",
      "Compute UNC without ob uncertainty\n",
      "UNC = 0.2469 (with self.ob_uncertainty=False)\n",
      "BSS = 0.1764\n",
      "\n",
      "PROBS 1pc:       XES = 1.1661 ====================\n",
      "rel=array([6.644, 0.304, 0.322, 1.   , 1.737, 0.322, 0.152, 0.014]), REL = 1.1661\n",
      "dsc=array([1.17 , 1.696, 0.848, 0.848, 0.848, 1.17 , 1.17 , 1.17 ]), DSC = 0.9911\n",
      "UNC = 0.9911\n",
      "XES from comp= 1.1661\n",
      "\n",
      "UNC = 0.9911\n",
      "XESS = -0.1766\n",
      "\n",
      "BOTH 1pc FALSE:       BS = 0.1970 ====================\n",
      "rel=array([0.96 , 0.016, 0.036, 0.24 , 0.476, 0.036, 0.008, 0.   ]), REL = 0.1970\n",
      "dsc=array([0.296, 0.379, 0.19 , 0.19 , 0.19 , 0.296, 0.296, 0.296]), DSC = 0.2371\n",
      "Compute UNC without ob uncertainty\n",
      "UNC = 0.2470 (with self.ob_uncertainty=False)\n",
      "BS from comp = 0.2069\n",
      "\n",
      "Compute UNC without ob uncertainty\n",
      "UNC = 0.2470 (with self.ob_uncertainty=False)\n",
      "BSS = 0.2025\n",
      "\n",
      "BOTH 1pc FALSE:       XES = 1.1798 ====================\n",
      "rel=array([6.497, 0.206, 0.261, 0.919, 1.644, 0.261, 0.103, 0.   ]), REL = 1.0990\n",
      "dsc=array([1.082, 1.546, 0.773, 0.773, 0.773, 1.082, 1.082, 1.082]), DSC = 0.9106\n",
      "UNC = 0.9914\n",
      "XES from comp= 1.1798\n",
      "\n",
      "UNC = 0.9914\n",
      "XESS = -0.1900\n",
      "\n",
      "OBS 0.1pc:       BS = 0.2027 ====================\n",
      "rel=array([0.978, 0.02 , 0.04 , 0.249, 0.489, 0.04 , 0.01 , 0.   ]), REL = 0.2027\n",
      "dsc=array([0.307, 0.393, 0.197, 0.197, 0.197, 0.307, 0.307, 0.307]), DSC = 0.2459\n",
      "Compute UNC without ob uncertainty\n",
      "UNC = 0.2469 (with self.ob_uncertainty=True)\n",
      "BS from comp = 0.2037\n",
      "\n",
      "Compute UNC without ob uncertainty\n",
      "UNC = 0.2469 (with self.ob_uncertainty=True)\n",
      "BSS = 0.1791\n",
      "\n",
      "OBS 0.1pc:       XES = 1.1675 ====================\n",
      "rel=array([6.626, 0.288, 0.313, 0.989, 1.724, 0.313, 0.144, 0.01 ]), REL = 1.1561\n",
      "dsc=array([1.158, 1.674, 0.837, 0.837, 0.837, 1.158, 1.158, 1.158]), DSC = 0.9797\n",
      "UNC = 0.9911\n",
      "XES from comp= 1.1675\n",
      "\n",
      "UNC = 0.9911\n",
      "XESS = -0.1780\n",
      "\n",
      "IGNORE NULL       BS = 0.2031 ====================\n",
      "rel=array([0.978, 0.02 , 0.04 , 0.25 , 0.49 , 0.04 , 0.01 , 0.   ]), REL = 0.2031\n",
      "dsc=array([0.308, 0.394, 0.197, 0.197, 0.197, 0.308, 0.308, 0.308]), DSC = 0.2464\n",
      "Compute UNC without ob uncertainty\n",
      "UNC = 0.2469 (with self.ob_uncertainty=True)\n",
      "BS from comp = 0.2035\n",
      "\n",
      "Compute UNC without ob uncertainty\n",
      "UNC = 0.2469 (with self.ob_uncertainty=True)\n",
      "BSS = 0.1774\n",
      "\n",
      "IGNORE NULL       XES = 1.1667 ====================\n",
      "rel=array([6.626, 0.304, 0.322, 1.   , 1.737, 0.313, 0.144, 0.01 ]), REL = 1.1616\n",
      "dsc=array([1.16 , 1.694, 0.847, 0.847, 0.847, 1.16 , 1.16 , 1.16 ]), DSC = 0.9859\n",
      "UNC = 0.9909\n",
      "XES from comp= 1.1667\n",
      "\n",
      "UNC = 0.9909\n",
      "XESS = -0.1774\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-5-2cb4a79a9e3c>:85: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  BSS = (BS-BS_UNC)/(0-BS_UNC)\n"
     ]
    }
   ],
   "source": [
    "do_brier = True\n",
    "do_xes = True\n",
    "def print_tests(f,o,title,ob_unc):\n",
    "    if do_brier:\n",
    "        brier = BrierScore(f,o,ob_uncertainty=ob_unc)\n",
    "        print(title,(\" \"*5),f\"BS = {brier.compute_BS():.4f}\",(\"=\"*20),)\n",
    "        print(f\"BS from comp = {brier.compute_BS(from_components=True):.4f}\\n\")\n",
    "        print(f\"BSS = {brier.compute_BSS():.4f}\\n\")\n",
    "    if do_xes:\n",
    "        xes = CrossEntropy(f,o)\n",
    "        print(title,(\" \"*5),\n",
    "            f\"XES = {xes.compute_XES(from_components=False):.4f}\",\n",
    "                        (\"=\"*20),)\n",
    "        print(f\"XES from comp= {xes.compute_XES(from_components=True):.4f}\\n\")\n",
    "        print(f\"XESS = {xes.compute_XESS():.4f}\\n\")\n",
    "    return\n",
    "\n",
    "\n",
    "# Testing for one forecast\n",
    "\n",
    "# Probs are unbounded\n",
    "# Obs are binary\n",
    "f = [0.4,]\n",
    "o = [1,]\n",
    "print_tests(f,o,\"SINGLE:\",False)\n",
    "\n",
    "###########################################################\n",
    "\n",
    "# Testing for multiple forecasts\n",
    "\n",
    "# Probs are unbounded\n",
    "# Obs are binary\n",
    "f = [0.1,0.5,0.9,0.8,0.7,0.2,1,0,0.1]\n",
    "o = [0,0,1,1,0,0,1,1,0]\n",
    "print_tests(f,o,\"RAW:\",False)\n",
    "\n",
    "# Probs are bounded to [0.01,0.99]\n",
    "# Obs are binary\n",
    "f = [0.1,0.5,0.9,0.8,0.7,0.2,0.99,0.01,0.1]\n",
    "o = [0,0,1,1,0,0,1,1,0]\n",
    "print_tests(f,o,\"PROBS 1pc:\",False)\n",
    "\n",
    "# Probs are bounded to [0.01,0.99]\n",
    "# Obs are from {0.01,0.99} (same certainty)\n",
    "f = [0.1,0.5,0.9,0.8,0.7,0.2,0.99,0.01,0.1]\n",
    "o = [0.01,0.01,0.99,0.99,0.01,0.01,0.99,0.99,0.01]\n",
    "print_tests(f,o,\"BOTH 1pc FALSE:\",False)\n",
    "\n",
    "# Probs are bounded to [0.01,0.99]\n",
    "# Obs are from {0.001,0.999} (high certainty)\n",
    "f = [0.1,0.5,0.9,0.8,0.7,0.2,0.99,0.01,0.1]\n",
    "o = [0.001,0.001,0.999,0.999,0.001,0.001,0.999,0.999,0.001]\n",
    "print_tests(f,o,\"OBS 0.1pc:\",True)\n",
    "\n",
    "# Testing for multiple forecasts\n",
    "# Probs are bounded to [0.01,0.99]\n",
    "# Obs are from {0,0.999} (ignore null forecast error)\n",
    "f = [0.1,0.5,0.9,0.8,0.7,0.2,0.99,0.01,0.1]\n",
    "o = [0,0,0.999,0.999,0,0,0.999,0.999,0]\n",
    "print_tests(f,o,\"IGNORE NULL\",True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compute UNC without ob uncertainty\n",
      "UNC = 0.0900 (with self.ob_uncertainty=True)\n",
      "Compute UNC without ob uncertainty\n",
      "UNC = 0.0900 (with self.ob_uncertainty=True)\n",
      "-1.6631018518518523\n",
      "[ 0.9    0.91   0.92   0.929  0.937  0.946  0.953  0.96   0.966  0.972\n",
      "  0.978  0.982  0.986  0.99   0.993  0.996  0.997  0.999  1.     1.\n",
      "  1.     0.999  0.997  0.996  0.993  0.99   0.986  0.982  0.978  0.972\n",
      "  0.966  0.96   0.953  0.946  0.937  0.929  0.92   0.91   0.9    0.889\n",
      "  0.877  0.866  0.853  0.84   0.826  0.812  0.797  0.782  0.766  0.75\n",
      "  0.733  0.716  0.698  0.679  0.66   0.64   0.62   0.599  0.578  0.556\n",
      "  0.533  0.51   0.486  0.462  0.437  0.412  0.386  0.36   0.333  0.306\n",
      "  0.277  0.249  0.22   0.19   0.16   0.129  0.097  0.066  0.033 -0.\n",
      " -0.034 -0.068 -0.103 -0.138 -0.174 -0.21  -0.247 -0.284 -0.322 -0.361\n",
      " -0.4   -0.44  -0.48  -0.521 -0.563 -0.604 -0.647 -0.69  -0.734 -0.778\n",
      " -0.823 -0.868 -0.914 -0.96  -1.007 -1.054 -1.103 -1.151 -1.2   -1.25\n",
      " -1.3   -1.351 -1.403 -1.454 -1.507 -1.56  -1.614 -1.668 -1.722 -1.778\n",
      " -1.834 -1.89  -1.947 -2.004 -2.062 -2.121 -2.18  -2.24  -2.3   -2.361\n",
      " -2.423 -2.484 -2.547 -2.61  -2.674 -2.738 -2.803 -2.868 -2.934 -3.\n",
      " -3.067 -3.134 -3.202 -3.271 -3.34  -3.41  -3.48  -3.551 -3.623 -3.694\n",
      " -3.767 -3.84  -3.914 -3.988 -4.063 -4.138 -4.214 -4.29  -4.367 -4.444\n",
      " -4.523 -4.601 -4.68  -4.76  -4.84  -4.921 -5.003 -5.084 -5.167 -5.25\n",
      " -5.334 -5.418 -5.503 -5.588 -5.674 -5.76  -5.847 -5.934 -6.023 -6.111\n",
      " -6.2   -6.29  -6.38  -6.471 -6.563 -6.654 -6.747 -6.84  -6.934 -7.028\n",
      " -7.123 -7.218 -7.314 -7.41  -7.507 -7.604 -7.703 -7.801]\n",
      "UNC = 0.4690\n",
      "UNC = 0.4690\n",
      "-1.9776780071582463\n",
      "[ -0.644  -0.444  -0.334  -0.259  -0.205  -0.163  -0.13   -0.103  -0.081\n",
      "  -0.064  -0.049  -0.037  -0.027  -0.019  -0.013  -0.008  -0.004  -0.002\n",
      "  -0.     -0.     -0.     -0.002  -0.004  -0.006  -0.009  -0.013  -0.017\n",
      "  -0.022  -0.028  -0.034  -0.04   -0.046  -0.053  -0.061  -0.069  -0.077\n",
      "  -0.085  -0.094  -0.103  -0.113  -0.123  -0.133  -0.143  -0.154  -0.165\n",
      "  -0.176  -0.187  -0.199  -0.211  -0.223  -0.235  -0.248  -0.261  -0.274\n",
      "  -0.287  -0.301  -0.315  -0.329  -0.343  -0.358  -0.373  -0.388  -0.403\n",
      "  -0.418  -0.434  -0.45   -0.466  -0.482  -0.499  -0.516  -0.533  -0.55\n",
      "  -0.567  -0.585  -0.603  -0.621  -0.639  -0.658  -0.677  -0.696  -0.715\n",
      "  -0.735  -0.755  -0.775  -0.795  -0.816  -0.837  -0.858  -0.879  -0.901\n",
      "  -0.923  -0.945  -0.967  -0.99   -1.013  -1.036  -1.06   -1.084  -1.108\n",
      "  -1.132  -1.157  -1.182  -1.207  -1.233  -1.259  -1.286  -1.312  -1.339\n",
      "  -1.367  -1.395  -1.423  -1.451  -1.48   -1.509  -1.539  -1.569  -1.6\n",
      "  -1.631  -1.662  -1.694  -1.726  -1.759  -1.792  -1.826  -1.86   -1.895\n",
      "  -1.93   -1.966  -2.002  -2.039  -2.076  -2.115  -2.153  -2.193  -2.233\n",
      "  -2.273  -2.315  -2.357  -2.399  -2.443  -2.487  -2.532  -2.578  -2.625\n",
      "  -2.673  -2.722  -2.771  -2.822  -2.874  -2.926  -2.98   -3.035  -3.092\n",
      "  -3.149  -3.208  -3.268  -3.33   -3.393  -3.458  -3.524  -3.593  -3.663\n",
      "  -3.735  -3.809  -3.885  -3.963  -4.044  -4.127  -4.213  -4.302  -4.394\n",
      "  -4.49   -4.589  -4.691  -4.798  -4.909  -5.025  -5.147  -5.274  -5.407\n",
      "  -5.547  -5.695  -5.852  -6.018  -6.195  -6.385  -6.588  -6.808  -7.047\n",
      "  -7.31   -7.6    -7.924  -8.292  -8.717  -9.221  -9.837 -10.632 -11.753]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-6-8a770ebb756c>:69: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.\n",
      "  XESS = (XES-UNC)/(0-UNC)\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "x and y must have same first dimension, but have shapes (198,) and (1,)",
     "output_type": "error",
     "traceback": [
      "\u001B[0;31m---------------------------------------------------------------------------\u001B[0m",
      "\u001B[0;31mValueError\u001B[0m                                Traceback (most recent call last)",
      "\u001B[0;32m<ipython-input-8-896996694f2c>\u001B[0m in \u001B[0;36m<module>\u001B[0;34m\u001B[0m\n\u001B[1;32m     45\u001B[0m     \u001B[0mprint\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mXESS_all\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m     46\u001B[0m \u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0;32m---> 47\u001B[0;31m     \u001B[0mplot_BSSvXESS\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mf\u001B[0m\u001B[0;34m,\u001B[0m\u001B[0mBSS\u001B[0m\u001B[0;34m,\u001B[0m\u001B[0mXESS\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0m\u001B[1;32m     48\u001B[0m \u001B[0;34m\u001B[0m\u001B[0m\n",
      "\u001B[0;32m<ipython-input-8-896996694f2c>\u001B[0m in \u001B[0;36mplot_BSSvXESS\u001B[0;34m(f, BSS, XESS)\u001B[0m\n\u001B[1;32m      2\u001B[0m     \u001B[0mfig\u001B[0m\u001B[0;34m,\u001B[0m\u001B[0max\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0mplt\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0msubplots\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0;36m1\u001B[0m\u001B[0;34m,\u001B[0m\u001B[0mfigsize\u001B[0m\u001B[0;34m=\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0;36m7\u001B[0m\u001B[0;34m,\u001B[0m\u001B[0;36m4\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[1;32m      3\u001B[0m \u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0;32m----> 4\u001B[0;31m     \u001B[0max\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mplot\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mf\u001B[0m\u001B[0;34m,\u001B[0m\u001B[0mBSS\u001B[0m\u001B[0;34m,\u001B[0m\u001B[0mlabel\u001B[0m\u001B[0;34m=\u001B[0m\u001B[0;34m'BSS'\u001B[0m\u001B[0;34m,\u001B[0m\u001B[0mzorder\u001B[0m\u001B[0;34m=\u001B[0m\u001B[0;36m3\u001B[0m\u001B[0;34m,\u001B[0m\u001B[0mcolor\u001B[0m\u001B[0;34m=\u001B[0m\u001B[0;34m\"red\"\u001B[0m\u001B[0;34m,\u001B[0m\u001B[0mlinestyle\u001B[0m\u001B[0;34m=\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[0m\u001B[1;32m      5\u001B[0m     \u001B[0max\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mplot\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mf\u001B[0m\u001B[0;34m,\u001B[0m\u001B[0mXESS\u001B[0m\u001B[0;34m,\u001B[0m\u001B[0mlabel\u001B[0m\u001B[0;34m=\u001B[0m\u001B[0;34m'XESS'\u001B[0m\u001B[0;34m,\u001B[0m\u001B[0mzorder\u001B[0m\u001B[0;34m=\u001B[0m\u001B[0;36m3\u001B[0m\u001B[0;34m,\u001B[0m\u001B[0mcolor\u001B[0m\u001B[0;34m=\u001B[0m\u001B[0;34m\"green\"\u001B[0m\u001B[0;34m,\u001B[0m\u001B[0mlinestyle\u001B[0m\u001B[0;34m=\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[1;32m      6\u001B[0m \u001B[0;34m\u001B[0m\u001B[0m\n",
      "\u001B[0;32m~/anaconda3/envs/research/lib/python3.8/site-packages/matplotlib/axes/_axes.py\u001B[0m in \u001B[0;36mplot\u001B[0;34m(self, scalex, scaley, data, *args, **kwargs)\u001B[0m\n\u001B[1;32m   1741\u001B[0m         \"\"\"\n\u001B[1;32m   1742\u001B[0m         \u001B[0mkwargs\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0mcbook\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mnormalize_kwargs\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mkwargs\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mmlines\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mLine2D\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[0;32m-> 1743\u001B[0;31m         \u001B[0mlines\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0;34m[\u001B[0m\u001B[0;34m*\u001B[0m\u001B[0mself\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0m_get_lines\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0;34m*\u001B[0m\u001B[0margs\u001B[0m\u001B[0;34m,\u001B[0m \u001B[0mdata\u001B[0m\u001B[0;34m=\u001B[0m\u001B[0mdata\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[0m\u001B[1;32m   1744\u001B[0m         \u001B[0;32mfor\u001B[0m \u001B[0mline\u001B[0m \u001B[0;32min\u001B[0m \u001B[0mlines\u001B[0m\u001B[0;34m:\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m   1745\u001B[0m             \u001B[0mself\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0madd_line\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mline\u001B[0m\u001B[0;34m)\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n",
      "\u001B[0;32m~/anaconda3/envs/research/lib/python3.8/site-packages/matplotlib/axes/_base.py\u001B[0m in \u001B[0;36m__call__\u001B[0;34m(self, data, *args, **kwargs)\u001B[0m\n\u001B[1;32m    271\u001B[0m                 \u001B[0mthis\u001B[0m \u001B[0;34m+=\u001B[0m \u001B[0margs\u001B[0m\u001B[0;34m[\u001B[0m\u001B[0;36m0\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[1;32m    272\u001B[0m                 \u001B[0margs\u001B[0m \u001B[0;34m=\u001B[0m \u001B[0margs\u001B[0m\u001B[0;34m[\u001B[0m\u001B[0;36m1\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--> 273\u001B[0;31m             \u001B[0;32myield\u001B[0m \u001B[0;32mfrom\u001B[0m \u001B[0mself\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0m_plot_args\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mthis\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[0m\n\u001B[0m\u001B[1;32m    274\u001B[0m \u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m    275\u001B[0m     \u001B[0;32mdef\u001B[0m \u001B[0mget_next_color\u001B[0m\u001B[0;34m(\u001B[0m\u001B[0mself\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~/anaconda3/envs/research/lib/python3.8/site-packages/matplotlib/axes/_base.py\u001B[0m in \u001B[0;36m_plot_args\u001B[0;34m(self, tup, kwargs)\u001B[0m\n\u001B[1;32m    397\u001B[0m \u001B[0;34m\u001B[0m\u001B[0m\n\u001B[1;32m    398\u001B[0m         \u001B[0;32mif\u001B[0m \u001B[0mx\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mshape\u001B[0m\u001B[0;34m[\u001B[0m\u001B[0;36m0\u001B[0m\u001B[0;34m]\u001B[0m \u001B[0;34m!=\u001B[0m \u001B[0my\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mshape\u001B[0m\u001B[0;34m[\u001B[0m\u001B[0;36m0\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--> 399\u001B[0;31m             raise ValueError(f\"x and y must have same first dimension, but \"\n\u001B[0m\u001B[1;32m    400\u001B[0m                              f\"have shapes {x.shape} and {y.shape}\")\n\u001B[1;32m    401\u001B[0m         \u001B[0;32mif\u001B[0m \u001B[0mx\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mndim\u001B[0m \u001B[0;34m>\u001B[0m \u001B[0;36m2\u001B[0m \u001B[0;32mor\u001B[0m \u001B[0my\u001B[0m\u001B[0;34m.\u001B[0m\u001B[0mndim\u001B[0m \u001B[0;34m>\u001B[0m \u001B[0;36m2\u001B[0m\u001B[0;34m:\u001B[0m\u001B[0;34m\u001B[0m\u001B[0;34m\u001B[0m\u001B[0m\n",
      "\u001B[0;31mValueError\u001B[0m: x and y must have same first dimension, but have shapes (198,) and (1,)"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbQAAAD8CAYAAAAfSFHzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAANnUlEQVR4nO3dX4hc93mH8edbKYbmT2MTb0IqOUQtchxd2MXeOKE0rdPQRnIvRMAXdkJMTUCYxiGXNoUmF75pLgohxIkQRpjcRBeNSZTixBRK4oLrVCvwP8XYbGVqb2WwHIcUHKiR/fZipmWy3t05Gs+uVq+eDwzsmfOb3ZcfKx6d3dFRqgpJki52v3OhB5AkaR4MmiSpBYMmSWrBoEmSWjBokqQWDJokqYWpQUtyNMnLSZ5e53ySfDPJcpInk1w//zElSdrYkCu0B4D9G5w/AOwdPw4B33n7Y0mSdH6mBq2qHgFe3WDJQeC7NfIYcHmSD85rQEmShtg5h8+xC3hx4nhl/NxLqxcmOcToKo53vetdN1xzzTVz+PKSpC5Onjz5SlUtzPLaeQQtazy35v20quoIcARgcXGxlpaW5vDlJUldJPnPWV87j3c5rgBXTRzvBs7M4fNKkjTYPIJ2HLh9/G7HTwC/rqq3/LhRkqTNNPVHjkm+B9wEXJlkBfga8A6AqjoMPATcDCwDvwHu2KxhJUlaz9SgVdVtU84X8KW5TSRJ0gy8U4gkqQWDJklqwaBJklowaJKkFgyaJKkFgyZJasGgSZJaMGiSpBYMmiSpBYMmSWrBoEmSWjBokqQWDJokqQWDJklqwaBJklowaJKkFgyaJKkFgyZJasGgSZJaMGiSpBYMmiSpBYMmSWrBoEmSWjBokqQWDJokqQWDJklqwaBJklowaJKkFgyaJKkFgyZJasGgSZJaMGiSpBYMmiSpBYMmSWrBoEmSWhgUtCT7kzybZDnJPWucf2+SHyV5IsmpJHfMf1RJktY3NWhJdgD3AQeAfcBtSfatWvYl4BdVdR1wE/APSS6b86ySJK1ryBXajcByVZ2uqteBY8DBVWsKeE+SAO8GXgXOzXVSSZI2MCRou4AXJ45Xxs9N+hbwUeAM8BTwlap6c/UnSnIoyVKSpbNnz844siRJbzUkaFnjuVp1/BngceD3gT8CvpXk997yoqojVbVYVYsLCwvnOaokSesbErQV4KqJ492MrsQm3QE8WCPLwPPANfMZUZKk6YYE7QSwN8me8Rs9bgWOr1rzAvBpgCQfAD4CnJ7noJIkbWTntAVVdS7JXcDDwA7gaFWdSnLn+Pxh4F7ggSRPMfoR5d1V9comzi1J0m+ZGjSAqnoIeGjVc4cnPj4D/OV8R5MkaTjvFCJJasGgSZJaMGiSpBYMmiSpBYMmSWrBoEmSWjBokqQWDJokqQWDJklqwaBJklowaJKkFgyaJKkFgyZJasGgSZJaMGiSpBYMmiSpBYMmSWrBoEmSWjBokqQWDJokqQWDJklqwaBJklowaJKkFgyaJKkFgyZJasGgSZJaMGiSpBYMmiSpBYMmSWrBoEmSWjBokqQWDJokqQWDJklqwaBJkloYFLQk+5M8m2Q5yT3rrLkpyeNJTiX52XzHlCRpYzunLUiyA7gP+AtgBTiR5HhV/WJizeXAt4H9VfVCkvdv0rySJK1pyBXajcByVZ2uqteBY8DBVWs+BzxYVS8AVNXL8x1TkqSNDQnaLuDFieOV8XOTrgauSPLTJCeT3L7WJ0pyKMlSkqWzZ8/ONrEkSWsYErSs8VytOt4J3AD8FfAZ4O+SXP2WF1UdqarFqlpcWFg472ElSVrP1N+hMboiu2rieDdwZo01r1TVa8BrSR4BrgOem8uUkiRNMeQK7QSwN8meJJcBtwLHV635IfDJJDuTvBP4OPDMfEeVJGl9U6/QqupckruAh4EdwNGqOpXkzvH5w1X1TJKfAE8CbwL3V9XTmzm4JEmTUrX612FbY3FxsZaWli7I15YkbU9JTlbV4iyv9U4hkqQWDJokqQWDJklqwaBJklowaJKkFgyaJKkFgyZJasGgSZJaMGiSpBYMmiSpBYMmSWrBoEmSWjBokqQWDJokqQWDJklqwaBJklowaJKkFgyaJKkFgyZJasGgSZJaMGiSpBYMmiSpBYMmSWrBoEmSWjBokqQWDJokqQWDJklqwaBJklowaJKkFgyaJKkFgyZJasGgSZJaMGiSpBYMmiSpBYMmSWphUNCS7E/ybJLlJPdssO5jSd5Icsv8RpQkabqpQUuyA7gPOADsA25Lsm+ddV8HHp73kJIkTTPkCu1GYLmqTlfV68Ax4OAa674MfB94eY7zSZI0yJCg7QJenDheGT/3/5LsAj4LHN7oEyU5lGQpydLZs2fPd1ZJktY1JGhZ47ladfwN4O6qemOjT1RVR6pqsaoWFxYWBo4oSdJ0OwesWQGumjjeDZxZtWYROJYE4Erg5iTnquoH8xhSkqRphgTtBLA3yR7gv4Bbgc9NLqiqPf/3cZIHgH8yZpKkrTQ1aFV1LsldjN69uAM4WlWnktw5Pr/h780kSdoKQ67QqKqHgIdWPbdmyKrqr9/+WJIknR/vFCJJasGgSZJaMGiSpBYMmiSpBYMmSWrBoEmSWjBokqQWDJokqQWDJklqwaBJklowaJKkFgyaJKkFgyZJasGgSZJaMGiSpBYMmiSpBYMmSWrBoEmSWjBokqQWDJokqQWDJklqwaBJklowaJKkFgyaJKkFgyZJasGgSZJaMGiSpBYMmiSpBYMmSWrBoEmSWjBokqQWDJokqQWDJklqwaBJkloYFLQk+5M8m2Q5yT1rnP98kifHj0eTXDf/USVJWt/UoCXZAdwHHAD2Abcl2bdq2fPAn1XVtcC9wJF5DypJ0kaGXKHdCCxX1emqeh04BhycXFBVj1bVr8aHjwG75zumJEkbGxK0XcCLE8cr4+fW80Xgx2udSHIoyVKSpbNnzw6fUpKkKYYELWs8V2suTD7FKGh3r3W+qo5U1WJVLS4sLAyfUpKkKXYOWLMCXDVxvBs4s3pRkmuB+4EDVfXL+YwnSdIwQ67QTgB7k+xJchlwK3B8ckGSDwEPAl+oqufmP6YkSRubeoVWVeeS3AU8DOwAjlbVqSR3js8fBr4KvA/4dhKAc1W1uHljS5L021K15q/DNt3i4mItLS1dkK8tSdqekpyc9YLIO4VIklowaJKkFgyaJKkFgyZJasGgSZJaMGiSpBYMmiSpBYMmSWrBoEmSWjBokqQWDJokqQWDJklqwaBJklowaJKkFgyaJKkFgyZJasGgSZJaMGiSpBYMmiSpBYMmSWrBoEmSWjBokqQWDJokqQWDJklqwaBJklowaJKkFgyaJKkFgyZJasGgSZJaMGiSpBYMmiSpBYMmSWrBoEmSWjBokqQWDJokqYVBQUuyP8mzSZaT3LPG+ST55vj8k0mun/+okiStb2rQkuwA7gMOAPuA25LsW7XsALB3/DgEfGfOc0qStKEhV2g3AstVdbqqXgeOAQdXrTkIfLdGHgMuT/LBOc8qSdK6dg5Yswt4ceJ4Bfj4gDW7gJcmFyU5xOgKDuB/kjx9XtMK4ErglQs9xEXIfZudezcb9202H5n1hUOCljWeqxnWUFVHgCMASZaqanHA19cE92027tvs3LvZuG+zSbI062uH/MhxBbhq4ng3cGaGNZIkbZohQTsB7E2yJ8llwK3A8VVrjgO3j9/t+Ang11X10upPJEnSZpn6I8eqOpfkLuBhYAdwtKpOJblzfP4w8BBwM7AM/Aa4Y8DXPjLz1Jc292027tvs3LvZuG+zmXnfUvWWX3VJknTR8U4hkqQWDJokqYVND5q3zZrNgH37/Hi/nkzyaJLrLsSc2820fZtY97EkbyS5ZSvn266G7FuSm5I8nuRUkp9t9Yzb0YA/p+9N8qMkT4z3bcj7C9pLcjTJy+v9W+SZu1BVm/Zg9CaS/wD+ALgMeALYt2rNzcCPGf1btk8AP9/MmS6Gx8B9+2PgivHHB9y3Yfs2se5fGL2Z6ZYLPfeFfgz8frsc+AXwofHx+y/03Bf6MXDf/hb4+vjjBeBV4LILPfuFfgB/ClwPPL3O+Zm6sNlXaN42azZT962qHq2qX40PH2P0b/8udUO+3wC+DHwfeHkrh9vGhuzb54AHq+oFgKpy74btWwHvSRLg3YyCdm5rx9x+quoRRnuxnpm6sNlBW++WWOe75lJzvnvyRUZ/m7nUTd23JLuAzwKHt3Cu7W7I99vVwBVJfprkZJLbt2y67WvIvn0L+CijG008BXylqt7cmvEuajN1Ycitr96Oud026xIzeE+SfIpR0P5kUye6OAzZt28Ad1fVG6O/NIth+7YTuAH4NPC7wL8leayqntvs4baxIfv2GeBx4M+BPwT+Ocm/VtV/b/JsF7uZurDZQfO2WbMZtCdJrgXuBw5U1S+3aLbtbMi+LQLHxjG7Erg5ybmq+sGWTLg9Df1z+kpVvQa8luQR4DrgUg7akH27A/j7Gv1iaDnJ88A1wL9vzYgXrZm6sNk/cvS2WbOZum9JPgQ8CHzhEv9b8qSp+1ZVe6rqw1X1YeAfgb+5xGMGw/6c/hD4ZJKdSd7J6H/ceGaL59xuhuzbC4yuaknyAUZ3kj+9pVNenGbqwqZeodXm3TartYH79lXgfcC3x1cb5+oSv7P3wH3TKkP2raqeSfIT4EngTeD+qrqk//ungd9v9wIPJHmK0Y/R7q6qS/6/lEnyPeAm4MokK8DXgHfA2+uCt76SJLXgnUIkSS0YNElSCwZNktSCQZMktWDQJEktGDRJUgsGTZLUwv8CYt4dDJmp3HkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_BSSvXESS(f,BSS,XESS):\n",
    "    fig,ax = plt.subplots(1,figsize=(7,4))\n",
    "    \n",
    "    ax.plot(f,BSS,label='BSS',zorder=3,color=\"red\",linestyle='-')\n",
    "    ax.plot(f,XESS,label='XESS',zorder=3,color=\"green\",linestyle=\"-\")\n",
    "    \n",
    "    ax.legend()\n",
    "    ax.hlines(y=0,xmin=f.min(),xmax=f.max(),color='k',\n",
    "              zorder=2)\n",
    "    ax.vlines(x=o,ymin=-2,ymax=1,color='k',zorder=2)\n",
    "    ax.set_xlabel('Forecast probability',**ff)\n",
    "    ax.set_ylabel('Skill score',**ff)\n",
    "    ax.set_title(f\"Skill scores for observed event of {int(o*100):d}% frequency\",**ff)\n",
    "    ax.set_yticks([-1,0,1])\n",
    "    ax.set_xticks(np.arange(0.0,1.1,0.1))\n",
    "\n",
    "    ax.fill_between(x=[0,0.05],y1=[-2,-2],y2=[1,1],\n",
    "                    zorder=1,color='lime',alpha=0.2)\n",
    "    ax.fill_between(x=[0.95,1],y1=[-2,-2],y2=[1,1],\n",
    "                    zorder=1,color='lime',alpha=0.2)\n",
    "    display(fig)\n",
    "    plt.close(fig)\n",
    "\n",
    "f = np.arange(0.005,0.995,0.005)\n",
    "\n",
    "for o in (0.1,0.3,0.5):\n",
    "    o = np.array([o,]*len(f))\n",
    "    assert o.size == f.size\n",
    "    bs = BrierScore(f,o)\n",
    "    BS = bs.compute_BS()\n",
    "    BS_UNC = bs.compute_UNC()\n",
    "    BSS, BSS_all = bs.compute_BSS(return_all=True) \n",
    "    BSS2 = (BS-BS_UNC)/(0-BS_UNC)\n",
    "    assert BSS == BSS2\n",
    "    print(BSS)\n",
    "    print(BSS_all)\n",
    "\n",
    "    xes = CrossEntropy(f,o)\n",
    "    XES = xes.compute_XES()\n",
    "    XES_UNC = xes.compute_UNC()\n",
    "    XESS, XESS_all = xes.compute_XESS(return_all=True)\n",
    "    XESS2 = (XES-XES_UNC)/(0-XES_UNC)\n",
    "    assert XESS == XESS2\n",
    "    print(XESS)\n",
    "    print(XESS_all)\n",
    "    \n",
    "    plot_BSSvXESS(f,BSS,XESS)\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The lowest 5% forecast probabilities are highlighted in green. We see the slopes of the curves diverge markedly at the extreme probabilities. (i.e., the change in skill estimation with absolute probability error for a very rare or very common event). Also, BSS goes into negative (skill is lost) at different locations to XESS, meaning the two scores are not interchangeable at the extremes. __So how do we pick one?__ \n",
    "\n",
    "# Derivation of cross-entropy components\n",
    "\n",
    "So let's look more in depth at how the score behaves for predictions of rare, intermittent events forecast by an ensemble forecasting system (EFS) comprising $N_e$ members. Probabilities will be generated from the multiple integrations with small perturbations in the initial conditions (i.e., we assume a perfect model.) by performing the mean over the vector $\\mathbf{f}$ of Trues and Falses in the forecast at each $\\alpha$. Just like the Brier Score (ref), the XES can be decomposed (e.g., ref) into four different quantities as: \n",
    "\n",
    "$$\n",
    "\\textrm{XES} = \\underbrace{\\frac{1}{N_\\alpha} \\sum^{N_k}_{k=1} n_k D_\\textrm{KL} (\\mathbf{\\hat{o}_k} \\| \\mathbf{f_k}) }_\\textrm{Reliability (REL)} - \n",
    "\\underbrace{ \\frac{1}{N_\\alpha} \\sum^{N_k}_{k=1} n_k D_\\textrm{KL} (\\mathbf{\\hat{o}_k} \\| \\mathbf{\\hat{o}}) }_\\textrm{Discrimination (DSC)} + \n",
    "\\underbrace{ \\frac{1}{N_\\alpha} \\sum^{N_\\alpha}_{\\alpha=1} D_\\textrm{KL} (\\mathbf{o_\\alpha} \\| \\mathbf{\\hat{o}}) }_\\textrm{Uncertainty (UNC)} + \n",
    "\\underbrace{ \\frac{1}{N_\\alpha} \\sum^{N_\\alpha}_{\\alpha=1} H(\\mathbf{o_\\alpha}) }_\\textrm{Obs. Uncertainty} \\\\ \\hspace{1cm} \\\\\n",
    "\\forall k \\in \\{ \\underline{p},\\frac{1}{N_e},\\frac{2}{N_e},...,\\frac{N_e-1}{N_e},\\bar{p} \\}\n",
    "$$\n",
    "\n",
    "where $k$ is each probability level. The components DSC, REL, and UNC represent discrimination, reliability, and uncertainty:\n",
    "\n",
    "* Discrimination, DSC, is a measure of useful information gained, sometimes termed \"inherent goodness\". Also known as mutual information between forecasts and observations. Higher values are better.\n",
    "* Reliability, REL, is a measure of the \"probabilistic false alarm\", or the amount of information that proved detrimental to reducing surprise. In EFSs, this stems from, e.g., issuing a 30% chance of rain before 20% of the observed rain events. Lower values are better.\n",
    "* Uncertainty, UNC, is the portion of missing information that is inherent and inextricable to the system. This itself comprises both (1) observational error (2) and inherent uncertainty in the system regardless of observational error from chaotic error growth (exponential information loss). Lower values indicate higher predictability (lower inherent uncertainty).\n",
    "\n",
    "Both DSC and REL can be divided by UNC to create skill-score components - like the XESS, this is neither a proper score (ref) nor measured in bits, but it gives a normalised approximation. Unlike XESS, DSCSS and RELSS are positive values.\n",
    "\n",
    "The Brier Score can be written in components analogously as:\n",
    "\n",
    "$$\n",
    "\\textrm{BS} = \\underbrace{\\frac{1}{N_\\alpha}\\sum\\limits _{k=1}^{N_k}{n_{k}(\\mathbf{f_{k}}-\\mathbf{\\bar{o}}_{\\mathbf{k}})}^{2}}_\\textrm{Reliability (REL)} - \\underbrace{\\frac{1}{N_\\alpha}\\sum\\limits _{k=1}^{N_k}{n_{k}(\\mathbf{\\bar{o}_{k}}-\\bar{\\mathbf{o}})}^{2}}_\\textrm{Discrimination (DSC)} + \\underbrace{\\mathbf{\\bar{o}}\\left({1-\\mathbf{\\bar{o}}}\\right)}_\\textrm{Uncertainty (UNC)} + \\textrm{Obs Unc.}\n",
    "$$\n",
    "\n",
    "A rarer phenomenon has less average surprise (the assumption of \"no event\" is usually correct), but a much larger surprise when it occurs without being forecast. Hence, there is less uncertainty to resolve, averaged over the dataset, for a rarer event. We must therefore normalise the estimate of skill by the entropy of the observations, or the UNC component. We name this the Cross-Entropy Skill Score in analogy to (Weijs et al ref):\n",
    "\n",
    "$$\n",
    "\\textrm{XESS} = \\frac{\\textrm{DSC}}{\\textrm{UNC}} - \\frac{\\textrm{REL}}{\\textrm{UNC}} \\\\\n",
    "= \\frac{\\textrm{XES}-\\textrm{UNC}}{0-\\textrm{UNC}} \n",
    "$$\n",
    "\n",
    "To recapitulate: $H_\\times$ measures how many bits of missing information remain after a forecast is obtained, but XESS is a normalised measure whereby positive values indicate forecast skill that outperforms the use of the climatological mean. We might preserve the time-series of individual $D_\\textrm{KL}$ values in the above component equations. For rare and impactful events, the event itself may contribute very little to values averaged over a long time. This can be considered by attributing larger value to the most impactful events (i.e., where the event was observed; this mimicks societal preference for false alarms over missed hits). Here is an information-theoretical analogue of the 2x2 contingency table (confusion matrix) developed by (Green and Swets), where each box indicates the gain or loss of information in terms of surprise. We convert the probabilistic forecast into binary by setting $o$ to True if greater than 0.5, and False otherwise; a forecast is True if $f$ is greater than the end-user's critical probability $p_\\textrm{cric}$, above which the user takes action such as damage mitigation:\n",
    "\n",
    "| | $f > p_{\\textrm{crit}}$ | $f < p_{\\textrm{crit}}$ |\n",
    "| --- | --- | --- |\n",
    "| $o > 0.5$ | $+S$ | $-S(p)$ |\n",
    "| $o < 0.5$ | 0 | $+S(1-p)$ |\n",
    "\n",
    "Insufficiently accurate initial conditions, and erroneous models, mean that we will waste resources past a quantisation (resolution) level where observational error becomes substantial. This is the curse of dimensionality. Probability space is also insufficiently sampled: there will always be regions of phase space not explored by the EFS, plus an operational model is not perfect. (Herein, we use a perfect L63 model, but perturbations can be added during the integration to simulate model-error drift.) Further, smaller spatial and time scales require more information to determine its state. Hence, at finer spatial resolutions (which come with finer temporal resolutions due to the CFL criterion), more members are required to capture the larger information dimension of the pdf, or a kernel dressing should be applied to interpolate the pdf between members (refs). \n",
    "\n",
    "Next, we:\n",
    "* Generate example time-series observations and forecasts using the Lorenz-63 model (ref) to analyse information gain estimates with respect to:\n",
    "    - Intermittency (\"burstiness\" of the rare event)\n",
    "    - Ensemble membership\n",
    "    - Temporal quantisation\n",
    "    - Initial-condition uncertainty of the EFS\n",
    "\n",
    "We deploy Lorenz's 1963 three-variable model of convection (L63;ref) in an intermittent regime (i.e., extreme events come in infrequent bursts). We will create three ensembles, each varyign the $\\rho$ parameter from L63. We will also consider different levels of exceedence in the L63 z-variable time series, affecting the event frequency.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# double approx help from github user kendixon\n",
    "\n",
    "class L63:\n",
    "    def __init__(self,x0,y0,z0,rho,\n",
    "                 sigma=10.0,beta=2.667,dt=1e-3):\n",
    "        \"\"\"Chop off spin-up time?\n",
    "        \n",
    "        TODO:\n",
    "            * Remove spin-up time (transients)\n",
    "            * Add perturbations (with random seed allowed as arg)\n",
    "        the paper values use the intermittent regime\n",
    "        \n",
    "        rho\n",
    "        166.08 is intermittent\n",
    "        166.1 is the (more intermittent) variation\n",
    "        28.0 is the lorenz default\n",
    "        \"\"\"\n",
    "        self.sigma = sigma\n",
    "        self.rho = rho\n",
    "        self.beta = beta\n",
    "        self.dt =dt\n",
    "        self.t = 0.0\n",
    "\n",
    "        # The data array in x,y,z. Time appended in axis=3\n",
    "        self.output = np.zeros((3,1))\n",
    "        self.output[:,0] = [x0,y0,z0]\n",
    "        # self.\n",
    "\n",
    "    def dxdt(self,x,y):\n",
    "        return self.sigma * (y-x)\n",
    "\n",
    "    def dydt(self,x,y,z):\n",
    "        return x * (self.rho - z) - y\n",
    "\n",
    "    def dzdt(self,x,y,z):\n",
    "        return (x*y) - (self.beta * z)\n",
    "    \n",
    "    def double_approx(self,x,y,z):\n",
    "        dx1 = self.dxdt(x,y)*self.dt + x\n",
    "        dy1 = self.dydt(x,y,z)*self.dt + y\n",
    "        dz1 = self.dzdt(x,y,z)*self.dt + z\n",
    "        \n",
    "        dx2 = self.dxdt(dx1,dy1)\n",
    "        dy2 = self.dydt(dx1,dy1,dz1)\n",
    "        dz2 = self.dzdt(dx1,dy1,dz1)\n",
    "        \n",
    "        x2 = self.dt*(0.5*dx1 + 0.5*dx2) + x\n",
    "        y2 = self.dt*(0.5*dy1 + 0.5*dy2) + y\n",
    "        z2 = self.dt*(0.5*dz1 + 0.5*dz2) + z\n",
    "        \n",
    "        return x2,y2,z2\n",
    "        \n",
    "    def integrate_once(self,):\n",
    "        \"\"\"\n",
    "        Args:\n",
    "            t0    : Initial time\n",
    "            y0    : Intiial state\n",
    "        \"\"\"\n",
    "        x = self.output[0,-1]\n",
    "        y = self.output[1,-1]\n",
    "        z = self.output[2,-1]\n",
    "        # print(x,y,z)\n",
    "        ret = self.double_approx(x,y,z)\n",
    "        # print(ret)\n",
    "        new_time = np.expand_dims(np.array(ret),axis=1)\n",
    "        return new_time\n",
    "        \n",
    "    def integrate(self,n,clip=None):\n",
    "        # Time 0 is already done\n",
    "        for n in range(n-1):\n",
    "            next_data = self.integrate_once()\n",
    "            # print(self.output.shape,next_data.shape)\n",
    "            self.output = np.concatenate((self.output,next_data),axis=1)\n",
    "            # print(self.output)\n",
    "        #print(self.output)\n",
    "        if clip is not None:\n",
    "            self.output = self.output[:,clip:]\n",
    "            # cut_nt = nt - clip\n",
    "        return \n",
    "    \n",
    "    def get_z(self,):\n",
    "        return self.output[2,:]\n",
    "\n",
    "    def get_exceedence_ts(self,pc,vrbl='z'):\n",
    "        data = self.get_z()\n",
    "        exceed = self.get_pc_exceed(data,pc)\n",
    "        return exceed\n",
    "    \n",
    "    @staticmethod\n",
    "    def get_pc_exceed(ts,pc):\n",
    "        return ts > np.percentile(ts,pc)\n",
    "    \n",
    "    @staticmethod\n",
    "    def do_windowing(ts,wsize):\n",
    "        \"\"\"\n",
    "        Args:\n",
    "        ts    : time series of raw data.\n",
    "        wsize : epoch window size in time steps\n",
    "        \"\"\"\n",
    "        nt = ts.size\n",
    "        epochs = np.zeros(int(nt/wsize)+1).astype(bool)\n",
    "        widx = np.arange(0,nt,wsize)\n",
    "        windows = np.zeros_like(ts).astype(bool)\n",
    "        epoch_idx = np.arange(int(wsize/2),nt,wsize)\n",
    "        for nc,cidx in enumerate(widx):\n",
    "            idx0 = int(nc*wsize)\n",
    "            idx1 = int(idx0 + wsize)\n",
    "            this_window = ts[idx0:idx1]\n",
    "            epochs[nc] = True in this_window\n",
    "            windows[idx0:idx1] = bool(epochs[nc])    \n",
    "        return windows        \n",
    "    \n",
    "    @classmethod\n",
    "    def mask_and_window(cls,arr,pc,wsize):\n",
    "        if arr.ndim == 1:\n",
    "            nt = arr.size\n",
    "            arr = arr.reshape(1,nt)\n",
    "        pc_arr = cls.get_pc_exceed(arr,pc)\n",
    "        window_ts = np.zeros_like(pc_arr)\n",
    "        for ne in range(pc_arr.shape[0]):\n",
    "            window_ts[ne,:] = cls.do_windowing(pc_arr[ne,:],wsize)\n",
    "        # window_ts = cls.do_windowing(pc_ts,wsize)\n",
    "        return window_ts\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the z-component (in L63, this represents rising motion in convective cells) as our time series. We convert to a binary forecast by identifying where the time series exceeds the given percentile. We then quantise in time to both introduce tolerance and mimic EFS real-world output.\n",
    "\n",
    "* The parameter $\\rho$ controls the intermittency (burstiness), which through chaotic growth (see Lyapunov exponents) destroys information about the state of the system until no further information can be obtained. This asymptote is an analogue of Lorenz's predictability limit after error saturation.\n",
    "* Window size $w$ mimicks the time- and spatial-scale-aware quantisation methods used by verification scores such as Fractions Skill Score (refs). Within each overlapping window segment, moved to centre on each data point from the raw time series, we assign True for that $\\alpha$ if the percentile is exceeded for that member or observation; False otherwise. This avoids verifying the noisy sample-by-sample evaluation, which is analgous to evaluating a weather model at the truncated scale (i.e., gridpoint-to-gridpoint). It also reflects the tolerance for false alarms over missed hits. The larger the value of $w$, the more likely a forecast captures the event, but the less information is transferred about the precise time of occurrence when interpreted in real-time $\\tau$. Hence we may expect an optimised sweet-spot.\n",
    "* The larger the percentile, the rarer the event. Coupled with the variation in regime via $\\rho$, this generally represents the forecast challenge of rare, intermittent (bursty) events such as supercell outbreaks.\n",
    "* We run 200 ensemble members ($N_e = 200$), but subsample randomly to create smaller ensembles of various sizes. This will assess the benefit of increasing the required computer resources linearly despite diminishing returns: each additional member contains more mutual information with the rest of the EFS members, meaning increasingly less novel (correct) information is acquired), and increasingly more incorrect information. Hence we may expect a sweet-spot.\n",
    "* We perturb initial conditions with a range of small perturbations from a uniform distribution to reflect initial condition uncertainty. As these perturbvations increase in magnitude, predictability is lost more quickly (information gain cumulates over time logarithmically). However, no perturbations at all would defeat the purpose of an EFS. Hence, we expect another sweet-spot.\n",
    "\n",
    "The high number of permutations of the above parameters mean we aggregate over all other parameters for each given parameters: this is not rigorous but will show the variability of information gain (and BSS versus XESS) as a function of $\\rho$, $w$, percetile, and $N_e$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_z_ts(x0,y0,z0,rho,nt_all,cutoff):\n",
    "    lorenz = L63(x0,y0,z0,rho=rho)\n",
    "    lorenz.integrate(nt_all,clip=cutoff)\n",
    "    return lorenz.get_z()\n",
    "\n",
    "def make_itr(tweak_max,Ne):\n",
    "    tweaks = np.random.uniform(-tweak_max,tweak_max,Ne)\n",
    "    for ne, tweak in enumerate(tweaks):\n",
    "        yield x0+tweak, y0+tweak, z0+tweak, rho\n",
    "\n",
    "def do_run(i):\n",
    "    x0,y0,z0,rho = i\n",
    "    lorenz = L63(x0,y0,z0,rho)\n",
    "    lorenz.integrate(nt_all,clip=cutoff)\n",
    "    return lorenz.get_z()\n",
    "\n",
    "def generate_ensemble(x0,y0,z0,Ne,rho,tweak_max,nt_all,cutoff):\n",
    "    itr = make_itr(tweak_max,Ne)\n",
    "    if ncpus == 1:\n",
    "        ensemble = np.zeros([Ne,nt_all-cutoff])\n",
    "        for n,i in enumerate(itr):\n",
    "            print(f\"Computing member {n} of {Ne}\")\n",
    "            ensemble[n,:] = do_run(i)\n",
    "    else:\n",
    "        with multiprocessing.Pool(ncpus) as pool:\n",
    "            returns = pool.map(do_run,itr)\n",
    "        ensemble = np.array(returns)\n",
    "    return ensemble\n",
    "                      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Want truth to be window size 100?\n",
    "# Then test with 200,500 window size\n",
    "\n",
    "# Catalogue of runs\n",
    "# first idx is obs, rest is fcst \n",
    "\n",
    "\n",
    "# Do numerous \"days\" for a longer experiment?\n",
    "\n",
    "# Parallelisation?\n",
    "# Initial conditions\n",
    "x0 = 1.0\n",
    "y0 = 1.0\n",
    "z0 = 1.0\n",
    "\n",
    "# Total integration steps\n",
    "nt_all = 25000\n",
    "\n",
    "# Spin-up period to be removed\n",
    "cutoff = 15000\n",
    "\n",
    "# Window size for scale-aware\n",
    "# Can do later when computing scores\n",
    "# wsize = 500 \n",
    "\n",
    "# Number of ensemble members\n",
    "# Can subset later\n",
    "Ne = 200\n",
    "\n",
    "# Number of time steps in raw data after spin-up is cut off\n",
    "nt = nt_all - cutoff\n",
    "\n",
    "# Measure of intermittency\n",
    "rho_list = [166.08,166.09,166.1]\n",
    "# percentiles for thresholding\n",
    "pc_list = [93.5,95.5,97.5,98.5,99.5,99.9]\n",
    "# IC perturbations for ensembles\n",
    "tweak_max_list = [1e-8,1e-7,5e-7,1e-6,1e-5] \n",
    "\n",
    "# wsize = 100 will be minimum (i.e. truth)\n",
    "wsize_list = [25,50,100,250,500,1000]\n",
    "\n",
    "# Raw\n",
    "OBS = {r:{} for r in rho_list}\n",
    "FCST = {r:{t:{} for t in tweak_max_list} for r in rho_list}\n",
    "\n",
    "# Windowed and percentile-masked\n",
    "OBS_PW = {r:{pc:{w:0 for w in wsize_list} for pc in pc_list} \n",
    "                for r in rho_list}\n",
    "FCST_PW = {r:{t:{pc:{w:0 for w in wsize_list} for pc in pc_list}\n",
    "                for t in tweak_max_list} for r in rho_list}\n",
    "    \n",
    "for rho, tweak_max in itertools.product(rho_list,tweak_max_list):\n",
    "    obs_file = f'data/obs_{int(rho*1000)}.npy'\n",
    "    fcst_file = f'data/fcst_{int(rho*1000)}_{int(tweak_max*1e8)}.npy'\n",
    "    print(fcst_file)\n",
    "\n",
    "    # Create observation data (only once per IC perturbation value (tweak_max))\n",
    "    if tweak_max == tweak_max_list[0]:\n",
    "        if os.path.exists(obs_file):\n",
    "            OBS[rho] = np.load(obs_file)\n",
    "        else:\n",
    "            obs = generate_z_ts(x0,y0,z0,rho,nt_all,cutoff)\n",
    "            OBS[rho] = obs.reshape(1,obs.size)\n",
    "            np.save(arr=obs,file=obs_file)\n",
    "\n",
    "        for pc, wsize in itertools.product(pc_list,wsize_list):\n",
    "            ts_pw = L63.mask_and_window(OBS[rho],pc,wsize)\n",
    "            OBS_PW[rho][pc][wsize] = ts_pw\n",
    "    \n",
    "    # Create forecast data\n",
    "    if os.path.exists(fcst_file):\n",
    "        FCST[rho][tweak_max] = np.load(fcst_file)\n",
    "    else:\n",
    "        if __name__ == \"__main__\":\n",
    "            fcst = generate_ensemble(x0,y0,z0,Ne,rho,tweak_max,nt_all,cutoff)\n",
    "        FCST[rho][tweak_max] = fcst\n",
    "        np.save(arr=fcst,file=fcst_file)\n",
    "            \n",
    "    for pc, wsize in itertools.product(pc_list,wsize_list):\n",
    "        ts_pw = L63.mask_and_window(FCST[rho][tweak_max],pc,wsize)\n",
    "        FCST_PW[rho][tweak_max][pc][wsize] = ts_pw\n",
    "        \n",
    "# [Ne x Nt]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below is an example of the first two EFS members for $\\rho = 0$, $\\textrm{ws} = 0$, exceeding a XXth percentile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_example(raw_obs,obs_pw,raw_fcst,fcst_pw,pc,wsize,title,nnn):\n",
    "    # Plotting\n",
    "    fig,axes = plt.subplots(figsize=(8,8),ncols=1,nrows=4)\n",
    "\n",
    "    nt = raw_obs.size\n",
    "    raw_xlocs = np.arange(nt)\n",
    "    window_locs = np.arange(wsize,nt+wsize,wsize)\n",
    "    x_lines = np.arange(0,nt+wsize,wsize)\n",
    "\n",
    "    axes.flat[0].plot(raw_xlocs,raw_obs,color='b')\n",
    "    axes.flat[0].set_title(\"Observations\",**ff)\n",
    "    axes.flat[0].hlines(y=np.percentile(raw_obs,pc),\n",
    "                           xmin=0,xmax=nt,color='k')\n",
    "\n",
    "    # axes.flat[1].fill(raw_xlocs,obs_pw,color='m')\n",
    "    # axes.flat[1].plot(raw_xlocs,obs_pw,color='m')\n",
    "    axes.flat[1].fill_between(x=window_locs,y1=obs_pw[::wsize],\n",
    "                                y2=0,step='pre',color='b')\n",
    "    axes.flat[1].vlines(x=x_lines,ymin=0,ymax=1,color='k')\n",
    "    \n",
    "    axes.flat[2].set_title(f\"Forecast member {nnn}\",**ff)\n",
    "    axes.flat[2].plot(raw_xlocs,raw_fcst,color='r')\n",
    "    axes.flat[2].hlines(y=np.percentile(raw_fcst,pc),\n",
    "                           xmin=0,xmax=nt,color='k')\n",
    "        \n",
    "    axes.flat[3].fill_between(x=window_locs,y1=fcst_pw[::wsize],\n",
    "                                y2=0,step='pre',color='r')\n",
    "\n",
    "    axes.flat[3].vlines(x=x_lines,ymin=0,ymax=1,color='k')\n",
    "    # axes.flat[3].fill(raw_xlocs,fcst_pw,color='k')   \n",
    "    \n",
    "    # make prettier\n",
    "    for n in (1,3):\n",
    "        axes.flat[n].hlines(y=(0,1),xmin=0,xmax=nt,color='k')\n",
    "    for n in range(4):\n",
    "        axes.flat[n].tick_params(axis='both',which='both',bottom=False,left=False,\n",
    "                                    top=False,labelbottom=False,labelleft=False)  \n",
    "        axes.flat[n].axis('off')\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rho = 166.08\n",
    "pc = 98.5\n",
    "wsize = 500\n",
    "tweak_max = 5e-7\n",
    "raw_fcst = FCST[rho][tweak_max][0,:]\n",
    "fcst_pw = FCST_PW[rho][tweak_max][pc][wsize][0,:]\n",
    "\n",
    "# Weirdness with dimensions\n",
    "raw_obs = OBS[rho]\n",
    "nnn = 0\n",
    "if raw_obs.ndim == 2:\n",
    "    raw_obs = raw_obs[nnn,:]\n",
    "obs_pw = OBS_PW[rho][pc][wsize][nnn,:]\n",
    "\n",
    "plot_example(raw_obs,obs_pw,raw_fcst,fcst_pw,\n",
    "             pc,wsize,f\"{rho}\",nnn=nnn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nnn = 1\n",
    "raw_fcst = FCST[rho][tweak_max][nnn,:]\n",
    "fcst_pw = FCST_PW[rho][tweak_max][pc][wsize][nnn,:]\n",
    "\n",
    "plot_example(raw_obs,obs_pw,raw_fcst,fcst_pw,\n",
    "             pc,wsize,f\"{rho}\",nnn=nnn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We note the increasing divergence over time between the three time series. Let's compare the verification of BS and XESS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def subsample(ensemble,Ne):\n",
    "    mems = random.sample(range(200),Ne)\n",
    "    return ensemble[mems,:]\n",
    "\n",
    "score_components = [\"DSC\",\"DSCSS\",\"REL\",\"RELSS\",\"H\",\"UNC\",\"XES\",\"COMPSUM_XES_noH\",\"COMPSUM_XES\",\n",
    "                        \"BS\",\"XESS\",\"COMPSUM_XESS_noH\",\"COMPSUM_XESS\",\"BSS\",]\n",
    "feature_list = ['rho','tweak_max','pc','wsize','Ne',]\n",
    "col_names = feature_list + score_components\n",
    "\n",
    "# scores = pd.DataFrame()#columns=col_names)\n",
    "scores = None\n",
    "\n",
    "Ne_list = [5,10,20,30,40,50,100,150,200]\n",
    "# Verification\n",
    "for rho, tweak_max, pc, wsize, Ne in itertools.product(\n",
    "        rho_list, tweak_max_list, pc_list, wsize_list, Ne_list):\n",
    "    \n",
    "    # Randomly select members\n",
    "    raw_f = FCST_PW[rho][tweak_max][pc][wsize][:,::wsize]\n",
    "    raw_p = subsample(raw_f,Ne)\n",
    "    probs = InfoGain.bound(np.mean(raw_p,axis=0),Ne=Ne)\n",
    "    # obs = OBS_PW[rho][pc][wsize][0,::wsize].astype(int)\n",
    "    # obs = InfoGain.bound(OBS_PW[rho][pc][wsize][0,::wsize].astype(int),thresh=0.0005)\n",
    "    obs = InfoGain.bound(OBS_PW[rho][pc][wsize][0,::wsize].astype(float),thresh=0.001)\n",
    "    # pdb.set_trace()\n",
    "    IG = InfoGain(f=probs,o=obs)\n",
    "    BrS = BrierScore(f=probs,o=obs)\n",
    "    # XES, XES_all = IG.compute_DS()\n",
    "    # XES, XES_all = IG.compute_XES()\n",
    "    XES = IG.compute_XES()\n",
    "    # XESS = IG.compute_DSS()\n",
    "    XESS = IG.compute_XESS()\n",
    "    \n",
    "    DSC = IG.compute_DSC()\n",
    "    REL = IG.compute_REL()\n",
    "    UNC = IG.compute_UNC()\n",
    "    DSCSS = DSC/UNC\n",
    "    RELSS = REL/UNC\n",
    "    H = IG.compute_H()\n",
    "\n",
    "    # Compare with adding components?    \n",
    "    COMPSUM_XES = IG.compute_XES()\n",
    "    COMPSUM_XES_noH = IG.compute_XES(with_H=False,from_components=True)\n",
    "\n",
    "    COMPSUM_XESS = IG.compute_XESS()\n",
    "    COMPSUM_XESS_noH = IG.compute_XESS(with_H=False,from_components=True)\n",
    "    \n",
    "    # BS = compute_BS(f=probs,o=obs)\n",
    "    BS = BrS.compute_BS()\n",
    "    # BSS = compute_BSS(f=probs,o=obs)\n",
    "    BSS = BrS.compute_BSS()\n",
    "\n",
    "    # pd.DataFrame({'f':probs,'o':obs,'DS_t':DS_all}).head(10)\n",
    "    # new_df = pd.DataFrame.from_records(list(\n",
    "    new_df = pd.DataFrame.from_dict(\n",
    "                    {'rho':rho,'tweak_max':tweak_max*1e8,\n",
    "                          'pc':pc,'wsize':wsize,'Ne':Ne,'COMPSUM_XES':COMPSUM_XES,\n",
    "                            'XES':XES,'XESS':XESS,\"H\":H,'UNC':UNC,'REL':REL,\n",
    "                           'DSC':DSC,'DSCSS':DSCSS,'RELSS':RELSS,\n",
    "                            'BS':BS,'BSS':BSS,\"COMPSUM_XESS\":COMPSUM_XESS,\n",
    "                           \"COMPSUM_XES_noH\":COMPSUM_XES_noH,\n",
    "                           \"COMPSUM_XESS_noH\":COMPSUM_XESS_noH},orient=\"index\").T\n",
    "    if scores is None:\n",
    "        scores = new_df\n",
    "    else:\n",
    "        scores = pd.concat([scores, new_df],axis=0,ignore_index=True)\n",
    "    \n",
    "scores['wsize'] = scores['wsize'].astype('int16')\n",
    "scores['tweak_max'] = scores['tweak_max'].astype('int16')\n",
    "scores['Ne'] = scores['Ne'].astype('int16')\n",
    "# display(scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def visualise_scores(ax,scores,diag,score,title_cap):\n",
    "    cl = copy.copy(col_names)\n",
    "    cl.remove(diag)\n",
    "    diag_vals = scores[diag].unique()\n",
    "    data = np.zeros([len(diag_vals),len(rho_list)])\n",
    "    for ndv,dv in enumerate(diag_vals):\n",
    "        for nrho, rho in enumerate(rho_list):\n",
    "            sub_df = scores[(scores[diag]==dv) & (scores['rho']==rho)]\n",
    "            data[ndv,nrho] = sub_df[score].mean() \n",
    "    \n",
    "    # fig,ax = plt.subplots(1,figsize=(7,7))\n",
    "    x = np.arange(len(diag_vals))\n",
    "    ax.bar(x-0.3,data[:,0],width=0.1,color='r',align='center',label=rho_list[0])\n",
    "    ax.bar(x-0.1,data[:,1],width=0.1,color='g',align='center',label=rho_list[1])\n",
    "    ax.bar(x+0.1,data[:,2],width=0.1,color='b',align='center',label=rho_list[2])\n",
    "    #ax.bar(x+0.3,data[:,3],width=0.1,color='r',align='center',label=rho_list[2])\n",
    "    \n",
    "    ax.plot(x-0.3,data[:,0],color='r')\n",
    "    ax.plot(x-0.1,data[:,1],color='g')\n",
    "    ax.plot(x+0.1,data[:,2],color='b')\n",
    "    #ax.plot(x+0.3,data[:,3],color='r')\n",
    "    \n",
    "    ax.hlines(y=0,xmin=0,xmax=x[-1],color='k')    \n",
    "    ax.set_xticks(x)\n",
    "    if diag.endswith(\"SS\"):\n",
    "        ylab = \"Skill Score\"\n",
    "    else:\n",
    "        ylab = \"Entropy (bits)\"\n",
    "    ax.set_ylabel(ylab)\n",
    "    ax.set_xticklabels(diag_vals)\n",
    "    ax.legend(facecolor='lightblue')\n",
    "    \n",
    "    ax.set_title(f\"({title_cap}) {score} with respect to {diag}\",**ff)\n",
    "    # display(fig)\n",
    "    # pdb.set_trace()\n",
    "    return ax\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for bar_score in score_components:#bar_score = \"DSC\"\n",
    "def plot_component(bar_score):\n",
    "    fig,axes = plt.subplots(ncols=2,nrows=2,figsize=(10,8))\n",
    "    _ = visualise_scores(axes.flat[0],scores,\"Ne\",bar_score,\"a\")\n",
    "    _ = visualise_scores(axes.flat[1],scores,\"tweak_max\",bar_score,\"b\")\n",
    "    _ = visualise_scores(axes.flat[2],scores,\"wsize\",bar_score,\"c\")\n",
    "    _ = visualise_scores(axes.flat[3],scores,\"pc\",bar_score,\"d\")\n",
    "    return\n",
    "\n",
    "score_itr = iter(score_components)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Explain DSCSS and RELSS, Toth paper?)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_component(next(score_itr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_component(next(score_itr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As ensemble membership increases, increasing useful information (DSC, DSCSS) is acquired by the EFS. We note the diminishing returns as membership increases. As IC error increases, DSCSS understandably decreases. As larger windowing is used DSC is larger. As the rarity of the event increases, the useful information acquired markedly drops. However, the rarer the event, the more sensitive the score to an event, and the larger the range of information loss/gain (bimodal)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_component(next(score_itr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_component(next(score_itr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reliability error (REL) increases with ensemble membership, just as DSC, as more members increase the chance the individual ensemble member has delivered erroneous information. If DSC > REL (or DSCSS > RELSS), erroneous information does not completely erode the useful information acquired, and hence information is gained.\n",
    "\n",
    "Reliability error increases with IC error, though it varies highly with regime ($\\rho$). This is seen in RELSS, and hence is independent of inherent uncertainty (UNC), and suggests the calibration error is a function of the flow. The lowest REL error occurs at smaller window sizes, though the $\\rho = 166.08$ and $= 166.09$ regimes show a minimum near 50 and 100 unit times, respectively. This is analgous to identifying the smallest skillful scale: evaluating data without time--space tolerance will appear poor due to a multitude of phase errors. The variation of optimal window size with regime reiterates the need to evaluate on numerous time--space scales.\n",
    "\n",
    "Finally, REL error also has a minimum around the 99th percentile. For events rarer than this, error increases rapidly, again potentially due to the low sample size for this percentile. A rarer event requires more data points to evaluate average surprise appropriately. It may therefore be useful to preserve surprise or XES/XESS per epoch $\\alpha$ in a time series to identify events that were best or worst for removing surprise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_component(next(score_itr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_component(next(score_itr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_component(next(score_itr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_component(next(score_itr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_component(next(score_itr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_component(next(score_itr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_component(next(score_itr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_component(next(score_itr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_component(next(score_itr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_component(next(score_itr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparing BSS and XESS, we see that XESS is more sensitive to poor probabilities than BSS for a small number of ensemble members. The two scores are similar for evaluating skill with respect to IC error and window size. But at extreme percentiles (d), the XESS punishes poor forecasts more severely in the top percentile. It is this non-negligible discrepancy, coupled with the mathematical support of information theory, that supports use of XES and XESS especially in applications where high value is placed on useful information about rare events (i.e., usually high-impact hazards). For instance, forecasting extreme events such as extreme percentile values for rotating supercells may be more accurately evaluated with information theory.\n",
    "\n",
    "Hence, information gain is the difference between prior (naïve) and posterior uncertainty. The prior uncertainty may simply be UNC; alternatively, we discuss how a prior estimate of uncertainty can be computed from a time series to form a rate of information loss (Kolmogorov-Sinai Entropy or $H_\\mu$; ref). Hence, we can form an information \"stream\" from naivity to the end-user's decision, whereby the EFS, human forecasters, emergency managers, etc represent \"information dams\". At each dam, the goal is to reduce uncertainty and maximise information gain. The components can be conserved throughout the flow of information."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summary of repercussions\n",
    "Here:\n",
    "\n",
    "- Scores based on BS, like FSS/eFSS, has the same problem as BS. Can just switch the scoring rule from MSE to XES\n",
    "- Next\n",
    "    \n",
    "Blah."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Some definitions:\n",
    "\n",
    "| Symbol| Meaning |\n",
    "| ------ | ------- |\n",
    "| $N_e$    | Ensemble membership |\n",
    "| $N_\\alpha$ | Number of \"epochs\" or \"unit time steps\" |\n",
    "| $\\tau$ | Total time (sec) |\n",
    "| $H$ | Average entropy or average surprise |\n",
    "| $S$ | Self-information or event surprise (bits) |\n",
    "| $H_\\times$ | Cross-entropy (bits) |\n",
    "| $\\textrm{IG}$ | Information gain (bits) |\n",
    "| $I(x{;}y)$ | Mutual information between $x$ and $y$ |\n",
    "| $$~D_{KL}(x\\|y)$$ | Kullback-Liebler Divergence between $x$ and $y$ |\n",
    "| $f$ | Forecasted probability of event |\n",
    "| $o$ | Probability that event occurred (i.e., observation) |\n",
    "| $\\mathbf{f}$ | 2-D array of forecast probabilities: $\\alpha$ by $i$. |\n",
    "| $\\mathbf{o}$ | 2-D array of probability of occurrence: $\\alpha$ by $i$. |\n",
    "| $\\bar{o}$ | Frequency of observed event |\n",
    "| REL | Reliability of the forecast set|\n",
    "| DSC | Discrimination (or Brier-type resolution) |\n",
    "| UNC | Uncertainty of the observed frequency |\n",
    "\n",
    "Selected mathematical symbols:\n",
    "\n",
    "| Symbol| Meaning |\n",
    "| ------ | ------- |\n",
    "| $\\in$   | Belonging to the set of |\n",
    "| $\\mathbb{Z}^+$ | The set of non-zero positive integers "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Bibliography (also see paper)\n",
    "\n",
    "Here."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}