function [Model]=TrainMESOM(Samples,Type,NumRowsMap,NumColsMap,NumSteps)
% Train a M-estimator SOM model (MESOM)
% Inputs:
%   Samples=Input samples (one sample per column)
%   Type=M-estimator type (see valid types below), as discussed in:
%       Zhang, Zhengyou (1997). Parameter estimation techniques: a tutorial 
%       with application to conic fitting. Image and Vision Computing 15,
%       59-76.
%   NumRowsMap,NumColsMap=Size of the self-organizing map
%   NumSteps=Number of time steps
% Output:
%   Model=Trained MESOM model

if isempty(Samples),
    Model = [];
else
    
    [Dimension,NumSamples]=size(Samples);

    switch Type
        case 'L2'
            Psi=@(x) x; % L2 (original Kohonen's SOFM)
        case 'L1'
            Psi=@(x) sign(x); % L1
        case 'L1-L2'
            Psi=@(x) x./sqrt(1+0.5*x.^2); % L1-L2
        case 'Lp'
            Psi=@(x) sign(x).*(abs(x).^(1.2-1)); % Lp with \nu=1.2
        case 'Fair'
            Psi=@(x) x./(1+abs(x)/1.3998); % 'Fair' with c=1.3998
        case 'Huber'
            Psi=@(x) x.*(abs(x)<=1.345)+1.345*sign(x).*(abs(x)>1.345); % Huber with k=1.345
        case 'Cauchy'
            Psi=@(x) x./(1+(x/2.3849).^2); % Cauchy with c=2.3849
        case 'German-McClure'
            Psi=@(x) x./((1+x.^2).^2); % German-McClure
        case 'Welsch'
            Psi=@(x) x.*exp(-(x/2.9846).^2); % Welsch with c=2.9846
        case 'Tukey'
            Psi=@(x) (abs(x)<=4.6851).*x.*((1-(x/4.6851).^2).^2); % Tukey's biweight with c=4.6851 (MAL)
        otherwise
            error('Unknown M-estimator type')
    end


    % Initialization
    fprintf('Initializing MESOM')
    NumNeuro=NumRowsMap*NumColsMap;
    Model.NumColsMap=NumColsMap;
    Model.NumRowsMap=NumRowsMap;
    Model.Dimension=Dimension;
    Model.Prototypes=zeros(Dimension,NumRowsMap,NumColsMap);
    Model.Psi=Psi;
    Model.Type=Type;
    NumPatIni=max([Dimension+1,ceil(NumSamples/(NumRowsMap*NumColsMap))]);

    for NdxRow=1:NumRowsMap
        for NdxCol=1:NumColsMap
            MySamples=Samples(:,ceil(NumSamples*rand(1,NumPatIni)));
            Model.Prototypes(:,NdxRow,NdxCol)=mean(MySamples,2);        
            fprintf('.')
        end
    end
    fprintf('\n');

    % Compute scaling constant
    TestCentroids=Samples(:,ceil(NumSamples*rand(1,NumNeuro)));
    DistCentroids=distance(Samples,TestCentroids);
    MinimalDists=min(DistCentroids,[],2);
    Model.Sigma=median(MinimalDists);
    Model.InvSigma=1/Model.Sigma;

    % Precompute topological distances
    [AllXCoords,AllYCoords]=ind2sub([NumRowsMap NumColsMap],1:NumNeuro);
    AllCoords(1,:)=AllXCoords;
    AllCoords(2,:)=AllYCoords;
    TopolDist=cell(NumNeuro,1);
    for NdxNeuro=1:NumNeuro    
        TopolDist{NdxNeuro}=sum((repmat(AllCoords(:,NdxNeuro),1,NumNeuro)-AllCoords).^2,1);
    end

    % Training
    fprintf('Training MESOM\n')
    MaxRadius=(NumRowsMap+NumColsMap)/8;
    for NdxStep=1:NumSteps
        MySample=Samples(:,ceil(NumSamples*rand(1)));
        if NdxStep<0.5*NumSteps   
            % Ordering phase: linear decay
            LearningRate=0.4*(1-NdxStep/NumSteps);
            MyRadius=MaxRadius*(1-(NdxStep-1)/NumSteps);
        else
            % Convergence phase: keep constant
            LearningRate=0.01;
            MyRadius=0.1;
        end

        % Compute squared Euclidean distances to the input sample and winner
        % neuron
        Diffs=repmat(MySample,1,NumNeuro)-Model.Prototypes(:,:);
        SquaredDists=sum(Diffs.^2,1);
        [Minimum,NdxWinner]=min(SquaredDists);
        Dists=sqrt(SquaredDists);

        % Compute update coefficients
        Coef=Model.Sigma*repmat(LearningRate*exp(-TopolDist{NdxWinner}/(MyRadius^2)).*Psi(Model.InvSigma*Dists),Dimension,1);

        % Compute normalized difference vectors
        NormDiffs=Diffs./repmat(Dists,[Dimension 1]);

        % Update the neurons
        if ~isnan(Coef),            
            Model.Prototypes(:,:)=Model.Prototypes(:,:)+...
                Coef.*NormDiffs;
        end

    end

    fprintf('Training finished\n')
end
    
    
        
