merge optimize branch

Author: TomDonoghue

specparam/core/funcs.py

@@ -32,9 +32,7 @@ def gaussian_function(xs, *params):
 
     ys = np.zeros_like(xs)
 
-    for ii in range(0, len(params), 3):
-
-        ctr, hgt, wid = params[ii:ii+3]
+    for ctr, hgt, wid in zip(*[iter(params)] * 3):
 
         ys = ys + hgt * np.exp(-(xs-ctr)**2 / (2*wid**2))
 
@@ -60,11 +58,8 @@ def expo_function(xs, *params):
         Output values for exponential function.
     """
 
-    ys = np.zeros_like(xs)
-
     offset, knee, exp = params
-
-    ys = ys + offset - np.log10(knee + xs**exp)
+    ys = offset - np.log10(knee + xs**exp)
 
     return ys
 
@@ -88,11 +83,8 @@ def expo_nk_function(xs, *params):
         Output values for exponential function, without a knee.
     """
 
-    ys = np.zeros_like(xs)
-
     offset, exp = params
-
-    ys = ys + offset - np.log10(xs**exp)
+    ys = offset - np.log10(xs**exp)
 
     return ys
 
@@ -113,11 +105,8 @@ def linear_function(xs, *params):
         Output values for linear function.
     """
 
-    ys = np.zeros_like(xs)
-
     offset, slope = params
-
-    ys = ys + offset + (xs*slope)
+    ys = offset + (xs*slope)
 
     return ys
 
@@ -138,11 +127,8 @@ def quadratic_function(xs, *params):
         Output values for quadratic function.
     """
 
-    ys = np.zeros_like(xs)
-
     offset, slope, curve = params
-
-    ys = ys + offset + (xs*slope) + ((xs**2)*curve)
+    ys = offset + (xs*slope) + ((xs**2)*curve)
 
     return ys
 

specparam/core/jacobians.py

@@ -0,0 +1,103 @@
+""""Functions for computing Jacobian matrices to be used during fitting.
+
+Notes
+-----
+These functions line up with those in `funcs`.
+The parameters in these functions are labeled {a, b, c, ...}, but follow the order in `funcs`.
+These functions are designed to be passed into `curve_fit` to provide a computed Jacobian.
+"""
+
+import numpy as np
+
+###################################################################################################
+###################################################################################################
+
+## Periodic Jacobian functions
+
+def jacobian_gauss(xs, *params):
+    """Create the Jacobian matrix for the Gaussian function.
+
+    Parameters
+    ----------
+    xs : 1d array
+        Input x-axis values.
+    *params : float
+        Parameters for the function.
+
+    Returns
+    -------
+    jacobian : 2d array
+        Jacobian matrix, with shape [len(xs), n_params].
+    """
+
+    jacobian = np.zeros((len(xs), len(params)))
+
+    for i, (a, b, c) in enumerate(zip(*[iter(params)] * 3)):
+
+        ax = -a + xs
+        ax2 = ax**2
+
+        c2 = c**2
+        c3 = c**3
+
+        exp = np.exp(-ax2 / (2 * c2))
+        exp_b = exp * b
+
+        ii = i * 3
+        jacobian[:, ii] = (exp_b * ax) / c2
+        jacobian[:, ii+1] = exp
+        jacobian[:, ii+2] = (exp_b * ax2) / c3
+
+    return jacobian
+
+
+## Aperiodic Jacobian functions
+
+def jacobian_expo(xs, *params):
+    """Create the Jacobian matrix for the exponential function.
+
+    Parameters
+    ----------
+    xs : 1d array
+        Input x-axis values.
+    *params : float
+        Parameters for the function.
+
+    Returns
+    -------
+    jacobian : 2d array
+        Jacobian matrix, with shape [len(xs), n_params].
+    """
+
+    a, b, c = params
+
+    xs_c = xs**c
+    b_xs_c = xs_c + b
+
+    jacobian = np.ones((len(xs), len(params)))
+    jacobian[:, 1] = -1 / b_xs_c
+    jacobian[:, 2] = -(xs_c * np.log10(xs)) / b_xs_c
+
+    return jacobian
+
+
+def jacobian_expo_nk(xs, *params):
+    """Create the Jacobian matrix for the exponential no-knee function.
+
+    Parameters
+    ----------
+    xs : 1d array
+        Input x-axis values.
+    *params : float
+        Parameters for the function.
+
+    Returns
+    -------
+    jacobian : 2d array
+        Jacobian matrix, with shape [len(xs), n_params].
+    """
+
+    jacobian = np.ones((len(xs), len(params)))
+    jacobian[:, 1] = -np.log10(xs)
+
+    return jacobian

specparam/objs/fit.py

@@ -71,6 +71,7 @@
 from specparam.core.modutils import copy_doc_func_to_method
 from specparam.core.utils import group_three, check_array_dim
 from specparam.core.funcs import gaussian_function, get_ap_func, infer_ap_func
+from specparam.core.jacobians import jacobian_gauss
 from specparam.core.errors import (FitError, NoModelError, DataError,
                                    NoDataError, InconsistentDataError)
 from specparam.core.strings import (gen_settings_str, gen_model_results_str,
@@ -191,12 +192,17 @@ def __init__(self, peak_width_limits=(0.5, 12.0), max_n_peaks=np.inf, min_peak_h
         self._gauss_overlap_thresh = 0.75
         # Parameter bounds for center frequency when fitting gaussians, in terms of +/- std dev
         self._cf_bound = 1.5
-        # The maximum number of calls to the curve fitting function
-        self._maxfev = 5000
         # The error metric to calculate, post model fitting. See `_calc_error` for options
         #   Note: this is for checking error post fitting, not an objective function for fitting
         self._error_metric = 'MAE'
 
+        ## PRIVATE CURVE_FIT SETTINGS
+        # The maximum number of calls to the curve fitting function
+        self._maxfev = 5000
+        # The tolerance setting for curve fitting (see scipy.curve_fit - ftol / xtol / gtol)
+        #   Here reduce tolerance to speed fitting. Set value to 1e-8 to match curve_fit default
+        self._tol = 0.00001
+
         ## RUN MODES
         # Set default debug mode - controls if an error is raised if model fitting is unsuccessful
         self._debug = False
@@ -944,7 +950,9 @@ def _simple_ap_fit(self, freqs, power_spectrum):
                 warnings.simplefilter("ignore")
                 aperiodic_params, _ = curve_fit(get_ap_func(self.aperiodic_mode),
                                                 freqs, power_spectrum, p0=guess,
-                                                maxfev=self._maxfev, bounds=ap_bounds)
+                                                maxfev=self._maxfev, bounds=ap_bounds,
+                                                ftol=self._tol, xtol=self._tol, gtol=self._tol,
+                                                check_finite=False)
         except RuntimeError as excp:
             error_msg = ("Model fitting failed due to not finding parameters in "
                          "the simple aperiodic component fit.")
@@ -1001,7 +1009,9 @@ def _robust_ap_fit(self, freqs, power_spectrum):
                 warnings.simplefilter("ignore")
                 aperiodic_params, _ = curve_fit(get_ap_func(self.aperiodic_mode),
                                                 freqs_ignore, spectrum_ignore, p0=popt,
-                                                maxfev=self._maxfev, bounds=ap_bounds)
+                                                maxfev=self._maxfev, bounds=ap_bounds,
+                                                ftol=self._tol, xtol=self._tol, gtol=self._tol,
+                                                check_finite=False)
         except RuntimeError as excp:
             error_msg = ("Model fitting failed due to not finding "
                          "parameters in the robust aperiodic fit.")
@@ -1147,7 +1157,9 @@ def _fit_peak_guess(self, guess):
         # Fit the peaks
         try:
             gaussian_params, _ = curve_fit(gaussian_function, self.freqs, self._spectrum_flat,
-                                           p0=guess, maxfev=self._maxfev, bounds=gaus_param_bounds)
+                                           p0=guess, maxfev=self._maxfev, bounds=gaus_param_bounds,
+                                           ftol=self._tol, xtol=self._tol, gtol=self._tol,
+                                           check_finite=False, jac=jacobian_gauss)
         except RuntimeError as excp:
             error_msg = ("Model fitting failed due to not finding "
                          "parameters in the peak component fit.")

specparam/plts/spectra.py

@@ -60,8 +60,8 @@ def plot_spectra(freqs, power_spectra, log_freqs=False, log_powers=False, freq_r
         freq_range = np.log10(freq_range) if log_freqs else freq_range
 
     # Make inputs iterable if need to be passed multiple times to plot each spectrum
-    plt_powers = np.reshape(power_spectra, (1, -1)) if np.ndim(power_spectra) == 1 else \
-        power_spectra
+    plt_powers = np.reshape(power_spectra, (1, -1)) if isinstance(freqs, np.ndarray) and \
+        np.ndim(power_spectra) == 1 else power_spectra
     plt_freqs = repeat(freqs) if isinstance(freqs, np.ndarray) and freqs.ndim == 1 else freqs
 
     # Set labels
@@ -131,6 +131,10 @@ def plot_spectra_shading(freqs, power_spectra, shades, shade_colors='r',
                         plot_kwargs.get('log_powers', False))
 
 
+# Alias `plot_spectrum_shading` to `plot_spectra_shading` for backwards compatibility
+plot_spectrum_shading = plot_spectra_shading
+
+
 @savefig
 @style_plot
 @check_dependency(plt, 'matplotlib')

specparam/tests/core/test_jacobians.py

@@ -0,0 +1,33 @@
+"""Tests for fooof.core.jacobians."""
+
+from fooof.core.jacobians import *
+
+###################################################################################################
+###################################################################################################
+
+def test_jacobian_gauss():
+
+    xs = np.arange(1, 100)
+    ctr, hgt, wid = 50, 5, 10
+
+    jacobian = jacobian_gauss(xs, ctr, hgt, wid)
+    assert isinstance(jacobian, np.ndarray)
+    assert jacobian.shape == (len(xs), 3)
+
+def test_jacobian_expo():
+
+    xs = np.arange(1, 100)
+    off, knee, exp = 10, 5, 2
+
+    jacobian = jacobian_expo(xs, off, knee, exp)
+    assert isinstance(jacobian, np.ndarray)
+    assert jacobian.shape == (len(xs), 3)
+
+def test_jacobian_expo_nk():
+
+    xs = np.arange(1, 100)
+    off, exp = 10, 2
+
+    jacobian = jacobian_expo_nk(xs, off, exp)
+    assert isinstance(jacobian, np.ndarray)
+    assert jacobian.shape == (len(xs), 2)

specparam/tests/objs/test_fit.py

@@ -391,7 +391,7 @@ def test_fit_failure():
 
     ## Induce a runtime error, and check it runs through
     tfm = SpectralModel(verbose=False)
-    tfm._maxfev = 5
+    tfm._maxfev = 2
 
     tfm.fit(*sim_power_spectrum([3, 50], [50, 2], [10, 0.5, 2, 20, 0.3, 4]))
 
@@ -417,7 +417,7 @@ def test_debug():
     """Test model object in debug mode, including with fit failures."""
 
     tfm = SpectralModel(verbose=False)
-    tfm._maxfev = 5
+    tfm._maxfev = 2
 
     tfm.set_debug_mode(True)
     assert tfm._debug is True

specparam/tests/plts/test_spectra.py

@@ -15,18 +15,22 @@
 @plot_test
 def test_plot_spectra(tfm, tfg, skip_if_no_mpl):
 
-    # Test with 1d inputs - 1d freq array and list of 1d power spectra
+    # Test with 1d inputs - 1d freq array & list of 1d power spectra
     plot_spectra(tfm.freqs, tfm.power_spectrum,
                  file_path=TEST_PLOTS_PATH, file_name='test_plot_spectra_1d.png')
 
-    # Test with 1d inputs - 1d freq array and list of 1d power spectra
+    # Test with 1d inputs - 1d freq array & list of 1d power spectra
     plot_spectra(tfg.freqs, [tfg.power_spectra[0, :], tfg.power_spectra[1, :]],
                  file_path=TEST_PLOTS_PATH, file_name='test_plot_spectra_list_1d.png')
 
     # Test with multiple freq inputs - list of 1d freq array and list of 1d power spectra
     plot_spectra([tfg.freqs, tfg.freqs], [tfg.power_spectra[0, :], tfg.power_spectra[1, :]],
-                 file_path=TEST_PLOTS_PATH,
-                 file_name='test_plot_spectra_lists_1d.png')
+                 file_path=TEST_PLOTS_PATH, file_name='test_plot_spectra_list_1d_freqs.png')
+
+    # Test with multiple lists - list of 1d freqs & list of 1d power spectra (different f ranges)
+    plot_spectra([tfg.freqs, tfg.freqs[:-5]],
+                 [tfg.power_spectra[0, :], tfg.power_spectra[1, :-5]],
+                 file_path=TEST_PLOTS_PATH, file_name='test_plot_spectra_lists_1d.png')
 
     # Test with 2d array inputs
     plot_spectra(np.vstack([tfg.freqs, tfg.freqs]),