commit fb337ef6236d964991f8fd43c274a56f7405385f
parent b6cb3e32af3963739a363c265f196fda4f494748
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Thu, 7 Jul 2016 10:15:47 +0200
Add uniform disk distribution.
Diffstat:
3 files changed, 116 insertions(+), 0 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -156,6 +156,7 @@ if(NOT NO_TEST)
new_test(test_ssp_rng_proxy)
new_test(test_ssp_ran_gaussian)
new_test(test_ssp_ran_piecewise_linear)
+ new_test(test_ssp_ran_uniform_disk)
endif()
################################################################################
diff --git a/src/ssp.h b/src/ssp.h
@@ -643,6 +643,26 @@ ssp_ran_piecewise_linear_get
struct ssp_rng* rng,
double *res);
+/*******************************************************************************
+ * Uniform disk distribution
+ ******************************************************************************/
+
+static FINLINE res_T
+ssp_ran_uniform_disk
+ (struct ssp_rng* rng,
+ double radius,
+ double pt[2])
+{
+ double theta, r;
+ if(!rng || !pt || radius <= 0) return RES_BAD_ARG;
+
+ theta = ssp_rng_uniform_double(rng, 0, 2 * PI);
+ r = radius * sqrt(ssp_rng_canonical(rng));
+ pt[0] = r * cos(theta);
+ pt[1] = r * sin(theta);
+ return RES_OK;
+}
+
END_DECLS
#endif /* SSP_H */
diff --git a/src/test_ssp_ran_uniform_disk.c b/src/test_ssp_ran_uniform_disk.c
@@ -0,0 +1,95 @@
+/* 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/math.h>
+
+#define NBS 1000000
+
+int
+main(int argc, char** argv)
+{
+ struct ssp_rng_proxy* proxy;
+ struct ssp_rng* rng;
+ struct mem_allocator allocator;
+ int i;
+ unsigned r, c, nb = 0;
+ double pt[2];
+ double x = 0, x2 = 0;
+ double mean, std;
+ double exp_mean = 20 * 20 * NBS / (PI * 100 * 100 ), exp_std = 0;
+ unsigned counts[10][10] = {0};
+ (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_uniform_disk(NULL, 100, pt), RES_BAD_ARG);
+ CHECK(ssp_ran_uniform_disk(rng, -100, pt), RES_BAD_ARG);
+ CHECK(ssp_ran_uniform_disk(rng, 100, pt), RES_OK);
+
+ for (i = 0; i < NBS; i++) {
+ CHECK(ssp_ran_uniform_disk(rng, 100, pt), RES_OK);
+ /* locate pt in a 10x10 grid */
+ r = (unsigned)((100 + pt[0]) / 20);
+ c = (unsigned)((100 + pt[1]) / 20);
+ ++counts[r][c];
+ }
+
+ for (r = 0; r < 10; r++) {
+ unsigned _x = (r >= 5 ? r - 4 : r - 5) * 20;
+ for (c = 0; c < 10; c++) {
+ unsigned _y = (c >= 5 ? c - 4 : c - 5) * 20;
+ unsigned r2 = _x * _x + _y * _y;
+ if (r2 > 100 * 100)
+ /* this square is not (fully) in the disk */
+ continue;
+ ++nb;
+ x += counts[r][c];
+ x2 += counts[r][c] * counts[r][c];
+ }
+ }
+ mean = x/nb;
+ std = sqrt(x2/nb - x/nb*x/nb);
+ printf("%g %g\n", mean, std);
+ CHECK(fabs(mean - exp_mean) < 10, 1);
+ CHECK(fabs(std - exp_std) < 100, 1);
+
+ 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;
+}
