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