schiff

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

commit 936a8b08f826d6d5f6313d6fe7de96596aa12de3
parent 75cb616b0a705844cbb55e3be5c64dcc77eb6763
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Thu, 10 Mar 2016 16:29:36 +0100

Parse the YAML sphere and parameter distributions

Diffstat:
Mcmake/CMakeLists.txt | 2--
Msrc/schiff_args.c | 748+++++++++++++++++++++++++++++++++++++++++++++----------------------------------
Msrc/schiff_args.h | 2+-
Msrc/schiff_geometry.c | 53+++++++++++++++++++++++++++++++++++++++++++++++++++--
Msrc/schiff_geometry.h | 3+--
Dsrc/schiff_histogram.c | 302------------------------------------------------------------------------------
Dsrc/schiff_histogram.h | 56--------------------------------------------------------
7 files changed, 479 insertions(+), 687 deletions(-)

diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -71,13 +71,11 @@ set(SCHIFF_FILES_SRC schiff.c schiff_args.c schiff_geometry.c - schiff_histogram.c schiff_mesh.c schiff_optical_properties.c) set(SCHIFF_FILES_INC schiff_args.h schiff_geometry.h - schiff_histogram.h schiff_mesh.h schiff_optical_properties.h schiff_streamline.h) diff --git a/src/schiff_args.c b/src/schiff_args.c @@ -180,415 +180,527 @@ print_help(const char* binary) } static res_T -parse_param_distribution(const char* str, struct schiff_param* param) +parse_wavelengths(const char* str, struct schiff_args* args) { - struct str buf; - char* tk; + size_t len; + size_t i; res_T res = RES_OK; - ASSERT(str && param); - param->distribution = SCHIFF_PARAM_NONE; + ASSERT(args && str); - str_init(&mem_default_allocator, &buf); - res = str_set(&buf, str); - if(res != RES_OK) { - fprintf(stderr, - "Not enough memory. Couldn't parse the parameter distribution `%s'.\n", - str); - goto error; - } + /* How many wavelengths are submitted */ + res = cstr_to_list_double(str, NULL, &len, 0); + if(res != RES_OK) goto error; - tk = strtok(str_get(&buf), "="); - if(!strcmp(tk, "c")) { /* Constant */ - param->distribution = SCHIFF_PARAM_CONSTANT; - res = cstr_to_double(strtok(NULL, "\0"), &param->data.constant); - } else if(!strcmp(tk, "h")) { /* Histogram */ - param->distribution = SCHIFF_PARAM_HISTOGRAM; - darray_hentry_init(&mem_default_allocator, &param->data.histogram); - res = schiff_histogram_load(&param->data.histogram, strtok(NULL, "\0")); - } else if(!strcmp(tk, "l")) { /* lognormal */ - double list[2]; - size_t len; - param->distribution = SCHIFF_PARAM_LOGNORMAL; - res = cstr_to_list_double(strtok(NULL, "\0"), list, &len, 2); - if(res == RES_OK) { - if(len != 2) { - res = RES_BAD_ARG; - } else { - param->data.lognormal.zeta = list[0]; - param->data.lognormal.sigma = list[1]; - } - } - } else { - res = RES_BAD_ARG; - } + /* Reserve the wavelengths memory space */ + sa_clear(args->wavelengths); + args->wavelengths = sa_add(args->wavelengths, len); - if(res != RES_OK) { - fprintf(stderr, "Invalid parameter distribution `%s'.\n", str); - goto error; - } + /* Read the wavelengths */ + res = cstr_to_list_double(optarg, args->wavelengths, NULL, len); + if(res != RES_OK) goto error; + /* Check the validity of read wavelengths */ + FOR_EACH(i, 0, len) { + if(args->wavelengths[i] < 0.0) { + res = RES_BAD_ARG; + goto error; + } + } exit: - str_release(&buf); return res; error: - if(param->distribution == SCHIFF_PARAM_HISTOGRAM) - darray_hentry_release(&param->data.histogram); goto exit; } +static INLINE void +log_err + (const char* filename, + const yaml_node_t* node, + const char* fmt, + ...) +{ + va_list vargs_list; + ASSERT(node && fmt); + + fprintf(stderr, "%s:%lu:%lu: ", + filename, + (unsigned long)node->start_mark.line+1, + (unsigned long)node->start_mark.column+1); + va_start(vargs_list, fmt); + vfprintf(stderr, fmt, vargs_list); + va_end(vargs_list); +} + static res_T -parse_cylinder_distribution(const char* str, struct schiff_args* args) +parse_yaml_double(const char* filename, const yaml_node_t* node, double* value) { - struct str buf; - char* tk; - char* tk_ctx; res_T res = RES_OK; - ASSERT(args && str); - - str_init(&mem_default_allocator, &buf); + ASSERT(filename && node && value); - /* The distribution was already set to a non sphere distribution. */ - if(args->geom.type != SCHIFF_NONE && args->geom.type != SCHIFF_CYLINDER) { - res = RES_BAD_ARG; - goto error; + if(node->type != YAML_SCALAR_NODE) { + log_err(filename, node, "expecting a floating point number.\n"); + return RES_BAD_ARG; } - - res = str_set(&buf, str); + res = cstr_to_double((char*)node->data.scalar.value, value); if(res != RES_OK) { - fprintf(stderr, - "Not enough memory. Couldn't parse the cylinder distribution `%s'.\n", - str); - goto error; + log_err(filename, node, "invalid floating point number `%s'.\n", + node->data.scalar.value); + return RES_BAD_ARG; } - /* Distribution of the cylinder radius */ - tk = strtok_r(str_get(&buf), ",", &tk_ctx); - res = parse_param_distribution(tk, &args->geom.data.cylinder.radius); - if(res != RES_OK) goto error; - - /* Distribution of the cylinder height */ - tk = strtok_r(NULL, ",", &tk_ctx); - res = parse_param_distribution(tk, &args->geom.data.cylinder.height); - if(res != RES_OK) goto error; - - /* Cylinder discretisation */ - tk = strtok_r(NULL, "\0", &tk_ctx); - if(!tk) { - args->geom.data.cylinder.nslices = SCHIFF_CYLINDER_DEFAULT.nslices; - } else { - res = cstr_to_uint(tk, &args->geom.data.cylinder.nslices); - if(res != RES_OK) goto error; - } - args->geom.type = SCHIFF_CYLINDER; - -exit: - str_release(&buf); - return res; -error: - goto exit; + return RES_OK; } static res_T -parse_cylinder_as_sphere_distribution(const char* str, struct schiff_args* args) +parse_yaml_param_lognormal + (const char* filename, + yaml_document_t* doc, + const yaml_node_t* lognormal, + struct schiff_param* param) { - struct str buf; - char* tk; - char* tk_ctx; + enum { + ZETA = BIT(0), + SIGMA = BIT(1) + }; + size_t nattrs; + size_t i; + int mask = 0; /* Register the parsed lognormal attributes */ res_T res = RES_OK; - ASSERT(args && str); - - str_init(&mem_default_allocator, &buf); + ASSERT(filename && doc && lognormal && param); - /* The distribution was already set to a non sphere distribution. */ - if(args->geom.type != SCHIFF_NONE && args->geom.type != SCHIFF_CYLINDER) { - res = RES_BAD_ARG; - goto error; + if(lognormal->type != YAML_MAPPING_NODE) { + log_err(filename, lognormal, + "expecting a mapping of the lognormal attributes.\n"); + return RES_BAD_ARG; } - res = str_set(&buf, str); - if(res != RES_OK) { - fprintf(stderr, - "Not enough memory. Couldn't parse the cylinder distribution `%s'.\n", - str); - goto error; - } + param->distribution = SCHIFF_PARAM_LOGNORMAL; - /* Cylinder aspect ratio */ - tk = strtok_r(str_get(&buf), ",", &tk_ctx); - res = cstr_to_double(tk, &args->geom.data.cylinder.aspect_ratio); - if(res != RES_OK) goto error; + /* Parse the log normal attributes */ + nattrs = (size_t) + (lognormal->data.mapping.pairs.top - lognormal->data.mapping.pairs.start); + FOR_EACH(i, 0, nattrs) { + yaml_node_t* key, *val; + + key=yaml_document_get_node(doc, lognormal->data.mapping.pairs.start[i].key); + val=yaml_document_get_node(doc, lognormal->data.mapping.pairs.start[i].value); + ASSERT(key->type == YAML_SCALAR_NODE); + + /* lognormal zeta attribute */ + if(!strcmp((char*)key->data.scalar.value, "zeta")) { + if(mask & ZETA) { + log_err(filename, key, + "the `zeta' lognormal attribute is already defined.\n"); + return RES_BAD_ARG; + } + mask |= ZETA; + res = parse_yaml_double(filename, val, &param->data.lognormal.zeta); + + /* lognormal sigma attribute */ + } else if(!strcmp((char*)key->data.scalar.value, "sigma")) { + if(mask & SIGMA) { + log_err(filename, key, + "the `sigma' lognormal attribute is already defined.\n"); + return RES_BAD_ARG; + } + mask |= SIGMA; + res = parse_yaml_double(filename, val, &param->data.lognormal.sigma); - /* Distribution of the radius of the equivalent sphere */ - tk = strtok_r(NULL, ",", &tk_ctx); - res = parse_param_distribution(tk, &args->geom.data.cylinder.radius); - if(res != RES_OK) goto error; + /* unknown attribute */ + } else { + log_err(filename, key, "unknown lognormal attribute `%s'.\n", + key->data.scalar.value); + return RES_BAD_ARG; + } + if(res != RES_OK) return res; + } - /* Cylinder discretisation */ - tk = strtok_r(NULL, "\0", &tk_ctx); - if(!tk) { - args->geom.data.cylinder.nslices = SCHIFF_CYLINDER_DEFAULT.nslices; - } else { - res = cstr_to_uint(tk, &args->geom.data.cylinder.nslices); - if(res != RES_OK) goto error; + /* Ensure that the lognormal attributes are all parsed */ + if(!(mask & ZETA)) { + log_err(filename, lognormal, "missing the zeta lognormal attribute.\n"); + return RES_BAD_ARG; + } else if(!(mask & SIGMA)) { + log_err(filename, lognormal, "missing the sigma lognormal attribute.\n"); + return RES_BAD_ARG; } - args->geom.type = SCHIFF_CYLINDER_AS_SPHERE; -exit: - str_release(&buf); - return res; -error: - goto exit; + return RES_OK; } - static res_T -parse_sphere_distribution(const char* str, struct schiff_args* args) +parse_yaml_param_histogram + (const char* filename, + yaml_document_t* doc, + const yaml_node_t* histo, + struct schiff_param* param) { - struct str buf; - char* tk; - char* tk_ctx; + enum { + LOWER = BIT(0), + UPPER = BIT(1), + PROBAS = BIT(2) + }; + size_t ientry, nentries, nattrs; + size_t i; + int mask = 0; /* Register the parsed histogram attributes */ + double accum_proba; res_T res = RES_OK; - ASSERT(args && str); + ASSERT(filename && doc && histo && param); - str_init(&mem_default_allocator, &buf); + param->distribution = SCHIFF_PARAM_HISTOGRAM; + param->data.histogram.entries = NULL; - /* The distribution was already set to a non sphere distribution. */ - if(args->geom.type != SCHIFF_NONE && args->geom.type != SCHIFF_SPHERE) { + if(histo->type != YAML_MAPPING_NODE) { + log_err(filename, histo, "expecting a mapping of the histogram data.\n"); res = RES_BAD_ARG; goto error; } - res = str_set(&buf, str); - if(res != RES_OK) { - fprintf(stderr, - "Not enough memory. Couldn't parse the sphere distribution `%s'.\n", - str); - goto error; - } + /* Parse the histogram data */ + nattrs = (size_t) + (histo->data.mapping.pairs.top - histo->data.mapping.pairs.start); + FOR_EACH(i, 0, nattrs) { + yaml_node_t *key, *val; + + key = yaml_document_get_node(doc, histo->data.mapping.pairs.start[i].key); + val = yaml_document_get_node(doc, histo->data.mapping.pairs.start[i].value); + ASSERT(key->type == YAML_SCALAR_NODE); + + /* Histogram lower bound */ + if(!strcmp((char*)key->data.scalar.value, "lower")) { + if(mask & LOWER) { + log_err(filename, key, + "the histogram lower bound is already defined.\n"); + res = RES_BAD_ARG; + goto error; + } + mask |= LOWER; + res = parse_yaml_double(filename, val, &param->data.histogram.lower); + if(res != RES_OK) goto error; + + /* Histogram upper bound */ + } else if(!strcmp((char*)key->data.scalar.value, "upper")) { + if(mask & UPPER) { + log_err(filename, key, + "the histogram upper bound is already defined.\n"); + res = RES_BAD_ARG; + goto error; + } + mask |= UPPER; + res = parse_yaml_double(filename, val, &param->data.histogram.upper); + if(res != RES_OK) goto error; + + /* Histogram entries */ + } else if(!strcmp((char*)key->data.scalar.value, "probabilities")) { + if(mask & PROBAS) { + log_err(filename, key, "the histogram data is already defined.\n"); + res = RES_BAD_ARG; + goto error; + } + mask |= PROBAS; - /* Distribution of the sphere radius */ - tk = strtok_r(str_get(&buf), ",", &tk_ctx); - res = parse_param_distribution(tk, &args->geom.data.sphere.radius); - if(res != RES_OK) goto error; + if(val->type != YAML_SEQUENCE_NODE) { + log_err(filename, val, + "expecting a sequence of floating point numbers.\n"); + res = RES_BAD_ARG; + goto error; + } - /* Sphere discretisation */ - tk = strtok_r(NULL, "\0", &tk_ctx); - if(!tk) { - args->geom.data.sphere.nslices = SCHIFF_SPHERE_DEFAULT.nslices; - } else { - res = cstr_to_uint(tk, &args->geom.data.sphere.nslices); - if(res != RES_OK) goto error; - } - args->geom.type = SCHIFF_SPHERE; + nentries = (size_t) + (val->data.sequence.items.top - val->data.sequence.items.start); + if(!sa_add(param->data.histogram.entries, nentries)) { + log_err(filename, val, + "couldn't allocate an histogram with %lu entries.\n", + (unsigned long)nentries); + res = RES_MEM_ERR; + goto error; + } -exit: - str_release(&buf); - return res; -error: - goto exit; -} + /* Parse histogram entries */ + accum_proba = 0; + FOR_EACH(ientry, 0, nentries) { + double proba; + yaml_node_t* node; + node = yaml_document_get_node(doc, val->data.sequence.items.start[ientry]); -static res_T -parse_super_formula(const char* str, struct schiff_param formula[6]) -{ - struct str buf; - char* tk; - char* tk_ctx = NULL; - int i; - res_T res = RES_OK; - ASSERT(formula); + res = parse_yaml_double(filename, node, &proba); + if(res != RES_OK) goto error; - str_init(&mem_default_allocator, &buf); + accum_proba += proba; + param->data.histogram.entries[ientry] = accum_proba; + } - if(!str) { - fprintf(stderr, "Missing super formula.\n"); - res = RES_BAD_ARG; - goto error; - } + /* Normalize the histogram entries */ + FOR_EACH(ientry, 0, nentries) + param->data.histogram.entries[ientry] /= accum_proba; - res = str_set(&buf, str); - if(res != RES_OK) { - fprintf(stderr, - "Not enough memory. Couldn't parse the super formula `%s'.\n", - str); - goto error; + } else { + log_err(filename, key, "unknown histogram data `%s'.\n", + key->data.scalar.value); + res = RES_BAD_ARG; + goto error; + } } - /* Avoid an overzealos GCC warn that assumes that str_get(&buf) may be NULL */ - tk_ctx = NULL; - - tk = strtok_r(str_get(&buf), ",", &tk_ctx); - for(i=0; tk && i < 6; ++i ) { - res = parse_param_distribution(tk, formula + i); - if(res != RES_OK) goto error; - tk = strtok_r(NULL, ",", &tk_ctx); + /* Ensure that the histogram data are all parsed. */ + if(!(mask & LOWER)) { + log_err(filename, histo, "missing the histogram lower parameter.\n"); + res = RES_BAD_ARG; + goto error; + } else if(!(mask & UPPER)) { + log_err(filename, histo, "missing the histogram upper parameter.\n"); + res = RES_BAD_ARG; + goto error; + } else if(!(mask & PROBAS)) { + log_err(filename, histo, "missing the histogram probabilities.\n"); + res = RES_BAD_ARG; + goto error; } - if(i != 6 || tk) { + /* Check the histogram interval */ + if(param->data.histogram.upper <= param->data.histogram.lower) { + log_err(filename, histo, "invalid histogram interval [%g, %g].\n", + param->data.histogram.lower, param->data.histogram.upper); res = RES_BAD_ARG; goto error; } exit: - str_release(&buf); return res; error: - fprintf(stderr, "Invalid super formula `%s'.\n", str); + if(param->data.histogram.entries) { + sa_release(param->data.histogram.entries); + param->data.histogram.entries = NULL; + } goto exit; } static res_T -parse_super_shape_distribution(const char* str, struct schiff_args* args) +parse_yaml_param_distribution + (const char* filename, + yaml_document_t* doc, + const yaml_node_t* node, + struct schiff_param* param) { - struct str buf; - char* tk; - char* tk_ctx; res_T res = RES_OK; - ASSERT(args && str); + ASSERT(filename && doc && node && param); - str_init(&mem_default_allocator, &buf); + if(node->type == YAML_SCALAR_NODE) { /* Floating point constant */ + param->distribution = SCHIFF_PARAM_CONSTANT; + res = parse_yaml_double(filename, node, &param->data.constant); + if(res != RES_OK) return res; - /* The distribution was already set to a non sphere distribution. */ - if(args->geom.type != SCHIFF_NONE && args->geom.type != SCHIFF_SUPER_SHAPE) { - res = RES_BAD_ARG; - goto error; - } + } else if(node->type == YAML_MAPPING_NODE) { + yaml_node_t* key, *val; - res = str_set(&buf, str); - if(res != RES_OK) { - fprintf(stderr, - "Not enough memory. Couldn't parse the super shape distribution `%s'.\n", - str); - goto error; - } + if(node->data.mapping.pairs.top - node->data.mapping.pairs.start != 1) { + log_err(filename, node, + "expecting a mapping from a parameter distribution to its attributes.\n"); + return RES_BAD_ARG; + } - /* Distribution of the super formula 0 parameters */ - tk = strtok_r(str_get(&buf), "#", &tk_ctx); - res = parse_super_formula(tk, args->geom.data.super_shape.formulas[0]); - if(res != RES_OK) goto error; + key = yaml_document_get_node(doc, node->data.mapping.pairs.start[0].key); + val = yaml_document_get_node(doc, node->data.mapping.pairs.start[0].value); + ASSERT(key->type == YAML_SCALAR_NODE); - /* Distribution of the super formula 1 parameters */ - tk = strtok_r(NULL, "#", &tk_ctx); - res = parse_super_formula(tk, args->geom.data.super_shape.formulas[1]); - if(res != RES_OK) goto error; + if(!strcmp((char*)key->data.scalar.value, "lognormal")) { + res = parse_yaml_param_lognormal(filename, doc, val, param); + if(res != RES_OK) return res; + } else if(!strcmp((char*)key->data.scalar.value, "histogram")) { + res = parse_yaml_param_histogram(filename, doc, val, param); + if(res != RES_OK) return res; + } else { + log_err(filename, key, "unknown parameter distribution `%s'.\n", + key->data.scalar.value); + return RES_BAD_ARG; + } - /* Super shape discretisation */ - tk = strtok_r(NULL, "\0", &tk_ctx); - if(!tk) { - args->geom.data.super_shape.nslices = SCHIFF_SUPER_SHAPE_DEFAULT.nslices; } else { - res = cstr_to_uint(tk, &args->geom.data.super_shape.nslices); - if(res != RES_OK) goto error; + log_err(filename, node, "unexpected YAML node.\n"); + return RES_BAD_ARG; } - args->geom.type = SCHIFF_SUPER_SHAPE; -exit: - str_release(&buf); - return res; -error: - goto exit; + return RES_OK; } static res_T -parse_wavelengths(const char* str, struct schiff_args* args) +parse_yaml_cylinder + (const char* filename, + yaml_document_t* doc, + const yaml_node_t* node, + struct schiff_geometry* geom) { - size_t len; + enum { + PROBA = BIT(0), + RADIUS = BIT(1), + HEIGHT = BIT(2), + ASPECT_RATIO = BIT(3) + }; + int mask = 0; /* Register the parsed attributes */ + size_t nattrs; size_t i; res_T res = RES_OK; - ASSERT(args && str); - - /* How many wavelengths are submitted */ - res = cstr_to_list_double(str, NULL, &len, 0); - if(res != RES_OK) goto error; - - /* Reserve the wavelengths memory space */ - sa_clear(args->wavelengths); - args->wavelengths = sa_add(args->wavelengths, len); + ASSERT(filename && doc && node && geom); - /* Read the wavelengths */ - res = cstr_to_list_double(optarg, args->wavelengths, NULL, len); - if(res != RES_OK) goto error; + if(node->type != YAML_MAPPING_NODE) { + log_err(filename, node, "expecting a mapping of cylinder attributes.\n"); + return RES_BAD_ARG; + } - /* Check the validity of read wavelengths */ - FOR_EACH(i, 0, len) { - if(args->wavelengths[i] < 0.0) { + nattrs = (size_t) + (node->data.mapping.pairs.top - node->data.mapping.pairs.start); + + FOR_EACH(i, 0, nattrs) { + yaml_node_t* key; + yaml_node_t* val; + + key = yaml_document_get_node(doc, node->data.mapping.pairs.start[i].key); + val = yaml_document_get_node(doc, node->data.mapping.pairs.start[i].value); + ASSERT(key->type == YAML_SCALAR_NODE); + + /* Distribution probability */ + if(!strcmp((char*)key->data.scalar.value, "proba")) { + double proba; /* TODO use this proba */ + if(mask & PROBA) { + log_err(filename, key, "the cylinder proba is already defined.\n"); + return RES_BAD_ARG; + } + mask |= PROBA; + res = parse_yaml_double(filename, val, &proba); + + /* Cylinder radius */ + } else if(!strcmp((char*)key->data.scalar.value, "radius")) { + if(mask & RADIUS) { + log_err(filename, key, "the cylinder radius is already defined .\n"); + return RES_BAD_ARG; + } + mask |= RADIUS; + res = parse_yaml_param_distribution + (filename, doc, val, &geom->data.cylinder.radius); + + /* Cylinder height */ + } else if(!strcmp((char*)key->data.scalar.value, "height")) { + if(mask & HEIGHT) { + log_err(filename, key, "the cylinder height is already defined .\n"); + return RES_BAD_ARG; + } + mask |= HEIGHT; + res = parse_yaml_param_distribution + (filename, doc, val, &geom->data.cylinder.height); + } else if(!strcmp((char*)key->data.scalar.value, "aspect_ratio")) { + if(mask & ASPECT_RATIO) { + log_err(filename, key, "the aspect_ratio is already defined.\n"); + return RES_BAD_ARG; + } + mask |= ASPECT_RATIO; + res = parse_yaml_double + (filename, val, &geom->data.cylinder.aspect_ratio); + } else { + log_err(filename, key, "unkown cylinder attribute `%s'.\n", + key->data.scalar.value); res = RES_BAD_ARG; - goto error; } + if(res != RES_OK) return res; + } + + /* Ensure the completness of the cylinder distribution */ + if(!(mask & RADIUS)) { + log_err(filename, node, "missing the radius attribute.\n"); + return RES_BAD_ARG; + } else if(!(mask & PROBA)) { + log_err(filename, node, "missing the proba attribute.\n"); + return RES_BAD_ARG; + } else if(!(mask & (HEIGHT|ASPECT_RATIO))) { + log_err(filename, node, "missing the height or aspect_ratio attribute.\n"); + return RES_BAD_ARG; + } else if((mask & HEIGHT) && (mask & ASPECT_RATIO)) { + log_err(filename, node, + "the height and the aspect_ratio attributes cannot be defined " + "simultanously.\n"); + return RES_BAD_ARG; + } + + if(mask & HEIGHT) { + geom->type = SCHIFF_CYLINDER; + } else { + ASSERT(mask & ASPECT_RATIO); + geom->type = SCHIFF_CYLINDER_AS_SPHERE; } -exit: - return res; -error: - goto exit; -} - -static INLINE void -log_yaml_node_error - (const char* filename, - const yaml_node_t* node, - const char* fmt, - ...) -{ - va_list vargs_list; - ASSERT(node && fmt); - - fprintf(stderr, "%s:%lu:%lu: ", - filename, - (unsigned long)node->start_mark.line+1, - (unsigned long)node->start_mark.column+1); - va_start(vargs_list, fmt); - vfprintf(stderr, fmt, vargs_list); - va_end(vargs_list); + return RES_OK; } static res_T -parse_yaml_cylinder +parse_yaml_sphere (const char* filename, yaml_document_t* doc, - const yaml_node_t* cylinder) + const yaml_node_t* node, + struct schiff_geometry* geom) { - size_t iattr, nattrs; + enum { + PROBA = BIT(0), + RADIUS = BIT(1) + }; + int mask = 0; /* Register the parsed attributes */ + size_t nattrs; + size_t i; res_T res = RES_OK; - ASSERT(doc && cylinder); + ASSERT(filename && doc && node && geom); - if(cylinder->type != YAML_MAPPING_NODE) { - log_yaml_node_error(filename, cylinder, - "expecting a YAML mapping of cylinder attributes.\n"); - res = RES_BAD_ARG; - goto error; + if(node->type != YAML_MAPPING_NODE) { + log_err(filename, node, "expecting a mapping of sphere attributes.\n"); + return RES_BAD_ARG; } nattrs = (size_t) - (cylinder->data.mapping.pairs.top - cylinder->data.mapping.pairs.start); + (node->data.mapping.pairs.top - node->data.mapping.pairs.start); - FOR_EACH(iattr, 0, nattrs) { - yaml_node_t* attr; + FOR_EACH(i, 0, nattrs) { + yaml_node_t* key; + yaml_node_t* val; - attr = yaml_document_get_node - (doc, cylinder->data.mapping.pairs.start[iattr].key); - ASSERT(attr->type == YAML_SCALAR_NODE); + key = yaml_document_get_node(doc, node->data.mapping.pairs.start[i].key); + val = yaml_document_get_node(doc, node->data.mapping.pairs.start[i].value); + ASSERT(key->type == YAML_SCALAR_NODE); - if(!strcmp((char*)attr->data.scalar.value, "proba")) { - } else if(!strcmp((char*)attr->data.scalar.value, "radius")) { - } else if(!strcmp((char*)attr->data.scalar.value, "height")) { - } else if(!strcmp((char*)attr->data.scalar.value, "aspect_ratio")) { + if(!strcmp((char*)key->data.scalar.value, "proba")) { + double proba; /* TODO use this proba */ + if(mask & PROBA) { + log_err(filename, node, "the sphere proba is already defined.\n"); + return RES_BAD_ARG; + } + mask |= PROBA; + res = parse_yaml_double(filename, val, &proba); + } else if(!strcmp((char*)key->data.scalar.value, "radius")) { + if(mask & RADIUS) { + log_err(filename, key, "the sphere radius is already defined.\n"); + return RES_BAD_ARG; + } + mask |= RADIUS; + res = parse_yaml_param_distribution + (filename, doc, val, &geom->data.sphere.radius); } else { - log_yaml_node_error(filename, attr, "unkown cylinder attribute `%s'.\n", - attr->data.scalar.value); - res = RES_BAD_ARG; - goto error; + log_err(filename, key, "unkown sphere attribute `%s'.\n", + key->data.scalar.value); + return RES_BAD_ARG; } + if(res != RES_OK) return res; } -exit: - return res; -error: - goto exit; + /* Ensure the completness of the spherical distribution */ + if(!(mask & RADIUS)) { + log_err(filename, node, "missing the radius attribute.\n"); + return RES_BAD_ARG; + } else if(!(mask & PROBA)) { + log_err(filename, node, "missing the proba attribute.\n"); + return RES_BAD_ARG; + } + + return RES_OK; } static res_T parse_yaml(const char* filename) { + struct schiff_geometry geom; yaml_parser_t parser; yaml_document_t doc; yaml_node_t* root; @@ -633,36 +745,33 @@ parse_yaml(const char* filename) ndistribs = (size_t) (root->data.sequence.items.top - root->data.sequence.items.start); FOR_EACH(idistrib, 0, ndistribs) { - yaml_node_t* distrib, *name, *data; + yaml_node_t* distrib, *key, *val; distrib = yaml_document_get_node (&doc, root->data.sequence.items.start[idistrib]); - if(distrib->data.mapping.pairs.top-distrib->data.mapping.pairs.start > 1) { - log_yaml_node_error(filename, distrib, "expecting only one mapping pair\n"); - res = RES_BAD_ARG; - goto error; - } - if(distrib->type != YAML_MAPPING_NODE) { - log_yaml_node_error(filename, distrib, "expecting a YAML mapping.\n"); + if(distrib->type != YAML_MAPPING_NODE + || distrib->data.mapping.pairs.top - distrib->data.mapping.pairs.start > 1) { + log_err(filename, distrib, + "expecting a mapping of the geometry distribution to its parameters\n"); res = RES_BAD_ARG; goto error; } - name = yaml_document_get_node(&doc, distrib->data.mapping.pairs.start[0].key); - data = yaml_document_get_node(&doc, distrib->data.mapping.pairs.start[0].value); - ASSERT(name->type == YAML_SCALAR_NODE); + key = yaml_document_get_node(&doc, distrib->data.mapping.pairs.start[0].key); + val = yaml_document_get_node(&doc, distrib->data.mapping.pairs.start[0].value); + ASSERT(key->type == YAML_SCALAR_NODE); - if(!strcmp((const char*)name->data.scalar.value, "cylinder")) { - res = parse_yaml_cylinder(filename, &doc, data); - /* TODO parse cylinder */ - } else if(!strcmp((const char*)name->data.scalar.value, "sphere")) { - /* TODO parse sphere */ + if(!strcmp((char*)key->data.scalar.value, "cylinder")) { + res = parse_yaml_cylinder(filename, &doc, val, &geom); + } else if(!strcmp((char*)key->data.scalar.value, "sphere")) { + res = parse_yaml_sphere(filename, &doc, val, &geom); } else { - log_yaml_node_error(filename, name, "unknown distribution type `%s'.\n", - name->data.scalar.value); + log_err(filename, key, "unknown distribution `%s'.\n", + key->data.scalar.value); res = RES_BAD_ARG; goto error; } + if(res != RES_OK) goto error; } exit: @@ -689,10 +798,8 @@ schiff_args_init ASSERT(argc && argv && args); *args = SCHIFF_ARGS_NULL; - while((opt = getopt(argc, argv, "c:C:d:g:G:hi:n:o:qs:u:w:")) != -1) { + while((opt = getopt(argc, argv, "d:g:G:hi:n:o:qw:")) != -1) { switch(opt) { - case 'c': res = parse_cylinder_distribution(optarg, args); break; - case 'C': res = parse_cylinder_as_sphere_distribution(optarg, args); break; case 'd': res = cstr_to_uint(optarg, &args->ndirs); break; case 'g': res = cstr_to_uint(optarg, &args->ngeoms); break; case 'G': res = cstr_to_uint(optarg, &args->ngeoms_dump); break; @@ -710,8 +817,6 @@ schiff_args_init break; case 'o': args->output_filename = optarg; break; case 'q': quiet = 1; break; - case 's': res = parse_sphere_distribution(optarg, args); break; - case 'u': res = parse_super_shape_distribution(optarg, args); break; case 'w': res = parse_wavelengths(optarg, args); break; default: res = RES_BAD_ARG; break; } @@ -752,7 +857,6 @@ schiff_args_init fprintf(stderr, "Enter the optical properties. Type ^D (i.e. CTRL+d) to stop:\n"); #endif - } res = schiff_optical_properties_load_stream(&args->properties, stdin, "stdin"); } else { diff --git a/src/schiff_args.h b/src/schiff_args.h @@ -48,7 +48,7 @@ static const struct schiff_args SCHIFF_ARGS_NULL = { NULL, /* Output filename */ NULL, /* List of optical properties */ NULL, /* List of wavelength to integrate */ - SCHIFF_GEOMETRY_NULL__, + SCHIFF_GEOMETRY_NULL__, /* List of Schiff geometries */ 0, /* # Dumped geometries */ 1000, /* # Sampled geometries */ 100, /* # Sampled directions per gemetry */ diff --git a/src/schiff_geometry.c b/src/schiff_geometry.c @@ -30,6 +30,9 @@ #include "schiff_mesh.h" #include "schiff_optical_properties.h" +#include <rsys/algorithm.h> +#include <rsys/stretchy_array.h> + #include <star/s3d.h> #include <star/ssp.h> #include <star/sschiff.h> @@ -79,6 +82,48 @@ struct geometry_distribution_context { /******************************************************************************* * Helper functions ******************************************************************************/ +static int +cmp_double(const void* op0, const void* op1) +{ + const double* a = (const double*)op0; + const double* b = (const double*)op1; + return *a < *b ? -1 : (*a > *b ? 1 : 0); +} + +static double +histogram_sample + (const double* entries, + const size_t nentries, + const double lower, + const double upper, + const double u) +{ + const double* find; + double step; + double from, to; + double v; + ASSERT(entries && nentries && lower < upper && u >= 0 && u < 1.0); + + find = search_lower_bound(&u, entries, nentries, sizeof(double), cmp_double); + + if(!find) /* Handle numerical issue */ + find = entries + (nentries - 1); + + step = (upper - lower) / (double)nentries; /* Size of an histogram interval */ + from = lower + (double)(find - entries) * step; + to = from + step; + + /* Linear interpolation */ + v = (u - from) / step; + if(eq_eps(v, 0, 1.e-6)) { + return from; + } else if(eq_eps(v, 1, 1.e-6)) { + return to; + } else { + return v*to + (1 - v)*from; + } +} + static INLINE double eval_param(const struct schiff_param* param, struct ssp_rng* rng) { @@ -94,8 +139,12 @@ eval_param(const struct schiff_param* param, struct ssp_rng* rng) log(param->data.lognormal.sigma)); break; case SCHIFF_PARAM_HISTOGRAM: - val = schiff_histogram_sample - (&param->data.histogram, ssp_rng_canonical(rng)); + val = histogram_sample + (param->data.histogram.entries, + sa_size(param->data.histogram.entries), + param->data.histogram.lower, + param->data.histogram.upper, + ssp_rng_canonical(rng)); break; default: FATAL("Unreachable code.\n"); break; } diff --git a/src/schiff_geometry.h b/src/schiff_geometry.h @@ -30,7 +30,6 @@ #define SCHIFF_GEOMETRY_H #include <rsys/rsys.h> -#include "schiff_histogram.h" enum schiff_param_distribution { SCHIFF_PARAM_CONSTANT, @@ -44,7 +43,7 @@ struct schiff_param { union { double constant; struct { double zeta, sigma; } lognormal; - struct darray_hentry histogram; + struct { double *entries, lower, upper; } histogram; } data; }; #define SCHIFF_PARAM_DEFAULT__ {SCHIFF_PARAM_CONSTANT, {1.0}} diff --git a/src/schiff_histogram.c b/src/schiff_histogram.c @@ -1,302 +0,0 @@ -/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com) - * - * 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 "schiff_histogram.h" -#include "schiff_streamline.h" - -#include <rsys/algorithm.h> -#include <rsys/str.h> -#include <rsys/cstr.h> -#include <rsys/stretchy_array.h> - -/******************************************************************************* - * Helper functions - ******************************************************************************/ -static res_T -parse_histogram_proba - (const char* str, - const char* filename, - const size_t iline, - double* proba) -{ - char buf[128]; - double d; - char* tk; - res_T res = RES_OK; - ASSERT(proba && filename && str); - - if(strlen(str) >= sizeof(buf) - 1/*NULL char*/) { - fprintf(stderr, - "%s:%lu: not enough memory: cannot parse the schiff histogram value `%s'.\n", - filename, (unsigned long)iline, str); - return RES_MEM_ERR; - } - - strncpy(buf, str, sizeof(buf)); - 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(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; -} - -static INLINE int -cmp_entry(const void* a, const void* b) -{ - const struct schiff_hentry* key = a; - const struct schiff_hentry* entry = b; - ASSERT(a && b); - return key->accum_proba < entry->accum_proba - ? -1 : key->accum_proba > 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; - goto error; - } - *size = (unsigned long)l; - -exit: - str_release(&buf); - return res; -error: - goto exit; -} - -/******************************************************************************* - * Local functions - ******************************************************************************/ -res_T -schiff_histogram_load(struct darray_hentry* histo, const char* filename) -{ - FILE* fp = NULL; - res_T res = RES_OK; - ASSERT(histo && filename); - - fp = fopen(filename, "r"); - if(!fp) { - fprintf(stderr, "Cannot open the file of optical properties `%s'.\n", - filename); - res = RES_IO_ERR; - goto error; - } - - res = schiff_histogram_load_stream(histo, fp, filename); - if(res != RES_OK) goto error; - -exit: - if(fp) fclose(fp); - return res; -error: - goto exit; -} - -res_T -schiff_histogram_load_stream - (struct darray_hentry* histo, - FILE* stream, - const char* stream_name) -{ - struct schiff_streamline streamline; - size_t iline; - size_t i; - char* line; - double bounds[2]; - double step; - double val; - double accum_proba; - unsigned long nentries; - res_T res = RES_OK; - ASSERT(histo && stream_name && stream_name); - - darray_hentry_clear(histo); - schiff_streamline_init(&streamline); - - iline = 1; - - /* 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_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; - ASSERT(step >= 0); - val = bounds[0]; - accum_proba = 0; - for(i=0;(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; - - 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(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_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; -error: - darray_hentry_clear(histo); - goto exit; -} - -double -schiff_histogram_sample(const struct darray_hentry* histo, const double u) -{ - const struct schiff_hentry* entries; - const struct schiff_hentry* find; - struct schiff_hentry samp; - size_t nentries; - double v; - ASSERT(histo && u >= 0.0 && u < 1.0); - - entries = darray_hentry_cdata_get(histo); - nentries = darray_hentry_size_get(histo); - ASSERT(nentries); - - samp.accum_proba = u; - find = search_lower_bound(&samp, entries, nentries, - sizeof(struct schiff_hentry), cmp_entry); - - 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 @@ -1,56 +0,0 @@ -/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com) - * - * 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. */ - -#ifndef SCHIFF_HISTOGRAM_H -#define SCHIFF_HISTOGRAM_H - -#include <rsys/dynamic_array.h> - -struct schiff_hentry { double value, accum_proba; }; -#define DARRAY_NAME hentry -#define DARRAY_DATA struct schiff_hentry -#include <rsys/dynamic_array.h> - -extern LOCAL_SYM res_T -schiff_histogram_load - (struct darray_hentry* histo, - const char* filename); - -extern LOCAL_SYM res_T -schiff_histogram_load_stream - (struct darray_hentry* histo, - FILE* stream, - const char* stream_name); - -extern LOCAL_SYM double -schiff_histogram_sample - (const struct darray_hentry* histo, - const double u); /* in [0, 1[ */ - -#endif /* SCHIFF_HISTOGRAM_H */ -