1 classdef EI < Greedy.Plugin.EICommon
2 % Plugin
for the Greedy.Algorithm
class generating a collateral reduced basis
3 % space plus interpolation DOFs and a local grid.
5 % This can be used
for the empirical interpolation of parametrized functions
6 % or
operator evaluations.
8 properties (Hidden=
true)
9 % temporary matrix generated by prepare_reduced_data()
14 generated_basis_type = 'ei';
17 properties (Constant, Hidden=true)
19 % cell array of field names that shall be copied to the generated
20 % Greedy.DataTree.Detailed.EILeafNode instance.
21 info_fields = {
'stop_Mmax',
'ei_target_error',
'noof_extensions', ...
22 'minimum_residual',
'enable_regularisation' };
26 % maximum number of generated collateral reduced basis vectors.
29 % error mode
for the empirical interpolation error
31 % This can be either
'interpol' or
'approx' for the
true interpolation
32 % error or the projection error of the
operator interpolation error in the
33 % collateral reduced basis space.
34 ei_target_error =
'approx';
36 %
boolean flag indicating that the error indicator which is assumed to be
37 % minimized is not equal to the ones returned by error_indicators()
39 % This is
set to
'true' by the Greedy.Plugin.PODEI
class.
40 target_error_external =
false;
42 % specifies how many collateral reduced basis functions shall be added in
43 % each extension step.
45 % This makes only sense to be
set greater than
'1' for trajectories, i.e.
46 % when more than one possible candidates
for basis extension exist
for one
50 %
if the residual is below
this barrier, the basis extension step is skipped.
52 % This option is only used in
case of #target_error_external .
53 minimum_residual = 1e-6;
55 %
if this regularisation flag is
set to
true, the collateral reduced basis
56 % spaces are forced to have zero mean, enforcing global conservation in
57 %
case of
operator interpolation.
58 enable_regularisation =
false;
62 function eide =
EI(generator)
63 % constructor storing the snapshot generator
66 % generator: snapshot generator of type SnapshotsGenerator.Cached.
67 eide = eide@Greedy.Plugin.
EICommon(generator);
70 function set.ei_target_error(
this, value)
71 %
if isequal(value,
'linfty-interpol')
72 % this.use_l2_error = 0;
74 % this.use_l2_error = 1;
76 this.ei_target_error = value;
79 function detailed_data = init_basis(
this, rmodel, model_data, M_train)
80 %
function detailed_data = init_basis(
this, rmodel, model_data, M_train)
81 % @copybrief .Greedy.Plugin.Interface.init_basis()
86 % rmodel: of type .Greedy.User.IReducedModel
89 % detailed_data: of type .Greedy.DataTree.Detailed.EILeafNode
91 % This method initializes the detailed data node.
92 if isa(model_data,
'Greedy.DataTree.Detailed.INode')
93 detailed_data = model_data;
94 model_data = detailed_data.model_data;
96 detailed_data = Greedy.DataTree.Detailed.EILeafNode(this, model_data);
97 set_fields(detailed_data, this, Greedy.Plugin.
EI.info_fields);
99 ophandle = this.generator.ophandle;
100 W = ophandle.inner_product_matrix(model_data);
102 % determine initial size of vectors:
103 detailed_data.QM = zeros(size(W,1),0);
104 detailed_data.TM = zeros(0,1);
105 detailed_data.BM = zeros(0,0);
109 function prepare_reduced_data(this, rmodel, detailed_data)
110 % function prepare_reduced_data(this, rmodel, detailed_data)
111 % @copybrief Greedy.Plugin::
Interfaceprepare_reduced_data()
113 % @copydetails Greedy.Plugin::
Interfaceprepare_reduced_data()
117 % detailed_data: of type .Greedy.DataTree.Detailed.EILeafNode
119 % This method pre-computes the normal solution matrix for
120 % `L^2`-approximation in the collateral reduced space, in case
121 % 'ei_target_error' is set to 'approx'.
122 ophandle = this.generator.ophandle;
123 if isequal(get_field(detailed_data, 'ei_target_error'), 'approx')
124 W = ophandle.inner_product_matrix(detailed_data.model_data);
125 QM = detailed_data.QM;
126 this.A = (QM' * W * QM)^(-1) * QM' * W;
132 function Uapprox = generate_reduced(this, rmodel, reduced_data, detailed_data, U)
133 % function Uapprox = generate_reduced(this, [], [], detailed_data, U)
134 % @copybrief Greedy.Plugin::
Interfacegenerate_reduced()
136 % @copydetails Greedy.Plugin::
Interfacegenerate_reduced()
143 switch get_field(detailed_data, 'ei_target_error')
146 Uapprox = detailed_data.QM*P;
147 case {
'interpol',
'linfty-interpol'}
148 coeff = detailed_data.BM \ U(detailed_data.TM,:);
149 Uapprox = detailed_data.QM * coeff;
151 error(
'ei_target_error unknown!');
155 function [breakloop, reason] = pre_check_for_end(
this, rmodel, detailed_data)
156 %
function [breakloop, reason] = pre_check_for_end(
this, rmodel, detailed_data)
162 % rmodel: of type .Greedy.User.IReducedModel
163 % detailed_data: of type .Greedy.DataTree.Detailed.EILeafNode
166 dd_leaf = get_active_leaf(detailed_data, rmodel);
168 if length(detailed_data.TM) >= get_field(dd_leaf,
'stop_Mmax')
169 set_stop_flag(dd_leaf, 'stopped_on_max_basis_size');
171 reason = '
EI stopped on maximum basis size';
172 elseif stopped_on_any_child(dd_leaf, 'stopped_on_empty_extension')
174 reason = '
EI stopped on empty extension';
175 elseif stopped_on_all_leafs(dd_leaf, 'stopped_on_epsilon')
177 reason = '
EI stopped on epsilon';
179 max_errs = get_field(dd_leaf, 'max_err_sequence');
180 % check for monotonicity, but not in case of external error estimation,
181 % because this does not make any sense then.
182 if length(max_errs) > 1 && (max_errs(end)>max_errs(end-1)) ...
183 && isequal(get_field(dd_leaf, 'ei_target_error'), 'approx') ...
184 && ~this.target_error_external
185 set_stop_flag(dd_leaf, 'stopped_on_monotonicity');
187 reason = '
EI stopped on monotonicity';
192 function detailed_data = basis_extension(this, rmodel, detailed_data, max_err_seq, mu)
193 % function detailed_data = basis_extension(this, rmodel, detailed_data, max_err_seq, mu)
194 % @copybrief Greedy.Plugin::
Interfacebasis_extension()
196 % @copydetails Greedy.Plugin::
Interfacebasis_extension()
200 % detailed_data: of type .Greedy.DataTree.Detailed.EILeafNode
203 % detailed_data: of type .Greedy.DataTree.Detailed.EILeafNode
205 dmodel = rmodel.detailed_model;
207 LU = generate(this.generator, dmodel, detailed_data.model_data);
209 for ext_num = 1:get_field(detailed_data, 'noof_extensions')
212 prepare_reduced_data(this, rmodel, detailed_data);
213 max_err_seq = compute_error(this, rmodel, [], detailed_data);
216 [dummy, max_i] = max(max_err_seq);
217 zeta_M = LU(:,max_i(1));
222 QM = detailed_data.QM;
223 TM = detailed_data.TM;
224 BM = detailed_data.BM;
227 %% compute interpolant
230 % get target values in interpolation points
231 zeta_M_part = zeta_M(TM);
233 % get interpolation coefficients
234 %sigma = bicgstab(BM,zeta_M_part,[],1000);
235 sigma = BM \ zeta_M_part;
238 r_M = zeta_M - QM * sigma;
240 disp(['positive minus negative contributions of residual (should be zero:)', num2str(sum(r_M))]);
241 disp(['positive minus negative contributions of zeta_M (should be zero:)', num2str(sum(zeta_M))]);
243 % if(sum(r_M) > 1e-4)
247 max_r = max(abs(r_M));
248 if (this.target_error_external || ext_num >=2)...
249 && max_r < get_field(detailed_data, 'minimum_residual') * 10^(ext_num-1)
251 disp(['skipping basis extension after ', num2str(ext_num), 'extension steps...']);
252 ext_num = ext_num - 1;
255 t_M = find(abs(r_M) >= max_r - 1e-12, 1);
258 if ~isempty(find(t_M(1)==TM, 1))
259 % if target error is external, this error is not immediately fatal!
260 if this.target_error_external
261 disp('skipping EIDetailed basis extension because a magic point is already selected');
262 disp('This might happen because of very bad parameter choice for ei extension.');
263 disp(['Error at twice selected point is: ', num2str(max_r)]);
264 disp('The error should be very low!');
265 ext_num = ext_num - 1;
268 error('empirical interpolation selected magic point twice!!');
271 q_M = r_M / r_M(t_M(1)); % then sign is respected
short version
272 % (after verifying zero-1 structure...)
273 % BM = [BM, zeros(m-1,1) ; ...
274 % QM(t_M,:), 1 ]; % interpolation matrix mit BM(i,j) = q_j(t_i)
275 % TM = [TM; t_M]; % new vector of interpolation DOFs
276 % QM = [QM, q_M]; % new matrix of basis functions
278 disp(['positive minus negative contributions of new basis function (should be zero:)', num2str(sum(q_M))]);
280 if abs(sum(q_M)) > 1e-12 && get_field(detailed_data, 'enable_regularisation')
281 subset = abs(q_M) > 1e-8;
282 if length(subset) > 10
283 q_M(subset) = q_M(subset) - sum(q_M)/length(q_M(subset));
286 disp('WARNING: applied regularisation!');
287 disp(['positive minus negative contributions of new basis function (should be zero:)', num2str(sum(q_M))]);
292 BM = [BM, q_M(TM) ; ...
293 QM(t_M,:), q_M(t_M) ]; % interpolation matrix mit BM(i,j) = q_j(t_i)
294 TM = [TM; t_M]; % new vector of interpolation DOFs
295 QM = [QM, q_M]; % new matrix of basis functions
297 detailed_data.BM = BM;
298 detailed_data.TM = TM;
299 detailed_data.QM = QM;
303 set_stop_flag(detailed_data, 'stopped_on_empty_extension');
309 append_field(detailed_data, 'ei_ext_count', ext_num);
314 function cop = copy(this)
315 % function cop = copy(this)
316 % deep copy of Greedy.Plugin.
EI object
319 % cop: copied
object of type Greedy.Plugin.
EI
320 cop =
EI(this.generator);
322 fnames = fieldnames(this);
323 for i=1:length(fnames)
324 cop.(fnames{i}) = this.(fnames{i});