From 3b02a5e509134722803c7f89eb093e10110d7e1d Mon Sep 17 00:00:00 2001 From: CommonClimate Date: Fri, 17 May 2024 15:14:56 -0700 Subject: [PATCH 1/5] 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/5] 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(): From 143ca6928cd86528a4bbdc7854dbede1aeb3debb Mon Sep 17 00:00:00 2001 From: CommonClimate Date: Sat, 18 May 2024 12:00:25 -0700 Subject: [PATCH 3/5] bug fixes when number =1 --- pyleoclim/core/surrogateseries.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/pyleoclim/core/surrogateseries.py b/pyleoclim/core/surrogateseries.py index 5b69a3ab..f2472fcd 100644 --- a/pyleoclim/core/surrogateseries.py +++ b/pyleoclim/core/surrogateseries.py @@ -263,8 +263,11 @@ def from_param(self, param, length=50, time_pattern = 'even', settings=None): if nparam >1: raise ValueError(f"1 parameter is needed for this model; only the first of the provided {nparam} will be used") y_surr = np.empty((length,self.number)) - for i in range(self.number): - y_surr[:,i] = tsmodel.colored_noise(alpha=param[0],t=times[:,i]) + if self.number > 1: + for i in range(self.number): + y_surr[:,i] = tsmodel.colored_noise(alpha=param[0],t=times[:,i]) + else: + y_surr = tsmodel.colored_noise(alpha=param[0],t=times) elif self.method == 'phaseran': raise ValueError("Phase-randomization is only available in from_series().") @@ -278,7 +281,7 @@ def from_param(self, param, length=50, time_pattern = 'even', settings=None): verbose=False, auto_time_params=True) s_list.append(ts) else: - ts = Series(time=times, value=y_surr, + ts = Series(time=np.squeeze(times), value=np.squeeze(y_surr), label = str(self.label or '') + " #`", verbose=False, auto_time_params=True) s_list.append(ts) From 388773ea3ba1fee3501c66d5599742357adc3c6f Mon Sep 17 00:00:00 2001 From: Alexander James Date: Mon, 20 May 2024 11:50:34 -0700 Subject: [PATCH 4/5] surrogate series tests update --- pyleoclim/tests/test_core_SurrogateSeries.py | 131 +++++++++++-------- 1 file changed, 74 insertions(+), 57 deletions(-) diff --git a/pyleoclim/tests/test_core_SurrogateSeries.py b/pyleoclim/tests/test_core_SurrogateSeries.py index 51af8689..06dfdd63 100644 --- a/pyleoclim/tests/test_core_SurrogateSeries.py +++ b/pyleoclim/tests/test_core_SurrogateSeries.py @@ -19,17 +19,66 @@ from numpy.testing import assert_array_equal from scipy import stats -class TestUISurrogatesSeries: - ''' Tests for pyleoclim.core.ui.SurrogateSeries - ''' +class TestSurrogateSeriesFromSeries: + '''Tests for SurrogateSeries.from_series() method''' + @pytest.mark.parametrize('method', ['uar1','ar1sim','phaseran']) + @pytest.mark.parametrize('number', [1,10]) + @pytest.mark.parametrize('seed', [None,42]) + def test_fromseries_t0(self,gen_ts,method,number,seed): + '''Test from_series method''' + series = gen_ts + surr = pyleo.SurrogateSeries(method=method,number=number,seed=seed) + surr.from_series(target_series=series,method=method,number=number,seed=seed) + + @pytest.mark.parametrize('method',['ar1sim','uar1','phaseran','CN']) + @pytest.mark.parametrize('number',[1,5]) + def test_from_series_match(self, method, number): + ''' Test that from_series() works with all methods, matches the original + time axis AND that number can be varied + ''' + t, v = pyleo.utils.gen_ts(model='colored_noise', nt=100) + ts = pyleo.Series(time = t, value = v, verbose=False, auto_time_params=True) + surr = pyleo.SurrogateSeries(method=method, number=number) + surr.from_series(ts) + assert len(surr.series_list) == number # test right number + for s in surr.series_list: + assert_array_equal(s.time, ts.time) # test time axis match + +class TestSurrogateSeriesFromParam: + '''Tests for SurrogateSeries.from_param() method''' + @pytest.mark.parametrize('method, noise_param', [('uar1',[1,5]), + ('CN',[1]), + ('ar1sim',[1,5]),]) + @pytest.mark.parametrize('number', [1,10]) + @pytest.mark.parametrize('seed', [None,42]) + @pytest.mark.parametrize('time_pattern,settings',[('even',None), + ('random',None), + ('specified',{'time':np.arange(50)})]) + def test_from_param_t0(self,method,noise_param,number,seed,time_pattern,settings): + '''Test from_series method''' + surr = pyleo.SurrogateSeries(method=method,number=number,seed=seed) + surr.from_param(noise_param=noise_param,time_pattern=time_pattern,settings=settings) + + @pytest.mark.parametrize('delta_t',[1,3]) + @pytest.mark.parametrize('length',[10,20]) + def test_from_param_even(self, delta_t, length, number=2): + ''' Test from_param() with even time axes with varying resolution + ''' + surr = pyleo.SurrogateSeries(method='ar1sim', number=number) + surr.from_param(param = [5,1], time_pattern='even', length= length, + settings={"delta_t" :delta_t}) + for ts in surr.series_list: + assert(np.allclose(tsmodel.inverse_cumsum(ts.time),delta_t)) + assert len(ts.time) == length + @pytest.mark.parametrize('method',['ar1sim','uar1']) @pytest.mark.parametrize('param',[[5,1],[2,2]]) # stacking pytests is legal! - def test_surrogates_ar1(self, method, param, nsim=10, eps=1, seed = 108): + def test_from_param_ar1(self, method, param, nsim=10, eps=1, seed = 108): ''' Generate AR(1) surrogates based on a AR(1) series with certain parameters, estimate the parameters from the surrogates series and verify accuracy ''' - ar1 = pyleo.SurrogateSeries(method='ar1sim', number=nsim, seed=seed) + ar1 = pyleo.SurrogateSeries(method=method, number=nsim, seed=seed) ar1.from_param(param = param) tau = np.empty((nsim)) sigma_2 = np.empty_like(tau) @@ -37,9 +86,9 @@ def test_surrogates_ar1(self, method, param, nsim=10, eps=1, seed = 108): tau[i], sigma_2[i] = tsmodel.uar1_fit(ts.value, ts.time) assert np.abs(tau.mean()-param[0]) < eps assert np.abs(sigma_2.mean()-param[1]) < eps - + @pytest.mark.parametrize('param',[1.0,[1.0]]) - def test_surrogates_CN(self, param, nsim=10, eps=.1, seed = 108): + def test_from_param_CN_t0(self, param, nsim=10, eps=.1, seed = 108): ''' Generate colored noise surrogates based known parameters, estimate the parameters from the surrogates series and verify accuracy ''' @@ -50,41 +99,27 @@ def test_surrogates_CN(self, param, nsim=10, eps=.1, seed = 108): beta[i] = ts.interp().spectral(method='cwt').beta_est().beta_est_res['beta'] assert np.abs(beta.mean()-1.0) < eps - def test_surrogates_CN2(self, nsim=10): + def test_from_param_CN_t1(self, nsim=10): ''' Generate colored noise surrogates with too many params; check error msg ''' CN = pyleo.SurrogateSeries(method='CN', number=nsim) with pytest.raises(ValueError): CN.from_param(param = [2.0, 3.0], length=200) # check that this returns an error - - @pytest.mark.parametrize('method',['ar1sim','uar1','phaseran','CN']) - @pytest.mark.parametrize('number',[1,5]) - def test_surrogates_match(self, method, number): - ''' Test that from_series() works with all methods, matches the original - time axis AND that number can be varied - ''' - t, v = pyleo.utils.gen_ts(model='colored_noise', nt=100) - ts = pyleo.Series(time = t, value = v, verbose=False, auto_time_params=True) - surr = pyleo.SurrogateSeries(method=method, number=number) - surr.from_series(ts) - assert len(surr.series_list) == number # test right number - for s in surr.series_list: - assert_array_equal(s.time, ts.time) # test time axis match - - @pytest.mark.parametrize('delta_t',[1,3]) - @pytest.mark.parametrize('length',[10,20]) - def test_surrogates_from_param_even(self, delta_t, length, number=2): - ''' Test from_param() with even time axes with varying resolution + + @pytest.mark.parametrize('length', [10, 100]) + def test_from_param_specified_time(self, length, number=2): + ''' Test AR(1) model with specified time axis ''' - surr = pyleo.SurrogateSeries(method='ar1sim', number=number) - surr.from_param(param = [5,1], time_pattern='even', length= length, - settings={"delta_t" :delta_t}) + ti = tsmodel.random_time_axis(n=length) + surr = pyleo.SurrogateSeries(method='ar1sim', number=number, seed=108) + surr.from_param(param = [5,1], time_pattern='specified', + settings={"time":ti}) + for ts in surr.series_list: - assert(np.allclose(tsmodel.inverse_cumsum(ts.time),delta_t)) - assert len(ts.time) == length - + assert_array_equal(ts.time, ti) # test time axis match + @pytest.mark.parametrize('dist', ["exponential", "poisson"]) - def test_surrogates_random_time_t0(self, dist, number=2, tol = 3, param=[1]): + def test_from_param_random_time_t0(self, dist, number=2, tol = 3, param=[1]): ''' Test from_param() with random time axes with two 1-parameter distributions ''' surr = pyleo.SurrogateSeries(method='ar1sim', number=number, seed=108) @@ -109,7 +144,7 @@ def test_surrogates_random_time_t0(self, dist, number=2, tol = 3, param=[1]): l2_norm = np.linalg.norm(empirical_cdf - theoretical_cdf) assert(l2_norm Date: Mon, 20 May 2024 12:01:41 -0700 Subject: [PATCH 5/5] surrogate_series test fixes --- pyleoclim/tests/test_core_SurrogateSeries.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyleoclim/tests/test_core_SurrogateSeries.py b/pyleoclim/tests/test_core_SurrogateSeries.py index 06dfdd63..61444802 100644 --- a/pyleoclim/tests/test_core_SurrogateSeries.py +++ b/pyleoclim/tests/test_core_SurrogateSeries.py @@ -28,7 +28,7 @@ def test_fromseries_t0(self,gen_ts,method,number,seed): '''Test from_series method''' series = gen_ts surr = pyleo.SurrogateSeries(method=method,number=number,seed=seed) - surr.from_series(target_series=series,method=method,number=number,seed=seed) + surr.from_series(target_series=series) @pytest.mark.parametrize('method',['ar1sim','uar1','phaseran','CN']) @pytest.mark.parametrize('number',[1,5]) @@ -57,7 +57,7 @@ class TestSurrogateSeriesFromParam: def test_from_param_t0(self,method,noise_param,number,seed,time_pattern,settings): '''Test from_series method''' surr = pyleo.SurrogateSeries(method=method,number=number,seed=seed) - surr.from_param(noise_param=noise_param,time_pattern=time_pattern,settings=settings) + surr.from_param(param=noise_param,time_pattern=time_pattern,settings=settings) @pytest.mark.parametrize('delta_t',[1,3]) @pytest.mark.parametrize('length',[10,20])