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