function [check,S] = mgacmp(X,Xref,varargin) 
% function [check,S] = mgacmp(X,Xref,nodes,tol,tolV,verbose) 
% MGA results comparisons

 %% Argument Parsing
 p = inputParser;
 p.addRequired('X');
 p.addRequired('Xref');
 p.addOptional('nodes','all',@(x) ischar(x) || iscellstr(x));
 p.addOptional('tol'       ,  0.002);
 p.addOptional('tolV'       ,  0.01);
 p.addOptional('verb'      ,false);
 
 p.parse(X,Xref,varargin{:});
 
 X        = p.Results.X;
 Xref       = p.Results.Xref;
 nodes = p.Results.nodes;
 tol        = p.Results.tol;
 tolV        = p.Results.tolV;
 verb       = p.Results.verb;
 
 %% Function body
 if strcmp(nodes,'all')
  nodes = {'AAASW' 'G1SW' 'G2SW' 'G3SW' 'AMLSW' 'AIPOFT' 'BZEROFT' 'OHCUR' 'TOFT' 'ZEREF' 'CAYREF' 'EFCUR' 'EFVOLT' 'EFWAVE'};
 else
  nodes = cellstr(nodes);
 end
 
 % Compare
 check = true;
 
 % Resample if needed
 if numel(X.TOH) == numel(Xref.TOH) && max(abs(X.TOH - Xref.TOH)) < eps('single')
  t = Xref.TOH;
 else
  fprintv(verb,'... resampling ...\n');
  dt = (X.TOH(end) - X.TOH(1))/numel(X.TOH);
  dtref = (Xref.TOH(end) - Xref.TOH(1))/numel(Xref.TOH);
  t = max(X.TOH(1),Xref.TOH(1)):max(dt,dtref):min(X.TOH(end),Xref.TOH(end));
  assert(~isempty(t),'No common timebase between X and Xref');
 end
 
 for nodec = nodes
  node = nodec{1}; S.(node) = false;
  if strcmp(node,'EFVOLT'), rtol = tolV;
  else,                                  rtol = tol;
  end
  switch node
   case {'AAASW' 'AMLSW'}
    S.(node) = matcmp(node,Xref.(node)(:,:,1),X.(node)(:,:,1),rtol,verb);
   case {'G1SW' 'G2SW' 'G3SW'}
    check = check && size(Xref.(node),3) == size(X.(node),3);
    nsw = size(X.(node),3);
    S.(node) = matcmp(node,permute(Xref.(node)(:,:,1:nsw),[2 1 3]),permute(X.(node)(:,:,1:nsw),[2 1 3]),rtol,verb);
   case {'AIPOFT' 'BZEROFT' 'NTOFT' 'OHCUR' 'TOFT'}
    S.(node) = sigcmp(node,Xref.(node),X.(node),rtol,verb);
   case {'ZEREF' 'CAYREF' 'EFCUR' 'EFVOLT' 'EFWAVE'}
    [x,y] = resamp(Xref.(node),Xref.TOH,X.(node),X.TOH,t);
    S.(node) = sigcmp(node,x(2:end-1,:),y(2:end-1,:),rtol,verb); % Exclude first and last time points
  end
  check = check && all(abs(S.(node)(:)) < rtol);
 end
 
end

function s = sigcmp(n,x,y,tol,verb)
 kend = min(size(x,2),size(y,2));
 s = NaN(1,kend);
 for k = 1:kend
  if all(x(:,k) == 0)
   if all(y(:,k) == 0)
    e = 0;
   else
    e = Inf;
   end
  else
   e = norm(x(:,k)-y(:,k),2)./norm(x(:,k),2);
  end
  s(k) = e;
  if ~(e <= tol)
   fprintv(verb,'%s(:,%d) %5.2f%%\n',n,k,e*100)
  end
 end
end

function s = matcmp(n,X,Y,tol,verb)
 dim = [size(X) 1];
 ktlast = dim(3);
 %  for kt = ktlast-1:-1:1
 %   if isequal(X(:,:,kt),X(:,:,ktlast)) && isequal(Y(:,:,kt),Y(:,:,ktlast))
 %    ktlast = kt;
 %   end
 %  end
 s = zeros(dim(1),ktlast);
 for kt = 1:ktlast
  for k = 1:dim(1)
   if ~isequal(X(k,:,kt),Y(k,:,kt))
    if any(isnan(Y(k,:,kt)))
     l = find(isnan(Y(k,:,kt)));
     fprintv(verb,'%s(%d,%s,%d) = NaN\n',n,k,ind2str(l),kt)
    else
     if all(Y(k,:,kt) == 0), e = Inf;
     else,                   e = norm(Y(k,:,kt)' - X(k,:,kt)',2)./norm(X(k,:,kt)',2); end
     %     e = max([e sqrt(sum((Y(k,:,kt)-X(k,:,kt)).^2)/sum(Y(k,:,kt).^2))]);
     s(k,kt) = e;
     if ~(e <= tol)
      l = find(X(k,:,kt) ~= 0 | Y(k,:,kt) ~= 0);
      fprintv(verb,'%s(%d,%s,%d) %5.2f%% ^Xref vX\n',n,k,ind2str(l),kt,e*100)
      if verb, disp([Y(k,l,kt);X(k,l,kt)]); end
     end
    end
   end
  end
 end
 
end

function s = ind2str(l)
 switch length(l)
  case 1, s = num2str(l);
  case {0 2}, s = mat2str(l);
  otherwise
   k = [0 find(diff(l) ~= 1) length(l)];
   s = '[';
   for i = 1:length(k)-1
    s = [s num2str(l(k(i)+1))];
    if k(i)+1 ~= k(i+1)
     s = [s ':' num2str(l(k(i+1)))];
    end
    s = [s ' '];
   end
   s(end) = ']';
 end
end

function fprintv(verb,varargin)
% verbose printing
if verb
  fprintf(varargin{:});
end
end
