commit e875fb97b81112c7279344947a232466e8944d15
parent f187fa13f89458d34d7675d91250c9c04daf6bcf
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 30 Jul 2018 13:58:03 +0200
Implement and test htgop_spectral_interval_sample_quadrature
Diffstat:
5 files changed, 218 insertions(+), 130 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -98,10 +98,10 @@ if(NOT NO_TEST)
new_test(test_htgop)
build_test(test_htgop_load)
- build_test(test_htgop_sw_sample)
+ build_test(test_htgop_sample)
add_test(test_htgop_load test_htgop_load ${_etc_dst}/ecrad_opt_prop.txt)
- add_test(test_htgop_sw_sample test_htgop_sw_sample ${_etc_dst}/ecrad_opt_prop.txt)
+ add_test(test_htgop_sample test_htgop_sample ${_etc_dst}/ecrad_opt_prop.txt)
endif()
################################################################################
diff --git a/src/htgop.c b/src/htgop.c
@@ -13,6 +13,8 @@
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>. */
+#define _POSIX_C_SOURCE 200112L /* nextafter support */
+
#include "htgop.h"
#include "htgop_c.h"
#include "htgop_reader.h"
@@ -239,6 +241,7 @@ sample_sw_spectral_interval_CIE
const double r, /* Canonical number in [0, 1[ */
size_t* id)
{
+ double r_next = nextafter(r, DBL_MAX);
double* find;
size_t i;
ASSERT(htgop && cdf && cdf_length);
@@ -249,17 +252,18 @@ sample_sw_spectral_interval_CIE
return RES_BAD_ARG;
}
- find = search_lower_bound(&r, cdf, cdf_length, sizeof(double), cmp_dbl);
- if(!find) {
+ /* Use r_next rather than r in order to find the first entry that is not less
+ * that *or equal* to r */
+ find = search_lower_bound(&r_next, cdf, cdf_length, sizeof(double), cmp_dbl);
+ if(!find) { /* <=> CIE XYZ band does not intersect the loaded intervals */
log_err(htgop,
"%s: could not sample a spectral interval in the CIE XYZ spectral band.\n",
caller);
return RES_BAD_OP;
}
- ASSERT(find);
i = (size_t)(find - cdf);
- ASSERT(i < cdf_length && cdf[i] > r && (!i || cdf[i-1] < r));
+ ASSERT(i < cdf_length && cdf[i] > r && (!i || cdf[i-1] <= r));
*id = i;
return RES_OK;
}
@@ -813,6 +817,7 @@ htgop_get_lw_spectral_interval
interval->quadrature_pdf = darray_double_cdata_get(&quad_pdfs[ispecint]);
interval->quadrature_cdf = darray_double_cdata_get(&quad_cdfs[ispecint]);
interval->quadrature_length = darray_double_size_get(&quad_pdfs[ispecint]);
+ interval->data__ = htgop;
return RES_OK;
}
@@ -841,6 +846,40 @@ htgop_get_sw_spectral_interval
interval->quadrature_pdf = darray_double_cdata_get(&quad_pdfs[ispecint]);
interval->quadrature_cdf = darray_double_cdata_get(&quad_cdfs[ispecint]);
interval->quadrature_length = darray_double_size_get(&quad_pdfs[ispecint]);
+ interval->data__ = htgop;
+ return RES_OK;
+}
+
+res_T
+htgop_spectral_interval_sample_quadrature
+ (const struct htgop_spectral_interval* specint,
+ const double r, /* Canonical number in [0, 1[ */
+ size_t* iquad_point) /* Id of the sample quadrature point */
+{
+ const struct htgop* htgop;
+ double r_next = nextafter(r, DBL_MAX);
+ double* find;
+ size_t i;
+
+ if(!specint || !iquad_point) return RES_BAD_ARG;
+
+ htgop = specint->data__;
+ if(r < 0 || r >= 1) {
+ log_err(htgop, "%s: invalid canonical random number `%g'.\n", FUNC_NAME, r);
+ return RES_BAD_ARG;
+ }
+
+
+ /* Use r_next rather than r in order to find the first entry that is not less
+ * that *or equal* to r */
+ find = search_lower_bound(&r_next, specint->quadrature_cdf,
+ specint->quadrature_length, sizeof(double), cmp_dbl);
+ ASSERT(find);
+ i = (size_t)(find - specint->quadrature_cdf);
+ ASSERT(i < specint->quadrature_length);
+ ASSERT(specint->quadrature_cdf[i] > r);
+ ASSERT(!i || specint->quadrature_cdf[i-1] <= r);
+ *iquad_point = i;
return RES_OK;
}
diff --git a/src/htgop.h b/src/htgop.h
@@ -54,6 +54,7 @@ struct htgop_spectral_interval {
const double* quadrature_pdf; /* PDF of the quadrature */
const double* quadrature_cdf; /* CDF of the quadrature */
size_t quadrature_length; /* #quadrature points of the interval */
+ const void* data__; /* Internal data */
};
struct htgop_level {
@@ -203,11 +204,10 @@ htgop_get_sw_spectral_interval
struct htgop_spectral_interval* interval);
HTGOP_API res_T
-htgop_sw_spectral_interval_sample_quadrature
- (const struct htgop* htgop,
- const size_t ispectral_interval,
+htgop_spectral_interval_sample_quadrature
+ (const struct htgop_spectral_interval* interval,
const double r, /* Canonical random number in [0, 1[ */
- size_t* iquad); /* Id of the sampled interval */
+ size_t* iquad_point); /* Id of the sample quadrature point */
HTGOP_API res_T
htgop_layer_get_lw_spectral_interval
diff --git a/src/test_htgop_sample.c b/src/test_htgop_sample.c
@@ -0,0 +1,169 @@
+/* Copyright (C) |Meso|Star> 2018 (contact@meso-star.com)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "htgop.h"
+#include "test_htgop_utils.h"
+
+#include <rsys/math.h>
+#include <string.h>
+
+#define N 100000
+#define CIE_XYZ_WAVENUMBER_MIN (1.0e7/780.0)
+#define CIE_XYZ_WAVENUMBER_MAX (1.0e7/380.0)
+
+static void
+check_sw_sample_specint
+ (struct htgop* htgop,
+ res_T (*sample_func)(const struct htgop*, const double, size_t*),
+ int* hist)
+{
+ const double* wnums;
+ size_t nspecints;
+ size_t ispecint;
+ size_t i;
+
+ CHK(sample_func && htgop && hist);
+ CHK(htgop_get_sw_spectral_intervals_count(htgop, &nspecints) == RES_OK);
+ CHK(htgop_get_sw_spectral_intervals_wave_numbers(htgop, &wnums) == RES_OK);
+
+ CHK(sample_func(NULL, 0, &ispecint) == RES_BAD_ARG);
+ CHK(sample_func(htgop, 1, &ispecint) == RES_BAD_ARG);
+ CHK(sample_func(htgop, 0, NULL) == RES_BAD_ARG);
+
+ memset(hist, 0, sizeof(*hist)*nspecints);
+ FOR_EACH(i, 0, N) {
+ const double r = rand_canonic();
+ CHK(sample_func(htgop, r, &ispecint) == RES_OK);
+ CHK(ispecint < nspecints);
+ CHK(wnums[ispecint+0] <= CIE_XYZ_WAVENUMBER_MAX);
+ CHK(wnums[ispecint+1] >= CIE_XYZ_WAVENUMBER_MIN);
+ hist[ispecint] += 1;
+ }
+}
+
+static void
+check_sample_quadrature
+ (struct htgop* htgop,
+ res_T (*get_nspecints)(const struct htgop*, size_t*),
+ res_T (*get_specint)
+ (const struct htgop*, const size_t, struct htgop_spectral_interval*))
+{
+ struct htgop_spectral_interval specint;
+ size_t nspecints;
+ size_t iquadpt;
+ size_t i;
+
+ CHK(htgop && get_specint);
+
+ CHK(get_nspecints(htgop, &nspecints) == RES_OK);
+ CHK(nspecints);
+
+ CHK(get_specint(htgop, 0, &specint) == RES_OK);
+ CHK(htgop_spectral_interval_sample_quadrature(NULL, 0, &iquadpt) == RES_BAD_ARG);
+ CHK(htgop_spectral_interval_sample_quadrature(&specint, 1, &iquadpt) == RES_BAD_ARG);
+ CHK(htgop_spectral_interval_sample_quadrature(&specint, 1, NULL) == RES_BAD_ARG);
+ CHK(htgop_spectral_interval_sample_quadrature(&specint, 1, NULL) == RES_BAD_ARG);
+
+ FOR_EACH(i, 0, 10) {
+ int* hist;
+ size_t ispecint;
+ size_t j;
+
+ ispecint = (size_t)(rand_canonic() * (double)nspecints);
+ CHK(get_specint(htgop, ispecint, &specint) == RES_OK);
+
+ CHK(hist = mem_calloc(specint.quadrature_length, sizeof(*hist)));
+ FOR_EACH(j, 0, N) {
+ const double r = rand_canonic();
+ CHK(htgop_spectral_interval_sample_quadrature(&specint, r, &iquadpt) == RES_OK);
+ CHK(iquadpt < specint.quadrature_length);
+ CHK(specint.quadrature_cdf[iquadpt] > r);
+ CHK(!iquadpt || specint.quadrature_cdf[iquadpt-1] <= r);
+
+ hist[iquadpt] += 1;
+ }
+
+ FOR_EACH(iquadpt, 0, specint.quadrature_length) {
+ if(!eq_eps(specint.quadrature_pdf[iquadpt], 0, 1.e-4))
+ CHK(hist[iquadpt] > 0);
+ }
+
+ mem_rm(hist);
+ }
+}
+
+
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct htgop* htgop;
+ const double* wnums;
+ int* hist_X;
+ int* hist_Y;
+ int* hist_Z;
+ size_t nspecints;
+ size_t ispecint;
+
+ if(argc < 2) {
+ fprintf(stderr, "Usage: %s FILENAME\n", argv[0]);
+ return 1;
+ }
+
+ CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
+
+ CHK(htgop_create(NULL, &allocator, 1, &htgop) == RES_OK);
+ CHK(htgop_load(htgop, argv[1]) == RES_OK);
+ CHK(htgop_get_sw_spectral_intervals_count(htgop, &nspecints) == RES_OK);
+ CHK(nspecints > 0);
+
+ CHK(htgop_get_sw_spectral_intervals_wave_numbers(htgop, &wnums) == RES_OK);
+
+ CHK(hist_X = MEM_CALLOC(&allocator, nspecints, sizeof(*hist_X)));
+ CHK(hist_Y = MEM_CALLOC(&allocator, nspecints, sizeof(*hist_Y)));
+ CHK(hist_Z = MEM_CALLOC(&allocator, nspecints, sizeof(*hist_Z)));
+
+ check_sw_sample_specint(htgop, htgop_sample_sw_spectral_interval_CIE_X, hist_X);
+ check_sw_sample_specint(htgop, htgop_sample_sw_spectral_interval_CIE_Y, hist_Y);
+ check_sw_sample_specint(htgop, htgop_sample_sw_spectral_interval_CIE_Z, hist_Z);
+
+ FOR_EACH(ispecint, 0, nspecints) {
+ if(wnums[ispecint+0] > CIE_XYZ_WAVENUMBER_MAX
+ || wnums[ispecint+1] < CIE_XYZ_WAVENUMBER_MIN) {
+ CHK(hist_X[ispecint] == 0);
+ CHK(hist_Y[ispecint] == 0);
+ CHK(hist_Z[ispecint] == 0);
+ } else {
+ CHK(!hist_X[ispecint] || hist_X[ispecint] != hist_Y[ispecint]);
+ CHK(!hist_X[ispecint] || hist_X[ispecint] != hist_Z[ispecint]);
+ CHK(!hist_Y[ispecint] || hist_Y[ispecint] != hist_Z[ispecint]);
+ }
+ }
+
+ MEM_RM(&allocator, hist_X);
+ MEM_RM(&allocator, hist_Y);
+ MEM_RM(&allocator, hist_Z);
+
+ check_sample_quadrature(htgop, htgop_get_sw_spectral_intervals_count,
+ htgop_get_sw_spectral_interval);
+ check_sample_quadrature(htgop, htgop_get_lw_spectral_intervals_count,
+ htgop_get_lw_spectral_interval);
+
+ CHK(htgop_ref_put(htgop) == RES_OK);
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
diff --git a/src/test_htgop_sw_sample.c b/src/test_htgop_sw_sample.c
@@ -1,120 +0,0 @@
-/* Copyright (C) |Meso|Star> 2018 (contact@meso-star.com)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#include "htgop.h"
-#include "test_htgop_utils.h"
-
-#include <string.h>
-
-#define N 100000
-#define CIE_XYZ_WAVENUMBER_MIN (1.0e7/780.0)
-#define CIE_XYZ_WAVENUMBER_MAX (1.0e7/380.0)
-
-static void
-run_sample_func
- (struct htgop* htgop,
- res_T (*sample_func)(const struct htgop*, const double, size_t*),
- int* hist)
-{
- const double* wnums;
- size_t nspecints;
- size_t ispecint;
- size_t i;
-
- CHK(sample_func && htgop && hist);
- CHK(htgop_get_sw_spectral_intervals_count(htgop, &nspecints) == RES_OK);
- CHK(htgop_get_sw_spectral_intervals_wave_numbers(htgop, &wnums) == RES_OK);
-
- memset(hist, 0, sizeof(*hist)*nspecints);
- FOR_EACH(i, 0, N) {
- const double r = rand_canonic();
- CHK(sample_func(htgop, r, &ispecint) == RES_OK);
- CHK(ispecint < nspecints);
- CHK(wnums[ispecint+0] <= CIE_XYZ_WAVENUMBER_MAX);
- CHK(wnums[ispecint+1] >= CIE_XYZ_WAVENUMBER_MIN);
- hist[ispecint] += 1;
- }
-}
-
-int
-main(int argc, char** argv)
-{
- struct mem_allocator allocator;
- struct htgop* htgop;
- const double* wnums;
- int* hist_X;
- int* hist_Y;
- int* hist_Z;
- size_t nspecints;
- size_t ispecint;
-
- if(argc < 2) {
- fprintf(stderr, "Usage: %s FILENAME\n", argv[0]);
- return 1;
- }
-
- CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK);
-
- CHK(htgop_create(NULL, &allocator, 1, &htgop) == RES_OK);
- CHK(htgop_load(htgop, argv[1]) == RES_OK);
- CHK(htgop_get_sw_spectral_intervals_count(htgop, &nspecints) == RES_OK);
- CHK(nspecints > 0);
-
- CHK(htgop_get_sw_spectral_intervals_wave_numbers(htgop, &wnums) == RES_OK);
-
- hist_X = MEM_CALLOC(&allocator, nspecints, sizeof(*hist_X));
- hist_Y = MEM_CALLOC(&allocator, nspecints, sizeof(*hist_Y));
- hist_Z = MEM_CALLOC(&allocator, nspecints, sizeof(*hist_Z));
- CHK(hist_X && hist_Y && hist_Z);
-
- CHK(htgop_sample_sw_spectral_interval_CIE_X(NULL, 0, &ispecint) == RES_BAD_ARG);
- CHK(htgop_sample_sw_spectral_interval_CIE_X(htgop, 1, &ispecint) == RES_BAD_ARG);
- CHK(htgop_sample_sw_spectral_interval_CIE_X(htgop, 0, NULL) == RES_BAD_ARG);
-
- CHK(htgop_sample_sw_spectral_interval_CIE_Y(NULL, 0, &ispecint) == RES_BAD_ARG);
- CHK(htgop_sample_sw_spectral_interval_CIE_Y(htgop, 1, &ispecint) == RES_BAD_ARG);
- CHK(htgop_sample_sw_spectral_interval_CIE_Y(htgop, 0, NULL) == RES_BAD_ARG);
-
- CHK(htgop_sample_sw_spectral_interval_CIE_Z(NULL, 0, &ispecint) == RES_BAD_ARG);
- CHK(htgop_sample_sw_spectral_interval_CIE_Z(htgop, 1, &ispecint) == RES_BAD_ARG);
- CHK(htgop_sample_sw_spectral_interval_CIE_Z(htgop, 0, NULL) == RES_BAD_ARG);
-
- run_sample_func(htgop, htgop_sample_sw_spectral_interval_CIE_X, hist_X);
- run_sample_func(htgop, htgop_sample_sw_spectral_interval_CIE_Y, hist_Y);
- run_sample_func(htgop, htgop_sample_sw_spectral_interval_CIE_Z, hist_Z);
-
- FOR_EACH(ispecint, 0, nspecints) {
- if(wnums[ispecint+0] > CIE_XYZ_WAVENUMBER_MAX
- || wnums[ispecint+1] < CIE_XYZ_WAVENUMBER_MIN) {
- CHK(hist_X[ispecint] == 0);
- CHK(hist_Y[ispecint] == 0);
- CHK(hist_Z[ispecint] == 0);
- } else {
- CHK(!hist_X[ispecint] || hist_X[ispecint] != hist_Y[ispecint]);
- CHK(!hist_X[ispecint] || hist_X[ispecint] != hist_Z[ispecint]);
- CHK(!hist_Y[ispecint] || hist_Y[ispecint] != hist_Z[ispecint]);
- }
- }
-
- MEM_RM(&allocator, hist_X);
- MEM_RM(&allocator, hist_Y);
- MEM_RM(&allocator, hist_Z);
-
- CHK(htgop_ref_put(htgop) == RES_OK);
- check_memory_allocator(&allocator);
- mem_shutdown_proxy_allocator(&allocator);
- CHK(mem_allocated_size() == 0);
- return 0;
-}
