% ELEMENTOP applies the given operator to the data.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: ELEMENTOP applies the given operator to the data.
%
% CALL:        a = elementOp(callerIsMethod, @getInfo, @getDefaultPlist, op, opname, opsym, aosNames, varargin(:));
%
% VERSION:     $Id: elementOp.m,v 1.45 2011/04/18 16:53:12 ingo Exp $
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function res = elementOp(varargin)
 
  
  import utils.const.*
  
  % Settings
  callerIsMethod = varargin{1};
  getInfo = varargin{2};
  getDefaultPlist = varargin{3};

  if callerIsMethod
    dpl = [];
    infoObj = [];
  else    
    % Check if this is a call for parameters
    if utils.helper.isinfocall(varargin{end}{:})
      res = getInfo(varargin{end}{3});
      return
    end
    
    dpl = getDefaultPlist();
    infoObj = getInfo('None');
        
  end
  
  op     = varargin{4};
  opname = varargin{5};
  opsym  = varargin{6};
    
  % variable names
  varnames = varargin{7};
  
  % Collect AO inputs but preserve the element shapes
  % ... also collect numeric terms and preserve input names
  argsin = varargin{8};
  
  plin = [];
  aos         = {};
  aosVarNames = {};
  
  if numel(argsin) == 1
    for kk=1:numel(argsin{1})
      aos = [aos {argsin{1}(kk)}];
      if ~callerIsMethod
        aosVarNames = [aosVarNames varnames(1)];
      end
    end
  else
    for kk=1:numel(argsin)
      if isa(argsin{kk}, 'ao')
        aos = [aos argsin(kk)];
        if ~callerIsMethod
          aosVarNames = [aosVarNames varnames(kk)];
        end
      elseif isnumeric(argsin{kk}) || islogical(argsin{kk})
        % When promoting the number to an AO, we have to be sure to call
        % the fromVals and allow it to add history.
        a = fromVals(ao, plist('vals', argsin{kk}), false);
        aos = [aos {a}];
        if all(size(argsin{kk}) == [1 1])
          aosVarNames = [aosVarNames {num2str(argsin{kk})}];
        elseif any(size(argsin{kk}) == [1 1])
          aosVarNames = [aosVarNames 'vector'];
        else
          aosVarNames = [aosVarNames 'matrix'];
        end
      elseif isa(argsin{kk}, 'plist')
        if isempty(plin)
          plin = argsin{kk};
        else
          plin = combine(plin, argsin{kk});
        end
      end
    end
  end
  
  % Combine input PLIST with default PLIST
  if callerIsMethod
    axis = 'y';
    pl = [];
  else    
    pl = applyDefaults(dpl, plin);
    axis = pl.find('axis');
    % operate at least the y-values
    if isempty(axis)
      axis = 'y';
    end
  end  
  
  if numel(aos) < 2
    error('### A %s operator requires at least two AO inputs.', opname)
  end
  
%   utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), op);
  
  if numel(aos) > 2
    
    % we recursively pass back to this method
    res     = copy(aos{1}, 1);
    resName = aosVarNames{1};
    for kk=2:numel(aos)
      res = ao.elementOp(callerIsMethod, getInfo, getDefaultPlist, op, opname, opsym, {resName, aosVarNames{kk}}, {res, aos{kk}});
      resName = res.name;
    end
    
  else % args == 2
    
    a1 = aos{1};
    a2 = aos{2};
    
    %%%%%%%%%%   Rule 3:
    if numel(a1) > 1 && numel(a2) > 1
      if isVector(a1) && isVector(a2) && numel(a1) ~= numel(a2)
        error('### It is not possible to %s two AO vectors of different lengths', opname);
      end
    end
    
    %%%%%%%%%%   Rule 8
    if ismatrix(a1) && isVector(a2)
      if nrows(a2) > 1 && nrows(a2) ~= nrows(a1)
        error('### The number of rows in AO matrix should match the number of rows in the column vector.');
      elseif ncols(a2)>1 && ncols(a2) ~= ncols(a1)
        error('### The number of cols in AO matrix should match the number of cols in the row vector.');
      end
    end
    if ismatrix(a2) && isVector(a1)
      if nrows(a1) > 1 && nrows(a1) ~= nrows(a2)
        error('### The number of rows in AO matrix should match the number of rows in the column vector.');
      elseif ncols(a1)>1 && ncols(a1) ~= ncols(a2)
        error('### The number of cols in AO matrix should match the number of cols in the row vector.');
      end
    end
    
    %%%%%%%%%%   Rule 9
    if ismatrix(a1) && ismatrix(a2)
      if ~isequal(size(a1), size(a2))
        error('### Two AO matrices must be the same size to %s them.', opname);
      end
    end
    
    %------------- Now perform operation
    
    if numel(a1) == 1 && numel(a2) == 1
      
      %%%%%%%%%%   Rule 1: single AO + single AO
      res = ao.initObjectWithSize(1,1);
      operateSingleObject(res, a1, [], a2, []);
      
    elseif isVector(a1) && numel(a2) == 1
      
      %%%%%%%%%%   Rule 2a: vector + single AO
      res = ao.initObjectWithSize(size(a1));
      
      for ii =1:numel(a1);
        operateSingleObject(res(ii), a1(ii), ii, a2, []);
      end
      
    elseif numel(a1) == 1 && isVector(a2)
      
      %%%%%%%%%%   Rule 2b: single AO + vector
      res = ao.initObjectWithSize(size(a2));
      
      for ii =1:numel(a2);
        operateSingleObject(res(ii), a1, [], a2(ii), ii);
      end
      
    elseif isVector(a1) && isVector(a2) && numel(a1) == numel(a2)
      
      %%%%%%%%%%   Rule 4: vector + vector
      res = ao.initObjectWithSize(size(a1));
      
      for ii =1:numel(a1);
        operateSingleObject(res(ii), a1(ii), ii, a2(ii), ii);
      end
      
    elseif ismatrix(a1) && numel(a2) == 1
      
      %%%%%%%%%%   Rule 5a: matrix + single AO
      res = ao.initObjectWithSize(size(a1));
      
      for ii =1:numel(a1);
        operateSingleObject(res(ii), a1(ii), ii, a2, []);
      end
      
    elseif numel(a1) == 1 && ismatrix(a2)
      
      %%%%%%%%%%   Rule 5b: single AO + matrix
      res = ao.initObjectWithSize(size(a2));
      
      for ii =1:numel(a2);
        operateSingleObject(res(ii), a1, [], a2(ii), ii);
      end
      
    elseif ismatrix(a1) && isVector(a2) && size(a1,1) == length(a2)
      
      %%%%%%%%%%   Rule 6a: matrix NP + vector N
      res = ao.initObjectWithSize(size(a1));
      
      for nn = 1:size(a1,1)
        for pp = 1:size(a1,2)
          operateSingleObject(res(nn,pp), a1(nn,pp), [nn pp], a2(nn), nn);
        end
      end
      
    elseif isVector(a1) && ismatrix(a2) && size(a2,1) == length(a1)
      
      %%%%%%%%%%   Rule 6b: vector N + matrix NP
      res = ao.initObjectWithSize(size(a2));
      
      for nn = 1:size(a2,1)
        for pp = 1:size(a2,2)
          operateSingleObject(res(nn,pp), a1(nn), nn, a2(nn,pp), [nn pp]);
        end
      end
      
    elseif ismatrix(a1) && isVector(a2) && size(a1,2) == length(a2)
      
      %%%%%%%%%%   Rule 7a: matrix NP + vector P
      res = ao.initObjectWithSize(size(a1));
      
      for nn = 1:size(a1,1)
        for pp = 1:size(a1,2)
          operateSingleObject(res(nn,pp), a1(nn,pp), [nn pp], a2(pp), pp);
        end
      end
      
    elseif isVector(a1) && ismatrix(a2) && size(a2,2) == length(a1)
      
      %%%%%%%%%%   Rule 7b: vector P + matrix NP
      res = ao.initObjectWithSize(size(a2));
      
      for nn = 1:size(a2,1)
        for pp = 1:size(a2,2)
          operateSingleObject(res(nn,pp), a1(pp), pp, a2(nn,pp), [nn pp]);
        end
      end
      
    elseif ismatrix(a1) && ismatrix(a2) && size(a1,1) == size(a2,1) && size(a1,2) == size(a2,2)
      
      %%%%%%%%%%   Rule 10: matrix NP + matrix NP
      res = ao.initObjectWithSize(size(a2));
      
      for nn = 1:size(a1,1)
        for pp = 1:size(a1,2)
          operateSingleObject(res(nn,pp), a1(nn,pp), [nn pp], a2(nn,pp), [nn pp]);
        end
      end
      
    else
      error('### Should not happen.')
    end
    
  end
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %
  % DESCRIPTION: Applies the given operator to single input objects.
  %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  function operateSingleObject(res, a1, a1Idx, a2, a2Idx)
    
    % Set data object
    res.data = operateData(a1, a2, op, axis);
        
    if callerIsMethod
      % do nothing
    else
      % Set name
      n1 = getName(a1.name, aosVarNames{1}, a1Idx);
      n2 = getName(a2.name, aosVarNames{2}, a2Idx);
      res.name = sprintf('(%s %s %s)', n1, opsym, n2);
      % Set description
      if ~isempty(a1.description) || ~isempty(a2.description)
        if isempty(a1.description)
          res.description = a2.description;
        elseif isempty(a2.description)
          res.description = a1.description;
        else
          res.description = strtrim([a1.description, ', ', a2.description]);
        end
      end
      % Set plotinfo
      if ~isempty(a1.plotinfo) || ~isempty(a2.plotinfo)
        res.plotinfo = combine(a2.plotinfo, a1.plotinfo);
      end
      res.addHistory(infoObj, pl, {n1, n2}, [a1.hist a2.hist]);
    end
    
  end
  
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               Local Functions                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Applies the given operator to the data object.
%
function data = operateData(a1, a2, op, axis)
  
  if isDataCompatible(a1, a2, op)
    
    data = getDataObject(a1.data, a2.data, op);
    
    operateValues(data, a1.data, a2.data, op, axis);
    operateError(data, a1.data, a2.data, op, axis);
    operateUnits(data, a1.data, a2.data, op, axis);
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Checks if the data objects are compatible
%
function res = isDataCompatible(a1, a2, op)
  
  d1 = a1.data;
  d2 = a2.data;
  %%%%% check: Data types
  if (isa(d1, 'fsdata') && isa(d2, 'tsdata')) || ...
      isa(d2, 'fsdata') && isa(d1, 'tsdata')
    error('### Can not operate time-series data to frequency-series data for the %s operator.', opname);
  end
  %%%%% check: Frequency for tsdata
  if  isa(d1, 'tsdata') && isa(d2, 'tsdata') && ...
      d1.isprop('fs') && d2.isprop('fs') && ...
      ~isempty(d1.fs) && ~isempty(d2.fs) && ...
      ~isnan(d1.fs) && ~isnan(d2.fs) && ...
      d1.fs ~= d2.fs
    error('### The objects have different sample rates. Please resample one of the objects.')
  end
  if any(strcmpi(op, {'plus', 'minus'}))
    %%%%% check: Y units
    if ~isempty(d1.yunits.strs) && ~isempty(d2.yunits.strs)
      if d1.yunits ~= d2.yunits
        error('### Y units should be equal for the %s operator %s <-> %s', op, char(a1.yunits), char(a2.yunits));
      end
    end
  end
  %%%%% check: X units for all data types
  if ~isa(d1, 'cdata') && ~isa(d2, 'cdata')
    if ~isempty(d1.xunits.strs) && ~isempty(d2.xunits.strs)
      if d1.xunits ~= d2.xunits
        error('### X units should be equal for the %s operator', op);
      end
    end
  end
  %%%%% check x base of the fsdata objects
  if isa(d1, 'fsdata') && isa(d2, 'fsdata') && any(abs(a1.x - a2.x) > 2*eps(a1.x))
    error('### It is not possible to make any operation on frequency-series data if the x values are not the same.');
  end
  
% Remove this condition because I haven't asked the other developers.
%   %%%%% check: t0 of the tsdata objects
%   if isa(d1, 'tsdata') && isa(d2, 'tsdata') && ~eq(d1.t0, d2.t0)
%     error('### It is not possible to make any operation on time-series data if t0 is not the same.');
%   end
  res = true;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Decide which data object should be used as the output object.
%
function dout = getDataObject(d1, d2, op)
  
  % The output data object is always a copy.
  if utils.helper.ismember(op, {'or', 'and', 'xor'})
    dout = cdata.initObjectWithSize(1,1);
  else
    if isa(d1, 'data2D') && isa(d2, 'data2D')
      if numel(d1.getY) > 1
        dout = copy(d1, 1);
      elseif numel(d2.getY) > 1
        dout = copy(d2, 1);
      else
        dout = copy(d1, 1);
      end
    elseif isa(d1, 'data2D') && isa(d2, 'cdata')
      dout = copy(d1, 1);
    elseif isa(d1, 'cdata') && isa(d2, 'data2D')
      dout = copy(d2, 1);
    else
      dout = copy(d1, 1);
    end
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Evaluate the output values.
%
function operateValues(dout, d1, d2, op, axis)
  
  if strcmp(op, 'mtimes') || strcmp(op, 'mrdivide')
    y = feval(op, d1.y, d2.y);
  else
    y = feval(op, d1.getY, d2.getY);
  end
  
  if isa(dout, 'cdata')
    if any(find(axis == 'y'))
      dout.setY(y);
    else
      error('cdata objects only have y axis');
    end
  else
    if any(find(axis == 'y'))
      dout.setY(y);
    end
    if any(find(axis == 'x'))
      dout.setX(feval(op, d1.getX, d2.getX));
    end
  end
  
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Evaluate the errors
%
function operateError(dout, d1, d2, op, axis)
  
  % Define function for the errors
  switch op
    case {'plus', 'minus'}
      err = @(err1, err2, val1, val2) sqrt(err1 .^2 + err2.^2);
    case {'times', 'mtimes'}
      err = @(err1, err2, val1, val2) sqrt( (err1./val1).^2 + (err2./val2).^2 ) .* abs(val1.*val2);
    case {'rdivide', 'mrdivide'}
      err = @(err1, err2, val1, val2) sqrt( (err1./val1).^2 + (err2./val2).^2 ) .* abs(val1./val2);
    otherwise
      err = @(err1, err2, val1, val2) [];
  end
  
  % Compute the error for the y-axis
  if isa(dout, 'cdata')
    if any(find(axis == 'y'))
      if ~isempty(d1.dy) || ~isempty(d2.dy)
        
        dy1 = d1.getDy;
        dy2 = d2.getDy;
        
        if isempty(dy1)
          dy1 = zeros(size(d1.getY));
        end
        if isempty(dy2)
          dy2 = zeros(size(d2.getY));
        end
        
        dy = err(dy1, dy2, d1.getY, d2.getY);
      else
        dy = [];
      end
      dout.setDy(dy);
    else
      warning('!!! The output data object is a ''cdata'' object and you operate only on the x-axis but this axis doesn''t exist on a cdata object.');
    end
    
  else
    
    % Compute the error for the y-axis
    if any(find(axis == 'y'))
      if ~isempty(d1.dy) || ~isempty(d2.dy)
        
        dy1 = d1.getDy;
        dy2 = d2.getDy;
        
        if isempty(dy1)
          dy1 = zeros(size(d1.getY));
        end
        if isempty(dy2)
          dy2 = zeros(size(d2.getY));
        end
        
        dy = err(dy1, dy2, d1.getY, d2.getY);
      else
        dy = [];
      end
      dout.setDy(dy);
    end
    
%     % Compute the error for the x-axis
%     if any(find(axis == 'x'))
%       if ~isempty(d1.dx) || ~isempty(d2.dx)
%         
%         dx1 = d1.getDx;
%         dx2 = d2.getDx;
%         
%         if isempty(dx1)
%           dx1 = zeros(size(d1.getX));
%         end
%         if isempty(dx2)
%           dx2 = zeros(size(d2.getX));
%         end
%         
%         dx = err(dx1, dx2, d1.getX, d2.getX);
%       else
%         dx = [];
%       end
%       dout.setDx(dx);
%     end
    
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Evaluate the units
%
function operateUnits(data, d1, d2, op, axis)
  if utils.helper.ismember(op, {'or', 'and', 'xor'})
    data.setYunits('');
  elseif any(strcmpi(op, {'plus', 'minus'}))
    % return the first non-empty
    if ~isempty(d1.yunits.strs)
      data.setYunits(d1.yunits);
    else
      data.setYunits(d2.yunits);
    end
  else
    % For other operators we need to apply the operator
    data.setYunits(feval(op, d1.yunits, d2.yunits));
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
function name = getName(objName, varName, idx)
  
  if strcmpi(objName, 'none') || isempty(objName)
    if ~isempty(varName)
      useName = varName;
    else
      useName = objName;
    end
    % Set the name depending to the index
    if isempty(idx)
      name = useName;
    elseif numel(idx) == 1
      name = sprintf('%s(%d)', useName, idx(1));
    else
      name = sprintf('%s(%d,%d)', useName, idx(1), idx(2));
    end
  else
    name = objName;
  end
  
  if isempty(name)
    name = '?';
  end
  
end

function res = isVector(a)
  res = any(size(a) > 1) && any(size(a) == 1);
end

function res = ismatrix(a)
  res = all(size(a) > 1);
end

%-------------------------------------
% Return number of rows in the array
function r = nrows(a)
  r = size(a,1);
end

%-------------------------------------
% Return number of cols in the array
function r = ncols(a)
  r = size(a,2);
end