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: NonDBayesCalibration
10 //- Description: Base class for generic Bayesian inference
11 //- Owner: Laura Swiler
12 //- Checked by:
13 //- Version:
14
15 #ifndef NOND_BAYES_CALIBRATION_H
16 #define NOND_BAYES_CALIBRATION_H
17
18 #include "NonDCalibration.hpp"
19 #include "MarginalsCorrDistribution.hpp"
20 #include "InvGammaRandomVariable.hpp"
21 #include "GaussianKDE.hpp"
22 #include "ANN/ANN.h"
23
24 //#define DEBUG
25
26 namespace Dakota {
27
28
29 /// Base class for Bayesian inference: generates posterior
30 /// distribution on model parameters given experimental data
31
32 /** This class will eventually provide a general-purpose framework for
33 Bayesian inference. In the short term, it only collects shared
34 code between QUESO and GPMSA implementations. */
35
36 class NonDBayesCalibration: public NonDCalibration
37 {
38 public:
39
40 //
41 //- Heading: Constructors and destructor
42 //
43
44 /// standard constructor
45 NonDBayesCalibration(ProblemDescDB& problem_db, Model& model);
46 /// destructor
47 ~NonDBayesCalibration();
48
49 //
50 //- Heading: Member functions
51 //
52
53 /// compute the prior PDF for a particular MCMC sample
54 template <typename VectorType>
55 Real prior_density(const VectorType& vec);
56 /// compute the log prior PDF for a particular MCMC sample
57 template <typename VectorType>
58 Real log_prior_density(const VectorType& vec);
59
60 /// draw a multivariate sample from the prior distribution
61 template <typename Engine>
62 void prior_sample(Engine& rng, RealVector& prior_samples);
63
64 /// return the mean of the prior distribution
65 template <typename VectorType>
66 void prior_mean(VectorType& mean_vec) const;
67
68 /// return the covariance of the prior distribution
69 template <typename MatrixType>
70 void prior_variance(MatrixType& var_mat) const;
71
72 /**
73 * \brief Compute the proposal covariance C based on low-rank approximation
74 * to the prior-preconditioned misfit Hessian.
75 */
76 static void get_positive_definite_covariance_from_hessian(
77 const RealSymMatrix &hessian, const RealMatrix& prior_chol_fact,
78 RealSymMatrix &covariance, short output_lev);
79
80 // compute information metrics
81 static Real knn_kl_div(RealMatrix& distX_samples, RealMatrix& distY_samples,
82 size_t dim);
83 static Real knn_mutual_info(RealMatrix& Xmatrix, int dimX, int dimY,
84 unsigned short alg);
85
86 protected:
87
88 //
89 //- Heading: Virtual function redefinitions
90 //
91
92 void pre_run();
93 void core_run();
94
95 void derived_init_communicators(ParLevLIter pl_iter);
96 void derived_set_communicators(ParLevLIter pl_iter);
97 void derived_free_communicators(ParLevLIter pl_iter);
98
99 void print_results(std::ostream& s, short results_state = FINAL_RESULTS);
100
101 const Model& algorithm_space_model() const;
102
103 //
104 //- Heading: New virtual functions
105 //
106
107 /// Perform Bayesian calibration (all derived classes must implement)
108 virtual void calibrate() = 0;
109
110 //
111 //- Heading: Member functions
112 //
113
114 /// construct mcmcModel (no emulation, GP, PCE, or SC) that wraps
115 /// inbound Model
116 void construct_mcmc_model();
117
118 /// initialize the hyper-parameter priors
119 void init_hyper_parameters();
120
121 /// initialize the MAP optimizer selection
122 void init_map_optimizer();
123 /// construct the negative log posterior RecastModel (negLogPostModel)
124 void construct_map_model();
125 /// construct the MAP optimizer for minimization of negLogPostModel
126 void construct_map_optimizer();
127
128 /// initialize emulator model and probability space transformations
129 void initialize_model();
130
131 /// calibrate the model to a high-fidelity data source, using mutual
132 /// information-guided design of experiments (adaptive experimental
133 /// design)
134 void calibrate_to_hifi();
135 /// evaluate stopping criteria for calibrate_to_hifi
136 bool eval_hi2lo_stop(bool stop_metric, double prev_MI, double max_MI,
137 int num_it, int num_hifi, int max_hifi,
138 int num_candidates);
139 /// print calibrate_to_hifi progress to file
140 void print_hi2lo_file(std::ostream& out_file, int num_it, int batchEvals,
141 RealMatrix& optimal_config_matrix, const RealVector&
142 MI_vec, int max_hifi, RealMatrix& resp_matrix, const
143 RealVector& optimal_config, double max_MI);
144 /// print calibrate_to_hifi progress
145 void print_hi2lo_begin(int num_it);
146 void print_hi2lo_status(int num_it, int i, const RealVector& xi_i, double MI);
147 void print_hi2lo_batch_status(int num_it, int batch_n, int batchEvals,
148 const RealVector& optimal_config, double max_MI);
149 void print_hi2lo_selected(int num_it, int batchEvals, RealMatrix&
150 optimal_config_matrix, const RealVector&
151 optimal_config, double max_MI);
152 /// supplement high-fidelity data with LHS samples
153 void add_lhs_hifi_data();
154 /// apply simulation error vector
155 void apply_error_vec(const RealVector& error_vec, int &seed, int experiment);
156 /// build matrix of errors
157 void build_error_matrix(const RealVector& sim_error_vec,
158 RealMatrix& sim_error_matrix, int &seed);
159 /// build matrix of candidate points
160 void build_designs(RealMatrix& design_matrix);
161 /// build matrix to calculate mutual information for calibrate_to_hifi
162 void build_hi2lo_xmatrix(RealMatrix& Xmatrix, int i, const RealMatrix&
163 mi_chain, RealMatrix& sim_error_matrix);
164 /// run high-fidelity model at several configs and add to experiment data
165 void run_hifi(RealMatrix& optimal_config_matrix, RealMatrix& resp_matrix);
166
167 /// calculate model discrepancy with respect to experimental data
168 void build_model_discrepancy();
169 void build_scalar_discrepancy();
170 void build_field_discrepancy();
171 void build_GP_field(const RealMatrix& t, RealMatrix& t_pred,
172 const RealVector& concat_disc, RealVector& disc_pred,
173 RealVector& disc_var);
174
175
176 /// calculate a Kernel Density Estimate (KDE) for the posterior samples
177 void calculate_kde();
178 /// calculate the model evidence
179 void calculate_evidence();
180
181 void extract_selected_posterior_samples(const std::vector<int> &points_to_keep,
182 const RealMatrix &samples_for_posterior_eval,
183 const RealVector &posterior_density, RealMatrix &posterior_data ) const;
184
185 void export_posterior_samples_to_file(const std::string filename,
186 const RealMatrix &posterior_data ) const;
187
188 /// compute the (approximate) gradient of the negative log posterior by
189 /// augmenting the (approximate) gradient of the negative log likelihood
190 /// with the gradient of the negative log prior
191 template <typename VectorType1, typename VectorType2>
192 void augment_gradient_with_log_prior(VectorType1& log_grad,
193 const VectorType2& vec);
194 /// compute the (approximate) Hessian of the negative log posterior by
195 /// augmenting the (approximate) Hessian of the negative log likelihood
196 /// with the Hessian of the negative log prior
197 template <typename MatrixType, typename VectorType>
198 void augment_hessian_with_log_prior(MatrixType& log_hess,
199 const VectorType& vec);
200
201 /// calculate log-likelihood from the passed residuals (assuming
202 /// they are already sized and scaled by covariance / hyperparams...
203 Real
204 log_likelihood(const RealVector& residuals, const RealVector& hyper_params);
205
206 /// compute priorCovCholFactor based on prior distributions for random
207 /// variables and any hyperparameters
208 void prior_cholesky_factorization();
209
210 /// member version forwards member data to static function
211 void get_positive_definite_covariance_from_hessian(
212 const RealSymMatrix &hessian, RealSymMatrix &covariance);
213
214 /// static function passed by pointer to negLogPostModel recast model
215 static void neg_log_post_resp_mapping(const Variables& model_vars,
216 const Variables& nlpost_vars,
217 const Response& model_resp,
218 Response& nlpost_resp);
219
220 /// Wrap iteratedModel in a RecastModel that performs response scaling
221 void scale_model();
222
223 /// Wrap iteratedModel in a RecastModel that weights the residuals
224 void weight_model();
225
226 //
227 //- Heading: Data
228 //
229
230 String scalarDataFilename;
231
232 // technically doesn't apply to GPMSA, but leaving here for now
233 /// the emulator type: NO_EMULATOR, GP_EMULATOR, PCE_EMULATOR,
234 /// SC_EMULATOR, ML_PCE_EMULATOR, MF_PCE_EMULATOR, or MF_SC_EMULATOR
235 short emulatorType;
236 /// Model instance employed in the likelihood function; provides response
237 /// function values from Gaussian processes, stochastic expansions (PCE/SC),
238 /// or direct access to simulations (no surrogate option)
239 Model mcmcModel;
240
241 /// whether the MCMC Model is a surrogate, or a thin transformation
242 /// around a surrogate, so can be cheaply re-evaluated in chain recovery
243 bool mcmcModelHasSurrogate;
244
245 /// DataTransformModel wrapping the mcmcModel
246 Model residualModel;
247
248 /// SQP or NIP optimizer for pre-solving for the MAP point prior to MCMC.
249 /// This is restricted to emulator cases for now, but as for derivative
250 /// preconditioning, could be activated for no-emulator cases with a
251 /// specification option (not active by default).
252 Iterator mapOptimizer;
253 /// RecastModel for solving for MAP: reduces residualModel to scalar
254 /// negative log posterior
255 Model negLogPostModel;
256 /// user setting of the MAP optimization algorithm type
257 unsigned short mapOptAlgOverride;
258
259 /// NonDPolynomialChaos or NonDStochCollocation instance for defining a
260 /// PCE/SC-based mcmcModel
261 Iterator stochExpIterator;
262
263 /// number of samples in the chain (e.g. number of MCMC samples);
264 /// for iterative update cycles, number of samples per update cycle
265 int chainSamples;
266 /// random seed for MCMC process
267 int randomSeed;
268
269 /// order of derivatives used in MCMC process (bitwise like ASV)
270 short mcmcDerivOrder;
271
272 /// solution for most recent MAP pre-solve; also serves as initial guess
273 /// for initializing the first solve and warm-starting the next solve
274 /// (posterior emulator refinement)
275 RealVector mapSoln;
276
277 // settings specific to adaptive DOE
278
279 /// whether to perform iterative design of experiments with
280 /// high-fidelity model
281 bool adaptExpDesign;
282 /// number of candidate designs for adaptive Bayesian experimental design
283 size_t numCandidates;
284 /// whether to import candidate design points for adaptive Bayesian
285 /// experimental design
286 String importCandPtsFile;
287 /// tabular format for the candidate design points import file
288 unsigned short importCandFormat;
289 /// maximum number of high-fidelity model runs to be used for adaptive
290 /// Bayesian experimental design
291 int maxHifiEvals;
292 /// number of optimal designs selected per iteration of experimental design
293 /// algorithm
294 int batchEvals;
295 /// algorithm to employ in calculating mutual information
296 unsigned short mutualInfoAlg;
297
298 // settings specific to model discrepancy
299
300 /// need field coordinates for model discrepancy
301 bool readFieldCoords;
302 /// flag whether to calculate model discrepancy
303 bool calModelDiscrepancy;
304 /// set discrepancy type
305 String discrepancyType;
306 /// filename for corrected model (model+discrepancy) calculations output
307 String exportCorrModelFile;
308 /// filename for discrepancy calculations output
309 String exportDiscrepFile;
310 /// filename for corrected model variance calculations
311 String exportCorrVarFile;
312 /// format options for corrected model output
313 unsigned short exportCorrModelFormat;
314 /// format options for discrepancy output
315 unsigned short exportDiscrepFormat;
316 /// format options for corrected model variance output
317 unsigned short exportCorrVarFormat;
318 /// specify polynomial or trend order for discrepancy correction
319 short approxCorrectionOrder;
320 /// number of prediction configurations at which to calculate model
321 /// discrepancy
322 size_t numPredConfigs;
323 /// lower bounds on configuration domain
324 RealVector configLowerBnds;
325 /// upper bounds on configuration domain
326 RealVector configUpperBnds;
327 /// array containing predicted of model+discrepancy
328 ResponseArray discrepancyResponses;
329 /// array containing predicted of model+discrepancy
330 ResponseArray correctedResponses;
331 /// matrix containing variances of model+discrepancy
332 RealMatrix correctedVariances;
333 /// list of prediction configurations at which to calculate model discrepancy
334 RealVector predictionConfigList;
335 /// whether to import prediction configurations at which to calculate model
336 /// discrepancy
337 String importPredConfigs;
338 /// tabular format for prediction configurations import file
339 unsigned short importPredConfigFormat;
340 /// print tabular files containing model+discrepancy responses and variances
341 void export_discrepancy(RealMatrix& pred_config_mat);
342 /// print tabular files containing model+discrepancy responses and variances
343 /// for field responses
344 void export_field_discrepancy(RealMatrix& pred_vars_mat);
345 /// array containing predicted of model+discrepancy
346 RealVector discrepancyFieldResponses;
347 /// array containing predicted of model+discrepancy
348 RealVector correctedFieldResponses;
349 /// matrix containing variances of model+discrepancy
350 RealVector correctedFieldVariances;
351
352 /// a high-fidelity model data source (given by pointer in input)
353 Model hifiModel;
354 /// initial high-fidelity model samples
355 int initHifiSamples;
356 /// LHS iterator to generate hi-fi model data
357 Iterator hifiSampler;
358
359 /// the Cholesky factor of the prior covariance
360 RealMatrix priorCovCholFactor;
361
362 /// mode for number of observation error multipliers to calibrate
363 /// (default none)
364 unsigned short obsErrorMultiplierMode;
365 /// calculated number of hyperparameters augmenting the calibration
366 /// parameter set, e.g., due to calibrate observation error multipliers
367 int numHyperparams;
368 /// alphas for inverse gamma distribution on hyper-params
369 RealVector invGammaAlphas;
370 /// alphas for inverse gamma distribution on hyper-params
371 RealVector invGammaBetas;
372 /// distributions for hyper-params
373 std::vector<Pecos::RandomVariable> invGammaDists;
374
375 /// flag indicating use of a variable transformation to standardized
376 /// probability space for the model or emulator
377 bool standardizedSpace;
378 /// flag indicating the calculation of KL divergence between prior and
379 /// posterior
380 bool posteriorStatsKL;
381 /// flag indicating the calculation of mutual information between prior and
382 /// posterior
383 bool posteriorStatsMutual;
384 /// flag indicating the calculation of the kernel density estimate of the
385 /// posteriors
386 bool posteriorStatsKDE;
387 /// flag indicating calculation of chain diagnostics
388 bool chainDiagnostics;
389 /// flag indicating calculation of confidence intervals as a chain
390 /// diagnositc
391 bool chainDiagnosticsCI;
392 /// flag indicating calculation of the evidence of the model
393 bool calModelEvidence;
394 /// flag indicating use of Monte Carlo approximation to calculate evidence
395 bool calModelEvidMC;
396 /// flag indicating use of Laplace approximation to calculate evidence
397 bool calModelEvidLaplace;
398 /// number of samples to be used in model evidence calculation
399 int evidenceSamples;
400 /// flag indicating usage of adaptive posterior refinement; currently makes
401 /// sense for unstructured grids in GP and PCE least squares/CS
402 bool adaptPosteriorRefine;
403
404 /// approach for defining proposal covariance
405 String proposalCovarType;
406 /// data from user input of proposal covariance
407 RealVector proposalCovarData;
408 /// filename for user-specified proposal covariance
409 String proposalCovarFilename;
410 /// approach for defining proposal covariance
411 String proposalCovarInputType;
412
413 /// Post-processing-related controls
414
415 /// accumulation of acceptance chain across restarts (stored in
416 /// user-space) TO DO: retire once restarts are retired; optimize to
417 /// convert to user-space only in final results
418 RealMatrix acceptanceChain;
419 /// cached function values corresponding to acceptanceChain for
420 /// final statistics reporting
421 RealMatrix acceptedFnVals;
422
423 /// container for managing best MCMC samples (points and associated
424 /// log posterior) collected across multiple (restarted) chains
425 std::/*multi*/map<Real, RealVector> bestSamples;
426
427 /// number of MCMC samples to discard from acceptance chain
428 int burnInSamples;
429 /// period or skip in post-processing the acceptance chain
430 int subSamplingPeriod;
431
432 /// Pointer to current class instance for use in static callback functions
433 static NonDBayesCalibration* nonDBayesInstance;
434
435 /// Compute final stats for MCMC chains
436 void compute_statistics();
437 RealMatrix chainStats;
438 RealMatrix fnStats;
439
440 /// export the acceptance chain in user space
441 void export_chain();
442
443 /// Print filtered posterior and function values (later: credibility
444 /// and prediction intervals)
445 void export_chain(RealMatrix& filtered_chain, RealMatrix& filtered_fn_vals);
446
447 /// Perform chain filtering based on target chain length
448 void filter_chain(const RealMatrix& acceptance_chain, RealMatrix& filtered_chain,
449 int target_length);
450 /// Perform chain filtering with burn-in and sub-sampling
451 void filter_chain(const RealMatrix& acceptance_chain, RealMatrix& filtered_chain);
452 void filter_fnvals(const RealMatrix& accepted_fn_vals, RealMatrix& filtered_fn_vals);
453
454 /// return a newly allocated filtered matrix including start_index and
455 /// every stride-th index after; for burn-in cases, start_index is the
456 /// number of burn-in discards
457 void filter_matrix_cols(const RealMatrix& orig_matrix, int start_index,
458 int stride, RealMatrix& filtered_matrix);
459
460 /// Compute credibility and prediction intervals of final chain
461 RealMatrix predVals;
462 /// cached filtered function values for printing (may be a view of acceptedFnVals)
463 RealMatrix filteredFnVals;
464 void compute_intervals();
465 void compute_prediction_vals(RealMatrix& filtered_fn_vals,
466 RealMatrix& PredVals, int num_filtered,
467 size_t num_exp, size_t num_concatenated);
468 void print_intervals_file(std::ostream& stream, RealMatrix& functionvalsT,
469 RealMatrix& predvalsT, int length,
470 size_t aug_length);
471 void print_intervals_screen(std::ostream& stream, RealMatrix& functionvalsT,
472 RealMatrix& predvalsT, int length);
473 /// output filename for the MCMC chain
474 String exportMCMCFilename;
475 // BMA TODO: user control of filtered file name and format? Or use
476 // a base name + extensions?
477 /// output formatting options for MCMC export
478 short exportMCMCFormat;
479 short filteredMCMCFormat;
480
481 /// Compute information metrics
482 void kl_post_prior(RealMatrix& acceptanceChain);
483 void prior_sample_matrix(RealMatrix& prior_dist_samples);
484 void mutual_info_buildX();
485 static void ann_dist(const ANNpointArray matrix1, const ANNpointArray matrix2,
486 RealVector& distances, int NX, int NY, int dim2, IntVector& k,
487 double eps);
488 static void ann_dist(const ANNpointArray matrix1,
489 const ANNpointArray matrix2, RealVector& distances,
490 Int2DArray& indices, int NX, int NY, int dim2,
491 IntVector& k, double eps);
492 Real kl_est;
493 void print_kl(std::ostream& stream);
494 void print_chain_diagnostics(std::ostream& s);
495 void print_batch_means_intervals(std::ostream& s);
496
497 /// whether response scaling is active
498 bool scaleFlag;
499 // /// Shallow copy of the scaling transformation model, when present
500 // /// (cached in case further wrapped by other transformations)
501 // Model scalingModel;
502
503 /// whether weight scaling is active
504 bool weightFlag;
505
506 private:
507
508 //
509 //- Heading: Data
510 //
511
512 };
513
514
algorithm_space_model() const515 inline const Model& NonDBayesCalibration::algorithm_space_model() const
516 { return mcmcModel; }
517
518
519 template <typename VectorType>
prior_density(const VectorType & vec)520 Real NonDBayesCalibration::prior_density(const VectorType& vec)
521 {
522 // template supports QUESO::GslVector which has to use long form
523
524 // TO DO: consider QUESO-based approach for this using priorRv.pdf(),
525 // which may in turn call back to our GenericVectorRV prior plug-in
526
527 const Pecos::MultivariateDistribution& mv_dist
528 = (standardizedSpace) ? mcmcModel.multivariate_distribution()
529 : iteratedModel.multivariate_distribution();
530 if (mv_dist.correlation()) {
531 Cerr << "Error: prior_density() uses a product of marginal densities\n"
532 << " and can only be used for independent random variables."
533 << std::endl;
534 abort_handler(METHOD_ERROR);
535 }
536
537 // Don't need to map indices so long as activeVars corresponds to vec
538 //const SharedVariablesData& svd
539 // = iteratedModel.current_variables().shared_data();
540
541 const BitArray& active_vars = mv_dist.active_variables();
542 bool no_mask = active_vars.empty();
543 size_t v, num_rv = mv_dist.random_variables().size(),
544 active_rv = (no_mask) ? num_rv : active_vars.count();
545 if (active_rv != numContinuousVars) { // avoid call to .length()/.size()
546 PCerr << "Error: active variable size mismatch in NonDBayesCalibration::"
547 << "prior_density(): " << active_rv << " expected, "
548 << numContinuousVars << " provided." << std::endl;
549 abort_handler(METHOD_ERROR);
550 }
551 Real pdf = 1.;
552 if (no_mask)
553 for (v=0; v<num_rv; ++v) {
554 #ifdef DEBUG
555 Cout << "Variable " << v+1 << " = " << vec[v] << " prior density = "
556 << mv_dist.pdf(vec[v], v) << std::endl;
557 #endif // DEBUG
558 pdf *= mv_dist.pdf(vec[v], v);
559 }
560 else {
561 size_t av_cntr = 0;
562 for (v=0; v<num_rv; ++v)
563 if (active_vars[v]) {
564 #ifdef DEBUG
565 Cout << "Variable " << v+1 << " = " << vec[av_cntr]<<" prior density = "
566 << mv_dist.pdf(vec[av_cntr], v) << std::endl;
567 #endif // DEBUG
568 pdf *= mv_dist.pdf(vec[av_cntr++], v);
569 }
570 }
571
572 // the estimated param is mult^2 ~ invgamma(alpha,beta)
573 for (size_t i=0; i<numHyperparams; ++i)
574 pdf *= invGammaDists[i].pdf(vec[numContinuousVars + i]);
575
576 return pdf;
577 }
578
579
580 template <>
prior_density(const RealVector & vec)581 inline Real NonDBayesCalibration::prior_density(const RealVector& vec)
582 {
583 // template specialization supports RealVector (e.g., for MAP pre-solve)
584
585 const Pecos::MultivariateDistribution& mv_dist
586 = (standardizedSpace) ? mcmcModel.multivariate_distribution() // u_dist
587 : iteratedModel.multivariate_distribution(); // x_dist
588
589 Real pdf;
590 if (numHyperparams) {
591 // evaluate PDF for random variables
592 RealVector rv_vec(Teuchos::View, vec.values(), numContinuousVars);
593 pdf = mv_dist.pdf(rv_vec);
594 // augment w/ hyperparams: estimated param is mult^2 ~ invgamma(alpha,beta)
595 for (size_t i=0; i<numHyperparams; ++i)
596 pdf *= invGammaDists[i].pdf(vec[numContinuousVars + i]);
597 }
598 else
599 pdf = mv_dist.pdf(vec);
600
601 return pdf;
602 }
603
604
605 template <typename VectorType>
log_prior_density(const VectorType & vec)606 Real NonDBayesCalibration::log_prior_density(const VectorType& vec)
607 {
608 const Pecos::MultivariateDistribution& mv_dist
609 = (standardizedSpace) ? mcmcModel.multivariate_distribution()
610 : iteratedModel.multivariate_distribution();
611 if (mv_dist.correlation()) {
612 Cerr << "Error: log_prior_density() uses a sum of log marginal densities\n"
613 << " and can only be used for independent random variables."
614 << std::endl;
615 abort_handler(METHOD_ERROR);
616 }
617
618 // Don't need to map indices so long as activeVars corresponds to vec
619 //const SharedVariablesData& svd
620 // = iteratedModel.current_variables().shared_data();
621
622 const BitArray& active_vars = mv_dist.active_variables();
623 bool no_mask = active_vars.empty();
624 size_t v, num_rv = mv_dist.random_variables().size(),
625 active_rv = (no_mask) ? num_rv : active_vars.count();
626 if (active_rv != numContinuousVars) { // avoid call to .length()/.size()
627 PCerr << "Error: active variable size mismatch in NonDBayesCalibration::"
628 << "log_prior_density(): " << active_rv << " expected, "
629 << numContinuousVars << " provided." << std::endl;
630 abort_handler(METHOD_ERROR);
631 }
632 Real log_pdf = 0.;
633 if (no_mask)
634 for (v=0; v<num_rv; ++v) {
635 #ifdef DEBUG
636 Cout << "Variable " << v+1 << " = " << vec[v] << " prior density = "
637 << mv_dist.pdf(vec[v], v) << " log prior density = "
638 << mv_dist.log_pdf(vec[v], v) << std::endl;
639 #endif // DEBUG
640 log_pdf += mv_dist.log_pdf(vec[v], v);
641 }
642 else {
643 size_t av_cntr = 0;
644 for (v=0; v<num_rv; ++v)
645 if (active_vars[v]) {
646 #ifdef DEBUG
647 Cout << "Variable " << v+1 << " = " << vec[av_cntr]
648 << " prior density = " << mv_dist.pdf(vec[av_cntr], v)
649 << " log prior density = " << mv_dist.log_pdf(vec[av_cntr], v)
650 << std::endl;
651 #endif // DEBUG
652 log_pdf += mv_dist.log_pdf(vec[av_cntr++], v);
653 }
654 }
655
656 // the estimated param is mult^2 ~ invgamma(alpha,beta)
657 for (size_t i=0; i<numHyperparams; ++i)
658 log_pdf += invGammaDists[i].log_pdf(vec[numContinuousVars + i]);
659
660 return log_pdf;
661 }
662
663
664 template <>
log_prior_density(const RealVector & vec)665 inline Real NonDBayesCalibration::log_prior_density(const RealVector& vec)
666 {
667 const Pecos::MultivariateDistribution& mv_dist
668 = (standardizedSpace) ? mcmcModel.multivariate_distribution() // u_dist
669 : iteratedModel.multivariate_distribution(); // x_dist
670
671 Real log_pdf;
672 if (numHyperparams) {
673 // evaluate PDF for random variables
674 RealVector rv_vec(Teuchos::View, vec.values(), numContinuousVars);
675 log_pdf = mv_dist.log_pdf(rv_vec);
676 // augment w/ hyperparams: estimated param is mult^2 ~ invgamma(alpha,beta)
677 for (size_t i=0; i<numHyperparams; ++i)
678 log_pdf += invGammaDists[i].log_pdf(vec[numContinuousVars + i]);
679 }
680 else
681 log_pdf = mv_dist.log_pdf(vec);
682
683 return log_pdf;
684 }
685
686
687 template <typename Engine>
prior_sample(Engine & rng,RealVector & prior_samples)688 void NonDBayesCalibration::prior_sample(Engine& rng, RealVector& prior_samples)
689 {
690 if (prior_samples.empty())
691 prior_samples.sizeUninitialized(numContinuousVars + numHyperparams);
692
693 const Pecos::MultivariateDistribution& mv_dist
694 = (standardizedSpace) ? mcmcModel.multivariate_distribution()
695 : iteratedModel.multivariate_distribution();
696 const std::shared_ptr<Pecos::MarginalsCorrDistribution> mv_dist_rep =
697 std::static_pointer_cast<Pecos::MarginalsCorrDistribution>
698 (mv_dist.multivar_dist_rep());
699 const SharedVariablesData& svd
700 = iteratedModel.current_variables().shared_data();
701 if (mv_dist_rep->correlation()) {
702 Cerr << "Error: prior_sample() does not support correlated prior samples."
703 << std::endl;
704 abort_handler(METHOD_ERROR);
705 }
706 for (size_t i=0; i<numContinuousVars; ++i)
707 prior_samples[i]
708 = mv_dist_rep->draw_sample(svd.cv_index_to_all_index(i), rng);
709
710 // the estimated param is mult^2 ~ invgamma(alpha,beta)
711 for (size_t i=0; i<numHyperparams; ++i)
712 prior_samples[numContinuousVars + i] = invGammaDists[i].draw_sample(rng);
713 }
714
715
716 /** Assume the target mean_vec is sized by client */
717 template <typename VectorType>
prior_mean(VectorType & mean_vec) const718 void NonDBayesCalibration::prior_mean(VectorType& mean_vec) const
719 {
720 // SVD index conversion is more general, but not required for current uses
721 //const SharedVariablesData& svd
722 // = iteratedModel.current_variables().shared_data();
723
724 // returns set of means that correspond to activeVars
725 RealVector mvd_means
726 = (standardizedSpace) ? mcmcModel.multivariate_distribution().means()
727 : iteratedModel.multivariate_distribution().means();
728
729 for (size_t i=0; i<numContinuousVars; ++i)
730 mean_vec[i] = mvd_means[i];//[svd.cv_index_to_active_index(i)];
731
732 for (size_t i=0; i<numHyperparams; ++i)
733 mean_vec[numContinuousVars + i] = invGammaDists[i].mean();
734 }
735
736
737 /** Assumes the target var_mat is sized by client */
738 template <typename MatrixType>
prior_variance(MatrixType & var_mat) const739 void NonDBayesCalibration::prior_variance(MatrixType& var_mat) const
740 {
741 // SVD index conversion is more general, but not required for current uses
742 //const SharedVariablesData& svd
743 // = iteratedModel.current_variables().shared_data();
744
745 if (standardizedSpace) {
746 // returns set of RV variances that correspond to activeVars
747 RealVector u_var = mcmcModel.multivariate_distribution().variances();
748 for (size_t i=0; i<numContinuousVars; ++i)
749 var_mat(i,i) = u_var[i];//[svd.cv_index_to_active_index(i)];
750 }
751 else {
752 const Pecos::MultivariateDistribution& x_dist
753 = iteratedModel.multivariate_distribution();
754 if (x_dist.correlation()) {
755 // returns set of RV stdevs that correspond to activeVars
756 RealVector x_std = x_dist.std_deviations();
757 const RealSymMatrix& x_correl = x_dist.correlation_matrix();
758 size_t i, j;//, rv_i, rv_j;
759 for (i=0; i<numContinuousVars; ++i) {
760 //rv_i = svd.cv_index_to_active_index(i);
761 var_mat(i,i) = x_std[i] * x_std[i];//x_std[rv_i] * x_std[rv_i];
762 for (j=0; j<i; ++j) {
763 //rv_j = svd.cv_index_to_active_index(j);
764 var_mat(i,j) = var_mat(j,i) = x_correl(i,j) * x_std[i] * x_std[j];
765 //var_mat(i,j) = var_mat(j,i)
766 // = x_correl(rv_i,rv_j) * x_std[rv_i] * x_std[rv_j];
767 }
768 }
769 }
770 else {
771 // returns set of RV variances that correspond to activeVars
772 RealVector x_var = x_dist.variances();
773 for (size_t i=0; i<numContinuousVars; ++i)
774 var_mat(i,i) = x_var[i];//[svd.cv_index_to_active_index(i)];
775 }
776 }
777
778 for (size_t i=0; i<numHyperparams; ++i)
779 var_mat(numContinuousVars + i, numContinuousVars + i) =
780 invGammaDists[i].variance();
781 }
782
783
784 template <typename VectorType1, typename VectorType2>
785 void NonDBayesCalibration::
augment_gradient_with_log_prior(VectorType1 & log_grad,const VectorType2 & vec)786 augment_gradient_with_log_prior(VectorType1& log_grad, const VectorType2& vec)
787 {
788 // neg log posterior = neg log likelihood + neg log prior = misfit - log prior
789 // --> gradient of neg log posterior = misfit gradient - log prior gradient
790 const Pecos::MultivariateDistribution& mv_dist
791 = (standardizedSpace) ? mcmcModel.multivariate_distribution()
792 : iteratedModel.multivariate_distribution();
793 const SharedVariablesData& svd
794 = iteratedModel.current_variables().shared_data();
795
796 for (size_t i=0; i<numContinuousVars; ++i)
797 log_grad[i] -=
798 mv_dist.log_pdf_gradient(vec[i], svd.cv_index_to_all_index(i));
799 }
800
801
802 template <typename MatrixType, typename VectorType>
803 void NonDBayesCalibration::
augment_hessian_with_log_prior(MatrixType & log_hess,const VectorType & vec)804 augment_hessian_with_log_prior(MatrixType& log_hess, const VectorType& vec)
805 {
806 // neg log posterior = neg log likelihood + neg log prior = misfit - log prior
807 // --> Hessian of neg log posterior = misfit Hessian - log prior Hessian
808 const Pecos::MultivariateDistribution& mv_dist
809 = (standardizedSpace) ? mcmcModel.multivariate_distribution()
810 : iteratedModel.multivariate_distribution();
811 const SharedVariablesData& svd
812 = iteratedModel.current_variables().shared_data();
813
814 for (size_t i=0; i<numContinuousVars; ++i)
815 log_hess(i, i) -=
816 mv_dist.log_pdf_hessian(vec[i], svd.cv_index_to_all_index(i));
817 }
818
819
820 /*
821 template <>
822 inline void NonDBayesCalibration::
823 augment_gradient_with_log_prior(RealVector& log_grad, const RealVector& vec)
824 {
825 // neg log posterior = neg log likelihood + neg log prior = misfit - log prior
826 // --> Hessian of neg log posterior = misfit Hessian - log prior Hessian
827 const Pecos::MultivariateDistribution& mv_dist
828 = (standardizedSpace) ? mcmcModel.multivariate_distribution()
829 : iteratedModel.multivariate_distribution();
830
831 log_grad -= mv_dist.log_pdf_gradient(vec);
832 }
833
834
835 template <>
836 inline void NonDBayesCalibration::
837 augment_hessian_with_log_prior(RealSymMatrix& log_hess, const RealVector& vec)
838 {
839 // neg log posterior = neg log likelihood + neg log prior = misfit - log prior
840 // --> Hessian of neg log posterior = misfit Hessian - log prior Hessian
841 const Pecos::MultivariateDistribution& mv_dist
842 = (standardizedSpace) ? mcmcModel.multivariate_distribution()
843 : iteratedModel.multivariate_distribution();
844
845 log_hess -= mv_dist.log_pdf_hessian(vec); // only need diagonal correction...
846 }
847 */
848
849
850 inline void NonDBayesCalibration::
get_positive_definite_covariance_from_hessian(const RealSymMatrix & hessian,RealSymMatrix & covariance)851 get_positive_definite_covariance_from_hessian( const RealSymMatrix &hessian,
852 RealSymMatrix &covariance )
853 {
854 get_positive_definite_covariance_from_hessian(hessian, priorCovCholFactor,
855 covariance, outputLevel);
856 }
857
858 } // namespace Dakota
859
860 #endif
861