function [ber]=multiuser(SNR,method,K,Nc,seed,runname)
% MULTIUSER  CDMA multiuser detection
%
%    BER=MULTIUSER(SNR,METHOD,K,NC,SEED,RUNNAME) calculates 
%    bit error rates BER for the specified signal-to-ratios 
%    SNR=-10*log10(2*sigma^2), where sigma^2 is the noise 
%    variance on the signal. SNR can be a vector of SNR 
%    values for instance SNR=0:25 as used in Refs. [1,2]. 
%    Other inputs:
%
%    METHOD can one of the following (optional, default is 
%    adaptiveTAP):
%
%	adaptiveTAP		adaptive TAP   
%	exact			exact calculation of 
%				individual optimal 	
%				detector 
%	naive			naive mean field aka 
%				substractive interference
%				cancellation with tanh 
%				non-linearity
%	selfaveragingTAP	self-averaging TAP
%
%    K is the number of users (optional, default is 8)
%    NC is the spreading factor (optional, default is 16)
%    SEED is the seed for the random generators 
%    (optional, default is 0).
%    RUNNAME is filename to which results are written
%    (optional, default is no output written to file).
%
%    The program uses a number of other paramters set in
%    the top of the program:  
%
%    Mean field annealing [1] is used with temperatures
%    from Tmax=1 down to Tmin=sigma^2 with NT=10 either 
%    logaritmically equally spaced temperatures or equally 
%    spaced temperatures, linearTscale=1. Default is set 
%    to equally spaced temperatures.   
%
%    The maximal number of iterations at each temperature
%    max_ite=100 and for adaptive TAP a minimum of 
%    min_ite=3 iterations is set. Convergence is tested in 
%    terms of the summed squared deviation between the lhs 
%    and rhs of the mean field equations for the posterior 
%    mean. The fault tolerance for this quantity is set to 
%    ftol=10^-6. 
%
%    The number of bit errors simulated is set to 
%    Nerrors=prc^-2=800 with prc=0.05/sqrt(2) or to at most
%    Nframes=1e12 have be tested.


% References [1,2]:
%
%   Analysis of Mean Field Annealing in Substactive Interference Cancellation
%   Thomas Fabricius and Ole Winther
%   Submitted to IEEE Transactions on Communications (2002).
%
%   Correcting the Bias of Subtractive Interference Cancellation in CDMA:
%   Advanced Mean Field Theory
%   Thomas Fabricius and Ole Winther
%   Submitted to IEEE Transactions on Information Theory (2003).

global max_ite;

% set default parameter
ftol=10^-6;
Nframes=1e12; 
max_ite=100;
sigma_all=(10.^(-SNR/20))/sqrt(2);
prc=.05/sqrt(2); %Accuracy of BER curve -- Nerrors=prc^-2=800 with this choice
Nerrors=prc^-2;
Tmax=1; % should in more generality be maximum user power max(A_k^2) 
NT=10;  % number of temperatures
linearTscale=1;

% set input parameters
try K; catch K=8; end
try Nc; catch Nc=16; end
try method; catch method='adaptiveTAP'; end
try seed; catch seed=0; end
exactmean = 0; 
switch method
  case 'adaptiveTAP'
    disp(' %%%%%%% Runs adaptive TAP mean field annealing %%%%%%%  '); 
    al_cmd=sprintf( '[M_est,V,dMdM,ite]=adaptiveTAP(-R/T,Z/T,M_est,ones(K,1)/T,V,ftol);');     
  case 'exact'
    disp(' %%%%%%% Runs exact mean values %%%%%%% ');  
    exactmean = 1;  
  case 'naive'
    disp(' %%%%%%% Runs naive mean field annealing %%%%%%% ');
    al_cmd=sprintf('[M_est,dMdM,ite]=naive(-R/T,Z/T,M_est,ftol);');    
  case 'selfaveragingTAP'
    disp(' %%%%%%% Runs self-averaging TAP mean field annealing %%%%%%%  ');
    al_cmd=sprintf('[M_est,dMdM,ite]=selfaveragingTAP(-R/T,Z/T,M_est,ftol,Nc,T);');    
end

berT=zeros(NT,length(SNR));
ber=zeros(1,length(SNR));

j=0;							
for sigma=sigma_all					% run over all SNRs
  tic
  j=j+1;
  currentSNR=SNR(j)
  Tmin=sigma^2; 
  if linearTscale
    Tvals=fliplr(linspace(Tmin,Tmax,NT));
  else % logaritmic T scale
    Tvals=fliplr(exp(linspace(log(Tmin),log(Tmax),NT)));
  end
  i=0; count=0; tot_ite=0; Sum_er=0; Sum_erT=zeros(NT,1); non_converged=0;
  while and(Sum_er<(prc^-2),i<Nframes)			% loop over frames
     i=i+1;

     rand('seed',seed+i);randn('seed',seed+i);          % set seed 
     S=sign(rand(K,Nc)-.5)/sqrt(Nc);			% generate codes
     R=S*S';                                           	% code correlation matrix
     [Z_single,Z,B]=generate_data(R,1,K,sigma,eye(K)); 	% generate data
     							% Z is received signal, 
							% B is transmitted bits

     R=R-eye(K);					% substract diagonal
      
     if exactmean					% run exact
       M_est=exact(Z,R,sigma); 
       this_error=sum(sum((1-sign(M_est).*B)/2));
     else 						% run mean field algorithm
       M_est=zeros(K,1); V=zeros(K,1);
       Tindex=NT; Told=Tmax;  
       for T=Tvals				      	% annealing loop 
         eval(al_cmd);
         count=count + 1; tot_ite=tot_ite + ite; 
         if ite==max_ite				% check convergence
           non_converged=non_converged+1;
           if Tindex==NT  
             disp('didn''t convergence at highest T');
             this_error=sum(sum((1-sign(M_est).*B)/2));
           end 
           break; % break if not converged!!!
         end
         this_error=sum(sum((1-sign(M_est).*B)/2));
         Sum_erT(Tindex)=this_error+Sum_erT(Tindex);
         Tindex=Tindex-1; Told=T;
       end 
       if Tindex>0 % use last converged solution as prediction for remaining Ts
         Sum_erT(Tindex:-1:1)=Sum_erT(Tindex:-1:1)+this_error;
       end 
     end   
   
     Sum_er=this_error+Sum_er;
  end

  							% output to screen:
  ber(j)=(Sum_er/K/i)                       		% ber
  time_total=toc                			% total time 
  total_number_of_frames=i             			% total number of frames 
  time_per_frame=time_total/i       			% time per frame

  if not(exactmean)
    berTe4(:,j)=Sum_erT/K/i*1e4                   	% ber times 1e4 all T
    Tvals=(fliplr(Tvals))'				% temperatures used:
    average_ite=tot_ite/count 				% average number of iterations
    non_converged                  			% number non-converged
  end

  if nargin > 5     
    cmd=sprintf('save %s',runname); 
    eval(cmd);
  end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                 generate data                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Z_single,Z,B]=generate_data(R,N,K,sigma,A)
% Z is received signal -- M_est=Z is the conventional detector
% B is transmitted bits
% R is is code correlation
% sigma noise std., so noise cov = sigma*R
B=sign(rand(K,N)-0.5);
Noise=randn(K,N);
[U,D]=eig(R);
D=diag(D);
D(find(D<0))=0;
Z=R*A*B+sigma*U*diag(D.^.5)*U'*Noise;
Z_single=A*B+sigma*Noise;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                  adaptive TAP                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [M,V,dMdM,ite]=adaptiveTAP(J,Z,M,V_0,V,ftol)
global max_ite;
min_ite=3;
K=length(M);
alpha=M; 
ite=0; dMdM=Inf; Hnegcount=0;
Omega=zeros(K,1); dV=zeros(K,1); dM=zeros(K,1); dMdg=zeros(K,1); eta=0.9; etaV=0.5; 
H=J;
while (dMdM>ftol & ite<max_ite) | ite<min_ite 
  ite=ite+1;
  indx=randperm(K);
  for sindx=1:K
    cindx=indx(sindx);

    % update hcav and V
    hpred=H(cindx,:)*alpha;
    Hc=H(cindx,cindx);
    dV(cindx)=Hc./(1+Hc.*Omega(cindx))-V(cindx);
    V(cindx)=V(cindx)+etaV*dV(cindx);
    V(cindx)=V(cindx).*(V(cindx)>0);  
    hcav=hpred-V(cindx)*(Omega(cindx)*hpred+alpha(cindx)); 
       
    % update S and derivative
    M_old=M(cindx);
    M(cindx)=eta*tanh(hcav+Z(cindx))+(1-eta)*M(cindx); dMdg(cindx)=1-M(cindx)^2;
    dM(cindx)=M(cindx)-M_old;

    % update alpha and Omega
    alpha_old=alpha(cindx); alpha(cindx)=(M(cindx)-dMdg(cindx)*hcav)./(1+dMdg(cindx)*V(cindx));
    Omega_old=Omega(cindx); dOmega=dMdg(cindx)/(1+dMdg(cindx)*V(cindx))-Omega(cindx);
    Omega(cindx)=Omega(cindx)+dOmega;
          
    % update H - make sure that diag(M)>0 to ensure non-negative V 
    if any((diag(H)+dOmega/(1-dOmega*H(cindx,cindx)-eps)*H(:,cindx).^2)<zeros(K,1)) 
      Hnegcount=Hnegcount+1;
      dOmega=0; Omega(cindx)=Omega_old+dOmega;
      alpha(cindx)=alpha_old; 
    end
    if (dOmega~=0) % use Sherman-Woodbury
      b=H(:,cindx); H=H+dOmega/(1-dOmega*H(cindx,cindx))*b*b'; 
    end

  end % over variables - sindx  

  dMdM=dM'*dM;

end; % epoch loop - ite

%dVdV=dV'*dV

% calculate chi - covariance matrix
%diagOmega=diag(Omega);
%chi = diagOmega+diagOmega*H*diagOmega
 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%              exact posterior mean                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function M=exact(Z,R,sigma)
[K,N]=size(Z);
Bs=(de2bi(0:2^K-1)'-.5)*2;
E=repmat(diag(Bs'*R*Bs),[1,N])-2*Bs'*Z;
W=exp(-(E-repmat(min(E),[2^K,1]))/(2*sigma^2));
isin=isinf(W);
if any(isin)           %if any weights are inf
   [I,J]=find(E==min(E));
   %[I,J]=find(isin);   %make them delta function
   W(:,J)=isin(:,J);   %and the rest zero
end
W=W/max(max(W));
M=real((Bs*W./repmat(sum(W,1),[K,1])));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%               naive mean field                   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [M,dMdM,ite]=naive(J,Z,M,ftol)
global max_ite;
K=length(M); 
ite=0; dMdM=Inf;
dM=zeros(K,1);
while  dMdM>ftol & ite<max_ite
  ite=ite+1;
  indx=randperm(K);
  for sindx=1:K
    cindx=indx(sindx);  
    M_old=M(cindx);
    M(cindx)=tanh(J(cindx,:)*M+Z(cindx)); 
    dM(cindx)=M(cindx)-M_old;
  end      
  dMdM=dM'*dM;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                self-averaging TAP                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [M,dMdM,ite]=selfaveragingTAP(J,Z,M,ftol,Nc,T,A)
global max_ite;
try A=A; catch A=1; end
K=length(M);
ite=0; dMdM=Inf;
dM=zeros(K,1);
q=M'*M/K;
V=0;
while  dMdM>ftol & ite<max_ite
  ite=ite+1;
  indx=randperm(K);
  for sindx=1:K
    cindx=indx(sindx);
    M_old=M(cindx);
    V=K/Nc*A^4/T^2*(1-q)/(1+K/Nc*A^2/T*(1-q)); 	% CDMA-model
    %V=J(cindx,:).^2*(1-M.^2);         		% SK-model
    M(cindx)=tanh(J(cindx,:)*M+Z(cindx)-V*M(cindx));
    dM(cindx)=M(cindx)-M_old;
    q=q+(M(cindx)^2-M_old^2)/K; 
  end
  dMdM=dM'*dM;
end