Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2021 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 "BasisFunctions.h"
24 :
25 : #include "core/ActionRegister.h"
26 : #include "lepton/Lepton.h"
27 :
28 :
29 : namespace PLMD {
30 : namespace ves {
31 :
32 : //+PLUMEDOC VES_BASISF BF_CUSTOM
33 : /*
34 : Basis functions given by arbitrary mathematical expressions.
35 :
36 : This allows you to define basis functions using arbitrary mathematical expressions
37 : that are parsed using the lepton library.
38 : The basis functions
39 : $f_{i}(x)$ are given in mathematical expressions with _x_ as a variable using
40 : the numbered FUNC keywords that start from
41 : FUNC1. Consistent with other basis functions is $f_{0}(x)=1$ defined as
42 : the constant. The interval on which the basis functions are defined is
43 : given using the MINIMUM and MAXIMUM keywords.
44 :
45 : Using the TRANSFORM keyword it is possible to define a function $x(t)$ that
46 : is used to transform the argument before calculating the basis functions
47 : values. The variables _min_ and _max_ can be used to indicate the minimum
48 : and the maximum of the interval. By default the arguments are not transformed,
49 : i.e. $x(t)=t$.
50 :
51 : For periodic basis functions you should use the PERIODIC flag to indicate
52 : that they are periodic.
53 :
54 : The basis functions $f_{i}(x)$ and the transform function $x(t)$ need
55 : to be well behaved in the interval on which the basis functions are defined,
56 : e.g. not result in a not a number (nan) or infinity (inf).
57 : The code will not perform checks to make sure that this is the case unless the
58 : flag CHECK_NAN_INF is enabled.
59 :
60 : ## Examples
61 :
62 : Defining Legendre polynomial basis functions of order 6 using BF_CUSTOM
63 : where the appropriate transform function is given by the TRANSFORM keyword.
64 : This is just an example of what can be done, in practice you should use
65 : [BF_LEGENDRE](BF_LEGENDRE.md) for Legendre polynomial basis functions.
66 :
67 : ```plumed
68 : BF_CUSTOM ...
69 : TRANSFORM=(t-(min+max)/2)/((max-min)/2)
70 : FUNC1=x
71 : FUNC2=(1/2)*(3*x^2-1)
72 : FUNC3=(1/2)*(5*x^3-3*x)
73 : FUNC4=(1/8)*(35*x^4-30*x^2+3)
74 : FUNC5=(1/8)*(63*x^5-70*x^3+15*x)
75 : FUNC6=(1/16)*(231*x^6-315*x^4+105*x^2-5)
76 : MINIMUM=-4.0
77 : MAXIMUM=4.0
78 : LABEL=bf1
79 : ... BF_CUSTOM
80 : ```
81 :
82 : Defining Fourier basis functions of order 3 using BF_CUSTOM where the
83 : periodicity is indicated using the PERIODIC flag. This is just an example
84 : of what can be done, in practice you should use [BF_FOURIER](BF_FOURIER.md)
85 : for Fourier basis functions.
86 :
87 : ```plumed
88 : BF_CUSTOM ...
89 : FUNC1=cos(x)
90 : FUNC2=sin(x)
91 : FUNC3=cos(2*x)
92 : FUNC4=sin(2*x)
93 : FUNC5=cos(3*x)
94 : FUNC6=sin(3*x)
95 : MINIMUM=-pi
96 : MAXIMUM=+pi
97 : LABEL=bf1
98 : PERIODIC
99 : ... BF_CUSTOM
100 : ```
101 :
102 :
103 : */
104 : //+ENDPLUMEDOC
105 :
106 : class BF_Custom : public BasisFunctions {
107 : private:
108 : lepton::CompiledExpression transf_value_expression_;
109 : lepton::CompiledExpression transf_deriv_expression_;
110 : double* transf_value_lepton_ref_;
111 : double* transf_deriv_lepton_ref_;
112 : std::vector<lepton::CompiledExpression> bf_values_expressions_;
113 : std::vector<lepton::CompiledExpression> bf_derivs_expressions_;
114 : std::vector<double*> bf_values_lepton_ref_;
115 : std::vector<double*> bf_derivs_lepton_ref_;
116 : std::string variable_str_;
117 : std::string transf_variable_str_;
118 : bool do_transf_;
119 : bool check_nan_inf_;
120 : public:
121 : static void registerKeywords( Keywords&);
122 : explicit BF_Custom(const ActionOptions&);
123 30 : ~BF_Custom() {};
124 : void getAllValues(const double, double&, bool&, std::vector<double>&, std::vector<double>&) const override;
125 : };
126 :
127 : PLUMED_REGISTER_ACTION(BF_Custom,"BF_CUSTOM")
128 :
129 12 : void BF_Custom::registerKeywords(Keywords& keys) {
130 12 : BasisFunctions::registerKeywords(keys);
131 24 : keys.remove("ORDER");
132 12 : keys.add("numbered","FUNC","The basis functions f_i(x) given in mathematical expressions using _x_ as a variable.");
133 12 : keys.add("optional","TRANSFORM","An optional function that can be used to transform the argument before calculating the basis function values. You should use _t_ as a variable. You can use the variables _min_ and _max_ to give the minimum and the maximum of the interval.");
134 12 : keys.addFlag("PERIODIC",false,"Indicate that the basis functions are periodic.");
135 12 : keys.addFlag("CHECK_NAN_INF",false,"Check that the basis functions do not result in a not a number (nan) or infinity (inf).");
136 12 : keys.remove("NUMERICAL_INTEGRALS");
137 12 : }
138 :
139 :
140 10 : BF_Custom::BF_Custom(const ActionOptions&ao):
141 : PLUMED_VES_BASISFUNCTIONS_INIT(ao),
142 10 : transf_value_lepton_ref_(nullptr),
143 10 : transf_deriv_lepton_ref_(nullptr),
144 10 : bf_values_expressions_(0),
145 10 : bf_derivs_expressions_(0),
146 10 : bf_values_lepton_ref_(0,nullptr),
147 10 : bf_derivs_lepton_ref_(0,nullptr),
148 10 : variable_str_("x"),
149 10 : transf_variable_str_("t"),
150 10 : do_transf_(false),
151 20 : check_nan_inf_(false) {
152 : std::vector<std::string> bf_str;
153 10 : std::string str_t1="1";
154 10 : bf_str.push_back(str_t1);
155 10 : for(int i=1;; i++) {
156 : std::string str_t2;
157 112 : if(!parseNumbered("FUNC",i,str_t2)) {
158 : break;
159 : }
160 : std::string is;
161 46 : Tools::convert(i,is);
162 92 : addKeywordToList("FUNC"+is,str_t2);
163 46 : bf_str.push_back(str_t2);
164 46 : }
165 : //
166 10 : if(bf_str.size()==1) {
167 0 : plumed_merror(getName()+" with label "+getLabel()+": No FUNC keywords given");
168 : }
169 :
170 10 : setOrder(bf_str.size()-1);
171 10 : setNumberOfBasisFunctions(getOrder()+1);
172 10 : setIntrinsicInterval(intervalMin(),intervalMax());
173 10 : bool periodic = false;
174 10 : parseFlag("PERIODIC",periodic);
175 10 : addKeywordToList("PERIODIC",periodic);
176 10 : if(periodic) {
177 : setPeriodic();
178 : } else {
179 : setNonPeriodic();
180 : }
181 : setIntervalBounded();
182 10 : setType("custom_functions");
183 30 : setDescription("Custom Functions");
184 : //
185 10 : std::vector<std::string> bf_values_parsed(getNumberOfBasisFunctions());
186 10 : std::vector<std::string> bf_derivs_parsed(getNumberOfBasisFunctions());
187 : bf_values_parsed[0] = "1";
188 : bf_derivs_parsed[0] = "0";
189 : //
190 10 : bf_values_expressions_.resize(getNumberOfBasisFunctions());
191 10 : bf_derivs_expressions_.resize(getNumberOfBasisFunctions());
192 10 : bf_values_lepton_ref_.resize(getNumberOfBasisFunctions());
193 10 : bf_derivs_lepton_ref_.resize(getNumberOfBasisFunctions());
194 : //
195 56 : for(unsigned int i=1; i<getNumberOfBasisFunctions(); i++) {
196 : std::string is;
197 46 : Tools::convert(i,is);
198 : try {
199 46 : lepton::ParsedExpression pe_value = lepton::Parser::parse(bf_str[i]).optimize(lepton::Constants());
200 46 : std::ostringstream tmp_stream;
201 46 : tmp_stream << pe_value;
202 46 : bf_values_parsed[i] = tmp_stream.str();
203 46 : bf_values_expressions_[i] = pe_value.createCompiledExpression();
204 46 : } catch(PLMD::lepton::Exception& exc) {
205 0 : plumed_merror("There was some problem in parsing the function "+bf_str[i]+" given in FUNC"+is + " with lepton");
206 0 : }
207 :
208 : std::vector<std::string> var_str;
209 92 : for(auto &p: bf_values_expressions_[i].getVariables()) {
210 46 : var_str.push_back(p);
211 : }
212 46 : if(var_str.size()!=1) {
213 0 : plumed_merror("Problem with function "+bf_str[i]+" given in FUNC"+is+": there should only be one variable");
214 : }
215 46 : if(var_str[0]!=variable_str_) {
216 0 : plumed_merror("Problem with function "+bf_str[i]+" given in FUNC"+is+": you should use "+variable_str_+" as a variable");
217 : }
218 :
219 : try {
220 92 : lepton::ParsedExpression pe_deriv = lepton::Parser::parse(bf_str[i]).differentiate(variable_str_).optimize(lepton::Constants());
221 46 : std::ostringstream tmp_stream2;
222 46 : tmp_stream2 << pe_deriv;
223 46 : bf_derivs_parsed[i] = tmp_stream2.str();
224 46 : bf_derivs_expressions_[i] = pe_deriv.createCompiledExpression();
225 46 : } catch(PLMD::lepton::Exception& exc) {
226 0 : plumed_merror("There was some problem in parsing the derivative of the function "+bf_str[i]+" given in FUNC"+is + " with lepton");
227 0 : }
228 :
229 : try {
230 46 : bf_values_lepton_ref_[i] = &bf_values_expressions_[i].getVariableReference(variable_str_);
231 0 : } catch(PLMD::lepton::Exception& exc) {}
232 :
233 : try {
234 46 : bf_derivs_lepton_ref_[i] = &bf_derivs_expressions_[i].getVariableReference(variable_str_);
235 8 : } catch(PLMD::lepton::Exception& exc) {}
236 :
237 46 : }
238 :
239 : std::string transf_value_parsed;
240 : std::string transf_deriv_parsed;
241 : std::string transf_str;
242 20 : parse("TRANSFORM",transf_str);
243 10 : if(transf_str.size()>0) {
244 6 : do_transf_ = true;
245 12 : addKeywordToList("TRANSFORM",transf_str);
246 : for(unsigned int k=0;; k++) {
247 10 : if(transf_str.find("min")!=std::string::npos) {
248 8 : transf_str.replace(transf_str.find("min"), std::string("min").length(),intervalMinStr());
249 : } else {
250 : break;
251 : }
252 : }
253 : for(unsigned int k=0;; k++) {
254 10 : if(transf_str.find("max")!=std::string::npos) {
255 8 : transf_str.replace(transf_str.find("max"), std::string("max").length(),intervalMaxStr());
256 : } else {
257 : break;
258 : }
259 : }
260 :
261 : try {
262 6 : lepton::ParsedExpression pe_value = lepton::Parser::parse(transf_str).optimize(lepton::Constants());;
263 6 : std::ostringstream tmp_stream;
264 6 : tmp_stream << pe_value;
265 6 : transf_value_parsed = tmp_stream.str();
266 6 : transf_value_expression_ = pe_value.createCompiledExpression();
267 6 : } catch(PLMD::lepton::Exception& exc) {
268 0 : plumed_merror("There was some problem in parsing the function "+transf_str+" given in TRANSFORM with lepton");
269 0 : }
270 :
271 : std::vector<std::string> var_str;
272 12 : for(auto &p: transf_value_expression_.getVariables()) {
273 6 : var_str.push_back(p);
274 : }
275 6 : if(var_str.size()!=1) {
276 0 : plumed_merror("Problem with function "+transf_str+" given in TRANSFORM: there should only be one variable");
277 : }
278 6 : if(var_str[0]!=transf_variable_str_) {
279 0 : plumed_merror("Problem with function "+transf_str+" given in TRANSFORM: you should use "+transf_variable_str_+" as a variable");
280 : }
281 :
282 : try {
283 12 : lepton::ParsedExpression pe_deriv = lepton::Parser::parse(transf_str).differentiate(transf_variable_str_).optimize(lepton::Constants());;
284 6 : std::ostringstream tmp_stream2;
285 6 : tmp_stream2 << pe_deriv;
286 6 : transf_deriv_parsed = tmp_stream2.str();
287 6 : transf_deriv_expression_ = pe_deriv.createCompiledExpression();
288 6 : } catch(PLMD::lepton::Exception& exc) {
289 0 : plumed_merror("There was some problem in parsing the derivative of the function "+transf_str+" given in TRANSFORM with lepton");
290 0 : }
291 :
292 : try {
293 6 : transf_value_lepton_ref_ = &transf_value_expression_.getVariableReference(transf_variable_str_);
294 0 : } catch(PLMD::lepton::Exception& exc) {}
295 :
296 : try {
297 6 : transf_deriv_lepton_ref_ = &transf_deriv_expression_.getVariableReference(transf_variable_str_);
298 3 : } catch(PLMD::lepton::Exception& exc) {}
299 :
300 6 : }
301 : //
302 10 : log.printf(" Using the following functions [lepton parsed function and derivative]:\n");
303 66 : for(unsigned int i=0; i<getNumberOfBasisFunctions(); i++) {
304 56 : log.printf(" %u: %s [ %s | %s ] \n",i,bf_str[i].c_str(),bf_values_parsed[i].c_str(),bf_derivs_parsed[i].c_str());
305 :
306 : }
307 : //
308 10 : if(do_transf_) {
309 6 : log.printf(" Arguments are transformed using the following function [lepton parsed function and derivative]:\n");
310 6 : log.printf(" %s [ %s | %s ] \n",transf_str.c_str(),transf_value_parsed.c_str(),transf_deriv_parsed.c_str());
311 : } else {
312 4 : log.printf(" Arguments are not transformed\n");
313 : }
314 : //
315 :
316 10 : parseFlag("CHECK_NAN_INF",check_nan_inf_);
317 10 : addKeywordToList("CHECK_NAN_INF",check_nan_inf_);
318 10 : if(check_nan_inf_) {
319 6 : log.printf(" The code will check that values given are numercially stable, e.g. do not result in a not a number (nan) or infinity (inf).\n");
320 : } else {
321 4 : log.printf(" The code will NOT check that values given are numercially stable, e.g. do not result in a not a number (nan) or infinity (inf).\n");
322 : }
323 :
324 10 : setupBF();
325 10 : checkRead();
326 20 : }
327 :
328 :
329 24204 : void BF_Custom::getAllValues(const double arg, double& argT, bool& inside_range, std::vector<double>& values, std::vector<double>& derivs) const {
330 24204 : inside_range=true;
331 24204 : argT=checkIfArgumentInsideInterval(arg,inside_range);
332 24204 : double transf_derivf=1.0;
333 : //
334 24204 : if(do_transf_) {
335 :
336 14062 : if(transf_value_lepton_ref_) {
337 14062 : *transf_value_lepton_ref_ = argT;
338 : }
339 14062 : if(transf_deriv_lepton_ref_) {
340 6125 : *transf_deriv_lepton_ref_ = argT;
341 : }
342 :
343 14062 : argT = transf_value_expression_.evaluate();
344 14062 : transf_derivf = transf_deriv_expression_.evaluate();
345 :
346 14062 : if(check_nan_inf_ && (std::isnan(argT) || std::isinf(argT)) ) {
347 : std::string vs;
348 0 : Tools::convert(argT,vs);
349 0 : plumed_merror(getName()+" with label "+getLabel()+": problem with the transform function, it gives " + vs);
350 : }
351 :
352 14062 : if(check_nan_inf_ && (std::isnan(transf_derivf) || std::isinf(transf_derivf)) ) {
353 : std::string vs;
354 0 : Tools::convert(transf_derivf,vs);
355 0 : plumed_merror(getName()+" with label "+getLabel()+": problem with the transform function, its derivative gives " + vs);
356 : }
357 : }
358 : //
359 24204 : values[0]=1.0;
360 24204 : derivs[0]=0.0;
361 137488 : for(unsigned int i=1; i < getNumberOfBasisFunctions(); i++) {
362 :
363 113284 : if(bf_values_lepton_ref_[i]) {
364 113284 : *bf_values_lepton_ref_[i] = argT;
365 : }
366 113284 : if(bf_derivs_lepton_ref_[i]) {
367 93594 : *bf_derivs_lepton_ref_[i] = argT;
368 : }
369 :
370 113284 : values[i] = bf_values_expressions_[i].evaluate();
371 113284 : derivs[i] = bf_derivs_expressions_[i].evaluate();
372 :
373 113284 : if(do_transf_) {
374 68306 : derivs[i]*=transf_derivf;
375 : }
376 : // NaN checks
377 113284 : if(check_nan_inf_ && (std::isnan(values[i]) || std::isinf(values[i])) ) {
378 : std::string vs;
379 0 : Tools::convert(values[i],vs);
380 : std::string is;
381 0 : Tools::convert(i,is);
382 0 : plumed_merror(getName()+" with label "+getLabel()+": problem with the basis function given in FUNC"+is+", it gives "+vs);
383 : }
384 : //
385 113284 : if(check_nan_inf_ && (std::isnan(derivs[i])|| std::isinf(derivs[i])) ) {
386 : std::string vs;
387 0 : Tools::convert(derivs[i],vs);
388 : std::string is;
389 0 : Tools::convert(i,is);
390 0 : plumed_merror(getName()+" with label "+getLabel()+": problem with derivative of the basis function given in FUNC"+is+", it gives "+vs);
391 : }
392 : }
393 24204 : if(!inside_range) {
394 2014 : for(unsigned int i=0; i<derivs.size(); i++) {
395 1705 : derivs[i]=0.0;
396 : }
397 : }
398 24204 : }
399 :
400 :
401 :
402 :
403 : }
404 : }
|
