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