function [items_to_take, max_fitted_value, max_fitted_volume] = backpack_metropolis( backpack_size, volumes, values )
	% Solving the backpack problem with the simulated annealing (aka Metropolis) algorithm.
	% Returns list of items which fit in backpack and have maximum total value.
	% backpack_size - the total volume of backpack
	% volumes - the vector of items volumes
	% values - the vector of items values


	N=length(volumes);

	% x is input vector with the  no take/take flag/bit (i.e, 0/1 values) for each item
	function xnew=change_x(xold)
		% new vector will be generated via random mutation 
		% of every take/no take flag
		flip_probability = 1./N; % in average 1 bit will be flipped 
		bits_to_flip = (rand(1,N) < flip_probability );
		xnew=xold;
		xnew(bits_to_flip)=xor( xold(bits_to_flip) , 1 ); % xor operator flip the chosen flag
		if ( any( xnew ~= xold) )
		       % at least 1 flag is flipped we are good to return
		       return;
	       else
		       % none of the flags is flipped, so we try again
		       xnew=change_x(xold); % recursive call to itself
	       end
	end

	function [E,  items_value, items_volume] =backpack_merit(x, backpack_size, volumes, values)
		% calculating merit function for backpack
		items_volume=sum(volumes .* x);
		items_value=sum(values  .* x);
		% for Metropolis algorithm to work E we need to look for the minimum of E
		E= - items_value; 

		if ( (items_volume > backpack_size) )
			% Items do not fit and backpack, bad choice of input vector.
			% We need to add penalty.
			penalty=(items_volume-backpack_size); % overfill penalty
			% The penalty coefficient (mu) must be quite big,
			% but not too big or we will stack in a local minimum.
			% Choosing this coefficient require a little tweaking and 
			% depends on size of backpack, values and volumes vectors
			mu=100; 
			E=E+mu*penalty;
		end
	end


	%% Initialization
	xold=zeros(1,N);
	xbest=xold;
	[Eold, max_fitted_value, max_fitted_volume]=backpack_merit(xold, backpack_size, volumes, values);

	Ncycles=1000;
	kT=max(values)*5; % should be large enough to permit even large and non optimal merit values
	kTmin=min(values)/5; % should be smaller then the smallest step in energy
	annealing_coef= power(kTmin/kT, 1/Ncycles); % kTmin=kT*annealing_coef^Ncycles 

	dbg_energy=zeros(1,Ncycles); % this array is used for debugging or illustrations

	% the main annealing cycle 
	for c=1:Ncycles
		xnew=change_x(xold);
		[Enew, items_value_new, items_volume_new] = ...
	       		backpack_merit(xnew, backpack_size, volumes, values);

		prob=rand(1,1);
		if ( (Enew < Eold) || ( prob < exp(-(Enew-Eold)/kT) ) )
			% Either this point has smaller energy
			%  and we go there without thinking
			% or 
			%  according to the Metropolis algorithm
			%  there is the probability exp(-dE/kT) to move away from the current optimum
			xold = xnew;
			Eold = Enew;
			max_fitted_value=items_value_new;
			max_fitted_volume=items_volume_new;
		end
		% anneal or cool the temperature
		kT=annealing_coef*kT;

		dbg_energy(c)=Eold;
	end
	plot(1:Ncycles, dbg_energy);  % debug/demo plot
	xlabel('cycle number');
	ylabel('merit function');

	% the Metropolis algorithm can return non valid solution
	% but no checks are done for simplicity
	indexes=1:N;
	items_to_take=indexes(xold==1); % think about this for a while
end
			
