commit 52687f24428cca216e2031bbe473e648e0e96c89
parent 66cb2687eda5448174d0af13f3ccd3a20cf4eb3d
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 15 Sep 2021 11:17:27 +0200
Upd the way the compute_trans_sensib process is written
Diffstat:
1 file changed, 117 insertions(+), 79 deletions(-)
diff --git a/src/sgs_compute_trans_sensib.c b/src/sgs_compute_trans_sensib.c
@@ -26,16 +26,18 @@
#include <star/smc.h>
#include <star/ssp.h>
+/* Sugar syntax */
enum {X, Y, Z};
+enum {WEIGHT, SENSIB, WEIGHTS_COUNT__};
+/* FIXME this should be variables */
static const double RECV_MIN[2] = {0.5, 1.5};
static const double RECV_MAX[2] = {1.5, 2.5};
-static const double EMIT_THRESHOLD = 2;
-static const double V[3] = {0, 0, 1};
-
-/* FIXME this should be a variable */
-static const double EMIT_SZ[3] = {4, 4, 0};
+static const double EMIT_E_THRESHOLD = 2;
static const double EMIT_E_SZ[3] = {0, 4, 2};
+static const double EMIT_S_SZ[3] = {4, 4, 0};
+
+static const double V[3] = {0, 0, 1};
#if 0
static const double RECV_MIN[2] = {0.125, 0.375};
@@ -48,36 +50,43 @@ static const double V[3] = {0, 0, 1};
* Helper functions
******************************************************************************/
static res_T
-realisation(void* weight, struct ssp_rng* rng, const unsigned ithread, void* ctx)
+realisation
+ (void* weights,
+ struct ssp_rng* rng,
+ const unsigned ithread,
+ void* ctx)
{
+ struct smc_doubleN_context* dblN_ctx = ctx;
+ struct sgs* sgs = dblN_ctx->integrand_data;
struct sgs_hit hit = SGS_HIT_NULL;
-
struct sgs_fragment frag = SGS_FRAGMENT_NULL;
- struct sgs* sgs = ctx;
- double dir_emit[3] = {0, 0, 0};
- double dir_spec[3] = {0, 0, 0};
- double pos_emit[3] = {0, 0, 0};
- double pos_recv[3] = {0, 0, 0};
- double pos_emit_e[3] = {0, 0, 0};
- double dir_spec_e[3] = {0, 0, 0};
+ double normal_e[3];
+ double normal_s[3];
- double range[2] = {0, DBL_MAX};
+ double dir_emit_s[3];
+ double dir_spec_e[3];
+ double dir_spec_s[3];
+ double pos_emit_e[3];
+ double pos_emit_s[3];
+ double pos_recv[3];
- double u_emit[3];
- double u_spec[3];
+ double range[2];
+
+ double u_emit_s[3];
+ double u_spec_s[3];
double u_emit_e[3];
double u_spec_e[3];
- double tmp[3];
- double beta_emit;
- double beta_spec;
- double alpha_emit;
- double alpha_spec;
- double beta_emit_e;
- double beta_spec_e;
+
double alpha_emit_e;
double alpha_spec_e;
-
+ double alpha_emit_s;
+ double alpha_spec_s;
+ double beta_emit_e;
+ double beta_spec_e;
+ double beta_emit_s;
+ double beta_spec_s;
+
double rho;
double grad_rho[3];
double I;
@@ -86,81 +95,103 @@ realisation(void* weight, struct ssp_rng* rng, const unsigned ithread, void* ctx
double S;
double w = 0;
+ double s = 0;
- enum sgs_surface_type surf_emit = SGS_SURFACE_NONE;
+ enum sgs_surface_type surf_emit_e;
+ enum sgs_surface_type surf_emit_s;
res_T res = RES_OK;
(void)ithread;
+ range[0] = 0, range[1] = INF;
+
/* Sample the sensitivity emissive surface */
sgs_geometry_sample(sgs->geom, rng, &frag);
- surf_emit = frag.surface;
+ pos_emit_s[X] = frag.position[X];
+ pos_emit_s[Y] = frag.position[Y];
+ pos_emit_s[Z] = frag.position[Z];
+ normal_s[X] = frag.normal[X];
+ normal_s[Y] = frag.normal[Y];
+ normal_s[Z] = frag.normal[Z];
+ surf_emit_s = frag.surface;
/* Sample the cosine weighted sampling of the emissive direction */
- ssp_ran_hemisphere_cos(rng, frag.normal, dir_emit, NULL);
+ ssp_ran_hemisphere_cos(rng, frag.normal, dir_emit_s, NULL);
- pos_emit[X] = frag.position[X];
- pos_emit[Y] = frag.position[Y];
- pos_emit[Z] = frag.position[Z];
-
- range[0] = 0;
- range[1] = INF;
-
- sgs_geometry_trace_ray(sgs->geom, pos_emit, dir_emit, range, surf_emit, &hit);
+ /* Trace the ray from the sensitivity emissive surface */
+ sgs_geometry_trace_ray
+ (sgs->geom, pos_emit_s, dir_emit_s, range, surf_emit_s, &hit);
if(SGS_HIT_NONE(&hit)) {
res = RES_OK;
goto error;
}
- pos_recv[X] = pos_emit[X] + hit.distance*dir_emit[X];
- pos_recv[Y] = pos_emit[Y] + hit.distance*dir_emit[Y];
- pos_recv[Z] = pos_emit[Z] + hit.distance*dir_emit[Z];
+ pos_recv[X] = pos_emit_s[X] + hit.distance*dir_emit_s[X];
+ pos_recv[Y] = pos_emit_s[Y] + hit.distance*dir_emit_s[Y];
+ pos_recv[Z] = pos_emit_s[Z] + hit.distance*dir_emit_s[Z];
- /* Does not reach the receiver */
+ /* The ray does not reach the receiver, the MC weight is NULL */
if(hit.surface != SGS_SURFACE_Z_NEG
- || RECV_MIN[X] > pos_recv[X] || pos_recv[X] > RECV_MAX[X]
- || RECV_MIN[Y] > pos_recv[Y] || pos_recv[Y] > RECV_MAX[Y]) {
+ || pos_recv[X] < RECV_MIN[X] || RECV_MAX[X] < pos_recv[X]
+ || pos_recv[Y] < RECV_MIN[Y] || RECV_MAX[Y] < pos_recv[Y]) {
goto exit;
}
- reflect(dir_spec, dir_emit, frag.normal);
- sgs_geometry_trace_ray(sgs->geom, pos_emit, dir_spec, range, surf_emit, &hit);
+ /* Find the reflected position */
+ reflect(dir_spec_s, dir_emit_s, frag.normal);
+ sgs_geometry_trace_ray
+ (sgs->geom, pos_emit_s, dir_spec_s, range, surf_emit_s, &hit);
if(SGS_HIT_NONE(&hit)) {
res = RES_OK;
goto error;
}
- pos_emit_e[X] = pos_emit[X] + hit.distance*dir_spec[X];
- pos_emit_e[Y] = pos_emit[Y] + hit.distance*dir_spec[Y];
- pos_emit_e[Z] = pos_emit[Z] + hit.distance*dir_spec[Z];
-
- if(hit.surface != SGS_SURFACE_X_POS || pos_emit_e[Z] > EMIT_THRESHOLD) {
+ pos_emit_e[X] = pos_emit_s[X] + hit.distance*dir_spec_s[X];
+ pos_emit_e[Y] = pos_emit_s[Y] + hit.distance*dir_spec_s[Y];
+ pos_emit_e[Z] = pos_emit_s[Z] + hit.distance*dir_spec_s[Z];
+ normal_e[X] = hit.normal[X];
+ normal_e[Y] = hit.normal[Y];
+ normal_e[Z] = hit.normal[Z];
+ surf_emit_e = hit.surface;
+
+ /* The reflected position is not an the emitter, the MC weight is null */
+ if(surf_emit_e != SGS_SURFACE_X_POS || pos_emit_e[Z] > EMIT_E_THRESHOLD) {
goto exit;
}
- alpha_emit = d3_dot(V, frag.normal) / d3_dot(dir_emit, frag.normal);
- alpha_spec = d3_dot(V, frag.normal) / d3_dot(dir_spec, frag.normal);
- d3_sub(u_emit, V, d3_muld(tmp, dir_emit, alpha_emit));
- d3_sub(u_spec, V, d3_muld(tmp, dir_spec, alpha_spec));
- beta_emit = d3_normalize(u_emit, u_emit);
- beta_spec = d3_normalize(u_spec, u_spec);
-
- d3_minus(dir_spec_e, dir_spec);
-
- alpha_emit_e = d3_dot(u_emit, hit.normal) / d3_dot(dir_spec_e, hit.normal);
- alpha_spec_e = d3_dot(u_spec, hit.normal) / d3_dot(dir_spec_e, hit.normal);
- d3_sub(u_emit_e, u_emit, d3_muld(tmp, dir_spec_e, alpha_emit_e));
- d3_sub(u_spec_e, u_spec, d3_muld(tmp, dir_spec_e, alpha_spec_e));
+ alpha_emit_s = d3_dot(V, normal_s) / d3_dot(dir_emit_s, normal_s);
+ alpha_spec_s = d3_dot(V, normal_s) / d3_dot(dir_spec_s, normal_s);
+ u_emit_s[X] = V[X] - alpha_emit_s*dir_emit_s[X];
+ u_emit_s[Y] = V[Y] - alpha_emit_s*dir_emit_s[Y];
+ u_emit_s[Z] = V[Z] - alpha_emit_s*dir_emit_s[Z];
+ u_spec_s[X] = V[X] - alpha_spec_s*dir_spec_s[X];
+ u_spec_s[Y] = V[Y] - alpha_spec_s*dir_spec_s[Y];
+ u_spec_s[Z] = V[Z] - alpha_spec_s*dir_spec_s[Z];
+ beta_emit_s = d3_normalize(u_emit_s, u_emit_s);
+ beta_spec_s = d3_normalize(u_spec_s, u_spec_s);
+
+ dir_spec_e[X] = -dir_spec_s[X];
+ dir_spec_e[Y] = -dir_spec_s[Y];
+ dir_spec_e[Z] = -dir_spec_s[Z];
+
+ alpha_emit_e = d3_dot(u_emit_s, normal_e) / d3_dot(dir_spec_e, normal_e);
+ alpha_spec_e = d3_dot(u_spec_s, normal_e) / d3_dot(dir_spec_e, normal_e);
+ u_emit_e[X] = u_emit_s[X] - alpha_emit_e*dir_spec_e[X];
+ u_emit_e[Y] = u_emit_s[Y] - alpha_emit_e*dir_spec_e[Y];
+ u_emit_e[Z] = u_emit_s[Z] - alpha_emit_e*dir_spec_e[Z];
+ u_spec_e[X] = u_spec_s[X] - alpha_spec_e*dir_spec_e[X];
+ u_spec_e[Y] = u_spec_s[Y] - alpha_spec_e*dir_spec_e[Y];
+ u_spec_e[Z] = u_spec_s[Z] - alpha_spec_e*dir_spec_e[Z];
beta_emit_e = d3_normalize(u_emit_e, u_emit_e);
beta_spec_e = d3_normalize(u_spec_e, u_spec_e);
rho = 0.25
- * (1 - cos(2*PI*pos_emit[X]/EMIT_SZ[X]))
- * (1 - cos(2*PI*pos_emit[Y]/EMIT_SZ[Y]));
+ * (1 - cos(2*PI*pos_emit_s[X]/EMIT_S_SZ[X]))
+ * (1 - cos(2*PI*pos_emit_s[Y]/EMIT_S_SZ[Y]));
grad_rho[X] = 0.25
- * (((2*PI)/EMIT_SZ[X])*sin(2*PI*pos_emit[X]/EMIT_SZ[X]))
- * (1 - cos(2*PI*pos_emit[Y]/EMIT_SZ[Y]));
+ * (((2*PI)/EMIT_S_SZ[X])*sin(2*PI*pos_emit_s[X]/EMIT_S_SZ[X]))
+ * (1 - cos(2*PI*pos_emit_s[Y]/EMIT_S_SZ[Y]));
grad_rho[Y] = 0.25
- * (((2*PI)/EMIT_SZ[Y])*sin(2*PI*pos_emit[Y]/EMIT_SZ[Y]))
- * (1 - cos(2*PI*pos_emit[X]/EMIT_SZ[X]));
+ * (((2*PI)/EMIT_S_SZ[Y])*sin(2*PI*pos_emit_s[Y]/EMIT_S_SZ[Y]))
+ * (1 - cos(2*PI*pos_emit_s[X]/EMIT_S_SZ[X]));
grad_rho[Z] = 0;
I =
@@ -175,15 +206,16 @@ realisation(void* weight, struct ssp_rng* rng, const unsigned ithread, void* ctx
* (1 - cos(2*PI*pos_emit_e[Y]/EMIT_E_SZ[Y]));
S =
- - beta_emit * d3_dot(grad_rho, u_emit) * I
- - rho * beta_emit * beta_emit_e * d3_dot(grad_I, u_emit_e)
- + rho * beta_spec * beta_spec_e * d3_dot(grad_I, u_spec_e);
+ - beta_emit_s * d3_dot(grad_rho, u_emit_s) * I
+ - rho * beta_emit_s * beta_emit_e * d3_dot(grad_I, u_emit_e)
+ + rho * beta_spec_s * beta_spec_e * d3_dot(grad_I, u_spec_e);
- /* w = I*EMIT_SZ[X]*EMIT_SZ[Y]*PI;*/ /* Weight */
- w = S*EMIT_SZ[X]*EMIT_SZ[Y]*PI; /* Sensib */
+ w = I*EMIT_S_SZ[X]*EMIT_S_SZ[Y]*PI; /* Weight */
+ s = S*EMIT_S_SZ[X]*EMIT_S_SZ[Y]*PI; /* Sensib */
exit:
- SMC_DOUBLE(weight) = w;
+ SMC_DOUBLEN(weights)[WEIGHT] = w;
+ SMC_DOUBLEN(weights)[SENSIB] = s;
return res;
error:
goto exit;
@@ -199,6 +231,7 @@ compute_trans_sensib(struct sgs* sgs)
struct smc_estimator_status status;
struct smc_device* smc = NULL;
struct smc_estimator* estimator = NULL;
+ struct smc_doubleN_context ctx;
res_T res = RES_OK;
ASSERT(sgs);
@@ -213,12 +246,14 @@ compute_trans_sensib(struct sgs* sgs)
/* Setup the integrator */
integrator.integrand = realisation;
- integrator.type = &smc_double;
+ integrator.type = &smc_doubleN;
integrator.max_realisations = sgs->nrealisations;
integrator.max_failures = (size_t)((double)sgs->nrealisations*0.01);
+ ctx.count = WEIGHTS_COUNT__;
+ ctx.integrand_data = sgs;
/* Run Monte-Carlo integration */
- res = smc_solve(smc, &integrator, sgs, &estimator);
+ res = smc_solve(smc, &integrator, &ctx, &estimator);
if(res != RES_OK) {
sgs_log_err(sgs, "Error while solving 4V/S -- %s.\n",
res_to_cstr(res));
@@ -234,9 +269,12 @@ compute_trans_sensib(struct sgs* sgs)
}
/* Print the result */
- sgs_log(sgs, "%g +/- %g; #failures = %lu/%lu\n",
- SMC_DOUBLE(status.E),
- SMC_DOUBLE(status.SE),
+ sgs_log(sgs, "estim ~ %g +/- %g; sensib ~ %g +/- %g\n",
+ SMC_DOUBLEN(status.E)[WEIGHT],
+ SMC_DOUBLEN(status.SE)[WEIGHT],
+ SMC_DOUBLEN(status.E)[SENSIB],
+ SMC_DOUBLEN(status.SE)[SENSIB]);
+ sgs_log(sgs, "#failures = %lu/%lu\n",
(unsigned long)status.NF,
(unsigned long)sgs->nrealisations);
