Skip to content

Commit

Permalink
a94: gammainc patch, tfce1d for csv
Browse files Browse the repository at this point in the history
  • Loading branch information
@andersonwinkler
andersonwinkler committed Apr 29, 2016
1 parent af179de commit 1423abb
Show file tree
Hide file tree
Showing 8 changed files with 116 additions and 66 deletions.
3 changes: 1 addition & 2 deletions palm_core.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3735,8 +3735,7 @@ function savedof(df1,df2,fname)
% Chi^2 staistic, p-value, and z-score:
X2 = (abs(Oconf-Econf)-.5).^2./Econf;
X2 = sum(sum(X2,1),3);
df = ones(size(X2))*size(X,2);
P = gammainc(X2/2,df/2,'upper'); % division of P by 2 omitted.
P = palm_gammainc(X2/2,size(X,2)/2,'upper'); % division of P by 2 omitted.
Z = sqrt(2)*erfcinv(P); % multiplication of P by 2 omitted.

% ==============================================================
Expand Down
1 change: 1 addition & 0 deletions palm_defaults.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,6 +74,7 @@
opts.tfce.E = 0.5; % TFCE E parameter
opts.tfce.conn = 6; % TFCE connectivity neighbourhood
opts.tfce.deltah = 'auto'; % Delta-h of the TFCE equation
opts.tfce.tfce1d = false; % Has the option -tfce1d been given?
opts.useniiclass = true; % Use the NIFTI class (saves memory)
opts.inormal = false; % Do an inverse-normal transformation?
opts.inormal_meth = 'Waerden'; % Method for the inverse-normal transformation.
Expand Down
2 changes: 1 addition & 1 deletion palm_gamma.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,7 @@
tail = 'lower';
end
end
pvals = gammainc((G-cpar)./tpar,kpar,tail);
pvals = palm_gammainc((G-cpar)./tpar,kpar,tail);

% Deal with imaginary parts.
if ~ isreal(pvals),
Expand Down
47 changes: 47 additions & 0 deletions palm_gammainc.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,47 @@
function Y = palm_gammainc(X,A,tail)
% Quick patch until the "gammainc" function in Octave is fixed
% http://savannah.gnu.org/bugs/?47800
%
% Inputs and outputs are as in "gammainc", except
% that A must be scalar.
%
% _____________________________________
% Anderson Winkler
% FMRIB / University of Oxford
% Apr/2016
% http://brainder.org

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% PALM -- Permutation Analysis of Linear Models
% Copyright (C) 2016 Anderson M. Winkler
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

% For most cases, this will be enough:
Y = gammainc(X,A,tail);

% The patch is really only for Octave:
if palm_isoctave,
idx = Y == 0;
if any(idx),
gfun = @(t)exp(-t).*t.^(A-1);
if strcmpi(tail,'upper'),
Y(idx) = arrayfun(@(x)quad(gfun,x,Inf),X(idx));
else
Y(idx) = arrayfun(@(x)quad(gfun,0,x),X(idx));
end
Y(idx) = bsxfun(@rdivide,Y(idx),gamma(A));
end
end
4 changes: 2 additions & 2 deletions palm_gcdf.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,7 @@
elseif df1 < 0,

% Chi^2, via lower Gamma incomplete for precision and speed
df2 = bsxfun(@times,ones(size(G)),df2);
gcdf = gammainc(G/2,df2/2);
%df2 = bsxfun(@times,ones(size(G)),df2);
gcdf = palm_gammainc(G/2,df2/2,'lower');

end
6 changes: 2 additions & 4 deletions palm_gpval.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,9 +78,7 @@
elseif df1 < 0,

% Chi^2, via upper Gamma incomplete for precision and speed
df2 = bsxfun(@times,ones(size(G)),df2);
pvals = gammainc(G/2,df2/2,'upper');
%df2 = bsxfun(@times,ones(size(G)),df2);
pvals = palm_gammainc(G/2,df2/2,'upper');

end


117 changes: 61 additions & 56 deletions palm_takeargs.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -404,6 +404,7 @@
opts.tfce.H = 2;
opts.tfce.E = 2;
opts.tfce.conn = 6;
opts.tfce.tfce1d = true;
a = a + 1;

case '-tfce2D', % basic
Expand Down Expand Up @@ -1596,66 +1597,70 @@
% statistics has been invoked, check if surfaces are available
% and with compatible size, then compute the area (dpv or dpf).
% Also take this opportunity to compute the adjacency matrix.
if opts.spatial.do && ...
any(strcmpi(plm.masks{i}.readwith,{'load','fs_load_mgh'})),
if Ns == 0,
error([ ...
'To use spatial statistics with vertexwise or facewise data it is\n'...
'necessary to provide the surface files (with the option "-s").%s'],'');
elseif Ns == 1,
s = 1;
else
s = i;
end

% String defining the types, for the filenames and other tasks.
if any(size(plm.srf{s}.data.vtx,1) == ...
size(plm.masks{i}.data));
plm.Yissrf(i) = true;
plm.Yisvtx(i) = true;
plm.Yisfac(i) = false;
plm.Ykindstr{i} = '_dpv';
elseif any(size(plm.srf{s}.data.fac,1) == ...
size(plm.masks{i}.data));
plm.Yissrf(i) = true;
plm.Yisvtx(i) = false;
plm.Yisfac(i) = true;
plm.Ykindstr{i} = '_dpf';
else
error([...
'Surface file does not match the input data:\n' ...
'- Surface file %d has %d vertices and %d faces (%s)\n' ...
'- Input data file %d has %d points (%s)'],...
s,size(plm.srf{s}.data.vtx,1),size(plm.srf{s}.data.fac,1),opts.s{s},...
i,max(size(plm.masks{i}.data)),opts.i{i});
end

% Surface area, to be used by the spatial statistics
if isempty(plm.srfarea{s}.data),

% No area file given, use the actual surface area
plm.Yarea{i} = palm_calcarea(plm.srf{s}.data,plm.Yisvtx(i));

elseif numel(plm.srfarea{s}.data) == 1,

% A weight given (such as 1): use that for each vertex or face,
% treating as if all had the same area.
if plm.Yisvtx(i),
plm.Yarea{i} = plm.srfarea{s}.data .* ...
ones(size(plm.srf{s}.data.vtx,1),1);
elseif plm.Yisfac(i),
plm.Yarea{i} = plm.srfarea{s}.data .* ...
ones(size(plm.srf{s}.data.fac,1),1);
if opts.tfce.tfce1d,
plm.Yisvol(i) = true;
else
if opts.spatial.do && ...
any(strcmpi(plm.masks{i}.readwith,{'load','fs_load_mgh'})),
if Ns == 0,
error([ ...
'To use spatial statistics with vertexwise or facewise data it is\n'...
'necessary to provide the surface files (with the option "-s").%s'],'');
elseif Ns == 1,
s = 1;
else
s = i;
end

else
% String defining the types, for the filenames and other tasks.
if any(size(plm.srf{s}.data.vtx,1) == ...
size(plm.masks{i}.data));
plm.Yissrf(i) = true;
plm.Yisvtx(i) = true;
plm.Yisfac(i) = false;
plm.Ykindstr{i} = '_dpv';
elseif any(size(plm.srf{s}.data.fac,1) == ...
size(plm.masks{i}.data));
plm.Yissrf(i) = true;
plm.Yisvtx(i) = false;
plm.Yisfac(i) = true;
plm.Ykindstr{i} = '_dpf';
else
error([...
'Surface file does not match the input data:\n' ...
'- Surface file %d has %d vertices and %d faces (%s)\n' ...
'- Input data file %d has %d points (%s)'],...
s,size(plm.srf{s}.data.vtx,1),size(plm.srf{s}.data.fac,1),opts.s{s},...
i,max(size(plm.masks{i}.data)),opts.i{i});
end

% Surface area, to be used by the spatial statistics
if isempty(plm.srfarea{s}.data),

% No area file given, use the actual surface area
plm.Yarea{i} = palm_calcarea(plm.srf{s}.data,plm.Yisvtx(i));

elseif numel(plm.srfarea{s}.data) == 1,

% A weight given (such as 1): use that for each vertex or face,
% treating as if all had the same area.
if plm.Yisvtx(i),
plm.Yarea{i} = plm.srfarea{s}.data .* ...
ones(size(plm.srf{s}.data.vtx,1),1);
elseif plm.Yisfac(i),
plm.Yarea{i} = plm.srfarea{s}.data .* ...
ones(size(plm.srf{s}.data.fac,1),1);
end

else

% Otherwise, just use the data from the file (already loaded).
plm.Yarea{i} = plm.srfarea{s}.data;
end

% Otherwise, just use the data from the file (already loaded).
plm.Yarea{i} = plm.srfarea{s}.data;
% Compute the adjacency matrix
plm.Yadjacency{i} = palm_adjacency(plm.srf{s}.data.fac,plm.Yisvtx(i));
end

% Compute the adjacency matrix
plm.Yadjacency{i} = palm_adjacency(plm.srf{s}.data.fac,plm.Yisvtx(i));
end
end
plm.nmasks = numel(plm.masks);
Expand Down
2 changes: 1 addition & 1 deletion palm_version.txt
View file
Original file line number Diff line number Diff line change
@@ -1 +1 @@
Apr/2016 (version alpha93)
Apr/2016 (version alpha94)

0 comments on commit 1423abb

Please sign in to comment.