Question

The Solution from the EDIT is now posted as answer.

Old Question

I want to include spline interpolation of Eigen::Spline into a larger formula and want to determine the derivative of this formula with the help of Eigen::AutoDiff.

I tried the following code:

#include <iostream>
#include <Eigen/Eigen>
#include <unsupported/Eigen/AutoDiff>
#include <unsupported/Eigen/Splines>

double array[] = {
    1.0, 2.0,
    3.0, 4.0,
    5.0, 1.0,
    6.0, 2.0
};

constexpr size_t nCols=2;
constexpr size_t rowSize=sizeof(array)/sizeof(double)/nCols;
constexpr size_t derOrder=1;

typedef Eigen::Map<Eigen::Matrix<double,1,rowSize>,Eigen::Unaligned,Eigen::InnerStride<nCols> > Map;

typedef const Eigen::Spline<double,nCols-1> Spline;

constexpr size_t interpolOrder=3;
constexpr double& xStop=array[nCols*(rowSize-1)];

Spline spline=Eigen::SplineFitting<Spline>::Interpolate(Map(array+1),interpolOrder,Map(static_cast<double*>(array))/xStop);

typedef Eigen::AutoDiffScalar<Eigen::Matrix<double,1,2> > DerType;

DerType f(const DerType& x){
    DerType ret;
    auto y = spline.derivatives<1>(x.value()/xStop);
//*** Compilation is okay if previous line is substituted with:
// Eigen::Array<double,1,2> y; y << 1.0, 2.0;
    ret.value() = y(0,0);
    ret.derivatives() = x.derivatives()*y(0,0);
    return ret;
}


int main() {
    DerType x(1.0,DerType::DerType(1.0,1.0));

    auto y = f(x);

    std::cout << "\nValue=" << y.value() << "\nDer=" << y.derivatives() << '\n';

    return 0;
}


/*
    Local Variables:
    compile-command: "g++ -g -std=c++11 -I/usr/local/include/eigen3 eigenInterpolAD.cc -o a.exe && (echo \"Running\"; ./a.exe);"
    End:
*/

Pityingly, the compilation of the code gives the following error message:

g++ -g -std=c++11 -I/usr/local/include/eigen3 eigenInterpolAD.cc -o a.exe && (echo "Running"; ./a.exe);
In file included from /usr/local/include/eigen3/Eigen/Core:254:0,
                 from /usr/local/include/eigen3/Eigen/Dense:1,
                 from /usr/local/include/eigen3/Eigen/Eigen:1,
                 from eigenInterpolAD.cc:5:
/usr/local/include/eigen3/Eigen/src/Core/PlainObjectBase.h: In instantiation of ‘static void Eigen::PlainObjectBase<Derived>::_check_template_params() [with Derived = Eigen::Array<double, 1, -1, 0, 1, 2>]’:
/usr/local/include/eigen3/Eigen/src/Core/Array.h:195:36:   required from ‘Eigen::Array<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::Array(const Eigen::Array<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>&) [with _Scalar = double; int _Rows = 1; int _Cols = -1; int _Options = 0; int _MaxRows = 1; int _MaxCols = 2]’
eigenInterpolAD.cc:33:48:   required from here
/usr/local/include/eigen3/Eigen/src/Core/util/StaticAssert.h:32:40: error: static assertion failed: INVALID_MATRIX_TEMPLATE_PARAMETERS
     #define EIGEN_STATIC_ASSERT(X,MSG) static_assert(X,#MSG);
                                        ^
/usr/local/include/eigen3/Eigen/src/Core/PlainObjectBase.h:657:7: note: in expansion of macro ‘EIGEN_STATIC_ASSERT’
       EIGEN_STATIC_ASSERT((EIGEN_IMPLIES(MaxRowsAtCompileTime==1 && MaxColsAtCompileTime!=1, (Options&RowMajor)==RowMajor)
       ^

How can I avoid the error? Thanks in advance for any helpful hint.

The eigen-version is 3.2.1. Compilation also fails with versions 3.1 and 3.0 of Eigen.

The compiler version is:

gcc version 4.8.1 (Ubuntu/Linaro 4.8.1-10ubuntu9)

Was it helpful?

Solution

As user Anycorn pointed out the source code in Section "Old Question" compiles and runs with the very newest version of Eigen on Bitbucket:

$hg showconfig
bundle.mainreporoot=/usr/local/eigen
paths.default=https://bitbucket.org/eigen/eigen
/usr/local/eigen
$hg identify
fabd880592ac tip
/usr/local/eigen

I've tested the principle with the following code:

/**
     Spline interpolation with AD.
 */
#include <iostream>
#include <fstream>
#include <Eigen/Eigen>
#include <unsupported/Eigen/AutoDiff>
#include <unsupported/Eigen/Splines>

constexpr const double array[] = {
    1.0, 2.0,
    3.0, 4.0,
    5.0, 1.0,
    6.0, 2.0
};

constexpr size_t nCols=2;
constexpr size_t rowSize=sizeof(array)/sizeof(double)/nCols;
constexpr size_t derOrder=1;

typedef Eigen::Map<Eigen::Matrix<double,1,rowSize>,Eigen::Unaligned,Eigen::InnerStride<nCols> > Map;

typedef Eigen::Spline<double,nCols-1> Spline;

constexpr size_t interpolOrder=3;
constexpr double xStart=array[0];
constexpr double xStop=array[nCols*(rowSize-1)];
constexpr double xDelta=xStop-xStart;

const Spline spline=Eigen::SplineFitting<Spline>::Interpolate(Map(
        const_cast<double*>(array)+1),
    interpolOrder,
    (Map(const_cast<double*>(array)).array()-xStart)/xDelta);

typedef Eigen::AutoDiffScalar<Eigen::Matrix<double,1,2> > DerType;

DerType f(const DerType& x){
    DerType ret;
    auto y = spline.derivatives<1>((x.value()-xStart)/xDelta);
    ret.value() = y(0,0);
    ret.derivatives() = x.derivatives()*y(0,1)/xDelta;
    return ret;
}


int main() {
    std::ofstream of("/temp/test.dat");

    constexpr size_t n=101;
    const double xStart=0.0;
    const double xStop=6.0;
    const double dx=(xStop-xStart)/(n-1);

    double x=xStart;

    for(size_t i=0; i!=n; i++, x+=dx) {
        DerType xAD(x,DerType::DerType(1.0,1.0));
        auto yAD = f(xAD);

        of << x << ' ' << yAD.value() << ' ' << yAD.derivatives()[0] << '\n';
    }

    of.close();

    return 0;
}


/*
    Local Variables:
    compile-command: "g++ -g -std=c++11 -I/usr/local/include/eigen3 eigenInterpolAD.cc -o a.exe && (echo \"Running\"; ./a.exe);"
    End:
*/

The following plot of the generated data shows that the code works fine now with the hg-version of Eigen. The red curve is the interpolating spline and the green curve is its derivative.

enter image description here

Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top