41 #include <boost/lexical_cast.hpp>
42 #include <boost/algorithm/string/trim.hpp>
43 #include <boost/random/detail/seed.hpp>
46 #include <boost/preprocessor/cat.hpp>
50 #define NEKTAR_MATH_NAME(x) BOOST_PP_CAT(_, x)
52 #define NEKTAR_MATH_NAME(x) x
55 #if( BOOST_VERSION / 100 % 1000 >= 36 )
56 using namespace boost::spirit::classic;
67 namespace LibUtilities
97 NekDouble
sign(NekDouble arg)
99 return (arg > 0.0) - (arg < 0.0);
106 NekDouble
awgn(NekDouble sigma)
109 boost::variate_generator<
111 boost::normal_distribution<>
112 > _normal(rng, boost::normal_distribution<>(0, sigma) );
154 template <
typename ScannerT>
157 expression = logical_or;
159 logical_or = logical_and >> *( (root_node_d[str_p(
"||")] >> logical_and) );
161 logical_and = equality >> *( (root_node_d[str_p(
"&&")] >> equality) );
163 equality = lt_gt >> *( (root_node_d[str_p(
"==")] >> lt_gt) );
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)
172 add_sub = mult_div >>
173 *( (root_node_d[ch_p(
'+')] >> mult_div)
174 | (root_node_d[ch_p(
'-')] >> mult_div)
177 mult_div = exponential >>
178 *( (root_node_d[ch_p(
'*')] >> exponential)
179 | (root_node_d[ch_p(
'/')] >> exponential)
182 exponential = factor >>
183 *( (root_node_d[ch_p(
'^')] >> factor) );
190 | inner_node_d[ch_p(
'(') >> expression >> ch_p(
')')]
191 | (root_node_d[ch_p(
'-')] >> factor);
193 parameter = leaf_node_d[ lexeme_d[
194 (alpha_p |
'_' |
'$') >> *(alnum_p |
'_' |
'$')
198 infix_node_d[inner_node_d[ch_p(
'(') >> expression >> *(
',' >> expression) >> ch_p(
')')]];
200 variable = leaf_node_d[ lexeme_d[
204 number = leaf_node_d[ lexeme_d[
208 constant = leaf_node_d[ lexeme_d[
212 op = eps_p( end_p |
"||" |
"&&" |
"==" |
"<=" |
">=" |
'<' |
'>' |
'+' |
'-' |
'*' |
'/' |
'^' |
')' |
',' );
220 AnalyticExpressionEvaluator::AnalyticExpressionEvaluator():
315 for (std::map<std::string, NekDouble>::const_iterator it = constants.begin(); it != constants.end(); ++it)
337 std::string errormsg(
"Attempt to add numerically different constants under the same name: ");
339 std::cout << errormsg << std::endl;
350 ASSERTL1(value != NULL,
"Constant variable not found: " + name);
357 for (std::map<std::string, NekDouble>::const_iterator it = params.begin(); it != params.end(); it++)
420 std::vector<std::string> variableNames;
421 parse((
char*) vlist.c_str(), ( *space_p >>
423 +(+graph_p)[push_back_a(variableNames)]
433 node_val_data_factory<NekDouble>,
434 std::string::const_iterator,
438 (expr.begin(), expr.end(), myGrammar, space_p);
440 ASSERTL1(parseInfo.full !=
false,
"Unable to fully parse function. Stopped just before: "
441 + std::string(parseInfo.stop, parseInfo.stop + 15));
444 throw std::runtime_error(
"Unable to fully parse function at: "
445 + std::string(parseInfo.stop, parseInfo.stop + 15));
460 for (
int i = 0; i < variableNames.size(); i++)
462 variableMap[variableNames[i]] = i;
473 ASSERTL1(stack.size() == 0,
"Constant expression yeilds non-empty execution stack. Bug in PrepareExecutionAsYouParse()");
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 ) );
495 ASSERTL1(
m_executionStack.size() > expression_id,
"unknown analytic expression, it must first be defined with DefineFunction(...)");
500 for (
int i = 0; i < stack.size(); i++)
502 (*stack[i]).run_once();
512 const int expression_id,
520 ASSERTL1(
m_executionStack.size() > expression_id,
"unknown analytic expression, it must first be defined with DefineFunction(...)");
541 for (
int i = 0; i < stack.size(); i++)
543 (*stack[i]).run_once();
556 ASSERTL1(
m_executionStack.size() > expression_id,
"unknown analytic expression, it must first be defined with DefineFunction(...)");
561 ASSERTL1(point.size() == variableMap.size(),
"The number of variables used to define this expression should match the point dimensionality.");
566 VariableMap::const_iterator it;
567 for (it = variableMap.begin(); it != variableMap.end(); ++it)
572 for (
int i = 0; i < stack.size(); i++)
574 (*stack[i]).run_once();
586 const int expression_id,
595 std::vector<Array<OneD, const NekDouble> >points;
602 Evaluate(expression_id,points,result);
610 const int expression_id,
616 const int num_points = points[0].num_elements();
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");
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);
642 if (result.num_elements() < num_points)
648 int work_left = num_points;
651 const int this_chunk_size = (std::min)(work_left, 1024);
652 for (
int i = 0; i < this_chunk_size; i++)
654 for(
int j = 0; j < nvals; ++j)
656 m_variable[i+this_chunk_size*j] = points[j][offset + i];
659 for (
int i = 0; i < stack.size(); i++)
661 (*stack[i]).run_many(this_chunk_size);
663 for (
int i = 0; i < this_chunk_size; i++)
665 result[offset + i] =
m_state[i];
667 work_left -= this_chunk_size;
668 offset += this_chunk_size;
681 std::string valueStr(location->value.begin(), location->value.end());
682 boost::algorithm::trim(valueStr);
684 const parser_id parserID = location->value.id();
685 #if defined(NEKTAR_DEBUG) || defined(NEKTAR_FULLDEBUG)
686 const int num_children = location->children.size();
691 ASSERTL1(num_children == 0,
"Illegal children under constant node: " + valueStr);
696 return std::make_pair(
true,
m_constant[it->second]);
700 ASSERTL1(num_children == 0,
"Illegal children under number node: " + valueStr);
701 return std::make_pair(
true, boost::lexical_cast<NekDouble>(valueStr.c_str()) );
705 ASSERTL1(num_children == 0,
"Illegal children under variable node: " + valueStr);
707 VariableMap::const_iterator it = variableMap.find(valueStr);
708 ASSERTL1(it != variableMap.end(),
"Unknown variable parsed: " + valueStr);
711 stack.push_back ( makeStep<StoreVar>( stateIndex, it->second ) );
712 return std::make_pair(
false, 0);
716 ASSERTL1(num_children == 0,
"Illegal children under parameter node: " + valueStr);
722 stack.push_back ( makeStep<StorePrm>( stateIndex, it->second ) );
723 return std::make_pair(
false, 0);
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");
731 if (location->children.size() == 1)
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);
746 return std::make_pair(
true,
m_function[it->second](v.second) );
755 if (
true == v1.first &&
true == v2.first)
757 return std::make_pair(
true,
m_function2[it->second](v1.second, v2.second) );
766 stack.push_back ( makeStep<EvalAbs>( stateIndex, stateIndex) );
767 return std::make_pair(
false,0);
769 stack.push_back ( makeStep<EvalAsin>( stateIndex, stateIndex) );
770 return std::make_pair(
false,0);
772 stack.push_back ( makeStep<EvalAcos>( stateIndex, stateIndex ) );
773 return std::make_pair(
false,0);
775 stack.push_back ( makeStep<EvalAtan>( stateIndex, stateIndex ) );
776 return std::make_pair(
false,0);
778 stack.push_back ( makeStep<EvalAtan2>( stateIndex, stateIndex, stateIndex+1 ) );
779 return std::make_pair(
false,0);
781 stack.push_back ( makeStep<EvalAng>( stateIndex, stateIndex, stateIndex+1 ) );
782 return std::make_pair(
false,0);
784 stack.push_back ( makeStep<EvalCeil>( stateIndex, stateIndex ) );
785 return std::make_pair(
false,0);
787 stack.push_back ( makeStep<EvalCos>( stateIndex, stateIndex ) );
788 return std::make_pair(
false,0);
790 stack.push_back ( makeStep<EvalCosh>( stateIndex, stateIndex ) );
791 return std::make_pair(
false,0);
793 stack.push_back ( makeStep<EvalExp>( stateIndex, stateIndex ) );
794 return std::make_pair(
false,0);
796 stack.push_back ( makeStep<EvalFabs>( stateIndex, stateIndex ) );
797 return std::make_pair(
false,0);
799 stack.push_back ( makeStep<EvalFloor>( stateIndex, stateIndex ) );
800 return std::make_pair(
false,0);
802 stack.push_back ( makeStep<EvalLog>( stateIndex, stateIndex ) );
803 return std::make_pair(
false,0);
805 stack.push_back ( makeStep<EvalLog10>( stateIndex, stateIndex ) );
806 return std::make_pair(
false,0);
808 stack.push_back ( makeStep<EvalRad>( stateIndex, stateIndex, stateIndex+1 ) );
809 return std::make_pair(
false,0);
811 stack.push_back ( makeStep<EvalSin>( stateIndex, stateIndex ) );
812 return std::make_pair(
false,0);
814 stack.push_back ( makeStep<EvalSinh>( stateIndex, stateIndex ) );
815 return std::make_pair(
false,0);
817 stack.push_back ( makeStep<EvalSqrt>( stateIndex, stateIndex ) );
818 return std::make_pair(
false,0);
820 stack.push_back ( makeStep<EvalTan>( stateIndex, stateIndex ) );
821 return std::make_pair(
false,0);
823 stack.push_back ( makeStep<EvalTanh>( stateIndex, stateIndex ) );
824 return std::make_pair(
false,0);
826 stack.push_back ( makeStep<EvalSign>( stateIndex, stateIndex ) );
827 return std::make_pair(
false,0);
829 ASSERTL0(
false,
"Evaluation of " + valueStr +
" is not implemented yet");
831 return std::make_pair(
false,0);
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);
843 return std::make_pair(
true, - v.second );
845 stack.push_back (makeStep<EvalNeg>( stateIndex, stateIndex) );
846 return std::make_pair(
false,0);
850 ASSERTL1(num_children == 2,
"Too few or too many arguments for mathematical operator: " + valueStr);
856 if ((
true == left.first) && (
true == right.first))
858 switch(*valueStr.begin())
861 return std::make_pair(
true, left.second + right.second );
863 return std::make_pair(
true, left.second - right.second );
865 return std::make_pair(
true, left.second * right.second );
867 return std::make_pair(
true, left.second / right.second );
869 return std::make_pair(
true, std::pow(left.second, right.second) );
871 return std::make_pair(
true, left.second == right.second );
873 if (*(valueStr.end()-1) ==
'=')
875 return std::make_pair(
true, left.second <= right.second );
879 return std::make_pair(
true, left.second < right.second );
881 return std::make_pair(
false,0);
883 if (*(valueStr.end()-1) ==
'=')
885 return std::make_pair(
true, left.second >= right.second );
889 return std::make_pair(
true, left.second > right.second );
891 return std::make_pair(
false,0);
893 ASSERTL0(
false,
"Invalid operator encountered: " + valueStr);
895 return std::make_pair(
false,0);
900 if (
true == left.first)
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 ) );
905 if (
true == right.first)
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 ) );
912 switch(*valueStr.begin())
915 stack.push_back (makeStep<EvalSum>( stateIndex, stateIndex, stateIndex+1 ) );
916 return std::make_pair(
false,0);
918 stack.push_back (makeStep<EvalSub> (stateIndex, stateIndex, stateIndex+1 ) );
919 return std::make_pair(
false,0);
921 stack.push_back (makeStep<EvalMul>( stateIndex, stateIndex, stateIndex+1 ) );
922 return std::make_pair(
false,0);
924 stack.push_back (makeStep<EvalDiv>( stateIndex, stateIndex, stateIndex+1 ) );
925 return std::make_pair(
false,0);
927 stack.push_back (makeStep<EvalPow>( stateIndex, stateIndex, stateIndex+1 ) );
928 return std::make_pair(
false,0);
930 stack.push_back (makeStep<EvalLogicalEqual>( stateIndex, stateIndex, stateIndex+1 ) );
931 return std::make_pair(
false,0);
933 if (*(valueStr.end()-1) ==
'=')
935 stack.push_back (makeStep<EvalLogicalLeq>( stateIndex, stateIndex, stateIndex+1 ) );
939 stack.push_back (makeStep<EvalLogicalLess>( stateIndex, stateIndex, stateIndex+1 ) );
941 return std::make_pair(
false,0);
944 if (*(valueStr.end()-1) ==
'=')
946 stack.push_back (makeStep<EvalLogicalGeq>( stateIndex, stateIndex, stateIndex+1 ) );
950 stack.push_back (makeStep<EvalLogicalGreater>( stateIndex, stateIndex, stateIndex+1 ) );
952 return std::make_pair(
false,0);
955 ASSERTL0(
false,
"Invalid operator encountered: " + valueStr);
957 return std::make_pair(
false,0);
959 ASSERTL0(
false,
"Illegal expression encountered: " + valueStr);
960 return std::make_pair(
false,0);
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...
FunctionNameMap m_functionMapNameToInstanceType
boost_spirit::tree_parse_info< std::string::const_iterator, boost_spirit::node_val_data_factory< NekDouble > > ParsedTreeInfo
static const int constantID
#define ASSERTL0(condition, msg)
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.
NekDouble m_total_eval_time
std::vector< VariableMap > m_stackVariableMap
Keeping map of variables individually per each analytic expression allows correctly handling expressi...
std::vector< NekDouble > m_parameter
#define sign(a, b)
return the sign(b)*a
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...
std::map< std::string, int > VariableMap
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...
static const int parameterID
static const int variableID
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.
boost_spirit::symbols< NekDouble > m_constantsParser
RandomGeneratorType m_generator
std::map< int, TwoArgFunc > m_function2
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 SetRandomSeed(unsigned int seed=123u)
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)
static const int functionID
std::map< int, OneArgFunc > m_function
Timer m_timer
Timer and sum of evaluation times.
static const int operatorID
static const int factorID
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
std::vector< NekDouble > m_variable
static NekDouble rad(NekDouble x, NekDouble y)
std::pair< bool, NekDouble > PrecomputedValue
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< ExecutionStack > m_executionStack
std::vector< int > m_state_sizes
Vector of state sizes per each.
std::vector< EvaluationStep * > ExecutionStack
std::vector< NekDouble > m_constant
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)
~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.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
ConstantMap m_constantMapNameToId
void SetParameters(std::map< std::string, NekDouble > const ¶ms)
Parameters are like constants, but they are inserted into the function at the time Evaluate is called...
boost::mt19937 RandomGeneratorType
static const int numberID