1 /* _______________________________________________________________________
2
3 PECOS: Parallel Environment for Creation Of Stochastics
4 Copyright (c) 2011, Sandia National Laboratories.
5 This software is distributed under the GNU Lesser General Public License.
6 For more information, see the README file in the top Pecos directory.
7 _______________________________________________________________________ */
8
9 #ifndef MARGINALS_CORR_DISTRIBUTION_HPP
10 #define MARGINALS_CORR_DISTRIBUTION_HPP
11
12 #include "MultivariateDistribution.hpp"
13 #include "RandomVariable.hpp"
14
15
16 namespace Pecos {
17
18
19 /// Class for multivariate distribution based on 1D marginals + correlations
20
21 /** The MarginalsCorrDistribution models a multivariate random variable
22 distribution using 1D marginals plus a correlation matrix. */
23
24 class MarginalsCorrDistribution: public MultivariateDistribution
25 {
26 public:
27
28 //
29 //- Heading: Constructors and destructor
30 //
31
32 MarginalsCorrDistribution(); ///< constructor
33 ~MarginalsCorrDistribution(); ///< destructor
34
35 //
36 //- Heading: Virtual function redefinitions
37 //
38
39 /// return randomVars
40 const std::vector<RandomVariable>& random_variables() const;
41 /// return randomVars
42 std::vector<RandomVariable>& random_variables();
43
44 /// return randomVars[i]
45 const RandomVariable& random_variable(size_t i) const;
46 /// return randomVars[i]
47 RandomVariable& random_variable(size_t i);
48
49 /// return ranVarTypes
50 const ShortArray& random_variable_types() const;
51 /// set ranVarTypes
52 void random_variable_types(const ShortArray& rv_types);
53
54 /// return ranVarTypes[i]
55 short random_variable_type(size_t i) const;
56 /// set ranVarTypes[i]
57 void random_variable_type(short rv_type, size_t i);
58
59 // BMA TODO: Review why these don't follow typical letter/envelope
60 // pattern. Left them for now...
61
62 /// pull non-standardized distribution parameters from pull_mvd
63 void pull_distribution_parameters(const MultivariateDistribution& pull_mvd);
64 /// pull non-standardized distribution parameters from pull_mvd_rep
65 void pull_distribution_parameters
66 (const std::shared_ptr<MultivariateDistribution> pull_mvd_rep);
67
68 /// pull non-standardized distribution parameters from pull_mvd, aligning
69 /// variables based on label lookups
70 void pull_distribution_parameters(const MultivariateDistribution& pull_mvd,
71 const StringArray& pull_labels, const StringArray& push_labels);
72 /// pull non-standardized distribution parameters from pull_mvd_rep, aligning
73 /// variables based on label lookups
74 void pull_distribution_parameters
75 (const std::shared_ptr<MultivariateDistribution> pull_mvd_rep,
76 const StringArray& pull_labels, const StringArray& push_labels);
77
78 /// pull non-standardized distribution parameters for a particular
79 /// random variable from pull_mvd
80 void pull_distribution_parameters(const MultivariateDistribution& pull_mvd,
81 size_t pull_index, size_t push_index);
82 /// pull non-standardized distribution parameters for a particular
83 /// random variable from pull_mvd_rep
84 void pull_distribution_parameters
85 (const std::shared_ptr<MultivariateDistribution> pull_mvd_rep,
86 size_t pull_index, size_t push_index);
87
88 /// return activeVars
89 const BitArray& active_variables() const;
90 /// set activeVars
91 void active_variables(const BitArray& );
92 /// return activeCorr
93 const BitArray& active_correlations() const;
94 /// set activeCorr
95 void active_correlations(const BitArray& );
96
97 /// return corrMatrix
98 const RealSymMatrix& correlation_matrix() const;
99 /// set corrMatrix
100 void correlation_matrix(const RealSymMatrix& corr);
101
102 /// assemble means and standard deviations from RandomVariable::moments()
103 RealRealPairArray moments() const;
104 /// assemble means from RandomVariable::mean()
105 RealVector means() const;
106 /// assemble standard deviations from RandomVariable::standard_deviation()
107 RealVector std_deviations() const;
108 /// assemble variances from RandomVariable::variance()
109 RealVector variances() const;
110
111 /// assemble distribution lower and upper bounds from RandomVariable::bounds()
112 RealRealPairArray distribution_bounds() const;
113 /// assemble distribution lower bounds from RandomVariable::bounds()
114 RealVector distribution_lower_bounds() const;
115 /// assemble distribution upper bounds from RandomVariable::bounds()
116 RealVector distribution_upper_bounds() const;
117
118 // these bounds are distinct from distribution bounds: for Dakota, they are
119 // "global" bounds; in Pecos, they include RangeVariable bounds
120
121 bool global_bounds() const;
122
123 /// check incoming vec for length equal to number of active random variables
124 template <typename OrdinalType, typename ScalarType>
125 void check_active_length(
126 const Teuchos::SerialDenseVector<OrdinalType, ScalarType>& vec,
127 const BitArray& mask) const;
128
129 /// set real lower bounds for variable ranges
130 void lower_bounds(const RealVector& l_bnds,
131 const BitArray& mask = BitArray());
132 /// set integer lower bounds for variable ranges
133 void lower_bounds(const IntVector& l_bnds,
134 const BitArray& mask = BitArray());
135 /// set real lower bound for a variable range
136 void lower_bound(Real l_bnd, size_t rv_index);
137 /// set int lower bound for a variable range
138 void lower_bound(int l_bnd, size_t rv_index);
139
140 /// set real upper bounds for variable ranges
141 void upper_bounds(const RealVector& u_bnds,
142 const BitArray& mask = BitArray());
143 /// set int upper bounds for variable ranges
144 void upper_bounds(const IntVector& u_bnds,
145 const BitArray& mask = BitArray());
146 /// set real upper bound for a variable range
147 void upper_bound(Real u_bnd, size_t rv_index);
148 /// set int upper bound for a variable range
149 void upper_bound(int u_bnd, size_t rv_index);
150
151 /* This requires clients to bypass the base class and cast a rep pointer
152 template <typename OrdinalType, typename ScalarType>
153 void lower_bounds(
154 const Teuchos::SerialDenseVector<OrdinalType, ScalarType>& l_bnds,
155 const BitArray& mask = BitArray());
156 template <typename ScalarType>
157 void lower_bound(ScalarType l_bnd, size_t rv_index);
158
159 template <typename OrdinalType, typename ScalarType>
160 void upper_bounds(
161 const Teuchos::SerialDenseVector<OrdinalType, ScalarType>& u_bnds,
162 const BitArray& mask = BitArray());
163 template <typename ScalarType>
164 void upper_bound(ScalarType u_bnd, size_t rv_index);
165 */
166
167 //
168 //- Heading: Member function definitions
169 //
170
171 /// set ranVarTypes and initialize randomVars
172 void initialize_types(const ShortArray& rv_types,
173 const BitArray& active_vars = BitArray());
174 /// assigns corrMatrix and activeCorr; invokes initialize_correlations()
175 void initialize_correlations(const RealSymMatrix& corr,
176 const BitArray& active_corr = BitArray());
177 /// initializes correlationFlag and performs sanity checks
178 void initialize_correlations();
179
180 /// update a scalar distribution parameter within randomVars[v]
181 template <typename ValueType>
182 void push_parameter(size_t v, short dist_param, ValueType value);
183 /// update values for one distribution parameter across a sequence
184 /// of random variables
185 template <typename OrdinalType, typename ScalarType>
186 void push_parameters(size_t start_v, size_t num_v, short dist_param,
187 const Teuchos::SerialDenseVector<OrdinalType, ScalarType>& values);
188 /// update values for one distribution parameter across the set
189 /// of random variables with matching RV type
190 template <typename OrdinalType, typename ScalarType>
191 void push_parameters(short rv_type, short dist_param,
192 const Teuchos::SerialDenseVector<OrdinalType, ScalarType>& values);
193 /// update values for one distribution parameter across a sequence
194 /// of random variables
195 template <typename ValueType>
196 void push_parameters(size_t start_v, size_t num_v, short dist_param,
197 const std::vector<ValueType>& values);
198 /// update values for one distribution parameter across the set
199 /// of random variables with matching RV type
200 template <typename ValueType>
201 void push_parameters(short rv_type, short dist_param,
202 const std::vector<ValueType>& values);
203
204 /// update val from a scalar distribution parameter from randomVars[v]
205 template <typename ValueType>
206 void pull_parameter(size_t v, short dist_param, ValueType& val) const;
207 /// define array of values using the identified distribution parameter
208 /// across the specified range of random variables
209 template <typename ValueType>
210 void pull_parameters(size_t start_v, size_t num_v, short dist_param,
211 std::vector<ValueType>& values) const;
212 /// define array of values using the identified distribution parameter
213 /// across the specified range of random variables
214 template <typename ValueType>
215 void pull_parameters(short rv_type, short dist_param,
216 std::vector<ValueType>& values) const;
217 /// update values for one distribution parameter across the set
218 /// of random variables with matching RV type
219 template <typename OrdinalType, typename ScalarType>
220 void pull_parameters(short rv_type, short dist_param,
221 Teuchos::SerialDenseVector<OrdinalType, ScalarType>& values) const;
222
223 /// return a distribution parameter from randomVars[v]
224 template <typename ValueType>
225 ValueType pull_parameter(size_t v, short dist_param) const;
226 /// return array of values for one distribution parameter across the set
227 /// of random variables with matching RV type
228 template <typename ValueType>
229 std::vector<ValueType> pull_parameters(short rv_type, short dist_param) const;
230
231 /// return the size of (non-scalar) distribution data from randomVars[v]
232 template <typename ValueType>
233 size_t pull_parameter_size(size_t v, short dist_param) const;
234
235 /*
236 /// expand corrMatrix from probabilistic variables to combined variables
237 void expand_correlation_matrix(size_t num_lead_v, size_t num_prob_v,
238 size_t num_trail_v);
239 /// contract corrMatrix from combined variables to probabilistic variables
240 void contract_correlation_matrix(size_t num_lead_v, size_t num_prob_v,
241 size_t num_trail_v);
242 */
243
244 /// set globalBndsFlag based on ranVarTypes
245 void check_global_bounds();
246 /// verify that randomVars[i].type() equals rv_type
247 void check_random_variable_type(size_t i, short rv_type) const;
248
249 /// draw a sample from the i-th RandomVariable
250 template <typename Engine>
251 Real draw_sample(size_t i, Engine& rng) const;
252 /// draw a sample from the i-th standardized RandomVariable
253 template <typename Engine>
254 Real draw_standard_sample(size_t i, Engine& rng) const;
255
256 protected:
257
258 //
259 //- Heading: Virtual function redefinitions
260 //
261
262 Real pdf(const RealVector& pt) const;
263 Real log_pdf(const RealVector& pt) const;
264
265 Real pdf(Real val, size_t rv_index) const;
266 Real pdf_gradient(Real val, size_t rv_index) const;
267 Real pdf_hessian(Real val, size_t rv_index) const;
268 Real log_pdf(Real val, size_t rv_index) const;
269 Real log_pdf_gradient(Real val, size_t rv_index) const;
270 Real log_pdf_hessian(Real val, size_t rv_index) const;
271
272 /// copy marginals + correlation data between representations
273 void copy_rep(std::shared_ptr<MultivariateDistribution> mvd_rep);
274
275 //
276 //- Heading: Data
277 //
278
279 /// vector of types of each u-space standardized uncertain variable to
280 /// which each x-space variable is transformed
281 ShortArray ranVarTypes;
282 /// vector of random variables encapsulating distribution parameters and
283 /// statistical functions (pdf, cdf, etc.)
284 std::vector<RandomVariable> randomVars;
285 /// subset of randomVars that are currently active (if empty, then
286 /// no subset: all variables are active)
287 BitArray activeVars;
288
289 /// matrix of random variable correlation coefficients
290 RealSymMatrix corrMatrix;
291 /// subset of randomVars to which the corrMatrix refers (if empty,
292 /// then no subset: correlations are provided for all variables)
293 BitArray activeCorr;
294
295 /// precompute the presence of random variables that support global bounds,
296 /// for fto accelerate run-time operations
297 bool globalBndsFlag;
298
299 private:
300
301 //
302 //- Heading: Data
303 //
304
305 };
306
307
MarginalsCorrDistribution()308 inline MarginalsCorrDistribution::MarginalsCorrDistribution():
309 MultivariateDistribution(BaseConstructor()), globalBndsFlag(false)
310 { }
311
312
~MarginalsCorrDistribution()313 inline MarginalsCorrDistribution::~MarginalsCorrDistribution()
314 { }
315
316
317 inline const std::vector<RandomVariable>& MarginalsCorrDistribution::
random_variables() const318 random_variables() const
319 { return randomVars; }
320
321
322 inline std::vector<RandomVariable>& MarginalsCorrDistribution::
random_variables()323 random_variables()
324 { return randomVars; }
325
326
327 inline const RandomVariable& MarginalsCorrDistribution::
random_variable(size_t rv_index) const328 random_variable(size_t rv_index) const
329 { return randomVars[rv_index]; }
330
331
332 inline RandomVariable& MarginalsCorrDistribution::
random_variable(size_t rv_index)333 random_variable(size_t rv_index)
334 { return randomVars[rv_index]; }
335
336
check_global_bounds()337 inline void MarginalsCorrDistribution::check_global_bounds()
338 {
339 // As a first cut, identify range variables (design, state) as having
340 // "global" bounds, which need to synchronize from global parameter space
341 // updates, separate from distribution parameter updates
342 globalBndsFlag = false;
343 size_t i, num_rv = ranVarTypes.size();
344 for (i=0; i<num_rv; ++i)
345 if (ranVarTypes[i] == CONTINUOUS_RANGE || ranVarTypes[i] == DISCRETE_RANGE)
346 { globalBndsFlag = true; break; }
347 }
348
349
350 inline void MarginalsCorrDistribution::
check_random_variable_type(size_t rv_index,short rv_type) const351 check_random_variable_type(size_t rv_index, short rv_type) const
352 {
353 if (randomVars[rv_index].type() != rv_type) {
354 PCerr << "Error: inconsistent random variable type in MarginalsCorr"
355 << "Distribution::check_random_variable_type()." << std::endl;
356 abort_handler(-1);
357 }
358 }
359
360
361 inline const ShortArray& MarginalsCorrDistribution::
random_variable_types() const362 random_variable_types() const
363 { return ranVarTypes; }
364
365
366 inline void MarginalsCorrDistribution::
random_variable_types(const ShortArray & rv_types)367 random_variable_types(const ShortArray& rv_types)
368 {
369 ranVarTypes = rv_types;
370 check_global_bounds();
371 }
372
373
374 inline short MarginalsCorrDistribution::
random_variable_type(size_t rv_index) const375 random_variable_type(size_t rv_index) const
376 { return ranVarTypes[rv_index]; }
377
378
379 inline void MarginalsCorrDistribution::
random_variable_type(short rv_type,size_t rv_index)380 random_variable_type(short rv_type, size_t rv_index)
381 {
382 bool new_rv_global = (rv_type == CONTINUOUS_RANGE ||
383 rv_type == DISCRETE_RANGE);
384 if (globalBndsFlag) {
385 short& curr_rv_type = ranVarTypes[rv_index];
386 bool curr_rv_global = (curr_rv_type == CONTINUOUS_RANGE ||
387 curr_rv_type == DISCRETE_RANGE);
388 curr_rv_type = rv_type;
389 if (curr_rv_global && !new_rv_global) // previous global type removed
390 check_global_bounds(); // --> check for flag update
391 }
392 else { // no current types define global bnds --> update globalBndsFlag
393 // to true if incoming rv_type defines global bnds
394 ranVarTypes[rv_index] = rv_type;
395 globalBndsFlag = new_rv_global;
396 }
397 }
398
399
active_variables() const400 inline const BitArray& MarginalsCorrDistribution::active_variables() const
401 { return activeVars; }
402
403
404 inline void MarginalsCorrDistribution::
active_variables(const BitArray & active_vars)405 active_variables(const BitArray& active_vars)
406 { activeVars = active_vars; }
407
408
active_correlations() const409 inline const BitArray& MarginalsCorrDistribution::active_correlations() const
410 { return activeCorr; }
411
412
413 inline void MarginalsCorrDistribution::
active_correlations(const BitArray & active_corr)414 active_correlations(const BitArray& active_corr)
415 { activeCorr = active_corr; }
416
417
418 inline const RealSymMatrix& MarginalsCorrDistribution::
correlation_matrix() const419 correlation_matrix() const
420 { return corrMatrix; }
421
422
423 inline void MarginalsCorrDistribution::
correlation_matrix(const RealSymMatrix & corr)424 correlation_matrix(const RealSymMatrix& corr)
425 { corrMatrix = corr; }
426
427
428 //inline const RealMatrix& MarginalsCorrDistribution::correlation_factor() const
429 //{ return corrCholeskyFactor; }
430
431
432 /** For consistent random variable ordering. */
433 inline void MarginalsCorrDistribution::
pull_distribution_parameters(const std::shared_ptr<MultivariateDistribution> pull_mvd_rep)434 pull_distribution_parameters
435 (const std::shared_ptr<MultivariateDistribution> pull_mvd_rep)
436 {
437 size_t v, num_rv = ranVarTypes.size();
438 for (v=0; v<num_rv; ++v)
439 pull_distribution_parameters(pull_mvd_rep, v, v);
440 }
441
442
443 inline void MarginalsCorrDistribution::
pull_distribution_parameters(const MultivariateDistribution & pull_mvd)444 pull_distribution_parameters(const MultivariateDistribution& pull_mvd)
445 { pull_distribution_parameters(pull_mvd.multivar_dist_rep()); }
446
447
448 /** For potentially inconsistent random variable ordering that requires
449 a lookup. */
450 inline void MarginalsCorrDistribution::
pull_distribution_parameters(const std::shared_ptr<MultivariateDistribution> pull_mvd_rep,const StringArray & pull_labels,const StringArray & push_labels)451 pull_distribution_parameters
452 (const std::shared_ptr<MultivariateDistribution> pull_mvd_rep,
453 const StringArray& pull_labels, const StringArray& push_labels)
454 {
455 size_t v, num_rv = ranVarTypes.size();
456 for (v=0; v<num_rv; ++v) {
457 size_t push_index = find_index(push_labels, pull_labels[v]);
458 if (push_index != _NPOS)
459 pull_distribution_parameters(pull_mvd_rep, v, push_index);
460 }
461 }
462
463
464 inline void MarginalsCorrDistribution::
pull_distribution_parameters(const MultivariateDistribution & pull_mvd,const StringArray & pull_labels,const StringArray & push_labels)465 pull_distribution_parameters(const MultivariateDistribution& pull_mvd,
466 const StringArray& pull_labels,
467 const StringArray& push_labels)
468 {
469 pull_distribution_parameters(pull_mvd.multivar_dist_rep(),
470 pull_labels, push_labels);
471 }
472
473
474 inline void MarginalsCorrDistribution::
pull_distribution_parameters(const MultivariateDistribution & pull_mvd,size_t pull_index,size_t push_index)475 pull_distribution_parameters(const MultivariateDistribution& pull_mvd,
476 size_t pull_index, size_t push_index)
477 {
478 pull_distribution_parameters(pull_mvd.multivar_dist_rep(),
479 pull_index, push_index);
480 }
481
482
483 template <typename ValueType>
484 void MarginalsCorrDistribution::
push_parameter(size_t v,short dist_param,ValueType value)485 push_parameter(size_t v, short dist_param, ValueType value)
486 { randomVars[v].push_parameter(dist_param, value); }
487
488
489 template <typename OrdinalType, typename ScalarType>
490 void MarginalsCorrDistribution::
push_parameters(size_t start_v,size_t num_v,short dist_param,const Teuchos::SerialDenseVector<OrdinalType,ScalarType> & values)491 push_parameters(size_t start_v, size_t num_v, short dist_param,
492 const Teuchos::SerialDenseVector<OrdinalType, ScalarType>& values)
493 {
494 // set one type of distribution parameter for a range of random variables
495 // TO DO: would like to retire this version if Dakota migrates from Teuchos
496 // to std::vector for dist params
497
498 size_t i, v, num_updates = std::min((size_t)values.length(), num_v);
499 for (i=0, v=start_v; i<num_updates; ++i, ++v)
500 randomVars[v].push_parameter(dist_param, values[i]);
501 }
502
503
504 template <typename ValueType>
505 void MarginalsCorrDistribution::
push_parameters(size_t start_v,size_t num_v,short dist_param,const std::vector<ValueType> & values)506 push_parameters(size_t start_v, size_t num_v, short dist_param,
507 const std::vector<ValueType>& values)
508 {
509 // set one distribution parameter type for a range of random variables
510
511 size_t i, v, num_updates = std::min(values.size(), num_v);
512 for (i=0, v=start_v; i<num_updates; ++i, ++v)
513 randomVars[v].push_parameter(dist_param, values[i]);
514 }
515
516
517 template <typename OrdinalType, typename ScalarType>
518 void MarginalsCorrDistribution::
push_parameters(short rv_type,short dist_param,const Teuchos::SerialDenseVector<OrdinalType,ScalarType> & values)519 push_parameters(short rv_type, short dist_param,
520 const Teuchos::SerialDenseVector<OrdinalType, ScalarType>& values)
521 {
522 // rv_type eliminates need to check for dist_param support
523 // TO DO: would like to retire this version if Dakota migrates from Teuchos
524 // to std::vector for dist params
525
526 size_t v, num_v = ranVarTypes.size(), cntr = 0, num_vals = values.length();
527 for (v=0; v < num_v && cntr < num_vals; ++v)
528 if (ranVarTypes[v] == rv_type)
529 randomVars[v].push_parameter(dist_param, values[cntr++]);
530 }
531
532
533 template <typename ValueType>
534 void MarginalsCorrDistribution::
push_parameters(short rv_type,short dist_param,const std::vector<ValueType> & values)535 push_parameters(short rv_type, short dist_param,
536 const std::vector<ValueType>& values)
537 {
538 // rv_type eliminates need to check for dist_param support
539
540 size_t v, num_v = ranVarTypes.size(), cntr = 0, num_vals = values.size();
541 for (v=0; v < num_v && cntr < num_vals; ++v)
542 if (ranVarTypes[v] == rv_type)
543 randomVars[v].push_parameter(dist_param, values[cntr++]);
544 }
545
546
547 template <typename ValueType>
548 void MarginalsCorrDistribution::
pull_parameter(size_t v,short dist_param,ValueType & val) const549 pull_parameter(size_t v, short dist_param, ValueType& val) const
550 { randomVars[v].pull_parameter(dist_param, val); }
551
552
553 template <typename ValueType>
554 void MarginalsCorrDistribution::
pull_parameters(size_t start_v,size_t num_v,short dist_param,std::vector<ValueType> & values) const555 pull_parameters(size_t start_v, size_t num_v, short dist_param,
556 std::vector<ValueType>& values) const
557 {
558 // set one distribution parameter type for a range of random variables
559
560 size_t i, v;
561 values.resize(num_v);
562 for (i=0, v=start_v; i<num_v; ++i, ++v)
563 randomVars[v].pull_parameter(dist_param, values[i]);
564 }
565
566
567 template <typename ValueType>
568 void MarginalsCorrDistribution::
pull_parameters(short rv_type,short dist_param,std::vector<ValueType> & values) const569 pull_parameters(short rv_type, short dist_param,
570 std::vector<ValueType>& values) const
571 {
572 // rv_type eliminates need to check for dist_param support
573
574 values.resize(std::count(ranVarTypes.begin(), ranVarTypes.end(), rv_type));
575
576 size_t v, num_v = ranVarTypes.size(), cntr = 0;
577 for (v=0; v<num_v; ++v)
578 if (ranVarTypes[v] == rv_type)
579 randomVars[v].pull_parameter(dist_param, values[cntr++]);
580 }
581
582
583 template <typename OrdinalType, typename ScalarType>
584 void MarginalsCorrDistribution::
pull_parameters(short rv_type,short dist_param,Teuchos::SerialDenseVector<OrdinalType,ScalarType> & values) const585 pull_parameters(short rv_type, short dist_param,
586 Teuchos::SerialDenseVector<OrdinalType, ScalarType>& values) const
587 {
588 // rv_type eliminates need to check for dist_param support
589
590 size_t v, num_v = ranVarTypes.size(), cntr = 0;
591 values.sizeUninitialized(
592 std::count(ranVarTypes.begin(), ranVarTypes.end(), rv_type));
593
594 for (v=0; v<num_v; ++v)
595 if (ranVarTypes[v] == rv_type)
596 randomVars[v].pull_parameter(dist_param, values[cntr++]);
597 }
598
599
600 template <typename ValueType>
601 ValueType MarginalsCorrDistribution::
pull_parameter(size_t v,short dist_param) const602 pull_parameter(size_t v, short dist_param) const
603 {
604 ValueType val;
605 randomVars[v].pull_parameter(dist_param, val);
606 return val;
607 }
608
609
610 template <typename ValueType>
611 std::vector<ValueType> MarginalsCorrDistribution::
pull_parameters(short rv_type,short dist_param) const612 pull_parameters(short rv_type, short dist_param) const
613 {
614 std::vector<ValueType> vals;
615 pull_parameters<ValueType>(rv_type, dist_param, vals);
616 return vals;
617 }
618
619
620 template <typename ValueType>
621 size_t MarginalsCorrDistribution::
pull_parameter_size(size_t v,short dist_param) const622 pull_parameter_size(size_t v, short dist_param) const
623 {
624 // overhead of one array copy instead of two or more
625 ValueType val;
626 randomVars[v].pull_parameter(dist_param, val);
627 return val.size();
628 }
629
630
631 /* These APIs are not currently needed but could be useful in the future:
632 template <typename VectorType>
633 void MarginalsCorrDistribution::
634 push_parameters(size_t v, const ShortArray& dist_params,
635 const VectorType& values)
636 {
637 // set multiple distribution parameters for a single variable
638
639 RandomVariable& random_var = randomVars[v];
640 size_t i, num_params = std::min(dist_params.size(), values.length());
641 for (i=0; i<num_params; ++i)
642 random_var.push_parameter(dist_params[i], values[i]);
643 }
644
645
646 template <typename VectorType>
647 void MarginalsCorrDistribution::
648 push_parameters(short dist_param, const VectorType& values)
649 {
650 // Without rv_type, query RV for dist_param support to match values to RV
651
652 size_t v, num_v = randomVars.size(), cntr = 0, num_vals = values.length();
653 for (v=0; v < num_v && cntr < num_vals; ++v)
654 if (randomVars[v].supports(dist_param))
655 randomVars[v].push_parameter(dist_param, values[cntr++]);
656 }
657 */
658
659
moments() const660 inline RealRealPairArray MarginalsCorrDistribution::moments() const
661 {
662 size_t v, num_v = randomVars.size();
663 RealRealPairArray mom;
664 if (activeVars.empty()) {
665 mom.resize(num_v);
666 for (v=0; v<num_v; ++v)
667 mom[v] = randomVars[v].moments();
668 }
669 else {
670 mom.resize(activeVars.count());
671 size_t av_cntr = 0;
672 for (v=0; v<num_v; ++v)
673 if (activeVars[v])
674 mom[av_cntr++] = randomVars[v].moments();
675 }
676 return mom;
677 }
678
679
means() const680 inline RealVector MarginalsCorrDistribution::means() const
681 {
682 size_t v, num_v = randomVars.size();
683 RealVector means;
684 if (activeVars.empty()) {
685 means.sizeUninitialized(num_v);
686 for (v=0; v<num_v; ++v)
687 means[v] = randomVars[v].mean();
688 }
689 else {
690 means.sizeUninitialized(activeVars.count());
691 size_t av_cntr = 0;
692 for (v=0; v<num_v; ++v)
693 if (activeVars[v])
694 means[av_cntr++] = randomVars[v].mean();
695 }
696 return means;
697 }
698
699
std_deviations() const700 inline RealVector MarginalsCorrDistribution::std_deviations() const
701 {
702 size_t v, num_v = randomVars.size();
703 RealVector std_devs;
704 if (activeVars.empty()) {
705 std_devs.sizeUninitialized(num_v);
706 for (v=0; v<num_v; ++v)
707 std_devs[v] = randomVars[v].standard_deviation();
708 }
709 else {
710 std_devs.sizeUninitialized(activeVars.count());
711 size_t av_cntr = 0;
712 for (v=0; v<num_v; ++v)
713 if (activeVars[v])
714 std_devs[av_cntr++] = randomVars[v].standard_deviation();
715 }
716 return std_devs;
717 }
718
719
variances() const720 inline RealVector MarginalsCorrDistribution::variances() const
721 {
722 size_t v, num_v = randomVars.size();
723 RealVector vars;
724 if (activeVars.empty()) {
725 vars.sizeUninitialized(num_v);
726 for (v=0; v<num_v; ++v)
727 vars[v] = randomVars[v].variance();
728 }
729 else {
730 vars.sizeUninitialized(activeVars.count());
731 size_t av_cntr = 0;
732 for (v=0; v<num_v; ++v)
733 if (activeVars[v])
734 vars[av_cntr++] = randomVars[v].variance();
735 }
736 return vars;
737 }
738
739
distribution_bounds() const740 inline RealRealPairArray MarginalsCorrDistribution::distribution_bounds() const
741 {
742 size_t v, num_v = randomVars.size();
743 RealRealPairArray bnds;
744 if (activeVars.empty()) {
745 bnds.resize(num_v);
746 for (v=0; v<num_v; ++v)
747 bnds[v] = randomVars[v].distribution_bounds();
748 }
749 else {
750 bnds.resize(activeVars.count());
751 size_t av_cntr = 0;
752 for (v=0; v<num_v; ++v)
753 if (activeVars[v])
754 bnds[av_cntr++] = randomVars[v].distribution_bounds();
755 }
756 return bnds;
757 }
758
759
distribution_lower_bounds() const760 inline RealVector MarginalsCorrDistribution::distribution_lower_bounds() const
761 {
762 size_t v, num_v = randomVars.size();
763 RealVector lwr_bnds;
764 if (activeVars.empty()) {
765 lwr_bnds.sizeUninitialized(num_v);
766 for (v=0; v<num_v; ++v)
767 lwr_bnds[v] = randomVars[v].distribution_bounds().first;
768 }
769 else {
770 lwr_bnds.sizeUninitialized(activeVars.count());
771 size_t av_cntr = 0;
772 for (v=0; v<num_v; ++v)
773 if (activeVars[v])
774 lwr_bnds[av_cntr++] = randomVars[v].distribution_bounds().first;
775 }
776 return lwr_bnds;
777 }
778
779
distribution_upper_bounds() const780 inline RealVector MarginalsCorrDistribution::distribution_upper_bounds() const
781 {
782 size_t v, num_v = randomVars.size();
783 RealVector upr_bnds;
784 if (activeVars.empty()) {
785 upr_bnds.sizeUninitialized(num_v);
786 for (v=0; v<num_v; ++v)
787 upr_bnds[v] = randomVars[v].distribution_bounds().second;
788 }
789 else {
790 upr_bnds.sizeUninitialized(activeVars.count());
791 size_t av_cntr = 0;
792 for (v=0; v<num_v; ++v)
793 if (activeVars[v])
794 upr_bnds[av_cntr++] = randomVars[v].distribution_bounds().second;
795 }
796 return upr_bnds;
797 }
798
799
global_bounds() const800 inline bool MarginalsCorrDistribution::global_bounds() const
801 { return globalBndsFlag; }
802
803
804 template <typename OrdinalType, typename ScalarType>
check_active_length(const Teuchos::SerialDenseVector<OrdinalType,ScalarType> & vec,const BitArray & mask) const805 void MarginalsCorrDistribution::check_active_length(
806 const Teuchos::SerialDenseVector<OrdinalType, ScalarType>& vec,
807 const BitArray& mask) const
808 {
809 size_t vec_len = vec.length(),
810 expect_len = (mask.empty()) ? randomVars.size() : mask.count();
811 if (vec_len != expect_len) {
812 PCerr << "Error: bad active vector length (" << vec_len << "); "
813 << expect_len << "expected." << std::endl;
814 abort_handler(-1);
815 }
816 }
817
818
819 inline void MarginalsCorrDistribution::
lower_bounds(const RealVector & l_bnds,const BitArray & mask)820 lower_bounds(const RealVector& l_bnds, const BitArray& mask)
821 {
822 check_active_length(l_bnds, mask);
823 size_t v, num_v = randomVars.size();
824 if (mask.empty())
825 for (v=0; v<num_v; ++v)
826 randomVars[v].lower_bound(l_bnds[v]);
827 else {
828 size_t av_cntr = 0;
829 for (v=0; v<num_v; ++v)
830 if (mask[v])
831 randomVars[v].lower_bound(l_bnds[av_cntr++]);
832 }
833 }
834
835
836 inline void MarginalsCorrDistribution::
lower_bounds(const IntVector & l_bnds,const BitArray & mask)837 lower_bounds(const IntVector& l_bnds, const BitArray& mask)
838 {
839 check_active_length(l_bnds, mask);
840 size_t v, num_v = randomVars.size();
841 if (mask.empty())
842 for (v=0; v<num_v; ++v)
843 randomVars[v].lower_bound(l_bnds[v]);
844 else {
845 size_t av_cntr = 0;
846 for (v=0; v<num_v; ++v)
847 if (mask[v])
848 randomVars[v].lower_bound(l_bnds[av_cntr++]);
849 }
850 }
851
852
lower_bound(Real l_bnd,size_t rv_index)853 inline void MarginalsCorrDistribution::lower_bound(Real l_bnd, size_t rv_index)
854 {
855 if (rv_index < randomVars.size())
856 randomVars[rv_index].lower_bound(l_bnd);
857 else {
858 PCerr << "Error: rv_index (" << rv_index << ") out of range in Marginals"
859 << "CorrDistribution::lower_bound(Real, size_t)" << std::endl;
860 abort_handler(-1);
861 }
862 }
863
864
lower_bound(int l_bnd,size_t rv_index)865 inline void MarginalsCorrDistribution::lower_bound(int l_bnd, size_t rv_index)
866 {
867 if (rv_index < randomVars.size())
868 randomVars[rv_index].lower_bound(l_bnd);
869 else {
870 PCerr << "Error: rv_index (" << rv_index << ") out of range in Marginals"
871 << "CorrDistribution::lower_bound(int, size_t)" << std::endl;
872 abort_handler(-1);
873 }
874 }
875
876
877 inline void MarginalsCorrDistribution::
upper_bounds(const RealVector & u_bnds,const BitArray & mask)878 upper_bounds(const RealVector& u_bnds, const BitArray& mask)
879 {
880 check_active_length(u_bnds, mask);
881 size_t v, num_v = randomVars.size();
882 if (mask.empty())
883 for (v=0; v<num_v; ++v)
884 randomVars[v].upper_bound(u_bnds[v]);
885 else {
886 size_t av_cntr = 0;
887 for (v=0; v<num_v; ++v)
888 if (mask[v])
889 randomVars[v].upper_bound(u_bnds[av_cntr++]);
890 }
891 }
892
893
894 inline void MarginalsCorrDistribution::
upper_bounds(const IntVector & u_bnds,const BitArray & mask)895 upper_bounds(const IntVector& u_bnds, const BitArray& mask)
896 {
897 check_active_length(u_bnds, mask);
898 size_t v, num_v = randomVars.size();
899 if (mask.empty())
900 for (v=0; v<num_v; ++v)
901 randomVars[v].upper_bound(u_bnds[v]);
902 else {
903 size_t av_cntr = 0;
904 for (v=0; v<num_v; ++v)
905 if (mask[v])
906 randomVars[v].upper_bound(u_bnds[av_cntr++]);
907 }
908 }
909
910
upper_bound(Real u_bnd,size_t rv_index)911 inline void MarginalsCorrDistribution::upper_bound(Real u_bnd, size_t rv_index)
912 {
913 if (rv_index < randomVars.size())
914 randomVars[rv_index].upper_bound(u_bnd);
915 else {
916 PCerr << "Error: rv_index (" << rv_index << ") out of range in Marginals"
917 << "CorrDistribution::upper_bound(Real, size_t)" << std::endl;
918 abort_handler(-1);
919 }
920 }
921
922
upper_bound(int u_bnd,size_t rv_index)923 inline void MarginalsCorrDistribution::upper_bound(int u_bnd, size_t rv_index)
924 {
925 if (rv_index < randomVars.size())
926 randomVars[rv_index].upper_bound(u_bnd);
927 else {
928 PCerr << "Error: rv_index (" << rv_index << ") out of range in Marginals"
929 << "CorrDistribution::upper_bound(int, size_t)" << std::endl;
930 abort_handler(-1);
931 }
932 }
933
934
pdf(Real val,size_t rv_index) const935 inline Real MarginalsCorrDistribution::pdf(Real val, size_t rv_index) const
936 { return randomVars[rv_index].pdf(val); }
937
938
939 inline Real MarginalsCorrDistribution::
pdf_gradient(Real val,size_t rv_index) const940 pdf_gradient(Real val, size_t rv_index) const
941 { return randomVars[rv_index].pdf_gradient(val); }
942
943
944 inline Real MarginalsCorrDistribution::
pdf_hessian(Real val,size_t rv_index) const945 pdf_hessian(Real val, size_t rv_index) const
946 { return randomVars[rv_index].pdf_hessian(val); }
947
948
log_pdf(Real val,size_t rv_index) const949 inline Real MarginalsCorrDistribution::log_pdf(Real val, size_t rv_index) const
950 { return randomVars[rv_index].log_pdf(val); }
951
952
953 inline Real MarginalsCorrDistribution::
log_pdf_gradient(Real val,size_t rv_index) const954 log_pdf_gradient(Real val, size_t rv_index) const
955 { return randomVars[rv_index].log_pdf_gradient(val); }
956
957
958 inline Real MarginalsCorrDistribution::
log_pdf_hessian(Real val,size_t rv_index) const959 log_pdf_hessian(Real val, size_t rv_index) const
960 { return randomVars[rv_index].log_pdf_hessian(val); }
961
962
pdf(const RealVector & pt) const963 inline Real MarginalsCorrDistribution::pdf(const RealVector& pt) const
964 {
965 // correlated density handled via upstream transform to standardized space
966 if (correlationFlag) {
967 PCerr << "Error: MarginalsCorrDistribution::pdf() currently uses a "
968 << "product of marginal densities\n and can only be used for "
969 << "independent random variables." << std::endl;
970 abort_handler(-1);
971 }
972
973 check_active_length(pt, activeVars);
974 size_t v, num_v = randomVars.size();
975 Real density = 1.;
976 if (activeVars.empty())
977 for (v=0; v<num_v; ++v)
978 density *= pdf(pt[v], v);
979 else {
980 size_t av_cntr = 0;
981 for (v=0; v<num_v; ++v)
982 if (activeVars[v])
983 density *= pdf(pt[av_cntr++], v);
984 }
985 return density;
986 }
987
988
log_pdf(const RealVector & pt) const989 inline Real MarginalsCorrDistribution::log_pdf(const RealVector& pt) const
990 {
991 // correlated density handled via upstream transform to standardized space
992 if (correlationFlag) {
993 PCerr << "Error: MarginalsCorrDistribution::log_pdf() currently uses a "
994 << "sum of log marginal densities\n and can only be used for "
995 << "independent random variables." << std::endl;
996 abort_handler(-1);
997 }
998
999 check_active_length(pt, activeVars);
1000 size_t v, num_v = randomVars.size();
1001 Real log_density = 0.;
1002 if (activeVars.empty())
1003 for (v=0; v<num_v; ++v)
1004 log_density += log_pdf(pt[v], v);
1005 else {
1006 size_t av_cntr = 0;
1007 for (v=0; v<num_v; ++v)
1008 if (activeVars[v])
1009 log_density += log_pdf(pt[av_cntr++], v);
1010 }
1011 return log_density;
1012 }
1013
1014
1015 template <typename Engine>
draw_sample(size_t rv_index,Engine & rng) const1016 Real MarginalsCorrDistribution::draw_sample(size_t rv_index, Engine& rng) const
1017 { return randomVars[rv_index].draw_sample(rng); }
1018
1019
1020 template <typename Engine>
1021 Real MarginalsCorrDistribution::
draw_standard_sample(size_t rv_index,Engine & rng) const1022 draw_standard_sample(size_t rv_index, Engine& rng) const
1023 { return randomVars[rv_index].draw_standard_sample(rng); }
1024
1025 } // namespace Pecos
1026
1027 #endif
1028