function colors = bsizer(Y,T,nu_0,sig_0,lam_0,meth,alpha,Sflag,l_vec,ns)
%
% function MAP = bsizer(Y,T,nu_0,sig_0,lam_0,meth,alpha,Sflag,l_vec,alpha,ns)
%
% Bayesian Sizer - fixed lambda_0
%                 
%
% 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 - scale parameter sigma_0 for the prior of sigma^2 (NOTE: not a squared sigma_0!)
% lam_0 - scalar lambda_0 (fixed)
% 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)]
% ns    - number of samples in simulation (required in meth={1,2})
%         [OPTIONAL: default 1000]
%
% 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(5,10,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 lam_0<0 
  error('Parameter lam_0 (lambda_0) must be non-negative')
end;

if nargin <= 5 ;    %  
  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 <= 6 ;    
  alpha_u = 0.8 ;
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 <= 7 ;    
  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 <= 8 ;    %  
  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 <= 9 ;    %  
  ns_u = 3000 ;         
else ;              %  ns was specified, use that
  ns_u = ns ;
end ;
ns=ns_u;
clear ns_u;

if ns<=1 | (ns-ceil(ns)>0)
  error('number of samples must be positive integer')
end;  


% Let 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 
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);

S_0=pinv(eye(n)+lam_0*C'*C);
const=(nu_0*sig_0^2+Y'*Y-Y'*S_0*Y)/(nu);
Sigma_0=const*S_0;





if (meth ==1 | meth==2)  % Simulate only in the case of simultaneous maps
    % Perform whitening for the t-sampler (that needs diagonal scale matrix)
    mean_0_T=Sigma_0^(-.5)*(S_0*Y);         % Whiten the mean
    sample_0=t_rand(nu,mean_0_T,eye(n),ns);
    sample_0=Sigma_0^(.5)*sample_0;         % Do inverse whitening
    meanS=mean(sample_0');
    meanS=real(meanS);                      % This is for some versions of Matlab (*)
    covS=cov(sample_0');                    
    covS=real(covS);                        % This is for some versions of Matlab (*)
end;  

% (*) Some versions of the Matlab (e.g. Matlab 6.0) seems to behave numerically
% in an unsatisfactory way; suggest a imaginary matrix when the truth is real

disp('     Calculating, please wait ....')
%wb = waitbar(0,'Progress');

 for smth_par=1:n_lvalues    % loop over smoothing parameters 
    disp( smth_par /n_lvalues )
    %waitbar( smth_par /n_lvalues);

    S_l=pinv(eye(n)+l_vec(smth_par)*C'*C);

    mu_bar=S_l*S_0*Y;
    Sigma=S_l*Sigma_0*S_l;

    % Next calculate the parameters of the joint difference of difference quotients
    Upsilon=D*Sigma*D';
    g_eta=diff(mu_bar)./diff(T');               % Posterior mean of difference quotients
    g_Upsilon=Gamma'*Upsilon*Gamma;             % Posterior covariance matrix of diff. quotients.
    

    quant=zeros(n-1,1);
    
    if (meth ==1 | meth==2)  % Do this only in the case of simultaneous maps
    
      mean_dif=diag(diff(T))^(-1)*D*pinv(eye(n)+l_vec(smth_par)*C'*C)*meanS'; 
      cov_mat=diag(diff(T))^(-1)*D*pinv(eye(n)+l_vec(smth_par)*C'*C)*covS*...
	      pinv(eye(n)+l_vec(smth_par)*C'*C)*D'*diag(diff(T))^(-1);
      std_dif=sqrt(diag(cov_mat));
      sample=diag(diff(T))^(-1)*D*pinv(eye(n)+l_vec(smth_par)*C'*C)*sample_0;
      sample=real(sample);
    
      % Quantities calculated from t-distribution
      z = betainv(2*(1-alpha),nu/2,0.5);
      dist_from_zero = sqrt(nu.*(std_dif.^2) ./ z - nu.*(std_dif.^2));
      quant(mean_dif<0)=mean_dif(mean_dif<0)+dist_from_zero(mean_dif<0);
      quant(mean_dif>0)=mean_dif(mean_dif>0)-dist_from_zero(mean_dif>0);
     
      
      % The following normal approximations can also be used 
      %quant(mean_dif<0)= sqrt(2) * std_dif(mean_dif<0) .* erfinv(2 * alpha - 1) + mean_dif(mean_dif<0);
      %quant(mean_dif>0)=-sqrt(2) * std_dif(mean_dif>0) .* erfinv(2 * alpha - 1) + mean_dif(mean_dif>0);
    
    else    

      % Quantities calculated from t-distribution
      z = betainv(2*(1-alpha),nu/2,0.5);
      dist_from_zero = sqrt(nu.*diag(g_Upsilon) ./ z - nu.*diag(g_Upsilon));
      quant(g_eta<0)=g_eta(g_eta<0)+dist_from_zero(g_eta<0);
      quant(g_eta>0)=g_eta(g_eta>0)-dist_from_zero(g_eta>0);
    
      % The following normal approximations can be used 
      %quant(g_eta<0)= sqrt(2) * std_dif(g_eta<0) .* erfinv(2 * alpha - 1) + g_eta(g_eta<0);
      %quant(g_eta>0)=-sqrt(2) * std_dif(g_eta>0) .* erfinv(2 * alpha - 1) + g_eta(g_eta>0);

    end;
    
   
   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 (g_eta(r)>0 & quant(r)>0)
	     colors(smth_par,r,1)=0;
	     colors(smth_par,r,2)=0;
	     colors(smth_par,r,3)=1;
 	   end;
	  if (g_eta(r)<0 & quant(r)<0)
	     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 (g_eta(r)>0 & quant(r)>0)
             z = 1 ./ (1 + g_eta(r).^2/(nu*g_Upsilon(r,r)));
             p(r) = 1-betainc(z, nu/2, 0.5)/2;
	     colors(smth_par,r,1)=((p(r)-alpha)/(1-alpha))^ce*.6;
	     colors(smth_par,r,2)=((p(r)-alpha)/(1-alpha))^ce*.6;
	     colors(smth_par,r,3)=1;
 	  end;
	  if (g_eta(r)<0 & quant(r)<0)
             z = 1 ./ (1 + g_eta(r).^2/(nu*g_Upsilon(r,r)));
             p(r) = 1-betainc(z, nu/2, 0.5)/2;
	     colors(smth_par,r,1)=1; 
	     colors(smth_par,r,2)=((p(r)-alpha)/(1-alpha))^ce*.6;
	     colors(smth_par,r,3)=((p(r)-alpha)/(1-alpha))^ce*.6;
	  end;  
       end;	
     
    end;
       
   else                        % simultaneous coloring scheme
      
      
       %%%%%%%%%%%%%% Simultaneous inference %%%%%%%%%%%%%%%
	
    if meth==1               % Simultaneous credible intervals I
	 
         %sample=diag(diff(T))^(-1)*D*pinv(eye(n)+l_vec(smth_par)*C'*C)*sample_0;
	 stded_sample=abs(((sample-repmat(mean_dif,1,ns)))./repmat(std_dif,1,ns));
         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)))/ns)>alpha
	      break
	   end;
	 end; 
	  
	 mean_dif=mean_dif';
	 
	 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
 	        z = nu*std_dif(r).^2 ./ (nu*std_dif(r).^2+mean_dif(r).^2);
                temp = 1 - betainc(z, nu/2, 0.5);
                p=1-(1-temp)/2;
	        colors(smth_par,r,1)=((p-alpha)/(1-alpha))^ce*.6;
	        colors(smth_par,r,2)=((p-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	 
		 z = nu*std_dif(r).^2 ./ (nu*std_dif(r).^2+mean_dif(r).^2);
                 temp = 1 - betainc(z, nu/2, 0.5);
                 p=1-(1-temp)/2;
		 colors(smth_par,r,1)=1; 
	         colors(smth_par,r,2)=((p-alpha)/(1-alpha))^ce*.6;
	         colors(smth_par,r,3)=((p-alpha)/(1-alpha))^ce*.6;
	       end;	 
	   end;  
     	   
         end;	 
	 
    end;	 % meth==1

    
    
    if meth==2              % highest pointwise probability method

     %%%%%%%%%%%%%%%%%%%%%%% Highest pointwise probability method  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      
         %sample=diag(diff(T))^(-1)*D*pinv(eye(n)+l_vec(smth_par)*C'*C)*sample_0;
         mean_dif=mean_dif';
         p(1:(n-1))=zeros(1,(n-1));

         for r=1:(n-1)         % loop over measure points
              
	    if ( mean_dif(r) > 0 & quant(r) > 0 )
                z = nu*std_dif(r).^2 ./ (nu*std_dif(r).^2+mean_dif(r).^2);
                temp = 1 - betainc(z, nu/2, 0.5);
                p(r)=1-(1-temp)/2;
 	    end;
	     
	    if ( mean_dif(r) < 0 & quant(r) < 0 )
                z = nu*std_dif(r).^2 ./ (nu*std_dif(r).^2+mean_dif(r).^2);
                temp = 1 - betainc(z, nu/2, 0.5);
                p(r)=1-(1-temp)/2;
	    end;

        end;                  % of for (loop over measure points)
       

       [b,k]=sort(-p);
       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:ns)
	     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/ns)<=alpha)
	      break;                     
          else
	     if (quant(k(r)) > 0 & mean_dif(k(r)) > 0 )           
	        if Sflag==0
		  colors(smth_par,k(r),1)=0;
	          colors(smth_par,k(r),2)=0;
	          colors(smth_par,k(r),3)=1;
	        else
 	          colors(smth_par,k(r),1)=((p(k(r))-alpha)/(1-alpha))^ce*.6;
	          colors(smth_par,k(r),2)=((p(k(r))-alpha)/(1-alpha))^ce*.6;
	          colors(smth_par,k(r),3)=1;
	        end;	 
             end;

             if (quant(k(r)) < 0 & mean_dif(k(r)) < 0 )         
                if Sflag==0
	          colors(smth_par,k(r),1)=1;
	          colors(smth_par,k(r),2)=0;
	          colors(smth_par,k(r),3)=0;
	        else
	          colors(smth_par,k(r),1)=1; 
	          colors(smth_par,k(r),2)=((p(k(r))-alpha)/(1-alpha))^ce*.6;
	          colors(smth_par,k(r),3)=((p(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')
  %set(gca,'XTick',[0:12:length(Y)+3])
  %set(gca,'XTickLabel',['1980';'    ';'1982';'    ';'1984';'    ';'1986';'    ';'1988';'    ';'1990';'    ';'1992';'    '; ])
  ylabel('log_{10}(\lambda)','FontSize',13);
  xlabel('Data index')
  if (meth==0)
    title(['Bayesian SiZer map - Pointwise credible intervals. \alpha= ',num2str(alpha), ' \lambda_0 =',num2str(lam_0)])
  end;
  if (meth==1)
    title(['Bayesian SiZer map - Simultaneous credible intervals I. \alpha= ',num2str(alpha), ' \lambda_0 =',num2str(lam_0)])
  end;
  if (meth==2)
    title(['Bayesian SiZer map - Simultaneous credible intervals II. \alpha= ',num2str(alpha), ' \lambda_0 =',num2str(lam_0)])
  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), ' \lambda_0 =',num2str(lam_0)])
  end; 
  if (meth==1)
         title(['Bayesian SiZer map - Conservative simultaneous credible intervals. \alpha= ',num2str(alpha), ' \lambda_0 =',num2str(lam_0)])
  end;
  if (meth==2)
         title(['Bayesian SiZer map - Greedy simultaneous credible intervals II. \alpha= ',num2str(alpha), ' \lambda_0 =',num2str(lam_0)])
  end;  
  
end;

 
 
 
%close(wb);
disp('     Finished!')