rjmcmc_sinu3(y , hyper , x);

function [out , param ]  = rjmcmc_sinu3(y , hyper , x);

%
%[out , param ]  = rjmcmc_sinu3(y , hyper , x);
%

[N , c]                = size(y);

if ( c ~= 1)
    
    error('y must be (N x 1) ');
    
end

if (nargin < 2 || isempty(hyper))
    
    hyper.T        = 5000;
    
    hyper.burnin   = hyper.T/3;
    
    hyper.kmax     = 12;
    
    hyper.c        = 0.5;
    
    hyper.O        = 12;
    
    
    hyper.upsilon  = 0.1;
    
    hyper.gamma    = 0.1;
    
    hyper.epsi1    = 10e-4;
    
    hyper.epsi2    = 10e-4;
    
    hyper.alpha    = 2;
    
    hyper.beta     = 10;
    
    hyper.lambda   = 0.2;
    
    hyper.merge    = 0.1;
    
    hyper.sRW      = (1/(5*N))^2;
    
    hyper.walk     = 0.0;
    
end

if (nargin < 3 || isempty(x))
    
    x                      = (1 : N)';
    
end


%% ---- Outputs memory allocation -------%%%%

k                      = zeros(1 , hyper.T);

wk                     = zeros(hyper.kmax , hyper.T);

ak                     = zeros(hyper.kmax , hyper.T);


nabla                  = zeros(1 , hyper.T);       % Poisson parameter.

delta                  = zeros(1 , hyper.T);       % Regularisation parameter.

sigma                  = zeros(1 , hyper.T);       % Output noise variance.


out.arwupdate          = 0;

out.rrwupdate          = 0;

out.abirth             = 0;

out.rbirth             = 0;

out.adeath             = 0;

out.rdeath             = 0;

out.amerge             = 0;

out.rmerge             = 0;

out.asplit             = 0;

out.rsplit             = 0;

%% ------ Usefull matrices & vectors --------------%%%%

IN                     = eye(N);

yt                     = y';


hvN                    = (hyper.upsilon + N)*0.5;

epsi1                  = 0.5 + hyper.epsi1;

epsi2                  = 1 + hyper.epsi2;

walkint                = pi*(1 + 2*hyper.walk);

S                      = hyper.O*N;


%%--- Compute periodogram with eventually zero-padding -----%%

temp_fft_y             = abs(fft(y , S)).^2;

perio_y                = temp_fft_y( 1 : S/2 );

cdfperio_y             = cumsum(perio_y/sum(perio_y));    %Cummulative density function of the normalized periodogram 

cdfperio_y             = [0 ; cdfperio_y(1:end - 1)];

w                      = (0 : S/2 - 1)'*(2*pi/S);          % build the frequency vector


%%------ Initial draw for hyperparameters  -----------%%%

sigma(1)               = 1./gammarnd(hyper.upsilon/2 , hyper.gamma/2);

delta(1)               = 1./gammarnd(hyper.alpha , hyper.beta);

nabla(1)               = gammarnd(epsi1 , epsi2 );




k(1)                   = hyper.kmax;

ind_k                  = (1 : k(1));

I1                     = eye(k(1));

decision               = rand(1 , hyper.T - 1);



%%----------Initial drawn for amplitude (supposed centered) ---%%

ak(ind_k , 1)         = chol(sigma(1)*I1).'*randn(k(1) , 1);

%%----------Initial drawn for pulsations [0 , pi] ---%%

tp_wk                 = pi*rand(k(1) , 1);

wk(ind_k , 1)         = sort(tp_wk);


%%------ Regression matrix (N x 2k(t)) -------%%

D                     = bsinu(wk(ind_k , 1) , x);

%%------ Initialization -------%%


Dt                    = D';

DtD                   = Dt*D;

M                     = inv(DtD + (1/delta(: , 1))*I1);

MDt                   = M*Dt;

ytPy                  = (yt*(IN - D*MDt))*y;


for t = 1 : hyper.T - 1
    
    %disp(sprintf('iterations = %1.0f/%1.0f' , t , hyper.T));
    
    birth    = hyper.c*min(1 , (nabla(t)/(k(t) + 1)));         % p(k + 1)/p(k) = nabla(t)/(k(t) + 1)
    
    death    = hyper.c*min(1 , ((k(t) + 1)/nabla(t)));
    
    %%--------- Detection part => Reverse Jump for estimate k --%%
    
    if ( (decision(t) <= birth) && (k(t) < hyper.kmax) )
        
        [k1 , tp_wk , D , DtD , M , MDt , ytPy , out]  = birth_move(y , yt , x , tp_wk , D , DtD , M , MDt , ytPy , delta , t , hyper , k  , IN , out);
        
    elseif ( (decision(t) <= (birth + death)) && (k(t) > 0) )
        
        [k1 , tp_wk , D , DtD , M , MDt , ytPy , out ] = death_move(y , yt , x , tp_wk , D , DtD , M , MDt , ytPy , delta , t , hyper , k  , IN , out);
        
        
    elseif ( (decision(t) <= 2*birth + death) && (k(t) < hyper.kmax) && ( k(t) > 1 ) )
        
        [k1 , tp_wk , D , DtD , M , MDt , ytPy , out ] = split_move(y , yt , x , tp_wk , D , DtD , M , MDt , ytPy , delta , walkint , t , hyper , k  , IN , out);        
        
    elseif ( ( decision(t) <= 2*(birth + death)) && ( k(t) > 1 ) )
        
        [k1 , tp_wk , D , DtD , M , MDt , ytPy , out ] = merge_move(y , yt, x , tp_wk , D , DtD , M , MDt , ytPy , delta  , t , hyper , k  , IN , out);
        
        
    else
        
        %%--------- Estimation part a) => Estimation of wk --%%
        
        
        %%========= Metropolis Hasting with 2 steps ============%% 
        
        [k1 , tp_wk , D , DtD , M , MDt , ytPy , out] = rwupdate_move(y , yt , x , cdfperio_y , w , tp_wk , D , DtD , M , MDt , ytPy , delta , t , hyper , k , IN , out);            
                
    end
        
    k(t + 1)           = k1;
    
    ind_k              = (1 : k1);
        
    wk(ind_k , t + 1)  = sort(tp_wk);
    
    sigma(t + 1)       = 1/gammarnd( hvN , 0.5*(hyper.gamma + ytPy)  );        % (eq 17) %
    
    %%----------- Estimation part b) => Estimation of ak ----------------------%%
    
    ak(ind_k , t + 1)  = MDt*y + chol(sigma(t + 1)*M).'*randn(k1 , 1);         % (eq 17) %
    
    %%============ Update hyper-parameters with Gibbs sampler ==============%%
    
    delta(t + 1)       = 1/gammarnd(k1 + hyper.alpha , ((ak(ind_k , t + 1)'*DtD)*ak(ind_k , t + 1))/(2*sigma(t + 1)) + hyper.beta );
    
    nabla(t + 1)       = gammarnd(epsi1 + k1 , epsi2 );                        %eq (24)
    
    
end

%%============= Outputs =============%%

[maxi , k_est]            = max(histc(k(k < hyper.burnin) , 0:hyper.kmax));

k_est                     = k_est - 1;

ind_k_est                 = find( k == k_est );

ind_k_est((ind_k_est < hyper.burnin)) = [];


wk_est                    = mean(wk(: , ind_k_est) , 2);

wk_est((wk_est == 0))     = [];

ak_est                    = mean(ak(: , ind_k_est) , 2);

ak_est((ak_est == 0))     = [];

sig_est                   =  mean(sigma(ind_k_est));


out.k_est                 = k_est;

out.wk_est                = wk_est;

out.ak_est                = ak_est;

out.sig_est               = sqrt(sig_est);


param.k                   = k;

param.ak                  = ak;

param.wk                  = wk;

param.delta               = delta;

param.sigma               = sigma;

param.nabla               = nabla;