function [x,dist,A]=anneal(A,b,beta0,faktor,steg,iter);
% [x,dist]=anneal(A,b,beta0,faktor,steg,iter);
% Simlated annealing : Binary equation system Ax=b
%
% Input:
% A=A-matrix
% b=b-vector 
%
% NOTE: Take for example b=rand(n,1)<0.5;
% and A=rand(n,m)<delta; where delta 0.05-0.25 (larger values are more difficult)
% m can be e g 10-50. True solutions for m<=15 can be investigated using procedure solcheck.
%
% SUGGESTED VALUES : n=100 and m=50 and delta=0.1
%
% beta0=starting value for beta (1/temperature) (SUGGESTED VALUE 0.01)
% faktor= increase in beta-value in each step (SUGGESTED VALUE 1.2)
% steg= number of increases in beta-value (Number of temperatures)
%          SUGGESTED VALUE 40
% iter=numer of simulations for each beta-value (SUGGESTED VALUE 100-500) 
%
% Output:
% x=final value of x-vector
% dist=vector of successive values of |Ax-b|
%
% Note: The beta-step and last value of dist is reported 
% for each beta-value. A plot is also produced. 
% TURN OFF PLOTS BEFORE STARTING PROCEDURE!!!


% start at a random x-vector
[n m]=size(A);
xstart=rand(m,1)<0.5;
xold=xstart;
Ax=A*xstart;
% transform Ax using modulo 2-calculation
Ax=(2*floor(Ax/2)~=Ax);
% distance from target-vector b
diffold=sum(abs(Ax-b));
dist=diffold;
beta=beta0;
% start iterations for beta
for j=1:steg
% increase beta
beta=faktor*beta;
j
% start iteration for current beta-value
for i=1:iter
	xprop=xold;
% generate a uniformly distributed neighbour to current x
	prop=unidrnd(m);
% flip 0 to 1
	if xold(prop)==0, 
		xprop(prop)=1;
% calculate new Ax-value
		Axny=Ax+A(:,prop);
% transform Ax using modulo 2-calculation
		Axny=(2*floor(Axny/2)~=Axny);		
	else, 
% flip 1 to 0
		xprop(prop)=0;
% calculate new Ax-value
		Axny=Ax+2-A(:,prop);
% transform Ax using modulo 2-calculation
		Axny=(2*floor(Axny/2)~=Axny);		
	end;
% calculate distance for proposed vector
	diff=sum(abs(Axny-b));
% go to proposed vector if it is better than old vector
	if diff-diffold<=0
		Ax=Axny;
		x=xprop;
		diffold=diff;
% calculate propability if new vector is worse than old vector
	elseif rand(1,1)<exp(-beta*(diff-diffold)) 
% "worse" proposal is accepted 
			Ax=Axny;
			x=xprop;
			diffold=diff;
% "worse" proposal is not accepted
	else x=xold;
	end
	xold=x;
% update distance vector
	dist=[dist diffold];
end
% display current distance
	diffold
% display current beta-value
	beta
% display current number of non-zero entries in x-vector
	xsum=sum(x)
% make new plot plot
	plot(dist)
	pause(2)
end;
% display final distance
q=size(dist);
dist(q(2))
