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: NonD
10 //- Description: Base class for NonDeterministic methods
11 //- Owner: Mike Eldred
12 //- Checked by:
13 //- Version:
14
15 #include "DakotaResponse.hpp"
16 #include "DakotaNonD.hpp"
17 #include "NonDLHSSampling.hpp"
18 #include "ProblemDescDB.hpp"
19 #include "dakota_tabular_io.hpp"
20 #include "NormalRandomVariable.hpp"
21 #include "ParallelLibrary.hpp"
22 #include "dakota_stat_util.hpp"
23
24 //#define DEBUG
25
26 static const char rcsId[]="@(#) $Id: DakotaNonD.cpp 7024 2010-10-16 01:24:42Z mseldre $";
27
28 namespace Dakota {
29
30 // initialization of statics
31 NonD* NonD::nondInstance(NULL);
32
33
NonD(ProblemDescDB & problem_db,Model & model)34 NonD::NonD(ProblemDescDB& problem_db, Model& model):
35 Analyzer(problem_db, model),
36 respLevelTarget(problem_db.get_short("method.nond.response_level_target")),
37 respLevelTargetReduce(
38 problem_db.get_short("method.nond.response_level_target_reduce")),
39 requestedRespLevels(problem_db.get_rva("method.nond.response_levels")),
40 requestedProbLevels(problem_db.get_rva("method.nond.probability_levels")),
41 requestedRelLevels(problem_db.get_rva("method.nond.reliability_levels")),
42 requestedGenRelLevels(
43 problem_db.get_rva("method.nond.gen_reliability_levels")),
44 totalLevelRequests(0),
45 cdfFlag(problem_db.get_short("method.nond.distribution") != COMPLEMENTARY),
46 pdfOutput(false),
47 finalMomentsType(problem_db.get_short("method.nond.final_moments"))
48 {
49 initialize_counts();
50
51 // When specifying z/p/beta/beta* levels, a spec with an index key (number of
52 // levels = list of ints) will result in multiple vectors of levels, one for
53 // each response fn. A spec without an index key will result in a single
54 // vector of levels. This is as much logic as the parser can support since it
55 // cannot access the number of response fns. Here, support the shorthand spec
56 // where there are multiple response fns but only one vector of levels (no
57 // index key provided), which are to be evenly distributed among the response
58 // fns. This provides some measure of backwards compatibility.
59 distribute_levels(requestedRespLevels); // always ascending
60 distribute_levels(requestedProbLevels, cdfFlag); // ascending if cumulative
61 distribute_levels(requestedRelLevels, !cdfFlag);// descending if cumulative
62 distribute_levels(requestedGenRelLevels, !cdfFlag);// descending if cumulative
63
64 for (size_t i=0; i<numFunctions; i++)
65 totalLevelRequests += requestedRespLevels[i].length() +
66 requestedProbLevels[i].length() + requestedRelLevels[i].length() +
67 requestedGenRelLevels[i].length();
68
69 // simple PDF output control for now (suppressed if quiet/silent output level)
70 if (totalLevelRequests && outputLevel >= NORMAL_OUTPUT)
71 pdfOutput = true;
72 }
73
74
NonD(unsigned short method_name,Model & model)75 NonD::NonD(unsigned short method_name, Model& model):
76 Analyzer(method_name, model), totalLevelRequests(0), cdfFlag(true),
77 pdfOutput(false), finalMomentsType(STANDARD_MOMENTS)
78 {
79 // NonDEvidence and NonDAdaptImpSampling use this ctor
80
81 initialize_counts();
82
83 // current set of statistics is mean, standard deviation, and
84 // probability of failure for each response function
85 //ShortArray asv(3*numFunctions, 1);
86 //finalStatistics = Response(numUncertainVars, asv);
87 }
88
89
NonD(unsigned short method_name,const RealVector & lower_bnds,const RealVector & upper_bnds)90 NonD::NonD(unsigned short method_name, const RealVector& lower_bnds,
91 const RealVector& upper_bnds):
92 Analyzer(method_name), epistemicStats(false), totalLevelRequests(0),
93 cdfFlag(true), pdfOutput(false), finalMomentsType(STANDARD_MOMENTS)
94 {
95 // ConcurrentStrategy uses this ctor for design opt, either for multi-start
96 // initial points or multibjective weight sets.
97
98 startCAUV = 0;
99 numCAUV = numContinuousVars
100 = std::min(lower_bnds.length(), upper_bnds.length());
101 numDiscreteIntVars = numDiscreteStringVars = numDiscreteRealVars = 0;
102 }
103
104
initialize_counts()105 void NonD::initialize_counts()
106 {
107 const Variables& vars = iteratedModel.current_variables();
108 //short active_view = vars.view().first;
109 const SizetArray& ac_totals = vars.shared_data().active_components_totals();
110 // convenience looping bounds
111 startCAUV = ac_totals[TOTAL_CDV]; numCAUV = ac_totals[TOTAL_CAUV];
112 // stats type
113 epistemicStats = (ac_totals[TOTAL_CEUV] || ac_totals[TOTAL_DEUIV] ||
114 ac_totals[TOTAL_DEUSV] || ac_totals[TOTAL_DEURV]);
115 }
116
117
resize()118 bool NonD::resize()
119 {
120 bool parent_reinit_comms = Analyzer::resize();
121
122 initialize_counts(); // sufficient at this time
123
124 return parent_reinit_comms;
125 }
126
127
derived_set_communicators(ParLevLIter pl_iter)128 void NonD::derived_set_communicators(ParLevLIter pl_iter)
129 {
130 miPLIndex = methodPCIter->mi_parallel_level_index(pl_iter);
131 iteratedModel.set_communicators(pl_iter, maxEvalConcurrency);
132 }
133
134
135 void NonD::
requested_levels(const RealVectorArray & req_resp_levels,const RealVectorArray & req_prob_levels,const RealVectorArray & req_rel_levels,const RealVectorArray & req_gen_rel_levels,short resp_lev_tgt,short resp_lev_tgt_reduce,bool cdf_flag,bool pdf_output)136 requested_levels(const RealVectorArray& req_resp_levels,
137 const RealVectorArray& req_prob_levels,
138 const RealVectorArray& req_rel_levels,
139 const RealVectorArray& req_gen_rel_levels,
140 short resp_lev_tgt, short resp_lev_tgt_reduce,
141 bool cdf_flag, bool pdf_output)
142 {
143 respLevelTarget = resp_lev_tgt;
144 respLevelTargetReduce = resp_lev_tgt_reduce;
145 cdfFlag = cdf_flag;
146
147 size_t i;
148 totalLevelRequests = 0;
149 if (req_resp_levels.empty())
150 requestedRespLevels.resize(numFunctions); // array of empty vectors
151 else {
152 requestedRespLevels = req_resp_levels;
153 // In current usage, incoming levels are already distributed.
154 //distribute_levels(requestedRespLevels);
155 for (i=0; i<numFunctions; ++i)
156 totalLevelRequests += requestedRespLevels[i].length();
157 }
158 if (req_prob_levels.empty())
159 requestedProbLevels.resize(numFunctions); // array of empty vectors
160 else {
161 requestedProbLevels = req_prob_levels;
162 // In current usage, incoming levels are already distributed.
163 //distribute_levels(requestedProbLevels, cdfFlag);
164 for (i=0; i<numFunctions; ++i)
165 totalLevelRequests += requestedProbLevels[i].length();
166 }
167 if (req_rel_levels.empty())
168 requestedRelLevels.resize(numFunctions); // array of empty vectors
169 else {
170 requestedRelLevels = req_rel_levels;
171 // In current usage, incoming levels are already distributed.
172 //distribute_levels(requestedRelLevels, !cdfFlag);
173 for (i=0; i<numFunctions; ++i)
174 totalLevelRequests += requestedRelLevels[i].length();
175 }
176 if (req_gen_rel_levels.empty())
177 requestedGenRelLevels.resize(numFunctions); // array of empty vectors
178 else {
179 requestedGenRelLevels = req_gen_rel_levels;
180 // In current usage, incoming levels are already distributed.
181 //distribute_levels(requestedGenRelLevels, !cdfFlag);
182 for (i=0; i<numFunctions; ++i)
183 totalLevelRequests += requestedGenRelLevels[i].length();
184 }
185
186 if (totalLevelRequests && pdf_output)
187 pdfOutput = true;
188
189 initialize_final_statistics();
190 initialize_response_covariance();
191 }
192
193
distribute_levels(RealVectorArray & levels,bool ascending)194 void NonD::distribute_levels(RealVectorArray& levels, bool ascending)
195 {
196 size_t i, j, num_level_arrays = levels.size();
197 // when not already performed by the parser (i.e., a num_levels key
198 // was not provided), distribute levels among an array of vectors
199 if (num_level_arrays != numFunctions) {
200 if (num_level_arrays == 0) // create array of empty vectors
201 levels.resize(numFunctions);
202 // NOTE: old response_thresholds default was a single 0 per resp fn. New
203 // default is no output level calculations if no requested input levels.
204 else if (num_level_arrays == 1) { // evenly distribute among all resp fns
205 RealVector level_array0 = levels[0];
206 size_t total_len = level_array0.length();
207 // check for divisibility
208 if (total_len%numFunctions) {
209 Cerr << "\nError: number of levels not evenly divisible by the number "
210 << "of response functions." << std::endl;
211 abort_handler(-1);
212 }
213 size_t new_len = total_len/numFunctions;
214 levels.resize(numFunctions);
215 for (i=0; i<numFunctions; ++i) {
216 RealVector& vec_i = levels[i];
217 vec_i.sizeUninitialized(new_len);
218 for (j=0; j<new_len; ++j)
219 vec_i[j] = level_array0[i*new_len+j];
220 }
221 }
222 else {
223 Cerr << "\nError: num_levels specification differs from the number of "
224 << "response functions." << std::endl;
225 abort_handler(-1);
226 }
227 }
228 // now ensure that each vector is sorted (ascending order is default)
229 for (i=0; i<numFunctions; i++) {
230 size_t len_i = levels[i].length();
231 if (len_i > 1) {
232 Real* vec_i = levels[i].values();
233 if (ascending)
234 std::sort(vec_i, vec_i + len_i); // ascending sort by default
235 else
236 std::sort(vec_i, vec_i + len_i, std::greater<Real>());// descending sort
237 }
238 }
239 }
240
241
242 void NonD::
construct_lhs(Iterator & u_space_sampler,Model & u_model,unsigned short sample_type,int num_samples,int seed,const String & rng,bool vary_pattern,short sampling_vars_mode)243 construct_lhs(Iterator& u_space_sampler, Model& u_model,
244 unsigned short sample_type, int num_samples, int seed,
245 const String& rng, bool vary_pattern, short sampling_vars_mode)
246 {
247 // sanity checks
248 if (num_samples <= 0) {
249 Cerr << "Error: bad samples specification (" << num_samples << ") in "
250 << "NonD::construct_lhs()." << std::endl;
251 abort_handler(-1);
252 }
253
254 // construct NonDLHSSampling with default sampling_vars_mode (ACTIVE)
255 u_space_sampler.assign_rep(std::make_shared<NonDLHSSampling>(u_model,
256 sample_type, num_samples, seed, rng, vary_pattern, sampling_vars_mode));
257 }
258
259
initialize_response_covariance()260 void NonD::initialize_response_covariance()
261 { } // default is no-op
262
263
264 /** Default definition of virtual function (used by sampling, reliability,
265 and stochastic expansion methods) defines the set of statistical results
266 to include the first two moments and level mappings for each QoI. */
initialize_final_statistics()267 void NonD::initialize_final_statistics()
268 {
269 size_t i, j, num_levels, cntr = 0, rl_len = 0, num_final_stats,
270 num_active_vars = iteratedModel.cv();
271 if (epistemicStats)
272 num_final_stats = 2*numFunctions;
273 else { // aleatory UQ
274 num_final_stats = (finalMomentsType) ? 2*numFunctions : 0;
275 num_final_stats += totalLevelRequests;
276 if (respLevelTargetReduce) {
277 rl_len = requestedRespLevels[0].length();
278 for (i=1; i<numFunctions; ++i)
279 if (requestedRespLevels[i].length() != rl_len) {
280 Cerr << "Error: system metric aggregation from component metrics "
281 << "requires\n consistency in length of response_levels "
282 << "among response function set." << std::endl;
283 abort_handler(-1);
284 }
285 num_final_stats += rl_len; // 1 aggregated system metric per resp level
286 }
287 }
288 // Instantiate finalStatistics:
289 // > default response ASV/DVV may be overridden by NestedModel update
290 // in subIterator.response_results_active_set(sub_iterator_set)
291 // > inactive views are not set until NestedModel ctor defines them, and
292 // subIterator construction follows in NestedModel::derived_init_comms()
293 // --> invocation of this fn from NonD ctors should have inactive view
294 ActiveSet stats_set(num_final_stats);//, num_active_vars); // default RV = 1
295 stats_set.derivative_vector(iteratedModel.inactive_continuous_variable_ids());
296 finalStatistics = Response(SIMULATION_RESPONSE, stats_set);
297
298 // Assign meaningful labels to finalStatistics (appear in NestedModel output)
299 char resp_tag[10];
300 StringArray stats_labels(num_final_stats);
301 if (epistemicStats) { // epistemic & mixed aleatory/epistemic
302 for (i=0; i<numFunctions; ++i) {
303 std::sprintf(resp_tag, "_r%zu", i+1);
304 stats_labels[cntr++] = String("z_lo") + String(resp_tag);
305 stats_labels[cntr++] = String("z_up") + String(resp_tag);
306 }
307 }
308 else { // aleatory
309 char lev_tag[10];
310 String dist_tag = (cdfFlag) ? String("cdf") : String("ccdf");
311 for (i=0; i<numFunctions; ++i) {
312 std::sprintf(resp_tag, "_r%zu", i+1);
313 if (finalMomentsType) {
314 stats_labels[cntr++] = String("mean") + String(resp_tag);
315 stats_labels[cntr++] = (finalMomentsType == CENTRAL_MOMENTS) ?
316 String("variance") + String(resp_tag) :
317 String("std_dev") + String(resp_tag);
318 }
319 num_levels = requestedRespLevels[i].length();
320 for (j=0; j<num_levels; ++j, ++cntr) {
321 switch (respLevelTarget) {
322 case PROBABILITIES: std::sprintf(lev_tag, "_plev%zu", j+1); break;
323 case RELIABILITIES: std::sprintf(lev_tag, "_blev%zu", j+1); break;
324 case GEN_RELIABILITIES: std::sprintf(lev_tag, "_b*lev%zu", j+1); break;
325 }
326 stats_labels[cntr] = dist_tag + String(lev_tag) + String(resp_tag);
327 }
328 num_levels = requestedProbLevels[i].length() +
329 requestedRelLevels[i].length() + requestedGenRelLevels[i].length();
330 for (j=0; j<num_levels; ++j, ++cntr) {
331 std::sprintf(lev_tag, "_zlev%zu", j+1);
332 stats_labels[cntr] = dist_tag + String(lev_tag) + String(resp_tag);
333 }
334 }
335 if (respLevelTargetReduce) {
336 String sys_tag("_sys");
337 for (j=0; j<rl_len; ++j, ++cntr) {
338 switch (respLevelTarget) {
339 case PROBABILITIES: std::sprintf(lev_tag, "_plev%zu", j+1); break;
340 case RELIABILITIES: std::sprintf(lev_tag, "_blev%zu", j+1); break;
341 case GEN_RELIABILITIES: std::sprintf(lev_tag, "_b*lev%zu", j+1); break;
342 }
343 stats_labels[cntr] = dist_tag + String(lev_tag) + sys_tag;
344 }
345 }
346 }
347
348 finalStatistics.function_labels(stats_labels);
349 }
350
351
pull_level_mappings(RealVector & level_maps,size_t offset)352 void NonD::pull_level_mappings(RealVector& level_maps, size_t offset)
353 {
354 if (level_maps.length() < offset + totalLevelRequests)
355 level_maps.resize(totalLevelRequests);
356
357 size_t i, j, num_lev, cntr = offset;
358 for (i=0; i<numFunctions; ++i) {
359 num_lev = requestedRespLevels[i].length();
360 switch (respLevelTarget) {
361 case PROBABILITIES:
362 for (j=0; j<num_lev; ++j, ++cntr)
363 level_maps[cntr] = computedProbLevels[i][j];
364 break;
365 case RELIABILITIES:
366 for (j=0; j<num_lev; ++j, ++cntr)
367 level_maps[cntr] = computedRelLevels[i][j];
368 break;
369 case GEN_RELIABILITIES:
370 for (j=0; j<num_lev; ++j, ++cntr)
371 level_maps[cntr] = computedGenRelLevels[i][j];
372 break;
373 }
374 num_lev = requestedProbLevels[i].length() + requestedRelLevels[i].length()
375 + requestedGenRelLevels[i].length();
376 for (j=0; j<num_lev; ++j, ++cntr)
377 level_maps[cntr] = computedRespLevels[i][j];
378 }
379 }
380
381
push_level_mappings(const RealVector & level_maps,size_t offset)382 void NonD::push_level_mappings(const RealVector& level_maps, size_t offset)
383 {
384 if (level_maps.length() < offset + totalLevelRequests) {
385 Cerr << "Error: insufficient vector length in NonD::push_level_mappings()"
386 << std::endl;
387 abort_handler(METHOD_ERROR);
388 }
389
390 size_t i, j, num_lev, cntr = offset;
391 for (i=0; i<numFunctions; ++i) {
392 num_lev = requestedRespLevels[i].length();
393 switch (respLevelTarget) {
394 case PROBABILITIES:
395 for (j=0; j<num_lev; ++j, ++cntr)
396 computedProbLevels[i][j] = level_maps[cntr];
397 break;
398 case RELIABILITIES:
399 for (j=0; j<num_lev; ++j, ++cntr)
400 computedRelLevels[i][j] = level_maps[cntr];
401 break;
402 case GEN_RELIABILITIES:
403 for (j=0; j<num_lev; ++j, ++cntr)
404 computedGenRelLevels[i][j] = level_maps[cntr];
405 break;
406 }
407 num_lev = requestedProbLevels[i].length() + requestedRelLevels[i].length()
408 + requestedGenRelLevels[i].length();
409 for (j=0; j<num_lev; ++j, ++cntr)
410 computedRespLevels[i][j] = level_maps[cntr];
411 }
412 }
413
414
415 void NonD::
load_pilot_sample(const SizetArray & pilot_spec,SizetArray & delta_N_l)416 load_pilot_sample(const SizetArray& pilot_spec, SizetArray& delta_N_l)
417 {
418 size_t pilot_size = pilot_spec.size(), delta_size = delta_N_l.size();
419 if (delta_size == pilot_size)
420 delta_N_l = pilot_spec;
421 else if (pilot_size <= 1) {
422 size_t num_samp = (pilot_size) ? pilot_spec[0] : 100;
423 delta_N_l.assign(delta_size, num_samp);
424 }
425 else {
426 Cerr << "Error: inconsistent pilot sample size (" << pilot_size
427 << ") in load_pilot_sample(SizetArray). " << delta_size
428 << " expected." << std::endl;
429 abort_handler(METHOD_ERROR);
430 }
431
432 Cout << "\nPilot sample:\n" << delta_N_l << std::endl;
433 }
434
435
436 void NonD::
load_pilot_sample(const SizetArray & pilot_spec,const Sizet3DArray & N_l,SizetArray & delta_N_l)437 load_pilot_sample(const SizetArray& pilot_spec, const Sizet3DArray& N_l,
438 SizetArray& delta_N_l)
439 {
440 size_t num_mf = N_l.size(), pilot_size = pilot_spec.size(), delta_size;
441
442 if (num_mf > 1) { // CV only case
443 delta_size = num_mf;
444 for (size_t i=0; i<num_mf; ++i)
445 if (N_l[i].size() > 1) {
446 Cerr << "Error: multidimensional N_l not expected in 1-dimensional "
447 << "load_pilot_sample(SizetArray)" << std::endl;
448 abort_handler(METHOD_ERROR);
449 }
450 Cout << "\nMultifidelity pilot sample:\n";
451 }
452 else { // ML only case
453 delta_size = N_l[0].size();
454 Cout << "\nMultilevel pilot sample:\n";
455 }
456
457 if (delta_size == pilot_size)
458 delta_N_l = pilot_spec;
459 else if (pilot_size <= 1) {
460 size_t num_samp = (pilot_size) ? pilot_spec[0] : 100;
461 delta_N_l.assign(delta_size, num_samp);
462 }
463 else {
464 Cerr << "Error: inconsistent pilot sample size (" << pilot_size
465 << ") in load_pilot_sample(SizetArray). " << delta_size
466 << " expected." << std::endl;
467 abort_handler(METHOD_ERROR);
468 }
469
470 Cout << delta_N_l << std::endl;
471 }
472
473
474 void NonD::
load_pilot_sample(const SizetArray & pilot_spec,const Sizet3DArray & N_l,Sizet2DArray & delta_N_l)475 load_pilot_sample(const SizetArray& pilot_spec, const Sizet3DArray& N_l,
476 Sizet2DArray& delta_N_l)
477 {
478 size_t i, num_samp, pilot_size = pilot_spec.size(), num_mf = N_l.size();
479 delta_N_l.resize(num_mf);
480
481 // allow several different pilot sample specifications
482 if (pilot_size <= 1) {
483 num_samp = (pilot_size) ? pilot_spec[0] : 100;
484 for (i=0; i<num_mf; ++i)
485 delta_N_l[i].assign(N_l[i].size(), num_samp);
486 }
487 else {
488 size_t j, num_lev, num_prev_lev, num_total_lev = 0;
489 bool same_lev = true;
490
491 for (i=0; i<num_mf; ++i) {
492 // for now, only SimulationModel supports solution_levels()
493 num_lev = N_l[i].size();
494 delta_N_l[i].resize(num_lev);
495 if (i && num_lev != num_prev_lev) same_lev = false;
496 num_total_lev += num_lev; num_prev_lev = num_lev;
497 }
498
499 if (same_lev && pilot_size == num_lev)
500 for (j=0; j<num_lev; ++j) {
501 num_samp = pilot_spec[j];
502 for (i=0; i<num_mf; ++i)
503 delta_N_l[i][j] = num_samp;
504 }
505 else if (pilot_size == num_total_lev) {
506 size_t cntr = 0;
507 for (i=0; i<num_mf; ++i) {
508 SizetArray& delta_N_li = delta_N_l[i]; num_lev = delta_N_li.size();
509 for (j=0; j<num_lev; ++j, ++cntr)
510 delta_N_li[j] = pilot_spec[cntr];
511 }
512 }
513 else {
514 Cerr << "Error: inconsistent pilot sample size (" << pilot_size
515 << ") in load_pilot_sample(Sizet2DArray)." << std::endl;
516 abort_handler(METHOD_ERROR);
517 }
518 }
519
520 Cout << "\nMultilevel-multifidelity pilot sample:\n";
521 print_multilevel_evaluation_summary(Cout, delta_N_l);
522 }
523
524
resize_final_statistics_gradients()525 void NonD::resize_final_statistics_gradients()
526 {
527 if (finalStatistics.is_null()) // not all ctor chains track final stats
528 return;
529
530 const ShortArray& final_asv = finalStatistics.active_set_request_vector();
531 const SizetArray& final_dvv = finalStatistics.active_set_derivative_vector();
532 size_t i, num_final_stats = final_asv.size();
533 bool final_grad_flag = false;
534 for (i=0; i<num_final_stats; i++)
535 if (final_asv[i] & 2) // no need to distinguish moment/level mapping grads
536 { final_grad_flag = true; break; }
537 finalStatistics.reshape(num_final_stats, final_dvv.size(),
538 final_grad_flag, false); // no final Hessians
539 }
540
541
update_final_statistics()542 void NonD::update_final_statistics()
543 {
544 if (finalStatistics.is_null()) // some ctor chains do not track final stats
545 return;
546
547 // this default implementation gets overridden/augmented in derived classes
548 update_aleatory_final_statistics();
549
550 if (respLevelTargetReduce) {
551 update_system_final_statistics();
552 update_system_final_statistics_gradients();
553 }
554 }
555
556
update_aleatory_final_statistics()557 void NonD::update_aleatory_final_statistics()
558 {
559 // update finalStatistics from computed{Resp,Prob,Rel,GenRel}Levels
560 size_t i, j, cntr = 0, rl_len, pl_bl_gl_len;
561 for (i=0; i<numFunctions; ++i) {
562 // final stats from compute_moments()
563 if (finalMomentsType) {
564 if (momentStats.empty())
565 cntr += 2;
566 else {
567 const Real* mom_i = momentStats[i];
568 finalStatistics.function_value(mom_i[0], cntr++); // mean
569 finalStatistics.function_value(mom_i[1], cntr++); // stdev or var
570 }
571 }
572 // final stats from compute_level_mappings()
573 rl_len = requestedRespLevels[i].length();
574 switch (respLevelTarget) { // individual component p/beta/beta*
575 case PROBABILITIES:
576 for (j=0; j<rl_len; ++j, ++cntr)
577 finalStatistics.function_value(computedProbLevels[i][j], cntr);
578 break;
579 case RELIABILITIES:
580 for (j=0; j<rl_len; ++j, ++cntr)
581 finalStatistics.function_value(computedRelLevels[i][j], cntr);
582 break;
583 case GEN_RELIABILITIES:
584 for (j=0; j<rl_len; ++j, ++cntr)
585 finalStatistics.function_value(computedGenRelLevels[i][j], cntr);
586 break;
587 }
588 pl_bl_gl_len = requestedProbLevels[i].length()
589 + requestedRelLevels[i].length() + requestedGenRelLevels[i].length();
590 for (j=0; j<pl_bl_gl_len; ++j, ++cntr)
591 finalStatistics.function_value(computedRespLevels[i][j], cntr);
592 }
593 }
594
595
update_system_final_statistics()596 void NonD::update_system_final_statistics()
597 {
598 // same rl_len enforced for all resp fns in initialize_final_statistics()
599 size_t i, j, rl_len = requestedRespLevels[0].length(),
600 cntr = totalLevelRequests;
601 if (finalMomentsType) cntr += 2*numFunctions;
602 for (j=0; j<rl_len; ++j, ++cntr) {
603 // compute system probability
604 Real system_p = 1.;
605 switch (respLevelTargetReduce) {
606 case SYSTEM_SERIES: // system p_success = product of component p_success
607 for (i=0; i<numFunctions; ++i)
608 system_p *= (1.-computedProbLevels[i][j]);
609 system_p = 1. - system_p; // convert back to p_fail
610 break;
611 case SYSTEM_PARALLEL: // system p_fail = product of component p_fail
612 for (i=0; i<numFunctions; ++i)
613 system_p *= computedProbLevels[i][j];
614 break;
615 }
616 #ifdef DEBUG
617 Cout << "\nSystem p = " << system_p << " from component p:\n";
618 for (i=0; i<numFunctions; ++i)
619 Cout << " " << computedProbLevels[i][j] << '\n';
620 Cout << '\n';
621 #endif // DEBUG
622 // convert system probability to desired system metric
623 switch (respLevelTarget) {
624 case PROBABILITIES:
625 finalStatistics.function_value(system_p, cntr); break;
626 case RELIABILITIES: case GEN_RELIABILITIES:
627 finalStatistics.function_value(
628 -Pecos::NormalRandomVariable::inverse_std_cdf(system_p), cntr);
629 break;
630 }
631 }
632 }
633
634
update_system_final_statistics_gradients()635 void NonD::update_system_final_statistics_gradients()
636 {
637 const ShortArray& final_asv = finalStatistics.active_set_request_vector();
638 const SizetArray& final_dvv = finalStatistics.active_set_derivative_vector();
639 // same rl_len enforced for all resp fns in initialize_final_statistics()
640 size_t l, v, s, p, rl_len = requestedRespLevels[0].length(),
641 num_deriv_vars = final_dvv.size(), cntr = totalLevelRequests,
642 moment_offset = (finalMomentsType) ? 2 : 0;
643 if (finalMomentsType) cntr += 2*numFunctions;
644 RealVector final_stat_grad(num_deriv_vars, false);
645 RealVectorArray component_grad(numFunctions);
646 Real prod;
647 for (l=0; l<rl_len; ++l, ++cntr) {
648 if (final_asv[cntr] & 2) {
649 // Retrieve component probability gradients from finalStatistics
650 size_t index = 0;
651 for (s=0; s<numFunctions; ++s) {
652 index += moment_offset;
653 if (respLevelTarget == PROBABILITIES)
654 component_grad[s] = finalStatistics.function_gradient_view(index+l);
655 else {
656 component_grad[s] = finalStatistics.function_gradient_copy(index+l);
657 Real component_beta = finalStatistics.function_value(index+l);
658 component_grad[s].scale(
659 -Pecos::NormalRandomVariable::std_pdf(-component_beta));
660 }
661 index += rl_len + requestedProbLevels[s].length() +
662 requestedRelLevels[s].length() + requestedGenRelLevels[s].length();
663 #ifdef DEBUG
664 Cout << "component gradient " << s+1 << ":\n"
665 << component_grad[s] << '\n';
666 #endif // DEBUG
667 }
668 // Compute system probability
669 for (v=0; v<num_deriv_vars; ++v) {
670 // apply product rule over n factors
671 Real& sum = final_stat_grad[v]; sum = 0.;
672 for (s=0; s<numFunctions; ++s) {
673 prod = 1.;
674 switch (respLevelTargetReduce) {
675 case SYSTEM_SERIES:// system p_success = prod of component p_success
676 for (p=0; p<numFunctions; ++p)
677 prod *= (p == s) ? -component_grad[p][v] :
678 1.-computedProbLevels[p][l];
679 break;
680 case SYSTEM_PARALLEL: // system p_fail = product of component p_fail
681 for (p=0; p<numFunctions; ++p)
682 prod *= (p == s) ? component_grad[p][v] :
683 computedProbLevels[p][l];
684 break;
685 }
686 sum += prod;
687 }
688 }
689 Real factor = 1.; bool scale = false;
690 // negate gradient if converting system p_success to system p_fail
691 if (respLevelTargetReduce == SYSTEM_SERIES)
692 { factor *= -1.; scale = true; }
693 // define any scaling for system metric type
694 if (respLevelTarget != PROBABILITIES) {
695 Real sys_beta = finalStatistics.function_value(cntr);
696 factor /= -Pecos::NormalRandomVariable::std_pdf(-sys_beta);
697 scale = true;
698 }
699 if (scale) final_stat_grad.scale(factor);
700 finalStatistics.function_gradient(final_stat_grad, cntr);
701 }
702 }
703 }
704
705
initialize_level_mappings()706 void NonD::initialize_level_mappings()
707 {
708 // Default sizing assumes no distinction between requested and achieved
709 // levels for the same measure (the request is always achieved) and assumes
710 // probability (e.g., computed by binning) and reliability (e.g., computed
711 // by moment projection) are not collapsible.
712 if (computedRespLevels.empty() || computedProbLevels.empty() ||
713 computedRelLevels.empty() || computedGenRelLevels.empty()) {
714 computedRespLevels.resize(numFunctions);
715 computedProbLevels.resize(numFunctions);
716 computedRelLevels.resize(numFunctions);
717 computedGenRelLevels.resize(numFunctions);
718 for (size_t i=0; i<numFunctions; ++i) {
719 switch (respLevelTarget) {
720 case PROBABILITIES:
721 computedProbLevels[i].resize(requestedRespLevels[i].length()); break;
722 case RELIABILITIES:
723 computedRelLevels[i].resize(requestedRespLevels[i].length()); break;
724 case GEN_RELIABILITIES:
725 computedGenRelLevels[i].resize(requestedRespLevels[i].length()); break;
726 }
727 computedRespLevels[i].resize(requestedProbLevels[i].length() +
728 requestedRelLevels[i].length() +
729 requestedGenRelLevels[i].length());
730 }
731 }
732 }
733
734
735 /** This function infers PDFs from the CDF/CCDF level mappings, in order
736 to enable PDF computation after CDF/CCDF probability level refinement
737 (e.g., from importance sampling).
738
739 prob_refinement alerts the routine to exclude inverse mappings from
740 the PDF, since refinement only applies to z->p forward mappings and
741 mixing refined and unrefined probability mappings results in an
742 inconsistency (potentially manifesting as negative density values).
743
744 all_levels_computed is an option used by reliability methods where
745 computed*Levels are defined across the union of all requested levels. */
746 void NonD::
compute_densities(const RealRealPairArray & min_max_fns,bool prob_refinement,bool all_levels_computed)747 compute_densities(const RealRealPairArray& min_max_fns,
748 bool prob_refinement, bool all_levels_computed)
749 {
750 if (!pdfOutput)
751 return;
752
753 computedPDFAbscissas.resize(numFunctions);
754 computedPDFOrdinates.resize(numFunctions);
755 archive_allocate_pdf();
756
757 size_t i, j, cntr, core_pdf_bins, pdf_size, offset;
758 Real z, min, max, prev_r, prev_p, new_r, new_p, last_r;
759 std::map<Real, Real> cdf_map;
760 std::map<Real, Real>::iterator it, it_last;
761 pdfComputed.resize(numFunctions, false);
762 for (i=0; i<numFunctions; ++i) {
763
764 // CDF/CCDF mappings: z -> p/beta/beta* and p/beta/beta* -> z
765 const RealVector& req_rlev_i = requestedRespLevels[i];
766 const RealVector& req_plev_i = requestedProbLevels[i];
767 const RealVector& req_glev_i = requestedGenRelLevels[i];
768 const RealVector& comp_rlev_i = computedRespLevels[i];
769 const RealVector& comp_plev_i = computedProbLevels[i];
770 const RealVector& comp_glev_i = computedGenRelLevels[i];
771 size_t rl_len = req_rlev_i.length(), pl_len = req_plev_i.length(),
772 bl_len = requestedRelLevels[i].length(), gl_len = req_glev_i.length();
773
774 // Define a unique sorted map of response levels -> cdf probabilities.
775 // (refer to NonD::initialize_level_mappings() for default indexing;
776 cdf_map.clear(); min = min_max_fns[i].first; max = min_max_fns[i].second;
777 if (all_levels_computed) {
778 for (j=0; j<rl_len; ++j) {
779 z = comp_rlev_i[j]; // request may be outside extreme values
780 cdf_map[z] = (cdfFlag) ? comp_plev_i[j] : 1.-comp_plev_i[j];
781 }
782 if (!prob_refinement || !rl_len) { // don't combine refined/unrefined
783 size_t total_len = rl_len + pl_len + bl_len + gl_len;
784 for (j=rl_len; j<total_len; ++j) {
785 z = comp_rlev_i[j];
786 //if (z >= min && z <= max) // exclude any extrapolations outside bnds
787 cdf_map[z] = (cdfFlag) ? comp_plev_i[j] : 1.-comp_plev_i[j];
788 }
789 }
790 }
791 else {
792 switch (respLevelTarget) {
793 case PROBABILITIES:
794 for (j=0; j<rl_len; ++j) {
795 z = req_rlev_i[j]; // request may be outside extreme values
796 cdf_map[z] = (cdfFlag) ? comp_plev_i[j] : 1.-comp_plev_i[j];
797 }
798 break;
799 //case RELIABILITIES: // exclude reliability level mappings from cdf_map
800 case GEN_RELIABILITIES:
801 for (j=0; j<rl_len; ++j) {
802 z = req_rlev_i[j]; // request may be outside extreme values
803 Real g_cdf = (cdfFlag) ? comp_glev_i[j] : -comp_glev_i[j];
804 cdf_map[z] = Pecos::NormalRandomVariable::std_cdf(-g_cdf);
805 }
806 break;
807 }
808 if (!prob_refinement || !rl_len) { // don't combine refined/unrefined
809 for (j=0; j<pl_len; ++j) {
810 z = comp_rlev_i[j];
811 //if (z >= min && z <= max) // exclude any extrapolations outside bnds
812 cdf_map[z] = (cdfFlag) ? req_plev_i[j] : 1.-req_plev_i[j];
813 }
814 // exclude reliability level mappings from cdf_map
815 for (j=0, cntr=pl_len+bl_len; j<gl_len; ++j, ++cntr) {
816 z = comp_rlev_i[cntr];
817 //if (z >= min && z <= max) {//exclude any extrapolations outside bnds
818 Real g_cdf = (cdfFlag) ? req_glev_i[j] : -req_glev_i[j];
819 cdf_map[z] = Pecos::NormalRandomVariable::std_cdf(-g_cdf);
820 //}
821 }
822 }
823 }
824
825 if (!cdf_map.empty()) {
826 it = cdf_map.begin(); prev_r = it->first;
827 it_last = --cdf_map.end(); last_r = it_last->first;
828 // compute computedPDF{Abscissas,Ordinates} from bin counts and widths
829 core_pdf_bins = cdf_map.size()-1; pdf_size = core_pdf_bins;
830 if (min < prev_r) ++pdf_size;
831 if (max > last_r) ++pdf_size;
832 RealVector& abs_i = computedPDFAbscissas[i]; abs_i.resize(pdf_size+1);
833 RealVector& ord_i = computedPDFOrdinates[i]; ord_i.resize(pdf_size);
834 if (min < prev_r) { // first bin accumulates p0 over [min,lev0]
835 offset = 1; prev_p = it->second;
836 abs_i[0] = min; ord_i[0] = prev_p/(prev_r - min);
837 }
838 else { // first bin accumulates p0+p1 over [lev0,lev1]
839 offset = 0; prev_p = 0.;
840 }
841 for (j=0; j<core_pdf_bins; ++j) {
842 ++it; new_r = it->first; new_p = it->second;
843 abs_i[j+offset] = prev_r;
844 ord_i[j+offset] = (new_p - prev_p) / (new_r - prev_r);
845 prev_r = new_r; prev_p = new_p;
846 }
847 if (max > last_r) {
848 abs_i[pdf_size-1] = last_r;
849 ord_i[pdf_size-1] = (1. - it_last->second)/(max - last_r);
850 abs_i[pdf_size] = max; // no ordinate for final abscissa
851 }
852 else
853 abs_i[pdf_size] = last_r; // no ordinate for final abscissa
854 pdfComputed[i] = true;
855 }
856 }
857 }
858
859
860 /** Print distribution mappings, including to file per response. */
861 void NonD::
print_level_mappings(std::ostream & s,String qoi_type,const StringArray & qoi_labels) const862 print_level_mappings(std::ostream& s, String qoi_type,
863 const StringArray& qoi_labels) const
864 {
865 print_densities(s, qoi_type, qoi_labels);
866
867 // output CDF/CCDF probabilities resulting from binning or CDF/CCDF
868 // reliabilities resulting from number of std devs separating mean & target
869 s << std::scientific << std::setprecision(write_precision)
870 << "\nLevel mappings for each " << qoi_type << ":\n";
871 for (size_t i=0; i<numFunctions; ++i)
872 if (!requestedRespLevels[i].empty() || !requestedProbLevels[i].empty() ||
873 !requestedRelLevels[i].empty() || !requestedGenRelLevels[i].empty()) {
874 print_level_map(s, i, qoi_labels[i]);
875 // optionally write the distribution mapping to a .dist file
876 if (outputLevel >= VERBOSE_OUTPUT)
877 level_mappings_file(i, qoi_labels[i]);
878 }
879 }
880
881
882 /** This version differs in its use of a concatenated vector of level
883 mappings, rathen than computed{Resp,Prob,Real,GenRel}Levels. */
884 void NonD::
print_level_mappings(std::ostream & s,const RealVector & level_maps,bool moment_offset,const String & prepend)885 print_level_mappings(std::ostream& s, const RealVector& level_maps,
886 bool moment_offset, const String& prepend)
887 {
888 if (level_maps.empty()) return;
889
890 if (prepend.empty()) s << "\nLevel mappings for each response function:\n";
891 else s << '\n' << prepend << " level mappings for each response function:\n";
892
893 size_t i, j, cntr,
894 width = write_precision+7, w2p2 = 2*width+2, w3p4 = 3*width+4;
895 const StringArray& qoi_labels = iteratedModel.response_labels();
896 for (i=0, cntr=0; i<numFunctions; ++i) {
897 if (moment_offset) cntr += 2; // skip over moments, if present
898 if (cdfFlag) s << "Cumulative Distribution Function (CDF) for ";
899 else s << "Complementary Cumulative Distribution Function (CCDF) for ";
900 s << qoi_labels[i] << ":\n Response Level Probability Level "
901 << "Reliability Index General Rel Index\n -------------- "
902 << "----------------- ----------------- -----------------\n";
903 const RealVector& req_z_levs = requestedRespLevels[i];
904 size_t num_z_levs = req_z_levs.length();
905 for (j=0; j<num_z_levs; ++j, ++cntr) {
906 s << " " << std::setw(width) << req_z_levs[j] << " ";
907 switch (respLevelTarget) {
908 case PROBABILITIES:
909 s << std::setw(width) << level_maps[cntr] << '\n'; break;
910 case RELIABILITIES:
911 s << std::setw(w2p2) << level_maps[cntr] << '\n'; break;
912 case GEN_RELIABILITIES:
913 s << std::setw(w3p4) << level_maps[cntr] << '\n'; break;
914 }
915 }
916 const RealVector& req_p_levs = requestedProbLevels[i];
917 size_t num_p_levs = req_p_levs.length();
918 for (j=0; j<num_p_levs; ++j, ++cntr)
919 s << " " << std::setw(width) << level_maps[cntr]
920 << " " << std::setw(width) << req_p_levs[j] << '\n';
921 const RealVector& req_b_levs = requestedRelLevels[i];
922 size_t num_b_levs = req_b_levs.length();
923 for (j=0; j<num_b_levs; ++j, ++cntr)
924 s << " " << std::setw(width) << level_maps[cntr]
925 << " " << std::setw(w2p2) << req_b_levs[j] << '\n';
926 const RealVector& req_g_levs = requestedGenRelLevels[i];
927 size_t num_g_levs = req_g_levs.length();
928 for (j=0; j<num_g_levs; ++j, ++cntr)
929 s << " " << std::setw(width) << level_maps[cntr]
930 << " " << std::setw(w3p4) << req_g_levs[j] << '\n';
931 }
932 }
933
934
935 void NonD::
print_densities(std::ostream & s,String qoi_type,const StringArray & pdf_labels) const936 print_densities(std::ostream& s, String qoi_type,
937 const StringArray& pdf_labels) const
938 {
939 if (!pdfOutput)
940 return;
941
942 size_t i, j, wpp7 = write_precision+7, num_qoi = computedPDFOrdinates.size();
943 if (num_qoi)
944 s << std::scientific << std::setprecision(write_precision)
945 << "\nProbability Density Function (PDF) histograms for each "
946 << qoi_type << ":\n";
947 for (i=0; i<num_qoi; ++i) {
948 const RealVector& ord_i = computedPDFOrdinates[i];
949 const RealVector& abs_i = computedPDFAbscissas[i];
950 size_t pdf_len = ord_i.length();
951 if (pdf_len) {
952 s << "PDF for " << pdf_labels[i] << ":\n"
953 << " Bin Lower Bin Upper Density Value\n"
954 << " --------- --------- -------------\n";
955 for (j=0; j<pdf_len; ++j)
956 s << " " << std::setw(wpp7) << abs_i[j]
957 << " " << std::setw(wpp7) << abs_i[j+1]
958 << " " << std::setw(wpp7) << ord_i[j] << '\n';
959 }
960 }
961 }
962
963
964 /// Write distribution mappings to a file for a single response
level_mappings_file(size_t fn_index,const String & qoi_label) const965 void NonD::level_mappings_file(size_t fn_index, const String& qoi_label) const
966 {
967 std::string dist_filename = qoi_label + ".dist";
968 std::ofstream dist_file;
969 TabularIO::open_file(dist_file, dist_filename, "Distribution Map Output");
970 dist_file << std::scientific << std::setprecision(write_precision);
971 print_level_map(dist_file, fn_index, qoi_label);
972 }
973
974
975 /** Print the distribution mapping for a single response function to
976 the passed output stream. This base class version maps from one
977 requested level type to one computed level type; some derived
978 class implementations (e.g., local and global reliability) output
979 multiple computed level types. */
980 void NonD::
print_level_map(std::ostream & s,size_t fn_index,const String & qoi_label) const981 print_level_map(std::ostream& s, size_t fn_index, const String& qoi_label) const
982 {
983 size_t j, width = write_precision+7, w2p2 = 2*width+2, w3p4 = 3*width+4;
984
985 if (cdfFlag)
986 s << "Cumulative Distribution Function (CDF) for ";
987 else
988 s << "Complementary Cumulative Distribution Function (CCDF) for ";
989 s << qoi_label << ":\n Response Level Probability Level "
990 << "Reliability Index General Rel Index\n -------------- "
991 << "----------------- ----------------- -----------------\n";
992 size_t num_resp_levels = requestedRespLevels[fn_index].length();
993 for (j=0; j<num_resp_levels; j++) { // map from 1 requested to 1 computed
994 s << " " << std::setw(width) << requestedRespLevels[fn_index][j] << " ";
995 switch (respLevelTarget) {
996 case PROBABILITIES:
997 s << std::setw(width) << computedProbLevels[fn_index][j] << '\n'; break;
998 case RELIABILITIES:
999 s << std::setw(w2p2) << computedRelLevels[fn_index][j] << '\n'; break;
1000 case GEN_RELIABILITIES:
1001 s << std::setw(w3p4) << computedGenRelLevels[fn_index][j] << '\n'; break;
1002 }
1003 }
1004 size_t num_prob_levels = requestedProbLevels[fn_index].length();
1005 for (j=0; j<num_prob_levels; j++) // map from 1 requested to 1 computed
1006 s << " " << std::setw(width) << computedRespLevels[fn_index][j]
1007 << " " << std::setw(width) << requestedProbLevels[fn_index][j] << '\n';
1008 size_t num_rel_levels = requestedRelLevels[fn_index].length(),
1009 offset = num_prob_levels;
1010 for (j=0; j<num_rel_levels; j++) // map from 1 requested to 1 computed
1011 s << " " << std::setw(width) << computedRespLevels[fn_index][j+offset]
1012 << " " << std::setw(w2p2) << requestedRelLevels[fn_index][j] << '\n';
1013 size_t num_gen_rel_levels = requestedGenRelLevels[fn_index].length();
1014 offset += num_rel_levels;
1015 for (j=0; j<num_gen_rel_levels; j++) // map from 1 requested to 1 computed
1016 s << " " << std::setw(width) << computedRespLevels[fn_index][j+offset]
1017 << " " << std::setw(w3p4) << requestedGenRelLevels[fn_index][j] << '\n';
1018 }
1019
1020
print_system_mappings(std::ostream & s) const1021 void NonD::print_system_mappings(std::ostream& s) const
1022 {
1023 size_t rl_len = requestedRespLevels[0].length();
1024 if (respLevelTargetReduce && rl_len) {
1025 size_t i, width = write_precision+7, g_width = 2*width+2,
1026 cntr = totalLevelRequests;
1027 if (finalMomentsType) cntr += 2*numFunctions;
1028 const RealVector& final_stats = finalStatistics.function_values();
1029 s << std::scientific << std::setprecision(write_precision)
1030 << "\nSystem response level mappings:\n";
1031 if (cdfFlag) s << "Cumulative distribution metrics ";
1032 else s << "Complementary cumulative distribution metrics ";
1033 if (respLevelTargetReduce == SYSTEM_SERIES) s << "for series ";
1034 else if (respLevelTargetReduce == SYSTEM_PARALLEL) s << "for parallel ";
1035 s << "system:\n Resp Level Set Probability Level Reliability Index "
1036 << "General Rel Index\n -------------- ----------------- "
1037 << "----------------- -----------------\n";
1038 Real prob, gen_rel;
1039 for (i=0; i<rl_len; ++i, ++cntr) {
1040 switch (respLevelTarget) {
1041 case PROBABILITIES:
1042 prob = final_stats[cntr];
1043 gen_rel = -Pecos::NormalRandomVariable::inverse_std_cdf(prob); break;
1044 default:
1045 gen_rel = final_stats[cntr];
1046 prob = Pecos::NormalRandomVariable::std_cdf(-gen_rel); break;
1047 }
1048 s << " " << std::setw(width) << i+1 << " " << std::setw(width) << prob
1049 << " " << std::setw(g_width) << gen_rel << '\n';
1050 }
1051 }
1052 }
1053
1054
1055 void NonD::
print_multilevel_evaluation_summary(std::ostream & s,const SizetArray & N_samp)1056 print_multilevel_evaluation_summary(std::ostream& s, const SizetArray& N_samp)
1057 {
1058 size_t j, width = write_precision+7, num_lev = N_samp.size();
1059 for (j=0; j<num_lev; ++j)
1060 s << " " << std::setw(width) << N_samp[j] << '\n';
1061 }
1062
1063
1064 void NonD::
print_multilevel_evaluation_summary(std::ostream & s,const Sizet2DArray & N_samp)1065 print_multilevel_evaluation_summary(std::ostream& s, const Sizet2DArray& N_samp)
1066 {
1067 size_t j, width = write_precision+7, num_lev = N_samp.size();
1068 for (j=0; j<num_lev; ++j) {
1069 const SizetArray& N_j = N_samp[j];
1070 s << " " << std::setw(width) << N_j[0];
1071 if (!homogeneous(N_j)) { // print all counts in this 1D array
1072 size_t q, num_q = N_j.size();
1073 for (size_t q=1; q<num_q; ++q)
1074 s << ' ' << N_j[q];
1075 }
1076 s << '\n';
1077 }
1078 }
1079
1080
1081 void NonD::
print_multilevel_evaluation_summary(std::ostream & s,const Sizet3DArray & N_samp)1082 print_multilevel_evaluation_summary(std::ostream& s, const Sizet3DArray& N_samp)
1083 {
1084 size_t i, j, num_mf = N_samp.size(), width = write_precision+7;
1085 if (num_mf == 1) s << "<<<<< Final samples per level:\n";
1086 else s << "<<<<< Final samples per model form:\n";
1087 for (i=0; i<num_mf; ++i) {
1088 if (num_mf > 1) s << " Model Form " << i+1 << ":\n";
1089 print_multilevel_evaluation_summary(s, N_samp[i]);
1090 }
1091 }
1092
1093
archive_allocate_mappings()1094 void NonD::archive_allocate_mappings()
1095 {
1096 if (!resultsDB.active()) return;
1097
1098 bool req_resp = false, req_prob = false, req_rel = false, req_gen = false;
1099 for (size_t i=0; i<numFunctions; ++i) {
1100 if (requestedRespLevels[i].length() > 0) req_resp = true;
1101 if (requestedProbLevels[i].length() > 0) req_prob = true;
1102 if (requestedRelLevels[i].length() > 0) req_rel = true;
1103 if (requestedGenRelLevels[i].length() > 0) req_gen = true;
1104 }
1105
1106 if (req_resp) {
1107 std::string resp_target, data_name;
1108 switch (respLevelTarget) {
1109 case PROBABILITIES:
1110 resp_target = "Probability";
1111 data_name = resultsNames.map_resp_prob;
1112 break;
1113 case RELIABILITIES:
1114 resp_target = "Reliability";
1115 data_name = resultsNames.map_resp_rel;
1116 break;
1117 case GEN_RELIABILITIES:
1118 resp_target = "Generalized Reliability";
1119 data_name = resultsNames.map_resp_genrel;
1120 break;
1121 }
1122
1123 // mapping per function, possibly empty
1124 MetaDataType md;
1125 md["Array Spans"] = make_metadatavalue("Response Functions");
1126 md["Column Labels"] =
1127 make_metadatavalue("Response Level", resp_target + " Level");
1128 resultsDB.array_allocate<RealMatrix>
1129 (run_identifier(), data_name, numFunctions, md);
1130 }
1131 if (req_prob) {
1132 // mapping per function, possibly empty
1133 MetaDataType md;
1134 md["Array Spans"] = make_metadatavalue("Response Functions");
1135 md["Column Labels"] =
1136 make_metadatavalue("Probability Level", "Response Level");
1137 resultsDB.array_allocate<RealMatrix>
1138 (run_identifier(), resultsNames.map_prob_resp, numFunctions, md);
1139 }
1140 if (req_rel) {
1141 // mapping per function, possibly empty
1142 MetaDataType md;
1143 md["Array Spans"] = make_metadatavalue("Response Functions");
1144 md["Column Labels"] =
1145 make_metadatavalue("Reliability Level", "Response Level");
1146 resultsDB.array_allocate<RealMatrix>
1147 (run_identifier(), resultsNames.map_rel_resp, numFunctions, md);
1148 }
1149 if (req_gen) {
1150 // mapping per function, possibly empty
1151 MetaDataType md;
1152 md["Array Spans"] = make_metadatavalue("Response Functions");
1153 md["Column Labels"] =
1154 make_metadatavalue("Generalized Reliability Level", "Response Level");
1155 resultsDB.array_allocate<RealMatrix>
1156 (run_identifier(), resultsNames.map_genrel_resp, numFunctions, md);
1157 }
1158 }
1159
1160
1161 // archive the mappings from response levels
archive_from_resp(size_t i,size_t inc_id)1162 void NonD::archive_from_resp(size_t i, size_t inc_id)
1163 {
1164 // only insert if active and response levels specified
1165 size_t num_resp_levels = requestedRespLevels[i].length();
1166 if (!resultsDB.active() || num_resp_levels == 0) return;
1167
1168 size_t j;
1169 std::string data_name;
1170
1171 RealMatrix mapping(num_resp_levels, 2);
1172
1173 DimScaleMap scale;
1174 scale.emplace(0, RealScale("response_levels", requestedRespLevels[i]));
1175 const StringArray &labels = iteratedModel.response_labels();
1176 RealVector *result;
1177
1178 // TODO: could use SetCol?
1179 switch (respLevelTarget) {
1180 case PROBABILITIES:
1181 data_name = resultsNames.map_resp_prob;
1182 for (j=0; j<num_resp_levels; ++j) {
1183 mapping(j, 0) = requestedRespLevels[i][j];
1184 mapping(j, 1) = computedProbLevels[i][j];
1185 }
1186 result = &computedProbLevels[i];
1187 break;
1188 case RELIABILITIES:
1189 data_name = resultsNames.map_resp_rel;
1190 for (j=0; j<num_resp_levels; ++j) {
1191 mapping(j, 0) = requestedRespLevels[i][j];
1192 mapping(j, 1) = computedRelLevels[i][j];
1193 }
1194 result = &computedRelLevels[i];
1195 break;
1196 case GEN_RELIABILITIES:
1197 data_name = resultsNames.map_resp_genrel;
1198 for (j=0; j<num_resp_levels; ++j) {
1199 mapping(j, 0) = requestedRespLevels[i][j];
1200 mapping(j, 1) = computedGenRelLevels[i][j];
1201 }
1202 result = &computedGenRelLevels[i];
1203 break;
1204 }
1205
1206 StringArray location;
1207 if(inc_id) location.push_back(String("increment:") + std::to_string(inc_id));
1208 location.push_back(String("response_levels"));
1209 location.push_back(labels[i]);
1210
1211 resultsDB.insert(run_identifier(), location, *result, scale);
1212
1213 resultsDB.array_insert<RealMatrix>(run_identifier(), data_name, i, mapping);
1214 }
1215
1216
1217 // archive the mappings to response levels
archive_to_resp(size_t i,size_t inc_id)1218 void NonD::archive_to_resp(size_t i, size_t inc_id)
1219 {
1220 if (!resultsDB.active()) return;
1221
1222 DimScaleMap scale;
1223 const StringArray &labels = iteratedModel.response_labels();
1224 StringArray location;
1225 size_t r_index = 0;
1226 if(inc_id) {
1227 location.push_back(String("increment:") + std::to_string(inc_id));
1228 r_index = 1;
1229 }
1230 location.push_back(String(""));
1231 location.push_back(labels[i]);
1232 size_t j;
1233 size_t num_prob_levels = requestedProbLevels[i].length();
1234 if (num_prob_levels > 0) {
1235 RealMatrix mapping(num_prob_levels, 2);
1236 for (j=0; j<num_prob_levels; j++) {
1237 mapping(j, 0) = requestedProbLevels[i][j];
1238 mapping(j, 1) = computedRespLevels[i][j];
1239 }
1240 resultsDB.
1241 array_insert<RealMatrix>(run_identifier(),
1242 resultsNames.map_prob_resp, i, mapping);
1243 location[r_index] = String("probability_levels");
1244 scale.emplace(0, RealScale("probability_levels", requestedProbLevels[i]));
1245 RealVector result(Teuchos::View, &computedRespLevels[i][0], num_prob_levels);
1246 resultsDB.insert(run_identifier(), location, result, scale);
1247 }
1248 size_t num_rel_levels = requestedRelLevels[i].length();
1249 size_t offset = num_prob_levels;
1250 if (num_rel_levels > 0) {
1251 RealMatrix mapping(num_rel_levels, 2);
1252 for (j=0; j<num_rel_levels; j++) {
1253 mapping(j, 0) = requestedRelLevels[i][j];
1254 mapping(j, 1) = computedRespLevels[i][j+offset];
1255 }
1256 resultsDB.
1257 array_insert<RealMatrix>(run_identifier(),
1258 resultsNames.map_rel_resp, i, mapping);
1259
1260 scale.emplace(0, RealScale("reliability_levels", requestedRelLevels[i]));
1261 RealVector result(Teuchos::View, &computedRespLevels[i][0] + offset, num_rel_levels);
1262 location[r_index] = String("reliability_levels");
1263 resultsDB.insert(run_identifier(), location, result, scale);
1264 }
1265 size_t num_gen_rel_levels = requestedGenRelLevels[i].length();
1266 offset += num_rel_levels;
1267 if (num_gen_rel_levels > 0) {
1268 RealMatrix mapping(num_gen_rel_levels, 2);
1269 for (j=0; j<num_gen_rel_levels; j++) {
1270 mapping(j, 0) = requestedGenRelLevels[i][j];
1271 mapping(j, 1) = computedRespLevels[i][j+offset];
1272 }
1273 resultsDB.
1274 array_insert<RealMatrix>(run_identifier(),
1275 resultsNames.map_genrel_resp, i, mapping);
1276
1277 scale.emplace(0, RealScale("gen_reliability_levels", requestedGenRelLevels[i]));
1278 RealVector result(Teuchos::View, &computedRespLevels[i][0] + offset, num_gen_rel_levels);
1279 location[r_index] = String("gen_reliability_levels");
1280 resultsDB.insert(run_identifier(), location, result, scale);
1281
1282 }
1283 }
1284
1285
archive_allocate_pdf()1286 void NonD::archive_allocate_pdf() // const
1287 {
1288 if (!resultsDB.active()) return;
1289
1290 // pdf per function, possibly empty
1291 MetaDataType md;
1292 md["Array Spans"] = make_metadatavalue("Response Functions");
1293 md["Row Labels"] =
1294 make_metadatavalue("Bin Lower", "Bin Upper", "Density Value");
1295 resultsDB.array_allocate<RealMatrix>
1296 (run_identifier(), resultsNames.pdf_histograms, numFunctions, md);
1297 }
1298
1299
archive_pdf(size_t i,size_t inc_id)1300 void NonD::archive_pdf(size_t i, size_t inc_id) // const
1301 {
1302 if (!resultsDB.active() || !pdfOutput ) return;
1303
1304 size_t pdf_len = computedPDFOrdinates[i].length();
1305 RealMatrix pdf(3, pdf_len);
1306 for (size_t j=0; j<pdf_len; ++j) {
1307 pdf(0, j) = computedPDFAbscissas[i][j];
1308 pdf(1, j) = computedPDFAbscissas[i][j+1];
1309 pdf(2, j) = computedPDFOrdinates[i][j];
1310 }
1311
1312 resultsDB.array_insert<RealMatrix>
1313 (run_identifier(), resultsNames.pdf_histograms, i, pdf);
1314
1315 const StringArray &labels = iteratedModel.response_labels();
1316 StringArray location;
1317 if(inc_id) location.push_back(String("increment:") + std::to_string(inc_id));
1318 location.push_back("probability_density");
1319 location.push_back(labels[i]);
1320 DimScaleMap scales;
1321 scales.emplace(0, RealScale("lower_bounds", &computedPDFAbscissas[i][0], pdf_len, ScaleScope::UNSHARED));
1322 scales.emplace(0, RealScale("upper_bounds", &computedPDFAbscissas[i][1], pdf_len, ScaleScope::UNSHARED));
1323 resultsDB.insert(run_identifier(),location, computedPDFOrdinates[i], scales);
1324 }
1325
archive_equiv_hf_evals(const Real equiv_hf_evals)1326 void NonD::archive_equiv_hf_evals(const Real equiv_hf_evals) {
1327 if (!resultsDB.active()) return;
1328 resultsDB.add_metadata_to_execution(run_identifier(),
1329 {ResultAttribute<Real>("equiv_hf_evals",equiv_hf_evals)});
1330 }
1331
1332 } // namespace Dakota
1333