fmda/moisture_models.py

   1 import numpy as np
   2 import math
   3 import matplotlib.pyplot as plt
   4 import copy
   5
   6 # ODE + Augmented Kalman Filter Code
   7 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   8
   9 def model_decay(m0,E,partials=0,T1=0.1,tlen=1):
  10     # Arguments:
  11     #   m0          fuel moisture content at start dimensionless, unit (1)
  12     #   E           fuel moisture eqilibrium (1)
  13     #   partials=0: return m1 = fuel moisture contents after time tlen (1)
  14     #           =1: return m1, dm0/dm0
  15     #           =2: return m1, dm1/dm0, dm1/dE
  16     #           =3: return m1, dm1/dm0, dm1/dE dm1/dT1
  17     #   T1          1/T, where T is the time constant approaching the equilibrium
  18     #               default 0.1/hour
  19     #   tlen        the time interval length, default 1 hour
  20
  21     exp_t = np.exp(-tlen*T1)                  # compute this subexpression only once
  22     m1 = E + (m0 - E)*exp_t                   # the solution at end
  23     if partials==0:
  24         return m1
  25     dm1_dm0 = exp_t
  26     if partials==1:
  27         return m1, dm1_dm0          # return value and Jacobian
  28     dm1_dE = 1 - exp_t
  29     if partials==2:
  30         return m1, dm1_dm0, dm1_dE
  31     dm1_dT1 = -(m0 - E)*tlen*exp_t            # partial derivative dm1 / dT1
  32     if partials==3:
  33         return m1, dm1_dm0, dm1_dE, dm1_dT1       # return value and all partial derivatives wrt m1 and parameters
  34     raise('Bad arg partials')
  35
  36
  37 def ext_kf(u,P,F,Q=0,d=None,H=None,R=None):
  38     """
  39     One step of the extended Kalman filter.
  40     If there is no data, only advance in time.
  41     :param u:   the state vector, shape n
  42     :param P:   the state covariance, shape (n,n)
  43     :param F:   the model function, args vector u, returns F(u) and Jacobian J(u)
  44     :param Q:   the process model noise covariance, shape (n,n)
  45     :param d:   data vector, shape (m). If none, only advance in time
  46     :param H:   observation matrix, shape (m,n)
  47     :param R:   data error covariance, shape (n,n)
  48     :return ua: the analysis state vector, shape (n)
  49     :return Pa: the analysis covariance matrix, shape (n,n)
  50     """
  51     def d2(a):
  52         return np.atleast_2d(a) # convert to at least 2d array
  53
  54     def d1(a):
  55         return np.atleast_1d(a) # convert to at least 1d array
  56
  57     # forecast
  58     uf, J  = F(u)          # advance the model state in time and get the Jacobian
  59     uf = d1(uf)            # if scalar, make state a 1D array
  60     J = d2(J)              # if scalar, make jacobian a 2D array
  61     P = d2(P)              # if scalar, make Jacobian as 2D array
  62     Pf  = d2(J.T @ P) @ J + Q  # advance the state covariance Pf = J' * P * J + Q
  63     # analysis
  64     if d is None or not d.size :  # no data, no analysis
  65         return uf, Pf
  66     # K = P H' * inverse(H * P * H' + R) = (inverse(H * P * H' + R)*(H P))'
  67     H = d2(H)
  68     HP  = d2(H @ P)            # precompute a part used twice
  69     K   = d2(np.linalg.solve( d2(HP @ H.T) + R, HP)).T  # Kalman gain
  70     # print('H',H)
  71     # print('K',K)
  72     res = d1(H @ d1(uf) - d)          # res = H*uf - d
  73     ua = uf - K @ res # analysis mean uf - K*res
  74     Pa = Pf - K @ d2(H @ P)        # analysis covariance
  75     return ua, d2(Pa)
  76
  77 ### Define model function with drying, wetting, and rain equilibria
  78
  79 # Parameters
  80 r0 = 0.05                                   # threshold rainfall [mm/h]
  81 rs = 8.0                                    # saturation rain intensity [mm/h]
  82 Tr = 14.0                                   # time constant for rain wetting model [h]
  83 S = 250                                     # saturation intensity [dimensionless]
  84 T = 10.0                                    # time constant for wetting/drying
  85
  86 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  87
  88 def model_moisture(m0,Eqd,Eqw,r,t=None,partials=0,T=10.0,tlen=1.0):
  89     # arguments:
  90     # m0         starting fuel moistureb (%s
  91     # Eqd        drying equilibrium      (%)
  92     # Eqw        wetting equilibrium     (%)
  93     # r          rain intensity          (mm/h)
  94     # t          time
  95     # partials = 0, 1, 2
  96     # returns: same as model_decay
  97     #   if partials==0: m1 = fuel moisture contents after time 1 hour
  98     #              ==1: m1, dm1/dm0
  99     #              ==2: m1, dm1/dm0, dm1/dE
 100
 101     if r > r0:
 102         # print('raining')
 103         E = S
 104         T1 =  (1.0 - np.exp(- (r - r0) / rs)) / Tr
 105     elif m0 <= Eqw:
 106         # print('wetting')
 107         E=Eqw
 108         T1 = 1.0/T
 109     elif m0 >= Eqd:
 110         # print('drying')
 111         E=Eqd
 112         T1 = 1.0/T
 113     else: # no change'
 114         E = m0
 115         T1=0.0
 116     exp_t = np.exp(-tlen*T1)
 117     m1 = E + (m0 - E)*exp_t
 118     dm1_dm0 = exp_t
 119     dm1_dE = 1 - exp_t
 120     #if t>=933 and t < 940:
 121     #  print('t,Eqw,Eqd,r,T1,E,m0,m1,dm1_dm0,dm1_dE',
 122     #        t,Eqw,Eqd,r,T1,E,m0,m1,dm1_dm0,dm1_dE)
 123     if partials==0:
 124         return m1
 125     if partials==1:
 126         return m1, dm1_dm0
 127     if partials==2:
 128         return m1, dm1_dm0, dm1_dE
 129     raise('bad partials')
 130
 131 def model_augmented(u0,Ed,Ew,r,t):
 132     # state u is the vector [m,dE] with dE correction to equilibria Ed and Ew at t
 133     #
 134     m0, Ec = u0  # decompose state u0
 135     # reuse model_moisture(m0,Eqd,Eqw,r,partials=0):
 136     # arguments:
 137     # m0         starting fuel moistureb (1)
 138     # Ed         drying equilibrium      (1)
 139     # Ew         wetting equilibrium     (1)
 140     # r          rain intensity          (mm/h)
 141     # partials = 0, 1, 2
 142     # returns: same as model_decay
 143     #   if partials==0: m1 = fuel moisture contents after time 1 hour
 144     #              ==1: m1, dm0/dm0
 145     #              ==2: m1, dm1/dm0, dm1/dE
 146     m1, dm1_dm0, dm1_dE  = model_moisture(m0,Ed + Ec, Ew + Ec, r, t, partials=2)
 147     u1 = np.array([m1,Ec])   # dE is just copied
 148     J =  np.array([[dm1_dm0, dm1_dE],
 149                    [0.     ,     1.]])
 150     return u1, J
 151
 152
 153 ### Default Uncertainty Matrices
 154 Q = np.array([[1e-3, 0.],
 155             [0,  1e-3]]) # process noise covariance
 156 H = np.array([[1., 0.]])  # first component observed
 157 R = np.array([1e-3]) # data variance
 158
 159 def run_augmented_kf(dat0,h2=None,hours=None, H=H, Q=Q, R=R):
 160     dat = copy.deepcopy(dat0)
 161
 162     if h2 is None:
 163         h2 = int(dat['h2'])
 164     if hours is None:
 165         hours = int(dat['hours'])
 166
 167     d = dat['y']
 168     feats = dat['features_list']
 169     Ed = dat['X'][:,feats.index('Ed')]
 170     Ew = dat['X'][:,feats.index('Ew')]
 171     rain = dat['X'][:,feats.index('rain')]
 172
 173     u = np.zeros((2,hours))
 174     u[:,0]=[0.1,0.0]       # initialize,background state
 175     P = np.zeros((2,2,hours))
 176     P[:,:,0] = np.array([[1e-3, 0.],
 177                       [0.,  1e-3]]) # background state covariance
 178     # Q = np.array([[1e-3, 0.],
 179     #             [0,  1e-3]]) # process noise covariance
 180     # H = np.array([[1., 0.]])  # first component observed
 181     # R = np.array([1e-3]) # data variance
 182
 183     for t in range(1,h2):
 184       # use lambda construction to pass additional arguments to the model
 185         u[:,t],P[:,:,t] = ext_kf(u[:,t-1],P[:,:,t-1],
 186                                   lambda uu: model_augmented(uu,Ed[t],Ew[t],rain[t],t),
 187                                   Q,d[t],H=H,R=R)
 188       # print('time',t,'data',d[t],'filtered',u[0,t],'Ec',u[1,t])
 189     for t in range(h2,hours):
 190         u[:,t],P[:,:,t] = ext_kf(u[:,t-1],P[:,:,t-1],
 191                                   lambda uu: model_augmented(uu,Ed[t],Ew[t],rain[t],t),
 192                                   Q*0.0)
 193       # print('time',t,'data',d[t],'forecast',u[0,t],'Ec',u[1,t])
 194     return u
 195
 196 # General Machine Learning Models
 197 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 198
 199 class MLModel(ABC):
 200     def __init__(self, params: dict):
 201         self.params = params
 202         if type(self) is MLModel:
 203             raise TypeError("MLModel is an abstract class and cannot be instantiated directly")
 204         super().__init__()
 205
 206     @abstractmethod
 207     def fit(self, X_train, y_train, weights=None):
 208         pass
 209
 210     @abstractmethod
 211     def predict(self, X):
 212         pass
 213
 214     def eval(self, X_test, y_test):
 215         preds = self.predict(X_test)
 216         rmse = np.sqrt(mean_squared_error(y_test, preds))
 217         rmse_ros = np.sqrt(mean_squared_error(ros_3wind(y_test), ros_3wind(preds)))
 218         print(f"Test RMSE: {rmse}")
 219         print(f"Test RMSE (ROS): {rmse_ros}")
 220         return rmse, rmse_ros
 221
 222 class XGB(MLModel):
 223     def __init__(self, params: dict):
 224         super().__init__(params)
 225         self.model = XGBRegressor(**self.params)
 226
 227     def fit(self, X_train, y_train, weights=None):
 228         print(f"Training XGB with params: {self.params}")
 229         self.model.fit(X_train, y_train, sample_weight=weights)
 230
 231     def predict(self, X):
 232         print("Predicting with XGB")
 233         preds = self.model.predict(X)
 234         return preds
 235
 236 class RF(MLModel):
 237     def __init__(self, params: dict):
 238         super().__init__(params)
 239         self.model = RandomForestRegressor(**self.params)
 240
 241     def fit(self, X_train, y_train, weights=None):
 242         print(f"Training RF with params: {self.params}")
 243         self.model.fit(X_train, y_train, sample_weight=weights)
 244
 245     def predict(self, X):
 246         print("Predicting with RF")
 247         preds = self.model.predict(X)
 248         return preds
 249
 250 class LM(MLModel):
 251     def __init__(self, params: dict):
 252         super().__init__(params)
 253         self.model = LinearRegression(**self.params)
 254
 255     def fit(self, X_train, y_train, weights=None):
 256         self.model.fit(X_train, y_train, sample_weight=weights)
 257         print(f"Training LM with params: {self.params}")
 258
 259     def predict(self, X):
 260         print("Predicting with LM")
 261         preds = self.model.predict(X)
 262         return preds
 263
 264
 265
 266
 267
 268
 269
 270
 271
 272
 273
 274
 275
 276
 277
 278
 279
 280
 281