function [J, e] = x0estTimeErrorJacForDTModel(A, B, C, D, yu, x0sim, ...
   FreeIndex, maxsize, sqrlam)
% Error/Jacobain calculation for initial state estimation in time domain.
% Return J, e for solving J*x0 = e. Note: e is not output error!
%
% Inputs:
%  A, B, C, D: state-space matrices
%  yu: [y, u]: matrix of output/input data for one experiment.
%  x0sim: initial state structure created using linsimstate. Must contain
%         relevant input lag info. For passing in negative delay
%         information, add a "InputDelay" field to x0sim returned by
%         linsimstate.
%  FreeIndex: free elements of x0; logical vector.
%  maxsize: options.Advanced.MaxSize: size of largest matrix allowed to be
%           formed (default: 250e3)
%  sqrlam: sqrtm(OutputWeight)

% Copyright 1986-2022 The MathWorks, Inc.

nx = size(A,1); [ny, nu] = size(D); 
Ncap = size(yu,1); Ncap2 = Ncap;
nxf = sum(FreeIndex);
M = max(nxf+1,floor(maxsize/(nxf+nu+ny)));
HasNegativeDelay = isfield(x0sim,'InputDelay') && any(x0sim.InputDelay<0);
if HasNegativeDelay
   id = x0sim.InputDelay;
   simsys = ltipack.ssdata(A,B,C,D,[],1);
   simsys.Delay.Input = id;
   x0 = x0sim.State; x0_ = x0;
   Ncap2 = Ncap-max(-id);
   if Ncap2<=0
      ctrlMsgUtils.error('Ident:estimation:x0estNegDelStructInit')
   end
end

% Simulate using initial guess of x0; that way J\e gives delta_x0.
if ny>1 || M<Ncap
   if HasNegativeDelay
      psisys = ltipack.ssdata(A,zeros(nx,1),C,zeros(ny,1),[],1);
   end
   R1 = [];
   for kc = 1:M:Ncap
      jj = (kc:min(Ncap,kc-1+M));
      if jj(length(jj))<Ncap
         jjz = [jj,jj(length(jj))+1];
      else
         jjz = jj;
      end
      
      psi = zeros(ny*length(jj),nxf);
      if ~HasNegativeDelay
         [yh,x0sim] = idsssim(A,B,C,D,yu(jjz,ny+1:end).',x0sim);
         yh = yh.';
      else
         t = (0:1:(numel(jjz)-1))';
         [yh,xh] = idpack.ltidata.idlsim(...
            simsys,yu(jjz,ny+1:end),t,x0_,'zoh');
         x0_ = xh(end,:).';
      end
      
      e = (yu(jj,1:ny)-yh(1:length(jj),:))*sqrlam;
      evec = e(:);
      
      kl = 1;
      x0dum = x0sim;
      for kx = 1:nx
         if ~FreeIndex(kx), continue; end
         if ~HasNegativeDelay
            if kc==1
               x0dum.State = zeros(nx,1); x0dum.State(kx) = 1;
               x0dum.Input = {zeros(0,1)};
            else
               x0dum.State = X00(:,kl); x0dum.Input = {zeros(0,1)};
            end
            [psitemp, xf] = idsssim(A,zeros(nx,1),C,zeros(ny,1),zeros(1,length(jjz)),x0dum);
            psitemp = psitemp.';
            X00(:,kl) = xf.State; %#ok<AGROW>
         else
            % gradient calculation for negative delay case
            if kc==1
               x0psi_ = zeros(nx,1); x0psi_(kx) = 1;
            else
               x0psi_ = xh(end,:).';
            end
            t = (0:1:length(jjz)-1)';
            [psitemp, xh] = idpack.ltidata.idlsim(psisys,zeros(length(jjz),1),t,x0psi_,'zoh');
         end
         psitemp = psitemp(1:length(jj),:)*sqrlam;
         psi(:,kl) = psitemp(:); kl = kl+1;
      end
      if ~isempty(R1)
         if size(R1,1)<nxf+1
            ctrlMsgUtils.error('Ident:estimation:x0estSmallMaxSize')
         end
         R1 = R1(1:nxf+1,:);
      end
      H1 = [R1;[psi,evec] ];
      if HasNegativeDelay
         % cut out the last max(-id) samples because those were obtained
         % assuming the input was zero.
         H1 = H1(1:min(length(jj),Ncap2),:);
      end
      R1 = triu(qr(H1));
   end
   
   [nRr,nRc] = size(R1);
   R1 = R1(1:min(nRr, nRc),:);
   J = R1(:,1:end-1); e = R1(:,end);
else
   % Direct estimation without segmentation
   % RE: State trajectory returned by IDSSSIM and TIMERESP is "one step
   % behind" relative to value obtained as unit impulse response of
   % sensitivity system   
   if ~HasNegativeDelay
      y0 = idsssim(A,B,C,D,yu(:,2:end).',x0sim).';
      x0sim.Input = {zeros(0,1)};
      [~,~,psix0] = idsssim(A.',C.',zeros(1,nx),zeros(1,ny),[1,zeros(1,Ncap)],x0sim);
      psix0 = psix0(:,2:end).';
   else
      syspsi = ltipack.ssdata(A',C',zeros(0,nx),zeros(0,ny),[],1);
      t = (0:1:size(yu,1))';
      y0 = idpack.ltidata.idlsim(simsys,yu(:,2:end),t(1:end-1),x0,'zoh');
      [~,~,~,psix0] = timeresp(syspsi,'impulse',t,RespConfig(InitialState=x0*0));
      psix0 = psix0(2:end,:);
   end
   J = psix0(:,FreeIndex)*sqrlam; e = (yu(:,1)-y0)*sqrlam;
   if HasNegativeDelay
      % Cut out the last max(-id) samples because those were obtained
      % assuming the input was zero.
      J = J(1:Ncap2,:); e = e(1:Ncap2);
   end
end
