Merge pull request #296 from fooof-tools/getmodel

[ENH] - Add getter functions for data & model components

Author: TomDonoghue

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.
 #

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.
 

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).
 

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]

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."""