commit b7e07af739449f04789ade3a82adf8acc83f73a4
parent d5a0eb6f6f9ee0b1a994e1839056d0cb57f9f9ea
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 4 Dec 2015 12:10:38 +0100
Add and test the discrete random variate
Tests are actually really basics.
Diffstat:
4 files changed, 326 insertions(+), 1 deletion(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -70,7 +70,7 @@ set(VERSION_MINOR 3)
set(VERSION_PATCH 0)
set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
-set(SSP_FILES_SRC ssp_rng.c ssp_rng_proxy.c)
+set(SSP_FILES_SRC ssp_ran_discrete.c ssp_rng.c ssp_rng_proxy.c)
set(SSP_FILES_INC ssp_rng_c.h)
set(SSP_FILES_INC_API ssp.h)
set(SSP_FILES_DOC COPYING.en COPYING.fr README.md)
@@ -139,6 +139,7 @@ if(NOT NO_TEST)
if(MSVC)
register_test(test_ssp_rng_aes test_ssp_rng aes)
endif()
+ new_test(test_ssp_ran_discrete)
new_test(test_ssp_ran_hemisphere ${MATH_LIB})
new_test(test_ssp_ran_sphere ${MATH_LIB})
new_test(test_ssp_ran_triangle ${MATH_LIB})
diff --git a/src/ssp.h b/src/ssp.h
@@ -57,6 +57,7 @@
/* Forward declaration of opaque types */
struct ssp_rng;
struct ssp_rng_proxy;
+struct ssp_ran_discrete;
/* Forward declaration of external types */
struct mem_allocator;
@@ -245,6 +246,39 @@ ssp_rng_proxy_get_type
struct ssp_rng_type* type);
/*******************************************************************************
+ * General discrete distribution
+ ******************************************************************************/
+/* Create a discrete random variate data structure from a list of weights.
+ * `weights' contain the weights of `nweights' discrete events. Its elements
+ * must be positive but they needn't add up to one. */
+SSP_API res_T
+ssp_ran_discrete_create
+ (struct mem_allocator* allocator, /* May be NULL <=> Use default allocator */
+ const double* weights,
+ const size_t nweights,
+ struct ssp_ran_discrete** ran);
+
+SSP_API res_T
+ssp_ran_discrete_ref_get
+ (struct ssp_ran_discrete* ran);
+
+SSP_API res_T
+ssp_ran_discrete_ref_put
+ (struct ssp_ran_discrete* ran);
+
+/* Return the index of the sampled discret event. */
+SSP_API size_t
+ssp_ran_discrete
+ (struct ssp_rng* rng,
+ struct ssp_ran_discrete* ran);
+
+/* Return the PDF of the discret event `i'. */
+SSP_API double
+ssp_ran_discrete_pdf
+ (const size_t i,
+ struct ssp_ran_discrete* ran);
+
+/*******************************************************************************
* Miscellaneous distributions
******************************************************************************/
/* Random variate from the exponential distribution with mean `mu':
diff --git a/src/ssp_ran_discrete.c b/src/ssp_ran_discrete.c
@@ -0,0 +1,182 @@
+/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com)
+ *
+ * This software is a library whose purpose is to generate [pseudo] random
+ * numbers and random variates.
+ *
+ * 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 <rsys/algorithm.h>
+#include <rsys/dynamic_array_double.h>
+#include <rsys/mem_allocator.h>
+#include <rsys/ref_count.h>
+
+struct ssp_ran_discrete {
+ struct darray_double cdf;
+ struct darray_double pdf;
+ ref_T ref;
+ struct mem_allocator* allocator;
+};
+
+/*******************************************************************************
+ * Helper function
+ ******************************************************************************/
+static FINLINE int
+cmp_double(const void* key, const void* val)
+{
+ const double* a = (const double*)key;
+ const double* b = (const double*)val;
+ ASSERT(a && b);
+ return *a < *b ? -1 : *a > *b ? 1 : 0;
+}
+
+static void
+ran_discrete_release(ref_T* ref)
+{
+ struct ssp_ran_discrete* ran;
+ ASSERT(ref);
+ ran = CONTAINER_OF(ref, struct ssp_ran_discrete, ref);
+ darray_double_release(&ran->cdf);
+ darray_double_release(&ran->pdf);
+ MEM_RM(ran->allocator, ran);
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+ssp_ran_discrete_create
+ (struct mem_allocator* mem_allocator,
+ const double* weights,
+ const size_t nweights,
+ struct ssp_ran_discrete** out_ran)
+{
+ struct mem_allocator* allocator = NULL;
+ struct ssp_ran_discrete* ran = NULL;
+ double* cdf;
+ double* pdf;
+ size_t i;
+ res_T res = RES_OK;
+
+ if(!weights || !nweights || !out_ran) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ allocator = mem_allocator ? mem_allocator : &mem_default_allocator;
+ ran = (struct ssp_ran_discrete*)
+ MEM_CALLOC(allocator, 1, sizeof(struct ssp_ran_discrete));
+ if(!ran) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ ref_init(&ran->ref);
+ ran->allocator = allocator;
+ darray_double_init(allocator, &ran->cdf);
+ darray_double_init(allocator, &ran->pdf);
+
+ res = darray_double_resize(&ran->cdf, nweights);
+ if(res != RES_OK) goto error;
+ res = darray_double_resize(&ran->pdf, nweights);
+ if(res != RES_OK) goto error;
+
+ cdf = darray_double_data_get(&ran->cdf);
+ pdf = darray_double_data_get(&ran->pdf);
+
+ if(weights[0] < 0.f) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ pdf[0] = weights[0];
+ cdf[0] = weights[0];
+ FOR_EACH(i, 1, nweights) {
+ if(weights[i] < 0.f) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ pdf[i] = weights[i];
+ cdf[i] = weights[i] + cdf[i-1];
+ }
+ if(!eq_eps(weights[nweights-1], 1.0, 1.e-8)) {
+ /* Normalize the weights */
+ const double len = cdf[nweights-1];
+ FOR_EACH(i, 0, nweights) {
+ pdf[i] /= len;
+ cdf[i] /= len;
+ }
+ }
+
+exit:
+ if(out_ran) *out_ran = ran;
+ return res;
+error:
+ if(ran) {
+ SSP(ran_discrete_ref_put(ran));
+ ran = NULL;
+ }
+ goto exit;
+}
+
+res_T
+ssp_ran_discrete_ref_get(struct ssp_ran_discrete* ran)
+{
+ if(!ran) return RES_BAD_ARG;
+ ref_get(&ran->ref);
+ return RES_OK;
+}
+
+res_T
+ssp_ran_discrete_ref_put(struct ssp_ran_discrete* ran)
+{
+ if(!ran) return RES_BAD_ARG;
+ ref_put(&ran->ref, ran_discrete_release);
+ return RES_OK;
+}
+
+size_t
+ssp_ran_discrete(struct ssp_rng* rng, struct ssp_ran_discrete* ran)
+{
+ const double* found;
+ double r;
+ ASSERT(rng && ran);
+
+ r = ssp_rng_canonical(rng);
+ found = (const double*)search_lower_bound(&r,
+ darray_double_cdata_get(&ran->cdf),
+ darray_double_size_get(&ran->cdf),
+ sizeof(double), cmp_double);
+ ASSERT(found && *found >= r);
+ return found - darray_double_cdata_get(&ran->cdf);
+}
+
+double
+ssp_ran_discrete_pdf(const size_t i, struct ssp_ran_discrete* ran)
+{
+ ASSERT(ran && i < darray_double_size_get(&ran->pdf));
+ return darray_double_cdata_get(&ran->pdf)[i];
+}
+
diff --git a/src/test_ssp_ran_discrete.c b/src/test_ssp_ran_discrete.c
@@ -0,0 +1,108 @@
+/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com)
+ *
+ * This software is a library whose purpose is to generate [pseudo] random
+ * numbers and random variates.
+ *
+ * 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"
+
+#include <string.h>
+
+#define NSAMPS 1024
+
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct ssp_ran_discrete* ran;
+ struct ssp_rng* rng;
+ const double weights[] = { 0.5, 0.1, 0.2, 0.05, 0.15 };
+ const size_t nweights = sizeof(weights)/sizeof(double);
+ size_t tmp[sizeof(weights)/sizeof(double)];
+ double accum;
+ size_t i;
+ (void)argc, (void)argv;
+
+ mem_init_proxy_allocator(&allocator, &mem_default_allocator);
+ CHECK(ssp_rng_create(&allocator, &ssp_rng_threefry, &rng), RES_OK);
+
+ CHECK(ssp_ran_discrete_create(NULL, NULL, 0, NULL), RES_BAD_ARG);
+ CHECK(ssp_ran_discrete_create(NULL, weights, 0, NULL), RES_BAD_ARG);
+ CHECK(ssp_ran_discrete_create(NULL, NULL, nweights, NULL), RES_BAD_ARG);
+ CHECK(ssp_ran_discrete_create(NULL, weights, nweights, NULL), RES_BAD_ARG);
+ CHECK(ssp_ran_discrete_create(NULL, NULL, 0, &ran), RES_BAD_ARG);
+ CHECK(ssp_ran_discrete_create(NULL, weights, 0, &ran), RES_BAD_ARG);
+ CHECK(ssp_ran_discrete_create(NULL, NULL, nweights, &ran), RES_BAD_ARG);
+ CHECK(ssp_ran_discrete_create(NULL, weights, nweights, &ran), RES_OK);
+
+ CHECK(ssp_ran_discrete_ref_get(NULL), RES_BAD_ARG);
+ CHECK(ssp_ran_discrete_ref_get(ran), RES_OK);
+ CHECK(ssp_ran_discrete_ref_put(NULL), RES_BAD_ARG);
+ CHECK(ssp_ran_discrete_ref_put(ran), RES_OK);
+ CHECK(ssp_ran_discrete_ref_put(ran), RES_OK);
+
+ CHECK(ssp_ran_discrete_create(&allocator, weights, nweights, &ran), RES_OK);
+
+ memset(tmp, 0, sizeof(tmp));
+ FOR_EACH(i, 0, NSAMPS) {
+ const size_t k = ssp_ran_discrete(rng, ran);
+ double pdf;
+ CHECK(k < nweights, 1);
+ pdf = ssp_ran_discrete_pdf(k, ran);
+ ++tmp[k];
+ CHECK(pdf, weights[k]);
+ CHECK(pdf >= 0.f, 1);
+ CHECK(pdf <= 1.f, 1);
+ }
+ FOR_EACH(i, 0, nweights) {
+ NCHECK(tmp[i], 0);
+ }
+ CHECK(ssp_ran_discrete_ref_put(ran), RES_OK);
+
+ CHECK(ssp_ran_discrete_create(&allocator, weights, nweights-1, &ran), RES_OK);
+ FOR_EACH(i, 0, NSAMPS) {
+ const size_t k = ssp_ran_discrete(rng, ran);
+ double pdf;
+ CHECK(k < nweights-1, 1);
+ pdf = ssp_ran_discrete_pdf(k, ran);
+ CHECK(pdf >= 0.f, 1);
+ CHECK(pdf <= 1.f, 1);
+ }
+ accum = 0;
+ FOR_EACH(i, 0, nweights-1) accum += ssp_ran_discrete_pdf(i, ran);
+ CHECK(eq_eps(accum, 1, 1.e-8), 1);
+ CHECK(ssp_ran_discrete_ref_put(ran), RES_OK);
+
+ CHECK(ssp_rng_ref_put(rng), RES_OK);
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHECK(mem_allocated_size(), 0);
+ return 0;
+}
+
