mfbox_som.m

toolbox for spm 5 for data, model free analysis

function [grid,cost]=mfbox_som(grid,gdim,data,w,varargin)
% usage: [grid,cost]=mfbox_som(grid,gdim,data,w,varargin)
%
% grid - neuron positions to be learned (dim1 x dim2 x ... x dimn x n) or (gnum x n)
% gdim - topology of the grid (dim1,dim2,...) with dim1*dim2* ... = gnum 
% data - data postitions (num x n)
% w    - data values (num)
%
% Copyright by Peter Gruber
% 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

error(nargchk(3,Inf,nargin));

sd = size(data);
sg = size(grid);
lg = prod(gdim);

if ((nargin<4)||(isempty(w))), w = ones(sd(1:(end-1))); end

maxiter = 500;
nbhfunc = 'quad';
nhdeg = 1;
verbose = 0;
plotting = 0;

l = nargin-4;
for i=1:2:l
    op = varargin{i};
    pa = varargin{i+1};
    switch lower(op)
        case 'maxiter'
            maxiter = pa;
        case 'nbhfunc'
            nbhfunc = pa;
        case 'nhdeg'
            nhdeg = pa;
        case 'plotting'
           p = find(strcmp(pa,{'off','on'}));
           plotting = p(1)-1;
        case 'verbose'
           p = find(strcmp(pa,{'off','on'}));
           verbose = p(1)-1;
       otherwise
           warning(sprintf('Illegal parameter %s ',op));
    end
end

if (length(sd)~=2), error('wrong data'); end
if (sd(2)~=sg(end)), error('grid and data incompatible'); end
if (prod(sg(1:(end-1)))~=lg), error('grid incompatible with gdim'); end

g = mfbox_mkgrid(gdim);
sumg = sum(g.^2,2);

grid = reshape(grid,prod(sg(1:end-1)),sg(end));

distfunc = str2func(sprintf('%sdist',nbhfunc));

t = 1;
lr = lrfunc(maxiter,nhdeg);
cost = [];
while (t<maxiter)
    oldgrid = grid;
    bmu = mfbox_voronoimap(data,grid);
    s = zeros(lg,sg(end));
    for j=1:lg, s(j,:) = sum(data(bmu==j,:).*repmat(w(bmu==j),1,sd(2)),1); end
    nv = zeros(lg,1);
    for j=1:lg, nv(j) = sum(w(bmu==j)); end
    delta = zeros(lg,sg(end));
    for i=1:lg,
        d = distfunc(sqrt(sumg(i)+sumg-2*g*g(i,:)'),lr(t));
        dnv = sum(d.*nv);
        if (dnv>eps), grid(i,:) = sum(repmat(d,1,sg(end)).*s,1)/dnv; end
    end
    if (plotting>0), scplot(grid'); drawnow; end
    c = sum(sum((grid-oldgrid).^2,2));
    if (c<sqrt(eps)), break; end
    cost(t) = c;
    if ((verbose>0)&&(mod(t,floor(maxiter/20))==0)), plot(cost); drawnow; end
    t = t+1;
end

return


function r=lrfunc(m,d)
if (d==1), r = linspace(1,0,m);
elseif (d==Inf), r = exp(1:m)/exp(1);
else r = (1:m).^(-d);
end


function d=expdist(d,rho)
d = exp((-d.^2)/(2*rho^2));


function d=chardist(d,rho)
d = 1*(abs(d)<rho);


function d=quaddist(d,rho)
d = max(0,1-(rho-d).^2);


function d=lindist(d,rho)
d = max(0,1-abs(rho-d));