function dydx = utEvalStateJacobian(x,par,type)
%UTEVALSTATEJACOBIAN Gateway function for jacobian computation (dy/dx).
%
%  dydx = utEvalStateJacobian(x,par,type)
%   Output:
%       dydx: derivative of output w.r.t regressors (x).
%   Inputs:
%       x: regressors
%       par: parameter struct (may contain additional info, such as NumberOfUnits)
%       type: class name of nonlinearity estimator for which jacobian is to
%       be computed.
%
%  This function is called from S function blocks during linearization, as
%  a mexCallMATLAB call.
%
%  See also idnlfun/setpar, idnlfun/getJacobian, utEvalCustomReg.

%   Written by: Rajiv Singh
%   Copyright 2007-2023 The MathWorks, Inc.

SetParFlag = true;
nu = length(x);

try
   switch lower(type)
      case {'idtreepartition','treepartition'}
         nlobj = idTreePartition;
         rI = par.RelativeNonlinearIndex;
         rJ = par.RelativeLinearIndex;
         In = "u"+(1:nu);
         nlobj = initializeMultipleIOs(nlobj,In,"y",In(rJ),"y",In(rI),"y");
         nlobj.NonlinearFcn.NumberOfUnits = par.NumberOfUnits;
      case {'idlinear','linear'}
         nlobj = idLinear;
         q = size(par.LinearSubspace,2);
         In = "u"+(1:nu);
         nlobj = initializeIOs(nlobj,In,"y");
         nlobj = initializeProjectionDim(nlobj,q);
      case {'idwaveletnetwork','wavenet'}
         nlobj = idWaveletNetwork;
         lq = size(par.LinearSubspace,2);
         nlq = size(par.NonlinearSubspace,2);
         ns = length(par.ScalingCoefficient);
         nw = length(par.WaveletCoefficient);
         I1 = par.RelativeLinearIndex;
         I2 = par.RelativeNonlinearIndex;
         In = "u"+(1:nu);
         nlobj = initializeMultipleIOs(nlobj,In,"y",In(I1),"y",In(I2),"y");
         nlobj = initializeMultipleProjectionDim(nlobj,lq,nlq);
         nlobj.NonlinearFcn = initializeNumBasis(nlobj.NonlinearFcn, ns, nw);
      case {'idsigmoidnetwork','sigmoidnet','idgaussianprocess'}
         IsGP = strcmpi(type,'idgaussianprocess');
         if IsGP
            nlobj = idGaussianProcess;
         else
            nlobj = idSigmoidNetwork;
         end
         lq = size(par.LinearSubspace,2);
         nlq = size(par.NonlinearSubspace,2);
         I1 = par.RelativeLinearIndex;
         I2 = par.RelativeNonlinearIndex;
         In = "u"+(1:nu);
         nlobj = initializeMultipleIOs(nlobj,In,"y",In(I1),"y",In(I2),"y");
         nlobj = initializeMultipleProjectionDim(nlobj,lq,nlq);
         if ~IsGP
            nlobj.NonlinearFcn.NumberOfUnits = par.NumberOfUnits;
         end
      case {'idpiecewiselinear','pwlinear'}
         nlobj = idPiecewiseLinear;
      case {'idfeedforwardnetwork','idcustomnetwork','neuralnet','customnet',...
            'idtreeensemble','idsupportvectormachine'}
         nlobj = par; % par itself is the object
         SetParFlag = false;
      case {'idpolynomial1d','poly1d'}
         nlobj = idPolynomial1D;
      case {'idsaturation','saturation'}
         nlobj = idSaturation;
      case {'iddeadzone','deadzone'}
         nlobj = idDeadZone;
      case {'idunitgain','unitgain'}
         nlobj = idUnitGain;
      otherwise
         error(message('Ident:utility:JacobianCalculationFailure1',upper(type)))
   end
   
   if SetParFlag
      nlobj = setpar(nlobj,par);
   end
   [~,~,dydx] = getJacobian(nlobj,x,false);
catch
   warning(message('Ident:utility:JacobianCalculationFailure2'))
   dydx = zeros(size(x));
end
