From 3b02a5e509134722803c7f89eb093e10110d7e1d Mon Sep 17 00:00:00 2001 From: CommonClimate Date: Fri, 17 May 2024 15:14:56 -0700 Subject: [PATCH 1/2] bug fixes including #544 --- pyleoclim/core/ensembleseries.py | 4 +++- pyleoclim/core/surrogateseries.py | 36 +++++++++++++++++++------------ pyleoclim/utils/tsmodel.py | 2 +- 3 files changed, 26 insertions(+), 16 deletions(-) diff --git a/pyleoclim/core/ensembleseries.py b/pyleoclim/core/ensembleseries.py index 569daa43..c5352d97 100644 --- a/pyleoclim/core/ensembleseries.py +++ b/pyleoclim/core/ensembleseries.py @@ -681,7 +681,9 @@ def plot_traces(self, figsize=[10, 4], xlabel=None, ylabel=None, title=None, num if title is not None: ax.set_title(title) - + else: + ax.set_title(self.label) + if plot_legend: lgd_args = {'frameon': False} lgd_args.update(lgd_kwargs) diff --git a/pyleoclim/core/surrogateseries.py b/pyleoclim/core/surrogateseries.py index 57a75c33..0b951e92 100644 --- a/pyleoclim/core/surrogateseries.py +++ b/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}") @@ -265,7 +273,7 @@ def from_param(self, param, length=50, time_pattern = 'even', settings=None): 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), + label = str(self.label or '') + " #" + str(i+1), verbose=False, auto_time_params=True) s_list.append(ts) diff --git a/pyleoclim/utils/tsmodel.py b/pyleoclim/utils/tsmodel.py index c91b80b4..62cdfe6a 100644 --- a/pyleoclim/utils/tsmodel.py +++ b/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 From f812b3c0fc6cb217907a854a5647d841c6d0ac39 Mon Sep 17 00:00:00 2001 From: CommonClimate Date: Fri, 17 May 2024 16:16:47 -0700 Subject: [PATCH 2/2] fix plot labeling --- pyleoclim/core/ensembleseries.py | 6 ++++-- pyleoclim/core/multipleseries.py | 2 +- pyleoclim/core/surrogateseries.py | 14 ++++++++++---- pyleoclim/tests/test_core_EnsembleSeries.py | 9 +++++---- 4 files changed, 20 insertions(+), 11 deletions(-) diff --git a/pyleoclim/core/ensembleseries.py b/pyleoclim/core/ensembleseries.py index c5352d97..bbcffb46 100644 --- a/pyleoclim/core/ensembleseries.py +++ b/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): @@ -682,7 +683,8 @@ def plot_traces(self, figsize=[10, 4], xlabel=None, ylabel=None, title=None, num if title is not None: ax.set_title(title) else: - ax.set_title(self.label) + if self.label is not None: + ax.set_title(self.label) if plot_legend: lgd_args = {'frameon': False} diff --git a/pyleoclim/core/multipleseries.py b/pyleoclim/core/multipleseries.py index e3c1aabb..26efb28f 100644 --- a/pyleoclim/core/multipleseries.py +++ b/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 diff --git a/pyleoclim/core/surrogateseries.py b/pyleoclim/core/surrogateseries.py index 0b951e92..5b69a3ab 100644 --- a/pyleoclim/core/surrogateseries.py +++ b/pyleoclim/core/surrogateseries.py @@ -80,7 +80,7 @@ def __init__(self, series_list=[], number=1, method=None, label=None, param=None elif method == 'phaseran': self.label = "Phase-randomized surrogates" elif method == 'uar1': - self.label = "AR(1) surrogates(MLE)" + self.label = "AR(1) surrogates (MLE)" elif method == 'CN': self.label = r'Power-law surrogates ($S(f) \propto f^{-\beta}$)' else: @@ -271,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 '') + " #" + 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) diff --git a/pyleoclim/tests/test_core_EnsembleSeries.py b/pyleoclim/tests/test_core_EnsembleSeries.py index 33b39bd9..6fd72a16 100644 --- a/pyleoclim/tests/test_core_EnsembleSeries.py +++ b/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():