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);
        w = mrf.imagesize(2)
          
        for y = 1:h
              for x = 1:w
                curr_class_label = mrf.classmask(y, x);
                
                match = 0;
                no_match = 0;
                
                for i1 = y-1 : y+1
                    for i2 = x-1 : x+1
                        if(i1 >0 && i1 < h && i2 >0 && i2 < w)
                            if(curr_class_label == mrf.classmask(i1, i2))
                                match = match + 1;
                            else
                                no_match = no_match + 1;
                            end
                        end
                    end
                end
                
                post_prob = mrf.logProbs{curr_class_label}(y, x);
                pri_prob = (-1 * beta * match) + (beta* no_match);
                rand_class = ceil(rand * cnum);
                
                while rand_class == curr_class_label
                    rand_class = ceil(rand * cnum);
                end
                
                match = 0;
                no_match = 0;
                
                for i1 = y-1 : y+1
                    for i2 = x-1 : x+1
                        if(i1 >0 && i1 < h && i2 >0 && i2 < w)
                            if(rand_class == mrf.classmask(i1, i2))
                                match = match + 1;
                            else
                                no_match = no_match + 1;
                            end
                        end
                    end
                end
                
                new_post_prob = mrf.logProbs{rand_class}(y, x);
                new_pri_prob = (-1 * beta * match) + (beta* no_match);
            
                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
