Roadmap

[regeirk]

Purpose

I'm thinking about refactoring this module and I invite everyone who is interested in sharing their thoughts.

General goals

1. Easy and user-friendly framework for wavelet analysis and statistical significance tests;
2. Handling of pandas.DataFrame();
3. Return result data class instead of individual vectors and matrices;
4. Add plotting tools;
5. Interfaces for interactive development (e.g. IPython, Jupyter)

[klapo]

Hello @regeirk! Thanks for taking up the tasks of modernizing and maintaining this library. I have some time I can dedicate to helping, although not a lot. I'm splitting my responses up for each point.

[regeirk]

Hi @klapo, thank you for your sharing your thoughts :)

[klapo]

I hope they were useful and not just ramblings!

[klapo]

(1) Yes I fully agree and support this! I built some helper functions to handle the set-up and processing of the pycwt output since I had to do these operations every time.

In regards to sig. testing:
The current wavelet significance testing is somewhat rigid and can lead to misleading results. I show this for cwt here and the xwt and coherence here. My wish list would include:

• Ability to specify the red noise coefficient directly for cwt significance testing (not being able to leads to very misleading results).
• Improving the Monte-Carlo significance test for biwavelet signficance. It is painfully slow (slower than performing significance testing directly using Monte-Carlo!) and does not agree with the actual results from a Monte-Carlo test of the significance.

(2) I think handling step 1 easily addresses step 2, no?

[klapo]

I need to use wavelets again for my current work, so I was thinking of taking a stab at building a class for the pycwt to make my work easier. Do you have any particular items you would like me to include when creating this prototype class?

So far, the class will ingest:

• The actual signal
• Spacing of the signal, dx
• Coordinate to perform along, x
• Mother wavelet
• Octaves (with some reasonable defaults)
• Scales to build an average wavelet power (default = None)
• Keywords for how/whether to normalize or detrend the signal -> potentially move this to a static helper function for the class.
• Keywords specifying whether to normalize the global wavelet power (default = True)
• Keywords specifying whether to rectify (default = True)
• The possibility to specify the data variance (e.g., if the data were normalized previously or derived through a separate method)
• Keywords specifying whether to perform a significance test (default = True)
• Keywords for calculating the significance (i.e., significance level, red noise or white noise coefficients)

and the class would have the following accessible attributes:

• The CWT periods
• coi
• Wavelet power
• Global significance curve
• Global power
• Some variable specifying portions of the wavelet periodogram that were found to be signifiant (i.e., an index)
• Scale averaged power
• Scale averaged significance
• The scale periods

The current work returns all of these quantities as numpy arrays.

[klapo]

I made a (still very messy) driver class with a sklearn-like .fit syntax. You can take a look at the direction I'm going on my branch. Is this close to what you had in mind?

The general workflow would look like the below, which generates a plot that looks like the attached (note that it includes statistical testing).

# Generate the driver object.
cwt = pycwt.wvlt()

# Fit the data from `signal` using some defined normalization procedure
# and specified significance testing.
cwt = cwt.fit(
    signal,
    dx,
    x,
    mother,
    norm_kwargs=norm_kwargs,
    variance=variance,
    sig_kwargs={
        'alpha': alpha,
        'sig_lvl': sig_lvl,
    }
)

# Send it the fit object to a plotter
pycwt.plotters.plot_wv_power(
    cwt,
)

[klapo]

(3) I think this makes a lot of sense and would love to help here. I would caution to make sure that we don't shoehorn a data model onto users that isn't deeply flexible. Maybe as a point of reference we can look at other python libraries with complex fitting operations like arviz and pymc. They use xarray objects under the hood which are able to flexibly handle multidimensional data. Using a data model that can robustly handle multidimensinal data makes a lot of sense to me since wavelets expand data to include at least one new dimension (frequency) and potentially several more (scale averaged quantities).

I have examples of how xarray could be integrated into a pycwt data model if you are interested in pursuing this (note: not a public repository).

[klapo]

(4) Again I absolutely agree. I have some plotting tools already and you are welcome to use any components from them that you think would be useful. I think point 1 would go nicely improving the capabilities of pycwt and make the plotting simpler to call and adapt.

[klapo]

(5) Could you elaborate what your intention is here?

[regeirk]

I thought of adding rich display methods as described here.

[klapo]

Oh ok! Like for displaying the pycwt object. Thanks for clarifying, I thought you meant for graphing.

I think we could circle back to this idea after tackling items 1-3. How does this sound for a rough priority?

[klapo]

We should also consider converting from a travis CI to using github actions and enforcing a code style. I vote for using the black code style.

[regeirk]

After more than a year, I have finally managed to resume the refactoring of this package.

I want to make things more user friendly and use more default parameters. I thought that we should be able to run a basic and simple wavelet analysis using something like this:

import pycwt
tseries = pycwt.load_dataset("sst_nino3")  # Load sample time-series.
wavelet = pycwt.Morlet(w0=6)               # Instantiate mother wavelet.
result = wavelet.run(tseries)              # Run wavelet analysis.
result.plot()

Here, it is expected that the tseries is of type pandas.Series, has a constant sampling frequency dt, no data gaps (NaN values) and that time is given by the series index .

Parameters -- The default scale and significance parameters for this analysis are:

• Starting scale: s0 = 2 * dt
• Spacing between discrete scales: dj = 1 / 12 -- twelve sub-octaves per octave.
• Number of scales: J = 7 / dj -- seven powers of two with dj sub-octaves.
• Lag-1 autocorrelation for red noise: alpha = 1.
• Significance level: significance_level = 0.95.

[klapo]

Hi @regeirk,

I am happy this is back on our radar! Here are my thoughts:

Does it make sense to force a user to use pandas.Series? What about just accepting anything array-like so that someone could e.g. pass in an xarray DataArray and making the user pass in an explicit dt? It would be in line with how other packages require users pass in data to their decomposition (e.g., PyDMD). I think it would work either way, I just personally don't use pandas.Series as my base object of data manipulation.

I would also suggest we use the sklearn .fit syntax instead of .run.

As for the significance testing, I think it can be exceptionally misleading if a user gets out a significance measure that does not represent the actual data they are testing, e.g. if my system has an alpha of 0.2 but I use the defaults I get out a significance test that doesn't make sense. I would argue the user should have to pass the sig testing parameters explicitly to avoid spurious detection.