function [mst,tg] = poisson_process(x, tao, prs, t, ...
    nodes, nfig, mthd, bsz)
% mthd:0 -> Quadrature rule
%      1 -> Importance sampling 
%      2 -> Importance sampling with Qr proposals
%      3-> Grid method

if nargin < 5, nodes = 200; end % numb. of nodes or samples
if nargin < 6, nfig = 1; end
if nargin < 7, mthd = 1; end 
if nargin < 8, bsz = 1; end 
wshift = 1;

pi2 = 2*pi;
ind_upd = 0;

[n,d] = size(x);
td = t(2)-t(1);
A = eye(1);
C= eye(1);
mst0 = prs.m0;
vst0 = prs.v0;
Lt0 = inv(vst0);

ns = 1;
% ng = floor(n/bsz);
ng = 1 + floor( (n-bsz)/wshift );
mst = zeros(ns, ng);
vst = zeros(ns,ns,ng);
vst(:,:,1) = 1e3*eye(ns);

muo0=zeros(1,1); 
muot=zeros(1,ng);
ho0=1e-8*eye(1);
hot = 1e2*ones(1,1,ng);
hot(:,:,1)=ho0;

global wqr prnodes tsp;
tsp = 0;% 1:sample from t distri; 0: sample from normal
if mthd == 0
    [ndps, w] = quadrl(nodes); % ndps is row vec.
elseif mthd == 1
    ndps =[];
    w = ones(nodes,1); % for the directly smapling method
else
    prnodes  = 7;
    [ndps, wqr] = quadrl(prnodes); % ndps is row vec.
    w = ones(nodes,1); % for the directly smapling method
end

tg = 0;
g = 0;
tg(1) = t(1);
for i = 1:wshift:(ng-1)*wshift +1
    g = g+1;
    tg(g) = t(i+floor(bsz/2));
  
    if mod(i,200)== 0, fprintf('\n process %d data point', i); end 
     % step A: compute forward message
    if g > 1 
         mst(:,g) =  A*mst(:,g-1);
         vst(:,:,g) = A*vst(:,:,g-1)*A' + tao;
     else
         mst(:,g) =  mst0;
         vst(:,:,g) = vst0;  
     end
    mstold = mst(:,g);
    vstold = vst(:,:,g);
    % Step B.C.D: compute the moments based on new observation
    % Only good for bsz =1's case   
    g=g
    [mstn, vstn] = upd_moments_psn_obs(mstold, ...
     vstold, mst(:,g), vst(:,:,g), x(g), nodes);
    
     % step E and F: compute the observation message and new belief
     [muot(:,g),hot(:,:,g),mst(:,g), vst(:,:,g)]...
        =upd_mssg_belief( mstold, ...
        vstold, mstn, vstn, muo0, ho0,1);  
        
     Lt(:,:,g) = inv(A*vst(:,:,g)*A' + tao);
end
g_old = g;

if nfig == 1    
    figure(1),
    plot(tg,mst);
    hold on;
end     

niter = 5;

for k = 2:niter
    
    mdfall = 0;
    % step G and H: kalman smoothing and obs. message updating
    g = g_old - 1;
    for i = ((ng-2)*wshift+1):-wshift:1  
       if mod(i,200)== 0, fprintf('\n k:%d  i:%d', k,i); end 
    
       mstki =  mst(:,g);
       % compute the difference 
       mdfall = max( abs(mst(:,g) - mstki)) + mdfall; 
               
       % G step: Kalman smoothing
       vsAt = vst(:,:,g)*A';
       Jt = vsAt * Lt(:,:,g);    
        mst(:,g) = mst(:,g) + Jt*(mst(:,g+1) - A*mst(:,g));
       % vst(:,:,g) = vst(:,:,g)*(ee+ (Lt(:,:,g)*vst(:,:,i+1) - ee)*Jt');
       vst(:,:,g) = vst(:,:,g) + Jt*( vst(:,:,g+1)*Jt' ...
           - vsAt');
         
        %  H step 
        % first remove the term from q to avoid double counts
        ivst = inv(vst(:,:, g));
        vstold = inv(ivst - hot(:,:,g));
        mstold = vstold * ( ivst*mst(:,g) - muot(:,g) );
        
%        if isposdef( vstold ) == 0,  
        if  min(eig(vstold)) <= 0
           wstr = ['vstold: not positive definite at backward iter.', num2str(k), 'and ', num2str(i) 'th dt point'];
           warning(wstr);
           g = g-1;
           continue;
        end 
        
       % ADF update: compute the moments based on observation message
       [mstn, vstn] = upd_moments_psn_obs(mstold, ...
         vstold, mst(:,g), vst(:,:,g), x(g), nodes);
     
       % Update the obs message and the new belief       
        muot_old = muot(:,g);
        hot_old = hot(:,:,g);
        [muot(:,g), hot(:,:,g), mst(:,g),vst(:,:,g)]= ...
        upd_mssg_belief( mstold, vstold, ...
        mstn, vstn,muot_old,hot_old);  
      
        g = g-1;
   end % end of F and G steps: smoothing and message updating
   
   if nfig == 1    
       figure(1),
       plot(tg,mst);
       hold on;
   end     
   
%   if mdfall/n < 1e-5, 
%        fprintf('\nep converge at iteration %d', k); 
%       break; 
%   end
   if k == niter, break; end
   
    % Step I: Kalman filtering
    % vst(:,:,1) = inv( Lt0 + C'*hot(:,:,1)*C);
    % mst(:,1) =  vst(:,:,1)*( C'*muot(:,1) + Lt0*A*mst0);     
        
    for g = 2:g_old  
       Lt(:,:,g-1) = inv(A*vst(:,:,g-1)*A' + tao);
       mstold = mst(:,g);
       vstold = vst(:,:,g);
       vst(:,:,g) = inv( Lt(:,:,g-1) + C'*hot(:,:,g)*C);
%       K = ee - vst(:,:,g)*Lt(:,:,i-1);
       mst(:,g) = vst(:,:,g)*(  C'*muot(:,g) + Lt(:,:,g-1)*A*mst(:,g-1) );
       if min(eig(vst(:,:,g)) ) <= 0
           mst(:,g)= mstold;
           vst(:,:,g) = vstold;                   
           wstr = ['vstnew:not positive definite at forward iter.', num2str(k), 'and ', num2str(g) 'th dt point'];
           warning(wstr);
       end
    end       
end    
