mfbox_findpermutations.m

toolbox for spm 5 for data, model free analysis

function [G,costcp,costp]=mfbox_findpermutations(m,M,epsilon,maxiter,delta,verbose)
% usage: [G,costcp,costp] = mfbox_findpermutations(m,M,epsilon,maxiter,delta,verbose)
%
% Copyright by Peter Gruber and Fabian J. Theis
% Signal Processing & Information Theory group
% Institute of Biophysics, University of Regensburg, Germany
% Homepage: http://research.fabian.theis.name
%           http://www-aglang.uni-regensburg.de
%
% This file is free software, subject to the 
% GNU GENERAL PUBLIC LICENSE, see gpl.txt

if (nargin<3) epsilon = 0.1; end
if (nargin<4) maxiter = 500; end
if (nargin<5) delta = 10^(-3); end
if (nargin<6) verbose = 0; end
if (~iscell(M)), M = {M}; end

G = {};
lM = size(M);
G{1} = eye(size(M{1,1}));
for k=1:lM(1)
    for l=(k+1):lM(2)
        if (m(l,k))
            tt = G{k}*M{k,l};
            t = abs(tt);
            v = zeros(size(M{l,l}));
            for j=1:size(v,1)
                [x,y] = mmax(t);
                v(x,y) = sign(tt(x,y));
                t(x,:) = 0;
                t(:,y) = 0;
            end
            G{l} = v;
        else
            G{l} = eye(size(M{l,l}));
        end
    end
end
if (epsilon<0)
    return
end

delta = epsilon*delta;
costcp = [];
costp = [];

lM = size(M);

i = 1;
while (i<maxiter)
    c = 0;
    for k=1:lM(1)
        for l=1:lM(2)
            if ((k~=l)&&(m(k,l)))
                [dV,dW] = permgrad(G{k},M{k,l},G{l});
                oGv = G{k}; oGw = G{l};
                G{k} = G{k}-epsilon*dV;
                G{k} = (G{k}*G{k}')^(-1/2)*G{k};
                G{l} = G{l}-epsilon*dW;
                G{l} = (G{l}*G{l}')^(-1/2)*G{l};
                c = c+norm(G{k}-oGv,'fro')+norm(G{l}-oGw,'fro');
            end
        end
    end
    
    costcp(i) = c;
    if (i==1)
        maxc = c;
    end
    
    if ((verbose>0)&&(mod(i,floor(maxiter/20))==0))
        plot(gca,1:i,costcp); set(gca,'YScale','log'); drawnow;
        title(gca,'Change per iteration');
    end
    if (c<delta*maxc), break; end
    i = i+1;
end

if (i==maxiter)
    warning(sprintf('No convergence of cross permutations after %d iterations',i));
end

i = 1;
while (i<maxiter)
    c = 0;
    for k=1:lM(1)
        if (m(k,k))
            dV = offgrad(M{k,k},G{k});
            oG = G{k};
            G{k} = G{k}-epsilon*dV;
            G{k} = (G{k}*G{k}')^(-1/2)*G{k};
            c = c+norm(G{k}-oG,'fro');
        end
    end
    costp(i) = c;

    if ((verbose>0)&&(mod(i,floor(maxiter/20))==0))
        plot(gca,1:i,costp); set(gca,'YScale','log'); drawnow; 
        title(gca,'Change per iteration');
    end
    if (c<delta), break; end
    i = i+1;
end

if (i==maxiter), warning(sprintf('No convergence after %d iterations',i)); end

for k=1:lM(1)
    t = G{k};
    tt = zeros(size(t));
    [d,p] = max(abs(t));
    for i=1:size(t,1)
        tt(i,p(i)) = sign(t(i,p(i)));
    end
    G{k} = tt;
end

return


function [VC,WC]=permgrad(V,C,W)

M = diag(diag(V*C*W'));
WC = M*V*C;
WC = WC'-WC;
VC = M*W*C';
VC = VC'-VC;

return

function V=offgrad(C,V)

M = diag(diag(V*C*V'));
V = 2*M*V*C;
V = V'-V;

return

function [x,y]=mmax(M)

[d,f] = max(M);
[l,y] = max(d);
x = f(y);

return