|
1 | | -# FOOOF.jl |
| 1 | +# PyFOOOF.jl |
2 | 2 | Julia interface to [FOOOF](https://github.com/fooof-tools/fooof) |
| 3 | + |
| 4 | +[![Build Status][build-img]][build-url] [![CodeCov][codecov-img]][codecov-url] |
| 5 | + |
| 6 | +[build-img]: https://github.com/beacon-biosignals/PyFOOOF.jl/workflows/CI/badge.svg |
| 7 | +[build-url]: https://github.com/beacon-biosignals/PyFOOOF.jl/actions |
| 8 | +[codecov-img]: https://codecov.io/github/beacon-biosignals/PyFOOOF.jl/badge.svg?branch=main |
| 9 | +[codecov-url]: https://codecov.io/github/beacon-biosignals/PyFOOOF.jl?branch=main |
| 10 | + |
| 11 | + |
| 12 | +## Installation |
| 13 | +This package uses [`PyCall`](https://github.com/JuliaPy/PyCall.jl/) to make |
| 14 | +[FOOOF](https://fooof-tools.github.io/fooof/index.html) available from within Julia. |
| 15 | +Unsurprisingly, FOOOF and its dependencies need to be installed in order for this to work |
| 16 | +and PyFOOOF will attempt to install when the package is built. |
| 17 | + |
| 18 | +By default, this installation happens in the "global" path for the Python used |
| 19 | +by PyCall. If you're using PyCall via its hidden Miniconda install, your own |
| 20 | +Anaconda environment, or a Python virtual environment, this is what you want. |
| 21 | +(The "global" path is sandboxed to the Conda/virtual environment.) If you are |
| 22 | +however using system Python, then you should set `ENV["PIPFLAGS"] = "--user"` |
| 23 | +before `add`ing / `build`ing the package. By default, PyFOOOF will use the latest |
| 24 | +FOOOF 1.x release available on [PyPI](https://pypi.org/project/FOOOF/), but this can also |
| 25 | +be changed via the `ENV["FOOOFVERSION"] = version_number` for your preferred |
| 26 | +`version_number`. |
| 27 | + |
| 28 | +Note that FOOOF uses [`matplotlib`](https://matplotlib.org/) for plotting, but does not install it automatically as a dependency. |
| 29 | +The Julia package [`PyPlot`](https://github.com/JuliaPy/PyPlot.jl), which provides a Julia interface to `matplotlib`, is useful for installing `matplotlib` and manipulating the rendered plots. |
| 30 | + |
| 31 | +FOOOF can also be installed manually ahead of time. |
| 32 | +From the shell, use `python -m pip install fooof` for the latest stable release |
| 33 | +or `python -m pip install fooof==version_number` for a given `version_number`, |
| 34 | +ensuring that `python` is the same one that PyCall is using. Alternatively, |
| 35 | +you can run this from within Julia: |
| 36 | +```julia |
| 37 | +julia> using PyCall |
| 38 | +julia> pip = pyimport("pip"); |
| 39 | + |
| 40 | +julia> pip.main(["install", "fooof==version_number"]); # specific version |
| 41 | +``` |
| 42 | + |
| 43 | +If you do not specify a version via `==version`, then the latest versions will be |
| 44 | +installed. If you wish to upgrade versions, you can use |
| 45 | +`python -m pip install --upgrade fooof` or |
| 46 | +```julia |
| 47 | +julia> using PyCall |
| 48 | + |
| 49 | +julia> pip = pyimport("pip"); |
| 50 | + |
| 51 | +julia> pip.main(["install", "--upgrade", "FOOOF"]); |
| 52 | +``` |
| 53 | + |
| 54 | +You can test your setup with `using PyCall; pyimport("fooof")`. |
| 55 | + |
| 56 | +## Usage |
| 57 | + |
| 58 | +In the same philosophy as PyCall, this allows for the transparent use of |
| 59 | +FOOOF from within Julia. |
| 60 | +The major things the package does are wrap the installation of FOOOF in the |
| 61 | +package `build` step, load all the FOOOF functionality into the module namespace, |
| 62 | +and provide a few accessors. |
| 63 | + |
| 64 | + |
| 65 | +### Exposing FOOOF in Julia |
| 66 | + |
| 67 | +For example, in Python you can create a new FOOOF model like this: |
| 68 | + |
| 69 | +```python |
| 70 | +>>> import fooof |
| 71 | + |
| 72 | +>>> fm = fooof.FOOOF() |
| 73 | +``` |
| 74 | + |
| 75 | +With PyFOOOF, you can do this from within Julia. |
| 76 | + |
| 77 | +```julia |
| 78 | +julia> using PyFOOOF |
| 79 | + |
| 80 | +julia> fm = PyFOOOF.FOOOF(); |
| 81 | +``` |
| 82 | + |
| 83 | +The PyCall infrastructure also means that Python docstrings are available |
| 84 | +in Julia: |
| 85 | + |
| 86 | +```julia |
| 87 | +help?> PyFOOOF.FOOOF |
| 88 | +Model a physiological power spectrum as a combination of aperiodic and periodic components. |
| 89 | + |
| 90 | + WARNING: FOOOF expects frequency and power values in linear space. |
| 91 | + |
| 92 | + Passing in logged frequencies and/or power spectra is not detected, |
| 93 | + and will silently produce incorrect results. |
| 94 | + |
| 95 | + Parameters |
| 96 | + ---------- |
| 97 | + peak_width_limits : tuple of (float, float), optional, default: (0.5, 12.0) |
| 98 | + Limits on possible peak width, in Hz, as (lower_bound, upper_bound). |
| 99 | + max_n_peaks : int, optional, default: inf |
| 100 | + Maximum number of peaks to fit. |
| 101 | + min_peak_height : float, optional, default: 0 |
| 102 | + Absolute threshold for detecting peaks, in units of the input data. |
| 103 | + peak_threshold : float, optional, default: 2.0 |
| 104 | + Relative threshold for detecting peaks, in units of standard deviation of the input data. |
| 105 | + aperiodic_mode : {'fixed', 'knee'} |
| 106 | + Which approach to take for fitting the aperiodic component. |
| 107 | + verbose : bool, optional, default: True |
| 108 | + Verbosity mode. If True, prints out warnings and general status updates. |
| 109 | + |
| 110 | + Attributes |
| 111 | + ---------- |
| 112 | + freqs : 1d array |
| 113 | + Frequency values for the power spectrum. |
| 114 | + power_spectrum : 1d array |
| 115 | + Power values, stored internally in log10 scale. |
| 116 | + freq_range : list of [float, float] |
| 117 | + Frequency range of the power spectrum, as [lowest_freq, highest_freq]. |
| 118 | + freq_res : float |
| 119 | + Frequency resolution of the power spectrum. |
| 120 | + fooofed_spectrum_ : 1d array |
| 121 | + The full model fit of the power spectrum, in log10 scale. |
| 122 | + aperiodic_params_ : 1d array |
| 123 | + Parameters that define the aperiodic fit. As [Offset, (Knee), Exponent]. |
| 124 | + The knee parameter is only included if aperiodic component is fit with a knee. |
| 125 | + peak_params_ : 2d array |
| 126 | + Fitted parameter values for the peaks. Each row is a peak, as [CF, PW, BW]. |
| 127 | + gaussian_params_ : 2d array |
| 128 | + Parameters that define the gaussian fit(s). |
| 129 | + Each row is a gaussian, as [mean, height, standard deviation]. |
| 130 | + r_squared_ : float |
| 131 | + R-squared of the fit between the input power spectrum and the full model fit. |
| 132 | + error_ : float |
| 133 | + Error of the full model fit. |
| 134 | + n_peaks_ : int |
| 135 | + The number of peaks fit in the model. |
| 136 | + has_data : bool |
| 137 | + Whether data is loaded to the object. |
| 138 | + has_model : bool |
| 139 | + Whether model results are available in the object. |
| 140 | + |
| 141 | + Notes |
| 142 | + ----- |
| 143 | + - Commonly used abbreviations used in this module include: |
| 144 | + CF: center frequency, PW: power, BW: Bandwidth, AP: aperiodic |
| 145 | + - Input power spectra must be provided in linear scale. |
| 146 | + Internally they are stored in log10 scale, as this is what the model operates upon. |
| 147 | + - Input power spectra should be smooth, as overly noisy power spectra may lead to bad fits. |
| 148 | + For example, raw FFT inputs are not appropriate. Where possible and appropriate, use |
| 149 | + longer time segments for power spectrum calculation to get smoother power spectra, |
| 150 | + as this will give better model fits. |
| 151 | + - The gaussian params are those that define the gaussian of the fit, where as the peak |
| 152 | + params are a modified version, in which the CF of the peak is the mean of the gaussian, |
| 153 | + the PW of the peak is the height of the gaussian over and above the aperiodic component, |
| 154 | + and the BW of the peak, is 2*std of the gaussian (as 'two sided' bandwidth). |
| 155 | +``` |
| 156 | + |
| 157 | +### Helping with type conversions |
| 158 | + |
| 159 | +PyCall can be rather aggressive in converting standard types, such as |
| 160 | +dictionaries, to their native Julia equivalents, but this can create problems |
| 161 | +due to differences in the way inheritance is traditionally used between |
| 162 | +languages. |
| 163 | +In particular, Julia arrays are converted to NumPy arrays and *not* Python lists. |
| 164 | +This conversion creates problems where FOOOF expects a list and not an array, for example in the `freq_range` keyword argument: |
| 165 | + |
| 166 | +```julia |
| 167 | +julia> fm = FOOOF(; peak_width_limits=[1,8]) |
| 168 | +ERROR: PyError ($(Expr(:escape, :(ccall(#= /home/ubuntu/.julia/packages/PyCall/BD546/src/pyfncall.jl:43 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw))))) <class 'ValueError'> |
| 169 | +ValueError('The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()') |
| 170 | + File "/home/ubuntu/anaconda3/lib/python3.8/site-packages/fooof/objs/fit.py", line 193, in __init__ |
| 171 | + self._reset_internal_settings() |
| 172 | + File "/home/ubuntu/anaconda3/lib/python3.8/site-packages/fooof/objs/fit.py", line 236, in _reset_internal_settings |
| 173 | + if self.peak_width_limits: |
| 174 | +... |
| 175 | +``` |
| 176 | +(The particular problem arises here because FOOOF is depending on the Python's automatic conversion of `None` and empty lists to `False` and non-empty lists to `True`.) |
| 177 | + |
| 178 | +Note that simply wrapping the array as a Python literal (`py"[1,8]"`) does not suffice because this is converted to a Julia vector and thus then to a NumPy array when passed back to Python. Instead, we have to force PyCall to not convert the resulting Python object with the `o` suffix: |
| 179 | + |
| 180 | +```julia |
| 181 | +julia> fm = FOOOF(; peak_width_limits=py"[1,8]"o) |
| 182 | +PyObject <fooof.objs.fit.FOOOF object at 0x7fea38b5b040> |
| 183 | +``` |
| 184 | + |
| 185 | +Another conversion problem arises in cases where nesting lists and eltypes creates problems. |
| 186 | +For example, the [FOOOF tutorial "Tuning and Troubleshooting"](https://fooof-tools.github.io/fooof/auto_tutorials/plot_07-TroubleShooting.html) includes this statement |
| 187 | +```python |
| 188 | +>>> gauss_params = [[10, 1.0, 2.5], [20, 0.8, 2], [32, 0.6, 1]] |
| 189 | +``` |
| 190 | +When executed via PyCall, this results in a Julia `Matrix` with eltype `Real`. |
| 191 | + |
| 192 | +```julia |
| 193 | +julia> py"$([[10, 1.0, 2.5], [20, 0.8, 2], [32, 0.6, 1]])" |
| 194 | +PyObject [array([10. , 1. , 2.5]), array([20. , 0.8, 2. ]), array([32. , 0.6, 1. ])] |
| 195 | +julia> py"$([[10, 1.0, 2.5], [20, 0.8, 2], [32, 0.6, 1]])" |
| 196 | +3-element Vector{Vector{Float64}}: |
| 197 | + [10.0, 1.0, 2.5] |
| 198 | + [20.0, 0.8, 2.0] |
| 199 | + [32.0, 0.6, 1.0] |
| 200 | +``` |
| 201 | +This results in an error: |
| 202 | +```julia |
| 203 | +julia> gen_power_spectrum = fooof.sim.gen.gen_power_spectrum; |
| 204 | + |
| 205 | +julia> f_range = [1, 50]; |
| 206 | + |
| 207 | +julia> ap_params = [20, 2]; |
| 208 | + |
| 209 | +julia> nlv = 0.1; |
| 210 | + |
| 211 | +julia> gauss_params = [[10, 1.0, 2.5], [20, 0.8, 2], [32, 0.6, 1]]; |
| 212 | + |
| 213 | +julia> freqs, spectrum = gen_power_spectrum(f_range, ap_params, gauss_params, nlv) |
| 214 | +ERROR: PyError ($(Expr(:escape, :(ccall(#= /home/ubuntu/.julia/packages/PyCall/BD546/src/pyfncall.jl:43 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw))))) <class 'ValueError'> |
| 215 | +ValueError('operands could not be broadcast together with shapes (99,) (3,) ') |
| 216 | +... |
| 217 | +``` |
| 218 | + |
| 219 | +When the statement is executed in Python (`py"[]"`) and then roundtripped through Julia, PyCall converts the Python return value to a `Matrix`, which works in the subsequent function call: |
| 220 | +```julia |
| 221 | +julia> gauss_params = py"[[10, 1.0, 2.5], [20, 0.8, 2], [32, 0.6, 1]]" |
| 222 | +3×3 Matrix{Real}: |
| 223 | + 10 1.0 2.5 |
| 224 | + 20 0.8 2 |
| 225 | + 32 0.6 1 |
| 226 | +julia> freqs, spectrum = gen_power_spectrum(f_range, ap_params, gauss_params, nlv) |
| 227 | +([1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5 … 45.5, 46.0, 46.5, 47.0, 47.5, 48.0, 48.5, 49.0, 49.5, 50.0], [9.112713501760112e19, 6.707288550822094e19, 3.4304395055235047e19, 1.6048034860916263e19, 1.3121468876633584e19, 9.23648446980319e18, 7.068034503219047e18, 7.474675398285033e18, 5.682794734823231e18, 6.002884162025267e18 … 5.494028369147603e16, 5.044411758605143e16, 4.528833513498138e16, 4.080554951080287e16, 4.064069219484658e16, 3.9731296024126536e16, 3.21719026879766e16, 4.828351597256686e16, 4.441592192173848e16, 4.129641670786365e16]) |
| 228 | +``` |
| 229 | + |
| 230 | +However, naively using a true 2d array (`Matrix`) in Julia also results in error: |
| 231 | +```julia |
| 232 | +julia> gauss_params = [10 1.0 2.5; 20 0.8 2; 32 0.6 1] |
| 233 | +3×3 Matrix{Float64}: |
| 234 | + 10.0 1.0 2.5 |
| 235 | + 20.0 0.8 2.0 |
| 236 | + 32.0 0.6 1.0 |
| 237 | +julia> freqs, spectrum = gen_power_spectrum(f_range, ap_params, gauss_params, nlv) |
| 238 | +ERROR: PyError ($(Expr(:escape, :(ccall(#= /home/ubuntu/.julia/packages/PyCall/BD546/src/pyfncall.jl:43 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw))))) <class 'ValueError'> |
| 239 | +ValueError('operands could not be broadcast together with shapes (99,) (3,) ') |
| 240 | +... |
| 241 | +``` |
| 242 | + |
| 243 | +The problem is in the eltype: we need to force it to `Real` so that the integers are preserved as integers when passed to Python: |
| 244 | + |
| 245 | +```julia |
| 246 | +julia> gauss_params = Real[10 1.0 2.5; 20 0.8 2; 32 0.6 1] |
| 247 | +3×3 Matrix{Real}: |
| 248 | + 10 1.0 2.5 |
| 249 | + 20 0.8 2 |
| 250 | + 32 0.6 1 |
| 251 | +``` |
| 252 | + |
| 253 | +If other automatic type conversions are found to be problematic or there are |
| 254 | +particular FOOOF functions that don't play nice via the default PyCall mechanisms, |
| 255 | +then issues and pull requests are welcome. |
| 256 | + |
| 257 | +Many of these problematic conversions can be fixed with relatively straightforward (and backward compatible) changes to FOOOF; we are in the process of opening PRs for this purpose. |
0 commit comments