commit 665323415ddd2fb32796cb7c0e1604e2586e7cb0
parent 52687f24428cca216e2031bbe473e648e0e96c89
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 15 Sep 2021 11:45:56 +0200
Add the sgs_geometry_get_type and sgs_geometry_type_to_cstr functions
Diffstat:
7 files changed, 326 insertions(+), 290 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -55,7 +55,7 @@ set(SGS_FILES_SRC
sgs_args.c
sgs.c
sgs_compute_4v_s.c
- sgs_compute_trans_sensib.c
+ sgs_compute_sensitivity_translation.c
sgs_geometry.c
sgs_geometry_box.c
sgs_geometry_slope.c
diff --git a/src/sgs.c b/src/sgs.c
@@ -145,7 +145,7 @@ sgs_run(struct sgs* sgs)
if(res != RES_OK) goto error;
} else {
#if 1
- res = compute_trans_sensib(sgs);
+ res = compute_sensitivity_translation(sgs);
#else
res = compute_4v_s(sgs);
#endif
diff --git a/src/sgs_c.h b/src/sgs_c.h
@@ -62,7 +62,7 @@ compute_4v_s
(struct sgs* sgs);
extern LOCAL_SYM res_T
-compute_trans_sensib
+compute_sensitivity_translation
(struct sgs* sgs);
#endif /* SGS_C_H */
diff --git a/src/sgs_compute_sensitivity_translation.c b/src/sgs_compute_sensitivity_translation.c
@@ -0,0 +1,298 @@
+/* Copyright (C) 2021
+ * CNRS/RAPSODEE,
+ * CNRS/LMAP,
+ * |Meso|Star> (contact@meso-star.com),
+ * UPS/Laplace.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "sgs_c.h"
+#include "sgs_geometry.h"
+#include "sgs_log.h"
+
+#include <rsys/cstr.h>
+
+#include <star/smc.h>
+#include <star/ssp.h>
+
+/* Sugar syntax */
+enum {X, Y, Z};
+enum {WEIGHT, SENSIT, 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_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};
+static const double RECV_MAX[2] = {0.375, 0.625};
+static const double EMIT_THRESHOLD = 0.5;
+static const double V[3] = {0, 0, 1};
+#endif
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static res_T
+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;
+
+ double normal_e[3];
+ double normal_s[3];
+
+ 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 range[2];
+
+ double u_emit_s[3];
+ double u_spec_s[3];
+ double u_emit_e[3];
+ double u_spec_e[3];
+
+ 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;
+ double grad_I[3];
+
+ double S;
+
+ double w = 0;
+ double s = 0;
+
+ 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);
+ 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_s, NULL);
+
+ /* 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_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];
+
+ /* The ray does not reach the receiver, the MC weight is NULL */
+ if(hit.surface != SGS_SURFACE_Z_NEG
+ || 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;
+ }
+
+ /* 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_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_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_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_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_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 =
+ (1 - cos(2*PI*pos_emit_e[Y]/EMIT_E_SZ[Y]))
+ * (1 - cos(2*PI*pos_emit_e[Z]/EMIT_E_SZ[Z]));
+ grad_I[X] = 0;
+ grad_I[Y] =
+ (((2*PI)/EMIT_E_SZ[Y])*sin(2*PI*pos_emit_e[Y]/EMIT_E_SZ[Y]))
+ * (1 - cos(2*PI*pos_emit_e[Z]/EMIT_E_SZ[Z]));
+ grad_I[Z] =
+ (((2*PI)/EMIT_E_SZ[Z])*sin(2*PI*pos_emit_e[Z]/EMIT_E_SZ[Z]))
+ * (1 - cos(2*PI*pos_emit_e[Y]/EMIT_E_SZ[Y]));
+
+ S =
+ - 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_S_SZ[X]*EMIT_S_SZ[Y]*PI; /* Weight */
+ s = S*EMIT_S_SZ[X]*EMIT_S_SZ[Y]*PI; /* Sensib */
+
+exit:
+ SMC_DOUBLEN(weights)[WEIGHT] = w;
+ SMC_DOUBLEN(weights)[SENSIT] = s;
+ return res;
+error:
+ goto exit;
+}
+
+/*******************************************************************************
+ * Local function
+ ******************************************************************************/
+res_T
+compute_sensitivity_translation(struct sgs* sgs)
+{
+ struct smc_integrator integrator;
+ struct smc_estimator_status status;
+ struct smc_device* smc = NULL;
+ struct smc_estimator* estimator = NULL;
+ struct smc_doubleN_context ctx;
+ enum sgs_geometry_type type;
+ res_T res = RES_OK;
+ ASSERT(sgs);
+
+ type = sgs_geometry_get_type(sgs->geom);
+ if(type != SGS_GEOMETRY_BOX) {
+ sgs_log_err(sgs,
+ "Invalid `%s' geometry. The box geometry is the only one that is "
+ "currently supported by the translation sensitivity computation.\n",
+ sgs_geometry_type_to_cstr(type));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Create the Star-MonteCarlo device */
+ res = smc_device_create
+ (&sgs->logger, sgs->allocator, sgs->nthreads, &ssp_rng_mt19937_64, &smc);
+ if(res != RES_OK) {
+ sgs_log_err(sgs, "Could not create the Star-MonteCarlo device -- %s.\n",
+ res_to_cstr(res));
+ goto error;
+ }
+
+ /* Setup the integrator */
+ integrator.integrand = realisation;
+ 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, &ctx, &estimator);
+ if(res != RES_OK) {
+ sgs_log_err(sgs, "Error while solving 4V/S -- %s.\n",
+ res_to_cstr(res));
+ goto error;
+ }
+
+ /* Retrieve the estimated value */
+ SMC(estimator_get_status(estimator, &status));
+ if(status.NF >= integrator.max_failures) {
+ sgs_log_err(sgs, "Too many failures: %lu\n", (unsigned long)status.NF);
+ res = RES_BAD_OP;
+ goto error;
+ }
+
+ /* Print the result */
+ sgs_log(sgs, "estim ~ %g +/- %g; sensib ~ %g +/- %g\n",
+ SMC_DOUBLEN(status.E)[WEIGHT],
+ SMC_DOUBLEN(status.SE)[WEIGHT],
+ SMC_DOUBLEN(status.E)[SENSIT],
+ SMC_DOUBLEN(status.SE)[SENSIT]);
+ sgs_log(sgs, "#failures = %lu/%lu\n",
+ (unsigned long)status.NF,
+ (unsigned long)sgs->nrealisations);
+
+exit:
+ if(smc) SMC(device_ref_put(smc));
+ if(estimator) SMC(estimator_ref_put(estimator));
+ return res;
+error:
+ goto exit;
+}
diff --git a/src/sgs_compute_trans_sensib.c b/src/sgs_compute_trans_sensib.c
@@ -1,287 +0,0 @@
-/* Copyright (C) 2021
- * CNRS/RAPSODEE,
- * CNRS/LMAP,
- * |Meso|Star> (contact@meso-star.com),
- * UPS/Laplace.
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>. */
-
-#include "sgs_c.h"
-#include "sgs_geometry.h"
-#include "sgs_log.h"
-
-#include <rsys/cstr.h>
-
-#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_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};
-static const double RECV_MAX[2] = {0.375, 0.625};
-static const double EMIT_THRESHOLD = 0.5;
-static const double V[3] = {0, 0, 1};
-#endif
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-static res_T
-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;
-
- double normal_e[3];
- double normal_s[3];
-
- 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 range[2];
-
- double u_emit_s[3];
- double u_spec_s[3];
- double u_emit_e[3];
- double u_spec_e[3];
-
- 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;
- double grad_I[3];
-
- double S;
-
- double w = 0;
- double s = 0;
-
- 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);
- 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_s, NULL);
-
- /* 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_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];
-
- /* The ray does not reach the receiver, the MC weight is NULL */
- if(hit.surface != SGS_SURFACE_Z_NEG
- || 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;
- }
-
- /* 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_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_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_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_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_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 =
- (1 - cos(2*PI*pos_emit_e[Y]/EMIT_E_SZ[Y]))
- * (1 - cos(2*PI*pos_emit_e[Z]/EMIT_E_SZ[Z]));
- grad_I[X] = 0;
- grad_I[Y] =
- (((2*PI)/EMIT_E_SZ[Y])*sin(2*PI*pos_emit_e[Y]/EMIT_E_SZ[Y]))
- * (1 - cos(2*PI*pos_emit_e[Z]/EMIT_E_SZ[Z]));
- grad_I[Z] =
- (((2*PI)/EMIT_E_SZ[Z])*sin(2*PI*pos_emit_e[Z]/EMIT_E_SZ[Z]))
- * (1 - cos(2*PI*pos_emit_e[Y]/EMIT_E_SZ[Y]));
-
- S =
- - 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_S_SZ[X]*EMIT_S_SZ[Y]*PI; /* Weight */
- s = S*EMIT_S_SZ[X]*EMIT_S_SZ[Y]*PI; /* Sensib */
-
-exit:
- SMC_DOUBLEN(weights)[WEIGHT] = w;
- SMC_DOUBLEN(weights)[SENSIB] = s;
- return res;
-error:
- goto exit;
-}
-
-/*******************************************************************************
- * Local function
- ******************************************************************************/
-res_T
-compute_trans_sensib(struct sgs* sgs)
-{
- struct smc_integrator integrator;
- 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);
-
- /* Create the Star-MonteCarlo device */
- res = smc_device_create
- (&sgs->logger, sgs->allocator, sgs->nthreads, &ssp_rng_mt19937_64, &smc);
- if(res != RES_OK) {
- sgs_log_err(sgs, "Could not create the Star-MonteCarlo device -- %s.\n",
- res_to_cstr(res));
- goto error;
- }
-
- /* Setup the integrator */
- integrator.integrand = realisation;
- 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, &ctx, &estimator);
- if(res != RES_OK) {
- sgs_log_err(sgs, "Error while solving 4V/S -- %s.\n",
- res_to_cstr(res));
- goto error;
- }
-
- /* Retrieve the estimated value */
- SMC(estimator_get_status(estimator, &status));
- if(status.NF >= integrator.max_failures) {
- sgs_log_err(sgs, "Too many failures: %lu\n", (unsigned long)status.NF);
- res = RES_BAD_OP;
- goto error;
- }
-
- /* Print the result */
- 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);
-
-exit:
- if(smc) SMC(device_ref_put(smc));
- if(estimator) SMC(estimator_ref_put(estimator));
- return res;
-error:
- goto exit;
-}
diff --git a/src/sgs_geometry.c b/src/sgs_geometry.c
@@ -182,6 +182,13 @@ sgs_geometry_ref_put(struct sgs_geometry* geom)
ref_put(&geom->ref, release_geometry);
}
+enum sgs_geometry_type
+sgs_geometry_get_type(const struct sgs_geometry* geom)
+{
+ ASSERT(geom);
+ return geom->type;
+}
+
int
sgs_geometry_get_sampling_mask(const struct sgs_geometry* geom)
{
diff --git a/src/sgs_geometry.h b/src/sgs_geometry.h
@@ -143,6 +143,20 @@ struct sgs;
struct sgs_geometry;
struct ssp_rng;
+static INLINE const char*
+sgs_geometry_type_to_cstr(const enum sgs_geometry_type type)
+{
+ const char* str = "unknown";
+ switch(type) {
+ case SGS_GEOMETRY_BOX: str = "box"; break;
+ case SGS_GEOMETRY_SLOPE: str = "slope"; break;
+ case SGS_GEOMETRY_STEP: str = "step"; break;
+ case SGS_GEOMETRY_NONE: str = "none"; break;
+ default: FATAL("Unreachable code.\n"); break;
+ }
+ return str;
+}
+
/*******************************************************************************
* Geometry API
******************************************************************************/
@@ -172,6 +186,10 @@ extern LOCAL_SYM void
sgs_geometry_ref_put
(struct sgs_geometry* geom);
+extern LOCAL_SYM enum sgs_geometry_type
+sgs_geometry_get_type
+ (const struct sgs_geometry* geom);
+
extern LOCAL_SYM int
sgs_geometry_get_sampling_mask
(const struct sgs_geometry* geom);
