Add support for inference with uncorrelated and correlated gaussianb …

…noise
JBjoernskov · Jan 9, 2024 · 71fc2b7 · 71fc2b7
1 parent 99fb356
commit 71fc2b7
Show file tree

Hide file tree

Showing 2 changed files with 17 additions and 13 deletions.
diff --git a/twin4build/estimator/tests/test_load_emcee_chain.py b/twin4build/estimator/tests/test_load_emcee_chain.py
@@ -48,13 +48,13 @@ def test_load_emcee_chain():
     sky_blue = colors[9]
     # plot.load_params()
 
-    do_iac_plot = True
-    do_logl_plot = True
-    do_trace_plot = True
-    do_swap_plot = True
-    do_jump_plot = True
-    do_corner_plot = True
-    do_inference = False
+    do_iac_plot = False
+    do_logl_plot = False
+    do_trace_plot = False
+    do_swap_plot = False
+    do_jump_plot = False
+    do_corner_plot = False
+    do_inference = True
 
     # loaddir = os.path.join(uppath(os.path.abspath(__file__), 2), "chain_logs", "20230829_155706_chain_log.pickle")
     # loaddir = os.path.join(uppath(os.path.abspath(__file__), 2), "chain_logs", "20230830_194210_chain_log.pickle")
@@ -121,7 +121,7 @@ def test_load_emcee_chain():
         result["chain.T"] = 1/result["chain.betas"] ##################################
 
 
-
+    burnin = int(result["chain.x"].shape[0])-500 #800
 #########################################
     list_ = ["integratedAutoCorrelatedTime", "chain.jumps_accepted", "chain.jumps_proposed", "chain.swaps_accepted", "chain.swaps_proposed"]
     for key in list_:
@@ -216,11 +216,11 @@ def test_load_emcee_chain():
         model.chain_log["n_par"] = n_par
         model.chain_log["n_par_map"] = n_par_map
         parameter_chain = result["chain.x"][burnin:,0,:,:]
-        del result
-        del model.chain_log["chain.x"]
+        # del result
+        # del model.chain_log["chain.x"]
         # parameter_chain = result["chain.x"][-1:,0,:,:] #[-1:,0,:,:]
         parameter_chain = parameter_chain.reshape((parameter_chain.shape[0]*parameter_chain.shape[1], parameter_chain.shape[2]))
-        fig, axes = simulator.run_emcee_inference(model, parameter_chain, targetParameters, targetMeasuringDevices, startTime, endTime, stepSize)
+        fig, axes = simulator.run_emcee_inference(model, parameter_chain, targetParameters, targetMeasuringDevices, startTime, endTime, stepSize, assume_uncorrelated_noise=False)
         ylabels = [r"$u_v [1]$", r"$T_{c,w,in} [^\circ\!C]$", r"$T_{c,w,out} [^\circ\!C]$", r"$T_{c,a,out} [^\circ\!C]$", r"$\dot{P}_f [W]$"]
         fig.subplots_adjust(hspace=0.3)
         fig.set_size_inches((15,10))
@@ -263,7 +263,7 @@ def test_load_emcee_chain():
     nrows = math.ceil(nparam/ncols)
 
 
-    burnin = int(result["chain.x"].shape[0])-500 #800
+
     # cm = plt.get_cmap('RdYlBu', ntemps)
     # cm_sb = sns.color_palette("vlag_r", n_colors=ntemps, center="dark") #vlag_r
     cm_sb = sns.diverging_palette(210, 0, s=50, l=50, n=ntemps, center="dark") #vlag_r

diff --git a/twin4build/simulator/simulator.py b/twin4build/simulator/simulator.py
@@ -314,11 +314,15 @@ def _sim_func_gaussian_process(self, model, theta, stepSize, startTime, endTime,
                             show_progress_bar=False)
             simulation_readings_train = np.zeros((len(self.dateTimeSteps), len(self.targetMeasuringDevices)))
             actual_readings_train = np.zeros((len(self.dateTimeSteps), len(self.targetMeasuringDevices)))
+            x_train = {}
             # t_train = self.secondTimeSteps
             # x_train = self.get_actual_readings(startTime=model.chain_log["startTime_train"], endTime=model.chain_log["endTime_train"], stepSize=model.chain_log["stepSize_train"], reading_type="input").to_numpy()
             for j, measuring_device in enumerate(self.targetMeasuringDevices):
                 simulation_readings_train[:,j] = np.array(next(iter(measuring_device.savedInput.values())))
                 actual_readings_train[:,j] = df_actual_readings_train[measuring_device.id].to_numpy()
+                source_component = [cp.connectsSystemThrough.connectsSystem for cp in measuring_device.connectsAt][0]
+                x = np.array(list(source_component.savedInput.values())).transpose()
+                x_train[measuring_device.id] = x
 
             self.simulate(model,
                             stepSize=stepSize,
@@ -364,7 +368,7 @@ def _sim_func_gaussian_process(self, model, theta, stepSize, startTime, endTime,
                 # kernel = kernels.Matern32Kernel(metric=scale_lengths, ndim=scale_lengths.size)
                 kernel = kernels.ExpSquaredKernel(metric=scale_lengths, ndim=scale_lengths.size)
                 gp = george.GP(a*kernel)
-                gp.compute(x, standardDeviation[j])
+                gp.compute(x_train[measuring_device.id], standardDeviation[j])
                 y_model[:,j] = simulation_readings
                 y[:,j] = np.mean(gp.sample_conditional(actual_readings_train[:,j]-simulation_readings_train[:,j], x, n_samples), axis=0) + simulation_readings
                 n_prev = n