line source
+ − % 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