Merge pull request #240 from stan-dev/feature/issue-202-vectorize-all

Fixes #202. Feature/issue 202 vectorize all

stan/math/fwd/mat/vectorize/apply_scalar_unary.hpp

@@ -0,0 +1,42 @@
+#ifndef STAN_MATH_FWD_MAT_VECTORIZE_APPLY_UNARY_SCALAR_HPP
+#define STAN_MATH_FWD_MAT_VECTORIZE_APPLY_UNARY_SCALAR_HPP
+
+#include <stan/math/prim/mat/vectorize/apply_scalar_unary.hpp>
+#include <stan/math/fwd/core/fvar.hpp>
+
+namespace stan {
+
+  namespace math {
+
+    /**
+     * Template specialization to fvar for vectorizing a unary scalar
+     * function.  This is a base scalar specialization.  It applies
+     * the function specified by the template parameter to the
+     * argument.  
+     *
+     * @tparam F Type of function to apply.
+     * @tparam T Value and tangent type for for forward-mode
+     * autodiff variable.
+     */
+    template <typename F, typename T>
+    struct apply_scalar_unary<F, stan::math::fvar<T> > {
+      /**
+       * Function return type, which is same as the argument type for
+       * the function, <code>fvar&lt;T&gt;</code>.
+       */
+      typedef stan::math::fvar<T> return_t;
+
+      /**
+       * Apply the function specified by F to the specified argument.
+       *
+       * @param x Argument variable.
+       * @return Function applied to the variable.
+       */
+      static inline return_t apply(const stan::math::fvar<T>& x) {
+        return F::fun(x);
+      }
+    };
+
+  }
+}
+#endif

stan/math/prim/mat/vectorize/apply_scalar_unary.hpp

@@ -0,0 +1,160 @@
+#ifndef STAN_MATH_PRIM_MAT_VECTORIZE_APPLY_UNARY_SCALAR_HPP
+#define STAN_MATH_PRIM_MAT_VECTORIZE_APPLY_UNARY_SCALAR_HPP
+
+#include <Eigen/Dense>
+#include <vector>
+
+namespace stan {
+
+  namespace math {
+
+    /**
+     * Base template class for vectorization of unary scalar functions
+     * defined by a template class <code>F</code> to a scalar,
+     * standard library vector, or Eigen dense matrix expression
+     * template.  
+     *
+     * <p>The base class applies to any Eigen dense matrix expression
+     * template.  Specializations define applications to scalars
+     * (primitive or autodiff in the corresponding autodiff library
+     * directories) or to standard library vectors of vectorizable
+     * types (primitives, Eigen dense matrix expressions, or further
+     * standard vectors).
+     *
+     * <p>Each specialization must define the typedef
+     * <code>return_t</code> for the vectorized return type and the
+     * function <code>apply</code> which defines the vectorization or
+     * base application of the function defined statically by the
+     * class F.  The function definition class F defines a static
+     * function <code>fun()</code>, which defines the function's
+     * behavior on scalars.
+     *
+     * @tparam F Type of function to apply.
+     * @tparam T Type of argument to which function is applied.
+     */
+    template <typename F, typename T>
+    struct apply_scalar_unary {
+      /**
+       * Type of underlying scalar for the matrix type T.
+       */
+      typedef typename Eigen::internal::traits<T>::Scalar scalar_t;
+
+      /**
+       * Return type for applying the function elementwise to a matrix
+       * expression template of type T.
+       */
+      typedef Eigen::Matrix<scalar_t, T::RowsAtCompileTime,
+                            T::ColsAtCompileTime>
+      return_t;
+
+      /**
+       * Return the result of applying the function defined by the
+       * template parameter F to the specified matrix argument. 
+       *
+       * @param x Matrix to which operation is applied.
+       * @return Componentwise application of the function specified
+       * by F to the specified matrix.
+       */
+      static inline return_t apply(const T& x) {
+        return_t result(x.rows(), x.cols());
+        for (int j = 0; j < x.cols(); ++j)
+          for (int i = 0; i < x.rows(); ++i)
+            result(i, j) = apply_scalar_unary<F, scalar_t>::apply(x(i, j));
+        return result;
+      }
+    };
+
+    /**
+     * Template specialization for vectorized functions applying to
+     * double arguments.
+     *
+     * @tparam F Type of function defining static apply function.
+     */
+    template <typename F>
+    struct apply_scalar_unary<F, double> {
+      /**
+       * The return type, double.
+       */
+      typedef double return_t;
+
+      /**
+       * Apply the function specified by F to the specified argument.
+       * This is defined through a direct application of
+       * <code>F::fun()</code>, which must be defined for double
+       * arguments. 
+       *
+       * @param x Argument scalar.
+       * @return Result of applying F to the scalar.
+       */
+      static inline return_t apply(double x) {
+        return F::fun(x);
+      }
+    };
+
+    /**
+     * Template specialization for vectorized functions applying to
+     * integer arguments.  Although the argument is integer, the
+     * return type is specified as double.  This allows promotion of
+     * integers to doubles in vectorized functions, or in containers.  
+     *
+     * @tparam F Type of function defining static apply function.
+     */
+    template <typename F>
+    struct apply_scalar_unary<F, int> {
+      /**
+       * The return type, double.
+       */
+      typedef double return_t;
+
+      /**
+       * Apply the function specified by F to the specified argument.
+       * This is defined through a direct application of
+       * <code>F::fun()</code>, which must be defined for double
+       * arguments. 
+       *
+       * @param x Argument scalar.
+       * @return Result of applying F to the scalar.
+       */
+      static inline return_t apply(int x) {
+        return F::fun(static_cast<double>(x));
+      }
+    };
+
+    /**
+     * Template specialization for vectorized functions applying to
+     * standard vector containers.  The lowest-level scalar type of
+     * the argument will determine the return type.  Integers are
+     * promoted to double values.
+     *
+     * @tparam F Class defining a static apply function.
+     * @tparam T Type of element contained in standard vector.
+     */
+    template <typename F, typename T>
+    struct apply_scalar_unary<F, std::vector<T> > {
+      /**
+       * Return type, which is calculated recursively as a standard
+       * vector of the return type of the contained type T.
+       */
+      typedef typename std::vector<typename apply_scalar_unary<F, T>::return_t>
+      return_t;
+
+      /**
+       * Apply the function specified by F elementwise to the
+       * specified argument.  This is defined recursively through this
+       * class applied to elements of type T.
+       *
+       * @param x Argument container.
+       * @return Elementwise application of F to the elements of the
+       * container.
+       */
+      static inline return_t apply(const std::vector<T>& x) {
+        return_t fx(x.size());
+        for (size_t i = 0; i < x.size(); ++i)
+          fx[i] = apply_scalar_unary<F, T>::apply(x[i]);
+        return fx;
+      }
+    };
+
+  }
+}
+#endif

stan/math/rev/mat/vectorize/apply_scalar_unary.hpp

@@ -0,0 +1,39 @@
+#ifndef STAN_MATH_REV_MAT_VECTORIZE_APPLY_UNARY_SCALAR_HPP
+#define STAN_MATH_REV_MAT_VECTORIZE_APPLY_UNARY_SCALAR_HPP
+
+#include <stan/math/prim/mat/vectorize/apply_scalar_unary.hpp>
+#include <stan/math/rev/core/var.hpp>
+
+namespace stan {
+
+  namespace math {
+
+    /**
+     * Template specialization to var for vectorizing a unary scalar
+     * function.  This is a base scalar specialization.  It applies
+     * the function specified by the template parameter to the
+     * argument.  
+     *
+     * @tparam F Type of function to apply.
+     */
+    template <typename F>
+    struct apply_scalar_unary<F, stan::math::var> {
+      /**
+       * Function return type, which is <code>var</code>.
+       */
+      typedef stan::math::var return_t;
+
+      /**
+       * Apply the function specified by F to the specified argument.  
+       * 
+       * @param x Argument variable.
+       * @return Function applied to the variable.
+       */
+      static inline return_t apply(const stan::math::var& x) {
+        return F::fun(x);
+      }
+    };
+
+  }
+}
+#endif

test/unit/math/fwd/mat/vectorize/build_fwd_matrix.hpp

@@ -0,0 +1,30 @@
+#ifndef TEST_UNIT_MATH_FWD_MAT_VECTORIZE_BUILD_FWD_MATRIX_HPP
+#define TEST_UNIT_MATH_FWD_MAT_VECTORIZE_BUILD_FWD_MATRIX_HPP
+
+#include <stan/math/fwd/mat.hpp>
+#include <Eigen/Dense>
+#include <vector>
+#include <test/unit/math/fwd/mat/vectorize/build_fwd_vector.hpp>
+
+template <typename F, typename T, int R, int C>
+static inline Eigen::Matrix<T, R, C>
+build_fwd_matrix(const Eigen::Matrix<T, R, C>& x, int seed_index = -1) {
+  using Eigen::Matrix;
+  using std::vector;
+
+  Eigen::Matrix<T, R, C> fvar_matrix(x.rows(), x.cols());
+  size_t num_inputs = F::valid_inputs().size();
+  //Fills matrix with copies of valid_input values
+  for (int i = 0; i < x.size(); ++i) {
+    std::vector<T> inputs;
+    if (seed_index == i)
+      inputs
+      = build_fwd_vector<F>(std::vector<T>(), seed_index % num_inputs);
+    else
+      inputs = build_fwd_vector<F>(std::vector<T>()); 
+    fvar_matrix(i) = inputs[i % num_inputs];
+  }
+  return fvar_matrix;
+}
+
+#endif

test/unit/math/fwd/mat/vectorize/build_fwd_vector.hpp

@@ -0,0 +1,35 @@
+#ifndef TEST_UNIT_MATH_FWD_MAT_VECTORIZE_BUILD_FWD_VECTOR_HPP
+#define TEST_UNIT_MATH_FWD_MAT_VECTORIZE_BUILD_FWD_VECTOR_HPP
+
+#include <stan/math/fwd/mat.hpp>
+#include <vector>
+
+template <typename F>
+static inline std::vector<double>
+build_fwd_vector(std::vector<double> double_vector,
+                          int seed_index = -1) { 
+  return F::valid_inputs();
+}
+
+template <typename F, typename T>
+static inline std::vector<stan::math::fvar<T> >
+build_fwd_vector(std::vector<stan::math::fvar<T> > fvar_vector,
+                        int seed_index = -1) { 
+  using std::vector;
+  using stan::math::fvar;
+
+  vector<T> template_v = build_fwd_vector<F>(vector<T>(), seed_index);
+
+  for (size_t i = 0; i < template_v.size(); ++i) {
+
+    // For fvar<fvar<double> >, this will fill in 
+    // all four components
+    if (seed_index == static_cast<int>(i))
+      fvar_vector.push_back(fvar<T>(template_v[i], template_v[i]));
+    else
+      fvar_vector.push_back(fvar<T>(template_v[i]));
+  }
+  return fvar_vector;
+}
+
+#endif

test/unit/math/fwd/mat/vectorize/expect_fwd_errors.hpp

@@ -0,0 +1,26 @@
+#ifndef TEST_UNIT_MATH_FWD_MAT_VECTORIZE_EXPECT_FWD_ERRORS_HPP
+#define TEST_UNIT_MATH_FWD_MAT_VECTORIZE_EXPECT_FWD_ERRORS_HPP
+
+#include <stan/math/fwd/mat.hpp>
+#include <test/unit/math/prim/mat/vectorize/expect_scalar_error.hpp>
+#include <test/unit/math/prim/mat/vectorize/expect_std_vector_error.hpp>
+#include <test/unit/math/prim/mat/vectorize/expect_matrix_error.hpp>
+#include <test/unit/math/prim/mat/vectorize/expect_vector_error.hpp>
+#include <test/unit/math/prim/mat/vectorize/expect_row_vector_error.hpp>
+
+template <typename F>
+void expect_fwd_errors() {
+  using stan::math::fvar;
+  expect_scalar_error<F, fvar<double> >();
+  expect_scalar_error<F, fvar<fvar<double> > >();
+  expect_std_vector_error<F, fvar<double> >();
+  expect_std_vector_error<F, fvar<fvar<double> > >();
+  expect_matrix_error<F, fvar<double> >();
+  expect_matrix_error<F, fvar<fvar<double> > >();
+  expect_vector_error<F, fvar<double> >();
+  expect_vector_error<F, fvar<fvar<double> > >();
+  expect_row_vector_error<F, fvar<double> >();
+  expect_row_vector_error<F, fvar<fvar<double> > >();
+}
+
+#endif

test/unit/math/fwd/mat/vectorize/expect_fwd_matrix_value.hpp

@@ -0,0 +1,49 @@
+#ifndef TEST_UNIT_MATH_FWD_MAT_VECTORIZE_EXPECT_FWD_MATRIX_VALUE_HPP
+#define TEST_UNIT_MATH_FWD_MAT_VECTORIZE_EXPECT_FWD_MATRIX_VALUE_HPP
+
+#include <stan/math/fwd/mat.hpP>
+#include <vector>
+#include <Eigen/Dense>
+#include <test/unit/math/fwd/mat/vectorize/build_fwd_matrix.hpp>
+#include <test/unit/math/fwd/mat/vectorize/expect_val_deriv_eq.hpp>
+
+template <typename F, typename T>
+void expect_fwd_matrix_value() {
+  using stan::math::fvar;
+  using std::vector;
+  typedef Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> matrix_t;
+
+  int num_inputs = F::valid_inputs().size();
+  int num_cols = 3;
+  matrix_t template_m(num_inputs, num_cols);
+
+  for (int i = 0; i < template_m.size(); ++i) {
+    matrix_t a = build_fwd_matrix<F>(template_m, i);
+    matrix_t fa = F::template apply<matrix_t>(a);
+    EXPECT_EQ(a.size(), fa.size());
+    expect_val_deriv_eq(F::apply_base(a(i)), fa(i));
+  }
+
+  size_t vector_matrix_size = 2;
+  for (size_t i = 0; i < vector_matrix_size; ++i) {
+    for (int j = 0; j < template_m.size(); ++j) {
+      vector<matrix_t> b;
+      for (size_t k = 0; k < vector_matrix_size; ++k)
+        if (k == i)
+          b.push_back(build_fwd_matrix<F>(template_m, j));
+        else
+          b.push_back(build_fwd_matrix<F>(template_m));
+      vector<matrix_t> fb = F::template apply<vector<matrix_t> >(b);
+      EXPECT_EQ(b.size(), fb.size());
+      EXPECT_EQ(b[i].size(), fb[i].size());
+      expect_val_deriv_eq(F::apply_base(b[i](j)), fb[i](j));
+    }
+  }
+
+  int seed_i = num_inputs + 1;
+  matrix_t a = build_fwd_matrix<F>(template_m, seed_i);
+  matrix_t fab = F::template apply<matrix_t>(a.block(1, 1, 1, 1));
+  expect_val_deriv_eq(F::apply_base(a(1,1)), fab(0,0));
+}
+
+#endif