Skip to content

Commit

Permalink
Add dynamic expressions
Browse files Browse the repository at this point in the history
(squashed commit)

Add dynamic expressions.
    - Add Proxy, Dynamic, RefProxy, and related classes in `dynamic` module.
    - Add tests and benchmaks
    - Add page to docs
    - Add dynamic forward and reverse Jacobian evaluators
    - Add MatrixMap container for returning Jacobians

Change storage. Expressions hold references to other expressions if their
StoreByRef trait is true (leaves are no longer special), unless they were
received as rvalues.

    - simplify arg_t and wave_ref_sel_t: store lightweight expressions by ref,
      unless T&&
    - Consistently use traits<T>::RhsDerived instead of T::RhsDerived
    - Make traits<T>::RhsDerived include the &&
    - Rework PrepareExpr (no recursion if returning a reference)
    - PreparedType trait is now T& for "no change / PrepareExpr::run returns
      reference" and T&& for "it is a new temporary expression" (e.g. for
      Convert)

Minor fixes. Update some benchmarks and CMake.
  • Loading branch information
leokoppel committed Oct 20, 2018
1 parent d0fd696 commit c90e633
Show file tree
Hide file tree
Showing 61 changed files with 2,164 additions and 517 deletions.
5 changes: 4 additions & 1 deletion CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -2,13 +2,16 @@

## Unreleased
### New features
- New `dynamic` module for heap-allocated, dynamically-composable `Proxy` expressions
- New documentation built with Sphinx

### Backward-incompatible API changes
- C++14 is now required

### Fixes and minor changes
- Fixed finding googletest source package on Ubuntu bionic
- Fix finding googletest source package on Ubuntu bionic
- Fix (trivial) reverse-mode AD on a single leaf
- Move numerical Jacobian evaluator into `core` module

## [0.3.0](https://github.com/wavelab/wave_geometry/compare/0.2.0...0.3.0) (2018-08-19)
### New features
Expand Down
1 change: 0 additions & 1 deletion benchmarks/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
find_package(benchmark REQUIRED)

set(CMAKE_VERBOSE_MAKEFILE ON)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=native")

wave_geometry_add_benchmark(util_cross_matrix_bench util_cross_matrix_bench.cpp)
wave_geometry_add_benchmark(util_identity_bench util_identity_bench.cpp)
Expand Down
5 changes: 3 additions & 2 deletions benchmarks/bechmark_helpers.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -26,9 +26,10 @@
/** Return a vector of random objects of type T, using setRandom() */
template <typename T>
std::vector<T, Eigen::aligned_allocator<T>> randomMatrices(int N) {
std::vector<T, Eigen::aligned_allocator<T>> v(N);
std::vector<T, Eigen::aligned_allocator<T>> v;
v.reserve(N);
for (auto i = N; i--;) {
v[i] = T::Random();
v.push_back(T::Random());
}
return v;
}
Expand Down
20 changes: 18 additions & 2 deletions benchmarks/rotate_chain/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -3,14 +3,30 @@
#
# Each file tests a different implementation (ceres autodiff, hand, wave) for different N.

find_package(GTSAM QUIET)
if(GTSAM_FOUND)
add_executable(rotate_chain_gtsam_bench rotate_chain_gtsam_bench.cpp)
target_link_libraries(rotate_chain_gtsam_bench gtsam benchmark::benchmark GTest::GTest)
set_target_properties(rotate_chain_gtsam_bench PROPERTIES
RUNTIME_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/benchmarks)
target_include_directories(rotate_chain_gtsam_bench PRIVATE ${PROJECT_SOURCE_DIR}/include ${GTSAM_INCLUDE_DIR})
endif()

# Use C++17 for structured bindings
# Rotate_chain uses C++17
list (FIND CMAKE_CXX_KNOWN_FEATURES "cxx_std_17" index)
list (FIND CMAKE_CXX_COMPILE_FEATURES "cxx_std_17" index)
if (${index} GREATER -1)
set(CMAKE_CXX_STANDARD 17)
else()
message(STATUS "Not configuring rotate_chain benchmarks because C++17 is not known to this version of CMake")
message(STATUS "Not configuring rotate_chain benchmarks because (according to CMake) compiler does not support C++17")
return()
endif()

wave_geometry_add_benchmark(rotate_chain_hand_bench rotate_chain_hand_bench.cpp)
wave_geometry_add_benchmark(rotate_chain_wave_bench rotate_chain_wave_bench.cpp)
wave_geometry_add_benchmark(rotate_chain_wave_untyped_bench rotate_chain_wave_untyped_bench.cpp)
wave_geometry_add_benchmark(rotate_chain_wave_reverse_bench rotate_chain_wave_reverse_bench.cpp)
wave_geometry_add_benchmark(rotate_chain_wave_dynamic_bench rotate_chain_wave_dynamic_bench.cpp)


wave_geometry_add_benchmark(imu_preint imu_preint.cpp)
291 changes: 62 additions & 229 deletions benchmarks/rotate_chain/rotate_chain_gtsam_bench.cpp
Original file line number Diff line number Diff line change
@@ -1,253 +1,86 @@
#include <benchmark/benchmark.h>
#include "../bechmark_helpers.hpp"
#include <Eigen/Core>
// Some adaptions for gtsam's bundled Eigen being old
#if !(EIGEN_VERSION_AT_LEAST(3, 3, 0))
#define EIGEN_DEVICE_FUNC
namespace Eigen {
using Index = EIGEN_DEFAULT_DENSE_INDEX_TYPE;
}
#endif

#include <gtsam/slam/expressions.h>
#include <benchmark/benchmark.h>
#include <gtsam/nonlinear/expressionTesting.h>


#include <gtsam/slam/expressions.h>
#include "wave/geometry/src/util/math/math.hpp"


template <typename T>
using EigenVector = std::vector<T, Eigen::aligned_allocator<T>>;


using gtsam::Rot3;
using gtsam::Point3;
using gtsam::Expression;
using gtsam::Point3;
using gtsam::Rot3;

/** Return a vector of random gtsam rotations */
EigenVector<Rot3> randomRot3s(int N) {
EigenVector<Rot3> v;
v.reserve(N);
for (auto i = N; i--;) {
const auto &q = wave::randomQuaternion<double>();
v.emplace_back(q);
}
return v;
}

inline Eigen::Matrix3d crossMatrix(const Eigen::Vector3d &v) {
Eigen::Matrix3d m;
m << 0.0, -v.z(), v.y(), //
v.z(), 0.0, -v.x(), //
-v.y(), v.x(), 0.0;
return m;
};

class RotateChain : public benchmark::Fixture {
protected:
const int N = 1000;
std::vector<gtsam::Values> vv;
const EigenVector<Rot3> R1 = randomRot3s(N);
const EigenVector<Rot3> R2 = randomRot3s(N);
const EigenVector<Rot3> R3 = randomRot3s(N);
const EigenVector<Rot3> R4 = randomRot3s(N);
const EigenVector<Rot3> R5 = randomRot3s(N);
const EigenVector<Rot3> R6 = randomRot3s(N);
const EigenVector<Rot3> R7 = randomRot3s(N);
const EigenVector<Rot3> R8 = randomRot3s(N);
const EigenVector<Rot3> R9 = randomRot3s(N);
const EigenVector<Rot3> R10 = randomRot3s(N);
const EigenVector<Point3> v1 = randomMatrices<Point3>(N);

RotateChain() {
for (auto i = 0; i < N; ++i) {
gtsam::Values values{};
values.insert(gtsam::Symbol{'R', 1}, Rot3{wave::randomQuaternion<double>()});
values.insert(gtsam::Symbol{'R', 2}, Rot3{wave::randomQuaternion<double>()});
values.insert(gtsam::Symbol{'R', 3}, Rot3{wave::randomQuaternion<double>()});
values.insert(gtsam::Symbol{'R', 4}, Rot3{wave::randomQuaternion<double>()});
values.insert(gtsam::Symbol{'R', 5}, Rot3{wave::randomQuaternion<double>()});
values.insert(gtsam::Symbol{'R', 6}, Rot3{wave::randomQuaternion<double>()});
values.insert(gtsam::Symbol{'R', 7}, Rot3{wave::randomQuaternion<double>()});
values.insert(gtsam::Symbol{'R', 8}, Rot3{wave::randomQuaternion<double>()});
values.insert(gtsam::Symbol{'R', 9}, Rot3{wave::randomQuaternion<double>()});
values.insert(gtsam::Symbol{'R', 10}, Rot3{wave::randomQuaternion<double>()});
values.insert(gtsam::Symbol{'p', 1}, Point3{Eigen::Vector3d::Random()});
vv.push_back(values);
}
}
};

BENCHMARK_F(RotateChain, Gtsam1)(benchmark::State &state) {
for (auto _ : state) {
Expression<Rot3> R1_{'R', 1};
Expression<Point3> p1_{'p', 1};
Expression<Point3> p2_ = rotate(R1_, p1_);
std::vector<gtsam::Matrix> jacs(2);

for (auto i = N; i-- > 0;) {
Point3 p2 = p2_.value(vv[i], jacs);

benchmark::DoNotOptimize(p2);
benchmark::DoNotOptimize(jacs.data());
}
}
}

BENCHMARK_F(RotateChain, Gtsam2)(benchmark::State &state) {
const Expression<Rot3> R1_{Rot3{}}, R2_{'R', 2};
const Expression<Point3> p1_{'p', 1};
const Expression<Point3> p2_ = rotate(R1_ * R2_, p1_);

std::vector<gtsam::Matrix> jacs(3);

for (auto _ : state) {
for (auto i = N; i-- > 0;) {
Point3 p2 = p2_.value(vv[i], jacs);

benchmark::DoNotOptimize(p2);
benchmark::DoNotOptimize(jacs);
}
}
}

BENCHMARK_F(RotateChain, Gtsam3)(benchmark::State &state) {
const Expression<Rot3> R1_{'R', 1}, R2_{'R', 2}, R3_{'R', 3};
const Expression<Point3> p1_{'p', 1};
const Expression<Point3> p2_ = rotate(R1_ * R2_ * R3_, p1_);

std::vector<gtsam::Matrix> jacs(4);

for (auto _ : state) {
for (auto i = N; i-- > 0;) {
Point3 p2 = p2_.value(vv[i], jacs);

benchmark::DoNotOptimize(p2);
benchmark::DoNotOptimize(jacs);
}
}
}

BENCHMARK_F(RotateChain, Gtsam4)(benchmark::State &state) {
const Expression<Rot3> R1_{'R', 1}, R2_{'R', 2}, R3_{'R', 3}, R4_{'R', 4};
const Expression<Point3> p1_{'p', 1};
const Expression<Point3> p2_ = rotate(R1_ * R2_ * R3_ * R4_, p1_);

std::vector<gtsam::Matrix> jacs(5);

for (auto _ : state) {
for (auto i = N; i-- > 0;) {
Point3 p2 = p2_.value(vv[i], jacs);

benchmark::DoNotOptimize(p2);
benchmark::DoNotOptimize(jacs);
}
}
}

BENCHMARK_F(RotateChain, Gtsam5)(benchmark::State &state) {
const Expression<Rot3> R1_{'R', 1}, R2_{'R', 2}, R3_{'R', 3}, R4_{'R', 4},
R5_{'R', 5};
const Expression<Point3> p1_{'p', 1};
const Expression<Point3> p2_ = rotate(R1_ * R2_ * R3_ * R4_ * R5_, p1_);

std::vector<gtsam::Matrix> jacs(6);

for (auto _ : state) {
for (auto i = N; i-- > 0;) {
Point3 p2 = p2_.value(vv[i], jacs);

benchmark::DoNotOptimize(p2);
benchmark::DoNotOptimize(jacs);
}
}
}

BENCHMARK_F(RotateChain, Gtsam6)(benchmark::State &state) {
const Expression<Rot3> R1_{'R', 1}, R2_{'R', 2}, R3_{'R', 3}, R4_{'R', 4},
R5_{'R', 5}, R6_{'R', 6};
const Expression<Point3> p1_{'p', 1};
const Expression<Point3> p2_ = rotate(R1_ * R2_ * R3_ * R4_ * R5_ * R6_, p1_);

std::vector<gtsam::Matrix> jacs(7);

for (auto _ : state) {
for (auto i = N; i-- > 0;) {
Point3 p2 = p2_.value(vv[i], jacs);

benchmark::DoNotOptimize(p2);
benchmark::DoNotOptimize(jacs);
}
}
}

BENCHMARK_F(RotateChain, Gtsam7)(benchmark::State &state) {
const Expression<Rot3> R1_{'R', 1}, R2_{'R', 2}, R3_{'R', 3}, R4_{'R', 4},
R5_{'R', 5}, R6_{'R', 6}, R7_{'R', 7};
const Expression<Point3> p1_{'p', 1};
const Expression<Point3> p2_ = rotate(R1_ * R2_ * R3_ * R4_ * R5_ * R6_ * R7_, p1_);

std::vector<gtsam::Matrix> jacs(8);

for (auto _ : state) {
for (auto i = N; i-- > 0;) {
Point3 p2 = p2_.value(vv[i], jacs);

benchmark::DoNotOptimize(p2);
benchmark::DoNotOptimize(jacs);
}
void BM_gtsamAll(benchmark::State &state) {
const auto N = state.range(0);
state.SetComplexityN(N);
// Produce the expression tree
gtsam::Values values;
auto p1 = Expression<Point3>{'p', 1};
auto expr = p1;
values.insert(gtsam::Symbol{'p', 1}, Point3{Eigen::Vector3d::Random()});

for (auto i = 0u; i < N; ++i) {
expr = Expression<Point3>{rotate(Expression<Rot3>{'R', i}, expr)};
values.insert(gtsam::Symbol{'R', i}, Rot3{wave::randomQuaternion<double>()});
}
}

BENCHMARK_F(RotateChain, Gtsam8)(benchmark::State &state) {
const Expression<Rot3> R1_{'R', 1}, R2_{'R', 2}, R3_{'R', 3}, R4_{'R', 4},
R5_{'R', 5}, R6_{'R', 6}, R7_{'R', 7}, R8_{'R', 8};
const Expression<Point3> p1_{'p', 1};
const Expression<Point3> p2_ =
rotate(R1_ * R2_ * R3_ * R4_ * R5_ * R6_ * R7_ * R8_, p1_);

std::vector<gtsam::Matrix> jacs(9);
std::vector<gtsam::Matrix> jacs(N + 1);

for (auto _ : state) {
for (auto i = N; i-- > 0;) {
Point3 p2 = p2_.value(vv[i], jacs);
Point3 p2 = expr.value(values, jacs);

benchmark::DoNotOptimize(p2);
benchmark::DoNotOptimize(jacs);
}
benchmark::DoNotOptimize(p2);
benchmark::DoNotOptimize(jacs);
}
}


BENCHMARK_F(RotateChain, Gtsam9)(benchmark::State &state) {
const Expression<Rot3> R1_{'R', 1}, R2_{'R', 2}, R3_{'R', 3}, R4_{'R', 4},
R5_{'R', 5}, R6_{'R', 6}, R7_{'R', 7}, R8_{'R', 8}, R9_{'R', 9};
const Expression<Point3> p1_{'p', 1};
const Expression<Point3> p2_ =
rotate(R1_ * R2_ * R3_ * R4_ * R5_ * R6_ * R7_ * R8_ * R9_, p1_);

std::vector<gtsam::Matrix> jacs(10);

for (auto _ : state) {
for (auto i = N; i-- > 0;) {
Point3 p2 = p2_.value(vv[i], jacs);

benchmark::DoNotOptimize(p2);
benchmark::DoNotOptimize(jacs);
}
// This benchmark tests gtsam's valueAndJacobianMap() function, which is normally private
#ifdef LK_CUSTOM_GTSAM_PUBLIC
void BM_gtsamPrivate(benchmark::State &state) {
const auto N = state.range(0);
state.SetComplexityN(N);
// Produce the expression tree
gtsam::Values values;
auto p1 = Expression<Point3>{'p', 1};
auto expr = p1;
values.insert(gtsam::Symbol{'p', 1}, Point3{Eigen::Vector3d::Random()});

for (auto i = 0u; i < N; ++i) {
expr = Expression<Point3>{rotate(Expression<Rot3>{'R', i}, expr)};
values.insert(gtsam::Symbol{'R', i}, Rot3{wave::randomQuaternion<double>()});
}
}


BENCHMARK_F(RotateChain, Gtsam10)(benchmark::State &state) {
const Expression<Rot3> R1_{'R', 1}, R2_{'R', 2}, R3_{'R', 3}, R4_{'R', 4},
R5_{'R', 5}, R6_{'R', 6}, R7_{'R', 7}, R8_{'R', 8}, R9_{'R', 9}, R10_{'R', 10};
const Expression<Point3> p1_{'p', 1};
const Expression<Point3> p2_ =
rotate(R1_ * R2_ * R3_ * R4_ * R5_ * R6_ * R7_ * R8_ * R9_ * R10_, p1_);

std::vector<gtsam::Matrix> jacs(11);
std::vector<gtsam::Matrix> jacs(N + 1);

for (auto _ : state) {
for (auto i = N; i-- > 0;) {
Point3 p2 = p2_.value(vv[i], jacs);

benchmark::DoNotOptimize(p2);
benchmark::DoNotOptimize(jacs);
}
// Pre-allocate and zero VerticalBlockMatrix
gtsam::KeyVector keys;
gtsam::FastVector<int> dims;
boost::tie(keys, dims) = expr.keysAndDims();
static const int Dim = gtsam::traits<Point3>::dimension;
gtsam::VerticalBlockMatrix Ab(dims, Dim);
Ab.matrix().setZero();
gtsam::internal::JacobianMap jacobianMap(keys, Ab);

// Call unsafe version
Point3 p2 = expr.valueAndJacobianMap(values, jacobianMap);

benchmark::DoNotOptimize(p2);
benchmark::DoNotOptimize(values);
benchmark::DoNotOptimize(jacobianMap);
}
}
BENCHMARK(BM_gtsamPrivate)->RangeMultiplier(2)->Range(1, 1 << 14)->Complexity();
#endif

WAVE_BENCHMARK_MAIN()
BENCHMARK(BM_gtsamAll)->RangeMultiplier(2)->Range(1, 1 << 14)->Complexity();
BENCHMARK_MAIN();
2 changes: 1 addition & 1 deletion benchmarks/rotate_chain/rotate_chain_hand_bench.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ using EigenVector = std::vector<T, Eigen::aligned_allocator<T>>;
template <class RotType>
class RotateChain : public benchmark::Fixture {
protected:
const int N = 1000;
const int N = 1;
// Use wave::RotationM just for the Random() method, which produces valid SO(3)
using waveRot = wave::MatrixRotation<RotType>;
const EigenVector<waveRot> R1 = randomMatrices<waveRot>(N);
Expand Down
Loading

0 comments on commit c90e633

Please sign in to comment.