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