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: NonDGlobalReliability
10 //- Description: Implementation code for NonDGlobalReliability class
11 //- Owner: Barron Bichon, Mike Eldred
12 //- Checked by:
13 //- Version:
14
15 #include "dakota_system_defs.hpp"
16 #include "dakota_data_io.hpp"
17 #include "NonDGlobalReliability.hpp"
18 #include "NonDLHSSampling.hpp"
19 //#include "DDACEDesignCompExp.hpp"
20 //#include "FSUDesignCompExp.hpp"
21 #include "NonDAdaptImpSampling.hpp"
22 //#ifdef HAVE_ACRO
23 //#include "COLINOptimizer.hpp"
24 //#endif
25 #ifdef HAVE_NCSU
26 #include "NCSUOptimizer.hpp"
27 #endif
28 #include "ProbabilityTransformModel.hpp"
29 #include "DataFitSurrModel.hpp"
30 #include "DakotaApproximation.hpp"
31 #include "ProblemDescDB.hpp"
32 #include "NormalRandomVariable.hpp"
33
34 //#define DEBUG
35 //#define DEGUG_PLOTS
36
37 static const char rcsId[] = "@(#) $Id: NonDGlobalReliability.cpp 4058 2006-10-25 01:39:40Z mseldre $";
38
39 extern "C" {
40 #ifdef DAKOTA_F90
41 #if defined(HAVE_CONFIG_H) && !defined(DISABLE_DAKOTA_CONFIG_H)
42
43 // Deprecated; continue to support legacy, clashing macros for ONE RELEASE
44 #define BVLS_WRAPPER_FC FC_FUNC_(bvls_wrapper,BVLS_WRAPPER)
45 void BVLS_WRAPPER_FC( Dakota::Real* a, int& m, int& n, Dakota::Real* b,
46 Dakota::Real* bnd, Dakota::Real* x, Dakota::Real& rnorm,
47 int& nsetp, Dakota::Real* w, int* index, int& ierr );
48 #else
49
50 // Use the CMake-generated fortran name mangling macros (eliminate warnings)
51 #include "dak_f90_config.h"
52 #define BVLS_WRAPPER_FC DAK_F90_GLOBAL_(bvls_wrapper,BVLS_WRAPPER)
53 void BVLS_WRAPPER_FC( Dakota::Real* a, int& m, int& n, Dakota::Real* b,
54 Dakota::Real* bnd, Dakota::Real* x, Dakota::Real& rnorm,
55 int& nsetp, Dakota::Real* w, int* index, int& ierr );
56 #endif // HAVE_CONFIG_H and !DISABLE_DAKOTA_CONFIG_H
57 #endif // DAKOTA_F90
58 }
59
60 namespace Dakota {
61
62 // initialization of statics
63 NonDGlobalReliability* NonDGlobalReliability::nondGlobRelInstance(NULL);
64
65
66 NonDGlobalReliability::
NonDGlobalReliability(ProblemDescDB & problem_db,Model & model)67 NonDGlobalReliability(ProblemDescDB& problem_db, Model& model):
68 NonDReliability(problem_db, model),
69 meritFunctionType(AUGMENTED_LAGRANGIAN_MERIT), dataOrder(1)
70 {
71 if (mppSearchType < EGRA_X) {
72 Cerr << "Error: only x-space and u-space EGRA are currently supported in "
73 << "global_reliability."<< std::endl;
74 abort_handler(-1);
75 }
76
77 // standard reliability indices are not defined and should be precluded
78 // via the input spec. requestedRelLevels is default sized in NonD and
79 // cannot be used for the empty() test below.
80 if (!probDescDB.get_rva("method.nond.reliability_levels").empty() ||
81 respLevelTarget == RELIABILITIES) {
82 Cerr << "Error: reliability indices are not defined for global reliability "
83 << "methods. Use generalized reliability instead." << std::endl;
84 abort_handler(-1);
85 }
86
87 // requestedProbLevels & requestedGenRelLevels are not yet supported
88 // since PMA EGRA requires additional R&D. Note: requestedProbLevels and
89 // requestedGenRelLevels are default sized in NonD and cannot be used for
90 // the empty() tests below.
91 if (!probDescDB.get_rva("method.nond.probability_levels").empty() ||
92 !probDescDB.get_rva("method.nond.gen_reliability_levels").empty()) {
93 Cerr << "Error: Inverse reliability mappings not currently supported in "
94 << "global_reliability."<< std::endl;
95 abort_handler(-1);
96 }
97
98 #ifndef DAKOTA_F90
99 if (meritFunctionType == LAGRANGIAN_MERIT) {
100 Cerr << "Error: F90 required for standard Lagrangian merit function in "
101 << "global_reliability."<< std::endl;
102 abort_handler(-1);
103 }
104 #endif // DAKOTA_F90
105
106 // Size the output arrays, augmenting sizing in NonDReliability. Relative to
107 // other NonD methods, the output storage for reliability methods is greater
108 // since there may be differences between requested and computed levels for
109 // the same level type (the request is not always achieved) and since
110 // probability and reliability are carried along in parallel (due to their
111 // direct correspondence). Relative to NonDLocalReliability, RelLevels are
112 // ignored since probabilities are estimated directly from MMAIS.
113 size_t i;
114 for (i=0; i<numFunctions; i++) {
115 size_t num_levels = requestedRespLevels[i].length() +
116 requestedProbLevels[i].length() + requestedGenRelLevels[i].length();
117 computedRespLevels[i].resize(num_levels);
118 computedProbLevels[i].resize(num_levels);
119 computedGenRelLevels[i].resize(num_levels);
120 }
121
122 // The Gaussian process model of the limit state in u-space [G-hat(u)] is
123 // constructed here one time.
124
125 // Always build a global Gaussian process model. No correction is needed.
126 String approx_type = "global_kriging";
127 if (probDescDB.get_short("method.nond.emulator") == GP_EMULATOR)
128 approx_type = "global_gaussian";
129 else if (probDescDB.get_short("method.nond.emulator") == EXPGP_EMULATOR)
130 approx_type = "global_exp_gauss_proc";
131
132 unsigned short sample_type = SUBMETHOD_DEFAULT;
133 UShortArray approx_order; // not used for GP/kriging
134 short corr_order = -1, corr_type = NO_CORRECTION,
135 active_view = iteratedModel.current_variables().view().first;
136 if (probDescDB.get_bool("method.derivative_usage")) {
137 if (approx_type == "global_gaussian") {
138 Cerr << "\nError: efficient_global does not support gaussian_process "
139 << "when derivatives present; use kriging instead." << std::endl;
140 abort_handler(-1);
141 }
142 if (iteratedModel.gradient_type() != "none") dataOrder |= 2;
143 if (iteratedModel.hessian_type() != "none") dataOrder |= 4;
144 }
145 String sample_reuse
146 = (active_view == RELAXED_ALL || active_view == MIXED_ALL) ? "all" : "none";
147
148 int db_samples = probDescDB.get_int("method.samples");
149 int samples = (db_samples > 0) ? db_samples :
150 (numContinuousVars+1)*(numContinuousVars+2)/2;
151
152 int lhs_seed = probDescDB.get_int("method.random_seed");
153 const String& rng = probDescDB.get_string("method.random_number_generator");
154 bool vary_pattern = false; // for consistency across outer loop invocations
155 // get point samples file
156 const String& import_pts_file
157 = probDescDB.get_string("method.import_build_points_file");
158 if (!import_pts_file.empty())
159 { samples = 0; sample_reuse = "all"; }
160
161 //int symbols = samples; // symbols needed for DDACE
162 Iterator dace_iterator;
163 // instantiate the Nataf ProbabilityTransform and GP DataFit recursions
164 if (mppSearchType == EGRA_X) { // Recast( DataFit( iteratedModel ) )
165
166 // The following uses on the fly derived ctor:
167 auto lhs_sampler_rep = std::make_shared<NonDLHSSampling>
168 (iteratedModel, sample_type, samples,
169 lhs_seed, rng, vary_pattern, ACTIVE_UNIFORM);
170 //unsigned short dace_method = SUBMETHOD_LHS; // submethod enum
171 //lhs_sampler_rep = new DDACEDesignCompExp(iteratedModel, samples, symbols,
172 // lhs_seed, dace_method);
173 //unsigned short dace_method = FSU_HAMMERSLEY;
174 //lhs_sampler_rep = new FSUDesignCompExp(iteratedModel, samples, lhs_seed,
175 // dace_method);
176 dace_iterator.assign_rep(lhs_sampler_rep);
177
178 // Construct g-hat(x) using a GP approximation over the active/uncertain
179 // vars (same view as iteratedModel: not the typical All view for DACE).
180 // For RBDO with GP over {u}+{d}, view set to All from all_variables spec.
181 Model g_hat_x_model;
182 IntSet surr_fn_indices; // Only want functions with requested levels
183 ActiveSet set = iteratedModel.current_response().active_set();// copy
184 set.request_values(0);
185 for (i=0; i<numFunctions; ++i)
186 if (!computedRespLevels[i].empty()) // sized to total req levels above
187 { set.request_value(dataOrder, i); surr_fn_indices.insert(i); }
188 dace_iterator.active_set(set);
189 //const Variables& curr_vars = iteratedModel.current_variables();
190 ActiveSet gp_set = iteratedModel.current_response().active_set(); // copy
191 gp_set.request_values(1);// no surr deriv evals, but GP may be grad-enhanced
192 g_hat_x_model.assign_rep(std::make_shared<DataFitSurrModel>
193 (dace_iterator, iteratedModel,
194 gp_set, approx_type, approx_order, corr_type, corr_order, dataOrder,
195 outputLevel, sample_reuse, import_pts_file,
196 probDescDB.get_ushort("method.import_build_format"),
197 probDescDB.get_bool("method.import_build_active_only"),
198 probDescDB.get_string("method.export_approx_points_file"),
199 probDescDB.get_ushort("method.export_approx_format")));
200 g_hat_x_model.surrogate_function_indices(surr_fn_indices);
201
202 // Recast g-hat(x) to G-hat(u); truncate dist bnds
203 uSpaceModel.assign_rep(std::make_shared<ProbabilityTransformModel>
204 (g_hat_x_model, STD_NORMAL_U, true, 5.));
205 }
206 else { // DataFit( Recast( iteratedModel ) )
207
208 // Recast g(x) to G(u); truncate dist bnds
209 Model g_u_model;
210 g_u_model.assign_rep(std::make_shared<ProbabilityTransformModel>
211 (iteratedModel, STD_NORMAL_U, true, 5.));
212
213 // For additional generality, could develop on the fly envelope ctor:
214 //Iterator dace_iterator(g_u_model, dace_method, ...);
215
216 // The following use on-the-fly derived ctors:
217 auto lhs_sampler_rep = std::make_shared<NonDLHSSampling>
218 (g_u_model, sample_type,
219 samples, lhs_seed, rng, vary_pattern, ACTIVE_UNIFORM);
220 //unsigned short dace_method = SUBMETHOD_LHS; // submethod enum
221 //lhs_sampler_rep = new DDACEDesignCompExp(g_u_model, samples, symbols,
222 // lhs_seed, dace_method);
223 //unsigned short dace_method = FSU_HAMMERSLEY;
224 //lhs_sampler_rep = new FSUDesignCompExp(g_u_model, samples, lhs_seed,
225 // dace_method);
226 dace_iterator.assign_rep(lhs_sampler_rep);
227
228 // Construct G-hat(u) using a GP approximation over the active/uncertain
229 // variables (using the same view as iteratedModel/g_u_model: not the
230 // typical All view for DACE).
231 // For RBDO with GP over {u}+{d}, view set to All from all_variables spec.
232 IntSet surr_fn_indices; // Only want functions with requested levels
233 ActiveSet set = iteratedModel.current_response().active_set();// copy
234 set.request_values(0);
235 for (i=0; i<numFunctions; ++i)
236 if (!computedRespLevels[i].empty()) // sized to total req levels above
237 { set.request_value(dataOrder, i); surr_fn_indices.insert(i); }
238 dace_iterator.active_set(set);
239
240 //const Variables& g_u_vars = g_u_model.current_variables();
241 ActiveSet gp_set = g_u_model.current_response().active_set(); // copy
242 gp_set.request_values(1);// no surr deriv evals, but GP may be grad-enhanced
243 uSpaceModel.assign_rep(std::make_shared<DataFitSurrModel>
244 (dace_iterator, g_u_model,
245 gp_set, approx_type, approx_order, corr_type, corr_order, dataOrder,
246 outputLevel, sample_reuse, import_pts_file,
247 probDescDB.get_ushort("method.import_build_format"),
248 probDescDB.get_bool("method.import_build_active_only"),
249 probDescDB.get_string("method.export_approx_points_file"),
250 probDescDB.get_ushort("method.export_approx_format")));
251 uSpaceModel.surrogate_function_indices(surr_fn_indices);
252 }
253
254 // Following this ctor, IteratorScheduler::init_iterator() initializes the
255 // parallel configuration for NonDGlobalReliability + iteratedModel using
256 // NonDGlobalReliability's maxEvalConcurrency. During uSpaceModel construction
257 // above, DataFitSurrModel::derived_init_communicators() initializes the
258 // parallel configuration for dace_iterator + iteratedModel using
259 // dace_iterator's maxEvalConcurrency. The only iteratedModel concurrency
260 // currently exercised is that used by dace_iterator within the initial GP
261 // construction, but the NonDGlobalReliability maxEvalConcurrency must still
262 // be set so as to avoid parallel configuration errors resulting from
263 // avail_procs > max_concurrency within IteratorScheduler::init_iterator().
264 // A max of the local derivative concurrency and the DACE concurrency is used
265 // for this purpose.
266 maxEvalConcurrency = std::max(maxEvalConcurrency,
267 dace_iterator.maximum_evaluation_concurrency());
268
269 // Configure a RecastModel with one objective and no constraints using the
270 // alternate minimalist constructor. The RIA/PMA expected improvement/
271 // expected feasibility formulations may vary with the level requests, so
272 // the recast fn pointers are reset for each level within the run fn.
273 SizetArray recast_vars_comps_total; // default: empty; no change in size
274 BitArray all_relax_di, all_relax_dr; // default: empty; no discrete relaxation
275 short recast_resp_order = 1; // nongradient-based optimizers
276 mppModel.assign_rep(std::make_shared<RecastModel>
277 (uSpaceModel, recast_vars_comps_total, all_relax_di,
278 all_relax_dr, 1, 0, 0, recast_resp_order));
279
280 // For formulations with one objective and one equality constraint,
281 // use the following instead:
282 //mppModel.assign_rep(new RecastModel(uSpaceModel, ..., 1, 1, 0, ...), false);
283 //RealVector nln_eq_targets(1, 0.);
284 //mppModel.nonlinear_eq_constraint_targets(nln_eq_targets);
285
286 // must use alternate NoDB ctor chain
287 int max_iter = 1000, max_eval = 10000;
288 double min_box_size = 1.e-15, vol_box_size = 1.e-15;
289 #ifdef HAVE_NCSU
290 mppOptimizer.assign_rep(std::make_shared<NCSUOptimizer>
291 (mppModel, max_iter, max_eval,
292 min_box_size, vol_box_size));
293 //#ifdef HAVE_ACRO
294 //int coliny_seed = 0; // system-generated, for now
295 //mppOptimizer.assign_rep(new
296 // COLINOptimizer<coliny::DIRECT>(mppModel, coliny_seed), false);
297 //mppOptimizer.assign_rep(new
298 // COLINOptimizer<coliny::EAminlp>(mppModel, coliny_seed), false);
299 //#endif
300 #else
301 Cerr << "NCSU DIRECT Optimizer is not available to use in the MPP search "
302 << "in global reliability optimization: aborting process." << std::endl;
303 abort_handler(-1);
304 #endif //HAVE_NCSU
305
306 // The importance sampler uses uSpaceModel (without additional recasting)
307 // and may be constructed here. Thus, NonDGlobal applies integration
308 // refinement to the G-hat(u) surrogate. Behavior needs to be repeatable
309 // and AIS is not part of the EGRA spec: either reuse lhs_seed or hardwire.
310 int refine_samples = 1000, refine_seed = 123457;
311 // we pass a u-space model and enforce the EGRA GP bounds on the samples;
312 // extreme values are needed to define bounds for outer PDF bins.
313 bool x_model_flag = false, use_model_bounds = true, track_extreme = pdfOutput;
314 integrationRefinement = MMAIS; vary_pattern = true;
315
316 auto importance_sampler_rep = std::make_shared<NonDAdaptImpSampling>
317 (uSpaceModel, sample_type, refine_samples, refine_seed,
318 rng, vary_pattern, integrationRefinement, cdfFlag,
319 x_model_flag, use_model_bounds, track_extreme);
320 importanceSampler.assign_rep(importance_sampler_rep);
321 }
322
323
~NonDGlobalReliability()324 NonDGlobalReliability::~NonDGlobalReliability()
325 { }
326
327
resize()328 bool NonDGlobalReliability::resize()
329 {
330 bool parent_reinit_comms = NonDReliability::resize();
331
332 Cerr << "\nError: Resizing is not yet supported in method "
333 << method_enum_to_string(methodName) << "." << std::endl;
334 abort_handler(METHOD_ERROR);
335
336 return parent_reinit_comms;
337 }
338
339
derived_init_communicators(ParLevLIter pl_iter)340 void NonDGlobalReliability::derived_init_communicators(ParLevLIter pl_iter)
341 {
342 iteratedModel.init_communicators(pl_iter, maxEvalConcurrency);
343
344 // mppModel.init_communicators() recursion is currently sufficient for
345 // uSpaceModel. An additional uSpaceModel.init_communicators() call would be
346 // motivated by special parallel usage of uSpaceModel below that is not
347 // otherwise covered by the recursion.
348 //uSpaceMaxConcurrency = maxEvalConcurrency; // local derivative concurrency
349 //uSpaceModel.init_communicators(pl_iter, uSpaceMaxConcurrency);
350
351 // mppOptimizer and importanceSampler use NoDBBaseConstructor, so no
352 // need to manage DB list nodes at this level
353 mppOptimizer.init_communicators(pl_iter);
354 importanceSampler.init_communicators(pl_iter);
355 }
356
357
derived_set_communicators(ParLevLIter pl_iter)358 void NonDGlobalReliability::derived_set_communicators(ParLevLIter pl_iter)
359 {
360 NonD::derived_set_communicators(pl_iter);
361
362 //uSpaceMaxConcurrency = maxEvalConcurrency; // local derivative concurrency
363 //uSpaceModel.set_communicators(pl_iter, uSpaceMaxConcurrency);
364
365 // mppOptimizer and importanceSampler use NoDBBaseConstructor, so no
366 // need to manage DB list nodes at this level
367 mppOptimizer.set_communicators(pl_iter);
368 importanceSampler.set_communicators(pl_iter);
369 }
370
371
derived_free_communicators(ParLevLIter pl_iter)372 void NonDGlobalReliability::derived_free_communicators(ParLevLIter pl_iter)
373 {
374 // deallocate communicators for MMAIS on uSpaceModel
375 importanceSampler.free_communicators(pl_iter);
376
377 // deallocate communicators for DIRECT on mppModel
378 mppOptimizer.free_communicators(pl_iter);
379
380 //uSpaceMaxConcurrency = maxEvalConcurrency; // local derivative concurrency
381 //uSpaceModel.free_communicators(pl_iter, uSpaceMaxConcurrency);
382
383 iteratedModel.free_communicators(pl_iter, maxEvalConcurrency);
384 }
385
386
pre_run()387 void NonDGlobalReliability::pre_run()
388 {
389 Analyzer::pre_run();
390
391 // IteratorScheduler::run_iterator() + Analyzer::initialize_run() ensure
392 // initialization of Model mappings for iteratedModel, but local recursions
393 // are not visible -> recur DataFitSurr + ProbabilityTransform if needed.
394 // > Note: part of this occurs at DataFit build time. Therefore, take
395 // care to avoid redundancy using mappingInitialized flag.
396 if (!mppModel.mapping_initialized()) {
397 ParLevLIter pl_iter = methodPCIter->mi_parallel_level_iterator(miPLIndex);
398 /*bool var_size_changed =*/ mppModel.initialize_mapping(pl_iter);
399 //if (var_size_changed) resize();
400 }
401
402 // now that vars/labels/bounds/targets have flowed down at run-time from
403 // any higher level recursions, propagate them up local Model recursions
404 // so that they are correct when they propagate back down. There is no
405 // need to recur below iteratedModel.
406 size_t layers = 3; // DFSModel + PTModel + RModel
407 mppModel.update_from_subordinate_model(layers-1); // recur twice
408 }
409
410
core_run()411 void NonDGlobalReliability::core_run()
412 {
413 // set the object instance pointer for use within static member functions
414 NonDGlobalReliability* prev_grel_instance = nondGlobRelInstance;
415 nondGlobRelInstance = this;
416
417 // Optimize the GP for all levels and then use MAIS for _all_ levels
418 optimize_gaussian_process();
419 importance_sampling();
420
421 // restore in case of recursion
422 nondGlobRelInstance = prev_grel_instance;
423 }
424
425
optimize_gaussian_process()426 void NonDGlobalReliability::optimize_gaussian_process()
427 {
428 if (mppSearchType == EGRA_X) {
429 // Assign non-default global variable bounds for use in PStudyDACE methods
430 // that require a bounded region (NIDR default truncates infinite/semi-
431 // infinite tails using 3-sigma by default to allow use outside NonD). As
432 // defined for u-space in the ctor, we use a 5-sigma truncation for EGRA.
433 // Note: This does not affect any uncertain variable distribution bounds
434 // or how NonD methods sample. However, PStudyDACE will sample uniformly
435 // in this x-space interval, which could result in very irregular coverage
436 // in u-space (it would be preferable to sample uniformly in u-space).
437 // Note: this cannot currently be defined in the constructor. EGRA_U does
438 // this by defining bounds truncation in ProbabilityTransformModel, but
439 // EGRA_X needs trans_U_to_X() which isn't available until dist parameters
440 // are updated at run time. Therefore, the code below sets the bounds
441 // for the DataFitSurrModel, which then propagates to actualModel in
442 // DataFitSurrModel::update_actual_model()
443 // Note: since u-space type is STD_NORMAL_U, u-space aleatory variables
444 // are all unbounded and global bounds are set to +/-5.
445 Pecos::ProbabilityTransformation& nataf
446 = uSpaceModel.probability_transformation();
447 RealVector x_l_bnds, x_u_bnds;
448 nataf.trans_U_to_X(uSpaceModel.continuous_lower_bounds(), x_l_bnds);
449 nataf.trans_U_to_X(uSpaceModel.continuous_upper_bounds(), x_u_bnds);
450 Model& g_hat_x_model = uSpaceModel.subordinate_model();
451 g_hat_x_model.continuous_lower_bounds(x_l_bnds);
452 g_hat_x_model.continuous_upper_bounds(x_u_bnds);
453 }
454
455 // Build initial GP once for all response functions
456 uSpaceModel.build_approximation();
457
458 // Loop over each response function in the responses specification. It is
459 // important to note that the MPP iteration is different for each response
460 // function, and it is not possible to combine the model evaluations for
461 // multiple response functions.
462 ParLevLIter pl_iter = methodPCIter->mi_parallel_level_iterator(miPLIndex);
463 for (respFnCount=0; respFnCount<numFunctions; respFnCount++) {
464
465 // The most general case is to allow a combination of response, probability,
466 // and reliability level specifications for each response function.
467 size_t rl_len = requestedRespLevels[respFnCount].length(),
468 pl_len = requestedProbLevels[respFnCount].length(),
469 gl_len = requestedGenRelLevels[respFnCount].length(),
470 num_levels = rl_len + pl_len + gl_len;
471
472 // Loop over response/probability/reliability levels
473 for (levelCount=0; levelCount<num_levels; levelCount++) {
474
475 // The rl_len response levels are performed first using the RIA
476 // formulation, followed by the pl_len probability levels and the
477 // gl_len generalized reliability levels using the PMA formulation.
478 bool ria_flag = (levelCount < rl_len) ? true : false;
479 if (ria_flag) {
480 requestedTargetLevel = requestedRespLevels[respFnCount][levelCount];
481 Cout << "\n>>>>> Reliability Index Approach (RIA) for response level "
482 << levelCount+1 << " = " << requestedTargetLevel << '\n';
483 }
484 else if (levelCount < rl_len + pl_len) {
485 size_t index = levelCount-rl_len;
486 Real p = requestedProbLevels[respFnCount][index];
487 Cout << "\n>>>>> Performance Measure Approach (PMA) for probability "
488 << "level " << index+1 << " = " << p << '\n';
489 requestedTargetLevel = -Pecos::NormalRandomVariable::inverse_std_cdf(p);
490 Real gen_beta_cdf = (cdfFlag) ?
491 requestedTargetLevel : -requestedTargetLevel;
492 pmaMaximizeG = (gen_beta_cdf < 0.);
493 }
494 else {
495 size_t index = levelCount-rl_len-pl_len;
496 requestedTargetLevel = requestedGenRelLevels[respFnCount][index];
497 Cout << "\n>>>>> Performance Measure Approach (PMA) for reliability "
498 << "level " << index+1 << " = " << requestedTargetLevel << '\n';
499 Real gen_beta_cdf = (cdfFlag) ?
500 requestedTargetLevel : -requestedTargetLevel;
501 pmaMaximizeG = (gen_beta_cdf < 0.);
502 }
503
504 bool pma_aug_lag_flag
505 = (!ria_flag && meritFunctionType == AUGMENTED_LAGRANGIAN_MERIT);
506 if (pma_aug_lag_flag) {
507 augLagrangeMult = 0.; penaltyParameter = 1.;
508 lastConstraintViolation = DBL_MAX; lastIterateAccepted = false;
509 }
510
511 // Iterate until EGRA converges
512 approxIters = 0;
513 approxConverged = false;
514 while (!approxConverged) {
515
516 approxIters++;
517
518 // construct global optimizer and its EI/EF recast model
519 Sizet2DArray vars_map, primary_resp_map, secondary_resp_map;
520 BoolDequeArray nonlinear_resp_map(1, BoolDeque(1, true));
521 std::shared_ptr<RecastModel> mpp_model_rep =
522 std::static_pointer_cast<RecastModel>(mppModel.model_rep());
523 if (ria_flag) {
524 // Standard RIA : min u'u s.t. g = z_bar
525 // use RIA evaluators to recast g into global opt subproblem
526 //primary_resp_map.reshape(1); // one objective, no contributors
527 //secondary_resp_map.reshape(1); // one constraint, one contributor
528 //secondary_resp_map[0].reshape(1);
529 //secondary_resp_map[0][0] = respFnCount;
530 //BoolDequeArray nonlinear_resp_map(2);
531 //nonlinear_resp_map[1] = BoolDeque(1, false);
532 //mpp_model_rep->init_maps(vars_map, false, NULL, NULL,
533 // primary_resp_map, secondary_resp_map, nonlinear_resp_map,
534 // RIA_objective_eval, RIA_constraint_eval);
535
536 // EFF formulation : max EFF s.t. bound constraints
537 // use EFF evaluators to recast g into global opt subproblem
538 primary_resp_map.resize(1);
539 primary_resp_map[0].resize(1);
540 primary_resp_map[0][0] = respFnCount;
541 mpp_model_rep->init_maps(vars_map, false, NULL, NULL,
542 primary_resp_map, secondary_resp_map, nonlinear_resp_map,
543 EFF_objective_eval, NULL);
544 }
545 else {
546 // Standard PMA : min/max g s.t. u'u = beta_bar^2
547 // use PMA evaluators to recast g into global opt subproblem
548 //void (*set_map) (const ShortArray& recast_asv,
549 // ShortArray& sub_model_asv) =
550 // (integrationOrder == 2) ? PMA2_set_mapping : NULL;
551 //primary_resp_map.reshape(1); // one objective, one contributor
552 //primary_resp_map[0].reshape(1);
553 //primary_resp_map[0][0] = respFnCount;
554 //secondary_resp_map.reshape(1); // one constraint, no contributors
555 //BoolDequeArray nonlinear_resp_map(2);
556 //nonlinear_resp_map[0] = BoolDeque(1, false);
557 //mpp_model_rep->init_maps(vars_map, false, NULL, set_map,
558 // primary_resp_map, secondary_resp_map, nonlinear_resp_map,
559 // PMA_objective_eval, PMA_constraint_eval);
560
561 // EIF formulation : max Phi(EIF) [merit function on EIF]
562 // determine fnStar from among sample data
563 get_best_sample();
564 // use EIF evaluators to recast g into global opt subproblem
565 primary_resp_map.resize(1);
566 primary_resp_map[0].resize(1);
567 primary_resp_map[0][0] = respFnCount;
568 mpp_model_rep->init_maps(vars_map, false, NULL, NULL,
569 primary_resp_map, secondary_resp_map, nonlinear_resp_map,
570 EIF_objective_eval, NULL);
571 }
572
573 // Execute GLOBAL search and retrieve u-space results
574 Cout << "\n>>>>> Initiating global reliability optimization\n";
575 mppOptimizer.run(pl_iter);
576 // Use these two lines for COLINY optimizers
577 //const VariablesArray& vars_star
578 // = mppOptimizer.variables_array_results();
579 //const RealVector& c_vars_u = vars_star[0].continuous_variables();
580 // Use these two lines for NCSU DIRECT
581 const Variables& vars_star = mppOptimizer.variables_results();
582 const RealVector& c_vars_u = vars_star.continuous_variables();
583
584 // Get expected value at u* for output
585 uSpaceModel.continuous_variables(c_vars_u);
586 uSpaceModel.evaluate();
587 const RealVector& g_hat_fns
588 = uSpaceModel.current_response().function_values();
589
590 // Re-evaluate the expected improvement/feasibility at vars_star
591 Real beta_star = 0.,/* TO DO */ exp_fns_star = (ria_flag) ?
592 expected_feasibility(g_hat_fns, vars_star) :
593 expected_improvement(g_hat_fns, vars_star);
594
595 Cout << "\nResults of EGRA iteration:\nFinal point (u-space) =\n"
596 << c_vars_u;
597 size_t wpp7 = write_precision+7;
598 if (ria_flag) {
599 Cout << "Expected Feasibility =\n "
600 << std::setw(wpp7) << -exp_fns_star << "\n "
601 << std::setw(wpp7) << g_hat_fns[respFnCount]-requestedTargetLevel
602 << " [G_hat(u) - z]\n";
603 // << " " << std::setw(wpp7) << beta_star
604 // << " [beta*]\n";
605 //Cout << "RIA optimum =\n "
606 // << std::setw(wpp7) << exp_fns_star << " [u'u]\n"
607 // << " " << std::setw(wpp7) << exp_fns_star[1]
608 // << " [G(u) - z]\n";
609 }
610 else {
611 // Calculate beta^2 for output (and aug_lag update)
612 Cout << "Expected Improvement =\n "
613 << std::setw(wpp7) << -exp_fns_star << "\n "
614 << std::setw(wpp7) << beta_star - requestedTargetLevel
615 << " [beta* - bar-beta*]\n "
616 << std::setw(wpp7) << g_hat_fns[respFnCount] << " [G_hat(u)]\n";
617 //Cout << "PMA optimum =\n "
618 // << std::setw(wpp7) << exp_fns_star << " [";
619 //if (pmaMaximizeG) Cout << '-';
620 //Cout << "G(u)]\n " << std::setw(wpp7)
621 // << exp_fns_star[1] << " [u'u - B^2]\n";
622 }
623
624 // Update parameters for the augmented Lagrangian merit function
625 if (pma_aug_lag_flag) {
626 // currently only used for PMA with EIF
627 Real c_violation = beta_star - requestedTargetLevel;
628 if (c_violation < lastConstraintViolation)
629 lastIterateAccepted = true;
630 lastConstraintViolation = c_violation;
631 }
632
633 // Check for convergence based on max EIF/EFF
634 // BMA: was previously hard-wired: convergenceTol = .001;
635 if (maxIterations < 0)
636 maxIterations = 25*numContinuousVars;
637 if (approxIters >= maxIterations || -exp_fns_star < convergenceTol)
638 approxConverged = true;
639 else if (mppSearchType == EGRA_X) {
640 // Evaluate response_star_truth in x-space
641 x_truth_evaluation(c_vars_u, dataOrder);
642 // Update the GP approximation in x-space
643 IntResponsePair resp_star_truth(iteratedModel.evaluation_id(),
644 iteratedModel.current_response());
645 uSpaceModel.append_approximation(iteratedModel.current_variables(),
646 resp_star_truth, true);
647 }
648 else {
649 // Evaluate response_star_truth in u-space
650 u_truth_evaluation(c_vars_u, dataOrder);
651 // Update the GP approximation in u-space
652 IntResponsePair resp_star_truth(uSpaceModel.evaluation_id(),
653 uSpaceModel.current_response());
654 uSpaceModel.append_approximation(vars_star, resp_star_truth, true);
655 }
656 } // end approx convergence while loop
657
658 if (ria_flag)
659 Cout << "\n<<<<< GP model has converged for response level ";
660 else
661 Cout << "\n<<<<< GP model has converged for response level ";
662 Cout << levelCount+1 << " = " << requestedTargetLevel << '\n';
663
664 } // end loop over levels
665
666 if (num_levels)
667 Cout << "\n<<<<< GP model has converged for all requested levels of "
668 << "response function " << respFnCount+1 << '\n';
669
670 #ifdef DEBUG
671 // DEBUG - output set of samples used to build the GP
672 // If problem is 2d, output a grid of points on the GP, variance,
673 // eff, and truth (if requested)
674 if (num_levels) { // can only plot if GP was built!
675 std::string samsfile("egra_sams");
676 std::string tag = "_" + std::to_string(respFnCount+1) + ".out";
677 samsfile += tag;
678 std::ofstream samsOut(samsfile.c_str(),std::ios::out);
679 samsOut << std::scientific;
680 const Pecos::SurrogateData& gp_data
681 = uSpaceModel.approximation_data(respFnCount);
682 size_t num_data_pts = gp_data.size(), num_vars = uSpaceModel.cv();
683 for (size_t i=0; i<num_data_pts; ++i) {
684 const RealVector& sams = gp_data.continuous_variables(i); // view
685 Real true_fn = gp_data.response_function(i);
686
687 if (mppSearchType == EGRA_X) {
688 RealVector sams_u(num_vars);
689 uSpaceModel.probability_transformation().trans_X_to_U(sams,sams_u);
690
691 samsOut << '\n';
692 for (size_t j=0; j<num_vars; j++)
693 samsOut << std::setw(13) << sams_u[j] << ' ';
694 samsOut<< std::setw(13) << true_fn;
695 }
696 else {
697 samsOut << '\n';
698 for (size_t j=0; j<num_vars; j++)
699 samsOut << std::setw(13) << sams[j] << ' ';
700 samsOut<< std::setw(13) << true_fn;
701 }
702 }
703 samsOut << std::endl;
704
705 #ifdef DEBUG_PLOTS
706 // Plotting the GP, etc over a grid is intended for visualization and
707 // is therefore only available for 2D problems
708 if (num_vars==2) {
709 bool true_plot = false;
710 std::string truefile("egra_true"), gpfile("egra_gp"),
711 varfile("egra_var"), efffile("egra_eff");
712 truefile += tag; gpfile += tag; varfile += tag; efffile += tag;
713 std::ofstream trueOut(truefile.c_str(), std::ios::out),
714 gpOut(gpfile.c_str(), std::ios::out),
715 varOut(varfile.c_str(), std::ios::out),
716 effOut(efffile.c_str(), std::ios::out);
717 gpOut << std::scientific; varOut << std::scientific;
718 effOut << std::scientific; if (true_plot) trueOut << std::scientific;
719 RealVector u_pt(2), x_pt(2);
720 Real lbnd = -5., ubnd = 5.;
721 //Real lbnd = -10., ubnd = 10.;
722 Real interval = (ubnd-lbnd)/100.;
723 for (size_t i=0; i<101; i++){
724 u_pt[0] = lbnd + float(i)*interval;
725 for (size_t j=0; j<101; j++){
726 u_pt[1] = lbnd + float(j)*interval;
727
728 u_evaluation(u_pt, 1);
729 const Response& gp_resp = uSpaceModel.current_response();
730 const RealVector& gp_fns = gp_resp.function_values();
731 gpOut << '\n' << std::setw(13) << u_pt[0] << ' ' << std::setw(13)
732 << u_pt[1] << ' ' << std::setw(13) << gp_fns[respFnCount];
733
734 RealVector variance;
735 if (mppSearchType == EGRA_X) { // Recast( DataFit( iteratedModel ) )
736 // RecastModel::derived_evaluate() propagates u_pt to x_pt
737 Model& dfs_model = uSpaceModel.subordinate_model();
738 variance = dfs_model.approximation_variances(
739 dfs_model.current_variables()); // x_pt
740 }
741 else // EGRA_U: DataFit( Recast( iteratedModel ) )
742 variance = uSpaceModel.approximation_variances(
743 uSpaceModel.current_variables()); // u_pt
744
745 varOut << '\n' << std::setw(13) << u_pt[0] << ' ' << std::setw(13)
746 << u_pt[1] << ' ' << std::setw(13) << variance[respFnCount];
747
748 Real eff = expected_feasibility(gp_fns, u_pt);
749
750 effOut << '\n' << std::setw(13) << u_pt[0] << ' ' << std::setw(13)
751 << u_pt[1] << ' ' << std::setw(13) << -eff;
752
753 // plotting the true function can be expensive, but is available
754 if (true_plot) {
755 x_truth_evaluation(u_pt, 1);
756 const Response& x_resp = iteratedModel.current_response();
757 trueOut << '\n' << std::setw(13) << u_pt[0] << ' '
758 << std::setw(13) << u_pt[1] << ' ' << std::setw(13)
759 << x_resp.function_value(respFnCount);
760 }
761 }
762 gpOut << std::endl; varOut << std::endl; effOut << std::endl;
763 if (true_plot) trueOut << std::endl;
764 }
765 }
766 #endif //DEBUG_PLOTS
767 }
768 #endif //DEBUG
769
770 } // end loop over response fns
771
772 // (conditionally) export final surrogates after all responses built
773 // User might expect x-space, but will get u-space if they requested it
774 if (mppSearchType == EGRA_X) {
775 // uSpaceModel = Recast(DataFit(iteratedModel))
776 Model& dfs_model = uSpaceModel.subordinate_model();
777 export_final_surrogates(dfs_model);
778 }
779 else
780 export_final_surrogates(uSpaceModel);
781 }
782
783
importance_sampling()784 void NonDGlobalReliability::importance_sampling()
785 {
786 bool x_data_flag = (mppSearchType == EGRA_X);
787 size_t i;
788 statCount = 0;
789 const ShortArray& final_res_asv = finalStatistics.active_set_request_vector();
790 ParLevLIter pl_iter = methodPCIter->mi_parallel_level_iterator(miPLIndex);
791 // rep needed for access to functions not mapped to Iterator level
792 std::shared_ptr<NonDAdaptImpSampling> importance_sampler_rep =
793 std::static_pointer_cast<NonDAdaptImpSampling>
794 (importanceSampler.iterator_rep());
795
796 for (respFnCount=0; respFnCount<numFunctions; respFnCount++) {
797
798 // The most general case is to allow a combination of response, probability,
799 // and reliability level specifications for each response function.
800 size_t rl_len = requestedRespLevels[respFnCount].length(),
801 pl_len = requestedProbLevels[respFnCount].length(),
802 gl_len = requestedGenRelLevels[respFnCount].length(),
803 num_levels = rl_len + pl_len + gl_len;
804
805 RealVectorArray gp_inputs;
806 if (num_levels==0) {
807 x_truth_evaluation(1); // c_vars_u?
808 const Response& true_resp = iteratedModel.current_response();
809 finalStatistics.function_value(
810 true_resp.function_value(respFnCount), statCount);
811 }
812 else {
813 // extract the approximation data from the surrogate model.
814 // Note 1: this data is either x-space (EGRA_X) or u-space (EGRA_U).
815 // Note 2: this returns _all_ truth model data for this response fn,
816 // including initial DACE and EGRA added points for _all_ levels.
817 // TO DO: likely need to remove DACE and partition remaining data among
818 // levels for importance sampling efficiency.
819 const Pecos::SurrogateData& gp_data
820 = uSpaceModel.approximation_data(respFnCount);
821 const Pecos::SDVArray& sdv_array = gp_data.variables_data();
822 size_t num_data_pts = sdv_array.size();
823 gp_inputs.resize(num_data_pts);
824 for (i=0; i<num_data_pts; ++i)
825 gp_inputs[i] = sdv_array[i].continuous_variables(); // view OK
826 }
827 statCount++;
828
829 // Standard deviations & sensitivities not available, skip
830 statCount++;
831
832 // Loop over response/probability/reliability levels
833 for (levelCount=0; levelCount<num_levels; levelCount++) {
834
835 Cout << "\n<<<<< Performing importance sampling for response function "
836 << respFnCount+1 << " level " << levelCount+1 << '\n';
837
838 bool ria_flag = (levelCount < rl_len) ? true : false;
839 if (ria_flag) {
840 Real z = computedRespLevels[respFnCount][levelCount]
841 = requestedRespLevels[respFnCount][levelCount];
842 importance_sampler_rep->
843 initialize(gp_inputs, x_data_flag, respFnCount, 0., z);
844 }
845 else // not operational (see error traps in ctor)
846 importance_sampler_rep->initialize(gp_inputs, x_data_flag,
847 respFnCount, 0., computedRespLevels[respFnCount][levelCount]);
848
849 importanceSampler.run(pl_iter);
850
851 Real p = importance_sampler_rep->final_probability();
852 #ifdef DEBUG
853 Cout << "\np = " << p << std::endl;
854 #endif // DEBUG
855 // RIA z-bar -> p
856 computedProbLevels[respFnCount][levelCount] = p;
857 // RIA z-bar -> generalized beta
858 Real gen_beta = -Pecos::NormalRandomVariable::inverse_std_cdf(p);
859 computedGenRelLevels[respFnCount][levelCount] = gen_beta;
860 switch (respLevelTarget) {
861 case PROBABILITIES:
862 finalStatistics.function_value(p, statCount++); break;
863 case GEN_RELIABILITIES:
864 finalStatistics.function_value(gen_beta, statCount++); break;
865 }
866 }
867 }
868 // post-process level mappings to define PDFs (using prob_refined and
869 // all_levels_computed modes)
870 if (pdfOutput)
871 compute_densities(importance_sampler_rep->extreme_values(), true, true);
872 }
873
874
875 void NonDGlobalReliability::
EIF_objective_eval(const Variables & sub_model_vars,const Variables & recast_vars,const Response & sub_model_response,Response & recast_response)876 EIF_objective_eval(const Variables& sub_model_vars,
877 const Variables& recast_vars,
878 const Response& sub_model_response,
879 Response& recast_response)
880 {
881 const ShortArray& recast_asv = recast_response.active_set_request_vector();
882 if (recast_asv[0] & 1) {
883 Real ei = nondGlobRelInstance->expected_improvement(
884 sub_model_response.function_values(), recast_vars);
885 recast_response.function_value(ei, 0);
886 }
887 }
888
889
890 void NonDGlobalReliability::
EFF_objective_eval(const Variables & sub_model_vars,const Variables & recast_vars,const Response & sub_model_response,Response & recast_response)891 EFF_objective_eval(const Variables& sub_model_vars,
892 const Variables& recast_vars,
893 const Response& sub_model_response,
894 Response& recast_response)
895 {
896 const ShortArray& recast_asv = recast_response.active_set_request_vector();
897 if (recast_asv[0] & 1) {
898 Real ef = nondGlobRelInstance->expected_feasibility(
899 sub_model_response.function_values(), recast_vars);
900 recast_response.function_value(ef, 0);
901 }
902 }
903
904
905 Real NonDGlobalReliability::
expected_improvement(const RealVector & expected_values,const Variables & recast_vars)906 expected_improvement(const RealVector& expected_values,
907 const Variables& recast_vars)
908 {
909 // Get variance from the GP; Expected values are passed in
910 // If GP built in x-space, transform input point to x-space to get variance
911 RealVector variances;
912 if (mppSearchType == EGRA_X) { // uSpaceModel = Recast(DataFit(iteratedModel))
913 Model& dfs_model = uSpaceModel.subordinate_model();
914 // assume recast_vars have been propagated to GP just prior to call
915 variances
916 = dfs_model.approximation_variances(dfs_model.current_variables());
917 }
918 else // EGRA_U: uSpaceModel = DataFit(Recast(iteratedModel))
919 variances = uSpaceModel.approximation_variances(recast_vars);
920
921 Real mean = expected_values[respFnCount];
922 Real stdv = std::sqrt(variances[respFnCount]);
923
924 // Calculate and apply penalty to the mean
925 Real beta_star = 0.; // TO DO
926 // calculate the equality constraint: beta* - bar-beta*
927 Real cfn = beta_star - requestedTargetLevel;
928 Real penalty = constraint_penalty(cfn, recast_vars.continuous_variables());
929 // Calculation of EI will have different form for the CDF/CCDF cases
930 Real penalized_mean = (pmaMaximizeG) ? mean - penalty : mean + penalty;
931
932 // Calculate the expected improvement
933 // Note: This is independent of the CDF/CCDF check
934
935 Real ei, cdf, pdf;
936 Real snv = (fnStar-penalized_mean); // not normalized yet
937 if(std::fabs(snv)>=std::fabs(stdv)*50.0) {
938 //this will trap the denominator=0.0 case even if numerator=0.0
939 pdf = 0.;
940 cdf = (snv > 0.) ? 1. : 0.;
941 }
942 else{
943 snv /= stdv; // now snv is the standard normal variate
944 cdf = Pecos::NormalRandomVariable::std_cdf(snv);
945 pdf = Pecos::NormalRandomVariable::std_pdf(snv);
946 }
947 ei = (pmaMaximizeG) ? (penalized_mean - fnStar)*(1.-cdf) + stdv*pdf
948 : (fnStar - penalized_mean)*cdf + stdv*pdf;
949 return -ei; // return -EI because we are maximizing EI
950 }
951
952
953 Real NonDGlobalReliability::
expected_feasibility(const RealVector & expected_values,const Variables & recast_vars)954 expected_feasibility(const RealVector& expected_values,
955 const Variables& recast_vars)
956 {
957 // Get variance from the GP; Expected values are passed in
958 // If GP built in x-space, transform input point to x-space to get variance
959 RealVector variances;
960 if (mppSearchType == EGRA_X) { // uSpaceModel = Recast(DataFit(iteratedModel))
961 Model& dfs_model = uSpaceModel.subordinate_model();
962 // assume recast_vars have been propagated to GP just prior to call
963 variances
964 = dfs_model.approximation_variances(dfs_model.current_variables());
965 }
966 else // EGRA_U: uSpaceModel = DataFit(Recast(iteratedModel))
967 variances = uSpaceModel.approximation_variances(recast_vars);
968
969 Real mean = expected_values[respFnCount],
970 stdv = std::sqrt(variances[respFnCount]),
971 zbar = requestedTargetLevel,
972 alpha = 2.; // may want to try values other than 2
973
974 // calculate standard normal variate +/- alpha
975
976 Real cdfz, pdfz, cdfp, pdfp, cdfm, pdfm;
977 Real snvz = (zbar - mean);
978 if (std::fabs(snvz) >= std::fabs(stdv)*50.) {
979 pdfm = pdfp = pdfz = 0.;
980 cdfm = cdfp = cdfz = (snvz > 0.) ? 1. : 0.;
981 }
982 else {
983 snvz /= stdv; Real snvp = snvz + alpha, snvm = snvz - alpha;
984 pdfz = Pecos::NormalRandomVariable::std_pdf(snvz);
985 cdfz = Pecos::NormalRandomVariable::std_cdf(snvz);
986 pdfp = Pecos::NormalRandomVariable::std_pdf(snvp);
987 cdfp = Pecos::NormalRandomVariable::std_cdf(snvp);
988 pdfm = Pecos::NormalRandomVariable::std_pdf(snvm);
989 cdfm = Pecos::NormalRandomVariable::std_cdf(snvm);
990 }
991 // calculate expected feasibility function
992 Real ef = (mean - zbar)*(2.*cdfz - cdfm - cdfp) //exploit
993 - stdv*(2.*pdfz - pdfm - pdfp //explore
994 - alpha*cdfp + alpha*cdfm);
995
996 return -ef; // return -EF because we are maximizing
997 }
998
999
get_best_sample()1000 void NonDGlobalReliability::get_best_sample()
1001 {
1002 // Pull the samples and responses from data used to build latest GP
1003 // and apply any penalties to calculate fnStar for use in the
1004 // expected improvement function
1005 // This is only done for PMA - there is no "best solution" for
1006 // the expected feasibility function used in RIA
1007
1008 Iterator& dace_iterator = uSpaceModel.subordinate_iterator();
1009 const RealMatrix& true_vars_x = dace_iterator.all_samples();
1010 const IntResponseMap& true_responses = dace_iterator.all_responses();
1011 size_t i, j, num_samples = true_vars_x.numCols(),
1012 num_vars = true_vars_x.numRows();
1013
1014 // If GP built in x-space, transform true_vars_x to u-space to calculate beta
1015 RealVectorArray true_c_vars_u(num_samples); RealVector true_vars_x_cv;
1016 for (i=0; i<num_samples; i++) {
1017 true_vars_x_cv = Teuchos::getCol(Teuchos::View,
1018 const_cast<RealMatrix&>(true_vars_x), (int)i);
1019 if (mppSearchType == EGRA_X)
1020 uSpaceModel.probability_transformation().trans_X_to_U(true_vars_x_cv,
1021 true_c_vars_u[i]);
1022 else
1023 true_c_vars_u[i] = true_vars_x_cv; // view OK
1024 }
1025
1026 // Calculation of fnStar will have different form for CDF/CCDF cases
1027 fnStar = (pmaMaximizeG) ? -DBL_MAX : DBL_MAX; IntRespMCIter it;
1028 for (i=0, it=true_responses.begin(); i<num_samples; i++, ++it) {
1029 // calculate the reliability index (beta)
1030 Real beta_star = 0., penalized_response; // TO DO
1031 // calculate the equality constraint: u'u - beta_target^2
1032 Real cfn = beta_star - requestedTargetLevel;
1033 Real penalty = constraint_penalty(cfn, true_c_vars_u[i]);
1034 penalized_response = (pmaMaximizeG) ?
1035 it->second.function_value(0) - penalty :
1036 it->second.function_value(0) + penalty;
1037 if ( ( pmaMaximizeG && penalized_response > fnStar) ||
1038 (!pmaMaximizeG && penalized_response < fnStar) )
1039 fnStar = penalized_response;
1040 }
1041 }
1042
1043
1044 Real NonDGlobalReliability::
constraint_penalty(const Real & c_viol,const RealVector & u)1045 constraint_penalty(const Real& c_viol, const RealVector& u)
1046 {
1047 if (meritFunctionType == PENALTY_MERIT)
1048 return exp(approxIters/10.)*c_viol*c_viol; // try other schedules
1049 else if (meritFunctionType == AUGMENTED_LAGRANGIAN_MERIT) {
1050 if (lastIterateAccepted)
1051 augLagrangeMult += 2.*penaltyParameter*c_viol;
1052 else
1053 penaltyParameter *= 2.;
1054 return augLagrangeMult*c_viol + penaltyParameter*c_viol*c_viol;
1055 }
1056 else if (meritFunctionType == LAGRANGIAN_MERIT) {
1057 #ifdef DAKOTA_F90
1058 // form [A] = grad[u'u - beta^2]
1059 int m = u.length();
1060 RealVector A(m, false);
1061 for (size_t i=0; i<m; i++)
1062 A[i] = 2.*u[i];
1063
1064 // form -{grad_f} = m_grad_f = -grad[G_hat(u)]
1065 uSpaceModel.continuous_variables(u);
1066 uSpaceModel.evaluate();
1067 RealVector m_grad_f
1068 = uSpaceModel.current_response().function_gradient_copy(0);
1069 m_grad_f.scale(-1.);
1070
1071 // solve for lambda : [A]{lambda} = {m_grad_f}
1072 int ierr, nsetp, n = 1; Real res_norm;
1073 IntVector index(1); RealVector lambda(1), w(1), bnd(2);
1074 // lawson_hanson2.f90: BVLS ignore bounds based on huge(), so +/-DBL_MAX
1075 // is sufficient here
1076 bnd[0] = -DBL_MAX; bnd[1] = DBL_MAX;
1077 BVLS_WRAPPER_FC(A.values(), m, n, m_grad_f.values(), bnd.values(),
1078 lambda.values(), res_norm, nsetp, w.values(),
1079 index.values(), ierr);
1080 if (ierr) {
1081 Cerr << "\nError: BVLS failed in constraint_penalty() in NonDGR"
1082 << std::endl;
1083 abort_handler(-1);
1084 }
1085
1086 lagrangeMult = lambda[0];
1087 return lagrangeMult*c_viol;
1088 #endif // DAKOTA_F90
1089 }
1090 else
1091 return 0.; // add NO_PENALTY to the enum instead?
1092 }
1093
1094
print_results(std::ostream & s,short results_state)1095 void NonDGlobalReliability::print_results(std::ostream& s, short results_state)
1096 {
1097 size_t i, j, wpp7 = write_precision + 7;
1098 const StringArray& fn_labels = iteratedModel.response_labels();
1099 s << "-----------------------------------------------------------------------"
1100 << "------";
1101
1102 print_densities(s);
1103
1104 // output CDF/CCDF level mappings (replaces NonD::print_level_mappings())
1105 s << std::scientific << std::setprecision(write_precision)
1106 << "\nLevel mappings for each response function:\n";
1107 for (i=0; i<numFunctions; i++) {
1108
1109 size_t num_levels = computedRespLevels[i].length();
1110 if (num_levels) {
1111 if (cdfFlag)
1112 s << "Cumulative Distribution Function (CDF) for ";
1113 else
1114 s << "Complementary Cumulative Distribution Function (CCDF) for ";
1115 s << fn_labels[i] << ":\n Response Level Probability Level "
1116 << "Reliability Index General Rel Index\n -------------- "
1117 << "----------------- ----------------- -----------------\n";
1118 for (j=0; j<num_levels; j++)
1119 s << " " << std::setw(wpp7) << computedRespLevels[i][j]
1120 << " " << std::setw(wpp7) << computedProbLevels[i][j]
1121 << std::setw(2*write_precision+18) << computedGenRelLevels[i][j]
1122 << '\n';
1123 }
1124 }
1125
1126 s << "-----------------------------------------------------------------------"
1127 << "------" << std::endl;
1128 }
1129
1130 } // namespace Dakota
1131