04a-data-verifying_punctual_metric(slow).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Importing modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\SERRA_R\\AppData\\Local\\Continum\\anaconda3\\lib\\site-packages\\statsmodels\\compat\\pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.\n",
      "  from pandas.core import datetools\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a49f64fef0>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a49f64fef0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import time\n",
    "import warnings\n",
    "import itertools\n",
    "from scipy import stats as st\n",
    "import statsmodels as sm\n",
    "import statsmodels.api as sma\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('fivethirtyeight')\n",
    "import h5py\n",
    "from statsmodels.distributions.empirical_distribution import ECDF\n",
    "from scipy.stats.kde import gaussian_kde\n",
    "from scipy.stats import norm\n",
    "import statisticalTools as sT\n",
    "from sklearn.neighbors.kde import KernelDensity\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "%matplotlib inline \n",
    "#%matplotlib \n",
    "plt.figure(figsize=(15, 6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Defining functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a4a37bdac8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline \n",
    "#%matplotlib \n",
    "\n",
    "plt.figure(figsize=(15, 6))\n",
    "\n",
    "\n",
    "def find_nearest(array, value):\n",
    "    '''Finds the closest value to a given value in an array. Returns closest value'''\n",
    "    array = np.asarray(array)\n",
    "    idx = (np.abs(array - value)).argmin()\n",
    "    return array[idx]\n",
    "\n",
    "def mid(x):\n",
    "    '''Returns a vector with the mean of two consecutive values.'''\n",
    "    middle = np.zeros(len(x)-1)\n",
    "    for i in range(len(x)-1):\n",
    "        middle[i] = (x[i+1]+x[i])/2\n",
    "    return middle\n",
    "\n",
    "def give_me(x0,x1,prob, tol = 0.03,est = 0):\n",
    "    '''Gives the probability of the values to me inside a zone [va,vb], the most probable [vstar], \n",
    "    and the interval length [abs(va-vb)]. Requires the data of a histogram.'''\n",
    "    \n",
    "    #The sum of a distribution is <= 1\n",
    "    if prob > 1:\n",
    "        return print('Put a value <= than 1')\n",
    "    if est == 0:\n",
    "        est = np.diff(x1)*x0 #finds the probabilities of each interval\n",
    "    else:\n",
    "        est = x1*x0\n",
    "    x_max = np.argmax(est) #finds the biggest probability in the interval\n",
    "    probb = est[x_max]\n",
    "    inds = [x_max]\n",
    "    x_act_min = x_max\n",
    "    x_act_max = x_max\n",
    "    count = 0\n",
    "    vstar = (x1[x_max]+x1[x_max+1])/2\n",
    "    \n",
    "    while probb + tol < prob:\n",
    "        #print(x_act_min)\n",
    "        count = count + 1\n",
    "        if x_act_min  > 0:\n",
    "            neig0 = est[x_act_min -1]\n",
    "        else:\n",
    "            neig0 = -1\n",
    "            \n",
    "        if x_act_max + 1 < len(x0):\n",
    "            neig1 = est[x_act_max +1]\n",
    "        else:\n",
    "            neig1 = -1\n",
    "        if (neig0 > neig1):\n",
    "            #print(neig0,neig1)\n",
    "            probb = probb + neig0\n",
    "            x_act_min = x_act_min - 1\n",
    "        elif (neig0 < neig1):\n",
    "            #print(neig0,neig1)\n",
    "            probb = probb + neig1\n",
    "            x_act_max = x_act_max + 1\n",
    "        else:\n",
    "            #print(neig0,neig1)\n",
    "            est[x_max] = 0\n",
    "            if x_act_min  > 0 and x_act_max + 1 < len(x0):\n",
    "                x_mmin = np.argmax(est[0:x_act_min])\n",
    "                x_mmax = np.argmax(est[x_act_max:])\n",
    "                if est[x_mmin] > est[x_mmax+x_act_max]:\n",
    "                    probb = probb + neig0\n",
    "                    x_act_min = x_act_min - 1 \n",
    "                else:\n",
    "                    probb = probb + neig1\n",
    "                    x_act_max = x_act_max + 1\n",
    "            elif x_act_min > 0:\n",
    "                probb = probb + neig0\n",
    "                x_act_min = x_act_min - 1 \n",
    "            elif x_act_max + 1 < len(x0):\n",
    "                probb = probb + neig1\n",
    "                x_act_max = x_act_max + 1\n",
    "            else:\n",
    "                #print(x_act_min,x_act_max)\n",
    "                warnings.warn(\"Warning: Verify the function give_me or try reducing the probability to 0.99\")\n",
    "                va = x1[x_act_min]\n",
    "                vb = x1[1 +x_act_max]\n",
    "                return probb, va, vb, abs(va-vb), vstar\n",
    "                #continue\n",
    "                \n",
    "        if count > len(x0):\n",
    "            return print('Impossible to finish')\n",
    "        \n",
    "    #print(x_act_max)\n",
    "    va = x1[x_act_min]\n",
    "    vb = x1[1 +x_act_max]\n",
    "    return probb, va, vb, abs(va-vb), vstar\n",
    "\n",
    "def proba(x0,x1):\n",
    "    '''Finds the probabilities of each interval'''\n",
    "    return np.diff(x1)*x0\n",
    "\n",
    "def verify_res(nforc,method,probability = 1, tole = 'None'):\n",
    "    '''Verify the interval of a distribution given the ensemble forecast data and a proability.'''\n",
    "    if method == 'hist10':\n",
    "        ai0, ai1 = np.histogram(nforc, bins = 10, density = 'True')\n",
    "        if tole == 'None':\n",
    "            tole = 0.05\n",
    "        prob, va, vb, dist, vstar = give_me(ai0,ai1,probability,tol = tole)\n",
    "    elif method == 'kde':\n",
    "        KDEpdf = gaussian_kde(nforc)\n",
    "        mean = np.abs(np.mean(nforc))\n",
    "        xmin = np.min(nforc)- .25*mean\n",
    "        xmax = np.max(nforc)+ .25*mean\n",
    "        x = np.linspace(xmin,xmax,100)\n",
    "        xprob = KDEpdf(x)\n",
    "        xprob = (xprob[:-1]+xprob[1:])/2\n",
    "        if tole == 'None':\n",
    "            tole = 0.0001\n",
    "        prob, va, vb, dist, vstar = give_me(xprob,x,probability,tol = tole)\n",
    "    elif method == 'kde_modif':\n",
    "        KDEpdf = gaussian_kde(nforc)\n",
    "        samp = KDEpdf.resample([100])\n",
    "        ai0, ai1 = np.histogram(samp, bins = 10, density = 'True')\n",
    "        if tole == 'None':\n",
    "            tole = 0.05\n",
    "        prob, va, vb, dist, vstar = give_me(ai0,ai1,probability,tol = tole)\n",
    "    elif method  == 'best_dist':\n",
    "        name, param = sT.best_fit_distribution(nforc)\n",
    "        result = getattr(st, name)\n",
    "        prob = probability\n",
    "        va,vb = result.interval(prob,*param)\n",
    "    elif method == 'norm':\n",
    "        param = norm.fit(nforc)\n",
    "        prob = probability\n",
    "        va,vb = norm.interval(prob,*param)\n",
    "    else:\n",
    "        return print('Method not known')\n",
    "    return prob, va, vb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reading data and extracting information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Read data\n",
    "file_name = 'data/GENS_test_AllLevels3'+'.hdf'\n",
    "dt = pd.read_hdf(file_name, 'Wind_vector')\n",
    "#Sort the dataset in function of date\n",
    "dt = dt.set_index('Timestamp')\n",
    "dt = dt.sort_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Read data\n",
    "file_name = 'data/gens-a_3'+'.hdf'\n",
    "DT = pd.read_hdf(file_name, 'Wind_vector')\n",
    "DT = DT.set_index('Timestamp')\n",
    "DT = DT.sort_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Selects the time steps presents on both datasets.\n",
    "ind0 = dt.index.unique()\n",
    "ind1 = DT.index.unique()\n",
    "ind = set(ind0) & set(ind1)\n",
    "\n",
    "dtt = dt.loc[ind]\n",
    "DTT = DT.loc[ind]\n",
    "\n",
    "dtc = dtt.where(dtt.Forecast == '000').where(dtt.NForecast == 999).dropna(how='all') #Observation - Most probable forecast\n",
    "dtd = dtt.where(dtt.Forecast == '006').where(dtt.NForecast == 999).dropna(how='all') #Forecast - Most probable forecast"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Selects the latitude, longitudes and levels steps presents on both datasets.\n",
    "lats0 = dtt.where(dtt.NForecast == 999).dropna(how = 'all').lat.unique()\n",
    "longs0 = dtt.where(dtt.NForecast == 999).dropna(how = 'all').long.unique()\n",
    "levels0 = dtt.where(dtt.NForecast == 999).dropna(how = 'all').level.unique()\n",
    "\n",
    "lats1 = DTT.lat.unique()\n",
    "longs1 = DTT.long.unique()\n",
    "levels1 = DTT.level.unique()\n",
    "\n",
    "lats = list(set(lats0) & set(lats1))\n",
    "longs = list(set(longs0) & set(longs1))\n",
    "levels = list(set(levels0) & set(levels1))\n",
    "\n",
    "#shows the latitudes, longitudes and levels presents in both datasets\n",
    "print(lats)\n",
    "print(longs)\n",
    "print(levels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Selecting data to compare\n",
    "\n",
    "dx = dt\n",
    "dt = dtt\n",
    "dt = dt.sort_index()\n",
    "dt = dt.where(dt.Forecast == '000')\n",
    "dt = dt.where(dt.NForecast == 999)\n",
    "dy = dt.dropna()\n",
    "#dt = dt.ix[dt.index[1]]\n",
    "\n",
    "dt = dtt\n",
    "dt = dt.sort_index()\n",
    "dt = dt.where(dt.Forecast == '006')\n",
    "dt = dt.where(dt.NForecast == 999)\n",
    "dyf = dt.dropna()\n",
    "\n",
    "dt = dx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "oi = 0\n",
    "\n",
    "dfp = DTT\n",
    "dyp = dy\n",
    "dyfp = dyf\n",
    "\n",
    "dfp = dfp.where(dfp.lat == lats[oi] )\n",
    "dfp = dfp.where(dfp.long == longs[oi] )\n",
    "dfp = dfp.where(dfp.level == levels[oi] )\n",
    "dfp = dfp.dropna(how = 'all')\n",
    "\n",
    "dyp = dyp.where(dyp.lat == lats[oi] )\n",
    "dyp = dyp.where(dyp.long == longs[oi]) \n",
    "dyp = dyp.where(dyp.level == levels[oi]) \n",
    "dyp = dyp.dropna(how='all')\n",
    "\n",
    "dyfp = dyfp.where(dyfp.lat == lats[oi]) \n",
    "dyfp = dyfp.where(dyfp.long == longs[oi]) \n",
    "dyfp = dyfp.where(dyfp.level == levels[oi]) \n",
    "dyfp = dyfp.dropna(how='all')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Verify results using punctual metric"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#defining variables and parameters\n",
    "count_loop = 0\n",
    "\n",
    "param = 'Wy'\n",
    "\n",
    "count = 0\n",
    "count0 = 0\n",
    "sizet = 0\n",
    "evalu = 0\n",
    "\n",
    "ind0 = dtc.index.unique()\n",
    "ind1 = DTT.index.unique()\n",
    "ind = set(ind0) & set(ind1)\n",
    "\n",
    "dtc = dtc.loc[ind]\n",
    "DTT = DTT.loc[ind]\n",
    "dtd = dtd.loc[ind]\n",
    "\n",
    "lats = DTT.lat.unique()\n",
    "longs = DTT.long.unique()\n",
    "levels = DTT.level.unique()\n",
    "\n",
    "#if one wants to test fast if it works\n",
    "#levels = levels[6]\n",
    "#lats = lats[0:1]\n",
    "#longs = longs[0:1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_t = time.time()\n",
    "for i in levels:\n",
    "    count_loop += 1\n",
    "    print(count_loop,'of',len(levels), end = ' ||| ')\n",
    "    for j in range(len(lats)):\n",
    "        for k in range(len(longs)):\n",
    "            temp = dtc.where(dtc.level == i)\n",
    "            temp = temp.where(temp.lat == lats[j]).where(temp.long == longs[k]).dropna(how='all')\n",
    "                \n",
    "            #temp0 = DTT.where(DTT.level == i)\n",
    "            #temp0 = temp0.where(temp0.lat == lats[j]).where(temp0.long == longs[k]).dropna(how='all')\n",
    "            \n",
    "            temp1 = dtd.where(dtd.level == i)\n",
    "            temp1 = temp1.where(temp1.lat == lats[j]).where(temp1.long == longs[k]).dropna(how='all')\n",
    "            \n",
    "            for timer in ind:\n",
    "                \n",
    "                Obs = temp.loc[timer]\n",
    "                #Nforc = temp0.loc[timer]\n",
    "                ONEforc = temp1.loc[timer]\n",
    "\n",
    "                #param = 'Wy'\n",
    "                #nforc = Nforc[param]\n",
    "                oneforc = ONEforc[param]\n",
    "                Obs = Obs[param]\n",
    "                \n",
    "                probability= .7\n",
    "                prob,va,vb = verify_res(nforc,'norm',probability,tole = 'None')\n",
    "                \n",
    "                Result_max = (Obs < vb)\n",
    "                Result_min = (Obs > va)\n",
    "                Result = (Result_max & Result_min)\n",
    "                \n",
    "                count += Result\n",
    "                sizet += 1\n",
    "                evalu += prob\n",
    "                \n",
    "                \n",
    "    print(str(time.time()-t_t)+ \" have passed since the start of the calcul.\")  \n",
    "    print(str( (time.time()-t_t)*(len(levels)-count_loop)/count_loop)+ \"s to finish\")  \n",
    "print(str(time.time()-t_t)+ \" seconds to calculate\")            \n",
    "evalu = evalu/sizet\n",
    "\n",
    "print('#','prob == ', probability, 'val =', count/sizet*100,'%','evalu=',evalu*100,'%')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 1\n",
    "prop = evalu\n",
    "print('Considered proportion:',prop)\n",
    "zval,pval = sm.stats.proportion.proportions_ztest(count,sizet*alpha,prop)\n",
    "print(count,'occurrences amid a total of',sizet*alpha)\n",
    "print('P-value of the test:',pval,'meaning the test is:',pval>0.05)"
   ]
  }
 ],
 "metadata": {
  "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}