update_simuilation_design

Publication/src/AUVSimulator/AUVSimulator.py

@@ -28,6 +28,7 @@ class AUVSimulator:
     def __init__(self, random_seed: int = 0, sigma: float = 1.,
                  loc_start: np.ndarray = np.array([0, 0]), temporal_truth: bool = True) -> None:
         self.temporal_truth = temporal_truth
+        print("Temporal truth is on.") if self.temporal_truth else print("Temporal truth is off.")
         self.ctd = CTDSimulator(random_seed=random_seed, sigma=sigma)
         self.messenger = Messenger()
         self.__loc = loc_start
@@ -110,7 +111,6 @@ def __collect_data(self) -> None:
             y_path = np.linspace(y_start, y_end, N)
             loc = np.stack((x_path, y_path), axis=1)
             if self.temporal_truth:
-                print("Temporal truth is on.")
                 sal = self.ctd.get_salinity_at_dt_loc(dt, loc)
             else:
                 raise NotImplementedError

Publication/src/AUVSimulator/CTDSimulator.py

@@ -25,6 +25,7 @@ def __init__(self, random_seed: int = 0,
         """
         Set up the CTD simulated truth field.
         """
+        # Load SINMOD data from a specific file
         np.random.seed(random_seed)
         filepath_sinmod = filepath
         datestring = filepath_sinmod.split("/")[-1].split("_")[-1][:-3].replace('.', '-') + " 10:00:00"
@@ -53,70 +54,27 @@ def __init__(self, random_seed: int = 0,
         self.L = np.linalg.cholesky(cov)
         print("Cholesky decomposition takes: ", time() - t0)
 
+        self.mu_truth = np.empty([0, len(self.grid_sinmod)])
         self.construct_ground_truth_field()
 
     def construct_ground_truth_field(self) -> None:
         """
         This function constructs the ground truth field given all the timestamps and locations.
         """
-        print("Start constructing ground truth field!")
         x0 = (self.salinity_sinmod[0, :, :].flatten() +
              (self.L @ np.random.randn(len(self.L)).reshape(-1, 1)).flatten())
         xt_1 = x0
         xt = x0
-        self.mu_truth = np.empty([0, len(xt)])
         self.mu_truth = np.vstack((self.mu_truth, xt))
 
-        import matplotlib.pyplot as plt
-        from matplotlib.pyplot import get_cmap
-        from matplotlib.gridspec import GridSpec
-        figpath = os.getcwd() + "/../../../../OneDrive - NTNU/MASCOT_PhD/Projects" \
-                                "/GOOGLE/Docs/fig/Sim_2DNidelva/Simulator/groundtruth/"
         t0 = time()
         for i in range(1, len(self.timestamp_sinmod)):
-            print("Time string: ", datetime.fromtimestamp(self.timestamp_sinmod[i]).strftime("%Y-%m-%d %H:%M:%S"))
-            fig = plt.figure(figsize=(35, 10))
-            gs = GridSpec(1, 3, figure=fig)
-            ax = fig.add_subplot(gs[0])
-            im = ax.scatter(self.grid_sinmod[:, 1], self.grid_sinmod[:, 0], c=xt, cmap=get_cmap("BrBG", 10),
-                        vmin=10, vmax=30)
-            plt.colorbar(im)
-            ax.set_title("Ground truth at time: " +
-                         datetime.fromtimestamp(self.timestamp_sinmod[i]).strftime("%Y-%m-%d %H:%M:%S"))
-            ax.set_xlabel("East (m)")
-            ax.set_ylabel("North (m)")
-
-            ax = fig.add_subplot(gs[1])
-            im = ax.scatter(self.grid_sinmod[:, 1], self.grid_sinmod[:, 0], c=self.salinity_sinmod[i, :, :].flatten(),
-                        cmap=get_cmap("BrBG", 10), vmin=10, vmax=30)
-            plt.colorbar(im)
-            ax.set_title("SINMOD at time: " +
-                         datetime.fromtimestamp(self.timestamp_sinmod[i]).strftime("%Y-%m-%d %H:%M:%S"))
-            ax.set_xlabel("East (m)")
-            ax.set_ylabel("North (m)")
-
-            ax = fig.add_subplot(gs[2])
-            im = ax.scatter(self.grid_sinmod[:, 1], self.grid_sinmod[:, 0],
-                            c=xt - self.salinity_sinmod[i, :, :].flatten(),
-                            cmap=get_cmap("RdBu", 10), vmin=-5, vmax=5)
-            plt.colorbar(im)
-            ax.set_title("Difference at time: " +
-                            datetime.fromtimestamp(self.timestamp_sinmod[i]).strftime("%Y-%m-%d %H:%M:%S"))
-            ax.set_xlabel("East (m)")
-            ax.set_ylabel("North (m)")
-
-            plt.savefig(figpath + "P_{:03d}.png".format(i))
-            plt.close("all")
-
             xt = self.salinity_sinmod[i, :, :].flatten() + \
                  self.ar1_corr * (xt_1 - self.salinity_sinmod[i-1, :, :].flatten()) + \
                  (np.sqrt(1 - self.ar1_corr ** 2) * (self.L @ np.random.randn(len(self.L)).reshape(-1, 1))).flatten()
             xt_1 = xt
             self.mu_truth = np.vstack((self.mu_truth, xt))
-
-        self.mu_truth
-        print("Time consumed: ", time() - t0)
-        xt_1
+        print("Constructing ground truth field takes: ", time() - t0)
 
     def get_salinity_at_dt_loc(self, dt: float, loc: np.ndarray) -> Union[np.ndarray, None]:
         """
@@ -131,12 +89,11 @@ def get_salinity_at_dt_loc(self, dt: float, loc: np.ndarray) -> Union[np.ndarray
         ts = np.array([self.timestamp])
         dist, ind_time = self.timestamp_sinmod_tree.query(ts)
         t1 = time()
-        sorted_salinity = self.salinity_sinmod[ind_time, :, :].flatten()
+        sorted_salinity = self.mu_truth[ind_time, :].flatten()
+        # sorted_salinity = self.salinity_sinmod[ind_time, :, :].flatten()
         dist, ind_loc = self.grid_sinmod_tree.query(loc)
         print("Query salinity at timestamp and location takes: ", time() - t1)
-
-
-        return sorted_salinity[ind_loc] + (L @ np.random.randn(len(L)).reshape(-1, 1)).flatten()
+        return sorted_salinity[ind_loc]
 
 
 if __name__ == "__main__":

Publication/src/GRF/GRF.py

@@ -9,6 +9,7 @@
 Date: 2023-08-22
 """
 from Field import Field
+from SINMOD import SINMOD
 from usr_func.vectorize import vectorize
 from usr_func.checkfolder import checkfolder
 from usr_func.normalize import normalize
@@ -21,14 +22,16 @@
 from pykdtree.kdtree import KDTree
 import time
 import pandas as pd
+import os
 
 
 class GRF:
     """
     GRF kernel
     """
     def __init__(self, sigma: float = 1., nugget: float = .4,
-                 approximate_eibv: bool = False, fast_eibv: bool = True) -> None:
+                 approximate_eibv: bool = False, fast_eibv: bool = True,
+                 filepath_prior: str = os.getcwd() + "/../sinmod/samples_2022.05.10.nc") -> None:
         """ Initializes the parameters in GRF kernel. """
         self.__approximate_eibv = approximate_eibv
         self.__fast_eibv = fast_eibv
@@ -78,6 +81,8 @@ def __init__(self, sigma: float = 1., nugget: float = .4,
         self.__construct_grf_field()
 
         # s1: update prior mean
+        self.__sinmod = SINMOD(filepath_prior)
+        self.__salinity_sinmod = self.__sinmod.get_salinity()[:, 0, :, :]
         self.__construct_prior_mean()
         self.__mu_prior = self.__mu
         self.__Sigma_prior = self.__Sigma
@@ -92,15 +97,14 @@ def __construct_grf_field(self) -> None:
                                             np.exp(-self.__eta * self.__distance_matrix))
 
     def __construct_prior_mean(self) -> None:
-        # s0: get delft3d dataset
-        dataset_sinmod = pd.read_csv("./../prior/sinmod.csv").to_numpy()
-        grid_sinmod = dataset_sinmod[:, :2]
-        sal_sinmod = dataset_sinmod[:, -1]
+        # dataset_sinmod = pd.read_csv("./../prior/sinmod.csv").to_numpy()
+        # grid_sinmod = dataset_sinmod[:, :2]
+        # sal_sinmod = dataset_sinmod[:, -1]
 
         # s1: interpolate onto grid.
-        dm_grid_sinmod = cdist(self.grid, grid_sinmod)
-        ind_close = np.argmin(dm_grid_sinmod, axis=1)
-        self.__mu = vectorize(sal_sinmod[ind_close])
+        # dm_grid_sinmod = cdist(self.grid, grid_sinmod)
+        # ind_close = np.argmin(dm_grid_sinmod, axis=1)
+        # self.__mu = vectorize(sal_sinmod[ind_close])
 
     def __load_cdf_table(self) -> None:
         """

Publication/src/tests/test_auv_simulator.py

@@ -23,7 +23,7 @@ def value(x, y):
 class TestAUVSimulator(TestCase):
 
     def setUp(self) -> None:
-        self.auv = AUVSimulator()
+        self.auv = AUVSimulator(random_seed=0, sigma=1., loc_start=np.array([0, 0]), temporal_truth=True)
         self.polygon_border_wgs = pd.read_csv(os.getcwd() + "/csv/polygon_border.csv").to_numpy()
         x, y = WGS.latlon2xy(self.polygon_border_wgs[:, 0], self.polygon_border_wgs[:, 1])
         self.polygon_border = np.stack((x, y), axis=1)

Publication/src/tests/test_ctd_simulator.py

@@ -1,4 +1,8 @@
 """ Unit test for CTD simulator
+
+Author: Yaolin Ge
+Email: [email protected]
+Date: 2023-08-30
 """
 
 from unittest import TestCase
@@ -17,14 +21,23 @@
 class TestCTDSimulator(TestCase):
 
     def setUp(self) -> None:
-        self.ctd = CTDSimulator(random_seed=13)
         filepath_sinmod = os.getcwd() + "/../sinmod/samples_2022.05.11.nc"
+        self.ctd = CTDSimulator(random_seed=14, filepath=filepath_sinmod, sigma=1.)
         self.sinmod = SINMOD(filepath_sinmod)
 
     def test_get_salinity_at_time_loc(self) -> None:
         # c1, whole region
-        # x, y, depth = self.sinmod.get_coordinates()
-        # loc = np.vstack((x.flatten(), y.flatten())).T
+        x, y, depth = self.sinmod.get_coordinates()
+        loc = np.vstack((x.flatten(), y.flatten())).T
+
+        dt = 1200
+        for i in range(10):
+            salinity = self.ctd.get_salinity_at_dt_loc(dt, loc)
+            plt.figure()
+            plt.scatter(loc[:, 1], loc[:, 0], c=salinity, cmap=get_cmap("BrBG", 10), vmin=10, vmax=30)
+            plt.colorbar()
+            plt.title("Time: " + datetime.fromtimestamp(self.ctd.timestamp).strftime("%Y-%m-%d %H:%M:%S"))
+            plt.show()
 
         # c2, regional locations
         xx = np.linspace(2000, 4000, 100)