# Path: pyleoclim/utils/tsmodel.py
# def ar1_fit_evenly(y):
#     ''' Returns the lag-1 autocorrelation from AR(1) fit.
# 
#     Parameters
#     ----------
#     y : array
#         vector of (float) numbers as a time series
# 
#     Returns
#     -------
#     g : float
#         lag-1 autocorrelation coefficient
# 
#     '''
#     # syntax compatible with statsmodels v0.11.1
#     #ar1_mod = sm.tsa.ARMA(y, (1, 0), missing='drop').fit(trend='nc', disp=0)
#     # syntax compatible with statsmodels v0.12
#     ar1_mod = ARIMA(y, order = (1, 0, 0), missing='drop',trend='ct').fit()
#     g = ar1_mod.params[2]
# 
#     if g > 1:
#         print('Warning: AR(1) fitted autocorrelation greater than 1; setting to 1-eps^{1/4}')
#         eps = np.spacing(1.0)
#         g = 1.0 - eps**(1/4)
# 
#     return g
# 
# Path: pyleoclim/utils/correlation.py
# def sm_ar1_sim(n, p, g, sig):
#     ''' Produce p realizations of an AR1 process of length n with lag-1 autocorrelation g using statsmodels
# 
#     Parameters
#     ----------
# 
#     n : int
#         row dimensions
#     p : int
#         column dimensions
# 
#     g : float
#         lag-1 autocorrelation coefficient
#         
#     sig : float
#         the standard deviation of the original time series
# 
#     Returns
#     -------
# 
#     red : numpy matrix
#         n rows by p columns matrix of an AR1 process
#     '''
#     # specify model parameters (statsmodel wants lag0 coefficents 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 red to achieve the same variance as X
# 
#     red = np.empty(shape=(n, p)) # declare array
# 
#     # simulate AR(1) model for each column
#     for i in np.arange(p):
#         red[:, i] = sm.tsa.arma_generate_sample(ar=ar, ma=ma, nsample=n, burnin=50, scale=sig_n)
# 
#     return red
# 
# def phaseran(recblk, nsurr):
#     ''' Simultaneous phase randomization of a set of time series
#     
#     It creates blocks of surrogate data with the same second order properties as the original
#     time series dataset by transforming the oriinal data into the frequency domain, randomizing the
#     phases simultaneoulsy across the time series and converting the data back into the time domain. 
#     
#     Written by Carlos Gias for MATLAB
# 
#     http://www.mathworks.nl/matlabcentral/fileexchange/32621-phase-randomization/content/phaseran.m
# 
#     Parameters
#     ----------
# 
#     recblk : numpy array
#         2D array , Row: time sample. Column: recording.
#         An odd number of time samples (height) is expected.
#         If that is not the case, recblock is reduced by 1 sample before the surrogate data is created.
#         The class must be double and it must be nonsparse.
#     
#     nsurr : int
#         is the number of image block surrogates that you want to generate.
# 
#     Returns
#     -------
# 
#     surrblk : numpy array
#         3D multidimensional array image block with the surrogate datasets along the third dimension
# 
#     References
#     ----------
# 
#     - Prichard, D., Theiler, J. Generating Surrogate Data for Time Series with Several Simultaneously Measured Variables (1994)
#     Physical Review Letters, Vol 73, Number 7
#     
#     - Carlos Gias (2020). Phase randomization, MATLAB Central File Exchange
#     '''
#     # Get parameters
#     nfrms = recblk.shape[0]
# 
#     if nfrms % 2 == 0:
#         nfrms = nfrms-1
#         recblk = recblk[0:nfrms]
# 
#     len_ser = int((nfrms-1)/2)
#     interv1 = np.arange(1, len_ser+1)
#     interv2 = np.arange(len_ser+1, nfrms)
# 
#     # Fourier transform of the original dataset
#     fft_recblk = np.fft.fft(recblk)
# 
#     surrblk = np.zeros((nfrms, nsurr))
# 
#     #  for k in tqdm(np.arange(nsurr)):
#     for k in np.arange(nsurr):
#         ph_rnd = np.random.rand(len_ser)
# 
#         # Create the random phases for all the time series
#         ph_interv1 = np.exp(2*np.pi*1j*ph_rnd)
#         ph_interv2 = np.conj(np.flipud(ph_interv1))
# 
#         # Randomize all the time series simultaneously
#         fft_recblk_surr = np.copy(fft_recblk)
#         fft_recblk_surr[interv1] = fft_recblk[interv1] * ph_interv1
#         fft_recblk_surr[interv2] = fft_recblk[interv2] * ph_interv2
# 
#         # Inverse transform
#         surrblk[:, k] = np.real(np.fft.ifft(fft_recblk_surr))
# 
#     return surrblk


import numpy as np
from statsmodels.tsa.stattools import grangercausalitytests
from tqdm import tqdm
from .tsmodel import ar1_fit_evenly
from .correlation import sm_ar1_sim, phaseran
from scipy.stats.mstats import mquantiles
                   **kwargs):
    ''' significance test with surrogates with randomized phases

    Parameters
    ----------

    y1, y2 : array
            vectors of (real) numbers with identical length, no NaNs allowed
    method : {'liang'}
            estimates for the Liang method
    npt : int>=1
         time advance in performing Euler forward differencing,
         e.g., 1, 2. Unless the series are generated with a highly chaotic deterministic system,
         npt=1 should be used.
    nsim : int
          the number of surrogates for significance test
    qs : list
        the quantiles for significance test

    Returns
    -------

    res_dict : dict
        A dictionary with the following information:
          T21_noise_qs : list
                        the quantiles of the information flow from noise2 to noise1 for significance testing
          tau21_noise_qs : list
                          the quantiles of the standardized information flow from noise2 to noise1 for significance testing
    '''
    

# Path: pyleoclim/utils/tsmodel.py
# def ar1_fit_evenly(y):
#     ''' Returns the lag-1 autocorrelation from AR(1) fit.
# 
#     Parameters
#     ----------
#     y : array
#         vector of (float) numbers as a time series
# 
#     Returns
#     -------
#     g : float
#         lag-1 autocorrelation coefficient
# 
#     '''
#     # syntax compatible with statsmodels v0.11.1
#     #ar1_mod = sm.tsa.ARMA(y, (1, 0), missing='drop').fit(trend='nc', disp=0)
#     # syntax compatible with statsmodels v0.12
#     ar1_mod = ARIMA(y, order = (1, 0, 0), missing='drop',trend='ct').fit()
#     g = ar1_mod.params[2]
# 
#     if g > 1:
#         print('Warning: AR(1) fitted autocorrelation greater than 1; setting to 1-eps^{1/4}')
#         eps = np.spacing(1.0)
#         g = 1.0 - eps**(1/4)
# 
#     return g
# 
# Path: pyleoclim/utils/correlation.py
# def sm_ar1_sim(n, p, g, sig):
#     ''' Produce p realizations of an AR1 process of length n with lag-1 autocorrelation g using statsmodels
# 
#     Parameters
#     ----------
# 
#     n : int
#         row dimensions
#     p : int
#         column dimensions
# 
#     g : float
#         lag-1 autocorrelation coefficient
#         
#     sig : float
#         the standard deviation of the original time series
# 
#     Returns
#     -------
# 
#     red : numpy matrix
#         n rows by p columns matrix of an AR1 process
#     '''
#     # specify model parameters (statsmodel wants lag0 coefficents 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 red to achieve the same variance as X
# 
#     red = np.empty(shape=(n, p)) # declare array
# 
#     # simulate AR(1) model for each column
#     for i in np.arange(p):
#         red[:, i] = sm.tsa.arma_generate_sample(ar=ar, ma=ma, nsample=n, burnin=50, scale=sig_n)
# 
#     return red
# 
# def phaseran(recblk, nsurr):
#     ''' Simultaneous phase randomization of a set of time series
#     
#     It creates blocks of surrogate data with the same second order properties as the original
#     time series dataset by transforming the oriinal data into the frequency domain, randomizing the
#     phases simultaneoulsy across the time series and converting the data back into the time domain. 
#     
#     Written by Carlos Gias for MATLAB
# 
#     http://www.mathworks.nl/matlabcentral/fileexchange/32621-phase-randomization/content/phaseran.m
# 
#     Parameters
#     ----------
# 
#     recblk : numpy array
#         2D array , Row: time sample. Column: recording.
#         An odd number of time samples (height) is expected.
#         If that is not the case, recblock is reduced by 1 sample before the surrogate data is created.
#         The class must be double and it must be nonsparse.
#     
#     nsurr : int
#         is the number of image block surrogates that you want to generate.
# 
#     Returns
#     -------
# 
#     surrblk : numpy array
#         3D multidimensional array image block with the surrogate datasets along the third dimension
# 
#     References
#     ----------
# 
#     - Prichard, D., Theiler, J. Generating Surrogate Data for Time Series with Several Simultaneously Measured Variables (1994)
#     Physical Review Letters, Vol 73, Number 7
#     
#     - Carlos Gias (2020). Phase randomization, MATLAB Central File Exchange
#     '''
#     # Get parameters
#     nfrms = recblk.shape[0]
# 
#     if nfrms % 2 == 0:
#         nfrms = nfrms-1
#         recblk = recblk[0:nfrms]
# 
#     len_ser = int((nfrms-1)/2)
#     interv1 = np.arange(1, len_ser+1)
#     interv2 = np.arange(len_ser+1, nfrms)
# 
#     # Fourier transform of the original dataset
#     fft_recblk = np.fft.fft(recblk)
# 
#     surrblk = np.zeros((nfrms, nsurr))
# 
#     #  for k in tqdm(np.arange(nsurr)):
#     for k in np.arange(nsurr):
#         ph_rnd = np.random.rand(len_ser)
# 
#         # Create the random phases for all the time series
#         ph_interv1 = np.exp(2*np.pi*1j*ph_rnd)
#         ph_interv2 = np.conj(np.flipud(ph_interv1))
# 
#         # Randomize all the time series simultaneously
#         fft_recblk_surr = np.copy(fft_recblk)
#         fft_recblk_surr[interv1] = fft_recblk[interv1] * ph_interv1
#         fft_recblk_surr[interv2] = fft_recblk[interv2] * ph_interv2
# 
#         # Inverse transform
#         surrblk[:, k] = np.real(np.fft.ifft(fft_recblk_surr))
# 
#     return surrblk


# Path: pyleoclim/utils/causality.py
import numpy as np
from statsmodels.tsa.stattools import grangercausalitytests
from tqdm import tqdm
from .tsmodel import ar1_fit_evenly
from .correlation import sm_ar1_sim, phaseran
from scipy.stats.mstats import mquantiles

                   **kwargs):
    ''' significance test with surrogates with randomized phases

    Parameters
    ----------

    y1, y2 : array
            vectors of (real) numbers with identical length, no NaNs allowed
    method : {'liang'}
            estimates for the Liang method
    npt : int>=1
         time advance in performing Euler forward differencing,
         e.g., 1, 2. Unless the series are generated with a highly chaotic deterministic system,
         npt=1 should be used.
    nsim : int
          the number of surrogates for significance test
    qs : list
        the quantiles for significance test

    Returns
    -------

    res_dict : dict
        A dictionary with the following information:
          T21_noise_qs : list
                        the quantiles of the information flow from noise2 to noise1 for significance testing
          tau21_noise_qs : list
                          the quantiles of the standardized information flow from noise2 to noise1 for significance testing
    '''
    

    noise1 = phaseran(y1, nsim)