From ae33787f3e4632f58848b51f1a0a9cdc2beb5acf Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Fri, 6 Jul 2018 16:13:15 -0400 Subject: [PATCH 1/9] Updating .gitignore --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index c0df14bafc0..3dde639abb1 100644 --- a/.gitignore +++ b/.gitignore @@ -56,3 +56,5 @@ stan.kdev4 # local make include /make/local + +*.gch \ No newline at end of file From 59601381d2e524793faa76aeec13542167372914 Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Fri, 6 Jul 2018 20:34:01 -0400 Subject: [PATCH 2/9] updating whitespace --- src/stan/services/util/initialize.hpp | 370 +++++++++++++------------- 1 file changed, 185 insertions(+), 185 deletions(-) diff --git a/src/stan/services/util/initialize.hpp b/src/stan/services/util/initialize.hpp index 74613ee6b6b..383dd28533d 100644 --- a/src/stan/services/util/initialize.hpp +++ b/src/stan/services/util/initialize.hpp @@ -13,202 +13,202 @@ #include namespace stan { - namespace services { - namespace util { +namespace services { +namespace util { - /** - * Returns a valid initial value of the parameters of the model - * on the unconstrained scale. - * - * For identical inputs (model, init, rng, init_radius), this - * function will produce the same initialization. - * - * Initialization first tries to use the provided - * stan::io::var_context, then it will generate - * random uniform values from -init_radius to +init_radius for missing - * parameters. - * - * When the var_context provides all variables or - * the init_radius is 0, this function will only evaluate the - * log probability of the model with the unconstrained - * parameters once to see if it's valid. - * - * When at least some of the initialization is random, it will - * randomly initialize until it finds a set of unconstrained - * parameters that are valid or it hits MAX_INIT_TRIES = - * 100 (hard-coded). - * - * Valid initialization is defined as a finite, non-NaN value - * for the evaluation of the log probability and all its - * gradients. - * - * @param[in] model the model - * @param[in] init a var_context with initial values - * @param[in,out] rng random number generator - * @param[in] init_radius the radius for generating random values. - * A value of 0 indicates that the unconstrained parameters (not - * provided by init) should be initialized with 0. - * @param[in] print_timing indicates whether a timing message should - * be printed to the logger - * @param[in,out] logger logger for messages - * @param[in,out] init_writer init writer (on the unconstrained scale) - * @throws exception passed through from the model if the model has a - * fatal error (not a std::domain_error) - * @throws std::domain_error if the model can not be initialized and - * the model does not have a fatal error (only allows for - * std::domain_error) - * @return valid unconstrained parameters for the model - */ - template - std::vector 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) { - std::vector unconstrained; - std::vector disc_vector; +/** + * Returns a valid initial value of the parameters of the model + * on the unconstrained scale. + * + * For identical inputs (model, init, rng, init_radius), this + * function will produce the same initialization. + * + * Initialization first tries to use the provided + * stan::io::var_context, then it will generate + * random uniform values from -init_radius to +init_radius for missing + * parameters. + * + * When the var_context provides all variables or + * the init_radius is 0, this function will only evaluate the + * log probability of the model with the unconstrained + * parameters once to see if it's valid. + * + * When at least some of the initialization is random, it will + * randomly initialize until it finds a set of unconstrained + * parameters that are valid or it hits MAX_INIT_TRIES = + * 100 (hard-coded). + * + * Valid initialization is defined as a finite, non-NaN value + * for the evaluation of the log probability and all its + * gradients. + * + * @param[in] model the model + * @param[in] init a var_context with initial values + * @param[in,out] rng random number generator + * @param[in] init_radius the radius for generating random values. + * A value of 0 indicates that the unconstrained parameters (not + * provided by init) should be initialized with 0. + * @param[in] print_timing indicates whether a timing message should + * be printed to the logger + * @param[in,out] logger logger for messages + * @param[in,out] init_writer init writer (on the unconstrained scale) + * @throws exception passed through from the model if the model has a + * fatal error (not a std::domain_error) + * @throws std::domain_error if the model can not be initialized and + * the model does not have a fatal error (only allows for + * std::domain_error) + * @return valid unconstrained parameters for the model + */ +template +std::vector 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) { + std::vector unconstrained; + std::vector disc_vector; - bool is_fully_initialized = true; - bool any_initialized = false; - std::vector param_names; - model.get_param_names(param_names); - for (size_t n = 0; n < param_names.size(); n++) { - is_fully_initialized &= init.contains_r(param_names[n]); - any_initialized |= init.contains_r(param_names[n]); - } - - bool init_zero = init_radius <= std::numeric_limits::min(); - - int MAX_INIT_TRIES = is_fully_initialized || init_zero ? 1 : 100; - int num_init_tries = 0; - for (; num_init_tries < MAX_INIT_TRIES; num_init_tries++) { - stan::io::random_var_context - random_context(model, rng, init_radius, init_zero); + bool is_fully_initialized = true; + bool any_initialized = false; + std::vector param_names; + model.get_param_names(param_names); + for (size_t n = 0; n < param_names.size(); n++) { + is_fully_initialized &= init.contains_r(param_names[n]); + any_initialized |= init.contains_r(param_names[n]); + } - if (!any_initialized) { - unconstrained = random_context.get_unconstrained(); - } else { - stan::io::chained_var_context context(init, random_context); + bool init_zero = init_radius <= std::numeric_limits::min(); - std::stringstream msg; - try { - model.transform_inits(context, - disc_vector, - unconstrained, - &msg); - } catch (const std::exception& e) { - if (msg.str().length() > 0) - logger.info(msg); - logger.info(e.what()); - throw; - } - if (msg.str().length() > 0) - logger.info(msg); - } - double log_prob(0); - std::stringstream msg; - try { - log_prob = model.template log_prob - (unconstrained, disc_vector, &msg); - if (msg.str().length() > 0) - logger.info(msg); - } catch (std::domain_error& e) { - if (msg.str().length() > 0) - logger.info(msg); - logger.info("Rejecting initial value:"); - logger.info(" Error evaluating the log probability" - " at the initial value."); - logger.info(e.what()); - continue; - } catch (std::exception& e) { - if (msg.str().length() > 0) - logger.info(msg); - logger.info("Unrecoverable error evaluating the log probability" - " at the initial value."); - logger.info(e.what()); - throw; - } - if (!boost::math::isfinite(log_prob)) { - logger.info("Rejecting initial value:"); - logger.info(" Log probability evaluates to log(0)," - " i.e. negative infinity."); - logger.info(" Stan can't start sampling from this" - " initial value."); - continue; - } - std::stringstream log_prob_msg; - std::vector gradient; - bool gradient_ok = true; - clock_t start_check = clock(); - try { - log_prob = stan::model::log_prob_grad - (model, unconstrained, disc_vector, - gradient, &log_prob_msg); - } catch (const std::exception& e) { - if (log_prob_msg.str().length() > 0) - logger.info(log_prob_msg); - logger.info(e.what()); - throw; - } - clock_t end_check = clock(); - double deltaT = static_cast(end_check - start_check) - / CLOCKS_PER_SEC; - if (log_prob_msg.str().length() > 0) - logger.info(log_prob_msg); + int MAX_INIT_TRIES = is_fully_initialized || init_zero ? 1 : 100; + int num_init_tries = 0; + for (; num_init_tries < MAX_INIT_TRIES; num_init_tries++) { + stan::io::random_var_context + random_context(model, rng, init_radius, init_zero); - for (size_t i = 0; i < gradient.size(); ++i) { - if (gradient_ok && !boost::math::isfinite(gradient[i])) { - logger.info("Rejecting initial value:"); - logger.info(" Gradient evaluated at the initial value" - " is not finite."); - logger.info(" Stan can't start sampling from this" - " initial value."); - gradient_ok = false; - } - } - if (gradient_ok && print_timing) { - logger.info(""); - std::stringstream msg1; - msg1 << "Gradient evaluation took " << deltaT << " seconds"; - logger.info(msg1); + if (!any_initialized) { + unconstrained = random_context.get_unconstrained(); + } else { + stan::io::chained_var_context context(init, random_context); - std::stringstream msg2; - msg2 << "1000 transitions using 10 leapfrog steps" - << " per transition would take" - << " " << 1e4 * deltaT << " seconds."; - logger.info(msg2); + std::stringstream msg; + try { + model.transform_inits(context, + disc_vector, + unconstrained, + &msg); + } catch (const std::exception& e) { + if (msg.str().length() > 0) + logger.info(msg); + logger.info(e.what()); + throw; + } + if (msg.str().length() > 0) + logger.info(msg); + } + double log_prob(0); + std::stringstream msg; + try { + log_prob = model.template log_prob + (unconstrained, disc_vector, &msg); + if (msg.str().length() > 0) + logger.info(msg); + } catch (std::domain_error& e) { + if (msg.str().length() > 0) + logger.info(msg); + logger.info("Rejecting initial value:"); + logger.info(" Error evaluating the log probability" + " at the initial value."); + logger.info(e.what()); + continue; + } catch (std::exception& e) { + if (msg.str().length() > 0) + logger.info(msg); + logger.info("Unrecoverable error evaluating the log probability" + " at the initial value."); + logger.info(e.what()); + throw; + } + if (!boost::math::isfinite(log_prob)) { + logger.info("Rejecting initial value:"); + logger.info(" Log probability evaluates to log(0)," + " i.e. negative infinity."); + logger.info(" Stan can't start sampling from this" + " initial value."); + continue; + } + std::stringstream log_prob_msg; + std::vector gradient; + bool gradient_ok = true; + clock_t start_check = clock(); + try { + log_prob = stan::model::log_prob_grad + (model, unconstrained, disc_vector, + gradient, &log_prob_msg); + } catch (const std::exception& e) { + if (log_prob_msg.str().length() > 0) + logger.info(log_prob_msg); + logger.info(e.what()); + throw; + } + clock_t end_check = clock(); + double deltaT = static_cast(end_check - start_check) + / CLOCKS_PER_SEC; + if (log_prob_msg.str().length() > 0) + logger.info(log_prob_msg); - logger.info("Adjust your expectations accordingly!"); - logger.info(""); - logger.info(""); - } - if (gradient_ok) - break; - } + for (size_t i = 0; i < gradient.size(); ++i) { + if (gradient_ok && !boost::math::isfinite(gradient[i])) { + logger.info("Rejecting initial value:"); + logger.info(" Gradient evaluated at the initial value" + " is not finite."); + logger.info(" Stan can't start sampling from this" + " initial value."); + gradient_ok = false; + } + } + if (gradient_ok && print_timing) { + logger.info(""); + std::stringstream msg1; + msg1 << "Gradient evaluation took " << deltaT << " seconds"; + logger.info(msg1); - if (num_init_tries == MAX_INIT_TRIES) { - if (!is_fully_initialized && !init_zero) { - logger.info(""); - std::stringstream msg; - msg << "Initialization between (-" << init_radius - << ", " << init_radius << ") failed after" - << " " << MAX_INIT_TRIES << " attempts. "; - logger.info(msg); - logger.info(" Try specifying initial values," - " reducing ranges of constrained values," - " or reparameterizing the model."); - } - throw std::domain_error("Initialization failed."); - } + std::stringstream msg2; + msg2 << "1000 transitions using 10 leapfrog steps" + << " per transition would take" + << " " << 1e4 * deltaT << " seconds."; + logger.info(msg2); - init_writer(unconstrained); - return unconstrained; - } + logger.info("Adjust your expectations accordingly!"); + logger.info(""); + logger.info(""); + } + if (gradient_ok) + break; + } + if (num_init_tries == MAX_INIT_TRIES) { + if (!is_fully_initialized && !init_zero) { + logger.info(""); + std::stringstream msg; + msg << "Initialization between (-" << init_radius + << ", " << init_radius << ") failed after" + << " " << MAX_INIT_TRIES << " attempts. "; + logger.info(msg); + logger.info(" Try specifying initial values," + " reducing ranges of constrained values," + " or reparameterizing the model."); } + throw std::domain_error("Initialization failed."); } + + init_writer(unconstrained); + return unconstrained; +} + +} +} } #endif From 344f82c3a909b90d9b2e63f61974aa18170324b3 Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Fri, 6 Jul 2018 16:14:08 -0400 Subject: [PATCH 3/9] adding failing test for issue #2562 --- .../unit/services/util/initialize_test.cpp | 427 ++++++++++++------ 1 file changed, 284 insertions(+), 143 deletions(-) diff --git a/src/test/unit/services/util/initialize_test.cpp b/src/test/unit/services/util/initialize_test.cpp index 63c8fe0ad8d..736e604f3a7 100644 --- a/src/test/unit/services/util/initialize_test.cpp +++ b/src/test/unit/services/util/initialize_test.cpp @@ -11,11 +11,11 @@ #include class ServicesUtilInitialize : public testing::Test { -public: + public: ServicesUtilInitialize() - : model(empty_context, 12345, &model_ss), - message(message_ss), - rng(stan::services::util::create_rng(0, 1)) {} + : model(empty_context, 12345, &model_ss), + message(message_ss), + rng(stan::services::util::create_rng(0, 1)) {} stan_model model; stan::io::empty_var_context empty_context; @@ -36,7 +36,7 @@ TEST_F(ServicesUtilInitialize, radius_zero__print_false) { init_radius, print_timing, logger, init); ASSERT_EQ(model.num_params_r(), params.size()) - << "2 parameters"; + << "2 parameters"; EXPECT_FLOAT_EQ(0, params[0]); EXPECT_FLOAT_EQ(0, params[1]); @@ -56,7 +56,7 @@ TEST_F(ServicesUtilInitialize, radius_zero__print_true) { init_radius, print_timing, logger, init); ASSERT_EQ(model.num_params_r(), params.size()) - << "2 parameters"; + << "2 parameters"; EXPECT_FLOAT_EQ(0, params[0]); EXPECT_FLOAT_EQ(0, params[1]); @@ -78,7 +78,7 @@ TEST_F(ServicesUtilInitialize, radius_two__print_false) { init_radius, print_timing, logger, init); ASSERT_EQ(model.num_params_r(), params.size()) - << "2 parameters"; + << "2 parameters"; EXPECT_GT(params[0], -init_radius); EXPECT_LT(params[0], init_radius); EXPECT_GT(params[1], -init_radius); @@ -100,7 +100,7 @@ TEST_F(ServicesUtilInitialize, radius_two__print_true) { init_radius, print_timing, logger, init); ASSERT_EQ(model.num_params_r(), params.size()) - << "2 parameters"; + << "2 parameters"; EXPECT_GT(params[0], -init_radius); EXPECT_LT(params[0], init_radius); EXPECT_GT(params[1], -init_radius); @@ -136,7 +136,7 @@ TEST_F(ServicesUtilInitialize, full_init__print_false) { init_radius, print_timing, logger, init); ASSERT_EQ(model.num_params_r(), params.size()) - << "2 parameters"; + << "2 parameters"; EXPECT_FLOAT_EQ(1.5, params[0]); EXPECT_FLOAT_EQ(-0.5, params[1]); @@ -168,7 +168,7 @@ TEST_F(ServicesUtilInitialize, full_init__print_true) { init_radius, print_timing, logger, init); ASSERT_EQ(model.num_params_r(), params.size()) - << "2 parameters"; + << "2 parameters"; EXPECT_FLOAT_EQ(1.5, params[0]); EXPECT_FLOAT_EQ(-0.5, params[1]); @@ -182,88 +182,89 @@ TEST_F(ServicesUtilInitialize, full_init__print_true) { } namespace test { - // Mock Throwing Model throws exception - class mock_throwing_model: public stan::model::prob_grad { - public: +// Mock Throwing Model throws exception +class mock_throwing_model: public stan::model::prob_grad { + public: - mock_throwing_model(): + mock_throwing_model(): stan::model::prob_grad(1), templated_log_prob_calls(0), transform_inits_calls(0), write_array_calls(0), log_prob_return_value(0.0) { } - void reset() { - templated_log_prob_calls = 0; - transform_inits_calls = 0; - write_array_calls = 0; - log_prob_return_value = 0.0; + void reset() { + templated_log_prob_calls = 0; + transform_inits_calls = 0; + write_array_calls = 0; + log_prob_return_value = 0.0; + } + + template + T__ log_prob(std::vector& params_r__, + std::vector& params_i__, + std::ostream* pstream__ = 0) const { + ++templated_log_prob_calls; + throw std::domain_error("throwing within log_prob"); + return log_prob_return_value; + } + + void transform_inits(const stan::io::var_context& context__, + std::vector& params_i__, + std::vector& params_r__, + std::ostream* pstream__) const { + ++transform_inits_calls; + for (size_t n = 0; n < params_r__.size(); ++n) { + params_r__[n] = n; } + } - template - T__ log_prob(std::vector& params_r__, - std::vector& params_i__, - std::ostream* pstream__ = 0) const { - ++templated_log_prob_calls; - throw std::domain_error("throwing within log_prob"); - return log_prob_return_value; - } + void get_dims(std::vector >& dimss__) const { + dimss__.resize(0); + std::vector scalar_dim; + dimss__.push_back(scalar_dim); + } - void transform_inits(const stan::io::var_context& context__, - std::vector& params_i__, - std::vector& params_r__, - std::ostream* pstream__) const { - ++transform_inits_calls; - for (size_t n = 0; n < params_r__.size(); ++n) { - params_r__[n] = n; - } - } + void constrained_param_names(std::vector& param_names__, + bool include_tparams__ = true, + bool include_gqs__ = true) const { + param_names__.push_back("theta"); + } - void get_dims(std::vector >& dimss__) const { - dimss__.resize(0); - std::vector scalar_dim; - dimss__.push_back(scalar_dim); - } + void get_param_names(std::vector& names) const { + constrained_param_names(names); + } - void constrained_param_names(std::vector& param_names__, + void unconstrained_param_names(std::vector& param_names__, bool include_tparams__ = true, bool include_gqs__ = true) const { - param_names__.push_back("theta"); - } - - void get_param_names(std::vector& names) const { - constrained_param_names(names); - } - - void unconstrained_param_names(std::vector& param_names__, - bool include_tparams__ = true, - bool include_gqs__ = true) const { - param_names__.clear(); - for (size_t n = 0; n < num_params_r__; ++n) { - std::stringstream param_name; - param_name << "param_" << n; - param_names__.push_back(param_name.str()); - } - } - template - void write_array(RNG& base_rng__, - std::vector& params_r__, - std::vector& params_i__, - std::vector& vars__, - bool include_tparams__ = true, - bool include_gqs__ = true, - std::ostream* pstream__ = 0) const { - ++write_array_calls; - vars__.resize(0); - for (size_t i = 0; i < params_r__.size(); ++i) - vars__.push_back(params_r__[i]); + param_names__.clear(); + for (size_t n = 0; n < num_params_r__; ++n) { + std::stringstream param_name; + param_name << "param_" << n; + param_names__.push_back(param_name.str()); } + } + template + void write_array(RNG& base_rng__, + std::vector& params_r__, + std::vector& params_i__, + std::vector& vars__, + bool include_tparams__ = true, + bool include_gqs__ = true, + std::ostream* pstream__ = 0) const { + ++write_array_calls; + vars__.resize(0); + for (size_t i = 0; i < params_r__.size(); ++i) + vars__.push_back(params_r__[i]); + } + + mutable int templated_log_prob_calls; + mutable int transform_inits_calls; + mutable int write_array_calls; + double log_prob_return_value; +}; - mutable int templated_log_prob_calls; - mutable int transform_inits_calls; - mutable int write_array_calls; - double log_prob_return_value; - }; } TEST_F(ServicesUtilInitialize, model_throws__radius_zero) { @@ -322,88 +323,88 @@ TEST_F(ServicesUtilInitialize, model_throws__full_init) { namespace test { - // Mock Throwing Model throws exception - class mock_error_model: public stan::model::prob_grad { - public: +// Mock Throwing Model throws exception +class mock_error_model: public stan::model::prob_grad { + public: - mock_error_model(): + mock_error_model(): stan::model::prob_grad(1), templated_log_prob_calls(0), transform_inits_calls(0), write_array_calls(0), log_prob_return_value(0.0) { } - void reset() { - templated_log_prob_calls = 0; - transform_inits_calls = 0; - write_array_calls = 0; - log_prob_return_value = 0.0; + void reset() { + templated_log_prob_calls = 0; + transform_inits_calls = 0; + write_array_calls = 0; + log_prob_return_value = 0.0; + } + + template + T__ log_prob(std::vector& params_r__, + std::vector& params_i__, + std::ostream* pstream__ = 0) const { + ++templated_log_prob_calls; + throw std::out_of_range("out_of_range error in log_prob"); + return log_prob_return_value; + } + + void transform_inits(const stan::io::var_context& context__, + std::vector& params_i__, + std::vector& params_r__, + std::ostream* pstream__) const { + ++transform_inits_calls; + for (size_t n = 0; n < params_r__.size(); ++n) { + params_r__[n] = n; } + } - template - T__ log_prob(std::vector& params_r__, - std::vector& params_i__, - std::ostream* pstream__ = 0) const { - ++templated_log_prob_calls; - throw std::out_of_range("out_of_range error in log_prob"); - return log_prob_return_value; - } + void get_dims(std::vector >& dimss__) const { + dimss__.resize(0); + std::vector scalar_dim; + dimss__.push_back(scalar_dim); + } - void transform_inits(const stan::io::var_context& context__, - std::vector& params_i__, - std::vector& params_r__, - std::ostream* pstream__) const { - ++transform_inits_calls; - for (size_t n = 0; n < params_r__.size(); ++n) { - params_r__[n] = n; - } - } + void constrained_param_names(std::vector& param_names__, + bool include_tparams__ = true, + bool include_gqs__ = true) const { + param_names__.push_back("theta"); + } - void get_dims(std::vector >& dimss__) const { - dimss__.resize(0); - std::vector scalar_dim; - dimss__.push_back(scalar_dim); - } + void get_param_names(std::vector& names) const { + constrained_param_names(names); + } - void constrained_param_names(std::vector& param_names__, + void unconstrained_param_names(std::vector& param_names__, bool include_tparams__ = true, bool include_gqs__ = true) const { - param_names__.push_back("theta"); - } - - void get_param_names(std::vector& names) const { - constrained_param_names(names); - } - - void unconstrained_param_names(std::vector& param_names__, - bool include_tparams__ = true, - bool include_gqs__ = true) const { - param_names__.clear(); - for (size_t n = 0; n < num_params_r__; ++n) { - std::stringstream param_name; - param_name << "param_" << n; - param_names__.push_back(param_name.str()); - } + param_names__.clear(); + for (size_t n = 0; n < num_params_r__; ++n) { + std::stringstream param_name; + param_name << "param_" << n; + param_names__.push_back(param_name.str()); } - template - void write_array(RNG& base_rng__, - std::vector& params_r__, - std::vector& params_i__, - std::vector& vars__, - bool include_tparams__ = true, - bool include_gqs__ = true, - std::ostream* pstream__ = 0) const { - ++write_array_calls; - vars__.resize(0); - for (size_t i = 0; i < params_r__.size(); ++i) - vars__.push_back(params_r__[i]); - } - - mutable int templated_log_prob_calls; - mutable int transform_inits_calls; - mutable int write_array_calls; - double log_prob_return_value; - }; + } + template + void write_array(RNG& base_rng__, + std::vector& params_r__, + std::vector& params_i__, + std::vector& vars__, + bool include_tparams__ = true, + bool include_gqs__ = true, + std::ostream* pstream__ = 0) const { + ++write_array_calls; + vars__.resize(0); + for (size_t i = 0; i < params_r__.size(); ++i) + vars__.push_back(params_r__[i]); + } + + mutable int templated_log_prob_calls; + mutable int transform_inits_calls; + mutable int write_array_calls; + double log_prob_return_value; +}; } @@ -463,3 +464,143 @@ TEST_F(ServicesUtilInitialize, model_errors__full_init) { EXPECT_EQ(2, logger.call_count_info()); EXPECT_EQ(1, logger.find_info("out_of_range error in log_prob")); } + +namespace test { +// mock_throwing_model_in_write_array throws exception in the write_array() +// method +class mock_throwing_model_in_write_array: public stan::model::prob_grad { + public: + + mock_throwing_model_in_write_array(): + stan::model::prob_grad(1), + templated_log_prob_calls(0), + transform_inits_calls(0), + write_array_calls(0), + log_prob_return_value(0.0) { } + + void reset() { + templated_log_prob_calls = 0; + transform_inits_calls = 0; + write_array_calls = 0; + log_prob_return_value = 0.0; + } + + template + T__ log_prob(std::vector& params_r__, + std::vector& params_i__, + std::ostream* pstream__ = 0) const { + ++templated_log_prob_calls; + return log_prob_return_value; + } + + void transform_inits(const stan::io::var_context& context__, + std::vector& params_i__, + std::vector& params_r__, + std::ostream* pstream__) const { + ++transform_inits_calls; + for (size_t n = 0; n < params_r__.size(); ++n) { + params_r__[n] = n; + } + } + + void get_dims(std::vector >& dimss__) const { + dimss__.resize(0); + std::vector scalar_dim; + dimss__.push_back(scalar_dim); + } + + void constrained_param_names(std::vector& param_names__, + bool include_tparams__ = true, + bool include_gqs__ = true) const { + param_names__.push_back("theta"); + } + + void get_param_names(std::vector& names) const { + constrained_param_names(names); + } + + void unconstrained_param_names(std::vector& param_names__, + bool include_tparams__ = true, + bool include_gqs__ = true) const { + param_names__.clear(); + for (size_t n = 0; n < num_params_r__; ++n) { + std::stringstream param_name; + param_name << "param_" << n; + param_names__.push_back(param_name.str()); + } + } + template + void write_array(RNG& base_rng__, + std::vector& params_r__, + std::vector& params_i__, + std::vector& vars__, + bool include_tparams__ = true, + bool include_gqs__ = true, + std::ostream* pstream__ = 0) const { + ++write_array_calls; + throw std::domain_error("throwing within write_array"); + vars__.resize(0); + for (size_t i = 0; i < params_r__.size(); ++i) + vars__.push_back(params_r__[i]); + } + + mutable int templated_log_prob_calls; + mutable int transform_inits_calls; + mutable int write_array_calls; + double log_prob_return_value; +}; +} + +TEST_F(ServicesUtilInitialize, model_throws_in_write_array__radius_zero) { + test::mock_throwing_model_in_write_array throwing_model; + + double init_radius = 0; + bool print_timing = false; + EXPECT_THROW(stan::services::util::initialize(throwing_model, empty_context, rng, + init_radius, print_timing, + logger, init), + std::domain_error); + + EXPECT_EQ(3, logger.call_count()); + EXPECT_EQ(3, logger.call_count_info()); + EXPECT_EQ(1, logger.find_info("throwing within write_array")); +} + +TEST_F(ServicesUtilInitialize, model_throws_in_write_array__radius_two) { + test::mock_throwing_model_in_write_array throwing_model; + + double init_radius = 2; + bool print_timing = false; + EXPECT_THROW(stan::services::util::initialize(throwing_model, empty_context, rng, + init_radius, print_timing, + logger, init), + std::domain_error); + EXPECT_EQ(303, logger.call_count()); + EXPECT_EQ(303, logger.call_count_info()); + EXPECT_EQ(100, logger.find_info("throwing within write_array")); +} + +TEST_F(ServicesUtilInitialize, model_throws_in_write_array__full_init) { + std::vector names_r; + std::vector values_r; + std::vector > dim_r; + names_r.push_back("y"); + values_r.push_back(6.35149); // 1.5 unconstrained: -10 + 20 * inv.logit(1.5) + values_r.push_back(-2.449187); // -0.5 unconstrained + std::vector d; + d.push_back(2); + dim_r.push_back(d); + stan::io::array_var_context init_context(names_r, values_r, dim_r); + + test::mock_throwing_model_in_write_array throwing_model; + + double init_radius = 2; + bool print_timing = false; + EXPECT_THROW(stan::services::util::initialize(throwing_model, init_context, rng, + init_radius, print_timing, + logger, init), + std::domain_error); + EXPECT_EQ(303, logger.call_count()); + EXPECT_EQ(303, logger.call_count_info()); + EXPECT_EQ(100, logger.find_info("throwing within write_array")); +} From 0e26281a4c761d2088a4d14ee2e0534b81d5d446 Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Fri, 6 Jul 2018 16:25:07 -0400 Subject: [PATCH 4/9] Adding additional test for random_var_context that shows the current behavior --- src/test/unit/io/random_var_context_test.cpp | 104 ++++++++++++++++++- 1 file changed, 100 insertions(+), 4 deletions(-) diff --git a/src/test/unit/io/random_var_context_test.cpp b/src/test/unit/io/random_var_context_test.cpp index 5311da855de..60e9b308c23 100644 --- a/src/test/unit/io/random_var_context_test.cpp +++ b/src/test/unit/io/random_var_context_test.cpp @@ -4,17 +4,107 @@ #include // L'Ecuyer RNG #include #include +#include + +namespace test { +// mock_throwing_model_in_write_array throws exception in the write_array() +// method +class mock_throwing_model_in_write_array: public stan::model::prob_grad { + public: + + mock_throwing_model_in_write_array(): + stan::model::prob_grad(1), + templated_log_prob_calls(0), + transform_inits_calls(0), + write_array_calls(0), + log_prob_return_value(0.0) { } + + void reset() { + templated_log_prob_calls = 0; + transform_inits_calls = 0; + write_array_calls = 0; + log_prob_return_value = 0.0; + } + + template + T__ log_prob(std::vector& params_r__, + std::vector& params_i__, + std::ostream* pstream__ = 0) const { + ++templated_log_prob_calls; + return log_prob_return_value; + } + + void transform_inits(const stan::io::var_context& context__, + std::vector& params_i__, + std::vector& params_r__, + std::ostream* pstream__) const { + ++transform_inits_calls; + for (size_t n = 0; n < params_r__.size(); ++n) { + params_r__[n] = n; + } + } + + void get_dims(std::vector >& dimss__) const { + dimss__.resize(0); + std::vector scalar_dim; + dimss__.push_back(scalar_dim); + } + + void constrained_param_names(std::vector& param_names__, + bool include_tparams__ = true, + bool include_gqs__ = true) const { + param_names__.push_back("theta"); + } + + void get_param_names(std::vector& names) const { + constrained_param_names(names); + } + + void unconstrained_param_names(std::vector& param_names__, + bool include_tparams__ = true, + bool include_gqs__ = true) const { + param_names__.clear(); + for (size_t n = 0; n < num_params_r__; ++n) { + std::stringstream param_name; + param_name << "param_" << n; + param_names__.push_back(param_name.str()); + } + } + template + void write_array(RNG& base_rng__, + std::vector& params_r__, + std::vector& params_i__, + std::vector& vars__, + bool include_tparams__ = true, + bool include_gqs__ = true, + std::ostream* pstream__ = 0) const { + ++write_array_calls; + throw std::domain_error("throwing within write_array"); + vars__.resize(0); + for (size_t i = 0; i < params_r__.size(); ++i) + vars__.push_back(params_r__[i]); + } + + mutable int templated_log_prob_calls; + mutable int transform_inits_calls; + mutable int write_array_calls; + double log_prob_return_value; +}; +} + class random_var_context : public testing::Test { -public: + public: random_var_context() - : empty_context(), - model(empty_context, static_cast(0)), - rng(0) { } + : empty_context(), + model(empty_context, static_cast(0)), + rng(0), + throwing_model() { } stan::io::empty_var_context empty_context; stan_model model; boost::ecuyer1988 rng; + test::mock_throwing_model_in_write_array throwing_model; }; TEST_F(random_var_context, contains_r) { @@ -81,3 +171,9 @@ TEST_F(random_var_context, names_i) { EXPECT_NO_THROW(context.names_i(names_i)); EXPECT_EQ(0, names_i.size()); } + +TEST_F(random_var_context, construct) { + EXPECT_THROW_MSG(stan::io::random_var_context(throwing_model, rng, 2, false), + std::domain_error, + "throwing within write_array"); +} From b456e0e5e21050b5b3482075387091e7d3e22034 Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Fri, 6 Jul 2018 20:46:20 -0400 Subject: [PATCH 5/9] Fixes #2562. Updating src/stan/services/util/initialize.hpp to handle exception behavior in new generated code. --- src/stan/services/util/initialize.hpp | 37 +++++++++++++++++---------- 1 file changed, 23 insertions(+), 14 deletions(-) diff --git a/src/stan/services/util/initialize.hpp b/src/stan/services/util/initialize.hpp index 383dd28533d..34e6917af67 100644 --- a/src/stan/services/util/initialize.hpp +++ b/src/stan/services/util/initialize.hpp @@ -85,31 +85,40 @@ std::vector initialize(Model& model, int MAX_INIT_TRIES = is_fully_initialized || init_zero ? 1 : 100; int num_init_tries = 0; for (; num_init_tries < MAX_INIT_TRIES; num_init_tries++) { - stan::io::random_var_context - random_context(model, rng, init_radius, init_zero); + std::stringstream msg; + try { + stan::io::random_var_context + random_context(model, rng, init_radius, init_zero); - if (!any_initialized) { - unconstrained = random_context.get_unconstrained(); - } else { - stan::io::chained_var_context context(init, random_context); + if (!any_initialized) { + unconstrained = random_context.get_unconstrained(); + } else { + stan::io::chained_var_context context(init, random_context); - std::stringstream msg; - try { model.transform_inits(context, disc_vector, unconstrained, &msg); - } catch (const std::exception& e) { - if (msg.str().length() > 0) - logger.info(msg); - logger.info(e.what()); - throw; } + } catch (std::domain_error& e) { + if (msg.str().length() > 0) + logger.info(msg); + logger.info("Rejecting initial value:"); + logger.info(" Error evaluating the log probability" + " at the initial value."); + logger.info(e.what()); + continue; + } catch (std::exception& e) { if (msg.str().length() > 0) logger.info(msg); + logger.info("Unrecoverable error evaluating the log probability" + " at the initial value."); + logger.info(e.what()); + throw; } + + msg.str(""); double log_prob(0); - std::stringstream msg; try { log_prob = model.template log_prob (unconstrained, disc_vector, &msg); From 813c89e94e1a0863ddfcf889908bee8f34912512 Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Tue, 10 Jul 2018 10:14:25 -0400 Subject: [PATCH 6/9] adding newline --- .gitignore | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index 3dde639abb1..7ce6ffac48b 100644 --- a/.gitignore +++ b/.gitignore @@ -57,4 +57,4 @@ stan.kdev4 # local make include /make/local -*.gch \ No newline at end of file +*.gch From 055f60f19c9abbccaec83a49cb9a915ada6974e5 Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Tue, 10 Jul 2018 10:14:33 -0400 Subject: [PATCH 7/9] updating initialization code --- src/stan/services/optimize/bfgs.hpp | 4 ++-- src/stan/services/optimize/lbfgs.hpp | 4 ++-- src/stan/services/optimize/newton.hpp | 4 ++-- src/stan/services/util/initialize.hpp | 29 ++++++++++++++++++--------- 4 files changed, 25 insertions(+), 16 deletions(-) diff --git a/src/stan/services/optimize/bfgs.hpp b/src/stan/services/optimize/bfgs.hpp index d14b0639351..5829ce718ca 100644 --- a/src/stan/services/optimize/bfgs.hpp +++ b/src/stan/services/optimize/bfgs.hpp @@ -65,8 +65,8 @@ namespace stan { std::vector disc_vector; std::vector cont_vector - = util::initialize(model, init, rng, init_radius, false, - logger, init_writer); + = util::initialize(model, init, rng, init_radius, false, + logger, init_writer); std::stringstream bfgs_ss; typedef stan::optimization::BFGSLineSearch diff --git a/src/stan/services/optimize/lbfgs.hpp b/src/stan/services/optimize/lbfgs.hpp index 48a307cf70c..042bd91730a 100644 --- a/src/stan/services/optimize/lbfgs.hpp +++ b/src/stan/services/optimize/lbfgs.hpp @@ -66,8 +66,8 @@ namespace stan { std::vector disc_vector; std::vector cont_vector - = util::initialize(model, init, rng, init_radius, false, - logger, init_writer); + = util::initialize(model, init, rng, init_radius, false, + logger, init_writer); std::stringstream lbfgs_ss; typedef stan::optimization::BFGSLineSearch diff --git a/src/stan/services/optimize/newton.hpp b/src/stan/services/optimize/newton.hpp index 07611ebbc37..573484c9c5f 100644 --- a/src/stan/services/optimize/newton.hpp +++ b/src/stan/services/optimize/newton.hpp @@ -51,8 +51,8 @@ namespace stan { std::vector disc_vector; std::vector cont_vector - = util::initialize(model, init, rng, init_radius, false, - logger, init_writer); + = util::initialize(model, init, rng, init_radius, false, + logger, init_writer); double lp(0); diff --git a/src/stan/services/util/initialize.hpp b/src/stan/services/util/initialize.hpp index 34e6917af67..cbf32c8e1d2 100644 --- a/src/stan/services/util/initialize.hpp +++ b/src/stan/services/util/initialize.hpp @@ -7,7 +7,6 @@ #include #include #include -#include #include #include #include @@ -38,10 +37,15 @@ namespace util { * parameters that are valid or it hits MAX_INIT_TRIES = * 100 (hard-coded). * - * Valid initialization is defined as a finite, non-NaN value - * for the evaluation of the log probability and all its + * Valid initialization is defined as a finite, non-NaN value for the + * evaluation of the log probability density function and all its * gradients. * + * @tparam Jacobian indicates whether to include the Jacobian term when + * evaluating the log density function + * @tparam Model the type of the model class + * @tparam RNG the type of the random number generator + * * @param[in] model the model * @param[in] init a var_context with initial values * @param[in,out] rng random number generator @@ -59,7 +63,7 @@ namespace util { * std::domain_error) * @return valid unconstrained parameters for the model */ -template +template std::vector initialize(Model& model, stan::io::var_context& init, RNG& rng, @@ -80,15 +84,15 @@ std::vector initialize(Model& model, any_initialized |= init.contains_r(param_names[n]); } - bool init_zero = init_radius <= std::numeric_limits::min(); + bool is_initialized_with_zero = init_radius == 0.0; - int MAX_INIT_TRIES = is_fully_initialized || init_zero ? 1 : 100; + int MAX_INIT_TRIES = is_fully_initialized || is_initialized_with_zero ? 1 : 100; int num_init_tries = 0; for (; num_init_tries < MAX_INIT_TRIES; num_init_tries++) { std::stringstream msg; try { stan::io::random_var_context - random_context(model, rng, init_radius, init_zero); + random_context(model, rng, init_radius, is_initialized_with_zero); if (!any_initialized) { unconstrained = random_context.get_unconstrained(); @@ -120,7 +124,10 @@ std::vector initialize(Model& model, msg.str(""); double log_prob(0); try { - log_prob = model.template log_prob + // we evaluate the log_prob function with propto=false + // because we're evaluating with `double` as the type of + // the parameters. + log_prob = model.template log_prob (unconstrained, disc_vector, &msg); if (msg.str().length() > 0) logger.info(msg); @@ -153,7 +160,9 @@ std::vector initialize(Model& model, bool gradient_ok = true; clock_t start_check = clock(); try { - log_prob = stan::model::log_prob_grad + // we evaluate this with propto=true since we're + // evaluating with autodiff variables + log_prob = stan::model::log_prob_grad (model, unconstrained, disc_vector, gradient, &log_prob_msg); } catch (const std::exception& e) { @@ -199,7 +208,7 @@ std::vector initialize(Model& model, } if (num_init_tries == MAX_INIT_TRIES) { - if (!is_fully_initialized && !init_zero) { + if (!is_fully_initialized && !is_initialized_with_zero) { logger.info(""); std::stringstream msg; msg << "Initialization between (-" << init_radius From cacf1ff470761a51268f82dedadc948d5d50783d Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Tue, 10 Jul 2018 10:52:00 -0400 Subject: [PATCH 8/9] simplifying logic --- src/stan/services/util/initialize.hpp | 52 +++++++++---------- .../unit/services/util/initialize_test.cpp | 14 +++-- 2 files changed, 30 insertions(+), 36 deletions(-) diff --git a/src/stan/services/util/initialize.hpp b/src/stan/services/util/initialize.hpp index cbf32c8e1d2..cfaf910f038 100644 --- a/src/stan/services/util/initialize.hpp +++ b/src/stan/services/util/initialize.hpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include #include @@ -157,7 +158,6 @@ std::vector initialize(Model& model, } std::stringstream log_prob_msg; std::vector gradient; - bool gradient_ok = true; clock_t start_check = clock(); try { // we evaluate this with propto=true since we're @@ -177,15 +177,14 @@ std::vector initialize(Model& model, if (log_prob_msg.str().length() > 0) logger.info(log_prob_msg); - for (size_t i = 0; i < gradient.size(); ++i) { - if (gradient_ok && !boost::math::isfinite(gradient[i])) { - logger.info("Rejecting initial value:"); - logger.info(" Gradient evaluated at the initial value" - " is not finite."); - logger.info(" Stan can't start sampling from this" - " initial value."); - gradient_ok = false; - } + bool gradient_ok = boost::math::isfinite(stan::math::sum(gradient)); + + if (!gradient_ok) { + logger.info("Rejecting initial value:"); + logger.info(" Gradient evaluated at the initial value" + " is not finite."); + logger.info(" Stan can't start sampling from this" + " initial value."); } if (gradient_ok && print_timing) { logger.info(""); @@ -203,27 +202,24 @@ std::vector initialize(Model& model, logger.info(""); logger.info(""); } - if (gradient_ok) - break; - } - - if (num_init_tries == MAX_INIT_TRIES) { - if (!is_fully_initialized && !is_initialized_with_zero) { - logger.info(""); - std::stringstream msg; - msg << "Initialization between (-" << init_radius - << ", " << init_radius << ") failed after" - << " " << MAX_INIT_TRIES << " attempts. "; - logger.info(msg); - logger.info(" Try specifying initial values," - " reducing ranges of constrained values," - " or reparameterizing the model."); + if (gradient_ok) { + init_writer(unconstrained); + return unconstrained; } - throw std::domain_error("Initialization failed."); } - init_writer(unconstrained); - return unconstrained; + if (!is_initialized_with_zero) { + logger.info(""); + std::stringstream msg; + msg << "Initialization between (-" << init_radius + << ", " << init_radius << ") failed after" + << " " << MAX_INIT_TRIES << " attempts. "; + logger.info(msg); + logger.info(" Try specifying initial values," + " reducing ranges of constrained values," + " or reparameterizing the model."); + } + throw std::domain_error("Initialization failed."); } } diff --git a/src/test/unit/services/util/initialize_test.cpp b/src/test/unit/services/util/initialize_test.cpp index 736e604f3a7..522f31aaf50 100644 --- a/src/test/unit/services/util/initialize_test.cpp +++ b/src/test/unit/services/util/initialize_test.cpp @@ -47,22 +47,20 @@ TEST_F(ServicesUtilInitialize, radius_zero__print_false) { EXPECT_EQ(params[1], init.vector_double_values()[0][1]); } -TEST_F(ServicesUtilInitialize, radius_zero__print_true) { +TEST_F(ServicesUtilInitialize, radius_zero__initialize_with_Jacobian) { std::vector params; double init_radius = 0; - bool print_timing = true; - params = stan::services::util::initialize(model, empty_context, rng, - init_radius, print_timing, - logger, init); + bool print_timing = false; + params = stan::services::util::initialize(model, empty_context, rng, + init_radius, print_timing, + logger, init); ASSERT_EQ(model.num_params_r(), params.size()) << "2 parameters"; EXPECT_FLOAT_EQ(0, params[0]); EXPECT_FLOAT_EQ(0, params[1]); - EXPECT_EQ(6, logger.call_count()); - EXPECT_EQ(6, logger.call_count_info()); - EXPECT_EQ(1, logger.find_info("Gradient evaluation")); + EXPECT_EQ(0, logger.call_count()); ASSERT_EQ(1, init.vector_double_values().size()); ASSERT_EQ(2, init.vector_double_values()[0].size()); EXPECT_EQ(params[0], init.vector_double_values()[0][0]); From 83b3f78e893a0e58859fe58e96cb9fe548f882bd Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Tue, 10 Jul 2018 11:06:07 -0400 Subject: [PATCH 9/9] cpplint --- src/stan/services/util/initialize.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/stan/services/util/initialize.hpp b/src/stan/services/util/initialize.hpp index cfaf910f038..a971afef2cc 100644 --- a/src/stan/services/util/initialize.hpp +++ b/src/stan/services/util/initialize.hpp @@ -87,7 +87,8 @@ std::vector initialize(Model& model, bool is_initialized_with_zero = init_radius == 0.0; - int MAX_INIT_TRIES = is_fully_initialized || is_initialized_with_zero ? 1 : 100; + int MAX_INIT_TRIES = is_fully_initialized || is_initialized_with_zero + ? 1 : 100; int num_init_tries = 0; for (; num_init_tries < MAX_INIT_TRIES; num_init_tries++) { std::stringstream msg;