--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@ao/melementOp.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,404 @@
+% MELEMENTOP applies the given matrix operator to the data.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% DESCRIPTION: MELEMENTOP applies the given matrix operator to the data.
+%
+% CALL:
+%              a = melementOp(callerIsMethod, op, opname, opsym, minfo, pl, a1, a2,...)
+%
+%
+% VERSION:     $Id: melementOp.m,v 1.12 2011/04/18 16:55:43 ingo Exp $
+%
+% HISTORY: 01-02-07 M Hewitson
+%             Creation
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function varargout = melementOp(varargin)
+  
+  import utils.const.*
+  
+  % Settings
+  callerIsMethod = varargin{1};
+  op     = varargin{2};
+  opname = varargin{3};
+  opsym  = varargin{4};
+  % Info to pass to history
+  iobj = varargin{5};
+  pl = varargin{6};
+  
+  % variable names
+  varnames = varargin{8};
+  
+  % Collect AO inputs but preserve the element shapes
+  % ... also collect numeric terms and preserve input names
+  argsin = varargin{7};
+  args = {};
+  in_names = {};
+  for kk=1:numel(argsin)
+    if isa(argsin{kk}, 'ao')
+      args = [args argsin(kk)];
+      in_names = [in_names varnames(kk)];
+    elseif isnumeric(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}), 0);
+      args = [args {a}];
+      if all(size(argsin{kk}) == [1 1])
+        in_names = [in_names num2str(argsin{kk})];
+      elseif any(size(argsin{kk}) == [1 1])
+        in_names = [in_names 'vector'];
+      else
+        in_names = [in_names 'matrix'];
+      end
+    end
+  end
+  
+  if numel(args) < 2
+    error('### %s operator requires at least two AO inputs.', opname)
+  end
+  
+  if numel(args) == 2
+    
+    % get the two arrays
+    a1 = args{1};
+    a2 = args{2};
+    
+    % check the data
+    for kk=1:numel(a1)
+      if ~isa(a1(kk).data, 'ltpda_data')
+        error('### one of the input AOs has an empty data field');
+      end
+    end
+    for kk=1:numel(a2)
+      if ~isa(a2(kk).data, 'ltpda_data')
+        error('### one of the input AOs has an empty data field');
+      end
+    end
+    
+    % Here we operate on two AO arrays according to the rules
+    
+    %---------- Deal with error cases first
+    r1 = size(a1,1);
+    c1 = size(a1,2);
+    r2 = size(a2,1);
+    c2 = size(a2,2);
+    
+    %== Rule 4: [1xN] */ [Nx1]
+    if r1 == 1 && r2 == 1 && c1==c2 && c1>1
+      error('### It is not possible to %s two AO vectors of size [1xN]', opname);
+    end
+    
+    %== Rule 6: [Nx1] */ [Nx1]
+    if r1 == r2 && c1==1 && c2==1 && r1>1
+      error('### It is not possible to %s two AO vectors of the size [Nx1]', opname);
+    end
+    
+    %== Rule 7: [NxP] */ [Nx1]
+    if r1 == r2 && c1>1 && c2==1 && c1~=r1 && r1>1
+      error('### It is not possible to %s [NxP] and [Nx1]', opname);
+    end
+    
+    %== Rule 8: [NxP] */ [Px1]
+    if c1 == c2 && r1>1 && r2==1 && c1>1
+      error('### It is not possible to %s [NxP] and [1xP]', opname);
+    end
+    
+    %== Rule 9: [NxP] */ [NxP]
+    if isequal(size(a1), size(a2)) && r1>1 && c1>1
+      if size(a1,1) ~= size(a1,2)
+        error('### It is not possible to %s [NxP] and [NxP]', opname);
+      end
+    end
+    
+    
+    %------------- Now perform operation
+    if numel(a1)==1 || numel(a2)==1
+      
+      % Rules 1,2,5
+      if isvector(a1) || isvector(a2) || ismatrix(a1) || ismatrix(a2)
+        % Rule 2,5: vector or matrix + single AO
+        if isvector(a1) || ismatrix(a1)
+          res = copy(a1,1);
+          for ee=1:numel(res)
+            res(ee).data = compatibleData(res(ee),a2);
+            res(ee).data.setY(operate(a1(ee), a2));
+            res(ee).data.setDy(operateError(a1(ee), a2));
+            % set history and name
+            if ~callerIsMethod
+              names = getNames(in_names, res(ee), ee, a2, []);
+              res(ee).addHistory(iobj, pl, names(1:2), [res(ee).hist a2.hist]);
+              res(ee).name = names{3};
+            end
+            res(ee).data.setYunits(getYunits(a1(ee), a2));
+          end
+        else
+          res = copy(a2,1);
+          for ee=1:numel(res)
+            res(ee).data = compatibleData(res(ee),a1);
+            res(ee).data.setY(operate(a2(ee), a1));
+            res(ee).data.setDy(operateError(a2(ee), a1));
+            % set history and name
+            if ~callerIsMethod
+              names = getNames(in_names, a1, [], res(ee), ee);
+              res(ee).addHistory(iobj, pl, names(1:2), [a1.hist res(ee).hist]);
+              res(ee).name = names{3};
+            end
+            res(ee).data.setYunits(getYunits(a1, a2(ee)));
+          end
+        end
+      else
+        % Rule 1: [1x1] */ [1x1]
+        res = copy(a1,1);
+        res.data = compatibleData(res,a2);
+        res.data.setY(operate(a1, a2));
+        res.data.setDy(operateError(a1, a2));
+        % set history and name
+        if ~callerIsMethod
+          names = getNames(in_names, res, [], a2, []);
+          res.addHistory(iobj, pl, names(1:2), [res.hist a2.hist]);
+          res.name = names{3};
+        end
+        res.data.setYunits(getYunits(a1, a2));
+      end
+    elseif isvector(a1) && isvector(a2) && r1==1 && c2==1 && r2==c1
+      % Rule 3: [1xN] */ [Nx1]
+      if strcmp(op, 'mrdivide')
+        error('### It is not possible to divide two matrices with different sizes');
+      end
+      res = [];
+      if strcmp(op, 'mtimes')
+        inner = 'times';
+      else
+        inner = 'rdivide';
+      end
+      
+      for ee=1:numel(a1)
+        if isempty(res)
+          res = feval(inner,a1(ee),a2(ee));
+        else
+          res = res + feval(inner,a1(ee),a2(ee));
+        end
+      end
+    elseif isvector(a1) && isvector(a2) && r1>1 && c1==1 && r2==1 && c2>1
+      % Rule 5: [Nx1] */ [1xM]
+      res(r1,c2) = ao();
+      for kk=1:r1
+        for ll=1:c2
+          res(kk,ll) = feval(op,a1(kk),a2(ll));
+        end
+      end
+    elseif ismatrix(a1) && (ismatrix(a2) || isvector(a2))
+      if strcmp(op, 'mrdivide') && ~isequal(size(a1),size(a2))
+        error('### Can only divide matrices of the same size');
+      end
+      % Rule 10: matrix */ matrix
+      res(r1,c2) = ao;
+      for kk=1:r1
+        for ll=1:c2
+          res(kk,ll) = feval(op,a1(kk,:),a2(:,ll));
+        end
+      end
+    else
+      error('### The inputs were not properly handled. This shouldn''t happen.');
+    end
+    
+    % Did something go wrong?
+    if isempty(res)
+      error('### The inputs were not properly handled. This shouldn''t happen.');
+    end
+    
+  else
+    % we recursively pass back to this method
+    res = copy(args{1}, 1);
+    for kk=2:numel(args)
+      res = feval(op, res, args{kk});
+    end
+  end
+  
+  % Set output
+  varargout{1} = res;
+  
+  %---------- nested functions
+  
+  %-------------------------------------------------
+  % Check the two inputs have compatible data types
+  function dout = compatibleData(a1,a2)
+    %== Data types
+    if (isa(a1.data, 'fsdata') && isa(a2.data, 'tsdata')) || ...
+        isa(a2.data, 'fsdata') && isa(a1.data, 'tsdata')
+      error('### Can not %s time-series data to frequency-series data.', opname);
+    end
+    % check X units for all data types
+    if ~isa(a1.data, 'cdata') && ~isa(a2.data, 'cdata')
+      if ~isempty(a1.data.xunits.strs) && ~isempty(a2.data.xunits.strs)
+        if a1.data.xunits ~= a2.data.xunits
+          error('### X units should be equal for the %s operator', op);
+        end
+      end
+    end
+    
+    % determine output data type
+    d1 = copy(a1.data,1);
+    d2 = copy(a2.data,1);
+    
+    if isa(d1, 'data2D') && isa(d2, 'data2D')
+      if numel(d1.y) > 1
+        dout = d1;
+      elseif numel(d2.y) > 1
+        dout = d2;
+      else
+        dout = d1;
+      end
+    elseif isa(d1, 'data2D') && isa(d2, 'cdata')
+      dout = d1;
+    elseif isa(d1, 'cdata') && isa(d2, 'data2D')
+      dout = d2;
+    else
+      dout = d1;
+    end
+    
+  end
+  
+  function uo = getYunits(a1, a2)
+    % For other operators we need to apply the operator
+    uo = feval(op, a1.data.yunits, a2.data.yunits);
+  end
+  
+  % Perform the desired operation on the data
+  function y = operate(a1, a2)
+    y = feval(op, a1.data.y, a2.data.y);
+  end
+  
+  % Perform the desired operation on the data uncertainty
+  function dy = operateError(a1, a2)
+    
+    if ~isempty(a1.dy) || ~isempty(a2.dy)
+      
+      da1 = a1.dy;
+      da2 = a2.dy;
+      
+      if isempty(da1)
+        da1 = zeros(size(a1.y));
+      end
+      if isempty(da2)
+        da2 = zeros(size(a2.y));
+      end
+      
+      switch op
+        case {'plus', 'minus'}
+          dy = sqrt(da1.^2 + da2.^2);
+        case {'times', 'mtimes'}
+          dy = sqrt( (da1./a1.y).^2 + (da2./a2.y).^2 ) .* abs(a1.y.*a2.y);
+        case {'rdivide', 'mrdivide'}
+          dy = sqrt( (da1./a1.y).^2 + (da2./a2.y).^2 ) .* abs(a1.y./a2.y);
+        otherwise
+          dy = [];
+      end
+      
+    else
+      dy = [];
+    end
+    
+  end
+  
+  %-----------------------------------------------
+  % Get two new AO names from the input var names,
+  % the input AO names, and the indices.
+  function names = getNames(in_names, a1, jj, a2, kk)
+    
+    % First variable name
+    if isempty(a1.name)  && ~isempty(in_names{1})
+      if ~isempty(jj)
+        if numel(jj) == 1
+          names{1} = sprintf('%s(%d)', in_names{1}, jj);
+        else
+          names{1} = sprintf('%s(%d,%d)', in_names{1}, jj(1), jj(2));
+        end
+      else
+        names{1} = in_names{1};
+      end
+    else
+      if ~isempty(jj)
+        if numel(jj) == 1
+          %           names{1} = sprintf('%s(%d)', a1.name, jj);
+          names{1} = sprintf('%s', a1.name);
+        else
+          %           names{1} = sprintf('%s(%d,%d)', a1.name, jj(1), jj(2));
+          names{1} = sprintf('%s', a1.name);
+        end
+      else
+        names{1} = a1.name;
+      end
+    end
+    % Second variable name
+    if isempty(a2.name) && ~isempty(in_names{2})
+      if isempty(in_names{2})
+        in_names{2} = a2.name;
+      end
+      if ~isempty(kk)
+        if numel(kk) == 1
+          names{2} = sprintf('%s(%d)', in_names{2}, kk);
+        else
+          names{1} = sprintf('%s(%d,%d)', in_names{2}, kk(1), kk(2));
+        end
+      else
+        names{2} = in_names{2};
+      end
+    else
+      names{2} = a2.name;
+      if ~isempty(kk)
+        if numel(kk) == 1
+          %           names{2} = sprintf('%s(%d)', a2.name, kk);
+          names{2} = sprintf('%s', a2.name);
+        else
+          %           names{2} = sprintf('%s(%d,%d)', a2.name, kk(1), kk(2));
+          names{2} = sprintf('%s', a2.name);
+        end
+      else
+        names{2} = a2.name;
+      end
+    end
+    
+    % The output AO name
+    names{3} = sprintf('(%s%s%s)', names{1}, opsym, names{2});
+  end
+  
+  %-------------------------------------
+  % Return true if the input is a matrix
+  function r = ismatrix(a)
+    if nrows(a) > 1 && ncols(a) > 1
+      r = true;
+    else
+      r = false;
+    end
+  end
+  
+  %-------------------------------------
+  % Return true if the input is a vector
+  function r = isvector(a)
+    if (nrows(a)==1 && ncols(a)>1) || (ncols(a)==1 && nrows(a)>1)
+      r = true;
+    else
+      r = false;
+    end
+  end
+  
+  %-------------------------------------
+  % Return numnber of rows in the array
+  function r = nrows(a)
+    r = size(a,1);
+  end
+  
+  %-------------------------------------
+  % Return numnber of cols in the array
+  function r = ncols(a)
+    r = size(a,2);
+  end
+  
+  
+end % End of add
+
+
+
+% END