fmda/moisture_models.py

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