classdef LinearODE < pmodel.ODE
   % Linear ODE model structure.
   % This represents a state space model with customized parameterization.
   
   %   Author(s): Rajiv Singh
   %   Copyright 2010-2017 The MathWorks, Inc.
   
   properties (Dependent)
      % The sample time handling behavior of the ODE function.
      % 
      % The "FunctionType" property specifies how the sample time affects the
      % state-space output arguments returned by the ODE function. It takes
      % one the following values:
      %   'c': Denotes that the ODE function returns data corresponding to
      %        a continuous-time state-space description, regardless of the
      %        value of the sample time (Ts). Default. 
      %   'd': Denotes that the ODE function returns data corresponding to
      %        a discrete-time state-space description, regardless of the
      %        value of the sample time (Ts). The output arguments should
      %        not depend upon the value of Ts.
      %  'cd': Denotes that the ODE function returns continuous-time
      %        state-space matrices when the sample time is zero and
      %        discrete-time matrices when the sample time is non-zero.
      %        The continuous- and discrete-time matrix sizes must match.
      %        Note that the sample time (Ts) is an input argument to the
      %        ODE function.
      FunctionType
   end
   
   properties (Hidden)
      FcnType % Alias for "FunctionType" property
   end   
      
   methods (Static, Hidden)
      function D = array(size)
         % Create an array of a given size
         D = pmodel.LinearODE.newarray(size);
      end
   end
   
   methods (Access = protected)
      
      % PARAMETER SERIALIZATION INTERFACE (PMODEL.GENERIC)
      function ps = getParamSet(M)
         ps = M.Parameters;
      end
      
      function M = setParamSet(M,ps)
         M.Parameters = ps;
      end
      
      function dispScalarModel(M)
         % Display of a scalar linear ODE object.
         Fcn0 = M.Function;
         if isempty(Fcn0)
            disp(getString(message('Control:pmodel:msgEmptyODE')))
            return
         elseif ischar(Fcn0) || isstring(Fcn0)
            [~, Fcn] = fileparts(Fcn0);
         else
            Fcn = func2str(Fcn0);
         end
         FcnType = M.FcnType;
         S = struct('Function',Fcn0,'FunctionType',FcnType);
         S.Parameters = M.Parameters;
         S.ExtraArguments = M.ExtraArgs;
         S.StateName = M.StateName;
         S.StateUnit = M.StateUnit;
         disp(S)
         
         disp(ctrlMsgUtils.message('Control:pmodel:msgLinearODEDisp1',Fcn,numel(M.Parameters)))         
         switch FcnType
            case 'c'
               disp(ctrlMsgUtils.message('Control:pmodel:msgLinearODEDisp2',Fcn))
            case 'd'
               disp(ctrlMsgUtils.message('Control:pmodel:msgLinearODEDisp3',Fcn))
            case 'cd'
               disp(ctrlMsgUtils.message('Control:pmodel:msgLinearODEDisp4',Fcn))
         end
         disp(' ')
         
      end
      
      function dispModelType(~,Size)
         % Display model type information.
         if prod(Size)~=1
            disp(ctrlMsgUtils.message('Control:pmodel:msgLinearODEDispArray',mat2str(Size)))
         end
      end
   end
   
   methods
      function M = LinearODE(varargin)
         % PMODEL.LinearODE Constructor.
         % Syntax:
         %   M = LinearODE(Orders, ODEFun, FcnType, Parameters, ExtraArgs)
         %
         % Orders: [ny,nu,nx]
         % ODEfun: name of ODE file or function handle
         % FcnType: 'c', 'd', or 'cd'
         % Parameters: Cell array of parameters, each element a numeric
         %             array
         % ExtraArgs: Cell array of arbitrary values to be passed to the
         %            ODE function as extra input arguments. Optional.
         
         % Currently, number of states nx must be fixed.
         ni = nargin;         
         if ni>0
            FcnType = varargin{3};
            varargin = varargin([1 2 4 5]);
         else
            FcnType = 'c';
         end
         M = M@pmodel.ODE(varargin{:});
         M.FunctionType = FcnType;
      end
      
      function Value = get.FunctionType(M)
         %GET method for "FunctionType" property
         Value = M.FcnType;
      end
      
      function M = set.FunctionType(M, Value)
         %SET method for "FunctionType" property
         M.FcnType = char(Value);
      end
      
      function [A,B,C,D,K,X0,varargout] = getABCDKX0(M,Ts,Nargout)
         % Get state space data.
         % Syntax:
         %     [A,B,C,D,K,X0,varargout] = getABCDKX0(M,Ts);
         Fcn = getFunctionHandle(M);
         Par = getParValue(M);
         Extras = M.ExtraArgs;
         Input = [{Fcn},Par,{Ts},Extras];
         if nargout<5
            % ODE function must return at least four output arguments.
            [A,B,C,D] = feval(Input{:});
         else
            if nargin<3 || isempty(Nargout)
               [~,A,B,C,D,K,X0] = getFcnNargout(M, Ts);
            elseif Nargout==4
               [A,B,C,D] = feval(Input{:});
               K = NaN(size(C,2),size(C,1));
               X0 = NaN(size(A,1),1);
            elseif Nargout==5
               [A,B,C,D,K] = feval(Input{:});
               X0 = NaN(size(A,1),1);
            elseif Nargout>=6
               [A,B,C,D,K,X0,varargout{1:nargout-6}] = feval(Input{:});
            else
               % Nargout is unknown (-1)
               [~,A,B,C,D,K,X0] = getFcnNargout(M, Ts);
            end
         end
      end
      
      function [N,A,B,C,D,K,X0] = getFcnNargout(M, Ts)
         % Determine number of output arguments of ODE function.
         Fcn = getFunctionHandle(M);
         Par = getParValue(M);
         Extras = M.ExtraArgs;
         Input = [{Fcn},Par,{Ts},Extras];
         
         try
            N = nargout(Fcn);
            if N>=4
               [varg{1:N}] = feval(Input{:});
               [A,B,C,D] = deal(varg{1:4});
               sz = size(C);
               K = nan(sz(2),sz(1));
               X0 = nan(sz(2),1);
               if N>4
                  K = varg{5};
                  if N>5
                     X0 = varg{6};
                  end
               end
               return
            elseif N>=0
               error(message('Control:pmodel:linearODEChk2'));
            end
         catch E
            throw(E)
         end
         
         % this code applies only if N==-1
         try
            [A,B,C,D,K,X0] = feval(Input{:});
            N = 6;
         catch E
            expids = {'MATLAB:TooManyOutputs','MATLAB:unassignedOutputs'};
            if any(strcmp(E.identifier,expids))
               try
                  [A,B,C,D,K] = feval(Input{:});
                  N = 5;
                  X0 = nan(size(A,1),1);
               catch E2
                  if any(strcmp(E2.identifier,expids))
                     [A,B,C,D] = feval(Input{:});
                     N = 4; 
                     K = nan(size(C,2),size(C,1));
                     X0 = nan(size(A,1),1);
                  else
                     throw(E2)
                  end
               end
            else
               throw(E)
            end
         end
      end      
      
   end
   
   methods (Hidden)
      function M = checkDataType(M, Ts)
         % Validate model properties data type.
         M = checkDataType@pmodel.ODE(M);
         F = M.FcnType;
         if ~any(strcmpi(F,{'c','d','cd'}))
            ctrlMsgUtils.error('Control:pmodel:linearODEChk1')
         else
            M.FcnType = lower(F);
         end
         
         Nargout = getFcnNargout(M, Ts);
         if Nargout>0 && Nargout<4
            ctrlMsgUtils.error('Control:pmodel:linearODEChk2')
         end
         
         if (Ts==0 && strcmp(M.FcnType,'d')) || (Ts~=0 && strcmp(M.FcnType,'c'))
            ctrlMsgUtils.error('Control:pmodel:idgreyFcnTypeTsCompat')
         end
         
         M = checkOutputArgConsistency(M, Ts, Nargout);
      end
      
      function M = checkOutputArgConsistency(M, Ts, varargin)
         % Check that output arguments are dimensionally consistent.
         % varargin for "Nargout" specification 
         try
            [A,B,C,D,K,X0] = getABCDKX0(M,Ts,varargin{:});
         catch E
            Fcn = M.Function;
            if isa(Fcn,'function_handle'), Fcn = func2str(Fcn); end
            ctrlMsgUtils.error('Control:pmodel:odeChk3',Fcn,E.message)
         end
         N = M.Orders;
         ny = N(1); nu = N(2); nx = N(3);
         % Support for multiexperiment X0 is disabled here.
         if ~isequal(size(A),[nx nx]) || ~isequal(size(B), [nx, nu]) ||...
               ~isequal(size(C),[ny nx]) || ~isequal(size(D),[ny nu]) ||...
               ~isequal(size(K),[nx ny]) || ~isequal(size(X0),[nx 1])
            ctrlMsgUtils.error('Control:pmodel:linearODEChk3')
         end
      end
      
   end   
   
end
