%MCMC knapsack
%finds a solution to 0/1 knapsack problem with 
%random weights and values

clf reset; clear all;

b=30; %weight bound
n=100; %no. of items
noite=100000; %no of iterations

b=3; n=20;noite=50000;

w=rand(1,n); %weights
v=rand(1,n); %values

hold on;
[v, iv]=sort(v); w=w(iv);
bar([1:n],v);
% plot([1:n],w, '-gs', 'LineStyle','none',...
%     'MarkerSize',5,...
%     'MarkerEdgeColor','red',...
%     'MarkerFaceColor',[1,0,0]);
bar([1:n],w, 0.2);
ylim([0,1.1]);

waitforbuttonpress;
hold off;


omg=zeros(1,n);  
weight=0; value=0; 
 

evolvalue=zeros(1,noite); 
evolweight=zeros(1,noite); 
 
gam=0.6;
T=10;
str = sprintf('$T=%1.4f, v=%1.4f, w=%1.4f$',T, value, weight);
h1=text(0.01,1.05,str, 'Interpreter', 'Latex');


for ite=1:noite
     
    if (mod(ite,10)==0) T=10*((10)^gam)/(((10)^gam)+(ite^gam));
    end
    
     ind=randi([1,n]);
    if (omg(ind)==1)
        prob=exp(-v(ind)/T);
        if (rand<prob)
            omg(ind)=0;
            value=value-v(ind);
            weight=weight-w(ind);
        end
    else
        newweight=weight+w(ind);
        if (newweight<=b) %feasibility
           omg(ind)=1;
           weight=newweight; 
           value=value+v(ind);
        end
    end
    
%     ind1=randi([1,n]);ind2=randi([1,n]);
%     if ((omg(ind1)==1)&&(omg(ind2)==0)&&(weight-w(ind1)+w(ind2)<=b)&&(v(ind2)>v(ind1)))
%         omg(ind1)=0; omg(ind2)=1; 
%         value=value-v(ind1)+v(ind2);
%         weight=weight-w(ind1)+w(ind2);
%     end
    if (mod(ite,100)==0) 
        fprintf('%d: %f %f %f\n',ite, T, value, weight);
    end
    evolvalue(ite)=value;
    evolweight(ite)=weight; %sum(omg);
    

    if (mod(ite,500)==0)
    hold on;
    delete(h1);
    bar([1:n],v, 'FaceColor',[0 0 0.5]);
    ylim([0,1.1]);
    bar(v.*omg, 'FaceColor',[0 1 0]);
    bar([1:n],w, 0.2,'FaceColor',[1 0 0]);
    str = sprintf('$T=%1.4f, v=%1.4f, w=%1.4f$',T, value, weight);
    h1=text(0.01,1.05,str, 'Interpreter', 'Latex');
    end
     
    %waitforbuttonpress;
    drawnow;
    
    
    hold off
    %waitforbuttonpress;
    
end

figure;
plot([1:noite],evolvalue, 'LineWidth',1, 'Color', 'red' );
hold on;
 
%plot([1:noite], evolweight, 'LineWidth',2, 'Color', 'blue' );
xlabel('time')
ylabel('total value');
hold off


