Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 
267 
268  m_function[ E_ABS ] = std::abs;
269  m_function[ E_ASIN ] = asin;
270  m_function[ E_ACOS ] = acos;
271  m_function[ E_ATAN ] = atan;
272  m_function[ E_CEIL ] = ceil;
273  m_function[ E_COS ] = cos;
274  m_function[ E_COSH ] = cosh;
275  m_function[ E_EXP ] = exp;
276  m_function[ E_FABS ] = fabs;
277  m_function[ E_FLOOR] = floor;
278  m_function[ E_LOG ] = log;
279  m_function[ E_LOG10] = log10;
280  m_function[ E_SIN ] = sin;
281  m_function[ E_SINH ] = sinh;
282  m_function[ E_SQRT ] = sqrt;
283  m_function[ E_TAN ] = tan;
284  m_function[ E_TANH ] = tanh;
285  m_function[ E_SIGN ] = sign;
286  m_function2[E_ATAN2] = atan2;
287  m_function2[E_ANG] = ang;
288  m_function2[E_RAD] = rad;
289  m_function2[E_BESSEL] = boost::math::cyl_bessel_j;
290 
291  // there is no entry to m_function that correspond to awgn function.
292  // this is made in purpose. This function need not be pre-evaluated once!
293  }
294 
295 
297  {
298  for (std::vector<ExecutionStack>::iterator it_es = m_executionStack.begin(); it_es != m_executionStack.end(); ++it_es)
299  {
300  for (std::vector<EvaluationStep*>::iterator it = (*it_es).begin(); it != (*it_es).end(); ++it)
301  {
302  delete *it;
303  }
304  (*it_es).clear();
305  }
306  m_executionStack.clear();
307  }
308 
309 
311  {
312  m_generator.seed(seed);
313  }
314 
315 
316  void AnalyticExpressionEvaluator::AddConstants(std::map<std::string, NekDouble> const& constants)
317  {
318  for (std::map<std::string, NekDouble>::const_iterator it = constants.begin(); it != constants.end(); ++it)
319  {
320  AddConstant(it->first, it->second);
321  }
322  }
323 
324  int AnalyticExpressionEvaluator::AddConstant(std::string const& name, NekDouble value)
325  {
326  ConstantMap::const_iterator it = m_constantMapNameToId.find(name);
327  if (it == m_constantMapNameToId.end())
328  {
329  // we are trying to avoid duplicating entries in m_constantParser and m_constants
330  m_constantsParser.add(name.c_str(), value);
331  int index = m_constant.size();
332  m_constantMapNameToId[name] = index;
333  m_constant.push_back(value);
334  return index;
335  }
336  else
337  {
338  if(m_constant[it->second] != value)
339  {
340  std::string errormsg("Attempt to add numerically different constants under the same name: ");
341  errormsg += name;
342  std::cout << errormsg << std::endl;
343  }
344  //ASSERTL1(m_constant[it->second] == value, "Attempt to add numerically different constants under the same name: " + name);
345  }
346  return it->second;
347  }
348 
349  NekDouble AnalyticExpressionEvaluator::GetConstant(std::string const& name)
350  {
351  NekDouble* value = find(m_constantsParser, name.c_str());
352 
353  ASSERTL1(value != NULL, "Constant variable not found: " + name);
354 
355  return *value;
356  }
357 
358  void AnalyticExpressionEvaluator::SetParameters(std::map<std::string, NekDouble> const& params)
359  {
360  for (std::map<std::string, NekDouble>::const_iterator it = params.begin(); it != params.end(); it++)
361  {
362  SetParameter(it->first, it->second);
363  }
364  }
365 
366  void AnalyticExpressionEvaluator::SetParameter(std::string const& name, NekDouble value)
367  {
368  ParameterMap::const_iterator it = m_parameterMapNameToId.find(name);
369  if (it == m_parameterMapNameToId.end())
370  {
371  m_parameterMapNameToId[name] = m_parameter.size();
372  m_parameter.push_back(value);
373  }
374  else
375  {
376  // if parameter is known, change its value
377  m_parameter[ it->second ] = value;
378  }
379  }
380 
381 
382  NekDouble AnalyticExpressionEvaluator::GetParameter(std::string const& name)
383  {
384  ParameterMap::const_iterator it = m_parameterMapNameToId.find(name);
385 
386  ASSERTL1(it != m_parameterMapNameToId.end(), "Parameter not found: " + name);
387 
388  return m_parameter[ it->second ];
389  }
390 
391 
393  {
394  return m_total_eval_time;
395  }
396 
397 
398  // ======================================================
399  // Public evaluate methods
400  // ======================================================
401 
402  int AnalyticExpressionEvaluator::DefineFunction(const std::string& vlist, const std::string& expr)
403  {
404  // Find the previous parsing.
405  ExpressionMap::const_iterator it = m_parsedMapExprToExecStackId.find(expr);
406  if (it != m_parsedMapExprToExecStackId.end())
407  {
408  // if this function is already defined, don't do anything but
409  // return its ID.
410  return it->second;
411  }
412 
413  // ----------------------------------------------
414  // Prepare an iterator that allows to walk along
415  // the string representing an analytic expression in the order
416  // that respects its recursive structure (thanks to boost::spirit).
417  // ----------------------------------------------
418 
419  // Parse the vlist input and separate the variables into ordered entries
420  // in a vector<string> object. These need to be ordered because this is
421  // the order the variables will get assigned to in the Map when Evaluate(...)
422  // is called.
423  std::vector<std::string> variableNames;
424  parse((char*) vlist.c_str(), ( *space_p >>
425  *(
426  +(+graph_p)[push_back_a(variableNames)]
427  >> +space_p
428  )
429  )
430  );
431  // Set up our grammar
432  AnalyticExpression myGrammar(&m_constantsParser, variableNames);
433 
434  // Do the actual parsing with boost::spirit and alert the user if there was an error with an exception.
435  ParsedTreeInfo parseInfo = ast_parse<
436  node_val_data_factory<NekDouble>,
437  std::string::const_iterator,
439  space_parser
440  >
441  (expr.begin(), expr.end(), myGrammar, space_p);
442 
443  ASSERTL1(parseInfo.full != false, "Unable to fully parse function. Stopped just before: "
444  + std::string(parseInfo.stop, parseInfo.stop + 15));
445  if (!parseInfo.full)
446  {
447  throw std::runtime_error("Unable to fully parse function at: "
448  + std::string(parseInfo.stop, parseInfo.stop + 15));
449  }
450 
451 
452  // ----------------------------------------------
453  // Data parsed, start setting up internal data structures.
454  // ----------------------------------------------
455 
456  ExecutionStack stack;
457  VariableMap variableMap;
458 
459  int stackId = m_executionStack.size();
460  m_state_size = 1;
461 
462  // register all variables declared in the expression
463  for (int i = 0; i < variableNames.size(); i++)
464  {
465  variableMap[variableNames[i]] = i;
466  }
467 
468  // then prepare an execution stack.
469  // this method also calculates a length of internal
470  // state storage (m_state_size) for this function.
471  PrecomputedValue v = PrepareExecutionAsYouParse(parseInfo.trees.begin(), stack, variableMap, 0);
472 
473  // constant expression, fully evaluated
474  if (true == v.first)
475  {
476  ASSERTL1(stack.size() == 0, "Constant expression yeilds non-empty execution stack. Bug in PrepareExecutionAsYouParse()");
477 
478  int const_index = AddConstant(std::string("EXPRESSION_") + boost::lexical_cast<std::string>(stackId), v.second);
479  stack.push_back ( makeStep<StoreConst>( 0, const_index ) );
480  }
481 
482  m_parsedMapExprToExecStackId[expr] = stackId;
483 
484  // the execution stack and its corresponding variable index map are
485  // two parallel std::vectors that share their ids. This split helps
486  // to achieve some performance improvement.
487  m_executionStack.push_back(stack);
488  m_stackVariableMap.push_back(variableMap);
489  m_state_sizes.push_back(m_state_size);
490  return stackId;
491  }
492 
493 
494  NekDouble AnalyticExpressionEvaluator::Evaluate(const int expression_id)
495  {
496  m_timer.Start();
497 
498  ASSERTL1(m_executionStack.size() > expression_id, "unknown analytic expression, it must first be defined with DefineFunction(...)");
499 
500  ExecutionStack &stack = m_executionStack[expression_id];
501 
502  m_state.resize(m_state_sizes[expression_id]);
503  for (int i = 0; i < stack.size(); i++)
504  {
505  (*stack[i]).run_once();
506  }
507 
508  m_timer.Stop();
510 
511  return m_state[0];
512  }
513 
515  const int expression_id,
516  const NekDouble x,
517  const NekDouble y,
518  const NekDouble z,
519  const NekDouble t)
520  {
521  m_timer.Start();
522 
523  ASSERTL1(m_executionStack.size() > expression_id, "unknown analytic expression, it must first be defined with DefineFunction(...)");
524 
525  ExecutionStack &stack = m_executionStack[expression_id];
526 
527  // initialise internal vector of variable values
528  m_state.resize(m_state_sizes[expression_id]);
529 
530  if (m_variable.size() < 4)
531  {
532  m_variable.resize(4);
533  }
534 
535  // no flexibility, no change of variable ordering in m_variable
536  // container depending on their names ordering in the input vlist
537  // argument of DefineFunction. Ordering convention (x,y,z,t) is assumed.
538  m_variable[0] = x;
539  m_variable[1] = y;
540  m_variable[2] = z;
541  m_variable[3] = t;
542 
543  // main execution cycle is hidden here
544  for (int i = 0; i < stack.size(); i++)
545  {
546  (*stack[i]).run_once();
547  }
548 
549  m_timer.Stop();
551 
552  return m_state[0];
553  }
554 
555  NekDouble AnalyticExpressionEvaluator::EvaluateAtPoint(const int expression_id, const std::vector<NekDouble> point)
556  {
557  m_timer.Start();
558 
559  ASSERTL1(m_executionStack.size() > expression_id, "unknown analytic expression, it must first be defined with DefineFunction(...)");
560 
561  ExecutionStack& stack = m_executionStack[expression_id];
562  VariableMap& variableMap = m_stackVariableMap[expression_id];
563 
564  ASSERTL1(point.size() == variableMap.size(), "The number of variables used to define this expression should match the point dimensionality.");
565 
566  // initialise internal vector of variable values
567  m_state.resize(m_state_sizes[expression_id]);
568  m_variable.resize(point.size());
569  VariableMap::const_iterator it;
570  for (it = variableMap.begin(); it != variableMap.end(); ++it)
571  {
572  m_variable[it->second] = point[it->second];
573  }
574  // main execution cycle is hidden here
575  for (int i = 0; i < stack.size(); i++)
576  {
577  (*stack[i]).run_once();
578  }
579 
580  m_timer.Stop();
582 
583  return m_state[0];
584  }
585 
586 
587  // wrapper function call
589  const int expression_id,
594  Array<OneD, NekDouble>& result)
595  {
596  m_timer.Start();
597 
598  std::vector<Array<OneD, const NekDouble> >points;
599 
600  points.push_back(x);
601  points.push_back(y);
602  points.push_back(z);
603  points.push_back(t);
604 
605  Evaluate(expression_id,points,result);
606 
607  m_timer.Stop();
609  }
610 
611 
613  const int expression_id,
614  const std::vector<Array<OneD, const NekDouble> > points,
615  Array<OneD, NekDouble>& result)
616  {
617  m_timer.Start();
618 
619  const int num_points = points[0].num_elements();
620  ASSERTL1(m_executionStack.size() > expression_id,
621  "unknown analytic expression, it must first be defined "
622  "with DefineFunction(...)");
623  ASSERTL1(result.num_elements() >= num_points,
624  "destination array must have enough capacity to store "
625  "expression values at each given point");
626 
627  ExecutionStack &stack = m_executionStack[expression_id];
628 
629  /// If number of points tends to 10^6, one may end up
630  /// with up to ~0.5Gb data allocated for m_state only.
631  /// Lets split the work into cache-sized chunks.
632  /// Ahtung, magic constant!
633  const int max_chunk_size = 1024;
634  const int nvals = points.size();
635  const int chunk_size = (std::min)(max_chunk_size, num_points);
636 
637  if (m_state.size() < chunk_size * m_state_sizes[expression_id])
638  {
639  m_state.resize(m_state_sizes[expression_id] * chunk_size, 0.0);
640  }
641  if (m_variable.size() < nvals * chunk_size)
642  {
643  m_variable.resize( nvals * chunk_size, 0.0);
644  }
645  if (result.num_elements() < num_points)
646  {
647  result = Array<OneD, NekDouble>(num_points, 0.0);
648  }
649 
650  int offset = 0;
651  int work_left = num_points;
652  while(work_left > 0)
653  {
654  const int this_chunk_size = (std::min)(work_left, 1024);
655  for (int i = 0; i < this_chunk_size; i++)
656  {
657  for(int j = 0; j < nvals; ++j)
658  {
659  m_variable[i+this_chunk_size*j] = points[j][offset + i];
660  }
661  }
662  for (int i = 0; i < stack.size(); i++)
663  {
664  (*stack[i]).run_many(this_chunk_size);
665  }
666  for (int i = 0; i < this_chunk_size; i++)
667  {
668  result[offset + i] = m_state[i];
669  }
670  work_left -= this_chunk_size;
671  offset += this_chunk_size;
672  }
673  m_timer.Stop();
675  }
676 
677 
679  const ParsedTreeIterator& location,
680  ExecutionStack& stack,
681  VariableMap& variableMap,
682  int stateIndex)
683  {
684  std::string valueStr(location->value.begin(), location->value.end());
685  boost::algorithm::trim(valueStr);
686 
687  const parser_id parserID = location->value.id();
688 #if defined(NEKTAR_DEBUG) || defined(NEKTAR_FULLDEBUG)
689  const int num_children = location->children.size();
690 #endif
691 
692  if (parserID == AnalyticExpression::constantID)
693  {
694  ASSERTL1(num_children == 0, "Illegal children under constant node: " + valueStr);
695 
696  ConstantMap::const_iterator it = m_constantMapNameToId.find(valueStr);
697  ASSERTL1(it != m_constantMapNameToId.end(), "Cannot find the value for the specified constant: " + valueStr);
698 
699  return std::make_pair(true, m_constant[it->second]);
700  }
701  else if (parserID == AnalyticExpression::numberID)
702  {
703  ASSERTL1(num_children == 0, "Illegal children under number node: " + valueStr);
704  return std::make_pair(true, boost::lexical_cast<NekDouble>(valueStr.c_str()) );
705  }
706  else if (parserID == AnalyticExpression::variableID)
707  {
708  ASSERTL1(num_children == 0, "Illegal children under variable node: " + valueStr);
709 
710  VariableMap::const_iterator it = variableMap.find(valueStr);
711  ASSERTL1(it != variableMap.end(), "Unknown variable parsed: " + valueStr);
712 
713  // Variables are not defined at the time of this parse.
714  stack.push_back ( makeStep<StoreVar>( stateIndex, it->second ) );
715  return std::make_pair(false, 0);
716  }
717  else if (parserID == AnalyticExpression::parameterID)
718  {
719  ASSERTL1(num_children == 0, "Illegal children under parameter node: " + valueStr);
720 
721  ParameterMap::const_iterator it = m_parameterMapNameToId.find(valueStr);
722  ASSERTL1(it != m_parameterMapNameToId.end(), "Unknown parameter parsed: " + valueStr);
723 
724  // Parameters may change in between of evalutions.
725  stack.push_back ( makeStep<StorePrm>( stateIndex, it->second ) );
726  return std::make_pair(false, 0);
727  }
728  else if (parserID == AnalyticExpression::functionID)
729  {
730  FunctionNameMap::const_iterator it = m_functionMapNameToInstanceType.find(valueStr);
731  ASSERTL1(it != m_functionMapNameToInstanceType.end(), "Invalid function specified: " + valueStr);
732  ASSERTL1(num_children == 1 || num_children == 2, "Function " + valueStr + " would like to have too few or too many arguments. This is not implemented yet");
733 
734  if (location->children.size() == 1)
735  {
736  PrecomputedValue v = PrepareExecutionAsYouParse(location->children.begin(), stack, variableMap, stateIndex);
737 
738  // additive white gaussian noise function
739  if (it->second == E_AWGN)
740  {
741  int const_index = AddConstant(std::string("SUB_EXPR_") + boost::lexical_cast<std::string>(m_constant.size()), v.second);
742  stack.push_back ( makeStep<StoreConst>( stateIndex, const_index ) );
743  stack.push_back ( makeStep<EvalAWGN>( stateIndex, stateIndex ) );
744  return std::make_pair(false,0);
745  }
746 
747  if (true == v.first)
748  {
749  return std::make_pair( true, m_function[it->second](v.second) );
750  }
751  }
752  else
753  {
754  PrecomputedValue v1 = PrepareExecutionAsYouParse(location->children.begin()+0, stack, variableMap, stateIndex);
755  PrecomputedValue v2 = PrepareExecutionAsYouParse(location->children.begin()+1, stack, variableMap, stateIndex+1);
756  m_state_size++;
757 
758  if (true == v1.first && true == v2.first)
759  {
760  return std::make_pair( true, m_function2[it->second](v1.second, v2.second) );
761  }
762  }
763 
764  // if somewhere down the parse tree there is a variable or parameter, set up an
765  // evaluation sequence.
766  switch (it->second)
767  {
768  case E_ABS:
769  stack.push_back ( makeStep<EvalAbs>( stateIndex, stateIndex) );
770  return std::make_pair(false,0);
771  case E_ASIN:
772  stack.push_back ( makeStep<EvalAsin>( stateIndex, stateIndex) );
773  return std::make_pair(false,0);
774  case E_ACOS:
775  stack.push_back ( makeStep<EvalAcos>( stateIndex, stateIndex ) );
776  return std::make_pair(false,0);
777  case E_ATAN:
778  stack.push_back ( makeStep<EvalAtan>( stateIndex, stateIndex ) );
779  return std::make_pair(false,0);
780  case E_ATAN2:
781  stack.push_back ( makeStep<EvalAtan2>( stateIndex, stateIndex, stateIndex+1 ) );
782  return std::make_pair(false,0);
783  case E_ANG:
784  stack.push_back ( makeStep<EvalAng>( stateIndex, stateIndex, stateIndex+1 ) );
785  return std::make_pair(false,0);
786  case E_BESSEL:
787  stack.push_back ( makeStep<EvalBessel>( stateIndex, stateIndex, stateIndex+1 ) );
788  return std::make_pair(false,0);
789  case E_CEIL:
790  stack.push_back ( makeStep<EvalCeil>( stateIndex, stateIndex ) );
791  return std::make_pair(false,0);
792  case E_COS:
793  stack.push_back ( makeStep<EvalCos>( stateIndex, stateIndex ) );
794  return std::make_pair(false,0);
795  case E_COSH:
796  stack.push_back ( makeStep<EvalCosh>( stateIndex, stateIndex ) );
797  return std::make_pair(false,0);
798  case E_EXP:
799  stack.push_back ( makeStep<EvalExp>( stateIndex, stateIndex ) );
800  return std::make_pair(false,0);
801  case E_FABS:
802  stack.push_back ( makeStep<EvalFabs>( stateIndex, stateIndex ) );
803  return std::make_pair(false,0);
804  case E_FLOOR:
805  stack.push_back ( makeStep<EvalFloor>( stateIndex, stateIndex ) );
806  return std::make_pair(false,0);
807  case E_LOG:
808  stack.push_back ( makeStep<EvalLog>( stateIndex, stateIndex ) );
809  return std::make_pair(false,0);
810  case E_LOG10:
811  stack.push_back ( makeStep<EvalLog10>( stateIndex, stateIndex ) );
812  return std::make_pair(false,0);
813  case E_RAD:
814  stack.push_back ( makeStep<EvalRad>( stateIndex, stateIndex, stateIndex+1 ) );
815  return std::make_pair(false,0);
816  case E_SIN:
817  stack.push_back ( makeStep<EvalSin>( stateIndex, stateIndex ) );
818  return std::make_pair(false,0);
819  case E_SINH:
820  stack.push_back ( makeStep<EvalSinh>( stateIndex, stateIndex ) );
821  return std::make_pair(false,0);
822  case E_SQRT:
823  stack.push_back ( makeStep<EvalSqrt>( stateIndex, stateIndex ) );
824  return std::make_pair(false,0);
825  case E_TAN:
826  stack.push_back ( makeStep<EvalTan>( stateIndex, stateIndex ) );
827  return std::make_pair(false,0);
828  case E_TANH:
829  stack.push_back ( makeStep<EvalTanh>( stateIndex, stateIndex ) );
830  return std::make_pair(false,0);
831  case E_SIGN:
832  stack.push_back ( makeStep<EvalSign>( stateIndex, stateIndex ) );
833  return std::make_pair(false,0);
834  default:
835  ASSERTL0(false, "Evaluation of " + valueStr + " is not implemented yet");
836  }
837  return std::make_pair(false,0);
838  }
839  else if (parserID == AnalyticExpression::factorID)
840  {
841  ASSERTL1(*valueStr.begin() == '-', "Illegal factor - it can only be '-' and it was: " + valueStr);
842  ASSERTL1(num_children == 1, "Illegal number of children under factor node: " + valueStr);
843 
844  PrecomputedValue v = PrepareExecutionAsYouParse(location->children.begin(), stack, variableMap, stateIndex);
845 
846  // if precomputed value is valid, process it further.
847  if (true == v.first)
848  {
849  return std::make_pair( true, - v.second );
850  }
851  stack.push_back (makeStep<EvalNeg>( stateIndex, stateIndex) );
852  return std::make_pair(false,0);
853  }
854  else if (parserID == AnalyticExpression::operatorID)
855  {
856  ASSERTL1(num_children == 2, "Too few or too many arguments for mathematical operator: " + valueStr);
857  PrecomputedValue left = PrepareExecutionAsYouParse(location->children.begin()+0, stack, variableMap, stateIndex);
858  PrecomputedValue right = PrepareExecutionAsYouParse(location->children.begin()+1, stack, variableMap, stateIndex+1);
859  m_state_size++;
860 
861  // if both precomputed values are valid, process them further.
862  if ((true == left.first) && (true == right.first))
863  {
864  switch(*valueStr.begin())
865  {
866  case '+':
867  return std::make_pair( true, left.second + right.second );
868  case '-':
869  return std::make_pair( true, left.second - right.second );
870  case '*':
871  return std::make_pair( true, left.second * right.second );
872  case '/':
873  return std::make_pair( true, left.second / right.second );
874  case '^':
875  return std::make_pair( true, std::pow(left.second, right.second) );
876  case '=':
877  return std::make_pair( true, left.second == right.second );
878  case '<':
879  if (*(valueStr.end()-1) == '=')
880  {
881  return std::make_pair( true, left.second <= right.second );
882  }
883  else
884  {
885  return std::make_pair( true, left.second < right.second );
886  }
887  return std::make_pair(false,0);
888  case '>':
889  if (*(valueStr.end()-1) == '=')
890  {
891  return std::make_pair( true, left.second >= right.second );
892  }
893  else
894  {
895  return std::make_pair( true, left.second > right.second );
896  }
897  return std::make_pair(false,0);
898  default:
899  ASSERTL0(false, "Invalid operator encountered: " + valueStr);
900  }
901  return std::make_pair(false,0);
902  }
903 
904  // either operator argument is not fully evaluated
905  // add pre-evaluated value to the contaner of constants
906  if (true == left.first)
907  {
908  int const_index = AddConstant(std::string("SUB_EXPR_") + boost::lexical_cast<std::string>(m_constant.size()), left.second);
909  stack.push_back ( makeStep<StoreConst>( stateIndex, const_index ) );
910  }
911  if (true == right.first)
912  {
913  int const_index = AddConstant(std::string("SUB_EXPR_") + boost::lexical_cast<std::string>(m_constant.size()), right.second);
914  stack.push_back ( makeStep<StoreConst>( stateIndex+1, const_index ) );
915  }
916 
917 
918  switch(*valueStr.begin())
919  {
920  case '+':
921  stack.push_back (makeStep<EvalSum>( stateIndex, stateIndex, stateIndex+1 ) );
922  return std::make_pair(false,0);
923  case '-':
924  stack.push_back (makeStep<EvalSub> (stateIndex, stateIndex, stateIndex+1 ) );
925  return std::make_pair(false,0);
926  case '*':
927  stack.push_back (makeStep<EvalMul>( stateIndex, stateIndex, stateIndex+1 ) );
928  return std::make_pair(false,0);
929  case '/':
930  stack.push_back (makeStep<EvalDiv>( stateIndex, stateIndex, stateIndex+1 ) );
931  return std::make_pair(false,0);
932  case '^':
933  stack.push_back (makeStep<EvalPow>( stateIndex, stateIndex, stateIndex+1 ) );
934  return std::make_pair(false,0);
935  case '=':
936  stack.push_back (makeStep<EvalLogicalEqual>( stateIndex, stateIndex, stateIndex+1 ) );
937  return std::make_pair(false,0);
938  case '<':
939  if (*(valueStr.end()-1) == '=')
940  {
941  stack.push_back (makeStep<EvalLogicalLeq>( stateIndex, stateIndex, stateIndex+1 ) );
942  }
943  else
944  {
945  stack.push_back (makeStep<EvalLogicalLess>( stateIndex, stateIndex, stateIndex+1 ) );
946  }
947  return std::make_pair(false,0);
948 
949  case '>':
950  if (*(valueStr.end()-1) == '=')
951  {
952  stack.push_back (makeStep<EvalLogicalGeq>( stateIndex, stateIndex, stateIndex+1 ) );
953  }
954  else
955  {
956  stack.push_back (makeStep<EvalLogicalGreater>( stateIndex, stateIndex, stateIndex+1 ) );
957  }
958  return std::make_pair(false,0);
959 
960  default:
961  ASSERTL0(false, "Invalid operator encountered: " + valueStr);
962  }
963  return std::make_pair(false,0);
964  }
965  ASSERTL0(false, "Illegal expression encountered: " + valueStr);
966  return std::make_pair(false,0);
967  }
968 
969  };
970 };
static NekDouble ang(NekDouble x, NekDouble y)
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:198
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:27
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(* PFD1)(NekDouble)
NekDouble(* PFD4)(NekDouble, NekDouble, NekDouble, NekDouble)
double NekDouble
bg::model::point< double, 3, bg::cs::cartesian > point
Definition: BLMesh.cpp:53
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...
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...
NekDouble(* PFD3)(NekDouble, NekDouble, NekDouble)
NekDouble(* PFD2)(NekDouble, NekDouble)
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
Definition: StdRegions.hpp:316
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:228
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...