diff --git a/examples/hwpv/pv1_poly.py b/examples/hwpv/pv1_poly.py index 3077e7b..3996ef6 100644 --- a/examples/hwpv/pv1_poly.py +++ b/examples/hwpv/pv1_poly.py @@ -18,15 +18,20 @@ model_folder = r'./models' idx_in = [0,1,2,3,4,5,6] -idx_out = [7,8,9] +idx_out = [7,8,9,10] +idx_out1 = idx_out[:2] +idx_out2 = idx_out[2:] +idx_out3 = idx_out[3:] class pv1(): def __init__(self, training_config=None, sim_config=None): self.init_to_none() if training_config is not None: self.load_training_config (training_config) + print('train_config') if sim_config is not None: self.load_sim_config (sim_config) + print('sim_config') def init_to_none(self): self.lr = None @@ -58,6 +63,7 @@ def init_to_none(self): def load_training_config(self, filename): fp = open (filename, 'r') + print(filename) config = json.load (fp) fp.close() self.lr = config['lr'] @@ -101,7 +107,7 @@ def read_lti(self, config): dict['b_coeff'][i,j,:] = torch.Tensor(b) block.load_state_dict (dict) -# print ('state dict', block.state_dict()) + print ('state dict', block.state_dict()) return block def read_net(self, config): @@ -140,10 +146,10 @@ def load_sim_config(self, filename, model_only=True): # print ('COL_U', self.COL_U) # print ('COL_Y', self.COL_Y) # print ('t_step', self.t_step) -# print ('F1', self.F1) -# print ('H1', self.H1) -# print ('F2', self.F2) -# print ('normfacs', self.normfacs) + # print ('F1', self.F1) + # print ('H1', self.H1) + # print ('F2', self.F2) + # print ('normfacs', self.normfacs) #------------------------------------- def append_net(self, model, label, F): @@ -264,6 +270,7 @@ def loadNormalization(self, filename): def loadAndApplyNormalization(self, filename): self.loadNormalization(filename) idx = 0 + #for c in self.COL_U : for c in self.COL_U + self.COL_Y: dmean = self.normfacs[c]['offset'] drange = self.normfacs[c]['scale'] @@ -282,8 +289,16 @@ def initializeModelStructure(self): self.y0 = torch.zeros((self.batch_size, self.na), dtype=torch.float) self.u0 = torch.zeros((self.batch_size, self.nb), dtype=torch.float) self.F2 = MimoStaticNonLinearity(in_channels=len(idx_out), out_channels=len(idx_out), n_hidden=self.nh2, activation=self.activation) - - def trainModelCoefficients(self): + self.H2 = MimoLinearDynamicalOperator(in_channels=len(idx_out2), out_channels=len(idx_out3), n_b=3, n_a=1, n_k=0) + + H2a_coeff = torch.tensor([[[0.00],[0.00 ]]]) + self.H2.a_coeff = nn.Parameter(H2a_coeff) + H2b_coeff = torch.tensor([[[-0.955e-2,0.955e-2,0],[1.0095,-1.91e-2,0.955e-2]]]) + self.H2.b_coeff = nn.Parameter(H2b_coeff) + print ('H2.a', self.H2.a_coeff) + print ('H2.b', self.H2.b_coeff) + + def trainModelCoefficients(self, bMAE = False): self.optimizer = torch.optim.Adam([ {'params': self.F1.parameters(), 'lr': self.lr}, {'params': self.H1.parameters(), 'lr': self.lr}, @@ -303,25 +318,62 @@ def trainModelCoefficients(self): y_non = self.F1 (ub) y_lin = self.H1 (y_non, self.y0, self.u0) y_hat = self.F2 (y_lin) + y_hat1 = y_hat[:,:,0:2] + y_hat2 = y_hat[:,:,2:4] + + # print('ub', ub.size()) + # print('y_non', y_non.size()) + # print('y_lin', y_lin.size()) + # print('y_hat', y_hat.size()) + # print('y_hat1', y_hat1.size()) + # print('y_hat2', y_hat2.size()) + + + y_lin2 = self.H2 (y_hat2, self.y0, self.u0) + #print('y_lin2', y_lin2.size()) + #print('yb', yb.size()) + y_hat3 = torch.cat((y_hat1,y_lin2), dim = 2) + # print('y_hat3', y_hat2.size()) # Compute fit loss - err_fit = yb - y_hat - loss_fit = torch.sum(torch.abs(err_fit)) + err_fit0 = yb[:,:,0]*1 - y_hat3[:,:,0]*1 + err_fit1 = yb[:,:,1]*1 - y_hat3[:,:,1]*1 + err_fit2 = yb[:,:,2]/2 - y_hat3[:,:,2]/2 + + err_fit = torch.cat((err_fit0,err_fit1,err_fit2)) + if bMAE: + loss_fit = torch.sum(torch.abs(err_fit)) + #print('mae error') + else: + loss_fit = torch.mean(err_fit**2) + + + loss = loss_fit + LOSS.append(loss.item()) if itr % self.print_freq == 0: - print('Iter {:4d} of {:4d} | Loss {:12.4f}'.format (itr, self.num_iter, loss_fit)) + print('Iter {:4d} of {:4d} | Loss {:12.4f}'.format (itr, self.num_iter, loss)) # Optimize loss.backward() self.optimizer.step() train_time = time.time() - start_time + + # print ('F1', self.F1) + print ('H2.a', self.H2.a_coeff) + print ('H2.b', self.H2.b_coeff) + # print ('H1.out', self.H1.out_channels) + # print ('F2', self.F2) + return train_time, LOSS def saveModelCoefficients(self, model_folder): torch.save(self.F1.state_dict(), os.path.join(model_folder, "F1.pkl")) torch.save(self.H1.state_dict(), os.path.join(model_folder, "H1.pkl")) torch.save(self.F2.state_dict(), os.path.join(model_folder, "F2.pkl")) + torch.save(self.H2.state_dict(), os.path.join(model_folder, "H2.pkl")) + def loadModelCoefficients(self, model_folder): B1 = torch.load(os.path.join(model_folder, "F1.pkl")) @@ -330,6 +382,10 @@ def loadModelCoefficients(self, model_folder): self.H1.load_state_dict(B2) B3 = torch.load(os.path.join(model_folder, "F2.pkl")) self.F2.load_state_dict(B3) + B4 = torch.load(os.path.join(model_folder, "H2.pkl")) + self.H2.load_state_dict(B4) + + def exportModel(self, filename): config = {'name':'PV1', 'type':'F1+H1+F2', 't_step': self.t_step} @@ -372,6 +428,8 @@ def printStateDicts(self): print ('F2', self.F2.state_dict()) print ('H1', self.H1.in_channels, self.H1.out_channels, self.H1.n_a, self.H1.n_b, self.H1.n_k) print (self.H1.state_dict()) + print ('H2', self.H2.in_channels, self.H2.out_channels, self.H2.n_a, self.H2.n_b, self.H2.n_k) + print (self.H2.state_dict()) def testOneCase(self, case_idx): case_data = self.data_train[[case_idx],:,:] @@ -379,14 +437,22 @@ def testOneCase(self, case_idx): y_non = self.F1 (ub) y_lin = self.H1 (y_non, self.y0, self.u0) y_hat = self.F2 (y_lin) + y_hat1 = y_hat[:,:,0:2] + y_hat2 = y_hat[:,:,2:4] + y_lin2 = self.H2 (y_hat2, self.y0, self.u0) + y_hat3 = torch.cat((y_hat1,y_lin2), dim = 2) + print (ub.shape, y_non.shape, y_lin.shape, y_hat.shape) + print ('H2', self.H2.in_channels, self.H2.out_channels, self.H2.n_a, self.H2.n_b, self.H2.n_k) + print (self.H2.state_dict()) # self.printStateDicts() # print (y_lin) - y_hat = y_hat.detach().numpy()[[0], :, :] - y_true = np.transpose(case_data[0,:,idx_out]) - rmse = dynonet.metrics.error_rmse(y_true, y_hat[0]) - return rmse, y_hat, y_true, np.transpose(case_data[0,:,idx_in]) + y_hat3 = y_hat3.detach().numpy()[[0], :, :] + y_hat2 = y_hat2.detach().numpy()[[0], :, :] + y_true = np.transpose(case_data[0,:,idx_out[:3]]) + rmse = dynonet.metrics.error_rmse(y_true, y_hat3[0]) + return rmse, y_hat3, y_true, np.transpose(case_data[0,:,idx_in]), y_hat2 def stepOneCase(self, case_idx): case_data = self.data_train[case_idx,:,:] @@ -419,9 +485,13 @@ def trainingErrors(self, bByCase=False): y_non = self.F1 (ub) y_lin = self.H1 (y_non, self.y0, self.u0) y_hat = self.F2 (y_lin) + y_hat1 = y_hat[:,:,0:2] + y_hat2 = y_hat[:,:,2:4] + y_lin2 = self.H2 (y_hat2, self.y0, self.u0) + y_hat3 = torch.cat((y_hat1,y_lin2), dim = 2) - y_hat = y_hat.detach().numpy() - y_true = self.data_train[:,:,idx_out] + y_hat3 = y_hat3.detach().numpy() + y_true = self.data_train[:,:,idx_out[:3]] self.n_cases = self.data_train.shape[0] icol = 0 @@ -429,18 +499,21 @@ def trainingErrors(self, bByCase=False): total_mae = {} if bByCase: case_rmse = lst1 = [dict() for i in range(self.n_cases)] + case_mae = lst2 = [dict() for i in range(self.n_cases)] else: case_rmse = None + case_mae = None for col in self.COL_Y: SUMSQ = 0.0 MAE = 0.0 for icase in range(self.n_cases): y1 = y_true[icase,:,icol] - y2 = y_hat[icase,:,icol] + y2 = y_hat3[icase,:,icol] colmae = dynonet.metrics.error_mae(y1, y2) colrms = dynonet.metrics.error_rmse(y1, y2) if bByCase: case_rmse[icase][col] = colrms + case_mae[icase][col] = colmae MAE += colmae SUMSQ += (colrms*colrms) MAE /= self.n_cases @@ -449,7 +522,7 @@ def trainingErrors(self, bByCase=False): total_mae[col] = MAE total_rmse[col] = RMSE icol += 1 - return total_rmse, total_mae, case_rmse + return total_rmse, total_mae, case_rmse, case_mae def set_LCL_filter(self, Lf, Cf, Lc): self.Lf = Lf diff --git a/examples/hwpv/pv1_test.py b/examples/hwpv/pv1_test.py index d4bf769..f724b06 100644 --- a/examples/hwpv/pv1_test.py +++ b/examples/hwpv/pv1_test.py @@ -7,18 +7,27 @@ #import pv1_fhf as pv1_model #import pv1_model #import pv1_feedback as pv1_model +os.environ["KMP_DUPLICATE_LIB_OK"]="TRUE" data_path = r'./data/pv1.hdf5' model_folder = r'./models' def plot_case(model, idx): - rmse, y_hat, y_true, u = model.testOneCase(idx) + rmse, y_hat, y_true, u, VsIs = model.testOneCase(idx) # rmse, y_hat, y_true, u = model.stepOneCase(idx) print ('column', model.COL_Y, 'RMS errors', rmse) valstr = ' '.join('{:.4f}'.format(rms) for rms in rmse) # print ('y_hat shape', y_hat.shape) # print ('y_true shape', y_true.shape) # print ('u shape', u.shape) + + fig, ax = plt.subplots (1, 2, sharex = 'col', figsize=(7,4), constrained_layout=True) + ax[0].plot (model.t, VsIs[0,:,0]*290.83447265625 + 203.58287048339844, label='Vs') + ax[0].legend() + ax[1].plot (model.t, VsIs[0,:,1]*66.92162322998047 + 30.236408233642578, label='Is') + ax[1].legend() + + plt.show() fig, ax = plt.subplots (2, 5, sharex = 'col', figsize=(15,8), constrained_layout=True) fig.suptitle ('Case {:d} Simulation; Output RMSE = {:s}'.format(idx, valstr)) @@ -43,6 +52,10 @@ def plot_case(model, idx): if bNormalized: scale = 1.0 offset = 0.0 + print('offset',offset) + print('scale',scale) + #scale = 1.0 + #offset = 0.0 ax[1,j].set_title ('Output {:s}'.format (col)) ax[1,j].plot (model.t, y_true[:,j]*scale + offset, label='y') ax[1,j].plot (model.t, y_hat[0,:,j]*scale + offset, label='y_hat') @@ -50,6 +63,21 @@ def plot_case(model, idx): ax[1,j].legend() j += 1 plt.show() + + plt.figure() + #plt.plot (model.t, y_true[:,2], label='Irms_true') + plt.plot (model.t, y_hat[0,:,2], label='Irms_hat') + plt.plot (model.t, VsIs[0,:,1], label='Is') + plt.legend() + plt.show() + + plt.figure() + #plt.plot (model.t, y_true[:,2], label='Irms_true') + #plt.plot (model.t, y_hat[0,:,2], label='Irms_hat') + plt.plot (model.t, VsIs[0,:,1]-y_hat[0,:,2], label='Is - Irms') + plt.legend() + plt.show() + if __name__ == '__main__': @@ -72,4 +100,4 @@ def plot_case(model, idx): plot_case (model, idx) else: plot_case (model, case_idx) - + \ No newline at end of file diff --git a/examples/hwpv/pv1_training.py b/examples/hwpv/pv1_training.py index b9a485a..49f9d98 100644 --- a/examples/hwpv/pv1_training.py +++ b/examples/hwpv/pv1_training.py @@ -5,6 +5,7 @@ #import pv1_fhf as pv1_model #import pv1_model #import pv1_feedback as pv1_model +os.environ["KMP_DUPLICATE_LIB_OK"]="TRUE" data_path = r'./data/pv1.hdf5' model_folder = r'./models' @@ -15,18 +16,20 @@ model.loadTrainingData(data_path) model.applyAndSaveNormalization(model_folder) model.initializeModelStructure() - train_time, LOSS = model.trainModelCoefficients() + train_time, LOSS = model.trainModelCoefficients(bMAE=False) model.saveModelCoefficients(model_folder) # quit() - rmse, mae, case_rmse = model.trainingErrors(False) + rmse, mae, case_rmse, case_mae = model.trainingErrors(False) nlookback = 10 * int(model.n_cases / model.batch_size) recent_loss = LOSS[len(LOSS)-nlookback:] # print (nlookback, recent_loss) - print ('COL_Y', model.COL_Y) +# print ('COL_Y', model.COL_Y) valstr = ' '.join('{:.4f}'.format(rmse[col]) for col in model.COL_Y) print ('Train time: {:.2f}, Recent loss: {:.2f}, RMS Errors: {:s}'.format (train_time, np.mean(recent_loss), valstr)) + valstr = ' '.join('{:.4f}'.format(mae[col]) for col in model.COL_Y) + print (' MAE Errors: {:s}'.format (valstr)) plt.figure() plt.plot(LOSS)