commit 5b4bfe8c675f554652e9becc7a7c88a05262b772
parent 4932c7c325e283cb6c25eacc3c2fa1a7f5fc4e4a
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Wed, 25 Oct 2017 15:17:28 +0200
Add some missing pdf functions.
Diffstat:
3 files changed, 113 insertions(+), 2 deletions(-)
diff --git a/src/ssp.h b/src/ssp.h
@@ -815,6 +815,16 @@ ssp_ranst_piecewise_linear_float_get
(const struct ssp_ranst_piecewise_linear_float *ran,
struct ssp_rng* rng);
+SSP_API double
+ssp_ranst_piecewise_linear_pdf
+ (const struct ssp_ranst_piecewise_linear *ran,
+ double x);
+
+SSP_API float
+ssp_ranst_piecewise_linear_float_pdf
+ (const struct ssp_ranst_piecewise_linear_float *ran,
+ float x);
+
/*******************************************************************************
* Uniform disk distribution.
*
@@ -833,6 +843,18 @@ ssp_ran_uniform_disk_float
const float radius,
float pt[2]);
+static INLINE double
+ssp_ran_uniform_disk_pdf(const double radius)
+{
+ return 1 / (radius * radius);
+}
+
+static INLINE float
+ssp_ran_uniform_disk_float_pdf(const float radius)
+{
+ return 1 / (radius * radius);
+}
+
END_DECLS
#endif /* SSP_H */
diff --git a/src/ssp_ranst_piecewise_linear.c b/src/ssp_ranst_piecewise_linear.c
@@ -34,6 +34,8 @@
#include <rsys/mem_allocator.h>
#include <rsys/ref_count.h>
+#include <algorithm>
+
using piecewise_dist = RAN_NAMESPACE::piecewise_linear_distribution<double>;
using piecewise_dist_float = RAN_NAMESPACE::piecewise_linear_distribution<float>;
@@ -176,9 +178,16 @@ ssp_ranst_piecewise_linear_setup
const double* weights,
const size_t size)
{
+ size_t i;
if(!ran || !intervals || !weights || size < 2)
return RES_BAD_ARG;
+ /* Checking param validity to avoid an assert when using ran */
+ for(i=0; i < size-1; i++) {
+ if(weights[i] < 0) return RES_BAD_ARG;
+ if(intervals[i+1] <= intervals[i]) return RES_BAD_ARG;
+ }
+ if (intervals[size-1] - intervals[size-2] <= 0) return RES_BAD_ARG;
piecewise_dist::param_type p{intervals, intervals + size, weights};
ran->distrib->param(p);
return RES_OK;
@@ -191,9 +200,16 @@ ssp_ranst_piecewise_linear_float_setup
const float* weights,
const size_t size)
{
+ size_t i;
if(!ran || !intervals || !weights || size < 2)
return RES_BAD_ARG;
+ /* Checking param validity to avoid an assert when using ran */
+ for(i=0; i < size-1; i++) {
+ if(weights[i] < 0) return RES_BAD_ARG;
+ if(intervals[i+1] <= intervals[i]) return RES_BAD_ARG;
+ }
+ if (weights[size-1] < 0) return RES_BAD_ARG;
piecewise_dist_float::param_type p{ intervals, intervals + size, weights };
ran->distrib->param(p);
return RES_OK;
@@ -253,3 +269,42 @@ ssp_ranst_piecewise_linear_float_get
return ran->distrib->operator()(r);
}
+double
+ssp_ranst_piecewise_linear_pdf
+ (const struct ssp_ranst_piecewise_linear *ran,
+ double x)
+{
+ ASSERT(ran);
+ if (x<ran->distrib->min() || x>ran->distrib->max())
+ return 0;
+
+ const auto& inter = ran->distrib->intervals();
+ const auto& dens = ran->distrib->densities();
+ auto b = std::lower_bound(inter.begin(), inter.end(), x);
+ size_t idx = b - inter.begin();
+ if (x == *b) return dens[idx];
+ idx--;
+ ASSERT(idx < inter.size() - 1);
+ return (dens[idx+1] * (x - inter[idx]) + dens[idx] * (inter[idx+1] - x))
+ / (inter[idx+1]- inter[idx]);
+}
+
+float
+ssp_ranst_piecewise_linear_float_pdf
+ (const struct ssp_ranst_piecewise_linear_float *ran,
+ float x)
+{
+ ASSERT(ran);
+ if (x<ran->distrib->min() || x>ran->distrib->max())
+ return 0;
+
+ const auto& inter = ran->distrib->intervals();
+ const auto& dens = ran->distrib->densities();
+ auto b = std::lower_bound(inter.begin(), inter.end(), x);
+ size_t idx = b - inter.begin();
+ if (x == *b) return dens[idx];
+ idx--;
+ ASSERT(idx < inter.size() - 1);
+ return (dens[idx+1] * (x - inter[idx]) + dens[idx] * (inter[idx+1] - x))
+ / (inter[idx+1] - inter[idx]);
+}
+\ No newline at end of file
diff --git a/src/test_ssp_ran_piecewise_linear.h b/src/test_ssp_ran_piecewise_linear.h
@@ -45,9 +45,11 @@
#define RANST_PIECEWISE_LINEAR_CREATE ssp_ranst_piecewise_linear_create
#define RANST_PIECEWISE_LINEAR_SETUP ssp_ranst_piecewise_linear_setup
#define RANST_PIECEWISE_LINEAR_GET ssp_ranst_piecewise_linear_get
+#define RANST_PIECEWISE_LINEAR_PDF ssp_ranst_piecewise_linear_pdf
#define RANST_PIECEWISE_LINEAR_REF_GET ssp_ranst_piecewise_linear_ref_get
#define RANST_PIECEWISE_LINEAR_REF_PUT ssp_ranst_piecewise_linear_ref_put
#define EQ_EPS_R eq_eps
+#define EPS_R DBL_EPSILON
#define SQRT sqrt
#elif TYPE_FLOAT==1
@@ -57,9 +59,11 @@
#define RANST_PIECEWISE_LINEAR_CREATE ssp_ranst_piecewise_linear_float_create
#define RANST_PIECEWISE_LINEAR_SETUP ssp_ranst_piecewise_linear_float_setup
#define RANST_PIECEWISE_LINEAR_GET ssp_ranst_piecewise_linear_float_get
+#define RANST_PIECEWISE_LINEAR_PDF ssp_ranst_piecewise_linear_float_pdf
#define RANST_PIECEWISE_LINEAR_REF_GET ssp_ranst_piecewise_linear_float_ref_get
#define RANST_PIECEWISE_LINEAR_REF_PUT ssp_ranst_piecewise_linear_float_ref_put
#define EQ_EPS_R eq_epsf
+#define EPS_R FLT_EPSILON
#define SQRT sqrtf
#else
#error "TYPE_SUFFIX must be defined either 0 or 1"
@@ -75,8 +79,8 @@ TEST()
REAL exp_mean = 5, mean;
REAL exp_std = 10 / SQRT(12) /*sqrt((b - a)² / 12) */, std;
REAL x = 0, x2 = 0;
- const REAL intervals[] = { 0, 1, 3, 5, 7, 8, 10 };
- const REAL weights[] = { 1, 1, 1, 1, 1, 1, 1 };
+ REAL intervals[] = { 0, 1, 3, 5, 7, 8, 10 };
+ REAL weights[] = { 1, 1, 1, 1, 1, 1, 1 };
mem_init_proxy_allocator(&allocator, &mem_default_allocator);
@@ -100,6 +104,27 @@ TEST()
CHECK(RANST_PIECEWISE_LINEAR_SETUP
(pwl, intervals, weights, sizeof(intervals)/sizeof(REAL)), RES_OK);
+ CHECK(RANST_PIECEWISE_LINEAR_REF_PUT(pwl), RES_OK);
+ CHECK(RANST_PIECEWISE_LINEAR_CREATE(&allocator, &pwl), RES_OK);
+
+ weights[1] = -1;
+ CHECK(RANST_PIECEWISE_LINEAR_SETUP
+ (pwl, intervals, weights, sizeof(intervals) / sizeof(REAL)), RES_BAD_ARG);
+ weights[1] = 1;
+
+ intervals[1] = 4;
+ CHECK(RANST_PIECEWISE_LINEAR_SETUP
+ (pwl, intervals, weights, sizeof(intervals) / sizeof(REAL)), RES_BAD_ARG);
+ intervals[1] = 1;
+
+ intervals[1] = 3;
+ CHECK(RANST_PIECEWISE_LINEAR_SETUP
+ (pwl, intervals, weights, sizeof(intervals) / sizeof(REAL)), RES_BAD_ARG);
+ intervals[1] = 1;
+
+ CHECK(RANST_PIECEWISE_LINEAR_SETUP
+ (pwl, intervals, weights, sizeof(intervals) / sizeof(REAL)), RES_OK);
+
CHECK(RANST_PIECEWISE_LINEAR_REF_GET(NULL), RES_BAD_ARG);
CHECK(RANST_PIECEWISE_LINEAR_REF_GET(pwl), RES_OK);
@@ -109,9 +134,15 @@ TEST()
FOR_EACH(i, 0, NBS) {
REAL r = RANST_PIECEWISE_LINEAR_GET(pwl, rng);
CHECK(0 <= r && r <= 10, 1);
+ REAL pdf = RANST_PIECEWISE_LINEAR_PDF(pwl, r);
+ CHECK(EQ_EPS_R(pdf, (REAL)0.1, EPS_R), 1);
x += r;
x2 += r * r;
}
+ CHECK(EQ_EPS_R(RANST_PIECEWISE_LINEAR_PDF(pwl, 0), (REAL)0.1, EPS_R), 1);
+ CHECK(EQ_EPS_R(RANST_PIECEWISE_LINEAR_PDF(pwl, 10), (REAL)0.1, EPS_R), 1);
+ CHECK(RANST_PIECEWISE_LINEAR_PDF(pwl, -1), 0);
+ CHECK(RANST_PIECEWISE_LINEAR_PDF(pwl, 11), 0);
mean = x/NBS;
std = SQRT(x2/NBS - mean*mean);
@@ -134,8 +165,10 @@ TEST()
#undef RANST_PIECEWISE_LINEAR_CREATE
#undef RANST_PIECEWISE_LINEAR_SETUP
#undef RANST_PIECEWISE_LINEAR_GET
+#undef RANST_PIECEWISE_LINEAR_PDF
#undef RANST_PIECEWISE_LINEAR_REF_GET
#undef RANST_PIECEWISE_LINEAR_REF_PUT
#undef EQ_EPS_R
+#undef EPS_R
#undef SQRT
#undef TYPE_FLOAT
