1 #ifndef STAN_SERVICES_UTIL_INITIALIZE_HPP
2 #define STAN_SERVICES_UTIL_INITIALIZE_HPP
3
4 #include <stan/callbacks/logger.hpp>
5 #include <stan/callbacks/writer.hpp>
6 #include <stan/io/var_context.hpp>
7 #include <stan/io/random_var_context.hpp>
8 #include <stan/io/chained_var_context.hpp>
9 #include <stan/model/log_prob_grad.hpp>
10 #include <stan/math/prim.hpp>
11 #include <chrono>
12 #include <sstream>
13 #include <string>
14 #include <vector>
15
16 namespace stan {
17 namespace services {
18 namespace util {
19
20 /**
21 * Returns a valid initial value of the parameters of the model
22 * on the unconstrained scale.
23 *
24 * For identical inputs (model, init, rng, init_radius), this
25 * function will produce the same initialization.
26 *
27 * Initialization first tries to use the provided
28 * <code>stan::io::var_context</code>, then it will generate
29 * random uniform values from -init_radius to +init_radius for missing
30 * parameters.
31 *
32 * When the <code>var_context</code> provides all variables or
33 * the init_radius is 0, this function will only evaluate the
34 * log probability of the model with the unconstrained
35 * parameters once to see if it's valid.
36 *
37 * When at least some of the initialization is random, it will
38 * randomly initialize until it finds a set of unconstrained
39 * parameters that are valid or it hits <code>MAX_INIT_TRIES =
40 * 100</code> (hard-coded).
41 *
42 * Valid initialization is defined as a finite, non-NaN value for the
43 * evaluation of the log probability density function and all its
44 * gradients.
45 *
46 * @tparam Jacobian indicates whether to include the Jacobian term when
47 * evaluating the log density function
48 * @tparam Model the type of the model class
49 * @tparam RNG the type of the random number generator
50 *
51 * @param[in] model the model
52 * @param[in] init a var_context with initial values
53 * @param[in,out] rng random number generator
54 * @param[in] init_radius the radius for generating random values.
55 * A value of 0 indicates that the unconstrained parameters (not
56 * provided by init) should be initialized with 0.
57 * @param[in] print_timing indicates whether a timing message should
58 * be printed to the logger
59 * @param[in,out] logger logger for messages
60 * @param[in,out] init_writer init writer (on the unconstrained scale)
61 * @throws exception passed through from the model if the model has a
62 * fatal error (not a std::domain_error)
63 * @throws std::domain_error if the model can not be initialized and
64 * the model does not have a fatal error (only allows for
65 * std::domain_error)
66 * @return valid unconstrained parameters for the model
67 */
68 template <bool Jacobian = true, typename Model, typename InitContext,
69 typename RNG>
initialize(Model & model,const InitContext & init,RNG & rng,double init_radius,bool print_timing,stan::callbacks::logger & logger,stan::callbacks::writer & init_writer)70 std::vector<double> initialize(Model& model, const InitContext& init, RNG& rng,
71 double init_radius, bool print_timing,
72 stan::callbacks::logger& logger,
73 stan::callbacks::writer& init_writer) {
74 std::vector<double> unconstrained;
75 std::vector<int> disc_vector;
76
77 bool is_fully_initialized = true;
78 bool any_initialized = false;
79 std::vector<std::string> param_names;
80 model.get_param_names(param_names);
81 for (size_t n = 0; n < param_names.size(); n++) {
82 is_fully_initialized &= init.contains_r(param_names[n]);
83 any_initialized |= init.contains_r(param_names[n]);
84 }
85
86 bool is_initialized_with_zero = init_radius == 0.0;
87
88 int MAX_INIT_TRIES
89 = is_fully_initialized || is_initialized_with_zero ? 1 : 100;
90 int num_init_tries = 0;
91 for (; num_init_tries < MAX_INIT_TRIES; num_init_tries++) {
92 std::stringstream msg;
93 try {
94 stan::io::random_var_context random_context(model, rng, init_radius,
95 is_initialized_with_zero);
96
97 if (!any_initialized) {
98 unconstrained = random_context.get_unconstrained();
99 } else {
100 stan::io::chained_var_context context(init, random_context);
101
102 model.transform_inits(context, disc_vector, unconstrained, &msg);
103 }
104 } catch (std::domain_error& e) {
105 if (msg.str().length() > 0)
106 logger.info(msg);
107 logger.info("Rejecting initial value:");
108 logger.info(
109 " Error evaluating the log probability"
110 " at the initial value.");
111 logger.info(e.what());
112 continue;
113 } catch (std::exception& e) {
114 if (msg.str().length() > 0)
115 logger.info(msg);
116 logger.info(
117 "Unrecoverable error evaluating the log probability"
118 " at the initial value.");
119 logger.info(e.what());
120 throw;
121 }
122
123 msg.str("");
124 double log_prob(0);
125 try {
126 // we evaluate the log_prob function with propto=false
127 // because we're evaluating with `double` as the type of
128 // the parameters.
129 log_prob = model.template log_prob<false, Jacobian>(unconstrained,
130 disc_vector, &msg);
131 if (msg.str().length() > 0)
132 logger.info(msg);
133 } catch (std::domain_error& e) {
134 if (msg.str().length() > 0)
135 logger.info(msg);
136 logger.info("Rejecting initial value:");
137 logger.info(
138 " Error evaluating the log probability"
139 " at the initial value.");
140 logger.info(e.what());
141 continue;
142 } catch (std::exception& e) {
143 if (msg.str().length() > 0)
144 logger.info(msg);
145 logger.info(
146 "Unrecoverable error evaluating the log probability"
147 " at the initial value.");
148 logger.info(e.what());
149 throw;
150 }
151 if (!std::isfinite(log_prob)) {
152 logger.info("Rejecting initial value:");
153 logger.info(
154 " Log probability evaluates to log(0),"
155 " i.e. negative infinity.");
156 logger.info(
157 " Stan can't start sampling from this"
158 " initial value.");
159 continue;
160 }
161 std::stringstream log_prob_msg;
162 std::vector<double> gradient;
163 auto start = std::chrono::steady_clock::now();
164 try {
165 // we evaluate this with propto=true since we're
166 // evaluating with autodiff variables
167 log_prob = stan::model::log_prob_grad<true, Jacobian>(
168 model, unconstrained, disc_vector, gradient, &log_prob_msg);
169 } catch (const std::exception& e) {
170 if (log_prob_msg.str().length() > 0)
171 logger.info(log_prob_msg);
172 logger.info(e.what());
173 throw;
174 }
175 auto end = std::chrono::steady_clock::now();
176 double deltaT
177 = std::chrono::duration_cast<std::chrono::microseconds>(end - start)
178 .count()
179 / 1000000.0;
180 if (log_prob_msg.str().length() > 0)
181 logger.info(log_prob_msg);
182
183 bool gradient_ok = std::isfinite(stan::math::sum(gradient));
184
185 if (!gradient_ok) {
186 logger.info("Rejecting initial value:");
187 logger.info(
188 " Gradient evaluated at the initial value"
189 " is not finite.");
190 logger.info(
191 " Stan can't start sampling from this"
192 " initial value.");
193 }
194 if (gradient_ok && print_timing) {
195 logger.info("");
196 std::stringstream msg1;
197 msg1 << "Gradient evaluation took " << deltaT << " seconds";
198 logger.info(msg1);
199
200 std::stringstream msg2;
201 msg2 << "1000 transitions using 10 leapfrog steps"
202 << " per transition would take"
203 << " " << 1e4 * deltaT << " seconds.";
204 logger.info(msg2);
205
206 logger.info("Adjust your expectations accordingly!");
207 logger.info("");
208 logger.info("");
209 }
210 if (gradient_ok) {
211 init_writer(unconstrained);
212 return unconstrained;
213 }
214 }
215
216 if (!is_initialized_with_zero) {
217 logger.info("");
218 std::stringstream msg;
219 msg << "Initialization between (-" << init_radius << ", " << init_radius
220 << ") failed after"
221 << " " << MAX_INIT_TRIES << " attempts. ";
222 logger.info(msg);
223 logger.info(
224 " Try specifying initial values,"
225 " reducing ranges of constrained values,"
226 " or reparameterizing the model.");
227 }
228 throw std::domain_error("Initialization failed.");
229 }
230
231 } // namespace util
232 } // namespace services
233 } // namespace stan
234 #endif
235