classdef NCFBalancedTruncation < mor.GenericBTSpec & matlab.mixin.CustomDisplay
   % Specification for reducing LTI model using Balanced Truncation of
   % Normalized Coprime Factors.
   %
   %   R = REDUCESPEC(SYS,'ncf') creates a mor.NCFBalancedTruncation spec R
   %   for reducing the LTI model SYS by applying Balanced Truncation to 
   %   the left normalized coprime factors [Ml,Nl] of SYS (see LNCF).
   %
   %   R = PROCESS(R) processes the specification R and collects the
   %   following information in R:
   %     Sigma     Hankel singular values of [Ml,Nl].
   %     Energy    Normalized energy of each state, measured as
   %               (Sigma(k)/Sigma(1))^2.
   %     Error     Bound on absolute approximation error in the frequency
   %               domain. Error(k) is the error when setting the order to
   %               r=k-1.
   %     LNCF      State-space realization (A,B,C,D) of [Ml Nl].
   %     TL,TR     Left/right transformations.
   %     Lr,Lo     Cholesky factors of controllability and observability 
   %               gramians for (TL'*A*TR,TL'*B,C*TR).
   %
   %   The VIEW method offers the following plots to help decide which states
   %   to keep:
   %     'sigma'    Bar chart of Hankel singular values and associated error
   %                bound.
   %     'energy'   Bar chart of normalized energies.
   %   You can customize these plots using Name/Value pairs, for example,
   %      view(R,'sigma',Axes=gca,Yscale='linear')
   %   See HSVOPTIONS for available options besides 'Axes'.
   %
   %   The GETROM method lets you pick the order with one of the following
   %   criteria:
   %     Order      Desired order (scalar or vector for batch selection).
   %     MaxError   Maximum approximation error. This picks the lowest order
   %                for which the error stays below this value.
   %     MinEnergy  States with lower energy are discarded.
   %
   %   Example: Reduce a 30th-order state-space model using balanced 
   %   truncation of normalized coprime factors.
   %      rng(0), G = rss(30,2,3);
   %      R = reducespec(G,'ncf');
   %      % View Hankel singular values
   %      view(R)
   %      % Get lowest-order model with approximation error below 1%
   %      rsys = getrom(R,MaxError=0.01)
   %      sigma(G,G-Gr)
   %
   %   See also REDUCESPEC, MOR.PROCESS, MOR.VIEW, MOR.GETROM, LNCF.

   % Copyright 2020-2023 The MathWorks, Inc.
   properties (SetAccess=protected, Dependent)
      % [Ml,Nl] (model being reduced)
      LNCF
      TL
      TR
   end

   properties (Access=protected)
      LNCF_  % LNCF [Ml,Nl]
   end

   methods

      function R = NCFBalancedTruncation(sys)
         % Convert to @ss
         if hasdelay(sys)
            error(message('Control:transformation:ROM9'))
         end
         try
            R.Model_ = ss(sys);
         catch
            error(message('Control:transformation:ROM1'))
         end
      end

      function value = get.LNCF(R)
         if isempty(R.Cache_)
            value = [];
         else
            value = R.LNCF_;
         end
      end      

      function value = get.TL(R)
         % Note: Cache uses TL*A*TR but BalancedTruncation uses TL'*A*TR
         % for consistency with ModalTruncation.
         if isempty(R.Cache_)
            value = [];
         else
            value = R.Cache_.TL';
         end
      end

      function value = get.TR(R)
         if isempty(R.Cache_)
            value = [];
         else
            value = R.Cache_.TR;
         end
      end

      function R = process(R)
         % Note: No advantage in calling it separately (no options)
         if isempty(R.Cache_)
            try
               % Note: GL is stable and all-pass (GL*GL~ = I)
               GL = lncf(R.Model_);
               R.LNCF_ = GL;
            catch ME
               error(message('Robust:transformation:ncfmr4',ME.message))
            end
            try
               R.Cache_ = runBT(GL,mor.BalancedTruncationOptions);
            catch ME
               throw(ME)
            end
         end
      end

      function [rsys,INFO] = getrom(R,arg2,optArgs)
         %GETROM  Obtain reduced-order model (ROM).
         %
         %   RSYS = GETROM(R,Name=Value) returns a reduced-order model RSYS
         %   based on one of the following order selection options:
         %     Order      Desired order (scalar or vector for batch selection).
         %     MaxError   Maximum approximation error. This picks the lowest order
         %                for which the error does not exceed this value.
         %     MinEnergy  This discards all states with lower normalized energy.
         %
         %   [RSYS,INFO] = GETROM(R,...) also returns a struct INFO containing
         %   the approximation RLNCF=[Mlr,Nlr] of the left normalized coprime 
         %   factorization [Ml,Nl] (see R.LNCF). The reduced-order model is 
         %      RSYS = Mlr \ Nlr.
         %
         %   See also REDUCESPEC.
         arguments
            R
            arg2 = '';
            optArgs.?mor.getrom.NCFBalancedTruncationOptions;
         end

         % Handle help request getrom(R,'-help')
         if nargin==2
            if nargout>0 || ~strcmp(arg2,'-help')
               error(message('Control:transformation:getrom2'))
            end
            help mor.NCFBalancedTruncation.getrom
            return
         end

         % Update cache, grab options, and resolve order
         try
            R = process(R);
            ROMopt = initOptions(mor.getrom.NCFBalancedTruncationOptions,namedargs2cell(optArgs));
         catch ME
            throw(ME)
         end
         orders = resolveOrder(ROMopt,R.Cache_);

         % Reduce LNCF
         sys = R.Model_;
         [ny,nu] = size(sys);
         GLr = reduceBT(R.LNCF_,orders,R.Cache_,false);
         try
            rsys = -imp2exp(GLr,1:ny,ny+1:ny+nu,'min');
         catch
            error(message('Robust:transformation:ncfmr7'))
         end
         INFO = struct('RLNCF',GLr);
         % Clear notes, userdata, etc
         rsys.InputName = sys.InputName;
         rsys.InputUnit = sys.InputUnit;
         rsys.InputGroup = sys.InputGroup;
         rsys.OutputName = sys.OutputName;
         rsys.OutputUnit = sys.OutputUnit;
         rsys.OutputGroup = sys.OutputGroup;
         rsys.Name = '';  rsys.Notes = {};  rsys.UserData = [];
      end

      % NCFMR support
      function Info = getNCFInfo(R)
         % Get data structure formerly returned by NCFMR 
         Info = rctutil.ncfmrInfo(R.Cache_);
         Info.GL = R.LNCF_;
      end

      function R = setNCFInfo(R,Info)
         % Info = rctutil.ncfmrInfo
         R.Cache_ = strip(Info);
         R.LNCF_ = Info.GL;
      end

   end

   methods (Access = protected)
      % Reorder properties
      function propgrp = getPropertyGroups(~)
         proplist = {'Sigma','Energy','Error','LNCF','TL','TR','Lr','Lo','Options'};
         propgrp = matlab.mixin.util.PropertyGroup(proplist);
      end
   end

end

