palm_fliptree.m

function Pset = palm_fliptree(Ptree,nP,cmc,idxout,maxP)
% Return a set of permutations from a permutation tree.
% 
% Usage:
% Sset = palm_fliptree(Ptree,nS,cmc,idxout,maxS)
% 
% Inputs:
% - Ptree  : Tree with the dependence structure between
%            observations, as generated by 'palm_tree'.
% - nS     : Number of permutations. Use 0 for exhaustive.
% - cmc    : A boolean indicating whether conditional
%            Monte Carlo should be used or not. If not used,
%            there is a possibility of having repeated
%            permutations. The more possible permutations,
%            the less likely to find repetitions.
% - idxout : (Optional) is supplied, Pset is an array of indices
%            rather than a cell array with sparse matrices.
% - maxS   : (Optional) Maximum number of possible sign flips.
%            If not supplied, it's calculated internally. If
%            supplied, it's not calculated internally and some
%            warnings that could be printed are omitted.
%            Also, this automatically allows nS>maxS (via CMC).
%
% Outputs:
% - Sset   : A cell array of size nP by 1 containing sparse
%            sign-flipping matrices.
% 
% Reference:
% * Winkler AM, Webster MA, Vidaurre D, Nichols TE, Smith SM.
%   Multi-level block permutation. Neuroimage. 2015;123:253-68.
%
% _____________________________________
% Anderson M. Winkler
% FMRIB / University of Oxford
% Nov/2013
% http://brainder.org

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% PALM -- Permutation Analysis of Linear Models
% Copyright (C) 2015 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/>.
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

% Note that the varnames follow the pattern of the similar
% function palm_permtree.m, for easier readability.

% Get the number of possible sign-flips.
if nargin < 5,
    maxP = palm_maxshuf(Ptree,'flips');
    if nP > maxP,
        nP = maxP; % the cap is only imposed if maxP isn't supplied
    end
end
if nargin < 4,
    idxout = false;
end

% Sign-flip #1 is no flip, regardless.
P = cell2mat(pickflip(Ptree,{},ones(size(Ptree,1)))');
P = horzcat(P,zeros(length(P),nP-1));

% All other sign flips up to nP
if nP == 0 || nP == maxP,
    
    % This will compute exhaustively all possible sign flips,
    % one branch at a time. If nP is too large print a warning.
    if nP > 1e5 && nargin <= 3;
        warning([...
            'Number of possible sign-flips is %g.\n' ...
            '         Performing all exhaustively.'],maxP);
    end
    for p = 2:maxP,
        Ptree  = nextflip(Ptree);
        P(:,p) = cell2mat(pickflip(Ptree,{},ones(size(Ptree,1)))');
    end
    
elseif cmc || nP > maxP,

    % Conditional Monte Carlo. Repeated sign flips allowed.
    for p = 2:nP,
        Ptree  = randomflip(Ptree);
        P(:,p) = cell2mat(pickflip(Ptree,{},ones(size(Ptree,1)))');
    end
    
else
    
    % Otherwise, repeated sign-flips are not allowed.
    % For this to work, maxP needs to be reasonably larger than
    % nP, otherwise it will take forever to run, so print a
    % warning.
    if nP > maxP/2 && nargin <= 3,
        warning([
            'The maximum number of sign flips (%g) is not much larger than\n' ...
            'the number you chose to run (%d). This means it may take a while (from\n' ...
            'a few seconds to several minutes) to find non-repeated sign flips.\n' ...
            'Consider instead running exhaustively all possible' ...
            'flips. It may be faster.'],maxP,nP);
    end
    
    % For each flip, keeps trying to find a new, non-repeated instance
    for p = 2:nP,
        whiletest = true;
        while whiletest,
            Ptree  = randomflip(Ptree);
            P(:,p) = cell2mat(pickflip(Ptree,{},ones(size(Ptree,1)))');
            whiletest = any(all(bsxfun(@eq,P(:,p),P(:,1:p-1))));
        end
    end
end

% Sort correctly rows using the 1st permutation
[~,idx] = palm_permtree(Ptree,1,false);
P = P(idx,:);

% For compatibility, convert each permutaion to a sparse permutation
% matrix. This section may be removed in the future if the
% remaining of the code is modified.
if idxout,
    Pset = P;
else
    Pset = cell(nP,1);
    for p = 1:nP,
        Pset{p} = sparse(diag(P(:,p)));
    end
end

% ==============================================================
function [Ptree,incremented] = nextflip(Ptree)
% Make the next sign flip of tree branches, and returns
% the shuffled tree. This can be used to compute exhaustively
% all possible sign flippings.

% Some vars for later
nU = size(Ptree,1);

% For each branch of the current node
for u = 1:nU,
    if isempty(Ptree{u,2}),
        
        % If the branches at this node cannot be considered for
        % flipping, go to the deeper levels, if they exist.
        if size(Ptree{u,3},2) > 1,
            [Ptree{u,3},incremented] = nextflip(Ptree{u,3});
            if incremented,
                if u > 1,
                    Ptree(1:u-1,:) = resetflips(Ptree(1:u-1,:));
                end
                break;
            end
        end
        
    else
        % If the branches at this node are to be considered
        % for sign-flippings (already being done or not)
        
        if sum(Ptree{u,2}) < numel(Ptree{u,2}),
            % If the current branch can be flipped, but haven't
            % reached the last possibility yet, flip and break
            % the loop.
            Ptree{u,2} = palm_incrbin(Ptree{u,2});
            incremented = true;
            if u > 1,
                Ptree(1:u-1,:) = resetflips(Ptree(1:u-1,:));
            end
            break;
        else
            % If the current branch could be flipped, but
            % it's the last possibility, reset it and 
            % don't break the loop.
            incremented = false;
        end
    end
end

% ==============================================================
function Ptree = resetflips(Ptree)
% Recursively reset all flips of a permutation tree
% back to the original state

for u = 1:size(Ptree,1),
    if isempty(Ptree{u,2}) && size(Ptree{u,3},2) > 1,
        Ptree{u,3} = resetflips(Ptree{u,3});
    else
        Ptree{u,2} = false(size(Ptree{u,2}));
    end
end

% ==============================================================
function Ptree = randomflip(Ptree)
% Make the a random sign-flip of all branches in the tree.

% For each branch of the current node
nU = size(Ptree,1);
for u = 1:nU,
    if isempty(Ptree{u,2}) == 1 && size(Ptree{u,3},2) > 1,
        % Go down more levels
        Ptree{u,3} = randomflip(Ptree{u,3});
    else
        % Or make a random flip if no deeper to go
        Ptree{u,2} = rand(size(Ptree{u,2})) > .5;
    end
end

% ==============================================================
function P = pickflip(Ptree,P,sgn)
% Take a tree in a given state and return the sign flip. This
% won't flip, only return the indices for the already flipped
% tree. This function is recursive and for the 1st iteration,
% P = {}, i.e., a 0x0 cell.

nU = size(Ptree,1);
if size(Ptree,2) == 3,
    for u = 1:nU,
        if isempty(Ptree{u,2}),
            bidx = sgn(u)*ones(size(Ptree{u,3},1),1);
        else
            bidx = double(~Ptree{u,2});
            bidx(Ptree{u,2}) = -1;
        end
        P = pickflip(Ptree{u,3},P,bidx);
    end
elseif size(Ptree,2) == 1,
    for u = 1:nU,
        if numel(sgn) == 1,
            v = 1;
        else
            v = u;
        end
        P{numel(P)+1} = sgn(v)*ones(size(Ptree{u}));
    end
end