function maxb = palm_maxshuf(Ptree,stype,uselog)
% Computes the maximum number of possible permutations given
% a tree that specifies the depencence between the observations.
%
% Usage:
% maxb = palm_maxshuf(Ptree,ptype,uselog)
% 
% - Ptree  : Permutation tree, generated by palm_tree.
% - stype  : Shuffling type to count. It can be one of:
%            - 'perms' for permutations.
%            - 'flips' for sign-flips
%            - 'both' for permutations with sign-flips.
% - uselog : A true/false indicating whether compute in logs.
%            Default is false.
% - maxb   : Maximum number of possible shufflings.
%
% _____________________________________
% Anderson M. Winkler
% FMRIB / University of Oxford
% Oct/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/>.
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

if nargin == 1,
    stype  = 'perms';
    uselog = false;
elseif nargin == 2,
    uselog = false;
end

if uselog,
    switch lower(stype),
        case 'perms',
            maxb = lmaxpermnode(Ptree,0);
        case 'flips',
            maxb = lmaxflipnode(Ptree,0);
            maxb = maxb/log2(exp(1));
        case 'both',
            maxp = lmaxpermnode(Ptree,0);
            maxs = lmaxflipnode(Ptree,0);
            maxs = maxs/log2(exp(1));
            maxb = maxp + maxs;
    end
else
    switch lower(stype),
        case 'perms',
            maxb = maxpermnode(Ptree,1);
        case 'flips',
            maxb = maxflipnode(Ptree,1);
        case 'both',
            maxp = maxpermnode(Ptree,1);
            maxs = maxflipnode(Ptree,1);
            maxb = maxp * double(maxs);
    end
end

% ==============================================================
function np = maxpermnode(Ptree,np)
% Number of permutations per node, recursive and
% incremental.
for u = 1:size(Ptree,1),
    np = np * seq2np(Ptree{u,1}(:,1));
    if size(Ptree{u,3},2) > 1,
        np = maxpermnode(Ptree{u,3},np);
    end
end

% ==============================================================
function np = seq2np(S)
% Takes a sequence of integers and computes the 
% number of possible permutations.
U   = unique(S);
nU  = numel(U);
cnt = zeros(size(U));
for u = 1:nU,
    cnt(u) = sum(S == U(u));
end
np = factorial(numel(S))/prod(factorial(cnt));

% ==============================================================
function ns = maxflipnode(Ptree,ns)
% Number of sign-flips per node, recursive and
% incremental.
for u = 1:size(Ptree,1),
    if size(Ptree{u,3},2) > 1,
        ns = maxflipnode(Ptree{u,3},ns);
    end
    ns = ns * 2^length(Ptree{u,2});
end

% ==============================================================
function np = lmaxpermnode(Ptree,np)
% Number of permutations per node, recursive and
% incremental.
for u = 1:size(Ptree,1),
    np = np + lseq2np(Ptree{u,1}(:,1));
    if size(Ptree{u,3},2) > 1,
        np = lmaxpermnode(Ptree{u,3},np);
    end
end

% ==============================================================
function np = lseq2np(S)
% Takes a sequence of integers and computes the 
% number of possible permutations.
nS  = numel(S);
U   = unique(S);
nU  = numel(U);
cnt = zeros(size(U));
for u = 1:nU,
    cnt(u) = sum(S == U(u));
end
lfac = palm_factorial(nS);
np   = lfac(nS+1) - sum(lfac(cnt+1));

% ==============================================================
function ns = lmaxflipnode(Ptree,ns)
% Number of sign-flips per node, recursive and
% incremental. Note the in/output are base2 logarithm.
for u = 1:size(Ptree,1),
    if size(Ptree{u,3},2) > 1,
        ns = lmaxflipnode(Ptree{u,3},ns);
    end
    ns = ns + length(Ptree{u,2});
end