Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2018 The VES code team
3 : (see the PEOPLE-VES file at the root of this folder for a list of names)
4 :
5 : See http://www.ves-code.org for more information.
6 :
7 : This file is part of VES code module.
8 :
9 : The VES code module is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : The VES code module is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with the VES code module. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 :
23 : #include "Optimizer.h"
24 : #include "CoeffsVector.h"
25 : #include "CoeffsMatrix.h"
26 : #include "VesBias.h"
27 : #include "VesTools.h"
28 :
29 : #include "tools/Exception.h"
30 : #include "core/PlumedMain.h"
31 : #include "core/ActionSet.h"
32 : #include "tools/Communicator.h"
33 : #include "tools/File.h"
34 : #include "tools/FileBase.h"
35 :
36 : namespace PLMD {
37 : namespace ves {
38 :
39 71 : Optimizer::Optimizer(const ActionOptions&ao):
40 : Action(ao),
41 : ActionPilot(ao),
42 : ActionWithValue(ao),
43 : description_("Undefined"),
44 : type_("Undefined"),
45 : stepsizes_(0),
46 : current_stepsizes(0),
47 : fixed_stepsize_(true),
48 : iter_counter(0),
49 : use_hessian_(false),
50 : diagonal_hessian_(true),
51 : monitor_instantaneous_gradient_(false),
52 : use_mwalkers_mpi_(false),
53 : mwalkers_mpi_single_files_(true),
54 : dynamic_targetdists_(0),
55 : ustride_targetdist_(0),
56 : ustride_reweightfactor_(0),
57 : coeffssetid_prefix_(""),
58 : coeffs_wstride_(100),
59 : coeffsOFiles_(0),
60 : coeffs_output_fmt_(""),
61 : gradient_wstride_(100),
62 : gradientOFiles_(0),
63 : gradient_output_fmt_(""),
64 : hessian_wstride_(100),
65 : hessianOFiles_(0),
66 : hessian_output_fmt_(""),
67 : targetdist_averages_wstride_(0),
68 : targetdist_averagesOFiles_(0),
69 : targetdist_averages_output_fmt_(""),
70 : nbiases_(0),
71 : bias_pntrs_(0),
72 : ncoeffssets_(0),
73 : coeffs_pntrs_(0),
74 : aux_coeffs_pntrs_(0),
75 : gradient_pntrs_(0),
76 : aver_gradient_pntrs_(0),
77 : hessian_pntrs_(0),
78 : coeffs_mask_pntrs_(0),
79 : targetdist_averages_pntrs_(0),
80 : identical_coeffs_shape_(true),
81 : bias_output_active_(false),
82 : bias_output_stride_(0),
83 : fes_output_active_(false),
84 : fes_output_stride_(0),
85 : fesproj_output_active_(false),
86 : fesproj_output_stride_(0),
87 : targetdist_output_active_(false),
88 : targetdist_output_stride_(0),
89 : targetdist_proj_output_active_(false),
90 : targetdist_proj_output_stride_(0),
91 142 : isFirstStep(true)
92 : {
93 142 : std::vector<std::string> bias_labels(0);
94 142 : parseVector("BIAS",bias_labels);
95 71 : plumed_massert(bias_labels.size()>0,"problem with BIAS keyword");
96 71 : nbiases_ = bias_labels.size();
97 : //
98 71 : std::string error_msg = "";
99 213 : bias_pntrs_ = VesTools::getPointersFromLabels<VesBias*>(bias_labels,plumed.getActionSet(),error_msg);
100 71 : if(error_msg.size()>0) {plumed_merror("Error in keyword BIAS of "+getName()+": "+error_msg);}
101 :
102 373 : for(unsigned int i=0; i<bias_pntrs_.size(); i++) {
103 77 : bias_pntrs_[i]->linkOptimizer(this);
104 : //
105 77 : std::vector<CoeffsVector*> pntrs_coeffs = bias_pntrs_[i]->getCoeffsPntrs();
106 77 : std::vector<CoeffsVector*> pntrs_gradient = bias_pntrs_[i]->getGradientPntrs();
107 77 : std::vector<CoeffsVector*> pntrs_targetdist_averages = bias_pntrs_[i]->getTargetDistAveragesPntrs();
108 77 : plumed_massert(pntrs_coeffs.size()==pntrs_gradient.size(),"something wrong in the coefficients and gradient passed from VES bias");
109 77 : plumed_massert(pntrs_coeffs.size()==pntrs_targetdist_averages.size(),"something wrong in the coefficients and target distribution averages passed from VES bias");
110 385 : for(unsigned int k=0; k<pntrs_coeffs.size(); k++) {
111 77 : plumed_massert(pntrs_coeffs[k] != NULL,"some coefficient is not linked correctly");
112 77 : plumed_massert(pntrs_gradient[k] != NULL,"some gradient is not linked correctly");
113 77 : plumed_massert(pntrs_targetdist_averages[k] != NULL,"some target distribution average is not linked correctly");
114 : pntrs_coeffs[k]->turnOnIterationCounter();
115 77 : coeffs_pntrs_.push_back(pntrs_coeffs[k]);
116 77 : pntrs_gradient[k]->turnOnIterationCounter();
117 77 : gradient_pntrs_.push_back(pntrs_gradient[k]);
118 77 : pntrs_targetdist_averages[k]->turnOnIterationCounter();
119 77 : targetdist_averages_pntrs_.push_back(pntrs_targetdist_averages[k]);
120 : //
121 77 : CoeffsVector* aux_coeffs_tmp = new CoeffsVector(*pntrs_coeffs[k]);
122 77 : std::string aux_label = pntrs_coeffs[k]->getLabel();
123 77 : if(aux_label.find("coeffs")!=std::string::npos) {
124 231 : aux_label.replace(aux_label.find("coeffs"), std::string("coeffs").length(), "aux_coeffs");
125 : }
126 : else {
127 : aux_label += "_aux";
128 : }
129 77 : aux_coeffs_tmp->setLabels(aux_label);
130 77 : aux_coeffs_pntrs_.push_back(aux_coeffs_tmp);
131 77 : AuxCoeffs(i).setValues( Coeffs(i) );
132 : }
133 : }
134 71 : ncoeffssets_ = coeffs_pntrs_.size();
135 71 : plumed_massert(aux_coeffs_pntrs_.size()==ncoeffssets_,"problems in linking aux coefficients");
136 71 : plumed_massert(gradient_pntrs_.size()==ncoeffssets_,"problems in linking gradients");
137 71 : setAllCoeffsSetIterationCounters();
138 :
139 :
140 : //
141 71 : identical_coeffs_shape_ = true;
142 83 : for(unsigned int i=1; i<ncoeffssets_; i++) {
143 12 : if(!coeffs_pntrs_[0]->sameShape(*coeffs_pntrs_[i])) {
144 0 : identical_coeffs_shape_ = false;
145 0 : break;
146 : }
147 : }
148 : //
149 142 : if(keywords.exists("STEPSIZE")) {
150 140 : plumed_assert(!keywords.exists("INITIAL_STEPSIZE"));
151 70 : fixed_stepsize_=true;
152 140 : parseMultipleValues("STEPSIZE",stepsizes_);
153 70 : setCurrentStepSizes(stepsizes_);
154 : }
155 142 : if(keywords.exists("INITIAL_STEPSIZE")) {
156 0 : plumed_assert(!keywords.exists("STEPSIZE"));
157 0 : fixed_stepsize_=false;
158 0 : parseMultipleValues("INITIAL_STEPSIZE",stepsizes_);
159 0 : setCurrentStepSizes(stepsizes_);
160 : }
161 : //
162 71 : if(ncoeffssets_==1) {
163 204 : log.printf(" optimizing VES bias %s with label %s: \n",bias_pntrs_[0]->getName().c_str(),bias_pntrs_[0]->getLabel().c_str());
164 204 : log.printf(" KbT: %f\n",bias_pntrs_[0]->getKbT());
165 204 : log.printf(" number of coefficients: %zu\n",coeffs_pntrs_[0]->numberOfCoeffs());
166 68 : if(stepsizes_.size()>0) {
167 67 : if(fixed_stepsize_) {log.printf(" using a constant step size of %f\n",stepsizes_[0]);}
168 0 : else {log.printf(" using an initial step size of %f\n",stepsizes_[0]);}
169 : }
170 : }
171 : else {
172 3 : log.printf(" optimizing %u coefficent sets from following %u VES biases:\n",ncoeffssets_,nbiases_);
173 21 : for(unsigned int i=0; i<nbiases_; i++) {
174 36 : log.printf(" %s of type %s (KbT: %f) \n",bias_pntrs_[i]->getLabel().c_str(),bias_pntrs_[i]->getName().c_str(),bias_pntrs_[i]->getKbT());
175 : }
176 : size_t tot_ncoeffs = 0;
177 21 : for(unsigned int i=0; i<ncoeffssets_; i++) {
178 9 : log.printf(" coefficient set %u: \n",i);
179 27 : log.printf(" used in bias %s (type %s)\n",coeffs_pntrs_[i]->getPntrToAction()->getLabel().c_str(),coeffs_pntrs_[i]->getPntrToAction()->getName().c_str());
180 27 : log.printf(" number of coefficients: %zu\n",coeffs_pntrs_[i]->numberOfCoeffs());
181 9 : if(stepsizes_.size()>0) {
182 18 : if(fixed_stepsize_) {log.printf(" using a constant step size of %f\n",stepsizes_[i]);}
183 0 : else {log.printf(" using an initial step size of %f\n",stepsizes_[i]);}
184 : }
185 18 : tot_ncoeffs += coeffs_pntrs_[i]->numberOfCoeffs();
186 : }
187 3 : log.printf(" total number of coefficients: %zu\n",tot_ncoeffs);
188 3 : if(identical_coeffs_shape_) {
189 3 : log.printf(" the indices shape is identical for all coefficient sets\n");
190 : }
191 : else {
192 0 : log.printf(" the indices shape differs between coefficient sets\n");
193 : }
194 : }
195 :
196 : //
197 142 : if(keywords.exists("FULL_HESSIAN")) {
198 0 : bool full_hessian=false;
199 0 : parseFlag("FULL_HESSIAN",full_hessian);
200 0 : diagonal_hessian_ = !full_hessian;
201 : }
202 : //
203 : bool mw_single_files = false;
204 142 : if(keywords.exists("MULTIPLE_WALKERS")) {
205 142 : parseFlag("MULTIPLE_WALKERS",use_mwalkers_mpi_);
206 71 : if(use_mwalkers_mpi_) {
207 : mw_single_files=true;
208 : }
209 : }
210 :
211 71 : int numwalkers=1;
212 71 : int walker_rank=0;
213 71 : if(comm.Get_rank()==0) {
214 62 : numwalkers = multi_sim_comm.Get_size();
215 62 : walker_rank = multi_sim_comm.Get_rank();
216 : }
217 71 : comm.Bcast(numwalkers,0);
218 71 : comm.Bcast(walker_rank,0);
219 71 : if(use_mwalkers_mpi_ && numwalkers==1) {
220 0 : plumed_merror("using the MULTIPLE_WALKERS keyword does not make sense if running the MD code with a single replica");
221 : }
222 71 : if(use_mwalkers_mpi_) {
223 12 : log.printf(" optimization performed using multiple walkers connected via MPI:\n");
224 12 : log.printf(" number of walkers: %d\n",numwalkers);
225 12 : log.printf(" walker number: %d\n",walker_rank);
226 12 : log.printf(" please see and cite ");
227 36 : log << plumed.cite("Raiteri, Laio, Gervasio, Micheletti, and Parrinello, J. Phys. Chem. B 110, 3533 (2006)");
228 12 : log.printf("\n");
229 : }
230 :
231 71 : dynamic_targetdists_.resize(nbiases_,false);
232 142 : if(keywords.exists("TARGETDIST_STRIDE")) {
233 : bool need_stride = false;
234 222 : for(unsigned int i=0; i<nbiases_; i++) {
235 152 : dynamic_targetdists_[i] = bias_pntrs_[i]->dynamicTargetDistribution();
236 76 : if(dynamic_targetdists_[i]) {need_stride = true;}
237 : }
238 140 : parse("TARGETDIST_STRIDE",ustride_targetdist_);
239 70 : if(need_stride && ustride_targetdist_==0) {
240 0 : plumed_merror("one of the biases has a dynamic target distribution so you need to give stride for updating it by using the TARGETDIST_STRIDE keyword");
241 : }
242 70 : if(!need_stride && ustride_targetdist_!=0) {
243 0 : plumed_merror("using the TARGETDIST_STRIDE keyword doesn't make sense as there is no dynamic target distribution to update");
244 : }
245 70 : if(ustride_targetdist_>0) {
246 33 : if(nbiases_==1) {
247 33 : log.printf(" the target distribution will be updated very %u coefficent iterations\n",ustride_targetdist_);
248 : }
249 : else {
250 0 : log.printf(" the target distribution will be updated very %u coefficent iterations for the following biases\n ",ustride_targetdist_);
251 0 : for(unsigned int i=0; i<nbiases_; i++) {
252 0 : log.printf("%s ",bias_pntrs_[i]->getLabel().c_str());
253 : }
254 0 : log.printf("\n");
255 : }
256 33 : log.printf(" See and cite ");
257 99 : log << plumed.cite("Valsson and Parrinello, J. Chem. Theory Comput. 11, 1996-2002 (2015)");
258 33 : log.printf("\n");
259 : }
260 : }
261 :
262 142 : if(keywords.exists("REWEIGHT_FACTOR_STRIDE")) {
263 : bool reweightfactor_calculated = false;
264 0 : for(unsigned int i=0; i<nbiases_; i++) {
265 0 : reweightfactor_calculated = bias_pntrs_[i]->isReweightFactorCalculated();
266 : }
267 0 : parse("REWEIGHT_FACTOR_STRIDE",ustride_reweightfactor_);
268 0 : if(ustride_reweightfactor_==0 && reweightfactor_calculated) {
269 0 : plumed_merror("the calculation of the reweight factor is enabled, You need to use the REWEIGHT_FACTOR_STRIDE keyword to specfiy how often it should be updated.");
270 : }
271 0 : if(ustride_reweightfactor_>0) {
272 0 : if(!reweightfactor_calculated) {
273 0 : plumed_merror("In order to use the REWEIGHT_FACTOR_STRIDE keyword you need to enable the calculation of the reweight factor in the VES bias by using the CALC_REWEIGHT_FACTOR flag.");
274 : }
275 0 : log.printf(" the reweight factor c(t) will be updated very %u coefficent iterations\n",ustride_reweightfactor_);
276 : }
277 : }
278 :
279 142 : if(keywords.exists("MONITOR_INSTANTANEOUS_GRADIENT")) {
280 142 : parseFlag("MONITOR_INSTANTANEOUS_GRADIENT",monitor_instantaneous_gradient_);
281 : }
282 :
283 142 : if(keywords.exists("MONITOR_AVERAGE_GRADIENT")) {
284 71 : bool monitor_aver_gradient = false;
285 142 : parseFlag("MONITOR_AVERAGE_GRADIENT",monitor_aver_gradient);
286 71 : if(monitor_aver_gradient) {
287 2 : unsigned int averaging_exp_decay=0;
288 4 : parse("MONITOR_AVERAGES_GRADIENT_EXP_DECAY",averaging_exp_decay);
289 : aver_gradient_pntrs_.clear();
290 6 : for(unsigned int i=0; i<ncoeffssets_; i++) {
291 4 : CoeffsVector* aver_gradient_tmp = new CoeffsVector(*gradient_pntrs_[i]);
292 2 : aver_gradient_tmp->clear();
293 2 : std::string aver_grad_label = aver_gradient_tmp->getLabel();
294 2 : if(aver_grad_label.find("gradient")!=std::string::npos) {
295 6 : aver_grad_label.replace(aver_grad_label.find("gradient"), std::string("gradient").length(), "aver_gradient");
296 : }
297 : else {
298 : aver_grad_label += "_aver";
299 : }
300 2 : aver_gradient_tmp->setLabels(aver_grad_label);
301 2 : if(averaging_exp_decay>0) {
302 1 : aver_gradient_tmp->setupExponentiallyDecayingAveraging(averaging_exp_decay);
303 : }
304 2 : aver_gradient_pntrs_.push_back(aver_gradient_tmp);
305 : }
306 : }
307 : }
308 :
309 :
310 71 : if(ncoeffssets_>1) {
311 : coeffssetid_prefix_="c-";
312 6 : if(keywords.exists("COEFFS_SET_ID_PREFIX")) {
313 6 : parse("COEFFS_SET_ID_PREFIX",coeffssetid_prefix_);
314 : }
315 : }
316 : else {
317 : coeffssetid_prefix_="";
318 136 : if(keywords.exists("COEFFS_SET_ID_PREFIX")) {
319 136 : parse("COEFFS_SET_ID_PREFIX",coeffssetid_prefix_);
320 : }
321 68 : if(coeffssetid_prefix_.size()>0) {
322 0 : plumed_merror("COEFFS_SET_ID_PREFIX should only be given if optimizing multiple coefficent sets");
323 : }
324 : }
325 :
326 142 : if(keywords.exists("INITIAL_COEFFS")) {
327 71 : std::vector<std::string> initial_coeffs_fnames;
328 142 : parseFilenames("INITIAL_COEFFS",initial_coeffs_fnames);
329 71 : if(initial_coeffs_fnames.size()>0) {
330 1 : readCoeffsFromFiles(initial_coeffs_fnames,false);
331 1 : comm.Barrier();
332 1 : if(comm.Get_rank()==0 && use_mwalkers_mpi_) {
333 0 : multi_sim_comm.Barrier();
334 : }
335 1 : setAllCoeffsSetIterationCounters();
336 : }
337 : }
338 : //
339 :
340 71 : std::vector<std::string> coeffs_fnames;
341 142 : if(keywords.exists("COEFFS_FILE")) {
342 213 : parseFilenames("COEFFS_FILE",coeffs_fnames,"coeffs.data");
343 71 : bool start_opt_afresh=false;
344 142 : if(keywords.exists("START_OPTIMIZATION_AFRESH")) {
345 140 : parseFlag("START_OPTIMIZATION_AFRESH",start_opt_afresh);
346 71 : if(start_opt_afresh && !getRestart()) {
347 0 : plumed_merror("the START_OPTIMIZATION_AFRESH keyword should only be used when a restart has been triggered by the RESTART keyword or the MD code");
348 : }
349 : }
350 142 : if(getRestart()) {
351 71 : for(unsigned int i=0; i<coeffs_fnames.size(); i++) {
352 30 : IFile ifile;
353 15 : ifile.link(*this);
354 19 : if(use_mwalkers_mpi_) {ifile.enforceSuffix("");}
355 15 : bool file_exist = ifile.FileExist(coeffs_fnames[i]);
356 15 : if(!file_exist) {
357 0 : std::string fname = FileBase::appendSuffix(coeffs_fnames[i],ifile.getSuffix());
358 0 : plumed_merror("Cannot find coefficient file " + fname + " when trying to restart an optimzation. If you don't want to restart the optimzation please remove the RESTART keyword or use the RESTART=NO within the "+getName()+" action to locally disable the restart.");
359 : }
360 : }
361 13 : readCoeffsFromFiles(coeffs_fnames,true);
362 13 : comm.Barrier();
363 13 : if(comm.Get_rank()==0 && use_mwalkers_mpi_) {
364 4 : multi_sim_comm.Barrier();
365 : }
366 13 : unsigned int iter_opt_tmp = coeffs_pntrs_[0]->getIterationCounter();
367 17 : for(unsigned int i=1; i<ncoeffssets_; i++) {
368 6 : plumed_massert(coeffs_pntrs_[i]->getIterationCounter()==iter_opt_tmp,"the iteraton counter should be the same for all files when restarting from previous coefficient files\n");
369 : }
370 13 : if(start_opt_afresh) {
371 : setIterationCounter(0);
372 1 : log.printf(" Optimization started afresh at iteration %u\n",getIterationCounter());
373 3 : for(unsigned int i=0; i<ncoeffssets_; i++) {
374 1 : AuxCoeffs(i).setValues( Coeffs(i) );
375 : }
376 : }
377 : else {
378 : setIterationCounter(coeffs_pntrs_[0]->getIterationCounter());
379 12 : log.printf(" Optimization restarted at iteration %u\n",getIterationCounter());
380 : }
381 13 : setAllCoeffsSetIterationCounters();
382 : }
383 :
384 71 : std::string coeffs_wstride_tmpstr="";
385 142 : parse("COEFFS_OUTPUT",coeffs_wstride_tmpstr);
386 142 : if(coeffs_wstride_tmpstr!="OFF" && coeffs_wstride_tmpstr.size()>0) {
387 71 : Tools::convert(coeffs_wstride_tmpstr,coeffs_wstride_);
388 : }
389 71 : if(coeffs_wstride_tmpstr=="OFF") {
390 0 : coeffs_fnames.clear();
391 : }
392 71 : setupOFiles(coeffs_fnames,coeffsOFiles_,mw_single_files);
393 142 : parse("COEFFS_FMT",coeffs_output_fmt_);
394 71 : if(coeffs_output_fmt_.size()>0) {
395 222 : for(unsigned int i=0; i<ncoeffssets_; i++) {
396 152 : coeffs_pntrs_[i]->setOutputFmt(coeffs_output_fmt_);
397 : }
398 : }
399 142 : if(!getRestart()) {
400 302 : for(unsigned int i=0; i<coeffsOFiles_.size(); i++) {
401 186 : coeffs_pntrs_[i]->writeToFile(*coeffsOFiles_[i],aux_coeffs_pntrs_[i],false);
402 : }
403 : }
404 71 : if(coeffs_fnames.size()>0) {
405 71 : if(ncoeffssets_==1) {
406 272 : log.printf(" Coefficients will be written out to file %s every %u iterations\n",coeffsOFiles_[0]->getPath().c_str(),coeffs_wstride_);
407 : }
408 : else {
409 3 : log.printf(" Coefficients will be written out to the following files every %u iterations:\n",coeffs_wstride_);
410 33 : for(unsigned int i=0; i<coeffs_fnames.size(); i++) {
411 36 : log.printf(" coefficient set %u: %s\n",i,coeffsOFiles_[i]->getPath().c_str());
412 : }
413 : }
414 : }
415 : else {
416 0 : log.printf(" Output of coefficients to file has been disabled\n");
417 : }
418 : }
419 :
420 71 : std::vector<std::string> gradient_fnames;
421 142 : if(keywords.exists("GRADIENT_FILE")) {
422 142 : parseFilenames("GRADIENT_FILE",gradient_fnames);
423 142 : parse("GRADIENT_OUTPUT",gradient_wstride_);
424 :
425 71 : if(coeffs_fnames.size()>0) {
426 352 : for(unsigned int i=0; i<gradient_fnames.size(); i++) {
427 70 : plumed_massert(gradient_fnames[i]!=coeffs_fnames[i],"COEFFS_FILE and GRADIENT_FILE cannot be the same");
428 : }
429 : }
430 71 : setupOFiles(gradient_fnames,gradientOFiles_,mw_single_files);
431 142 : parse("GRADIENT_FMT",gradient_output_fmt_);
432 71 : if(gradient_output_fmt_.size()>0) {
433 204 : for(unsigned int i=0; i<ncoeffssets_; i++) {
434 140 : gradient_pntrs_[i]->setOutputFmt(gradient_output_fmt_);
435 : }
436 : }
437 :
438 71 : if(gradient_fnames.size()>0) {
439 64 : if(ncoeffssets_==1) {
440 244 : log.printf(" Gradient will be written out to file %s every %u iterations\n",gradientOFiles_[0]->getPath().c_str(),gradient_wstride_);
441 : }
442 : else {
443 3 : log.printf(" Gradient will be written out to the following files every %u iterations:\n",gradient_wstride_);
444 33 : for(unsigned int i=0; i<gradient_fnames.size(); i++) {
445 36 : log.printf(" coefficient set %u: %s\n",i,gradientOFiles_[i]->getPath().c_str());
446 : }
447 : }
448 : }
449 : }
450 :
451 71 : std::vector<std::string> hessian_fnames;
452 142 : if(keywords.exists("HESSIAN_FILE")) {
453 142 : parseFilenames("HESSIAN_FILE",hessian_fnames);
454 142 : parse("HESSIAN_OUTPUT",hessian_wstride_);
455 :
456 71 : if(coeffs_fnames.size()>0) {
457 337 : for(unsigned int i=0; i<hessian_fnames.size(); i++) {
458 65 : plumed_massert(hessian_fnames[i]!=coeffs_fnames[i],"COEFFS_FILE and HESSIAN_FILE cannot be the same");
459 : }
460 : }
461 71 : if(gradient_fnames.size()>0) {
462 323 : for(unsigned int i=0; i<hessian_fnames.size(); i++) {
463 65 : plumed_massert(hessian_fnames[i]!=gradient_fnames[i],"GRADIENT_FILE and HESSIAN_FILE cannot be the same");
464 : }
465 : }
466 71 : setupOFiles(hessian_fnames,hessianOFiles_,mw_single_files);
467 142 : parse("HESSIAN_FMT",hessian_output_fmt_);
468 :
469 71 : if(hessian_fnames.size()>0) {
470 61 : if(ncoeffssets_==1) {
471 236 : log.printf(" Hessian will be written out to file %s every %u iterations\n",hessianOFiles_[0]->getPath().c_str(),hessian_wstride_);
472 : }
473 : else {
474 2 : log.printf(" Hessian will be written out to the following files every %u iterations:\n",hessian_wstride_);
475 22 : for(unsigned int i=0; i<hessian_fnames.size(); i++) {
476 24 : log.printf(" coefficient set %u: %s\n",i,hessianOFiles_[i]->getPath().c_str());
477 : }
478 : }
479 : }
480 : }
481 :
482 :
483 : //
484 142 : if(keywords.exists("MASK_FILE")) {
485 70 : std::vector<std::string> mask_fnames_in;
486 140 : parseVector("MASK_FILE",mask_fnames_in);
487 70 : if(mask_fnames_in.size()==1 && ncoeffssets_>1) {
488 0 : if(identical_coeffs_shape_) {mask_fnames_in.resize(ncoeffssets_,mask_fnames_in[0]);}
489 0 : else {plumed_merror("the coefficients indices shape differs between biases so you need to give a separate file for each coefficient set\n");}
490 : }
491 70 : if(mask_fnames_in.size()>0 && mask_fnames_in.size()!=ncoeffssets_) {
492 0 : plumed_merror("Error in MASK_FILE keyword: either give one value for all biases or a separate value for each coefficient set");
493 : }
494 :
495 70 : coeffs_mask_pntrs_.resize(ncoeffssets_);
496 222 : for(unsigned int i=0; i<ncoeffssets_; i++) {
497 228 : coeffs_mask_pntrs_[i] = new CoeffsVector(*coeffs_pntrs_[i]);
498 152 : coeffs_mask_pntrs_[i]->setLabels("mask");
499 76 : coeffs_mask_pntrs_[i]->setValues(1.0);
500 152 : coeffs_mask_pntrs_[i]->setOutputFmt("%f");
501 : }
502 :
503 70 : if(mask_fnames_in.size()>0) {
504 1 : if(ncoeffssets_==1) {
505 1 : size_t nread = coeffs_mask_pntrs_[0]->readFromFile(mask_fnames_in[0],true,true);
506 2 : log.printf(" read %zu values from mask file %s\n",nread,mask_fnames_in[0].c_str());
507 1 : size_t ndeactived = coeffs_mask_pntrs_[0]->countValues(0.0);
508 1 : log.printf(" deactived optimization of %zu coefficients\n",ndeactived);
509 : }
510 : else {
511 0 : for(unsigned int i=0; i<ncoeffssets_; i++) {
512 0 : size_t nread = coeffs_mask_pntrs_[i]->readFromFile(mask_fnames_in[i],true,true);
513 0 : log.printf(" mask for coefficent set %u:\n",i);
514 0 : log.printf(" read %zu values from file %s\n",nread,mask_fnames_in[i].c_str());
515 0 : size_t ndeactived = coeffs_mask_pntrs_[0]->countValues(0.0);
516 0 : log.printf(" deactived optimization of %zu coefficients\n",ndeactived);
517 : }
518 : }
519 : }
520 :
521 70 : std::vector<std::string> mask_fnames_out;
522 140 : parseFilenames("OUTPUT_MASK_FILE",mask_fnames_out);
523 :
524 143 : for(unsigned int i=0; i<mask_fnames_out.size(); i++) {
525 1 : if(mask_fnames_in.size()>0) {
526 1 : plumed_massert(mask_fnames_out[i]!=mask_fnames_in[i],"MASK_FILE and OUTPUT_MASK_FILE cannot be the same");
527 : }
528 2 : OFile maskOFile;
529 1 : maskOFile.link(*this);
530 1 : maskOFile.enforceBackup();
531 1 : if(use_mwalkers_mpi_ && mwalkers_mpi_single_files_) {
532 0 : unsigned int r=0;
533 0 : if(comm.Get_rank()==0) {r=multi_sim_comm.Get_rank();}
534 0 : comm.Bcast(r,0);
535 0 : if(r>0) {mask_fnames_out[i]="/dev/null";}
536 0 : maskOFile.enforceSuffix("");
537 : }
538 1 : maskOFile.open(mask_fnames_out[i]);
539 1 : coeffs_mask_pntrs_[i]->writeToFile(maskOFile,true);
540 1 : maskOFile.close();
541 : }
542 : }
543 :
544 142 : if(getRestart() && ustride_targetdist_>0) {
545 24 : for(unsigned int i=0; i<nbiases_; i++) {
546 16 : if(dynamic_targetdists_[i]) {
547 8 : bias_pntrs_[i]->restartTargetDistributions();
548 : }
549 : }
550 : }
551 :
552 :
553 71 : std::vector<std::string> targetdist_averages_fnames;
554 142 : if(keywords.exists("TARGETDIST_AVERAGES_FILE")) {
555 213 : parseFilenames("TARGETDIST_AVERAGES_FILE",targetdist_averages_fnames,"targetdist-averages.data");
556 142 : parse("TARGETDIST_AVERAGES_OUTPUT",targetdist_averages_wstride_);
557 :
558 71 : if(coeffs_fnames.size()>0) {
559 373 : for(unsigned int i=0; i<targetdist_averages_fnames.size(); i++) {
560 77 : plumed_massert(targetdist_averages_fnames[i]!=coeffs_fnames[i],"COEFFS_FILE and TARGETDIST_AVERAGES_FILE cannot be the same");
561 : }
562 : }
563 71 : if(gradient_fnames.size()>0) {
564 338 : for(unsigned int i=0; i<targetdist_averages_fnames.size(); i++) {
565 70 : plumed_massert(targetdist_averages_fnames[i]!=gradient_fnames[i],"GRADIENT_FILE and TARGETDIST_AVERAGES_FILE cannot be the same");
566 : }
567 : }
568 71 : if(hessian_fnames.size()>0) {
569 317 : for(unsigned int i=0; i<targetdist_averages_fnames.size(); i++) {
570 65 : plumed_massert(targetdist_averages_fnames[i]!=hessian_fnames[i],"HESSIAN_FILE and TARGETDIST_AVERAGES_FILE cannot be the same");
571 : }
572 : }
573 71 : setupOFiles(targetdist_averages_fnames,targetdist_averagesOFiles_,mw_single_files);
574 142 : parse("TARGETDIST_AVERAGES_FMT",targetdist_averages_output_fmt_);
575 71 : if(targetdist_averages_output_fmt_.size()>0) {
576 225 : for(unsigned int i=0; i<ncoeffssets_; i++) {
577 154 : targetdist_averages_pntrs_[i]->setOutputFmt(targetdist_averages_output_fmt_);
578 : }
579 : }
580 :
581 373 : for(unsigned int i=0; i<targetdist_averagesOFiles_.size(); i++) {
582 154 : targetdist_averages_pntrs_[i]->writeToFile(*targetdist_averagesOFiles_[i]);
583 : }
584 :
585 71 : if(targetdist_averages_wstride_==0) {
586 208 : for(unsigned int i=0; i<targetdist_averagesOFiles_.size(); i++) {
587 44 : targetdist_averagesOFiles_[i]->close();
588 44 : delete targetdist_averagesOFiles_[i];
589 : }
590 : targetdist_averagesOFiles_.clear();
591 : }
592 :
593 71 : if(targetdist_averages_fnames.size()>0 && targetdist_averages_wstride_ > 0) {
594 33 : if(ncoeffssets_==1) {
595 132 : log.printf(" Target distribution averages will be written out to file %s every %u iterations\n",targetdist_averagesOFiles_[0]->getPath().c_str(),targetdist_averages_wstride_);
596 : }
597 : else {
598 0 : log.printf(" Target distribution averages will be written out to the following files every %u iterations:\n",targetdist_averages_wstride_);
599 0 : for(unsigned int i=0; i<targetdist_averages_fnames.size(); i++) {
600 0 : log.printf(" coefficient set %u: %s\n",i,targetdist_averagesOFiles_[i]->getPath().c_str());
601 : }
602 : }
603 : }
604 : }
605 :
606 :
607 142 : if(keywords.exists("BIAS_OUTPUT")) {
608 142 : parse("BIAS_OUTPUT",bias_output_stride_);
609 71 : if(bias_output_stride_>0) {
610 70 : bias_output_active_=true;
611 222 : for(unsigned int i=0; i<nbiases_; i++) {
612 152 : bias_pntrs_[i]->enableBiasFileOutput();
613 76 : bias_pntrs_[i]->setupBiasFileOutput();
614 76 : bias_pntrs_[i]->writeBiasToFile();
615 : }
616 : }
617 : else {
618 1 : bias_output_active_=false;
619 1 : bias_output_stride_=1000;
620 : }
621 : }
622 :
623 142 : if(keywords.exists("FES_OUTPUT")) {
624 142 : parse("FES_OUTPUT",fes_output_stride_);
625 71 : if(fes_output_stride_>0) {
626 70 : fes_output_active_=true;
627 222 : for(unsigned int i=0; i<nbiases_; i++) {
628 152 : bias_pntrs_[i]->enableFesFileOutput();
629 76 : bias_pntrs_[i]->setupFesFileOutput();
630 76 : bias_pntrs_[i]->writeFesToFile();
631 : }
632 : }
633 : else {
634 1 : fes_output_active_=false;
635 1 : fes_output_stride_=1000;
636 : }
637 : }
638 :
639 142 : if(keywords.exists("FES_PROJ_OUTPUT")) {
640 142 : parse("FES_PROJ_OUTPUT",fesproj_output_stride_);
641 71 : if(fesproj_output_stride_>0) {
642 16 : fesproj_output_active_=true;
643 48 : for(unsigned int i=0; i<nbiases_; i++) {
644 32 : bias_pntrs_[i]->enableFesProjFileOutput();
645 16 : bias_pntrs_[i]->setupFesProjFileOutput();
646 16 : bias_pntrs_[i]->writeFesProjToFile();
647 : }
648 : }
649 : else {
650 55 : fesproj_output_active_=false;
651 55 : fesproj_output_stride_=1000;
652 : }
653 : }
654 :
655 225 : for(unsigned int i=0; i<nbiases_; i++) {
656 242 : if(!dynamic_targetdists_[i] && bias_pntrs_[i]->isStaticTargetDistFileOutputActive()) {
657 4 : bias_pntrs_[i]->setupTargetDistFileOutput();
658 4 : bias_pntrs_[i]->writeTargetDistToFile();
659 4 : bias_pntrs_[i]->setupTargetDistProjFileOutput();
660 4 : bias_pntrs_[i]->writeTargetDistProjToFile();
661 : }
662 : }
663 :
664 142 : if(keywords.exists("TARGETDIST_OUTPUT")) {
665 140 : parse("TARGETDIST_OUTPUT",targetdist_output_stride_);
666 70 : if(targetdist_output_stride_>0) {
667 32 : if(ustride_targetdist_==0) {
668 0 : plumed_merror("it doesn't make sense to use the TARGETDIST_OUTPUT keyword if you don't have a target distribution that needs to be updated");
669 : }
670 32 : if(targetdist_output_stride_%ustride_targetdist_!=0) {
671 0 : plumed_merror("the value given in TARGETDIST_OUTPUT doesn't make sense, it should be multiple of TARGETDIST_STRIDE");
672 : }
673 32 : if(targetdist_output_stride_%coeffs_wstride_!=0) {
674 0 : plumed_merror("the value given in TARGETDIST_OUTPUT doesn't make sense, it should be multiple of COEFFS_OUTPUT");
675 : }
676 :
677 32 : targetdist_output_active_=true;
678 96 : for(unsigned int i=0; i<nbiases_; i++) {
679 64 : if(dynamic_targetdists_[i]) {
680 32 : bias_pntrs_[i]->enableDynamicTargetDistFileOutput();
681 32 : bias_pntrs_[i]->setupTargetDistFileOutput();
682 32 : bias_pntrs_[i]->writeTargetDistToFile();
683 : }
684 : }
685 : }
686 : else {
687 38 : targetdist_output_active_=false;
688 38 : targetdist_output_stride_=1000;
689 : }
690 : }
691 :
692 142 : if(keywords.exists("TARGETDIST_PROJ_OUTPUT")) {
693 140 : parse("TARGETDIST_PROJ_OUTPUT",targetdist_proj_output_stride_);
694 70 : if(targetdist_proj_output_stride_>0) {
695 3 : if(ustride_targetdist_==0) {
696 0 : plumed_merror("it doesn't make sense to use the TARGETDIST_PROJ_OUTPUT keyword if you don't have a target distribution that needs to be updated");
697 : }
698 3 : if(targetdist_proj_output_stride_%ustride_targetdist_!=0) {
699 0 : plumed_merror("the value given in TARGETDIST_PROJ_OUTPUT doesn't make sense, it should be multiple of TARGETDIST_STRIDE");
700 : }
701 :
702 3 : targetdist_proj_output_active_=true;
703 9 : for(unsigned int i=0; i<nbiases_; i++) {
704 6 : if(dynamic_targetdists_[i]) {
705 3 : bias_pntrs_[i]->enableDynamicTargetDistFileOutput();
706 3 : bias_pntrs_[i]->setupTargetDistProjFileOutput();
707 3 : bias_pntrs_[i]->writeTargetDistProjToFile();
708 : }
709 : }
710 : }
711 : else {
712 67 : targetdist_proj_output_active_=false;
713 67 : targetdist_proj_output_stride_=1000;
714 : }
715 : }
716 :
717 71 : if(ncoeffssets_==1) {
718 68 : log.printf(" Output Components:\n");
719 68 : log.printf(" ");
720 68 : if(monitor_instantaneous_gradient_) {
721 6 : addComponent("gradrms"); componentIsNotPeriodic("gradrms");
722 2 : log.printf(" ");
723 6 : addComponent("gradmax"); componentIsNotPeriodic("gradmax");
724 : }
725 68 : if(aver_gradient_pntrs_.size()>0) {
726 2 : log.printf(" ");
727 6 : addComponent("avergradrms"); componentIsNotPeriodic("avergradrms");
728 2 : log.printf(" ");
729 6 : addComponent("avergradmax"); componentIsNotPeriodic("avergradmax");
730 : }
731 68 : if(!fixed_stepsize_) {
732 0 : log.printf(" ");
733 0 : addComponent("stepsize"); componentIsNotPeriodic("stepsize");
734 0 : getPntrToComponent("stepsize")->set( getCurrentStepSize(0) );
735 : }
736 : }
737 : else {
738 21 : for(unsigned int i=0; i<ncoeffssets_; i++) {
739 9 : log.printf(" Output Components for coefficent set %u:\n",i);
740 27 : std::string is=""; Tools::convert(i,is); is = "_" + coeffssetid_prefix_ + is;
741 9 : log.printf(" ");
742 9 : if(monitor_instantaneous_gradient_) {
743 0 : addComponent("gradrms"+is); componentIsNotPeriodic("gradrms"+is);
744 0 : log.printf(" ");
745 0 : addComponent("gradmax"+is); componentIsNotPeriodic("gradmax"+is);
746 : }
747 9 : if(aver_gradient_pntrs_.size()>0) {
748 0 : log.printf(" ");
749 0 : addComponent("avergradrms"+is); componentIsNotPeriodic("avergradrms"+is);
750 0 : log.printf(" ");
751 0 : addComponent("avergradmax"+is); componentIsNotPeriodic("avergradmax"+is);
752 : }
753 9 : if(!fixed_stepsize_) {
754 0 : log.printf(" ");
755 0 : addComponent("stepsize"+is); componentIsNotPeriodic("stepsize"+is);
756 0 : getPntrToComponent("stepsize"+is)->set( getCurrentStepSize(i) );
757 : }
758 : }
759 : }
760 :
761 71 : }
762 :
763 :
764 142 : Optimizer::~Optimizer() {
765 : //
766 225 : for(unsigned int i=0; i<ncoeffssets_; i++) {
767 153 : if(coeffsOFiles_.size()>0 && getIterationCounter()%coeffs_wstride_!=0) {
768 8 : coeffs_pntrs_[i]->writeToFile(*coeffsOFiles_[i],aux_coeffs_pntrs_[i],false);
769 : }
770 77 : if(targetdist_averagesOFiles_.size()>0 && iter_counter%targetdist_averages_wstride_!=0) {
771 3 : targetdist_averages_pntrs_[i]->writeToFile(*targetdist_averagesOFiles_[i]);
772 : }
773 : }
774 : //
775 71 : if(!isTargetDistOutputActive()) {
776 129 : for(unsigned int i=0; i<nbiases_; i++) {
777 90 : if(dynamic_targetdists_[i]) {
778 1 : bias_pntrs_[i]->enableDynamicTargetDistFileOutput();
779 1 : bias_pntrs_[i]->setupTargetDistFileOutput();
780 1 : bias_pntrs_[i]->writeTargetDistToFile();
781 : }
782 : }
783 : }
784 : //
785 141 : if(isBiasOutputActive() && getIterationCounter()%getBiasOutputStride()!=0) {
786 1 : writeBiasOutputFiles();
787 : }
788 141 : if(isFesOutputActive() && getIterationCounter()%getFesOutputStride()!=0) {
789 1 : writeFesOutputFiles();
790 : }
791 87 : if(isFesProjOutputActive() && getIterationCounter()%getFesProjOutputStride()!=0) {
792 1 : writeFesProjOutputFiles();
793 : }
794 103 : if(isTargetDistOutputActive() && getIterationCounter()%getTargetDistOutputStride()!=0) {
795 2 : writeTargetDistOutputFiles();
796 : }
797 74 : if(isTargetDistProjOutputActive() && getIterationCounter()%getTargetDistProjOutputStride()!=0) {
798 1 : writeTargetDistProjOutputFiles();
799 : }
800 : //
801 373 : for(unsigned int i=0; i<aux_coeffs_pntrs_.size(); i++) {
802 77 : delete aux_coeffs_pntrs_[i];
803 : }
804 : aux_coeffs_pntrs_.clear();
805 : //
806 148 : for(unsigned int i=0; i<aver_gradient_pntrs_.size(); i++) {
807 2 : delete aver_gradient_pntrs_[i];
808 : }
809 : aver_gradient_pntrs_.clear();
810 : //
811 370 : for(unsigned int i=0; i<coeffs_mask_pntrs_.size(); i++) {
812 76 : delete coeffs_mask_pntrs_[i];
813 : }
814 : coeffs_mask_pntrs_.clear();
815 : //
816 370 : for(unsigned int i=0; i<coeffsOFiles_.size(); i++) {
817 76 : coeffsOFiles_[i]->close();
818 76 : delete coeffsOFiles_[i];
819 : }
820 : coeffsOFiles_.clear();
821 352 : for(unsigned int i=0; i<gradientOFiles_.size(); i++) {
822 70 : gradientOFiles_[i]->close();
823 70 : delete gradientOFiles_[i];
824 : }
825 : gradientOFiles_.clear();
826 337 : for(unsigned int i=0; i<hessianOFiles_.size(); i++) {
827 65 : hessianOFiles_[i]->close();
828 65 : delete hessianOFiles_[i];
829 : }
830 : hessianOFiles_.clear();
831 241 : for(unsigned int i=0; i<targetdist_averagesOFiles_.size(); i++) {
832 33 : targetdist_averagesOFiles_[i]->close();
833 33 : delete targetdist_averagesOFiles_[i];
834 : }
835 : targetdist_averagesOFiles_.clear();
836 71 : }
837 :
838 :
839 73 : void Optimizer::registerKeywords( Keywords& keys ) {
840 73 : Action::registerKeywords(keys);
841 73 : ActionPilot::registerKeywords(keys);
842 73 : ActionWithValue::registerKeywords(keys);
843 : //
844 146 : keys.remove("NUMERICAL_DERIVATIVES");
845 : // Default always active keywords
846 292 : keys.add("compulsory","BIAS","the label of the VES bias to be optimized");
847 292 : keys.add("compulsory","STRIDE","the frequency of updating the coefficients given in the number of MD steps.");
848 365 : keys.add("compulsory","COEFFS_FILE","coeffs.data","the name of output file for the coefficients");
849 365 : keys.add("compulsory","COEFFS_OUTPUT","100","how often the coefficients should be written to file. This parameter is given as the number of iterations.");
850 292 : keys.add("optional","COEFFS_FMT","specify format for coefficient file(s) (useful for decrease the number of digits in regtests)");
851 292 : keys.add("optional","COEFFS_SET_ID_PREFIX","suffix to add to the filename given in FILE to identfy the bias, should only be given if a single filename is given in FILE when optimizing multiple biases.");
852 : //
853 292 : keys.add("optional","INITIAL_COEFFS","the name(s) of file(s) with the initial coefficents");
854 : // Hidden keywords to output the gradient to a file.
855 292 : keys.add("hidden","GRADIENT_FILE","the name of output file for the gradient");
856 292 : keys.add("hidden","GRADIENT_OUTPUT","how often the gradient should be written to file. This parameter is given as the number of bias iterations. It is by default 100 if GRADIENT_FILE is specficed");
857 292 : keys.add("hidden","GRADIENT_FMT","specify format for gradient file(s) (useful for decrease the number of digits in regtests)");
858 : // Either use a fixed stepsize (useFixedStepSizeKeywords) or changing stepsize (useDynamicsStepSizeKeywords)
859 292 : keys.reserve("compulsory","STEPSIZE","the step size used for the optimization");
860 292 : keys.reserve("compulsory","INITIAL_STEPSIZE","the initial step size used for the optimization");
861 : // Keywords related to the Hessian, actived with the useHessianKeywords function
862 219 : keys.reserveFlag("FULL_HESSIAN",false,"if the full Hessian matrix should be used for the optimization, otherwise only the diagonal part of the Hessian is used");
863 292 : keys.reserve("hidden","HESSIAN_FILE","the name of output file for the Hessian");
864 292 : keys.reserve("hidden","HESSIAN_OUTPUT","how often the Hessian should be written to file. This parameter is given as the number of bias iterations. It is by default 100 if HESSIAN_FILE is specficed");
865 292 : keys.reserve("hidden","HESSIAN_FMT","specify format for hessian file(s) (useful for decrease the number of digits in regtests)");
866 : // Keywords related to the multiple walkers, actived with the useMultipleWalkersKeywords function
867 219 : keys.reserveFlag("MULTIPLE_WALKERS",false,"if optimization is to be performed using multiple walkers connected via MPI");
868 : // Keywords related to the mask file, actived with the useMaskKeywords function
869 292 : keys.reserve("optional","MASK_FILE","read in a mask file which allows one to employ different step sizes for different coefficents and/or deactive the optimization of certain coefficients (by putting values of 0.0). One can write out the resulting mask by using the OUTPUT_MASK_FILE keyword.");
870 292 : keys.reserve("optional","OUTPUT_MASK_FILE","Name of the file to write out the mask resulting from using the MASK_FILE keyword. Can also be used to generate a template mask file.");
871 : //
872 219 : keys.reserveFlag("START_OPTIMIZATION_AFRESH",false,"if the iterations should be started afresh when a restart has been triggered by the RESTART keyword or the MD code.");
873 : //
874 219 : keys.addFlag("MONITOR_INSTANTANEOUS_GRADIENT",false,"if quantities related to the instantaneous gradient should be outputted.");
875 : //
876 219 : keys.reserveFlag("MONITOR_AVERAGE_GRADIENT",false,"if the averaged gradient should be monitored and quantities related to it should be outputted.");
877 292 : keys.reserve("optional","MONITOR_AVERAGES_GRADIENT_EXP_DECAY","use an exponentially decaying averaging with a given time constant when monitoring the averaged gradient");
878 : //
879 292 : keys.reserve("optional","TARGETDIST_STRIDE","stride for updating a target distribution that is iteratively updated during the optimization. Note that the value is given in terms of coefficent iterations.");
880 292 : keys.reserve("optional","TARGETDIST_OUTPUT","how often the dynamic target distribution(s) should be written out to file. Note that the value is given in terms of coefficent iterations.");
881 292 : keys.reserve("optional","TARGETDIST_PROJ_OUTPUT","how often the projections of the dynamic target distribution(s) should be written out to file. Note that the value is given in terms of coefficent iterations.");
882 : //
883 292 : keys.add("optional","TARGETDIST_AVERAGES_FILE","the name of output file for the target distribution averages. By default it is targetdist-averages.data.");
884 292 : keys.add("optional","TARGETDIST_AVERAGES_OUTPUT","how often the target distribution averages should be written out to file. Note that the value is given in terms of coefficent iterations. If no value is given are the averages only written at the begining of the optimization");
885 292 : keys.add("hidden","TARGETDIST_AVERAGES_FMT","specify format for target distribution averages file(s) (useful for decrease the number of digits in regtests)");
886 : //
887 292 : keys.add("optional","BIAS_OUTPUT","how often the bias(es) should be written out to file. Note that the value is given in terms of coefficent iterations.");
888 292 : keys.add("optional","FES_OUTPUT","how often the FES(s) should be written out to file. Note that the value is given in terms of coefficent iterations.");
889 292 : keys.add("optional","FES_PROJ_OUTPUT","how often the projections of the FES(s) should be written out to file. Note that the value is given in terms of coefficent iterations.");
890 : //
891 292 : keys.reserve("optional","REWEIGHT_FACTOR_STRIDE","stride for updating the reweighting factor c(t). Note that the value is given in terms of coefficent iterations.");
892 : //
893 146 : keys.use("RESTART");
894 : //
895 146 : keys.use("UPDATE_FROM");
896 146 : keys.use("UPDATE_UNTIL");
897 : // Components that are always active
898 292 : keys.addOutputComponent("gradrms","MONITOR_INSTANTANEOUS_GRADIENT","the root mean square value of the coefficent gradient. For multiple biases this component is labeled using the number of the bias as gradrms-#.");
899 292 : keys.addOutputComponent("gradmax","MONITOR_INSTANTANEOUS_GRADIENT","the largest absolute value of the coefficent gradient. For multiple biases this component is labeled using the number of the bias as gradmax-#.");
900 73 : ActionWithValue::useCustomisableComponents(keys);
901 : // keys.addOutputComponent("gradmaxidx","default","the index of the maximum absolute value of the gradient");
902 :
903 73 : }
904 :
905 :
906 73 : void Optimizer::useHessianKeywords(Keywords& keys) {
907 : // keys.use("FULL_HESSIAN");
908 146 : keys.use("HESSIAN_FILE");
909 146 : keys.use("HESSIAN_OUTPUT");
910 146 : keys.use("HESSIAN_FMT");
911 73 : }
912 :
913 :
914 73 : void Optimizer::useMultipleWalkersKeywords(Keywords& keys) {
915 146 : keys.use("MULTIPLE_WALKERS");
916 73 : }
917 :
918 :
919 71 : void Optimizer::useFixedStepSizeKeywords(Keywords& keys) {
920 142 : keys.use("STEPSIZE");
921 71 : }
922 :
923 :
924 0 : void Optimizer::useDynamicStepSizeKeywords(Keywords& keys) {
925 0 : keys.use("INITIAL_STEPSIZE");
926 0 : keys.addOutputComponent("stepsize","default","the current value of step size used to update the coefficients. For multiple biases this component is labeled using the number of the bias as stepsize-#.");
927 0 : }
928 :
929 :
930 71 : void Optimizer::useMaskKeywords(Keywords& keys) {
931 142 : keys.use("MASK_FILE");
932 142 : keys.use("OUTPUT_MASK_FILE");
933 71 : }
934 :
935 :
936 71 : void Optimizer::useRestartKeywords(Keywords& keys) {
937 142 : keys.use("START_OPTIMIZATION_AFRESH");
938 71 : }
939 :
940 :
941 73 : void Optimizer::useMonitorAverageGradientKeywords(Keywords& keys) {
942 146 : keys.use("MONITOR_AVERAGE_GRADIENT");
943 146 : keys.use("MONITOR_AVERAGES_GRADIENT_EXP_DECAY");
944 292 : keys.addOutputComponent("avergradrms","MONITOR_AVERAGE_GRADIENT","the root mean square value of the averaged coefficent gradient. For multiple biases this component is labeled using the number of the bias as gradrms-#.");
945 292 : keys.addOutputComponent("avergradmax","MONITOR_AVERAGE_GRADIENT","the largest absolute value of the averaged coefficent gradient. For multiple biases this component is labeled using the number of the bias as gradmax-#.");
946 73 : }
947 :
948 :
949 71 : void Optimizer::useDynamicTargetDistributionKeywords(Keywords& keys) {
950 142 : keys.use("TARGETDIST_STRIDE");
951 142 : keys.use("TARGETDIST_OUTPUT");
952 142 : keys.use("TARGETDIST_PROJ_OUTPUT");
953 71 : }
954 :
955 :
956 0 : void Optimizer::useReweightFactorKeywords(Keywords& keys) {
957 0 : keys.use("REWEIGHT_FACTOR_STRIDE");
958 0 : }
959 :
960 :
961 71 : void Optimizer::turnOnHessian() {
962 71 : plumed_massert(hessian_pntrs_.size()==0,"turnOnHessian() should only be run during initialization");
963 71 : use_hessian_=true;
964 71 : hessian_pntrs_.clear();
965 225 : for(unsigned int i=0; i<nbiases_; i++) {
966 154 : std::vector<CoeffsMatrix*> pntrs_hessian = enableHessian(bias_pntrs_[i],diagonal_hessian_);
967 385 : for(unsigned int k=0; k<pntrs_hessian.size(); k++) {
968 77 : pntrs_hessian[k]->turnOnIterationCounter();
969 77 : pntrs_hessian[k]->setIterationCounterAndTime(getIterationCounter(),getTime());
970 77 : hessian_pntrs_.push_back(pntrs_hessian[k]);
971 : }
972 : }
973 71 : plumed_massert(hessian_pntrs_.size()==ncoeffssets_,"problems in linking Hessians");
974 71 : if(diagonal_hessian_) {
975 71 : log.printf(" Optimization performed using diagonal Hessian matrix\n");
976 : }
977 : else {
978 0 : log.printf(" Optimization performed using full Hessian matrix\n");
979 : }
980 : //
981 71 : if(hessian_output_fmt_.size()>0) {
982 191 : for(unsigned int i=0; i<ncoeffssets_; i++) {
983 130 : hessian_pntrs_[i]->setOutputFmt(hessian_output_fmt_);
984 : }
985 : }
986 :
987 71 : }
988 :
989 :
990 0 : void Optimizer::turnOffHessian() {
991 0 : use_hessian_=false;
992 0 : for(unsigned int i=0; i<nbiases_; i++) {
993 0 : bias_pntrs_[i]->disableHessian();
994 : }
995 : hessian_pntrs_.clear();
996 0 : for(unsigned int i=0; i<hessianOFiles_.size(); i++) {
997 0 : hessianOFiles_[i]->close();
998 0 : delete hessianOFiles_[i];
999 : }
1000 : hessianOFiles_.clear();
1001 0 : }
1002 :
1003 :
1004 77 : std::vector<CoeffsMatrix*> Optimizer::enableHessian(VesBias* bias_pntr_in, const bool diagonal_hessian) {
1005 77 : plumed_massert(use_hessian_,"the Hessian should not be used");
1006 77 : bias_pntr_in->enableHessian(diagonal_hessian);
1007 : std::vector<CoeffsMatrix*> hessian_pntrs_out = bias_pntr_in->getHessianPntrs();
1008 385 : for(unsigned int k=0; k<hessian_pntrs_out.size(); k++) {
1009 77 : plumed_massert(hessian_pntrs_out[k] != NULL,"Hessian is needed but not linked correctly");
1010 : }
1011 77 : return hessian_pntrs_out;
1012 : }
1013 :
1014 :
1015 : // CoeffsMatrix* Optimizer::switchToDiagonalHessian(VesBias* bias_pntr_in) {
1016 : // plumed_massert(use_hessian_,"it does not make sense to switch to diagonal Hessian if it Hessian is not used");
1017 : // diagonal_hessian_=true;
1018 : // bias_pntr_in->enableHessian(diagonal_hessian_);
1019 : // CoeffsMatrix* hessian_pntr_out = bias_pntr_in->getHessianPntr();
1020 : // plumed_massert(hessian_pntr_out != NULL,"Hessian is needed but not linked correctly");
1021 : // //
1022 : // log.printf(" %s (with label %s): switching to a diagonal Hessian for VES bias %s (with label %s) at time %f\n",getName().c_str(),getLabel().c_str(),bias_pntr_in->getName().c_str(),bias_pntr_in->getLabel().c_str(),getTime());
1023 : // return hessian_pntr_out;
1024 : // }
1025 :
1026 :
1027 : // CoeffsMatrix* Optimizer::switchToFullHessian(VesBias* bias_pntr_in) {
1028 : // plumed_massert(use_hessian_,"it does not make sense to switch to diagonal Hessian if it Hessian is not used");
1029 : // diagonal_hessian_=false;
1030 : // bias_pntr_in->enableHessian(diagonal_hessian_);
1031 : // CoeffsMatrix* hessian_pntr_out = bias_pntr_in->getHessianPntr();
1032 : // plumed_massert(hessian_pntr_out != NULL,"Hessian is needed but not linked correctly");
1033 : // //
1034 : // log.printf(" %s (with label %s): switching to a diagonal Hessian for VES bias %s (with label %s) at time %f\n",getName().c_str(),getLabel().c_str(),bias_pntr_in->getName().c_str(),bias_pntr_in->getLabel().c_str(),getTime());
1035 : // return hessian_pntr_out;
1036 : // }
1037 :
1038 :
1039 821 : void Optimizer::update() {
1040 821 : if(onStep() && !isFirstStep) {
1041 2250 : for(unsigned int i=0; i<nbiases_; i++) {
1042 1540 : bias_pntrs_[i]->updateGradientAndHessian(use_mwalkers_mpi_);
1043 : }
1044 2250 : for(unsigned int i=0; i<ncoeffssets_; i++) {
1045 2310 : if(gradient_pntrs_[i]->isActive()) {coeffsUpdate(i);}
1046 : else {
1047 360 : std::string msg = "iteration " + getIterationCounterStr(+1) +
1048 270 : " for " + bias_pntrs_[i]->getLabel() +
1049 90 : " - the coefficients are not updated as CV values are outside the bias intervals";
1050 90 : warning(msg);
1051 : }
1052 :
1053 : // +1 as this is done before increaseIterationCounter() is used
1054 770 : unsigned int curr_iter = getIterationCounter()+1;
1055 770 : double curr_time = getTime();
1056 770 : coeffs_pntrs_[i]->setIterationCounterAndTime(curr_iter,curr_time);
1057 770 : aux_coeffs_pntrs_[i]->setIterationCounterAndTime(curr_iter,curr_time);
1058 770 : gradient_pntrs_[i]->setIterationCounterAndTime(curr_iter,curr_time);
1059 770 : targetdist_averages_pntrs_[i]->setIterationCounterAndTime(curr_iter,curr_time);
1060 770 : if(use_hessian_) {
1061 770 : hessian_pntrs_[i]->setIterationCounterAndTime(curr_iter,curr_time);
1062 : }
1063 770 : if(aver_gradient_pntrs_.size()>0) {
1064 20 : aver_gradient_pntrs_[i]->setIterationCounterAndTime(curr_iter,curr_time);
1065 40 : aver_gradient_pntrs_[i]->addToAverage(*gradient_pntrs_[i]);
1066 : }
1067 : }
1068 : increaseIterationCounter();
1069 :
1070 710 : if(ustride_targetdist_>0 && getIterationCounter()%ustride_targetdist_==0) {
1071 990 : for(unsigned int i=0; i<nbiases_; i++) {
1072 660 : if(dynamic_targetdists_[i]) {
1073 330 : bias_pntrs_[i]->updateTargetDistributions();
1074 : }
1075 : }
1076 : }
1077 710 : if(ustride_reweightfactor_>0 && getIterationCounter()%ustride_reweightfactor_==0) {
1078 0 : for(unsigned int i=0; i<nbiases_; i++) {
1079 0 : bias_pntrs_[i]->updateReweightFactor();
1080 : }
1081 : }
1082 :
1083 710 : updateOutputComponents();
1084 2250 : for(unsigned int i=0; i<ncoeffssets_; i++) {
1085 770 : writeOutputFiles(i);
1086 : }
1087 : //
1088 1410 : if(isBiasOutputActive() && getIterationCounter()%getBiasOutputStride()==0) {
1089 70 : writeBiasOutputFiles();
1090 : }
1091 1410 : if(isFesOutputActive() && getIterationCounter()%getFesOutputStride()==0) {
1092 70 : writeFesOutputFiles();
1093 : }
1094 870 : if(isFesProjOutputActive() && getIterationCounter()%getFesProjOutputStride()==0) {
1095 16 : writeFesProjOutputFiles();
1096 : }
1097 1030 : if(isTargetDistOutputActive() && getIterationCounter()%getTargetDistOutputStride()==0) {
1098 31 : writeTargetDistOutputFiles();
1099 : }
1100 740 : if(isTargetDistProjOutputActive() && getIterationCounter()%getTargetDistProjOutputStride()==0) {
1101 3 : writeTargetDistProjOutputFiles();
1102 : }
1103 : }
1104 : else {
1105 111 : isFirstStep=false;
1106 : }
1107 821 : }
1108 :
1109 :
1110 710 : void Optimizer::updateOutputComponents() {
1111 710 : if(ncoeffssets_==1) {
1112 680 : if(!fixed_stepsize_) {
1113 0 : getPntrToComponent("stepsize")->set( getCurrentStepSize(0) );
1114 : }
1115 680 : if(monitor_instantaneous_gradient_) {
1116 60 : getPntrToComponent("gradrms")->set( gradient_pntrs_[0]->getRMS() );
1117 20 : size_t gradient_maxabs_idx=0;
1118 60 : getPntrToComponent("gradmax")->set( gradient_pntrs_[0]->getMaxAbsValue(gradient_maxabs_idx) );
1119 : }
1120 680 : if(aver_gradient_pntrs_.size()>0) {
1121 60 : getPntrToComponent("avergradrms")->set( aver_gradient_pntrs_[0]->getRMS() );
1122 20 : size_t avergradient_maxabs_idx=0;
1123 60 : getPntrToComponent("avergradmax")->set( aver_gradient_pntrs_[0]->getMaxAbsValue(avergradient_maxabs_idx) );
1124 : }
1125 : }
1126 : else {
1127 210 : for(unsigned int i=0; i<ncoeffssets_; i++) {
1128 270 : std::string is=""; Tools::convert(i,is); is = "_" + coeffssetid_prefix_ + is;
1129 90 : if(!fixed_stepsize_) {
1130 0 : getPntrToComponent("stepsize"+is)->set( getCurrentStepSize(i) );
1131 : }
1132 90 : if(monitor_instantaneous_gradient_) {
1133 0 : getPntrToComponent("gradrms"+is)->set( gradient_pntrs_[i]->getRMS() );
1134 0 : size_t gradient_maxabs_idx=0;
1135 0 : getPntrToComponent("gradmax"+is)->set( gradient_pntrs_[i]->getMaxAbsValue(gradient_maxabs_idx) );
1136 : }
1137 90 : if(aver_gradient_pntrs_.size()>0) {
1138 0 : getPntrToComponent("avergradrms"+is)->set( aver_gradient_pntrs_[0]->getRMS() );
1139 0 : size_t avergradient_maxabs_idx=0;
1140 0 : getPntrToComponent("avergradmax"+is)->set( aver_gradient_pntrs_[0]->getMaxAbsValue(avergradient_maxabs_idx) );
1141 : }
1142 : }
1143 : }
1144 710 : }
1145 :
1146 :
1147 1 : void Optimizer::turnOffCoeffsOutputFiles() {
1148 5 : for(unsigned int i=0; i<coeffsOFiles_.size(); i++) {
1149 1 : coeffsOFiles_[i]->close();
1150 1 : delete coeffsOFiles_[i];
1151 : }
1152 : coeffsOFiles_.clear();
1153 1 : }
1154 :
1155 :
1156 770 : void Optimizer::writeOutputFiles(const unsigned int coeffs_id) {
1157 770 : if(coeffsOFiles_.size()>0 && iter_counter%coeffs_wstride_==0) {
1158 2972 : coeffs_pntrs_[coeffs_id]->writeToFile(*coeffsOFiles_[coeffs_id],aux_coeffs_pntrs_[coeffs_id],false);
1159 : }
1160 770 : if(gradientOFiles_.size()>0 && iter_counter%gradient_wstride_==0) {
1161 700 : if(aver_gradient_pntrs_.size()==0) {
1162 2040 : gradient_pntrs_[coeffs_id]->writeToFile(*gradientOFiles_[coeffs_id],false);
1163 : }
1164 : else {
1165 80 : gradient_pntrs_[coeffs_id]->writeToFile(*gradientOFiles_[coeffs_id],aver_gradient_pntrs_[coeffs_id],false);
1166 : }
1167 : }
1168 770 : if(hessianOFiles_.size()>0 && iter_counter%hessian_wstride_==0) {
1169 1950 : hessian_pntrs_[coeffs_id]->writeToFile(*hessianOFiles_[coeffs_id]);
1170 : }
1171 770 : if(targetdist_averagesOFiles_.size()>0 && iter_counter%targetdist_averages_wstride_==0) {
1172 969 : targetdist_averages_pntrs_[coeffs_id]->writeToFile(*targetdist_averagesOFiles_[coeffs_id]);
1173 : }
1174 770 : }
1175 :
1176 :
1177 354 : void Optimizer::setupOFiles(std::vector<std::string>& fnames, std::vector<OFile*>& OFiles, const bool multi_sim_single_files) {
1178 354 : plumed_assert(ncoeffssets_>0);
1179 708 : OFiles.resize(fnames.size(),NULL);
1180 1593 : for(unsigned int i=0; i<fnames.size(); i++) {
1181 590 : OFiles[i] = new OFile();
1182 295 : OFiles[i]->link(*this);
1183 295 : if(multi_sim_single_files) {
1184 48 : unsigned int r=0;
1185 48 : if(comm.Get_rank()==0) {r=multi_sim_comm.Get_rank();}
1186 48 : comm.Bcast(r,0);
1187 48 : if(r>0) {fnames[i]="/dev/null";}
1188 96 : OFiles[i]->enforceSuffix("");
1189 : }
1190 590 : OFiles[i]->open(fnames[i]);
1191 295 : OFiles[i]->setHeavyFlush();
1192 : }
1193 354 : }
1194 :
1195 :
1196 14 : void Optimizer::readCoeffsFromFiles(const std::vector<std::string>& fnames, const bool read_aux_coeffs) {
1197 14 : plumed_assert(ncoeffssets_>0);
1198 14 : plumed_assert(fnames.size()==ncoeffssets_);
1199 14 : if(ncoeffssets_==1) {
1200 13 : log.printf(" Read in coefficents from file ");
1201 : }
1202 : else {
1203 1 : log.printf(" Read in coefficents from files:\n");
1204 : }
1205 46 : for(unsigned int i=0; i<ncoeffssets_; i++) {
1206 32 : IFile ifile;
1207 16 : ifile.link(*this);
1208 16 : if(use_mwalkers_mpi_ && mwalkers_mpi_single_files_) {
1209 8 : ifile.enforceSuffix("");
1210 : }
1211 32 : ifile.open(fnames[i]);
1212 48 : if(!ifile.FieldExist(coeffs_pntrs_[i]->getDataLabel())) {
1213 0 : std::string error_msg = "Problem with reading coefficents from file " + ifile.getPath() + ": no field with name " + coeffs_pntrs_[i]->getDataLabel() + "\n";
1214 0 : plumed_merror(error_msg);
1215 : }
1216 16 : size_t ncoeffs_read = coeffs_pntrs_[i]->readFromFile(ifile,false,false);
1217 16 : if(ncoeffssets_==1) {
1218 52 : log.printf("%s (read %zu of %zu values)\n", ifile.getPath().c_str(),ncoeffs_read,coeffs_pntrs_[i]->numberOfCoeffs());
1219 : }
1220 : else {
1221 12 : log.printf(" coefficent set %u: %s (read %zu of %zu values)\n",i,ifile.getPath().c_str(),ncoeffs_read,coeffs_pntrs_[i]->numberOfCoeffs());
1222 : }
1223 16 : ifile.close();
1224 16 : if(read_aux_coeffs) {
1225 15 : ifile.open(fnames[i]);
1226 45 : if(!ifile.FieldExist(aux_coeffs_pntrs_[i]->getDataLabel())) {
1227 0 : std::string error_msg = "Problem with reading coefficents from file " + ifile.getPath() + ": no field with name " + aux_coeffs_pntrs_[i]->getDataLabel() + "\n";
1228 0 : plumed_merror(error_msg);
1229 : }
1230 15 : aux_coeffs_pntrs_[i]->readFromFile(ifile,false,false);
1231 15 : ifile.close();
1232 : }
1233 : else {
1234 1 : AuxCoeffs(i).setValues( Coeffs(i) );
1235 : }
1236 : }
1237 14 : }
1238 :
1239 :
1240 77 : void Optimizer::addCoeffsSetIDsToFilenames(std::vector<std::string>& fnames, std::string& coeffssetid_prefix) {
1241 77 : if(ncoeffssets_==1) {return;}
1242 : //
1243 9 : if(fnames.size()==1) {
1244 0 : fnames.resize(ncoeffssets_,fnames[0]);
1245 : }
1246 9 : plumed_assert(fnames.size()==ncoeffssets_);
1247 : //
1248 63 : for(unsigned int i=0; i<ncoeffssets_; i++) {
1249 27 : std::string is=""; Tools::convert(i,is);
1250 162 : fnames[i] = FileBase::appendSuffix(fnames[i],"."+coeffssetid_prefix_+is);
1251 : }
1252 : }
1253 :
1254 :
1255 85 : void Optimizer::setAllCoeffsSetIterationCounters() {
1256 271 : for(unsigned int i=0; i<ncoeffssets_; i++) {
1257 186 : coeffs_pntrs_[i]->setIterationCounterAndTime(getIterationCounter(),getTime());
1258 93 : aux_coeffs_pntrs_[i]->setIterationCounterAndTime(getIterationCounter(),getTime());
1259 93 : gradient_pntrs_[i]->setIterationCounterAndTime(getIterationCounter(),getTime());
1260 93 : targetdist_averages_pntrs_[i]->setIterationCounterAndTime(getIterationCounter(),getTime());
1261 93 : if(use_hessian_) {
1262 0 : hessian_pntrs_[i]->setIterationCounterAndTime(getIterationCounter(),getTime());
1263 : }
1264 : }
1265 85 : }
1266 :
1267 :
1268 90 : std::string Optimizer::getIterationCounterStr(const int offset) const {
1269 : std::string s;
1270 90 : Tools::convert(getIterationCounter()+offset,s);
1271 90 : return s;
1272 : }
1273 :
1274 :
1275 71 : void Optimizer::writeBiasOutputFiles() const {
1276 225 : for(unsigned int i=0; i<nbiases_; i++) {
1277 154 : bias_pntrs_[i]->writeBiasToFile();
1278 : }
1279 71 : }
1280 :
1281 :
1282 71 : void Optimizer::writeFesOutputFiles() const {
1283 225 : for(unsigned int i=0; i<nbiases_; i++) {
1284 154 : bias_pntrs_[i]->writeFesToFile();
1285 : }
1286 71 : }
1287 :
1288 :
1289 17 : void Optimizer::writeFesProjOutputFiles() const {
1290 51 : for(unsigned int i=0; i<nbiases_; i++) {
1291 34 : bias_pntrs_[i]->writeFesProjToFile();
1292 : }
1293 17 : }
1294 :
1295 :
1296 33 : void Optimizer::writeTargetDistOutputFiles() const {
1297 99 : for(unsigned int i=0; i<nbiases_; i++) {
1298 66 : if(dynamic_targetdists_[i]) {
1299 33 : bias_pntrs_[i]->writeTargetDistToFile();
1300 : }
1301 : }
1302 33 : }
1303 :
1304 :
1305 4 : void Optimizer::writeTargetDistProjOutputFiles() const {
1306 12 : for(unsigned int i=0; i<nbiases_; i++) {
1307 8 : if(dynamic_targetdists_[i]) {
1308 4 : bias_pntrs_[i]->writeTargetDistProjToFile();
1309 : }
1310 : }
1311 4 : }
1312 :
1313 :
1314 :
1315 : }
1316 4839 : }
|
