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