map_rect.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_MAT_FUNCTOR_MAP_RECT_HPP
2 #define STAN_MATH_PRIM_MAT_FUNCTOR_MAP_RECT_HPP
3 
7 
8 #define STAN_REGISTER_MAP_RECT(CALLID, FUNCTOR)
9 
10 #ifdef STAN_MPI
12 #else
14 #endif
15 
16 #include <vector>
17 
18 namespace stan {
19 namespace math {
20 
21 /**
22  * Map N function evaluations to parameters and data which are in
23  * rectangular format. Each function evaluation may return a column
24  * vector of different sizes and the output is the concatenated vector
25  * from all evaluations.
26  *
27  * In addition to job specific parameters, real and int data, a shared
28  * parameter vector is repeated in all evaluations. All input
29  * parameters are stored as vectors whereas data is stored as arrays.
30  *
31  * For N jobs the output of this function is
32  *
33  * [
34  * f(shared_params, job_params[1], x_r[1], x_i[1]),
35  * f(shared_params, job_params[2], x_r[2], x_i[2]),
36  * ... ]'.
37  *
38  * The function is implemented with serial execution and with
39  * parallelism using threading or MPI (TODO). The threading version is
40  * used if the compiler flag STAN_THREADS is set during compilation
41  * while the MPI version is only available if STAN_MPI is defined. The
42  * MPI parallelism takes precedence over serial or threading execution
43  * of the function.
44  *
45  * For the threaded parallelism the N jobs are chunked into T blocks
46  * which are executed asynchronously using the async C++11
47  * facility. This ensure that at most T threads are used, but the
48  * actual number of threads is controlled by the implementation of
49  * async provided by the compiler. Note that nested calls of map_rect
50  * will lead to a multiplicative increase in the number of job chunks
51  * generated. The number of threads T is controlled at runtime via the
52  * STAN_NUM_threads environment variable, see the get_num_threads
53  * function for details.
54  *
55  * For the MPI version to work this function has these special
56  * non-standard conventions:
57  *
58  * - The call_id template parameter is considered as a label for the
59  * functor F and data combination. Since MPI communication is
60  * expensive, the real and int data is transmitted only a single
61  * time per functor F / call_id combination to the workers.
62  * - The MPI implementation requires that the functor type fully
63  * specifies the functor and hence requires a default constructible
64  * function object. This choice reduces the need for communication
65  * across MPI as the type is sufficient and the state of the
66  * functor is not necessary to transmit. Thus, the functor is
67  * specified as template argument only.
68  * - The size of the returned vector of each job must stay consistent
69  * when performing repeated calls to map_rect (which usually vary
70  * the values of the parameters).
71  * - To achieve the exact same results between the serial and the MPI
72  * evaluation scheme both variants work in exactly the same way
73  * which is to use nested AD calls and pre-computing all gradients
74  * when the function gets called. The usual evaluation scheme would
75  * build an AD graph instead of doing on-the-spot gradient
76  * evaluation. For large problems this results in speedups for the
77  * serial version even on a single core due to smaller AD graph
78  * sizes.
79  * - In MPI operation mode, no outputs from the workers are streamed
80  * back to the root.
81  * - Finally, each map_rect call must be registered with the
82  * STAN_REGISTER_MAP_RECT macro. This is required to enable
83  * communication between processes for the MPI case. The macro
84  * definition is empty if MPI is not enabled.
85  *
86  * The functor F is expected to have the usual operator() function
87  * with signature
88  *
89  * template <typename T1, typename T2>
90  * Eigen::Matrix<typename stan::return_type<T1, T2>::type, Eigen::Dynamic, 1>
91  * operator()(const Eigen::Matrix<T1, Eigen::Dynamic, 1>& eta,
92  * const Eigen::Matrix<T2, Eigen::Dynamic, 1>& theta,
93  * const std::vector<double>& x_r, const std::vector<int>& x_i,
94  * std::ostream* msgs = 0) const { ... }
95  *
96  * WARNING: For the MPI case, the data arguments are NOT checked if
97  * they are unchanged between repeated evaluations for a given
98  * call_id/functor F pair. This is silently assumed to be immutable
99  * between evaluations.
100  *
101  * @tparam T_shared_param Type of shared parameters.
102  * @tparam T_job_param Type of job specific parameters.
103  * @param shared_params shared parameter vector passed as first
104  * argument to functor for all jobs
105  * @param job_params Array of job specific parameter vectors. All job
106  * specific parameters must have matching sizes.
107  * @param x_r Array of real arrays for each job. The first dimension
108  * must match the number of jobs (which is determined by the first
109  * dimension of job_params) and each entry must have the same size.
110  * @param x_i Array of int data with the same conventions as x_r.
111  * @param msgs Output stream for messages.
112  * @tparam call_id Label for functor/data combination. See above for
113  * details.
114  * @tparam F Functor which is applied to all job specific parameters
115  * with conventions described.
116  * @return concatenated results from all jobs
117  */
118 
119 template <int call_id, typename F, typename T_shared_param,
120  typename T_job_param>
122  Eigen::Dynamic, 1>
123 map_rect(const Eigen::Matrix<T_shared_param, Eigen::Dynamic, 1>& shared_params,
124  const std::vector<Eigen::Matrix<T_job_param, Eigen::Dynamic, 1>>&
125  job_params,
126  const std::vector<std::vector<double>>& x_r,
127  const std::vector<std::vector<int>>& x_i,
128  std::ostream* msgs = nullptr) {
129  static const char* function = "map_rect";
130  typedef Eigen::Matrix<
132  Eigen::Dynamic, 1>
133  return_t;
134 
135  check_matching_sizes(function, "job parameters", job_params, "real data",
136  x_r);
137  check_matching_sizes(function, "job parameters", job_params, "int data", x_i);
138 
139  // check size consistency of inputs per job
140  const std::vector<int> job_params_dims = dims(job_params);
141  const int size_job_params = job_params_dims[1];
142  const int size_x_r = dims(x_r)[1];
143  const int size_x_i = dims(x_i)[1];
144  for (int i = 1; i < job_params_dims[0]; i++) {
145  check_size_match(function,
146  "Size of one of the vectors of "
147  "the job specific parameters",
148  job_params[i].size(),
149  "size of another vector of the "
150  "job specifc parameters",
151  size_job_params);
152  check_size_match(function,
153  "Size of one of the arrays of "
154  "the job specific real data",
155  x_r[i].size(),
156  "size of another array of the "
157  "job specifc real data",
158  size_x_r);
159  check_size_match(function,
160  "Size of one of the arrays of "
161  "the job specific int data",
162  x_i[i].size(),
163  "size of another array of the "
164  "job specifc int data",
165  size_x_i);
166  }
167 
168  if (job_params_dims[0] == 0)
169  return return_t();
170 
171 #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);
174 #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);
177 #endif
178 }
179 
180 } // namespace math
181 } // namespace stan
182 
183 #endif
boost::math::tools::promote_args< typename scalar_type< T1 >::type, typename scalar_type< T2 >::type, typename scalar_type< T3 >::type, typename scalar_type< T4 >::type, typename scalar_type< T5 >::type, typename scalar_type< T6 >::type >::type type
Definition: return_type.hpp:20
#define F(x, y, z)
void check_size_match(const char *function, const char *name_i, T_size1 i, const char *name_j, T_size2 j)
void dims(const T &x, std::vector< int > &result)
Definition: dims.hpp:11
int size(const std::vector< T > &x)
Definition: size.hpp:17
void check_matching_sizes(const char *function, const char *name1, const T_y1 &y1, const char *name2, const T_y2 &y2)
Eigen::Matrix< typename stan::return_type< T_shared_param, T_job_param >::type, Eigen::Dynamic, 1 > map_rect(const Eigen::Matrix< T_shared_param, Eigen::Dynamic, 1 > &shared_params, const std::vector< Eigen::Matrix< T_job_param, Eigen::Dynamic, 1 >> &job_params, const std::vector< std::vector< double >> &x_r, const std::vector< std::vector< int >> &x_i, std::ostream *msgs=nullptr)
Definition: map_rect.hpp:123