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