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