function h = compareplot(varargin)
%COMPAREPLOT Compares system response with the measured data.
%
%   COMPAREPLOT(DATA, SYS, M)
%     compares the response of the system SYS to the validation data DATA.
%     It plots the response of SYS superimposed over the response data in
%     DATA and displays the normalized root mean square (NRMSE) measure of
%     the goodness of fit.
%
%     DATA is validation data (an IDDATA or IDFRD object) for which the
%     comparison is made. For linear systems, DATA can be either in time-
%     or frequency-domain. For nonlinear systems (such as IDNLARX), DATA
%     must be in the time domain.
%
%     SYS is a dynamic system (such as IDTF, IDSS, FRD, TF, ZPK object) or
%     an IDDATA object.
%      * When SYS is an IDDATA object, DATA must be an IDDATA object too
%        with matching domains, number of experiments and time or frequency
%        vectors.
%      * When SYS is an FRD model (objects of class FRD or IDFRD), DATA
%        must be an FRD model too. If the frequency vectors in DATA and
%        SYS do not match, DATA is interpolated to the frequency values
%        in SYS (no extrapolation) and the fit is computed over the
%        overlapping frequency regions.
%      * When time and/or frequency units of DATA do not match those of SYS,
%        the system is rescaled to the units of DATA.
%
%     M is the prediction horizon. Measured output values in DATA up to
%     time t-M are used to predict the output of SYS at time t. All
%     relevant inputs are used. M = Inf gives a pure simulation of the
%     system. (Default: M = Inf). M is ignored when SYS is an IDDATA
%     object, an FRD model or a dynamic system with no noise component
%     (such as systems belonging to Control System Toolbox).
%
%     The fit is calculated as:
%         FIT = 100 * (1-norm(Y-YHAT)/norm(Y-mean(Y))) (in %)
%     where Y is the output of the validation data and YHAT is the test
%     system output. The matching of input/output channels in DATA and SYS
%     is based on the channel names. It is thus possible to evaluate models
%     that do not use all input channels available in DATA.
%
%   COMPAREPLOT(DATA, SYS, STYLE, M)
%     allows specification of line style, color or marker for the plot of
%     the response of model SYS. STYLE should be a string of applicable
%     characters, as available for the PLOT command. STYLE, if used, should
%     precede the specification of prediction horizon M.
%
%   COMPAREPLOT(DATA, SYS1, SYS2, ..., SYSn, M)
%     compares several models. You can also specify a color, line style,
%     and marker for each system. If both plot styles and prediction
%     horizon are specified, the prediction horizon must be specified as
%     the last argument (following the list of models and their plot
%     styles), as in:
%      COMPAREPLOT(DATA, sys1, 'r', sys2, 'y--', sys3, 'gx', M).
%
%   COMPAREPLOT(DATA, SYS1, SYS2, ..., SYSn, M, OPTIONS)
%     gives access to additional options for configuring the handling of
%     initial conditions, data offsets and other settings affecting the
%     computation of model response and fit values. See COMPAREOPTIONS for
%     more information.
%
%   COMPAREPLOT(AX, ...) plots into the axes with handle AX.
%
%   See also COMPARE, COMPAREOPTIONS, IDDATA/PLOT, DYNAMICSYSTEM/BODE,
%   IDLTI/SPECTRUM, PREDICT, IDPARAMETRIC/FORECAST.

%   Author(s): Rajiv Singh
%   Copyright 2011-2023 The MathWorks, Inc.

%#function idguis.internal.plots.compare
%#function iduis.utSpectrumData
%#function iduis.utMagPhaseIOData

ni = nargin;
narginchk(1,Inf)
ArgNames = cell(ni,1);
for ct = 1:ni
   ArgNames(ct) = {inputname(ct)};
end

% Parse input list
% Check for axes argument
if ishghandle(varargin{1},'axes')
   ax = varargin{1};
   varargin(1) = [];
   ArgNames(1) = [];
else
   ax = [];
end

% Parse input list
try
   [refDataStruct, sysList, Extras, PVPairs] = idmodel.parseCompareResidPredictInputs(varargin,ArgNames,'compare',true);
   [refData, sysList, Options, ynames, unames, Horizon] = ...
      idmodel.validateCompareInputs(refDataStruct, sysList,Extras,PVPairs,'compare');
catch E
   throw(E)
end

if isempty(refData.Name)
   refData.Name = getString(message('Ident:plots:msgValiDat'));
end

IsIddata = isa(refData,'iddata'); % to distinguish between IDDATA, FRD

TD = IsIddata && strcmpi(refData.Domain,'time');
% IsReal: determines whether real and imaginary parts have to be shown
% separately (in case curves are complex valued) or not.
IsReal = TD && isreal(refData);
PSmod = isa(refData,'FRDModel') && size(refData,2)==0;
nsys = numel(sysList);

if IsReal
   for ksys = 1:nsys
      IsReal = isreal(sysList(ksys).System) && localIsreal(sysList(ksys).IC);
      if ~IsReal, break; end
   end
end
src = repmat(handle(NaN),nsys,1);
for ct = 1:length(sysList)
   src(ct) = plotpack.idsimsource(sysList(ct).System,'Name',sysList(ct).Name);
   if ~IsIddata
      if PSmod
         yname = sysList(ct).System.OutputName;
         if ~isempty(yname) && ~all(strcmp(yname,''))
            cnames = cellfun(@(x)['e@',x],yname,'uni',false);
         else
            cnames = repmat({''},[length(yname),1]);
         end
      else
         cnames = sysList(ct).System.InputName;
      end
      uOrd = sysList(ct).DataInputOrder;
      I = ~isnan(uOrd);
      cnames = cnames(I);
      src(ct).ColumnNames = cnames;
   elseif IsIddata && TD && ~IsReal
      src(ct).ColumnNames = {'Real','Imag'};
   else
      src(ct).ColumnNames = {''};
   end
end

yname = refData.OutputName;
for ky = 1:numel(yname)
   if isempty(yname{ky}), yname{ky} = ['y',int2str(ky)]; end
end
refData.OutputName = yname;

if PSmod
   uname = cellfun(@(x)['e@',x],yname,'uni',false);
else
   uname = refData.InputName;
end

for ku = 1:numel(uname)
   if isempty(uname{ku}), uname{ku} = ['u',int2str(ku)]; end
end

if ~PSmod
   refData.InputName = uname;
end

% Create plot (visibility ='off')
if isempty(ax)
   ax = gca;
   % REVISIT: cla(ax,'reset')
end

if ~isempty(ancestor(ax,'matlab.graphics.layout.Layout'))
   %error(message('MATLAB:handle_graphics:Layout:UnsupportedChart'))
   clf
   ax = gca;
end

h = localCompare(ax, refData, src, IsReal, Options, ynames, unames, Horizon, sysList);
if nargout==0, clear h, end

%--------------------------------------------------------------------------
function h = localCompare(ax, refData, src, IsReal, Options, OutputNames, ...
   InputNames, OriginalHorizon, sysList)
% Generate compare plot for single-experiment IDDATA reference.
% IsReal: TRUE only for real time domain data

% Get plot object if exists
h = gcr(ax);
% Get plot add/replace status
NewPlot = strcmp(get(ax,'NextPlot'),'replace');
NewRespPlot = NewPlot || isempty(h);
IsIddata = isa(refData,'iddata');
IsFRD = isa(refData,'frd');
TD = IsIddata && strcmpi(refData.Domain,'time');
FD = IsIddata && ~TD;
Spect = IsFRD && size(refData,2)==0;

% todo
if TD
   PlotOptions = ltiplotoption('lsim',[],cstprefs.tbxprefs,NewRespPlot,h);
elseif Spect
   PlotOptions = ltiplotoption('noisespectrum',[],cstprefs.tbxprefs,NewRespPlot,h);
else
   PlotOptions = ltiplotoption('bode',[],cstprefs.tbxprefs,NewRespPlot,h);
end

if IsIddata
   %todo: When PlotOptions are exposed, make default InputVisible 'off'
   PlotOptions.InputVisible = {'off'};
end

% Check hold state
% Used to see if grid is on
if ~NewPlot && isempty(h) && ~isempty(findall(ax,'Tag','CSTgridLines'))
   PlotOptions.Grid = 'on';
end

% Clear and reset axes if new plot
if NewPlot

   % Notify the live editor of newplot
   matlab.graphics.internal.clearNotify(ax);

   % Clear any existing response plot upfront (otherwise style
   % settings below get erased by CLA in respplot/check_hold)
   if ~isempty(h)
      cla(h.AxesGrid,handle(ax))  % speed optimization
   end

   % Release manual limits and hide axis for optimal performance
   % RE: Manual negative Xlim can cause warning for BODE (not reset by clear)
   set(ax,'Visible','off','XlimMode','auto','YlimMode','auto')
end

if NewRespPlot
   % Ax can be a numeric so cast it as a handle
   axHndl = handle(ax);
   
   set(ax,...
      'XGrid',      PlotOptions.Grid,...
      'YGrid',      PlotOptions.Grid,...
      'FontSize',   PlotOptions.TickLabel.FontSize,...
      'FontWeight', PlotOptions.TickLabel.FontWeight,...
      'FontAngle',  PlotOptions.TickLabel.FontAngle,...
      'Selected',   'off')

   set(get(ax,'Title'),...
      'FontSize',  PlotOptions.Title.FontSize,...
      'FontWeight',PlotOptions.Title.FontWeight,...
      'FontAngle', PlotOptions.Title.FontAngle)

   set(get(ax,'XLabel'),...
      'FontSize',  PlotOptions.XLabel.FontSize,...
      'FontWeight',PlotOptions.XLabel.FontWeight,...
      'FontAngle', PlotOptions.XLabel.FontAngle)

   set(get(ax,'YLabel'),...
      'FontSize',  PlotOptions.YLabel.FontSize,...
      'FontWeight',PlotOptions.YLabel.FontWeight,...
      'FontAngle', PlotOptions.YLabel.FontAngle)

   if strcmp(PlotOptions.ColorMode.TickLabel,"manual")
        set(ax,'XColor',PlotOptions.TickLabel.Color,'YColor',PlotOptions.TickLabel.Color);
    end

    if strcmp(PlotOptions.ColorMode.XLabel,"manual")
        axHndl.XLabel.Color = PlotOptions.XLabel.Color;
    end

    if strcmp(PlotOptions.ColorMode.YLabel,"manual")
        axHndl.YLabel.Color = PlotOptions.YLabel.Color;
    end
end

% Create plot
nu = length(InputNames);
ny = length(OutputNames);
iosize = [ny nu];

p = PlotOptions;
if TD
   if isinf(OriginalHorizon)
      p.Title.String = getString(message('Ident:plots:msgSimComp'));
   else
      p.Title.String = getString(message('Ident:plots:msgPredComp',OriginalHorizon));
   end
   h = plotpack.tsimplot(ax, iosize, IsReal, 'OutputName', OutputNames, ...
      'InputName', InputNames, 'Tag','compareTD', 'PlotType', 'compare');
   setoptions(h,p)
elseif Spect
   InputNames = cellfun(@(x)['e@',x],OutputNames,'uni',false);
   p.Title.String = getString(message('Ident:plots:msgOutputSpecComp'));
   h = plotpack.fspectrumplot(ax, [iosize(1),iosize(1)], 'OutputName', ...
      OutputNames, 'InputName', InputNames, 'Tag','compareSpect',...
      'PlotType', 'compare');
   % TODO: when options are exposed, handle this by explicit
   % re-interpretation method for units = "auto"
   if strcmp(p.FreqUnits,'auto')
      tu = refData.TimeUnit;
      if strcmp(tu,'seconds')
         tu = 's';
      else
         tu = tu(1:end-1);
      end
      Unit = strrep(refData.FrequencyUnit,'TimeUnit',tu);
      p.FreqUnits = Unit;
   end
   setoptions(h,p);
else
   p.PhaseWrapping = 'on';
   if IsIddata
      p.Title.String = getString(message('Ident:plots:msgFFTComp'));
   else
      p.Title.String = getString(message('Ident:plots:msgFRComp'));
   end
   h = plotpack.fsimplot(ax, iosize, IsIddata, 'OutputName', ...
      OutputNames, 'InputName', InputNames, 'Tag','compareFD',...
      'PlotType', 'compare');
   % TODO: when options are exposed, handle this by explicit
   % re-interpretation method for units = "auto"
   if strcmp(p.FreqUnits,'auto')
      tu = refData.TimeUnit;
      if strcmp(tu,'seconds')
         tu = 's';
      else
         tu = tu(1:end-1);
      end
      Unit = strrep(refData.FrequencyUnit,'TimeUnit',tu);
      p.FreqUnits = Unit;
   end
   setoptions(h,p);
end

% Delete datatips when the axis is clicked
% check if figure is in Live Editor, if so remove th second
% toolbar (g3034662)
allaxes_h = allaxes(h);
if matlab.internal.editor.figure.FigureUtils.isEditorFigure(ancestor(ax,'figure'))
   if nu>0 && (TD || FD)
      for ku = 1:nu
         allaxes_h(ny+ku).Toolbar = [];
      end
   end   
end
set(allaxes_h,'ButtonDownFcn',{@LocalAxesButtonDownFcn h}) % Temporary workaround

%setoptions(h,PlotOptions); % REVISIT: create plot options!

if Spect
   localConfigureSpectPlot(h,src,Options,refData,sysList);
elseif FD || IsFRD
   localConfigureFDPlot(h,src,Options,refData,sysList);
else
   localConfigureTDPlot(h,src,Options,refData,sysList);
end

% Control cursor and datatip popups over characteristic markers
% REVISIT: remove this code when MouseEntered/Exited event available
Fig = ancestor(h.AxesGrid.Parent,'figure');
if isempty(get(Fig,'WindowButtonMotionFcn'))
   set(Fig,'WindowButtonMotionFcn',@(x,y) hoverfig(Fig))
   % Customize datacursor to use datatip style and not
   % cursor window
   hTool = datacursormode(Fig);
   set(hTool,'DisplayStyle','datatip');
end

if ~Spect && (FD || IsFRD)
   p.PhaseWrapping = 'off';
   setoptions(h,p);
end

% Turn off confidence region menu option for finite horizon plot
Menu = handle(findobj(h.AxesGrid.UIContextMenu,'Tag','ConfidenceRegion'));
if idIsValidHandle(Menu)
   if isinf(OriginalHorizon)
      Menu.Visible = 'on';
   else
      Menu.Visible = 'off';
   end
end

h.RefData.Context.ExpNo = 1;
h.draw;
h.Visible = 'on'; % required to compute response curves (before fit table data is computed)

%--------------------------------------------------------------------------
function localConfigureTDPlot(h,src,Options,refData,sysList)
% prepare the time domain compare plot

% Set global time focus for user-defined range/vector (sets preferred X limits)

TimeUnits = refData.TimeUnit;
t = pvget(refData,'SamplingInstants'); nexp = numel(t);
t = t{1}; % Exp 1 is default selection
h.setTimeFocus(t,TimeUnits)

vdatasrc = getDataPlotSource(refData,refData.Name);
vdatasrc.IsReal = h.IsReal;

% Configure validation data
rv = h.addRefDataWave(vdatasrc);
% add mean value and peak value characteristics
DefinedCharacteristics = vdatasrc.getCharacteristics('time');
rv.setCharacteristics(DefinedCharacteristics);
rv.DataFcn = {'getData', vdatasrc, rv};
if nexp>1
   rv.Name = sprintf('%s:%s',refData.Name,refData.ExperimentName{1});
end

% Configure model responses
for ksys = 1:length(src)
   % Link each response to system source
   Context = idguis.internal.plots.utGenDataModelPlotContext(sysList(ksys), nexp, size(refData,1), 'compare');
   if ~isempty(Context.DataOutputOrder) && all(isnan(Context.DataOutputOrder))
      continue;
   end
   if h.IsReal
      ColIndex = 1;
   else
      ColIndex = 1:2;
   end
   r = idguis.internal.plots.addDataModelWaveform(h, src(ksys), Context, ColIndex);
   DefinedCharacteristics = src(ksys).getCharacteristics('timeplot');
   r.setCharacteristics(DefinedCharacteristics);
   r.DataFcn = {'idguis.internal.plots.compare', r, vdatasrc};
   r.Context.Options = Options;
   r.Context.ExpNo = 1; % used for selecting proper IC
   % Styles and options
   initsysresp(r,'idsim',h.Options,r.Context.Style)
   for k = 1:numel(r.Data)
      r.Data(k).TimeUnits = TimeUnits;
   end
end

plotpack.simplotmenu(h,'time',size(refData,'Ne')>1);

%--------------------------------------------------------------------------
function localConfigureFDPlot(h,src,Options,refData,sysList)
% prepare the frequency domain compare plot
% comparison against FD iddata or frd

IsIddata = isa(refData,'iddata');
h.IsIddata = IsIddata;
if IsIddata
   nexp = size(refData,4);
   if nexp==1
      freq = refData.Frequency;
      fUnits = refData.FrequencyUnit;
   else
      % use first experiment data for instantiation
      freq = refData.Frequency{1};
      fUnits = refData.FrequencyUnit{1};
   end
   ncaps = size(refData,1);
else
   nexp = 1;
   freq = refData.Frequency;
   fUnits = refData.FrequencyUnit;
   ncaps = length(freq);
end
freq = freq(freq>=0);
fUnits = strrep(fUnits,'TimeUnit',refData.TimeUnit);
if numel(freq)==1
   h.setFreqFocus({0 max(freq,1)},fUnits)
else
   h.setFreqFocus(freq,fUnits)
end

vdatasrc = getDataPlotSource(refData,refData.Name);
% Configure reference data
rv = h.addRefDataWave(vdatasrc, IsIddata);
if IsIddata && nexp>1
   rv.Name = sprintf('%s:%s',refData.Name,refData.ExperimentName{1});
end

if IsIddata
   DefinedCharacteristics = vdatasrc.getCharacteristics('frequency');
else %frd
   DefinedCharacteristics = vdatasrc.getCharacteristics('freqplot');
end
rv.setCharacteristics(DefinedCharacteristics);
if IsIddata
   rv.DataFcn = {'iduis.utMagPhaseIOData', vdatasrc, rv};
else
   fac = funitconv(refData.FrequencyUnit,'rad/TimeUnit',refData.TimeUnit);
   rv.DataFcn = {'magphaseresp' vdatasrc 'bode' rv fac*refData.Frequency};
end

% Configure responses
for ksys = 1:length(src)
   % Link each response to system source
   Context = idguis.internal.plots.utGenDataModelPlotContext(sysList(ksys), nexp, ncaps, 'compare');
   if (~isempty(Context.DataInputOrder) && all(isnan(Context.DataInputOrder))) ||...
         (~isempty(Context.DataOutputOrder) && all(isnan(Context.DataOutputOrder)))
      continue;
   end
   if IsIddata
      r = idguis.internal.plots.addDataModelWaveform(h, src(ksys), Context, 1);
   else
      r = idguis.internal.plots.addDataModelWaveform(h, src(ksys), Context);
   end
   DefinedCharacteristics = src(ksys).getCharacteristics('freqplot');
   r.setCharacteristics(DefinedCharacteristics);
   r.DataFcn = {'idguis.internal.plots.compare', r, vdatasrc};
   r.Context.Options = Options;
   r.Context.ExpNo = 1;
   r.setstyle(r.Context.Style);
end

plotpack.simplotmenu(h,'frequency',IsIddata && size(refData,'Ne')>1);

%--------------------------------------------------------------------------
function localConfigureSpectPlot(h,src,Options,refData,sysList)
% Prepare the spectrum compare plot
fUnits = strrep(refData.FrequencyUnit,'TimeUnit',refData.TimeUnit);
freq = refData.Frequency;
freq = freq(freq>=0);
if numel(freq)==1
   h.setFreqFocus({0 max(freq,1)},fUnits)
else
   h.setFreqFocus(freq,fUnits)
end

vdatasrc = getDataPlotSource(refData,refData.Name);
rv = h.addRefDataWave(vdatasrc);
DefinedCharacteristics = vdatasrc.getCharacteristics('spectrumcompareplot');
rv.setCharacteristics(DefinedCharacteristics);
rv.DataFcn = {'iduis.utSpectrumData', vdatasrc, rv};

% Configure responses
for ksys = 1:length(src)
   % Link each response to system source
   Context = idguis.internal.plots.utGenDataModelPlotContext(sysList(ksys), 1, ...
      length(freq), 'compare');
   if (~isempty(Context.DataInputOrder) && all(isnan(Context.DataInputOrder))) ||...
         (~isempty(Context.DataOutputOrder) && all(isnan(Context.DataOutputOrder)))
      continue;
   end
   r = idguis.internal.plots.addDataModelWaveform(h, src(ksys), Context);
   DefinedCharacteristics = src(ksys).getCharacteristics('spectrumplot');
   r.setCharacteristics(DefinedCharacteristics);
   r.DataFcn = {'idguis.internal.plots.compare', r, vdatasrc};
   r.Context.Options = Options;
   r.Context.ExpNo = 1;
   r.setstyle(r.Context.Style);
end

plotpack.simplotmenu(h,'spectrum',false);

%--------------------------------------------------------------------------
function LocalAxesButtonDownFcn(EventSrc,~,~)
% Axes ButtonDown function to delete datatips when clicked
% Temporary workaround

RespPlot = gcr(EventSrc);
if ishandle(RespPlot)
   switch get(ancestor(EventSrc,'figure'),'SelectionType')
      case 'normal'
         PropEdit = PropEditor(RespPlot,'current');  % handle of (unique) property editor
         if ~isempty(PropEdit) && isvalid(PropEdit) && PropEdit.IsVisible
            % Left-click & property editor open: quick target change
            PropEdit.setTarget(RespPlot);
         end
         % Get the cursor mode object
         hTool = datacursormode(ancestor(EventSrc,'figure'));
         % Clear all data tips
         target = handle(EventSrc);
         if ishghandle(target,'axes')
            removeAllDataCursors(hTool,target);
         end
      case 'open'
         % Double-click: open editor
         PropEdit = PropEditor(RespPlot);
         PropEdit.setTarget(RespPlot);
         PropEdit.show();
   end
end


%--------------------------------------------------------------------------
function boo = localIsreal(IC)
% Check if initial conditions are real

boo = true;
if isnumeric(IC) || isa(IC,'param.Continuous')
   boo = isreal(IC);
elseif isstruct(IC) && isfield(IC,'Input') && isfield(IC,'Output')
   for ct = 1:numel(IC)
      boo = boo && isreal(IC(ct).Input) && isreal(IC(ct).Output);
   end
end
