1 classdef Assessment < handle
2 % a
class used to compute reduced several reduced simulations over a huge
3 % parameter sample extracting useful information
5 % The
class generates one or two dimensional cell arrays of these data fields
6 % storing information
for a variation of some attributes of the reduced
7 % model, which can be freely chosen by specification of #plot_fields .
9 properties(Access=
public)
10 % is a vector of field names with at most 2 elements over which the error
11 % landscape is computed or a cell with empty entries.
13 % A reasonable choice would be
'{ \'N\', \'M\' }'. The error landscape is
14 % plotted over a variation of these fields specified by samples.
15 % In
case of an empty cell entry,
'N' and
'M' are coupled by the fixed
16 % ratio #M_by_N_ratio and samples should be a vector of coupling constant
17 % between 0 and 1. See also: update_rmodel()
22 % is a cell array of scalar row vectors. The scalar values are
'plot_field'
23 % values
for which the error is computed and tested.
25 % @
default =
'{[0.1:0.1:1]}'
26 samples = {(0.1:0.1:1)};
28 % an parameter sampling
object of type ParameterSampling.Interface
30 % @
default is a random sampling with 10 elements and seed 654321;
32 % See also: init_random_sampling()
35 % is a cell array of description texts for the plot fields. If this field
36 % is not set it is set to 'plot_fields'.
37 plot_field_descr = {};
40 % character
string specifying
this test run. The
string should be unique
41 %
for every parameter combination.
43 % @
default rmodel.name + '_stoachastic_assessment';
46 %
boolean value determining whether error estimates shall be computed
47 compute_estimates =
false;
49 %
boolean value determining whether the
"true" error `\|u_h(\mu) -
50 % u_{\text{red}}(\mu)\|_{{\cal W}_h}` shall be computed
51 compute_errors =
false;
53 %
boolean value determining whether the condition numbers of the system
54 % matrices shall be computed
55 compute_conditions =
false;
57 %
boolean value indicating whether the refinement steps of the reduced
58 % basis generation by Greedy.TrainingSetAdaptation shall be taken into account
60 % This option is only useful
if
61 % -# #M_test is set to the training sample of the parameter space
63 % -# the #plot_fields include the basis dimensions
64 % In this case the parameter sampling used during basis generation for the
65 % given basis dimension is selected.
66 follow_refinement_steps = false;
68 %
if M and N are couple
this value specifies the constant
for the coupling.
70 % A value of zero means `M_{\max} / N_{\max}`
83 % the
"real" collection of vectors used
for the #plot_fields variation
88 % the underlying detailed data of type Greedy.DataTree.Detailed.INode
90 % @todo make
this cached...
96 function dd =
get.detailed_data(
this)
97 dd = this.cache_object.
obj.detailed_data;
100 function rsamples =
get.rsamples(
this)
101 if length(this.samples) == 2
102 rsamples = combvec(this.samples{1},this.samples{2});
104 rsamples = this.samples{1};
108 function sta =
Assessment(matfile, rmodel, mu_set_size, seed)
109 %
function sta =
Assessment(matfile, rmodel, mu_set_size, seed)
113 % matfile: name of result file where a
IReducedModel object and a
114 % Greedy.DataTree.Detailed.INode
object must be stored.
115 % mu_set_size: Optional argument determining the number of random
116 % parameters in the validation sample. (Default = 10)
117 % seed: a random seed used
for initialization of
118 % ParameterSampling.Random
object for the #M_test parameter
119 % sample
set. (Default = 654321)
121 if nargin == 1 && isa(matfile,
'Postprocess.StochasticAssessment.Assessment')
124 fns = m.Properties(cellfun(@(x) ~isequal(x.SetAccess, 'none') && ~x.Dependent, m.Properties));
125 fns = cellfun(@(x) x.Name, fns, 'UniformOutput', false);
126 fns = setdiff(fns, {
'cache_object'});
127 for i = 1:length(fns);
128 sta.(fns{i}) = cp.(fns{i});
130 sta.cache_object =
CacheableObject(cp.cache_object.matfile,
'detailed_data');
131 sta.rmodel = copy(cp.rmodel);
133 sta.rmodel = copy(rmodel);
141 init_random_sampling(sta, mu_set_size, seed);
143 if isempty(sta.plot_field_descr)
144 sta.plot_field_descr = sta.plot_fields;
146 assert(length(sta.plot_field_descr) == length(sta.plot_fields));
149 sta.run_name = [rmodel.descr.name, '_stochastic_assessment'];
153 function output=compute(this)
154 % function compute(this)
155 % stochastic estimation of error between reduced and detailed simulation over a
156 % test sample of `\mu` vectors
158 % This function stochastically estimates the error between reduced and detailed
159 % simulations for given `\mu`-vectors and various reduced and collateral basis
160 % sizes. The results are visualized in a surface plot for problems with CRB.
162 % If required by the users, averaged time measurements for the reduced and the
163 % detailed simulations are computed, too.
166 % output: an
Output object
168 %if isequal(params.mode, 'error_ei')
169 % model.detailed_simulation = model.detailed_ei_simulation;
172 %if isequal(params.mode, 'error_ei_rb_proj')
173 % model.detailed_simulation = model.detailed_ei_rb_proj_simulation;
176 reduced_data = gen_reduced_data(this.rmodel, this.detailed_data);
178 if this.compute_estimates
179 assert(this.rmodel.Mstrich > 0);
180 assert(this.rmodel.enable_error_estimator);
182 assert(~this.rmodel.enable_error_estimator);
183 assert(this.rmodel.Mstrich == 0);
186 basis_gen = this.detailed_data.bg_algorithm;
187 while ~isa(basis_gen, 'Greedy.Algorithm')
188 switch(class(basis_gen))
189 case 'Greedy.Combined'
190 basis_gen = get_rb_basis_generator(basis_gen);
191 case 'Greedy.TrainingSetAdaptation'
192 basis_gen = basis_gen.child;
197 generator = basis_gen.detailed_gen;
199 if ~generator.is_valid
200 generator = SnapshotsGenerator.Trajectories(this.rmodel.detailed_model);
203 % default value (only need for mode == 'check' when no 'mu_set_size' is given)
205 savefile = [this.run_name];
206 tmpfile = [savefile,'_tmp.mat'];
207 savefile = [savefile,'.mat'];
209 rps_size = length(this.rsamples);
211 complete_out = cell(1, rps_size);
213 if this.compute_conditions
214 rd_conds = cell(1, rps_size);
219 for comb = 1:rps_size
221 pf_val=this.rsamples(:,comb);
223 for i=1:length(pf_val)
224 update_rmodel(this, this.plot_fields{i}, pf_val(i));
226 temp_rd = extract_reduced_data_subset(this.rmodel, reduced_data);
228 if this.compute_conditions
229 rd_conds{comb} = get_conds(temp_rd);
232 disp([
'Processing combination ', num2str(comb),
'/', num2str(rps_size)]);
234 M = this.get_test_space_sample;
235 generator.prepare(this.rmodel.detailed_model,
this.detailed_data.model_data, M);
238 part_out =
struct(
'max_Delta' , cell(1, mu_size), ...
239 'max_Delta_ind', cell(1, mu_size), ...
240 'max_err' , cell(1, mu_size), ...
241 'max_err_ind' , cell(1, mu_size), ...
242 'rtime' , cell(1, mu_size), ...
243 'rrtime' , cell(1, mu_size), ...
244 'dtime' , cell(1, mu_size), ...
245 'rconds' , cell(1, mu_size), ...
246 'dconds' , cell(1, mu_size));
247 %
for mu_ind = 1:mu_size
248 parfor mu_ind = 1:mu_size
251 tdm = trm.detailed_model;
255 disp([
'processing parameter vector ',num2str(mu_ind),...
256 '/',num2str(size(M,1))]);
258 set_mu(trm, M(mu_ind,:));
259 set_mu(tdm, M(mu_ind,:));
261 compute_out = compute_impl(tthis, tgen, trm, trd);
263 part_out(mu_ind) = compute_out;
266 complete_out{comb} = part_out;
267 % disp([
'saving temporary workspace in ', tmpfile]);
268 % rmodel = this.rmodel;
269 % detailed_data = this.detailed_data;
270 % save(fullfile(rbmatlabtemp,tmpfile),
'complete_out',
'part_out',
'this',
'rmodel',
'detailed_data',
'reduced_data');
274 if length(this.samples) == 2
275 ns1 = length(this.samples{1});
276 ns2 = rps_size / ns1;
277 complete_out = reshape(complete_out, ns1, ns2);
278 if this.compute_conditions
279 rd_conds = reshape(rd_conds, ns1, ns2);
283 output = Postprocess.StochasticAssessment.Output(complete_out,
this);
284 output.rd_conds = rd_conds;
286 save(fullfile(rbmatlabresult, savefile));
290 function init_random_sampling(
this, mu_set_size, seed)
291 %
function init_random_sampling(
this, mu_set_size, seed)
292 % generates the initial random sample
set #M_test
298 detailed_model = this.rmodel.detailed_model;
299 this.M_test = ParameterSampling.Random(mu_set_size, seed);
300 if this.M_test.init_required()
301 init_sample(this.M_test, detailed_model);
305 function other = copy(
this)
306 %
function other = copy(
this)
307 % deep copy of
this object
311 other = Postprocess.StochasticAssessment.
Assessment(
this);
314 function equal_distribution_samples(
this, min, max, sample_size)
315 %
function equal_distribution_samples(min, max, sample_size)
316 % adapts the #samples attribute such that they are equally distributed in
320 % min: is a vector of minimum values of
'plot_fields' variables.
321 % max: is a vector
for maximum values of
'plot_fields'
323 % sample_size: specifies the number of sample values between
'min' and
324 %
'max' for which the error is computed and plotted. The
325 %
default value is
'max-min'.
328 sample_size = max-min;
332 assert(all(max > min));
334 this.samples = cell(1,length(max));
335 for i = 1:length(max)
336 this.samples{i} = min(i) + ...
337 round((0:sample_size(i)-1) .* ...
338 (max(i)-min(i))./(sample_size(i)-1));
342 function update_rmodel(
this, plot_field, new_value)
343 %
function update_rmodel(
this, plot_field, new_value)
344 % updates the reduced model, i.e. applies the changes defined by the
345 % #plot_fields attribute
348 % plot_field: an entry in #plot_fields,
if it is empty, a coupling of
349 %
'N' and
'M' is assured.
350 % new_value: the
new value it should be
set to.
352 if ~isempty(plot_field)
353 this.rmodel.(plot_field) = new_value;
358 rmodel = this.rmodel;
359 disp(['N = ', num2str(rmodel.N), ' M = ', num2str(rmodel.M), ' Mstrich = ', num2str(rmodel.Mstrich)]);
368 function [rb_sim_data, tictoc] = rb_simulation_tictoc(varargin)
369 % function [rb_sim_data, tictoc] = rb_simulation_tictoc(varargin)
370 % a wrapper around the
IReducedModel.rb_simulation() measuring the
374 % rb_sim_data: the reduced simulation data
375 % tictoc: execution time
376 [rb_sim_data, tictoc] = tictoc_wrapper(@rb_simulation, varargin{:});
382 methods(Access=
private)
384 function output = compute_impl(this, generator, rmodel, reduced_data)
388 model_data = this.detailed_data.model_data;
399 rmodel.compute_conditions = this.compute_conditions;
401 [rb_sim_data,tictoc] = this.rb_simulation_tictoc(rmodel, reduced_data);
404 if this.compute_estimates
405 Delta = rb_sim_data.Delta;
406 [max_Delta, max_Delta_ind] = max(Delta);
411 if this.compute_conditions
412 rconds = rb_sim_data.conds;
414 % condition computations may make time measurements bogus.
418 if this.compute_errors
419 dmodel = rmodel.detailed_model;
420 if this.compute_conditions
422 dmodel.compute_conditions =
true;
426 [U_H, dtic_toc, opts] = generate(generator, dmodel, model_data, fields);
428 if this.compute_conditions
431 % condition computations may make time measurements bogus.
436 rb_sim_data = rb_reconstruction(rmodel, this.detailed_data.datatree, rb_sim_data);
440 errors = dmodel.l2_error_sequence_algorithm(U_H, rb_sim_data.U, ...
442 [max_err,max_err_ind] = max(errors);
448 warning([
'catched an error: ', ME.message]);
449 disp(
'setting error to NaN');
452 output.max_Delta = max_Delta;
453 output.max_Delta_ind = max_Delta_ind;
454 output.max_err = max_err;
455 output.max_err_ind = max_err_ind;
456 output.rtime = rtime;
457 output.rrtime = rrtime;
458 output.dtime = dtime;
459 output.rconds = rconds;
460 output.dconds = dconds;
466 % modifies the reduced basis size fields of
'model' by a single variable.
469 % c: sets the reduced basis sizes to
470 % `(N,M) = c (N_{\mbox{max}}, c_{MbyN} N_{\mbox{max}})`. For time-adaptive
471 % schemes with several collateral reduced basis spaces, the following
472 % formula is applied: `(N,M^k) = c (N_{\mbox{max}}), c_{MbyN}
473 % \frac{M^k_{\mbox{max}}}{\mbox{max}_{k=0,...,K} M^k_{\mbox{max}}}
476 dd = get_by_description(this.detailed_data.datatree,
'rb');
477 Nmax = get_rb_size(dd);
478 Mmax = get_ei_size(dd);
479 maxMmax = max(Mmax)-this.rmodel.Mstrich;
480 this.rmodel.N = max(round(c*Nmax),1);
481 if this.M_by_N_ratio == 0
482 M_by_N_r = maxMmax / Nmax;
483 elseif this.M_by_N_ratio == -1
484 leaf_dd = get_active_leaf(this.detailed_data, this.rmodel);
485 assert(isa(leaf_dd,
'Greedy.DataTree.Detailed.PODEI'));
486 ext_ind = this.rmodel.N;
487 ext_counts = get_field(leaf_dd,
'ei_ext_count');
488 ei_ext_steps = get_field(leaf_dd,
'ei_ext_steps');
489 rb_real_ext_steps = get_field(leaf_dd,
'rb_real_ext_steps');
490 if isempty(ext_counts)
491 n_ei_exts = length(ei_ext_steps);
492 ext_counts = ones(size(ei_ext_steps))*maxMmax/n_ei_exts;
494 Mcomp = round(sum((ei_ext_steps <= rb_real_ext_steps(ext_ind)).*ext_counts));
495 this.rmodel.M = min(maxMmax, max(Mcomp, 1));
498 M_by_N_r = this.M_by_N_ratio;
500 this.rmodel.M = arrayfun(@(x) min(x, max(round(c*x/maxMmax*M_by_N_r*Nmax),1)), maxMmax);
503 function M = get_test_space_sample(this)
504 if this.follow_refinement_steps
505 info = active_info(this.detailed_data);
506 ref_level = get_ref_level(info, this.rmodel);
507 M = get_sample_at_ref_level(this.M_test, ref_level);
509 M = this.M_test.sample;