Skip to content

Commit

Permalink
many relevant changes
Browse files Browse the repository at this point in the history
  • Loading branch information
@andersonwinkler
andersonwinkler committed Oct 28, 2013
1 parent 7b5ced3 commit 3160a8b
Show file tree
Hide file tree
Showing 15 changed files with 839 additions and 14 deletions.
1 change: 1 addition & 0 deletions palm
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
palm.m
180 changes: 180 additions & 0 deletions palm.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,180 @@
% #!/usr/bin/octave -q
function palm(varargin)
% ===========================================
% PALM: Permutation Analysis of Linear Models
% ===========================================
%
% PALM performs permutation inference for the general linear
% models (GLMs) of arbitrary complexity, taking as inputs data
% in various formats, and being able to take into account certain
% cases of well-structured non-independence.
%
% The options are:
%
% -i <file> : Input(s). More than one can be specified.
% -m <file> : Mask(s). Either one for all inputs, or one per input,
% supplied in the same order as the respective -i appear.
% -s <file> : Surface file(s), for the inputs that are scalars over
% a surface geometry. Each -s should be supplied in the
% same order as the respective -i.
% -d <file> : Design matrix. It can be in .csv format, or in FSL's
% "VEST" format.
% -t <file> : t-contrasts file, in VEST format.
% -f <file> : F-contrasts file, in VEST format, just as in FSL.
% -eb <file> : Exchangeability blocks file, in VEST format.
% -vg <file> : Variance groups file, in VEST format.
% -o <string> : Output prefix. It may itself be prefixed by a path.
% -n <integer> : Number of permutations. Use -n 0 for exhaustive.
% -c <real> : Theshold for cluster extent for t contrasts.
% -C <real> : Theshold for cluster mass for t contrasts.
% -F <positive> : Theshold for cluster extent for F contrasts.
% -S <positive> : Theshold for cluster mass for F contrasts.
% -T : Enables TFCE inference for 3D (volume) data, i.e.,
% with H=2, E=0.5, C=6.
% -T2 : Enables TFCE inference for 2D (surface, TBSS) data,
% i.e., H=2, E=1, C=26.
% -cnpc <real> : Threshold for NPC cluster extent (z-stat).
% -Cnpc <real> : Threshold for NPC cluster mass (z-stat).
% -Tnpc : Enables TFCE inference for NPC on 3D data,
% with H=2, E=0.5, C=6.
% -T2npc : Enables TFCE inference for NPC on 2D data,
% i.e., H=2, E=1, C=26.
% -tfce_H <real> : Set the TFCE H (height power) parameter.
% -tfce_E <real> : Set the TFCE E (extent power) parameter.
% -tfce_C <integer> : Set the TFCE C (connectivity) parameter.
% -pb : Permute blocks as a whole.
% -ise : Assume independent and symmetric errors, to
% allow sign-flipping.
% -ee : Assume exchangeable errors, to allow
% permutations.
% -pmethod <string> : Method to partition the model. It can be one of:
% 'Guttman', 'Beckmann' or 'Ridgway'.
% Default is 'Beckmann'.
% -corrmod : Applies FWER-correction of p-values over multiple
% modalities.
% -corrcon : Applies FWER-correction of p-values over multiple
% contrasts.
% -saveparametric : Saves also the parametric p-values.
% -savemask : Save the effective masks used for each modality.
% -rmethod <string> : Method for regression/permutation. It can be one of:
% 'Draper-Stoneman','Still-White', 'Freedman-Lane',
% 'terBraak', 'Kennedy', 'Manly', 'Huh-Jhun' or 'Smith'.
% Default is 'Freedman-Lane'.
% -npc : Use non-parametric combination (NPC).
% -cmethod <string> : Method for combination in the NPC. It can be one
% of: 'Tippett', 'Fisher', 'Pearson-David', 'Stouffer',
% 'Wilkinson', 'Winer', 'Edgington', 'Mudholkar-George',
% 'Friston', 'Darlington-Hayes', 'Zaykin',
% 'Dudbridge-Koeleman', 'Dudbridge-Koeleman2',
% 'Nichols', 'Taylor-Tibshirani', 'Jiang'.
% Default is 'Tippett'.
% -draft : Run in the "draft mode". No FWER correction is
% possible, only FDR-adjustment.
% -fdr : Produces FDR-adjusted p-values.
%
% _____________________________________
% Anderson M. Winkler
% FMRIB / Univ. of Oxford
% Sep/2013
% http://brainder.org

% If Octave
if palm_isoctave,

% Disable memory dump on SIGTERM
sigterm_dumps_octave_core(0);

% If running as a script, take the input arguments
cmdname = program_invocation_name();
if ~strcmpi(cmdname(end-5:end),'octave'),
varargin = argv();
end
end

% This is probably redundant but fix a bug in an old version
nargin = numel(varargin);

% Print usage if no inputs are given
if nargin == 0 || strcmp(varargin{1},'-q'),
fprintf('===========================================\n');
fprintf('PALM: Permutation Analysis of Linear Models\n');
fprintf('===========================================\n');
fprintf('\n');
fprintf('PALM performs permutation inference for the general linear\n');
fprintf('models (GLMs) of arbitrary complexity, taking as inputs data\n');
fprintf('in various formats, and being able to take into account certain\n');
fprintf('cases of well-structured non-independence.\n');
fprintf('\n');
fprintf('The options are:\n');
fprintf('\n');
fprintf('-i <file> : Input(s). More than one can be specified.\n');
fprintf('-m <file> : Mask(s). Either one for all inputs, or one per input,\n');
fprintf(' supplied in the same order as the respective -i appear.\n');
fprintf('-s <file> : Surface file(s), for the inputs that are scalars over\n');
fprintf(' a surface geometry. Each -s should be supplied in the\n');
fprintf(' same order as the respective -i.\n');
fprintf('-d <file> : Design matrix. It can be in .csv format, or in FSL''s\n');
fprintf(' "VEST" format.\n');
fprintf('-t <file> : t-contrasts file, in VEST format.\n');
fprintf('-f <file> : F-contrasts file, in VEST format, just as in FSL.\n');
fprintf('-eb <file> : Exchangeability blocks file, in VEST format.\n');
fprintf('-vg <file> : Variance groups file, in VEST format.\n');
fprintf('-o <string> : Output prefix. It may itself be prefixed by a path.\n');
fprintf('-n <integer> : Number of permutations. Use -n 0 for exhaustive.\n');
fprintf('-c <real> : Theshold for cluster extent for t contrasts.\n');
fprintf('-C <real> : Theshold for cluster mass for t contrasts.\n');
fprintf('-F <positive> : Theshold for cluster extent for F contrasts.\n');
fprintf('-S <positive> : Theshold for cluster mass for F contrasts.\n');
fprintf('-T : Enables TFCE inference for 3D (volume) data, i.e.,\n');
fprintf(' with H=2, E=0.5, C=6.\n');
fprintf('-T2 : Enables TFCE inference for 2D (surface, TBSS) data,\n');
fprintf(' i.e., H=2, E=1, C=26.\n');
fprintf('-cnpc <real> : Threshold for NPC cluster extent (z-stat).\n');
fprintf('-Cnpc <real> : Threshold for NPC cluster mass (z-stat).\n');
fprintf('-Tnpc : Enables TFCE inference for NPC on 3D data,\n');
fprintf(' with H=2, E=0.5, C=6.\n');
fprintf('-T2npc : Enables TFCE inference for NPC on 2D data,\n');
fprintf(' i.e., H=2, E=1, C=26.\n');
fprintf('-tfce_H <real> : Set the TFCE H (height power) parameter.\n');
fprintf('-tfce_E <real> : Set the TFCE E (extent power) parameter.\n');
fprintf('-tfce_C <integer> : Set the TFCE C (connectivity) parameter.\n');
fprintf('-pb : Permute blocks as a whole.\n');
fprintf('-ise : Assume independent and symmetric errors, to\n');
fprintf(' allow sign-flipping.\n');
fprintf('-ee : Assume exchangeable errors, to allow\n');
fprintf(' permutations.\n');
fprintf('-pmethod <string> : Method to partition the model. It can be one of:\n');
fprintf(' ''Guttman'', ''Beckmann'' or ''Ridgway''.\n');
fprintf(' Default is ''Beckmann''.\n');
fprintf('-corrmod : Applies FWER-correction of p-values over multiple\n');
fprintf(' modalities.\n');
fprintf('-corrcon : Applies FWER-correction of p-values over multiple\n');
fprintf(' contrasts.\n');
fprintf('-saveparametric : Saves also the parametric p-values.\n');
fprintf('-savemask : Save the effective masks used for each modality.\n');
fprintf('-rmethod <string> : Method for regression/permutation. It can be one of:\n');
fprintf(' ''Draper-Stoneman'',''Still-White'', ''Freedman-Lane'',\n');
fprintf(' ''terBraak'', ''Kennedy'', ''Manly'', ''Huh-Jhun'' or ''Smith''.\n');
fprintf(' Default is ''Freedman-Lane''.\n');
fprintf('-npc : Use non-parametric combination (NPC).\n');
fprintf('-cmethod <string> : Method for combination in the NPC. It can be one\n');
fprintf(' of: ''Tippett'', ''Fisher'', ''Pearson-David'', ''Stouffer'',\n');
fprintf(' ''Wilkinson'', ''Winer'', ''Edgington'', ''Mudholkar-George'',\n');
fprintf(' ''Friston'', ''Darlington-Hayes'', ''Zaykin'',\n');
fprintf(' ''Dudbridge-Koeleman'', ''Dudbridge-Koeleman2'',\n');
fprintf(' ''Nichols'', ''Taylor-Tibshirani'', ''Jiang''.\n');
fprintf(' Default is ''Tippett''.\n');
fprintf('-draft : Run in the "draft mode". No FWER correction is\n');
fprintf(' possible, only FDR-adjustment.\n');
fprintf('-fdr : Produces FDR-adjusted p-values.\n');
fprintf('\n');
fprintf('_____________________________________\n');
fprintf('Anderson M. Winkler\n');
fprintf('FMRIB / Univ. of Oxford\n');
fprintf('Sep/2013\n');
fprintf('http://brainder.org\n');
return;
end

% Now run what matters
palm_backend(varargin);
96 changes: 96 additions & 0 deletions palm_algol.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,96 @@
function P = palm_algol(varargin)
% Given a sequence of integers S (stored as a vector of length n),
% returns a set P of permutations. P can be an array n by nP of
% permutation indices, or a cell-array with nP elements, each
% being a sparse permutation matrix.
%
% Usage:
% P = palm_algol(S,flag)
%
% S : Sequence of n integers (row or column vector).
% flag : Logical true or false. If true, P is a cell-array
% with nP elements, each a sparse permutation matrix.
% If false, P is a n by nP array with permutation
% indices. Default is false (faster and requiring
% less memory).
% P : Output permutation matrix.
%
% This function is an implementation of the "Algorithm L",
% published by D. Knuth in his "The Art of Computer Programming",
% Volume 4, Fascicle 2: Generating All Tuples and Permutations.
%
% _____________________________________
% Anderson M. Winkler
% FMRIB / University of Oxford
% Feb/2012 (first version with Algorithm L)
% Oct/2013 (this version)
% http://brainder.org

% Compute the largest number of unique permutations
S = varargin{1};
n = length(S);
U = unique(S);
nU = zeros(size(U));
for u = 1:numel(U),
nU(u) = sum(S == U(u));
end
maxP = factorial(n)/prod(factorial(nU));

% Algorithm L
% The 1st permutation is the sorted
[a,tmp] = sort(S); % sort, and keep indices to permute back
[~,idxback] = sort(tmp); % to 'unsort' the permutations
p = 1; % current permutation number
P = zeros(n,maxP); % init the array to store all perms
P(:,p) = a; % 1st permutation is the identity (Step L1)

% Step L2
j = n - 1;
while j > 0 && a(j) >= a(j+1),
j = j - 1;
end

% Interrupt if j == 0
while j > 0,

% Step L3
l = n;
while a(j) >= a(l),
l = l - 1;
end
tmp = a(j);
a(j) = a(l);
a(l) = tmp;

% Step L4
k = j + 1;
l = n;
while k < l,
tmp = a(k);
a(k) = a(l);
a(l) = tmp;
k = k + 1;
l = l - 1;
end
p = p + 1;
P(:,p) = a;

% Step L2 again
j = n - 1;
while j > 0 && a(j) >= a(j+1),
j = j - 1;
end
end

% 'Unsort' the permutations
P = P(idxback,:);

% Convert to sparse permutation matrices
if nargin > 1 && varargin{2},
[~,tmp] = sort(P);
[~,idx] = sort(tmp);
P = cell(maxP,1);
for p = 1:maxP,
P{p} = palm_idx2perm(idx(:,p));
end
end
2 changes: 1 addition & 1 deletion palm_backend.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ function palm_backend(varargin)

% Start by taking what matters from the arguments
global plm opts; % uncomment for debugging
[opts,plm] = palm_takeargs(varargin);
[opts,plm] = palm_takeargs(varargin{1});

% To store the statistic name for each contrast, to be used later when
% saving the statistic image to a file
Expand Down
5 changes: 2 additions & 3 deletions palm_mat2seq.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,6 @@
S = zeros(N,1);
U = unique(M,'rows');
for u = 1:size(U,1),
uidx = all(repmat(U(u,:),[N 1]) == M,2);
uidx = all(bsxfun(@eq,U(u,:),M),2);
S(uidx) = u;
end

end
8 changes: 4 additions & 4 deletions palm_miscread.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,7 +54,7 @@

% Read an FSL "VEST" file.
X.readwith = 'vestread';
[X.data,X.extra.PPH] = vestread(X.filename);
[X.data,X.extra.PPH] = palm_vestread(X.filename);

case '.gz',

Expand Down Expand Up @@ -131,7 +131,7 @@
'.pial','.sulc','.thickness','.volume'},

% Read a FreeSurfer curvature file
extern = checkprogs;
extern = palm_checkprogs;
if extern.fs,
X.readwith = 'fs_read_curv';
[X.data,X.extra.fnum] = read_curv(X.filename);
Expand All @@ -149,7 +149,7 @@
'smoothwm.nofix','sphere','sphere.reg','white'},

% Read a FreeSurfer surface file
extern = checkprogs;
extern = palm_checkprogs;
if extern.fs,
X.readwith = 'fs_read_surf';
[X.data.vtx,X.data.fac] = read_surf(X.filename);
Expand All @@ -165,7 +165,7 @@
case {'.mgh','.mgz'},

% Read a FreeSurfer MGH/MGZ file
extern = checkprogs;
extern = palm_checkprogs;
if extern.fs,
X.readwith = 'fs_load_mgh';
[X.data,X.extra.M,X.extra.mr_parms,X.extra.volsz] = load_mgh(X.filename);
Expand Down
Loading

0 comments on commit 3160a8b

Please sign in to comment.