1 #ifndef STAN_SERVICES_INIT_INITIALIZE_STATE_HPP
2 #define STAN_SERVICES_INIT_INITIALIZE_STATE_HPP
4 #include <boost/lexical_cast.hpp>
5 #include <boost/random/additive_combine.hpp>
6 #include <boost/random/uniform_real_distribution.hpp>
7 #include <boost/random/variate_generator.hpp>
15 #include <stan/math/prim/scal/fun/is_inf.hpp>
16 #include <stan/math/prim/scal/fun/is_nan.hpp>
17 #include <stan/math/prim/mat/fun/Eigen.hpp>
46 void rm_indices_from_name(std::string& name) {
47 size_t x = name.find_first_of(
".");
48 if (std::string::npos == x)
return;
57 void rm_indices_from_name(std::vector<std::string>& names) {
58 for (
size_t i = 0; i < names.size(); i++)
59 rm_indices_from_name(names[i]);
60 std::vector<std::string>::iterator it;
61 it = std::unique(names.begin(), names.end());
62 names.resize(std::distance(names.begin(), it));
78 template <
class Model>
79 bool are_all_pars_initialized(
const Model& model,
81 std::vector<std::string> names;
82 model.constrained_param_names(names,
false,
false);
83 rm_indices_from_name(names);
84 for (
size_t i = 0; i < names.size(); i++)
85 if (!context.
contains_r(names[i]))
return false;
89 template <
class Model>
90 bool validate_unconstrained_initialization(Eigen::VectorXd& cont_params,
92 for (
int n = 0; n < cont_params.size(); n++) {
93 if (stan::math::is_inf(cont_params[n])
94 || stan::math::is_nan(cont_params[n])) {
95 std::vector<std::string> param_names;
96 model.unconstrained_param_names(param_names);
98 std::stringstream msg;
99 msg << param_names[n] <<
" initialized to invalid value ("
100 << cont_params[n] <<
")";
102 throw std::invalid_argument(msg.str());
118 template <
class Model>
124 validate_unconstrained_initialization(cont_params, model);
125 }
catch (
const std::exception& e) {
129 double init_log_prob;
130 Eigen::VectorXd init_grad = Eigen::VectorXd::Zero(model.num_params_r());
134 }
catch (
const std::exception& e) {
137 writer(
"Rejecting initial value:");
138 writer(
" Error evaluating the log probability "
139 "at the initial value.");
142 if (!boost::math::isfinite(init_log_prob)) {
143 writer(
"Rejecting initial value:");
144 writer(
" Log probability evaluates to log(0), "
145 "i.e. negative infinity.");
146 writer(
" Stan can't start sampling from this initial value.");
149 for (
int i = 0; i < init_grad.size(); ++i) {
150 if (!boost::math::isfinite(init_grad(i))) {
151 writer(
"Rejecting initial value:");
152 writer(
" Gradient evaluated at the initial value "
154 writer(
" Stan can't start sampling from this initial value.");
171 template <
class Model>
175 cont_params.setZero();
193 template <
class Model,
class RNG>
195 Eigen::VectorXd& cont_params,
199 int num_init_tries = -1;
201 boost::random::uniform_real_distribution<double>
202 init_range_distribution(-R, R);
204 boost::variate_generator
205 <RNG&, boost::random::uniform_real_distribution<double> >
206 init_rng(base_rng, init_range_distribution);
209 static int MAX_INIT_TRIES = 100;
211 for (num_init_tries = 1; num_init_tries <= MAX_INIT_TRIES;
213 for (
int i = 0; i < cont_params.size(); ++i)
214 cont_params(i) = init_rng();
219 if (num_init_tries > MAX_INIT_TRIES) {
220 std::stringstream R_ss, MAX_INIT_TRIES_ss;
222 MAX_INIT_TRIES_ss << MAX_INIT_TRIES;
226 writer(
"Initialization between (-" + R_ss.str() +
", " + R_ss.str()
228 + MAX_INIT_TRIES_ss.str() +
" attempts. ");
229 writer(
" Try specifying initial values,"
230 " reducing ranges of constrained values,"
231 " or reparameterizing the model.");
256 template <
class ContextFactory,
class Model,
class RNG>
259 Eigen::VectorXd& cont_params,
263 ContextFactory& context_factory) {
265 boost::random::uniform_real_distribution<double>
266 init_range_distribution(-R, R);
267 boost::variate_generator
268 <RNG&, boost::random::uniform_real_distribution<double> >
269 init_rng(base_rng, init_range_distribution);
271 static int MAX_INIT_TRIES = 100;
273 int num_init_tries = -1;
274 std::vector<std::string> cont_names;
275 model.constrained_param_names(cont_names,
false,
false);
276 rm_indices_from_name(cont_names);
277 for (num_init_tries = 1; num_init_tries <= MAX_INIT_TRIES;
279 std::vector<double> cont_vecs(cont_params.size());
280 for (
int i = 0; i < cont_params.size(); ++i) {
281 cont_vecs[i] = init_rng();
282 cont_params[i] = cont_vecs[i];
284 typename ContextFactory::var_context_t context
285 = context_factory(source);
286 std::vector<double> cont_vecs_constrained;
287 std::vector<int> int_vecs;
288 model.write_array(base_rng, cont_vecs, int_vecs,
289 cont_vecs_constrained,
false,
false, 0);
290 std::vector<std::vector<size_t> > dims;
291 model.get_dims(dims);
293 cont_vecs_constrained,
296 model.transform_inits(cvc, cont_params, 0);
301 if (num_init_tries > MAX_INIT_TRIES) {
302 std::stringstream R_ss, MAX_INIT_TRIES_ss;
304 MAX_INIT_TRIES_ss << MAX_INIT_TRIES;
308 writer(
"Initialization between (-" + R_ss.str() +
", " + R_ss.str()
310 + MAX_INIT_TRIES_ss.str() +
" attempts. ");
311 writer(
" Try specifying initial values,"
312 " reducing ranges of constrained values,"
313 " or reparameterizing the model.");
316 }
catch(
const std::exception& e) {
317 writer(
"Initialization partially from source failed.");
345 template <
class ContextFactory,
class Model,
class RNG>
347 Eigen::VectorXd& cont_params,
352 ContextFactory& context_factory,
353 bool enable_random_init =
false,
356 typename ContextFactory::var_context_t context
357 = context_factory(source);
359 if (enable_random_init && !are_all_pars_initialized(model, context)) {
368 model.transform_inits(context, cont_params, 0);
369 }
catch(
const std::exception& e) {
370 writer(
"Initialization from source failed.");
388 val = boost::lexical_cast<
double>(s);
389 }
catch (
const boost::bad_lexical_cast& e) {
390 val = std::numeric_limits<double>::quiet_NaN();
416 template <
class ContextFactory,
class Model,
class RNG>
418 Eigen::VectorXd& cont_params,
422 ContextFactory& context_factory,
423 bool enable_random_init =
false,
void write_error_msg(interface_callbacks::writer::base_writer &writer, const std::exception &e)
Writes a Metropolis rejection message.
bool initialize_state_values(Eigen::VectorXd &cont_params, Model &model, interface_callbacks::writer::base_writer &writer)
bool initialize_state_source(const std::string source, Eigen::VectorXd &cont_params, Model &model, RNG &base_rng, interface_callbacks::writer::base_writer &writer, ContextFactory &context_factory, bool enable_random_init=false, double R=2)
Creates the initial state using the source parameter.
Probability, optimization and sampling library.
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)
virtual bool contains_r(const std::string &name) const =0
Return true if the specified variable name is defined.
bool initialize_state(const std::string &init, Eigen::VectorXd &cont_params, Model &model, RNG &base_rng, interface_callbacks::writer::base_writer &writer, ContextFactory &context_factory, bool enable_random_init=false, double init_r=2)
Creates the initial state.
bool get_double_from_string(const std::string &s, double &val)
Converts string to double.
A var_reader reads array variables of integer and floating point type by name and dimension...
bool initialize_state_random(const double R, Eigen::VectorXd &cont_params, Model &model, RNG &base_rng, interface_callbacks::writer::base_writer &writer)
Initializes state to random uniform values within range.
base_writer is an abstract base class defining the interface for Stan writer callbacks.
bool initialize_state_zero(Eigen::VectorXd &cont_params, Model &model, interface_callbacks::writer::base_writer &writer)
Sets initial state to zero.
bool initialize_state_source_and_random(const std::string &source, double R, Eigen::VectorXd &cont_params, Model &model, RNG &base_rng, interface_callbacks::writer::base_writer &writer, ContextFactory &context_factory)
Creates the initial state.
A chained_var_context object represents two objects of var_context as one.
An array_var_context object represents a named arrays with dimensions constructed from an array...