commit bc749a762177ed5c0f5f7dd9cbb46ec51652667e
parent 5e365c2c362833832829289c92a23fa25b41deb9
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 23 Feb 2016 14:21:32 +0100
First draft of the estimation of the phase function connector values
Diffstat:
2 files changed, 158 insertions(+), 74 deletions(-)
diff --git a/src/sschiff_estimator.c b/src/sschiff_estimator.c
@@ -47,6 +47,7 @@
#include <omp.h>
#include <gsl/gsl_sf.h>
+#include <gsl/gsl_multiroots.h>
#define MAX_FAILURES 10
@@ -191,6 +192,19 @@ build_transform
transform[11] = -f3_dot(basis+6, pos);
}
+static void
+get_mc_value
+ (const struct mc_data* data,
+ const size_t nrealisations,
+ struct sschiff_state* value)
+{
+ ASSERT(data && value);
+ value->E = data->weight / (double)nrealisations;
+ value->V = data->sqr_weight / (double)nrealisations
+ - (data->weight * data->weight) / (double)(nrealisations * nrealisations);
+ value->SE = sqrt(value->V / (double)nrealisations);
+}
+
/* Dump a thumbnail of the sampled geometry in given sampled direction */
static INLINE void
draw_thumbnail
@@ -447,102 +461,178 @@ accum_differential_cross_section
}
}
-#if 0
-struct genevieve_arg {
- double limit_angle;
+struct function_arg {
double scattering_cross_section;
- double limit_differential_cross_section;
- double limit_differential_cross_section_cumulative;
+ double angle;
+ double differential_cross_section;
+ double differential_cross_section_cumulative;
+
+ /* Precomputed values */
+ double cos_2_angle;
+ double cos_angle;
+ double cos_half_angle;
+ double sin_half_angle;
+ double sqr_cos_angle;
+ double sqr_sin_half_angle;
};
-int
-genevieve(const gsl_vector* x, void* params, gsl_vector* f)
+static int
+function_fdf(const gsl_vector* vec, void* params, gsl_vector* res, gsl_matrix* J)
{
- struct genevieve_arg* arg = params;
- const double x = gsl_vector_get(x, 0);
- double f0;
- double square_cos_limit_angle;
- double square_half_sin_limit_angle;
-
- if(x==2 || x==4 || x==6)
+ struct function_arg* arg = params;
+ const double n = gsl_vector_get(vec, 0);
+ double u, u_prime;
+ double v, v_prime;
+ double coef;
+ double f, df;
+
+ if(n==2 || n==4 || n==6)
return GSL_EDOM;
- square_cos_limit_angle = cos(limit_angle);
- square_cos_limit_angle *= square_cos_limit_angle;
+ coef = (4*PI * arg->differential_cross_section) / (1+arg->sqr_cos_angle);
+
+ u = ( arg->sqr_sin_half_angle
+ * ( (n*n - 6*n + 8) * arg->cos_2_angle
+ + 3*n*n
+ - 8*(n - 2) * arg->cos_angle
+ - 26*n + 72 )
+ - pow(arg->sin_half_angle, n) * (4*n*n - 24*n + 40) );
+ v = ((n - 2) * (n - 4) * (n - 6));
- square_half_sin_limit_angle = sin(limit_angle) / 2;
- square_half_sin_limit_angle *= square_half_sin_limit_angle;
+ u_prime = arg->sqr_sin_half_angle
+ * ((2*n - 6) * arg->cos_2_angle + 6*n - 8*arg->cos_angle - 26)
+ - n/2.0 * arg->cos_half_angle * pow(arg->sin_half_angle, n-1)
+ * (4*n*n - 24*n + 40) - pow(arg->sin_half_angle, n) * (8*n - 24);
+ v_prime = (n-4)*(n-6) + (n-2)*(n-6) + (n-2)*(n-4);
- f0 = arg->limit_differential_cross_section_cumulative
+ f = arg->differential_cross_section_cumulative
- arg->scattering_cross_section
- + (4*PI * arg->limit_differential_cross_section) / (1+square_cos_limit_angle)
- * (square_half_sin_limit_angle * ((x*x - 6*x + 8) * cos(2*arg->limit_angle) + 3*x*x - 8 * (x - 2) * cos(arg->limit_angle) - 26*x + 72) - pow(sin(arg->limit_angle/2.0), x) * (4*x*x - 24*x + 40))
- / ((x-2)*(x-4)*(x-6));
+ + coef * u / v;
+ df = (v*u_prime - u*v_prime) / (v*v);
+ df = df * coef;
+ gsl_vector_set(res, 0, f);
+ gsl_matrix_set(J, 0, 0, df);
return GSL_SUCCESS;
}
-static void
-phase_function_post_process()
+static res_T
+phase_function_post_process
+ (struct sschiff_device* dev,
+ const size_t iwavelength,
+ const struct geometry_context* ctx,
+ struct sschiff_estimator* estimator)
{
- const gsl_multiroot_fsolver* solver = NULL;
- const gsl_multiroot_function func;
- double n;
- double r;
- struct genevieve_param* arg;
- struct gsl_vector* x = NULL;
- struct gsl_vector* root;
+ gsl_multiroot_fdfsolver* solver = NULL;
+ gsl_multiroot_function_fdf func;
+ gsl_vector* init_val = NULL;
+ gsl_vector* root = NULL;
+ struct mc_accum* accums;
+ struct function_arg arg;
+ double n, r; /* Estimated value */
+ size_t ilimit_angle; /* Index of the limit angle of the current wavelength */
+ int status;
+ int iter;
+ const int max_iter = 1000; /* Maximum number of Newton iterations */
res_T res = RES_OK;
+ ASSERT(dev);
- solver = gsl_multiroot_fsolver_alloc(gsl_multiroot_fsolver_dnewton, 1);
+ /* Create a Newton-Raphson nonlinear solver without derivatives */
+ solver = gsl_multiroot_fdfsolver_alloc(gsl_multiroot_fdfsolver_newton, 1);
if(!solver) {
log_error(dev,
- "Not enough memory. Couldn't allocate the GSL Newton-Raphson solver.\n");
+ "Not enough memory. Couldn't allocate the GSL Newton solver.\n");
res = RES_MEM_ERR;
goto error;
}
- x = gsl_vector_alloc(1);
- if(!x) {
+ /* Allocate the result vector. */
+ init_val = gsl_vector_alloc(1);
+ if(!init_val) {
log_error(dev,
- "Not enough memory. Couldn't allocat the initial GSL vector.\n");
+ "Not enough memory. Couldn't allocate the GSL init vector.\n");
res = RES_MEM_ERR;
goto error;
}
- /* TODO fill arg */
- func.f = genevieve;
- func.n = 1;
- func.params = &arg;
+ /* Fetch attribute */
+ accums = estimator->accums + iwavelength;
+ ilimit_angle = ctx->limit_angles[iwavelength]-1;
+
+ /* Fill system arguments */
+ arg.angle = estimator->angles[ilimit_angle];
+ arg.scattering_cross_section =
+ accums->scattering_cross_section.weight
+ / (double)estimator->nrealisations;
+ arg.differential_cross_section =
+ accums->differential_cross_section[ilimit_angle].weight
+ / (double)estimator->nrealisations;
+ arg.differential_cross_section_cumulative =
+ accums->differential_cross_section_cumulative[ilimit_angle].weight
+ / (double)estimator->nrealisations;
- gsl_vector_set(x, 0, 3.0);
- gsl_multiroot_fsolver_set(solver, f, x);
-
- FOR_EACH(iwlen, 0, nwavelengths) {
- int status;
- int iter;
- do {
- status = gsl_multiroot_fsolver_iterate(solver);
- if(states != GSL_SUCCESS) break; /* ERROR */
+ {
+ struct sschiff_state state;
+ get_mc_value
+ (&accums->differential_cross_section[ilimit_angle],
+ estimator->nrealisations,
+ &state);
+ }
+ /* Precompute some values to speed up Newton iterations */
+ arg.cos_2_angle = cos(2.0 * arg.angle);
+ arg.cos_angle = cos(arg.angle);
+ arg.cos_half_angle = cos(arg.angle * 0.5);
+ arg.sin_half_angle = sin(arg.angle * 0.5);
+ arg.sqr_cos_angle = arg.cos_angle * arg.cos_angle;
+ arg.sqr_sin_half_angle = arg.sin_half_angle * arg.sin_half_angle;
+
+ /* Setup system functions */
+ func.f = NULL;
+ func.df = NULL;
+ func.fdf = function_fdf;
+ func.n = 1; /* Number of equations */
+ func.params = &arg; /* System parameters */
+
+ /* Initialise the solver */
+ gsl_vector_set(init_val, 0, 3.0); /* n */
+ gsl_multiroot_fdfsolver_set(solver, &func, init_val);
+ gsl_set_error_handler_off();
+
+ /* Launch the Newton estimation. Iterate up to `max_iter'or until the
+ * estimation is "good enough". */
+ iter = 0;
+ do {
+ status = gsl_multiroot_fdfsolver_iterate(solver);
+ if(status == GSL_SUCCESS) {
status = gsl_multiroot_test_residual(solver->f, 1.e-3);
- } while(status == GSL_CONTINUE && iter < 1000);
+ }
+ } while(status == GSL_CONTINUE && ++iter < max_iter);
+
+ if(status != GSL_SUCCESS) {
+ /* TODO enhance the error message. */
+ log_error(dev, "Error in Newton estimation - GSL error: %s.\n",
+ gsl_strerror(status));
+ res = RES_BAD_OP;
+ goto error;
}
- root = gsl_multiroot_fsolver_root(solver->x);
+ /* Retrieve the estimated "n" value */
+ root = gsl_multiroot_fdfsolver_root(solver);
n = gsl_vector_get(root, 0);
+ printf("%g\n", n);
- r = 2 * limit_differential_cross_section * pow(sin(limit_angle/2), n)
- / (1 + pow(cos(limit_angle), 2));
+ /* Compute r from the estimated n */
+ r = 2 * arg.differential_cross_section * pow(arg.sin_half_angle, n)
+ / (1 + arg.sqr_cos_angle);
exit:
- if(solver) gsl_multiroot_fsolver_free(solver);
- if(x) gsl_vector_free(x);
+ if(solver) gsl_multiroot_fdfsolver_free(solver);
+ if(init_val) gsl_vector_free(init_val);
return res;
error:
goto exit;
}
-#endif
static res_T
accum_monte_carlo_weight
@@ -763,19 +853,6 @@ error:
goto exit;
}
-static void
-get_mc_value
- (const struct mc_data* data,
- const size_t nrealisations,
- struct sschiff_state* value)
-{
- ASSERT(data && value);
- value->E = data->weight / (double)nrealisations;
- value->V = data->sqr_weight / (double)nrealisations
- - (data->weight * data->weight) / (double)(nrealisations * nrealisations);
- value->SE = sqrt(value->V / (double)nrealisations);
-}
-
static char
check_distribution(struct sschiff_geometry_distribution* distrib)
{
@@ -1125,6 +1202,7 @@ sschiff_integrate
struct ssp_rng** rngs = NULL;
struct ssp_rng_proxy* rng_proxy = NULL;
size_t i;
+ size_t iwlen;
int igeom;
ATOMIC res = (ATOMIC)RES_OK;
memset(&geom_ctx, 0, sizeof(geom_ctx));
@@ -1211,7 +1289,6 @@ sschiff_integrate
/* Merge the per thread integration ressults */
FOR_EACH(i, 0, dev->nthreads) {
- size_t iwlen;
FOR_EACH(iwlen, 0, nwavelengths) {
size_t iangle;
#define MC_ACCUM(Data) { \
@@ -1232,6 +1309,13 @@ sschiff_integrate
}
}
+#if 0
+ FOR_EACH(i, 0, nwavelengths) {
+ res = phase_function_post_process(dev, i, &geom_ctx, estimator);
+ if(res != RES_OK) goto error;
+ }
+#endif
+
exit:
geometry_context_release(&geom_ctx);
if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
diff --git a/src/test_sschiff_estimator_cylinder.c b/src/test_sschiff_estimator_cylinder.c
@@ -553,9 +553,9 @@ main(int argc, char** argv)
2.010619438e+2, 2.018997018e+2, 2.027374599e+2
};
const size_t nx = sizeof(x)/sizeof(double);
- const size_t nscatt_angles = 3;
- const size_t ngeoms = 100;
- const size_t ndirs = 10;
+ const size_t nscatt_angles = 1000;
+ const size_t ngeoms = 10000;
+ const size_t ndirs = 100;
size_t i;
(void)argc, (void)argv;
@@ -567,7 +567,7 @@ main(int argc, char** argv)
cylinder_init(&sampler_ctx.geometry, 64);
- FOR_EACH(i, 0, nx) {
+ FOR_EACH(i, 10, nx) {
const double wavelength = 0.6; /* In micron */
struct sschiff_cross_section cross_section;
double interval[2];
