star-sp

Random number generators and distributions
git clone git://git.meso-star.fr/star-sp.git
Log | Files | Refs | README | LICENSE

commit b6cb3e32af3963739a363c265f196fda4f494748
parent 5a90e5071b7981ee82146989272cd1c382c5431f
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date:   Thu,  7 Jul 2016 10:14:06 +0200

Some refactoring; mainly to decrease code duplication.

Diffstat:
Mcmake/CMakeLists.txt | 3+--
Msrc/ssp.h | 25++++++++++++-------------
Msrc/ssp_distributions.cpp | 127++++++++++++++++++++++++++++++++++++++++++++++---------------------------------
Dsrc/ssp_distributions.h | 27---------------------------
Msrc/test_ssp_ran_gaussian.c | 2+-
Msrc/test_ssp_ran_piecewise_linear.c | 2+-
6 files changed, 89 insertions(+), 97 deletions(-)

diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -76,8 +76,7 @@ set(SSP_FILES_SRC ssp_rng_proxy.c ssp_distributions.cpp) set(SSP_FILES_INC - ssp_rng_c.h - ssp_distributions.h) + ssp_rng_c.h) set(SSP_FILES_INC_API ssp.h) set(SSP_FILES_DOC diff --git a/src/ssp.h b/src/ssp.h @@ -58,8 +58,7 @@ struct ssp_rng; struct ssp_rng_proxy; struct ssp_ran_discrete; -struct ssp_ran_gaussian; -struct ssp_ran_piecewise_linear; +struct ssp_ran_type; /* Forward declaration of external types */ struct mem_allocator; @@ -586,31 +585,31 @@ ssp_ran_sphere_hg SSP_API res_T ssp_ran_gaussian_create (struct mem_allocator* allocator, /* May be NULL <=> Use default allocator */ - struct ssp_ran_gaussian** ran); + struct ssp_ran_type** ran); SSP_API res_T ssp_ran_gaussian_init - (struct ssp_ran_gaussian* ran, + (struct ssp_ran_type* ran, double mu, double sigma); SSP_API res_T ssp_ran_gaussian_ref_get - (struct ssp_ran_gaussian* ran); + (struct ssp_ran_type* ran); SSP_API res_T ssp_ran_gaussian_ref_put - (struct ssp_ran_gaussian* ran); + (struct ssp_ran_type* ran); SSP_API res_T ssp_ran_gaussian_get - (const struct ssp_ran_gaussian* ran, + (const struct ssp_ran_type* ran, struct ssp_rng* rng, double* res); SSP_API res_T ssp_distribution_gaussian_pdf - (const struct ssp_ran_gaussian* ran, + (const struct ssp_ran_type* ran, double x, double* res); @@ -621,26 +620,26 @@ ssp_distribution_gaussian_pdf 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); + struct ssp_ran_type **ran); SSP_API res_T ssp_ran_piecewise_linear_init - (struct ssp_ran_piecewise_linear *ran, + (struct ssp_ran_type *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); + (struct ssp_ran_type* ran); SSP_API res_T ssp_ran_piecewise_linear_ref_put - (struct ssp_ran_piecewise_linear* ran); + (struct ssp_ran_type* ran); SSP_API res_T ssp_ran_piecewise_linear_get - (const struct ssp_ran_piecewise_linear *ran, + (const struct ssp_ran_type *ran, struct ssp_rng* rng, double *res); diff --git a/src/ssp_distributions.cpp b/src/ssp_distributions.cpp @@ -1,9 +1,16 @@ -#include "ssp_distributions.h" - #include <rsys/mem_allocator.h> +#include <rsys/ref_count.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 ******************************************************************************/ @@ -12,7 +19,7 @@ using normal_dist=RAN_NAMESPACE::normal_distribution<double>; struct ran_gaussian_state { - normal_dist *distrib; + normal_dist* distrib; double mu; double K1; // 1.0 / sigma; double K2; // 1.0 / (sigma * SQRT_2_PI) @@ -21,13 +28,15 @@ struct ran_gaussian_state static void gaussian_release(ref_T* ref) { - struct ssp_ran_gaussian* ran; + ssp_ran_type* ran; + ran_gaussian_state* state; ASSERT(ref); - ran = CONTAINER_OF(ref, struct ssp_ran_gaussian, ref); - if(ran->state) { - ran->state->distrib->~normal_dist(); - MEM_RM(ran->allocator, ran->state->distrib); - MEM_RM(ran->allocator, ran->state); + 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); } MEM_RM(ran->allocator, ran); } @@ -40,19 +49,21 @@ using piecewise_dist=RAN_NAMESPACE::piecewise_linear_distribution<double>; struct ran_piecewise_linear_state { - piecewise_dist *distrib; + piecewise_dist* distrib; }; static void piecewise_release(ref_T* ref) { - ssp_ran_piecewise_linear* ran; + ssp_ran_type* ran; + ran_piecewise_linear_state* state; 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); + 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); } MEM_RM(ran->allocator, ran); } @@ -63,10 +74,11 @@ piecewise_release(ref_T* ref) res_T ssp_ran_gaussian_create (struct mem_allocator* allocator, - struct ssp_ran_gaussian** out_ran) + struct ssp_ran_type** out_ran) { - struct ssp_ran_gaussian *ran = nullptr; - void *tmp = nullptr; + ssp_ran_type* ran = nullptr; + ran_gaussian_state* state; + void* tmp = nullptr; res_T res = RES_OK; if (!out_ran) @@ -74,16 +86,16 @@ ssp_ran_gaussian_create allocator = allocator ? allocator : &mem_default_allocator; - ran = static_cast<struct ssp_ran_gaussian*>( - MEM_CALLOC(allocator, 1, sizeof(struct ssp_ran_gaussian))); + ran = static_cast<ssp_ran_type*>(MEM_CALLOC(allocator, 1, sizeof(ssp_ran_type))); if (!ran) { res = RES_MEM_ERR; goto err; } ref_init(&ran->ref); - ran->state = static_cast<ran_gaussian_state*>( + state = static_cast<ran_gaussian_state*>( MEM_CALLOC(allocator, 1, sizeof(ran_gaussian_state))); - if (!ran->state) { + ran->state = state; + if (!state) { res = RES_MEM_ERR; goto err; } @@ -93,10 +105,10 @@ ssp_ran_gaussian_create goto err; } ran->allocator = allocator; - ran->state->distrib = static_cast<normal_dist*>(new(tmp) normal_dist); - ran->state->mu = 0; - ran->state->K1 = 1; - ran->state->K2 = 1 / SQRT_2_PI; + state->distrib = static_cast<normal_dist*>(new(tmp) normal_dist); + state->mu = 0; + state->K1 = 1; + state->K2 = 1 / SQRT_2_PI; if (out_ran) *out_ran = ran; end: return res; @@ -110,25 +122,27 @@ ssp_ran_gaussian_create res_T ssp_ran_gaussian_init - (struct ssp_ran_gaussian* ran, + (struct ssp_ran_type* ran, double mu, double sigma) { + ran_gaussian_state* state; if (!ran || sigma < 0) return RES_BAD_ARG; normal_dist::param_type p{mu, sigma}; - ran->state->distrib->param(p); - ran->state->mu = mu; - ran->state->K1 = 1.0 / sigma; - ran->state->K2 = 1.0 / (sigma * SQRT_2_PI); + 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); return RES_OK; } res_T ssp_ran_gaussian_ref_get - (struct ssp_ran_gaussian* ran) + (struct ssp_ran_type* ran) { if(!ran) return RES_BAD_ARG; ref_get(&ran->ref); @@ -137,7 +151,7 @@ ssp_ran_gaussian_ref_get res_T ssp_ran_gaussian_ref_put - (struct ssp_ran_gaussian* ran) + (struct ssp_ran_type* ran) { if(!ran) return RES_BAD_ARG; ref_put(&ran->ref, gaussian_release); @@ -146,21 +160,23 @@ ssp_ran_gaussian_ref_put res_T ssp_ran_gaussian_get - (const struct ssp_ran_gaussian *ran, struct ssp_rng* rng, double* res) + (const struct ssp_ran_type* ran, struct ssp_rng* rng, double* res) { if(!ran || !rng || !res) return RES_BAD_ARG; class rng_cxx r(*rng); - *res = ran->state->distrib->operator()(r); + ran_gaussian_state* state = static_cast<ran_gaussian_state*>(ran->state); + *res = state->distrib->operator()(r); return RES_OK; } res_T ssp_ran_gaussian_pdf - (const struct ssp_ran_gaussian *ran, double x, double* res) + (const struct ssp_ran_type* ran, double x, double* res) { if(!ran || !res) return RES_BAD_ARG; - const double tmp = (x - ran->state->mu) * ran->state->K1; - *res = ran->state->K2 * exp(-0.5 * tmp * tmp); + ran_gaussian_state* state = static_cast<ran_gaussian_state*>(ran->state); + const double tmp = (x - state->mu) * state->K1; + *res = state->K2 * exp(-0.5 * tmp * tmp); return RES_OK; } @@ -171,10 +187,11 @@ ssp_ran_gaussian_pdf res_T ssp_ran_piecewise_linear_create (struct mem_allocator* allocator, - struct ssp_ran_piecewise_linear **out_ran) + struct ssp_ran_type** out_ran) { - ssp_ran_piecewise_linear *ran = nullptr; - void *tmp = nullptr; + ssp_ran_type* ran = nullptr; + ran_piecewise_linear_state* state; + void* tmp = nullptr; res_T res = RES_OK; if (!out_ran) @@ -182,16 +199,16 @@ ssp_ran_piecewise_linear_create allocator = allocator ? allocator : &mem_default_allocator; - ran = static_cast<ssp_ran_piecewise_linear*>( - MEM_CALLOC(allocator, 1, sizeof(ssp_ran_piecewise_linear))); + ran = static_cast<ssp_ran_type*>(MEM_CALLOC(allocator, 1, sizeof(ssp_ran_type))); if (!ran) { res = RES_MEM_ERR; goto err; } ref_init(&ran->ref); - ran->state = static_cast<ran_piecewise_linear_state*>( + state = static_cast<ran_piecewise_linear_state*>( MEM_CALLOC(allocator, 1, sizeof(ran_piecewise_linear_state))); - if (!ran->state) { + ran->state = state; + if (!state) { res = RES_MEM_ERR; goto err; } @@ -201,7 +218,7 @@ ssp_ran_piecewise_linear_create goto err; } ran->allocator = allocator; - ran->state->distrib = static_cast<piecewise_dist*>(new(tmp) piecewise_dist); + state->distrib = static_cast<piecewise_dist*>(new(tmp) piecewise_dist); if (out_ran) *out_ran = ran; end: return res; @@ -215,7 +232,7 @@ ssp_ran_piecewise_linear_create res_T ssp_ran_piecewise_linear_init - (struct ssp_ran_piecewise_linear *ran, + (struct ssp_ran_type* ran, const double* intervals, const double* weights, size_t size) @@ -224,14 +241,16 @@ ssp_ran_piecewise_linear_init return RES_BAD_ARG; piecewise_dist::param_type p{intervals, intervals + size, weights}; - ran->state->distrib->param(p); + ran_piecewise_linear_state* state + = static_cast<ran_piecewise_linear_state*>(ran->state); + state->distrib->param(p); return RES_OK; } res_T ssp_ran_piecewise_linear_ref_get - (struct ssp_ran_piecewise_linear* ran) + (struct ssp_ran_type* ran) { if(!ran) return RES_BAD_ARG; ref_get(&ran->ref); @@ -240,7 +259,7 @@ ssp_ran_piecewise_linear_ref_get res_T ssp_ran_piecewise_linear_ref_put - (struct ssp_ran_piecewise_linear* ran) + (struct ssp_ran_type* ran) { if(!ran) return RES_BAD_ARG; ref_put(&ran->ref, piecewise_release); @@ -249,12 +268,14 @@ ssp_ran_piecewise_linear_ref_put res_T ssp_ran_piecewise_linear_get - (const struct ssp_ran_piecewise_linear *ran, + (const struct ssp_ran_type* ran, struct ssp_rng* rng, double* res) { if(!ran || !rng || !res) return RES_BAD_ARG; class rng_cxx r(*rng); - *res = ran->state->distrib->operator()(r); + ran_piecewise_linear_state* state + = static_cast<ran_piecewise_linear_state*>(ran->state); + *res = state->distrib->operator()(r); return RES_OK; } diff --git a/src/ssp_distributions.h b/src/ssp_distributions.h @@ -1,27 +0,0 @@ -#ifndef __SSP_DISTRIB_GAUSSIAN_H__ -#define __SSP_DISTRIB_GAUSSIAN_H__ - -#include <rsys/rsys.h> -#include <rsys/ref_count.h> - -/* Forward declaration of opaque types */ -struct ran_gaussian_state; -struct ran_piecewise_linear_state; - -/* Forward declaration of external types */ -struct ssp_rng; -struct mem_allocator; - -struct ssp_ran_gaussian { - ref_T ref; - struct mem_allocator* allocator; - struct ran_gaussian_state *state; -}; - -struct ssp_ran_piecewise_linear { - ref_T ref; - struct mem_allocator* allocator; - struct ran_piecewise_linear_state *state; -}; - -#endif 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_gaussian *gaussian; + struct ssp_ran_type *gaussian; int i; double x = 0, x2 = 0; double exp_mean = 10, mean; 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_piecewise_linear *pwl; + struct ssp_ran_type *pwl; int i; double exp_mean = 5, mean; double exp_std = 10 / sqrt(12) /*sqrt((b - a)² / 12) */, std;