diff --git a/examples/models/plot_data_components.py b/examples/models/plot_data_components.py index 768434e4..ba0e72c4 100644 --- a/examples/models/plot_data_components.py +++ b/examples/models/plot_data_components.py @@ -29,8 +29,8 @@ fm.fit(freqs, powers) ################################################################################################### -# Data Components -# ~~~~~~~~~~~~~~~ +# Data & Model Components +# ----------------------- # # The model fit process includes procedures for isolating aperiodic and periodic components in # the data, fitting each of these components separately, and then combining the model components @@ -39,8 +39,13 @@ # In doing this process, the model fit procedure computes and stores isolated data components, # which are available in the model. # -# Before diving into the isolated data components, let's check the data (`power_spectrum`) -# and full model fit of a model object (`fooofed_spectrum`). + +################################################################################################### +# Full Data & Model Components +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Before diving into the isolated data components, let's check 'full' components, including +# the data (`power_spectrum`) and full model fit of a model object (`fooofed_spectrum`). # ################################################################################################### @@ -53,25 +58,39 @@ # Plot the power spectrum model from the object plot_spectra(fm.freqs, fm.fooofed_spectrum_, color='red') +################################################################################################### +# Isolated Components +# ------------------- +# +# As well as the 'full' data & model components above, the model fitting procedure includes +# steps that result in isolated periodic and aperiodic components, in both the +# data and model. These isolated components are stored internally in the model. +# +# To access these components, we can use the following `getter` methods: +# +# - :meth:`~fooof.FOOOF.get_data`: allows for accessing data components +# - :meth:`~fooof.FOOOF.get_model`: allows for accessing model components +# + ################################################################################################### # Aperiodic Component # ~~~~~~~~~~~~~~~~~~~ # # To fit the aperiodic component, the model fit procedure includes a peak removal process. # -# The resulting 'peak-removed' data component is stored in the model object, in the -# `_spectrum_peak_rm` attribute. +# The resulting 'peak-removed' data component is stored in the model object, as well as the +# isolated aperiodic component model fit. # ################################################################################################### # Plot the peak removed spectrum data component -plot_spectra(fm.freqs, fm._spectrum_peak_rm, color='black') +plot_spectra(fm.freqs, fm.get_data('aperiodic'), color='black') ################################################################################################### # Plot the peak removed spectrum, with the model aperiodic fit -plot_spectra(fm.freqs, [fm._spectrum_peak_rm, fm._ap_fit], +plot_spectra(fm.freqs, [fm.get_data('aperiodic'), fm.get_model('aperiodic')], colors=['black', 'blue'], linestyle=['-', '--']) ################################################################################################### @@ -81,19 +100,20 @@ # To fit the periodic component, the model fit procedure removes the fit peaks from the power # spectrum. # -# The resulting 'flattened' data component is stored in the model object, in the -# `_spectrum_flat` attribute. +# The resulting 'flattened' data component is stored in the model object, as well as the +# isolated periodic component model fit. # ################################################################################################### # Plot the flattened spectrum data component -plot_spectra(fm.freqs, fm._spectrum_flat, color='black') +plot_spectra(fm.freqs, fm.get_data('peak'), color='black') ################################################################################################### # Plot the flattened spectrum data with the model peak fit -plot_spectra(fm.freqs, [fm._spectrum_flat, fm._peak_fit], colors=['black', 'green']) +plot_spectra(fm.freqs, [fm.get_data('peak'), fm.get_model('peak')], + colors=['black', 'green'], linestyle=['-', '--']) ################################################################################################### # Full Model Fit @@ -106,18 +126,88 @@ ################################################################################################### # Plot the full model fit, as the combination of the aperiodic and peak model components -plot_spectra(fm.freqs, [fm._ap_fit + fm._peak_fit], color='red') +plot_spectra(fm.freqs, [fm.get_model('aperiodic') + fm.get_model('peak')], color='red') ################################################################################################### -# Notes on Analyzing Data Components -# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# Linear vs Log Spacing +# --------------------- # # The above shows data components as they are available on the model object, and used in -# the fitting process. Some analyses may aim to use these isolated components to compute -# certain measures of interest on the data. Note that these data components are stored in -# 'private' attributes (indicated by a leading underscore), meaning in normal function they -# are not expected to be accessed by the user, but as we've seen above they can still be accessed. -# However, analyses derived from these isolated data components is not currently officially -# supported by the module, and so users who wish to do so should consider the benefits and -# limitations of any such analyses. +# the fitting process - notable, in log10 spacing. +# +# Some analyses may aim to use these isolated components to compute certain measures of +# interest on the data. However, when doing so, one may often want the linear power +# representations of these components. +# +# Both the `get_data` and `get_model` methods accept a 'space' argument, whereby the user +# can specify whether the return the components in log10 or linear spacing. + + + +################################################################################################### +# Aperiodic Components in Linear Space +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# First we can examine the aperiodic data & model components, in linear space. +# + +################################################################################################### + +# Plot the peak removed spectrum, with the model aperiodic fit +plot_spectra(fm.freqs, [fm.get_data('aperiodic', 'linear'), fm.get_model('aperiodic', 'linear')], + colors=['black', 'blue'], linestyle=['-', '--']) + +################################################################################################### +# Peak Component in Linear Space +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Next, we can examine the peak data & model components, in linear space. +# + +################################################################################################### + +# Plot the flattened spectrum data with the model peak fit +plot_spectra(fm.freqs, [fm.get_data('peak', 'linear'), fm.get_model('peak', 'linear')], + colors=['black', 'green'], linestyle=['-', '--']) + +################################################################################################### +# Linear Space Additive Model +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Note that specifying 'linear' does not simply unlog the data components to return them +# in linear space, but instead defines the space of the additive data definition such that +# `power_spectrum = aperiodic_component + peak_component` (for data and/or model). +# +# We can see this by plotting the linear space data (or model) with the corresponding +# aperiodic and periodic components summed together. Note that if you simply unlog +# the components and sum them, they does not add up to reflecting the full data / model. +# + +################################################################################################### + +# Plot the linear data, showing the combination of peak + aperiodic matches the full data +plot_spectra(fm.freqs, + [fm.get_data('full', 'linear'), + fm.get_data('aperiodic', 'linear') + fm.get_data('peak', 'linear')], + linestyle=['-', 'dashed'], colors=['black', 'red'], alpha=[0.3, 0.75]) + +################################################################################################### + +# Plot the linear model, showing the combination of peak + aperiodic matches the full model +plot_spectra(fm.freqs, + [fm.get_model('full', 'linear'), + fm.get_model('aperiodic', 'linear') + fm.get_model('peak', 'linear')], + linestyle=['-', 'dashed'], colors=['black', 'red'], alpha=[0.3, 0.75]) + +################################################################################################### +# Notes on Analyzing Data & Model Components +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# The functionality here allows for accessing the model components in log space (as used by +# the model for fitting), as well as recomputing in linear space. +# +# If you are aiming to analyze these components, it is important to consider which version of +# the data you should analyze for the question at hand, as there are key differences to the +# different representations. Users who wish to do so post-hoc analyses of these data and model +# components should consider the benefits and limitations the different representations. # diff --git a/fooof/core/utils.py b/fooof/core/utils.py index 3b50e44c..29dd2103 100644 --- a/fooof/core/utils.py +++ b/fooof/core/utils.py @@ -8,6 +8,25 @@ ################################################################################################### ################################################################################################### +def unlog(arr, base=10): + """Helper function to unlog an array. + + Parameters + ---------- + arr : ndarray + Array. + base : float + Base of the log to undo. + + Returns + ------- + ndarray + Unlogged array. + """ + + return np.power(base, arr) + + def group_three(vec): """Group an array of values into threes. diff --git a/fooof/objs/fit.py b/fooof/objs/fit.py index 1648192a..89932505 100644 --- a/fooof/objs/fit.py +++ b/fooof/objs/fit.py @@ -63,6 +63,7 @@ from numpy.linalg import LinAlgError from scipy.optimize import curve_fit +from fooof.core.utils import unlog from fooof.core.items import OBJ_DESC from fooof.core.info import get_indices from fooof.core.io import save_fm, load_json @@ -596,6 +597,95 @@ def get_meta_data(self): for key in OBJ_DESC['meta_data']}) + def get_data(self, component='full', space='log'): + """Get a data component. + + Parameters + ---------- + component : {'full', 'aperiodic', 'peak'} + Which data component to return. + 'full' - full power spectrum + 'aperiodic' - isolated aperiodic data component + 'peak' - isolated peak data component + space : {'log', 'linear'} + Which space to return the data component in. + 'log' - returns in log10 space. + 'linear' - returns in linear space. + + Returns + ------- + output : 1d array + Specified data component, in specified spacing. + + Notes + ----- + The 'space' parameter doesn't just define the spacing of the data component + values, but rather defines the space of the additive data definition such that + `power_spectrum = aperiodic_component + peak_component`. + With space set as 'log', this combination holds in log space. + With space set as 'linear', this combination holds in linear space. + """ + + assert space in ['linear', 'log'], "Input for 'space' invalid." + + if component == 'full': + output = self.power_spectrum if space == 'log' else unlog(self.power_spectrum) + elif component == 'aperiodic': + output = self._spectrum_peak_rm if space == 'log' else \ + unlog(self.power_spectrum) / unlog(self._peak_fit) + elif component == 'peak': + output = self._spectrum_flat if space == 'log' else \ + unlog(self.power_spectrum) - unlog(self._ap_fit) + else: + raise ValueError('Input for component invalid.') + + return output + + + def get_model(self, component='full', space='log'): + """Get a model component. + + Parameters + ---------- + component : {'full', 'aperiodic', 'peak'} + Which model component to return. + 'full' - full model + 'aperiodic' - isolated aperiodic model component + 'peak' - isolated peak model component + space : {'log', 'linear'} + Which space to return the model component in. + 'log' - returns in log10 space. + 'linear' - returns in linear space. + + Returns + ------- + output : 1d array + Specified model component, in specified spacing. + + Notes + ----- + The 'space' parameter doesn't just define the spacing of the model component + values, but rather defines the space of the additive model such that + `model = aperiodic_component + peak_component`. + With space set as 'log', this combination holds in log space. + With space set as 'linear', this combination holds in linear space. + """ + + assert space in ['linear', 'log'], "Input for 'space' invalid." + + if component == 'full': + output = self.fooofed_spectrum_ if space == 'log' else unlog(self.fooofed_spectrum_) + elif component == 'aperiodic': + output = self._ap_fit if space == 'log' else unlog(self._ap_fit) + elif component == 'peak': + output = self._peak_fit if space == 'log' else \ + unlog(self.fooofed_spectrum_) - unlog(self._ap_fit) + else: + raise ValueError('Input for component invalid.') + + return output + + def get_params(self, name, col=None): """Return model fit parameters for specified feature(s). diff --git a/fooof/tests/core/test_utils.py b/fooof/tests/core/test_utils.py index 2dbce06d..01115eaa 100644 --- a/fooof/tests/core/test_utils.py +++ b/fooof/tests/core/test_utils.py @@ -12,6 +12,13 @@ ################################################################################################### ################################################################################################### +def test_unlog(): + + orig = np.array([1, 2, 3, 4]) + logged = np.log10(orig) + unlogged = unlog(logged) + assert np.array_equal(orig, unlogged) + def test_group_three(): dat = [0, 1, 2, 3, 4, 5] diff --git a/fooof/tests/objs/test_fit.py b/fooof/tests/objs/test_fit.py index 4add3fe7..23a49432 100644 --- a/fooof/tests/objs/test_fit.py +++ b/fooof/tests/objs/test_fit.py @@ -311,6 +311,17 @@ def test_obj_gets(tfm): results = tfm.get_results() assert isinstance(results, FOOOFResults) +def test_get_components(tfm): + + # Make sure test object has been fit + tfm.fit() + + # Test get data & model components + for comp in ['full', 'aperiodic', 'peak']: + for space in ['log', 'linear']: + assert isinstance(tfm.get_data(comp, space), np.ndarray) + assert isinstance(tfm.get_model(comp, space), np.ndarray) + def test_get_params(tfm): """Test the get_params method."""