Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Reduce volume memory usage #168

Merged
merged 18 commits into from
Mar 17, 2020
Merged

Reduce volume memory usage #168

Show file tree
Hide file tree
Changes from 12 commits
Commits
Show all changes
18 commits
Select commit Hold shift + click to select a range
348c75d
Reduce volume memory usage
glasserm Mar 3, 2020
8ad2ca5
More logical and efficient sequence of steps for the the hp>0 condition
mharms Mar 4, 2020
b7303a8
Forgot to modify the fslmaths -bptf call to operate on the original f…
mharms Mar 4, 2020
e2a0f24
Change shell rm to matlab delete
coalsont Mar 4, 2020
9ac0eaa
change single(zeros(... to zeros(..., 'single')
coalsont Mar 4, 2020
97761a9
Add a couple linebreaks just to match other code
mharms Mar 5, 2020
624c8d1
MATLAB compile: functionhighpassandvariancenormalize
hodgem Mar 9, 2020
2430741
Update icaDim.m
glasserm Mar 10, 2020
bad46a5
Update icaDim.m
glasserm Mar 10, 2020
c05996a
Update icaDim.m
glasserm Mar 10, 2020
6edc725
Update icaDim.m
glasserm Mar 11, 2020
d122f4c
Some minor code restructuring to hopefully faciliate readability
mharms Mar 11, 2020
e0c8ca2
Decrease memory usage of functionhighpass by using 'range' instead of…
mharms Mar 12, 2020
a5b0cd2
Decrease memory usage of icaDim by using 'range' instead of 'std'; Ad…
mharms Mar 12, 2020
f1b2c5f
Delete 'data' once done with it
mharms Mar 12, 2020
9f08366
Fix some whitespace (tab vs space) misalignment
mharms Mar 12, 2020
ca736b1
Add reporting of percentage of masked voxels
mharms Mar 16, 2020
a75919f
MATLAB compile: functionhighpassandvariancenormalize
hodgem Mar 16, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Binary file modified ...cripts/Compiled_functionhighpassandvariancenormalize/functionhighpassandvariancenormalize
View file
Open in desktop
Binary file not shown.
83 changes: 57 additions & 26 deletions ICAFIX/scripts/functionhighpassandvariancenormalize.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -93,13 +93,6 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin)
hpstring = ['_hp' pdstring num2str(hp)];
end

%% Load volume time series
if dovol > 0
cts=single(read_avw([fmri '.nii.gz']));
ctsX=size(cts,1); ctsY=size(cts,2); ctsZ=size(cts,3); ctsT=size(cts,4);
cts=reshape(cts,ctsX*ctsY*ctsZ,ctsT);
end

%% Load the motion confounds, and the CIFTI (if hp>=0) (don't need either if hp<0)
if hp>=0
confounds=load([fmri hpstring '.ica/mc/prefiltered_func_data_mcf.par']);
Expand All @@ -114,13 +107,28 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin)
%% Apply hp filtering of the motion confounds, volume (if requested), and CIFTI
if pdflag % polynomial detrend case
if dovol > 0
% Save and filter confounds, as a NIFTI, as expected by FIX
save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf'],'f',[1 1 1 TR]);
confounds=detrendpoly(confounds,hp);
save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf_hp'],'f',[1 1 1 TR]);


% Load volume time series and reduce to just the non-zero voxels (for memory efficiency)
ctsfull=single(read_avw([fmri '.nii.gz']));
ctsX=size(ctsfull,1); ctsY=size(ctsfull,2); ctsZ=size(ctsfull,3); ctsT=size(ctsfull,4);
ctsfull=reshape(ctsfull,ctsX*ctsY*ctsZ,ctsT);
ctsmask=std(ctsfull,[],2)>0;
cts=ctsfull(ctsmask,:);
clear ctsfull;

% Polynomial detrend
cts=detrendpoly(cts',hp)';
save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri hpstring '.nii.gz'],'f',[1 1 1 TR]);
% Use -d flag in fslcpgeom (even if not technically necessary) to reduce possibility of mistakes when editing script

% Write out result, restoring to original size
ctsfull=zeros(ctsX*ctsY*ctsZ,ctsT, 'single');
ctsfull(ctsmask,:)=cts;
save_avw(reshape(ctsfull,ctsX,ctsY,ctsZ,ctsT),[fmri hpstring '.nii.gz'],'f',[1 1 1 TR]);
clear ctsfull;
% Use -d flag in fslcpgeom (even if not technically necessary) to reduce possibility of mistakes when editing script
call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri hpstring '.nii.gz -d']);
end

Expand All @@ -129,26 +137,42 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin)

elseif hp>0 % "fslmaths -bptf" based filtering
if dovol > 0
% Save and filter confounds, as a NIFTI, as expected by FIX
save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf'],'f',[1 1 1 TR]);
call_fsl(sprintf(['fslmaths ' fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf -bptf %f -1 ' fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf_hp'],0.5*hp/TR));

save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri hpstring '.nii.gz'],'f',[1 1 1 TR]);
call_fsl(['fslmaths ' fmri hpstring '.nii.gz -bptf ' num2str(0.5*hp/TR) ' -1 ' fmri hpstring '.nii.gz']);
cts=single(read_avw([fmri hpstring '.nii.gz']));
cts=reshape(cts,ctsX*ctsY*ctsZ,ctsT);
call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri hpstring '.nii.gz -d']);
% bptf filtering; no masking here, so this is probably memory inefficient
call_fsl(['fslmaths ' fmri '.nii.gz -bptf ' num2str(0.5*hp/TR) ' -1 ' fmri hpstring '.nii.gz']);
call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri hpstring '.nii.gz -d']);

% Load in result and reduce to just the non-zero voxels (for memory efficiency going forward)
ctsfull=single(read_avw([fmri hpstring '.nii.gz']));
ctsX=size(ctsfull,1); ctsY=size(ctsfull,2); ctsZ=size(ctsfull,3); ctsT=size(ctsfull,4);
ctsfull=reshape(ctsfull,ctsX*ctsY*ctsZ,ctsT);
ctsmask=std(ctsfull,[],2)>0;
cts=ctsfull(ctsmask,:);
clear ctsfull;
end

BOdimX=size(BO.cdata,1); BOdimZnew=ceil(BOdimX/100); BOdimT=size(BO.cdata,2);
save_avw(reshape([BO.cdata ; zeros(100*BOdimZnew-BOdimX,BOdimT)],10,10,BOdimZnew,BOdimT),'Atlas','f',[1 1 1 TR]);
call_fsl(sprintf('fslmaths Atlas -bptf %f -1 Atlas',0.5*hp/TR));
grot=reshape(single(read_avw('Atlas')),100*BOdimZnew,BOdimT); BO.cdata=grot(1:BOdimX,:); clear grot;
call_fsl('rm Atlas.nii.gz');
ciftisave(BO,[fmri '_Atlas' regstring hpstring '.dtseries.nii'],WBC);
delete('Atlas.nii.gz');

elseif hp<0 % If no hp filtering, still need to at least demean the volumetric time series
if dovol > 0
cts=demean(cts')';

% Load volume time series and reduce to just the non-zero voxels (for memory efficiency)
ctsfull=single(read_avw([fmri '.nii.gz']));
ctsX=size(ctsfull,1); ctsY=size(ctsfull,2); ctsZ=size(ctsfull,3); ctsT=size(ctsfull,4);
ctsfull=reshape(ctsfull,ctsX*ctsY*ctsZ,ctsT);
ctsmask=std(ctsfull,[],2)>0;
cts=ctsfull(ctsmask,:);
clear ctsfull;

cts=demean(cts')';
end

end
Expand Down Expand Up @@ -176,14 +200,18 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin)
% (Only need these for the individual runs, which, per above, are expected to be passed in with hp requested).
if hp>=0
if dovol > 0
fname=[fmri hpstring '_vn.nii.gz'];
save_avw(reshape(Outcts.noise_unst_std,ctsX,ctsY,ctsZ,1),fname,'f',[1 1 1 1]);
call_fsl(['fslcpgeom ' fmri '_mean.nii.gz ' fname ' -d']);
fname=[fmri hpstring '_vn.nii.gz'];
vnfull=zeros(ctsX*ctsY*ctsZ,1, 'single');
vnfull(ctsmask)=Outcts.noise_unst_std;

save_avw(reshape(vnfull,ctsX,ctsY,ctsZ,1),fname,'f',[1 1 1 1]);
clear vnfull;
call_fsl(['fslcpgeom ' fmri '_mean.nii.gz ' fname ' -d']);
end

VN=BO;
VN.cdata=OutBO.noise_unst_std;
fname=[fmri '_Atlas' regstring hpstring '_vn.dscalar.nii'];
fname=[fmri '_Atlas' regstring hpstring '_vn.dscalar.nii'];
disp(['saving ' fname]);
ciftisavereset(VN,fname,WBC);
end
Expand All @@ -193,11 +221,14 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin)
% But CIFTI version of TCS only saved if hp>0
if dovol > 0
cts=cts./repmat(Outcts.noise_unst_std,1,ctsT);
% Use '_vnts' (volume normalized time series) as the suffix for the volumetric VN'ed TCS
fname=[fmri hpstring '_vnts.nii.gz'];
save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),fname,'f',[1 1 1 1]);
% N.B. Version of 'fslcpgeom' in FSL 6.0.0 requires a patch because it doesn't copy both the qform and sform faithfully
call_fsl(['fslcpgeom ' fmri '.nii.gz ' fname ' -d']);
% Use '_vnts' (volume normalized time series) as the suffix for the volumetric VN'ed TCS
fname=[fmri hpstring '_vnts.nii.gz'];
ctsfull=zeros(ctsX*ctsY*ctsZ,ctsT, 'single');
ctsfull(ctsmask,:)=cts;
save_avw(reshape(ctsfull,ctsX,ctsY,ctsZ,ctsT),fname,'f',[1 1 1 1]);
clear ctsfull;
% N.B. Version of 'fslcpgeom' in FSL 6.0.0 requires a patch because it doesn't copy both the qform and sform faithfully
call_fsl(['fslcpgeom ' fmri '.nii.gz ' fname ' -d']);
end
% For CIFTI, we can use the extension to distinguish between VN maps (.dscalar) and VN'ed time series (.dtseries)
if hp>=0
Expand Down
147 changes: 91 additions & 56 deletions ICAFIX/scripts/icaDim.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
function [Out] = icaDim(Origdata,DEMDT,VN,Iterate,NDist)
function [Out] = icaDim(data,DEMDT,VN,Iterate,NDist)
% Matthew F. Glasser, Chad Donahue, Steve Smith, Christian Beckmann

%%%%%%%%%%%%%%%%%%%%
Expand All @@ -19,60 +19,74 @@

%%%%%%%%%%%%%%%%%%%%

Nvox = size(data,1);
Ntp = size(data,2);

% Mask of non-zero voxels
mask = std(data,[],2)>0;
% Apply mask, reusing 'data' variable for memory efficiency
data = data(mask,:);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Input variables don't count as used memory until they are assigned to, nor does a trivial copy (no indexing expression). I prefer not overwriting variables with modified versions, and in this case and data_detrend, it shouldn't improve memory usage versus my version.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'll add the usemask stuff, but I don't see the disadvantage of overwriting variables with a modified version, esp. since it is a misnomer to name the variable data_detrend since whether it is actually detrended depends on the value of DEMDT.

Nmask = size(data,1);

%Remove Constant Timeseries and Detrend Data
% Again, reuse 'data' variable for memory efficiency
if DEMDT==1
data = Origdata(std(Origdata,[],2)>0,:);
% In this case, preserve the trend for adding back later
data_detrend = detrend(data')';
data_trend = data - data_detrend;
data_trend = data - data_detrend;
data = data_detrend;
clear data_detrend;
elseif DEMDT==0
data = Origdata(std(Origdata,[],2)>0,:);
data_detrend = demean(data')';
data_trend = single(zeros(size(data,1),size(data,2)));
% In this case, preserve the mean (over time) for adding back later
data_mean = mean(data')';
data = demean(data')';
elseif DEMDT==-1
data = Origdata(std(Origdata,[],2)>0,:);
data_detrend = data;
data_trend = single(zeros(size(data,1),size(data,2)));
% No demeaning or detrending; no additional prep necessary
end

%Precompute PCA reconstruction
%octave doesn't have "rng", but it does accept "rand('twister', 0);" for matlab compatibility

[u,EigS,v]=nets_svds(data_detrend',0); %Eigenvalues of detrended data
[u,EigS,v]=nets_svds(data',0); %Eigenvalues of data
Out.DOF=sum(diag(EigS)>(std(diag(EigS))*0.1)); %Compute degrees of freedom

u(isnan(u))=0; v(isnan(v))=0;

c=1;
stabCount = 0;
while stabCount < stabThresh %Loop until dim output is stable

while stabCount < stabThresh
%Loop until dim output is stable
% Note: within the while loop, 'data_vn' gets reused, but 'data' stays
% constant (albeit, possibly demeaned and/or detrended previously, per above)
c
%Variance normalization via PCA reconstruction: Isolate unstructured noise
if VN~=0
noise_unst = (u(:,Out.VNDIM(c):Out.DOF)*EigS(Out.VNDIM(c):Out.DOF,Out.VNDIM(c):Out.DOF)*v(:,Out.VNDIM(c):Out.DOF)')';
Out.noise_unst_std = max(std(noise_unst,[],2),0.001);
data_detrend_vn = data_detrend ./ repmat(Out.noise_unst_std,1,size(data_detrend,2));
Out.noise_unst_std = max(std(noise_unst,[],2),0.001); clear noise_unst;
data_vn = data ./ repmat(Out.noise_unst_std,1,Ntp);
elseif VN==0
data_detrend_vn = data_detrend;
Out.noise_unst_std = single(ones(size(data_detrend,1),1));
data_vn = data;
Out.noise_unst_std = ones(Nmask,1,'single');
end
if size(data_detrend_vn,1)<size(data_detrend_vn,2)
if size(data_vn,1)<size(data_vn,2)
if DEMDT~=-1
d_pcaNorm=eig(cov(data_detrend_vn'));
d_pcaNorm=eig(cov(data_vn'));
else
d_pcaNorm=eig(data_detrend_vn*data_detrend_vn');
d_pcaNorm=eig(data_vn*data_vn');
end
else
if DEMDT~=-1
d_pcaNorm=eig(cov(data_detrend_vn));
d_pcaNorm=eig(cov(data_vn));
else
d_pcaNorm=eig(data_detrend_vn'*data_detrend_vn);
d_pcaNorm=eig(data_vn'*data_vn);
end
end
Out.lambda_pcaNorm=flipud(d_pcaNorm);

%Enable fitting multiple noise distributions
DOF=Out.DOF;
%DOF=sum(Out.lambda_pcaNorm>std(Out.lambda_pcaNorm)*0.1); %Recompute DOF because of demeaning in cov
MaxX=length(data_detrend_vn);
MaxX=length(data_vn);
if DEMDT==-1
MaxX=2000000;
end
Expand All @@ -82,11 +96,11 @@
for i=1:NDist
[x(i),en] = FitWishart(lnb,S,DOF,MaxX,lambda(1:DOF));
lambda=lambda(1:DOF)-en(1:DOF);
en=[en;single(zeros(Out.DOF-length(en)+1,1))];
en=[en; zeros(Out.DOF-length(en)+1,1,'single')];
EN(:,i)=en(1:Out.DOF);
MaxX=x(i);
DOF=min(find(lambda <= 0));
lambda=[lambda(1:DOF-1);single(zeros(Out.DOF-DOF+1,1))];
lambda=[lambda(1:DOF-1); zeros(Out.DOF-DOF+1,1,'single')];
lambda=lambda(1:Out.DOF);
end

Expand All @@ -105,7 +119,7 @@
%MaxS=5;
%DOF=Out.DOF;
%MaxX=round(Out.x(c));
%MaxX=length(data_detrend_vn);
%MaxX=length(data_vn);
%lambda=Out.lambda_pcaNorm;
%[Out.x(c),Out.EN,Out.s(c)] = SmoothEst(lnb,MaxS,DOF,MaxX,lambda);

Expand Down Expand Up @@ -169,20 +183,20 @@
disp([' dimavg: ' mat2str(priordimavg)]);
end
end %End while loop for dim calcs

if Iterate < 0
Out.calcDim=round(priordimavg);
end
if Iterate > 3
Out.calcDim=round(mean(Out.VNDIM(4:end)));
end

if DEMDT~=-1
[u,EigS,v]=nets_svds(demean(data_detrend_vn)',0);
else
[u,EigS,v]=nets_svds(data_detrend_vn',0);
end
[u,EigS,v]=nets_svds(data_vn',0);

clear data_vn;

u(isnan(u))=0; v(isnan(v))=0;
Out.EigSAdj=single(zeros(length(EigS),1));
Out.EigSAdj=zeros(length(EigS),1,'single');
Out.grot_one=diag(EigS(1:length(Out.EN),1:length(Out.EN)));
Out.grot_two=Out.grot_one.^2;
Out.grot_three=(Out.grot_two./max(Out.grot_two)).*max(Out.lambda_pcaNorm);
Expand Down Expand Up @@ -219,15 +233,36 @@

Out.EigSAdj(1:length(Out.EN))=sqrt(Out.grot_six);

Out.data=single(zeros(size(Origdata,1),size(Origdata,2)));
if DEMDT~=-1
Out.data(std(Origdata,[],2)>0,:)=data_trend + (((u*diag(Out.EigSAdj)*v')'+repmat(mean(data_detrend_vn),size(data_detrend_vn,1),1)) .* repmat(Out.noise_unst_std,1,size(data_detrend_vn,2)));
else
Out.data(std(Origdata,[],2)>0,:)=data_trend + (((u*diag(Out.EigSAdj)*v')') .* repmat(Out.noise_unst_std,1,size(data_detrend_vn,2)));
% Form the masked voxel version of the output data
tmp=(u*diag(Out.EigSAdj)*v')'; clear v;
tmp = tmp .* repmat(Out.noise_unst_std,1,Ntp);
if DEMDT==1
% Add the trend back in
tmp = tmp + data_trend;
clear data_trend;
elseif DEMDT==0
% Add the mean back in
tmp = tmp + repmat(data_mean,1,Ntp);
elseif DEMDT==-1
% No demeaning or detrending performed; nothing to add back in
end
temp=single(zeros(size(Origdata,1),1)); temp(std(Origdata,[],2)>0,:)=Out.noise_unst_std; Out.noise_unst_std=max(temp,0.001); clear temp;

end %End function
% Create the fully-formed (non-masked) version of the output data
Out.data=zeros(Nvox,Ntp,'single');
Out.data(mask,:)=tmp;
clear tmp;

% Create fully-formed (non-masked) version of Out.noise_unst_std, and
% make sure value is at least 0.001Make sure Out.noise_unst_std is at least 0.001
temp=zeros(Nvox,1,'single');
temp(mask,:)=Out.noise_unst_std;
Out.noise_unst_std=max(temp,0.001);
clear temp;

end %End main function


%% ----------- Internal Helper Functions ------------ %%

function [out] = lpdist(in)
out=pdist(log(in));
Expand Down Expand Up @@ -387,27 +422,27 @@
end

function [O] = Smooth(M,S)
if S==0
if S==0
O=M;
else
sigma = S/(2*sqrt(2*log(2)));
%sz = S*6; % length of gaussFilter vector
%x = linspace(-sz / 2, sz / 2, sz);
%gaussFilter = exp(-x .^ 2 / (2 * sigma ^ 2));
%gaussFilter = gaussFilter / sum (gaussFilter); % normalize
else
sigma = S/(2*sqrt(2*log(2)));
%sz = S*6; % length of gaussFilter vector
%x = linspace(-sz / 2, sz / 2, sz);
%gaussFilter = exp(-x .^ 2 / (2 * sigma ^ 2));
%gaussFilter = gaussFilter / sum (gaussFilter); % normalize


width = round((6*sigma - 1)/2);
support = (-width:width);
gaussFilter = exp( -(support).^2 ./ (2*sigma^2) );
gaussFilter = gaussFilter/ sum(gaussFilter);
width = round((6*sigma - 1)/2);
support = (-width:width);
gaussFilter = exp( -(support).^2 ./ (2*sigma^2) );
gaussFilter = gaussFilter/ sum(gaussFilter);

O=single(zeros(size(M,1),size(M,2)));
O=zeros(size(M,1),size(M,2),'single');

for i=1:size(M,2)
op
O(:,i) = conv (M(:,i), gaussFilter, 'same');
%O(:,i) = filter (gaussFilter,1,M(:,i));
end
end
for i=1:size(M,2)
op
O(:,i) = conv (M(:,i), gaussFilter, 'same');
%O(:,i) = filter (gaussFilter,1,M(:,i));
end
end
end