function [smll,out] = ekfit(par,i,y,m0,P0,skip,numNoise,funa,funh)
% Minus Log Liklihood of measurments given parameters fit0 data y stim i
% function smll = hhfilter(par,i,y,m0,P0,skip)
% 
% Model: x_k = funa(x_{k-1}, q_{k-1});  q_{k-1} ~ N(0,Q);    
%        y_k = funh(x_k,r_k);           r_k     ~ N(0,R); 

% Simplest noise model
Q = eye(numNoise(1));
R = eye(numNoise(2));

% Zero vectors for differentiating wrt the other operand 
zeroq = zeros(1,numNoise(1));
zeror = zeros(1,numNoise(2));

% Zero vectors for differentiating wrt the 0 operand (AD) 
zeroq_AD = myAD(zeroq);
zeror_AD = myAD(zeror);

% Initialize the minus-log-likelihood (mll)
mll = zeros(1,size(i,1));

% Main Filtering Loop
for k = 2:size(y,1)

  % Evaluate derivatives for prediction
  result = funa(m0,i(k,:),zeroq_AD,par);
  Wm = getderivs(result); % Wm = d/dq a(m0,i_k,0_q,p)
  m0_AD = myAD(m0);
  result = funa(m0_AD,i(k,:),zeroq,par);
  Am = getderivs(result); % Am = d/dx a(m0,i_k,0_q,p)

  % Prediction
  mb = getvalue(result);  % mb = a(m0,i_k,0_q,p)
  Pb = Wm*Q*Wm' + Am*P0*Am';

  % Evaluate derivatives for update
  mb_AD = myAD(mb);
  result = funh(mb,i(k,:),zeror_AD,par);
  V = getderivs(result);  % V = d/dr h(mb,0_r);
  result = funh(mb_AD,i(k,:),zeror,par);
  H = getderivs(result);  % H = d/dx h(mb,0_r);

  % Update
  h = getvalue(result);
  v = y(k) - h;
  S = H*Pb*H' + V*R*V';
  K = Pb*H'/S;
  m0 = mb + K*v;
  P0 = Pb - K*S*K';

  % Marginal Measurement Likelihood
  mll(k) = log(det(S)) + v'/S*v;
end

smll = sum(mll(skip:end));

