function [mean_smooths,colors] = mw_bsizer(X,T,meth,alpha,Sflag,l_vec)

% Model within BSiZer

% INPUTS:
% X     - sampled realizations of the signal (n*nsM -matrix [nsM=number of samples])
% T     - sampled realizations of the explanatory variable (n*nsM -matrix).
%         If fixed, give a repeated matrix (e.g. repmat([1:n]',1,nsM))
% meth  - 0 independent features
%         1 simultaneous features (conservative version)
%         2 simultaneous features using highest pointwise probabilities
%         3 twice simultaneous (simultaneous in x and lambda scales)
%         [OPTIONAL: default 0]
% alpha - Confidence level (0.5 < alpha < 1)
%         [OPTIONAL: default 0.8]
% Sflag - 0 - binary coloring
%         1 - color maps with shades based on posterior probabilities from independent features
%         [OPTIONAL: default 0]
% l_vec - column vector containing the smoothing parameters in increasing order
%         [OPTIONAL: default logspace(log10(.5),6,80)]
%
% OUTPUTS:
% colors - Bayesian Sizer map
% mean_smooths - Posterior means of scaled mu (mean smooths with different levels of smoothing)



[n,nsM]=size(X);
[nT,nsT]=size(T);



error(nargchk(2,8,nargin))

if n~=nT | nsM~=nsM
    error('X and T must be of the same size')
end;

if nsM<100
    disp('Sample size seems to be quite small. Are you sure you want to continue?')
    disp('Press a key to continue. Ctrl-C to cancel')
end;


if nargin <= 2 ;    %
    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 <= 3;
    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 <= 4 ;
    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 <= 5 ;
    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;


T_mean=mean(T');

n_lvalues=length(l_vec);

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
colors(:,:,2)=.42*ones(n_lvalues,(n-1)); % Initialize green channel
colors(:,:,3)=.59*ones(n_lvalues,(n-1)); % Initialize blue channel


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) ] ]';





for i=1:nsM
    disp(i)
    K1=diag(diff(T(2:end,i)).^(-1));
    K2=diag(diff(T(1:(end-1),i)).^(-1));
    C=(K1*A-K2*B);

    for smth_par=1:n_lvalues
        sample_0(:,i,smth_par)=diag(diff(T(:,i)))^(-1)*D*pinv(eye(n)+l_vec(smth_par)*C'*C)*X(:,i);
        smooth(:,i,smth_par)=pinv(eye(n)+l_vec(smth_par)*C'*C)*X(:,i);
    end;
end;
mean_smooths=mean(smooth,2);

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

for smth_par=1:n_lvalues    % loop over smoothing parameters

    disp( smth_par /n_lvalues )

    p_vec=zeros(n-1,1);
    sample=sample_0(:,:,smth_par);
    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  % 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)=1;
                    colors(smth_par,r,2)=0;
                    colors(smth_par,r,3)=0;
                end;
                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;
            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)=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;
                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;
            end;

        end;

    else                        % simultaneous coloring scheme

        %%%%%%%%%%%%%% Simultaneous inference  %%%%%%%%%%%%%%%

        if meth==1               % Simultaneous credible intervals I (conservative version)

            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)=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;
                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;

            end;

        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,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_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;

                    if (p_vec(k(r)) >=alpha & 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_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;
                end; % ((cum_sum/ns)<=alpha)

            end;                  % for r=1:(n-1)

        end;        % of method 2

    end;   % coloring scheme selection (0,1,2)

end; % smoothing parameter

if meth==3

    p_vec=zeros((n-1)*n_lvalues,1);

    for smth_par=1:n_lvalues    % loop over smoothing parameters
        sample=sample_0(:,:,smth_par);
        mean_dif=mean(sample');

        p_vec((smth_par-1)*(n-1)+find(mean_dif<0))=sum(sample(find(mean_dif<0),:)'<0)./nsM;
        p_vec((smth_par-1)*(n-1)+find(mean_dif>0))=sum(sample(find(mean_dif>0),:)'>0)./nsM;
        p_vec=p_vec.*(p_vec>=alpha);

        sample_vec(((smth_par-1)*(n-1)+1):((smth_par)*(n-1)),:)=sample;
        mean_dif_vec(((smth_par-1)*(n-1)+1):((smth_par)*(n-1)))=mean(sample');

    end;

    [b,k]=sort(-p_vec);

    for r=1:((n-1)*n_lvalues)

        cum_sum=0;
        ind_vec=zeros(length(k(1:r)),1);            % Initialize the event
        ind_vec=mean_dif_vec(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_vec(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_vec(k(r)) > 0 )
                if Sflag==0
                    colors(floor((k(r)-1)/(n-1))+1,mod(k(r)-1,(n-1))+1,1)=1;
                    colors(floor((k(r)-1)/(n-1))+1,mod(k(r)-1,(n-1))+1,2)=0;
                    colors(floor((k(r)-1)/(n-1))+1,mod(k(r)-1,(n-1))+1,3)=0;
                else
                    colors(floor((k(r)-1)/(n-1))+1,mod(k(r)-1,(n-1))+1,1)=1;
                    colors(floor((k(r)-1)/(n-1))+1,mod(k(r)-1,(n-1))+1,2)=((p_vec(k(r))-alpha)/(1-alpha))^ce*.6;
                    colors(floor((k(r)-1)/(n-1))+1,mod(k(r)-1,(n-1))+1,3)=((p_vec(k(r))-alpha)/(1-alpha))^ce*.6;
                end;
            end;

            if (p_vec(k(r)) >=alpha & mean_dif_vec(k(r)) < 0 )
                if Sflag==0
                    colors(floor((k(r)-1)/(n-1))+1,mod(k(r)-1,(n-1))+1,1)=0;
                    colors(floor((k(r)-1)/(n-1))+1,mod(k(r)-1,(n-1))+1,2)=0;
                    colors(floor((k(r)-1)/(n-1))+1,mod(k(r)-1,(n-1))+1,3)=1;
                else
                    colors(floor((k(r)-1)/(n-1))+1,mod(k(r)-1,(n-1))+1,1)=((p_vec(k(r))-alpha)/(1-alpha))^ce*.6;
                    colors(floor((k(r)-1)/(n-1))+1,mod(k(r)-1,(n-1))+1,2)=((p_vec(k(r))-alpha)/(1-alpha))^ce*.6;
                    colors(floor((k(r)-1)/(n-1))+1,mod(k(r)-1,(n-1))+1,3)=1;
                end;
            end;

        end; % ((cum_sum/ns)<=alpha)

    end;                  % for r=1:(n-1)


end;  % meth==2


% Drawing the output

for smth_par=1:n_lvalues    % loop over smoothing parameters
    xl(smth_par)=n/2-4*sqrt(sqrt(l_vec(smth_par)));
    xr(smth_par)=n/2+4*sqrt(sqrt(l_vec(smth_par)));
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_mean,[log10(min(l_vec)) log10(max(l_vec))],colors)
    set(gca,'YDir','normal')
    ylabel('log_{10}(\lambda)','FontSize',13);
    xlabel('Data index')
    if (meth==0)
        title(['Pointwise credible intervals. \alpha= ',num2str(alpha)])
    end;
    if (meth==1)
        title(['Simultaneous credible intervals (conservative). \alpha= ',num2str(alpha)])
    end;
    if (meth==2)
        title(['Simultaneous credible intervals. \alpha= ',num2str(alpha)])
    end;
    if (meth==3)
        title(['Twice simultaneous credible intervals. \alpha= ',num2str(alpha)])
    end;
    hold on
    plot(xl*(max(mean(T'))-min(mean(T')))/n,log10(l_vec),'w',xr*(max(mean(T'))-min(mean(T')))/n,log10(l_vec),'w');
    hold off



    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_mean,[log10(min(l_vec)) log10(max(l_vec))],colors)
    set(gca,'YDir','normal')
    ylabel('log_{10}(\lambda)','FontSize',13);
    xlabel('Years before present')
    if (meth==0)
        title(['Pointwise credible intervals. \alpha= ',num2str(alpha)])
    end;
    if (meth==1)
        title(['Simultaneous credible intervals (conservative). \alpha= ',num2str(alpha)])
    end;
    if (meth==2)
        title(['Simultaneous credible intervals. \alpha= ',num2str(alpha)])
    end;
    if (meth==3)
        title(['Twice simultaneous credible intervals. \alpha= ',num2str(alpha)])
    end;
    hold on
    plot(xl*(max(mean(T'))-min(mean(T')))/n,log10(l_vec),'w',xr*(max(mean(T'))-min(mean(T')))/n,log10(l_vec),'w');
    hold off

end;



disp('     Finished!')