Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2017-2023 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_MetainferenceBase_h
23 : #define __PLUMED_isdb_MetainferenceBase_h
24 :
25 : #include "core/ActionWithValue.h"
26 : #include "core/ActionAtomistic.h"
27 : #include "core/ActionWithArguments.h"
28 : #include "tools/Communicator.h"
29 : #include "core/PlumedMain.h"
30 : #include "tools/Random.h"
31 : #include "tools/OpenMP.h"
32 :
33 : #define PLUMED_METAINF_INIT(ao) Action(ao),MetainferenceBase(ao)
34 :
35 : namespace PLMD {
36 : namespace isdb {
37 :
38 : /**
39 : \ingroup INHERIT
40 : This is the abstract base class to use for implementing new ISDB Metainference actions, within it there is
41 : information as to how to go about implementing a new Metainference action.
42 : */
43 :
44 : class MetainferenceBase :
45 : public ActionAtomistic,
46 : public ActionWithArguments,
47 : public ActionWithValue
48 : {
49 : private:
50 : std::vector<double> forces;
51 : std::vector<double> forcesToApply;
52 :
53 : // activate metainference
54 : bool doscore_;
55 : unsigned write_stride_;
56 : // number of experimental data
57 : unsigned narg;
58 : // experimental data
59 : std::vector<double> parameters;
60 : // metainference derivatives
61 : std::vector<double> metader_;
62 : // vector of back-calculated experimental data
63 : std::vector<double> calc_data_;
64 :
65 : // noise type
66 : unsigned noise_type_;
67 : enum { GAUSS, MGAUSS, OUTLIERS, MOUTLIERS, GENERIC };
68 : unsigned gen_likelihood_;
69 : enum { LIKE_GAUSS, LIKE_LOGN };
70 : bool doscale_;
71 : unsigned scale_prior_;
72 : enum { SC_GAUSS, SC_FLAT };
73 : double scale_;
74 : double scale_mu_;
75 : double scale_min_;
76 : double scale_max_;
77 : double Dscale_;
78 : // scale is data scaling factor
79 : // noise type
80 : unsigned offset_prior_;
81 : bool dooffset_;
82 : double offset_;
83 : double offset_mu_;
84 : double offset_min_;
85 : double offset_max_;
86 : double Doffset_;
87 : // scale and offset regression
88 : bool doregres_zero_;
89 : int nregres_zero_;
90 : // sigma is data uncertainty
91 : std::vector<double> sigma_;
92 : std::vector<double> sigma_min_;
93 : std::vector<double> sigma_max_;
94 : std::vector<double> Dsigma_;
95 : // sigma_mean is uncertainty in the mean estimate
96 : std::vector<double> sigma_mean2_;
97 : // this is the estimator of the mean value per replica for generic metainference
98 : std::vector<double> ftilde_;
99 : double Dftilde_;
100 :
101 : // temperature in kbt
102 : double kbt_;
103 :
104 : // Monte Carlo stuff
105 : std::vector<Random> random;
106 : unsigned MCsteps_;
107 : long long unsigned MCaccept_;
108 : long long unsigned MCacceptScale_;
109 : long long unsigned MCacceptFT_;
110 : long long unsigned MCtrial_;
111 : unsigned MCchunksize_;
112 :
113 : // output
114 : Value* valueScore;
115 : Value* valueScale;
116 : Value* valueOffset;
117 : Value* valueAccept;
118 : Value* valueAcceptScale;
119 : Value* valueAcceptFT;
120 : std::vector<Value*> valueSigma;
121 : std::vector<Value*> valueSigmaMean;
122 : std::vector<Value*> valueFtilde;
123 :
124 : // restart
125 : std::string status_file_name_;
126 : OFile sfile_;
127 :
128 : // others
129 : bool firstTime;
130 : std::vector<bool> firstTimeW;
131 : bool master;
132 : bool do_reweight_;
133 : unsigned do_optsigmamean_;
134 : unsigned nrep_;
135 : unsigned replica_;
136 :
137 : // selector
138 : unsigned nsel_;
139 : std::string selector_;
140 : unsigned iselect;
141 :
142 : // optimize sigma mean
143 : std::vector< std::vector < std::vector <double> > > sigma_mean2_last_;
144 : unsigned optsigmamean_stride_;
145 : // optimize sigma max
146 : unsigned N_optimized_step_;
147 : unsigned optimized_step_;
148 : bool sigmamax_opt_done_;
149 : std::vector<double> sigma_max_est_;
150 :
151 : // average weights
152 : double decay_w_;
153 : std::vector< std::vector <double> > average_weights_;
154 :
155 : double getEnergyMIGEN(const std::vector<double> &mean, const std::vector<double> &ftilde, const std::vector<double> &sigma,
156 : const double scale, const double offset);
157 : double getEnergySP(const std::vector<double> &mean, const std::vector<double> &sigma,
158 : const double scale, const double offset);
159 : double getEnergySPE(const std::vector<double> &mean, const std::vector<double> &sigma,
160 : const double scale, const double offset);
161 : double getEnergyGJ(const std::vector<double> &mean, const std::vector<double> &sigma,
162 : const double scale, const double offset);
163 : double getEnergyGJE(const std::vector<double> &mean, const std::vector<double> &sigma,
164 : const double scale, const double offset);
165 : void setMetaDer(const unsigned index, const double der);
166 : void getEnergyForceSP(const std::vector<double> &mean, const std::vector<double> &dmean_x, const std::vector<double> &dmean_b);
167 : void getEnergyForceSPE(const std::vector<double> &mean, const std::vector<double> &dmean_x, const std::vector<double> &dmean_b);
168 : void getEnergyForceGJ(const std::vector<double> &mean, const std::vector<double> &dmean_x, const std::vector<double> &dmean_b);
169 : void getEnergyForceGJE(const std::vector<double> &mean, const std::vector<double> &dmean_x, const std::vector<double> &dmean_b);
170 : void getEnergyForceMIGEN(const std::vector<double> &mean, const std::vector<double> &dmean_x, const std::vector<double> &dmean_b);
171 : double getCalcData(const unsigned index);
172 : void get_weights(double &weight, double &norm, double &neff);
173 : void replica_averaging(const double weight, const double norm, std::vector<double> &mean, std::vector<double> &dmean_b);
174 : void get_sigma_mean(const double weight, const double norm, const double neff, const std::vector<double> &mean);
175 : void do_regression_zero(const std::vector<double> &mean);
176 : void moveTilde(const std::vector<double> &mean_, double &old_energy);
177 : void moveScaleOffset(const std::vector<double> &mean_, double &old_energy);
178 : void moveSigmas(const std::vector<double> &mean_, double &old_energy, const unsigned i, const std::vector<unsigned> &indices, bool &breaknow);
179 : double doMonteCarlo(const std::vector<double> &mean);
180 :
181 : public:
182 : static void registerKeywords( Keywords& keys );
183 : explicit MetainferenceBase(const ActionOptions&);
184 : ~MetainferenceBase();
185 : void Initialise(const unsigned input);
186 : void Selector();
187 : unsigned getNarg();
188 : void setNarg(const unsigned input);
189 : void setParameters(const std::vector<double>& input);
190 : void setParameter(const double input);
191 : void setCalcData(const unsigned index, const double datum);
192 : void setCalcData(const std::vector<double>& data);
193 : bool getDoScore();
194 : unsigned getWstride();
195 : double getScore();
196 : void setScore(const double score);
197 : void setDerivatives();
198 : double getMetaDer(const unsigned index);
199 : void writeStatus();
200 : void turnOnDerivatives() override;
201 : unsigned getNumberOfDerivatives() override;
202 : void lockRequests() override;
203 : void unlockRequests() override;
204 : void calculateNumericalDerivatives( ActionWithValue* a ) override;
205 : void apply() override;
206 : void setArgDerivatives(Value *v, const double &d);
207 : void setAtomsDerivatives(Value*v, const unsigned i, const Vector&d);
208 : void setBoxDerivatives(Value*v, const Tensor&d);
209 : };
210 :
211 : inline
212 : void MetainferenceBase::setNarg(const unsigned input)
213 : {
214 31 : narg = input;
215 : }
216 :
217 : inline
218 : bool MetainferenceBase::getDoScore()
219 : {
220 22634 : return doscore_;
221 : }
222 :
223 : inline
224 : unsigned MetainferenceBase::getWstride()
225 : {
226 1434 : return write_stride_;
227 : }
228 :
229 : inline
230 : unsigned MetainferenceBase::getNarg()
231 : {
232 2225 : return narg;
233 : }
234 :
235 : inline
236 : void MetainferenceBase::setMetaDer(const unsigned index, const double der)
237 : {
238 7146 : metader_[index] = der;
239 : }
240 :
241 : inline
242 : double MetainferenceBase::getMetaDer(const unsigned index)
243 : {
244 909788 : return metader_[index];
245 : }
246 :
247 : inline
248 : double MetainferenceBase::getCalcData(const unsigned index)
249 : {
250 : return calc_data_[index];
251 : }
252 :
253 : inline
254 : void MetainferenceBase::setCalcData(const unsigned index, const double datum)
255 : {
256 3868 : calc_data_[index] = datum;
257 3868 : }
258 :
259 : inline
260 : void MetainferenceBase::setCalcData(const std::vector<double>& data)
261 : {
262 : for(unsigned i=0; i<data.size(); i++) calc_data_[i] = data[i];
263 : }
264 :
265 : inline
266 27 : void MetainferenceBase::setParameters(const std::vector<double>& input) {
267 168 : for(unsigned i=0; i<input.size(); i++) parameters.push_back(input[i]);
268 27 : }
269 :
270 : inline
271 : void MetainferenceBase::setParameter(const double input) {
272 2356 : parameters.push_back(input);
273 2356 : }
274 :
275 : inline
276 2225 : void MetainferenceBase::setScore(const double score) {
277 2225 : valueScore->set(score);
278 2225 : }
279 :
280 : inline
281 112 : void MetainferenceBase::setDerivatives() {
282 : // Get appropriate number of derivatives
283 : // Derivatives are first for arguments and then for atoms
284 : unsigned nder;
285 112 : if( getNumberOfAtoms()>0 ) {
286 112 : nder = 3*getNumberOfAtoms() + 9 + getNumberOfArguments();
287 : } else {
288 0 : nder = getNumberOfArguments();
289 : }
290 :
291 : // Resize all derivative arrays
292 112 : forces.resize( nder ); forcesToApply.resize( nder );
293 21460 : for(int i=0; i<getNumberOfComponents(); ++i) getPntrToComponent(i)->resizeDerivatives(nder);
294 112 : }
295 :
296 : inline
297 2747 : void MetainferenceBase::turnOnDerivatives() {
298 2747 : ActionWithValue::turnOnDerivatives();
299 2747 : }
300 :
301 : inline
302 2828290 : unsigned MetainferenceBase::getNumberOfDerivatives() {
303 2828290 : if( getNumberOfAtoms()>0 ) {
304 2828290 : return 3*getNumberOfAtoms() + 9 + getNumberOfArguments();
305 : }
306 0 : return getNumberOfArguments();
307 : }
308 :
309 : inline
310 717 : void MetainferenceBase::lockRequests() {
311 : ActionAtomistic::lockRequests();
312 : ActionWithArguments::lockRequests();
313 717 : }
314 :
315 : inline
316 717 : void MetainferenceBase::unlockRequests() {
317 : ActionAtomistic::unlockRequests();
318 : ActionWithArguments::unlockRequests();
319 717 : }
320 :
321 : inline
322 75 : void MetainferenceBase::calculateNumericalDerivatives( ActionWithValue* a=NULL ) {
323 75 : if( getNumberOfArguments()>0 ) {
324 48 : ActionWithArguments::calculateNumericalDerivatives( a );
325 : }
326 75 : if( getNumberOfAtoms()>0 ) {
327 75 : Matrix<double> save_derivatives( getNumberOfComponents(), getNumberOfArguments() );
328 1293 : for(int j=0; j<getNumberOfComponents(); ++j) {
329 2322 : for(unsigned i=0; i<getNumberOfArguments(); ++i) if(getPntrToComponent(j)->hasDerivatives()) save_derivatives(j,i)=getPntrToComponent(j)->getDerivative(i);
330 : }
331 75 : calculateAtomicNumericalDerivatives( a, getNumberOfArguments() );
332 1293 : for(int j=0; j<getNumberOfComponents(); ++j) {
333 2322 : for(unsigned i=0; i<getNumberOfArguments(); ++i) if(getPntrToComponent(j)->hasDerivatives()) getPntrToComponent(j)->addDerivative( i, save_derivatives(j,i) );
334 : }
335 : }
336 75 : }
337 :
338 : inline
339 717 : void MetainferenceBase::apply() {
340 717 : bool wasforced=false; forcesToApply.assign(forcesToApply.size(),0.0);
341 32982 : for(int i=0; i<getNumberOfComponents(); ++i) {
342 32265 : if( getPntrToComponent(i)->applyForce( forces ) ) {
343 : wasforced=true;
344 41664408 : for(unsigned i=0; i<forces.size(); ++i) forcesToApply[i]+=forces[i];
345 : }
346 : }
347 717 : if( wasforced ) {
348 350 : unsigned ind=0; addForcesOnArguments( 0, forcesToApply, ind, getLabel() );
349 350 : if( getNumberOfAtoms()>0 ) setForcesOnAtoms( forcesToApply, ind );
350 : }
351 717 : }
352 :
353 : inline
354 : void MetainferenceBase::setArgDerivatives(Value *v, const double &d) {
355 160 : v->addDerivative(0,d);
356 : }
357 :
358 : inline
359 2771792 : void MetainferenceBase::setAtomsDerivatives(Value*v, const unsigned i, const Vector&d) {
360 2771792 : const unsigned noa=getNumberOfArguments();
361 2771792 : v->addDerivative(noa+3*i+0,d[0]);
362 2771792 : v->addDerivative(noa+3*i+1,d[1]);
363 2771792 : v->addDerivative(noa+3*i+2,d[2]);
364 2771792 : }
365 :
366 : inline
367 12016 : void MetainferenceBase::setBoxDerivatives(Value* v,const Tensor&d) {
368 12016 : const unsigned noa=getNumberOfArguments();
369 : const unsigned nat=getNumberOfAtoms();
370 12016 : v->addDerivative(noa+3*nat+0,d(0,0));
371 12016 : v->addDerivative(noa+3*nat+1,d(0,1));
372 12016 : v->addDerivative(noa+3*nat+2,d(0,2));
373 12016 : v->addDerivative(noa+3*nat+3,d(1,0));
374 12016 : v->addDerivative(noa+3*nat+4,d(1,1));
375 12016 : v->addDerivative(noa+3*nat+5,d(1,2));
376 12016 : v->addDerivative(noa+3*nat+6,d(2,0));
377 12016 : v->addDerivative(noa+3*nat+7,d(2,1));
378 12016 : v->addDerivative(noa+3*nat+8,d(2,2));
379 12016 : }
380 :
381 :
382 : }
383 : }
384 :
385 : #endif
386 :
|
