Unit Test TestSolvePolynomialQuadraticSmoothed disabled

Summary

The (currently disabled) unit test PPSolveFactoryTestFixture.TestSolvePolynomialQuadraticSmoothed is crashing.

Steps to reproduce

Enable the TestSolvePolynomialQuadraticSmoothed test in PPSolveFactoryTest.cpp and run the test.

What is the current bug behavior?

Crash

What is the expected correct behavior?

Successful test

Relevant logs and/or screenshots

[ RUN      ] PPSolveFactoryTestFixture.TestSolvePolynomialQuadraticSmoothed
gsl: ../gsl/gsl_vector_double.h:206: ERROR: index out of range
Default GSL error handler invoked.

Program received signal SIGABRT, Aborted.
0x0000003da42324f5 in raise () from /lib64/libc.so.6
Missing separate debuginfos, use: debuginfo-install glibc-2.12-1.212.el6.x86_64 infinipath-psm-3.0.1-115.1015_open.2.el6.x86_64 libnl-1.1.4-2.el6.x86_64 zlib-1.2.3-29.el6.x86_64
(gdb) bt
#0  0x0000003da42324f5 in raise () from /lib64/libc.so.6
#1  0x0000003da4233cd5 in abort () from /lib64/libc.so.6
#2  0x00007ffff7bce65d in gsl_error () at /var/tmp/gsell/gsl-2.5/src/err/error.c:47
#3  0x00007ffff7d34669 in gsl_vector_ptr () from /opt/psi/Compiler/gsl/2.5/gcc/7.3.0/lib64/libgsl.so.23
#4  0x0000000000b1a010 in interpolation::PPSolveFactory::getDerivs(interpolation::Mesh::Iterator) ()
    at /home/scratch/OPAL/OPAL-fork/OPAL-src/src/Classic/Fields/Interpolation/PPSolveFactory.cpp:289
#5  0x0000000000b22a6c in interpolation::PPSolveFactory::solve() ()
    at /home/scratch/OPAL/OPAL-fork/OPAL-src/src/Classic/Fields/Interpolation/PPSolveFactory.cpp:322
#6  0x00000000006be2f0 in PPSolveFactoryTestFixture_TestSolvePolynomialQuadraticSmoothed_Test::TestBody() ()
    at /home/scratch/OPAL/OPAL-fork/OPAL-src/tests/classic_src/Fields/Interpolation/PPSolveFactoryTest.cpp:194

Possible fixes

This is due to these lines (https://gitlab.psi.ch/OPAL/src/blob/master/src/Classic/Fields/Interpolation/PPSolveFactory.cpp#L284):

            MMatrix<double> coeffs =
                  polynomials_m[nearest.toInteger()]->GetCoefficientsAsMatrix();
            MVector<double> values = coeffs*derivPolyVec_m[i];
            derivValues_m[i] = std::vector<double>(polyDim_m);
            for(int j = 0; j < posDim; ++j) {
                derivValues_m[i][j] = values(j+1);
            }

posDim is 3, while values and derivValues_m[i] have only length 2. Presumably the loop should be until polyDim_m? However while this change runs through without crash, the unit test fails.

Todo List

• Add a method to SquarePolynomialVector that analytically calculates the derivative of the polynomial vector at some point;
  • add appropriate test.
• Add a test to SolveFactory;
  • check that SolveFactory does indeed correctly solve when given derivatives as well as values and returns a correct polynomial.
• Fix the test in PPSolveFactory; in principle this is just handing off values to SolveFactory (it is a wrapper to handle building a number of polynomials from a grid, e.g. a field map).
• Modify the boundary condition handling to improve quality of the fit when doing smoothing.
• Add notes on interpolation to manual or wiki