From de07000cae8faa2d864ba53f44dd05d7475274bd Mon Sep 17 00:00:00 2001 From: Jesse Mapel Date: Mon, 3 Dec 2018 17:26:03 -0700 Subject: [PATCH] Added position interpolation (#41) * added environment file * Added GSL * Fixed travis.yml: * Maybed actually fixed travis * More travis fixes * interpolation now compiling * Added tests * added interp tests * now building * Added better cubci spline test * Added exceptions * Check for openblas error on linux * Force openblas update? * Trying to get a better openblas version * Checking conda version * Specify openblas version? * Added tests to travois * Moved interpolation into separate methods * Added better interpolation tests --- .travis.yml | 1 + CMakeLists.txt | 11 ++-- environment.yml | 2 + include/eal.h | 57 +++++++++-------- src/eal.cpp | 124 ++++++++++++++++++++++++++++++------- tests/CMakeLists.txt | 10 ++- tests/EalTest.cpp | 144 +++++++++++++++++++++++++++++++++++++++++++ tests/Test.cpp | 8 --- 8 files changed, 291 insertions(+), 66 deletions(-) create mode 100644 tests/EalTest.cpp delete mode 100644 tests/Test.cpp diff --git a/.travis.yml b/.travis.yml index 5e262d1..894f33f 100644 --- a/.travis.yml +++ b/.travis.yml @@ -47,6 +47,7 @@ script: - cd build - cmake -DCOVERAGE=ON .. - cmake --build . + - ctest -VV after_success: - | diff --git a/CMakeLists.txt b/CMakeLists.txt index 49fc833..507499d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -22,20 +22,17 @@ find_package(GSL) # Library setup add_library(eal SHARED - src/eal.cpp -) + ${CMAKE_CURRENT_SOURCE_DIR}/src/eal.cpp) set_target_properties(eal PROPERTIES VERSION ${PROJECT_VERSION} - SOVERSION 1 -) + SOVERSION 1) set(EAL_INCLUDE_DIRS "${CMAKE_CURRENT_SOURCE_DIR}/include/" - "${CMAKE_CURRENT_SOURCE_DIR}/include/json") + "${CMAKE_CURRENT_SOURCE_DIR}/include/json") target_include_directories(eal PRIVATE ${GSL_INCLUDE_DIRS} PUBLIC - ${EAL_INCLUDE_DIRS} -) + ${EAL_INCLUDE_DIRS}) # Setup for GoogleTest find_package (Threads) diff --git a/environment.yml b/environment.yml index 6452d76..430f7de 100644 --- a/environment.yml +++ b/environment.yml @@ -2,7 +2,9 @@ name: eal channels: - conda-forge + - default dependencies: - cmake>=3.10 - gsl + - openblas>=0.3.0 diff --git a/include/eal.h b/include/eal.h index e5cdc4d..c33665a 100644 --- a/include/eal.h +++ b/include/eal.h @@ -5,32 +5,41 @@ #include #include -using json = nlohmann::json; - -using namespace std; - namespace eal { - json constructStateFromIsd(const string positionRotationData); - - vector getPosition(vector> coords, vector times, - string interp, double time); - vector getVelocity(vector> coords, vector times, - string interp, double time, bool interpolation); - - vector getPosition(vector coeffs, string interp, double time); - vector getVelocity(vector coeffs, string interp, double time); - - vector getRotation(string from, string to, vector> rotations, - vector times, string interp, double time); - vector getAngularVelocity(string from, string to, vector> rotations, - vector times, string interp, double time, bool interpolation); - - vector getRotation(string from, string to, vector coefficients, - string interp, double time); - vector getAngularVelocity(string from, string to, vector coefficients, - string interp, double time); - + enum interpolation { + linear, + spline + }; + + nlohmann::json constructStateFromIsd(const std::string positionRotationData); + + std::vector getPosition(std::vector> coords, + std::vector times, + interpolation interp, double time); + std::vector getVelocity(std::vector> coords, + std::vector times, + interpolation interp, double time); + + std::vector getPosition(std::vector coeffs, double time); + std::vector getVelocity(std::vector coeffs, double time); + + std::vector getRotation(std::string from, std::string to, + std::vector> rotations, + std::vector times, + interpolation interp, double time); + std::vector getAngularVelocity(std::string from, std::string to, + std::vector> rotations, + std::vector times, + interpolation interp, double time); + + std::vector getRotation(std::string from, std::string to, + std::vector coefficients, double time); + std::vector getAngularVelocity(std::string from, std::string to, + std::vector coefficients, double time); + + double linearInterpolate(std::vector points, std::vector times, double time); + double splineInterpolate(std::vector points, std::vector times, double time); } #endif // EAL_H diff --git a/src/eal.cpp b/src/eal.cpp index 32dafd1..8da683f 100644 --- a/src/eal.cpp +++ b/src/eal.cpp @@ -2,51 +2,78 @@ #include +#include #include #include +#include using json = nlohmann::json; +using namespace std; namespace eal { // Parsing the JSON - json constructStateFromIsd(const std::string positionRotationData) { - // Parse the position and rotation data from isd - json isd = json::parse(positionRotationData); - json state; + json constructStateFromIsd(const string positionRotationData) { + // Parse the position and rotation data from isd + json isd = json::parse(positionRotationData); + json state; - state["m_w"] = isd.at("w"); - state["m_x"] = isd.at("x"); - state["m_y"] = isd.at("y"); - state["m_z"] = isd.at("z"); + state["m_w"] = isd.at("w"); + state["m_x"] = isd.at("x"); + state["m_y"] = isd.at("y"); + state["m_z"] = isd.at("z"); - return state; - } + return state; + } // Positional Functions // Position Data Functions vector getPosition(vector> coords, vector times, - string interp, double time) { + interpolation interp, double time) { + // Check that all of the data sizes are okay + // TODO is there a cleaner way to do this? We're going to have to do this a lot. + if (coords.size() != 3) { + throw invalid_argument("Invalid input positions, expected three vectors."); + } + + // GSL setup vector coordinate = {0.0, 0.0, 0.0}; + switch(interp) { + case linear: + coordinate = { linearInterpolate(coords[0], times, time), + linearInterpolate(coords[1], times, time), + linearInterpolate(coords[2], times, time) }; + break; + + case spline: + coordinate = { splineInterpolate(coords[0], times, time), + splineInterpolate(coords[1], times, time), + splineInterpolate(coords[2], times, time) }; + break; + + default: + throw invalid_argument("Invalid interpolation method."); + break; + } + return coordinate; } vector getVelocity(vector> coords, vector times, - string interp, double time, bool interpolation) { + interpolation interp, double time) { vector coordinate = {0.0, 0.0, 0.0}; return coordinate; } // Postion Function Functions - vector getPosition(vector coeffs, string interp, double time) { + vector getPosition(vector coeffs, double time) { vector coordinate = {0.0, 0.0, 0.0}; return coordinate; } - vector getVelocity(vector coeffs, string interp, double time, - bool interpolation) { + vector getVelocity(vector coeffs, double time) { vector coordinate = {0.0, 0.0, 0.0}; return coordinate; } @@ -55,28 +82,83 @@ namespace eal { // Rotation Data Functions vector getRotation(string from, string to, vector> rotations, - vector times, string interp, double time) { + vector times, interpolation interp, double time) { vector coordinate = {0.0, 0.0, 0.0}; return coordinate; } vector getAngularVelocity(string from, string to, vector> rotations, - vector times, string interp, double time, bool interpolation) { + vector times, interpolation interp, double time) { vector coordinate = {0.0, 0.0, 0.0}; return coordinate; } // Rotation Function Functions - vector getRotation(string from, string to, vector coefficients, - string interp, double time) { + vector getRotation(string from, string to, + vector coefficients, double time) { vector coordinate = {0.0, 0.0, 0.0}; return coordinate; } - vector getAngularVelocity(string from, string to, vector coefficients, - string interp, double time) { + vector getAngularVelocity(string from, string to, + vector coefficients, double time) { vector coordinate = {0.0, 0.0, 0.0}; return coordinate; } + // Interpolation helper functions + double linearInterpolate(vector points, vector times, double time) { + size_t numPoints = points.size(); + if (numPoints < 2) { + throw invalid_argument("At least two points must be input to interpolate over."); + } + if (points.size() != times.size()) { + throw invalid_argument("Invalid interpolation data, must have the same number of points as times."); + } + if (time < times.front() || time > times.back()) { + throw invalid_argument("Invalid interpolation time, outside of input times."); + } + + // GSL setup + gsl_interp *interpolator = gsl_interp_alloc(gsl_interp_linear, numPoints); + gsl_interp_init(interpolator, ×[0], &points[0], numPoints); + gsl_interp_accel *acc = gsl_interp_accel_alloc(); + + // GSL evaluate + double result = gsl_interp_eval(interpolator, ×[0], &points[0], time, acc); + + // GSL clean up + gsl_interp_free(interpolator); + gsl_interp_accel_free(acc); + + return result; + } + + double splineInterpolate(vector points, vector times, double time) { + size_t numPoints = points.size(); + if (numPoints < 2) { + throw invalid_argument("Invalid interpolation data, at least two points are required to interpolate."); + } + if (points.size() != times.size()) { + throw invalid_argument("Invalid interpolation data, must have the same number of points as times."); + } + if (time < times.front() || time > times.back()) { + throw invalid_argument("Invalid interpolation time, outside of input times."); + } + + // GSL setup + gsl_interp *interpolator = gsl_interp_alloc(gsl_interp_cspline, numPoints); + gsl_interp_init(interpolator, ×[0], &points[0], numPoints); + gsl_interp_accel *acc = gsl_interp_accel_alloc(); + + // GSL evaluate + double result = gsl_interp_eval(interpolator, ×[0], &points[0], time, acc); + + // GSL clean up + gsl_interp_free(interpolator); + gsl_interp_accel_free(acc); + + return result; + } + } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index b74a3d3..4755fb7 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -1,12 +1,10 @@ cmake_minimum_required(VERSION 3.10) +# collect all of the test sources file(GLOB test_source "${CMAKE_SOURCE_DIR}/tests/*.cpp") -# Link runEALTests with what we want to test and the GTest and pthread library -add_executable(runEALTests - TestMain.cpp - ${test_source}) - -target_link_libraries(runEALTests eal ${ALLLIBS} ${GTEST_LIBRARIES} ${GTEST_MAIN_LIBRARIES} pthread) +# setup test executable +add_executable(runEALTests ${test_source}) +target_link_libraries(runEALTests eal ${GTEST_LIBRARIES} ${GTEST_MAIN_LIBRARIES} pthread) gtest_discover_tests(runEALTests WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/tests) diff --git a/tests/EalTest.cpp b/tests/EalTest.cpp new file mode 100644 index 0000000..7747784 --- /dev/null +++ b/tests/EalTest.cpp @@ -0,0 +1,144 @@ +#include "gtest/gtest.h" + +#include "eal.h" + +#include + +using namespace std; + +TEST(PositionInterpTest, LinearInterp) { + vector times = { -3, -2, -1, 0, 1, 2}; + vector> data = {{ -3, -2, -1, 0, 1, 2}, + { 9, 4, 1, 0, 1, 4}, + {-27, -8, -1, 0, 1, 8}}; + + vector coordinate = eal::getPosition(data, times, eal::linear, -1.5); + + ASSERT_EQ(3, coordinate.size()); + EXPECT_DOUBLE_EQ(-1.5, coordinate[0]); + EXPECT_DOUBLE_EQ(2.5, coordinate[1]); + EXPECT_DOUBLE_EQ(-4.5, coordinate[2]); +} + +TEST(PositionInterpTest, SplineInterp) { + vector times = {0, 1, 2, 3}; + vector> data = {{0, 0, 0, 0}, + {0, 1, 2, 3}, + {0, 2, 1, 0}}; + + vector coordinate = eal::getPosition(data, times, eal::spline, 0.5); + + ASSERT_EQ(3, coordinate.size()); + EXPECT_DOUBLE_EQ(0, coordinate[0]); + EXPECT_DOUBLE_EQ(0.5, coordinate[1]); + EXPECT_DOUBLE_EQ(2.8 * 0.5 - 0.8 * 0.125, coordinate[2]); +} + +TEST(PositionInterpTest, FourCoordinates) { + vector times = { -3, -2, -1, 0, 1, 2}; + vector> data = {{ -3, -2, -1, 0, 1, 2}, + { 9, 4, 1, 0, 1, 4}, + {-27, -8, -1, 0, 1, 8}, + { 25, 0, -5, 25, 3, 6}}; + + EXPECT_THROW(eal::getPosition(data, times, eal::linear, 0.0), + invalid_argument); +} + +TEST(PositionInterpTest, BadInterpolation) { + vector times = { -3, -2, -1, 0, 1, 2}; + vector> data = {{ -3, -2, -1, 0, 1, 2}, + { 9, 4, 1, 0, 1, 4}, + {-27, -8, -1, 0, 1, 8}}; + + EXPECT_THROW(eal::getPosition(data, times, (eal::interpolation)1000, 0.0), + invalid_argument); +} + +TEST(LinearInterpTest, ExampleInterpolation) { + vector times = {0, 1, 2, 3}; + vector data = {0, 2, 1, 0}; + + EXPECT_DOUBLE_EQ(0.0, eal::linearInterpolate(data, times, 0.0)); + EXPECT_DOUBLE_EQ(1.0, eal::linearInterpolate(data, times, 0.5)); + EXPECT_DOUBLE_EQ(2.0, eal::linearInterpolate(data, times, 1.0)); + EXPECT_DOUBLE_EQ(1.5, eal::linearInterpolate(data, times, 1.5)); + EXPECT_DOUBLE_EQ(1.0, eal::linearInterpolate(data, times, 2.0)); + EXPECT_DOUBLE_EQ(0.5, eal::linearInterpolate(data, times, 2.5)); + EXPECT_DOUBLE_EQ(0.0, eal::linearInterpolate(data, times, 3.0)); +} + +TEST(LinearInterpTest, NoPoints) { + vector times = {}; + vector data = {}; + + EXPECT_THROW(eal::linearInterpolate(data, times, 0.0), + invalid_argument); +} + +TEST(LinearInterpTest, DifferentCounts) { + vector times = { -3, -2, -1, 0, 2}; + vector data = { -3, -2, 1, 2}; + + EXPECT_THROW(eal::linearInterpolate(data, times, 0.0), + invalid_argument); +} + +TEST(LinearInterpTest, Extrapolate) { + vector times = {0, 1, 2, 3}; + vector data = {0, 2, 1, 0}; + + EXPECT_THROW(eal::linearInterpolate(data, times, -1.0), + invalid_argument); + EXPECT_THROW(eal::linearInterpolate(data, times, 4.0), + invalid_argument); +} + +TEST(SplineInterpTest, ExampleInterpolation) { + // From http://www.maths.nuigalway.ie/~niall/teaching/Archive/1617/MA378/2-2-CubicSplines.pdf + vector times = {0, 1, 2, 3}; + vector data = {0, 2, 1, 0}; + // Spline functions is: + // 2.8x - 0.8x^3, x in [0, 1] + // S(x) = x^3 - 5.4x^2 + 8.2x - 1.8, x in [1, 2] + // -0.2x^3 + 1.8x^2 - 6.2x + 7.8, x in [2, 3] + + // The spline interpolation is only ~1e-10 so we have to define a tolerance + double tolerance = 1e-10; + EXPECT_NEAR(0.0, eal::splineInterpolate(data, times, 0.0), tolerance); + EXPECT_NEAR(2.8 * 0.5 - 0.8 * 0.125, + eal::splineInterpolate(data, times, 0.5), tolerance); + EXPECT_NEAR(2.0, eal::splineInterpolate(data, times, 1.0), tolerance); + EXPECT_NEAR(3.375 - 5.4 * 2.25 + 8.2 * 1.5 - 1.8, + eal::splineInterpolate(data, times, 1.5), tolerance); + EXPECT_NEAR(1.0, eal::splineInterpolate(data, times, 2.0), tolerance); + EXPECT_NEAR(-0.2 * 15.625 + 1.8 * 6.25 - 6.2 * 2.5 + 7.8, + eal::splineInterpolate(data, times, 2.5), tolerance); + EXPECT_NEAR(0.0, eal::splineInterpolate(data, times, 3.0), tolerance); +} + +TEST(SplineInterpTest, NoPoints) { + vector times = {}; + vector data = {}; + + EXPECT_THROW(eal::splineInterpolate(data, times, 0.0), + invalid_argument); +} + +TEST(SplineInterpTest, DifferentCounts) { + vector times = { -3, -2, -1, 0, 2}; + vector data = { -3, -2, 1, 2}; + + EXPECT_THROW(eal::splineInterpolate(data, times, 0.0), + invalid_argument); +} + +TEST(SplineInterpTest, Extrapolate) { + vector times = {0, 1, 2, 3}; + vector data = {0, 2, 1, 0}; + + EXPECT_THROW(eal::splineInterpolate(data, times, -1.0), + invalid_argument); + EXPECT_THROW(eal::splineInterpolate(data, times, 4.0), + invalid_argument); +} diff --git a/tests/Test.cpp b/tests/Test.cpp deleted file mode 100644 index 7dad61b..0000000 --- a/tests/Test.cpp +++ /dev/null @@ -1,8 +0,0 @@ -#include "gtest/gtest.h" - - -TEST (Test, SimpleTest) { - int test = 0; - - ASSERT_EQ(test, 0); -} -- GitLab