Loading sails/modelfit.py +1 −2 Original line number Diff line number Diff line Loading @@ -192,7 +192,6 @@ class AbstractLinearModel(AbstractAnam): return resid def simulate_data(self, num_samples=1000, num_realisations=1, use_cov=True): num_sources = self.nsignals # Preallocate output Loading sails/stft.py +45 −42 Original line number Diff line number Diff line Loading @@ -45,7 +45,6 @@ import warnings import fractions from dataclasses import dataclass from functools import wraps from copy import deepcopy import numpy as np from scipy import fft as sp_fft Loading Loading @@ -529,11 +528,14 @@ def _proc_spectrum_scaling(pxx, scale, side, mode, nfft): def _proc_get_freq_inds(freqvals, fmin, fmax): """Get indices of specified frequencies within range.""" fidx = (freqvals >= fmin) & \ (freqvals <= fmax) return fidx def _proc_get_time_range(nperseg, nstep, fs, xlen): """Get time indices for a given STFT.""" # Create time window vector noverlap = nperseg - nstep return np.arange(nperseg/2, xlen - nperseg/2 + 1, nperseg - noverlap)/float(fs) Loading @@ -560,7 +562,7 @@ def _proc_nan_freq_range(result, freqvals, fmin=None, fmax=None, axis=-1): """ if fmin >= fmax: logging.error('Selected fmin ({}) is larger than fmax ({})'.format(fmin,)) logging.error('Selected fmin ({}) is larger than fmax ({})'.format(fmin, fmax)) logging.info('Setting data within range {0} - {1} to NaN'.format(fmin, fmax)) fidx = _proc_get_freq_inds(freqvals, fmin, fmax) Loading Loading @@ -597,10 +599,10 @@ def _proc_trim_freq_range(result, freqvals, fmin=None, fmax=None, axis=-1): """ if fmin is None and fmax is None: # Just passing through return freqs, result return freqvals, result if fmin >= fmax: logging.error('Selected fmin ({}) is larger than fmax ({})'.format(fmin,)) logging.error('Selected fmin ({}) is larger than fmax ({})'.format(fmin, fmax)) logging.info('Trimming freq axis to range {0} - {1}'.format(fmin, fmax)) fmin = freqvals[0] if fmin is None else fmin Loading @@ -621,7 +623,7 @@ def _proc_apply_average(p, average, axis=0, keepdims=False): p : array_like Array to compute average across average : {'mean', 'median', 'median_bias', 'min'} Method of Method of averaging axis : int Index of dimension to average across (default = 0) keepdims : bool, optional Loading Loading @@ -652,7 +654,7 @@ def _proc_apply_average(p, average, axis=0, keepdims=False): p = np.nanmedian(p, axis=axis, keepdims=keepdims) / bias elif average == 'min': p = np.nanmin(p, axis=axis, keepdims=keepdims) elif average == None: elif average is None: pass else: msg = "'average' value of '{0}' not recognised - please use 'mean' or 'median'" Loading Loading @@ -1125,6 +1127,7 @@ class GLMIRASAConfig(GLMPeriodogramConfig, IRASAConfig): # These functions take input data, run the option handling and execute whatever # computations are needed @dataclass class SpectrumResult: Loading Loading @@ -1659,6 +1662,7 @@ def _glm_fit_glmtools(pxx, reg_categorical, reg_ztrans, reg_unitmax, DC.add_regressor(name='Constant', rtype='Constant') for key in reg_categorical.keys(): logging.debug('Adding Condition : {0}'.format(key)) # TODO: Need to review this # DC.add_regressor(name=key, rtype='Categorical', datainfo=key, codes=[1]) DC.add_regressor(name=key, rtype='Parametric', datainfo=key, preproc=None) for key in reg_ztrans.keys(): Loading Loading @@ -1800,17 +1804,9 @@ def glm_multitaper(X, reg_categorical=None, reg_ztrans=None, reg_unitmax=None, if axis == -1: axis = X.ndim - 1 if X.ndim != 1 and fit_method in ['pinv', 'lstsq']: msg = "Data input should be vector for 'pinv' and 'lstsq' fits - data shape {0} was passed in" if axis == -1: axis = X.ndim - 1 if X.ndim != 1 and fit_method in ['pinv', 'lstsq']: msg = "Data input should be vector for 'pinv' and 'lstsq' fits - data shape {0} was passed in" logging.error(msg.format(X.shape)) logging.error("Use fit_method='glmtools' for multdimensional data") raise ValueError("Fit methods 'pinv' and 'lstsq' not implemented for multidimensional data") # Set configuration logging.info('Setting config options') config = MultiTaperConfig(X.shape[axis], Loading Loading @@ -1851,28 +1847,36 @@ def glm_multitaper(X, reg_categorical=None, reg_ztrans=None, reg_unitmax=None, # 3. Quick # ISSUES WITH IRASA # 1. Resampling pairs do nothing & geometric mean within pairs does very little, better that each contributes to a robust median # 1. Resampling pairs do nothing & geometric mean within pairs does very little, # better that each contributes to a robust median # 2. Fraction inversion resampling factors are more effective on the high end than the low end. # They systematically leave a residue on the low frequency end of broad peaks # 2b. Though the factors are symmetric in the time domain, oscillatory peaks are irregularly spaced in the frequency domain # 3. IRASA estimates the aperodic component whilst the oscillatory component is a residual and is therefore mixed with noise. # 2b. Though the factors are symmetric in the time domain, # oscillatory peaks are irregularly spaced in the frequency domain # 3. IRASA estimates the aperodic component whilst the oscillatory component is a # residual and is therefore mixed with noise. # 4. Oscillatory, but not aperiodic, component can be negative - neither should ever be negative. # Solutions # 1. Resampling is carried out one at a time and all resamples contribute to median where they benefit robust estimation # 2. Resampling factors are logarithmically spaced - ensuring that oscillatory peaks are (nearly) uniformly spaced in freq-domain # 3. The final aperiodic component is ap(f) = min(ap(f), full(f)) - this forces aperiodic and oscillatory components to be +ve # 1. Resampling is carried out one at a time and all resamples contribute # to median where they benefit robust estimation # 2. Resampling factors are logarithmically spaced - ensuring that oscillatory # peaks are (nearly) uniformly spaced in freq-domain # 3. The final aperiodic component is ap(f) = min(ap(f), full(f)) - # this forces aperiodic and oscillatory components to be +ve # (This constraint should not be enforced for GLM-IRASA) # ISSUE WITH 1/f fitting # 1. Unreasonable to assume that one single slope fits whole spectrum # 1b. This relates to the necessity for a 'knee' in fooof and spectral-plateau as well as lots of empirical results # 1c. Also need to acknoweldge fundamental ambiguity in solution, particularly where there are broad spans of oscillatory peaks # 1c. Also need to acknoweldge fundamental ambiguity in solution, # particularly where there are broad spans of oscillatory peaks # Solutions # 1. 1/f slope is parameterised for purpose such as group/condition comparison or correlation # 2. Estimate slope in flexible manner (IRASA) that allows for some flexiblity in range of slopes (not perfect though) # 3. Compute contrasts or correlations with GLM-Spectrum in freq-by-freq approach and visualise whole spectrum effect within aperiodic component. # 3. Compute contrasts or correlations with GLM-Spectrum in freq-by-freq # approach and visualise whole spectrum effect within aperiodic component. # IDEAS # 1. IRASA with multitaper would also expect slope to be self-similar across smoothness parameters? Loading @@ -1884,7 +1888,6 @@ def glm_multitaper(X, reg_categorical=None, reg_ztrans=None, reg_unitmax=None, @set_verbose def irasa(x, method='original', resample_factors=None, aperiodic_average='median', #bootstrap_osc=None, bs_average='mean', # Periodogram args average='median', # STFT window args Loading sails/_docstring_utils.py +2 −2 File changed.Contains only whitespace changes. Show changes Loading
sails/modelfit.py +1 −2 Original line number Diff line number Diff line Loading @@ -192,7 +192,6 @@ class AbstractLinearModel(AbstractAnam): return resid def simulate_data(self, num_samples=1000, num_realisations=1, use_cov=True): num_sources = self.nsignals # Preallocate output Loading
sails/stft.py +45 −42 Original line number Diff line number Diff line Loading @@ -45,7 +45,6 @@ import warnings import fractions from dataclasses import dataclass from functools import wraps from copy import deepcopy import numpy as np from scipy import fft as sp_fft Loading Loading @@ -529,11 +528,14 @@ def _proc_spectrum_scaling(pxx, scale, side, mode, nfft): def _proc_get_freq_inds(freqvals, fmin, fmax): """Get indices of specified frequencies within range.""" fidx = (freqvals >= fmin) & \ (freqvals <= fmax) return fidx def _proc_get_time_range(nperseg, nstep, fs, xlen): """Get time indices for a given STFT.""" # Create time window vector noverlap = nperseg - nstep return np.arange(nperseg/2, xlen - nperseg/2 + 1, nperseg - noverlap)/float(fs) Loading @@ -560,7 +562,7 @@ def _proc_nan_freq_range(result, freqvals, fmin=None, fmax=None, axis=-1): """ if fmin >= fmax: logging.error('Selected fmin ({}) is larger than fmax ({})'.format(fmin,)) logging.error('Selected fmin ({}) is larger than fmax ({})'.format(fmin, fmax)) logging.info('Setting data within range {0} - {1} to NaN'.format(fmin, fmax)) fidx = _proc_get_freq_inds(freqvals, fmin, fmax) Loading Loading @@ -597,10 +599,10 @@ def _proc_trim_freq_range(result, freqvals, fmin=None, fmax=None, axis=-1): """ if fmin is None and fmax is None: # Just passing through return freqs, result return freqvals, result if fmin >= fmax: logging.error('Selected fmin ({}) is larger than fmax ({})'.format(fmin,)) logging.error('Selected fmin ({}) is larger than fmax ({})'.format(fmin, fmax)) logging.info('Trimming freq axis to range {0} - {1}'.format(fmin, fmax)) fmin = freqvals[0] if fmin is None else fmin Loading @@ -621,7 +623,7 @@ def _proc_apply_average(p, average, axis=0, keepdims=False): p : array_like Array to compute average across average : {'mean', 'median', 'median_bias', 'min'} Method of Method of averaging axis : int Index of dimension to average across (default = 0) keepdims : bool, optional Loading Loading @@ -652,7 +654,7 @@ def _proc_apply_average(p, average, axis=0, keepdims=False): p = np.nanmedian(p, axis=axis, keepdims=keepdims) / bias elif average == 'min': p = np.nanmin(p, axis=axis, keepdims=keepdims) elif average == None: elif average is None: pass else: msg = "'average' value of '{0}' not recognised - please use 'mean' or 'median'" Loading Loading @@ -1125,6 +1127,7 @@ class GLMIRASAConfig(GLMPeriodogramConfig, IRASAConfig): # These functions take input data, run the option handling and execute whatever # computations are needed @dataclass class SpectrumResult: Loading Loading @@ -1659,6 +1662,7 @@ def _glm_fit_glmtools(pxx, reg_categorical, reg_ztrans, reg_unitmax, DC.add_regressor(name='Constant', rtype='Constant') for key in reg_categorical.keys(): logging.debug('Adding Condition : {0}'.format(key)) # TODO: Need to review this # DC.add_regressor(name=key, rtype='Categorical', datainfo=key, codes=[1]) DC.add_regressor(name=key, rtype='Parametric', datainfo=key, preproc=None) for key in reg_ztrans.keys(): Loading Loading @@ -1800,17 +1804,9 @@ def glm_multitaper(X, reg_categorical=None, reg_ztrans=None, reg_unitmax=None, if axis == -1: axis = X.ndim - 1 if X.ndim != 1 and fit_method in ['pinv', 'lstsq']: msg = "Data input should be vector for 'pinv' and 'lstsq' fits - data shape {0} was passed in" if axis == -1: axis = X.ndim - 1 if X.ndim != 1 and fit_method in ['pinv', 'lstsq']: msg = "Data input should be vector for 'pinv' and 'lstsq' fits - data shape {0} was passed in" logging.error(msg.format(X.shape)) logging.error("Use fit_method='glmtools' for multdimensional data") raise ValueError("Fit methods 'pinv' and 'lstsq' not implemented for multidimensional data") # Set configuration logging.info('Setting config options') config = MultiTaperConfig(X.shape[axis], Loading Loading @@ -1851,28 +1847,36 @@ def glm_multitaper(X, reg_categorical=None, reg_ztrans=None, reg_unitmax=None, # 3. Quick # ISSUES WITH IRASA # 1. Resampling pairs do nothing & geometric mean within pairs does very little, better that each contributes to a robust median # 1. Resampling pairs do nothing & geometric mean within pairs does very little, # better that each contributes to a robust median # 2. Fraction inversion resampling factors are more effective on the high end than the low end. # They systematically leave a residue on the low frequency end of broad peaks # 2b. Though the factors are symmetric in the time domain, oscillatory peaks are irregularly spaced in the frequency domain # 3. IRASA estimates the aperodic component whilst the oscillatory component is a residual and is therefore mixed with noise. # 2b. Though the factors are symmetric in the time domain, # oscillatory peaks are irregularly spaced in the frequency domain # 3. IRASA estimates the aperodic component whilst the oscillatory component is a # residual and is therefore mixed with noise. # 4. Oscillatory, but not aperiodic, component can be negative - neither should ever be negative. # Solutions # 1. Resampling is carried out one at a time and all resamples contribute to median where they benefit robust estimation # 2. Resampling factors are logarithmically spaced - ensuring that oscillatory peaks are (nearly) uniformly spaced in freq-domain # 3. The final aperiodic component is ap(f) = min(ap(f), full(f)) - this forces aperiodic and oscillatory components to be +ve # 1. Resampling is carried out one at a time and all resamples contribute # to median where they benefit robust estimation # 2. Resampling factors are logarithmically spaced - ensuring that oscillatory # peaks are (nearly) uniformly spaced in freq-domain # 3. The final aperiodic component is ap(f) = min(ap(f), full(f)) - # this forces aperiodic and oscillatory components to be +ve # (This constraint should not be enforced for GLM-IRASA) # ISSUE WITH 1/f fitting # 1. Unreasonable to assume that one single slope fits whole spectrum # 1b. This relates to the necessity for a 'knee' in fooof and spectral-plateau as well as lots of empirical results # 1c. Also need to acknoweldge fundamental ambiguity in solution, particularly where there are broad spans of oscillatory peaks # 1c. Also need to acknoweldge fundamental ambiguity in solution, # particularly where there are broad spans of oscillatory peaks # Solutions # 1. 1/f slope is parameterised for purpose such as group/condition comparison or correlation # 2. Estimate slope in flexible manner (IRASA) that allows for some flexiblity in range of slopes (not perfect though) # 3. Compute contrasts or correlations with GLM-Spectrum in freq-by-freq approach and visualise whole spectrum effect within aperiodic component. # 3. Compute contrasts or correlations with GLM-Spectrum in freq-by-freq # approach and visualise whole spectrum effect within aperiodic component. # IDEAS # 1. IRASA with multitaper would also expect slope to be self-similar across smoothness parameters? Loading @@ -1884,7 +1888,6 @@ def glm_multitaper(X, reg_categorical=None, reg_ztrans=None, reg_unitmax=None, @set_verbose def irasa(x, method='original', resample_factors=None, aperiodic_average='median', #bootstrap_osc=None, bs_average='mean', # Periodogram args average='median', # STFT window args Loading