From 1541ef8eeafc473ec95f126ae850bb71431bf3f3 Mon Sep 17 00:00:00 2001 From: Christian Kothe Date: Thu, 31 Oct 2013 01:27:19 -0700 Subject: [PATCH] Minor additions and improvements. --- RELEASE NOTES.TXT | 30 ++++++++++++++++++++--- code/dataset_editing/set_insert_markers.m | 5 +++- code/filters/flt_ica.m | 9 +++++-- code/helpers/hlp_diskcache.m | 11 +++++---- code/misc/window_func.m | 13 +++++----- code/paradigms/ParadigmMTCSP.m | 10 +++++--- code/utils/utl_calc_crossspec.m | 4 +-- code/utils/utl_partition_bundle.m | 2 +- dependencies/fast_decode.m | 1 + dependencies/fast_encode.m | 1 + 10 files changed, 62 insertions(+), 24 deletions(-) diff --git a/RELEASE NOTES.TXT b/RELEASE NOTES.TXT index 67678627..50087f5c 100644 --- a/RELEASE NOTES.TXT +++ b/RELEASE NOTES.TXT @@ -99,6 +99,27 @@ If you are making extensive use of the underlying EEGLAB tools, please consider +=== Changes since 1.1-beta === + +Major New Features + +Minor New Features +- set_insert_markers: now supports negative upper bounds in the absoluterange segment specification +- window_func: now supports some more alternative names for window functions and one or another func +- hlp_diskcache: now has a tad less overhead +- ParadigmMTCSP: now supports LogTransform option, also a bit more efficient +- flt_ica: now allows to retain channel labels when TransformData is set + +Major Changes + +Minor Changes +- partitioning of continuous data during a parameter search will now cache the data set rather than + reload it every time + +Serious Fixes + +Minor Fixes +- flt_ica, robust_sphere option: was calling an obsolete function === Changes since 1.01 === @@ -252,7 +273,7 @@ Minor New Features listing data sets to process (using rdir() rather than dir()), such as /data/mystudy/**/*.set - hlp_serialize: now serializes gpuArray data (as double arrays) -Major Changes: +Major Changes - arg_tovals: now by default sets the arg_direct flag to false rather than true, meaning that some existing functions taking outputs generated via arg_tovals will likely run the slow path unless the flag is set to true explicitly (e.g., via arg_setdirect) @@ -262,7 +283,7 @@ Major Changes: - some methods are now marked as “experimental” and not shown in the GUI by default (can be toggled at startup or in the bcilab_config.m) -Minor Changes: +Minor Changes - utl_fileinfo: now returns the version identifier string for the given file if one is present as the first argument instead of the code hash (otherwise returns the hash as before); the 5th output always contains the code hash @@ -300,13 +321,14 @@ Minor Changes: - bci_batchtrain: was previously running locally by default; now the default matches the current global setting - fit_eeg_distribution: now by default fits an ensemble of truncated distributions and returns a quantile thereof - flt_clean_windows: updated defaults in accordance with fit_eeg_distribution +- ParadigmDALERP: now uses a FIR filter by default instead of IIR -Serious Fixes: +Serious Fixes - ml_calcloss: fixed calculation of FP and FN (Dan Roberts) - utl_crossval: calculated mean AUC loss incorrectly, usually pessimistically (per-fold AUC’s were correct) -Minor Fixes: +Minor Fixes - arg_subswitch: fixed an error message - arg_subtoggle: now properly supports ‘off’ / ‘on’ strings in the defaults - env_showmenu: now highlights the menu when it’s already open diff --git a/code/dataset_editing/set_insert_markers.m b/code/dataset_editing/set_insert_markers.m index 6f19219d..0765cd1f 100644 --- a/code/dataset_editing/set_insert_markers.m +++ b/code/dataset_editing/set_insert_markers.m @@ -129,7 +129,7 @@ arg_subswitch({'segmentspec','SegmentSpec','segment'},'absoluterange', ... {'absoluterange',{ ... arg({'lo','BeginOffset'},0,[],'Lower bound of insertion interval. Events will be inserted beginning from this time point, in seconds.'), ... - arg({'hi','EndOffset'},Inf,[],'Upper bound of insertion interval. Events will be inserted beginning from this time point, in seconds.')}, ... + arg({'hi','EndOffset'},Inf,[],'Upper bound of insertion interval. Events will be inserted up to this time point, in seconds. If this is negative, it counts from the end of the recording.')}, ... 'relativerange',{ ... arg({'event','EventType'},'event1',[],'Reference event type. New events will be inserted in a range around each event of this type.'), ... arg({'lo','BeginOffset'},-0.5,[],'Lower bound relative to event. This is the lower boundary of insertion intervals relative to the reference events. In seconds.'), ... @@ -179,6 +179,9 @@ opts.segmentspec.lo = signal.xmin*signal.srate; end if opts.segmentspec.hi == Inf opts.segmentspec.hi = signal.xmax*signal.srate; end + if opts.segmentspec.hi < 0 + opts.segmentspec.hi = signal.xmax*signal.srate - opts.segmentspec.hi; end + % inject using absolute latencies if isempty(opts.event) error('an event type must be specified'); end diff --git a/code/filters/flt_ica.m b/code/filters/flt_ica.m index 4ffdcb99..f390e302 100644 --- a/code/filters/flt_ica.m +++ b/code/filters/flt_ica.m @@ -346,6 +346,7 @@ }, 'ICA variant. AMICA is the highest quality (but slowest, except if run on a cluster), Infomax is second-highest quality, FastICA is fastest (but can fail to converge and gives poorer results), KernelICA is experimental.'), ... arg_sub({'data_cleaning','DataCleaning','CleaningLevel','clean'},{}, @flt_clean_settings,'Optional data cleaning prior to running an ICA. The computed ICA solution will be applied to the original uncleaned data.'), ... arg({'do_transform','TransformData','transform'},false,[],'Transform the data rather than annotate. By default, ICA decompositions are added as annotations to the data set.'),... + arg({'retain_labels','RetainLabels'},true,[],'Retain labels when transforming. If this is false the channel labels will be replaced by 1:k for k components, if TransformData is checked.'),... arg({'clear_after_trans','ClearAfterTransform'},true,[],'Clear .icaweights after transform. This is so that later functions do not attempt to transform the already transformed data.'),... arg({'do_calcact','CalculateActivation'},false,[],'Calculate component activations. If true, the .icaact field will be populated.'),... arg({'cleaned_data','OutputCleanedData'},false,[],'Emit cleaned data. Whether the cleaned data, instead of the original data should be output (note: this is not applicable for online use, since most cleaning filters cannot be run online).'),... @@ -916,7 +917,7 @@ pre.icasphere = inv(sqrtm(cov(pre.data'))); pre.icaweights = eye(size(pre.data,1)); case 'robust_sphere' - pre.icasphere = inv(sqrtm(hlp_diskcache('icaweights',@cov_robust,pre.data'))); + pre.icasphere = inv(sqrtm(hlp_diskcache('icaweights',@cov_blockgeom,pre.data'))); pre.icaweights = eye(size(pre.data,1)); otherwise % let pop_runica handle all the rest @@ -989,7 +990,11 @@ if isfield(signal.etc,'amica') && size(signal.etc.amica.W,3) > 1 warn_once('Note: The signal will only be transformed according to the 1st amica model.'); end signal.data = (signal.icaweights*signal.icasphere)*signal.data(signal.icachansind,:); - signal.chanlocs = struct('labels',cellfun(@num2str,num2cell(1:length(signal.icachansind),1),'UniformOutput',false)); + if retain_labels && nnz(signal.icachansind) == size(signal.data,1) + signal.chanlocs = signal.chanlocs(signal.icachansind); + else + signal.chanlocs = struct('labels',cellfun(@num2str,num2cell(1:length(signal.icachansind),1),'UniformOutput',false)); + end signal.nbchan = size(signal.data,1); if clear_after_trans signal.icaweights = []; diff --git a/code/helpers/hlp_diskcache.m b/code/helpers/hlp_diskcache.m index e5ad9848..59683f24 100644 --- a/code/helpers/hlp_diskcache.m +++ b/code/helpers/hlp_diskcache.m @@ -150,6 +150,9 @@ archive_version = 1.0; % version of the archive format persistent settings; +persistent have_translatepath; +if isempty(have_translatepath) + have_translatepath = exist('env_translatepath','file'); end % parse options if iscell(options) || isstruct(options) @@ -180,7 +183,7 @@ end % make path platform-specific -if ~exist('env_translatepath','file') +if ~have_translatepath options.folder = strrep(strrep(options.folder,'\',filesep),'/',filesep); else options.folder = env_translatepath(options.folder); @@ -384,11 +387,9 @@ if isfield(x,'tracking') && isfield(x.tracking,'expression') x = trim_expression(x.tracking.expression); elseif iscell(x) - for i=1:length(x) - x{i} = trim_expression(x{i}); end + x = cellfun(@trim_expression,x,'UniformOutput',false); elseif isfield(x,{'head','parts'}) - for i=1:length(x.parts) - x.parts{i} = trim_expression(x.parts{i}); end + x.parts = trim_expression(x.parts); end diff --git a/code/misc/window_func.m b/code/misc/window_func.m index b3189a83..125a749c 100644 --- a/code/misc/window_func.m +++ b/code/misc/window_func.m @@ -41,20 +41,21 @@ % write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 % USA - p = (0:(m-1))/(m-1); switch name case 'bartlett' w = 1 - abs(((0:(m-1)) - (m-1)/2)/((m-1)/2)); - case 'barthann' + case {'barthann','barthannwin'} w = 0.62 - 0.48*abs(p-0.5) - 0.38*cos(2*pi*p); case 'blackman' w = 0.42-0.5*cos(2*pi*p) + 0.08*cos(4*pi*p); case 'blackmanharris' w = 0.35875 - 0.48829*cos(2*pi*p) + 0.14128*cos(4*pi*p) - 0.01168*cos(6*pi*p); - case 'flattop' + case {'bohman','bohmanwin'} + w = (1-abs(p*2-1)).*cos(pi*abs(p*2-1)) + (1/pi)*sin(pi*abs(p*2-1)); + case {'flattop','flattopwin'} w = 0.2157 - 0.4163*cos(2*pi*p) + 0.2783*cos(4*pi*p) - 0.0837*cos(6*pi*p) + 0.0060*cos(8*pi*p); - case 'gauss' + case {'gauss','gausswin'} if nargin < 3 param = 2.5; end w = exp(-0.5*(param*2*(p-0.5)).^2); @@ -68,9 +69,9 @@ w = besseli(0,param*sqrt(1-(2*p-1).^2))/besseli(0,param); case 'lanczos' w = sin(pi*(2*p-1))./(pi*(2*p-1)); w(isnan(w)) = 1; - case 'nuttall' + case {'nuttall','nuttallwin'} w = 0.3635819 - 0.4891775*cos(2*pi*p) + 0.1365995*cos(4*pi*p) - 0.0106411*cos(6*pi*p); - case 'rect' + case {'rect','rectwin'} w = ones(1,m); case 'triang' w = 1 - abs(((0:(m-1)) - (m-1)/2)/((m+1)/2)); diff --git a/code/paradigms/ParadigmMTCSP.m b/code/paradigms/ParadigmMTCSP.m index 7ad64924..5785ee23 100644 --- a/code/paradigms/ParadigmMTCSP.m +++ b/code/paradigms/ParadigmMTCSP.m @@ -37,6 +37,7 @@ arg({'patterns','PatternPairs'},uint32(3),[],'CSP patterns per frequency (times two).'), ... arg({'whycsp','SkipCSP','WhyCSP'},false,[],'Classify cross-spectrum directly. This results in much higher-dimensional features, but can be approached with appropriately regularized classifiers (see ml_trainproximal).'), ... arg({'normalize_spectrum','NormalizeSpectrum'},false,[],'Normalize the spectrum. Recommended if using sophisticated regularized classifiers.'), ... + arg({'logtransform','LogTransform'},false,[],'Log-transform output. Log-transformed spectra are more likely to be separable by a linear classifier.'), ... arg({'vectorize_features','VectorizeFeatures'},true,[],'Vectorize the features. For compatibility with basic classifiers.')); if args.signal.nbchan == 1 @@ -45,6 +46,7 @@ error('Multi-taper CSP prefers to work on at least as many channels as you request output patterns. Please reduce the number of pattern pairs.'); end if isempty(args.timewnds) args.timewnds = struct(); end + [C,dummy,T] = size(args.signal.data); %#ok if args.whycsp % shortcut for w = size(args.timewnds,1):-1:1 @@ -62,7 +64,7 @@ mean_covar{w}(~isfinite(mean_covar{w})) = 0; weighted_covar{w}(~isfinite(weighted_covar{w})) = 0; % calculate spatial filters for each frequency for f=size(mean_covar{w},1):-1:1 - [V,D] = eig(squeeze(weighted_covar{w}(f,:,:)),squeeze(mean_covar{w,1}(f,:,:))); %#ok + [V,D] = eig(reshape(weighted_covar{w}(f,:,:),C,C),reshape(mean_covar{w,1}(f,:,:),C,C)); %#ok P = inv(V); % retain k best filters/patterns at both ends of the eigenvalue spectrum filters(w,f,:,:) = real(V(:,[1:args.patterns end-args.patterns+1:end])); @@ -80,14 +82,14 @@ end % solve a CSP instance for each frequency for f=size(covar{w,1},1):-1:1 - [V,D] = eig(squeeze(covar{w,1}(f,:,:)),squeeze(covar{w,1}(f,:,:)+covar{w,2}(f,:,:))); %#ok + [V,D] = eig(reshape(covar{w,1}(f,:,:),C,C),reshape(covar{w,1}(f,:,:),C,C)+reshape(covar{w,2}(f,:,:),C,C)); %#ok P = inv(V); filters(w,f,:,:) = real(V(:,[1:args.patterns end-args.patterns+1:end])); patterns(w,f,:,:) = real(P([1:args.patterns end-args.patterns+1:end],:))'; end end end - model = struct('filters',{filters},'patterns',{patterns},'time_args',{time_args},'spec_args',{args.spectral_estimation}, 'covar',{covar}, 'mean_covar',{mean_covar}, 'weighted_covar',{weighted_covar}, 'chanlocs',{args.signal.chanlocs},'vectorize_features',{args.vectorize_features},'whycsp',{args.whycsp},'normalize_spectrum',{args.normalize_spectrum}); + model = struct('filters',{filters},'patterns',{patterns},'time_args',{time_args},'spec_args',{args.spectral_estimation}, 'covar',{covar}, 'mean_covar',{mean_covar}, 'weighted_covar',{weighted_covar}, 'chanlocs',{args.signal.chanlocs},'vectorize_features',{args.vectorize_features},'whycsp',{args.whycsp},'normalize_spectrum',{args.normalize_spectrum},'logtransform',{args.logtransform}); end global tracking; %#ok tracking.inspection.signal = args.signal; @@ -121,6 +123,8 @@ freqs = freqs(1):(freqs(2)-freqs(1))/(nfreqs-1):freqs(2); features = bsxfun(@times,features,max(1,1./freqs')); end + if featuremodel.logtransform + features = log(features); end end end % apply minimal conditioning to features diff --git a/code/utils/utl_calc_crossspec.m b/code/utils/utl_calc_crossspec.m index ab79ff80..8ea2e102 100644 --- a/code/utils/utl_calc_crossspec.m +++ b/code/utils/utl_calc_crossspec.m @@ -47,7 +47,7 @@ arg({'bandwidth','TimeBandwidth'},5,[],'Spectral smoothing. Controls the bias vs. variance of the spectral estimation. Reasonable values are 1 to 3 (1 being fairly noisy, and 3 being fairly smooth but 5x slower)'), ... arg({'tapers','Tapers'},[],[],'Number of tapers. Should be an integer smaller than 2*TimeBandwith; default 2*TimeBandwidth-1'), ... arg({'padding','Padding'},0,[],'FFT padding factor. Controls the oversampling of the spectrum; 0 is the next largest power of two, 1 is 2x as much, etc.'), ... - arg({'subsample_spectrum','SubsampleSpectrum'},1,[],'Sub-sample the spectrum. This is according to the TimeBandwidth parameter.'), ... + arg({'subsample_spectrum','SubsampleSpectrum'},1,[],'Sub-sample the spectrum. This is the sub-sampling factor.'), ... arg({'robust_estimation','RobustEstimation'},false,[],'Robust cross-spectral estimation. Whether cross-spectral matrices should be aggregated across trials in a robust manner.'), ... arg({'sum_weights','SumWeights'}, [],[],'Weights for weighted averaging. If passed, the weighted average cross-spectrum will be returned as a second output.'), ... arg({'filtered_subsampling','FilteredSubsampling'},false,[],'Use a filter prior to sub-sampling. Slower but yields a better spectral estimate.'), ... @@ -98,7 +98,7 @@ for t = 1:T cs = abs(CrossSpecMatc(signal.data(:,:,t)',signal.pnts/signal.srate,struct('tapers',[2*args.bandwidth args.tapers],'pad',args.padding,'Fs',signal.srate,'fpass',args.freqwnd))); % filter and sub-sample - if args.subsample_spectrum > 1 && args. filtered_subsampling + if args.subsample_spectrum > 1 && args.filtered_subsampling cs = filter(wnd,1,cs,1); end cs = cs(1:args.subsample_spectrum:end,:,:); if isempty(meanspec) diff --git a/code/utils/utl_partition_bundle.m b/code/utils/utl_partition_bundle.m index d7008c45..0f4f1064 100644 --- a/code/utils/utl_partition_bundle.m +++ b/code/utils/utl_partition_bundle.m @@ -19,7 +19,7 @@ if isempty(inds) % compute index set size (from first stream) - res = exp_eval(set_partition(bundle.streams{1},[],varargin{:})); + res = exp_eval_optimized(set_partition(bundle.streams{1},[],varargin{:})); else % partition the streams (symbolically) for s=1:length(bundle.streams) diff --git a/dependencies/fast_decode.m b/dependencies/fast_decode.m index 95bcec14..22470628 100644 --- a/dependencies/fast_decode.m +++ b/dependencies/fast_decode.m @@ -1,3 +1,4 @@ function bytes = fast_decode(str) +% decode a string into a uint8 vector bytes = uint8(str)-uint8(32); bytes = (bytes(1:end/2) + bitshift(bytes((end/2+1):end),uint8(4)))'; diff --git a/dependencies/fast_encode.m b/dependencies/fast_encode.m index 20c27c94..aa760656 100644 --- a/dependencies/fast_encode.m +++ b/dependencies/fast_encode.m @@ -1,3 +1,4 @@ function str = fast_encode(bytes) +% encode a uint8 vector into a string bytes = uint8(32) + [bitand(bytes,uint8(15)) bitshift(bitand(bytes,uint8(15*16)),-4)]; str = char(bytes(:)');