Serialization for sum_to_zero

src/stan/io/deserializer.hpp

@@ -581,6 +581,64 @@ class deserializer {
     return ret;
   }
 
+  /**
+   * Return the next zero sum vector of the specified size (using one fewer
+   * unconstrained scalars), incrementing the specified reference with the
+   * log absolute Jacobian determinant (no adjustment, in this case).
+   *
+   * <p>See <code>stan::math::sum_to_zero_constrain(Eigen::Matrix,T&)</code>.
+   *
+   * @tparam Ret The type to return.
+   * @tparam Jacobian Whether to increment the log of the absolute Jacobian
+   * determinant of the transform.
+   * @tparam LP Type of log probability.
+   * @tparam Sizes A parameter pack of integral types.
+   * @param lp The reference to the variable holding the log
+   * @param size Number of cells in zero sum vector to generate.
+   * @return The next zero sum of the specified size.
+   * @throws std::invalid_argument if number of dimensions (`k`) is zero
+   */
+  template <typename Ret, bool Jacobian, typename LP,
+            require_not_std_vector_t<Ret>* = nullptr>
+  inline auto read_constrain_sum_to_zero(LP& lp, size_t size) {
+    stan::math::check_positive("read_sum_to_zero", "size", size);
+    return stan::math::sum_to_zero_constrain<Jacobian>(
+        this->read<Ret>(size - 1), lp);
+  }
+
+  /**
+   * Return the next zero sum vector of the specified size (using one fewer
+   * unconstrained scalars), incrementing the specified reference with the
+   * log absolute Jacobian determinant (no adjustment, in this case).
+   *
+   * <p>See <code>stan::math::sum_to_zero_constrain(Eigen::Matrix,T&)</code>.
+   *
+   * @tparam Ret The type to return.
+   * @tparam Jacobian Whether to increment the log of the absolute Jacobian
+   * determinant of the transform.
+   * @tparam LP Type of log probability.
+   * @tparam Sizes A parameter pack of integral types.
+   * @param lp The reference to the variable holding the log
+   * probability to increment.
+   * @param vecsize The size of the return vector.
+   * @param sizes Pack of integrals to use to construct the return's type.
+   * @return The next zero sum of the specified size.
+   * @throws std::invalid_argument if number of dimensions (`k`) is zero
+   */
+  template <typename Ret, bool Jacobian, typename LP, typename... Sizes,
+            require_std_vector_t<Ret>* = nullptr>
+  inline auto read_constrain_sum_to_zero(LP& lp, const size_t vecsize,
+                                         Sizes... sizes) {
+    std::decay_t<Ret> ret;
+    ret.reserve(vecsize);
+    for (size_t i = 0; i < vecsize; ++i) {
+      ret.emplace_back(
+          this->read_constrain_sum_to_zero<value_type_t<Ret>, Jacobian>(
+              lp, sizes...));
+    }
+    return ret;
+  }
+
   /**
    * Return the next ordered vector of the specified
    * size, incrementing the specified reference with the log
@@ -933,7 +991,7 @@ class deserializer {
    * @param lp The reference to the variable holding the log
    * probability to increment.
    * @param rows Rows of matrix
-   * @param cows Cows of matrix
+   * @param cols Cols of matrix
    */
   template <typename Ret, bool Jacobian, typename LP,
             require_not_std_vector_t<Ret>* = nullptr,
@@ -989,7 +1047,7 @@ class deserializer {
    * @param lp The reference to the variable holding the log
    * probability to increment.
    * @param rows Rows of matrix
-   * @param cows Cows of matrix
+   * @param cols Cols of matrix
    */
   template <typename Ret, bool Jacobian, typename LP,
             require_not_std_vector_t<Ret>* = nullptr,
@@ -1160,6 +1218,37 @@ class deserializer {
     return ret;
   }
 
+  /**
+   * Read a serialized sum_to_zero vector and unconstrain it
+   *
+   * @tparam Ret Type of output
+   * @return Unconstrained vector
+   */
+  template <typename Ret, require_not_std_vector_t<Ret>* = nullptr>
+  inline auto read_free_sum_to_zero(size_t size) {
+    return stan::math::sum_to_zero_free(this->read<Ret>(size));
+  }
+
+  /**
+   * Read serialized zero-sum vectors and unconstrain them
+   *
+   * @tparam Ret Type of output
+   * @tparam Sizes Types of dimensions of output
+   * @param vecsize Vector size
+   * @param sizes dimensions
+   * @return Unconstrained vectors
+   */
+  template <typename Ret, typename... Sizes,
+            require_std_vector_t<Ret>* = nullptr>
+  inline auto read_free_sum_to_zero(size_t vecsize, Sizes... sizes) {
+    std::decay_t<Ret> ret;
+    ret.reserve(vecsize);
+    for (size_t i = 0; i < vecsize; ++i) {
+      ret.emplace_back(read_free_sum_to_zero<value_type_t<Ret>>(sizes...));
+    }
+    return ret;
+  }
+
   /**
    * Read a serialized ordered and unconstrain it
    *
@@ -1358,7 +1447,7 @@ class deserializer {
    *
    * @tparam Ret Type of output
    * @param rows Rows of matrix
-   * @param cows Cows of matrix
+   * @param cols Cols of matrix
    * @return Unconstrained matrix
    */
   template <typename Ret, require_not_std_vector_t<Ret>* = nullptr>
@@ -1392,7 +1481,7 @@ class deserializer {
    *
    * @tparam Ret Type of output
    * @param rows Rows of matrix
-   * @param cows Cows of matrix
+   * @param cols Cols of matrix
    * @return Unconstrained matrix
    */
   template <typename Ret, require_not_std_vector_t<Ret>* = nullptr>

src/stan/io/serializer.hpp

@@ -347,6 +347,31 @@ class serializer {
     }
   }
 
+  /**
+   * Write a serialized sum-to-zero vector and unconstrain it
+   *
+   * @tparam Vec An Eigen type with either fixed rows or columns at compile
+   * time.
+   * @param x The vector to read from.
+   */
+  template <typename Vec, require_not_std_vector_t<Vec>* = nullptr>
+  inline void write_free_sum_to_zero(const Vec& x) {
+    this->write(stan::math::sum_to_zero_free(x));
+  }
+
+  /**
+   * Write serialized zero sum vectors and unconstrain them
+   *
+   * @tparam StdVec A `std:vector`
+   * @param x An std vector.
+   */
+  template <typename StdVec, require_std_vector_t<StdVec>* = nullptr>
+  inline void write_free_sum_to_zero(const StdVec& x) {
+    for (const auto& ret_i : x) {
+      this->write_free_sum_to_zero(ret_i);
+    }
+  }
+
   /**
    * Write a serialized ordered and unconstrain it
    *

src/test/unit/io/deserializer_free_test.cpp

@@ -248,6 +248,30 @@ TEST(deserializer_vector, read_free_simplex) {
                     SimplexConstrain>(std::make_tuple(3, 2, 4));
 }
 
+// sum_to_zero
+template <typename Ret>
+struct SumToZeroConstrain {
+  static auto read() {
+    return [](stan::io::deserializer<double>& deserializer, auto&&... args) {
+      return deserializer.read_constrain_sum_to_zero<Ret, false>(args...);
+    };
+  }
+  static auto free() {
+    return [](stan::io::deserializer<double>& deserializer, auto&&... args) {
+      return deserializer.read_free_sum_to_zero<Ret>(args...);
+    };
+  }
+};
+
+TEST(deserializer_vector, read_free_sum_to_zero) {
+  using stan::test::deserializer_test;
+  deserializer_test<Eigen::VectorXd, SumToZeroConstrain>(std::make_tuple(4));
+  deserializer_test<std::vector<Eigen::VectorXd>, SumToZeroConstrain>(
+      std::make_tuple(2, 4));
+  deserializer_test<std::vector<std::vector<Eigen::VectorXd>>,
+                    SumToZeroConstrain>(std::make_tuple(3, 2, 4));
+}
+
 // ordered
 template <typename Ret>
 struct OrderedConstrain {

src/test/unit/io/deserializer_stdvector_test.cpp

@@ -129,6 +129,52 @@ TEST(deserializer_array, simplex) {
   }
 }
 
+// sum_to_zero
+
+TEST(deserializer_array, sum_to_zero) {
+  std::vector<int> theta_i;
+  std::vector<double> theta;
+  for (size_t i = 0; i < 100U; ++i)
+    theta.push_back(static_cast<double>(i));
+
+  stan::io::deserializer<double> deserializer1(theta, theta_i);
+  stan::io::deserializer<double> deserializer2(theta, theta_i);
+
+  // no jac
+  {
+    double lp_ref = 0.0;
+    double lp = 0.0;
+    auto y
+        = deserializer1
+              .read_constrain_sum_to_zero<std::vector<Eigen::VectorXd>, false>(
+                  lp, 4, 3);
+    for (size_t i = 0; i < 4; ++i) {
+      stan::test::expect_near_rel(
+          "test_std_vector_deserializer", y[i],
+          deserializer2.read_constrain_sum_to_zero<Eigen::VectorXd, false>(
+              lp_ref, 3));
+    }
+    EXPECT_FLOAT_EQ(lp_ref, lp);
+  }
+
+  // jac
+  {
+    double lp_ref = 0.0;
+    double lp = 0.0;
+    auto y
+        = deserializer1
+              .read_constrain_sum_to_zero<std::vector<Eigen::VectorXd>, true>(
+                  lp, 4, 3);
+    for (size_t i = 0; i < 4; ++i) {
+      stan::test::expect_near_rel(
+          "test_std_vector_deserializer", y[i],
+          deserializer2.read_constrain_sum_to_zero<Eigen::VectorXd, true>(
+              lp_ref, 3));
+    }
+    EXPECT_FLOAT_EQ(lp_ref, lp);
+  }
+}
+
 // ordered
 
 TEST(deserializer_array, ordered) {

src/test/unit/io/deserializer_test.cpp

@@ -672,6 +672,48 @@ TEST(deserializer_vector, simplex_jacobian) {
   EXPECT_FLOAT_EQ(lp_ref, lp);
 }
 
+// sum_to_zero
+
+TEST(deserializer_vector, sum_to_zero_constrain) {
+  std::vector<int> theta_i;
+  std::vector<double> theta;
+  theta.push_back(3.0);
+  theta.push_back(-1.0);
+  theta.push_back(-2.0);
+  theta.push_back(0.0);
+  stan::io::deserializer<double> deserializer(theta, theta_i);
+  double lp = 0;
+  Eigen::VectorXd reference
+      = stan::math::sum_to_zero_constrain(stan::math::to_vector(theta));
+  Eigen::VectorXd phi(
+      deserializer.read_constrain_sum_to_zero<Eigen::VectorXd, false>(
+          lp, theta.size() + 1));
+  for (size_t i = 0; i < phi.size(); ++i) {
+    EXPECT_FLOAT_EQ(reference(i), phi[i]);
+  }
+}
+
+TEST(deserializer_vector, sum_to_zero_jacobian) {
+  std::vector<int> theta_i;
+  std::vector<double> theta;
+  theta.push_back(3.0);
+  theta.push_back(-1.0);
+  theta.push_back(-2.0);
+  theta.push_back(0.0);
+  stan::io::deserializer<double> deserializer(theta, theta_i);
+  double lp = 0.0;
+  double lp_ref = 0.0;
+  Eigen::VectorXd reference
+      = stan::math::sum_to_zero_constrain(stan::math::to_vector(theta), lp_ref);
+  Eigen::VectorXd phi(
+      deserializer.read_constrain_sum_to_zero<Eigen::VectorXd, true>(
+          lp, theta.size() + 1));
+  for (size_t i = 0; i < phi.size(); ++i) {
+    EXPECT_FLOAT_EQ(reference(i), phi[i]);
+  }
+  EXPECT_FLOAT_EQ(lp_ref, lp);
+}
+
 // ordered
 
 TEST(deserializer_vector, ordered_constrain) {

src/test/unit/io/deserializer_var_test.cpp

@@ -597,6 +597,49 @@ TEST(deserializer_vector, simplex_jacobian) {
   EXPECT_FLOAT_EQ(lp_ref.val(), lp.val());
 }
 
+// sum_to_zero
+
+TEST(deserializer_vector, sum_to_zero_constrain) {
+  std::vector<int> theta_i;
+  std::vector<stan::math::var> theta;
+  theta.push_back(3.0);
+  theta.push_back(-1.0);
+  theta.push_back(-2.0);
+  theta.push_back(0.0);
+  stan::io::deserializer<stan::math::var> deserializer(theta, theta_i);
+  stan::math::var lp = 0;
+  Eigen::Matrix<stan::math::var, Eigen::Dynamic, 1> reference
+      = stan::math::sum_to_zero_constrain(stan::math::to_vector(theta));
+  Eigen::Matrix<stan::math::var, Eigen::Dynamic, 1> phi(
+      deserializer.read_constrain_sum_to_zero<
+          Eigen::Matrix<stan::math::var, Eigen::Dynamic, 1>, false>(
+          lp, theta.size() + 1));
+  for (size_t i = 0; i < phi.size(); ++i) {
+    EXPECT_FLOAT_EQ(reference(i).val(), phi[i].val());
+  }
+}
+
+TEST(deserializer_vector, sum_to_zero_jacobian) {
+  std::vector<int> theta_i;
+  std::vector<stan::math::var> theta;
+  theta.push_back(3.0);
+  theta.push_back(-1.0);
+  theta.push_back(-2.0);
+  theta.push_back(0.0);
+  stan::io::deserializer<stan::math::var> deserializer(theta, theta_i);
+  stan::math::var lp = 0.0;
+  stan::math::var lp_ref = 0.0;
+  Eigen::Matrix<stan::math::var, Eigen::Dynamic, 1> reference
+      = stan::math::sum_to_zero_constrain(stan::math::to_vector(theta), lp_ref);
+  Eigen::Matrix<stan::math::var, Eigen::Dynamic, 1> phi(
+      deserializer.read_constrain_sum_to_zero<
+          Eigen::Matrix<stan::math::var, Eigen::Dynamic, 1>, true>(
+          lp, theta.size() + 1));
+  for (size_t i = 0; i < phi.size(); ++i) {
+    EXPECT_FLOAT_EQ(reference(i).val(), phi[i].val());
+  }
+  EXPECT_FLOAT_EQ(lp_ref.val(), lp.val());
+}
 // ordered
 
 TEST(deserializer_vector, ordered_constrain) {

src/test/unit/io/deserializer_varmat_test.cpp

@@ -325,6 +325,44 @@ TEST(deserializer, read_constrain_simplex_jacobian) {
   EXPECT_FLOAT_EQ(lp_ref.val(), lp.val());
 }
 
+// sum_to_zero
+
+TEST(deserializer, read_constrain_sum_to_zero_constrain) {
+  std::vector<int> theta_i;
+  std::vector<stan::math::var> theta;
+  theta.push_back(-2.0);
+  theta.push_back(3.0);
+  theta.push_back(-1.0);
+  theta.push_back(0.0);
+  stan::io::deserializer<stan::math::var> deserializer(theta, theta_i);
+  stan::math::var lp = 0.0;
+  auto reference
+      = stan::math::sum_to_zero_constrain(stan::math::to_vector(theta));
+  auto y = deserializer.read_constrain_sum_to_zero<var_vector_t, false>(
+      lp, theta.size() + 1);
+  EXPECT_TRUE((std::is_same<var_vector_t, decltype(y)>::value));
+  stan::test::expect_near_rel("deserializer tests", reference.val(), y.val());
+}
+
+TEST(deserializer, read_constrain_sum_to_zero_jacobian) {
+  std::vector<int> theta_i;
+  std::vector<stan::math::var> theta;
+  theta.push_back(-2.0);
+  theta.push_back(3.0);
+  theta.push_back(-1.0);
+  theta.push_back(0.0);
+  stan::io::deserializer<stan::math::var> deserializer(theta, theta_i);
+  stan::math::var lp_ref = 0.0;
+  stan::math::var lp = 0.0;
+  auto reference
+      = stan::math::sum_to_zero_constrain(stan::math::to_vector(theta), lp_ref);
+  auto y = deserializer.read_constrain_sum_to_zero<var_vector_t, true>(
+      lp, theta.size() + 1);
+  EXPECT_TRUE((std::is_same<var_vector_t, decltype(y)>::value));
+  stan::test::expect_near_rel("deserializer tests", reference.val(), y.val());
+  EXPECT_FLOAT_EQ(lp_ref.val(), lp.val());
+}
+
 // ordered
 
 TEST(deserializer, read_constrain_ordered_constrain) {