Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2017-2026 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed 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 : plumed 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 plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #ifndef __PLUMED_isdb_MetainferenceCore_h
23 : #define __PLUMED_isdb_MetainferenceCore_h
24 :
25 : #include "core/Action.h"
26 : #include "core/ActionWithValue.h"
27 : #include "tools/File.h"
28 : #include "tools/Random.h"
29 : #include "tools/OpenMP.h"
30 : #include <functional>
31 : #include <string>
32 : #include <vector>
33 :
34 : namespace PLMD {
35 :
36 : class Communicator;
37 : class Keywords;
38 : class Log;
39 : class Value;
40 :
41 : namespace isdb {
42 :
43 : /// History of the estimates of the squared standard error of the mean of one datum:
44 : /// a circular buffer of fixed capacity that keeps track of its maximum
45 12070 : class SigmaMeanHistory {
46 : std::vector<double> buf_;
47 : unsigned first_=0;
48 : unsigned n_=0;
49 : double max_=0.;
50 1162 : void recomputeMax() {
51 1162 : max_=get(0);
52 1410 : for(unsigned k=1; k<n_; ++k) {
53 248 : if(max_<get(k)) {
54 100 : max_=get(k);
55 : }
56 : }
57 1162 : }
58 : public:
59 192 : explicit SigmaMeanHistory(const unsigned capacity=1): buf_(capacity>0?capacity:1) {}
60 : unsigned size() const {
61 7464 : return n_;
62 : }
63 : /// the maximum of the stored values (the history must not be empty)
64 : double max() const {
65 9751 : return max_;
66 : }
67 : /// the k-th value, from the oldest
68 : double get(const unsigned k) const {
69 21719 : return buf_[(first_+k)%buf_.size()];
70 : }
71 : void clear() {
72 448 : first_=0;
73 448 : n_=0;
74 : }
75 : /// add a value, removing the oldest one if the history is full
76 8787 : void push(const double v) {
77 8787 : if(n_==buf_.size()) {
78 600 : const double old=buf_[first_];
79 600 : buf_[first_]=v;
80 600 : first_=(first_+1)%buf_.size();
81 600 : if(max_<v) {
82 164 : max_=v;
83 436 : } else if(old==max_) {
84 100 : recomputeMax();
85 : }
86 : } else {
87 8187 : buf_[(first_+n_)%buf_.size()]=v;
88 8187 : n_++;
89 8187 : if(n_==1||max_<v) {
90 7387 : max_=v;
91 : }
92 : }
93 8787 : }
94 : /// replace the oldest value
95 : void setFirst(const double v) {
96 1062 : buf_[first_]=v;
97 1062 : recomputeMax();
98 288 : }
99 : };
100 :
101 : /**
102 : The Metainference algorithm shared by the METAINFERENCE bias and by the
103 : ISDB collective variables that inherit from MetainferenceBase.
104 :
105 : This is not an Action: it is owned by an Action, which it uses for input
106 : parsing, communication, log, restart and status files. The owner adds its
107 : own value (the bias or the score) and calls setup() once the number of data
108 : and the experimental values are known; setup() adds the components of
109 : Metainference to the owner. At every step the owner fills data() with the
110 : back-calculated values, calls calculate() and applies the derivatives of the
111 : energy with respect to the data (getMetaDer()) and to the bias used for
112 : reweighting (getBiasDer()).
113 : */
114 : class MetainferenceCore {
115 : public:
116 : /// get a pointer to a component of the owner action, given its name
117 : typedef std::function<Value*(const std::string&)> ComponentGetter;
118 : private:
119 : Action& action_;
120 : Communicator& comm;
121 : Communicator& multi_sim_comm;
122 : Log& log;
123 :
124 : // true if Metainference is active (always for the bias, DOSCORE for the CVs)
125 : bool active_;
126 : unsigned write_stride_;
127 : // number of experimental data
128 : unsigned narg;
129 : // experimental data
130 : std::vector<double> parameters;
131 : // metainference derivatives
132 : std::vector<double> metader_;
133 : // vector of back-calculated experimental data
134 : std::vector<double> calc_data_;
135 :
136 : // noise type
137 : unsigned noise_type_;
138 : enum { GAUSS, MGAUSS, OUTLIERS, MOUTLIERS, GENERIC };
139 : unsigned gen_likelihood_;
140 : enum { LIKE_GAUSS, LIKE_LOGN };
141 : bool doscale_;
142 : unsigned scale_prior_;
143 : enum { SC_GAUSS, SC_FLAT };
144 : double scale_;
145 : double scale_mu_;
146 : double scale_min_;
147 : double scale_max_;
148 : double Dscale_;
149 : // scale is data scaling factor
150 : // noise type
151 : unsigned offset_prior_;
152 : bool dooffset_;
153 : double offset_;
154 : double offset_mu_;
155 : double offset_min_;
156 : double offset_max_;
157 : double Doffset_;
158 : // scale and offset regression
159 : bool doregres_zero_;
160 : int nregres_zero_;
161 : // sigma is data uncertainty
162 : std::vector<double> sigma_;
163 : std::vector<double> sigma_min_;
164 : std::vector<double> sigma_max_;
165 : std::vector<double> Dsigma_;
166 : // sigma_mean is uncertainty in the mean estimate
167 : std::vector<double> sigma_mean2_;
168 : // this is the estimator of the mean value per replica for generic metainference
169 : std::vector<double> ftilde_;
170 : double Dftilde_;
171 :
172 : // temperature in kbt
173 : double kbt_;
174 :
175 : // Monte Carlo stuff
176 : std::vector<Random> random;
177 : unsigned MCsteps_;
178 : long long unsigned MCaccept_;
179 : long long unsigned MCacceptScale_;
180 : long long unsigned MCacceptFT_;
181 : long long unsigned MCtrial_;
182 : unsigned MCchunksize_;
183 :
184 : // output
185 : Value* valueScale=nullptr;
186 : Value* valueOffset=nullptr;
187 : Value* valueAccept=nullptr;
188 : Value* valueAcceptScale=nullptr;
189 : Value* valueAcceptFT=nullptr;
190 : Value* valueBiasDer_=nullptr;
191 : Value* valueWeight_=nullptr;
192 : Value* valueNeff_=nullptr;
193 : std::vector<Value*> valueSigma;
194 : std::vector<Value*> valueSigmaMean;
195 : std::vector<Value*> valueFtilde;
196 :
197 : // restart
198 : std::string status_file_name_;
199 : OFile sfile_;
200 : std::string fmt_;
201 :
202 : // others
203 : bool firstTime;
204 : std::vector<bool> firstTimeW;
205 : bool master;
206 : bool do_reweight_;
207 : unsigned do_optsigmamean_;
208 : unsigned nrep_;
209 : unsigned replica_;
210 :
211 : // selector
212 : unsigned nsel_;
213 : std::string selector_;
214 : unsigned iselect;
215 :
216 : // optimize sigma mean
217 : std::vector< std::vector <SigmaMeanHistory> > sigma_mean2_last_;
218 : unsigned optsigmamean_stride_;
219 : // optimize sigma max
220 : unsigned N_optimized_step_;
221 : unsigned optimized_step_;
222 : bool sigmamax_opt_done_;
223 : std::vector<double> sigma_max_est_;
224 : // GENERIC: sigma mean (the minimum of sigma) has exceeded SIGMA_MAX
225 : bool warned_sigma_max_=false;
226 :
227 : // average weights
228 : double decay_w_;
229 : std::vector< std::vector <double> > average_weights_;
230 :
231 : // derivative of the energy with respect to the bias used for reweighting
232 : double biasDer_;
233 : // value of the bias used for reweighting at this step
234 : double reweight_bias_;
235 : // the sigmas moved by the last MC move (MC_CHUNKSIZE)
236 : std::vector<unsigned> moved_;
237 : // the energy of the MC has been updated with differences
238 : bool energy_from_differences_=false;
239 :
240 : double getEnergyMIGEN(const std::vector<double> &mean, const std::vector<double> &ftilde, const std::vector<double> &sigma,
241 : const double scale, const double offset);
242 : double getEnergySP(const std::vector<double> &mean, const std::vector<double> &sigma,
243 : const double scale, const double offset);
244 : double getEnergySPE(const std::vector<double> &mean, const std::vector<double> &sigma,
245 : const double scale, const double offset);
246 : double getEnergyGJ(const std::vector<double> &mean, const std::vector<double> &sigma,
247 : const double scale, const double offset);
248 : double getEnergyGJE(const std::vector<double> &mean, const std::vector<double> &sigma,
249 : const double scale, const double offset);
250 : double getEnergyTerm(const unsigned i, const double mean_i, const double sigma_i) const;
251 : void getEnergyForceSP(const std::vector<double> &mean, const std::vector<double> &dmean_x, const std::vector<double> &dmean_b);
252 : void getEnergyForceSPE(const std::vector<double> &mean, const std::vector<double> &dmean_x, const std::vector<double> &dmean_b);
253 : void getEnergyForceGJ(const std::vector<double> &mean, const std::vector<double> &dmean_x, const std::vector<double> &dmean_b);
254 : void getEnergyForceGJE(const std::vector<double> &mean, const std::vector<double> &dmean_x, const std::vector<double> &dmean_b);
255 : void getEnergyForceMIGEN(const std::vector<double> &mean, const std::vector<double> &dmean_x, const std::vector<double> &dmean_b);
256 : void get_weights(double &weight, double &norm, double &neff);
257 : void replica_averaging(const double weight, const double norm, std::vector<double> &mean, std::vector<double> &dmean_b);
258 : void get_sigma_mean(const double weight, const double norm, const double neff, const std::vector<double> &mean);
259 : void do_regression_zero(const std::vector<double> &mean);
260 : void moveTilde(const std::vector<double> &mean_, double &old_energy);
261 : void moveScaleOffset(const std::vector<double> &mean_, double &old_energy);
262 : void moveSigmas(const std::vector<double> &mean_, double &old_energy, const unsigned i, const std::vector<unsigned> &indices, bool &breaknow);
263 : double doMonteCarlo(const std::vector<double> &mean);
264 :
265 : public:
266 : /// register the keywords and the components shared by all the Metainference actions
267 : static void registerKeywords( Keywords& keys );
268 : /// read the input of the owner action; active is false for a CV without DOSCORE
269 : MetainferenceCore(Action& action, bool active);
270 : ~MetainferenceCore();
271 : MetainferenceCore(const MetainferenceCore&) = delete;
272 : MetainferenceCore& operator=(const MetainferenceCore&) = delete;
273 : /// set the experimental values, read the status file at restart and add the components to the owner
274 : void setup(const std::vector<double>& params, ActionWithValue& av, const ComponentGetter& getComponent);
275 : /// the back-calculated data, to be set before calculate()
276 : std::vector<double>& data() {
277 : return calc_data_;
278 : }
279 : /// run Metainference with the current data and return the energy summed over the replicas
280 : /// isel is the selector index, reweight_bias the value of the bias used with REWEIGHT
281 : double calculate(unsigned isel, double reweight_bias);
282 : /// derivative of the energy with respect to the back-calculated datum i
283 : double getMetaDer(const unsigned i) const {
284 17048 : return metader_[i];
285 : }
286 : /// derivative of the energy with respect to the bias used with REWEIGHT
287 : double getBiasDer() const {
288 2084 : return biasDer_;
289 : }
290 : /// the selector index (0 without SELECTOR)
291 : unsigned getSelector() const;
292 : bool doReweight() const {
293 6373 : return do_reweight_;
294 : }
295 : bool isActive() const {
296 : return active_;
297 : }
298 : unsigned getWriteStride() const {
299 : return write_stride_;
300 : }
301 : unsigned getNarg() const {
302 : return narg;
303 : }
304 : /// write the status file if this is a step at which it should be written
305 : void update();
306 : void writeStatus();
307 : };
308 :
309 : }
310 : }
311 :
312 : #endif
|
