commit 47102f3c8e342e35326a71b6bec2030008fdf995
parent 7988e07b4898f422c200216a9870942fcf4e9572
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Thu, 7 Jul 2016 16:55:08 +0200
API break: change naming pattern for state-based distributions.
discrete distribution now uses C++11 std::discrete_distribution, though.
Diffstat:
7 files changed, 326 insertions(+), 397 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -71,7 +71,6 @@ set(VERSION_PATCH 0)
set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
set(SSP_FILES_SRC
- ssp_ran_discrete.c
ssp_rng.c
ssp_rng_proxy.c
ssp_distributions.cpp)
diff --git a/src/ssp.h b/src/ssp.h
@@ -57,8 +57,9 @@
/* Forward declaration of opaque types */
struct ssp_rng;
struct ssp_rng_proxy;
-struct ssp_ran_discrete;
-struct ssp_ran_type;
+struct ssp_ranst_discrete;
+struct ssp_ranst_gaussian;
+struct ssp_ranst_piecewise_linear;
/* Forward declaration of external types */
struct mem_allocator;
@@ -252,41 +253,41 @@ ssp_rng_proxy_get_type
struct ssp_rng_type* type);
/*******************************************************************************
- * General discrete distribution
+ * General discrete distribution with state
******************************************************************************/
/* 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
+ssp_ranst_discrete_create
(struct mem_allocator* allocator, /* May be NULL <=> Use default allocator */
- struct ssp_ran_discrete** ran);
+ struct ssp_ranst_discrete** ran);
SSP_API res_T
-ssp_ran_discrete_setup
- (struct ssp_ran_discrete* discrete,
+ssp_ranst_discrete_setup
+ (struct ssp_ranst_discrete* discrete,
const double* weights,
const size_t nweights);
SSP_API res_T
-ssp_ran_discrete_ref_get
- (struct ssp_ran_discrete* ran);
+ssp_ranst_discrete_ref_get
+ (struct ssp_ranst_discrete* ran);
SSP_API res_T
-ssp_ran_discrete_ref_put
- (struct ssp_ran_discrete* ran);
+ssp_ranst_discrete_ref_put
+ (struct ssp_ranst_discrete* ran);
/* Return the index of the sampled discret event. */
SSP_API size_t
-ssp_ran_discrete
+ssp_ranst_discrete_get
(struct ssp_rng* rng,
- struct ssp_ran_discrete* ran);
+ const struct ssp_ranst_discrete* ran);
/* Return the PDF of the discret event `i'. */
SSP_API double
-ssp_ran_discrete_pdf
+ssp_ranst_discrete_pdf
(const size_t i,
- struct ssp_ran_discrete* ran);
+ const struct ssp_ranst_discrete* ran);
/*******************************************************************************
* Miscellaneous distributions
@@ -303,7 +304,8 @@ ssp_ran_exp_pdf
(const double x,
const double mu);
-/* Random variate from the gaussian (or normal) distribution with mean `mu':
+/* Stateless random variate from the gaussian (or normal) distribution
+ * with mean `mu':
* P(x) dx = 1 / (sigma*sqrt(2*PI)) * exp(1/2*((x-mu)/sigma)^2) dx */
SSP_API double
ssp_ran_gaussian
@@ -579,65 +581,65 @@ ssp_ran_sphere_hg
}
/*******************************************************************************
- * Gaussian distribution
+ * Gaussian distribution with state
******************************************************************************/
SSP_API res_T
-ssp_ran_gaussian_create
+ssp_ranst_gaussian_create
(struct mem_allocator* allocator, /* May be NULL <=> Use default allocator */
- struct ssp_ran_type** ran);
+ struct ssp_ranst_gaussian** ran);
SSP_API res_T
-ssp_ran_gaussian_init
- (struct ssp_ran_type* ran,
+ssp_ranst_gaussian_setup
+ (struct ssp_ranst_gaussian* ran,
double mu,
double sigma);
SSP_API res_T
-ssp_ran_gaussian_ref_get
- (struct ssp_ran_type* ran);
+ssp_ranst_gaussian_ref_get
+ (struct ssp_ranst_gaussian* ran);
SSP_API res_T
-ssp_ran_gaussian_ref_put
- (struct ssp_ran_type* ran);
+ssp_ranst_gaussian_ref_put
+ (struct ssp_ranst_gaussian* ran);
SSP_API double
-ssp_ran_gaussian_get
- (const struct ssp_ran_type* ran,
+ssp_ranst_gaussian_get
+ (const struct ssp_ranst_gaussian* ran,
struct ssp_rng* rng);
SSP_API double
-ssp_distribution_gaussian_pdf
- (const struct ssp_ran_type* ran,
+ssp_ranst_gaussian_pdf
+ (const struct ssp_ranst_gaussian* ran,
double x);
/*******************************************************************************
- * Piecewise linear distribution
+ * Piecewise linear distribution with state
******************************************************************************/
SSP_API res_T
-ssp_ran_piecewise_linear_create
+ssp_ranst_piecewise_linear_create
(struct mem_allocator* allocator, /* May be NULL <=> Use default allocator */
- struct ssp_ran_type **ran);
+ struct ssp_ranst_piecewise_linear **ran);
SSP_API res_T
-ssp_ran_piecewise_linear_init
- (struct ssp_ran_type *ran,
+ssp_ranst_piecewise_linear_setup
+ (struct ssp_ranst_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_type* ran);
+ssp_ranst_piecewise_linear_ref_get
+ (struct ssp_ranst_piecewise_linear* ran);
SSP_API res_T
-ssp_ran_piecewise_linear_ref_put
- (struct ssp_ran_type* ran);
+ssp_ranst_piecewise_linear_ref_put
+ (struct ssp_ranst_piecewise_linear* ran);
SSP_API double
-ssp_ran_piecewise_linear_get
- (const struct ssp_ran_type *ran,
+ssp_ranst_piecewise_linear_get
+ (const struct ssp_ranst_piecewise_linear *ran,
struct ssp_rng* rng);
/*******************************************************************************
diff --git a/src/ssp_distributions.cpp b/src/ssp_distributions.cpp
@@ -1,15 +1,10 @@
#include <rsys/mem_allocator.h>
#include <rsys/ref_count.h>
+#include <rsys/dynamic_array_double.h>
#include <rsys/math.h>
+
#include "ssp_rng_c.h"
-/* one single type for all distributions
- * only the state type depends on the distribution type */
-struct ssp_ran_type {
- ref_T ref;
- struct mem_allocator* allocator;
- void* state;
-};
/*******************************************************************************
* Gaussian distribution utilities
@@ -17,8 +12,10 @@ struct ssp_ran_type {
using normal_dist=RAN_NAMESPACE::normal_distribution<double>;
-struct ran_gaussian_state
+struct ssp_ranst_gaussian
{
+ ref_T ref;
+ struct mem_allocator* allocator;
normal_dist* distrib;
double mu;
double K1; // 1.0 / sigma;
@@ -28,15 +25,12 @@ struct ran_gaussian_state
static void
gaussian_release(ref_T* ref)
{
- ssp_ran_type* ran;
- ran_gaussian_state* state;
+ ssp_ranst_gaussian* ran;
ASSERT(ref);
- ran = CONTAINER_OF(ref, ssp_ran_type, ref);
- state = static_cast<ran_gaussian_state*>(ran->state);
- if(state) {
- state->distrib->~normal_dist();
- MEM_RM(ran->allocator, state->distrib);
- MEM_RM(ran->allocator, state);
+ ran = CONTAINER_OF(ref, ssp_ranst_gaussian, ref);
+ if(ran->distrib) {
+ ran->distrib->~normal_dist();
+ MEM_RM(ran->allocator, ran->distrib);
}
MEM_RM(ran->allocator, ran);
}
@@ -47,38 +41,194 @@ gaussian_release(ref_T* ref)
using piecewise_dist=RAN_NAMESPACE::piecewise_linear_distribution<double>;
-struct ran_piecewise_linear_state
+struct ssp_ranst_piecewise_linear
{
+ ref_T ref;
+ struct mem_allocator* allocator;
piecewise_dist* distrib;
};
static void
piecewise_release(ref_T* ref)
{
- ssp_ran_type* ran;
- ran_piecewise_linear_state* state;
+ ssp_ranst_piecewise_linear* ran;
+ ASSERT(ref);
+ ran = CONTAINER_OF(ref, ssp_ranst_piecewise_linear, ref);
+ if(ran->distrib) {
+ ran->distrib->~piecewise_dist();
+ MEM_RM(ran->allocator, ran->distrib);
+ }
+ MEM_RM(ran->allocator, ran);
+}
+
+/*******************************************************************************
+ * Discrete distribution utilities
+ ******************************************************************************/
+
+using discrete_dist=RAN_NAMESPACE::discrete_distribution<size_t>;
+
+struct ssp_ranst_discrete {
+ struct darray_double pdf;
+ ref_T ref;
+ struct mem_allocator* allocator;
+ discrete_dist* distrib;
+};
+
+static void
+ran_discrete_release(ref_T* ref)
+{
+ struct ssp_ranst_discrete* ran;
ASSERT(ref);
- ran = CONTAINER_OF(ref, ssp_ran_type, ref);
- state = static_cast<ran_piecewise_linear_state*>(ran->state);
- if(state) {
- state->distrib->~piecewise_dist();
- MEM_RM(ran->allocator, state->distrib);
- MEM_RM(ran->allocator, state);
+ ran = CONTAINER_OF(ref, struct ssp_ranst_discrete, ref);
+ if(ran->distrib) {
+ ran->distrib->~discrete_dist();
+ MEM_RM(ran->allocator, ran->distrib);
}
+ darray_double_release(&ran->pdf);
MEM_RM(ran->allocator, ran);
}
/*******************************************************************************
+ * Discrete distribution exported functions
+ ******************************************************************************/
+
+res_T
+ssp_ranst_discrete_create
+ (struct mem_allocator* allocator,
+ struct ssp_ranst_discrete** out_ran)
+{
+ ssp_ranst_discrete* ran = nullptr;
+ void* mem = nullptr;
+ res_T res = RES_OK;
+
+ if(!out_ran) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ allocator = allocator ? allocator : &mem_default_allocator;
+
+ ran = static_cast<ssp_ranst_discrete*>(
+ MEM_CALLOC(allocator, 1, sizeof(ssp_ranst_discrete)));
+ if(!ran) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ ref_init(&ran->ref);
+
+ mem = MEM_ALLOC(allocator, sizeof(discrete_dist));
+ if (!mem) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ ran->allocator = allocator;
+ ran->distrib = static_cast<discrete_dist*>(new(mem) discrete_dist);
+ darray_double_init(allocator, &ran->pdf);
+
+exit:
+ if(out_ran) *out_ran = ran;
+ return res;
+error:
+ if(ran) {
+ SSP(ranst_discrete_ref_put(ran));
+ ran = NULL;
+ }
+ goto exit;
+}
+
+res_T
+ssp_ranst_discrete_setup
+ (struct ssp_ranst_discrete* ran,
+ const double* weights,
+ const size_t nweights)
+{
+ double* pdf;
+ double sum = 0;
+ size_t i;
+ res_T res = RES_OK;
+
+ if(!ran || !weights || !nweights) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ res = darray_double_resize(&ran->pdf, nweights);
+ if(res != RES_OK) goto error;
+
+ pdf = darray_double_data_get(&ran->pdf);
+
+ for(i = 0; i < nweights; i++) {
+ if(weights[i] < 0.f) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ sum += weights[i];
+ pdf[i] = weights[i];
+ }
+ if (sum == 0) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ sum = 1 / sum;
+ for(i = 0; i < nweights; i++) {
+ pdf[i] *= sum;
+ }
+
+ {
+ discrete_dist::param_type p{weights, weights + nweights};
+ ran->distrib->param(p);
+ }
+
+exit:
+ return res;
+error:
+ if(ran) {
+ darray_double_clear(&ran->pdf);
+ }
+ goto exit;
+}
+
+res_T
+ssp_ranst_discrete_ref_get(struct ssp_ranst_discrete* ran)
+{
+ if(!ran) return RES_BAD_ARG;
+ ref_get(&ran->ref);
+ return RES_OK;
+}
+
+res_T
+ssp_ranst_discrete_ref_put(struct ssp_ranst_discrete* ran)
+{
+ if(!ran) return RES_BAD_ARG;
+ ref_put(&ran->ref, ran_discrete_release);
+ return RES_OK;
+}
+
+size_t
+ssp_ranst_discrete_get(struct ssp_rng* rng, const struct ssp_ranst_discrete* ran)
+{
+ class rng_cxx r(*rng);
+ return ran->distrib->operator()(r);
+}
+
+double
+ssp_ranst_discrete_pdf(const size_t i, const struct ssp_ranst_discrete* ran)
+{
+ ASSERT(ran && i < darray_double_size_get(&ran->pdf));
+ return darray_double_cdata_get(&ran->pdf)[i];
+}
+
+/*******************************************************************************
* Gaussian distribution exported functions
******************************************************************************/
+
res_T
-ssp_ran_gaussian_create
+ssp_ranst_gaussian_create
(struct mem_allocator* allocator,
- struct ssp_ran_type** out_ran)
+ struct ssp_ranst_gaussian** out_ran)
{
- ssp_ran_type* ran = nullptr;
- ran_gaussian_state* state;
- void* tmp = nullptr;
+ ssp_ranst_gaussian* ran = nullptr;
+ void* mem = nullptr;
res_T res = RES_OK;
if (!out_ran)
@@ -86,61 +236,54 @@ ssp_ran_gaussian_create
allocator = allocator ? allocator : &mem_default_allocator;
- ran = static_cast<ssp_ran_type*>(MEM_CALLOC(allocator, 1, sizeof(ssp_ran_type)));
+ ran = static_cast<ssp_ranst_gaussian*>(
+ MEM_CALLOC(allocator, 1, sizeof(ssp_ranst_gaussian)));
if (!ran) {
res = RES_MEM_ERR;
goto err;
}
ref_init(&ran->ref);
- state = static_cast<ran_gaussian_state*>(
- MEM_CALLOC(allocator, 1, sizeof(ran_gaussian_state)));
- ran->state = state;
- if (!state) {
- res = RES_MEM_ERR;
- goto err;
- }
- tmp = MEM_ALLOC(allocator, sizeof(normal_dist));
- if (!tmp) {
+
+ mem = MEM_ALLOC(allocator, sizeof(normal_dist));
+ if (!mem) {
res = RES_MEM_ERR;
goto err;
}
ran->allocator = allocator;
- state->distrib = static_cast<normal_dist*>(new(tmp) normal_dist);
- state->K1 = -1; /* invalid */
+ ran->distrib = static_cast<normal_dist*>(new(mem) normal_dist);
+ ran->K1 = -1; /* invalid */
if (out_ran) *out_ran = ran;
end:
return res;
err:
if (ran) {
- ref_put(&ran->ref, gaussian_release);
+ SSP(ranst_gaussian_ref_put(ran));
ran = nullptr;
}
goto end;
}
res_T
-ssp_ran_gaussian_init
- (struct ssp_ran_type* ran,
+ssp_ranst_gaussian_setup
+ (struct ssp_ranst_gaussian* ran,
double mu,
double sigma)
{
- ran_gaussian_state* state;
if (!ran || sigma < 0)
return RES_BAD_ARG;
normal_dist::param_type p{mu, sigma};
- state = static_cast<ran_gaussian_state*>(ran->state);
- state->distrib->param(p);
- state->mu = mu;
- state->K1 = 1.0 / sigma;
- state->K2 = 1.0 / (sigma * SQRT_2_PI);
+ ran->distrib->param(p);
+ ran->mu = mu;
+ ran->K1 = 1.0 / sigma;
+ ran->K2 = 1.0 / (sigma * SQRT_2_PI);
return RES_OK;
}
res_T
-ssp_ran_gaussian_ref_get
- (struct ssp_ran_type* ran)
+ssp_ranst_gaussian_ref_get
+ (struct ssp_ranst_gaussian* ran)
{
if(!ran) return RES_BAD_ARG;
ref_get(&ran->ref);
@@ -148,8 +291,8 @@ ssp_ran_gaussian_ref_get
}
res_T
-ssp_ran_gaussian_ref_put
- (struct ssp_ran_type* ran)
+ssp_ranst_gaussian_ref_put
+ (struct ssp_ranst_gaussian* ran)
{
if(!ran) return RES_BAD_ARG;
ref_put(&ran->ref, gaussian_release);
@@ -157,23 +300,21 @@ ssp_ran_gaussian_ref_put
}
double
-ssp_ran_gaussian_get
- (const struct ssp_ran_type* ran, struct ssp_rng* rng)
+ssp_ranst_gaussian_get
+ (const struct ssp_ranst_gaussian* ran, struct ssp_rng* rng)
{
class rng_cxx r(*rng);
- ran_gaussian_state* state = static_cast<ran_gaussian_state*>(ran->state);
- ASSERT(state->K1 > 0);
- return state->distrib->operator()(r);
+ ASSERT(ran->K1 > 0);
+ return ran->distrib->operator()(r);
}
double
-ssp_ran_gaussian_pdf
- (const struct ssp_ran_type* ran, double x)
+ssp_ranst_gaussian_pdf
+ (const struct ssp_ranst_gaussian* ran, double x)
{
- ran_gaussian_state* state = static_cast<ran_gaussian_state*>(ran->state);
- const double tmp = (x - state->mu) * state->K1;
- ASSERT(state->K1 > 0);
- return state->K2 * exp(-0.5 * tmp * tmp);
+ const double tmp = (x - ran->mu) * ran->K1;
+ ASSERT(ran->K1 > 0);
+ return ran->K2 * exp(-0.5 * tmp * tmp);
}
/*******************************************************************************
@@ -181,13 +322,12 @@ ssp_ran_gaussian_pdf
******************************************************************************/
res_T
-ssp_ran_piecewise_linear_create
+ssp_ranst_piecewise_linear_create
(struct mem_allocator* allocator,
- struct ssp_ran_type** out_ran)
+ struct ssp_ranst_piecewise_linear** out_ran)
{
- ssp_ran_type* ran = nullptr;
- ran_piecewise_linear_state* state;
- void* tmp = nullptr;
+ ssp_ranst_piecewise_linear* ran = nullptr;
+ void* mem = nullptr;
res_T res = RES_OK;
if (!out_ran)
@@ -195,40 +335,35 @@ ssp_ran_piecewise_linear_create
allocator = allocator ? allocator : &mem_default_allocator;
- ran = static_cast<ssp_ran_type*>(MEM_CALLOC(allocator, 1, sizeof(ssp_ran_type)));
+ ran = static_cast<ssp_ranst_piecewise_linear*>(
+ MEM_CALLOC(allocator, 1, sizeof(ssp_ranst_piecewise_linear)));
if (!ran) {
res = RES_MEM_ERR;
goto err;
}
ref_init(&ran->ref);
- state = static_cast<ran_piecewise_linear_state*>(
- MEM_CALLOC(allocator, 1, sizeof(ran_piecewise_linear_state)));
- ran->state = state;
- if (!state) {
- res = RES_MEM_ERR;
- goto err;
- }
- tmp = MEM_ALLOC(allocator, sizeof(piecewise_dist));
- if (!tmp) {
+
+ mem = MEM_ALLOC(allocator, sizeof(piecewise_dist));
+ if (!mem) {
res = RES_MEM_ERR;
goto err;
}
ran->allocator = allocator;
- state->distrib = static_cast<piecewise_dist*>(new(tmp) piecewise_dist);
+ ran->distrib = static_cast<piecewise_dist*>(new(mem) piecewise_dist);
if (out_ran) *out_ran = ran;
end:
return res;
err:
if (ran) {
- ref_put(&ran->ref, piecewise_release);
+ SSP(ranst_piecewise_linear_ref_put(ran));
ran = nullptr;
}
goto end;
}
res_T
-ssp_ran_piecewise_linear_init
- (struct ssp_ran_type* ran,
+ssp_ranst_piecewise_linear_setup
+ (struct ssp_ranst_piecewise_linear* ran,
const double* intervals,
const double* weights,
size_t size)
@@ -237,16 +372,14 @@ ssp_ran_piecewise_linear_init
return RES_BAD_ARG;
piecewise_dist::param_type p{intervals, intervals + size, weights};
- ran_piecewise_linear_state* state
- = static_cast<ran_piecewise_linear_state*>(ran->state);
- state->distrib->param(p);
+ ran->distrib->param(p);
return RES_OK;
}
res_T
-ssp_ran_piecewise_linear_ref_get
- (struct ssp_ran_type* ran)
+ssp_ranst_piecewise_linear_ref_get
+ (struct ssp_ranst_piecewise_linear* ran)
{
if(!ran) return RES_BAD_ARG;
ref_get(&ran->ref);
@@ -254,8 +387,8 @@ ssp_ran_piecewise_linear_ref_get
}
res_T
-ssp_ran_piecewise_linear_ref_put
- (struct ssp_ran_type* ran)
+ssp_ranst_piecewise_linear_ref_put
+ (struct ssp_ranst_piecewise_linear* ran)
{
if(!ran) return RES_BAD_ARG;
ref_put(&ran->ref, piecewise_release);
@@ -263,12 +396,10 @@ ssp_ran_piecewise_linear_ref_put
}
double
-ssp_ran_piecewise_linear_get
- (const struct ssp_ran_type* ran,
+ssp_ranst_piecewise_linear_get
+ (const struct ssp_ranst_piecewise_linear* ran,
struct ssp_rng* rng)
{
class rng_cxx r(*rng);
- ran_piecewise_linear_state* state
- = static_cast<ran_piecewise_linear_state*>(ran->state);
- return state->distrib->operator()(r);
+ return ran->distrib->operator()(r);
}
diff --git a/src/ssp_ran_discrete.c b/src/ssp_ran_discrete.c
@@ -1,203 +0,0 @@
-/* Copyright (C) |Meso|Star> 2015-2016 (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,
- struct ssp_ran_discrete** out_ran)
-{
- struct mem_allocator* allocator = NULL;
- struct ssp_ran_discrete* ran = NULL;
- res_T res = RES_OK;
-
- if(!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);
-
-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_setup
- (struct ssp_ran_discrete* ran,
- const double* weights,
- const size_t nweights)
-{
- double* cdf;
- double* pdf;
- size_t i;
- res_T res = RES_OK;
-
- if(!ran || !weights || !nweights) {
- res = RES_BAD_ARG;
- goto error;
- }
-
- 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:
- return res;
-error:
- if(ran) {
- darray_double_clear(&ran->cdf);
- darray_double_clear(&ran->pdf);
- }
- 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 && darray_double_size_get(&ran->cdf) != 0);
-
- 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
@@ -40,7 +40,7 @@ int
main(int argc, char** argv)
{
struct mem_allocator allocator;
- struct ssp_ran_discrete* ran;
+ struct ssp_ranst_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);
@@ -52,32 +52,32 @@ main(int argc, char** 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), RES_BAD_ARG);
- CHECK(ssp_ran_discrete_create(NULL, &ran), RES_OK);
+ CHECK(ssp_ranst_discrete_create(NULL, NULL), RES_BAD_ARG);
+ CHECK(ssp_ranst_discrete_create(NULL, &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_ranst_discrete_ref_get(NULL), RES_BAD_ARG);
+ CHECK(ssp_ranst_discrete_ref_get(ran), RES_OK);
+ CHECK(ssp_ranst_discrete_ref_put(NULL), RES_BAD_ARG);
+ CHECK(ssp_ranst_discrete_ref_put(ran), RES_OK);
+ CHECK(ssp_ranst_discrete_ref_put(ran), RES_OK);
- CHECK(ssp_ran_discrete_create(&allocator, &ran), RES_OK);
+ CHECK(ssp_ranst_discrete_create(&allocator, &ran), RES_OK);
- CHECK(ssp_ran_discrete_setup(NULL, NULL, 0), RES_BAD_ARG);
- CHECK(ssp_ran_discrete_setup(ran, NULL, 0), RES_BAD_ARG);
- CHECK(ssp_ran_discrete_setup(NULL, weights, 0), RES_BAD_ARG);
- CHECK(ssp_ran_discrete_setup(ran, weights, 0), RES_BAD_ARG);
- CHECK(ssp_ran_discrete_setup(NULL, NULL, nweights), RES_BAD_ARG);
- CHECK(ssp_ran_discrete_setup(ran, NULL, nweights), RES_BAD_ARG);
- CHECK(ssp_ran_discrete_setup(NULL, weights, nweights), RES_BAD_ARG);
- CHECK(ssp_ran_discrete_setup(ran, weights, nweights), RES_OK);
+ CHECK(ssp_ranst_discrete_setup(NULL, NULL, 0), RES_BAD_ARG);
+ CHECK(ssp_ranst_discrete_setup(ran, NULL, 0), RES_BAD_ARG);
+ CHECK(ssp_ranst_discrete_setup(NULL, weights, 0), RES_BAD_ARG);
+ CHECK(ssp_ranst_discrete_setup(ran, weights, 0), RES_BAD_ARG);
+ CHECK(ssp_ranst_discrete_setup(NULL, NULL, nweights), RES_BAD_ARG);
+ CHECK(ssp_ranst_discrete_setup(ran, NULL, nweights), RES_BAD_ARG);
+ CHECK(ssp_ranst_discrete_setup(NULL, weights, nweights), RES_BAD_ARG);
+ CHECK(ssp_ranst_discrete_setup(ran, weights, nweights), RES_OK);
memset(tmp, 0, sizeof(tmp));
FOR_EACH(i, 0, NSAMPS) {
- const size_t k = ssp_ran_discrete(rng, ran);
+ const size_t k = ssp_ranst_discrete_get(rng, ran);
double pdf;
CHECK(k < nweights, 1);
- pdf = ssp_ran_discrete_pdf(k, ran);
+ pdf = ssp_ranst_discrete_pdf(k, ran);
++tmp[k];
CHECK(pdf, weights[k]);
CHECK(pdf >= 0.f, 1);
@@ -87,19 +87,19 @@ main(int argc, char** argv)
NCHECK(tmp[i], 0);
}
- CHECK(ssp_ran_discrete_setup(ran, weights, nweights-1), RES_OK);
+ CHECK(ssp_ranst_discrete_setup(ran, weights, nweights-1), RES_OK);
FOR_EACH(i, 0, NSAMPS) {
- const size_t k = ssp_ran_discrete(rng, ran);
+ const size_t k = ssp_ranst_discrete_get(rng, ran);
double pdf;
CHECK(k < nweights-1, 1);
- pdf = ssp_ran_discrete_pdf(k, ran);
+ pdf = ssp_ranst_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);
+ FOR_EACH(i, 0, nweights-1) accum += ssp_ranst_discrete_pdf(i, ran);
CHECK(eq_eps(accum, 1, 1.e-8), 1);
- CHECK(ssp_ran_discrete_ref_put(ran), RES_OK);
+ CHECK(ssp_ranst_discrete_ref_put(ran), RES_OK);
CHECK(ssp_rng_ref_put(rng), RES_OK);
check_memory_allocator(&allocator);
diff --git a/src/test_ssp_ran_gaussian.c b/src/test_ssp_ran_gaussian.c
@@ -40,7 +40,7 @@ main(int argc, char** argv)
struct ssp_rng_proxy* proxy;
struct ssp_rng* rng;
struct mem_allocator allocator;
- struct ssp_ran_type *gaussian;
+ struct ssp_ranst_gaussian *gaussian;
int i;
double x = 0, x2 = 0;
double exp_mean = 10, mean;
@@ -54,26 +54,26 @@ main(int argc, char** argv)
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_gaussian_create(NULL, NULL), RES_BAD_ARG);
- CHECK(ssp_ran_gaussian_create(NULL, &gaussian), RES_OK);
- CHECK(ssp_ran_gaussian_ref_put(gaussian), RES_OK);
+ CHECK(ssp_ranst_gaussian_create(NULL, NULL), RES_BAD_ARG);
+ CHECK(ssp_ranst_gaussian_create(NULL, &gaussian), RES_OK);
+ CHECK(ssp_ranst_gaussian_ref_put(gaussian), RES_OK);
- CHECK(ssp_ran_gaussian_create(&allocator, NULL), RES_BAD_ARG);
- CHECK(ssp_ran_gaussian_create(&allocator, &gaussian), RES_OK);
+ CHECK(ssp_ranst_gaussian_create(&allocator, NULL), RES_BAD_ARG);
+ CHECK(ssp_ranst_gaussian_create(&allocator, &gaussian), RES_OK);
- CHECK(ssp_ran_gaussian_init(NULL, exp_mean, exp_std), RES_BAD_ARG);
- CHECK(ssp_ran_gaussian_init(gaussian, exp_mean, -1), RES_BAD_ARG);
- CHECK(ssp_ran_gaussian_init(gaussian, exp_mean, exp_std), RES_OK);
+ CHECK(ssp_ranst_gaussian_setup(NULL, exp_mean, exp_std), RES_BAD_ARG);
+ CHECK(ssp_ranst_gaussian_setup(gaussian, exp_mean, -1), RES_BAD_ARG);
+ CHECK(ssp_ranst_gaussian_setup(gaussian, exp_mean, exp_std), RES_OK);
- CHECK(ssp_ran_gaussian_ref_get(NULL), RES_BAD_ARG);
- CHECK(ssp_ran_gaussian_ref_get(gaussian), RES_OK);
+ CHECK(ssp_ranst_gaussian_ref_get(NULL), RES_BAD_ARG);
+ CHECK(ssp_ranst_gaussian_ref_get(gaussian), RES_OK);
- CHECK(ssp_ran_gaussian_ref_put(NULL), RES_BAD_ARG);
- CHECK(ssp_ran_gaussian_ref_put(gaussian), RES_OK);
+ CHECK(ssp_ranst_gaussian_ref_put(NULL), RES_BAD_ARG);
+ CHECK(ssp_ranst_gaussian_ref_put(gaussian), RES_OK);
time_current(&start);
for (i = 0; i < NBS; i++) {
- double _x = ssp_ran_gaussian_get(gaussian, rng);
+ double _x = ssp_ranst_gaussian_get(gaussian, rng);
x += _x;
x2 += _x * _x;
}
@@ -109,7 +109,7 @@ main(int argc, char** argv)
CHECK(fabs(mean - exp_mean) < 0.003, 1);
CHECK(fabs(std - exp_std) < 0.001, 1);
- CHECK(ssp_ran_gaussian_ref_put(gaussian), RES_OK);
+ CHECK(ssp_ranst_gaussian_ref_put(gaussian), RES_OK);
CHECK(ssp_rng_ref_put(rng), RES_OK);
CHECK(ssp_rng_proxy_ref_put(proxy), RES_OK);
diff --git a/src/test_ssp_ran_piecewise_linear.c b/src/test_ssp_ran_piecewise_linear.c
@@ -39,7 +39,7 @@ main(int argc, char** argv)
struct ssp_rng_proxy* proxy;
struct ssp_rng* rng;
struct mem_allocator allocator;
- struct ssp_ran_type *pwl;
+ struct ssp_ranst_piecewise_linear *pwl;
int i;
double exp_mean = 5, mean;
double exp_std = 10 / sqrt(12) /*sqrt((b - a)² / 12) */, std;
@@ -53,27 +53,27 @@ main(int argc, char** argv)
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, NULL), RES_BAD_ARG);
- CHECK(ssp_ran_piecewise_linear_create(NULL, &pwl), RES_OK);
- CHECK(ssp_ran_piecewise_linear_ref_put(pwl), RES_OK);
+ CHECK(ssp_ranst_piecewise_linear_create(NULL, NULL), RES_BAD_ARG);
+ CHECK(ssp_ranst_piecewise_linear_create(NULL, &pwl), RES_OK);
+ CHECK(ssp_ranst_piecewise_linear_ref_put(pwl), RES_OK);
- CHECK(ssp_ran_piecewise_linear_create(&allocator, NULL), RES_BAD_ARG);
- CHECK(ssp_ran_piecewise_linear_create(&allocator, &pwl), RES_OK);
+ CHECK(ssp_ranst_piecewise_linear_create(&allocator, NULL), RES_BAD_ARG);
+ CHECK(ssp_ranst_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_ranst_piecewise_linear_setup(NULL, intervals, weights, sizeof(intervals)/sizeof(double)), RES_BAD_ARG);
+ CHECK(ssp_ranst_piecewise_linear_setup(pwl, NULL, weights, sizeof(intervals)/sizeof(double)), RES_BAD_ARG);
+ CHECK(ssp_ranst_piecewise_linear_setup(pwl, intervals, NULL, sizeof(intervals)/sizeof(double)), RES_BAD_ARG);
+ CHECK(ssp_ranst_piecewise_linear_setup(pwl, intervals, weights, 1), RES_BAD_ARG);
+ CHECK(ssp_ranst_piecewise_linear_setup(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_ranst_piecewise_linear_ref_get(NULL), RES_BAD_ARG);
+ CHECK(ssp_ranst_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_ranst_piecewise_linear_ref_put(NULL), RES_BAD_ARG);
+ CHECK(ssp_ranst_piecewise_linear_ref_put(pwl), RES_OK);
for (i = 0; i < NBS; i++) {
- double _x = ssp_ran_piecewise_linear_get(pwl, rng);
+ double _x = ssp_ranst_piecewise_linear_get(pwl, rng);
CHECK(0 <= _x && _x <= 10, 1);
x += _x;
x2 += _x * _x;
@@ -85,7 +85,7 @@ main(int argc, char** argv)
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_ranst_piecewise_linear_ref_put(pwl), RES_OK);
CHECK(ssp_rng_ref_put(rng), RES_OK);
CHECK(ssp_rng_proxy_ref_put(proxy), RES_OK);
