commit 8dc79ed97b2cb6fb26a5916c8b16a95455fe58f1
parent ac25ce9e4c2a482bad8eec78acd616263d27490f
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 23 Nov 2015 12:26:33 +0100
Update the histogram distribution
Replace the histogram of discreet probabilities by an histogram of
probabilities to sample an interval. Values inside the interval are
retrieved by a linear interpolation of interval boundaries.
Diffstat:
3 files changed, 63 insertions(+), 26 deletions(-)
diff --git a/src/schiff_args.c b/src/schiff_args.c
@@ -79,7 +79,11 @@ print_help(const char* binary)
" h=HISTOGRAM parameter follows the distribution whose the normalized\n"
" histogram is described by the HISTOGRAM file. Each line of the\n"
" HISTOGRAM file must be formatted as \"VAL PROBA\", where PROBA is\n"
-" the probability of the VAL value.\n");
+" the probability to sample the interval whose upper bound is VAL\n"
+" and lower bound is the VAL defined by the previous line.\n");
+ printf(
+" Consquently the lines must be sorted in ascending order\n"
+" according to VAL while the PROBA of the first line must be 0.\n");
printf(
" l=ZETA:SIGMA parameter is distributed with respect to the lognormal\n"
" distribution:\n"
@@ -96,14 +100,18 @@ print_help(const char* binary)
SCHIFF_CYLINDER_DEFAULT.nslices);
printf(
" -C ASPECT_RATIO,RDISTRIB[,N]\n"
-" the soft particles are cylindrical meshes discretized in N\n"
-" slices. By default N is %u. The ratio between their height\n"
-" and their radius is fixed by the ASPECT_RATIO floating point\n"
-" parameter. Their volume is defined according to the volume of\n"
-" a sphere whose radius is distributed with respect to the\n"
-" RDISTRIB distribution.\n",
+" the soft particles are cylindrical meshes discretized in N\n"
+" slices. By default N is %u. The volume of the cylinders is\n"
+" defined according to the volume of a sphere whose radius \"R\" is\n"
+" distributed with respect to the RDISTRIB distribution. The\n"
+" ratio between the cylinder radius \"r\" and its height \"h\" is\n"
+" fixed by the ASPECT_RATIO floating point parameter.\n",
SCHIFF_CYLINDER_DEFAULT.nslices);
printf(
+" r = R * ((2*ASPECT_RATIO)/3)^1/3\n"
+" h = r * 2 / ASPECT_RATIO\n");
+
+ printf(
" -d DIRS number of sampled directions for each geometry. Default is %u.\n",
SCHIFF_ARGS_NULL.ndirs);
printf(
@@ -134,7 +142,7 @@ print_help(const char* binary)
" the soft particles are 3D super shapes discretized in N slices\n"
" along 2PI. By default N is %u. The A, B, N0 and N1 parameters\n"
" of the super formulas are fixed to 1 while the MDISTRIBG and\n"
-" the N2DISTRIB control de distribution of the M, and N2\n"
+" the N2DISTRIB control the distribution of the M, and N2\n"
" supershape parameters, respectively.\n",
SCHIFF_SUPER_SHAPE_DEFAULT.nslices);
printf(
diff --git a/src/schiff_histogram.c b/src/schiff_histogram.c
@@ -79,7 +79,12 @@ parse_schiff_hentry
return RES_BAD_ARG;
}
entry->value = tmp[0];
- entry->cdf = tmp[1];
+ entry->accum_proba = tmp[1];
+ if(entry->accum_proba < 0.0) {
+ fprintf(stderr, "%s:%lu: wrong accumulated probability `%g'.\n",
+ filename, (unsigned long)iline, entry->accum_proba);
+ return RES_BAD_ARG;
+ }
return RES_OK;
}
@@ -90,7 +95,7 @@ cmp_entry(const void* a, const void* b)
const struct schiff_hentry* entry = b;
ASSERT(a && b);
val = *(const double*)a;
- return val < entry->cdf ? -1 : val > entry->cdf ? 1 : 0;
+ return val < entry->accum_proba ? -1 : val > entry->accum_proba ? 1 : 0;
}
/*******************************************************************************
@@ -130,7 +135,7 @@ schiff_histogram_load_stream
struct schiff_streamline streamline;
size_t iline;
char* line;
- double cdf = 0.0;
+ double accum_proba = 0.0;
res_T res = RES_OK;
ASSERT(histo && stream_name && stream_name);
@@ -145,27 +150,40 @@ schiff_histogram_load_stream
continue;
res = parse_schiff_hentry(&entry, stream_name, iline, strtok(line, "#"));
- if(res == RES_OK) { /* *NO* error */
- cdf += entry.cdf;
- entry.cdf = cdf;
- res = darray_hentry_push_back(histo, &entry);
- if(res != RES_OK) {
- fprintf(stderr, "%s:%lu: couldn't store the histogram entry.\n",
- stream_name, (unsigned long)iline);
- goto error;
- }
+ if(res != RES_OK) goto error;
+
+ /* On first line the accum_proba must be 0 */
+ if(iline == 1 && entry.accum_proba != 0.0) {
+ res = RES_BAD_ARG;
+ fprintf(stderr,
+"%s:%lu: invalid entry `%g %g'. The accumulated probability of the first entry must be 0.\n",
+ stream_name, (unsigned long)iline, entry.value, entry.accum_proba);
+ goto error;
+ } else if(darray_hentry_cdata_get(histo)[iline-2].value >= entry.value) {
+ res = RES_BAD_ARG;
+ fprintf(stderr,
+"%s:%lu: invalid entry `%g %g'. Histogram value must be sorted in ascending order.\n",
+ stream_name, (unsigned long)iline, entry.value, entry.accum_proba);
+ goto error;
+ }
+ accum_proba += entry.accum_proba;
+ entry.accum_proba = accum_proba;
+ res = darray_hentry_push_back(histo, &entry);
+ if(res != RES_OK) {
+ fprintf(stderr, "%s:%lu: couldn't store the histogram entry.\n",
+ stream_name, (unsigned long)iline);
+ goto error;
}
++iline;
}
if(res != RES_EOF) goto error;
res = RES_OK;
- if(!eq_eps(cdf, 1, 1.e-6)) {
+ if(!eq_eps(accum_proba, 1, 1.e-6)) {
fprintf(stderr, "%s: unormalized histogram.\n", stream_name);
res = RES_BAD_ARG;
goto error;
}
-
exit:
schiff_streamline_release(&streamline);
return res;
@@ -180,6 +198,7 @@ schiff_histogram_sample(const struct darray_hentry* histo, const double u)
const struct schiff_hentry* entries;
const struct schiff_hentry* find;
size_t nentries;
+ double v;
ASSERT(histo && u >= 0.0 && u < 1.0);
entries = darray_hentry_cdata_get(histo);
@@ -188,10 +207,20 @@ schiff_histogram_sample(const struct darray_hentry* histo, const double u)
find = search_lower_bound
(&u, entries, nentries, sizeof(struct schiff_hentry), cmp_entry);
- if(!find) { /* Numerical issue */
- return entries[nentries - 1].value;
- } else {
+
+ if(find == entries)
return find->value;
+
+ if(!find) find = entries + (nentries - 1); /* Handle numerical issue */
+
+ /* Linear interpolation */
+ v = (u - find[-1].accum_proba) / (find[0].accum_proba - find[-1].accum_proba);
+ if(eq_eps(v, 0, 1.e-6)) {
+ return find[-1].value;
+ } else if(eq_eps(v, 1, 1.e-6)) {
+ return find[0].value;
+ } else {
+ return v * (find[0].value - find[-1].value) + find[-1].value;
}
}
diff --git a/src/schiff_histogram.h b/src/schiff_histogram.h
@@ -31,7 +31,7 @@
#include <rsys/dynamic_array.h>
-struct schiff_hentry { double value, cdf; };
+struct schiff_hentry { double value, accum_proba; };
#define DARRAY_NAME hentry
#define DARRAY_DATA struct schiff_hentry
#include <rsys/dynamic_array.h>
