Merge conflict

Author: andrjohns

stan/math/opencl/kernel_generator/select.hpp

@@ -3,6 +3,7 @@
 #ifdef STAN_OPENCL
 
 #include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/fun/select.hpp>
 #include <stan/math/opencl/matrix_cl_view.hpp>
 #include <stan/math/opencl/kernel_generator/type_str.hpp>
 #include <stan/math/opencl/kernel_generator/name_generator.hpp>
@@ -150,22 +151,6 @@ select(T_condition&& condition, T_then&& then, T_else&& els) {  // NOLINT
           as_operation_cl(std::forward<T_else>(els))};
 }
 
-/**
- * Scalar overload of the selection operation.
- * @tparam T_then type of then scalar
- * @tparam T_else type of else scalar
- * @param condition condition
- * @param then then result
- * @param els else result
- * @return `condition ? then : els`
- */
-template <typename T_then, typename T_else,
-          require_all_arithmetic_t<T_then, T_else>* = nullptr>
-inline std::common_type_t<T_then, T_else> select(bool condition, T_then then,
-                                                 T_else els) {
-  return condition ? then : els;
-}
-
 /** @}*/
 }  // namespace math
 }  // namespace stan

stan/math/prim/fun.hpp

@@ -7,6 +7,8 @@
 #include <stan/math/prim/fun/acosh.hpp>
 #include <stan/math/prim/fun/add.hpp>
 #include <stan/math/prim/fun/add_diag.hpp>
+#include <stan/math/prim/fun/all.hpp>
+#include <stan/math/prim/fun/any.hpp>
 #include <stan/math/prim/fun/append_array.hpp>
 #include <stan/math/prim/fun/append_col.hpp>
 #include <stan/math/prim/fun/append_row.hpp>
@@ -62,6 +64,7 @@
 #include <stan/math/prim/fun/cov_matrix_free.hpp>
 #include <stan/math/prim/fun/cov_matrix_free_lkj.hpp>
 #include <stan/math/prim/fun/crossprod.hpp>
+#include <stan/math/prim/fun/csr_extract.hpp>
 #include <stan/math/prim/fun/csr_extract_u.hpp>
 #include <stan/math/prim/fun/csr_extract_v.hpp>
 #include <stan/math/prim/fun/csr_extract_w.hpp>
@@ -82,6 +85,8 @@
 #include <stan/math/prim/fun/dot_product.hpp>
 #include <stan/math/prim/fun/dot_self.hpp>
 #include <stan/math/prim/fun/eigen_comparisons.hpp>
+#include <stan/math/prim/fun/eigendecompose.hpp>
+#include <stan/math/prim/fun/eigendecompose_sym.hpp>
 #include <stan/math/prim/fun/eigenvalues.hpp>
 #include <stan/math/prim/fun/eigenvalues_sym.hpp>
 #include <stan/math/prim/fun/eigenvectors.hpp>
@@ -275,6 +280,7 @@
 #include <stan/math/prim/fun/qr.hpp>
 #include <stan/math/prim/fun/qr_Q.hpp>
 #include <stan/math/prim/fun/qr_R.hpp>
+#include <stan/math/prim/fun/qr_thin.hpp>
 #include <stan/math/prim/fun/qr_thin_Q.hpp>
 #include <stan/math/prim/fun/qr_thin_R.hpp>
 #include <stan/math/prim/fun/quad_form.hpp>
@@ -306,6 +312,7 @@
 #include <stan/math/prim/fun/sd.hpp>
 #include <stan/math/prim/fun/segment.hpp>
 #include <stan/math/prim/fun/serializer.hpp>
+#include <stan/math/prim/fun/select.hpp>
 #include <stan/math/prim/fun/sign.hpp>
 #include <stan/math/prim/fun/signbit.hpp>
 #include <stan/math/prim/fun/simplex_constrain.hpp>
@@ -331,6 +338,7 @@
 #include <stan/math/prim/fun/sub_row.hpp>
 #include <stan/math/prim/fun/subtract.hpp>
 #include <stan/math/prim/fun/sum.hpp>
+#include <stan/math/prim/fun/svd.hpp>
 #include <stan/math/prim/fun/svd_U.hpp>
 #include <stan/math/prim/fun/svd_V.hpp>
 #include <stan/math/prim/fun/symmetrize_from_lower_tri.hpp>

stan/math/prim/fun/all.hpp

@@ -0,0 +1,85 @@
+#ifndef STAN_MATH_PRIM_FUN_ALL_HPP
+#define STAN_MATH_PRIM_FUN_ALL_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/functor/for_each.hpp>
+#include <algorithm>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return true if all values in the input are true.
+ *
+ * Overload for a single integral input
+ *
+ * @tparam T Any type convertible to `bool`
+ * @param x integral input
+ * @return The input unchanged
+ */
+template <typename T, require_t<std::is_convertible<T, bool>>* = nullptr>
+constexpr inline bool all(T x) {
+  return x;
+}
+
+/**
+ * Return true if all values in the input are true.
+ *
+ * Overload for Eigen types
+ *
+ * @tparam ContainerT A type derived from `Eigen::EigenBase` that has an
+ *                      `integral` scalar type
+ * @param x Eigen object of boolean inputs
+ * @return Boolean indicating whether all elements are true
+ */
+template <typename ContainerT,
+          require_eigen_st<std::is_integral, ContainerT>* = nullptr>
+inline bool all(const ContainerT& x) {
+  return x.all();
+}
+
+// Forward-declaration for correct resolution of all(std::vector<std::tuple>)
+template <typename... Types>
+inline bool all(const std::tuple<Types...>& x);
+
+/**
+ * Return true if all values in the input are true.
+ *
+ * Overload for a std::vector/nested inputs. The Eigen::Map/apply_vector_unary
+ * approach cannot be used as std::vector<bool> types do not have a .data()
+ * member and are not always stored contiguously.
+ *
+ * @tparam InnerT Type within std::vector
+ * @param x Nested container of boolean inputs
+ * @return Boolean indicating whether all elements are true
+ */
+template <typename InnerT>
+inline bool all(const std::vector<InnerT>& x) {
+  return std::all_of(x.begin(), x.end(), [](const auto& i) { return all(i); });
+}
+
+/**
+ * Return true if all values in the input are true.
+ *
+ * Overload for a tuple input.
+ *
+ * @tparam Types of items within tuple
+ * @param x Tuple of boolean scalar-type elements
+ * @return Boolean indicating whether all elements are true
+ */
+template <typename... Types>
+inline bool all(const std::tuple<Types...>& x) {
+  bool all_true = true;
+  math::for_each(
+      [&all_true](const auto& i) {
+        all_true = all_true && all(i);
+        return;
+      },
+      x);
+  return all_true;
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif

stan/math/prim/fun/any.hpp

@@ -0,0 +1,85 @@
+#ifndef STAN_MATH_PRIM_FUN_ANY_HPP
+#define STAN_MATH_PRIM_FUN_ANY_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/functor/for_each.hpp>
+#include <algorithm>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return true if any values in the input are true.
+ *
+ * Overload for a single boolean input
+ *
+ * @tparam T Any type convertible to `bool`
+ * @param x boolean input
+ * @return The input unchanged
+ */
+template <typename T, require_t<std::is_convertible<T, bool>>* = nullptr>
+constexpr inline bool any(T x) {
+  return x;
+}
+
+/**
+ * Return true if any values in the input are true.
+ *
+ * Overload for Eigen types
+ *
+ * @tparam ContainerT A type derived from `Eigen::EigenBase` that has an
+ *                      `integral` scalar type
+ * @param x Eigen object of boolean inputs
+ * @return Boolean indicating whether any elements are true
+ */
+template <typename ContainerT,
+          require_eigen_st<std::is_integral, ContainerT>* = nullptr>
+inline bool any(const ContainerT& x) {
+  return x.any();
+}
+
+// Forward-declaration for correct resolution of any(std::vector<std::tuple>)
+template <typename... Types>
+inline bool any(const std::tuple<Types...>& x);
+
+/**
+ * Return true if any values in the input are true.
+ *
+ * Overload for a std::vector/nested inputs. The Eigen::Map/apply_vector_unary
+ * approach cannot be used as std::vector<bool> types do not have a .data()
+ * member and are not always stored contiguously.
+ *
+ * @tparam InnerT Type within std::vector
+ * @param x Nested container of boolean inputs
+ * @return Boolean indicating whether any elements are true
+ */
+template <typename InnerT>
+inline bool any(const std::vector<InnerT>& x) {
+  return std::any_of(x.begin(), x.end(), [](const auto& i) { return any(i); });
+}
+
+/**
+ * Return true if any values in the input are true.
+ *
+ * Overload for a tuple input.
+ *
+ * @tparam Types of items within tuple
+ * @param x Tuple of boolean scalar-type elements
+ * @return Boolean indicating whether any elements are true
+ */
+template <typename... Types>
+inline bool any(const std::tuple<Types...>& x) {
+  bool any_true = false;
+  math::for_each(
+      [&any_true](const auto& i) {
+        any_true = any_true || any(i);
+        return;
+      },
+      x);
+  return any_true;
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif

stan/math/prim/fun/complex_schur_decompose.hpp

@@ -27,13 +27,13 @@ namespace math {
 template <typename M, require_eigen_dense_dynamic_t<M>* = nullptr>
 inline Eigen::Matrix<complex_return_t<scalar_type_t<M>>, -1, -1>
 complex_schur_decompose_u(const M& m) {
-  if (m.size() == 0)
+  if (unlikely(m.size() == 0)) {
     return m;
+  }
   check_square("complex_schur_decompose_u", "m", m);
   using MatType = Eigen::Matrix<scalar_type_t<M>, -1, -1>;
   // copy because ComplexSchur requires Eigen::Matrix type
-  MatType mv = m;
-  Eigen::ComplexSchur<MatType> cs(mv);
+  Eigen::ComplexSchur<MatType> cs{MatType(m)};
   return cs.matrixU();
 }
 
@@ -51,16 +51,46 @@ complex_schur_decompose_u(const M& m) {
 template <typename M, require_eigen_dense_dynamic_t<M>* = nullptr>
 inline Eigen::Matrix<complex_return_t<scalar_type_t<M>>, -1, -1>
 complex_schur_decompose_t(const M& m) {
-  if (m.size() == 0)
+  if (unlikely(m.size() == 0)) {
     return m;
+  }
   check_square("complex_schur_decompose_t", "m", m);
   using MatType = Eigen::Matrix<scalar_type_t<M>, -1, -1>;
   // copy because ComplexSchur requires Eigen::Matrix type
-  MatType mv = m;
-  Eigen::ComplexSchur<MatType> cs(mv, false);
+  Eigen::ComplexSchur<MatType> cs{MatType(m), false};
   return cs.matrixT();
 }
 
+/**
+ * Return the complex Schur decomposition of the
+ * specified square matrix.
+ *
+ * The complex Schur decomposition of a square matrix `A` produces a
+ * complex unitary matrix `U` and a complex upper-triangular Schur
+ * form matrix `T` such that `A = U * T * inv(U)`.  Further, the
+ * unitary matrix's inverse is equal to its conjugate transpose,
+ * `inv(U) = U*`, where `U*(i, j) = conj(U(j, i))`
+ *
+ * @tparam M type of matrix
+ * @param m real matrix to decompose
+ * @return a tuple (U,T) where U is the complex unitary matrix of the complex
+ * Schur decomposition of `m` and T is the Schur form matrix of
+ * the complex Schur decomposition of `m`
+ */
+template <typename M, require_eigen_dense_dynamic_t<M>* = nullptr>
+inline std::tuple<Eigen::Matrix<complex_return_t<scalar_type_t<M>>, -1, -1>,
+                  Eigen::Matrix<complex_return_t<scalar_type_t<M>>, -1, -1>>
+complex_schur_decompose(const M& m) {
+  if (unlikely(m.size() == 0)) {
+    return std::make_tuple(m, m);
+  }
+  check_square("complex_schur_decompose", "m", m);
+  using MatType = Eigen::Matrix<scalar_type_t<M>, -1, -1>;
+  // copy because ComplexSchur requires Eigen::Matrix type
+  Eigen::ComplexSchur<MatType> cs{MatType(m)};
+  return std::make_tuple(std::move(cs.matrixU()), std::move(cs.matrixT()));
+}
+
 }  // namespace math
 }  // namespace stan
 #endif

stan/math/prim/fun/csr_extract.hpp

@@ -0,0 +1,66 @@
+#ifndef STAN_MATH_PRIM_FUN_CSR_EXTRACT_HPP
+#define STAN_MATH_PRIM_FUN_CSR_EXTRACT_HPP
+
+#include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
+
+namespace stan {
+namespace math {
+
+/** \addtogroup csr_format
+ *  @{
+ */
+
+/**
+ * Extract the non-zero values, column indexes for non-zero values, and
+ * the NZE index for each entry from a sparse matrix.
+ *
+ * @tparam T type of elements in the matrix
+ * @param[in] A sparse matrix.
+ * @return a tuple W,V,U.
+ */
+template <typename T>
+const std::tuple<Eigen::Matrix<T, Eigen::Dynamic, 1>, std::vector<int>,
+                 std::vector<int>>
+csr_extract(const Eigen::SparseMatrix<T, Eigen::RowMajor>& A) {
+  auto a_nonzeros = A.nonZeros();
+  Eigen::Matrix<T, Eigen::Dynamic, 1> w
+      = Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(a_nonzeros);
+  std::vector<int> v(a_nonzeros);
+  for (int nze = 0; nze < a_nonzeros; ++nze) {
+    w[nze] = *(A.valuePtr() + nze);
+    v[nze] = *(A.innerIndexPtr() + nze) + stan::error_index::value;
+  }
+  std::vector<int> u(A.outerSize() + 1);  // last entry is garbage.
+  for (int nze = 0; nze <= A.outerSize(); ++nze) {
+    u[nze] = *(A.outerIndexPtr() + nze) + stan::error_index::value;
+  }
+  return std::make_tuple(std::move(w), std::move(v), std::move(u));
+}
+
+/* Extract the non-zero values from a dense matrix by converting
+ * to sparse and calling the sparse matrix extractor.
+ *
+ * @tparam T type of elements in the matrix
+ * @tparam R number of rows, can be Eigen::Dynamic
+ * @tparam C number of columns, can be Eigen::Dynamic
+ *
+ * @param[in] A dense matrix.
+ * @return a tuple W,V,U.
+ */
+template <typename T, require_eigen_dense_base_t<T>* = nullptr>
+const std::tuple<Eigen::Matrix<scalar_type_t<T>, Eigen::Dynamic, 1>,
+                 std::vector<int>, std::vector<int>>
+csr_extract(const T& A) {
+  // conversion to sparse seems to touch data twice, so we need to call to_ref
+  Eigen::SparseMatrix<scalar_type_t<T>, Eigen::RowMajor> B
+      = to_ref(A).sparseView();
+  return csr_extract(B);
+}
+
+/** @} */  // end of csr_format group
+
+}  // namespace math
+}  // namespace stan
+
+#endif