commit 2f167ea14d8d18716a0f35b1687438f57843c266
parent 4a0eadada0928281eb8a9bcb5d47037f8ecc9c8f
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 23 Mar 2016 10:38:45 +0100
Add the gaussian parameter distribution
Diffstat:
3 files changed, 103 insertions(+), 90 deletions(-)
diff --git a/src/schiff_args.c b/src/schiff_args.c
@@ -240,79 +240,81 @@ parse_yaml_uint
}
static res_T
-parse_yaml_param_lognormal
+parse_yaml_param_mu_sigma
(const char* filename,
yaml_document_t* doc,
- const yaml_node_t* lognormal,
+ const yaml_node_t* distrib,
+ const char* distrib_name,
const double min_val, /* Minimum valid value fo the parameter */
const double max_val, /* Maximum valid value of the parameter */
- struct schiff_param* param)
+ double* mu,
+ double* sigma)
{
enum {
- ZETA = BIT(0),
+ MU = BIT(0),
SIGMA = BIT(1)
};
size_t nattrs;
size_t i;
- int mask = 0; /* Register the parsed lognormal attributes */
+ int mask = 0; /* Register the parsed attributes */
res_T res = RES_OK;
- ASSERT(filename && doc && lognormal && param && min_val < max_val);
+ ASSERT(filename && doc && distrib && distrib_name);
+ ASSERT(min_val < max_val && mu && sigma);
- if(lognormal->type != YAML_MAPPING_NODE) {
- log_err(filename, lognormal,
- "expecting a mapping of the lognormal attributes.\n");
+ if(distrib->type != YAML_MAPPING_NODE) {
+ log_err(filename, distrib,
+ "expecting a mapping of the %s attributes.\n", distrib_name);
return RES_BAD_ARG;
}
- param->distribution = SCHIFF_PARAM_LOGNORMAL;
-
/* Parse the log normal attributes */
nattrs = (size_t)
- (lognormal->data.mapping.pairs.top - lognormal->data.mapping.pairs.start);
+ (distrib->data.mapping.pairs.top - distrib->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);
+ key=yaml_document_get_node(doc, distrib->data.mapping.pairs.start[i].key);
+ val=yaml_document_get_node(doc, distrib->data.mapping.pairs.start[i].value);
ASSERT(key->type == YAML_SCALAR_NODE);
- /* lognormal mu attribute */
+ #define SETUP_MASK(Flag, Name) { \
+ if(mask & Flag) { \
+ log_err(filename, key, "the "Name" %s attribute is already defined.\n",\
+ distrib_name); \
+ return RES_BAD_ARG; \
+ } \
+ mask |= Flag; \
+ } (void)0
+
+ /* mu attribute */
if(!strcmp((char*)key->data.scalar.value, "mu")) {
- if(mask & ZETA) {
- log_err(filename, key,
- "the `mu' lognormal attribute is already defined.\n");
- return RES_BAD_ARG;
- }
- mask |= ZETA;
- res = parse_yaml_double /* mean value */
- (filename, val, min_val, max_val, ¶m->data.lognormal.mu);
+ SETUP_MASK(MU, "`mu'");
+ res = parse_yaml_double(filename, val, min_val, max_val, mu);
- /* lognormal sigma attribute */
+ /* 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, DBL_MIN, DBL_MAX, ¶m->data.lognormal.sigma);
+ SETUP_MASK(SIGMA, "`sigma'");
+ res = parse_yaml_double(filename, val, DBL_MIN, DBL_MAX, sigma);
/* unknown attribute */
} else {
- log_err(filename, key, "unknown lognormal attribute `%s'.\n",
- key->data.scalar.value);
+ log_err(filename, key, "unknown %s attribute `%s'.\n",
+ distrib_name, key->data.scalar.value);
return RES_BAD_ARG;
}
if(res != RES_OK) return res;
+
+ #undef SETUP_MASK
}
- /* Ensure that the lognormal attributes are all parsed */
- if(!(mask & ZETA)) {
- log_err(filename, lognormal, "missing the mu lognormal attribute.\n");
+ /* Ensure that the attributes are all parsed */
+ if(!(mask & MU)) {
+ log_err(filename, distrib, "missing the `mu' %s attribute.\n",
+ distrib_name);
return RES_BAD_ARG;
} else if(!(mask & SIGMA)) {
- log_err(filename, lognormal, "missing the sigma lognormal attribute.\n");
+ log_err(filename, distrib, "missing the `sigma' %s attribute.\n",
+ distrib_name);
return RES_BAD_ARG;
}
@@ -359,40 +361,31 @@ parse_yaml_param_histogram
val = yaml_document_get_node(doc, histo->data.mapping.pairs.start[i].value);
ASSERT(key->type == YAML_SCALAR_NODE);
+ #define SETUP_MASK(Flag, Name) { \
+ if(mask & Flag) { \
+ log_err(filename, key, "the "Name" is already defined.\n"); \
+ return RES_BAD_ARG; \
+ } \
+ mask |= Flag; \
+ } (void)0
+
/* 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;
+ SETUP_MASK(LOWER, "histogram lower bound");
res = parse_yaml_double
(filename, val, min_val, max_val, ¶m->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;
+ SETUP_MASK(UPPER, "histogram upper bound");
res = parse_yaml_double
(filename, val, min_val, max_val, ¶m->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;
+ SETUP_MASK(PROBAS, "histogram data");
if(val->type != YAML_SEQUENCE_NODE) {
log_err(filename, val,
@@ -435,6 +428,8 @@ parse_yaml_param_histogram
res = RES_BAD_ARG;
goto error;
}
+
+ #undef SETUP_MASK
}
/* Ensure that the histogram data are all parsed. */
@@ -501,14 +496,27 @@ parse_yaml_param_distribution
val = yaml_document_get_node(doc, node->data.mapping.pairs.start[0].value);
ASSERT(key->type == YAML_SCALAR_NODE);
+ /* Lognormal distribution */
if(!strcmp((char*)key->data.scalar.value, "lognormal")) {
- res = parse_yaml_param_lognormal
- (filename, doc, val, min_val, max_val, param);
+ res = parse_yaml_param_mu_sigma(filename, doc, val, "lognormal", min_val,
+ max_val, ¶m->data.lognormal.mu, ¶m->data.lognormal.sigma);
+ if(res != RES_OK) return res;
+ param->distribution = SCHIFF_PARAM_LOGNORMAL;
+
+ /* Gaussian distribution */
+ } else if(!strcmp((char*)key->data.scalar.value, "gaussian")) {
+ res = parse_yaml_param_mu_sigma(filename, doc, val, "gaussian", min_val,
+ max_val, ¶m->data.gaussian.mu, ¶m->data.lognormal.sigma);
if(res != RES_OK) return res;
+ param->distribution = SCHIFF_PARAM_GAUSSIAN;
+
+ /* Histogram distribution */
} else if(!strcmp((char*)key->data.scalar.value, "histogram")) {
res = parse_yaml_param_histogram
(filename, doc, val, min_val, max_val, param);
if(res != RES_OK) return res;
+
+ /* Unknown distribution */
} else {
log_err(filename, key, "unknown parameter distribution `%s'.\n",
key->data.scalar.value);
@@ -563,7 +571,7 @@ parse_yaml_ellipsoid
#define SETUP_MASK(Flag, Name) { \
if(mask & Flag) { \
- log_err(filename, node, "the "Name" is already defined.\n"); \
+ log_err(filename, key, "the "Name" is already defined.\n"); \
return RES_BAD_ARG; \
} \
mask |= Flag; \
@@ -610,13 +618,14 @@ parse_yaml_ellipsoid
return RES_BAD_ARG;
}
if(res != RES_OK) return res;
+
+ #undef SETUP_MASK
}
- #undef SETUP_MASK
/* Ensure the completeness of the ellipsoid distribution */
if((mask & (LENGTH_a | LENGTH_c))) {
if(mask & (RADIUS | ASPECT_RATIO)) {
- log_err(filename, node,
+ log_err(filename, node,
"the ellipsoid semi-principal axis and its equivalent sphere cannot "
"be defined simultaneously.\n");
return RES_BAD_ARG;
@@ -632,7 +641,7 @@ parse_yaml_ellipsoid
geom->type = SCHIFF_ELLIPSOID;
} else if (mask & (RADIUS | ASPECT_RATIO)) {
if(mask & (LENGTH_a | LENGTH_c)) {
- log_err(filename, node,
+ log_err(filename, node,
"the ellipsoid semi-principal axis and its equivalent sphere cannot "
"be defined simultaneously.\n");
return RES_BAD_ARG;
@@ -692,53 +701,40 @@ parse_yaml_cylinder
val = yaml_document_get_node(doc, node->data.mapping.pairs.start[i].value);
ASSERT(key->type == YAML_SCALAR_NODE);
+ #define SETUP_MASK(Flag, Name) { \
+ if(mask & Flag) { \
+ log_err(filename, key, "the "Name" is already defined.\n"); \
+ return RES_BAD_ARG; \
+ } \
+ mask |= Flag; \
+ } (void)0
+
/* Distribution probability */
if(!strcmp((char*)key->data.scalar.value, "proba")) {
- if(mask & PROBA) {
- log_err(filename, key, "the cylinder proba is already defined.\n");
- return RES_BAD_ARG;
- }
- mask |= PROBA;
+ SETUP_MASK(PROBA, "cylinder proba");
res = parse_yaml_double(filename, val, DBL_MIN, DBL_MAX, geom_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;
+ SETUP_MASK(RADIUS, "cylinder radius");
res = parse_yaml_param_distribution
(filename, doc, val, DBL_MIN, DBL_MAX, &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;
+ SETUP_MASK(HEIGHT, "cylinder height");
res = parse_yaml_param_distribution
(filename, doc, val, DBL_MIN, DBL_MAX, &geom->data.cylinder.height);
/* Cylinder aspect ratio */
} 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;
+ SETUP_MASK(ASPECT_RATIO, "cylinder height/radius aspect_ratio");
res = parse_yaml_double
(filename, val, DBL_MIN, DBL_MAX, &geom->data.cylinder.aspect_ratio);
/* # slices used to discretized the triangular cylinder */
} else if(!strcmp((char*)key->data.scalar.value, "slices")) {
- if(mask & SLICES) {
- log_err(filename, key,
- "the cylinder number of slices is already defined.\n");
- return RES_BAD_ARG;
- }
- mask |= SLICES;
+ SETUP_MASK(SLICES, "cylinder number of slices");
res = parse_yaml_uint
(filename, val, 4, 32768, &geom->data.cylinder.nslices);
@@ -749,6 +745,8 @@ parse_yaml_cylinder
res = RES_BAD_ARG;
}
if(res != RES_OK) return res;
+
+ #undef SETUP_MASK
}
/* Ensure the completness of the cylinder distribution */
diff --git a/src/schiff_geometry.c b/src/schiff_geometry.c
@@ -128,6 +128,15 @@ histogram_sample
}
}
+static FINLINE double
+ran_positive_gaussian(struct ssp_rng* rng, const double mu, const double sigma)
+{
+ double val;
+ ASSERT(rng && mu > 0);
+ do { val = ssp_ran_gaussian(rng, mu, sigma); } while(val <= 0);
+ return val;
+}
+
static INLINE double
eval_param(const struct schiff_param* param, struct ssp_rng* rng)
{
@@ -142,6 +151,10 @@ eval_param(const struct schiff_param* param, struct ssp_rng* rng)
val = ssp_ran_lognormal(rng, log(param->data.lognormal.mu),
log(param->data.lognormal.sigma));
break;
+ case SCHIFF_PARAM_GAUSSIAN:
+ val = ran_positive_gaussian
+ (rng, param->data.gaussian.mu, param->data.gaussian.sigma);
+ break;
case SCHIFF_PARAM_HISTOGRAM:
val = histogram_sample
(param->data.histogram.entries,
diff --git a/src/schiff_geometry.h b/src/schiff_geometry.h
@@ -22,6 +22,7 @@
enum schiff_param_distribution {
SCHIFF_PARAM_CONSTANT,
SCHIFF_PARAM_LOGNORMAL,
+ SCHIFF_PARAM_GAUSSIAN,
SCHIFF_PARAM_HISTOGRAM,
SCHIFF_PARAM_NONE
};
@@ -31,6 +32,7 @@ struct schiff_param {
union {
double constant;
struct { double mu, sigma; } lognormal;
+ struct { double mu, sigma; } gaussian;
struct { double *entries, lower, upper; } histogram;
} data;
};
