update_threshold

Nidelva3D/src/Agents/Agent_adaptive.py

@@ -59,6 +59,8 @@ def __init__(self, kernel: str = "GMRF", num_steps: int = 5,
         self.__ibv = np.zeros([self.__num_steps, 1])
         self.__vr = np.zeros([self.__num_steps, 1])
         self.__rmse = np.zeros([self.__num_steps, 1])
+        self.__ce = np.zeros([self.__num_steps, 1])
+        self.__corr_coef = np.zeros([self.__num_steps, 1])
         self.update_metrics()
 
         self.debug = debug
@@ -186,6 +188,35 @@ def update_metrics(self) -> None:
 
         mu[mu < 0] = 0
         self.__rmse[self.__counter] = mean_squared_error(self.mu_truth.flatten(), mu.flatten(), squared=False)
+        self.__corr_coef[self.__counter] = self.__cal_corr_coef(self.mu_truth.flatten(), mu.flatten())
+        self.__ce[self.__counter] = self.__cal_classification_error(self.mu_truth.flatten(),
+                                                                    mu.flatten(), self.__threshold)
+
+    @staticmethod
+    def __cal_corr_coef(x, y) -> float:
+        """
+        Calculate the correlation coefficient between two vectors.
+
+        Method:
+        1. Calculate the covariance matrix between two vectors.
+        2. Calculate the correlation coefficient.
+
+        """
+        d1 = x - np.mean(x)
+        d2 = y - np.mean(y)
+        cov = d1.T @ d2
+        corr = cov / np.linalg.norm(d1) / np.linalg.norm(d2)
+        return corr
+
+    @staticmethod
+    def __cal_classification_error(x, y, threshold) -> float:
+        """
+        Calculate the classification error between two vectors.
+        """
+        X = np.where(x > threshold, 1, 0)
+        Y = np.where(y > threshold, 1, 0)
+        CE = np.sum(np.abs(X - Y)) / len(X)
+        return CE
 
     @staticmethod
     def __get_ibv(threshold: float, mu: np.ndarray, sigma_diag: np.ndarray) -> np.ndarray:
@@ -204,7 +235,7 @@ def get_counter(self):
         return self.__counter
 
     def get_metrics(self) -> tuple:
-        return self.__ibv, self.__vr, self.__rmse
+        return self.__ibv, self.__vr, self.__rmse, self.__corr_coef, self.__ce
 
 
 if __name__ == "__main__":

Nidelva3D/src/GMRF/spde.py

@@ -122,7 +122,7 @@ def setThreshold(self):
         """Set threshold for Excursion set
         """
         ind = np.load(FILEPATH + 'models/boundary.npy')
-        self.threshold = 23 # self.mu[ind].mean()
+        self.threshold = self.mu[ind].mean()
         print('Treshold is set to %.2f'%(self.threshold))
 
     def getThreshold(self) -> np.ndarray:

Nidelva3D/src/Simulators/Simulator.py

@@ -29,10 +29,12 @@ def __init__(self, num_steps: int = 5, random_seed: int = 0, temporal_truth: boo
     def run(self) -> 'pd.DataFrame':
         self.grf_agent.run()
         self.gmrf_agent.run()
-        self.ibv_grf, self.vr_grf, self.rmse_grf = self.grf_agent.get_metrics()
-        self.ibv_gmrf, self.vr_gmrf, self.rmse_gmrf = self.gmrf_agent.get_metrics()
-        df = np.hstack((self.ibv_grf, self.vr_grf, self.rmse_grf, self.ibv_gmrf, self.vr_gmrf, self.rmse_gmrf))
-        df = pd.DataFrame(df, columns=["ibv_grf", "vr_grf", "rmse_grf", "ibv_gmrf", "vr_gmrf", "rmse_gmrf"])
+        self.ibv_grf, self.vr_grf, self.rmse_grf, self.corr_grf, self.ce_grf = self.grf_agent.get_metrics()
+        self.ibv_gmrf, self.vr_gmrf, self.rmse_gmrf, self.corr_gmrf, self.ce_gmrf = self.gmrf_agent.get_metrics()
+        df = np.hstack((self.ibv_grf, self.vr_grf, self.rmse_grf, self.corr_grf, self.ce_grf,
+                        self.ibv_gmrf, self.vr_gmrf, self.rmse_gmrf, self.corr_gmrf, self.ce_gmrf))
+        df = pd.DataFrame(df, columns=["ibv_grf", "vr_grf", "rmse_grf", "corr_grf", "ce_grf",
+                                       "ibv_gmrf", "vr_gmrf", "rmse_gmrf", "corr_gmrf", "ce_gmrf"])
         return df
 
 

Nidelva3D/src/tests/test_simulator.py

@@ -14,8 +14,8 @@
 class TestSimulator(TestCase):
 
     def setUp(self) -> None:
-        num_steps = 10
-        random_seed = 18
+        num_steps = 20
+        random_seed = 19
         debug = True
         temporal_truth = True
         self.simulator = Simulator(num_steps=num_steps, random_seed=random_seed,
@@ -26,29 +26,41 @@ def test_run_simulator(self) -> None:
         ibv = np.stack((df["ibv_grf"], df["ibv_gmrf"]), axis=1)
         vr = np.stack((df["vr_grf"], df["vr_gmrf"]), axis=1)
         rmse = np.stack((df["rmse_grf"], df["rmse_gmrf"]), axis=1)
+        corr = np.stack((df["corr_grf"], df["corr_gmrf"]), axis=1)
+        ce = np.stack((df["ce_grf"], df["ce_gmrf"]), axis=1)
 
-        plt.figure(figsize=(25, 8))
-        plt.subplot(131)
+        plt.figure(figsize=(40, 8))
+        plt.subplot(151)
         plt.plot(ibv[:, 0], label="GRF")
         plt.plot(ibv[:, 1], label="GMRF")
         # plt.plot(ibv[:, 0]/np.amax(ibv[:, 0]), label="GRF")
         # plt.plot(ibv[:, 1]/np.amax(ibv[:, 1]), label="GMRF")
         plt.legend()
         plt.title("IBV")
-        plt.subplot(132)
+        plt.subplot(152)
         plt.plot(vr[:, 0], label="GRF")
         plt.plot(vr[:, 1], label="GMRF")
         # plt.plot(vr[:, 0]/np.amax(vr[:, 0]), label="GRF")
         # plt.plot(vr[:, 1]/np.amax(vr[:, 1]), label="GMRF")
         plt.legend()
         plt.title("VR")
-        plt.subplot(133)
+        plt.subplot(153)
         plt.plot(rmse[:, 0], label="GRF")
         plt.plot(rmse[:, 1], label="GMRF")
         # plt.plot(rmse[:, 0]/np.amax(rmse[:, 0]), label="GRF")
         # plt.plot(rmse[:, 1]/np.amax(rmse[:, 1]), label="GMRF")
         plt.legend()
         plt.title("RMSE")
+        plt.subplot(154)
+        plt.plot(corr[:, 0], label="GRF")
+        plt.plot(corr[:, 1], label="GMRF")
+        plt.legend()
+        plt.title("Correlation")
+        plt.subplot(155)
+        plt.plot(ce[:, 0], label="GRF")
+        plt.plot(ce[:, 1], label="GMRF")
+        plt.legend()
+        plt.title("CE")
         plt.show()
         print("stop")