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