Nektar++
AnalyticExpressionEvaluator.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File AnalyticExpressionEvaluator.hpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
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
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR 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 OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Parser and evaluator of analytic expressions.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 
38 #include <stdexcept>
39 #include <algorithm>
40 #include <iostream>
41 #include <boost/lexical_cast.hpp>
42 #include <boost/algorithm/string/trim.hpp>
43 #include <boost/random/detail/seed.hpp>
44 
45 #ifdef _MSC_VER
46 #include <boost/preprocessor/cat.hpp>
47 #endif //MSC_VER
48 
49 #ifdef _MSC_VER
50 #define NEKTAR_MATH_NAME(x) BOOST_PP_CAT(_, x)
51 #else
52 #define NEKTAR_MATH_NAME(x) x
53 #endif
54 
55 #if( BOOST_VERSION / 100 % 1000 >= 36 )
56 using namespace boost::spirit::classic;
57 #else
58 using namespace boost::spirit;
59 #endif
60 
61 // trying to avoid incompatibility between standart <algorithm> header and
62 // windows.h header which defines max and min macros.
63 using std::min;
64 
65 namespace Nektar
66 {
67  namespace LibUtilities
68  {
69 
70  // =========================================================================
71  // AnalyticExpression definitions for Spirit Parser
72  // =========================================================================
73 
74  typedef NekDouble (*PFD)();
75  typedef NekDouble (*PFD1)(NekDouble);
79  struct func
80  {
81  func(PFD1 p) : func1(p), size(1) {};
82  func(PFD2 p) : func2(p), size(2) {};
83  func(PFD3 p) : func3(p), size(3) {};
84  func(PFD4 p) : func4(p), size(4) {};
85 
86  union // Pointer to a function
87  {
92  };
93  size_t size;
94  };
95 
96  // signum function
97  NekDouble sign(NekDouble arg)
98  {
99  return (arg > 0.0) - (arg < 0.0);
100  }
101 
102  // Additive white Gaussian noise function.
103  // Arg: sigma of the zero-mean gaussian distribution
104  // Attention: this function is not actually used for
105  // evaluation purposes.
106  NekDouble awgn(NekDouble sigma)
107  {
109  boost::variate_generator<
111  boost::normal_distribution<>
112  > _normal(rng, boost::normal_distribution<>(0, sigma) );
113  return _normal();
114  }
115 
116  /** This struct creates a parser that matches the function
117  definitions from math.h. All of the functions accept one
118  of more NekDoubles as arguments and returns a NekDouble. **/
119  static struct functions : symbols<func>
120  {
122  {
123  // Add all of the functions from math.h
124  add
125  ("abs", std::abs)
126  ("asin", asin)
127  ("acos", acos)
128  ("atan", atan)
129  ("atan2", atan2)
130  ("ang", ang)
131  ("ceil", ceil)
132  ("cos", cos)
133  ("cosh", cosh)
134  ("exp", exp)
135  ("fabs", fabs) // same as abs
136  ("floor", floor)
137  ("log", log)
138  ("log10", log10)
139  ("rad", rad)
140  ("sin", sin)
141  ("sinh", sinh)
142  ("sqrt", sqrt)
143  ("tan", tan)
144  ("tanh", tanh)
145  // and few more custom functions
146  ("sign", sign)
147  ("awgn", awgn)
148  ;
149  }
150  } functions_p;
151 
152 
153  /** This function specifies the grammar of the MathAnalyticExpression parser. **/
154  template <typename ScannerT>
156  {
157  expression = logical_or;
158 
159  logical_or = logical_and >> *( (root_node_d[str_p("||")] >> logical_and) );
160 
161  logical_and = equality >> *( (root_node_d[str_p("&&")] >> equality) );
162 
163  equality = lt_gt >> *( (root_node_d[str_p("==")] >> lt_gt) );
164 
165  lt_gt = add_sub >>
166  *( (root_node_d[str_p("<=")] >> add_sub)
167  | (root_node_d[str_p(">=")] >> add_sub)
168  | (root_node_d[ch_p('<')] >> add_sub)
169  | (root_node_d[ch_p('>')] >> add_sub)
170  );
171 
172  add_sub = mult_div >>
173  *( (root_node_d[ch_p('+')] >> mult_div)
174  | (root_node_d[ch_p('-')] >> mult_div)
175  );
176 
177  mult_div = exponential >>
178  *( (root_node_d[ch_p('*')] >> exponential)
179  | (root_node_d[ch_p('/')] >> exponential)
180  );
181 
182  exponential = factor >>
183  *( (root_node_d[ch_p('^')] >> factor) );
184 
185  factor = number
186  | function
187  | variable
188  | constant
189  | parameter
190  | inner_node_d[ch_p('(') >> expression >> ch_p(')')]
191  | (root_node_d[ch_p('-')] >> factor);
192 
193  parameter = leaf_node_d[ lexeme_d[
194  (alpha_p | '_' | '$') >> *(alnum_p | '_' | '$')
195  ] ] >> op;
196 
197  function = root_node_d[functions_p] >>
198  infix_node_d[inner_node_d[ch_p('(') >> expression >> *(',' >> expression) >> ch_p(')')]];
199 
200  variable = leaf_node_d[ lexeme_d[
201  self.variables_p
202  ] ] >> op;
203 
204  number = leaf_node_d[ lexeme_d[
205  real_p
206  ] ] >> op;
207 
208  constant = leaf_node_d[ lexeme_d[
209  *self.constants_p
210  ] ] >> op;
211 
212  op = eps_p( end_p | "||" | "&&" | "==" | "<=" | ">=" | '<' | '>' | '+' | '-' | '*' | '/' | '^' | ')' | ',' );
213  }
214 
215  // =========================================================================
216  // AnalyticExpressionEvaluator constructor and setting up methods
217  // =========================================================================
218 
219  // \brief Initializes the evaluator. Call DefineFunction(...) next.
220  AnalyticExpressionEvaluator::AnalyticExpressionEvaluator():
221  m_timer(),
222  m_total_eval_time(0)
223  {
224  m_state_size = 1;
225 
226  AddConstant("MEANINGLESS", 0.0);
227  AddConstant("E", 2.71828182845904523536); // Natural logarithm
228  AddConstant("LOG2E", 1.4426950408889634074); // log_2 e
229  AddConstant("LOG10E", 0.43429448190325182765); // log_10 e
230  AddConstant("LN2", 0.69314718055994530942); // log_e 2
231  AddConstant("LN10", 2.30258509299404568402); // log_e 10
232  AddConstant("PI", 3.14159265358979323846); // pi
233  AddConstant("PI_2", 1.57079632679489661923); // pi/2
234  AddConstant("PI_4", 0.78539816339744830962); // pi/4
235  AddConstant("1_PI", 0.31830988618379067154); // 1/pi
236  AddConstant("2_PI", 0.63661977236758134308); // 2/pi
237  AddConstant("2_SQRTPI", 1.12837916709551257390); // 2/sqrt(pi)
238  AddConstant("SQRT2", 1.41421356237309504880); // sqrt(2)
239  AddConstant("SQRT1_2", 0.70710678118654752440); // 1/sqrt(2)
240  AddConstant("GAMMA", 0.57721566490153286060); // Euler
241  AddConstant("DEG", 57.2957795130823208768); // deg/radian
242  AddConstant("PHI", 1.61803398874989484820); // golden ratio
243 
266 
267  m_function[ E_ABS ] = std::abs;
268  m_function[ E_ASIN ] = asin;
269  m_function[ E_ACOS ] = acos;
270  m_function[ E_ATAN ] = atan;
271  m_function[ E_CEIL ] = ceil;
272  m_function[ E_COS ] = cos;
273  m_function[ E_COSH ] = cosh;
274  m_function[ E_EXP ] = exp;
275  m_function[ E_FABS ] = fabs;
276  m_function[ E_FLOOR] = floor;
277  m_function[ E_LOG ] = log;
278  m_function[ E_LOG10] = log10;
279  m_function[ E_SIN ] = sin;
280  m_function[ E_SINH ] = sinh;
281  m_function[ E_SQRT ] = sqrt;
282  m_function[ E_TAN ] = tan;
283  m_function[ E_TANH ] = tanh;
284  m_function[ E_SIGN ] = sign;
285  m_function2[E_ATAN2] = atan2;
286  m_function2[E_ANG] = ang;
287  m_function2[E_RAD] = rad;
288  // there is no entry to m_function that correspond to awgn function.
289  // this is made in purpose. This function need not be pre-evaluated once!
290  }
291 
292 
294  {
295  for (std::vector<ExecutionStack>::iterator it_es = m_executionStack.begin(); it_es != m_executionStack.end(); ++it_es)
296  {
297  for (std::vector<EvaluationStep*>::iterator it = (*it_es).begin(); it != (*it_es).end(); ++it)
298  {
299  delete *it;
300  }
301  (*it_es).clear();
302  }
303  m_executionStack.clear();
304  }
305 
306 
308  {
309  m_generator.seed(seed);
310  }
311 
312 
313  void AnalyticExpressionEvaluator::AddConstants(std::map<std::string, NekDouble> const& constants)
314  {
315  for (std::map<std::string, NekDouble>::const_iterator it = constants.begin(); it != constants.end(); ++it)
316  {
317  AddConstant(it->first, it->second);
318  }
319  }
320 
321  int AnalyticExpressionEvaluator::AddConstant(std::string const& name, NekDouble value)
322  {
323  ConstantMap::const_iterator it = m_constantMapNameToId.find(name);
324  if (it == m_constantMapNameToId.end())
325  {
326  // we are trying to avoid duplicating entries in m_constantParser and m_constants
327  m_constantsParser.add(name.c_str(), value);
328  int index = m_constant.size();
329  m_constantMapNameToId[name] = index;
330  m_constant.push_back(value);
331  return index;
332  }
333  else
334  {
335  if(m_constant[it->second] != value)
336  {
337  std::string errormsg("Attempt to add numerically different constants under the same name: ");
338  errormsg += name;
339  std::cout << errormsg << std::endl;
340  }
341  //ASSERTL1(m_constant[it->second] == value, "Attempt to add numerically different constants under the same name: " + name);
342  }
343  return it->second;
344  }
345 
346  NekDouble AnalyticExpressionEvaluator::GetConstant(std::string const& name)
347  {
348  NekDouble* value = find(m_constantsParser, name.c_str());
349 
350  ASSERTL1(value != NULL, "Constant variable not found: " + name);
351 
352  return *value;
353  }
354 
355  void AnalyticExpressionEvaluator::SetParameters(std::map<std::string, NekDouble> const& params)
356  {
357  for (std::map<std::string, NekDouble>::const_iterator it = params.begin(); it != params.end(); it++)
358  {
359  SetParameter(it->first, it->second);
360  }
361  }
362 
363  void AnalyticExpressionEvaluator::SetParameter(std::string const& name, NekDouble value)
364  {
365  ParameterMap::const_iterator it = m_parameterMapNameToId.find(name);
366  if (it == m_parameterMapNameToId.end())
367  {
368  m_parameterMapNameToId[name] = m_parameter.size();
369  m_parameter.push_back(value);
370  }
371  else
372  {
373  // if parameter is known, change its value
374  m_parameter[ it->second ] = value;
375  }
376  }
377 
378 
379  NekDouble AnalyticExpressionEvaluator::GetParameter(std::string const& name)
380  {
381  ParameterMap::const_iterator it = m_parameterMapNameToId.find(name);
382 
383  ASSERTL1(it != m_parameterMapNameToId.end(), "Parameter not found: " + name);
384 
385  return m_parameter[ it->second ];
386  }
387 
388 
390  {
391  return m_total_eval_time;
392  }
393 
394 
395  // ======================================================
396  // Public evaluate methods
397  // ======================================================
398 
399  int AnalyticExpressionEvaluator::DefineFunction(const std::string& vlist, const std::string& expr)
400  {
401  // Find the previous parsing.
402  ExpressionMap::const_iterator it = m_parsedMapExprToExecStackId.find(expr);
403  if (it != m_parsedMapExprToExecStackId.end())
404  {
405  // if this function is already defined, don't do anything but
406  // return its ID.
407  return it->second;
408  }
409 
410  // ----------------------------------------------
411  // Prepare an iterator that allows to walk along
412  // the string representing an analytic expression in the order
413  // that respects its recursive structure (thanks to boost::spirit).
414  // ----------------------------------------------
415 
416  // Parse the vlist input and separate the variables into ordered entries
417  // in a vector<string> object. These need to be ordered because this is
418  // the order the variables will get assigned to in the Map when Evaluate(...)
419  // is called.
420  std::vector<std::string> variableNames;
421  parse((char*) vlist.c_str(), ( *space_p >>
422  *(
423  +(+graph_p)[push_back_a(variableNames)]
424  >> +space_p
425  )
426  )
427  );
428  // Set up our grammar
429  AnalyticExpression myGrammar(&m_constantsParser, variableNames);
430 
431  // Do the actual parsing with boost::spirit and alert the user if there was an error with an exception.
432  ParsedTreeInfo parseInfo = ast_parse<
433  node_val_data_factory<NekDouble>,
434  std::string::const_iterator,
436  space_parser
437  >
438  (expr.begin(), expr.end(), myGrammar, space_p);
439 
440  ASSERTL1(parseInfo.full != false, "Unable to fully parse function. Stopped just before: "
441  + std::string(parseInfo.stop, parseInfo.stop + 15));
442  if (!parseInfo.full)
443  {
444  throw std::runtime_error("Unable to fully parse function at: "
445  + std::string(parseInfo.stop, parseInfo.stop + 15));
446  }
447 
448 
449  // ----------------------------------------------
450  // Data parsed, start setting up internal data structures.
451  // ----------------------------------------------
452 
453  ExecutionStack stack;
454  VariableMap variableMap;
455 
456  int stackId = m_executionStack.size();
457  m_state_size = 1;
458 
459  // register all variables declared in the expression
460  for (int i = 0; i < variableNames.size(); i++)
461  {
462  variableMap[variableNames[i]] = i;
463  }
464 
465  // then prepare an execution stack.
466  // this method also calculates a length of internal
467  // state storage (m_state_size) for this function.
468  PrecomputedValue v = PrepareExecutionAsYouParse(parseInfo.trees.begin(), stack, variableMap, 0);
469 
470  // constant expression, fully evaluated
471  if (true == v.first)
472  {
473  ASSERTL1(stack.size() == 0, "Constant expression yeilds non-empty execution stack. Bug in PrepareExecutionAsYouParse()");
474 
475  int const_index = AddConstant(std::string("EXPRESSION_") + boost::lexical_cast<std::string>(stackId), v.second);
476  stack.push_back ( makeStep<StoreConst>( 0, const_index ) );
477  }
478 
479  m_parsedMapExprToExecStackId[expr] = stackId;
480 
481  // the execution stack and its corresponding variable index map are
482  // two parallel std::vectors that share their ids. This split helps
483  // to achieve some performance improvement.
484  m_executionStack.push_back(stack);
485  m_stackVariableMap.push_back(variableMap);
486  m_state_sizes.push_back(m_state_size);
487  return stackId;
488  }
489 
490 
491  NekDouble AnalyticExpressionEvaluator::Evaluate(const int expression_id)
492  {
493  m_timer.Start();
494 
495  ASSERTL1(m_executionStack.size() > expression_id, "unknown analytic expression, it must first be defined with DefineFunction(...)");
496 
497  ExecutionStack &stack = m_executionStack[expression_id];
498 
499  m_state.resize(m_state_sizes[expression_id]);
500  for (int i = 0; i < stack.size(); i++)
501  {
502  (*stack[i]).run_once();
503  }
504 
505  m_timer.Stop();
507 
508  return m_state[0];
509  }
510 
512  const int expression_id,
513  const NekDouble x,
514  const NekDouble y,
515  const NekDouble z,
516  const NekDouble t)
517  {
518  m_timer.Start();
519 
520  ASSERTL1(m_executionStack.size() > expression_id, "unknown analytic expression, it must first be defined with DefineFunction(...)");
521 
522  ExecutionStack &stack = m_executionStack[expression_id];
523 
524  // initialise internal vector of variable values
525  m_state.resize(m_state_sizes[expression_id]);
526 
527  if (m_variable.size() < 4)
528  {
529  m_variable.resize(4);
530  }
531 
532  // no flexibility, no change of variable ordering in m_variable
533  // container depending on their names ordering in the input vlist
534  // argument of DefineFunction. Ordering convention (x,y,z,t) is assumed.
535  m_variable[0] = x;
536  m_variable[1] = y;
537  m_variable[2] = z;
538  m_variable[3] = t;
539 
540  // main execution cycle is hidden here
541  for (int i = 0; i < stack.size(); i++)
542  {
543  (*stack[i]).run_once();
544  }
545 
546  m_timer.Stop();
548 
549  return m_state[0];
550  }
551 
552  NekDouble AnalyticExpressionEvaluator::EvaluateAtPoint(const int expression_id, const std::vector<NekDouble> point)
553  {
554  m_timer.Start();
555 
556  ASSERTL1(m_executionStack.size() > expression_id, "unknown analytic expression, it must first be defined with DefineFunction(...)");
557 
558  ExecutionStack& stack = m_executionStack[expression_id];
559  VariableMap& variableMap = m_stackVariableMap[expression_id];
560 
561  ASSERTL1(point.size() == variableMap.size(), "The number of variables used to define this expression should match the point dimensionality.");
562 
563  // initialise internal vector of variable values
564  m_state.resize(m_state_sizes[expression_id]);
565  m_variable.resize(point.size());
566  VariableMap::const_iterator it;
567  for (it = variableMap.begin(); it != variableMap.end(); ++it)
568  {
569  m_variable[it->second] = point[it->second];
570  }
571  // main execution cycle is hidden here
572  for (int i = 0; i < stack.size(); i++)
573  {
574  (*stack[i]).run_once();
575  }
576 
577  m_timer.Stop();
579 
580  return m_state[0];
581  }
582 
583 
584  // wrapper function call
586  const int expression_id,
591  Array<OneD, NekDouble>& result)
592  {
593  m_timer.Start();
594 
595  std::vector<Array<OneD, const NekDouble> >points;
596 
597  points.push_back(x);
598  points.push_back(y);
599  points.push_back(z);
600  points.push_back(t);
601 
602  Evaluate(expression_id,points,result);
603 
604  m_timer.Stop();
606  }
607 
608 
610  const int expression_id,
611  const std::vector<Array<OneD, const NekDouble> > points,
612  Array<OneD, NekDouble>& result)
613  {
614  m_timer.Start();
615 
616  const int num_points = points[0].num_elements();
617  ASSERTL1(m_executionStack.size() > expression_id,
618  "unknown analytic expression, it must first be defined "
619  "with DefineFunction(...)");
620  ASSERTL1(result.num_elements() >= num_points,
621  "destination array must have enough capacity to store "
622  "expression values at each given point");
623 
624  ExecutionStack &stack = m_executionStack[expression_id];
625 
626  /// If number of points tends to 10^6, one may end up
627  /// with up to ~0.5Gb data allocated for m_state only.
628  /// Lets split the work into cache-sized chunks.
629  /// Ahtung, magic constant!
630  const int max_chunk_size = 1024;
631  const int nvals = points.size();
632  const int chunk_size = (std::min)(max_chunk_size, num_points);
633 
634  if (m_state.size() < chunk_size * m_state_sizes[expression_id])
635  {
636  m_state.resize(m_state_sizes[expression_id] * chunk_size, 0.0);
637  }
638  if (m_variable.size() < nvals * chunk_size)
639  {
640  m_variable.resize( nvals * chunk_size, 0.0);
641  }
642  if (result.num_elements() < num_points)
643  {
644  result = Array<OneD, NekDouble>(num_points, 0.0);
645  }
646 
647  int offset = 0;
648  int work_left = num_points;
649  while(work_left > 0)
650  {
651  const int this_chunk_size = (std::min)(work_left, 1024);
652  for (int i = 0; i < this_chunk_size; i++)
653  {
654  for(int j = 0; j < nvals; ++j)
655  {
656  m_variable[i+this_chunk_size*j] = points[j][offset + i];
657  }
658  }
659  for (int i = 0; i < stack.size(); i++)
660  {
661  (*stack[i]).run_many(this_chunk_size);
662  }
663  for (int i = 0; i < this_chunk_size; i++)
664  {
665  result[offset + i] = m_state[i];
666  }
667  work_left -= this_chunk_size;
668  offset += this_chunk_size;
669  }
670  m_timer.Stop();
672  }
673 
674 
676  const ParsedTreeIterator& location,
677  ExecutionStack& stack,
678  VariableMap& variableMap,
679  int stateIndex)
680  {
681  std::string valueStr(location->value.begin(), location->value.end());
682  boost::algorithm::trim(valueStr);
683 
684  const parser_id parserID = location->value.id();
685 #if defined(NEKTAR_DEBUG) || defined(NEKTAR_FULLDEBUG)
686  const int num_children = location->children.size();
687 #endif
688 
689  if (parserID == AnalyticExpression::constantID)
690  {
691  ASSERTL1(num_children == 0, "Illegal children under constant node: " + valueStr);
692 
693  ConstantMap::const_iterator it = m_constantMapNameToId.find(valueStr);
694  ASSERTL1(it != m_constantMapNameToId.end(), "Cannot find the value for the specified constant: " + valueStr);
695 
696  return std::make_pair(true, m_constant[it->second]);
697  }
698  else if (parserID == AnalyticExpression::numberID)
699  {
700  ASSERTL1(num_children == 0, "Illegal children under number node: " + valueStr);
701  return std::make_pair(true, boost::lexical_cast<NekDouble>(valueStr.c_str()) );
702  }
703  else if (parserID == AnalyticExpression::variableID)
704  {
705  ASSERTL1(num_children == 0, "Illegal children under variable node: " + valueStr);
706 
707  VariableMap::const_iterator it = variableMap.find(valueStr);
708  ASSERTL1(it != variableMap.end(), "Unknown variable parsed: " + valueStr);
709 
710  // Variables are not defined at the time of this parse.
711  stack.push_back ( makeStep<StoreVar>( stateIndex, it->second ) );
712  return std::make_pair(false, 0);
713  }
714  else if (parserID == AnalyticExpression::parameterID)
715  {
716  ASSERTL1(num_children == 0, "Illegal children under parameter node: " + valueStr);
717 
718  ParameterMap::const_iterator it = m_parameterMapNameToId.find(valueStr);
719  ASSERTL1(it != m_parameterMapNameToId.end(), "Unknown parameter parsed: " + valueStr);
720 
721  // Parameters may change in between of evalutions.
722  stack.push_back ( makeStep<StorePrm>( stateIndex, it->second ) );
723  return std::make_pair(false, 0);
724  }
725  else if (parserID == AnalyticExpression::functionID)
726  {
727  FunctionNameMap::const_iterator it = m_functionMapNameToInstanceType.find(valueStr);
728  ASSERTL1(it != m_functionMapNameToInstanceType.end(), "Invalid function specified: " + valueStr);
729  ASSERTL1(num_children == 1 || num_children == 2, "Function " + valueStr + " would like to have too few or too many arguments. This is not implemented yet");
730 
731  if (location->children.size() == 1)
732  {
733  PrecomputedValue v = PrepareExecutionAsYouParse(location->children.begin(), stack, variableMap, stateIndex);
734 
735  // additive white gaussian noise function
736  if (it->second == E_AWGN)
737  {
738  int const_index = AddConstant(std::string("SUB_EXPR_") + boost::lexical_cast<std::string>(m_constant.size()), v.second);
739  stack.push_back ( makeStep<StoreConst>( stateIndex, const_index ) );
740  stack.push_back ( makeStep<EvalAWGN>( stateIndex, stateIndex ) );
741  return std::make_pair(false,0);
742  }
743 
744  if (true == v.first)
745  {
746  return std::make_pair( true, m_function[it->second](v.second) );
747  }
748  }
749  else
750  {
751  PrecomputedValue v1 = PrepareExecutionAsYouParse(location->children.begin()+0, stack, variableMap, stateIndex);
752  PrecomputedValue v2 = PrepareExecutionAsYouParse(location->children.begin()+1, stack, variableMap, stateIndex+1);
753  m_state_size++;
754 
755  if (true == v1.first && true == v2.first)
756  {
757  return std::make_pair( true, m_function2[it->second](v1.second, v2.second) );
758  }
759  }
760 
761  // if somewhere down the parse tree there is a variable or parameter, set up an
762  // evaluation sequence.
763  switch (it->second)
764  {
765  case E_ABS:
766  stack.push_back ( makeStep<EvalAbs>( stateIndex, stateIndex) );
767  return std::make_pair(false,0);
768  case E_ASIN:
769  stack.push_back ( makeStep<EvalAsin>( stateIndex, stateIndex) );
770  return std::make_pair(false,0);
771  case E_ACOS:
772  stack.push_back ( makeStep<EvalAcos>( stateIndex, stateIndex ) );
773  return std::make_pair(false,0);
774  case E_ATAN:
775  stack.push_back ( makeStep<EvalAtan>( stateIndex, stateIndex ) );
776  return std::make_pair(false,0);
777  case E_ATAN2:
778  stack.push_back ( makeStep<EvalAtan2>( stateIndex, stateIndex, stateIndex+1 ) );
779  return std::make_pair(false,0);
780  case E_ANG:
781  stack.push_back ( makeStep<EvalAng>( stateIndex, stateIndex, stateIndex+1 ) );
782  return std::make_pair(false,0);
783  case E_CEIL:
784  stack.push_back ( makeStep<EvalCeil>( stateIndex, stateIndex ) );
785  return std::make_pair(false,0);
786  case E_COS:
787  stack.push_back ( makeStep<EvalCos>( stateIndex, stateIndex ) );
788  return std::make_pair(false,0);
789  case E_COSH:
790  stack.push_back ( makeStep<EvalCosh>( stateIndex, stateIndex ) );
791  return std::make_pair(false,0);
792  case E_EXP:
793  stack.push_back ( makeStep<EvalExp>( stateIndex, stateIndex ) );
794  return std::make_pair(false,0);
795  case E_FABS:
796  stack.push_back ( makeStep<EvalFabs>( stateIndex, stateIndex ) );
797  return std::make_pair(false,0);
798  case E_FLOOR:
799  stack.push_back ( makeStep<EvalFloor>( stateIndex, stateIndex ) );
800  return std::make_pair(false,0);
801  case E_LOG:
802  stack.push_back ( makeStep<EvalLog>( stateIndex, stateIndex ) );
803  return std::make_pair(false,0);
804  case E_LOG10:
805  stack.push_back ( makeStep<EvalLog10>( stateIndex, stateIndex ) );
806  return std::make_pair(false,0);
807  case E_RAD:
808  stack.push_back ( makeStep<EvalRad>( stateIndex, stateIndex, stateIndex+1 ) );
809  return std::make_pair(false,0);
810  case E_SIN:
811  stack.push_back ( makeStep<EvalSin>( stateIndex, stateIndex ) );
812  return std::make_pair(false,0);
813  case E_SINH:
814  stack.push_back ( makeStep<EvalSinh>( stateIndex, stateIndex ) );
815  return std::make_pair(false,0);
816  case E_SQRT:
817  stack.push_back ( makeStep<EvalSqrt>( stateIndex, stateIndex ) );
818  return std::make_pair(false,0);
819  case E_TAN:
820  stack.push_back ( makeStep<EvalTan>( stateIndex, stateIndex ) );
821  return std::make_pair(false,0);
822  case E_TANH:
823  stack.push_back ( makeStep<EvalTanh>( stateIndex, stateIndex ) );
824  return std::make_pair(false,0);
825  case E_SIGN:
826  stack.push_back ( makeStep<EvalSign>( stateIndex, stateIndex ) );
827  return std::make_pair(false,0);
828  default:
829  ASSERTL0(false, "Evaluation of " + valueStr + " is not implemented yet");
830  }
831  return std::make_pair(false,0);
832  }
833  else if (parserID == AnalyticExpression::factorID)
834  {
835  ASSERTL1(*valueStr.begin() == '-', "Illegal factor - it can only be '-' and it was: " + valueStr);
836  ASSERTL1(num_children == 1, "Illegal number of children under factor node: " + valueStr);
837 
838  PrecomputedValue v = PrepareExecutionAsYouParse(location->children.begin(), stack, variableMap, stateIndex);
839 
840  // if precomputed value is valid, process it further.
841  if (true == v.first)
842  {
843  return std::make_pair( true, - v.second );
844  }
845  stack.push_back (makeStep<EvalNeg>( stateIndex, stateIndex) );
846  return std::make_pair(false,0);
847  }
848  else if (parserID == AnalyticExpression::operatorID)
849  {
850  ASSERTL1(num_children == 2, "Too few or too many arguments for mathematical operator: " + valueStr);
851  PrecomputedValue left = PrepareExecutionAsYouParse(location->children.begin()+0, stack, variableMap, stateIndex);
852  PrecomputedValue right = PrepareExecutionAsYouParse(location->children.begin()+1, stack, variableMap, stateIndex+1);
853  m_state_size++;
854 
855  // if both precomputed values are valid, process them further.
856  if ((true == left.first) && (true == right.first))
857  {
858  switch(*valueStr.begin())
859  {
860  case '+':
861  return std::make_pair( true, left.second + right.second );
862  case '-':
863  return std::make_pair( true, left.second - right.second );
864  case '*':
865  return std::make_pair( true, left.second * right.second );
866  case '/':
867  return std::make_pair( true, left.second / right.second );
868  case '^':
869  return std::make_pair( true, std::pow(left.second, right.second) );
870  case '=':
871  return std::make_pair( true, left.second == right.second );
872  case '<':
873  if (*(valueStr.end()-1) == '=')
874  {
875  return std::make_pair( true, left.second <= right.second );
876  }
877  else
878  {
879  return std::make_pair( true, left.second < right.second );
880  }
881  return std::make_pair(false,0);
882  case '>':
883  if (*(valueStr.end()-1) == '=')
884  {
885  return std::make_pair( true, left.second >= right.second );
886  }
887  else
888  {
889  return std::make_pair( true, left.second > right.second );
890  }
891  return std::make_pair(false,0);
892  default:
893  ASSERTL0(false, "Invalid operator encountered: " + valueStr);
894  }
895  return std::make_pair(false,0);
896  }
897 
898  // either operator argument is not fully evaluated
899  // add pre-evaluated value to the contaner of constants
900  if (true == left.first)
901  {
902  int const_index = AddConstant(std::string("SUB_EXPR_") + boost::lexical_cast<std::string>(m_constant.size()), left.second);
903  stack.push_back ( makeStep<StoreConst>( stateIndex, const_index ) );
904  }
905  if (true == right.first)
906  {
907  int const_index = AddConstant(std::string("SUB_EXPR_") + boost::lexical_cast<std::string>(m_constant.size()), right.second);
908  stack.push_back ( makeStep<StoreConst>( stateIndex+1, const_index ) );
909  }
910 
911 
912  switch(*valueStr.begin())
913  {
914  case '+':
915  stack.push_back (makeStep<EvalSum>( stateIndex, stateIndex, stateIndex+1 ) );
916  return std::make_pair(false,0);
917  case '-':
918  stack.push_back (makeStep<EvalSub> (stateIndex, stateIndex, stateIndex+1 ) );
919  return std::make_pair(false,0);
920  case '*':
921  stack.push_back (makeStep<EvalMul>( stateIndex, stateIndex, stateIndex+1 ) );
922  return std::make_pair(false,0);
923  case '/':
924  stack.push_back (makeStep<EvalDiv>( stateIndex, stateIndex, stateIndex+1 ) );
925  return std::make_pair(false,0);
926  case '^':
927  stack.push_back (makeStep<EvalPow>( stateIndex, stateIndex, stateIndex+1 ) );
928  return std::make_pair(false,0);
929  case '=':
930  stack.push_back (makeStep<EvalLogicalEqual>( stateIndex, stateIndex, stateIndex+1 ) );
931  return std::make_pair(false,0);
932  case '<':
933  if (*(valueStr.end()-1) == '=')
934  {
935  stack.push_back (makeStep<EvalLogicalLeq>( stateIndex, stateIndex, stateIndex+1 ) );
936  }
937  else
938  {
939  stack.push_back (makeStep<EvalLogicalLess>( stateIndex, stateIndex, stateIndex+1 ) );
940  }
941  return std::make_pair(false,0);
942 
943  case '>':
944  if (*(valueStr.end()-1) == '=')
945  {
946  stack.push_back (makeStep<EvalLogicalGeq>( stateIndex, stateIndex, stateIndex+1 ) );
947  }
948  else
949  {
950  stack.push_back (makeStep<EvalLogicalGreater>( stateIndex, stateIndex, stateIndex+1 ) );
951  }
952  return std::make_pair(false,0);
953 
954  default:
955  ASSERTL0(false, "Invalid operator encountered: " + valueStr);
956  }
957  return std::make_pair(false,0);
958  }
959  ASSERTL0(false, "Illegal expression encountered: " + valueStr);
960  return std::make_pair(false,0);
961  }
962 
963  };
964 };
static NekDouble ang(NekDouble x, NekDouble y)
NekDouble(* PFD3)(NekDouble, NekDouble, NekDouble)
NekDouble GetTime() const
Returns the total time spent in evaluation procedures, seconds.
void AddConstants(std::map< std::string, NekDouble > const &constants)
Constants are evaluated and inserted into the function at the time it is parsed when calling the Defi...
boost_spirit::tree_parse_info< std::string::const_iterator, boost_spirit::node_val_data_factory< NekDouble > > ParsedTreeInfo
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
ExpressionMap m_parsedMapExprToExecStackId
These vector and map store pre-processed evaluation sequences for the analytic expressions. Each ExecutionStack is an ordered container of steps of sequential execution process which evaluates an analytic expression.
std::vector< VariableMap > m_stackVariableMap
Keeping map of variables individually per each analytic expression allows correctly handling expressi...
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:22
void SetParameter(std::string const &name, NekDouble value)
This function behaves in the same way as SetParameters, but it only adds one parameter and it does no...
Nektar::LibUtilities::functions functions_p
int AddConstant(std::string const &name, NekDouble value)
This function behaves in the same way as AddConstants, but it only adds one constant at a time...
NekDouble Evaluate(const int AnalyticExpression_id)
Evaluation method for expressions depending on parameters only.
std::vector< NekDouble > m_state
This vector stores the execution state (memory) used by the sequential execution process.
NekDouble sign(NekDouble arg)
NekDouble EvaluateAtPoint(const int AnalyticExpression_id, std::vector< NekDouble > point)
Evaluation method for expressions depending on unspecified number of variables. This suitable for exp...
void Stop()
Definition: Timer.cpp:62
PrecomputedValue PrepareExecutionAsYouParse(const ParsedTreeIterator &root, ExecutionStack &stack, VariableMap &varMap, int stateIndex)
This method prepares the execution stack (an ordered sequence of operators that perform the evaluatio...
NekDouble GetConstant(std::string const &name)
If a constant with the specified name exists, it returns the NekDouble value that the constant stores...
int m_state_size
This counter is used by PrepareExecutionAsYouParse for finding the minimal state size necessary for e...
ParameterMap m_parameterMapNameToId
The following data structures hold input data to be used on evaluation stage. There are three types o...
NekDouble awgn(NekDouble sigma)
NekDouble(* PFD4)(NekDouble, NekDouble, NekDouble, NekDouble)
double NekDouble
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
static NekDouble rad(NekDouble x, NekDouble y)
int DefineFunction(const std::string &vlist, const std::string &function)
This function allows one to define a function to evaluate. The first argument (vlist) is a list of va...
NekDouble(* PFD1)(NekDouble)
std::vector< int > m_state_sizes
Vector of state sizes per each.
NekDouble GetParameter(std::string const &name)
If a parameter with the specified name exists, it returns the NekDouble value that the parameter stor...
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
Definition: StdRegions.hpp:312
void Start()
Definition: Timer.cpp:51
~AnalyticExpressionEvaluator(void)
Destroys the execution stack.
boost_spirit::tree_match< std::string::const_iterator, boost_spirit::node_val_data_factory< NekDouble > >::tree_iterator ParsedTreeIterator
NekDouble TimePerTest(unsigned int n)
Returns amount of seconds per iteration in a test with n iterations.
Definition: Timer.cpp:108
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
void SetParameters(std::map< std::string, NekDouble > const &params)
Parameters are like constants, but they are inserted into the function at the time Evaluate is called...
NekDouble(* PFD2)(NekDouble, NekDouble)