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 80 : 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 160 : isFirstStep(true)
92 : {
93 80 : std::vector<std::string> bias_labels(0);
94 160 : parseVector("BIAS",bias_labels);
95 80 : plumed_massert(bias_labels.size()>0,"problem with BIAS keyword");
96 80 : nbiases_ = bias_labels.size();
97 : //
98 80 : std::string error_msg = "";
99 240 : bias_pntrs_ = VesTools::getPointersFromLabels<VesBias*>(bias_labels,plumed.getActionSet(),error_msg);
100 80 : if(error_msg.size()>0) {plumed_merror("Error in keyword BIAS of "+getName()+": "+error_msg);}
101 :
102 252 : for(unsigned int i=0; i<bias_pntrs_.size(); i++) {
103 86 : bias_pntrs_[i]->linkOptimizer(this);
104 : //
105 86 : std::vector<CoeffsVector*> pntrs_coeffs = bias_pntrs_[i]->getCoeffsPntrs();
106 86 : std::vector<CoeffsVector*> pntrs_gradient = bias_pntrs_[i]->getGradientPntrs();
107 86 : std::vector<CoeffsVector*> pntrs_targetdist_averages = bias_pntrs_[i]->getTargetDistAveragesPntrs();
108 86 : plumed_massert(pntrs_coeffs.size()==pntrs_gradient.size(),"something wrong in the coefficients and gradient passed from VES bias");
109 86 : plumed_massert(pntrs_coeffs.size()==pntrs_targetdist_averages.size(),"something wrong in the coefficients and target distribution averages passed from VES bias");
110 258 : for(unsigned int k=0; k<pntrs_coeffs.size(); k++) {
111 86 : plumed_massert(pntrs_coeffs[k] != NULL,"some coefficient is not linked correctly");
112 86 : plumed_massert(pntrs_gradient[k] != NULL,"some gradient is not linked correctly");
113 86 : plumed_massert(pntrs_targetdist_averages[k] != NULL,"some target distribution average is not linked correctly");
114 : pntrs_coeffs[k]->turnOnIterationCounter();
115 86 : coeffs_pntrs_.push_back(pntrs_coeffs[k]);
116 86 : pntrs_gradient[k]->turnOnIterationCounter();
117 86 : gradient_pntrs_.push_back(pntrs_gradient[k]);
118 86 : pntrs_targetdist_averages[k]->turnOnIterationCounter();
119 86 : targetdist_averages_pntrs_.push_back(pntrs_targetdist_averages[k]);
120 : //
121 86 : CoeffsVector* aux_coeffs_tmp = new CoeffsVector(*pntrs_coeffs[k]);
122 86 : std::string aux_label = pntrs_coeffs[k]->getLabel();
123 86 : if(aux_label.find("coeffs")!=std::string::npos) {
124 258 : aux_label.replace(aux_label.find("coeffs"), std::string("coeffs").length(), "aux_coeffs");
125 : }
126 : else {
127 : aux_label += "_aux";
128 : }
129 86 : aux_coeffs_tmp->setLabels(aux_label);
130 86 : aux_coeffs_pntrs_.push_back(aux_coeffs_tmp);
131 86 : AuxCoeffs(i).setValues( Coeffs(i) );
132 : }
133 : }
134 80 : ncoeffssets_ = coeffs_pntrs_.size();
135 80 : plumed_massert(aux_coeffs_pntrs_.size()==ncoeffssets_,"problems in linking aux coefficients");
136 80 : plumed_massert(gradient_pntrs_.size()==ncoeffssets_,"problems in linking gradients");
137 80 : setAllCoeffsSetIterationCounters();
138 :
139 :
140 : //
141 80 : identical_coeffs_shape_ = true;
142 86 : 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 160 : if(keywords.exists("STEPSIZE")) {
150 152 : plumed_assert(!keywords.exists("INITIAL_STEPSIZE"));
151 76 : fixed_stepsize_=true;
152 152 : parseMultipleValues("STEPSIZE",stepsizes_);
153 76 : setCurrentStepSizes(stepsizes_);
154 : }
155 160 : if(keywords.exists("INITIAL_STEPSIZE")) {
156 6 : plumed_assert(!keywords.exists("STEPSIZE"));
157 3 : fixed_stepsize_=false;
158 6 : parseMultipleValues("INITIAL_STEPSIZE",stepsizes_);
159 3 : setCurrentStepSizes(stepsizes_);
160 : }
161 : //
162 80 : if(ncoeffssets_==1) {
163 154 : log.printf(" optimizing VES bias %s with label %s: \n",bias_pntrs_[0]->getName().c_str(),bias_pntrs_[0]->getLabel().c_str());
164 154 : log.printf(" KbT: %f\n",bias_pntrs_[0]->getKbT());
165 154 : log.printf(" number of coefficients: %zu\n",coeffs_pntrs_[0]->numberOfCoeffs());
166 77 : if(stepsizes_.size()>0) {
167 76 : if(fixed_stepsize_) {log.printf(" using a constant step size of %f\n",stepsizes_[0]);}
168 3 : else {log.printf(" using an initial step size of %f\n",stepsizes_[0]);}
169 : }
170 : }
171 : else {
172 3 : log.printf(" optimizing %u coefficient sets from following %u VES biases:\n",ncoeffssets_,nbiases_);
173 9 : 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 9 : 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 18 : 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 160 : 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 160 : if(keywords.exists("MULTIPLE_WALKERS")) {
205 160 : parseFlag("MULTIPLE_WALKERS",use_mwalkers_mpi_);
206 80 : if(use_mwalkers_mpi_) {
207 : mw_single_files=true;
208 : }
209 : }
210 :
211 80 : int numwalkers=1;
212 80 : int walker_rank=0;
213 80 : if(comm.Get_rank()==0) {
214 71 : numwalkers = multi_sim_comm.Get_size();
215 71 : walker_rank = multi_sim_comm.Get_rank();
216 : }
217 80 : comm.Bcast(numwalkers,0);
218 80 : comm.Bcast(walker_rank,0);
219 80 : 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 80 : 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 80 : dynamic_targetdists_.resize(nbiases_,false);
232 160 : if(keywords.exists("TARGETDIST_STRIDE")) {
233 : bool need_stride = false;
234 83 : for(unsigned int i=0; i<nbiases_; i++) {
235 166 : dynamic_targetdists_[i] = bias_pntrs_[i]->dynamicTargetDistribution();
236 83 : if(dynamic_targetdists_[i]) {need_stride = true;}
237 : }
238 154 : parse("TARGETDIST_STRIDE",ustride_targetdist_);
239 77 : 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 77 : 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 77 : if(ustride_targetdist_>0) {
246 37 : if(nbiases_==1) {
247 37 : log.printf(" the target distribution will be updated very %u coefficient iterations\n",ustride_targetdist_);
248 : }
249 : else {
250 0 : log.printf(" the target distribution will be updated very %u coefficient 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 37 : log.printf(" See and cite ");
257 111 : log << plumed.cite("Valsson and Parrinello, J. Chem. Theory Comput. 11, 1996-2002 (2015)");
258 37 : log.printf("\n");
259 : }
260 : }
261 :
262 160 : 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 coefficient iterations\n",ustride_reweightfactor_);
276 : }
277 : }
278 :
279 160 : if(keywords.exists("MONITOR_INSTANTANEOUS_GRADIENT")) {
280 160 : parseFlag("MONITOR_INSTANTANEOUS_GRADIENT",monitor_instantaneous_gradient_);
281 : }
282 :
283 160 : if(keywords.exists("MONITOR_AVERAGE_GRADIENT")) {
284 75 : bool monitor_aver_gradient = false;
285 150 : parseFlag("MONITOR_AVERAGE_GRADIENT",monitor_aver_gradient);
286 75 : 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 4 : 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 80 : 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 154 : if(keywords.exists("COEFFS_SET_ID_PREFIX")) {
319 154 : parse("COEFFS_SET_ID_PREFIX",coeffssetid_prefix_);
320 : }
321 77 : if(coeffssetid_prefix_.size()>0) {
322 0 : plumed_merror("COEFFS_SET_ID_PREFIX should only be given if optimizing multiple coefficient sets");
323 : }
324 : }
325 :
326 160 : if(keywords.exists("INITIAL_COEFFS")) {
327 : std::vector<std::string> initial_coeffs_fnames;
328 160 : parseFilenames("INITIAL_COEFFS",initial_coeffs_fnames);
329 80 : 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 80 : }
337 : }
338 : //
339 :
340 80 : std::vector<std::string> coeffs_fnames;
341 160 : if(keywords.exists("COEFFS_FILE")) {
342 240 : parseFilenames("COEFFS_FILE",coeffs_fnames,"coeffs.data");
343 80 : bool start_opt_afresh=false;
344 160 : if(keywords.exists("START_OPTIMIZATION_AFRESH")) {
345 154 : parseFlag("START_OPTIMIZATION_AFRESH",start_opt_afresh);
346 78 : 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 160 : if(getRestart()) {
351 43 : for(unsigned int i=0; i<coeffs_fnames.size(); i++) {
352 15 : 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 15 : }
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 15 : 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 1 : 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 80 : std::string coeffs_wstride_tmpstr="";
385 160 : parse("COEFFS_OUTPUT",coeffs_wstride_tmpstr);
386 160 : if(coeffs_wstride_tmpstr!="OFF" && coeffs_wstride_tmpstr.size()>0) {
387 80 : Tools::convert(coeffs_wstride_tmpstr,coeffs_wstride_);
388 : }
389 80 : if(coeffs_wstride_tmpstr=="OFF") {
390 0 : coeffs_fnames.clear();
391 : }
392 80 : setupOFiles(coeffs_fnames,coeffsOFiles_,mw_single_files);
393 160 : parse("COEFFS_FMT",coeffs_output_fmt_);
394 80 : if(coeffs_output_fmt_.size()>0) {
395 85 : for(unsigned int i=0; i<ncoeffssets_; i++) {
396 170 : coeffs_pntrs_[i]->setOutputFmt(coeffs_output_fmt_);
397 : }
398 : }
399 160 : if(!getRestart()) {
400 209 : for(unsigned int i=0; i<coeffsOFiles_.size(); i++) {
401 213 : coeffs_pntrs_[i]->writeToFile(*coeffsOFiles_[i],aux_coeffs_pntrs_[i],false);
402 : }
403 : }
404 80 : if(coeffs_fnames.size()>0) {
405 80 : if(ncoeffssets_==1) {
406 308 : 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 21 : for(unsigned int i=0; i<coeffs_fnames.size(); i++) {
411 27 : 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 80 : std::vector<std::string> gradient_fnames;
421 160 : if(keywords.exists("GRADIENT_FILE")) {
422 160 : parseFilenames("GRADIENT_FILE",gradient_fnames);
423 160 : parse("GRADIENT_OUTPUT",gradient_wstride_);
424 :
425 80 : if(coeffs_fnames.size()>0) {
426 228 : for(unsigned int i=0; i<gradient_fnames.size(); i++) {
427 74 : plumed_massert(gradient_fnames[i]!=coeffs_fnames[i],"COEFFS_FILE and GRADIENT_FILE cannot be the same");
428 : }
429 : }
430 80 : setupOFiles(gradient_fnames,gradientOFiles_,mw_single_files);
431 160 : parse("GRADIENT_FMT",gradient_output_fmt_);
432 80 : if(gradient_output_fmt_.size()>0) {
433 74 : for(unsigned int i=0; i<ncoeffssets_; i++) {
434 148 : gradient_pntrs_[i]->setOutputFmt(gradient_output_fmt_);
435 : }
436 : }
437 :
438 80 : if(gradient_fnames.size()>0) {
439 68 : if(ncoeffssets_==1) {
440 260 : 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 21 : for(unsigned int i=0; i<gradient_fnames.size(); i++) {
445 27 : log.printf(" coefficient set %u: %s\n",i,gradientOFiles_[i]->getPath().c_str());
446 : }
447 : }
448 : }
449 : }
450 :
451 80 : std::vector<std::string> hessian_fnames;
452 160 : if(keywords.exists("HESSIAN_FILE")) {
453 150 : parseFilenames("HESSIAN_FILE",hessian_fnames);
454 150 : parse("HESSIAN_OUTPUT",hessian_wstride_);
455 :
456 75 : if(coeffs_fnames.size()>0) {
457 207 : for(unsigned int i=0; i<hessian_fnames.size(); i++) {
458 66 : plumed_massert(hessian_fnames[i]!=coeffs_fnames[i],"COEFFS_FILE and HESSIAN_FILE cannot be the same");
459 : }
460 : }
461 75 : if(gradient_fnames.size()>0) {
462 197 : for(unsigned int i=0; i<hessian_fnames.size(); i++) {
463 66 : plumed_massert(hessian_fnames[i]!=gradient_fnames[i],"GRADIENT_FILE and HESSIAN_FILE cannot be the same");
464 : }
465 : }
466 75 : setupOFiles(hessian_fnames,hessianOFiles_,mw_single_files);
467 150 : parse("HESSIAN_FMT",hessian_output_fmt_);
468 :
469 75 : if(hessian_fnames.size()>0) {
470 62 : if(ncoeffssets_==1) {
471 240 : 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 14 : for(unsigned int i=0; i<hessian_fnames.size(); i++) {
476 18 : log.printf(" coefficient set %u: %s\n",i,hessianOFiles_[i]->getPath().c_str());
477 : }
478 : }
479 : }
480 : }
481 :
482 :
483 : //
484 160 : if(keywords.exists("MASK_FILE")) {
485 : std::vector<std::string> mask_fnames_in;
486 158 : parseVector("MASK_FILE",mask_fnames_in);
487 79 : 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 82 : 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 79 : coeffs_mask_pntrs_.resize(ncoeffssets_);
496 85 : for(unsigned int i=0; i<ncoeffssets_; i++) {
497 255 : coeffs_mask_pntrs_[i] = new CoeffsVector(*coeffs_pntrs_[i]);
498 255 : coeffs_mask_pntrs_[i]->setLabels("mask");
499 85 : coeffs_mask_pntrs_[i]->setValues(1.0);
500 255 : coeffs_mask_pntrs_[i]->setOutputFmt("%f");
501 : }
502 :
503 79 : if(mask_fnames_in.size()>0) {
504 3 : if(ncoeffssets_==1) {
505 3 : size_t nread = coeffs_mask_pntrs_[0]->readFromFile(mask_fnames_in[0],true,true);
506 3 : log.printf(" read %zu values from mask file %s\n",nread,mask_fnames_in[0].c_str());
507 3 : size_t ndeactived = coeffs_mask_pntrs_[0]->countValues(0.0);
508 3 : 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 coefficient 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 79 : std::vector<std::string> mask_fnames_out;
522 158 : parseFilenames("OUTPUT_MASK_FILE",mask_fnames_out);
523 :
524 164 : for(unsigned int i=0; i<mask_fnames_out.size(); i++) {
525 3 : if(mask_fnames_in.size()>0) {
526 3 : plumed_massert(mask_fnames_out[i]!=mask_fnames_in[i],"MASK_FILE and OUTPUT_MASK_FILE cannot be the same");
527 : }
528 3 : OFile maskOFile;
529 3 : maskOFile.link(*this);
530 3 : maskOFile.enforceBackup();
531 3 : 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 3 : maskOFile.open(mask_fnames_out[i]);
539 3 : coeffs_mask_pntrs_[i]->writeToFile(maskOFile,true);
540 3 : maskOFile.close();
541 82 : }
542 : }
543 :
544 160 : if(getRestart() && ustride_targetdist_>0) {
545 8 : 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 80 : std::vector<std::string> targetdist_averages_fnames;
554 160 : if(keywords.exists("TARGETDIST_AVERAGES_FILE")) {
555 240 : parseFilenames("TARGETDIST_AVERAGES_FILE",targetdist_averages_fnames,"targetdist-averages.data");
556 160 : parse("TARGETDIST_AVERAGES_OUTPUT",targetdist_averages_wstride_);
557 :
558 80 : if(coeffs_fnames.size()>0) {
559 252 : for(unsigned int i=0; i<targetdist_averages_fnames.size(); i++) {
560 86 : plumed_massert(targetdist_averages_fnames[i]!=coeffs_fnames[i],"COEFFS_FILE and TARGETDIST_AVERAGES_FILE cannot be the same");
561 : }
562 : }
563 80 : if(gradient_fnames.size()>0) {
564 216 : for(unsigned int i=0; i<targetdist_averages_fnames.size(); i++) {
565 74 : plumed_massert(targetdist_averages_fnames[i]!=gradient_fnames[i],"GRADIENT_FILE and TARGETDIST_AVERAGES_FILE cannot be the same");
566 : }
567 : }
568 80 : if(hessian_fnames.size()>0) {
569 194 : for(unsigned int i=0; i<targetdist_averages_fnames.size(); i++) {
570 66 : plumed_massert(targetdist_averages_fnames[i]!=hessian_fnames[i],"HESSIAN_FILE and TARGETDIST_AVERAGES_FILE cannot be the same");
571 : }
572 : }
573 80 : setupOFiles(targetdist_averages_fnames,targetdist_averagesOFiles_,mw_single_files);
574 160 : parse("TARGETDIST_AVERAGES_FMT",targetdist_averages_output_fmt_);
575 80 : if(targetdist_averages_output_fmt_.size()>0) {
576 82 : for(unsigned int i=0; i<ncoeffssets_; i++) {
577 164 : targetdist_averages_pntrs_[i]->setOutputFmt(targetdist_averages_output_fmt_);
578 : }
579 : }
580 :
581 252 : for(unsigned int i=0; i<targetdist_averagesOFiles_.size(); i++) {
582 172 : targetdist_averages_pntrs_[i]->writeToFile(*targetdist_averagesOFiles_[i]);
583 : }
584 :
585 80 : if(targetdist_averages_wstride_==0) {
586 153 : for(unsigned int i=0; i<targetdist_averagesOFiles_.size(); i++) {
587 53 : targetdist_averagesOFiles_[i]->close();
588 53 : delete targetdist_averagesOFiles_[i];
589 : }
590 : targetdist_averagesOFiles_.clear();
591 : }
592 :
593 80 : 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 160 : if(keywords.exists("BIAS_OUTPUT")) {
608 160 : parse("BIAS_OUTPUT",bias_output_stride_);
609 80 : if(bias_output_stride_>0) {
610 79 : bias_output_active_=true;
611 164 : for(unsigned int i=0; i<nbiases_; i++) {
612 170 : bias_pntrs_[i]->enableBiasFileOutput();
613 85 : bias_pntrs_[i]->setupBiasFileOutput();
614 85 : bias_pntrs_[i]->writeBiasToFile();
615 : }
616 : }
617 : else {
618 1 : bias_output_active_=false;
619 1 : bias_output_stride_=1000;
620 : }
621 : }
622 :
623 160 : if(keywords.exists("FES_OUTPUT")) {
624 160 : parse("FES_OUTPUT",fes_output_stride_);
625 80 : if(fes_output_stride_>0) {
626 79 : fes_output_active_=true;
627 164 : for(unsigned int i=0; i<nbiases_; i++) {
628 170 : bias_pntrs_[i]->enableFesFileOutput();
629 85 : bias_pntrs_[i]->setupFesFileOutput();
630 85 : bias_pntrs_[i]->writeFesToFile();
631 : }
632 : }
633 : else {
634 1 : fes_output_active_=false;
635 1 : fes_output_stride_=1000;
636 : }
637 : }
638 :
639 160 : if(keywords.exists("FES_PROJ_OUTPUT")) {
640 160 : parse("FES_PROJ_OUTPUT",fesproj_output_stride_);
641 80 : if(fesproj_output_stride_>0) {
642 16 : fesproj_output_active_=true;
643 32 : 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 64 : fesproj_output_active_=false;
651 64 : fesproj_output_stride_=1000;
652 : }
653 : }
654 :
655 86 : for(unsigned int i=0; i<nbiases_; i++) {
656 270 : 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 160 : if(keywords.exists("TARGETDIST_OUTPUT")) {
665 154 : parse("TARGETDIST_OUTPUT",targetdist_output_stride_);
666 77 : if(targetdist_output_stride_>0) {
667 36 : 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 36 : 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 :
674 36 : targetdist_output_active_=true;
675 72 : for(unsigned int i=0; i<nbiases_; i++) {
676 72 : if(dynamic_targetdists_[i]) {
677 36 : bias_pntrs_[i]->enableDynamicTargetDistFileOutput();
678 36 : bias_pntrs_[i]->setupTargetDistFileOutput();
679 36 : bias_pntrs_[i]->writeTargetDistToFile();
680 : }
681 : }
682 : }
683 : else {
684 41 : targetdist_output_active_=false;
685 41 : targetdist_output_stride_=1000;
686 : }
687 : }
688 :
689 160 : if(keywords.exists("TARGETDIST_PROJ_OUTPUT")) {
690 154 : parse("TARGETDIST_PROJ_OUTPUT",targetdist_proj_output_stride_);
691 77 : if(targetdist_proj_output_stride_>0) {
692 3 : if(ustride_targetdist_==0) {
693 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");
694 : }
695 3 : if(targetdist_proj_output_stride_%ustride_targetdist_!=0) {
696 0 : plumed_merror("the value given in TARGETDIST_PROJ_OUTPUT doesn't make sense, it should be multiple of TARGETDIST_STRIDE");
697 : }
698 :
699 3 : targetdist_proj_output_active_=true;
700 6 : for(unsigned int i=0; i<nbiases_; i++) {
701 6 : if(dynamic_targetdists_[i]) {
702 3 : bias_pntrs_[i]->enableDynamicTargetDistFileOutput();
703 3 : bias_pntrs_[i]->setupTargetDistProjFileOutput();
704 3 : bias_pntrs_[i]->writeTargetDistProjToFile();
705 : }
706 : }
707 : }
708 : else {
709 74 : targetdist_proj_output_active_=false;
710 74 : targetdist_proj_output_stride_=1000;
711 : }
712 : }
713 :
714 80 : if(ncoeffssets_==1) {
715 77 : log.printf(" Output Components:\n");
716 77 : log.printf(" ");
717 77 : if(monitor_instantaneous_gradient_) {
718 6 : addComponent("gradrms"); componentIsNotPeriodic("gradrms");
719 2 : log.printf(" ");
720 6 : addComponent("gradmax"); componentIsNotPeriodic("gradmax");
721 : }
722 77 : if(aver_gradient_pntrs_.size()>0) {
723 2 : log.printf(" ");
724 6 : addComponent("avergradrms"); componentIsNotPeriodic("avergradrms");
725 2 : log.printf(" ");
726 6 : addComponent("avergradmax"); componentIsNotPeriodic("avergradmax");
727 : }
728 77 : if(!fixed_stepsize_) {
729 3 : log.printf(" ");
730 9 : addComponent("stepsize"); componentIsNotPeriodic("stepsize");
731 6 : getPntrToComponent("stepsize")->set( getCurrentStepSize(0) );
732 : }
733 : }
734 : else {
735 9 : for(unsigned int i=0; i<ncoeffssets_; i++) {
736 9 : log.printf(" Output Components for coefficient set %u:\n",i);
737 27 : std::string is=""; Tools::convert(i,is); is = "-" + coeffssetid_prefix_ + is;
738 9 : log.printf(" ");
739 9 : if(monitor_instantaneous_gradient_) {
740 9 : addComponent("gradrms"+is); componentIsNotPeriodic("gradrms"+is);
741 3 : log.printf(" ");
742 9 : addComponent("gradmax"+is); componentIsNotPeriodic("gradmax"+is);
743 : }
744 9 : if(aver_gradient_pntrs_.size()>0) {
745 0 : log.printf(" ");
746 0 : addComponent("avergradrms"+is); componentIsNotPeriodic("avergradrms"+is);
747 0 : log.printf(" ");
748 0 : addComponent("avergradmax"+is); componentIsNotPeriodic("avergradmax"+is);
749 : }
750 9 : if(!fixed_stepsize_) {
751 0 : log.printf(" ");
752 0 : addComponent("stepsize"+is); componentIsNotPeriodic("stepsize"+is);
753 0 : getPntrToComponent("stepsize"+is)->set( getCurrentStepSize(i) );
754 : }
755 : }
756 80 : }
757 :
758 80 : }
759 :
760 :
761 160 : Optimizer::~Optimizer() {
762 : //
763 80 : if(!isTargetDistOutputActive()) {
764 50 : for(unsigned int i=0; i<nbiases_; i++) {
765 100 : if(dynamic_targetdists_[i]) {
766 1 : bias_pntrs_[i]->enableDynamicTargetDistFileOutput();
767 1 : bias_pntrs_[i]->setupTargetDistFileOutput();
768 1 : bias_pntrs_[i]->writeTargetDistToFile();
769 : }
770 : }
771 : }
772 72 : else if(isTargetDistOutputActive() && getIterationCounter()%getTargetDistOutputStride()!=0) {
773 1 : writeTargetDistOutputFiles();
774 : }
775 : //
776 252 : for(unsigned int i=0; i<aux_coeffs_pntrs_.size(); i++) {
777 86 : delete aux_coeffs_pntrs_[i];
778 : }
779 : aux_coeffs_pntrs_.clear();
780 : //
781 164 : for(unsigned int i=0; i<aver_gradient_pntrs_.size(); i++) {
782 2 : delete aver_gradient_pntrs_[i];
783 : }
784 : aver_gradient_pntrs_.clear();
785 : //
786 330 : for(unsigned int i=0; i<coeffs_mask_pntrs_.size(); i++) {
787 85 : delete coeffs_mask_pntrs_[i];
788 : }
789 : coeffs_mask_pntrs_.clear();
790 : //
791 330 : for(unsigned int i=0; i<coeffsOFiles_.size(); i++) {
792 85 : coeffsOFiles_[i]->close();
793 85 : delete coeffsOFiles_[i];
794 : }
795 : coeffsOFiles_.clear();
796 308 : for(unsigned int i=0; i<gradientOFiles_.size(); i++) {
797 74 : gradientOFiles_[i]->close();
798 74 : delete gradientOFiles_[i];
799 : }
800 : gradientOFiles_.clear();
801 292 : for(unsigned int i=0; i<hessianOFiles_.size(); i++) {
802 66 : hessianOFiles_[i]->close();
803 66 : delete hessianOFiles_[i];
804 : }
805 : hessianOFiles_.clear();
806 226 : for(unsigned int i=0; i<targetdist_averagesOFiles_.size(); i++) {
807 33 : targetdist_averagesOFiles_[i]->close();
808 33 : delete targetdist_averagesOFiles_[i];
809 : }
810 : targetdist_averagesOFiles_.clear();
811 80 : }
812 :
813 :
814 84 : void Optimizer::registerKeywords( Keywords& keys ) {
815 84 : Action::registerKeywords(keys);
816 84 : ActionPilot::registerKeywords(keys);
817 84 : ActionWithValue::registerKeywords(keys);
818 : //
819 168 : keys.remove("NUMERICAL_DERIVATIVES");
820 : // Default always active keywords
821 336 : keys.add("compulsory","BIAS","the label of the VES bias to be optimized");
822 336 : keys.add("compulsory","STRIDE","the frequency of updating the coefficients given in the number of MD steps.");
823 420 : keys.add("compulsory","COEFFS_FILE","coeffs.data","the name of output file for the coefficients");
824 420 : keys.add("compulsory","COEFFS_OUTPUT","100","how often the coefficients should be written to file. This parameter is given as the number of iterations.");
825 336 : keys.add("optional","COEFFS_FMT","specify format for coefficient file(s) (useful for decrease the number of digits in regtests)");
826 336 : 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.");
827 : //
828 336 : keys.add("optional","INITIAL_COEFFS","the name(s) of file(s) with the initial coefficients");
829 : // Hidden keywords to output the gradient to a file.
830 336 : keys.add("hidden","GRADIENT_FILE","the name of output file for the gradient");
831 336 : 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");
832 336 : keys.add("hidden","GRADIENT_FMT","specify format for gradient file(s) (useful for decrease the number of digits in regtests)");
833 : // Either use a fixed stepsize (useFixedStepSizeKeywords) or changing stepsize (useDynamicsStepSizeKeywords)
834 336 : keys.reserve("compulsory","STEPSIZE","the step size used for the optimization");
835 336 : keys.reserve("compulsory","INITIAL_STEPSIZE","the initial step size used for the optimization");
836 : // Keywords related to the Hessian, actived with the useHessianKeywords function
837 252 : 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");
838 336 : keys.reserve("hidden","HESSIAN_FILE","the name of output file for the Hessian");
839 336 : 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");
840 336 : keys.reserve("hidden","HESSIAN_FMT","specify format for hessian file(s) (useful for decrease the number of digits in regtests)");
841 : // Keywords related to the multiple walkers, actived with the useMultipleWalkersKeywords function
842 252 : keys.reserveFlag("MULTIPLE_WALKERS",false,"if optimization is to be performed using multiple walkers connected via MPI");
843 : // Keywords related to the mask file, actived with the useMaskKeywords function
844 336 : 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.");
845 336 : 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.");
846 : //
847 252 : 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.");
848 : //
849 252 : keys.addFlag("MONITOR_INSTANTANEOUS_GRADIENT",false,"if quantities related to the instantaneous gradient should be outputted.");
850 : //
851 252 : keys.reserveFlag("MONITOR_AVERAGE_GRADIENT",false,"if the averaged gradient should be monitored and quantities related to it should be outputted.");
852 336 : keys.reserve("optional","MONITOR_AVERAGES_GRADIENT_EXP_DECAY","use an exponentially decaying averaging with a given time constant when monitoring the averaged gradient");
853 : //
854 336 : 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.");
855 336 : 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.");
856 336 : 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.");
857 : //
858 336 : keys.add("optional","TARGETDIST_AVERAGES_FILE","the name of output file for the target distribution averages. By default it is targetdist-averages.data.");
859 336 : 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");
860 336 : keys.add("hidden","TARGETDIST_AVERAGES_FMT","specify format for target distribution averages file(s) (useful for decrease the number of digits in regtests)");
861 : //
862 336 : 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.");
863 336 : 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.");
864 336 : 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.");
865 : //
866 336 : 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.");
867 : //
868 168 : keys.use("RESTART");
869 : //
870 168 : keys.use("UPDATE_FROM");
871 168 : keys.use("UPDATE_UNTIL");
872 : // Components that are always active
873 336 : 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-#.");
874 336 : 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-#.");
875 : // keys.addOutputComponent("gradmaxidx","default","the index of the maximum absolute value of the gradient");
876 :
877 84 : }
878 :
879 :
880 77 : void Optimizer::useHessianKeywords(Keywords& keys) {
881 : // keys.use("FULL_HESSIAN");
882 154 : keys.use("HESSIAN_FILE");
883 154 : keys.use("HESSIAN_OUTPUT");
884 154 : keys.use("HESSIAN_FMT");
885 77 : }
886 :
887 :
888 84 : void Optimizer::useMultipleWalkersKeywords(Keywords& keys) {
889 168 : keys.use("MULTIPLE_WALKERS");
890 84 : }
891 :
892 :
893 78 : void Optimizer::useFixedStepSizeKeywords(Keywords& keys) {
894 156 : keys.use("STEPSIZE");
895 78 : }
896 :
897 :
898 4 : void Optimizer::useDynamicStepSizeKeywords(Keywords& keys) {
899 8 : keys.use("INITIAL_STEPSIZE");
900 16 : 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-#.");
901 4 : }
902 :
903 :
904 82 : void Optimizer::useMaskKeywords(Keywords& keys) {
905 164 : keys.use("MASK_FILE");
906 164 : keys.use("OUTPUT_MASK_FILE");
907 82 : }
908 :
909 :
910 79 : void Optimizer::useRestartKeywords(Keywords& keys) {
911 158 : keys.use("START_OPTIMIZATION_AFRESH");
912 79 : }
913 :
914 :
915 77 : void Optimizer::useMonitorAverageGradientKeywords(Keywords& keys) {
916 154 : keys.use("MONITOR_AVERAGE_GRADIENT");
917 154 : keys.use("MONITOR_AVERAGES_GRADIENT_EXP_DECAY");
918 308 : 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-#.");
919 308 : 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-#.");
920 77 : }
921 :
922 :
923 79 : void Optimizer::useDynamicTargetDistributionKeywords(Keywords& keys) {
924 158 : keys.use("TARGETDIST_STRIDE");
925 158 : keys.use("TARGETDIST_OUTPUT");
926 158 : keys.use("TARGETDIST_PROJ_OUTPUT");
927 79 : }
928 :
929 :
930 0 : void Optimizer::useReweightFactorKeywords(Keywords& keys) {
931 0 : keys.use("REWEIGHT_FACTOR_STRIDE");
932 0 : }
933 :
934 :
935 75 : void Optimizer::turnOnHessian() {
936 75 : plumed_massert(hessian_pntrs_.size()==0,"turnOnHessian() should only be run during initialization");
937 75 : use_hessian_=true;
938 : hessian_pntrs_.clear();
939 156 : for(unsigned int i=0; i<nbiases_; i++) {
940 162 : std::vector<CoeffsMatrix*> pntrs_hessian = enableHessian(bias_pntrs_[i],diagonal_hessian_);
941 243 : for(unsigned int k=0; k<pntrs_hessian.size(); k++) {
942 81 : pntrs_hessian[k]->turnOnIterationCounter();
943 162 : pntrs_hessian[k]->setIterationCounterAndTime(getIterationCounter(),getTime());
944 81 : hessian_pntrs_.push_back(pntrs_hessian[k]);
945 : }
946 : }
947 75 : plumed_massert(hessian_pntrs_.size()==ncoeffssets_,"problems in linking Hessians");
948 75 : if(diagonal_hessian_) {
949 75 : log.printf(" Optimization performed using diagonal Hessian matrix\n");
950 : }
951 : else {
952 0 : log.printf(" Optimization performed using full Hessian matrix\n");
953 : }
954 : //
955 75 : if(hessian_output_fmt_.size()>0) {
956 66 : for(unsigned int i=0; i<ncoeffssets_; i++) {
957 132 : hessian_pntrs_[i]->setOutputFmt(hessian_output_fmt_);
958 : }
959 : }
960 :
961 75 : }
962 :
963 :
964 0 : void Optimizer::turnOffHessian() {
965 0 : use_hessian_=false;
966 0 : for(unsigned int i=0; i<nbiases_; i++) {
967 0 : bias_pntrs_[i]->disableHessian();
968 : }
969 : hessian_pntrs_.clear();
970 0 : for(unsigned int i=0; i<hessianOFiles_.size(); i++) {
971 0 : hessianOFiles_[i]->close();
972 0 : delete hessianOFiles_[i];
973 : }
974 : hessianOFiles_.clear();
975 0 : }
976 :
977 :
978 81 : std::vector<CoeffsMatrix*> Optimizer::enableHessian(VesBias* bias_pntr_in, const bool diagonal_hessian) {
979 81 : plumed_massert(use_hessian_,"the Hessian should not be used");
980 81 : bias_pntr_in->enableHessian(diagonal_hessian);
981 : std::vector<CoeffsMatrix*> hessian_pntrs_out = bias_pntr_in->getHessianPntrs();
982 324 : for(unsigned int k=0; k<hessian_pntrs_out.size(); k++) {
983 81 : plumed_massert(hessian_pntrs_out[k] != NULL,"Hessian is needed but not linked correctly");
984 : }
985 81 : return hessian_pntrs_out;
986 : }
987 :
988 :
989 : // CoeffsMatrix* Optimizer::switchToDiagonalHessian(VesBias* bias_pntr_in) {
990 : // plumed_massert(use_hessian_,"it does not make sense to switch to diagonal Hessian if it Hessian is not used");
991 : // diagonal_hessian_=true;
992 : // bias_pntr_in->enableHessian(diagonal_hessian_);
993 : // CoeffsMatrix* hessian_pntr_out = bias_pntr_in->getHessianPntr();
994 : // plumed_massert(hessian_pntr_out != NULL,"Hessian is needed but not linked correctly");
995 : // //
996 : // 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());
997 : // return hessian_pntr_out;
998 : // }
999 :
1000 :
1001 : // CoeffsMatrix* Optimizer::switchToFullHessian(VesBias* bias_pntr_in) {
1002 : // plumed_massert(use_hessian_,"it does not make sense to switch to diagonal Hessian if it Hessian is not used");
1003 : // diagonal_hessian_=false;
1004 : // bias_pntr_in->enableHessian(diagonal_hessian_);
1005 : // CoeffsMatrix* hessian_pntr_out = bias_pntr_in->getHessianPntr();
1006 : // plumed_massert(hessian_pntr_out != NULL,"Hessian is needed but not linked correctly");
1007 : // //
1008 : // 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());
1009 : // return hessian_pntr_out;
1010 : // }
1011 :
1012 :
1013 40910 : void Optimizer::update() {
1014 40910 : if(onStep() && !isFirstStep) {
1015 40820 : for(unsigned int i=0; i<nbiases_; i++) {
1016 81640 : bias_pntrs_[i]->updateGradientAndHessian(use_mwalkers_mpi_);
1017 : }
1018 40820 : for(unsigned int i=0; i<ncoeffssets_; i++) {
1019 122460 : if(gradient_pntrs_[i]->isActive()) {coeffsUpdate(i);}
1020 : else {
1021 500 : std::string msg = "iteration " + getIterationCounterStr(+1) +
1022 250 : " for " + bias_pntrs_[i]->getLabel() +
1023 : " - the coefficients are not updated as CV values are outside the bias intervals";
1024 125 : warning(msg);
1025 : }
1026 :
1027 : // +1 as this is done before increaseIterationCounter() is used
1028 40820 : unsigned int curr_iter = getIterationCounter()+1;
1029 40820 : double curr_time = getTime();
1030 40820 : coeffs_pntrs_[i]->setIterationCounterAndTime(curr_iter,curr_time);
1031 40820 : aux_coeffs_pntrs_[i]->setIterationCounterAndTime(curr_iter,curr_time);
1032 40820 : gradient_pntrs_[i]->setIterationCounterAndTime(curr_iter,curr_time);
1033 40820 : targetdist_averages_pntrs_[i]->setIterationCounterAndTime(curr_iter,curr_time);
1034 40820 : if(use_hessian_) {
1035 40770 : hessian_pntrs_[i]->setIterationCounterAndTime(curr_iter,curr_time);
1036 : }
1037 40820 : if(aver_gradient_pntrs_.size()>0) {
1038 20 : aver_gradient_pntrs_[i]->setIterationCounterAndTime(curr_iter,curr_time);
1039 40 : aver_gradient_pntrs_[i]->addToAverage(*gradient_pntrs_[i]);
1040 : }
1041 : }
1042 : increaseIterationCounter();
1043 40760 : updateOutputComponents();
1044 81580 : for(unsigned int i=0; i<ncoeffssets_; i++) {
1045 40820 : writeOutputFiles(i);
1046 : }
1047 81090 : if(ustride_targetdist_>0 && getIterationCounter()%ustride_targetdist_==0) {
1048 346 : for(unsigned int i=0; i<nbiases_; i++) {
1049 692 : if(dynamic_targetdists_[i]) {
1050 346 : bias_pntrs_[i]->updateTargetDistributions();
1051 : }
1052 : }
1053 : }
1054 40760 : if(ustride_reweightfactor_>0 && getIterationCounter()%ustride_reweightfactor_==0) {
1055 0 : for(unsigned int i=0; i<nbiases_; i++) {
1056 0 : bias_pntrs_[i]->updateReweightFactor();
1057 : }
1058 : }
1059 :
1060 :
1061 : //
1062 81510 : if(isBiasOutputActive() && getIterationCounter()%getBiasOutputStride()==0) {
1063 79 : writeBiasOutputFiles();
1064 : }
1065 81510 : if(isFesOutputActive() && getIterationCounter()%getFesOutputStride()==0) {
1066 79 : writeFesOutputFiles();
1067 : }
1068 40920 : if(isFesProjOutputActive() && getIterationCounter()%getFesProjOutputStride()==0) {
1069 16 : writeFesProjOutputFiles();
1070 : }
1071 81080 : if(isTargetDistOutputActive() && getIterationCounter()%getTargetDistOutputStride()==0) {
1072 35 : writeTargetDistOutputFiles();
1073 : }
1074 40790 : if(isTargetDistProjOutputActive() && getIterationCounter()%getTargetDistProjOutputStride()==0) {
1075 3 : writeTargetDistProjOutputFiles();
1076 : }
1077 : }
1078 : else {
1079 150 : isFirstStep=false;
1080 : }
1081 40910 : }
1082 :
1083 :
1084 40760 : void Optimizer::updateOutputComponents() {
1085 40760 : if(ncoeffssets_==1) {
1086 40730 : if(!fixed_stepsize_) {
1087 60 : getPntrToComponent("stepsize")->set( getCurrentStepSize(0) );
1088 : }
1089 40730 : if(monitor_instantaneous_gradient_) {
1090 40 : getPntrToComponent("gradrms")->set( gradient_pntrs_[0]->getRMS() );
1091 20 : size_t gradient_maxabs_idx=0;
1092 40 : getPntrToComponent("gradmax")->set( gradient_pntrs_[0]->getMaxAbsValue(gradient_maxabs_idx) );
1093 : }
1094 40730 : if(aver_gradient_pntrs_.size()>0) {
1095 40 : getPntrToComponent("avergradrms")->set( aver_gradient_pntrs_[0]->getRMS() );
1096 20 : size_t avergradient_maxabs_idx=0;
1097 40 : getPntrToComponent("avergradmax")->set( aver_gradient_pntrs_[0]->getMaxAbsValue(avergradient_maxabs_idx) );
1098 : }
1099 : }
1100 : else {
1101 90 : for(unsigned int i=0; i<ncoeffssets_; i++) {
1102 270 : std::string is=""; Tools::convert(i,is); is = "-" + coeffssetid_prefix_ + is;
1103 90 : if(!fixed_stepsize_) {
1104 0 : getPntrToComponent("stepsize"+is)->set( getCurrentStepSize(i) );
1105 : }
1106 90 : if(monitor_instantaneous_gradient_) {
1107 90 : getPntrToComponent("gradrms"+is)->set( gradient_pntrs_[i]->getRMS() );
1108 30 : size_t gradient_maxabs_idx=0;
1109 60 : getPntrToComponent("gradmax"+is)->set( gradient_pntrs_[i]->getMaxAbsValue(gradient_maxabs_idx) );
1110 : }
1111 90 : if(aver_gradient_pntrs_.size()>0) {
1112 0 : getPntrToComponent("avergradrms"+is)->set( aver_gradient_pntrs_[0]->getRMS() );
1113 0 : size_t avergradient_maxabs_idx=0;
1114 0 : getPntrToComponent("avergradmax"+is)->set( aver_gradient_pntrs_[0]->getMaxAbsValue(avergradient_maxabs_idx) );
1115 : }
1116 : }
1117 : }
1118 40760 : }
1119 :
1120 :
1121 1 : void Optimizer::turnOffCoeffsOutputFiles() {
1122 4 : for(unsigned int i=0; i<coeffsOFiles_.size(); i++) {
1123 1 : coeffsOFiles_[i]->close();
1124 1 : delete coeffsOFiles_[i];
1125 : }
1126 : coeffsOFiles_.clear();
1127 1 : }
1128 :
1129 :
1130 40820 : void Optimizer::writeOutputFiles(const unsigned int coeffs_id) {
1131 40820 : if(coeffsOFiles_.size()>0 && iter_counter%coeffs_wstride_==0) {
1132 3232 : coeffs_pntrs_[coeffs_id]->writeToFile(*coeffsOFiles_[coeffs_id],aux_coeffs_pntrs_[coeffs_id],false);
1133 : }
1134 40820 : if(gradientOFiles_.size()>0 && iter_counter%gradient_wstride_==0) {
1135 740 : if(aver_gradient_pntrs_.size()==0) {
1136 2160 : gradient_pntrs_[coeffs_id]->writeToFile(*gradientOFiles_[coeffs_id],false);
1137 : }
1138 : else {
1139 80 : gradient_pntrs_[coeffs_id]->writeToFile(*gradientOFiles_[coeffs_id],aver_gradient_pntrs_[coeffs_id],false);
1140 : }
1141 : }
1142 40820 : if(hessianOFiles_.size()>0 && iter_counter%hessian_wstride_==0) {
1143 1980 : hessian_pntrs_[coeffs_id]->writeToFile(*hessianOFiles_[coeffs_id]);
1144 : }
1145 40820 : if(targetdist_averagesOFiles_.size()>0 && iter_counter%targetdist_averages_wstride_==0) {
1146 990 : targetdist_averages_pntrs_[coeffs_id]->writeToFile(*targetdist_averagesOFiles_[coeffs_id]);
1147 : }
1148 40820 : }
1149 :
1150 :
1151 389 : void Optimizer::setupOFiles(std::vector<std::string>& fnames, std::vector<OFile*>& OFiles, const bool multi_sim_single_files) {
1152 389 : plumed_assert(ncoeffssets_>0);
1153 778 : OFiles.resize(fnames.size(),NULL);
1154 1414 : for(unsigned int i=0; i<fnames.size(); i++) {
1155 636 : OFiles[i] = new OFile();
1156 636 : OFiles[i]->link(*this);
1157 318 : if(multi_sim_single_files) {
1158 48 : unsigned int r=0;
1159 48 : if(comm.Get_rank()==0) {r=multi_sim_comm.Get_rank();}
1160 48 : comm.Bcast(r,0);
1161 48 : if(r>0) {fnames[i]="/dev/null";}
1162 144 : OFiles[i]->enforceSuffix("");
1163 : }
1164 636 : OFiles[i]->open(fnames[i]);
1165 318 : OFiles[i]->setHeavyFlush();
1166 : }
1167 389 : }
1168 :
1169 :
1170 14 : void Optimizer::readCoeffsFromFiles(const std::vector<std::string>& fnames, const bool read_aux_coeffs) {
1171 14 : plumed_assert(ncoeffssets_>0);
1172 14 : plumed_assert(fnames.size()==ncoeffssets_);
1173 14 : if(ncoeffssets_==1) {
1174 13 : log.printf(" Read in coefficients from file ");
1175 : }
1176 : else {
1177 1 : log.printf(" Read in coefficients from files:\n");
1178 : }
1179 16 : for(unsigned int i=0; i<ncoeffssets_; i++) {
1180 16 : IFile ifile;
1181 16 : ifile.link(*this);
1182 16 : if(use_mwalkers_mpi_ && mwalkers_mpi_single_files_) {
1183 8 : ifile.enforceSuffix("");
1184 : }
1185 32 : ifile.open(fnames[i]);
1186 48 : if(!ifile.FieldExist(coeffs_pntrs_[i]->getDataLabel())) {
1187 0 : std::string error_msg = "Problem with reading coefficients from file " + ifile.getPath() + ": no field with name " + coeffs_pntrs_[i]->getDataLabel() + "\n";
1188 0 : plumed_merror(error_msg);
1189 : }
1190 16 : size_t ncoeffs_read = coeffs_pntrs_[i]->readFromFile(ifile,false,false);
1191 16 : if(ncoeffssets_==1) {
1192 39 : log.printf("%s (read %zu of %zu values)\n", ifile.getPath().c_str(),ncoeffs_read,coeffs_pntrs_[i]->numberOfCoeffs());
1193 : }
1194 : else {
1195 9 : log.printf(" coefficient set %u: %s (read %zu of %zu values)\n",i,ifile.getPath().c_str(),ncoeffs_read,coeffs_pntrs_[i]->numberOfCoeffs());
1196 : }
1197 16 : ifile.close();
1198 16 : if(read_aux_coeffs) {
1199 15 : ifile.open(fnames[i]);
1200 45 : if(!ifile.FieldExist(aux_coeffs_pntrs_[i]->getDataLabel())) {
1201 0 : std::string error_msg = "Problem with reading coefficients from file " + ifile.getPath() + ": no field with name " + aux_coeffs_pntrs_[i]->getDataLabel() + "\n";
1202 0 : plumed_merror(error_msg);
1203 : }
1204 15 : aux_coeffs_pntrs_[i]->readFromFile(ifile,false,false);
1205 15 : ifile.close();
1206 : }
1207 : else {
1208 1 : AuxCoeffs(i).setValues( Coeffs(i) );
1209 : }
1210 16 : }
1211 14 : }
1212 :
1213 :
1214 86 : void Optimizer::addCoeffsSetIDsToFilenames(std::vector<std::string>& fnames, std::string& coeffssetid_prefix) {
1215 172 : if(ncoeffssets_==1) {return;}
1216 : //
1217 9 : if(fnames.size()==1) {
1218 0 : fnames.resize(ncoeffssets_,fnames[0]);
1219 : }
1220 9 : plumed_assert(fnames.size()==ncoeffssets_);
1221 : //
1222 27 : for(unsigned int i=0; i<ncoeffssets_; i++) {
1223 27 : std::string is=""; Tools::convert(i,is);
1224 162 : fnames[i] = FileBase::appendSuffix(fnames[i],"."+coeffssetid_prefix_+is);
1225 : }
1226 : }
1227 :
1228 :
1229 94 : void Optimizer::setAllCoeffsSetIterationCounters() {
1230 196 : for(unsigned int i=0; i<ncoeffssets_; i++) {
1231 306 : coeffs_pntrs_[i]->setIterationCounterAndTime(getIterationCounter(),getTime());
1232 204 : aux_coeffs_pntrs_[i]->setIterationCounterAndTime(getIterationCounter(),getTime());
1233 204 : gradient_pntrs_[i]->setIterationCounterAndTime(getIterationCounter(),getTime());
1234 204 : targetdist_averages_pntrs_[i]->setIterationCounterAndTime(getIterationCounter(),getTime());
1235 102 : if(use_hessian_) {
1236 0 : hessian_pntrs_[i]->setIterationCounterAndTime(getIterationCounter(),getTime());
1237 : }
1238 : }
1239 94 : }
1240 :
1241 :
1242 125 : std::string Optimizer::getIterationCounterStr(const int offset) const {
1243 : std::string s;
1244 125 : Tools::convert(getIterationCounter()+offset,s);
1245 125 : return s;
1246 : }
1247 :
1248 :
1249 79 : void Optimizer::writeBiasOutputFiles() const {
1250 164 : for(unsigned int i=0; i<nbiases_; i++) {
1251 170 : bias_pntrs_[i]->writeBiasToFile();
1252 : }
1253 79 : }
1254 :
1255 :
1256 79 : void Optimizer::writeFesOutputFiles() const {
1257 164 : for(unsigned int i=0; i<nbiases_; i++) {
1258 170 : bias_pntrs_[i]->writeFesToFile();
1259 : }
1260 79 : }
1261 :
1262 :
1263 16 : void Optimizer::writeFesProjOutputFiles() const {
1264 32 : for(unsigned int i=0; i<nbiases_; i++) {
1265 32 : bias_pntrs_[i]->writeFesProjToFile();
1266 : }
1267 16 : }
1268 :
1269 :
1270 36 : void Optimizer::writeTargetDistOutputFiles() const {
1271 72 : for(unsigned int i=0; i<nbiases_; i++) {
1272 72 : if(dynamic_targetdists_[i]) {
1273 36 : bias_pntrs_[i]->writeTargetDistToFile();
1274 : }
1275 : }
1276 36 : }
1277 :
1278 :
1279 3 : void Optimizer::writeTargetDistProjOutputFiles() const {
1280 6 : for(unsigned int i=0; i<nbiases_; i++) {
1281 6 : if(dynamic_targetdists_[i]) {
1282 3 : bias_pntrs_[i]->writeTargetDistProjToFile();
1283 : }
1284 : }
1285 3 : }
1286 :
1287 :
1288 :
1289 : }
1290 5874 : }
|
