function colors = bsizer_rand_l(Y,T,nu_0,sig_0,a_l0,b_l0,meth,alpha,Sflag,l_vec,nsB,nsM)
%
% function MAP = bsizer_rand_l(Y,T,nu_0,sig_0,a_l0,b_l0,meth,alpha,Sflag,l_vec,nB,nsM)
%
% Bayesian Sizer - Gamma prior for lambda_0 (Uses Gibbs sampling)
%                 
%
% INPUTS:
% Y     - data points (n*1 -vector)
% T     - measure points, e.g. time points (n*1 -vector)
% nu_0  - scalar parameter nu_0 for the prior of sigma^2 (deg. of freedom)
% sig_0 - scalar parameter sigma_0 for the prior of sigma^2 (NOTE: not a squared sigma_0!)
% a_l0  - prior shape parameter for lambda_0
% b_l0  - prior scale parameter for lambda_0
% meth  - 0 independent maps, 
%         1 conservative simultaneous, 
%         2 simultaneous using highest pointwise probabilities
%         [OPTIONAL: default 0]
% alpha - Confidence level (0.5 < alpha < 1)
%         [OPTIONAL: default 0.8]
% Sflag - 0 - binary coloring 
%         1 - color maps with shades (works only with independent and MP maps)
%         [OPTIONAL: default 0]          
% l_vec - column vector containing the smoothing parameters in increasing order
%         [OPTIONAL: default logspace(log10(.5),6,80)]
% nsB   - Length of burnin in MCMC simulation
%         [OPTIONAL: default 1000]
% nsM   - Length of monitor phase in MCMC simulation
%         [OPTIONAL: default 3000]
%
% OUTPUT: Bayesian Sizer map 
%
% The routine uses Statistics toolbox (to generate random numbers)
%


[rY,cY]=size(Y);
[rT,cT]=size(T);
T=T';

error(nargchk(6,12,nargin))

if cY>1 | cT>1
  error('Y and T must be column vectors')
end;

if rY~=rT
  error('Y and T must be of same length')
end;

if nu_0<=2 
   disp('Note: Parameter nu_0 must be larger than 2 to get finite mean')
end;

if sig_0<=0 
  error('Parameter sig_0 (sigma_0) must be positive')
end;

if a_l0<=0 
  error('Shape parameter a_l0 must be positive')
end;

if b_l0<=0 
  error('Scale parameter b_l0 must be positive')
end;

if nargin <= 6 ;    %  
  meth_u = 0 ;
else ;              %  meth was specified, use that
  if length(meth)>0
     meth_u = meth ;
  else
    meth_u = 0 ;
  end;  
end ;
meth=meth_u;
clear meth_u;

if ~(ismember(meth,[0 1 2]))
 error('meth must be either 0, 1 or 2') 
end;

if nargin <= 7;    
  alpha_u = 0.80 ;
else ;              %  alpha was specified, use that
  if length(alpha)>0
    alpha_u = alpha ;
  else
    alpha_u = 0.8;
  end; 
end ;
alpha=alpha_u;
clear alpha_u;

if ~(alpha>0.5 & alpha<1)
  error('alpha must be between 0.5 and 1')
end;

if nargin <= 8 ;    
  Sflag_u = 0;    
else ;              %  Sflag was specified, use that
  Sflag_u = Sflag ;
end ;
Sflag=Sflag_u;
clear Sflag_u;

if ~(ismember(Sflag,[0 1]))
 error('Sflag must be either 0 or 1') 
end;

if nargin <= 9 ;    
  l_vec_u = [0 logspace(log10(.5),6,80)];    
else               %  l_vec was specified, use that
  if length(l_vec)>0  % The input was not empty matrix
     l_vec_u = l_vec ;
  else
     l_vec_u = [0 logspace(log10(.5),6,80)];
  end;  
end ;
l_vec=l_vec_u;
clear l_vec_u;

[sv,si]=sort(l_vec);
if min(l_vec)<0 | sum(abs(si-[1:length(l_vec)]))>1
  error('elements of l_vec should be non negative and in increasing order')
end;

if nargin <= 10 ;      
  nsB_u = 1000 ;  %1000
else ;              %  nsB was specified, use that
  nsB_u = nsB ;
end ;
nsB=nsB_u;
clear nsB_u;

if nsB<0 | (nsB-ceil(nsB)>0)
  error('length of burnin phase must be non negative integer')
end;  

if nargin <= 11 ;      
  nsM_u = 3000 ;  
else ;              %  nsM was specified, use that
  nsM_u = nsM ;
end ;
nsM=nsM_u;
clear nsM_u;

if nsM<=100 | (nsM-ceil(nsM)>0)
  error('length of monitor phase must be positive integer and larger than 100')
end;  




% Lets us first define some variables that are needed later
Gamma=diag(diff(T));
n_lvalues=length(l_vec);
n=length(Y);
nu=nu_0+n-2;
D=diag(-1*ones(n,1))+diag(ones(n-1,1),1);
D=D(1:n-1,:);
ce=5; % Exponent used in color mapping (in case of maps with color shades)

% Background coloring of the maps
colors(:,:,1)=.29*ones(n_lvalues,(n-1)); % Initialize red channel   %.5
colors(:,:,2)=.42*ones(n_lvalues,(n-1)); % Initialize green channel %.2 
colors(:,:,3)=.59*ones(n_lvalues,(n-1)); % Initialize blue channel  %.5

% Set some initial values (for MCMC sampler)
K1=diag(diff(T(2:end)).^(-1));
K2=diag(diff(T(1:(end-1))).^(-1));
A=diag(-1*ones(n-3,1),1)+diag(ones(n-4,1),2);
A=[A' ; [zeros(1,n-4) 1 -1]; [zeros(1,n-3) 1] ]'; 
B=diag(-1*ones(n-2,1))+ diag(ones(n-3,1),1);
B=[B' ; [zeros(1,n-3)  1]; [zeros(1,n-2) ] ]';
C=(K1*A-K2*B);

sigma2=nu_0/(nu_0-2)*sig_0^2;
lambda_0=a_l0/b_l0;
mu=pinv(eye(n)+lambda_0*C'*C)*Y;

disp('     Generating samples, please wait ....')
 for s=1:(nsM+nsB)

  if (mod(s,(nsM+nsB)/10)==0) 
    fprintf(1,'%2.0f percent ready \n',(s/(nsM+nsB)*100)); 
  end;
  b_s2=(nu_0*sig_0^2+(mu-Y)'*(mu-Y)+lambda_0*mu'*C'*C*mu)/2;
  sigma2=1/gamrnd((nu_0+2*n-2)/2,1/b_s2);
 
  lambda_0=gamrnd(a_l0+(n-2)/2,1/(b_l0+1/2*mu'*C'*C*mu/sigma2));
  mu=mvnrnd((pinv(eye(n)+lambda_0*C'*C)*Y)', sigma2*pinv(eye(n)+lambda_0*C'*C),1)';
  if (s>nsB) sample_0(s-nsB,:)=mu';  end;  
  
 end;
disp('     Sample generation finished. Now calculating credible regions with different values of lambda. Please wait ....')



sample_0=sample_0';

%wb = waitbar(0,'Progress');

 for smth_par=1:n_lvalues    % loop over smoothing parameters 

    disp( smth_par /n_lvalues )
    %waitbar( smth_par /n_lvalues);
    
    p_vec=zeros(n-1,1);
    
    sample=diag(diff(T))^(-1)*D*pinv(eye(n)+l_vec(smth_par)*C'*C)*sample_0;
    mean_dif=mean(sample');
    
    p_vec(find(mean_dif<0))=sum(sample(find(mean_dif<0),:)'<0)./nsM;
    p_vec(find(mean_dif>0))=sum(sample(find(mean_dif>0),:)'>0)./nsM;
    p_vec=p_vec.*(p_vec>=alpha);
    
    if meth==0                  % Independent coloring scheme

       %%%%%%%%%%%% Independent coloring scheme %%%%%%%%%%%%%%%  
      
     if Sflag==0  % The binary coloring using only two colors
       
       for r=1:(n-1)            % loop over measure points
          if (mean_dif(r)>0 & p_vec(r)>=alpha)
	       colors(smth_par,r,1)=0;
	       colors(smth_par,r,2)=0;
	       colors(smth_par,r,3)=1;
 	     end;
	    if (mean_dif(r)<0 & p_vec(r)>=alpha)
           colors(smth_par,r,1)=1;
	       colors(smth_par,r,2)=0;
	       colors(smth_par,r,3)=0;
	    end;  
       end;	
       
      else  % Coloring for maps using blue and red palette
       
        for r=1:(n-1)            % loop over measure points
          if (mean_dif(r)>0 & p_vec(r)>=alpha)
	        colors(smth_par,r,1)=((p_vec(r)-alpha)/(1-alpha))^ce*.6;
	        colors(smth_par,r,2)=((p_vec(r)-alpha)/(1-alpha))^ce*.6;
	        colors(smth_par,r,3)=1;
 	     end;
	     if (mean_dif(r)<0 & p_vec(r)>=alpha)
	       colors(smth_par,r,1)=1; 
	       colors(smth_par,r,2)=((p_vec(r)-alpha)/(1-alpha))^ce*.6;
	       colors(smth_par,r,3)=((p_vec(r)-alpha)/(1-alpha))^ce*.6;
	     end;  
       end;	

     end;   
       
   else                        % simultaneous coloring scheme
      
        %%%%%%%%%%%%%% Simultaneous inference  %%%%%%%%%%%%%%%
	
      if meth==1               % Simultaneous credible intervals I
	
	 mean_dif=mean_dif';
	 std_dif=sqrt(diag(cov(sample')));
	 stded_sample=abs(((sample-repmat(mean_dif,1,nsM)))./repmat(std_dif,1,nsM));
         d=[floor(min(max(stded_sample'))*10)/10:.025:floor(min(max(stded_sample'))*10)/10+200*.025];
     for d_ind=1:length(d);
	   if (sum(min(stded_sample<d(d_ind)))/nsM)>alpha
	      break
	   end;
	 end; 
	  
	 
	 for r=1:(n-1)         % loop over measure points
	   
	   if (mean_dif(r)-d(d_ind)*std_dif(r)>0)
	      if Sflag==0
	         colors(smth_par,r,1)=0;
	         colors(smth_par,r,2)=0;
	         colors(smth_par,r,3)=1;    
	      else       
                 colors(smth_par,r,1)=((p_vec(r)-alpha)/(1-alpha))^ce*.6;
	         colors(smth_par,r,2)=((p_vec(r)-alpha)/(1-alpha))^ce*.6;
	         colors(smth_par,r,3)=1;
	      end; 
	   end;
	   if (mean_dif(r)+d(d_ind)*std_dif(r)<0)
	      if Sflag==0
   	  	 colors(smth_par,r,1)=1;
	         colors(smth_par,r,2)=0;
	         colors(smth_par,r,3)=0;
	      else	 
  		 colors(smth_par,r,1)=1; 
	         colors(smth_par,r,2)=((p_vec(r)-alpha)/(1-alpha))^ce*.6;
	         colors(smth_par,r,3)=((p_vec(r)-alpha)/(1-alpha))^ce*.6;
	      end;	 
	    end;  
     	   
          end;	  % for 
	 
    end;	 % meth==1

    
    
    if meth==2              % highest pointwise probability method

     %%%%%%%%%%%%%%%%%%%%%%% Highest pointwise probability method  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      
 

    [b,k]=sort(-p_vec);
    for r=1:(n-1)     
   	  cum_sum=0;
	  ind_vec=zeros(length(k(1:r)),1);          % Initialize the event
          ind_vec=mean_dif(k(1:r))>0;               % Define the event; positive slopes are coded with ones, 
	                                            % negative with zeros
	  for (s_ind=1:nsM)
	    cum_sum=cum_sum+isequal(sample(k(1:r),s_ind)>0,ind_vec'); % Count the cumulative sum of favourable cases 
	                                                              % (event takes place)
          end;
	   
	  if ((cum_sum/nsM)<=alpha)
            break;                     
          else
	     if (p_vec(k(r)) >= alpha & mean_dif(k(r)) > 0 )           
	        if Sflag==0
		   colors(smth_par,r,1)=0;
	           colors(smth_par,r,2)=0;
	           colors(smth_par,r,3)=1;
	        else
 	           colors(smth_par,k(r),1)=((p_vec(k(r))-alpha)/(1-alpha))^ce*.6;
	           colors(smth_par,k(r),2)=((p_vec(k(r))-alpha)/(1-alpha))^ce*.6;
	           colors(smth_par,k(r),3)=1;
	        end;	 
             end;

             if (p_vec(k(r)) >=alpha & mean_dif(k(r)) < 0 )         
               if Sflag==0
		  colors(smth_par,r,1)=1;
	          colors(smth_par,r,2)=0;
	          colors(smth_par,r,3)=0;
	        else
	          colors(smth_par,k(r),1)=1; 
	          colors(smth_par,k(r),2)=((p_vec(k(r))-alpha)/(1-alpha))^ce*.6;
	          colors(smth_par,k(r),3)=((p_vec(k(r))-alpha)/(1-alpha))^ce*.6;
	        end; 
             end;    
	 end; % ((cum_sum/ns)<=alpha) 
	
    end;                  % for r=1:(n-1)
      
  end;        % of method 2  
      
 end;   % coloring scheme selection (independent/simultaneous)

 end; % smoothing parameter

 
 % Drawing the output

clf;
if min(l_vec)==0;
    temp=find(l_vec>0);
    min_l=l_vec(temp(1));
 else
    min_l=min(l_vec);
end;

if Sflag==1

  col_inf_R(:,:,1)=ones(1,101);
  col_inf_R(:,:,2)=[0:.01:1].^ce*.6;
  col_inf_R(:,:,3)=[0:.01:1].^ce*.6;
  col_inf_B(:,:,1)=[0:.01:1].^ce*.6;
  col_inf_B(:,:,2)=[0:.01:1].^ce*.6;
  col_inf_B(:,:,3)=ones(1,101);
  
  subplot('Position',[.1 .25 .8 .65])
  image(T,[log10(min_l) log10(max(l_vec))],colors)
  set(gca,'YDir','normal')
  ylabel('log_{10}(\lambda)','FontSize',13);
  xlabel('Data index')
  if (meth==0)
    title(['Bayesian SiZer map - Pointwise credible intervals. \alpha= ',num2str(alpha)])
  end;
  if (meth==1)
    title(['Bayesian SiZer map - Simultaneous credible intervals I. \alpha= ',num2str(alpha)])
  end;
  if (meth==2)
    title(['Bayesian SiZer map - Simultaneous credible intervals II. \alpha= ',num2str(alpha)])
  end;  
  

  subplot('Position',[.52 .07 .37 .03])
  image(col_inf_R);
  set(gca,'YTick',[])
  ticks=[1 20 40 60 80 101];
  set(gca,'XTick',ticks);
  set(gca,'XTickLabel',[alpha+[0:.2:1].*(1-alpha)])
  
  subplot('Position',[.1 .07 .37 .03])
  image(col_inf_B);
  set(gca,'YTick',[])
  set(gca,'XTick',ticks);
  set(gca,'XTickLabel',[alpha+[0:.2:1].*(1-alpha)])
  
else
  
  image(T,[log10(min_l) log10(max(l_vec))],colors)
  set(gca,'YDir','normal')
  ylabel('log_{10}(\lambda)','FontSize',13);
  xlabel('Data index')
  if (meth==0)
         title(['Bayesian SiZer map - Pointwise credible intervals. \alpha= ',num2str(alpha)])
  end; 
  if (meth==1)
         title(['Bayesian SiZer map - Conservative simultaneous credible intervals. \alpha= ',num2str(alpha)])
  end;
  if (meth==2)
         title(['Bayesian SiZer map - Greedy simultaneous credible intervals. \alpha= ',num2str(alpha)])
  end;  
  
end;
 
 
%close(wb);
disp('     Finished!')