Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
40 changes: 17 additions & 23 deletions src/stan/services/optimize/laplace_sample.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -61,67 +61,62 @@ 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
Eigen::MatrixXd hessian;
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<char>::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();
Eigen::MatrixXd inv_sqrt_neg_hessian = L_neg_hessian.inverse().transpose();
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);
Comment thread
mitzimorris marked this conversation as resolved.
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<char>::eof())
logger.info(write_array_msgs);
// output draw, log_p, log_q
std::vector<double> 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<double> 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

/**
Expand Down Expand Up @@ -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;
Expand Down
6 changes: 6 additions & 0 deletions src/test/test-models/good/optimization/constrain_sigma.stan
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
parameters {
real<lower=0> sigma;
}
model {
sigma ~ normal(3, 1);
}
75 changes: 75 additions & 0 deletions src/test/unit/services/optimize/laplace_jacobian_test.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,75 @@
#include <boost/algorithm/string.hpp>
#include <gtest/gtest.h>
#include <stan/callbacks/stream_logger.hpp>
#include <stan/callbacks/stream_writer.hpp>
#include <stan/io/dump.hpp>
#include <stan/io/empty_var_context.hpp>
#include <stan/io/stan_csv_reader.hpp>
#include <stan/services/error_codes.hpp>
#include <stan/services/optimize/laplace_sample.hpp>
#include <test/test-models/good/optimization/constrain_sigma.hpp>
#include <test/unit/services/instrumented_callbacks.hpp>
#include <test/unit/util.hpp>
#include <cmath>
#include <iostream>
#include <vector>

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<true>(
*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<false>(*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));
}
26 changes: 25 additions & 1 deletion src/test/unit/services/optimize/laplace_sample_test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<true>(
Expand Down Expand Up @@ -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<true>(
*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));
}