Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

uar1 surrogates #537

Merged
merged 10 commits into from
Apr 22, 2024
19 changes: 7 additions & 12 deletions pyleoclim/core/coherence.py
Original file line number Diff line number Diff line change
Expand Up @@ -692,17 +692,17 @@ def signif_test(self, number=200, method='ar1sim', seed=None, qs=[0.95], setting

coh.signif_test(method='phaseran').plot()
'''

from ..core.surrogateseries import SurrogateSeries

if number == 0:
return self

new = self.copy()
surr1 = self.timeseries1.surrogates(
number=number, seed=seed, method=method, settings=settings
)
surr2 = self.timeseries2.surrogates(
number=number, seed=seed, method=method, settings=settings
)

surr1 = SurrogateSeries(method=method,number=number, seed=seed)
surr1.from_series(self.timeseries1)
surr2 = SurrogateSeries(method=method,number=number, seed=seed)
surr2.from_series(self.timeseries2)

# adjust time axis

Expand Down Expand Up @@ -730,11 +730,6 @@ def signif_test(self, number=200, method='ar1sim', seed=None, qs=[0.95], setting
wtc_qs = np.ndarray(shape=(nq, nf, nt))
xwt_qs = np.empty_like(wtc_qs)

# for i in range(nf):
# for j in range(nt):
# wtc_qs[:,i,j] = mquantiles(wtcs[:,i,j], qs)
# xwt_qs[:,i,j] = mquantiles(xwts[:,i,j], qs)

# extract quantiles and reshape
wtc_qs = mquantiles(wtcs_r, qs, axis=0)
wtc_qs = np.reshape(wtc_qs, (nq, nf, nt))
Expand Down
6 changes: 3 additions & 3 deletions pyleoclim/core/ensembleseries.py
Original file line number Diff line number Diff line change
Expand Up @@ -22,8 +22,6 @@
import matplotlib.transforms as transforms
import matplotlib as mpl
from tqdm import tqdm
#import warnings
#from scipy.stats.mstats import mquantiles

class EnsembleSeries(MultipleSeries):
''' EnsembleSeries object
Expand Down Expand Up @@ -524,6 +522,8 @@ def plot_traces(self, figsize=[10, 4], xlabel=None, ylabel=None, title=None, num
'''
savefig_settings = {} if savefig_settings is None else savefig_settings.copy()
lgd_kwargs = {} if lgd_kwargs is None else lgd_kwargs.copy()

num_traces = min(num_traces, len(self.series_list)) # restrict to the smaller of the two

# generate default axis labels
time_label, value_label = self.make_labels()
Expand Down Expand Up @@ -972,7 +972,7 @@ def histplot(self, figsize=[10, 4], title=None, savefig_settings=None,
ax=None, ylabel='KDE', vertical=False, edgecolor='w', **plot_kwargs):
""" Plots the distribution of the timeseries across ensembles

Reuses seaborn [histplot](https://seaborn.pydata.org/generated/seaborn.histplot.html) function.
Reuses the seaborn [histplot](https://seaborn.pydata.org/generated/seaborn.histplot.html) function.

Parameters
----------
Expand Down
34 changes: 4 additions & 30 deletions pyleoclim/core/psds.py
Original file line number Diff line number Diff line change
@@ -1,10 +1,6 @@
from ..utils import plotting
from ..utils import wavelet as waveutils
from ..utils import spectral as specutils
#from ..utils import tsbase

#from ..core.multiplepsd import *
#import ..core.multiplepsd


import matplotlib.pyplot as plt
Expand All @@ -16,28 +12,6 @@

from matplotlib.ticker import ScalarFormatter, FormatStrFormatter


# def infer_period_unit_from_time_unit(time_unit):
# ''' infer a period unit based on the given time unit

# '''
# if time_unit is None:
# period_unit = None
# else:
# unit_group = lipdutils.timeUnitsCheck(time_unit)
# if unit_group != 'unknown':
# if unit_group == 'kage_units':
# period_unit = 'kyrs'
# else:
# period_unit = 'yrs'
# else:
# if time_unit[-1] == 's':
# period_unit = time_unit
# else:
# period_unit = f'{time_unit}s'

# return period_unit

def infer_period_unit_from_time_unit(time_unit):
''' infer a period unit based on the given time unit

Expand Down Expand Up @@ -273,6 +247,7 @@ def signif_test(self, method='ar1sim', number=None, seed=None, qs=[0.95],
pyleoclim.core.series.Series.wavelet : Performs wavelet analysis on Pyleoclim Series

'''
from ..core.surrogateseries import SurrogateSeries

if self.spec_method == 'wwz' and method == 'ar1asym':
raise ValueError('Asymptotic solution is not supported for the wwz method')
Expand Down Expand Up @@ -300,10 +275,9 @@ def signif_test(self, method='ar1sim', number=None, seed=None, qs=[0.95],
return self

new = self.copy()
surr = self.timeseries.surrogates(
number=number, seed=seed, method=method, settings=settings
)

surr = SurrogateSeries(method=method,number=number, seed=seed)
surr.from_series(self.timeseries)

if signif_scals:
surr_psd = surr.spectral(
method=self.spec_method, settings=self.spec_args, scalogram_list=signif_scals
Expand Down
10 changes: 6 additions & 4 deletions pyleoclim/core/scalograms.py
Original file line number Diff line number Diff line change
Expand Up @@ -482,6 +482,7 @@ def signif_test(self, method='ar1sim', number=None, seed=None, qs=[0.95],
scalogram = ts.wavelet().signif_test(number=2, export_scal=True)

'''

if self.wave_method == 'wwz' and method == 'ar1asym':
raise ValueError('Asymptotic solution is not supported for the wwz method')

Expand All @@ -490,6 +491,7 @@ def signif_test(self, method='ar1sim', number=None, seed=None, qs=[0.95],
raise ValueError(f"The available methods are: {acceptable_methods}")

if method in ['ar1sim','uar1'] :
from ..core.surrogateseries import SurrogateSeries

if hasattr(self,'signif_scals'):
signif_scals = self.signif_scals
Expand Down Expand Up @@ -517,14 +519,14 @@ def signif_test(self, method='ar1sim', number=None, seed=None, qs=[0.95],
number -= len(scalogram_list)
surr_scal_tmp = []
surr_scal_tmp.extend(scalogram_list)
surr = self.timeseries.surrogates(number=number, seed=seed,
method=method, settings=settings)
surr = SurrogateSeries(method=method,number=number, seed=seed)
surr.from_series(self.timeseries)

surr_scal_tmp.extend(surr.wavelet(method=self.wave_method, settings=self.wave_args).scalogram_list)
surr_scal = MultipleScalogram(scalogram_list=surr_scal_tmp)
else:
surr = self.timeseries.surrogates(number=number, seed=seed,
method=method, settings=settings)
surr = SurrogateSeries(method=method,number=number, seed=seed)
surr.from_series(self.timeseries)
surr_scal = surr.wavelet(method=self.wave_method, settings=self.wave_args)

if type(qs) is not list:
Expand Down
Loading
Loading