Skip to content

Commit

Permalink
changed test
Browse files Browse the repository at this point in the history
  • Loading branch information
@Franzi2114
Franzi2114 committed Aug 1, 2023
1 parent 0d188e9 commit 24759b4
Show file tree
Hide file tree
Showing 3 changed files with 43 additions and 47 deletions.
25 changes: 11 additions & 14 deletions stan/math/prim/functor/hcubature.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -192,10 +192,9 @@ template <typename F, typename ParsPairT>
std::pair<double, double> gauss_kronrod(const F& integrand, const double a,
const double b,
const ParsPairT& pars_pair) {
std::vector<double> c(1, 0);
std::vector<double> cp(1, 0);
std::vector<double> cm(1, 0);
c[0] = 0.5 * (a + b);
Eigen::VectorXd c {{0.5 * (a + b)}};
Eigen::VectorXd cp(1);
Eigen::VectorXd cm(1);
double delta = 0.5 * (b - a);
double f0 = math::apply(
[&integrand, &c](auto&&... args) { return integrand(c, args...); },
Expand Down Expand Up @@ -308,10 +307,10 @@ std::tuple<double, double, int> integrate_GenzMalik(
double twelvef1 = 12 * f1;

double maxdivdiff = 0.0;
std::vector<double> divdiff(dim);
std::vector<double> p2(dim);
std::vector<double> p3(dim);
std::vector<double> cc(dim, 0);
Eigen::VectorXd divdiff(dim);
Eigen::VectorXd p2(dim);
Eigen::VectorXd p3(dim);
Eigen::VectorXd cc(dim);

for (std::size_t i = 0; i != dim; i++) {
for (std::size_t j = 0; j != dim; j++) {
Expand Down Expand Up @@ -352,7 +351,7 @@ std::tuple<double, double, int> integrate_GenzMalik(
f3 += f3i;
divdiff[i] = fabs(f3i + twelvef1 - 7 * f2i);
}
std::vector<double> p4(dim);
Eigen::VectorXd p4(dim);
double f4 = 0.0;
for (std::size_t i = 0; i != points[2].cols(); i++) {
for (std::size_t j = 0; j != dim; j++) {
Expand All @@ -367,7 +366,7 @@ std::tuple<double, double, int> integrate_GenzMalik(
f4 += temp;
}
double f5 = 0.0;
std::vector<double> p5(dim);
Eigen::VectorXd p5(dim);
for (std::size_t i = 0; i != points[3].cols(); i++) {
for (std::size_t j = 0; j != dim; j++) {
p5[j] = deltac[j] * points[3](j, i);
Expand Down Expand Up @@ -521,12 +520,10 @@ double hcubature(const F& integrand, const ParsTuple& pars, const int dim,
auto&& box = ms[err_idx];

double w = (box.b_[box.kdiv_] - box.a_[box.kdiv_]) / 2;
Eigen::VectorXd ma
= Eigen::Map<const Eigen::VectorXd>(box.a_.data(), box.a_.size());
Eigen::VectorXd ma = Eigen::Map<const Eigen::VectorXd>(box.a_.data(), box.a_.size());

ma[box.kdiv_] += w;
Eigen::VectorXd mb
= Eigen::Map<const Eigen::VectorXd>(box.b_.data(), box.b_.size());
Eigen::VectorXd mb = Eigen::Map<const Eigen::VectorXd>(box.b_.data(), box.b_.size());
mb[box.kdiv_] -= w;

double result_1;
Expand Down
3 changes: 2 additions & 1 deletion stan/math/prim/prob/wiener_full_lpdf.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -288,7 +288,7 @@ inline ReturnT wiener_full_lpdf(const T_y& y, const T_a& a, const T_t0& t0,
T_st0>::value) {
return 0;
}

using T_y_ref = ref_type_t<T_y>;
using T_a_ref = ref_type_t<T_a>;
using T_v_ref = ref_type_t<T_v>;
Expand Down Expand Up @@ -392,6 +392,7 @@ inline ReturnT wiener_full_lpdf(const T_y& y, const T_a& a, const T_t0& t0,
}
}


const double log_error_density = log(1e-6); // precision for density
const double error_bound = precision_derivatives; // precision for
// derivatives (controllable by user)
Expand Down
62 changes: 30 additions & 32 deletions test/unit/math/prim/functor/hcubature_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -105,8 +105,8 @@ stan::return_type_t<T_x, double> f7(const T_x& x, double a) {

template <typename F, typename ArgsTupleT>
void test_integration(const F& f, const ArgsTupleT& pars, int dim,
std::vector<double> a, std::vector<double> b, int maxEval,
double reqAbsError, std::vector<double> reqRelError,
const Eigen::VectorXd& a, const Eigen::VectorXd& b, int maxEval,
double reqAbsError, const Eigen::VectorXd& reqRelError,
double val) {
using stan::math::hcubature;

Expand All @@ -122,46 +122,44 @@ TEST(StanMath_hcubature_prim, test1) {
// https://www.quantargo.com/help/r/latest/packages/cubature/2.0.4.1/hcubature

int dim = 1;
std::vector<double> a = {0.0};
std::vector<double> b = {1.0};
std::vector<double> reqRelError = {1e-4, 1e-6, 1e-7};
test_integration(hcubature_test::f1<std::vector<double>>, std::make_tuple(),
const Eigen::VectorXd a {{0.0}};
const Eigen::VectorXd b {{1.0}};
const Eigen::VectorXd reqRelError {{1e-4, 1e-6, 1e-7}};
test_integration(hcubature_test::f1<const Eigen::VectorXd>, std::make_tuple(),
dim, a, b, 6000, 0.0, reqRelError, 0.841471);

dim = 2;
a = {0.0, 0.0};
b = {1.0, 1.0};
reqRelError = {1e-4, 1e-6, 1e-7};
test_integration(hcubature_test::f2<std::vector<double>>, std::make_tuple(),
dim, a, b, 6000, 0.0, reqRelError, 0.7080734);
const Eigen::VectorXd a_2 {{0.0, 0.0}};
const Eigen::VectorXd b_2 {{1.0, 1.0}};
test_integration(hcubature_test::f2<const Eigen::VectorXd>, std::make_tuple(),
dim, a_2, b_2, 6000, 0.0, reqRelError, 0.7080734);

reqRelError = {1e-4};
test_integration(hcubature_test::f3<std::vector<double>>,
std::make_tuple(0.50124145262344534123412), dim, a, b, 10000,
0.0, reqRelError, 0.1972807);
const Eigen::VectorXd reqRelError_2 {{1e-4}};
test_integration(hcubature_test::f3<const Eigen::VectorXd>,
std::make_tuple(0.50124145262344534123412), dim, a_2, b_2, 10000,
0.0, reqRelError_2, 0.1972807);

// (Gaussian centered at 1/2)
reqRelError = {1e-4, 1e-6, 1e-7};
test_integration(hcubature_test::f4<std::vector<double>>,
std::make_tuple(0.1), dim, a, b, 6000, 0.0, reqRelError, 1);
test_integration(hcubature_test::f4<const Eigen::VectorXd>,
std::make_tuple(0.1), dim, a_2, b_2, 6000, 0.0, reqRelError, 1);

dim = 3;
a = {0.0, 0.0, 0.0};
b = {1.0, 1.0, 1.0};
reqRelError = {1e-4, 1e-6};
test_integration(hcubature_test::f5<std::vector<double>>, std::make_tuple(),
dim, a, b, 6000, 0.0, reqRelError, 1.00001);
const Eigen::VectorXd a_3 {{0.0, 0.0, 0.0}};
const Eigen::VectorXd b_3 {{1.0, 1.0, 1.0}};
const Eigen::VectorXd reqRelError_3 {{1e-4, 1e-6}};
test_integration(hcubature_test::f5<const Eigen::VectorXd>, std::make_tuple(),
dim, a_3, b_3, 6000, 0.0, reqRelError_3, 1.00001);

reqRelError = {1e-4, 1e-6, 1e-8};
test_integration(hcubature_test::f6<std::vector<double>>, std::make_tuple(),
dim, a, b, 6000, 0.0, reqRelError, 1);
const Eigen::VectorXd reqRelError_4 {{1e-4, 1e-6, 1e-8}};
test_integration(hcubature_test::f6<const Eigen::VectorXd>, std::make_tuple(),
dim, a_3, b_3, 6000, 0.0, reqRelError_4, 1);

// (Tsuda's example)
dim = 4;
a = {0.0, 0.0, 0.0, 0.0};
b = {1.0, 1.0, 1.0, 1.0};
reqRelError = {1e-4, 1e-6};
test_integration(hcubature_test::f7<std::vector<double>>,
std::make_tuple((1 + sqrt(10.0)) / 9.0), dim, a, b, 20000,
0.0, reqRelError, 0.999998);
const Eigen::VectorXd a_4 {{0.0, 0.0, 0.0, 0.0}};
const Eigen::VectorXd b_4 {{1.0, 1.0, 1.0, 1.0}};
test_integration(hcubature_test::f7<const Eigen::VectorXd>,
std::make_tuple((1 + sqrt(10.0)) / 9.0), dim, a_4, b_4, 20000,
0.0, reqRelError_3, 0.999998);

}

0 comments on commit 24759b4

Please sign in to comment.