Mercurial > hg > ltpda

comparison m-toolbox/classes/@pest/eval.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 % EVAL evaluate a pest object
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: EVAL evaluate a pest model.
5 %
6 % CALL: b = eval(p, pl)
7 % b = eval(p, x, pl)
8 % b = eval(p, x1, ... , xN, pl)
9 % b = eval(p, [x1 ... xN], pl)
10 %
11 % INPUTS: p - input pest(s) containing parameter values.
12 % xi - input ao(s) containing x values (as x or y fields, depending on the 'xfield' parameter)
13 % pl - parameter list (see below)
14 %
15 % OUTPUTs: b - an AO containing the model evaluated at the given X
16 % values, with the given parameter values.
17 %
18 % <a href="matlab:utils.helper.displayMethodInfo('pest', 'eval')">Parameters Description</a>
19 %
20 % VERSION: $Id: eval.m,v 1.28 2011/05/15 22:47:14 mauro Exp $
21 %
22 % EXAMPLES:
23 %
24 % % 1)
25 % % Prepare the symbolic model
26 % mdl = smodel(plist('expression', 'a1.*x + a2.*x.^2 + a0', 'xvar', 'x', 'yunits', 'V'));
27 %
28 % % Prepare the pest object
29 %
30 % p = pest(plist('paramnames', {'a0','a1','a2'}, 'y', [1 2 3], 'models', mdl));
31 %
32 % % Evaluate the object
33 % a1 = eval(p, plist('xdata', ao([1:10])))
34 % a2 = eval(p, ao([1:10]))
35 %
36 % % 2)
37 % % Prepare the symbolic model
38 % mdl = smodel(plist('expression', 'a1.*x1 + a2.*x2 + a0', 'xvar', {'x1', 'x2'}, 'yunits', 'm', 'xunits', {'T', 'K'}));
39 %
40 % % Prepare the pest object
41 %
42 % p = pest(plist('paramnames', {'a0','a1','a2'}, 'y', [1 2 3], 'yunits', {'m', 'T/m', 'K/m'}, 'models', mdl));
43 %
44 % % Evaluate the object
45 % x1 = ao(plist('yvals', [1:10], 'fs', 1, 'yunits', 'T'));
46 % x2 = ao(plist('yvals', [1:10], 'fs', 1, 'yunits', 'K'));
47 % a1 = eval(p, plist('xdata', [x1 x2]))
48 % a2 = eval(p, [x1 x2])
49 % a3 = eval(p, x1, x2)
50 %
51 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
52
53 function varargout = eval(varargin)
54
55 % Check if this is a call for parameters
56 if utils.helper.isinfocall(varargin{:})
57 varargout{1} = getInfo(varargin{3});
58 return
59 end
60
61 import utils.const.*
62 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
63
64 % Collect input variable names
65 in_names = cell(size(varargin));
66 for ii = 1:nargin,in_names{ii} = inputname(ii);end
67
68 % Collect all AOs and plists
69 [psts, pst_invars, rest] = utils.helper.collect_objects(varargin(:), 'pest', in_names);
70 [pl, pl_invars, rest] = utils.helper.collect_objects(rest, 'plist', in_names);
71 [as, as_invars, rest] = utils.helper.collect_objects(rest, 'ao', in_names);
72 [c, c_invars] = utils.helper.collect_objects(rest, 'cell', in_names);
73
74 if nargout == 0
75 error('### eval can not be used as a modifier method. Please give at least one output');
76 end
77
78 % combine plists
79 pl = parse(pl, getDefaultPlist());
80
81 % Extract necessary parameters
82 index = find(pl, 'index');
83 x = find(pl, 'Xdata');
84 if ~isempty(as)
85 x = as;
86 end
87 % I don't know how to deal with the history of a cell array of aos
88 if ~isempty(c)
89 error('Please pass the arguments in a vector or a list')
90 end
91
92 % Extract the information about the x field, if necessary
93 xfield = pl.find('xfield');
94 switch xfield
95 case 'x'
96 op = str2func('x');
97 case 'y'
98 op = str2func('y');
99 otherwise
100 error('oops')
101 end
102
103 % Extract the xvals for the smodel, and the output x for the ao
104 switch class(x)
105 case 'cell'
106 switch class(x{1})
107 case 'ao'
108 data_type = find(pl, 'type');
109 if isempty(data_type)
110 % Nothing to do, output data type will be inherited from the first AO
111 if ~isa(x{1}.data, 'cdata')
112 out_x = x{1};
113 else
114 out_x = [];
115 end
116 out_xunits = '';
117 else
118 % In this case I have to extract the vector with the output x and xunits
119 out_x = feval(op, x{1});
120 out_xunits = feval(str2func([xfield 'units']), x{1});
121 end
122 case 'double'
123 data_type = find(pl, 'type');
124 xvals = x;
125 switch data_type
126 case {'tsdata', 'fsdata', 'xydata'}
127 out_x = x;
128 out_xunits = find(pl, 'xunits');
129 case {'cdata', ''}
130 out_x = [];
131 out_xunits = '';
132 otherwise
133 error('LTPDA error')
134 end
135 otherwise
136 end
137 case 'ao'
138 data_type = find(pl, 'type');
139 if isempty(data_type)
140 % Nothing to do, output data type will be inherited from the first AO
141 if ~isa(x(1).data, 'cdata')
142 out_x = x(1);
143 else
144 out_x = [];
145 end
146 out_xunits = '';
147 else
148 % In this case I have to extract the vector with the output x and xunits
149 out_x = feval(op, x(1));
150 out_xunits = feval(str2func([xfield 'units']), x(1));
151 end
152 case 'double'
153 data_type = find(pl, 'type');
154 xvals = x;
155 switch data_type
156 case {'tsdata', 'fsdata', 'xydata'}
157 out_x = x;
158 case {'cdata', ''}
159 out_x = [];
160 otherwise
161 error('LTPDA error')
162 end
163 out_xunits = find(pl, 'xunits');
164 otherwise
165 end
166
167 % If the user wants to override the pest/smodel yunits, let's get them
168 % This works only if the user sets them to something not empty
169 yunits = find(pl, 'yunits');
170
171 % If we have AOs in a cell....
172 if iscell(x) && all(cellfun(@(x)isa(x, 'ao'), x))
173 % if we have multiple x, we need to convert the y values into a cell array
174 xvals = cellfun(op, x, 'UniformOutput', false);
175 elseif isa(x, 'ao')
176 % we put the y values in to a cell array
177 if numel(x) > 1
178 xvals = cellfun(op, num2cell(x), 'UniformOutput', false);
179 else
180 % ... we take the x values, as per the help
181 xvals = feval(op, x);
182 end
183 end
184
185 % Loop over input objects
186 for jj = 1:numel(psts)
187
188 pst = psts(jj);
189 % evaluate models
190 m = copy(pst.models, true);
191 switch class(m)
192 case 'smodel'
193 % Make sure the smodel parameters are named the same as the pest
194 m(index).setParams(pst.names, pst.y);
195 % If the user provided the x vector(s), override the smodel x with these
196 if ~isempty(xvals)
197 m(index).setXvals(xvals);
198 end
199 % Go for the model evaluation
200 out(jj) = eval(m(index), plist(...
201 'output type', data_type, 'output x', out_x, 'output xunits', out_xunits));
202
203 % Setting the units of the evaluated model
204 if ~isempty(yunits)
205 out(jj).setYunits(yunits);
206 end
207
208 case 'ao'
209 % do linear combination: using lincom
210 out(jj) = lincom(m, pst);
211 out(jj).simplifyYunits;
212
213
214 case 'matrix'
215 % check objects of the matrix and switch
216 switch class(m(1).objs)
217 case 'smodel'
218 % Make sure the smodel parameters are named the same as the pest
219 for ii = 1:numel(m.objs)
220 m.objs(ii).setParams(pst.names, pst.y);
221 end
222 % If the user provided the x vector(s), override the smodel x with these
223 if ~isempty(x)
224 for ii = 1:numel(m.objs)
225 m.objs(ii).setXvals(x);
226 end
227 end
228 % Go for the model evaluation
229 tout = ao.initObjectWithSize(size(m.objs,1),size(m.objs,2));
230 for ii=1:size(m.objs,1)
231 for kk=1:size(m.objs,2)
232 tout(ii,kk) = eval(m.objs(ii,kk), plist(...
233 'output type', data_type, 'output x', out_x));
234 % Setting the units of the evaluated model
235 if ~isempty(yunits)
236 tout(ii,kk).setYunits(yunits);
237 end
238 end
239 end
240 out(jj) = matrix(tout);
241 case 'ao'
242 % get params from the pest object
243 prms = pst.y;
244 % build cdata aos
245 prmsao = ao.initObjectWithSize(numel(prms),1);
246 for ii = 1:numel(prms)
247 prmsao(ii) = ao(cdata(prms(ii)));
248 prmsao(ii).setYunits(pst.yunits(ii));
249 end
250 % build matrix for parameters
251 prm = matrix(prmsao);
252 % build matrix for the model
253 mm = ao.initObjectWithSize(numel(m(1).objs),numel(prms));
254 for ii = 1:numel(m(1).objs)
255 for kk = 1:numel(prms)
256 mm(ii,kk) = m(kk).getObjectAtIndex(ii);
257 end
258 end
259 mmat = matrix(mm);
260 % eval model
261 tout = mmat*prm;
262 out(jj) = tout;
263 end
264 otherwise
265 error('### current version of pest/eval needs the ''models'' field to be a smodel')
266 end
267
268
269
270 % uncertainties for the evaluated model: calculate them from covariance matrix
271 if ~isempty(pst.cov) && utils.prog.yes2true(pl.find('errors'));
272 switch class(m)
273 case 'smodel'
274 C = pst.cov;
275 p = pst.names;
276 % here we need a matrix of "functions" which are the derivatives wrt parameters,
277 % evaluated at each point x:
278
279 F = [];
280 for kk = 1:length(p)
281 md = eval(diff(m(index), plist('var', p{kk})));
282 F = [F md.y];
283 end
284 % The formula is:
285 % D = F * C * F';
286 % and then we need to take
287 % dy = sqrt(diag(D))
288 if size(md.y, 1) > 1
289 % Make sure we work with columns
290 out(jj).setDy(sqrt(sum((F * C)' .* F'))');
291 else
292 out(jj).setDy(sqrt(sum((F' * C)' .* F)));
293 end
294
295 otherwise
296 warning('Propagation of the errors on the model not yet implemented')
297 end
298 end
299
300 % Set output AO name
301 name = sprintf('eval(%s,', pst.name);
302 for kk = 2:numel(pst)
303 name = [name pst(kk).name ','];
304 end
305 name = [name(1:end-1) ')'];
306 out(jj).name = name;
307 % Add history
308 if isempty(as)
309 out(jj).addHistory(getInfo('None'), pl, pst_invars, pst(:).hist);
310 else
311 out(jj).addHistory(getInfo('None'), pl, {pst_invars as_invars}, [pst(:).hist as(:).hist]);
312 end
313 end
314
315 % Set output
316 varargout = utils.helper.setoutputs(nargout, out);
317
318 end
319
320
321 %--------------------------------------------------------------------------
322 % Get Info Object
323 %--------------------------------------------------------------------------
324 function ii = getInfo(varargin)
325 if nargin == 1 && strcmpi(varargin{1}, 'None')
326 sets = {};
327 pl = [];
328 else
329 sets = {'Default'};
330 pl = getDefaultPlist();
331 end
332 % Build info object
333 ii = minfo(mfilename, 'pest', 'ltpda', utils.const.categories.sigproc, '$Id: eval.m,v 1.28 2011/05/15 22:47:14 mauro Exp $', sets, pl);
334 end
335
336 %--------------------------------------------------------------------------
337 % Get Default Plist
338 %--------------------------------------------------------------------------
339 function plout = getDefaultPlist()
340 persistent pl;
341 if ~exist('pl', 'var') || isempty(pl)
342 pl = buildplist();
343 end
344 plout = pl;
345 end
346
347 function pl = buildplist()
348
349 pl = plist();
350
351 % INDEX
352 p = param({'index', 'Select which model must be evaluated if more than one.'}, paramValue.DOUBLE_VALUE(1));
353 pl.append(p);
354
355 % XDATA
356 p = param({'Xdata', ['The X values to evaluate the model at. This can be:<ul>'...
357 '<li>a double vector </li>' ...
358 '<li>a cell array of double vectors</li>' ...
359 '<li>a single AO (from which the Y data will be extracted)</li>' ...
360 '<li>a cell array of AOs (from which the Y data will be extracted)</li></ul>' ...
361 ]}, paramValue.EMPTY_DOUBLE);
362 pl.append(p);
363
364 % XFIELD
365 p = param({'xfield', 'Choose the field to extract the x values from when inputting AOs for parameter ''xdata''.'}, {2, {'x', 'y'}, paramValue.SINGLE});
366 pl.append(p);
367
368 % TYPE
369 pv = paramValue.DATA_TYPES;
370 % Add an 'empty' on top of the list
371 pv{2} = [{''} pv{2}];
372 p = param({'type', ['Choose the data type for the output ao.<br>'...
373 'If empty, and if the user input AOs as ''XDATA'', the type will be inherited.']}, pv);
374 p.val.setValIndex(1);
375 pl.append(p);
376
377 % YUNITS
378 p = param({'yunits','Unit on Y axis.'}, paramValue.STRING_VALUE(''));
379 pl.append(p);
380
381 % XUNITS
382 p = param({'xunits','Unit on X axis.'}, paramValue.STRING_VALUE(''));
383 pl.append(p);
384
385 % ERRORS
386 p = param({'errors', ['Estimate the uncertainty of the output values based <br>' ...
387 'on the parameters covariance matrix']}, paramValue.TRUE_FALSE);
388 pl.append(p);
389
390 end
391 % END