1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 * -------------------------------------------------------------------------- *
3 * Lepton *
4 * -------------------------------------------------------------------------- *
5 * This is part of the Lepton expression parser originating from *
6 * Simbios, the NIH National Center for Physics-Based Simulation of *
7 * Biological Structures at Stanford, funded under the NIH Roadmap for *
8 * Medical Research, grant U54 GM072970. See https://simtk.org. *
9 * *
10 * Portions copyright (c) 2013-2016 Stanford University and the Authors. *
11 * Authors: Peter Eastman *
12 * Contributors: *
13 * *
14 * Permission is hereby granted, free of charge, to any person obtaining a *
15 * copy of this software and associated documentation files (the "Software"), *
16 * to deal in the Software without restriction, including without limitation *
17 * the rights to use, copy, modify, merge, publish, distribute, sublicense, *
18 * and/or sell copies of the Software, and to permit persons to whom the *
19 * Software is furnished to do so, subject to the following conditions: *
20 * *
21 * The above copyright notice and this permission notice shall be included in *
22 * all copies or substantial portions of the Software. *
23 * *
24 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
25 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
26 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
27 * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
28 * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
29 * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
30 * USE OR OTHER DEALINGS IN THE SOFTWARE. *
31 * -------------------------------------------------------------------------- *
32 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
33 /* -------------------------------------------------------------------------- *
34 * lepton *
35 * -------------------------------------------------------------------------- *
36 * This is part of the lepton expression parser originating from *
37 * Simbios, the NIH National Center for Physics-Based Simulation of *
38 * Biological Structures at Stanford, funded under the NIH Roadmap for *
39 * Medical Research, grant U54 GM072970. See https://simtk.org. *
40 * *
41 * Portions copyright (c) 2009-2015 Stanford University and the Authors. *
42 * Authors: Peter Eastman *
43 * Contributors: *
44 * *
45 * Permission is hereby granted, free of charge, to any person obtaining a *
46 * copy of this software and associated documentation files (the "Software"), *
47 * to deal in the Software without restriction, including without limitation *
48 * the rights to use, copy, modify, merge, publish, distribute, sublicense, *
49 * and/or sell copies of the Software, and to permit persons to whom the *
50 * Software is furnished to do so, subject to the following conditions: *
51 * *
52 * The above copyright notice and this permission notice shall be included in *
53 * all copies or substantial portions of the Software. *
54 * *
55 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
56 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
57 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
58 * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
59 * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
60 * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
61 * USE OR OTHER DEALINGS IN THE SOFTWARE. *
62 * -------------------------------------------------------------------------- */
63
64 #include "Parser.h"
65 #include "CustomFunction.h"
66 #include "Exception.h"
67 #include "ExpressionTreeNode.h"
68 #include "Operation.h"
69 #include "ParsedExpression.h"
70 #include <cctype>
71 #include <iostream>
72
73 namespace PLMD {
74 using namespace lepton;
75 using namespace std;
76
77 namespace lepton {
78
79 static const string Digits = "0123456789";
80 static const string Operators = "+-*/^";
81 static const bool LeftAssociative[] = {true, true, true, true, false};
82 static const int Precedence[] = {0, 0, 1, 1, 3};
83 static const Operation::Id OperationId[] = {Operation::ADD, Operation::SUBTRACT, Operation::MULTIPLY, Operation::DIVIDE, Operation::POWER};
84
Constants()85 const std::map<std::string, double> & Constants() {
86 static const std::map<std::string, double> constants = {
87 {"e", std::exp(1.0)},
88 {"log2e", 1.0/std::log(2.0)},
89 {"log10e", 1.0/std::log(10.0)},
90 {"ln2", std::log(2.0)},
91 {"ln10", std::log(10.0)},
92 {"pi", 3.14159265358979323844},
93 {"pi_2", 3.14159265358979323844*0.5},
94 {"pi_4", 3.14159265358979323844*0.25},
95 // {"1_pi", 1.0/pi},
96 // {"2_pi", 2.0/pi},
97 // {"2_sqrtpi", 2.0/std::sqrt(pi)},
98 {"sqrt2", std::sqrt(2.0)},
99 {"sqrt1_2", std::sqrt(0.5)}
100 };
101 return constants;
102 }
103
104 }
105
106
107 class lepton::ParseToken {
108 public:
109 enum Type {Number, Operator, Variable, Function, LeftParen, RightParen, Comma, Whitespace};
110
ParseToken(string text,Type type)111 ParseToken(string text, Type type) : text(text), type(type) {
112 }
getText() const113 const string& getText() const {
114 return text;
115 }
getType() const116 Type getType() const {
117 return type;
118 }
119 private:
120 string text;
121 Type type;
122 };
123
trim(const string & expression)124 string Parser::trim(const string& expression) {
125 // Remove leading and trailing spaces.
126
127 int start, end;
128 for (start = 0; start < (int) expression.size() && isspace(expression[start]); start++)
129 ;
130 for (end = (int) expression.size()-1; end > start && isspace(expression[end]); end--)
131 ;
132 if (start == end && isspace(expression[end]))
133 return "";
134 return expression.substr(start, end-start+1);
135 }
136
getNextToken(const string & expression,int start)137 ParseToken Parser::getNextToken(const string& expression, int start) {
138 char c = expression[start];
139 if (c == '(')
140 return ParseToken("(", ParseToken::LeftParen);
141 if (c == ')')
142 return ParseToken(")", ParseToken::RightParen);
143 if (c == ',')
144 return ParseToken(",", ParseToken::Comma);
145 if (Operators.find(c) != string::npos)
146 return ParseToken(string(1, c), ParseToken::Operator);
147 if (isspace(c)) {
148 // White space
149
150 for (int pos = start+1; pos < (int) expression.size(); pos++) {
151 if (!isspace(expression[pos]))
152 return ParseToken(expression.substr(start, pos-start), ParseToken::Whitespace);
153 }
154 return ParseToken(expression.substr(start, string::npos), ParseToken::Whitespace);
155 }
156 if (c == '.' || Digits.find(c) != string::npos) {
157 // A number
158
159 bool foundDecimal = (c == '.');
160 bool foundExp = false;
161 int pos;
162 for (pos = start+1; pos < (int) expression.size(); pos++) {
163 c = expression[pos];
164 if (Digits.find(c) != string::npos)
165 continue;
166 if (c == '.' && !foundDecimal) {
167 foundDecimal = true;
168 continue;
169 }
170 if ((c == 'e' || c == 'E') && !foundExp) {
171 foundExp = true;
172 if (pos < (int) expression.size()-1 && (expression[pos+1] == '-' || expression[pos+1] == '+'))
173 pos++;
174 continue;
175 }
176 break;
177 }
178 return ParseToken(expression.substr(start, pos-start), ParseToken::Number);
179 }
180
181 // A variable, function, or left parenthesis
182
183 for (int pos = start; pos < (int) expression.size(); pos++) {
184 c = expression[pos];
185 if (c == '(')
186 return ParseToken(expression.substr(start, pos-start+1), ParseToken::Function);
187 if (Operators.find(c) != string::npos || c == ',' || c == ')' || isspace(c))
188 return ParseToken(expression.substr(start, pos-start), ParseToken::Variable);
189 }
190 return ParseToken(expression.substr(start, string::npos), ParseToken::Variable);
191 }
192
tokenize(const string & expression)193 vector<ParseToken> Parser::tokenize(const string& expression) {
194 vector<ParseToken> tokens;
195 int pos = 0;
196 while (pos < (int) expression.size()) {
197 ParseToken token = getNextToken(expression, pos);
198 if (token.getType() != ParseToken::Whitespace)
199 tokens.push_back(token);
200 pos += (int) token.getText().size();
201 }
202 return tokens;
203 }
204
parse(const string & expression)205 ParsedExpression Parser::parse(const string& expression) {
206 return parse(expression, map<string, CustomFunction*>());
207 }
208
parse(const string & expression,const map<string,CustomFunction * > & customFunctions)209 ParsedExpression Parser::parse(const string& expression, const map<string, CustomFunction*>& customFunctions) {
210 try {
211 // First split the expression into subexpressions.
212
213 string primaryExpression = expression;
214 vector<string> subexpressions;
215 while (true) {
216 string::size_type pos = primaryExpression.find_last_of(';');
217 if (pos == string::npos)
218 break;
219 string sub = trim(primaryExpression.substr(pos+1));
220 if (sub.size() > 0)
221 subexpressions.push_back(sub);
222 primaryExpression = primaryExpression.substr(0, pos);
223 }
224
225 // Parse the subexpressions.
226
227 map<string, ExpressionTreeNode> subexpDefs;
228 for (int i = 0; i < (int) subexpressions.size(); i++) {
229 string::size_type equalsPos = subexpressions[i].find('=');
230 if (equalsPos == string::npos)
231 throw Exception("subexpression does not specify a name");
232 string name = trim(subexpressions[i].substr(0, equalsPos));
233 if (name.size() == 0)
234 throw Exception("subexpression does not specify a name");
235 vector<ParseToken> tokens = tokenize(subexpressions[i].substr(equalsPos+1));
236 int pos = 0;
237 subexpDefs[name] = parsePrecedence(tokens, pos, customFunctions, subexpDefs, 0);
238 if (pos != tokens.size())
239 throw Exception("unexpected text at end of subexpression: "+tokens[pos].getText());
240 }
241
242 // Now parse the primary expression.
243
244 vector<ParseToken> tokens = tokenize(primaryExpression);
245 int pos = 0;
246 ExpressionTreeNode result = parsePrecedence(tokens, pos, customFunctions, subexpDefs, 0);
247 if (pos != tokens.size())
248 throw Exception("unexpected text at end of expression: "+tokens[pos].getText());
249 return ParsedExpression(result);
250 }
251 catch (Exception& ex) {
252 throw Exception("Parse error in expression \""+expression+"\": "+ex.what());
253 }
254 }
255
parsePrecedence(const vector<ParseToken> & tokens,int & pos,const map<string,CustomFunction * > & customFunctions,const map<string,ExpressionTreeNode> & subexpressionDefs,int precedence)256 ExpressionTreeNode Parser::parsePrecedence(const vector<ParseToken>& tokens, int& pos, const map<string, CustomFunction*>& customFunctions,
257 const map<string, ExpressionTreeNode>& subexpressionDefs, int precedence) {
258 if (pos == tokens.size())
259 throw Exception("unexpected end of expression");
260
261 // Parse the next value (number, variable, function, parenthesized expression)
262
263 ParseToken token = tokens[pos];
264 ExpressionTreeNode result;
265 if (token.getType() == ParseToken::Number) {
266 double value;
267 stringstream(token.getText()) >> value;
268 result = ExpressionTreeNode(new Operation::Constant(value));
269 pos++;
270 }
271 else if (token.getType() == ParseToken::Variable) {
272 map<string, ExpressionTreeNode>::const_iterator subexp = subexpressionDefs.find(token.getText());
273 if (subexp == subexpressionDefs.end()) {
274 Operation* op = new Operation::Variable(token.getText());
275 result = ExpressionTreeNode(op);
276 }
277 else
278 result = subexp->second;
279 pos++;
280 }
281 else if (token.getType() == ParseToken::LeftParen) {
282 pos++;
283 result = parsePrecedence(tokens, pos, customFunctions, subexpressionDefs, 0);
284 if (pos == tokens.size() || tokens[pos].getType() != ParseToken::RightParen)
285 throw Exception("unbalanced parentheses");
286 pos++;
287 }
288 else if (token.getType() == ParseToken::Function) {
289 pos++;
290 vector<ExpressionTreeNode> args;
291 bool moreArgs;
292 do {
293 args.push_back(parsePrecedence(tokens, pos, customFunctions, subexpressionDefs, 0));
294 moreArgs = (pos < (int) tokens.size() && tokens[pos].getType() == ParseToken::Comma);
295 if (moreArgs)
296 pos++;
297 } while (moreArgs);
298 if (pos == tokens.size() || tokens[pos].getType() != ParseToken::RightParen)
299 throw Exception("unbalanced parentheses");
300 pos++;
301 Operation* op = getFunctionOperation(token.getText(), customFunctions);
302 try {
303 result = ExpressionTreeNode(op, args);
304 }
305 catch (...) {
306 delete op;
307 throw;
308 }
309 }
310 else if (token.getType() == ParseToken::Operator && token.getText() == "-") {
311 pos++;
312 ExpressionTreeNode toNegate = parsePrecedence(tokens, pos, customFunctions, subexpressionDefs, 2);
313 result = ExpressionTreeNode(new Operation::Negate(), toNegate);
314 }
315 else
316 throw Exception("unexpected token: "+token.getText());
317
318 // Now deal with the next binary operator.
319
320 while (pos < (int) tokens.size() && tokens[pos].getType() == ParseToken::Operator) {
321 token = tokens[pos];
322 int opIndex = (int) Operators.find(token.getText());
323 int opPrecedence = Precedence[opIndex];
324 if (opPrecedence < precedence)
325 return result;
326 pos++;
327 ExpressionTreeNode arg = parsePrecedence(tokens, pos, customFunctions, subexpressionDefs, LeftAssociative[opIndex] ? opPrecedence+1 : opPrecedence);
328 Operation* op = getOperatorOperation(token.getText());
329 try {
330 result = ExpressionTreeNode(op, result, arg);
331 }
332 catch (...) {
333 delete op;
334 throw;
335 }
336 }
337 return result;
338 }
339
getOperatorOperation(const std::string & name)340 Operation* Parser::getOperatorOperation(const std::string& name) {
341 switch (OperationId[Operators.find(name)]) {
342 case Operation::ADD:
343 return new Operation::Add();
344 case Operation::SUBTRACT:
345 return new Operation::Subtract();
346 case Operation::MULTIPLY:
347 return new Operation::Multiply();
348 case Operation::DIVIDE:
349 return new Operation::Divide();
350 case Operation::POWER:
351 return new Operation::Power();
352 default:
353 throw Exception("unknown operator");
354 }
355 }
356
getFunctionOperation(const std::string & name,const map<string,CustomFunction * > & customFunctions)357 Operation* Parser::getFunctionOperation(const std::string& name, const map<string, CustomFunction*>& customFunctions) {
358
359 const static map<string, Operation::Id> opMap ={
360 { "sqrt" , Operation::SQRT },
361 { "exp" , Operation::EXP },
362 { "log" , Operation::LOG },
363 { "sin" , Operation::SIN },
364 { "cos" , Operation::COS },
365 { "sec" , Operation::SEC },
366 { "csc" , Operation::CSC },
367 { "tan" , Operation::TAN },
368 { "cot" , Operation::COT },
369 { "asin" , Operation::ASIN },
370 { "acos" , Operation::ACOS },
371 { "atan" , Operation::ATAN },
372 { "sinh" , Operation::SINH },
373 { "cosh" , Operation::COSH },
374 { "tanh" , Operation::TANH },
375 { "erf" , Operation::ERF },
376 { "erfc" , Operation::ERFC },
377 { "step" , Operation::STEP },
378 { "delta" , Operation::DELTA },
379 { "nandelta" , Operation::NANDELTA },
380 { "square" , Operation::SQUARE },
381 { "cube", Operation::CUBE },
382 { "recip" , Operation::RECIPROCAL },
383 { "min" , Operation::MIN },
384 { "max" , Operation::MAX },
385 { "abs" , Operation::ABS },
386 { "floor" , Operation::FLOOR },
387 { "ceil" , Operation::CEIL },
388 { "select" , Operation::SELECT },
389 { "acot" , Operation::ACOT },
390 { "asec" , Operation::ASEC },
391 { "acsc" , Operation::ACSC },
392 { "coth" , Operation::COTH },
393 { "sech" , Operation::SECH },
394 { "csch" , Operation::CSCH },
395 { "asinh" , Operation::ASINH },
396 { "acosh" , Operation::ACOSH },
397 { "atanh" , Operation::ATANH },
398 { "acoth" , Operation::ACOTH },
399 { "asech" , Operation::ASECH },
400 { "acsch" , Operation::ACSCH },
401 { "atan2" , Operation::ATAN2 },
402 };
403 string trimmed = name.substr(0, name.size()-1);
404
405 // First check custom functions.
406
407 map<string, CustomFunction*>::const_iterator custom = customFunctions.find(trimmed);
408 if (custom != customFunctions.end())
409 return new Operation::Custom(trimmed, custom->second->clone());
410
411 // Now try standard functions.
412
413 map<string, Operation::Id>::const_iterator iter = opMap.find(trimmed);
414 if (iter == opMap.end())
415 throw Exception("unknown function: "+trimmed);
416 switch (iter->second) {
417 case Operation::SQRT:
418 return new Operation::Sqrt();
419 case Operation::EXP:
420 return new Operation::Exp();
421 case Operation::LOG:
422 return new Operation::Log();
423 case Operation::SIN:
424 return new Operation::Sin();
425 case Operation::COS:
426 return new Operation::Cos();
427 case Operation::SEC:
428 return new Operation::Sec();
429 case Operation::CSC:
430 return new Operation::Csc();
431 case Operation::TAN:
432 return new Operation::Tan();
433 case Operation::COT:
434 return new Operation::Cot();
435 case Operation::ASIN:
436 return new Operation::Asin();
437 case Operation::ACOS:
438 return new Operation::Acos();
439 case Operation::ATAN:
440 return new Operation::Atan();
441 case Operation::SINH:
442 return new Operation::Sinh();
443 case Operation::COSH:
444 return new Operation::Cosh();
445 case Operation::TANH:
446 return new Operation::Tanh();
447 case Operation::ERF:
448 return new Operation::Erf();
449 case Operation::ERFC:
450 return new Operation::Erfc();
451 case Operation::STEP:
452 return new Operation::Step();
453 case Operation::DELTA:
454 return new Operation::Delta();
455 case Operation::NANDELTA:
456 return new Operation::Nandelta();
457 case Operation::SQUARE:
458 return new Operation::Square();
459 case Operation::CUBE:
460 return new Operation::Cube();
461 case Operation::RECIPROCAL:
462 return new Operation::Reciprocal();
463 case Operation::MIN:
464 return new Operation::Min();
465 case Operation::MAX:
466 return new Operation::Max();
467 case Operation::ABS:
468 return new Operation::Abs();
469 case Operation::FLOOR:
470 return new Operation::Floor();
471 case Operation::CEIL:
472 return new Operation::Ceil();
473 case Operation::SELECT:
474 return new Operation::Select();
475 case Operation::ACOT:
476 return new Operation::Acot();
477 case Operation::ASEC:
478 return new Operation::Asec();
479 case Operation::ACSC:
480 return new Operation::Acsc();
481 case Operation::COTH:
482 return new Operation::Coth();
483 case Operation::SECH:
484 return new Operation::Sech();
485 case Operation::CSCH:
486 return new Operation::Csch();
487 case Operation::ASINH:
488 return new Operation::Asinh();
489 case Operation::ACOSH:
490 return new Operation::Acosh();
491 case Operation::ATANH:
492 return new Operation::Atanh();
493 case Operation::ACOTH:
494 return new Operation::Acoth();
495 case Operation::ASECH:
496 return new Operation::Asech();
497 case Operation::ACSCH:
498 return new Operation::Acsch();
499 case Operation::ATAN2:
500 return new Operation::Atan2();
501 default:
502 throw Exception("unknown function");
503 }
504 }
505 }
506