diff --git a/src/stan/services/optimize/laplace_sample.hpp b/src/stan/services/optimize/laplace_sample.hpp index 81755bcdfe1..ae1d8b3f452 100644 --- a/src/stan/services/optimize/laplace_sample.hpp +++ b/src/stan/services/optimize/laplace_sample.hpp @@ -61,7 +61,7 @@ void laplace_sample(const Model& model, const Eigen::VectorXd& theta_hat, // calculate inverse negative Hessian's Cholesky factor if (refresh > 0) { - logger.info("Calculating Hessian\n"); + logger.info("Calculating Hessian"); } double log_p; // dummy Eigen::VectorXd grad; // dummy @@ -69,16 +69,13 @@ void laplace_sample(const Model& model, const Eigen::VectorXd& theta_hat, interrupt(); math::internal::finite_diff_hessian_auto(log_density_fun, theta_hat, log_p, grad, hessian); - if (refresh > 0) { + if (refresh > 0 && log_density_msgs.peek() != std::char_traits::eof()) logger.info(log_density_msgs); - logger.info("\n"); - } // calculate Cholesky factor and inverse interrupt(); if (refresh > 0) { logger.info("Calculating inverse of Cholesky factor"); - logger.info("\n"); } Eigen::MatrixXd L_neg_hessian = (-hessian).llt().matrixL(); interrupt(); @@ -86,42 +83,40 @@ void laplace_sample(const Model& model, const Eigen::VectorXd& theta_hat, interrupt(); Eigen::MatrixXd half_hessian = 0.5 * hessian; - // generate draws and output to sample writer if (refresh > 0) { logger.info("Generating draws"); - logger.info("\n"); } + // generate draws + std::stringstream refresh_msg; boost::ecuyer1988 rng = util::create_rng(random_seed, 0); Eigen::VectorXd draw_vec; // declare draw_vec, msgs here to avoid re-alloc for (int m = 0; m < draws; ++m) { interrupt(); // allow interpution each iteration if (refresh > 0 && m % refresh == 0) { - logger.info("iteration: "); - logger.info(std::to_string(m)); - logger.info("\n"); + refresh_msg << "iteration: " << std::to_string(m); + logger.info(refresh_msg); + refresh_msg.str(std::string()); } Eigen::VectorXd z(num_unc_params); for (int n = 0; n < num_unc_params; ++n) { z(n) = math::std_normal_rng(rng); } - Eigen::VectorXd unc_draw = theta_hat + inv_sqrt_neg_hessian * z; + std::stringstream write_array_msgs; + model.write_array(rng, unc_draw, draw_vec, include_tp, include_gq, + &write_array_msgs); + if (refresh > 0 && write_array_msgs.peek() != std::char_traits::eof()) + logger.info(write_array_msgs); + // output draw, log_p, log_q + std::vector draw(&draw_vec(0), &draw_vec(0) + draw_size); double log_p = log_density_fun(unc_draw).val(); + draw.push_back(log_p); Eigen::VectorXd diff = unc_draw - theta_hat; double log_q = diff.transpose() * half_hessian * diff; - - std::stringstream msgs; - model.write_array(rng, unc_draw, draw_vec, include_tp, include_gq, &msgs); - if (refresh > 0) { - logger.info(msgs); - logger.info("\n"); - } - std::vector draw(&draw_vec(0), &draw_vec(0) + draw_size); - draw.push_back(log_p); draw.push_back(log_q); sample_writer(draw); } -} +} // namespace internal } // namespace internal /** @@ -168,11 +163,10 @@ int laplace_sample(const Model& model, const Eigen::VectorXd& theta_hat, } catch (const std::exception& e) { if (refresh >= 0) { logger.error(e.what()); - logger.error("\n"); } } catch (...) { if (refresh >= 0) { - logger.error("unknown exception during execution\n"); + logger.error("unknown exception during execution"); } } return error_codes::DATAERR; diff --git a/src/test/test-models/good/optimization/constrain_sigma.stan b/src/test/test-models/good/optimization/constrain_sigma.stan new file mode 100644 index 00000000000..cf6ce126096 --- /dev/null +++ b/src/test/test-models/good/optimization/constrain_sigma.stan @@ -0,0 +1,6 @@ +parameters { + real sigma; +} +model { + sigma ~ normal(3, 1); +} diff --git a/src/test/unit/services/optimize/laplace_jacobian_test.cpp b/src/test/unit/services/optimize/laplace_jacobian_test.cpp new file mode 100644 index 00000000000..105058d5d7e --- /dev/null +++ b/src/test/unit/services/optimize/laplace_jacobian_test.cpp @@ -0,0 +1,75 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +class ServicesLaplaceJacobian : public ::testing::Test { + public: + ServicesLaplaceJacobian() : logger(msgs, msgs, msgs, msgs, msgs) {} + + void SetUp() { + stan::io::empty_var_context var_context; + model = new stan_model(var_context); // typedef from multi_normal.hpp + } + + void TearDown() { delete model; } + + stan_model* model; + std::stringstream msgs; + stan::callbacks::stream_logger logger; + stan::test::unit::instrumented_interrupt interrupt; +}; + +TEST_F(ServicesLaplaceJacobian, laplace_jacobian_adjust) { + Eigen::VectorXd theta_hat(1); + theta_hat << 1.09; // test for mode sigma at 3.1, take log + int draws = 5; + unsigned int seed = 1234; + int refresh = 1; + std::stringstream sample1_ss; + stan::callbacks::stream_writer sample_writer(sample1_ss, ""); + int return_code = stan::services::laplace_sample( + *model, theta_hat, draws, seed, refresh, interrupt, logger, + sample_writer); + EXPECT_EQ(stan::services::error_codes::OK, return_code); + std::stringstream out; + stan::io::stan_csv draws1 + = stan::io::stan_csv_reader::parse(sample1_ss, &out); + + std::stringstream sample2_ss; + stan::callbacks::stream_writer sample_writer2(sample2_ss, ""); + return_code = stan::services::laplace_sample(*model, theta_hat, draws, + seed, refresh, interrupt, + logger, sample_writer2); + EXPECT_EQ(stan::services::error_codes::OK, return_code); + stan::io::stan_csv draws2 + = stan::io::stan_csv_reader::parse(sample2_ss, &out); + + EXPECT_EQ(3, draws1.header.size()); + EXPECT_EQ("sigma", draws1.header[0]); + EXPECT_EQ("log_p", draws1.header[1]); + EXPECT_EQ("log_q", draws1.header[2]); + + EXPECT_EQ(3, draws2.header.size()); + EXPECT_EQ("sigma", draws2.header[0]); + EXPECT_EQ("log_p", draws2.header[1]); + EXPECT_EQ("log_q", draws2.header[2]); + + Eigen::MatrixXd sample1 = draws1.samples; + Eigen::MatrixXd sample2 = draws2.samples; + + EXPECT_EQ(sample1.coeff(0, 0), sample2.coeff(0, 0)); + EXPECT_NE(sample1.coeff(0, 1), sample2.coeff(0, 1)); + EXPECT_EQ(sample1.coeff(0, 2), sample2.coeff(0, 2)); +} diff --git a/src/test/unit/services/optimize/laplace_sample_test.cpp b/src/test/unit/services/optimize/laplace_sample_test.cpp index f5467bb842f..3737f0d8a88 100644 --- a/src/test/unit/services/optimize/laplace_sample_test.cpp +++ b/src/test/unit/services/optimize/laplace_sample_test.cpp @@ -36,7 +36,7 @@ TEST_F(ServicesLaplaceSample, values) { theta_hat << 2, 3; int draws = 50000; // big to enable mean, var & covar test precision unsigned int seed = 1234; - int refresh = 1; + int refresh = 100; std::stringstream sample_ss; stan::callbacks::stream_writer sample_writer(sample_ss, ""); int return_code = stan::services::laplace_sample( @@ -126,3 +126,27 @@ TEST_F(ServicesLaplaceSample, nonPositiveDrawsError) { sample_writer); EXPECT_EQ(stan::services::error_codes::DATAERR, RC); } + +TEST_F(ServicesLaplaceSample, consoleOutput) { + Eigen::VectorXd theta_hat(2); + theta_hat << 2, 3; + int draws = 10; + unsigned int seed = 1234; + int refresh = 1; + std::stringstream sample_ss; + stan::callbacks::stream_writer sample_writer(sample_ss, ""); + std::stringstream logger_ss; + stan::callbacks::stream_logger sample_logger(logger_ss, logger_ss, logger_ss, + logger_ss, logger_ss); + int return_code = stan::services::laplace_sample( + *model, theta_hat, draws, seed, refresh, interrupt, sample_logger, + sample_writer); + EXPECT_EQ(stan::services::error_codes::OK, return_code); + std::string console_str = logger_ss.str(); + EXPECT_EQ(1, + count_matches( + "Calculating Hessian\nCalculating inverse of Cholesky factor\n", + console_str)); + EXPECT_EQ(1, count_matches("Generating draws\niteration: 0\niteration: 1", + console_str)); +}