This question is a spin-off of another question I had concerning boost::fusion. The idea is to use boost::fusion to iterate over a large C-style struct containing N-dimension arrays. Computation on these arrays is done by Eigen. By using boost::fusion, one can apply simple arithmetic operations on the whole C struct, e.g. scalar multiplication or vector addition.
When dealing with binary operations, I use boost::fusion::zip to form a single sequence, and boost::fusion::for_each to iterate other that sequence.
The problem with boost::fusion::zip is that it constructs const sequences while I need to modify one of the values (e.g. return value of the addition). Thus, I end up using const_cast to modify that value (Eigen vector), but for some reason I cannot use a result_ref in the add() function. Why is that?
Moreover, is there any better (or simpler) way to achieve what I am trying to do? boost::fusion::zip may not be the best fit, but I could not find any other simple way to do that.
#include <iostream>
#include <boost/fusion/adapted/struct/adapt_struct.hpp>
#include <boost/fusion/include/adapt_struct.hpp>
#include <boost/fusion/algorithm/iteration/for_each.hpp>
#include <boost/fusion/include/for_each.hpp>
#include <boost/fusion/algorithm/transformation/zip.hpp>
#include <boost/fusion/include/zip.hpp>
#include <boost/bind.hpp>
#include <boost/fusion/container/vector/vector30.hpp>
#include <boost/fusion/include/vector30.hpp>
#include <boost/fusion/sequence/intrinsic/at_c.hpp>
#include <boost/fusion/include/at_c.hpp>
#include <boost/type_traits/remove_const.hpp>
#include <boost/type_traits/remove_reference.hpp>
#include <Eigen/Core>
template <class type_const_ref>
struct remove_const_ref
{
typedef typename boost::remove_reference <type_const_ref>::type type_const;
typedef typename boost::remove_const <type_const_ref>::type type_ref;
typedef typename boost::remove_const <type_const >::type type;
};
namespace demo
{
template<typename T, int SIZE1, int SIZE2>
struct data
{
T ar1[SIZE1][SIZE2];
T ar2[SIZE1][SIZE2];
};
template<typename T>
struct EigenMap
{
typedef Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, 1> > type;
};
template<typename T>
struct data_eigen
{
template <int SIZE1, int SIZE2>
data_eigen(data<T,SIZE1,SIZE2>& src)
: ar1(typename EigenMap<T>::type(&src.ar1[0][0], SIZE1*SIZE2)),
ar2(typename EigenMap<T>::type(&src.ar2[0][0], SIZE1*SIZE2))
{
}
typename EigenMap<T>::type ar1;
typename EigenMap<T>::type ar2;
};
struct print
{
template<typename T>
void operator()(const Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, 1> >& t) const
{
std::cout << t.transpose() << std::endl;
}
};
struct scalarMult
{
template<typename T, typename U>
void operator()(T& t, U& u) const
{
t *= u;
}
};
template <typename T>
struct add
{
template <typename U>
void operator() (const boost::fusion::vector3<U,U,U>& t) const
{
typedef typename remove_const_ref<U>::type_ref vector_ref;
typedef typename remove_const_ref<U>::type vector_type;
// FIXME: find why we cannot use vector_ref
vector_type result_ref = const_cast<vector_ref>(boost::fusion::at_c<2>(t));
result_ref = boost::fusion::at_c<0>(t) + boost::fusion::at_c<1>(t);
}
};
}
BOOST_FUSION_ADAPT_TPL_STRUCT
(
(T),
(demo::data_eigen) (T),
(typename demo::EigenMap<T>::type, ar1)
(typename demo::EigenMap<T>::type, ar2)
)
int main()
{
typedef float REALTYPE;
const int SIZE1 = 2;
const int SIZE2 = 2;
// Basic data structure with multidimensional arrays
demo::data<REALTYPE, SIZE1, SIZE2> d1;
for (unsigned int i = 0; i < SIZE1; ++i)
for (unsigned int j = 0; j < SIZE2; ++j)
{
d1.ar1[i][j] = (i+1)*(j+1);
d1.ar2[i][j] = i + j;
}
demo::data<REALTYPE, SIZE1, SIZE2> d2;
demo::data<REALTYPE, SIZE1, SIZE2> d3;
memset(&d3, 0, sizeof(demo::data<REALTYPE, SIZE1, SIZE2>));
for (unsigned int i = 0; i < SIZE1; ++i)
for (unsigned int j = 0; j < SIZE2; ++j)
{
d2.ar1[i][j] = 1.0;
d2.ar2[i][j] = 2.0;
}
// Eigen::Map + BOOST_FUSION_ADAPT_TPL_STRUCT
demo::data_eigen<REALTYPE> eig_d1(d1);
demo::data_eigen<REALTYPE> eig_d2(d2);
demo::data_eigen<REALTYPE> eig_d3(d3);
std::cout << "d1:" << std::endl;
boost::fusion::for_each(eig_d1, demo::print());
std::cout << std::endl;
std::cout << "d2:" << std::endl;
boost::fusion::for_each(eig_d2, demo::print());
std::cout << std::endl;
boost::fusion::for_each(eig_d1, boost::bind<void>(demo::scalarMult(), _1, 2.0));
std::cout << "d1 = 2 * d1:" << std::endl;
boost::fusion::for_each(eig_d1, demo::print());
std::cout << std::endl;
boost::fusion::for_each(boost::fusion::zip(eig_d1, eig_d2, eig_d3),
demo::add<REALTYPE>());
std::cout << "d3 = d1 + d2:" << std::endl;
boost::fusion::for_each(eig_d3, demo::print());
std::cout << std::endl;
return EXIT_SUCCESS;
}
I believe what you need to use is zip_view.
Your for_each invocation would be:
typedef demo::data_eigen<REALTYPE>& vector_ref;
typedef boost::fusion::vector<vector_ref,vector_ref,vector_ref> my_zip;
boost::fusion::for_each(boost::fusion::zip_view<my_zip>(my_zip(eig_d1, eig_d2, eig_d3)), demo::add());
and your add functor would simply be:
struct add
{
template <typename ZipView>
void operator() (const ZipView& t) const //CHANGED
{
boost::fusion::at_c<2>(t) = boost::fusion::at_c<0>(t) + boost::fusion::at_c<1>(t);
}
};
The following code has been tested with g++ 4.8.0 and outputs the same result as yours:
#include <iostream>
#include <boost/bind.hpp>
#include <boost/fusion/include/at_c.hpp>
#include <boost/fusion/include/adapt_struct.hpp>
#include <boost/fusion/include/for_each.hpp>
#include <boost/fusion/include/zip_view.hpp>
#include <boost/fusion/include/vector.hpp>
#include <Eigen/Core>
namespace demo
{
template<typename T, int SIZE1, int SIZE2>
struct data
{
T ar1[SIZE1][SIZE2];
T ar2[SIZE1][SIZE2];
};
template<typename T>
struct EigenMap
{
typedef Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, 1> > type;
};
template<typename T>
struct data_eigen
{
template <int SIZE1, int SIZE2>
data_eigen(data<T,SIZE1,SIZE2>& src)
: ar1(typename EigenMap<T>::type(&src.ar1[0][0], SIZE1*SIZE2)),
ar2(typename EigenMap<T>::type(&src.ar2[0][0], SIZE1*SIZE2))
{
}
typename EigenMap<T>::type ar1;
typename EigenMap<T>::type ar2;
};
struct print
{
template<typename T>
void operator()(const Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, 1> >& t) const
{
std::cout << t.transpose() << std::endl;
}
};
struct scalarMult
{
template<typename T, typename U>
void operator()(T& t, U& u) const
{
t *= u;
}
};
struct add
{
template <typename ZipView>
void operator() (const ZipView& t) const //CHANGED
{
boost::fusion::at_c<2>(t) = boost::fusion::at_c<0>(t) + boost::fusion::at_c<1>(t);
}
};
}
BOOST_FUSION_ADAPT_TPL_STRUCT
(
(T),
(demo::data_eigen) (T),
(typename demo::EigenMap<T>::type, ar1)
(typename demo::EigenMap<T>::type, ar2)
)
int main()
{
typedef float REALTYPE;
const int SIZE1 = 2;
const int SIZE2 = 2;
// Basic data structure with multidimensional arrays
demo::data<REALTYPE, SIZE1, SIZE2> d1;
for (unsigned int i = 0; i < SIZE1; ++i)
for (unsigned int j = 0; j < SIZE2; ++j)
{
d1.ar1[i][j] = (i+1)*(j+1);
d1.ar2[i][j] = i + j;
}
demo::data<REALTYPE, SIZE1, SIZE2> d2;
demo::data<REALTYPE, SIZE1, SIZE2> d3;
memset(&d3, 0, sizeof(demo::data<REALTYPE, SIZE1, SIZE2>));
for (unsigned int i = 0; i < SIZE1; ++i)
for (unsigned int j = 0; j < SIZE2; ++j)
{
d2.ar1[i][j] = 1.0;
d2.ar2[i][j] = 2.0;
}
// Eigen::Map + BOOST_FUSION_ADAPT_TPL_STRUCT
demo::data_eigen<REALTYPE> eig_d1(d1);
demo::data_eigen<REALTYPE> eig_d2(d2);
demo::data_eigen<REALTYPE> eig_d3(d3);
std::cout << "d1:" << std::endl;
boost::fusion::for_each(eig_d1, demo::print());
std::cout << std::endl;
std::cout << "d2:" << std::endl;
boost::fusion::for_each(eig_d2, demo::print());
std::cout << std::endl;
boost::fusion::for_each(eig_d1, boost::bind<void>(demo::scalarMult(), _1, 2.0));
std::cout << "d1 = 2 * d1:" << std::endl;
boost::fusion::for_each(eig_d1, demo::print());
std::cout << std::endl;
typedef demo::data_eigen<REALTYPE>& vector_ref; //ADDITION
typedef boost::fusion::vector<vector_ref,vector_ref,vector_ref> my_zip; //ADDITION
boost::fusion::for_each(boost::fusion::zip_view<my_zip>(my_zip(eig_d1, eig_d2, eig_d3)), //CHANGED
demo::add());
std::cout << "d3 = d1 + d2:" << std::endl;
boost::fusion::for_each(eig_d3, demo::print());
std::cout << std::endl;
return EXIT_SUCCESS;
}