vectorization template metaprograms & problematic inclusion order

[Bob Carpenter]

Here's what I've come up with to deal with template metaprogramming
for vectorization. It all works as far as I can tell, and will
allow vectorization to any depth. But, it doesn't play nicely
with our dependency order (Eigen::Matrix -> std::vector -> scalar),
instead making it (scalar -> Eigen::Matrix, std::vector -> Eigen::matrix).

I really don't like that the base class involves Eigen, but I can't
figure out any other way to get it to work with
all the Eigen matrix types. That is, I can't find an Eigen type
I can drop in as a replacement.

I very much want to be able to write a single templated
version of the functor (example here is exp_fun) and have it
generate all of the instances (double, double[], double[,], ...,
vector, vector[], ..., matrix, matrix[], ...). And have it be
possible to write a more specialized version for var and fvar.

The problem is that the base class for double depends on the base
template, which involves Eigen::Matrix. And I think it reintroduces
order dependency because each class needs to come after the base
template.

Any suggestions on how to fix this would be most appreciated. I'm
stumped right now in terms of how to get it into our directory structure
without violating our dependency rules.

I already wrote a simplified query to the Eigen list, but haven't
gotten a response.

- Bob


template <typename F, typename T>
struct apply_fun {
typedef typename Eigen::internal::traits<T>::Scalar scalar_t;

typedef Eigen::Matrix<scalar_t, T::RowsAtCompileTime, T::ColsAtCompileTime>
return_t;

static inline return_t apply(const T& x) {
return_t result(x.rows(), x.cols());
for (int i = 0; i < x.size(); ++i)
result(i) = apply_fun<F, scalar_t>::apply(x(i));
return result;
}
};

template <typename F>
struct apply_fun<F, double> {
typedef double return_t;

static inline return_t apply(double x) {
return F::fun(x);
}
};

template <typename F>
struct apply_fun<F, stan::math::var> {
typedef stan::math::var return_t;

static inline return_t apply(const stan::math::var& x) {
return F::fun(x);
}
};

template <typename F, typename T>
struct apply_fun<F, stan::math::fvar<T> > {
typedef stan::math::fvar<T> return_t;

static inline return_t apply(const stan::math::fvar<T>& x) {
return F::fun(x);
}
};

template <typename F, typename T>
struct apply_fun<F, std::vector<T> > {
typedef typename std::vector<typename apply_fun<F, T>::return_t> return_t;

static inline return_t apply(const std::vector<T>& x) {
std::vector<T> fx(x.size());
for (size_t i = 0; i < x.size(); ++i)
fx[i] = apply_fun<F, T>::apply(x[i]);
return fx;
}
};



struct exp_fun {
template <typename T>
static inline T fun(const T& x) {
using std::exp;
return exp(x);
}
};


template <typename T>
inline typename apply_fun<exp_fun, T>::return_t
my_exp(const T& x) {
return apply_fun<exp_fun, T>::apply(x);
}

[Michael Betancourt]

Is there no way to do partial template specialization and split
up the specializations for double and arrays/vectors/matrices
into their respective directories?

> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Daniel Lee]

No partial template specification on functions. Only on classes.

Bob, I'm going to think about it tonight and try to split it up reasonably.

Daniel

[Bob Carpenter]

Michael --- yes, they can be split up, but the problem is
that they all depend on the base class template, which
involves Eigen types. So there's no way to split them up
with inducing a dependency from /prim to /mat.

Barring a better solution, I'll just have to drop it all into
/mat and define the templated versions of functions there. It makes
sense in that it's where it would have to be defined anyway because
the top-level function depends on everything. What it disallows
is having a std::vector only version of the vectorization, but given
our main use case in Stan, that's no big deal.

Now, onto the C++ tutorial part of the email. Sorry in advance
if it's obvious---this took me a while to figure out.

You're both right. Yes, Michael, those are partial class template
specializations. Yes, Daniel, you can't do that with functions,
which is why I did it with classes and a simple templated wrapper
function.

That is, instead of a template function foo

template <typename T0, ..., TN>
T0 foo(T1, ..., TN);

which can be called directly as foo<T0,...,TN>(x1,...,xN),
I create a template class

template <typename T0, ..., TN>
struct foo_fun {
static T0 apply(T1, ..., TN);
}

and call foo_fun<T0,...,TN>::apply(x1,...,XN).

The upside is that it is order independent because of
the way class template specializations are resolved. The
downside is that you need to specify all those parameters
in the call and can't do type inference.

The latter problem's taken care of with the generic
wrapper function which passes the types of the arguments
(and function class name type parameter) through.

- Bob