Skip to content

Commit 43f3fea

Browse files
author
Nicusor Serban
committed
Merge branch 'develop' into release/v2.36.0
2 parents 13b0c5c + f4fb7cb commit 43f3fea

File tree

6 files changed

+334
-12
lines changed

6 files changed

+334
-12
lines changed

‎src/stan/io/stan_csv_reader.hpp‎

Lines changed: 9 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -250,6 +250,10 @@ class stan_csv_reader {
250250

251251
int rows = lines - 3;
252252
int cols = std::count(line.begin(), line.end(), ',') + 1;
253+
if (cols == 1) {
254+
// model has no parameters
255+
return;
256+
}
253257
adaptation.metric.resize(rows, cols);
254258
char comment; // Buffer for comment indicator, #
255259

@@ -330,7 +334,11 @@ class stan_csv_reader {
330334
for (int col = 0; col < cols; col++) {
331335
std::getline(ls, line, ',');
332336
boost::trim(line);
333-
std::stringstream(line) >> samples(row, col);
337+
try {
338+
samples(row, col) = static_cast<double>(std::stold(line));
339+
} catch (const std::out_of_range& e) {
340+
samples(row, col) = std::numeric_limits<double>::quiet_NaN();
341+
}
334342
}
335343
}
336344
}

‎src/stan/mcmc/chainset.hpp‎

Lines changed: 32 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -253,7 +253,7 @@ class chainset {
253253
double median(const std::string& name) const { return median(index(name)); }
254254

255255
/**
256-
* Compute maximum absolute deviation (mad) for specified parameter.
256+
* Compute median absolute deviation (mad) for specified parameter.
257257
*
258258
* Follows R implementation: constant * median(abs(x - center))
259259
* where the value of center is median(x) and the constant is 1.4826,
@@ -263,16 +263,20 @@ class chainset {
263263
* @param index parameter index
264264
* @return sample mad
265265
*/
266-
double max_abs_deviation(const int index) const {
266+
double med_abs_deviation(const int index) const {
267267
Eigen::MatrixXd draws = samples(index);
268268
auto center = median(index);
269269
Eigen::MatrixXd abs_dev = (draws.array() - center).abs();
270270
Eigen::Map<Eigen::VectorXd> map(abs_dev.data(), abs_dev.size());
271-
return 1.4826 * stan::math::quantile(map, 0.5);
271+
try {
272+
return 1.4826 * stan::math::quantile(map, 0.5);
273+
} catch (const std::logic_error& e) {
274+
return std::numeric_limits<double>::quiet_NaN();
275+
}
272276
}
273277

274278
/**
275-
* Compute maximum absolute deviation (mad) for specified parameter.
279+
* Compute median absolute deviation (mad) for specified parameter.
276280
*
277281
* Follows R implementation: constant * median(abs(x - center))
278282
* where the value of center is median(x) and the constant is 1.4826,
@@ -282,27 +286,33 @@ class chainset {
282286
* @param name parameter name
283287
* @return sample mad
284288
*/
285-
double max_abs_deviation(const std::string& name) const {
286-
return max_abs_deviation(index(name));
289+
double med_abs_deviation(const std::string& name) const {
290+
return med_abs_deviation(index(name));
287291
}
288292

289293
/**
290294
* Compute the quantile value of the specified parameter
291-
* at the specified probability.
295+
* at the specified probability via a call to stan::math::quantile.
292296
*
293297
* Calls stan::math::quantile which throws
294298
* std::invalid_argument If any element of samples_vec is NaN, or size 0.
295299
* and std::domain_error If `p<0` or `p>1`.
300+
* If this happens, error will be caught, quantile value is NaN.
296301
*
297302
* @param index parameter index
298303
* @param prob probability
299304
* @return parameter value at quantile
300305
*/
301306
double quantile(const int index, const double prob) const {
302-
// Ensure the probability is within [0, 1]
303307
Eigen::MatrixXd draws = samples(index);
304308
Eigen::Map<Eigen::VectorXd> map(draws.data(), draws.size());
305-
return stan::math::quantile(map, prob);
309+
double result;
310+
try {
311+
result = stan::math::quantile(map, prob);
312+
} catch (const std::logic_error& e) {
313+
return std::numeric_limits<double>::quiet_NaN();
314+
}
315+
return result;
306316
}
307317

308318
/**
@@ -321,6 +331,11 @@ class chainset {
321331
* Compute the quantile values of the specified parameter
322332
* for a set of specified probabilities.
323333
*
334+
* Calls stan::math::quantile which throws
335+
* std::invalid_argument If any element of samples_vec is NaN, or size 0.
336+
* and std::domain_error If `p<0` or `p>1`.
337+
* If this happens, error will be caught, quantile value is NaN.
338+
*
324339
* @param index parameter index
325340
* @param probs vector of probabilities
326341
* @return vector of parameter values for quantiles
@@ -332,7 +347,14 @@ class chainset {
332347
Eigen::MatrixXd draws = samples(index);
333348
Eigen::Map<Eigen::VectorXd> map(draws.data(), draws.size());
334349
std::vector<double> probs_vec(probs.data(), probs.data() + probs.size());
335-
std::vector<double> quantiles = stan::math::quantile(map, probs_vec);
350+
std::vector<double> quantiles;
351+
try {
352+
quantiles = stan::math::quantile(map, probs_vec);
353+
} catch (const std::logic_error& e) {
354+
Eigen::VectorXd nans(probs.size());
355+
nans.setConstant(std::numeric_limits<double>::quiet_NaN());
356+
return nans;
357+
}
336358
return Eigen::Map<Eigen::VectorXd>(quantiles.data(), quantiles.size());
337359
}
338360

‎src/test/unit/io/stan_csv_reader_test.cpp‎

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,8 @@ class StanIoStanCsvReader : public testing::Test {
3030
"src/test/unit/io/test_csv_files/bernoulli_warmup.csv");
3131
fixed_param_stream.open(
3232
"src/test/unit/io/test_csv_files/fixed_param_output.csv");
33+
no_params_stream.open(
34+
"src/test/unit/io/test_csv_files/no_parameters_hmc.csv");
3335
}
3436

3537
void TearDown() {
@@ -46,6 +48,7 @@ class StanIoStanCsvReader : public testing::Test {
4648
bernoulli_thin_stream.close();
4749
bernoulli_warmup_stream.close();
4850
fixed_param_stream.close();
51+
no_params_stream.close();
4952
}
5053

5154
std::ifstream blocker0_stream, epil0_stream;
@@ -58,6 +61,7 @@ class StanIoStanCsvReader : public testing::Test {
5861
std::ifstream bernoulli_thin_stream;
5962
std::ifstream bernoulli_warmup_stream;
6063
std::ifstream fixed_param_stream;
64+
std::ifstream no_params_stream;
6165
};
6266

6367
TEST_F(StanIoStanCsvReader, read_metadata1) {
@@ -570,6 +574,14 @@ TEST_F(StanIoStanCsvReader, fixed_param) {
570574
ASSERT_EQ(10, fixed_param.samples.rows());
571575
}
572576

577+
TEST_F(StanIoStanCsvReader, no_parameters) {
578+
stan::io::stan_csv no_parameters_hmc;
579+
std::stringstream out;
580+
no_parameters_hmc = stan::io::stan_csv_reader::parse(no_params_stream, &out);
581+
ASSERT_EQ(100, no_parameters_hmc.samples.rows());
582+
ASSERT_EQ(0, no_parameters_hmc.adaptation.metric.size());
583+
}
584+
573585
TEST_F(StanIoStanCsvReader, no_samples) {
574586
std::ifstream no_samples_stream;
575587
no_samples_stream.open(
@@ -590,4 +602,16 @@ TEST_F(StanIoStanCsvReader, variational) {
590602
= stan::io::stan_csv_reader::parse(variational_stream, &out);
591603
variational_stream.close();
592604
ASSERT_EQ(1000, variational.metadata.num_samples);
605+
ASSERT_EQ(0, variational.adaptation.metric.size());
606+
}
607+
608+
TEST_F(StanIoStanCsvReader, read_nans) {
609+
std::ifstream datagen_stream;
610+
datagen_stream.open("src/test/unit/mcmc/test_csv_files/datagen_output.csv",
611+
std::ifstream::in);
612+
std::stringstream out;
613+
stan::io::stan_csv datagen
614+
= stan::io::stan_csv_reader::parse(datagen_stream, &out);
615+
datagen_stream.close();
616+
ASSERT_TRUE(std::isnan(datagen.samples(0, 2)));
593617
}

‎src/test/unit/io/test_csv_files/no_parameters_hmc.csv‎

Lines changed: 157 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,157 @@
1+
# stan_version_major = 2
2+
# stan_version_minor = 35
3+
# stan_version_patch = 0
4+
# model = sim_model
5+
# start_datetime = 2024-11-25 16:05:51 UTC
6+
# method = sample (Default)
7+
# sample
8+
# num_samples = 100
9+
# num_warmup = 100
10+
# save_warmup = false (Default)
11+
# thin = 1 (Default)
12+
# adapt
13+
# engaged = true (Default)
14+
# gamma = 0.05 (Default)
15+
# delta = 0.8 (Default)
16+
# kappa = 0.75 (Default)
17+
# t0 = 10 (Default)
18+
# init_buffer = 75 (Default)
19+
# term_buffer = 50 (Default)
20+
# window = 25 (Default)
21+
# save_metric = false (Default)
22+
# algorithm = hmc (Default)
23+
# hmc
24+
# engine = nuts (Default)
25+
# nuts
26+
# max_depth = 10 (Default)
27+
# metric = diag_e (Default)
28+
# metric_file = (Default)
29+
# stepsize = 1 (Default)
30+
# stepsize_jitter = 0 (Default)
31+
# num_chains = 1 (Default)
32+
# id = 1 (Default)
33+
# data
34+
# file = (Default)
35+
# init = 2 (Default)
36+
# random
37+
# seed = 1799096017 (Default)
38+
# output
39+
# file = output.csv (Default)
40+
# diagnostic_file = (Default)
41+
# refresh = 100 (Default)
42+
# sig_figs = -1 (Default)
43+
# profile_file = profile.csv (Default)
44+
# save_cmdstan_config = false (Default)
45+
# num_threads = 1 (Default)
46+
# stanc_version = %%NAME%%3 %%VERSION%%
47+
# stancflags =
48+
lp__,accept_stat__,stepsize__,treedepth__,n_leapfrog__,divergent__,energy__,x
49+
# Adaptation terminated
50+
# Step size = 1.35996e+14
51+
# Diagonal elements of inverse mass matrix:
52+
#
53+
0,1,1.35996e+14,1,1,0,0,8.57152
54+
0,1,1.35996e+14,1,1,0,0,-4.56101
55+
0,1,1.35996e+14,1,1,0,0,9.30019
56+
0,1,1.35996e+14,1,1,0,0,-14.6254
57+
0,1,1.35996e+14,1,1,0,0,8.30123
58+
0,1,1.35996e+14,1,1,0,0,19.2067
59+
0,1,1.35996e+14,1,1,0,0,12.6027
60+
0,1,1.35996e+14,1,1,0,0,-5.70173
61+
0,1,1.35996e+14,1,1,0,0,3.04002
62+
0,1,1.35996e+14,1,1,0,0,-8.28215
63+
0,1,1.35996e+14,1,1,0,0,3.20386
64+
0,1,1.35996e+14,1,1,0,0,-6.45707
65+
0,1,1.35996e+14,1,1,0,0,-5.29269
66+
0,1,1.35996e+14,1,1,0,0,5.1527
67+
0,1,1.35996e+14,1,1,0,0,-3.69775
68+
0,1,1.35996e+14,1,1,0,0,-1.77925
69+
0,1,1.35996e+14,1,1,0,0,-13.7562
70+
0,1,1.35996e+14,1,1,0,0,-0.789554
71+
0,1,1.35996e+14,1,1,0,0,-14.0911
72+
0,1,1.35996e+14,1,1,0,0,-4.71911
73+
0,1,1.35996e+14,1,1,0,0,-12.7874
74+
0,1,1.35996e+14,1,1,0,0,15.029
75+
0,1,1.35996e+14,1,1,0,0,-6.30709
76+
0,1,1.35996e+14,1,1,0,0,3.50805
77+
0,1,1.35996e+14,1,1,0,0,9.01495
78+
0,1,1.35996e+14,1,1,0,0,5.9809
79+
0,1,1.35996e+14,1,1,0,0,11.02
80+
0,1,1.35996e+14,1,1,0,0,9.68621
81+
0,1,1.35996e+14,1,1,0,0,-4.20659
82+
0,1,1.35996e+14,1,1,0,0,-7.71128
83+
0,1,1.35996e+14,1,1,0,0,-13.6125
84+
0,1,1.35996e+14,1,1,0,0,0.541567
85+
0,1,1.35996e+14,1,1,0,0,-6.96969
86+
0,1,1.35996e+14,1,1,0,0,2.49902
87+
0,1,1.35996e+14,1,1,0,0,8.31074
88+
0,1,1.35996e+14,1,1,0,0,13.4856
89+
0,1,1.35996e+14,1,1,0,0,-4.60591
90+
0,1,1.35996e+14,1,1,0,0,6.84733
91+
0,1,1.35996e+14,1,1,0,0,-16.8676
92+
0,1,1.35996e+14,1,1,0,0,4.43181
93+
0,1,1.35996e+14,1,1,0,0,-10.6785
94+
0,1,1.35996e+14,1,1,0,0,-5.56113
95+
0,1,1.35996e+14,1,1,0,0,-0.95401
96+
0,1,1.35996e+14,1,1,0,0,11.2198
97+
0,1,1.35996e+14,1,1,0,0,3.43417
98+
0,1,1.35996e+14,1,1,0,0,-11.2942
99+
0,1,1.35996e+14,1,1,0,0,-14.3029
100+
0,1,1.35996e+14,1,1,0,0,3.69492
101+
0,1,1.35996e+14,1,1,0,0,0.319324
102+
0,1,1.35996e+14,1,1,0,0,-5.95097
103+
0,1,1.35996e+14,1,1,0,0,5.99333
104+
0,1,1.35996e+14,1,1,0,0,-6.59629
105+
0,1,1.35996e+14,1,1,0,0,14.1795
106+
0,1,1.35996e+14,1,1,0,0,-7.58818
107+
0,1,1.35996e+14,1,1,0,0,4.89377
108+
0,1,1.35996e+14,1,1,0,0,4.63195
109+
0,1,1.35996e+14,1,1,0,0,-4.62905
110+
0,1,1.35996e+14,1,1,0,0,-11.4145
111+
0,1,1.35996e+14,1,1,0,0,4.03017
112+
0,1,1.35996e+14,1,1,0,0,-10.0459
113+
0,1,1.35996e+14,1,1,0,0,-11.8674
114+
0,1,1.35996e+14,1,1,0,0,-0.161997
115+
0,1,1.35996e+14,1,1,0,0,-5.75037
116+
0,1,1.35996e+14,1,1,0,0,-13.3027
117+
0,1,1.35996e+14,1,1,0,0,4.86817
118+
0,1,1.35996e+14,1,1,0,0,11.1937
119+
0,1,1.35996e+14,1,1,0,0,13.918
120+
0,1,1.35996e+14,1,1,0,0,-12.2423
121+
0,1,1.35996e+14,1,1,0,0,22.3588
122+
0,1,1.35996e+14,1,1,0,0,-8.30628
123+
0,1,1.35996e+14,1,1,0,0,-6.87127
124+
0,1,1.35996e+14,1,1,0,0,12.721
125+
0,1,1.35996e+14,1,1,0,0,-7.86135
126+
0,1,1.35996e+14,1,1,0,0,-8.56196
127+
0,1,1.35996e+14,1,1,0,0,-4.04709
128+
0,1,1.35996e+14,1,1,0,0,-21.6766
129+
0,1,1.35996e+14,1,1,0,0,-19.6485
130+
0,1,1.35996e+14,1,1,0,0,1.99421
131+
0,1,1.35996e+14,1,1,0,0,11.2645
132+
0,1,1.35996e+14,1,1,0,0,-9.35154
133+
0,1,1.35996e+14,1,1,0,0,-3.37081
134+
0,1,1.35996e+14,1,1,0,0,2.46874
135+
0,1,1.35996e+14,1,1,0,0,7.28248
136+
0,1,1.35996e+14,1,1,0,0,19.3846
137+
0,1,1.35996e+14,1,1,0,0,-4.92502
138+
0,1,1.35996e+14,1,1,0,0,-11.4687
139+
0,1,1.35996e+14,1,1,0,0,3.12818
140+
0,1,1.35996e+14,1,1,0,0,-5.79636
141+
0,1,1.35996e+14,1,1,0,0,12.2196
142+
0,1,1.35996e+14,1,1,0,0,5.59427
143+
0,1,1.35996e+14,1,1,0,0,-22.9084
144+
0,1,1.35996e+14,1,1,0,0,2.70951
145+
0,1,1.35996e+14,1,1,0,0,-0.604509
146+
0,1,1.35996e+14,1,1,0,0,1.06201
147+
0,1,1.35996e+14,1,1,0,0,-3.10986
148+
0,1,1.35996e+14,1,1,0,0,6.57804
149+
0,1,1.35996e+14,1,1,0,0,9.0426
150+
0,1,1.35996e+14,1,1,0,0,3.30854
151+
0,1,1.35996e+14,1,1,0,0,8.85378
152+
0,1,1.35996e+14,1,1,0,0,-1.71584
153+
#
154+
# Elapsed Time: 0 seconds (Warm-up)
155+
# 0 seconds (Sampling)
156+
# 0 seconds (Total)
157+
#

‎src/test/unit/mcmc/chainset_test.cpp‎

Lines changed: 30 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -170,7 +170,7 @@ TEST_F(McmcChains, summary_stats) {
170170
EXPECT_NEAR(theta_sd_expect, bern_chains.sd("theta"), 1e-5);
171171

172172
double theta_mad_expect = 0.12309;
173-
EXPECT_NEAR(theta_mad_expect, bern_chains.max_abs_deviation("theta"), 1e-5);
173+
EXPECT_NEAR(theta_mad_expect, bern_chains.med_abs_deviation("theta"), 1e-5);
174174

175175
double theta_mcse_mean_expect = 0.003234;
176176
EXPECT_NEAR(theta_mcse_mean_expect, bern_chains.mcse_mean("theta"), 1e-4);
@@ -207,3 +207,32 @@ TEST_F(McmcChains, summary_stats) {
207207
EXPECT_NEAR(theta_ac(i), theta_ac_expect(i), 0.0005);
208208
}
209209
}
210+
211+
TEST_F(McmcChains, quantile_tests) {
212+
std::ifstream datagen_stream;
213+
datagen_stream.open("src/test/unit/mcmc/test_csv_files/datagen_output.csv",
214+
std::ifstream::in);
215+
stan::io::stan_csv datagen_csv
216+
= stan::io::stan_csv_reader::parse(datagen_stream, &out);
217+
datagen_stream.close();
218+
stan::mcmc::chainset datagen_chains(datagen_csv);
219+
220+
Eigen::VectorXd probs(6);
221+
probs << 0.0, 0.01, 0.05, 0.95, 0.99, 1.0;
222+
Eigen::VectorXd stepsize_quantiles
223+
= datagen_chains.quantiles("stepsize__", probs);
224+
for (size_t i = 0; i < probs.size(); ++i) {
225+
EXPECT_TRUE(std::isnan(stepsize_quantiles(i)));
226+
}
227+
228+
double stepsize_MAD = datagen_chains.med_abs_deviation("stepsize__");
229+
EXPECT_TRUE(std::isnan(stepsize_MAD));
230+
231+
Eigen::VectorXd bad_probs(3);
232+
bad_probs << 5, 50, 95;
233+
Eigen::VectorXd y_sim_quantiles
234+
= datagen_chains.quantiles("y_sim[1]", bad_probs);
235+
for (size_t i = 0; i < bad_probs.size(); ++i) {
236+
EXPECT_TRUE(std::isnan(y_sim_quantiles(i)));
237+
}
238+
}

0 commit comments

Comments
 (0)