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, Brian Adams 12 //- Checked by: 13 //- Version: 14 15 #ifndef NOND_QUESO_BAYES_CALIBRATION_H 16 #define NOND_QUESO_BAYES_CALIBRATION_H 17 18 #include "NonDBayesCalibration.hpp" 19 20 // forward declare here to isolate QUESO includes to Dakota .cpp files 21 namespace QUESO { 22 class GslVector; 23 class GslMatrix; 24 class EnvOptionsValues; 25 class FullEnvironment; 26 template<class V, class M> class VectorSpace; 27 template<class V, class M> class BoxSubset; 28 template<class V, class M> class GenericScalarFunction; 29 template<class V, class M> class BaseVectorRV; 30 template<class V, class M> class GenericVectorRV; 31 template<class V, class M> class StatisticalInverseProblem; 32 class SipOptionsValues; 33 class MhOptionsValues; 34 } 35 36 namespace Dakota { 37 38 // forward declarations of Dakota specializations 39 template <class V, class M> class DerivInformedPropCovTK; 40 template <class V, class M> class DerivInformedPropCovLogitTK; 41 42 43 /// Bayesian inference using the QUESO library from UT Austin 44 45 /** This class wraps the Quantification of Uncertainty for Estimation, 46 Simulation, and Optimization (QUESO) library, developed as part of 47 the Predictive Science Academic Alliance Program (PSAAP)-funded 48 Predictive Engineering and Computational Sciences (PECOS) Center 49 at UT Austin. */ 50 class NonDQUESOBayesCalibration: public NonDBayesCalibration 51 { 52 /// Random walk transition kernel needs callback access to QUESO details 53 friend class DerivInformedPropCovTK<QUESO::GslVector, QUESO::GslMatrix>; 54 /// Logit random walk transition kernel needs callback access to QUESO details 55 friend class DerivInformedPropCovLogitTK<QUESO::GslVector, QUESO::GslMatrix>; 56 57 public: 58 59 // 60 //- Heading: Constructors and destructor 61 // 62 63 /// standard constructor 64 NonDQUESOBayesCalibration(ProblemDescDB& problem_db, Model& model); 65 /// destructor 66 ~NonDQUESOBayesCalibration(); 67 68 // 69 //- Heading: public member functions 70 // 71 72 protected: 73 74 // 75 //- Heading: Virtual function redefinitions 76 // 77 78 void calibrate(); 79 void print_results(std::ostream& s, short results_state = FINAL_RESULTS); 80 81 /// convenience function to print calibration parameters, e.g., for 82 /// MAP / best parameters 83 void print_variables(std::ostream& s, const RealVector& c_vars); 84 85 /// initialize the QUESO FullEnvironment on the Dakota MPIComm 86 void init_queso_environment(); 87 88 /// initialize the ASV value for preconditioned cases 89 void init_precond_request_value(); 90 91 /// define solver options, likelihood callback, posterior RV, and 92 /// inverse problem 93 virtual void init_queso_solver(); 94 95 /// use derivative information from the emulator to define the proposal 96 /// covariance (inverse of misfit Hessian) 97 void precondition_proposal(unsigned int chain_index); 98 99 /// perform the MCMC process 100 void run_queso_solver(); 101 102 void map_pre_solve(); 103 104 /// short term option to restart the MCMC chain with updated proposal 105 /// density computed from the emulator at a new starting point 106 void run_chain(); 107 108 /// cache the chain to acceptanceChain and acceptedFnVals 109 void cache_chain(); 110 111 /// log at most batchSize best chain points into bestSamples 112 void log_best(); 113 114 /// extract batchSize points from the MCMC chain and store final 115 /// aggregated set within allSamples; unique points with best 116 /// conditioning are selected, as determined by pivoted LU 117 void filter_chain_by_conditioning(); 118 119 /// copy bestSamples to allSamples to use in surrogate update 120 void best_to_all(); 121 122 /// evaluates allSamples on iteratedModel and update the mcmcModel emulator 123 /// with all{Samples,Responses} 124 void update_model(); 125 126 /// compute the L2 norm of the change in emulator coefficients 127 Real assess_emulator_convergence(); 128 129 /// intialize the QUESO parameter space, min, max, initial, and domain 130 void init_parameter_domain(); 131 132 void init_proposal_covariance(); 133 134 /// use covariance of prior distribution for setting proposal covariance 135 void prior_proposal_covariance(); 136 137 /// set proposal covariance from user-provided diagonal or matrix 138 void user_proposal_covariance(const String& input_fmt, 139 const RealVector& cov_data, 140 const String& cov_filename); 141 142 // perform sanity checks on proposalCovMatrix 143 void validate_proposal(); 144 145 /// set inverse problem options calIpOptionsValues common to all solvers 146 void set_ip_options(); 147 /// set MH-specific inverse problem options calIpMhOptionsValues 148 void set_mh_options(); 149 /// update MH-specific inverse problem options calIpMhOptionsValues 150 void update_chain_size(unsigned int size); 151 152 //The likelihood routine is in the format that QUESO requires, 153 //with a particular argument list that QUESO expects. 154 //We are not using all of these arguments but may in the future. 155 /// Log Likelihood function for call-back from QUESO to DAKOTA for evaluation 156 static double dakotaLogLikelihood(const QUESO::GslVector& paramValues, 157 const QUESO::GslVector* paramDirection, 158 const void* functionDataPtr, 159 QUESO::GslVector* gradVector, 160 QUESO::GslMatrix* hessianMatrix, 161 QUESO::GslVector* hessianEffect); 162 163 /// local copy_data utility from GslVector to RealVector 164 void copy_gsl(const QUESO::GslVector& qv, RealVector& rv); 165 /// local copy_data utility from RealVector to GslVector 166 void copy_gsl(const RealVector& rv, QUESO::GslVector& qv); 167 168 /// local copy_data utility from portion of GslVector to RealVector 169 void copy_gsl_partial(const QUESO::GslVector& qv, size_t start, 170 RealVector& rv); 171 /// local copy_data utility from RealVector to portion of GslVector 172 void copy_gsl_partial(const RealVector& rv, QUESO::GslVector& qv, 173 size_t start); 174 175 /// local copy_data utility from GslVector to column in RealMatrix 176 void copy_gsl(const QUESO::GslVector& qv, RealMatrix& rm, int i); 177 178 /// equality tester for two GslVectors 179 bool equal_gsl(const QUESO::GslVector& qv1, const QUESO::GslVector& qv2); 180 181 // 182 //- Heading: Data 183 // 184 185 /// MCMC type ("dram" or "delayed_rejection" or "adaptive_metropolis" 186 /// or "metropolis_hastings" or "multilevel", within QUESO) 187 String mcmcType; 188 /// period (number of accepted chain samples) for proposal covariance update 189 int propCovUpdatePeriod; 190 /// number of points to add to surrogate at each iteration 191 unsigned int batchSize; 192 /// the active set request value to use in proposal preconditioning 193 short precondRequestValue; 194 /// flag indicating user activation of logit transform option 195 /** this option is useful for preventing rejection or resampling for 196 out-of-bounds samples by transforming bounded domains to [-inf,inf]. */ 197 bool logitTransform; 198 199 // the following QUESO objects listed in order of construction; 200 // scoped_ptr more appropriate, but don't want to include QUESO 201 // headers here (would be needed for checked delete on scoped_ptr) 202 203 // TODO: see if this can be a local withing init function 204 /// options for setting up the QUESO Environment 205 std::shared_ptr<QUESO::EnvOptionsValues> envOptionsValues; 206 207 /// top-level QUESO Environment 208 std::shared_ptr<QUESO::FullEnvironment> quesoEnv; 209 210 /// QUESO parameter space based on number of calibrated parameters 211 std::shared_ptr<QUESO::VectorSpace<QUESO::GslVector,QUESO::GslMatrix> > 212 paramSpace; 213 214 /// QUESO parameter domain: hypercube based on min/max values 215 std::shared_ptr<QUESO::BoxSubset<QUESO::GslVector,QUESO::GslMatrix> > 216 paramDomain; 217 218 /// initial parameter values at which to start chain 219 std::shared_ptr<QUESO::GslVector> paramInitials; 220 221 /// random variable for the prior 222 std::shared_ptr<QUESO::BaseVectorRV<QUESO::GslVector,QUESO::GslMatrix> > 223 priorRv; 224 225 /// proposal covariance for DRAM 226 std::shared_ptr<QUESO::GslMatrix> proposalCovMatrix; 227 228 /// optional multiplier to scale prior-based proposal covariance 229 double priorPropCovMult; 230 231 /// general inverse problem options 232 std::shared_ptr<QUESO::SipOptionsValues> calIpOptionsValues; 233 234 /// MH-specific inverse problem options 235 std::shared_ptr<QUESO::MhOptionsValues> calIpMhOptionsValues; 236 237 std::shared_ptr<QUESO::GenericScalarFunction<QUESO::GslVector, 238 QUESO::GslMatrix> > likelihoodFunctionObj; 239 240 /// random variable for the posterior 241 std::shared_ptr<QUESO::GenericVectorRV<QUESO::GslVector,QUESO::GslMatrix> > 242 postRv; 243 244 /// QUESO inverse problem solver 245 std::shared_ptr<QUESO::StatisticalInverseProblem<QUESO::GslVector, 246 QUESO::GslMatrix> > inverseProb; 247 248 /// advanced options file name (GPMSA only); settings from this file 249 /// override any C++ / Dakota input file settings 250 String advancedOptionsFile; 251 252 /// Pointer to current class instance for use in static callback functions 253 static NonDQUESOBayesCalibration* nonDQUESOInstance; 254 255 private: 256 257 // 258 // - Heading: Data 259 // 260 261 /// cache previous expansion coefficients for assessing convergence of 262 /// emulator refinement process 263 RealVectorArray prevCoeffs; 264 }; 265 266 } // namespace Dakota 267 268 #endif 269