32 #include <boost/math/special_functions/gamma.hpp> 33 #include <boost/math/special_functions/polygamma.hpp> 40 void ScalarEvalOp::eval(
double fv[],
size_t n,
const double sv[])
46 assert(
size_t(out)<n);
47 fv[out]=evaluate(0,0);
50 assert(out+in1.size()<=n);
51 for (
unsigned i=0; i<in1.size(); ++i)
52 fv[out+i]=evaluate(flow1? fv[in1[i]]: sv[in1[i]], 0);
55 assert(out+in1.size()<=n);
56 for (
unsigned i=0; i<in1.size(); ++i)
59 const double* v=flow2? fv: sv;
61 x2+=j.weight*v[j.idx];
62 fv[out+i]=evaluate(flow1? fv[in1[i]]: sv[in1[i]], x2);
93 void ScalarEvalOp::deriv(
double df[],
size_t n,
const double ds[],
94 const double sv[],
const double fv[])
96 assert(out>=0 &&
size_t(out)<n);
104 assert((flow1 && in1[0]<ValueVector::flowVars.size()) ||
105 (!flow1 && in1[0]<ValueVector::stockVars.size()));
106 const double x1=flow1? fv[in1[0]]: sv[in1[0]];
107 const double dx1=flow1? df[in1[0]]: ds[in1[0]];
108 df[out] = dx1!=0? dx1 * d1(x1,0): 0;
113 assert((flow1 && in1[0]<ValueVector::flowVars.size()) ||
114 (!flow1 && in1[0]<ValueVector::stockVars.size()));
115 assert((flow2 && in2[0][0].idx<ValueVector::flowVars.size()) ||
116 (!flow2 && in2[0][0].idx<ValueVector::stockVars.size()));
117 const double x1=flow1? fv[in1[0]]: sv[in1[0]];
118 const double x2=flow2? fv[in2[0][0].idx]: sv[in2[0][0].idx];
119 const double dx1=flow1? df[in1[0]]: ds[in1[0]];
120 const double dx2=flow2? df[in2[0][0].idx]: ds[in2[0][0].idx];
121 df[out] = (dx1!=0? dx1 * d1(x1,x2): 0) +
122 (dx2!=0? dx2 *
d2(x1,x2): 0);
127 throw error(
"Invalid operation detected on a %s operation",
128 OperationBase::typeName(type()).c_str());
131 double ConstantEvalOp::evaluate(
double in1,
double in2)
const 143 double EvalOpBase::t;
158 {
return 2.71828182845904523536028747135266249775724709369995;}
168 {
return 3.14159265358979323846264338327950288419716939937510;}
198 {
return numeric_limits<double>::max();}
228 {
throw error(
"evaluate for integral operator should not be called");}
238 {
throw error(
"evaluate for derivative operator should not be called");}
251 {
return state?
dynamic_cast<DataOp&
>(*state).
deriv(x1): 0;}
271 {
return 0.5/
::sqrt(fabs(x1));}
301 {
return 1/(x1*
::log(x2));}
321 {
throw error(
"lt cannot be used with an implicit method");}
331 {
throw error(
"le cannot be used with an implicit method");}
341 {
throw error(
"eq cannot be used with an implicit method");}
351 {
throw error(
"user functions cannot be used with an implicit method");}
354 {
throw error(
"user functions cannot be used with an implicit method");}
359 if (
isnan(in1))
return in2;
360 if (
isnan(in2))
return in1;
361 return std::min(in1,in2);
373 if (
isnan(in1))
return in2;
374 if (
isnan(in2))
return in1;
375 return std::max(in1,in2);
386 {
return in1>0.5 && in2>0.5;}
389 {
throw error(
"and cannot be used with an implicit method");}
396 {
return in1>0.5 || in2>0.5;}
399 {
throw error(
"or cannot be used with an implicit method");}
409 {
throw error(
"not cannot be used with an implicit method");}
439 {
return 1/
sqr(::
cos(x1));}
446 {return ::asin(in1);}
456 {return ::acos(in1);}
466 {return ::atan(in1);}
469 {
return 1/(1+
sqr(x1));}
496 {return ::tanh(in1);}
506 {return ::fabs(in1);}
509 {
return (x1<0)? -1: 1;}
516 {return ::floor(in1);}
519 {
throw error(
"floor cannot be used with an implicit method");}
526 {
return in1-::floor(in1);}
529 {
throw error(
"frac cannot be used with an implicit method");}
536 {
return in1 > 0? boost::math::tgamma(in1) : numeric_limits<double>::max();}
539 {
return boost::math::digamma(fabs(x1))*boost::math::tgamma(fabs(x1));}
556 {
return in1 > -1? boost::math::tgamma(in1+1) : 1;}
602 {
return -x1/(x2*x2);}
607 {
throw error(
"calling evaluate() on EvalOp<numOps> invalid");}
610 {
throw error(
"calling d1 on EvalOp<numOps> invalid");}
613 {
throw error(
"calling d2() on EvalOp<numOps> invalid");}
619 switch (classify(
op))
632 if (
auto f=dynamic_cast<UserFunction*>(state.get()))
651 const std::shared_ptr<VariableValue>& to,
655 auto t=ScalarEvalOp::create(
op,state);
658 assert(t->numArgs()==0 || (from1.
idx()>=0 && (t->numArgs()==1 || from2.
idx()>=0)));
660 if (
auto f=dynamic_cast<UserFunction*>(state.get()))
664 t->in1.push_back(from1.
idx());
668 if (to->idx()==-1) to->allocValue();
Expr polygamma(const Expr &x, const Expr &y)
represents the operation when evaluating the equations
bool isFlowVar() const
returns true if variable's data is allocated on the flowVariables vector
UnitsExpressionWalker pow(const UnitsExpressionWalker &x, const UnitsExpressionWalker &y)
double interpolate(double) const
interpolates y data between x values bounding the argument
std::shared_ptr< Item > ItemPtr
Legacy EvalOp base interface.
Creation and access to the minskyTCL_obj object, which has code to record whenever Minsky's state cha...
double isfinite(double x)
factory class template, for creating objects of type T<O>. N is the maximum value of O ...
double evaluate(double x, double y)
UnitsExpressionWalker timeUnit
legacy data importer object
OperationFactory< ScalarEvalOp, EvalOp, OperationType::sum-1 > evalOpFactory
double deriv(double) const
derivative of the interpolate function. At the data points, the derivative is defined as the weighted...
float d2(float x0, float y0, float x1, float y1)
1.8.14