function [count] = KRS(n,prob_add,prob_del,prob_swap,interval,total_count)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program solves Challenge 18.4 (b)
%   in Homework #2
%   SUNGWOO PARK
%   KRS.m OCT 2006
%   Minor changes: DPO 10/2006 (to comments and to normalization)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [count] = KRS(n,prob_add,prob_del,prob_swap,interval,total_count)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For an n by n lattice, we consider covering it with a mixture of
% monomers and dimers.  This function estimates the number of distinct 
% arrangements of k dimers using the KRS algorithm. 
% This algorithm is Monte Carlo based, and forms a sequence of
% dimer arrangements.  At each step, we either add a dimer, delete 
% a dimer, swap a dimer, or make a null move.
%
% This function KRS has six input parameters
%
%   n - size of lattice
%   prob_add - a probability to add a new dimer.
%   prob_del - a probability to delete an existent dimer.
%   prob_swap - a probability to swap a dimer with a new one.
%   interval - how many steps to change arrangement of dimers.
%              before counting a new arrangement
%   total_count 
%
% There is one output parameter:
%
% The array count stores the desired coefficients:
%    count(j) is the estimate of the number of distinct arrangements
%    of j-1 dimers in the nxn lattice, j=1,...,n(n-1)+1.
%
% Reference:
%  C. Kenyon, D. Randall, and A. Sinclair, 
%  Approximating the number of monomer-dimer coverings of a lattice
%  J. Stat. Phys. 83, 637 (1996)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (1) Initialization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Construct n by n array called lattice corresponding to the nodes
% of the monomer-dimer lattice.
%  The element is 0 if it is a monomer,
%                 1 if it is a dimer connected to its right neighbor,
%                -1 if it is a dimer connected to its left  neighbor,
%                 2 if it is a dimer connected to its upper neighbor,
%                -2 if it is a dimer connected to its lower neighbor,
% Initially, the lattice is filled with monomers.

lattice = zeros(n,n);
% E = number of edges.
E = 2* n*(n-1);
% upper bound of # of dimers
max_dimer = E/2;

% Initialize the counts to zero.

count = zeros(max_dimer+1,1);

numOfDimer = 0;


for i=1:total_count
    % We will rearrange a dimer interval times after counting
    for j=1:interval,
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % (1) Choose an adjacent pair of sites
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        if(rand(1,1)>0.5)
           % create horizontal edge
            x = ceil(n * rand(1,1));
            y = ceil((n-1)* rand(1,1));
            p1 = [x y];
            p2 = [x y+1];
        else
            % create vertical edge
            x = ceil((n-1)* rand(1,1));
            y = ceil(n * rand(1,1));
            p1 = [x y];
            p2 = [x+1 y];
        end
        
        statusP1 = lattice(p1(1),p1(2));
        statusP2 = lattice(p2(1),p2(2));

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % (2) - (a) when both sites have monomers
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        if statusP1 == 0 && statusP2 == 0
            if(rand(1) < prob_add)
                % decide to add with prob_add
                numOfDimer = numOfDimer +1;

                if p1(1) == p2(1) 
                    %if a new edge is horizontal
                    lattice(p1(1),p1(2)) = 1;
                    lattice(p2(1),p2(2)) = -1;
                else
                    %if a new edge is vertical
                    lattice(p1(1),p1(2)) = 2;
                    lattice(p2(1),p2(2)) = -2;
                end
            end
            
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % (2) - (b) when the sites are occupied by a single dimer
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        elseif ((statusP1 == 1 && statusP2 == -1 && p1(1)==p2(1)) ||(statusP1 == 2 && statusP2 == -2 && p1(2)==p2(2)) )
            if(rand(1) < prob_del)
                % decide to delete with prob_del
                numOfDimer = numOfDimer -1;
                lattice(p1(1),p1(2)) = 0;
                lattice(p2(1),p2(2)) = 0;
            end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % (2) - (c) when one site is occupied by a dimer
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        elseif statusP1 ==0 || statusP2 ==0
            if(rand(1) < prob_swap)
                % decide to swap with prob_swap

                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                % * delete the overlapped dimer
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                if statusP2==0
                    switch(statusP1)
                        case 1
                            lattice(p1(1),p1(2)+1)=0;
                        case -1
                            lattice(p1(1),p1(2)-1)=0;
                        case 2
                            lattice(p1(1)+1,p1(2))=0;
                        case -2
                            lattice(p1(1)-1,p1(2))=0;
                    end
                elseif statusP1 ==0
                    switch(statusP2)
                        case 1
                            lattice(p2(1),p2(2)+1)=0;
                        case -1
                            lattice(p2(1),p2(2)-1)=0;
                        case 2
                            lattice(p2(1)+1,p2(2))=0;
                        case -2
                            lattice(p2(1)-1,p2(2))=0;
                    end
                end
                
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                % * add a new dimer
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                if p1(1) == p2(1) 
                    %if a new edge is horizontal
                    lattice(p1(1),p1(2)) = 1;
                    lattice(p2(1),p2(2)) = -1;
                else
                    %if a new edge is vertical
                    lattice(p1(1),p1(2)) = 2;
                    lattice(p2(1),p2(2)) = -2;
                end

            end
        end
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % (3) count the current configuration
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    count(numOfDimer+1) = count(numOfDimer+1)+1;
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (4) DPO: Normalize the counts using the fact 
%          that C(1) = n(n-1)*2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if count(2)~=0
    count = n*(n-1)*2*floor( count./count(2) );
end
sum =0;

% Print the result and plot the function C(k).

for i=0:max_dimer
    fprintf('C(%d)=\t%d\n',i,count(i+1));
    plot(i,count(i+1),'r+');
    sum = sum + count(i+1);
    hold on ;
end
fprintf('Sum of C(k) = %d\n',sum);
xlabel('k');
ylabel('C(k)');
title('KRS algorithm');