Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/melementOp.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 % MELEMENTOP applies the given matrix operator to the data.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: MELEMENTOP applies the given matrix operator to the data.
5 %
6 % CALL:
7 % a = melementOp(callerIsMethod, op, opname, opsym, minfo, pl, a1, a2,...)
8 %
9 %
10 % VERSION: $Id: melementOp.m,v 1.12 2011/04/18 16:55:43 ingo Exp $
11 %
12 % HISTORY: 01-02-07 M Hewitson
13 % Creation
14 %
15 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
16
17 function varargout = melementOp(varargin)
18
19 import utils.const.*
20
21 % Settings
22 callerIsMethod = varargin{1};
23 op = varargin{2};
24 opname = varargin{3};
25 opsym = varargin{4};
26 % Info to pass to history
27 iobj = varargin{5};
28 pl = varargin{6};
29
30 % variable names
31 varnames = varargin{8};
32
33 % Collect AO inputs but preserve the element shapes
34 % ... also collect numeric terms and preserve input names
35 argsin = varargin{7};
36 args = {};
37 in_names = {};
38 for kk=1:numel(argsin)
39 if isa(argsin{kk}, 'ao')
40 args = [args argsin(kk)];
41 in_names = [in_names varnames(kk)];
42 elseif isnumeric(argsin{kk})
43 % When promoting the number to an AO, we have to be sure to call
44 % the fromVals and allow it to add history.
45 a = fromVals(ao, plist('vals', argsin{kk}), 0);
46 args = [args {a}];
47 if all(size(argsin{kk}) == [1 1])
48 in_names = [in_names num2str(argsin{kk})];
49 elseif any(size(argsin{kk}) == [1 1])
50 in_names = [in_names 'vector'];
51 else
52 in_names = [in_names 'matrix'];
53 end
54 end
55 end
56
57 if numel(args) < 2
58 error('### %s operator requires at least two AO inputs.', opname)
59 end
60
61 if numel(args) == 2
62
63 % get the two arrays
64 a1 = args{1};
65 a2 = args{2};
66
67 % check the data
68 for kk=1:numel(a1)
69 if ~isa(a1(kk).data, 'ltpda_data')
70 error('### one of the input AOs has an empty data field');
71 end
72 end
73 for kk=1:numel(a2)
74 if ~isa(a2(kk).data, 'ltpda_data')
75 error('### one of the input AOs has an empty data field');
76 end
77 end
78
79 % Here we operate on two AO arrays according to the rules
80
81 %---------- Deal with error cases first
82 r1 = size(a1,1);
83 c1 = size(a1,2);
84 r2 = size(a2,1);
85 c2 = size(a2,2);
86
87 %== Rule 4: [1xN] */ [Nx1]
88 if r1 == 1 && r2 == 1 && c1==c2 && c1>1
89 error('### It is not possible to %s two AO vectors of size [1xN]', opname);
90 end
91
92 %== Rule 6: [Nx1] */ [Nx1]
93 if r1 == r2 && c1==1 && c2==1 && r1>1
94 error('### It is not possible to %s two AO vectors of the size [Nx1]', opname);
95 end
96
97 %== Rule 7: [NxP] */ [Nx1]
98 if r1 == r2 && c1>1 && c2==1 && c1~=r1 && r1>1
99 error('### It is not possible to %s [NxP] and [Nx1]', opname);
100 end
101
102 %== Rule 8: [NxP] */ [Px1]
103 if c1 == c2 && r1>1 && r2==1 && c1>1
104 error('### It is not possible to %s [NxP] and [1xP]', opname);
105 end
106
107 %== Rule 9: [NxP] */ [NxP]
108 if isequal(size(a1), size(a2)) && r1>1 && c1>1
109 if size(a1,1) ~= size(a1,2)
110 error('### It is not possible to %s [NxP] and [NxP]', opname);
111 end
112 end
113
114
115 %------------- Now perform operation
116 if numel(a1)==1 || numel(a2)==1
117
118 % Rules 1,2,5
119 if isvector(a1) || isvector(a2) || ismatrix(a1) || ismatrix(a2)
120 % Rule 2,5: vector or matrix + single AO
121 if isvector(a1) || ismatrix(a1)
122 res = copy(a1,1);
123 for ee=1:numel(res)
124 res(ee).data = compatibleData(res(ee),a2);
125 res(ee).data.setY(operate(a1(ee), a2));
126 res(ee).data.setDy(operateError(a1(ee), a2));
127 % set history and name
128 if ~callerIsMethod
129 names = getNames(in_names, res(ee), ee, a2, []);
130 res(ee).addHistory(iobj, pl, names(1:2), [res(ee).hist a2.hist]);
131 res(ee).name = names{3};
132 end
133 res(ee).data.setYunits(getYunits(a1(ee), a2));
134 end
135 else
136 res = copy(a2,1);
137 for ee=1:numel(res)
138 res(ee).data = compatibleData(res(ee),a1);
139 res(ee).data.setY(operate(a2(ee), a1));
140 res(ee).data.setDy(operateError(a2(ee), a1));
141 % set history and name
142 if ~callerIsMethod
143 names = getNames(in_names, a1, [], res(ee), ee);
144 res(ee).addHistory(iobj, pl, names(1:2), [a1.hist res(ee).hist]);
145 res(ee).name = names{3};
146 end
147 res(ee).data.setYunits(getYunits(a1, a2(ee)));
148 end
149 end
150 else
151 % Rule 1: [1x1] */ [1x1]
152 res = copy(a1,1);
153 res.data = compatibleData(res,a2);
154 res.data.setY(operate(a1, a2));
155 res.data.setDy(operateError(a1, a2));
156 % set history and name
157 if ~callerIsMethod
158 names = getNames(in_names, res, [], a2, []);
159 res.addHistory(iobj, pl, names(1:2), [res.hist a2.hist]);
160 res.name = names{3};
161 end
162 res.data.setYunits(getYunits(a1, a2));
163 end
164 elseif isvector(a1) && isvector(a2) && r1==1 && c2==1 && r2==c1
165 % Rule 3: [1xN] */ [Nx1]
166 if strcmp(op, 'mrdivide')
167 error('### It is not possible to divide two matrices with different sizes');
168 end
169 res = [];
170 if strcmp(op, 'mtimes')
171 inner = 'times';
172 else
173 inner = 'rdivide';
174 end
175
176 for ee=1:numel(a1)
177 if isempty(res)
178 res = feval(inner,a1(ee),a2(ee));
179 else
180 res = res + feval(inner,a1(ee),a2(ee));
181 end
182 end
183 elseif isvector(a1) && isvector(a2) && r1>1 && c1==1 && r2==1 && c2>1
184 % Rule 5: [Nx1] */ [1xM]
185 res(r1,c2) = ao();
186 for kk=1:r1
187 for ll=1:c2
188 res(kk,ll) = feval(op,a1(kk),a2(ll));
189 end
190 end
191 elseif ismatrix(a1) && (ismatrix(a2) || isvector(a2))
192 if strcmp(op, 'mrdivide') && ~isequal(size(a1),size(a2))
193 error('### Can only divide matrices of the same size');
194 end
195 % Rule 10: matrix */ matrix
196 res(r1,c2) = ao;
197 for kk=1:r1
198 for ll=1:c2
199 res(kk,ll) = feval(op,a1(kk,:),a2(:,ll));
200 end
201 end
202 else
203 error('### The inputs were not properly handled. This shouldn''t happen.');
204 end
205
206 % Did something go wrong?
207 if isempty(res)
208 error('### The inputs were not properly handled. This shouldn''t happen.');
209 end
210
211 else
212 % we recursively pass back to this method
213 res = copy(args{1}, 1);
214 for kk=2:numel(args)
215 res = feval(op, res, args{kk});
216 end
217 end
218
219 % Set output
220 varargout{1} = res;
221
222 %---------- nested functions
223
224 %-------------------------------------------------
225 % Check the two inputs have compatible data types
226 function dout = compatibleData(a1,a2)
227 %== Data types
228 if (isa(a1.data, 'fsdata') && isa(a2.data, 'tsdata')) || ...
229 isa(a2.data, 'fsdata') && isa(a1.data, 'tsdata')
230 error('### Can not %s time-series data to frequency-series data.', opname);
231 end
232 % check X units for all data types
233 if ~isa(a1.data, 'cdata') && ~isa(a2.data, 'cdata')
234 if ~isempty(a1.data.xunits.strs) && ~isempty(a2.data.xunits.strs)
235 if a1.data.xunits ~= a2.data.xunits
236 error('### X units should be equal for the %s operator', op);
237 end
238 end
239 end
240
241 % determine output data type
242 d1 = copy(a1.data,1);
243 d2 = copy(a2.data,1);
244
245 if isa(d1, 'data2D') && isa(d2, 'data2D')
246 if numel(d1.y) > 1
247 dout = d1;
248 elseif numel(d2.y) > 1
249 dout = d2;
250 else
251 dout = d1;
252 end
253 elseif isa(d1, 'data2D') && isa(d2, 'cdata')
254 dout = d1;
255 elseif isa(d1, 'cdata') && isa(d2, 'data2D')
256 dout = d2;
257 else
258 dout = d1;
259 end
260
261 end
262
263 function uo = getYunits(a1, a2)
264 % For other operators we need to apply the operator
265 uo = feval(op, a1.data.yunits, a2.data.yunits);
266 end
267
268 % Perform the desired operation on the data
269 function y = operate(a1, a2)
270 y = feval(op, a1.data.y, a2.data.y);
271 end
272
273 % Perform the desired operation on the data uncertainty
274 function dy = operateError(a1, a2)
275
276 if ~isempty(a1.dy) || ~isempty(a2.dy)
277
278 da1 = a1.dy;
279 da2 = a2.dy;
280
281 if isempty(da1)
282 da1 = zeros(size(a1.y));
283 end
284 if isempty(da2)
285 da2 = zeros(size(a2.y));
286 end
287
288 switch op
289 case {'plus', 'minus'}
290 dy = sqrt(da1.^2 + da2.^2);
291 case {'times', 'mtimes'}
292 dy = sqrt( (da1./a1.y).^2 + (da2./a2.y).^2 ) .* abs(a1.y.*a2.y);
293 case {'rdivide', 'mrdivide'}
294 dy = sqrt( (da1./a1.y).^2 + (da2./a2.y).^2 ) .* abs(a1.y./a2.y);
295 otherwise
296 dy = [];
297 end
298
299 else
300 dy = [];
301 end
302
303 end
304
305 %-----------------------------------------------
306 % Get two new AO names from the input var names,
307 % the input AO names, and the indices.
308 function names = getNames(in_names, a1, jj, a2, kk)
309
310 % First variable name
311 if isempty(a1.name) && ~isempty(in_names{1})
312 if ~isempty(jj)
313 if numel(jj) == 1
314 names{1} = sprintf('%s(%d)', in_names{1}, jj);
315 else
316 names{1} = sprintf('%s(%d,%d)', in_names{1}, jj(1), jj(2));
317 end
318 else
319 names{1} = in_names{1};
320 end
321 else
322 if ~isempty(jj)
323 if numel(jj) == 1
324 % names{1} = sprintf('%s(%d)', a1.name, jj);
325 names{1} = sprintf('%s', a1.name);
326 else
327 % names{1} = sprintf('%s(%d,%d)', a1.name, jj(1), jj(2));
328 names{1} = sprintf('%s', a1.name);
329 end
330 else
331 names{1} = a1.name;
332 end
333 end
334 % Second variable name
335 if isempty(a2.name) && ~isempty(in_names{2})
336 if isempty(in_names{2})
337 in_names{2} = a2.name;
338 end
339 if ~isempty(kk)
340 if numel(kk) == 1
341 names{2} = sprintf('%s(%d)', in_names{2}, kk);
342 else
343 names{1} = sprintf('%s(%d,%d)', in_names{2}, kk(1), kk(2));
344 end
345 else
346 names{2} = in_names{2};
347 end
348 else
349 names{2} = a2.name;
350 if ~isempty(kk)
351 if numel(kk) == 1
352 % names{2} = sprintf('%s(%d)', a2.name, kk);
353 names{2} = sprintf('%s', a2.name);
354 else
355 % names{2} = sprintf('%s(%d,%d)', a2.name, kk(1), kk(2));
356 names{2} = sprintf('%s', a2.name);
357 end
358 else
359 names{2} = a2.name;
360 end
361 end
362
363 % The output AO name
364 names{3} = sprintf('(%s%s%s)', names{1}, opsym, names{2});
365 end
366
367 %-------------------------------------
368 % Return true if the input is a matrix
369 function r = ismatrix(a)
370 if nrows(a) > 1 && ncols(a) > 1
371 r = true;
372 else
373 r = false;
374 end
375 end
376
377 %-------------------------------------
378 % Return true if the input is a vector
379 function r = isvector(a)
380 if (nrows(a)==1 && ncols(a)>1) || (ncols(a)==1 && nrows(a)>1)
381 r = true;
382 else
383 r = false;
384 end
385 end
386
387 %-------------------------------------
388 % Return numnber of rows in the array
389 function r = nrows(a)
390 r = size(a,1);
391 end
392
393 %-------------------------------------
394 % Return numnber of cols in the array
395 function r = ncols(a)
396 r = size(a,2);
397 end
398
399
400 end % End of add
401
402
403
404 % END