allow colored_noise to output normalized series and use it to generat…

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 = []

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
 

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,23 +149,23 @@ 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]
 
         ts_M = pyleo.MultipleSeries(serieslist)
 
         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)

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()

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
 

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)
+    tsutils.interp(unevenly_spaced_series.time,unevenly_spaced_series.value,start=start,stop=stop,step=step,step_style=step_style)

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,13 +380,20 @@ 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
 
     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