Merge pull request #180 from slacgismo/dilation-dev

Development of dilation module. Fixes a bug on NaN handling and provides options for dilating both raw and filled matrices. Bypassing branch protection because tests all passed.
slacgismo · Jan 22, 2025 · d5aa80d · d5aa80d
2 parents 6055471 + d6a8d4b
commit d5aa80d
Show file tree

Hide file tree

Showing 3 changed files with 332 additions and 228 deletions.
diff --git a/notebooks/archive/Signal_dilatation.ipynb b/notebooks/archive/Signal_dilatation.ipynb
diff --git a/notebooks/archive/Signal_dilation.ipynb b/notebooks/archive/Signal_dilation.ipynb
diff --git a/solardatatools/algorithms/dilation.py b/solardatatools/algorithms/dilation.py
@@ -10,22 +10,28 @@
 
 DEFAULT = {
     "nvals_dil": 101,
+    "matrix": "raw",
 }
 
 
 class Dilation:
     def __init__(self, data_handler, **config):
-        self.dh = data_handler
-        self.nvals_ori = data_handler.raw_data_matrix.shape[0]
-        self.ndays = data_handler.raw_data_matrix.shape[1]
-        self.idx_ori = None
-        self.idx_dil = None
-        self.signal_ori = data_handler.raw_data_matrix.ravel(order="F")
-        self.signal_dil = None
         if len(config) == 0:
             self.config = DEFAULT
         else:
             self.config = config
+        if self.config["matrix"] == "raw":
+            mat = data_handler.raw_data_matrix
+        elif self.config["matrix"] == "filled":
+            mat = data_handler.filled_data_matrix
+        else:
+            raise ValueError("Invalid value for matrix. Choose from: ['raw', 'filled']")
+        self.dh = data_handler
+        self.nvals_ori, self.ndays = mat.shape
+        self.idx_ori = None
+        self.idx_dil = None
+        self.signal_ori = mat.ravel(order="F")
+        self.signal_dil = None
         self.nvals_dil = self.config["nvals_dil"]
         self.run()
 
@@ -278,17 +284,18 @@ def compute_integral_interpolation(ttnew, xxnew, new_indices):
     """
     # replace NaNs with zeros
     xxnew_filled = np.nan_to_num(xxnew, nan=0)
-    # compute the piecewise than cumulative integral of the signal
     piecewise_integrals = np.diff(ttnew) * xxnew_filled[:-1]
     cumulative_integrals = np.zeros(ttnew.shape[0])
     # Add the initial zero as a baseline for the cumulative sum
     cumulative_integrals[1:] = np.cumsum(piecewise_integrals)
+    # the new value at each new time point is the difference between the cumulative integrals
+    # at this point and the following one
+    y = np.diff(cumulative_integrals[new_indices])
     # set NaNs back to Na
     was_nan = np.isnan(xxnew)
     # the last point is not used in the interpolation, it shouldn't propagate NaNs
     was_nan[-1] = False
-    cumulative_integrals[was_nan] = np.nan
-    # the new value at each new time point is the difference between the cumulative integrals
-    # at this point and the following one
-    y = np.diff(cumulative_integrals[new_indices])
+    subarrays = np.split(was_nan, new_indices)
+    nan_flags = np.array([np.any(subarray) for subarray in subarrays])[1:-1]
+    y[nan_flags] = np.nan
     return y