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: NonDStochCollocation
10 //- Description: Implementation code for NonDStochCollocation class
11 //- Owner: Mike Eldred
12
13 #include "dakota_system_defs.hpp"
14 #include "NonDStochCollocation.hpp"
15 #include "NonDSparseGrid.hpp"
16 #include "DakotaResponse.hpp"
17 #include "ProblemDescDB.hpp"
18 #include "DataFitSurrModel.hpp"
19 #include "ProbabilityTransformModel.hpp"
20 #include "SharedPecosApproxData.hpp"
21 #include "PecosApproximation.hpp"
22 #include "SharedInterpPolyApproxData.hpp"
23
24 //#define ALLOW_GLOBAL_HERMITE_INTERPOLATION
25 //#define DEBUG
26
27
28 namespace Dakota {
29
30 /** This constructor is called for a standard letter-envelope iterator
31 instantiation using the ProblemDescDB. */
32 NonDStochCollocation::
NonDStochCollocation(ProblemDescDB & problem_db,Model & model)33 NonDStochCollocation(ProblemDescDB& problem_db, Model& model):
34 NonDExpansion(problem_db, model)
35 {
36 // ----------------
37 // Resolve settings
38 // ----------------
39 short data_order,
40 u_space_type = probDescDB.get_short("method.nond.expansion_type");
41 resolve_inputs(u_space_type, data_order);
42
43 // -------------------
44 // Recast g(x) to G(u)
45 // -------------------
46 Model g_u_model;
47 g_u_model.assign_rep(std::make_shared<ProbabilityTransformModel>
48 (iteratedModel, u_space_type)); // retain dist bounds
49
50 // -------------------------
51 // Construct u_space_sampler
52 // -------------------------
53 // LHS/Incremental LHS/Quadrature/SparseGrid samples in u-space
54 // generated using active sampling view:
55 Iterator u_space_sampler;
56 config_integration(probDescDB.get_ushort("method.nond.quadrature_order"),
57 probDescDB.get_ushort("method.nond.sparse_grid_level"),
58 probDescDB.get_rv("method.nond.dimension_preference"),
59 u_space_type, u_space_sampler, g_u_model);
60 String pt_reuse, approx_type;
61 config_approximation_type(approx_type);
62
63 // --------------------------------
64 // Construct G-hat(u) = uSpaceModel
65 // --------------------------------
66 // G-hat(u) uses an orthogonal polynomial approximation over the
67 // active/uncertain variables (using same view as iteratedModel/g_u_model:
68 // not the typical All view for DACE). No correction is employed.
69 // *** Note: for SCBDO with polynomials over {u}+{d}, change view to All.
70 short corr_order = -1, corr_type = NO_CORRECTION;
71 UShortArray approx_order; // empty
72 const ActiveSet& recast_set = g_u_model.current_response().active_set();
73 // DFSModel: consume any QoI aggregation; support surrogate grad evals at most
74 ShortArray asv(g_u_model.qoi(), 3); // for stand alone mode
75 ActiveSet sc_set(asv, recast_set.derivative_vector());
76 String empty_str; // build data import not supported for structured grids
77 uSpaceModel.assign_rep(std::make_shared<DataFitSurrModel>
78 (u_space_sampler, g_u_model,
79 sc_set, approx_type, approx_order, corr_type, corr_order, data_order,
80 outputLevel, pt_reuse, empty_str, TABULAR_ANNOTATED, false,
81 probDescDB.get_string("method.export_approx_points_file"),
82 probDescDB.get_ushort("method.export_approx_format")));
83 initialize_u_space_model();
84
85 // -------------------------------
86 // Construct expSampler, if needed
87 // -------------------------------
88 construct_expansion_sampler(problem_db.get_ushort("method.sample_type"),
89 problem_db.get_string("method.random_number_generator"),
90 problem_db.get_ushort("method.nond.integration_refinement"),
91 problem_db.get_iv("method.nond.refinement_samples"),
92 probDescDB.get_string("method.import_approx_points_file"),
93 probDescDB.get_ushort("method.import_approx_format"),
94 probDescDB.get_bool("method.import_approx_active_only"));
95
96 if (parallelLib.command_line_check())
97 Cout << "\nStochastic collocation construction completed: initial grid "
98 << "size of " << numSamplesOnModel << " evaluations to be performed."
99 << std::endl;
100 }
101
102
103 /** This constructor is used for helper iterator instantiation on the fly. */
104 NonDStochCollocation::
NonDStochCollocation(Model & model,short exp_coeffs_approach,unsigned short num_int,const RealVector & dim_pref,short u_space_type,short refine_type,short refine_control,short covar_control,short rule_nest,short rule_growth,bool piecewise_basis,bool use_derivs)105 NonDStochCollocation(Model& model, short exp_coeffs_approach,
106 unsigned short num_int, const RealVector& dim_pref,
107 short u_space_type, short refine_type,
108 short refine_control, short covar_control,
109 short rule_nest, short rule_growth,
110 bool piecewise_basis, bool use_derivs):
111 NonDExpansion(STOCH_COLLOCATION, model, exp_coeffs_approach, dim_pref, 0,
112 refine_type, refine_control, covar_control, 0., rule_nest,
113 rule_growth, piecewise_basis, use_derivs)
114 // Note: non-zero seed would be needed for expansionSampler, if defined
115 {
116 // ----------------
117 // Resolve settings
118 // ----------------
119 short data_order;
120 resolve_inputs(u_space_type, data_order);
121
122 // -------------------
123 // Recast g(x) to G(u)
124 // -------------------
125 Model g_u_model;
126 g_u_model.assign_rep(std::make_shared<ProbabilityTransformModel>
127 (iteratedModel, u_space_type)); // retain dist bounds
128
129 // -------------------------
130 // Construct u_space_sampler
131 // -------------------------
132 // LHS/Incremental LHS/Quadrature/SparseGrid samples in u-space
133 // generated using active sampling view:
134 Iterator u_space_sampler;
135 config_integration(exp_coeffs_approach, num_int, dim_pref, u_space_sampler,
136 g_u_model);
137 String pt_reuse, approx_type;
138 config_approximation_type(approx_type);
139
140 // --------------------------------
141 // Construct G-hat(u) = uSpaceModel
142 // --------------------------------
143 // G-hat(u) uses an orthogonal polynomial approximation over the
144 // active/uncertain variables (using same view as iteratedModel/g_u_model:
145 // not the typical All view for DACE). No correction is employed.
146 // *** Note: for SCBDO with polynomials over {u}+{d}, change view to All.
147 short corr_order = -1, corr_type = NO_CORRECTION;
148 UShortArray approx_order; // empty
149 const ActiveSet& recast_set = g_u_model.current_response().active_set();
150 // DFSModel: consume any QoI aggregation.
151 // TO DO: support surrogate Hessians in helper mode.
152 ShortArray asv(g_u_model.qoi(), 3); // TO DO: consider passing in data_mode
153 ActiveSet sc_set(asv, recast_set.derivative_vector());
154 uSpaceModel.assign_rep(std::make_shared<DataFitSurrModel>
155 (u_space_sampler, g_u_model,
156 sc_set, approx_type, approx_order, corr_type, corr_order, data_order,
157 outputLevel, pt_reuse));
158 initialize_u_space_model();
159
160 // no expansionSampler, no numSamplesOnExpansion
161 }
162
163
164 /** This constructor is called from derived class constructors that
165 customize the object construction. */
166 NonDStochCollocation::
NonDStochCollocation(unsigned short method_name,ProblemDescDB & problem_db,Model & model)167 NonDStochCollocation(unsigned short method_name, ProblemDescDB& problem_db,
168 Model& model):
169 NonDExpansion(problem_db, model)
170 {
171 // Logic delegated to derived class constructor...
172 }
173
174
175 /** This constructor is called from derived class constructors that
176 customize the object construction. */
177 NonDStochCollocation::
NonDStochCollocation(unsigned short method_name,Model & model,short exp_coeffs_approach,const RealVector & dim_pref,short refine_type,short refine_control,short covar_control,short ml_alloc_control,short ml_discrep,short rule_nest,short rule_growth,bool piecewise_basis,bool use_derivs)178 NonDStochCollocation(unsigned short method_name, Model& model,
179 short exp_coeffs_approach, const RealVector& dim_pref,
180 short refine_type, short refine_control,
181 short covar_control, short ml_alloc_control,
182 short ml_discrep, short rule_nest, short rule_growth,
183 bool piecewise_basis, bool use_derivs):
184 NonDExpansion(method_name, model, exp_coeffs_approach, dim_pref, 0,
185 refine_type, refine_control, covar_control, 0., rule_nest,
186 rule_growth, piecewise_basis, use_derivs)
187 {
188 multilevAllocControl = ml_alloc_control;
189 multilevDiscrepEmulation = ml_discrep;
190
191 // Logic delegated to derived class constructor...
192 }
193
194
~NonDStochCollocation()195 NonDStochCollocation::~NonDStochCollocation()
196 { }
197
198
199 void NonDStochCollocation::
config_integration(unsigned short quad_order,unsigned short ssg_level,const RealVector & dim_pref,short u_space_type,Iterator & u_space_sampler,Model & g_u_model)200 config_integration(unsigned short quad_order, unsigned short ssg_level,
201 const RealVector& dim_pref, short u_space_type,
202 Iterator& u_space_sampler, Model& g_u_model)
203 {
204 // -------------------------
205 // Construct u_space_sampler
206 // -------------------------
207 if (quad_order != USHRT_MAX) {
208 expansionCoeffsApproach = Pecos::QUADRATURE;
209 // TensorProductDriver does not currently support a hierarchical grid via
210 // increment_grid(), etc., although default rule type is set to nested
211 // within NonDExpansion::construct_quadrature() for refinement cases
212 expansionBasisType = Pecos::NODAL_INTERPOLANT;
213 construct_quadrature(u_space_sampler, g_u_model, quad_order, dim_pref);
214 }
215 else if (ssg_level != USHRT_MAX) {
216 switch (expansionBasisType) {
217 case Pecos::HIERARCHICAL_INTERPOLANT:
218 // Note: nestedRules not defined until construct_sparse_grid()
219 if (ruleNestingOverride == Pecos::NON_NESTED) {
220 Cerr << "Error: hierarchical interpolants currently require nested "
221 << "rules. Please remove \"non_nested\" override." << std::endl;
222 abort_handler(-1);
223 }
224 expansionCoeffsApproach = Pecos::HIERARCHICAL_SPARSE_GRID;
225 break;
226 case Pecos::NODAL_INTERPOLANT:
227 expansionCoeffsApproach = (refineControl) ?
228 Pecos::INCREMENTAL_SPARSE_GRID : Pecos::COMBINED_SPARSE_GRID;
229 break;
230 case Pecos::DEFAULT_BASIS:
231 // hierarchical interpolation currently requires nested rules, which
232 // are fairly limited outside of CC, F2, Gauss-Patterson --> use Nodal
233 // as default unless conditions are just right.
234 // Note: nestedRules not defined until construct_sparse_grid()
235 if (u_space_type == STD_UNIFORM_U && refineControl &&
236 ruleNestingOverride != Pecos::NON_NESTED) {//TO DO: rm this constraint
237 expansionCoeffsApproach = Pecos::HIERARCHICAL_SPARSE_GRID;
238 expansionBasisType = Pecos::HIERARCHICAL_INTERPOLANT;
239 }
240 else {
241 expansionCoeffsApproach = (refineControl) ?
242 Pecos::INCREMENTAL_SPARSE_GRID : Pecos::COMBINED_SPARSE_GRID;
243 expansionBasisType = Pecos::NODAL_INTERPOLANT;
244 }
245 break;
246 }
247 /*
248 if (refineControl == Pecos::LOCAL_ADAPTIVE_CONTROL) {
249 if (!piecewiseBasis ||
250 expansionBasisType != Pecos::HIERARCHICAL_INTERPOLANT) {
251 // TO DO: promote this error check to resolve_inputs()
252 PCerr << "Warning: overriding basis type to local hierarchical\n.";
253 piecewiseBasis = true;
254 expansionBasisType = Pecos::HIERARCHICAL_INTERPOLANT;
255 }
256 expansionCoeffsApproach = Pecos::HIERARCHICAL_SPARSE_GRID;
257 }
258 */
259 construct_sparse_grid(u_space_sampler, g_u_model, ssg_level, dim_pref);
260 }
261 }
262
263
264 void NonDStochCollocation::
config_integration(short exp_coeffs_approach,unsigned short num_int,const RealVector & dim_pref,Iterator & u_space_sampler,Model & g_u_model)265 config_integration(short exp_coeffs_approach, unsigned short num_int,
266 const RealVector& dim_pref, Iterator& u_space_sampler,
267 Model& g_u_model)
268 {
269 // -------------------------
270 // Construct u_space_sampler
271 // -------------------------
272 switch (expansionCoeffsApproach) {
273 case Pecos::QUADRATURE:
274 expansionBasisType = Pecos::NODAL_INTERPOLANT;
275 construct_quadrature(u_space_sampler, g_u_model, num_int, dim_pref);
276 break;
277 case Pecos::COMBINED_SPARSE_GRID: case Pecos::INCREMENTAL_SPARSE_GRID:
278 expansionBasisType = Pecos::NODAL_INTERPOLANT;
279 construct_sparse_grid(u_space_sampler, g_u_model, num_int, dim_pref);
280 break;
281 case Pecos::HIERARCHICAL_SPARSE_GRID:
282 expansionBasisType = Pecos::HIERARCHICAL_INTERPOLANT;
283 construct_sparse_grid(u_space_sampler, g_u_model, num_int, dim_pref);
284 break;
285 }
286 }
287
288
config_approximation_type(String & approx_type)289 void NonDStochCollocation::config_approximation_type(String& approx_type)
290 {
291 if (piecewiseBasis)
292 approx_type = (expansionBasisType == Pecos::HIERARCHICAL_INTERPOLANT) ?
293 "piecewise_hierarchical_interpolation_polynomial" :
294 "piecewise_nodal_interpolation_polynomial";
295 else
296 approx_type = (expansionBasisType == Pecos::HIERARCHICAL_INTERPOLANT) ?
297 "global_hierarchical_interpolation_polynomial" :
298 "global_nodal_interpolation_polynomial";
299 }
300
301
resize()302 bool NonDStochCollocation::resize()
303 {
304 bool parent_reinit_comms = NonDExpansion::resize();
305
306 Cerr << "\nError: Resizing is not yet supported in method "
307 << method_enum_to_string(methodName) << "." << std::endl;
308 abort_handler(METHOD_ERROR);
309
310 return parent_reinit_comms;
311 }
312
313
314 void NonDStochCollocation::
resolve_inputs(short & u_space_type,short & data_order)315 resolve_inputs(short& u_space_type, short& data_order)
316 {
317 // perform first due to piecewiseBasis overrides
318 NonDExpansion::resolve_inputs(u_space_type, data_order);
319
320 // There are two derivative cases of interest: (1) derivatives used as
321 // additional data for forming the approximation (derivatives w.r.t. the
322 // expansion variables), and (2) derivatives that will be approximated
323 // separately (derivatives w.r.t. auxilliary variables). The data_order
324 // passed through the DataFitSurrModel defines Approximation::buildDataOrder,
325 // which is restricted to managing the former case. If we need to manage the
326 // latter case in the future, we do not have a good way to detect this state
327 // at construct time, as neither the finalStats ASV/DVV nor subIteratorFlag
328 // have yet been defined. One indicator that is defined is the presence of
329 // inactive continuous vars, since the subModel inactive view is updated
330 // within the NestedModel ctor prior to subIterator instantiation.
331 data_order = 1;
332 if (useDerivs) { // input specification
333 if (iteratedModel.gradient_type() != "none") data_order |= 2;
334 //if (iteratedModel.hessian_type() != "none") data_order |= 4; // not yet
335 #ifdef ALLOW_GLOBAL_HERMITE_INTERPOLATION
336 if (data_order == 1)
337 Cerr << "\nWarning: use_derivatives option in stoch_collocation "
338 << "requires a response\n gradient specification. "
339 << "Option will be ignored.\n" << std::endl;
340 #else
341 if (piecewiseBasis) {
342 if (data_order == 1)
343 Cerr << "\nWarning: use_derivatives option in stoch_collocation "
344 << "requires a response\n gradient specification. "
345 << "Option will be ignored.\n" << std::endl;
346 }
347 else {
348 Cerr << "\nWarning: use of global gradient-enhanced interpolants is "
349 << "disallowed in production\n executables. To activate "
350 << "this research capability, define\n ALLOW_GLOBAL_HERMITE_"
351 << "INTERPOLATION in Dakota::NonDStochCollocation and recompile.\n"
352 << std::endl;
353 data_order = 1;
354 }
355 #endif
356 }
357 useDerivs = (data_order > 1); // override input specification
358
359 // override u_space_type to STD_UNIFORM_U for global Hermite interpolation
360 if (useDerivs && !piecewiseBasis) {
361
362 switch (u_space_type) {
363 //case EXTENDED_U: // default; not user-selectable -> quiet default reassign
364 // break;
365 case ASKEY_U: case PARTIAL_ASKEY_U: // non-default
366 Cerr << "\nWarning: overriding transformation from ASKEY to STD_UNIFORM "
367 << "for Hermite interpolation.\n" << std::endl;
368 break;
369 case STD_NORMAL_U: // non-default
370 Cerr << "\nWarning: overriding transformation from WIENER to STD_UNIFORM "
371 << "for Hermite interpolation.\n" << std::endl;
372 break;
373 }
374
375 u_space_type = STD_UNIFORM_U;
376 }
377 }
378
379
initialize_u_space_model()380 void NonDStochCollocation::initialize_u_space_model()
381 {
382 NonDExpansion::initialize_u_space_model();
383 configure_pecos_options(); // pulled out of base because C3 does not use it
384
385 // initialize product accumulators with PolynomialApproximation pointers
386 // used in covariance calculations
387 if ( expansionBasisType == Pecos::HIERARCHICAL_INTERPOLANT &&
388 refineControl && ( refineMetric == Pecos::COVARIANCE_METRIC ||
389 refineMetric == Pecos::MIXED_STATS_METRIC ) )
390 initialize_covariance();
391
392 // Precedes construct_basis() since basis is stored in Pecos driver
393 SharedApproxData& shared_data = uSpaceModel.shared_approximation();
394 shared_data.integration_iterator(uSpaceModel.subordinate_iterator());
395
396 // DataFitSurrModel copies u-space mvDist from ProbabilityTransformModel
397 const Pecos::MultivariateDistribution& u_mvd
398 = uSpaceModel.multivariate_distribution();
399 // construct the polynomial basis (shared by integration drivers)
400 shared_data.construct_basis(u_mvd);
401 // mainly a run-time requirement, but also needed at construct time
402 // (e.g., to initialize NumericGenOrthogPolynomial::distributionType)
403 //shared_data.update_basis_distribution_parameters(u_mvd);
404
405 initialize_u_space_grid();
406 }
407
408
initialize_covariance()409 void NonDStochCollocation::initialize_covariance()
410 {
411 std::vector<Approximation>& poly_approxs = uSpaceModel.approximations();
412 size_t i, j;
413 for (i=0; i<numFunctions; ++i) {
414 std::shared_ptr<PecosApproximation> pa_rep_i =
415 std::static_pointer_cast<PecosApproximation>(
416 poly_approxs[i].approx_rep());
417 pa_rep_i->clear_covariance_pointers();
418 for (j=0; j<=i; ++j)
419 pa_rep_i->initialize_covariance(poly_approxs[j]);
420 }
421 }
422
423
compute_delta_mean(bool update_ref)424 void NonDStochCollocation::compute_delta_mean(bool update_ref)
425 {
426 std::vector<Approximation>& poly_approxs = uSpaceModel.approximations();
427 bool warn_flag = false,
428 combined_stats = (statsMetricMode == Pecos::COMBINED_EXPANSION_STATS);
429
430 if (deltaRespMean.empty()) deltaRespMean.sizeUninitialized(numFunctions);
431 for (size_t i=0; i<numFunctions; ++i) {
432 std::shared_ptr<PecosApproximation> pa_rep_i =
433 std::static_pointer_cast<PecosApproximation>(
434 poly_approxs[i].approx_rep());
435 if (pa_rep_i->expansion_coefficient_flag()) {
436 if (combined_stats)
437 // refinement assessed for impact on combined expansion from roll up
438 deltaRespMean[i] = (allVars) ?
439 pa_rep_i->delta_combined_mean(initialPtU) :
440 pa_rep_i->delta_combined_mean();
441 else // refinement assessed for impact on the current expansion
442 deltaRespMean[i] = (allVars) ?
443 pa_rep_i->delta_mean(initialPtU) : pa_rep_i->delta_mean();
444
445 if (update_ref) {
446 if (combined_stats) {
447 Real new_mean = pa_rep_i->combined_moment(0) + deltaRespMean[i];
448 pa_rep_i->combined_moment(new_mean, 0);
449 }
450 else {
451 Real new_mean = pa_rep_i->moment(0) + deltaRespMean[i];
452 pa_rep_i->moment(new_mean, 0);
453 }
454 }
455 }
456 else
457 { warn_flag = true; deltaRespMean[i] = 0.; }
458 }
459
460 // no current need for printing mean values by themselves:
461 //if (print_metric) print_mean(Cout, deltaRespMean, "Change in");
462 // mean values not tracked outside PolynomialApproximation:
463 //if (update_ref) respMean += deltaRespMean;
464 if (warn_flag)
465 Cerr << "Warning: expansion coefficients unavailable in NonD"
466 << "StochCollocation::compute_delta_mean().\n "
467 << "Zeroing affected deltaRespMean terms." << std::endl;
468 }
469
470
471 void NonDStochCollocation::
compute_delta_variance(bool update_ref,bool print_metric)472 compute_delta_variance(bool update_ref, bool print_metric)
473 {
474 std::vector<Approximation>& poly_approxs = uSpaceModel.approximations();
475 bool warn_flag = false,
476 combined_stats = (statsMetricMode == Pecos::COMBINED_EXPANSION_STATS);
477
478 if (deltaRespVariance.empty())
479 deltaRespVariance.sizeUninitialized(numFunctions);
480 for (size_t i=0; i<numFunctions; ++i) {
481 std::shared_ptr<PecosApproximation> pa_rep_i =
482 std::static_pointer_cast<PecosApproximation>
483 (poly_approxs[i].approx_rep());
484 Real& delta = deltaRespVariance[i];
485 if (pa_rep_i->expansion_coefficient_flag()) {
486 if (combined_stats)
487 // refinement assessed for impact on combined expansion from roll up
488 delta = (allVars) ? pa_rep_i->delta_combined_variance(initialPtU) :
489 pa_rep_i->delta_combined_variance();
490 else // refinement assessed for impact on the current expansion
491 delta = (allVars) ? pa_rep_i->delta_variance(initialPtU) :
492 pa_rep_i->delta_variance();
493
494 if (update_ref) {
495 respVariance[i] += delta;
496 if (combined_stats) pa_rep_i->combined_moment(respVariance[i], 1);
497 else pa_rep_i->moment(respVariance[i], 1);
498 }
499 }
500 else
501 { warn_flag = true; delta = 0.; }
502 }
503
504 if (print_metric) print_variance(Cout, deltaRespVariance, "Change in");
505 if (warn_flag)
506 Cerr << "Warning: expansion coefficients unavailable in NonDStoch"
507 << "Collocation::compute_delta_variance().\n "
508 << "Zeroing affected deltaRespVariance terms." << std::endl;
509 }
510
511
512 void NonDStochCollocation::
compute_delta_covariance(bool update_ref,bool print_metric)513 compute_delta_covariance(bool update_ref, bool print_metric)
514 {
515 std::vector<Approximation>& poly_approxs = uSpaceModel.approximations();
516 bool warn_flag = false,
517 combined_stats = (statsMetricMode == Pecos::COMBINED_EXPANSION_STATS);
518 size_t i, j;
519
520 if (deltaRespCovariance.empty())
521 deltaRespCovariance.shapeUninitialized(numFunctions);
522 for (i=0; i<numFunctions; ++i) {
523 std::shared_ptr<PecosApproximation> pa_rep_i =
524 std::static_pointer_cast<PecosApproximation>
525 (poly_approxs[i].approx_rep());
526 if (pa_rep_i->expansion_coefficient_flag()) {
527 for (j=0; j<=i; ++j) {
528 Approximation& approx_j = poly_approxs[j];
529 Real& delta = deltaRespCovariance(i,j);
530 if (approx_j.expansion_coefficient_flag()) {
531 if (combined_stats)
532 // refinement assessed for impact on combined exp from roll up
533 delta = (allVars) ?
534 pa_rep_i->delta_combined_covariance(initialPtU, approx_j) :
535 pa_rep_i->delta_combined_covariance(approx_j);
536 else // refinement assessed for impact on the current expansion
537 delta = (allVars) ?
538 pa_rep_i->delta_covariance(initialPtU, approx_j) :
539 pa_rep_i->delta_covariance(approx_j);
540
541 if (update_ref) {
542 respCovariance(i,j) += delta;
543 if (i == j) {
544 if (combined_stats)
545 pa_rep_i->combined_moment(respCovariance(i,i), 1);
546 else
547 pa_rep_i->moment(respCovariance(i,i), 1);
548 }
549 }
550 }
551 else
552 { warn_flag = true; delta = 0.; }
553 }
554 }
555 else {
556 warn_flag = true;
557 for (j=0; j<=i; ++j)
558 deltaRespCovariance(i,j) = 0.;
559 }
560 }
561
562 if (print_metric) print_covariance(Cout, deltaRespCovariance, "Change in");
563 if (warn_flag)
564 Cerr << "Warning: expansion coefficients unavailable in NonDStoch"
565 << "Collocation::compute_delta_covariance().\n "
566 << "Zeroing affected deltaRespCovariance terms." << std::endl;
567 }
568
569
570 Real NonDStochCollocation::
compute_covariance_metric(bool revert,bool print_metric)571 compute_covariance_metric(bool revert, bool print_metric)
572 {
573 if (expansionBasisType == Pecos::HIERARCHICAL_INTERPOLANT) {
574 bool update_ref = !revert;
575
576 // augment delta covar -> ref covar with delta_mean -> ref mean;
577 // > supports base compute/print requirements without incurring overhead
578 // of metric_roll_up(), avoiding need to more broadly redefine
579 // compute_statistics(), print_results(), pull_*(), push_*() to
580 // exclude compute_moments()
581 compute_delta_mean(update_ref);
582
583 // Metric scale is determined from reference covariance. While defining
584 // the scale from an updated covariance would eliminate problems with zero
585 // covariance for adaptations from level 0, different refinement candidates
586 // would score equally at 1 (induced 100% of change in updated covariance)
587 // in this initial set of candidates. Therefore, use reference covariance
588 // as the scale and mitigate underflow of its norm.
589 Real scale, delta_norm;
590 switch (covarianceControl) {
591 case DIAGONAL_COVARIANCE:
592 if (relativeMetric) // norm of reference variance, bounded from zero
593 scale = std::max(Pecos::SMALL_NUMBER, respVariance.normFrobenius());
594 compute_delta_variance(update_ref, print_metric);
595 delta_norm = deltaRespVariance.normFrobenius();
596 break;
597 case FULL_COVARIANCE:
598 if (relativeMetric) // norm of reference covariance, bounded from zero
599 scale = std::max(Pecos::SMALL_NUMBER, respCovariance.normFrobenius());
600 compute_delta_covariance(update_ref, print_metric);
601 delta_norm = deltaRespCovariance.normFrobenius();
602 break;
603 }
604
605 return (relativeMetric) ? delta_norm / scale : delta_norm;
606 }
607 else // use default implementation
608 return NonDExpansion::compute_covariance_metric(revert, print_metric);
609 }
610
611
612 Real NonDStochCollocation::
compute_level_mappings_metric(bool revert,bool print_metric)613 compute_level_mappings_metric(bool revert, bool print_metric)
614 {
615 // Focus only on level mappings and neglect moment deltas
616
617 // combine delta_beta() and delta_z() from HierarchInterpPolyApproximation
618 // with default definition of delta-{p,beta*}
619
620 if (expansionBasisType == Pecos::HIERARCHICAL_INTERPOLANT) {
621
622 // ensure moment updates for mixed stats:
623 if (refineMetric == Pecos::MIXED_STATS_METRIC) {
624 bool update_ref = !revert;
625 compute_delta_mean(update_ref);
626 switch (covarianceControl) {
627 case DIAGONAL_COVARIANCE:
628 compute_delta_variance(update_ref, false); break;
629 case FULL_COVARIANCE:
630 compute_delta_covariance(update_ref, false); break;
631 }
632 }
633
634 // Note: it would be desirable to include support for all active statistics,
635 // including delta_mean() and delta_std_deviation(). With access to nested
636 // response mappings passed down from an outer context, a more comprehensive
637 // set of stats could be supported in the logic below.
638
639 bool beta_map = false, numerical_map = false; size_t i, j, cntr;
640 for (i=0; i<numFunctions; ++i) {
641 if (!requestedRelLevels[i].empty()) beta_map = true;
642 if (!requestedProbLevels[i].empty() || !requestedGenRelLevels[i].empty())
643 numerical_map = true;
644 if (!requestedRespLevels[i].empty()) {
645 if (respLevelTarget == RELIABILITIES) beta_map = true;
646 else numerical_map = true;
647 }
648 }
649 if (beta_map) { // hierarchical increments in beta-bar->z and z-bar->beta
650
651 size_t offset = 0;
652 RealVector level_maps_ref, level_maps_new;
653 pull_level_mappings(level_maps_ref, offset);
654 if (numerical_map) { // merge in z-bar->p,beta* & p-bar,beta*-bar->z
655 //metric_roll_up(REFINEMENT_RESULTS); // TO DO: support combined exp in numerical stats
656 compute_numerical_level_mappings();
657 pull_level_mappings(level_maps_new, offset);// analytic overlaid at end
658 deltaLevelMaps = level_maps_new; deltaLevelMaps -= level_maps_ref;
659 }
660 else {
661 deltaLevelMaps.size(totalLevelRequests); // init to 0
662 if (!revert) level_maps_new.size(totalLevelRequests); // init to 0
663 }
664
665 bool warn_flag = false,
666 combined_stats = (statsMetricMode == Pecos::COMBINED_EXPANSION_STATS);
667 std::vector<Approximation>& poly_approxs = uSpaceModel.approximations();
668 Real delta, ref, sum_sq = 0., scale_sq = 0., z_bar, beta_bar;
669 for (i=0, cntr=0; i<numFunctions; ++i) {
670 size_t rl_len = requestedRespLevels[i].length(),
671 pl_len = requestedProbLevels[i].length(),
672 bl_len = requestedRelLevels[i].length(),
673 gl_len = requestedGenRelLevels[i].length();
674 std::shared_ptr<PecosApproximation> pa_rep_i =
675 std::static_pointer_cast<PecosApproximation>
676 (poly_approxs[i].approx_rep());
677 if (pa_rep_i->expansion_coefficient_flag()) {
678 if (respLevelTarget == RELIABILITIES)
679 for (j=0; j<rl_len; ++j, ++cntr) {
680 z_bar = requestedRespLevels[i][j];
681 if (combined_stats)
682 delta = deltaLevelMaps[cntr] = (allVars) ?
683 pa_rep_i->delta_combined_beta(initialPtU, cdfFlag, z_bar) :
684 pa_rep_i->delta_combined_beta(cdfFlag, z_bar);
685 else
686 delta = deltaLevelMaps[cntr] = (allVars) ?
687 pa_rep_i->delta_beta(initialPtU, cdfFlag, z_bar) :
688 pa_rep_i->delta_beta(cdfFlag, z_bar);
689 sum_sq += delta * delta;
690 // avoid proliferating numerical exception:
691 ref = level_maps_ref[cntr];
692 if (std::abs(ref) == Pecos::LARGE_NUMBER) {
693 // ref is undefined and delta neglects term;
694 // recompute new w/o delta in analytic_delta_level_mappings()
695
696 // do not increment scale for dummy beta value --> may result
697 // in SMALL_NUMBER scaling below if no meaningful refs exist
698 }
699 else if (relativeMetric)
700 scale_sq += ref * ref;// ref,delta are valid --> update scale
701 }
702 else
703 for (j=0; j<rl_len; ++j, ++cntr) {
704 delta = deltaLevelMaps[cntr]; sum_sq += delta * delta;
705 if (relativeMetric)
706 { ref = level_maps_ref[cntr]; scale_sq += ref * ref; }
707 }
708 for (j=0; j<pl_len; ++j, ++cntr) {
709 delta = deltaLevelMaps[cntr]; sum_sq += delta * delta;
710 if (relativeMetric)
711 { ref = level_maps_ref[cntr]; scale_sq += ref * ref; }
712 }
713 for (j=0; j<bl_len; ++j, ++cntr) {
714 beta_bar = requestedRelLevels[i][j];
715 if (combined_stats)
716 delta = deltaLevelMaps[cntr] = (allVars) ?
717 pa_rep_i->delta_combined_z(initialPtU, cdfFlag, beta_bar) :
718 pa_rep_i->delta_combined_z(cdfFlag, beta_bar);
719 else
720 delta = deltaLevelMaps[cntr] = (allVars) ?
721 pa_rep_i->delta_z(initialPtU, cdfFlag, beta_bar) :
722 pa_rep_i->delta_z(cdfFlag, beta_bar);
723 sum_sq += delta * delta;
724 ref = level_maps_ref[cntr];
725 if (relativeMetric) scale_sq += ref * ref;
726 }
727 for (j=0; j<gl_len; ++j, ++cntr) {
728 delta = deltaLevelMaps[cntr]; sum_sq += delta * delta;
729 if (relativeMetric)
730 { ref = level_maps_ref[cntr]; scale_sq += ref * ref; }
731 }
732 }
733 else {
734 warn_flag = true;
735 cntr += rl_len + pl_len + bl_len + gl_len;
736 }
737 }
738 if (warn_flag)
739 Cerr << "Warning: expansion coefficients unavailable in "
740 << "NonDStochCollocation::compute_level_mappings_metric().\n"
741 << " Omitting affected level mappings." << std::endl;
742 // As for compute_delta_covariance(), print level mapping deltas
743 // (without a moment offset as for final_stats):
744 if (print_metric)
745 print_level_mappings(Cout, deltaLevelMaps, false, "Change in");
746
747 // Level mappings: promote to new or revert to previous (if required)
748 if (!revert) { // retain updated values
749 analytic_delta_level_mappings(level_maps_ref, level_maps_new);
750 push_level_mappings(level_maps_new, offset);
751 }
752 else if (numerical_map) // restore ref values that were overwritten
753 push_level_mappings(level_maps_ref, offset); // restore reference
754 //else deltaLevelMaps does not impact existing level mappings
755
756 // Metric scale is determined from reference stats, not updated stats,
757 // as consistent with compute_covariance_metric() above.
758 if (relativeMetric) {
759 Real scale = std::max(Pecos::SMALL_NUMBER, std::sqrt(scale_sq));
760 return std::sqrt(sum_sq) / scale;
761 }
762 else
763 return std::sqrt(sum_sq);
764 }
765 else // use default implementation if no beta-mapping increments
766 return
767 NonDExpansion::compute_level_mappings_metric(revert, print_metric);
768 }
769 else // use default implementation for Nodal
770 return NonDExpansion::compute_level_mappings_metric(revert, print_metric);
771 }
772
773
774 /*
775 Real NonDStochCollocation::
776 compute_final_statistics_metric(bool revert, bool print_metric)
777 {
778 // Focus only on level mappings and neglect moment deltas
779
780 // combine delta_beta() and delta_z() from HierarchInterpPolyApproximation
781 // with default definition of delta-{p,beta*}
782
783 if (expansionBasisType == Pecos::HIERARCHICAL_INTERPOLANT) {
784 // Note: it would be desirable to include support for all active statistics,
785 // including delta_mean() and delta_std_deviation(). With access to nested
786 // response mappings passed down from an outer context, a more comprehensive
787 // set of stats could be supported in the logic below.
788 bool beta_map = false, numerical_map = false,
789 combined_stats = (statsMetricMode == Pecos::COMBINED_EXPANSION_STATS);
790 size_t i, j, cntr;
791 for (i=0; i<numFunctions; ++i) {
792 if (!requestedRelLevels[i].empty()) beta_map = true;
793 if (!requestedProbLevels[i].empty() || !requestedGenRelLevels[i].empty())
794 numerical_map = true;
795 if (!requestedRespLevels[i].empty()) {
796 if (respLevelTarget == RELIABILITIES) beta_map = true;
797 else numerical_map = true;
798 }
799 }
800 if (beta_map) { // hierarchical increments in beta-bar->z and z-bar->beta
801
802 RealVector delta_final_stats, final_stats_new,
803 final_stats_ref = finalStatistics.function_values();
804 // *** Note: this requires that the reference includes FINAL_RESULTS,
805 // *** which is not currently true (only INTERMEDIATE_RESULTS)
806 if (numerical_map) { // merge in z-bar->p,beta* & p-bar,beta*-bar->z
807 compute_statistics(FINAL_RESULTS);//no finalStats for REFINEMENT_RESULTS
808 delta_final_stats = final_stats_new = finalStatistics.function_values();
809 delta_final_stats -= final_stats_ref; // compute delta
810 }
811 else {
812 size_t num_stats = finalStatistics.num_functions();
813 delta_final_stats.size(num_stats); // init to 0
814 if (!revert) final_stats_new.size(num_stats); // init to 0
815 }
816
817 bool warn_flag = false;
818 std::vector<Approximation>& poly_approxs = uSpaceModel.approximations();
819 Real delta, ref, sum_sq = 0., scale_sq = 0., z_bar, beta_bar;
820 for (i=0, cntr=0; i<numFunctions; ++i) {
821 size_t rl_len = requestedRespLevels[i].length(),
822 pl_len = requestedProbLevels[i].length(),
823 bl_len = requestedRelLevels[i].length(),
824 gl_len = requestedGenRelLevels[i].length();
825 std::shared_ptr<PecosApproximation> pa_rep_i =
826 std::static_pointer_cast<PecosApproximation>
827 (poly_approxs[i].approx_rep());
828 cntr += moment_offset; // skip moments if present
829 if (pa_rep_i->expansion_coefficient_flag()) {
830 if (respLevelTarget == RELIABILITIES)
831 for (j=0; j<rl_len; ++j, ++cntr) {
832 z_bar = requestedRespLevels[i][j];
833 if (combined_stats)
834 delta = delta_final_stats[cntr] = (allVars) ?
835 pa_rep_i->delta_combined_beta(initialPtU, cdfFlag, z_bar) :
836 pa_rep_i->delta_combined_beta(cdfFlag, z_bar);
837 else
838 delta = delta_final_stats[cntr] = (allVars) ?
839 pa_rep_i->delta_beta(initialPtU, cdfFlag, z_bar) :
840 pa_rep_i->delta_beta(cdfFlag, z_bar);
841 sum_sq += delta * delta;
842 ref = final_stats_ref[cntr];
843 // Note: this captures the more likely of the Pecos::
844 // HierarchInterpPolyApproximation::delta_beta_map() exceptions
845 // (sigma_ref = 0), but not rare case of sigma_new = 0 by itself.
846 if (std::abs(ref) == Pecos::LARGE_NUMBER) {
847 // ref is undefined and delta neglects term; must compute new
848 if (!revert) {
849 if (combined_stats)
850 final_stats_new[cntr] = (allVars) ?
851 pa_rep_i->combined_beta(initialPtU, cdfFlag, z_bar) :
852 pa_rep_i->combined_beta(cdfFlag, z_bar);
853 else
854 final_stats_new[cntr] = (allVars) ?
855 pa_rep_i->beta(initialPtU, cdfFlag, z_bar) :
856 pa_rep_i->beta(cdfFlag, z_bar);
857 }
858 // do not increment scale for dummy beta value --> may result
859 // in SMALL_NUMBER scaling below if no meaningful refs exist
860 }
861 else { // ref and delta are valid --> update scale and new
862 if (relativeMetric) scale_sq += ref * ref;
863 if (!revert) final_stats_new[cntr] = ref + delta;
864 }
865 }
866 else
867 for (j=0; j<rl_len; ++j, ++cntr) {
868 delta = delta_final_stats[cntr]; sum_sq += delta * delta;
869 if (relativeMetric)
870 { ref = final_stats_ref[cntr]; scale_sq += ref * ref; }
871 }
872 for (j=0; j<pl_len; ++j, ++cntr) {
873 delta = delta_final_stats[cntr]; sum_sq += delta * delta;
874 if (relativeMetric)
875 { ref = final_stats_ref[cntr]; scale_sq += ref * ref; }
876 }
877 for (j=0; j<bl_len; ++j, ++cntr) {
878 beta_bar = requestedRelLevels[i][j];
879 if (combined_stats)
880 delta = delta_final_stats[cntr] = (allVars) ?
881 pa_rep_i->delta_combined_z(initialPtU, cdfFlag, beta_bar) :
882 pa_rep_i->delta_combined_z(cdfFlag, beta_bar);
883 else
884 delta = delta_final_stats[cntr] = (allVars) ?
885 pa_rep_i->delta_z(initialPtU, cdfFlag, beta_bar) :
886 pa_rep_i->delta_z(cdfFlag, beta_bar);
887 sum_sq += delta * delta;
888 ref = final_stats_ref[cntr];
889 if (relativeMetric) scale_sq += ref * ref;
890 if (!revert) final_stats_new[cntr] = ref + delta;
891 }
892 for (j=0; j<gl_len; ++j, ++cntr) {
893 delta = delta_final_stats[cntr]; sum_sq += delta * delta;
894 if (relativeMetric)
895 { ref = final_stats_ref[cntr]; scale_sq += ref * ref; }
896 }
897 }
898 else {
899 warn_flag = true;
900 cntr += rl_len + pl_len + bl_len + gl_len;
901 }
902 }
903 if (warn_flag)
904 Cerr << "Warning: expansion coefficients unavailable in "
905 << "NonDStochCollocation::compute_final_statistics_metric().\n"
906 << " Omitting affected final statistics." << std::endl;
907
908 // As for compute_delta_covariance(), print level mapping deltas:
909 if (print_metric) {
910 bool moments = (finalMomentsType > NO_MOMENTS);
911 print_level_mappings(Cout, delta_final_stats, moments, "Change in");
912 }
913 // Final stats: revert to previous or promote to new
914 if (!revert)
915 finalStatistics.function_values(final_stats_new);
916 else if (numerical_map) // revert from final_stats_new to final_stats_ref
917 finalStatistics.function_values(final_stats_ref);
918
919 // Metric scale is determined from reference stats, not updated stats,
920 // as consistent with compute_covariance_metric() above.
921 if (relativeMetric) {
922 Real scale = std::max(Pecos::SMALL_NUMBER, std::sqrt(scale_sq));
923 return std::sqrt(sum_sq) / scale;
924 }
925 else
926 return std::sqrt(sum_sq);
927 }
928 else // use default implementation if no beta-mapping increments
929 return
930 NonDExpansion::compute_final_statistics_metric(revert, print_metric);
931 }
932 else // use default implementation for Nodal
933 return NonDExpansion::compute_final_statistics_metric(revert, print_metric);
934 }
935 */
936
937
938 /* Calculate analytic and numerical statistics from the expansion and
939 log results within final_stats for use in OUU.
940 void NonDStochCollocation::compute_statistics(short results_state)
941 {
942 if (expansionBasisType == Pecos::HIERARCHICAL_INTERPOLANT)
943 // specialize base implementation since metric_roll_up() is not used and
944 // extra stats that require roll-up no longer come for free
945 switch (results_state) {
946 case REFINEMENT_RESULTS:
947 // compute_{covariance,level_mapping,final_statistics}_metric() performs
948 // only the necessary computations for resolving delta.norm()
949
950 // augment delta covar -> ref covar with delta_mean -> ref mean;
951 // > supports base compute/print requirements without incurring overhead
952 // of metric_roll_up(), avoiding need to more broadly redefine
953 // compute_statistics(), print_results(), pull_*(), push_*() to
954 // exclude compute_moments()
955 compute_delta_mean(true); // update ref
956 break;
957 //case INTERMEDIATE_RESULTS: case FINAL_RESULTS: break;
958 }
959
960 NonDExpansion::compute_statistics(results_state);
961 }
962
963
964 void NonDStochCollocation::pull_candidate(RealVector& stats_star)
965 {
966 if (expansionBasisType == Pecos::HIERARCHICAL_INTERPOLANT)
967 // If pulling updated values as in NonDExpansion
968 switch (refineMetric) {
969 case Pecos::COVARIANCE_METRIC: {
970 std::vector<Approximation>& poly_approxs = uSpaceModel.approximations();
971 std::shared_ptr<PecosApproximation> poly_approx_rep;
972 bool full_covar = (covarianceControl == FULL_COVARIANCE),
973 combined_stats = (statsMetricMode == Pecos::COMBINED_EXPANSION_STATS);
974 size_t vec_len = (full_covar) ?
975 (numFunctions*(numFunctions + 3))/2 : 2*numFunctions;
976 if (stats_star.length() != vec_len) stats_star.sizeUninitialized(vec_len);
977 // pull means
978 for (size_t i=0; i<numFunctions; ++i) {
979 poly_approx_rep =
980 std::static_pointer_cast<PecosApproximation>(
981 poly_approxs[i].approx_rep());
982 stats_star[i] = (combined_stats) ?
983 poly_approx_rep->combined_moment(0) + deltaRespMean[i] :
984 poly_approx_rep->moment(0) + deltaRespMean[i];
985 }
986 // pull resp{V,Cov}ariance
987 if (full_covar) {
988 RealSymMatrix new_covar(respCovariance);
989 new_covar += deltaRespCovariance;
990 pull_lower_triangle(new_covar, stats_star, numFunctions);
991 }
992 else
993 for (size_t i=0; i<numFunctions; ++i)
994 stats_star[i+numFunctions] = respVariance[i] + deltaRespVariance[i];
995 break;
996 }
997 default:
998 // define an offset for MIXED_STATS_METRIC
999 pull_level_mappings(stats_star, offset); // pull updated numerical stats
1000 analytic_delta_level_mappings(stats_star, stats_star); // update analytic
1001 // (stats_star provides ref and becomes new for these entries)
1002 break;
1003 }
1004 else
1005 NonDExpansion::pull_candidate(stats_star);
1006 }
1007 */
1008
1009
1010 /** In this function, we leave numerical stats alone, updating
1011 analytic level stats either using ref+delta or, if ref is invalid,
1012 though recomputation. */
1013 void NonDStochCollocation::
analytic_delta_level_mappings(const RealVector & level_maps_ref,RealVector & level_maps_new)1014 analytic_delta_level_mappings(const RealVector& level_maps_ref,
1015 RealVector& level_maps_new)
1016 {
1017 // preserve existing as this may be an overlay
1018 if (level_maps_new.length() != totalLevelRequests)
1019 level_maps_new.resize(totalLevelRequests);
1020
1021 size_t i, j, cntr, rl_len, pl_len, bl_len, gl_len, pl_bl_gl_len;
1022 std::vector<Approximation>& poly_approxs = uSpaceModel.approximations();
1023 Real delta, ref, sum_sq = 0., scale_sq = 0., z_bar, beta_bar;
1024 bool combined_stats = (statsMetricMode == Pecos::COMBINED_EXPANSION_STATS);
1025 for (i=0, cntr=0; i<numFunctions; ++i) {
1026 rl_len = requestedRespLevels[i].length();
1027 pl_len = requestedProbLevels[i].length();
1028 bl_len = requestedRelLevels[i].length();
1029 gl_len = requestedGenRelLevels[i].length();
1030 pl_bl_gl_len = pl_len+bl_len+gl_len;
1031 std::shared_ptr<PecosApproximation> pa_rep_i =
1032 std::static_pointer_cast<PecosApproximation>(
1033 poly_approxs[i].approx_rep());
1034 if (respLevelTarget == RELIABILITIES)
1035 for (j=0; j<rl_len; ++j, ++cntr) {
1036 // Note: this captures the more likely of the Pecos::
1037 // HierarchInterpPolyApproximation::delta_beta_map() exceptions
1038 // (sigma_ref = 0), but not rare case of sigma_new = 0 by itself.
1039 ref = level_maps_ref[cntr];
1040 if (std::abs(ref) == Pecos::LARGE_NUMBER) {
1041 // ref is undefined and delta neglects term; must compute new
1042 z_bar = requestedRespLevels[i][j];
1043 if (combined_stats)
1044 level_maps_new[cntr] = (allVars) ?
1045 pa_rep_i->combined_beta(initialPtU, cdfFlag, z_bar) :
1046 pa_rep_i->combined_beta(cdfFlag, z_bar);
1047 else
1048 level_maps_new[cntr] = (allVars) ?
1049 pa_rep_i->beta(initialPtU, cdfFlag, z_bar) :
1050 pa_rep_i->beta(cdfFlag, z_bar);
1051 }
1052 else // ref and delta are valid
1053 level_maps_new[cntr] = ref + deltaLevelMaps[cntr];
1054 }
1055 else
1056 cntr += rl_len;
1057 cntr += pl_len;
1058 for (j=0; j<bl_len; ++j)
1059 level_maps_new[cntr] = level_maps_ref[cntr] + deltaLevelMaps[cntr];
1060 cntr += gl_len;
1061 }
1062 }
1063
1064 } // namespace Dakota
1065