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 12 : keys.remove("ORDER");
131 24 : keys.add("numbered","FUNC","The basis functions f_i(x) given in mathematical expressions using _x_ as a variable.");
132 24 : 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 24 : keys.addFlag("PERIODIC",false,"Indicate that the basis functions are periodic.");
134 24 : 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 : {
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)) {break;}
158 46 : std::string is; Tools::convert(i,is);
159 92 : addKeywordToList("FUNC"+is,str_t2);
160 46 : bf_str.push_back(str_t2);
161 46 : }
162 : //
163 10 : if(bf_str.size()==1) {plumed_merror(getName()+" with label "+getLabel()+": No FUNC keywords given");}
164 :
165 10 : setOrder(bf_str.size()-1);
166 10 : setNumberOfBasisFunctions(getOrder()+1);
167 10 : setIntrinsicInterval(intervalMin(),intervalMax());
168 10 : bool periodic = false;
169 20 : parseFlag("PERIODIC",periodic); addKeywordToList("PERIODIC",periodic);
170 10 : if(periodic) {setPeriodic();}
171 : else {setNonPeriodic();}
172 : setIntervalBounded();
173 10 : setType("custom_functions");
174 20 : setDescription("Custom Functions");
175 : //
176 10 : std::vector<std::string> bf_values_parsed(getNumberOfBasisFunctions());
177 10 : std::vector<std::string> bf_derivs_parsed(getNumberOfBasisFunctions());
178 : bf_values_parsed[0] = "1";
179 : bf_derivs_parsed[0] = "0";
180 : //
181 10 : bf_values_expressions_.resize(getNumberOfBasisFunctions());
182 10 : bf_derivs_expressions_.resize(getNumberOfBasisFunctions());
183 10 : bf_values_lepton_ref_.resize(getNumberOfBasisFunctions());
184 10 : bf_derivs_lepton_ref_.resize(getNumberOfBasisFunctions());
185 : //
186 56 : for(unsigned int i=1; i<getNumberOfBasisFunctions(); i++) {
187 46 : std::string is; Tools::convert(i,is);
188 : try {
189 46 : lepton::ParsedExpression pe_value = lepton::Parser::parse(bf_str[i]).optimize(lepton::Constants());
190 46 : std::ostringstream tmp_stream; tmp_stream << pe_value;
191 46 : bf_values_parsed[i] = tmp_stream.str();
192 46 : bf_values_expressions_[i] = pe_value.createCompiledExpression();
193 46 : }
194 0 : catch(PLMD::lepton::Exception& exc) {
195 0 : plumed_merror("There was some problem in parsing the function "+bf_str[i]+" given in FUNC"+is + " with lepton");
196 0 : }
197 :
198 : std::vector<std::string> var_str;
199 92 : for(auto &p: bf_values_expressions_[i].getVariables()) {
200 46 : var_str.push_back(p);
201 : }
202 46 : if(var_str.size()!=1) {
203 0 : plumed_merror("Problem with function "+bf_str[i]+" given in FUNC"+is+": there should only be one variable");
204 : }
205 46 : if(var_str[0]!=variable_str_) {
206 0 : plumed_merror("Problem with function "+bf_str[i]+" given in FUNC"+is+": you should use "+variable_str_+" as a variable");
207 : }
208 :
209 : try {
210 92 : lepton::ParsedExpression pe_deriv = lepton::Parser::parse(bf_str[i]).differentiate(variable_str_).optimize(lepton::Constants());
211 46 : std::ostringstream tmp_stream2; tmp_stream2 << pe_deriv;
212 46 : bf_derivs_parsed[i] = tmp_stream2.str();
213 46 : bf_derivs_expressions_[i] = pe_deriv.createCompiledExpression();
214 46 : }
215 0 : catch(PLMD::lepton::Exception& exc) {
216 0 : plumed_merror("There was some problem in parsing the derivative of the function "+bf_str[i]+" given in FUNC"+is + " with lepton");
217 0 : }
218 :
219 : try {
220 46 : bf_values_lepton_ref_[i] = &bf_values_expressions_[i].getVariableReference(variable_str_);
221 0 : } catch(PLMD::lepton::Exception& exc) {}
222 :
223 : try {
224 46 : bf_derivs_lepton_ref_[i] = &bf_derivs_expressions_[i].getVariableReference(variable_str_);
225 8 : } catch(PLMD::lepton::Exception& exc) {}
226 :
227 46 : }
228 :
229 : std::string transf_value_parsed;
230 : std::string transf_deriv_parsed;
231 : std::string transf_str;
232 20 : parse("TRANSFORM",transf_str);
233 10 : if(transf_str.size()>0) {
234 6 : do_transf_ = true;
235 12 : addKeywordToList("TRANSFORM",transf_str);
236 : for(unsigned int k=0;; k++) {
237 14 : if(transf_str.find("min")!=std::string::npos) {transf_str.replace(transf_str.find("min"), std::string("min").length(),intervalMinStr());}
238 : else {break;}
239 : }
240 : for(unsigned int k=0;; k++) {
241 14 : if(transf_str.find("max")!=std::string::npos) {transf_str.replace(transf_str.find("max"), std::string("max").length(),intervalMaxStr());}
242 : else {break;}
243 : }
244 :
245 : try {
246 6 : lepton::ParsedExpression pe_value = lepton::Parser::parse(transf_str).optimize(lepton::Constants());;
247 6 : std::ostringstream tmp_stream; tmp_stream << pe_value;
248 6 : transf_value_parsed = tmp_stream.str();
249 6 : transf_value_expression_ = pe_value.createCompiledExpression();
250 6 : }
251 0 : catch(PLMD::lepton::Exception& exc) {
252 0 : plumed_merror("There was some problem in parsing the function "+transf_str+" given in TRANSFORM with lepton");
253 0 : }
254 :
255 : std::vector<std::string> var_str;
256 12 : for(auto &p: transf_value_expression_.getVariables()) {
257 6 : var_str.push_back(p);
258 : }
259 6 : if(var_str.size()!=1) {
260 0 : plumed_merror("Problem with function "+transf_str+" given in TRANSFORM: there should only be one variable");
261 : }
262 6 : if(var_str[0]!=transf_variable_str_) {
263 0 : plumed_merror("Problem with function "+transf_str+" given in TRANSFORM: you should use "+transf_variable_str_+" as a variable");
264 : }
265 :
266 : try {
267 12 : lepton::ParsedExpression pe_deriv = lepton::Parser::parse(transf_str).differentiate(transf_variable_str_).optimize(lepton::Constants());;
268 6 : std::ostringstream tmp_stream2; tmp_stream2 << pe_deriv;
269 6 : transf_deriv_parsed = tmp_stream2.str();
270 6 : transf_deriv_expression_ = pe_deriv.createCompiledExpression();
271 6 : }
272 0 : catch(PLMD::lepton::Exception& exc) {
273 0 : plumed_merror("There was some problem in parsing the derivative of the function "+transf_str+" given in TRANSFORM with lepton");
274 0 : }
275 :
276 : try {
277 6 : transf_value_lepton_ref_ = &transf_value_expression_.getVariableReference(transf_variable_str_);
278 0 : } catch(PLMD::lepton::Exception& exc) {}
279 :
280 : try {
281 6 : transf_deriv_lepton_ref_ = &transf_deriv_expression_.getVariableReference(transf_variable_str_);
282 3 : } catch(PLMD::lepton::Exception& exc) {}
283 :
284 6 : }
285 : //
286 10 : log.printf(" Using the following functions [lepton parsed function and derivative]:\n");
287 66 : for(unsigned int i=0; i<getNumberOfBasisFunctions(); i++) {
288 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());
289 :
290 : }
291 : //
292 10 : if(do_transf_) {
293 6 : log.printf(" Arguments are transformed using the following function [lepton parsed function and derivative]:\n");
294 6 : log.printf(" %s [ %s | %s ] \n",transf_str.c_str(),transf_value_parsed.c_str(),transf_deriv_parsed.c_str());
295 : }
296 : else {
297 4 : log.printf(" Arguments are not transformed\n");
298 : }
299 : //
300 :
301 20 : parseFlag("CHECK_NAN_INF",check_nan_inf_); addKeywordToList("CHECK_NAN_INF",check_nan_inf_);
302 10 : if(check_nan_inf_) {
303 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");
304 : }
305 : else {
306 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");
307 : }
308 :
309 10 : setupBF();
310 10 : checkRead();
311 20 : }
312 :
313 :
314 24204 : void BF_Custom::getAllValues(const double arg, double& argT, bool& inside_range, std::vector<double>& values, std::vector<double>& derivs) const {
315 24204 : inside_range=true;
316 24204 : argT=checkIfArgumentInsideInterval(arg,inside_range);
317 24204 : double transf_derivf=1.0;
318 : //
319 24204 : if(do_transf_) {
320 :
321 14062 : if(transf_value_lepton_ref_) {*transf_value_lepton_ref_ = argT;}
322 14062 : if(transf_deriv_lepton_ref_) {*transf_deriv_lepton_ref_ = argT;}
323 :
324 14062 : argT = transf_value_expression_.evaluate();
325 14062 : transf_derivf = transf_deriv_expression_.evaluate();
326 :
327 14062 : if(check_nan_inf_ && (std::isnan(argT) || std::isinf(argT)) ) {
328 0 : std::string vs; Tools::convert(argT,vs);
329 0 : plumed_merror(getName()+" with label "+getLabel()+": problem with the transform function, it gives " + vs);
330 : }
331 :
332 14062 : if(check_nan_inf_ && (std::isnan(transf_derivf) || std::isinf(transf_derivf)) ) {
333 0 : std::string vs; Tools::convert(transf_derivf,vs);
334 0 : plumed_merror(getName()+" with label "+getLabel()+": problem with the transform function, its derivative gives " + vs);
335 : }
336 : }
337 : //
338 24204 : values[0]=1.0;
339 24204 : derivs[0]=0.0;
340 137488 : for(unsigned int i=1; i < getNumberOfBasisFunctions(); i++) {
341 :
342 113284 : if(bf_values_lepton_ref_[i]) {*bf_values_lepton_ref_[i] = argT;}
343 113284 : if(bf_derivs_lepton_ref_[i]) {*bf_derivs_lepton_ref_[i] = argT;}
344 :
345 113284 : values[i] = bf_values_expressions_[i].evaluate();
346 113284 : derivs[i] = bf_derivs_expressions_[i].evaluate();
347 :
348 113284 : if(do_transf_) {derivs[i]*=transf_derivf;}
349 : // NaN checks
350 113284 : if(check_nan_inf_ && (std::isnan(values[i]) || std::isinf(values[i])) ) {
351 0 : std::string vs; Tools::convert(values[i],vs);
352 0 : std::string is; Tools::convert(i,is);
353 0 : plumed_merror(getName()+" with label "+getLabel()+": problem with the basis function given in FUNC"+is+", it gives "+vs);
354 : }
355 : //
356 113284 : if(check_nan_inf_ && (std::isnan(derivs[i])|| std::isinf(derivs[i])) ) {
357 : std::string vs; Tools::convert(derivs[i],vs);
358 0 : std::string is; Tools::convert(i,is);
359 0 : plumed_merror(getName()+" with label "+getLabel()+": problem with derivative of the basis function given in FUNC"+is+", it gives "+vs);
360 : }
361 : }
362 26218 : if(!inside_range) {for(unsigned int i=0; i<derivs.size(); i++) {derivs[i]=0.0;}}
363 24204 : }
364 :
365 :
366 :
367 :
368 : }
369 : }
|
