function mrf=gmrf_doMMD(mrf)

         cmap = load('MRF_colormap.mat'); % the colormap
            h = mrf.imagesize(1);         % height of the image
            w = mrf.imagesize(2);         % width of the image
         cnum = mrf.classnum;             % number of classes
         beta = mrf.Beta;                 % value of parameter beta
    DeltaUmin = mrf.DeltaUmin;            % value of minimal necessary energy change
            T = mrf.T0;                   % temperature at the begining
            c = mrf.c;                    % the c constant parameter
     

           cycle = 0;
    summa_deltaE = 2 * DeltaUmin; % the first iteration is guaranteed

    while summa_deltaE > DeltaUmin 
        
        % ====================================== %
        %                                        %
        %    Please, put your implementation     %
        %             BELOW THIS LINE            %
        %                                        %
        % ====================================== %
        
        summa_deltaE=0;
        cycle = cycle + 1;
        h = mrf.imagesize(1);         % height of the image
        w = mrf.imagesize(2);         % width of the image

        for y=1:h
              for x=1:w
                curr_class_label = mrf.classmask(y, x);
                mecsi_mecsi = 0;
                nem_mecsi = 0;
                for y_ = y-1 : y+1
                    for x_ = x-1 : x+1
                        if(y_ >0 && y_ < h && x_ >0 && x_ < w)
                            if(curr_class_label == mrf.classmask(y_, x_))
                                mecsi_mecsi = mecsi_mecsi + 1;
                            else
                                nem_mecsi = nem_mecsi + 1;
                            end
                        end
                    end
                end
                post_prob = mrf.logProbs{curr_class_label}(y, x);
                pri_prob = (-1 * beta * mecsi_mecsi) + (beta* nem_mecsi);
                rand_class = ceil(rand * cnum);
                while rand_class == curr_class_label
                    rand_class = ceil(rand * cnum);
                end
                mecsi_mecsi = 0;
                nem_mecsi = 0;
                for y_ = y-1 : y+1
                    for x_ = x-1 : x+1
                        if(y_ >0 && y_ < h && x_ >0 && x_ < w)
                            if(rand_class == mrf.classmask(y_, x_))
                                mecsi_mecsi = mecsi_mecsi + 1;
                            else
                                nem_mecsi = nem_mecsi + 1;
                            end
                        end
                    end
                end
                new_post_prob = mrf.logProbs{rand_class}(y, x);
                new_pri_prob = (-1 * beta * mecsi_mecsi) + (beta* nem_mecsi);
            
                U = post_prob + pri_prob;
                new_U = new_post_prob + new_pri_prob;
            
                dU = new_U - U;
                if (rand < exp(-dU/T) && dU >0 || dU < 0 ) 
                    summa_deltaE = summa_deltaE + abs(dU);
                    mrf.classmask(y, x) = rand_class;
                end
              end
        end
        T = T*c;
        % ====================================== %
        %                                        %
        %    Please, put your implementation     %
        %             ABOVE THIS LINE            %
        %                                        %
        % ====================================== %    
        
        imshow(uint8(255*reshape(cmap.color(mrf.classmask,:), h, w, 3)));
        %fprintf('Iteration #%i\n', cycle);
        title(['Class map in cycle ', num2str(cycle)]);
        drawnow;
    end
end
