function [DCM,BMR,BMA] = spm_dcm_bmr_all(DCM,field,OPT)
% Bayesian model reduction of all permutations of model parameters
% FORMAT [DCM,BMR,BMA] = spm_dcm_bmr_all(DCM,field,OPT)
%
% DCM      - A single estimated DCM (or PEB) structure:
%
%  DCM.M.pE  - prior expectation
%  DCM.M.pC  - prior covariance
%  DCM.Ep    - posterior expectation
%  DCM.Cp    - posterior covariances
%  DCM.beta  - prior expectation of reduced parameters (default: 0)
%              NB: beta  = 'full' uses full prior expectations
%  DCM.gamma - prior variance of reduced parameters (default: 0)
%              NB:  multiplies the prior variance from the full model
%              NB2: gamma = 'full' uses full prior variances
%
% field      - parameter fields in DCM{i}.Ep to optimise [default: {'A','B'}]
%             'All' will invoke all fields (i.e. random effects)
%             If Ep is not a structure, all parameters will be considered
%
% OPT        - Bayesian model selection or averaging: 'BMS' or 'BMA'
%              [default: 'BMA']
%
% Returns:
%
% DCM - Bayesian Model Average (BMA) over models in the final iteration of 
%       the search:
%       DCM.M.pE  - reduced prior expectation
%       DCM.M.pC  - reduced prior covariance
%       DCM.Ep    - reduced (BMA/BMS) posterior expectation
%       DCM.Cp    - reduced (BMA/BMS) posterior covariance
%       DCM.Pp    - Model posterior over parameters (with and without)
%
% BMR -  (Nsub) summary structure reporting the model space from the last
%        iteration of the search:
%
%        BMR.name - character/cell array of parameter names
%        BMR.F    - free energies (relative to full model)
%        BMR.P    - and posterior (model) probabilities
%        BMR.K    - [models x parameters] model space (1 = off, 0 = on)
%
% BMA - Baysian model average (over reduced models; see spm_dcm_bma)
%
%--------------------------------------------------------------------------
% This routine searches over reduced (nested) models of a full model (DCM) 
% using Bayesian model reduction and performs Bayesian Model Averaging.
% 'Reduced' means some free parameters (parameters with a non-
% zero prior covariance) are switched off by fixing their prior variance 
% to zero. 
%
% If there are fewer than nmax = 8 free parameters, all permutations of 
% switching off parameters will be tested. Otherwise, this routine 
% implements the following greedy search procedure. The nmax parameters 
% are identified which, when switched off individually, produce the least 
% reduction (greatest increase) in model evidence. All permutations of 
% switching off these parameters are then evaluated and the best 
% permutation is retained. This procedure is repeated until all nmax
% parameters are retained or there are no more parameters to consider. 
% Finally, BMA is performed on the models from the last iteration.
% 
% NB: The full model should be estimated prior to running this function. A
% summary of the reduced model is plotted when the number of output
% arguments is greater than one.
%
% See also: spm_dcm_post_hoc - this routine is essentially a simplified
% version of spm_dcm_post_hoc
%__________________________________________________________________________

% Karl Friston, Peter Zeidman
% Copyright (C) 2010-2022 Wellcome Centre for Human Neuroimaging


%-specification of null prior covariance
%--------------------------------------------------------------------------
if isfield(DCM,'beta'),  beta  = DCM.beta;  else, beta  = 0; end
if isfield(DCM,'gamma'), gamma = DCM.gamma; else, gamma = 0; end

%-check for preference not to produce output
%--------------------------------------------------------------------------
noplot  = spm_get_defaults('cmdline');
verbose = spm_get_defaults('dcm.verbose');

%-Check fields of parameter structure (and options)
%--------------------------------------------------------------------------
if nargin < 3
    OPT   = 'BMA';
end
if nargin < 2 || isempty(field)
    field = {'A','B'};
end
if ischar(field)
    if strcmpi(field,'All') && isstruct(DCM.Ep)
        field = fieldnames(DCM.Ep);
    else
        field = {field};
    end
end

%-deal with filenames structure
%--------------------------------------------------------------------------
if ischar(DCM)
    DCM = load(DCM,'DCM');
    DCM = DCM.DCM;
end

% Get prior covariances
%--------------------------------------------------------------------------
if isstruct(DCM.M.pC), DCM.M.pC = diag(spm_vec(DCM.M.pC)); end
if spm_length(DCM.M.pE) ~= size(DCM.M.pC,1)
    DCM.M.pC = diag(spm_vec(DCM.M.pC));
end

% Get priors and posteriors
%--------------------------------------------------------------------------
qE  = DCM.Ep;
qC  = DCM.Cp;
pE  = DCM.M.pE;
pC  = DCM.M.pC;

% Remove (a priori) null space
%--------------------------------------------------------------------------
U   = spm_svd(pC);
qE  = U'*spm_vec(qE);
pE  = U'*spm_vec(pE);
qC  = U'*qC*U;
pC  = U'*pC*U;


%-Greedy search (GS) - eliminating parameters in a top down fashion
%==========================================================================

% Accumulated reduction vector (C)
%--------------------------------------------------------------------------
q   = diag(DCM.M.pC);
if sum(q < 1024)
    C   = double(q > mean(q(q < 1024))/1024);
else
    C   = double(q > 0);
end
GS  = 1;
while GS
    
    %-Find free coupling parameters
    %----------------------------------------------------------------------
    if isstruct(DCM.Ep)
        k = spm_fieldindices(DCM.Ep,field{:});
    else
        k = 1:spm_length(DCM.Ep);
    end
    k = k(find(C(k))); %#ok<FNDSB>
    
    % If there are too many parameters find those with the least evidence
    %----------------------------------------------------------------------
    nparam = length(k);
    nmax   = fix(max(nparam/4,8));
    if nparam > nmax
        
        % Model search over new prior without the i-th parameter
        %------------------------------------------------------------------
        Z     = zeros(1,nparam);
        for i = 1:nparam
            
            % Identify parameters to retain r and to remove s
            %--------------------------------------------------------------
            r   = C; r(k(i)) = 0; s = 1 - r;

            % Create reduced prior covariance matrix
            %--------------------------------------------------------------
            if isnumeric(gamma)
                R   = U'*diag(r + s*gamma)*U;
                rC  = R*pC*R;
            else
                rC  = pC;
            end
            
            % Create reduced prior means
            %--------------------------------------------------------------
            if isnumeric(beta)
                S  = U'*diag(r)*U;
                rE = S*pE + U'*s*beta;
            else
                rE = pE;
            end
            
            Z(i) = spm_log_evidence(qE,qC,pE,pC,rE,rC);
        end
        
        % Find parameters with the least evidence
        %------------------------------------------------------------------
        [z,i] = sort(-Z);
        k     = k(i(1:nmax));
        
        % Flag a greedy search
        %------------------------------------------------------------------
        GS = 1;
        
    elseif isempty(k)
        fprintf('\nThere are no free parameters in this model.\n')
        return
    else
        GS = 0;
    end
    
    
    % compare models
    %======================================================================
    for j = 1:2
        
        if j == 1
            % compare models with and without nmax parameters first
            %--------------------------------------------------------------
            K = repmat(logical([1;0]),1,numel(k));
        else
            % compare all combinations
            %--------------------------------------------------------------
            k = k(1:min(8,end));
            K = spm_perm_mtx(numel(k));
        end
        
        % Model search over new prior (covariance)
        %------------------------------------------------------------------
        nK    = size(K,1);
        G     = zeros(1,nK);
        for i = 1:nK
            
            % Identify parameters to retain (r) and to remove (s)
            %--------------------------------------------------------------
            r    = C; r(k(K(i,:))) = 0; s = 1 - r;
            
            % Create reduced prior covariance matrix
            %--------------------------------------------------------------
            if isnumeric(gamma)
                R   = U'*diag(r + s*gamma)*U;
                rC  = R*pC*R;
            else
                rC  = pC;
            end
            
            % Create reduced prior means
            %--------------------------------------------------------------
            if isnumeric(beta)
                S  = U'*diag(r)*U;
                rE = S*pE + U'*s*beta;
            else
                rE = pE;
            end
            
            G(i) = spm_log_evidence(qE,qC,pE,pC,rE,rC);
        end
        
        % if sufficient complexity reduction then omit combinations
        %------------------------------------------------------------------
        if G(1) - G(end) > nmax && nparam > nmax
            break;
        else
            nmax = 8;
        end
    end
    
    % posterior probability
    %----------------------------------------------------------------------
    p            = spm_softmax(G(:));
    
    %-Get selected model and prune redundant parameters
    %======================================================================
    [z,i]        = max(p);
    C(k(K(i,:))) = 0;
    
    % Continue greedy search if any parameters have been eliminated
    %----------------------------------------------------------------------
    nelim  = full(sum(K(i,:)));
    GS     = GS & nelim;
    
    % Show results
    % --------------------------------------------------------------------- 
    if verbose
        fprintf('%i out of %i free parameters removed \n',nelim,nparam)
    end
    
    if ~noplot && nmax <= 8
            spm_figure('Getwin','BMR - all'); clf
            subplot(3,2,1)
            if numel(G) > 32, plot(G,'k'), else, bar(G,'c'), end
            title('log-posterior','FontSize',16)
            xlabel('model','FontSize',12)
            ylabel('log-probability','FontSize',12)
            axis square

            subplot(3,2,2)
            if numel(G) > 32, plot(p,'k'), else, bar(p,'r'), end
            title('model posterior','FontSize',16)
            xlabel('model','FontSize',12)
            ylabel('probability','FontSize',12)
            axis square
            drawnow
     end
    
end


%-Inference over families (one family per coupling parameter)
%==========================================================================
Pk = zeros(2,length(k));
for i = 1:length(k)
    Pk(1,i) = mean(p(~K(:,i))); % mean prob. of models with parameter i on
    Pk(2,i) = mean(p( K(:,i))); % mean prob. of models with parameter i off
end

% Normalise probabilities
Pk = Pk(1,:) ./ sum(Pk);

% Set within overall parameter vector
Pp    = C;
Pp(k) = Pk;

%-Bayesian model selection or average
%==========================================================================
qE    = DCM.Ep;
qC    = DCM.Cp;
pE    = DCM.M.pE;
pC    = DCM.M.pC;
pE    = spm_vec(pE);

switch OPT
    
    case('BMA')
        % Bayesian model averaging
        %------------------------------------------------------------------
        Gmax     = max(G);
        
    case('BMS')
        % Bayesian model selection (place winning G outside Occam's window
        %------------------------------------------------------------------
        [Gmax,i] = max(G);
        G(i)     = G(i) + 16;
        Gmax     = Gmax + 16;
        
    otherwise
end

BMA   = {};
for i = 1:length(K)
    
    % if this mdel is in Occam's window, include in BMA
    %----------------------------------------------------------------------
    if G(i) > (Gmax - 8)
        
        % reduced model
        %------------------------------------------------------------------
        r            = C;
        r(k(K(i,:))) = 0;
        s            = 1 - r;
        if isnumeric(gamma)
            R  = diag(r + s*gamma);
            rC = R*pC*R;
        else
            rC = pC;
        end
        S      = diag(r);
        if isnumeric(beta)
            rE = S*spm_vec(pE) + s*beta;
        else
            rE = pE;
        end
        
        % BMR
        %------------------------------------------------------------------
        [F,Ep,Cp]    = spm_log_evidence_reduce(qE,qC,pE,pC,rE,rC);
        BMA{end + 1} = struct('Ep',Ep,'Cp',Cp,'F',F);
    end
end

switch OPT
    case('BMA')
        
        % Bayesian model averaging
        %------------------------------------------------------------------
        BMA   = spm_dcm_bma(BMA);
        Ep    = BMA.Ep;
        Cp    = BMA.Cp;
        
    case('BMS')
        
        % Bayesian model selection
        %------------------------------------------------------------------
        BMA   = BMA{1};
        Ep    = BMA.Ep;
        Cp    = BMA.Cp;
        
    otherwise
end

% convert BMA variances into a covariance matrix, if necessary
%--------------------------------------------------------------------------
if isstruct(Cp) || (spm_length(Cp) == spm_length(Ep))
    Cp = diag(spm_vec(Cp));
end

% Show full and reduced conditional estimates (for Bayesian average)
%--------------------------------------------------------------------------
spm_figure('Getwin','BMR - all');

if isstruct(DCM.Ep)
    i  = spm_find_pC(pC,DCM.Ep,field);
else
    i  = 1:spm_length(DCM.Ep);
end
qE     = spm_vec(qE);
Ep     = spm_vec(Ep);

j = i(ismember(i,1:length(spm_vec(Ep))));

% BMR summary and plotting, if requested
%--------------------------------------------------------------------------
if nargout > 1
    try
        Pnames     = spm_fieldindices(DCM.Ep,k);
    catch
        try
            Np     = numel(DCM.Pnames);
            Pnames = DCM.Pnames(rem(k - 1,Np) + 1);
        catch
            Pnames = 'all parameters';
        end
    end
    BMR.name = Pnames;
    BMR.F    = G;
    BMR.P    = p;
    BMR.K    = K;
    BMR.k    = k;
    
    if ~noplot
        subplot(3,2,3), spm_plot_ci(qE(i),qC(i,i))
        title('MAP (full)','FontSize',16)
        axis square, a = axis;

        subplot(3,2,4), spm_plot_ci(Ep(j),abs(Cp(j,j)))
        title('MAP (reduced)','FontSize',16), axis square, axis(a)

        subplot(3,2,5), imagesc(1 - K')
        xlabel('model'), ylabel('parameter'), title('model space','FontSize',16)
        set(gca,'YTickLabel',BMR.name);
        axis tight, axis square

        subplot(3,2,6)
        bar(Pp(i))
        xlabel('parameter'), title(' posterior','FontSize',16)
        axis square
        drawnow, axis([0 (numel(i) + 1) 0 1])
    end
    
end


%-Save Bayesian parameter average (Ep,Cp) and family-wise inference (Pp)
%==========================================================================
if isstruct(DCM.Ep)
    if numel(Pp) < 32
        legend(spm_fieldindices(DCM.Ep,i))
    end
    Pp    = spm_unvec(Pp,DCM.Ep);
    Ep    = spm_unvec(Ep,DCM.Ep);
end

DCM.M.pE  = rE;        % reduced prior expectation
DCM.M.pC  = rC;        % reduced prior covariance
DCM.Ep    = Ep;        % Bayesian model averages
DCM.Cp    = Cp;        % Bayesian model variance
DCM.Pp    = Pp;        % Model posterior over parameters (with and without)

% Clear free energy if supplied (which is no longer meaningful)
%--------------------------------------------------------------------------
if isfield(DCM,'F')
    DCM = rmfield(DCM,'F');
end
