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