1 /* _______________________________________________________________________
2
3 DAKOTA: Design Analysis Kit for Optimization and Terascale Applications
4 Copyright 2014-2020 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
5 This software is distributed under the GNU Lesser General Public License.
6 For more information, see the README file in the top Dakota directory.
7 _______________________________________________________________________ */
8
9 //- Class: NonDMultilevelSampling
10 //- Description: Implementation code for NonDMultilevelSampling class
11 //- Owner: Mike Eldred
12 //- Checked by:
13 //- Version:
14
15 #include "dakota_system_defs.hpp"
16 #include "dakota_data_io.hpp"
17 #include "dakota_tabular_io.hpp"
18 #include "DakotaModel.hpp"
19 #include "DakotaResponse.hpp"
20 #include "NonDMultilevelSampling.hpp"
21 #include "ProblemDescDB.hpp"
22 #include "DiscrepancyCalculator.hpp"
23 #include "DakotaIterator.hpp"
24
25 #ifdef HAVE_NPSOL
26 #include "NPSOLOptimizer.hpp"
27 #elif HAVE_OPTPP
28 #include "SNLLOptimizer.hpp"
29 #endif
30
31 static const char rcsId[]="@(#) $Id: NonDMultilevelSampling.cpp 7035 2010-10-22 21:45:39Z mseldre $";
32
33
34 namespace Dakota {
35
36 /** This constructor is called for a standard letter-envelope iterator
37 instantiation. In this case, set_db_list_nodes has been called and
38 probDescDB can be queried for settings from the method specification. */
39 NonDMultilevelSampling::
NonDMultilevelSampling(ProblemDescDB & problem_db,Model & model)40 NonDMultilevelSampling(ProblemDescDB& problem_db, Model& model):
41 NonDSampling(problem_db, model),
42 pilotSamples(problem_db.get_sza("method.nond.pilot_samples")),
43 randomSeedSeqSpec(problem_db.get_sza("method.random_seed_sequence")),
44 mlmfIter(0),
45 allocationTarget(problem_db.get_short("method.nond.allocation_target")),
46 qoiAggregation(problem_db.get_short("method.nond.qoi_aggregation")),
47 useTargetVarianceOptimizationFlag(problem_db.get_bool("method.nond.allocation_target.variance.optimization")),
48 finalCVRefinement(true),
49 exportSampleSets(problem_db.get_bool("method.nond.export_sample_sequence")),
50 exportSamplesFormat(
51 problem_db.get_ushort("method.nond.export_samples_format"))
52 {
53 // initialize scalars from sequence
54 seedSpec = randomSeed = random_seed(0);
55
56 // Support multilevel LHS as a specification override. The estimator variance
57 // is known/correct for MC and an assumption/approximation for LHS. To get an
58 // accurate LHS estimator variance, one would need:
59 // (a) assumptions about separability -> analytic variance reduction by a
60 // constant factor
61 // (b) similarly, assumptions about the form relative to MC (e.g., a constant
62 // factor largely cancels out within the relative sample allocation.)
63 // (c) numerically-generated estimator variance (from, e.g., replicated LHS)
64 if (!sampleType) // SUBMETHOD_DEFAULT
65 sampleType = SUBMETHOD_RANDOM;
66
67 // check iteratedModel for model form hierarchy and/or discretization levels;
68 // set initial response mode for set_communicators() (precedes core_run()).
69 if (iteratedModel.surrogate_type() == "hierarchical")
70 aggregated_models_mode();
71 else {
72 Cerr << "Error: Multilevel Monte Carlo requires a hierarchical "
73 << "surrogate model specification." << std::endl;
74 abort_handler(METHOD_ERROR);
75 }
76
77 ModelList& ordered_models = iteratedModel.subordinate_models(false);
78 size_t i, j, num_mf = ordered_models.size(), num_lev,
79 prev_lev = std::numeric_limits<size_t>::max(),
80 pilot_size = pilotSamples.size();
81 ModelLRevIter ml_rit; bool err_flag = false;
82 NLev.resize(num_mf);
83 for (i=num_mf-1, ml_rit=ordered_models.rbegin();
84 ml_rit!=ordered_models.rend(); --i, ++ml_rit) { // high fid to low fid
85 // for now, only SimulationModel supports solution_{levels,costs}()
86 num_lev = ml_rit->solution_levels(); // lower bound is 1 soln level
87
88 if (num_lev > prev_lev) {
89 Cerr << "\nWarning: unused solution levels in multilevel sampling for "
90 << "model " << ml_rit->model_id() << ".\n Ignoring "
91 << num_lev - prev_lev << " of " << num_lev << " levels."<< std::endl;
92 num_lev = prev_lev;
93 }
94
95 // Ensure there is consistent cost data available as SimulationModel must
96 // be allowed to have empty solnCntlCostMap (when optional solution control
97 // is not specified). Passing false bypasses lower bound of 1.
98 if (num_lev > ml_rit->solution_levels(false)) { // default is 0 soln costs
99 Cerr << "Error: insufficient cost data provided for multilevel sampling."
100 << "\n Please provide solution_level_cost estimates for model "
101 << ml_rit->model_id() << '.' << std::endl;
102 err_flag = true;
103 }
104
105 //Sizet2DArray& Nl_i = NLev[i];
106 NLev[i].resize(num_lev); //Nl_i.resize(num_lev);
107 //for (j=0; j<num_lev; ++j)
108 // Nl_i[j].resize(numFunctions); // defer to pre_run()
109
110 prev_lev = num_lev;
111 }
112 if (err_flag)
113 abort_handler(METHOD_ERROR);
114
115 switch (pilot_size) {
116 case 0: maxEvalConcurrency *= 100; break;
117 case 1: maxEvalConcurrency *= pilotSamples[0]; break;
118 default: {
119 size_t max_ps = 0;
120 for (i=0; i<pilot_size; ++i)
121 if (pilotSamples[i] > max_ps)
122 max_ps = pilotSamples[i];
123 if (max_ps)
124 maxEvalConcurrency *= max_ps;
125 break;
126 }
127 }
128
129 // For testing multilevel_mc_Qsum():
130 //subIteratorFlag = true;
131 }
132
133
~NonDMultilevelSampling()134 NonDMultilevelSampling::~NonDMultilevelSampling()
135 { }
136
137
resize()138 bool NonDMultilevelSampling::resize()
139 {
140 bool parent_reinit_comms = NonDSampling::resize();
141
142 Cerr << "\nError: Resizing is not yet supported in method "
143 << method_enum_to_string(methodName) << "." << std::endl;
144 abort_handler(METHOD_ERROR);
145
146 return parent_reinit_comms;
147 }
148
149
pre_run()150 void NonDMultilevelSampling::pre_run()
151 {
152 NonDSampling::pre_run();
153
154 // reset sample counters to 0
155 size_t i, j, num_mf = NLev.size(), num_lev;
156 for (i=0; i<num_mf; ++i) {
157 Sizet2DArray& Nl_i = NLev[i];
158 num_lev = Nl_i.size();
159 for (j=0; j<num_lev; ++j)
160 Nl_i[j].assign(numFunctions, 0);
161 }
162 }
163
164
165 /** The primary run function manages the general case: a hierarchy of model
166 forms (from the ordered model fidelities within a HierarchSurrModel),
167 each of which may contain multiple discretization levels. */
core_run()168 void NonDMultilevelSampling::core_run()
169 {
170 //model,
171 // surrogate hierarchical
172 // ordered_model_fidelities = 'LF' 'MF' 'HF'
173 //
174 // Future: include peer alternatives (1D list --> matrix)
175 // For MLMC, could seek adaptive selection of most correlated
176 // alternative (or a convex combination of alternatives).
177
178 // TO DO: hierarchy incl peers (not peers each optionally incl hierarchy)
179 // num_mf = iteratedModel.model_hierarchy_depth();
180 // num_peer_i = iteratedModel.model_peer_breadth(i);
181
182 // TO DO: this initial logic is limiting:
183 // > allow MLMC and CVMC for either model forms or discretization levels
184 // > separate method specs that both map to NonDMultilevelSampling ???
185
186 // TO DO: following pilot sample across levels and fidelities in mixed case,
187 // could pair models for CVMC based on estimation of rho2_LH.
188
189 // For two-model control variate methods, select lowest,highest fidelities
190 size_t num_mf = NLev.size();
191 unsigned short lf_form = 0, hf_form = num_mf - 1; // ordered_models = low:high
192 if (num_mf > 1) {
193 size_t num_hf_lev = NLev.back().size();
194 if (num_hf_lev > 1) { // ML performed on HF with CV using available LF
195 // multiple model forms + multiple solutions levels --> perform MLMC on
196 // HF model and bind 1:min(num_hf,num_lf) LF control variates starting
197 // at coarsest level (TO DO: validate case of unequal levels)
198 if (true) // reformulated approach using 1 new QoI correlation per level
199 multilevel_control_variate_mc_Qcorr(lf_form, hf_form);
200 else // original approach using 1 discrepancy correlation per level
201 multilevel_control_variate_mc_Ycorr(lf_form, hf_form);
202 }
203 else { // multiple model forms (only) --> CVMC
204 // use nominal value from user input, ignoring solution_level_control
205 UShortArray hf_lf_key; unsigned short lev = USHRT_MAX;
206 Pecos::DiscrepancyCalculator::
207 form_key(0, hf_form, lev, lf_form, lev, hf_lf_key);
208 control_variate_mc(hf_lf_key);
209 }
210 }
211 else { // multiple solutions levels (only) --> traditional ML-MC
212 if (true) //(subIteratorFlag)
213 multilevel_mc_Qsum(hf_form); // w/ error est, unbiased central moments
214 else
215 multilevel_mc_Ysum(hf_form); // lighter weight
216 }
217 }
218
219
220 /** This function performs "geometrical" MLMC on a single model form
221 with multiple discretization levels. */
multilevel_mc_Ysum(unsigned short model_form)222 void NonDMultilevelSampling::multilevel_mc_Ysum(unsigned short model_form)
223 {
224 // Formulate as a coordinated progression towards convergence, where, e.g.,
225 // time step is inferred from the spatial discretization (NOT an additional
226 // solution control) based on stability criteria, e.g. CFL condition. Can
227 // we reliably capture runtime estimates as part of pilot run w/i Dakota?
228 // Ultimately seems desirable to support either online or offline cost
229 // estimates, to allow more accurate resource allocation when possible
230 // or necessary (e.g., combustion processes with expense that is highly
231 // parameter dependent).
232 // model
233 // id_model = 'LF'
234 // simulation
235 // # point to state vars; ordered based on set values for h, delta-t
236 // solution_level_control = 'dssiv1'
237 // # relative cost estimates in same order as state set values
238 // # --> re-sort into map keyed by increasing cost
239 // solution_level_cost = 10 2 200
240
241 // How to manage the set of MLMC statistics:
242 // 1. Simplest: proposal is to use the mean estimator to drive the algorithm,
243 // but carry along other estimates.
244 // 2. Later: could consider a refinement for converging the estimator of the
245 // variance after convergence of the mean estimator.
246
247 // How to manage the vector of QoI:
248 // 1. Worst case: select N_l based only on QoI w/ highest total variance
249 // from pilot run --> fix for all levels and don't allow switching
250 // across major iterations (possible oscillation? Or simple overlay
251 // of resolution reqmts?)
252 // 2. Better: select N_l based on convergence in aggregated variance.
253
254 // assign truth model form (solution level assignment is deferred until loop)
255 UShortArray truth_key;
256 unsigned short seq_index = 2, lev = USHRT_MAX; // lev updated in loop below
257 Pecos::DiscrepancyCalculator::form_key(0, model_form, lev, truth_key);
258 iteratedModel.active_model_key(truth_key);
259 Model& truth_model = iteratedModel.truth_model();
260
261 size_t qoi, num_steps = truth_model.solution_levels();// 1 model form
262 unsigned short& step = (true) ? lev : model_form; // option not active
263 size_t max_iter = (maxIterations < 0) ? 25 : maxIterations; // default = -1
264 Real eps_sq_div_2, sum_sqrt_var_cost, estimator_var0 = 0., lev_cost;
265 // retrieve cost estimates across soln levels for a particular model form
266 RealVector cost = truth_model.solution_level_costs(), agg_var(num_steps);
267 // For moment estimation, we accumulate telescoping sums for Q^i using
268 // discrepancies Yi = Q^i_{lev} - Q^i_{lev-1} (sum_Y[i] for i=1:4).
269 // For computing N_l from estimator variance, we accumulate square of Y1
270 // estimator (YY[i] = (Y^i)^2 for i=1).
271 IntRealMatrixMap sum_Y; RealMatrix sum_YY(numFunctions, num_steps);
272 initialize_ml_Ysums(sum_Y, num_steps);
273
274 // Initialize for pilot sample
275 SizetArray delta_N_l;
276 load_pilot_sample(pilotSamples, NLev, delta_N_l);
277
278 // raw eval counts are accumulation of allSamples irrespective of resp faults
279 SizetArray raw_N_l(num_steps, 0);
280 RealVectorArray mu_hat(num_steps);
281 Sizet2DArray& N_l = NLev[model_form];
282
283 // now converge on sample counts per level (N_l)
284 mlmfIter = 0;
285 while (Pecos::l1_norm(delta_N_l) && mlmfIter <= max_iter) {
286
287 sum_sqrt_var_cost = 0.;
288 for (step=0; step<num_steps; ++step) { // step is reference to lev
289
290 configure_indices(step, model_form, lev, seq_index);
291 lev_cost = level_cost(cost, step);
292
293 // set the number of current samples from the defined increment
294 numSamples = delta_N_l[step];
295
296 // aggregate variances across QoI for estimating N_l (justification:
297 // for independent QoI, sum of QoI variances = variance of QoI sum)
298 Real& agg_var_l = agg_var[step]; // carried over from prev iter if no samp
299 if (numSamples) {
300
301 // advance any sequence specifications (seed_sequence)
302 assign_specification_sequence(step);
303 // generate new MC parameter sets
304 get_parameter_sets(iteratedModel);// pull dist params from any model
305
306 // export separate output files for each data set. truth_model()
307 // has the correct data when in bypass-surrogate mode.
308 if (exportSampleSets)
309 export_all_samples("ml_", iteratedModel.truth_model(), mlmfIter, step);
310
311 // compute allResponses from allVariables using hierarchical model
312 evaluate_parameter_sets(iteratedModel, true, false);
313
314 // process allResponses: accumulate new samples for each qoi and
315 // update number of successful samples for each QoI
316 //if (mlmfIter == 0) accumulate_offsets(mu_hat[lev]);
317 accumulate_ml_Ysums(sum_Y, sum_YY, lev, mu_hat[step], N_l[step]);
318 if (outputLevel == DEBUG_OUTPUT)
319 Cout << "Accumulated sums (Y1, Y2, Y3, Y4, Y1sq):\n" << sum_Y[1]
320 << sum_Y[2] << sum_Y[3] << sum_Y[4] << sum_YY << std::endl;
321 // update raw evaluation counts
322 raw_N_l[step] += numSamples;
323
324 // compute estimator variance from current sample accumulation:
325 if (outputLevel >= DEBUG_OUTPUT)
326 Cout << "variance of Y[" << step << "]: ";
327 agg_var_l
328 = aggregate_variance_Ysum(sum_Y[1][step], sum_YY[step], N_l[step]);
329 }
330
331 sum_sqrt_var_cost += std::sqrt(agg_var_l * lev_cost);
332 // MSE reference is MLMC with pilot sample, prior to any N_l adaptation:
333 if (mlmfIter == 0)
334 estimator_var0
335 += aggregate_mse_Ysum(sum_Y[1][step], sum_YY[step], N_l[step]);
336 }
337 // compute epsilon target based on relative tolerance: total MSE = eps^2
338 // which is equally apportioned (eps^2 / 2) among discretization MSE and
339 // estimator variance (\Sum var_Y_l / N_l). Since we do not know the
340 // discretization error, we compute an initial estimator variance from MLMC
341 // on the pilot sample and then seek to reduce it by a relative_factor <= 1.
342 if (mlmfIter == 0) { // eps^2 / 2 = var * relative factor
343 eps_sq_div_2 = estimator_var0 * convergenceTol;
344 if (outputLevel == DEBUG_OUTPUT)
345 Cout << "Epsilon squared target = " << eps_sq_div_2 << std::endl;
346 }
347
348 // update targets based on variance estimates
349 Real fact = sum_sqrt_var_cost / eps_sq_div_2, N_target;
350 for (step=0; step<num_steps; ++step) {
351 // Equation 3.9 in CTR Annual Research Briefs:
352 // "A multifidelity control variate approach for the multilevel Monte
353 // Carlo technique," Geraci, Eldred, Iaccarino, 2015.
354 N_target = std::sqrt(agg_var[step] / level_cost(cost, step)) * fact;
355 delta_N_l[step] = one_sided_delta(average(N_l[step]), N_target);
356 }
357 ++mlmfIter;
358 Cout << "\nMLMC iteration " << mlmfIter << " sample increments:\n"
359 << delta_N_l << std::endl;
360 }
361
362 // aggregate expected value of estimators for Y, Y^2, Y^3, Y^4. Final expected
363 // value is sum of expected values from telescopic sum. Note: raw moments
364 // have no bias correction (no additional variance from an estimated center).
365 RealMatrix Q_raw_mom(numFunctions, 4);
366 RealMatrix &sum_Y1 = sum_Y[1], &sum_Y2 = sum_Y[2],
367 &sum_Y3 = sum_Y[3], &sum_Y4 = sum_Y[4];
368 for (qoi=0; qoi<numFunctions; ++qoi) {
369 for (step=0; step<num_steps; ++step) {
370 size_t Nlq = N_l[lev][qoi];
371 Q_raw_mom(qoi,0) += sum_Y1(qoi,step) / Nlq;
372 Q_raw_mom(qoi,1) += sum_Y2(qoi,step) / Nlq;
373 Q_raw_mom(qoi,2) += sum_Y3(qoi,step) / Nlq;
374 Q_raw_mom(qoi,3) += sum_Y4(qoi,step) / Nlq;
375 }
376 }
377 // Convert uncentered raw moment estimates to final moments (central or std)
378 convert_moments(Q_raw_mom, momentStats);
379
380 // compute the equivalent number of HF evaluations (includes any sim faults)
381 equivHFEvals = raw_N_l[0] * cost[0]; // first level is single eval
382 for (step=1; step<num_steps; ++step) // subsequent levels incur 2 model costs
383 equivHFEvals += raw_N_l[step] * (cost[step] + cost[step-1]);
384 equivHFEvals /= cost[num_steps-1]; // normalize into equivalent HF evals
385 }
386
387
388 static RealVector *static_lev_cost_vec(NULL);
389 static size_t *static_qoi(NULL);
390 static Real *static_eps_sq_div_2(NULL);
391 static RealVector *static_Nlq_pilot(NULL);
392
393 static IntRealMatrixMap *static_sum_Ql(NULL);
394 static IntRealMatrixMap *static_sum_Qlm1(NULL);
395 static IntIntPairRealMatrixMap *static_sum_QlQlm1(NULL);
396
target_var_objective_eval_npsol(int & mode,int & n,double * x,double & f,double * gradf,int & nstate)397 void NonDMultilevelSampling::target_var_objective_eval_npsol(int &mode, int &n, double *x, double &f, double *gradf,
398 int &nstate) {
399 RealVector optpp_x;
400 RealVector optpp_grad_f;
401 optpp_x.size(n);
402 optpp_grad_f.size(n);
403
404 f = -1; //Dummy value
405
406 for (size_t i = 0; i < n; ++i) {
407 optpp_x[i] = x[i];
408 }
409 target_var_objective_eval_optpp(mode, n, optpp_x, f, optpp_grad_f, nstate);
410
411 for (size_t i = 0; i < n && mode; ++i) {
412 gradf[i] = optpp_grad_f[i];
413 }
414
415 }
416
417 void
target_var_constraint_eval_npsol(int & mode,int & m,int & n,int & ldJ,int * needc,double * x,double * g,double * grad_g,int & nstate)418 NonDMultilevelSampling::target_var_constraint_eval_npsol(int &mode, int &m, int &n, int &ldJ, int *needc, double *x,
419 double *g, double *grad_g, int &nstate) {
420 RealVector optpp_x;
421 RealVector optpp_g;
422 RealMatrix optpp_grad_g(1, n);
423 optpp_x.size(n);
424 optpp_g.size(n);
425
426 for (size_t i = 0; i < n; ++i) {
427 optpp_x[i] = x[i];
428 }
429
430 target_var_constraint_eval_optpp(mode, n, optpp_x, optpp_g, optpp_grad_g, nstate);
431
432 g[0] = optpp_g[0];
433 for (size_t i = 0; i < n && mode; ++i) {
434 grad_g[i] = optpp_grad_g[0][i];
435 }
436 }
437
438
target_var_constraint_eval_logscale_npsol(int & mode,int & m,int & n,int & ldJ,int * needc,double * x,double * g,double * grad_g,int & nstate)439 void NonDMultilevelSampling::target_var_constraint_eval_logscale_npsol(int& mode, int& m, int& n, int& ldJ, int* needc, double* x, double* g, double* grad_g, int& nstate){
440 RealVector optpp_x;
441 RealVector optpp_g;
442 RealMatrix optpp_grad_g(1, n);
443 optpp_x.size(n);
444 optpp_g.size(n);
445
446 for (size_t i = 0; i < n; ++i) {
447 optpp_x[i] = x[i];
448 }
449
450 target_var_constraint_eval_logscale_optpp(mode, n, optpp_x, optpp_g, optpp_grad_g, nstate);
451
452 g[0] = optpp_g[0];
453 for (size_t i = 0; i < n && mode; ++i) {
454 grad_g[i] = optpp_grad_g[0][i];
455 }
456 }
457
target_var_objective_eval_optpp(int mode,int n,const RealVector & x,double & f,RealVector & grad_f,int & result_mode)458 void NonDMultilevelSampling::target_var_objective_eval_optpp(int mode, int n, const RealVector &x, double &f,
459 RealVector &grad_f, int &result_mode) {
460 f = 0;
461
462 #ifdef HAVE_NPSOL
463 #elif HAVE_OPTPP
464 if(mode & OPTPP::NLPFunction){
465 result_mode = OPTPP::NLPFunction;
466 #endif
467 for (int i = 0; i < n; ++i) {
468 f += x[i] * (*static_lev_cost_vec)[i];
469 }
470 #ifdef HAVE_NPSOL
471 #elif HAVE_OPTPP
472 }
473 #endif
474
475 #ifdef HAVE_NPSOL
476 #elif HAVE_OPTPP
477 if(mode & OPTPP::NLPGradient){
478 result_mode = OPTPP::NLPGradient;
479 #endif
480 for (int i = 0; i < n; ++i) {
481 grad_f[i] = (*static_lev_cost_vec)[i];
482 }
483 #ifdef HAVE_NPSOL
484 #elif HAVE_OPTPP
485 }
486 #endif
487
488 }
489
target_var_constraint_eval_logscale_optpp(int mode,int n,const RealVector & x,RealVector & g,RealMatrix & grad_g,int & result_mode)490 void NonDMultilevelSampling::target_var_constraint_eval_logscale_optpp(int mode, int n, const RealVector& x, RealVector& g,
491 RealMatrix& grad_g, int& result_mode){
492 Real agg_estim_var = 0;
493 size_t num_lev = n;
494 target_var_constraint_eval_optpp(mode, n, x, g, grad_g, result_mode);
495 agg_estim_var = g[0];
496 #ifdef HAVE_NPSOL
497 #elif HAVE_OPTPP
498 if(mode & OPTPP::NLPFunction) {
499 #endif
500 g[0] = std::log(g[0]); // - (*static_eps_sq_div_2);
501 #ifdef HAVE_NPSOL
502 #elif HAVE_OPTPP
503 }
504 if(mode & OPTPP::NLPGradient) {
505 #endif
506 for (size_t lev = 0; lev < num_lev; ++lev) {
507 grad_g[0][lev] = grad_g[0][lev]/agg_estim_var;
508 }
509 #ifdef HAVE_NPSOL
510 #elif HAVE_OPTPP
511 }
512 #endif
513 }
514
515
target_var_constraint_eval_optpp(int mode,int n,const RealVector & x,RealVector & g,RealMatrix & grad_g,int & result_mode)516 void NonDMultilevelSampling::target_var_constraint_eval_optpp(int mode, int n, const RealVector &x, RealVector &g,
517 RealMatrix &grad_g, int &result_mode) {
518
519 bool compute_gradient = false;
520 #ifdef HAVE_NPSOL
521 compute_gradient = mode; //if mode == 0, NPSOL ignores gradients
522 #elif HAVE_OPTPP
523 if(mode & OPTPP::NLPFunction) {
524 result_mode = OPTPP::NLPFunction;
525 }
526 if(mode & OPTPP::NLPGradient){
527 compute_gradient = true;
528 result_mode = OPTPP::NLPGradient;
529 }
530 #endif
531
532 // std error in variance or std deviation estimate
533 size_t lev = 0;
534 Real Nlq = x[lev];
535 size_t Nlq_pilot = (*static_Nlq_pilot)[lev];
536 size_t qoi = *static_qoi;
537 size_t num_lev = n;
538 RealVector agg_estim_var_l;
539 agg_estim_var_l.size(num_lev);
540 Real agg_estim_var = 0;
541
542 agg_estim_var_l[0] = var_of_var_ml_l0(*static_sum_Ql, *static_sum_Qlm1, *static_sum_QlQlm1, Nlq_pilot, Nlq, qoi,
543 compute_gradient, grad_g[0][0]);
544 agg_estim_var += agg_estim_var_l[0];
545 for (lev = 1; lev < num_lev; ++lev) {
546 Nlq = x[lev];
547 Nlq_pilot = (*static_Nlq_pilot)[lev];
548
549 agg_estim_var_l[lev] = var_of_var_ml_l(*static_sum_Ql, *static_sum_Qlm1, *static_sum_QlQlm1, Nlq_pilot, Nlq, qoi, lev,
550 compute_gradient, grad_g[0][lev]);
551 agg_estim_var += agg_estim_var_l[lev];
552 }
553
554 g[0] = agg_estim_var; // - (*static_eps_sq_div_2);
555 }
556
557 /** This function performs "geometrical" MLMC on a single model form
558 with multiple discretization levels. */
multilevel_mc_Qsum(unsigned short model_form)559 void NonDMultilevelSampling::multilevel_mc_Qsum(unsigned short model_form)
560 {
561 // assign truth model form (solution level assignment is deferred until loop)
562 UShortArray truth_key;
563 unsigned short seq_index = 2, lev = USHRT_MAX; // lev updated in loop below
564 Pecos::DiscrepancyCalculator::form_key(0, model_form, lev, truth_key);
565 iteratedModel.active_model_key(truth_key);
566 Model& truth_model = iteratedModel.truth_model();
567
568 size_t qoi, num_steps = truth_model.solution_levels();//1 model form
569 unsigned short& step = (true) ? lev : model_form; // option not active
570
571 size_t max_iter = (maxIterations < 0) ? 25 : maxIterations; // default = -1
572 Real eps_sq_div_2, sum_sqrt_var_cost, estimator_var0 = 0., lev_cost, place_holder;
573 // retrieve cost estimates across soln levels for a particular model form
574 RealVector cost = truth_model.solution_level_costs(), agg_var(num_steps), agg_var_of_var(num_steps),
575 estimator_var0_qoi(numFunctions), eps_sq_div_2_qoi(numFunctions),
576 sum_sqrt_var_cost_qoi(numFunctions), sum_sqrt_var_cost_var_qoi(numFunctions), sum_sqrt_var_cost_mean_qoi(numFunctions),
577 agg_var_l_qoi(numFunctions), level_cost_vec(num_steps);
578 RealMatrix agg_var_qoi(numFunctions, num_steps), agg_var_mean_qoi(numFunctions, num_steps),
579 agg_var_var_qoi(numFunctions, num_steps), N_target_qoi(numFunctions, num_steps);
580 RealMatrix N_target_qoi_FN(numFunctions, num_steps);
581
582 // For moment estimation, we accumulate telescoping sums for Q^i using
583 // discrepancies Yi = Q^i_{lev} - Q^i_{lev-1} (Y_diff_Qpow[i] for i=1:4).
584 // For computing N_l from estimator variance, we accumulate square of Y1
585 // estimator (YY[i] = (Y^i)^2 for i=1).
586 IntRealMatrixMap sum_Ql, sum_Qlm1;
587 IntIntPairRealMatrixMap sum_QlQlm1;
588 initialize_ml_Qsums(sum_Ql, sum_Qlm1, sum_QlQlm1, num_steps);
589 IntIntPair pr11(1, 1);
590
591 bool have_npsol = false, have_optpp = false;
592 #ifdef HAVE_NPSOL
593 have_npsol = true;
594 #endif
595 #ifdef HAVE_OPTPP
596 have_optpp = true;
597 #endif
598
599 // Initialize for pilot sample
600 SizetArray delta_N_l;
601 IntMatrix delta_N_l_qoi(numFunctions, num_steps);
602 load_pilot_sample(pilotSamples, NLev, delta_N_l);
603
604 // raw eval counts are accumulation of allSamples irrespective of resp faults
605 SizetArray raw_N_l(num_steps, 0);
606 RealVectorArray mu_hat(num_steps);
607 Sizet2DArray &N_l = NLev[model_form];
608
609 // Safe guard.
610 for(qoi = 0; qoi < numFunctions; ++qoi) {
611 for (step = 0; step < num_steps; ++step) {
612 N_target_qoi(qoi, step) = -1;
613 }
614 }
615
616 // now converge on sample counts per level (N_l)
617 mlmfIter = 0;
618 while (Pecos::l1_norm(delta_N_l) && mlmfIter <= max_iter) {
619 Real underrelaxation_factor = static_cast<Real>(mlmfIter + 1)/static_cast<Real>(max_iter + 1);
620
621 sum_sqrt_var_cost = 0.;
622 for (qoi = 0; qoi < numFunctions; ++qoi) {
623 sum_sqrt_var_cost_qoi[qoi] = 0.;
624 sum_sqrt_var_cost_mean_qoi[qoi] = 0.;
625 sum_sqrt_var_cost_var_qoi[qoi] = 0.;
626 }
627 for (step=0; step<num_steps; ++step) {
628
629 configure_indices(step, model_form, lev, seq_index);
630 lev_cost = level_cost(cost, step);
631 level_cost_vec[step] = lev_cost;
632
633 // set the number of current samples from the defined increment
634 numSamples = delta_N_l[step];
635
636 // aggregate variances across QoI for estimating N_l (justification:
637 // for independent QoI, sum of QoI variances = variance of QoI sum)
638 Real &agg_var_l = agg_var[step]; // carried over from prev iter if no samp
639
640 if (numSamples) {
641
642 // advance any sequence specifications (seed_sequence)
643 assign_specification_sequence(step);
644 // generate new MC parameter sets
645 get_parameter_sets(iteratedModel);// pull dist params from any model
646
647 // export separate output files for each data set. truth_model()
648 // has the correct data when in bypass-surrogate mode.
649 if (exportSampleSets)
650 export_all_samples("ml_", iteratedModel.truth_model(),
651 mlmfIter, step);
652
653 // compute allResponses from allVariables using hierarchical model
654 evaluate_parameter_sets(iteratedModel, true, false);
655
656 // process allResponses: accumulate new samples for each qoi and
657 // update number of successful samples for each QoI
658 //if (mlmfIter == 0) accumulate_offsets(mu_hat[step]);
659
660 accumulate_ml_Qsums(sum_Ql, sum_Qlm1, sum_QlQlm1, step,
661 mu_hat[step], N_l[step]);
662 if (outputLevel == DEBUG_OUTPUT)
663 Cout << "Accumulated sums (Ql[1,2], Qlm1[1,2]):\n" << sum_Ql[1]
664 << sum_Ql[2] << sum_Qlm1[1] << sum_Qlm1[2] << std::endl;
665 // update raw evaluation counts
666 raw_N_l[step] += numSamples;
667
668 // compute estimator variance from current sample accumulation:
669 if (outputLevel >= DEBUG_OUTPUT)
670 Cout << "variance of Y[" << step << "]: ";
671 //if (target_mean)
672 agg_var_l = 0;
673 if (qoiAggregation==QOI_AGGREGATION_SUM) {
674 if (allocationTarget == TARGET_MEAN) {
675 agg_var_l = aggregate_variance_Qsum(sum_Ql[1][step], sum_Qlm1[1][step],
676 sum_Ql[2][step], sum_QlQlm1[pr11][step], sum_Qlm1[2][step],
677 N_l[step], step);
678 } else if (allocationTarget == TARGET_VARIANCE) {
679 for (qoi = 0; qoi < numFunctions; ++qoi) {
680 agg_var_l += ((step == 0) ? var_of_var_ml_l0(sum_Ql, sum_Qlm1, sum_QlQlm1, N_l[step][qoi],
681 N_l[step][qoi], qoi, false, place_holder)
682 : var_of_var_ml_l(sum_Ql, sum_Qlm1, sum_QlQlm1, N_l[step][qoi],
683 N_l[step][qoi], qoi, step, false, place_holder)) *
684 N_l[step][qoi]; //As described in the paper by Krumscheid, Pisaroni, Nobile
685 }
686 } else{
687 Cout << "NonDMultilevelSampling::multilevel_mc_Qsum: allocationTarget is not implemented.\n";
688 abort_handler(INTERFACE_ERROR);
689 }
690 check_negative(agg_var_l);
691 }else if (qoiAggregation==QOI_AGGREGATION_MAX) {
692 for (qoi = 0; qoi < numFunctions; ++qoi) {
693 agg_var_mean_qoi(qoi, step) = aggregate_variance_Qsum(sum_Ql[1][step], sum_Qlm1[1][step],
694 sum_Ql[2][step], sum_QlQlm1[pr11][step],
695 sum_Qlm1[2][step],
696 N_l[step], step, qoi);
697 agg_var_var_qoi(qoi, step) = ((step == 0) ? var_of_var_ml_l0(sum_Ql, sum_Qlm1, sum_QlQlm1, N_l[step][qoi],
698 N_l[step][qoi], qoi, false, place_holder)
699 : var_of_var_ml_l(sum_Ql, sum_Qlm1, sum_QlQlm1, N_l[step][qoi],
700 N_l[step][qoi], qoi, step, false, place_holder)) *
701 N_l[step][qoi]; //As described in the paper by Krumscheid, Pisaroni, Nobile
702 agg_var_qoi(qoi, step) = (allocationTarget == TARGET_MEAN) ? agg_var_mean_qoi(qoi, step) : agg_var_var_qoi(qoi, step);
703 check_negative(agg_var_qoi(qoi, step));
704 }
705 }else{
706 Cout << "NonDMultilevelSampling::multilevel_mc_Qsum: qoiAggregation is not known.\n";
707 abort_handler(INTERFACE_ERROR);
708 }
709 }
710
711 sum_sqrt_var_cost += std::sqrt(agg_var_l * lev_cost);
712 for (qoi = 0; qoi < numFunctions && qoiAggregation==QOI_AGGREGATION_MAX; ++qoi) {
713 sum_sqrt_var_cost_qoi[qoi] += std::sqrt(agg_var_qoi(qoi, step) * lev_cost);
714 }
715 // MSE reference is MC applied to HF:
716 if (mlmfIter == 0) {
717 estimator_var0 += aggregate_mse_Qsum(sum_Ql[1][step], sum_Qlm1[1][step],
718 sum_Ql[2][step], sum_QlQlm1[pr11][step], sum_Qlm1[2][step],
719 N_l[step], step);
720 for (qoi = 0; qoi < numFunctions; ++qoi) {
721 estimator_var0_qoi[qoi] += aggregate_mse_Qsum(sum_Ql[1][step], sum_Qlm1[1][step],
722 sum_Ql[2][step], sum_QlQlm1[pr11][step], sum_Qlm1[2][step],
723 N_l[step], step, qoi);
724 }
725 }
726 }
727 // compute epsilon target based on relative tolerance: total MSE = eps^2
728 // which is equally apportioned (eps^2 / 2) among discretization MSE and
729 // estimator variance (\Sum var_Y_l / N_l). Since we do not know the
730 // discretization error, we compute an initial estimator variance and
731 // then seek to reduce it by a relative_factor <= 1.
732 if (mlmfIter == 0) { // eps^2 / 2 = var * relative factor
733 //if(target_mean)
734 eps_sq_div_2 = estimator_var0 * convergenceTol;
735 for (qoi = 0; qoi < numFunctions; ++qoi) {
736 eps_sq_div_2_qoi[qoi] = estimator_var0_qoi[qoi] * convergenceTol; //1.389824213484928e-7; //2.23214285714257e-5; //estimator_var0_qoi[qoi] * convergenceTol;
737 }
738 if (outputLevel == DEBUG_OUTPUT) {
739 Cout << "Epsilon squared target = " << eps_sq_div_2 << std::endl;
740 Cout << "Epsilon squared target per QoI = " << eps_sq_div_2_qoi << std::endl;
741 }
742 }
743
744 // update targets based on variance estimates
745 //if(target_mean){
746 Real fact = sum_sqrt_var_cost / eps_sq_div_2, fact_qoi, N_target;
747 if (outputLevel == DEBUG_OUTPUT) Cout << "N_target: " << std::endl;
748 for (step = 0; step < num_steps && (qoiAggregation==QOI_AGGREGATION_SUM); ++step) {
749 // Equation 3.9 in CTR Annual Research Briefs:
750 // "A multifidelity control variate approach for the multilevel Monte
751 // Carlo technique," Geraci, Eldred, Iaccarino, 2015.
752 N_target = std::sqrt(agg_var[step] / level_cost_vec[step]) * fact;
753 for (qoi = 0; qoi < numFunctions; ++qoi) N_target_qoi(qoi, step) = N_target;
754 delta_N_l[step] = one_sided_delta(average(N_l[step]), N_target);
755
756 if (outputLevel == DEBUG_OUTPUT) Cout << N_target << " ";
757 }
758
759 for (qoi = 0; qoi < numFunctions && (qoiAggregation==QOI_AGGREGATION_MAX); ++qoi) {
760 fact_qoi = sum_sqrt_var_cost_qoi[qoi]/eps_sq_div_2_qoi[qoi];
761 if (outputLevel == DEBUG_OUTPUT)
762 Cout << "\n\tN_target for Qoi: " << qoi << ", with lagrange: " << fact_qoi << std::endl;
763 for (step = 0; step < num_steps; ++step) {
764 N_target_qoi(qoi, step) = std::sqrt(agg_var_qoi(qoi, step) / level_cost_vec[step]) * fact_qoi;
765 //N_target_qoi_FN(qoi, step) = N_target_qoi(qoi, step);
766 if (outputLevel == DEBUG_OUTPUT) {
767 Cout << "\t\tVar of target: " << agg_var_qoi(qoi, step) << std::endl;
768 Cout << "\t\tCost: " << level_cost_vec[step] << "\n";
769 Cout << "\t\tN_target_qoi: " << N_target_qoi(qoi, step) << "\n";
770 }
771 }
772 }
773 if (outputLevel == DEBUG_OUTPUT) Cout << "\n";
774
775 if( allocationTarget == TARGET_VARIANCE && (have_npsol || have_optpp) && useTargetVarianceOptimizationFlag){
776
777 Cout << "Numerical Optimization for sample allocation targeting variance using " << (have_npsol ? "NPSOL" : "OPTPP") << std::endl;
778 for (qoi = 0; qoi < numFunctions; ++qoi) {
779 RealVector initial_point, pilot_samples;
780 initial_point.size(N_l.size());
781 pilot_samples.size(N_l.size());
782 for (step = 0; step < N_l.size(); ++step) {
783 pilot_samples[step] = N_l[step][qoi];
784 initial_point[step] = 8. > N_target_qoi(qoi, step) ? 8 : N_target_qoi(qoi, step);
785 }
786
787 RealVector var_lower_bnds, var_upper_bnds, lin_ineq_lower_bnds, lin_ineq_upper_bnds, lin_eq_targets,
788 nonlin_ineq_lower_bnds, nonlin_ineq_upper_bnds, nonlin_eq_targets;
789 RealMatrix lin_ineq_coeffs, lin_eq_coeffs;
790
791 //Bound constraints only allowing positive values for Nlq
792 var_lower_bnds.size(num_steps); //init to 0
793 for (step = 0; step < N_l.size(); ++step) {
794 var_lower_bnds[step] = 6.;
795 }
796
797 var_upper_bnds.size(num_steps); //init to 0
798 var_upper_bnds.putScalar(1e10); //Set to high upper bound
799
800 //Number of linear inequality constraints = 0
801 lin_ineq_coeffs.shape(0, 0);
802 lin_ineq_lower_bnds.size(0);
803 lin_ineq_upper_bnds.size(0);
804
805 //Number of linear equality constraints = 0
806 lin_eq_coeffs.shape(0, 0);
807 lin_eq_targets.size(0);
808 //Number of nonlinear inequality bound constraints = 0
809 nonlin_ineq_lower_bnds.size(0);
810 nonlin_ineq_upper_bnds.size(0);
811 //Number of nonlinear equality constraints = 1, s.t. c_eq: c_1(Nlq) = convergenceTol;
812 nonlin_eq_targets.size(1); //init to 0
813 nonlin_eq_targets[0] = eps_sq_div_2_qoi[qoi]; //convergenceTol;
814
815 assign_static_member(nonlin_eq_targets[0], qoi, level_cost_vec, sum_Ql, sum_Qlm1, sum_QlQlm1, pilot_samples);
816
817 std::unique_ptr<Iterator> optimizer;
818 #ifdef HAVE_NPSOL
819 optimizer.reset(new NPSOLOptimizer(initial_point,
820 var_lower_bnds, var_upper_bnds,
821 lin_ineq_coeffs, lin_ineq_lower_bnds,
822 lin_ineq_upper_bnds, lin_eq_coeffs,
823 lin_eq_targets, nonlin_ineq_lower_bnds,
824 nonlin_ineq_upper_bnds, nonlin_eq_targets,
825 &target_var_objective_eval_npsol,
826 &target_var_constraint_eval_npsol,
827 3, 1e-15) //derivative_level = 3 means user_supplied gradients
828 );
829 #elif HAVE_OPTPP
830 optimizer.reset(new SNLLOptimizer(initial_point,
831 var_lower_bnds, var_upper_bnds,
832 lin_ineq_coeffs, lin_ineq_lower_bnds,
833 lin_ineq_upper_bnds, lin_eq_coeffs,
834 lin_eq_targets,nonlin_ineq_lower_bnds,
835 nonlin_ineq_upper_bnds, nonlin_eq_targets,
836 &target_var_objective_eval_optpp,
837 &target_var_constraint_eval_optpp,
838 100000, 100000, 1.e-14,
839 1.e-14, 100000)
840 );
841 #endif
842 optimizer->output_level(DEBUG_OUTPUT);
843 optimizer->run();
844
845
846 //Cout << optimizer->all_variables() << std::endl;
847 if (outputLevel == DEBUG_OUTPUT) {
848 Cout << "Optimization Run: Initial point: \n";
849 for (int i = 0; i < initial_point.length(); ++i) {
850 Cout << initial_point[i] << " ";
851 }
852 Cout << "\nOptimization Run. Best point: \n";
853 Cout << optimizer->variables_results().continuous_variables() << std::endl;
854 Cout << "Objective: " << optimizer->response_results().function_value(0) << std::endl;
855 Cout << "Constraint: " << optimizer->response_results().function_value(1) << std::endl;
856 Cout << "Relative Constraint violation: " << std::abs(1 - optimizer->response_results().function_value(1)/nonlin_eq_targets[0]) << std::endl;
857 Cout << "\n";
858 }
859
860 if(std::abs(1. - optimizer->response_results().function_value(1)/nonlin_eq_targets[0]) > 1.0e-5){
861 Cout << "Relative Constraint violation violated: Switching to log scale " << std::endl;
862 for (step = 0; step < N_l.size(); ++step) {
863 initial_point[step] = 8. > N_target_qoi(qoi, step) ? 8 : N_target_qoi(qoi, step); //optimizer->variables_results().continuous_variable(step) > pilot_samples[step] ? optimizer->variables_results().continuous_variable(step) : pilot_samples[step];
864 }
865 nonlin_eq_targets[0] = std::log(eps_sq_div_2_qoi[qoi]); //std::log(convergenceTol);
866 #ifdef HAVE_NPSOL
867 optimizer.reset(new NPSOLOptimizer(initial_point,
868 var_lower_bnds, var_upper_bnds,
869 lin_ineq_coeffs, lin_ineq_lower_bnds,
870 lin_ineq_upper_bnds, lin_eq_coeffs,
871 lin_eq_targets, nonlin_ineq_lower_bnds,
872 nonlin_ineq_upper_bnds, nonlin_eq_targets,
873 &target_var_objective_eval_npsol,
874 &target_var_constraint_eval_logscale_npsol,
875 3, 1e-15)
876 ); //derivative_level = 3 means user_supplied gradients
877 #elif HAVE_OPTPP
878 optimizer.reset(new SNLLOptimizer(initial_point,
879 var_lower_bnds, var_upper_bnds,
880 lin_ineq_coeffs, lin_ineq_lower_bnds,
881 lin_ineq_upper_bnds, lin_eq_coeffs,
882 lin_eq_targets, nonlin_ineq_lower_bnds,
883 nonlin_ineq_upper_bnds, nonlin_eq_targets,
884 &target_var_objective_eval_optpp,
885 &target_var_constraint_eval_logscale_optpp,
886 100000, 100000, 1.e-14,
887 1.e-14, 100000));
888 #endif
889 optimizer->run();
890
891 }
892
893 for (step=0; step<num_steps; ++step) {
894 N_target_qoi(qoi, step) = optimizer->variables_results().continuous_variable(step);
895 }
896
897 }
898
899 if (outputLevel == DEBUG_OUTPUT) {
900 Cout << "Final Optimization results: \n";
901 Cout << N_target_qoi << std::endl<< std::endl;
902 }
903 }
904
905 if(allocationTarget == TARGET_VARIANCE && qoiAggregation==QOI_AGGREGATION_SUM){
906 for (step = 0; step < num_steps; ++step) {
907
908 Real N_l_step_avg = 0;
909 Real N_l_target_step_avg = 0;
910 for (qoi = 0; qoi < numFunctions; ++qoi) {
911 N_l_step_avg += N_l[step][qoi];
912 N_l_target_step_avg += N_target_qoi(qoi, step);
913 }
914 N_l_step_avg /= numFunctions;
915 N_l_target_step_avg /= numFunctions;
916
917 delta_N_l[step] = std::ceil(std::min<Real>(N_l_step_avg*2, one_sided_delta( N_l_step_avg, N_l_target_step_avg)));
918 }
919 }
920 if(qoiAggregation==QOI_AGGREGATION_MAX) {
921 for (qoi = 0; qoi < numFunctions; ++qoi) {
922 for (step=0; step<num_steps; ++step) {
923 if(allocationTarget == TARGET_MEAN){
924 delta_N_l_qoi(qoi, step) = one_sided_delta(N_l[step][qoi], N_target_qoi(qoi, step));
925 }else if (allocationTarget == TARGET_VARIANCE){
926 delta_N_l_qoi(qoi, step) = std::min(N_l[step][qoi]*2, one_sided_delta(N_l[step][qoi], N_target_qoi(qoi, step)));
927 }else{
928 Cout << "NonDMultilevelSampling::multilevel_mc_Qsum: allocationTarget is not implemented.\n";
929 abort_handler(INTERFACE_ERROR);
930 }
931 }
932 }
933 Real max_qoi_idx = -1, max_cost = -1, cur_cost = 0;
934 for (step=0; step<num_steps; ++step) {
935 max_qoi_idx = 0;
936 for (qoi = 1; qoi < numFunctions; ++qoi) {
937 max_qoi_idx = delta_N_l_qoi(qoi, step) > delta_N_l_qoi(max_qoi_idx, step) ? qoi : max_qoi_idx;
938 }
939 delta_N_l[step] = delta_N_l_qoi(max_qoi_idx, step);
940 }
941 }
942
943 ++mlmfIter;
944 Cout << "\nMLMC iteration " << mlmfIter << " sample increments:\n"
945 << delta_N_l << std::endl;
946
947 }
948
949 // Roll up expected value estimators for central moments. Final expected
950 // value is sum of expected values from telescopic sum. Note: raw moments
951 // have no bias correction (no additional variance from an estimated center).
952 //RealMatrix Q_raw_mom(numFunctions, 4);
953 RealMatrix &sum_Q1l = sum_Ql[1], &sum_Q2l = sum_Ql[2],
954 &sum_Q3l = sum_Ql[3], &sum_Q4l = sum_Ql[4],
955 &sum_Q1lm1 = sum_Qlm1[1], &sum_Q2lm1 = sum_Qlm1[2],
956 &sum_Q3lm1 = sum_Qlm1[3], &sum_Q4lm1 = sum_Qlm1[4];
957 Real cm1, cm2, cm3, cm4, cm1l, cm2l, cm3l, cm4l;
958 if (momentStats.empty())
959 momentStats.shapeUninitialized(4, numFunctions);
960 for (qoi = 0; qoi < numFunctions; ++qoi) {
961 cm1 = cm2 = cm3 = cm4 = 0.;
962 for (step=0; step<num_steps; ++step) {
963 size_t Nlq = N_l[step][qoi];
964 // roll up unbiased moments centered on level mean
965 uncentered_to_centered(sum_Q1l(qoi, step) / Nlq, sum_Q2l(qoi, step) / Nlq,
966 sum_Q3l(qoi, step) / Nlq, sum_Q4l(qoi, step) / Nlq,
967 cm1l, cm2l, cm3l, cm4l, Nlq);
968 cm1 += cm1l;
969 cm2 += cm2l;
970 cm3 += cm3l;
971 cm4 += cm4l;
972 if (outputLevel == DEBUG_OUTPUT)
973 Cout << "CM_l for level " << step << ": "
974 << cm1l << ' ' << cm2l << ' ' << cm3l << ' ' << cm4l << '\n';
975 if (step) {
976 uncentered_to_centered(sum_Q1lm1(qoi, step) / Nlq, sum_Q2lm1(qoi, step) / Nlq,
977 sum_Q3lm1(qoi, step) / Nlq, sum_Q4lm1(qoi, step) / Nlq,
978 cm1l, cm2l, cm3l, cm4l, Nlq);
979 cm1 -= cm1l;
980 cm2 -= cm2l;
981 cm3 -= cm3l;
982 cm4 -= cm4l;
983 if (outputLevel == DEBUG_OUTPUT)
984 Cout << "CM_lm1 for level " << step << ": "
985 << cm1l << ' ' << cm2l << ' ' << cm3l << ' ' << cm4l << '\n';
986 }
987 }
988 check_negative(cm2);
989 check_negative(cm4);
990 Real *mom_q = momentStats[qoi];
991 if (finalMomentsType == CENTRAL_MOMENTS) {
992 mom_q[0] = cm1;
993 mom_q[1] = cm2;
994 mom_q[2] = cm3;
995 mom_q[3] = cm4;
996 } else
997 centered_to_standard(cm1, cm2, cm3, cm4,
998 mom_q[0], mom_q[1], mom_q[2], mom_q[3]);
999 }
1000
1001 // populate finalStatErrors
1002 compute_error_estimates(sum_Ql, sum_Qlm1, sum_QlQlm1, N_l);
1003
1004 // compute the equivalent number of HF evaluations (includes any sim faults)
1005 equivHFEvals = raw_N_l[0] * cost[0]; // first level is single eval
1006 for (step=1; step<num_steps; ++step)// subsequent levels incur 2 model costs
1007 equivHFEvals += raw_N_l[step] * (cost[step] + cost[step - 1]);
1008 equivHFEvals /= cost[num_steps - 1]; // normalize into equivalent HF evals
1009 }
1010
1011
assign_static_member(Real & conv_tol,size_t & qoi,RealVector & level_cost_vec,IntRealMatrixMap & sum_Ql,IntRealMatrixMap & sum_Qlm1,IntIntPairRealMatrixMap & sum_QlQlm1,RealVector & pilot_samples) const1012 void NonDMultilevelSampling::assign_static_member(Real &conv_tol, size_t &qoi, RealVector &level_cost_vec,
1013 IntRealMatrixMap &sum_Ql, IntRealMatrixMap &sum_Qlm1,
1014 IntIntPairRealMatrixMap &sum_QlQlm1,
1015 RealVector &pilot_samples) const
1016 {
1017 static_lev_cost_vec= &level_cost_vec;
1018 static_qoi = &qoi;
1019 static_sum_Ql = &sum_Ql;
1020 static_sum_Qlm1 = &sum_Qlm1;
1021 static_sum_QlQlm1 = &sum_QlQlm1;
1022 static_eps_sq_div_2 = &conv_tol;
1023 static_Nlq_pilot = &pilot_samples;
1024 }
1025
1026
1027 /** This function performs control variate MC across two combinations of
1028 model form and discretization level. */
1029 void NonDMultilevelSampling::
control_variate_mc(const UShortArray & active_key)1030 control_variate_mc(const UShortArray& active_key)
1031 {
1032 // Current implementation performs pilot + shared increment + LF increment,
1033 // where these increments are targeting a prescribed MSE reduction.
1034 // **********
1035 // *** TO DO: should CV MC iterate (new shared + LF increments)
1036 // *** until MSE target is met?
1037 // **********
1038
1039 aggregated_models_mode();
1040 iteratedModel.active_model_key(active_key); // data group 0
1041 Model& truth_model = iteratedModel.truth_model();
1042 Model& surr_model = iteratedModel.surrogate_model();
1043
1044 // retrieve active index
1045 //unsigned short lf_lev_index = surr_model.solution_level_index(),
1046 // hf_lev_index = truth_model.solution_level_index();
1047 // retrieve cost estimates across model forms for a particular soln level
1048 Real lf_cost = surr_model.solution_level_cost(),
1049 hf_cost = truth_model.solution_level_cost(),
1050 cost_ratio = hf_cost / lf_cost, avg_eval_ratio, avg_mse_ratio;
1051
1052 IntRealVectorMap sum_L_shared, sum_L_refined, sum_H, sum_LL, sum_LH;
1053 initialize_cv_sums(sum_L_shared, sum_L_refined, sum_H, sum_LL, sum_LH);
1054 RealVector sum_HH(numFunctions), var_H(numFunctions, false),
1055 rho2_LH(numFunctions, false);
1056
1057 // Initialize for pilot sample
1058 SizetArray delta_N_l;
1059 load_pilot_sample(pilotSamples, NLev, delta_N_l);
1060
1061 // NLev allocations currently enforce truncation to #HF levels (1)
1062 UShortArray hf_key, lf_key;
1063 Pecos::DiscrepancyCalculator::extract_keys(active_key, hf_key, lf_key);
1064 unsigned short hf_model_form = hf_key[1], lf_model_form = lf_key[1];
1065 SizetArray& N_lf = NLev[lf_model_form][0];//[lf_lev_index];
1066 SizetArray& N_hf = NLev[hf_model_form][0];//[hf_lev_index];
1067 size_t raw_N_lf = 0, raw_N_hf = 0;
1068 RealVector mu_hat;
1069
1070 // ---------------------
1071 // Compute Pilot Samples
1072 // ---------------------
1073
1074 mlmfIter = 0;
1075
1076 // Initialize for pilot sample (shared sample count discarding any excess)
1077 numSamples = std::min(delta_N_l[lf_model_form], delta_N_l[hf_model_form]);
1078 shared_increment(mlmfIter, 0);
1079 accumulate_cv_sums(sum_L_shared, sum_L_refined, sum_H, sum_LL, sum_LH,
1080 sum_HH, mu_hat, N_lf, N_hf);
1081 raw_N_lf += numSamples; raw_N_hf += numSamples;
1082
1083 // Compute the LF/HF evaluation ratio, averaged over the QoI.
1084 // This includes updating var_H and rho2_LH.
1085 avg_eval_ratio = eval_ratio(sum_L_shared[1], sum_H[1], sum_LL[1], sum_LH[1],
1086 sum_HH, cost_ratio, N_hf, var_H, rho2_LH);
1087 // compute the ratio of MC and CVMC mean squared errors (controls convergence)
1088 avg_mse_ratio = MSE_ratio(avg_eval_ratio, var_H, rho2_LH, mlmfIter, N_hf);
1089
1090 // ----------------------------------------------------------
1091 // Compute shared increment targeting specified MSE reduction
1092 // ----------------------------------------------------------
1093
1094 // bypass refinement if maxIterations == 0 or convergenceTol already
1095 // satisfied by pilot sample
1096 if (maxIterations && avg_mse_ratio > convergenceTol) {
1097
1098 // Assuming rho_AB, evaluation_ratio and var_H to be relatively invariant,
1099 // we seek a relative reduction in MSE using the convergence tol spec:
1100 // convTol = CV_mse / MC^0_mse = mse_ratio * N0 / N
1101 // delta_N = mse_ratio*N0/convTol - N0 = (mse_ratio/convTol - 1) * N0
1102 Real incr = (avg_mse_ratio / convergenceTol - 1.) * numSamples;
1103 numSamples = (size_t)std::floor(incr + .5); // round
1104
1105 if (numSamples) {
1106 shared_increment(++mlmfIter, 0);
1107 accumulate_cv_sums(sum_L_shared, sum_L_refined, sum_H, sum_LL, sum_LH,
1108 sum_HH, mu_hat, N_lf, N_hf);
1109 raw_N_lf += numSamples; raw_N_hf += numSamples;
1110 // update ratios:
1111 avg_eval_ratio = eval_ratio(sum_L_shared[1], sum_H[1], sum_LL[1],
1112 sum_LH[1], sum_HH, cost_ratio, N_hf,
1113 var_H, rho2_LH);
1114 avg_mse_ratio = MSE_ratio(avg_eval_ratio, var_H, rho2_LH,mlmfIter,N_hf);
1115 }
1116 }
1117
1118 // --------------------------------------------------
1119 // Compute LF increment based on the evaluation ratio
1120 // --------------------------------------------------
1121 uncorrected_surrogate_mode(); // also needed for assignment of lf_key below
1122 // Group id in lf_key is not currently important, since no SurrogateData
1123 // (correlations are computed based on the paired LF/HF data group, prior
1124 // to the augmentation, which could imply a future group segregation)
1125 iteratedModel.active_model_key(lf_key); // sets activeKey and surrModelKey
1126 if (lf_increment(avg_eval_ratio, N_lf, N_hf, ++mlmfIter, 0)) { // level 0
1127 accumulate_cv_sums(sum_L_refined, mu_hat, N_lf);
1128 raw_N_lf += numSamples;
1129 }
1130
1131 // Compute/apply control variate parameter to estimate uncentered raw moments
1132 RealMatrix H_raw_mom(numFunctions, 4);
1133 cv_raw_moments(sum_L_shared, sum_H, sum_LL, sum_LH, N_hf, sum_L_refined, N_lf,
1134 rho2_LH, H_raw_mom);
1135 // Convert uncentered raw moment estimates to final moments (central or std)
1136 convert_moments(H_raw_mom, momentStats);
1137
1138 // compute the equivalent number of HF evaluations
1139 equivHFEvals = raw_N_hf + (Real)raw_N_lf / cost_ratio;
1140 }
1141
1142
1143 /** This function performs "geometrical" MLMC across discretization
1144 levels for the high fidelity model form where CVMC si employed
1145 across two model forms to exploit correlation in the discrepancies
1146 at each level (Y_l). */
1147 void NonDMultilevelSampling::
multilevel_control_variate_mc_Ycorr(unsigned short lf_model_form,unsigned short hf_model_form)1148 multilevel_control_variate_mc_Ycorr(unsigned short lf_model_form,
1149 unsigned short hf_model_form)
1150 {
1151 // assign model forms (solution level assignments are deferred until loop)
1152 UShortArray active_key;
1153 unsigned short seq_index = 2, lev = USHRT_MAX; // lev updated in loop below
1154 Pecos::DiscrepancyCalculator::
1155 form_key(0, hf_model_form, lev, lf_model_form, lev, active_key);
1156 iteratedModel.active_model_key(active_key);
1157 Model& truth_model = iteratedModel.truth_model();
1158 Model& surr_model = iteratedModel.surrogate_model();
1159
1160 size_t qoi, num_hf_lev = truth_model.solution_levels(),
1161 num_cv_lev = std::min(num_hf_lev, surr_model.solution_levels());
1162 unsigned short& group = lev; // no alias switch for this algorithm
1163 size_t max_iter = (maxIterations < 0) ? 25 : maxIterations; // default = -1
1164 Real avg_eval_ratio, eps_sq_div_2, sum_sqrt_var_cost, estimator_var0 = 0.,
1165 lf_lev_cost, hf_lev_cost;
1166 // retrieve cost estimates across solution levels for HF model
1167 RealVector hf_cost = truth_model.solution_level_costs(),
1168 lf_cost = surr_model.solution_level_costs(), agg_var_hf(num_hf_lev),
1169 avg_eval_ratios(num_cv_lev);
1170 // For moment estimation, we accumulate telescoping sums for Q^i using
1171 // discrepancies Yi = Q^i_{lev} - Q^i_{lev-1} (Y_diff_Qpow[i] for i=1:4).
1172 // For computing N_l from estimator variance, we accumulate square of Y1
1173 // estimator (YY[1] = (Y^i)^2 for i=1).
1174 IntRealMatrixMap sum_L_refined, sum_L_shared, sum_H, sum_LL, sum_LH, sum_HH;
1175 initialize_mlcv_sums(sum_L_shared, sum_L_refined, sum_H, sum_LL, sum_LH,
1176 sum_HH, num_hf_lev, num_cv_lev);
1177 RealMatrix var_H(numFunctions, num_cv_lev, false),
1178 rho2_LH(numFunctions, num_cv_lev, false);
1179 RealVector Lambda(num_cv_lev, false), avg_rho2_LH(num_cv_lev, false);
1180
1181 // Initialize for pilot sample
1182 Sizet2DArray& N_lf = NLev[lf_model_form];
1183 Sizet2DArray& N_hf = NLev[hf_model_form];
1184 Sizet2DArray delta_N_l; load_pilot_sample(pilotSamples, NLev, delta_N_l);
1185 //SizetArray& delta_N_lf = delta_N_l[lf_model_form];
1186 SizetArray& delta_N_hf = delta_N_l[hf_model_form];
1187
1188 // raw eval counts are accumulation of allSamples irrespective of resp faults
1189 SizetArray raw_N_lf(num_cv_lev, 0), raw_N_hf(num_hf_lev, 0);
1190 RealVector mu_L_hat, mu_H_hat;
1191
1192 // now converge on sample counts per level (N_hf)
1193 mlmfIter = 0;
1194 while (Pecos::l1_norm(delta_N_hf) && mlmfIter <= max_iter) {
1195
1196 sum_sqrt_var_cost = 0.;
1197 for (lev=0; lev<num_hf_lev; ++lev) {
1198
1199 configure_indices(group, hf_model_form, lev, seq_index);
1200 hf_lev_cost = level_cost(hf_cost, lev);
1201
1202 // set the number of current samples from the defined increment
1203 numSamples = delta_N_hf[lev];
1204
1205 // aggregate variances across QoI for estimating N_hf (justification:
1206 // for independent QoI, sum of QoI variances = variance of QoI sum)
1207 Real& agg_var_hf_l = agg_var_hf[lev];//carried over from prev iter if!samp
1208 if (numSamples) {
1209
1210 // advance any sequence specifications (seed_sequence)
1211 assign_specification_sequence(lev);
1212 // generate new MC parameter sets
1213 get_parameter_sets(iteratedModel);// pull dist params from any model
1214
1215 // export separate output files for each data set. Note that
1216 // truth_model() is indexed with hf_model_form at this stage for
1217 // all levels. The exported discretization level (e.g., state variable
1218 // value) can't capture a level discrepancy for lev>0 and will reflect
1219 // the most recent evaluation state.
1220 if (exportSampleSets)
1221 export_all_samples("ml_", iteratedModel.truth_model(), mlmfIter, lev);
1222
1223 // compute allResponses from allVariables using hierarchical model
1224 evaluate_parameter_sets(iteratedModel, true, false);
1225
1226 // control variate betwen LF and HF for this discretization level:
1227 // if unequal number of levels, loop over all HF levels for MLMC and
1228 // apply CVMC when LF levels are available. LF levels are assigned as
1229 // control variates to the leading set of HF levels, since these will
1230 // tend to have larger variance.
1231 if (lev < num_cv_lev) {
1232
1233 // store allResponses used for sum_H (and sum_HH)
1234 IntResponseMap hf_resp = allResponses; // shallow copy
1235 // activate LF response (lev 0) or LF response discrepancy (lev > 0)
1236 // within the hierarchical surrogate model. Level indices & surrogate
1237 // response mode are same as HF above, only the model form changes.
1238 // However, we must pass the unchanged level index to update the
1239 // corresponding variable values for the new model form.
1240 configure_indices(group, lf_model_form, lev, seq_index);
1241 lf_lev_cost = level_cost(lf_cost, lev);
1242 // compute allResp w/ LF model form reusing allVars from MLMC step
1243 evaluate_parameter_sets(iteratedModel, true, false);
1244 // process previous and new set of allResponses for CV sums
1245 accumulate_mlcv_Ysums(allResponses, hf_resp, sum_L_shared,
1246 sum_L_refined, sum_H, sum_LL, sum_LH,
1247 sum_HH, lev, mu_L_hat, mu_H_hat,
1248 N_lf[lev], N_hf[lev]);
1249 if (outputLevel == DEBUG_OUTPUT)
1250 Cout << "Accumulated sums (L_shared[1,2], L_refined[1,2], LH[1,2])"
1251 << ":\n" << sum_L_shared[1] << sum_L_shared[2]
1252 << sum_L_refined[1] << sum_L_refined[2]
1253 << sum_LH[1] << sum_LH[2];
1254 // update raw evaluation counts
1255 raw_N_lf[lev] += numSamples; raw_N_hf[lev] += numSamples;
1256
1257 // compute the average evaluation ratio and Lambda factor
1258 avg_eval_ratio = avg_eval_ratios[lev] =
1259 eval_ratio(sum_L_shared[1], sum_H[1], sum_LL[1], sum_LH[1],
1260 sum_HH[1], hf_lev_cost/lf_lev_cost, lev, N_hf[lev],
1261 var_H, rho2_LH);
1262 avg_rho2_LH[lev] = average(rho2_LH[lev], numFunctions);
1263 Lambda[lev] = 1. - avg_rho2_LH[lev]
1264 * (avg_eval_ratio - 1.) / avg_eval_ratio;
1265 agg_var_hf_l = sum(var_H[lev], numFunctions);
1266 }
1267 else { // no LF model for this level; accumulate only multilevel sums
1268 RealMatrix& sum_HH1 = sum_HH[1];
1269 accumulate_ml_Ysums(sum_H, sum_HH1, lev, mu_H_hat,
1270 N_hf[lev]); // sum_Y for lev>0
1271 if (outputLevel == DEBUG_OUTPUT)
1272 Cout << "Accumulated sums (H[1], H[2], HH):\n"
1273 << sum_H[1] << sum_H[2] << sum_HH1;
1274 raw_N_hf[lev] += numSamples;
1275 // aggregate Y variances across QoI for this level
1276 if (outputLevel >= DEBUG_OUTPUT)
1277 Cout << "variance of Y[" << lev << "]: ";
1278 agg_var_hf_l
1279 = aggregate_variance_Ysum(sum_H[1][lev], sum_HH1[lev], N_hf[lev]);
1280 }
1281 }
1282
1283 // accumulate sum of sqrt's of estimator var * cost used in N_target
1284 sum_sqrt_var_cost += (lev < num_cv_lev) ?
1285 std::sqrt(agg_var_hf_l * hf_lev_cost / (1. - avg_rho2_LH[lev]))
1286 * Lambda[lev] : std::sqrt(agg_var_hf_l * hf_lev_cost);
1287 // MSE reference is MLMF MC applied to {HF,LF} pilot sample aggregated
1288 // across qoi. Note: if the pilot sample for LF is not shaped, then r=1
1289 // will result in no additional variance reduction beyond MLMC.
1290 if (mlmfIter == 0)
1291 estimator_var0 += (lev < num_cv_lev) ?
1292 aggregate_mse_Yvar(var_H[lev], N_hf[lev]) :
1293 aggregate_mse_Ysum(sum_H[1][lev], sum_HH[1][lev], N_hf[lev]);
1294 }
1295 // compute epsilon target based on relative tolerance: total MSE = eps^2
1296 // which is equally apportioned (eps^2 / 2) among discretization MSE and
1297 // estimator variance (\Sum var_Y_l / N_l). Since we do not know the
1298 // discretization error, we compute an initial estimator variance and
1299 // then seek to reduce it by a relative_factor <= 1.
1300 if (mlmfIter == 0) { // eps^2 / 2 = var * relative factor
1301 eps_sq_div_2 = estimator_var0 * convergenceTol;
1302 if (outputLevel == DEBUG_OUTPUT)
1303 Cout << "Epsilon squared target = " << eps_sq_div_2 << std::endl;
1304 }
1305
1306 // All CV lf_increment() calls now follow all ML level evals:
1307 for (lev=0; lev<num_cv_lev; ++lev) {
1308 if (delta_N_hf[lev]) {
1309 configure_indices(group, lf_model_form, lev, seq_index);//augment LF grp
1310
1311 // execute additional LF sample increment, if needed
1312 if (lf_increment(avg_eval_ratios[lev], N_lf[lev], N_hf[lev],
1313 mlmfIter, lev)) {
1314 accumulate_mlcv_Ysums(sum_L_refined, lev, mu_L_hat, N_lf[lev]);
1315 raw_N_lf[lev] += numSamples;
1316 if (outputLevel == DEBUG_OUTPUT)
1317 Cout << "Accumulated sums (L_refined[1,2]):\n"
1318 << sum_L_refined[1] << sum_L_refined[2];
1319 }
1320 }
1321 }
1322
1323 // update targets based on variance estimates
1324 Real fact = sum_sqrt_var_cost / eps_sq_div_2, N_target;
1325 for (lev=0; lev<num_hf_lev; ++lev) {
1326 hf_lev_cost = (lev) ? hf_cost[lev] + hf_cost[lev-1] : hf_cost[lev];
1327 N_target = (lev < num_cv_lev) ? fact *
1328 std::sqrt(agg_var_hf[lev] / hf_lev_cost * (1. - avg_rho2_LH[lev])) :
1329 fact * std::sqrt(agg_var_hf[lev] / hf_lev_cost);
1330 delta_N_hf[lev] = one_sided_delta(average(N_hf[lev]), N_target);
1331 }
1332 ++mlmfIter;
1333 Cout << "\nMLCVMC iteration " << mlmfIter << " sample increments:\n"
1334 << delta_N_hf << std::endl;
1335 }
1336
1337 // Iteration complete. Now roll up raw moments from CVMC and MLMC estimators.
1338 RealMatrix Y_mlmc_mom(numFunctions, 4), Y_cvmc_mom(numFunctions, 4, false);
1339 for (lev=0; lev<num_cv_lev; ++lev) {
1340 cv_raw_moments(sum_L_shared, sum_H, sum_LL, sum_LH, N_hf[lev],
1341 sum_L_refined, N_lf[lev], rho2_LH, lev, Y_cvmc_mom);
1342 Y_mlmc_mom += Y_cvmc_mom;
1343 }
1344 if (num_hf_lev > num_cv_lev) {
1345 RealMatrix &sum_H1 = sum_H[1], &sum_H2 = sum_H[2],
1346 &sum_H3 = sum_H[3], &sum_H4 = sum_H[4];
1347 for (qoi=0; qoi<numFunctions; ++qoi) {
1348 for (lev=num_cv_lev; lev<num_hf_lev; ++lev) {
1349 size_t Nlq = N_hf[lev][qoi];
1350 Y_mlmc_mom(qoi,0) += sum_H1(qoi,lev) / Nlq;
1351 Y_mlmc_mom(qoi,1) += sum_H2(qoi,lev) / Nlq;
1352 Y_mlmc_mom(qoi,2) += sum_H3(qoi,lev) / Nlq;
1353 Y_mlmc_mom(qoi,3) += sum_H4(qoi,lev) / Nlq;
1354 }
1355 }
1356 }
1357 // Convert uncentered raw moment estimates to final moments (central or std)
1358 convert_moments(Y_mlmc_mom, momentStats);
1359
1360 // compute the equivalent number of HF evaluations
1361 equivHFEvals = raw_N_hf[0] * hf_cost[0] + raw_N_lf[0] * lf_cost[0]; // 1st lev
1362 for (lev=1; lev<num_hf_lev; ++lev) // subsequent levels incur 2 model costs
1363 equivHFEvals += raw_N_hf[lev] * (hf_cost[lev] + hf_cost[lev-1]);
1364 for (lev=1; lev<num_cv_lev; ++lev) // subsequent levels incur 2 model costs
1365 equivHFEvals += raw_N_lf[lev] * (lf_cost[lev] + lf_cost[lev-1]);
1366 equivHFEvals /= hf_cost[num_hf_lev-1]; // normalize into equivalent HF evals
1367 }
1368
1369
1370 /** This function performs "geometrical" MLMC across discretization
1371 levels for the high fidelity model form where CVMC is employed
1372 across two model forms. It generalizes the Y_l correlation case
1373 to separately target correlations for each QoI level embedded
1374 within the level discrepancies. */
1375 void NonDMultilevelSampling::
multilevel_control_variate_mc_Qcorr(unsigned short lf_model_form,unsigned short hf_model_form)1376 multilevel_control_variate_mc_Qcorr(unsigned short lf_model_form,
1377 unsigned short hf_model_form)
1378 {
1379 // assign model forms (solution level assignments are deferred until loop)
1380 UShortArray active_key;
1381 unsigned short seq_index = 2, lev = USHRT_MAX; // lev updated in loop below
1382 Pecos::DiscrepancyCalculator::
1383 form_key(0, hf_model_form, lev, lf_model_form, lev, active_key);
1384 iteratedModel.active_model_key(active_key);
1385 Model& truth_model = iteratedModel.truth_model();
1386 Model& surr_model = iteratedModel.surrogate_model();
1387
1388 size_t qoi, num_hf_lev = truth_model.solution_levels(),
1389 num_cv_lev = std::min(num_hf_lev, surr_model.solution_levels());
1390 unsigned short& group = lev; // no alias switch for this algorithm
1391 size_t max_iter = (maxIterations < 0) ? 25 : maxIterations; // default = -1
1392 Real avg_eval_ratio, eps_sq_div_2, sum_sqrt_var_cost, estimator_var0 = 0.,
1393 lf_lev_cost, hf_lev_cost;
1394 // retrieve cost estimates across solution levels for HF model
1395 RealVector hf_cost = truth_model.solution_level_costs(),
1396 lf_cost = surr_model.solution_level_costs(), agg_var_hf(num_hf_lev),
1397 avg_eval_ratios(num_cv_lev);
1398
1399 // CV requires cross-level covariance combinations in Qcorr approach
1400 IntRealMatrixMap sum_Ll, sum_Llm1,
1401 sum_Ll_refined, sum_Llm1_refined, sum_Hl, sum_Hlm1,
1402 sum_Ll_Ll, sum_Ll_Llm1, // for Var(Q_l^L), Covar(Q_l^L,Q_lm1^L)
1403 sum_Llm1_Llm1, sum_Hl_Ll, // for Var(Q_lm1^L), Covar(Q_l^H,Q_l^L)
1404 sum_Hl_Llm1, sum_Hlm1_Ll, // for Covar(Q_l^H,Q_lm1^L), Covar(Q_lm1^H,Q_l^L)
1405 sum_Hlm1_Llm1, // for Covar(Q_lm1^H,Q_lm1^L)
1406 sum_Hl_Hl, // for Var(Q_l^H)
1407 sum_Hl_Hlm1, // for Covar(Q_l^H,Q_lm1^H)
1408 sum_Hlm1_Hlm1; // for Var(Q_lm1^H)
1409 // Initialize accumulators and related arrays/maps, allowing for different
1410 // number of ML and CV levels (num_hf_lev & num_cv_lev, respectively).
1411 initialize_mlcv_sums(sum_Ll, sum_Llm1, sum_Ll_refined, sum_Llm1_refined,
1412 sum_Hl, sum_Hlm1, sum_Ll_Ll, sum_Ll_Llm1, sum_Llm1_Llm1,
1413 sum_Hl_Ll, sum_Hl_Llm1, sum_Hlm1_Ll, sum_Hlm1_Llm1,
1414 sum_Hl_Hl, sum_Hl_Hlm1, sum_Hlm1_Hlm1, num_hf_lev,
1415 num_cv_lev);
1416 RealMatrix var_Yl(numFunctions, num_cv_lev, false),
1417 rho_dot2_LH(numFunctions, num_cv_lev, false);
1418 RealVector Lambda(num_cv_lev, false), avg_rho_dot2_LH(num_cv_lev, false);
1419
1420 // Initialize for pilot sample
1421 Sizet2DArray& N_lf = NLev[lf_model_form];
1422 Sizet2DArray& N_hf = NLev[hf_model_form];
1423 Sizet2DArray delta_N_l; load_pilot_sample(pilotSamples, NLev, delta_N_l);
1424 //SizetArray& delta_N_lf = delta_N_l[lf_model_form];
1425 SizetArray& delta_N_hf = delta_N_l[hf_model_form];
1426
1427 // raw eval counts are accumulation of allSamples irrespective of resp faults
1428 SizetArray raw_N_lf(num_cv_lev, 0), raw_N_hf(num_hf_lev, 0);
1429 RealVector mu_L_hat, mu_H_hat;
1430
1431 // now converge on sample counts per level (N_hf)
1432 mlmfIter = 0;
1433 while (Pecos::l1_norm(delta_N_hf) && mlmfIter <= max_iter) {
1434
1435 sum_sqrt_var_cost = 0.;
1436 for (lev=0; lev<num_hf_lev; ++lev) {
1437
1438 configure_indices(group, hf_model_form, lev, seq_index);
1439 hf_lev_cost = level_cost(hf_cost, lev);
1440
1441 // set the number of current samples from the defined increment
1442 numSamples = delta_N_hf[lev];
1443
1444 // aggregate variances across QoI for estimating N_hf (justification:
1445 // for independent QoI, sum of QoI variances = variance of QoI sum)
1446 Real& agg_var_hf_l = agg_var_hf[lev];//carried over from prev iter if!samp
1447 if (numSamples) {
1448
1449 // advance any sequence specifications (seed_sequence)
1450 assign_specification_sequence(lev);
1451 // generate new MC parameter sets
1452 get_parameter_sets(iteratedModel);// pull dist params from any model
1453
1454 // export separate output files for each data set. Note that
1455 // truth_model() is indexed with hf_model_form at this stage for
1456 // all levels. The exported discretization level (e.g., state variable
1457 // value) can't capture a level discrepancy for lev>0 and will reflect
1458 // the most recent evaluation state.
1459 if (exportSampleSets)
1460 export_all_samples("ml_", iteratedModel.truth_model(), mlmfIter, lev);
1461
1462 // compute allResponses from allVariables using hierarchical model
1463 evaluate_parameter_sets(iteratedModel, true, false);
1464
1465 // control variate betwen LF and HF for this discretization level:
1466 // if unequal number of levels, loop over all HF levels for MLMC and
1467 // apply CVMC when LF levels are available. LF levels are assigned as
1468 // control variates to the leading set of HF levels, since these will
1469 // tend to have larger variance.
1470 if (lev < num_cv_lev) {
1471
1472 // store allResponses used for sum_H (and sum_HH)
1473 IntResponseMap hf_resp = allResponses; // shallow copy is sufficient
1474 // activate LF response (lev 0) or LF response discrepancy (lev > 0)
1475 // within the hierarchical surrogate model. Level indices & surrogate
1476 // response mode are same as HF above, only the model form changes.
1477 // However, we must pass the unchanged level index to update the
1478 // corresponding variable values for the new model form.
1479 configure_indices(group, lf_model_form, lev, seq_index);
1480 lf_lev_cost = level_cost(lf_cost, lev);
1481 // eval allResp w/ LF model reusing allVars from ML step above
1482 evaluate_parameter_sets(iteratedModel, true, false);
1483 // process previous and new set of allResponses for MLCV sums;
1484 accumulate_mlcv_Qsums(allResponses, hf_resp, sum_Ll, sum_Llm1,
1485 sum_Ll_refined, sum_Llm1_refined, sum_Hl,
1486 sum_Hlm1, sum_Ll_Ll, sum_Ll_Llm1, sum_Llm1_Llm1,
1487 sum_Hl_Ll, sum_Hl_Llm1, sum_Hlm1_Ll,
1488 sum_Hlm1_Llm1, sum_Hl_Hl, sum_Hl_Hlm1,
1489 sum_Hlm1_Hlm1, lev, mu_L_hat, mu_H_hat,
1490 N_lf[lev], N_hf[lev]);
1491 if (outputLevel == DEBUG_OUTPUT)
1492 Cout << "Accumulated sums (Ll[1,2], L_refined[1,2], Hl[1,2]):\n"
1493 << sum_Ll[1] << sum_Ll[2] << sum_Ll_refined[1]
1494 << sum_Ll_refined[2] << sum_Hl[1] << sum_Hl[2];
1495 // update raw evaluation counts
1496 raw_N_lf[lev] += numSamples; raw_N_hf[lev] += numSamples;
1497
1498 // compute the average evaluation ratio and Lambda factor
1499 avg_eval_ratio = avg_eval_ratios[lev] =
1500 eval_ratio(sum_Ll[1], sum_Llm1[1], sum_Hl[1], sum_Hlm1[1],
1501 sum_Ll_Ll[1], sum_Ll_Llm1[1], sum_Llm1_Llm1[1],
1502 sum_Hl_Ll[1], sum_Hl_Llm1[1], sum_Hlm1_Ll[1],
1503 sum_Hlm1_Llm1[1], sum_Hl_Hl[1], sum_Hl_Hlm1[1],
1504 sum_Hlm1_Hlm1[1], hf_lev_cost/lf_lev_cost, lev,
1505 N_hf[lev], var_Yl, rho_dot2_LH);
1506 avg_rho_dot2_LH[lev] = average(rho_dot2_LH[lev], numFunctions);
1507 Lambda[lev] = 1. - avg_rho_dot2_LH[lev]
1508 * (avg_eval_ratio - 1.) / avg_eval_ratio;
1509 agg_var_hf_l = sum(var_Yl[lev], numFunctions);
1510 }
1511 else { // no LF model for this level; accumulate only multilevel
1512 // discrepancy sums (Hl is Yl) as in standard MLMC
1513 RealMatrix& sum_HH1 = sum_Hl_Hl[1];
1514 accumulate_ml_Ysums(sum_Hl, sum_HH1, lev, mu_H_hat,
1515 N_hf[lev]); // sum_Y for lev>0
1516 if (outputLevel == DEBUG_OUTPUT)
1517 Cout << "Accumulated sums (H[1], H[2], HH[1]):\n"
1518 << sum_Hl[1] << sum_Hl[2] << sum_HH1;
1519 raw_N_hf[lev] += numSamples;
1520 // aggregate Y variances across QoI for this level
1521 if (outputLevel >= DEBUG_OUTPUT)
1522 Cout << "variance of Y[" << lev << "]: ";
1523 agg_var_hf_l
1524 = aggregate_variance_Ysum(sum_Hl[1][lev], sum_HH1[lev], N_hf[lev]);
1525 }
1526 }
1527
1528 // accumulate sum of sqrt's of estimator var * cost used in N_target
1529 sum_sqrt_var_cost += (lev < num_cv_lev) ?
1530 std::sqrt(agg_var_hf_l * hf_lev_cost / (1. - avg_rho_dot2_LH[lev]))
1531 * Lambda[lev] : std::sqrt(agg_var_hf_l * hf_lev_cost);
1532 // MSE reference is MLMF MC applied to {HF,LF} pilot sample aggregated
1533 // across qoi. Note: if the pilot sample for LF is not shaped, then r=1
1534 // will result in no additional variance reduction beyond MLMC.
1535 if (mlmfIter == 0)
1536 estimator_var0 += (lev < num_cv_lev) ?
1537 aggregate_mse_Yvar(var_Yl[lev], N_hf[lev]) :
1538 aggregate_mse_Ysum(sum_Hl[1][lev], sum_Hl_Hl[1][lev], N_hf[lev]);
1539 }
1540 // compute epsilon target based on relative tolerance: total MSE = eps^2
1541 // which is equally apportioned (eps^2 / 2) among discretization MSE and
1542 // estimator variance (\Sum var_Y_l / N_l). Since we do not know the
1543 // discretization error, we compute an initial estimator variance and
1544 // then seek to reduce it by a relative_factor <= 1.
1545 if (mlmfIter == 0) { // eps^2 / 2 = var * relative factor
1546 eps_sq_div_2 = estimator_var0 * convergenceTol;
1547 if (outputLevel == DEBUG_OUTPUT)
1548 Cout << "Epsilon squared target = " << eps_sq_div_2 << std::endl;
1549 }
1550
1551 // All CV lf_increment() calls now follow all ML level evals:
1552 // > Provides separation of pilot sample from refinements (simplifying
1553 // offline execution with data importing w/o undesirable seed progression)
1554 // > Improves application of max_iterations control in general: user
1555 // specification results in consistent count for ML and CV refinements
1556 // > Incurs a bit more overhead: avg_eval_ratios array, mode resetting
1557 // > Could potentially have parallel scheduling benefits by grouping
1558 // similar Model eval sets for aggregated scheduling
1559 for (lev=0; lev<num_cv_lev; ++lev) {
1560 if (delta_N_hf[lev]) {
1561 configure_indices(group, lf_model_form, lev, seq_index);//augment LF grp
1562
1563 // now execute additional LF sample increment, if needed
1564 if (lf_increment(avg_eval_ratios[lev], N_lf[lev], N_hf[lev],
1565 mlmfIter, lev)) {
1566 accumulate_mlcv_Qsums(sum_Ll_refined, sum_Llm1_refined, lev, mu_L_hat,
1567 N_lf[lev]);
1568 raw_N_lf[lev] += numSamples;
1569 if (outputLevel == DEBUG_OUTPUT)
1570 Cout << "Accumulated sums (L_refined[1,2]):\n"
1571 << sum_Ll_refined[1] << sum_Ll_refined[2];
1572 }
1573 }
1574 }
1575
1576 // update targets based on variance estimates
1577 Real fact = sum_sqrt_var_cost / eps_sq_div_2, N_target;
1578 for (lev=0; lev<num_hf_lev; ++lev) {
1579 hf_lev_cost = (lev) ? hf_cost[lev] + hf_cost[lev-1] : hf_cost[lev];
1580 N_target = (lev < num_cv_lev) ? fact *
1581 std::sqrt(agg_var_hf[lev] / hf_lev_cost * (1. - avg_rho_dot2_LH[lev])) :
1582 fact * std::sqrt(agg_var_hf[lev] / hf_lev_cost);
1583 delta_N_hf[lev] = one_sided_delta(average(N_hf[lev]), N_target);
1584 }
1585 ++mlmfIter;
1586 Cout << "\nMLCVMC iteration " << mlmfIter << " sample increments:\n"
1587 << delta_N_hf << std::endl;
1588 }
1589
1590 // Iteration complete. Now roll up raw moments from CVMC and MLMC estimators.
1591 RealMatrix Y_mlmc_mom(numFunctions, 4), Y_cvmc_mom(numFunctions, 4, false);
1592 for (lev=0; lev<num_cv_lev; ++lev) {
1593 cv_raw_moments(sum_Ll, sum_Llm1, sum_Hl, sum_Hlm1, sum_Ll_Ll, sum_Ll_Llm1,
1594 sum_Llm1_Llm1, sum_Hl_Ll, sum_Hl_Llm1, sum_Hlm1_Ll,
1595 sum_Hlm1_Llm1, sum_Hl_Hl, sum_Hl_Hlm1, sum_Hlm1_Hlm1,
1596 N_hf[lev], sum_Ll_refined, sum_Llm1_refined, N_lf[lev],
1597 rho_dot2_LH, lev, Y_cvmc_mom);
1598 Y_mlmc_mom += Y_cvmc_mom;
1599 }
1600 if (num_hf_lev > num_cv_lev) {
1601 // MLMC without CV: sum_H = HF Q sums for lev = 0 and HF Y sums for lev > 0
1602 RealMatrix &sum_H1 = sum_Hl[1], &sum_H2 = sum_Hl[2],
1603 &sum_H3 = sum_Hl[3], &sum_H4 = sum_Hl[4];
1604 for (qoi=0; qoi<numFunctions; ++qoi) {
1605 for (lev=num_cv_lev; lev<num_hf_lev; ++lev) {
1606 size_t Nlq = N_hf[lev][qoi];
1607 Y_mlmc_mom(qoi,0) += sum_H1(qoi,lev) / Nlq;
1608 Y_mlmc_mom(qoi,1) += sum_H2(qoi,lev) / Nlq;
1609 Y_mlmc_mom(qoi,2) += sum_H3(qoi,lev) / Nlq;
1610 Y_mlmc_mom(qoi,3) += sum_H4(qoi,lev) / Nlq;
1611 }
1612 }
1613 }
1614 // Convert uncentered raw moment estimates to final moments (central or std)
1615 convert_moments(Y_mlmc_mom, momentStats);
1616
1617 // compute the equivalent number of HF evaluations
1618 equivHFEvals = raw_N_hf[0] * hf_cost[0] + raw_N_lf[0] * lf_cost[0]; // 1st lev
1619 for (lev=1; lev<num_hf_lev; ++lev) // subsequent levels incur 2 model costs
1620 equivHFEvals += raw_N_hf[lev] * (hf_cost[lev] + hf_cost[lev-1]);
1621 for (lev=1; lev<num_cv_lev; ++lev) // subsequent levels incur 2 model costs
1622 equivHFEvals += raw_N_lf[lev] * (lf_cost[lev] + lf_cost[lev-1]);
1623 equivHFEvals /= hf_cost[num_hf_lev-1]; // normalize into equivalent HF evals
1624 }
1625
1626
1627 void NonDMultilevelSampling::
configure_indices(unsigned short group,unsigned short form,unsigned short lev,unsigned short s_index)1628 configure_indices(unsigned short group, unsigned short form,
1629 unsigned short lev, unsigned short s_index)
1630 {
1631 // Notes:
1632 // > could consolidate with NonDExpansion::configure_indices() with a passed
1633 // model and virtual *_mode() assignments. Leaving separate for now...
1634 // > group index is assigned based on step in model form/resolution sequence
1635 // > CVMC does not use this helper; it requires uncorrected_surrogate_mode()
1636
1637 UShortArray hf_key;
1638 Pecos::DiscrepancyCalculator::form_key(group, form, lev, hf_key);
1639
1640 if (hf_key[s_index] == 0) { // step 0 in the sequence
1641 bypass_surrogate_mode();
1642 iteratedModel.active_model_key(hf_key); // one active fidelity
1643 }
1644 else { //if (multilevDiscrepEmulation == DISTINCT_EMULATION) {
1645 aggregated_models_mode();
1646
1647 UShortArray lf_key(hf_key), aggregate_key;
1648 Pecos::DiscrepancyCalculator::decrement_key(lf_key, s_index);
1649 Pecos::DiscrepancyCalculator::aggregate_keys(hf_key, lf_key, aggregate_key);
1650 iteratedModel.active_model_key(aggregate_key); // two active fidelities
1651 }
1652 }
1653
1654
assign_specification_sequence(size_t index)1655 void NonDMultilevelSampling::assign_specification_sequence(size_t index)
1656 {
1657 // Note: seedSpec/randomSeed initialized from randomSeedSeqSpec in ctor
1658
1659 // advance any sequence specifications, as admissible
1660 // Note: no colloc pts sequence as load_pilot_sample() handles this separately
1661 int seed_i = random_seed(index);
1662 if (seed_i) randomSeed = seed_i;// propagate to NonDSampling::initialize_lhs()
1663 // else previous value will allow existing RNG to continue for varyPattern
1664 }
1665
1666
1667 void NonDMultilevelSampling::
initialize_ml_Ysums(IntRealMatrixMap & sum_Y,size_t num_lev)1668 initialize_ml_Ysums(IntRealMatrixMap& sum_Y, size_t num_lev)
1669 {
1670 // sum_* are running sums across all increments
1671 std::pair<int, RealMatrix> empty_pr;
1672 for (int i=1; i<=4; ++i) {
1673 empty_pr.first = i;
1674 // std::map::insert() returns std::pair<IntRMMIter, bool>:
1675 // use iterator to shape RealMatrix in place and init sums to 0
1676 sum_Y.insert(empty_pr).first->second.shape(numFunctions, num_lev);
1677 }
1678 }
1679
1680
1681 void NonDMultilevelSampling::
initialize_ml_Qsums(IntRealMatrixMap & sum_Ql,IntRealMatrixMap & sum_Qlm1,IntIntPairRealMatrixMap & sum_QlQlm1,size_t num_lev)1682 initialize_ml_Qsums(IntRealMatrixMap& sum_Ql, IntRealMatrixMap& sum_Qlm1,
1683 IntIntPairRealMatrixMap& sum_QlQlm1, size_t num_lev)
1684 {
1685 // sum_* are running sums across all increments
1686 std::pair<int, RealMatrix> empty_irm_pr; int i, j;
1687 for (i=1; i<=4; ++i) {
1688 empty_irm_pr.first = i;
1689 // std::map::insert() returns std::pair<IntRMMIter, bool>:
1690 // use iterator to shape RealMatrix in place and init sums to 0
1691 sum_Ql.insert(empty_irm_pr).first->second.shape(numFunctions, num_lev);
1692 sum_Qlm1.insert(empty_irm_pr).first->second.shape(numFunctions, num_lev);
1693 }
1694 std::pair<IntIntPair, RealMatrix> empty_iirm_pr;
1695 IntIntPair& ii = empty_iirm_pr.first;
1696 for (i=1; i<=2; ++i)
1697 for (j=1; j<=2; ++j) {
1698 ii.first = i; ii.second = j;
1699 // std::map::insert() returns std::pair<IIPRMMap::iterator, bool>
1700 sum_QlQlm1.insert(empty_iirm_pr).first->
1701 second.shape(numFunctions, num_lev);
1702 }
1703 }
1704
1705
1706 void NonDMultilevelSampling::
initialize_cv_sums(IntRealVectorMap & sum_L_shared,IntRealVectorMap & sum_L_refined,IntRealVectorMap & sum_H,IntRealVectorMap & sum_LL,IntRealVectorMap & sum_LH)1707 initialize_cv_sums(IntRealVectorMap& sum_L_shared,
1708 IntRealVectorMap& sum_L_refined, IntRealVectorMap& sum_H,
1709 IntRealVectorMap& sum_LL, //IntRealVectorMap& sum_HH,
1710 IntRealVectorMap& sum_LH)
1711 {
1712 // sum_* are running sums across all increments
1713 std::pair<int, RealVector> empty_pr;
1714 for (int i=1; i<=4; ++i) {
1715 empty_pr.first = i;
1716 // std::map::insert() returns std::pair<IntRVMIter, bool>:
1717 // use iterator to size RealVector in place and init sums to 0
1718 sum_L_shared.insert(empty_pr).first->second.size(numFunctions);
1719 sum_L_refined.insert(empty_pr).first->second.size(numFunctions);
1720 sum_H.insert(empty_pr).first->second.size(numFunctions);
1721 sum_LL.insert(empty_pr).first->second.size(numFunctions);
1722 //sum_HH.insert(empty_pr).first->second.size(numFunctions);
1723 sum_LH.insert(empty_pr).first->second.size(numFunctions);
1724 }
1725 }
1726
1727
1728 void NonDMultilevelSampling::
initialize_mlcv_sums(IntRealMatrixMap & sum_L_shared,IntRealMatrixMap & sum_L_refined,IntRealMatrixMap & sum_H,IntRealMatrixMap & sum_LL,IntRealMatrixMap & sum_LH,IntRealMatrixMap & sum_HH,size_t num_ml_lev,size_t num_cv_lev)1729 initialize_mlcv_sums(IntRealMatrixMap& sum_L_shared,
1730 IntRealMatrixMap& sum_L_refined, IntRealMatrixMap& sum_H,
1731 IntRealMatrixMap& sum_LL, IntRealMatrixMap& sum_LH,
1732 IntRealMatrixMap& sum_HH, size_t num_ml_lev,
1733 size_t num_cv_lev)
1734 {
1735 // sum_* are running sums across all increments
1736 std::pair<int, RealMatrix> empty_pr;
1737 for (int i=1; i<=4; ++i) {
1738 empty_pr.first = i;
1739 // std::map::insert() returns std::pair<IntRVMIter, bool>:
1740 // use iterator to shape RealMatrix in place and init sums to 0
1741
1742 // num_cv_lev:
1743 sum_L_shared.insert(empty_pr).first->second.shape(numFunctions, num_cv_lev);
1744 sum_L_refined.insert(empty_pr).first->second.shape(numFunctions,num_cv_lev);
1745 sum_LL.insert(empty_pr).first->second.shape(numFunctions, num_cv_lev);
1746 sum_LH.insert(empty_pr).first->second.shape(numFunctions, num_cv_lev);
1747
1748 // num_ml_lev:
1749 sum_H.insert(empty_pr).first->second.shape(numFunctions, num_ml_lev);
1750 }
1751
1752 // Only need order 1 accumulation for HH
1753 empty_pr.first = 1;
1754 sum_HH.insert(empty_pr).first->second.shape(numFunctions, num_ml_lev);
1755 }
1756
1757
1758 void NonDMultilevelSampling::
initialize_mlcv_sums(IntRealMatrixMap & sum_Ll,IntRealMatrixMap & sum_Llm1,IntRealMatrixMap & sum_Ll_refined,IntRealMatrixMap & sum_Llm1_refined,IntRealMatrixMap & sum_Hl,IntRealMatrixMap & sum_Hlm1,IntRealMatrixMap & sum_Ll_Ll,IntRealMatrixMap & sum_Ll_Llm1,IntRealMatrixMap & sum_Llm1_Llm1,IntRealMatrixMap & sum_Hl_Ll,IntRealMatrixMap & sum_Hl_Llm1,IntRealMatrixMap & sum_Hlm1_Ll,IntRealMatrixMap & sum_Hlm1_Llm1,IntRealMatrixMap & sum_Hl_Hl,IntRealMatrixMap & sum_Hl_Hlm1,IntRealMatrixMap & sum_Hlm1_Hlm1,size_t num_ml_lev,size_t num_cv_lev)1759 initialize_mlcv_sums(IntRealMatrixMap& sum_Ll, IntRealMatrixMap& sum_Llm1,
1760 IntRealMatrixMap& sum_Ll_refined,
1761 IntRealMatrixMap& sum_Llm1_refined,
1762 IntRealMatrixMap& sum_Hl, IntRealMatrixMap& sum_Hlm1,
1763 IntRealMatrixMap& sum_Ll_Ll, IntRealMatrixMap& sum_Ll_Llm1,
1764 IntRealMatrixMap& sum_Llm1_Llm1,
1765 IntRealMatrixMap& sum_Hl_Ll, IntRealMatrixMap& sum_Hl_Llm1,
1766 IntRealMatrixMap& sum_Hlm1_Ll,
1767 IntRealMatrixMap& sum_Hlm1_Llm1,
1768 IntRealMatrixMap& sum_Hl_Hl, IntRealMatrixMap& sum_Hl_Hlm1,
1769 IntRealMatrixMap& sum_Hlm1_Hlm1, size_t num_ml_lev,
1770 size_t num_cv_lev)
1771 {
1772 // sum_* are running sums across all increments
1773 std::pair<int, RealMatrix> empty_pr;
1774 for (int i=1; i<=4; ++i) {
1775 empty_pr.first = i;
1776 // std::map::insert() returns std::pair<IntRVMIter, bool>:
1777 // use iterator to shape RealMatrix in place and init sums to 0
1778
1779 // num_cv_lev:
1780 sum_Ll.insert(empty_pr).first->second.shape(numFunctions, num_cv_lev);
1781 sum_Llm1.insert(empty_pr).first->second.shape(numFunctions, num_cv_lev);
1782 sum_Ll_refined.insert(empty_pr).first->
1783 second.shape(numFunctions, num_cv_lev);
1784 sum_Llm1_refined.insert(empty_pr).first->
1785 second.shape(numFunctions, num_cv_lev);
1786 sum_Hlm1.insert(empty_pr).first->second.shape(numFunctions, num_cv_lev);
1787 sum_Ll_Ll.insert(empty_pr).first->second.shape(numFunctions, num_cv_lev);
1788 sum_Ll_Llm1.insert(empty_pr).first->second.shape(numFunctions, num_cv_lev);
1789 sum_Llm1_Llm1.insert(empty_pr).first->second.shape(numFunctions,num_cv_lev);
1790 sum_Hl_Ll.insert(empty_pr).first->second.shape(numFunctions, num_cv_lev);
1791 sum_Hl_Llm1.insert(empty_pr).first->second.shape(numFunctions, num_cv_lev);
1792 sum_Hlm1_Ll.insert(empty_pr).first->second.shape(numFunctions, num_cv_lev);
1793 sum_Hlm1_Llm1.insert(empty_pr).first->second.shape(numFunctions,num_cv_lev);
1794 // num_ml_lev:
1795 sum_Hl.insert(empty_pr).first->second.shape(numFunctions, num_ml_lev);
1796 sum_Hl_Hl.insert(empty_pr).first->second.shape(numFunctions, num_ml_lev);
1797 sum_Hl_Hlm1.insert(empty_pr).first->second.shape(numFunctions, num_ml_lev);
1798 sum_Hlm1_Hlm1.insert(empty_pr).first->second.shape(numFunctions,num_ml_lev);
1799 }
1800 }
1801
1802
1803 void NonDMultilevelSampling::
accumulate_ml_Qsums(IntRealMatrixMap & sum_Q,size_t lev,const RealVector & offset,SizetArray & num_Q)1804 accumulate_ml_Qsums(IntRealMatrixMap& sum_Q, size_t lev,
1805 const RealVector& offset, SizetArray& num_Q)
1806 {
1807 using std::isfinite;
1808 Real q_l, q_l_prod;
1809 int ord, active_ord; size_t qoi;
1810 IntRespMCIter r_it; IntRMMIter q_it;
1811 bool os = !offset.empty();
1812
1813 for (r_it=allResponses.begin(); r_it!=allResponses.end(); ++r_it) {
1814 const RealVector& fn_vals = r_it->second.function_values();
1815
1816 for (qoi=0; qoi<numFunctions; ++qoi) {
1817 q_l_prod = q_l = (os) ? fn_vals[qoi] - offset[qoi] : fn_vals[qoi];
1818
1819 if (isfinite(q_l)) { // neither NaN nor +/-Inf
1820 q_it = sum_Q.begin(); ord = q_it->first;
1821 active_ord = 1;
1822 while (q_it!=sum_Q.end()) {
1823
1824 if (ord == active_ord) {
1825 q_it->second(qoi,lev) += q_l_prod; ++q_it;
1826 ord = (q_it == sum_Q.end()) ? 0 : q_it->first;
1827 }
1828
1829 q_l_prod *= q_l; ++active_ord;
1830 }
1831 ++num_Q[qoi];
1832 }
1833 }
1834 //Cout << r_it->first << ": " << sum_Q[1];
1835 }
1836 }
1837
1838
1839 void NonDMultilevelSampling::
accumulate_ml_Qsums(IntRealMatrixMap & sum_Ql,IntRealMatrixMap & sum_Qlm1,IntIntPairRealMatrixMap & sum_QlQlm1,size_t lev,const RealVector & offset,SizetArray & num_Q)1840 accumulate_ml_Qsums(IntRealMatrixMap& sum_Ql, IntRealMatrixMap& sum_Qlm1,
1841 IntIntPairRealMatrixMap& sum_QlQlm1, size_t lev,
1842 const RealVector& offset, SizetArray& num_Q)
1843 {
1844 if (lev == 0)
1845 accumulate_ml_Qsums(sum_Ql, lev, offset, num_Q);
1846 else {
1847 using std::isfinite;
1848 Real q_l, q_lm1, q_l_prod, q_lm1_prod, qq_prod;
1849 int l1_ord, l2_ord, active_ord; size_t qoi;
1850 IntRespMCIter r_it; IntRMMIter l1_it, l2_it; IntIntPair pr;
1851 bool os = !offset.empty();
1852
1853 for (r_it=allResponses.begin(); r_it!=allResponses.end(); ++r_it) {
1854 const RealVector& fn_vals = r_it->second.function_values();
1855
1856 for (qoi=0; qoi<numFunctions; ++qoi) {
1857 // response mode AGGREGATED_MODELS orders HF (active model key)
1858 // followed by LF (previous/decremented model key)
1859 q_l_prod = q_l = (os) ? fn_vals[qoi] - offset[qoi] : fn_vals[qoi];
1860 q_lm1_prod = q_lm1 = (os) ?
1861 fn_vals[qoi+numFunctions] - offset[qoi+numFunctions] :
1862 fn_vals[qoi+numFunctions];
1863
1864 // sync sample counts for Ql and Qlm1
1865 if (isfinite(q_l) && isfinite(q_lm1)) { // neither NaN nor +/-Inf
1866
1867 // covariance terms: products of q_l and q_lm1
1868 pr.first = pr.second = 1;
1869 qq_prod = q_l_prod * q_lm1_prod;
1870 sum_QlQlm1[pr](qoi,lev) += qq_prod;
1871 pr.second = 2;
1872 sum_QlQlm1[pr](qoi,lev) += qq_prod * q_lm1_prod;
1873 pr.first = 2; pr.second = 1;
1874 qq_prod *= q_l_prod;
1875 sum_QlQlm1[pr](qoi,lev) += qq_prod;
1876 pr.second = 2;
1877 sum_QlQlm1[pr](qoi,lev) += qq_prod * q_lm1_prod;
1878
1879 // mean,variance terms: products of q_l or products of q_lm1
1880 l1_it = sum_Ql.begin();
1881 l2_it = sum_Qlm1.begin();
1882 l1_ord = (l1_it == sum_Ql.end()) ? 0 : l1_it->first;
1883 l2_ord = (l2_it == sum_Qlm1.end()) ? 0 : l2_it->first;
1884 active_ord = 1;
1885 while (l1_it != sum_Ql.end() || l2_it != sum_Qlm1.end()) {
1886
1887 // Low: Ll, Llm1
1888 if (l1_ord == active_ord) {
1889 l1_it->second(qoi,lev) += q_l_prod; ++l1_it;
1890 l1_ord = (l1_it == sum_Ql.end()) ? 0 : l1_it->first;
1891 }
1892 if (l2_ord == active_ord) {
1893 l2_it->second(qoi,lev) += q_lm1_prod; ++l2_it;
1894 l2_ord = (l2_it == sum_Qlm1.end()) ? 0 : l2_it->first;
1895 }
1896
1897 q_l_prod *= q_l; q_lm1_prod *= q_lm1; ++active_ord;
1898 }
1899 ++num_Q[qoi];
1900 }
1901 }
1902 //Cout << r_it->first << ": " << sum_Ql[1] << sum_Qlm1[1];
1903 }
1904 }
1905 }
1906
1907
1908 void NonDMultilevelSampling::
accumulate_ml_Ysums(IntRealMatrixMap & sum_Y,RealMatrix & sum_YY,size_t lev,const RealVector & offset,SizetArray & num_Y)1909 accumulate_ml_Ysums(IntRealMatrixMap& sum_Y, RealMatrix& sum_YY, size_t lev,
1910 const RealVector& offset, SizetArray& num_Y)
1911 {
1912 using std::isfinite;
1913 Real lf_fn, lf_prod;
1914 int y_ord, active_ord; size_t qoi;
1915 IntRespMCIter r_it; IntRMMIter y_it;
1916 bool os = !offset.empty();
1917
1918 if (lev == 0) {
1919 for (r_it=allResponses.begin(); r_it!=allResponses.end(); ++r_it) {
1920 const RealVector& fn_vals = r_it->second.function_values();
1921 for (qoi=0; qoi<numFunctions; ++qoi) {
1922
1923 lf_prod = lf_fn = (os) ? fn_vals[qoi] - offset[qoi] : fn_vals[qoi];
1924 if (isfinite(lf_fn)) { // neither NaN nor +/-Inf
1925 // add to sum_YY: running sums across all sample increments
1926 sum_YY(qoi,lev) += lf_prod * lf_prod;
1927
1928 // add to sum_Y: running sums across all sample increments
1929 y_it = sum_Y.begin(); y_ord = y_it->first;
1930 active_ord = 1;
1931 while (y_it!=sum_Y.end() || active_ord <= 1) {
1932 if (y_ord == active_ord) {
1933 y_it->second(qoi,lev) += lf_prod; ++y_it;
1934 y_ord = (y_it == sum_Y.end()) ? 0 : y_it->first;
1935 }
1936 lf_prod *= lf_fn; ++active_ord;
1937 }
1938 ++num_Y[qoi];
1939 }
1940 }
1941 }
1942 }
1943 else {
1944 Real hf_fn, hf_prod;
1945 for (r_it=allResponses.begin(); r_it!=allResponses.end(); ++r_it) {
1946 const RealVector& fn_vals = r_it->second.function_values();
1947 for (qoi=0; qoi<numFunctions; ++qoi) {
1948
1949 // response mode AGGREGATED_MODELS orders HF (active model key)
1950 // followed by LF (previous/decremented model key)
1951 hf_prod = hf_fn = (os) ? fn_vals[qoi] - offset[qoi] : fn_vals[qoi];
1952 lf_prod = lf_fn = (os) ?
1953 fn_vals[qoi+numFunctions] - offset[qoi+numFunctions] :
1954 fn_vals[qoi+numFunctions];
1955 if (isfinite(lf_fn) && isfinite(hf_fn)) { // neither NaN nor +/-Inf
1956
1957 // add to sum_YY: running sums across all sample increments
1958 Real delta_prod = hf_prod - lf_prod;
1959 sum_YY(qoi,lev) += delta_prod * delta_prod; // (HF^p-LF^p)^2 for p=1
1960
1961 // add to sum_Y: running sums across all sample increments
1962 y_it = sum_Y.begin(); y_ord = y_it->first;
1963 active_ord = 1;
1964 while (y_it!=sum_Y.end() || active_ord <= 1) {
1965 if (y_ord == active_ord) {
1966 y_it->second(qoi,lev) += hf_prod - lf_prod; // HF^p-LF^p
1967 ++y_it; y_ord = (y_it == sum_Y.end()) ? 0 : y_it->first;
1968 }
1969 hf_prod *= hf_fn; lf_prod *= lf_fn; ++active_ord;
1970 }
1971 ++num_Y[qoi];
1972 }
1973 }
1974 }
1975 }
1976 }
1977
1978
1979 void NonDMultilevelSampling::
accumulate_cv_sums(IntRealVectorMap & sum_L,const RealVector & offset,SizetArray & num_L)1980 accumulate_cv_sums(IntRealVectorMap& sum_L, const RealVector& offset,
1981 SizetArray& num_L)
1982 {
1983 // uses one set of allResponses in BYPASS_SURROGATE mode
1984 // IntRealVectorMap is not a multilevel case --> no discrepancies
1985
1986 using std::isfinite;
1987 Real fn_val, prod;
1988 int ord, active_ord; size_t qoi;
1989 IntRespMCIter r_it; IntRVMIter l_it;
1990 bool os = !offset.empty();
1991
1992 for (r_it=allResponses.begin(); r_it!=allResponses.end(); ++r_it) {
1993 const RealVector& fn_vals = r_it->second.function_values();
1994
1995 for (qoi=0; qoi<numFunctions; ++qoi) {
1996 prod = fn_val = (os) ? fn_vals[qoi] - offset[qoi] : fn_vals[qoi];
1997
1998 if (isfinite(fn_val)) { // neither NaN nor +/-Inf
1999 l_it = sum_L.begin(); ord = l_it->first; active_ord = 1;
2000 while (l_it!=sum_L.end()) {
2001
2002 if (ord == active_ord) {
2003 l_it->second[qoi] += prod; ++l_it;
2004 ord = (l_it == sum_L.end()) ? 0 : l_it->first;
2005 }
2006
2007 prod *= fn_val; ++active_ord;
2008 }
2009 ++num_L[qoi];
2010 }
2011 }
2012 }
2013 }
2014
2015
2016 void NonDMultilevelSampling::
accumulate_cv_sums(IntRealVectorMap & sum_L_shared,IntRealVectorMap & sum_L_refined,IntRealVectorMap & sum_H,IntRealVectorMap & sum_LL,IntRealVectorMap & sum_LH,RealVector & sum_HH,const RealVector & offset,SizetArray & num_L,SizetArray & num_H)2017 accumulate_cv_sums(IntRealVectorMap& sum_L_shared,
2018 IntRealVectorMap& sum_L_refined, IntRealVectorMap& sum_H,
2019 IntRealVectorMap& sum_LL, IntRealVectorMap& sum_LH,
2020 RealVector& sum_HH, const RealVector& offset,
2021 SizetArray& num_L, SizetArray& num_H)
2022 {
2023 // uses one set of allResponses in AGGREGATED_MODELS mode
2024 // IntRealVectorMap is not a multilevel case so no discrepancies
2025
2026 using std::isfinite;
2027 Real lf_fn, hf_fn, lf_prod, hf_prod;
2028 IntRespMCIter r_it; IntRVMIter ls_it, lr_it, h_it, ll_it, lh_it;
2029 int ls_ord, lr_ord, h_ord, ll_ord, lh_ord, active_ord; size_t qoi;
2030 bool os = !offset.empty();
2031
2032 for (r_it=allResponses.begin(); r_it!=allResponses.end(); ++r_it) {
2033 const RealVector& fn_vals = r_it->second.function_values();
2034
2035 for (qoi=0; qoi<numFunctions; ++qoi) {
2036
2037 // response mode AGGREGATED_MODELS orders HF (active model key)
2038 // followed by LF (previous/decremented model key)
2039 hf_prod = hf_fn = (os) ? fn_vals[qoi] - offset[qoi] : fn_vals[qoi];
2040 lf_prod = lf_fn = (os) ?
2041 fn_vals[qoi+numFunctions] - offset[qoi+numFunctions] :
2042 fn_vals[qoi+numFunctions];
2043
2044 // sync sample counts for all L and H interactions at this level
2045 if (isfinite(lf_fn) && isfinite(hf_fn)) { // neither NaN nor +/-Inf
2046
2047 // High-High
2048 sum_HH[qoi] += hf_prod * hf_prod;
2049
2050 ls_it = sum_L_shared.begin(); lr_it = sum_L_refined.begin();
2051 h_it = sum_H.begin(); ll_it = sum_LL.begin(); lh_it = sum_LH.begin();
2052 ls_ord = /*(ls_it == sum_L_shared.end()) ? 0 :*/ ls_it->first;
2053 lr_ord = /*(lr_it == sum_L_refined.end()) ? 0 :*/ lr_it->first;
2054 h_ord = /*(h_it == sum_H.end()) ? 0 :*/ h_it->first;
2055 ll_ord = /*(ll_it == sum_LL.end()) ? 0 :*/ ll_it->first;
2056 lh_ord = /*(lh_it == sum_LH.end()) ? 0 :*/ lh_it->first;
2057 active_ord = 1;
2058 while (ls_it!=sum_L_shared.end() || lr_it!=sum_L_refined.end() ||
2059 h_it!=sum_H.end() || ll_it!=sum_LL.end() ||
2060 lh_it!=sum_LH.end() || active_ord <= 1) {
2061
2062 // Low shared
2063 if (ls_ord == active_ord) {
2064 ls_it->second[qoi] += lf_prod;
2065 ++ls_it; ls_ord = (ls_it == sum_L_shared.end()) ? 0 : ls_it->first;
2066 }
2067 // Low refined
2068 if (lr_ord == active_ord) {
2069 lr_it->second[qoi] += lf_prod;
2070 ++lr_it; lr_ord = (lr_it == sum_L_refined.end()) ? 0 : lr_it->first;
2071 }
2072 // High
2073 if (h_ord == active_ord) {
2074 h_it->second[qoi] += hf_prod;
2075 ++h_it; h_ord = (h_it == sum_H.end()) ? 0 : h_it->first;
2076 }
2077 // Low-Low
2078 if (ll_ord == active_ord) {
2079 ll_it->second[qoi] += lf_prod * lf_prod;
2080 ++ll_it; ll_ord = (ll_it == sum_LL.end()) ? 0 : ll_it->first;
2081 }
2082 // Low-High
2083 if (lh_ord == active_ord) {
2084 lh_it->second[qoi] += lf_prod * hf_prod;
2085 ++lh_it; lh_ord = (lh_it == sum_LH.end()) ? 0 : lh_it->first;
2086 }
2087
2088 if (ls_ord || lr_ord || ll_ord || lh_ord) lf_prod *= lf_fn;
2089 if (h_ord || lh_ord) hf_prod *= hf_fn;
2090 ++active_ord;
2091 }
2092 ++num_L[qoi]; ++num_H[qoi];
2093 }
2094 }
2095 }
2096 }
2097
2098
2099 void NonDMultilevelSampling::
accumulate_mlcv_Qsums(IntRealMatrixMap & sum_Ql,IntRealMatrixMap & sum_Qlm1,size_t lev,const RealVector & offset,SizetArray & num_Q)2100 accumulate_mlcv_Qsums(IntRealMatrixMap& sum_Ql, IntRealMatrixMap& sum_Qlm1,
2101 size_t lev, const RealVector& offset, SizetArray& num_Q)
2102 {
2103 if (lev == 0)
2104 accumulate_ml_Qsums(sum_Ql, lev, offset, num_Q);
2105 else {
2106 using std::isfinite;
2107 Real q_l, q_l_prod, q_lm1_prod, q_lm1;
2108 int l1_ord, l2_ord, active_ord; size_t qoi;
2109 IntRespMCIter r_it; IntRMMIter l1_it, l2_it;
2110 bool os = !offset.empty();
2111
2112 for (r_it=allResponses.begin(); r_it!=allResponses.end(); ++r_it) {
2113 const RealVector& fn_vals = r_it->second.function_values();
2114
2115 for (qoi=0; qoi<numFunctions; ++qoi) {
2116 // response mode AGGREGATED_MODELS orders HF (active model key)
2117 // followed by LF (previous/decremented model key)
2118 q_l_prod = q_l = (os) ? fn_vals[qoi] - offset[qoi] : fn_vals[qoi];
2119 q_lm1_prod = q_lm1 = (os) ?
2120 fn_vals[qoi+numFunctions] - offset[qoi+numFunctions] :
2121 fn_vals[qoi+numFunctions];
2122
2123 // sync sample counts for Ql and Qlm1
2124 if (isfinite(q_l) && isfinite(q_lm1)) { // neither NaN nor +/-Inf
2125 l1_it = sum_Ql.begin();
2126 l2_it = sum_Qlm1.begin();
2127 l1_ord = (l1_it == sum_Ql.end()) ? 0 : l1_it->first;
2128 l2_ord = (l2_it == sum_Qlm1.end()) ? 0 : l2_it->first;
2129
2130 active_ord = 1;
2131 while (l1_it != sum_Ql.end() || l2_it != sum_Qlm1.end()) {
2132
2133 // Low: Ll, Llm1
2134 if (l1_ord == active_ord) {
2135 l1_it->second(qoi,lev) += q_l_prod; ++l1_it;
2136 l1_ord = (l1_it == sum_Ql.end()) ? 0 : l1_it->first;
2137 }
2138 if (l2_ord == active_ord) {
2139 l2_it->second(qoi,lev) += q_lm1_prod; ++l2_it;
2140 l2_ord = (l2_it == sum_Qlm1.end()) ? 0 : l2_it->first;
2141 }
2142
2143 q_l_prod *= q_l; q_lm1_prod *= q_lm1; ++active_ord;
2144 }
2145 ++num_Q[qoi];
2146 }
2147 }
2148 //Cout << r_it->first << ": " << sum_Ql[1] << sum_Qlm1[1];
2149 }
2150 }
2151 }
2152
2153
2154 void NonDMultilevelSampling::
accumulate_mlcv_Ysums(IntRealMatrixMap & sum_Y,size_t lev,const RealVector & offset,SizetArray & num_Y)2155 accumulate_mlcv_Ysums(IntRealMatrixMap& sum_Y, size_t lev,
2156 const RealVector& offset, SizetArray& num_Y)
2157 {
2158 // uses one set of allResponses in BYPASS_SURROGATE (level 0) or
2159 // AGGREGATED_MODELS (lev > 0) modes. IntRealMatrixMap is a multilevel
2160 // case with discrepancies, indexed by level.
2161
2162 if (lev == 0)
2163 accumulate_ml_Qsums(sum_Y, lev, offset, num_Y);
2164 else { // AGGREGATED_MODELS -> 2 sets of qoi per response map
2165 using std::isfinite;
2166 Real fn_l, prod_l, fn_lm1, prod_lm1;
2167 int ord, active_ord; size_t qoi;
2168 IntRespMCIter r_it; IntRMMIter y_it;
2169 bool os = !offset.empty();
2170
2171 for (r_it=allResponses.begin(); r_it!=allResponses.end(); ++r_it) {
2172 const RealVector& fn_vals = r_it->second.function_values();
2173
2174 for (qoi=0; qoi<numFunctions; ++qoi) {
2175 // response mode AGGREGATED_MODELS orders HF (active model key)
2176 // followed by LF (previous/decremented model key)
2177 prod_l = fn_l = (os) ? fn_vals[qoi] - offset[qoi] : fn_vals[qoi];
2178 prod_lm1 = fn_lm1 = (os) ?
2179 fn_vals[qoi+numFunctions] - offset[qoi+numFunctions] :
2180 fn_vals[qoi+numFunctions];
2181
2182 if (isfinite(fn_l) && isfinite(fn_lm1)) { // neither NaN nor +/-Inf
2183 y_it = sum_Y.begin(); ord = y_it->first; active_ord = 1;
2184 while (y_it!=sum_Y.end()) {
2185
2186 if (ord == active_ord) {
2187 y_it->second(qoi,lev) += prod_l - prod_lm1; ++y_it;
2188 ord = (y_it == sum_Y.end()) ? 0 : y_it->first;
2189 }
2190
2191 prod_l *= fn_l; prod_lm1 *= fn_lm1;
2192 ++active_ord;
2193 }
2194 ++num_Y[qoi];
2195 }
2196 }
2197 //Cout << r_it->first << ": " << sum_Y[1];
2198 }
2199 }
2200 }
2201
2202
2203 void NonDMultilevelSampling::
accumulate_mlcv_Qsums(const IntResponseMap & lf_resp_map,const IntResponseMap & hf_resp_map,IntRealMatrixMap & sum_L_shared,IntRealMatrixMap & sum_L_refined,IntRealMatrixMap & sum_H,IntRealMatrixMap & sum_LL,IntRealMatrixMap & sum_LH,IntRealMatrixMap & sum_HH,size_t lev,const RealVector & lf_offset,const RealVector & hf_offset,SizetArray & num_L,SizetArray & num_H)2204 accumulate_mlcv_Qsums(const IntResponseMap& lf_resp_map,
2205 const IntResponseMap& hf_resp_map,
2206 IntRealMatrixMap& sum_L_shared,
2207 IntRealMatrixMap& sum_L_refined, IntRealMatrixMap& sum_H,
2208 IntRealMatrixMap& sum_LL, IntRealMatrixMap& sum_LH,
2209 IntRealMatrixMap& sum_HH, size_t lev,
2210 const RealVector& lf_offset, const RealVector& hf_offset,
2211 SizetArray& num_L, SizetArray& num_H)
2212 {
2213 using std::isfinite;
2214 Real lf_l, hf_l, lf_l_prod, hf_l_prod;
2215 IntRespMCIter lf_r_it, hf_r_it;
2216 IntRMMIter ls_it, lr_it, h_it, ll_it, lh_it, hh_it;
2217 int ls_ord, lr_ord, h_ord, ll_ord, lh_ord, hh_ord, active_ord;
2218 size_t qoi;
2219 bool lfos = !lf_offset.empty(), hfos = !hf_offset.empty();
2220
2221 for (lf_r_it =lf_resp_map.begin(), hf_r_it =hf_resp_map.begin();
2222 lf_r_it!=lf_resp_map.end() && hf_r_it!=hf_resp_map.end();
2223 ++lf_r_it, ++hf_r_it) {
2224 const RealVector& lf_fn_vals = lf_r_it->second.function_values();
2225 const RealVector& hf_fn_vals = hf_r_it->second.function_values();
2226
2227 for (qoi=0; qoi<numFunctions; ++qoi) {
2228
2229 lf_l_prod = lf_l = (lfos) ?
2230 lf_fn_vals[qoi] - lf_offset[qoi] : lf_fn_vals[qoi];
2231 hf_l_prod = hf_l = (hfos) ?
2232 hf_fn_vals[qoi] - hf_offset[qoi] : hf_fn_vals[qoi];
2233
2234 // sync sample counts for all L and H interactions at this level
2235 if (isfinite(lf_l) && isfinite(hf_l)) { // neither NaN nor +/-Inf
2236 ls_it = sum_L_shared.begin(); lr_it = sum_L_refined.begin();
2237 h_it = sum_H.begin(); ll_it = sum_LL.begin();
2238 lh_it = sum_LH.begin(); hh_it = sum_HH.begin();
2239 ls_ord = (ls_it == sum_L_shared.end()) ? 0 : ls_it->first;
2240 lr_ord = (lr_it == sum_L_refined.end()) ? 0 : lr_it->first;
2241 h_ord = (h_it == sum_H.end()) ? 0 : h_it->first;
2242 ll_ord = (ll_it == sum_LL.end()) ? 0 : ll_it->first;
2243 lh_ord = (lh_it == sum_LH.end()) ? 0 : lh_it->first;
2244 hh_ord = (hh_it == sum_HH.end()) ? 0 : hh_it->first;
2245 active_ord = 1;
2246
2247 while (ls_it!=sum_L_shared.end() || lr_it!=sum_L_refined.end() ||
2248 h_it !=sum_H.end() || ll_it!=sum_LL.end() ||
2249 lh_it!=sum_LH.end() || hh_it!=sum_HH.end()) {
2250
2251 // Low shared
2252 if (ls_ord == active_ord) {
2253 ls_it->second(qoi,lev) += lf_l_prod; ++ls_it;
2254 ls_ord = (ls_it == sum_L_shared.end()) ? 0 : ls_it->first;
2255 }
2256 // Low refined
2257 if (lr_ord == active_ord) {
2258 lr_it->second(qoi,lev) += lf_l_prod; ++lr_it;
2259 lr_ord = (lr_it == sum_L_refined.end()) ? 0 : lr_it->first;
2260 }
2261 // High
2262 if (h_ord == active_ord) {
2263 h_it->second(qoi,lev) += hf_l_prod; ++h_it;
2264 h_ord = (h_it == sum_H.end()) ? 0 : h_it->first;
2265 }
2266 // Low-Low
2267 if (ll_ord == active_ord) {
2268 ll_it->second(qoi,lev) += lf_l_prod * lf_l_prod; ++ll_it;
2269 ll_ord = (ll_it == sum_LL.end()) ? 0 : ll_it->first;
2270 }
2271 // Low-High
2272 if (lh_ord == active_ord) {
2273 lh_it->second(qoi,lev) += lf_l_prod * hf_l_prod; ++lh_it;
2274 lh_ord = (lh_it == sum_LH.end()) ? 0 : lh_it->first;
2275 }
2276 // High-High (no map, only a single matrix for order 1)
2277 if (hh_ord == active_ord) {
2278 hh_it->second(qoi,lev) += hf_l_prod * hf_l_prod; ++hh_it;
2279 hh_ord = (hh_it == sum_HH.end()) ? 0 : hh_it->first;
2280 }
2281
2282 if (ls_ord || lr_ord || ll_ord || lh_ord) lf_l_prod *= lf_l;
2283 if (h_ord || lh_ord || hh_ord) hf_l_prod *= hf_l;
2284 ++active_ord;
2285 }
2286 ++num_L[qoi]; ++num_H[qoi];
2287 }
2288 }
2289 }
2290 }
2291
2292
2293 void NonDMultilevelSampling::
accumulate_mlcv_Ysums(const IntResponseMap & lf_resp_map,const IntResponseMap & hf_resp_map,IntRealMatrixMap & sum_L_shared,IntRealMatrixMap & sum_L_refined,IntRealMatrixMap & sum_H,IntRealMatrixMap & sum_LL,IntRealMatrixMap & sum_LH,IntRealMatrixMap & sum_HH,size_t lev,const RealVector & lf_offset,const RealVector & hf_offset,SizetArray & num_L,SizetArray & num_H)2294 accumulate_mlcv_Ysums(const IntResponseMap& lf_resp_map,
2295 const IntResponseMap& hf_resp_map,
2296 IntRealMatrixMap& sum_L_shared,
2297 IntRealMatrixMap& sum_L_refined, IntRealMatrixMap& sum_H,
2298 IntRealMatrixMap& sum_LL, IntRealMatrixMap& sum_LH,
2299 IntRealMatrixMap& sum_HH, size_t lev,
2300 const RealVector& lf_offset, const RealVector& hf_offset,
2301 SizetArray& num_L, SizetArray& num_H)
2302 {
2303 // uses two sets of responses (LF & HF) in BYPASS_SURROGATE (level 0) or
2304 // AGGREGATED_MODELS (lev > 0) modes. IntRealMatrixMap are for multilevel
2305 // case with discrepancies, indexed by level.
2306
2307 if (lev == 0) // BYPASS_SURROGATE -> 1 set of qoi per response map
2308 accumulate_mlcv_Qsums(lf_resp_map, hf_resp_map, sum_L_shared,
2309 sum_L_refined, sum_H, sum_LL, sum_LH, sum_HH,
2310 lev, lf_offset, hf_offset, num_L, num_H);
2311 else { // AGGREGATED_MODELS -> 2 sets of qoi per response map
2312 using std::isfinite;
2313 Real lf_l, lf_l_prod, lf_lm1, lf_lm1_prod,
2314 hf_l, hf_l_prod, hf_lm1, hf_lm1_prod;
2315 IntRespMCIter lf_r_it, hf_r_it;
2316 IntRMMIter ls_it, lr_it, h_it, ll_it, lh_it, hh_it;
2317 int ls_ord, lr_ord, h_ord, ll_ord, lh_ord, hh_ord, active_ord;
2318 size_t qoi;
2319 bool lfos = !lf_offset.empty(), hfos = !hf_offset.empty();
2320
2321 for (lf_r_it =lf_resp_map.begin(), hf_r_it =hf_resp_map.begin();
2322 lf_r_it!=lf_resp_map.end() && hf_r_it!=hf_resp_map.end();
2323 ++lf_r_it, ++hf_r_it) {
2324 const RealVector& lf_fn_vals = lf_r_it->second.function_values();
2325 const RealVector& hf_fn_vals = hf_r_it->second.function_values();
2326
2327 for (qoi=0; qoi<numFunctions; ++qoi) {
2328
2329 // response mode AGGREGATED_MODELS orders level l (active model key)
2330 // followed by level l-1 (previous/decremented model key)
2331 lf_l_prod = lf_l = (lfos) ?
2332 lf_fn_vals[qoi] - lf_offset[qoi] : lf_fn_vals[qoi];
2333 lf_lm1_prod = lf_lm1 = (lfos) ?
2334 lf_fn_vals[qoi+numFunctions] - lf_offset[qoi+numFunctions] :
2335 lf_fn_vals[qoi+numFunctions];
2336 hf_l_prod = hf_l = (hfos) ?
2337 hf_fn_vals[qoi] - hf_offset[qoi] : hf_fn_vals[qoi];
2338 hf_lm1_prod = hf_lm1 = (hfos) ?
2339 hf_fn_vals[qoi+numFunctions] - hf_offset[qoi+numFunctions] :
2340 hf_fn_vals[qoi+numFunctions];
2341
2342 // sync sample counts for all L and H interactions at this level
2343 if (isfinite(lf_l) && isfinite(lf_lm1) &&
2344 isfinite(hf_l) && isfinite(hf_lm1)) { // neither NaN nor +/-Inf
2345 ls_it = sum_L_shared.begin(); lr_it = sum_L_refined.begin();
2346 h_it = sum_H.begin(); ll_it = sum_LL.begin();
2347 lh_it = sum_LH.begin(); hh_it = sum_HH.begin();
2348 ls_ord = (ls_it == sum_L_shared.end()) ? 0 : ls_it->first;
2349 lr_ord = (lr_it == sum_L_refined.end()) ? 0 : lr_it->first;
2350 h_ord = (h_it == sum_H.end()) ? 0 : h_it->first;
2351 ll_ord = (ll_it == sum_LL.end()) ? 0 : ll_it->first;
2352 lh_ord = (lh_it == sum_LH.end()) ? 0 : lh_it->first;
2353 hh_ord = (hh_it == sum_HH.end()) ? 0 : hh_it->first;
2354 active_ord = 1;
2355
2356 while (ls_it!=sum_L_shared.end() || lr_it!=sum_L_refined.end() ||
2357 h_it !=sum_H.end() || ll_it!=sum_LL.end() ||
2358 lh_it!=sum_LH.end() || hh_it!=sum_HH.end()) {
2359
2360 // Low shared
2361 if (ls_ord == active_ord) {
2362 ls_it->second(qoi,lev) += lf_l_prod - lf_lm1_prod; ++ls_it;
2363 ls_ord = (ls_it == sum_L_shared.end()) ? 0 : ls_it->first;
2364 }
2365 // Low refined
2366 if (lr_ord == active_ord) {
2367 lr_it->second(qoi,lev) += lf_l_prod - lf_lm1_prod; ++lr_it;
2368 lr_ord = (lr_it == sum_L_refined.end()) ? 0 : lr_it->first;
2369 }
2370 // High
2371 if (h_ord == active_ord) {
2372 h_it->second(qoi,lev) += hf_l_prod - hf_lm1_prod; ++h_it;
2373 h_ord = (h_it == sum_H.end()) ? 0 : h_it->first;
2374 }
2375 // Low-Low
2376 if (ll_ord == active_ord) {
2377 Real delta_prod = lf_l_prod - lf_lm1_prod;
2378 ll_it->second(qoi,lev) += delta_prod * delta_prod;
2379 ++ll_it; ll_ord = (ll_it == sum_LL.end()) ? 0 : ll_it->first;
2380 }
2381 // Low-High
2382 if (lh_ord == active_ord) {
2383 lh_it->second(qoi,lev) += (lf_l_prod - lf_lm1_prod)
2384 * (hf_l_prod - hf_lm1_prod);
2385 ++lh_it; lh_ord = (lh_it == sum_LH.end()) ? 0 : lh_it->first;
2386 }
2387 // High-High (map only contains order 1 in some contexts)
2388 if (hh_ord == active_ord) {
2389 Real delta_prod = hf_l_prod - hf_lm1_prod;
2390 hh_it->second(qoi,lev) += delta_prod * delta_prod;
2391 ++hh_it; hh_ord = (hh_it == sum_HH.end()) ? 0 : hh_it->first;
2392 }
2393
2394 if (ls_ord || lr_ord || ll_ord || lh_ord)
2395 { lf_l_prod *= lf_l; lf_lm1_prod *= lf_lm1; }
2396 if (h_ord || lh_ord || hh_ord)
2397 { hf_l_prod *= hf_l; hf_lm1_prod *= hf_lm1; }
2398 ++active_ord;
2399 }
2400 ++num_L[qoi]; ++num_H[qoi];
2401 }
2402 }
2403 }
2404 }
2405 }
2406
2407
2408 void NonDMultilevelSampling::
accumulate_mlcv_Qsums(const IntResponseMap & lf_resp_map,const IntResponseMap & hf_resp_map,IntRealMatrixMap & sum_Ll,IntRealMatrixMap & sum_Llm1,IntRealMatrixMap & sum_Ll_refined,IntRealMatrixMap & sum_Llm1_refined,IntRealMatrixMap & sum_Hl,IntRealMatrixMap & sum_Hlm1,IntRealMatrixMap & sum_Ll_Ll,IntRealMatrixMap & sum_Ll_Llm1,IntRealMatrixMap & sum_Llm1_Llm1,IntRealMatrixMap & sum_Hl_Ll,IntRealMatrixMap & sum_Hl_Llm1,IntRealMatrixMap & sum_Hlm1_Ll,IntRealMatrixMap & sum_Hlm1_Llm1,IntRealMatrixMap & sum_Hl_Hl,IntRealMatrixMap & sum_Hl_Hlm1,IntRealMatrixMap & sum_Hlm1_Hlm1,size_t lev,const RealVector & lf_offset,const RealVector & hf_offset,SizetArray & num_L,SizetArray & num_H)2409 accumulate_mlcv_Qsums(const IntResponseMap& lf_resp_map,
2410 const IntResponseMap& hf_resp_map,
2411 IntRealMatrixMap& sum_Ll, IntRealMatrixMap& sum_Llm1,
2412 IntRealMatrixMap& sum_Ll_refined,
2413 IntRealMatrixMap& sum_Llm1_refined,
2414 IntRealMatrixMap& sum_Hl, IntRealMatrixMap& sum_Hlm1,
2415 IntRealMatrixMap& sum_Ll_Ll,
2416 IntRealMatrixMap& sum_Ll_Llm1,
2417 IntRealMatrixMap& sum_Llm1_Llm1,
2418 IntRealMatrixMap& sum_Hl_Ll,
2419 IntRealMatrixMap& sum_Hl_Llm1,
2420 IntRealMatrixMap& sum_Hlm1_Ll,
2421 IntRealMatrixMap& sum_Hlm1_Llm1,
2422 IntRealMatrixMap& sum_Hl_Hl,
2423 IntRealMatrixMap& sum_Hl_Hlm1,
2424 IntRealMatrixMap& sum_Hlm1_Hlm1, size_t lev,
2425 const RealVector& lf_offset, const RealVector& hf_offset,
2426 SizetArray& num_L, SizetArray& num_H)
2427 {
2428 // uses two sets of responses (LF & HF) in BYPASS_SURROGATE (level 0) or
2429 // AGGREGATED_MODELS (lev > 0) modes. IntRealMatrixMap are for multilevel
2430 // case with discrepancies, indexed by level.
2431
2432 if (lev == 0) // level lm1 not available; accumulate only level l
2433 accumulate_mlcv_Qsums(lf_resp_map, hf_resp_map, sum_Ll, sum_Ll_refined,
2434 sum_Hl, sum_Ll_Ll, sum_Hl_Ll, sum_Hl_Hl, lev,
2435 lf_offset, hf_offset, num_L, num_H);
2436 else {
2437 using std::isfinite;
2438 Real lf_l_prod, lf_l, lf_lm1_prod, lf_lm1,
2439 hf_l_prod, hf_l, hf_lm1_prod, hf_lm1;
2440 IntRespMCIter lf_r_it, hf_r_it;
2441 IntRMMIter l1_it, l2_it, lr1_it, lr2_it, h1_it, h2_it, ll1_it, ll2_it,
2442 ll3_it, lh1_it, lh2_it, lh3_it, lh4_it, hh1_it, hh2_it, hh3_it;
2443 int active_ord, l1_ord, l2_ord, lr1_ord, lr2_ord, h1_ord, h2_ord,
2444 ll1_ord, ll2_ord, ll3_ord, lh1_ord, lh2_ord, lh3_ord, lh4_ord,
2445 hh1_ord, hh2_ord, hh3_ord;
2446 size_t qoi;
2447 bool lfos = !lf_offset.empty(), hfos = !hf_offset.empty();
2448
2449 for (lf_r_it =lf_resp_map.begin(), hf_r_it =hf_resp_map.begin();
2450 lf_r_it!=lf_resp_map.end() && hf_r_it!=hf_resp_map.end();
2451 ++lf_r_it, ++hf_r_it) {
2452 const RealVector& lf_fn_vals = lf_r_it->second.function_values();
2453 const RealVector& hf_fn_vals = hf_r_it->second.function_values();
2454
2455 for (qoi=0; qoi<numFunctions; ++qoi) {
2456
2457 // response mode AGGREGATED_MODELS orders level l (active model key)
2458 // followed by level l-1 (previous/decremented model key)
2459 lf_l_prod = lf_l = (lfos) ?
2460 lf_fn_vals[qoi] - lf_offset[qoi] : lf_fn_vals[qoi];
2461 lf_lm1_prod = lf_lm1 = (lfos) ?
2462 lf_fn_vals[qoi+numFunctions] - lf_offset[qoi+numFunctions] :
2463 lf_fn_vals[qoi+numFunctions];
2464 hf_l_prod = hf_l = (hfos) ?
2465 hf_fn_vals[qoi] - hf_offset[qoi] : hf_fn_vals[qoi];
2466 hf_lm1_prod = hf_lm1 = (hfos) ?
2467 hf_fn_vals[qoi+numFunctions] - hf_offset[qoi+numFunctions] :
2468 hf_fn_vals[qoi+numFunctions];
2469
2470 // sync sample counts for all L and H interactions at this level
2471 if (isfinite(lf_l) && isfinite(lf_lm1) &&
2472 isfinite(hf_l) && isfinite(hf_lm1)) { // neither NaN nor +/-Inf
2473 // Low: Ll, Llm1, Ll_refined, Llm1_refined
2474 l1_it = sum_Ll.begin(); l2_it = sum_Llm1.begin();
2475 lr1_it = sum_Ll_refined.begin(); lr2_it = sum_Llm1_refined.begin();
2476 l1_ord = (l1_it == sum_Ll.end()) ? 0 : l1_it->first;
2477 l2_ord = (l2_it == sum_Llm1.end()) ? 0 : l2_it->first;
2478 lr1_ord = (lr1_it == sum_Ll_refined.end()) ? 0 : lr1_it->first;
2479 lr2_ord = (lr2_it == sum_Llm1_refined.end()) ? 0 : lr2_it->first;
2480 // High: Hl, Hlm1
2481 h1_it = sum_Hl.begin();
2482 h2_it = sum_Hlm1.begin();
2483 h1_ord = (h1_it == sum_Hl.end()) ? 0 : h1_it->first;
2484 h2_ord = (h2_it == sum_Hlm1.end()) ? 0 : h2_it->first;
2485 // Low-Low: Ll_Ll, Ll_Llm1, Llm1_Llm1
2486 ll1_it = sum_Ll_Ll.begin(); ll2_it = sum_Ll_Llm1.begin();
2487 ll3_it = sum_Llm1_Llm1.begin();
2488 ll1_ord = (ll1_it == sum_Ll_Ll.end()) ? 0 : ll1_it->first;
2489 ll2_ord = (ll2_it == sum_Ll_Llm1.end()) ? 0 : ll2_it->first;
2490 ll3_ord = (ll3_it == sum_Llm1_Llm1.end()) ? 0 : ll3_it->first;
2491 // Low-High: Hl_Ll, Hl_Llm1, Hlm1_Ll, Hlm1_Llm1
2492 lh1_it = sum_Hl_Ll.begin(); lh2_it = sum_Hl_Llm1.begin();
2493 lh3_it = sum_Hlm1_Ll.begin(); lh4_it = sum_Hlm1_Llm1.begin();
2494 lh1_ord = (lh1_it == sum_Hl_Ll.end()) ? 0 : lh1_it->first;
2495 lh2_ord = (lh2_it == sum_Hl_Llm1.end()) ? 0 : lh2_it->first;
2496 lh3_ord = (lh3_it == sum_Hlm1_Ll.end()) ? 0 : lh3_it->first;
2497 lh4_ord = (lh4_it == sum_Hlm1_Llm1.end()) ? 0 : lh4_it->first;
2498 // High-High: Hl_Hl, Hl_Hlm1, Hlm1_Hlm1
2499 hh1_it = sum_Hl_Hl.begin(); hh2_it = sum_Hl_Hlm1.begin();
2500 hh3_it = sum_Hlm1_Hlm1.begin();
2501 hh1_ord = (hh1_it == sum_Hl_Hl.end()) ? 0 : hh1_it->first;
2502 hh2_ord = (hh2_it == sum_Hl_Hlm1.end()) ? 0 : hh2_it->first;
2503 hh3_ord = (hh3_it == sum_Hlm1_Hlm1.end()) ? 0 : hh3_it->first;
2504
2505 active_ord = 1;
2506
2507 while (l1_it !=sum_Ll.end() || l2_it !=sum_Llm1.end() ||
2508 lr1_it!=sum_Ll_refined.end() ||
2509 lr2_it!=sum_Llm1_refined.end() || h1_it !=sum_Hl.end() ||
2510 h2_it !=sum_Hlm1.end() || ll1_it!=sum_Ll_Ll.end() ||
2511 ll2_it!=sum_Ll_Llm1.end() || ll3_it!=sum_Llm1_Llm1.end() ||
2512 lh1_it!=sum_Hl_Ll.end() || lh2_it!=sum_Hl_Llm1.end() ||
2513 lh3_it!=sum_Hlm1_Ll.end() || lh4_it!=sum_Hlm1_Llm1.end() ||
2514 hh1_it!=sum_Hl_Hl.end() || hh2_it!=sum_Hl_Hlm1.end() ||
2515 hh3_it!=sum_Hlm1_Hlm1.end()) {
2516
2517 // Low: Ll, Llm1, Ll_refined, Llm1_refined
2518 if (l1_ord == active_ord) {
2519 l1_it->second(qoi,lev) += lf_l_prod; ++l1_it;
2520 l1_ord = (l1_it == sum_Ll.end()) ? 0 : l1_it->first;
2521 }
2522 if (l2_ord == active_ord) {
2523 l2_it->second(qoi,lev) += lf_lm1_prod; ++l2_it;
2524 l2_ord = (l2_it == sum_Llm1.end()) ? 0 : l2_it->first;
2525 }
2526 if (lr1_ord == active_ord) {
2527 lr1_it->second(qoi,lev) += lf_l_prod; ++lr1_it;
2528 lr1_ord = (lr1_it == sum_Ll_refined.end()) ? 0 : lr1_it->first;
2529 }
2530 if (lr2_ord == active_ord) {
2531 lr2_it->second(qoi,lev) += lf_lm1_prod; ++lr2_it;
2532 lr2_ord = (lr2_it == sum_Llm1_refined.end()) ? 0 : lr2_it->first;
2533 }
2534 // High: Hl, Hlm1
2535 if (h1_ord == active_ord) {
2536 h1_it->second(qoi,lev) += hf_l_prod; ++h1_it;
2537 h1_ord = (h1_it == sum_Hl.end()) ? 0 : h1_it->first;
2538 }
2539 if (h2_ord == active_ord) {
2540 h2_it->second(qoi,lev) += hf_lm1_prod; ++h2_it;
2541 h2_ord = (h2_it == sum_Hlm1.end()) ? 0 : h2_it->first;
2542 }
2543 // Low-Low: Ll_Ll, Ll_Llm1, Llm1_Llm1
2544 if (ll1_ord == active_ord) {
2545 ll1_it->second(qoi,lev) += lf_l_prod * lf_l_prod; ++ll1_it;
2546 ll1_ord = (ll1_it == sum_Ll_Ll.end()) ? 0 : ll1_it->first;
2547 }
2548 if (ll2_ord == active_ord) {
2549 ll2_it->second(qoi,lev) += lf_l_prod * lf_lm1_prod; ++ll2_it;
2550 ll2_ord = (ll2_it == sum_Ll_Llm1.end()) ? 0 : ll2_it->first;
2551 }
2552 if (ll3_ord == active_ord) {
2553 ll3_it->second(qoi,lev) += lf_lm1_prod * lf_lm1_prod; ++ll3_it;
2554 ll3_ord = (ll3_it == sum_Llm1_Llm1.end()) ? 0 : ll3_it->first;
2555 }
2556 // Low-High: Hl_Ll, Hl_Llm1, Hlm1_Ll, Hlm1_Llm1
2557 if (lh1_ord == active_ord) {
2558 lh1_it->second(qoi,lev) += hf_l_prod * lf_l_prod; ++lh1_it;
2559 lh1_ord = (lh1_it == sum_Hl_Ll.end()) ? 0 : lh1_it->first;
2560 }
2561 if (lh2_ord == active_ord) {
2562 lh2_it->second(qoi,lev) += hf_l_prod * lf_lm1_prod; ++lh2_it;
2563 lh2_ord = (lh2_it == sum_Hl_Llm1.end()) ? 0 : lh2_it->first;
2564 }
2565 if (lh3_ord == active_ord) {
2566 lh3_it->second(qoi,lev) += hf_lm1_prod * lf_l_prod; ++lh3_it;
2567 lh3_ord = (lh3_it == sum_Hlm1_Ll.end()) ? 0 : lh3_it->first;
2568 }
2569 if (lh4_ord == active_ord) {
2570 lh4_it->second(qoi,lev) += hf_lm1_prod * lf_lm1_prod; ++lh4_it;
2571 lh4_ord = (lh4_it == sum_Hlm1_Llm1.end()) ? 0 : lh4_it->first;
2572 }
2573 // High-High: Hl_Hl, Hl_Hlm1, Hlm1_Hlm1
2574 if (hh1_ord == active_ord) {
2575 hh1_it->second(qoi,lev) += hf_l_prod * hf_l_prod; ++hh1_it;
2576 hh1_ord = (hh1_it == sum_Hl_Hl.end()) ? 0 : hh1_it->first;
2577 }
2578 if (hh2_ord == active_ord) {
2579 hh2_it->second(qoi,lev) += hf_l_prod * hf_lm1_prod; ++hh2_it;
2580 hh2_ord = (hh2_it == sum_Hl_Hlm1.end()) ? 0 : hh2_it->first;
2581 }
2582 if (hh3_ord == active_ord) {
2583 hh3_it->second(qoi,lev) += hf_lm1_prod * hf_lm1_prod; ++hh3_it;
2584 hh3_ord = (hh3_it == sum_Hlm1_Hlm1.end()) ? 0 : hh3_it->first;
2585 }
2586
2587 if (l1_ord || lr1_ord || ll1_ord || ll2_ord || lh1_ord || lh3_ord)
2588 lf_l_prod *= lf_l;
2589 if (l2_ord || lr2_ord || ll2_ord || ll3_ord || lh2_ord || lh4_ord)
2590 lf_lm1_prod *= lf_lm1;
2591 if (h1_ord || lh1_ord || lh2_ord || hh1_ord || hh2_ord)
2592 hf_l_prod *= hf_l;
2593 if (h2_ord || lh3_ord || lh4_ord || hh2_ord || hh3_ord)
2594 hf_lm1_prod *= hf_lm1;
2595 ++active_ord;
2596 }
2597 ++num_L[qoi]; ++num_H[qoi];
2598 }
2599 }
2600 }
2601 }
2602 }
2603
2604
shared_increment(size_t iter,size_t lev)2605 void NonDMultilevelSampling::shared_increment(size_t iter, size_t lev)
2606 {
2607 if (iter == _NPOS) Cout << "\nCVMC sample increments: ";
2608 else if (iter == 0) Cout << "\nCVMC pilot sample: ";
2609 else Cout << "\nCVMC iteration " << iter << " sample increments: ";
2610 Cout << "LF = " << numSamples << " HF = " << numSamples << '\n';
2611
2612 //aggregated_models_mode(); // set at calling level for CV
2613
2614 // generate new MC parameter sets
2615 get_parameter_sets(iteratedModel);// pull dist params from any model
2616
2617 // export separate output files for each data set:
2618 if (exportSampleSets) // for HF+LF models, use the HF tags
2619 export_all_samples("cv_", iteratedModel.truth_model(), iter, lev);
2620
2621 // compute allResponses from allVariables using hierarchical model
2622 evaluate_parameter_sets(iteratedModel, true, false);
2623 }
2624
2625
2626 bool NonDMultilevelSampling::
lf_increment(Real avg_eval_ratio,const SizetArray & N_lf,const SizetArray & N_hf,size_t iter,size_t lev)2627 lf_increment(Real avg_eval_ratio, const SizetArray& N_lf,
2628 const SizetArray& N_hf, size_t iter, size_t lev)
2629 {
2630 // ----------------------------------------------
2631 // Compute Final LF increment for control variate
2632 // ----------------------------------------------
2633
2634 // update LF samples based on evaluation ratio
2635 // r = m/n -> m = r*n -> delta = m-n = (r-1)*n
2636 // or with inverse r -> delta = m-n = n/inverse_r - n
2637
2638 numSamples = one_sided_delta(average(N_lf), average(N_hf) * avg_eval_ratio);
2639 if (numSamples) {
2640 Cout << "\nCVMC LF sample increment = " << numSamples;
2641 if (outputLevel >= DEBUG_OUTPUT)
2642 Cout << " from avg LF = " << average(N_lf) << ", avg HF = "
2643 << average(N_hf) << ", avg eval_ratio = " << avg_eval_ratio;
2644 Cout << std::endl;
2645
2646 // generate new MC parameter sets
2647 get_parameter_sets(iteratedModel);// pull dist params from any model
2648 // export separate output files for each data set:
2649 if (exportSampleSets)
2650 export_all_samples("cv_", iteratedModel.surrogate_model(), iter, lev);
2651
2652 // Iteration 0 is defined as the pilot sample, and each subsequent iter
2653 // can be defined as a CV increment followed by an ML increment. In this
2654 // case, terminating based on max_iterations results in a final ML increment
2655 // without a corresponding CV refinement; thus the number of ML and CV
2656 // refinements is consistent although the final sample profile is not
2657 // self-consistent -- to override this and finish with a final CV increment
2658 // corresponding to the final ML increment, the finalCVRefinement flag can
2659 // be set. Note: termination based on delta_N_hf=0 has a final ML increment
2660 // of zero and corresponding final CV increment of zero. Therefore, this
2661 // iteration completes on the previous CV increment and is more consistent
2662 // with finalCVRefinement=true.
2663 size_t max_iter = (maxIterations < 0) ? 25 : maxIterations; // default = -1
2664 if (iter < max_iter || finalCVRefinement) {
2665 // hierarchical surrogate mode could be BYPASS_SURROGATE for CV or
2666 // BYPASS_SURROGATE/AGGREGATED_MODELS for ML-CV
2667 //bypass_surrogate_mode(); // set at calling level for CV or ML-CV
2668
2669 // compute allResponses from allVariables using hierarchical model
2670 evaluate_parameter_sets(iteratedModel, true, false);
2671 return true;
2672 }
2673 else
2674 return false;
2675 }
2676 else {
2677 Cout << "\nNo CVMC LF sample increment";
2678 if (outputLevel >= DEBUG_OUTPUT)
2679 Cout << " from avg LF = " << average(N_lf) << ", avg HF = "
2680 << average(N_hf) << ", avg eval_ratio = " << avg_eval_ratio;
2681 Cout << std::endl;
2682 return false;
2683 }
2684 }
2685
2686
2687 Real NonDMultilevelSampling::
eval_ratio(const RealVector & sum_L_shared,const RealVector & sum_H,const RealVector & sum_LL,const RealVector & sum_LH,const RealVector & sum_HH,Real cost_ratio,const SizetArray & N_shared,RealVector & var_H,RealVector & rho2_LH)2688 eval_ratio(const RealVector& sum_L_shared, const RealVector& sum_H,
2689 const RealVector& sum_LL, const RealVector& sum_LH,
2690 const RealVector& sum_HH, Real cost_ratio,
2691 const SizetArray& N_shared, RealVector& var_H, RealVector& rho2_LH)
2692 {
2693 Real eval_ratio, avg_eval_ratio = 0.; size_t num_avg = 0;
2694 for (size_t qoi=0; qoi<numFunctions; ++qoi) {
2695
2696 Real& rho_sq = rho2_LH[qoi];
2697 compute_control(sum_L_shared[qoi], sum_H[qoi], sum_LL[qoi], sum_LH[qoi],
2698 sum_HH[qoi], N_shared[qoi], var_H[qoi], rho_sq);
2699
2700 // compute evaluation ratio which determines increment for LF samples
2701 // > the sample increment optimizes the total computational budget and is
2702 // not treated as a worst case accuracy reqmt --> use the QoI average
2703 // > refinement based only on QoI mean statistics
2704 // Given use of 1/r in MSE_ratio, one approach would average 1/r, but
2705 // this does not seem to behave as well in limited numerical experience.
2706 //if (rho_sq > Pecos::SMALL_NUMBER) {
2707 // avg_inv_eval_ratio += std::sqrt((1. - rho_sq)/(cost_ratio * rho_sq));
2708 if (rho_sq < 1.) { // protect against division by 0
2709 eval_ratio = std::sqrt(cost_ratio * rho_sq / (1. - rho_sq));
2710 if (outputLevel >= DEBUG_OUTPUT)
2711 Cout << "eval_ratio() QoI " << qoi+1 << ": cost_ratio = " << cost_ratio
2712 << " rho_sq = " << rho_sq << " eval_ratio = " << eval_ratio
2713 << std::endl;
2714 avg_eval_ratio += eval_ratio;
2715 ++num_avg;
2716 }
2717 }
2718 if (outputLevel >= DEBUG_OUTPUT)
2719 Cout << "variance of HF Q:\n" << var_H;
2720
2721 if (num_avg) avg_eval_ratio /= num_avg;
2722 else // should not happen, but provide a reasonable upper bound
2723 avg_eval_ratio = (Real)maxFunctionEvals / average(N_shared);
2724
2725 return avg_eval_ratio;
2726 }
2727
2728
2729 Real NonDMultilevelSampling::
eval_ratio(RealMatrix & sum_L_shared,RealMatrix & sum_H,RealMatrix & sum_LL,RealMatrix & sum_LH,RealMatrix & sum_HH,Real cost_ratio,size_t lev,const SizetArray & N_shared,RealMatrix & var_H,RealMatrix & rho2_LH)2730 eval_ratio(RealMatrix& sum_L_shared, RealMatrix& sum_H, RealMatrix& sum_LL,
2731 RealMatrix& sum_LH, RealMatrix& sum_HH, Real cost_ratio, size_t lev,
2732 const SizetArray& N_shared, RealMatrix& var_H, RealMatrix& rho2_LH)
2733 {
2734 Real eval_ratio, avg_eval_ratio = 0.; size_t num_avg = 0;
2735 for (size_t qoi=0; qoi<numFunctions; ++qoi) {
2736
2737 Real& rho_sq = rho2_LH(qoi,lev);
2738 compute_control(sum_L_shared(qoi,lev), sum_H(qoi,lev), sum_LL(qoi,lev),
2739 sum_LH(qoi,lev), sum_HH(qoi,lev), N_shared[qoi],
2740 var_H(qoi,lev), rho_sq);
2741
2742 if (rho_sq < 1.) { // protect against division by 0
2743 eval_ratio = std::sqrt(cost_ratio * rho_sq / (1. - rho_sq));
2744 if (outputLevel >= DEBUG_OUTPUT)
2745 Cout << "eval_ratio() QoI " << qoi+1 << ": cost_ratio = " << cost_ratio
2746 << " rho_sq = " << rho_sq << " eval_ratio = " << eval_ratio
2747 << std::endl;
2748 avg_eval_ratio += eval_ratio;
2749 ++num_avg;
2750 }
2751 }
2752 if (outputLevel >= DEBUG_OUTPUT) {
2753 Cout << "variance of HF Q[" << lev << "]:\n";
2754 write_col_vector_trans(Cout, (int)lev, (int)numFunctions, var_H);
2755 }
2756
2757 if (num_avg) avg_eval_ratio /= num_avg;
2758 else // should not happen, but provide a reasonable upper bound
2759 avg_eval_ratio = (Real)maxFunctionEvals / average(N_shared);
2760
2761 return avg_eval_ratio;
2762 }
2763
2764
2765 Real NonDMultilevelSampling::
eval_ratio(RealMatrix & sum_Ll,RealMatrix & sum_Llm1,RealMatrix & sum_Hl,RealMatrix & sum_Hlm1,RealMatrix & sum_Ll_Ll,RealMatrix & sum_Ll_Llm1,RealMatrix & sum_Llm1_Llm1,RealMatrix & sum_Hl_Ll,RealMatrix & sum_Hl_Llm1,RealMatrix & sum_Hlm1_Ll,RealMatrix & sum_Hlm1_Llm1,RealMatrix & sum_Hl_Hl,RealMatrix & sum_Hl_Hlm1,RealMatrix & sum_Hlm1_Hlm1,Real cost_ratio,size_t lev,const SizetArray & N_shared,RealMatrix & var_YHl,RealMatrix & rho_dot2_LH)2766 eval_ratio(RealMatrix& sum_Ll, RealMatrix& sum_Llm1, RealMatrix& sum_Hl,
2767 RealMatrix& sum_Hlm1, RealMatrix& sum_Ll_Ll, RealMatrix& sum_Ll_Llm1,
2768 RealMatrix& sum_Llm1_Llm1, RealMatrix& sum_Hl_Ll,
2769 RealMatrix& sum_Hl_Llm1, RealMatrix& sum_Hlm1_Ll,
2770 RealMatrix& sum_Hlm1_Llm1, RealMatrix& sum_Hl_Hl,
2771 RealMatrix& sum_Hl_Hlm1, RealMatrix& sum_Hlm1_Hlm1,
2772 Real cost_ratio, size_t lev, const SizetArray& N_shared,
2773 RealMatrix& var_YHl, RealMatrix& rho_dot2_LH)
2774 {
2775 if (lev == 0)
2776 return eval_ratio(sum_Ll, sum_Hl, sum_Ll_Ll, sum_Hl_Ll, sum_Hl_Hl,
2777 cost_ratio, lev, N_shared, var_YHl, rho_dot2_LH);
2778 else {
2779 Real beta_dot, gamma, eval_ratio, avg_eval_ratio = 0.;
2780 size_t qoi, num_avg = 0;
2781 for (qoi=0; qoi<numFunctions; ++qoi) {
2782 Real& rho_dot_sq = rho_dot2_LH(qoi,lev);
2783 compute_control(sum_Ll(qoi,lev), sum_Llm1(qoi,lev), sum_Hl(qoi,lev),
2784 sum_Hlm1(qoi,lev), sum_Ll_Ll(qoi,lev),
2785 sum_Ll_Llm1(qoi,lev), sum_Llm1_Llm1(qoi,lev),
2786 sum_Hl_Ll(qoi,lev), sum_Hl_Llm1(qoi,lev),
2787 sum_Hlm1_Ll(qoi,lev), sum_Hlm1_Llm1(qoi,lev),
2788 sum_Hl_Hl(qoi,lev), sum_Hl_Hlm1(qoi,lev),
2789 sum_Hlm1_Hlm1(qoi,lev), N_shared[qoi], var_YHl(qoi,lev),
2790 rho_dot_sq, beta_dot, gamma);
2791
2792 if (rho_dot_sq < 1.) { // protect against division by 0
2793 eval_ratio = std::sqrt(cost_ratio * rho_dot_sq / (1.-rho_dot_sq));
2794 if (outputLevel >= DEBUG_OUTPUT)
2795 Cout << "eval_ratio() QoI " << qoi+1 << ": cost_ratio = "
2796 << cost_ratio << " rho_dot_sq = " << rho_dot_sq
2797 << " eval_ratio = " << eval_ratio << std::endl;
2798 avg_eval_ratio += eval_ratio;
2799 ++num_avg;
2800 }
2801 }
2802 if (outputLevel >= DEBUG_OUTPUT) {
2803 Cout << "variance of HF Y[" << lev << "]:\n";
2804 write_col_vector_trans(Cout, (int)lev, (int)numFunctions, var_YHl);
2805 }
2806
2807 if (num_avg) avg_eval_ratio /= num_avg;
2808 else // should not happen, but provide a reasonable upper bound
2809 avg_eval_ratio = (Real)maxFunctionEvals / average(N_shared);
2810
2811 return avg_eval_ratio;
2812 }
2813 }
2814
2815
2816 Real NonDMultilevelSampling::
MSE_ratio(Real avg_eval_ratio,const RealVector & var_H,const RealVector & rho2_LH,size_t iter,const SizetArray & N_hf)2817 MSE_ratio(Real avg_eval_ratio, const RealVector& var_H,
2818 const RealVector& rho2_LH, size_t iter, const SizetArray& N_hf)
2819 {
2820 if (iter == 0) mcMSEIter0.sizeUninitialized(numFunctions);
2821
2822 Real mc_mse, cvmc_mse, mse_ratio, avg_mse_ratio = 0.;//,avg_mse_iter_ratio=0.;
2823 for (size_t qoi=0; qoi<numFunctions; ++qoi) {
2824 // Compute ratio of MSE for high fidelity MC and multifidelity CVMC
2825 mse_ratio = 1. - rho2_LH[qoi] * (1. - 1. / avg_eval_ratio); // Ng 2014
2826 mc_mse = var_H[qoi] / N_hf[qoi]; cvmc_mse = mc_mse * mse_ratio;
2827 Cout << "Mean square error for QoI " << qoi+1 << " reduced from " << mc_mse
2828 << " (MC) to " << cvmc_mse << " (CV); factor = " << mse_ratio << '\n';
2829
2830 if (iter == 0)
2831 { mcMSEIter0[qoi] = mc_mse; avg_mse_ratio += mse_ratio; }
2832 else
2833 avg_mse_ratio += cvmc_mse / mcMSEIter0[qoi];
2834 }
2835 //avg_mse_iter_ratio /= numFunctions;
2836 avg_mse_ratio /= numFunctions;
2837 Cout //<< "Average MSE reduction factor from CV for iteration "
2838 //<< std::setw(4) << iter << " = " << avg_mse_iter_ratio << '\n'
2839 << "Average MSE reduction factor since pilot MC = " << avg_mse_ratio
2840 << " targeting convergence tol = " << convergenceTol << "\n\n";
2841 return avg_mse_ratio;
2842 }
2843
2844
2845 void NonDMultilevelSampling::
cv_raw_moments(IntRealVectorMap & sum_L_shared,IntRealVectorMap & sum_H,IntRealVectorMap & sum_LL,IntRealVectorMap & sum_LH,const SizetArray & N_shared,IntRealVectorMap & sum_L_refined,const SizetArray & N_refined,const RealVector & rho2_LH,RealMatrix & H_raw_mom)2846 cv_raw_moments(IntRealVectorMap& sum_L_shared, IntRealVectorMap& sum_H,
2847 IntRealVectorMap& sum_LL, IntRealVectorMap& sum_LH,
2848 const SizetArray& N_shared, IntRealVectorMap& sum_L_refined,
2849 const SizetArray& N_refined, const RealVector& rho2_LH,
2850 RealMatrix& H_raw_mom)
2851 {
2852 if (H_raw_mom.empty()) H_raw_mom.shapeUninitialized(numFunctions, 4);
2853 RealVector beta(numFunctions, false);
2854
2855 // rho2_LH not stored for i > 1
2856 for (size_t qoi=0; qoi<numFunctions; ++qoi)
2857 Cout << "rho_LH (Pearson correlation) for QoI " << qoi+1 << " = "
2858 << std::setw(9) << std::sqrt(rho2_LH[qoi]) << '\n';
2859 //<< ", effectiveness ratio = " << std::setw(9) << rho2_LH[qoi] * cr1;
2860
2861 for (int i=1; i<=4; ++i) {
2862 compute_control(sum_L_shared[i], sum_H[i], sum_LL[i], sum_LH[i], N_shared,
2863 beta);
2864 Cout << "Moment " << i << ":\n";
2865 RealVector H_rm_col(Teuchos::View, H_raw_mom[i-1], numFunctions);
2866 apply_control(sum_H[i], sum_L_shared[i], N_shared, sum_L_refined[i],
2867 N_refined, beta, H_rm_col);
2868 }
2869 }
2870
2871
2872 void NonDMultilevelSampling::
cv_raw_moments(IntRealMatrixMap & sum_L_shared,IntRealMatrixMap & sum_H,IntRealMatrixMap & sum_LL,IntRealMatrixMap & sum_LH,const SizetArray & N_shared,IntRealMatrixMap & sum_L_refined,const SizetArray & N_refined,const RealMatrix & rho2_LH,size_t lev,RealMatrix & H_raw_mom)2873 cv_raw_moments(IntRealMatrixMap& sum_L_shared, IntRealMatrixMap& sum_H,
2874 IntRealMatrixMap& sum_LL, IntRealMatrixMap& sum_LH,
2875 const SizetArray& N_shared, IntRealMatrixMap& sum_L_refined,
2876 const SizetArray& N_refined, const RealMatrix& rho2_LH,
2877 size_t lev, RealMatrix& H_raw_mom)
2878 {
2879 if (H_raw_mom.empty()) H_raw_mom.shapeUninitialized(numFunctions, 4);
2880 RealVector beta(numFunctions, false);
2881
2882 // rho2_LH not stored for i > 1
2883 for (size_t qoi=0; qoi<numFunctions; ++qoi)
2884 Cout << "rho_LH (Pearson correlation) for QoI " << qoi+1 << " = "
2885 << std::setw(9) << std::sqrt(rho2_LH(qoi,lev)) << '\n';
2886 //<< ", effectiveness ratio = " << std::setw(9) << rho2_LH(qoi,lev)*cr1;
2887
2888 for (int i=1; i<=4; ++i) {
2889 compute_control(sum_L_shared[i], sum_H[i], sum_LL[i], sum_LH[i], N_shared,
2890 lev, beta);
2891 Cout << "Moment " << i << ":\n";
2892 RealVector H_rm_col(Teuchos::View, H_raw_mom[i-1], numFunctions);
2893 apply_control(sum_H[i], sum_L_shared[i], N_shared, sum_L_refined[i],
2894 N_refined, lev, beta, H_rm_col);
2895 }
2896 Cout << '\n'; // for loop over levels
2897 }
2898
2899
2900 void NonDMultilevelSampling::
cv_raw_moments(IntRealMatrixMap & sum_Ll,IntRealMatrixMap & sum_Llm1,IntRealMatrixMap & sum_Hl,IntRealMatrixMap & sum_Hlm1,IntRealMatrixMap & sum_Ll_Ll,IntRealMatrixMap & sum_Ll_Llm1,IntRealMatrixMap & sum_Llm1_Llm1,IntRealMatrixMap & sum_Hl_Ll,IntRealMatrixMap & sum_Hl_Llm1,IntRealMatrixMap & sum_Hlm1_Ll,IntRealMatrixMap & sum_Hlm1_Llm1,IntRealMatrixMap & sum_Hl_Hl,IntRealMatrixMap & sum_Hl_Hlm1,IntRealMatrixMap & sum_Hlm1_Hlm1,const SizetArray & N_shared,IntRealMatrixMap & sum_Ll_refined,IntRealMatrixMap & sum_Llm1_refined,const SizetArray & N_refined,const RealMatrix & rho_dot2_LH,size_t lev,RealMatrix & H_raw_mom)2901 cv_raw_moments(IntRealMatrixMap& sum_Ll, IntRealMatrixMap& sum_Llm1,
2902 IntRealMatrixMap& sum_Hl, IntRealMatrixMap& sum_Hlm1,
2903 IntRealMatrixMap& sum_Ll_Ll, IntRealMatrixMap& sum_Ll_Llm1,
2904 IntRealMatrixMap& sum_Llm1_Llm1, IntRealMatrixMap& sum_Hl_Ll,
2905 IntRealMatrixMap& sum_Hl_Llm1, IntRealMatrixMap& sum_Hlm1_Ll,
2906 IntRealMatrixMap& sum_Hlm1_Llm1, IntRealMatrixMap& sum_Hl_Hl,
2907 IntRealMatrixMap& sum_Hl_Hlm1, IntRealMatrixMap& sum_Hlm1_Hlm1,
2908 const SizetArray& N_shared, IntRealMatrixMap& sum_Ll_refined,
2909 IntRealMatrixMap& sum_Llm1_refined, const SizetArray& N_refined,
2910 const RealMatrix& rho_dot2_LH, size_t lev, RealMatrix& H_raw_mom)
2911 {
2912 if (lev == 0)
2913 cv_raw_moments(sum_Ll, sum_Hl, sum_Ll_Ll, sum_Hl_Ll, N_shared,
2914 sum_Ll_refined, N_refined, rho_dot2_LH, lev, H_raw_mom);
2915 else {
2916 if (H_raw_mom.empty()) H_raw_mom.shapeUninitialized(numFunctions, 4);
2917 RealVector beta_dot(numFunctions, false), gamma(numFunctions, false);
2918
2919 // rho_dot2_LH not stored for i > 1
2920 for (size_t qoi=0; qoi<numFunctions; ++qoi)
2921 Cout << "rho_dot_LH (Pearson correlation) for QoI " << qoi+1 << " = "
2922 << std::setw(9) << std::sqrt(rho_dot2_LH(qoi,lev)) << '\n';
2923 //<< ", effectiveness ratio = " << rho_dot2_LH(qoi,lev) * cr1;
2924
2925 // aggregate expected value of estimators for E[Y] for Y=LF^k or Y=HF^k
2926 for (int i=1; i<=4; ++i) {
2927 compute_control(sum_Ll[i], sum_Llm1[i], sum_Hl[i], sum_Hlm1[i],
2928 sum_Ll_Ll[i], sum_Ll_Llm1[i], sum_Llm1_Llm1[i],
2929 sum_Hl_Ll[i], sum_Hl_Llm1[i], sum_Hlm1_Ll[i],
2930 sum_Hlm1_Llm1[i], sum_Hl_Hl[i], sum_Hl_Hlm1[i],
2931 sum_Hlm1_Hlm1[i], N_shared, lev, beta_dot, gamma);
2932 Cout << "Moment " << i << ":\n";
2933 RealVector H_rm_col(Teuchos::View, H_raw_mom[i-1], numFunctions);
2934 apply_control(sum_Hl[i], sum_Hlm1[i], sum_Ll[i], sum_Llm1[i], N_shared,
2935 sum_Ll_refined[i], sum_Llm1_refined[i], N_refined, lev,
2936 beta_dot, gamma, H_rm_col);
2937 }
2938 Cout << '\n'; // for loop over levels
2939 }
2940 }
2941
2942
2943 void NonDMultilevelSampling::
compute_control(Real sum_L,Real sum_H,Real sum_LL,Real sum_LH,size_t N_shared,Real & beta)2944 compute_control(Real sum_L, Real sum_H, Real sum_LL, Real sum_LH,
2945 size_t N_shared, Real& beta)
2946 {
2947 // unbiased mean estimator X-bar = 1/N * sum
2948 // unbiased sample variance estimator = 1/(N-1) sum[(X_i - X-bar)^2]
2949 // = 1/(N-1) [ N Raw_X - N X-bar^2 ] = bessel * [Raw_X - X-bar^2]
2950 //Real mu_L = sum_L / N_shared, mu_H = sum_H / N_shared;
2951 //Real var_L = (sum_LL / N_shared - mu_L * mu_L) * bessel_corr,
2952 // cov_LH = (sum_LH / N_shared - mu_L * mu_H) * bessel_corr;
2953
2954 // Cancel repeated N_shared and bessel_corr within beta = cov_LH / var_L:
2955 beta = (sum_LH - sum_L * sum_H / N_shared)
2956 / (sum_LL - sum_L * sum_L / N_shared);
2957 }
2958
2959
2960 void NonDMultilevelSampling::
compute_control(Real sum_L,Real sum_H,Real sum_LL,Real sum_LH,Real sum_HH,size_t N_shared,Real & var_H,Real & rho2_LH)2961 compute_control(Real sum_L, Real sum_H, Real sum_LL, Real sum_LH, Real sum_HH,
2962 size_t N_shared, Real& var_H, Real& rho2_LH)
2963 {
2964 Real bessel_corr = (Real)N_shared / (Real)(N_shared - 1);
2965
2966 // unbiased mean estimator X-bar = 1/N * sum
2967 Real mu_L = sum_L / N_shared, mu_H = sum_H / N_shared;
2968 // unbiased sample variance estimator = 1/(N-1) sum[(X_i - X-bar)^2]
2969 // = 1/(N-1) [ N Raw_X - N X-bar^2 ] = bessel * [Raw_X - X-bar^2]
2970 Real var_L = (sum_LL / N_shared - mu_L * mu_L) * bessel_corr,
2971 cov_LH = (sum_LH / N_shared - mu_L * mu_H) * bessel_corr;
2972 var_H = (sum_HH / N_shared - mu_H * mu_H) * bessel_corr;
2973
2974 //beta = cov_LH / var_L;
2975 rho2_LH = cov_LH / var_L * cov_LH / var_H;
2976 }
2977
2978
2979 void NonDMultilevelSampling::
compute_control(Real sum_Ll,Real sum_Llm1,Real sum_Hl,Real sum_Hlm1,Real sum_Ll_Ll,Real sum_Ll_Llm1,Real sum_Llm1_Llm1,Real sum_Hl_Ll,Real sum_Hl_Llm1,Real sum_Hlm1_Ll,Real sum_Hlm1_Llm1,Real sum_Hl_Hl,Real sum_Hl_Hlm1,Real sum_Hlm1_Hlm1,size_t N_shared,Real & var_YH,Real & rho_dot2_LH,Real & beta_dot,Real & gamma)2980 compute_control(Real sum_Ll, Real sum_Llm1, Real sum_Hl, Real sum_Hlm1,
2981 Real sum_Ll_Ll, Real sum_Ll_Llm1, Real sum_Llm1_Llm1,
2982 Real sum_Hl_Ll, Real sum_Hl_Llm1, Real sum_Hlm1_Ll,
2983 Real sum_Hlm1_Llm1, Real sum_Hl_Hl, Real sum_Hl_Hlm1,
2984 Real sum_Hlm1_Hlm1, size_t N_shared, Real& var_YH,
2985 Real& rho_dot2_LH, Real& beta_dot, Real& gamma)
2986 {
2987 Real bessel_corr = (Real)N_shared / (Real)(N_shared - 1);
2988
2989 // means, variances, covariances for Q
2990 // Note: sum_*[i][lm1] is not the same as sum_*lm1[i][lev] due to
2991 // discrepancy evaluations with different sample sets!
2992 Real mu_Ll = sum_Ll / N_shared, mu_Llm1 = sum_Llm1 / N_shared;
2993 Real mu_Hl = sum_Hl / N_shared, mu_Hlm1 = sum_Hlm1 / N_shared;
2994
2995 Real var_Ll = (sum_Ll_Ll / N_shared - mu_Ll * mu_Ll) * bessel_corr;
2996 Real var_Llm1 = (sum_Llm1_Llm1 / N_shared - mu_Llm1 * mu_Llm1) * bessel_corr;
2997 Real var_Hl = (sum_Hl_Hl / N_shared - mu_Hl * mu_Hl) * bessel_corr;
2998 Real var_Hlm1 = (sum_Hlm1_Hlm1 / N_shared - mu_Hlm1 * mu_Hlm1) * bessel_corr;
2999
3000 Real cov_Hl_Ll = (sum_Hl_Ll / N_shared - mu_Hl * mu_Ll) * bessel_corr;
3001 Real cov_Hl_Llm1 = (sum_Hl_Llm1 / N_shared - mu_Hl * mu_Llm1) * bessel_corr;
3002 Real cov_Hlm1_Ll = (sum_Hlm1_Ll / N_shared - mu_Hlm1 * mu_Ll) * bessel_corr;
3003 Real cov_Hlm1_Llm1
3004 = (sum_Hlm1_Llm1 / N_shared - mu_Hlm1 * mu_Llm1) * bessel_corr;
3005
3006 Real cov_Ll_Llm1 = (sum_Ll_Llm1 / N_shared - mu_Ll * mu_Llm1) * bessel_corr;
3007 Real cov_Hl_Hlm1 = (sum_Hl_Hlm1 / N_shared - mu_Hl * mu_Hlm1) * bessel_corr;
3008
3009 // quantities derived from Q moments
3010 // gamma:
3011 Real cov_YHl_Ll = cov_Hl_Ll - cov_Hlm1_Ll;
3012 Real cov_YHl_Llm1 = cov_Hl_Llm1 - cov_Hlm1_Llm1;
3013 gamma = (cov_YHl_Llm1 * cov_Ll_Llm1 - var_Llm1 * cov_YHl_Ll)
3014 / (var_Ll * cov_YHl_Llm1 - cov_YHl_Ll * cov_Ll_Llm1);
3015 // theta, tau, beta:
3016 Real cov_YHl_YLldot = gamma * (cov_Hl_Ll - cov_Hlm1_Ll)
3017 - cov_Hl_Llm1 + cov_Hlm1_Llm1;
3018 Real cov_YHl_YLl = cov_Hl_Ll - cov_Hlm1_Ll - cov_Hl_Llm1 + cov_Hlm1_Llm1;
3019 Real var_YLldot = gamma * (gamma * var_Ll - 2. * cov_Ll_Llm1) + var_Llm1;
3020 Real var_YHl = var_Hl - 2. * cov_Hl_Hlm1 + var_Hlm1,
3021 var_YLl = var_Ll - 2. * cov_Ll_Llm1 + var_Llm1,
3022 theta = cov_YHl_YLldot / cov_YHl_YLl, tau = var_YLldot / var_YLl;
3023 // carry forwards:
3024 var_YH = var_Hl - 2. * cov_Hl_Hlm1 + var_Hlm1; // var(H_l - H_lm1)
3025 beta_dot = cov_YHl_YLldot / var_YLldot;
3026
3027 // compute evaluation ratio which determines increment for LF samples
3028 // > the sample increment optimizes the total computational budget and is
3029 // not treated as a worst case accuracy reqmt --> use the QoI average
3030 // > refinement based only on QoI mean statistics
3031 // Given use of 1/r in MSE_ratio, one approach would average 1/r, but
3032 // this does not seem to behave as well in limited numerical experience.
3033 Real rho2_LH = cov_YHl_YLl / var_YHl * cov_YHl_YLl / var_YLl,
3034 ratio = theta * theta / tau;
3035 // variance reduction test
3036 //if (ratio < 1.) ...switch to Ycorr-based control...
3037 rho_dot2_LH = rho2_LH * ratio;
3038 //rho_dot2_LH = cov_YHl_YLldot / var_YHl * cov_YHl_YLldot / var_YLldot;
3039
3040 if (outputLevel == DEBUG_OUTPUT)
3041 Cout << "compute_control(): var reduce ratio = " << ratio << " rho2_LH = "
3042 << rho2_LH << " rho_dot2_LH = " << rho_dot2_LH << std::endl;
3043 }
3044
3045
3046 void NonDMultilevelSampling::
apply_control(Real sum_H,Real sum_L_shared,size_t N_shared,Real sum_L_refined,size_t N_refined,Real beta,Real & H_raw_mom)3047 apply_control(Real sum_H, Real sum_L_shared, size_t N_shared,
3048 Real sum_L_refined, size_t N_refined, Real beta, Real& H_raw_mom)
3049 {
3050 // apply control for HF uncentered raw moment estimates:
3051 H_raw_mom = sum_H / N_shared // mu_H from shared samples
3052 - beta * (sum_L_shared / N_shared - // mu_L from shared samples
3053 sum_L_refined / N_refined); // refined_mu_L incl increment
3054 }
3055
3056
3057 void NonDMultilevelSampling::
apply_control(Real sum_Hl,Real sum_Hlm1,Real sum_Ll,Real sum_Llm1,size_t N_shared,Real sum_Ll_refined,Real sum_Llm1_refined,size_t N_refined,Real beta_dot,Real gamma,Real & H_raw_mom)3058 apply_control(Real sum_Hl, Real sum_Hlm1, Real sum_Ll, Real sum_Llm1,
3059 size_t N_shared, Real sum_Ll_refined, Real sum_Llm1_refined,
3060 size_t N_refined, Real beta_dot, Real gamma, Real& H_raw_mom)
3061 {
3062 // updated LF expectations following final sample increment:
3063 Real mu_Hl = sum_Hl / N_shared, mu_Hlm1 = sum_Hlm1 / N_shared,
3064 mu_Ll = sum_Ll / N_shared, mu_Llm1 = sum_Llm1 / N_shared;
3065 Real refined_mu_Ll = sum_Ll_refined / N_refined;
3066 Real refined_mu_Llm1 = sum_Llm1_refined / N_refined;
3067
3068 // apply control for HF uncentered raw moment estimates:
3069 Real mu_YH = mu_Hl - mu_Hlm1;
3070 Real mu_YLdot = gamma * mu_Ll - mu_Llm1;
3071 Real refined_mu_YLdot = gamma * refined_mu_Ll - refined_mu_Llm1;
3072 H_raw_mom = mu_YH - beta_dot * (mu_YLdot - refined_mu_YLdot);
3073 }
3074
3075
3076 void NonDMultilevelSampling::
convert_moments(const RealMatrix & raw_mom,RealMatrix & final_mom)3077 convert_moments(const RealMatrix& raw_mom, RealMatrix& final_mom)
3078 {
3079 // Note: raw_mom is numFunctions x 4 and final_mom is the transpose
3080 if (final_mom.empty())
3081 final_mom.shapeUninitialized(4, numFunctions);
3082
3083 // Convert uncentered raw moment estimates to central moments
3084 if (finalMomentsType == CENTRAL_MOMENTS) {
3085 for (size_t qoi=0; qoi<numFunctions; ++qoi)
3086 uncentered_to_centered(raw_mom(qoi,0), raw_mom(qoi,1), raw_mom(qoi,2),
3087 raw_mom(qoi,3), final_mom(0,qoi), final_mom(1,qoi),
3088 final_mom(2,qoi), final_mom(3,qoi));
3089 }
3090 // Convert uncentered raw moment estimates to standardized moments
3091 else { //if (finalMomentsType == STANDARD_MOMENTS) {
3092 Real cm1, cm2, cm3, cm4;
3093 for (size_t qoi=0; qoi<numFunctions; ++qoi) {
3094 uncentered_to_centered(raw_mom(qoi,0), raw_mom(qoi,1), raw_mom(qoi,2),
3095 raw_mom(qoi,3), cm1, cm2, cm3, cm4);
3096 centered_to_standard(cm1, cm2, cm3, cm4, final_mom(0,qoi),
3097 final_mom(1,qoi), final_mom(2,qoi),
3098 final_mom(3,qoi));
3099 }
3100 }
3101
3102 if (outputLevel >= DEBUG_OUTPUT)
3103 for (size_t qoi=0; qoi<numFunctions; ++qoi)
3104 Cout << "raw mom 1 = " << raw_mom(qoi,0)
3105 << " final mom 1 = " << final_mom(0,qoi) << '\n'
3106 << "raw mom 2 = " << raw_mom(qoi,1)
3107 << " final mom 2 = " << final_mom(1,qoi) << '\n'
3108 << "raw mom 3 = " << raw_mom(qoi,2)
3109 << " final mom 3 = " << final_mom(2,qoi) << '\n'
3110 << "raw mom 4 = " << raw_mom(qoi,3)
3111 << " final mom 4 = " << final_mom(3,qoi) << "\n\n";
3112 }
3113
3114
3115 void NonDMultilevelSampling::
compute_error_estimates(IntRealMatrixMap & sum_Ql,IntRealMatrixMap & sum_Qlm1,IntIntPairRealMatrixMap & sum_QlQlm1,Sizet2DArray & num_Q)3116 compute_error_estimates(IntRealMatrixMap &sum_Ql, IntRealMatrixMap &sum_Qlm1,
3117 IntIntPairRealMatrixMap &sum_QlQlm1,
3118 Sizet2DArray &num_Q) {
3119 if (!finalMomentsType)
3120 return;
3121
3122 if (finalStatErrors.empty())
3123 finalStatErrors.size(finalStatistics.num_functions()); // init to 0.
3124
3125 Real agg_estim_var, var_Yl, cm1l, cm2l, cm3l, cm4l, cm1lm1, cm2lm1,
3126 cm3lm1, cm4lm1, cm1l_sq, cm1lm1_sq, cm2l_sq, cm2lm1_sq, var_Ql, var_Qlm1,
3127 mu_Q2l, mu_Q2lm1, mu_Q2lQ2lm1,
3128 mu_Q1lm1_mu_Q2lQ1lm1, mu_Q1lm1_mu_Q1lm1_muQ2l, mu_Q1l_mu_Q1lQ2lm1, mu_Q1l_mu_Q1l_mu_Q2lm1,
3129 mu_Q1l_mu_Qlm1_mu_Q1lQ1lm1, mu_Q1l_mu_Q1l_mu_Q1lm1_muQ1lm1, mu_Q2l_muQ2lm1, mu_Q1lQ1lm1_mu_Q1lQ1lm1,
3130 mu_P2lP2lm1, var_P2l, var_P2lm1, covar_P2lP2lm1, term, bessel_corr;
3131 size_t lev, qoi, cntr = 0, Nlq,
3132 num_lev = iteratedModel.truth_model().solution_levels();
3133 IntIntPair pr11(1, 1), pr12(1, 2), pr21(2, 1), pr22(2, 2);
3134 RealMatrix &sum_Q1l = sum_Ql[1], &sum_Q1lm1 = sum_Qlm1[1],
3135 &sum_Q2l = sum_Ql[2], &sum_Q2lm1 = sum_Qlm1[2],
3136 &sum_Q3l = sum_Ql[3], &sum_Q3lm1 = sum_Qlm1[3],
3137 &sum_Q4l = sum_Ql[4], &sum_Q4lm1 = sum_Qlm1[4],
3138 &sum_Q1lQ1lm1 = sum_QlQlm1[pr11], &sum_Q1lQ2lm1 = sum_QlQlm1[pr12],
3139 &sum_Q2lQ1lm1 = sum_QlQlm1[pr21], &sum_Q2lQ2lm1 = sum_QlQlm1[pr22];
3140 for (qoi = 0; qoi < numFunctions; ++qoi) {
3141
3142 // std error in mean estimate
3143 lev = 0;
3144 Nlq = num_Q[lev][qoi];
3145
3146 bessel_corr = (Real) Nlq / (Real) (Nlq - 1);
3147 cm1l = sum_Q1l(qoi, lev) / Nlq;
3148 var_Yl = (sum_Q2l(qoi, lev) / Nlq - cm1l * cm1l) * bessel_corr; // var_Ql
3149 agg_estim_var = var_Yl / Nlq;
3150 for (lev = 1; lev < num_lev; ++lev) {
3151 Nlq = num_Q[lev][qoi];
3152 cm1l = sum_Q1l(qoi, lev) / Nlq;
3153 cm1lm1 = sum_Q1lm1(qoi, lev) / Nlq;
3154 //var_Yl = var_Ql - 2.* covar_QlQlm1 + var_Qlm1;
3155 var_Yl = (sum_Q2l(qoi, lev) / Nlq - cm1l * cm1l
3156 - 2. * (sum_Q1lQ1lm1(qoi, lev) / Nlq - cm1l * cm1lm1) +
3157 sum_Q2lm1(qoi, lev) / Nlq - cm1lm1 * cm1lm1) * bessel_corr;
3158 agg_estim_var += var_Yl / Nlq;
3159 }
3160 check_negative(agg_estim_var);
3161
3162 finalStatErrors[cntr++] = std::sqrt(agg_estim_var); // std error
3163 if (outputLevel >= DEBUG_OUTPUT) {
3164 Cout << "Estimator variance for mean = " << agg_estim_var << "\n";
3165 Cout << "Estimator SE for mean = " << finalStatErrors[cntr - 1] << "\n";
3166 }
3167 // std error in variance or std deviation estimate
3168 lev = 0;
3169 Nlq = num_Q[lev][qoi];
3170 uncentered_to_centered(sum_Q1l(qoi, lev) / Nlq, sum_Q2l(qoi, lev) / Nlq,
3171 sum_Q3l(qoi, lev) / Nlq, sum_Q4l(qoi, lev) / Nlq,
3172 cm1l, cm2l, cm3l, cm4l, Nlq);
3173 cm2l_sq = cm2l * cm2l;
3174
3175 //[fm] bias correction for var_P2l
3176 var_P2l = Nlq * (Nlq - 1.) / (Nlq * Nlq - 2. * Nlq + 3.) * (cm4l - (Nlq - 3.) / (Nlq - 1.) * cm2l_sq);
3177 agg_estim_var = var_P2l / Nlq;
3178 for (lev = 1; lev < num_lev; ++lev) {
3179 Nlq = num_Q[lev][qoi];
3180 mu_Q2l = sum_Q2l(qoi, lev) / Nlq;
3181 uncentered_to_centered(sum_Q1l(qoi, lev) / Nlq, mu_Q2l,
3182 sum_Q3l(qoi, lev) / Nlq, sum_Q4l(qoi, lev) / Nlq,
3183 cm1l, cm2l, cm3l, cm4l, Nlq);
3184 mu_Q2lm1 = sum_Q2lm1(qoi, lev) / Nlq;
3185 uncentered_to_centered(sum_Q1lm1(qoi, lev) / Nlq, mu_Q2lm1,
3186 sum_Q3lm1(qoi, lev) / Nlq, sum_Q4lm1(qoi, lev) / Nlq,
3187 cm1lm1, cm2lm1, cm3lm1, cm4lm1, Nlq);
3188 cm2l_sq = cm2l * cm2l;
3189 cm2lm1_sq = cm2lm1 * cm2lm1;
3190
3191 //[fm] unbiased products of mean
3192 mu_Q2lQ2lm1 = sum_Q2lQ2lm1(qoi, lev) / Nlq;
3193 mu_Q1lm1_mu_Q2lQ1lm1 = unbiased_mean_product_pair(sum_Q1lm1(qoi, lev), sum_Q2lQ1lm1(qoi, lev),
3194 sum_Q2lQ2lm1(qoi, lev), Nlq);
3195 mu_Q1lm1_mu_Q1lm1_muQ2l = unbiased_mean_product_triplet(sum_Q1lm1(qoi, lev), sum_Q1lm1(qoi, lev),
3196 sum_Q2l(qoi, lev),
3197 sum_Q2lm1(qoi, lev), sum_Q2lQ1lm1(qoi, lev),
3198 sum_Q2lQ1lm1(qoi, lev),
3199 sum_Q2lQ2lm1(qoi, lev), Nlq);
3200 mu_Q1l_mu_Q1lQ2lm1 = unbiased_mean_product_pair(sum_Q1l(qoi, lev), sum_Q1lQ2lm1(qoi, lev),
3201 sum_Q2lQ2lm1(qoi, lev), Nlq);
3202 mu_Q1l_mu_Q1l_mu_Q2lm1 = unbiased_mean_product_triplet(sum_Q1l(qoi, lev), sum_Q1l(qoi, lev),
3203 sum_Q2lm1(qoi, lev),
3204 sum_Q2l(qoi, lev), sum_Q1lQ2lm1(qoi, lev),
3205 sum_Q1lQ2lm1(qoi, lev),
3206 sum_Q2lQ2lm1(qoi, lev), Nlq);
3207 mu_Q1l_mu_Qlm1_mu_Q1lQ1lm1 = unbiased_mean_product_triplet(sum_Q1l(qoi, lev), sum_Q1lm1(qoi, lev),
3208 sum_Q1lQ1lm1(qoi, lev),
3209 sum_Q1lQ1lm1(qoi, lev), sum_Q2lQ1lm1(qoi, lev),
3210 sum_Q1lQ2lm1(qoi, lev),
3211 sum_Q2lQ2lm1(qoi, lev), Nlq);
3212 mu_Q1l_mu_Q1l_mu_Q1lm1_muQ1lm1 = unbiased_mean_product_pairpair(sum_Q1l(qoi, lev), sum_Q1lm1(qoi, lev),
3213 sum_Q1lQ1lm1(qoi, lev),
3214 sum_Q2l(qoi, lev), sum_Q2lm1(qoi, lev),
3215 sum_Q2lQ1lm1(qoi, lev), sum_Q1lQ2lm1(qoi, lev),
3216 sum_Q2lQ2lm1(qoi, lev), Nlq);
3217 mu_Q2l_muQ2lm1 = unbiased_mean_product_pair(sum_Q2l(qoi, lev), sum_Q2lm1(qoi, lev), sum_Q2lQ2lm1(qoi, lev),
3218 Nlq);
3219 mu_P2lP2lm1 = mu_Q2lQ2lm1 //E[QL2 Ql2]
3220 - 2. * mu_Q1lm1_mu_Q2lQ1lm1 //E[Ql] E[QL2Ql]
3221 + 2. * mu_Q1lm1_mu_Q1lm1_muQ2l //E[Ql]2 E[QL2]
3222 - 2. * mu_Q1l_mu_Q1lQ2lm1 //E[QL] E[QLQl2]
3223 + 2. * mu_Q1l_mu_Q1l_mu_Q2lm1 //E[QL]2 E[Ql2]
3224 + 4. * mu_Q1l_mu_Qlm1_mu_Q1lQ1lm1 //E[QL] E[Ql] E[QLQl]
3225 - 4. * mu_Q1l_mu_Q1l_mu_Q1lm1_muQ1lm1 //E[QL]2 E[Ql]2
3226 - mu_Q2l_muQ2lm1; //E[QL2] E[Ql2]
3227
3228 // [fm] bias correction for var_P2l and var_P2lm1
3229 var_P2l = Nlq * (Nlq - 1.) / (Nlq * Nlq - 2. * Nlq + 3.) * (cm4l - (Nlq - 3.) / (Nlq - 1.) * cm2l_sq);
3230 var_P2lm1 =
3231 Nlq * (Nlq - 1.) / (Nlq * Nlq - 2. * Nlq + 3.) * (cm4lm1 - (Nlq - 3.) / (Nlq - 1.) * cm2lm1_sq);
3232
3233 // [fm] unbiased by opening up the square and compute three different term
3234 mu_Q1lQ1lm1_mu_Q1lQ1lm1 = unbiased_mean_product_pair(sum_Q1lQ1lm1(qoi, lev), sum_Q1lQ1lm1(qoi, lev),
3235 sum_Q2lQ2lm1(qoi, lev), Nlq);
3236 term = mu_Q1lQ1lm1_mu_Q1lQ1lm1 - 2. * mu_Q1l_mu_Qlm1_mu_Q1lQ1lm1 + mu_Q1l_mu_Q1l_mu_Q1lm1_muQ1lm1;
3237
3238 //[fm] Using an unbiased estimator we include the var_Ql * var_Qlm1 term in mu_P2lP2lm1
3239 // and term is already squared out
3240 covar_P2lP2lm1
3241 = mu_P2lP2lm1 + term / (Nlq - 1.);
3242 agg_estim_var += (var_P2l + var_P2lm1 - 2. * covar_P2lP2lm1) / Nlq;
3243
3244 if (outputLevel >= DEBUG_OUTPUT) {
3245 Cout << "Estimator for covariance for variance = " << covar_P2lP2lm1 << "\n";
3246 Cout << "Estimator for covariance for variance first term = " << mu_P2lP2lm1 << "\n";
3247 Cout << "Estimator for covariance for variance second term = " << term << "\n";
3248 }
3249 }
3250 check_negative(agg_estim_var);
3251 if (outputLevel >= DEBUG_OUTPUT)
3252 Cout << "Estimator SE for variance = " << sqrt(agg_estim_var) << "\n";
3253
3254 Real mom2 = momentStats(1, qoi);
3255 if (finalMomentsType == STANDARD_MOMENTS && mom2 > 0.) {
3256 // std error of std deviation estimator
3257
3258 // An approximation for std error of a fn of another std error estimator
3259 // = derivative of function * std error of the other estimator -->
3260 // d/dtheta of sqrt( variance(theta) ) = 1/2 variance^{-1/2} = 1/(2 stdev)
3261 // Note: this approx. assumes normality in the estimator distribution.
3262 // Harding et al. 2014 assumes normality in the QoI distribution and has
3263 // been observed to contain bias in numerical experiments, whereas bias
3264 // in the derivative approx goes to zero asymptotically.
3265 finalStatErrors[cntr] = std::sqrt(agg_estim_var) / (2. * mom2);
3266 ++cntr;
3267 if (outputLevel >= DEBUG_OUTPUT)
3268 Cout << "Estimator SE for stddev = " << finalStatErrors[cntr - 1] << "\n\n";
3269 } else // std error of variance estimator
3270 finalStatErrors[cntr++] = std::sqrt(agg_estim_var);
3271
3272 // level mapping errors not implemented at this time
3273 cntr +=
3274 requestedRespLevels[qoi].length() + requestedProbLevels[qoi].length() +
3275 requestedRelLevels[qoi].length() + requestedGenRelLevels[qoi].length();
3276 }
3277 }
3278
3279 void NonDMultilevelSampling::
export_all_samples(String root_prepend,const Model & model,size_t iter,size_t lev)3280 export_all_samples(String root_prepend, const Model& model, size_t iter,
3281 size_t lev)
3282 {
3283 String tabular_filename(root_prepend);
3284 const String& iface_id = model.interface_id();
3285 size_t i, num_samp = allSamples.numCols();
3286 if (iface_id.empty()) tabular_filename += "NO_ID_i";
3287 else tabular_filename += iface_id + "_i";
3288 tabular_filename += std::to_string(iter) + "_l"
3289 + std::to_string(lev) + '_'
3290 + std::to_string(num_samp) + ".dat";
3291 Variables vars(model.current_variables().copy());
3292
3293 String context_message("NonDMultilevelSampling::export_all_samples");
3294 StringArray no_resp_labels; String cntr_label("sample_id");
3295
3296 // Rather than hard override, rely on output_precision user spec
3297 //int save_wp = write_precision;
3298 //write_precision = 16; // override
3299 std::ofstream tabular_stream;
3300 TabularIO::open_file(tabular_stream, tabular_filename, context_message);
3301 TabularIO::write_header_tabular(tabular_stream, vars, no_resp_labels,
3302 cntr_label, exportSamplesFormat);
3303 for (i=0; i<num_samp; ++i) {
3304 sample_to_variables(allSamples[i], vars); // NonDSampling version
3305 TabularIO::write_data_tabular(tabular_stream, vars, iface_id, i+1,
3306 exportSamplesFormat);
3307 }
3308
3309 TabularIO::close_file(tabular_stream, tabular_filename, context_message);
3310 //write_precision = save_wp; // restore
3311 }
3312
3313
post_run(std::ostream & s)3314 void NonDMultilevelSampling::post_run(std::ostream& s)
3315 {
3316 // Final moments are generated within core_run() by convert_moments().
3317 // No addtional stats are currently supported.
3318 //if (statsFlag) // calculate statistics on allResponses
3319 // compute_statistics(allSamples, allResponses);
3320
3321 // NonD::update_aleatory_final_statistics() pushes momentStats into
3322 // finalStatistics
3323 update_final_statistics();
3324
3325 Analyzer::post_run(s);
3326 }
3327
3328
print_results(std::ostream & s,short results_state)3329 void NonDMultilevelSampling::print_results(std::ostream& s, short results_state)
3330 {
3331 if (statsFlag) {
3332 print_multilevel_evaluation_summary(s, NLev);
3333 s << "<<<<< Equivalent number of high fidelity evaluations: "
3334 << equivHFEvals << "\n\nStatistics based on multilevel sample set:\n";
3335
3336 //print_statistics(s);
3337 print_moments(s, "response function",
3338 iteratedModel.truth_model().response_labels());
3339 archive_moments();
3340 archive_equiv_hf_evals(equivHFEvals);
3341 }
3342 }
3343
3344 } // namespace Dakota
3345