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: ApproximationInterface
10 //- Description: Implementation of base class for approximation interfaces
11 //- Owner: Mike Eldred
12
13 #include "dakota_system_defs.hpp"
14 #include "dakota_tabular_io.hpp"
15 #include "ApproximationInterface.hpp"
16 #include "DakotaVariables.hpp"
17 #include "DakotaResponse.hpp"
18 #include "ProblemDescDB.hpp"
19 #include "PRPMultiIndex.hpp"
20
21 //#define DEBUG
22
23
24 namespace Dakota {
25
26 extern PRPCache data_pairs;
27
28 size_t ApproximationInterface::approxIdNum = 0;
29
30 ApproximationInterface::
ApproximationInterface(ProblemDescDB & problem_db,const Variables & am_vars,bool am_cache,const String & am_interface_id,const StringArray & fn_labels)31 ApproximationInterface(ProblemDescDB& problem_db, const Variables& am_vars,
32 bool am_cache, const String& am_interface_id,
33 const StringArray& fn_labels):
34 Interface(BaseConstructor(), problem_db),
35 approxFnIndices(problem_db.get_is("model.surrogate.function_indices")),
36 //graph3DFlag(problem_db.get_bool("environment.graphics")),
37 challengeFile(problem_db.get_string("model.surrogate.challenge_points_file")),
38 challengeFormat(
39 problem_db.get_ushort("model.surrogate.challenge_points_file_format")),
40 challengeUseVarLabels(
41 problem_db.get_bool("model.surrogate.challenge_use_variable_labels")),
42 challengeActiveOnly(
43 problem_db.get_bool("model.surrogate.challenge_points_file_active")),
44 actualModelVars(am_vars.copy()), actualModelCache(am_cache),
45 actualModelInterfaceId(am_interface_id)
46 {
47 // Some specification-based attributes inherited from Interface may
48 // be incorrect since there is no longer an approximation interface
49 // specification (assign_rep() is used from DataFitSurrModel).
50 // Override these inherited settings.
51 interfaceId = String("APPROX_INTERFACE_") + std::to_string(++approxIdNum);
52 interfaceType = APPROX_INTERFACE;
53 algebraicMappings = false; // for now; *** TO DO ***
54
55 // process approxFnIndices. IntSets are sorted and unique. Error checking
56 // is performed in SurrogateModel and does not need to be replicated here.
57 size_t num_fns = fn_labels.size();
58 if (approxFnIndices.empty()) // default: all fns are approximated
59 for (int i=0; i<num_fns; i++)
60 approxFnIndices.insert(i);
61
62 // This section moved to the constructor so that ApproxInterfaces can be
63 // queried for their state prior to build_approximation() (e.g., to configure
64 // parallelism w/ max_concurrency). Also, Approximation classes now
65 // use the problem_db, so their instantiation must occur in constructors to
66 // assure proper list node settings.
67 functionSurfaces.resize(num_fns);
68 // despite view mappings, x in map() always = size of active actualModelVars
69 size_t num_vars = actualModelVars.cv() + actualModelVars.div()
70 + actualModelVars.dsv() + actualModelVars.drv();
71 sharedData = SharedApproxData(problem_db, num_vars);
72 for (ISIter it=approxFnIndices.begin(); it!=approxFnIndices.end(); ++it)
73 functionSurfaces[*it] = Approximation(problem_db, sharedData,
74 fn_labels[*it]);
75
76 /*
77 // Old approach for scaling/offsetting approximation results read the data
78 // from a hardwired filename (approx_scale_offset.in). This capability had
79 // been used (~Spring 2000) to reuse existing surrogates with modified design
80 // goals, particularly modified constraint allowables via approximation
81 // offsets. The intent was to eventually make this a supported capability
82 // within the interface keyword of the input specification, e.g.:
83 // [approx_scale = <LISTof><REAL>]
84 // [approx_offset = <LISTof><REAL>]
85 // An outstanding question was which approximation types to support, i.e.
86 // (1) global: primary use, (2) local/multipoint: somewhat useful, and
87 // (3) hierarchical: not supported in HierarchSurrModel and not that useful
88 // since correction is required.
89 //
90 // Since that time, the capability for specifying inequality constraint bounds
91 // and equality constraint targets (~Spring 2001) has largely superceded the
92 // capability to scale/offset approximation results, since modifying the
93 // constraint bounds/targets for subsequent surrogate-based runs has the same
94 // effect. However, if the need to offset/scale surrogate objective functions
95 // or least squares terms arises, then this is not supported in the existing
96 // code and the block below (as well as corresponding lines in map()) could be
97 // reactivated.
98 ifstream scale_offset("approx_scale_offset.in");
99 if (scale_offset) {
100 approxScale.reshape(num_fns);
101 approxOffset.reshape(num_fns);
102 approxScale = 1.; // default
103 approxOffset = 0.; // default
104 int start_id, end_id;
105 Real scale, offset;
106 while (!scale_offset.eof()) {
107 // Read a record with format: <start id> <end id> <scale> <offset>
108 scale_offset >> start_id >> end_id >> scale >> offset;
109 Cout << "approx_scale_offset.in: " << start_id << ' ' << end_id << ' '
110 << scale << ' ' << offset << '\n';
111 // Check that start_id/end_id are reasonable
112 if ( start_id < 1 || start_id > num_fns || end_id < 1
113 || end_id > num_fns || start_id > end_id) {
114 Cerr << "Error: start/end id's from approx_scale_offset.in do not "
115 << "define a valid range\n within the " << num_fns
116 << " response functions." << std::endl;
117 abort_handler(-1);
118 }
119 // Set scale and offset for the start_id/end_id range
120 for (int i=start_id-1; i<=end_id-1; i++) { // 1-based
121 approxScale[i] = scale;
122 approxOffset[i] = offset;
123 }
124 }
125 }
126 */
127 }
128
129
130 ApproximationInterface::
ApproximationInterface(const String & approx_type,const UShortArray & approx_order,const Variables & am_vars,bool am_cache,const String & am_interface_id,size_t num_fns,short data_order,short output_level)131 ApproximationInterface(const String& approx_type,
132 const UShortArray& approx_order,
133 const Variables& am_vars, bool am_cache,
134 const String& am_interface_id, size_t num_fns,
135 short data_order, short output_level):
136 Interface(NoDBBaseConstructor(), num_fns, output_level), //graph3DFlag(false),
137 challengeFormat(TABULAR_ANNOTATED), challengeActiveOnly(false),
138 actualModelVars(am_vars.copy()),
139 actualModelCache(am_cache), actualModelInterfaceId(am_interface_id)
140 {
141 interfaceId = String("APPROX_INTERFACE_") + std::to_string(++approxIdNum);
142 interfaceType = APPROX_INTERFACE;
143
144 functionSurfaces.resize(num_fns);
145 // despite view mappings, x in map() always = size of active actualModelVars
146 size_t num_vars = actualModelVars.cv() + actualModelVars.div()
147 + actualModelVars.dsv() + actualModelVars.drv();
148 sharedData = SharedApproxData(approx_type, approx_order, num_vars,
149 data_order, output_level);
150 for (int i=0; i<num_fns; i++) {
151 approxFnIndices.insert(i);
152 functionSurfaces[i] = Approximation(sharedData);
153 }
154 }
155
156
157 void ApproximationInterface::
map(const Variables & vars,const ActiveSet & set,Response & response,bool asynch_flag)158 map(const Variables& vars, const ActiveSet& set, Response& response,
159 bool asynch_flag)
160 {
161 ++evalIdCntr; // all calls to map (used throughout as eval id)
162 ++newEvalIdCntr; // nonduplicate evaluations (used ONLY in fn. eval. summary)
163 if (fineGrainEvalCounters) { // detailed evaluation reporting
164 const ShortArray& asv = set.request_vector();
165 size_t i, num_fns = asv.size();
166 init_evaluation_counters(num_fns);
167 for (i=0; i<num_fns; ++i) {
168 short asv_val = asv[i];
169 if (asv_val & 1) { ++fnValCounter[i]; ++newFnValCounter[i]; }
170 if (asv_val & 2) { ++fnGradCounter[i]; ++newFnGradCounter[i]; }
171 if (asv_val & 4) { ++fnHessCounter[i]; ++newFnHessCounter[i]; }
172 }
173 if (fnLabels.empty())
174 fnLabels = response.function_labels();
175 }
176
177 // output the current parameter values prior to running the Analysis
178 if (outputLevel > NORMAL_OUTPUT) {
179 Cout << "\n------------------------------------\n"
180 << "Begin Approximate Fn Evaluation " << std::setw(4) << evalIdCntr;
181 // This may be more confusing than helpful:
182 //if (evalIdRefPt)
183 // Cout << " (local evaluation " << evalIdCntr - evalIdRefPt << ")";
184 Cout << "\n------------------------------------\nParameters for "
185 << "approximate fn evaluation " << evalIdCntr << ":\n" << vars << '\n';
186 }
187 else if (evalIdCntr == 1)
188 Cout << "Beginning Approximate Fn Evaluations..." << std::endl;
189
190 //if (asvControlFlag) // else leave response ActiveSet as initialized
191 response.active_set(set); // set the current ActiveSet within the response
192
193 // Subdivide set for algebraic_mappings() and derived_map()
194 Response algebraic_response, core_response; // empty handles
195 ActiveSet core_set;
196
197 if (algebraicMappings) {
198 if (evalIdCntr == 1)
199 init_algebraic_mappings(vars, response);
200 if (coreMappings) { // both mappings
201 ActiveSet algebraic_set;
202 asv_mapping(set, algebraic_set, core_set);
203 algebraic_response = Response(SIMULATION_RESPONSE, algebraic_set);
204 algebraic_mappings(vars, algebraic_set, algebraic_response);
205 // separate core_response from response
206 core_response = response.copy();
207 core_response.active_set(core_set);
208 }
209 else // algebraic mappings only
210 algebraic_mappings(vars, set, response);
211 }
212 else if (coreMappings) { // analysis_driver mappings only
213 core_set = set;
214 core_response = response; // shared rep
215 }
216
217 if (coreMappings) {
218
219 // Evaluate functionSurfaces at vars and populate core_response.
220 const ShortArray& core_asv = core_set.request_vector();
221 size_t i, num_core_fns = core_asv.size();
222 if ( num_core_fns != functionSurfaces.size() ) {
223 Cerr << "Error: mismatch in number of functions in ApproximationInterface"
224 << "::map()" << std::endl;
225 abort_handler(-1);
226 }
227 // perform any conversions necessary to map between the approximation and
228 // underlying actualModel variables views. Within SurrogateModel,
229 // check_submodel_compatibility() is responsible for verifying that the
230 // views are compatible and force_rebuild() is responsible for verifying
231 // that no unapproximated variables have changed (which would invalidate
232 // the approximation).
233 short approx_active_view = vars.view().first,
234 actual_active_view = actualModelVars.view().first;
235 bool am_vars = true;
236 if (approx_active_view == actual_active_view)
237 am_vars = false;
238 else if ( ( actual_active_view == RELAXED_ALL ||
239 actual_active_view == MIXED_ALL ) &&
240 approx_active_view >= RELAXED_DESIGN ) { // Distinct to All
241 if (vars.acv())
242 actualModelVars.continuous_variables(vars.all_continuous_variables());
243 if (vars.adiv())
244 actualModelVars.discrete_int_variables(
245 vars.all_discrete_int_variables());
246 if (vars.adsv())
247 actualModelVars.discrete_string_variables(
248 vars.all_discrete_string_variables());
249 if (vars.adrv())
250 actualModelVars.discrete_real_variables(
251 vars.all_discrete_real_variables());
252 }
253 else if ( ( approx_active_view == RELAXED_ALL ||
254 approx_active_view == MIXED_ALL ) &&
255 actual_active_view >= RELAXED_DESIGN) { // All to Distinct
256 if (vars.cv())
257 actualModelVars.all_continuous_variables(vars.continuous_variables());
258 if (vars.div())
259 actualModelVars.all_discrete_int_variables(
260 vars.discrete_int_variables());
261 if (vars.dsv())
262 actualModelVars.all_discrete_string_variables(
263 vars.discrete_string_variables());
264 if (vars.drv())
265 actualModelVars.all_discrete_real_variables(
266 vars.discrete_real_variables());
267 }
268 // TO DO: extend for aleatory/epistemic uncertain views
269 else {
270 Cerr << "Error: unsupported variable view differences in "
271 << "ApproximationInterface::map()" << std::endl;
272 abort_handler(-1);
273 }
274 // active subsets of actualModelVars are used in surrogate construction
275 // and evaluation
276 const Variables& surf_vars = (am_vars) ? actualModelVars : vars;
277
278 //size_t num_core_vars = x.length(),
279 //bool approx_scale_len = (approxScale.length()) ? true : false;
280 //bool approx_offset_len = (approxOffset.length()) ? true : false;
281 for (ISIter it=approxFnIndices.begin(); it!=approxFnIndices.end(); ++it) {
282 int index = *it;
283 if (core_asv[index] & 1) {
284 Real fn_val = functionSurfaces[index].value(surf_vars);
285 //if (approx_scale_len)
286 // fn_val *= approxScale[index];
287 //if (approx_offset_len)
288 // fn_val += approxOffset[index];
289 core_response.function_value(fn_val, index);
290 }
291 if (core_asv[index] & 2) { // TO DO: ACV->DVV
292 const RealVector& fn_grad = functionSurfaces[index].gradient(surf_vars);
293 //if (approx_scale_len)
294 // for (size_t j=0; j<num_core_vars; j++)
295 // fn_grad[j] *= approxScale[index];
296 core_response.function_gradient(fn_grad, index);
297 }
298 if (core_asv[index] & 4) { // TO DO: ACV->DVV
299 const RealSymMatrix& fn_hess
300 = functionSurfaces[index].hessian(surf_vars);
301 //if (approx_scale_len)
302 // for (size_t j=0; j<num_core_vars; j++)
303 // for (size_t k=0; k<num_core_vars; k++)
304 // fn_hess[j][k] *= approxScale[index];
305 core_response.function_hessian(fn_hess, index);
306 }
307 }
308 }
309
310 if (algebraicMappings && coreMappings)
311 response_mapping(algebraic_response, core_response, response);
312
313 if (outputLevel > NORMAL_OUTPUT)
314 Cout << "\nActive response data for approximate fn evaluation "
315 << evalIdCntr << ":\n" << response << '\n';
316 if (asynch_flag)
317 beforeSynchResponseMap[evalIdCntr] = response.copy();
318 }
319
320
321 // Little distinction between blocking and nonblocking synch since all
322 // responses are completed.
synchronize()323 const IntResponseMap& ApproximationInterface::synchronize()
324 {
325 // move data from beforeSynch map to completed map
326 rawResponseMap.clear();
327 std::swap(beforeSynchResponseMap, rawResponseMap);
328
329 // augment with any cached evals
330 rawResponseMap.insert(cachedResponseMap.begin(), cachedResponseMap.end());
331 cachedResponseMap.clear();
332
333 return rawResponseMap;
334 }
335
336
synchronize_nowait()337 const IntResponseMap& ApproximationInterface::synchronize_nowait()
338 {
339 // move data from beforeSynch map to completed map
340 rawResponseMap.clear();
341 std::swap(beforeSynchResponseMap, rawResponseMap);
342
343 // augment with any cached evals
344 rawResponseMap.insert(cachedResponseMap.begin(), cachedResponseMap.end());
345 cachedResponseMap.clear();
346
347 return rawResponseMap;
348 }
349
350
351 /** This function populates/replaces each Approximation::anchorPoint
352 with the incoming variables/response data point. */
353 void ApproximationInterface::
update_approximation(const Variables & vars,const IntResponsePair & response_pr)354 update_approximation(const Variables& vars, const IntResponsePair& response_pr)
355 {
356 // NOTE: variable sets passed in from DataFitSurrModel::build_approximation()
357 // correspond to the active continuous variables for either the top level
358 // model or sub-model (DataFitSurrModel::currentVariables or
359 // DataFitSurrModel::actualModel::currentVariables) since the view is the same
360 // in the local case. This can be inconsistent with the use of all continuous
361 // variables above (in constructor and map()) if there are inactive variables
362 // which are changing (e.g. OUU w/ surrogate at UQ level). Currently, a
363 // Taylor series does not compute response derivs w.r.t. inactive variables
364 // and the approximation therefore cannot capture any changes in the inactive
365 // variable values. For this reason, DataFitSurrModel::force_rebuild() forces
366 // a Taylor series rebuild whenever the inactive variable values change.
367
368 // add/replace SurrogateData::anchor{Vars,Resp}
369 if (actualModelCache) {
370 // anchor vars/resp are not sufficiently persistent for use in shallow
371 // copies. Therefore, use ordered id lookup in global PRPCache.
372 IntStringPair ids(response_pr.first, actualModelInterfaceId);
373 PRPCacheCIter p_it;
374 // sign indicates dataset source (see DataFitSurrModel::build_global()):
375 // eval id > 0 for unique evals from current execution (in data_pairs)
376 // eval id = 0 for evals from file import (not in data_pairs)
377 // eval id < 0 for non-unique evals from restart (in data_pairs)
378 if (response_pr.first > 0) // unique evals: current run
379 p_it = lookup_by_ids(data_pairs, ids);
380 else { // non-unique eval ids from restart/file import
381 // rather than resorting to lookup_by_val(), use a two-pass approach
382 // to process multiple returns from equal_range(search_ids)
383 if(actualModelInterfaceId.empty()) {
384 ParamResponsePair search_pr(vars, "NO_ID", response_pr.second);
385 p_it = lookup_by_ids(data_pairs, ids, search_pr);
386 } else {
387 ParamResponsePair search_pr(vars, actualModelInterfaceId,
388 response_pr.second);
389 p_it = lookup_by_ids(data_pairs, ids, search_pr);
390 }
391 }
392 if (p_it == data_pairs.end()) // deep response copies with vars sharing
393 mixed_add(vars, response_pr.second, true);
394 else // shallow copies of cached vars/resp data
395 shallow_add(p_it->variables(), p_it->response(), true);
396 }
397 else // deep response copies with vars sharing
398 mixed_add(vars, response_pr.second, true);
399
400 // reset active approxData key using SharedApproxData::approxDataKeys
401 restore_data_key();
402 }
403
404
405 /** This function populates/replaces each Approximation::currentPoints
406 with the incoming variables/response arrays. */
407 void ApproximationInterface::
update_approximation(const RealMatrix & samples,const IntResponseMap & resp_map)408 update_approximation(const RealMatrix& samples, const IntResponseMap& resp_map)
409 {
410 size_t i, num_pts = resp_map.size();
411 if (samples.numCols() != num_pts) {
412 Cerr << "Error: mismatch in variable and response set lengths in "
413 << "ApproximationInterface::update_approximation()." << std::endl;
414 abort_handler(-1);
415 }
416 // clear active SurrogateData::{vars,resp}Data
417 ISIter a_it; IntRespMCIter r_it;
418 for (a_it=approxFnIndices.begin(); a_it!=approxFnIndices.end(); ++a_it)
419 functionSurfaces[*a_it].clear_active_data();
420 // replace active SurrogateData::{vars,resp}Data
421 if (actualModelCache) {
422 PRPCacheCIter p_it; size_t num_cv = samples.numRows();
423 for (i=0, r_it=resp_map.begin(); i<num_pts; ++i, ++r_it) {
424 // allVariables/allResponses are not sufficiently persistent for use in
425 // shallow copies. Therefore, use ordered id lookup in global PRPCache.
426 IntStringPair ids(r_it->first, actualModelInterfaceId);
427 // sign indicates dataset source (see DataFitSurrModel::build_global()):
428 // eval id > 0 for unique evals from current execution (in data_pairs)
429 // eval id = 0 for evals from file import (not in data_pairs)
430 // eval id < 0 for non-unique evals from restart (in data_pairs)
431 if (r_it->first > 0) // unique evals: current run
432 p_it = lookup_by_ids(data_pairs, ids);
433 else { // nonunique eval ids from restart/file import
434 // rather than resorting to lookup_by_val(), use a two-pass approach
435 // to process multiple returns from equal_range(search_ids)
436 sample_to_variables(samples[i], num_cv, actualModelVars);
437 if (actualModelInterfaceId.empty()) {
438 ParamResponsePair search_pr(actualModelVars, "NO_ID", r_it->second);
439 p_it = lookup_by_ids(data_pairs, ids, search_pr);
440 } else {
441 ParamResponsePair search_pr(actualModelVars, actualModelInterfaceId,
442 r_it->second);
443 p_it = lookup_by_ids(data_pairs, ids, search_pr);
444 }
445 }
446 if (p_it == data_pairs.end()) // deep response copies with vars sharing
447 mixed_add(samples[i], r_it->second, false);
448 else // shallow copies of cached vars/resp data
449 shallow_add(p_it->variables(), p_it->response(), false);
450 }
451 }
452 else // deep response copies with vars sharing
453 for (i=0, r_it=resp_map.begin(); i<num_pts; ++i, ++r_it)
454 mixed_add(samples[i], r_it->second, false);
455
456 // reset active approxData key using SharedApproxData::approxDataKeys
457 restore_data_key();
458 }
459
460
461 /** This function populates/replaces each Approximation::currentPoints
462 with the incoming variables/response arrays. */
463 void ApproximationInterface::
update_approximation(const VariablesArray & vars_array,const IntResponseMap & resp_map)464 update_approximation(const VariablesArray& vars_array,
465 const IntResponseMap& resp_map)
466 {
467 size_t i, num_pts = resp_map.size();
468 if (vars_array.size() != num_pts) {
469 Cerr << "Error: mismatch in variable and response set lengths in "
470 << "ApproximationInterface::update_approximation()." << std::endl;
471 abort_handler(-1);
472 }
473 // clear active SurrogateData::{vars,resp}Data
474 ISIter a_it; IntRespMCIter r_it;
475 for (a_it=approxFnIndices.begin(); a_it!=approxFnIndices.end(); ++a_it)
476 functionSurfaces[*a_it].clear_active_data();
477 // replace active SurrogateData::{vars,resp}Data
478 if (actualModelCache) {
479 PRPCacheCIter p_it;
480 for (i=0, r_it=resp_map.begin(); i<num_pts; ++i, ++r_it) {
481 // allVariables/allResponses are not sufficiently persistent for use in
482 // shallow copies. Therefore, use ordered id lookup in global PRPCache.
483 IntStringPair ids(r_it->first, actualModelInterfaceId);
484 // sign indicates dataset source (see DataFitSurrModel::build_global()):
485 // eval id > 0 for unique evals from current execution (in data_pairs)
486 // eval id = 0 for evals from file import (not in data_pairs)
487 // eval id < 0 for non-unique evals from restart (in data_pairs)
488 if (r_it->first > 0) // unique evals: current run
489 p_it = lookup_by_ids(data_pairs, ids);
490 else { // nonunique eval ids from restart/file import
491 // rather than resorting to lookup_by_val(), use a two-pass approach
492 // to process multiple returns from equal_range(search_ids)
493 if(actualModelInterfaceId.empty()) {
494 ParamResponsePair search_pr(vars_array[i], "NO_ID", r_it->second);
495 p_it = lookup_by_ids(data_pairs, ids, search_pr);
496 } else {
497 ParamResponsePair search_pr(vars_array[i], actualModelInterfaceId,
498 r_it->second);
499 p_it = lookup_by_ids(data_pairs, ids, search_pr);
500 }
501 }
502 if (p_it == data_pairs.end()) // deep response copies with vars sharing
503 mixed_add(vars_array[i], r_it->second, false);
504 else // shallow copies of cached vars/resp data
505 shallow_add(p_it->variables(), p_it->response(), false);
506 }
507 }
508 else // deep response copies with vars sharing
509 for (i=0, r_it=resp_map.begin(); i<num_pts; ++i, ++r_it)
510 mixed_add(vars_array[i], r_it->second, false);
511
512 // reset active approxData key using SharedApproxData::approxDataKeys
513 restore_data_key();
514 }
515
516
517 /** This function appends to each Approximation::currentPoints with
518 one incoming variables/response data point. */
519 void ApproximationInterface::
append_approximation(const Variables & vars,const IntResponsePair & response_pr)520 append_approximation(const Variables& vars, const IntResponsePair& response_pr)
521 {
522 // append a single point to SurrogateData::{vars,resp}Data
523 if (actualModelCache) {
524 // anchor vars/resp are not sufficiently persistent for use in shallow
525 // copies. Therefore, use ordered id lookup in global PRPCache.
526 IntStringPair ids(response_pr.first, actualModelInterfaceId);
527 PRPCacheCIter p_it;
528 // sign indicates dataset source (see DataFitSurrModel::build_global()):
529 // eval id > 0 for unique evals from current execution (in data_pairs)
530 // eval id = 0 for evals from file import (not in data_pairs)
531 // eval id < 0 for non-unique evals from restart (in data_pairs)
532 if (response_pr.first > 0) // unique evals: current run
533 p_it = lookup_by_ids(data_pairs, ids);
534 else { // nonunique eval ids from restart/file import
535 // rather than resorting to lookup_by_val(), use a two-pass approach
536 // to process multiple returns from equal_range(search_ids)
537 if(actualModelInterfaceId.empty()) {
538 ParamResponsePair search_pr(vars, "NO_ID", response_pr.second);
539 p_it = lookup_by_ids(data_pairs, ids, search_pr);
540 } else {
541 ParamResponsePair search_pr(vars, actualModelInterfaceId,
542 response_pr.second);
543 p_it = lookup_by_ids(data_pairs, ids, search_pr);
544 }
545 }
546 if (p_it == data_pairs.end()) // deep response copies with vars sharing
547 mixed_add(vars, response_pr.second, false);
548 else // shallow copies of cached vars/resp data
549 shallow_add(p_it->variables(), p_it->response(), false);
550 }
551 else // deep response copies with vars sharing
552 mixed_add(vars, response_pr.second, false);
553
554 update_pop_counts(response_pr);
555
556 // reset active approxData key using SharedApproxData::approxDataKeys
557 restore_data_key();
558 }
559
560
561 /** This function appends to each Approximation::currentPoints with
562 multiple incoming variables/response data points. */
563 void ApproximationInterface::
append_approximation(const RealMatrix & samples,const IntResponseMap & resp_map)564 append_approximation(const RealMatrix& samples, const IntResponseMap& resp_map)
565 {
566 size_t i, num_pts = resp_map.size();
567 if (samples.numCols() != num_pts) {
568 Cerr << "Error: mismatch in variable and response set lengths in "
569 << "ApproximationInterface::append_approximation()." << std::endl;
570 abort_handler(-1);
571 }
572 // append multiple points to SurrogateData::{vars,resp}Data
573 IntRespMCIter r_it;
574 if (actualModelCache) {
575 PRPCacheCIter p_it; size_t num_cv = samples.numRows();
576 for (i=0, r_it=resp_map.begin(); i<num_pts; ++i, ++r_it) {
577 // allVariables/allResponses are not sufficiently persistent for use in
578 // shallow copies. Therefore, use ordered id lookup in global PRPCache.
579 IntStringPair ids(r_it->first, actualModelInterfaceId);
580 // sign indicates dataset source (see DataFitSurrModel::build_global()):
581 // eval id > 0 for unique evals from current execution (in data_pairs)
582 // eval id = 0 for evals from file import (not in data_pairs)
583 // eval id < 0 for non-unique evals from restart (in data_pairs)
584 if (r_it->first > 0) // unique evals: current run
585 p_it = lookup_by_ids(data_pairs, ids);
586 else { // nonunique eval ids from restart/file import
587 // rather than resorting to lookup_by_val(), use a two-pass approach
588 // to process multiple returns from equal_range(search_ids)
589 sample_to_variables(samples[i], num_cv, actualModelVars);
590 if(actualModelInterfaceId.empty()) {
591 ParamResponsePair search_pr(actualModelVars, "NO_ID", r_it->second);
592 p_it = lookup_by_ids(data_pairs, ids, search_pr);
593 } else {
594 ParamResponsePair search_pr(actualModelVars, actualModelInterfaceId,
595 r_it->second);
596 p_it = lookup_by_ids(data_pairs, ids, search_pr);
597 }
598 }
599 if (p_it == data_pairs.end()) // deep response copies with vars sharing
600 mixed_add(samples[i], r_it->second, false);
601 else // shallow copies of cached vars/resp data
602 shallow_add(p_it->variables(), p_it->response(), false);
603 }
604 }
605 else // deep response copies with vars sharing
606 for (i=0, r_it=resp_map.begin(); i<num_pts; ++i, ++r_it)
607 mixed_add(samples[i], r_it->second, false);
608
609 update_pop_counts(resp_map);
610
611 // reset active approxData key using SharedApproxData::approxDataKeys
612 restore_data_key();
613 }
614
615
616 /** This function appends to each Approximation::currentPoints with
617 multiple incoming variables/response data points. */
618 void ApproximationInterface::
append_approximation(const VariablesArray & vars_array,const IntResponseMap & resp_map)619 append_approximation(const VariablesArray& vars_array,
620 const IntResponseMap& resp_map)
621 {
622 size_t i, num_pts = resp_map.size();
623 if (vars_array.size() != num_pts) {
624 Cerr << "Error: mismatch in variable and response set lengths in "
625 << "ApproximationInterface::append_approximation()." << std::endl;
626 abort_handler(-1);
627 }
628 // append multiple points to SurrogateData::{vars,resp}Data
629 IntRespMCIter r_it;
630 if (actualModelCache) {
631 PRPCacheCIter p_it;
632 for (i=0, r_it=resp_map.begin(); i<num_pts; ++i, ++r_it) {
633 // allVariables/allResponses are not sufficiently persistent for use in
634 // shallow copies. Therefore, use ordered id lookup in global PRPCache.
635 IntStringPair ids(r_it->first, actualModelInterfaceId);
636 // sign indicates dataset source (see DataFitSurrModel::build_global()):
637 // eval id > 0 for unique evals from current execution (in data_pairs)
638 // eval id = 0 for evals from file import (not in data_pairs)
639 // eval id < 0 for non-unique evals from restart (in data_pairs)
640 if (r_it->first > 0) // unique evals: current run
641 p_it = lookup_by_ids(data_pairs, ids);
642 else { // nonunique eval ids from restart/file import
643 // rather than resorting to lookup_by_val(), use a two-pass approach
644 // to process multiple returns from equal_range(search_ids)
645 if(actualModelInterfaceId.empty()) {
646 ParamResponsePair search_pr(vars_array[i], "NO_ID", r_it->second);
647 p_it = lookup_by_ids(data_pairs, ids, search_pr);
648 } else {
649 ParamResponsePair search_pr(vars_array[i], actualModelInterfaceId,
650 r_it->second);
651 p_it = lookup_by_ids(data_pairs, ids, search_pr);
652 }
653 }
654 if (p_it == data_pairs.end()) // deep response copies with vars sharing
655 mixed_add(vars_array[i], r_it->second, false);
656 else // shallow copies of cached vars/resp data
657 shallow_add(p_it->variables(), p_it->response(), false);
658 }
659 }
660 else // deep response copies with vars sharing
661 for (i=0, r_it=resp_map.begin(); i<num_pts; ++i, ++r_it)
662 mixed_add(vars_array[i], r_it->second, false);
663
664 update_pop_counts(resp_map);
665
666 // reset active approxData key using SharedApproxData::approxDataKeys
667 restore_data_key();
668 }
669
670
671 /** This function finds the coefficients for each Approximation based
672 on the data passed through update_approximation() calls. The
673 bounds are used only for graphics visualization. */
674 void ApproximationInterface::
build_approximation(const RealVector & c_l_bnds,const RealVector & c_u_bnds,const IntVector & di_l_bnds,const IntVector & di_u_bnds,const RealVector & dr_l_bnds,const RealVector & dr_u_bnds)675 build_approximation(const RealVector& c_l_bnds, const RealVector& c_u_bnds,
676 const IntVector& di_l_bnds, const IntVector& di_u_bnds,
677 const RealVector& dr_l_bnds, const RealVector& dr_u_bnds)
678 {
679 // initialize the data shared among approximation instances
680 sharedData.set_bounds(c_l_bnds, c_u_bnds, di_l_bnds, di_u_bnds,
681 dr_l_bnds, dr_u_bnds);
682 sharedData.build();
683 // build the approximation surface instances
684 for (ISIter it=approxFnIndices.begin(); it!=approxFnIndices.end(); ++it) {
685 int fn_index = *it;
686 // construct the approximation
687 functionSurfaces[fn_index].build();
688
689 // manage diagnostics
690 if (functionSurfaces[fn_index].diagnostics_available()) {
691 // print default or user-requested metrics and cross-validation
692 functionSurfaces[fn_index].primary_diagnostics(fn_index);
693 // for user-provided challenge data, we assume there are
694 // function values for all functions in the analysis, not just
695 // the indices for which surrogates are being built
696
697 // BMA TODO: can this move to ctor?
698 if (!challengeFile.empty()) {
699 if (challengePoints.empty())
700 read_challenge_points();
701 functionSurfaces[fn_index].challenge_diagnostics(fn_index,
702 challengePoints, Teuchos::getCol(Teuchos::View, challengeResponses,
703 fn_index));
704 }
705 }
706 }
707
708 /* Old 3D graphics capability:
709 int fn_index = *approxFnIndices.begin();
710 // if graphics is on for 2 variables, plot first functionSurface in 3D
711 if (graph3DFlag && sharedData.num_variables() == 2) {
712 functionSurfaces[fn_index].draw_surface();
713 }
714 */
715 }
716
717
718 /** This function calls export on each approximation */
export_approximation()719 void ApproximationInterface::export_approximation()
720 {
721 for (ISIter it=approxFnIndices.begin(); it!=approxFnIndices.end(); ++it)
722 functionSurfaces[*it].export_model();
723 }
724
725
726 /** This function updates the coefficients for each Approximation based
727 on data increments provided by {update,append}_approximation(). */
rebuild_approximation(const BitArray & rebuild_fns)728 void ApproximationInterface::rebuild_approximation(const BitArray& rebuild_fns)
729 {
730 // rebuild data shared among approximation instances
731 sharedData.rebuild();
732 // rebuild the approximation surfaces
733 for (ISIter it=approxFnIndices.begin(); it!=approxFnIndices.end(); ++it) {
734 int fn_index = *it;
735 // check for rebuild request (defaults to true if no deque defined)
736 if (rebuild_fns.empty() || rebuild_fns[fn_index]) {
737 // approx bounds not updated as in build_approximation()
738
739 // invokes increment_coefficients()
740 functionSurfaces[fn_index].rebuild();
741
742 // diagnostics not currently active on rebuild
743 }
744 }
745 }
746
747
748 void ApproximationInterface::
mixed_add(const Variables & vars,const Response & response,bool anchor)749 mixed_add(const Variables& vars, const Response& response, bool anchor)
750 {
751 Pecos::SurrogateDataVars sdv; bool first_vars = true;
752 const ShortArray& asv = response.active_set_request_vector();
753 size_t i, fn_index, num_fns = functionSurfaces.size(), num_asv = asv.size(),
754 key_index;
755 for (ISIter it=approxFnIndices.begin(); it!=approxFnIndices.end(); ++it) {
756 fn_index = *it;
757 // asv may be larger than num_fns due to response aggregation modes
758 // (e.g., multifidelity) --> use num_fns as stride and support add()
759 // to vector of SurrogateData
760 for (i=fn_index, key_index=0; i<num_asv; i+=num_fns, ++key_index)
761 if (asv[i]) {
762 Approximation& fn_surf = functionSurfaces[fn_index];
763 // rather than unrolling the response (containing all response fns)
764 // into per-response function arrays for input to fn_surf, pass the
765 // complete response along with a response function index.
766 if (first_vars) {
767 fn_surf.add(vars, anchor, true, key_index); // deep
768 fn_surf.add(response, i, anchor, true, key_index); // deep
769 // carry newly added sdv over to other approx fn indices:
770 sdv = (anchor) ? fn_surf.surrogate_data().anchor_variables() :
771 fn_surf.surrogate_data().variables_data().back();
772 first_vars = false;
773 }
774 else {
775 fn_surf.add(sdv, anchor, false, key_index); // shallow
776 fn_surf.add(response, i, anchor, true, key_index); // deep
777 }
778 }
779 }
780 }
781
782
783 void ApproximationInterface::
mixed_add(const Real * c_vars,const Response & response,bool anchor)784 mixed_add(const Real* c_vars, const Response& response, bool anchor)
785 {
786 Pecos::SurrogateDataVars sdv; bool first_vars = true;
787 const ShortArray& asv = response.active_set_request_vector();
788 size_t i, fn_index, num_fns = functionSurfaces.size(), num_asv = asv.size(),
789 key_index;
790 for (ISIter it=approxFnIndices.begin(); it!=approxFnIndices.end(); ++it) {
791 fn_index = *it;
792 // asv may be larger than num_fns due to response aggregation modes
793 // (e.g., multifidelity) --> use num_fns as stride and support add()
794 // to vector of SurrogateData
795 for (i=fn_index, key_index=0; i<num_asv; i+=num_fns, ++key_index)
796 if (asv[i]) {
797 Approximation& fn_surf = functionSurfaces[fn_index];
798 // rather than unrolling the response (containing all response fns)
799 // into per-response function arrays for input to fn_surf, pass the
800 // complete response along with a response function index.
801 if (first_vars) {
802 fn_surf.add(c_vars, anchor, true, key_index); // deep
803 fn_surf.add(response, i, anchor, true, key_index); // deep
804 // carry newly added sdv over to other approx fn indices:
805 sdv = (anchor) ? fn_surf.surrogate_data().anchor_variables() :
806 fn_surf.surrogate_data().variables_data().back();
807 first_vars = false;
808 }
809 else {
810 fn_surf.add(sdv, anchor, false, key_index); // shallow
811 fn_surf.add(response, i, anchor, true, key_index); // deep
812 }
813 }
814 }
815 }
816
817
818 void ApproximationInterface::
shallow_add(const Variables & vars,const Response & response,bool anchor)819 shallow_add(const Variables& vars, const Response& response, bool anchor)
820 {
821 const ShortArray& asv = response.active_set_request_vector();
822 size_t i, fn_index, num_fns = functionSurfaces.size(), num_asv = asv.size(),
823 key_index;
824 for (ISIter it=approxFnIndices.begin(); it!=approxFnIndices.end(); ++it) {
825 fn_index = *it;
826 // asv may be larger than num_fns due to response aggregation modes
827 // (e.g., multifidelity) --> use num_fns as stride and support add()
828 // to vector of SurrogateData
829 for (i=fn_index, key_index=0; i<num_asv; i+=num_fns, ++key_index)
830 if (asv[i]) {
831 Approximation& fn_surf = functionSurfaces[fn_index];
832 fn_surf.add(vars, anchor, false, key_index); // shallow
833 fn_surf.add(response, i, anchor, false, key_index); // shallow
834 }
835 }
836 }
837
838
839 void ApproximationInterface::
update_pop_counts(const IntResponsePair & response_pr)840 update_pop_counts(const IntResponsePair& response_pr)
841 {
842 // update pop counts for data in response_pr
843 const ShortArray& asv = response_pr.second.active_set_request_vector();
844 size_t i, key_index, fn_index, num_fns = functionSurfaces.size(),
845 num_asv = asv.size(), count;
846 ISIter fn_it;
847 for (fn_it=approxFnIndices.begin(); fn_it!=approxFnIndices.end(); ++fn_it) {
848 fn_index = *fn_it;
849 // asv may be larger than num_fns due to response aggregation modes
850 for (i=fn_index, key_index=0; i<num_asv; i+=num_fns, ++key_index) {
851 count = (asv[i]) ? 1 : 0; // either one or no pts appended
852 functionSurfaces[fn_index].pop_count(count, key_index);
853 }
854 }
855 }
856
857
update_pop_counts(const IntResponseMap & resp_map)858 void ApproximationInterface::update_pop_counts(const IntResponseMap& resp_map)
859 {
860 ISIter fn_it; IntRespMCIter r_it = resp_map.begin();
861 size_t i, count, key_index, fn_index, num_fns = functionSurfaces.size(),
862 num_asv = r_it->second.active_set_request_vector().size();
863 for (fn_it=approxFnIndices.begin(); fn_it!=approxFnIndices.end(); ++fn_it) {
864 fn_index = *fn_it;
865 // asv may be larger than num_fns due to response aggregation modes
866 for (i=fn_index, key_index=0; i<num_asv; i+=num_fns, ++key_index) {
867 for (r_it=resp_map.begin(), count=0; r_it!=resp_map.end(); ++r_it)
868 if (r_it->second.active_set_request_vector()[i])
869 ++count;
870 functionSurfaces[fn_index].pop_count(count, key_index);
871 }
872 }
873 }
874
875
876 Real2DArray ApproximationInterface::
cv_diagnostics(const StringArray & metric_types,unsigned num_folds)877 cv_diagnostics(const StringArray& metric_types, unsigned num_folds)
878 {
879 Real2DArray cv_diags;
880 ISIter a_it = approxFnIndices.begin();
881 ISCIter a_end = approxFnIndices.end();
882 for ( ; a_it != a_end; ++a_it) {
883 size_t index = *a_it;
884 cv_diags.push_back(functionSurfaces[index].
885 cv_diagnostic(metric_types, num_folds));
886 }
887 return cv_diags;
888 }
889
890
891 Real2DArray ApproximationInterface::
challenge_diagnostics(const StringArray & metric_types,const RealMatrix & challenge_pts,const RealVector & challenge_resps)892 challenge_diagnostics(const StringArray& metric_types,
893 const RealMatrix& challenge_pts,
894 const RealVector& challenge_resps)
895 {
896 Real2DArray chall_diags;
897 ISIter a_it = approxFnIndices.begin();
898 ISCIter a_end = approxFnIndices.end();
899 for ( ; a_it != a_end; ++a_it) {
900 size_t index = *a_it;
901 chall_diags.push_back(functionSurfaces[index].
902 challenge_diagnostic(metric_types, challenge_pts,
903 challenge_resps));
904 }
905 return chall_diags;
906 }
907
908
909
910 const RealVectorArray& ApproximationInterface::
approximation_coefficients(bool normalized)911 approximation_coefficients(bool normalized)
912 {
913 // only assign the functionSurfaceCoeffs array if it's requested
914 // (i.e., do it here rather than in build/update functions above).
915 if (functionSurfaceCoeffs.empty())
916 functionSurfaceCoeffs.resize(functionSurfaces.size());
917 for (ISIter it=approxFnIndices.begin(); it!=approxFnIndices.end(); ++it) {
918 int index = *it;
919 functionSurfaceCoeffs[index]
920 = functionSurfaces[index].approximation_coefficients(normalized);
921 }
922 return functionSurfaceCoeffs;
923 }
924
925
926 void ApproximationInterface::
approximation_coefficients(const RealVectorArray & approx_coeffs,bool normalized)927 approximation_coefficients(const RealVectorArray& approx_coeffs,
928 bool normalized)
929 {
930 for (ISIter it=approxFnIndices.begin(); it!=approxFnIndices.end(); ++it) {
931 int index = *it;
932 functionSurfaces[index].approximation_coefficients(approx_coeffs[index],
933 normalized);
934 }
935 //functionSurfaceCoeffs = approx_coeffs;
936 }
937
938
939 const RealVector& ApproximationInterface::
approximation_variances(const Variables & vars)940 approximation_variances(const Variables& vars)
941 {
942 // only assign the functionSurfaceVariances array if it's requested
943 // (i.e., do it here rather than in build/update functions above).
944 if (functionSurfaceVariances.empty())
945 functionSurfaceVariances.sizeUninitialized(functionSurfaces.size());
946 for (ISIter it=approxFnIndices.begin(); it!=approxFnIndices.end(); ++it) {
947 int index = *it;
948 functionSurfaceVariances[index]
949 = functionSurfaces[index].prediction_variance(vars);
950 }
951 return functionSurfaceVariances;
952 }
953
954
955 // TODO: What does it even mean to challenge at index or string
956 // data?!? Is a Response object available too?
957 /** Challenge data defaults to active/inactive, but user can override
958 to active only. */
read_challenge_points()959 void ApproximationInterface::read_challenge_points()
960 {
961 size_t num_fns = functionSurfaces.size();
962 String context_msg = "Surrogate model, interface id '" + interface_id() +
963 "' import_challenge_points_file";
964 bool verbose = (outputLevel > NORMAL_OUTPUT);
965 // use a Variables object to easily read active vs. all
966 TabularIO::read_data_tabular(challengeFile, context_msg,
967 actualModelVars.copy(), num_fns, challengePoints,
968 challengeResponses, challengeFormat,
969 verbose, challengeUseVarLabels,
970 challengeActiveOnly);
971 }
972
973 } // namespace Dakota
974