function mst = poisson_process2(x, tao, prs, t, ...
    nodes, mthd, nfig, 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, mthd = 1; end 
if nargin < 7, nfig = 1; end
if nargin < 8, bsz = 1; end 
wshift = 1;

global scl;
pi2 = 2*pi;
ind_upd = 0;

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

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

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

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

% lnty =['b' 'g' 'r' 'm' 'c' 'k'];
lnty =['b' 'g' 'r' 'm' 'k'];
niter = 5;
% lr = 1;
for k = 1:niter
  % lr = linspace(ng/7,ng,ng)/ng/3; lr = lr.^.5;
  lr = .3*ones(ng,1);
  for g = 1:ng
    if mod(g,200)== 0, fprintf('\n process %d data point', g); end 
     % step A: compute forward message
     if g > 1 
         mst_pr = mst(:,g); 
         vst_pr = vst(:,:,g); 
         mst(:,g) =  A*mst(:,g-1);
         vst(:,:,g) = A*vst(:,:,g-1)*A' + tao;
     else
         if k == 1
             mst(:,g) =  mst0;
             vst(:,:,g) = vst0;  
         else
             continue;
         end
     end
     
     if k == 1
         ivst = inv(vst(:,:, g));
         vstold = inv(ivst - hot(:,:,g));
         mstold = vstold * ( ivst*mst(:,g) - muot(:,g) );
         if  min(eig(vstold)) <= 0
             wstr = ['vstold: not positive definite at forward iter.', num2str(k), 'and ', num2str(g) 'th dt point'];
             warning(wstr);
             mst(:,g) = mst_pr; 
             vst(:,:,g) = vst_pr; 
             continue;
         end 
        % Step B.C.D: compute the moments based on new observation
        [mstn, vstn] = upd_moments_psn_obs(mstold, ...
         vstold, mst(:,g), vst(:,:,g), x(:,g), nodes, ndps, w, mthd);
        
         % step E and F: compute the observation message and new belief
         muot_old = muot(:,g);
         hot_old = hot(:,:,g);
         while 1
            [muot(:,g),hot(:,:,g),mst(:,g), vst(:,:,g)]...
                =upd_mssg_belief( mstold, ...
            vstold, mstn, vstn, muot_old, hot_old,lr(g));   
            if min(eig(vst(:,:,g))) > 0, break; end
            lr(g) = lr(g)/2;
        end
        Lt(:,:,g) = inv(A*vst(:,:,g)*A' + tao);
     else
       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);
       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
   
   if 0
     for i =1:ng  
       vtmp(i) = 1/hot(:,:,i);
       mtmp(i) = vtmp(i)*muot(:,i);
       vh(i) = hot(:,:,i);
       vs(i) = vst(:,:,i);
    end
    figure(10),
    subplot(411), plot(vtmp),
    subplot(412),plot(mtmp),
    subplot(413), plot(muot),
    subplot(414),plot(vh),
   end
   
    if k == niter,
        break; 
    end
        
    % lr = linspace(ng,ng/7,ng)/ng/3;  lr = lr.^.5;
    lr = .3*ones(ng,1);    
    mdfall = 0;
    % step G and H: kalman smoothing and obs. message updating
    for g = ng-1:-1:1  
       if mod(g,200)== 0, fprintf('\n k:%d  i:%d',k,g); 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');
         
       % end
       % if 0
        %  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(g) 'th dt point'];
           warning(wstr);
           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, ndps, w, mthd);
     
       % Update the obs message and the new belief       
        muot_old = muot(:,g);
        hot_old = hot(:,:,g);
        while 1
            [muot(:,g), hot(:,:,g), mst(:,g),vst(:,:,g)]= ...
                upd_mssg_belief( mstold, vstold, ...
                mstn, vstn,muot_old,hot_old,lr(g));
            if min(eig(vst(:,:,g))) > 0, break; end    
            lr(g) = lr(g)/2;
        end
   end % end of F and G steps: smoothing and message updating
   
   if nfig == 1    
       figure(1),
       plot(t(1:end),exp(mst(1:end)*scl),lnty(k));
       hold on;
   end        
%   if mdfall/n < 1e-5, 
%        fprintf('\nep converge at iteration %d', k); 
%       break; 
%   end
end    

    mdfall = 0;
    % step G and H: kalman smoothing and obs. message updating
    for g = ng-1:-1:1  
       if mod(g,50)== 0, fprintf('\n k:%d  i:%d',k,g); 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');
         
   end
   
   if nfig == 1    
       figure(1),
       plot(t(1:end),exp(mst(1:end)*scl),'k');
       hold on;
   end