classdef ModalTruncation < mor.GenericSpec
   % Specification for reducing LTI model using Modal Truncation.
   %
   %   For non-sparse LTI models SYS, R = REDUCESPEC(SYS,'modal') creates a 
   %   mor.ModalTruncation object R for reducing SYS with the Modal Truncation
   %   method. This computes a modal decomposition (see MODSEP) and eliminates
   %   non-essential modes based on various criteria. Use R.Options to further
   %   configure this spec, type "help(R.Options)" for details.
   %
   %   R = PROCESS(R) processes the specification R and collects the
   %   following information in R:
   %     Mode       Average pole location in modal components.
   %     Damping    Damping ratio of average poles.
   %     Frequency  Natural frequency of average poles.
   %     Component  Modal components.
   %     DCGain     Normalized DC contribution of modal components (with
   %                optional input and output scaling).
   %     TL,TR      Left and right block-diagonalizing transformations 
   %                (TL'*A*TR is block diagonal).
   %
   %   Use VIEW to decide which modes to keep with one of these plots:
   %     'mode'     Mode locations.
   %     'damp'     Mode dampings and natural frequencies.
   %     'contrib'  DC contributions of modal components.
   %
   %   Use GETROM to obtain the modal form of SYS or a reduced-order model
   %   with one of these mode selection criteria:
   %     Selection  Selected modes (logical or index vector into R.Modes).
   %     Frequency  Frequency range [FMIN,FMAX] of interest. Modes outside
   %                this range are discarded.
   %     Damping    Damping range [ZMIN,ZMAX] of interest. Modes outside this
   %                range are discarded.
   %     MinDC      Modes with smaller normalized DC contribution are
   %                discarded.
   %     Method     Mode elimination method ('matchDC' or 'truncate').
   %                The "matchDC" method (default) preserves the DC gain
   %                (steady-state value of step response) while the "truncate"
   %                method preserves the feedthrough term.
   %
   %   Example: Keep only the modes with natural frequency in [0,2] and a DC
   %   contribution of at least 10%:
   %      rng(3), sys = rss(30);
   %      R = reducespec(sys,'modal');
   %      % Run algorithm once for efficiency (optional)
   %      R = process(R)
   %      % View relative DC contribution of each mode
   %      view(R,'contrib')
   %      % Get reduced-order model
   %      rsys = getrom(R,Frequency=[0 10],MinDC=0.1)
   %      figure(1), bode(sys,rsys)
   %      figure(2), step(sys,rsys)
   %
   %   See also REDUCESPEC, MOR.PROCESS, MOR.VIEW, MOR.GETROM, MODSEP.

   % Copyright 2020-2023 The MathWorks, Inc.
   properties (SetAccess=protected, Dependent)
      Mode
      Damping
      Frequency
      Component
      DCGain
      TL
      TR
   end

   methods

      function R = ModalTruncation(sys)
         % Convert to @ss
         if nargin>0  % for call by subclass
            try
               R.Model_ = ss(sys);
            catch
               error(message('Control:transformation:ROM1'))
            end
            R.Options_ = mor.ModalTruncationOptions;
         end
      end

      function value = get.Mode(this)
         if isempty(this.Cache_)
            value = [];
         else
            value = this.Cache_.Mode;
         end
      end

      function value = get.Damping(this)
         if isempty(this.Cache_)
            value = [];
         else
            value = this.Cache_.Damping;
         end
      end

      function value = get.Frequency(this)
         if isempty(this.Cache_)
            value = [];
         else
            value = this.Cache_.Frequency;
         end
      end

      function value = get.Component(this)
         if isempty(this.Cache_) || this.Options_.ModeOnly
            value = [];
         else
            value = ss.make(this.Cache_.Component);
         end
      end

      function value = get.DCGain(this)
         if isempty(this.Cache_) || this.Options_.ModeOnly
            value = [];
         else
            value = this.Cache_.DCGain;
         end
      end

      function value = get.TL(this)
         if isempty(this.Cache_) || this.Options_.ModeOnly
            value = [];
         else
            value = this.Cache_.TL;
         end
      end

      function value = get.TR(this)
         if isempty(this.Cache_) || this.Options_.ModeOnly
            value = [];
         else
            value = this.Cache_.TR;
         end
      end

      function R = process(R)
         % Recompute derived quantities
         if isempty(R.Cache_)
            try
               R.Cache_ = runMT(R.Model_,R.Options_);
            catch ME
               throw(ME)
            end
         end
      end

      function [rsys,varargout] = getrom(R,arg2,optArgs)
         %GETROM  Obtain reduced-order model (ROM).
         %
         %   RSYS = GETROM(R) returns the modal form of the original system.
         %
         %   RSYS = GETROM(R,Name1=Value1,Name2=Value2,...) returns a reduced-order
         %   model RSYS based on the specified mode selection criteria. Available
         %   criteria and options include:
         %     Selection  Selected modes (logical or index vector into R.Modes).
         %     Frequency  Frequency range [FMIN,FMAX] of interest. Modes outside
         %                this range are discarded.
         %     Damping    Damping range [ZMIN,ZMAX] of interest. Modes outside this
         %                range are discarded.
         %     MinDC      Modes with smaller normalized DC contribution are
         %                discarded.
         %     Method     Mode elimination method ('matchDC' or 'truncate').
         %                The "matchDC" method (default) preserves the DC gain
         %                (steady-state value of step response) while the "truncate"
         %                method preserves the feedthrough term.
         %
         %   [RSYS,INFO] = GETROM(R,...) also returns a struct INFO containing the
         %   left and right projectors PL and PR that map the original matrices
         %   (A,B,C,D) to the reduced-order matrices (PL'*A*PR,PL'*B,C*PR,D).
         %
         %   See also REDUCESPEC.
         arguments
            R
            arg2 = '';
            optArgs.?mor.getrom.ModalTruncationOptions;
         end
         %narginchk(2,Inf)

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

         % Update cache
         try
            R = process(R);
         catch ME
            throw(ME)
         end
         Cache = R.Cache_;

         % Resolve selection
         ROMopt = initOptions(mor.getrom.ModalTruncationOptions,namedargs2cell(optArgs));
         select = ROMopt.Selection(:);
         nModes = numel(Cache.Mode);
         if isempty(select)
            select = (Cache.Frequency>=ROMopt.Frequency(1) & Cache.Frequency<=ROMopt.Frequency(2) & ...
               Cache.Damping>=ROMopt.Damping(1) & Cache.Damping<=ROMopt.Damping(2));
            if ~R.Options_.ModeOnly
               select = (select & Cache.DCGain>=ROMopt.MinDC);
            end
         elseif isnumeric(select)
            % Index vector
            if any(select>nModes)
               error(message('Control:transformation:MODALROM19',nModes))
            end
            select = ismember((1:nModes)',select);
         else
            % Logical vector
            if numel(select)~=nModes
               error(message('Control:transformation:MODALROM2',nModes))
            end
         end

         % Get reduced model
         MatchDC = strcmp(ROMopt.Method,'matchDC');
         try
            [rsys,varargout{1:nargout-1}] = reduceMT(R.Model_,select,Cache,MatchDC);
         catch ME
            throw(ME)
         end
         % Clear notes, userdata, etc
         rsys.Name = '';  rsys.Notes = {};  rsys.UserData = [];
      end
      

      function varargout = view(R,plotType,optArgs)
         %VIEW  Graphically select order of reduced model.
         %
         %   VIEW(R,PLOTTYPE) shows the modal content of the original model
         %   and provides insight into which modes can be eliminated to reduce 
         %   the model order. Available plot types include:
         %     'mode'     Mode locations.
         %     'damp'     Mode dampings and natural frequencies.
         %     'contrib'  Bar chart of normalized DC contributions.
         %
         %   VIEW(R) is equivalent to VIEW(R,'mode').
         %
         %   VIEW(R,...,Axes=AX) plots into the axes with handle AX.
         %
         %   H = VIEW(R,...) returns the plot handle.
         %
         %   See also REDUCESPEC.
         arguments
            R
            plotType = 'mode';
            optArgs.?mor.view.ModalTruncationOptions;
         end
         nargoutchk(0,1)
         % Handle help request view(R,'-help')
         if strcmp(plotType,'-help')
            if nargout>0
               error(message('Control:transformation:view1'))
            end
            help mor.ModalTruncation.view
            return
         end
         % Resolve plot type
         plotType = ltipack.matchKey(plotType,{'mode','damp','contribution'});
         if isempty(plotType)
            error(message('Control:transformation:MODALROM3'))
         end
         % Update data
         try
            R = process(R);
         catch ME
            throw(ME)
         end
         % Create plot
         if isfield(optArgs,'Axes') && isgraphics(optArgs.Axes,'axes')
            ax = optArgs.Axes;
         else
            ax = gca;
         end
         switch plotType
            case 'mode'
               % Mode locations
               p = R.Mode;
               if R.Cache_.REAL && ~R.Options_.ModeOnly
                  p = [p ; conj(p(imag(p)>0,:))];
               end
               plot(ax,real(p),imag(p),'x',LineWidth=1);
               title(ax,getString(message('Control:transformation:MODALROM4')))
               xlabel(ax,getString(message('Controllib:plots:strRealAxis')))
               ylabel(ax,getString(message('Controllib:plots:strImaginaryAxis')))
               grid(ax,'on')
            case 'damp'
               % Damping vs natural frequency
               semilogy(ax,R.Damping,R.Frequency,'x',LineWidth=1)
               xlim = get(ax,'XLim');
               set(ax,'XLim',max(-1,min(1,xlim)))
               title(ax,getString(message('Control:transformation:MODALROM5')))
               xlabel(ax,getString(message('Controllib:plots:strDamping')))
               ylabel(ax,getString(message('Control:transformation:MODALROM6')))
               grid(ax,'on')
            case 'contribution'
               if R.Options_.ModeOnly
                  error(message('Control:transformation:MODALROM7'))
               end
               bar(ax,sort(R.DCGain,'descend'));
               set(ax,'XLim',[0.5 numel(R.Mode)+0.5],'Yscale','log')
               title(ax,getString(message('Control:transformation:MODALROM8')))
               xlabel(ax,getString(message('Control:transformation:MODALROM9')))
               ylabel(ax,getString(message('Control:transformation:MODALROM10')))
               grid(ax,'on')
         end
         % Return handle if requested
         if nargout>0
            varargout = {ax};
         end
      end

   end

end

