commit ae7e4d277f18694008077cb9a5230f5ab361c84b
parent 9d5495061132471fcece2919d05d993b5b51cb3a
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Mon, 4 Jul 2016 16:19:13 +0200
Starting version 0.4, that is mainly about distributions.
This first commit adds a gaussian/normal distribution with its test.
Diffstat:
7 files changed, 326 insertions(+), 41 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -66,14 +66,24 @@ rcmake_append_runtime_dirs(_runtime_dirs RSys ${Boost_LIBRARY_DIRS})
# Configure and define targets
################################################################################
set(VERSION_MAJOR 0)
-set(VERSION_MINOR 3)
+set(VERSION_MINOR 4)
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)
-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)
+set(SSP_FILES_SRC
+ ssp_ran_discrete.c
+ ssp_rng.c
+ ssp_rng_proxy.c
+ ssp_distributions.cpp)
+set(SSP_FILES_INC
+ ssp_rng_c.h
+ ssp_distributions.h)
+set(SSP_FILES_INC_API
+ ssp.h)
+set(SSP_FILES_DOC
+ COPYING.en
+ COPYING.fr
+ README.md)
# Prepend each file in the `SSP_FILES_<SRC|INC>' list by `SSP_SOURCE_DIR'
rcmake_prepend_path(SSP_FILES_SRC ${SSP_SOURCE_DIR})
@@ -108,7 +118,7 @@ if(NOT NO_TEST)
add_executable(${_name}
${SSP_SOURCE_DIR}/${_name}.c
${SSP_SOURCE_DIR}/test_ssp_utils.h)
- target_link_libraries(${_name} ssp RSys)
+ target_link_libraries(${_name} ssp RSys ${MATH_LIB})
set(_libraries ${ARGN})
foreach(_lib ${_libraries})
target_link_libraries(${_name} ${_lib})
@@ -145,6 +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)
endif()
################################################################################
diff --git a/src/ssp.h b/src/ssp.h
@@ -58,6 +58,7 @@
struct ssp_rng;
struct ssp_rng_proxy;
struct ssp_ran_discrete;
+struct ssp_distrib_gaussian;
/* Forward declaration of external types */
struct mem_allocator;
@@ -577,6 +578,33 @@ ssp_ran_sphere_hg
return f33_mulf3(sample, f33_basis(basis, up), sample_local);
}
+/*******************************************************************************
+ * Gaussian distribution
+ ******************************************************************************/
+
+SSP_API res_T
+ssp_distribution_gaussian_init
+(struct mem_allocator* allocator,
+ double mu,
+ double sigma,
+ struct ssp_distrib_gaussian **out_dist);
+
+SSP_API res_T
+ssp_distribution_gaussian_ref_get
+(struct ssp_distrib_gaussian* distrib);
+
+SSP_API res_T
+ssp_distribution_gaussian_ref_put
+(struct ssp_distrib_gaussian* distrib);
+
+SSP_API double
+ssp_distribution_gaussian_get
+(const struct ssp_distrib_gaussian *dist, struct ssp_rng* rng);
+
+SSP_API double
+ssp_distribution_gaussian_pdf
+(const struct ssp_distrib_gaussian *dist, double x);
+
END_DECLS
#endif /* SSP_H */
diff --git a/src/ssp_distributions.cpp b/src/ssp_distributions.cpp
@@ -0,0 +1,118 @@
+#include "ssp_distributions.h"
+
+#include <rsys/mem_allocator.h>
+#include <rsys/math.h>
+#include "ssp_rng_c.h"
+
+using normal_dist_double=RAN_NAMESPACE::normal_distribution<double>;
+
+struct ssp_distrib_gaussian_state
+{
+ normal_dist_double *distribution;
+ double mu;
+ double K1; // 1.0 / sigma;
+ double K2; // 1.0 / (sigma * SQRT_2_PI)
+};
+
+static void
+distrib_gaussian_release(ref_T* ref)
+{
+ ssp_distrib_gaussian* distrib;
+ 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);
+ }
+ MEM_RM(distrib->allocator, distrib);
+}
+
+/*******************************************************************************
+ * 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_distrib_gaussian *dist = 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) {
+ 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) {
+ res = RES_MEM_ERR;
+ goto err;
+ }
+ tmp = MEM_ALLOC(allocator, sizeof(normal_dist_double));
+ 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;
+ end:
+ return res;
+ err:
+ if (dist) {
+ ref_put(&dist->ref, distrib_gaussian_release);
+ dist = nullptr;
+ }
+ goto end;
+}
+
+res_T
+ssp_distribution_gaussian_ref_get
+(struct ssp_distrib_gaussian* distrib)
+{
+ if(!distrib) return RES_BAD_ARG;
+ ref_get(&distrib->ref);
+ return RES_OK;
+}
+
+res_T
+ssp_distribution_gaussian_ref_put
+(struct ssp_distrib_gaussian* distrib)
+{
+ if(!distrib) return RES_BAD_ARG;
+ ref_put(&distrib->ref, distrib_gaussian_release);
+ return RES_OK;
+}
+
+double
+ssp_distribution_gaussian_get
+(const struct ssp_distrib_gaussian *dist, struct ssp_rng* rng)
+{
+ class rng_cxx r(*rng);
+ return dist->state->distribution->operator()(r);
+}
+
+double
+ssp_distribution_gaussian_pdf
+(const struct ssp_distrib_gaussian *dist, double x)
+{
+ const double tmp = (x - dist->state->mu) * dist->state->K1;
+ return dist->state->K2 * exp(-0.5 * tmp * tmp);
+}
diff --git a/src/ssp_distributions.h b/src/ssp_distributions.h
@@ -0,0 +1,22 @@
+#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 ssp_distrib_gaussian_state;
+
+/* Forward declaration of external types */
+struct ssp_rng;
+struct mem_allocator;
+
+struct ssp_distrib_gaussian {
+ size_t sizeof_params;
+ size_t sizeof_state;
+ ref_T ref;
+ struct mem_allocator* allocator;
+ struct ssp_distrib_gaussian_state *state;
+};
+
+#endif
diff --git a/src/ssp_rng.c b/src/ssp_rng.c
@@ -34,25 +34,6 @@
#include <rsys/dynamic_array_char.h>
#include <rsys/mem_allocator.h>
-#ifdef COMPILER_CL
- /* Disable some warnings in boost includes */
- #pragma warning(push)
- #pragma warning(disable:4244) /* possible loss of data due to data conversion */
- #pragma warning(disable:4458) /* declaration of a variable hides class member */
-
- /* The random C++11 library is bugged on MSVC 12 & 14. Use the boost version */
- #include <boost/random.hpp>
- #include <boost/random/random_device.hpp>
- #define RAN_NAMESPACE boost::random
-
- #pragma warning(pop)
- #pragma warning(push)
- #pragma warning(disable:4706) /* Assignment in conditional expression */
-#else
- #include <random>
- #define RAN_NAMESPACE std
-#endif
-
/* disable some warnings in Random123 includes */
#ifdef COMPILER_GCC
#pragma GCC diagnostic push
@@ -81,22 +62,6 @@
#include <cstring>
#include <limits>
-#define SQRT_2_PI 2.50662827463100050240
-
-class rng_cxx /* Wrap a SSP random generator into a CXX object */
-{
-public:
- FINLINE rng_cxx(struct ssp_rng& rng) : rng(&rng) {}
- FINLINE uint64_t operator()() { return ssp_rng_get(rng); }
- FINLINE uint64_t min() const { return ssp_rng_min(rng); }
- FINLINE uint64_t max() const { return ssp_rng_max(rng); }
-
- typedef uint64_t result_type;
-
-private:
- ssp_rng* rng;
-};
-
/*******************************************************************************
* KISS PRNG
******************************************************************************/
diff --git a/src/ssp_rng_c.h b/src/ssp_rng_c.h
@@ -35,6 +35,27 @@
#include "ssp.h"
#include <rsys/ref_count.h>
+#define SQRT_2_PI 2.50662827463100050240
+
+#ifdef COMPILER_CL
+ /* Disable some warnings in boost includes */
+ #pragma warning(push)
+ #pragma warning(disable:4244) /* possible loss of data due to data conversion */
+ #pragma warning(disable:4458) /* declaration of a variable hides class member */
+
+ /* The random C++11 library is bugged on MSVC 12 & 14. Use the boost version */
+ #include <boost/random.hpp>
+ #include <boost/random/random_device.hpp>
+ #define RAN_NAMESPACE boost::random
+
+ #pragma warning(pop)
+ #pragma warning(push)
+ #pragma warning(disable:4706) /* Assignment in conditional expression */
+#else
+ #include <random>
+ #define RAN_NAMESPACE std
+#endif
+
struct ssp_rng {
struct ssp_rng_type type;
void* state;
@@ -42,5 +63,19 @@ struct ssp_rng {
ref_T ref;
};
+class rng_cxx /* Wrap a SSP random generator into a CXX object */
+{
+public:
+ FINLINE rng_cxx(struct ssp_rng& _rng) : rng(_rng) {}
+ FINLINE uint64_t operator()() { return rng.type.get(rng.state); }
+ FINLINE uint64_t min() const { return rng.type.min; }
+ FINLINE uint64_t max() const { return rng.type.max; }
+
+ typedef uint64_t result_type;
+
+private:
+ ssp_rng& rng;
+};
+
#endif /* SSP_RNG_C_H */
diff --git a/src/test_ssp_distrib_gaussian.c b/src/test_ssp_distrib_gaussian.c
@@ -0,0 +1,106 @@
+/* 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;
+}
