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: ScalingModel
10 //- Description: Implementation code for the ScalingModel class
11 //- Owner: Brian Adams
12 //- Checked by:
13
14 #include "ScalingModel.hpp"
15
16 static const char rcsId[]="@(#) $Id$";
17
18 namespace Dakota {
19
20 // ---------
21 // Scaling-related constants
22 // ---------
23
24 /// minimum value allowed for a characteristic value when scaling; ten
25 /// orders of magnitude greater than DBL_MIN
26 const double SCALING_MIN_SCALE = 1.0e10*DBL_MIN;
27 /// lower bound on domain of logarithm function when scaling
28 const double SCALING_MIN_LOG = SCALING_MIN_SCALE;
29 /// logarithm base to be used when scaling
30 const double SCALING_LOGBASE = 10.0;
31 /// ln(SCALING_LOGBASE); needed in transforming variables in several places
32 const double SCALING_LN_LOGBASE = std::log(SCALING_LOGBASE);
33
34 /// to indicate type of object being scaled
35 enum { CDV, LINEAR, NONLIN, FN_LSQ };
36 /// to restrict type of auto scaling allowed
37 enum { DISALLOW, TARGET, BOUNDS };
38
39
40 /// initialization of static needed by RecastModel
41 ScalingModel* ScalingModel::scaleModelInstance(NULL);
42
43
44 /** This constructor computes various indices and mappings, then
45 updates the properties of the RecastModel */
46 ScalingModel::
ScalingModel(Model & sub_model)47 ScalingModel(Model& sub_model):
48 // BMA TODO: should the BitArrays be empty or same as submodel?
49 // recast_secondary_offset is the index to the equality constraints within
50 // the secondary responses
51 RecastModel(sub_model, SizetArray(), BitArray(), BitArray(),
52 sub_model.num_primary_fns(), sub_model.num_secondary_fns(),
53 sub_model.num_nonlinear_ineq_constraints(),
54 response_order(sub_model))
55 {
56 if (outputLevel >= DEBUG_OUTPUT)
57 Cout << "Initializing scaling transformation" << std::endl;
58
59 modelId = RecastModel::recast_model_id(root_model_id(), "SCALING");
60 // RecastModel is constructed, then later initialized because scaled
61 // properties need to be set on the RecastModel, like bounds, but
62 // the nonlinearity of the mapping is determined by the scales
63 // themselves.
64
65 // initialize_scaling function needs to modify the iteratedModel
66 // compute needed class data...
67 initialize_scaling(sub_model);
68
69
70 // No change in sizes for scaling
71 size_t num_primary = sub_model.num_primary_fns(),
72 num_secondary = sub_model.num_secondary_fns(),
73 num_recast_fns = num_primary + num_secondary;
74
75 // the scaling transformation doesn't change any counts of variables
76 // or responses, but may require a nonlinear transformation of
77 // responses (grad, Hess) when variables are transformed
78
79 // ---
80 // Variables mapping (one-to-one)
81 // ---
82
83 // BMA TODO: Scaling only works with continuous variables, but
84 // verify this is the correct way to default map the other
85 // variables.
86
87 // We assume the mapping is for all active variables
88 size_t total_active_vars =
89 sub_model.cv() + sub_model.div() + sub_model.dsv() + sub_model.drv();
90 Sizet2DArray vars_map_indices(total_active_vars);
91 bool nonlinear_vars_mapping = false;
92 for (size_t i=0; i<total_active_vars; ++i) {
93 vars_map_indices[i].resize(1);
94 vars_map_indices[i][0] = i;
95 if (varsScaleFlag && cvScaleTypes[i] & SCALE_LOG)
96 nonlinear_vars_mapping = true;
97 }
98
99
100 // mappings from the submodel Response to residual Response
101 BoolDequeArray nonlinear_resp_mapping(num_recast_fns);
102
103 // ---
104 // Primary mapping
105 // ---
106 Sizet2DArray primary_resp_map_indices(num_primary);
107 for (size_t i=0; i<num_primary; i++) {
108 primary_resp_map_indices[i].resize(1);
109 primary_resp_map_indices[i][0] = i;
110 nonlinear_resp_mapping[i].resize(1);
111 nonlinear_resp_mapping[i][0] =
112 (primaryRespScaleFlag && responseScaleTypes[i] & SCALE_LOG);
113 }
114
115 // ---
116 // Secondary mapping (one-to-one)
117 // ---
118 Sizet2DArray secondary_resp_map_indices(num_secondary);
119 for (size_t i=0; i<num_secondary; i++) {
120 secondary_resp_map_indices[i].resize(1);
121 secondary_resp_map_indices[i][0] = num_primary + i;
122 nonlinear_resp_mapping[num_primary+i].resize(1);
123 nonlinear_resp_mapping[num_primary+i][0] = secondaryRespScaleFlag &&
124 responseScaleTypes[num_primary + i] & SCALE_LOG;
125 }
126
127 // callbacks for RecastModel transformations: default maps when not needed
128 void (*variables_map) (const Variables&, Variables&) = variables_scaler;
129 void (*set_map) (const Variables&, const ActiveSet&, ActiveSet&) = NULL;
130 // register primary response scaler if requested, or variables scaled
131 void (*primary_resp_map)
132 (const Variables&, const Variables&, const Response&, Response&) =
133 (primaryRespScaleFlag || varsScaleFlag) ? primary_resp_scaler : NULL;
134 // scale secondary response if requested, or variables scaled
135 void (*secondary_resp_map) (const Variables&, const Variables&,
136 const Response&, Response&) =
137 (secondaryRespScaleFlag || varsScaleFlag) ? secondary_resp_scaler : NULL;
138
139 RecastModel::
140 init_maps(vars_map_indices, nonlinear_vars_mapping, variables_map, set_map,
141 primary_resp_map_indices, secondary_resp_map_indices,
142 nonlinear_resp_mapping, primary_resp_map, secondary_resp_map);
143
144 // need inverse vars mapping for use with late updates from sub-model
145 inverse_mappings(variables_unscaler, NULL, NULL, NULL);
146
147 // Preserve weights through scaling transformation
148 primary_response_fn_weights(sub_model.primary_response_fn_weights());
149
150 // Preserve sense through scaling transformation
151 // Note: for a specification of negative scaling, we will assume that
152 // the user's intent is to overlay the scaling and sense as specified,
153 // such that we will not enforce a flip in sense for negative scaling.
154 primary_response_fn_sense(sub_model.primary_response_fn_sense());
155
156 // BMA TODO: consume scales so they aren't here anymore?
157 }
158
159
~ScalingModel()160 ScalingModel::~ScalingModel()
161 { /* empty dtor */}
162
163
164 /** Since this convenience function is public, it must have a
165 fall-through to return a copy for when this scaling type isn't
166 active. */
cv_scaled2native(const RealVector & scaled_cv) const167 RealVector ScalingModel::cv_scaled2native(const RealVector& scaled_cv) const
168 {
169 return (varsScaleFlag) ?
170 modify_s2n(scaled_cv, cvScaleTypes, cvScaleMultipliers, cvScaleOffsets) :
171 scaled_cv;
172 }
173
174
175 /** Since this convenience function is public, it must behave
176 correctly when this scale type isn't active. It does, because it
177 modifies in-place */
resp_scaled2native(const Variables & native_vars,Response & updated_resp) const178 void ScalingModel::resp_scaled2native(const Variables& native_vars,
179 Response& updated_resp) const
180 {
181 if (primaryRespScaleFlag || secondaryRespScaleFlag ||
182 // NOTE: formerly DakotaLeastSq::post_run didn't have this check:
183 need_resp_trans_byvars(updated_resp.active_set_request_vector(), 0,
184 num_primary_fns())){
185 size_t num_nln_cons =
186 num_nonlinear_ineq_constraints() + num_nonlinear_eq_constraints();
187 Response tmp_response = updated_resp.copy();
188 if (primaryRespScaleFlag ||
189 need_resp_trans_byvars(tmp_response.active_set_request_vector(), 0,
190 num_primary_fns())) {
191 response_modify_s2n(native_vars, updated_resp, tmp_response, 0,
192 num_primary_fns());
193 updated_resp.update_partial(0, num_primary_fns(), tmp_response, 0 );
194 }
195 if (secondaryRespScaleFlag ||
196 need_resp_trans_byvars(tmp_response.active_set_request_vector(),
197 num_primary_fns(), num_nln_cons)) {
198 response_modify_s2n(native_vars, updated_resp, tmp_response,
199 num_primary_fns(), num_nln_cons);
200 updated_resp.update_partial(num_primary_fns(), num_nln_cons,
201 tmp_response, num_primary_fns());
202 }
203 }
204 }
205
206
207 /** Since this convenience function is public, it must have a
208 fall-through to return a copy for when this scaling type isn't
209 active. */
210 void ScalingModel::
secondary_resp_scaled2native(const RealVector & scaled_nln_cons,const ShortArray & asv,RealVector & native_fns) const211 secondary_resp_scaled2native(const RealVector& scaled_nln_cons,
212 const ShortArray& asv,
213 RealVector& native_fns) const
214 {
215 size_t num_nln_cons =
216 num_nonlinear_ineq_constraints() + num_nonlinear_eq_constraints();
217 if (secondaryRespScaleFlag ||
218 need_resp_trans_byvars(asv, num_primary_fns(), num_nln_cons)) {
219 // scale all functions, but only copy constraints
220 copy_data_partial
221 (modify_s2n(scaled_nln_cons, responseScaleTypes, responseScaleMultipliers,
222 responseScaleOffsets),
223 num_primary_fns(), num_nln_cons, native_fns, num_primary_fns());
224 }
225 else
226 copy_data_partial(scaled_nln_cons, num_primary_fns(), num_nln_cons,
227 native_fns, num_primary_fns());
228 }
229
230
231 /** Initialize scaling types, multipliers, and offsets. Update the
232 iteratedModel appropriately */
initialize_scaling(Model & sub_model)233 void ScalingModel::initialize_scaling(Model& sub_model)
234 {
235 if (outputLevel > NORMAL_OUTPUT)
236 Cout << "\nScalingModel: Scaling enabled ('auto' scaling is reported as "
237 << "derived values)" << std::endl;
238 else if (outputLevel > SILENT_OUTPUT)
239 Cout << "\nScalingModel: Scaling enabled" << std::endl;
240
241 // in the scaled case, perform numerical derivatives at the RecastModel level
242 // (override the RecastModel default and the subModel default)
243 supportsEstimDerivs = true;
244 sub_model.supports_derivative_estimation(false);
245
246 size_t num_cv = cv(), num_primary = num_primary_fns(),
247 num_nln_ineq = num_nonlinear_ineq_constraints(),
248 num_nln_eq = num_nonlinear_eq_constraints(),
249 num_lin_ineq = num_linear_ineq_constraints(),
250 num_lin_eq = num_linear_eq_constraints();
251
252 // temporary arrays
253 UShortArray tmp_types;
254 RealVector tmp_multipliers, tmp_offsets;
255 RealVector lbs, ubs, targets;
256 //RealMatrix linear_constraint_coeffs;
257
258 // NOTE: When retrieving scaling vectors from database, excepting linear
259 // constraints, they've already been checked at input to have length 0,
260 // 1, or number of vars, constraints, etc.
261
262 // -----------------
263 // CONTINUOUS DESIGN
264 // -----------------
265 const UShortArray& cdv_spec_types = scalingOpts.cvScaleTypes;
266 const RealVector& cdv_scales = scalingOpts.cvScales;
267 varsScaleFlag = scaling_active(cdv_spec_types);
268
269 copy_data(sub_model.continuous_lower_bounds(), lbs); // view->copy
270 copy_data(sub_model.continuous_upper_bounds(), ubs); // view->copy
271
272
273 compute_scaling(CDV, BOUNDS, num_cv, lbs, ubs, targets,
274 cdv_spec_types, cdv_scales, cvScaleTypes,
275 cvScaleMultipliers, cvScaleOffsets);
276
277 continuous_lower_bounds(lbs);
278 continuous_upper_bounds(ubs);
279 continuous_variables(
280 modify_n2s(sub_model.continuous_variables(), cvScaleTypes,
281 cvScaleMultipliers, cvScaleOffsets) );
282
283 if (outputLevel > NORMAL_OUTPUT && varsScaleFlag) {
284 StringArray cv_labels;
285 copy_data(continuous_variable_labels(), cv_labels);
286 print_scaling("Continuous design variable scales", cvScaleTypes,
287 cvScaleMultipliers, cvScaleOffsets, cv_labels);
288 }
289
290 // each responseScale* = [fnScale*, nonlinearIneqScale*, nonlinearEqScale*]
291 // to make transformations faster at run time
292 // numFns should reflect size of user-space model
293 responseScaleTypes.resize(numFns);
294 responseScaleMultipliers.resize(numFns);
295 responseScaleOffsets.resize(numFns);
296
297 // -------------------------
298 // OBJECTIVE FNS / LSQ TERMS
299 // -------------------------
300 const UShortArray& primary_spec_types = scalingOpts.priScaleTypes;
301 const RealVector& primary_scales = scalingOpts.priScales;
302 primaryRespScaleFlag = scaling_active(primary_spec_types);
303
304 lbs.size(0); ubs.size(0);
305 compute_scaling(FN_LSQ, DISALLOW, num_primary, lbs, ubs, targets,
306 primary_spec_types, primary_scales, tmp_types,
307 tmp_multipliers, tmp_offsets);
308
309 for (int i=0; i<num_primary; ++i) {
310 responseScaleTypes[i] = tmp_types[i];
311 responseScaleMultipliers[i] = tmp_multipliers[i];
312 responseScaleOffsets[i] = 0;
313 }
314
315 // --------------------
316 // NONLINEAR INEQUALITY
317 // --------------------
318 const UShortArray& nln_ineq_spec_types = scalingOpts.nlnIneqScaleTypes;
319 const RealVector& nln_ineq_scales = scalingOpts.nlnIneqScales;
320 secondaryRespScaleFlag = scaling_active(nln_ineq_spec_types);
321
322 lbs = sub_model.nonlinear_ineq_constraint_lower_bounds();
323 ubs = sub_model.nonlinear_ineq_constraint_upper_bounds();
324
325 compute_scaling(NONLIN, BOUNDS, num_nln_ineq, lbs, ubs,
326 targets, nln_ineq_spec_types, nln_ineq_scales, tmp_types,
327 tmp_multipliers, tmp_offsets);
328
329 for (int i=0; i<num_nln_ineq; ++i) {
330 responseScaleTypes[num_primary+i] = tmp_types[i];
331 responseScaleMultipliers[num_primary+i] = tmp_multipliers[i];
332 responseScaleOffsets[num_primary+i] = tmp_offsets[i];
333 }
334
335 nonlinear_ineq_constraint_lower_bounds(lbs);
336 nonlinear_ineq_constraint_upper_bounds(ubs);
337
338 // --------------------
339 // NONLINEAR EQUALITY
340 // --------------------
341 const UShortArray& nln_eq_spec_types = scalingOpts.nlnEqScaleTypes;
342 const RealVector& nln_eq_scales = scalingOpts.nlnEqScales;
343 secondaryRespScaleFlag
344 = (secondaryRespScaleFlag || scaling_active(nln_eq_spec_types));
345
346 lbs.size(0); ubs.size(0);
347 targets = sub_model.nonlinear_eq_constraint_targets();
348 compute_scaling(NONLIN, TARGET, num_nln_eq,
349 lbs, ubs, targets, nln_eq_spec_types, nln_eq_scales,
350 tmp_types, tmp_multipliers, tmp_offsets);
351
352 // BMA TODO: use copy_data?
353 for (int i=0; i<num_nln_eq; ++i) {
354 responseScaleTypes[num_primary+num_nln_ineq+i]
355 = tmp_types[i];
356 responseScaleMultipliers[num_primary+num_nln_ineq+i]
357 = tmp_multipliers[i];
358 responseScaleOffsets[num_primary+num_nln_ineq+i]
359 = tmp_offsets[i];
360 }
361
362 nonlinear_eq_constraint_targets(targets);
363
364 if (outputLevel > NORMAL_OUTPUT &&
365 (primaryRespScaleFlag || secondaryRespScaleFlag) )
366 print_scaling("Response scales", responseScaleTypes,
367 responseScaleMultipliers, responseScaleOffsets,
368 sub_model.response_labels());
369
370 // LINEAR CONSTRAINT SCALING:
371 // computed scales account for scaling in CVs x
372 // updating constraint coefficient matrices now handled in derived classes
373 // NOTE: cannot use offsets since would be an affine scale
374 // ScaleOffsets in this case are only for applying to bounds
375
376 // -----------------
377 // LINEAR INEQUALITY
378 // -----------------
379 const UShortArray& lin_ineq_spec_types = scalingOpts.linIneqScaleTypes;
380 const RealVector& lin_ineq_scales = scalingOpts.linIneqScales;
381
382 if ( ( lin_ineq_spec_types.size() != 0 &&
383 lin_ineq_spec_types.size() != 1 &&
384 lin_ineq_spec_types.size() != num_lin_ineq ) ||
385 ( lin_ineq_scales.length() != 0 && lin_ineq_scales.length() != 1 &&
386 lin_ineq_scales.length() != num_lin_ineq ) ) {
387 Cerr << "Error: linear_inequality_scale specifications must have length 0, "
388 << "1, or " << num_lin_ineq << ".\n";
389 abort_handler(-1);
390 }
391
392 linearIneqScaleOffsets.resize(num_lin_ineq);
393
394 lbs = sub_model.linear_ineq_constraint_lower_bounds();
395 ubs = sub_model.linear_ineq_constraint_upper_bounds();
396 targets.size(0);
397
398 const RealMatrix& lin_ineq_coeffs
399 = sub_model.linear_ineq_constraint_coeffs();
400 for (int i=0; i<num_lin_ineq; ++i) {
401
402 // compute A_i*cvScaleOffset for current constraint -- discrete variables
403 // aren't scaled so don't contribute
404 linearIneqScaleOffsets[i] = 0.0;
405 for (int j=0; j<num_cv; ++j)
406 linearIneqScaleOffsets[i] += lin_ineq_coeffs(i,j)*cvScaleOffsets[j];
407
408 lbs[i] -= linearIneqScaleOffsets[i];
409 ubs[i] -= linearIneqScaleOffsets[i];
410
411 }
412 compute_scaling(LINEAR, BOUNDS, num_lin_ineq,
413 lbs, ubs, targets, lin_ineq_spec_types, lin_ineq_scales,
414 linearIneqScaleTypes, linearIneqScaleMultipliers,
415 tmp_offsets);
416
417 linear_ineq_constraint_lower_bounds(lbs);
418 linear_ineq_constraint_upper_bounds(ubs);
419 linear_ineq_constraint_coeffs(
420 lin_coeffs_modify_n2s(lin_ineq_coeffs, cvScaleMultipliers,
421 linearIneqScaleMultipliers) );
422
423 if (outputLevel > NORMAL_OUTPUT && num_lin_ineq > 0)
424 print_scaling("Linear inequality scales (incl. any variable scaling)",
425 linearIneqScaleTypes, linearIneqScaleMultipliers,
426 linearIneqScaleOffsets, StringArray());
427
428 // ---------------
429 // LINEAR EQUALITY
430 // ---------------
431 const UShortArray& lin_eq_spec_types = scalingOpts.linEqScaleTypes;
432 const RealVector& lin_eq_scales = scalingOpts.linEqScales;
433
434 if ( ( lin_eq_spec_types.size() != 0 && lin_eq_spec_types.size() != 1 &&
435 lin_eq_spec_types.size() != num_lin_eq ) ||
436 ( lin_eq_scales.length() != 0 && lin_eq_scales.length() != 1 &&
437 lin_eq_scales.length() != num_lin_eq ) ) {
438 Cerr << "Error: linear_equality_scale specifications must have length 0, "
439 << "1, or " << num_lin_eq << ".\n";
440 abort_handler(-1);
441 }
442
443 linearEqScaleOffsets.resize(num_lin_eq);
444
445 lbs.size(0); ubs.size(0);
446 targets = sub_model.linear_eq_constraint_targets();
447
448 const RealMatrix& lin_eq_coeffs
449 = sub_model.linear_eq_constraint_coeffs();
450 for (int i=0; i<num_lin_eq; ++i) {
451 // compute A_i*cvScaleOffset for current constraint
452 linearEqScaleOffsets[i] = 0.0;
453 for (int j=0; j<num_cv; ++j)
454 linearEqScaleOffsets[i] += lin_eq_coeffs(i,j)*cvScaleOffsets[j];
455
456 targets[i] -= linearEqScaleOffsets[i];
457 }
458 compute_scaling(LINEAR, TARGET, num_lin_eq,
459 lbs, ubs, targets, lin_eq_spec_types, lin_eq_scales,
460 linearEqScaleTypes, linearEqScaleMultipliers,
461 tmp_offsets);
462
463 linear_eq_constraint_targets(targets);
464 linear_eq_constraint_coeffs(
465 lin_coeffs_modify_n2s(lin_eq_coeffs, cvScaleMultipliers,
466 linearEqScaleMultipliers) );
467
468 if (outputLevel > NORMAL_OUTPUT && num_lin_eq > 0)
469 print_scaling("Linear equality scales (incl. any variable scaling)",
470 linearEqScaleTypes, linearEqScaleMultipliers,
471 linearEqScaleOffsets, StringArray());
472
473 if (outputLevel > NORMAL_OUTPUT)
474 Cout << std::endl;
475 }
476
477
scaling_active(const UShortArray & scale_types)478 bool ScalingModel::scaling_active(const UShortArray& scale_types)
479 {
480 for (const auto& sc_type : scale_types) {
481 if (sc_type > SCALE_NONE)
482 return true;
483 }
484 return false; // false if array empty or all types == SCALE_NONE
485 }
486
487
488 // compute_scaling will potentially modify lbs, ubs, and targets; will resize
489 // and set class data referenced by scale_types, scale_mults, and scale_offsets
490 void ScalingModel::
compute_scaling(int object_type,int auto_type,int num_vars,RealVector & lbs,RealVector & ubs,RealVector & targets,const UShortArray & spec_types,const RealVector & scales,UShortArray & scale_types,RealVector & scale_mults,RealVector & scale_offsets)491 compute_scaling(int object_type, // type of object being scaled
492 int auto_type, // option for auto scaling type
493 int num_vars, // length of object being scaled
494 RealVector& lbs, RealVector& ubs,
495 RealVector& targets, const UShortArray& spec_types,
496 const RealVector& scales, UShortArray& scale_types,
497 RealVector& scale_mults, RealVector& scale_offsets)
498 {
499 // temporary arrays
500 unsigned short tmp_scl_type;
501 Real tmp_bound, tmp_mult, tmp_offset;
502 //RealMatrix linear_constraint_coeffs;
503
504 const int num_scale_types = spec_types.size();
505 const int num_scales = scales.length();
506
507 scale_types.resize(num_vars);
508 scale_mults.resize(num_vars);
509 scale_offsets.resize(num_vars);
510
511 for (int i=0; i<num_vars; ++i) {
512
513 //set defaults
514 scale_types[i] = SCALE_NONE;
515 scale_mults[i] = 1.0;
516 scale_offsets[i] = 0.0;
517
518 // set the string for scale_type, depending on whether user sent
519 // 0, 1, or N scale_strings
520 tmp_scl_type = SCALE_NONE;
521 if (num_scale_types == 1)
522 tmp_scl_type = spec_types[0];
523 else if (num_scale_types > 1)
524 tmp_scl_type = spec_types[i];
525
526 if (tmp_scl_type > SCALE_NONE) {
527 size_t num_lin_cons =
528 num_linear_ineq_constraints() + num_linear_eq_constraints();
529 if (object_type == CDV && num_lin_cons > 0 && tmp_scl_type == SCALE_LOG) {
530 Cerr << "Error: Continuous design variables cannot be logarithmically "
531 << "scaled when linear\nconstraints are present.\n";
532 abort_handler(-1);
533 }
534 else if ( num_scales > 0 ) {
535
536 // process scale_values for all types of scaling
537 // indicate that scale values are active, update bounds, poss. negating
538 scale_types[i] |= SCALE_VALUE;
539 scale_mults[i] = (num_scales == 1) ? scales[0] : scales[i];
540 if (std::fabs(scale_mults[i]) < SCALING_MIN_SCALE)
541 Cout << "Warning: abs(scale) < " << SCALING_MIN_SCALE
542 << " provided; carefully verify results.\n";
543 // adjust bounds or targets
544 if (!lbs.empty()) {
545 // don't scale bounds if the user intended no bound
546 if (-BIG_REAL_BOUND < lbs[i])
547 lbs[i] /= scale_mults[i];
548 if (ubs[i] < BIG_REAL_BOUND)
549 ubs[i] /= scale_mults[i];
550 if (scale_mults[i] < 0) {
551 tmp_bound = lbs[i];
552 lbs[i] = ubs[i];
553 ubs[i] = tmp_bound;
554 }
555 }
556 else if (!targets.empty())
557 targets[i] /= scale_mults[i];
558 }
559
560 } // endif for generic scaling preprocessing
561
562 // At this point bounds/targets are scaled with user-provided values and
563 // scale_mults are set to user-provided values.
564 // Now auto or log scale as relevant and allowed:
565 if ( tmp_scl_type == SCALE_AUTO && auto_type > DISALLOW ) {
566 bool scale_flag = false; // will be true for valid auto-scaling
567 if ( auto_type == TARGET ) {
568 scale_flag = compute_scale_factor(targets[i], &tmp_mult);
569 tmp_offset = 0.0;
570 }
571 else if (auto_type == BOUNDS )
572 scale_flag = compute_scale_factor(lbs[i], ubs[i],
573 &tmp_mult, &tmp_offset);
574 if (scale_flag) {
575
576 scale_types[i] |= SCALE_VALUE;
577 // tmp_offset was calculated based on scaled bounds, so
578 // includes the effect of user scale values, so in computing
579 // the offset, need to include the effect of any user-supplied
580 // characteristic value scaling, then update multipliers
581 scale_offsets[i] += tmp_offset*scale_mults[i];
582 scale_mults[i] *= tmp_mult;
583
584 // necessary since the initial values may have already been value scaled
585 if (auto_type == BOUNDS) {
586 // don't scale bounds if the user intended no bound
587 if (-BIG_REAL_BOUND < lbs[i])
588 lbs[i] = (lbs[i] - tmp_offset)/tmp_mult;
589 if (ubs[i] < BIG_REAL_BOUND)
590 ubs[i] = (ubs[i] - tmp_offset)/tmp_mult;
591 }
592 else if (auto_type == TARGET)
593 targets[i] /= tmp_mult;
594
595 }
596 }
597 else if ( tmp_scl_type == SCALE_LOG ) {
598
599 scale_types[i] |= SCALE_LOG;
600 if (auto_type == BOUNDS) {
601 if (-BIG_REAL_BOUND < lbs[i]) {
602 if ( lbs[i] < SCALING_MIN_LOG )
603 Cout << "Warning: scale_type 'log' used without positive lower "
604 << "bound.\n";
605 lbs[i] = std::log(lbs[i])/SCALING_LN_LOGBASE;
606 }
607 if (ubs[i] < BIG_REAL_BOUND) {
608 if ( ubs[i] < SCALING_MIN_LOG )
609 Cout << "Warning: scale_type 'log' used without positive upper "
610 << "bound.\n";
611 ubs[i] = std::log(ubs[i])/SCALING_LN_LOGBASE;
612 }
613 }
614 else if (auto_type == TARGET) {
615 targets[i] = std::log(targets[i])/SCALING_LN_LOGBASE;
616 if ( targets[i] < SCALING_MIN_LOG )
617 Cout << "Warning: scale_type 'log' used without positive target.\n";
618 }
619 }
620
621 } // end for each variable
622 }
623
624
625 // automatically compute scaling factor
626 // bounds case allows for negative multipliers
627 // returns true if a valid scaling factor was computed
compute_scale_factor(const Real lower_bound,const Real upper_bound,Real * multiplier,Real * offset)628 bool ScalingModel::compute_scale_factor(const Real lower_bound,
629 const Real upper_bound,
630 Real *multiplier, Real *offset)
631 {
632 /* Compute scaleMultipliers for each design var, fn, constr, etc.
633 1. user-specified scaling was already detected at higher level
634 2. check for two-sided bounds
635 3. then check for one-sided bounds
636 4. else resort to no scaling
637
638 Auto-scaling is to [0,1] (affine scaling) in two sided case
639 or value of bound in one-sided case
640 */
641
642 bool lb_flag = false;
643 bool ub_flag = false;
644
645 if (-BIG_REAL_BOUND < lower_bound)
646 lb_flag = true;
647 if (upper_bound < BIG_REAL_BOUND)
648 ub_flag = true;
649
650 // process two-sided, then single-sided bounds
651 if ( lb_flag && ub_flag ) {
652 *multiplier = upper_bound - lower_bound;
653 *offset = lower_bound;
654 }
655 else if (lb_flag) {
656 *multiplier = lower_bound;
657 *offset = 0.0;
658 }
659 else if (ub_flag) {
660 *multiplier = upper_bound;
661 *offset = 0.0;
662 }
663 else {
664 Cout << "Warning: abs(bounds) > BIG_REAL_BOUND. Not auto-scaling "
665 << "component." << std::endl;
666 *multiplier = 1.0;
667 *offset = 0.0;
668 return(false);
669 }
670
671 if (std::fabs(*multiplier) < SCALING_MIN_SCALE) {
672 *multiplier = (*multiplier >= 0.0) ? SCALING_MIN_SCALE :
673 -(SCALING_MIN_SCALE);
674 Cout << "Warning: in auto-scaling abs(computed scale) < "
675 << SCALING_MIN_SCALE << "; resetting scale = "
676 << *multiplier << ".\n";
677 }
678
679 return(true);
680 }
681
682
683 // automatically compute scaling factor
684 // target case allows for negative multipliers
685 // returns true if a valid scaling factor was computed
compute_scale_factor(const Real target,Real * multiplier)686 bool ScalingModel::compute_scale_factor(const Real target, Real *multiplier)
687 {
688 if ( std::fabs(target) < BIG_REAL_BOUND )
689 *multiplier = target;
690 else {
691 Cout << "Automatic Scaling Warning: abs(target) > BIG_REAL_BOUND. "
692 << "Not scaling this component." << std::endl;
693 *multiplier = 1.0;
694 return(false);
695 }
696
697 if (std::fabs(*multiplier) < SCALING_MIN_SCALE) {
698 *multiplier = (*multiplier >= 0.0) ? SCALING_MIN_SCALE :
699 -(SCALING_MIN_SCALE);
700 Cout << "Warning: in auto-scaling abs(computed scale) < "
701 << SCALING_MIN_SCALE << "; resetting scale = "
702 << *multiplier << ".\n";
703 }
704
705 return(true);
706 }
707
708
709 /** compute scaled linear constraint matrix given design variable
710 multipliers and linear scaling multipliers. Only scales components
711 corresponding to continuous variables so for src_coeffs of size MxN,
712 lin_multipliers.size() <= M, cv_multipliers.size() <= N */
713 RealMatrix ScalingModel::
lin_coeffs_modify_n2s(const RealMatrix & src_coeffs,const RealVector & cv_multipliers,const RealVector & lin_multipliers) const714 lin_coeffs_modify_n2s(const RealMatrix& src_coeffs,
715 const RealVector& cv_multipliers,
716 const RealVector& lin_multipliers) const
717 {
718 RealMatrix dest_coeffs(src_coeffs);
719 for (int i=0; i<lin_multipliers.length(); ++i)
720 for (int j=0; j<cv_multipliers.length(); ++j)
721 dest_coeffs(i,j) *= cv_multipliers[j] / lin_multipliers[i];
722
723 return(dest_coeffs);
724 }
725
726
print_scaling(const String & info,const UShortArray & scale_types,const RealVector & scale_mults,const RealVector & scale_offsets,const StringArray & labels)727 void ScalingModel::print_scaling(const String& info, const UShortArray& scale_types,
728 const RealVector& scale_mults,
729 const RealVector& scale_offsets,
730 const StringArray& labels)
731 {
732 // labels will be empty for linear constraints
733 Cout << "\n" << info << ":\n";
734 Cout << "scale type " << std::setw(write_precision+7) << "multiplier" << " "
735 << std::setw(write_precision+7) << "offset"
736 << (labels.empty() ? " constraint number" : " label") << std::endl;
737 for (size_t i=0; i<scale_types.size(); ++i) {
738 switch (scale_types[i]) {
739 case SCALE_NONE:
740 Cout << "none ";
741 break;
742 case SCALE_VALUE:
743 Cout << "value ";
744 break;
745 case SCALE_LOG:
746 Cout << "log ";
747 break;
748 case (SCALE_VALUE | SCALE_LOG):
749 Cout << "value+log ";
750 break;
751 }
752 Cout << std::setw(write_precision+7) << scale_mults[i] << " "
753 << std::setw(write_precision+7) << scale_offsets[i] << " ";
754 if (labels.empty())
755 Cout << i << std::endl;
756 else
757 Cout << labels[i] << std::endl;
758 }
759 }
760
761
762 /** Variables map from iterator/scaled space to user/native space
763 using a RecastModel. */
764 void ScalingModel::
variables_scaler(const Variables & scaled_vars,Variables & native_vars)765 variables_scaler(const Variables& scaled_vars, Variables& native_vars)
766 {
767 if (scaleModelInstance->outputLevel > NORMAL_OUTPUT) {
768 Cout << "\n----------------------------------";
769 Cout << "\nPre-processing Function Evaluation";
770 Cout << "\n----------------------------------";
771 Cout << "\nVariables before unscaling transformation:\n";
772 write_data(Cout, scaled_vars.continuous_variables(),
773 scaled_vars.continuous_variable_labels());
774 Cout << std::endl;
775 }
776 native_vars.continuous_variables
777 (scaleModelInstance->modify_s2n(scaled_vars.continuous_variables(),
778 scaleModelInstance->cvScaleTypes,
779 scaleModelInstance->cvScaleMultipliers,
780 scaleModelInstance->cvScaleOffsets));
781 }
782
783 void ScalingModel::
variables_unscaler(const Variables & native_vars,Variables & scaled_vars)784 variables_unscaler(const Variables& native_vars, Variables& scaled_vars)
785 {
786 scaled_vars.continuous_variables
787 (scaleModelInstance->modify_n2s(native_vars.continuous_variables(),
788 scaleModelInstance->cvScaleTypes,
789 scaleModelInstance->cvScaleMultipliers,
790 scaleModelInstance->cvScaleOffsets));
791 }
792
793
794 void ScalingModel::
primary_resp_scaler(const Variables & native_vars,const Variables & scaled_vars,const Response & native_response,Response & iterator_response)795 primary_resp_scaler(const Variables& native_vars, const Variables& scaled_vars,
796 const Response& native_response, Response& iterator_response)
797 {
798 // scaling is always applied on a model with user's original size,
799 // so can query this object or the underlying model for sizes
800
801 // need to scale if primary responses are scaled or (variables are
802 // scaled and grad or hess requested)
803 size_t start_offset = 0;
804 size_t num_responses = scaleModelInstance->num_primary_fns();
805 bool scale_transform_needed =
806 scaleModelInstance->primaryRespScaleFlag ||
807 scaleModelInstance->need_resp_trans_byvars
808 (native_response.active_set_request_vector(), start_offset, num_responses);
809
810 if (scale_transform_needed) {
811 if (scaleModelInstance->outputLevel > NORMAL_OUTPUT) {
812 Cout << "\n--------------------------------------------";
813 Cout << "\nPost-processing Function Evaluation: Primary";
814 Cout << "\n--------------------------------------------" << std::endl;
815 }
816 scaleModelInstance->
817 response_modify_n2s(native_vars, native_response,
818 iterator_response, start_offset, num_responses);
819
820 }
821 else
822 // could reach this if variables are scaled and only functions are requested
823 iterator_response.update_partial(start_offset, num_responses,
824 native_response, start_offset);
825 }
826
827
828 /** Constraint function map from user/native space to iterator/scaled/combined
829 space using a RecastModel. */
830 void ScalingModel::
secondary_resp_scaler(const Variables & native_vars,const Variables & scaled_vars,const Response & native_response,Response & iterator_response)831 secondary_resp_scaler(const Variables& native_vars,
832 const Variables& scaled_vars,
833 const Response& native_response,
834 Response& iterator_response)
835 {
836 // scaling is always applied on a model with user's original size,
837 // so can query this object or the underlying model for sizes
838
839 // need to scale if secondary responses are scaled or (variables are
840 // scaled and grad or hess requested)
841 size_t start_offset = scaleModelInstance->num_primary_fns();
842 size_t num_nln_cons = scaleModelInstance->num_nonlinear_ineq_constraints() +
843 scaleModelInstance->num_nonlinear_eq_constraints();
844 bool scale_transform_needed =
845 scaleModelInstance->secondaryRespScaleFlag ||
846 scaleModelInstance->need_resp_trans_byvars
847 (native_response.active_set_request_vector(), start_offset, num_nln_cons);
848
849 if (scale_transform_needed) {
850 if (scaleModelInstance->outputLevel > NORMAL_OUTPUT) {
851 Cout << "\n----------------------------------------------";
852 Cout << "\nPost-processing Function Evaluation: Secondary";
853 Cout << "\n----------------------------------------------" << std::endl;
854 }
855 scaleModelInstance->
856 response_modify_n2s(native_vars, native_response,
857 iterator_response, start_offset, num_nln_cons);
858 }
859 else
860 // could reach this if variables are scaled and only functions are requested
861 iterator_response.update_partial(start_offset, num_nln_cons,
862 native_response, start_offset);
863 }
864
865
866 /** Determine if variable transformations present and derivatives
867 requested, which implies a response transformation is necessay */
868 bool ScalingModel::
need_resp_trans_byvars(const ShortArray & asv,int start_index,int num_resp) const869 need_resp_trans_byvars(const ShortArray& asv, int start_index,
870 int num_resp) const
871 {
872 if (varsScaleFlag)
873 for (size_t i=start_index; i<start_index+num_resp; ++i)
874 if (asv[i] & 2 || asv[i] & 4)
875 return(true);
876
877 return(false);
878 }
879
880 /** general RealVector mapping from native to scaled variables; loosely,
881 in greatest generality:
882 scaled_var = log( (native_var - offset) / multiplier ) */
883 RealVector ScalingModel::
modify_n2s(const RealVector & native_vars,const UShortArray & scale_types,const RealVector & multipliers,const RealVector & offsets) const884 modify_n2s(const RealVector& native_vars, const UShortArray& scale_types,
885 const RealVector& multipliers, const RealVector& offsets) const
886 {
887 RealVector scaled_vars(native_vars.length(), false);
888 for (int i=0; i<native_vars.length(); ++i) {
889
890 if (scale_types[i] & SCALE_VALUE)
891 scaled_vars[i] = (native_vars[i] - offsets[i]) / multipliers[i];
892 else
893 scaled_vars[i] = native_vars[i];
894
895 if (scale_types[i] & SCALE_LOG)
896 scaled_vars[i] = std::log(scaled_vars[i])/SCALING_LN_LOGBASE;
897
898 }
899
900 return(scaled_vars);
901 }
902
903 /** general RealVector mapping from scaled to native variables and/or vals;
904 loosely, in greatest generality:
905 scaled_var = (LOG_BASE^scaled_var) * multiplier + offset */
906 RealVector ScalingModel::
modify_s2n(const RealVector & scaled_vars,const UShortArray & scale_types,const RealVector & multipliers,const RealVector & offsets) const907 modify_s2n(const RealVector& scaled_vars, const UShortArray& scale_types,
908 const RealVector& multipliers, const RealVector& offsets) const
909 {
910 RealVector native_vars(scaled_vars.length(), false);
911 for (RealVector::ordinalType i=0; i<scaled_vars.length(); ++i) {
912
913 if (scale_types[i] & SCALE_LOG)
914 native_vars[i] = std::pow(SCALING_LOGBASE, scaled_vars[i]);
915 else
916 native_vars[i] = scaled_vars[i];
917
918 if (scale_types[i] & SCALE_VALUE)
919 native_vars[i] = native_vars[i]*multipliers[i] + offsets[i];
920
921 }
922
923 return(native_vars);
924 }
925
926
927 /** Scaling response mapping: modifies response from a model
928 (user/native) for use in iterators (scaled). Maps num_responses
929 starting at response_offset */
response_modify_n2s(const Variables & native_vars,const Response & native_response,Response & recast_response,int start_offset,int num_responses) const930 void ScalingModel::response_modify_n2s(const Variables& native_vars,
931 const Response& native_response,
932 Response& recast_response,
933 int start_offset,
934 int num_responses) const
935 {
936 int i, j, k;
937 int end_offset = start_offset + num_responses;
938 SizetMultiArray var_ids;
939 RealVector cdv;
940
941 // unroll the unscaled (native/user-space) response
942 const ActiveSet& set = native_response.active_set();
943 const ShortArray& asv = set.request_vector();
944 const SizetArray& dvv = set.derivative_vector();
945 size_t num_deriv_vars = dvv.size();
946
947 if (dvv == native_vars.continuous_variable_ids()) {
948 var_ids.resize(boost::extents[native_vars.cv()]);
949 var_ids = native_vars.continuous_variable_ids();
950 cdv = native_vars.continuous_variables(); // view OK
951 }
952 else if (dvv == native_vars.inactive_continuous_variable_ids()) {
953 var_ids.resize(boost::extents[native_vars.icv()]);
954 var_ids = native_vars.inactive_continuous_variable_ids();
955 cdv = native_vars.inactive_continuous_variables(); // view OK
956 }
957 else { // general derivatives
958 var_ids.resize(boost::extents[native_vars.acv()]);
959 var_ids = native_vars.all_continuous_variable_ids();
960 cdv = native_vars.all_continuous_variables(); // view OK
961 }
962
963 if (outputLevel > NORMAL_OUTPUT) {
964 if (start_offset < num_primary_fns())
965 Cout << "Primary response after scaling transformation:\n";
966 else
967 Cout << "Secondary response after scaling transformation:\n";
968 }
969
970 // scale functions and constraints
971 // assumes Multipliers and Offsets are 1 and 0 when not in use
972 // there's a tradeoff here between flops and logic simplicity
973 // (responseScaleOffsets may be nonzero for constraints)
974 // iterate over components of ASV-requested functions and scale when necessary
975 Real recast_val;
976 const RealVector& native_vals = native_response.function_values();
977 const StringArray& recast_labels = recast_response.function_labels();
978 for (i = start_offset; i < end_offset; ++i)
979 if (asv[i] & 1) {
980 // SCALE_LOG case here includes case of SCALE_LOG && SCALE_VALUE
981 if (responseScaleTypes[i] & SCALE_LOG)
982 recast_val = std::log( (native_vals[i] - responseScaleOffsets[i]) /
983 responseScaleMultipliers[i] )/SCALING_LN_LOGBASE;
984 else if (responseScaleTypes[i] & SCALE_VALUE)
985 recast_val = (native_vals[i] - responseScaleOffsets[i]) /
986 responseScaleMultipliers[i];
987 else
988 recast_val = native_vals[i];
989 recast_response.function_value(recast_val, i);
990 if (outputLevel > NORMAL_OUTPUT)
991 Cout << " " << std::setw(write_precision+7)
992 << recast_val << ' ' << recast_labels[i] << '\n';
993 }
994
995 // scale gradients
996 const RealMatrix& native_grads = native_response.function_gradients();
997 for (i = start_offset; i < end_offset; ++i)
998 if (asv[i] & 2) {
999
1000 Real fn_divisor;
1001 if (responseScaleTypes[i] & SCALE_LOG)
1002 fn_divisor = (native_vals[i] - responseScaleOffsets[i]) *
1003 SCALING_LN_LOGBASE;
1004 else if (responseScaleTypes[i] & SCALE_VALUE)
1005 fn_divisor = responseScaleMultipliers[i];
1006 else
1007 fn_divisor = 1.;
1008
1009 RealVector recast_grad
1010 = recast_response.function_gradient_view(i);
1011 copy_data(native_grads[i], (int)num_deriv_vars, recast_grad);
1012 for (j=0; j<num_deriv_vars; ++j) {
1013 size_t xj_index = find_index(var_ids, dvv[j]);
1014
1015 // first multiply by d(f_scaled)/d(f) based on scaling type
1016 recast_grad[xj_index] /= fn_divisor;
1017
1018 // now multiply by d(x)/d(x_scaled)
1019 if (cvScaleTypes[xj_index] & SCALE_LOG)
1020 recast_grad[xj_index] *=
1021 (cdv[xj_index] - cvScaleOffsets[xj_index]) * SCALING_LN_LOGBASE;
1022 else if (cvScaleTypes[xj_index] & SCALE_VALUE)
1023 recast_grad[xj_index] *= cvScaleMultipliers[xj_index];
1024 }
1025 if (outputLevel > NORMAL_OUTPUT) {
1026 write_col_vector_trans(Cout, i, recast_response.function_gradients(),
1027 true, true, false);
1028 Cout << recast_labels[i] << " gradient\n";
1029 }
1030 }
1031
1032 // scale hessians
1033 const RealSymMatrixArray& native_hessians
1034 = native_response.function_hessians();
1035 for (i = start_offset; i < end_offset; ++i)
1036 if (asv[i] & 4) {
1037 RealSymMatrix recast_hess
1038 = recast_response.function_hessian_view(i);
1039 recast_hess.assign(native_hessians[i]);
1040
1041 Real offset_fn = 1.;
1042 if (responseScaleTypes[i] & SCALE_LOG)
1043 offset_fn = native_vals[i] - responseScaleOffsets[i];
1044 for (j=0; j<num_deriv_vars; ++j) {
1045 size_t xj_index = find_index(var_ids, dvv[j]);
1046 for (k=0; k<=j; ++k) {
1047 size_t xk_index = find_index(var_ids, dvv[k]);
1048
1049 // first multiply by d(f_scaled)/d(f) based on scaling type
1050 if (responseScaleTypes[i] & SCALE_LOG) {
1051
1052 recast_hess(xj_index,xk_index) -=
1053 native_grads(xj_index,i)*native_grads(xk_index,i) / offset_fn;
1054
1055 recast_hess(xj_index,xk_index) /= offset_fn*SCALING_LN_LOGBASE;
1056
1057 }
1058 else if (responseScaleTypes[i] & SCALE_VALUE)
1059 recast_hess(xj_index,xk_index) /= responseScaleMultipliers[i];
1060
1061 // now multiply by d(x)/d(x_scaled) for j,k
1062 if (cvScaleTypes[xj_index] & SCALE_LOG)
1063 recast_hess(xj_index,xk_index) *=
1064 (cdv[xj_index] - cvScaleOffsets[xj_index]) * SCALING_LN_LOGBASE;
1065 else if (cvScaleTypes[xj_index] & SCALE_VALUE)
1066 recast_hess(xj_index,xk_index) *= cvScaleMultipliers[xj_index];
1067
1068 if (cvScaleTypes[xk_index] & SCALE_LOG)
1069 recast_hess(xj_index,xk_index) *=
1070 (cdv[xk_index] - cvScaleOffsets[xk_index]) * SCALING_LN_LOGBASE;
1071 else if (cvScaleTypes[xk_index] & SCALE_VALUE)
1072 recast_hess(xj_index,xk_index) *= cvScaleMultipliers[xk_index];
1073
1074 // need gradient term only for diagonal entries
1075 if (xj_index == xk_index && cvScaleTypes[xj_index] & SCALE_LOG) {
1076 if (responseScaleTypes[i] & SCALE_LOG)
1077 recast_hess(xj_index,xk_index) += native_grads(xj_index,i)*
1078 (cdv[xj_index] - cvScaleOffsets[xj_index]) *
1079 SCALING_LN_LOGBASE / offset_fn;
1080 else
1081 recast_hess(xj_index,xk_index) += native_grads(xj_index,i)*
1082 (cdv[xj_index] - cvScaleOffsets[xj_index]) * SCALING_LN_LOGBASE
1083 * SCALING_LN_LOGBASE / responseScaleMultipliers[i];
1084 }
1085
1086 }
1087 }
1088 if (outputLevel > NORMAL_OUTPUT) {
1089 write_data(Cout, recast_hess, true, true, false);
1090 Cout << recast_labels[i] << " Hessian\n";
1091 }
1092 }
1093
1094 if (outputLevel > NORMAL_OUTPUT)
1095 Cout << std::endl;
1096 }
1097
1098 /** Unscaling response mapping: modifies response from scaled
1099 (iterator) to native (user) space. Maps num_responses starting at
1100 response_offset. If response_unscale = false, only variables will
1101 be unscaled, and responses left in scaled space. */
response_modify_s2n(const Variables & native_vars,const Response & scaled_response,Response & native_response,int start_offset,int num_responses,bool unscale_resp) const1102 void ScalingModel::response_modify_s2n(const Variables& native_vars,
1103 const Response& scaled_response,
1104 Response& native_response,
1105 int start_offset,
1106 int num_responses,
1107 bool unscale_resp) const
1108 {
1109 using std::pow;
1110
1111 int i, j, k;
1112 int end_offset = start_offset + num_responses;
1113 SizetMultiArray var_ids;
1114 RealVector cdv;
1115
1116 // unroll the unscaled (native/user-space) response
1117 const ActiveSet& set = scaled_response.active_set();
1118 const ShortArray& asv = set.request_vector();
1119 const SizetArray& dvv = set.derivative_vector();
1120 size_t num_deriv_vars = dvv.size();
1121
1122 if (dvv == native_vars.continuous_variable_ids()) {
1123 var_ids.resize(boost::extents[native_vars.cv()]);
1124 var_ids = native_vars.continuous_variable_ids();
1125 cdv = native_vars.continuous_variables(); // view OK
1126 }
1127 else if (dvv == native_vars.inactive_continuous_variable_ids()) {
1128 var_ids.resize(boost::extents[native_vars.icv()]);
1129 var_ids = native_vars.inactive_continuous_variable_ids();
1130 cdv = native_vars.inactive_continuous_variables(); // view OK
1131 }
1132 else { // general derivatives
1133 var_ids.resize(boost::extents[native_vars.acv()]);
1134 var_ids = native_vars.all_continuous_variable_ids();
1135 cdv = native_vars.all_continuous_variables(); // view OK
1136 }
1137
1138 if (outputLevel > NORMAL_OUTPUT) {
1139 if (start_offset < num_primary_fns())
1140 Cout << "Primary response after unscaling transformation:\n";
1141 else
1142 Cout << "Secondary response after unscaling transformation:\n";
1143 }
1144
1145 // scale functions and constraints
1146 // assumes Multipliers and Offsets are 1 and 0 when not in use
1147 // there's a tradeoff here between flops and logic simplicity
1148 // (responseScaleOffsets may be nonzero for constraints)
1149 // iterate over components of ASV-requested functions and scale when necessary
1150 Real native_val;
1151 const RealVector& scaled_vals = scaled_response.function_values();
1152 const StringArray& native_labels = native_response.function_labels();
1153 for (i = start_offset; i < end_offset; ++i)
1154 if (asv[i] & 1) {
1155 // SCALE_LOG case here includes case of SCALE_LOG && SCALE_VALUE
1156 if (unscale_resp && responseScaleTypes[i] & SCALE_LOG)
1157 native_val = pow(SCALING_LOGBASE, scaled_vals[i]) *
1158 responseScaleMultipliers[i] + responseScaleOffsets[i];
1159
1160 else if (unscale_resp && responseScaleTypes[i] & SCALE_VALUE)
1161 native_val = scaled_vals[i]*responseScaleMultipliers[i] +
1162 responseScaleOffsets[i];
1163 else
1164 native_val = scaled_vals[i];
1165 native_response.function_value(native_val, i);
1166 if (outputLevel > NORMAL_OUTPUT)
1167 Cout << " " << std::setw(write_precision+7)
1168 << native_val << ' ' << native_labels[i] << '\n';
1169 }
1170
1171 // scale gradients
1172 Real df_dfscl;
1173 const RealMatrix& scaled_grads = scaled_response.function_gradients();
1174 for (i = start_offset; i < end_offset; ++i)
1175 if (asv[i] & 2) {
1176
1177 if (unscale_resp && responseScaleTypes[i] & SCALE_LOG)
1178 df_dfscl = pow(SCALING_LOGBASE, scaled_vals[i]) * SCALING_LN_LOGBASE *
1179 responseScaleMultipliers[i];
1180 else if (unscale_resp && responseScaleTypes[i] & SCALE_VALUE)
1181 df_dfscl = responseScaleMultipliers[i];
1182 else
1183 df_dfscl = 1.;
1184
1185 RealVector native_grad
1186 = native_response.function_gradient_view(i);
1187 copy_data(scaled_grads[i], (int)num_deriv_vars, native_grad);
1188 for (j=0; j<num_deriv_vars; ++j) {
1189 size_t xj_index = find_index(var_ids, dvv[j]);
1190
1191 // first multiply by d(f)/d(f_scaled) based on scaling type
1192 native_grad[xj_index] *= df_dfscl;
1193
1194 // now multiply by d(x_scaled)/d(x)
1195 if (cvScaleTypes[xj_index] & SCALE_LOG)
1196 native_grad[xj_index] /= (cdv[xj_index] - cvScaleOffsets[xj_index]) *
1197 SCALING_LN_LOGBASE;
1198 else if (cvScaleTypes[xj_index] & SCALE_VALUE)
1199 native_grad[xj_index] /= cvScaleMultipliers[xj_index];
1200 }
1201 if (outputLevel > NORMAL_OUTPUT) {
1202 write_col_vector_trans(Cout, i, native_response.function_gradients(),
1203 true, true, false);
1204 Cout << native_labels[i] << " gradient\n";
1205 }
1206 }
1207
1208 // scale hessians
1209 const RealSymMatrixArray& scaled_hessians
1210 = scaled_response.function_hessians();
1211 for (i = start_offset; i < end_offset; ++i)
1212 if (asv[i] & 4) {
1213
1214 if (unscale_resp && responseScaleTypes[i] & SCALE_LOG)
1215 df_dfscl = pow(SCALING_LOGBASE, scaled_vals[i]) * SCALING_LN_LOGBASE *
1216 responseScaleMultipliers[i];
1217 else if (unscale_resp && responseScaleTypes[i] & SCALE_VALUE)
1218 df_dfscl = responseScaleMultipliers[i];
1219 else
1220 df_dfscl = 1.;
1221
1222 RealSymMatrix native_hess
1223 = native_response.function_hessian_view(i);
1224 native_hess.assign(scaled_hessians[i]);
1225 for (j=0; j<num_deriv_vars; ++j) {
1226 size_t xj_index = find_index(var_ids, dvv[j]);
1227 for (k=0; k<=j; ++k) {
1228 size_t xk_index = find_index(var_ids, dvv[k]);
1229
1230 // first multiply by d(f_scaled)/d(f) based on scaling type
1231 // have to skip the addend when no response scaling...
1232 if (unscale_resp && responseScaleTypes[i] & SCALE_LOG) {
1233
1234 native_hess(xj_index,xk_index) += scaled_grads(xj_index,i) *
1235 scaled_grads(xk_index,i) * SCALING_LN_LOGBASE;
1236
1237 native_hess(xj_index,xk_index) *= df_dfscl;
1238
1239 }
1240 else if (unscale_resp && responseScaleTypes[i] & SCALE_VALUE)
1241 native_hess(xj_index,xk_index) *= df_dfscl;
1242
1243 // now multiply by d(x_scaled)/d(x) for j,k
1244 if (cvScaleTypes[xj_index] & SCALE_LOG)
1245 native_hess(xj_index,xk_index) /=
1246 (cdv[xj_index] - cvScaleOffsets[xj_index]) * SCALING_LN_LOGBASE;
1247 else if (cvScaleTypes[xj_index] & SCALE_VALUE)
1248 native_hess(xj_index,xk_index) /= cvScaleMultipliers[xj_index];
1249
1250 if (cvScaleTypes[xk_index] & SCALE_LOG)
1251 native_hess(xj_index,xk_index) /=
1252 (cdv[xk_index] - cvScaleOffsets[xk_index]) * SCALING_LN_LOGBASE;
1253 else if (cvScaleTypes[xk_index] & SCALE_VALUE)
1254 native_hess(xj_index,xk_index) /= cvScaleMultipliers[xk_index];
1255
1256 // need gradient term only for diagonal entries
1257 if (xj_index == xk_index && cvScaleTypes[xj_index] & SCALE_LOG)
1258 native_hess(xj_index,xk_index) -=
1259 df_dfscl * scaled_grads(xj_index,i) /
1260 (cdv[xj_index] - cvScaleOffsets[xj_index]) *
1261 (cdv[xj_index] - cvScaleOffsets[xj_index]) * SCALING_LN_LOGBASE;
1262
1263 }
1264 }
1265 if (outputLevel > NORMAL_OUTPUT) {
1266 write_data(Cout, native_hess, true, true, false);
1267 Cout << native_labels[i] << " Hessian\n";
1268 }
1269 }
1270
1271 if (outputLevel > NORMAL_OUTPUT)
1272 Cout << std::endl;
1273 }
1274
default_active_set()1275 ActiveSet ScalingModel::default_active_set() {
1276 // A ScalingModel has the same number of responses as its
1277 // submodel. It is also assumed to have supportEstimDerivs == true
1278 ActiveSet set;
1279 set.derivative_vector(currentVariables.all_continuous_variable_ids());
1280 bool has_deriv_vars = set.derivative_vector().size() != 0;
1281 // The ScalingModel can return at least everything that the submodel can.
1282 ShortArray asv(subModel.default_active_set().request_vector());
1283
1284 // In addition, if mixed or numerical gradients are active, the ScalingModel
1285 // can return gradients for all responses
1286 if(has_deriv_vars) {
1287 if(gradientType != "none")
1288 for(auto &a : asv)
1289 a |= 2;
1290 // Also, if mixed, numerical, or quasi hessians are active, the ScalingModel
1291 // can return hessians for all responses
1292 if(hessianType != "none")
1293 for(auto &a : asv)
1294 a |= 4;
1295 }
1296 set.request_vector(asv);
1297 return set;
1298 }
1299
1300 } // namespace Dakota
1301