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