Commit 544cbb0b authored by Niels Kloosterman's avatar Niels Kloosterman
Browse files

fixed point averaging coarse graining + minor changes

parents a602d86c e216a76a
...@@ -410,6 +410,131 @@ for s = 1:numel(timescales) % loop through timescales ...@@ -410,6 +410,131 @@ for s = 1:numel(timescales) % loop through timescales
fprintf('%d ', ic) fprintf('%d ', ic)
y_inds = transpose(chunk_borders(ic,1):chunk_borders(ic,2)); y_inds = transpose(chunk_borders(ic,1):chunk_borders(ic,2));
% do point skipping for scales > 1, non-HP option
cg_data = {};
switch coarsegrainmethod
case 'filtskip'
if strcmp(filtmethod, 'hp')
nloops = 1; % keep original sampling rate
stepSize = 1;
else
nloops = sc;
stepSize = sc;
end
cg_data = cell(nloops,1); % make cell: cg_data{istart}{trials}(chan-by-time)
resamp_x = data_sel.trial;
for is = 1:nloops % loop over starting points here!
cg_data{is} = cellfun(@(resamp_x) resamp_x(:, is:(stepSize-1+1):end), resamp_x, 'UniformOutput', false ); % downsample data
end
clear resamp_x;
case 'pointavg' % original point averaging coarse graining, no loop over starting points
if sc == 1 || strcmp(filtmethod, 'hp') % no coarse graining for native sampling rate or high-pass entropy
cg_data{1} = data_sel.trial; %only keep trial data
nloops = 1; % no loop across starting points
else % coarse-grain time series at this time scale
nloops = 1; % no loop across starting points
nchan = size(data_sel.trial{1},1);
for itrial = 1:length(data_sel.trial)
num_cg_tpts = floor(length(data_sel.trial{itrial})/sc); % number of coarse-grained time points
cg_data{1}{itrial} = zeros(nchan, num_cg_tpts); % preallocate cg_data matrix
for t = 1:num_cg_tpts % calculate coarse_grained time series
cg_data{1}{itrial}(:,t) = mean( data_sel.trial{itrial}(:, (t-1)*sc + [1:sc]) ,2);
end
end
end
end
% after coarsegraining, loop mse computation across starting points
allcont = zeros(sc, nchan, m+1); % start_chan_m
for istart = 1:nloops
if max(cellfun(@(x) size(x,2), cg_data{istart})) == m % TODO check this at start
fprintf('Coarse grained trials below %d + 1 samples, skipping remaining starting points\n', m)
break
end
% concatenate trials and convert to single
y = single(cell2mat(cg_data{istart}));
% collect trial bounds and create mask with valid time points for pats
trl_bounds = cumsum(cellfun(@(x) size(x,2), cg_data{istart}))';
trl_mask = true(size(y,2),1);
if allowgpu && gpuavailable
trl_mask = gpuArray(trl_mask);
end
trl_mask([trl_bounds-1; trl_bounds]) = false;
% break if n data points < 100 (See Grandy et al., 2016)
ndatapoints = length(trl_mask); % TODO check this at start
if ndatapoints < 100
fprintf('N data points < 100, breaking\n')
break
end
% Calculate sample entropy of coarse grained time series
if strcmp(recompute_r, 'perscale_toi_sp')
r_new = r * std(y,1,2);
end
% keep the estimated r parameter
r_estimate(:, s, itoi, istart) = r_new; % dimord chan nsc ntoi nstartingpts
% chunk y to keep memory usage in check
max_numel = mem_available/4; % single = 4 bytes
chunk_size = round(max_numel/numel(y));
n_chunks = ceil(size(y,2)/chunk_size);
temp = 1;
chunk_borders = zeros(n_chunks,2);
for ic = 1:n_chunks
chunk_borders(ic,:) = [temp temp+chunk_size];
temp = temp+chunk_size-1; % chunks need to overlap to avoid missing pats at chunk borders
end
chunk_borders(end) = size(y,2);
clear temp
%fprintf('starting point %d\n', istart)
cont = zeros(m+1, n_chunks, nchan);
y_chunk1 = shiftdim(y', -1 ); % insert singleton dim
r_new2 = shiftdim(r_new, -2);
if allowgpu && gpuavailable
cont = gpuArray(cont);
y_chunk1 = gpuArray(y_chunk1);
r_new2 = gpuArray(r_new2);
end
fprintf('%d chunks: ', n_chunks)
for ic = 1:n_chunks
fprintf('%d ', ic)
y_inds = transpose(chunk_borders(ic,1):chunk_borders(ic,2));
y_chunk2 = permute(y_chunk1(1,y_inds,:), [2 1 3]); % insert singleton dim
if allowgpu && gpuavailable
y_chunk2 = gpuArray(y_chunk2);
end
ymat = bsxfun(@le, abs(bsxfun(@minus, y_chunk1, y_chunk2 )), r_new2 );
for ichan=1:nchan % loop since triu only supports 2D
ymat(:,:,ichan) = triu(ymat(:,:,ichan), chunk_borders(ic,1));
end
for k = 1:m+1
if k >= m % TODO try for m > 2
cont(k,ic,:) = sum(reshape(ymat(trl_mask(y_inds(1:end-2)), trl_mask, :), [], nchan));
end
if k < m+1
ymat = ymat & circshift(ymat, [-1 -1 0]);
end
end
clear ymat y_inds y_chunk2
end
allcont(istart, :, :) = gather(squeeze(sum(cont,2))'); % sum over chunks. dimord: start_chan_m
%fprintf('\n')
end % cg starting points
allcont = sum(allcont,1); % sum counts over starting points
y_chunk2 = permute(y_chunk1(1,y_inds,:), [2 1 3]); % insert singleton dim y_chunk2 = permute(y_chunk1(1,y_inds,:), [2 1 3]); % insert singleton dim
if allowgpu && gpuavailable if allowgpu && gpuavailable
......
...@@ -246,18 +246,25 @@ for s = 1:numel(timescales) % loop through timescales ...@@ -246,18 +246,25 @@ for s = 1:numel(timescales) % loop through timescales
break % subsequent time points will also not work break % subsequent time points will also not work
end end
% do point skipping for scales > 1, non-HP option
cg_data = {}; cg_data = {};
switch coarsegrainmethod switch coarsegrainmethod
case 'filtskip' case 'filtskip'
nloops = sc; if strcmp(filtmethod, 'hp')
nloops = 1; % keep original sampling rate for hp option
stepSize = 1;
else
nloops = sc;
stepSize = sc;
end
cg_data = cell(nloops,1); % make cell: cg_data{istart}{trials}(chan-by-time) cg_data = cell(nloops,1); % make cell: cg_data{istart}{trials}(chan-by-time)
for is = 1:nloops % loop over starting points here! for is = 1:nloops % loop over starting points here!
resamp_x = data_sel.trial; resamp_x = data_sel.trial;
cg_data{is} = cellfun(@(resamp_x) resamp_x(:, is:(sc-1+1):end), resamp_x, 'UniformOutput', false ); % add padding% Filter cg_data{is} = cellfun(@(resamp_x) resamp_x(:, is:(stepSize-1+1):end), resamp_x, 'UniformOutput', false ); % add padding% Filter
end end
clear resamp_x; clear resamp_x;
case 'pointavg' % original point averaging coarse graining, no loop over starting points case 'pointavg' % original point averaging coarse graining, no loop over starting points
if sc == 1 % no coarse graining for native sampling rate if sc == 1 || strcmp(filtmethod, 'hp') % no coarse graining for native sampling rate or high-pass entropy
cg_data{1} = data_sel.trial; %only keep trial data cg_data{1} = data_sel.trial; %only keep trial data
nloops = 1; % no loop across starting points nloops = 1; % no loop across starting points
else % coarse-grain time series at this time scale else % coarse-grain time series at this time scale
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment