OilfieldToolsNet/IFEM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

added: template EvalFunction over a Scalar type

Browse Source 
rename to EvalFuncSpatial and make EvalFunction a type alias
for EvalFuncSpatial<Real> to avoid breaking existing code.

add deriv, dderiv, gradient and hessian specializations
for autodiff::var, thus adding automatic differentiation
This commit is contained in:
Arne Morten Kvarving 2023-09-23 18:14:02 +02:00
parent 264146b22e
commit 67a22decb0
3 changed files with 560 additions and 296 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

221
src/Utility/ExprFunctions.C
View File

@ -305,7 +305,9 @@ Real EvalFuncScalar<autodiff::var>::deriv (Real x) const
} }
EvalFunction::EvalFunction (const char* function, Real epsX, Real epsT) template<class Scalar>
EvalFuncSpatial<Scalar>::
EvalFuncSpatial (const char* function, Real epsX, Real epsT)
: dx(epsX), dt(epsT) : dx(epsX), dt(epsT)
{ {
try { try {
@ -354,27 +356,31 @@ EvalFunction::EvalFunction (const char* function, Real epsX, Real epsT)
} }
EvalFunction::~EvalFunction () = default; template<class Scalar>
EvalFuncSpatial<Scalar>::~EvalFuncSpatial () = default;
void EvalFunction::addDerivative (const std::string& function, template<class Scalar>
const std::string& variables, void EvalFuncSpatial<Scalar>::
int d1, int d2) addDerivative (const std::string& function,
const std::string& variables,
int d1, int d2)
{ {
if (d1 > 0 && d1 <= 4 && d2 < 1) // A first order derivative is specified if (d1 > 0 && d1 <= 4 && d2 < 1) // A first order derivative is specified
{ {
if (!derivative1[--d1]) if (!derivative1[--d1])
derivative1[d1] = std::make_unique<EvalFunction>((variables+function).c_str()); derivative1[d1] = std::make_unique<FuncType>((variables+function).c_str());
} }
else if ((d1 = voigtIdx(d1,d2)) >= 0) // A second order derivative is specified else if ((d1 = voigtIdx(d1,d2)) >= 0) // A second order derivative is specified
{ {
if (!derivative2[d1]) if (!derivative2[d1])
derivative2[d1] = std::make_unique<EvalFunction>((variables+function).c_str()); derivative2[d1] = std::make_unique<FuncType>((variables+function).c_str());
} }
} }
Real EvalFunction::evaluate (const Vec3& X) const template<class Scalar>
Real EvalFuncSpatial<Scalar>::evaluate (const Vec3& X) const
{ {
const Vec4* Xt = dynamic_cast<const Vec4*>(&X); const Vec4* Xt = dynamic_cast<const Vec4*>(&X);
Real result = Real(0); Real result = Real(0);
@ -390,7 +396,10 @@ Real EvalFunction::evaluate (const Vec3& X) const
*arg[i].y = X.y; *arg[i].y = X.y;
*arg[i].z = X.z; *arg[i].z = X.z;
*arg[i].t = Xt ? Xt->t : Real(0); *arg[i].t = Xt ? Xt->t : Real(0);
result = expr[i]->Evaluate(); if constexpr (std::is_same_v<Scalar,Real>)
result = expr[i]->Evaluate();
else
result = expr[i]->Evaluate().expr->val;
} }
catch (ExprEval::Exception& e) { catch (ExprEval::Exception& e) {
ExprException(e,"evaluating expression"); ExprException(e,"evaluating expression");
@ -400,7 +409,8 @@ Real EvalFunction::evaluate (const Vec3& X) const
} }
Real EvalFunction::deriv (const Vec3& X, int dir) const template<>
Real EvalFuncSpatial<Real>::deriv (const Vec3& X, int dir) const
{ {
if (dir < 1) if (dir < 1)
return Real(0); return Real(0);
@ -431,7 +441,31 @@ Real EvalFunction::deriv (const Vec3& X, int dir) const
} }
Real EvalFunction::dderiv (const Vec3& X, int d1, int d2) const template<>
Real EvalFuncSpatial<autodiff::var>::deriv (const Vec3& X, int dir) const
{
if (dir < 1 || dir > 4)
return Real(0);
size_t i = 0;
#ifdef USE_OPENMP
i = omp_get_thread_num();
#endif
const Vec4* Xt = dynamic_cast<const Vec4*>(&X);
*arg[i].x = X.x;
*arg[i].y = X.y;
*arg[i].z = X.z;
*arg[i].t = Xt ? Xt->t : Real(0);
// Evaluate spatial derivative using auto-diff
return derivativesx(expr[i]->Evaluate(),
autodiff::wrt(arg[i].get(dir)))[0].expr->val;
}
template<>
Real EvalFuncSpatial<Real>::dderiv (const Vec3& X, int d1, int d2) const
{ {
if ((d1 = voigtIdx(d1,d2)) < 0) if ((d1 = voigtIdx(d1,d2)) < 0)
return Real(0); return Real(0);
@ -440,10 +474,88 @@ Real EvalFunction::dderiv (const Vec3& X, int d1, int d2) const
} }
void EvalFunction::setParam (const std::string& name, double value) template<>
Real EvalFuncSpatial<autodiff::var>::dderiv (const Vec3& X, int d1, int d2) const
{
if (d1 < 1 || d1 > 3 ||
d2 < 1 || d2 > 3)
return Real(0);
size_t i = 0;
#ifdef USE_OPENMP
i = omp_get_thread_num();
#endif
const Vec4* Xt = dynamic_cast<const Vec4*>(&X);
*arg[i].x = X.x;
*arg[i].y = X.y;
*arg[i].z = X.z;
*arg[i].t = Xt ? Xt->t : Real(0);
return derivativesx(derivativesx(expr[i]->Evaluate(),
autodiff::wrt(arg[i].get(d1)))[0],
autodiff::wrt(arg[i].get(d2)))[0].expr->val;
}
template<>
Vec3 EvalFuncSpatial<autodiff::var>::gradient (const Vec3& X) const
{
size_t i = 0;
#ifdef USE_OPENMP
i = omp_get_thread_num();
#endif
const Vec4* Xt = dynamic_cast<const Vec4*>(&X);
double t;
*arg[i].x = X.x;
*arg[i].y = X.y;
*arg[i].z = X.z;
*arg[i].t = t = Xt ? Xt->t : Real(0);
const auto dx = derivativesx(expr[i]->Evaluate(),
autodiff::wrt(*arg[i].x, *arg[i].y, *arg[i].z));
return Vec4(dx[0].expr->val, dx[1].expr->val, dx[2].expr->val, t);
}
template<>
SymmTensor EvalFuncSpatial<autodiff::var>::hessian (const Vec3& X) const
{
size_t i = 0;
#ifdef USE_OPENMP
i = omp_get_thread_num();
#endif
const Vec4* Xt = dynamic_cast<const Vec4*>(&X);
*arg[i].x = X.x;
*arg[i].y = X.y;
*arg[i].z = X.z;
*arg[i].t = Xt ? Xt->t : Real(0);
const auto dx =
derivativesx(expr[i]->Evaluate(), autodiff::wrt(*arg[i].x, *arg[i].y, *arg[i].z));
const auto [uxx, uxy, uxz] =
derivativesx(dx[0], autodiff::wrt(*arg[i].x, *arg[i].y, *arg[i].z));
const auto [uyy, uyz] =
derivativesx(dx[1], autodiff::wrt(*arg[i].y, *arg[i].z));
const auto [uzz] =
derivativesx(dx[2], autodiff::wrt(*arg[i].z));
return SymmTensor({uxx.expr->val, uyy.expr->val, uzz.expr->val,
uxy.expr->val, uyz.expr->val, uxz.expr->val});
}
template<class Scalar>
void EvalFuncSpatial<Scalar>::setParam (const std::string& name, double value)
{ {
for (std::unique_ptr<ValueList>& v1 : v) { for (std::unique_ptr<ValueList>& v1 : v) {
double* address = v1->GetAddress(name); Scalar* address = v1->GetAddress(name);
if (!address) if (!address)
v1->Add(name,value,false); v1->Add(name,value,false);
else else
@ -452,21 +564,24 @@ void EvalFunction::setParam (const std::string& name, double value)
} }
EvalFunctions::EvalFunctions (const std::string& functions, template<class Scalar>
const std::string& variables, EvalFunctions<Scalar>::EvalFunctions (const std::string& functions,
const Real epsX, const Real epsT) const std::string& variables,
const Real epsX, const Real epsT)
{ {
std::vector<std::string> components = splitComps(functions,variables); std::vector<std::string> components = splitComps(functions,variables);
for (const std::string& comp : components) for (const std::string& comp : components)
p.emplace_back(std::make_unique<EvalFunction>(comp.c_str(),epsX,epsT)); p.emplace_back(std::make_unique<FuncType>(comp.c_str(),epsX,epsT));
} }
EvalFunctions::~EvalFunctions () = default; template<class Scalar>
EvalFunctions<Scalar>::~EvalFunctions () = default;
void EvalFunctions::addDerivative (const std::string& functions, template<class Scalar>
const std::string& variables, int d1, int d2) void EvalFunctions<Scalar>::addDerivative (const std::string& functions,
const std::string& variables, int d1, int d2)
{ {
std::vector<std::string> components = splitComps(functions,variables); std::vector<std::string> components = splitComps(functions,variables);
for (size_t i = 0; i < p.size() && i < components.size(); i++) for (size_t i = 0; i < p.size() && i < components.size(); i++)
@ -474,8 +589,9 @@ void EvalFunctions::addDerivative (const std::string& functions,
} }
template <class ParentFunc, class Ret> template <class ParentFunc, class Ret, class Scalar>
Ret EvalMultiFunction<ParentFunc,Ret>::evaluate (const Vec3& X) const Ret EvalMultiFunction<ParentFunc,Ret,Scalar>::
evaluate (const Vec3& X) const
{ {
std::vector<Real> res_array(this->p.size()); std::vector<Real> res_array(this->p.size());
for (size_t i = 0; i < this->p.size(); ++i) for (size_t i = 0; i < this->p.size(); ++i)
@ -485,44 +601,47 @@ Ret EvalMultiFunction<ParentFunc,Ret>::evaluate (const Vec3& X) const
} }
template <class ParentFunc, class Ret> template <class ParentFunc, class Ret, class Scalar>
void EvalMultiFunction<ParentFunc,Ret>::setNoDims () void EvalMultiFunction<ParentFunc,Ret,Scalar>::setNoDims ()
{ {
std::tie(nsd, this->ncmp) = getNoDims<Ret>(this->p.size()); std::tie(nsd, this->ncmp) = getNoDims<Ret>(this->p.size());
} }
template<class ParentFunc, class Ret> template<class ParentFunc, class Ret, class Scalar>
Ret EvalMultiFunction<ParentFunc,Ret>::deriv (const Vec3& X, int dir) const Ret EvalMultiFunction<ParentFunc,Ret,Scalar>::
deriv (const Vec3& X, int dir) const
{ {
std::vector<Real> tmp(p.size()); std::vector<Real> tmp(this->p.size());
for (size_t i = 0; i < p.size(); ++i) for (size_t i = 0; i < this->p.size(); ++i)
tmp[i] = p[i]->deriv(X,dir); tmp[i] = this->p[i]->deriv(X,dir);
return Ret(tmp); return Ret(tmp);
} }
template<class ParentFunc, class Ret> template<class ParentFunc, class Ret, class Scalar>
Ret EvalMultiFunction<ParentFunc,Ret>::dderiv (const Vec3& X, int d1, int d2) const Ret EvalMultiFunction<ParentFunc,Ret,Scalar>::
dderiv (const Vec3& X, int d1, int d2) const
{ {
std::vector<Real> tmp(p.size()); std::vector<Real> tmp(this->p.size());
for (size_t i = 0; i < p.size(); ++i) for (size_t i = 0; i < this->p.size(); ++i)
tmp[i] = p[i]->dderiv(X,d1,d2); tmp[i] = this->p[i]->dderiv(X,d1,d2);
return Ret(tmp); return Ret(tmp);
} }
template <class ParentFunc, class Ret> template <class ParentFunc, class Ret, class Scalar>
std::vector<Real> std::vector<Real>
EvalMultiFunction<ParentFunc, Ret>::evalGradient (const Vec3& X) const EvalMultiFunction<ParentFunc,Ret,Scalar>::
evalGradient (const Vec3& X) const
{ {
std::vector<Real> result; std::vector<Real> result;
result.reserve(this->ncmp*this->nsd); result.reserve(this->ncmp*this->nsd);
std::vector<Vec3> dx; std::vector<Vec3> dx;
dx.reserve(this->p.size()); dx.reserve(this->p.size());
for (const std::unique_ptr<EvalFunction>& f : this->p) for (const std::unique_ptr<FuncType>& f : this->p)
dx.push_back(f->gradient(X)); dx.push_back(f->gradient(X));
for (size_t d = 1; d <= this->nsd; ++d) for (size_t d = 1; d <= this->nsd; ++d)
@ -533,15 +652,16 @@ EvalMultiFunction<ParentFunc, Ret>::evalGradient (const Vec3& X) const
} }
template <class ParentFunc, class Ret> template <class ParentFunc, class Ret, class Scalar>
std::vector<Real> std::vector<Real>
EvalMultiFunction<ParentFunc,Ret>::evalHessian (const Vec3& X) const EvalMultiFunction<ParentFunc,Ret,Scalar>::
evalHessian (const Vec3& X) const
{ {
std::vector<Real> result; std::vector<Real> result;
result.reserve(this->p.size()*this->nsd*this->nsd); result.reserve(this->p.size()*this->nsd*this->nsd);
std::vector<SymmTensor> dx; std::vector<SymmTensor> dx;
dx.reserve(this->p.size()); dx.reserve(this->p.size());
for (const std::unique_ptr<EvalFunction>& f : this->p) for (const std::unique_ptr<FuncType>& f : this->p)
dx.push_back(f->hessian(X)); dx.push_back(f->hessian(X));
for (size_t d2 = 1; d2 <= this->nsd; ++d2) for (size_t d2 = 1; d2 <= this->nsd; ++d2)
@ -553,14 +673,14 @@ EvalMultiFunction<ParentFunc,Ret>::evalHessian (const Vec3& X) const
} }
template <class ParentFunc, class Ret> template <class ParentFunc, class Ret, class Scalar>
std::vector<Real> std::vector<Real>
EvalMultiFunction<ParentFunc, Ret>::evalTimeDerivative (const Vec3& X) const EvalMultiFunction<ParentFunc,Ret,Scalar>::evalTimeDerivative (const Vec3& X) const
{ {
std::vector<Real> result; std::vector<Real> result;
result.reserve(this->ncmp); result.reserve(this->ncmp);
for (const std::unique_ptr<EvalFunction>& f : this->p) for (const std::unique_ptr<FuncType>& f : this->p)
result.push_back(f->deriv(X,4)); result.push_back(f->timeDerivative(X));
return result; return result;
} }
@ -568,6 +688,13 @@ EvalMultiFunction<ParentFunc, Ret>::evalTimeDerivative (const Vec3& X) const
template class EvalFuncScalar<Real>; template class EvalFuncScalar<Real>;
template class EvalFuncScalar<autodiff::var>; template class EvalFuncScalar<autodiff::var>;
template class EvalMultiFunction<VecFunc,Vec3>; template class EvalFuncSpatial<Real>;
template class EvalMultiFunction<TensorFunc,Tensor>; template class EvalFuncSpatial<autodiff::var>;
template class EvalMultiFunction<STensorFunc,SymmTensor>; template class EvalFunctions<Real>;
template class EvalFunctions<autodiff::var>;
template class EvalMultiFunction<VecFunc,Vec3,Real>;
template class EvalMultiFunction<VecFunc,Vec3,autodiff::var>;
template class EvalMultiFunction<TensorFunc,Tensor,Real>;
template class EvalMultiFunction<TensorFunc,Tensor,autodiff::var>;
template class EvalMultiFunction<STensorFunc,SymmTensor,Real>;
template class EvalMultiFunction<STensorFunc,SymmTensor,autodiff::var>;

79
src/Utility/ExprFunctions.h
View File

@ -80,11 +80,13 @@ protected:
\brief A scalar-valued spatial function, general function expression. \brief A scalar-valued spatial function, general function expression.
*/ */
class EvalFunction : public RealFunc template<class Scalar>
class EvalFuncSpatial : public RealFunc
{ {
using Expression = ExprEval::Expression<Real>; //!< Type alias for expression tree using Expression = ExprEval::Expression<Scalar>; //!< Type alias for expression tree
using FunctionList = ExprEval::FunctionList<Real>; //!< Type alias for function list using FunctionList = ExprEval::FunctionList<Scalar>; //!< Type alias for function list
using ValueList = ExprEval::ValueList<Real>; //!< Type alias for value list using ValueList = ExprEval::ValueList<Scalar>; //!< Type alias for value list
using FuncType = EvalFuncSpatial<Scalar>; //!< Type alias for function
std::vector<std::unique_ptr<Expression>> expr; //!< Roots of the expression tree std::vector<std::unique_ptr<Expression>> expr; //!< Roots of the expression tree
std::vector<std::unique_ptr<FunctionList>> f; //!< Lists of functions std::vector<std::unique_ptr<FunctionList>> f; //!< Lists of functions
std::vector<std::unique_ptr<ValueList>> v; //!< Lists of variables and constants std::vector<std::unique_ptr<ValueList>> v; //!< Lists of variables and constants
@ -92,16 +94,27 @@ class EvalFunction : public RealFunc
//! \brief A struct representing a spatial function argument. //! \brief A struct representing a spatial function argument.
struct Arg struct Arg
{ {
Real* x; //!< X-coordinate Scalar* x; //!< X-coordinate
Real* y; //!< Y-coordinate Scalar* y; //!< Y-coordinate
Real* z; //!< Z-coordinate Scalar* z; //!< Z-coordinate
Real* t; //!< Time Scalar* t; //!< Time
const Scalar& get(int dir) const
{
switch (dir) {
default:
case 1: return *x;
case 2: return *y;
case 3: return *z;
case 4: return *t;
}
}
}; };
std::vector<Arg> arg; //!< Function argument values std::vector<Arg> arg; //!< Function argument values
std::array<std::unique_ptr<EvalFunction>,4> derivative1; //!< First order derivative expressions std::array<std::unique_ptr<FuncType>,4> derivative1; //!< First order derivative expressions
std::array<std::unique_ptr<EvalFunction>,6> derivative2; //!< Second order derivative expressions std::array<std::unique_ptr<FuncType>,6> derivative2; //!< Second order derivative expressions
bool IAmConstant; //!< Indicates whether the time coordinate is given or not bool IAmConstant; //!< Indicates whether the time coordinate is given or not
@ -110,12 +123,12 @@ class EvalFunction : public RealFunc
public: public:
//! \brief The constructor parses the expression string. //! \brief The constructor parses the expression string.
explicit EvalFunction(const char* function, explicit EvalFuncSpatial(const char* function,
Real epsX = Real(1.0e-8), Real epsT = Real(1.0e-12)); Real epsX = Real(1.0e-8), Real epsT = Real(1.0e-12));
//! \brief Defaulted destructor. //! \brief Defaulted destructor.
//! \details The implementation needs to be in compile unit so we have the //! \details The implementation needs to be in compile unit so we have the
//! definition for the types of the unique_ptr's. //! definition for the types of the unique_ptr's.
virtual ~EvalFunction(); virtual ~EvalFuncSpatial();
//! \brief Adds an expression function for a first or second derivative. //! \brief Adds an expression function for a first or second derivative.
void addDerivative(const std::string& function, const std::string& variables, void addDerivative(const std::string& function, const std::string& variables,
@ -132,11 +145,19 @@ public:
//! \brief Set an additional parameter in the function. //! \brief Set an additional parameter in the function.
void setParam(const std::string& name, double value); void setParam(const std::string& name, double value);
//! \brief Evaluates first derivatives of the function.
Vec3 gradient(const Vec3& X) const override
{
return this->RealFunc::gradient(X);
}
//! \brief Evaluates first derivatives of the function.
SymmTensor hessian(const Vec3& X) const override
{
return this->RealFunc::hessian(X);
}
protected: protected:
//! \brief Non-implemented copy constructor to disallow copying.
EvalFunction(const EvalFunction&) = delete;
//! \brief Non-implemented assignment operator to disallow copying.
EvalFunction& operator=(const EvalFunction&) = delete;
//! \brief Evaluates the function expression. //! \brief Evaluates the function expression.
Real evaluate(const Vec3& X) const override; Real evaluate(const Vec3& X) const override;
}; };
@ -146,9 +167,12 @@ protected:
\brief A base class for multi-component expression functions. \brief A base class for multi-component expression functions.
*/ */
template<class Scalar>
class EvalFunctions class EvalFunctions
{ {
protected: protected:
using FuncType = EvalFuncSpatial<Scalar>; //!< Type alias for function
//! \brief The constructor parses the expression string for each component. //! \brief The constructor parses the expression string for each component.
EvalFunctions(const std::string& functions, const std::string& variables, EvalFunctions(const std::string& functions, const std::string& variables,
const Real epsX, const Real epsT); const Real epsX, const Real epsT);
@ -163,7 +187,7 @@ public:
const std::string& variables, int d1, int d2 = 0); const std::string& variables, int d1, int d2 = 0);
protected: protected:
std::vector<std::unique_ptr<EvalFunction>> p; //!< Array of component expressions std::vector<std::unique_ptr<FuncType>> p; //!< Array of component expressions
}; };
@ -172,10 +196,11 @@ protected:
\details The function is implemented as an array of EvalFunction objects. \details The function is implemented as an array of EvalFunction objects.
*/ */
template <class ParentFunc, class Ret> template <class ParentFunc, class Ret, class Scalar>
class EvalMultiFunction : public ParentFunc, public EvalFunctions class EvalMultiFunction : public ParentFunc, public EvalFunctions<Scalar>
{ {
size_t nsd; //!< Number of spatial dimensions size_t nsd; //!< Number of spatial dimensions
using FuncType = typename EvalFunctions<Scalar>::FuncType; //!< Type alias for function
public: public:
//! \brief The constructor parses the expression string for each component. //! \brief The constructor parses the expression string for each component.
@ -183,7 +208,7 @@ public:
const std::string& variables = "", const std::string& variables = "",
const Real epsX = 1e-8, const Real epsX = 1e-8,
const Real epsT = 1e-12) const Real epsT = 1e-12)
: EvalFunctions(functions,variables,epsX,epsT), nsd(0) { this->setNoDims(); } : EvalFunctions<Scalar>(functions,variables,epsX,epsT), nsd(0) { this->setNoDims(); }
//! \brief Empty destructor. //! \brief Empty destructor.
virtual ~EvalMultiFunction() {} virtual ~EvalMultiFunction() {}
@ -191,7 +216,7 @@ public:
//! \brief Returns whether the function is time-independent or not. //! \brief Returns whether the function is time-independent or not.
bool isConstant() const override bool isConstant() const override
{ {
for (const std::unique_ptr<EvalFunction>& func : p) for (const std::unique_ptr<FuncType>& func : this->p)
if (!func->isConstant()) if (!func->isConstant())
return false; return false;
return true; return true;
@ -208,7 +233,7 @@ public:
//! \brief Set an additional parameter in the function. //! \brief Set an additional parameter in the function.
void setParam(const std::string& name, double value) void setParam(const std::string& name, double value)
{ {
for (std::unique_ptr<EvalFunction>& func : p) for (std::unique_ptr<FuncType>& func : this->p)
func->setParam(name, value); func->setParam(name, value);
} }
@ -234,11 +259,13 @@ protected:
//! Scalar-valued function expression //! Scalar-valued function expression
using EvalFunc = EvalFuncScalar<Real>; using EvalFunc = EvalFuncScalar<Real>;
//! Scalar-valued spatial function expression
using EvalFunction = EvalFuncSpatial<Real>;
//! Vector-valued function expression //! Vector-valued function expression
using VecFuncExpr = EvalMultiFunction<VecFunc,Vec3>; using VecFuncExpr = EvalMultiFunction<VecFunc,Vec3,Real>;
//! Tensor-valued function expression //! Tensor-valued function expression
using TensorFuncExpr = EvalMultiFunction<TensorFunc,Tensor>; using TensorFuncExpr = EvalMultiFunction<TensorFunc,Tensor,Real>;
//! Symmetric tensor-valued function expression //! Symmetric tensor-valued function expression
using STensorFuncExpr = EvalMultiFunction<STensorFunc,SymmTensor>; using STensorFuncExpr = EvalMultiFunction<STensorFunc,SymmTensor,Real>;
#endif #endif

556
src/Utility/Test/TestFunctions.C
View File

@ -23,6 +23,13 @@
#include <cmath> #include <cmath>
using EvalFuncAd = EvalFuncScalar<autodiff::var>;
using EvalFunctionAd = EvalFuncSpatial<autodiff::var>;
using VecFuncExprAd = EvalMultiFunction<VecFunc,Vec3,autodiff::var>;
using TensorFuncExprAd = EvalMultiFunction<TensorFunc,Tensor,autodiff::var>;
using STensorFuncExprAd = EvalMultiFunction<STensorFunc,SymmTensor,autodiff::var>;
TEST(TestScalarFunc, ParseDerivative) TEST(TestScalarFunc, ParseDerivative)
{ {
const char* func1 = "sin(1.5*t)*t"; const char* func1 = "sin(1.5*t)*t";
@ -31,7 +38,7 @@ TEST(TestScalarFunc, ParseDerivative)
ScalarFunc* f1 = utl::parseTimeFunc(func1,"expression"); ScalarFunc* f1 = utl::parseTimeFunc(func1,"expression");
ScalarFunc* f2 = utl::parseTimeFunc(func2,"expression"); ScalarFunc* f2 = utl::parseTimeFunc(func2,"expression");
EvalFuncScalar<autodiff::var> f3(func1,"t"); EvalFuncAd f3(func1,"t");
ASSERT_TRUE(f1 != nullptr); ASSERT_TRUE(f1 != nullptr);
ASSERT_TRUE(f2 != nullptr); ASSERT_TRUE(f2 != nullptr);
@ -56,17 +63,20 @@ TEST(TestScalarFunc, ParseDerivative)
TEST(TestRealFunc, Gradient) TEST(TestRealFunc, Gradient)
{ {
const char* f1 = "sin(x)*sin(y)*sin(z)"; const char* g = "sin(x)*sin(y)*sin(z)";
const char* d1 = "cos(x)*sin(y)*sin(z)"; const char* g_x = "cos(x)*sin(y)*sin(z)";
const char* d2 = "sin(x)*cos(y)*sin(z)"; const char* g_y = "sin(x)*cos(y)*sin(z)";
const char* d3 = "sin(x)*sin(y)*cos(z)"; const char* g_z = "sin(x)*sin(y)*cos(z)";
EvalFunction f(f1); EvalFunction f1(g);
f.addDerivative(d1, "", 1); f1.addDerivative(g_x, "", 1);
f.addDerivative(d2, "", 2); f1.addDerivative(g_y, "", 2);
f.addDerivative(d3, "", 3); f1.addDerivative(g_z, "", 3);
EXPECT_TRUE(f.isConstant()); EvalFunctionAd f2(g);
EXPECT_TRUE(f1.isConstant());
EXPECT_TRUE(f2.isConstant());
for (double x : {0.1, 0.2, 0.3}) for (double x : {0.1, 0.2, 0.3})
for (double y : {0.5, 0.6, 0.7}) for (double y : {0.5, 0.6, 0.7})
@ -76,10 +86,13 @@ TEST(TestRealFunc, Gradient)
sin(x)*cos(y)*sin(z), sin(x)*cos(y)*sin(z),
sin(x)*sin(y)*cos(z)); sin(x)*sin(y)*cos(z));
const Vec3 grad = f.gradient(X); for (const RealFunc* fp : {static_cast<const RealFunc*>(&f1),
for (size_t i = 1; i <= 3; ++i) { static_cast<const RealFunc*>(&f2)}) {
EXPECT_DOUBLE_EQ(f.deriv(X, i), r[i-1]); const Vec3 grad = fp->gradient(X);
EXPECT_DOUBLE_EQ(grad[i-1], r[i-1]); for (size_t i = 1; i <= 3; ++i) {
EXPECT_DOUBLE_EQ(fp->deriv(X, i), r[i-1]);
EXPECT_DOUBLE_EQ(grad[i-1], r[i-1]);
}
} }
} }
} }
@ -116,23 +129,26 @@ TEST(TestRealFunc, GradientFD)
TEST(TestRealFunc, Hessian) TEST(TestRealFunc, Hessian)
{ {
const char* f1 = "sin(x)*sin(y)*sin(z)"; const char* g = "sin(x)*sin(y)*sin(z)";
const char* d11 = "-sin(x)*sin(y)*sin(z)"; const char* g_xx = "-sin(x)*sin(y)*sin(z)";
const char* d22 = "-sin(x)*sin(y)*sin(z)"; const char* g_xy = "cos(x)*cos(y)*sin(z)";
const char* d33 = "-sin(x)*sin(y)*sin(z)"; const char* g_xz = "cos(x)*sin(y)*cos(z)";
const char* d12 = "cos(x)*cos(y)*sin(z)"; const char* g_yy = "-sin(x)*sin(y)*sin(z)";
const char* d13 = "cos(x)*sin(y)*cos(z)"; const char* g_zz = "-sin(x)*sin(y)*sin(z)";
const char* d23 = "sin(x)*cos(y)*cos(z)"; const char* g_yz = "sin(x)*cos(y)*cos(z)";
EvalFunction f(f1); EvalFunction f1(g);
f.addDerivative(d11,"",1,1); f1.addDerivative(g_xx,"",1,1);
f.addDerivative(d12,"",1,2); f1.addDerivative(g_xy,"",1,2);
f.addDerivative(d13,"",1,3); f1.addDerivative(g_xz,"",1,3);
f.addDerivative(d22,"",2,2); f1.addDerivative(g_yy,"",2,2);
f.addDerivative(d23,"",2,3); f1.addDerivative(g_yz,"",2,3);
f.addDerivative(d33,"",3,3); f1.addDerivative(g_zz,"",3,3);
EXPECT_TRUE(f.isConstant()); EvalFunctionAd f2(g);
EXPECT_TRUE(f1.isConstant());
EXPECT_TRUE(f2.isConstant());
for (double x : {0.1, 0.2, 0.3}) for (double x : {0.1, 0.2, 0.3})
for (double y : {0.5, 0.6, 0.7}) for (double y : {0.5, 0.6, 0.7})
@ -145,11 +161,14 @@ TEST(TestRealFunc, Hessian)
sin(x)*cos(y)*cos(z), sin(x)*cos(y)*cos(z),
cos(x)*sin(y)*cos(z)}); cos(x)*sin(y)*cos(z)});
const SymmTensor hess = f.hessian(X); for (const RealFunc* fp : {static_cast<const RealFunc*>(&f1),
for (size_t i = 1; i <= 3; ++i) static_cast<const RealFunc*>(&f2)}) {
for (size_t j = 1; j <= 3; ++j) { const SymmTensor hess = fp->hessian(X);
EXPECT_DOUBLE_EQ(f.dderiv(X,i,j), r(i,j)); for (size_t i = 1; i <= 3; ++i)
EXPECT_DOUBLE_EQ(hess(i,j), r(i,j)); for (size_t j = 1; j <= 3; ++j) {
EXPECT_DOUBLE_EQ(fp->dderiv(X,i,j), r(i,j));
EXPECT_DOUBLE_EQ(hess(i,j), r(i,j));
}
} }
} }
} }
@ -159,17 +178,21 @@ TEST(TestVecFunc, Evaluate)
{ {
const char* func = "sin(x) | cos (y) | exp(z)"; const char* func = "sin(x) | cos (y) | exp(z)";
VecFuncExpr f(func); VecFuncExpr f1(func);
VecFuncExprAd f2(func);
for (double x : {0.1, 0.2, 0.3}) for (double x : {0.1, 0.2, 0.3})
for (double y : {0.5, 0.6, 0.7}) for (double y : {0.5, 0.6, 0.7})
for (double z : {0.8, 0.9, 1.0}) { for (double z : {0.8, 0.9, 1.0}) {
const Vec3 X(x,y,z); const Vec3 X(x,y,z);
const Vec3 fx = f(X);
const Vec3 r(sin(x), cos(y), exp(z)); const Vec3 r(sin(x), cos(y), exp(z));
EXPECT_DOUBLE_EQ(fx.x, r.x); for (const VecFunc* fp : {static_cast<const VecFunc*>(&f1),
EXPECT_DOUBLE_EQ(fx.y, r.y); static_cast<const VecFunc*>(&f2)}) {
EXPECT_DOUBLE_EQ(fx.z, r.z); const Vec3 fx = (*fp)(X);
EXPECT_DOUBLE_EQ(fx.x, r.x);
EXPECT_DOUBLE_EQ(fx.y, r.y);
EXPECT_DOUBLE_EQ(fx.z, r.z);
}
} }
} }
@ -180,11 +203,14 @@ TEST(TestVecFunction, Gradient2D)
const char* g_x = "cos(x)*sin(y) | 2*x*y*y"; const char* g_x = "cos(x)*sin(y) | 2*x*y*y";
const char* g_y = "sin(x)*cos(y) | 2*x*x*y"; const char* g_y = "sin(x)*cos(y) | 2*x*x*y";
VecFuncExpr f(g); VecFuncExpr f1(g);
f.addDerivative(g_x,"",1); f1.addDerivative(g_x,"",1);
f.addDerivative(g_y,"",2); f1.addDerivative(g_y,"",2);
EXPECT_TRUE(f.isConstant()); VecFuncExprAd f2(g);
EXPECT_TRUE(f1.isConstant());
EXPECT_TRUE(f2.isConstant());
for (double x : {0.1, 0.2, 0.3}) for (double x : {0.1, 0.2, 0.3})
for (double y : {0.5, 0.6, 0.7}) { for (double y : {0.5, 0.6, 0.7}) {
@ -192,12 +218,15 @@ TEST(TestVecFunction, Gradient2D)
const Tensor r({cos(x)*sin(y), 2*x*y*y, const Tensor r({cos(x)*sin(y), 2*x*y*y,
sin(x)*cos(y), 2*x*x*y}); sin(x)*cos(y), 2*x*x*y});
const Tensor grad = f.gradient(X); for (const VecFunc* fp : {static_cast<const VecFunc*>(&f1),
for (size_t d = 1; d <= 2; ++d) { static_cast<const VecFunc*>(&f2)}) {
const Vec3 dx = f.deriv(X,d); const Tensor grad = fp->gradient(X);
for (size_t i = 1; i <= 2; ++i) { for (size_t d = 1; d <= 2; ++d) {
EXPECT_DOUBLE_EQ(dx[i-1], r(i,d)); const Vec3 dx = fp->deriv(X,d);
EXPECT_DOUBLE_EQ(grad(i,d), r(i,d)); for (size_t i = 1; i <= 2; ++i) {
EXPECT_DOUBLE_EQ(dx[i-1], r(i,d));
EXPECT_DOUBLE_EQ(grad(i,d), r(i,d));
}
} }
} }
} }
@ -243,12 +272,15 @@ TEST(TestVecFunction, Gradient3D)
const char* g_y = "sin(x)*cos(y)*sin(z) | 2*x*x*y*z*z*z | 2*exp(-x)*exp(2*y)*exp(z*z)"; const char* g_y = "sin(x)*cos(y)*sin(z) | 2*x*x*y*z*z*z | 2*exp(-x)*exp(2*y)*exp(z*z)";
const char* g_z = "sin(x)*sin(y)*cos(z) | 3*x*x*y*y*z*z | 2*z*exp(-x)*exp(2*y)*exp(z*z)"; const char* g_z = "sin(x)*sin(y)*cos(z) | 3*x*x*y*y*z*z | 2*z*exp(-x)*exp(2*y)*exp(z*z)";
VecFuncExpr f(g); VecFuncExpr f1(g);
f.addDerivative(g_x,"",1); f1.addDerivative(g_x,"",1);
f.addDerivative(g_y,"",2); f1.addDerivative(g_y,"",2);
f.addDerivative(g_z,"",3); f1.addDerivative(g_z,"",3);
EXPECT_TRUE(f.isConstant()); VecFuncExprAd f2(g);
EXPECT_TRUE(f1.isConstant());
EXPECT_TRUE(f2.isConstant());
for (double x : {0.1, 0.2, 0.3}) for (double x : {0.1, 0.2, 0.3})
for (double y : {0.5, 0.6, 0.7}) for (double y : {0.5, 0.6, 0.7})
@ -258,12 +290,15 @@ TEST(TestVecFunction, Gradient3D)
sin(x)*cos(y)*sin(z), 2*x*x*y*z*z*z, 2*exp(-x)*exp(2*y)*exp(z*z), sin(x)*cos(y)*sin(z), 2*x*x*y*z*z*z, 2*exp(-x)*exp(2*y)*exp(z*z),
sin(x)*sin(y)*cos(z), 3*x*x*y*y*z*z, 2*z*exp(-x)*exp(2*y)*exp(z*z)}); sin(x)*sin(y)*cos(z), 3*x*x*y*y*z*z, 2*z*exp(-x)*exp(2*y)*exp(z*z)});
const Tensor grad = f.gradient(X); for (const VecFunc* fp : {static_cast<const VecFunc*>(&f1),
for (size_t d = 1; d <= 3; ++d) { static_cast<const VecFunc*>(&f2)}) {
const Vec3 dx = f.deriv(X,d); const Tensor grad = fp->gradient(X);
for (size_t i = 1; i <= 3; ++i) { for (size_t d = 1; d <= 3; ++d) {
EXPECT_DOUBLE_EQ(dx[i-1], r(i,d)); const Vec3 dx = fp->deriv(X,d);
EXPECT_DOUBLE_EQ(grad(i,d), r(i,d)); for (size_t i = 1; i <= 3; ++i) {
EXPECT_DOUBLE_EQ(dx[i-1], r(i,d));
EXPECT_DOUBLE_EQ(grad(i,d), r(i,d));
}
} }
} }
} }
@ -314,12 +349,15 @@ TEST(TestVecFunction, Hessian2D)
const char* g_xy = "cos(x)*cos(y) | 2*x*2*y"; const char* g_xy = "cos(x)*cos(y) | 2*x*2*y";
const char* g_yy = "-sin(x)*sin(y) | 2*x*x"; const char* g_yy = "-sin(x)*sin(y) | 2*x*x";
VecFuncExpr f(g); VecFuncExpr f1(g);
f.addDerivative(g_xx,"",1,1); f1.addDerivative(g_xx,"",1,1);
f.addDerivative(g_yy,"",2,2); f1.addDerivative(g_yy,"",2,2);
f.addDerivative(g_xy,"",1,2); f1.addDerivative(g_xy,"",1,2);
EXPECT_TRUE(f.isConstant()); VecFuncExprAd f2(g);
EXPECT_TRUE(f1.isConstant());
EXPECT_TRUE(f2.isConstant());
for (double x : {0.1, 0.2, 0.3}) for (double x : {0.1, 0.2, 0.3})
for (double y : {0.5, 0.6, 0.7}) { for (double y : {0.5, 0.6, 0.7}) {
@ -330,15 +368,18 @@ TEST(TestVecFunction, Hessian2D)
cos(x)*cos(y), 2*x*2*y, cos(x)*cos(y), 2*x*2*y,
-sin(x)*sin(y), 2*x*x}.data()); -sin(x)*sin(y), 2*x*x}.data());
const utl::matrix3d<Real> hess = f.hessian(X); for (const VecFunc* fp : {static_cast<const VecFunc*>(&f1),
for (size_t d1 = 1; d1 <= 2; ++d1) static_cast<const VecFunc*>(&f2)}) {
for (size_t d2 = 1; d2 <= 2; ++d2) { const utl::matrix3d<Real> hess = fp->hessian(X);
const Vec3 d2x = f.dderiv(X,d1,d2); for (size_t d1 = 1; d1 <= 2; ++d1)
for (size_t i = 1; i <= 2; ++i) { for (size_t d2 = 1; d2 <= 2; ++d2) {
EXPECT_DOUBLE_EQ(hess(i,d1,d2), r(i,d1,d2)); const Vec3 d2x = fp->dderiv(X,d1,d2);
EXPECT_DOUBLE_EQ(d2x[i-1], r(i,d1,d2)); for (size_t i = 1; i <= 2; ++i) {
EXPECT_DOUBLE_EQ(hess(i,d1,d2), r(i,d1,d2));
EXPECT_DOUBLE_EQ(d2x[i-1], r(i,d1,d2));
}
} }
} }
} }
} }
@ -353,15 +394,18 @@ TEST(TestVecFunction, Hessian3D)
const char* g_yz = "sin(x)*cos(y)*cos(z) | x*x*2*y*2*z | exp(x)*2*exp(2*y)*3*exp(3*z)"; const char* g_yz = "sin(x)*cos(y)*cos(z) | x*x*2*y*2*z | exp(x)*2*exp(2*y)*3*exp(3*z)";
const char* g_zz = "-sin(x)*sin(y)*sin(z) | x*x*y*y*2 | exp(x)*exp(2*y)*9*exp(3*z)"; const char* g_zz = "-sin(x)*sin(y)*sin(z) | x*x*y*y*2 | exp(x)*exp(2*y)*9*exp(3*z)";
VecFuncExpr f(g); VecFuncExpr f1(g);
f.addDerivative(g_xx,"",1,1); f1.addDerivative(g_xx,"",1,1);
f.addDerivative(g_yy,"",2,2); f1.addDerivative(g_yy,"",2,2);
f.addDerivative(g_zz,"",3,3); f1.addDerivative(g_zz,"",3,3);
f.addDerivative(g_xy,"",1,2); f1.addDerivative(g_xy,"",1,2);
f.addDerivative(g_xz,"",1,3); f1.addDerivative(g_xz,"",1,3);
f.addDerivative(g_yz,"",2,3); f1.addDerivative(g_yz,"",2,3);
EXPECT_TRUE(f.isConstant()); VecFuncExprAd f2(g);
EXPECT_TRUE(f1.isConstant());
EXPECT_TRUE(f2.isConstant());
for (double x : {0.1, 0.2, 0.3}) for (double x : {0.1, 0.2, 0.3})
for (double y : {0.5, 0.6, 0.7}) for (double y : {0.5, 0.6, 0.7})
@ -380,15 +424,19 @@ TEST(TestVecFunction, Hessian3D)
sin(x)*cos(y)*cos(z), x*x*2*y*2*z, exp(x)*2*exp(2*y)*3*exp(3*z), sin(x)*cos(y)*cos(z), x*x*2*y*2*z, exp(x)*2*exp(2*y)*3*exp(3*z),
-sin(x)*sin(y)*sin(z), x*x*y*y*2, exp(x)*exp(2*y)*9*exp(3*z)}.data()); -sin(x)*sin(y)*sin(z), x*x*y*y*2, exp(x)*exp(2*y)*9*exp(3*z)}.data());
const utl::matrix3d<Real> hess = f.hessian(X);
for (size_t d1 = 1; d1 <= 3; ++d1) for (const VecFunc* fp : {static_cast<const VecFunc*>(&f1),
for (size_t d2 = 1; d2 <= 3; ++d2) { static_cast<const VecFunc*>(&f2)}) {
const Vec3 d2x = f.dderiv(X,d1,d2); const utl::matrix3d<Real> hess = fp->hessian(X);
for (size_t i = 1; i <= 3; ++i) { for (size_t d1 = 1; d1 <= 3; ++d1)
EXPECT_DOUBLE_EQ(d2x[i-1], r(i,d1,d2)); for (size_t d2 = 1; d2 <= 3; ++d2) {
EXPECT_DOUBLE_EQ(hess(i,d1,d2), r(i,d1,d2)); const Vec3 d2x = fp->dderiv(X,d1,d2);
for (size_t i = 1; i <= 3; ++i) {
EXPECT_DOUBLE_EQ(d2x[i-1], r(i,d1,d2));
EXPECT_DOUBLE_EQ(hess(i,d1,d2), r(i,d1,d2));
}
} }
} }
} }
} }
@ -470,17 +518,22 @@ TEST(TestTensorFunc, Evaluate)
{ {
const char* func = "sin(x) | cos (y) | exp(z) | sin(x)*cos(y)"; const char* func = "sin(x) | cos (y) | exp(z) | sin(x)*cos(y)";
TensorFuncExpr f(func); TensorFuncExpr f1(func);
TensorFuncExprAd f2(func);
for (double x : {0.1, 0.2, 0.3}) for (double x : {0.1, 0.2, 0.3})
for (double y : {0.5, 0.6, 0.7}) for (double y : {0.5, 0.6, 0.7})
for (double z : {0.8, 0.9, 1.0}) { for (double z : {0.8, 0.9, 1.0}) {
const Vec3 X(x,y,z); const Vec3 X(x,y,z);
const Tensor fx = f(X);
const Tensor r({sin(x), cos(y), exp(z), sin(x)*cos(y)}); const Tensor r({sin(x), cos(y), exp(z), sin(x)*cos(y)});
for (size_t i = 1; i <= 2; ++i)
for (size_t j = 1; j <= 2; ++j) for (const TensorFunc* fp : {static_cast<const TensorFunc*>(&f1),
EXPECT_DOUBLE_EQ(fx(i,j), r(i,j)); static_cast<const TensorFunc*>(&f2)}) {
const Tensor fx = (*fp)(X);
for (size_t i = 1; i <= 2; ++i)
for (size_t j = 1; j <= 2; ++j)
EXPECT_DOUBLE_EQ(fx(i,j), r(i,j));
}
} }
} }
@ -491,11 +544,13 @@ TEST(TestTensorFunction, Gradient2D)
const char* g_x = "cos(x)*sin(y) | 2*x*y*y | exp(x)*exp(2*y) | -2*exp(-2*x)*exp(y)"; const char* g_x = "cos(x)*sin(y) | 2*x*y*y | exp(x)*exp(2*y) | -2*exp(-2*x)*exp(y)";
const char* g_y = "sin(x)*cos(y) | 2*x*x*y | 2*exp(x)*exp(2*y) | exp(-2*x)*exp(y)"; const char* g_y = "sin(x)*cos(y) | 2*x*x*y | 2*exp(x)*exp(2*y) | exp(-2*x)*exp(y)";
TensorFuncExpr f(g); TensorFuncExpr f1(g);
f.addDerivative(g_x,"",1); f1.addDerivative(g_x,"",1);
f.addDerivative(g_y,"",2); f1.addDerivative(g_y,"",2);
EXPECT_TRUE(f.isConstant()); TensorFuncExprAd f2(g);
EXPECT_TRUE(f1.isConstant());
for (double x : {0.1, 0.2, 0.3}) for (double x : {0.1, 0.2, 0.3})
for (double y : {0.5, 0.6, 0.7}) { for (double y : {0.5, 0.6, 0.7}) {
@ -503,15 +558,17 @@ TEST(TestTensorFunction, Gradient2D)
utl::matrix3d<Real> r(2,2,2); utl::matrix3d<Real> r(2,2,2);
r.fill(std::array{cos(x)*sin(y), 2*x*y*y, exp(x)*exp(2*y), -2*exp(-2*x)*exp(y), r.fill(std::array{cos(x)*sin(y), 2*x*y*y, exp(x)*exp(2*y), -2*exp(-2*x)*exp(y),
sin(x)*cos(y), 2*x*x*y, 2*exp(x)*exp(2*y), exp(-2*x)*exp(y)}.data()); sin(x)*cos(y), 2*x*x*y, 2*exp(x)*exp(2*y), exp(-2*x)*exp(y)}.data());
for (const TensorFunc* fp : {static_cast<const TensorFunc*>(&f1),
const utl::matrix3d<Real> grad = f.gradient(X); static_cast<const TensorFunc*>(&f2)}) {
for (size_t d = 1; d <= 2; ++d) { const utl::matrix3d<Real> grad = fp->gradient(X);
const Tensor dx = f.deriv(X,d); for (size_t d = 1; d <= 2; ++d) {
for (size_t i = 1; i <= 2; ++i) const Tensor dx = fp->deriv(X,d);
for (size_t j = 1; j <= 2; ++j) { for (size_t i = 1; i <= 2; ++i)
EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d)); for (size_t j = 1; j <= 2; ++j) {
EXPECT_DOUBLE_EQ(grad(i,j,d), r(i,j,d)); EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d));
} EXPECT_DOUBLE_EQ(grad(i,j,d), r(i,j,d));
}
}
} }
} }
} }
@ -570,12 +627,15 @@ TEST(TestTensorFunction, Gradient3D)
"exp(-2*x)*exp(y) | x*y | x*y*2*z |" "exp(-2*x)*exp(y) | x*y | x*y*2*z |"
"0.0 | 0.0 | 1.0"; "0.0 | 0.0 | 1.0";
TensorFuncExpr f(g); TensorFuncExpr f1(g);
f.addDerivative(g_x,"",1); f1.addDerivative(g_x,"",1);
f.addDerivative(g_y,"",2); f1.addDerivative(g_y,"",2);
f.addDerivative(g_z,"",3); f1.addDerivative(g_z,"",3);
EXPECT_TRUE(f.isConstant()); TensorFuncExprAd f2(g);
EXPECT_TRUE(f1.isConstant());
EXPECT_TRUE(f2.isConstant());
for (double x : {0.1, 0.2, 0.3}) for (double x : {0.1, 0.2, 0.3})
for (double y : {0.5, 0.6, 0.7}) for (double y : {0.5, 0.6, 0.7})
@ -594,14 +654,17 @@ TEST(TestTensorFunction, Gradient3D)
exp(-2*x)*exp(y), x*y, x*y*2*z, exp(-2*x)*exp(y), x*y, x*y*2*z,
0.0, 0.0, 1.0}.data()); 0.0, 0.0, 1.0}.data());
const utl::matrix3d<Real> grad = f.gradient(X); for (const TensorFunc* fp : {static_cast<const TensorFunc*>(&f1),
for (size_t d = 1; d <= 3; ++d) { static_cast<const TensorFunc*>(&f2)}) {
const Tensor dx = f.deriv(X,d); const utl::matrix3d<Real> grad = fp->gradient(X);
for (size_t i = 1; i <= 3; ++i) for (size_t d = 1; d <= 3; ++d) {
for (size_t j = 1; j <= 3; ++j) { const Tensor dx = fp->deriv(X,d);
EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d)); for (size_t i = 1; i <= 3; ++i)
EXPECT_DOUBLE_EQ(grad(i,j,d), r(i,j,d)); for (size_t j = 1; j <= 3; ++j) {
} EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d));
EXPECT_DOUBLE_EQ(grad(i,j,d), r(i,j,d));
}
}
} }
} }
} }
@ -674,12 +737,15 @@ TEST(TestTensorFunction, Hessian2D)
const char* g_xy = "cos(x)*cos(y) | 2*x*2*y | exp(x)*2*exp(2*y) | -2*exp(-2*x)*exp(y)"; const char* g_xy = "cos(x)*cos(y) | 2*x*2*y | exp(x)*2*exp(2*y) | -2*exp(-2*x)*exp(y)";
const char* g_yy = "-sin(x)*sin(y) | 2*x*x | exp(x)*4*exp(2*y) | exp(-2*x)*exp(y)"; const char* g_yy = "-sin(x)*sin(y) | 2*x*x | exp(x)*4*exp(2*y) | exp(-2*x)*exp(y)";
TensorFuncExpr f(g); TensorFuncExpr f1(g);
f.addDerivative(g_xx,"",1,1); f1.addDerivative(g_xx,"",1,1);
f.addDerivative(g_yy,"",2,2); f1.addDerivative(g_yy,"",2,2);
f.addDerivative(g_xy,"",1,2); f1.addDerivative(g_xy,"",1,2);
EXPECT_TRUE(f.isConstant()); TensorFuncExprAd f2(g);
EXPECT_TRUE(f1.isConstant());
EXPECT_TRUE(f2.isConstant());
for (double x : {0.1, 0.2, 0.3}) for (double x : {0.1, 0.2, 0.3})
for (double y : {0.5, 0.6, 0.7}) { for (double y : {0.5, 0.6, 0.7}) {
@ -698,16 +764,19 @@ TEST(TestTensorFunction, Hessian2D)
-sin(x)*sin(y), 2*x*x, -sin(x)*sin(y), 2*x*x,
exp(x)*4*exp(2*y), exp(-2*x)*exp(y)}.data()); exp(x)*4*exp(2*y), exp(-2*x)*exp(y)}.data());
const utl::matrix4d<Real> hess = f.hessian(X); for (const TensorFunc* fp : {static_cast<const TensorFunc*>(&f1),
for (size_t d1 = 1; d1 <= 2; ++d1) static_cast<const TensorFunc*>(&f2)}) {
for (size_t d2 = 1; d2 <= 2; ++d2) { const utl::matrix4d<Real> hess = fp->hessian(X);
const Tensor dx = f.dderiv(X,d1,d2); for (size_t d1 = 1; d1 <= 2; ++d1)
for (size_t i = 1; i <= 2; ++i) for (size_t d2 = 1; d2 <= 2; ++d2) {
for (size_t j = 1; j <= 2; ++j) { const Tensor dx = fp->dderiv(X,d1,d2);
EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d1,d2)); for (size_t i = 1; i <= 2; ++i)
EXPECT_DOUBLE_EQ(hess(i,j,d1,d2), r(i,j,d1,d2)); for (size_t j = 1; j <= 2; ++j) {
} EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d1,d2));
} EXPECT_DOUBLE_EQ(hess(i,j,d1,d2), r(i,j,d1,d2));
}
}
}
} }
} }
@ -736,15 +805,18 @@ TEST(TestTensorFunction, Hessian3D)
"0.0 | 0.0 | 2.0*x*y |" "0.0 | 0.0 | 2.0*x*y |"
"0.0 | 0.0 | 2.0"; "0.0 | 0.0 | 2.0";
TensorFuncExpr f(g); TensorFuncExpr f1(g);
f.addDerivative(g_xx,"",1,1); f1.addDerivative(g_xx,"",1,1);
f.addDerivative(g_yy,"",2,2); f1.addDerivative(g_yy,"",2,2);
f.addDerivative(g_zz,"",3,3); f1.addDerivative(g_zz,"",3,3);
f.addDerivative(g_xy,"",1,2); f1.addDerivative(g_xy,"",1,2);
f.addDerivative(g_xz,"",1,3); f1.addDerivative(g_xz,"",1,3);
f.addDerivative(g_yz,"",2,3); f1.addDerivative(g_yz,"",2,3);
EXPECT_TRUE(f.isConstant()); TensorFuncExprAd f2(g);
EXPECT_TRUE(f1.isConstant());
EXPECT_TRUE(f2.isConstant());
for (double x : {0.1, 0.2, 0.3}) for (double x : {0.1, 0.2, 0.3})
for (double y : {0.5, 0.6, 0.7}) for (double y : {0.5, 0.6, 0.7})
@ -787,16 +859,19 @@ TEST(TestTensorFunction, Hessian3D)
0.0, 0.0, 2.0*x*y, 0.0, 0.0, 2.0*x*y,
0.0, 0.0, 2.0}.data()); 0.0, 0.0, 2.0}.data());
const utl::matrix4d<Real> hess = f.hessian(X); for (const TensorFunc* fp : {static_cast<const TensorFunc*>(&f1),
for (size_t d1 = 1; d1 <= 3; ++d1) static_cast<const TensorFunc*>(&f2)}) {
for (size_t d2 = 1; d2 <= 3; ++d2) { const utl::matrix4d<Real> hess = fp->hessian(X);
const Tensor dx = f.dderiv(X,d1,d2); for (size_t d1 = 1; d1 <= 3; ++d1)
for (size_t i = 1; i <= 3; ++i) for (size_t d2 = 1; d2 <= 3; ++d2) {
for (size_t j = 1; j <= 3; ++j) { const Tensor dx = fp->dderiv(X,d1,d2);
EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d1,d2)); for (size_t i = 1; i <= 3; ++i)
EXPECT_DOUBLE_EQ(hess(i,j,d1,d2), r(i,j,d1,d2)); for (size_t j = 1; j <= 3; ++j) {
} EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d1,d2));
} EXPECT_DOUBLE_EQ(hess(i,j,d1,d2), r(i,j,d1,d2));
}
}
}
} }
} }
@ -854,17 +929,24 @@ TEST(TestSTensorFunc, Evaluate2D)
{ {
const char* func = "sin(x) | cos (y) | exp(z)"; const char* func = "sin(x) | cos (y) | exp(z)";
STensorFuncExpr f(func); STensorFuncExpr f1(func);
STensorFuncExprAd f2(func);
EXPECT_TRUE(f1.isConstant());
EXPECT_TRUE(f2.isConstant());
for (double x : {0.1, 0.2, 0.3}) for (double x : {0.1, 0.2, 0.3})
for (double y : {0.5, 0.6, 0.7}) for (double y : {0.5, 0.6, 0.7})
for (double z : {0.8, 0.9, 1.0}) { for (double z : {0.8, 0.9, 1.0}) {
const Vec3 X(x,y,z); const Vec3 X(x,y,z);
const SymmTensor fx = f(X);
const Tensor r({sin(x), exp(z), exp(z), cos(y)}); const Tensor r({sin(x), exp(z), exp(z), cos(y)});
for (size_t i = 1; i <= 2; ++i) for (const STensorFunc* fp : {static_cast<const STensorFunc*>(&f1),
for (size_t j = 1; j <= 2; ++j) static_cast<const STensorFunc*>(&f2)}) {
EXPECT_DOUBLE_EQ(fx(i,j), r(i,j)); const SymmTensor fx = (*fp)(X);
for (size_t i = 1; i <= 2; ++i)
for (size_t j = 1; j <= 2; ++j)
EXPECT_DOUBLE_EQ(fx(i,j), r(i,j));
}
} }
} }
@ -894,19 +976,23 @@ TEST(TestSTensorFunc, Evaluate3D)
const char* func = "sin(x) | cos (y) | exp(z) |" const char* func = "sin(x) | cos (y) | exp(z) |"
"sin(x)*sin(y) | sin(x)*cos(y) | exp(x)*exp(y)"; "sin(x)*sin(y) | sin(x)*cos(y) | exp(x)*exp(y)";
STensorFuncExpr f(func); STensorFuncExpr f1(func);
STensorFuncExprAd f2(func);
for (double x : {0.1, 0.2, 0.3}) for (double x : {0.1, 0.2, 0.3})
for (double y : {0.5, 0.6, 0.7}) for (double y : {0.5, 0.6, 0.7})
for (double z : {0.8, 0.9, 1.0}) { for (double z : {0.8, 0.9, 1.0}) {
const Vec3 X(x,y,z); const Vec3 X(x,y,z);
const SymmTensor fx = f(X);
const Tensor r({sin(x), sin(x)*sin(y), exp(x)*exp(y), const Tensor r({sin(x), sin(x)*sin(y), exp(x)*exp(y),
sin(x)*sin(y), cos(y), sin(x)*cos(y), sin(x)*sin(y), cos(y), sin(x)*cos(y),
exp(x)*exp(y), sin(x)*cos(y), exp(z)}); exp(x)*exp(y), sin(x)*cos(y), exp(z)});
for (size_t i = 1; i <= 3; ++i) for (const STensorFunc* fp : {static_cast<const STensorFunc*>(&f1),
for (size_t j = 1; j <= 3; ++j) static_cast<const STensorFunc*>(&f2)}) {
EXPECT_DOUBLE_EQ(fx(i,j), r(i,j)); const SymmTensor fx = (*fp)(X);
for (size_t i = 1; i <= 3; ++i)
for (size_t j = 1; j <= 3; ++j)
EXPECT_DOUBLE_EQ(fx(i,j), r(i,j));
}
} }
} }
@ -917,11 +1003,14 @@ TEST(TestSTensorFunction, Gradient2D)
const char* g_x = "cos(x)*sin(y) | exp(x)*exp(2*y) | 2*x*y*y"; const char* g_x = "cos(x)*sin(y) | exp(x)*exp(2*y) | 2*x*y*y";
const char* g_y = "sin(x)*cos(y) | 2*exp(x)*exp(2*y) | 2*x*x*y"; const char* g_y = "sin(x)*cos(y) | 2*exp(x)*exp(2*y) | 2*x*x*y";
STensorFuncExpr f(g); STensorFuncExpr f1(g);
f.addDerivative(g_x,"",1); f1.addDerivative(g_x,"",1);
f.addDerivative(g_y,"",2); f1.addDerivative(g_y,"",2);
EXPECT_TRUE(f.isConstant()); STensorFuncExprAd f2(g);
EXPECT_TRUE(f1.isConstant());
EXPECT_TRUE(f2.isConstant());
for (double x : {0.1, 0.2, 0.3}) for (double x : {0.1, 0.2, 0.3})
for (double y : {0.5, 0.6, 0.7}) { for (double y : {0.5, 0.6, 0.7}) {
@ -930,14 +1019,17 @@ TEST(TestSTensorFunction, Gradient2D)
r.fill(std::array{cos(x)*sin(y), 2*x*y*y, 2*x*y*y, exp(x)*exp(2*y), r.fill(std::array{cos(x)*sin(y), 2*x*y*y, 2*x*y*y, exp(x)*exp(2*y),
sin(x)*cos(y), 2*x*x*y, 2*x*x*y, 2*exp(x)*exp(2*y)}.data()); sin(x)*cos(y), 2*x*x*y, 2*x*x*y, 2*exp(x)*exp(2*y)}.data());
const utl::matrix3d<Real> grad = f.gradient(X); for (const STensorFunc* fp : {static_cast<const STensorFunc*>(&f1),
for (size_t d = 1; d <= 2; ++d) { static_cast<const STensorFunc*>(&f2)}) {
const SymmTensor dx = f.deriv(X,d); const utl::matrix3d<Real> grad = fp->gradient(X);
for (size_t i = 1; i <= 2; ++i) for (size_t d = 1; d <= 2; ++d) {
for (size_t j = 1; j <= 2; ++j) { const SymmTensor dx = fp->deriv(X,d);
EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d)); for (size_t i = 1; i <= 2; ++i)
EXPECT_DOUBLE_EQ(grad(i,j,d), r(i,j,d)); for (size_t j = 1; j <= 2; ++j) {
} EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d));
EXPECT_DOUBLE_EQ(grad(i,j,d), r(i,j,d));
}
}
} }
} }
} }
@ -1022,12 +1114,15 @@ TEST(TestSTensorFunction, Gradient3D)
const char* g_z = "sin(x)*sin(y)*cos(z) | exp(x)*exp(2*y)*exp(z) | x*x*y*y*2*z |" const char* g_z = "sin(x)*sin(y)*cos(z) | exp(x)*exp(2*y)*exp(z) | x*x*y*y*2*z |"
"x*y | x*x*y | 1"; "x*y | x*x*y | 1";
STensorFuncExpr f(g); STensorFuncExpr f1(g);
f.addDerivative(g_x,"",1); f1.addDerivative(g_x,"",1);
f.addDerivative(g_y,"",2); f1.addDerivative(g_y,"",2);
f.addDerivative(g_z,"",3); f1.addDerivative(g_z,"",3);
EXPECT_TRUE(f.isConstant()); STensorFuncExprAd f2(g);
EXPECT_TRUE(f1.isConstant());
EXPECT_TRUE(f2.isConstant());
for (double x : {0.1, 0.2, 0.3}) for (double x : {0.1, 0.2, 0.3})
for (double y : {0.5, 0.6, 0.7}) for (double y : {0.5, 0.6, 0.7})
@ -1046,14 +1141,17 @@ TEST(TestSTensorFunction, Gradient3D)
x*y, exp(x)*exp(2*y)*exp(z), x*x*y, x*y, exp(x)*exp(2*y)*exp(z), x*x*y,
1.0, x*x*y, x*x*y*y*2*z}.data()); 1.0, x*x*y, x*x*y*y*2*z}.data());
const utl::matrix3d<Real> grad = f.gradient(X); for (const STensorFunc* fp : {static_cast<const STensorFunc*>(&f1),
for (size_t d = 1; d <= 3; ++d) { static_cast<const STensorFunc*>(&f2)}) {
const SymmTensor dx = f.deriv(X,d); const utl::matrix3d<Real> grad = fp->gradient(X);
for (size_t i = 1; i <= 3; ++i) for (size_t d = 1; d <= 3; ++d) {
for (size_t j = 1; j <= 3; ++j) { const SymmTensor dx = fp->deriv(X,d);
EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d)); for (size_t i = 1; i <= 3; ++i)
EXPECT_DOUBLE_EQ(grad(i,j,d), r(i,j,d)); for (size_t j = 1; j <= 3; ++j) {
} EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d));
EXPECT_DOUBLE_EQ(grad(i,j,d), r(i,j,d));
}
}
} }
} }
} }
@ -1145,12 +1243,15 @@ TEST(TestSTensorFunction, Hessian2D)
const char* g_xy = "cos(x)*cos(y) | exp(x)*2*exp(2*y) | 2*x*2*y"; const char* g_xy = "cos(x)*cos(y) | exp(x)*2*exp(2*y) | 2*x*2*y";
const char* g_yy = "-sin(x)*sin(y) | exp(x)*4*exp(2*y) | x*x*2"; const char* g_yy = "-sin(x)*sin(y) | exp(x)*4*exp(2*y) | x*x*2";
STensorFuncExpr f(g); STensorFuncExpr f1(g);
f.addDerivative(g_xx,"",1,1); f1.addDerivative(g_xx,"",1,1);
f.addDerivative(g_xy,"",1,2); f1.addDerivative(g_xy,"",1,2);
f.addDerivative(g_yy,"",2,2); f1.addDerivative(g_yy,"",2,2);
EXPECT_TRUE(f.isConstant()); STensorFuncExprAd f2(g);
EXPECT_TRUE(f1.isConstant());
EXPECT_TRUE(f2.isConstant());
for (double x : {0.1, 0.2, 0.3}) for (double x : {0.1, 0.2, 0.3})
for (double y : {0.5, 0.6, 0.7}) { for (double y : {0.5, 0.6, 0.7}) {
@ -1168,16 +1269,19 @@ TEST(TestSTensorFunction, Hessian2D)
-sin(x)*sin(y), x*x*2, -sin(x)*sin(y), x*x*2,
x*x*2, exp(x)*4*exp(2*y)}.data()); x*x*2, exp(x)*4*exp(2*y)}.data());
const utl::matrix4d<Real> hess = f.hessian(X); for (const STensorFunc* fp : {static_cast<const STensorFunc*>(&f1),
for (size_t d1 = 1; d1 <= 2; ++d1) static_cast<const STensorFunc*>(&f2)}) {
for (size_t d2 = 1; d2 <= 2; ++d2) { const utl::matrix4d<Real> hess = fp->hessian(X);
const SymmTensor dx = f.dderiv(X,d1,d2); for (size_t d1 = 1; d1 <= 2; ++d1)
for (size_t i = 1; i <= 2; ++i) for (size_t d2 = 1; d2 <= 2; ++d2) {
for (size_t j = 1; j <= 2; ++j) { const SymmTensor dx = fp->dderiv(X,d1,d2);
EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d1,d2)); for (size_t i = 1; i <= 2; ++i)
EXPECT_DOUBLE_EQ(hess(i,j,d1,d2), r(i,j,d1,d2)); for (size_t j = 1; j <= 2; ++j) {
} EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d1,d2));
} EXPECT_DOUBLE_EQ(hess(i,j,d1,d2), r(i,j,d1,d2));
}
}
}
} }
} }
@ -1200,15 +1304,18 @@ TEST(TestSTensorFunction, Hessian3D)
"0.0 | 0.0 | 0.0"; "0.0 | 0.0 | 0.0";
STensorFuncExpr f(g); STensorFuncExpr f1(g);
f.addDerivative(g_xx,"",1,1); f1.addDerivative(g_xx,"",1,1);
f.addDerivative(g_xy,"",1,2); f1.addDerivative(g_xy,"",1,2);
f.addDerivative(g_xz,"",1,3); f1.addDerivative(g_xz,"",1,3);
f.addDerivative(g_yy,"",2,2); f1.addDerivative(g_yy,"",2,2);
f.addDerivative(g_yz,"",2,3); f1.addDerivative(g_yz,"",2,3);
f.addDerivative(g_zz,"",3,3); f1.addDerivative(g_zz,"",3,3);
EXPECT_TRUE(f.isConstant()); STensorFuncExprAd f2(g);
EXPECT_TRUE(f1.isConstant());
EXPECT_TRUE(f2.isConstant());
for (double x : {0.1, 0.2, 0.3}) for (double x : {0.1, 0.2, 0.3})
for (double y : {0.5, 0.6, 0.7}) for (double y : {0.5, 0.6, 0.7})
@ -1251,17 +1358,20 @@ TEST(TestSTensorFunction, Hessian3D)
0.0, exp(x)*exp(2*y)*exp(z), 0.0, 0.0, exp(x)*exp(2*y)*exp(z), 0.0,
0.0, 0.0, x*x*y*y*2}.data()); 0.0, 0.0, x*x*y*y*2}.data());
const utl::matrix4d<Real> hess = f.hessian(X); for (const STensorFunc* fp : {static_cast<const STensorFunc*>(&f1),
for (size_t d1 = 1; d1 <= 3; ++d1) static_cast<const STensorFunc*>(&f2)}) {
for (size_t d2 = 1; d2 <= 3; ++d2) { const utl::matrix4d<Real> hess = fp->hessian(X);
const SymmTensor dx = f.dderiv(X,d1,d2); for (size_t d1 = 1; d1 <= 3; ++d1)
for (size_t i = 1; i <= 3; ++i) for (size_t d2 = 1; d2 <= 3; ++d2) {
for (size_t j = 1; j <= 3; ++j) { const SymmTensor dx = fp->dderiv(X,d1,d2);
EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d1,d2)); for (size_t i = 1; i <= 3; ++i)
EXPECT_DOUBLE_EQ(hess(i,j,d1,d2), r(i,j,d1,d2)); for (size_t j = 1; j <= 3; ++j) {
} EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d1,d2));
} EXPECT_DOUBLE_EQ(hess(i,j,d1,d2), r(i,j,d1,d2));
} }
}
}
}
} }
@ -1384,5 +1494,5 @@ TEST(TestEvalFunction, Derivatives)
EXPECT_DOUBLE_EQ(f.dderiv(X,3,1), cos(x)*sin(y)*cos(z)*sin(t)); EXPECT_DOUBLE_EQ(f.dderiv(X,3,1), cos(x)*sin(y)*cos(z)*sin(t));
EXPECT_DOUBLE_EQ(f.dderiv(X,3,2), sin(x)*cos(y)*cos(z)*sin(t)); EXPECT_DOUBLE_EQ(f.dderiv(X,3,2), sin(x)*cos(y)*cos(z)*sin(t));
EXPECT_DOUBLE_EQ(f.dderiv(X,3,3), -sin(x)*sin(y)*sin(z)*sin(t)); EXPECT_DOUBLE_EQ(f.dderiv(X,3,3), -sin(x)*sin(y)*sin(z)*sin(t));
} }
} }