Fix compilation for both ColMajor/RowMajor + provide a define.

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)

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<T>& 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_;
     };
 

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 <int SO>
+    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.).

include/roboptim/core/function.hh

@@ -27,6 +27,7 @@
 # include <utility>
 # include <vector>
 
+# include <boost/static_assert.hpp>
 # include <boost/algorithm/string/replace.hpp>
 # include <boost/tuple/tuple.hpp>
 # include <boost/preprocessor/punctuation/comma.hpp>
@@ -43,11 +44,21 @@
 # include <roboptim/core/indent.hh>
 # include <roboptim/core/portability.hh>
 
+# include <roboptim/core/detail/utility.hh>
+
+// 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##_t> NAME##_ref;			\
   typedef const Eigen::Ref<const NAME##_t>& 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<NAME##_t, 0, detail::row_vector_stride<StorageOrder>::type> NAME##_ref;			\
+  typedef const Eigen::Ref<const NAME##_t, 0, detail::row_vector_stride<StorageOrder>::type>& const_##NAME##_ref
+
 # define ROBOPTIM_GENERATE_TYPEDEFS_REF(NAME,TYPE)	\
   typedef TYPE NAME##_t;				\
   typedef NAME##_t& NAME##_ref;				\
@@ -68,6 +79,11 @@
   typedef typename GenericFunctionTraits<T>::NAME##_ref NAME##_ref;	\
   typedef typename GenericFunctionTraits<T>::const_##NAME##_ref const_##NAME##_ref
 
+# define ROBOPTIM_GENERATE_TRAITS_REFS_T(NAME,TRAITS)				\
+  typedef typename GenericFunctionTraits<TRAITS>::NAME##_t NAME##_t;		\
+  typedef typename GenericFunctionTraits<TRAITS>::NAME##_ref NAME##_ref;	\
+  typedef typename GenericFunctionTraits<TRAITS>::const_##NAME##_ref const_##NAME##_ref
+
 # define ROBOPTIM_FUNCTION_FWD_TYPEDEFS(PARENT)	\
   typedef PARENT parent_t;			\
   typedef parent_t::value_type value_type;	\
@@ -76,6 +92,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 +103,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 +171,7 @@ namespace roboptim
     /// used for computations.
     typedef typename GenericFunctionTraits<T>::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 +181,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,63 +587,91 @@ namespace roboptim
   struct GenericFunctionTraits<EigenMatrixDense>
   {
     /// \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<double BOOST_PP_COMMA()
+     Eigen::Matrix<value_type BOOST_PP_COMMA()
      Eigen::Dynamic BOOST_PP_COMMA()
      Eigen::Dynamic BOOST_PP_COMMA()
      StorageOrder>);
+
+    // Vector types
     ROBOPTIM_GENERATE_TYPEDEFS_EIGEN_REF
     (vector,
-     Eigen::Matrix<double BOOST_PP_COMMA()
+     Eigen::Matrix<value_type BOOST_PP_COMMA()
      Eigen::Dynamic BOOST_PP_COMMA()
      1>);
 
+    // Row vector types
+    ROBOPTIM_GENERATE_TYPEDEFS_EIGEN_REF
+    (rowVector,
+     Eigen::Matrix<value_type BOOST_PP_COMMA()
+     1 BOOST_PP_COMMA()
+     Eigen::Dynamic>);
+
     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);
+    ROBOPTIM_GENERATE_TYPEDEFS_EIGEN_REF(derivative,vector_t);
   };
 
   /// \brief Trait specializing GenericFunction for Eigen sparse matrices.
   template <>
   struct GenericFunctionTraits<EigenMatrixSparse>
   {
     /// \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<double BOOST_PP_COMMA() StorageOrder>);
+     Eigen::SparseMatrix<value_type BOOST_PP_COMMA() StorageOrder>);
 
+    // Vector types
     ROBOPTIM_GENERATE_TYPEDEFS_EIGEN_REF
     (vector,
-     Eigen::Matrix<double BOOST_PP_COMMA()
+     Eigen::Matrix<value_type BOOST_PP_COMMA()
      Eigen::Dynamic BOOST_PP_COMMA()
      1>);
 
+    // Row vector types
+    ROBOPTIM_GENERATE_TYPEDEFS_EIGEN_REF
+    (rowVector,
+     Eigen::Matrix<value_type BOOST_PP_COMMA()
+     1 BOOST_PP_COMMA()
+     Eigen::Dynamic>);
+
     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<double>);
+    ROBOPTIM_GENERATE_TYPEDEFS_REF(gradient,Eigen::SparseVector<value_type \
+                                   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);
   };
 
   /// @}

include/roboptim/core/n-times-derivable-function.hh

@@ -25,6 +25,14 @@
 # include <roboptim/core/fwd.hh>
 # include <roboptim/core/twice-differentiable-function.hh>
 
+# 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);