function [Err, T_estimate] = e_sifir_tst(x_plus, Eamp, T, Param, Trim, varargin)
%E_SIFIR_TST Test a non-linear FIR EMG-torque model.
%
% [Err T_estimate] = e_sifir_tst(x_plus, Eamp, T, [Q D Tol ii], Trim)
%
% Computes the estimated torque(s) from a previously fit non-linear FIR
% EMG-torque model [presumably fit using function e_sifir_trn()]
% and the error(s) between the estimated and actual
% torque(s).  The default error measure is the RMS.  For multiple-output
% systems, errors may be computed along each output component (default)
% or in a distance sense. Full definitions of the required function
% inputs are provided in e_sifir_trn().  x_plus holds
% the model coefficients generated by e_sifir_trn().  Eamp and T are
% the data from the test trial(s), and Param and Trim are the same
% parameters used in training.  The Tol value of Param(3) is not used.
%
% If Eamp and T are NOT cell arrays,
% this function computes the EMG-predicted torque, T_estimate, using
% the entire USABLE sample range of Eamp.  This range will omit the
% first Q torque values (due to the lags needed in the FIR model) and
% ii additional values due to estimation into the future.  However, it is
% useful for T_estimate to have the same length as T.  Thus, the first
% useable T value is pre-pended (Q+ii) times so that T_estimate is
% returned to the caller with the same length as T, and so that each
% sample in T_estimate is properly time-aligned with samples in Eamp
% and T.
%
% After T_estimate is computed, the error (default = RMS) between it and T
% is computed [ignoring the first Q+ii+Trim(1) values and also the last
% Trim(2) values].  This error value is returned as the scalar Err,
% and T_estimate is returned as a vector. When T is an array (i.e., when
% there are multiple outputs), the error can be computed in two alternative
% manners.  By default, error is computed along each output component.
% If variable argument 'Edist' is set to 'distance, then a distance error
% is first computed at
% each time instant (treating each output as a dimension in the distance
% computation) and then the RMS or MAV is computed across time.
%
% When Eamp and T ARE cell arrays, then output T_estimate will be
% a cell array of the same dimensions.  Each element of T_estimate
% holds the corresponding EMG-predicted torque associated with the
% corresponding Eamp{} and T{} elements.  Output Err will be a regular
% matrix of the same dimensions as Eamp and T.  Each element of vector Err
% holds the error associated with the corresponding Eamp{}
% and T{} elements.
% 
% EMG Amplitude Estimation Toolbox - Edward (Ted) A. Clancy - WPI

% Copyright (c) Edward A. Clancy, 2015.
% This work is licensed under the Aladdin free public license.
% For copying permissions see license.txt.
% email: ted@wpi.edu

% 22 June 2015.

%%%%%%%%%%%%%%%%%%%%%%%%%% Process Command Line %%%%%%%%%%%%%%%%%%%%%%%%%%
% A few overall checks.
if nargin<5, error('Need >= 5 input arguments.'); end
if iscell(Eamp) == 0, Eamp = {Eamp}; end  % Coerce to cell array.
if iscell(T)    == 0, T    = {T};    end  % Coerce to cell array.
if sum(size(T)==size(Eamp))~=2, error('T and Eamp dimensions differ.'); end

% Eamp and T elements.
for el=2:numel(Eamp)  % Are all elements the same size?
  if sum(size(Eamp{1})==size(Eamp{el}))~=2, error('{Eamp} sizes differ.'); end
  if sum(size(T{1})   ==size(T{el}))   ~=2, error('{T} sizes differ.');    end
end
if length(Eamp{1})~=length(T{1}), error('Eamp and T time durations differ.'); end
for el=1:numel(Eamp)  % Coerce each Eamp as Eamp(ci,m); T as T(m,co).
  if size(Eamp{el},1) > size(Eamp{el},2), Eamp{el} = Eamp{el}'; end
  if size(T{el},1)    < size(T{el},2),    T{el}    = T{el}';    end
end

% Parameters.
if size(Param,1)~=1 && size(Param,2~=1), error('Param not a vector.'); end
if length(Param)~=4, error('length(Param) ~= 4.'); end
Q = Param(1);  D = Param(2);  ii = Param(4);
if  Q<0, error('Q less than zero.');  end
if  D<1, error('D less than one.');   end
if ii<0, error('ii less than zero.'); end

% Trim.
if size(Trim,1)~=1 && size(Trim,2~=1), error('Trim not a vector.'); end
if length(Trim)~=2, error('length(Trim) ~= 2.'); end
if Trim(1)<0, error('Trim(1) must be non-negative.'); end
if Trim(2)<0, error('Trim(2) must be non-negative.'); end

% Process optional commands, if any.
if rem(length(varargin),2)==1, error('Need even number of optional args.'); end
Edist = 'Component'; % Default method for multiple-output data.
Emeth = 'RMS';  % Default error computation method.

for i = 1:2:length(varargin)
  switch varargin{i}
      
    case 'Edist'  % Specify the method for multiple outputs.
      switch varargin{i+1}
        case {'Component', 'Distance'}, Edist = varargin{i+1};
        otherwise, error(['Bogus "Edist": "' varargin{i+1} '".'])
      end
    
    case 'Emeth'  % Specify the method for error computation.
      switch varargin{i+1}
        case {'MAV', 'RMS'}, Emeth = varargin{i+1};
        otherwise, error(['Bogus "Emeth": "' varargin{i+1} '".'])
      end
    
    otherwise, error(['Bogus PropertyName "' varargin{i} '".'])
      
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%% Compute. %%%%%%%%%%%%%%%%%%%%%%%%%%%%
T_estimate = cell(size(Eamp));  % Pre-allocate.
Err = zeros(size(Eamp));  % Pre-allocate.

for m = 1:numel(Eamp)  % Loop over the test data sets.
  % Build the test data design matix A using ALL available data.
  [A, b1] = e_sifir_ab(Eamp(m), T(m), Param);  % A  is (Nr-Q-ii) by (Q+1)*ci*D.
                                               % b1 is (Nr-Q-ii)*Ns by Nco.
  
  % Estimate last (Nr-Q-ii) torques.
  T_estimate{m} = A*x_plus;  % T_estimate is (Nr-Q-ii)*Ns by co.

  % Pre-pend Q+ii torques ==> T, T_estimate are aligned + have same length.
  T_estimate{m} = [ones(Q+ii,1)*T_estimate{m}(1,:); T_estimate{m}];
  b = [ones(Q+ii,1)*b1(1,:); b1];

  % Compute the desired error.
  Range = 1+Q+ii+Trim(1) : length(T{1})-Trim(2); % Computation index range.
  switch Edist  % Component vs. distance error for multiple channels.
    case 'Component'
      switch Emeth
        case 'MAV'
          Err(m) = mean(mean( abs(T_estimate{m}(Range,:)-b(Range,:)) ));
        case 'RMS'
          Err(m) = sqrt(  mean(mean( (T_estimate{m}(Range,:)-b(Range,:)).^2 ))  );
      end
    case 'Distance'
      E2 = (T_estimate{m}(Range,:)-b(Range,:)) .^ 2; % Array/vector: squared errors.
      if size(E2,2)>1, E2 = sum(E2,2); end  % Sum if there are multiple outputs.
      switch Emeth  % E2 is now a row vector of squared distance errors.
        case 'MAV',  Err(m) = mean(sqrt(E2));
        case 'RMS',  Err(m) = sqrt(mean(E2));
      end
  end
end

if length(Eamp)==1, T_estimate = T_estimate{1}; end  % If one test set, no cell.

return