fmda/presentations/fmda_kf_rnn_orig.ipynb

   1 {
   2  "cells": [
   3   {
   4    "cell_type": "code",
   5    "execution_count": null,
   6    "metadata": {},
   7    "outputs": [],
   8    "source": [
   9     "import os\n",
  10     "os.chdir('..')\n",
  11     "import reproducibility"
  12    ]
  13   },
  14   {
  15    "cell_type": "markdown",
  16    "metadata": {
  17     "id": "khOTxJsYc91W"
  18    },
  19    "source": [
  20     "# Kalman Filtering and Recurrent Neural Networks for Fuel Moisture\n",
  21     "## Jan Mandel, University of Colorado Denver\n"
  22    ]
  23   },
  24   {
  25    "cell_type": "markdown",
  26    "metadata": {
  27     "id": "sXaqfI-EdCEk"
  28    },
  29    "source": [
  30     "## Abstract"
  31    ]
  32   },
  33   {
  34    "cell_type": "markdown",
  35    "metadata": {
  36     "id": "ZbtSiYm4dF7B"
  37    },
  38    "source": [
  39     "''Fuel moisture is an important factor of the spread of wildland fires. Some weather stations have fuel moisture sensors and data are available online. We review a simple model of fuel moisture from atmospheric conditions, and show how to adjust the model using the weather station data."
  40    ]
  41   },
  42   {
  43    "cell_type": "markdown",
  44    "metadata": {
  45     "id": "eZ6dfHlZ63j1"
  46    },
  47    "source": [
  48     "## Table of contents"
  49    ]
  50   },
  51   {
  52    "cell_type": "markdown",
  53    "metadata": {
  54     "id": "WHIgN2uZ689b"
  55    },
  56    "source": [
  57     "1 Introduction\n",
  58     "\n",
  59     "2 Background\n",
  60     "\n",
  61     "2.1 Imports\n",
  62     "\n",
  63     "2.2 Kalman filter\n",
  64     "\n",
  65     "2.2.1 Overview\n",
  66     "\n",
  67     "2.2.2 Formulation\n",
  68     "\n",
  69     "2.2.3 A Kalman filter tester\n",
  70     "\n",
  71     "2.3 Fuel moisture model\n",
  72     "\n",
  73     "2.3.1 A simple time lag model\n",
  74     "\n",
  75     "2.3.1 Fuel moisture model with drying equilibrium, wetting equilibrium, and rain\n",
  76     "\n",
  77     "3 Methods\n",
  78     "\n",
  79     "3.1 Kalman filter demonstration on the simple model\n",
  80     "\n",
  81     "3.1.1 Creating synthetic data\n",
  82     "\n",
  83     "3.1.2 Running the Kalman filter\n",
  84     "\n",
  85     "3.2 Acquisition and preprocessing of real data\n",
  86     "\n",
  87     "3.2.1 Acquisition of fuel moisture observations\n",
  88     "\n",
  89     "3.2.2 Acquisition of weather data\n",
  90     "\n",
  91     "3.2.3 Preprocessing and visualization of the weather data\n",
  92     "\n",
  93     "4 Results\n",
  94     "\n",
  95     "4.1 Kalman filter with fuel moisture observations, followed by forecasting\n",
  96     "\n",
  97     "4.2 Model with an augmented state\n",
  98     "\n",
  99     "4.3 Kalman filter on the augmented model\n",
 100     "\n",
 101     "4.4 A comment on the information flow in the Kalman filter and in neural networks\n",
 102     "\n",
 103     "5. Conclusion\n",
 104     "\n",
 105     "Contributions of Authors\n",
 106     "\n",
 107     "Acknowledgements\n",
 108     "\n",
 109     "References\n",
 110     "\n",
 111     "\n",
 112     "\n",
 113     "\n",
 114     "\n",
 115     "\n",
 116     "\n"
 117    ]
 118   },
 119   {
 120    "cell_type": "markdown",
 121    "metadata": {
 122     "id": "ZFafUPCTO1N1"
 123    },
 124    "source": [
 125     "## 1 Introduction"
 126    ]
 127   },
 128   {
 129    "cell_type": "markdown",
 130    "metadata": {
 131     "id": "4_RcdWybPFks"
 132    },
 133    "source": [
 134     "The Kalman filter is at the foundation of many technologies in daily use, from GPS to weather forecasting. No model is completely accurate. Think space navigation: the movement of a Apollo 13 between the moon and the earth, subject to gravitational forces and propulsion, with the position ascertained by visual measurements. No matter how accurate the model of spacecraft motion is, the measurements are always burdened with noise. The idea of Kalman filter is to evolve a quantification of the of the state (here, positin and velocity of the spacecraft) in the form of a covariance matrix, and, using an estimate of the uncertainty of the data, adjust the state to split the difference every time measurements are taken. \n",
 135     "\n",
 136     "Here, we use the Kalman filter to estimate the evolution of fuel (dead wood) moisture content from a simple theoretical model, adjusting the state of the model hourly for measurements from fuel moisture a sensor in a wood stick exposed to the elements. This is needed for forecasting of wildfire progress; for this purpose, we also want to have the filter adjust the model from the data, so that it gives more accurate data for future when we only have hourly weather forecast but no actual data - because the future has not happened yet. "
 137    ]
 138   },
 139   {
 140    "cell_type": "markdown",
 141    "metadata": {
 142     "id": "M2kbwDPBTB7A"
 143    },
 144    "source": [
 145     "## 2 Background"
 146    ]
 147   },
 148   {
 149    "cell_type": "markdown",
 150    "metadata": {
 151     "id": "ar1BbXac49hO"
 152    },
 153    "source": [
 154     "In this section, we take care of preliminaries: we install some packages we need, and then proceed with the Kalman filter."
 155    ]
 156   },
 157   {
 158    "cell_type": "markdown",
 159    "metadata": {
 160     "id": "_5F5CuRqc91X"
 161    },
 162    "source": [
 163     "### 2.1 Imports"
 164    ]
 165   },
 166   {
 167    "cell_type": "markdown",
 168    "metadata": {
 169     "id": "K6sWUMf0c91Y"
 170    },
 171    "source": [
 172     "We may need the pygrib package to read weather data, but pygrib requires current numpy while google colab is using an old numpy version for compatibility with tensorflow. We will upgrade numpy and restart the runtime then the notebook will need to be run again. If numpy is current, we just download and import packages we need."
 173    ]
 174   },
 175   {
 176    "cell_type": "markdown",
 177    "metadata": {
 178     "id": "X9rvlymMZdJg"
 179    },
 180    "source": [
 181     "### 2.2 Kalman filter"
 182    ]
 183   },
 184   {
 185    "cell_type": "markdown",
 186    "metadata": {
 187     "id": "x5E2UE3F5gf2"
 188    },
 189    "source": [
 190     "#### 2.2.1 Overview"
 191    ]
 192   },
 193   {
 194    "cell_type": "markdown",
 195    "metadata": {
 196     "id": "NPgTHlCLAlA-"
 197    },
 198    "source": [
 199     "The Kalman filter provides an estimate $u$ of the time evolution of some unknown process, called \"nature\" or \"truth\". We do not know with certainty what the nature is, but we can observe it at regular intervals (steps) with some error. In each step, model $F$ advances the model state $u$ in time, $ u \\leftarrow F(u)$, and attempts to reconcile the state with an observation $d$ of the true state, so $u \\approx d$. The filter modifies the model state $u$ to balance the uncertainty in the model and the data (this is called *analysis*) and the cycle continues. For that purpose, the filter evolves also an estimate of the uncertainly of the model.\n",
 200     "\n",
 201     "More generally, instead of $u \\approx d$, only a part of the state is observed, and $Hu \\approx d$ where $H$ is a matrix, or observation function. Basically, $Hu$ is what the data would be if the model was completely accurate. \n",
 202     "\n",
 203     "See Kalman (1960) for the original publication, Kalnay (2003) for a gentle introduction, and the [Wikipedia article](https://en.wikipedia.org/wiki/Extended_Kalman_filter)."
 204    ]
 205   },
 206   {
 207    "cell_type": "markdown",
 208    "metadata": {
 209     "id": "y6j34L5s5pEL"
 210    },
 211    "source": [
 212     "#### 2.2.2 Formulation\n",
 213     "\n",
 214     "---\n",
 215     "\n"
 216    ]
 217   },
 218   {
 219    "cell_type": "markdown",
 220    "metadata": {
 221     "id": "b3GZW5vP5_o8"
 222    },
 223    "source": [
 224     "We present the Kalman filter in perhaps the most used form, as extended to nonlinear models.\n",
 225     " Consider a discrete time model of some natural\n",
 226     "process. At time step $k$, the model has state $u_{k}\\in\\mathbb{R}^{n}$, which\n",
 227     "can be approximated from the previous step $u_{k-1}$ by applying the model\n",
 228     "$\\mathcal{M}$ to get a forecast $u_{k}^{f}=\\mathcal{M}\\left(  u_{k-1}\\right)\n",
 229     "$. We model uncertainty in the model itself by adding normally distributed\n",
 230     "noise with mean zero and covariance $Q$ to the uncertainty of $u_{k}^{f}$. We\n",
 231     "also need to estimate now the uncertainty in the previous state $u_{k-1}$\n",
 232     "propagates to the uncertainty of the forecast $u_{k}^{f}$. So, assume that the\n",
 233     "model is differentiable and quantify the uncertainty of the state by a\n",
 234     "covariance matrix. That is,  assume that at step $k-1$, the state has\n",
 235     "(approximately) normal distribution with mean $u_{k-1}$ and covariance\n",
 236     "$P_{k-1}$. Using the Taylor expansion of order $1$ of the model operator at\n",
 237     "$u_{k-1}$, $\\mathcal{M}\\left(  u\\right)  \\approx\\mathcal{M}\\left(\n",
 238     "u_{k-1}\\right)  +\\mathcal{M}^{\\prime}\\left(  u_{k-1}\\right)  \\left(\n",
 239     "u-u_{k-1}\\right)  $, where $\\mathcal{M}^{\\prime}\\left(  u_{k-1}\\right)  $ is\n",
 240     "the Jacobian matrix of $\\mathcal{M}$ at $u_{k-1}$. It can be shown that the\n",
 241     "forecast has then (approximately)\\ normal distribution with mean and\n",
 242     "covariance\n",
 243     "$$\n",
 244     "u_{k}^{f}=\\mathcal{M}\\left(  u_{k-1}\\right)  ,\\ P_{k}^{f}=\\mathcal{M}\\left(\n",
 245     "u_{k-1}\\right)  P_{k-1}\\mathcal{M}^{\\prime}\\left(  u_{k-1}\\right)  +Q\n",
 246     "$$\n",
 247     "At time $k$, we also have an observation $d_{k}\\approx Hu_{k}$, where $H$ is a\n",
 248     "given observation operator, and we want to find $u_{k}$ so that both\n",
 249     "$$\n",
 250     "u_{k}\\approx u_{k}^{f}\\text{ and }d_{k}\\approx Hu_{k}.\n",
 251     "$$\n",
 252     "We quantify the uncertainly of the error of observation $d_{k}$ by a covariance\n",
 253     "matrix $R$: assume that the observation error has normal probability\n",
 254     "distribution with a known covariance $R$. Then, the likelihood of state $u$ is\n",
 255     "proportional to $e^{-\\left\\Vert d_{k}-Hu\\right\\Vert _{R^{-1}}^{2}/2}$, where\n",
 256     "we used the notation for the norm $\\left\\Vert v\\right\\Vert _{A}%\n",
 257     "=\\left(v^{\\top}Av\\right)^{1/2}$ induced by a positive definite matrix $A$. Similarly, we quantify the\n",
 258     "uncertainty of the state by a covariance matrix $P_{k}$. That is, the forecast\n",
 259     "state has (approximately) normal distribution with mean $u_{k}^{f}$  and covariance\n",
 260     "$P_{k}^{f}$. From the Bayes theorem of statistics, the probability distribution\n",
 261     "of the state after taking the data into account has density\n",
 262     "$$\n",
 263     "p_{k}\\left(  u\\right) \\propto e^\\frac{-\\left\\Vert d_{k}\n",
 264     "-Hu\\right\\Vert_{R^{-1}}^{2}}{2}e^\\frac{-\\left\\Vert u-u_{k}^{f}\\right\\Vert _{\n",
 265     "{P_{k}^f}^{-1}  }^{2}}{2}%\n",
 266     "$$\n",
 267     "where $\\propto$ means proportional.\n",
 268     "Note that the probability density at $u$ is maximal when $\\left\\Vert\n",
 269     "d_{k}-Hu\\right\\Vert _{R^{-1}}^{2}+\\left\\Vert u-u_{k}\\right\\Vert _{{P_{k}^{f}}^{-1}}^{2}$\n",
 270     " is minimal, which quantifies the statement that $d_{k}\\approx\n",
 271     "Hu_{k}$ and $u\\approx u_{k}^{f}$.  By a direct computation completing the\n",
 272     "square and using the Sherman-Morrison-Woodbury formula, \n",
 273     "$$p_{k}\\left(\n",
 274     "\t\tu\n",
 275     "\t   \\right) \\propto \n",
 276     "e^{-\\frac{\n",
 277     "\t\\left\\Vert u-u_{k\n",
 278     "\t         }\n",
 279     "\t\\right\\Vert_\n",
 280     "\t\t{P_{k\n",
 281     "\t\t      }^{-1}\n",
 282     "\t\t}^{2}\n",
 283     "\t}\n",
 284     "\t{2}},\n",
 285     "$$ \n",
 286     "which is the density of the normal distribution with the mean\n",
 287     "$$\n",
 288     "u_{k}^{f}=u_{k}^{f}+K_{k}(d-Hu_{k}^{f}),\\ \\text{where }K_{k}=P_{k}%\n",
 289     "^{f}H^{\\mathrm{T}}(HP_{k}^{f}H^{\\mathrm{T}}+R)^{-1}%\n",
 290     "$$\n",
 291     "and covariance\n",
 292     "$$\n",
 293     "P_{k}=\\left(  \\left(  P_{k}^{f}\\right)  ^{-1}+H^{\\mathrm{T}}R^{-1}H\\right)\n",
 294     "^{-1}=(I-KH)P_{k}^{f}.\n",
 295     "$$\n",
 296     "\n",
 297     "These are the equations of the extended Kalman filter. The original Kalman (1960) filter was\n",
 298     "formulated for a linear process. The extension to the\n",
 299     "nonlinear case made broad array of applications possible, including the Apollo spacecraft naviation (McGee and Schmidt, 1966),  and is\n",
 300     "still a de-facto standard in navigation and GPS.\n"
 301    ]
 302   },
 303   {
 304    "cell_type": "code",
 305    "execution_count": null,
 306    "metadata": {
 307     "id": "-bvUtJ_OLwQA"
 308    },
 309    "outputs": [],
 310    "source": [
 311     "import numpy as np\n",
 312     "def ext_kf(u,P,F,Q=0,d=None,H=None,R=None):\n",
 313     "  \"\"\"\n",
 314     "  One step of the extended Kalman filter. \n",
 315     "  If there is no data, only advance in time.\n",
 316     "  :param u:   the state vector, shape n\n",
 317     "  :param P:   the state covariance, shape (n,n)\n",
 318     "  :param F:   the model function, args vector u, returns F(u) and Jacobian J(u)\n",
 319     "  :param Q:   the process model noise covariance, shape (n,n)\n",
 320     "  :param d:   data vector, shape (m). If none, only advance in time\n",
 321     "  :param H:   observation matrix, shape (m,n)\n",
 322     "  :param R:   data error covariance, shape (n,n)\n",
 323     "  :return ua: the analysis state vector, shape (n)\n",
 324     "  :return Pa: the analysis covariance matrix, shape (n,n)\n",
 325     "  \"\"\"\n",
 326     "  def d2(a):\n",
 327     "    return np.atleast_2d(a) # convert to at least 2d array\n",
 328     "\n",
 329     "  def d1(a):\n",
 330     "    return np.atleast_1d(a) # convert to at least 1d array\n",
 331     "\n",
 332     "  # forecast\n",
 333     "  uf, J  = F(u)          # advance the model state in time and get the Jacobian\n",
 334     "  uf = d1(uf)            # if scalar, make state a 1D array\n",
 335     "  J = d2(J)              # if scalar, make jacobian a 2D array\n",
 336     "  P = d2(P)              # if scalar, make Jacobian as 2D array\n",
 337     "  Pf  = d2(J.T @ P) @ J + Q  # advance the state covariance Pf = J' * P * J + Q\n",
 338     "  # analysis\n",
 339     "  if d is None or not d.size :  # no data, no analysis\n",
 340     "    return uf, Pf\n",
 341     "  # K = P H' * inverse(H * P * H' + R) = (inverse(H * P * H' + R)*(H P))'\n",
 342     "  H = d2(H)\n",
 343     "  HP  = d2(H @ P)            # precompute a part used twice  \n",
 344     "  K   = d2(np.linalg.solve( d2(HP @ H.T) + R, HP)).T  # Kalman gain\n",
 345     "  # print('H',H)\n",
 346     "  # print('K',K)\n",
 347     "  res = d1(H @ d1(uf) - d)          # res = H*uf - d\n",
 348     "  ua = uf - K @ res # analysis mean uf - K*res\n",
 349     "  Pa = Pf - K @ d2(H @ P)        # analysis covariance\n",
 350     "  return ua, d2(Pa)\n"
 351    ]
 352   },
 353   {
 354    "cell_type": "markdown",
 355    "metadata": {
 356     "id": "Uvsbbv2XZ2Hd"
 357    },
 358    "source": [
 359     "#### 2.2.3 A Kalman filter tester"
 360    ]
 361   },
 362   {
 363    "cell_type": "markdown",
 364    "metadata": {
 365     "id": "gcmGBqPOU1e5"
 366    },
 367    "source": [
 368     "It is a very good idea to make write a simple tester for every piece of code. How else would we know it actually works, and that something basic did not get broken inadvertently, perhaps as a side effect of changing something else? A simple tester may save a great deal of time trying to debug cryptic errors later. And, what better place for a tester that right after the code it is testing so that it gets run every time?"
 369    ]
 370   },
 371   {
 372    "cell_type": "code",
 373    "execution_count": null,
 374    "metadata": {
 375     "id": "OsOqvQk6ZXZV"
 376    },
 377    "outputs": [],
 378    "source": [
 379     "# a basic ext_kf test\n",
 380     "import numpy as np\n",
 381     "u = [1,\n",
 382     "     2]\n",
 383     "P = [[2 , -1],\n",
 384     "    [-1 , 2]]\n",
 385     "A = [ [1 ,2],\n",
 386     "      [3 ,4]]\n",
 387     "u = np.array(u)      \n",
 388     "Q = np.array([[1,0],[0,1]])\n",
 389     "A = np.array(A)\n",
 390     "def fun(u):\n",
 391     "  return A @ u, A\n",
 392     "F = lambda u: fun(u)\n",
 393     "H = [[1, 0],\n",
 394     "     [0, 1]]\n",
 395     "d = [2,\n",
 396     "    3]\n",
 397     "R = [[2, 0],\n",
 398     "    [0, 2]]\n",
 399     "H = np.array(H)      \n",
 400     "d = np.array(d)\n",
 401     "R = np.array(R)\n",
 402     "ua,Pa = ext_kf(u,P,F,Q)\n",
 403     "print('ua=',ua)\n",
 404     "print('Pa=',Pa)\n",
 405     "ua,Pa = ext_kf(u,P,F,Q,d,H,R)\n",
 406     "print('ua=',ua)\n",
 407     "print('Pa=',Pa)\n"
 408    ]
 409   },
 410   {
 411    "cell_type": "markdown",
 412    "metadata": {
 413     "id": "A9ZpmNcdRpmp"
 414    },
 415    "source": [
 416     "### 2.3  Fuel moisture models\n",
 417     "\n",
 418     "\n"
 419    ]
 420   },
 421   {
 422    "cell_type": "markdown",
 423    "metadata": {
 424     "id": "eZL8gN7ISGVh"
 425    },
 426    "source": [
 427     "#### 2.3.1 A simple fuel moisture model"
 428    ]
 429   },
 430   {
 431    "cell_type": "markdown",
 432    "metadata": {
 433     "id": "1XvOC4kYSQgH"
 434    },
 435    "source": [
 436     "First consider a simplified fuel moisture model without considering the effect of rain.\n",
 437     "The evolution of fuel moisture content $m(t)$ is modeled by the time-lag differential equation on interval $\\left[\n",
 438     "t_{0},t_{1}\\right]  $,\n",
 439     "$$\n",
 440     "\\frac{dm}{dt}=\\frac{E-m(t)}{T},\\quad m(t_{0})=m_{0}.\n",
 441     "$$\n",
 442     "where the initial fuel moisture content $m_{0}=m\\left(  t_{0}\\right)  $ is the\n",
 443     "input, and $m_{1}=m(t_{1})$ is the output. Tnus, $m_1=F(m_0)$. The parameters of the model are the\n",
 444     "fuel moisture equilibrium $E$, assumed to be constant over the interval $\\left[\n",
 445     "t_{0},t_{1}\\right]  $, NS the characteristic decay time $T$. \n",
 446     "\n",
 447     "We can build the general model later by calling this simple model with different\n",
 448     "equilibria and time constants (drying, wetting, rain).\n",
 449     "\n",
 450     "Since $E$ is constant in time, the solution can be found\n",
 451     "analytically,\n",
 452     "$$\n",
 453     "m\\left(  t\\right)  =E+\\left(  m_{0}-E\\right)  e^{-t/T}%\n",
 454     "$$\n",
 455     "For convenience, we use $T_{1}=1/T$ instead of $T$, and the model becomes\n",
 456     "$$\n",
 457     "m_{1}=E+\\left(  m_{0}-E\\right)  e^{-\\left(  t_{1}-t_{0}\\right)  T_{1}}%\n",
 458     "$$\n",
 459     "In the extended Kalman filter, we will need the partial derivatives of $m_{1}$\n",
 460     "with respect to the input and the parameters. Compute\n",
 461     "$$\n",
 462     "\\frac{dm_{1}}{d_{m0}}=e^{-\\left(  t_{1}-t_{0}\\right)  T_{1}}\n",
 463     "$$\n",
 464     "$$\n",
 465     "\\frac{dm_{1}}{dE}=1-e^{-\\left(  t_{1}-t_{0}\\right)  T_{1}}\n",
 466     "$$\n",
 467     "$$\n",
 468     "\\frac{dm_{1}}{dT_{1}}=-\\left(  m_{0}-E\\right)  \\left(  t_{1}-t_{0}\\right)\n",
 469     "e^{-\\left(  t_{1}-t_{0}\\right)  T_{1}}\n",
 470     "$$\n",
 471     "At the moment, we need only ${dm_{1}}/{dm_{0}}$ but we put in the code all partials for possible use in future.\n"
 472    ]
 473   },
 474   {
 475    "cell_type": "code",
 476    "execution_count": null,
 477    "metadata": {
 478     "id": "wuVIAGLiSeR8"
 479    },
 480    "outputs": [],
 481    "source": [
 482     "import numpy as np\n",
 483     "def model_decay(m0,E,partials=0,T1=0.1,tlen=1):  \n",
 484     "  # Arguments: \n",
 485     "  #   m0          fuel moisture content at start dimensionless, unit (1)\n",
 486     "  #   E           fuel moisture eqilibrium (1)\n",
 487     "  #   partials=0: return m1 = fuel moisture contents after time tlen (1)\n",
 488     "  #           =1: return m1, dm0/dm0 \n",
 489     "  #           =2: return m1, dm1/dm0, dm1/dE\n",
 490     "  #           =3: return m1, dm1/dm0, dm1/dE dm1/dT1   \n",
 491     "  #   T1          1/T, where T is the time constant approaching the equilibrium\n",
 492     "  #               default 0.1/hour\n",
 493     "  #   tlen        the time interval length, default 1 hour\n",
 494     "\n",
 495     "  exp_t = np.exp(-tlen*T1)                  # compute this subexpression only once\n",
 496     "  m1 = E + (m0 - E)*exp_t                   # the solution at end\n",
 497     "  if partials==0:\n",
 498     "    return m1\n",
 499     "  dm1_dm0 = exp_t\n",
 500     "  if partials==1:\n",
 501     "    return m1, dm1_dm0          # return value and Jacobian\n",
 502     "  dm1_dE = 1 - exp_t      \n",
 503     "  if partials==2:\n",
 504     "     return m1, dm1_dm0, dm1_dE \n",
 505     "  dm1_dT1 = -(m0 - E)*tlen*exp_t            # partial derivative dm1 / dT1\n",
 506     "  if partials==3:\n",
 507     "    return m1, dm1_dm0, dm1_dE, dm1_dT1       # return value and all partial derivatives wrt m1 and parameters\n",
 508     "  raise('Bad arg partials')\n",
 509     "  "
 510    ]
 511   },
 512   {
 513    "cell_type": "markdown",
 514    "metadata": {
 515     "id": "dOARZlj-RUCi"
 516    },
 517    "source": [
 518     "#### 2.3.2 Fuel moisture model with drying equilibrium, wetting equilibrium, and rain"
 519    ]
 520   },
 521   {
 522    "cell_type": "markdown",
 523    "metadata": {
 524     "id": "AJp6FTpTSx5B"
 525    },
 526    "source": [
 527     "Here is a little more realistic fuel moisture model from Mandel et al. (2004). A rain-wetting lag time $t_{\\mathrm{r}}$ is reached for heavy rain only\n",
 528     "asymptotically, when the rain intensity $r$ (mm/h) is\n",
 529     "large:\n",
 530     "$$\n",
 531     "\\frac{\\mathrm{d}m}{\\mathrm{d}t}=\\frac{S-m}{t_{\\mathrm{r}}}\\left(1-\\exp\\left(-\\frac{r-r_0}{r_{\\mathrm{s}}}\n",
 532     "\\right)  \\right),\\ \\text{if}\\ r>r_0, \n",
 533     "$$\n",
 534     "where $r_0$ is the threshold rain intensity below which no perceptible\n",
 535     "wetting occurs, and $r_{\\mathrm{s}}$ is the saturation rain\n",
 536     "intensity. At the saturation rain intensity, $1-1/e\\approx 0.63$ of\n",
 537     "the maximal rain-wetting rate is achieved. For 10h fuel, the model takes $S=250\\,{\\%}$,\n",
 538     "$t_{\\mathrm{r}}=14$h, $r_0=0.05$mm/h and\n",
 539     "$r_{\\mathrm{s}}=8$mm/h. "
 540    ]
 541   },
 542   {
 543    "cell_type": "code",
 544    "execution_count": null,
 545    "metadata": {
 546     "id": "ITsKE0psRblG"
 547    },
 548    "outputs": [],
 549    "source": [
 550     "### Define model function with drying, wetting, and rain equilibria\n",
 551     "\n",
 552     "# Parameters\n",
 553     "r0 = 0.05                                   # threshold rainfall [mm/h]\n",
 554     "rs = 8.0                                    # saturation rain intensity [mm/h]\n",
 555     "Tr = 14.0                                   # time constant for rain wetting model [h]\n",
 556     "S = 250                                     # saturation intensity [dimensionless]\n",
 557     "T = 10.0                                    # time constant for wetting/drying\n",
 558     "\n",
 559     "def model_moisture(m0,Eqd,Eqw,r,t,partials=0,T=10.0,tlen=1.0):\n",
 560     "    # arguments:\n",
 561     "    # m0         starting fuel moistureb (%s\n",
 562     "    # Eqd        drying equilibrium      (%) \n",
 563     "    # Eqw        wetting equilibrium     (%)\n",
 564     "    # r          rain intensity          (mm/h)\n",
 565     "    # t          time\n",
 566     "    # partials = 0, 1, 2\n",
 567     "    # returns: same as model_decay\n",
 568     "    #   if partials==0: m1 = fuel moisture contents after time 1 hour\n",
 569     "    #              ==1: m1, dm1/dm0 \n",
 570     "    #              ==2: m1, dm1/dm0, dm1/dE  \n",
 571     "    \n",
 572     "    if r > r0:\n",
 573     "        # print('raining')\n",
 574     "        E = S\n",
 575     "        T1 =  (1.0 - np.exp(- (r - r0) / rs)) / Tr\n",
 576     "    elif m0 <= Eqw: \n",
 577     "        # print('wetting')\n",
 578     "        E=Eqw\n",
 579     "        T1 = 1.0/T\n",
 580     "    elif m0 >= Eqd:\n",
 581     "        # print('drying')\n",
 582     "        E=Eqd\n",
 583     "        T1 = 1.0/T\n",
 584     "    else: # no change'\n",
 585     "        E = m0\n",
 586     "        T1=0.0\n",
 587     "    exp_t = np.exp(-tlen*T1)\n",
 588     "    m1 = E + (m0 - E)*exp_t  \n",
 589     "    dm1_dm0 = exp_t\n",
 590     "    dm1_dE = 1 - exp_t\n",
 591     "    #if t>=933 and t < 940:\n",
 592     "    #  print('t,Eqw,Eqd,r,T1,E,m0,m1,dm1_dm0,dm1_dE',\n",
 593     "    #        t,Eqw,Eqd,r,T1,E,m0,m1,dm1_dm0,dm1_dE)   \n",
 594     "    if partials==0: \n",
 595     "        return m1\n",
 596     "    if partials==1:\n",
 597     "        return m1, dm1_dm0\n",
 598     "    if partials==2:\n",
 599     "        return m1, dm1_dm0, dm1_dE\n",
 600     "    raise('bad partials')"
 601    ]
 602   },
 603   {
 604    "cell_type": "markdown",
 605    "metadata": {
 606     "id": "JDLU3B_jV42l"
 607    },
 608    "source": [
 609     "## 3. Methods"
 610    ]
 611   },
 612   {
 613    "cell_type": "markdown",
 614    "metadata": {
 615     "id": "hLPJT3FcA2a7"
 616    },
 617    "source": [
 618     "### 3.1 Kalman filter demonstration on the simple model"
 619    ]
 620   },
 621   {
 622    "cell_type": "markdown",
 623    "metadata": {
 624     "id": "kIA3X8vluFdd"
 625    },
 626    "source": [
 627     "We demonstrate the Kalman filter for this model on a simple artificial example. The model is solving the differential equation for one hour. The equilibrium $E$ is constant during the hour, but it changes over the day so that it is higher at night and lower during the day, with a 24-hour period.  First, we create the \"truth\" by choosing the equilibrium $E$ and solving the differential aquation every hour, with a small additive noise. The synthetic data is obtained as the values of the \"truth\", with random noise to simulate observation error."
 628    ]
 629   },
 630   {
 631    "cell_type": "markdown",
 632    "metadata": {
 633     "id": "bBv10PTiChhm"
 634    },
 635    "source": [
 636     "#### 3.1.1 Creating synthetic data"
 637    ]
 638   },
 639   {
 640    "cell_type": "code",
 641    "execution_count": null,
 642    "metadata": {
 643     "id": "my6nnrk1iQo8"
 644    },
 645    "outputs": [],
 646    "source": []
 647   },
 648   {
 649    "cell_type": "code",
 650    "execution_count": null,
 651    "metadata": {
 652     "id": "-_pz-wXnCMnP"
 653    },
 654    "outputs": [],
 655    "source": [
 656     "def create_synthetic_data(days=20,power=4,data_noise=0.02,process_noise=0.0,DeltaE=0.0):\n",
 657     "  import numpy as np, random\n",
 658     "  hours = days*24\n",
 659     "  h2 = int(hours/2)\n",
 660     "  hour = np.array(range(hours))\n",
 661     "  day = np.array(range(hours))/24.\n",
 662     "\n",
 663     "  # artificial equilibrium data\n",
 664     "  E = np.power(np.sin(np.pi*day),4) # diurnal curve\n",
 665     "  E = 0.05+0.25*E\n",
 666     "  # FMC free run\n",
 667     "  m_f = np.zeros(hours)\n",
 668     "  m_f[0] = 0.1         # initial FMC\n",
 669     "  process_noise=0.\n",
 670     "  for t in range(hours-1):\n",
 671     "    m_f[t+1] = max(0.,model_decay(m_f[t],E[t])  + random.gauss(0,process_noise) )\n",
 672     "  data = m_f + np.random.normal(loc=0,scale=data_noise,size=hours)\n",
 673     "  E = E + DeltaE    \n",
 674     "\n",
 675     "  %matplotlib inline\n",
 676     "  import matplotlib.pyplot as plt \n",
 677     "  # fig1, ax1 = plt.subplots()\n",
 678     "\n",
 679     "  plt.figure(figsize=(16,4))\n",
 680     "  plt.plot(hour,E,linestyle='--',c='r',label='Equilibrium')\n",
 681     "  plt.plot(hour,m_f,linestyle='-',c='k',label='10-h fuel truth')\n",
 682     "  plt.scatter(hour[:h2],data[:h2],c='b',label='10-h fuel data')\n",
 683     "  plt.title('Synthetic data')\n",
 684     "  plt.xlabel('Time (hours)')\n",
 685     "  plt.ylabel('Fuel moisture content (%)')\n",
 686     "  plt.legend()\n",
 687     "  return E,m_f,data,hour,h2,DeltaE\n"
 688    ]
 689   },
 690   {
 691    "cell_type": "code",
 692    "execution_count": null,
 693    "metadata": {
 694     "id": "GfRxLp4HkcVz"
 695    },
 696    "outputs": [],
 697    "source": [
 698     "E,m_f,data,hour,h2,DeltaE = create_synthetic_data(days=20,power=4,data_noise=0.01,process_noise=0.0,DeltaE=0.0) "
 699    ]
 700   },
 701   {
 702    "cell_type": "markdown",
 703    "metadata": {
 704     "id": "z-3WLAEpD2yJ"
 705    },
 706    "source": [
 707     "#### 3.1.2 Running the Kalman filter"
 708    ]
 709   },
 710   {
 711    "cell_type": "markdown",
 712    "metadata": {
 713     "id": "T4g-RrrYAlBD"
 714    },
 715    "source": [
 716     "We have used the same code for model and for the truth, and run the Kalman filter for 10 days. The graph below shows that the model state was remarkably close to the truth, even if the model is fed only noisy observations. This is because the dynamics of the model and of the truth are the same. After 10 days, we let the model continue without any new data to simulate forecasting the future, and the agreement with the truth was still very good."
 717    ]
 718   },
 719   {
 720    "cell_type": "code",
 721    "execution_count": null,
 722    "metadata": {
 723     "id": "_-CjONZkD18n"
 724    },
 725    "outputs": [],
 726    "source": [
 727     "import numpy as np\n",
 728     "import matplotlib.pyplot as plt \n",
 729     "\n",
 730     "# using global E, m_f\n",
 731     "\n",
 732     "def plot_m(m,Ec=None,title=None,):  # global hour\n",
 733     "  hours=hour.shape[0]\n",
 734     "  %matplotlib inline\n",
 735     "  plt.figure(figsize=(16,4))\n",
 736     "  plt.plot(hour,E,linestyle='--',c='r',label='E=Equilibrium data')\n",
 737     "  # print(len(hour),len(m_f))\n",
 738     "  plt.plot(hour,m_f,linestyle='-',c='b',label='m_f=10-h fuel truth')\n",
 739     "  plt.scatter(hour[:h2],data[:h2],c='b',label='data=10-h fuel data')\n",
 740     "  if m is not None:\n",
 741     "    plt.plot(hour[:h2],m[:h2],linestyle='-',c='k',label='m=filtered')\n",
 742     "    plt.plot(hour[h2:hours],m[h2:hours],linestyle='-',c='r',label='m=forecast')\n",
 743     "  if Ec is not None:\n",
 744     "    plt.plot(hour,Ec,linestyle='-',c='g',label='Ec=Equilibrium correction')\n",
 745     "  if title is not None:\n",
 746     "    plt.title(title) \n",
 747     "  else:\n",
 748     "    plt.title('Kalman filtering and forecast on artificial data')\n",
 749     "  plt.xlabel('Time (hours)') \n",
 750     "  plt.ylabel('Fuel moisture content (%)')\n",
 751     "  plt.legend()\n",
 752     "\n",
 753     "def kf_example(DeltaE):\n",
 754     "  hours=hour.shape[0]\n",
 755     "  m = np.zeros(hours)\n",
 756     "  m[0]=0.1             # background state  \n",
 757     "  P = np.zeros(hours)\n",
 758     "  P[0] = 0.03 # background state variance\n",
 759     "  Q = np.array([0.02]) # process noise variance\n",
 760     "  H = np.array([1.])   # all observed\n",
 761     "  R = np.array([0.02]) # data variance\n",
 762     "\n",
 763     "  for t in range(h2):\n",
 764     "    # use lambda construction to pass additional arguments to the model \n",
 765     "    m[t+1],P[t+1] = ext_kf(m[t],P[t],lambda u: model_decay(u,E[t]+DeltaE,partials=1),Q,\n",
 766     "                    d=data[t],H=H,R=R)\n",
 767     "  for t in range(h2,hours - 1):\n",
 768     "    m[t+1],P[t+1] = ext_kf(m[t],P[t],lambda u: model_decay(u,E[t]+DeltaE,partials=1))\n",
 769     "  return m, P\n",
 770     " "
 771    ]
 772   },
 773   {
 774    "cell_type": "code",
 775    "execution_count": null,
 776    "metadata": {
 777     "id": "d0EFhTPZAlBD",
 778     "scrolled": true
 779    },
 780    "outputs": [],
 781    "source": [
 782     "DeltaE = 0.0          # bias\n",
 783     "m, P = kf_example(DeltaE)\n",
 784     "plot_m(m)\n"
 785    ]
 786   },
 787   {
 788    "cell_type": "markdown",
 789    "metadata": {
 790     "id": "vqyB2Yz3uCsD"
 791    },
 792    "source": [
 793     "We have recovered the fuel moisture from data with random noise - we **filtered** the noise out. "
 794    ]
 795   },
 796   {
 797    "cell_type": "markdown",
 798    "metadata": {
 799     "id": "Dl7pBZ9B3Nox"
 800    },
 801    "source": [
 802     "Let's have a look at the evolution of the filter's estimate of the variance $P$. A common problem with the Kalman filter is when the variance converges to zero over time, then, since the filter trusts the model too much, it ignores the observations. Of course, once we switch to forecasting mode, the variance is not of interest. We could keep evolving the variance to bridge over periods when there are no observations, but not in this simplified version."
 803    ]
 804   },
 805   {
 806    "cell_type": "code",
 807    "execution_count": null,
 808    "metadata": {
 809     "id": "wRJgbmGLc91g"
 810    },
 811    "outputs": [],
 812    "source": [
 813     "%matplotlib inline\n",
 814     "plt.figure(figsize=(16,4))\n",
 815     "plt.plot(P,linestyle='-',c='b',label='Estimated state variance P')\n",
 816     "plt.title('Kalman filtering and forecast on artificial data')\n",
 817     "plt.xlabel('Time (hours)') \n",
 818     "plt.ylabel('Estimated variance of fuel moisture (%^2)')\n",
 819     "plt.legend()"
 820    ]
 821   },
 822   {
 823    "cell_type": "markdown",
 824    "metadata": {
 825     "id": "Ccr-uKbmAlBE"
 826    },
 827    "source": [
 828     "Now what if the model is wrong - different from nature? That is always so in reality. Now suppose that the model and the truth are not the same. That is always the case in reality.  Consider a simple case when the model thinks that the equilibrium $E$ is too high."
 829    ]
 830   },
 831   {
 832    "cell_type": "code",
 833    "execution_count": null,
 834    "metadata": {
 835     "id": "spMdGW8oAlBE"
 836    },
 837    "outputs": [],
 838    "source": [
 839     "DeltaE = -0.05\n",
 840     "m, P = kf_example(DeltaE)\n",
 841     "plot_m(m)"
 842    ]
 843   },
 844   {
 845    "cell_type": "code",
 846    "execution_count": null,
 847    "metadata": {
 848     "id": "Z_INvPNGCmM2"
 849    },
 850    "outputs": [],
 851    "source": [
 852     "DeltaE = 0.05\n",
 853     "m, P = kf_example(DeltaE)\n",
 854     "plot_m(m)"
 855    ]
 856   },
 857   {
 858    "cell_type": "markdown",
 859    "metadata": {
 860     "id": "DQeF7J8T4j2i"
 861    },
 862    "source": [
 863     "We have found a good estimate of the state $m$, while data is available. Also, the estimated state variance $P$ converges with time - we have *learned* the variance that balances the noise. But for forecasting fuel moisture, we need to continue the fuel moisture model into the future, and we can't have any measurements from future. We only have the equilibrium from weather forecast. And the forecast and the truth disagree - as soon as there is no data to attract the simulation, the model is doing its own thing."
 864    ]
 865   },
 866   {
 867    "cell_type": "markdown",
 868    "metadata": {
 869     "id": "-WMWCDz4DX45"
 870    },
 871    "source": [
 872     "#### 3.2 Model with an augmented state"
 873    ]
 874   },
 875   {
 876    "cell_type": "markdown",
 877    "metadata": {
 878     "id": "jivOYEhiXMi5"
 879    },
 880    "source": [
 881     "In reality, the equilibrium moisture $E$ computed from atmospheric conditions\n",
 882     "generally does not agree with the data. We want to add a correction $\\Delta\n",
 883     "E$ to $E$ constant in time, and identify the new parameter $\\Delta E$ from data. \n",
 884     "Because the Kalman filter identifies state, add the parameter to the state.\n",
 885     "Define augmented state $u=\\left[\n",
 886     "\\begin{array}\n",
 887     "[c]{c}\n",
 888     "m\\\\\n",
 889     "\\Delta E\n",
 890     "\\end{array}\n",
 891     "\\right]  .$ Since $\\Delta E$ is constant in time, it satisfies the\n",
 892     "differential equation $\\frac{d\\Delta E}{dt}=0.$ So, we want to estimate the\n",
 893     "state $u$ governed by the\n",
 894     "$$\n",
 895     "\\frac{d}{dt}\\left[\n",
 896     "\\begin{array}\n",
 897     "[c]{c}\n",
 898     "m\\\\\n",
 899     "\\Delta E\n",
 900     "\\end{array}\n",
 901     "\\right]  =\\left[\n",
 902     "\\begin{array}\n",
 903     "[c]{c}\n",
 904     "\\frac{E+\\Delta E-m(t)}{T}\\\\\n",
 905     "0\n",
 906     "\\end{array}\n",
 907     "\\right]  ,\n",
 908     "$$\n",
 909     "which we write as $\\frac{du}{dt}=F(u),$ where\n",
 910     "$$\n",
 911     "F(u)=\\left[\n",
 912     "\\begin{array}\n",
 913     "[c]{c}\n",
 914     "F_{1}\\left(  u\\right)  \\\\\n",
 915     "F_{2}\\left(  u\\right)\n",
 916     "\\end{array}\n",
 917     "\\right]  =F\\left(  \\left[\n",
 918     "\\begin{array}\n",
 919     "[c]{c}\n",
 920     "m\\\\\n",
 921     "\\Delta E\n",
 922     "\\end{array}\n",
 923     "\\right]  \\right)  =\\left[\n",
 924     "\\begin{array}\n",
 925     "[c]{c}\n",
 926     "\\left(  E+\\Delta E-m(t)\\right)  T_{1}\\\\\n",
 927     "0\n",
 928     "\\end{array}\n",
 929     "\\right]  ,\\quad T_{1}=\\frac{1}{T}.\n",
 930     "$$\n",
 931     "The Jacobian of $F$ is\n",
 932     "$$\n",
 933     "\\left[\n",
 934     "\\begin{array}\n",
 935     "[c]{cc}\n",
 936     "\\frac{\\partial F_{1}}{\\partial u_{1}} & \\frac{\\partial F_{1}}{\\partial u_{2}\n",
 937     "}\\\\\n",
 938     "\\frac{\\partial F_{2}}{\\partial u_{1}} & \\frac{\\partial F_{2}}{\\partial u_{2}}\n",
 939     "\\end{array}\n",
 940     "\\right]  =\\left[\n",
 941     "\\begin{array}\n",
 942     "[c]{cc}\n",
 943     "\\frac{\\partial m_{1}}{\\partial m_{0}} & \\frac{\\partial m_{1}}{\\partial E}\\\\\n",
 944     "\\frac{\\partial\\Delta E}{\\partial m_{0}} & \\frac{\\partial\\Delta E}\n",
 945     "{\\partial\\Delta E}\n",
 946     "\\end{array}\n",
 947     "\\right]  =\\left[\n",
 948     "\\begin{array}\n",
 949     "[c]{cc}\n",
 950     "\\frac{\\partial m_{1}}{\\partial m_{0}} & \\frac{\\partial m_{1}}{\\partial E}\\\\\n",
 951     "0 & 1\n",
 952     "\\end{array}\n",
 953     "\\right]\n",
 954     "$$\n",
 955     "Here is a function that implements the augmented model $F$. The input is\n",
 956     "$u_{0}$. The output is $u_{1}$ and the Jacobian $du_{1}/du_{0}$."
 957    ]
 958   },
 959   {
 960    "cell_type": "markdown",
 961    "metadata": {
 962     "id": "MJ1C_1Omc91s"
 963    },
 964    "source": [
 965     "\n",
 966     "Define augmented model function. Also, add use drying, wetting, and rain equilibria"
 967    ]
 968   },
 969   {
 970    "cell_type": "code",
 971    "execution_count": null,
 972    "metadata": {
 973     "id": "GHtAaGp9WSHT"
 974    },
 975    "outputs": [],
 976    "source": [
 977     "def model_augmented(u0,Ed,Ew,r,t):\n",
 978     "    # state u is the vector [m,dE] with dE correction to equilibria Ed and Ew at t\n",
 979     "    # \n",
 980     "    m0, Ec = u0  # decompose state u0\n",
 981     "    # reuse model_moisture(m0,Eqd,Eqw,r,partials=0):\n",
 982     "    # arguments:\n",
 983     "    # m0         starting fuel moistureb (1)\n",
 984     "    # Ed         drying equilibrium      (1) \n",
 985     "    # Ew         wetting equilibrium     (1)\n",
 986     "    # r          rain intensity          (mm/h)\n",
 987     "    # partials = 0, 1, 2\n",
 988     "    # returns: same as model_decay\n",
 989     "    #   if partials==0: m1 = fuel moisture contents after time 1 hour\n",
 990     "    #              ==1: m1, dm0/dm0 \n",
 991     "    #              ==2: m1, dm1/dm0, dm1/dE \n",
 992     "    m1, dm1_dm0, dm1_dE  = model_moisture(m0,Ed + Ec, Ew + Ec, r, t, partials=2)\n",
 993     "    u1 = np.array([m1,Ec])   # dE is just copied\n",
 994     "    J =  np.array([[dm1_dm0, dm1_dE],\n",
 995     "                   [0.     ,     1.]])\n",
 996     "    return u1, J"
 997    ]
 998   },
 999   {
1000    "cell_type": "code",
1001    "execution_count": null,
1002    "metadata": {
1003     "id": "1No3g6HyAEh_"
1004    },
1005    "outputs": [],
1006    "source": [
1007     "def run_augmented_kf(d,Ed,Ew,rain,h2,hours):\n",
1008     "  u = np.zeros((2,hours))\n",
1009     "  u[:,0]=[0.1,0.0]       # initialize,background state  \n",
1010     "  P = np.zeros((2,2,hours))\n",
1011     "  P[:,:,0] = np.array([[1e-3, 0.],\n",
1012     "                      [0.,  1e-3]]) # background state covariance\n",
1013     "  Q = np.array([[1e-3, 0.],\n",
1014     "                [0,  1e-3]]) # process noise covariance\n",
1015     "  H = np.array([[1., 0.]])  # first component observed\n",
1016     "  R = np.array([1e-3]) # data variance\n",
1017     "\n",
1018     "  # ext_kf(u,P,F,Q=0,d=None,H=None,R=None) returns ua, Pa\n",
1019     "\n",
1020     "  # print('initial u=',u,'P=',P)\n",
1021     "  # print('Q=',Q,'H=',H,'R=',R)\n",
1022     "\n",
1023     "  for t in range(1,h2):\n",
1024     "      # use lambda construction to pass additional arguments to the model \n",
1025     "      u[:,t],P[:,:,t] = ext_kf(u[:,t-1],P[:,:,t-1],\n",
1026     "                                  lambda uu: model_augmented(uu,Ed[t],Ew[t],rain[t],t),\n",
1027     "                                  Q,d[t],H=H,R=R)\n",
1028     "      # print('time',t,'data',d[t],'filtered',u[0,t],'Ec',u[1,t])\n",
1029     "  for t in range(h2,hours):\n",
1030     "      u[:,t],P[:,:,t] = ext_kf(u[:,t-1],P[:,:,t-1],\n",
1031     "                                  lambda uu: model_augmented(uu,Ed[t],Ew[t],rain[t],t),\n",
1032     "                                  Q*0.0)\n",
1033     "      # print('time',t,'data',d[t],'forecast',u[0,t],'Ec',u[1,t])\n",
1034     "  return u"
1035    ]
1036   },
1037   {
1038    "cell_type": "code",
1039    "execution_count": null,
1040    "metadata": {
1041     "id": "9vUv1Uc5OfPI"
1042    },
1043    "outputs": [],
1044    "source": [
1045     "def augmented_example(DeltaE):\n",
1046     "  hours=hour.shape[0]\n",
1047     "  h2 = int(hours/2)\n",
1048     "  m, Ec = run_augmented_kf(data,E+DeltaE,E+DeltaE,0*E,h2,hours)  # data, E, hours are global\n",
1049     "  return m, Ec"
1050    ]
1051   },
1052   {
1053    "cell_type": "code",
1054    "execution_count": null,
1055    "metadata": {
1056     "id": "QQnRlereDHbf"
1057    },
1058    "outputs": [],
1059    "source": [
1060     "m, Ec=augmented_example(0.1)\n",
1061     "plot_m(m, Ec)"
1062    ]
1063   },
1064   {
1065    "cell_type": "code",
1066    "execution_count": null,
1067    "metadata": {
1068     "id": "ZOC6bplsIRCx"
1069    },
1070    "outputs": [],
1071    "source": [
1072     "m, Ec=augmented_example(0.0)\n",
1073     "plot_m(m, Ec)"
1074    ]
1075   },
1076   {
1077    "cell_type": "code",
1078    "execution_count": null,
1079    "metadata": {
1080     "id": "xuMx41fXIejB"
1081    },
1082    "outputs": [],
1083    "source": [
1084     "m, Ec=augmented_example(-0.1)\n",
1085     "plot_m(m, Ec)"
1086    ]
1087   },
1088   {
1089    "cell_type": "markdown",
1090    "metadata": {
1091     "id": "ejlc1sbENm6R"
1092    },
1093    "source": [
1094     "## From Kalman filter to neural networks"
1095    ]
1096   },
1097   {
1098    "cell_type": "markdown",
1099    "metadata": {
1100     "id": "6zeGArbmNkAW"
1101    },
1102    "source": [
1103     "In the Kalman filter, at each time step $k$,\n",
1104     "\n",
1105     "* the input state is $u_{k-1}$ size $n$ and its covariance matrix $P_{k-1}$ size $n \\times n$.\n",
1106     "* the model is applied to external data $e_k$ and the input $u_{k-1},P_{k-1}$ produce the forecast $u_k^f$ and its covariance $P^f_k$\n",
1107     "* the new state $u_k$ is found by minimizing $|| u^f_k - u_k||^2_{P^f_k} + ||H u_k - d_k||^2_{R}$   \n",
1108     "* the new state covariance is $P_k = ( (P^f_k)^{-1} + H^\\top R^{-1} H)^{-1}$.\n",
1109     "\n",
1110     "Here, the state consists of \n",
1111     "* the fuel moisture and the adjustment to the equilibrium, dimension 2\n",
1112     "* the covariance matrix of vector of dimension 2, which is symmetric $2 \\times 2$ matrix, given by 3 numbers because it is symmetric\n",
1113     "Thus, the dimension of the state is 2 + 3 = 5. The first component of the state, the fuel moisture, is the quantity of interest, the rest are auxiliary.\n",
1114     "\n",
1115     "\n",
1116     "This can be understood as:\n",
1117     "\n",
1118     "* a mapping $M$ of the 5-dimensional hidden and external data state to a new hidden state:\n",
1119     "$$M:(u_{k-1},P_{k-1},e_k) \\mapsto (u_{k},P_{k})$$\n",
1120     "* retrieving the output (the quantity of interest) as the first component of the hiddent state\n",
1121     "* feeding the hiddent state back to the mapping $M$ for the next step $k+1$\n",
1122     "* training consists of fitting the hidden state to minimize a loss function\n",
1123     "$$\\ell(u_{k},P_{k},d_k,R_k) \\to \\min$$\n",
1124     "\n",
1125     "Note that in the augmented Kalman filter above, the mapping $M$ is fixed and it has a one component of the hidden state as a parameter. To get a better fit, we could increase the number of parameters, e.g., by modeling the moisture in multiple layers, as in van der Kamp et al. (2017) two-layer model.\n"
1126    ]
1127   },
1128   {
1129    "cell_type": "markdown",
1130    "metadata": {
1131     "id": "fk72YB2mjuGk"
1132    },
1133    "source": [
1134     "Building and evaluating RNN"
1135    ]
1136   },
1137   {
1138    "cell_type": "markdown",
1139    "metadata": {
1140     "id": "svt8wUAsSA67"
1141    },
1142    "source": [
1143     "A recurrent neural network (RNN) has a similar information flow but it can be more flexible and look for the best model automatically, i.e., build the model from data. \n",
1144     "\n",
1145     "We'll start by how to evaluate the map, then actually create it later.\n",
1146     "\n",
1147     "Some of the code is from https://machinelearningmastery.com/understanding-simple-recurrent-neural-networks-in-keras/\n"
1148    ]
1149   },
1150   {
1151    "cell_type": "code",
1152    "execution_count": null,
1153    "metadata": {
1154     "id": "H3RTQCDG9q-4"
1155    },
1156    "outputs": [],
1157    "source": [
1158     "import numpy as np\n",
1159     "import tensorflow as tf\n",
1160     "from keras.models import Sequential\n",
1161     "from keras.layers import Dense, SimpleRNN\n",
1162     "from keras.utils.vis_utils import plot_model\n",
1163     "from sklearn.preprocessing import MinMaxScaler\n",
1164     "from sklearn.metrics import mean_squared_error\n",
1165     "import math\n",
1166     "import matplotlib.pyplot as plt\n",
1167     "import tensorflow as tf\n",
1168     "import keras.backend as K"
1169    ]
1170   },
1171   {
1172    "cell_type": "code",
1173    "execution_count": null,
1174    "metadata": {
1175     "id": "pcIU5lWhVEAy"
1176    },
1177    "outputs": [],
1178    "source": [
1179     "def create_RNN(hidden_units, dense_units, input_shape, activation):\n",
1180     "    inputs = tf.keras.Input(shape=input_shape)\n",
1181     "    # https://stackoverflow.com/questions/43448029/how-can-i-print-the-values-of-keras-tensors\n",
1182     "    # inputs2 = K.print_tensor(inputs, message='inputs = ')  # change allso inputs to inputs2 below, must be used\n",
1183     "    x = tf.keras.layers.SimpleRNN(hidden_units, input_shape=input_shape,\n",
1184     "                        activation=activation[0])(inputs)\n",
1185     "    outputs = tf.keras.layers.Dense(dense_units, activation=activation[1])(x)\n",
1186     "    model = tf.keras.Model(inputs=inputs, outputs=outputs)\n",
1187     "    model.compile(loss='mean_squared_error', optimizer='adam')\n",
1188     "    return model"
1189    ]
1190   },
1191   {
1192    "cell_type": "code",
1193    "execution_count": null,
1194    "metadata": {
1195     "id": "lphzeucMfI7L"
1196    },
1197    "outputs": [],
1198    "source": [
1199     "# Demo example\n",
1200     "hidden=5\n",
1201     "features=2\n",
1202     "timesteps=3\n",
1203     "demo_model = create_RNN(hidden_units=hidden, dense_units=1, \n",
1204     "                        input_shape=(timesteps,features), \n",
1205     "                        activation=['linear', 'linear'])\n",
1206     "print(demo_model.summary())\n",
1207     "w = demo_model.get_weights()\n",
1208     "#print(len(w),' weight arrays:',w)\n",
1209     "wname=('wx','wh','bh','wy','by','wz','bz')\n",
1210     "for i in range(len(w)):\n",
1211     "  print(i,':',wname[i],'shape=',w[i].shape)\n",
1212     "wx, wh, bh, wy, by = w\n",
1213     "plot_model(demo_model, to_file='model_plot.png', \n",
1214     "  show_shapes=True, show_layer_names=True,\n",
1215     "  expand_nested=True,)"
1216    ]
1217   },
1218   {
1219    "cell_type": "markdown",
1220    "metadata": {
1221     "id": "uFVQdXL0fODX"
1222    },
1223    "source": [
1224     "The input layer here is just a formality. The input of the hidden layer `simple_rnn` consist of vector passed by the input layer, followed by its own output from the previous time step.\n",
1225     "\n",
1226     "Now let’s do a simple experiment to see how the layers from a SimpleRNN and Dense layer produce an output. Keep this figure in view.\n",
1227     "<img src=\"https://machinelearningmastery.com/wp-content/uploads/2021/09/rnnCode1.png\">"
1228    ]
1229   },
1230   {
1231    "cell_type": "markdown",
1232    "metadata": {
1233     "id": "TcxvQPo1hYip"
1234    },
1235    "source": [
1236     "We’ll input x for three time steps and let the network generate an output. The values of the hidden units at time steps 1, 2 and 3 will be computed. $h(0)$ is initialized to the zero vector. The output $o(3)$ is computed from $h(3)$ and $w(3)$. An activation function is linear, $f(x)=x$, so the update of  $h(k)$  and the output $o(k)$ are given by\n",
1237     "\\begin{align*}\n",
1238     "h\\left(  0\\right)  = &0  \\\\\n",
1239     "h\\left(  k+1\\right)  =& \n",
1240     "x\\left(  k\\right) w_{x}\n",
1241     "  +h(k) w_{h}  + b_{h}\\\\\n",
1242     "o(k+1)=& h(k+1)w_{y} + b_y\n",
1243     "\\end{align*}"
1244    ]
1245   },
1246   {
1247    "cell_type": "code",
1248    "execution_count": null,
1249    "metadata": {
1250     "id": "fqL5TEfpml7q"
1251    },
1252    "outputs": [],
1253    "source": [
1254     "# Reshape the input to sample_size x time_steps x features \n",
1255     "samples=4   # number of samples\n",
1256     "x = tf.reshape(tf.range(samples*timesteps*features),[samples,timesteps,features]) \n",
1257     "print('test input x=',x)\n",
1258     "print('model.predict start')\n",
1259     "y_pred_model = demo_model.predict(x)\n",
1260     "print('model.predict end')\n",
1261     "\n",
1262     "o3=np.zeros([samples,1])\n",
1263     "for i in range(samples):\n",
1264     "  h_0 = np.zeros(hidden)\n",
1265     "  h_1 = np.dot(x[i,0,:], wx) + np.dot(h_0,wh) + bh\n",
1266     "  h_2 = np.dot(x[i,1,:], wx) + np.dot(h_1,wh) + bh\n",
1267     "  h_3 = np.dot(x[i,2,:], wx) + np.dot(h_2,wh) + bh\n",
1268     "  o3[i,0] = np.dot(h_3, wy) + by\n",
1269     "#print('h1 = ', h_1,'h2 = ', h_2,'h3 = ', h_3)\n",
1270     "\n",
1271     "print(\"Prediction from network \", y_pred_model)\n",
1272     "print(\"Prediction from our computation \", o3)"
1273    ]
1274   },
1275   {
1276    "cell_type": "markdown",
1277    "metadata": {
1278     "id": "6qdqOFCvhQL1"
1279    },
1280    "source": [
1281     "The result is the same."
1282    ]
1283   },
1284   {
1285    "cell_type": "markdown",
1286    "metadata": {
1287     "id": "AkyiGlZF0WrM"
1288    },
1289    "source": [
1290     "#### Training and forecasting with the RNN"
1291    ]
1292   },
1293   {
1294    "cell_type": "markdown",
1295    "metadata": {
1296     "id": "e265QFwlw22_"
1297    },
1298    "source": [
1299     "We are given a sequence `x` of inputs size `[train_steps+forecast_steps,features]` and want to train a model so that at step `i` in `range(train_steps)`, the model returns close to `features[i,:]`. The trained model then returns for `i` in `range(train_steps,train_steps+forecast_steps)` a forecast `features[i,:]`."
1300    ]
1301   },
1302   {
1303    "cell_type": "code",
1304    "execution_count": null,
1305    "metadata": {
1306     "id": "owX7OTg-RnMA"
1307    },
1308    "outputs": [],
1309    "source": [
1310     "def staircase(x,y,timesteps,trainsteps,return_sequences=False):\n",
1311     "  # x [trainsteps+forecaststeps,features]    all inputs\n",
1312     "  # y [trainsteps,outputs]\n",
1313     "  # timesteps: split x and y into samples length timesteps, shifted by 1\n",
1314     "  # trainsteps: number of timesteps to use for training, no more than y.shape[0]\n",
1315     "  print('shape x = ',x.shape)\n",
1316     "  print('shape y = ',y.shape)\n",
1317     "  print('timesteps=',timesteps)\n",
1318     "  print('trainsteps=',trainsteps)\n",
1319     "  outputs = y.shape[1]\n",
1320     "  features = x.shape[1]\n",
1321     "  forecaststeps = x.shape[0]-trainsteps\n",
1322     "  samples = trainsteps-timesteps+1\n",
1323     "  print('staircase: samples=',samples,'timesteps=',timesteps,'features=',features)\n",
1324     "  x_train = np.empty([samples, timesteps, features])\n",
1325     "  print('return_sequences=',return_sequences)\n",
1326     "  if return_sequences:\n",
1327     "    print('returning all timesteps in a sample')\n",
1328     "    y_train = np.empty([samples, timesteps, outputs])  # all\n",
1329     "    for i in range(samples):\n",
1330     "      for k in range(timesteps):\n",
1331     "        for j in range(features):\n",
1332     "          x_train[i,k,j] = x[i+k,j]\n",
1333     "        for j in range(outputs):\n",
1334     "          y_train[i,k,j] = y[i+k,j]\n",
1335     "  else:\n",
1336     "    print('returning only the last timestep in a sample')\n",
1337     "    y_train = np.empty([samples, outputs])\n",
1338     "    for i in range(samples):\n",
1339     "      for j in range(features):\n",
1340     "        for k in range(timesteps):\n",
1341     "          x_train[i,k,j] = x[i+k,j]\n",
1342     "      for j in range(outputs):\n",
1343     "        y_train[i,j] = y[i+timesteps-1,j]\n",
1344     "\n",
1345     "  return x_train, y_train"
1346    ]
1347   },
1348   {
1349    "cell_type": "code",
1350    "execution_count": null,
1351    "metadata": {
1352     "id": "FzOotSFf-tPR"
1353    },
1354    "outputs": [],
1355    "source": [
1356     "def seq2batches(x,y,timesteps,trainsteps):\n",
1357     "  # x [trainsteps+forecaststeps,features]    all inputs\n",
1358     "  # y [trainsteps,outputs]\n",
1359     "  # timesteps: split x and y into samples length timesteps, shifted by 1\n",
1360     "  # trainsteps: number of timesteps to use for training, no more than y.shape[0]\n",
1361     "  print('shape x = ',x.shape)\n",
1362     "  print('shape y = ',y.shape)\n",
1363     "  print('timesteps=',timesteps)\n",
1364     "  print('trainsteps=',trainsteps)\n",
1365     "  outputs = y.shape[1]\n",
1366     "  features = x.shape[1]\n",
1367     "  samples= trainsteps - timesteps + 1\n",
1368     "  print('samples=',samples)\n",
1369     "  x_train = np.empty([samples, timesteps, features])\n",
1370     "  y_train = np.empty([samples, timesteps, outputs])  # only the last\n",
1371     "  print('samples=',samples,' timesteps=',timesteps,\n",
1372     "        ' features=',features,' outputs=',outputs)\n",
1373     "  for i in range(samples):\n",
1374     "    for k in range(timesteps):\n",
1375     "      for j in range(features):\n",
1376     "        x_train[i,k,j] = x[i+k,j]\n",
1377     "      for j in range(outputs):\n",
1378     "        y_train[i,k,j] = y[i+k,j]  # return sequences\n",
1379     "  return x_train, y_train"
1380    ]
1381   },
1382   {
1383    "cell_type": "code",
1384    "execution_count": null,
1385    "metadata": {
1386     "id": "Kg7wSrkk-HrE"
1387    },
1388    "outputs": [],
1389    "source": [
1390     "print('test preprocessing for RNN')\n",
1391     "trainsteps=5\n",
1392     "features=1\n",
1393     "outputs=1\n",
1394     "timesteps=3\n",
1395     "x = tf.reshape(tf.range(trainsteps*features),[trainsteps,features])\n",
1396     "y = tf.reshape(tf.range(trainsteps*outputs),[trainsteps,outputs])\n",
1397     "print('x=',x)\n",
1398     "print('y=',y)\n",
1399     "x_train, y_train = staircase(x,y,timesteps,trainsteps)\n",
1400     "print('x_train=',x_train)\n",
1401     "print('y_train=',y_train)\n",
1402     "x_train, y_train = seq2batches(x,y,timesteps,trainsteps)\n",
1403     "print('x_train=',x_train)\n",
1404     "print('y_train=',y_train)"
1405    ]
1406   },
1407   {
1408    "cell_type": "code",
1409    "execution_count": null,
1410    "metadata": {
1411     "id": "rHWxqARpSO_f"
1412    },
1413    "outputs": [],
1414    "source": [
1415     "E,m_f,data,hour,h2,DeltaE = create_synthetic_data(days=20,power=4,data_noise=0.01,process_noise=0.0,DeltaE=0.1) "
1416    ]
1417   },
1418   {
1419    "cell_type": "code",
1420    "execution_count": null,
1421    "metadata": {
1422     "id": "vfcxN9JCH5Ku"
1423    },
1424    "outputs": [],
1425    "source": [
1426     "scale=False\n",
1427     "# transform as 2D, (timesteps, features) and (timesteps, outputs)\n",
1428     "Et = np.reshape(E,[E.shape[0],1])\n",
1429     "datat = np.reshape(data,[data.shape[0],1])\n",
1430     "if scale:\n",
1431     "  scalerx = MinMaxScaler()\n",
1432     "  scalerx.fit(Et)\n",
1433     "  Et = scalerx.transform(Et)\n",
1434     "  scalery = MinMaxScaler()\n",
1435     "  scalery.fit(datat)\n",
1436     "  datat = scalery.transform(datat)"
1437    ]
1438   },
1439   {
1440    "cell_type": "code",
1441    "execution_count": null,
1442    "metadata": {
1443     "id": "PaHfJW7mSJE1"
1444    },
1445    "outputs": [],
1446    "source": [
1447     "def create_RNN_2(hidden_units, dense_units, activation, stateful=False, \n",
1448     "                 batch_shape=None, input_shape=None, dense_layers=1,\n",
1449     "                 rnn_layers=1,return_sequences=False,\n",
1450     "                 initial_state=None):\n",
1451     "    if stateful:\n",
1452     "      inputs = tf.keras.Input(batch_shape=batch_shape)\n",
1453     "    else:\n",
1454     "      inputs = tf.keras.Input(shape=input_shape)\n",
1455     "    # https://stackoverflow.com/questions/43448029/how-can-i-print-the-values-of-keras-tensors\n",
1456     "    # inputs2 = K.print_tensor(inputs, message='inputs = ')  # change allso inputs to inputs2 below, must be used\n",
1457     "    x = inputs\n",
1458     "    for i in range(rnn_layers):\n",
1459     "      x = tf.keras.layers.SimpleRNN(hidden_units,activation=activation[0],\n",
1460     "              stateful=stateful,return_sequences=return_sequences)(x\n",
1461     "              # ,initial_state=initial_state\n",
1462     "              )\n",
1463     "    # x = tf.keras.layers.Dense(hidden_units, activation=activation[1])(x)\n",
1464     "    for i in range(dense_layers):\n",
1465     "      x = tf.keras.layers.Dense(dense_units, activation=activation[1])(x)\n",
1466     "    model = tf.keras.Model(inputs=inputs, outputs=x)\n",
1467     "    model.compile(loss='mean_squared_error', optimizer='adam')\n",
1468     "    return model"
1469    ]
1470   },
1471   {
1472    "cell_type": "code",
1473    "execution_count": null,
1474    "metadata": {
1475     "id": "ZjJxHiEVL5sJ"
1476    },
1477    "outputs": [],
1478    "source": [
1479     "# split data\n",
1480     "return_sequences=False\n",
1481     "shift = 0.\n",
1482     "print('shifting inputs by',shift)\n",
1483     "x_train, y_train = staircase(Et+shift,datat+shift,timesteps=5,trainsteps=h2,\n",
1484     "                             return_sequences=return_sequences)\n",
1485     "print('x_train shape=',x_train.shape)\n",
1486     "samples, timesteps, features = x_train.shape\n",
1487     "print('y_train shape=',y_train.shape)\n",
1488     "# the simplest model possible\n",
1489     "activation=['linear','linear']\n",
1490     "hidden_units=1\n",
1491     "dense_units=1\n",
1492     "dense_layers=1\n",
1493     "features=1\n",
1494     "hours=Et.shape[0]\n",
1495     "h0 = tf.convert_to_tensor(datat[:samples],dtype=tf.float32)\n",
1496     "# print('initial state=',h0)\n",
1497     "# statefull model version for traning\n",
1498     "\n",
1499     "# model_fit=create_RNN_2(hidden_units=hidden_units, \n",
1500     "model_fit=create_RNN_2(hidden_units=hidden_units, \n",
1501     "                        dense_units=dense_units, \n",
1502     "                        batch_shape=(samples,timesteps,features),\n",
1503     "                        stateful=True,\n",
1504     "                        return_sequences=return_sequences,\n",
1505     "                        # initial_state=h0,\n",
1506     "                        activation=activation,\n",
1507     "                        dense_layers=dense_layers)\n",
1508     "# same model stateless for prediction on the entire dataset - to start onlg\n",
1509     "# the real application will switch to prediction after training data end\n",
1510     "# and start from the state there\n",
1511     "print('model_fit input shape',x_train.shape,'output shape',model_fit(x_train).shape)\n",
1512     "from keras.utils.vis_utils import plot_model\n",
1513     "plot_model(model_fit, to_file='model_plot.png', \n",
1514     "           show_shapes=True, show_layer_names=True)"
1515    ]
1516   },
1517   {
1518    "cell_type": "code",
1519    "execution_count": null,
1520    "metadata": {
1521     "id": "ipX9EJqz17Lr"
1522    },
1523    "outputs": [],
1524    "source": [
1525     "model_predict=create_RNN_2(hidden_units=hidden_units, dense_units=dense_units,  \n",
1526     "                        input_shape=(hours,features),stateful = False,\n",
1527     "                        return_sequences=True,\n",
1528     "                        activation=activation,dense_layers=dense_layers)\n",
1529     "# model_predict=create_RNN_sequences(hidden_units=1, dense_units=1, input_shape=(hours,1), \n",
1530     "#                        activation=['linear', 'linear'])\n",
1531     "print('model_predict input shape',Et.shape,'output shape',model_predict(Et).shape)\n",
1532     "print(model_predict.summary())\n",
1533     "from keras.utils.vis_utils import plot_model\n",
1534     "plot_model(model_predict, to_file='model_plot.png', \n",
1535     "           show_shapes=True, show_layer_names=True)"
1536    ]
1537   },
1538   {
1539    "cell_type": "code",
1540    "execution_count": null,
1541    "metadata": {
1542     "id": "dE0OHg0tGVE9"
1543    },
1544    "outputs": [],
1545    "source": [
1546     "# fitting\n",
1547     "w_exact=  [np.array([[1.-np.exp(-0.1)]]), np.array([[np.exp(-0.1)]]), np.array([0.]),np.array([[1.0]]),np.array([-1.*DeltaE])]\n",
1548     "w_initial=[np.array([[1.-np.exp(-0.1)]]), np.array([[np.exp(-0.1)]]), np.array([0.]),np.array([[1.0]]),np.array([0.*DeltaE])]\n",
1549     "model_fit.set_weights(w_initial)\n",
1550     "model_fit.fit(x_train, y_train, epochs=1000, verbose=0,batch_size=samples)\n",
1551     "w_fitted=model_fit.get_weights()\n",
1552     "for i in range(len(w)):\n",
1553     "  print('weight',i,' exact:',w_exact[i],':  initial:',w_initial[i],' fitted:',w_fitted[i])"
1554    ]
1555   },
1556   {
1557    "cell_type": "code",
1558    "execution_count": null,
1559    "metadata": {
1560     "id": "-T8lCS6nBHCj"
1561    },
1562    "outputs": [],
1563    "source": [
1564     "def model_eval(w,title):\n",
1565     "  # prediction on the entire dataset from zero state\n",
1566     "  model_predict.set_weights(w)\n",
1567     "  hours=Et.shape[0]\n",
1568     "  print('Et.shape=',Et.shape,'hours=',hours)\n",
1569     "  x_input=np.reshape(Et,(1, hours, 1))\n",
1570     "  y_output = model_predict.predict(x_input)\n",
1571     "  print('x_input.shape=',x_input.shape,'y_output.shape=',y_output.shape)\n",
1572     "  m = np.reshape(y_output,hours) - shift\n",
1573     "  print('weights=',w)\n",
1574     "  if scale:\n",
1575     "    print('scaling')\n",
1576     "    m = scalery.inverse_transform(m)\n",
1577     "  m = np.reshape(m,hours)\n",
1578     "  plot_m(m,title=title)\n",
1579     "  return m"
1580    ]
1581   },
1582   {
1583    "cell_type": "code",
1584    "execution_count": null,
1585    "metadata": {
1586     "id": "R2jkoZlAIaSb"
1587    },
1588    "outputs": [],
1589    "source": [
1590     "m_fitted=model_eval(w_fitted,'RNN prediction with fitted weights')"
1591    ]
1592   },
1593   {
1594    "cell_type": "code",
1595    "execution_count": null,
1596    "metadata": {
1597     "id": "bO1ewTj9gGvg"
1598    },
1599    "outputs": [],
1600    "source": [
1601     "m_exact=model_eval(w_exact,'RNN prediction with exact weights')"
1602    ]
1603   },
1604   {
1605    "cell_type": "code",
1606    "execution_count": null,
1607    "metadata": {
1608     "id": "henChC0cmbuy"
1609    },
1610    "outputs": [],
1611    "source": [
1612     "m_initial=model_eval(w_initial,'RNN prediction with initial weights')"
1613    ]
1614   },
1615   {
1616    "cell_type": "code",
1617    "execution_count": null,
1618    "metadata": {
1619     "id": "PZw7DNQD4Inr"
1620    },
1621    "outputs": [],
1622    "source": [
1623     "out = np.empty((hours,1))\n",
1624     "w=w_exact\n",
1625     "h=0\n",
1626     "for i in range(Et.shape[0]):\n",
1627     "  h=np.dot(Et[i,0],w[0])+np.dot(h,w[1]) + w[2]\n",
1628     "  out[i]=np.dot(h,w[3]) + w[4]\n",
1629     "if scale:\n",
1630     "  print('scaling')\n",
1631     "  out = scalery.inverse_transform(out)\n",
1632     "out=np.reshape(out,hours)\n",
1633     "print('max abs diff',np.max(np.abs(m_exact-out)))\n",
1634     "plot_m(out,title='Hand computed RNN prediction with exact weights')"
1635    ]
1636   },
1637   {
1638    "cell_type": "markdown",
1639    "metadata": {
1640     "id": "6uXVJj9koGF2"
1641    },
1642    "source": [
1643     "### 3.2 Acquisition and preprocessing of real data"
1644    ]
1645   },
1646   {
1647    "cell_type": "markdown",
1648    "metadata": {
1649     "id": "q3BpOBuzc91i"
1650    },
1651    "source": [
1652     "Data assimilation for fuel moisture from Remote Automated Weather Stations (RAWS) was developed in Vejmelka et al. (2016). First, they use regression from all RAWS in a given area to extend the data spatially from RAWS to a grid in the whole area, then they run the extended Kalman filter at each grid node. Here, we are interested in a simplified problem: estimate future fuel moisture at a single RAWS location from weather data.  "
1653    ]
1654   },
1655   {
1656    "cell_type": "markdown",
1657    "metadata": {
1658     "id": "c8Y6bL1Yc91i"
1659    },
1660    "source": [
1661     "#### 3.2.1 Acquisition of fuel moisture observations"
1662    ]
1663   },
1664   {
1665    "cell_type": "markdown",
1666    "metadata": {
1667     "id": "0CuXyWBFc91i"
1668    },
1669    "source": [
1670     "We try to load the data from a saved file first. If that fails, retrieve the fuel moisture data from sensors on weather stations in the Mesowest network. Get all stations with fuel moisture data in a spatial box within one hour, then pick one station and retrieve the whole time series."
1671    ]
1672   },
1673   {
1674    "cell_type": "code",
1675    "execution_count": null,
1676    "metadata": {},
1677    "outputs": [],
1678    "source": [
1679     "import json\n",
1680     "os.chdir('data')"
1681    ]
1682   },
1683   {
1684    "cell_type": "code",
1685    "execution_count": null,
1686    "metadata": {
1687     "id": "LFrlbbMmc91i"
1688    },
1689    "outputs": [],
1690    "source": [
1691     "jfile = 'raws.json'; vars='fuel_moisture'; case = 1\n",
1692     "# jfile = 'raws2.json'; vars='fuel_moisture,precip_accum_one_hour'; case = 2\n",
1693     "def json_w(j,f):\n",
1694     "  print('writing json file',f)\n",
1695     "  json.dump(j,open(f,'w'),indent=4)\n",
1696     "try:\n",
1697     "    #! wget --no-clobber http://math.ucdenver.edu/~jmandel/data/math4779f21/raws.json\n",
1698     "    j = json.load(open(jfile,'r'))\n",
1699     "    print('loaded from ',jfile)\n",
1700     "    # Take the first station in the boulding box that has data between time_start and time_s2.\n",
1701     "    # Then retrieve data for that station between time_start and time_end\n",
1702     "    time_start = j['time_start']      # start of data time series\n",
1703     "    # time_s2    = j['time_s2']         # end of segment to read coordinates\n",
1704     "    time_end  = j['time_end']         # end of data time series\n",
1705     "    meso_ts  = j['meso_ts']           # get meso observations time series\n",
1706     "    obs_lon =   j['obs_lon']          # where we retrieved observations\n",
1707     "    obs_lat =   j['obs_lat']\n",
1708     "except:\n",
1709     "    print(\"can't read\",jfile,', creating')\n",
1710     "    # set up bounds\n",
1711     "    time_start = \"201806010800\"  # June 1 2018 08:00 in format yyyymmddHHMM\n",
1712     "    time_s2    = \"201806010900\"  # June 1 2018 09:00 in format yyyymmddHHMM \n",
1713     "    time_end   = \"201907200900\"  # June 20 2018 09:00 in format yyyymmddHHMM \n",
1714     "    #time_start=  \"201810230100\"\n",
1715     "    #time_s2=  \"201810230300\"\n",
1716     "    #time_end  =  \"201806022300\"\n",
1717     "    !pip install MesoPy\n",
1718     "    from MesoPy import Meso\n",
1719     "    bounding_box = \"-115, 38, -110, 40\"  # min longtitude, latitude\n",
1720     "    meso_token=\"b40cb52cbdef43ef81329b84e8fd874f\"       # you should get your own if you do more of this\n",
1721     "    m = Meso(meso_token)# create a Meso object\n",
1722     "    print('reading MesoWest fuel moisture data')\n",
1723     "    json_w(m.variables(),'variables.json')\n",
1724     "    meso_obss = m.timeseries(time_start, time_s2, bbox=bounding_box, \n",
1725     "                             showemptystations = '0', vars=vars)   # ask the object for data\n",
1726     "    json_w(meso_obss,'meso_obss.json')                        \n",
1727     "    # pick one station and retrieve the whole time series.\n",
1728     "    station=meso_obss['STATION'][0]\n",
1729     "    json_w(station,'station.json')\n",
1730     "    lon,lat = (float(station['LONGITUDE']),float(station['LATITUDE']))\n",
1731     "    print(station['NAME'],'station',station['STID'],'at',lon,lat)\n",
1732     "    e = 0.01   # tolerance\n",
1733     "    bb = '%s, %s, %s, %s' % (lon - e, lat - e, lon + e, lat + e)\n",
1734     "    print('bounding box',bb)\n",
1735     "    meso_ts = m.timeseries(time_start, time_end, bbox=bb, showemptystations = '0', vars=vars)   # ask the object for data\n",
1736     "    json_w(meso_ts,'meso_ts.json')                        \n",
1737     "    obs_lon, obs_lat = (lon, lat)   # remember station coordinates for later\n",
1738     "    j={'time_start':time_start,'time_s2':time_s2,'time_end':time_end,\n",
1739     "       'meso_ts':meso_ts,'obs_lon':obs_lon,'obs_lat':obs_lat}\n",
1740     "    json_w(j,jfile)\n",
1741     "    print('done')"
1742    ]
1743   },
1744   {
1745    "cell_type": "code",
1746    "execution_count": null,
1747    "metadata": {
1748     "id": "3bXopS3btyz0",
1749     "scrolled": true
1750    },
1751    "outputs": [],
1752    "source": [
1753     "os.chdir('..')\n",
1754     "# process the data retrieved for this station\n",
1755     "# print(json.dumps(meso_ts['STATION'][0], indent=4))\n",
1756     "from datetime import datetime, timedelta, time\n",
1757     "import numpy as np\n",
1758     "import matplotlib.pyplot as plt\n",
1759     "import pytz\n",
1760     "station = meso_ts['STATION'][0]\n",
1761     "time_str  = station['OBSERVATIONS']['date_time']\n",
1762     "obs_time = [datetime.strptime(t, '%Y-%m-%dT%H:%M:%SZ').replace(tzinfo=pytz.UTC) for t in time_str]\n",
1763     "start_time = obs_time[0].replace(minute=0)     # remember obs_time and start_time for later\n",
1764     "end_time = obs_time[-1]\n",
1765     "obs_data = np.array(station['OBSERVATIONS'][\"fuel_moisture_set_1\"])\n",
1766     "# obs_data = np.array(station['OBSERVATIONS'][\"fuel_moisture\"])\n",
1767     "# display the data retrieved\n",
1768     "#for o_time,o_data in zip (obs_time,obs_data):\n",
1769     "#    print(o_time,o_data)\n",
1770     "%matplotlib inline\n",
1771     "plt.figure(figsize=(16,4))\n",
1772     "plt.plot(obs_data,linestyle='-',c='k',label='10-h fuel data')\n",
1773     "plt.title(station['STID'] + ' 10 h fuel moisture data')\n",
1774     "plt.xlabel('Time (hours)') \n",
1775     "plt.ylabel('Fuel moisture content (%)')\n",
1776     "plt.legend()\n",
1777     " "
1778    ]
1779   },
1780   {
1781    "cell_type": "markdown",
1782    "metadata": {
1783     "id": "pY4hPeATK9wZ"
1784    },
1785    "source": [
1786     "#### 3.2.2 Acquisition of weather data"
1787    ]
1788   },
1789   {
1790    "cell_type": "markdown",
1791    "metadata": {
1792     "id": "xhyjXqxVN6B2"
1793    },
1794    "source": [
1795     "Our weather data are results from atmospheric models, with assimilated observations from weather stations, satellites, radars, etc. The models can be run in reanalysis mode (for the past, with data for the period modeled)  or in forecast mode (for the future, with only past data assimilated - because future data are not here yet). We use the Real-Time Mesoscale Analysis ([RTMA](https://www.nco.ncep.noaa.gov/pmb/products/rtma/)) interpolated to the RAWS location. RTMA is a real-time product, posted hourly, and available only for few days in the past. We have our own collection of selected RAWS data over past few years, obtained as a side effect of running the fuel moisture modeling software [WRFXPY](https://github.com/openwfm/wrfxpy).\n",
1796     "\n",
1797     "First try to read the data already extracted for this RAWS and staged for download."
1798    ]
1799   },
1800   {
1801    "cell_type": "code",
1802    "execution_count": null,
1803    "metadata": {
1804     "id": "WlqJRP8Vc91o"
1805    },
1806    "outputs": [],
1807    "source": [
1808     "os.chdir('data')\n",
1809     "import json\n",
1810     "jfile = 'rtma.json'\n",
1811     "try:\n",
1812     "    ! wget --no-clobber http://math.ucdenver.edu/~jmandel/data/math4779f21/rtma.json\n",
1813     "    j = json.load(open(jfile,'r'))\n",
1814     "    print('loaded from ',jfile)\n",
1815     "    if j['obs_lat']!=obs_lat or j['obs_lon']!=obs_lon:\n",
1816     "      print('lon lat doesnot agree, need to load original RTMA files')\n",
1817     "      read_rtma=True\n",
1818     "    else:\n",
1819     "      read_rtma=False\n",
1820     "except:\n",
1821     "    print(\"can't read\",jfile,', creating')\n",
1822     "    read_rtma=True\n",
1823     "\n",
1824     "print('')"
1825    ]
1826   },
1827   {
1828    "cell_type": "markdown",
1829    "metadata": {
1830     "id": "THI6gElyHOOc"
1831    },
1832    "source": [
1833     "Next, functions to get the files, open as grib, and interpolate to the station coordinates"
1834    ]
1835   },
1836   {
1837    "cell_type": "markdown",
1838    "metadata": {
1839     "id": "2iBNHQg5hPxB"
1840    },
1841    "source": [
1842     "####<font color=red>Note: If read_rtma==True, the notebook will say it crashed when run the first time. This is because it needs to install different version of some python packages and restart runtime. Simply run it again.</fonr>"
1843    ]
1844   },
1845   {
1846    "cell_type": "code",
1847    "execution_count": null,
1848    "metadata": {
1849     "id": "mxZABVDxt0gd"
1850    },
1851    "outputs": [],
1852    "source": [
1853     "# Set up environment to read RTMA gribs\n",
1854     "# we will need current numpy for pygrib - needed on Colab, tensorflow is using numpy 1.19\\\n",
1855     "if read_rtma:\n",
1856     "  import subprocess,os\n",
1857     "  def load_rtma(path,file,reload=0):\n",
1858     "    url='http://math.ucdenver.edu/~jmandel/rtma/' + path \n",
1859     "    if os.path.exists(file):\n",
1860     "      if reload:\n",
1861     "        print(file + ' already exists, removing')\n",
1862     "        os.remove(file)\n",
1863     "      else:\n",
1864     "        print(file + ' already exists, exiting')\n",
1865     "        # add checking size here\n",
1866     "        return 0\n",
1867     "    try:\n",
1868     "      ret = subprocess.check_output(['wget','--no-clobber','--output-document='+ file, url,],stderr=subprocess.STDOUT).decode() # execute command from python strings\n",
1869     "      if os.path.exists(file):\n",
1870     "        print('loaded ' + url + ' as ' + file)\n",
1871     "        return 0\n",
1872     "      else: \n",
1873     "        print('file transfer completed, but the file is missing? ' + url)  \n",
1874     "      return 1\n",
1875     "    except:\n",
1876     "      print('file transfer failed: ' + url)\n",
1877     "      return 2\n"
1878    ]
1879   },
1880   {
1881    "cell_type": "markdown",
1882    "metadata": {
1883     "id": "dQ-uJI2sy6I3"
1884    },
1885    "source": [
1886     "Create a function to transfer RTMA files in GRIB2 format from the stash. The function returns zero if the file transfer succeeded. If the file is not available, it returns a nonzero value. Note: if needed, maybe in future add more sophisticated checks, check the return code of wget and if the file size is correct."
1887    ]
1888   },
1889   {
1890    "cell_type": "code",
1891    "execution_count": null,
1892    "metadata": {
1893     "id": "PL3gxK67AlBI"
1894    },
1895    "outputs": [],
1896    "source": [
1897     "if read_rtma:\n",
1898     "  def rtma_grib(t,var):\n",
1899     "    tpath = '%4i%02i%02i/%02i' % (t.year, t.month, t.day, t.hour)  # remote path on server\n",
1900     "    tstr  = '%4i%02i%02i%02i_' % (t.year, t.month, t.day, t.hour)  # time string for local path\n",
1901     "    gribfile = os.path.join('data',tstr + var + '.grib')\n",
1902     "    remote = tpath + '/' + var + '.grib'\n",
1903     "    if load_rtma(remote,gribfile):\n",
1904     "        print('cannot load remote file',remote,'as',gribfile)\n",
1905     "        return []\n",
1906     "    else:\n",
1907     "        try:\n",
1908     "            gf=GribFile(gribfile)\n",
1909     "            v = np.array(gf[1].values())\n",
1910     "        except:\n",
1911     "            print('cannot read grib file',gribfile)\n",
1912     "            return []\n",
1913     "        print('loaded ',gribfile,' containing array shape ',v.shape)\n",
1914     "        return gf[1]   # grib message\n"
1915    ]
1916   },
1917   {
1918    "cell_type": "code",
1919    "execution_count": null,
1920    "metadata": {
1921     "id": "OY1oTYKlfd17"
1922    },
1923    "outputs": [],
1924    "source": [
1925     "if read_rtma:\n",
1926     "    times = pd.date_range(start=time_start,end=time_end,freq='1H')\n",
1927     "    varnames=['temp','td','precipa']\n",
1928     "    j =    read_interp_rtma(varnames,times,obs_lat,obs_lon)      # temperature\n",
1929     "    for varname in varnames:\n",
1930     "        j[varname]=j[varname].tolist() \n",
1931     "    j['obs_lat']=obs_lat\n",
1932     "    j['obs_lon']=obs_lon\n",
1933     "    json.dump(j,open('rtma.json','w'),indent=4)\n",
1934     "    print('done')"
1935    ]
1936   },
1937   {
1938    "cell_type": "code",
1939    "execution_count": null,
1940    "metadata": {
1941     "id": "ccp10kurAlBI"
1942    },
1943    "outputs": [],
1944    "source": [
1945     "from scipy.interpolate import LinearNDInterpolator, interpn\n",
1946     "from scipy.optimize import root\n",
1947     "def interp_to_lat_lon_slow(lats,lons,v,lat,lon): \n",
1948     "    # on mesh with coordinates lats and lons interpolate v to given lat lon\n",
1949     "    interp=LinearNDInterpolator(list(zip(lats.flatten(),lons.flatten())),v.flatten())\n",
1950     "    return interp(lat,lon)\n",
1951     "def interp_to_lat_lon(lats,lons,v,lat,lon):\n",
1952     "    # on mesh with coordinates lats and lons interpolate v to given lat lon\n",
1953     "    points=(np.array(range(lats.shape[0]),float),np.array(range(lats.shape[1]),float))  # uniform mesh\n",
1954     "    def res(ij):  # interpolation of lons lats on the uniform mesh, to noninteger coordinates   \n",
1955     "       return np.hstack((interpn(points,lats,ij)-lat, interpn(points,lons,ij)-lon))\n",
1956     "    # solve for xi,xj such that lats(xi,xj)=lat lons(xi,xj)=lon, then interpolate to (xi, xj) on uniform grid \n",
1957     "    result = root(res,(0,0)) # solve res(ij) = 0\n",
1958     "    if not result.success:\n",
1959     "        print(result.message)\n",
1960     "        exit(1)\n",
1961     "    return interpn(points,v,result.x) \n"
1962    ]
1963   },
1964   {
1965    "cell_type": "markdown",
1966    "metadata": {
1967     "id": "jvnpq6S5AlBI"
1968    },
1969    "source": [
1970     "The interpolation function needs to  be tested."
1971    ]
1972   },
1973   {
1974    "cell_type": "code",
1975    "execution_count": null,
1976    "metadata": {
1977     "id": "NVMJBYI7AlBI"
1978    },
1979    "outputs": [],
1980    "source": [
1981     "def interp_to_lat_lon_test(lats,lons):\n",
1982     "    print('testing interp_to_lat_lon')\n",
1983     "    vx, vy = np.meshgrid(range(lats.shape[0]),range(lats.shape[1]),indexing='ij')\n",
1984     "    i, j = (1,2)\n",
1985     "    lat,lon = ((lats[i,j]+lats[i+1,j+1])/2,(lons[i,j]+lons[i+1,j+1])/2)\n",
1986     "    vi = interp_to_lat_lon(lats,lons,vx,lat,lon)\n",
1987     "    vj = interp_to_lat_lon(lats,lons,vy,lat,lon)\n",
1988     "    print(vi,vj,'should be about',i+0.5,j+0.5)\n",
1989     "    test_slow = 0\n",
1990     "    if test_slow:\n",
1991     "        print('Testing against the standard slow method scipy.interpolate.LinearNDInterpolator. Please wait...')\n",
1992     "        vi_slow = interp_to_lat_lon_slow(lats,lons,vx,lat,lon)\n",
1993     "        print(vi_slow)\n",
1994     "        vj_slow = interp_to_lat_lon_slow(lats,lons,vy,lat,lon)\n",
1995     "        print(vj_slow)\n",
1996     "        \n",
1997     "#gf = rtma_grib(start_time,'temp')      #  read the first grib file and use it to test interpolation\n",
1998     "#lats, lons = gf.latlons()\n",
1999     "#interp_to_lat_lon_test(lats,lons)\n"
2000    ]
2001   },
2002   {
2003    "cell_type": "code",
2004    "execution_count": null,
2005    "metadata": {
2006     "id": "vt-Mk8fIc91m"
2007    },
2008    "outputs": [],
2009    "source": [
2010     "#%debug\n"
2011    ]
2012   },
2013   {
2014    "cell_type": "markdown",
2015    "metadata": {
2016     "id": "LQbWB_3GAlBI"
2017    },
2018    "source": [
2019     "Now we are ready for a function to read the RTMA files and interpolate to the station coordinates"
2020    ]
2021   },
2022   {
2023    "cell_type": "code",
2024    "execution_count": null,
2025    "metadata": {
2026     "id": "b3JJH3XPAlBI"
2027    },
2028    "outputs": [],
2029    "source": [
2030     "if read_rtma:\n",
2031     "  import pandas as pd, json\n",
2032     "  def read_interp_rtma(varnames,times,lat,lon):\n",
2033     "    # read RTMA from start_time to end_time and interpolate to obs_lat obs_lon\n",
2034     "    ntimes = len(times)\n",
2035     "    time_str = 'time_str'\n",
2036     "    j={time_str:times.strftime('%Y-%m-%d %H:%M').tolist()}\n",
2037     "    for varname in varnames:\n",
2038     "        j[varname]=np.full(ntimes,np.nan)  # initialize array of nans as list\n",
2039     "    n=0\n",
2040     "    for t in times:\n",
2041     "        tim=t.strftime('%Y-%m-%d %H:%M')\n",
2042     "        should_be = j[time_str][n]\n",
2043     "        if tim != should_be:\n",
2044     "            print('n=',n,'time',tim,'expected',should_be)\n",
2045     "            raise 'Invalid time' \n",
2046     "        for varname in varnames:\n",
2047     "            gf = rtma_grib(t,varname)   # read and create grib object, download if needed\n",
2048     "            if gf:\n",
2049     "                lats,lons = gf.latlons()    # coordinates\n",
2050     "                v = gf.values()\n",
2051     "                vi=interp_to_lat_lon(lats,lons,v,lat,lon) # append to array\n",
2052     "                print(varname,'at',t,'interpolated to',lat,lon,' value ',vi)\n",
2053     "                j[varname][n] = vi\n",
2054     "            else:\n",
2055     "                print(varname,'at',t,' could not be loaded')\n",
2056     "        n = n+1\n",
2057     "    return j"
2058    ]
2059   },
2060   {
2061    "cell_type": "code",
2062    "execution_count": null,
2063    "metadata": {
2064     "id": "bMpYIZT6c91o"
2065    },
2066    "outputs": [],
2067    "source": [
2068     "# %debug\n"
2069    ]
2070   },
2071   {
2072    "cell_type": "markdown",
2073    "metadata": {
2074     "id": "KVXBjGA0CiXr"
2075    },
2076    "source": [
2077     "#### 3.2.3 Preprocessing and visualization of the weather data"
2078    ]
2079   },
2080   {
2081    "cell_type": "code",
2082    "execution_count": null,
2083    "metadata": {
2084     "id": "fNA3Vbo1c91o"
2085    },
2086    "outputs": [],
2087    "source": [
2088     "rtma = j\n",
2089     "td = np.array(rtma['td'])\n",
2090     "t2 = np.array(rtma['temp'])\n",
2091     "rain=np.array(rtma['precipa'])\n",
2092     "# compute relative humidity\n",
2093     "rh = 100*np.exp(17.625*243.04*(td - t2) / (243.04 + t2 - 273.15) / (243.0 + td - 273.15))\n",
2094     "Ed = 0.924*rh**0.679 + 0.000499*np.exp(0.1*rh) + 0.18*(21.1 + 273.15 - t2)*(1 - np.exp(-0.115*rh))\n",
2095     "Ew = 0.618*rh**0.753 + 0.000454*np.exp(0.1*rh) + 0.18*(21.1 + 273.15 - t2)*(1 - np.exp(-0.115*rh))"
2096    ]
2097   },
2098   {
2099    "cell_type": "code",
2100    "execution_count": null,
2101    "metadata": {
2102     "id": "tZIK59bJAlBJ"
2103    },
2104    "outputs": [],
2105    "source": [
2106     "%matplotlib inline\n",
2107     "plt.figure(figsize=(16,4))\n",
2108     "plt.plot(t2,linestyle='-',c='k',label='Temperature')\n",
2109     "plt.title(station['STID'] + ' Temperature')\n",
2110     "plt.xlabel('Time (hours)') \n",
2111     "plt.ylabel('Temperature (K)')\n",
2112     "plt.legend()"
2113    ]
2114   },
2115   {
2116    "cell_type": "code",
2117    "execution_count": null,
2118    "metadata": {
2119     "id": "LbyqcuXYAlBJ"
2120    },
2121    "outputs": [],
2122    "source": [
2123     "%matplotlib inline\n",
2124     "plt.figure(figsize=(16,4))\n",
2125     "plt.plot(td,linestyle='-',c='k',label='Dew point')\n",
2126     "plt.title(station['STID'] + ' Dew point (K)')\n",
2127     "plt.xlabel('Time (hours)') \n",
2128     "plt.ylabel('Dew point (K)')\n",
2129     "plt.legend()"
2130    ]
2131   },
2132   {
2133    "cell_type": "code",
2134    "execution_count": null,
2135    "metadata": {
2136     "id": "dfoOK2kSc91p"
2137    },
2138    "outputs": [],
2139    "source": [
2140     "%matplotlib inline\n",
2141     "plt.figure(figsize=(16,4))\n",
2142     "plt.plot(rh,linestyle='-',c='k',label='Dew point')\n",
2143     "plt.title(station['STID'] + ' relative humidity')\n",
2144     "plt.xlabel('Time (hours)') \n",
2145     "plt.ylabel('Relative humidity (%)')\n",
2146     "plt.legend()"
2147    ]
2148   },
2149   {
2150    "cell_type": "code",
2151    "execution_count": null,
2152    "metadata": {
2153     "id": "MWTJ5b2kc91p"
2154    },
2155    "outputs": [],
2156    "source": [
2157     "%matplotlib inline\n",
2158     "plt.figure(figsize=(16,4))\n",
2159     "plt.plot(Ed,linestyle='-',c='r',label='drying equilibrium')\n",
2160     "plt.plot(Ew,linestyle=':',c='b',label='wetting equilibrium')\n",
2161     "plt.title(station['STID'] + ' drying and wetting equilibria')\n",
2162     "plt.xlabel('Time (hours)') \n",
2163     "plt.ylabel('Fuel moisture contents (%)')\n",
2164     "plt.legend()"
2165    ]
2166   },
2167   {
2168    "cell_type": "markdown",
2169    "metadata": {
2170     "id": "jY3_eeBRc91p"
2171    },
2172    "source": [
2173     " "
2174    ]
2175   },
2176   {
2177    "cell_type": "code",
2178    "execution_count": null,
2179    "metadata": {
2180     "id": "PQKSRvRSAlBJ"
2181    },
2182    "outputs": [],
2183    "source": [
2184     "%matplotlib inline\n",
2185     "plt.figure(figsize=(16,4))\n",
2186     "plt.plot(rain,linestyle='-',c='k',label='Precipitation')\n",
2187     "plt.title(station['STID'] + ' Precipitation' )\n",
2188     "plt.xlabel('Time (hours)') \n",
2189     "plt.ylabel('Precipitation (mm/hour)')\n",
2190     "plt.legend()"
2191    ]
2192   },
2193   {
2194    "cell_type": "code",
2195    "execution_count": null,
2196    "metadata": {
2197     "id": "Dwbt4UXfro5x"
2198    },
2199    "outputs": [],
2200    "source": [
2201     "print(rain[1900:2000])"
2202    ]
2203   },
2204   {
2205    "cell_type": "markdown",
2206    "metadata": {
2207     "id": "_yRu_7WvHc6P"
2208    },
2209    "source": [
2210     "Precipitation from RTMA is in kg/m${}^2$. 1m water depth over 1m${}^2$ is 1m${}^3$ with mass 1000 kg thus 1 kg/m${}^2$ is the same as 1 mm of precipitation. RTMA values are accumulations over 1 h so these are values in mm/h. So 9999 mm/h = 10m/h makes no sense. Replace anything over 1m/h by nan and try again."
2211    ]
2212   },
2213   {
2214    "cell_type": "code",
2215    "execution_count": null,
2216    "metadata": {
2217     "id": "XPYO_Iuvc91q"
2218    },
2219    "outputs": [],
2220    "source": [
2221     "rain[rain > 1000] = np.NaN"
2222    ]
2223   },
2224   {
2225    "cell_type": "code",
2226    "execution_count": null,
2227    "metadata": {
2228     "id": "GYWTfbBBc91q",
2229     "scrolled": true
2230    },
2231    "outputs": [],
2232    "source": [
2233     "%matplotlib inline\n",
2234     "plt.figure(figsize=(16,4))\n",
2235     "plt.plot(rain,linestyle='-',c='k',label='Precipitation')\n",
2236     "plt.title(station['STID'] + ' Precipitation' )\n",
2237     "plt.xlabel('Time (hours)') \n",
2238     "plt.ylabel('Precipitation (mm/hour)')\n",
2239     "plt.legend()"
2240    ]
2241   },
2242   {
2243    "cell_type": "markdown",
2244    "metadata": {
2245     "id": "Q_L0R2Njc91q"
2246    },
2247    "source": [
2248     "Fix some missing data, then we can use the data for up to 1942 hours until a biger gap."
2249    ]
2250   },
2251   {
2252    "cell_type": "code",
2253    "execution_count": null,
2254    "metadata": {
2255     "id": "_tkU7UJic91q"
2256    },
2257    "outputs": [],
2258    "source": [
2259     "# fix isolated nans\n",
2260     "def fixnan(a,n):\n",
2261     "    for c in range(n):\n",
2262     "        for i in np.where(np.isnan(a)):\n",
2263     "            a[i]=0.5*(a[i-1]+a[i+1])\n",
2264     "        if not any(np.isnan(a)):\n",
2265     "            break\n",
2266     "    return a\n",
2267     "\n",
2268     "rain=fixnan(rain,2)\n",
2269     "t2=fixnan(t2,2)\n",
2270     "rh=fixnan(rh,2)\n",
2271     "obs_data=fixnan(obs_data,2)\n",
2272     "Ed=fixnan(Ed,2)\n",
2273     "Ew=fixnan(Ew,2)\n",
2274     "\n",
2275     "print(np.where(np.isnan(rain)))\n",
2276     "print(np.where(np.isnan(t2)))\n",
2277     "print(np.where(np.isnan(rh)))\n",
2278     "print(np.where(np.isnan(obs_data)))"
2279    ]
2280   },
2281   {
2282    "cell_type": "markdown",
2283    "metadata": {
2284     "id": "XqQYnyI9DIy1"
2285    },
2286    "source": [
2287     "## 4 Results"
2288    ]
2289   },
2290   {
2291    "cell_type": "markdown",
2292    "metadata": {
2293     "id": "2tIC_Tqnc91r"
2294    },
2295    "source": [
2296     "### 4.1 Kalman filter with fuel moisture observations, followed by forecasting\n",
2297     "We run the model first with Kalman filter for 150 hours. The observations are the RAWS data\n",
2298     "After 150 hours, we run in forecast mode - the RAWS data are no longer used, and we run the model from the weather data without the Kalman filter. The weather data are taken to be RTMA interpolated to one RAWS location.\n",
2299     "In a real forecasting application, the model would be run from weather forecast rather than data."
2300    ]
2301   },
2302   {
2303    "cell_type": "code",
2304    "execution_count": null,
2305    "metadata": {
2306     "id": "aXnSQM7wc91r"
2307    },
2308    "outputs": [],
2309    "source": [
2310     "# run KF on an initial data seqment\n",
2311     "import numpy as np\n",
2312     "import matplotlib.pyplot as plt \n",
2313     "\n",
2314     "hours=1200 # total simulation\n",
2315     "h2 = 300\n",
2316     "m = np.zeros(hours) # preallocate\n",
2317     "m[0]= obs_data[0]             # initial state  \n",
2318     "P = np.zeros(hours)\n",
2319     "P[0] = 1e-3 # background state variance\n",
2320     "H = np.array([1.])   # all oQ = np.array([0.02]) # process noise variancebserved\n",
2321     "Q = np.array([1e-3]) # process noise variance\n",
2322     "R = np.array([1e-3]) # data variance\n",
2323     "for t in range(hours-1):\n",
2324     "    # using lambda construction to pass additional arguments to the model \n",
2325     "    if t < h2 and not np.isnan(obs_data[t]) and not np.isnan(Ew[t]) and not np.isnan(rain[t]): # advance model and run KF\n",
2326     "        m[t+1],P[t+1] = ext_kf(m[t],P[t],lambda u: model_moisture(u,Ed[t],Ew[t],rain[t],t,partials=1),Q,\n",
2327     "                    d=obs_data[t],H=H,R=R)\n",
2328     "    else:  # just advance to next hour, no process noise\n",
2329     "        m[t+1],P[t+1] = ext_kf(m[t],P[t],lambda u: model_moisture(u,Ed[t],Ew[t],rain[t],t,partials=1),Q*0.0)"
2330    ]
2331   },
2332   {
2333    "cell_type": "code",
2334    "execution_count": null,
2335    "metadata": {
2336     "id": "peMi-OF3c91r",
2337     "scrolled": true
2338    },
2339    "outputs": [],
2340    "source": [
2341     "%matplotlib inline\n",
2342     "plt.figure(figsize=(16,4))\n",
2343     "plt.plot(Ed[:hours],linestyle='--',c='r',label='Drying Equilibrium')\n",
2344     "plt.plot(Ew[:hours],linestyle='--',c='b',label='Wetting Equilibrium')\n",
2345     "plt.plot(obs_data[:hours],linestyle=':',c='k',label='RAWS data')\n",
2346     "plt.plot(m[:h2],linestyle='-',c='k',label='filtered')\n",
2347     "plt.plot(range(h2,hours),m[h2:hours],linestyle='-',c='r',label='forecast')\n",
2348     "plt.title(station['STID'] + ' Kalman filtering and forecast with real data')\n",
2349     "plt.xlabel('Time (hours)') \n",
2350     "plt.ylabel('Fuel moisture content (%)')\n",
2351     "plt.legend()"
2352    ]
2353   },
2354   {
2355    "cell_type": "markdown",
2356    "metadata": {
2357     "id": "3TnwXYcLc91r"
2358    },
2359    "source": [
2360     "Clearly, there is a problem - the forecast fuel moisture is too high. We need to assimilate also some parameters of the model, not just its output state. "
2361    ]
2362   },
2363   {
2364    "cell_type": "markdown",
2365    "metadata": {
2366     "id": "8SuVNg8TsW4d"
2367    },
2368    "source": [
2369     "### 4.3 Kalman filter on the augmented model"
2370    ]
2371   },
2372   {
2373    "cell_type": "markdown",
2374    "metadata": {
2375     "id": "FYAbWNCfk2wD"
2376    },
2377    "source": [
2378     "Run augmented filter and plot the result:\n"
2379    ]
2380   },
2381   {
2382    "cell_type": "code",
2383    "execution_count": null,
2384    "metadata": {
2385     "id": "Q3NHr3oBsDg6"
2386    },
2387    "outputs": [],
2388    "source": [
2389     "m,Ec = run_augmented_kf(obs_data,Ed,Ew,rain,h2,hours)  # extract from state"
2390    ]
2391   },
2392   {
2393    "cell_type": "code",
2394    "execution_count": null,
2395    "metadata": {
2396     "id": "hlkby_oTlB_f"
2397    },
2398    "outputs": [],
2399    "source": [
2400     "title = station['STID'] +' Kalman filtering and forecast with augmented state, real data. Training 0:%i hmax' % h2\n",
2401     "def plot_moisture(hmin,hmax):\n",
2402     "  print('training from 0 to',h2,'plot from',hmin,'to',hmax)\n",
2403     "  plt.figure(figsize=(16,4))\n",
2404     "  plt.plot(range(hmin,hmax),Ed[hmin:hmax],linestyle='--',c='r',label='Drying Equilibrium (%)')\n",
2405     "  plt.plot(range(hmin,hmax),Ew[hmin:hmax],linestyle='--',c='b',label='Wetting Equilibrium (%)')\n",
2406     "  plt.plot(range(hmin,hmax),Ec[hmin:hmax],linestyle='--',c='g',label='Equilibrium Correction (%)')\n",
2407     "  plt.plot(range(hmin,hmax),m[hmin:hmax],linestyle='-',c='k',label='filtered')\n",
2408     "  plt.plot(range(hmin,hmax),obs_data[hmin:hmax],linestyle='-',c='b',label='RAWS data (%)')\n",
2409     "  plt.plot(range(hmin,hmax),rain[hmin:hmax],linestyle='-',c='b',label='RTMA rain (mm/h)')\n",
2410     "  plt.title(title)\n",
2411     "  if hmin>=h2:\n",
2412     "    plt.plot(m[hmin:h2],linestyle='-',c='k',label='Filtered')\n",
2413     "  h1 = np.maximum(hmin,h2)\n",
2414     "  plt.plot(range(h1,hmax),m[h1:hmax],linestyle='-',c='r',label='Forecast (%)')\n",
2415     "  plt.xlabel('Time (hours)') \n",
2416     "  plt.ylabel('Fuel moisture content (%)')\n",
2417     "  plt.legend()"
2418    ]
2419   },
2420   {
2421    "cell_type": "code",
2422    "execution_count": null,
2423    "metadata": {},
2424    "outputs": [],
2425    "source": [
2426     "os.chdir('..')\n",
2427     "from data_funcs import to_json, from_json"
2428    ]
2429   },
2430   {
2431    "cell_type": "code",
2432    "execution_count": null,
2433    "metadata": {},
2434    "outputs": [],
2435    "source": [
2436     "os.chdir('data')\n",
2437     "kf_orig={'title':title,'hours':hours,'h2':h2,'Ed':Ed,'Ew':Ew,'Ec':Ec,'rain':rain,\n",
2438     "            'fm':obs_data,'m':m,'note':'RAWS and RTMA data + m from augmented KF in fmda_kf_rnn_orig'}\n",
2439     "to_json(kf_orig,'kf_orig.json')"
2440    ]
2441   },
2442   {
2443    "cell_type": "code",
2444    "execution_count": null,
2445    "metadata": {
2446     "id": "q-h5omKgnow2"
2447    },
2448    "outputs": [],
2449    "source": [
2450     "plot_moisture(0,hours)"
2451    ]
2452   },
2453   {
2454    "cell_type": "markdown",
2455    "metadata": {
2456     "id": "0w0YtHtqnza5"
2457    },
2458    "source": [
2459     "A detailed view of transition from training to forecast:"
2460    ]
2461   },
2462   {
2463    "cell_type": "code",
2464    "execution_count": null,
2465    "metadata": {
2466     "id": "B7sXGUotc91s"
2467    },
2468    "outputs": [],
2469    "source": [
2470     "plot_moisture(0,600)\n",
2471     "\n"
2472    ]
2473   },
2474   {
2475    "cell_type": "code",
2476    "execution_count": null,
2477    "metadata": {
2478     "id": "xy7sIs0z_Kk6"
2479    },
2480    "outputs": [],
2481    "source": [
2482     "plot_moisture(300,500)"
2483    ]
2484   },
2485   {
2486    "cell_type": "code",
2487    "execution_count": null,
2488    "metadata": {
2489     "id": "y-C6IRFVxGUR"
2490    },
2491    "outputs": [],
2492    "source": [
2493     "plot_moisture(300,800)"
2494    ]
2495   },
2496   {
2497    "cell_type": "code",
2498    "execution_count": null,
2499    "metadata": {
2500     "id": "TvlCtT0X2ejp"
2501    },
2502    "outputs": [],
2503    "source": [
2504     "plot_moisture(800,1200)"
2505    ]
2506   },
2507   {
2508    "cell_type": "markdown",
2509    "metadata": {
2510     "id": "7W03QTo3c91t"
2511    },
2512    "source": [
2513     "Filtering by extended Kalman filter using RAWS data until 150 hours, then forecasting mode - running the model from interpolated RTMA only. For the first 60 hours the forecast is good, the equilibium correction made the model quite close to data. But then the big spike in equilibrium moisture around 230 hours attracted the solution, and it took a while for it to get back. The spike in the RAWS measurement is there but much smaller. The model becomes inaccurate during periods when the fuel moisture equilibrium is large.\n",
2514     "\n",
2515     "Possible reasons include: 1. There was something in the data we do not know about - maybe it rained but RTMA did not tell us. Try comparing with data from the RAWS itself? 2. The model is too simple, assumes the whole depth of the wood stick is wetting and drying at the same time. Perhaps the moisture got stored in the inside layers of the measurement stick. Try a two-layer model as in van der Kamp (2017) and make the state larger? "
2516    ]
2517   },
2518   {
2519    "cell_type": "markdown",
2520    "metadata": {
2521     "id": "owEI4EtTo7Ek"
2522    },
2523    "source": [
2524     "A detailed view of rain episode:"
2525    ]
2526   },
2527   {
2528    "cell_type": "code",
2529    "execution_count": null,
2530    "metadata": {
2531     "id": "C_hoDjgtpMEJ"
2532    },
2533    "outputs": [],
2534    "source": [
2535     "plot_moisture(900,1100)"
2536    ]
2537   },
2538   {
2539    "cell_type": "markdown",
2540    "metadata": {
2541     "id": "DRraWhwdpSkV"
2542    },
2543    "source": [
2544     "It seems there is some rain that the model does not know about."
2545    ]
2546   },
2547   {
2548    "cell_type": "markdown",
2549    "metadata": {
2550     "id": "1STfnlT40rPX"
2551    },
2552    "source": [
2553     "## RNN for real data, no rain yet"
2554    ]
2555   },
2556   {
2557    "cell_type": "markdown",
2558    "metadata": {
2559     "id": "3cwY43iSnQ0t"
2560    },
2561    "source": [
2562     "#### Linear modeling by RELU - potential for generalization"
2563    ]
2564   },
2565   {
2566    "cell_type": "code",
2567    "execution_count": null,
2568    "metadata": {
2569     "id": "MotzNBvOnFvC"
2570    },
2571    "outputs": [],
2572    "source": [
2573     "def RELU(x):\n",
2574     "  if x>0. :\n",
2575     "    return x\n",
2576     "  else:\n",
2577     "    return 0.\n",
2578     "\n",
2579     "# network computing z = a*x1 + b*x2 with offset c\n",
2580     "def linrelu(x,a,b,c):\n",
2581     "  y = np.dot(np.array([[a, b], [-a, -b] ]), x) + np.array([c, -c])\n",
2582     "  y[0]=RELU(y[0])\n",
2583     "  y[1]=RELU(y[1])\n",
2584     "  return(np.dot([1,-1],y))-c\n",
2585     "x = [1,2]\n",
2586     "a = 2\n",
2587     "b = 4\n",
2588     "c = 3\n",
2589     "print(a*x[0]+b*x[1])\n",
2590     "linrelu(x,a,b,c)"
2591    ]
2592   },
2593   {
2594    "cell_type": "markdown",
2595    "metadata": {
2596     "id": "-p6dcLua_udD"
2597    },
2598    "source": [
2599     "### Basic RNN on real data "
2600    ]
2601   },
2602   {
2603    "cell_type": "markdown",
2604    "metadata": {
2605     "id": "gSmbDPZIHbTr"
2606    },
2607    "source": [
2608     "Try with E average between drying and wetting"
2609    ]
2610   },
2611   {
2612    "cell_type": "code",
2613    "execution_count": null,
2614    "metadata": {
2615     "id": "ymhNMZkoHfCl"
2616    },
2617    "outputs": [],
2618    "source": [
2619     "E = (Ed + Ew)/2\n",
2620     "print(Ed.shape,Ew.shape,rain.shape)\n",
2621     "first_rain=np.nonzero(rain>0)[0][0]\n",
2622     "print(first_rain)\n",
2623     "hours=first_rain\n",
2624     "E=E[:hours]\n",
2625     "data=obs_data[:hours]\n",
2626     "scale=False\n",
2627     "\n",
2628     "# transform as 2D, (timesteps, features) and (timesteps, outputs)\n",
2629     "Et = np.reshape(E,[E.shape[0],1])\n",
2630     "datat = np.reshape(data,[data.shape[0],1])\n",
2631     "if scale:\n",
2632     "  scalerx = MinMaxScaler()\n",
2633     "  scalerx.fit(Et)\n",
2634     "  Et = scalerx.transform(Et)\n",
2635     "  scalery = MinMaxScaler()\n",
2636     "  scalery.fit(datat)\n",
2637     "  datat = scalery.transform(datat)"
2638    ]
2639   },
2640   {
2641    "cell_type": "markdown",
2642    "metadata": {
2643     "id": "DPcxv85XILdn"
2644    },
2645    "source": [
2646     "Create the model again"
2647    ]
2648   },
2649   {
2650    "cell_type": "code",
2651    "execution_count": null,
2652    "metadata": {
2653     "id": "gEkbHZSqIOq1"
2654    },
2655    "outputs": [],
2656    "source": [
2657     "# split data\n",
2658     "return_sequences=False\n",
2659     "x_train, y_train = staircase(Et,datat,timesteps=5,trainsteps=h2,\n",
2660     "                             return_sequences=return_sequences)\n",
2661     "print('x_train shape=',x_train.shape)\n",
2662     "samples, timesteps, features = x_train.shape\n",
2663     "print('y_train shape=',y_train.shape)\n",
2664     "# the simplest model possible\n",
2665     "activation=['linear','linear']\n",
2666     "hidden_units=3\n",
2667     "dense_units=1\n",
2668     "dense_layers=1\n",
2669     "features=1\n",
2670     "hours=Et.shape[0]\n",
2671     "h0 = tf.convert_to_tensor(datat[:samples],dtype=tf.float32)\n",
2672     "# print('initial state=',h0)\n",
2673     "# statefull model version for traning\n",
2674     "\n",
2675     "model_fit=create_RNN_2(hidden_units=hidden_units, \n",
2676     "                        dense_units=dense_units, \n",
2677     "                        batch_shape=(samples,timesteps,features),\n",
2678     "                        stateful=True,\n",
2679     "                        return_sequences=return_sequences,\n",
2680     "                        # initial_state=h0,\n",
2681     "                        activation=activation,\n",
2682     "                        dense_layers=dense_layers)\n",
2683     "# same model stateless for prediction on the entire dataset - to start onlg\n",
2684     "# the real application will switch to prediction after training data end\n",
2685     "# and start from the state there\n",
2686     "print('model_fit input shape',x_train.shape,'output shape',model_fit(x_train).shape)\n",
2687     "from keras.utils.vis_utils import plot_model\n",
2688     "plot_model(model_fit, to_file='model_plot.png', \n",
2689     "           show_shapes=True, show_layer_names=True)"
2690    ]
2691   },
2692   {
2693    "cell_type": "code",
2694    "execution_count": null,
2695    "metadata": {
2696     "id": "jtFJQu33NqfL"
2697    },
2698    "outputs": [],
2699    "source": [
2700     "model_predict=create_RNN_2(hidden_units=hidden_units, dense_units=dense_units,  \n",
2701     "                        input_shape=(hours,features),stateful = False,\n",
2702     "                        return_sequences=True,\n",
2703     "                        activation=activation,dense_layers=dense_layers)\n",
2704     "# model_predict=create_RNN_sequences(hidden_units=1, dense_units=1, input_shape=(hours,1), \n",
2705     "#                        activation=['linear', 'linear'])\n",
2706     "print('model_predict input shape',Et.shape,'output shape',model_predict(Et).shape)\n",
2707     "print(model_predict.summary())\n",
2708     "from keras.utils.vis_utils import plot_model\n",
2709     "plot_model(model_predict, to_file='model_plot.png', \n",
2710     "           show_shapes=True, show_layer_names=True)"
2711    ]
2712   },
2713   {
2714    "cell_type": "code",
2715    "execution_count": null,
2716    "metadata": {
2717     "id": "wuxh5pq0OMSa"
2718    },
2719    "outputs": [],
2720    "source": [
2721     "# fitting\n",
2722     "DeltaE = 0\n",
2723     "w_exact=  [np.array([[1.-np.exp(-0.1)]]), np.array([[np.exp(-0.1)]]), np.array([0.]),np.array([[1.0]]),np.array([-1.*DeltaE])]\n",
2724     "w_initial=[np.array([[1.-np.exp(-0.1)]]), np.array([[np.exp(-0.1)]]), np.array([0.]),np.array([[1.0]]),np.array([-1.0])]\n",
2725     "w=model_fit.get_weights()\n",
2726     "for i in range(len(w)):\n",
2727     "  print('weight',i,'shape',w[i].shape,'ndim',w[i].ndim,'given',w_initial[i].shape)\n",
2728     "  for j in range(w[i].shape[0]):\n",
2729     "    if w[i].ndim==2:\n",
2730     "      for k in range(w[i].shape[1]):\n",
2731     "        w[i][j][k]=w_initial[i][0][0]/w[i].shape[0]\n",
2732     "    else:\n",
2733     "      w[i][j]=w_initial[i][0]\n",
2734     "model_fit.set_weights(w)\n",
2735     "model_fit.fit(x_train, y_train, epochs=5000, verbose=0,batch_size=samples)\n",
2736     "w_fitted=model_fit.get_weights()\n",
2737     "for i in range(len(w)):\n",
2738     "  print('weight',i,' exact:',w_exact[i],':  initial:',w_initial[i],' fitted:',w_fitted[i])"
2739    ]
2740   },
2741   {
2742    "cell_type": "code",
2743    "execution_count": null,
2744    "metadata": {
2745     "tags": []
2746    },
2747    "outputs": [],
2748    "source": []
2749   },
2750   {
2751    "cell_type": "code",
2752    "execution_count": null,
2753    "metadata": {
2754     "id": "uJz1EgPyRTEH"
2755    },
2756    "outputs": [],
2757    "source": [
2758     "# evaluate model\n",
2759     "model_predict.set_weights(w_fitted)\n",
2760     "x_input=np.reshape(Et,(1, hours, 1))\n",
2761     "y_output = model_predict.predict(x_input)\n",
2762     "print('x_input.shape=',x_input.shape,'y_output.shape=',y_output.shape)\n",
2763     "print(shift)\n",
2764     "m = np.reshape(y_output,hours)\n",
2765     "print('weights=',w)\n",
2766     "if scale:\n",
2767     "    print('scaling')\n",
2768     "    m = scalery.inverse_transform(m)\n",
2769     "m = np.reshape(m,hours)\n",
2770     "hour=np.array(range(hours))\n",
2771     "title=\"First RNN forecast\"\n",
2772     "plt.figure(figsize=(16,4))\n",
2773     "plt.plot(hour,E,linestyle='--',c='r',label='E=Equilibrium data')\n",
2774     "# print(len(hour),len(m_f))\n",
2775     "plt.scatter(hour,data,c='b',label='data=10-h fuel data')\n",
2776     "if m is not None:\n",
2777     "    plt.plot(hour[:h2],m[:h2],linestyle='-',c='k',label='m=filtered')\n",
2778     "    plt.plot(hour[h2:hours],m[h2:hours],linestyle='-',c='r',label='m=forecast')\n",
2779     "plt.title(title) \n",
2780     "plt.legend()\n"
2781    ]
2782   },
2783   {
2784    "cell_type": "code",
2785    "execution_count": null,
2786    "metadata": {},
2787    "outputs": [],
2788    "source": []
2789   },
2790   {
2791    "cell_type": "code",
2792    "execution_count": null,
2793    "metadata": {
2794     "id": "VSwtgKPJPnH4"
2795    },
2796    "outputs": [],
2797    "source": [
2798     "# plot subinterval only\n",
2799     "def plot_int(lb=0,ub=hours,title=\"RNN forecast\"):\n",
2800     "  hour=np.array(range(hours))\n",
2801     "  plt.figure(figsize=(16,4))\n",
2802     "  plt.plot(hour[lb:ub],E[lb:ub],linestyle='--',c='r',label='Equilibrium data')\n",
2803     "  # plt.scatter(hour[lb:ub],data[lb:ub],c='b',label='data=10-h fuel data')\n",
2804     "  plt.plot(hour[lb:ub],m[lb:ub],linestyle='-',c='b',label='data=10-h fuel data')\n",
2805     "  if lb <= h2:\n",
2806     "    ub1 = min(h2,ub)\n",
2807     "    plt.plot(hour[lb:ub1],m[lb:ub1],linestyle='-',c='k',label='filtered')\n",
2808     "  if ub >= h2:\n",
2809     "    lb1 = max(h2,lb)\n",
2810     "    plt.plot(hour[lb1:ub],m[lb1:ub],linestyle='-',c='r',label='forecast')\n",
2811     "  plt.title(title) \n",
2812     "  plt.legend()"
2813    ]
2814   },
2815   {
2816    "cell_type": "code",
2817    "execution_count": null,
2818    "metadata": {
2819     "id": "vCjk9hZtkFym"
2820    },
2821    "outputs": [],
2822    "source": [
2823     "plot_int()"
2824    ]
2825   },
2826   {
2827    "cell_type": "code",
2828    "execution_count": null,
2829    "metadata": {
2830     "id": "Sd3fDOnvmmdp"
2831    },
2832    "outputs": [],
2833    "source": [
2834     "plot_int(0,300)"
2835    ]
2836   },
2837   {
2838    "cell_type": "code",
2839    "execution_count": null,
2840    "metadata": {
2841     "id": "vHkc4KHdkAJp"
2842    },
2843    "outputs": [],
2844    "source": [
2845     "plot_int(300,500)"
2846    ]
2847   },
2848   {
2849    "cell_type": "code",
2850    "execution_count": null,
2851    "metadata": {
2852     "id": "Km5VWhcJlyvV"
2853    },
2854    "outputs": [],
2855    "source": [
2856     "plot_int(500,800)"
2857    ]
2858   },
2859   {
2860    "cell_type": "markdown",
2861    "metadata": {
2862     "id": "TBayRudFcZWP"
2863    },
2864    "source": [
2865     "Next step: two features - drying and wetting equilibria"
2866    ]
2867   },
2868   {
2869    "cell_type": "code",
2870    "execution_count": null,
2871    "metadata": {
2872     "id": "SGbgxOm_kEc4"
2873    },
2874    "outputs": [],
2875    "source": [
2876     "print(Ed.shape,Ew.shape,rain.shape)\n",
2877     "first_rain=np.nonzero(rain>0)[0][0]\n",
2878     "print(first_rain)\n",
2879     "hours=first_rain\n",
2880     "Ed=Ed[:hours]\n",
2881     "Ew=Ew[:hours]\n",
2882     "h2 = 300\n",
2883     "# print(Ed.shape,Ew.shape)\n",
2884     "# (timesteps, features)\n",
2885     "Et = np.vstack((Ed, Ew)).T\n",
2886     "print(E.shape)\n",
2887     "data=obs_data[:hours]\n",
2888     "\n",
2889     "scale=False\n",
2890     "\n",
2891     "# transform as 2D, (timesteps, features) and (timesteps, outputs)\n",
2892     "datat = np.reshape(data,[data.shape[0],1])\n",
2893     "if scale:\n",
2894     "  scalerx = MinMaxScaler()\n",
2895     "  scalerx.fit(Et)\n",
2896     "  Et = scalerx.transform(Et)\n",
2897     "  scalery = MinMaxScaler()\n",
2898     "  scalery.fit(datat)\n",
2899     "  datat = scalery.transform(datat)"
2900    ]
2901   },
2902   {
2903    "cell_type": "code",
2904    "execution_count": null,
2905    "metadata": {},
2906    "outputs": [],
2907    "source": [
2908     "from utils import hash2"
2909    ]
2910   },
2911   {
2912    "cell_type": "code",
2913    "execution_count": null,
2914    "metadata": {},
2915    "outputs": [],
2916    "source": [
2917     "# Set seed for reproducibility\n",
2918     "reproducibility.set_seed()"
2919    ]
2920   },
2921   {
2922    "cell_type": "code",
2923    "execution_count": null,
2924    "metadata": {
2925     "id": "b6aJAvBEkEBl"
2926    },
2927    "outputs": [],
2928    "source": [
2929     "# split data\n",
2930     "return_sequences=False\n",
2931     "x_train, y_train = staircase(Et,datat,timesteps=5,trainsteps=h2,\n",
2932     "                             return_sequences=return_sequences)\n",
2933     "print('x_train shape=',x_train.shape)\n",
2934     "samples, timesteps, features = x_train.shape\n",
2935     "print('y_train shape=',y_train.shape)\n",
2936     "# the simplest model possible\n",
2937     "activation=['linear','linear']\n",
2938     "hidden_units=6\n",
2939     "dense_units=1\n",
2940     "dense_layers=1\n",
2941     "features=Et.shape[1]\n",
2942     "hours=Et.shape[0]\n",
2943     "h0 = tf.convert_to_tensor(datat[:samples],dtype=tf.float32)\n",
2944     "# print('initial state=',h0)\n",
2945     "# statefull model version for traning\n",
2946     "\n",
2947     "model_fit=create_RNN_2(hidden_units=hidden_units, \n",
2948     "                        dense_units=dense_units, \n",
2949     "                        batch_shape=(samples,timesteps,features),\n",
2950     "                        stateful=True,\n",
2951     "                        return_sequences=return_sequences,\n",
2952     "                        # initial_state=h0,\n",
2953     "                        activation=activation,\n",
2954     "                        dense_layers=dense_layers)\n",
2955     "# same model stateless for prediction on the entire dataset - to start onlg\n",
2956     "# the real application will switch to prediction after training data end\n",
2957     "# and start from the state there\n",
2958     "print('model_fit input shape',x_train.shape,'output shape',model_fit(x_train).shape)\n",
2959     "print('model_fit input shape',x_train.shape,'output shape',y_train.shape)\n",
2960     "from keras.utils.vis_utils import plot_model\n",
2961     "plot_model(model_fit, to_file='model_plot.png', \n",
2962     "           show_shapes=True, show_layer_names=True)"
2963    ]
2964   },
2965   {
2966    "cell_type": "code",
2967    "execution_count": null,
2968    "metadata": {},
2969    "outputs": [],
2970    "source": [
2971     "## Check 1: equilibrium input data the same\n",
2972     "\n",
2973     "print(hash2(Et))\n",
2974     "print(hash2(x_train))\n",
2975     "print(hash2(y_train))"
2976    ]
2977   },
2978   {
2979    "cell_type": "code",
2980    "execution_count": null,
2981    "metadata": {},
2982    "outputs": [],
2983    "source": [
2984     "## Check 2: Untrained RNN initialized with same weights\n",
2985     "\n",
2986     "hash2(model_fit.get_weights())"
2987    ]
2988   },
2989   {
2990    "cell_type": "code",
2991    "execution_count": null,
2992    "metadata": {
2993     "id": "ClBMYe8Lqr7P"
2994    },
2995    "outputs": [],
2996    "source": [
2997     "model_predict=create_RNN_2(hidden_units=hidden_units, dense_units=dense_units,  \n",
2998     "                        input_shape=(hours,features),stateful = False,\n",
2999     "                        return_sequences=True,\n",
3000     "                        activation=activation,dense_layers=dense_layers)\n",
3001     "# model_predict=create_RNN_sequences(hidden_units=1, dense_units=1, input_shape=(hours,1), \n",
3002     "#                        activation=['linear', 'linear'])\n",
3003     "# print('model_predict input shape',Et.shape,'output shape',model_predict(Et).shape)\n",
3004     "print(model_predict.summary())\n",
3005     "from keras.utils.vis_utils import plot_model\n",
3006     "plot_model(model_predict, to_file='model_plot.png', \n",
3007     "           show_shapes=True, show_layer_names=True)"
3008    ]
3009   },
3010   {
3011    "cell_type": "code",
3012    "execution_count": null,
3013    "metadata": {},
3014    "outputs": [],
3015    "source": [
3016     "## Check 3: Second model initialization same weights\n",
3017     "\n",
3018     "hash2(model_predict.get_weights())"
3019    ]
3020   },
3021   {
3022    "cell_type": "code",
3023    "execution_count": null,
3024    "metadata": {
3025     "id": "4U0kTEiksNZs"
3026    },
3027    "outputs": [],
3028    "source": [
3029     "w_initial=[np.array([[1.-np.exp(-0.1)]]), np.array([[np.exp(-0.1)]]), np.array([0.]),np.array([[1.0]]),np.array([-1.0])]\n",
3030     "w=model_fit.get_weights()\n",
3031     "for i in range(len(w)):\n",
3032     "  print('weight',i,'shape',w[i].shape,'ndim',w[i].ndim,'given',w_initial[i].shape)\n",
3033     "  for j in range(w[i].shape[0]):\n",
3034     "    if w[i].ndim==2:\n",
3035     "      for k in range(w[i].shape[1]):\n",
3036     "        w[i][j][k]=w_initial[i][0][0]/w[i].shape[0]\n",
3037     "    else:\n",
3038     "      w[i][j]=w_initial[i][0]\n",
3039     "model_fit.set_weights(w)"
3040    ]
3041   },
3042   {
3043    "cell_type": "code",
3044    "execution_count": null,
3045    "metadata": {},
3046    "outputs": [],
3047    "source": [
3048     "## Check 4: weights the same after this step \n",
3049     "\n",
3050     "print(hash2(model_fit.get_weights()))\n",
3051     "print(hash2(x_train))\n",
3052     "print(hash2(y_train))"
3053    ]
3054   },
3055   {
3056    "cell_type": "code",
3057    "execution_count": null,
3058    "metadata": {},
3059    "outputs": [],
3060    "source": [
3061     "model_fit.fit(x_train, y_train, epochs=5000, verbose=0,batch_size=samples)\n",
3062     "w_fitted=model_fit.get_weights()\n",
3063     "for i in range(len(w)):\n",
3064     "  print('weight',i,' exact:',w_exact[i],':  initial:',w_initial[i],' fitted:',w_fitted[i])"
3065    ]
3066   },
3067   {
3068    "cell_type": "code",
3069    "execution_count": null,
3070    "metadata": {},
3071    "outputs": [],
3072    "source": [
3073     "## Check 5: Weights the same after fitting\n",
3074     "\n",
3075     "hash2(model_fit.get_weights())"
3076    ]
3077   },
3078   {
3079    "cell_type": "code",
3080    "execution_count": null,
3081    "metadata": {
3082     "id": "o10lIOl4sndv"
3083    },
3084    "outputs": [],
3085    "source": [
3086     "# evaluate model\n",
3087     "model_predict.set_weights(w_fitted)\n",
3088     "x_input=np.reshape(Et,(1, hours, 2))\n",
3089     "y_output = model_predict.predict(x_input)\n",
3090     "print('x_input.shape=',x_input.shape,'y_output.shape=',y_output.shape)\n",
3091     "print(shift)\n",
3092     "m = np.reshape(y_output,hours)\n",
3093     "# print('weights=',w) ## Note for Jan - this should be w_fitted right?\n",
3094     "print('weights=',w_fitted) \n",
3095     "if scale:\n",
3096     "    print('scaling')\n",
3097     "    m = scalery.inverse_transform(m)\n",
3098     "m = np.reshape(m,hours)\n",
3099     "hour=np.array(range(hours))\n",
3100     "title=\"First RNN forecast\"\n",
3101     "plt.figure(figsize=(16,4))\n",
3102     "plt.plot(hour,Ed,linestyle='--',c='r',label='Drying equilibrium')\n",
3103     "plt.plot(hour,Ew,linestyle='--',c='b',label='Wetting equilibrium')\n",
3104     "# print(len(hour),len(m_f))\n",
3105     "plt.scatter(hour,data,c='b',label='data=10-h fuel data')\n",
3106     "if m is not None:\n",
3107     "    plt.plot(hour[:h2],m[:h2],linestyle='-',c='k',label='m=filtered')\n",
3108     "    plt.plot(hour[h2:hours],m[h2:hours],linestyle='-',c='r',label='m=forecast')\n",
3109     "plt.title(title) \n",
3110     "plt.legend()\n"
3111    ]
3112   },
3113   {
3114    "cell_type": "code",
3115    "execution_count": null,
3116    "metadata": {},
3117    "outputs": [],
3118    "source": [
3119     "rnn_orig={'title':'RNN fitting and prediction - original','hours':hours,'h2':h2,'Ed':Ed,'Ew':Ew,'rain':rain,\n",
3120     "            'fm':obs_data,'m':m}\n",
3121     "# 'w_exact':w_exact,'w_initial':w_initial,'w_fitted':w_fitted\n",
3122     "to_json(rnn_orig,'rnn_orig.json')"
3123    ]
3124   },
3125   {
3126    "cell_type": "code",
3127    "execution_count": null,
3128    "metadata": {},
3129    "outputs": [],
3130    "source": [
3131     "## Check 6: Final Model Hash\n",
3132     "\n",
3133     "hash2(rnn_orig['m'])"
3134    ]
3135   },
3136   {
3137    "cell_type": "code",
3138    "execution_count": null,
3139    "metadata": {
3140     "id": "mrWioCJVuU-G"
3141    },
3142    "outputs": [],
3143    "source": [
3144     "# plot subinterval only\n",
3145     "def plot_int(lb=0,ub=hours,title=\"RNN Prediction\"):\n",
3146     "  hour=np.array(range(hours))\n",
3147     "  plt.figure(figsize=(16,4))\n",
3148     "  plt.plot(hour[lb:ub],Ed[lb:ub],linestyle='--',c='r',label='Drying equilibrium')\n",
3149     "  plt.plot(hour[lb:ub],Ew[lb:ub],linestyle='--',c='b',label='Wetting equilibrium')\n",
3150     "  plt.plot(hour[lb:ub],data[lb:ub],linestyle='-',c='b',label='RAWS fuel moisture data')\n",
3151     "  if lb <= h2:\n",
3152     "    ub1 = min(h2,ub)\n",
3153     "    plt.plot(hour[lb:ub1],m[lb:ub1],linestyle='-',c='k',label='Fuel moisture fitted')\n",
3154     "  if ub >= h2:\n",
3155     "    lb1 = max(h2,lb)\n",
3156     "    plt.plot(hour[lb1:ub],m[lb1:ub],linestyle='-',c='r',label='Fuel moisture prediction')\n",
3157     "  plt.title(title) \n",
3158     "  plt.legend()"
3159    ]
3160   },
3161   {
3162    "cell_type": "code",
3163    "execution_count": null,
3164    "metadata": {
3165     "id": "qmGPeG61uqGI"
3166    },
3167    "outputs": [],
3168    "source": [
3169     "plot_int(0,600,title='RNN fitting and prediction')  # again the whole thing"
3170    ]
3171   },
3172   {
3173    "cell_type": "code",
3174    "execution_count": null,
3175    "metadata": {
3176     "id": "SwnOSJlOuvAA"
3177    },
3178    "outputs": [],
3179    "source": [
3180     "plot_int(0,300,title='RNN Fitting') "
3181    ]
3182   },
3183   {
3184    "cell_type": "code",
3185    "execution_count": null,
3186    "metadata": {
3187     "id": "EqCZD7uCvDrS"
3188    },
3189    "outputs": [],
3190    "source": [
3191     "plot_int(300,500)"
3192    ]
3193   },
3194   {
3195    "cell_type": "code",
3196    "execution_count": null,
3197    "metadata": {
3198     "id": "hYgLAXpUvSLo"
3199    },
3200    "outputs": [],
3201    "source": [
3202     "plot_int(500,800)"
3203    ]
3204   },
3205   {
3206    "cell_type": "markdown",
3207    "metadata": {
3208     "id": "gVQxv9Blc91t"
3209    },
3210    "source": [
3211     "### 4.4 A comment on the information flow in the Kalman filter and in neural networks"
3212    ]
3213   },
3214   {
3215    "cell_type": "markdown",
3216    "metadata": {
3217     "id": "_g_OTEg6ePb9"
3218    },
3219    "source": [
3220     "## 5. Conclusion"
3221    ]
3222   },
3223   {
3224    "cell_type": "markdown",
3225    "metadata": {
3226     "id": "aNxw7xI3FqFt"
3227    },
3228    "source": [
3229     "We have shown how to combine a model and data for improved forecasting of fuel moisture from weather forecast using the Kalman filter. Augmenting the filter state by a model parameter and joint estimation of augmented state resulted in an improvement of the forecast."
3230    ]
3231   },
3232   {
3233    "cell_type": "markdown",
3234    "metadata": {
3235     "id": "IWpmDwUPGElR"
3236    },
3237    "source": [
3238     "## Contributions of authors "
3239    ]
3240   },
3241   {
3242    "cell_type": "markdown",
3243    "metadata": {
3244     "id": "jujW1VFgGOCn"
3245    },
3246    "source": [
3247     "Not applicable."
3248    ]
3249   },
3250   {
3251    "cell_type": "markdown",
3252    "metadata": {
3253     "id": "HWslw7HmGZmP"
3254    },
3255    "source": [
3256     "## Acknowledgements"
3257    ]
3258   },
3259   {
3260    "cell_type": "markdown",
3261    "metadata": {
3262     "id": "xubqDAV2GjkZ"
3263    },
3264    "source": [
3265     "This Math Clinic was sponsored by the team of investigators of the NASA grant no. 80NSSC19K1091 *Coupled Interactive Forecasting of Weather, Fire Behavior, and Smoke Impact for Improved Wildland Fire Decision Making* under the NASA ROSES18 Disasters program. The author would like to thank Brian Zhang from the Math Clinic class for bringing the reference van der Kamp et al. (2017) to his attention."
3266    ]
3267   },
3268   {
3269    "cell_type": "markdown",
3270    "metadata": {
3271     "id": "ZsNZxOv7c91t"
3272    },
3273    "source": [
3274     "## References"
3275    ]
3276   },
3277   {
3278    "cell_type": "markdown",
3279    "metadata": {
3280     "id": "vFY-iS1Wc91t"
3281    },
3282    "source": [
3283     "J. Mandel, S. Amram, J. D. Beezley, G. Kelman, A. K. Kochanski, V. Y. Kondratenko, B. H. Lynn, B. Regev, and M. Vejmelka. *Recent advances and applications of WRF-SFIRE.* Natural Hazards and Earth System Science, 14(10):2829–2845, 2014. [doi:10.5194/nhessd-2-1759-2014](https://doi.org/10.5194/nhessd-2-1759-2014)\n",
3284     "\n",
3285     "R. E. Kalman. *A new approach to linear filtering and prediction problems.* Transactions of the ASME – Journal of Basic Engineering, Series D, 82:35–45, 1960. [doi:10.1115/1.3662552](https://doi.org/10.1115/1.3662552)\n",
3286     "\n",
3287     "E. Kalnay. *Atmospheric Modeling, Data Assimilation and Predictability.* Cambridge University Press, 2003. [doi:10.1017/CBO9780511802270](https://doi.org/10.1017/CBO9780511802270)\n",
3288     "\n",
3289     "D. W. van der Kamp, R. D. Moore, and I. G. McKendry. *A model for simulating the moisture content of standardized fuel sticks of various sizes.* Agricultural and Forest Meteorology, 236:123–134, 2017. [doi:10.1016/j.agrformet.2017.01.013](https://doi.org/10.1016/j.agrformet.2017.01.013)\n",
3290     "\n",
3291     "S. F. Schmidt. *Application of state-space methods to navigation problems.* volume 3 of Advances in Control Systems, C. T.  Leondes, ed., pages 293–340. Elsevier, 1966. [doi:10.1016/B978-1-4831-6716-9.50011-4](https://doi.org/10.1016/B978-1-4831-6716-9.50011-4)\n",
3292     "\n",
3293     "M. Vejmelka, A. K. Kochanski, and J. Mandel. *Data assimilation of dead fuel moisture observations from remote automatic weather stations.* International Journal of Wildland Fire, 25:558– 568, 2016. [doi:10.1071/WF14085](https://doi.org/10.1071/WF14085)\n"
3294    ]
3295   }
3296  ],
3297  "metadata": {
3298   "accelerator": "GPU",
3299   "colab": {
3300    "collapsed_sections": [],
3301    "provenance": []
3302   },
3303   "gpuClass": "standard",
3304   "kernelspec": {
3305    "display_name": "Python 3 (ipykernel)",
3306    "language": "python",
3307    "name": "python3"
3308   },
3309   "language_info": {
3310    "codemirror_mode": {
3311     "name": "ipython",
3312     "version": 3
3313    },
3314    "file_extension": ".py",
3315    "mimetype": "text/x-python",
3316    "name": "python",
3317    "nbconvert_exporter": "python",
3318    "pygments_lexer": "ipython3",
3319    "version": "3.10.9"
3320   }
3321  },
3322  "nbformat": 4,
3323  "nbformat_minor": 4
3324 }