# Path: pyleoclim/utils/tsmodel.py
# def ar1_sim(y, p, t=None):
#     ''' Produce p realizations of an AR(1) process of length n with lag-1 autocorrelation g calculated from `y` and (if provided) `t`
# 
#     Parameters
#     ----------
# 
#     y : array
#         a time series; NaNs not allowed
#     p : int
#         column dimension (number of surrogates)
#     t : array
#         the time axis of the series
# 
#     Returns
#     -------
# 
#     ysim : array
#         n by p matrix of simulated AR(1) vector
# 
#     See Also
#     --------
# 
#     pyleoclim.utils.tsmodel.ar1_model : Simulates a (possibly irregularly-sampled) AR(1) process with given decay constant tau, à la REDFIT.
# 
#     pyleoclim.utils.tsmodel.ar1_fit : Returns the lag-1 autocorrelation from AR(1) fit OR persistence from tauest.
# 
#     pyleoclim.utils.tsmodel.ar1_fit_evenly : Returns the lag-1 autocorrelation from AR(1) fit assuming even temporal spacing.
# 
#     pyleoclim.utils.tsmodel.tau_estimation : Estimates the  temporal decay scale of an (un)evenly spaced time series.
# 
#     '''
#     n = np.size(y)
#     ysim = np.empty(shape=(n, p))  # declare array
# 
#     sig = np.std(y)
#     if is_evenly_spaced(t):
#         g = ar1_fit_evenly(y)
# 
#         # specify model parameters (statmodel want lag0 coefficent as unity)
#         ar = np.r_[1, -g]  # AR model parameter
#         ma = np.r_[1, 0.0]  # MA model parameters
#         sig_n = sig*np.sqrt(1-g**2)  # theoretical noise variance for the process to achieve the same variance as y
# 
#         # simulate AR(1) model for each column
#         for i in np.arange(p):
#             #ysim[:, i] = sm.tsa.arma_generate_sample(ar=ar, ma=ma, nsample=n, burnin=50, sigma=sig_n) # old statsmodels syntax
#             #ysim[:, i] = sm.tsa.ArmaProcess(ar, ma).generate_sample(nsample=n, scale=sig_n, burnin=50) # statsmodels v0.11.1-?
#             ysim[:, i] = arma_generate_sample(ar, ma, nsample=n, scale=sig_n, burnin=50) # statsmodels v0.12+
#     else:
#         #  tau_est = ar1_fit(y, t=t, detrend=detrend, params=params)
#         tau_est = tau_estimation(y, t)
#         for i in np.arange(p):
#             # the output of ar1_model has unit variance,
#             # multiply by sig to be consistent with the original input timeseries
#             ysim[:, i] = ar1_model(t, tau_est, output_sigma=sig)
# 
#     if p == 1:
#         ysim = ysim[:, 0]
# 
#     return ysim
# 
# def colored_noise(alpha, t, f0=None, m=None, seed=None):
#     ''' Generate a colored noise timeseries
# 
#     Parameters
#     ----------
#     alpha : float
#         exponent of the 1/f^alpha noise
# 
#     t : float
#         time vector of the generated noise
# 
#     f0 : float
#         fundamental frequency
# 
#     m : int
#         maximum number of the waves, which determines the highest frequency of the components in the synthetic noise
# 
#     Returns
#     -------
# 
#     y : array
#         the generated 1/f^alpha noise
# 
#     References
#     ----------
# 
#     Eq. (15) in Kirchner, J. W. Aliasing in 1/f(alpha) noise spectra: origins, consequences, and remedies.
#         Phys Rev E Stat Nonlin Soft Matter Phys 71, 066110 (2005).
#     '''
#     n = np.size(t)  # number of time points
#     y = np.zeros(n)
# 
#     if f0 is None:
#         f0 = 1/n  # fundamental frequency
#     if m is None:
#         m = n//2
# 
#     k = np.arange(m) + 1  # wave numbers
# 
#     if seed is not None:
#         np.random.seed(seed)
# 
#     theta = np.random.rand(int(m))*2*np.pi  # random phase
#     for j in range(n):
#         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


import numpy as np
import pandas as pd
import pytest
import pyleoclim as pyleo
from numpy.testing import assert_array_equal
from pandas.testing import assert_frame_equal
from pyleoclim.utils.tsmodel import (
    ar1_sim,
    colored_noise,
)

Naming rules:
1. class: Test{filename}{Class}{method} with appropriate camel case
2. function: test_{method}_t{test_id}

Notes on how to test:
0. Make sure [pytest](https://docs.pytest.org) has been installed: `pip install pytest`
1. execute `pytest {directory_path}` in terminal to perform all tests in all testing files inside the specified directory
2. execute `pytest {file_path}` in terminal to perform all tests in the specified file
3. execute `pytest {file_path}::{TestClass}::{test_method}` in terminal to perform a specific test class/method inside the specified file
4. after `pip install pytest-xdist`, one may execute "pytest -n 4" to test in parallel with number of workers specified by `-n`
5. for more details, see https://docs.pytest.org/en/stable/usage.html
'''




# a collection of useful functions

def gen_normal(loc=0, scale=1, nt=100):
    ''' Generate random data with a Gaussian distribution
    '''
    t = np.arange(nt)
    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):
    ''' Generate colored noise
    '''
    t = np.arange(nt)

# Path: pyleoclim/utils/tsmodel.py
# def ar1_sim(y, p, t=None):
#     ''' Produce p realizations of an AR(1) process of length n with lag-1 autocorrelation g calculated from `y` and (if provided) `t`
# 
#     Parameters
#     ----------
# 
#     y : array
#         a time series; NaNs not allowed
#     p : int
#         column dimension (number of surrogates)
#     t : array
#         the time axis of the series
# 
#     Returns
#     -------
# 
#     ysim : array
#         n by p matrix of simulated AR(1) vector
# 
#     See Also
#     --------
# 
#     pyleoclim.utils.tsmodel.ar1_model : Simulates a (possibly irregularly-sampled) AR(1) process with given decay constant tau, à la REDFIT.
# 
#     pyleoclim.utils.tsmodel.ar1_fit : Returns the lag-1 autocorrelation from AR(1) fit OR persistence from tauest.
# 
#     pyleoclim.utils.tsmodel.ar1_fit_evenly : Returns the lag-1 autocorrelation from AR(1) fit assuming even temporal spacing.
# 
#     pyleoclim.utils.tsmodel.tau_estimation : Estimates the  temporal decay scale of an (un)evenly spaced time series.
# 
#     '''
#     n = np.size(y)
#     ysim = np.empty(shape=(n, p))  # declare array
# 
#     sig = np.std(y)
#     if is_evenly_spaced(t):
#         g = ar1_fit_evenly(y)
# 
#         # specify model parameters (statmodel want lag0 coefficent as unity)
#         ar = np.r_[1, -g]  # AR model parameter
#         ma = np.r_[1, 0.0]  # MA model parameters
#         sig_n = sig*np.sqrt(1-g**2)  # theoretical noise variance for the process to achieve the same variance as y
# 
#         # simulate AR(1) model for each column
#         for i in np.arange(p):
#             #ysim[:, i] = sm.tsa.arma_generate_sample(ar=ar, ma=ma, nsample=n, burnin=50, sigma=sig_n) # old statsmodels syntax
#             #ysim[:, i] = sm.tsa.ArmaProcess(ar, ma).generate_sample(nsample=n, scale=sig_n, burnin=50) # statsmodels v0.11.1-?
#             ysim[:, i] = arma_generate_sample(ar, ma, nsample=n, scale=sig_n, burnin=50) # statsmodels v0.12+
#     else:
#         #  tau_est = ar1_fit(y, t=t, detrend=detrend, params=params)
#         tau_est = tau_estimation(y, t)
#         for i in np.arange(p):
#             # the output of ar1_model has unit variance,
#             # multiply by sig to be consistent with the original input timeseries
#             ysim[:, i] = ar1_model(t, tau_est, output_sigma=sig)
# 
#     if p == 1:
#         ysim = ysim[:, 0]
# 
#     return ysim
# 
# def colored_noise(alpha, t, f0=None, m=None, seed=None):
#     ''' Generate a colored noise timeseries
# 
#     Parameters
#     ----------
#     alpha : float
#         exponent of the 1/f^alpha noise
# 
#     t : float
#         time vector of the generated noise
# 
#     f0 : float
#         fundamental frequency
# 
#     m : int
#         maximum number of the waves, which determines the highest frequency of the components in the synthetic noise
# 
#     Returns
#     -------
# 
#     y : array
#         the generated 1/f^alpha noise
# 
#     References
#     ----------
# 
#     Eq. (15) in Kirchner, J. W. Aliasing in 1/f(alpha) noise spectra: origins, consequences, and remedies.
#         Phys Rev E Stat Nonlin Soft Matter Phys 71, 066110 (2005).
#     '''
#     n = np.size(t)  # number of time points
#     y = np.zeros(n)
# 
#     if f0 is None:
#         f0 = 1/n  # fundamental frequency
#     if m is None:
#         m = n//2
# 
#     k = np.arange(m) + 1  # wave numbers
# 
#     if seed is not None:
#         np.random.seed(seed)
# 
#     theta = np.random.rand(int(m))*2*np.pi  # random phase
#     for j in range(n):
#         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


# Path: pyleoclim/tests/test_ui_PSD.py
import numpy as np
import pandas as pd
import pytest
import pyleoclim as pyleo
from numpy.testing import assert_array_equal
from pandas.testing import assert_frame_equal
from pyleoclim.utils.tsmodel import (
    ar1_sim,
    colored_noise,
)


Naming rules:
1. class: Test{filename}{Class}{method} with appropriate camel case
2. function: test_{method}_t{test_id}

Notes on how to test:
0. Make sure [pytest](https://docs.pytest.org) has been installed: `pip install pytest`
1. execute `pytest {directory_path}` in terminal to perform all tests in all testing files inside the specified directory
2. execute `pytest {file_path}` in terminal to perform all tests in the specified file
3. execute `pytest {file_path}::{TestClass}::{test_method}` in terminal to perform a specific test class/method inside the specified file
4. after `pip install pytest-xdist`, one may execute "pytest -n 4" to test in parallel with number of workers specified by `-n`
5. for more details, see https://docs.pytest.org/en/stable/usage.html
'''




# a collection of useful functions

def gen_normal(loc=0, scale=1, nt=100):
    ''' Generate random data with a Gaussian distribution
    '''
    t = np.arange(nt)
    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):
    ''' Generate colored noise
    '''
    t = np.arange(nt)

    v = colored_noise(alpha=alpha, t=t, f0=f0, m=m, seed=seed)