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: EffGlobalMinimizer
10 //- Description: Implementation code for the EffGlobalMinimizer class
11 //- Owner: Barron J Bichon, Vanderbilt University
12 //- Checked by:
13 //- Version:
14
15 //- Edited by: Anh Tran in 2020 for parallelization
16
17
18 #include "EffGlobalMinimizer.hpp"
19 #include "dakota_system_defs.hpp"
20 #include "dakota_data_io.hpp"
21 #include "NonDLHSSampling.hpp"
22 #include "RecastModel.hpp"
23 #include "DataFitSurrModel.hpp"
24 #include "DakotaApproximation.hpp"
25 #include "ProblemDescDB.hpp"
26 #include "DakotaGraphics.hpp"
27 #ifdef HAVE_NCSU
28 #include "NCSUOptimizer.hpp"
29 #endif
30 #include "DakotaModel.hpp"
31 #include "DakotaResponse.hpp"
32 #include "NormalRandomVariable.hpp"
33
34 //#define DEBUG
35 //#define DEBUG_PLOTS
36
37 namespace Dakota {
38
39 EffGlobalMinimizer* EffGlobalMinimizer::effGlobalInstance(NULL);
40
41
42 // This constructor accepts a Model
EffGlobalMinimizer(ProblemDescDB & problem_db,Model & model)43 EffGlobalMinimizer::EffGlobalMinimizer(ProblemDescDB& problem_db, Model& model):
44 SurrBasedMinimizer(problem_db, model, std::shared_ptr<TraitsBase>(new EffGlobalTraits())),
45 BatchSizeAcquisition(probDescDB.get_int("method.batch_size")),
46 BatchSizeExploration(probDescDB.get_int("method.batch_size.exploration")),
47 setUpType("model"), dataOrder(1)
48 {
49 // historical default convergence tolerances
50 if (convergenceTol < 0.0) convergenceTol = 1.0e-12;
51 distanceTol = probDescDB.get_real("method.x_conv_tol");
52 if (distanceTol < 0.0) distanceTol = 1.0e-8;
53
54 bestVariablesArray.push_back(iteratedModel.current_variables().copy());
55
56 // initialize augmented Lagrange multipliers
57 size_t num_multipliers = numNonlinearEqConstraints;
58 for (size_t i=0; i<numNonlinearIneqConstraints; i++) {
59 if (origNonlinIneqLowerBnds[i] > -bigRealBoundSize) // g has a lower bound
60 ++num_multipliers;
61 if (origNonlinIneqUpperBnds[i] < bigRealBoundSize) // g has an upper bound
62 ++num_multipliers;
63 }
64 augLagrangeMult.resize(num_multipliers);
65 augLagrangeMult = 0.;
66
67 truthFnStar.resize(numFunctions);
68
69 // Always build a global Gaussian process model. No correction is needed.
70 String approx_type = "global_kriging";
71 if (probDescDB.get_short("method.nond.emulator") == GP_EMULATOR)
72 approx_type = "global_gaussian";
73 else if (probDescDB.get_short("method.nond.emulator") == EXPGP_EMULATOR)
74 approx_type = "global_exp_gauss_proc";
75
76 String sample_reuse = "none"; // *** TO DO: allow reuse separate from import
77 UShortArray approx_order; // empty
78 short corr_order = -1, corr_type = NO_CORRECTION;
79 if (probDescDB.get_bool("method.derivative_usage")) {
80 if (approx_type == "global_gaussian") {
81 Cerr << "\nError: efficient_global does not support gaussian_process "
82 << "when derivatives present; use kriging instead." << std::endl;
83 abort_handler(METHOD_ERROR);
84 }
85 if (iteratedModel.gradient_type() != "none") dataOrder |= 2;
86 if (iteratedModel.hessian_type() != "none") dataOrder |= 4;
87 }
88 int db_samples = probDescDB.get_int("method.samples");
89 int samples = (db_samples > 0) ? db_samples :
90 (numContinuousVars+1)*(numContinuousVars+2)/2;
91 int lhs_seed = probDescDB.get_int("method.random_seed");
92 unsigned short sample_type = SUBMETHOD_DEFAULT;
93 String rng; // empty string: use default
94 //int symbols = samples; // symbols needed for DDACE
95 bool vary_pattern = false;// for consistency across any outer loop invocations
96 // get point samples file
97 const String& import_pts_file
98 = probDescDB.get_string("method.import_build_points_file");
99 if (!import_pts_file.empty()) // *** TO DO: allow reuse separate from import
100 { samples = 0; sample_reuse = "all"; }
101
102 Iterator dace_iterator;
103 // The following uses on the fly derived ctor:
104 dace_iterator.assign_rep(std::make_shared<NonDLHSSampling>(iteratedModel, sample_type,
105 samples, lhs_seed, rng, vary_pattern, ACTIVE_UNIFORM));
106 // only use derivatives if the user requested and they are available
107 dace_iterator.active_set_request_values(dataOrder);
108
109 // Construct f-hat (fHatModel) using a GP approximation for each response function over
110 // the active/design vars (same view as iteratedModel: not the typical All
111 // view for DACE).
112 //const Variables& curr_vars = iteratedModel.current_variables();
113 ActiveSet gp_set = iteratedModel.current_response().active_set(); // copy
114 gp_set.request_values(1); // no surr deriv evals, but GP may be grad-enhanced
115 fHatModel.assign_rep(std::make_shared<DataFitSurrModel>
116 (dace_iterator, iteratedModel,
117 gp_set, approx_type, approx_order, corr_type, corr_order, dataOrder,
118 outputLevel, sample_reuse, import_pts_file,
119 probDescDB.get_ushort("method.import_build_format"),
120 probDescDB.get_bool("method.import_build_active_only"),
121 probDescDB.get_string("method.export_approx_points_file"),
122 probDescDB.get_ushort("method.export_approx_format")));
123
124 // Following this ctor, IteratorScheduler::init_iterator() initializes the
125 // parallel configuration for EffGlobalMinimizer + iteratedModel using
126 // EffGlobalMinimizer's maxEvalConcurrency. During fHatModel construction
127 // above, DataFitSurrModel::derived_init_communicators() initializes the
128 // parallel config for dace_iterator + iteratedModel using dace_iterator's
129 // maxEvalConcurrency. The only iteratedModel concurrency currently exercised
130 // is that used by dace_iterator within the initial GP construction, but the
131 // EffGlobalMinimizer maxEvalConcurrency must still be set so as to avoid
132 // parallel config errors resulting from avail_procs > max_concurrency within
133 // IteratorScheduler::init_iterator(). A max of the local derivative
134 // concurrency and the DACE concurrency is used for this purpose.
135 maxEvalConcurrency = std::max(maxEvalConcurrency, dace_iterator.maximum_evaluation_concurrency());
136
137 // Configure a RecastModel with one objective and no constraints using the
138 // alternate minimalist constructor: the recast fn pointers are reset for
139 // each level within the run fn.
140 SizetArray recast_vars_comps_total; // default: empty; no change in size
141 BitArray all_relax_di, all_relax_dr; // default: empty; no discrete relaxation
142 short recast_resp_order = 1; // nongradient-based optimizers
143 eifModel.assign_rep(std::make_shared<RecastModel>
144 (fHatModel, recast_vars_comps_total, all_relax_di,
145 all_relax_dr, 1, 0, 0, recast_resp_order));
146
147 // must use alternate NoDB ctor chain
148 int max_iterations = 10000, max_fn_evals = 50000;
149 double min_box_size = 1.e-15, vol_box_size = 1.e-15;
150 #ifdef HAVE_NCSU
151 approxSubProbMinimizer.assign_rep(std::make_shared<NCSUOptimizer>(eifModel, max_iterations,
152 max_fn_evals, min_box_size, vol_box_size));
153 #else
154 Cerr << "NCSU DIRECT is not available to optimize the GP subproblems. "
155 << "Aborting process." << std::endl;
156 abort_handler(METHOD_ERROR);
157 #endif //HAVE_NCSU
158 }
159
160
161
~EffGlobalMinimizer()162 EffGlobalMinimizer::~EffGlobalMinimizer() {}
163
164
core_run()165 void EffGlobalMinimizer::core_run()
166 {
167 if (setUpType=="model")
168 minimize_surrogates_on_model();
169 else if (setUpType=="user_functions") {
170 //minimize_surrogates_on_user_functions();
171 Cerr << "Error: bad setUpType in EffGlobalMinimizer::core_run()."
172 << std::endl;
173 abort_handler(METHOD_ERROR);
174 }
175 else {
176 Cerr << "Error: bad setUpType in EffGlobalMinimizer::core_run()."
177 << std::endl;
178 abort_handler(METHOD_ERROR);
179 }
180 }
181
182
minimize_surrogates_on_model()183 void EffGlobalMinimizer::minimize_surrogates_on_model()
184 {
185 //Parallel EGO: --------------------------
186 // Solve the problem.
187 //Parallel EGO: --------------------------
188 EffGlobalMinimizer* prev_instance = effGlobalInstance;
189 // build initial GP once for all response functions: fHatModel.build_approximation()
190 initialize();
191
192 // Iterate until EGO converges
193 unsigned short eif_convergence_cntr = 0, dist_convergence_cntr = 0,
194 eif_convergence_limit = 2, dist_convergence_limit = 1;
195 globalIterCount = 0;
196 bool approx_converged = false;
197
198 // Decided for now (10-25-2013) to have EGO take the maxIterations
199 // as the default from minimizer, so it will be initialized as 100
200 // maxIterations = 25*numContinuousVars;
201 RealVector prev_cv_star;
202
203 while (!approx_converged) {
204
205 // Add safeguard: If model cannot run asynchronously
206 // then reset BatchSizeAcquisition = 1
207 // and throw warnings on screens
208 // and proceed with batchSize = 1
209
210 bool parallel_flag = false; // default: run sequential by default
211 // bool parallel_flag = true; // debug -- always run parallel
212
213 if (BatchSizeAcquisition > 1 || BatchSizeExploration > 1) {
214 if (iteratedModel.asynch_flag()) // change model.asynch_flag() to iteratedModel.asynch_flag()
215 parallel_flag = true; // turn on parallel_flag if the requirements are satisfied
216 else {
217 Cerr << "Warning: concurrent operations not supported by model. Batch size request ignored." << std::endl;
218 BatchSizeAcquisition = 1; // reverse to the default sequential
219 BatchSizeExploration = 0; // reverse to the default sequential
220 }
221 }
222
223 if (parallel_flag) { // begin if parallel_flag = true -- then run in parallel
224
225 ++globalIterCount;
226
227 // Reset the convergence counters
228 dist_convergence_cntr = 0; // reset distance convergence counters
229 dist_convergence_limit = BatchSizeAcquisition; // set convergence limit for parallel EGO
230
231 // Initialize the input array for the batch
232 VariablesArray input_array_batch_acquisition(BatchSizeAcquisition);
233
234 // Note: vars_star: input; resp_star: output (liar); resp_star_truth: output (true)
235 size_t numDataPts; // debug
236
237 // Construct the batch
238 Cout << "\n>>>>> Initiating global optimization\n";
239
240 for (int i_batch_acquisition = 0; i_batch_acquisition < BatchSizeAcquisition; i_batch_acquisition++) {
241
242 // Initialize EIF recast model
243 Sizet2DArray vars_map, primary_resp_map(1), secondary_resp_map;
244 BoolDequeArray nonlinear_resp_map(1, BoolDeque(numFunctions, true));
245 primary_resp_map[0].resize(numFunctions);
246 for (size_t i=0; i<numFunctions; i++)
247 primary_resp_map[0][i] = i;
248 std::shared_ptr<RecastModel> eif_model_rep =
249 std::static_pointer_cast<RecastModel>(eifModel.model_rep());
250 eif_model_rep->init_maps(vars_map, false, NULL, NULL,
251 primary_resp_map, secondary_resp_map, nonlinear_resp_map,
252 EIF_objective_eval, NULL);
253
254 // Determine fnStar from among sample data
255 get_best_sample();
256 bestVariablesArray.front().continuous_variables(varStar); // debug
257 bestResponseArray.front().function_values(truthFnStar); // debug
258
259 // Execute GLOBAL search and retrieve results
260 // Cout << "\n>>>>> Initiating global optimization\n";
261 ParLevLIter pl_iter = methodPCIter->mi_parallel_level_iterator(miPLIndex);
262 approxSubProbMinimizer.run(pl_iter); // maximize the EI acquisition fucntion
263 const Variables& vars_star = approxSubProbMinimizer.variables_results();
264 const RealVector& c_vars = vars_star.continuous_variables();
265 const Response& resp_star = approxSubProbMinimizer.response_results();
266 const Real& eif_star = resp_star.function_value(0);
267
268 // Cout << "c_vars(i_batch_acquisition = " << i_batch_acquisition << ") = " << c_vars << ".\n"; // debug
269
270 // Get expected value for output
271 fHatModel.continuous_variables(c_vars);
272 fHatModel.evaluate();
273 const Response& approx_response = fHatModel.current_response();
274
275 Real aug_lag = augmented_lagrangian_merit(approx_response.function_values(),
276 iteratedModel.primary_response_fn_sense(),
277 iteratedModel.primary_response_fn_weights(),
278 origNonlinIneqLowerBnds,
279 origNonlinIneqUpperBnds,
280 origNonlinEqTargets);
281
282 Cout << "\nResults of EGO iteration:\nFinal point =\n"
283 << c_vars << "Expected Improvement =\n "
284 << std::setw(write_precision+7) << -eif_star
285 << "\n " << std::setw(write_precision+7)
286 << aug_lag << " [merit]\n";
287
288 // Impose constant liar -- temporarily cast constant liar as observations
289 // const IntResponsePair resp_star_liar(iteratedModel.evaluation_id(), approx_response);
290 const IntResponsePair resp_star_liar(iteratedModel.evaluation_id() + i_batch_acquisition + 1, approx_response); // implement a liar counter
291
292 // Update GP
293 // debug
294 numDataPts = fHatModel.approximation_data(0).points(); // debug
295 Cout << "\nParallel EGO: Adding liar response...\n"; // debug
296 // Cout << "globalIterCount = " << globalIterCount << ".\n"; // debug
297 // Cout << "i_batch_acquisition = " << i_batch_acquisition << ".\n"; // debug
298
299 // Append constant liar to fHatModel (aka heuristic liar)
300 if (i_batch_acquisition < BatchSizeAcquisition - 1)
301 fHatModel.append_approximation(vars_star, resp_star_liar, true); // debug
302
303 // Cout << "fHatModel # of points = " << numDataPts << ".\n"; // debug
304 // Update constraints based on the constant liar
305 if (numNonlinearConstraints) {
306 const RealVector& fns_star_liar = resp_star_liar.second.function_values();
307 Real norm_cv_star = std::sqrt(constraint_violation(fns_star_liar, 0.));
308 if (norm_cv_star < etaSequence)
309 update_augmented_lagrange_multipliers(fns_star_liar);
310 else
311 update_penalty();
312 }
313 Cout << "Parallel EGO: Finished adding liar responses!\n"; // debug
314 // const RealVector& variances = fHatModel.approximation_variances(vars_star); // debug
315 // Cout << "debug: vars_star = " << vars_star; // debug
316 // Cout << "debug: variances(vars_star) = " << variances; // debug
317 // Copy vars_star (input) to the batch before querying
318 input_array_batch_acquisition[i_batch_acquisition] = vars_star.copy();
319 }
320
321 // Delete liar responses
322 for (int i_batch_acquisition = 0; i_batch_acquisition < BatchSizeAcquisition; i_batch_acquisition++) {
323 if (i_batch_acquisition < BatchSizeAcquisition - 1)
324 fHatModel.pop_approximation(false);
325
326 // debug
327 numDataPts = fHatModel.approximation_data(0).points(); // debug
328 Cout << "\nParallel EGO: Deleting liar response...\n"; // debug
329 // Cout << "globalIterCount = " << globalIterCount << ".\n"; // debug
330 // Cout << "i_batch_acquisition = " << i_batch_acquisition << ".\n"; // debug
331 // Cout << "fHatModel # of points = " << numDataPts << ".\n"; // debug
332 }
333 Cout << "\nParallel EGO: Finished deleting liar responses!\n"; // debug
334
335 // Query the batch
336 for (int i_batch_acquisition = 0; i_batch_acquisition < BatchSizeAcquisition; i_batch_acquisition++) {
337 fHatModel.component_parallel_mode(TRUTH_MODEL_MODE);
338 iteratedModel.active_variables(input_array_batch_acquisition[i_batch_acquisition]);
339 ActiveSet set = iteratedModel.current_response().active_set();
340 set.request_values(dataOrder);
341 iteratedModel.evaluate_nowait(set);
342 }
343
344 // Get true responses resp_star_truth
345 const IntResponseMap resp_star_truth = iteratedModel.synchronize();
346
347 // Update the GP approximation with batch results
348 Cout << "\nParallel EGO: Adding true response...\n"; // debug
349 fHatModel.append_approximation(input_array_batch_acquisition, resp_star_truth, true);
350 Cout << "\nParallel EGO: Finished adding true responses!\n"; // debug
351
352 // debug
353 numDataPts = fHatModel.approximation_data(0).points(); // debug
354 // Cout << "globalIterCount = " << globalIterCount << ".\n"; // debug
355 // Cout << "fHatModel # of points = " << numDataPts << ".\n"; // debug
356
357 // Update constraints
358 for (int i_batch_acquisition = 0; i_batch_acquisition < BatchSizeAcquisition; i_batch_acquisition++) {
359 const Variables& vars_star = approxSubProbMinimizer.variables_results();
360 const RealVector& c_vars = vars_star.continuous_variables();
361 const Response& resp_star = approxSubProbMinimizer.response_results();
362 const Real& eif_star = resp_star.function_value(0);
363
364 debug_print_values();
365 // check_convergence(eif_star, c_vars, prev_cv_star, eif_convergence_cntr, dist_convergence_cntr);
366 // Check for convergence based on max EIF
367 if ( -eif_star < convergenceTol )
368 ++eif_convergence_cntr;
369
370 // Check for convergence based in distance between successive points.
371 // If the dist between successive points is very small, then there is
372 // little value in updating the GP since the new training point will
373 // essentially be the previous optimal point.
374
375 Real distCStar = (prev_cv_star.empty()) ? DBL_MAX : rel_change_L2(c_vars, prev_cv_star);
376 // update prev_cv_star
377 copy_data(c_vars, prev_cv_star);
378
379 if (distCStar < distanceTol)
380 ++dist_convergence_cntr;
381
382 // If DIRECT failed to find a point with EIF>0, it returns the
383 // center point as the optimal solution. EGO may have converged,
384 // but DIRECT may have just failed to find a point with a good
385 // EIF value. Adding this midpoint can alter the GPs enough to
386 // to allow DIRECT to find something useful, so we force
387 // max(EIF)<tol twice to make sure. Note that we cannot make
388 // this check more than 2 because it would cause EGO to add
389 // the center point more than once, which will damage the GPs.
390 // Unfortunately, when it happens the second time, it may still
391 // be that DIRECT failed and not that EGO converged.
392 debug_print_counter(globalIterCount, eif_star, distCStar, dist_convergence_cntr);
393 }
394
395 // Check convergence
396 if ( dist_convergence_cntr >= dist_convergence_limit ||
397 eif_convergence_cntr >= eif_convergence_limit ||
398 globalIterCount >= maxIterations ) {
399
400 approx_converged = true;
401
402 // only verbose if debug
403 if (outputLevel > NORMAL_OUTPUT) {
404 // debug
405 if (dist_convergence_cntr >= dist_convergence_limit) {
406 Cout << "\nStopping criteria met: dist_convergence_cntr (="
407 << dist_convergence_cntr
408 << ") >= dist_convergence_limit (="
409 << dist_convergence_limit
410 << ").\n";
411 }
412 if (eif_convergence_cntr >= eif_convergence_limit) {
413 Cout << "\nStopping criteria met: eif_convergence_cntr (="
414 << eif_convergence_cntr
415 << ") >= eif_convergence_limit (="
416 << eif_convergence_limit
417 << ").\n";
418 }
419 if (globalIterCount >= maxIterations) {
420 Cout << "\nStopping criteria met: globalIterCount (="
421 << globalIterCount
422 << ") >= maxIterations (="
423 << maxIterations
424 << ").\n";
425 }
426 }
427 }
428 else {
429 approx_converged = false;
430
431 // only verbose if debug
432 if (outputLevel > NORMAL_OUTPUT) {
433 // debug
434 if (dist_convergence_cntr < dist_convergence_limit) {
435 Cout << "\nStopping criteria not met: dist_convergence_cntr (="
436 << dist_convergence_cntr
437 << ") < dist_convergence_limit (="
438 << dist_convergence_limit
439 << ").\n";
440 }
441 if (eif_convergence_cntr < eif_convergence_limit) {
442 Cout << "\nStopping criteria not met: eif_convergence_cntr (="
443 << eif_convergence_cntr
444 << ") < eif_convergence_limit (="
445 << eif_convergence_limit
446 << ").\n";
447 }
448 if (globalIterCount < maxIterations) {
449 Cout << "\nStopping criteria not met: globalIterCount (="
450 << globalIterCount
451 << ") < maxIterations (="
452 << maxIterations
453 << ").\n";
454 }
455 }
456 }
457
458
459 // Update constraints
460 if (numNonlinearConstraints) {
461 IntRespMCIter batch_response_it; // IntRMMIter (iterator)? IntRMMCIter (const_iterator)?
462 for (batch_response_it = resp_star_truth.begin(); batch_response_it != resp_star_truth.end(); batch_response_it++) {
463 // Update the merit function parameters
464 // Logic follows Conn, Gould, and Toint, section 14.4:
465 // const RealVector& fns_star_truth = resp_star_truth.second.function_values(); // old implementation
466 const RealVector& fns_star_truth = batch_response_it->second.function_values(); // current implementation
467 Real norm_cv_star = std::sqrt(constraint_violation(fns_star_truth, 0.));
468 if (norm_cv_star < etaSequence)
469 update_augmented_lagrange_multipliers(fns_star_truth);
470 else
471 update_penalty();
472 }
473 }
474 } // end if parallel_flag = true -- then run in parallel
475 else { // else begin parallel_flag = false -- then run sequentially (reinstate old sequential implementation)
476
477 ++globalIterCount;
478
479 // Initialize EIF recast model
480 Sizet2DArray vars_map, primary_resp_map(1), secondary_resp_map;
481 BoolDequeArray nonlinear_resp_map(1, BoolDeque(numFunctions, true));
482 primary_resp_map[0].resize(numFunctions);
483 for (size_t i=0; i<numFunctions; i++)
484 primary_resp_map[0][i] = i;
485 std::shared_ptr<RecastModel> eif_model_rep =
486 std::static_pointer_cast<RecastModel>(eifModel.model_rep());
487 eif_model_rep->init_maps(vars_map, false, NULL, NULL,
488 primary_resp_map, secondary_resp_map, nonlinear_resp_map,
489 EIF_objective_eval, NULL);
490
491 // Determine fnStar from among sample data
492 get_best_sample();
493
494 // Execute GLOBAL search and retrieve results
495 Cout << "\n>>>>> Initiating global optimization\n";
496 ParLevLIter pl_iter = methodPCIter->mi_parallel_level_iterator(miPLIndex);
497 approxSubProbMinimizer.run(pl_iter); // maximize the EI acquisition fucntion
498 const Variables& vars_star = approxSubProbMinimizer.variables_results();
499 const RealVector& c_vars = vars_star.continuous_variables();
500 const Response& resp_star = approxSubProbMinimizer.response_results();
501 const Real& eif_star = resp_star.function_value(0);
502
503 // Get expected value for output
504 fHatModel.continuous_variables(c_vars);
505 fHatModel.evaluate();
506 const Response& approx_response = fHatModel.current_response(); // .function_values();
507
508 Real aug_lag = augmented_lagrangian_merit(approx_response.function_values(),
509 iteratedModel.primary_response_fn_sense(),
510 iteratedModel.primary_response_fn_weights(),
511 origNonlinIneqLowerBnds,
512 origNonlinIneqUpperBnds,
513 origNonlinEqTargets);
514
515 Cout << "\nResults of EGO iteration:\nFinal point =\n"
516 << c_vars << "Expected Improvement =\n "
517 << std::setw(write_precision+7) << -eif_star
518 << "\n " << std::setw(write_precision+7)
519 << aug_lag << " [merit]\n";
520
521 debug_print_values();
522 // check_convergence(eif_star, c_vars, prev_cv_star, eif_convergence_cntr, dist_convergence_cntr);
523 // Check for convergence based on max EIF
524 if ( -eif_star < convergenceTol )
525 ++eif_convergence_cntr;
526
527 // Check for convergence based in distance between successive points.
528 // If the dist between successive points is very small, then there is
529 // little value in updating the GP since the new training point will
530 // essentially be the previous optimal point.
531
532 Real distCStar = (prev_cv_star.empty()) ? DBL_MAX : rel_change_L2(c_vars, prev_cv_star);
533 // update prev_cv_star
534 copy_data(c_vars, prev_cv_star);
535 if (distCStar < distanceTol)
536 ++dist_convergence_cntr;
537
538 // If DIRECT failed to find a point with EIF>0, it returns the
539 // center point as the optimal solution. EGO may have converged,
540 // but DIRECT may have just failed to find a point with a good
541 // EIF value. Adding this midpoint can alter the GPs enough to
542 // to allow DIRECT to find something useful, so we force
543 // max(EIF)<tol twice to make sure. Note that we cannot make
544 // this check more than 2 because it would cause EGO to add
545 // the center point more than once, which will damage the GPs.
546 // Unfortunately, when it happens the second time, it may still
547 // be that DIRECT failed and not that EGO converged.
548 debug_print_counter(globalIterCount, eif_star, distCStar, dist_convergence_cntr);
549
550 // Evaluate response_star_truth
551 fHatModel.component_parallel_mode(TRUTH_MODEL_MODE);
552 iteratedModel.continuous_variables(c_vars);
553 ActiveSet set = iteratedModel.current_response().active_set();
554 set.request_values(dataOrder);
555 iteratedModel.evaluate(set);
556 IntResponsePair resp_star_truth(iteratedModel.evaluation_id(), iteratedModel.current_response());
557
558 // Update the GP approximation
559 fHatModel.append_approximation(vars_star, resp_star_truth, true);
560
561 // Check convergence
562 if ( dist_convergence_cntr >= dist_convergence_limit ||
563 eif_convergence_cntr >= eif_convergence_limit ||
564 globalIterCount >= maxIterations )
565 approx_converged = true;
566 else {
567 // Update constraints
568 if (numNonlinearConstraints) {
569 // Update the merit function parameters
570 // Logic follows Conn, Gould, and Toint, section 14.4:
571 const RealVector& fns_star_truth = resp_star_truth.second.function_values();
572 Real norm_cv_star = std::sqrt(constraint_violation(fns_star_truth, 0.));
573 if (norm_cv_star < etaSequence)
574 update_augmented_lagrange_multipliers(fns_star_truth);
575 else
576 update_penalty();
577 }
578 }
579 } // end if parallel_flag = false -- then run sequentially
580 } // end approx convergence while loop
581
582
583 // (conditionally) export final surrogates
584 export_final_surrogates(fHatModel);
585
586 // Set best variables and response for use by strategy level.
587 // c_vars, fmin contain the optimal design
588 get_best_sample(); // pull optimal result from sample data
589 bestVariablesArray.front().continuous_variables(varStar);
590 bestResponseArray.front().function_values(truthFnStar);
591
592 // restore in case of recursion
593 effGlobalInstance = prev_instance;
594
595 debug_plots();
596
597 }
598
599
600 /** To maximize expected improvement, the approxSubProbMinimizer will minimize -(expected_improvement). */
EIF_objective_eval(const Variables & sub_model_vars,const Variables & recast_vars,const Response & sub_model_response,Response & recast_response)601 void EffGlobalMinimizer::EIF_objective_eval(const Variables& sub_model_vars,
602 const Variables& recast_vars,
603 const Response& sub_model_response,
604 Response& recast_response) {
605 // Means are passed in, but must retrieve variance from the GP
606 const RealVector& means = sub_model_response.function_values();
607 const RealVector& variances = effGlobalInstance->fHatModel.approximation_variances(recast_vars);
608 const ShortArray& recast_asv = recast_response.active_set_request_vector();
609
610 if (recast_asv[0] & 1) { // return -EI since we are maximizing
611 Real neg_ei = -effGlobalInstance->expected_improvement(means, variances);
612 recast_response.function_value(neg_ei, 0);
613 }
614 }
615
616
expected_improvement(const RealVector & means,const RealVector & variances)617 Real EffGlobalMinimizer::expected_improvement(const RealVector& means, const RealVector& variances)
618 {
619 // Objective calculation will incorporate any sense changes or
620 // weights, such that this is an objective to minimize.
621 Real mean = objective(means, iteratedModel.primary_response_fn_sense(),
622 iteratedModel.primary_response_fn_weights()), stdv;
623
624 if ( numNonlinearConstraints ) {
625 // mean_M = mean_f + lambda*EV + r_p*EV*EV
626 // stdv_M = stdv_f
627 const RealVector& ev = expected_violation(means, variances);
628 for (size_t i=0; i<numNonlinearConstraints; ++i)
629 mean += augLagrangeMult[i]*ev[i] + penaltyParameter*ev[i]*ev[i];
630 stdv = std::sqrt(variances[0]); // ***
631 }
632 else { // extend for NLS/MOO ***
633 // mean_M = M(mu_f)
634 // stdv_M = sqrt(var_f)
635 stdv = std::sqrt(variances[0]); // *** sqrt(sum(variances(1:nUsrPrimaryFns))
636 }
637
638 // Calculate the expected improvement
639 Real cdf, pdf;
640 Real snv = (meritFnStar - mean); // standard normal variate
641 if(std::fabs(snv) >= std::fabs(stdv)*50.0) {
642 //this will trap the denominator=0.0 case even if numerator=0.0
643 pdf=0.0;
644 cdf=(snv>0.0)?1.0:0.0;
645 }
646 else{
647 snv/=stdv;
648 cdf = Pecos::NormalRandomVariable::std_cdf(snv);
649 pdf = Pecos::NormalRandomVariable::std_pdf(snv);
650 }
651
652 Real ei = (meritFnStar - mean) * cdf + stdv * pdf;
653
654 return ei;
655 }
656
657
expected_violation(const RealVector & means,const RealVector & variances)658 RealVector EffGlobalMinimizer::expected_violation(const RealVector& means, const RealVector& variances)
659 {
660 RealVector ev(numNonlinearConstraints);
661
662 size_t i, cntr=0;
663 // inequality constraints
664 for (i=0; i<numNonlinearIneqConstraints; i++) {
665 const Real& mean = means[numUserPrimaryFns+i];
666 const Real& stdv = std::sqrt(variances[numUserPrimaryFns+i]);
667 const Real& lbnd = origNonlinIneqLowerBnds[i];
668 const Real& ubnd = origNonlinIneqUpperBnds[i];
669 if (lbnd > -bigRealBoundSize) {
670 Real cdf, pdf;
671 Real snv = (lbnd-mean);
672 if(std::fabs(snv)>=std::fabs(stdv)*50.0){
673 pdf=0.0;
674 cdf=(snv>0.0)?1.0:0.0;
675 }
676 else {
677 snv/=stdv; //now snv is the standard normal variate
678 cdf = Pecos::NormalRandomVariable::std_cdf(snv);
679 pdf = Pecos::NormalRandomVariable::std_pdf(snv);
680 }
681 ev[cntr++] = (lbnd-mean)*cdf + stdv*pdf;
682 }
683 if (ubnd < bigRealBoundSize) {
684 Real cdf, pdf;
685 Real snv = (ubnd-mean);
686 if(std::fabs(snv)>=std::fabs(stdv)*50.0){
687 pdf=0.0;
688 cdf=(snv>0.0)?1.0:0.0;
689 }
690 else{
691 snv/=stdv;
692 cdf = Pecos::NormalRandomVariable::std_cdf(snv);
693 pdf = Pecos::NormalRandomVariable::std_pdf(snv);
694 }
695 ev[cntr++] = (mean-ubnd)*(1.-cdf) + stdv*pdf;
696 }
697 }
698 // equality constraints
699 for (i=0; i<numNonlinearEqConstraints; i++) {
700 const Real& mean = means[numUserPrimaryFns+numNonlinearIneqConstraints+i];
701 const Real& stdv = std::sqrt(variances[numUserPrimaryFns+numNonlinearIneqConstraints+i]);
702 const Real& zbar = origNonlinEqTargets[i];
703 Real cdf, pdf;
704 Real snv = (zbar-mean);
705 if(std::fabs(snv)*50.0>=std::fabs(stdv)) {
706 pdf=0.0;
707 cdf=(snv>=0.0)?1.0:0.0;
708 }
709 else{
710 snv/=stdv;
711 cdf = Pecos::NormalRandomVariable::std_cdf(snv);
712 pdf = Pecos::NormalRandomVariable::std_pdf(snv);
713 }
714 ev[cntr++] = (zbar-mean)*(2.*cdf-1.) + 2.*stdv*pdf;
715 }
716
717 return ev;
718 }
719
720
get_best_sample()721 void EffGlobalMinimizer::get_best_sample()
722 {
723 // Pull the samples and responses from data used to build latest GP
724 // to determine fnStar for use in the expected improvement function
725
726 const Pecos::SurrogateData& gp_data_0 = fHatModel.approximation_data(0);
727 const Pecos::SDVArray& sdv_array = gp_data_0.variables_data();
728 const Pecos::SDRArray& sdr_array = gp_data_0.response_data();
729
730 size_t i, sam_star_idx = 0, num_data_pts = gp_data_0.points();
731 Real fn, fn_star = DBL_MAX;
732
733 for (i=0; i<num_data_pts; ++i) {
734 // Cout << "\nget_best_sample(): num_data_pts = " << num_data_pts << "\n" ; // debug
735 const RealVector& sams = sdv_array[i].continuous_variables();
736
737 fHatModel.continuous_variables(sams);
738 fHatModel.evaluate();
739 const RealVector& f_hat = fHatModel.current_response().function_values();
740 fn = augmented_lagrangian_merit(f_hat,
741 iteratedModel.primary_response_fn_sense(),
742 iteratedModel.primary_response_fn_weights(),
743 origNonlinIneqLowerBnds,
744 origNonlinIneqUpperBnds,
745 origNonlinEqTargets);
746
747 if (fn < fn_star) {
748 copy_data(sams, varStar);
749 sam_star_idx = i;
750 fn_star = fn;
751 meritFnStar = fn;
752 truthFnStar[0] = sdr_array[i].response_function();
753 }
754 }
755
756 // update truthFnStar with all additional primary/secondary fns corresponding
757 // to lowest merit function value
758 for (i=1; i<numFunctions; ++i) {
759 const Pecos::SDRArray& sdr_array = fHatModel.approximation_data(i).response_data();
760 truthFnStar[i] = sdr_array[sam_star_idx].response_function();
761 }
762 }
763
764
update_penalty()765 void EffGlobalMinimizer::update_penalty()
766 {
767 // Logic follows Conn, Gould, and Toint, section 14.4, step 3
768 // CGT use mu *= tau with tau = 0.01 -> r_p *= 50
769 // Rodriguez, Renaud, Watson: r_p *= 10
770 // Robinson, Willcox, Eldred, and Haimes: r_p *= 5
771 penaltyParameter *= 10.;
772 //penaltyParameter = std::min(penaltyParameter, 1.e+20); // cap the max penalty?
773 Real mu = 1./2./penaltyParameter; // conversion between r_p and mu penalties
774 etaSequence = eta * std::pow(mu, alphaEta);
775
776 #ifdef DEBUG
777 Cout << "Penalty updated: " << penaltyParameter << '\n'
778 << "eta updated: " << etaSequence << '\n'
779 << "Augmented Lagrange multipliers:\n" << augLagrangeMult;
780 #endif
781 }
782
783
784 // This override exists purely to prevent an optimizer/minimizer from declaring sources
785 // when it's being used to evaluate a user-defined function (e.g. finding the correlation
786 // lengths of Dakota's GP).
declare_sources()787 void EffGlobalMinimizer::declare_sources() {
788 if(setUpType == "user_functions")
789 return;
790 else
791 Iterator::declare_sources();
792 }
793
794
get_augmented_lagrangian(const RealVector & mean,const RealVector & c_vars,const Real & eif_star)795 Real EffGlobalMinimizer::get_augmented_lagrangian(const RealVector& mean,
796 const RealVector& c_vars,
797 const Real& eif_star) {
798 Real aug_lag = augmented_lagrangian_merit(mean,
799 iteratedModel.primary_response_fn_sense(),
800 iteratedModel.primary_response_fn_weights(), origNonlinIneqLowerBnds,
801 origNonlinIneqUpperBnds, origNonlinEqTargets);
802
803 // print out message
804 Cout << "\nResults of EGO iteration:\nFinal point =\n" << c_vars
805 << "Expected Improvement =\n "
806 << std::setw(write_precision+7) << -eif_star
807 << "\n " << std::setw(write_precision+7)
808 << aug_lag << " [merit]\n";
809 return aug_lag;
810 }
811
812
check_convergence(const Real & eif_star,const RealVector & c_vars,RealVector prev_cv_star,unsigned short eif_convergence_cntr,unsigned short dist_convergence_cntr)813 void EffGlobalMinimizer::check_convergence(const Real& eif_star,
814 const RealVector& c_vars,
815 RealVector prev_cv_star,
816 unsigned short eif_convergence_cntr,
817 unsigned short dist_convergence_cntr) {
818 // Check for convergence based on max EIF
819 if ( -eif_star < convergenceTol )
820 ++eif_convergence_cntr;
821
822 // Check for convergence based in distance between successive points.
823 // If the dist between successive points is very small, then there is
824 // little value in updating the GP since the new training point will
825 // essentially be the previous optimal point.
826
827 Real distCStar = (prev_cv_star.empty()) ? DBL_MAX : rel_change_L2(c_vars, prev_cv_star);
828 // update prev_cv_star
829 copy_data(c_vars, prev_cv_star);
830 if (distCStar < distanceTol)
831 ++dist_convergence_cntr;
832
833 // If DIRECT failed to find a point with EIF>0, it returns the
834 // center point as the optimal solution. EGO may have converged,
835 // but DIRECT may have just failed to find a point with a good
836 // EIF value. Adding this midpoint can alter the GPs enough to
837 // to allow DIRECT to find something useful, so we force
838 // max(EIF)<tol twice to make sure. Note that we cannot make
839 // this check more than 2 because it would cause EGO to add
840 // the center point more than once, which will damage the GPs.
841 // Unfortunately, when it happens the second time, it may still
842 // be that DIRECT failed and not that EGO converged.
843
844 return;
845
846 }
847
initialize()848 void EffGlobalMinimizer::initialize() {
849 // set the object instance pointers for use within the static member fns
850 effGlobalInstance = this;
851
852 // now that variables/labels/bounds/targets have flowed down at run-time from
853 // any higher level recursions, propagate them up the instantiate-on-the-fly
854 // Model recursion so that they are correct when they propagate back down.
855 eifModel.update_from_subordinate_model(); // depth = max
856
857 // (We might want a more selective update from submodel, or make a
858 // new EIFModel specialization of RecastModel.) Always want to
859 // minimize the negative expected improvement as posed in the
860 // eifModel, which consumes min/max sense and weights, and recasts
861 // nonlinear constraints, so we don't let these propagate to the
862 // approxSubproblemMinimizer.
863 eifModel.primary_response_fn_sense(BoolDeque());
864 eifModel.primary_response_fn_weights(RealVector(), false); // no recursion
865 eifModel.reshape_constraints(0,0, eifModel.num_linear_ineq_constraints(),
866 eifModel.num_linear_eq_constraints());
867
868 // Build initial GP once for all response functions
869 fHatModel.build_approximation();
870
871 return;
872 }
873
debug_print_values()874 void EffGlobalMinimizer::debug_print_values() {
875 #ifdef DEBUG
876 RealVector variance = fHatModel.approximation_variances(vars_star);
877 RealVector ev = expected_violation(mean,variance);
878 RealVector stdv(numFunctions);
879 for (size_t i=0; i<numFunctions; i++)
880 stdv[i] = std::sqrt(variance[i]);
881 Cout << "\nexpected values =\n" << mean
882 << "\nstandard deviation =\n" << stdv
883 << "\nexpected violation =\n" << ev << std::endl;
884 #endif //DEBUG
885 }
886
debug_print_counter(unsigned short globalIterCount,const Real & eif_star,Real distCStar,unsigned short dist_convergence_cntr)887 void EffGlobalMinimizer::debug_print_counter(unsigned short globalIterCount,
888 const Real& eif_star,
889 Real distCStar,
890 unsigned short dist_convergence_cntr) {
891 #ifdef DEBUG
892 Cout << "EGO Iteration " << globalIterCount << "\neif_star " << eif_star
893 << "\ndistCStar " << distCStar << "\ndist_convergence_cntr "
894 << dist_convergence_cntr << '\n';
895 #endif //DEBUG
896 }
897
debug_plots()898 void EffGlobalMinimizer::debug_plots() {
899
900 #ifdef DEBUG_PLOTS
901 // DEBUG - output set of samples used to build the GP
902 // If problem is 2d, output a grid of points on the GP
903 // and truth (if requested)
904 for (size_t i=0; i<numFunctions; i++) {
905 std::string samsfile("ego_sams");
906 std::string tag = "_" + std::to_string(i+1) + ".out";
907 samsfile += tag;
908 std::ofstream samsOut(samsfile.c_str(),std::ios::out);
909 samsOut << std::scientific;
910 const Pecos::SurrogateData& gp_data = fHatModel.approximation_data(i);
911 const Pecos::SDVArray& sdv_array = gp_data.variables_data();
912 const Pecos::SDRArray& sdr_array = gp_data.response_data();
913 size_t num_data_pts = gp_data.size(), num_vars = fHatModel.cv();
914 for (size_t j=0; j<num_data_pts; ++j) {
915 samsOut << '\n';
916 const RealVector& sams = sdv_array[j].continuous_variables();
917 for (size_t k=0; k<num_vars; k++)
918 samsOut << std::setw(13) << sams[k] << ' ';
919 samsOut << std::setw(13) << sdr_array[j].response_function();
920 }
921 samsOut << std::endl;
922
923 // Plotting the GP over a grid is intended for visualization and
924 // is therefore only available for 2D problems
925 if (num_vars==2) {
926 std::string gpfile("ego_gp");
927 std::string varfile("ego_var");
928 gpfile += tag;
929 varfile += tag;
930 std::ofstream gpOut(gpfile.c_str(), std::ios::out);
931 std::ofstream varOut(varfile.c_str(), std::ios::out);
932 std::ofstream eifOut("ego_eif.out", std::ios::out);
933 gpOut << std::scientific;
934 varOut << std::scientific;
935 eifOut << std::scientific;
936 RealVector test_pt(2);
937 const RealVector& lbnd = fHatModel.continuous_lower_bounds();
938 const RealVector& ubnd = fHatModel.continuous_upper_bounds();
939 Real interval0 = (ubnd[0] - lbnd[0])/100., interval1 = (ubnd[1] - lbnd[1])/100.;
940 for (size_t j=0; j<101; j++) {
941 test_pt[0] = lbnd[0] + float(j) * interval0;
942 for (size_t k=0; k<101; k++) {
943 test_pt[1] = lbnd[1] + float(k) * interval1;
944
945 fHatModel.continuous_variables(test_pt);
946 fHatModel.evaluate();
947 const Response& gp_resp = fHatModel.current_response();
948 const RealVector& gp_fn = gp_resp.function_values();
949
950 gpOut << '\n' << std::setw(13) << test_pt[0] << ' ' << std::setw(13)
951 << test_pt[1] << ' ' << std::setw(13) << gp_fn[i];
952
953 RealVector variances
954 = fHatModel.approximation_variances(fHatModel.current_variables());
955
956 varOut << '\n' << std::setw(13) << test_pt[0] << ' ' << std::setw(13)
957 << test_pt[1] << ' ' << std::setw(13) << variances[i];
958
959 if (i==numFunctions-1) {
960 Real m = augmented_lagrangian_merit(gp_fn,
961 iteratedModel.primary_response_fn_sense(),
962 iteratedModel.primary_response_fn_weights(),
963 origNonlinIneqLowerBnds, origNonlinIneqUpperBnds,
964 origNonlinEqTargets);
965 RealVector merit(1);
966 merit[0] = m;
967
968 Real ei = expected_improvement(merit, test_pt);
969
970 eifOut << '\n' << std::setw(13) << test_pt[0] << ' '
971 << std::setw(13) << test_pt[1] << ' ' << std::setw(13) << ei;
972 }
973 }
974 gpOut << std::endl;
975 varOut << std::endl;
976 if (i == numFunctions - 1)
977 eifOut << std::endl;
978 }
979 }
980 }
981 #endif //DEBUG_PLOTS
982 }
983
984 } // namespace Dakota
985