Surrogate bug fixes (#545)

* bug fixes

including #544

* fix plot labeling

pyleoclim/core/ensembleseries.py

@@ -36,8 +36,9 @@ class EnsembleSeries(MultipleSeries):
     and visualization (e.g., envelope plot) that are unavailable to other classes.
 
     '''
-    def __init__(self, series_list):
+    def __init__(self, series_list, label=None):
         self.series_list = series_list
+        self.label = label
 
     @classmethod
     def from_AgeEnsembleArray(self, series, age_array, value_depth = None, age_depth = None, extrapolate=True,verbose=True):
@@ -681,7 +682,10 @@ def plot_traces(self, figsize=[10, 4], xlabel=None, ylabel=None, title=None, num
 
         if title is not None:
             ax.set_title(title)
-
+        else:
+            if self.label is not None:
+                ax.set_title(self.label)
+
         if plot_legend:
             lgd_args = {'frameon': False}
             lgd_args.update(lgd_kwargs)

pyleoclim/core/multipleseries.py

@@ -1651,7 +1651,7 @@ def plot(self, figsize=[10, 4],
             fig, ax = plt.subplots(figsize=figsize)
 
         if title is None and self.label is not None:
-            ax.set_title(self.label, fontweight='bold')
+            ax.set_title(self.label)
 
         if ylabel is None:
             consistent_ylabels = True

pyleoclim/core/surrogateseries.py

@@ -76,13 +76,13 @@ def __init__(self, series_list=[], number=1, method=None, label=None, param=None
         # refine the display name
         if label is None:
             if method == 'ar1sim':
-                self.label =  "AR(1) (MoM)"
+                self.label =  "AR(1) surrogates (MoM)"
             elif method == 'phaseran':
-                self.label = "phase-randomized"
+                self.label = "Phase-randomized surrogates"
             elif method == 'uar1':
-                self.label = "AR(1) (MLE)"
+                self.label = "AR(1) surrogates (MLE)"
             elif method == 'CN':
-                self.label = r'$f^{-\beta}$'
+                self.label = r'Power-law surrogates ($S(f) \propto f^{-\beta}$)'
             else:
                 raise ValueError(f"Unknown method: {self.method}. Please use one of {supported_surrogates}")
 
@@ -131,7 +131,9 @@ def from_series(self, target_series):
 
         # apply surrogate method
         if self.method == 'ar1sim':
-                y_surr = tsmodel.ar1_sim(target_series.value, self.number, target_series.time)  
+            tsi = target_series if target_series.is_evenly_spaced() else target_series.interp()
+            mu = tsi.value.mean()
+            y_surr = tsmodel.ar1_sim(target_series.value, self.number, target_series.time) + mu  
 
         elif self.method == 'phaseran':
             if target_series.is_evenly_spaced():
@@ -142,32 +144,35 @@ def from_series(self, target_series):
         elif self.method == 'uar1':
             # estimate theta with MLE
             tau, sigma_2 = tsmodel.uar1_fit(target_series.value, target_series.time)
+            tsi = target_series if target_series.is_evenly_spaced() else target_series.interp()
+            mu = tsi.value.mean()
             self.param = [tau, sigma_2] # assign parameters for future use
             # generate time axes according to provided pattern
             times = np.squeeze(np.tile(target_series.time, (self.number, 1)).T)  
-            # generate matrix
-            y_surr = tsmodel.uar1_sim(t = times, tau=tau, sigma_2=sigma_2)
+            # generate matrix and add  the mean
+            y_surr = tsmodel.uar1_sim(t = times, tau=tau, sigma_2=sigma_2) + mu
 
         elif self.method == 'CN':
             tsi = target_series if target_series.is_evenly_spaced() else target_series.interp()
+            mu = tsi.value.mean()
             sigma = tsi.value.std()
             alpha = tsi.spectral(method='cwt').beta_est().beta_est_res['beta'] # fit the parameter using the smoothest spectral method
             self.param = [alpha]
             y_surr = np.empty((len(target_series.time),self.number))
             for i in range(self.number):
-                y_surr[:,i] = tsmodel.colored_noise(alpha=alpha,t=target_series.time, std = sigma)
+                y_surr[:,i] = tsmodel.colored_noise(alpha=alpha,t=target_series.time, std = sigma) + mu
 
         if self.number > 1:
             s_list = []
             for i, y in enumerate(y_surr.T):
                 ts = target_series.copy() # copy Series
-                ts.value = y # replace value with y_surr column
-                ts.label = str(target_series.label or '') + " surr #" + str(i+1)
+                ts.value = y 
+                ts.label = str(target_series.label or '') + " #" + str(i+1)
                 s_list.append(ts)
         else:
             ts = target_series.copy() # copy Series
-            ts.value = y_surr # replace value with y_surr column
-            ts.label = str(target_series.label or '') + " surr" 
+            ts.value = y_surr
+            ts.label = str(target_series.label or '') 
             s_list = [ts]
 
         self.series_list = s_list
@@ -238,7 +243,10 @@ def from_param(self, param, length=50, time_pattern = 'even', settings=None):
             if "time" not in settings:
                 raise ValueError("'time' not found in settings")
             else:
-                times =  np.tile(settings["time"], (self.number, 1)).T  
+                if self.number > 1:
+                    times =  np.tile(settings["time"], (self.number, 1)).T 
+                else:
+                    times = settings["time"]
         else:
             raise ValueError(f"Unknown time pattern: {time_pattern}")
 
@@ -263,9 +271,15 @@ def from_param(self, param, length=50, time_pattern = 'even', settings=None):
 
         # create the series_list    
         s_list = []
-        for i, (t, y) in enumerate(zip(times.T,y_surr.T)):
-            ts = Series(time=t, value=y,  
-                           label = str(self.label or '') + " surr #" + str(i+1),
+        if self.number > 1:
+            for i, (t, y) in enumerate(zip(times.T,y_surr.T)):
+                ts = Series(time=t, value=y,  
+                               label = str(self.label or '') + " #" + str(i+1),
+                               verbose=False, auto_time_params=True)
+                s_list.append(ts)
+        else:
+            ts = Series(time=times, value=y_surr,  
+                           label = str(self.label or '') + " #`",
                            verbose=False, auto_time_params=True)
             s_list.append(ts)
 

pyleoclim/tests/test_core_EnsembleSeries.py

@@ -239,8 +239,9 @@ def test_plot_envelope_t0(self):
 
         fig, ax = ts_ens.plot_envelope(curve_lw=1.5)
         pyleo.closefig(fig)
-
-    def test_plot_traces_t0(self):
+
+    @pytest.mark.parametrize('label', ['ensemble', None])
+    def test_plot_traces_t0(self,label):
         ''' Test EnsembleSeries.plot_traces() on a list of colored noise
         '''
         nn = 30 # number of noise realizations
@@ -254,9 +255,9 @@ def test_plot_traces_t0(self):
             ts = pyleo.Series(time=signal.time, value=signal.value+noise[:,idx])
             series_list.append(ts)
 
-        ts_ens = pyleo.EnsembleSeries(series_list)
+        ts_ens = pyleo.EnsembleSeries(series_list, label=label)
 
-        fig, ax = ts_ens.plot_traces(alpha=0.2, num_traces=8)
+        fig, ax = ts_ens.plot_traces(alpha=0.2, num_traces=8) # test transparency and num_traces at the same time
         pyleo.closefig(fig)
 
 class TestUIEnsembleSeriesQuantiles():

pyleoclim/utils/tsmodel.py

@@ -415,7 +415,7 @@ def colored_noise(alpha, t, std = 1.0, f0=None, m=None, seed=None):
     y = np.zeros(n)
 
     if f0 is None:
-        f0 = 1/n  # fundamental frequency
+        f0 = 1/np.ptp(t)  # fundamental frequency
     if m is None:
         m = n//2