Merge pull request #69 from brown-ccv/fix-step2-update-generalized_FT

chore: clean up generalized_FT

.github/workflows/test-B01_SL_load_single_file.yml

@@ -37,3 +37,4 @@ jobs:
         run: python src/icesat2_tracks/analysis_db/B05_define_angle.py SH_20190502_05180312 SH_testSLsinglefile2 True
       - name: Seventh step B06_correct_separate_var
         run: python src/icesat2_tracks/analysis_db/B06_correct_separate_var.py SH_20190502_05180312 SH_testSLsinglefile2 True
+

src/icesat2_tracks/ICEsat2_SI_tools/generalized_FT.py

@@ -1,8 +1,16 @@
-import numpy as np
+import copy
+import time
 
-import icesat2_tracks.ICEsat2_SI_tools.spectral_estimates as spec
-import icesat2_tracks.ICEsat2_SI_tools.lanczos as lanczos
+from numpy import linalg
+import numpy as np
+import pandas as pd
+import xarray as xr
 import matplotlib.pyplot as plt
+from scipy.signal import detrend
+import lmfit as LM
+
+from icesat2_tracks.ICEsat2_SI_tools import lanczos, spectral_estimates as spec
+import icesat2_tracks.local_modules.JONSWAP_gamma as spectal_models
 
 
 def rebin(data, dk):
@@ -69,15 +77,15 @@ def get_weights_from_data(
     pars = Spec_fft.set_parameters(flim=np.sqrt(9.81 * k[-1]) / 2 / np.pi)
     k_max = (pars["f_max"].value * 2 * np.pi) ** 2 / 9.81
 
-    if method is "gaussian":
+    if method == "gaussian":
         # simple gaussian weight
         def gaus(x, x_0, amp, sigma_g):
             return amp * np.exp(-0.5 * ((x - x_0) / sigma_g) ** 2)
 
         weight = gaus(k, k_max, 1, 0.02) ** (1 / 2)
         params = None
 
-    elif method is "parametric":
+    elif method == "parametric":
         # JONSWAP weight
         f = np.sqrt(9.81 * k) / (2 * np.pi)
         weight = Spec_fft.create_weight(freq=f, plot_flag=False, max_nfev=max_nfev)
@@ -121,8 +129,6 @@ def define_weight_shutter(weight, k, Ncut=3):
 
 
 def make_xarray_from_dict(D, name, dims, coords):
-    import xarray as xr
-
     D_return = dict()
     for xi, I in D.items():
         coords["x"] = xi
@@ -136,7 +142,7 @@ def define_weights(stancil, prior, x, y, dx, k, max_nfev, plot_flag=False):
     return weights normalized to 1, prior_pars used for the next iteration
     """
 
-    if (type(prior[0]) is bool) and not prior[0]:
+    if isinstance(prior[0], bool) and not prior[0]:
         # fit function to data
         weight, prior_pars = get_weights_from_data(
             x, y, dx, stancil, k, max_nfev, plot_flag=plot_flag, method="parametric"
@@ -157,8 +163,6 @@ def define_weights(stancil, prior, x, y, dx, k, max_nfev, plot_flag=False):
     weight = weight * define_weight_shutter(weight, k, Ncut=3)
 
     if plot_flag:
-        import matplotlib.pyplot as plt
-
         plt.plot(k, weight, zorder=12, c="darkgreen", linewidth=0.8, label=weight_name)
 
     # peak normlize weights by std of data
@@ -228,9 +232,6 @@ def cal_spectrogram(
         self.GG, params_dataframe
             params_dataframe is a pd.DataFrame that containes all the parameters of the fitting process (and may contain uncertainties too once they are calculated)
         """
-        import xarray as xr
-        import copy
-        import pandas as pd
 
         X = self.x if x is None else x  # all x positions
         DATA = self.data if data is None else data  # all data points
@@ -246,9 +247,6 @@ def calc_gFT_apply(stancil, prior):
             windows the data accoding to stencil and applies LS spectrogram
             returns: stancil center, spectrum for this stencil, number of datapoints in stancil
             """
-            from scipy.signal import detrend
-            import matplotlib.pyplot as plt
-            import time
 
             ta = time.perf_counter()
             x_mask = (stancil[0] <= X) & (X <= stancil[-1])
@@ -277,8 +275,6 @@ def calc_gFT_apply(stancil, prior):
             FT = generalized_Fourier(x, y, self.k)
 
             if plot_flag:
-                import matplotlib.pyplot as plt
-
                 plt.figure(figsize=(3.34, 1.8), dpi=300)
 
             # define weights. Weights are normalized to 1
@@ -338,10 +334,7 @@ def calc_gFT_apply(stancil, prior):
             inverse_stats = FT.get_stats(self.dk, Lpoints_full, print_flag=plot_flag)
             # add fitting parameters of Prior to stats dict
             for k, I in prior_pars.items():
-                try:
-                    inverse_stats[k] = I.value
-                except:
-                    inverse_stats[k] = np.nan
+                inverse_stats[k] = I.value if hasattr(I, "value") else np.nan
 
             print("compute time stats : ", time.perf_counter() - ta)
 
@@ -634,15 +627,11 @@ def calc_var(self):
 
     def parceval(self, add_attrs=True, weight_data=False):
         "test Parceval theorem"
-        import copy
 
         DATA = self.data
-        L = self.Lmeters
         X = self.x
 
         def get_stancil_var_apply(stancil):
-            from scipy.signal import detrend
-
             "returns the variance of yy for stancil"
             x_mask = (stancil[0] < X) & (X <= stancil[-1])
             idata = DATA[x_mask]
@@ -756,8 +745,6 @@ def __init__(self, x, ydata, k):
         """
         non_dimensionalize (bool, default=True) if True, then the data and R_data_uncertainty is non-dimensionalized by the std of the data
         """
-        import numpy as np
-        from numpy import linalg
 
         self.x, self.ydata, self.k = x, ydata, k
         self.M = self.k.size  # number of wavenumbers
@@ -769,7 +756,7 @@ def __init__(self, x, ydata, k):
             # test if the data is real, not nan and not inf
             assert np.isrealobj(self.ydata), "data is not real"
             assert np.isfinite(self.ydata).all(), "data is not finite"
-            assert np.isnan(self.ydata).all() == False, "data is not nan"
+            assert not np.isnan(self.ydata).all(), "data is not nan"
 
     # data matrix
     def get_H(self, xx=None):
@@ -792,8 +779,6 @@ def define_problem(self, P_weight, R_data_uncertainty):
         self.R_1d = R_data_uncertainty
 
     def solve(self):
-        from numpy import linalg
-
         inv = linalg.inv
         """ 
         solves the linear inverse problem, return hessian and p_hat
@@ -865,7 +850,6 @@ def parceval(self, dk, Nx_full):
     def get_stats(self, dk, Nx_full, print_flag=False):
         residual = self.ydata - self.model()
 
-        Lmeters = self.x[-1] - self.x[0]
         pars = {
             "data_var": self.ydata.var(),
             "model_var": self.model().var(),
@@ -896,8 +880,6 @@ def get_stats(self, dk, Nx_full, print_flag=False):
 
 class get_prior_spec:
     def __init__(self, freq, data):
-        import lmfit as LM
-
         self.LM = LM
         self.data = data
         self.freq = freq
@@ -919,7 +901,6 @@ def set_parameters(self, flim=None):
         self.params LMfit.parameters class needed for optimization
 
         """
-        import numpy as np
 
         params = self.LM.Parameters()
 
@@ -949,9 +930,6 @@ def model_func(self, f, params):
         )
 
     def non_dim_spec_model(self, f, f_max, amp, gamma=1, angle_rad=0):
-        import icesat2_tracks.local_modules.JONSWAP_gamma as spectal_models
-
-        U = 20  # results are incensitive to U
         f_true = f * np.cos(angle_rad)
         model = spectal_models.JONSWAP_default_alt(f_true, f_max, 20, gamma=gamma)
         model = amp * model / np.nanmean(model)
@@ -984,13 +962,9 @@ def optimize(self, fitting_args=None, method="dual_annealing", max_nfev=None):
         return self.fitter
 
     def plot_data(self):
-        import matplotlib.pyplot as plt
-
         plt.plot(self.freq, self.data, "k")
 
     def plot_model(self, pars):
-        import matplotlib.pyplot as plt
-
         plt.plot(self.freq, self.model_func(self.freq, pars), "b--")
 
     def runningmean(self, var, m, tailcopy=False):