MaCh3 2.2.1
Reference Guide
Loading...
Searching...
No Matches
AdaptiveMCMCHandler.cpp
Go to the documentation of this file.
2
3namespace adaptive_mcmc{
4
5// ********************************************
7// ********************************************
12 total_steps = 0;
13
14 par_means = {};
15 adaptive_covariance = nullptr;
16}
17
18// ********************************************
20// ********************************************
21 if(adaptive_covariance != nullptr) {
23 }
24}
25
26// ********************************************
27bool AdaptiveMCMCHandler::InitFromConfig(const YAML::Node& adapt_manager, const std::string& matrix_name_str,
28 const std::vector<double>* parameters, const std::vector<double>* fixed) {
29// ********************************************
30 /*
31 * HW: Idea is that adaption can simply read the YAML config
32 * Options :
33 * External Info:
34 * UseExternalMatrix [bool] : Use an external matrix
35 * ExternalMatrixFileName [str] : Name of file containing external info
36 * ExternalMatrixName [str] : Name of external Matrix
37 * ExternalMeansName [str] : Name of external means vector [for updates]
38 *
39 * General Info:
40 * DoAdaption [bool] : Do we want to do adaption?
41 * AdaptionStartThrow [int] : Step we start throwing adaptive matrix from
42 * AdaptionEndUpdate [int] : Step we stop updating adaptive matrix
43 * AdaptionStartUpdate [int] : Do we skip the first N steps?
44 * AdaptionUpdateStep [int] : Number of steps between matrix updates
45 * AdaptionSaveNIterations [int]: You don't have to save every adaptive stage so decide how often you want to save
46 * Adaption blocks [vector<vector<int>>] : Splits the throw matrix into several block matrices
47 * OuputFileName [std::string] : Name of the file that the adaptive matrices will be saved into
48 */
49
50 // setAdaptionDefaults();
51 if(!adapt_manager["AdaptionOptions"]["Covariance"][matrix_name_str]) {
52 MACH3LOG_WARN("Adaptive Settings not found for {}, this is fine if you don't want adaptive MCMC", matrix_name_str);
53 return false;
54 }
55
56 // We"re going to grab this info from the YAML manager
57 if(!GetFromManager<bool>(adapt_manager["AdaptionOptions"]["Covariance"][matrix_name_str]["DoAdaption"], false)) {
58 MACH3LOG_WARN("Not using adaption for {}", matrix_name_str);
59 return false;
60 }
61
62 if(!CheckNodeExists(adapt_manager, "AdaptionOptions", "Settings", "OutputFileName")) {
63 MACH3LOG_ERROR("No OutputFileName specified in AdaptionOptions::Settings into your config file");
64 MACH3LOG_ERROR("This is required if you are using adaptive MCMC");
65 throw MaCh3Exception(__FILE__, __LINE__);
66 }
67
68 start_adaptive_throw = GetFromManager<int>(adapt_manager["AdaptionOptions"]["Settings"]["StartThrow"], 10);
69 start_adaptive_update = GetFromManager<int>(adapt_manager["AdaptionOptions"]["Settings"]["StartUpdate"], 0);
70 end_adaptive_update = GetFromManager<int>(adapt_manager["AdaptionOptions"]["Settings"]["EndUpdate"], 10000);
71 adaptive_update_step = GetFromManager<int>(adapt_manager["AdaptionOptions"]["Settings"]["UpdateStep"], 100);
72 adaptive_save_n_iterations = GetFromManager<int>(adapt_manager["AdaptionOptions"]["Settings"]["SaveNIterations"], 1);
73 output_file_name = GetFromManager<std::string>(adapt_manager["AdaptionOptions"]["Settings"]["OutputFileName"], "");
74
75 // We also want to check for "blocks" by default all parameters "know" about each other
76 // but we can split the matrix into independent block matrices
77 SetParams(parameters);
78 SetFixed(fixed);
79
81 adaption_scale = 2.38*2.38/GetNumParams();
82
83 // We"ll set a dummy variable here
84 auto matrix_blocks = GetFromManager<std::vector<std::vector<int>>>(adapt_manager["AdaptionOptions"]["Covariance"][matrix_name_str]["MatrixBlocks"], {{}});
85
86 SetAdaptiveBlocks(matrix_blocks);
87 return true;
88}
89
90// ********************************************
92// ********************************************
93 adaptive_covariance = new TMatrixDSym(GetNumParams());
94 adaptive_covariance->Zero();
95 par_means = std::vector<double>(GetNumParams(), 0);
96}
97
98// ********************************************
99void AdaptiveMCMCHandler::SetAdaptiveBlocks(const std::vector<std::vector<int>>& block_indices) {
100// ********************************************
101 /*
102 * In order to adapt efficient we want to setup our throw matrix to be a serious of block-diagonal (ish) matrices
103 *
104 * To do this we set sub-block in the config by parameter index. For example having
105 * [[0,4],[4, 6]] in your config will set up two blocks one with all indices 0<=i<4 and the other with 4<=i<6
106 */
107 // Set up block regions
108 adapt_block_matrix_indices = std::vector<int>(GetNumParams(), 0);
109
110 // Should also make a matrix of block sizes
111 adapt_block_sizes = std::vector<int>(block_indices.size()+1, 0);
113
114 if(block_indices.size()==0 || block_indices[0].size()==0) return;
115
116 int block_size = static_cast<int>(block_indices.size());
117 // Now we loop over our blocks
118 for(int iblock=0; iblock < block_size; iblock++){
119 // Loop over blocks in the block
120 int sub_block_size = static_cast<int>(block_indices.size()-1);
121 for(int isubblock=0; isubblock < sub_block_size ; isubblock+=2){
122 int block_lb = block_indices[iblock][isubblock];
123 int block_ub = block_indices[iblock][isubblock+1];
124
125 if(block_lb > GetNumParams() || block_ub > GetNumParams()){
126 MACH3LOG_ERROR("Cannot set matrix block with edges {}, {} for matrix of size {}",
127 block_lb, block_ub, GetNumParams());
128 throw MaCh3Exception(__FILE__, __LINE__);;
129 }
130 for(int ipar = block_lb; ipar < block_ub; ipar++){
131 adapt_block_matrix_indices[ipar] = iblock+1;
132 adapt_block_sizes[iblock+1] += 1;
133 adapt_block_sizes[0] -= 1;
134 }
135 }
136 }
137}
138
139// ********************************************
140//HW: Truly adaptive MCMC!
141void AdaptiveMCMCHandler::SaveAdaptiveToFile(const std::string &outFileName,
142 const std::string &systematicName, bool is_final) {
143// ********************************************
145
146 TFile *outFile = new TFile(outFileName.c_str(), "UPDATE");
147 if (outFile->IsZombie()) {
148 MACH3LOG_ERROR("Couldn't find {}", outFileName);
149 throw MaCh3Exception(__FILE__, __LINE__);
150 }
151
152 TVectorD *outMeanVec = new TVectorD(int(par_means.size()));
153 for (int i = 0; i < int(par_means.size()); i++) {
154 (*outMeanVec)(i) = par_means[i];
155 }
156
157 std::string adaptive_cov_name = systematicName + "_posfit_matrix";
158 std::string mean_vec_name = systematicName + "_mean_vec";
159 if (!is_final) {
160 // Some string to make the name of the saved adaptive matrix clear
161 std::string total_steps_str =
162 std::to_string(total_steps);
163 std::string syst_name_str = systematicName;
164 adaptive_cov_name =
165 total_steps_str + '_' + syst_name_str + std::string("_throw_matrix");
166 mean_vec_name =
167 total_steps_str + '_' + syst_name_str + std::string("_mean_vec");
168 }
169
170 outFile->cd();
171 adaptive_covariance->Write(adaptive_cov_name.c_str());
172 outMeanVec->Write(mean_vec_name.c_str());
173 outFile->Close();
174 delete outMeanVec;
175 delete outFile;
176 }
177}
178
179// ********************************************
180// HW : I would like this to be less painful to use!
181// First things first we need setters
182void AdaptiveMCMCHandler::SetThrowMatrixFromFile(const std::string& matrix_file_name,
183 const std::string& matrix_name,
184 const std::string& means_name,
185 bool& use_adaptive) {
186// ********************************************
187 // Lets you set the throw matrix externally
188 // Open file
189 std::unique_ptr<TFile>matrix_file(new TFile(matrix_file_name.c_str()));
190 use_adaptive = true;
191
192 if(matrix_file->IsZombie()){
193 MACH3LOG_ERROR("Couldn't find {}", matrix_file_name);
194 throw MaCh3Exception(__FILE__ , __LINE__ );
195 }
196
197 // Next we grab our matrix
198 adaptive_covariance = static_cast<TMatrixDSym*>(matrix_file->Get(matrix_name.c_str()));
200 MACH3LOG_ERROR("Couldn't find {} in {}", matrix_name, matrix_file_name);
201 throw MaCh3Exception(__FILE__ , __LINE__ );
202 }
203
204 // Finally we grab the means vector
205 TVectorD* means_vector = static_cast<TVectorD*>(matrix_file->Get(means_name.c_str()));
206
207 // This is fine to not exist!
208 if(means_vector){
209 // Yay our vector exists! Let's loop and fill it
210 // Should check this is done
211 if(means_vector->GetNrows()){
212 MACH3LOG_ERROR("External means vec size ({}) != matrix size ({})", means_vector->GetNrows(), GetNumParams());
213 throw MaCh3Exception(__FILE__, __LINE__);
214 }
215
216 par_means = std::vector<double>(GetNumParams());
217 for(int i = 0; i < GetNumParams(); i++){
218 par_means[i] = (*means_vector)(i);
219 }
220 MACH3LOG_INFO("Found Means in External File, Will be able to adapt");
221 }
222 // Totally fine if it doesn't exist, we just can't do adaption
223 else{
224 // We don't need a means vector, set the adaption=false
225 MACH3LOG_WARN("Cannot find means vector in {}, therefore I will not be able to adapt!", matrix_file_name);
226 use_adaptive = false;
227 }
228
229 matrix_file->Close();
230 MACH3LOG_INFO("Set up matrix from external file");
231}
232
233// ********************************************
235// ********************************************
236 std::vector<double> par_means_prev = par_means;
237 int steps_post_burn = total_steps - start_adaptive_update;
238
239 // Step 1: Update means and compute deviations
240 for (int i = 0; i < GetNumParams(); ++i) {
241 if (IsFixed(i)) continue;
242
243 par_means[i] = (CurrVal(i) + par_means_prev[i]*steps_post_burn)/(steps_post_burn+1);
244 // Left over from cyclic means
245 }
246
247 // Step 2: Update covariance
248 #ifdef MULTITHREAD
249 #pragma omp parallel for
250 #endif
251 for (int i = 0; i < GetNumParams(); ++i) {
252 if (IsFixed(i)) {
253 (*adaptive_covariance)(i, i) = 1.0;
254 continue;
255 }
256
257 int block_i = adapt_block_matrix_indices[i];
258
259 for (int j = 0; j <= i; ++j) {
260 if (IsFixed(j) || adapt_block_matrix_indices[j] != block_i) {
261 (*adaptive_covariance)(i, j) = 0.0;
262 (*adaptive_covariance)(j, i) = 0.0;
263 continue;
264 }
265
266 double cov_prev = (*adaptive_covariance)(i, j);
267 double cov_updated = 0.0;
268
269 if (steps_post_burn > 0) {
270 // Haario-style update
271 double cov_t = cov_prev * (steps_post_burn - 1) / steps_post_burn;
272 double prev_means_t = steps_post_burn * par_means_prev[i] * par_means_prev[j];
273 double curr_means_t = (steps_post_burn + 1) * par_means[i] * par_means[j];
274 double curr_step_t = CurrVal(i) * CurrVal(j);
275
276 cov_updated = cov_t + adaption_scale * (prev_means_t - curr_means_t + curr_step_t) / steps_post_burn;
277 }
278
279 (*adaptive_covariance)(i, j) = cov_updated;
280 (*adaptive_covariance)(j, i) = cov_updated;
281 }
282 }
283}
284
285// ********************************************
287// ********************************************
288 if(total_steps == start_adaptive_throw) return true;
289 else return false;
290}
291
292// ********************************************
294// ********************************************
296 // Check whether the number of steps is divisible by the adaptive update step
297 // e.g. if adaptive_update_step = 1000 and (total_step - start_adpative_throw) is 5000 then this is true
299 return true;
300 }
301 else return false;
302}
303
304// ********************************************
306// ********************************************
309 else return false;
310}
311
312// ********************************************
314// ********************************************
315 if(total_steps <= start_adaptive_throw) return true;
316 else return false;
317}
318
319// ********************************************
321// ********************************************
322 MACH3LOG_INFO("Adaptive MCMC Info:");
323 MACH3LOG_INFO("Throwing from New Matrix from Step : {}", start_adaptive_throw);
324 MACH3LOG_INFO("Adaption Matrix Start Update : {}", start_adaptive_update);
325 MACH3LOG_INFO("Adaption Matrix Ending Updates : {}", end_adaptive_update);
326 MACH3LOG_INFO("Steps Between Updates : {}", adaptive_update_step);
327 MACH3LOG_INFO("Saving matrices to file : {}", output_file_name);
328 MACH3LOG_INFO("Will only save every {} iterations" , adaptive_save_n_iterations);
329}
330
331double AdaptiveMCMCHandler::CurrVal(const int par_index){
334 return (*_fCurrVal)[par_index];
335}
336
337} //end adaptive_mcmc
#define MACH3LOG_ERROR
Definition: MaCh3Logger.h:25
#define MACH3LOG_INFO
Definition: MaCh3Logger.h:23
#define MACH3LOG_WARN
Definition: MaCh3Logger.h:24
bool CheckNodeExists(const YAML::Node &node, Args... args)
KS: Wrapper function to call the recursive helper.
Definition: YamlHelper.h:54
Custom exception class for MaCh3 errors.
std::vector< double > par_means
Mean values for all parameters.
bool IsFixed(const int ipar) const
Check if a parameter is fixed.
bool UpdateMatrixAdapt()
Tell whether matrix should be updated.
std::vector< int > adapt_block_matrix_indices
Indices for block-matrix adaption.
double CurrVal(const int par_index)
Get Current value of parameter.
bool InitFromConfig(const YAML::Node &adapt_manager, const std::string &matrix_name_str, const std::vector< double > *parameters, const std::vector< double > *fixed)
Read initial values from config file.
std::string output_file_name
Name of the file to save the adaptive matrices into.
bool AdaptionUpdate()
To be fair not a clue...
void SaveAdaptiveToFile(const std::string &outFileName, const std::string &systematicName, const bool is_final=false)
HW: Save adaptive throw matrix to file.
void SetAdaptiveBlocks(const std::vector< std::vector< int > > &block_indices)
HW: sets adaptive block matrix.
const std::vector< double > * _fCurrVal
Current values of parameters.
void Print()
Print all class members.
int start_adaptive_update
When do we stop update the adaptive matrix.
std::vector< int > adapt_block_sizes
Size of blocks for adaption.
bool IndivStepScaleAdapt()
Tell whether we want reset step scale or not.
void SetFixed(const std::vector< double > *fix)
Set the fixed parameters.
void SetParams(const std::vector< double > *params)
Set the current values of the parameters.
void CreateNewAdaptiveCovariance()
If we don't have a covariance matrix to start from for adaptive tune we need to make one!
int end_adaptive_update
Steps between changing throw matrix.
int total_steps
Total number of MCMC steps.
void UpdateAdaptiveCovariance()
Method to update adaptive MCMC .
int adaptive_update_step
Steps between changing throw matrix.
void SetThrowMatrixFromFile(const std::string &matrix_file_name, const std::string &matrix_name, const std::string &means_name, bool &use_adaptive)
sets throw matrix from a file
int GetNumParams() const
Get the current values of the parameters.
bool SkipAdaption()
Tell if we are Skipping Adaption.
TMatrixDSym * adaptive_covariance
Full adaptive covariance matrix.