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: Specialization of RecastModel to transform a sub-model to
11 //-              u-space.
12 //- Owner:       Brian Adams
13 //- Checked by:
14 //- Version: $Id$
15 
16 #ifndef PROBABILITY_TRANSFORM_MODEL_H
17 #define PROBABILITY_TRANSFORM_MODEL_H
18 
19 #include "RecastModel.hpp"
20 #include "ProbabilityTransformation.hpp"
21 
22 namespace Dakota
23 {
24 
25 
26 /// Probability transformation specialization of RecastModel
27 
28 /** Specialization of RecastModel to transform a sub-model to u-space. */
29 class ProbabilityTransformModel: public RecastModel
30 {
31 public:
32 
33   //
34   //- Heading: Constructor and destructor
35   //
36 
37   /// standard constructor
38   ProbabilityTransformModel(const Model& sub_model, short u_space_type,
39                             bool truncate_bnds = false, Real bnd = 10.);
40 
41   /// destructor
42   ~ProbabilityTransformModel();
43 
44   //
45   //- Heading: Member functions
46   //
47 
48   /// initialize transformed distribution types and instantiate mvDist
49   static void initialize_distribution_types(short u_space_type,
50     const Pecos::MultivariateDistribution& x_dist,
51     Pecos::MultivariateDistribution& u_dist);
52 
53 protected:
54 
55   //
56   //- Heading: Virtual function redefinitions
57   //
58 
59   Pecos::ProbabilityTransformation& probability_transformation();
60 
61   //bool initialize_mapping(ParLevLIter pl_iter);
62   //bool finalize_mapping();
63   bool resize_pending() const;
64   void update_from_subordinate_model(size_t depth =
65 				     std::numeric_limits<size_t>::max());
66 
67   /// set primaryACVarMapIndices and secondaryACVarMapTargets (only, for now)
68   void nested_variable_mappings(const SizetArray& c_index1,
69 				const SizetArray& di_index1,
70 				const SizetArray& ds_index1,
71 				const SizetArray& dr_index1,
72 				const ShortArray& c_target2,
73 				const ShortArray& di_target2,
74 				const ShortArray& ds_target2,
75 				const ShortArray& dr_target2);
76   /// return primaryACVarMapIndices
77   const SizetArray& nested_acv1_indices() const;
78   /// return secondaryACVarMapTargets
79   const ShortArray& nested_acv2_targets() const;
80 
81   /// calculate and return potential state of distribution parameter
82   /// derivatives, but do not cache value in distParamDerivs
83   short query_distribution_parameter_derivatives() const;
84   /// activate distParamDerivs to {NO,MIXED,ALL}_DERIVS
85   void activate_distribution_parameter_derivatives();
86   /// reset distParamDerivs to NO_DERIVS
87   void deactivate_distribution_parameter_derivatives();
88 
89   void assign_instance();
90 
91   void trans_grad_X_to_U(const RealVector& fn_grad_x, RealVector& fn_grad_u,
92 			 const RealVector& x_vars);
93   void trans_grad_U_to_X(const RealVector& fn_grad_u, RealVector& fn_grad_x,
94 			 const RealVector& x_vars);
95   void trans_grad_X_to_S(const RealVector& fn_grad_x, RealVector& fn_grad_s,
96 			 const RealVector& x_vars);
97   void trans_hess_X_to_U(const RealSymMatrix& fn_hess_x,
98 			 RealSymMatrix& fn_hess_u, const RealVector& x_vars,
99 			 const RealVector& fn_grad_x);
100 
101   //
102   //- Heading: Member functions
103   //
104 
105   /// initialize transformed distribution types and natafTransform
106   /// (construct time)
107   void initialize_transformation(short u_space_type);
108   /// update with latest distribution data (run time)
109   void update_transformation();
110 
111   /// instantiate and initialize natafTransform
112   void initialize_nataf();
113 
114   // x-space correlations assigned in Model and u-space is uncorrelated
115   //void update_distribution_correlations();
116   /// verify that correlation warping is supported by Nataf for given
117   /// variable types
118   void verify_correlation_support(short u_space_type);
119 
120   /// initialize the continuous/discrete variable types using u-space types
121   /// (converted from Pecos to Dakota)
122   void initialize_dakota_variable_types();
123   /// update model bounds using u-space (truncated) distribution bounds
124   void update_model_bounds(bool truncate_bnds, Real bnd);
125 
126   /// detect when the variables transformation is nonlinear
127   bool nonlinear_variables_mapping(
128     const Pecos::MultivariateDistribution& x_dist,
129     const Pecos::MultivariateDistribution& u_dist) const;
130 
131   /// convert vector<RandomVariable> index to active correlation index
132   size_t rv_index_to_corr_index(size_t rv_index);
133   /// convert allContinuousVars index to active correlation index
134   size_t acv_index_to_corr_index(size_t acv_index);
135 
136   /// convert from Pecos To Dakota variable enumeration type for continuous
137   /// aleatory uncertain variables used in variable transformations
138   unsigned short pecos_to_dakota_variable_type(unsigned short pecos_var_type,
139 					       size_t rv_index);
140 
141   /// static function for RecastModels used for forward mapping of u-space
142   /// variables from NonD Iterators to x-space variables for Model evaluations
143   static void vars_u_to_x_mapping(const Variables& u_vars, Variables& x_vars);
144   /// static function for RecastModels used for inverse mapping of x-space
145   /// variables from data import to u-space variables for NonD Iterators
146   static void vars_x_to_u_mapping(const Variables& x_vars, Variables& u_vars);
147 
148   /// static function for RecastModels used to map u-space ActiveSets
149   /// from NonD Iterators to x-space ActiveSets for Model evaluations
150   static void set_u_to_x_mapping(const Variables& u_vars,
151                                  const ActiveSet& u_set, ActiveSet& x_set);
152 
153   /// static function for RecastModels used to map x-space responses from
154   /// Model evaluations to u-space responses for return to NonD Iterator.
155   static void resp_x_to_u_mapping(const Variables& x_vars,
156                                   const Variables& u_vars,
157                                   const Response& x_response,
158                                   Response& u_response);
159 
160 private:
161 
162   /// Nonlinear variable transformation that encapsulates the required
163   /// data for performing transformations from X -> Z -> U and back.
164   Pecos::ProbabilityTransformation natafTransform;
165 
166   /// indicates state of derivatives of final results with respect to
167   /// distribution parameters s within resp_x_to_u_mapping() using the chain
168   /// rule df/dx dx/ds.  The default is to calculate derivatives with respect
169   /// to standard random variables u using the chain rule df/dx dx/du.
170   short distParamDerivs;
171 
172   /// boolean flag indicating use of distribution truncation for
173   /// defining global model bounds
174   bool truncatedBounds;
175   /// number of +/- standard deviations used for defining bounds truncation
176   Real boundVal;
177 
178   /// "primary" all continuous variable mapping indices flowed down
179   /// from higher level iteration
180   SizetArray primaryACVarMapIndices;
181   // "primary" all discrete int variable mapping indices flowed down from
182   // higher level iteration
183   //SizetArray primaryADIVarMapIndices;
184   // "primary" all discrete string variable mapping indices flowed down from
185   // higher level iteration
186   //SizetArray primaryADSVarMapIndices;
187   // "primary" all discrete real variable mapping indices flowed down from
188   // higher level iteration
189   //SizetArray primaryADRVarMapIndices;
190   /// "secondary" all continuous variable mapping targets flowed down
191   /// from higher level iteration
192   ShortArray secondaryACVarMapTargets;
193   // "secondary" all discrete int variable mapping targets flowed down
194   // from higher level iteration
195   //ShortArray secondaryADIVarMapTargets;
196   // "secondary" all discrete string variable mapping targets flowed down
197   // from higher level iteration
198   //ShortArray secondaryADSVarMapTargets;
199   // "secondary" all discrete real variable mapping targets flowed down
200   // from higher level iteration
201   //ShortArray secondaryADRVarMapTargets;
202 
203   /// static pointer to this class for use in static callbacks
204   static ProbabilityTransformModel* ptmInstance;
205 };
206 
207 
208 inline void ProbabilityTransformModel::
nested_variable_mappings(const SizetArray & c_index1,const SizetArray & di_index1,const SizetArray & ds_index1,const SizetArray & dr_index1,const ShortArray & c_target2,const ShortArray & di_target2,const ShortArray & ds_target2,const ShortArray & dr_target2)209 nested_variable_mappings(const SizetArray& c_index1,
210 			 const SizetArray& di_index1,
211 			 const SizetArray& ds_index1,
212 			 const SizetArray& dr_index1,
213 			 const ShortArray& c_target2,
214 			 const ShortArray& di_target2,
215 			 const ShortArray& ds_target2,
216 			 const ShortArray& dr_target2)
217 {
218   primaryACVarMapIndices      = c_index1;
219   //primaryADIVarMapIndices   = di_index1;
220   //primaryADSVarMapIndices   = ds_index1;
221   //primaryADRVarMapIndices   = dr_index1;
222   secondaryACVarMapTargets    = c_target2;
223   //secondaryADIVarMapTargets = di_target2;
224   //secondaryADSVarMapTargets = ds_target2;
225   //secondaryADRVarMapTargets = dr_target2;
226 }
227 
228 
229 inline short ProbabilityTransformModel::
query_distribution_parameter_derivatives() const230 query_distribution_parameter_derivatives() const
231 {
232   short dist_param_derivs = NO_DERIVS;
233   size_t i, num_outer_cv = secondaryACVarMapTargets.size();
234   if (num_outer_cv) {
235     bool tgt = false, no_tgt = false;
236     for (i=0; i<num_outer_cv; ++i)
237       if (secondaryACVarMapTargets[i] == Pecos::NO_TARGET) no_tgt = true;
238       else                                                    tgt = true;
239     if (tgt && no_tgt) dist_param_derivs = MIXED_DERIVS;
240     else if (tgt)      dist_param_derivs =   ALL_DERIVS;
241   }
242   return dist_param_derivs;
243 }
244 
245 
246 inline void ProbabilityTransformModel::
activate_distribution_parameter_derivatives()247 activate_distribution_parameter_derivatives()
248 { distParamDerivs = query_distribution_parameter_derivatives(); }
249 
250 
251 inline void ProbabilityTransformModel::
deactivate_distribution_parameter_derivatives()252 deactivate_distribution_parameter_derivatives()
253 { distParamDerivs = NO_DERIVS; }
254 
255 
nested_acv1_indices() const256 inline const SizetArray& ProbabilityTransformModel::nested_acv1_indices() const
257 { return primaryACVarMapIndices; }
258 
259 
nested_acv2_targets() const260 inline const ShortArray& ProbabilityTransformModel::nested_acv2_targets() const
261 { return secondaryACVarMapTargets; }
262 
263 
resize_pending() const264 inline bool ProbabilityTransformModel::resize_pending() const
265 { return subModel.resize_pending(); }
266 
267 
initialize_nataf()268 inline void ProbabilityTransformModel::initialize_nataf()
269 {
270   if (natafTransform.is_null()) {
271     natafTransform = Pecos::ProbabilityTransformation("nataf"); // for now
272     // shallow copies
273     natafTransform.x_distribution(subModel.multivariate_distribution());
274     natafTransform.u_distribution(mvDist);
275   }
276 }
277 
278 
update_transformation()279 inline void ProbabilityTransformModel::update_transformation()
280 {
281   mvDist.pull_distribution_parameters(subModel.multivariate_distribution());
282   // x-space correlations assigned in Model and u-space is uncorrelated
283   //update_distribution_correlations();
284 
285   // Modify the correlation matrix (Nataf) and compute its Cholesky factor.
286   // Since the uncertain variable distributions (means, std devs, correlations)
287   // may change, update of correlation warpings is performed regularly.
288   natafTransform.transform_correlations();
289 
290   update_model_bounds(truncatedBounds, boundVal);
291 }
292 
293 
294 inline void ProbabilityTransformModel::
initialize_transformation(short u_space_type)295 initialize_transformation(short u_space_type)
296 {
297   if (mvDist.is_null()) // already initialized: no current reason to update
298     mvDist = Pecos::MultivariateDistribution(Pecos::MARGINALS_CORRELATIONS);
299 
300   const Pecos::MultivariateDistribution& x_dist
301     = subModel.multivariate_distribution();
302   initialize_distribution_types(u_space_type, x_dist, mvDist);
303   initialize_nataf();
304   initialize_dakota_variable_types();
305   verify_correlation_support(u_space_type);
306 
307   // also perform run-time update as some construct-time operations require it
308   // (e.g. grid_size() for maxConcurrency, inverse vars transform in ctor).
309   update_transformation();
310 }
311 
312 
313 inline void ProbabilityTransformModel::
update_from_subordinate_model(size_t depth)314 update_from_subordinate_model(size_t depth)
315 {
316   // standard updates for RecastModels, including subModel recursion
317   //RecastModel::update_from_subordinate_model(depth);
318   // ordering problem with invMapping dependence on dist params
319 
320   // data flows from the bottom-up, so recurse first
321   if (depth == std::numeric_limits<size_t>::max())
322     subModel.update_from_subordinate_model(depth); // retain special value (inf)
323   else if (depth)
324     subModel.update_from_subordinate_model(depth - 1); // decrement
325   //else depth exhausted --> update this level only
326 
327   // propagate any subModel parameter updates to mvDist
328   update_transformation();
329 
330   // now pull additional updates from subModel (requires latest dist params)
331   RecastModel::update_from_model(subModel);
332 }
333 
334 
335 inline bool ProbabilityTransformModel::
nonlinear_variables_mapping(const Pecos::MultivariateDistribution & x_dist,const Pecos::MultivariateDistribution & u_dist) const336 nonlinear_variables_mapping(const Pecos::MultivariateDistribution& x_dist,
337 			    const Pecos::MultivariateDistribution& u_dist) const
338 {
339   bool nln_vars_map = false;
340   const ShortArray& x_types = x_dist.random_variable_types();
341   const ShortArray& u_types = u_dist.random_variable_types();
342   size_t i, num_types = std::min(x_types.size(), u_types.size());
343   const BitArray& active_v = x_dist.active_variables();
344   for (i=0; i<num_types; ++i)
345     if (active_v[i]) {
346       switch (u_types[i]) {
347       case Pecos::STD_NORMAL:
348 	nln_vars_map = (x_types[i] != Pecos::NORMAL); break;
349       case Pecos::STD_UNIFORM:
350 	switch (x_types[i]) {
351 	case Pecos::CONTINUOUS_RANGE: case Pecos::UNIFORM:
352 	case Pecos::HISTOGRAM_BIN:    case Pecos::CONTINUOUS_INTERVAL_UNCERTAIN:
353 	  break;
354 	default:  nln_vars_map = true;  break;
355 	}
356 	break;
357       case Pecos::STD_EXPONENTIAL:
358 	nln_vars_map = (x_types[i] != Pecos::EXPONENTIAL); break;
359       case Pecos::STD_BETA:  nln_vars_map = (x_types[i] != Pecos::BETA);  break;
360       case Pecos::STD_GAMMA: nln_vars_map = (x_types[i] != Pecos::GAMMA); break;
361       default:               nln_vars_map = (x_types[i] != u_types[i]);   break;
362       }
363 
364       if (nln_vars_map) break;
365     }
366 
367   return nln_vars_map;
368 }
369 
370 
371 inline Pecos::ProbabilityTransformation& ProbabilityTransformModel::
probability_transformation()372 probability_transformation()
373 { return natafTransform; }
374 
375 
assign_instance()376 inline void ProbabilityTransformModel::assign_instance()
377 { ptmInstance = this; }
378 
379 
380 /** Map the variables from iterator space (u) to simulation space (x). */
381 inline void ProbabilityTransformModel::
vars_u_to_x_mapping(const Variables & u_vars,Variables & x_vars)382 vars_u_to_x_mapping(const Variables& u_vars, Variables& x_vars)
383 {
384   ptmInstance->natafTransform.trans_U_to_X(u_vars.continuous_variables(),
385 					   x_vars.continuous_variables_view());
386   // *** TO DO: active discrete {int,string,real}
387 }
388 
389 
390 /** Map the variables from simulation space (x) to iterator space (u). */
391 inline void ProbabilityTransformModel::
vars_x_to_u_mapping(const Variables & x_vars,Variables & u_vars)392 vars_x_to_u_mapping(const Variables& x_vars, Variables& u_vars)
393 {
394   ptmInstance->natafTransform.trans_X_to_U(x_vars.continuous_variables(),
395 					   u_vars.continuous_variables_view());
396   // *** TO DO: active discrete {int,string,real}
397 }
398 
399 
400 inline void ProbabilityTransformModel::
trans_grad_X_to_U(const RealVector & fn_grad_x,RealVector & fn_grad_u,const RealVector & x_vars)401 trans_grad_X_to_U(const RealVector& fn_grad_x, RealVector& fn_grad_u,
402 		  const RealVector& x_vars)
403 {
404   SizetMultiArrayConstView cv_ids = currentVariables.continuous_variable_ids();
405   SizetArray x_dvv; copy_data(cv_ids, x_dvv);
406   natafTransform.trans_grad_X_to_U(fn_grad_x, fn_grad_u, x_vars, x_dvv, cv_ids);
407 }
408 
409 
410 inline void ProbabilityTransformModel::
trans_grad_U_to_X(const RealVector & fn_grad_u,RealVector & fn_grad_x,const RealVector & x_vars)411 trans_grad_U_to_X(const RealVector& fn_grad_u, RealVector& fn_grad_x,
412 		  const RealVector& x_vars)
413 {
414   SizetMultiArrayConstView cv_ids = currentVariables.continuous_variable_ids();
415   SizetArray x_dvv; copy_data(cv_ids, x_dvv);
416   natafTransform.trans_grad_U_to_X(fn_grad_u, fn_grad_x, x_vars, x_dvv, cv_ids);
417 }
418 
419 
420 inline void ProbabilityTransformModel::
trans_grad_X_to_S(const RealVector & fn_grad_x,RealVector & fn_grad_s,const RealVector & x_vars)421 trans_grad_X_to_S(const RealVector& fn_grad_x, RealVector& fn_grad_s,
422 		  const RealVector& x_vars)
423 {
424   SizetMultiArrayConstView cv_ids = currentVariables.continuous_variable_ids();
425   SizetArray x_dvv; copy_data(cv_ids, x_dvv);
426   natafTransform.trans_grad_X_to_S(fn_grad_x, fn_grad_s, x_vars, x_dvv, cv_ids,
427     currentVariables.all_continuous_variable_ids(), primaryACVarMapIndices,
428     secondaryACVarMapTargets);
429 }
430 
431 
432 inline void ProbabilityTransformModel::
trans_hess_X_to_U(const RealSymMatrix & fn_hess_x,RealSymMatrix & fn_hess_u,const RealVector & x_vars,const RealVector & fn_grad_x)433 trans_hess_X_to_U(const RealSymMatrix& fn_hess_x, RealSymMatrix& fn_hess_u,
434 		  const RealVector& x_vars, const RealVector& fn_grad_x)
435 {
436   SizetMultiArrayConstView cv_ids = currentVariables.continuous_variable_ids();
437   SizetArray x_dvv; copy_data(cv_ids, x_dvv);
438   natafTransform.trans_hess_X_to_U(fn_hess_x, fn_hess_u, x_vars, fn_grad_x,
439 				   x_dvv, cv_ids);
440 }
441 
442 
443 inline size_t ProbabilityTransformModel::
rv_index_to_corr_index(size_t rv_index)444 rv_index_to_corr_index(size_t rv_index)
445 {
446   const Pecos::BitArray& active_corr
447     = subModel.multivariate_distribution().active_correlations();
448   if (active_corr.empty())
449     return rv_index; // no mask
450   else if (active_corr[rv_index]) { // offset RV index to account for mask
451     size_t i, corr_index = 0;
452     for (i=0; i<rv_index; ++i)
453       if (active_corr[i])
454 	++corr_index;
455     return corr_index;
456   }
457   else // RV not active in correlated variable subset
458     return _NPOS;
459 }
460 
461 
462 inline size_t ProbabilityTransformModel::
acv_index_to_corr_index(size_t acv_index)463 acv_index_to_corr_index(size_t acv_index)
464 {
465   const SharedVariablesData& svd = subModel.current_variables().shared_data();
466   return rv_index_to_corr_index(svd.acv_index_to_all_index(acv_index));
467 }
468 
469 } // namespace Dakota
470 
471 #endif
472