Skip to content

Commit 2d99d55

Browse files
bbbales2bbbales2
committed
Merge remote-tracking branch 'origin/develop' into feature/a_few_varmats
2 parents ec4b37e + dcc2bb1 commit 2d99d55

32 files changed

+438
-2275
lines changed

‎make/tests

+9
Original file line numberDiff line numberDiff line change
@@ -35,6 +35,15 @@ test/unit/libmultiple.so : test/unit/multiple_translation_units1.o test/unit/mul
3535

3636
test/unit/multiple_translation_units_test$(EXE): test/unit/libmultiple.so
3737

38+
############################################################
39+
#
40+
# Expression tests
41+
##
42+
43+
EXPRESSION_TESTS := $(subst .cpp,$(EXE),$(call findfiles,test/expressions,*_test.cpp))
44+
$(EXPRESSION_TESTS) : $(LIBSUNDIALS)
45+
46+
3847
############################################################
3948
#
4049
# CVODES tests

‎stan/math/fwd/fun/multiply_lower_tri_self_transpose.hpp

+1-7
Original file line numberDiff line numberDiff line change
@@ -16,14 +16,8 @@ multiply_lower_tri_self_transpose(const EigMat& m) {
1616
if (m.rows() == 0) {
1717
return {};
1818
}
19-
plain_type_t<EigMat> L(m.rows(), m.cols());
20-
L.setZero();
19+
plain_type_t<EigMat> L = m.template triangularView<Eigen::Lower>();
2120

22-
for (size_type i = 0; i < m.rows(); i++) {
23-
for (size_type j = 0; (j < i + 1) && (j < m.cols()); j++) {
24-
L.coeffRef(i, j) = m.coeff(i, j);
25-
}
26-
}
2721
return multiply(L, L.transpose());
2822
}
2923

‎stan/math/prim/fun/block.hpp

+3-4
Original file line numberDiff line numberDiff line change
@@ -18,10 +18,9 @@ namespace math {
1818
* @param ncols Number of columns in block.
1919
* @throw std::out_of_range if either index is out of range.
2020
*/
21-
template <typename T>
22-
inline Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> block(
23-
const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& m, size_t i,
24-
size_t j, size_t nrows, size_t ncols) {
21+
template <typename T, require_eigen_t<T>* = nullptr>
22+
inline Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic> block(
23+
const T& m, size_t i, size_t j, size_t nrows, size_t ncols) {
2524
check_row_index("block", "i", m, i);
2625
check_row_index("block", "i+nrows-1", m, i + nrows - 1);
2726
check_column_index("block", "j", m, j);

‎stan/math/prim/fun/csr_extract_u.hpp

+6-3
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33

44
#include <stan/math/prim/meta.hpp>
55
#include <stan/math/prim/fun/Eigen.hpp>
6+
#include <stan/math/prim/fun/to_ref.hpp>
67
#include <vector>
78
#include <numeric>
89

@@ -39,9 +40,11 @@ const std::vector<int> csr_extract_u(
3940
* @param A Dense matrix.
4041
* @return Array of indexes into non-zero entries of A.
4142
*/
42-
template <typename T, int R, int C>
43-
const std::vector<int> csr_extract_u(const Eigen::Matrix<T, R, C>& A) {
44-
Eigen::SparseMatrix<T, Eigen::RowMajor> B = A.sparseView();
43+
template <typename T, require_eigen_dense_base_t<T>* = nullptr>
44+
const std::vector<int> csr_extract_u(const T& A) {
45+
// conversion to sparse seems to touch data twice, so we need to call to_ref
46+
Eigen::SparseMatrix<scalar_type_t<T>, Eigen::RowMajor> B
47+
= to_ref(A).sparseView();
4548
std::vector<int> u = csr_extract_u(B);
4649
return u;
4750
}

‎stan/math/prim/fun/csr_extract_v.hpp

+6-3
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33

44
#include <stan/math/prim/meta.hpp>
55
#include <stan/math/prim/fun/Eigen.hpp>
6+
#include <stan/math/prim/fun/to_ref.hpp>
67
#include <vector>
78
#include <numeric>
89

@@ -43,9 +44,11 @@ const std::vector<int> csr_extract_v(
4344
* @param[in] A dense matrix.
4445
* @return Array of column indexes to non-zero entries of A.
4546
*/
46-
template <typename T, int R, int C>
47-
const std::vector<int> csr_extract_v(const Eigen::Matrix<T, R, C>& A) {
48-
Eigen::SparseMatrix<T, Eigen::RowMajor> B = A.sparseView();
47+
template <typename T, require_eigen_dense_base_t<T>* = nullptr>
48+
const std::vector<int> csr_extract_v(const T& A) {
49+
// conversion to sparse seems to touch data twice, so we need to call to_ref
50+
Eigen::SparseMatrix<scalar_type_t<T>, Eigen::RowMajor> B
51+
= to_ref(A).sparseView();
4952
std::vector<int> v = csr_extract_v(B);
5053
return v;
5154
}

‎stan/math/prim/fun/csr_extract_w.hpp

+7-4
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@
22
#define STAN_MATH_PRIM_FUN_CSR_EXTRACT_W_HPP
33

44
#include <stan/math/prim/fun/Eigen.hpp>
5+
#include <stan/math/prim/fun/to_ref.hpp>
56

67
namespace stan {
78
namespace math {
@@ -37,10 +38,12 @@ const Eigen::Matrix<T, Eigen::Dynamic, 1> csr_extract_w(
3738
* @param[in] A dense matrix.
3839
* @return Vector of non-zero entries of A.
3940
*/
40-
template <typename T, int R, int C>
41-
const Eigen::Matrix<T, Eigen::Dynamic, 1> csr_extract_w(
42-
const Eigen::Matrix<T, R, C>& A) {
43-
Eigen::SparseMatrix<T, Eigen::RowMajor> B = A.sparseView();
41+
template <typename T, require_eigen_dense_base_t<T>* = nullptr>
42+
const Eigen::Matrix<scalar_type_t<T>, Eigen::Dynamic, 1> csr_extract_w(
43+
const T& A) {
44+
// conversion to sparse seems to touch data twice, so we need to call to_ref
45+
Eigen::SparseMatrix<scalar_type_t<T>, Eigen::RowMajor> B
46+
= to_ref(A).sparseView();
4447
return csr_extract_w(B);
4548
}
4649

‎stan/math/prim/fun/csr_matrix_times_vector.hpp

+6-8
Original file line numberDiff line numberDiff line change
@@ -49,8 +49,8 @@ namespace math {
4949
* column index of each value), the integer array u (containing
5050
* one-based indexes of where each row starts in w).
5151
*
52-
* @tparam T1 type of elements in the sparse matrix
53-
* @tparam T2 type of elements in the dense vector
52+
* @tparam T1 type of the sparse matrix
53+
* @tparam T2 type of the dense vector
5454
* @param m Number of rows in matrix.
5555
* @param n Number of columns in matrix.
5656
* @param w Vector of non-zero values in matrix.
@@ -71,10 +71,8 @@ namespace math {
7171
*/
7272
template <typename T1, typename T2>
7373
inline Eigen::Matrix<return_type_t<T1, T2>, Eigen::Dynamic, 1>
74-
csr_matrix_times_vector(int m, int n,
75-
const Eigen::Matrix<T1, Eigen::Dynamic, 1>& w,
76-
const std::vector<int>& v, const std::vector<int>& u,
77-
const Eigen::Matrix<T2, Eigen::Dynamic, 1>& b) {
74+
csr_matrix_times_vector(int m, int n, const T1& w, const std::vector<int>& v,
75+
const std::vector<int>& u, const T2& b) {
7876
using result_t = return_type_t<T1, T2>;
7977

8078
check_positive("csr_matrix_times_vector", "m", m);
@@ -99,9 +97,9 @@ csr_matrix_times_vector(int m, int n,
9997
for (int nze = u[row] - stan::error_index::value; nze < row_end_in_w;
10098
++nze, ++i) {
10199
check_range("csr_matrix_times_vector", "j", n, v[nze]);
102-
b_sub.coeffRef(i) = b.coeffRef(v[nze] - stan::error_index::value);
100+
b_sub.coeffRef(i) = b.coeff(v[nze] - stan::error_index::value);
103101
}
104-
Eigen::Matrix<T1, Eigen::Dynamic, 1> w_sub(
102+
Eigen::Matrix<value_type_t<T1>, Eigen::Dynamic, 1> w_sub(
105103
w.segment(u[row] - stan::error_index::value, idx));
106104
result.coeffRef(row) = dot_product(w_sub, b_sub);
107105
}

‎stan/math/prim/fun/csr_to_dense_matrix.hpp

+8-7
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44
#include <stan/math/prim/err.hpp>
55
#include <stan/math/prim/fun/csr_u_to_z.hpp>
66
#include <stan/math/prim/fun/Eigen.hpp>
7+
#include <stan/math/prim/fun/to_ref.hpp>
78
#include <vector>
89

910
namespace stan {
@@ -15,7 +16,7 @@ namespace math {
1516
/**
1617
* Construct a dense Eigen matrix from the CSR format components.
1718
*
18-
* @tparam T type of elements in the matrix
19+
* @tparam T type of the matrix
1920
* @param[in] m Number of matrix rows.
2021
* @param[in] n Number of matrix columns.
2122
* @param[in] w Values of non-zero matrix entries.
@@ -30,9 +31,9 @@ namespace math {
3031
* @throw std::out_of_range if any of the indices are out of range.
3132
*/
3233
template <typename T>
33-
inline Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> csr_to_dense_matrix(
34-
int m, int n, const Eigen::Matrix<T, Eigen::Dynamic, 1>& w,
35-
const std::vector<int>& v, const std::vector<int>& u) {
34+
inline Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic>
35+
csr_to_dense_matrix(int m, int n, const T& w, const std::vector<int>& v,
36+
const std::vector<int>& u) {
3637
using Eigen::Dynamic;
3738
using Eigen::Matrix;
3839

@@ -45,8 +46,8 @@ inline Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> csr_to_dense_matrix(
4546
for (int i : v) {
4647
check_range("csr_to_dense_matrix", "v[]", n, i);
4748
}
48-
49-
Matrix<T, Dynamic, Dynamic> result(m, n);
49+
const auto& w_ref = to_ref(w);
50+
Matrix<value_type_t<T>, Dynamic, Dynamic> result(m, n);
5051
result.setZero();
5152
for (int row = 0; row < m; ++row) {
5253
int row_end_in_w = (u[row] - stan::error_index::value) + csr_u_to_z(u, row);
@@ -55,7 +56,7 @@ inline Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> csr_to_dense_matrix(
5556
++nze) {
5657
// row is row index, v[nze] is column index. w[nze] is entry value.
5758
check_range("csr_to_dense_matrix", "j", n, v[nze]);
58-
result(row, v[nze] - stan::error_index::value) = w(nze);
59+
result(row, v[nze] - stan::error_index::value) = w_ref.coeff(nze);
5960
}
6061
}
6162
return result;

‎stan/math/prim/functor/map_rect.hpp

+15-6
Original file line numberDiff line numberDiff line change
@@ -94,6 +94,9 @@ namespace math {
9494
* const std::vector<double>& x_r, const std::vector<int>& x_i,
9595
* std::ostream* msgs = 0) const { ... }
9696
*
97+
* If an expression is passed as `shared_params`, the functor needs to accept
98+
* expression.
99+
*
97100
* WARNING: For the MPI case, the data arguments are NOT checked if
98101
* they are unchanged between repeated evaluations for a given
99102
* call_id/functor F pair. This is silently assumed to be immutable
@@ -118,9 +121,10 @@ namespace math {
118121
*/
119122

120123
template <int call_id, typename F, typename T_shared_param,
121-
typename T_job_param>
124+
typename T_job_param,
125+
require_eigen_col_vector_t<T_shared_param>* = nullptr>
122126
Eigen::Matrix<return_type_t<T_shared_param, T_job_param>, Eigen::Dynamic, 1>
123-
map_rect(const Eigen::Matrix<T_shared_param, Eigen::Dynamic, 1>& shared_params,
127+
map_rect(const T_shared_param& shared_params,
124128
const std::vector<Eigen::Matrix<T_job_param, Eigen::Dynamic, 1>>&
125129
job_params,
126130
const std::vector<std::vector<double>>& x_r,
@@ -169,11 +173,16 @@ map_rect(const Eigen::Matrix<T_shared_param, Eigen::Dynamic, 1>& shared_params,
169173
}
170174

171175
#ifdef STAN_MPI
172-
return internal::map_rect_mpi<call_id, F, T_shared_param, T_job_param>(
173-
shared_params, job_params, x_r, x_i, msgs);
176+
using T_plain_shared_param = plain_type_t<T_shared_param>;
177+
T_plain_shared_param shared_params_eval = shared_params;
178+
return internal::map_rect_mpi<call_id, F, T_plain_shared_param, T_job_param>(
179+
shared_params_eval, job_params, x_r, x_i, msgs);
174180
#else
175-
return internal::map_rect_concurrent<call_id, F, T_shared_param, T_job_param>(
176-
shared_params, job_params, x_r, x_i, msgs);
181+
using T_shared_param_ref = ref_type_t<T_shared_param>;
182+
T_shared_param_ref shared_params_ref = shared_params;
183+
return internal::map_rect_concurrent<call_id, F, T_shared_param_ref,
184+
T_job_param>(shared_params_ref,
185+
job_params, x_r, x_i, msgs);
177186
#endif
178187
}
179188

‎stan/math/prim/functor/map_rect_combine.hpp

+4-3
Original file line numberDiff line numberDiff line change
@@ -34,10 +34,11 @@ namespace internal {
3434
* @tparam T_shared_param type of shared parameters
3535
* @tparam T_job_param type of job specific parameters
3636
*/
37-
template <typename F, typename T_shared_param, typename T_job_param>
37+
template <typename F, typename T_shared_param, typename T_job_param,
38+
require_eigen_col_vector_t<T_shared_param>* = nullptr>
3839
class map_rect_combine {
3940
using ops_partials_t
40-
= operands_and_partials<Eigen::Matrix<T_shared_param, Eigen::Dynamic, 1>,
41+
= operands_and_partials<T_shared_param,
4142
Eigen::Matrix<T_job_param, Eigen::Dynamic, 1>>;
4243
std::vector<ops_partials_t> ops_partials_;
4344

@@ -51,7 +52,7 @@ class map_rect_combine {
5152
map_rect_combine()
5253
: ops_partials_(), num_shared_operands_(0), num_job_operands_(0) {}
5354
map_rect_combine(
54-
const Eigen::Matrix<T_shared_param, Eigen::Dynamic, 1>& shared_params,
55+
const T_shared_param& shared_params,
5556
const std::vector<Eigen::Matrix<T_job_param, Eigen::Dynamic, 1>>&
5657
job_params)
5758
: ops_partials_(),

‎stan/math/prim/functor/map_rect_concurrent.hpp

+3-2
Original file line numberDiff line numberDiff line change
@@ -10,10 +10,11 @@ namespace math {
1010
namespace internal {
1111

1212
template <int call_id, typename F, typename T_shared_param,
13-
typename T_job_param>
13+
typename T_job_param,
14+
require_eigen_col_vector_t<T_shared_param>* = nullptr>
1415
Eigen::Matrix<return_type_t<T_shared_param, T_job_param>, Eigen::Dynamic, 1>
1516
map_rect_concurrent(
16-
const Eigen::Matrix<T_shared_param, Eigen::Dynamic, 1>& shared_params,
17+
const T_shared_param& shared_params,
1718
const std::vector<Eigen::Matrix<T_job_param, Eigen::Dynamic, 1>>&
1819
job_params,
1920
const std::vector<std::vector<double>>& x_r,

‎stan/math/prim/functor/map_rect_mpi.hpp

+10-8
Original file line numberDiff line numberDiff line change
@@ -18,13 +18,14 @@ namespace internal {
1818
template <int call_id, typename F, typename T_shared_param,
1919
typename T_job_param>
2020
Eigen::Matrix<return_type_t<T_shared_param, T_job_param>, Eigen::Dynamic, 1>
21-
map_rect_mpi(
22-
const Eigen::Matrix<T_shared_param, Eigen::Dynamic, 1>& shared_params,
23-
const std::vector<Eigen::Matrix<T_job_param, Eigen::Dynamic, 1>>&
24-
job_params,
25-
const std::vector<std::vector<double>>& x_r,
26-
const std::vector<std::vector<int>>& x_i, std::ostream* msgs = nullptr) {
27-
using ReduceF = internal::map_rect_reduce<F, T_shared_param, T_job_param>;
21+
map_rect_mpi(const T_shared_param& shared_params,
22+
const std::vector<Eigen::Matrix<T_job_param, Eigen::Dynamic, 1>>&
23+
job_params,
24+
const std::vector<std::vector<double>>& x_r,
25+
const std::vector<std::vector<int>>& x_i,
26+
std::ostream* msgs = nullptr) {
27+
using ReduceF = internal::map_rect_reduce<F, scalar_type_t<T_shared_param>,
28+
T_job_param>;
2829
using CombineF = internal::map_rect_combine<F, T_shared_param, T_job_param>;
2930

3031
// whenever the cluster is already busy with some command we fall
@@ -51,7 +52,8 @@ map_rect_mpi(
5152
typedef FUNCTOR mpi_mr_##CALLID##_##SHARED##_##JOB##_; \
5253
typedef map_rect_reduce<mpi_mr_##CALLID##_##SHARED##_##JOB##_, SHARED, JOB> \
5354
mpi_mr_##CALLID##_##SHARED##_##JOB##_red_; \
54-
typedef map_rect_combine<mpi_mr_##CALLID##_##SHARED##_##JOB##_, SHARED, JOB> \
55+
typedef map_rect_combine<mpi_mr_##CALLID##_##SHARED##_##JOB##_, \
56+
Eigen::Matrix<SHARED, Eigen::Dynamic, 1>, JOB> \
5557
mpi_mr_##CALLID##_##SHARED##_##JOB##_comb_; \
5658
typedef mpi_parallel_call<CALLID, mpi_mr_##CALLID##_##SHARED##_##JOB##_red_, \
5759
mpi_mr_##CALLID##_##SHARED##_##JOB##_comb_> \

‎stan/math/prim/functor/mpi_parallel_call.hpp

+1-1
Original file line numberDiff line numberDiff line change
@@ -202,7 +202,7 @@ class mpi_parallel_call {
202202
*/
203203
template <typename T_shared_param, typename T_job_param>
204204
mpi_parallel_call(
205-
const Eigen::Matrix<T_shared_param, Eigen::Dynamic, 1>& shared_params,
205+
const T_shared_param& shared_params,
206206
const std::vector<Eigen::Matrix<T_job_param, Eigen::Dynamic, 1>>&
207207
job_params,
208208
const std::vector<std::vector<double>>& x_r,

0 commit comments

Comments
 (0)