function [items_to_take, max_fitted_value, max_fitted_volume] = backpack_binary( backpack_size, volumes, values )
	% returns list of items which fit in backpack and have maximum total value
	% backpack_size - total volume of backpack
	% volumes - vector of items volumes
	% values - vector of items values


	N=length(volumes);

	% x is input vectors with items take/ no take (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 bit
		flip_probability = 1./N; % at average 1 bit will be flipped 
		bits_to_flip = (rand(1,N) < flip_probability );
		xnew=xold;
		xnew(bits_to_flip)=xor( xnew(bits_to_flip) , 1 ); % xor operator flip the choice value
		if ( any( xnew ~= xold) )
		       % at least 1 bit flipped we are good to return
		       return;
	       else
		       % no bit flips so try again
		       xnew=change_x(xold);
	       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);
		items_max_value=sum(values); % maximum items value for very large backpack
		% for Metropolis algorithm to work E need to be positive, 
		% but smaller values are better (we are minimizing merit function)
		%E=items_max_value - items_value; 
		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); % oversize penalty
			mu=100; % penalty coefficient must be quite big
			E=E+mu*penalty;
		end
	end


	max_fitted_value=-Inf; % the maximum value which fits in the backpack
	xbest=NaN;
	xold=zeros(1,N);
	Eold=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)/50; % 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 debug 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);

		if (Enew < Eold) 
			% this point is closer to optimum
			% we go there without thinking
			xold = xnew;
			Eold = Enew;
			max_fitted_value=items_value_new;
			max_fitted_volume=items_volume_new;
		else
			% according to Metropolis algorithm
			% there is a probability exp(-E/kT) to move away from optimal
			prob=rand(1,1);
			if ( prob < exp(-(Enew-Eold)/kT) )
				xold = xnew;
				Eold = Enew;
				max_fitted_value=items_value_new;
				max_fitted_volume=items_volume_new;
			end
		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');

	% 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
			
