schiff

Estimate the radiative properties of soft particless
git clone git://git.meso-star.com/schiff.git
Log | Files | Refs | README | LICENSE

commit ba84ffade367374cc94abaad6409096e6abf733e
parent 348cff3810baae116c46565b33881ded31c87ad2
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Mon, 23 Nov 2015 19:30:59 +0100

Update the histogram file format

Diffstat:
Msrc/schiff_args.c | 14++++++--------
Msrc/schiff_histogram.c | 202+++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------
2 files changed, 146 insertions(+), 70 deletions(-)

diff --git a/src/schiff_args.c b/src/schiff_args.c @@ -77,14 +77,12 @@ print_help(const char* binary) printf( " c=VAL parameter is the constant VAL.\n"); printf( -" 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 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"); +" h=HISTOGRAM The first valid line of the HISTOGRAM file should be formatted\n" +" as \"LOWER UPPER NINTERVALS\" where \"LOWER\" and \"UPPER\" are,\n" +" respectively, the lower and the upper bounds of the histogram,\n" +" and \"NINTERVALS\" the number of intervals. The remaining valid\n" +" lines should list the probability of the interval upper bound\n" +" until all intervals are defined.\n"); printf( " l=ZETA:SIGMA parameter is distributed with respect to the lognormal\n" " distribution:\n" diff --git a/src/schiff_histogram.c b/src/schiff_histogram.c @@ -30,6 +30,7 @@ #include "schiff_streamline.h" #include <rsys/algorithm.h> +#include <rsys/str.h> #include <rsys/cstr.h> #include <rsys/stretchy_array.h> @@ -37,18 +38,17 @@ * Helper functions ******************************************************************************/ static res_T -parse_schiff_hentry - (struct schiff_hentry* entry, +parse_histogram_proba + (const char* str, const char* filename, const size_t iline, - const char* str) + double* proba) { char buf[128]; - double tmp[2]; + double d; char* tk; - int i; res_T res = RES_OK; - ASSERT(entry && filename && str); + ASSERT(proba && filename && str); if(strlen(str) >= sizeof(buf) - 1/*NULL char*/) { fprintf(stderr, @@ -58,33 +58,20 @@ parse_schiff_hentry } strncpy(buf, str, sizeof(buf)); - for(i=0, tk=strtok(buf, " \t"); tk && i < 2; ++i, tk=strtok(NULL, " \t")) { - res = cstr_to_double(tk, tmp + i); - if(res != RES_OK) { - fprintf(stderr, "%s:%lu: cannot read the histogram value `%s'.\n", - filename, (unsigned long)iline, tk); - return res; - } - } - if(i < 2) { - fprintf(stderr, - "%s:%lu: wrong histogram value `%s'. Expect <VALUE> <PROBABILITY>\n", - filename, (unsigned long)iline, str); - return RES_BAD_ARG; + tk = strtok(buf, "#"); + res = cstr_to_double(tk, &d); + if(res != RES_OK) { + fprintf(stderr, "%s:%lu: cannot read the histogram probability `%s'.\n", + filename, (unsigned long)iline, tk); + return res; } - if(tk) { - fprintf(stderr, "%s:%lu: wrong histogram entry `%s'.\n", - filename, (unsigned long)iline, str); - return RES_BAD_ARG; - } - entry->value = tmp[0]; - 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); + if(d < 0.0) { + fprintf(stderr, "%s:%lu: wrong probability `%g'.\n", + filename, (unsigned long)iline, d); return RES_BAD_ARG; } + *proba = d; return RES_OK; } @@ -98,6 +85,67 @@ cmp_entry(const void* a, const void* b) return val < entry->accum_proba ? -1 : val > entry->accum_proba ? 1 : 0; } +static res_T +parse_histogram_header + (const char* str, + const char* filename, + const size_t iline, + double bounds[2], + unsigned long* size) +{ + struct str buf; + char* tk; + long l; + res_T res = RES_OK; + ASSERT(str && filename && bounds && size); + + str_init(&mem_default_allocator, &buf); + res = str_set(&buf, str); + if(res != RES_OK) { + fprintf(stderr, + "Not enough memory. Couldn't parse the histogram header `%s'.\n", str); + goto error; + } + + tk = strtok(str_get(&buf), " \t"); + res = cstr_to_double(tk, bounds + 0); + if(res != RES_OK) { + fprintf(stderr, "%s:%lu: invalid histogram lower bound `%s'.\n", + filename, (unsigned long)iline, tk); + goto error; + } + + tk = strtok(NULL, " \t"); + res = cstr_to_double(tk, bounds + 1); + if(res != RES_OK) { + fprintf(stderr, "%s:%lu: invalid histogram upper bound `%s'.\n", + filename, (unsigned long)iline, tk); + goto error; + } + + if(bounds[0] > bounds[1]) { + fprintf(stderr, "%s:%lu: invalid histogram bounds `[%g, %g]'.\n", + filename, (unsigned long)iline, bounds[0], bounds[1]); + res = RES_BAD_ARG; + goto error; + } + + tk = strtok(NULL, "#"); + res = cstr_to_long(tk, &l); + if(res != RES_OK || l <= 0) { + fprintf(stderr, "%s:%lu: invalid number of histogram intervals `%s'.\n", + filename, (unsigned long)iline, tk); + if(res == RES_OK) res = RES_BAD_ARG; + } + *size = (size_t)l; + +exit: + str_release(&buf); + return res; +error: + goto exit; +} + /******************************************************************************* * Local functions ******************************************************************************/ @@ -134,8 +182,13 @@ schiff_histogram_load_stream { struct schiff_streamline streamline; size_t iline; + size_t i; char* line; - double accum_proba = 0.0; + double bounds[2]; + double step; + double val; + double accum_proba; + size_t nentries; res_T res = RES_OK; ASSERT(histo && stream_name && stream_name); @@ -143,47 +196,72 @@ schiff_histogram_load_stream schiff_streamline_init(&streamline); iline = 1; - while((res = schiff_streamline_read(&streamline, stream, &line)) != RES_EOF) { - struct schiff_hentry entry; - if(*line == '\0' /*Empty line*/ || *line == '#' /* Comment */) - continue; + /* Skip empty line and comments */ + while((res = schiff_streamline_read(&streamline, stream, &line)) != RES_EOF + && strcspn(line, "# \t") == 0) { + ++iline; + } + if(res == RES_EOF) { + fprintf(stderr, "%s: invalid empty histogram.\n", stream_name); + goto error; + } - res = parse_schiff_hentry(&entry, stream_name, iline, strtok(line, "#")); + res = parse_histogram_header + (strtok(line, "#"), stream_name, iline, bounds, &nentries); + if(res != RES_OK) goto error; + ++iline; + + res = darray_hentry_resize(histo, (size_t)nentries); + if(res != RES_OK) { + fprintf(stderr, "%s:%lu: couldn't allocate the histogram entries.\n", + stream_name, (unsigned long)iline); + goto error; + } + + step = (bounds[1] - bounds[0])/(double)nentries; + val = bounds[0]; + accum_proba = 0; + for(;(res = schiff_streamline_read(&streamline, stream, &line))!=RES_EOF + && i<nentries; ++iline) { + double proba; + + if(strcspn(line, "# \t") == 0) continue; /* Empty line */ + res = parse_histogram_proba(strtok(line, "#"), stream_name, iline, &proba); 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; + accum_proba += proba; + darray_hentry_data_get(histo)[i].value = val; + darray_hentry_data_get(histo)[i].accum_proba = accum_proba; + ++i; + val += step; } - if(res != RES_EOF) goto error; - res = RES_OK; - if(!eq_eps(accum_proba, 1, 1.e-6)) { - fprintf(stderr, "%s: unormalized histogram.\n", stream_name); + + if(i < nentries) { + fprintf(stderr, +"%s: missing probabilities. Expecting %lu values while %lu is submitted.\n", + stream_name, (unsigned long)nentries, (unsigned long)i); res = RES_BAD_ARG; goto error; } + if(res == RES_OK) { + fprintf(stderr, "%s:%lu ignore subsequent lines.\n", + stream_name, (unsigned long)iline); + } else if(res == RES_EOF) { + res = RES_OK; + } + if(res != RES_OK) goto error; + + /* Normalize the histogram */ + FOR_EACH(i, 0, darray_hentry_size_get(histo) - 1) + darray_hentry_data_get(histo)[i].accum_proba /= accum_proba; + + /* handle numerical issue on the last entry */ + i = darray_hentry_size_get(histo)-1; + darray_hentry_data_get(histo)[i].value = bounds[1]; + darray_hentry_data_get(histo)[i].accum_proba = 1.0; + exit: schiff_streamline_release(&streamline); return res;