function r = cons00(A, x0, eta, strategy, T, optional_param)
% r = cons00(A, x0, eta, strategy, T)
% r = cons00(A, x0, eta, strategy, T, optional_param)
%
% A:   adjacency matrix (n x n matrix)
% x0:  initial state values (n x 1 vector)
%  
%
% strategy: "*" denotes original contribution 
%   1    simple round()
%   2    probabilistic quantization
%   3  * probabilistic quantization with compensation
%   4  * probabilistic quantization with compensation, limited to 0-1
%   5  * deterministic
%   6    no quantization (y=x)
%   8    Carli et al
%   9    optional_param is bias
%  19    limited to 0-1
%
% T: number of stages to simulate
%
%
% Example use:
%
%   r = cons00(generate_connected(10), (1:10)', 0.02, 5, 1000)
	
	
	
	r = struct;
	r.strategy = strategy;
	r.n = size(A,1);
	
	if not(  (size(x0,1) == r.n) && (size(x0,2) == 1) )   
		fprintf('Size of x0 must be %d x 1, instead of %d x %d \n', r.n, size(x0,1),size(x0,2))
		error
	end
	
	if not( (eta > 0) && (eta < 1) )
		error(sprintf('eta out of bounds: %f',eta));
	end
	
	r.T = T;
	
	r.A = A;
	r.D = degree_matrix(r.A);
	r.Delta = max(max(r.D));
	r.L = r.D - r.A;

	r.eta = eta;
	r.epsilon = r.eta / r.Delta;

	r.P = eye(r.n) - r.epsilon * r.L;
	

	r.x = zeros(r.n,T) * NaN;
	r.y = zeros(r.n,T) * NaN;
	r.x(:,1) = x0;
	
	if size(x0,1) == 0
		r.x(:,1) = 10.3 + 5 * 2*(0.5-rand(r.n,1));
	end

	r.goal = mean(r.x(:,1));

	% For communication strategy 3
	r.compensate = zeros(r.n,T);
	r.error = zeros(r.n,T);

	
	strategy_carli = 8;

	for k=1:T
		progress(k,T);
	
		for i=1:r.n
	
			switch strategy
				% No discretization limits
				case 6
					r.y(i,k) = r.x(i,k);

				case {1, strategy_carli}
					r.y(i,k) = round(r.x(i,k));
					
				case 2			
					m = floor(r.x(i,k));
					p = r.x(i,k) - m;
		
					if rand < p            % with probability p
						r.y(i,k) = m + 1;
					else
						r.y(i,k) = m ;
					end

					r.error(i,k) = (r.y(i,k)-r.x(i,k));

				case 3
					value = r.x(i,k) - r.compensate(i,k);
				
					m = floor(value);
					p = value - m;

					if rand < p            % with probability p
						r.y(i,k) = m + 1;
					else
						r.y(i,k) = m ;
					end

					r.error(i,k) = (r.y(i,k)-r.x(i,k));
					r.compensate(i,k+1) = r.compensate(i,k) + r.error(i,k);

				case 4
					
					value = r.x(i,k) - r.compensate(i,k);
			
					if value > 1
						r.y(i,k) = 1;
					else
						if value < 0
							r.y(i,k) = 0;
						else
							p = value;
							if rand < p            % with probability p
								r.y(i,k) = 1;
							else
								r.y(i,k) = 0 ;
							end					
						end
					end
						
					
					r.error(i,k) = (r.y(i,k)-r.x(i,k));
					r.compensate(i,k+1) = r.compensate(i,k) + r.error(i,k);

				case 5
					value = r.x(i,k) - r.compensate(i,k);

					r.y(i,k) = round(value);

					r.error(i,k) = (r.y(i,k)-r.x(i,k));
					r.compensate(i,k+1) = r.compensate(i,k) + r.error(i,k);

				case 9
					value = r.x(i,k) - r.compensate(i,k);

					m = floor(value);
					p = optional_param;

					if rand < p            % with probability p
						r.y(i,k) = m + 1;
					else
						r.y(i,k) = m ;
					end

					r.error(i,k) = (r.y(i,k)-r.x(i,k));
					r.compensate(i,k+1) = r.compensate(i,k) + r.error(i,k);

				case 10

					value = r.x(i,k) - r.compensate(i,k);

					if value > 1
						r.y(i,k) = 1;
					else
						if value < 0
							r.y(i,k) = 0;
						else
							p = optional_param;
							if rand < p            % with probability p
								r.y(i,k) = 1;
							else
								r.y(i,k) = 0 ;
							end					
						end
					end

					r.error(i,k) = (r.y(i,k)-r.x(i,k));
					r.compensate(i,k+1) = r.compensate(i,k) + r.error(i,k);

			end
		end
	
		for i=1:r.n
			
			
			switch strategy
				case {1,2,3,4,5,6,9,10}
					s = 0;
					for j=1:r.n
						s = s + r.A(i,j) * (r.y(j,k) - r.x(i,k));
					end
		
					r.x(i,k+1) = r.x(i,k) + r.epsilon * s;
				case strategy_carli

					s = 0;
					for j=1:r.n
						s = s + r.P(i,j) * (r.y(j,k));
					end

					r.x(i,k+1) = r.x(i,k) - r.y(i,k) + s;
				
			end
		end
	end

	T = size(r.x,2);

	for k = 1:T
		x = r.x(:,k);
		r.disagreement(k) = x' * r.L * x;
		r.drift(k) = mean(r.x(:,k)) - r.goal;
		r.mse(k) = 0;
		for i=1:r.n
			r.mse(k)=r.mse(k)+ (1/r.n) * (r.x(i,k) - r.goal)^2;
		end
	end