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: TestDriverInterface
10 //- Description: Class implementation
11 //- Owner: Mike Eldred, Brian Adams
12
13 #include "TestDriverInterface.hpp"
14 #include "ParallelLibrary.hpp"
15 #include "DataMethod.hpp" // for output levels
16 //#include <thread> // for sleep_for
17 #ifdef DAKOTA_MODELCENTER
18 #include "PHXCppApi.h"
19 #endif
20 #include <boost/tokenizer.hpp>
21 #include "dakota_mersenne_twister.hpp"
22 #include <boost/random/normal_distribution.hpp>
23 #include <boost/assign.hpp>
24 #include <vector>
25 #include "Teuchos_SerialDenseHelpers.hpp"
26 #include "NonDLHSSampling.hpp"
27 #include "spectral_diffusion.hpp"
28 #include "predator_prey.hpp"
29
30 namespace Dakota {
31
32 /// offset used text_book exponent: 1.0 is nominal, 1.4 used for B&B testing
33 const double POW_VAL = 1.0 ;
34 /// levenshtein_distance computes the distance between its argument and this
35 // reference string
36 const String LEV_REF = "Dakota";
37
38 StringRealMap TestDriverInterface::levenshteinDistanceCache;
39
40 #ifdef DAKOTA_SALINAS
41 /// subroutine interface to SALINAS simulation code
42 int salinas_main(int argc, char *argv[], MPI_Comm* comm);
43 #endif // DAKOTA_SALINAS
44
45
TestDriverInterface(const ProblemDescDB & problem_db)46 TestDriverInterface::TestDriverInterface(const ProblemDescDB& problem_db)
47 : DirectApplicInterface(problem_db)
48 {
49 // register this class' analysis driver types with the string to enum map
50 // at the base class
51 driverTypeMap["cantilever"] = CANTILEVER_BEAM;
52 driverTypeMap["mod_cantilever"] = MOD_CANTILEVER_BEAM;
53 driverTypeMap["cantilever_ml"] = CANTILEVER_BEAM_ML;
54 driverTypeMap["cyl_head"] = CYLINDER_HEAD;
55 driverTypeMap["extended_rosenbrock"] = EXTENDED_ROSENBROCK;
56 driverTypeMap["generalized_rosenbrock"] = GENERALIZED_ROSENBROCK;
57 driverTypeMap["lf_rosenbrock"] = LF_ROSENBROCK;
58 driverTypeMap["extra_lf_rosenbrock"] = EXTRA_LF_ROSENBROCK;
59 driverTypeMap["mf_rosenbrock"] = MF_ROSENBROCK;
60 driverTypeMap["rosenbrock"] = ROSENBROCK;
61 driverTypeMap["modified_rosenbrock"] = MODIFIED_ROSENBROCK;
62 driverTypeMap["lf_poly_prod"] = LF_POLY_PROD;
63 driverTypeMap["poly_prod"] = POLY_PROD;
64 driverTypeMap["gerstner"] = GERSTNER;
65 driverTypeMap["scalable_gerstner"] = SCALABLE_GERSTNER;
66 driverTypeMap["log_ratio"] = LOGNORMAL_RATIO;
67 driverTypeMap["multimodal"] = MULTIMODAL;
68 driverTypeMap["lf_short_column"] = LF_SHORT_COLUMN;
69 driverTypeMap["mf_short_column"] = MF_SHORT_COLUMN;
70 driverTypeMap["short_column"] = SHORT_COLUMN;
71 driverTypeMap["side_impact_cost"] = SIDE_IMPACT_COST;
72 driverTypeMap["side_impact_perf"] = SIDE_IMPACT_PERFORMANCE;
73 driverTypeMap["sobol_rational"] = SOBOL_RATIONAL;
74 driverTypeMap["sobol_g_function"] = SOBOL_G_FUNCTION;
75 driverTypeMap["sobol_ishigami"] = SOBOL_ISHIGAMI;
76 driverTypeMap["steel_column_cost"] = STEEL_COLUMN_COST;
77 driverTypeMap["steel_column_perf"] = STEEL_COLUMN_PERFORMANCE;
78 driverTypeMap["text_book"] = TEXT_BOOK;
79 driverTypeMap["text_book1"] = TEXT_BOOK1;
80 driverTypeMap["text_book2"] = TEXT_BOOK2;
81 driverTypeMap["text_book3"] = TEXT_BOOK3;
82 driverTypeMap["text_book_ouu"] = TEXT_BOOK_OUU;
83 driverTypeMap["scalable_text_book"] = SCALABLE_TEXT_BOOK;
84 driverTypeMap["scalable_monomials"] = SCALABLE_MONOMIALS;
85 driverTypeMap["mogatest1"] = MOGATEST1;
86 driverTypeMap["mogatest2"] = MOGATEST2;
87 driverTypeMap["mogatest3"] = MOGATEST3;
88 driverTypeMap["illumination"] = ILLUMINATION;
89 driverTypeMap["barnes"] = BARNES;
90 driverTypeMap["barnes_lf"] = BARNES_LF;
91 driverTypeMap["herbie"] = HERBIE;
92 driverTypeMap["smooth_herbie"] = SMOOTH_HERBIE;
93 driverTypeMap["shubert"] = SHUBERT;
94 driverTypeMap["salinas"] = SALINAS;
95 driverTypeMap["mc_api_run"] = MODELCENTER;
96 driverTypeMap["modelcenter"] = MODELCENTER;
97 driverTypeMap["genz"] = GENZ;
98 driverTypeMap["damped_oscillator"] = DAMPED_OSCILLATOR;
99 driverTypeMap["steady_state_diffusion_1d"] = STEADY_STATE_DIFFUSION_1D;
100 driverTypeMap["ss_diffusion_discrepancy"] = SS_DIFFUSION_DISCREPANCY;
101 driverTypeMap["transient_diffusion_1d"] = TRANSIENT_DIFFUSION_1D;
102 driverTypeMap["predator_prey"] = PREDATOR_PREY;
103 driverTypeMap["aniso_quad_form"] = ANISOTROPIC_QUADRATIC_FORM;
104 driverTypeMap["bayes_linear"] = BAYES_LINEAR;
105 driverTypeMap["problem18"] = PROBLEM18;
106
107 // convert strings to enums for analysisDriverTypes, iFilterType, oFilterType
108 analysisDriverTypes.resize(numAnalysisDrivers);
109 std::map<String, driver_t>::iterator sd_iter;
110 for (size_t i=0; i<numAnalysisDrivers; ++i) {
111 sd_iter = driverTypeMap.find(analysisDrivers[i]);//toLower(Drivers[i]));
112 if (sd_iter == driverTypeMap.end()) {
113 if (outputLevel > NORMAL_OUTPUT)
114 Cerr << "Warning: analysis_driver \"" << analysisDrivers[i] << "\" not "
115 << "available at construct time in TestDriverInterface.\n "
116 << " Subsequent interface plug-in may resolve." << std::endl;
117 analysisDriverTypes[i] = NO_DRIVER;
118 }
119 else
120 analysisDriverTypes[i] = sd_iter->second;
121 }
122
123 sd_iter = driverTypeMap.find(iFilterName); //toLower(iFilterName));
124 if (sd_iter == driverTypeMap.end()) {
125 if (outputLevel > NORMAL_OUTPUT)
126 Cerr << "Warning: input filter \"" << iFilterName << "\" not available at"
127 << " construct time in TestDriverInterface.\n Subsequent "
128 << "interface plug-in may resolve." << std::endl;
129 iFilterType = NO_DRIVER;
130 }
131 else
132 iFilterType = sd_iter->second;
133
134 sd_iter = driverTypeMap.find(oFilterName); //toLower(oFilterName));
135 if (sd_iter == driverTypeMap.end()) {
136 if (outputLevel > NORMAL_OUTPUT)
137 Cerr << "Warning: output filter \"" << oFilterName << "\" not available "
138 << "at construct time in TestDriverInterface.\n Subsequent"
139 << " interface plug-in may resolve." << std::endl;
140 oFilterType = NO_DRIVER;
141 }
142 else
143 oFilterType = sd_iter->second;
144
145 // define localDataView from analysisDriverTypes,
146 // overriding any base class constructor setting
147 localDataView = 0;
148 for (size_t i=0; i<numAnalysisDrivers; ++i)
149 switch (analysisDriverTypes[i]) {
150 case CANTILEVER_BEAM: case MOD_CANTILEVER_BEAM: case CANTILEVER_BEAM_ML:
151 case ROSENBROCK: case LF_ROSENBROCK: case EXTRA_LF_ROSENBROCK:
152 case MF_ROSENBROCK: case MODIFIED_ROSENBROCK: case PROBLEM18:
153 case SHORT_COLUMN: case LF_SHORT_COLUMN: case MF_SHORT_COLUMN:
154 case SOBOL_ISHIGAMI: case STEEL_COLUMN_COST: case STEEL_COLUMN_PERFORMANCE:
155 localDataView |= VARIABLES_MAP; break;
156 case NO_DRIVER: // assume VARIABLES_VECTOR approach for plug-ins for now
157 case CYLINDER_HEAD: case LOGNORMAL_RATIO: case MULTIMODAL:
158 case GERSTNER: case SCALABLE_GERSTNER:
159 case EXTENDED_ROSENBROCK: case GENERALIZED_ROSENBROCK:
160 case SIDE_IMPACT_COST: case SIDE_IMPACT_PERFORMANCE:
161 case SOBOL_G_FUNCTION: case SOBOL_RATIONAL:
162 case TEXT_BOOK: case TEXT_BOOK1: case TEXT_BOOK2: case TEXT_BOOK3:
163 case TEXT_BOOK_OUU: case SCALABLE_TEXT_BOOK: case SCALABLE_MONOMIALS:
164 case MOGATEST1: case MOGATEST2: case MOGATEST3:
165 case ILLUMINATION: case LF_POLY_PROD: case POLY_PROD:
166 case BARNES: case BARNES_LF:
167 case HERBIE: case SMOOTH_HERBIE: case SHUBERT:
168 case SALINAS: case MODELCENTER:
169 case GENZ: case DAMPED_OSCILLATOR:
170 case STEADY_STATE_DIFFUSION_1D: case SS_DIFFUSION_DISCREPANCY:
171 case TRANSIENT_DIFFUSION_1D: case PREDATOR_PREY:
172 case ANISOTROPIC_QUADRATIC_FORM: case BAYES_LINEAR:
173 localDataView |= VARIABLES_VECTOR; break;
174 }
175
176 // define varTypeMap for analysis drivers based on xCM/XDM maps
177 if (localDataView & VARIABLES_MAP) {
178 // define the string to enumeration map
179 //switch (ac_name) {
180 //case ROSENBROCK: case LF_ROSENBROCK: case MF_ROSENBROCK:
181 //case SOBOL_ISHIGAMI:
182 varTypeMap["x1"] = VAR_x1; varTypeMap["x2"] = VAR_x2;
183 varTypeMap["x3"] = VAR_x3; //varTypeMap["x4"] = VAR_x4;
184 //varTypeMap["x5"] = VAR_x5; varTypeMap["x6"] = VAR_x6;
185 //varTypeMap["x7"] = VAR_x7; varTypeMap["x8"] = VAR_x8;
186 //varTypeMap["x9"] = VAR_x9; varTypeMap["x10"] = VAR_x10; break;
187 //case SHORT_COLUMN: case LF_SHORT_COLUMN: case MF_SHORT_COLUMN:
188 varTypeMap["b"] = VAR_b; varTypeMap["h"] = VAR_h;
189 varTypeMap["P"] = VAR_P; varTypeMap["M"] = VAR_M; varTypeMap["Y"] = VAR_Y;
190 //break;
191 //case MF_ROSENBROCK: case MF_SHORT_COLUMN: (add to previous)
192 varTypeMap["ModelForm"] = VAR_MForm;
193 //case CANTILEVER_BEAM: case MOD_CANTILEVER_BEAM:
194 varTypeMap["w"] = VAR_w; varTypeMap["t"] = VAR_t; varTypeMap["R"] = VAR_R;
195 varTypeMap["E"] = VAR_E; varTypeMap["X"] = VAR_X; varTypeMap["area_type"] = VAR_area_type;
196 //varTypeMap["Y"] = VAR_Y; break;
197 //case STEEL_COLUMN:
198 varTypeMap["Fs"] = VAR_Fs; varTypeMap["P1"] = VAR_P1;
199 varTypeMap["P2"] = VAR_P2; varTypeMap["P3"] = VAR_P3;
200 varTypeMap["B"] = VAR_B; varTypeMap["D"] = VAR_D;
201 varTypeMap["H"] = VAR_H; //varTypeMap["b"] = VAR_b;
202 varTypeMap["d"] = VAR_d; //varTypeMap["h"] = VAR_h;
203 varTypeMap["F0"] = VAR_F0; //varTypeMap["E"] = VAR_E; break;
204 //case PROBLEM18:
205 varTypeMap["x"] = VAR_x; varTypeMap["xi"] = VAR_xi;
206 varTypeMap["Af"] = VAR_Af; varTypeMap["Ac"] = VAR_Ac;
207 //}
208 }
209 }
210
211
~TestDriverInterface()212 TestDriverInterface::~TestDriverInterface()
213 {
214 // No tear-down required for now
215 }
216
217
218 /** Derived map to evaluate a particular built-in test analysis function */
derived_map_ac(const String & ac_name)219 int TestDriverInterface::derived_map_ac(const String& ac_name)
220 {
221 // NOTE: a Factory pattern might be appropriate in the future to manage the
222 // conditional presence of linked subroutines in TestDriverInterface.
223
224 #ifdef MPI_DEBUG
225 Cout << "analysis server " << analysisServerId << " invoking " << ac_name
226 << " within TestDriverInterface." << std::endl;
227 #endif // MPI_DEBUG
228 int fail_code = 0;
229 std::map<String, driver_t>::iterator sd_iter = driverTypeMap.find(ac_name);
230 driver_t ac_type
231 = (sd_iter!=driverTypeMap.end()) ? sd_iter->second : NO_DRIVER;
232 switch (ac_type) {
233 case CANTILEVER_BEAM:
234 fail_code = cantilever(); break;
235 case MOD_CANTILEVER_BEAM:
236 fail_code = mod_cantilever(); break;
237 case CANTILEVER_BEAM_ML:
238 fail_code = cantilever_ml(); break;
239 case CYLINDER_HEAD:
240 fail_code = cyl_head(); break;
241 case ROSENBROCK:
242 fail_code = rosenbrock(); break;
243 case MODIFIED_ROSENBROCK:
244 fail_code = modified_rosenbrock(); break;
245 case GENERALIZED_ROSENBROCK:
246 fail_code = generalized_rosenbrock(); break;
247 case EXTENDED_ROSENBROCK:
248 fail_code = extended_rosenbrock(); break;
249 case LF_ROSENBROCK:
250 fail_code = lf_rosenbrock(); break;
251 case EXTRA_LF_ROSENBROCK:
252 fail_code = extra_lf_rosenbrock(); break;
253 case MF_ROSENBROCK:
254 fail_code = mf_rosenbrock(); break;
255 case LF_POLY_PROD:
256 fail_code = lf_poly_prod(); break;
257 case POLY_PROD:
258 fail_code = poly_prod(); break;
259 case GERSTNER:
260 fail_code = gerstner(); break;
261 case SCALABLE_GERSTNER:
262 fail_code = scalable_gerstner(); break;
263 case LOGNORMAL_RATIO:
264 fail_code = log_ratio(); break;
265 case MULTIMODAL:
266 fail_code = multimodal(); break;
267 case SHORT_COLUMN:
268 fail_code = short_column(); break;
269 case LF_SHORT_COLUMN:
270 fail_code = lf_short_column(); break;
271 case MF_SHORT_COLUMN:
272 fail_code = mf_short_column(); break;
273 case SIDE_IMPACT_COST:
274 fail_code = side_impact_cost(); break;
275 case SIDE_IMPACT_PERFORMANCE:
276 fail_code = side_impact_perf(); break;
277 case SOBOL_RATIONAL:
278 fail_code = sobol_rational(); break;
279 case SOBOL_G_FUNCTION:
280 fail_code = sobol_g_function(); break;
281 case SOBOL_ISHIGAMI:
282 fail_code = sobol_ishigami(); break;
283 case STEEL_COLUMN_COST:
284 fail_code = steel_column_cost(); break;
285 case STEEL_COLUMN_PERFORMANCE:
286 fail_code = steel_column_perf(); break;
287 case TEXT_BOOK:
288 fail_code = text_book(); break;
289 case TEXT_BOOK1: // for testing concurrent analyses
290 fail_code = text_book1(); break;
291 case TEXT_BOOK2: // for testing concurrent analyses
292 fail_code = text_book2(); break;
293 case TEXT_BOOK3: // for testing concurrent analyses
294 fail_code = text_book3(); break;
295 case TEXT_BOOK_OUU:
296 fail_code = text_book_ouu(); break;
297 case SCALABLE_TEXT_BOOK:
298 fail_code = scalable_text_book(); break;
299 case SCALABLE_MONOMIALS:
300 fail_code = scalable_monomials(); break;
301 case MOGATEST1:
302 fail_code = mogatest1(); break;
303 case MOGATEST2:
304 fail_code = mogatest2(); break;
305 case MOGATEST3:
306 fail_code = mogatest3(); break;
307 case ILLUMINATION:
308 fail_code = illumination(); break;
309 case BARNES:
310 fail_code = barnes(); break;
311 case BARNES_LF:
312 fail_code = barnes_lf(); break;
313 case HERBIE:
314 fail_code = herbie(); break;
315 case SMOOTH_HERBIE:
316 fail_code = smooth_herbie(); break;
317 case SHUBERT:
318 fail_code = shubert(); break;
319 #ifdef DAKOTA_SALINAS
320 case SALINAS:
321 fail_code = salinas(); break;
322 #endif // DAKOTA_SALINAS
323 #ifdef DAKOTA_MODELCENTER
324 case MODELCENTER: //case MC_API_RUN:
325 fail_code = mc_api_run(); break;
326 #endif // DAKOTA_MODELCENTER
327 case GENZ:
328 fail_code = genz(); break;
329 case DAMPED_OSCILLATOR:
330 fail_code = damped_oscillator(); break;
331 case STEADY_STATE_DIFFUSION_1D:
332 fail_code = steady_state_diffusion_1d(); break;
333 case SS_DIFFUSION_DISCREPANCY:
334 fail_code = ss_diffusion_discrepancy(); break;
335 case TRANSIENT_DIFFUSION_1D:
336 fail_code = transient_diffusion_1d(); break;
337 case PREDATOR_PREY:
338 fail_code = predator_prey(); break;
339 case ANISOTROPIC_QUADRATIC_FORM:
340 fail_code = aniso_quad_form(); break;
341 case BAYES_LINEAR:
342 fail_code = bayes_linear(); break;
343 case PROBLEM18:
344 fail_code = problem18(); break;
345 default: {
346 Cerr << "Error: analysis_driver '" << ac_name << "' is not available in "
347 << "the direct interface." << std::endl;
348 abort_handler(INTERFACE_ERROR);
349 }
350 }
351
352 // Failure capturing
353 if (fail_code) {
354 std::string err_msg("Error evaluating direct analysis_driver ");
355 err_msg += ac_name;
356 throw FunctionEvalFailure(err_msg);
357 }
358
359 return 0;
360 }
361
362
363 // -----------------------------------------
364 // Begin direct interfaces to test functions
365 // -----------------------------------------
cantilever()366 int TestDriverInterface::cantilever()
367 {
368 using std::pow;
369
370 if (multiProcAnalysisFlag) {
371 Cerr << "Error: cantilever direct fn does not support multiprocessor "
372 << "analyses." << std::endl;
373 abort_handler(-1);
374 }
375 // cantilever normally has 6 variables: 2 design + 4 uncertain
376 // If, however, design variables are _inserted_ into the uncertain variable
377 // distribution parameters (e.g., dakota_rbdo_cantilever_mapvars.in) instead
378 // of augmenting the uncertain variables, then the number of variables is 4.
379 // Design gradients are not supported for the case of design var insertion.
380 if ( (numVars != 4 && numVars != 6) || numADIV || numADRV ||//var count, no dv
381 (gradFlag && numVars == 4 && numDerivVars != 4) ) { // design insertion
382 Cerr << "Error: Bad number of variables in cantilever direct fn."
383 << std::endl;
384 abort_handler(INTERFACE_ERROR);
385 }
386 if (numFns < 2 || numFns > 3) {
387 Cerr << "Error: Bad number of functions in cantilever direct fn."
388 << std::endl;
389 abort_handler(INTERFACE_ERROR);
390 }
391
392 // Compute the cross-sectional area, stress, and displacement of the
393 // cantilever beam. This simulator is unusual in that it must support both
394 // the case of design variable insertion and the case of design variable
395 // augmentation. It does not support mixed insertion/augmentation. In
396 // the 6 variable case, w,t,R,E,X,Y are all passed in; in the 4 variable
397 // case, w,t assume local values.
398 std::map<var_t, Real>::iterator m_iter = xCM.find(VAR_w);
399 Real w = (m_iter == xCM.end()) ? 2.5 : m_iter->second; // beam width
400 m_iter = xCM.find(VAR_t);
401 Real t = (m_iter == xCM.end()) ? 2.5 : m_iter->second; // beam thickness
402 Real R = xCM[VAR_R], // yield strength
403 E = xCM[VAR_E], // Young's modulus
404 X = xCM[VAR_X], // horizontal load
405 Y = xCM[VAR_Y]; // vertical load
406
407 // allow f,c1,c2 (optimization) or just c1,c2 (calibration)
408 bool objective; size_t c1i, c2i;
409 if (numFns == 2) { objective = false; c1i = 0; c2i = 1; }
410 else { objective = true; c1i = 1; c2i = 2; }
411
412 // optimization inequality constraint: <= 0 and scaled O(1)
413 //Real g_stress = stress/R - 1.0;
414 //Real g_disp = disp/D0 - 1.0;
415
416 Real D0 = 2.2535, L = 100., area = w*t, w_sq = w*w, t_sq = t*t,
417 R_sq = R*R, X_sq = X*X, Y_sq = Y*Y;
418 Real stress = 600.*Y/w/t_sq + 600.*X/w_sq/t;
419 Real D1 = 4.*pow(L,3)/E/area, D2 = pow(Y/t_sq, 2)+pow(X/w_sq, 2),
420 D3 = D1/std::sqrt(D2)/D0, D4 = D1*std::sqrt(D2)/D0;
421
422 // **** f:
423 if (objective && (directFnASV[0] & 1))
424 fnVals[0] = area;
425
426 // **** c1:
427 if (directFnASV[c1i] & 1)
428 fnVals[c1i] = stress/R - 1.;
429
430 // **** c2:
431 if (directFnASV[c2i] & 1)
432 fnVals[c2i] = D4 - 1.;
433
434 // **** df/dx:
435 if (objective && (directFnASV[0] & 2))
436 for (size_t i=0; i<numDerivVars; ++i)
437 switch (varTypeDVV[i]) {
438 case VAR_w: fnGrads[0][i] = t; break; // design var derivative
439 case VAR_t: fnGrads[0][i] = w; break; // design var derivative
440 default: fnGrads[0][i] = 0.; break; // uncertain var derivative
441 }
442
443 // **** dc1/dx:
444 if (directFnASV[c1i] & 2)
445 for (size_t i=0; i<numDerivVars; ++i)
446 switch (varTypeDVV[i]) {
447 case VAR_w: fnGrads[c1i][i] = -600.*(Y/t + 2.*X/w)/w_sq/t/R; break;// dv
448 case VAR_t: fnGrads[c1i][i] = -600.*(2.*Y/t + X/w)/w/t_sq/R; break;// dv
449 case VAR_R: fnGrads[c1i][i] = -stress/R_sq; break; // uncertain var deriv
450 case VAR_E: fnGrads[c1i][i] = 0.; break; // uncertain var deriv
451 case VAR_X: fnGrads[c1i][i] = 600./w_sq/t/R; break; // uncertain var deriv
452 case VAR_Y: fnGrads[c1i][i] = 600./w/t_sq/R; break; // uncertain var deriv
453 }
454
455 // **** dc2/dx:
456 if (directFnASV[c2i] & 2)
457 for (size_t i=0; i<numDerivVars; ++i)
458 switch (varTypeDVV[i]) {
459 case VAR_w: fnGrads[c2i][i] = -D3*2.*X_sq/w_sq/w_sq/w - D4/w; break;// dv
460 case VAR_t: fnGrads[c2i][i] = -D3*2.*Y_sq/t_sq/t_sq/t - D4/t; break;// dv
461 case VAR_R: fnGrads[c2i][i] = 0.; break; // unc var deriv
462 case VAR_E: fnGrads[c2i][i] = -D4/E; break; // unc var deriv
463 case VAR_X: fnGrads[c2i][i] = D3*X/w_sq/w_sq; break; // unc var deriv
464 case VAR_Y: fnGrads[c2i][i] = D3*Y/t_sq/t_sq; break; // unc var deriv
465 }
466
467 // **** d^2f/dx^2:
468 if (objective && (directFnASV[0] & 4))
469 for (size_t i=0; i<numDerivVars; ++i)
470 for (size_t j=0; j<=i; ++j)
471 fnHessians[0](i,j)
472 = ( (varTypeDVV[i] == VAR_w && varTypeDVV[j] == VAR_t) ||
473 (varTypeDVV[i] == VAR_t && varTypeDVV[j] == VAR_w) ) ? 1. : 0.;
474
475 // **** d^2c1/dx^2:
476 if (directFnASV[c1i] & 4) {
477 for (size_t i=0; i<numDerivVars; ++i)
478 for (size_t j=0; j<=i; ++j)
479 if (varTypeDVV[i] == VAR_w && varTypeDVV[j] == VAR_w)
480 fnHessians[c1i](i,j) = 1200.*(Y/t + 3.*X/w)/w_sq/area/R;
481 else if (varTypeDVV[i] == VAR_t && varTypeDVV[j] == VAR_t)
482 fnHessians[c1i](i,j) = 1200.*(3.*Y/t + X/w)/t_sq/area/R;
483 else if (varTypeDVV[i] == VAR_R && varTypeDVV[j] == VAR_R)
484 fnHessians[c1i](i,j) = 2.*stress/pow(R, 3);
485 else if ( (varTypeDVV[i] == VAR_w && varTypeDVV[j] == VAR_t) ||
486 (varTypeDVV[i] == VAR_t && varTypeDVV[j] == VAR_w) )
487 fnHessians[c1i](i,j) = 1200.*(Y/t + X/w)/w_sq/t_sq/R;
488 else if ( (varTypeDVV[i] == VAR_w && varTypeDVV[j] == VAR_R) ||
489 (varTypeDVV[i] == VAR_R && varTypeDVV[j] == VAR_w) )
490 fnHessians[c1i](i,j) = 600.*(Y/t + 2.*X/w)/w_sq/t/R_sq;
491 else if ( (varTypeDVV[i] == VAR_w && varTypeDVV[j] == VAR_X) ||
492 (varTypeDVV[i] == VAR_X && varTypeDVV[j] == VAR_w) )
493 fnHessians[c1i](i,j) = -1200./w_sq/w/t/R;
494 else if ( (varTypeDVV[i] == VAR_w && varTypeDVV[j] == VAR_Y) ||
495 (varTypeDVV[i] == VAR_Y && varTypeDVV[j] == VAR_w) )
496 fnHessians[c1i](i,j) = -600./w_sq/t_sq/R;
497 else if ( (varTypeDVV[i] == VAR_t && varTypeDVV[j] == VAR_R) ||
498 (varTypeDVV[i] == VAR_R && varTypeDVV[j] == VAR_t) )
499 fnHessians[c1i](i,j) = 600.*(2.*Y/t + X/w)/w/t_sq/R_sq;
500 else if ( (varTypeDVV[i] == VAR_t && varTypeDVV[j] == VAR_X) ||
501 (varTypeDVV[i] == VAR_X && varTypeDVV[j] == VAR_t) )
502 fnHessians[c1i](i,j) = -600./w_sq/t_sq/R;
503 else if ( (varTypeDVV[i] == VAR_t && varTypeDVV[j] == VAR_Y) ||
504 (varTypeDVV[i] == VAR_Y && varTypeDVV[j] == VAR_t) )
505 fnHessians[c1i](i,j) = -1200./w/t_sq/t/R;
506 else if ( (varTypeDVV[i] == VAR_R && varTypeDVV[j] == VAR_X) ||
507 (varTypeDVV[i] == VAR_X && varTypeDVV[j] == VAR_R) )
508 fnHessians[c1i](i,j) = -600./w_sq/t/R_sq;
509 else if ( (varTypeDVV[i] == VAR_R && varTypeDVV[j] == VAR_Y) ||
510 (varTypeDVV[i] == VAR_Y && varTypeDVV[j] == VAR_R) )
511 fnHessians[c1i](i,j) = -600./w/t_sq/R_sq;
512 else
513 fnHessians[c1i](i,j) = 0.;
514 }
515
516 // **** d^2c2/dx^2:
517 if (directFnASV[c2i] & 4) {
518 Real D5 = 1./std::sqrt(D2)/D0, D6 = -D1/2./D0/pow(D2,1.5);
519 Real D7 = std::sqrt(D2)/D0, D8 = D1/2./D0/std::sqrt(D2);
520 Real dD2_dX = 2.*X/w_sq/w_sq, dD3_dX = D6*dD2_dX, dD4_dX = D8*dD2_dX;
521 Real dD2_dY = 2.*Y/t_sq/t_sq, dD3_dY = D6*dD2_dY, dD4_dY = D8*dD2_dY;
522 Real dD1_dw = -D1/w, dD2_dw = -4.*X_sq/w_sq/w_sq/w,
523 dD3_dw = D5*dD1_dw + D6*dD2_dw, dD4_dw = D7*dD1_dw + D8*dD2_dw;
524 Real dD1_dt = -D1/t, dD2_dt = -4.*Y_sq/t_sq/t_sq/t,
525 dD3_dt = D5*dD1_dt + D6*dD2_dt, dD4_dt = D7*dD1_dt + D8*dD2_dt;
526 for (size_t i=0; i<numDerivVars; ++i)
527 for (size_t j=0; j<=i; ++j)
528 if (varTypeDVV[i] == VAR_w && varTypeDVV[j] == VAR_w)
529 fnHessians[c2i](i,j) = D3*10.*X_sq/pow(w_sq,3)
530 - 2.*X_sq/w_sq/w_sq/w*dD3_dw + D4/w_sq - dD4_dw/w;
531 else if (varTypeDVV[i] == VAR_t && varTypeDVV[j] == VAR_t)
532 fnHessians[c2i](i,j) = D3*10.*Y_sq/pow(t_sq,3)
533 - 2.*Y_sq/t_sq/t_sq/t*dD3_dt + D4/t_sq - dD4_dt/t;
534 else if (varTypeDVV[i] == VAR_E && varTypeDVV[j] == VAR_E) {
535 Real dD1_dE = -D1/E, dD4_dE = D7*dD1_dE;
536 fnHessians[c2i](i,j) = D4/E/E - dD4_dE/E;
537 }
538 else if (varTypeDVV[i] == VAR_X && varTypeDVV[j] == VAR_X)
539 fnHessians[c2i](i,j) = D3/w_sq/w_sq + X/w_sq/w_sq*dD3_dX;
540 else if (varTypeDVV[i] == VAR_Y && varTypeDVV[j] == VAR_Y)
541 fnHessians[c2i](i,j) = D3/t_sq/t_sq + Y/t_sq/t_sq*dD3_dY;
542 else if ( (varTypeDVV[i] == VAR_w && varTypeDVV[j] == VAR_t) ||
543 (varTypeDVV[i] == VAR_t && varTypeDVV[j] == VAR_w) )
544 fnHessians[c2i](i,j) = -2.*X_sq/w_sq/w_sq/w*dD3_dt - dD4_dt/w;
545 else if ( (varTypeDVV[i] == VAR_w && varTypeDVV[j] == VAR_E) ||
546 (varTypeDVV[i] == VAR_E && varTypeDVV[j] == VAR_w) )
547 fnHessians[c2i](i,j) = -dD4_dw/E;
548 else if ( (varTypeDVV[i] == VAR_w && varTypeDVV[j] == VAR_X) ||
549 (varTypeDVV[i] == VAR_X && varTypeDVV[j] == VAR_w) )
550 fnHessians[c2i](i,j) = -4.*X*D3/w_sq/w_sq/w + X/w_sq/w_sq*dD3_dw;
551 else if ( (varTypeDVV[i] == VAR_w && varTypeDVV[j] == VAR_Y) ||
552 (varTypeDVV[i] == VAR_Y && varTypeDVV[j] == VAR_w) )
553 fnHessians[c2i](i,j) = Y/t_sq/t_sq*dD3_dw;
554 else if ( (varTypeDVV[i] == VAR_t && varTypeDVV[j] == VAR_E) ||
555 (varTypeDVV[i] == VAR_E && varTypeDVV[j] == VAR_t) )
556 fnHessians[c2i](i,j) = -dD4_dt/E;
557 else if ( (varTypeDVV[i] == VAR_t && varTypeDVV[j] == VAR_X) ||
558 (varTypeDVV[i] == VAR_X && varTypeDVV[j] == VAR_t) )
559 fnHessians[c2i](i,j) = X/w_sq/w_sq*dD3_dt;
560 else if ( (varTypeDVV[i] == VAR_t && varTypeDVV[j] == VAR_Y) ||
561 (varTypeDVV[i] == VAR_Y && varTypeDVV[j] == VAR_t) )
562 fnHessians[c2i](i,j) = -4.*Y*D3/t_sq/t_sq/t + Y/t_sq/t_sq*dD3_dt;
563 else if ( (varTypeDVV[i] == VAR_E && varTypeDVV[j] == VAR_X) ||
564 (varTypeDVV[i] == VAR_X && varTypeDVV[j] == VAR_E) )
565 fnHessians[c2i](i,j) = -dD4_dX/E;
566 else if ( (varTypeDVV[i] == VAR_E && varTypeDVV[j] == VAR_Y) ||
567 (varTypeDVV[i] == VAR_Y && varTypeDVV[j] == VAR_E) )
568 fnHessians[c2i](i,j) = -dD4_dY/E;
569 else if ( (varTypeDVV[i] == VAR_X && varTypeDVV[j] == VAR_Y) ||
570 (varTypeDVV[i] == VAR_Y && varTypeDVV[j] == VAR_X) )
571 fnHessians[c2i](i,j) = X/w_sq/w_sq*dD3_dY;
572 else
573 fnHessians[c2i](i,j) = 0.;
574 }
575
576 return 0; // no failure
577 }
578
579
mod_cantilever()580 int TestDriverInterface::mod_cantilever()
581 {
582 using std::pow;
583
584 if (multiProcAnalysisFlag) {
585 Cerr << "Error: cantilever direct fn does not support multiprocessor "
586 << "analyses." << std::endl;
587 abort_handler(-1);
588 }
589 // cantilever normally has 6 variables: 2 design + 4 uncertain
590 // If, however, design variables are _inserted_ into the uncertain variable
591 // distribution parameters (e.g., dakota_rbdo_cantilever_mapvars.in) instead
592 // of augmenting the uncertain variables, then the number of variables is 4.
593 // Design gradients are not supported for the case of design var insertion.
594 if ( (numVars != 4 && numVars != 6) || numADIV || numADRV ||//var count, no dv
595 (gradFlag && numVars == 4 && numDerivVars != 4) ) { // design insertion
596 Cerr << "Error: Bad number of variables in cantilever direct fn."
597 << std::endl;
598 abort_handler(INTERFACE_ERROR);
599 }
600 if (numFns < 2 || numFns > 3) {
601 Cerr << "Error: Bad number of functions in mod_cantilever direct fn."
602 << std::endl;
603 abort_handler(INTERFACE_ERROR);
604 }
605
606 // Compute the cross-sectional area, stress, and displacement of the
607 // cantilever beam. This simulator is unusual in that it must support both
608 // the case of design variable insertion and the case of design variable
609 // augmentation. It does not support mixed insertion/augmentation. In
610 // the 6 variable case, w,t,R,E,X,Y are all passed in; in the 4 variable
611 // case, w,t assume local values.
612 std::map<var_t, Real>::iterator m_iter = xCM.find(VAR_w);
613 Real w = (m_iter == xCM.end()) ? 2.5 : m_iter->second; // beam width
614 m_iter = xCM.find(VAR_t);
615 Real t = (m_iter == xCM.end()) ? 2.5 : m_iter->second; // beam thickness
616 Real R = xCM[VAR_R], // yield strength
617 E = xCM[VAR_E], // Young's modulus
618 X = xCM[VAR_X], // horizontal load
619 Y = xCM[VAR_Y]; // vertical load
620
621 // allow f,c1,c2 (optimization) or just c1,c2 (calibration)
622 bool objective; size_t c1i, c2i;
623 if (numFns == 2) { objective = false; c1i = 0; c2i = 1; }
624 else { objective = true; c1i = 1; c2i = 2; }
625
626 // UQ limit state <= 0: don't scale stress by random variable r
627 //double g_stress = stress - r;
628 //double g_disp = displ - D0;
629
630 Real D0 = 2.2535, L = 100., area = w*t, w_sq = w*w, t_sq = t*t,
631 R_sq = R*R, X_sq = X*X, Y_sq = Y*Y;
632 Real stress = 600.*Y/w/t_sq + 600.*X/w_sq/t;
633 Real D1 = 4.*pow(L,3)/E/area, D2 = pow(Y/t_sq, 2)+pow(X/w_sq, 2),
634 D3 = D1/std::sqrt(D2), displ = D1*std::sqrt(D2);
635
636 // **** f:
637 if (objective && (directFnASV[0] & 1))
638 fnVals[0] = area;
639
640 // **** c1:
641 if (directFnASV[c1i] & 1)
642 fnVals[c1i] = stress - R;
643
644 // **** c2:
645 if (directFnASV[c2i] & 1)
646 fnVals[c2i] = displ - D0;
647
648 // **** df/dx:
649 if (objective && (directFnASV[0] & 2))
650 for (size_t i=0; i<numDerivVars; ++i)
651 switch (varTypeDVV[i]) {
652 case VAR_w: fnGrads[0][i] = t; break; // design var derivative
653 case VAR_t: fnGrads[0][i] = w; break; // design var derivative
654 default: fnGrads[0][i] = 0.; break; // uncertain var derivative
655 }
656
657 // **** dc1/dx:
658 if (directFnASV[c1i] & 2)
659 for (size_t i=0; i<numDerivVars; ++i)
660 switch (varTypeDVV[i]) {
661 case VAR_w: fnGrads[c1i][i] = -600.*(Y/t + 2.*X/w)/w_sq/t; break;//des var
662 case VAR_t: fnGrads[c1i][i] = -600.*(2.*Y/t + X/w)/w/t_sq; break;//des var
663 case VAR_R: fnGrads[c1i][i] = -1.; break; // uncertain var deriv
664 case VAR_E: fnGrads[c1i][i] = 0.; break; // uncertain var deriv
665 case VAR_X: fnGrads[c1i][i] = 600./w_sq/t; break; // uncertain var deriv
666 case VAR_Y: fnGrads[c1i][i] = 600./w/t_sq; break; // uncertain var deriv
667 }
668
669 // **** dc2/dx:
670 if (directFnASV[c2i] & 2)
671 for (size_t i=0; i<numDerivVars; ++i)
672 switch (varTypeDVV[i]) {
673 case VAR_w: fnGrads[c2i][i] = -D3*2.*X_sq/w_sq/w_sq/w - displ/w; break;
674 case VAR_t: fnGrads[c2i][i] = -D3*2.*Y_sq/t_sq/t_sq/t - displ/t; break;
675 case VAR_R: fnGrads[c2i][i] = 0.; break; // unc var deriv
676 case VAR_E: fnGrads[c2i][i] = -displ/E; break; // unc var deriv
677 case VAR_X: fnGrads[c2i][i] = D3*X/w_sq/w_sq; break; // unc var deriv
678 case VAR_Y: fnGrads[c2i][i] = D3*Y/t_sq/t_sq; break; // unc var deriv
679 }
680
681 /* Alternative modification: take E out of displ denominator to remove
682 singularity in tail (at 20 std deviations). In PCE/SC testing, this
683 had minimal impact and did not justify the nonstandard form.
684
685 Real D0 = 2.2535, L = 100., area = w*t, w_sq = w*w, t_sq = t*t,
686 R_sq = R*R, X_sq = X*X, Y_sq = Y*Y;
687 Real stress = 600.*Y/w/t_sq + 600.*X/w_sq/t;
688 Real D1 = 4.*pow(L,3)/area, D2 = pow(Y/t_sq, 2)+pow(X/w_sq, 2),
689 D3 = D1/std::sqrt(D2), D4 = D1*std::sqrt(D2);
690
691 // **** c2:
692 if (directFnASV[c2i] & 1)
693 fnVals[c2i] = D4 - D0*E;
694
695 // **** dc2/dx:
696 if (directFnASV[c2i] & 2)
697 for (size_t i=0; i<numDerivVars; ++i)
698 switch (varTypeDVV[i]) {
699 case VAR_w: fnGrads[c2i][i] = -D3*2.*X_sq/w_sq/w_sq/w - D4/w; break;// des
700 case VAR_t: fnGrads[c2i][i] = -D3*2.*Y_sq/t_sq/t_sq/t - D4/t; break;// des
701 case VAR_R: fnGrads[c2i][i] = 0.; break; // unc var deriv
702 case VAR_E: fnGrads[c2i][i] = -D0; break; // unc var deriv
703 case VAR_X: fnGrads[c2i][i] = D3*X/w_sq/w_sq; break; // unc var deriv
704 case VAR_Y: fnGrads[c2i][i] = D3*Y/t_sq/t_sq; break; // unc var deriv
705 }
706 */
707
708 return 0; // no failure
709 }
710
cantilever_ml()711 int TestDriverInterface::cantilever_ml()
712 {
713 using std::pow;
714
715 if (multiProcAnalysisFlag) {
716 Cerr << "Error: cantilever direct fn does not support multiprocessor "
717 << "analyses." << std::endl;
718 abort_handler(-1);
719 }
720 // cantilever normally has 6 variables: 2 design + 4 uncertain
721 // If, however, design variables are _inserted_ into the uncertain variable
722 // distribution parameters (e.g., dakota_rbdo_cantilever_mapvars.in) instead
723 // of augmenting the uncertain variables, then the number of variables is 4.
724 // Design gradients are not supported for the case of design var insertion.
725 if ( (numVars != 5 && numVars != 7) || (numADIV != 1) || numADRV ||//var count, no dv
726 (gradFlag && numVars == 5 && numDerivVars != 4) ) { // design insertion
727 Cerr << "Error: Bad number of variables in cantilever direct fn."
728 << std::endl;
729 Cerr << "Num vars:" << numVars << ", " << numDerivVars
730 << std::endl;
731 abort_handler(INTERFACE_ERROR);
732 }
733 if (numFns < 2 || numFns > 3) {
734 Cerr << "Error: Bad number of functions in mod_cantilever direct fn."
735 << std::endl;
736 abort_handler(INTERFACE_ERROR);
737 }
738
739 // Compute the cross-sectional area, stress, and displacement of the
740 // cantilever beam. This simulator is unusual in that it must support both
741 // the case of design variable insertion and the case of design variable
742 // augmentation. It does not support mixed insertion/augmentation. In
743 // the 6 variable case, w,t,R,E,X,Y are all passed in; in the 4 variable
744 // case, w,t assume local values.
745 std::map<var_t, Real>::iterator m_iter = xCM.find(VAR_w);
746 Real w = (m_iter == xCM.end()) ? 2.5 : m_iter->second; // beam width
747 m_iter = xCM.find(VAR_t);
748 Real t = (m_iter == xCM.end()) ? 2.5 : m_iter->second; // beam thickness
749 Real R = xCM[VAR_R], // yield strength
750 E = xCM[VAR_E], // Young's modulus
751 X = xCM[VAR_X], // horizontal load
752 Y = xCM[VAR_Y]; // vertical load
753
754 // allow f,c1,c2 (optimization) or just c1,c2 (calibration)
755 bool objective; size_t c1i, c2i;
756 if (numFns == 2) { objective = false; c1i = 0; c2i = 1; }
757 else { objective = true; c1i = 1; c2i = 2; }
758
759 // UQ limit state <= 0: don't scale stress by random variable r
760 //double g_stress = stress - r;
761 //double g_disp = displ - D0;
762
763 std::map<var_t, int>::iterator area_type_iter = xDIM.find(VAR_area_type);
764 int area_type = (area_type_iter == xDIM.end()) ? 1. : area_type_iter->second; // Correlation Af for objective
765
766 Real D0 = 2.2535, L = 100., area, w_sq, t_sq, R_sq, X_sq, Y_sq;
767 Real stress;
768 Real D1, D2, D3, displ;
769 area = w*t;
770 if(area_type == 1){// Rectangle
771 w_sq = w*w; t_sq = t*t;
772 R_sq = R*R; X_sq = X*X; Y_sq = Y*Y;
773
774 stress = 6.*L*Y/w/t_sq + 6.*L*X/w_sq/t;
775
776 D1 = 4.*pow(L, 3)/E/area;
777 D2 = pow(Y/t_sq, 2)+pow(X/w_sq, 2);
778 D3 = D1/std::sqrt(D2);
779 displ = D1*std::sqrt(D2);
780 }else if(area_type == 2){// Ellipse
781 const Real m_pi = 3.14159265358979323846;
782 Real a = t/2. * 4./m_pi;
783 Real b = w/2.;
784
785 stress = (4.*L)/(m_pi*a*b) * std::sqrt(pow(Y/a, 2) + pow(X/b, 2));
786
787 Real I_x = (m_pi * b * pow(a, 3))/4.;
788 Real I_y = (m_pi * pow(b, 3) * a)/4.;
789 displ = std::sqrt(
790 pow((pow(L, 3) * X)/(3.*E*I_y), 2) +
791 pow((pow(L, 3) * Y)/(3.*E*I_x), 2)
792 );
793 }else{
794 Cout << "TestDriverInterface::mod_cantilever_ml(): wrong area type.\n";
795 abort_handler(INTERFACE_ERROR);
796 }
797 // **** f:
798 if (objective && (directFnASV[0] & 1))
799 fnVals[0] = area;
800
801 // **** c1:
802 if (directFnASV[c1i] & 1)
803 fnVals[c1i] = stress/R - 1.;
804
805 // **** c2:
806 if (directFnASV[c2i] & 1)
807 fnVals[c2i] = displ/D0 - 1.;
808
809 // **** df/dx:
810 if (objective && (directFnASV[0] & 2))
811 for (size_t i=0; i<numDerivVars && area_type == 1; ++i)
812 switch (varTypeDVV[i]) {
813 case VAR_w: fnGrads[0][i] = t; break; // design var derivative
814 case VAR_t: fnGrads[0][i] = w; break; // design var derivative
815 default: fnGrads[0][i] = 0.; break; // uncertain var derivative
816 }
817
818 // **** dc1/dx:
819 if (directFnASV[c1i] & 2)
820 for (size_t i=0; i<numDerivVars && area_type == 1; ++i)
821 switch (varTypeDVV[i]) {
822 case VAR_w: fnGrads[c1i][i] = -600.*(Y/t + 2.*X/w)/w_sq/t; break;//des var
823 case VAR_t: fnGrads[c1i][i] = -600.*(2.*Y/t + X/w)/w/t_sq; break;//des var
824 case VAR_R: fnGrads[c1i][i] = -1.; break; // uncertain var deriv
825 case VAR_E: fnGrads[c1i][i] = 0.; break; // uncertain var deriv
826 case VAR_X: fnGrads[c1i][i] = 600./w_sq/t; break; // uncertain var deriv
827 case VAR_Y: fnGrads[c1i][i] = 600./w/t_sq; break; // uncertain var deriv
828 }
829
830 // **** dc2/dx:
831 if (directFnASV[c2i] & 2)
832 for (size_t i=0; i<numDerivVars && area_type == 1; ++i)
833 switch (varTypeDVV[i]) {
834 case VAR_w: fnGrads[c2i][i] = -D3*2.*X_sq/w_sq/w_sq/w - displ/w; break;
835 case VAR_t: fnGrads[c2i][i] = -D3*2.*Y_sq/t_sq/t_sq/t - displ/t; break;
836 case VAR_R: fnGrads[c2i][i] = 0.; break; // unc var deriv
837 case VAR_E: fnGrads[c2i][i] = -displ/E; break; // unc var deriv
838 case VAR_X: fnGrads[c2i][i] = D3*X/w_sq/w_sq; break; // unc var deriv
839 case VAR_Y: fnGrads[c2i][i] = D3*Y/t_sq/t_sq; break; // unc var deriv
840 }
841
842 /* Alternative modification: take E out of displ denominator to remove
843 singularity in tail (at 20 std deviations). In PCE/SC testing, this
844 had minimal impact and did not justify the nonstandard form.
845
846 Real D0 = 2.2535, L = 100., area = w*t, w_sq = w*w, t_sq = t*t,
847 R_sq = R*R, X_sq = X*X, Y_sq = Y*Y;
848 Real stress = 600.*Y/w/t_sq + 600.*X/w_sq/t;
849 Real D1 = 4.*pow(L,3)/area, D2 = pow(Y/t_sq, 2)+pow(X/w_sq, 2),
850 D3 = D1/std::sqrt(D2), D4 = D1*std::sqrt(D2);
851
852 // **** c2:
853 if (directFnASV[c2i] & 1)
854 fnVals[c2i] = D4 - D0*E;
855
856 // **** dc2/dx:
857 if (directFnASV[c2i] & 2)
858 for (size_t i=0; i<numDerivVars; ++i)
859 switch (varTypeDVV[i]) {
860 case VAR_w: fnGrads[c2i][i] = -D3*2.*X_sq/w_sq/w_sq/w - D4/w; break;// des
861 case VAR_t: fnGrads[c2i][i] = -D3*2.*Y_sq/t_sq/t_sq/t - D4/t; break;// des
862 case VAR_R: fnGrads[c2i][i] = 0.; break; // unc var deriv
863 case VAR_E: fnGrads[c2i][i] = -D0; break; // unc var deriv
864 case VAR_X: fnGrads[c2i][i] = D3*X/w_sq/w_sq; break; // unc var deriv
865 case VAR_Y: fnGrads[c2i][i] = D3*Y/t_sq/t_sq; break; // unc var deriv
866 }
867 */
868
869 return 0; // no failure
870 }
871
872
873
cyl_head()874 int TestDriverInterface::cyl_head()
875 {
876 if (multiProcAnalysisFlag) {
877 Cerr << "Error: cyl_head direct fn does not yet support "
878 << "multiprocessor analyses." << std::endl;
879 abort_handler(-1);
880 }
881 if (numVars != 2 || numADIV || numADRV || (gradFlag && numDerivVars != 2)) {
882 Cerr << "Error: Bad number of variables in cyl_head direct fn." <<std::endl;
883 abort_handler(INTERFACE_ERROR);
884 }
885 if (numFns != 4) {
886 Cerr << "Error: Bad number of functions in cyl_head direct fn." <<std::endl;
887 abort_handler(INTERFACE_ERROR);
888 }
889 if (hessFlag) {
890 Cerr << "Error: Hessians not supported in cyl_head direct fn." << std::endl;
891 abort_handler(INTERFACE_ERROR);
892 }
893
894 Real exhaust_offset = 1.34;
895 Real exhaust_dia = 1.556;
896 Real intake_offset = 3.25;
897 // Use nondimensional xC[1]:
898 // (0. <= nondimensional <= 4.), (0. in <= dimensional <= 0.004 in)
899 Real warranty = 100000. + 15000. * (4. - xC[1]);
900 Real cycle_time = 45. + 4.5*std::pow(4. - xC[1], 1.5);
901 Real wall_thickness = intake_offset - exhaust_offset - (xC[0]+exhaust_dia)/2.;
902 Real horse_power = 250.+200.*(xC[0]/1.833-1.);
903 Real max_stress = 750. + std::pow(std::fabs(wall_thickness),-2.5);
904
905 // **** f:
906 if (directFnASV[0] & 1)
907 fnVals[0] = -1.*(horse_power/250.+warranty/100000.);
908
909 // **** c1:
910 if (directFnASV[1] & 1)
911 fnVals[1] = max_stress/1500.-1.;
912
913 // **** c2:
914 if (directFnASV[2] & 1)
915 fnVals[2] = 1.-warranty/100000.;
916
917 // **** c3:
918 if (directFnASV[3] & 1)
919 fnVals[3] = cycle_time/60. - 1.;
920
921 // **** c4: (Unnecessary if intake_dia upper bound reduced to 2.164)
922 //if (directFnASV[4] & 1)
923 // fnVals[4] = 1.-20.*wall_thickness;
924
925 // **** df/dx:
926 if (directFnASV[0] & 2) {
927 fnGrads[0][0] = -.8/1.833;
928 fnGrads[0][1] = 0.15;
929 }
930
931 // **** dc1/dx:
932 if (directFnASV[1] & 2) {
933 fnGrads[1][0] = 1.25/1500*std::pow(wall_thickness, -3.5);
934 fnGrads[1][1] = 0.;
935 }
936
937 // **** dc2/dx:
938 if (directFnASV[2] & 2) {
939 fnGrads[2][0] = 0.;
940 fnGrads[2][1] = 0.15;
941 }
942
943 // **** dc3/dx:
944 if (directFnASV[3] & 2) {
945 fnGrads[3][0] = 0.;
946 fnGrads[3][1] = -0.1125*std::sqrt(4. - xC[1]);
947 }
948
949 return 0; // no failure
950 }
951
952
multimodal()953 int TestDriverInterface::multimodal()
954 {
955 if (multiProcAnalysisFlag) {
956 Cerr << "Error: multimodal direct fn does not support multiprocessor "
957 << "analyses." << std::endl;
958 abort_handler(-1);
959 }
960 if ( numVars != 2 || numADIV || numADRV ||
961 ( ( gradFlag || hessFlag ) && numDerivVars != 2 ) ) {
962 Cerr << "Error: Bad number of variables in multimodal direct fn."
963 << std::endl;
964 abort_handler(INTERFACE_ERROR);
965 }
966 if (numFns != 1) {
967 Cerr << "Error: Bad number of functions in multimodal direct fn."
968 << std::endl;
969 abort_handler(INTERFACE_ERROR);
970 }
971
972 // **** f:
973 if (directFnASV[0] & 1)
974 fnVals[0] = (xC[0]*xC[0]+4)*(xC[1]-1)/20 - std::sin(5*xC[0]/2) - 2;
975
976 // **** df/dx:
977 if (directFnASV[0] & 2) {
978 fnGrads[0][0] = xC[0]*(xC[1]-1)/10 - (5/2)*std::cos(5*xC[0]/2);
979 fnGrads[0][1] = (xC[0]*xC[0]+4)/20;
980 }
981
982 // **** d^2f/dx^2:
983 if (directFnASV[0] & 4) {
984 fnHessians[0](0,0) = (xC[1]-1)/10 + (25/4)*std::sin(5*xC[0]/2);
985 fnHessians[0](0,1) = fnHessians[0](1,0) = xC[0]/10;
986 fnHessians[0](1,1) = 0.0;
987 }
988
989 return 0; // no failure
990 }
991
rosenbrock()992 int TestDriverInterface::rosenbrock()
993 {
994 if (multiProcAnalysisFlag) {
995 Cerr << "Error: rosenbrock direct fn does not yet support multiprocessor "
996 << "analyses." << std::endl;
997 abort_handler(-1);
998 }
999 if (numACV != 2 || numADIV > 1 || numADRV) { // allow ModelForm discrete int
1000 Cerr << "Error: Bad number of variables in rosenbrock direct fn."
1001 << std::endl;
1002 abort_handler(INTERFACE_ERROR);
1003 }
1004 if (numFns > 2) { // 1 fn -> opt, 2 fns -> least sq
1005 Cerr << "Error: Bad number of functions in rosenbrock direct fn."
1006 << std::endl;
1007 abort_handler(INTERFACE_ERROR);
1008 }
1009
1010 bool least_sq_flag = (numFns > 1);
1011 Real x1 = xCM[VAR_x1], x2 = xCM[VAR_x2], f1 = x2-x1*x1, f2 = 1.-x1;
1012
1013 if (least_sq_flag) {
1014 // **** Residual R1:
1015 if (directFnASV[0] & 1)
1016 fnVals[0] = 10.*f1;
1017 // **** Residual R2:
1018 if (directFnASV[1] & 1)
1019 fnVals[1] = f2;
1020
1021 // **** dR1/dx:
1022 if (directFnASV[0] & 2)
1023 for (size_t i=0; i<numDerivVars; ++i)
1024 switch (varTypeDVV[i]) {
1025 case VAR_x1: fnGrads[0][i] = -20.*x1; break;
1026 case VAR_x2: fnGrads[0][i] = 10.; break;
1027 }
1028 // **** dR2/dx:
1029 if (directFnASV[1] & 2)
1030 for (size_t i=0; i<numDerivVars; ++i)
1031 switch (varTypeDVV[i]) {
1032 case VAR_x1: fnGrads[1][i] = -1.; break;
1033 case VAR_x2: fnGrads[1][i] = 0.; break;
1034 }
1035
1036 // **** d^2R1/dx^2:
1037 if (directFnASV[0] & 4)
1038 for (size_t i=0; i<numDerivVars; ++i)
1039 for (size_t j=0; j<=i; ++j)
1040 if (varTypeDVV[i] == VAR_x1 && varTypeDVV[j] == VAR_x1)
1041 fnHessians[0](i,j) = -20.;
1042 else
1043 fnHessians[0](i,j) = 0.;
1044 // **** d^2R2/dx^2:
1045 if (directFnASV[1] & 4)
1046 fnHessians[1] = 0.;
1047 }
1048 else {
1049 // **** f:
1050 if (directFnASV[0] & 1)
1051 fnVals[0] = 100.*f1*f1+f2*f2;
1052
1053 // **** df/dx:
1054 if (directFnASV[0] & 2)
1055 for (size_t i=0; i<numDerivVars; ++i)
1056 switch (varTypeDVV[i]) {
1057 case VAR_x1: fnGrads[0][i] = -400.*f1*x1 - 2.*f2; break;
1058 case VAR_x2: fnGrads[0][i] = 200.*f1; break;
1059 }
1060
1061 // **** d^2f/dx^2:
1062 if (directFnASV[0] & 4)
1063 for (size_t i=0; i<numDerivVars; ++i)
1064 for (size_t j=0; j<=i; ++j)
1065 if (varTypeDVV[i] == VAR_x1 && varTypeDVV[j] == VAR_x1)
1066 fnHessians[0](i,j) = -400.*(x2 - 3.*x1*x1) + 2.;
1067 else if ( (varTypeDVV[i] == VAR_x1 && varTypeDVV[j] == VAR_x2) ||
1068 (varTypeDVV[i] == VAR_x2 && varTypeDVV[j] == VAR_x1) )
1069 fnHessians[0](i,j) = -400.*x1;
1070 else if (varTypeDVV[i] == VAR_x2 && varTypeDVV[j] == VAR_x2)
1071 fnHessians[0](i,j) = 200.;
1072 }
1073
1074 return 0; // no failure
1075 }
1076
modified_rosenbrock()1077 int TestDriverInterface::modified_rosenbrock()
1078 {
1079 if (multiProcAnalysisFlag) {
1080 Cerr << "Error: modified rosenbrock direct fn does not yet support multiprocessor "
1081 << "analyses." << std::endl;
1082 abort_handler(-1);
1083 }
1084 if (numACV != 2 || numADIV > 1 || numADRV) { // allow ModelForm discrete int
1085 Cerr << "Error: Bad number of variables in modified rosenbrock direct fn."
1086 << std::endl;
1087 abort_handler(INTERFACE_ERROR);
1088 }
1089 if (numFns > 3) { // 1 fn -> opt, 2 fns -> least sq
1090 Cerr << "Error: Bad number of functions in modified rosenbrock direct fn."
1091 << std::endl;
1092 abort_handler(INTERFACE_ERROR);
1093 }
1094
1095 bool least_sq_flag = (numFns > 1);
1096 Real x1 = xCM[VAR_x1], x2 = xCM[VAR_x2], f1 = x2-x1*x1, f2 = 1.-x1,
1097 f3 = std::sin(2.*PI*x1+x2);
1098
1099 if (least_sq_flag) {
1100 // **** Residual R1:
1101 if (directFnASV[0] & 1)
1102 fnVals[0] = 10.*f1;
1103 // **** Residual R2:
1104 if (directFnASV[1] & 1)
1105 fnVals[1] = f2;
1106 // **** Residual R3:
1107 if (directFnASV[2] & 1)
1108 fnVals[2] = f3;
1109
1110 // **** dR1/dx:
1111 if (directFnASV[0] & 2)
1112 for (size_t i=0; i<numDerivVars; ++i)
1113 switch (varTypeDVV[i]) {
1114 case VAR_x1: fnGrads[0][i] = -20.*x1; break;
1115 case VAR_x2: fnGrads[0][i] = 10.; break;
1116 }
1117 // **** dR2/dx:
1118 if (directFnASV[1] & 2)
1119 for (size_t i=0; i<numDerivVars; ++i)
1120 switch (varTypeDVV[i]) {
1121 case VAR_x1: fnGrads[1][i] = -1.; break;
1122 case VAR_x2: fnGrads[1][i] = 0.; break;
1123 }
1124 // **** dR3/dx:
1125 if (directFnASV[2] & 2)
1126 for (size_t i=0; i<numDerivVars; ++i)
1127 switch (varTypeDVV[i]) {
1128 case VAR_x1: fnGrads[2][i] = 2.*PI*std::cos(2.*PI*x1+x2);break;
1129 case VAR_x2: fnGrads[2][i] = std::cos(2.*PI*x1+x2); break;
1130 }
1131
1132 // **** d^2R1/dx^2:
1133 if (directFnASV[0] & 4)
1134 for (size_t i=0; i<numDerivVars; ++i)
1135 for (size_t j=0; j<=i; ++j)
1136 if (varTypeDVV[i] == VAR_x1 && varTypeDVV[j] == VAR_x1)
1137 fnHessians[0](i,j) = -20.;
1138 else
1139 fnHessians[0](i,j) = 0.;
1140 // **** d^2R2/dx^2:
1141 if (directFnASV[1] & 4)
1142 fnHessians[1] = 0.;
1143 }
1144 else {
1145 // **** f:
1146 if (directFnASV[0] & 1)
1147 fnVals[0] = 100.*f1*f1+f2*f2+f3*f3;
1148
1149 // **** df/dx:
1150 if (directFnASV[0] & 2)
1151 for (size_t i=0; i<numDerivVars; ++i)
1152 switch (varTypeDVV[i]) {
1153 case VAR_x1: fnGrads[0][i] = -400.*f1*x1 - 2.*f2+
1154 2.*PI*std::sin(2.*(2.*PI*x1+x2));
1155 break;
1156 case VAR_x2: fnGrads[0][i] = 200.*f1+std::sin(2.*(2.*PI*x1+x2));
1157 break;
1158 }
1159
1160 // **** d^2f/dx^2:
1161 if (directFnASV[0] & 4)
1162 for (size_t i=0; i<numDerivVars; ++i)
1163 for (size_t j=0; j<=i; ++j)
1164 if (varTypeDVV[i] == VAR_x1 && varTypeDVV[j] == VAR_x1)
1165 fnHessians[0](i,j) = -400.*(x2 - 3.*x1*x1) + 2.+
1166 8.*PI*PI*std::cos(2.*(2.*PI*x1+x2));
1167 else if ( (varTypeDVV[i] == VAR_x1 && varTypeDVV[j] == VAR_x2) ||
1168 (varTypeDVV[i] == VAR_x2 && varTypeDVV[j] == VAR_x1) )
1169 fnHessians[0](i,j) = -400.*x1+4.*PI*std::cos(2.*(2.*PI*x1+x2));
1170 else if (varTypeDVV[i] == VAR_x2 && varTypeDVV[j] == VAR_x2)
1171 fnHessians[0](i,j) = 200.+2.*std::cos(2.*(2.*PI*x1+x2));
1172 }
1173
1174 return 0; // no failure
1175 }
1176
generalized_rosenbrock()1177 int TestDriverInterface::generalized_rosenbrock()
1178 {
1179 if (multiProcAnalysisFlag) {
1180 Cerr << "Error: generalized_rosenbrock direct fn does not support "
1181 << "multiprocessor analyses." << std::endl;
1182 abort_handler(-1);
1183 }
1184 if (numADIV || numADRV) {
1185 Cerr << "Error: discrete variables not supported in generalized_rosenbrock "
1186 << "direct fn." << std::endl;
1187 abort_handler(INTERFACE_ERROR);
1188 }
1189 if ( (directFnASV[0] & 6) && numVars != numDerivVars ) {
1190 Cerr << "Error: DVV subsets not supported in generalized_rosenbrock direct "
1191 << "fn." << std::endl;
1192 abort_handler(INTERFACE_ERROR);
1193 }
1194 if (numFns != 1 && numFns != 2*(numVars-1)) {
1195 Cerr << "Error: Bad number of functions in extended_rosenbrock direct fn."
1196 << std::endl;
1197 abort_handler(INTERFACE_ERROR);
1198 }
1199
1200 bool least_sq_flag = (numFns > 1);
1201
1202 // This multidimensional extension results from the overlay of numVars-1
1203 // coupled 2D Rosenbrock problems. When defining as a NLS problem,
1204 // residuals are paired and # residuals == 2 * (# vars - 1)
1205
1206 for (size_t i=1; i<numVars; ++i) {
1207 size_t index_ip1 = i, index_i = i-1; // offset by 1
1208 const Real& x_ip1 = xC[index_ip1];
1209 const Real& x_i = xC[index_i];
1210 Real f1 = x_ip1 - x_i*x_i, f2 = 1. - x_i;
1211
1212 if (least_sq_flag) {
1213 size_t rindex_2i = 2*i-1, rindex_2im1 = 2*i-2;
1214
1215 // **** R_2im1:
1216 if (directFnASV[rindex_2im1] & 1)
1217 fnVals[rindex_2im1] = 10.*f1;
1218 // **** R_2i:
1219 if (directFnASV[rindex_2i] & 1)
1220 fnVals[rindex_2i] = f2;
1221
1222 // *** dR_2im1/dx:
1223 if (directFnASV[rindex_2im1] & 2) { // define non-zeros (set_local_data)
1224 Real* grad = fnGrads[rindex_2im1];
1225 grad[index_i] = -20.*x_i;
1226 grad[index_ip1] = 10.;
1227 }
1228 // **** dR_2i/dx:
1229 if (directFnASV[rindex_2i] & 2) { // define non-zeros (set_local_data)
1230 Real* grad = fnGrads[rindex_2i];
1231 grad[index_i] = -1.;
1232 //grad[index_ip1] = 0.;
1233 }
1234
1235 // **** d^2R_2im1/dx^2:
1236 if (directFnASV[rindex_2im1] & 4) // define non-zeros (set_local_data)
1237 fnHessians[rindex_2im1](index_i,index_i) = -20.;
1238 // **** d^2R_2i/dx^2:
1239 if (directFnASV[rindex_2i] & 4)
1240 fnHessians[rindex_2i] = 0.;
1241
1242 }
1243 else {
1244
1245 // **** f:
1246 if (directFnASV[0] & 1)
1247 fnVals[0] += 100.*f1*f1 + f2*f2;
1248
1249 // **** df/dx:
1250 if (directFnASV[0] & 2) {
1251 fnGrads[0][index_i] += -400.*f1*x_i - 2.*f2;
1252 fnGrads[0][index_ip1] += 200.*f1;
1253 }
1254
1255 // **** d^2f/dx^2:
1256 if (directFnASV[0] & 4) {
1257 Real fx = x_ip1 - 3.*x_i*x_i;
1258 fnHessians[0](index_i,index_i) += -400.*fx + 2.0;
1259 fnHessians[0](index_i,index_ip1) += -400.*x_i;
1260 fnHessians[0](index_ip1,index_i) += -400.*x_i;
1261 fnHessians[0](index_ip1,index_ip1) += 200.;
1262 }
1263 }
1264 }
1265
1266 return 0; // no failure
1267 }
1268
1269
extended_rosenbrock()1270 int TestDriverInterface::extended_rosenbrock()
1271 {
1272 if (multiProcAnalysisFlag) {
1273 Cerr << "Error: extended_rosenbrock direct fn does not support "
1274 << "multiprocessor analyses." << std::endl;
1275 abort_handler(-1);
1276 }
1277 if (numADIV || numADRV) {
1278 Cerr << "Error: discrete variables not supported in extended_rosenbrock "
1279 << "direct fn." << std::endl;
1280 abort_handler(INTERFACE_ERROR);
1281 }
1282 if ( (directFnASV[0] & 6) && numVars != numDerivVars ) {
1283 Cerr << "Error: DVV subsets not supported in extended_rosenbrock direct fn."
1284 << std::endl;
1285 abort_handler(INTERFACE_ERROR);
1286 }
1287 if (numVars % 2) {
1288 Cerr << "Error: Bad number of variables in extended_rosenbrock direct fn."
1289 << std::endl;
1290 abort_handler(INTERFACE_ERROR);
1291 }
1292 if (numFns != 1 && numFns != numVars) {
1293 Cerr << "Error: Bad number of functions in extended_rosenbrock direct fn."
1294 << std::endl;
1295 abort_handler(INTERFACE_ERROR);
1296 }
1297
1298 // This multidimensional extension results from the sum of numVars/2
1299 // uncoupled 2D Rosenbrock problems. When defining as a NLS problem,
1300 // residuals are paired and # residuals == # vars
1301
1302 bool least_sq_flag = (numFns > 1);
1303 Real alpha = 100., sqrt_alpha = 10.;
1304 for (size_t i=1; i<=numVars/2; ++i) {
1305 size_t index_2i = 2*i-1, index_2im1 = 2*i-2; // offset by 1
1306 const Real& x_2i = xC[index_2i];
1307 const Real& x_2im1 = xC[index_2im1];
1308 Real f1 = x_2i - x_2im1*x_2im1, f2 = 1. - x_2im1;
1309
1310 if (least_sq_flag) {
1311
1312 // **** Residual R_2im1:
1313 if (directFnASV[index_2im1] & 1)
1314 fnVals[index_2im1] = sqrt_alpha*f1;
1315 // **** Residual R_2i:
1316 if (directFnASV[index_2i] & 1)
1317 fnVals[index_2i] = f2;
1318
1319 // *** dR_2im1/dx:
1320 if (directFnASV[index_2im1] & 2) { // define non-zeros (set_local_data)
1321 Real* grad = fnGrads[index_2im1];
1322 grad[index_2im1] = -2.*sqrt_alpha*x_2im1;
1323 grad[index_2i] = sqrt_alpha;
1324 }
1325 // **** dR_2i/dx:
1326 if (directFnASV[index_2i] & 2) { // define non-zeros (set_local_data)
1327 Real* grad = fnGrads[index_2i];
1328 grad[index_2im1] = -1.;
1329 //grad[index_2i] = 0.;
1330 }
1331
1332 // **** d^2R_2im1/dx^2:
1333 if (directFnASV[index_2im1] & 4) // define non-zeros (set_local_data)
1334 fnHessians[index_2im1](index_2im1,index_2im1) = -2.*sqrt_alpha;
1335 // **** d^2R_2i/dx^2:
1336 if (directFnASV[index_2i] & 4)
1337 fnHessians[index_2i] = 0.;
1338
1339 }
1340 else {
1341 // **** f:
1342 if (directFnASV[0] & 1)
1343 fnVals[0] += alpha*f1*f1 + f2*f2;
1344
1345 // **** df/dx:
1346 if (directFnASV[0] & 2) {
1347 fnGrads[0][index_2im1] += -4.*alpha*f1*x_2im1 - 2.*f2;
1348 fnGrads[0][index_2i] += 2.*alpha*f1;
1349 }
1350
1351 // **** d^2f/dx^2:
1352 if (directFnASV[0] & 4) {
1353 Real fx = x_2i - 3.*x_2im1*x_2im1;
1354 fnHessians[0](index_2im1,index_2im1) += -4.*alpha*fx + 2.0;
1355 fnHessians[0](index_2im1,index_2i) += -4.*alpha*x_2im1;
1356 fnHessians[0](index_2i,index_2im1) += -4.*alpha*x_2im1;
1357 fnHessians[0](index_2i,index_2i) += 2.*alpha;
1358 }
1359 }
1360 }
1361
1362 return 0; // no failure
1363 }
1364
1365
lf_rosenbrock()1366 int TestDriverInterface::lf_rosenbrock()
1367 {
1368 if (multiProcAnalysisFlag) {
1369 Cerr << "Error: lf_rosenbrock direct fn does not support "
1370 << "multiprocessor analyses." << std::endl;
1371 abort_handler(-1);
1372 }
1373 if (numACV != 2 || numADIV > 1 || numADRV) { // allow ModelForm discrete int
1374 Cerr << "Error: Bad number of variables in lf_rosenbrock direct fn."
1375 << std::endl;
1376 abort_handler(INTERFACE_ERROR);
1377 }
1378 if (numFns > 1) {
1379 Cerr << "Error: Bad number of functions in lf_rosenbrock direct fn."
1380 << std::endl;
1381 abort_handler(INTERFACE_ERROR);
1382 }
1383
1384 Real x1 = xCM[VAR_x1], x2 = xCM[VAR_x2],
1385 f1 = x2 - x1*x1 + .2, f2 = 0.8 - x1; // terms offset by +/- .2
1386
1387 // **** f:
1388 if (directFnASV[0] & 1)
1389 fnVals[0] = 100.*f1*f1+f2*f2;
1390
1391 // **** df/dx:
1392 if (directFnASV[0] & 2)
1393 for (size_t i=0; i<numDerivVars; ++i)
1394 switch (varTypeDVV[i]) {
1395 case VAR_x1: fnGrads[0][i] = -400.*f1*x1 - 2.*f2; break;
1396 case VAR_x2: fnGrads[0][i] = 200.*f1; break;
1397 }
1398
1399 // **** d^2f/dx^2:
1400 if (directFnASV[0] & 4)
1401 for (size_t i=0; i<numDerivVars; ++i)
1402 for (size_t j=0; j<=i; ++j)
1403 if (varTypeDVV[i] == VAR_x1 && varTypeDVV[j] == VAR_x1)
1404 fnHessians[0](i,j) = -400.*(x2 - 3.*x1*x1 + .2) + 2.;
1405 else if ( (varTypeDVV[i] == VAR_x1 && varTypeDVV[j] == VAR_x2) ||
1406 (varTypeDVV[i] == VAR_x2 && varTypeDVV[j] == VAR_x1) )
1407 fnHessians[0](i,j) = -400.*x1;
1408 else if (varTypeDVV[i] == VAR_x2 && varTypeDVV[j] == VAR_x2)
1409 fnHessians[0](i,j) = 200.;
1410
1411 return 0; // no failure
1412 }
1413
1414 // extra low fidelity rosenbrock
extra_lf_rosenbrock()1415 int TestDriverInterface::extra_lf_rosenbrock()
1416 {
1417 if (multiProcAnalysisFlag) {
1418 Cerr << "Error: extra_lf_rosenbrock direct fn does not support "
1419 << "multiprocessor analyses." << std::endl;
1420 abort_handler(-1);
1421 }
1422 if (numACV != 2 || numADIV > 1 || numADRV) { // allow ModelForm discrete int
1423 Cerr << "Error: Bad number of variables in lf_rosenbrock direct fn."
1424 << std::endl;
1425 abort_handler(INTERFACE_ERROR);
1426 }
1427 if (numFns > 1) {
1428 Cerr << "Error: Bad number of functions in extra_lf_rosenbrock direct fn."
1429 << std::endl;
1430 abort_handler(INTERFACE_ERROR);
1431 }
1432
1433 Real a = 0.92;
1434 Real b = -1.0;
1435 Real c = 1.1;
1436 Real d = 1.1;
1437
1438 Real x1 = xCM[VAR_x1], x2 = xCM[VAR_x2],
1439 f1 = x2 - a*x1*x1 + b, f2 = c - d*x1;
1440
1441 // **** f:
1442 if (directFnASV[0] & 1)
1443 fnVals[0] = 100.*f1*f1+f2*f2;
1444
1445 // **** df/dx:
1446 if (directFnASV[0] & 2)
1447 for (size_t i=0; i<numDerivVars; ++i)
1448 switch (varTypeDVV[i]) {
1449 case VAR_x1: fnGrads[0][i] = -400.*a*f1*x1 - 2.*d*f2; break;
1450 case VAR_x2: fnGrads[0][i] = 200.*f1; break;
1451 }
1452
1453 // **** d^2f/dx^2:
1454 if (directFnASV[0] & 4)
1455 for (size_t i=0; i<numDerivVars; ++i)
1456 for (size_t j=0; j<=i; ++j)
1457 if (varTypeDVV[i] == VAR_x1 && varTypeDVV[j] == VAR_x1)
1458 fnHessians[0](i,j) = -400.*a*f1 + 400.*a*a*x1*x1 + 2.*d*d;
1459 else if ( (varTypeDVV[i] == VAR_x1 && varTypeDVV[j] == VAR_x2) ||
1460 (varTypeDVV[i] == VAR_x2 && varTypeDVV[j] == VAR_x1) )
1461 fnHessians[0](i,j) = -400.*a*x1;
1462 else if (varTypeDVV[i] == VAR_x2 && varTypeDVV[j] == VAR_x2)
1463 fnHessians[0](i,j) = 200.;
1464
1465 return 0; // no failure
1466 }
1467
1468
mf_rosenbrock()1469 int TestDriverInterface::mf_rosenbrock()
1470 {
1471 if (multiProcAnalysisFlag) {
1472 Cerr << "Error: mf_rosenbrock direct fn does not support "
1473 << "multiprocessor analyses." << std::endl;
1474 abort_handler(-1);
1475 }
1476 if (numVars != 3 || numADRV) {
1477 Cerr << "Error: Bad number of variables in mf_rosenbrock direct fn."
1478 << std::endl;
1479 abort_handler(INTERFACE_ERROR);
1480 }
1481 if (numFns > 1) {
1482 Cerr << "Error: Bad number of functions in mf_rosenbrock direct fn."
1483 << std::endl;
1484 abort_handler(INTERFACE_ERROR);
1485 }
1486
1487 switch (xDIM[VAR_MForm]) {
1488 case 1: rosenbrock(); break;
1489 case 2: lf_rosenbrock(); break;
1490 default: return 1; break;
1491 }
1492
1493 return 0; // no failure
1494 }
1495
1496 /// modified lo-fi Rosenbrock to test SBO with hierarchical approximations
lf_poly_prod()1497 int TestDriverInterface::lf_poly_prod()
1498 {
1499 if (multiProcAnalysisFlag) {
1500 Cerr << "Error: lf_poly_prod direct fn does not yet support multiprocessor "
1501 << "analyses." << std::endl;
1502 abort_handler(-1);
1503 }
1504
1505 if ( (gradFlag || hessFlag) && (numADIV || numADRV) ) {
1506 Cerr << "Error: lf_poly_prod direct fn assumes no discrete variables in "
1507 << "derivative or hessian mode." << std::endl;
1508 abort_handler(INTERFACE_ERROR);
1509 }
1510
1511 if ( numACV != 2) {
1512 Cerr << "Error: Bad number of variables in lf_poly_prod direct fn."
1513 << std::endl;
1514 abort_handler(INTERFACE_ERROR);
1515 }
1516
1517 if ( numFns != 1 ) {
1518 Cerr << "Error: Bad number of functions in lf_poly_prod direct fn."
1519 << std::endl;
1520 abort_handler(INTERFACE_ERROR);
1521 }
1522
1523 if (directFnASV[0] & 1)
1524 fnVals[0] = xC[0]*xC[0] - 0.5*xC[1];
1525
1526 if (directFnASV[0] & 2) {
1527 fnGrads[0][0] = 2.0*xC[0];
1528 fnGrads[0][1] = -0.5;
1529 }
1530
1531 if (directFnASV[0] & 4)
1532 fnHessians[0](0,0) = 2.0;
1533
1534 return 0;
1535 }
1536
1537 /// modified lo-fi Rosenbrock to test SBO with hierarchical approximations
poly_prod()1538 int TestDriverInterface::poly_prod()
1539 {
1540 if (multiProcAnalysisFlag) {
1541 Cerr << "Error: poly_prod direct fn does not yet support multiprocessor "
1542 << "analyses." << std::endl;
1543 abort_handler(-1);
1544 }
1545
1546 if ( (gradFlag || hessFlag) && (numADIV || numADRV) ) {
1547 Cerr << "Error: poly_prod direct fn assumes no discrete variables in "
1548 << "derivative or hessian mode." << std::endl;
1549 abort_handler(INTERFACE_ERROR);
1550 }
1551
1552 if ( numACV != 2) {
1553 Cerr << "Error: Bad number of variables in poly_prod direct fn."
1554 << std::endl;
1555 abort_handler(INTERFACE_ERROR);
1556 }
1557
1558 if ( numFns != 1 ) {
1559 Cerr << "Error: Bad number of functions in poly_prod direct fn."
1560 << std::endl;
1561 abort_handler(INTERFACE_ERROR);
1562 }
1563
1564 // Compute and output responses
1565 double x0_sq = xC[0]*xC[0], x1_sq = xC[1]*xC[1];
1566 double t2 = x0_sq - 0.5*xC[1];
1567 double t1 = xC[0] + x1_sq/2.;
1568
1569 if (directFnASV[0] & 1)
1570 fnVals[0] = t1*t2;
1571
1572 if (directFnASV[0] & 2) {
1573 fnGrads[0][0] = 2.0*xC[0]*t1 + t2;
1574 fnGrads[0][1] = t2*xC[1] - t1/2.0;
1575 }
1576
1577 if (directFnASV[0] & 4) {
1578 fnHessians[0](0,0) = 4.0*xC[0] + 2.*t1;
1579 fnHessians[0](1,1) = t2 - xC[1];
1580 fnHessians[0](1,0) = 2.*xC[0]*xC[1] - 0.5;
1581 }
1582 return 0;
1583 }
1584
1585
gerstner()1586 int TestDriverInterface::gerstner()
1587 {
1588 if (multiProcAnalysisFlag) {
1589 Cerr << "Error: gerstner direct fn does not support multiprocessor "
1590 << "analyses." << std::endl;
1591 abort_handler(-1);
1592 }
1593 if (numVars != 2 || numADIV || numADRV || (gradFlag && numDerivVars != 2)) {
1594 Cerr << "Error: Bad number of variables in gerstner direct fn."<< std::endl;
1595 abort_handler(INTERFACE_ERROR);
1596 }
1597 if (numFns != 1) {
1598 Cerr << "Error: Bad number of functions in gerstner direct fn."<< std::endl;
1599 abort_handler(INTERFACE_ERROR);
1600 }
1601 if (hessFlag) {
1602 Cerr << "Error: Hessians not supported in gerstner direct fn." << std::endl;
1603 abort_handler(INTERFACE_ERROR);
1604 }
1605
1606 const Real& x = xC[0]; const Real& y = xC[1];
1607 String an_comp = (!analysisComponents.empty() &&
1608 !analysisComponents[analysisDriverIndex].empty()) ?
1609 analysisComponents[analysisDriverIndex][0] : "iso1";
1610 short test_fn; Real x_coeff, y_coeff, xy_coeff;
1611 if (an_comp == "iso1")
1612 { test_fn = 1; x_coeff = y_coeff = 10.; }
1613 else if (an_comp == "iso2")
1614 { test_fn = 2; x_coeff = y_coeff = xy_coeff = 1.; }
1615 else if (an_comp == "iso3")
1616 { test_fn = 3; x_coeff = y_coeff = 10.; }
1617 else if (an_comp == "aniso1")
1618 { test_fn = 1; x_coeff = 1.; y_coeff = 10.; }
1619 else if (an_comp == "aniso2")
1620 { test_fn = 2; x_coeff = 1.; y_coeff = xy_coeff = 10.; }
1621 else if (an_comp == "aniso3")
1622 { test_fn = 3; x_coeff = 10.; y_coeff = 5.; }
1623 else {
1624 Cerr << "Error: analysis component specification required in gerstner "
1625 << "direct fn." << std::endl;
1626 abort_handler(INTERFACE_ERROR);
1627 }
1628
1629 // **** f:
1630 if (directFnASV[0] & 1) {
1631 switch (test_fn) {
1632 case 1:
1633 fnVals[0] = x_coeff*std::exp(-x*x) + y_coeff*std::exp(-y*y); break;
1634 case 2:
1635 fnVals[0] = x_coeff*std::exp(x) + y_coeff*std::exp(y)
1636 + xy_coeff*std::exp(x*y); break;
1637 case 3:
1638 fnVals[0] = std::exp(-x_coeff*x*x - y_coeff*y*y); break;
1639 }
1640 }
1641
1642 // **** df/dx:
1643 if (directFnASV[0] & 2) {
1644 Real val;
1645 switch (test_fn) {
1646 case 1:
1647 fnGrads[0][0] = -2.*x*x_coeff*std::exp(-x*x);
1648 fnGrads[0][1] = -2.*y*y_coeff*std::exp(-y*y); break;
1649 case 2:
1650 val = xy_coeff*std::exp(x*y);
1651 fnGrads[0][0] = x_coeff*std::exp(x) + val*y;
1652 fnGrads[0][1] = y_coeff*std::exp(y) + val*x; break;
1653 case 3:
1654 val = std::exp(-x_coeff*x*x - y_coeff*y*y);
1655 fnGrads[0][0] = -2.*x*x_coeff*val;
1656 fnGrads[0][1] = -2.*y*y_coeff*val; break;
1657 }
1658 }
1659
1660 return 0; // no failure
1661 }
1662
scalable_gerstner()1663 int TestDriverInterface::scalable_gerstner()
1664 {
1665 if (multiProcAnalysisFlag) {
1666 Cerr << "Error: scalable_gerstner direct fn does not support "
1667 << "multiprocessor analyses." << std::endl;
1668 abort_handler(-1);
1669 }
1670 if (numADIV || numADRV) {
1671 Cerr << "Error: Bad variable types in scalable_gerstner direct fn."
1672 << std::endl;
1673 abort_handler(INTERFACE_ERROR);
1674 }
1675 if (numFns != 1) {
1676 Cerr << "Error: Bad number of functions in scalable_gerstner direct fn."
1677 << std::endl;
1678 abort_handler(INTERFACE_ERROR);
1679 }
1680 if (hessFlag) {
1681 Cerr << "Error: Hessians not supported in scalable_gerstner direct fn."
1682 << std::endl;
1683 abort_handler(INTERFACE_ERROR);
1684 }
1685
1686 String an_comp = (!analysisComponents.empty() &&
1687 !analysisComponents[analysisDriverIndex].empty()) ?
1688 analysisComponents[analysisDriverIndex][0] : "iso1";
1689 short test_fn; Real even_coeff, odd_coeff, inter_coeff;
1690 if (an_comp == "iso1")
1691 { test_fn = 1; even_coeff = odd_coeff = 10.; }
1692 else if (an_comp == "iso2")
1693 { test_fn = 2; even_coeff = odd_coeff = inter_coeff = 1.; }
1694 else if (an_comp == "iso3")
1695 { test_fn = 3; even_coeff = odd_coeff = 10.; }
1696 else if (an_comp == "aniso1")
1697 { test_fn = 1; even_coeff = 1.; odd_coeff = 10.; }
1698 else if (an_comp == "aniso2")
1699 { test_fn = 2; even_coeff = 1.; odd_coeff = inter_coeff = 10.; }
1700 else if (an_comp == "aniso3")
1701 { test_fn = 3; even_coeff = 10.; odd_coeff = 5.; }
1702 else {
1703 Cerr << "Error: analysis component specification required in gerstner "
1704 << "direct fn." << std::endl;
1705 abort_handler(INTERFACE_ERROR);
1706 }
1707
1708 // **** f:
1709 if (directFnASV[0] & 1) {
1710 switch (test_fn) {
1711 case 1:
1712 fnVals[0] = 0.;
1713 for (size_t i=0; i<numVars; ++i)
1714 fnVals[0] += (i%2) ? odd_coeff*std::exp(-xC[i]*xC[i]) :
1715 even_coeff*std::exp(-xC[i]*xC[i]); break;
1716 case 2:
1717 fnVals[0] = 0.;
1718 for (size_t i=0; i<numVars; ++i)
1719 if (i%2)
1720 fnVals[0] += odd_coeff*std::exp(xC[i])
1721 + inter_coeff*std::exp(xC[i]*xC[i-1]);
1722 else
1723 fnVals[0] += even_coeff*std::exp(xC[i]);
1724 break;
1725 case 3: {
1726 Real sum = 0;
1727 for (size_t i=0; i<numVars; ++i)
1728 sum -= (i%2) ? odd_coeff*xC[i]*xC[i] : even_coeff*xC[i]*xC[i];
1729 fnVals[0] = std::exp(sum); break;
1730 }
1731 }
1732 }
1733
1734 // **** df/dx:
1735 if (directFnASV[0] & 2) {
1736 Real val;
1737 switch (test_fn) {
1738 case 1:
1739 for (size_t i=0; i<numVars; ++i)
1740 fnGrads[0][i] = (i%2) ? -2.*xC[i]* odd_coeff*std::exp(-xC[i]*xC[i])
1741 : -2.*xC[i]*even_coeff*std::exp(-xC[i]*xC[i]);
1742 break;
1743 case 2:
1744 for (size_t i=0; i<numVars; ++i)
1745 {
1746 if (i%2)
1747 fnGrads[0][i] = odd_coeff*std::exp(xC[i])
1748 + inter_coeff*xC[i-1]*std::exp(xC[i]*xC[i-1]);
1749 else
1750 {
1751 fnGrads[0][i] = even_coeff*std::exp(xC[i]);
1752 if ( i+1 < numVars )
1753 fnGrads[0][i] += inter_coeff*xC[i+1]*std::exp(xC[i]*xC[i+1]);
1754 }
1755 }
1756 break;
1757 case 3:
1758 if (directFnASV[0] & 1)
1759 val = fnVals[0];
1760 else {
1761 Real sum = 0;
1762 for (size_t i=0; i<numVars; ++i)
1763 sum -= (i%2) ? odd_coeff*xC[i]*xC[i] : even_coeff*xC[i]*xC[i];
1764 val = std::exp(sum);
1765 }
1766 for (size_t i=0; i<numVars; ++i)
1767 fnGrads[0][i] = (i%2) ? -2.*xC[i]* odd_coeff*val
1768 : -2.*xC[i]*even_coeff*val;
1769 break;
1770 }
1771 }
1772
1773 return 0; // no failure
1774 }
1775
1776
1777 void TestDriverInterface::
get_genz_coefficients(int num_dims,Real factor,int c_type,RealVector & c,RealVector & w)1778 get_genz_coefficients( int num_dims, Real factor, int c_type,
1779 RealVector &c, RealVector &w )
1780 {
1781 c.resize( num_dims );
1782 w.resize( num_dims );
1783 switch ( c_type )
1784 {
1785 case 0:
1786 {
1787 Real csum = 0.0;
1788 for ( int d = 0; d < num_dims; d++ )
1789 {
1790 w[d] = 0.0;
1791 c[d] = ( (Real)d + 0.5 ) / (Real)num_dims;
1792 csum += c[d];
1793 }
1794 for ( int d = 0; d < num_dims; d++ )
1795 {
1796 c[d] *= ( factor / csum );
1797 }
1798 break;
1799 }
1800 case 1:
1801 {
1802 Real csum = 0.0;
1803 for ( int d = 0; d < num_dims; d++ )
1804 {
1805 w[d] = 0.0;
1806 c[d] = 1.0 / (Real)( ( d + 1 ) * ( d + 1 ) );
1807 csum += c[d];
1808 }
1809 for ( int d = 0; d < num_dims; d++ )
1810 {
1811 c[d] *= ( factor / csum );
1812 }
1813 break;
1814 }
1815 case 2:
1816 {
1817 Real csum = 0.0;
1818 for ( int d = 0; d < num_dims; d++ )
1819 {
1820 w[d] = 0.;
1821 c[d] = std::exp( (Real)(d+1)*std::log( 1.e-8 ) / (Real)num_dims );
1822 csum += c[d];
1823 }
1824 for ( int d = 0; d < num_dims; d++ )
1825 {
1826 c[d] *= ( factor / csum );
1827 }
1828 break;
1829 }
1830 default:
1831 throw(std::runtime_error("GetCoefficients() ensure type in [0,1]"));
1832 }
1833 }
1834
1835
1836 void TestDriverInterface::
steady_state_diffusion_core(SpectralDiffusionModel & model,RealVector & domain_limits)1837 steady_state_diffusion_core(SpectralDiffusionModel& model,
1838 RealVector& domain_limits)
1839 {
1840 // ------------------------------------------------------------- //
1841 // Pre-processing
1842 // ------------------------------------------------------------- //
1843
1844 bool err_flag = false;
1845 if (multiProcAnalysisFlag) {
1846 Cerr << "Error: steady_state_diffusion_1d direct fn does not support "
1847 << "multiprocessor analyses." << std::endl;
1848 err_flag = true;
1849 }
1850 if ( ( numVars < 1 ) || ( numADIV > 1 ) ) {
1851 Cerr << "Error: Bad variable types in steady_state_diffusion_1d direct fn."
1852 << std::endl;
1853 err_flag = true;
1854 }
1855 if (numFns < 1) {
1856 Cerr << "Error: Bad number of functions in steady_state_diffusion_1d "
1857 << "direct fn." << std::endl;
1858 err_flag = true;
1859 }
1860 if (hessFlag||gradFlag) {
1861 Cerr << "Error: Gradients and Hessians are not supported in "
1862 << "steady_state_diffusion_1d direct fn." << std::endl;
1863 err_flag = true;
1864 }
1865 if (err_flag)
1866 abort_handler(INTERFACE_ERROR);
1867
1868 // ------------------------------------------------------------- //
1869 // Read parameters from discrete state variables
1870 // ------------------------------------------------------------- //
1871
1872 // Get the positivity flag from the discrete string variables
1873 size_t pos_index = find_index(xDSLabels, "positivity");
1874 String pos_string = ( pos_index == _NPOS ) ? "on" : xDS[pos_index];
1875 bool positivity = ( pos_string == "on" );
1876
1877 // Get the field mean from the discrete real variables
1878 size_t mean_index = find_index(xDRLabels, "field_mean");
1879 Real field_mean = ( mean_index == _NPOS ) ? 1.0 : xDR[mean_index];
1880
1881 // Get the field mean from the discrete real variables
1882 size_t std_dev_index = find_index(xDRLabels, "field_std_dev");
1883 Real field_std_dev = ( std_dev_index == _NPOS ) ? 1.0 : xDR[std_dev_index];
1884
1885 // Get the kernel order from the discrete real variables
1886 size_t kern_ord_index = find_index(xDRLabels, "kernel_order");
1887 Real kernel_order = ( kern_ord_index == _NPOS ) ? 1.0 : xDR[kern_ord_index];
1888
1889 // Get the kernel length from the discrete real variables
1890 size_t kern_len_index = find_index(xDRLabels, "kernel_length");
1891 Real kernel_length = ( kern_len_index == _NPOS ) ? 1.0 : xDR[kern_len_index];
1892
1893 // Compute default QoI coordinates:
1894 RealVector qoi_coords( numFns, false );
1895 if (numFns > 1) {
1896 Real range = domain_limits[1]-domain_limits[0];
1897 Real h = (range*0.9) / (Real)(numFns-1);
1898 for (int i=0; i<numFns; i++)
1899 qoi_coords[i] = domain_limits[0]+range*0.05+(Real)i*h;
1900 }
1901 else
1902 qoi_coords[0] = (domain_limits[1]+domain_limits[0])/2.;
1903
1904 // If QoI coordinates provided through discrete real variables, overwrite
1905 // defaults:
1906 for (int i=0; i<numFns; i++) {
1907 size_t coord_index
1908 = find_index(xDRLabels, "coord_" + std::to_string(i));
1909 if ( coord_index != _NPOS )
1910 qoi_coords[i] = xDR[coord_index];
1911 }
1912
1913 // ------------------------------------------------------------- //
1914 // Initialize model
1915 // ------------------------------------------------------------- //
1916
1917 model.set_num_qoi( numFns );
1918 model.set_qoi_coords( qoi_coords );
1919 model.set_field_mean( field_mean );
1920 model.set_field_std_dev( field_std_dev );
1921
1922 // These are ignored if kernel is not exponential
1923 model.set_positivity( positivity );
1924 model.set_kernel_order( kernel_order );
1925 model.set_kernel_length( kernel_length );
1926 }
1927
1928
1929 /** \brief Solve the 1D diffusion equation with an uncertain variable
1930 * coefficient using the spectral Chebyshev collocation method.
1931 *
1932 * del(k del(u) ) = f on [0,1] subject to u(0) = 0 u(1) = 0
1933 *
1934 * Here we set f = -1 and
1935 * k = 1+4.*sum_d [cos(2*pi*x)/(pi*d)^2*z[d]] d=1,...,num_dims
1936 * where z_d are random variables, typically i.i.d uniform[-1,1]
1937 */
steady_state_diffusion_1d()1938 int TestDriverInterface::steady_state_diffusion_1d()
1939 {
1940 // Initialize domain and boundary conditions:
1941 RealVector bndry_conds(2), domain_limits(2); // initialize to zero
1942 domain_limits[1] = 1.;
1943
1944 SpectralDiffusionModel model;
1945 steady_state_diffusion_core(model, domain_limits);
1946
1947 // ------------------------------------------------------------- //
1948 // Read parameters from discrete state variables
1949 // ------------------------------------------------------------- //
1950
1951 // Get the mesh resolution from the first discrete integer variable
1952 size_t mesh_size_index = find_index(xDILabels, "mesh_size");
1953 int order = ( mesh_size_index == _NPOS ) ? 20 : xDI[mesh_size_index];
1954
1955 // Get the kernel specification from the discrete string variables
1956 size_t kernel_index = find_index(xDSLabels, "kernel_type");
1957 String kernel = ( kernel_index == _NPOS ) ? "default" : xDS[kernel_index];
1958
1959 if (order % 2) {
1960 Cerr << "Error: Mesh size must be even." << std::endl;
1961 abort_handler(INTERFACE_ERROR);
1962 }
1963 if (order + 1 < xC.length() && kernel == "exponential") {
1964 Cerr << "Error: Mesh size must be greater than or equal "
1965 << "to the number of random variables + 1 when using "
1966 << "the exponential kernel." << std::endl;
1967 abort_handler(INTERFACE_ERROR);
1968 }
1969
1970 // ------------------------------------------------------------- //
1971 // Initialize and evaluate model
1972 // ------------------------------------------------------------- //
1973
1974 model.initialize( order, kernel, bndry_conds, domain_limits );
1975 model.evaluate( xC, fnVals );
1976 return 0;
1977 }
1978
1979
ss_diffusion_discrepancy()1980 int TestDriverInterface::ss_diffusion_discrepancy()
1981 {
1982 // Initialize domain and boundary conditions:
1983 RealVector bndry_conds(2), domain_limits(2); // initialize to zero
1984 domain_limits[1] = 1.;
1985
1986 SpectralDiffusionModel model;
1987 steady_state_diffusion_core(model, domain_limits);
1988
1989 // ------------------------------------------------------------- //
1990 // Read parameters from discrete state variables
1991 // ------------------------------------------------------------- //
1992
1993 // Get the mesh resolution from the first discrete integer variable
1994 size_t mesh_size_index = find_index(xDILabels, "mesh_size");
1995 int order_l = ( mesh_size_index == _NPOS ) ? 20 : xDI[mesh_size_index];
1996 int order_lm1 = order_l / 2;
1997 bool err_flag = false;
1998 if (order_l % 2)
1999 { Cerr << "Error: mesh size must be even." << std::endl; err_flag = true; }
2000 else if (order_l < 4) {
2001 Cerr << "Error: mesh size must be at least 4 at level l for even mesh "
2002 << "size and level l-1." << std::endl;
2003 err_flag = true;
2004 }
2005
2006 // Get the kernel specification from the discrete string variables
2007 size_t kernel_index = find_index(xDSLabels, "kernel_type");
2008 String kernel = ( kernel_index == _NPOS ) ? "default" : xDS[kernel_index];
2009
2010 if (order_lm1 + 1 < xC.length() && kernel == "exponential") {
2011 Cerr << "Error: mesh size must be >= the number of random variables + 1 "
2012 << "when using the exponential kernel." << std::endl;
2013 err_flag = true;
2014 }
2015
2016 if (err_flag)
2017 abort_handler(INTERFACE_ERROR);
2018
2019 // ------------------------------------------------------------- //
2020 // Evaluate model twice to compute discrepancy across consecutive
2021 // mesh resolutions
2022 // ------------------------------------------------------------- //
2023
2024 model.initialize( order_l, kernel, bndry_conds, domain_limits );
2025 model.evaluate( xC, fnVals );
2026
2027 RealVector q_lm1(numFns, false);
2028 model.initialize( order_lm1, kernel, bndry_conds, domain_limits );
2029 model.evaluate( xC, q_lm1 );
2030 fnVals -= q_lm1;
2031
2032 return 0;
2033 }
2034
2035
transient_diffusion_1d()2036 int TestDriverInterface::transient_diffusion_1d()
2037 {
2038 if (multiProcAnalysisFlag) {
2039 Cerr << "Error: transient_diffusion_1d direct fn does not support "
2040 << "multiprocessor analyses." << std::endl;
2041 abort_handler(-1);
2042 }
2043 if (numACV != 7 || numADIV > 2) {
2044 Cerr << "Error: unsupported variable counts in transient_diffusion_1d "
2045 << "direct fn." << std::endl;
2046 abort_handler(INTERFACE_ERROR);
2047 }
2048 if (numFns > 1) {
2049 Cerr << "Error: unsupported function counts in transient_diffusion_1d "
2050 << "direct fn." << std::endl;
2051 abort_handler(INTERFACE_ERROR);
2052 }
2053 if (hessFlag || gradFlag) {
2054 Cerr << "Error: gradients and Hessians are not supported in "
2055 << "transient_diffusion_1d direct fn." << std::endl;
2056 abort_handler(INTERFACE_ERROR);
2057 }
2058 // Get the number of spatial discretization points and the number of
2059 // Fourier solution modes from the discrete integer variables
2060 size_t nx_index = find_index(xDILabels, "N_x"),
2061 nm_index = find_index(xDILabels, "N_mod");
2062 int N_x = ( nx_index == _NPOS ) ? 200 : xDI[nx_index];
2063 int N_mod = ( nm_index == _NPOS ) ? 21 : xDI[nm_index];
2064
2065 RealVector x(N_x+1, false), u_n(N_x+1, false), f_tilde(N_x+1, false),
2066 A(N_x, false);
2067 Real Pi = 4. * std::atan(1.);
2068
2069 // scaling from [-1,1] to [xi_min, xi_max]
2070 // [-1,1] -> [l,u]: xi = l + (u-l) * (xC+1) / 2
2071 RealVector xi(numVars);
2072 xi[0] = Pi*xC[0]; xi[1] = Pi*xC[1]; xi[2] = Pi*xC[2]; // [-pi,pi]
2073 xi[3] = 0.005 + 0.004*xC[3]; // [.001, .009]
2074 xi[4] = (xC[4]+1.)/2.; xi[5] = (xC[5]+1.)/2.; xi[6] = (xC[6]+1.)/2.; // [0,1]
2075
2076 // spatial discretization
2077 size_t i, n;
2078 Real dx = 1. / N_x;
2079 x[0] = 0.;
2080 for (i=1; i<=N_x; ++i)
2081 x[i] = x[i-1] + dx;
2082
2083 Real gamma = 1., T_f = 0.5, sum_A, int_u_n = 0.,
2084 f_t_1, f_t_2, x_i, xpi, npi_g, u_exp, xi2_sq = xi[2]*xi[2],
2085 sin_xi0 = std::sin(xi[0]), sin_xi1 = std::sin(xi[1]),
2086 i_fn = 3.5*(sin_xi0 + 7.*sin_xi1*sin_xi1 + xi2_sq*xi2_sq*sin_xi0/10.);
2087 Real g_fn = 50., a_i = -0.5;
2088 for (i=4; i<=6; ++i)
2089 g_fn *= (std::abs(4.*xi[i] - 2.) + a_i) / (1.+a_i);
2090
2091 f_tilde[0] = 0.;
2092 for (n=0; n<N_mod; ++n) {
2093 sum_A = 0.; u_n[0] = 0.; npi_g = n * Pi / gamma;
2094 u_exp = std::exp(-xi[3] * npi_g * npi_g * T_f);
2095 for (i=1; i<=N_x; ++i) {
2096 x_i = x[i]; xpi = Pi*x_i; f_t_1 = std::sin(xpi);
2097 f_t_2 = std::sin(2.*xpi) + std::sin(3.*xpi)
2098 + 50.* ( std::sin(9.*xpi) + std::sin(21.*xpi) );
2099 // This term is the product of (a realization of) the initial solution
2100 // and the n_th modal term
2101 f_tilde[i] = std::sin(npi_g * x_i) * (g_fn * f_t_1 + i_fn * f_t_2);
2102
2103 // n_th coeff of the modal expansion (integrated over a single interval)
2104 A[i-1] = (f_tilde[i] + f_tilde[i-1]) * dx / 2.; // area (trapezoidal rule)
2105 sum_A += A[i-1];
2106
2107 u_n[i] = std::sin(npi_g * x_i) * u_exp;
2108 }
2109
2110 // n_th coefficient of the modal expansion (...sum of the integrals)
2111 u_n.scale(2. * sum_A);
2112 for (i=1; i<=N_x; ++i)
2113 int_u_n += u_n[i-1] + u_n[i]; // ends counted once, interior counted twice
2114 }
2115
2116 // QoI is the integral of the solution over the physical space
2117 fnVals[0] = int_u_n * dx / 2.; // trapezoidal rule
2118
2119 return 0;
2120 }
2121
predator_prey()2122 int TestDriverInterface::predator_prey()
2123 {
2124 // ------------------------------------------------------------- //
2125 // Pre-processing
2126 // ------------------------------------------------------------- //
2127
2128 if (multiProcAnalysisFlag) {
2129 Cerr << "Error: predator_prey direct fn does not support "
2130 << "multiprocessor analyses." << std::endl;
2131 abort_handler(-1);
2132 }
2133 if ( numVars < 1 || numADIV > 1 || numADRV > 1 ) {
2134 Cerr << "Error: Bad variable types in predator_prey direct fn."<< std::endl;
2135 abort_handler(INTERFACE_ERROR);
2136 }
2137 if (numFns != 3) {
2138 Cerr << "Error: Bad number of functions in predator_prey direct fn."
2139 << std::endl;
2140 abort_handler(INTERFACE_ERROR);
2141 }
2142 if (hessFlag || gradFlag) {
2143 Cerr << "Error: Gradients and Hessians are not supported in "
2144 << "predator_prey direct fn." << std::endl;
2145 abort_handler(INTERFACE_ERROR);
2146 }
2147
2148 // ------------------------------------------------------------- //
2149 // Read parameters from discrete state variables
2150 // ------------------------------------------------------------- //
2151
2152 // Get the mesh resolution from the first discrete integer variable
2153 size_t step_index = find_index(xDILabels, "time_steps");
2154 int num_tsteps = ( step_index == _NPOS ) ? 101 : xDI[step_index];
2155 if (num_tsteps % 2 != 1) {
2156 Cerr << "Error: Number of time steps must be odd" << std::endl;
2157 abort_handler(INTERFACE_ERROR);
2158 }
2159
2160 size_t tf_index = find_index(xDRLabels, "final_time");
2161 Real final_time = ( tf_index == _NPOS ) ? 10. : xDR[tf_index];
2162
2163 RealVector init_conditions(3);
2164 init_conditions[0] = 0.7; init_conditions[1] = 0.5;
2165 init_conditions[2] = 0.2; // these could be made random variables
2166
2167 // ------------------------------------------------------------- //
2168 // Initialize and evaluate model
2169 // ------------------------------------------------------------- //
2170
2171 PredatorPreyModel model;
2172 model.set_initial_conditions( init_conditions );
2173 model.set_time(final_time, final_time/((double)num_tsteps-1.));
2174
2175 model.evaluate( xC, fnVals );
2176
2177 return 0;
2178 }
2179
genz()2180 int TestDriverInterface::genz()
2181 {
2182 if (multiProcAnalysisFlag) {
2183 Cerr << "Error: genz direct fn does not support "
2184 << "multiprocessor analyses." << std::endl;
2185 abort_handler(-1);
2186 }
2187 if (numADIV || numADRV) {
2188 Cerr << "Error: Bad variable types in genz direct fn."
2189 << std::endl;
2190 abort_handler(INTERFACE_ERROR);
2191 }
2192 if (numFns != 1) {
2193 Cerr << "Error: Bad number of functions in genz direct fn."
2194 << std::endl;
2195 abort_handler(INTERFACE_ERROR);
2196 }
2197 if (hessFlag) {
2198 Cerr << "Error: Hessians not supported in genz direct fn."
2199 << std::endl;
2200 abort_handler(INTERFACE_ERROR);
2201 }
2202
2203 String test;//, resp;
2204 if (analysisComponents.empty() ||
2205 analysisComponents[analysisDriverIndex].empty())
2206 { test = "os1"; /*resp = "generic";*/ }
2207 else {
2208 StringArray& an_comps = analysisComponents[analysisDriverIndex];
2209 test = an_comps[0];
2210 //if (an_comps.size() > 1) resp = an_comps[1];
2211 }
2212
2213 int coeff_type, fn_type;
2214 Real factor;
2215 if (test == "os1")
2216 { coeff_type = 0; fn_type = 0; factor = 4.5; }
2217 else if (test == "os2")
2218 { coeff_type = 1; fn_type = 0; factor = 4.5; }
2219 else if (test == "os3")
2220 { coeff_type = 2; fn_type = 0; factor = 4.5; }
2221 else if (test == "cp1")
2222 { coeff_type = 0; fn_type = 1; factor = .25; }
2223 else if (test == "cp2")
2224 { coeff_type = 1; fn_type = 1; factor = .25; }
2225 else if (test == "cp3")
2226 { coeff_type = 2; fn_type = 1; factor = .25; }
2227 else {
2228 Cerr << "Error: analysis component specification required in genz "
2229 << "direct fn." << std::endl;
2230 abort_handler(INTERFACE_ERROR);
2231 }
2232
2233 RealVector c, w;
2234 get_genz_coefficients( numVars, factor, coeff_type, c, w );
2235 Real pi = 4.0 * std::atan( 1.0 );
2236
2237 // **** f:
2238 if (directFnASV[0] & 1) {
2239 switch (fn_type) {
2240 case 0:
2241 fnVals[0] = 2.0 * pi * w[0];
2242 for ( int d = 0; d < numVars; d++ )
2243 fnVals[0] += c[d] * xC[d];
2244 fnVals[0] = std::cos( fnVals[0] );
2245 break;
2246 case 1:
2247 fnVals[0] = 1.0;
2248 for ( int d = 0; d < numVars; d++ )
2249 fnVals[0] += c[d]* xC[d];
2250 Real expo = -(Real)(numVars+1);
2251 //if (resp == "residual") expo /= 2.;// reproduce Genz after sum res^2
2252 fnVals[0] = std::pow( fnVals[0], expo );
2253 break;
2254 }
2255 }
2256
2257 return 0; // no failure
2258 }
2259
2260
damped_oscillator()2261 int TestDriverInterface::damped_oscillator()
2262 {
2263 if (multiProcAnalysisFlag) {
2264 Cerr << "Error: damped oscillator direct fn does not support "
2265 << "multiprocessor analyses." << std::endl;
2266 abort_handler(-1);
2267 }
2268 if (numVars < 1 || numVars > 6 || numADIV || numADRV) {
2269 Cerr << "Error: Bad variable types in damped oscillator direct fn."
2270 << std::endl;
2271 abort_handler(INTERFACE_ERROR);
2272 }
2273 if (numFns < 1) {
2274 Cerr << "Error: Bad number of functions in damped oscillator direct fn."
2275 << std::endl;
2276 abort_handler(INTERFACE_ERROR);
2277 }
2278 if (hessFlag || gradFlag) {
2279 Cerr << "Error: Gradients and Hessians not supported in damped oscillator "
2280 << "direct fn." << std::endl;
2281 abort_handler(INTERFACE_ERROR);
2282 }
2283
2284 Real initial_time = 0., final_time = 20., //delta_t = 0.3;
2285 delta_t = (final_time - initial_time) / numFns;
2286 //int num_time_steps = numFns;
2287
2288 Real pi = 4.0 * std::atan( 1.0 );
2289
2290 Real b = xC[0], k = 0.035, F = 0.1, w = 1.0, x0 = 0.5, v0 = 0.;
2291 if ( numVars >= 2 ) k = xC[1];
2292 if ( numVars >= 3 ) F = xC[2];
2293 if ( numVars >= 4 ) w = xC[3];
2294 if ( numVars >= 5 ) x0 = xC[4];
2295 if ( numVars >= 6 ) v0 = xC[5];
2296
2297 Real kw = ( k-w*w ), bw = b*w, g = b / 2.;
2298 Real zeta2 = bw*bw + kw*kw, zeta = std::sqrt(zeta2);
2299 Real phi = std::atan( -bw / kw );
2300 Real sqrtk = std::sqrt( k );
2301 Real wd = sqrtk*std::sqrt( 1.-g*g / k );
2302 if ( kw / zeta2 < 0 ) phi += pi;
2303 Real B1 = -F*( bw ) / zeta2, B2 = F*kw / zeta2,
2304 A1 = x0-B1, A2 = ( v0+g*A1-w*B2 ) / wd;
2305
2306 if ( sqrtk <= g ) {
2307 Cerr << "Error: damped_oscillator parameters do not result in under-damped "
2308 << "solution." << std::endl;
2309 abort_handler(INTERFACE_ERROR);
2310 }
2311
2312 // always include response at final_time plus additional time steps;
2313 // response at initial time isn't included as it isn't a fn of some params.
2314 Real time = initial_time, y_stead, y_trans;
2315 for (size_t i=0; i<numFns; ++i) {
2316 time += delta_t;
2317 if (directFnASV[i] & 1) {
2318 // Steady state solution (y_stead) for rhs = 0
2319 y_stead = F * std::sin( w*time + phi ) / zeta;
2320 // Compute transient (y_trans) component of solution
2321 y_trans = std::exp( -g*time )*( A1 * std::cos( wd*time ) +
2322 A2 * std::sin( wd*time ) );
2323 fnVals[i] = y_stead + y_trans;
2324 }
2325 }
2326
2327 return 0; // no failure
2328 }
2329
2330
log_ratio()2331 int TestDriverInterface::log_ratio()
2332 {
2333 if (multiProcAnalysisFlag) {
2334 Cerr << "Error: log_ratio direct fn does not support multiprocessor "
2335 << "analyses." << std::endl;
2336 abort_handler(-1);
2337 }
2338 if ( numVars != 2 || numADIV || numADRV ||
2339 ( ( gradFlag || hessFlag ) && numDerivVars != 2 ) ) {
2340 Cerr << "Error: Bad number of variables in log_ratio direct fn."
2341 << std::endl;
2342 abort_handler(INTERFACE_ERROR);
2343 }
2344 if (numFns != 1) {
2345 Cerr << "Error: Bad number of functions in log_ratio direct fn."
2346 << std::endl;
2347 abort_handler(INTERFACE_ERROR);
2348 }
2349
2350 const Real& x1 = xC[0]; const Real& x2 = xC[1];
2351
2352 // **** f:
2353 if (directFnASV[0] & 1)
2354 fnVals[0] = x1/x2;
2355
2356 // **** df/dx:
2357 if (directFnASV[0] & 2) {
2358 fnGrads[0][0] = 1./x2;
2359 fnGrads[0][1] = -x1/(x2*x2);
2360 }
2361
2362 // **** d^2f/dx^2:
2363 if (directFnASV[0] & 4) {
2364 fnHessians[0](0,0) = 0.0;
2365 fnHessians[0](0,1) = fnHessians[0](1,0) = -1./(x2*x2);
2366 fnHessians[0](1,1) = 2.*x1/std::pow(x2,3);
2367 }
2368
2369 return 0; // no failure
2370 }
2371
2372
short_column()2373 int TestDriverInterface::short_column()
2374 {
2375 if (multiProcAnalysisFlag) {
2376 Cerr << "Error: short_column direct fn does not support multiprocessor "
2377 << "analyses." << std::endl;
2378 abort_handler(-1);
2379 }
2380 if (numACV != 5 || numADIV > 1 || numADRV) { // allow ModelForm discrete int
2381 Cerr << "Error: Bad number of variables in short_column direct fn."
2382 << std::endl;
2383 abort_handler(INTERFACE_ERROR);
2384 }
2385 size_t ai, lsi;
2386 if (numFns == 1) // option for limit state only
2387 lsi = 0;
2388 else if (numFns == 2) // option for area + limit state
2389 { ai = 0; lsi = 1; }
2390 else {
2391 Cerr << "Error: Bad number of functions in short_column direct fn."
2392 << std::endl;
2393 abort_handler(INTERFACE_ERROR);
2394 }
2395
2396 // b = xC[0] = column base (design var.)
2397 // h = xC[1] = column height (design var.)
2398 // P = xC[2] (normal uncertain var.)
2399 // M = xC[3] (normal uncertain var.)
2400 // Y = xC[4] (lognormal uncertain var.)
2401 Real b = xCM[VAR_b], h = xCM[VAR_h], P = xCM[VAR_P], M = xCM[VAR_M],
2402 Y = xCM[VAR_Y], b_sq = b*b, h_sq = h*h, P_sq = P*P, Y_sq = Y*Y;
2403
2404 // **** f (objective = bh = cross sectional area):
2405 if (numFns > 1 && (directFnASV[ai] & 1))
2406 fnVals[ai] = b*h;
2407
2408 // **** g (limit state = short column response):
2409 if (directFnASV[lsi] & 1)
2410 fnVals[lsi] = 1. - 4.*M/(b*h_sq*Y) - P_sq/(b_sq*h_sq*Y_sq);
2411
2412 // **** df/dx (w.r.t. active variables):
2413 if (numFns > 1 && (directFnASV[ai] & 2))
2414 for (size_t i=0; i<numDerivVars; ++i)
2415 switch (varTypeDVV[i]) {
2416 case VAR_b: // design variable derivative
2417 fnGrads[ai][i] = h; break;
2418 case VAR_h: // design variable derivative
2419 fnGrads[ai][i] = b; break;
2420 default: // uncertain variable derivative
2421 fnGrads[ai][i] = 0.; break;
2422 }
2423
2424 // **** dg/dx (w.r.t. active variables):
2425 if (directFnASV[lsi] & 2)
2426 for (size_t i=0; i<numDerivVars; ++i)
2427 switch (varTypeDVV[i]) {
2428 case VAR_b: // design variable derivative
2429 fnGrads[lsi][i] = 4.*M/(b_sq*h_sq*Y) + 2.*P_sq/(b_sq*b*h_sq*Y_sq);break;
2430 case VAR_h: // design variable derivative
2431 fnGrads[lsi][i] = 8.*M/(b*h_sq*h*Y) + 2.*P_sq/(b_sq*h_sq*h*Y_sq);break;
2432 case VAR_P: // uncertain variable derivative
2433 fnGrads[lsi][i] = -2.*P/(b_sq*h_sq*Y_sq); break;
2434 case VAR_M: // uncertain variable derivative
2435 fnGrads[lsi][i] = -4./(b*h_sq*Y); break;
2436 case VAR_Y: // uncertain variable derivative
2437 fnGrads[lsi][i] = 4.*M/(b*h_sq*Y_sq) + 2.*P_sq/(b_sq*h_sq*Y_sq*Y);break;
2438 }
2439
2440 // **** d^2f/dx^2: (SORM)
2441 if (numFns > 1 && (directFnASV[ai] & 4))
2442 for (size_t i=0; i<numDerivVars; ++i)
2443 for (size_t j=0; j<=i; ++j)
2444 fnHessians[ai](i,j)
2445 = ( (varTypeDVV[i] == VAR_b && varTypeDVV[j] == VAR_h) ||
2446 (varTypeDVV[i] == VAR_h && varTypeDVV[j] == VAR_b) ) ? 1. : 0.;
2447
2448 // **** d^2g/dx^2: (SORM)
2449 if (directFnASV[lsi] & 4)
2450 for (size_t i=0; i<numDerivVars; ++i)
2451 for (size_t j=0; j<=i; ++j)
2452 if (varTypeDVV[i] == VAR_b && varTypeDVV[j] == VAR_b)
2453 fnHessians[lsi](i,j) = -8.*M/(b_sq*b*h_sq*Y)
2454 - 6.*P_sq/(b_sq*b_sq*h_sq*Y_sq);
2455 else if ( (varTypeDVV[i] == VAR_b && varTypeDVV[j] == VAR_h) ||
2456 (varTypeDVV[i] == VAR_h && varTypeDVV[j] == VAR_b) )
2457 fnHessians[lsi](i,j) = -8.*M/(b_sq*h_sq*h*Y)
2458 - 4.*P_sq/(b_sq*b*h_sq*h*Y_sq);
2459 else if (varTypeDVV[i] == VAR_h && varTypeDVV[j] == VAR_h)
2460 fnHessians[lsi](i,j) = -24.*M/(b*h_sq*h_sq*Y)
2461 - 6.*P_sq/(b_sq*h_sq*h_sq*Y_sq);
2462 else if (varTypeDVV[i] == VAR_P && varTypeDVV[j] == VAR_P)
2463 fnHessians[lsi](i,j) = -2./(b_sq*h_sq*Y_sq);
2464 else if ( (varTypeDVV[i] == VAR_P && varTypeDVV[j] == VAR_M) ||
2465 (varTypeDVV[i] == VAR_M && varTypeDVV[j] == VAR_P) )
2466 fnHessians[lsi](i,j) = 0.;
2467 else if ( (varTypeDVV[i] == VAR_P && varTypeDVV[j] == VAR_Y) ||
2468 (varTypeDVV[i] == VAR_Y && varTypeDVV[j] == VAR_P) )
2469 fnHessians[lsi](i,j) = 4.*P/(b_sq*h_sq*Y_sq*Y);
2470 else if (varTypeDVV[i] == VAR_M && varTypeDVV[j] == VAR_M)
2471 fnHessians[lsi](i,j) = 0.;
2472 else if ( (varTypeDVV[i] == VAR_M && varTypeDVV[j] == VAR_Y) ||
2473 (varTypeDVV[i] == VAR_Y && varTypeDVV[j] == VAR_M) )
2474 fnHessians[lsi](i,j) = 4./(b*h_sq*Y_sq);
2475 else if (varTypeDVV[i] == VAR_Y && varTypeDVV[j] == VAR_Y)
2476 fnHessians[lsi](i,j) = -8.*M/(b*h_sq*Y_sq*Y)
2477 - 6.*P_sq/(b_sq*h_sq*Y_sq*Y_sq);
2478 else { // unsupported cross-derivative
2479 Cerr << "Error: unsupported Hessian cross term in short_column."
2480 << std::endl;
2481 abort_handler(INTERFACE_ERROR);
2482 }
2483
2484 return 0; // no failure
2485 }
2486
2487
lf_short_column()2488 int TestDriverInterface::lf_short_column()
2489 {
2490 if (multiProcAnalysisFlag) {
2491 Cerr << "Error: lf_short_column direct fn does not support multiprocessor "
2492 << "analyses." << std::endl;
2493 abort_handler(-1);
2494 }
2495 if (numVars != 5 || numADIV || numADRV) {
2496 Cerr << "Error: Bad number of variables in lf_short_column direct fn."
2497 << std::endl;
2498 abort_handler(INTERFACE_ERROR);
2499 }
2500
2501 short form = 2; // high fidelity case is form 1
2502 if (!analysisComponents.empty() &&
2503 !analysisComponents[analysisDriverIndex].empty()) {
2504 const String& an_comp = analysisComponents[analysisDriverIndex][0];
2505 if (an_comp == "lf1") form = 2;
2506 else if (an_comp == "lf2") form = 3;
2507 else if (an_comp == "lf3") form = 4;
2508 }
2509 return alternate_short_column_forms(form);
2510 }
2511
2512
mf_short_column()2513 int TestDriverInterface::mf_short_column()
2514 {
2515 if (multiProcAnalysisFlag) {
2516 Cerr << "Error: mf_short_column direct fn does not support "
2517 << "multiprocessor analyses." << std::endl;
2518 abort_handler(-1);
2519 }
2520 if (numACV != 5 || numADIV > 1 || numADRV) {
2521 Cerr << "Error: Bad number of variables in mf_short_column direct fn."
2522 << std::endl;
2523 abort_handler(INTERFACE_ERROR);
2524 }
2525 if (numFns > 2) {
2526 Cerr << "Error: Bad number of functions in mf_short_column direct fn."
2527 << std::endl;
2528 abort_handler(INTERFACE_ERROR);
2529 }
2530
2531 int form = xDIM[VAR_MForm];
2532 switch (form) {
2533 case 1: return short_column(); break;
2534 default: return alternate_short_column_forms(form); break;
2535 // monotonic mean: 3 1 2 4
2536 //case 1: return alternate_short_column_forms(3); break;
2537 //case 2: return short_column(); break;
2538 //case 3: return alternate_short_column_forms(2); break;
2539 //case 4: return alternate_short_column_forms(4); break;
2540 // monotonic std dev: 2 1 4 3
2541 //case 1: return alternate_short_column_forms(2); break;
2542 //case 2: return short_column(); break;
2543 //case 3: return alternate_short_column_forms(4); break;
2544 //case 4: return alternate_short_column_forms(3); break;
2545 }
2546 }
2547
2548
alternate_short_column_forms(int form)2549 int TestDriverInterface::alternate_short_column_forms(int form)
2550 {
2551 size_t ai, lsi;
2552 if (numFns == 1) // option for limit state only
2553 lsi = 0;
2554 else if (numFns == 2) // option for area + limit state
2555 { ai = 0; lsi = 1; }
2556 else {
2557 Cerr << "Error: Bad number of functions in alternate_short_column_forms "
2558 << "direct fn." << std::endl;
2559 abort_handler(INTERFACE_ERROR);
2560 }
2561
2562 // b = xC[0] = column base (design var.)
2563 // h = xC[1] = column height (design var.)
2564 // P = xC[2] (normal uncertain var.)
2565 // M = xC[3] (normal uncertain var.)
2566 // Y = xC[4] (lognormal uncertain var.)
2567 Real b = xCM[VAR_b], h = xCM[VAR_h], P = xCM[VAR_P], M = xCM[VAR_M],
2568 Y = xCM[VAR_Y], b_sq = b*b, h_sq = h*h, P_sq = P*P, M_sq = M*M,
2569 Y_sq = Y*Y;
2570
2571 // **** f (objective = bh = cross sectional area):
2572 if (numFns > 1 && (directFnASV[ai] & 1))
2573 fnVals[ai] = b*h;
2574
2575 // **** g (limit state = short column response):
2576 if (directFnASV[lsi] & 1) {
2577 switch (form) {
2578 //case 1: short_column(); // original high fidelity case
2579 case 2:
2580 fnVals[lsi] = 1. - 4.*P/(b*h_sq*Y) - P_sq/(b_sq*h_sq*Y_sq); break;
2581 case 3:
2582 fnVals[lsi] = 1. - 4.*M/(b*h_sq*Y) - M_sq/(b_sq*h_sq*Y_sq); break;
2583 case 4:
2584 fnVals[lsi] = 1. - 4.*M/(b*h_sq*Y) - P_sq/(b_sq*h_sq*Y_sq)
2585 - 4.*(P - M)/(b*h*Y); break;
2586 default: return 1; break;
2587 }
2588 }
2589
2590 return 0; // no failure
2591 }
2592
2593
side_impact_cost()2594 int TestDriverInterface::side_impact_cost()
2595 {
2596 if (numVars != 7 || numFns != 1) {
2597 Cerr << "Error: wrong number of inputs/outputs in side_impact_cost."
2598 << std::endl;
2599 abort_handler(INTERFACE_ERROR);
2600 }
2601
2602 // **** f:
2603 if (directFnASV[0] & 1)
2604 fnVals[0] = 1.98 + 4.9*xC[0] + 6.67*xC[1] + 6.98*xC[2] + 4.01*xC[3]
2605 + 1.78*xC[4] + 2.73*xC[6];
2606
2607 // **** df/dx:
2608 if (directFnASV[0] & 2) {
2609 Real* fn_grad = fnGrads[0];
2610 fn_grad[0] = 4.9; fn_grad[1] = 6.67; fn_grad[2] = 6.98;
2611 fn_grad[3] = 4.01; fn_grad[4] = 1.78; fn_grad[5] = 0.;
2612 fn_grad[6] = 2.73;
2613 }
2614
2615 // **** d^2f/dx^2:
2616 if (directFnASV[0] & 4)
2617 fnHessians[0] = 0.;
2618
2619 return 0; // no failure
2620 }
2621
2622
side_impact_perf()2623 int TestDriverInterface::side_impact_perf()
2624 {
2625 if (numVars != 11 || numFns != 10) {
2626 Cerr << "Error: wrong number of inputs/outputs in side_impact_perf."
2627 << std::endl;
2628 abort_handler(INTERFACE_ERROR);
2629 }
2630
2631 // **** f:
2632 if (directFnASV[0] & 1)
2633 fnVals[0] = 1.16 - 0.3717*xC[1]*xC[3] - 0.00931*xC[1]*xC[9]
2634 - 0.484*xC[2]*xC[8] + 0.01343*xC[5]*xC[9];
2635 if (directFnASV[1] & 1)
2636 fnVals[1] = 28.98 + 3.818*xC[2] - 4.2*xC[0]*xC[1] + 0.0207*xC[4]*xC[9]
2637 + 6.63*xC[5]*xC[8] - 7.7*xC[6]*xC[7] + 0.32*xC[8]*xC[9];
2638 if (directFnASV[2] & 1)
2639 fnVals[2] = 33.86 + 2.95*xC[2] + 0.1792*xC[9] - 5.057*xC[0]*xC[1]
2640 - 11*xC[1]*xC[7] - 0.0215*xC[4]*xC[9] - 9.98*xC[6]*xC[7]
2641 + 22*xC[7]*xC[8];
2642 if (directFnASV[3] & 1)
2643 fnVals[3] = 46.36 - 9.9*xC[1] - 12.9*xC[0]*xC[7] + 0.1107*xC[2]*xC[9];
2644 if (directFnASV[4] & 1)
2645 fnVals[4] = 0.261 - 0.0159*xC[0]*xC[1] - 0.188*xC[0]*xC[7]
2646 - 0.019*xC[1]*xC[6] + 0.0144*xC[2]*xC[4] + 0.0008757*xC[4]*xC[9]
2647 + 0.08045*xC[5]*xC[8] + 0.00139*xC[7]*xC[10]
2648 + 0.00001575*xC[9]*xC[10];
2649 if (directFnASV[5] & 1)
2650 fnVals[5] = 0.214 + 0.00817*xC[4] - 0.131*xC[0]*xC[7] - 0.0704*xC[0]*xC[8]
2651 + 0.03099*xC[1]*xC[5] - 0.018*xC[1]*xC[6] + 0.0208*xC[2]*xC[7]
2652 + 0.121*xC[2]*xC[8] - 0.00364*xC[4]*xC[5] + 0.0007715*xC[4]*xC[9]
2653 - 0.0005354*xC[5]*xC[9] + 0.00121*xC[7]*xC[10];
2654 if (directFnASV[6] & 1)
2655 fnVals[6] = 0.74 - 0.61*xC[1] - 0.163*xC[2]*xC[7] + 0.001232*xC[2]*xC[9]
2656 - 0.166*xC[6]*xC[8] + 0.227*xC[1]*xC[1];
2657 if (directFnASV[7] & 1)
2658 fnVals[7] = 4.72 - 0.5*xC[3] - 0.19*xC[1]*xC[2] - 0.0122*xC[3]*xC[9]
2659 + 0.009325*xC[5]*xC[9] + 0.000191*xC[10]*xC[10];
2660 if (directFnASV[8] & 1)
2661 fnVals[8] = 10.58 - 0.674*xC[0]*xC[1] - 1.95*xC[1]*xC[7]
2662 + 0.02054*xC[2]*xC[9] - 0.0198*xC[3]*xC[9] + 0.028*xC[5]*xC[9];
2663 if (directFnASV[9] & 1)
2664 fnVals[9] = 16.45 - 0.489*xC[2]*xC[6] - 0.843*xC[4]*xC[5]
2665 + 0.0432*xC[8]*xC[9] - 0.0556*xC[8]*xC[10]
2666 - 0.000786*xC[10]*xC[10];
2667
2668 // **** df/dx and d^2f/dx^2:
2669 bool grad_flag = false, hess_flag = false;
2670 for (size_t i=0; i<numFns; ++i) {
2671 if (directFnASV[i] & 2) grad_flag = true;
2672 if (directFnASV[i] & 4) hess_flag = true;
2673 }
2674 if (grad_flag)
2675 Cerr << "Error: gradients not currently supported in side_impact_perf()."
2676 << std::endl;
2677 if (hess_flag)
2678 Cerr << "Error: Hessians not currently supported in side_impact_perf()."
2679 << std::endl;
2680 if (grad_flag || hess_flag)
2681 abort_handler(INTERFACE_ERROR);
2682
2683 return 0; // no failure
2684 }
2685
2686
steel_column_cost()2687 int TestDriverInterface::steel_column_cost()
2688 {
2689 if (numVars != 3 || numFns != 1) {
2690 Cerr << "Error: wrong number of inputs/outputs in steel_column_cost."
2691 << std::endl;
2692 abort_handler(INTERFACE_ERROR);
2693 }
2694
2695 // In the steel column description in Kuschel & Rackwitz, Cost is _not_
2696 // defined as a random variable. That is Cost is not a fn(B, D, H), but
2697 // is rather defined as a fn(b, d, h). Since dCost/dX|_{X=mean} is not the
2698 // same as dCost/dmean for non-normal X (jacobian_dX_dS is not 1), dCost/dX
2699 // may not be used and an optional interface must be defined for Cost.
2700
2701 // b = xC[0] = flange breadth (lognormal unc. var., mean is design var.)
2702 // d = xC[1] = flange thickness (lognormal unc. var., mean is design var.)
2703 // h = xC[2] = profile height (lognormal unc. var., mean is design var.)
2704
2705 Real b = xCM[VAR_b], d = xCM[VAR_d], h = xCM[VAR_h];
2706
2707 // **** f (objective = bd + 5h = cost of column):
2708 if (directFnASV[0] & 1)
2709 fnVals[0] = b*d + 5.*h;
2710
2711 // **** df/dx:
2712 if (directFnASV[0] & 2)
2713 for (size_t i=0; i<numDerivVars; ++i)
2714 switch (varTypeDVV[i]) {
2715 case VAR_b:
2716 fnGrads[0][i] = d; break;
2717 case VAR_d:
2718 fnGrads[0][i] = b; break;
2719 case VAR_h:
2720 fnGrads[0][i] = 5.; break;
2721 }
2722
2723 // **** d^2f/dx^2:
2724 if (directFnASV[0] & 4)
2725 for (size_t i=0; i<numDerivVars; ++i)
2726 for (size_t j=0; j<=i; ++j)
2727 fnHessians[0](i,j)
2728 = ( (varTypeDVV[i] == VAR_b && varTypeDVV[j] == VAR_d) ||
2729 (varTypeDVV[i] == VAR_d && varTypeDVV[j] == VAR_b) ) ? 1. : 0.;
2730
2731 return 0; // no failure
2732 }
2733
2734
steel_column_perf()2735 int TestDriverInterface::steel_column_perf()
2736 {
2737 if (numVars != 9 || numFns != 1) {
2738 Cerr << "Error: wrong number of inputs/outputs in steel_column_perf."
2739 << std::endl;
2740 abort_handler(INTERFACE_ERROR);
2741 }
2742
2743 // In the steel column description in Kuschel & Rackwitz, Cost is _not_
2744 // defined as a random variable. That is Cost is not a fn(B, D, H), but
2745 // is rather defined as a fn(b, d, h). Since dCost/dX|_{X=mean} is not the
2746 // same as dCost/dmean for non-normal X (jacobian_dX_dS is not 1), dCost/dX
2747 // may not be used and an optional interface must be defined for Cost.
2748
2749 // set effective length s based on assumed boundary conditions
2750 // actual length of the column is 7500 mm
2751 Real s = 7500;
2752
2753 // Fs = yield stress (lognormal unc. var.)
2754 // P1 = dead weight load (normal unc. var.)
2755 // P2 = variable load (gumbel unc. var.)
2756 // P3 = variable load (gumbel unc. var.)
2757 // B = flange breadth (lognormal unc. var., mean is design var.)
2758 // D = flange thickness (lognormal unc. var., mean is design var.)
2759 // H = profile height (lognormal unc. var., mean is design var.)
2760 // F0 = init. deflection (normal unc. var.)
2761 // E = elastic modulus (weibull unc. var.)
2762
2763 Real F0 = xCM[VAR_F0], B = xCM[VAR_B], D = xCM[VAR_D], H = xCM[VAR_H],
2764 Fs = xCM[VAR_Fs], E = xCM[VAR_E], P = xCM[VAR_P1]+xCM[VAR_P2]+xCM[VAR_P3],
2765 Pi = 3.1415926535897932385, Pi2 = Pi*Pi, Pi4 = Pi2*Pi2, Pi6 = Pi2*Pi4,
2766 B2 = B*B, D2 = D*D, H2 = H*H, H3 = H*H2, H5 = H2*H3, E2 = E*E, E3 = E*E2,
2767 s2 = s*s, X = Pi2*E*B*D*H2 - 2.*s2*P, X2 = X*X, X3 = X*X2;
2768
2769 // **** g (limit state):
2770 if (directFnASV[0] & 1)
2771 fnVals[0] = Fs - P*(1./2./B/D + Pi2*F0*E*H/X);
2772
2773 // **** dg/dx (w.r.t. active/uncertain variables):
2774 if (directFnASV[0] & 2)
2775 for (size_t i=0; i<numDerivVars; ++i)
2776 switch (varTypeDVV[i]) {
2777 case VAR_F0: // df/dF0
2778 fnGrads[0][i] = -E*H*P*Pi2/X; break;
2779 case VAR_P1: case VAR_P2: case VAR_P3: // df/dP1, df/dP2, df/dP3
2780 fnGrads[0][i] = -1./2./B/D - B*D*E2*F0*H3*Pi4/X2; break;
2781 case VAR_B: // df/dB
2782 fnGrads[0][i] = P*(1./2./B2/D + D*E2*F0*H3*Pi4/X2); break;
2783 case VAR_D: // df/dD
2784 fnGrads[0][i] = P*(1./2./B/D2 + B*E2*F0*H3*Pi4/X2); break;
2785 case VAR_H: // df/dH
2786 fnGrads[0][i] = E*F0*P*Pi2*(X + 4.*P*s2)/X2; break;
2787 case VAR_Fs: // df/dFs
2788 fnGrads[0][i] = 1.; break;
2789 case VAR_E: // df/dE
2790 fnGrads[0][i] = 2.*F0*H*P*P*Pi2*s2/X2; break;
2791 }
2792
2793 // **** d^2g/dx^2: (SORM)
2794 if (directFnASV[0] & 4)
2795 for (size_t i=0; i<numDerivVars; ++i)
2796 for (size_t j=0; j<=i; ++j)
2797 if (varTypeDVV[i] == VAR_Fs || varTypeDVV[j] == VAR_Fs)
2798 fnHessians[0](i,j) = 0.;
2799 else if ( (varTypeDVV[i] == VAR_P1 && varTypeDVV[j] == VAR_P1) ||
2800 (varTypeDVV[i] == VAR_P1 && varTypeDVV[j] == VAR_P2) ||
2801 (varTypeDVV[i] == VAR_P1 && varTypeDVV[j] == VAR_P3) ||
2802 (varTypeDVV[i] == VAR_P2 && varTypeDVV[j] == VAR_P1) ||
2803 (varTypeDVV[i] == VAR_P2 && varTypeDVV[j] == VAR_P2) ||
2804 (varTypeDVV[i] == VAR_P2 && varTypeDVV[j] == VAR_P3) ||
2805 (varTypeDVV[i] == VAR_P3 && varTypeDVV[j] == VAR_P1) ||
2806 (varTypeDVV[i] == VAR_P3 && varTypeDVV[j] == VAR_P2) ||
2807 (varTypeDVV[i] == VAR_P3 && varTypeDVV[j] == VAR_P3) )
2808 fnHessians[0](i,j) = -4.*B*D*E2*F0*H3*Pi4*s2/X3;
2809 else if ( (varTypeDVV[i] == VAR_P1 && varTypeDVV[j] == VAR_B) ||
2810 (varTypeDVV[i] == VAR_P2 && varTypeDVV[j] == VAR_B) ||
2811 (varTypeDVV[i] == VAR_P3 && varTypeDVV[j] == VAR_B) ||
2812 (varTypeDVV[i] == VAR_B && varTypeDVV[j] == VAR_P1) ||
2813 (varTypeDVV[i] == VAR_B && varTypeDVV[j] == VAR_P2) ||
2814 (varTypeDVV[i] == VAR_B && varTypeDVV[j] == VAR_P3) )
2815 fnHessians[0](i,j)
2816 = 1./2./B2/D + D*E2*F0*H3*Pi4/X2*(2.*B*D*E*H2*Pi2/X - 1.);
2817 else if ( (varTypeDVV[i] == VAR_P1 && varTypeDVV[j] == VAR_D) ||
2818 (varTypeDVV[i] == VAR_P2 && varTypeDVV[j] == VAR_D) ||
2819 (varTypeDVV[i] == VAR_P3 && varTypeDVV[j] == VAR_D) ||
2820 (varTypeDVV[i] == VAR_D && varTypeDVV[j] == VAR_P1) ||
2821 (varTypeDVV[i] == VAR_D && varTypeDVV[j] == VAR_P2) ||
2822 (varTypeDVV[i] == VAR_D && varTypeDVV[j] == VAR_P3) )
2823 fnHessians[0](i,j)
2824 = 1./2./B/D2 + B*E2*F0*H3*Pi4/X2*(2.*B*D*E*H2*Pi2/X - 1.);
2825 else if ( (varTypeDVV[i] == VAR_P1 && varTypeDVV[j] == VAR_H) ||
2826 (varTypeDVV[i] == VAR_P2 && varTypeDVV[j] == VAR_H) ||
2827 (varTypeDVV[i] == VAR_P3 && varTypeDVV[j] == VAR_H) ||
2828 (varTypeDVV[i] == VAR_H && varTypeDVV[j] == VAR_P1) ||
2829 (varTypeDVV[i] == VAR_H && varTypeDVV[j] == VAR_P2) ||
2830 (varTypeDVV[i] == VAR_H && varTypeDVV[j] == VAR_P3) )
2831 fnHessians[0](i,j) = B*D*E2*F0*H2*Pi4*(X+8.*P*s2)/X3;
2832 else if ( (varTypeDVV[i] == VAR_P1 && varTypeDVV[j] == VAR_F0) ||
2833 (varTypeDVV[i] == VAR_P2 && varTypeDVV[j] == VAR_F0) ||
2834 (varTypeDVV[i] == VAR_P3 && varTypeDVV[j] == VAR_F0) ||
2835 (varTypeDVV[i] == VAR_F0 && varTypeDVV[j] == VAR_P1) ||
2836 (varTypeDVV[i] == VAR_F0 && varTypeDVV[j] == VAR_P2) ||
2837 (varTypeDVV[i] == VAR_F0 && varTypeDVV[j] == VAR_P3) )
2838 fnHessians[0](i,j) = -B*D*E2*H3*Pi4/X2;
2839 else if ( (varTypeDVV[i] == VAR_P1 && varTypeDVV[j] == VAR_E) ||
2840 (varTypeDVV[i] == VAR_P2 && varTypeDVV[j] == VAR_E) ||
2841 (varTypeDVV[i] == VAR_P3 && varTypeDVV[j] == VAR_E) ||
2842 (varTypeDVV[i] == VAR_E && varTypeDVV[j] == VAR_P1) ||
2843 (varTypeDVV[i] == VAR_E && varTypeDVV[j] == VAR_P2) ||
2844 (varTypeDVV[i] == VAR_E && varTypeDVV[j] == VAR_P3) )
2845 fnHessians[0](i,j) = 4.*B*D*E*F0*H3*P*Pi4*s2/X3;
2846 else if (varTypeDVV[i] == VAR_B && varTypeDVV[j] == VAR_B)
2847 fnHessians[0](i,j) = -P*(1./B2/B/D + 2.*D2*E3*F0*H5*Pi6/X3);
2848 else if ( (varTypeDVV[i] == VAR_B && varTypeDVV[j] == VAR_D) ||
2849 (varTypeDVV[i] == VAR_D && varTypeDVV[j] == VAR_B) )
2850 fnHessians[0](i,j)
2851 = -P*(1./2./B2/D2 + E2*F0*H3*Pi4/X2*(2.*B*D*E*H2*Pi2/X - 1.));
2852 else if ( (varTypeDVV[i] == VAR_B && varTypeDVV[j] == VAR_H) ||
2853 (varTypeDVV[i] == VAR_H && varTypeDVV[j] == VAR_B) )
2854 fnHessians[0](i,j) = -D*E2*F0*H2*P*Pi4*(X + 8.*P*s2)/X3;
2855 else if ( (varTypeDVV[i] == VAR_F0 && varTypeDVV[j] == VAR_B) ||
2856 (varTypeDVV[i] == VAR_B && varTypeDVV[j] == VAR_F0) )
2857 fnHessians[0](i,j) = D*E2*H3*P*Pi4/X2;
2858 else if ( (varTypeDVV[i] == VAR_B && varTypeDVV[j] == VAR_E) ||
2859 (varTypeDVV[i] == VAR_E && varTypeDVV[j] == VAR_B) )
2860 fnHessians[0](i,j) = -4.*D*E*F0*H3*P*P*Pi4*s2/X3;
2861 else if (varTypeDVV[i] == VAR_D && varTypeDVV[j] == VAR_D)
2862 fnHessians[0](i,j) = -P*(1./B/D2/D + 2.*B2*E3*F0*H5*Pi6/X3);
2863 else if ( (varTypeDVV[i] == VAR_D && varTypeDVV[j] == VAR_H) ||
2864 (varTypeDVV[i] == VAR_H && varTypeDVV[j] == VAR_D) )
2865 fnHessians[0](i,j) = -B*E2*F0*H2*P*Pi4*(X + 8.*P*s2)/X3;
2866 else if ( (varTypeDVV[i] == VAR_F0 && varTypeDVV[j] == VAR_D) ||
2867 (varTypeDVV[i] == VAR_D && varTypeDVV[j] == VAR_F0) )
2868 fnHessians[0](i,j) = B*E2*H3*P*Pi4/X2;
2869 else if ( (varTypeDVV[i] == VAR_D && varTypeDVV[j] == VAR_E) ||
2870 (varTypeDVV[i] == VAR_E && varTypeDVV[j] == VAR_D) )
2871 fnHessians[0](i,j) = -4.*B*E*F0*H3*P*P*Pi4*s2/X3;
2872 else if (varTypeDVV[i] == VAR_H && varTypeDVV[j] == VAR_H)
2873 fnHessians[0](i,j) = -2.*B*D*E2*F0*H*P*Pi4*(X + 8.*P*s2)/X3;
2874 else if ( (varTypeDVV[i] == VAR_F0 && varTypeDVV[j] == VAR_H) ||
2875 (varTypeDVV[i] == VAR_H && varTypeDVV[j] == VAR_F0) )
2876 fnHessians[0](i,j) = E*P*Pi2*(X + 4.*P*s2)/X2;
2877 else if ( (varTypeDVV[i] == VAR_H && varTypeDVV[j] == VAR_E) ||
2878 (varTypeDVV[i] == VAR_E && varTypeDVV[j] == VAR_H) )
2879 fnHessians[0](i,j) = -2.*F0*P*P*Pi2*s2*(3.*X + 8.*P*s2)/X3;
2880 else if (varTypeDVV[i] == VAR_F0 && varTypeDVV[j] == VAR_F0)
2881 fnHessians[0](i,j) = 0.;
2882 else if ( (varTypeDVV[i] == VAR_F0 && varTypeDVV[j] == VAR_E) ||
2883 (varTypeDVV[i] == VAR_E && varTypeDVV[j] == VAR_F0) )
2884 fnHessians[0](i,j) = 2.*H*P*P*Pi2*s2/X2;
2885 else if (varTypeDVV[i] == VAR_E && varTypeDVV[j] == VAR_E)
2886 fnHessians[0](i,j) = -4.*B*D*F0*H3*P*P*Pi4*s2/X3;
2887 else { // unsupported derivative
2888 Cerr << "Error: unsupported Hessian cross term in steel_column."
2889 << std::endl;
2890 abort_handler(INTERFACE_ERROR);
2891 }
2892
2893 return 0; // no failure
2894 }
2895
2896
sobol_rational()2897 int TestDriverInterface::sobol_rational()
2898 {
2899 if (multiProcAnalysisFlag) {
2900 Cerr << "Error: sobol_rational direct fn does not support multiprocessor "
2901 << "analyses." << std::endl;
2902 abort_handler(-1);
2903 }
2904 if (numVars != 2 || numFns != 1) {
2905 Cerr << "Error: Bad number of inputs/outputs in sobol_rational direct fn."
2906 << std::endl;
2907 abort_handler(INTERFACE_ERROR);
2908 }
2909
2910 // f = (x2 + 0.5)^4 / (x1 + 0.5)^2
2911 // See Storlie et al. SAND2008-6570
2912
2913 const Real& x1 = xC[0]; const Real& x2 = xC[1];
2914
2915 // **** f:
2916 if (directFnASV[0] & 1)
2917 fnVals[0] = std::pow((x2 + 0.5), 4.) / std::pow((x1 + 0.5), 2.);
2918
2919 // **** df/dx:
2920 if (directFnASV[0] & 2)
2921 for (size_t i=0; i<numDerivVars; ++i)
2922 switch (directFnDVV[i]) {
2923 case 1: // x1
2924 fnGrads[0][i] = -2.*std::pow(x2 + 0.5, 4.)/std::pow(x1 + 0.5, 3.);
2925 break;
2926 case 2: // x2
2927 fnGrads[0][i] = 4.*std::pow(x2 + 0.5, 3.)/std::pow(x1 + 0.5, 2.);
2928 break;
2929 }
2930
2931 return 0; // no failure
2932 }
2933
2934
sobol_g_function()2935 int TestDriverInterface::sobol_g_function()
2936 {
2937 if (multiProcAnalysisFlag) {
2938 Cerr << "Error: sobol_g_function direct fn does not support multiprocessor "
2939 << "analyses." << std::endl;
2940 abort_handler(-1);
2941 }
2942 if (numVars < 1 || numVars > 10 || numFns != 1) {
2943 Cerr << "Error: Bad number of inputs/outputs in sobol_g_function direct fn."
2944 << std::endl;
2945 abort_handler(INTERFACE_ERROR);
2946 }
2947
2948 // Sobol g-Function: see Storlie et al. SAND2008-6570
2949 int a[] = {0,1,2,4,8,99,99,99,99,99};
2950
2951 // **** f:
2952 if (directFnASV[0] & 1) {
2953 fnVals[0] = 2.;
2954 for (int i=0; i<numVars; ++i)
2955 fnVals[0] *= ( std::abs(4.*xC[i] - 2.) + a[i] ) / ( 1. + a[i] );
2956 }
2957
2958 // **** df/dx:
2959 if (directFnASV[0] & 2) {
2960 for (size_t i=0; i<numDerivVars; ++i) {
2961 size_t var_index = directFnDVV[i] - 1;
2962 Real& fn_grad_0i = fnGrads[0][i];
2963 if (4.*xC[var_index] == 2.) // zero gradient assumed at discontinuity
2964 fn_grad_0i = 0.;
2965 else {
2966 fn_grad_0i = (4.*xC[var_index] > 2.) ?
2967 8. / ( 1. + a[var_index] ) : -8. / ( 1. + a[var_index] );
2968 for (int j=0; j<numVars; ++j)
2969 if (j != var_index)
2970 fn_grad_0i *= ( std::abs(4.*xC[j] - 2.) + a[j] ) / ( 1. + a[j] );
2971 }
2972 }
2973 }
2974
2975 return 0; // no failure
2976 }
2977
2978
sobol_ishigami()2979 int TestDriverInterface::sobol_ishigami()
2980 {
2981 using std::pow;
2982 using std::sin;
2983 using std::cos;
2984
2985 if (multiProcAnalysisFlag) {
2986 Cerr << "Error: sobol_ishigami direct fn does not support multiprocessor "
2987 << "analyses." << std::endl;
2988 abort_handler(-1);
2989 }
2990 if (numVars != 3 || numFns != 1) {
2991 Cerr << "Error: Bad number of inputs/outputs in sobol_ishigami direct fn."
2992 << std::endl;
2993 abort_handler(INTERFACE_ERROR);
2994 }
2995
2996 // Ishigami Function: see Storlie et al. SAND2008-6570
2997 const Real pi = 3.14159265358979324;
2998 Real x1 = xCM[VAR_x1], x2 = xCM[VAR_x2], x3 = xCM[VAR_x3];
2999
3000 // **** f:
3001 if (directFnASV[0] & 1)
3002 fnVals[0] = ( 1. + 0.1 * pow(2.*pi*x3 - pi, 4.0) ) *
3003 sin(2.*pi*x1 - pi) + 7. * pow(sin(2*pi*x2 - pi), 2.0);
3004
3005 // **** df/dx:
3006 if (directFnASV[0] & 2)
3007 for (size_t i=0; i<numDerivVars; ++i)
3008 switch (varTypeDVV[i]) {
3009 case VAR_x1:
3010 fnGrads[0][i]
3011 = 2.*pi * ( 1. + 0.1*pow(2.*pi*x3 - pi, 4.) ) * cos(2.*pi*x1 - pi);
3012 break;
3013 case VAR_x2:
3014 fnGrads[0][i] = 28.*pi * sin(2.*pi*x2 - pi) * cos(2.*pi*x2 - pi);
3015 break;
3016 case VAR_x3:
3017 fnGrads[0][i] = 0.8 * pow(2.*pi*x3 - pi, 3.) * sin(2.*pi*x1 - pi);
3018 break;
3019 }
3020
3021 return 0; // no failure
3022 }
3023
3024
text_book()3025 int TestDriverInterface::text_book()
3026 {
3027 // This version is used when multiple analysis drivers are not employed.
3028 // In this case, execute text_book1/2/3 sequentially.
3029
3030 if (numFns > 3) {
3031 Cerr << "Error: Bad number of functions in text_book direct fn."
3032 << std::endl;
3033 abort_handler(INTERFACE_ERROR);
3034 }
3035 // The presence of discrete variables can cause offsets in directFnDVV which
3036 // the text_book derivative logic does not currently account for.
3037 if ( (gradFlag || hessFlag) && (numADIV || numADRV || numADSV) ) {
3038 Cerr << "Error: text_book direct fn assumes no discrete variables in "
3039 << "derivative mode." << std::endl;
3040 abort_handler(INTERFACE_ERROR);
3041 }
3042
3043 text_book1(); // objective fn val/grad/Hessian
3044 if (numFns > 1)
3045 text_book2(); // constraint 1 val/grad/Hessian
3046 if (numFns > 2)
3047 text_book3(); // constraint 2 val/grad/Hessian
3048
3049 // Test failure capturing for Direct case
3050 //int r = rand();
3051 //if (r < 10000) // RAND_MAX = 32767
3052 // return 1; // failure
3053
3054 // for faking a more expensive evaluation:
3055 //std::this_thread::sleep_for(std::chrono::seconds(5));
3056 return 0; // no failure
3057 }
3058
3059
3060 // text_book1/2/3 are used when evalComm is split into multiple analysis
3061 // servers. In this case, the 3 portions are executed in parallel.
text_book1()3062 int TestDriverInterface::text_book1()
3063 {
3064 // **********************************
3065 // **** f: sum (x[i] - POWVAL)^4 ****
3066 // **********************************
3067 size_t i;
3068 if (directFnASV[0] & 1) {
3069 Real local_val = 0.0;
3070 for (i=analysisCommRank; i<numVars; i+=analysisCommSize) {
3071 // orders all continuous vars followed by all discrete vars. This is
3072 // fine in the direct case so long as everything is self-consistent.
3073 Real x_i;
3074 if (i<numACV)
3075 x_i = xC[i];
3076 else if (i<numACV+numADIV)
3077 x_i = (Real)xDI[i-numACV];
3078 else if (i<numACV+numADIV+numADRV)
3079 x_i = xDR[i-numACV-numADIV];
3080 else
3081 x_i = levenshtein_distance(xDS[i-numACV-numADIV-numADRV]);
3082 local_val += std::pow(x_i-POW_VAL, 4);
3083 #ifdef TB_EXPENSIVE
3084 // MDO98/WCSMO99 TFLOP/NOW timings: j<=15000 used for fnVals[0] only
3085 // (text_book1 only)
3086 // StrOpt_2002 TFLOP timings: j<=5000 used for all fnVals, fnGrads, and
3087 // fnHessians that are present (text_book1/2/3)
3088 for (size_t j=1; j<=5000; ++j)
3089 local_val += 1./(std::pow(x_i-POW_VAL,4)+j/100.)
3090 /(std::pow(x_i-POW_VAL,4)+j/100.);
3091 #endif // TB_EXPENSIVE
3092 }
3093
3094 if (multiProcAnalysisFlag) {
3095 Real global_val = 0.0;
3096 parallelLib.reduce_sum_a(&local_val, &global_val, 1);
3097 // only analysisCommRank 0 has the correct sum. This is OK (MPI_Allreduce
3098 // not needed) since only analysisCommRank 0 updates response for
3099 // evalCommRank 0 in overlay_response. evalCommRank 0 then returns the
3100 // results to the iterator in ApplicationInterface::serve_evaluations().
3101 if (analysisCommRank == 0)
3102 fnVals[0] = global_val;
3103 }
3104 else
3105 fnVals[0] = local_val;
3106 }
3107
3108 // ****************
3109 // **** df/dx: ****
3110 // ****************
3111 if (directFnASV[0] & 2) {
3112 //for (i=0; i<numDerivVars; ++i)
3113 // fnGrads[0][i] = 4.*pow(xC[i]-POW_VAL,3);
3114 std::fill_n(fnGrads[0], fnGrads.numRows(), 0.);
3115 for (i=analysisCommRank; i<numDerivVars; i+=analysisCommSize) {
3116 size_t var_index = directFnDVV[i] - 1; // assumes no discrete vars
3117 Real x_i = xC[var_index]; // assumes no discrete vars
3118 fnGrads[0][i] = 4.*std::pow(x_i-POW_VAL,3);
3119 #ifdef TB_EXPENSIVE
3120 for (size_t j=1; j<=5000; ++j)
3121 fnGrads[0][i] += 1./(std::pow(x_i-POW_VAL,3)+j/100.)
3122 /(std::pow(x_i-POW_VAL,3)+j/100.);
3123 #endif // TB_EXPENSIVE
3124 }
3125
3126 if (multiProcAnalysisFlag) {
3127 Real* sum_fns = (analysisCommRank) ? NULL : new Real [numDerivVars];
3128 parallelLib.reduce_sum_a((Real*)fnGrads[0], sum_fns,
3129 numDerivVars);
3130 if (analysisCommRank == 0) {
3131 RealVector fn_grad_col_vec = Teuchos::getCol(Teuchos::View, fnGrads, 0);
3132 copy_data(sum_fns, (int)numDerivVars, fn_grad_col_vec);
3133 delete [] sum_fns;
3134 }
3135 }
3136 }
3137
3138 // ********************
3139 // **** d^2f/dx^2: ****
3140 // ********************
3141 if (directFnASV[0] & 4) {
3142 fnHessians[0] = 0.;
3143 //for (i=0; i<numDerivVars; ++i)
3144 // fnHessians[0][i][i] = 12.*pow(xC[i]-POW_VAL,2);
3145 for (i=analysisCommRank; i<numDerivVars; i+=analysisCommSize) {
3146 size_t var_index = directFnDVV[i] - 1; // assumes no discrete vars
3147 Real x_i = xC[var_index]; // assumes no discrete vars
3148 fnHessians[0](i,i) = 12.*std::pow(x_i-POW_VAL,2);
3149 #ifdef TB_EXPENSIVE
3150 for (size_t j=0; j<numDerivVars; ++j)
3151 for (size_t k=1; k<=5000; ++k)
3152 fnHessians[0](i,j) += 1./(std::pow(x_i-POW_VAL,2)+k/100.)
3153 /(std::pow(x_i-POW_VAL,2)+k/100.);
3154 #endif // TB_EXPENSIVE
3155 }
3156
3157 if (multiProcAnalysisFlag) {
3158 int num_reals = numDerivVars * numDerivVars;
3159 Real* local_fns = new Real [num_reals];
3160 std::copy(fnHessians[0].values(), fnHessians[0].values() + num_reals,
3161 local_fns);
3162 Real* sum_fns = (analysisCommRank) ? NULL : new Real [num_reals];
3163 parallelLib.reduce_sum_a(local_fns, sum_fns, num_reals);
3164 delete [] local_fns;
3165 if (analysisCommRank == 0) {
3166 std::copy(sum_fns, sum_fns + num_reals, fnHessians[0].values());
3167 delete [] sum_fns;
3168 }
3169 }
3170 }
3171
3172 // for faking a more expensive evaluation:
3173 //std::this_thread::sleep_for(std::chrono::seconds(1));
3174
3175 return 0;
3176 }
3177
3178
text_book2()3179 int TestDriverInterface::text_book2()
3180 {
3181 // **********************************
3182 // **** c1: x[0]*x[0] - 0.5*x[1] ****
3183 // **********************************
3184 size_t i;
3185 if (directFnASV[1] & 1) {
3186 Real local_val = 0.0;
3187 // Definitely not the most efficient way to do this, but the point is to
3188 // demonstrate Comm communication.
3189 for (i=analysisCommRank; i<numVars; i+=analysisCommSize) {
3190 // orders all continuous vars followed by all discrete vars. This is
3191 // fine in the direct case so long as everything is self-consistent.
3192 Real x_i;
3193 if (i<numACV)
3194 x_i = xC[i];
3195 else if (i<numACV+numADIV)
3196 x_i = (Real)xDI[i-numACV];
3197 else if (i<numACV+numADIV+numADRV)
3198 x_i = xDR[i-numACV-numADIV];
3199 else
3200 x_i = levenshtein_distance(xDS[i-numACV-numADIV-numADRV]);
3201 if (i==0) // could be changed to i % 2 == 0 to get even vars.
3202 local_val += x_i*x_i;
3203 else if (i==1) // could be changed to i % 2 == 1 to get odd vars
3204 local_val -= 0.5*x_i;
3205 #ifdef TB_EXPENSIVE
3206 for (size_t j=1; j<=5000; ++j)
3207 local_val += 1./(std::pow(x_i-POW_VAL,4)+j/100.)
3208 /(std::pow(x_i-POW_VAL,4)+j/100.);
3209 #endif // TB_EXPENSIVE
3210 }
3211
3212 if (multiProcAnalysisFlag) {
3213 Real global_val = 0.0;
3214 parallelLib.reduce_sum_a(&local_val, &global_val, 1);
3215 // only analysisCommRank 0 has the correct sum. This is OK (MPI_Allreduce
3216 // not needed) since only analysisCommRank 0 updates response for
3217 // evalCommRank 0 in overlay_response. evalCommRank 0 then returns the
3218 // results to the iterator in ApplicationInterface::serve_evaluations().
3219 if (analysisCommRank == 0)
3220 fnVals[1] = global_val;
3221 }
3222 else
3223 fnVals[1] = local_val;
3224 }
3225
3226 // *****************
3227 // **** dc1/dx: ****
3228 // *****************
3229 if (directFnASV[1] & 2) {
3230 std::fill_n(fnGrads[1], fnGrads.numRows(), 0.);
3231
3232 //fnGrads[1][0] = 2.*xC[0];
3233 //fnGrads[1][1] = -0.5;
3234 for (i=analysisCommRank; i<numDerivVars; i+=analysisCommSize) {
3235 size_t var_index = directFnDVV[i] - 1; // assumes no discrete vars
3236 if (var_index == 0)
3237 fnGrads[1][i] = 2.*xC[0];
3238 else if (var_index == 1)
3239 fnGrads[1][i] = -0.5;
3240 #ifdef TB_EXPENSIVE
3241 Real x_i = xC[var_index];
3242 for (size_t j=1; j<=5000; ++j)
3243 fnGrads[1][i] += 1./(std::pow(x_i-POW_VAL,3)+j/100.)
3244 /(std::pow(x_i-POW_VAL,3)+j/100.);
3245 #endif // TB_EXPENSIVE
3246 }
3247
3248 if (multiProcAnalysisFlag) {
3249 Real* sum_fns = (analysisCommRank) ? NULL : new Real [numDerivVars];
3250 parallelLib.reduce_sum_a((Real*)fnGrads[1], sum_fns,
3251 numDerivVars);
3252 if (analysisCommRank == 0) {
3253 RealVector fn_grad_col_vec = Teuchos::getCol(Teuchos::View, fnGrads, 1);
3254 copy_data(sum_fns, (int)numDerivVars, fn_grad_col_vec);
3255 delete [] sum_fns;
3256 }
3257 }
3258 }
3259
3260 // *********************
3261 // **** d^2c1/dx^2: ****
3262 // *********************
3263 if (directFnASV[1] & 4) {
3264 fnHessians[1] = 0.;
3265 //fnHessians[1][0][0] = 2.;
3266 for (i=analysisCommRank; i<numDerivVars; i+=analysisCommSize) {
3267 size_t var_index = directFnDVV[i] - 1; // assumes no discrete vars
3268 if (var_index == 0)
3269 fnHessians[1](i,i) = 2.;
3270 #ifdef TB_EXPENSIVE
3271 Real x_i = xC[var_index];
3272 for (size_t j=0; j<numDerivVars; ++j)
3273 for (size_t k=1; k<=5000; ++k)
3274 fnHessians[1](i,j) += 1./(std::pow(x_i-POW_VAL,2)+k/100.)
3275 /(std::pow(x_i-POW_VAL,2)+k/100.);
3276 #endif // TB_EXPENSIVE
3277 }
3278
3279 if (multiProcAnalysisFlag) {
3280 int num_reals = numDerivVars * numDerivVars;
3281 Real* local_fns = new Real [num_reals];
3282 std::copy(fnHessians[1].values(), fnHessians[1].values() + num_reals,
3283 local_fns);
3284 Real* sum_fns = (analysisCommRank) ? NULL : new Real [num_reals];
3285 parallelLib.reduce_sum_a(local_fns, sum_fns, num_reals);
3286 delete [] local_fns;
3287 if (analysisCommRank == 0) {
3288 std::copy( sum_fns, sum_fns + num_reals, fnHessians[1].values() );
3289 delete [] sum_fns;
3290 }
3291 }
3292 }
3293
3294 // for faking a more expensive evaluation:
3295 //std::this_thread::sleep_for(std::chrono::seconds(1));
3296
3297 return 0;
3298 }
3299
3300
text_book3()3301 int TestDriverInterface::text_book3()
3302 {
3303 // **********************************
3304 // **** c2: x[1]*x[1] - 0.5*x[0] ****
3305 // **********************************
3306 size_t i;
3307 if (directFnASV[2] & 1) {
3308 Real local_val = 0.0;
3309 // Definitely not the most efficient way to do this, but the point is to
3310 // demonstrate Comm communication.
3311 for (i=analysisCommRank; i<numVars; i+=analysisCommSize) {
3312 // orders all continuous vars followed by all discrete vars. This is
3313 // fine in the direct case so long as everything is self-consistent.
3314 Real x_i;
3315 if (i<numACV)
3316 x_i = xC[i];
3317 else if (i<numACV+numADIV)
3318 x_i = (Real)xDI[i-numACV];
3319 else if (i<numACV+numADIV+numADRV)
3320 x_i = xDR[i-numACV-numADIV];
3321 else
3322 x_i = levenshtein_distance(xDS[i-numACV-numADIV-numADRV]);
3323 if (i==0) // could be changed to i % 2 == 0 to get even vars.
3324 local_val -= 0.5*x_i;
3325 else if (i==1) // could be changed to i % 2 == 1 to get odd vars
3326 local_val += x_i*x_i;
3327 #ifdef TB_EXPENSIVE
3328 for (size_t j=1; j<=5000; ++j)
3329 local_val += 1./(std::pow(x_i-POW_VAL,4)+j/100.)
3330 /(std::pow(x_i-POW_VAL,4)+j/100.);
3331 #endif // TB_EXPENSIVE
3332 }
3333
3334 if (multiProcAnalysisFlag) {
3335 Real global_val = 0.0;
3336 parallelLib.reduce_sum_a(&local_val, &global_val, 1);
3337 // only analysisCommRank 0 has the correct sum. This is OK (MPI_Allreduce
3338 // not needed) since only analysisCommRank 0 updates response for
3339 // evalCommRank 0 in overlay_response. evalCommRank 0 then returns the
3340 // results to the iterator in ApplicationInterface::serve_evaluations().
3341 if (analysisCommRank == 0)
3342 fnVals[2] = global_val;
3343 }
3344 else
3345 fnVals[2] = local_val;
3346 }
3347
3348 // *****************
3349 // **** dc2/dx: ****
3350 // *****************
3351 if (directFnASV[2] & 2) {
3352 std::fill_n(fnGrads[2], fnGrads.numRows(), 0.);
3353
3354 //fnGrads[2][0] = -0.5;
3355 //fnGrads[2][1] = 2.*xC[1];
3356 for (i=analysisCommRank; i<numDerivVars; i+=analysisCommSize) {
3357 size_t var_index = directFnDVV[i] - 1; // assumes no discrete vars
3358 if (var_index == 0)
3359 fnGrads[2][i] = -0.5;
3360 else if (var_index == 1)
3361 fnGrads[2][i] = 2.*xC[1];
3362 #ifdef TB_EXPENSIVE
3363 Real x_i = xC[var_index];
3364 for (size_t j=1; j<=5000; ++j)
3365 fnGrads[2][i] += 1./(std::pow(x_i-POW_VAL,3)+j/100.)
3366 /(std::pow(x_i-POW_VAL,3)+j/100.);
3367 #endif // TB_EXPENSIVE
3368 }
3369
3370 if (multiProcAnalysisFlag) {
3371 Real* sum_fns = (analysisCommRank) ? NULL : new Real [numDerivVars];
3372 parallelLib.reduce_sum_a((Real*)fnGrads[2], sum_fns,
3373 numDerivVars);
3374 if (analysisCommRank == 0) {
3375 RealVector fn_grad_col_vec = Teuchos::getCol(Teuchos::View, fnGrads, 2);
3376 copy_data(sum_fns, (int)numDerivVars, fn_grad_col_vec);
3377 delete [] sum_fns;
3378 }
3379 }
3380 }
3381
3382 // *********************
3383 // **** d^2c2/dx^2: ****
3384 // *********************
3385 if (directFnASV[2] & 4) {
3386 fnHessians[2] = 0.;
3387 //fnHessians[2][1][1] = 2.;
3388 for (i=analysisCommRank; i<numDerivVars; i+=analysisCommSize) {
3389 size_t var_index = directFnDVV[i] - 1; // assumes no discrete vars
3390 if (var_index == 1)
3391 fnHessians[2](i,i) = 2.;
3392 #ifdef TB_EXPENSIVE
3393 Real x_i = xC[var_index];
3394 for (size_t j=0; j<numDerivVars; ++j)
3395 for (size_t k=1; k<=5000; ++k)
3396 fnHessians[2](i,j) += 1./(std::pow(x_i-POW_VAL,2)+k/100.)
3397 /(std::pow(x_i-POW_VAL,2)+k/100.);
3398 #endif // TB_EXPENSIVE
3399 }
3400
3401 if (multiProcAnalysisFlag) {
3402 int num_reals = numDerivVars * numDerivVars;
3403 Real* local_fns = new Real [num_reals];
3404 std::copy(fnHessians[2].values(), fnHessians[2].values() + num_reals,
3405 local_fns);
3406 Real* sum_fns = (analysisCommRank) ? NULL : new Real [num_reals];
3407 parallelLib.reduce_sum_a(local_fns, sum_fns, num_reals);
3408 delete [] local_fns;
3409 if (analysisCommRank == 0) {
3410 std::copy(sum_fns, sum_fns + num_reals, fnHessians[2].values());
3411 delete [] sum_fns;
3412 }
3413 }
3414 }
3415
3416 // for faking a more expensive evaluation:
3417 //std::this_thread::sleep_for(std::chrono::seconds(1));
3418
3419 return 0;
3420 }
3421
3422
text_book_ouu()3423 int TestDriverInterface::text_book_ouu()
3424 {
3425 if (multiProcAnalysisFlag) {
3426 Cerr << "Error: text_book_ouu direct fn does not support multiprocessor "
3427 << "analyses." << std::endl;
3428 abort_handler(-1);
3429 }
3430 // typical usage is 2 design vars + 6 uncertain variables, although the
3431 // number of uncertain variables can be any factor of two.
3432 if (numVars < 4 || numVars % 2 || numADIV || numADRV) {
3433 Cerr << "Error: Bad number of variables in text_book_ouu direct fn."
3434 << std::endl;
3435 abort_handler(INTERFACE_ERROR);
3436 }
3437 if (numFns > 3) {
3438 Cerr << "Error: Bad number of functions in text_book_ouu direct fn."
3439 << std::endl;
3440 abort_handler(INTERFACE_ERROR);
3441 }
3442 if (hessFlag) {
3443 Cerr << "Error: Hessians not supported in text_book_ouu direct fn."
3444 << std::endl;
3445 abort_handler(INTERFACE_ERROR);
3446 }
3447
3448 size_t i, split = 2 + (numVars - 2)/2; // split the uncertain vars among d1,d2
3449 // xC[0], xC[1] = design
3450 // xC[2], xC[3],... = uncertain
3451
3452 // **** f:
3453 if (directFnASV[0] & 1) {
3454 Real f = 0.;
3455 for(i=2; i<split; ++i)
3456 f += std::pow(xC[i]-10.*xC[0], 4.0);
3457 for(i=split; i<numVars; ++i)
3458 f += std::pow(xC[i]-10.*xC[1], 4.0);
3459 fnVals[0] = f;
3460 }
3461
3462 // **** c1:
3463 if (numFns>1 && (directFnASV[1] & 1))
3464 fnVals[1] = xC[0]*(xC[2]*xC[2] - 0.5*xC[3]);
3465
3466 // **** c2:
3467 if (numFns>2 && (directFnASV[2] & 1))
3468 fnVals[2] = xC[1]*(xC[3]*xC[3] - 0.5*xC[2]);
3469
3470 // **** df/dx:
3471 if (directFnASV[0] & 2)
3472 for (i=0; i<numDerivVars; ++i)
3473 switch (directFnDVV[i]) {
3474 case 1: { // design variable derivative
3475 Real f0 = 0., xC0 = xC[0];
3476 for (size_t j=2; j<split; ++j)
3477 f0 += -40.*std::pow(xC[j]-10.*xC0, 3.0);
3478 fnGrads[0][i] = f0;
3479 break;
3480 }
3481 case 2: { // design variable derivative
3482 Real f1 = 0., xC1 = xC[1];
3483 for (size_t j=split; j<numVars; ++j)
3484 f1 += -40.*std::pow(xC[j]-10.*xC1, 3.0);
3485 fnGrads[0][i] = f1;
3486 break;
3487 }
3488 default: { // uncertain variable derivative
3489 size_t var_index = directFnDVV[i] - 1;
3490 Real xCvi = xC[var_index];
3491 fnGrads[0][i] = (var_index<split) ? 4.*std::pow(xCvi-10.*xC[0], 3)
3492 : 4.*std::pow(xCvi-10.*xC[1], 3);
3493 break;
3494 }
3495 }
3496
3497 // **** dc1/dx:
3498 if (numFns>1 && (directFnASV[1] & 2))
3499 for (i=0; i<numDerivVars; ++i)
3500 switch (directFnDVV[i]) {
3501 case 1: // design variable derivative
3502 fnGrads[1][i] = xC[2]*xC[2] - 0.5*xC[3]; break;
3503 case 3: // uncertain variable derivative
3504 fnGrads[1][i] = 2*xC[0]*xC[2]; break;
3505 case 4: // uncertain variable derivative
3506 fnGrads[1][i] = -0.5*xC[0]; break;
3507 default: // all other derivatives
3508 fnGrads[1][i] = 0.; break;
3509 }
3510
3511 // **** dc2/dx:
3512 if (numFns>2 && (directFnASV[2] & 2))
3513 for (i=0; i<numDerivVars; ++i)
3514 switch (directFnDVV[i]) {
3515 case 2: // design variable derivative
3516 fnGrads[2][i] = xC[3]*xC[3] - 0.5*xC[2]; break;
3517 case 3: // uncertain variable derivative
3518 fnGrads[2][i] = -0.5*xC[1]; break;
3519 case 4: // uncertain variable derivative
3520 fnGrads[2][i] = 2*xC[1]*xC[3]; break;
3521 default: // all other derivative
3522 fnGrads[2][i] = 0.; break;
3523 }
3524
3525 // for faking a more expensive evaluation:
3526 //std::this_thread::sleep_for(std::chrono::seconds(5));
3527
3528 return 0; // no failure
3529 }
3530
3531
scalable_text_book()3532 int TestDriverInterface::scalable_text_book()
3533 {
3534 if (numADIV || numADRV) {
3535 Cerr << "Error: scalable_text_book direct fn does not support discrete "
3536 << "variables." << std::endl;
3537 abort_handler(INTERFACE_ERROR);
3538 }
3539
3540 // **********************************
3541 // **** f: sum (x[i] - POWVAL)^4 ****
3542 // **********************************
3543 size_t i, j;
3544 if (directFnASV[0] & 1) {
3545 fnVals[0] = 0.;
3546 for (i=0; i<numACV; ++i)
3547 fnVals[0] += std::pow(xC[i]-POW_VAL, 4);
3548 }
3549
3550 // ****************
3551 // **** df/dx: ****
3552 // ****************
3553 if (directFnASV[0] & 2) {
3554 std::fill_n(fnGrads[0], fnGrads.numRows(), 0.);
3555 for (i=0; i<numDerivVars; ++i) {
3556 size_t var_index = directFnDVV[i] - 1;
3557 fnGrads[0][i] = 4.*std::pow(xC[var_index]-POW_VAL,3);
3558 }
3559 }
3560
3561 // ********************
3562 // **** d^2f/dx^2: ****
3563 // ********************
3564 if (directFnASV[0] & 4) {
3565 fnHessians[0] = 0.;
3566 for (i=0; i<numDerivVars; ++i) {
3567 size_t var_index = directFnDVV[i] - 1;
3568 fnHessians[0](i,i) = 12.*std::pow(xC[var_index]-POW_VAL,2);
3569 }
3570 }
3571
3572 // *********************
3573 // **** constraints ****
3574 // *********************
3575 // "symmetric" constraint pairs are defined from pairs of variables
3576 // (although odd constraint or variable counts are also allowable):
3577 // for i=1:num_fns-1, c[i] = x[i-1]^2 - x[i]/2 for odd i
3578 // c[i] = x[i-1]^2 - x[i-2]/2 for even i
3579 for (i=1; i<numFns; ++i) {
3580 // ************
3581 // **** c: ****
3582 // ************
3583 if (directFnASV[i] & 1) {
3584 fnVals[i] = (i-1 < numACV) ? xC[i-1]*xC[i-1] : 0.;
3585 if (i%2) // odd constraint
3586 { if (i < numACV) fnVals[i] -= xC[i]/2.; }
3587 else // even constraint
3588 { if (i-2 < numACV) fnVals[i] -= xC[i-2]/2.; }
3589 }
3590 // ****************
3591 // **** dc/dx: ****
3592 // ****************
3593 if (directFnASV[i] & 2) {
3594 std::fill_n(fnGrads[i], fnGrads.numRows(), 0.);
3595 for (j=0; j<numDerivVars; ++j) {
3596 size_t var_index = directFnDVV[j] - 1;
3597 if (i-1 < numACV && var_index == i-1) // both constraints
3598 fnGrads[i][j] = 2.*xC[i-1];
3599 else if (i%2) // odd constraint
3600 { if (i < numACV && var_index == i) fnGrads[i][j] = -0.5; }
3601 else // even constraint
3602 { if (i-2 < numACV && var_index == i-2) fnGrads[i][j] = -0.5; }
3603 }
3604 }
3605 // ********************
3606 // **** d^2c/dx^2: ****
3607 // ********************
3608 if (directFnASV[i] & 4) {
3609 fnHessians[i] = 0.;
3610 if (i-1 < numACV)
3611 for (j=0; j<numDerivVars; ++j)
3612 if (directFnDVV[j] == i) // both constraints
3613 fnHessians[i](j,j) = 2.;
3614 }
3615 }
3616
3617 // for faking a more expensive evaluation:
3618 //std::this_thread::sleep_for(std::chrono::seconds(5));
3619
3620 return 0; // no failure
3621 }
3622
3623
scalable_monomials()3624 int TestDriverInterface::scalable_monomials()
3625 {
3626 if (numADIV || numADRV) {
3627 Cerr << "Error: scalable_monomials direct fn does not support discrete "
3628 << "variables." << std::endl;
3629 abort_handler(INTERFACE_ERROR);
3630 }
3631 if (numFns != 1) {
3632 Cerr << "Error: Bad number of functions in scalable_monomials direct fn."
3633 << std::endl;
3634 abort_handler(INTERFACE_ERROR);
3635 }
3636
3637 // get power of monomial from analysis components, if available (default to 1)
3638 int power = 1;
3639 if (!analysisComponents.empty() &&
3640 !analysisComponents[analysisDriverIndex].empty())
3641 power = std::atoi(analysisComponents[analysisDriverIndex][0].c_str());
3642
3643 // ***************************
3644 // **** f: sum x[i]^power ****
3645 // ***************************
3646 if (directFnASV[0] & 1) {
3647 fnVals[0] = 0.;
3648 for (size_t i=0; i<numACV; ++i)
3649 fnVals[0] += std::pow(xC[i], power);
3650 }
3651
3652 // ****************
3653 // **** df/dx: ****
3654 // ****************
3655 if (directFnASV[0] & 2) {
3656 std::fill_n(fnGrads[0], fnGrads.numRows(), 0.);
3657 for (size_t i=0; i<numDerivVars; ++i) {
3658 size_t var_index = directFnDVV[i] - 1;
3659 switch (power) {
3660 case 0: fnGrads[0][i] = 0; break;
3661 default: fnGrads[0][i] = power*std::pow(xC[var_index], power-1); break;
3662 }
3663 }
3664 }
3665
3666 // ********************
3667 // **** d^2f/dx^2: ****
3668 // ********************
3669 if (directFnASV[0] & 4) {
3670 fnHessians[0] = 0.;
3671 for (size_t i=0; i<numDerivVars; ++i) {
3672 size_t var_index = directFnDVV[i] - 1;
3673 switch (power) {
3674 case 0: case 1:
3675 fnHessians[0](i,i) = 0; break;
3676 default:
3677 fnHessians[0](i,i) = power*(power-1)*std::pow(xC[var_index], power-2);
3678 break;
3679 }
3680 }
3681 }
3682
3683 return 0; // no failure
3684 }
3685
3686
3687
3688
3689
mogatest1()3690 int TestDriverInterface::mogatest1()
3691 {
3692 if (multiProcAnalysisFlag) {
3693 Cerr << "Error: mogatest1 direct fn does not yet support multiprocessor "
3694 << "analyses." << std::endl;
3695 abort_handler(-1);
3696 }
3697 if ( (numACV + numADIV + numADRV) != 3) {
3698 Cerr << "Error: Bad number of variables in mogatest1 direct fn."
3699 << std::endl;
3700 abort_handler(INTERFACE_ERROR);
3701 }
3702 if (numFns != 2) {
3703 Cerr << "Error: Bad number of functions in mogatest1 direct fn."
3704 << std::endl;
3705 abort_handler(INTERFACE_ERROR);
3706 }
3707
3708 // kernel shared with standlone driver in test/mogatest1.cpp
3709 double f0 = 0.0;
3710 double f1 = 0.0;
3711 for (size_t i=0; i<numVars; i++) {
3712 // orders all continuous vars followed by all discrete vars. This is
3713 // fine in the direct case so long as everything is self-consistent.
3714 Real x_i;
3715 if (i<numACV)
3716 x_i = xC[i];
3717 else if (i<numACV+numADIV)
3718 x_i = (Real)xDI[i-numACV];
3719 else
3720 x_i = xDR[i-numACV-numADIV];
3721
3722 f0 = pow(x_i-(1/sqrt(3.0)),2) + f0;
3723 f1 = pow(x_i+(1/sqrt(3.0)),2) + f1;
3724 }
3725 f0 = 1-exp(-1*f0);
3726 f1 = 1-exp(-1*f1);
3727
3728 if (directFnASV[0] & 1)
3729 fnVals[0] = f0;
3730 if (directFnASV[1] & 1)
3731 fnVals[1] = f1;
3732 if (directFnASV[0] & 2 || directFnASV[1] & 2) {
3733 Cerr << "Error: Analytic gradients not supported in mogatest1."
3734 << std::endl;
3735 abort_handler(INTERFACE_ERROR);
3736 }
3737 if (directFnASV[0] & 4 || directFnASV[1] & 4) {
3738 Cerr << "Error: Analytic Hessians not supported in mogatest1."
3739 << std::endl;
3740 abort_handler(INTERFACE_ERROR);
3741 }
3742
3743 return 0; // no failure
3744 }
3745
3746
mogatest2()3747 int TestDriverInterface::mogatest2()
3748 {
3749 if (multiProcAnalysisFlag) {
3750 Cerr << "Error: mogatest2 direct fn does not yet support multiprocessor "
3751 << "analyses." << std::endl;
3752 abort_handler(-1);
3753 }
3754 if (numACV != 2 || numADIV > 0 || numADRV > 0) {
3755 // TODO: allow discrete variables
3756 Cerr << "Error: Bad number of variables in mogatest2 direct fn."
3757 << std::endl;
3758 abort_handler(INTERFACE_ERROR);
3759 }
3760 if (numFns != 2) {
3761 Cerr << "Error: Bad number of functions in mogatest2 direct fn."
3762 << std::endl;
3763 abort_handler(INTERFACE_ERROR);
3764 }
3765
3766 // kernel shared with standlone driver in test/mogatest2.cpp
3767 double f0=0;
3768 double f1=0;
3769 const double PI = 3.14159265358979;
3770 f0 = xC[0];
3771 f1 = sin(2*PI*4*xC[0]);
3772 f1 = f1*(-1*xC[0]/(1+10*xC[1]));
3773 f1 = f1+1-pow((xC[0]/(1+10*xC[1])),2);
3774 f1 = f1*(1+10*xC[1]);
3775
3776 if (directFnASV[0] & 1)
3777 fnVals[0] = f0;
3778 if (directFnASV[1] & 1)
3779 fnVals[1] = f1;
3780 if (directFnASV[0] & 2 || directFnASV[1] & 2) {
3781 Cerr << "Error: Analytic gradients not supported in mogatest2."
3782 << std::endl;
3783 abort_handler(INTERFACE_ERROR);
3784 }
3785 if (directFnASV[0] & 4 || directFnASV[1] & 4) {
3786 Cerr << "Error: Analytic Hessians not supported in mogatest2."
3787 << std::endl;
3788 abort_handler(INTERFACE_ERROR);
3789 }
3790
3791 return 0; // no failure
3792 }
3793
mogatest3()3794 int TestDriverInterface::mogatest3()
3795 {
3796 if (multiProcAnalysisFlag) {
3797 Cerr << "Error: mogatest3 direct fn does not yet support multiprocessor "
3798 << "analyses." << std::endl;
3799 abort_handler(-1);
3800 }
3801 if (numACV != 2 || numADIV > 0 || numADRV > 0) {
3802 // TODO: allow discrete variables
3803 Cerr << "Error: Bad number of variables in mogatest3 direct fn."
3804 << std::endl;
3805 abort_handler(INTERFACE_ERROR);
3806 }
3807 if (numFns != 4) {
3808 Cerr << "Error: Bad number of functions in mogatest3 direct fn."
3809 << std::endl;
3810 abort_handler(INTERFACE_ERROR);
3811 }
3812
3813 // kernel shared with standlone driver in test/mogatest3.cpp
3814 double f0=0;
3815 double f1=0;
3816 double g0=0;
3817 double g1=0;
3818 f0 = pow(xC[0]-2,2)+pow(xC[1]-1,2)+2;
3819 f1 = 9*xC[0]-pow(xC[1]-1,2);
3820 g0 = (xC[0]*xC[0])+(xC[1]*xC[1])-225;
3821 g1 = xC[0]-3*xC[1]+10;
3822
3823 if (directFnASV[0] & 1)
3824 fnVals[0] = f0;
3825 if (directFnASV[1] & 1)
3826 fnVals[1] = f1;
3827 if (directFnASV[2] & 1)
3828 fnVals[2] = g0;
3829 if (directFnASV[3] & 1)
3830 fnVals[3] = g1;
3831 if (directFnASV[0] & 2 || directFnASV[1] & 2 ||
3832 directFnASV[2] & 2 || directFnASV[3] & 2) {
3833 Cerr << "Error: Analytic gradients not supported in mogatest3."
3834 << std::endl;
3835 abort_handler(INTERFACE_ERROR);
3836 }
3837 if (directFnASV[0] & 4 || directFnASV[1] & 4||
3838 directFnASV[2] & 4 || directFnASV[3] & 4) {
3839 Cerr << "Error: Analytic Hessians not supported in mogatest3."
3840 << std::endl;
3841 abort_handler(INTERFACE_ERROR);
3842 }
3843
3844 return 0; // no failure
3845 }
3846
illumination()3847 int TestDriverInterface::illumination()
3848 {
3849 if (multiProcAnalysisFlag) {
3850 Cerr << "Error: illumination direct fn does not yet support multiprocessor "
3851 << "analyses." << std::endl;
3852 abort_handler(-1);
3853 }
3854
3855 if ( (gradFlag || hessFlag) && (numADIV || numADRV) ) {
3856 Cerr << "Error: illumination direct fn assumes no discrete variables in "
3857 << "derivative or hessian mode." << std::endl;
3858 abort_handler(INTERFACE_ERROR);
3859 }
3860
3861 size_t const num_vars = numACV; // keep everything consistent with standalone
3862
3863 if ( num_vars != 7) {
3864 Cerr << "Error: Bad number of variables in illumination direct fn."
3865 << std::endl;
3866 abort_handler(INTERFACE_ERROR);
3867 }
3868
3869 if ( numFns != 1 ) {
3870 Cerr << "Error: Bad number of functions in illumination direct fn."
3871 << std::endl;
3872 abort_handler(INTERFACE_ERROR);
3873 }
3874
3875 // compute function and gradient values
3876 double A[11][7] ={
3877 { 0.347392, 0.205329, 0.191987, 0.077192, 0.004561, 0.024003, 0.000000},
3878 { 0.486058, 0.289069, 0.379202, 0.117711, 0.006667, 0.032256, 0.000000},
3879 { 0.752511, 0.611283, 2.417907, 0.701700, 0.473047, 0.285597, 0.319187},
3880 { 0.303582, 0.364620, 1.898185, 0.693173, 0.607718, 0.328582, 0.437394},
3881 { 0.540946, 0.411549, 1.696545, 0.391735, 0.177832, 0.110119, 0.083817},
3882 { 0.651840, 0.540687, 3.208793, 0.639020, 0.293811, 0.156842, 0.128499},
3883 { 0.098008, 0.245771, 0.742564, 0.807976, 0.929739, 0.435144, 0.669797},
3884 { 0.000000, 0.026963, 0.000000, 0.246606, 0.414657, 0.231777, 0.372202},
3885 { 0.285597, 0.320457, 0.851227, 0.584677, 0.616436, 0.341447, 0.477329},
3886 { 0.324622, 0.306394, 0.991904, 0.477744, 0.376266, 0.158288, 0.198745},
3887 { 0.000000, 0.050361, 0.000000, 0.212042, 0.434397, 0.286455, 0.462731} };
3888
3889 // KRD found bugs in previous computation; this Hessian no longer used:
3890 double harray[7][7] ={
3891 { 1.929437, 1.572662, 6.294004, 1.852205, 1.222324, 0.692036, 0.768564},
3892 { 1.572662, 1.354287, 5.511537, 1.787932, 1.320048, 0.724301, 0.870382},
3893 { 6.294004, 5.511537, 25.064512, 7.358494, 5.133563, 2.791970, 3.257364},
3894 { 1.852205, 1.787932, 7.358494, 2.883178, 2.497491, 1.321922, 1.747230},
3895 { 1.222324, 1.320048, 5.133563, 2.497491, 2.457733, 1.295927, 1.816568},
3896 { 0.692036, 0.724301, 2.791970, 1.321922, 1.295927, 0.694642, 0.968982},
3897 { 0.768564, 0.870382, 3.257364, 1.747230, 1.816568, 0.968982, 1.385357} };
3898
3899
3900 // Calculation of grad(f) = df/dx_I and hess(f) = d^2f/(dx_I dx_J):
3901 //
3902 // U = sum_{i=0:10}{ ( 1 - sum_{j=0:6}{ A[i][j]*x[j] } )^2 }
3903 // dU/dx_I = sum_{i=0:10}{ 2.0*( 1 - sum_{j=0:6}{ A[i][j]*x[j] } )*-A[i][I] }
3904 // d^2U/(dx_I dx_J) = sum_{i=0:10){ 2.0 * A[i][I] * A[i][J] }
3905 //
3906 // f = sqrt(U)
3907 // df/dx_I = 0.5*(dU/dx_I)/sqrt(U)
3908 // = 0.5*(dU/dx_I)/f
3909 // d^2f/(dx_I dx_J)
3910 // = -0.25 * U^(-1.5) * dU/dx_I * dU/dx_J + 0.5/sqrt(U) * d2U/(dx_I dx_J)
3911 // = ( -df/dx_I * df/dx_J + 0.5 * d2U/(dx_I dx_J) ) / f
3912 //
3913 // so to (efficiently) compute grad(f) we precompute f
3914 // and to (efficiently) compute hess(f) we precompute f and grad(f)
3915
3916 size_t Ider, Jder;
3917 double grad[7];
3918 for(Ider=0; Ider<num_vars; ++Ider)
3919 grad[Ider] = 0.0;
3920
3921 // **** f: (f is required to calculate any derivative of f; perform always)
3922 double U = 0.0;
3923 for (size_t i=0; i<11; i++) {
3924 double dtmp = 0.0;
3925 for (size_t j =0; j<num_vars; j++)
3926 dtmp += A[i][j] * xC[j];
3927 dtmp = 1.0 - dtmp;
3928 U += dtmp*dtmp;
3929
3930 // precompute grad(U) unconditionally as it might be needed (cheap)
3931 for(Ider=0; Ider<num_vars; ++Ider)
3932 grad[Ider] -= dtmp*2.0*A[i][Ider]; //this is grad(U)
3933 }
3934 double fx = sqrt(U);
3935 if (directFnASV[0] & 1)
3936 fnVals[0] = fx;
3937
3938 // **** df/dx:
3939 if (directFnASV[0] & 6) { // gradient required for itself and to calcuate the Hessian
3940 for(Ider=0; Ider<num_vars; ++Ider)
3941 grad[Ider] *= (0.5/fx); // this is now updated to be grad(f)
3942
3943 if (directFnASV[0] & 2) { // if ASV demands gradient, print it
3944 for (int Ider=0; Ider<num_vars; ++Ider)
3945 fnGrads[0][Ider] = grad[Ider];
3946 }
3947 }
3948
3949 // **** the hessian ddf/dxIdxJ
3950 if (directFnASV[0] & 4) {
3951 for(Ider=0; Ider<num_vars; ++Ider) {
3952 for(Jder=Ider; Jder<num_vars; ++Jder) {
3953 // define triangular part of hess(U)
3954 for(size_t i=0; i<11; ++i)
3955 fnHessians[0](Ider,Jder) += A[i][Ider]*A[i][Jder]; // this is 0.5*hess(U)
3956
3957 fnHessians[0](Jder,Ider) = fnHessians[0](Ider,Jder) // this is hess(f)
3958 = (fnHessians[0](Ider,Jder)- grad[Ider]*grad[Jder])/fx;
3959 }
3960 }
3961 }
3962 return 0;
3963 }
3964
3965 /// barnes test for SBO perforamnce from Rodriguez, Perez, Renaud, et al.
3966 /// Modified 3/7/18 to incorporate random a[].
barnes()3967 int TestDriverInterface::barnes()
3968 {
3969 if (multiProcAnalysisFlag) {
3970 Cerr << "Error: barnes direct fn does not yet support multiprocessor "
3971 << "analyses." << std::endl;
3972 abort_handler(-1);
3973 }
3974
3975 if ( hessFlag ) {
3976 Cerr << "Error: barnes direct fn does not yet support analytic Hessians."
3977 << std::endl;
3978 abort_handler(INTERFACE_ERROR);
3979 }
3980
3981 if ( gradFlag && (numADIV || numADRV) ) {
3982 Cerr << "Error: barnes direct fn assumes no discrete variables in "
3983 << "derivative mode." << std::endl;
3984 abort_handler(INTERFACE_ERROR);
3985 }
3986
3987 if ( numACV < 2 || numACV > 23) {
3988 Cerr << "Error: Bad number of variables in barnes direct fn."
3989 << std::endl;
3990 abort_handler(INTERFACE_ERROR);
3991 }
3992
3993 if ( numFns != 4 ) {
3994 Cerr << "Error: Bad number of functions in barnes direct fn."
3995 << std::endl;
3996 abort_handler(INTERFACE_ERROR);
3997 }
3998
3999
4000 // Verification test for SBO performance.
4001 // Taken from Rodriguez, Perez, Renaud, et al.
4002 // Constraints g >= 0.
4003
4004 double a[] = { 75.196, -3.8112, 0.12694, -2.0567e-3, 1.0345e-5,
4005 -6.8306, 0.030234, -1.28134e-3, 3.5256e-5, -2.266e-7,
4006 0.25645, -3.4604e-3, 1.3514e-5, -28.106, -5.2375e-6,
4007 -6.3e-8, 7.0e-10, 3.4054e-4, -1.6638e-6, -2.8673,
4008 0.0005};
4009 double x1 = xC[0], x2 = xC[1], x1x2 = x1*x2, x2_sq = x2*x2, x1_sq = x1*x1;
4010
4011 // modify trailing a[] coeffs, according to number of variables
4012 // Assume design followed by uncertain:
4013 // numACV = 2: 2 design, no uncertain
4014 // numACV = 5: 2 design, a[18]--a[20] are uncertain
4015 // numACV = 23: 2 design, a[0] --a[20] are uncertain
4016 size_t i, cntr;
4017 for (i=23-numACV, cntr=2; i<21; ++i, ++cntr)
4018 a[i] = xC[cntr];
4019
4020 // **** f
4021 if (directFnASV[0] & 1) {
4022 double f = a[0] + a[1]*x1 + a[2]*x1_sq + a[3]*x1_sq*x1 + a[4]*x1_sq*x1_sq
4023 + a[5]*x2 + a[6]*x1x2 + a[7]*x1*x1x2 + a[8]*x1x2*x1_sq
4024 + a[9]*x2*x1_sq*x1_sq + a[10]*x2_sq + a[11]*x2*x2_sq + a[12]*x2_sq*x2_sq
4025 + a[13]/(x2+1.) + a[14]*x2_sq*x1_sq + a[15]*x1*x1_sq*x2_sq
4026 + a[16]*x1x2*x2_sq*x1_sq + a[17]*x1*x2_sq + a[18]*x1x2*x2_sq
4027 + a[19]*exp(a[20]*x1x2);
4028 fnVals[0] = f;
4029 }
4030
4031 // **** g1
4032 if (directFnASV[1] & 1)
4033 fnVals[1] = x1x2/700. - 1.;
4034
4035 // **** g2
4036 if (directFnASV[2] & 1)
4037 fnVals[2] = x2/5. - x1_sq/625.;
4038
4039 // **** g3
4040 if (directFnASV[3] & 1)
4041 fnVals[3] = pow(x2/50. - 1., 2.) - x1/500. + 0.11;
4042
4043 // **** df/dx
4044 if (directFnASV[0] & 2) {
4045 for (size_t i=0; i<numDerivVars; i++) {
4046 int var_index = directFnDVV[i] - 1;
4047 switch (var_index) {
4048 case 0:
4049 fnGrads[0][i] = a[1] + 2.*a[2]*x1 + 3.*a[3]*x1_sq + 4.*a[4]*x1_sq*x1
4050 + a[6]*x2 + 2.*a[7]*x1x2 + 3.*a[8]*x2*x1_sq + 4.*a[9]*x1x2*x1_sq
4051 + 2.*a[14]*x2_sq*x1 + 3.*a[15]*x1_sq*x2_sq + 3.*a[16]*x2*x2_sq*x1_sq
4052 + a[17]*x2_sq + a[18]*x2*x2_sq + a[19]*a[20]*x2*exp(a[20]*x1x2);
4053 break;
4054 case 1:
4055 fnGrads[0][i] = a[5] + a[6]*x1 + a[7]*x1_sq + a[8]*x1*x1_sq
4056 + a[9]*x1_sq*x1_sq + 2.*a[10]*x2 + 3.*a[11]*x2_sq + 4.*a[12]*x2*x2_sq
4057 - a[13]/pow(x2+1., 2.) + 2.*a[14]*x2*x1_sq + 2.*a[15]*x1*x1_sq*x2
4058 + 3.*a[16]*x1*x2_sq*x1_sq + 2.*a[17]*x1x2 + 3.*a[18]*x1*x2_sq
4059 + a[19]*a[20]*x1*exp(a[20]*x1x2);
4060 break;
4061 }
4062 }
4063 }
4064
4065 // **** dg1/dx
4066 if (directFnASV[1] & 2) {
4067 for (size_t i=0; i<numDerivVars; i++) {
4068 int var_index = directFnDVV[i] - 1;
4069 switch (var_index) {
4070 case 0:
4071 fnGrads[1][i] = x2/700.;
4072 break;
4073 case 1:
4074 fnGrads[1][i] = x1/700.;
4075 break;
4076 }
4077 }
4078 }
4079
4080 // **** dg2/dx
4081 if (directFnASV[2] & 2) {
4082 for (size_t i=0; i<numDerivVars; i++) {
4083 int var_index = directFnDVV[i] - 1;
4084 switch (var_index) {
4085 case 0:
4086 fnGrads[2][i] = -2.*x1/625.;
4087 break;
4088 case 1:
4089 fnGrads[2][i] = 0.2;
4090 break;
4091 }
4092 }
4093 }
4094
4095 // **** dg3/dx
4096 if (directFnASV[3] & 2) {
4097 for (size_t i=0; i<numDerivVars; i++) {
4098 int var_index = directFnDVV[i] - 1;
4099 switch (var_index) {
4100 case 0:
4101 fnGrads[3][i] = -1./500.;
4102 break;
4103 case 1:
4104 fnGrads[3][i] = 2.*(x2/50. - 1.)/50.;
4105 break;
4106 }
4107 }
4108 }
4109 return 0;
4110 }
4111
4112 /// lo-fi barnes test for SBO perforamnce from Rodriguez, Perez, Renaud, et al.
barnes_lf()4113 int TestDriverInterface::barnes_lf()
4114 {
4115 if (multiProcAnalysisFlag) {
4116 Cerr << "Error: barnes_lf direct fn does not yet support multiprocessor "
4117 << "analyses." << std::endl;
4118 abort_handler(-1);
4119 }
4120
4121 if ( hessFlag ) {
4122 Cerr << "Error: barnes_lf direct fn does not yet support analytic Hessians."
4123 << std::endl;
4124 abort_handler(INTERFACE_ERROR);
4125 }
4126
4127 if ( gradFlag && (numADIV || numADRV) ) {
4128 Cerr << "Error: barnes_lf direct fn assumes no discrete variables in "
4129 << "derivative mode." << std::endl;
4130 abort_handler(INTERFACE_ERROR);
4131 }
4132
4133 if ( numACV != 2) {
4134 Cerr << "Error: Bad number of variables in barnes_lf direct fn."
4135 << std::endl;
4136 abort_handler(INTERFACE_ERROR);
4137 }
4138
4139 if ( numFns != 4 ) {
4140 Cerr << "Error: Bad number of functions in barnes_lf direct fn."
4141 << std::endl;
4142 abort_handler(INTERFACE_ERROR);
4143 }
4144
4145 // Verification test for SBO performance.
4146 // Taken from Rodriguez, Perez, Renaud, et al.
4147 // Constraints g >= 0.
4148
4149 double a[] = { 75.196, -3.8112, 0.12694, -2.0567e-3, 1.0345e-5,
4150 -6.8306, 0.030234, -1.28134e-3, 3.5256e-5, -2.266e-7,
4151 0.25645, -3.4604e-3, 1.3514e-5, -28.106, -5.2375e-6,
4152 -6.3e-8, 7.0e-10, 3.4054e-4, -1.6638e-6, -2.8673,
4153 0.0005};
4154 // Taylor series of Barnes function about the point (p1,p2)
4155 double p1 = 30.0;
4156 double p2 = 40.0;
4157 double x1 = xC[0]-p1;
4158 double x12 = x1*x1;
4159 double x13 = x12*x1;
4160 double x2 = xC[1]-p2;
4161 double x22 = x2*x2;
4162 double x23 = x22*x2;
4163
4164 // **** f
4165 if (directFnASV[0] & 1) {
4166 fnVals[0] =
4167 - 2.74465943148169
4168 + 0.01213957527281*x1
4169 + 0.00995748775273*x12
4170 - 5.557060816484793e-04*x13
4171 + (1.15084419109172+0.00947331101091*x1+2.994070392732408e-05*x12)*x2
4172 + (-0.02997939337414-1.676054720545071e-04*x1)*x22
4173 - 0.00132216646850*x23;
4174 }
4175
4176 // **** g1
4177 if (directFnASV[1] & 1)
4178 fnVals[1] = (xC[0]+xC[1]-50.)/10.;
4179
4180 // **** g2
4181 if (directFnASV[2] & 1)
4182 fnVals[2] = (-0.64*xC[0]+xC[1])/6.;
4183
4184 // **** g3
4185 if (directFnASV[3] & 1) {
4186 if (xC[1] > 50)
4187 fnVals[3] =
4188 -0.00599508167546*xC[0]
4189 + 0.0134054101569*xC[1]
4190 - 0.34054101569933;
4191 else
4192 fnVals[3] =
4193 -0.00599508167546*xC[0]
4194 - 0.01340541015699*xC[1]
4195 + 1.;
4196 }
4197
4198 // **** df/dx
4199 if (directFnASV[0] & 2) {
4200 for (size_t i=0; i<numDerivVars; i++) {
4201 int var_index = directFnDVV[i] - 1;
4202 switch (var_index) {
4203 case 0:
4204 fnGrads[0][i] =
4205 - 0.58530968989099+0.01991497550546*xC[0]-0.00166711824495*x12
4206 + (0.00767686877527+ 5.988140785464816e-05*xC[0])*x2
4207 - 1.676054720545071e-04*x22;
4208 break;
4209 case 1:
4210 fnGrads[0][i] =
4211 0.86664486076442+0.00947331101091*xC[0]+2.994070392732408e-05*x12
4212 + 2*(-0.02495122921250-1.676054720545071e-04*xC[0])*x2
4213 - 0.00396649940550*x22;
4214 break;
4215 }
4216 }
4217 }
4218
4219 // **** dg1/dx
4220 if (directFnASV[1] & 2) {
4221 for (size_t i=0; i<numDerivVars; i++) {
4222 int var_index = directFnDVV[i] - 1;
4223 switch (var_index) {
4224 case 0:
4225 fnGrads[1][i] = 1./10.;
4226 break;
4227 case 1:
4228 fnGrads[1][i] = 1./10.;
4229 break;
4230 }
4231 }
4232 }
4233
4234 // **** dg2/dx
4235 if (directFnASV[2] & 2) {
4236 for (size_t i=0; i<numDerivVars; i++) {
4237 int var_index = directFnDVV[i] - 1;
4238 switch (var_index) {
4239 case 0:
4240 fnGrads[2][i] = -0.64/6.;
4241 break;
4242 case 1:
4243 fnGrads[2][i] = 1./6.;
4244 break;
4245 }
4246 }
4247 }
4248
4249 // **** dg3/dx
4250 if (directFnASV[3] & 2) {
4251 for (size_t i=0; i<numDerivVars; i++) {
4252 int var_index = directFnDVV[i] - 1;
4253 switch (var_index) {
4254 case 0:
4255 fnGrads[3][i] = -0.00599508167546;
4256 break;
4257 case 1:
4258 if (xC[1] > 50)
4259 fnGrads[3][i] = 0.01340541015692;
4260 else
4261 fnGrads[3][i] = -0.01340541015692;
4262 break;
4263 }
4264 }
4265 }
4266 return 0;
4267 }
4268
4269
4270 /// 1D Herbie function and its derivatives (apart from a multiplicative factor)
4271 void TestDriverInterface::
herbie1D(size_t der_mode,Real xc_loc,std::vector<Real> & w_and_ders)4272 herbie1D(size_t der_mode, Real xc_loc, std::vector<Real>& w_and_ders)
4273 {
4274 w_and_ders[0]=w_and_ders[1]=w_and_ders[2]=0.0;
4275
4276 Real rtemp1=xc_loc-1.0;
4277 Real rtemp1_sq=rtemp1*rtemp1;
4278 Real rtemp2=xc_loc+1.0;
4279 Real rtemp2_sq=rtemp2*rtemp2;
4280 Real rtemp3=8.0*(xc_loc+0.1);
4281
4282 if(der_mode & 1) //1=2^0: the 0th derivative of the response (the response itself)
4283 w_and_ders[0]=
4284 std::exp(-rtemp1_sq)
4285 +std::exp(-0.8*rtemp2_sq)
4286 -0.05*std::sin(rtemp3);
4287 if(der_mode & 2) //2=2^1: the 1st derivative of the response
4288 w_and_ders[1]=
4289 -2.0*rtemp1*std::exp(-rtemp1_sq)
4290 -1.6*rtemp2*std::exp(-0.8*rtemp2_sq)
4291 -0.4*std::cos(rtemp3);
4292 if(der_mode & 4) //4=2^2: the 2nd derivative of the response
4293 w_and_ders[2]=
4294 (-2.0+4.0*rtemp1_sq)*std::exp(-rtemp1_sq)
4295 +(-1.6+2.56*rtemp2_sq)*std::exp(-0.8*rtemp2_sq)
4296 +3.2*std::sin(rtemp3);
4297 if(der_mode > 7) {
4298 Cerr << "only 0th through 2nd derivatives are implemented for herbie1D()\n";
4299 assert(false); //should throw an exception get brian to help
4300 }
4301 }
4302
4303
4304 /// 1D Smoothed Herbie= 1DHerbie minus the high frequency sine term, and its derivatives (apart from a multiplicative factor)
4305 void TestDriverInterface::
smooth_herbie1D(size_t der_mode,Real xc_loc,std::vector<Real> & w_and_ders)4306 smooth_herbie1D(size_t der_mode, Real xc_loc, std::vector<Real>& w_and_ders)
4307 {
4308 w_and_ders[0]=w_and_ders[1]=w_and_ders[2]=0.0;
4309
4310 Real rtemp1=xc_loc-1.0;
4311 Real rtemp1_sq=rtemp1*rtemp1;
4312 Real rtemp2=xc_loc+1.0;
4313 Real rtemp2_sq=rtemp2*rtemp2;
4314
4315 if(der_mode & 1) //1=2^0: the 0th derivative of the response (the response itself)
4316 w_and_ders[0]=
4317 std::exp(-rtemp1_sq)
4318 +std::exp(-0.8*rtemp2_sq);
4319 if(der_mode & 2) //2=2^1: the 1st derivative of the response
4320 w_and_ders[1]=
4321 -2.0*rtemp1*std::exp(-rtemp1_sq)
4322 -1.6*rtemp2*std::exp(-0.8*rtemp2_sq);
4323 if(der_mode & 4) //4=2^2: the 2nd derivative of the response
4324 w_and_ders[2]=
4325 (-2.0+4.0*rtemp1_sq)*std::exp(-rtemp1_sq)
4326 +(-1.6+2.56*rtemp2_sq)*std::exp(-0.8*rtemp2_sq);
4327 if(der_mode > 7) {
4328 Cerr << "only 0th through 2nd derivatives are implemented for smooth_herbie1D()\n";
4329 assert(false); //should throw an exception get brian to help
4330 }
4331 }
4332
4333
4334 /// 1D Shubert function and its derivatives (apart from a multiplicative factor)
4335 void TestDriverInterface::
shubert1D(size_t der_mode,Real xc_loc,std::vector<Real> & w_and_ders)4336 shubert1D(size_t der_mode, Real xc_loc, std::vector<Real>& w_and_ders)
4337 {
4338 w_and_ders[0]=w_and_ders[1]=w_and_ders[2]=0.0;
4339
4340 size_t k;
4341 Real k_real;
4342
4343 if(der_mode & 1) {
4344 for (k=1; k<=5; ++k) {
4345 k_real=static_cast<Real>(k);
4346 w_and_ders[0]+=k_real*std::cos(xc_loc*(k_real+1.0)+k_real);
4347 }
4348 }
4349 if(der_mode & 2) {
4350 for (k=1; k<=5; ++k) {
4351 k_real=static_cast<Real>(k);
4352 w_and_ders[1]+=k_real*(k_real+1.0)*-std::sin(xc_loc*(k_real+1.0)+k_real);
4353 }
4354 }
4355 if(der_mode & 4) {
4356 for (k=1; k<=5; ++k) {
4357 k_real=static_cast<Real>(k);
4358 w_and_ders[2]+=k_real*(k_real+1.0)*(k_real+1.0)*-std::cos(xc_loc*(k_real+1.0)+k_real);
4359 }
4360 }
4361 if(der_mode > 7) {
4362 Cerr << "only 0th through 2nd derivatives are implemented for shubert1D()\n";
4363 assert(false); //should throw an exception get brian to help
4364 }
4365 }
4366
4367
4368 /// N-D Herbie function and its derivatives
herbie()4369 int TestDriverInterface::herbie()
4370 {
4371 size_t i;
4372 std::vector<size_t> der_mode(numVars);
4373 for (i=0; i<numVars; ++i)
4374 der_mode[i]=1;
4375 if(directFnASV[0] >= 2)
4376 for (i=0; i<numDerivVars; ++i)
4377 der_mode[directFnDVV[i]-1]+=2;
4378 if(directFnASV[0] >= 4)
4379 for (i=0; i<numDerivVars; ++i)
4380 der_mode[directFnDVV[i]-1]+=4;
4381 std::vector<Real> w(numVars);
4382 std::vector<Real> d1w(numVars);
4383 std::vector<Real> d2w(numVars);
4384 std::vector<Real> w_and_ders(3);
4385
4386 for(i=0; i<numVars; ++i) {
4387 herbie1D(der_mode[i],xC[i],w_and_ders);
4388 w[i] =w_and_ders[0];
4389 d1w[i]=w_and_ders[1];
4390 d2w[i]=w_and_ders[2];
4391 }
4392
4393 separable_combine(-1.0,w,d1w,d2w);
4394 return 0;
4395 }
4396
4397 /// N-D Smoothed Herbie function and its derivatives
smooth_herbie()4398 int TestDriverInterface::smooth_herbie()
4399 {
4400 size_t i;
4401 std::vector<size_t> der_mode(numVars);
4402 for (i=0; i<numVars; ++i)
4403 der_mode[i]=1;
4404 if(directFnASV[0] >= 2)
4405 for (i=0; i<numDerivVars; ++i)
4406 der_mode[directFnDVV[i]-1]+=2;
4407 if(directFnASV[0] >= 4)
4408 for (i=0; i<numDerivVars; ++i)
4409 der_mode[directFnDVV[i]-1]+=4;
4410 std::vector<Real> w(numVars);
4411 std::vector<Real> d1w(numVars);
4412 std::vector<Real> d2w(numVars);
4413 std::vector<Real> w_and_ders(3);
4414
4415 for(i=0; i<numVars; ++i) {
4416 smooth_herbie1D(der_mode[i], xC[i], w_and_ders);
4417 w[i] =w_and_ders[0];
4418 d1w[i]=w_and_ders[1];
4419 d2w[i]=w_and_ders[2];
4420 }
4421
4422 separable_combine(-1.0,w,d1w,d2w);
4423 return 0;
4424 }
4425
shubert()4426 int TestDriverInterface::shubert()
4427 {
4428 size_t i;
4429 std::vector<size_t> der_mode(numVars);
4430 for (i=0; i<numVars; ++i)
4431 der_mode[i]=1;
4432 if(directFnASV[0] >= 2)
4433 for (i=0; i<numDerivVars; ++i)
4434 der_mode[directFnDVV[i]-1]+=2;
4435 if(directFnASV[0] >= 4)
4436 for (i=0; i<numDerivVars; ++i)
4437 der_mode[directFnDVV[i]-1]+=4;
4438 std::vector<Real> w(numVars);
4439 std::vector<Real> d1w(numVars);
4440 std::vector<Real> d2w(numVars);
4441 std::vector<Real> w_and_ders(3);
4442
4443 for(i=0; i<numVars; ++i) {
4444 shubert1D(der_mode[i],xC[i],w_and_ders);
4445 w[i] =w_and_ders[0];
4446 d1w[i]=w_and_ders[1];
4447 d2w[i]=w_and_ders[2];
4448 }
4449
4450 separable_combine(1.0,w,d1w,d2w);
4451 return 0;
4452 }
4453
bayes_linear()4454 int TestDriverInterface::bayes_linear()
4455 {
4456 // This test driver implements the linear verification example in the
4457 // document "User Guidelines and Best Practices for CASL VUQ Analysis
4458 // using Dakota", CASL-U-2014-0038-000. The example is in Appendix A
4459 // and Section 6.2.6. It is of the form:
4460 // y = g(x)*Beta + epsilon, where epsilon ~ N(0,1/lamda * R)
4461 // Beta is a set of unknown regression coefficients we are trying to infer,
4462 // epsilon is the vector of observational errors having variance
4463 // 1/lambda = sigma^2
4464 // and R is a fixed correlation function.
4465
4466 if (multiProcAnalysisFlag) {
4467 Cerr << "Error: bayes_linear direct fn does not support "
4468 << "multiprocessor analyses." << std::endl;
4469 abort_handler(-1);
4470 }
4471 if (numVars < 1 || numVars > 500 || numADIV || numADRV) {
4472 Cerr << "Error: Bad variable types in Bayes linear fn."
4473 << std::endl;
4474 abort_handler(INTERFACE_ERROR);
4475 }
4476 if (numFns < 1) {
4477 Cerr << "Error: Bad number of functions in Bayes linear direct fn."
4478 << std::endl;
4479 abort_handler(INTERFACE_ERROR);
4480 }
4481 if (hessFlag || gradFlag) {
4482 Cerr << "Error: Gradients and Hessians not supported in Bayes linear "
4483 << "direct fn." << std::endl;
4484 abort_handler(INTERFACE_ERROR);
4485 }
4486
4487 /*const Real pi = 3.14159265358979324;
4488
4489 Real mean_pred = 0.4; int i;
4490 for (i=1; i<numVars; i++) {
4491 mean_pred += 0.4 + 0.5*std::sin(2*pi*i/numVars);
4492 }
4493 RealVector n_means, n_std_devs, n_l_bnds, n_u_bnds;
4494 n_means.resize(numFns); n_std_devs.resize(numFns);
4495 n_l_bnds.resize(numFns); n_u_bnds.resize(numFns);
4496 Real dbl_inf = std::numeric_limits<Real>::infinity();
4497 for (i=0; i<numFns; i++) {
4498 n_means(i)=mean_pred;
4499 n_std_devs(i)=0.0316; // initial lambda = 1000.
4500 n_l_bnds(i)=-dbl_inf;
4501 n_u_bnds(i)= dbl_inf;
4502 }
4503
4504 NonDLHSSampling* normal_sampler;
4505 String rng("mt19937");
4506 unsigned short sample_type = SUBMETHOD_LHS;
4507 int seed = 1;
4508 seed += (int)clock();
4509 Cout << "seed " << seed;
4510 RealSymMatrix correl_matrix(numFns);
4511 correl_matrix = 0.8;
4512 for (i=0; i<numFns; i++)
4513 correl_matrix(i,i)=1.0;
4514
4515 Cout << "correl matrix " << correl_matrix << '\n';
4516
4517 normal_sampler = new NonDLHSSampling(sample_type, 1, seed,
4518 rng, n_means, n_std_devs,
4519 n_l_bnds, n_u_bnds, correl_matrix);
4520
4521 Cout << "normal samples " << normal_sampler->all_samples() << '\n';
4522 const RealMatrix& lhs_samples = normal_sampler->all_samples();
4523 RealVector temp_cvars = Teuchos::getCol(Teuchos::View, const_cast<RealMatrix&>(lhs_samples), 0);
4524
4525 for (i=0; i<numFns; i++) {
4526 fnVals[i]=temp_cvars[i];
4527 }
4528 */
4529
4530 Real pred=0.0;
4531 for (int i=0; i<numVars; i++)
4532 pred += xC[i];
4533 fnVals[0]= pred;
4534
4535 return 0;
4536 }
4537
problem18()4538 int TestDriverInterface::problem18(){
4539 // This test driver implements the 1D optimization benchmark problem 18 from
4540 // http://infinity77.net/global_optimization/test_functions_1d.html
4541 // as a benchmark function to test the MLMC approach with SNOWPAC
4542
4543 if (multiProcAnalysisFlag) {
4544 Cerr << "Error: problem18 direct fn does not support "
4545 << "multiprocessor analyses." << std::endl;
4546 abort_handler(-1);
4547 }
4548 //if (numVars < 1 || numVars > 500 || numADIV || numADRV) {
4549 // Cerr << "Error: Bad variable types in problem18 fn."
4550 // << std::endl;
4551 // abort_handler(INTERFACE_ERROR);
4552 //} Unsure about this yet
4553 if (numFns < 1) {
4554 Cerr << "Error: Bad number of functions in problem18 direct fn."
4555 << std::endl;
4556 abort_handler(INTERFACE_ERROR);
4557 }
4558 if (hessFlag || gradFlag) {
4559 Cerr << "Error: Gradients and Hessians not supported in problem18 "
4560 << "direct fn." << std::endl;
4561 abort_handler(INTERFACE_ERROR);
4562 }
4563
4564 std::map<var_t, Real>::iterator x_iter = xCM.find(VAR_x);
4565 Real x = (x_iter == xCM.end()) ? 0.5 : x_iter->second; // x
4566
4567 std::map<var_t, Real>::iterator xi_iter = xCM.find(VAR_xi);
4568 Real xi = (xi_iter == xCM.end()) ? 0. : xi_iter->second; // RV xi
4569
4570 std::map<var_t, Real>::iterator Af_iter = xDRM.find(VAR_Af);
4571 Real Af = (Af_iter == xDRM.end()) ? 1. : Af_iter->second; // Correlation Af for objective
4572
4573 std::map<var_t, Real>::iterator Ac_iter = xDRM.find(VAR_Ac);
4574 Real Ac = (Ac_iter == xDRM.end()) ? 1. : Ac_iter->second; // Correlation Ac for constraint
4575
4576 if(Af < 0){
4577 Af = problem18_Ax(Af, x);
4578 }
4579 if(Ac < 0){
4580 Ac = problem18_Ax(Ac, x);
4581 }
4582
4583 fnVals[0] = problem18_f(x) + Af * xi * xi * xi;
4584 fnVals[1] = problem18_g(x) - problem18_f(x) + Ac * xi * xi * xi;
4585
4586 //std::cout << "Input parameters:" << std::endl;
4587 //std::cout << "x: " << x << std::endl;
4588 //std::cout << "xi: " << xi << std::endl;
4589 //std::cout << "Af: " << Af << std::endl;
4590 //std::cout << "Ac: " << Ac << std::endl;
4591 //std::cout << "x reached end: " << (x_iter == xCM.end()) << std::endl;
4592 //std::cout << "xi reached end: " << (xi_iter == xCM.end()) << std::endl;
4593 //std::cout << "Af reached end: " << (Af_iter == xDRM.end()) << std::endl;
4594 //std::cout << "Ac reached end: " << (Ac_iter == xDRM.end()) << std::endl;
4595
4596 return 0;
4597 }
4598
problem18_f(const double & x)4599 double TestDriverInterface::problem18_f(const double &x){
4600 return x <= 3. ?
4601 (x - 2.) * (x - 2.) :
4602 2. * std::log(x - 2.) + 1;
4603 }
4604
problem18_g(const double & x)4605 double TestDriverInterface::problem18_g(const double &x){
4606 return (2. * std::log(3.5 - 2))/(3.5 - 1) * x + 1 - (2. *std::log(3.5 - 2))/(3.5 - 1);
4607 }
4608
problem18_Ax(const double & A,const double & x)4609 double TestDriverInterface::problem18_Ax(const double &A, const double &x){
4610 if(A == -1)
4611 return 0.5/6. * x + 0.4;
4612 else if(A == -2)
4613 return 0.5/6. * sin(x) + 0.4;
4614 else if(A == -3)
4615 return 0.5/6. * log(x) + 0.4;
4616 else if(A == -4)
4617 return 0.69*1./exp(2.*x)+0.3;
4618 else
4619 throw INTERFACE_ERROR;
4620 }
4621
4622
4623 /// this function combines N 1D functions and their derivatives to compute a N-D separable function and its derivatives, logic is general enough to support different 1D functions in different dimensions (can mix and match)
separable_combine(Real mult_scale_factor,std::vector<Real> & w,std::vector<Real> & d1w,std::vector<Real> & d2w)4624 void TestDriverInterface::separable_combine(Real mult_scale_factor, std::vector<Real>& w, std::vector<Real>& d1w, std::vector<Real>& d2w)
4625 {
4626 // *************************************************************
4627 // **** now that w(x_i), dw(x_i)/dx_i, and d^2w(x_i)/dx_i^2 ****
4628 // **** are defined we can calculate the response, gradient ****
4629 // **** of the response, and Hessian of the response in an ****
4630 // **** identical fashion ****
4631 // *************************************************************
4632
4633 Real local_val;
4634 size_t i, j, k, i_var_index, j_var_index;
4635
4636 // ****************************************
4637 // **** response ****
4638 // **** f=\prod_{i=1}^{numVars} w(x_i) ****
4639 // ****************************************
4640 if (directFnASV[0] & 1) {
4641 local_val=mult_scale_factor;
4642 for (i=0; i<numVars; ++i)
4643 local_val*=w[i];
4644 fnVals[0]=local_val;
4645 }
4646
4647 // **************************************************
4648 // **** gradient of response ****
4649 // **** df/dx_i=(\prod_{j=1}^{i-1} w(x_j)) ... ****
4650 // **** *(dw(x_i)/dx_i) ... ****
4651 // **** *(\prod_{j=i+1}^{numVars} w(x_j)) ****
4652 // **************************************************
4653 if (directFnASV[0] & 2) {
4654 std::fill_n(fnGrads[0], fnGrads.numRows(), 0.);
4655 for (i=0; i<numDerivVars; ++i) {
4656 i_var_index = directFnDVV[i] - 1;
4657 local_val=mult_scale_factor*d1w[i_var_index];
4658 for (j=0; j<i_var_index; ++j)
4659 local_val*=w[j];
4660 for (j=i_var_index+1; j<numVars; ++j)
4661 local_val*=w[j];
4662 fnGrads[0][i]=local_val;
4663 }
4664 }
4665
4666 // ***********************************************************
4667 // **** Hessian of response ****
4668 // **** if(i==j) ****
4669 // **** d^2f/dx_i^2=(\prod_{k=1}^{i-1} w(x_k)) ... ****
4670 // **** *(d^2w(x_i)/dx_i^2) ... ****
4671 // **** *(\prod_{k=i+1}^{numVars} w(x_k)) ****
4672 // **** ****
4673 // **** if(i<j) ****
4674 // **** d^2f/(dx_i*dx_j)=(\prod_{k=1}^{i-1} w(x_k)) ... ****
4675 // **** *(dw(x_i)/dx_i) ... ****
4676 // **** *(\prod_{k=i+1}^{j-1} w(x_k)) ... ****
4677 // **** *(dw(x_j)/dx_j) ... ****
4678 // **** *(\prod_{k=j+1}^{numVars} w(x_j)) ****
4679 // ***********************************************************
4680 if (directFnASV[0] & 4) {
4681 fnHessians[0] = 0.; //what does this do? I think it's vestigal
4682 for (i=0; i<numDerivVars; ++i) {
4683 i_var_index = directFnDVV[i] - 1;
4684 for (j=0; j<numDerivVars; ++j) {
4685 j_var_index = directFnDVV[j] - 1;
4686 if (i_var_index==j_var_index )
4687 local_val = mult_scale_factor*d2w[i_var_index];
4688 else
4689 local_val = mult_scale_factor*d1w[i_var_index]*d1w[j_var_index];
4690 for (k=0; k<numVars; ++k)
4691 if( (k!=i_var_index) && (k!=j_var_index) )
4692 local_val*=w[k];
4693 fnHessians[0](i,j) =local_val;
4694 }
4695 }
4696 }
4697 }
4698
levenshtein_distance(const String & v)4699 Real TestDriverInterface::levenshtein_distance(const String &v) {
4700 /// Levenshtein distance is the number of changes (single character
4701 // deletions, additions, or changes) to convert one string (v) to another
4702 // (levenshteinReference). This nice implementation shamelessly stolen/adapted
4703 // from Wikipedia, which attributes it to Wager and Fischer,
4704 // doi:10.1145/321796.321811.
4705 // Results are stored in levenshteinDistanceCache to avoid needless repeated
4706 // calcuations
4707 SRMCIter d_match;
4708 d_match = levenshteinDistanceCache.find(v);
4709 if (d_match != levenshteinDistanceCache.end())
4710 return d_match->second;
4711 const String &w = LEV_REF;
4712 size_t v_len = v.size(), w_len = w.size();
4713 IntMatrix d(v_len+1, w_len+1);
4714 size_t i, j;
4715 for(i=0; i<=v_len; ++i)
4716 d(i,0) = i;
4717 for(j=0; j<=w_len; ++j)
4718 d(0,j) = j;
4719 for(j=1; j<=w_len; ++j) {
4720 for(i=1; i<=v_len; ++i) {
4721 if(v[i-1] == w[j-1])
4722 d(i,j) = d(i-1,j-1);
4723 else
4724 d(i,j) = std::min( d(i-1,j)+1, std::min(d(i ,j-1) + 1,
4725 d(i-1,j-1) + 1));
4726 }
4727 }
4728 levenshteinDistanceCache[v] = d(v_len,w_len);
4729 return levenshteinDistanceCache[v];
4730 }
4731
4732 #ifdef DAKOTA_SALINAS
salinas()4733 int TestDriverInterface::salinas()
4734 {
4735 if (numFns < 1 || numFns > 3) {
4736 Cerr << "Error: Bad number of functions in salinas direct fn." << std::endl;
4737 abort_handler(INTERFACE_ERROR);
4738 }
4739 if (numVars < 1) {
4740 Cerr << "Error: Bad number of variables in salinas direct fn." << std::endl;
4741 abort_handler(INTERFACE_ERROR);
4742 }
4743 if (gradFlag || hessFlag) {
4744 Cerr << "Error: analytic derivatives not yet supported in salinas direct "
4745 << "fn." <<std::endl;
4746 abort_handler(INTERFACE_ERROR);
4747 }
4748
4749 // ------------------------
4750 // Salinas input processing
4751 // ------------------------
4752
4753 // Set up dummy argc and argv with name of modified input file.
4754 // NOTE: for concurrent analyses within each fn eval, may want something like
4755 // salinas.<eval>.<analysis>.inp (e.g., salinas.3.2.inp)
4756 int argc = 2;
4757 char* argv[3];
4758 argv[0] = "salinas"; // should be ignored
4759 char si[32];
4760 std::sprintf(si,"salinas%d.inp", evalServerId); // tag root name in root.inp
4761 argv[1] = si;
4762 argv[2] = NULL; // standard requires this, see Kern.&Ritchie, p. 115
4763
4764 // Insert vars into Salinas input file (Exodus model file avoided for now).
4765 // Set up loop to process input file and match tokens to variable tags. The
4766 // Salinas parser is not dependent on new lines, so don't worry about
4767 // formatting.
4768 std::ifstream fin("salinas.inp.template");
4769 std::ofstream fout(argv[1]);
4770 std::string token;
4771 while (fin) {
4772 fin >> token;
4773 if (token=="//")
4774 fin.ignore(256, '\n'); // comments will not be replicated in fout
4775 else if (token=="'./tagged.exo'") {
4776 // **** Issues with the Exodus input file. Exodus input files must be
4777 // tagged because the Exodus output uses the same name and must be
4778 // protected from conflict with other concurrent simulations. This
4779 // requires the creation of these tagged files by DAKOTA or their
4780 // existence a priori (which is a problem when tagging with open ended
4781 // indices like evaluation id). A few approaches to handling this:
4782 // 1.) could use system("cp root.exo root#.exo"), but no good on TFLOP
4783 // 2.) could tag w/ evalServerId & put sal[0-9].exo out before launching,
4784 // but Salinas must overwrite properly (it does) & data loss must be OK
4785 // 3.) could modify salinas to use (tagged) root.inp i/o root.exo in
4786 // creating root-out.exo, thereby removing the need to tag Exodus input
4787 char se[32];
4788 std::sprintf(se,"'./salinas%d.exo' ", evalServerId); // tag root in root.exo
4789 fout << se;
4790 }
4791 else if (localDataView & VARIABLES_MAP) {
4792 // Use a map-based lookup if tokens of interest are added to varTypeMap
4793 std::map<String, var_t>::iterator v_it = varTypeMap.find(token);
4794 if (v_it == varTypeMap.end())
4795 fout << token << ' ';
4796 else {
4797 var_t vt = v_it->second;
4798 std::map<var_t, Real>::iterator xc_it = xCM.find(vt),
4799 xdr_it = xDRM.find(vt);
4800 std::map<var_t, int>::iterator xdi_it = xDIM.find(vt);
4801 if (xc_it != xCM.end())
4802 fout << xc_it->second << ' ';
4803 else if (xdi_it != xDIM.end())
4804 fout << xdi_it->second << ' ';
4805 else if (xdr_it != xDRM.end())
4806 fout << xdr_it->second << ' ';
4807 elseint TestDriverInterface::scalable_text_book()
4808
4809 fout << token << ' ';
4810 }
4811 }
4812 else {
4813 bool found = false;
4814 size_t i;
4815 for (i=0; i<numACV; ++i) // loop to remove any order dependency
4816 if (token == xCLabels[i]) // detect variable label
4817 { fout << xC[i] << ' '; found = true; break; }
4818 if (!found)
4819 for (i=0; i<numADIV; ++i) // loop to remove any order dependency
4820 if (token == xDILabels[i]) // detect variable label
4821 { fout << xDI[i] << ' '; found = true; break; }
4822 if (!found)
4823 for (i=0; i<numADRV; ++i) // loop to remove any order dependency
4824 if (token == xDRLabels[i]) // detect variable label
4825 { fout << xDR[i] << ' '; found = true; break; }
4826 if (!found)
4827 fout << token << ' ';
4828 }
4829 }
4830 fout << std::flush;
4831
4832 // -----------------
4833 // Salinas execution
4834 // -----------------
4835
4836 // salinas_main is a bare bones wrapper for salinas. It is provided to
4837 // permit calling salinas as a subroutine.
4838
4839 // analysis_comm may be invalid if multiProcAnalysisFlag is not true!
4840 const ParallelLevel& ea_pl
4841 = parallelLib.parallel_configuration().ea_parallel_level();
4842 MPI_Comm analysis_comm = ea_pl.server_intra_communicator();
4843 int fail_code = salinas_main(argc, argv, &analysis_comm);
4844 if (fail_code)
4845 return fail_code;
4846
4847 // ---------------------------
4848 // Salinas response processing
4849 // ---------------------------
4850
4851 // Compute margins and return lowest margin as objective function to be
4852 // _maximized_ (minimize opposite sign). Constrain mass to be:
4853 // mass_low_bnd <= mass <= mass_upp_bnd
4854 //Real min_margin = 0.;
4855 Real lambda1, mass, mass_low_bnd=3., mass_upp_bnd=6.; // 4.608 nominal mass
4856
4857 // Call EXODUS utilities to retrieve stress & g data
4858
4859 // Retrieve data from salinas#.rslt
4860 char so[32];
4861 std::sprintf(so,"salinas%d.rslt",evalServerId);
4862 std::ifstream f2in(so);
4863 while (f2in) {
4864 f2in >> token;
4865 if (token=="Total") {
4866 for (size_t i=0; i<4; ++i) // After "Total", parse "Mass of Structure is"
4867 f2in >> token;
4868 f2in >> mass;
4869 }
4870 else if (token=="1:") {
4871 f2in >> lambda1;
4872 }
4873 else
4874 f2in.ignore(256, '\n');
4875 }
4876
4877 // **** f:
4878 if (directFnASV[0] & 1)
4879 fnVals[0] = -lambda1; // max fundamental freq. s.t. mass bounds
4880 //fnVals[0] = -min_margin;
4881
4882 // **** c1:
4883 if (numFns > 1 && (directFnASV[1] & 1))
4884 fnVals[1] = (mass_low_bnd - mass)/mass_low_bnd;
4885
4886 // **** c2:
4887 if (numFns > 2 && (directFnASV[2] & 1))
4888 fnVals[2] = (mass - mass_upp_bnd)/mass_upp_bnd;
4889
4890 // **** df/dx:
4891 //if (directFnASV[0] & 2) {
4892 //}
4893
4894 // **** dc1/dx:
4895 //if (numFns > 1 && (directFnASV[1] & 2)) {
4896 //}
4897
4898 // **** dc2/dx:
4899 //if (numFns > 2 && (directFnASV[2] & 2)) {
4900 //}
4901
4902 return 0;
4903 }
4904 #endif // DAKOTA_SALINAS
4905
4906
4907 #ifdef DAKOTA_MODELCENTER
4908 /** The ModelCenter interface doesn't have any specific construct
4909 vs. run time functions. For now, we manage it along with the
4910 integrated test drivers */
mc_api_run()4911 int TestDriverInterface::mc_api_run()
4912 {
4913 // ModelCenter interface through direct Dakota interface, HKIM 4/3/03
4914 // Modified to replace pxcFile with analysisComponents, MSE 6/20/05
4915
4916 if (multiProcAnalysisFlag) {
4917 Cerr << "Error: mc_api_run direct fn does not yet support multiprocessor "
4918 << "analyses." << std::endl;
4919 abort_handler(-1);
4920 }
4921
4922 int i, j, ireturn, iprint = 1;
4923
4924 if(!mc_ptr_int) { // If null, instantiate ModelCenter
4925 // the pxcFile (Phoenix configuration file) is passed through the
4926 // analysis_components specification
4927 if (!analysisComponents.empty() &&
4928 !analysisComponents[analysisDriverIndex].empty())
4929 mc_load_model(ireturn,iprint,mc_ptr_int,
4930 analysisComponents[analysisDriverIndex][0].c_str());
4931 else
4932 ireturn = -1;
4933 if(ireturn == -1 || mc_ptr_int == 0) abort_handler(INTERFACE_ERROR);
4934 }
4935
4936 // continuous variables
4937 for(i=0; i<numACV; ++i) {
4938 const char* inStr = xCLabels[i].c_str();
4939 mc_set_value(ireturn,iprint,mc_ptr_int,xC[i],inStr);
4940 if(ireturn == -1) abort_handler(INTERFACE_ERROR);
4941 }
4942
4943 // discrete, integer-valued variables (actual values sent, not indices)
4944 for(i=0; i<numADIV; ++i) {
4945 const char* inStr = xDILabels[i].c_str();
4946 mc_set_value(ireturn,iprint,mc_ptr_int,xDI[i],inStr);
4947 if(ireturn == -1) abort_handler(INTERFACE_ERROR);
4948 }
4949
4950 // discrete, real-valued variables (actual values sent, not indices)
4951 for(i=0; i<numADRV; ++i) {
4952 const char* inStr = xDRLabels[i].c_str();
4953 mc_set_value(ireturn,iprint,mc_ptr_int,xDR[i],inStr);
4954 if(ireturn == -1) abort_handler(INTERFACE_ERROR);
4955 }
4956
4957 int out_var_act_len = fnLabels.size();
4958 if (out_var_act_len != numFns) {
4959 Cerr << "Error: Mismatch in the number of responses in mc_api_run."
4960 << std::endl;
4961 abort_handler(INTERFACE_ERROR);
4962 }
4963
4964 for (i=0; i<out_var_act_len; ++i) {
4965 // **** f:
4966 if (directFnASV[i] & 1) {
4967 const char* outStr = fnLabels[i].c_str();
4968 mc_get_value(ireturn,iprint,mc_ptr_int,fnVals[i],outStr);
4969 if(ireturn == -1) {
4970 // Assume this is a failed function evaluation
4971 // TODO: check correctness / other possible return codes
4972 return(-1);
4973 }
4974
4975 }
4976 // **** df/dx:
4977 if (directFnASV[i] & 2) {
4978 Cerr << "Error: Analytic gradients not supported in mc_api_run."
4979 << std::endl;
4980 abort_handler(INTERFACE_ERROR);
4981 }
4982 // **** d^2f/dx^2:
4983 if (directFnASV[i] & 4) {
4984 Cerr << "Error: Analytic Hessians not supported in mc_api_run."
4985 << std::endl;
4986 abort_handler(INTERFACE_ERROR);
4987 }
4988
4989 }
4990
4991 if(dc_ptr_int) {
4992 mc_store_current_design_point(ireturn,iprint,dc_ptr_int);
4993 // ignore ireturn for now.
4994 }
4995
4996 return 0; // no failure
4997 }
4998 #endif // DAKOTA_MODELCENTER
4999
aniso_quad_form()5000 int TestDriverInterface::aniso_quad_form()
5001 {
5002 if (multiProcAnalysisFlag) {
5003 Cerr << "Error: aniso_quad_form direct fn does not yet support "
5004 << "multiprocessor analyses." << std::endl;
5005 abort_handler(-1);
5006 }
5007 if (numADIV || numADRV || (gradFlag && numDerivVars != numVars)) {
5008 Cerr << "Error: Bad number of variables in aniso_quad_form direct fn."
5009 << std::endl;
5010 abort_handler(INTERFACE_ERROR);
5011 }
5012 if (numFns != 1) {
5013 Cerr << "Error: Bad number of functions in aniso_quad_form direct fn."
5014 << std::endl;
5015 abort_handler(INTERFACE_ERROR);
5016 }
5017 if (hessFlag) {
5018 Cerr << "Error: Hessians not supported in aniso_quad_form direct fn."
5019 << std::endl;
5020 abort_handler(INTERFACE_ERROR);
5021 }
5022
5023 static bool initialized = false;
5024
5025 static RealMatrix q_mat(numVars, numVars);
5026
5027 if(!initialized)
5028 {
5029 size_t seed = std::time(NULL);
5030 std::vector<RealMatrix::scalarType> eigenvals =
5031 boost::assign::list_of(7.0)(4.0)(1.0);
5032
5033 if (!analysisComponents.empty() &&
5034 !analysisComponents[analysisDriverIndex].empty())
5035 {
5036 typedef boost::char_separator<char> sepT;
5037 typedef boost::tokenizer<sepT> tokenT;
5038 sepT sep(" :");
5039
5040 const StringArray& anal_comps = analysisComponents[analysisDriverIndex];
5041 for(StringArray::const_iterator comp = anal_comps.begin();
5042 comp != anal_comps.end(); ++comp)
5043 {
5044 tokenT tokens(*comp, sep);
5045 tokenT::iterator tok = tokens.begin();
5046 if(*tok == "seed")
5047 {
5048 if(++tok == tokens.end())
5049 {
5050 Cerr << "Seed not specified for aniso_quad_form." << std::endl;
5051 abort_handler(INTERFACE_ERROR);
5052 }
5053
5054 seed = static_cast<size_t>(std::stoull(*tok));
5055
5056 if(++tok != tokens.end())
5057 {
5058 Cerr << "Multiple fields given as a seed in analysis_components "
5059 "for aniso_quad_form." << std::endl;
5060 abort_handler(INTERFACE_ERROR);
5061 }
5062 }
5063 else if(*tok == "eigenvals")
5064 {
5065 eigenvals.clear();
5066
5067 ++tok;
5068 for(; tok != tokens.end(); ++tok)
5069 {
5070 eigenvals.push_back(stod(*tok));
5071 }
5072
5073 if(eigenvals.size() == 0)
5074 {
5075 Cerr << "No eigenvalues specified in analysis_components "
5076 "for aniso_quad_form." << std::endl;
5077 abort_handler(INTERFACE_ERROR);
5078 }
5079 else if(eigenvals.size() > numVars)
5080 {
5081 Cerr << "Too many eigenvalues specified in analysis_components "
5082 "for aniso_quad_form." << std::endl;
5083 abort_handler(INTERFACE_ERROR);
5084 }
5085 }
5086 else
5087 {
5088 Cerr << "Aniso_quad_form received invalid analysis_components."
5089 << std::endl;
5090 abort_handler(INTERFACE_ERROR);
5091 }
5092 }
5093 }
5094
5095 boost::random::mt19937 vec_RNG(seed);
5096 boost::random::normal_distribution<> sampler;
5097
5098 size_t num_prin_direc = eigenvals.size();
5099
5100 RealMatrix prin_vecs(numVars, num_prin_direc);
5101 size_t numElem = numVars * num_prin_direc;
5102
5103 RealMatrix::scalarType* vals = prin_vecs.values();
5104 for(size_t i = 0; i < numElem; ++i)
5105 {
5106 vals[i] = sampler(vec_RNG);
5107 }
5108
5109 for(int i = 0; i < num_prin_direc; ++i)
5110 {
5111 RealVector prev_ortho_vec = getCol(Teuchos::View, prin_vecs, i);
5112 prev_ortho_vec *= 1.0/prev_ortho_vec.normFrobenius();
5113
5114 for(int j = i + 1; j < num_prin_direc; ++j)
5115 {
5116 RealVector::scalarType proj_val = getCol(Teuchos::View, prin_vecs,
5117 j).dot(prev_ortho_vec);
5118
5119 RealVector prev_ortho_vec_copy(Teuchos::Copy,
5120 prev_ortho_vec.values(),
5121 prev_ortho_vec.length());
5122 prev_ortho_vec_copy *= proj_val;
5123
5124 getCol(Teuchos::View, prin_vecs, j) -= prev_ortho_vec_copy;
5125 }
5126
5127 q_mat.multiply(Teuchos::NO_TRANS, Teuchos::TRANS, eigenvals[i],
5128 prev_ortho_vec, prev_ortho_vec, 1.0);
5129 }
5130
5131 initialized = true;
5132 }
5133
5134 RealMatrix quad_prod(numVars, 1);
5135 quad_prod.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, q_mat, xC,
5136 0.0);
5137
5138 if (directFnASV[0] & 1) // **** f:
5139 {
5140 fnVals[0] = getCol(Teuchos::View, quad_prod, 0).dot(xC);
5141 }
5142
5143 if (directFnASV[0] & 2) // **** gradf
5144 {
5145 quad_prod *= 2.0;
5146 fnGrads = quad_prod;
5147 }
5148
5149 return 0; // no failure
5150 }
5151
5152 } // namespace Dakota
5153