From 6d7d12e0b733703bb2922800aa1fc045fcd13eb3 Mon Sep 17 00:00:00 2001 From: Jan Date: Thu, 5 Jan 2023 18:08:55 -0700 Subject: [PATCH] splitting moisture_rnn.py off moisture_models.py --- fmda/fmda_rnn_rain.ipynb | 2 +- fmda/moisture_models.py | 216 --------------------------- fmda/{moisture_models.py => moisture_rnn.py} | 179 ---------------------- 3 files changed, 1 insertion(+), 396 deletions(-) copy fmda/{moisture_models.py => moisture_rnn.py} (55%) diff --git a/fmda/fmda_rnn_rain.ipynb b/fmda/fmda_rnn_rain.ipynb index 5153006..917dfc9 100644 --- a/fmda/fmda_rnn_rain.ipynb +++ b/fmda/fmda_rnn_rain.ipynb @@ -33,7 +33,7 @@ "import data_funcs as datf\n", "from data_funcs import format_raws, retrieve_raws, format_precip, fixnan\n", "import moisture_models as mod\n", - "from moisture_models import create_RNN, create_RNN_2, staircase, seq2batches, create_rnn_data, train_rnn, rnn_predict\n", + "from moisture_rnn import create_RNN, create_RNN_2, staircase, seq2batches, create_rnn_data, train_rnn, rnn_predict\n", "\n", "meso_token=\"4192c18707b848299783d59a9317c6e1\"\n", "m=Meso(meso_token)" diff --git a/fmda/moisture_models.py b/fmda/moisture_models.py index 3655dbc..93b3e53 100644 --- a/fmda/moisture_models.py +++ b/fmda/moisture_models.py @@ -195,219 +195,3 @@ def run_augmented_kf(d,Ed,Ew,rain,h2,hours, H=H, Q=Q, R=R): # print('time',t,'data',d[t],'forecast',u[0,t],'Ec',u[1,t]) return u -#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -## RNN Model Funcs - -def create_RNN(hidden_units, dense_units, input_shape, activation): - inputs = tf.keras.Input(shape=input_shape) - # https://stackoverflow.com/questions/43448029/how-can-i-print-the-values-of-keras-tensors - # inputs2 = K.print_tensor(inputs, message='inputs = ') # change allso inputs to inputs2 below, must be used - x = tf.keras.layers.SimpleRNN(hidden_units, input_shape=input_shape, - activation=activation[0])(inputs) - outputs = tf.keras.layers.Dense(dense_units, activation=activation[1])(x) - model = tf.keras.Model(inputs=inputs, outputs=outputs) - model.compile(loss='mean_squared_error', optimizer='adam') - return model - -def staircase(x,y,timesteps,trainsteps,return_sequences=False, verbose = False): - # x [trainsteps+forecaststeps,features] all inputs - # y [trainsteps,outputs] - # timesteps: split x and y into samples length timesteps, shifted by 1 - # trainsteps: number of timesteps to use for training, no more than y.shape[0] - vprint('shape x = ',x.shape) - vprint('shape y = ',y.shape) - vprint('timesteps=',timesteps) - vprint('trainsteps=',trainsteps) - outputs = y.shape[1] - features = x.shape[1] - forecaststeps = x.shape[0]-trainsteps - samples = trainsteps-timesteps+1 - vprint('staircase: samples=',samples,'timesteps=',timesteps,'features=',features) - x_train = np.empty([samples, timesteps, features]) - vprint('return_sequences=',return_sequences) - if return_sequences: - vprint('returning all timesteps in a sample') - y_train = np.empty([samples, timesteps, outputs]) # all - for i in range(samples): - for k in range(timesteps): - for j in range(features): - x_train[i,k,j] = x[i+k,j] - for j in range(outputs): - y_train[i,k,j] = y[i+k,j] - else: - vprint('returning only the last timestep in a sample') - y_train = np.empty([samples, outputs]) - for i in range(samples): - for j in range(features): - for k in range(timesteps): - x_train[i,k,j] = x[i+k,j] - for j in range(outputs): - y_train[i,j] = y[i+timesteps-1,j] - - return x_train, y_train - - -def seq2batches(x,y,timesteps,trainsteps): - # x [trainsteps+forecaststeps,features] all inputs - # y [trainsteps,outputs] - # timesteps: split x and y into samples length timesteps, shifted by 1 - # trainsteps: number of timesteps to use for training, no more than y.shape[0] - print('shape x = ',x.shape) - print('shape y = ',y.shape) - print('timesteps=',timesteps) - print('trainsteps=',trainsteps) - outputs = y.shape[1] - features = x.shape[1] - samples= trainsteps - timesteps + 1 - print('samples=',samples) - x_train = np.empty([samples, timesteps, features]) - y_train = np.empty([samples, timesteps, outputs]) # only the last - print('samples=',samples,' timesteps=',timesteps, - ' features=',features,' outputs=',outputs) - for i in range(samples): - for k in range(timesteps): - for j in range(features): - x_train[i,k,j] = x[i+k,j] - for j in range(outputs): - y_train[i,k,j] = y[i+k,j] # return sequences - return x_train, y_train - -def create_RNN_2(hidden_units, dense_units, activation, stateful=False, - batch_shape=None, input_shape=None, dense_layers=1, - rnn_layers=1,return_sequences=False, - initial_state=None, verbose = True): - if stateful: - inputs = tf.keras.Input(batch_shape=batch_shape) - else: - inputs = tf.keras.Input(shape=input_shape) - # https://stackoverflow.com/questions/43448029/how-can-i-print-the-values-of-keras-tensors - # inputs2 = K.print_tensor(inputs, message='inputs = ') # change allso inputs to inputs2 below, must be used - x = inputs - for i in range(rnn_layers): - x = tf.keras.layers.SimpleRNN(hidden_units,activation=activation[0], - stateful=stateful,return_sequences=return_sequences)(x - # ,initial_state=initial_state - ) - # x = tf.keras.layers.Dense(hidden_units, activation=activation[1])(x) - for i in range(dense_layers): - x = tf.keras.layers.Dense(dense_units, activation=activation[1])(x) - model = tf.keras.Model(inputs=inputs, outputs=x) - model.compile(loss='mean_squared_error', optimizer='adam') - return model - -def create_rnn_data(dat, hours, h2, scale = False, verbose = False): - Ew = dat['Ew'] - Ed = dat['Ed'] - rain = dat['rain'] - fm = dat['fm'] - temp = dat['temp'] - - # Average Equilibrium - E = (Ed + Ew)/2 - - # transform as 2D, (timesteps, features) and (timesteps, outputs) - Et = np.reshape(E,[E.shape[0],1]) - datat = np.reshape(fm,[fm.shape[0],1]) - - # Scale Data if required - scale=False - if scale: - scalerx = MinMaxScaler() - scalerx.fit(Et) - Et = scalerx.transform(Et) - scalery = MinMaxScaler() - scalery.fit(datat) - datat = scalery.transform(datat) - - # split data - x_train, y_train = staircase(Et,datat,timesteps=5,trainsteps=h2, - return_sequences=False, verbose = verbose) - vprint('x_train shape=',x_train.shape) - samples, timesteps, features = x_train.shape - vprint('y_train shape=',y_train.shape) - - h0 = tf.convert_to_tensor(datat[:samples],dtype=tf.float32) - - # Set up return dictionary - - rnn_dat = { - 'x_train': x_train, - 'y_train': y_train, - 'Et': Et, - 'samples': samples, - 'timesteps': timesteps, - 'features': features, - 'h0': h0 - } - - return rnn_dat - -def train_rnn(rnn_dat, hours, activation, hidden_units, dense_units, dense_layers, verbose = False): - - samples = rnn_dat['samples'] - features = rnn_dat['features'] - timesteps = rnn_dat['timesteps'] - - model_fit=create_RNN_2(hidden_units=hidden_units, - dense_units=dense_units, - batch_shape=(samples,timesteps,features), - stateful=True, - return_sequences=False, - # initial_state=h0, - activation=activation, - dense_layers=dense_layers) - - Et = rnn_dat['Et'] - model_predict=create_RNN_2(hidden_units=hidden_units, dense_units=dense_units, - input_shape=(hours,features),stateful = False, - return_sequences=True, - activation=activation,dense_layers=dense_layers) - - vprint('model_predict input shape',Et.shape,'output shape',model_predict(Et).shape) - if verbose: print(model_predict.summary()) - - x_train = rnn_dat['x_train'] - y_train = rnn_dat['y_train'] - - # fitting - DeltaE = 0 - 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])] - 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])] - w=model_fit.get_weights() - for i in range(len(w)): - vprint('weight',i,'shape',w[i].shape,'ndim',w[i].ndim,'given',w_initial[i].shape) - for j in range(w[i].shape[0]): - if w[i].ndim==2: - for k in range(w[i].shape[1]): - w[i][j][k]=w_initial[i][0][0]/w[i].shape[0] - else: - w[i][j]=w_initial[i][0] - model_fit.set_weights(w) - model_fit.fit(x_train, y_train, epochs=5000,batch_size=samples, verbose=0) - w_fitted=model_fit.get_weights() - for i in range(len(w)): - vprint('weight',i,' exact:',w_exact[i],': initial:',w_initial[i],' fitted:',w_fitted[i]) - - model_predict.set_weights(w_fitted) - - return model_predict - - -def rnn_predict(model, rnn_dat, hours, scale = False, verbose = False): - scale = False - # model.set_weights(w_fitted) - x_input=np.reshape(rnn_dat['Et'],(1, hours, 1)) - y_output = model.predict(x_input, verbose = verbose) - - vprint('x_input.shape=',x_input.shape,'y_output.shape=',y_output.shape) - - m = np.reshape(y_output,hours) - # print('weights=',w) - if scale: - vprint('scaling') - m = scalery.inverse_transform(m) - m = np.reshape(m,hours) - - return m - diff --git a/fmda/moisture_models.py b/fmda/moisture_rnn.py similarity index 55% copy from fmda/moisture_models.py copy to fmda/moisture_rnn.py index 3655dbc..a7fb756 100644 --- a/fmda/moisture_models.py +++ b/fmda/moisture_rnn.py @@ -18,185 +18,6 @@ def vprint(*args): print(s, end=' ') print(args[-1]) - -def model_decay(m0,E,partials=0,T1=0.1,tlen=1): - # Arguments: - # m0 fuel moisture content at start dimensionless, unit (1) - # E fuel moisture eqilibrium (1) - # partials=0: return m1 = fuel moisture contents after time tlen (1) - # =1: return m1, dm0/dm0 - # =2: return m1, dm1/dm0, dm1/dE - # =3: return m1, dm1/dm0, dm1/dE dm1/dT1 - # T1 1/T, where T is the time constant approaching the equilibrium - # default 0.1/hour - # tlen the time interval length, default 1 hour - - exp_t = np.exp(-tlen*T1) # compute this subexpression only once - m1 = E + (m0 - E)*exp_t # the solution at end - if partials==0: - return m1 - dm1_dm0 = exp_t - if partials==1: - return m1, dm1_dm0 # return value and Jacobian - dm1_dE = 1 - exp_t - if partials==2: - return m1, dm1_dm0, dm1_dE - dm1_dT1 = -(m0 - E)*tlen*exp_t # partial derivative dm1 / dT1 - if partials==3: - return m1, dm1_dm0, dm1_dE, dm1_dT1 # return value and all partial derivatives wrt m1 and parameters - raise('Bad arg partials') - - -def ext_kf(u,P,F,Q=0,d=None,H=None,R=None): - """ - One step of the extended Kalman filter. - If there is no data, only advance in time. - :param u: the state vector, shape n - :param P: the state covariance, shape (n,n) - :param F: the model function, args vector u, returns F(u) and Jacobian J(u) - :param Q: the process model noise covariance, shape (n,n) - :param d: data vector, shape (m). If none, only advance in time - :param H: observation matrix, shape (m,n) - :param R: data error covariance, shape (n,n) - :return ua: the analysis state vector, shape (n) - :return Pa: the analysis covariance matrix, shape (n,n) - """ - def d2(a): - return np.atleast_2d(a) # convert to at least 2d array - - def d1(a): - return np.atleast_1d(a) # convert to at least 1d array - - # forecast - uf, J = F(u) # advance the model state in time and get the Jacobian - uf = d1(uf) # if scalar, make state a 1D array - J = d2(J) # if scalar, make jacobian a 2D array - P = d2(P) # if scalar, make Jacobian as 2D array - Pf = d2(J.T @ P) @ J + Q # advance the state covariance Pf = J' * P * J + Q - # analysis - if d is None or not d.size : # no data, no analysis - return uf, Pf - # K = P H' * inverse(H * P * H' + R) = (inverse(H * P * H' + R)*(H P))' - H = d2(H) - HP = d2(H @ P) # precompute a part used twice - K = d2(np.linalg.solve( d2(HP @ H.T) + R, HP)).T # Kalman gain - # print('H',H) - # print('K',K) - res = d1(H @ d1(uf) - d) # res = H*uf - d - ua = uf - K @ res # analysis mean uf - K*res - Pa = Pf - K @ d2(H @ P) # analysis covariance - return ua, d2(Pa) - -### Define model function with drying, wetting, and rain equilibria - -# Parameters -r0 = 0.05 # threshold rainfall [mm/h] -rs = 8.0 # saturation rain intensity [mm/h] -Tr = 14.0 # time constant for rain wetting model [h] -S = 250 # saturation intensity [dimensionless] -T = 10.0 # time constant for wetting/drying - -#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -def model_moisture(m0,Eqd,Eqw,r,t,partials=0,T=10.0,tlen=1.0): - # arguments: - # m0 starting fuel moistureb (%s - # Eqd drying equilibrium (%) - # Eqw wetting equilibrium (%) - # r rain intensity (mm/h) - # t time - # partials = 0, 1, 2 - # returns: same as model_decay - # if partials==0: m1 = fuel moisture contents after time 1 hour - # ==1: m1, dm1/dm0 - # ==2: m1, dm1/dm0, dm1/dE - - if r > r0: - # print('raining') - E = S - T1 = (1.0 - np.exp(- (r - r0) / rs)) / Tr - elif m0 <= Eqw: - # print('wetting') - E=Eqw - T1 = 1.0/T - elif m0 >= Eqd: - # print('drying') - E=Eqd - T1 = 1.0/T - else: # no change' - E = m0 - T1=0.0 - exp_t = np.exp(-tlen*T1) - m1 = E + (m0 - E)*exp_t - dm1_dm0 = exp_t - dm1_dE = 1 - exp_t - #if t>=933 and t < 940: - # print('t,Eqw,Eqd,r,T1,E,m0,m1,dm1_dm0,dm1_dE', - # t,Eqw,Eqd,r,T1,E,m0,m1,dm1_dm0,dm1_dE) - if partials==0: - return m1 - if partials==1: - return m1, dm1_dm0 - if partials==2: - return m1, dm1_dm0, dm1_dE - raise('bad partials') - -#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -def model_augmented(u0,Ed,Ew,r,t): - # state u is the vector [m,dE] with dE correction to equilibria Ed and Ew at t - # - m0, Ec = u0 # decompose state u0 - # reuse model_moisture(m0,Eqd,Eqw,r,partials=0): - # arguments: - # m0 starting fuel moistureb (1) - # Ed drying equilibrium (1) - # Ew wetting equilibrium (1) - # r rain intensity (mm/h) - # partials = 0, 1, 2 - # returns: same as model_decay - # if partials==0: m1 = fuel moisture contents after time 1 hour - # ==1: m1, dm0/dm0 - # ==2: m1, dm1/dm0, dm1/dE - m1, dm1_dm0, dm1_dE = model_moisture(m0,Ed + Ec, Ew + Ec, r, t, partials=2) - u1 = np.array([m1,Ec]) # dE is just copied - J = np.array([[dm1_dm0, dm1_dE], - [0. , 1.]]) - return u1, J - -#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -### Default Uncertainty Matrices -Q = np.array([[1e-3, 0.], - [0, 1e-3]]) # process noise covariance -H = np.array([[1., 0.]]) # first component observed -R = np.array([1e-3]) # data variance - -def run_augmented_kf(d,Ed,Ew,rain,h2,hours, H=H, Q=Q, R=R): - u = np.zeros((2,hours)) - u[:,0]=[0.1,0.0] # initialize,background state - P = np.zeros((2,2,hours)) - P[:,:,0] = np.array([[1e-3, 0.], - [0., 1e-3]]) # background state covariance - # Q = np.array([[1e-3, 0.], - # [0, 1e-3]]) # process noise covariance - # H = np.array([[1., 0.]]) # first component observed - # R = np.array([1e-3]) # data variance - - for t in range(1,h2): - # use lambda construction to pass additional arguments to the model - u[:,t],P[:,:,t] = ext_kf(u[:,t-1],P[:,:,t-1], - lambda uu: model_augmented(uu,Ed[t],Ew[t],rain[t],t), - Q,d[t],H=H,R=R) - # print('time',t,'data',d[t],'filtered',u[0,t],'Ec',u[1,t]) - for t in range(h2,hours): - u[:,t],P[:,:,t] = ext_kf(u[:,t-1],P[:,:,t-1], - lambda uu: model_augmented(uu,Ed[t],Ew[t],rain[t],t), - Q*0.0) - # print('time',t,'data',d[t],'forecast',u[0,t],'Ec',u[1,t]) - return u - -#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - ## RNN Model Funcs def create_RNN(hidden_units, dense_units, input_shape, activation): -- 2.11.4.GIT