function [NewXMap,NewYMap,WeightEvolution] = LearnBDK(Xtr,Ytr,XMap,YMap,MLKP);

if (MLKP.MaxIter < 1)
    NewXMap=XMap;
    NewYMap=YMap;
    WeightEvolution=0;
    return
end

OldRadiusX=floor(realmax);
OldRadiusY=floor(realmax);
[Nobjects,NvarX] = size(Xtr);
[NvarY,Nunits]=size(YMap);
WeightValues = DetermineWeights(Ytr,MLKP);
DistWXYAlpha = DetermineAlphaValues(MLKP);
if (upper(MLKP.MLXModeObjectWise) == 'N')
    if (upper(MLKP.MLXModeWeightsReverse) == 'Y')
        DistWYXAlpha=fliplr(DistWXYAlpha);
    else
        DistWYXAlpha=DistWXYAlpha;
    end
end

Accept=0;
indx=0;
for iter=1:MLKP.MaxIter
    % create neighbourhood and weights
    [AlphaX,AlphaY,RadiusX,RadiusY,DecayFract] = AdaptLearnParams(iter,MLKP);
    if (RadiusX < OldRadiusX)
        [ActUnitsX, WeightsX] = MakeWinnerTable(RadiusX,MLKP);
        OldRadiusX=RadiusX;
    end
    if (RadiusY < OldRadiusY)
        [ActUnitsY, WeightsY] = MakeWinnerTable(RadiusY,MLKP);
        OldRadiusY=RadiusY;
    end
    if (upper(MLKP.WeightDecay) == 'Y' | upper(MLKP.AddNoise) == 'Y')
        [XMap] = WeightDecayAddNoise(XMap, DecayFract, MLKP);
        [YMap] = WeightDecayAddNoise(YMap, DecayFract, MLKP);
    end   
    if (upper(MLKP.FastUpdate) == 'Y' & RadiusX == MLKP.FastUpdateSize...
                                      & RadiusY == MLKP.FastUpdateSize)
        [NewXMap,NewYMap] = FastUpdateMapXY(Xtr,XMap,Ytr,YMap,MLKP,iter);
        return;
    end
    
    OldXMap=XMap;
    OldYMap=YMap;
    perm=randperm(Nobjects);
    % fast mode, update maps per object
    if (upper(MLKP.MLXModeObjectWise) == 'Y')
        for iobj=1:Nobjects
            X=Xtr(perm(iobj),:);
            Y=Ytr(perm(iobj),:);
            [ValueX, ValueY, WinnerX, WinnerY] = DSWFast(X,XMap,Y,YMap,MLKP,DistWXYAlpha(iter));
            if (Boltzmann((ValueX+ValueY)/2,iter,MLKP) == 'Y')
                Accept=Accept+1;
                AF = WeightValues(perm(iobj))*ValueX;
                [NewXMap] = AdaptMap(X,XMap,WinnerX,AlphaX*AF,ActUnitsX,WeightsX);
                XMap=NewXMap;
                AF = WeightValues(perm(iobj))*ValueY;
                [NewYMap] = AdaptMap(Y,YMap,WinnerY,AlphaY*AF,ActUnitsY,WeightsY);
                YMap=NewYMap;
            end
        end
    % slow mode, update maps per whole training set
    else
        if (mod(iter,2)==1)
            for iobj=1:Nobjects
                X=Xtr(perm(iobj),:);
                Y=Ytr(perm(iobj),:);
                [Value, Winner] = DetermineSimultaneousWinner(X,XMap,Y,YMap,...
                                  MLKP,DistWYXAlpha(iter));
                if (Boltzmann(Value,iter,MLKP) == 'Y')
                    Accept=Accept+1;
                    AF = WeightValues(perm(iobj))*Value;
                    [NewXMap] = AdaptMap(X,XMap,Winner,AlphaX*AF,ActUnitsX,WeightsX);
                    XMap=NewXMap;
                end
            end
        else
            for iobj=1:Nobjects
                X=Xtr(perm(iobj),:);
                Y=Ytr(perm(iobj),:);
                [Value, Winner] = DetermineSimultaneousWinner(X,XMap,Y,YMap,...
                                  MLKP,DistWXYAlpha(iter));
                if (Boltzmann(Value,iter,MLKP) == 'Y')
                    Accept=Accept+1;
                    AF = WeightValues(perm(iobj))*Value;
                    [NewYMap] = AdaptMap(Y,YMap,Winner,AlphaY*AF,ActUnitsY,WeightsY);
                    YMap=NewYMap;
                end
            end
        end
    end
    WeightEvolution(1,iter) = sum(sum(abs(XMap-OldXMap)))/AlphaX;
    WeightEvolution(2,iter) = sum(sum(abs(YMap-OldYMap)))/AlphaY;
    if (iter > 1 & CheckConvergence(OldXMap, XMap, MLKP) == 'Y')
        return;
    end
end

if (upper(MLKP.BoltzmannMode) == 'THR' | upper(MLKP.BoltzmannMode) == 'EXP')
    Message=sprintf('MLX: %d updates out of %d trials accepted by Boltzmann algorithm',...
                    Accept, MLKP.MaxIter*Nobjects);
    disp(Message);
end