BIE
/home/weinberg/src/BIE/include/Simulation.h
00001 // This is really -*- C++ -*- 00002 00003 00004 #ifndef Simulation_h 00005 #define Simulation_h 00006 00007 #include "Serializable.h" 00008 00009 00010 #include <iostream> 00011 #include <iomanip> 00012 #include <fstream> 00013 #include <cmath> 00014 #include <vector> 00015 00016 using namespace std; 00017 00018 #include "ACG.h" 00019 #include "Uniform.h" 00020 #include "Normal.h" 00021 00022 #include "BIEconfig.h" 00023 #include "BaseDataTree.h" 00024 #include "Frontier.h" 00025 #include "Distribution.h" 00026 #include "MHWidth.h" 00027 #include "MixturePrior.h" 00028 #include "Model.h" 00029 #include "Converge.h" 00030 #include "Integration.h" 00031 #include "Chain.h" 00032 #include "MCAlgorithm.h" 00033 #include "LikelihoodFunction.h" 00034 00035 namespace BIE { 00036 00037 // forward declaration 00038 class LikelikehoodCompuation; 00039 class MCAlgorithm; 00040 00041 //+ CLICLASS Simulation 00046 class Simulation: public Serializable { 00047 friend class LikelihoodComputation; 00048 friend class MCAlgorithm; 00049 friend class MetropolisHastings; 00050 friend class ReversibleJump; 00051 00052 private: 00053 00054 bool first_state; 00055 int mstat1, mstat0; 00056 00057 protected: 00059 bool first_log; 00061 int M; 00063 int Ndim; 00064 00066 MHWidth* _mhwidth; 00068 Simulation* _last; 00070 BaseDataTree* _dist; 00072 Integration* _intgr; 00073 // Distribution* @autopersist(_stat); 00074 // Distribution* @autopersist(_dataStat); 00076 MixturePrior* _prior; 00078 Model* _model; 00080 Converge* _conv; 00082 LikelihoodComputation* _likelihoodComputation; 00084 MCAlgorithm* _mca; 00085 00087 Frontier* _currentFrontier; 00089 Frontier* _priorFrontier; 00090 00095 bool internalLP; 00096 00103 bool serialLP; 00104 00105 00107 bool state_not_initialized; 00108 00110 00111 00112 Chain ch, *chfid; 00114 vector<Chain> chains; 00116 int Mchains; 00118 00120 int nbad; 00121 00123 int ngood; 00124 00126 virtual void MCMethod(); 00127 00129 void print_bin_frequency(void); 00130 00132 double threshval; 00133 00139 double BinnedLikeProb(vector<double> &z, SampleDistribution* sd, 00140 double norm, Tile *t, State *s); 00141 00143 double PointLikeProb(vector<double> &z, SampleDistribution* sd, 00144 double norm, Tile *t, State *s); 00145 00147 bool debug_state; 00148 00150 string debug_tag; 00151 00153 int debug_cnt; 00154 00156 string algoname; 00157 00158 public: 00159 00161 static unsigned int state_diag; 00162 00164 static unsigned int output_prec; 00165 00190 //+ CLICONSTR int int MHWidth* BaseDataTree* Model* Integration* Converge* MixturePrior* LikelihoodComputation* MCAlgorithm* 00192 //+ CLICONSTR int int MHWidth* BaseDataTree* Model* Integration* Converge* MixturePrior* LikelihoodComputation* MCAlgorithm* Simulation* 00194 Simulation(int M, int Ndim, MHWidth* mhwidth, 00195 BaseDataTree* d, Model *m, 00196 Integration* i, Converge* c, MixturePrior* p, 00197 LikelihoodComputation *l, MCAlgorithm *mca, 00198 Simulation* last=NULL); 00199 00200 //+ CLICONSTR int int BaseDataTree* Model* Integration* Converge* MixturePrior* LikelihoodComputation* MCAlgorithm* 00202 //+ CLICONSTR int int BaseDataTree* Model* Integration* Converge* MixturePrior* LikelihoodComputation* MCAlgorithm* Simulation* 00204 Simulation(int M, int Ndim, 00205 BaseDataTree* d, Model *m, 00206 Integration* i, Converge* c, MixturePrior* p, 00207 LikelihoodComputation *l, MCAlgorithm *mca, 00208 Simulation* last=NULL); 00210 00212 virtual ~Simulation(); 00213 00215 virtual void Reinitialize(MHWidth* width, MixturePrior *p); 00216 00218 void OneStep(); 00219 00221 Converge *getConverge() const { return _conv; } 00222 00224 bool Convergence() { return _conv->Converged(); } 00225 00227 int BurnIn() { return _conv->ConvergedIndex(); } 00228 00230 00231 // Single chain 00232 bool AccumConverge(double value, State& state) 00233 { return _conv->AccumData(value, state); } 00234 00235 // Parallel chains 00236 bool AccumConverge(vector<double>& values, vector<State>& states) 00237 { return _conv->AccumData(values, states); } 00239 00241 void NewConvergence(int m, Ensemble* d, string id); 00242 00244 State Prior(); 00245 00247 virtual double Likelihood(State&, int indx); 00248 00250 double PriorValue(State&); 00251 00253 virtual State GetState() { return ch.GetState0(); } 00254 00257 virtual void PrintStepDiagnostic(); 00258 00261 virtual void PrintStateDiagnostic(); 00262 00264 00265 virtual void ReportState(State& s); 00266 virtual void ReportState(); 00268 00270 00271 virtual void LogState(string& logfile); 00272 virtual void LogState(int level, int iterno, string& outfile); 00274 00276 virtual double GetValue() { return ch.value; } 00277 00279 virtual double GetPrior() { return ch.prior_value; } 00280 00282 virtual double GetLikelihood() { return ch.like_value; } 00283 00284 //+ CLIMETHOD void SetThreshold double 00286 void SetThreshold(double v) { LP->threshval = v; } 00287 00288 //+ CLIMETHOD void SetUserLikelihood LikelihoodFunction* 00290 void SetUserLikelihood(LikelihoodFunction *fct); 00291 00293 Simulation* GetPreviousSimulation() { return _last; } 00294 00296 void ResetFrontier() { 00297 delete _currentFrontier; 00298 _currentFrontier = _dist->GetDefaultFrontier()->Copy(); 00299 } 00300 00302 void SetPriorFrontier(Frontier *f) { _priorFrontier = f->Copy(); } 00303 00305 int GetMmax() { return M; } 00306 00308 int GetMcur() { return chfid->Mcur; } 00309 00311 int GetNdim() { return Ndim; } 00312 00314 BaseDataTree* GetDataTree() { return _dist; } 00315 00317 LikelihoodComputation* GetLike() { return _likelihoodComputation; } 00318 00320 MCAlgorithm* GetMCA() { return _mca; } 00321 00323 MixturePrior* GetMixturePrior() { return _prior; } 00324 00326 MHWidth* GetMHWidth() { return _mhwidth; } 00327 00328 //+ CLIMETHOD Frontier* currentFrontier 00330 Frontier *currentFrontier() { return _currentFrontier; } 00331 00333 Frontier *priorFrontier() { return _priorFrontier; } 00334 00336 00337 virtual void NewState(int Mcur, vector<double>& w, vector< vector<double> >& p); 00338 virtual void NewState(Chain& ch); 00339 virtual void NewState(State& s); 00341 00342 00344 00345 State Mix_to_State(int Mcur, vector<double>& w, vector< vector<double> >& p); 00346 void State_to_Mix(int& Mcur, vector<double>& w, vector< vector<double> >& p, State& s); 00347 00348 virtual void PrintState() { chfid->PrintState(); } 00349 00350 bool Serial() { return serialLP; } 00351 bool Parallel() { return !serialLP; } 00353 00355 virtual void AdditionalInfo() {} 00356 00363 LikelihoodFunction * LP; 00364 00365 //+ CLIMETHOD void DebugState string 00367 void DebugState(string s) { 00368 debug_state = true; 00369 debug_tag = nametag + "." + s; 00370 } 00371 00372 private: 00373 00374 // Save the RNG 00375 // 00376 template<class Archive> 00377 void post_save(Archive &ar, const unsigned int file_version) const { 00378 vector<Serializable*> s; 00379 s.push_back(BIEgen); 00380 ar << BOOST_SERIALIZATION_NVP(s); 00381 } 00382 // Reassign pointer to fiducial chain and restore the RNG 00383 // 00384 template<class Archive> 00385 void post_load(Archive &ar, const unsigned int file_version) { 00386 vector<Serializable *>s; 00387 ar >> BOOST_SERIALIZATION_NVP(s); 00388 BIEgen = dynamic_cast<BIEACG*>(s[0]); 00389 chfid = &ch; 00390 } 00391 00392 // AUTO GENERATED BY ../persistence/autopersist.py 00393 protected: 00394 Simulation() {} 00395 private: 00396 friend class boost::serialization::access; 00397 BOOST_SERIALIZATION_SPLIT_MEMBER(); 00398 00399 template<class Archive> 00400 void save(Archive & ar, const unsigned int version) const { 00401 this->pre_save(ar, version); 00402 try { 00403 ar << BOOST_SERIALIZATION_BASE_OBJECT_NVP(Serializable); 00404 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00405 } 00406 try { 00407 ar << BOOST_SERIALIZATION_NVP(first_state); 00408 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00409 } 00410 try { 00411 ar << BOOST_SERIALIZATION_NVP(mstat1); 00412 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00413 } 00414 try { 00415 ar << BOOST_SERIALIZATION_NVP(mstat0); 00416 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00417 } 00418 try { 00419 ar << BOOST_SERIALIZATION_NVP(first_log); 00420 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00421 } 00422 try { 00423 ar << BOOST_SERIALIZATION_NVP(M); 00424 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00425 } 00426 try { 00427 ar << BOOST_SERIALIZATION_NVP(Ndim); 00428 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00429 } 00430 try { 00431 ar << BOOST_SERIALIZATION_NVP(_mhwidth); 00432 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00433 } 00434 try { 00435 ar << BOOST_SERIALIZATION_NVP(_last); 00436 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00437 } 00438 try { 00439 ar << BOOST_SERIALIZATION_NVP(_dist); 00440 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00441 } 00442 try { 00443 ar << BOOST_SERIALIZATION_NVP(_intgr); 00444 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00445 } 00446 try { 00447 ar << BOOST_SERIALIZATION_NVP(_prior); 00448 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00449 } 00450 try { 00451 ar << BOOST_SERIALIZATION_NVP(_model); 00452 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00453 } 00454 try { 00455 ar << BOOST_SERIALIZATION_NVP(_conv); 00456 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00457 } 00458 try { 00459 ar << BOOST_SERIALIZATION_NVP(_likelihoodComputation); 00460 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00461 } 00462 try { 00463 ar << BOOST_SERIALIZATION_NVP(_mca); 00464 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00465 } 00466 try { 00467 ar << BOOST_SERIALIZATION_NVP(_currentFrontier); 00468 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00469 } 00470 try { 00471 ar << BOOST_SERIALIZATION_NVP(_priorFrontier); 00472 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00473 } 00474 try { 00475 ar << BOOST_SERIALIZATION_NVP(internalLP); 00476 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00477 } 00478 try { 00479 ar << BOOST_SERIALIZATION_NVP(serialLP); 00480 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00481 } 00482 try { 00483 ar << BOOST_SERIALIZATION_NVP(state_not_initialized); 00484 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00485 } 00486 try { 00487 ar << BOOST_SERIALIZATION_NVP(ch); 00488 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00489 } 00490 try { 00491 ar << BOOST_SERIALIZATION_NVP(chains); 00492 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00493 } 00494 try { 00495 ar << BOOST_SERIALIZATION_NVP(Mchains); 00496 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00497 } 00498 try { 00499 ar << BOOST_SERIALIZATION_NVP(nbad); 00500 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00501 } 00502 try { 00503 ar << BOOST_SERIALIZATION_NVP(ngood); 00504 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00505 } 00506 try { 00507 ar << BOOST_SERIALIZATION_NVP(threshval); 00508 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00509 } 00510 try { 00511 ar << BOOST_SERIALIZATION_NVP(debug_state); 00512 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00513 } 00514 try { 00515 ar << BOOST_SERIALIZATION_NVP(debug_tag); 00516 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00517 } 00518 try { 00519 ar << BOOST_SERIALIZATION_NVP(debug_cnt); 00520 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00521 } 00522 try { 00523 ar << BOOST_SERIALIZATION_NVP(algoname); 00524 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00525 } 00526 try { 00527 ar << BOOST_SERIALIZATION_NVP(state_diag); 00528 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00529 } 00530 try { 00531 ar << BOOST_SERIALIZATION_NVP(output_prec); 00532 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00533 } 00534 try { 00535 ar << BOOST_SERIALIZATION_NVP(LP); 00536 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00537 } 00538 this->post_save(ar, version); 00539 } 00540 00541 template<class Archive> 00542 void load(Archive & ar, const unsigned int version) { 00543 this->pre_load(ar, version); 00544 try { 00545 ar >> BOOST_SERIALIZATION_BASE_OBJECT_NVP(Serializable); 00546 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00547 } 00548 try { 00549 ar >> BOOST_SERIALIZATION_NVP(first_state); 00550 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00551 } 00552 try { 00553 ar >> BOOST_SERIALIZATION_NVP(mstat1); 00554 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00555 } 00556 try { 00557 ar >> BOOST_SERIALIZATION_NVP(mstat0); 00558 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00559 } 00560 try { 00561 ar >> BOOST_SERIALIZATION_NVP(first_log); 00562 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00563 } 00564 try { 00565 ar >> BOOST_SERIALIZATION_NVP(M); 00566 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00567 } 00568 try { 00569 ar >> BOOST_SERIALIZATION_NVP(Ndim); 00570 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00571 } 00572 try { 00573 ar >> BOOST_SERIALIZATION_NVP(_mhwidth); 00574 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00575 } 00576 try { 00577 ar >> BOOST_SERIALIZATION_NVP(_last); 00578 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00579 } 00580 try { 00581 ar >> BOOST_SERIALIZATION_NVP(_dist); 00582 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00583 } 00584 try { 00585 ar >> BOOST_SERIALIZATION_NVP(_intgr); 00586 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00587 } 00588 try { 00589 ar >> BOOST_SERIALIZATION_NVP(_prior); 00590 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00591 } 00592 try { 00593 ar >> BOOST_SERIALIZATION_NVP(_model); 00594 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00595 } 00596 try { 00597 ar >> BOOST_SERIALIZATION_NVP(_conv); 00598 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00599 } 00600 try { 00601 ar >> BOOST_SERIALIZATION_NVP(_likelihoodComputation); 00602 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00603 } 00604 try { 00605 ar >> BOOST_SERIALIZATION_NVP(_mca); 00606 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00607 } 00608 try { 00609 ar >> BOOST_SERIALIZATION_NVP(_currentFrontier); 00610 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00611 } 00612 try { 00613 ar >> BOOST_SERIALIZATION_NVP(_priorFrontier); 00614 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00615 } 00616 try { 00617 ar >> BOOST_SERIALIZATION_NVP(internalLP); 00618 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00619 } 00620 try { 00621 ar >> BOOST_SERIALIZATION_NVP(serialLP); 00622 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00623 } 00624 try { 00625 ar >> BOOST_SERIALIZATION_NVP(state_not_initialized); 00626 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00627 } 00628 try { 00629 ar >> BOOST_SERIALIZATION_NVP(ch); 00630 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00631 } 00632 try { 00633 ar >> BOOST_SERIALIZATION_NVP(chains); 00634 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00635 } 00636 try { 00637 ar >> BOOST_SERIALIZATION_NVP(Mchains); 00638 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00639 } 00640 try { 00641 ar >> BOOST_SERIALIZATION_NVP(nbad); 00642 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00643 } 00644 try { 00645 ar >> BOOST_SERIALIZATION_NVP(ngood); 00646 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00647 } 00648 try { 00649 ar >> BOOST_SERIALIZATION_NVP(threshval); 00650 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00651 } 00652 try { 00653 ar >> BOOST_SERIALIZATION_NVP(debug_state); 00654 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00655 } 00656 try { 00657 ar >> BOOST_SERIALIZATION_NVP(debug_tag); 00658 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00659 } 00660 try { 00661 ar >> BOOST_SERIALIZATION_NVP(debug_cnt); 00662 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00663 } 00664 try { 00665 ar >> BOOST_SERIALIZATION_NVP(algoname); 00666 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00667 } 00668 try { 00669 ar >> BOOST_SERIALIZATION_NVP(state_diag); 00670 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00671 } 00672 try { 00673 ar >> BOOST_SERIALIZATION_NVP(output_prec); 00674 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00675 } 00676 try { 00677 ar >> BOOST_SERIALIZATION_NVP(LP); 00678 BIE_CATCH_BOOST_SERIALIZATION_EXCEPTION; 00679 } 00680 this->post_load(ar, version); 00681 } 00682 00683 }; 00684 } 00685 00686 BIE_CLASS_TYPE_INFO(BIE::Simulation) 00687 BIE_CLASS_EXPORT_KEY(BIE::Simulation) 00688 #endif
Send suggestions, questions, and feedback to WEINBERG at ASTRO dot UMASS dot EDU.
Documentation generated at Fri Mar 26 00:35:10 2010 by
