% STD_READPAC - read phase-amplitude correlation
%
% Usage:
% >> [STUDY, clustinfo] = std_readpac(STUDY, ALLEEG);
% >> [STUDY, clustinfo] = std_readpac(STUDY, ALLEEG, ...
% 'key', 'val');
% Inputs:
% STUDY - studyset structure containing some or all files in ALLEEG
% ALLEEG - vector of loaded EEG datasets
%
% Optional inputs:
% 'channels' - [cell] list of channels to import {default: all}
% 'clusters' - [integer] list of clusters to import {[]|default: all but
% the parent cluster (1) and any 'NotClust' clusters}
% 'freqrange' - [min max] frequency range {default: whole measure range}
% 'timerange' - [min max] time range {default: whole measure epoch}
%
% Output:
% STUDY - (possibly) updated STUDY structure
% clustinfo - structure of specified cluster information.
%
% Author: Arnaud Delorme, CERCO, 2009-
% Copyright (C) Arnaud Delorme, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.
function [STUDY, clustinfo] = std_readpac(STUDY, ALLEEG, varargin);
if nargin < 2
help std_readpac;
return;
end
[opt moreopts] = finputcheck( varargin, { ...
'condition' 'cell' [] {};
'channels1' 'cell' [] {};
'clusters1' 'integer' [] [];
'channels2' 'cell' [] {};
'clusters2' 'integer' [] [];
'onepersubj' 'string' { 'on','off' } 'off';
'forceread' 'string' { 'on','off' } 'off';
'recompute' 'string' { 'on','off' } 'off';
'freqrange' 'real' [] [];
'timerange' 'real' [] [] }, ...
'std_readpac', 'ignore');
if ischar(opt), error(opt); end
%STUDY = pop_pacparams(STUDY, 'default');
%if isempty(opt.timerange), opt.timerange = STUDY.etc.pacparams.timerange; end
%if isempty(opt.freqrange), opt.freqrange = STUDY.etc.pacparams.freqrange; end
nc = max(length(STUDY.condition),1);
ng = max(length(STUDY.group),1);
% find channel indices
% --------------------
if ~isempty(opt.channels1)
len1 = length(opt.channels1);
len2 = length(opt.channels2);
opt.indices1 = std_chaninds(STUDY, opt.channels1);
opt.indices2 = std_chaninds(STUDY, opt.channels2);
else
len1 = length(opt.clusters1);
len2 = length(opt.clusters2);
opt.indices1 = opt.clusters1;
opt.indices2 = opt.clusters2;
end
STUDY = std_convertoldsetformat(STUDY); %XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX REMOVE WHEN READY TO GET RID OF OLD FORMAT
for ind1 = 1:len1 % usually only one channel/component
for ind2 = 1:len2 % usually one channel/component
% find indices
% ------------
if ~isempty(opt.channels1)
tmpstruct1 = STUDY.changrp(opt.indices1(ind1));
tmpstruct2 = STUDY.changrp(opt.indices2(ind2));
else
tmpstruct1 = STUDY.cluster(opt.indices1(ind1));
tmpstruct2 = STUDY.cluster(opt.indices2(ind2));
end
allinds1 = tmpstruct1.allinds;
setinds1 = tmpstruct1.setinds;
allinds2 = tmpstruct2.allinds;
setinds2 = tmpstruct2.setinds;
% check if data is already here
% -----------------------------
dataread = 0;
if isfield(tmpstruct1, 'pacdata') && strcmpi(opt.forceread, 'off') && strcmpi(opt.recompute, 'off')
if ~isempty(tmpstruct1.pacdata) && iscell(tmpstruct1.pacdata) && length(tmpstruct1.pacdata) >= opt.indices2(ind2)
if ~isempty(tmpstruct1.pacdata{opt.indices2(ind2)})
%if isequal( STUDY.etc.pacparams.timerange, opt.timerange) && ...
% isequal( STUDY.etc.pacparams.freqrange, opt.freqrange) && ~isempty(tmpstruct.pacdata)
dataread = 1;
end
end
end
if ~dataread
% reserve arrays
% --------------
% pacarray = cell( max(length(STUDY.condition),1), max(length(STUDY.group),1) );
% tmpind1 = 1; while(isempty(setinds{tmpind1})), tmpind1 = tmpind1+1; end
% tmpind2 = 1; while(isempty(setinds{tmpind2})), tmpind2 = tmpind2+1; end
% if ~isempty(opt.channels1)
% [ tmp allfreqs alltimes ] = std_readpac( ALLEEG, 'channels1' , setinds1{tmpind}(1), 'channels2' , setinds2{tmpind}(1), 'timerange', opt.timerange, 'freqrange', opt.freqrange);
% else [ tmp allfreqs alltimes ] = std_readpac( ALLEEG, 'components1', setinds1{tmpind}(1), 'components2', setinds2{tmpind}(1), 'timerange', opt.timerange, 'freqrange', opt.freqrange);
% end
% for c = 1:nc
% for g = 1:ng
% pacarray{c, g} = repmat(zero, [length(alltimes), length(allfreqs), length(allinds1{c,g}) ]);
% end
% end
% read the data and select channels
% ---------------------------------
fprintf('Reading all PAC data...\n');
for c = 1:nc
for g = 1:ng
% scan all subjects
count = 1;
for subj = 1:length(STUDY.subject)
% get dataset indices for this subject
[inds1 inds2] = getsubjcomps(STUDY, subj, setinds1{c,g}, setinds2{c,g});
if setinds1{c,g}(inds1) ~= setinds2{c,g}(inds2), error('Wrong subject index'); end
if ~strcmpi(ALLEEG(setinds1{c,g}(inds1)).subject, STUDY.subject(subj)), error('Wrong subject index'); end
if ~isempty(inds1) && ~isempty(inds2)
if ~isempty(opt.channels1)
[pacarraytmp allfreqs alltimes] = std_pac( ALLEEG(setinds1{c,g}(subj)), 'channels1' , allinds1{c,g}(inds1), 'channels2', allinds2{c,g}(inds2), 'timerange', opt.timerange, 'freqrange', opt.freqrange, 'recompute', opt.recompute, moreopts{:});
else [pacarraytmp allfreqs alltimes] = std_pac( ALLEEG(setinds1{c,g}(subj)), 'components1', allinds1{c,g}(inds1), 'components2', allinds2{c,g}(inds2), 'timerange', opt.timerange, 'freqrange', opt.freqrange, 'recompute', opt.recompute, moreopts{:});
end
% collapse first 2 dimensions (comps x comps)
if ndims(pacarraytmp) == 4
pacarraytmp = reshape(pacarraytmp, size(pacarraytmp,1)*size(pacarraytmp,2), size(pacarraytmp,3), size(pacarraytmp,4));
else pacarraytmp = reshape(pacarraytmp, 1, size(pacarraytmp,1),size(pacarraytmp,2));
end
if strcmpi(opt.onepersubj, 'on')
pacarray{c, g}(:,:,count) = squeeze(mean(pacarraytmp,1));
count = count+1;
else
for tmpi = 1:size(pacarraytmp,1)
pacarray{c, g}(:,:,count) = pacarraytmp(tmpi,:,:);
count = count+1;
end
end
end
end
end
end
% copy data to structure
% ----------------------
if ~isempty(opt.channels1)
STUDY.changrp(opt.indices1(ind1)).pacfreqs = allfreqs;
STUDY.changrp(opt.indices1(ind1)).pactimes = alltimes;
STUDY.changrp(opt.indices1(ind1)).pacdata{opt.indices2(ind2)} = pacarray;
else STUDY.cluster(opt.indices1(ind1)).pacfreqs = allfreqs;
STUDY.cluster(opt.indices1(ind1)).pactimes = alltimes;
STUDY.cluster(opt.indices1(ind1)).pacdata{opt.indices2(ind2)} = pacarray;
end
end
end
end
% return structure
% ----------------
if ~isempty(opt.channels1)
clustinfo = STUDY.changrp(opt.indices1);
else clustinfo = STUDY.cluster(opt.indices1);
end
% get components common to a given subject
% ----------------------------------------
function [inds1 inds2] = getsubjcomps(STUDY, subj, setlist1, setlist2, complist1, complist2)
inds1 = [];
inds2 = [];
datasets = strmatch(STUDY.subject{subj}, { STUDY.datasetinfo.subject } ); % all datasets of subject
[tmp1] = intersect_bc(setlist1, datasets);
[tmp2] = intersect_bc(setlist2, datasets);
if length(tmp1) > 1, error('This function does not support sessions for subjects'); end
if length(tmp2) > 1, error('This function does not support sessions for subjects'); end
if tmp1 ~= tmp2, error('Different datasets while it should be the same'); end
if ~isempty(tmp1), inds1 = find(setlist1 == tmp1); end
if ~isempty(tmp2), inds2 = find(setlist2 == tmp2); end