From 5597f65ad582b27ff09bb581b9fbcd7aae22470f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Benjamin=20Chr=C3=A9tien?= Date: Tue, 7 Apr 2015 01:03:08 +0900 Subject: [PATCH 1/3] Fix compilation for both ColMajor/RowMajor + provide a define. Users can use the CMake option "STORAGE_ORDER" to change the storage order used by RobOptim. It defaults to ColMajor (Eigen's default), and RowMajor can be used for hypothetical increased performance in some scenarios (e.g. when mostly using impl_gradient to fill the Jacobian matrices). --- CMakeLists.txt | 13 +++ include/roboptim/core/detail/utility.hh | 14 ++++ include/roboptim/core/function.hh | 79 ++++++++++++++++--- include/roboptim/core/operator/concatenate.hh | 27 +++++++ .../roboptim/core/operator/concatenate.hxx | 51 +++++++++++- include/roboptim/core/operator/product.hxx | 13 ++- src/visualization/gnuplot-matrix.cc | 33 +------- tests/derivable-function.cc | 22 +++--- tests/detail-utility.cc | 44 +++++++++++ tests/function.cc | 3 + tests/operator-scalar.cc | 12 +-- ...ion-gnuplot-differentiable-function.stdout | 14 ++-- 12 files changed, 249 insertions(+), 76 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 5053cd0cf..a1b65f722 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -181,6 +181,19 @@ PKG_CONFIG_APPEND_LIBS(ltdl) PKG_CONFIG_APPEND_CFLAGS (-DEIGEN_RUNTIME_NO_MALLOC) HEADER_INSTALL("${HEADERS}") +# Set the default storage order for Eigen matrices (ColMajor or RowMajor). +# Note: RowMajor may yield better results for the dense Jacobian filling process +# if most functions only implement impl_gradient, but ColMajor is the usual +# default storage order for this kind of scenario. +SET(STORAGE_ORDER "ColMajor" CACHE + STRING "Storage order for Eigen matrices (ColMajor or RowMajor") +IF(NOT (STORAGE_ORDER STREQUAL "ColMajor" OR STORAGE_ORDER STREQUAL "RowMajor")) + MESSAGE(FATAL_ERROR "Storage order should be either ColMajor or RowMajor.") +ENDIF() + +SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DROBOPTIM_STORAGE_ORDER=${STORAGE_ORDER}") +PKG_CONFIG_APPEND_CFLAGS (-DROBOPTIM_STORAGE_ORDER=${STORAGE_ORDER}) + OPTION(DISABLE_TESTS "Disable test programs" OFF) ADD_SUBDIRECTORY(src) diff --git a/include/roboptim/core/detail/utility.hh b/include/roboptim/core/detail/utility.hh index 0435847df..10aee91ac 100644 --- a/include/roboptim/core/detail/utility.hh +++ b/include/roboptim/core/detail/utility.hh @@ -136,6 +136,20 @@ namespace roboptim }; + /// \brief Get the matrix stride type for a row vector, given a matrix + /// storage order. + /// + /// This solves compiler errors in Eigen when dealing with Eigen::Refs, + /// row vectors and different storage orders. + /// + /// \tparam SO storage order (Eigen::ColMajor or Eigen::RowMajor). + template + struct row_vector_stride + { + typedef Eigen::InnerStride<(SO == Eigen::RowMajor)? 1:-1> type; + }; + + /// \brief Converts CLIST to a boost::mpl::vector to ensure a similar /// behavior for codes using different random access sequences (vector, /// list, etc.). diff --git a/include/roboptim/core/function.hh b/include/roboptim/core/function.hh index 14f7888fe..59750e1bc 100644 --- a/include/roboptim/core/function.hh +++ b/include/roboptim/core/function.hh @@ -27,6 +27,7 @@ # include # include +# include # include # include # include @@ -43,11 +44,21 @@ # include # include +# include + +// Generate Eigen::Ref types for a given data type. # define ROBOPTIM_GENERATE_TYPEDEFS_EIGEN_REF(NAME,TYPE) \ typedef TYPE NAME##_t; \ typedef Eigen::Ref NAME##_ref; \ typedef const Eigen::Ref& const_##NAME##_ref +// Generate Eigen::Ref types for a given data type. Specialized version for row +// vectors with the proper stride, no matter what the storage order is. +# define ROBOPTIM_GENERATE_TYPEDEFS_EIGEN_REF_VEC(NAME,TYPE) \ + typedef TYPE NAME##_t; \ + typedef Eigen::Ref::type> NAME##_ref; \ + typedef const Eigen::Ref::type>& const_##NAME##_ref + # define ROBOPTIM_GENERATE_TYPEDEFS_REF(NAME,TYPE) \ typedef TYPE NAME##_t; \ typedef NAME##_t& NAME##_ref; \ @@ -76,6 +87,7 @@ ROBOPTIM_GENERATE_FWD_REFS(argument); \ ROBOPTIM_GENERATE_FWD_REFS(result); \ ROBOPTIM_GENERATE_FWD_REFS(vector); \ + ROBOPTIM_GENERATE_FWD_REFS(rowVector); \ ROBOPTIM_GENERATE_FWD_REFS(matrix) # define ROBOPTIM_FUNCTION_FWD_TYPEDEFS_(PARENT) \ @@ -86,13 +98,26 @@ ROBOPTIM_GENERATE_FWD_REFS_(argument); \ ROBOPTIM_GENERATE_FWD_REFS_(result); \ ROBOPTIM_GENERATE_FWD_REFS_(vector); \ + ROBOPTIM_GENERATE_FWD_REFS_(rowVector); \ ROBOPTIM_GENERATE_FWD_REFS_(matrix) +// Default storage order = ColMajor +# ifndef ROBOPTIM_STORAGE_ORDER +# define ROBOPTIM_STORAGE_ORDER ColMajor +# endif //! ROBOPTIM_STORAGE_ORDER + +BOOST_STATIC_ASSERT_MSG (Eigen::ROBOPTIM_STORAGE_ORDER == Eigen::ColMajor \ + || Eigen::ROBOPTIM_STORAGE_ORDER == Eigen::RowMajor, \ + "Wrong storage order provided by ROBOPTIM_STORAGE_ORDER."); + namespace roboptim { /// \addtogroup roboptim_meta_function /// @{ + /// \brief Default matrix storage order. + static const int StorageOrder = Eigen::ROBOPTIM_STORAGE_ORDER; + /// \brief GenericFunction traits /// /// This helper class is used to define which types a GenericFunction @@ -141,7 +166,7 @@ namespace roboptim /// used for computations. typedef typename GenericFunctionTraits::value_type value_type; - /// \brief Basic vector type. + /// \brief Basic (column) vector type. /// /// This basic vector type is used each time a vector of values /// is required. @@ -151,6 +176,12 @@ namespace roboptim /// implementation. ROBOPTIM_GENERATE_TRAITS_REFS_(vector); + /// \brief Row vector type. + /// + /// This basic vector type is used each time a row vector of values + /// is required (e.g. gradients). + ROBOPTIM_GENERATE_TRAITS_REFS_(rowVector); + /// \brief Basic matrix type. /// /// This basic matrix type is used each time a two dimensional @@ -551,30 +582,43 @@ namespace roboptim struct GenericFunctionTraits { /// \brief Matrix storage order. - static const int StorageOrder = Eigen::RowMajor; + static const int StorageOrder = roboptim::StorageOrder; + + /// \brief Value type. + typedef double value_type; // For each type, we have: // - type_t: the type itself // - type_ref: reference to type object // - const_type_ref: const reference to type object + + // Matrix types ROBOPTIM_GENERATE_TYPEDEFS_EIGEN_REF (matrix, - Eigen::Matrix); + + // Vector types ROBOPTIM_GENERATE_TYPEDEFS_EIGEN_REF (vector, - Eigen::Matrix); + // Row vector types + ROBOPTIM_GENERATE_TYPEDEFS_EIGEN_REF + (rowVector, + Eigen::Matrix); + typedef matrix_t::Index size_type; - typedef matrix_t::Scalar value_type; ROBOPTIM_GENERATE_TYPEDEFS_EIGEN_REF(result,vector_t); ROBOPTIM_GENERATE_TYPEDEFS_EIGEN_REF(argument,vector_t); - ROBOPTIM_GENERATE_TYPEDEFS_EIGEN_REF(gradient,vector_t); + ROBOPTIM_GENERATE_TYPEDEFS_EIGEN_REF_VEC(gradient,rowVector_t); ROBOPTIM_GENERATE_TYPEDEFS_EIGEN_REF(jacobian,matrix_t); ROBOPTIM_GENERATE_TYPEDEFS_EIGEN_REF(hessian,matrix_t); }; @@ -584,28 +628,41 @@ namespace roboptim struct GenericFunctionTraits { /// \brief Matrix storage order. - static const int StorageOrder = Eigen::RowMajor; + static const int StorageOrder = roboptim::StorageOrder; + + /// \brief Value type. + typedef double value_type; // For each type, we have: // - type_t: the type itself // - type_ref: reference to type object // - const_type_ref: const reference to type object + + // Matrix types ROBOPTIM_GENERATE_TYPEDEFS_REF (matrix, - Eigen::SparseMatrix); + Eigen::SparseMatrix); + // Vector types ROBOPTIM_GENERATE_TYPEDEFS_EIGEN_REF (vector, - Eigen::Matrix); + // Row vector types + ROBOPTIM_GENERATE_TYPEDEFS_EIGEN_REF + (rowVector, + Eigen::Matrix); + typedef matrix_t::Index size_type; - typedef matrix_t::Scalar value_type; ROBOPTIM_GENERATE_TYPEDEFS_EIGEN_REF(result,vector_t); ROBOPTIM_GENERATE_TYPEDEFS_EIGEN_REF(argument,vector_t); - ROBOPTIM_GENERATE_TYPEDEFS_REF(gradient,Eigen::SparseVector); + ROBOPTIM_GENERATE_TYPEDEFS_REF(gradient,Eigen::SparseVector); ROBOPTIM_GENERATE_TYPEDEFS_REF(jacobian,matrix_t); ROBOPTIM_GENERATE_TYPEDEFS_REF(hessian,matrix_t); }; diff --git a/include/roboptim/core/operator/concatenate.hh b/include/roboptim/core/operator/concatenate.hh index 47b1a24b3..0e554a4f3 100644 --- a/include/roboptim/core/operator/concatenate.hh +++ b/include/roboptim/core/operator/concatenate.hh @@ -19,7 +19,10 @@ # define ROBOPTIM_CORE_OPERATOR_CONCATENATE_HH # include # include + # include +# include +# include # include # include @@ -27,6 +30,12 @@ namespace roboptim { + namespace detail + { + template + struct ConcatenateTypes; + } // end of namespace detail + /// \addtogroup roboptim_operator /// @{ @@ -39,6 +48,9 @@ namespace roboptim typedef typename detail::AutopromoteTrait::T_type parentType_t; ROBOPTIM_DIFFERENTIABLE_FUNCTION_FWD_TYPEDEFS_ (parentType_t); + /// \brief Traits type. + typedef typename parentType_t::traits_t traits_t; + typedef boost::shared_ptr ConcatenateShPtr_t; explicit Concatenate (boost::shared_ptr left, @@ -76,6 +88,21 @@ namespace roboptim void impl_jacobian (jacobian_ref jacobian, const_argument_ref arg) const; + + private: + + template + void concatenateJacobian (jacobian_ref jacobian, + const_argument_ref argument, + typename detail::ConcatenateTypes::isDense_t::type* = 0) + const; + + template + void concatenateJacobian (jacobian_ref jacobian, + const_argument_ref argument, + typename detail::ConcatenateTypes::isNotDense_t::type* = 0) + const; + private: boost::shared_ptr left_; boost::shared_ptr right_; diff --git a/include/roboptim/core/operator/concatenate.hxx b/include/roboptim/core/operator/concatenate.hxx index 2eaba78fa..2dfc1ebe2 100644 --- a/include/roboptim/core/operator/concatenate.hxx +++ b/include/roboptim/core/operator/concatenate.hxx @@ -17,10 +17,31 @@ #ifndef ROBOPTIM_CORE_OPERATOR_CONCATENATE_HXX # define ROBOPTIM_CORE_OPERATOR_CONCATENATE_HXX + # include +# include + +# include namespace roboptim { + namespace detail + { + template + struct ConcatenateTypes + { + typedef typename T::traits_t traits_t; + + typedef typename boost::enable_if< + boost::is_same > + isDense_t; + + typedef typename boost::disable_if< + boost::is_same > + isNotDense_t; + }; + } // end of namespace detail + template Concatenate::Concatenate (boost::shared_ptr left, @@ -80,9 +101,10 @@ namespace roboptim } template - void - Concatenate::impl_jacobian (jacobian_ref jacobian, - const_argument_ref x) + template + void Concatenate::concatenateJacobian (jacobian_ref jacobian, + const_argument_ref x, + typename detail::ConcatenateTypes::isDense_t::type*) const { left_->jacobian (jacobianLeft_, x); @@ -92,6 +114,29 @@ namespace roboptim jacobian.middleRows (left_->outputSize (), right_->outputSize ()) = jacobianRight_; } + + template + template + void Concatenate::concatenateJacobian (jacobian_ref jacobian, + const_argument_ref x, + typename detail::ConcatenateTypes::isNotDense_t::type*) + const + { + left_->jacobian (jacobianLeft_, x); + right_->jacobian (jacobianRight_, x); + + copySparseBlock (jacobian, jacobianLeft_, 0, 0); + copySparseBlock (jacobian, jacobianRight_, left_->outputSize (), 0); + } + + template + void + Concatenate::impl_jacobian (jacobian_ref jacobian, + const_argument_ref x) + const + { + concatenateJacobian (jacobian, x); + } } // end of namespace roboptim. #endif //! ROBOPTIM_CORE_OPERATOR_CONCATENATE_HXX diff --git a/include/roboptim/core/operator/product.hxx b/include/roboptim/core/operator/product.hxx index 23fa6644c..357d377d7 100644 --- a/include/roboptim/core/operator/product.hxx +++ b/include/roboptim/core/operator/product.hxx @@ -48,6 +48,11 @@ namespace roboptim typedef typename U::vector_ref vectorU_ref; typedef typename V::vector_ref vectorV_ref; + typedef typename U::rowVector_t rowVectorU_t; + typedef typename V::rowVector_t rowVectorV_t; + typedef typename U::rowVector_ref rowVectorU_ref; + typedef typename V::rowVector_ref rowVectorV_ref; + typedef typename U::gradient_t gradientU_t; typedef typename V::gradient_t gradientV_t; typedef typename Product::gradient_t gradient_t; @@ -68,8 +73,8 @@ namespace roboptim template static void gradient (typename Types::gradient_ref grad_uv, - const typename Types::vectorU_ref u, - const typename Types::vectorV_ref v, + const typename Types::rowVectorU_ref u, + const typename Types::rowVectorV_ref v, const typename Types::gradientU_ref grad_u, const typename Types::gradientV_ref grad_v, typename boost::enable_if::fullDense_t>::type* = 0) @@ -84,8 +89,8 @@ namespace roboptim template static void gradient (typename Types::gradient_ref grad_uv, - const typename Types::vectorU_ref u, - const typename Types::vectorV_ref v, + const typename Types::rowVectorU_ref u, + const typename Types::rowVectorV_ref v, const typename Types::gradientU_ref grad_u, const typename Types::gradientV_ref grad_v, typename boost::disable_if::fullDense_t>::type* = 0) diff --git a/src/visualization/gnuplot-matrix.cc b/src/visualization/gnuplot-matrix.cc index 548b13a81..fc375eb6f 100644 --- a/src/visualization/gnuplot-matrix.cc +++ b/src/visualization/gnuplot-matrix.cc @@ -86,38 +86,9 @@ namespace roboptim std::string sparse_matrix_to_gnuplot (GenericFunctionTraits::const_matrix_ref mat) { - typedef GenericFunctionTraits::matrix_t matrix_t; + typedef GenericFunctionTraits::matrix_t denseMatrix_t; - std::string str = ""; - - // Set the header of the Gnuplot output (range, etc.) - set_matrix_header (str, mat); - - // Since Gnuplot does not support sparse matrices, we will need to - // plot all the zeros of the sparse matrices. - // Outer dimension - for (int k = 0; k < mat.outerSize (); ++k) - { - // Inner dimension - matrix_t::InnerIterator it (mat,k); - for (int col_it = 0; col_it < mat.cols (); ++col_it) - { - // Sparse nonzero: return 1 - if (col_it == it.col ()) - { - str += "1"; - ++it; - } - // Sparse zero: return 0 - else str += "0"; - - if (col_it < mat.cols () - 1) str += " "; - else str += "\n"; - } - } - str += "e\n"; - str += "e\n"; - return str; + return dense_matrix_to_gnuplot (denseMatrix_t (mat)); } } // end of namespace detail. diff --git a/tests/derivable-function.cc b/tests/derivable-function.cc index 77d8755e5..1bc0c5a42 100644 --- a/tests/derivable-function.cc +++ b/tests/derivable-function.cc @@ -224,13 +224,15 @@ BOOST_AUTO_TEST_CASE_TEMPLATE (jacobian_check, T, functionTypes_t) x[2] = 7.; x[3] = 2.; - BOOST_REQUIRE_EQUAL (f.jacobian (x).row (0).size (), + // TODO: use row() once Eigen is fixed + // cf. http://eigen.tuxfamily.org/bz/show_bug.cgi?id=980 + BOOST_REQUIRE_EQUAL (f.jacobian (x).block (0,0,1,4).size (), f.gradient (x, 0).size ()); - BOOST_REQUIRE_EQUAL (f.jacobian (x).row (1).size (), + BOOST_REQUIRE_EQUAL (f.jacobian (x).block (1,0,1,4).size (), f.gradient (x, 1).size ()); - BOOST_CHECK (f.jacobian (x).row (0).isApprox (f.gradient (x, 0).adjoint ())); - BOOST_CHECK (f.jacobian (x).row (1).isApprox (f.gradient (x, 1).adjoint ())); + BOOST_CHECK (f.jacobian (x).block (0,0,1,4).isApprox (f.gradient (x, 0))); + BOOST_CHECK (f.jacobian (x).block (1,0,1,4).isApprox (f.gradient (x, 1))); typename F::jacobian_t jac = f.jacobian (x); typename F::gradient_t grad0 = f.gradient (x, 0); @@ -239,8 +241,8 @@ BOOST_AUTO_TEST_CASE_TEMPLATE (jacobian_check, T, functionTypes_t) // Call f with a map to a C-style array. typename F::value_type ar[4] = {10., -1., 7., 2.}; Eigen::Map::argument_t> map_x (ar, 4); - BOOST_CHECK (f.jacobian (map_x).row (0).isApprox (f.gradient (map_x, 0).adjoint ())); - BOOST_CHECK (f.jacobian (map_x).row (1).isApprox (f.gradient (map_x, 1).adjoint ())); + BOOST_CHECK (f.jacobian (map_x).block (0,0,1,4).isApprox (f.gradient (map_x, 0))); + BOOST_CHECK (f.jacobian (map_x).block (1,0,1,4).isApprox (f.gradient (map_x, 1))); BOOST_CHECK (f.jacobian (map_x).isApprox (jac)); BOOST_CHECK (f.gradient (map_x, 0).isApprox (grad0)); BOOST_CHECK (f.gradient (map_x, 1).isApprox (grad1)); @@ -250,10 +252,10 @@ BOOST_AUTO_TEST_CASE_TEMPLATE (jacobian_check, T, functionTypes_t) typename F::argument_t full_x (8); full_x.setZero (); full_x.segment (2,4) = x; - BOOST_CHECK (f.jacobian (full_x.segment (2,4)).row (0).isApprox - (f.gradient (full_x.segment (2,4), 0).adjoint ())); - BOOST_CHECK (f.jacobian (full_x.segment (2,4)).row (1).isApprox - (f.gradient (full_x.segment (2,4), 1).adjoint ())); + BOOST_CHECK (f.jacobian (full_x.segment (2,4)).block (0,0,1,4).isApprox + (f.gradient (full_x.segment (2,4), 0))); + BOOST_CHECK (f.jacobian (full_x.segment (2,4)).block (1,0,1,4).isApprox + (f.gradient (full_x.segment (2,4), 1))); BOOST_CHECK (f.jacobian (full_x.segment (2,4)).isApprox (jac)); BOOST_CHECK (f.gradient (full_x.segment (2,4), 0).isApprox (grad0)); BOOST_CHECK (f.gradient (full_x.segment (2,4), 1).isApprox (grad1)); diff --git a/tests/detail-utility.cc b/tests/detail-utility.cc index 83a66b6b1..d624855d2 100644 --- a/tests/detail-utility.cc +++ b/tests/detail-utility.cc @@ -15,6 +15,9 @@ // You should have received a copy of the GNU Lesser General Public License // along with roboptim. If not, see . +#define EIGEN_RUNTIME_NO_MALLOC +#include + #include "shared-tests/fixture.hh" #include @@ -33,6 +36,38 @@ BOOST_FIXTURE_TEST_SUITE (core, TestSuiteConfiguration) BOOST_MPL_ASSERT ((PRED)); \ BOOST_CHECK ((PRED::type::value)); +namespace detail +{ + using namespace Eigen; + using namespace ::roboptim::detail; + + template + struct TestRowVectorStride + { + typedef Matrix matrix_t; + typedef Matrix rowVector_t; + typedef Ref::type> rowVector_ref; + + void foo (rowVector_ref v) + { + v.fill (0.); + } + + void run () + { + matrix_t m (2, 2); + m << 1., 2., 3., 4.; + + internal::set_is_malloc_allowed (false); + foo (m.row (1)); + internal::set_is_malloc_allowed (true); + + BOOST_CHECK (!m.row (0).isZero ()); + BOOST_CHECK (m.row (1).isZero ()); + } + }; +} // end of namespace detail + BOOST_AUTO_TEST_CASE (detail_utility) { @@ -91,6 +126,15 @@ BOOST_AUTO_TEST_CASE (detail_utility) TEST_PREDICATE (predicate3_t); } + // Test row_vector_stride. + { + detail::TestRowVectorStride testCol; + testCol.run (); + + detail::TestRowVectorStride testRow; + testRow.run (); + } + // Test contains_base_of. { typedef vector testSeq_t; diff --git a/tests/function.cc b/tests/function.cc index 56304bfe6..b38c5c104 100644 --- a/tests/function.cc +++ b/tests/function.cc @@ -72,6 +72,9 @@ BOOST_AUTO_TEST_CASE_TEMPLATE (null_function, T, functionTypes_t) boost::shared_ptr output = retrievePattern ("function"); + // Check storage order + BOOST_CHECK (roboptim::StorageOrder == Eigen::ROBOPTIM_STORAGE_ORDER); + Null null; NoTitle notitle; diff --git a/tests/operator-scalar.cc b/tests/operator-scalar.cc index 74cb87499..e9c2b5546 100644 --- a/tests/operator-scalar.cc +++ b/tests/operator-scalar.cc @@ -25,11 +25,6 @@ #include #include -#include -#include -#include - -#include #include using namespace roboptim; @@ -43,15 +38,12 @@ BOOST_FIXTURE_TEST_SUITE (core, TestSuiteConfiguration) BOOST_AUTO_TEST_CASE_TEMPLATE (scalar_test, T, functionTypes_t) { typename GenericIdentityFunction::result_t offset (5); - offset.setZero (); + offset << 1., 2., 0., 4., 5.; boost::shared_ptr > identity = boost::make_shared > (offset); - boost::shared_ptr > constant = - boost::make_shared > (offset); - boost::shared_ptr > - fct = 2. * identity * constant + constant - constant; + boost::shared_ptr > fct = 2. * identity; typename GenericIdentityFunction::vector_t x (5); x.setZero (); diff --git a/tests/visualization-gnuplot-differentiable-function.stdout b/tests/visualization-gnuplot-differentiable-function.stdout index 02dbbc514..09d453d70 100644 --- a/tests/visualization-gnuplot-differentiable-function.stdout +++ b/tests/visualization-gnuplot-differentiable-function.stdout @@ -30,13 +30,13 @@ set yrange [0:7] reverse set size ratio -1 unset colorbox plot '-' using ($1+0.5):($2+0.5):($3 == 0 ? 0 : 1) matrix with image notitle -1 0 0 0 1 1 0 -1 0 1 0 0 0 1 -1 0 1 0 0 0 1 -1 1 1 0 0 1 0 -0 0 1 0 1 0 0 -0 0 1 0 1 0 0 -0 0 1 0 1 1 1 +1.00000000 0.00000000 0.00000000 0.00000000 1.00000000 1.00000000 0.00000000 +1.00000000 0.00000000 1.00000000 0.00000000 0.00000000 0.00000000 1.00000000 +1.00000000 0.00000000 1.00000000 0.00000000 0.00000000 0.00000000 1.00000000 +1.00000000 1.00000000 1.00000000 0.00000000 0.00000000 1.00000000 0.00000000 +0.00000000 0.00000000 1.00000000 0.00000000 1.00000000 0.00000000 0.00000000 +0.00000000 0.00000000 1.00000000 0.00000000 1.00000000 0.00000000 0.00000000 +0.00000000 0.00000000 1.00000000 0.00000000 1.00000000 1.00000000 1.00000000 e e From d74356f364ddea8984e00d9ee282234a595ac80c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Benjamin=20Chr=C3=A9tien?= Date: Wed, 8 Apr 2015 11:22:42 +0900 Subject: [PATCH 2/3] Fix some types. --- .../core/decorator/finite-difference-gradient.hh | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/include/roboptim/core/decorator/finite-difference-gradient.hh b/include/roboptim/core/decorator/finite-difference-gradient.hh index 74442bf5b..bf493cbbb 100644 --- a/include/roboptim/core/decorator/finite-difference-gradient.hh +++ b/include/roboptim/core/decorator/finite-difference-gradient.hh @@ -55,9 +55,8 @@ namespace roboptim /// \return output stream virtual std::ostream& print (std::ostream& o) const; - /// \brief Gradient has been computed for this point. - vector_t x_; + argument_t x_; /// \brief Analytical gradient. gradient_t analyticalGradient_; @@ -106,9 +105,8 @@ namespace roboptim /// \return output stream virtual std::ostream& print (std::ostream& o) const; - /// \brief Jacobian has been computed for this point. - vector_t x_; + argument_t x_; /// \brief Analytical Jacobian. gradient_t analyticalJacobian_; @@ -183,10 +181,10 @@ namespace roboptim /// \brief Wrapped function. const GenericFunction& adaptee_; - /// \brief Vector storing temporary Jacobian row. - mutable gradient_t column_; - /// \brief Vector storing temporary Jacobian column. + mutable vector_t column_; + + /// \brief Vector storing temporary Jacobian row. mutable gradient_t gradient_; }; From a4dcac74c0b4d082b54572e122b7528ca56022ce Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Benjamin=20Chr=C3=A9tien?= Date: Wed, 8 Apr 2015 11:23:35 +0900 Subject: [PATCH 3/3] Add new derivative_t type. gradient_t was previously used for this, but gradient_t describes row vectors while derivative_t describes column vectors. --- include/roboptim/core/function.hh | 7 ++++ .../core/n-times-derivable-function.hh | 36 ++++++++++++++----- tests/n-times-derivable-function.cc | 4 +-- 3 files changed, 37 insertions(+), 10 deletions(-) diff --git a/include/roboptim/core/function.hh b/include/roboptim/core/function.hh index 59750e1bc..ab3929d12 100644 --- a/include/roboptim/core/function.hh +++ b/include/roboptim/core/function.hh @@ -79,6 +79,11 @@ typedef typename GenericFunctionTraits::NAME##_ref NAME##_ref; \ typedef typename GenericFunctionTraits::const_##NAME##_ref const_##NAME##_ref +# define ROBOPTIM_GENERATE_TRAITS_REFS_T(NAME,TRAITS) \ + typedef typename GenericFunctionTraits::NAME##_t NAME##_t; \ + typedef typename GenericFunctionTraits::NAME##_ref NAME##_ref; \ + typedef typename GenericFunctionTraits::const_##NAME##_ref const_##NAME##_ref + # define ROBOPTIM_FUNCTION_FWD_TYPEDEFS(PARENT) \ typedef PARENT parent_t; \ typedef parent_t::value_type value_type; \ @@ -621,6 +626,7 @@ namespace roboptim ROBOPTIM_GENERATE_TYPEDEFS_EIGEN_REF_VEC(gradient,rowVector_t); ROBOPTIM_GENERATE_TYPEDEFS_EIGEN_REF(jacobian,matrix_t); ROBOPTIM_GENERATE_TYPEDEFS_EIGEN_REF(hessian,matrix_t); + ROBOPTIM_GENERATE_TYPEDEFS_EIGEN_REF(derivative,vector_t); }; /// \brief Trait specializing GenericFunction for Eigen sparse matrices. @@ -665,6 +671,7 @@ namespace roboptim BOOST_PP_COMMA() Eigen::RowMajor>); ROBOPTIM_GENERATE_TYPEDEFS_REF(jacobian,matrix_t); ROBOPTIM_GENERATE_TYPEDEFS_REF(hessian,matrix_t); + ROBOPTIM_GENERATE_TYPEDEFS_EIGEN_REF(derivative,vector_t); }; /// @} diff --git a/include/roboptim/core/n-times-derivable-function.hh b/include/roboptim/core/n-times-derivable-function.hh index 2de9974d2..384cf1e41 100644 --- a/include/roboptim/core/n-times-derivable-function.hh +++ b/include/roboptim/core/n-times-derivable-function.hh @@ -25,6 +25,14 @@ # include # include +# define ROBOPTIM_NTIMES_DERIVABLE_FUNCTION_FWD_TYPEDEFS(PARENT) \ + ROBOPTIM_TWICE_DIFFERENTIABLE_FUNCTION_FWD_TYPEDEFS(PARENT); \ + ROBOPTIM_GENERATE_FWD_REFS (derivative) + +# define ROBOPTIM_NTIMES_DERIVABLE_FUNCTION_FWD_TYPEDEFS_(PARENT) \ + ROBOPTIM_TWICE_DIFFERENTIABLE_FUNCTION_FWD_TYPEDEFS_(PARENT); \ + ROBOPTIM_GENERATE_FWD_REFS_ (derivative) + namespace roboptim { /// \addtogroup roboptim_function @@ -42,6 +50,18 @@ namespace roboptim class NTimesDerivableFunction<2> : public TwiceDifferentiableFunction { public: + + /// \brief Parent type. + typedef TwiceDifferentiableFunction parent_t; + + /// \brief Traits type. + typedef parent_t::traits_t traits_t; + + /// \brief Derivative type. + /// + /// Derivatives are column vectors. + ROBOPTIM_GENERATE_TRAITS_REFS_T(derivative,traits_t); + /// \brief Function derivability order. One static const variable per class /// in inheritance structure. static const size_type derivabilityOrder = 2; @@ -65,7 +85,7 @@ namespace roboptim /// \brief Check if a derivative is valid (check sizes). /// \param derivative derivative vector to be checked /// \return true if valid, false if not - bool isValidDerivative (const_gradient_ref derivative) const + bool isValidDerivative (const_derivative_ref derivative) const { return derivative.size () == this->derivativeSize (); } @@ -106,9 +126,9 @@ namespace roboptim /// \param argument point at which the derivative will be computed /// \param order derivative order (if 0 then function is evaluated) /// \return derivative vector - gradient_t derivative (value_type argument, size_type order = 1) const + derivative_t derivative (value_type argument, size_type order = 1) const { - gradient_t derivative (derivativeSize ()); + derivative_t derivative (derivativeSize ()); derivative.setZero (); this->derivative (derivative, argument, order); return derivative; @@ -120,7 +140,7 @@ namespace roboptim /// \param derivative derivative will be stored in this vector /// \param argument point at which the derivative will be computed /// \param order derivative order (if 0 then function is evaluated) - void derivative (gradient_ref derivative, + void derivative (derivative_ref derivative, value_type argument, size_type order = 1) const { @@ -184,7 +204,7 @@ namespace roboptim /// The gradient is computed for a specific sub-function which id /// is passed through the functionId argument. /// \warning Do not call this function directly, call #gradient - //// or #derivative instead. + /// or #derivative instead. /// \param gradient gradient will be store in this argument /// \param argument point where the gradient will be computed /// \param functionId evaluated function id in the split representation @@ -196,7 +216,7 @@ namespace roboptim Eigen::internal::set_is_malloc_allowed (true); #endif //! ROBOPTIM_DO_NOT_CHECK_ALLOCATION - gradient_t derivative (derivativeSize ()); + derivative_t derivative (derivativeSize ()); derivative.setZero (); this->derivative (derivative, argument[0], 1); @@ -211,7 +231,7 @@ namespace roboptim /// \param derivative derivative will be store in this argument /// \param argument point where the gradient will be computed /// \param order derivative order (if 0 evaluates the function) - virtual void impl_derivative (gradient_ref derivative, + virtual void impl_derivative (derivative_ref derivative, value_type argument, size_type order = 1) const = 0; @@ -235,7 +255,7 @@ namespace roboptim assert (functionId == 0); - gradient_t derivative (derivativeSize ()); + derivative_t derivative (derivativeSize ()); derivative.setZero (); this->derivative (derivative, argument[0], 2); diff --git a/tests/n-times-derivable-function.cc b/tests/n-times-derivable-function.cc index e5a187c46..c014f35e2 100644 --- a/tests/n-times-derivable-function.cc +++ b/tests/n-times-derivable-function.cc @@ -36,7 +36,7 @@ struct F2 : public NTimesDerivableFunction<2> result.setZero (); } - virtual void impl_derivative (gradient_ref derivative, + virtual void impl_derivative (derivative_ref derivative, double, size_type ROBOPTIM_DEBUG_ONLY (order = 1)) const { @@ -58,7 +58,7 @@ struct F10 : public NTimesDerivableFunction<10> result.setZero (); } - virtual void impl_derivative (gradient_ref derivative, + virtual void impl_derivative (derivative_ref derivative, double, size_type ROBOPTIM_DEBUG_ONLY (order = 1)) const {