diff --git a/pyleoclim/core/surrogateseries.py b/pyleoclim/core/surrogateseries.py index 060407fd..57a75c33 100644 --- a/pyleoclim/core/surrogateseries.py +++ b/pyleoclim/core/surrogateseries.py @@ -150,11 +150,12 @@ def from_series(self, target_series): elif self.method == 'CN': tsi = target_series if target_series.is_evenly_spaced() else target_series.interp() + 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 + 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) + y_surr[:,i] = tsmodel.colored_noise(alpha=alpha,t=target_series.time, std = sigma) if self.number > 1: s_list = [] diff --git a/pyleoclim/tests/conftest.py b/pyleoclim/tests/conftest.py index e00242e9..2fe420c0 100644 --- a/pyleoclim/tests/conftest.py +++ b/pyleoclim/tests/conftest.py @@ -57,7 +57,7 @@ def unevenly_spaced_series_nans(): t = np.linspace(1,length,length) ** 2 v = np.ones(length) v[2:4] = np.nan - series = pyleo.Series(time=t, value =v, dropna=False, verbose=False) + series = pyleo.Series(time=t, value =v, dropna=False, verbose=False, auto_time_params=True) return series @pytest.fixture @@ -66,15 +66,16 @@ def evenly_spaced_series(): length = 10 t = np.linspace(1,length,length) v = np.cos(2*np.pi*t/10) - series = pyleo.Series(time=t, value=v, verbose=False) + series = pyleo.Series(time=t, value=v, verbose=False, auto_time_params=True) series.label = 'cosine' return series @pytest.fixture def pinkseries(): """Pyleoclim geoseries with 1/f spectrum """ - t,v = pyleo.utils.gen_ts(model='colored_noise',alpha=1.0, nt=100, seed=251) - ts = pyleo.Series(t,v, sort_ts='none', verbose=False) + t = np.arange(100) + v = pyleo.utils.colored_noise(alpha=1.0, t=t, seed=251, std=2.5) + ts = pyleo.Series(t,v, verbose=False, auto_time_params=True) ts.label = 'pink noise' return ts diff --git a/pyleoclim/tests/test_core_MultipleSeries.py b/pyleoclim/tests/test_core_MultipleSeries.py index 6725ff7c..f8d1b612 100644 --- a/pyleoclim/tests/test_core_MultipleSeries.py +++ b/pyleoclim/tests/test_core_MultipleSeries.py @@ -21,10 +21,7 @@ import pytest import pyleoclim as pyleo -from pyleoclim.utils.tsmodel import ( - ar1_sim, - colored_noise, -) +from pyleoclim.utils.tsmodel import colored_noise # a collection of useful functions @@ -41,11 +38,11 @@ def gen_normal(loc=0, scale=1, nt=100): v = np.random.normal(loc=loc, scale=scale, size=nt) return t, v -def gen_colored_noise(alpha=1, nt=100, f0=None, m=None, seed=None): +def gen_colored_noise(alpha=1, nt=100, std=1.0, f0=None, m=None, seed=None): ''' Generate colored noise ''' t = np.arange(nt) - v = colored_noise(alpha=alpha, t=t, f0=f0, m=m, seed=seed) + v = colored_noise(alpha=alpha, t=t, std=std, f0=f0, m=m, seed=seed) return t, v def load_data(): @@ -98,8 +95,8 @@ def test_plot(self): t_1, v_1 = gen_normal() #Create series objects - ts_0 = pyleo.Series(time = t_0, value = v_0) - ts_1 = pyleo.Series(time = t_1, value = v_1) + ts_0 = pyleo.Series(time = t_0, value = v_0, auto_time_params=True, verbose=False) + ts_1 = pyleo.Series(time = t_1, value = v_1, auto_time_params=True, verbose=False) #Create a list of series objects serieslist = [ts_0, ts_1] @@ -136,8 +133,8 @@ def test_stripes(self): t_1, v_1 = gen_normal() #Create series objects - ts_0 = pyleo.Series(time = t_0, value = v_0, label='series 1') - ts_1 = pyleo.Series(time = t_1, value = v_1, label='series 2') + ts_0 = pyleo.Series(time = t_0, value = v_0, label='series 1', auto_time_params=True, verbose=False) + ts_1 = pyleo.Series(time = t_1, value = v_1, label='series 2', auto_time_params=True, verbose=False) #Turn this list into a multiple series object ts_M = ts_0 & ts_1 @@ -152,11 +149,11 @@ class TestMultipleSeriesStandardize: only now we are running the test on series in a MultipleSeries object''' def test_standardize(self): - t_0, v_0 = gen_colored_noise() - t_1, v_1 = gen_colored_noise() + t_0, v_0 = gen_colored_noise(std = 10) + t_1, v_1 = gen_colored_noise(std= 20) - ts_0 = pyleo.Series(time = t_0, value = v_0) - ts_1 = pyleo.Series(time = t_1, value = v_1) + ts_0 = pyleo.Series(time = t_0, value = v_0, auto_time_params=True, verbose=False) + ts_1 = pyleo.Series(time = t_1, value = v_1, auto_time_params=True, verbose=False) serieslist = [ts_0, ts_1] @@ -164,11 +161,11 @@ def test_standardize(self): ts_M_std = ts_M.standardize() - x_axis_0 = ts_M_std.series_list[0].__dict__['time'] - x_axis_1 = ts_M_std.series_list[1].__dict__['time'] + x_axis_0 = ts_M_std.series_list[0].time + x_axis_1 = ts_M_std.series_list[1].time - y_axis_0 = ts_M_std.series_list[0].__dict__['value'] - y_axis_1 = ts_M_std.series_list[1].__dict__['value'] + y_axis_0 = ts_M_std.series_list[0].value + y_axis_1 = ts_M_std.series_list[1].value assert_array_equal(x_axis_0, t_0) assert_array_equal(x_axis_1, t_1) diff --git a/pyleoclim/tests/test_core_Series.py b/pyleoclim/tests/test_core_Series.py index d96ff6d8..bc664fea 100644 --- a/pyleoclim/tests/test_core_Series.py +++ b/pyleoclim/tests/test_core_Series.py @@ -378,18 +378,10 @@ class TestUISeriesStandardize: Standardize normalizes the series object, so we'll simply test maximum and minimum values''' - def test_standardize(self, pinkseries): + def test_standardize(self, pinkseries, eps=0.1): ts = pinkseries - - #Call function to be tested - ts_std = ts.standardize(keep_log=True) - - #Compare maximum and minimum values - value = ts.__dict__['value'] - value_std = ts_std.__dict__['value'] - - assert max(value) > max(value_std) - assert min(value) < min(value_std) + ts_std = ts.standardize(keep_log=False) + assert np.abs(ts_std.value.std()-1) < eps # current class TestUISeriesClean: '''Test for Series.clean() diff --git a/pyleoclim/tests/test_utils_tsmodel.py b/pyleoclim/tests/test_utils_tsmodel.py index a7e250f2..f2144f24 100644 --- a/pyleoclim/tests/test_utils_tsmodel.py +++ b/pyleoclim/tests/test_utils_tsmodel.py @@ -67,5 +67,11 @@ def test_uar1_fit(p, evenly_spaced, tol = 0.45): # test that assert np.abs(theta_hat_bar[0]-tau) < tol assert np.abs(theta_hat_bar[1]-sigma_2) < tol - + +@pytest.mark.parametrize('std',[None,2]) +def test_colored_noise(std, nt=100, eps = 0.1): + t = np.arange(nt) + v = tsmodel.colored_noise(alpha=1.0, t=t, std = std) + if std is not None: + assert np.abs(v.std() - std) < eps diff --git a/pyleoclim/tests/test_utils_tsutils.py b/pyleoclim/tests/test_utils_tsutils.py index 57a1d0cb..650e0d52 100644 --- a/pyleoclim/tests/test_utils_tsutils.py +++ b/pyleoclim/tests/test_utils_tsutils.py @@ -71,4 +71,4 @@ def test_interp_t1(unevenly_spaced_series): @pytest.mark.parametrize('step',[None,20]) @pytest.mark.parametrize('step_style',[None,'median']) def test_interp_t2(unevenly_spaced_series,start,stop,step,step_style): - tsutils.interp(unevenly_spaced_series.time,unevenly_spaced_series.value,start=start,stop=stop,step=step,step_style=step_style) \ No newline at end of file + tsutils.interp(unevenly_spaced_series.time,unevenly_spaced_series.value,start=start,stop=stop,step=step,step_style=step_style) diff --git a/pyleoclim/utils/tsmodel.py b/pyleoclim/utils/tsmodel.py index 9be5860e..c91b80b4 100644 --- a/pyleoclim/utils/tsmodel.py +++ b/pyleoclim/utils/tsmodel.py @@ -6,7 +6,7 @@ from statsmodels.tsa.arima_process import arma_generate_sample from statsmodels.tsa.arima.model import ARIMA # from tqdm import tqdm -# from .tsutils import standardize +from .tsutils import standardize # from stochastic.processes.noise import ColoredNoise # from stochastic.processes.noise import FractionalGaussianNoise @@ -370,7 +370,7 @@ def sm_ar1_sim(n, p, g, sig): return red -def colored_noise(alpha, t, f0=None, m=None, seed=None): +def colored_noise(alpha, t, std = 1.0, f0=None, m=None, seed=None): ''' Generate a colored noise timeseries Parameters @@ -380,6 +380,9 @@ def colored_noise(alpha, t, f0=None, m=None, seed=None): t : float time vector of the generated noise + + std : float + standard deviation of the series. defaults to 1.0 f0 : float fundamental frequency @@ -387,6 +390,10 @@ def colored_noise(alpha, t, f0=None, m=None, seed=None): m : int maximum number of the waves, which determines the highest frequency of the components in the synthetic noise + seed : int + seed for the random number generator + + Returns ------- @@ -422,8 +429,13 @@ def colored_noise(alpha, t, f0=None, m=None, seed=None): coeff = (k*f0)**(-alpha/2) sin_func = np.sin(2*np.pi*k*f0*t[j] + theta) y[j] = np.sum(coeff*sin_func) - - return y + + if std is not None: + ys, _, _ = standardize(y,scale=std) # rescale + else: + ys = y + + return ys def colored_noise_2regimes(alpha1, alpha2, f_break, t, f0=None, m=None, seed=None): ''' Generate a colored noise timeseries with two regimes