commit 99643458cecabf51beaf15721ba7d104b121b9ac
parent ae7e4d277f18694008077cb9a5230f5ab361c84b
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Mon, 4 Jul 2016 17:57:29 +0200
Rename *_distribution_* stuff into *_ran_*.
Split gaussian distribution creation into create and init too.
Diffstat:
6 files changed, 211 insertions(+), 177 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -155,7 +155,7 @@ if(NOT NO_TEST)
new_test(test_ssp_ran_hg ${MATH_LIB})
new_test(test_ssp_ran_triangle ${MATH_LIB})
new_test(test_ssp_rng_proxy)
- new_test(test_ssp_distrib_gaussian)
+ new_test(test_ssp_ran_gaussian)
endif()
################################################################################
diff --git a/src/ssp.h b/src/ssp.h
@@ -58,7 +58,7 @@
struct ssp_rng;
struct ssp_rng_proxy;
struct ssp_ran_discrete;
-struct ssp_distrib_gaussian;
+struct ssp_ran_gaussian;
/* Forward declaration of external types */
struct mem_allocator;
@@ -583,27 +583,31 @@ ssp_ran_sphere_hg
******************************************************************************/
SSP_API res_T
-ssp_distribution_gaussian_init
-(struct mem_allocator* allocator,
- double mu,
- double sigma,
- struct ssp_distrib_gaussian **out_dist);
+ssp_ran_gaussian_create
+ (struct mem_allocator* allocator, /* May be NULL <=> Use default allocator */
+ struct ssp_ran_gaussian** ran);
+
+SSP_API res_T
+ssp_ran_gaussian_init
+ (struct ssp_ran_gaussian* ran,
+ double mu,
+ double sigma);
SSP_API res_T
-ssp_distribution_gaussian_ref_get
-(struct ssp_distrib_gaussian* distrib);
+ssp_ran_gaussian_ref_get
+ (struct ssp_ran_gaussian* ran);
SSP_API res_T
-ssp_distribution_gaussian_ref_put
-(struct ssp_distrib_gaussian* distrib);
+ssp_ran_gaussian_ref_put
+ (struct ssp_ran_gaussian* ran);
SSP_API double
-ssp_distribution_gaussian_get
-(const struct ssp_distrib_gaussian *dist, struct ssp_rng* rng);
+ssp_ran_gaussian_get
+ (const struct ssp_ran_gaussian* ran, struct ssp_rng* rng);
SSP_API double
ssp_distribution_gaussian_pdf
-(const struct ssp_distrib_gaussian *dist, double x);
+ (const struct ssp_ran_gaussian* ran, double x);
END_DECLS
diff --git a/src/ssp_distributions.cpp b/src/ssp_distributions.cpp
@@ -4,115 +4,130 @@
#include <rsys/math.h>
#include "ssp_rng_c.h"
-using normal_dist_double=RAN_NAMESPACE::normal_distribution<double>;
+/*******************************************************************************
+ * Gaussian distribution utilities
+ ******************************************************************************/
-struct ssp_distrib_gaussian_state
+using normal_dist=RAN_NAMESPACE::normal_distribution<double>;
+
+struct ran_gaussian_state
{
- normal_dist_double *distribution;
+ normal_dist *distrib;
double mu;
double K1; // 1.0 / sigma;
double K2; // 1.0 / (sigma * SQRT_2_PI)
};
static void
-distrib_gaussian_release(ref_T* ref)
+gaussian_release(ref_T* ref)
{
- ssp_distrib_gaussian* distrib;
+ struct ssp_ran_gaussian* ran;
ASSERT(ref);
- distrib = CONTAINER_OF(ref, struct ssp_distrib_gaussian, ref);
- if(distrib->state) {
- distrib->state->distribution->~normal_distribution<double>();
- MEM_RM(distrib->allocator, distrib->state->distribution);
- MEM_RM(distrib->allocator, distrib->state);
+ 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);
}
- MEM_RM(distrib->allocator, distrib);
+ MEM_RM(ran->allocator, ran);
}
/*******************************************************************************
* Gaussian distribution exported functions
******************************************************************************/
-
res_T
-ssp_distribution_gaussian_init
-(struct mem_allocator* allocator,
- double mu,
- double sigma,
- struct ssp_distrib_gaussian **out_dist)
+ssp_ran_gaussian_create
+ (struct mem_allocator* allocator,
+ struct ssp_ran_gaussian** out_ran)
{
- ssp_distrib_gaussian *dist = nullptr;
+ struct ssp_ran_gaussian *ran = nullptr;
void *tmp = nullptr;
res_T res = RES_OK;
- if (sigma < 0)
- return RES_BAD_ARG;
-
allocator = allocator ? allocator : &mem_default_allocator;
- dist = static_cast<ssp_distrib_gaussian*>(
- MEM_CALLOC(allocator, 1, sizeof(ssp_distrib_gaussian)));
- if (!dist) {
+ ran = static_cast<struct ssp_ran_gaussian*>(
+ MEM_CALLOC(allocator, 1, sizeof(struct ssp_ran_gaussian)));
+ if (!ran) {
res = RES_MEM_ERR;
goto err;
}
- ref_init(&dist->ref);
- dist->state = static_cast<ssp_distrib_gaussian_state*>(
- MEM_CALLOC(allocator, 1, sizeof(ssp_distrib_gaussian_state)));
- if (!dist->state) {
+ ref_init(&ran->ref);
+ ran->state = static_cast<ran_gaussian_state*>(
+ MEM_CALLOC(allocator, 1, sizeof(ran_gaussian_state)));
+ if (!ran->state) {
res = RES_MEM_ERR;
goto err;
}
- tmp = MEM_ALLOC(allocator, sizeof(normal_dist_double));
+ tmp = MEM_ALLOC(allocator, sizeof(normal_dist));
if (!tmp) {
res = RES_MEM_ERR;
goto err;
}
- dist->allocator = allocator;
- dist->state->distribution = static_cast<normal_dist_double*>(
- new(tmp) normal_dist_double(mu, sigma));
- dist->state->mu = mu;
- dist->state->K1 = 1.0 / sigma;
- dist->state->K2 = 1.0 / (sigma * SQRT_2_PI);
- if (out_dist) *out_dist = dist;
+ 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;
+ if (out_ran) *out_ran = ran;
end:
return res;
err:
- if (dist) {
- ref_put(&dist->ref, distrib_gaussian_release);
- dist = nullptr;
+ if (ran) {
+ ref_put(&ran->ref, gaussian_release);
+ ran = nullptr;
}
goto end;
}
res_T
-ssp_distribution_gaussian_ref_get
-(struct ssp_distrib_gaussian* distrib)
+ssp_ran_gaussian_init
+ (struct ssp_ran_gaussian* ran,
+ double mu,
+ double sigma)
+{
+ 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);
+
+ return RES_OK;
+}
+
+res_T
+ssp_ran_gaussian_ref_get
+ (struct ssp_ran_gaussian* ran)
{
- if(!distrib) return RES_BAD_ARG;
- ref_get(&distrib->ref);
+ if(!ran) return RES_BAD_ARG;
+ ref_get(&ran->ref);
return RES_OK;
}
res_T
-ssp_distribution_gaussian_ref_put
-(struct ssp_distrib_gaussian* distrib)
+ssp_ran_gaussian_ref_put
+ (struct ssp_ran_gaussian* ran)
{
- if(!distrib) return RES_BAD_ARG;
- ref_put(&distrib->ref, distrib_gaussian_release);
+ if(!ran) return RES_BAD_ARG;
+ ref_put(&ran->ref, gaussian_release);
return RES_OK;
}
double
-ssp_distribution_gaussian_get
-(const struct ssp_distrib_gaussian *dist, struct ssp_rng* rng)
+ssp_ran_gaussian_get
+ (const struct ssp_ran_gaussian *ran, struct ssp_rng* rng)
{
class rng_cxx r(*rng);
- return dist->state->distribution->operator()(r);
+ return ran->state->distrib->operator()(r);
}
double
-ssp_distribution_gaussian_pdf
-(const struct ssp_distrib_gaussian *dist, double x)
+ssp_ran_gaussian_pdf
+ (const struct ssp_ran_gaussian *ran, double x)
{
- const double tmp = (x - dist->state->mu) * dist->state->K1;
- return dist->state->K2 * exp(-0.5 * tmp * tmp);
+ const double tmp = (x - ran->state->mu) * ran->state->K1;
+ return ran->state->K2 * exp(-0.5 * tmp * tmp);
}
diff --git a/src/ssp_distributions.h b/src/ssp_distributions.h
@@ -5,18 +5,18 @@
#include <rsys/ref_count.h>
/* Forward declaration of opaque types */
-struct ssp_distrib_gaussian_state;
+struct ran_gaussian_state;
/* Forward declaration of external types */
struct ssp_rng;
struct mem_allocator;
-struct ssp_distrib_gaussian {
+struct ssp_ran_gaussian {
size_t sizeof_params;
size_t sizeof_state;
ref_T ref;
struct mem_allocator* allocator;
- struct ssp_distrib_gaussian_state *state;
+ struct ran_gaussian_state *state;
};
#endif
diff --git a/src/test_ssp_distrib_gaussian.c b/src/test_ssp_distrib_gaussian.c
@@ -1,106 +0,0 @@
-/* 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"
-#include <rsys/clock_time.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_distrib_gaussian *gaussian;
- int i;
- double x = 0, x2 = 0;
- double exp_mean = 10, mean;
- double exp_std = 3, std;
- struct time start, end, res;
- char dump[512];
- (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_distribution_gaussian_init(&allocator, exp_mean, exp_std, &gaussian), RES_OK);
-
- time_current(&start);
- for (i = 0; i < NBS; i++) {
- double _x = ssp_distribution_gaussian_get(gaussian, rng);
- x += _x;
- x2 += _x * _x;
- }
- time_current(&end);
- time_sub(&res, &end, &start);
- time_dump(&res, TIME_SEC|TIME_MSEC|TIME_USEC, NULL, dump, sizeof(dump));
- printf("%s--\n", dump);
-
- mean = x/NBS;
- std = sqrt(x2/NBS - x/NBS*x/NBS);
- printf("%g %g\n", mean, std);
- CHECK(fabs(mean - exp_mean) < 0.003, 1);
- CHECK(fabs(std - exp_std) < 0.001, 1);
-
- /* same test with inline gaussian generation */
- x = 0;
- x2 = 0;
-
- time_current(&start);
- for (i = 0; i < NBS; i++) {
- double _x = ssp_ran_gaussian(rng, exp_mean, exp_std);
- x += _x;
- x2 += _x * _x;
- }
- time_current(&end);
- time_sub(&res, &end, &start);
- time_dump(&res, TIME_SEC|TIME_MSEC|TIME_USEC, NULL, dump, sizeof(dump));
- printf("%s--\n", dump);
-
- mean = x/NBS;
- std = sqrt(x2/NBS - x/NBS*x/NBS);
- printf("%g %g\n", mean, std);
- CHECK(fabs(mean - exp_mean) < 0.003, 1);
- CHECK(fabs(std - exp_std) < 0.001, 1);
-
- CHECK(ssp_distribution_gaussian_ref_put(gaussian), 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;
-}
diff --git a/src/test_ssp_ran_gaussian.c b/src/test_ssp_ran_gaussian.c
@@ -0,0 +1,121 @@
+/* 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"
+#include <rsys/clock_time.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_gaussian *gaussian;
+ int i;
+ double x = 0, x2 = 0;
+ double exp_mean = 10, mean;
+ double exp_std = 3, std;
+ struct time start, end, res;
+ char dump[512];
+ (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_gaussian_create(NULL, &gaussian), RES_BAD_ARG);
+ CHECK(ssp_ran_gaussian_create(&allocator, NULL), RES_BAD_ARG);
+ CHECK(ssp_ran_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_ran_gaussian_ref_get(NULL), RES_BAD_ARG);
+ CHECK(ssp_ran_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_ran_gaussian_get(NULL, rng), RES_BAD_ARG);
+ CHECK(ssp_ran_gaussian_get(gaussian, NULL), RES_BAD_ARG);
+
+ time_current(&start);
+ for (i = 0; i < NBS; i++) {
+ double _x = ssp_ran_gaussian_get(gaussian, rng);
+ x += _x;
+ x2 += _x * _x;
+ }
+ time_current(&end);
+ time_sub(&res, &end, &start);
+ time_dump(&res, TIME_SEC|TIME_MSEC|TIME_USEC, NULL, dump, sizeof(dump));
+ printf("%s--\n", dump);
+
+ mean = x/NBS;
+ std = sqrt(x2/NBS - x/NBS*x/NBS);
+ printf("%g %g\n", mean, std);
+ CHECK(fabs(mean - exp_mean) < 0.003, 1);
+ CHECK(fabs(std - exp_std) < 0.001, 1);
+
+ /* same test with inline gaussian generation */
+ x = 0;
+ x2 = 0;
+
+ time_current(&start);
+ for (i = 0; i < NBS; i++) {
+ double _x = ssp_ran_gaussian(rng, exp_mean, exp_std);
+ x += _x;
+ x2 += _x * _x;
+ }
+ time_current(&end);
+ time_sub(&res, &end, &start);
+ time_dump(&res, TIME_SEC|TIME_MSEC|TIME_USEC, NULL, dump, sizeof(dump));
+ printf("%s--\n", dump);
+
+ mean = x/NBS;
+ std = sqrt(x2/NBS - x/NBS*x/NBS);
+ printf("%g %g\n", mean, std);
+ 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_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;
+}
