fmda/moisture_models.py

   1 import numpy as np
   2 import tensorflow as tf
   3 from keras.models import Sequential
   4 from keras.layers import Dense, SimpleRNN
   5 from keras.utils.vis_utils import plot_model
   6 from sklearn.preprocessing import MinMaxScaler
   7 from sklearn.metrics import mean_squared_error
   8 import math
   9 import matplotlib.pyplot as plt
  10 import tensorflow as tf
  11 import keras.backend as K
  12 import tensorflow as tf
  13 from utils import vprint
  14
  15
  16 def model_decay(m0,E,partials=0,T1=0.1,tlen=1):
  17     # Arguments:
  18     #   m0          fuel moisture content at start dimensionless, unit (1)
  19     #   E           fuel moisture eqilibrium (1)
  20     #   partials=0: return m1 = fuel moisture contents after time tlen (1)
  21     #           =1: return m1, dm0/dm0
  22     #           =2: return m1, dm1/dm0, dm1/dE
  23     #           =3: return m1, dm1/dm0, dm1/dE dm1/dT1
  24     #   T1          1/T, where T is the time constant approaching the equilibrium
  25     #               default 0.1/hour
  26     #   tlen        the time interval length, default 1 hour
  27
  28     exp_t = np.exp(-tlen*T1)                  # compute this subexpression only once
  29     m1 = E + (m0 - E)*exp_t                   # the solution at end
  30     if partials==0:
  31         return m1
  32     dm1_dm0 = exp_t
  33     if partials==1:
  34         return m1, dm1_dm0          # return value and Jacobian
  35     dm1_dE = 1 - exp_t
  36     if partials==2:
  37         return m1, dm1_dm0, dm1_dE
  38     dm1_dT1 = -(m0 - E)*tlen*exp_t            # partial derivative dm1 / dT1
  39     if partials==3:
  40         return m1, dm1_dm0, dm1_dE, dm1_dT1       # return value and all partial derivatives wrt m1 and parameters
  41     raise('Bad arg partials')
  42
  43
  44 def ext_kf(u,P,F,Q=0,d=None,H=None,R=None):
  45     """
  46     One step of the extended Kalman filter.
  47     If there is no data, only advance in time.
  48     :param u:   the state vector, shape n
  49     :param P:   the state covariance, shape (n,n)
  50     :param F:   the model function, args vector u, returns F(u) and Jacobian J(u)
  51     :param Q:   the process model noise covariance, shape (n,n)
  52     :param d:   data vector, shape (m). If none, only advance in time
  53     :param H:   observation matrix, shape (m,n)
  54     :param R:   data error covariance, shape (n,n)
  55     :return ua: the analysis state vector, shape (n)
  56     :return Pa: the analysis covariance matrix, shape (n,n)
  57     """
  58     def d2(a):
  59         return np.atleast_2d(a) # convert to at least 2d array
  60
  61     def d1(a):
  62         return np.atleast_1d(a) # convert to at least 1d array
  63
  64     # forecast
  65     uf, J  = F(u)          # advance the model state in time and get the Jacobian
  66     uf = d1(uf)            # if scalar, make state a 1D array
  67     J = d2(J)              # if scalar, make jacobian a 2D array
  68     P = d2(P)              # if scalar, make Jacobian as 2D array
  69     Pf  = d2(J.T @ P) @ J + Q  # advance the state covariance Pf = J' * P * J + Q
  70     # analysis
  71     if d is None or not d.size :  # no data, no analysis
  72         return uf, Pf
  73     # K = P H' * inverse(H * P * H' + R) = (inverse(H * P * H' + R)*(H P))'
  74     H = d2(H)
  75     HP  = d2(H @ P)            # precompute a part used twice
  76     K   = d2(np.linalg.solve( d2(HP @ H.T) + R, HP)).T  # Kalman gain
  77     # print('H',H)
  78     # print('K',K)
  79     res = d1(H @ d1(uf) - d)          # res = H*uf - d
  80     ua = uf - K @ res # analysis mean uf - K*res
  81     Pa = Pf - K @ d2(H @ P)        # analysis covariance
  82     return ua, d2(Pa)
  83
  84 ### Define model function with drying, wetting, and rain equilibria
  85
  86 # Parameters
  87 r0 = 0.05                                   # threshold rainfall [mm/h]
  88 rs = 8.0                                    # saturation rain intensity [mm/h]
  89 Tr = 14.0                                   # time constant for rain wetting model [h]
  90 S = 250                                     # saturation intensity [dimensionless]
  91 T = 10.0                                    # time constant for wetting/drying
  92
  93 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  94
  95 def model_moisture(m0,Eqd,Eqw,r,t=None,partials=0,T=10.0,tlen=1.0):
  96     # arguments:
  97     # m0         starting fuel moistureb (%s
  98     # Eqd        drying equilibrium      (%)
  99     # Eqw        wetting equilibrium     (%)
 100     # r          rain intensity          (mm/h)
 101     # t          time
 102     # partials = 0, 1, 2
 103     # returns: same as model_decay
 104     #   if partials==0: m1 = fuel moisture contents after time 1 hour
 105     #              ==1: m1, dm1/dm0
 106     #              ==2: m1, dm1/dm0, dm1/dE
 107
 108     if r > r0:
 109         # print('raining')
 110         E = S
 111         T1 =  (1.0 - np.exp(- (r - r0) / rs)) / Tr
 112     elif m0 <= Eqw:
 113         # print('wetting')
 114         E=Eqw
 115         T1 = 1.0/T
 116     elif m0 >= Eqd:
 117         # print('drying')
 118         E=Eqd
 119         T1 = 1.0/T
 120     else: # no change'
 121         E = m0
 122         T1=0.0
 123     exp_t = np.exp(-tlen*T1)
 124     m1 = E + (m0 - E)*exp_t
 125     dm1_dm0 = exp_t
 126     dm1_dE = 1 - exp_t
 127     #if t>=933 and t < 940:
 128     #  print('t,Eqw,Eqd,r,T1,E,m0,m1,dm1_dm0,dm1_dE',
 129     #        t,Eqw,Eqd,r,T1,E,m0,m1,dm1_dm0,dm1_dE)
 130     if partials==0:
 131         return m1
 132     if partials==1:
 133         return m1, dm1_dm0
 134     if partials==2:
 135         return m1, dm1_dm0, dm1_dE
 136     raise('bad partials')
 137
 138 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 139
 140 ##  NOT TESTED
 141 def model_moisture_run(Eqd,Eqw,r,hours=None,T=10.0,tlen=1.0):
 142 # for arrays of FMC model input run the fuel moisture model
 143     if hours is None:
 144         hours = min(len(Eqd),len(Eqw),len(r))
 145     m = np.zeros(hours)
 146     m[0]=(Eqd[0]+Eqw[0])/2
 147     for k in range(hours-1):
 148         m[k+1]=model_moisture(m[k],Eqd[k],Eqw[k],r[k],T=T,tlen=tlen)
 149     return m
 150
 151 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 152
 153 def model_augmented(u0,Ed,Ew,r,t):
 154     # state u is the vector [m,dE] with dE correction to equilibria Ed and Ew at t
 155     #
 156     m0, Ec = u0  # decompose state u0
 157     # reuse model_moisture(m0,Eqd,Eqw,r,partials=0):
 158     # arguments:
 159     # m0         starting fuel moistureb (1)
 160     # Ed         drying equilibrium      (1)
 161     # Ew         wetting equilibrium     (1)
 162     # r          rain intensity          (mm/h)
 163     # partials = 0, 1, 2
 164     # returns: same as model_decay
 165     #   if partials==0: m1 = fuel moisture contents after time 1 hour
 166     #              ==1: m1, dm0/dm0
 167     #              ==2: m1, dm1/dm0, dm1/dE
 168     m1, dm1_dm0, dm1_dE  = model_moisture(m0,Ed + Ec, Ew + Ec, r, t, partials=2)
 169     u1 = np.array([m1,Ec])   # dE is just copied
 170     J =  np.array([[dm1_dm0, dm1_dE],
 171                    [0.     ,     1.]])
 172     return u1, J
 173
 174 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 175 ### Default Uncertainty Matrices
 176 Q = np.array([[1e-3, 0.],
 177             [0,  1e-3]]) # process noise covariance
 178 H = np.array([[1., 0.]])  # first component observed
 179 R = np.array([1e-3]) # data variance
 180
 181 def run_augmented_kf(dat,h2=None,hours=None, H=H, Q=Q, R=R):
 182     d = dat['fm']
 183     Ed = dat['Ed']
 184     Ew = dat['Ew']
 185     rain = dat['rain']
 186     if h2 is None:
 187         h2 = int(dat['h2'])
 188     if hours is None:
 189         hours = int(dat['hours'])
 190     u = np.zeros((2,hours))
 191     u[:,0]=[0.1,0.0]       # initialize,background state
 192     P = np.zeros((2,2,hours))
 193     P[:,:,0] = np.array([[1e-3, 0.],
 194                       [0.,  1e-3]]) # background state covariance
 195     # Q = np.array([[1e-3, 0.],
 196     #             [0,  1e-3]]) # process noise covariance
 197     # H = np.array([[1., 0.]])  # first component observed
 198     # R = np.array([1e-3]) # data variance
 199
 200     for t in range(1,h2):
 201       # use lambda construction to pass additional arguments to the model
 202         u[:,t],P[:,:,t] = ext_kf(u[:,t-1],P[:,:,t-1],
 203                                   lambda uu: model_augmented(uu,Ed[t],Ew[t],rain[t],t),
 204                                   Q,d[t],H=H,R=R)
 205       # print('time',t,'data',d[t],'filtered',u[0,t],'Ec',u[1,t])
 206     for t in range(h2,hours):
 207         u[:,t],P[:,:,t] = ext_kf(u[:,t-1],P[:,:,t-1],
 208                                   lambda uu: model_augmented(uu,Ed[t],Ew[t],rain[t],t),
 209                                   Q*0.0)
 210       # print('time',t,'data',d[t],'forecast',u[0,t],'Ec',u[1,t])
 211     return u
 212