Commit c6ce8c93 authored by Thibault Lejemble's avatar Thibault Lejemble

AutoDiff for weighting functions test

parent 87112240
......@@ -12,9 +12,84 @@
#include "../common/testing.h"
#include "../common/testUtils.h"
#include <unsupported/Eigen/AutoDiff>
using namespace std;
using namespace Grenaille;
template<typename DataPoint, typename WeightKernel>
void testFunctionAutoDiff()
{
// Define related structure
typedef typename WeightKernel::Scalar ScalarDiff;
typedef typename ScalarDiff::Scalar Scalar;
typedef typename DataPoint::VectorType VectorTypeDiff;
typedef typename DataPoint::MatrixType MatrixTypeDiff;
Scalar epsilon = testEpsilon<Scalar>();
ScalarDiff t = Eigen::internal::random<Scalar>(0.10, 10.0);
ScalarDiff tDiff = ScalarDiff(t.value(), DataPoint::Dim+1, 0);
DistWeightFunc<DataPoint, WeightKernel> wfunc(t);
DistWeightFunc<DataPoint, WeightKernel> wfuncDiff(tDiff);
// rejection sampling inside a sphere of radius t
VectorTypeDiff x = t*VectorTypeDiff::Random();
while ( x.norm() < 0.001*t || 0.999*t < x.norm())
x = t*VectorTypeDiff::Random();
VectorTypeDiff xDiff;
for(int i=0; i<DataPoint::Dim; ++i) {
xDiff[i] = ScalarDiff(x[i].value(), DataPoint::Dim+1, i+1);
}
DataPoint dummy;
ScalarDiff w;
VectorTypeDiff dw;
// reference values
VectorTypeDiff dx_w = wfunc.spacedw(x, dummy);
MatrixTypeDiff d2x_w = wfunc.spaced2w(x, dummy);
ScalarDiff dt_w = wfunc.scaledw(x, dummy);
ScalarDiff d2t_w = wfunc.scaled2w(x, dummy);
VectorTypeDiff d2tx_w = wfunc.scaleSpaced2w(x, dummy);
// auto diff values
VectorTypeDiff dx_w_;
MatrixTypeDiff d2x_w_;
ScalarDiff dt_w_;
ScalarDiff d2t_w_;
VectorTypeDiff d2tx_w_;
// 1st order space derivative
w = wfunc.w(xDiff, dummy);
dx_w_ = w.derivatives().template tail<DataPoint::Dim>();
// 2nd order space derivative
dw = wfunc.spacedw(xDiff, dummy);
for(int i=0; i<DataPoint::Dim; ++i) {
d2x_w_.row(i) = dw[i].derivatives().template tail<DataPoint::Dim>();
}
// 1st order scale derivative
w = wfuncDiff.w(x, dummy);
dt_w_ = w.derivatives()[0];
// 2nd order scale derivative
w = wfuncDiff.scaledw(x, dummy);
d2t_w_ = w.derivatives()[0];
// cross derivative
w = wfuncDiff.scaledw(xDiff, dummy);
d2tx_w_ = w.derivatives().template tail<DataPoint::Dim>();
VERIFY( (dx_w-dx_w_).array().abs().maxCoeff() < epsilon );
VERIFY( (d2x_w-d2x_w_).array().abs().maxCoeff() < epsilon );
VERIFY( (d2tx_w-d2tx_w_).array().abs().maxCoeff() < epsilon );
VERIFY( std::abs((dt_w-dt_w_).value()) < epsilon );
VERIFY( std::abs((d2t_w-d2t_w_).value()) < epsilon );
}
template<typename DataPoint, typename WeightKernel>
void testFunction()
{
......@@ -122,10 +197,20 @@ void callSubTests()
typedef SmoothWeightKernel<Scalar> SmoothKernel;
typedef ConstantWeightKernel<Scalar> ConstantKernel;
typedef Eigen::AutoDiffScalar<Eigen::Matrix<Scalar,Dim+1,1>> ScalarDiff;
typedef PointPositionNormal<ScalarDiff, Dim> DataPointDiff;
typedef SmoothWeightKernel<ScalarDiff> SmoothKernelDiff;
typedef ConstantWeightKernel<ScalarDiff> ConstantKernelDiff;
for(int i = 0; i < g_repeat; ++i)
{
CALL_SUBTEST(( testFunction<DataPoint, SmoothKernel>() ));
CALL_SUBTEST(( testFunction<DataPoint, ConstantKernel>() ));
CALL_SUBTEST(( testFunctionAutoDiff<DataPointDiff, SmoothKernelDiff>() ));
CALL_SUBTEST(( testFunctionAutoDiff<DataPointDiff, ConstantKernelDiff>() ));
}
cout << "ok" << endl;
}
......
......@@ -12,9 +12,40 @@
#include "../common/testing.h"
#include "../common/testUtils.h"
#include <unsupported/Eigen/AutoDiff>
using namespace std;
using namespace Grenaille;
template<class Kernel>
void testFunctionAutoDiff()
{
typedef typename Kernel::Scalar ScalarDiff;
typedef typename ScalarDiff::Scalar Scalar;
Scalar step = Scalar(0.05);
int n = Scalar(1.)/step;
Kernel k;
Scalar epsilon = testEpsilon<Scalar>();
// compare to automatic differentiation
for(int i=1; i<=n; ++i)
{
ScalarDiff a(i*step, 1, 0);
ScalarDiff f = k.f(a);
ScalarDiff df = k.df(a);
ScalarDiff ddf = k.ddf(a);
Scalar diff1 = std::abs( f.derivatives()[0] - df.value());
Scalar diff2 = std::abs(df.derivatives()[0] - ddf.value());
VERIFY(diff1 < epsilon);
VERIFY(diff2 < epsilon);
}
}
template<class Kernel>
void testFunction()
{
......@@ -53,8 +84,14 @@ void testFunction()
template<typename Scalar>
void callSubTests()
{
typedef Eigen::AutoDiffScalar<Eigen::Matrix<Scalar,1,1>> ScalarDiff;
typedef SmoothWeightKernel<Scalar> Kernel;
typedef SmoothWeightKernel<ScalarDiff> KernelAutoDiff;
CALL_SUBTEST(( testFunction<Kernel>() ));
CALL_SUBTEST(( testFunctionAutoDiff<KernelAutoDiff>() ));
cout << "ok" << endl;
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment