commit 8675beba45401991b00c05ac7fd9ba498575f7fb
parent 99643458cecabf51beaf15721ba7d104b121b9ac
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Mon, 4 Jul 2016 18:44:18 +0200
Add piecewise linear distribution
Diffstat:
5 files changed, 258 insertions(+), 0 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -156,6 +156,7 @@ if(NOT NO_TEST)
new_test(test_ssp_ran_triangle ${MATH_LIB})
new_test(test_ssp_rng_proxy)
new_test(test_ssp_ran_gaussian)
+ new_test(test_ssp_ran_piecewise_linear)
endif()
################################################################################
diff --git a/src/ssp.h b/src/ssp.h
@@ -59,6 +59,7 @@ struct ssp_rng;
struct ssp_rng_proxy;
struct ssp_ran_discrete;
struct ssp_ran_gaussian;
+struct ssp_ran_piecewise_linear;
/* Forward declaration of external types */
struct mem_allocator;
@@ -609,6 +610,34 @@ SSP_API double
ssp_distribution_gaussian_pdf
(const struct ssp_ran_gaussian* ran, double x);
+/*******************************************************************************
+ * Piecewise linear distribution
+ ******************************************************************************/
+
+SSP_API res_T
+ssp_ran_piecewise_linear_create
+ (struct mem_allocator* allocator, /* May be NULL <=> Use default allocator */
+ struct ssp_ran_piecewise_linear **ran);
+
+SSP_API res_T
+ssp_ran_piecewise_linear_init
+ (struct ssp_ran_piecewise_linear *ran,
+ const double* intervals,
+ const double* weights,
+ size_t size);
+
+SSP_API res_T
+ssp_ran_piecewise_linear_ref_get
+ (struct ssp_ran_piecewise_linear* ran);
+
+SSP_API res_T
+ssp_ran_piecewise_linear_ref_put
+ (struct ssp_ran_piecewise_linear* ran);
+
+SSP_API double
+ssp_ran_piecewise_linear_get
+ (const struct ssp_ran_piecewise_linear *ran, struct ssp_rng* rng);
+
END_DECLS
#endif /* SSP_H */
diff --git a/src/ssp_distributions.cpp b/src/ssp_distributions.cpp
@@ -33,6 +33,31 @@ gaussian_release(ref_T* ref)
}
/*******************************************************************************
+ * Piecewise linear distribution utilities
+ ******************************************************************************/
+
+using piecewise_dist=RAN_NAMESPACE::piecewise_linear_distribution<double>;
+
+struct ran_piecewise_linear_state
+{
+ piecewise_dist *distrib;
+};
+
+static void
+piecewise_release(ref_T* ref)
+{
+ ssp_ran_piecewise_linear* ran;
+ ASSERT(ref);
+ ran = CONTAINER_OF(ref, struct ssp_ran_piecewise_linear, ref);
+ if(ran->state) {
+ ran->state->distrib->~piecewise_dist();
+ MEM_RM(ran->allocator, ran->state->distrib);
+ MEM_RM(ran->allocator, ran->state);
+ }
+ MEM_RM(ran->allocator, ran);
+}
+
+/*******************************************************************************
* Gaussian distribution exported functions
******************************************************************************/
res_T
@@ -44,6 +69,9 @@ ssp_ran_gaussian_create
void *tmp = nullptr;
res_T res = RES_OK;
+ if (!allocator || !out_ran)
+ return RES_BAD_ARG;
+
allocator = allocator ? allocator : &mem_default_allocator;
ran = static_cast<struct ssp_ran_gaussian*>(
@@ -120,6 +148,7 @@ double
ssp_ran_gaussian_get
(const struct ssp_ran_gaussian *ran, struct ssp_rng* rng)
{
+ if(!ran || !rng) return RES_BAD_ARG;
class rng_cxx r(*rng);
return ran->state->distrib->operator()(r);
}
@@ -128,6 +157,99 @@ double
ssp_ran_gaussian_pdf
(const struct ssp_ran_gaussian *ran, double x)
{
+ if(!ran) return RES_BAD_ARG;
const double tmp = (x - ran->state->mu) * ran->state->K1;
return ran->state->K2 * exp(-0.5 * tmp * tmp);
}
+
+/*******************************************************************************
+ * Piecewise linear distribution exported functions
+ ******************************************************************************/
+
+res_T
+ssp_ran_piecewise_linear_create
+ (struct mem_allocator* allocator,
+ struct ssp_ran_piecewise_linear **out_ran)
+{
+ ssp_ran_piecewise_linear *ran = nullptr;
+ void *tmp = nullptr;
+ res_T res = RES_OK;
+
+ if (!allocator || !out_ran)
+ return RES_BAD_ARG;
+
+ allocator = allocator ? allocator : &mem_default_allocator;
+
+ ran = static_cast<ssp_ran_piecewise_linear*>(
+ MEM_CALLOC(allocator, 1, sizeof(ssp_ran_piecewise_linear)));
+ if (!ran) {
+ res = RES_MEM_ERR;
+ goto err;
+ }
+ ref_init(&ran->ref);
+ ran->state = static_cast<ran_piecewise_linear_state*>(
+ MEM_CALLOC(allocator, 1, sizeof(ran_piecewise_linear_state)));
+ if (!ran->state) {
+ res = RES_MEM_ERR;
+ goto err;
+ }
+ tmp = MEM_ALLOC(allocator, sizeof(piecewise_dist));
+ if (!tmp) {
+ res = RES_MEM_ERR;
+ goto err;
+ }
+ ran->allocator = allocator;
+ ran->state->distrib = static_cast<piecewise_dist*>(new(tmp) piecewise_dist);
+ if (out_ran) *out_ran = ran;
+ end:
+ return res;
+ err:
+ if (ran) {
+ ref_put(&ran->ref, piecewise_release);
+ ran = nullptr;
+ }
+ goto end;
+}
+
+res_T
+ssp_ran_piecewise_linear_init
+ (struct ssp_ran_piecewise_linear *ran,
+ const double* intervals,
+ const double* weights,
+ size_t size)
+{
+ if (!ran || !intervals || !weights || size < 2)
+ return RES_BAD_ARG;
+
+ piecewise_dist::param_type p{intervals, intervals + size, weights};
+ ran->state->distrib->param(p);
+
+ return RES_OK;
+}
+
+res_T
+ssp_ran_piecewise_linear_ref_get
+ (struct ssp_ran_piecewise_linear* ran)
+{
+ if(!ran) return RES_BAD_ARG;
+ ref_get(&ran->ref);
+ return RES_OK;
+}
+
+res_T
+ssp_ran_piecewise_linear_ref_put
+ (struct ssp_ran_piecewise_linear* ran)
+{
+ if(!ran) return RES_BAD_ARG;
+ ref_put(&ran->ref, piecewise_release);
+ return RES_OK;
+}
+
+double
+ssp_ran_piecewise_linear_get
+ (const struct ssp_ran_piecewise_linear *ran, struct ssp_rng* rng)
+{
+ if(!ran || !rng) return RES_BAD_ARG;
+ class rng_cxx r(*rng);
+ return ran->state->distrib->operator()(r);
+}
diff --git a/src/ssp_distributions.h b/src/ssp_distributions.h
@@ -6,6 +6,7 @@
/* Forward declaration of opaque types */
struct ran_gaussian_state;
+struct ran_piecewise_linear_state;
/* Forward declaration of external types */
struct ssp_rng;
@@ -19,4 +20,12 @@ struct ssp_ran_gaussian {
struct ran_gaussian_state *state;
};
+struct ssp_ran_piecewise_linear {
+ size_t sizeof_params;
+ size_t sizeof_state;
+ ref_T ref;
+ struct mem_allocator* allocator;
+ struct ran_piecewise_linear_state *state;
+};
+
#endif
diff --git a/src/test_ssp_ran_piecewise_linear.c b/src/test_ssp_ran_piecewise_linear.c
@@ -0,0 +1,97 @@
+/* Copyright (C) |Meso|Star> 2015-2016 (contact@meso-star.com)
+ *
+ * This software is a program whose purpose is to test the spp library.
+ *
+ * This software is governed by the CeCILL license under French law and
+ * abiding by the rules of distribution of free software. You can use,
+ * modify and/or redistribute the software under the terms of the CeCILL
+ * license as circulated by CEA, CNRS and INRIA at the following URL
+ * "http://www.cecill.info".
+ *
+ * As a counterpart to the access to the source code and rights to copy,
+ * modify and redistribute granted by the license, users are provided only
+ * with a limited warranty and the software's author, the holder of the
+ * economic rights, and the successive licensors have only limited
+ * liability.
+ *
+ * In this respect, the user's attention is drawn to the risks associated
+ * with loading, using, modifying and/or developing or reproducing the
+ * software by the user in light of its specific status of free software,
+ * that may mean that it is complicated to manipulate, and that also
+ * therefore means that it is reserved for developers and experienced
+ * professionals having in-depth computer knowledge. Users are therefore
+ * encouraged to load and test the software's suitability as regards their
+ * requirements in conditions enabling the security of their systems and/or
+ * data to be ensured and, more generally, to use and operate it in the
+ * same conditions as regards security.
+ *
+ * The fact that you are presently reading this means that you have had
+ * knowledge of the CeCILL license and that you accept its terms. */
+
+#include "ssp.h"
+#include "test_ssp_utils.h"
+
+#define NBS 1000000
+
+int
+main(int argc, char** argv)
+{
+ struct ssp_rng_proxy* proxy;
+ struct ssp_rng* rng;
+ struct mem_allocator allocator;
+ struct ssp_ran_piecewise_linear *pwl;
+ int i;
+ double exp_mean = 5, mean;
+ double exp_std = 10 / sqrt(12) /*sqrt((b - a)² / 12) */, std;
+ double x = 0, x2 = 0;
+ const double intervals[] = { 0, 1, 3, 5, 7, 8, 10 };
+ const double weights[] = { 1, 1, 1, 1, 1, 1, 1 };
+ (void)argc, (void)argv;
+
+ mem_init_proxy_allocator(&allocator, &mem_default_allocator);
+
+ CHECK(ssp_rng_proxy_create(&allocator, &ssp_rng_threefry, 1, &proxy), RES_OK);
+ CHECK(ssp_rng_proxy_create_rng(proxy, 0, &rng), RES_OK);
+
+ CHECK(ssp_ran_piecewise_linear_create(NULL, &pwl), RES_BAD_ARG);
+ CHECK(ssp_ran_piecewise_linear_create(&allocator, NULL), RES_BAD_ARG);
+ CHECK(ssp_ran_piecewise_linear_create(&allocator, &pwl), RES_OK);
+
+ CHECK(ssp_ran_piecewise_linear_init(NULL, intervals, weights, sizeof(intervals)/sizeof(double)), RES_BAD_ARG);
+ CHECK(ssp_ran_piecewise_linear_init(pwl, NULL, weights, sizeof(intervals)/sizeof(double)), RES_BAD_ARG);
+ CHECK(ssp_ran_piecewise_linear_init(pwl, intervals, NULL, sizeof(intervals)/sizeof(double)), RES_BAD_ARG);
+ CHECK(ssp_ran_piecewise_linear_init(pwl, intervals, weights, 1), RES_BAD_ARG);
+ CHECK(ssp_ran_piecewise_linear_init(pwl, intervals, weights, sizeof(intervals)/sizeof(double)), RES_OK);
+
+ CHECK(ssp_ran_piecewise_linear_ref_get(NULL), RES_BAD_ARG);
+ CHECK(ssp_ran_piecewise_linear_ref_get(pwl), RES_OK);
+
+ CHECK(ssp_ran_piecewise_linear_ref_put(NULL), RES_BAD_ARG);
+ CHECK(ssp_ran_piecewise_linear_ref_put(pwl), RES_OK);
+
+ CHECK(ssp_ran_piecewise_linear_get(NULL, rng), RES_BAD_ARG);
+ CHECK(ssp_ran_piecewise_linear_get(pwl, NULL), RES_BAD_ARG);
+
+ for (i = 0; i < NBS; i++) {
+ double _x = ssp_ran_piecewise_linear_get(pwl, rng);
+ CHECK(0 <= _x && _x <= 10, 1);
+ x += _x;
+ x2 += _x * _x;
+ }
+
+ mean = x/NBS;
+ std = sqrt(x2/NBS - x/NBS*x/NBS);
+ printf("%g %g\n", mean, std);
+ CHECK(fabs(mean - exp_mean) < 0.001, 1);
+ CHECK(fabs(std - exp_std) < 0.001, 1);
+
+ CHECK(ssp_ran_piecewise_linear_ref_put(pwl), RES_OK);
+
+ CHECK(ssp_rng_ref_put(rng), RES_OK);
+ CHECK(ssp_rng_proxy_ref_put(proxy), RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHECK(mem_allocated_size(), 0);
+ return 0;
+}
