fmda/presentations/fmda_rnn_test.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     "# Inside 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": "fk72YB2mjuGk"
  28    },
  29    "source": [
  30     "## Building and evaluating RNN"
  31    ]
  32   },
  33   {
  34    "cell_type": "markdown",
  35    "metadata": {
  36     "id": "svt8wUAsSA67"
  37    },
  38    "source": [
  39     "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",
  40     "\n",
  41     "We'll start by how to evaluate the map, then actually create it later.\n",
  42     "\n",
  43     "The following is based on https://machinelearningmastery.com/understanding-simple-recurrent-neural-networks-in-keras/\n"
  44    ]
  45   },
  46   {
  47    "cell_type": "code",
  48    "execution_count": null,
  49    "metadata": {
  50     "id": "H3RTQCDG9q-4"
  51    },
  52    "outputs": [],
  53    "source": [
  54     "import numpy as np\n",
  55     "import tensorflow as tf\n",
  56     "from keras.models import Sequential\n",
  57     "from keras.layers import Dense, SimpleRNN\n",
  58     "from keras.utils.vis_utils import plot_model\n",
  59     "from sklearn.preprocessing import MinMaxScaler\n",
  60     "from sklearn.metrics import mean_squared_error\n",
  61     "import math\n",
  62     "import matplotlib.pyplot as plt\n",
  63     "import tensorflow as tf\n",
  64     "import keras.backend as K"
  65    ]
  66   },
  67   {
  68    "cell_type": "code",
  69    "execution_count": null,
  70    "metadata": {
  71     "id": "pcIU5lWhVEAy"
  72    },
  73    "outputs": [],
  74    "source": [
  75     "def create_RNN(hidden_units, dense_units, input_shape, activation):\n",
  76     "    inputs = tf.keras.Input(shape=input_shape)\n",
  77     "    # https://stackoverflow.com/questions/43448029/how-can-i-print-the-values-of-keras-tensors\n",
  78     "    # inputs2 = K.print_tensor(inputs, message='inputs = ')  # change allso inputs to inputs2 below, must be used\n",
  79     "    x = tf.keras.layers.SimpleRNN(hidden_units, input_shape=input_shape,\n",
  80     "                        activation=activation[0])(inputs)\n",
  81     "    outputs = tf.keras.layers.Dense(dense_units, activation=activation[1])(x)\n",
  82     "    model = tf.keras.Model(inputs=inputs, outputs=outputs)\n",
  83     "    model.compile(loss='mean_squared_error', optimizer='adam')\n",
  84     "    return model"
  85    ]
  86   },
  87   {
  88    "cell_type": "code",
  89    "execution_count": null,
  90    "metadata": {
  91     "id": "lphzeucMfI7L"
  92    },
  93    "outputs": [],
  94    "source": [
  95     "# Demo example\n",
  96     "hidden=5\n",
  97     "features=2\n",
  98     "timesteps=3\n",
  99     "demo_model = create_RNN(hidden_units=hidden, dense_units=1, \n",
 100     "                        input_shape=(timesteps,features), \n",
 101     "                        activation=['linear', 'linear'])\n",
 102     "print(demo_model.summary())\n",
 103     "w = demo_model.get_weights()\n",
 104     "#print(len(w),' weight arrays:',w)\n",
 105     "wname=('wx','wh','bh','wy','by','wz','bz')\n",
 106     "for i in range(len(w)):\n",
 107     "  print(i,':',wname[i],'shape=',w[i].shape)\n",
 108     "wx, wh, bh, wy, by = w\n",
 109     "plot_model(demo_model, to_file='model_plot.png', \n",
 110     "  show_shapes=True, show_layer_names=True,\n",
 111     "  expand_nested=True,)"
 112    ]
 113   },
 114   {
 115    "cell_type": "markdown",
 116    "metadata": {
 117     "id": "uFVQdXL0fODX"
 118    },
 119    "source": [
 120     "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",
 121     "\n",
 122     "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",
 123     "<img src=\"https://machinelearningmastery.com/wp-content/uploads/2021/09/rnnCode1.png\">"
 124    ]
 125   },
 126   {
 127    "cell_type": "markdown",
 128    "metadata": {
 129     "id": "TcxvQPo1hYip"
 130    },
 131    "source": [
 132     "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",
 133     "\\begin{align*}\n",
 134     "h\\left(  0\\right)  = &0  \\\\\n",
 135     "h\\left(  k+1\\right)  =& \n",
 136     "x\\left(  k\\right) w_{x}\n",
 137     "  +h(k) w_{h}  + b_{h}\\\\\n",
 138     "o(k+1)=& h(k+1)w_{y} + b_y\n",
 139     "\\end{align*}"
 140    ]
 141   },
 142   {
 143    "cell_type": "code",
 144    "execution_count": null,
 145    "metadata": {
 146     "id": "fqL5TEfpml7q"
 147    },
 148    "outputs": [],
 149    "source": [
 150     "# Reshape the input to sample_size x time_steps x features \n",
 151     "samples=4   # number of samples\n",
 152     "outputs = 1 # number of outputs\n",
 153     "x_tf = tf.reshape(tf.range(samples*timesteps*features),[samples,timesteps,features]) \n",
 154     "# print('test input x=',x)\n",
 155     "print('model.predict start')\n",
 156     "y_pred_model = demo_model.predict(x_tf)\n",
 157     "print('model.predict end')\n",
 158     "x = x_tf.numpy()  # without converting x from tensor we get in loops zero\n",
 159     "\n",
 160     "o3=np.zeros([samples,outputs])\n",
 161     "h_3 = np.zeros(hidden)\n",
 162     "for i in range(samples):\n",
 163     "  h_0 = np.zeros(hidden)\n",
 164     "  h_1 = np.dot(x[i,0,:], wx) + np.dot(h_0,wh) + bh\n",
 165     "  h_2 = np.dot(x[i,1,:], wx) + np.dot(h_1,wh) + bh\n",
 166     "  h_3 = np.dot(x[i,2,:], wx) + np.dot(h_2,wh) + bh\n",
 167     "  o3[i,:] = np.dot(h_3, wy) + by\n",
 168     "# in loopa\n",
 169     "timesteps=3\n",
 170     "\n",
 171     "o = np.zeros([samples,outputs])\n",
 172     "for i in range(samples):\n",
 173     "    h = np.zeros(hidden)\n",
 174     "    for j in range(timesteps):\n",
 175     "        hout = np.zeros(hidden)\n",
 176     "        for k in range(hidden):\n",
 177     "            hout[k] = bh[k]\n",
 178     "            for l in range(features):\n",
 179     "                hout[k] += x[i,j,l]*wx[l,k] \n",
 180     "            for l in range(hidden):\n",
 181     "                hout[k] += h[l]*wh[l,k]\n",
 182     "        h = hout.copy()\n",
 183     "    o[i,:] = by\n",
 184     "    for l in range(outputs):\n",
 185     "         for k in range(hidden):\n",
 186     "            # print('i=',i,'k=',k,'l=',l)\n",
 187     "            o[i,l] += h[k]*wy[k,l]\n",
 188     "\n",
 189     "print(\"Prediction from network \", y_pred_model)\n",
 190     "print(\"Prediction from our computation \", o3)\n",
 191     "print(\"Prediction computed by loops \", o)\n",
 192     "if np.allclose(y_pred_model,o3) and np.allclose(y_pred_model,o):\n",
 193     "    print('The result is the same.')"
 194    ]
 195   },
 196   {
 197    "cell_type": "code",
 198    "execution_count": null,
 199    "metadata": {},
 200    "outputs": [],
 201    "source": [
 202     "input(\"Press Enter to continue...\")"
 203    ]
 204   },
 205   {
 206    "cell_type": "markdown",
 207    "metadata": {
 208     "id": "A9ZpmNcdRpmp"
 209    },
 210    "source": [
 211     "## Fuel moisture models\n",
 212     "\n",
 213     "\n"
 214    ]
 215   },
 216   {
 217    "cell_type": "markdown",
 218    "metadata": {
 219     "id": "eZL8gN7ISGVh"
 220    },
 221    "source": [
 222     "### A simple fuel moisture model"
 223    ]
 224   },
 225   {
 226    "cell_type": "markdown",
 227    "metadata": {
 228     "id": "1XvOC4kYSQgH"
 229    },
 230    "source": [
 231     "First consider a simplified fuel moisture model without considering the effect of rain.\n",
 232     "The evolution of fuel moisture content $m(t)$ is modeled by the time-lag differential equation on interval $\\left[\n",
 233     "t_{0},t_{1}\\right]  $,\n",
 234     "$$\n",
 235     "\\frac{dm}{dt}=\\frac{E-m(t)}{T},\\quad m(t_{0})=m_{0}.\n",
 236     "$$\n",
 237     "where the initial fuel moisture content $m_{0}=m\\left(  t_{0}\\right)  $ is the\n",
 238     "input, and $m_{1}=m(t_{1})$ is the output. Tnus, $m_1=F(m_0)$. The parameters of the model are the\n",
 239     "fuel moisture equilibrium $E$, assumed to be constant over the interval $\\left[\n",
 240     "t_{0},t_{1}\\right]  $, NS the characteristic decay time $T$. \n",
 241     "\n",
 242     "We can build the general model later by calling this simple model with different\n",
 243     "equilibria and time constants (drying, wetting, rain).\n",
 244     "\n",
 245     "Since $E$ is constant in time, the solution can be found\n",
 246     "analytically,\n",
 247     "$$\n",
 248     "m\\left(  t\\right)  =E+\\left(  m_{0}-E\\right)  e^{-t/T}%\n",
 249     "$$\n",
 250     "For convenience, we use $T_{1}=1/T$ instead of $T$, and the model becomes\n",
 251     "$$\n",
 252     "m_{1}=E+\\left(  m_{0}-E\\right)  e^{-\\left(  t_{1}-t_{0}\\right)  T_{1}}%\n",
 253     "$$\n",
 254     "In the extended Kalman filter, we will need the partial derivatives of $m_{1}$\n",
 255     "with respect to the input and the parameters. Compute\n",
 256     "$$\n",
 257     "\\frac{dm_{1}}{d_{m0}}=e^{-\\left(  t_{1}-t_{0}\\right)  T_{1}}\n",
 258     "$$\n",
 259     "$$\n",
 260     "\\frac{dm_{1}}{dE}=1-e^{-\\left(  t_{1}-t_{0}\\right)  T_{1}}\n",
 261     "$$\n",
 262     "$$\n",
 263     "\\frac{dm_{1}}{dT_{1}}=-\\left(  m_{0}-E\\right)  \\left(  t_{1}-t_{0}\\right)\n",
 264     "e^{-\\left(  t_{1}-t_{0}\\right)  T_{1}}\n",
 265     "$$\n",
 266     "At the moment, we need only ${dm_{1}}/{dm_{0}}$ but we put in the code all partials for possible use in future.\n"
 267    ]
 268   },
 269   {
 270    "cell_type": "code",
 271    "execution_count": null,
 272    "metadata": {
 273     "id": "wuVIAGLiSeR8"
 274    },
 275    "outputs": [],
 276    "source": [
 277     "import numpy as np\n",
 278     "def model_decay(m0,E,partials=0,T1=0.1,tlen=1):  \n",
 279     "  # Arguments: \n",
 280     "  #   m0          fuel moisture content at start dimensionless, unit (1)\n",
 281     "  #   E           fuel moisture eqilibrium (1)\n",
 282     "  #   partials=0: return m1 = fuel moisture contents after time tlen (1)\n",
 283     "  #           =1: return m1, dm0/dm0 \n",
 284     "  #           =2: return m1, dm1/dm0, dm1/dE\n",
 285     "  #           =3: return m1, dm1/dm0, dm1/dE dm1/dT1   \n",
 286     "  #   T1          1/T, where T is the time constant approaching the equilibrium\n",
 287     "  #               default 0.1/hour\n",
 288     "  #   tlen        the time interval length, default 1 hour\n",
 289     "\n",
 290     "  exp_t = np.exp(-tlen*T1)                  # compute this subexpression only once\n",
 291     "  m1 = E + (m0 - E)*exp_t                   # the solution at end\n",
 292     "  if partials==0:\n",
 293     "    return m1\n",
 294     "  dm1_dm0 = exp_t\n",
 295     "  if partials==1:\n",
 296     "    return m1, dm1_dm0          # return value and Jacobian\n",
 297     "  dm1_dE = 1 - exp_t      \n",
 298     "  if partials==2:\n",
 299     "     return m1, dm1_dm0, dm1_dE \n",
 300     "  dm1_dT1 = -(m0 - E)*tlen*exp_t            # partial derivative dm1 / dT1\n",
 301     "  if partials==3:\n",
 302     "    return m1, dm1_dm0, dm1_dE, dm1_dT1       # return value and all partial derivatives wrt m1 and parameters\n",
 303     "  raise('Bad arg partials')\n",
 304     "  "
 305    ]
 306   },
 307   {
 308    "cell_type": "markdown",
 309    "metadata": {
 310     "id": "dOARZlj-RUCi"
 311    },
 312    "source": [
 313     "### Fuel moisture model with drying equilibrium, wetting equilibrium, and rain"
 314    ]
 315   },
 316   {
 317    "cell_type": "markdown",
 318    "metadata": {
 319     "id": "AJp6FTpTSx5B"
 320    },
 321    "source": [
 322     "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",
 323     "asymptotically, when the rain intensity $r$ (mm/h) is\n",
 324     "large:\n",
 325     "$$\n",
 326     "\\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",
 327     "\\right)  \\right),\\ \\text{if}\\ r>r_0, \n",
 328     "$$\n",
 329     "where $r_0$ is the threshold rain intensity below which no perceptible\n",
 330     "wetting occurs, and $r_{\\mathrm{s}}$ is the saturation rain\n",
 331     "intensity. At the saturation rain intensity, $1-1/e\\approx 0.63$ of\n",
 332     "the maximal rain-wetting rate is achieved. For 10h fuel, the model takes $S=250\\,{\\%}$,\n",
 333     "$t_{\\mathrm{r}}=14$h, $r_0=0.05$mm/h and\n",
 334     "$r_{\\mathrm{s}}=8$mm/h. "
 335    ]
 336   },
 337   {
 338    "cell_type": "code",
 339    "execution_count": null,
 340    "metadata": {
 341     "id": "ITsKE0psRblG"
 342    },
 343    "outputs": [],
 344    "source": [
 345     "### Define model function with drying, wetting, and rain equilibria\n",
 346     "\n",
 347     "# Parameters\n",
 348     "r0 = 0.05                                   # threshold rainfall [mm/h]\n",
 349     "rs = 8.0                                    # saturation rain intensity [mm/h]\n",
 350     "Tr = 14.0                                   # time constant for rain wetting model [h]\n",
 351     "S = 250                                     # saturation intensity [dimensionless]\n",
 352     "T = 10.0                                    # time constant for wetting/drying\n",
 353     "\n",
 354     "def model_moisture(m0,Eqd,Eqw,r,t,partials=0,T=10.0,tlen=1.0):\n",
 355     "    # arguments:\n",
 356     "    # m0         starting fuel moistureb (%s\n",
 357     "    # Eqd        drying equilibrium      (%) \n",
 358     "    # Eqw        wetting equilibrium     (%)\n",
 359     "    # r          rain intensity          (mm/h)\n",
 360     "    # t          time\n",
 361     "    # partials = 0, 1, 2\n",
 362     "    # returns: same as model_decay\n",
 363     "    #   if partials==0: m1 = fuel moisture contents after time 1 hour\n",
 364     "    #              ==1: m1, dm1/dm0 \n",
 365     "    #              ==2: m1, dm1/dm0, dm1/dE  \n",
 366     "    \n",
 367     "    if r > r0:\n",
 368     "        # print('raining')\n",
 369     "        E = S\n",
 370     "        T1 =  (1.0 - np.exp(- (r - r0) / rs)) / Tr\n",
 371     "    elif m0 <= Eqw: \n",
 372     "        # print('wetting')\n",
 373     "        E=Eqw\n",
 374     "        T1 = 1.0/T\n",
 375     "    elif m0 >= Eqd:\n",
 376     "        # print('drying')\n",
 377     "        E=Eqd\n",
 378     "        T1 = 1.0/T\n",
 379     "    else: # no change'\n",
 380     "        E = m0\n",
 381     "        T1=0.0\n",
 382     "    exp_t = np.exp(-tlen*T1)\n",
 383     "    m1 = E + (m0 - E)*exp_t  \n",
 384     "    dm1_dm0 = exp_t\n",
 385     "    dm1_dE = 1 - exp_t\n",
 386     "    #if t>=933 and t < 940:\n",
 387     "    #  print('t,Eqw,Eqd,r,T1,E,m0,m1,dm1_dm0,dm1_dE',\n",
 388     "    #        t,Eqw,Eqd,r,T1,E,m0,m1,dm1_dm0,dm1_dE)   \n",
 389     "    if partials==0: \n",
 390     "        return m1\n",
 391     "    if partials==1:\n",
 392     "        return m1, dm1_dm0\n",
 393     "    if partials==2:\n",
 394     "        return m1, dm1_dm0, dm1_dE\n",
 395     "    raise('bad partials')"
 396    ]
 397   },
 398   {
 399    "cell_type": "markdown",
 400    "metadata": {
 401     "id": "AkyiGlZF0WrM"
 402    },
 403    "source": [
 404     "#### Training and forecasting with the RNN"
 405    ]
 406   },
 407   {
 408    "cell_type": "markdown",
 409    "metadata": {
 410     "id": "e265QFwlw22_"
 411    },
 412    "source": [
 413     "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,:]`."
 414    ]
 415   },
 416   {
 417    "cell_type": "code",
 418    "execution_count": null,
 419    "metadata": {
 420     "id": "owX7OTg-RnMA"
 421    },
 422    "outputs": [],
 423    "source": [
 424     "def staircase(x,y,timesteps,trainsteps,return_sequences=False):\n",
 425     "  # x [trainsteps+forecaststeps,features]    all inputs\n",
 426     "  # y [trainsteps,outputs]\n",
 427     "  # timesteps: split x and y into samples length timesteps, shifted by 1\n",
 428     "  # trainsteps: number of timesteps to use for training, no more than y.shape[0]\n",
 429     "  print('shape x = ',x.shape)\n",
 430     "  print('shape y = ',y.shape)\n",
 431     "  print('timesteps=',timesteps)\n",
 432     "  print('trainsteps=',trainsteps)\n",
 433     "  outputs = y.shape[1]\n",
 434     "  features = x.shape[1]\n",
 435     "  forecaststeps = x.shape[0]-trainsteps\n",
 436     "  samples = trainsteps-timesteps+1\n",
 437     "  print('staircase: samples=',samples,'timesteps=',timesteps,'features=',features)\n",
 438     "  x_train = np.empty([samples, timesteps, features])\n",
 439     "  print('return_sequences=',return_sequences)\n",
 440     "  if return_sequences:\n",
 441     "    print('returning all timesteps in a sample')\n",
 442     "    y_train = np.empty([samples, timesteps, outputs])  # all\n",
 443     "    for i in range(samples):\n",
 444     "      for k in range(timesteps):\n",
 445     "        for j in range(features):\n",
 446     "          x_train[i,k,j] = x[i+k,j]\n",
 447     "        for j in range(outputs):\n",
 448     "          y_train[i,k,j] = y[i+k,j]\n",
 449     "  else:\n",
 450     "    print('returning only the last timestep in a sample')\n",
 451     "    y_train = np.empty([samples, outputs])\n",
 452     "    for i in range(samples):\n",
 453     "      for j in range(features):\n",
 454     "        for k in range(timesteps):\n",
 455     "          x_train[i,k,j] = x[i+k,j]\n",
 456     "      for j in range(outputs):\n",
 457     "        y_train[i,j] = y[i+timesteps-1,j]\n",
 458     "\n",
 459     "  return x_train, y_train"
 460    ]
 461   },
 462   {
 463    "cell_type": "code",
 464    "execution_count": null,
 465    "metadata": {
 466     "id": "FzOotSFf-tPR"
 467    },
 468    "outputs": [],
 469    "source": [
 470     "def seq2batches(x,y,timesteps,trainsteps):\n",
 471     "  # x [trainsteps+forecaststeps,features]    all inputs\n",
 472     "  # y [trainsteps,outputs]\n",
 473     "  # timesteps: split x and y into samples length timesteps, shifted by 1\n",
 474     "  # trainsteps: number of timesteps to use for training, no more than y.shape[0]\n",
 475     "  print('shape x = ',x.shape)\n",
 476     "  print('shape y = ',y.shape)\n",
 477     "  print('timesteps=',timesteps)\n",
 478     "  print('trainsteps=',trainsteps)\n",
 479     "  outputs = y.shape[1]\n",
 480     "  features = x.shape[1]\n",
 481     "  samples= trainsteps - timesteps + 1\n",
 482     "  print('samples=',samples)\n",
 483     "  x_train = np.empty([samples, timesteps, features])\n",
 484     "  y_train = np.empty([samples, timesteps, outputs])  # only the last\n",
 485     "  print('samples=',samples,' timesteps=',timesteps,\n",
 486     "        ' features=',features,' outputs=',outputs)\n",
 487     "  for i in range(samples):\n",
 488     "    for k in range(timesteps):\n",
 489     "      for j in range(features):\n",
 490     "        x_train[i,k,j] = x[i+k,j]\n",
 491     "      for j in range(outputs):\n",
 492     "        y_train[i,k,j] = y[i+k,j]  # return sequences\n",
 493     "  return x_train, y_train"
 494    ]
 495   },
 496   {
 497    "cell_type": "code",
 498    "execution_count": null,
 499    "metadata": {
 500     "id": "Kg7wSrkk-HrE"
 501    },
 502    "outputs": [],
 503    "source": [
 504     "print('test preprocessing for RNN')\n",
 505     "trainsteps=5\n",
 506     "features=1\n",
 507     "outputs=1\n",
 508     "timesteps=3\n",
 509     "x = tf.reshape(tf.range(trainsteps*features),[trainsteps,features])\n",
 510     "y = tf.reshape(tf.range(trainsteps*outputs),[trainsteps,outputs])\n",
 511     "print('x=',x)\n",
 512     "print('y=',y)\n",
 513     "x_train, y_train = staircase(x,y,timesteps,trainsteps)\n",
 514     "print('x_train=',x_train)\n",
 515     "print('y_train=',y_train)\n",
 516     "x_train, y_train = seq2batches(x,y,timesteps,trainsteps)\n",
 517     "print('x_train=',x_train)\n",
 518     "print('y_train=',y_train)"
 519    ]
 520   },
 521   {
 522    "cell_type": "code",
 523    "execution_count": null,
 524    "metadata": {
 525     "id": "rHWxqARpSO_f"
 526    },
 527    "outputs": [],
 528    "source": [
 529     "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) "
 530    ]
 531   },
 532   {
 533    "cell_type": "code",
 534    "execution_count": null,
 535    "metadata": {
 536     "id": "vfcxN9JCH5Ku"
 537    },
 538    "outputs": [],
 539    "source": [
 540     "scale=False\n",
 541     "# transform as 2D, (timesteps, features) and (timesteps, outputs)\n",
 542     "Et = np.reshape(E,[E.shape[0],1])\n",
 543     "datat = np.reshape(data,[data.shape[0],1])\n",
 544     "if scale:\n",
 545     "  scalerx = MinMaxScaler()\n",
 546     "  scalerx.fit(Et)\n",
 547     "  Et = scalerx.transform(Et)\n",
 548     "  scalery = MinMaxScaler()\n",
 549     "  scalery.fit(datat)\n",
 550     "  datat = scalery.transform(datat)"
 551    ]
 552   },
 553   {
 554    "cell_type": "code",
 555    "execution_count": null,
 556    "metadata": {
 557     "id": "PaHfJW7mSJE1"
 558    },
 559    "outputs": [],
 560    "source": [
 561     "def create_RNN_2(hidden_units, dense_units, activation, stateful=False, \n",
 562     "                 batch_shape=None, input_shape=None, dense_layers=1,\n",
 563     "                 rnn_layers=1,return_sequences=False,\n",
 564     "                 initial_state=None):\n",
 565     "    if stateful:\n",
 566     "      inputs = tf.keras.Input(batch_shape=batch_shape)\n",
 567     "    else:\n",
 568     "      inputs = tf.keras.Input(shape=input_shape)\n",
 569     "    # https://stackoverflow.com/questions/43448029/how-can-i-print-the-values-of-keras-tensors\n",
 570     "    # inputs2 = K.print_tensor(inputs, message='inputs = ')  # change allso inputs to inputs2 below, must be used\n",
 571     "    x = inputs\n",
 572     "    for i in range(rnn_layers):\n",
 573     "      x = tf.keras.layers.SimpleRNN(hidden_units,activation=activation[0],\n",
 574     "              stateful=stateful,return_sequences=return_sequences)(x\n",
 575     "              # ,initial_state=initial_state\n",
 576     "              )\n",
 577     "    # x = tf.keras.layers.Dense(hidden_units, activation=activation[1])(x)\n",
 578     "    for i in range(dense_layers):\n",
 579     "      x = tf.keras.layers.Dense(dense_units, activation=activation[1])(x)\n",
 580     "    model = tf.keras.Model(inputs=inputs, outputs=x)\n",
 581     "    model.compile(loss='mean_squared_error', optimizer='adam')\n",
 582     "    return model"
 583    ]
 584   },
 585   {
 586    "cell_type": "code",
 587    "execution_count": null,
 588    "metadata": {
 589     "id": "ZjJxHiEVL5sJ"
 590    },
 591    "outputs": [],
 592    "source": [
 593     "# split data\n",
 594     "return_sequences=False\n",
 595     "shift = 0.\n",
 596     "print('shifting inputs by',shift)\n",
 597     "x_train, y_train = staircase(Et+shift,datat+shift,timesteps=5,trainsteps=h2,\n",
 598     "                             return_sequences=return_sequences)\n",
 599     "print('x_train shape=',x_train.shape)\n",
 600     "samples, timesteps, features = x_train.shape\n",
 601     "print('y_train shape=',y_train.shape)\n",
 602     "# the simplest model possible\n",
 603     "activation=['linear','linear']\n",
 604     "hidden_units=1\n",
 605     "dense_units=1\n",
 606     "dense_layers=1\n",
 607     "features=1\n",
 608     "hours=Et.shape[0]\n",
 609     "h0 = tf.convert_to_tensor(datat[:samples],dtype=tf.float32)\n",
 610     "# print('initial state=',h0)\n",
 611     "# statefull model version for traning\n",
 612     "import moisture_rnn\n",
 613     "# model_fit=create_RNN_2(hidden_units=hidden_units, \n",
 614     "model_fit=moisture_rnn.create_RNN_2(hidden_units=hidden_units, \n",
 615     "                        dense_units=dense_units, \n",
 616     "                        batch_shape=(samples,timesteps,features),\n",
 617     "                        stateful=True,\n",
 618     "                        return_sequences=return_sequences,\n",
 619     "                        # initial_state=h0,\n",
 620     "                        activation=activation,\n",
 621     "                        dense_layers=dense_layers)\n",
 622     "# same model stateless for prediction on the entire dataset - to start onlg\n",
 623     "# the real application will switch to prediction after training data end\n",
 624     "# and start from the state there\n",
 625     "print('model_fit input shape',x_train.shape,'output shape',model_fit(x_train).shape)\n",
 626     "from keras.utils.vis_utils import plot_model\n",
 627     "plot_model(model_fit, to_file='model_plot.png', \n",
 628     "           show_shapes=True, show_layer_names=True)"
 629    ]
 630   },
 631   {
 632    "cell_type": "code",
 633    "execution_count": null,
 634    "metadata": {
 635     "id": "ipX9EJqz17Lr"
 636    },
 637    "outputs": [],
 638    "source": [
 639     "model_predict=create_RNN_2(hidden_units=hidden_units, dense_units=dense_units,  \n",
 640     "                        input_shape=(hours,features),stateful = False,\n",
 641     "                        return_sequences=True,\n",
 642     "                        activation=activation,dense_layers=dense_layers)\n",
 643     "# model_predict=create_RNN_sequences(hidden_units=1, dense_units=1, input_shape=(hours,1), \n",
 644     "#                        activation=['linear', 'linear'])\n",
 645     "print('model_predict input shape',Et.shape,'output shape',model_predict(Et).shape)\n",
 646     "print(model_predict.summary())\n",
 647     "from keras.utils.vis_utils import plot_model\n",
 648     "plot_model(model_predict, to_file='model_plot.png', \n",
 649     "           show_shapes=True, show_layer_names=True)"
 650    ]
 651   },
 652   {
 653    "cell_type": "code",
 654    "execution_count": null,
 655    "metadata": {
 656     "id": "dE0OHg0tGVE9"
 657    },
 658    "outputs": [],
 659    "source": [
 660     "# fitting\n",
 661     "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",
 662     "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",
 663     "model_fit.set_weights(w_initial)\n",
 664     "model_fit.fit(x_train, y_train, epochs=1000, verbose=0,batch_size=samples)\n",
 665     "w_fitted=model_fit.get_weights()\n",
 666     "for i in range(len(w)):\n",
 667     "  print('weight',i,' exact:',w_exact[i],':  initial:',w_initial[i],' fitted:',w_fitted[i])"
 668    ]
 669   },
 670   {
 671    "cell_type": "code",
 672    "execution_count": null,
 673    "metadata": {
 674     "id": "-T8lCS6nBHCj"
 675    },
 676    "outputs": [],
 677    "source": [
 678     "def model_eval(w,title):\n",
 679     "  # prediction on the entire dataset from zero state\n",
 680     "  model_predict.set_weights(w)\n",
 681     "  hours=Et.shape[0]\n",
 682     "  print('Et.shape=',Et.shape,'hours=',hours)\n",
 683     "  x_input=np.reshape(Et,(1, hours, 1))\n",
 684     "  y_output = model_predict.predict(x_input)\n",
 685     "  print('x_input.shape=',x_input.shape,'y_output.shape=',y_output.shape)\n",
 686     "  m = np.reshape(y_output,hours) - shift\n",
 687     "  print('weights=',w)\n",
 688     "  if scale:\n",
 689     "    print('scaling')\n",
 690     "    m = scalery.inverse_transform(m)\n",
 691     "  m = np.reshape(m,hours)\n",
 692     "  plot_m(m,title=title)\n",
 693     "  return m"
 694    ]
 695   },
 696   {
 697    "cell_type": "code",
 698    "execution_count": null,
 699    "metadata": {
 700     "id": "R2jkoZlAIaSb"
 701    },
 702    "outputs": [],
 703    "source": [
 704     "m_fitted=model_eval(w_fitted,'RNN prediction with fitted weights')"
 705    ]
 706   },
 707   {
 708    "cell_type": "code",
 709    "execution_count": null,
 710    "metadata": {
 711     "id": "bO1ewTj9gGvg"
 712    },
 713    "outputs": [],
 714    "source": [
 715     "m_exact=model_eval(w_exact,'RNN prediction with exact weights')"
 716    ]
 717   },
 718   {
 719    "cell_type": "code",
 720    "execution_count": null,
 721    "metadata": {
 722     "id": "henChC0cmbuy"
 723    },
 724    "outputs": [],
 725    "source": [
 726     "m_initial=model_eval(w_initial,'RNN prediction with initial weights')"
 727    ]
 728   },
 729   {
 730    "cell_type": "code",
 731    "execution_count": null,
 732    "metadata": {
 733     "id": "PZw7DNQD4Inr"
 734    },
 735    "outputs": [],
 736    "source": [
 737     "out = np.empty((hours,1))\n",
 738     "w=w_exact\n",
 739     "h=0\n",
 740     "for i in range(Et.shape[0]):\n",
 741     "  h=np.dot(Et[i,0],w[0])+np.dot(h,w[1]) + w[2]\n",
 742     "  out[i]=np.dot(h,w[3]) + w[4]\n",
 743     "if scale:\n",
 744     "  print('scaling')\n",
 745     "  out = scalery.inverse_transform(out)\n",
 746     "out=np.reshape(out,hours)\n",
 747     "print('max abs diff',np.max(np.abs(m_exact-out)))\n",
 748     "plot_m(out,title='Hand computed RNN prediction with exact weights')"
 749    ]
 750   },
 751   {
 752    "cell_type": "markdown",
 753    "metadata": {
 754     "id": "6uXVJj9koGF2"
 755    },
 756    "source": [
 757     "### 3.2 Acquisition and preprocessing of real data"
 758    ]
 759   },
 760   {
 761    "cell_type": "markdown",
 762    "metadata": {
 763     "id": "q3BpOBuzc91i"
 764    },
 765    "source": [
 766     "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.  "
 767    ]
 768   },
 769   {
 770    "cell_type": "markdown",
 771    "metadata": {
 772     "id": "c8Y6bL1Yc91i"
 773    },
 774    "source": [
 775     "#### 3.2.1 Acquisition of fuel moisture observations"
 776    ]
 777   },
 778   {
 779    "cell_type": "markdown",
 780    "metadata": {
 781     "id": "0CuXyWBFc91i",
 782     "tags": []
 783    },
 784    "source": [
 785     "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."
 786    ]
 787   },
 788   {
 789    "cell_type": "code",
 790    "execution_count": null,
 791    "metadata": {},
 792    "outputs": [],
 793    "source": [
 794     "import json\n",
 795     "os.chdir('data')"
 796    ]
 797   },
 798   {
 799    "cell_type": "code",
 800    "execution_count": null,
 801    "metadata": {
 802     "id": "LFrlbbMmc91i"
 803    },
 804    "outputs": [],
 805    "source": [
 806     "jfile = 'raws.json'; vars='fuel_moisture'; case = 1\n",
 807     "# jfile = 'raws2.json'; vars='fuel_moisture,precip_accum_one_hour'; case = 2\n",
 808     "def json_w(j,f):\n",
 809     "  print('writing json file',f)\n",
 810     "  json.dump(j,open(f,'w'),indent=4)\n",
 811     "try:\n",
 812     "    #! wget --no-clobber http://math.ucdenver.edu/~jmandel/data/math4779f21/raws.json\n",
 813     "    j = json.load(open(jfile,'r'))\n",
 814     "    print('loaded from ',jfile)\n",
 815     "    # Take the first station in the boulding box that has data between time_start and time_s2.\n",
 816     "    # Then retrieve data for that station between time_start and time_end\n",
 817     "    time_start = j['time_start']      # start of data time series\n",
 818     "    # time_s2    = j['time_s2']         # end of segment to read coordinates\n",
 819     "    time_end  = j['time_end']         # end of data time series\n",
 820     "    meso_ts  = j['meso_ts']           # get meso observations time series\n",
 821     "    obs_lon =   j['obs_lon']          # where we retrieved observations\n",
 822     "    obs_lat =   j['obs_lat']\n",
 823     "except:\n",
 824     "    print(\"can't read\",jfile,', creating')\n",
 825     "    # set up bounds\n",
 826     "    time_start = \"201806010800\"  # June 1 2018 08:00 in format yyyymmddHHMM\n",
 827     "    time_s2    = \"201806010900\"  # June 1 2018 09:00 in format yyyymmddHHMM \n",
 828     "    time_end   = \"201907200900\"  # June 20 2018 09:00 in format yyyymmddHHMM \n",
 829     "    #time_start=  \"201810230100\"\n",
 830     "    #time_s2=  \"201810230300\"\n",
 831     "    #time_end  =  \"201806022300\"\n",
 832     "    !pip install MesoPy\n",
 833     "    from MesoPy import Meso\n",
 834     "    bounding_box = \"-115, 38, -110, 40\"  # min longtitude, latitude\n",
 835     "    meso_token=\"b40cb52cbdef43ef81329b84e8fd874f\"       # you should get your own if you do more of this\n",
 836     "    m = Meso(meso_token)# create a Meso object\n",
 837     "    print('reading MesoWest fuel moisture data')\n",
 838     "    json_w(m.variables(),'variables.json')\n",
 839     "    meso_obss = m.timeseries(time_start, time_s2, bbox=bounding_box, \n",
 840     "                             showemptystations = '0', vars=vars)   # ask the object for data\n",
 841     "    json_w(meso_obss,'meso_obss.json')                        \n",
 842     "    # pick one station and retrieve the whole time series.\n",
 843     "    station=meso_obss['STATION'][0]\n",
 844     "    json_w(station,'station.json')\n",
 845     "    lon,lat = (float(station['LONGITUDE']),float(station['LATITUDE']))\n",
 846     "    print(station['NAME'],'station',station['STID'],'at',lon,lat)\n",
 847     "    e = 0.01   # tolerance\n",
 848     "    bb = '%s, %s, %s, %s' % (lon - e, lat - e, lon + e, lat + e)\n",
 849     "    print('bounding box',bb)\n",
 850     "    meso_ts = m.timeseries(time_start, time_end, bbox=bb, showemptystations = '0', vars=vars)   # ask the object for data\n",
 851     "    json_w(meso_ts,'meso_ts.json')                        \n",
 852     "    obs_lon, obs_lat = (lon, lat)   # remember station coordinates for later\n",
 853     "    j={'time_start':time_start,'time_s2':time_s2,'time_end':time_end,\n",
 854     "       'meso_ts':meso_ts,'obs_lon':obs_lon,'obs_lat':obs_lat}\n",
 855     "    json_w(j,jfile)\n",
 856     "    print('done')"
 857    ]
 858   },
 859   {
 860    "cell_type": "code",
 861    "execution_count": null,
 862    "metadata": {
 863     "id": "3bXopS3btyz0"
 864    },
 865    "outputs": [],
 866    "source": [
 867     "os.chdir('..')\n",
 868     "# process the data retrieved for this station\n",
 869     "# print(json.dumps(meso_ts['STATION'][0], indent=4))\n",
 870     "from datetime import datetime, timedelta, time\n",
 871     "import numpy as np\n",
 872     "import matplotlib.pyplot as plt\n",
 873     "import pytz\n",
 874     "station = meso_ts['STATION'][0]\n",
 875     "time_str  = station['OBSERVATIONS']['date_time']\n",
 876     "obs_time = [datetime.strptime(t, '%Y-%m-%dT%H:%M:%SZ').replace(tzinfo=pytz.UTC) for t in time_str]\n",
 877     "start_time = obs_time[0].replace(minute=0)     # remember obs_time and start_time for later\n",
 878     "end_time = obs_time[-1]\n",
 879     "obs_data = np.array(station['OBSERVATIONS'][\"fuel_moisture_set_1\"])\n",
 880     "# obs_data = np.array(station['OBSERVATIONS'][\"fuel_moisture\"])\n",
 881     "# display the data retrieved\n",
 882     "#for o_time,o_data in zip (obs_time,obs_data):\n",
 883     "#    print(o_time,o_data)\n",
 884     "%matplotlib inline\n",
 885     "plt.figure(figsize=(16,4))\n",
 886     "plt.plot(obs_data,linestyle='-',c='k',label='10-h fuel data')\n",
 887     "plt.title(station['STID'] + ' 10 h fuel moisture data')\n",
 888     "plt.xlabel('Time (hours)') \n",
 889     "plt.ylabel('Fuel moisture content (%)')\n",
 890     "plt.legend()\n",
 891     " "
 892    ]
 893   },
 894   {
 895    "cell_type": "markdown",
 896    "metadata": {
 897     "id": "pY4hPeATK9wZ"
 898    },
 899    "source": [
 900     "#### 3.2.2 Acquisition of weather data"
 901    ]
 902   },
 903   {
 904    "cell_type": "markdown",
 905    "metadata": {
 906     "id": "xhyjXqxVN6B2"
 907    },
 908    "source": [
 909     "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",
 910     "\n",
 911     "First try to read the data already extracted for this RAWS and staged for download."
 912    ]
 913   },
 914   {
 915    "cell_type": "code",
 916    "execution_count": null,
 917    "metadata": {
 918     "id": "WlqJRP8Vc91o"
 919    },
 920    "outputs": [],
 921    "source": [
 922     "os.chdir('data')\n",
 923     "import json\n",
 924     "jfile = 'rtma.json'\n",
 925     "try:\n",
 926     "    ! wget --no-clobber http://math.ucdenver.edu/~jmandel/data/math4779f21/rtma.json\n",
 927     "    j = json.load(open(jfile,'r'))\n",
 928     "    print('loaded from ',jfile)\n",
 929     "    if j['obs_lat']!=obs_lat or j['obs_lon']!=obs_lon:\n",
 930     "      print('lon lat doesnot agree, need to load original RTMA files')\n",
 931     "      read_rtma=True\n",
 932     "    else:\n",
 933     "      read_rtma=False\n",
 934     "except:\n",
 935     "    print(\"can't read\",jfile,', creating')\n",
 936     "    read_rtma=True\n",
 937     "\n",
 938     "print('')"
 939    ]
 940   },
 941   {
 942    "cell_type": "markdown",
 943    "metadata": {
 944     "id": "THI6gElyHOOc"
 945    },
 946    "source": [
 947     "Next, functions to get the files, open as grib, and interpolate to the station coordinates"
 948    ]
 949   },
 950   {
 951    "cell_type": "markdown",
 952    "metadata": {
 953     "id": "2iBNHQg5hPxB"
 954    },
 955    "source": [
 956     "####<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>"
 957    ]
 958   },
 959   {
 960    "cell_type": "code",
 961    "execution_count": null,
 962    "metadata": {
 963     "id": "mxZABVDxt0gd"
 964    },
 965    "outputs": [],
 966    "source": [
 967     "# Set up environment to read RTMA gribs\n",
 968     "# we will need current numpy for pygrib - needed on Colab, tensorflow is using numpy 1.19\\\n",
 969     "if read_rtma:\n",
 970     "  import subprocess,os\n",
 971     "  def load_rtma(path,file,reload=0):\n",
 972     "    url='http://math.ucdenver.edu/~jmandel/rtma/' + path \n",
 973     "    if os.path.exists(file):\n",
 974     "      if reload:\n",
 975     "        print(file + ' already exists, removing')\n",
 976     "        os.remove(file)\n",
 977     "      else:\n",
 978     "        print(file + ' already exists, exiting')\n",
 979     "        # add checking size here\n",
 980     "        return 0\n",
 981     "    try:\n",
 982     "      ret = subprocess.check_output(['wget','--no-clobber','--output-document='+ file, url,],stderr=subprocess.STDOUT).decode() # execute command from python strings\n",
 983     "      if os.path.exists(file):\n",
 984     "        print('loaded ' + url + ' as ' + file)\n",
 985     "        return 0\n",
 986     "      else: \n",
 987     "        print('file transfer completed, but the file is missing? ' + url)  \n",
 988     "      return 1\n",
 989     "    except:\n",
 990     "      print('file transfer failed: ' + url)\n",
 991     "      return 2\n"
 992    ]
 993   },
 994   {
 995    "cell_type": "markdown",
 996    "metadata": {
 997     "id": "dQ-uJI2sy6I3"
 998    },
 999    "source": [
1000     "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."
1001    ]
1002   },
1003   {
1004    "cell_type": "code",
1005    "execution_count": null,
1006    "metadata": {
1007     "id": "PL3gxK67AlBI"
1008    },
1009    "outputs": [],
1010    "source": [
1011     "if read_rtma:\n",
1012     "  def rtma_grib(t,var):\n",
1013     "    tpath = '%4i%02i%02i/%02i' % (t.year, t.month, t.day, t.hour)  # remote path on server\n",
1014     "    tstr  = '%4i%02i%02i%02i_' % (t.year, t.month, t.day, t.hour)  # time string for local path\n",
1015     "    gribfile = os.path.join('data',tstr + var + '.grib')\n",
1016     "    remote = tpath + '/' + var + '.grib'\n",
1017     "    if load_rtma(remote,gribfile):\n",
1018     "        print('cannot load remote file',remote,'as',gribfile)\n",
1019     "        return []\n",
1020     "    else:\n",
1021     "        try:\n",
1022     "            gf=GribFile(gribfile)\n",
1023     "            v = np.array(gf[1].values())\n",
1024     "        except:\n",
1025     "            print('cannot read grib file',gribfile)\n",
1026     "            return []\n",
1027     "        print('loaded ',gribfile,' containing array shape ',v.shape)\n",
1028     "        return gf[1]   # grib message\n"
1029    ]
1030   },
1031   {
1032    "cell_type": "code",
1033    "execution_count": null,
1034    "metadata": {
1035     "id": "OY1oTYKlfd17"
1036    },
1037    "outputs": [],
1038    "source": [
1039     "if read_rtma:\n",
1040     "    times = pd.date_range(start=time_start,end=time_end,freq='1H')\n",
1041     "    varnames=['temp','td','precipa']\n",
1042     "    j =    read_interp_rtma(varnames,times,obs_lat,obs_lon)      # temperature\n",
1043     "    for varname in varnames:\n",
1044     "        j[varname]=j[varname].tolist() \n",
1045     "    j['obs_lat']=obs_lat\n",
1046     "    j['obs_lon']=obs_lon\n",
1047     "    json.dump(j,open('rtma.json','w'),indent=4)\n",
1048     "    print('done')"
1049    ]
1050   },
1051   {
1052    "cell_type": "code",
1053    "execution_count": null,
1054    "metadata": {
1055     "id": "ccp10kurAlBI"
1056    },
1057    "outputs": [],
1058    "source": [
1059     "from scipy.interpolate import LinearNDInterpolator, interpn\n",
1060     "from scipy.optimize import root\n",
1061     "def interp_to_lat_lon_slow(lats,lons,v,lat,lon): \n",
1062     "    # on mesh with coordinates lats and lons interpolate v to given lat lon\n",
1063     "    interp=LinearNDInterpolator(list(zip(lats.flatten(),lons.flatten())),v.flatten())\n",
1064     "    return interp(lat,lon)\n",
1065     "def interp_to_lat_lon(lats,lons,v,lat,lon):\n",
1066     "    # on mesh with coordinates lats and lons interpolate v to given lat lon\n",
1067     "    points=(np.array(range(lats.shape[0]),float),np.array(range(lats.shape[1]),float))  # uniform mesh\n",
1068     "    def res(ij):  # interpolation of lons lats on the uniform mesh, to noninteger coordinates   \n",
1069     "       return np.hstack((interpn(points,lats,ij)-lat, interpn(points,lons,ij)-lon))\n",
1070     "    # solve for xi,xj such that lats(xi,xj)=lat lons(xi,xj)=lon, then interpolate to (xi, xj) on uniform grid \n",
1071     "    result = root(res,(0,0)) # solve res(ij) = 0\n",
1072     "    if not result.success:\n",
1073     "        print(result.message)\n",
1074     "        exit(1)\n",
1075     "    return interpn(points,v,result.x) \n"
1076    ]
1077   },
1078   {
1079    "cell_type": "markdown",
1080    "metadata": {
1081     "id": "jvnpq6S5AlBI"
1082    },
1083    "source": [
1084     "The interpolation function needs to  be tested."
1085    ]
1086   },
1087   {
1088    "cell_type": "code",
1089    "execution_count": null,
1090    "metadata": {
1091     "id": "NVMJBYI7AlBI"
1092    },
1093    "outputs": [],
1094    "source": [
1095     "def interp_to_lat_lon_test(lats,lons):\n",
1096     "    print('testing interp_to_lat_lon')\n",
1097     "    vx, vy = np.meshgrid(range(lats.shape[0]),range(lats.shape[1]),indexing='ij')\n",
1098     "    i, j = (1,2)\n",
1099     "    lat,lon = ((lats[i,j]+lats[i+1,j+1])/2,(lons[i,j]+lons[i+1,j+1])/2)\n",
1100     "    vi = interp_to_lat_lon(lats,lons,vx,lat,lon)\n",
1101     "    vj = interp_to_lat_lon(lats,lons,vy,lat,lon)\n",
1102     "    print(vi,vj,'should be about',i+0.5,j+0.5)\n",
1103     "    test_slow = 0\n",
1104     "    if test_slow:\n",
1105     "        print('Testing against the standard slow method scipy.interpolate.LinearNDInterpolator. Please wait...')\n",
1106     "        vi_slow = interp_to_lat_lon_slow(lats,lons,vx,lat,lon)\n",
1107     "        print(vi_slow)\n",
1108     "        vj_slow = interp_to_lat_lon_slow(lats,lons,vy,lat,lon)\n",
1109     "        print(vj_slow)\n",
1110     "        \n",
1111     "#gf = rtma_grib(start_time,'temp')      #  read the first grib file and use it to test interpolation\n",
1112     "#lats, lons = gf.latlons()\n",
1113     "#interp_to_lat_lon_test(lats,lons)\n"
1114    ]
1115   },
1116   {
1117    "cell_type": "code",
1118    "execution_count": null,
1119    "metadata": {
1120     "id": "vt-Mk8fIc91m"
1121    },
1122    "outputs": [],
1123    "source": [
1124     "#%debug\n"
1125    ]
1126   },
1127   {
1128    "cell_type": "markdown",
1129    "metadata": {
1130     "id": "LQbWB_3GAlBI"
1131    },
1132    "source": [
1133     "Now we are ready for a function to read the RTMA files and interpolate to the station coordinates"
1134    ]
1135   },
1136   {
1137    "cell_type": "code",
1138    "execution_count": null,
1139    "metadata": {
1140     "id": "b3JJH3XPAlBI"
1141    },
1142    "outputs": [],
1143    "source": [
1144     "if read_rtma:\n",
1145     "  import pandas as pd, json\n",
1146     "  def read_interp_rtma(varnames,times,lat,lon):\n",
1147     "    # read RTMA from start_time to end_time and interpolate to obs_lat obs_lon\n",
1148     "    ntimes = len(times)\n",
1149     "    time_str = 'time_str'\n",
1150     "    j={time_str:times.strftime('%Y-%m-%d %H:%M').tolist()}\n",
1151     "    for varname in varnames:\n",
1152     "        j[varname]=np.full(ntimes,np.nan)  # initialize array of nans as list\n",
1153     "    n=0\n",
1154     "    for t in times:\n",
1155     "        tim=t.strftime('%Y-%m-%d %H:%M')\n",
1156     "        should_be = j[time_str][n]\n",
1157     "        if tim != should_be:\n",
1158     "            print('n=',n,'time',tim,'expected',should_be)\n",
1159     "            raise 'Invalid time' \n",
1160     "        for varname in varnames:\n",
1161     "            gf = rtma_grib(t,varname)   # read and create grib object, download if needed\n",
1162     "            if gf:\n",
1163     "                lats,lons = gf.latlons()    # coordinates\n",
1164     "                v = gf.values()\n",
1165     "                vi=interp_to_lat_lon(lats,lons,v,lat,lon) # append to array\n",
1166     "                print(varname,'at',t,'interpolated to',lat,lon,' value ',vi)\n",
1167     "                j[varname][n] = vi\n",
1168     "            else:\n",
1169     "                print(varname,'at',t,' could not be loaded')\n",
1170     "        n = n+1\n",
1171     "    return j"
1172    ]
1173   },
1174   {
1175    "cell_type": "code",
1176    "execution_count": null,
1177    "metadata": {
1178     "id": "bMpYIZT6c91o"
1179    },
1180    "outputs": [],
1181    "source": [
1182     "# %debug\n"
1183    ]
1184   },
1185   {
1186    "cell_type": "markdown",
1187    "metadata": {
1188     "id": "KVXBjGA0CiXr"
1189    },
1190    "source": [
1191     "#### 3.2.3 Preprocessing and visualization of the weather data"
1192    ]
1193   },
1194   {
1195    "cell_type": "code",
1196    "execution_count": null,
1197    "metadata": {
1198     "id": "fNA3Vbo1c91o",
1199     "jupyter": {
1200      "source_hidden": true
1201     },
1202     "tags": []
1203    },
1204    "outputs": [],
1205    "source": [
1206     "rtma = j\n",
1207     "td = np.array(rtma['td'])\n",
1208     "t2 = np.array(rtma['temp'])\n",
1209     "rain=np.array(rtma['precipa'])\n",
1210     "# compute relative humidity\n",
1211     "rh = 100*np.exp(17.625*243.04*(td - t2) / (243.04 + t2 - 273.15) / (243.0 + td - 273.15))\n",
1212     "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",
1213     "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))"
1214    ]
1215   },
1216   {
1217    "cell_type": "code",
1218    "execution_count": null,
1219    "metadata": {
1220     "id": "tZIK59bJAlBJ",
1221     "jupyter": {
1222      "source_hidden": true
1223     },
1224     "tags": []
1225    },
1226    "outputs": [],
1227    "source": [
1228     "%matplotlib inline\n",
1229     "plt.figure(figsize=(16,4))\n",
1230     "plt.plot(t2,linestyle='-',c='k',label='Temperature')\n",
1231     "plt.title(station['STID'] + ' Temperature')\n",
1232     "plt.xlabel('Time (hours)') \n",
1233     "plt.ylabel('Temperature (K)')\n",
1234     "plt.legend()"
1235    ]
1236   },
1237   {
1238    "cell_type": "code",
1239    "execution_count": null,
1240    "metadata": {
1241     "id": "LbyqcuXYAlBJ"
1242    },
1243    "outputs": [],
1244    "source": [
1245     "%matplotlib inline\n",
1246     "plt.figure(figsize=(16,4))\n",
1247     "plt.plot(td,linestyle='-',c='k',label='Dew point')\n",
1248     "plt.title(station['STID'] + ' Dew point (K)')\n",
1249     "plt.xlabel('Time (hours)') \n",
1250     "plt.ylabel('Dew point (K)')\n",
1251     "plt.legend()"
1252    ]
1253   },
1254   {
1255    "cell_type": "code",
1256    "execution_count": null,
1257    "metadata": {
1258     "id": "dfoOK2kSc91p"
1259    },
1260    "outputs": [],
1261    "source": [
1262     "%matplotlib inline\n",
1263     "plt.figure(figsize=(16,4))\n",
1264     "plt.plot(rh,linestyle='-',c='k',label='Dew point')\n",
1265     "plt.title(station['STID'] + ' relative humidity')\n",
1266     "plt.xlabel('Time (hours)') \n",
1267     "plt.ylabel('Relative humidity (%)')\n",
1268     "plt.legend()"
1269    ]
1270   },
1271   {
1272    "cell_type": "code",
1273    "execution_count": null,
1274    "metadata": {
1275     "id": "MWTJ5b2kc91p",
1276     "jupyter": {
1277      "source_hidden": true
1278     },
1279     "tags": []
1280    },
1281    "outputs": [],
1282    "source": [
1283     "%matplotlib inline\n",
1284     "plt.figure(figsize=(16,4))\n",
1285     "plt.plot(Ed,linestyle='-',c='r',label='drying equilibrium')\n",
1286     "plt.plot(Ew,linestyle=':',c='b',label='wetting equilibrium')\n",
1287     "plt.title(station['STID'] + ' drying and wetting equilibria')\n",
1288     "plt.xlabel('Time (hours)') \n",
1289     "plt.ylabel('Fuel moisture contents (%)')\n",
1290     "plt.legend()"
1291    ]
1292   },
1293   {
1294    "cell_type": "markdown",
1295    "metadata": {
1296     "id": "jY3_eeBRc91p"
1297    },
1298    "source": [
1299     " "
1300    ]
1301   },
1302   {
1303    "cell_type": "code",
1304    "execution_count": null,
1305    "metadata": {
1306     "id": "PQKSRvRSAlBJ"
1307    },
1308    "outputs": [],
1309    "source": [
1310     "%matplotlib inline\n",
1311     "plt.figure(figsize=(16,4))\n",
1312     "plt.plot(rain,linestyle='-',c='k',label='Precipitation')\n",
1313     "plt.title(station['STID'] + ' Precipitation' )\n",
1314     "plt.xlabel('Time (hours)') \n",
1315     "plt.ylabel('Precipitation (mm/hour)')\n",
1316     "plt.legend()"
1317    ]
1318   },
1319   {
1320    "cell_type": "code",
1321    "execution_count": null,
1322    "metadata": {
1323     "id": "Dwbt4UXfro5x"
1324    },
1325    "outputs": [],
1326    "source": [
1327     "print(rain[1900:2000])"
1328    ]
1329   },
1330   {
1331    "cell_type": "markdown",
1332    "metadata": {
1333     "id": "_yRu_7WvHc6P"
1334    },
1335    "source": [
1336     "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."
1337    ]
1338   },
1339   {
1340    "cell_type": "code",
1341    "execution_count": null,
1342    "metadata": {
1343     "id": "XPYO_Iuvc91q"
1344    },
1345    "outputs": [],
1346    "source": [
1347     "rain[rain > 1000] = np.NaN"
1348    ]
1349   },
1350   {
1351    "cell_type": "code",
1352    "execution_count": null,
1353    "metadata": {
1354     "id": "GYWTfbBBc91q"
1355    },
1356    "outputs": [],
1357    "source": [
1358     "%matplotlib inline\n",
1359     "plt.figure(figsize=(16,4))\n",
1360     "plt.plot(rain,linestyle='-',c='k',label='Precipitation')\n",
1361     "plt.title(station['STID'] + ' Precipitation' )\n",
1362     "plt.xlabel('Time (hours)') \n",
1363     "plt.ylabel('Precipitation (mm/hour)')\n",
1364     "plt.legend()"
1365    ]
1366   },
1367   {
1368    "cell_type": "markdown",
1369    "metadata": {
1370     "id": "Q_L0R2Njc91q"
1371    },
1372    "source": [
1373     "Fix some missing data, then we can use the data for up to 1942 hours until a biger gap."
1374    ]
1375   },
1376   {
1377    "cell_type": "code",
1378    "execution_count": null,
1379    "metadata": {
1380     "id": "_tkU7UJic91q"
1381    },
1382    "outputs": [],
1383    "source": [
1384     "# fix isolated nans\n",
1385     "def fixnan(a,n):\n",
1386     "    for c in range(n):\n",
1387     "        for i in np.where(np.isnan(a)):\n",
1388     "            a[i]=0.5*(a[i-1]+a[i+1])\n",
1389     "        if not any(np.isnan(a)):\n",
1390     "            break\n",
1391     "    return a\n",
1392     "\n",
1393     "rain=fixnan(rain,2)\n",
1394     "t2=fixnan(t2,2)\n",
1395     "rh=fixnan(rh,2)\n",
1396     "obs_data=fixnan(obs_data,2)\n",
1397     "Ed=fixnan(Ed,2)\n",
1398     "Ew=fixnan(Ew,2)\n",
1399     "\n",
1400     "print(np.where(np.isnan(rain)))\n",
1401     "print(np.where(np.isnan(t2)))\n",
1402     "print(np.where(np.isnan(rh)))\n",
1403     "print(np.where(np.isnan(obs_data)))"
1404    ]
1405   },
1406   {
1407    "cell_type": "markdown",
1408    "metadata": {
1409     "id": "XqQYnyI9DIy1"
1410    },
1411    "source": [
1412     "## 4 Results"
1413    ]
1414   },
1415   {
1416    "cell_type": "markdown",
1417    "metadata": {
1418     "id": "2tIC_Tqnc91r"
1419    },
1420    "source": [
1421     "### 4.1 Kalman filter with fuel moisture observations, followed by forecasting\n",
1422     "We run the model first with Kalman filter for 150 hours. The observations are the RAWS data\n",
1423     "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",
1424     "In a real forecasting application, the model would be run from weather forecast rather than data."
1425    ]
1426   },
1427   {
1428    "cell_type": "code",
1429    "execution_count": null,
1430    "metadata": {
1431     "id": "aXnSQM7wc91r"
1432    },
1433    "outputs": [],
1434    "source": [
1435     "# run KF on an initial data seqment\n",
1436     "import numpy as np\n",
1437     "import matplotlib.pyplot as plt \n",
1438     "\n",
1439     "hours=1200 # total simulation\n",
1440     "h2 = 300\n",
1441     "m = np.zeros(hours) # preallocate\n",
1442     "m[0]= obs_data[0]             # initial state  \n",
1443     "P = np.zeros(hours)\n",
1444     "P[0] = 1e-3 # background state variance\n",
1445     "H = np.array([1.])   # all oQ = np.array([0.02]) # process noise variancebserved\n",
1446     "Q = np.array([1e-3]) # process noise variance\n",
1447     "R = np.array([1e-3]) # data variance\n",
1448     "for t in range(hours-1):\n",
1449     "    # using lambda construction to pass additional arguments to the model \n",
1450     "    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",
1451     "        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",
1452     "                    d=obs_data[t],H=H,R=R)\n",
1453     "    else:  # just advance to next hour, no process noise\n",
1454     "        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)"
1455    ]
1456   },
1457   {
1458    "cell_type": "code",
1459    "execution_count": null,
1460    "metadata": {
1461     "id": "peMi-OF3c91r"
1462    },
1463    "outputs": [],
1464    "source": [
1465     "%matplotlib inline\n",
1466     "plt.figure(figsize=(16,4))\n",
1467     "plt.plot(Ed[:hours],linestyle='--',c='r',label='Drying Equilibrium')\n",
1468     "plt.plot(Ew[:hours],linestyle='--',c='b',label='Wetting Equilibrium')\n",
1469     "plt.plot(obs_data[:hours],linestyle=':',c='k',label='RAWS data')\n",
1470     "plt.plot(m[:h2],linestyle='-',c='k',label='filtered')\n",
1471     "plt.plot(range(h2,hours),m[h2:hours],linestyle='-',c='r',label='forecast')\n",
1472     "plt.title(station['STID'] + ' Kalman filtering and forecast with real data')\n",
1473     "plt.xlabel('Time (hours)') \n",
1474     "plt.ylabel('Fuel moisture content (%)')\n",
1475     "plt.legend()"
1476    ]
1477   },
1478   {
1479    "cell_type": "markdown",
1480    "metadata": {
1481     "id": "3TnwXYcLc91r"
1482    },
1483    "source": [
1484     "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. "
1485    ]
1486   },
1487   {
1488    "cell_type": "markdown",
1489    "metadata": {
1490     "id": "8SuVNg8TsW4d"
1491    },
1492    "source": [
1493     "### 4.3 Kalman filter on the augmented model"
1494    ]
1495   },
1496   {
1497    "cell_type": "markdown",
1498    "metadata": {
1499     "id": "FYAbWNCfk2wD"
1500    },
1501    "source": [
1502     "Run augmented filter and plot the result:\n"
1503    ]
1504   },
1505   {
1506    "cell_type": "code",
1507    "execution_count": null,
1508    "metadata": {
1509     "id": "Q3NHr3oBsDg6"
1510    },
1511    "outputs": [],
1512    "source": [
1513     "m,Ec = run_augmented_kf(obs_data,Ed,Ew,rain,h2,hours)  # extract from state"
1514    ]
1515   },
1516   {
1517    "cell_type": "code",
1518    "execution_count": null,
1519    "metadata": {
1520     "id": "hlkby_oTlB_f"
1521    },
1522    "outputs": [],
1523    "source": [
1524     "title = station['STID'] +' Kalman filtering and forecast with augmented state, real data. Training 0:%i hmax' % h2\n",
1525     "def plot_moisture(hmin,hmax):\n",
1526     "  print('training from 0 to',h2,'plot from',hmin,'to',hmax)\n",
1527     "  plt.figure(figsize=(16,4))\n",
1528     "  plt.plot(range(hmin,hmax),Ed[hmin:hmax],linestyle='--',c='r',label='Drying Equilibrium (%)')\n",
1529     "  plt.plot(range(hmin,hmax),Ew[hmin:hmax],linestyle='--',c='b',label='Wetting Equilibrium (%)')\n",
1530     "  plt.plot(range(hmin,hmax),Ec[hmin:hmax],linestyle='--',c='g',label='Equilibrium Correction (%)')\n",
1531     "  plt.plot(range(hmin,hmax),m[hmin:hmax],linestyle='-',c='k',label='filtered')\n",
1532     "  plt.plot(range(hmin,hmax),obs_data[hmin:hmax],linestyle='-',c='b',label='RAWS data (%)')\n",
1533     "  plt.plot(range(hmin,hmax),rain[hmin:hmax],linestyle='-',c='b',label='RTMA rain (mm/h)')\n",
1534     "  plt.title(title)\n",
1535     "  if hmin>=h2:\n",
1536     "    plt.plot(m[hmin:h2],linestyle='-',c='k',label='Filtered')\n",
1537     "  h1 = np.maximum(hmin,h2)\n",
1538     "  plt.plot(range(h1,hmax),m[h1:hmax],linestyle='-',c='r',label='Forecast (%)')\n",
1539     "  plt.xlabel('Time (hours)') \n",
1540     "  plt.ylabel('Fuel moisture content (%)')\n",
1541     "  plt.legend()"
1542    ]
1543   },
1544   {
1545    "cell_type": "code",
1546    "execution_count": null,
1547    "metadata": {},
1548    "outputs": [],
1549    "source": [
1550     "os.chdir('..')\n",
1551     "from data_funcs import to_json, from_json"
1552    ]
1553   },
1554   {
1555    "cell_type": "code",
1556    "execution_count": null,
1557    "metadata": {},
1558    "outputs": [],
1559    "source": [
1560     "os.chdir('data')\n",
1561     "kf_orig={'title':title,'hours':hours,'h2':h2,'Ed':Ed,'Ew':Ew,'Ec':Ec,'rain':rain,\n",
1562     "            'fm':obs_data,'m':m,'note':'RAWS and RTMA data + m from augmented KF in fmda_kf_rnn_orig'}\n",
1563     "to_json(kf_orig,'kf_orig.json')"
1564    ]
1565   },
1566   {
1567    "cell_type": "code",
1568    "execution_count": null,
1569    "metadata": {
1570     "id": "q-h5omKgnow2"
1571    },
1572    "outputs": [],
1573    "source": [
1574     "plot_moisture(0,hours)"
1575    ]
1576   },
1577   {
1578    "cell_type": "markdown",
1579    "metadata": {
1580     "id": "0w0YtHtqnza5"
1581    },
1582    "source": [
1583     "A detailed view of transition from training to forecast:"
1584    ]
1585   },
1586   {
1587    "cell_type": "code",
1588    "execution_count": null,
1589    "metadata": {
1590     "id": "B7sXGUotc91s"
1591    },
1592    "outputs": [],
1593    "source": [
1594     "plot_moisture(0,600)\n",
1595     "\n"
1596    ]
1597   },
1598   {
1599    "cell_type": "code",
1600    "execution_count": null,
1601    "metadata": {
1602     "id": "xy7sIs0z_Kk6"
1603    },
1604    "outputs": [],
1605    "source": [
1606     "plot_moisture(300,500)"
1607    ]
1608   },
1609   {
1610    "cell_type": "code",
1611    "execution_count": null,
1612    "metadata": {
1613     "id": "y-C6IRFVxGUR"
1614    },
1615    "outputs": [],
1616    "source": [
1617     "plot_moisture(300,800)"
1618    ]
1619   },
1620   {
1621    "cell_type": "code",
1622    "execution_count": null,
1623    "metadata": {
1624     "id": "TvlCtT0X2ejp"
1625    },
1626    "outputs": [],
1627    "source": [
1628     "plot_moisture(800,1200)"
1629    ]
1630   },
1631   {
1632    "cell_type": "markdown",
1633    "metadata": {
1634     "id": "7W03QTo3c91t"
1635    },
1636    "source": [
1637     "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",
1638     "\n",
1639     "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? "
1640    ]
1641   },
1642   {
1643    "cell_type": "markdown",
1644    "metadata": {
1645     "id": "owEI4EtTo7Ek"
1646    },
1647    "source": [
1648     "A detailed view of rain episode:"
1649    ]
1650   },
1651   {
1652    "cell_type": "code",
1653    "execution_count": null,
1654    "metadata": {
1655     "id": "C_hoDjgtpMEJ"
1656    },
1657    "outputs": [],
1658    "source": [
1659     "plot_moisture(900,1100)"
1660    ]
1661   },
1662   {
1663    "cell_type": "markdown",
1664    "metadata": {
1665     "id": "DRraWhwdpSkV"
1666    },
1667    "source": [
1668     "It seems there is some rain that the model does not know about."
1669    ]
1670   },
1671   {
1672    "cell_type": "markdown",
1673    "metadata": {
1674     "id": "1STfnlT40rPX"
1675    },
1676    "source": [
1677     "## RNN for real data, no rain yet"
1678    ]
1679   },
1680   {
1681    "cell_type": "markdown",
1682    "metadata": {
1683     "id": "3cwY43iSnQ0t"
1684    },
1685    "source": [
1686     "#### Linear modeling by RELU - potential for generalization"
1687    ]
1688   },
1689   {
1690    "cell_type": "code",
1691    "execution_count": null,
1692    "metadata": {
1693     "id": "MotzNBvOnFvC"
1694    },
1695    "outputs": [],
1696    "source": [
1697     "def RELU(x):\n",
1698     "  if x>0. :\n",
1699     "    return x\n",
1700     "  else:\n",
1701     "    return 0.\n",
1702     "\n",
1703     "# network computing z = a*x1 + b*x2 with offset c\n",
1704     "def linrelu(x,a,b,c):\n",
1705     "  y = np.dot(np.array([[a, b], [-a, -b] ]), x) + np.array([c, -c])\n",
1706     "  y[0]=RELU(y[0])\n",
1707     "  y[1]=RELU(y[1])\n",
1708     "  return(np.dot([1,-1],y))-c\n",
1709     "x = [1,2]\n",
1710     "a = 2\n",
1711     "b = 4\n",
1712     "c = 3\n",
1713     "print(a*x[0]+b*x[1])\n",
1714     "linrelu(x,a,b,c)"
1715    ]
1716   },
1717   {
1718    "cell_type": "markdown",
1719    "metadata": {
1720     "id": "-p6dcLua_udD"
1721    },
1722    "source": [
1723     "### Basic RNN on real data "
1724    ]
1725   },
1726   {
1727    "cell_type": "markdown",
1728    "metadata": {
1729     "id": "gSmbDPZIHbTr"
1730    },
1731    "source": [
1732     "Try with E average between drying and wetting"
1733    ]
1734   },
1735   {
1736    "cell_type": "code",
1737    "execution_count": null,
1738    "metadata": {
1739     "id": "ymhNMZkoHfCl"
1740    },
1741    "outputs": [],
1742    "source": [
1743     "E = (Ed + Ew)/2\n",
1744     "print(Ed.shape,Ew.shape,rain.shape)\n",
1745     "first_rain=np.nonzero(rain>0)[0][0]\n",
1746     "print(first_rain)\n",
1747     "hours=first_rain\n",
1748     "E=E[:hours]\n",
1749     "data=obs_data[:hours]\n",
1750     "scale=False\n",
1751     "\n",
1752     "# transform as 2D, (timesteps, features) and (timesteps, outputs)\n",
1753     "Et = np.reshape(E,[E.shape[0],1])\n",
1754     "datat = np.reshape(data,[data.shape[0],1])\n",
1755     "if scale:\n",
1756     "  scalerx = MinMaxScaler()\n",
1757     "  scalerx.fit(Et)\n",
1758     "  Et = scalerx.transform(Et)\n",
1759     "  scalery = MinMaxScaler()\n",
1760     "  scalery.fit(datat)\n",
1761     "  datat = scalery.transform(datat)"
1762    ]
1763   },
1764   {
1765    "cell_type": "markdown",
1766    "metadata": {
1767     "id": "DPcxv85XILdn"
1768    },
1769    "source": [
1770     "Create the model again"
1771    ]
1772   },
1773   {
1774    "cell_type": "code",
1775    "execution_count": null,
1776    "metadata": {
1777     "id": "gEkbHZSqIOq1"
1778    },
1779    "outputs": [],
1780    "source": [
1781     "# split data\n",
1782     "return_sequences=False\n",
1783     "x_train, y_train = staircase(Et,datat,timesteps=5,trainsteps=h2,\n",
1784     "                             return_sequences=return_sequences)\n",
1785     "print('x_train shape=',x_train.shape)\n",
1786     "samples, timesteps, features = x_train.shape\n",
1787     "print('y_train shape=',y_train.shape)\n",
1788     "# the simplest model possible\n",
1789     "activation=['linear','linear']\n",
1790     "hidden_units=3\n",
1791     "dense_units=1\n",
1792     "dense_layers=1\n",
1793     "features=1\n",
1794     "hours=Et.shape[0]\n",
1795     "h0 = tf.convert_to_tensor(datat[:samples],dtype=tf.float32)\n",
1796     "# print('initial state=',h0)\n",
1797     "# statefull model version for traning\n",
1798     "import moisture_rnn\n",
1799     "model_fit=moisture_rnn.create_RNN_2(hidden_units=hidden_units, \n",
1800     "                        dense_units=dense_units, \n",
1801     "                        batch_shape=(samples,timesteps,features),\n",
1802     "                        stateful=True,\n",
1803     "                        return_sequences=return_sequences,\n",
1804     "                        # initial_state=h0,\n",
1805     "                        activation=activation,\n",
1806     "                        dense_layers=dense_layers)\n",
1807     "# same model stateless for prediction on the entire dataset - to start onlg\n",
1808     "# the real application will switch to prediction after training data end\n",
1809     "# and start from the state there\n",
1810     "print('model_fit input shape',x_train.shape,'output shape',model_fit(x_train).shape)\n",
1811     "from keras.utils.vis_utils import plot_model\n",
1812     "plot_model(model_fit, to_file='model_plot.png', \n",
1813     "           show_shapes=True, show_layer_names=True)"
1814    ]
1815   },
1816   {
1817    "cell_type": "code",
1818    "execution_count": null,
1819    "metadata": {
1820     "id": "jtFJQu33NqfL"
1821    },
1822    "outputs": [],
1823    "source": [
1824     "model_predict=create_RNN_2(hidden_units=hidden_units, dense_units=dense_units,  \n",
1825     "                        input_shape=(hours,features),stateful = False,\n",
1826     "                        return_sequences=True,\n",
1827     "                        activation=activation,dense_layers=dense_layers)\n",
1828     "# model_predict=create_RNN_sequences(hidden_units=1, dense_units=1, input_shape=(hours,1), \n",
1829     "#                        activation=['linear', 'linear'])\n",
1830     "print('model_predict input shape',Et.shape,'output shape',model_predict(Et).shape)\n",
1831     "print(model_predict.summary())\n",
1832     "from keras.utils.vis_utils import plot_model\n",
1833     "plot_model(model_predict, to_file='model_plot.png', \n",
1834     "           show_shapes=True, show_layer_names=True)"
1835    ]
1836   },
1837   {
1838    "cell_type": "code",
1839    "execution_count": null,
1840    "metadata": {
1841     "id": "wuxh5pq0OMSa"
1842    },
1843    "outputs": [],
1844    "source": [
1845     "# fitting\n",
1846     "DeltaE = 0\n",
1847     "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",
1848     "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",
1849     "w=model_fit.get_weights()\n",
1850     "for i in range(len(w)):\n",
1851     "  print('weight',i,'shape',w[i].shape,'ndim',w[i].ndim,'given',w_initial[i].shape)\n",
1852     "  for j in range(w[i].shape[0]):\n",
1853     "    if w[i].ndim==2:\n",
1854     "      for k in range(w[i].shape[1]):\n",
1855     "        w[i][j][k]=w_initial[i][0][0]/w[i].shape[0]\n",
1856     "    else:\n",
1857     "      w[i][j]=w_initial[i][0]\n",
1858     "model_fit.set_weights(w)\n",
1859     "model_fit.fit(x_train, y_train, epochs=5000, verbose=0,batch_size=samples)\n",
1860     "w_fitted=model_fit.get_weights()\n",
1861     "for i in range(len(w)):\n",
1862     "  print('weight',i,' exact:',w_exact[i],':  initial:',w_initial[i],' fitted:',w_fitted[i])"
1863    ]
1864   },
1865   {
1866    "cell_type": "code",
1867    "execution_count": null,
1868    "metadata": {
1869     "tags": []
1870    },
1871    "outputs": [],
1872    "source": []
1873   },
1874   {
1875    "cell_type": "code",
1876    "execution_count": null,
1877    "metadata": {
1878     "id": "uJz1EgPyRTEH"
1879    },
1880    "outputs": [],
1881    "source": [
1882     "# evaluate model\n",
1883     "model_predict.set_weights(w_fitted)\n",
1884     "x_input=np.reshape(Et,(1, hours, 1))\n",
1885     "y_output = model_predict.predict(x_input)\n",
1886     "print('x_input.shape=',x_input.shape,'y_output.shape=',y_output.shape)\n",
1887     "print(shift)\n",
1888     "m = np.reshape(y_output,hours)\n",
1889     "print('weights=',w)\n",
1890     "if scale:\n",
1891     "    print('scaling')\n",
1892     "    m = scalery.inverse_transform(m)\n",
1893     "m = np.reshape(m,hours)\n",
1894     "hour=np.array(range(hours))\n",
1895     "title=\"First RNN forecast\"\n",
1896     "plt.figure(figsize=(16,4))\n",
1897     "plt.plot(hour,E,linestyle='--',c='r',label='E=Equilibrium data')\n",
1898     "# print(len(hour),len(m_f))\n",
1899     "plt.scatter(hour,data,c='b',label='data=10-h fuel data')\n",
1900     "if m is not None:\n",
1901     "    plt.plot(hour[:h2],m[:h2],linestyle='-',c='k',label='m=filtered')\n",
1902     "    plt.plot(hour[h2:hours],m[h2:hours],linestyle='-',c='r',label='m=forecast')\n",
1903     "plt.title(title) \n",
1904     "plt.legend()\n"
1905    ]
1906   },
1907   {
1908    "cell_type": "code",
1909    "execution_count": null,
1910    "metadata": {},
1911    "outputs": [],
1912    "source": []
1913   },
1914   {
1915    "cell_type": "code",
1916    "execution_count": null,
1917    "metadata": {
1918     "id": "VSwtgKPJPnH4"
1919    },
1920    "outputs": [],
1921    "source": [
1922     "# plot subinterval only\n",
1923     "def plot_int(lb=0,ub=hours,title=\"RNN forecast\"):\n",
1924     "  hour=np.array(range(hours))\n",
1925     "  plt.figure(figsize=(16,4))\n",
1926     "  plt.plot(hour[lb:ub],E[lb:ub],linestyle='--',c='r',label='Equilibrium data')\n",
1927     "  # plt.scatter(hour[lb:ub],data[lb:ub],c='b',label='data=10-h fuel data')\n",
1928     "  plt.plot(hour[lb:ub],m[lb:ub],linestyle='-',c='b',label='data=10-h fuel data')\n",
1929     "  if lb <= h2:\n",
1930     "    ub1 = min(h2,ub)\n",
1931     "    plt.plot(hour[lb:ub1],m[lb:ub1],linestyle='-',c='k',label='filtered')\n",
1932     "  if ub >= h2:\n",
1933     "    lb1 = max(h2,lb)\n",
1934     "    plt.plot(hour[lb1:ub],m[lb1:ub],linestyle='-',c='r',label='forecast')\n",
1935     "  plt.title(title) \n",
1936     "  plt.legend()"
1937    ]
1938   },
1939   {
1940    "cell_type": "code",
1941    "execution_count": null,
1942    "metadata": {
1943     "id": "vCjk9hZtkFym"
1944    },
1945    "outputs": [],
1946    "source": [
1947     "plot_int()"
1948    ]
1949   },
1950   {
1951    "cell_type": "code",
1952    "execution_count": null,
1953    "metadata": {
1954     "id": "Sd3fDOnvmmdp"
1955    },
1956    "outputs": [],
1957    "source": [
1958     "plot_int(0,300)"
1959    ]
1960   },
1961   {
1962    "cell_type": "code",
1963    "execution_count": null,
1964    "metadata": {
1965     "id": "vHkc4KHdkAJp"
1966    },
1967    "outputs": [],
1968    "source": [
1969     "plot_int(300,500)"
1970    ]
1971   },
1972   {
1973    "cell_type": "code",
1974    "execution_count": null,
1975    "metadata": {
1976     "id": "Km5VWhcJlyvV"
1977    },
1978    "outputs": [],
1979    "source": [
1980     "plot_int(500,800)"
1981    ]
1982   },
1983   {
1984    "cell_type": "markdown",
1985    "metadata": {
1986     "id": "TBayRudFcZWP"
1987    },
1988    "source": [
1989     "Next step: two features - drying and wetting equilibria"
1990    ]
1991   },
1992   {
1993    "cell_type": "code",
1994    "execution_count": null,
1995    "metadata": {
1996     "id": "SGbgxOm_kEc4"
1997    },
1998    "outputs": [],
1999    "source": [
2000     "print(Ed.shape,Ew.shape,rain.shape)\n",
2001     "first_rain=np.nonzero(rain>0)[0][0]\n",
2002     "print(first_rain)\n",
2003     "hours=first_rain\n",
2004     "Ed=Ed[:hours]\n",
2005     "Ew=Ew[:hours]\n",
2006     "h2 = 300\n",
2007     "# print(Ed.shape,Ew.shape)\n",
2008     "# (timesteps, features)\n",
2009     "Et = np.vstack((Ed, Ew)).T\n",
2010     "print(E.shape)\n",
2011     "data=obs_data[:hours]\n",
2012     "\n",
2013     "scale=False\n",
2014     "\n",
2015     "# transform as 2D, (timesteps, features) and (timesteps, outputs)\n",
2016     "datat = np.reshape(data,[data.shape[0],1])\n",
2017     "if scale:\n",
2018     "  scalerx = MinMaxScaler()\n",
2019     "  scalerx.fit(Et)\n",
2020     "  Et = scalerx.transform(Et)\n",
2021     "  scalery = MinMaxScaler()\n",
2022     "  scalery.fit(datat)\n",
2023     "  datat = scalery.transform(datat)"
2024    ]
2025   },
2026   {
2027    "cell_type": "code",
2028    "execution_count": null,
2029    "metadata": {},
2030    "outputs": [],
2031    "source": [
2032     "from utils import hash2"
2033    ]
2034   },
2035   {
2036    "cell_type": "code",
2037    "execution_count": null,
2038    "metadata": {},
2039    "outputs": [],
2040    "source": [
2041     "# Set seed for reproducibility\n",
2042     "reproducibility.set_seed()"
2043    ]
2044   },
2045   {
2046    "cell_type": "code",
2047    "execution_count": null,
2048    "metadata": {
2049     "id": "b6aJAvBEkEBl"
2050    },
2051    "outputs": [],
2052    "source": [
2053     "# split data\n",
2054     "return_sequences=False\n",
2055     "x_train, y_train = staircase(Et,datat,timesteps=5,trainsteps=h2,\n",
2056     "                             return_sequences=return_sequences)\n",
2057     "print('x_train shape=',x_train.shape)\n",
2058     "samples, timesteps, features = x_train.shape\n",
2059     "print('y_train shape=',y_train.shape)\n",
2060     "# the simplest model possible\n",
2061     "activation=['linear','linear']\n",
2062     "hidden_units=6\n",
2063     "dense_units=1\n",
2064     "dense_layers=1\n",
2065     "features=Et.shape[1]\n",
2066     "hours=Et.shape[0]\n",
2067     "h0 = tf.convert_to_tensor(datat[:samples],dtype=tf.float32)\n",
2068     "# print('initial state=',h0)\n",
2069     "# statefull model version for traning\n",
2070     "import moisture_rnn\n",
2071     "model_fit=moisture_rnn.create_RNN_2(hidden_units=hidden_units, \n",
2072     "                        dense_units=dense_units, \n",
2073     "                        batch_shape=(samples,timesteps,features),\n",
2074     "                        stateful=True,\n",
2075     "                        return_sequences=return_sequences,\n",
2076     "                        # initial_state=h0,\n",
2077     "                        activation=activation,\n",
2078     "                        dense_layers=dense_layers)\n",
2079     "# same model stateless for prediction on the entire dataset - to start onlg\n",
2080     "# the real application will switch to prediction after training data end\n",
2081     "# and start from the state there\n",
2082     "print('model_fit input shape',x_train.shape,'output shape',model_fit(x_train).shape)\n",
2083     "print('model_fit input shape',x_train.shape,'output shape',y_train.shape)\n",
2084     "from keras.utils.vis_utils import plot_model\n",
2085     "plot_model(model_fit, to_file='model_plot.png', \n",
2086     "           show_shapes=True, show_layer_names=True)"
2087    ]
2088   },
2089   {
2090    "cell_type": "code",
2091    "execution_count": null,
2092    "metadata": {},
2093    "outputs": [],
2094    "source": [
2095     "## Check 1: equilibrium input data the same\n",
2096     "\n",
2097     "print(hash2(Et))\n",
2098     "print(hash2(x_train))\n",
2099     "print(hash2(y_train))"
2100    ]
2101   },
2102   {
2103    "cell_type": "code",
2104    "execution_count": null,
2105    "metadata": {},
2106    "outputs": [],
2107    "source": [
2108     "## Check 2: Untrained RNN initialized with same weights\n",
2109     "\n",
2110     "hash2(model_fit.get_weights())"
2111    ]
2112   },
2113   {
2114    "cell_type": "code",
2115    "execution_count": null,
2116    "metadata": {
2117     "id": "ClBMYe8Lqr7P"
2118    },
2119    "outputs": [],
2120    "source": [
2121     "model_predict=moisture_rnn.create_RNN_2(hidden_units=hidden_units, dense_units=dense_units,  \n",
2122     "                        input_shape=(hours,features),stateful = False,\n",
2123     "                        return_sequences=True,\n",
2124     "                        activation=activation,dense_layers=dense_layers)\n",
2125     "# model_predict=create_RNN_sequences(hidden_units=1, dense_units=1, input_shape=(hours,1), \n",
2126     "#                        activation=['linear', 'linear'])\n",
2127     "# print('model_predict input shape',Et.shape,'output shape',model_predict(Et).shape)\n",
2128     "print(model_predict.summary())\n",
2129     "from keras.utils.vis_utils import plot_model\n",
2130     "plot_model(model_predict, to_file='model_plot.png', \n",
2131     "           show_shapes=True, show_layer_names=True)"
2132    ]
2133   },
2134   {
2135    "cell_type": "code",
2136    "execution_count": null,
2137    "metadata": {},
2138    "outputs": [],
2139    "source": [
2140     "## Check 3: Second model initialization same weights\n",
2141     "\n",
2142     "hash2(model_predict.get_weights())"
2143    ]
2144   },
2145   {
2146    "cell_type": "code",
2147    "execution_count": null,
2148    "metadata": {
2149     "id": "4U0kTEiksNZs"
2150    },
2151    "outputs": [],
2152    "source": [
2153     "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",
2154     "w=model_fit.get_weights()\n",
2155     "for i in range(len(w)):\n",
2156     "  print('weight',i,'shape',w[i].shape,'ndim',w[i].ndim,'given',w_initial[i].shape)\n",
2157     "  for j in range(w[i].shape[0]):\n",
2158     "    if w[i].ndim==2:\n",
2159     "      for k in range(w[i].shape[1]):\n",
2160     "        w[i][j][k]=w_initial[i][0][0]/w[i].shape[0]\n",
2161     "    else:\n",
2162     "      w[i][j]=w_initial[i][0]\n",
2163     "model_fit.set_weights(w)"
2164    ]
2165   },
2166   {
2167    "cell_type": "code",
2168    "execution_count": null,
2169    "metadata": {},
2170    "outputs": [],
2171    "source": [
2172     "## Check 4: weights the same after this step \n",
2173     "\n",
2174     "print(hash2(model_fit.get_weights()))\n",
2175     "print(hash2(x_train))\n",
2176     "print(hash2(y_train))"
2177    ]
2178   },
2179   {
2180    "cell_type": "code",
2181    "execution_count": null,
2182    "metadata": {},
2183    "outputs": [],
2184    "source": [
2185     "# print('model_fit input shape',x_train.shape,'output shape',y_train.shape)\n",
2186     "# print('x_train',x_train)\n",
2187     "# print('y_train',y_train)"
2188    ]
2189   },
2190   {
2191    "cell_type": "code",
2192    "execution_count": null,
2193    "metadata": {},
2194    "outputs": [],
2195    "source": [
2196     "# reproducibility.set_seed()"
2197    ]
2198   },
2199   {
2200    "cell_type": "code",
2201    "execution_count": null,
2202    "metadata": {},
2203    "outputs": [],
2204    "source": [
2205     "# model_fit.get_weights()"
2206    ]
2207   },
2208   {
2209    "cell_type": "code",
2210    "execution_count": null,
2211    "metadata": {},
2212    "outputs": [],
2213    "source": [
2214     "model_fit.fit(x_train, y_train, epochs=5000, verbose=0,batch_size=samples)\n",
2215     "w_fitted=model_fit.get_weights()\n",
2216     "for i in range(len(w)):\n",
2217     "  print('weight',i,' exact:',w_exact[i],':  initial:',w_initial[i],' fitted:',w_fitted[i])"
2218    ]
2219   },
2220   {
2221    "cell_type": "code",
2222    "execution_count": null,
2223    "metadata": {},
2224    "outputs": [],
2225    "source": [
2226     "## Check 5: Weights NOT the same after fitting\n",
2227     "\n",
2228     "hash2(model_fit.get_weights())"
2229    ]
2230   },
2231   {
2232    "cell_type": "code",
2233    "execution_count": null,
2234    "metadata": {
2235     "id": "o10lIOl4sndv"
2236    },
2237    "outputs": [],
2238    "source": [
2239     "# evaluate model\n",
2240     "model_predict.set_weights(w_fitted)\n",
2241     "x_input=np.reshape(Et,(1, hours, 2))\n",
2242     "y_output = model_predict.predict(x_input)\n",
2243     "print('x_input.shape=',x_input.shape,'y_output.shape=',y_output.shape)\n",
2244     "print(shift)\n",
2245     "m = np.reshape(y_output,hours)\n",
2246     "print('weights=',w)\n",
2247     "if scale:\n",
2248     "    print('scaling')\n",
2249     "    m = scalery.inverse_transform(m)\n",
2250     "m = np.reshape(m,hours)\n",
2251     "hour=np.array(range(hours))\n",
2252     "title=\"First RNN forecast\"\n",
2253     "plt.figure(figsize=(16,4))\n",
2254     "plt.plot(hour,Ed,linestyle='--',c='r',label='Drying equilibrium')\n",
2255     "plt.plot(hour,Ew,linestyle='--',c='b',label='Wetting equilibrium')\n",
2256     "# print(len(hour),len(m_f))\n",
2257     "plt.scatter(hour,data,c='b',label='data=10-h fuel data')\n",
2258     "if m is not None:\n",
2259     "    plt.plot(hour[:h2],m[:h2],linestyle='-',c='k',label='m=filtered')\n",
2260     "    plt.plot(hour[h2:hours],m[h2:hours],linestyle='-',c='r',label='m=forecast')\n",
2261     "plt.title(title) \n",
2262     "plt.legend()\n"
2263    ]
2264   },
2265   {
2266    "cell_type": "code",
2267    "execution_count": null,
2268    "metadata": {},
2269    "outputs": [],
2270    "source": [
2271     "rnn_orig={'title':'RNN fitting and prediction - original','hours':hours,'h2':h2,'Ed':Ed,'Ew':Ew,'rain':rain,\n",
2272     "            'fm':obs_data,'m':m}\n",
2273     "# 'w_exact':w_exact,'w_initial':w_initial,'w_fitted':w_fitted\n",
2274     "to_json(rnn_orig,'rnn_orig.json')"
2275    ]
2276   },
2277   {
2278    "cell_type": "code",
2279    "execution_count": null,
2280    "metadata": {},
2281    "outputs": [],
2282    "source": [
2283     "hash2(rnn_orig, ['Ed', 'Ew', 'fm', 'rain'])"
2284    ]
2285   },
2286   {
2287    "cell_type": "code",
2288    "execution_count": null,
2289    "metadata": {},
2290    "outputs": [],
2291    "source": [
2292     "import pandas as pd\n",
2293     "print(np.sum(pd.util.hash_array(rnn_orig['m'])))\n",
2294     "hash2(rnn_orig['m'])"
2295    ]
2296   },
2297   {
2298    "cell_type": "code",
2299    "execution_count": null,
2300    "metadata": {
2301     "tags": []
2302    },
2303    "outputs": [],
2304    "source": [
2305     "print(rnn_orig['m'][0])"
2306    ]
2307   },
2308   {
2309    "cell_type": "code",
2310    "execution_count": null,
2311    "metadata": {},
2312    "outputs": [],
2313    "source": [
2314     "print(hash2(rnn_orig['Ed']))\n",
2315     "print(rnn_orig['Ed'][0])"
2316    ]
2317   },
2318   {
2319    "cell_type": "code",
2320    "execution_count": null,
2321    "metadata": {
2322     "id": "mrWioCJVuU-G"
2323    },
2324    "outputs": [],
2325    "source": [
2326     "# plot subinterval only\n",
2327     "def plot_int(lb=0,ub=hours,title=\"RNN Prediction\"):\n",
2328     "  hour=np.array(range(hours))\n",
2329     "  plt.figure(figsize=(16,4))\n",
2330     "  plt.plot(hour[lb:ub],Ed[lb:ub],linestyle='--',c='r',label='Drying equilibrium')\n",
2331     "  plt.plot(hour[lb:ub],Ew[lb:ub],linestyle='--',c='b',label='Wetting equilibrium')\n",
2332     "  plt.plot(hour[lb:ub],data[lb:ub],linestyle='-',c='b',label='RAWS fuel moisture data')\n",
2333     "  if lb <= h2:\n",
2334     "    ub1 = min(h2,ub)\n",
2335     "    plt.plot(hour[lb:ub1],m[lb:ub1],linestyle='-',c='k',label='Fuel moisture fitted')\n",
2336     "  if ub >= h2:\n",
2337     "    lb1 = max(h2,lb)\n",
2338     "    plt.plot(hour[lb1:ub],m[lb1:ub],linestyle='-',c='r',label='Fuel moisture prediction')\n",
2339     "  plt.title(title) \n",
2340     "  plt.legend()"
2341    ]
2342   },
2343   {
2344    "cell_type": "code",
2345    "execution_count": null,
2346    "metadata": {
2347     "id": "qmGPeG61uqGI"
2348    },
2349    "outputs": [],
2350    "source": [
2351     "plot_int(0,600,title='RNN fitting and prediction')  # again the whole thing"
2352    ]
2353   },
2354   {
2355    "cell_type": "code",
2356    "execution_count": null,
2357    "metadata": {
2358     "id": "SwnOSJlOuvAA"
2359    },
2360    "outputs": [],
2361    "source": [
2362     "plot_int(0,300,title='RNN Fitting') "
2363    ]
2364   },
2365   {
2366    "cell_type": "code",
2367    "execution_count": null,
2368    "metadata": {
2369     "id": "EqCZD7uCvDrS"
2370    },
2371    "outputs": [],
2372    "source": [
2373     "plot_int(300,500)"
2374    ]
2375   },
2376   {
2377    "cell_type": "code",
2378    "execution_count": null,
2379    "metadata": {
2380     "id": "hYgLAXpUvSLo"
2381    },
2382    "outputs": [],
2383    "source": [
2384     "plot_int(500,800)"
2385    ]
2386   },
2387   {
2388    "cell_type": "code",
2389    "execution_count": null,
2390    "metadata": {},
2391    "outputs": [],
2392    "source": [
2393     "type(model_predict)"
2394    ]
2395   },
2396   {
2397    "cell_type": "markdown",
2398    "metadata": {
2399     "id": "gVQxv9Blc91t"
2400    },
2401    "source": [
2402     "### 4.4 A comment on the information flow in the Kalman filter and in neural networks"
2403    ]
2404   },
2405   {
2406    "cell_type": "markdown",
2407    "metadata": {
2408     "id": "CZmR4HPAc91t"
2409    },
2410    "source": []
2411   },
2412   {
2413    "cell_type": "markdown",
2414    "metadata": {
2415     "id": "_g_OTEg6ePb9"
2416    },
2417    "source": [
2418     "## 5. Conclusion"
2419    ]
2420   },
2421   {
2422    "cell_type": "markdown",
2423    "metadata": {
2424     "id": "aNxw7xI3FqFt"
2425    },
2426    "source": [
2427     "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."
2428    ]
2429   },
2430   {
2431    "cell_type": "markdown",
2432    "metadata": {
2433     "id": "IWpmDwUPGElR"
2434    },
2435    "source": [
2436     "## Contributions of authors "
2437    ]
2438   },
2439   {
2440    "cell_type": "markdown",
2441    "metadata": {
2442     "id": "jujW1VFgGOCn"
2443    },
2444    "source": [
2445     "Not applicable."
2446    ]
2447   },
2448   {
2449    "cell_type": "markdown",
2450    "metadata": {
2451     "id": "HWslw7HmGZmP"
2452    },
2453    "source": [
2454     "## Acknowledgements"
2455    ]
2456   },
2457   {
2458    "cell_type": "markdown",
2459    "metadata": {
2460     "id": "xubqDAV2GjkZ"
2461    },
2462    "source": [
2463     "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."
2464    ]
2465   },
2466   {
2467    "cell_type": "markdown",
2468    "metadata": {
2469     "id": "ZsNZxOv7c91t"
2470    },
2471    "source": [
2472     "## References"
2473    ]
2474   },
2475   {
2476    "cell_type": "markdown",
2477    "metadata": {
2478     "id": "vFY-iS1Wc91t"
2479    },
2480    "source": [
2481     "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",
2482     "\n",
2483     "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",
2484     "\n",
2485     "E. Kalnay. *Atmospheric Modeling, Data Assimilation and Predictability.* Cambridge University Press, 2003. [doi:10.1017/CBO9780511802270](https://doi.org/10.1017/CBO9780511802270)\n",
2486     "\n",
2487     "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",
2488     "\n",
2489     "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",
2490     "\n",
2491     "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"
2492    ]
2493   }
2494  ],
2495  "metadata": {
2496   "accelerator": "GPU",
2497   "colab": {
2498    "collapsed_sections": [],
2499    "provenance": []
2500   },
2501   "gpuClass": "standard",
2502   "kernelspec": {
2503    "display_name": "Python 3 (ipykernel)",
2504    "language": "python",
2505    "name": "python3"
2506   },
2507   "language_info": {
2508    "codemirror_mode": {
2509     "name": "ipython",
2510     "version": 3
2511    },
2512    "file_extension": ".py",
2513    "mimetype": "text/x-python",
2514    "name": "python",
2515    "nbconvert_exporter": "python",
2516    "pygments_lexer": "ipython3",
2517    "version": "3.10.9"
2518   }
2519  },
2520  "nbformat": 4,
2521  "nbformat_minor": 4
2522 }