Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/elementOp.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % ELEMENTOP applies the given operator to the data.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: ELEMENTOP applies the given operator to the data.
5 %
6 % CALL: a = elementOp(callerIsMethod, @getInfo, @getDefaultPlist, op, opname, opsym, aosNames, varargin(:));
7 %
8 % VERSION: $Id: elementOp.m,v 1.45 2011/04/18 16:53:12 ingo Exp $
9 %
10 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
11
12 function res = elementOp(varargin)
13
14
15 import utils.const.*
16
17 % Settings
18 callerIsMethod = varargin{1};
19 getInfo = varargin{2};
20 getDefaultPlist = varargin{3};
21
22 if callerIsMethod
23 dpl = [];
24 infoObj = [];
25 else
26 % Check if this is a call for parameters
27 if utils.helper.isinfocall(varargin{end}{:})
28 res = getInfo(varargin{end}{3});
29 return
30 end
31
32 dpl = getDefaultPlist();
33 infoObj = getInfo('None');
34
35 end
36
37 op = varargin{4};
38 opname = varargin{5};
39 opsym = varargin{6};
40
41 % variable names
42 varnames = varargin{7};
43
44 % Collect AO inputs but preserve the element shapes
45 % ... also collect numeric terms and preserve input names
46 argsin = varargin{8};
47
48 plin = [];
49 aos = {};
50 aosVarNames = {};
51
52 if numel(argsin) == 1
53 for kk=1:numel(argsin{1})
54 aos = [aos {argsin{1}(kk)}];
55 if ~callerIsMethod
56 aosVarNames = [aosVarNames varnames(1)];
57 end
58 end
59 else
60 for kk=1:numel(argsin)
61 if isa(argsin{kk}, 'ao')
62 aos = [aos argsin(kk)];
63 if ~callerIsMethod
64 aosVarNames = [aosVarNames varnames(kk)];
65 end
66 elseif isnumeric(argsin{kk}) || islogical(argsin{kk})
67 % When promoting the number to an AO, we have to be sure to call
68 % the fromVals and allow it to add history.
69 a = fromVals(ao, plist('vals', argsin{kk}), false);
70 aos = [aos {a}];
71 if all(size(argsin{kk}) == [1 1])
72 aosVarNames = [aosVarNames {num2str(argsin{kk})}];
73 elseif any(size(argsin{kk}) == [1 1])
74 aosVarNames = [aosVarNames 'vector'];
75 else
76 aosVarNames = [aosVarNames 'matrix'];
77 end
78 elseif isa(argsin{kk}, 'plist')
79 if isempty(plin)
80 plin = argsin{kk};
81 else
82 plin = combine(plin, argsin{kk});
83 end
84 end
85 end
86 end
87
88 % Combine input PLIST with default PLIST
89 if callerIsMethod
90 axis = 'y';
91 pl = [];
92 else
93 pl = applyDefaults(dpl, plin);
94 axis = pl.find('axis');
95 % operate at least the y-values
96 if isempty(axis)
97 axis = 'y';
98 end
99 end
100
101 if numel(aos) < 2
102 error('### A %s operator requires at least two AO inputs.', opname)
103 end
104
105 % utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), op);
106
107 if numel(aos) > 2
108
109 % we recursively pass back to this method
110 res = copy(aos{1}, 1);
111 resName = aosVarNames{1};
112 for kk=2:numel(aos)
113 res = ao.elementOp(callerIsMethod, getInfo, getDefaultPlist, op, opname, opsym, {resName, aosVarNames{kk}}, {res, aos{kk}});
114 resName = res.name;
115 end
116
117 else % args == 2
118
119 a1 = aos{1};
120 a2 = aos{2};
121
122 %%%%%%%%%% Rule 3:
123 if numel(a1) > 1 && numel(a2) > 1
124 if isVector(a1) && isVector(a2) && numel(a1) ~= numel(a2)
125 error('### It is not possible to %s two AO vectors of different lengths', opname);
126 end
127 end
128
129 %%%%%%%%%% Rule 8
130 if ismatrix(a1) && isVector(a2)
131 if nrows(a2) > 1 && nrows(a2) ~= nrows(a1)
132 error('### The number of rows in AO matrix should match the number of rows in the column vector.');
133 elseif ncols(a2)>1 && ncols(a2) ~= ncols(a1)
134 error('### The number of cols in AO matrix should match the number of cols in the row vector.');
135 end
136 end
137 if ismatrix(a2) && isVector(a1)
138 if nrows(a1) > 1 && nrows(a1) ~= nrows(a2)
139 error('### The number of rows in AO matrix should match the number of rows in the column vector.');
140 elseif ncols(a1)>1 && ncols(a1) ~= ncols(a2)
141 error('### The number of cols in AO matrix should match the number of cols in the row vector.');
142 end
143 end
144
145 %%%%%%%%%% Rule 9
146 if ismatrix(a1) && ismatrix(a2)
147 if ~isequal(size(a1), size(a2))
148 error('### Two AO matrices must be the same size to %s them.', opname);
149 end
150 end
151
152 %------------- Now perform operation
153
154 if numel(a1) == 1 && numel(a2) == 1
155
156 %%%%%%%%%% Rule 1: single AO + single AO
157 res = ao.initObjectWithSize(1,1);
158 operateSingleObject(res, a1, [], a2, []);
159
160 elseif isVector(a1) && numel(a2) == 1
161
162 %%%%%%%%%% Rule 2a: vector + single AO
163 res = ao.initObjectWithSize(size(a1));
164
165 for ii =1:numel(a1);
166 operateSingleObject(res(ii), a1(ii), ii, a2, []);
167 end
168
169 elseif numel(a1) == 1 && isVector(a2)
170
171 %%%%%%%%%% Rule 2b: single AO + vector
172 res = ao.initObjectWithSize(size(a2));
173
174 for ii =1:numel(a2);
175 operateSingleObject(res(ii), a1, [], a2(ii), ii);
176 end
177
178 elseif isVector(a1) && isVector(a2) && numel(a1) == numel(a2)
179
180 %%%%%%%%%% Rule 4: vector + vector
181 res = ao.initObjectWithSize(size(a1));
182
183 for ii =1:numel(a1);
184 operateSingleObject(res(ii), a1(ii), ii, a2(ii), ii);
185 end
186
187 elseif ismatrix(a1) && numel(a2) == 1
188
189 %%%%%%%%%% Rule 5a: matrix + single AO
190 res = ao.initObjectWithSize(size(a1));
191
192 for ii =1:numel(a1);
193 operateSingleObject(res(ii), a1(ii), ii, a2, []);
194 end
195
196 elseif numel(a1) == 1 && ismatrix(a2)
197
198 %%%%%%%%%% Rule 5b: single AO + matrix
199 res = ao.initObjectWithSize(size(a2));
200
201 for ii =1:numel(a2);
202 operateSingleObject(res(ii), a1, [], a2(ii), ii);
203 end
204
205 elseif ismatrix(a1) && isVector(a2) && size(a1,1) == length(a2)
206
207 %%%%%%%%%% Rule 6a: matrix NP + vector N
208 res = ao.initObjectWithSize(size(a1));
209
210 for nn = 1:size(a1,1)
211 for pp = 1:size(a1,2)
212 operateSingleObject(res(nn,pp), a1(nn,pp), [nn pp], a2(nn), nn);
213 end
214 end
215
216 elseif isVector(a1) && ismatrix(a2) && size(a2,1) == length(a1)
217
218 %%%%%%%%%% Rule 6b: vector N + matrix NP
219 res = ao.initObjectWithSize(size(a2));
220
221 for nn = 1:size(a2,1)
222 for pp = 1:size(a2,2)
223 operateSingleObject(res(nn,pp), a1(nn), nn, a2(nn,pp), [nn pp]);
224 end
225 end
226
227 elseif ismatrix(a1) && isVector(a2) && size(a1,2) == length(a2)
228
229 %%%%%%%%%% Rule 7a: matrix NP + vector P
230 res = ao.initObjectWithSize(size(a1));
231
232 for nn = 1:size(a1,1)
233 for pp = 1:size(a1,2)
234 operateSingleObject(res(nn,pp), a1(nn,pp), [nn pp], a2(pp), pp);
235 end
236 end
237
238 elseif isVector(a1) && ismatrix(a2) && size(a2,2) == length(a1)
239
240 %%%%%%%%%% Rule 7b: vector P + matrix NP
241 res = ao.initObjectWithSize(size(a2));
242
243 for nn = 1:size(a2,1)
244 for pp = 1:size(a2,2)
245 operateSingleObject(res(nn,pp), a1(pp), pp, a2(nn,pp), [nn pp]);
246 end
247 end
248
249 elseif ismatrix(a1) && ismatrix(a2) && size(a1,1) == size(a2,1) && size(a1,2) == size(a2,2)
250
251 %%%%%%%%%% Rule 10: matrix NP + matrix NP
252 res = ao.initObjectWithSize(size(a2));
253
254 for nn = 1:size(a1,1)
255 for pp = 1:size(a1,2)
256 operateSingleObject(res(nn,pp), a1(nn,pp), [nn pp], a2(nn,pp), [nn pp]);
257 end
258 end
259
260 else
261 error('### Should not happen.')
262 end
263
264 end
265
266 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
267 %
268 % DESCRIPTION: Applies the given operator to single input objects.
269 %
270 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
271
272 function operateSingleObject(res, a1, a1Idx, a2, a2Idx)
273
274 % Set data object
275 res.data = operateData(a1, a2, op, axis);
276
277 if callerIsMethod
278 % do nothing
279 else
280 % Set name
281 n1 = getName(a1.name, aosVarNames{1}, a1Idx);
282 n2 = getName(a2.name, aosVarNames{2}, a2Idx);
283 res.name = sprintf('(%s %s %s)', n1, opsym, n2);
284 % Set description
285 if ~isempty(a1.description) || ~isempty(a2.description)
286 if isempty(a1.description)
287 res.description = a2.description;
288 elseif isempty(a2.description)
289 res.description = a1.description;
290 else
291 res.description = strtrim([a1.description, ', ', a2.description]);
292 end
293 end
294 % Set plotinfo
295 if ~isempty(a1.plotinfo) || ~isempty(a2.plotinfo)
296 res.plotinfo = combine(a2.plotinfo, a1.plotinfo);
297 end
298 res.addHistory(infoObj, pl, {n1, n2}, [a1.hist a2.hist]);
299 end
300
301 end
302
303 end
304
305 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
306 % Local Functions %
307 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
308
309 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
310 %
311 % Applies the given operator to the data object.
312 %
313 function data = operateData(a1, a2, op, axis)
314
315 if isDataCompatible(a1, a2, op)
316
317 data = getDataObject(a1.data, a2.data, op);
318
319 operateValues(data, a1.data, a2.data, op, axis);
320 operateError(data, a1.data, a2.data, op, axis);
321 operateUnits(data, a1.data, a2.data, op, axis);
322 end
323 end
324
325 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
326 %
327 % Checks if the data objects are compatible
328 %
329 function res = isDataCompatible(a1, a2, op)
330
331 d1 = a1.data;
332 d2 = a2.data;
333 %%%%% check: Data types
334 if (isa(d1, 'fsdata') && isa(d2, 'tsdata')) || ...
335 isa(d2, 'fsdata') && isa(d1, 'tsdata')
336 error('### Can not operate time-series data to frequency-series data for the %s operator.', opname);
337 end
338 %%%%% check: Frequency for tsdata
339 if isa(d1, 'tsdata') && isa(d2, 'tsdata') && ...
340 d1.isprop('fs') && d2.isprop('fs') && ...
341 ~isempty(d1.fs) && ~isempty(d2.fs) && ...
342 ~isnan(d1.fs) && ~isnan(d2.fs) && ...
343 d1.fs ~= d2.fs
344 error('### The objects have different sample rates. Please resample one of the objects.')
345 end
346 if any(strcmpi(op, {'plus', 'minus'}))
347 %%%%% check: Y units
348 if ~isempty(d1.yunits.strs) && ~isempty(d2.yunits.strs)
349 if d1.yunits ~= d2.yunits
350 error('### Y units should be equal for the %s operator %s <-> %s', op, char(a1.yunits), char(a2.yunits));
351 end
352 end
353 end
354 %%%%% check: X units for all data types
355 if ~isa(d1, 'cdata') && ~isa(d2, 'cdata')
356 if ~isempty(d1.xunits.strs) && ~isempty(d2.xunits.strs)
357 if d1.xunits ~= d2.xunits
358 error('### X units should be equal for the %s operator', op);
359 end
360 end
361 end
362 %%%%% check x base of the fsdata objects
363 if isa(d1, 'fsdata') && isa(d2, 'fsdata') && any(abs(a1.x - a2.x) > 2*eps(a1.x))
364 error('### It is not possible to make any operation on frequency-series data if the x values are not the same.');
365 end
366
367 % Remove this condition because I haven't asked the other developers.
368 % %%%%% check: t0 of the tsdata objects
369 % if isa(d1, 'tsdata') && isa(d2, 'tsdata') && ~eq(d1.t0, d2.t0)
370 % error('### It is not possible to make any operation on time-series data if t0 is not the same.');
371 % end
372 res = true;
373 end
374
375 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
376 %
377 % Decide which data object should be used as the output object.
378 %
379 function dout = getDataObject(d1, d2, op)
380
381 % The output data object is always a copy.
382 if utils.helper.ismember(op, {'or', 'and', 'xor'})
383 dout = cdata.initObjectWithSize(1,1);
384 else
385 if isa(d1, 'data2D') && isa(d2, 'data2D')
386 if numel(d1.getY) > 1
387 dout = copy(d1, 1);
388 elseif numel(d2.getY) > 1
389 dout = copy(d2, 1);
390 else
391 dout = copy(d1, 1);
392 end
393 elseif isa(d1, 'data2D') && isa(d2, 'cdata')
394 dout = copy(d1, 1);
395 elseif isa(d1, 'cdata') && isa(d2, 'data2D')
396 dout = copy(d2, 1);
397 else
398 dout = copy(d1, 1);
399 end
400 end
401 end
402
403 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
404 %
405 % Evaluate the output values.
406 %
407 function operateValues(dout, d1, d2, op, axis)
408
409 if strcmp(op, 'mtimes') || strcmp(op, 'mrdivide')
410 y = feval(op, d1.y, d2.y);
411 else
412 y = feval(op, d1.getY, d2.getY);
413 end
414
415 if isa(dout, 'cdata')
416 if any(find(axis == 'y'))
417 dout.setY(y);
418 else
419 error('cdata objects only have y axis');
420 end
421 else
422 if any(find(axis == 'y'))
423 dout.setY(y);
424 end
425 if any(find(axis == 'x'))
426 dout.setX(feval(op, d1.getX, d2.getX));
427 end
428 end
429
430 end
431
432 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
433 %
434 % Evaluate the errors
435 %
436 function operateError(dout, d1, d2, op, axis)
437
438 % Define function for the errors
439 switch op
440 case {'plus', 'minus'}
441 err = @(err1, err2, val1, val2) sqrt(err1 .^2 + err2.^2);
442 case {'times', 'mtimes'}
443 err = @(err1, err2, val1, val2) sqrt( (err1./val1).^2 + (err2./val2).^2 ) .* abs(val1.*val2);
444 case {'rdivide', 'mrdivide'}
445 err = @(err1, err2, val1, val2) sqrt( (err1./val1).^2 + (err2./val2).^2 ) .* abs(val1./val2);
446 otherwise
447 err = @(err1, err2, val1, val2) [];
448 end
449
450 % Compute the error for the y-axis
451 if isa(dout, 'cdata')
452 if any(find(axis == 'y'))
453 if ~isempty(d1.dy) || ~isempty(d2.dy)
454
455 dy1 = d1.getDy;
456 dy2 = d2.getDy;
457
458 if isempty(dy1)
459 dy1 = zeros(size(d1.getY));
460 end
461 if isempty(dy2)
462 dy2 = zeros(size(d2.getY));
463 end
464
465 dy = err(dy1, dy2, d1.getY, d2.getY);
466 else
467 dy = [];
468 end
469 dout.setDy(dy);
470 else
471 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.');
472 end
473
474 else
475
476 % Compute the error for the y-axis
477 if any(find(axis == 'y'))
478 if ~isempty(d1.dy) || ~isempty(d2.dy)
479
480 dy1 = d1.getDy;
481 dy2 = d2.getDy;
482
483 if isempty(dy1)
484 dy1 = zeros(size(d1.getY));
485 end
486 if isempty(dy2)
487 dy2 = zeros(size(d2.getY));
488 end
489
490 dy = err(dy1, dy2, d1.getY, d2.getY);
491 else
492 dy = [];
493 end
494 dout.setDy(dy);
495 end
496
497 % % Compute the error for the x-axis
498 % if any(find(axis == 'x'))
499 % if ~isempty(d1.dx) || ~isempty(d2.dx)
500 %
501 % dx1 = d1.getDx;
502 % dx2 = d2.getDx;
503 %
504 % if isempty(dx1)
505 % dx1 = zeros(size(d1.getX));
506 % end
507 % if isempty(dx2)
508 % dx2 = zeros(size(d2.getX));
509 % end
510 %
511 % dx = err(dx1, dx2, d1.getX, d2.getX);
512 % else
513 % dx = [];
514 % end
515 % dout.setDx(dx);
516 % end
517
518 end
519 end
520
521 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
522 %
523 % Evaluate the units
524 %
525 function operateUnits(data, d1, d2, op, axis)
526 if utils.helper.ismember(op, {'or', 'and', 'xor'})
527 data.setYunits('');
528 elseif any(strcmpi(op, {'plus', 'minus'}))
529 % return the first non-empty
530 if ~isempty(d1.yunits.strs)
531 data.setYunits(d1.yunits);
532 else
533 data.setYunits(d2.yunits);
534 end
535 else
536 % For other operators we need to apply the operator
537 data.setYunits(feval(op, d1.yunits, d2.yunits));
538 end
539 end
540
541 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
542 %
543 function name = getName(objName, varName, idx)
544
545 if strcmpi(objName, 'none') || isempty(objName)
546 if ~isempty(varName)
547 useName = varName;
548 else
549 useName = objName;
550 end
551 % Set the name depending to the index
552 if isempty(idx)
553 name = useName;
554 elseif numel(idx) == 1
555 name = sprintf('%s(%d)', useName, idx(1));
556 else
557 name = sprintf('%s(%d,%d)', useName, idx(1), idx(2));
558 end
559 else
560 name = objName;
561 end
562
563 if isempty(name)
564 name = '?';
565 end
566
567 end
568
569 function res = isVector(a)
570 res = any(size(a) > 1) && any(size(a) == 1);
571 end
572
573 function res = ismatrix(a)
574 res = all(size(a) > 1);
575 end
576
577 %-------------------------------------
578 % Return number of rows in the array
579 function r = nrows(a)
580 r = size(a,1);
581 end
582
583 %-------------------------------------
584 % Return number of cols in the array
585 function r = ncols(a)
586 r = size(a,2);
587 end