clear
% path(path,'H:/Book/Programs')
% A matlab demonstration of the CE method for 
% solving a MAXCUT problem
% Problem is so that it is best to parition the first
% m and last m nodes.

% Written by Dirk Kroese 
% and Shie Mannor

% Setting parameters
m = 200; % (2m nodes total)
N = 1000; % sample size
rho = 0.1; 
maxIters=23; % Number of iterations

% set random seed
rand('state',1234);
% constructing the cost matrix
Z11 = rand(m);
for i = 1:m,
  for j =i:m,
 	Z11(j,i) = Z11(i,j);
  end;
  Z11(i,i) = 0;
end

Z22 = rand(m);
for i = 1:m,
  for j =i:m,
 	Z22(j,i) = Z22(i,j);
  end;
  Z22(i,i) = 0;
end

B12  = ones(m);
B21 = ones(m);

C = [Z11 B12; B21 Z22]; % cost matrix


% Calculating optimal score

%p = rand(1,2*m) % for testing
optp = [ones(1,m),zeros(1,m)]
y = rand(1,2*m);
x = (y <optp);  % generate cut vector, according to p
s = scut(C,x);
fprintf('Best score %f\n',s)

%Initializing
p = 1/2 * ones (1,2*m);
p(1) = 1;
gammas = zeros(1,maxIters);
curbests = zeros(1,maxIters);
ps = zeros(maxIters,2*m);
pdist = zeros(1,maxIters);

tic
curbest=0.0;
for j=(1:maxIters) % main CE loop
j % output interation number
% generate matrix X of cutvectors
     Y = rand(N,2*m);
     X = (Y < ones(N,1)*p);
     g = zeros(N,1);
     for i=1:N
           g(i,:) = scut(C,X(i,:));
     end
     [sortcut,sortidx] = sort(g); 
                         % sortidx(1) contains the index of the smallest
                         % cutvector
     % update gamma
     eidx = ceil((1-rho)*N);  % smallest index of best 
     gamma = sortcut(eidx);
     curbest = max(curbest,sortcut(end)) %keep track of best result
     gammas(j)=gamma;
	 curbests(j) = curbest;
     pdist(j) = min(norm(p -optp),norm(p - ~optp));

     % update p on basis of elite vectors, that is, the vectors
     % X(sortidx(eidx),:) to X(sortidx(N),:)
     for i=1:2*m
       p(i)= sum(X(sortidx(eidx:N),i))/(N-eidx+1);
     end
 	  ps(j,:) = p;   	
end
toc  %end timing

% gammas is the calculated gamma
% pdist is the distance from the optimal p
finalp = p

fprintf('Calculated best score %f\n',curbest);
figure(1);
plot(1:maxIters,gammas);
xlabel('t');ylabel('\gamma_t');
figure(4);
plot(1:maxIters,pdist);
xlabel('Number of iteration');ylabel('||p-p_o_p_t ||_2');
figure(2); 
subplot(10,1,1); bar(1/2*ones(1,2*m));
for j=1:9
   subplot(10,1,j+1); bar(ps(j,:));
end
figure(3); 
for j=1:10
   subplot(10,1,j); bar(ps(j+9,:));
end
% use plot editing to get y-axis limits in first plot from o to 1
% use the property editor> tools> alignment tool to get more space between
% the plots. Export as .eps
   
for i = 1:maxIters
  fprintf('%d   %f  %f  %f\n',i,gammas(i),curbests(i), pdist(i));
end
