diff --git a/stan/math/fwd/mat/meta/operands_and_partials.hpp b/stan/math/fwd/mat/meta/operands_and_partials.hpp
index ac6b1d2f9eb46e6e93f4118b679d39d17c02a254..695e26f3c7ce0ba5e2165d0281255afc495b5f3a 100644
--- a/stan/math/fwd/mat/meta/operands_and_partials.hpp
+++ b/stan/math/fwd/mat/meta/operands_and_partials.hpp
@@ -9,35 +9,34 @@
namespace stan {
namespace math {
namespace detail {
- template <typename ViewElt, typename Op, typename Dx>
- class ops_partials_edge_vec_fwd
- : public ops_partials_edge_mat_prim<ViewElt, Op, -1, 1> {
+ // Vectorized Univariate
+ template <typename ViewElt, typename Dx>
+ class ops_partials_edge<ViewElt, std::vector<fvar<Dx> > >
+ : public ops_partials_edge_mat_prim<ViewElt, std::vector<fvar<Dx> >, -1, 1> {
public:
- ops_partials_edge_vec_fwd(const Op& ops)
- : ops_partials_edge_mat_prim<ViewElt, Op, -1, 1>(ops) {}
+ explicit ops_partials_edge(const std::vector<fvar<Dx> >& ops)
+ : ops_partials_edge_mat_prim<ViewElt, std::vector<fvar<Dx> >, -1, 1>(ops) {}
Dx dx() {
Dx derivative(0);
for (int i = 0; i < this->size(); ++i) {
- derivative += this->partials[i] * this->operands[i].d_;
+ derivative += this->partials(i) * this->operands[i].d_;
}
return derivative;
}
};
-
- // Vectorized Univariate
- template <typename ViewElt, typename Dx>
- class ops_partials_edge<ViewElt, std::vector<fvar<Dx> > >
- : public ops_partials_edge_vec_fwd<ViewElt, std::vector<fvar<Dx> >, Dx> {
- public:
- explicit ops_partials_edge(const std::vector<fvar<Dx> >& ops)
- : ops_partials_edge_vec_fwd<ViewElt, std::vector<fvar<Dx> >, Dx>(ops) {}
- };
template <typename ViewElt, typename Dx, int R, int C>
class ops_partials_edge<ViewElt, Eigen::Matrix<fvar<Dx>, R, C> >
- : public ops_partials_edge_vec_fwd<ViewElt, Eigen::Matrix<fvar<Dx>, R, C>, Dx> {
+ : public ops_partials_edge_mat_prim<ViewElt, Eigen::Matrix<fvar<Dx>, R, C>, R, C> {
public:
explicit ops_partials_edge(const Eigen::Matrix<fvar<Dx>, R, C>& ops)
- : ops_partials_edge_vec_fwd<ViewElt, Eigen::Matrix<fvar<Dx>, R, C>, Dx>(ops) {}
+ : ops_partials_edge_mat_prim<ViewElt, Eigen::Matrix<fvar<Dx>, R, C>, R, C>(ops) {}
+ Dx dx() {
+ Dx derivative(0);
+ for (int i = 0; i < this->size(); ++i) {
+ derivative += this->partials(i) * this->operands(i).d_;
+ }
+ return derivative;
+ }
};
// Multivariate; vectors of eigen types
diff --git a/stan/math/fwd/scal/meta/OperandsAndPartials.hpp b/stan/math/fwd/scal/meta/OperandsAndPartials.hpp
deleted file mode 100644
index 61a269c0fd3e51e08e26757075cc7dc429c80555..0000000000000000000000000000000000000000
--- a/stan/math/fwd/scal/meta/OperandsAndPartials.hpp
+++ /dev/null
@@ -1,183 +0,0 @@
-#ifndef STAN_MATH_FWD_SCAL_META_OPERANDSANDPARTIALS_HPP
-#define STAN_MATH_FWD_SCAL_META_OPERANDSANDPARTIALS_HPP
-
-#include <stan/math/prim/scal/meta/is_constant_struct.hpp>
-#include <stan/math/prim/scal/meta/length.hpp>
-#include <stan/math/prim/scal/meta/OperandsAndPartials.hpp>
-#include <stan/math/fwd/core.hpp>
-
-namespace stan {
- namespace math {
-
- namespace {
- template <typename T_derivative,
- typename T,
- typename T_partials,
- bool is_vec = is_vector<T>::value,
- bool is_const = is_constant_struct<T>::value>
- struct increment_derivative {
- inline T_derivative operator()(const T& x,
- const T_partials& d_dx) {
- return 0;
- }
- };
-
- template <typename T_derivative,
- typename T,
- typename T_partials>
- struct increment_derivative<T_derivative, T, T_partials, false, false> {
- inline T_derivative operator()(const T& x,
- const T_partials& d_dx) {
- return d_dx[0] * x.d_;
- }
- };
-
- template <typename T_derivative,
- typename T,
- typename T_partials>
- struct increment_derivative<T_derivative, T, T_partials, true, false> {
- inline T_derivative operator()(const T& x,
- const T_partials& d_dx) {
- T_derivative derivative(0);
- for (size_t n = 0; n < length(x); n++)
- derivative += d_dx[n] * x[n].d_;
- return derivative;
- }
- };
-
- template <typename T,
- typename T1, typename D1, typename T2, typename D2,
- typename T3, typename D3, typename T4, typename D4,
- typename T5, typename D5, typename T6, typename D6>
- fvar<T> partials_to_fvar(T& logp,
- const T1& x1, D1& d_x1,
- const T2& x2, D2& d_x2,
- const T3& x3, D3& d_x3,
- const T4& x4, D4& d_x4,
- const T5& x5, D5& d_x5,
- const T6& x6, D6& d_x6) {
- T deriv = 0;
- if (!is_constant_struct<T1>::value)
- deriv += increment_derivative<T, T1, D1>()(x1, d_x1);
- if (!is_constant_struct<T2>::value)
- deriv += increment_derivative<T, T2, D2>()(x2, d_x2);
- if (!is_constant_struct<T3>::value)
- deriv += increment_derivative<T, T3, D3>()(x3, d_x3);
- if (!is_constant_struct<T4>::value)
- deriv += increment_derivative<T, T4, D4>()(x4, d_x4);
- if (!is_constant_struct<T5>::value)
- deriv += increment_derivative<T, T5, D5>()(x5, d_x5);
- if (!is_constant_struct<T6>::value)
- deriv += increment_derivative<T, T6, D6>()(x6, d_x6);
- return fvar<T>(logp, deriv);
- }
- }
-
- /**
- * This class builds partial derivatives with respect to a set of
- * operands. There are two reason for the generality of this
- * class. The first is to handle vector and scalar arguments
- * without needing to write additional code. The second is to use
- * this class for writing probability distributions that handle
- * primitives, reverse mode, and forward mode variables
- * seamlessly.
- *
- * This is the partial template specialization for when the return
- * type is fvar<T>.
- *
- * @tparam T1 First set of operands.
- * @tparam T2 Second set of operands.
- * @tparam T3 Third set of operands.
- * @tparam T4 Fourth set of operands.
- * @tparam T5 Fifth set of operands.
- * @tparam T6 Sixth set of operands.
- * @tparam T_return_type Return type of the expression. This defaults
- * to a template metaprogram that calculates the scalar promotion of
- * T1 -- T6.
- */
- template<typename T1, typename T2, typename T3,
- typename T4, typename T5, typename T6,
- typename T_partials_return>
- struct OperandsAndPartials<T1, T2, T3, T4, T5, T6,
- fvar<T_partials_return> > {
- typedef fvar<T_partials_return> T_return_type;
-
- const T1& x1_;
- const T2& x2_;
- const T3& x3_;
- const T4& x4_;
- const T5& x5_;
- const T6& x6_;
-
- size_t n_partials;
- T_partials_return* all_partials;
-
- VectorView<T_partials_return,
- is_vector<T1>::value,
- is_constant_struct<T1>::value> d_x1;
- VectorView<T_partials_return,
- is_vector<T2>::value,
- is_constant_struct<T2>::value> d_x2;
- VectorView<T_partials_return,
- is_vector<T3>::value,
- is_constant_struct<T3>::value> d_x3;
- VectorView<T_partials_return,
- is_vector<T4>::value,
- is_constant_struct<T4>::value> d_x4;
- VectorView<T_partials_return,
- is_vector<T5>::value,
- is_constant_struct<T5>::value> d_x5;
- VectorView<T_partials_return,
- is_vector<T6>::value,
- is_constant_struct<T6>::value> d_x6;
-
- OperandsAndPartials(const T1& x1 = 0, const T2& x2 = 0, const T3& x3 = 0,
- const T4& x4 = 0, const T5& x5 = 0, const T6& x6 = 0)
- : x1_(x1), x2_(x2), x3_(x3), x4_(x4), x5_(x5), x6_(x6),
- n_partials(!is_constant_struct<T1>::value * length(x1) +
- !is_constant_struct<T2>::value * length(x2) +
- !is_constant_struct<T3>::value * length(x3) +
- !is_constant_struct<T4>::value * length(x4) +
- !is_constant_struct<T5>::value * length(x5) +
- !is_constant_struct<T6>::value * length(x6)),
- all_partials(new T_partials_return[n_partials]),
- d_x1(all_partials),
- d_x2(all_partials
- + (!is_constant_struct<T1>::value) * length(x1)),
- d_x3(all_partials
- + (!is_constant_struct<T1>::value) * length(x1)
- + (!is_constant_struct<T2>::value) * length(x2)),
- d_x4(all_partials
- + (!is_constant_struct<T1>::value) * length(x1)
- + (!is_constant_struct<T2>::value) * length(x2)
- + (!is_constant_struct<T3>::value) * length(x3)),
- d_x5(all_partials
- + (!is_constant_struct<T1>::value) * length(x1)
- + (!is_constant_struct<T2>::value) * length(x2)
- + (!is_constant_struct<T3>::value) * length(x3)
- + (!is_constant_struct<T4>::value) * length(x4)),
- d_x6(all_partials
- + (!is_constant_struct<T1>::value) * length(x1)
- + (!is_constant_struct<T2>::value) * length(x2)
- + (!is_constant_struct<T3>::value) * length(x3)
- + (!is_constant_struct<T4>::value) * length(x4)
- + (!is_constant_struct<T5>::value) * length(x5)) {
- std::fill(all_partials, all_partials + n_partials, 0);
- }
-
- T_return_type
- value(T_partials_return value) {
- return partials_to_fvar(value,
- x1_, d_x1, x2_, d_x2,
- x3_, d_x3, x4_, d_x4,
- x5_, d_x4, x6_, d_x5);
- }
-
- ~OperandsAndPartials() {
- delete[] all_partials;
- }
- };
-
- }
-}
-#endif
diff --git a/stan/math/fwd/scal/meta/operands_and_partials.hpp b/stan/math/fwd/scal/meta/operands_and_partials.hpp
index 77fc9b0c286d8d156948765a572db7fbadb3f3fb..c9dc55bed7ee174ffea85c26846fabfdb81a342c 100644
--- a/stan/math/fwd/scal/meta/operands_and_partials.hpp
+++ b/stan/math/fwd/scal/meta/operands_and_partials.hpp
@@ -15,80 +15,40 @@ namespace stan {
ops_partials_edge(const fvar<Dx>& op)
: ops_partials_edge_singular<ViewElt, fvar<Dx> >(op) {}
Dx dx() {
- return this->partial * this->operand.d_;
+ return this->partials[0] * this->operand.d_;
}
};
+ } // end namespace detail
+ template <typename Op1, typename Op2, typename Op3, typename Op4,
+ typename Dx>
+ class operands_and_partials<Op1, Op2, Op3, Op4, fvar<Dx> > {
+ public:
+ // these are going to be stack local and get collected
+ detail::ops_partials_edge<Dx, Op1> edge1_;
+ detail::ops_partials_edge<Dx, Op2> edge2_;
+ detail::ops_partials_edge<Dx, Op3> edge3_;
+ detail::ops_partials_edge<Dx, Op4> edge4_;
+ typedef fvar<Dx> T_return_type;
+ operands_and_partials(const Op1& o1)
+ : edge1_(o1) { }
+ operands_and_partials(const Op1& o1, const Op2& o2)
+ : edge1_(o1), edge2_(o2) { }
+ operands_and_partials(const Op1& o1, const Op2& o2, const Op3& o3)
+ : edge1_(o1), edge2_(o2), edge3_(o3) { }
+ operands_and_partials(const Op1& o1, const Op2& o2, const Op3& o3, const Op4& o4)
+ : edge1_(o1), edge2_(o2), edge3_(o3), edge4_(o4) { }
- template <typename Op1, typename Op2, typename Op3, typename Op4,
- typename Dx>
- class operands_and_partials<Op1, Op2, Op3, Op4, fvar<Dx> > {
- public:
- typedef fvar<Dx> T_return_type;
- operands_and_partials(const Op1& o1)
- : edge1_(o1) { }
- operands_and_partials(const Op1& o1, const Op2& o2)
- : edge1_(o1), edge2_(o2) { }
- operands_and_partials(const Op1& o1, const Op2& o2, const Op3& o3)
- : edge1_(o1), edge2_(o2), edge3_(o3) { }
- operands_and_partials(const Op1& o1, const Op2& o2, const Op3& o3, const Op4& o4)
- : edge1_(o1), edge2_(o2), edge3_(o3), edge4_(o4) { }
-
- void increment_dx1(const int n, const Dx& adj) {
- edge1_.increment_dx(n, adj);
- }
- void increment_dx2(const int n, const Dx& adj) {
- edge2_.increment_dx(n, adj);
- }
- void increment_dx3(const int n, const Dx& adj) {
- edge3_.increment_dx(n, adj);
- }
- void increment_dx4(const int n, const Dx& adj) {
- edge4_.increment_dx(n, adj);
- }
-
- template<typename T>
- typename boost::enable_if_c<is_vector_like<Op1>::value
- && is_vector_like<T>::value, void>::type
- increment_dx1_vector(const int n, const T& adj) {
- edge1_.increment_dx_vector(n, adj);
- }
- template<typename T>
- typename boost::enable_if_c<is_vector_like<Op2>::value
- && is_vector_like<T>::value, void>::type
- increment_dx2_vector(const int n, const T& adj) {
- edge2_.increment_dx_vector(n, adj);
- }
- template<typename T>
- typename boost::enable_if_c<is_vector_like<Op3>::value
- && is_vector_like<T>::value, void>::type
- increment_dx3_vector(const int n, const T& adj) {
- edge3_.increment_dx_vector(n, adj);
- }
- template<typename T>
- typename boost::enable_if_c<is_vector_like<Op4>::value
- && is_vector_like<T>::value, void>::type
- increment_dx4_vector(const int n, const T& adj) {
- edge4_.increment_dx_vector(n, adj);
- }
-
- // this is what matters in terms of going on autodiff stack
- T_return_type build(Dx value) {
- // Just call chain here and end up with an fvar(val, adj)
- Dx deriv = 0;
- deriv += edge1_.dx();
- deriv += edge2_.dx();
- deriv += edge3_.dx();
- deriv += edge4_.dx();
- return T_return_type(value, deriv);
- };
- private:
- // these are going to be stack local and get collected
- ops_partials_edge<Dx, Op1> edge1_;
- ops_partials_edge<Dx, Op2> edge2_;
- ops_partials_edge<Dx, Op3> edge3_;
- ops_partials_edge<Dx, Op4> edge4_;
+ // this is what matters in terms of going on autodiff stack
+ T_return_type build(Dx value) {
+ // Just call chain here and end up with an fvar(val, adj)
+ Dx deriv = 0;
+ deriv += edge1_.dx();
+ deriv += edge2_.dx();
+ deriv += edge3_.dx();
+ deriv += edge4_.dx();
+ return T_return_type(value, deriv);
};
- } // end namespace detail
+ };
}
}
diff --git a/stan/math/prim/arr/meta/VectorView.hpp b/stan/math/prim/arr/meta/VectorView.hpp
deleted file mode 100644
index 4759217b2477129d1f7c31bc78785f986bd745f2..0000000000000000000000000000000000000000
--- a/stan/math/prim/arr/meta/VectorView.hpp
+++ /dev/null
@@ -1,42 +0,0 @@
-#ifndef STAN_MATH_ARR_SCAL_META_VECTORVIEW_HPP
-#define STAN_MATH_ARR_SCAL_META_VECTORVIEW_HPP
-
-#include <stan/math/prim/scal/meta/VectorView.hpp>
-#include <vector>
-
-namespace stan {
-
- template <typename T>
- class VectorView<std::vector<T>, true, false> {
- public:
- typedef typename scalar_type<T>::type scalar_t;
-
- template <typename X>
- explicit VectorView(X& x) : x_(&x[0]) { }
-
- scalar_t& operator[](int i) {
- return x_[i];
- }
-
- private:
- scalar_t* x_;
- };
-
- template <typename T>
- class VectorView<const std::vector<T>, true, false> {
- public:
- typedef typename boost::add_const<typename scalar_type<T>::type>::type
- scalar_t;
-
- template <typename X>
- explicit VectorView(X& x) : x_(&x[0]) { }
-
- scalar_t& operator[](int i) const {
- return x_[i];
- }
- private:
- scalar_t* x_;
- };
-
-}
-#endif
diff --git a/stan/math/prim/mat/meta/VectorView.hpp b/stan/math/prim/mat/meta/VectorView.hpp
deleted file mode 100644
index d1fe67d9e932f2cb2bc8f0accc45e895b226ad2d..0000000000000000000000000000000000000000
--- a/stan/math/prim/mat/meta/VectorView.hpp
+++ /dev/null
@@ -1,43 +0,0 @@
-#ifndef STAN_MATH_MAT_SCAL_META_VECTORVIEW_HPP
-#define STAN_MATH_MAT_SCAL_META_VECTORVIEW_HPP
-
-#include <stan/math/prim/mat/fun/Eigen.hpp>
-#include <stan/math/prim/mat/meta/scalar_type.hpp>
-#include <stan/math/prim/scal/meta/VectorView.hpp>
-#include <boost/type_traits.hpp>
-
-namespace stan {
-
- template <typename T, int R, int C>
- class VectorView<Eigen::Matrix<T, R, C>, true, false> {
- public:
- typedef typename scalar_type<T>::type scalar_t;
-
- template <typename X>
- explicit VectorView(X& x) : x_(x.data()) { }
-
- scalar_t& operator[](int i) {
- return x_[i];
- }
- private:
- scalar_t* x_;
- };
-
- template <typename T, int R, int C>
- class VectorView<const Eigen::Matrix<T, R, C>, true, false> {
- public:
- typedef typename boost::add_const<typename scalar_type<T>::type>::type
- scalar_t;
-
- template <typename X>
- explicit VectorView(X& x) : x_(x.data()) { }
-
- scalar_t& operator[](int i) const {
- return x_[i];
- }
- private:
- scalar_t* x_;
- };
-
-}
-#endif
diff --git a/stan/math/prim/mat/meta/operands_and_partials.hpp b/stan/math/prim/mat/meta/operands_and_partials.hpp
index 79106441f6d3b4078f3714a6976a0ba228fbf6b5..56c49c9143c8d66bd4087f475c552ecafa038025 100644
--- a/stan/math/prim/mat/meta/operands_and_partials.hpp
+++ b/stan/math/prim/mat/meta/operands_and_partials.hpp
@@ -3,30 +3,48 @@
#include <stan/math/prim/mat/fun/typedefs.hpp>
#include <stan/math/prim/scal/meta/operands_and_partials.hpp>
+#include <stan/math/prim/scal/meta/broadcast_array.hpp>
#include <vector>
namespace stan {
namespace math {
namespace detail {
+ template <typename T, typename Orig, int R, int C>
+ struct zero_vec_or_mat {
+ typedef Eigen::Matrix<T, R, C> ret;
+ static ret zero(const Orig& in) {
+ return ret::Zero(in.rows(), in.cols());
+ }
+ };
+ template <typename T, typename Orig>
+ struct zero_vec_or_mat<T, Orig, -1, 1> {
+ typedef Eigen::Matrix<T, -1, 1> ret;
+ static ret zero(const Orig& in) {
+ return ret::Zero(in.size(), 1);
+ }
+ };
+ template <typename T, typename Orig>
+ struct zero_vec_or_mat<T, Orig, 1, -1> {
+ typedef Eigen::Matrix<T, 1, -1> ret;
+ static ret zero(const Orig& in) {
+ return ret::Zero(1, in.size());
+ }
+ };
+
// Handles all the vectorized cases for "views of scalars":
// ViewElt = double, Arg = std::vector<var> d_ holds vector<double>
// ViewElt = double, Arg = std::vector<double> d_ holds dummy
// Op will always be some container of vars
template <typename ViewElt, typename Op, int R, int C>
class ops_partials_edge_mat_prim {
- protected:
+ public:
typedef Eigen::Matrix<ViewElt, R, C> partials_t;
const Op& operands;
partials_t partials;
- public:
+ broadcast_array<partials_t> partials_vec;
ops_partials_edge_mat_prim(const Op& ops)
- : operands(ops), partials(partials_t::Zero(ops.size())) {}
- void increment_dx(int n, const ViewElt& adj) {
- partials[n] += adj;
- }
- void increment_dx_vector(int /*n*/, const partials_t& adj) {
- partials += adj;
- }
+ : operands(ops), partials(zero_vec_or_mat<ViewElt, Op, R, C>::zero(ops)),
+ partials_vec(partials) {}
void dump_partials(double* partials) {
for (int i = 0; i < this->partials.size(); ++i) {
partials[i] = this->partials(i);
@@ -41,21 +59,17 @@ namespace stan {
class ops_partials_edge_multivariate_prim {
public:
typedef Eigen::Matrix<ViewElt, Eigen::Dynamic, Eigen::Dynamic> partial_t;
+ std::vector<partial_t> partials_vec;
+ const std::vector<Op>& operands;
ops_partials_edge_multivariate_prim(const std::vector<Op>& ops)
- : partials(ops.size()), operands(ops) {
+ : partials_vec(ops.size()), operands(ops) {
for (size_t i = 0; i < ops.size(); ++i) {
- partials[i] = partial_t::Zero(ops[i].rows(), ops[i].cols());
+ partials_vec[i] = partial_t::Zero(ops[i].rows(), ops[i].cols());
}
}
- void increment_dx_vector(int n, const partial_t& adj) {
- partials[n] += adj;
- }
int size() {
return this->operands.size() * this->operands[0].size();
}
- protected:
- std::vector<partial_t> partials;
- const std::vector<Op>& operands;
};
}
}
diff --git a/stan/math/prim/scal/meta/OperandsAndPartials.hpp b/stan/math/prim/scal/meta/OperandsAndPartials.hpp
deleted file mode 100644
index 0d5e684c364b2e75fa338c986710342e5d45b8d6..0000000000000000000000000000000000000000
--- a/stan/math/prim/scal/meta/OperandsAndPartials.hpp
+++ /dev/null
@@ -1,74 +0,0 @@
-#ifndef STAN_MATH_PRIM_SCAL_META_OPERANDSANDPARTIALS_HPP
-#define STAN_MATH_PRIM_SCAL_META_OPERANDSANDPARTIALS_HPP
-
-#include <stan/math/prim/scal/meta/return_type.hpp>
-#include <stan/math/prim/scal/meta/VectorView.hpp>
-
-namespace stan {
- namespace math {
-
- /**
- * This class builds partial derivatives with respect to a set of
- * operands. There are two reason for the generality of this
- * class. The first is to handle vector and scalar arguments
- * without needing to write additional code. The second is to use
- * this class for writing probability distributions that handle
- * primitives, reverse mode, and forward mode variables
- * seamlessly.
- *
- * The default template class handles the case where the arguments
- * are primitive. There are template specializations for reverse
- * mode and forward mode.
- *
- * @tparam T1 First set of operands.
- * @tparam T2 Second set of operands.
- * @tparam T3 Third set of operands.
- * @tparam T4 Fourth set of operands.
- * @tparam T5 Fifth set of operands.
- * @tparam T6 Sixth set of operands.
- * @tparam T_return_type Return type of the expression. This defaults
- * to a template metaprogram that calculates the scalar promotion of
- * T1 -- T6.
- */
- template<typename T1 = double, typename T2 = double, typename T3 = double,
- typename T4 = double, typename T5 = double, typename T6 = double,
- typename T_return_type
- = typename stan::return_type<T1, T2, T3, T4, T5, T6>::type>
- struct OperandsAndPartials {
- VectorView<T_return_type, false, true> d_x1;
- VectorView<T_return_type, false, true> d_x2;
- VectorView<T_return_type, false, true> d_x3;
- VectorView<T_return_type, false, true> d_x4;
- VectorView<T_return_type, false, true> d_x5;
- VectorView<T_return_type, false, true> d_x6;
-
- /**
- * Constructor.
- *
- * @param x1 first set of operands
- * @param x2 second set of operands
- * @param x3 third set of operands
- * @param x4 fourth set of operands
- * @param x5 fifth set of operands
- * @param x6 sixth set of operands
- */
- OperandsAndPartials(const T1& x1 = 0, const T2& x2 = 0,
- const T3& x3 = 0, const T4& x4 = 0,
- const T5& x5 = 0, const T6& x6 = 0) { }
-
- /**
- * Returns a T_return_type with the value specified with
- * the partial derivatves.
- *
- * @param[in] value Value of the variable
- * @returns a variable with the appropriate value
- */
- T_return_type
- value(double value) {
- return value;
- }
- };
-
- }
-}
-#endif
diff --git a/stan/math/prim/scal/meta/VectorView.hpp b/stan/math/prim/scal/meta/VectorView.hpp
deleted file mode 100644
index 1a34b3c1e5a5c9f514195657834763987b13e604..0000000000000000000000000000000000000000
--- a/stan/math/prim/scal/meta/VectorView.hpp
+++ /dev/null
@@ -1,144 +0,0 @@
-#ifndef STAN_MATH_PRIM_SCAL_META_VECTORVIEW_HPP
-#define STAN_MATH_PRIM_SCAL_META_VECTORVIEW_HPP
-
-#include <stan/math/prim/scal/meta/scalar_type.hpp>
-#include <stan/math/prim/scal/meta/is_vector_like.hpp>
-#include <boost/type_traits.hpp>
-#include <stdexcept>
-
-namespace stan {
-
- /**
- * VectorView is a template expression that is constructed with a
- * container or scalar, which it then allows to be used as an array
- * using <code>operator[]</code>.
- *
- * For a scalar value, any index returns the reference or pointer
- * used to construct the view.
- *
- * For a container, the index returns a reference to the position in
- * the underlying container used to construct the view. WARNING:
- * There is no bounds checking for container indices and they will
- * segfault if accessed beyond their boundaries.
- *
- * The first use is to read arguments to prob functions as vectors,
- * even if scalars, so they can be read by common code (and scalars
- * automatically broadcast up to behave like vectors) : VectorView
- * of immutable const array of double* (no allocation).
- *
- * The second use is to build up derivatives into common storage :
- * VectorView of mutable shared array (no allocation because
- * allocated on auto-diff arena memory).
- *
- * Because it deals with references to its inputs, it is up to the
- * client of VectorView to ensure that the container being wrapped
- * is not modified while the VectorView is in use in such a way as
- * to disrupt the indexing. Similarly, because it deals with
- * references, it cannot be constructed with a literal or expression.
- *
- * @tparam T Type of scalar or container being wrapped.
- * @tparam is_array True if underlying type T can be indexed with
- * operator[].
- * @tparam throw_if_accessed True if the behavior is to throw an
- * exception whenever <code>operator[]</code> is called.
- */
- template <typename T,
- bool is_array = stan::is_vector_like<T>::value,
- bool throw_if_accessed = false>
- class VectorView {
- public:
- typedef typename
- boost::conditional<boost::is_const<T>::value,
- typename boost::add_const<
- typename scalar_type<T>::type>::type,
- typename scalar_type<T>::type>::type scalar_t;
-
- template <typename X>
- explicit VectorView(X x) {
- throw std::logic_error("VectorView: the default template "
- "specialization not implemented");
- }
-
- scalar_t& operator[](int i) {
- throw std::logic_error("VectorView: the default template "
- "specialization not implemented");
- }
-
- scalar_t& operator[](int i) const {
- throw std::logic_error("VectorView: the default template "
- "specialization not implemented");
- }
- };
-
- template <typename T, bool is_array>
- class VectorView<T, is_array, true> {
- public:
- typedef typename
- boost::conditional<boost::is_const<T>::value,
- typename boost::add_const<
- typename scalar_type<T>::type>::type,
- typename scalar_type<T>::type>::type scalar_t;
- VectorView() { }
-
- template <typename X>
- explicit VectorView(X x) { }
-
- scalar_t& operator[](int i) {
- throw std::logic_error("VectorView: this cannot be accessed");
- }
-
- scalar_t& operator[](int i) const {
- throw std::logic_error("VectorView: this cannot be accessed");
- }
- };
-
- // this covers non-vectors: double
- template <typename T>
- class VectorView<T, false, false> {
- public:
- typedef typename
- boost::conditional<boost::is_const<T>::value,
- typename boost::add_const<
- typename scalar_type<T>::type>::type,
- typename scalar_type<T>::type>::type scalar_t;
-
- explicit VectorView(scalar_t& x) : x_(&x) { }
-
- explicit VectorView(scalar_t* x) : x_(x) { }
-
- scalar_t& operator[](int i) {
- return *x_;
- }
-
- scalar_t& operator[](int i) const {
- return *x_;
- }
- private:
- scalar_t* x_;
- };
-
- // this covers raw memory: double*
- template <typename T>
- class VectorView<T, true, false> {
- public:
- typedef typename
- boost::conditional<boost::is_const<T>::value,
- typename boost::add_const<
- typename scalar_type<T>::type>::type,
- typename scalar_type<T>::type>::type scalar_t;
-
- explicit VectorView(scalar_t* x) : x_(x) { }
-
- scalar_t& operator[](int i) {
- return x_[i];
- }
-
- scalar_t& operator[](int i) const {
- return x_[i];
- }
-
- private:
- scalar_t* x_;
- };
-}
-#endif
diff --git a/stan/math/prim/scal/meta/broadcast_array.hpp b/stan/math/prim/scal/meta/broadcast_array.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..6be607e4c38d3399c446dd15390cc72eb0949262
--- /dev/null
+++ b/stan/math/prim/scal/meta/broadcast_array.hpp
@@ -0,0 +1,34 @@
+#ifndef STAN_MATH_PRIM_SCAL_META_BROADCAST_ARRAY_HPP
+#define STAN_MATH_PRIM_SCAL_META_BROADCAST_ARRAY_HPP
+
+namespace stan {
+ namespace math {
+ namespace detail {
+ template <typename T>
+ class broadcast_array {
+ private:
+ T& prim_;
+
+ public:
+ explicit broadcast_array(T& prim): prim_(prim) {}
+
+ T& operator[] (int /*i*/) {
+ return prim_;
+ }
+ int size() { return 1; }
+ };
+
+ template <>
+ class broadcast_array<void> {
+ public:
+ explicit broadcast_array() {}
+ double& operator[] (int /*i*/) {
+ throw std::logic_error("Don't do this");
+ }
+ int size() { return 0; }
+ };
+ } // end namespace detail
+ } // end namespace math
+} // end namespace stan
+
+#endif
diff --git a/stan/math/prim/scal/meta/operands_and_partials.hpp b/stan/math/prim/scal/meta/operands_and_partials.hpp
index a14a0e98427f959391753f766e87dcc0fe36c606..b49195a9565073d0c0e70090af4dd20f899bb8e3 100644
--- a/stan/math/prim/scal/meta/operands_and_partials.hpp
+++ b/stan/math/prim/scal/meta/operands_and_partials.hpp
@@ -2,6 +2,7 @@
#define STAN_MATH_PRIM_SCAL_META_OPERANDS_AND_PARTIALS_HPP
#include <stan/math/prim/scal/meta/return_type.hpp>
+#include <stan/math/prim/scal/meta/broadcast_array.hpp>
namespace stan {
namespace math {
@@ -52,40 +53,10 @@ namespace stan {
public:
ops_partials_edge() { }
ops_partials_edge(const Op& /* a */){ }
- void increment_dx(int /* n */, const ViewElt& /* adj */) {}
void dump_partials(ViewElt* /* partials */) {} // used for vars
void dump_operands(void* /* operands */) {} // also used for reverse mode
double dx() const { return 0; } //used for fvars
- int size() {return 0; }
- };
-
- template <typename Op1 = double, typename Op2 = double,
- typename Op3 = double, typename Op4 = double,
- typename T_return_type = typename return_type<Op1, Op2, Op3, Op4>::type>
- class operands_and_partials {
- public:
- operands_and_partials(const Op1& op1) {}
- operands_and_partials(const Op1& op1, const Op2& op2) {}
- operands_and_partials(const Op1& op1, const Op2& op2, const Op3& op3) {}
- operands_and_partials(const Op1& op1, const Op2& op2, const Op3& op3, const Op4& op4) {}
-
- void increment_dx1(int, double) {}
- void increment_dx2(int, double) {}
- void increment_dx3(int, double) {}
- void increment_dx4(int, double) {}
-
- template <typename PartialVec>
- void increment_dx1_vector(int, PartialVec) {}
- template <typename PartialVec>
- void increment_dx2_vector(int, PartialVec) {}
- template <typename PartialVec>
- void increment_dx3_vector(int, PartialVec) {}
- template <typename PartialVec>
- void increment_dx4_vector(int, PartialVec) {}
-
- double build(const double value) {
- return value;
- }
+ int size() { return 0; }
};
// Base class shared between fvar and var specializations of uni- and
@@ -93,21 +64,45 @@ namespace stan {
// Singular version - underlying Var is the same as Op, so there's just one.
template <typename ViewElt, typename Op>
class ops_partials_edge_singular {
- protected:
+ public:
const Op& operand;
ViewElt partial;
-
- public:
- ops_partials_edge_singular(const Op& a) : operand(a), partial(0) { }
- void increment_dx(int /*n*/, const ViewElt& adj) {
- partial += adj;
- }
+ broadcast_array<ViewElt> partials;
+ ops_partials_edge_singular(const Op& a) : operand(a), partial(0), partials(partial) { }
// dump_operands implemented in specialization
- int size() {
- return 1;
- }
+ int size() { return 1; }
+ };
+
+ template<>
+ struct ops_partials_edge<double, double> {
+ ops_partials_edge() {}
+ ops_partials_edge(const double& /* a */) {}
+ broadcast_array<void> partials;
+ void dump_operands(void*) {}
+ void dump_partials(void*) {}
+ double dx() { return 0; }
+ int size() { return 0; }
};
} // end namespace detail
+ template <typename Op1 = double, typename Op2 = double,
+ typename Op3 = double, typename Op4 = double,
+ typename T_return_type = typename return_type<Op1, Op2, Op3, Op4>::type>
+ class operands_and_partials {
+ public:
+ operands_and_partials(const Op1& op1) {}
+ operands_and_partials(const Op1& op1, const Op2& op2) {}
+ operands_and_partials(const Op1& op1, const Op2& op2, const Op3& op3) {}
+ operands_and_partials(const Op1& op1, const Op2& op2, const Op3& op3, const Op4& op4) {}
+
+ double build(const double value) {
+ return value;
+ }
+
+ detail::ops_partials_edge<double, Op1> edge1_;
+ detail::ops_partials_edge<double, Op2> edge2_;
+ detail::ops_partials_edge<double, Op3> edge3_;
+ detail::ops_partials_edge<double, Op4> edge4_;
+ };
} // end namespace math
} // end namespace stan
#endif
diff --git a/stan/math/prim/scal/meta/partials_return_type.hpp b/stan/math/prim/scal/meta/partials_return_type.hpp
index a173265cd4eb42c707bf3169ade103f2e8fdd747..9d068f4c4d98cee17a3fa466d0a7d055ec60b29b 100644
--- a/stan/math/prim/scal/meta/partials_return_type.hpp
+++ b/stan/math/prim/scal/meta/partials_return_type.hpp
@@ -25,7 +25,5 @@ namespace stan {
::type
type;
};
-
}
#endif
-
diff --git a/stan/math/prim/scal/meta/partials_type.hpp b/stan/math/prim/scal/meta/partials_type.hpp
index 0a56a8418862ef8e235def667d01519615403dd4..dd23d4275cb5200ca4cb5a6b365ed43ab4a3ba5a 100644
--- a/stan/math/prim/scal/meta/partials_type.hpp
+++ b/stan/math/prim/scal/meta/partials_type.hpp
@@ -10,4 +10,3 @@ namespace stan {
}
#endif
-
diff --git a/stan/math/rev/mat/meta/operands_and_partials.hpp b/stan/math/rev/mat/meta/operands_and_partials.hpp
index 4322d5817af7b57af866724f08aaed4e9c7f2666..d12477e97b016a3d01b4a57dd121364e7320823d 100644
--- a/stan/math/rev/mat/meta/operands_and_partials.hpp
+++ b/stan/math/rev/mat/meta/operands_and_partials.hpp
@@ -9,51 +9,31 @@
namespace stan {
namespace math {
namespace detail {
- template <typename ViewElt, typename Op>
- class ops_partials_edge_vec_rev
- : public ops_partials_edge_mat_prim<ViewElt, Op, -1, 1> {
- public:
- ops_partials_edge_vec_rev(const Op& ops)
- : ops_partials_edge_mat_prim<ViewElt, Op, -1, 1>(ops) {}
- void dump_partials(double* partials) {
- for (int i = 0; i < this->size(); ++i) {
- partials[i] = this->partials[i];
- }
- }
- void dump_operands(vari** varis) {
- for (int i = 0; i < this->size(); ++i) {
- varis[i] = this->operands[i].vi_;
- }
- }
- };
-
+ // Vectorized Univariate
// ViewElt = double, Arg = std::vector<var> d_ holds vector<double>
template <typename ViewElt>
class ops_partials_edge<ViewElt, std::vector<var> >
- : public ops_partials_edge_vec_rev<ViewElt, std::vector<var> > {
+ : public ops_partials_edge_mat_prim<ViewElt, std::vector<var>, -1, 1> {
public:
explicit ops_partials_edge(const std::vector<var>& ops)
- : ops_partials_edge_vec_rev<ViewElt, std::vector<var> >(ops) {}
- };
-
- // ViewElt = double, Arg = vector_v d_ holds vector<double>
- // ViewElt = double, Arg = VectorXd d_ holds dummy
- template <typename ViewElt>
- class ops_partials_edge<ViewElt, vector_v >
- : public ops_partials_edge_vec_rev<ViewElt, vector_v > {
- public:
- explicit ops_partials_edge(const vector_v& ops)
- : ops_partials_edge_vec_rev<ViewElt, vector_v>(ops) {}
+ : ops_partials_edge_mat_prim<ViewElt, std::vector<var>, -1, 1>(ops) {}
+ void dump_operands(vari** varis) {
+ for (size_t i = 0; i < this->operands.size(); ++i) {
+ varis[i] = this->operands[i].vi_;
+ }
+ }
};
-
- // ViewElt = double, Arg = row_vector_v d_ holds = vector<double>
- // ViewElt = double, Arg = RowVectorXd, d_ holds dummy
- template <typename ViewElt>
- class ops_partials_edge<ViewElt, row_vector_v >
- : public ops_partials_edge_vec_rev<ViewElt, row_vector_v > {
+ template <typename ViewElt, int R, int C>
+ class ops_partials_edge<ViewElt, Eigen::Matrix<var, R, C> >
+ : public ops_partials_edge_mat_prim<ViewElt, Eigen::Matrix<var, R, C>, R, C> {
public:
- explicit ops_partials_edge(const row_vector_v& ops)
- : ops_partials_edge_vec_rev<ViewElt, row_vector_v>(ops) {}
+ explicit ops_partials_edge(const Eigen::Matrix<var, R, C>& ops)
+ : ops_partials_edge_mat_prim<ViewElt, Eigen::Matrix<var, R, C>, R, C>(ops) {}
+ void dump_operands(vari** varis) {
+ for (int i = 0; i < this->operands.size(); ++i) {
+ varis[i] = this->operands(i).vi_;
+ }
+ }
};
// MULTIVARIATE
@@ -70,9 +50,9 @@ namespace stan {
: ops_partials_edge_multivariate_prim<ViewElt, Op>(ops) {}
void dump_partials(double* partials) {
int p_i = 0;
- for (size_t r = 0; r < this->partials.size(); ++r) {
- for (int c = 0; c < this->partials[r].size(); ++c, ++p_i) {
- partials[p_i] = this->partials[r](c);
+ for (size_t r = 0; r < this->partials_vec.size(); ++r) {
+ for (int c = 0; c < this->partials_vec[r].size(); ++c, ++p_i) {
+ partials[p_i] = this->partials_vec[r](c);
}
}
}
@@ -102,38 +82,6 @@ namespace stan {
: ops_partials_edge_multi_rev<ViewElt, vector_v>(ops) {}
};
- // VIEWS of MATRICES (R = -1, C = -1)
- // =====================
- // ViewElt = MatrixXd, Arg = matrix_v, d_ holds MatrixXd
- // ViewElt = MatrixXd, Arg = MatrixXd d_ holds dummy
- template <typename ViewElt, int R, int C>
- class ops_partials_edge<ViewElt, Eigen::Matrix<var, R, C> > {
- public:
- typedef Eigen::Matrix<ViewElt, R, C> partials_t;
- typedef Eigen::Matrix<var, R, C> operands_t;
- ops_partials_edge(const operands_t& ops)
- : operands(ops), partials(partials_t::Zero(ops.rows(), ops.cols())) {}
- void increment_dx_vector(int /*n*/, const partials_t& adj) {
- partials += adj;
- }
- void dump_partials(double* partials) {
- for (int i = 0; i < this->partials.size(); ++i) {
- partials[i] = this->partials(i);
- }
- }
- void dump_operands(vari** varis) {
- for (int i = 0; i < this->operands.size(); ++i) {
- varis[i] = this->operands(i).vi_;
- }
- }
- int size() {
- return this->operands.size();
- }
- protected:
- const operands_t& operands;
- partials_t partials;
- };
-
// ViewElt = MatrixXd, Arg = vector<matrix_v> d_ holds vector<MatriXd>
// ViewElt = MatrixXd, Arg = vector<MatrixXd> d_ holds dummy
template <typename ViewElt>
diff --git a/stan/math/rev/scal/meta/OperandsAndPartials.hpp b/stan/math/rev/scal/meta/OperandsAndPartials.hpp
deleted file mode 100644
index d978dedec92b3b43a21899008ee28b323cbde523..0000000000000000000000000000000000000000
--- a/stan/math/rev/scal/meta/OperandsAndPartials.hpp
+++ /dev/null
@@ -1,191 +0,0 @@
-#ifndef STAN_MATH_REV_SCAL_META_OPERANDSANDPARTIALS_HPP
-#define STAN_MATH_REV_SCAL_META_OPERANDSANDPARTIALS_HPP
-
-#include <stan/math/prim/scal/meta/is_constant_struct.hpp>
-#include <stan/math/prim/scal/meta/length.hpp>
-#include <stan/math/prim/scal/meta/OperandsAndPartials.hpp>
-#include <stan/math/prim/scal/meta/VectorView.hpp>
-#include <stan/math/rev/core.hpp>
-
-namespace stan {
- namespace math {
-
- namespace {
- class partials_vari : public vari {
- private:
- const size_t N_;
- vari** operands_;
- double* partials_;
- public:
- partials_vari(double value,
- size_t N,
- vari** operands, double* partials)
- : vari(value),
- N_(N),
- operands_(operands),
- partials_(partials) { }
- void chain() {
- for (size_t n = 0; n < N_; ++n)
- operands_[n]->adj_ += adj_ * partials_[n];
- }
- };
-
- var partials_to_var(double logp, size_t nvaris,
- vari** all_varis,
- double* all_partials) {
- return var(new partials_vari(logp, nvaris, all_varis,
- all_partials));
- }
-
- template<typename T,
- bool is_vec = is_vector<T>::value,
- bool is_const = is_constant_struct<T>::value>
- struct set_varis {
- inline size_t set(vari** /*varis*/, const T& /*x*/) {
- return 0U;
- }
- };
- template<typename T>
- struct set_varis<T, true, false> {
- inline size_t set(vari** varis, const T& x) {
- for (size_t n = 0; n < length(x); n++)
- varis[n] = x[n].vi_;
- return length(x);
- }
- };
- template<>
- struct set_varis<var, false, false> {
- inline size_t set(vari** varis, const var& x) {
- varis[0] = x.vi_;
- return (1);
- }
- };
- }
-
- /**
- * This class builds partial derivatives with respect to a set of
- * operands. There are two reason for the generality of this
- * class. The first is to handle vector and scalar arguments
- * without needing to write additional code. The second is to use
- * this class for writing probability distributions that handle
- * primitives, reverse mode, and forward mode variables
- * seamlessly.
- *
- * This is the partial template specialization for when the return
- * type is stan::math::var.
- *
- * @tparam T1 First set of operands.
- * @tparam T2 Second set of operands.
- * @tparam T3 Third set of operands.
- * @tparam T4 Fourth set of operands.
- * @tparam T5 Fifth set of operands.
- * @tparam T6 Sixth set of operands.
- * @tparam T_return_type Return type of the expression. This defaults
- * to a template metaprogram that calculates the scalar promotion of
- * T1 -- T6.
- */
- template<typename T1, typename T2, typename T3,
- typename T4, typename T5, typename T6>
- struct OperandsAndPartials<T1, T2, T3, T4, T5, T6, stan::math::var> {
- size_t nvaris;
- vari** all_varis;
- double* all_partials;
-
- VectorView<double,
- is_vector<T1>::value,
- is_constant_struct<T1>::value> d_x1;
- VectorView<double,
- is_vector<T2>::value,
- is_constant_struct<T2>::value> d_x2;
- VectorView<double,
- is_vector<T3>::value,
- is_constant_struct<T3>::value> d_x3;
- VectorView<double,
- is_vector<T4>::value,
- is_constant_struct<T4>::value> d_x4;
- VectorView<double,
- is_vector<T5>::value,
- is_constant_struct<T5>::value> d_x5;
- VectorView<double,
- is_vector<T6>::value,
- is_constant_struct<T6>::value> d_x6;
-
- /**
- * Constructor.
- *
- * @param x1 first set of operands
- * @param x2 second set of operands
- * @param x3 third set of operands
- * @param x4 fourth set of operands
- * @param x5 fifth set of operands
- * @param x6 sixth set of operands
- */
- OperandsAndPartials(const T1& x1 = 0, const T2& x2 = 0, const T3& x3 = 0,
- const T4& x4 = 0, const T5& x5 = 0, const T6& x6 = 0)
- : nvaris(!is_constant_struct<T1>::value * length(x1) +
- !is_constant_struct<T2>::value * length(x2) +
- !is_constant_struct<T3>::value * length(x3) +
- !is_constant_struct<T4>::value * length(x4) +
- !is_constant_struct<T5>::value * length(x5) +
- !is_constant_struct<T6>::value * length(x6)),
- // TODO(carpenter): replace with array allocation fun
- all_varis(static_cast<vari**>
- (vari::operator new
- (sizeof(vari*) * nvaris))),
- all_partials(static_cast<double*>
- (vari::operator new
- (sizeof(double) * nvaris))),
- d_x1(all_partials),
- d_x2(all_partials
- + (!is_constant_struct<T1>::value) * length(x1)),
- d_x3(all_partials
- + (!is_constant_struct<T1>::value) * length(x1)
- + (!is_constant_struct<T2>::value) * length(x2)),
- d_x4(all_partials
- + (!is_constant_struct<T1>::value) * length(x1)
- + (!is_constant_struct<T2>::value) * length(x2)
- + (!is_constant_struct<T3>::value) * length(x3)),
- d_x5(all_partials
- + (!is_constant_struct<T1>::value) * length(x1)
- + (!is_constant_struct<T2>::value) * length(x2)
- + (!is_constant_struct<T3>::value) * length(x3)
- + (!is_constant_struct<T4>::value) * length(x4)),
- d_x6(all_partials
- + (!is_constant_struct<T1>::value) * length(x1)
- + (!is_constant_struct<T2>::value) * length(x2)
- + (!is_constant_struct<T3>::value) * length(x3)
- + (!is_constant_struct<T4>::value) * length(x4)
- + (!is_constant_struct<T5>::value) * length(x5)) {
- size_t base = 0;
- if (!is_constant_struct<T1>::value)
- base += set_varis<T1>().set(&all_varis[base], x1);
- if (!is_constant_struct<T2>::value)
- base += set_varis<T2>().set(&all_varis[base], x2);
- if (!is_constant_struct<T3>::value)
- base += set_varis<T3>().set(&all_varis[base], x3);
- if (!is_constant_struct<T4>::value)
- base += set_varis<T4>().set(&all_varis[base], x4);
- if (!is_constant_struct<T5>::value)
- base += set_varis<T5>().set(&all_varis[base], x5);
- if (!is_constant_struct<T6>::value)
- set_varis<T6>().set(&all_varis[base], x6);
- std::fill(all_partials, all_partials+nvaris, 0);
- }
-
- /**
- * Returns a T_return_type with the value specified with
- * the partial derivatves.
- *
- * @param[in] value Value of the variable
- * @returns a variable with the appropriate value and
- * the adjoints set for reverse mode autodiff
- */
- stan::math::var value(double value) {
- return partials_to_var(value, nvaris, all_varis,
- all_partials);
- }
- };
-
- }
-}
-#endif
diff --git a/stan/math/rev/scal/meta/operands_and_partials.hpp b/stan/math/rev/scal/meta/operands_and_partials.hpp
index 16d59ba07beb91523baf25fc21fb4eff6ee4a82b..64f02a596c24c040ed4fbc89422c3c86a5656b1d 100644
--- a/stan/math/rev/scal/meta/operands_and_partials.hpp
+++ b/stan/math/rev/scal/meta/operands_and_partials.hpp
@@ -19,83 +19,50 @@ namespace stan {
ops_partials_edge(const var& op)
: ops_partials_edge_singular<ViewElt, var>(op) {}
void dump_partials(ViewElt* partials) {
- *partials = this->partial;
+ *partials = this->partials[0];
}
void dump_operands(vari** varis) {
*varis = this->operand.vi_;
}
};
+ } // end namespace detail
+ template <typename Op1, typename Op2, typename Op3, typename Op4>
+ class operands_and_partials<Op1, Op2, Op3, Op4, var> {
+ public:
+ // these are going to be stack local and get collected
+ detail::ops_partials_edge<double, Op1> edge1_;
+ detail::ops_partials_edge<double, Op2> edge2_;
+ detail::ops_partials_edge<double, Op3> edge3_;
+ detail::ops_partials_edge<double, Op4> edge4_;
+ operands_and_partials(const Op1& o1)
+ : edge1_(o1) { }
+ operands_and_partials(const Op1& o1, const Op2& o2)
+ : edge1_(o1), edge2_(o2) { }
+ operands_and_partials(const Op1& o1, const Op2& o2, const Op3& o3)
+ : edge1_(o1), edge2_(o2), edge3_(o3) { }
+ operands_and_partials(const Op1& o1, const Op2& o2, const Op3& o3, const Op4& o4)
+ : edge1_(o1), edge2_(o2), edge3_(o3), edge4_(o4) { }
- template <typename Op1, typename Op2, typename Op3, typename Op4>
- class operands_and_partials<Op1, Op2, Op3, Op4, var> {
- public:
- operands_and_partials(const Op1& o1)
- : edge1_(o1) { }
- operands_and_partials(const Op1& o1, const Op2& o2)
- : edge1_(o1), edge2_(o2) { }
- operands_and_partials(const Op1& o1, const Op2& o2, const Op3& o3)
- : edge1_(o1), edge2_(o2), edge3_(o3) { }
- operands_and_partials(const Op1& o1, const Op2& o2, const Op3& o3, const Op4& o4)
- : edge1_(o1), edge2_(o2), edge3_(o3), edge4_(o4) { }
-
- void increment_dx1(const int n, const double adj) { edge1_.increment_dx(n, adj); }
- void increment_dx2(const int n, const double adj) { edge2_.increment_dx(n, adj); }
- void increment_dx3(const int n, const double adj) { edge3_.increment_dx(n, adj); }
- void increment_dx4(const int n, const double adj) { edge4_.increment_dx(n, adj); }
-
- template<typename T>
- typename boost::enable_if_c<is_vector_like<Op1>::value
- && is_vector_like<T>::value, void>::type
- increment_dx1_vector(const int n, const T& adj) {
- edge1_.increment_dx_vector(n, adj);
- }
- template<typename T>
- typename boost::enable_if_c<is_vector_like<Op2>::value
- && is_vector_like<T>::value, void>::type
- increment_dx2_vector(const int n, const T& adj) {
- edge2_.increment_dx_vector(n, adj);
- }
- template<typename T>
- typename boost::enable_if_c<is_vector_like<Op3>::value
- && is_vector_like<T>::value, void>::type
- increment_dx3_vector(const int n, const T& adj) {
- edge3_.increment_dx_vector(n, adj);
- }
- template<typename T>
- typename boost::enable_if_c<is_vector_like<Op4>::value
- && is_vector_like<T>::value, void>::type
- increment_dx4_vector(const int n, const T& adj) {
- edge4_.increment_dx_vector(n, adj);
- }
-
- // this is what matters in terms of going on autodiff stack
- var build(double value) {
- size_t size = edge1_.size() + edge2_.size() + edge3_.size() + edge4_.size();
- vari** varis = ChainableStack::memalloc_.alloc_array<vari*>(size);
- int idx = 0;
- edge1_.dump_operands(&varis[idx]);
- edge2_.dump_operands(&varis[idx += edge1_.size()]);
- edge3_.dump_operands(&varis[idx += edge2_.size()]);
- edge4_.dump_operands(&varis[idx += edge3_.size()]);
+ // this is what matters in terms of going on autodiff stack
+ var build(double value) {
+ size_t size = edge1_.size() + edge2_.size() + edge3_.size() + edge4_.size();
+ vari** varis = ChainableStack::memalloc_.alloc_array<vari*>(size);
+ int idx = 0;
+ edge1_.dump_operands(&varis[idx]);
+ edge2_.dump_operands(&varis[idx += edge1_.size()]);
+ edge3_.dump_operands(&varis[idx += edge2_.size()]);
+ edge4_.dump_operands(&varis[idx += edge3_.size()]);
- double* partials = ChainableStack::memalloc_.alloc_array<double>(size);
- idx = 0;
- edge1_.dump_partials(&partials[idx]);
- edge2_.dump_partials(&partials[idx += edge1_.size()]);
- edge3_.dump_partials(&partials[idx += edge2_.size()]);
- edge4_.dump_partials(&partials[idx += edge3_.size()]);
+ double* partials = ChainableStack::memalloc_.alloc_array<double>(size);
+ idx = 0;
+ edge1_.dump_partials(&partials[idx]);
+ edge2_.dump_partials(&partials[idx += edge1_.size()]);
+ edge3_.dump_partials(&partials[idx += edge2_.size()]);
+ edge4_.dump_partials(&partials[idx += edge3_.size()]);
- return var(new precomputed_gradients_vari(value, size, varis, partials));
- };
- private:
- // these are going to be stack local and get collected
- // TODO(sean): should we pass in ViewElt somehow? shape of it?
- ops_partials_edge<double, Op1> edge1_;
- ops_partials_edge<double, Op2> edge2_;
- ops_partials_edge<double, Op3> edge3_;
- ops_partials_edge<double, Op4> edge4_;
+ return var(new precomputed_gradients_vari(value, size, varis, partials));
};
- } // end namespace detail
+ };
}
}
#endif
diff --git a/test/foo_test.cpp b/test/foo_test.cpp
deleted file mode 100644
index b7368fbb9bfb5865df54aeba57b7b45ce3aa3875..0000000000000000000000000000000000000000
--- a/test/foo_test.cpp
+++ /dev/null
@@ -1,11 +0,0 @@
-#include <stan/math/mix/scal.hpp>
-#include <gtest/gtest.h>
-
-TEST(foo, bar) {
- double y = 0;
- double mu = 0;
- stan::math::fvar<stan::math::var> sigma = 1;
-
-
- stan::math::fvar<stan::math::var> lp = normal_lpdf(y, mu, sigma);
-}
diff --git a/test/unit/math/fwd/arr/meta/OperandsAndPartials_test.cpp b/test/unit/math/fwd/arr/meta/OperandsAndPartials_test.cpp
deleted file mode 100644
index 3cdc0d1406fd9c80d6c69c994b6fc313766b3dc8..0000000000000000000000000000000000000000
--- a/test/unit/math/fwd/arr/meta/OperandsAndPartials_test.cpp
+++ /dev/null
@@ -1,58 +0,0 @@
-#include <stan/math/fwd/arr.hpp>
-#include <gtest/gtest.h>
-
-TEST(AgradPartialsVari, OperandsAndPartialsFvarVec) {
- using stan::math::OperandsAndPartials;
- using stan::math::fvar;
-
- std::vector<fvar<double> > x1;
- x1.push_back(fvar<double>(2.0,2.0));
- x1.push_back(fvar<double>(1.0,3.0));
-
- fvar<double> x2 = 3.0;
- fvar<double> x3 = 5.0;
- x2.d_ = -1.0;
- x3.d_ = 4.0;
-
- OperandsAndPartials<std::vector<fvar<double> >,fvar<double>,fvar<double> > o(x1, x2, x3);
- o.d_x1[0] += 17.0;
- o.d_x1[1] += 13.0;
- o.d_x2[0] += 19.0;
- o.d_x2[0] += 19.0;
- o.d_x3[0] += 23.0;
- o.d_x3[0] += 23.0;
- fvar<double> y = o.value(-1.0);
-
- EXPECT_FLOAT_EQ(2*17 + 3*13 - 2*19 + 2*4*23,y.d_);
- EXPECT_FLOAT_EQ(-1,y.val_);
-}
-
-// TEST(AgradPartialsVari, incr_deriv_vec_fvar) {
-// using stan::VectorView;
-// using stan::math::incr_deriv;
-// using stan::is_vector;
-// using stan::is_constant_struct;
-// using stan::math::fvar;
-
-// fvar<double> c(1,1);
-
-// std::vector<fvar<double> > d;
-// d.push_back(c);
-// d.push_back(c);
-
-// std::vector<double> d_deriv;
-// d_deriv.push_back(3);
-// d_deriv.push_back(4);
-
-// VectorView<std::vector<double>,
-// stan::is_vector<std::vector<fvar<double> > >::value,
-// stan::is_constant_struct<std::vector<fvar<double> > >::value>
-// d_d(d_deriv);
-
-// double result3 = incr_deriv<VectorView<std::vector<double>,
-// is_vector<std::vector<fvar<double> > >::value,
-// is_constant_struct<std::vector<fvar<double> > >::value>,
-// std::vector<fvar<double> >,double>().incr(d_d,d);
-
-// EXPECT_FLOAT_EQ(7, result3);
-// }
diff --git a/test/unit/math/fwd/mat/meta/operands_and_partials_test.cpp b/test/unit/math/fwd/mat/meta/operands_and_partials_test.cpp
index 9e79c13daf2d5a41b37e6230c21d91e8809fbe35..66e1e65d19f9424ec292af473bcf4058d23a3de6 100644
--- a/test/unit/math/fwd/mat/meta/operands_and_partials_test.cpp
+++ b/test/unit/math/fwd/mat/meta/operands_and_partials_test.cpp
@@ -2,7 +2,7 @@
#include <gtest/gtest.h>
TEST(AgradPartialsVari, OperandsAndPartialsFvarVec) {
- using stan::math::detail::operands_and_partials;
+ using stan::math::operands_and_partials;
using stan::math::fvar;
std::vector<fvar<double> > x1;
@@ -18,11 +18,11 @@ TEST(AgradPartialsVari, OperandsAndPartialsFvarVec) {
dx1 << 17.0, 13.0;
operands_and_partials<std::vector<fvar<double> >,fvar<double>,fvar<double> > o(x1, x2, x3);
- o.increment_dx1_vector(0, dx1);
- o.increment_dx2(0, 19.0);
- o.increment_dx2(0, 19.0);
- o.increment_dx3(0, 23.0);
- o.increment_dx3(0, 23.0);
+ o.edge1_.partials_vec[0] += dx1;
+ o.edge2_.partials[0] += 19.0;
+ o.edge2_.partials[0] += 19.0;
+ o.edge3_.partials[0] += 23.0;
+ o.edge3_.partials[0] += 23.0;
fvar<double> y = o.build(-1.0);
EXPECT_FLOAT_EQ(2*17 + 3*13 - 2*19 + 2*4*23,y.d_);
diff --git a/test/unit/math/fwd/scal/meta/OperandsAndPartials_test.cpp b/test/unit/math/fwd/scal/meta/OperandsAndPartials_test.cpp
deleted file mode 100644
index 44150748c1f4c13b984ed612d315c86a212f999b..0000000000000000000000000000000000000000
--- a/test/unit/math/fwd/scal/meta/OperandsAndPartials_test.cpp
+++ /dev/null
@@ -1,48 +0,0 @@
-#include <stan/math/fwd/scal.hpp>
-#include <gtest/gtest.h>
-
-TEST(AgradPartialsVari, OperandsAndPartialsFvar) {
- using stan::math::OperandsAndPartials;
- using stan::math::fvar;
-
- fvar<double> x1 = 2.0;
- fvar<double> x2 = 3.0;
- fvar<double> x3 = 5.0;
- x1.d_ = 2.0;
- x2.d_ = -1.0;
- x3.d_ = 4.0;
-
- OperandsAndPartials<fvar<double>,fvar<double>,fvar<double> > o(x1, x2, x3);
- o.d_x1[0] += 17.0;
- o.d_x2[0] += 19.0;
- o.d_x3[0] += 23.0;
-
- fvar<double> y = o.value(-1.0);
- EXPECT_FLOAT_EQ(107,y.d_);
- EXPECT_FLOAT_EQ(-1,y.val_);
-}
-
-// TEST(AgradPartialsVari, incr_deriv_fvar) {
-// using stan::VectorView;
-// using stan::math::incr_deriv;
-// using stan::is_vector;
-// using stan::is_constant_struct;
-// using stan::math::fvar;
-
-// fvar<double> c;
-// c.val_ = 1.0;
-// c.d_ = 1.0;
-
-// double c_deriv = 2;
-// VectorView<const double,
-// stan::is_vector<fvar<double> >::value,
-// stan::is_constant_struct<fvar<double> >::value>
-// d_c(c_deriv);
-
-// double result2 = incr_deriv<VectorView<const double,
-// is_vector<fvar<double> >::value,
-// is_constant_struct<fvar<double> >::value>,
-// fvar<double>,double>().incr(d_c,c);
-
-// EXPECT_FLOAT_EQ(2, result2);
-// }
diff --git a/test/unit/math/fwd/scal/meta/operands_and_partials_test.cpp b/test/unit/math/fwd/scal/meta/operands_and_partials_test.cpp
index 24f22bb15bee61c3f230713ec495db826693511e..719fea335e273b892dbf16b87baaec177b06777c 100644
--- a/test/unit/math/fwd/scal/meta/operands_and_partials_test.cpp
+++ b/test/unit/math/fwd/scal/meta/operands_and_partials_test.cpp
@@ -2,7 +2,7 @@
#include <gtest/gtest.h>
TEST(AgradPartialsVari, OperandsAndPartialsFvar) {
- using stan::math::detail::operands_and_partials;
+ using stan::math::operands_and_partials;
using stan::math::fvar;
fvar<double> x1 = 2.0;
@@ -13,9 +13,9 @@ TEST(AgradPartialsVari, OperandsAndPartialsFvar) {
x3.d_ = 4.0;
operands_and_partials<fvar<double>,fvar<double>,fvar<double> > o(x1, x2, x3);
- o.increment_dx1(0, 17.0);
- o.increment_dx2(0, 19.0);
- o.increment_dx3(0, 23.0);
+ o.edge1_.partials[0] += 17.0;
+ o.edge2_.partials[0] += 19.0;
+ o.edge3_.partials[0] += 23.0;
fvar<double> y = o.build(-1.0);
EXPECT_FLOAT_EQ(107,y.d_);
diff --git a/test/unit/math/mix/scal/meta/broadcast_array_test.cpp b/test/unit/math/mix/scal/meta/broadcast_array_test.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ff6cfbff0e25453622192d2392ec9ff2c3dd9786
--- /dev/null
+++ b/test/unit/math/mix/scal/meta/broadcast_array_test.cpp
@@ -0,0 +1,15 @@
+#include <gtest/gtest.h>
+#include <stan/math/fwd/core/fvar.hpp>
+#include <stan/math/rev/core/var.hpp>
+#include <stan/math/prim/scal/meta/broadcast_array.hpp>
+
+TEST(foo, bar) {
+ using stan::math::detail::broadcast_array;
+ using stan::math::fvar;
+ using stan::math::var;
+
+ fvar<var> fv(0.0);
+ broadcast_array<fvar<var> > ba(fv);
+ EXPECT_EQ(0.0, ba[0].val_.vi_->val_);
+ EXPECT_EQ(0.0, ba[2].val_.vi_->val_);
+}
diff --git a/test/unit/math/prim/arr/meta/OperandsAndPartials_test.cpp b/test/unit/math/prim/arr/meta/OperandsAndPartials_test.cpp
deleted file mode 100644
index f49dcf8540e4257eae8fe9488a2871590a9573af..0000000000000000000000000000000000000000
--- a/test/unit/math/prim/arr/meta/OperandsAndPartials_test.cpp
+++ /dev/null
@@ -1,41 +0,0 @@
-#include <stan/math/prim/arr.hpp>
-#include <gtest/gtest.h>
-
-// TEST(AgradPartialsVari, incr_deriv_vec_double) {
-// using stan::VectorView;
-// using stan::math::incr_deriv;
-// using stan::is_vector;
-// using stan::is_constant_struct;
-
-// std::vector<double> b;
-// b.push_back(1);
-// b.push_back(1);
-
-// VectorView<const std::vector<double>,
-// stan::is_vector<std::vector<double> >::value,
-// stan::is_constant_struct<std::vector<double> >::value>
-// d_b(stan::length(b) * 0);
-
-// double result = incr_deriv<VectorView<const std::vector<double>,
-// is_vector<std::vector<double> >::value,
-// is_constant_struct<std::vector<double> >::value>,
-// std::vector<double>,double>().incr(d_b,b);
-
-// EXPECT_FLOAT_EQ(0, result);
-// }
-
-TEST(AgradPartialsVari, OperandsAndPartials_check_throw) {
- using stan::math::OperandsAndPartials;
- using std::vector;
-
- vector<double> D;
-
- OperandsAndPartials<vector<double>,vector<double>,vector<double>,
- vector<double>,vector<double>,vector<double> > o3(D,D,D,D,D,D);
- EXPECT_THROW(o3.d_x1[0], std::logic_error);
- EXPECT_THROW(o3.d_x2[0], std::logic_error);
- EXPECT_THROW(o3.d_x3[0], std::logic_error);
- EXPECT_THROW(o3.d_x4[0], std::logic_error);
- EXPECT_THROW(o3.d_x5[0], std::logic_error);
- EXPECT_THROW(o3.d_x6[0], std::logic_error);
-}
diff --git a/test/unit/math/prim/scal/meta/OperandsAndPartials_test.cpp b/test/unit/math/prim/scal/meta/OperandsAndPartials_test.cpp
deleted file mode 100644
index d37c4b96f15d0df143ec53d70496d9fe9da0e0c2..0000000000000000000000000000000000000000
--- a/test/unit/math/prim/scal/meta/OperandsAndPartials_test.cpp
+++ /dev/null
@@ -1,46 +0,0 @@
-#include <stan/math/prim/scal.hpp>
-#include <gtest/gtest.h>
-
-TEST(AgradPartialsVari, OperandsAndPartials) {
- using stan::math::OperandsAndPartials;
-
- OperandsAndPartials<double> o1;
- OperandsAndPartials<double,double,double,double> o2;
-
- SUCCEED() << "Construction should be ok.";
-}
-
-// TEST(AgradPartialsVari, incr_deriv_double) {
-// using stan::VectorView;
-// using stan::math::incr_deriv;
-// using stan::is_vector;
-// using stan::is_constant_struct;
-
-// double a = 1;
-
-// VectorView<const double,stan::is_vector<double>::value,
-// stan::is_constant_struct<double>::value> d_a(stan::length(a) * 0);
-
-
-
-// double result = incr_deriv<VectorView<const double,
-// is_vector<double>::value,
-// is_constant_struct<double>::value>,
-// double,double>().incr(d_a,a);
-// EXPECT_FLOAT_EQ(0, result);
-// }
-
-
-TEST(AgradPartialsVari, OperandsAndPartials_check_throw) {
- using stan::math::OperandsAndPartials;
-
- double d;
-
- OperandsAndPartials<> o1(d,d,d,d,d,d);
- EXPECT_THROW(o1.d_x1[0], std::logic_error);
- EXPECT_THROW(o1.d_x2[0], std::logic_error);
- EXPECT_THROW(o1.d_x3[0], std::logic_error);
- EXPECT_THROW(o1.d_x4[0], std::logic_error);
- EXPECT_THROW(o1.d_x5[0], std::logic_error);
- EXPECT_THROW(o1.d_x6[0], std::logic_error);
-}
diff --git a/test/unit/math/prim/scal/meta/operands_and_partials_test.cpp b/test/unit/math/prim/scal/meta/operands_and_partials_test.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d79c625c43fd9cea69c511b28b45d19f810a4db3
--- /dev/null
+++ b/test/unit/math/prim/scal/meta/operands_and_partials_test.cpp
@@ -0,0 +1,14 @@
+#include <stan/math/prim/scal.hpp>
+#include <gtest/gtest.h>
+
+TEST(AgradPartialsVari, OperandsAndPartials) {
+ using stan::math::OperandsAndPartials;
+
+ OperandsAndPartials<double> o1;
+ OperandsAndPartials<double,double,double,double> o2;
+
+ // TODO(Sean): Learn why it's size 6
+ EXPECT_GT(sizeof(&o2), sizeof(o2));
+
+ SUCCEED() << "Construction should be ok.";
+}
diff --git a/test/unit/math/rev/arr/meta/OperandsAndPartials_test.cpp b/test/unit/math/rev/arr/meta/OperandsAndPartials_test.cpp
deleted file mode 100644
index 7fe30453d6c9dbad245e0bf6b67eb05d25fce283..0000000000000000000000000000000000000000
--- a/test/unit/math/rev/arr/meta/OperandsAndPartials_test.cpp
+++ /dev/null
@@ -1,49 +0,0 @@
-#include <stan/math/rev/arr.hpp>
-#include <gtest/gtest.h>
-
-TEST(AgradPartialsVari, OperandsAndPartials) {
- using stan::math::OperandsAndPartials;
- using stan::math::var;
-
- std::vector<double> d_vec(4);
- OperandsAndPartials<std::vector<double> > o3(d_vec);
- std::vector<var> v_vec;
- v_vec.push_back(var(0.0));
- v_vec.push_back(var(1.0));
- v_vec.push_back(var(2.0));
- v_vec.push_back(var(3.0));
-
- std::vector<double> grad;
-
- OperandsAndPartials<std::vector<var> > o4(v_vec);
- o4.d_x1[0] = 10.0;
- o4.d_x1[1] = 20.0;
- o4.d_x1[2] = 30.0;
- o4.d_x1[3] = 40.0;
-
- var v = o4.value(10.0);
- v.grad(v_vec, grad);
- EXPECT_EQ(4U, o4.nvaris);
- EXPECT_FLOAT_EQ(10.0, v.val());
- EXPECT_FLOAT_EQ(10.0, grad[0]);
- EXPECT_FLOAT_EQ(20.0, grad[1]);
- EXPECT_FLOAT_EQ(30.0, grad[2]);
- EXPECT_FLOAT_EQ(40.0, grad[3]);
-}
-
-TEST(AgradPartialsVari, OperandsAndPartials_check_throw) {
- using stan::math::OperandsAndPartials;
- using stan::math::var;
- using std::vector;
-
- vector<var> V;
-
- OperandsAndPartials<vector<var>,vector<var>,vector<var>,
- vector<var>,vector<var>,vector<var> > o4(V,V,V,V,V,V);
- EXPECT_NO_THROW(o4.d_x1[0]);
- EXPECT_NO_THROW(o4.d_x2[0]);
- EXPECT_NO_THROW(o4.d_x3[0]);
- EXPECT_NO_THROW(o4.d_x4[0]);
- EXPECT_NO_THROW(o4.d_x5[0]);
- EXPECT_NO_THROW(o4.d_x6[0]);
-}
diff --git a/test/unit/math/rev/mat/meta/operands_and_partials_test.cpp b/test/unit/math/rev/mat/meta/operands_and_partials_test.cpp
index 491c1337672033121059f71bd27166bce9d4fc65..a3d4ff3b10130d136323ed6d3112bbddd83a2b86 100644
--- a/test/unit/math/rev/mat/meta/operands_and_partials_test.cpp
+++ b/test/unit/math/rev/mat/meta/operands_and_partials_test.cpp
@@ -5,14 +5,14 @@
#include <iostream>
TEST(AgradPartialsVari, OperandsAndPartialsVec) {
- using stan::math::detail::operands_and_partials;
+ using stan::math::operands_and_partials;
using stan::math::var;
using stan::math::vector_v;
using stan::math::vector_d;
vector_d d_vec(4);
operands_and_partials<vector_d > o3(d_vec);
- EXPECT_EQ(1, sizeof(o3));
+ EXPECT_EQ(4, sizeof(o3));
vector_v v_vec(4);
var v1 = var(0.0);
@@ -28,10 +28,10 @@ TEST(AgradPartialsVari, OperandsAndPartialsVec) {
v_stdvec.push_back(v4);
operands_and_partials<vector_v > o4(v_vec);
- o4.increment_dx1(0, 10.0);
- o4.increment_dx1(1, 20.0);
- o4.increment_dx1(2, 30.0);
- o4.increment_dx1(3, 40.0);
+ o4.edge1_.partials[0] += 10.0;
+ o4.edge1_.partials[1] += 20.0;
+ o4.edge1_.partials[2] += 30.0;
+ o4.edge1_.partials[3] += 40.0;
std::vector<double> grad;
var v = o4.build(10.0);
@@ -43,8 +43,43 @@ TEST(AgradPartialsVari, OperandsAndPartialsVec) {
EXPECT_FLOAT_EQ(40.0, grad[3]);
}
+TEST(AgradPartialsVari, OperandsAndPartialsStdVec) {
+ using stan::math::operands_and_partials;
+ using stan::math::var;
+
+ std::vector<double> d_vec(4);
+ operands_and_partials<std::vector<double> > o3(d_vec);
+ EXPECT_EQ(4, sizeof(o3));
+
+ std::vector<var> v_vec;
+ var v1 = var(0.0);
+ var v2 = var(1.0);
+ var v3 = var(2.0);
+ var v4 = var(3.0);
+ v_vec.push_back(v1);
+ v_vec.push_back(v2);
+ v_vec.push_back(v3);
+ v_vec.push_back(v4);
+
+ operands_and_partials<std::vector<var> > o4(v_vec);
+ o4.edge1_.partials[0] += 10.0;
+ o4.edge1_.partials[1] += 20.0;
+ o4.edge1_.partials[2] += 30.0;
+ o4.edge1_.partials[3] += 40.0;
+
+ std::vector<double> grad;
+ var v = o4.build(10.0);
+ v.grad(v_vec, grad);
+ EXPECT_FLOAT_EQ(10.0, v.val());
+ EXPECT_FLOAT_EQ(10.0, grad[0]);
+ EXPECT_FLOAT_EQ(20.0, grad[1]);
+ EXPECT_FLOAT_EQ(30.0, grad[2]);
+ EXPECT_FLOAT_EQ(40.0, grad[3]);
+}
+
+
TEST(AgradPartialsVari, OperandsAndPartialsMat) {
- using stan::math::detail::operands_and_partials;
+ using stan::math::operands_and_partials;
using stan::math::var;
using stan::math::matrix_v;
using stan::math::matrix_d;
@@ -53,7 +88,7 @@ TEST(AgradPartialsVari, OperandsAndPartialsMat) {
d_mat << 10.0, 20.0, 30.0, 40.0;
operands_and_partials<matrix_d > o3(d_mat);
- EXPECT_EQ(1, sizeof(o3));
+ EXPECT_EQ(4, sizeof(o3));
matrix_v v_mat(2, 2);
var v1 = var(0.0);
@@ -69,22 +104,22 @@ TEST(AgradPartialsVari, OperandsAndPartialsMat) {
v_stdvec.push_back(v4);
operands_and_partials<matrix_v > o4(v_mat);
- // TODO(Sean): Check assumption that Matrices are only used in the multivariate cases
- o4.increment_dx1_vector(1, d_mat);
- o4.increment_dx1_vector(27, d_mat); // Should affect the same vars as the call above
+ o4.edge1_.partials += d_mat;
+ o4.edge1_.partials_vec[1] += d_mat;
+ o4.edge1_.partials_vec[27] += d_mat; // Should affect the same vars as the call above
std::vector<double> grad;
var v = o4.build(10.0);
v.grad(v_stdvec, grad);
EXPECT_FLOAT_EQ(10.0, v.val());
- EXPECT_FLOAT_EQ(20.0, grad[0]);
- EXPECT_FLOAT_EQ(40.0, grad[1]);
- EXPECT_FLOAT_EQ(60.0, grad[2]);
- EXPECT_FLOAT_EQ(80.0, grad[3]);
+ EXPECT_FLOAT_EQ(30.0, grad[0]);
+ EXPECT_FLOAT_EQ(60.0, grad[1]);
+ EXPECT_FLOAT_EQ(90.0, grad[2]);
+ EXPECT_FLOAT_EQ(120.0, grad[3]);
}
TEST(AgradPartialsVari, OperandsAndPartialsMatMultivar) {
- using stan::math::detail::operands_and_partials;
+ using stan::math::operands_and_partials;
using stan::math::var;
using stan::math::matrix_v;
using stan::math::matrix_d;
@@ -95,7 +130,7 @@ TEST(AgradPartialsVari, OperandsAndPartialsMatMultivar) {
d_mat_vec.push_back(d_mat);
operands_and_partials<std::vector<matrix_d> > o3(d_mat_vec);
- EXPECT_EQ(1, sizeof(o3));
+ EXPECT_EQ(4, sizeof(o3));
matrix_v v_mat1(2, 2);
var v1 = var(0.0);
@@ -126,9 +161,9 @@ TEST(AgradPartialsVari, OperandsAndPartialsMatMultivar) {
v_stdvec.push_back(v8);
operands_and_partials<std::vector<matrix_v> > o4(v_mat_vec);
- o4.increment_dx1_vector(0, d_mat);
- o4.increment_dx1_vector(1, d_mat); // Should NOT affect the same vars as the call above
- o4.increment_dx1_vector(1, d_mat);
+ o4.edge1_.partials_vec[0] += d_mat;
+ o4.edge1_.partials_vec[1] += d_mat; // Should NOT affect the same vars as the call above
+ o4.edge1_.partials_vec[1] += d_mat;
std::vector<double> grad;
var v = o4.build(10.0);
@@ -145,7 +180,7 @@ TEST(AgradPartialsVari, OperandsAndPartialsMatMultivar) {
}
TEST(AgradPartialsVari, OperandsAndPartialsMultivar) {
- using stan::math::detail::operands_and_partials;
+ using stan::math::operands_and_partials;
using stan::math::var;
using stan::math::vector_v;
using stan::math::vector_d;
@@ -159,7 +194,7 @@ TEST(AgradPartialsVari, OperandsAndPartialsMultivar) {
d_vec_vec.push_back(d_vec2);
operands_and_partials<std::vector<vector_d> > o3(d_vec_vec);
- EXPECT_EQ(1, sizeof(o3));
+ EXPECT_EQ(4, sizeof(o3));
vector_v v_vec1(2);
var v1 = var(0.0);
@@ -180,8 +215,8 @@ TEST(AgradPartialsVari, OperandsAndPartialsMultivar) {
v_stdvec.push_back(v4);
operands_and_partials<std::vector<vector_v> > o4(v_vec);
- o4.increment_dx1_vector(0, d_vec1);
- o4.increment_dx1_vector(1, d_vec2);
+ o4.edge1_.partials_vec[0] += d_vec1;
+ o4.edge1_.partials_vec[1] += d_vec2;
std::vector<double> grad;
var v = o4.build(10.0);
@@ -195,7 +230,7 @@ TEST(AgradPartialsVari, OperandsAndPartialsMultivar) {
// XXX Test mixed - operands_and_partials<std::vector<matrix_v>, vector_d, vector_v>
TEST(AgradPartialsVari, OperandsAndPartialsMultivarMixed) {
- using stan::math::detail::operands_and_partials;
+ using stan::math::operands_and_partials;
using stan::math::var;
using stan::math::vector_v;
using stan::math::vector_d;
@@ -229,11 +264,11 @@ TEST(AgradPartialsVari, OperandsAndPartialsMultivarMixed) {
operands_and_partials<std::vector<vector_v>,
std::vector<vector_d>,
vector_v> o4(v_vec, d_vec_vec, v_vec2);
- o4.increment_dx1_vector(0, d_vec1);
- o4.increment_dx3_vector(0, d_vec2);
+ o4.edge1_.partials_vec[0] += d_vec1;
+ o4.edge3_.partials_vec[0] += d_vec2;
- // 2 partials stdvecs, 3 pointers to edges, 2 pointers to operands vecs
- EXPECT_EQ(2*sizeof(d_vec1) + 5*sizeof(&v_vec), sizeof(o4));
+ // 2 partials stdvecs, 4 pointers to edges, 2 pointers to operands vecs
+ EXPECT_EQ(2*sizeof(d_vec1) + 6*sizeof(&v_vec), sizeof(o4));
std::vector<double> grad;
var v = o4.build(10.0);
diff --git a/test/unit/math/rev/scal/meta/OperandsAndPartials_test.cpp b/test/unit/math/rev/scal/meta/OperandsAndPartials_test.cpp
deleted file mode 100644
index ef5345ad0818a706ff3a69f29c93f5d814da8266..0000000000000000000000000000000000000000
--- a/test/unit/math/rev/scal/meta/OperandsAndPartials_test.cpp
+++ /dev/null
@@ -1,73 +0,0 @@
-#include <stan/math/rev/scal.hpp>
-#include <gtest/gtest.h>
-
-TEST(AgradPartialsVari,OperandsAndPartials1) {
- using stan::math::vari;
- using stan::math::OperandsAndPartials;
- using stan::math::var;
-
- var x = 2.0;
- var z = -3.0 * x; // dz/dx = -3
- OperandsAndPartials<var> o(z);
- o.d_x1[0] += 5.0; // dy/dz = 5
-
- var y = o.value(-1.0);
-
- stan::math::grad(y.vi_);
- EXPECT_FLOAT_EQ(-15.0, x.adj()); // dy/dx = -15
-}
-TEST(AgradPartialsVari,OperandsAndPartials2) {
- using stan::math::OperandsAndPartials;
- using stan::math::vari;
- using stan::math::var;
-
- var x1 = 2.0;
- var x2 = 3.0;
- var z1 = -5.0 * x1; // dz1/dx1 = -5
- var z2 = -7.0 * x2; // dz2/dx2 = -7
- OperandsAndPartials<var,var> o(z1,z2);
- o.d_x1[0] += 11.0; // dy/dz1 = 11.0
- o.d_x2[0] += 13.0; // dy/dz2 = 13.0
- var y = o.value(-1.0);
-
- stan::math::grad(y.vi_);
- EXPECT_FLOAT_EQ(-55.0, x1.adj()); // dy/dx1 = -55
- EXPECT_FLOAT_EQ(-91.0, x2.adj()); // dy/dx2 = -91
-}
-TEST(AgradPartialsVari, OperandsAndPartials3) {
- using stan::math::OperandsAndPartials;
- using stan::math::vari;
- using stan::math::var;
-
- var x1 = 2.0;
- var x2 = 3.0;
- var x3 = 5.0;
- var z1 = -7.0 * x1; // dz1/dx1 = -5
- var z2 = -9.0 * x2; // dz2/dx2 = -7
- var z3 = -11.0 * x3; // dz3/dx3 = -11
-
- OperandsAndPartials<var,var,var> o(z1, z2, z3);
- o.d_x1[0] += 17.0; // dy/dz1 = 17.0
- o.d_x2[0] += 19.0; // dy/dz2 = 19.0
- o.d_x3[0] += 23.0; // dy/dz3 = 23.0
- var y = o.value(-1.0);
-
- stan::math::grad(y.vi_);
- EXPECT_FLOAT_EQ(-119.0, x1.adj()); // dy/dx1 = -119
- EXPECT_FLOAT_EQ(-171.0, x2.adj()); // dy/dx2 = -133
- EXPECT_FLOAT_EQ(-253.0, x3.adj()); // dy/dx2 = -253
-}
-TEST(AgradPartialsVari, OperandsAndPartials_check_throw) {
- using stan::math::OperandsAndPartials;
- using stan::math::var;
-
- var v;
-
- OperandsAndPartials<var,var,var,var,var,var> o2(v,v,v,v,v,v);
- EXPECT_NO_THROW(o2.d_x1[0]);
- EXPECT_NO_THROW(o2.d_x2[0]);
- EXPECT_NO_THROW(o2.d_x3[0]);
- EXPECT_NO_THROW(o2.d_x4[0]);
- EXPECT_NO_THROW(o2.d_x5[0]);
- EXPECT_NO_THROW(o2.d_x6[0]);
-}