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: NonDSampling
10 //- Description: Wrapper class for Fortran 90 LHS library
11 //- Owner: Mike Eldred
12 //- Checked by:
13 //- Version:
14
15 #ifndef NOND_MULTILEVEL_SAMPLING_H
16 #define NOND_MULTILEVEL_SAMPLING_H
17
18 #include "NonDSampling.hpp"
19 #include "DataMethod.hpp"
20
21
22 namespace Dakota {
23
24 /// Performs Multilevel Monte Carlo sampling for uncertainty quantification.
25
26 /** Multilevel Monte Carlo (MLMC) is a variance-reduction technique
27 that utilitizes lower fidelity simulations that have response QoI
28 that are correlated with the high-fidelity response QoI. */
29
30 class NonDMultilevelSampling: public NonDSampling
31 {
32 public:
33
34 //
35 //- Heading: Constructors and destructor
36 //
37
38 /// standard constructor
39 NonDMultilevelSampling(ProblemDescDB& problem_db, Model& model);
40 /// destructor
41 ~NonDMultilevelSampling();
42
43 //
44 //- Heading: Virtual function redefinitions
45 //
46
47 bool resize();
48
49 protected:
50
51 //
52 //- Heading: Virtual function redefinitions
53 //
54
55 void pre_run();
56 void core_run();
57 void post_run(std::ostream& s);
58 void print_results(std::ostream& s, short results_state = FINAL_RESULTS);
59
60 private:
61
62 //
63 //- Heading: Helper functions
64 //
65
66 /// Perform multilevel Monte Carlo across the discretization levels for a
67 /// particular model form using discrepancy accumulators (sum_Y)
68 void multilevel_mc_Ysum(unsigned short model_form);
69 /// Perform multilevel Monte Carlo across the discretization levels for a
70 /// particular model form using QoI accumulators (sum_Q)
71 void multilevel_mc_Qsum(unsigned short model_form);
72 /// Perform control variate Monte Carlo across two model forms
73 void control_variate_mc(const UShortArray& active_key);
74 /// Perform multilevel Monte Carlo across levels in combination with
75 /// control variate Monte Carlo across model forms at each level; CV
76 /// computes correlations for Y (LH correlations for level discrepancies)
77 void multilevel_control_variate_mc_Ycorr(unsigned short lf_model_form,
78 unsigned short hf_model_form);
79 /// Perform multilevel Monte Carlo across levels in combination with
80 /// control variate Monte Carlo across model forms at each level; CV
81 /// computes correlations for Q (LH correlations for QoI)
82 void multilevel_control_variate_mc_Qcorr(unsigned short lf_model_form,
83 unsigned short hf_model_form);
84
85 /// perform a shared increment of LF and HF samples for purposes of
86 /// computing/updating the evaluation ratio and the MSE ratio
87 void shared_increment(size_t iter, size_t lev);
88 /// perform final LF sample increment as indicated by the evaluation ratio
89 bool lf_increment(Real avg_eval_ratio, const SizetArray& N_lf,
90 const SizetArray& N_hf, size_t iter, size_t lev);
91
92 /// synchronize iteratedModel and activeSet on AGGREGATED_MODELS mode
93 void aggregated_models_mode();
94 /// synchronize iteratedModel and activeSet on BYPASS_SURROGATE mode
95 void bypass_surrogate_mode();
96 /// synchronize iteratedModel and activeSet on UNCORRECTED_SURROGATE mode
97 void uncorrected_surrogate_mode();
98
99 /// return (aggregate) level cost
100 Real level_cost(const RealVector& cost, unsigned short step);
101
102 /// advance any sequence specifications
103 void assign_specification_sequence(size_t index);
104 /// extract current random seed from randomSeedSeqSpec
105 int random_seed(size_t index) const;
106
107 /// manage response mode and active model key from {group,form,lev} triplet.
108 /// s_index is the sequence index that defines the active dimension for a
109 /// model sequence.
110 void configure_indices(unsigned short group, unsigned short form,
111 unsigned short lev, unsigned short s_index);
112
113 /// initialize the ML accumulators for computing means, variances, and
114 /// covariances across fidelity levels
115 void initialize_ml_Ysums(IntRealMatrixMap& sum_Y, size_t num_lev);
116 /// initialize the ML accumulators for computing means, variances, and
117 /// covariances across fidelity levels
118 void initialize_ml_Qsums(IntRealMatrixMap& sum_Ql, IntRealMatrixMap& sum_Qlm1,
119 IntIntPairRealMatrixMap& sum_QlQlm1, size_t num_lev);
120
121 /// initialize the CV accumulators for computing means, variances, and
122 /// covariances across fidelity levels
123 void initialize_cv_sums(IntRealVectorMap& sum_L_shared,
124 IntRealVectorMap& sum_L_refined,
125 IntRealVectorMap& sum_H, IntRealVectorMap& sum_LL,
126 IntRealVectorMap& sum_LH);
127
128 /// initialize the MLCV accumulators for computing means, variances, and
129 /// covariances across fidelity levels
130 void initialize_mlcv_sums(IntRealMatrixMap& sum_L_shared,
131 IntRealMatrixMap& sum_L_refined,
132 IntRealMatrixMap& sum_H,
133 IntRealMatrixMap& sum_LL, IntRealMatrixMap& sum_LH,
134 IntRealMatrixMap& sum_HH, size_t num_ml_lev,
135 size_t num_cv_lev);
136 /// initialize the MLCV accumulators for computing means, variances, and
137 /// covariances across fidelity levels
138 void initialize_mlcv_sums(IntRealMatrixMap& sum_Ll,
139 IntRealMatrixMap& sum_Llm1,
140 IntRealMatrixMap& sum_Ll_refined,
141 IntRealMatrixMap& sum_Llm1_refined,
142 IntRealMatrixMap& sum_Hl,
143 IntRealMatrixMap& sum_Hlm1,
144 IntRealMatrixMap& sum_Ll_Ll,
145 IntRealMatrixMap& sum_Ll_Llm1,
146 IntRealMatrixMap& sum_Llm1_Llm1,
147 IntRealMatrixMap& sum_Hl_Ll,
148 IntRealMatrixMap& sum_Hl_Llm1,
149 IntRealMatrixMap& sum_Hlm1_Ll,
150 IntRealMatrixMap& sum_Hlm1_Llm1,
151 IntRealMatrixMap& sum_Hl_Hl,
152 IntRealMatrixMap& sum_Hl_Hlm1,
153 IntRealMatrixMap& sum_Hlm1_Hlm1,
154 size_t num_ml_lev, size_t num_cv_lev);
155
156 /// accumulate initial approximation to mean vector, for use as offsets in
157 /// subsequent accumulations
158 void accumulate_offsets(RealVector& mu);
159
160 /// update running QoI sums for one model (sum_Q) using set of model
161 /// evaluations within allResponses; used for level 0 from other accumulators
162 void accumulate_ml_Qsums(IntRealMatrixMap& sum_Q, size_t lev,
163 const RealVector& offset, SizetArray& num_Q);
164 /// update accumulators for multilevel telescoping running sums
165 /// using set of model evaluations within allResponses
166 void accumulate_ml_Ysums(IntRealMatrixMap& sum_Y, RealMatrix& sum_YY,
167 size_t lev, const RealVector& offset,
168 SizetArray& num_Y);
169 /// update running QoI sums for two models (sum_Ql, sum_Qlm1) using set of
170 /// model evaluations within allResponses
171 void accumulate_ml_Qsums(IntRealMatrixMap& sum_Ql, IntRealMatrixMap& sum_Qlm1,
172 IntIntPairRealMatrixMap& sum_QlQlm1, size_t lev,
173 const RealVector& offset, SizetArray& num_Q);
174
175 /// update running sums for one model (sum_L) using set of model
176 /// evaluations within allResponses
177 void accumulate_cv_sums(IntRealVectorMap& sum_L, const RealVector& offset,
178 SizetArray& num_L);
179 /// update running sums for two models (sum_L, sum_H, and sum_LH)
180 /// from set of low/high fidelity model evaluations within allResponses
181 void accumulate_cv_sums(IntRealVectorMap& sum_L_shared,
182 IntRealVectorMap& sum_L_refined,
183 IntRealVectorMap& sum_H, IntRealVectorMap& sum_LL,
184 IntRealVectorMap& sum_LH, RealVector& sum_HH,
185 const RealVector& offset, SizetArray& num_L,
186 SizetArray& num_H);
187
188 /// update running QoI sums for one model at two levels (sum_Ql, sum_Qlm1)
189 /// using set of model evaluations within allResponses
190 void accumulate_mlcv_Qsums(IntRealMatrixMap& sum_Ql,
191 IntRealMatrixMap& sum_Qlm1, size_t lev,
192 const RealVector& offset, SizetArray& num_Q);
193 /// update running discrepancy sums for one model (sum_Y) using
194 /// set of model evaluations within allResponses
195 void accumulate_mlcv_Ysums(IntRealMatrixMap& sum_Y, size_t lev,
196 const RealVector& offset, SizetArray& num_Y);
197 /// update running QoI sums for two models (sum_L, sum_H, sum_LL, sum_LH,
198 /// and sum_HH) from set of low/high fidelity model evaluations within
199 /// {lf,hf}_resp_map; used for level 0 from other accumulators
200 void accumulate_mlcv_Qsums(const IntResponseMap& lf_resp_map,
201 const IntResponseMap& hf_resp_map,
202 IntRealMatrixMap& sum_L_shared,
203 IntRealMatrixMap& sum_L_refined,
204 IntRealMatrixMap& sum_H, IntRealMatrixMap& sum_LL,
205 IntRealMatrixMap& sum_LH, IntRealMatrixMap& sum_HH,
206 size_t lev, const RealVector& lf_offset,
207 const RealVector& hf_offset, SizetArray& num_L,
208 SizetArray& num_H);
209 /// update running two-level discrepancy sums for two models (sum_L,
210 /// sum_H, sum_LL, sum_LH, and sum_HH) from set of low/high fidelity
211 /// model evaluations within {lf,hf}resp_map
212 void accumulate_mlcv_Ysums(const IntResponseMap& lf_resp_map,
213 const IntResponseMap& hf_resp_map,
214 IntRealMatrixMap& sum_L_shared,
215 IntRealMatrixMap& sum_L_refined,
216 IntRealMatrixMap& sum_H, IntRealMatrixMap& sum_LL,
217 IntRealMatrixMap& sum_LH, IntRealMatrixMap& sum_HH,
218 size_t lev, const RealVector& lf_offset,
219 const RealVector& hf_offset,
220 SizetArray& num_L, SizetArray& num_H);
221 /// update running QoI sums for two models and two levels from set
222 /// of low/high fidelity model evaluations within {lf,hf}_resp_map
223 void accumulate_mlcv_Qsums(const IntResponseMap& lf_resp_map,
224 const IntResponseMap& hf_resp_map,
225 IntRealMatrixMap& sum_Ll,
226 IntRealMatrixMap& sum_Llm1,
227 IntRealMatrixMap& sum_Ll_refined,
228 IntRealMatrixMap& sum_Llm1_refined,
229 IntRealMatrixMap& sum_Hl,
230 IntRealMatrixMap& sum_Hlm1,
231 IntRealMatrixMap& sum_Ll_Ll,
232 IntRealMatrixMap& sum_Ll_Llm1,
233 IntRealMatrixMap& sum_Llm1_Llm1,
234 IntRealMatrixMap& sum_Hl_Ll,
235 IntRealMatrixMap& sum_Hl_Llm1,
236 IntRealMatrixMap& sum_Hlm1_Ll,
237 IntRealMatrixMap& sum_Hlm1_Llm1,
238 IntRealMatrixMap& sum_Hl_Hl,
239 IntRealMatrixMap& sum_Hl_Hlm1,
240 IntRealMatrixMap& sum_Hlm1_Hlm1, size_t lev,
241 const RealVector& lf_offset,
242 const RealVector& hf_offset,
243 SizetArray& num_L, SizetArray& num_H);
244
245 /// compute the LF/HF evaluation ratio, averaged over the QoI
246 Real eval_ratio(const RealVector& sum_L_shared, const RealVector& sum_H,
247 const RealVector& sum_LL, const RealVector& sum_LH,
248 const RealVector& sum_HH, Real cost_ratio,
249 const SizetArray& N_shared, RealVector& var_H,
250 RealVector& rho2_LH);
251 /// compute the LF/HF evaluation ratio, averaged over the QoI
252 Real eval_ratio(RealMatrix& sum_L_shared, RealMatrix& sum_H,
253 RealMatrix& sum_LL, RealMatrix& sum_LH, RealMatrix& sum_HH,
254 Real cost_ratio, size_t lev, const SizetArray& N_shared,
255 RealMatrix& var_H, RealMatrix& rho2_LH);
256 /// compute the LF/HF evaluation ratio, averaged over the QoI
257 Real eval_ratio(RealMatrix& sum_Ll, RealMatrix& sum_Llm1,
258 RealMatrix& sum_Hl, RealMatrix& sum_Hlm1,
259 RealMatrix& sum_Ll_Ll, RealMatrix& sum_Ll_Llm1,
260 RealMatrix& sum_Llm1_Llm1, RealMatrix& sum_Hl_Ll,
261 RealMatrix& sum_Hl_Llm1, RealMatrix& sum_Hlm1_Ll,
262 RealMatrix& sum_Hlm1_Llm1, RealMatrix& sum_Hl_Hl,
263 RealMatrix& sum_Hl_Hlm1, RealMatrix& sum_Hlm1_Hlm1,
264 Real cost_ratio, size_t lev, const SizetArray& N_shared,
265 RealMatrix& var_YHl, RealMatrix& rho_dot2_LH);
266
267 /// compute ratio of MC and CVMC mean squared errors, averaged over the QoI
268 Real MSE_ratio(Real avg_eval_ratio, const RealVector& var_H,
269 const RealVector& rho2_LH, size_t iter,
270 const SizetArray& N_hf);
271
272 /// compute control variate parameters for CVMC and estimate raw moments
273 void cv_raw_moments(IntRealVectorMap& sum_L_shared, IntRealVectorMap& sum_H,
274 IntRealVectorMap& sum_LL, IntRealVectorMap& sum_LH,
275 const SizetArray& N_shared,
276 IntRealVectorMap& sum_L_refined,
277 const SizetArray& N_refined, const RealVector& rho2_LH,
278 RealMatrix& H_raw_mom);
279 /// apply control variate parameters for MLCVMC to estimate raw
280 /// moment contributions
281 void cv_raw_moments(IntRealMatrixMap& sum_L_shared, IntRealMatrixMap& sum_H,
282 IntRealMatrixMap& sum_LL, IntRealMatrixMap& sum_LH,
283 const SizetArray& N_shared,
284 IntRealMatrixMap& sum_L_refined,
285 const SizetArray& N_refined, const RealMatrix& rho2_LH,
286 size_t lev, RealMatrix& H_raw_mom);
287 /// apply control variate parameters for MLCVMC to estimate raw
288 /// moment contributions
289 void cv_raw_moments(IntRealMatrixMap& sum_Ll, IntRealMatrixMap& sum_Llm1,
290 IntRealMatrixMap& sum_Hl, IntRealMatrixMap& sum_Hlm1,
291 IntRealMatrixMap& sum_Ll_Ll,
292 IntRealMatrixMap& sum_Ll_Llm1,
293 IntRealMatrixMap& sum_Llm1_Llm1,
294 IntRealMatrixMap& sum_Hl_Ll,
295 IntRealMatrixMap& sum_Hl_Llm1,
296 IntRealMatrixMap& sum_Hlm1_Ll,
297 IntRealMatrixMap& sum_Hlm1_Llm1,
298 IntRealMatrixMap& sum_Hl_Hl,
299 IntRealMatrixMap& sum_Hl_Hlm1,
300 IntRealMatrixMap& sum_Hlm1_Hlm1,
301 const SizetArray& N_shared,
302 IntRealMatrixMap& sum_Ll_refined,
303 IntRealMatrixMap& sum_Llm1_refined,
304 const SizetArray& N_refined,
305 const RealMatrix& rho_dot2_LH, size_t lev,
306 RealMatrix& H_raw_mom);
307
308 /// compute scalar control variate parameters
309 void compute_control(Real sum_L, Real sum_H, Real sum_LL, Real sum_LH,
310 size_t N_shared, Real& beta);
311 /// compute scalar variance and correlation parameters for control variates
312 void compute_control(Real sum_L, Real sum_H, Real sum_LL, Real sum_LH,
313 Real sum_HH, size_t N_shared, Real& var_H,
314 Real& rho2_LH);
315 /// compute scalar control variate parameters
316 void compute_control(Real sum_Ll, Real sum_Llm1, Real sum_Hl, Real sum_Hlm1,
317 Real sum_Ll_Ll, Real sum_Ll_Llm1, Real sum_Llm1_Llm1,
318 Real sum_Hl_Ll, Real sum_Hl_Llm1, Real sum_Hlm1_Ll,
319 Real sum_Hlm1_Llm1, Real sum_Hl_Hl, Real sum_Hl_Hlm1,
320 Real sum_Hlm1_Hlm1, size_t N_shared, Real& var_YH,
321 Real& rho_dot2_LH, Real& beta_dot, Real& gamma);
322 /// compute vector control variate parameters
323 void compute_control(const RealVector& sum_L, const RealVector& sum_H,
324 const RealVector& sum_LL, const RealVector& sum_LH,
325 const SizetArray& N_shared, RealVector& beta);
326 /// compute vector variance and correlation parameters for control variates
327 void compute_control(const RealVector& sum_L, const RealVector& sum_H,
328 const RealVector& sum_LL, const RealVector& sum_LH,
329 const RealVector& sum_HH, const SizetArray& N_shared,
330 RealVector& var_H, RealVector& rho2_LH);
331 /// compute matrix control variate parameters
332 void compute_control(const RealMatrix& sum_L, const RealMatrix& sum_H,
333 const RealMatrix& sum_LL, const RealMatrix& sum_LH,
334 const SizetArray& N_shared, size_t lev,
335 RealVector& beta);
336 /// compute matrix control variate parameters
337 void compute_control(const RealMatrix& sum_Ll, const RealMatrix& sum_Llm1,
338 const RealMatrix& sum_Hl, const RealMatrix& sum_Hlm1,
339 const RealMatrix& sum_Ll_Ll,
340 const RealMatrix& sum_Ll_Llm1,
341 const RealMatrix& sum_Llm1_Llm1,
342 const RealMatrix& sum_Hl_Ll,
343 const RealMatrix& sum_Hl_Llm1,
344 const RealMatrix& sum_Hlm1_Ll,
345 const RealMatrix& sum_Hlm1_Llm1,
346 const RealMatrix& sum_Hl_Hl,
347 const RealMatrix& sum_Hl_Hlm1,
348 const RealMatrix& sum_Hlm1_Hlm1,
349 const SizetArray& N_shared, size_t lev,
350 RealVector& beta_dot, RealVector& gamma);
351
352 /// apply scalar control variate parameter (beta) to approximate HF moment
353 void apply_control(Real sum_H, Real sum_L_shared, size_t N_shared,
354 Real sum_L_refined, size_t N_refined, Real beta,
355 Real& H_raw_mom);
356 /// apply scalar control variate parameter (beta) to approximate HF moment
357 void apply_control(Real sum_Hl, Real sum_Hlm1, Real sum_Ll, Real sum_Llm1,
358 size_t N_shared, Real sum_Ll_refined,
359 Real sum_Llm1_refined, size_t N_refined, Real beta_dot,
360 Real gamma, Real& H_raw_mom);
361 /// apply vector control variate parameter (beta) to approximate HF moment
362 void apply_control(const RealVector& sum_H, const RealVector& sum_L_shared,
363 const SizetArray& N_shared,
364 const RealVector& sum_L_refined,
365 const SizetArray& N_refined, const RealVector& beta,
366 RealVector& H_raw_mom);
367 /// apply matrix control variate parameter (beta) to approximate HF moment
368 void apply_control(const RealMatrix& sum_H, const RealMatrix& sum_L_shared,
369 const SizetArray& N_shared,
370 const RealMatrix& sum_L_refined,
371 const SizetArray& N_refined, size_t lev,
372 const RealVector& beta, RealVector& H_raw_mom);
373 /// apply matrix control variate parameter (beta) to approximate HF moment
374 void apply_control(const RealMatrix& sum_Hl, const RealMatrix& sum_Hlm1,
375 const RealMatrix& sum_Ll, const RealMatrix& sum_Llm1,
376 const SizetArray& N_shared,
377 const RealMatrix& sum_Ll_refined,
378 const RealMatrix& sum_Llm1_refined,
379 const SizetArray& N_refined, size_t lev,
380 const RealVector& beta_dot, const RealVector& gamma,
381 RealVector& H_raw_mom);
382
383 /// export allSamples to tagged tabular file
384 void export_all_samples(String root_prepend, const Model& model,
385 size_t iter, size_t lev);
386
387 /// convert uncentered raw moments (multilevel expectations) to
388 /// standardized moments
389 void convert_moments(const RealMatrix& raw_mom, RealMatrix& final_mom);
390
391 /// populate finalStatErrors for MLMC based on Q sums
392 void compute_error_estimates(IntRealMatrixMap& sum_Ql,
393 IntRealMatrixMap& sum_Qlm1,
394 IntIntPairRealMatrixMap& sum_QlQlm1,
395 Sizet2DArray& num_Q);
396
397 /// compute variance from sum accumulators
398 Real variance_Ysum(Real sum_Y, Real sum_YY, size_t Nlq);
399 /// compute variance from sum accumulators
400 Real variance_Qsum(Real sum_Ql, Real sum_Qlm1, Real sum_QlQl, Real sum_QlQlm1,
401 Real sum_Qlm1Qlm1, size_t Nlq);
402
403 /// sum up variances across QoI (using sum_YY with means from sum_Y)
404 Real aggregate_variance_Ysum(const Real* sum_Y, const Real* sum_YY,
405 const SizetArray& N_l);
406 /// sum up variances across QoI (using sum_YY with means from sum_Y)
407 Real aggregate_variance_Qsum(const Real* sum_Ql, const Real* sum_Qlm1,
408 const Real* sum_QlQl, const Real* sum_QlQlm1,
409 const Real* sum_Qlm1Qlm1, const SizetArray& N_l,
410 const size_t& lev);
411 Real aggregate_variance_Qsum(const Real* sum_Ql, const Real* sum_Qlm1,
412 const Real* sum_QlQl, const Real* sum_QlQlm1,
413 const Real* sum_Qlm1Qlm1, const SizetArray& N_l,
414 const size_t& lev, const size_t& qoi);
415
416 /// sum up Monte Carlo estimates for mean squared error (MSE) across
417 /// QoI using discrepancy variances
418 Real aggregate_mse_Yvar(const Real* var_Y, const SizetArray& N_l);
419 /// sum up Monte Carlo estimates for mean squared error (MSE) across
420 /// QoI using discrepancy sums
421 Real aggregate_mse_Ysum(const Real* sum_Y, const Real* sum_YY,
422 const SizetArray& N_l);
423 /// sum up Monte Carlo estimates for mean squared error (MSE) across
424 /// QoI using discrepancy sums
425 Real aggregate_mse_Qsum(const Real* sum_Ql, const Real* sum_Qlm1,
426 const Real* sum_QlQl, const Real* sum_QlQlm1,
427 const Real* sum_Qlm1Qlm1, const SizetArray& N_l,
428 const size_t& lev);
429
430 Real aggregate_mse_Qsum(const Real* sum_Ql, const Real* sum_Qlm1,
431 const Real* sum_QlQl, const Real* sum_QlQlm1,
432 const Real* sum_Qlm1Qlm1, const SizetArray& N_l,
433 const size_t& lev, const size_t& qoi);
434
435 /// convert uncentered (raw) moments to centered moments; biased estimators
436 static void uncentered_to_centered(Real rm1, Real rm2, Real rm3, Real rm4,
437 Real& cm1, Real& cm2, Real& cm3, Real& cm4,
438 size_t Nlq);
439 /// convert uncentered (raw) moments to centered moments; unbiased estimators
440 static void uncentered_to_centered(Real rm1, Real rm2, Real rm3, Real rm4,
441 Real& cm1, Real& cm2, Real& cm3, Real& cm4);
442 /// convert centered moments to standardized moments
443 static void centered_to_standard(Real cm1, Real cm2, Real cm3, Real cm4,
444 Real& sm1, Real& sm2, Real& sm3, Real& sm4);
445 /// detect, warn, and repair a negative central moment (for even orders)
446 static void check_negative(Real& cm);
447
448 /// compute sum of a set of observations
449 Real sum(const Real* vec, size_t vec_len) const;
450 /// compute average of a set of observations
451 Real average(const Real* vec, size_t vec_len) const;
452 /// compute average of a set of observations
453 Real average(const RealVector& vec) const;
454 /// compute average of a set of observations
455 Real average(const SizetArray& sa) const;
456
457 /// compute the unbiased product of two sampling means
458 static Real unbiased_mean_product_pair(const Real& sumQ1, const Real& sumQ2, const Real& sumQ1Q2, const size_t& Nlq);
459 /// compute the unbiased product of three sampling means
460 static Real unbiased_mean_product_triplet(const Real& sumQ1, const Real& sumQ2, const Real& sumQ3,
461 const Real& sumQ1Q2, const Real& sumQ1Q3, const Real& sumQ2Q3, const Real& sumQ1Q2Q3, const size_t& Nlq);
462 /// compute the unbiased product of two pairs of products of sampling means
463 static Real unbiased_mean_product_pairpair(const Real& sumQ1, const Real& sumQ2, const Real& sumQ1Q2,
464 const Real& sumQ1sq, const Real& sumQ2sq,
465 const Real& sumQ1sqQ2, const Real& sumQ1Q2sq, const Real& sumQ1sqQ2sq, const size_t& Nlq);
466
467 static Real var_of_var_ml_l0(IntRealMatrixMap sum_Ql, IntRealMatrixMap sum_Qlm1, IntIntPairRealMatrixMap sum_QlQlm1, const size_t& Nlq_pilot, const Real& Nlq, const size_t& qoi, const bool& compute_gradient, Real& grad_g);
468
469 static Real var_of_var_ml_lmax(IntRealMatrixMap sum_Ql, IntRealMatrixMap sum_Qlm1, IntIntPairRealMatrixMap sum_QlQlm1, const size_t& Nlq_pilot, const Real& Nlq, const size_t& qoi, const bool& compute_gradient, Real& grad_g);
470
471 static Real var_of_var_ml_l(IntRealMatrixMap sum_Ql, IntRealMatrixMap sum_Qlm1, IntIntPairRealMatrixMap sum_QlQlm1, const size_t& Nlq_pilot, const Real& Nlq, const size_t& qoi, const size_t& lev, const bool& compute_gradient, Real& grad_g);
472
473 ///OPTPP definition
474 static void target_var_objective_eval_optpp(int mode, int n, const RealVector& x, double& f,
475 RealVector& grad_f, int& result_mode);
476 static void target_var_constraint_eval_optpp(int mode, int n, const RealVector& x, RealVector& g,
477 RealMatrix& grad_g, int& result_mode);
478 static void target_var_constraint_eval_logscale_optpp(int mode, int n, const RealVector& x, RealVector& g,
479 RealMatrix& grad_g, int& result_mode);
480 /// NPSOL definition (Wrapper using OPTPP implementation above under the hood)
481 static void target_var_objective_eval_npsol(int& mode, int& n, double* x, double& f, double* gradf, int& nstate);
482 static 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);
483 static void 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);
484
485 void assign_static_member(Real &conv_tol, size_t &qoi, RealVector &level_cost_vec, IntRealMatrixMap &sum_Ql,
486 IntRealMatrixMap &sum_Qlm1, IntIntPairRealMatrixMap &sum_QlQlm1,
487 RealVector &pilot_samples) const;
488
489 void assign_static_member_problem18(Real &var_L_exact, Real &var_H_exact, Real &mu_four_L_exact, Real &mu_four_H_exact, Real &Ax, RealVector &level_cost_vec) const;
490
491 //
492 //- Heading: Data
493 //
494
495 /// total number of successful sample evaluations (excluding faults)
496 /// for each model form, discretization level, and QoI
497 Sizet3DArray NLev;
498
499 /// store the pilot_samples input specification, prior to run-time
500 /// invocation of load_pilot_sample()
501 SizetArray pilotSamples;
502
503 /// user specification for seed_sequence
504 SizetArray randomSeedSeqSpec;
505
506 /// major iteration counter
507 size_t mlmfIter;
508
509 /// store the allocation_target input specification, prior to run-time
510 /// Options right now:
511 /// - Mean = First moment (Mean)
512 /// - Variance = Second moment (Variance or standard deviation depending on moments central or standard)
513 short allocationTarget;
514
515 /// option to switch on numerical optimization for solution of sample alloation
516 /// of allocationTarget Variance
517 bool useTargetVarianceOptimizationFlag;
518
519 /// store the qoi_aggregation_norm input specification, prior to run-time
520 /// Options right now:
521 /// - sum = aggregate the variance over all QoIs, compute samples from that
522 /// - max = take maximum sample allocation over QoIs for each level
523 short qoiAggregation;
524
525 /// mean squared error of mean estimator from pilot sample MC on HF model
526 RealVector mcMSEIter0;
527 /// equivalent number of high fidelity evaluations accumulated using samples
528 /// across multiple model forms and/or discretization levels
529 Real equivHFEvals;
530
531 /// if defined, complete the final CV refinement when terminating MLCV based
532 /// on maxIterations (the total number of refinements beyond the pilot sample
533 /// will be one more for CV than for ML). This approach is consistent with
534 /// normal termination based on l1_norm(delta_N_hf) = 0.
535 bool finalCVRefinement;
536
537 /// if defined, export each of the sample increments in ML, CV, MLCV
538 /// using tagged tabular files
539 bool exportSampleSets;
540 /// format for exporting sample increments using tagged tabular files
541 unsigned short exportSamplesFormat;
542 };
543
544
aggregated_models_mode()545 inline void NonDMultilevelSampling::aggregated_models_mode()
546 {
547 if (iteratedModel.surrogate_response_mode() != AGGREGATED_MODELS) {
548 iteratedModel.surrogate_response_mode(AGGREGATED_MODELS); // set HF,LF
549 // synch activeSet with iteratedModel.response_size()
550 activeSet.reshape(2*numFunctions);
551 activeSet.request_values(1);
552 }
553 }
554
555
bypass_surrogate_mode()556 inline void NonDMultilevelSampling::bypass_surrogate_mode()
557 {
558 if (iteratedModel.surrogate_response_mode() != BYPASS_SURROGATE) {
559 iteratedModel.surrogate_response_mode(BYPASS_SURROGATE); // HF
560 activeSet.reshape(numFunctions);// synch with model.response_size()
561 }
562 }
563
564
uncorrected_surrogate_mode()565 inline void NonDMultilevelSampling::uncorrected_surrogate_mode()
566 {
567 if (iteratedModel.surrogate_response_mode() != UNCORRECTED_SURROGATE) {
568 iteratedModel.surrogate_response_mode(UNCORRECTED_SURROGATE); // LF
569 activeSet.reshape(numFunctions);// synch with model.response_size()
570 }
571 }
572
573
574 inline Real NonDMultilevelSampling::
level_cost(const RealVector & cost,unsigned short step)575 level_cost(const RealVector& cost, unsigned short step)
576 {
577 // discrepancies incur two level costs
578 return (step) ?
579 cost[step] + cost[step-1] : // aggregated {HF,LF} mode
580 cost[step]; // uncorrected LF mode
581 }
582
583
584 /** extract an active seed from a seed sequence */
random_seed(size_t index) const585 inline int NonDMultilevelSampling::random_seed(size_t index) const
586 {
587 // return 0 for cases where seed is undefined or will not be updated
588
589 if (randomSeedSeqSpec.empty()) return 0; // no spec -> non-repeatable samples
590 else if (!varyPattern) // continually reset seed to specified value
591 return (index < randomSeedSeqSpec.size()) ?
592 randomSeedSeqSpec[index] : randomSeedSeqSpec.back();
593 // only set sequence of seeds for first pass, then let RNG state continue
594 else if (mlmfIter == 0 && index < randomSeedSeqSpec.size()) // pilot iter only
595 return randomSeedSeqSpec[index];
596 else return 0; // seed sequence exhausted, do not update
597 }
598
599
600 inline void NonDMultilevelSampling::
compute_control(const RealVector & sum_L,const RealVector & sum_H,const RealVector & sum_LL,const RealVector & sum_LH,const SizetArray & N_shared,RealVector & beta)601 compute_control(const RealVector& sum_L, const RealVector& sum_H,
602 const RealVector& sum_LL, const RealVector& sum_LH,
603 const SizetArray& N_shared, RealVector& beta)
604 {
605 for (size_t qoi=0; qoi<numFunctions; ++qoi)
606 compute_control(sum_L[qoi], sum_H[qoi], sum_LL[qoi], sum_LH[qoi],
607 N_shared[qoi], beta[qoi]);
608 }
609
610
611 inline void NonDMultilevelSampling::
compute_control(const RealVector & sum_L,const RealVector & sum_H,const RealVector & sum_LL,const RealVector & sum_LH,const RealVector & sum_HH,const SizetArray & N_shared,RealVector & var_H,RealVector & rho2_LH)612 compute_control(const RealVector& sum_L, const RealVector& sum_H,
613 const RealVector& sum_LL, const RealVector& sum_LH,
614 const RealVector& sum_HH, const SizetArray& N_shared,
615 RealVector& var_H, RealVector& rho2_LH)
616 {
617 for (size_t qoi=0; qoi<numFunctions; ++qoi)
618 compute_control(sum_L[qoi], sum_H[qoi], sum_LL[qoi], sum_LH[qoi],
619 sum_HH[qoi], N_shared[qoi], var_H[qoi], rho2_LH[qoi]);
620 }
621
622
623 inline void NonDMultilevelSampling::
compute_control(const RealMatrix & sum_L,const RealMatrix & sum_H,const RealMatrix & sum_LL,const RealMatrix & sum_LH,const SizetArray & N_shared,size_t lev,RealVector & beta)624 compute_control(const RealMatrix& sum_L, const RealMatrix& sum_H,
625 const RealMatrix& sum_LL, const RealMatrix& sum_LH,
626 const SizetArray& N_shared, size_t lev, RealVector& beta)
627 {
628 for (size_t qoi=0; qoi<numFunctions; ++qoi)
629 compute_control(sum_L(qoi,lev), sum_H(qoi,lev), sum_LL(qoi,lev),
630 sum_LH(qoi,lev), N_shared[qoi], beta[qoi]);
631 }
632
633
634 inline void NonDMultilevelSampling::
compute_control(const RealMatrix & sum_Ll,const RealMatrix & sum_Llm1,const RealMatrix & sum_Hl,const RealMatrix & sum_Hlm1,const RealMatrix & sum_Ll_Ll,const RealMatrix & sum_Ll_Llm1,const RealMatrix & sum_Llm1_Llm1,const RealMatrix & sum_Hl_Ll,const RealMatrix & sum_Hl_Llm1,const RealMatrix & sum_Hlm1_Ll,const RealMatrix & sum_Hlm1_Llm1,const RealMatrix & sum_Hl_Hl,const RealMatrix & sum_Hl_Hlm1,const RealMatrix & sum_Hlm1_Hlm1,const SizetArray & N_shared,size_t lev,RealVector & beta_dot,RealVector & gamma)635 compute_control(const RealMatrix& sum_Ll, const RealMatrix& sum_Llm1,
636 const RealMatrix& sum_Hl, const RealMatrix& sum_Hlm1,
637 const RealMatrix& sum_Ll_Ll, const RealMatrix& sum_Ll_Llm1,
638 const RealMatrix& sum_Llm1_Llm1, const RealMatrix& sum_Hl_Ll,
639 const RealMatrix& sum_Hl_Llm1, const RealMatrix& sum_Hlm1_Ll,
640 const RealMatrix& sum_Hlm1_Llm1, const RealMatrix& sum_Hl_Hl,
641 const RealMatrix& sum_Hl_Hlm1, const RealMatrix& sum_Hlm1_Hlm1,
642 const SizetArray& N_shared, size_t lev, RealVector& beta_dot,
643 RealVector& gamma)
644 {
645 Real var_YH, rho_dot2_LH; // not needed for this context
646 for (size_t qoi=0; qoi<numFunctions; ++qoi)
647 compute_control(sum_Ll(qoi,lev), sum_Llm1(qoi,lev), sum_Hl(qoi,lev),
648 sum_Hlm1(qoi,lev), sum_Ll_Ll(qoi,lev), sum_Ll_Llm1(qoi,lev),
649 sum_Llm1_Llm1(qoi,lev), sum_Hl_Ll(qoi,lev),
650 sum_Hl_Llm1(qoi,lev), sum_Hlm1_Ll(qoi,lev),
651 sum_Hlm1_Llm1(qoi,lev), sum_Hl_Hl(qoi,lev),
652 sum_Hl_Hlm1(qoi,lev), sum_Hlm1_Hlm1(qoi,lev),
653 N_shared[qoi], var_YH, rho_dot2_LH,
654 beta_dot[qoi], gamma[qoi]);
655 }
656
657
658 inline void NonDMultilevelSampling::
apply_control(const RealVector & sum_H,const RealVector & sum_L_shared,const SizetArray & N_shared,const RealVector & sum_L_refined,const SizetArray & N_refined,const RealVector & beta,RealVector & H_raw_mom)659 apply_control(const RealVector& sum_H, const RealVector& sum_L_shared,
660 const SizetArray& N_shared, const RealVector& sum_L_refined,
661 const SizetArray& N_refined, const RealVector& beta,
662 RealVector& H_raw_mom)
663 {
664 for (size_t qoi=0; qoi<numFunctions; ++qoi) {
665 Cout << " QoI " << qoi+1 << ": control variate beta = "
666 << std::setw(9) << beta[qoi] << '\n';
667 apply_control(sum_H[qoi], sum_L_shared[qoi], N_shared[qoi],
668 sum_L_refined[qoi], N_refined[qoi], beta[qoi],
669 H_raw_mom[qoi]);
670 }
671 if (numFunctions > 1) Cout << '\n';
672 }
673
674
675 inline void NonDMultilevelSampling::
apply_control(const RealMatrix & sum_H,const RealMatrix & sum_L_shared,const SizetArray & N_shared,const RealMatrix & sum_L_refined,const SizetArray & N_refined,size_t lev,const RealVector & beta,RealVector & H_raw_mom)676 apply_control(const RealMatrix& sum_H, const RealMatrix& sum_L_shared,
677 const SizetArray& N_shared, const RealMatrix& sum_L_refined,
678 const SizetArray& N_refined, size_t lev, const RealVector& beta,
679 RealVector& H_raw_mom)
680 {
681 for (size_t qoi=0; qoi<numFunctions; ++qoi) {
682 Cout << " QoI " << qoi+1 << ": control variate beta = "
683 << std::setw(9) << beta[qoi] << '\n';
684 apply_control(sum_H(qoi,lev), sum_L_shared(qoi,lev), N_shared[qoi],
685 sum_L_refined(qoi,lev), N_refined[qoi], beta[qoi],
686 H_raw_mom[qoi]);
687 }
688 if (numFunctions > 1) Cout << '\n';
689 }
690
691
692 inline void NonDMultilevelSampling::
apply_control(const RealMatrix & sum_Hl,const RealMatrix & sum_Hlm1,const RealMatrix & sum_Ll,const RealMatrix & sum_Llm1,const SizetArray & N_shared,const RealMatrix & sum_Ll_refined,const RealMatrix & sum_Llm1_refined,const SizetArray & N_refined,size_t lev,const RealVector & beta_dot,const RealVector & gamma,RealVector & H_raw_mom)693 apply_control(const RealMatrix& sum_Hl, const RealMatrix& sum_Hlm1,
694 const RealMatrix& sum_Ll, const RealMatrix& sum_Llm1,
695 const SizetArray& N_shared, const RealMatrix& sum_Ll_refined,
696 const RealMatrix& sum_Llm1_refined, const SizetArray& N_refined,
697 size_t lev, const RealVector& beta_dot, const RealVector& gamma,
698 RealVector& H_raw_mom)
699 {
700 for (size_t qoi=0; qoi<numFunctions; ++qoi) {
701 Cout << " QoI " << qoi+1 << ": control variate beta_dot = "
702 << std::setw(9) << beta_dot[qoi] << '\n';
703 apply_control(sum_Hl(qoi,lev), sum_Hlm1(qoi,lev), sum_Ll(qoi,lev),
704 sum_Llm1(qoi,lev), N_shared[qoi], sum_Ll_refined(qoi,lev),
705 sum_Llm1_refined(qoi,lev), N_refined[qoi], beta_dot[qoi],
706 gamma[qoi], H_raw_mom[qoi]);
707 }
708 if (numFunctions > 1) Cout << '\n';
709 }
710
711
712 inline Real NonDMultilevelSampling::
variance_Ysum(Real sum_Y,Real sum_YY,size_t Nlq)713 variance_Ysum(Real sum_Y, Real sum_YY, /*Real offset,*/ size_t Nlq)
714 {
715 Real mu_Y = sum_Y / Nlq;
716 // Note: precision loss in variance is difficult to avoid without
717 // storing full sample history; must accumulate Y^2 across iterations
718 // instead of (Y-mean)^2 since mean is updated on each iteration.
719 Real var_Y = (sum_YY / Nlq - mu_Y * mu_Y)
720 * (Real)Nlq / (Real)(Nlq - 1); // Bessel's correction
721 return var_Y;
722
723 /*
724 Real new_mu_Y = mu_Y + offset;
725 return var_Y
726 // + offset * offset // uncenter from old mu_hat
727 // - new_mu_Y * new_mu_Y; // recenter with new_mu_Y
728 - mu_Y * mu_Y - 2. * mu_Y * offset; // cancel offset^2
729 */
730 }
731
732
733 inline Real NonDMultilevelSampling::
variance_Qsum(Real sum_Ql,Real sum_Qlm1,Real sum_QlQl,Real sum_QlQlm1,Real sum_Qlm1Qlm1,size_t Nlq)734 variance_Qsum(Real sum_Ql, Real sum_Qlm1, Real sum_QlQl, Real sum_QlQlm1,
735 Real sum_Qlm1Qlm1, size_t Nlq)
736 {
737 Real mu_Ql = sum_Ql / Nlq, mu_Qlm1 = sum_Qlm1 / Nlq;
738 //var_Y = var_Ql - 2.* covar_QlQlm1 + var_Qlm1;
739 return ( sum_QlQl / Nlq - mu_Ql * mu_Ql // var_Ql
740 - 2. * ( sum_QlQlm1 / Nlq - mu_Ql * mu_Qlm1 ) // covar_QlQlm1
741 + sum_Qlm1Qlm1 / Nlq - mu_Qlm1 * mu_Qlm1 ) // var_Qlm1
742 * (Real)Nlq / (Real)(Nlq - 1); // Bessel's correction
743 }
744
745
746 inline Real NonDMultilevelSampling::
aggregate_variance_Ysum(const Real * sum_Y,const Real * sum_YY,const SizetArray & N_l)747 aggregate_variance_Ysum(const Real* sum_Y, const Real* sum_YY,
748 const SizetArray& N_l)
749 {
750 Real agg_var_l = 0.;//, var_Y;
751 //if (outputLevel >= DEBUG_OUTPUT) Cout << "[ ";
752 for (size_t qoi=0; qoi<numFunctions; ++qoi) //{
753 agg_var_l += variance_Ysum(sum_Y[qoi], sum_YY[qoi], N_l[qoi]);
754 //if (outputLevel >= DEBUG_OUTPUT) Cout << var_Y << ' ';
755 //}
756 //if (outputLevel >= DEBUG_OUTPUT) Cout << "]\n";
757 return agg_var_l;
758 }
759
760
761 inline Real NonDMultilevelSampling::
aggregate_variance_Qsum(const Real * sum_Ql,const Real * sum_Qlm1,const Real * sum_QlQl,const Real * sum_QlQlm1,const Real * sum_Qlm1Qlm1,const SizetArray & N_l,const size_t & lev)762 aggregate_variance_Qsum(const Real* sum_Ql, const Real* sum_Qlm1,
763 const Real* sum_QlQl, const Real* sum_QlQlm1,
764 const Real* sum_Qlm1Qlm1, const SizetArray& N_l,
765 const size_t& lev)
766 {
767 Real agg_var_l = 0., var_Y;
768 //if (outputLevel >= DEBUG_OUTPUT) Cout << "[ ";
769 for (size_t qoi=0; qoi<numFunctions; ++qoi) //{
770 agg_var_l += (lev) ?
771 variance_Qsum(sum_Ql[qoi], sum_Qlm1[qoi], sum_QlQl[qoi], sum_QlQlm1[qoi],
772 sum_Qlm1Qlm1[qoi], N_l[qoi]) :
773 variance_Ysum(sum_Ql[qoi], sum_QlQl[qoi], N_l[qoi]);
774 //if (outputLevel >= DEBUG_OUTPUT) Cout << var_Y << ' ';
775 //}
776 //if (outputLevel >= DEBUG_OUTPUT) Cout << "]\n";
777 return agg_var_l;
778 }
779
780 inline Real NonDMultilevelSampling::
aggregate_variance_Qsum(const Real * sum_Ql,const Real * sum_Qlm1,const Real * sum_QlQl,const Real * sum_QlQlm1,const Real * sum_Qlm1Qlm1,const SizetArray & N_l,const size_t & lev,const size_t & qoi)781 aggregate_variance_Qsum(const Real* sum_Ql, const Real* sum_Qlm1,
782 const Real* sum_QlQl, const Real* sum_QlQlm1,
783 const Real* sum_Qlm1Qlm1, const SizetArray& N_l,
784 const size_t& lev, const size_t& qoi)
785 {
786 Real agg_var_l = 0., var_Y;
787 //if (outputLevel >= DEBUG_OUTPUT) Cout << "[ ";
788 //for (size_t qoi=0; qoi<numFunctions; ++qoi) //{
789 agg_var_l = (lev) ?
790 variance_Qsum(sum_Ql[qoi], sum_Qlm1[qoi], sum_QlQl[qoi], sum_QlQlm1[qoi],
791 sum_Qlm1Qlm1[qoi], N_l[qoi]) :
792 variance_Ysum(sum_Ql[qoi], sum_QlQl[qoi], N_l[qoi]);
793 //if (outputLevel >= DEBUG_OUTPUT) Cout << var_Y << ' ';
794 //}
795 //if (outputLevel >= DEBUG_OUTPUT) Cout << "]\n";
796 return agg_var_l;
797 }
798
799
800 inline Real NonDMultilevelSampling::
aggregate_mse_Yvar(const Real * var_Y,const SizetArray & N_l)801 aggregate_mse_Yvar(const Real* var_Y, const SizetArray& N_l)
802 {
803 Real agg_mse = 0.;
804 // aggregate MC estimator variance for each QoI
805 for (size_t qoi=0; qoi<numFunctions; ++qoi)
806 agg_mse += var_Y[qoi] / N_l[qoi];
807 return agg_mse;
808 }
809
810
811 inline Real NonDMultilevelSampling::
aggregate_mse_Ysum(const Real * sum_Y,const Real * sum_YY,const SizetArray & N_l)812 aggregate_mse_Ysum(const Real* sum_Y, const Real* sum_YY, const SizetArray& N_l)
813 {
814 Real agg_mse = 0.;//, var_Y;
815 // aggregate MC estimator variance for each QoI
816 for (size_t qoi=0; qoi<numFunctions; ++qoi)
817 agg_mse += variance_Ysum(sum_Y[qoi], sum_YY[qoi], N_l[qoi]) / N_l[qoi];
818 return agg_mse;
819 }
820
821
822 inline Real NonDMultilevelSampling::
aggregate_mse_Qsum(const Real * sum_Ql,const Real * sum_Qlm1,const Real * sum_QlQl,const Real * sum_QlQlm1,const Real * sum_Qlm1Qlm1,const SizetArray & N_l,const size_t & lev)823 aggregate_mse_Qsum(const Real* sum_Ql, const Real* sum_Qlm1,
824 const Real* sum_QlQl, const Real* sum_QlQlm1,
825 const Real* sum_Qlm1Qlm1, const SizetArray& N_l, const size_t& lev)
826 {
827 Real agg_mse = 0., mu_Ql, mu_Qlm1, var_Y; size_t Nlq;
828 for (size_t qoi=0; qoi<numFunctions; ++qoi) {
829 Nlq = N_l[qoi];
830 var_Y = (lev) ?
831 variance_Qsum(sum_Ql[qoi], sum_Qlm1[qoi], sum_QlQl[qoi], sum_QlQlm1[qoi],
832 sum_Qlm1Qlm1[qoi], Nlq) :
833 variance_Ysum(sum_Ql[qoi], sum_QlQl[qoi], Nlq);
834 agg_mse += var_Y / Nlq; // aggregate MC estimator variance for each QoI
835 }
836 return agg_mse;
837 }
838
839 inline Real NonDMultilevelSampling::
aggregate_mse_Qsum(const Real * sum_Ql,const Real * sum_Qlm1,const Real * sum_QlQl,const Real * sum_QlQlm1,const Real * sum_Qlm1Qlm1,const SizetArray & N_l,const size_t & lev,const size_t & qoi)840 aggregate_mse_Qsum(const Real* sum_Ql, const Real* sum_Qlm1,
841 const Real* sum_QlQl, const Real* sum_QlQlm1,
842 const Real* sum_Qlm1Qlm1, const SizetArray& N_l, const size_t& lev, const size_t& qoi)
843 {
844 Real agg_mse = 0., mu_Ql, mu_Qlm1, var_Y; size_t Nlq;
845 Nlq = N_l[qoi];
846 var_Y = (lev) ?
847 variance_Qsum(sum_Ql[qoi], sum_Qlm1[qoi], sum_QlQl[qoi], sum_QlQlm1[qoi],
848 sum_Qlm1Qlm1[qoi], Nlq) :
849 variance_Ysum(sum_Ql[qoi], sum_QlQl[qoi], Nlq);
850 agg_mse += var_Y / Nlq; // aggregate MC estimator variance for each QoI
851
852 return agg_mse;
853 }
854
855
accumulate_offsets(RealVector & mu)856 inline void NonDMultilevelSampling::accumulate_offsets(RealVector& mu)
857 {
858 IntRespMCIter r_it = allResponses.begin();
859 size_t qoi, num_samp, num_fns = r_it->second.num_functions();
860 mu.sizeUninitialized(num_fns);
861 Real q_l, sum;
862 for (qoi=0; qoi<num_fns; ++qoi) {
863 num_samp = 0; sum = 0.;
864 for (r_it=allResponses.begin(); r_it!=allResponses.end(); ++r_it) {
865 q_l = r_it->second.function_value(qoi);
866 if (std::isfinite(q_l)) // neither NaN nor +/-Inf
867 { sum += q_l; ++num_samp; }
868 }
869 mu[qoi] = sum / num_samp;
870 }
871 }
872
873
874 /** For single-level moment calculations with a scalar Nlq. */
875 inline void NonDMultilevelSampling::
uncentered_to_centered(Real rm1,Real rm2,Real rm3,Real rm4,Real & cm1,Real & cm2,Real & cm3,Real & cm4)876 uncentered_to_centered(Real rm1, Real rm2, Real rm3, Real rm4,
877 Real& cm1, Real& cm2, Real& cm3, Real& cm4)
878 {
879 // convert from uncentered ("raw") to centered moments for a single level
880
881 // For moments from sampling:
882 // > Raw moments are unbiased. Central moments are unbiased for an exact mean
883 // (i.e., the samples represent the full population).
884 // > For sampling a portion of the population, central moments {2,3,4} are
885 // biased estimators since the mean is approximated. The conversion to
886 // unbiased requires a correction based on the number of samples, as
887 // implemented in the subsequent function.
888
889 cm1 = rm1; // mean
890 cm2 = rm2 - cm1 * cm1; // variance
891
892 cm3 = rm3 - cm1 * (3. * cm2 + cm1 * cm1); // using cm
893 // = rm3 - cm1 * (3. * rm2 - 2. * cm1 * cm1); // using rm
894
895 // the 4th moment is the central moment (non-excess, not cumulant)
896 cm4 = rm4 - cm1 * (4. * cm3 + cm1 * (6. * cm2 + cm1 * cm1)); // using cm
897 // = rm4 - cm1 * (4. * rm3 - cm1 * (6. * rm2 - 3. * cm1 * cm1)); // using rm
898 }
899
900
901 /** For single-level moment calculations with a scalar Nlq. */
902 inline void NonDMultilevelSampling::
uncentered_to_centered(Real rm1,Real rm2,Real rm3,Real rm4,Real & cm1,Real & cm2,Real & cm3,Real & cm4,size_t Nlq)903 uncentered_to_centered(Real rm1, Real rm2, Real rm3, Real rm4, Real& cm1,
904 Real& cm2, Real& cm3, Real& cm4, size_t Nlq)
905 {
906 // convert from uncentered ("raw") to centered moments for a single level
907
908 // Biased central moment estimators:
909 uncentered_to_centered(rm1, rm2, rm3, rm4, cm1, cm2, cm3, cm4);
910
911 // Bias corrections for population-based estimators w/ estimated means:
912 if (Nlq > 3) {
913 Real cm1_sq = cm1 * cm1;
914 Real nm1 = Nlq - 1., nm2 = Nlq - 2., n_sq = Nlq * Nlq;
915 cm2 *= Nlq / nm1; // unbiased population variance from Bessel's correction
916 cm3 *= n_sq / (nm1 * nm2);
917 // From "Modeling with Data," Klemens 2009 (Appendix M).
918 // Notes:
919 // (1) this is 4th central moment (non-excess, unnormalized),
920 // which differs from the fourth cumulant (excess, unnormalized)
921 // (2) cm2 is now unbiased within following conversion:
922
923 //cm4 = ( n_sq * Nlq * cm4 / nm1 - (6. * Nlq - 9.) * cm2 * cm2 )
924 // / (n_sq - 3. * Nlq + 3.);
925 //[fm] account for bias correction due to cm2^2 term
926 cm4 = ( n_sq * Nlq * cm4 / nm1 - (6. * Nlq - 9.) * (n_sq - Nlq) / (n_sq - 2. * Nlq + 3) * cm2 * cm2 )
927 / ( (n_sq - 3. * Nlq + 3.) - (6. * Nlq - 9.) * (n_sq - Nlq) / (Nlq * (n_sq - 2. * Nlq + 3.)) );
928 }
929 else
930 Cerr << "Warning: due to small sample size, resorting to biased estimator "
931 << "conversion in NonDMultilevelSampling::uncentered_to_centered().\n";
932 }
933
934
935 inline void NonDMultilevelSampling::
centered_to_standard(Real cm1,Real cm2,Real cm3,Real cm4,Real & sm1,Real & sm2,Real & sm3,Real & sm4)936 centered_to_standard(Real cm1, Real cm2, Real cm3, Real cm4,
937 Real& sm1, Real& sm2, Real& sm3, Real& sm4)
938 {
939 // convert from centered to standardized moments
940 sm1 = cm1; // mean
941 if (cm2 > 0.) {
942 sm2 = std::sqrt(cm2); // std deviation
943 sm3 = cm3 / (cm2 * sm2); // skewness
944 sm4 = cm4 / (cm2 * cm2) - 3.; // excess kurtosis
945 }
946 else {
947 Cerr << "\nWarning: central to standard conversion failed due to "
948 << "non-positive\n variance. Retaining central moments.\n";
949 sm2 = 0.; sm3 = cm3; sm4 = cm4; // or assign NaN to sm{3,4}
950 }
951 }
952
953
check_negative(Real & cm)954 inline void NonDMultilevelSampling::check_negative(Real& cm)
955 {
956 if (cm < 0.) {
957 Cerr << "\nWarning: central moment less than zero (" << cm << "). "
958 << "Repairing to zero.\n";
959 cm = 0.;
960 // TO DO: consider hard error if COV < -tol (pass in mean and cm order)
961 }
962 }
963
964
sum(const Real * vec,size_t vec_len) const965 inline Real NonDMultilevelSampling::sum(const Real* vec, size_t vec_len) const
966 {
967 Real sum = 0.;
968 for (size_t i=0; i<vec_len; ++i)
969 sum += vec[i];
970 return sum;
971 }
972
973
974 inline Real NonDMultilevelSampling::
average(const Real * vec,size_t vec_len) const975 average(const Real* vec, size_t vec_len) const
976 { return sum(vec, vec_len) / (Real)vec_len; }
977
978
average(const RealVector & vec) const979 inline Real NonDMultilevelSampling::average(const RealVector& vec) const
980 { return average(vec.values(), vec.length()); }
981
982
average(const SizetArray & sa) const983 inline Real NonDMultilevelSampling::average(const SizetArray& sa) const
984 {
985 size_t i, len = sa.size(), sum = 0;
986 for (i=0; i<len; ++i)
987 sum += sa[i];
988 return (Real)sum / (Real)len;
989 }
990
var_of_var_ml_l0(IntRealMatrixMap sum_Ql,IntRealMatrixMap sum_Qlm1,IntIntPairRealMatrixMap sum_QlQlm1,const size_t & Nlq_pilot,const Real & Nlq,const size_t & qoi,const bool & compute_gradient,Real & grad_g)991 inline Real NonDMultilevelSampling::var_of_var_ml_l0(IntRealMatrixMap sum_Ql, IntRealMatrixMap sum_Qlm1, IntIntPairRealMatrixMap sum_QlQlm1,
992 const size_t& Nlq_pilot, const Real& Nlq, const size_t& qoi, const bool& compute_gradient, Real& grad_g){
993 Real cm1l, cm2l, cm3l, cm4l, cm2l_sq, var_of_var;
994
995 IntIntPair pr11(1, 1), pr12(1, 2), pr21(2, 1), pr22(2, 2);
996 RealMatrix &sum_Q1l = sum_Ql[1], &sum_Q1lm1 = sum_Qlm1[1],
997 &sum_Q2l = sum_Ql[2], &sum_Q2lm1 = sum_Qlm1[2],
998 &sum_Q3l = sum_Ql[3], &sum_Q3lm1 = sum_Qlm1[3],
999 &sum_Q4l = sum_Ql[4], &sum_Q4lm1 = sum_Qlm1[4],
1000 &sum_Q1lQ1lm1 = sum_QlQlm1[pr11], &sum_Q1lQ2lm1 = sum_QlQlm1[pr12],
1001 &sum_Q2lQ1lm1 = sum_QlQlm1[pr21], &sum_Q2lQ2lm1 = sum_QlQlm1[pr22];
1002
1003 uncentered_to_centered(sum_Q1l(qoi, 0) / Nlq_pilot, sum_Q2l(qoi, 0) / Nlq_pilot,
1004 sum_Q3l(qoi, 0) / Nlq_pilot, sum_Q4l(qoi, 0) / Nlq_pilot,
1005 cm1l, cm2l, cm3l, cm4l, Nlq_pilot);
1006
1007 cm2l_sq = cm2l * cm2l;
1008 var_of_var = (Nlq - 1.) / (Nlq * Nlq - 2. * Nlq + 3.) * (cm4l - (Nlq - 3.) / (Nlq - 1.) * cm2l_sq);
1009
1010 if(compute_gradient) {
1011 grad_g = ((Nlq * Nlq - 2. * Nlq + 3.) - (Nlq - 1.) * (2. * Nlq - 2.)) /
1012 ((Nlq * Nlq - 2. * Nlq + 3.) * (Nlq * Nlq - 2. * Nlq + 3.)) * cm4l
1013 - ((Nlq * Nlq - 2. * Nlq + 3.) - (Nlq - 3.) * (2. * Nlq - 2.)) /
1014 ((Nlq * Nlq - 2. * Nlq + 3.) * (Nlq * Nlq - 2. * Nlq + 3.)) * cm2l_sq;
1015 }
1016 //[fm] bias correction for var_P2l
1017 return var_of_var;
1018 }
1019
1020
var_of_var_ml_lmax(IntRealMatrixMap sum_Ql,IntRealMatrixMap sum_Qlm1,IntIntPairRealMatrixMap sum_QlQlm1,const size_t & Nlq_pilot,const Real & Nlq,const size_t & qoi,const bool & compute_gradient,Real & grad_g)1021 inline Real NonDMultilevelSampling::var_of_var_ml_lmax(IntRealMatrixMap sum_Ql, IntRealMatrixMap sum_Qlm1, IntIntPairRealMatrixMap sum_QlQlm1,
1022 const size_t& Nlq_pilot, const Real& Nlq, const size_t& qoi, const bool& compute_gradient, Real& grad_g){
1023 Real cm1l, cm2l, cm3l, cm4l, cm2l_sq, var_of_var;
1024
1025 IntIntPair pr11(1, 1), pr12(1, 2), pr21(2, 1), pr22(2, 2);
1026 RealMatrix &sum_Q1l = sum_Ql[1], &sum_Q1lm1 = sum_Qlm1[1],
1027 &sum_Q2l = sum_Ql[2], &sum_Q2lm1 = sum_Qlm1[2],
1028 &sum_Q3l = sum_Ql[3], &sum_Q3lm1 = sum_Qlm1[3],
1029 &sum_Q4l = sum_Ql[4], &sum_Q4lm1 = sum_Qlm1[4],
1030 &sum_Q1lQ1lm1 = sum_QlQlm1[pr11], &sum_Q1lQ2lm1 = sum_QlQlm1[pr12],
1031 &sum_Q2lQ1lm1 = sum_QlQlm1[pr21], &sum_Q2lQ2lm1 = sum_QlQlm1[pr22];
1032
1033 uncentered_to_centered(sum_Q1l(qoi, 1) / Nlq_pilot, sum_Q2l(qoi, 1) / Nlq_pilot,
1034 sum_Q3l(qoi, 1) / Nlq_pilot, sum_Q4l(qoi, 1) / Nlq_pilot,
1035 cm1l, cm2l, cm3l, cm4l, Nlq_pilot);
1036
1037 cm2l_sq = cm2l * cm2l;
1038 var_of_var = (Nlq - 1.) / (Nlq * Nlq - 2. * Nlq + 3.) * (cm4l - (Nlq - 3.) / (Nlq - 1.) * cm2l_sq);
1039
1040 if(compute_gradient) {
1041 grad_g = ((Nlq * Nlq - 2. * Nlq + 3.) - (Nlq - 1.) * (2. * Nlq - 2.)) /
1042 ((Nlq * Nlq - 2. * Nlq + 3.) * (Nlq * Nlq - 2. * Nlq + 3.)) * cm4l
1043 - ((Nlq * Nlq - 2. * Nlq + 3.) - (Nlq - 3.) * (2. * Nlq - 2.)) /
1044 ((Nlq * Nlq - 2. * Nlq + 3.) * (Nlq * Nlq - 2. * Nlq + 3.)) * cm2l_sq;
1045 }
1046 //[fm] bias correction for var_P2l
1047 return var_of_var;
1048 }
1049
var_of_var_ml_l(IntRealMatrixMap sum_Ql,IntRealMatrixMap sum_Qlm1,IntIntPairRealMatrixMap sum_QlQlm1,const size_t & Nlq_pilot,const Real & Nlq,const size_t & qoi,const size_t & lev,const bool & compute_gradient,Real & grad_g)1050 inline Real NonDMultilevelSampling::var_of_var_ml_l(IntRealMatrixMap sum_Ql, IntRealMatrixMap sum_Qlm1, IntIntPairRealMatrixMap sum_QlQlm1,
1051 const size_t& Nlq_pilot, const Real& Nlq, const size_t& qoi, const size_t& lev, const bool& compute_gradient, Real& grad_g){
1052 Real cm1l, cm2l, cm3l, cm4l, cm1lm1, cm2lm1,
1053 cm3lm1, cm4lm1, cm1l_sq, cm2l_sq, cm2lm1_sq,
1054 mu_Q2l, mu_Q2lm1, mu_Q2lQ2lm1,
1055 mu_Q1lm1_mu_Q2lQ1lm1, mu_Q1lm1_mu_Q1lm1_muQ2l, mu_Q1l_mu_Q1lQ2lm1, mu_Q1l_mu_Q1l_mu_Q2lm1,
1056 mu_Q1l_mu_Qlm1_mu_Q1lQ1lm1, mu_Q1l_mu_Q1l_mu_Q1lm1_muQ1lm1, mu_Q2l_muQ2lm1, mu_Q1lQ1lm1_mu_Q1lQ1lm1,
1057 mu_P2lP2lm1, var_P2l, var_P2lm1, covar_P2lP2lm1, term, var_of_var;
1058
1059 IntIntPair pr11(1, 1), pr12(1, 2), pr21(2, 1), pr22(2, 2);
1060 const RealMatrix &sum_Q1l = sum_Ql[1], &sum_Q1lm1 = sum_Qlm1[1],
1061 &sum_Q2l = sum_Ql[2], &sum_Q2lm1 = sum_Qlm1[2],
1062 &sum_Q3l = sum_Ql[3], &sum_Q3lm1 = sum_Qlm1[3],
1063 &sum_Q4l = sum_Ql[4], &sum_Q4lm1 = sum_Qlm1[4],
1064 &sum_Q1lQ1lm1 = sum_QlQlm1[pr11], &sum_Q1lQ2lm1 = sum_QlQlm1[pr12],
1065 &sum_Q2lQ1lm1 = sum_QlQlm1[pr21], &sum_Q2lQ2lm1 = sum_QlQlm1[pr22];
1066
1067 mu_Q2l = sum_Q2l(qoi, lev) / Nlq_pilot;
1068 uncentered_to_centered(sum_Q1l(qoi, lev) / Nlq_pilot, mu_Q2l,
1069 sum_Q3l(qoi, lev) / Nlq_pilot, sum_Q4l(qoi, lev) / Nlq_pilot,
1070 cm1l, cm2l, cm3l, cm4l, Nlq_pilot);
1071 mu_Q2lm1 = sum_Q2lm1(qoi, lev) / Nlq_pilot;
1072 uncentered_to_centered(sum_Q1lm1(qoi, lev) / Nlq_pilot, mu_Q2lm1,
1073 sum_Q3lm1(qoi, lev) / Nlq_pilot, sum_Q4lm1(qoi, lev) / Nlq_pilot,
1074 cm1lm1, cm2lm1, cm3lm1, cm4lm1, Nlq_pilot);
1075 cm2l_sq = cm2l * cm2l;
1076 cm2lm1_sq = cm2lm1 * cm2lm1;
1077
1078 // [fm] bias correction for var_P2l and var_P2lm1
1079 var_P2l = Nlq * (Nlq - 1.) / (Nlq * Nlq - 2. * Nlq + 3.) * (cm4l - (Nlq - 3.) / (Nlq - 1.) * cm2l_sq);
1080 var_P2lm1 =
1081 Nlq * (Nlq - 1.) / (Nlq * Nlq - 2. * Nlq + 3.) * (cm4lm1 - (Nlq - 3.) / (Nlq - 1.) * cm2lm1_sq);
1082
1083 //[fm] unbiased products of mean
1084 mu_Q2lQ2lm1 = sum_Q2lQ2lm1(qoi, lev) / Nlq_pilot;
1085 mu_Q1lm1_mu_Q2lQ1lm1 = unbiased_mean_product_pair(sum_Q1lm1(qoi, lev), sum_Q2lQ1lm1(qoi, lev),
1086 sum_Q2lQ2lm1(qoi, lev), Nlq_pilot);
1087 mu_Q1lm1_mu_Q1lm1_muQ2l = unbiased_mean_product_triplet(sum_Q1lm1(qoi, lev), sum_Q1lm1(qoi, lev),
1088 sum_Q2l(qoi, lev),
1089 sum_Q2lm1(qoi, lev), sum_Q2lQ1lm1(qoi, lev),
1090 sum_Q2lQ1lm1(qoi, lev),
1091 sum_Q2lQ2lm1(qoi, lev), Nlq_pilot);
1092 mu_Q1l_mu_Q1lQ2lm1 = unbiased_mean_product_pair(sum_Q1l(qoi, lev), sum_Q1lQ2lm1(qoi, lev), sum_Q2lQ2lm1(qoi, lev),
1093 Nlq_pilot);
1094 mu_Q1l_mu_Q1l_mu_Q2lm1 = unbiased_mean_product_triplet(sum_Q1l(qoi, lev), sum_Q1l(qoi, lev), sum_Q2lm1(qoi, lev),
1095 sum_Q2l(qoi, lev), sum_Q1lQ2lm1(qoi, lev),
1096 sum_Q1lQ2lm1(qoi, lev),
1097 sum_Q2lQ2lm1(qoi, lev), Nlq_pilot);
1098 mu_Q1l_mu_Qlm1_mu_Q1lQ1lm1 = unbiased_mean_product_triplet(sum_Q1l(qoi, lev), sum_Q1lm1(qoi, lev),
1099 sum_Q1lQ1lm1(qoi, lev),
1100 sum_Q1lQ1lm1(qoi, lev), sum_Q2lQ1lm1(qoi, lev),
1101 sum_Q1lQ2lm1(qoi, lev),
1102 sum_Q2lQ2lm1(qoi, lev), Nlq_pilot);
1103 mu_Q1l_mu_Q1l_mu_Q1lm1_muQ1lm1 = unbiased_mean_product_pairpair(sum_Q1l(qoi, lev), sum_Q1lm1(qoi, lev),
1104 sum_Q1lQ1lm1(qoi, lev),
1105 sum_Q2l(qoi, lev), sum_Q2lm1(qoi, lev),
1106 sum_Q2lQ1lm1(qoi, lev), sum_Q1lQ2lm1(qoi, lev),
1107 sum_Q2lQ2lm1(qoi, lev), Nlq_pilot);
1108 mu_Q2l_muQ2lm1 = unbiased_mean_product_pair(sum_Q2l(qoi, lev), sum_Q2lm1(qoi, lev), sum_Q2lQ2lm1(qoi, lev), Nlq_pilot);
1109 mu_P2lP2lm1 = mu_Q2lQ2lm1 //E[QL2 Ql2]
1110 - 2. * mu_Q1lm1_mu_Q2lQ1lm1 //E[Ql] E[QL2Ql]
1111 + 2. * mu_Q1lm1_mu_Q1lm1_muQ2l //E[Ql]2 E[QL2]
1112 - 2. * mu_Q1l_mu_Q1lQ2lm1 //E[QL] E[QLQl2]
1113 + 2. * mu_Q1l_mu_Q1l_mu_Q2lm1 //E[QL]2 E[Ql2]
1114 + 4. * mu_Q1l_mu_Qlm1_mu_Q1lQ1lm1 //E[QL] E[Ql] E[QLQl]
1115 - 4. * mu_Q1l_mu_Q1l_mu_Q1lm1_muQ1lm1 //E[QL]2 E[Ql]2
1116 - mu_Q2l_muQ2lm1; //E[QL2] E[Ql2]
1117
1118 // [fm] unbiased by opening up the square and compute three different term
1119 mu_Q1lQ1lm1_mu_Q1lQ1lm1 = unbiased_mean_product_pair(sum_Q1lQ1lm1(qoi, lev), sum_Q1lQ1lm1(qoi, lev),
1120 sum_Q2lQ2lm1(qoi, lev), Nlq_pilot);
1121 term = mu_Q1lQ1lm1_mu_Q1lQ1lm1 - 2. * mu_Q1l_mu_Qlm1_mu_Q1lQ1lm1 + mu_Q1l_mu_Q1l_mu_Q1lm1_muQ1lm1;
1122
1123 //[fm] Using only unbiased estimators the sum is also unbiased
1124 covar_P2lP2lm1
1125 = mu_P2lP2lm1 + term / (Nlq - 1.);
1126
1127 var_of_var = (var_P2l + var_P2lm1 - 2. * covar_P2lP2lm1) / Nlq;
1128
1129 if(compute_gradient) {
1130 grad_g = ((Nlq * Nlq - 2. * Nlq + 3.) - (Nlq - 1.) * (2. * Nlq - 2.)) /
1131 ((Nlq * Nlq - 2. * Nlq + 3.) * (Nlq * Nlq - 2. * Nlq + 3.)) * cm4l
1132 - ((Nlq * Nlq - 2. * Nlq + 3.) - (Nlq - 3.) * (2. * Nlq - 2.)) /
1133 ((Nlq * Nlq - 2. * Nlq + 3.) * (Nlq * Nlq - 2. * Nlq + 3.)) * cm2l_sq
1134 + ((Nlq * Nlq - 2. * Nlq + 3.) - (Nlq - 1.) * (2. * Nlq - 2.)) /
1135 ((Nlq * Nlq - 2. * Nlq + 3.) * (Nlq * Nlq - 2. * Nlq + 3.)) * cm4lm1
1136 - ((Nlq * Nlq - 2. * Nlq + 3.) - (Nlq - 3.) * (2. * Nlq - 2.)) /
1137 ((Nlq * Nlq - 2. * Nlq + 3.) * (Nlq * Nlq - 2. * Nlq + 3.)) * cm2lm1_sq
1138 - 2. * (-1. / (Nlq * Nlq) * mu_P2lP2lm1 +
1139 (-2. * Nlq + 1.) / ((Nlq * Nlq - Nlq) * (Nlq * Nlq - Nlq)) * term);
1140 }
1141
1142 return var_of_var;
1143 }
1144
unbiased_mean_product_pair(const Real & sumQ1,const Real & sumQ2,const Real & sumQ1Q2,const size_t & Nlq)1145 inline Real NonDMultilevelSampling::unbiased_mean_product_pair(const Real& sumQ1, const Real& sumQ2, const Real& sumQ1Q2, const size_t& Nlq)
1146 {
1147 Real mean1, mean2, bessel_corr1, bessel_corr2 = 0.;
1148
1149 mean1 = 1./Nlq * 1./Nlq * sumQ1 * sumQ2;
1150 mean2 = 1./Nlq * sumQ1Q2;
1151 bessel_corr1 = (Real)Nlq / ((Real)Nlq - 1.);
1152 bessel_corr2 = 1. / ((Real)Nlq - 1.);
1153
1154 return bessel_corr1*mean1 - bessel_corr2*mean2;
1155 }
1156
unbiased_mean_product_triplet(const Real & sumQ1,const Real & sumQ2,const Real & sumQ3,const Real & sumQ1Q2,const Real & sumQ1Q3,const Real & sumQ2Q3,const Real & sumQ1Q2Q3,const size_t & Nlq)1157 inline Real NonDMultilevelSampling::unbiased_mean_product_triplet(const Real& sumQ1, const Real& sumQ2, const Real& sumQ3,
1158 const Real& sumQ1Q2, const Real& sumQ1Q3, const Real& sumQ2Q3, const Real& sumQ1Q2Q3, const size_t& Nlq)
1159 {
1160 Real mean1, mean2, mean3, bessel_corr1, bessel_corr2, bessel_corr3 = 0.;
1161
1162 mean1 = 1./Nlq * 1./Nlq * 1./Nlq * sumQ1 * sumQ2 * sumQ3;
1163 mean2 = unbiased_mean_product_pair(sumQ1Q2, sumQ3, sumQ1Q2Q3, Nlq);
1164 mean2 += unbiased_mean_product_pair(sumQ2Q3, sumQ1, sumQ1Q2Q3, Nlq);
1165 mean2 += unbiased_mean_product_pair(sumQ1Q3, sumQ2, sumQ1Q2Q3, Nlq);
1166 mean3 = 1./((Real)Nlq) * sumQ1Q2Q3;
1167 bessel_corr1 = (Nlq * Nlq)/((Nlq - 1.)*(Nlq - 2.));
1168 bessel_corr2 = 1./(Nlq - 2.);
1169 bessel_corr3 = 1./((Nlq - 1.)*(Nlq - 2.));
1170
1171 return bessel_corr1 * mean1 - bessel_corr2 * mean2 - bessel_corr3 * mean3;
1172 }
1173
unbiased_mean_product_pairpair(const Real & sumQ1,const Real & sumQ2,const Real & sumQ1Q2,const Real & sumQ1sq,const Real & sumQ2sq,const Real & sumQ1sqQ2,const Real & sumQ1Q2sq,const Real & sumQ1sqQ2sq,const size_t & Nlq)1174 inline Real NonDMultilevelSampling::unbiased_mean_product_pairpair(const Real& sumQ1, const Real& sumQ2, const Real& sumQ1Q2,
1175 const Real& sumQ1sq, const Real& sumQ2sq,
1176 const Real& sumQ1sqQ2, const Real& sumQ1Q2sq, const Real& sumQ1sqQ2sq, const size_t& Nlq)
1177 {
1178 Real mean1, mean2, mean3, mean4, bessel_corr1, bessel_corr2, bessel_corr3, bessel_corr4 = 0.;
1179
1180 mean1 = 1./Nlq * 1./Nlq * 1./Nlq * 1./Nlq * sumQ1 * sumQ1 * sumQ2 * sumQ2;
1181
1182 mean2 = unbiased_mean_product_triplet(sumQ1sq, sumQ2, sumQ2, sumQ1sqQ2, sumQ1sqQ2, sumQ2sq, sumQ1sqQ2sq, Nlq);
1183 mean2 += 4. * unbiased_mean_product_triplet(sumQ1Q2, sumQ1, sumQ2, sumQ1sqQ2, sumQ1Q2sq, sumQ1Q2, sumQ1sqQ2sq, Nlq);
1184 mean2 += unbiased_mean_product_triplet(sumQ1, sumQ1, sumQ2sq, sumQ1sq, sumQ1Q2sq, sumQ1Q2sq, sumQ1sqQ2sq, Nlq);
1185
1186 mean3 = unbiased_mean_product_pair(sumQ1sq, sumQ2sq, sumQ1sqQ2sq, Nlq);
1187 mean3 += 2. * unbiased_mean_product_pair(sumQ1Q2, sumQ1Q2, sumQ1sqQ2sq, Nlq);
1188 mean3 += 2. * unbiased_mean_product_pair(sumQ1sqQ2, sumQ2, sumQ1sqQ2sq, Nlq);
1189 mean3 += 2. * unbiased_mean_product_pair(sumQ1, sumQ1Q2sq, sumQ1sqQ2sq, Nlq);
1190
1191 mean4 = 1./Nlq * sumQ1sqQ2sq;
1192
1193 bessel_corr1 = (Nlq * Nlq * Nlq)/((Nlq - 1.)*(Nlq - 2.)*(Nlq - 3.));
1194 bessel_corr2 = 1./(Nlq - 3.);
1195 bessel_corr3 = 1./((Nlq - 2.)*(Nlq - 3.));
1196 bessel_corr4 = 1./((Nlq - 1.)*(Nlq - 2.)*(Nlq - 3.));
1197
1198 return bessel_corr1 * mean1 - bessel_corr2 * mean2 - bessel_corr3 * mean3 - bessel_corr4 * mean4;
1199 }
1200
1201 } // namespace Dakota
1202
1203 #endif
1204