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: ProbabilityTransformModel
10 //- Description: Implementation code for the ProbabilityTransformModel class
11 //- Owner: Brian Adams
12 //- Checked by:
13
14 #include "ProbabilityTransformModel.hpp"
15 #include "MarginalsCorrDistribution.hpp"
16 #include "GammaRandomVariable.hpp"
17 #include "GumbelRandomVariable.hpp"
18 #include "FrechetRandomVariable.hpp"
19 #include "WeibullRandomVariable.hpp"
20
21 static const char rcsId[]="@(#) $Id$";
22
23 namespace Dakota
24 {
25
26 /// initialization of static needed by RecastModel
27 ProbabilityTransformModel* ProbabilityTransformModel::ptmInstance(NULL);
28
29
30 ProbabilityTransformModel::
ProbabilityTransformModel(const Model & x_model,short u_space_type,bool truncate_bnds,Real bnd)31 ProbabilityTransformModel(const Model& x_model, short u_space_type,
32 bool truncate_bnds, Real bnd) :
33 RecastModel(x_model), distParamDerivs(NO_DERIVS),
34 truncatedBounds(truncate_bnds), boundVal(bnd)
35 {
36 modelType = "probability_transform";
37 modelId = recast_model_id(root_model_id(), "PROBABILITY_TRANSFORM");
38
39 Sizet2DArray vars_map, primary_resp_map, secondary_resp_map;
40 SizetArray recast_vars_comps_total; // default: no change in cauv total
41 // we do not reorder the u-space variable types such that we preserve a
42 // 1-to-1 mapping with consistent ordering
43 Pecos::MultivariateDistribution& x_dist
44 = subModel.multivariate_distribution();
45 const BitArray& active_vars = x_dist.active_variables();
46 size_t i, num_active_rv = (active_vars.empty()) ?
47 x_dist.random_variables().size() : active_vars.count();
48 // *** TO DO: generalize vars_map across multiple active types
49 vars_map.resize(num_active_rv);
50 for (i=0; i<num_active_rv; ++i)
51 { vars_map[i].resize(1); vars_map[i][0] = i; }
52 primary_resp_map.resize(numFns);
53 for (i=0; i<numFns; ++i)
54 { primary_resp_map[i].resize(1); primary_resp_map[i][0] = i; }
55
56 // There is no additional response mapping beyond that required by the
57 // nonlinear variables mapping.
58 BoolDequeArray nonlinear_resp_map(numFns, BoolDeque(1, false));
59 const Response& x_resp = subModel.current_response();
60 const SharedVariablesData& svd = subModel.current_variables().shared_data();
61 const BitArray& all_relax_di = svd.all_relaxed_discrete_int();
62 const BitArray& all_relax_dr = svd.all_relaxed_discrete_real();
63 short recast_resp_order = 1; // recast resp order to be same as original resp
64 if (!x_resp.function_gradients().empty()) recast_resp_order |= 2;
65 if (!x_resp.function_hessians().empty()) recast_resp_order |= 4;
66
67 // initialize current{Variables,Response}, userDefinedConstraints
68 init_sizes(recast_vars_comps_total, all_relax_di, all_relax_dr, numFns,
69 0, 0, recast_resp_order);
70 // initialize invariant portions of probability transform within mvDist
71 // (requires currentVariables)
72 initialize_transformation(u_space_type);
73 // initialize Variables/Response/ActiveSet recastings (requires mvDist)
74 init_maps(vars_map, nonlinear_variables_mapping(x_dist, mvDist),
75 vars_u_to_x_mapping, set_u_to_x_mapping, primary_resp_map,
76 secondary_resp_map, nonlinear_resp_map, resp_x_to_u_mapping, NULL);
77 // publish inverse mappings for use in data imports. Since derivatives are
78 // not imported and response values are not transformed, an inverse variables
79 // transformation is sufficient for this purpose.
80 inverse_mappings(vars_x_to_u_mapping, NULL, NULL, NULL);
81 // initialize currentVariables based on subModel initial state
82 inverse_transform_variables(subModel.current_variables(), currentVariables);
83 }
84
85
~ProbabilityTransformModel()86 ProbabilityTransformModel::~ProbabilityTransformModel()
87 { }
88
89
initialize_dakota_variable_types()90 void ProbabilityTransformModel::initialize_dakota_variable_types()
91 {
92 // Note: ctor has called initialize_distribution_
93 // {transformation,types,correlations}(). Defining the transformation is
94 // deferred until Model::initialize_mapping() to allow for problem resizing.
95
96 const SharedVariablesData& x_svd = subModel.current_variables().shared_data();
97 bool cdv, ddv, cauv, dauv, ceuv, deuv, csv, dsv;
98 x_svd.active_subsets(cdv, ddv, cauv, dauv, ceuv, deuv, csv, dsv);
99 size_t i, num_cdv, num_ddiv, num_ddsv, num_ddrv, num_cauv, num_dauiv,
100 num_dausv, num_daurv, num_ceuv, num_deuiv, num_deusv, num_deurv,
101 num_csv, num_dsiv, num_dssv, num_dsrv, rv_cntr = 0, cv_cntr = 0,
102 div_cntr = 0, dsv_cntr = 0, drv_cntr = 0;
103 x_svd.design_counts(num_cdv, num_ddiv, num_ddsv, num_ddrv);
104 x_svd.aleatory_uncertain_counts(num_cauv, num_dauiv, num_dausv, num_daurv);
105 x_svd.epistemic_uncertain_counts(num_ceuv, num_deuiv, num_deusv, num_deurv);
106 x_svd.state_counts(num_csv, num_dsiv, num_dssv, num_dsrv);
107 const Pecos::ShortArray& u_types = mvDist.random_variable_types();
108
109 // Update active continuous/discrete variable types (needed for Model::
110 // continuous_{probability_density,distribution_bounds,distribution_moment}())
111 if (cdv)
112 for (i=0; i<num_cdv; ++i, ++rv_cntr, ++cv_cntr)
113 continuous_variable_type(
114 pecos_to_dakota_variable_type(u_types[rv_cntr], rv_cntr), cv_cntr);
115 else
116 rv_cntr += num_cdv;
117 if (ddv) {
118 for (i=0; i<num_ddiv; ++i, ++rv_cntr, ++div_cntr)
119 discrete_int_variable_type(
120 pecos_to_dakota_variable_type(u_types[rv_cntr], rv_cntr), div_cntr);
121 for (i=0; i<num_ddsv; ++i, ++rv_cntr, ++dsv_cntr)
122 discrete_string_variable_type(
123 pecos_to_dakota_variable_type(u_types[rv_cntr], rv_cntr), dsv_cntr);
124 for (i=0; i<num_ddrv; ++i, ++rv_cntr, ++drv_cntr)
125 discrete_real_variable_type(
126 pecos_to_dakota_variable_type(u_types[rv_cntr], rv_cntr), drv_cntr);
127 }
128 else
129 rv_cntr += num_ddiv + num_ddsv + num_ddrv;
130 if (cauv)
131 for (i=0; i<num_cauv; ++i, ++rv_cntr, ++cv_cntr)
132 continuous_variable_type(
133 pecos_to_dakota_variable_type(u_types[rv_cntr], rv_cntr), cv_cntr);
134 else
135 rv_cntr += num_cauv;
136 if (dauv) {
137 for (i=0; i<num_dauiv; ++i, ++rv_cntr, ++div_cntr)
138 discrete_int_variable_type(
139 pecos_to_dakota_variable_type(u_types[rv_cntr], rv_cntr), div_cntr);
140 for (i=0; i<num_dausv; ++i, ++rv_cntr, ++dsv_cntr)
141 discrete_string_variable_type(
142 pecos_to_dakota_variable_type(u_types[rv_cntr], rv_cntr), dsv_cntr);
143 for (i=0; i<num_daurv; ++i, ++rv_cntr, ++drv_cntr)
144 discrete_real_variable_type(
145 pecos_to_dakota_variable_type(u_types[rv_cntr], rv_cntr), drv_cntr);
146 }
147 else
148 rv_cntr += num_dauiv + num_dausv + num_daurv;
149 if (ceuv)
150 for (i=0; i<num_ceuv; ++i, ++rv_cntr, ++cv_cntr)
151 continuous_variable_type(
152 pecos_to_dakota_variable_type(u_types[rv_cntr], rv_cntr), cv_cntr);
153 else
154 rv_cntr += num_ceuv;
155 if (deuv) {
156 for (i=0; i<num_deuiv; ++i, ++rv_cntr, ++div_cntr)
157 discrete_int_variable_type(
158 pecos_to_dakota_variable_type(u_types[rv_cntr], rv_cntr), div_cntr);
159 for (i=0; i<num_deusv; ++i, ++rv_cntr, ++dsv_cntr)
160 discrete_string_variable_type(
161 pecos_to_dakota_variable_type(u_types[rv_cntr], rv_cntr), dsv_cntr);
162 for (i=0; i<num_deurv; ++i, ++rv_cntr, ++drv_cntr)
163 discrete_real_variable_type(
164 pecos_to_dakota_variable_type(u_types[rv_cntr], rv_cntr), drv_cntr);
165 }
166 else
167 rv_cntr += num_deuiv + num_deusv + num_deurv;
168 if (csv)
169 for (i=0; i<num_csv; ++i, ++rv_cntr, ++cv_cntr)
170 continuous_variable_type(
171 pecos_to_dakota_variable_type(u_types[rv_cntr], rv_cntr), cv_cntr);
172 else
173 rv_cntr += num_csv;
174 if (dsv) {
175 for (i=0; i<num_dsiv; ++i, ++rv_cntr, ++div_cntr)
176 discrete_int_variable_type(
177 pecos_to_dakota_variable_type(u_types[rv_cntr], rv_cntr), div_cntr);
178 for (i=0; i<num_dssv; ++i, ++rv_cntr, ++dsv_cntr)
179 discrete_string_variable_type(
180 pecos_to_dakota_variable_type(u_types[rv_cntr], rv_cntr), dsv_cntr);
181 for (i=0; i<num_dsrv; ++i, ++rv_cntr, ++drv_cntr)
182 discrete_real_variable_type(
183 pecos_to_dakota_variable_type(u_types[rv_cntr], rv_cntr), drv_cntr);
184 }
185 //else
186 // rv_cntr += num_dsiv + num_dssv + num_dsrv;
187 }
188
189
190 void ProbabilityTransformModel::
update_model_bounds(bool truncate_bnds,Real bnd)191 update_model_bounds(bool truncate_bnds, Real bnd)
192 {
193 // Here, we finesse the continuous "global" bounds for a Model (distinct from
194 // distribution bounds in Pecos::MultivariateDistribution), for benefit of
195 // methods that are not distribution-aware, instead operating on bounded
196 // domains (e.g., PStudyDACE). While the same concept could be extended to
197 // discrete distributions with semi-infinite support, we leave that exercise
198 // to be motivated by future use cases.
199
200 ////////////////////////////////////////////////////////////
201 // Populate continuous model bounds for transformed space //
202 ////////////////////////////////////////////////////////////
203
204 const Pecos::ShortArray& u_types = mvDist.random_variable_types();
205 const std::vector<Pecos::RandomVariable>& x_rv
206 = subModel.multivariate_distribution().random_variables();
207 size_t num_cv = currentVariables.cv(), num_rv = u_types.size();
208 // [-1,1] are u-space bounds for design, state, epistemic, uniform, & beta
209 RealVector c_l_bnds(num_cv, false); c_l_bnds = -1.;
210 RealVector c_u_bnds(num_cv, false); c_u_bnds = 1.;
211 Real dbl_inf = std::numeric_limits<Real>::infinity();
212
213 const SharedVariablesData& x_svd = subModel.current_variables().shared_data();
214 bool cdv, ddv, cauv, dauv, ceuv, deuv, csv, dsv;
215 x_svd.active_subsets(cdv, ddv, cauv, dauv, ceuv, deuv, csv, dsv);
216 size_t i, num_cdv, num_ddiv, num_ddsv, num_ddrv, num_cauv, num_dauiv,
217 num_dausv, num_daurv, num_ceuv, num_deuiv, num_deusv, num_deurv,
218 num_csv, num_dsiv, num_dssv, num_dsrv, cv_cntr = 0, rv_cntr = 0;
219 x_svd.design_counts(num_cdv, num_ddiv, num_ddsv, num_ddrv);
220 x_svd.aleatory_uncertain_counts(num_cauv, num_dauiv, num_dausv, num_daurv);
221 x_svd.epistemic_uncertain_counts(num_ceuv, num_deuiv, num_deusv, num_deurv);
222 x_svd.state_counts(num_csv, num_dsiv, num_dssv, num_dsrv);
223
224 if (truncate_bnds) {
225 // truncate unbounded distributions for approaches requiring bounds:
226 // > standard sampling modes: model bounds only used for design/state
227 // > *_UNIFORM modes: model bounds are used for all active variables
228
229 // all cdv are mapped to [-1,1]
230 if (cdv) cv_cntr += num_cdv;
231 rv_cntr += num_cdv + num_ddiv + num_ddsv + num_ddrv;
232
233 if (cauv) {
234 for (i=0; i<num_cauv; ++i, ++cv_cntr, ++rv_cntr) {
235 const Pecos::RandomVariable& rv_i = x_rv[rv_cntr];
236 switch (u_types[rv_cntr]) {
237 case Pecos::STD_NORMAL: // mean +/- bnd std devs
238 c_l_bnds[cv_cntr] = -bnd; c_u_bnds[cv_cntr] = bnd; break;
239 case Pecos::STD_EXPONENTIAL: // [0, mean + bnd std devs] for beta=1
240 c_l_bnds[cv_cntr] = 0.; c_u_bnds[cv_cntr] = 1.+bnd; break;
241 case Pecos::STD_GAMMA: {
242 Real alpha, mean, stdev;
243 rv_i.pull_parameter(Pecos::GA_ALPHA, alpha);
244 Pecos::GammaRandomVariable::
245 moments_from_params(alpha, 1., mean, stdev);
246 c_l_bnds[cv_cntr] = 0.;
247 c_u_bnds[cv_cntr] = mean + bnd * stdev; break;
248 }
249 case Pecos::BOUNDED_NORMAL: {
250 // Note: as for NIDR initialization, we use mean,std_dev parameters
251 // rather than computing actual mean,std_dev of bounded distribution
252 Real mean, stdev, l_bnd, u_bnd;
253 rv_i.pull_parameter(Pecos::N_MEAN, mean);
254 rv_i.pull_parameter(Pecos::N_STD_DEV, stdev);
255 rv_i.pull_parameter(Pecos::N_LWR_BND, l_bnd);
256 rv_i.pull_parameter(Pecos::N_UPR_BND, u_bnd);
257 c_l_bnds[cv_cntr] = (l_bnd > -dbl_inf) ? l_bnd : // use spec bound
258 mean - bnd * stdev; // infer bound
259 c_u_bnds[cv_cntr] = (u_bnd < dbl_inf) ? u_bnd : // use spec bound
260 mean + bnd * stdev; // infer bound
261 break;
262 }
263 case Pecos::LOGNORMAL: { // semi-bounded distribution
264 Real mean, stdev;
265 rv_i.pull_parameter(Pecos::LN_MEAN, mean);
266 rv_i.pull_parameter(Pecos::LN_STD_DEV, stdev);
267 c_l_bnds[cv_cntr] = 0.;
268 c_u_bnds[cv_cntr] = mean + bnd * stdev;
269 break;
270 }
271 case Pecos::BOUNDED_LOGNORMAL: {
272 Real l_bnd, u_bnd;
273 rv_i.pull_parameter(Pecos::LN_LWR_BND, l_bnd);
274 rv_i.pull_parameter(Pecos::LN_UPR_BND, u_bnd);
275 c_l_bnds[cv_cntr] = l_bnd; // spec or 0
276 if (u_bnd < dbl_inf)
277 c_u_bnds[cv_cntr] = u_bnd; // use spec bound
278 else { // infer bound
279 Real mean, stdev;
280 rv_i.pull_parameter(Pecos::LN_MEAN, mean);
281 rv_i.pull_parameter(Pecos::LN_STD_DEV, stdev);
282 // Note: as for NIDR initialization, we use mean,std_dev parameters
283 // rather than computing actual mean,std_dev of bounded distribution
284 c_u_bnds[cv_cntr] = mean + bnd * stdev;
285 }
286 break;
287 }
288 case Pecos::LOGUNIFORM: case Pecos::TRIANGULAR:
289 case Pecos::HISTOGRAM_BIN:
290 // bounded distributions: x-space has desired bounds
291 c_l_bnds[cv_cntr] = subModel.continuous_lower_bound(cv_cntr);
292 c_u_bnds[cv_cntr] = subModel.continuous_upper_bound(cv_cntr);
293 break;
294 // Note: Could use subModel bounds for the following cases as well
295 // except NIDR uses +/-3 sigma, whereas here we're using +/-10 sigma
296 case Pecos::GUMBEL: { // unbounded distribution
297 Real alpha, beta, mean, stdev;
298 rv_i.pull_parameter(Pecos::GU_ALPHA, alpha);
299 rv_i.pull_parameter(Pecos::GU_BETA, beta);
300 Pecos::GumbelRandomVariable::
301 moments_from_params(alpha, beta, mean, stdev);
302 c_l_bnds[cv_cntr] = mean - bnd * stdev;
303 c_u_bnds[cv_cntr] = mean + bnd * stdev;
304 break;
305 }
306 case Pecos::FRECHET: { // semibounded distribution
307 Real alpha, beta, mean, stdev;
308 rv_i.pull_parameter(Pecos::F_ALPHA, alpha);
309 rv_i.pull_parameter(Pecos::F_BETA, beta);
310 Pecos::FrechetRandomVariable::
311 moments_from_params(alpha, beta, mean, stdev);
312 c_l_bnds[cv_cntr] = 0.; c_u_bnds[cv_cntr] = mean + bnd * stdev; break;
313 }
314 case Pecos::WEIBULL: { // semibounded distribution
315 Real alpha, beta, mean, stdev;
316 rv_i.pull_parameter(Pecos::W_ALPHA, alpha);
317 rv_i.pull_parameter(Pecos::W_BETA, beta);
318 Pecos::WeibullRandomVariable::
319 moments_from_params(alpha, beta, mean, stdev);
320 c_l_bnds[cv_cntr] = 0.; c_u_bnds[cv_cntr] = mean + bnd * stdev; break;
321 }
322 }
323 }
324 }
325 else
326 rv_cntr += num_cauv;
327 rv_cntr += num_dauiv + num_dausv + num_daurv;
328
329 if (ceuv) {
330 for (i=0; i<num_ceuv; ++i, ++cv_cntr, ++rv_cntr) {
331 const Pecos::RandomVariable& rv_i = x_rv[rv_cntr];
332 switch (u_types[rv_cntr]) {
333 case Pecos::CONTINUOUS_INTERVAL_UNCERTAIN:
334 // bounded distributions: x-space has desired bounds
335 c_l_bnds[cv_cntr] = subModel.continuous_lower_bound(cv_cntr);
336 c_u_bnds[cv_cntr] = subModel.continuous_upper_bound(cv_cntr);
337 break;
338 }
339 }
340 }
341 //else
342 // rv_cntr += num_ceuv;
343 //rv_cntr += num_deuiv + num_deusv + num_deurv;
344
345 // all csv are mapped to [-1,1]
346 //if (csv) cv_cntr += num_csv;
347 //rv_cntr += num_csv + num_dsiv + num_dssv + num_dsrv;
348 }
349 else { // retain infinite model bounds where distributions are unbounded
350
351 // all cdv are mapped to [-1,1]
352 if (cdv) cv_cntr += num_cdv;
353 rv_cntr += num_cdv + num_ddiv + num_ddsv + num_ddrv;
354
355 if (cauv) {
356 for (i=0; i<num_cauv; ++i, ++cv_cntr, ++rv_cntr) {
357 const Pecos::RandomVariable& rv_i = x_rv[rv_cntr];
358 switch (u_types[rv_cntr]) {
359 case Pecos::STD_NORMAL: case Pecos::GUMBEL: // unbounded distributions
360 c_l_bnds[cv_cntr] = -dbl_inf; c_u_bnds[cv_cntr] = dbl_inf; break;
361 case Pecos::LOGNORMAL: case Pecos::STD_EXPONENTIAL:
362 case Pecos::STD_GAMMA: case Pecos::FRECHET:
363 case Pecos::WEIBULL: // semibounded distributions
364 c_l_bnds[cv_cntr] = 0.; c_u_bnds[cv_cntr] = dbl_inf; break;
365 case Pecos::BOUNDED_NORMAL:
366 // can't rely on subModel bounds since could be 1-sided
367 rv_i.pull_parameter(Pecos::N_LWR_BND, c_l_bnds[cv_cntr]);
368 rv_i.pull_parameter(Pecos::N_UPR_BND, c_u_bnds[cv_cntr]);
369 break;
370 case Pecos::BOUNDED_LOGNORMAL:
371 // can't rely on subModel bounds since could be 1-sided
372 rv_i.pull_parameter(Pecos::LN_LWR_BND, c_l_bnds[cv_cntr]);
373 rv_i.pull_parameter(Pecos::LN_UPR_BND, c_u_bnds[cv_cntr]);
374 break;
375 case Pecos::LOGUNIFORM: case Pecos::TRIANGULAR:
376 case Pecos::HISTOGRAM_BIN: // bounded distributions
377 // 2-sided: can rely on subModel bounds
378 c_l_bnds[cv_cntr] = subModel.continuous_lower_bound(cv_cntr);
379 c_u_bnds[cv_cntr] = subModel.continuous_upper_bound(cv_cntr); break;
380 }
381 }
382 }
383 else
384 rv_cntr += num_cauv;
385 rv_cntr += num_dauiv + num_dausv + num_daurv;
386
387 if (ceuv) {
388 for (i=0; i<num_ceuv; ++i, ++cv_cntr, ++rv_cntr) {
389 const Pecos::RandomVariable& rv_i = x_rv[rv_cntr];
390 switch (u_types[rv_cntr]) {
391 case Pecos::CONTINUOUS_INTERVAL_UNCERTAIN:
392 // bounded distributions: x-space has desired bounds
393 c_l_bnds[cv_cntr] = subModel.continuous_lower_bound(cv_cntr);
394 c_u_bnds[cv_cntr] = subModel.continuous_upper_bound(cv_cntr);
395 break;
396 }
397 }
398 }
399 //else
400 // rv_cntr += num_ceuv;
401 //rv_cntr += num_deuiv + num_deusv + num_deurv;
402
403 // all csv are mapped to [-1,1]
404 //if (csv) cv_cntr += num_csv;
405 //rv_cntr += num_csv + num_dsiv + num_dssv + num_dsrv;
406 }
407
408 continuous_lower_bounds(c_l_bnds); continuous_upper_bounds(c_u_bnds);
409 }
410
411
412 /** Build ProbabilityTransformation::ranVar arrays containing the
413 uncertain variable distribution types and their corresponding
414 means/standard deviations. This function is used when the Model
415 variables are in x-space. */
416 void ProbabilityTransformModel::
initialize_distribution_types(short u_space_type,const Pecos::MultivariateDistribution & x_dist,Pecos::MultivariateDistribution & u_dist)417 initialize_distribution_types(short u_space_type,
418 const Pecos::MultivariateDistribution& x_dist,
419 Pecos::MultivariateDistribution& u_dist)
420 {
421 // u_space_type is an enumeration for type of u-space transformation:
422 // > if STD_NORMAL_U (reliability, AIS, and Wiener PCE/SC), then u-space is
423 // defined exclusively with std normals;
424 // > if STD_UNIFORM_U (SC with local & global Hermite basis polynomials),
425 // then u-space is defined exclusively with std uniforms;
426 // > if PARTIAL_ASKEY_U (C3 with orthog polynomials), then u-space is defined
427 // by std normals and std uniforms;
428 // > if ASKEY_U (PCE/SC using Askey polynomials), then u-space is defined by
429 // std normals, std uniforms, std exponentials, std betas, and std gammas;
430 // > if EXTENDED_U (PCE/SC with Askey plus numerically-generated polynomials),
431 // then u-space involves at most linear scaling to std distributions.
432
433 const Pecos::ShortArray& x_types = x_dist.random_variable_types();
434 const Pecos::BitArray& active_rv = x_dist.active_variables();
435 size_t i, num_rv = x_types.size();
436 Pecos::ShortArray u_types(num_rv, Pecos::NO_TYPE);
437 bool err_flag = false;
438
439 for (i=0; i<num_rv; ++i)
440 if (active_rv[i])
441 switch (u_space_type) {
442 case STD_NORMAL_U: case STD_UNIFORM_U:
443 switch (x_types[i]) {
444 case Pecos::DISCRETE_RANGE: case Pecos::DISCRETE_SET_INT:
445 case Pecos::DISCRETE_SET_STRING: case Pecos::DISCRETE_SET_REAL:
446 case Pecos::POISSON: case Pecos::BINOMIAL:
447 case Pecos::NEGATIVE_BINOMIAL: case Pecos::GEOMETRIC:
448 case Pecos::HYPERGEOMETRIC: case Pecos::HISTOGRAM_PT_INT:
449 case Pecos::HISTOGRAM_PT_STRING: case Pecos::HISTOGRAM_PT_REAL:
450 case Pecos::DISCRETE_INTERVAL_UNCERTAIN:
451 case Pecos::DISCRETE_UNCERTAIN_SET_INT:
452 case Pecos::DISCRETE_UNCERTAIN_SET_STRING:
453 case Pecos::DISCRETE_UNCERTAIN_SET_REAL:
454 err_flag = true; break;
455 case Pecos::CONTINUOUS_RANGE: case Pecos::CONTINUOUS_INTERVAL_UNCERTAIN:
456 u_types[i] = Pecos::STD_UNIFORM; break;
457 default:
458 u_types[i] = (u_space_type == STD_UNIFORM_U)
459 ? Pecos::STD_UNIFORM : Pecos::STD_NORMAL; break;
460 }
461 break;
462 case PARTIAL_ASKEY_U: // used for cases with limited distrib support (C3)
463 switch (x_types[i]) {
464 case Pecos::NORMAL: case Pecos::BOUNDED_NORMAL:
465 case Pecos::LOGNORMAL: case Pecos::BOUNDED_LOGNORMAL:
466 case Pecos::EXPONENTIAL: case Pecos::GAMMA:
467 case Pecos::GUMBEL: case Pecos::FRECHET:
468 case Pecos::WEIBULL: // unbounded or semi-bounded dist; bounded N/logN
469 u_types[i] = Pecos::STD_NORMAL; break;
470 case Pecos::UNIFORM: case Pecos::LOGUNIFORM:
471 case Pecos::TRIANGULAR: case Pecos::BETA:
472 case Pecos::HISTOGRAM_BIN: case Pecos::CONTINUOUS_RANGE:
473 case Pecos::CONTINUOUS_INTERVAL_UNCERTAIN: // bounded
474 u_types[i] = Pecos::STD_UNIFORM; break;
475 // TO DO: discrete types ddv, dauv, deuv, dsv
476 default: err_flag = true; break;
477 }
478 break;
479 case ASKEY_U:
480 switch (x_types[i]) {
481 case Pecos::NORMAL: case Pecos::BOUNDED_NORMAL:
482 case Pecos::LOGNORMAL: case Pecos::BOUNDED_LOGNORMAL:
483 case Pecos::GUMBEL: case Pecos::FRECHET:
484 case Pecos::WEIBULL:
485 u_types[i] = Pecos::STD_NORMAL; break;
486 case Pecos::UNIFORM: case Pecos::LOGUNIFORM:
487 case Pecos::TRIANGULAR: case Pecos::HISTOGRAM_BIN:
488 case Pecos::CONTINUOUS_RANGE: case Pecos::CONTINUOUS_INTERVAL_UNCERTAIN:
489 u_types[i] = Pecos::STD_UNIFORM; break;
490 case Pecos::EXPONENTIAL: u_types[i] = Pecos::STD_EXPONENTIAL; break;
491 case Pecos::BETA: u_types[i] = Pecos::STD_BETA; break;
492 case Pecos::GAMMA: u_types[i] = Pecos::STD_GAMMA; break;
493 // TO DO: discrete types ddv, dauv, deuv, dsv
494 //case Pecos::POISSON: case Pecos::BINOMIAL:
495 //case Pecos::NEGATIVE_BINOMIAL: case Pecos::GEOMETRIC:
496 //case Pecos::HYPERGEOMETRIC:
497 default: err_flag = true; break;
498 }
499 break;
500 case EXTENDED_U:
501 switch (x_types[i]) {
502 case Pecos::CONTINUOUS_RANGE: case Pecos::UNIFORM:
503 u_types[i] = Pecos::STD_UNIFORM; break;
504 case Pecos::NORMAL: u_types[i] = Pecos::STD_NORMAL; break;
505 case Pecos::EXPONENTIAL: u_types[i] = Pecos::STD_EXPONENTIAL; break;
506 case Pecos::BETA: u_types[i] = Pecos::STD_BETA; break;
507 case Pecos::GAMMA: u_types[i] = Pecos::STD_GAMMA; break;
508 default: u_types[i] = x_types[i]; break;
509 }
510 break;
511 }
512 else // inactive variables are not transformed
513 u_types[i] = x_types[i];
514
515 if (err_flag) {
516 Cerr << "Error: unsupported mapping in ProbabilityTransformModel::"
517 << "initialize_distribution_transformation()." << std::endl;
518 abort_handler(MODEL_ERROR);
519 }
520
521 std::shared_ptr<Pecos::MarginalsCorrDistribution> u_dist_rep =
522 std::static_pointer_cast<Pecos::MarginalsCorrDistribution>
523 (u_dist.multivar_dist_rep());
524 u_dist_rep->initialize_types(u_types, active_rv);
525 }
526
527
verify_correlation_support(short u_space_type)528 void ProbabilityTransformModel::verify_correlation_support(short u_space_type)
529 {
530 Pecos::MultivariateDistribution& x_dist
531 = subModel.multivariate_distribution();
532 if (x_dist.correlation()) {
533 const Pecos::ShortArray& x_types = x_dist.random_variable_types();
534 const Pecos::ShortArray& u_types = mvDist.random_variable_types();
535 const Pecos::RealSymMatrix& x_corr = x_dist.correlation_matrix();
536 const Pecos::BitArray& active_corr = x_dist.active_correlations();
537 size_t i, j, corr_i, corr_j, num_rv = x_types.size();
538 bool no_mask = active_corr.empty();
539
540 // We can only decorrelate in std normal space; therefore, if a variable
541 // with a u_type other than STD_NORMAL is correlated with anything, change
542 // its u_type to STD_NORMAL.
543 if (u_space_type != STD_NORMAL_U) {
544 for (i=0, corr_i=0; i<num_rv; ++i)
545 if (no_mask || active_corr[i]) {
546 if (u_types[i] != Pecos::STD_NORMAL)
547 // since we don't check all rows, check all columns despite symmetry
548 for (j=0, corr_j=0; j<num_rv; ++j)
549 if (no_mask || active_corr[j]) {
550 if (i != j && std::abs(x_corr(corr_i, corr_j)) >
551 Pecos::SMALL_NUMBER) {
552 Cerr << "\nWarning: u-space type for random variable " << i+1
553 << " changed to\n STD_NORMAL due to "
554 << "decorrelation requirements.\n";
555 mvDist.random_variable_type(Pecos::STD_NORMAL, i);
556 break; // out of inner loop
557 }
558 ++corr_j;
559 }
560 ++corr_i;
561 }
562
563 /*
564 for (i=0; i<num_rv; ++i)
565 if (u_types[i] != Pecos::STD_NORMAL) {
566 // since we don't check all rows, check all columns despite symmetry
567 corr_i = rv_index_to_corr_index(i);
568 if (corr_i != _NPOS)
569 for (j=0, corr_j=0; j<num_rv; ++j)
570 if (i != j) {
571 corr_j = rv_index_to_corr_index(j);
572 if (corr_j != _NPOS &&
573 std::abs(x_corr(corr_i, corr_j)) > Pecos::SMALL_NUMBER) {
574 Cerr << "\nWarning: u-space type for random variable " << i+1
575 << " changed to\n STD_NORMAL due to "
576 << "decorrelation requirements.\n";
577 mvDist.random_variable_type(Pecos::STD_NORMAL, i);
578 break; // out of inner loop
579 }
580 }
581 }
582 */
583 }
584
585 // Check for correlations among variable types (bounded normal, bounded
586 // lognormal, loguniform, triangular, beta, and histogram) that are not
587 // supported by Der Kiureghian & Liu for correlation warping estimation
588 // when transforming to std normals.
589 bool err_flag = false;
590 for (i=0, corr_i=0; i<num_rv; ++i) {
591 bool distribution_error = false;
592 if (no_mask || active_corr[i]) {
593 short x_type = x_types[i];
594 if (x_type == Pecos::BOUNDED_NORMAL || x_type == Pecos::LOGUNIFORM ||
595 x_type == Pecos::BOUNDED_LOGNORMAL || x_type == Pecos::TRIANGULAR ||
596 x_type == Pecos::BETA || x_type == Pecos::HISTOGRAM_BIN)
597 // since we don't check all rows, check all columns despite symmetry
598 for (j=0, corr_j=0; j<num_rv; ++j)
599 if (no_mask || active_corr[j]) {
600 if (i != j &&
601 std::abs(x_corr(corr_i, corr_j)) > Pecos::SMALL_NUMBER)
602 { distribution_error = true; break; }
603 ++corr_j;
604 }
605 ++corr_i;
606 }
607 if (distribution_error) {
608 Cerr << "Error: correlation warping for Nataf variable transformation "
609 << "of bounded normal,\n bounded lognormal, loguniform, "
610 << "triangular, beta, and histogram bin\n distributions is "
611 << "not currently supported. Error detected for variable " << i+1
612 << "." << std::endl;
613 err_flag = true;
614 }
615 }
616 if (err_flag)
617 abort_handler(MODEL_ERROR);
618 }
619 }
620
621
622 unsigned short ProbabilityTransformModel::
pecos_to_dakota_variable_type(unsigned short pecos_var_type,size_t rv_index)623 pecos_to_dakota_variable_type(unsigned short pecos_var_type, size_t rv_index)
624 {
625 const SizetArray& vc_totals
626 = subModel.current_variables().shared_data().components_totals();
627 switch (pecos_var_type) {
628 case Pecos::CONTINUOUS_RANGE: // non-unique mapping
629 return (rv_index < vc_totals[TOTAL_CDV]) ? // not subject to active subsets
630 CONTINUOUS_DESIGN : CONTINUOUS_STATE; break;
631 case Pecos::DISCRETE_RANGE: // non-unique mapping
632 return (rv_index < vc_totals[TOTAL_CDV] + vc_totals[TOTAL_DDIV]) ?
633 DISCRETE_DESIGN_RANGE : DISCRETE_STATE_RANGE; break;
634 case Pecos::DISCRETE_SET_INT: // non-unique mapping
635 return (rv_index < vc_totals[TOTAL_CDV] + vc_totals[TOTAL_DDIV]) ?
636 DISCRETE_DESIGN_SET_INT : DISCRETE_STATE_SET_INT; break;
637 case Pecos::DISCRETE_SET_STRING: // non-unique mapping
638 return (rv_index < vc_totals[TOTAL_CDV] + vc_totals[TOTAL_DDIV] +
639 vc_totals[TOTAL_DDSV]) ?
640 DISCRETE_DESIGN_SET_STRING : DISCRETE_STATE_SET_STRING; break;
641 case Pecos::DISCRETE_SET_REAL: // non-unique mapping
642 return (rv_index < vc_totals[TOTAL_CDV] + vc_totals[TOTAL_DDIV] +
643 vc_totals[TOTAL_DDSV] + vc_totals[TOTAL_DDRV]) ?
644 DISCRETE_DESIGN_SET_REAL : DISCRETE_STATE_SET_REAL; break;
645 // continuous aleatory
646 case Pecos::STD_NORMAL: case Pecos::NORMAL: case Pecos::BOUNDED_NORMAL:
647 return NORMAL_UNCERTAIN; break;
648 case Pecos::LOGNORMAL: case Pecos::BOUNDED_LOGNORMAL:
649 return LOGNORMAL_UNCERTAIN; break;
650 case Pecos::STD_UNIFORM: case Pecos::UNIFORM:
651 return UNIFORM_UNCERTAIN; break;
652 case Pecos::LOGUNIFORM:
653 return LOGUNIFORM_UNCERTAIN; break;
654 case Pecos::TRIANGULAR:
655 return TRIANGULAR_UNCERTAIN; break;
656 case Pecos::STD_EXPONENTIAL: case Pecos::EXPONENTIAL:
657 return EXPONENTIAL_UNCERTAIN; break;
658 case Pecos::STD_BETA: case Pecos::BETA:
659 return BETA_UNCERTAIN; break;
660 case Pecos::STD_GAMMA: case Pecos::GAMMA:
661 return GAMMA_UNCERTAIN; break;
662 case Pecos::GUMBEL:
663 return GUMBEL_UNCERTAIN; break;
664 case Pecos::FRECHET:
665 return FRECHET_UNCERTAIN; break;
666 case Pecos::WEIBULL:
667 return WEIBULL_UNCERTAIN; break;
668 case Pecos::HISTOGRAM_BIN:
669 return HISTOGRAM_BIN_UNCERTAIN; break;
670 // discrete aleatory
671 case Pecos::POISSON:
672 return POISSON_UNCERTAIN; break;
673 case Pecos::BINOMIAL:
674 return BINOMIAL_UNCERTAIN; break;
675 case Pecos::NEGATIVE_BINOMIAL:
676 return NEGATIVE_BINOMIAL_UNCERTAIN; break;
677 case Pecos::GEOMETRIC:
678 return GEOMETRIC_UNCERTAIN; break;
679 case Pecos::HYPERGEOMETRIC:
680 return HYPERGEOMETRIC_UNCERTAIN; break;
681 case Pecos::HISTOGRAM_PT_INT:
682 return HISTOGRAM_POINT_UNCERTAIN_INT; break;
683 case Pecos::HISTOGRAM_PT_STRING:
684 return HISTOGRAM_POINT_UNCERTAIN_STRING; break;
685 case Pecos::HISTOGRAM_PT_REAL:
686 return HISTOGRAM_POINT_UNCERTAIN_REAL; break;
687 // continuous epistemic
688 case Pecos::CONTINUOUS_INTERVAL_UNCERTAIN:
689 return CONTINUOUS_INTERVAL_UNCERTAIN; break;
690 // discrete epistemic
691 case Pecos::DISCRETE_INTERVAL_UNCERTAIN:
692 return DISCRETE_INTERVAL_UNCERTAIN; break;
693 case Pecos:: DISCRETE_UNCERTAIN_SET_INT:
694 return DISCRETE_UNCERTAIN_SET_INT; break;
695 case Pecos::DISCRETE_UNCERTAIN_SET_STRING:
696 return DISCRETE_UNCERTAIN_SET_STRING; break;
697 case Pecos::DISCRETE_UNCERTAIN_SET_REAL:
698 return DISCRETE_UNCERTAIN_SET_REAL; break;
699 default:
700 Cerr << "Error: unsupported Pecos distribution type in "
701 << "pecos_to_dakota_variable_type()." << std::endl;
702 abort_handler(MODEL_ERROR); return 0; break;
703 }
704 }
705
706
707 void ProbabilityTransformModel::
resp_x_to_u_mapping(const Variables & x_vars,const Variables & u_vars,const Response & x_response,Response & u_response)708 resp_x_to_u_mapping(const Variables& x_vars, const Variables& u_vars,
709 const Response& x_response, Response& u_response)
710 {
711 const RealVector& x_cv = x_vars.continuous_variables();
712 SizetMultiArrayConstView x_cv_ids = x_vars.continuous_variable_ids();
713 SizetMultiArrayConstView x_acv_ids = x_vars.all_continuous_variable_ids();
714 const RealVector& x_fns = x_response.function_values();
715
716 // In this recasting, the inputs and outputs are mapped one-to-one, with no
717 // reordering. However, the x-space ASV may be augmented from the original
718 // u-space ASV due to nonlinear mapping logic in RecastModel::transform_set().
719 const ShortArray& u_asv = u_response.active_set_request_vector();
720 const SizetArray& u_dvv = u_response.active_set_derivative_vector();
721 const ShortArray& x_asv = x_response.active_set_request_vector();
722 const SizetArray& x_dvv = x_response.active_set_derivative_vector();
723 Pecos::MultivariateDistribution& x_dist
724 = ptmInstance->subModel.multivariate_distribution();
725 size_t i, j, num_fns = x_asv.size(), num_deriv_vars = x_dvv.size();
726 if (u_asv.size() != num_fns) {
727 Cerr << "Error: inconsistent response function definition in Probability"
728 << "TransformModel::resp_x_to_u_mapping().\n x-space response "
729 << "size = " << num_fns << ", u-space response size =\n"
730 << u_asv.size() << std::endl;
731 abort_handler(MODEL_ERROR);
732 }
733 if (!x_dist.correlation() && u_dvv != x_dvv) {
734 Cerr << "Error: inconsistent derivative component definition in Probability"
735 << "TransformModel::resp_x_to_u_mapping().\nx-space DVV =\n" << x_dvv
736 << "u-space DVV =\n" << u_dvv << std::endl;
737 abort_handler(MODEL_ERROR);
738 }
739 bool u_grad_flag = false, u_hess_flag = false;
740 for (i=0; i<num_fns; ++i) {
741 if (u_asv[i] & 2)
742 u_grad_flag = true;
743 if (u_asv[i] & 4)
744 u_hess_flag = true;
745 }
746
747 bool map_derivs = ( (u_grad_flag || u_hess_flag) &&
748 u_dvv != u_vars.inactive_continuous_variable_ids() );
749 RealVector fn_grad_x, fn_grad_us; RealSymMatrix fn_hess_us;
750 RealMatrix jacobian_xu, jacobian_xs; RealSymMatrixArray hessian_xu;
751
752 if (map_derivs) {
753 // The following transformation data is invariant w.r.t. the response fns
754 // and is computed outside of the num_fns loop
755 if (ptmInstance->distParamDerivs > NO_DERIVS)
756 ptmInstance->natafTransform.jacobian_dX_dS(x_cv, jacobian_xs,
757 x_cv_ids, x_acv_ids, ptmInstance->primaryACVarMapIndices,
758 ptmInstance->secondaryACVarMapTargets);
759 else {
760 if (u_grad_flag || u_hess_flag)
761 ptmInstance->natafTransform.jacobian_dX_dU(x_cv, jacobian_xu);
762 if (u_hess_flag && ptmInstance->nonlinearVarsMapping)
763 ptmInstance->natafTransform.hessian_d2X_dU2(x_cv, hessian_xu);
764 }
765 }
766
767 for (i=0; i<num_fns; ++i) {
768 short u_asv_val = u_asv[i]; // original request from iterator
769 short x_asv_val = x_asv[i]; // mapped request for sub-model
770
771 // map value g(x) to G(u)
772 if (u_asv_val & 1) {
773 if ( !(x_asv_val & 1) ) {
774 Cerr << "Error: missing required sub-model data in ProbabilityTransform"
775 << "Model::resp_x_to_u_mapping()" << std::endl;
776 abort_handler(MODEL_ERROR);
777 }
778 // no transformation: g(x) = G(u) by definition
779 u_response.function_value(x_fns[i], i);
780 }
781
782 // manage data requirements for derivative transformations: if fn_grad_x
783 // is needed for Hessian x-form (nonlinear I/O mapping), then x_asv has been
784 // augmented to include the gradient in RecastModel::transform_set().
785 if (x_asv_val & 2)
786 fn_grad_x = x_response.function_gradient_view(i);
787
788 // map gradient dg/dx to dG/du
789 if (u_asv_val & 2) {
790 if ( !(x_asv_val & 2) ) {
791 Cerr << "Error: missing required gradient sub-model data in Probability"
792 << "TransformModel::resp_x_to_u_mapping()" << std::endl;
793 abort_handler(MODEL_ERROR);
794 }
795 if (map_derivs) { // perform transformation
796 fn_grad_us = u_response.function_gradient_view(i);
797 if (ptmInstance->distParamDerivs > NO_DERIVS) // transform subset
798 ptmInstance->natafTransform.trans_grad_X_to_S(fn_grad_x,
799 fn_grad_us, jacobian_xs, x_dvv, x_cv_ids, x_acv_ids,
800 ptmInstance->primaryACVarMapIndices,
801 ptmInstance->secondaryACVarMapTargets);
802 else // transform subset of components
803 ptmInstance->natafTransform.trans_grad_X_to_U(fn_grad_x,
804 fn_grad_us, jacobian_xu, x_dvv, x_cv_ids);
805 }
806 else // no transformation: dg/dx = dG/du
807 u_response.function_gradient(fn_grad_x, i);
808 }
809
810 // map Hessian d^2g/dx^2 to d^2G/du^2
811 if (u_asv_val & 4) {
812 if ( !(x_asv_val & 4) || ( map_derivs &&
813 ptmInstance->nonlinearVarsMapping && !(x_asv_val & 2) ) ) {
814 Cerr << "Error: missing required sub-model data in Probability"
815 << "TransformModel::resp_x_to_u_mapping()" << std::endl;
816 abort_handler(MODEL_ERROR);
817 }
818 const RealSymMatrix& fn_hess_x = x_response.function_hessian(i);
819 if (map_derivs) { // perform transformation
820 fn_hess_us = u_response.function_hessian_view(i);
821 if (ptmInstance->distParamDerivs > NO_DERIVS) { // transform subset
822 Cerr << "Error: Hessians with respect to inserted variables not yet "
823 << "supported." << std::endl;
824 abort_handler(MODEL_ERROR);
825 //ptmInstance->natafTransform.trans_hess_X_to_S(fn_hess_x,
826 // fn_hess_us, jacobian_xs, hessian_xs, fn_grad_s, x_dvv,
827 // x_cv_ids, x_vars.all_continuous_variable_ids(),
828 // ptmInstance->primaryACVarMapIndices,
829 // ptmInstance->secondaryACVarMapTargets);
830 }
831 else // transform subset of components
832 ptmInstance->natafTransform.trans_hess_X_to_U(fn_hess_x, fn_hess_us,
833 jacobian_xu, hessian_xu, fn_grad_x, x_dvv, x_cv_ids);
834 }
835 else // no transformation: d^2g/dx^2 = d^2G/du^2
836 u_response.function_hessian(fn_hess_x, i);
837 }
838 }
839
840 #ifdef DEBUG
841 Cout << "\nx_response:\n" << x_response
842 << "\nu_response:\n" << u_response << std::endl;
843 #endif
844 }
845
846
847 /** Define the DVV for x-space derivative evaluations by augmenting
848 the iterator requests to account for correlations. */
849 void ProbabilityTransformModel::
set_u_to_x_mapping(const Variables & u_vars,const ActiveSet & u_set,ActiveSet & x_set)850 set_u_to_x_mapping(const Variables& u_vars, const ActiveSet& u_set,
851 ActiveSet& x_set)
852 {
853 Pecos::MultivariateDistribution& x_dist
854 = ptmInstance->subModel.multivariate_distribution();
855 //if (ptmInstance->distParamDerivs > NO_DERIVS) {
856 //}
857 //else
858 if (x_dist.correlation()) {
859 const SizetArray& u_dvv = u_set.derivative_vector();
860 bool std_dvv = (u_dvv == u_vars.continuous_variable_ids() ||
861 u_dvv == u_vars.inactive_continuous_variable_ids());
862 if (!std_dvv) { // partial random variable derivatives: check correlations
863 SizetArray x_dvv;
864 SizetMultiArrayConstView acv_ids = u_vars.all_continuous_variable_ids();
865 size_t i, j, corr_i, corr_j, acv_id_i, num_acv = acv_ids.size();
866 const RealSymMatrix& corr_x = x_dist.correlation_matrix();
867 for (i=0; i<num_acv; ++i) { // insert in ascending order
868 acv_id_i = acv_ids[i];
869 if (contains(u_dvv, acv_id_i))
870 x_dvv.push_back(acv_id_i);
871 else {
872 corr_i = ptmInstance->acv_index_to_corr_index(i);
873 if (corr_i != _NPOS)
874 for (j=0; j<num_acv; ++j)
875 if (j != i) {
876 corr_j = ptmInstance->acv_index_to_corr_index(j);
877 if (corr_j != _NPOS &&
878 std::abs(corr_x(corr_i, corr_j)) > Pecos::SMALL_NUMBER &&
879 contains(u_dvv, acv_ids[j]))
880 { x_dvv.push_back(acv_id_i); break; }
881 }
882 }
883 }
884 x_set.derivative_vector(x_dvv);
885 }
886 }
887 // else copying DVV in RecastModel::transform_set() is sufficient
888 }
889
890 } // namespace Dakota
891