Mercurial > hg > ltpda

comparison m-toolbox/classes/@pest/pest.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 a71a40911c27
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % PEST constructor for parameter estimates (pest) class.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: PEST constructor for parameter estimates (pest) class.
5 %
6 % CONSTRUCTOR:
7 %
8 % pe = pest()
9 % pe = pest(filename)
10 % pe = pest(plist-object)
11 % pe = pest(y, paramNames)
12 % pe = pest(y, paramNames, dy)
13 % pe = pest(y, paramNames, dy, cov)
14 % pe = pest(y, paramNames, dy, cov, chi2)
15 % pe = pest(y, paramNames, dy, cov, chi2, dof)
16 %
17 % INPUTS:
18 %
19 % y - best fit parameters
20 % paramNames - names of the parameters
21 % dy - standard errors of the parameters
22 % cov - covariance matrix of the parameters
23 % chi2 - reduced chi^2 of the final fit
24 % dof - degrees of freedom
25 %
26 % <a href="matlab:utils.helper.displayMethodInfo('pest', 'pest')">Parameters Description</a>
27 %
28 % VERSION: $Id: pest.m,v 1.33 2011/08/22 05:22:23 hewitson Exp $
29 %
30 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
31
32 classdef pest < ltpda_uoh
33
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35 % Property definition %
36 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
37
38 %---------- Protected read-only Properties ----------
39 properties (SetAccess = protected, Dependent = true)
40 end
41
42 %---------- Protected read-only Properties ----------
43 properties (GetAccess = public, SetAccess = protected)
44 dy % standard errors of the parameters.
45 y = []; % best fit parameters
46 names = {}; % names of the parameters, if any
47 yunits = unit.initObjectWithSize(1,0); % the units of each parameter
48 pdf = []; % posterior probability distribution of the parameters
49 cov = []; % covariance matrix of the parameters
50 corr = []; % correlation matrix of the parameters
51 chi2 = []; % reduced chi^2 of the final fit
52 dof = []; % degrees of freedom
53 chain = []; % monte carlo markov chain
54 models = []; % models that were fit
55 end
56
57 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
58 % Check property setting %
59 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
60
61 methods
62 % y
63 function set.y(obj, val)
64 try
65 % The genericSet method sets the value with a cell-array
66 val = [val{:}];
67 end
68 if ~isnumeric(val) || (~isvector(val) && ~isempty(val))
69 error('### The value for the property ''y'' must be a numeric vector');
70 end
71 if ~isempty(val)
72 obj.y = reshape(val, [], 1);
73 else
74 obj.y = val;
75 end
76 end
77 % dy
78 function set.dy(obj, val)
79 try
80 % The genericSet method sets the value with a cell-array
81 val = [val{:}];
82 end
83 if ~isnumeric(val) || (~isvector(val) && ~isempty(val))
84 error('### The value for the property ''dy'' must be a numeric vector');
85 end
86 if ~isempty(val)
87 obj.dy = reshape(val, [], 1);
88 else
89 obj.dy = val;
90 end
91 end
92 % names
93 function set.names(obj, val)
94 if numel(val) == 1 && iscell(val{1})
95 val = val{1};
96 end
97 if isempty(val)
98 val = {};
99 end
100 if ~(iscell(val) || ischar(val))
101 error('### The value for the property ''names'' must be a cell of parameter names.');
102 end
103 if iscell(val)
104 obj.names = val;
105 else
106 obj.names = cellstr(val);
107 end
108 end
109 % yunits
110 function set.yunits(obj, val)
111 if numel(val) == 1 && iscell(val) && iscell(val{1})
112 val = val{1};
113 end
114 if isempty(val)
115 obj.yunits = unit.initObjectWithSize(size(val,1), size(val,2));
116 elseif ischar(val)
117 obj.yunits = unit(val);
118 elseif isa(val, 'unit')
119 obj.yunits = copy(val,1);
120 elseif iscell(val)
121 obj.yunits = unit(val{:});
122 else
123 error('### The yunits value must be a vector of unit-objects or a string');
124 end
125 end
126 % pdf
127 function set.pdf(obj, val)
128 try
129 % The genericSet method sets the value with a cell-array
130 val = [val{:}];
131 end
132 if ~isnumeric(val) || ndims(val) ~= 2
133 error('### The value for the property ''pdf'' must be a numeric matrix');
134 end
135 obj.pdf = val;
136 end
137 % cov
138 function set.cov(obj, val)
139 try
140 % The genericSet method sets the value with a cell-array
141 val = [val{:}];
142 end
143 if ~isnumeric(val) || ndims(val) ~= 2
144 error('### The value for the property ''Cov'' must be a numeric matrix');
145 end
146 obj.cov = val;
147 end
148 % corr
149 function set.corr(obj, val)
150 try
151 % The genericSet method sets the value with a cell-array
152 val = [val{:}];
153 end
154 if ~isnumeric(val) || ndims(val) ~= 2
155 error('### The value for the property ''Corr'' must be a numeric matrix');
156 end
157 obj.corr = val;
158 end
159 % chi2
160 function set.chi2(obj, val)
161 try
162 % The genericSet method sets the value with a cell-array
163 val = [val{:}];
164 end
165 if ~isnumeric(val) || (any(size(val) ~= [1 1]) && ~isempty(val))
166 error('### The value for the property ''chi2'' must be a single double');
167 end
168 obj.chi2 = val;
169 end
170 % dof
171 function set.dof(obj, val)
172 try
173 % The genericSet method sets the value with a cell-array
174 val = [val{:}];
175 end
176 if ~isnumeric(val) || (any(size(val) ~= [1 1]) && ~isempty(val))
177 error('### The value for the property ''dof'' must be a single double');
178 end
179 obj.dof = val;
180 end
181 % dof
182 function set.chain(obj, val)
183 try
184 % The genericSet method sets the value with a cell-array
185 val = [val{:}];
186 end
187 if ~isnumeric(val)
188 error('### The value for the property ''chain'' must be a double');
189 end
190 obj.chain = val;
191 end
192
193 % models
194 function set.models(obj, val)
195 try
196 % The genericSet method sets the value with a cell-array
197 val = [val{:}];
198 end
199 if ~(isa(val, 'ltpda_uoh') || isempty(val))
200 error('### The value for the property ''models'' must be ltpda_uoh object');
201 end
202 if ~isempty(val)
203 obj.models = val;
204 else
205 obj.models = [];
206 end
207 end
208 end
209
210 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
211 % Compute the Dependent properties %
212 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
213
214 methods
215 % function val = get.dy(obj)
216 % val = sprintf('Please define a method which gets the dy from the object: %s', obj.name);
217 % end
218 end
219
220 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
221 % Constructor %
222 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
223
224 methods
225 function obj = pest(varargin)
226
227 import utils.const.*
228 utils.helper.msg(msg.OMNAME, 'running %s/%s', mfilename('class'), mfilename);
229
230 % Collect all pest objects
231 [objs, dummy, rest] = utils.helper.collect_objects(varargin(:), 'pest');
232
233 if isempty(rest) && ~isempty(objs)
234 % Do copy constructor and return
235 utils.helper.msg(msg.OPROC1, 'copy constructor');
236 obj = copy(objs, 1);
237 for kk=1:numel(obj)
238 obj(kk).addHistory(pest.getInfo('pest', 'None'), [], [], obj(kk).hist);
239 end
240 return
241 end
242
243 switch nargin
244 case 0
245 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
246 %%%%%%%%%%%%%%%%%%%%%%%%%%%% no input %%%%%%%%%%%%%%%%%%%%%%%%%%%*
247 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
248 utils.helper.msg(msg.OPROC1, 'empty constructor');
249 obj.addHistory(pest.getInfo('pest', 'None'), plist(), [], []);
250
251 case 1
252 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
253 %%%%%%%%%%%%%%%%%%%%%%%%%%% One input %%%%%%%%%%%%%%%%%%%%%%%%%%%
254 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
255
256 if ischar(varargin{1})
257 %%%%%%%%%% obj = pest('foo.mat') %%%%%%%%%%
258 %%%%%%%%%% obj = pest('foo.xml') %%%%%%%%%%
259 utils.helper.msg(msg.OPROC1, 'constructing from file %s', varargin{1});
260 obj = fromFile(obj, varargin{1});
261
262 elseif isstruct(varargin{1})
263 %%%%%%%%%% obj = pest(struct) %%%%%%%%%%
264 utils.helper.msg(msg.OPROC1, 'constructing from struct');
265 obj = obj.fromStruct(varargin{1});
266
267 elseif isa(varargin{1}, 'plist')
268 %%%%%%%%%% obj = pest(plist-object) %%%%%%%%%%
269
270 pl = varargin{1};
271
272 if pl.isparam('filename')
273 utils.helper.msg(msg.OPROC1, 'constructing from file %s', varargin{1});
274 obj = fromFile(obj, pl);
275
276 elseif pl.isparam('hostname') || pl.isparam('conn')
277 utils.helper.msg(msg.OPROC1, 'constructing from repository %s', pl.find('hostname'));
278 obj = obj.fromRepository(pl);
279
280 elseif pl.isparam('built-in')
281 utils.helper.msg(msg.OPROC1, 'constructing from built-in model');
282 obj = fromModel(obj, pl);
283
284 elseif pl.isparam('y')
285 utils.helper.msg(msg.OPROC1, 'constructing from values');
286 obj = obj.fromValues(pl);
287
288 else
289 obj.setProperties(pl);
290 obj.addHistory(pest.getInfo('pest', 'None'), pl, [], []);
291 end
292
293 elseif isnumeric(varargin{1})
294 %%%%%%%%%% obj = pest(y) %%%%%%%%%%
295 utils.helper.msg(msg.OPROC1, 'constructing from y');
296 obj.y = varargin{1};
297 pl = plist('y', varargin{1});
298 obj.fromValues(pl);
299
300 else
301 error('### Unknown single argument constructor.');
302 end
303
304 case 2
305 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
306 %%%%%%%%%%%%%%%%%%%%%%%%%%% two input %%%%%%%%%%%%%%%%%%%%%%%%%%%
307 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
308 if (isa(varargin{1}, 'database') || isa(varargin{1}, 'mpipeline.repository.RepositoryConnection')) && isnumeric(varargin{2})
309 %%%%%%%%%% obj = pest(<database-object>, [IDs]) %%%%%%%%%%
310 obj = obj.fromRepository(plist('conn', varargin{1}, 'id', varargin{2}));
311
312 elseif isnumeric(varargin{1}) && (iscell(varargin{2}) || ischar(varargin{2}))
313 %%%%%%%%%% obj = pest(y, paramNames) %%%%%%%%%%
314 utils.helper.msg(msg.OPROC1, 'constructing from y and param names');
315 pl = plist('y', varargin{1}, 'paramNames', varargin{2});
316 obj.fromValues(pl);
317
318 elseif isa(varargin{1}, 'pest') && isa(varargin{2}, 'plist') && isempty(varargin{2}.params)
319 %%%%%%%%%% obj = pest(pest, <empty-plist>) %%%%%%%%%%
320 obj = pest(varargin{1});
321
322 elseif isa(varargin{1}, 'org.apache.xerces.dom.DeferredElementImpl') && ...
323 isa(varargin{2}, 'history')
324 %%%%%%%%%% obj = pest(DOM node, history-objects) %%%%%%%%%%
325 obj = fromDom(obj, varargin{1}, varargin{2});
326
327 elseif isa(varargin{1}, 'ltpda_uoh') && isa(varargin{2}, 'plist')
328 %%%%%%%%%%% obj = pest(<ltpda_uoh>-object, plist-object) %%%%%%%%%%
329 % always recreate from plist
330
331 % If we are trying to load from file, and the file exists, do
332 % that. Otherwise, copy the input object.
333 if varargin{2}.isparam('filename')
334 if exist(fullfile('.', find(varargin{2}, 'filename')), 'file')==2
335 obj = pest(varargin{2});
336 else
337 obj = pest(varargin{1});
338 end
339 else
340 obj = pest(varargin{2});
341 end
342
343 else
344 error('### Unknown 2 argument constructor.');
345 end
346
347 case 3
348 if isnumeric(varargin{1}) && ...
349 (iscell(varargin{2}) || ischar(varargin{2})) && ...
350 isnumeric(varargin{3})
351 %%%%%%%%%% obj = pest(y, paramNames, dy) %%%%%%%%%%
352 utils.helper.msg(msg.OPROC1, 'constructing from y, param names and dy');
353 pl = plist(...
354 'y', varargin{1}, ...
355 'paramNames', varargin{2}, ...
356 'dy', varargin{3});
357 obj.fromValues(pl);
358 else
359 error('### Unknown 3 argument constructor.');
360 end
361
362 case 4
363 if isnumeric(varargin{1}) && ...
364 (iscell(varargin{2}) || ischar(varargin{2})) && ...
365 isnumeric(varargin{3}) && ...
366 isnumeric(varargin{4})
367 %%%%%%%%%% obj = pest(y, paramNames, dy) %%%%%%%%%%
368 utils.helper.msg(msg.OPROC1, 'constructing from y, param names, dy and cov');
369 pl = plist(...
370 'y', varargin{1}, ...
371 'paramNames', varargin{2}, ...
372 'dy', varargin{3}, ...
373 'cov', varargin{4});
374 obj.fromValues(pl);
375 else
376 error('### Unknown 4 argument constructor.');
377 end
378
379 case 5
380 if isnumeric(varargin{1}) && ...
381 (iscell(varargin{2}) || ischar(varargin{2})) && ...
382 isnumeric(varargin{3}) && ...
383 isnumeric(varargin{4}) && ...
384 isnumeric(varargin{5})
385 %%%%%%%%%% obj = pest(y, paramNames, dy) %%%%%%%%%%
386 utils.helper.msg(msg.OPROC1, 'constructing from y, param names, dy, cov and chi2');
387 pl = plist(...
388 'y', varargin{1}, ...
389 'paramNames', varargin{2}, ...
390 'dy', varargin{3}, ...
391 'cov', varargin{4}, ...
392 'chi2', varargin{5});
393 obj.fromValues(pl);
394 else
395 error('### Unknown 5 argument constructor.');
396 end
397
398 case 6
399 if isnumeric(varargin{1}) && ...
400 (iscell(varargin{2}) || ischar(varargin{2})) && ...
401 isnumeric(varargin{3}) && ...
402 isnumeric(varargin{4}) && ...
403 isnumeric(varargin{5}) && ...
404 isnumeric(varargin{6})
405 %%%%%%%%%% obj = pest(y, paramNames, dy) %%%%%%%%%%
406 utils.helper.msg(msg.OPROC1, 'constructing from y, param names, dy, cov, chi2 and dof');
407 pl = plist(...
408 'y', varargin{1}, ...
409 'paramNames', varargin{2}, ...
410 'dy', varargin{3}, ...
411 'cov', varargin{4}, ...
412 'chi2', varargin{5}, ...
413 'dof', varargin{6});
414 obj.fromValues(pl);
415 else
416 error('### Unknown 6 argument constructor.');
417 end
418
419 otherwise
420 error('### Unknown number of arguments.');
421 end
422
423 end % End constructor
424 end % End public methods
425
426 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
427 % Methods (Public, hidden) %
428 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
429
430 methods (Hidden = true)
431 varargout = attachToDom(varargin)
432 end
433
434 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
435 % Methods (protected) %
436 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
437
438 methods (Access = protected)
439 varargout = fromStruct(varargin)
440 varargout = fromDom(varargin)
441 end
442
443 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
444 % Methods (private) %
445 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
446
447 methods (Access = private)
448 % Constructors
449 varargout = fromValues(varargin)
450 varargout = genericSet(varargin)
451 end
452
453 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
454 % Methods (static) %
455 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
456 methods (Static)
457
458 function mdls = getBuiltInModels(varargin)
459 mdls = ltpda_uo.getBuiltInModels('pest');
460 end
461
462 function out = VEROUT()
463 out = '$Id: pest.m,v 1.33 2011/08/22 05:22:23 hewitson Exp $';
464 end
465
466 function ii = getInfo(varargin)
467 ii = utils.helper.generic_getInfo(varargin{:}, 'pest');
468 end
469
470 function out = SETS()
471 out = [SETS@ltpda_uoh, {'From Values'}];
472 end
473
474
475 function plout = getDefaultPlist(set)
476 persistent pl;
477 persistent lastset;
478 if exist('pl', 'var')==0 || isempty(pl) || ~strcmp(lastset, set)
479 pl = pest.buildplist(set);
480 lastset = set;
481 end
482 plout = pl;
483 end
484
485 function out = buildplist(set)
486
487 if ~utils.helper.ismember(lower(pest.SETS), lower(set))
488 error('### Unknown set [%s]', set);
489 end
490
491 out = plist();
492 out = pest.addGlobalKeys(out);
493 out = buildplist@ltpda_uoh(out, set);
494
495 switch lower(set)
496 case 'from values'
497 % y
498 p = param({'y','Best fit parameters'}, paramValue.EMPTY_DOUBLE);
499 out.append(p);
500 % parameter names
501 p = param({'paramNames','Names of the parameters'}, paramValue.EMPTY_CELL);
502 out.append(p);
503 % dy
504 p = param({'dy','Standard errors of the parameters'}, paramValue.EMPTY_DOUBLE);
505 out.append(p);
506 % cov
507 p = param({'cov','Covariance matrix of the parameters'}, paramValue.EMPTY_DOUBLE);
508 out.append(p);
509 % chi2
510 p = param({'chi2','Reduced chi^2 of the final fit'}, paramValue.EMPTY_DOUBLE);
511 out.append(p);
512 % dof
513 p = param({'dof','Degrees of freedom'}, paramValue.EMPTY_DOUBLE);
514 out.append(p);
515 % models
516 p = param({'models',['Models used for the fit in which the coefficients were calculated<br>' ...
517 'Please notice that the models need to be stored in a <tt>smodel</tt> object!']}, smodel);
518 out.append(p);
519 % pdf
520 p = param({'pdf','Probability density function, as output by MCMC methods'}, paramValue.EMPTY_DOUBLE);
521 out.append(p);
522 % yunits
523 p = param({'yunits','A cell-array of units for the parameters.'}, paramValue.EMPTY_CELL);
524 out.append(p);
525
526 end
527 end % function out = getDefaultPlist(varargin)
528
529 function obj = initObjectWithSize(n,m)
530 obj = pest.newarray([n m]);
531 end
532
533 end % End static methods
534
535 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
536 % Methods (static, private) %
537 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
538
539 methods (Static, Access = private)
540 end % End static, private methods
541
542 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
543 % Methods (static, hidden) %
544 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
545
546 methods (Static = true, Hidden = true)
547 varargout = loadobj(varargin)
548 varargout = update_struct(varargin);
549 end
550
551 end