commit a3377a8c8c785b7d1f5cebb286956bea98f380ca
parent d97612754c9ccd7592792a5d1d502b805778b7e2
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 1 Feb 2023 18:32:58 +0100
Remplace les sources C par leur version Noweb
Diffstat:
3 files changed, 14 insertions(+), 425 deletions(-)
diff --git a/Makefile b/Makefile
@@ -24,16 +24,19 @@
include config.mk
+DOC = sgs_compute_sensitivity_translation
+
+TMPC =\
+ src/sgs_compute_sensitivity_translation.c
+
SRC =\
src/sgs_args.c\
src/sgs.c\
- src/sgs_compute_sensitivity_translation.c\
src/sgs_geometry_box.c\
src/sgs_geometry.c\
src/sgs_log.c\
- src/sgs_main.c
-
-DOC = sgs_compute_sensitivity_translation
+ src/sgs_main.c\
+ $(TMPC)
NOWEB=\
src/$(DOC).nw
@@ -58,9 +61,12 @@ sgs: $(OBJ)
@echo "configure $(MAKE)"
@$(SHELL) configure.sh program && echo "config done" > $@ || exit 1
+$(TMPC): $(NOWEB) config.mk
+ @$(TANGLE) $(TANGLE_OPTS)\
+ -R$$(c="$@" && echo "$${c##*/}" | sed 's/_/\\_/g') $(NOWEB) | cpif $@
+
.SUFFIXES: .c .d .o
.c.d:
- @#echo "DEP $@"
@$(CC) $(CFLAGS) -MM -MT "$(@:.d=.o) $@" $< -MF $@
.c.o:
@@ -102,8 +108,8 @@ uninstall:
rm -f $(DESTDIR)$(PREFIX)/share/doc/star-gs/README.md
clean:
- @rm -f $(OBJ) $(TEX) sgs .config_program .config_document ./*.aux ./*.bbl\
- ./*.idx ./*.log ./*.dvi ./*.glo ./*.lof ./*.toc ./*.out ./*.blg\
+ @rm -f $(OBJ) $(TEX) $(TMPC) sgs .config_program .config_document ./*.aux\
+ ./*.bbl ./*.idx ./*.log ./*.dvi ./*.glo ./*.lof ./*.toc ./*.out ./*.blg\
./*.nav ./*.snm ./*.vrb weave* $(DOC).pdf
distclean: clean
diff --git a/config.mk b/config.mk
@@ -33,7 +33,7 @@ LATEX = pdflatex
TANGLE = notangle
WEAVE = noweave
-TANGLE_OPTS = #-L
+TANGLE_OPTS = -L
WEAVE_OPTS = -index -delay -v
################################################################################
diff --git a/src/sgs_compute_sensitivity_translation.c b/src/sgs_compute_sensitivity_translation.c
@@ -1,417 +0,0 @@
-/* Copyright (C) 2021-2023 Centre National de la Recherche Scientifique
- * Copyright (C) 2021-2023 Institut Mines Télécom Albi-Carmaux
- * Copyright (C) 2021-2023 |Méso|Star> (contact@meso-star.com)
- * Copyright (C) 2021-2023 Université Clermont Auvergne
- * Copyright (C) 2021-2023 Université de Lorraine
- * Copyright (C) 2021-2023 Université de Lyon
- * Copyright (C) 2021-2023 Université de Toulouse
- *
- * 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 */
-#if 0
-/* Il semblerait que le cube soit de taille {4, 4, 3} */
-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};
-#else
-static const double RECV_MIN[2] = {0.125/*<=>1/8*/, 0.375/*<=>3/8*/};
-static const double RECV_MAX[2] = {0.375/*<=>3/8*/, 0.625/*<=>5/8*/};
-static const double EMIT_E_THRESHOLD = 2.0/3.0;
-static const double EMIT_E_SZ[3] = {0, 1, 2.0/3.0};
-static const double EMIT_S_SZ[3] = {1, 1, 0};
-#endif
-
-/* Vecteur de déformation (à renommer en gamma ?) */
-static const double V[3] = {0, 0, 1};
-
-struct deformation_2d {
- double alpha;
- double beta;
- double u[3];
-};
-
-struct scene {
- /* Upper and lower limit of the scene */
- double box_low[3];
- double box_upp[3];
-
- double box_sz[3]; /* Size of the scene */
-
- /* Geometry of the receiver on the box bottom plane */
- double recv_min[3];
- double recv_max[3];
-
- double emit_e_threshold; /* Geometric upper bound of the radiative source */
- double emit_e_sz[3]; /* Size of the radiative source */
- double emit_s_sz[3]; /* Size of the specular surface */
-};
-
-/*******************************************************************************
- * Helper functions
- ******************************************************************************/
-static void
-setup_scene(const struct sgs_geometry* geometry, struct scene* scn)
-{
- ASSERT(geometry && scn);
-
- /* Compute the size of the geometry */
- sgs_geometry_get_aabb(geometry, scn->box_low, scn->box_upp);
- d3_sub(scn->box_sz, scn->box_upp, scn->box_low);
-
- /* Setup the geometry of the receiver */
- scn->recv_min[X] = RECV_MIN[X]*scn->box_sz[X] + scn->box_low[X];
- scn->recv_min[Y] = RECV_MIN[Y]*scn->box_sz[Y] + scn->box_low[Y];
- scn->recv_min[Z] = 0;
- scn->recv_max[X] = RECV_MAX[X]*scn->box_sz[X] + scn->box_low[X];
- scn->recv_max[Y] = RECV_MAX[Y]*scn->box_sz[Y] + scn->box_low[Y];
- scn->recv_max[Z] = 0;
-
- /* Setup the geometry of the radiative source */
- scn->emit_e_threshold = EMIT_E_THRESHOLD*scn->box_sz[Z] + scn->box_low[Z];
-
- /* Compute the size of the radiative source/specular surface */
- scn->emit_e_sz[X] = EMIT_E_SZ[X]*scn->box_sz[X] + scn->box_low[X];
- scn->emit_e_sz[Y] = EMIT_E_SZ[Y]*scn->box_sz[Y] + scn->box_low[Y];
- scn->emit_e_sz[Z] = EMIT_E_SZ[Z]*scn->box_sz[Z] + scn->box_low[Z];
- scn->emit_s_sz[X] = EMIT_S_SZ[X]*scn->box_sz[X] + scn->box_low[X];
- scn->emit_s_sz[Y] = EMIT_S_SZ[Y]*scn->box_sz[Y] + scn->box_low[Y];
- scn->emit_s_sz[Z] = EMIT_S_SZ[Z]*scn->box_sz[Z] + scn->box_low[Z];
-}
-
-/* Return 1 if the `hit' lies on the receiver and 0 otherwise */
-static int
-hit_receiver
- (const struct scene* scn,
- const double ray_org[3],
- const double ray_dir[3],
- const struct sgs_hit* hit,
- double pos_recv[3])
-{
- ASSERT(scn && ray_org && ray_dir && hit && pos_recv);
-
- pos_recv[X] = ray_org[X] + hit->distance*ray_dir[X];
- pos_recv[Y] = ray_org[Y] + hit->distance*ray_dir[Y];
- pos_recv[Z] = ray_org[Z] + hit->distance*ray_dir[Z];
-
- if(hit->surface != SGS_SURFACE_Z_MIN
- || pos_recv[X] < scn->recv_min[X] || scn->recv_max[X] < pos_recv[X]
- || pos_recv[Y] < scn->recv_min[Y] || scn->recv_max[Y] < pos_recv[Y]) {
- return 0; /* The ray does not intersect the receiver */
- } else {
- return 1;
- }
-}
-
-/* Return 1 if the `hit' lies on the source and 0 otherwise */
-static int
-hit_source
- (const struct scene* scn,
- const double ray_org[3],
- const double ray_dir[3],
- const struct sgs_hit* hit,
- double pos_emit_e[3])
-{
- pos_emit_e[X] = ray_org[X] + hit->distance*ray_dir[X];
- pos_emit_e[Y] = ray_org[Y] + hit->distance*ray_dir[Y];
- pos_emit_e[Z] = ray_org[Z] + hit->distance*ray_dir[Z];
-
- if(hit->surface != SGS_SURFACE_X_MAX
- || pos_emit_e[Z] > scn->emit_e_threshold) {
- return 0; /* The ray does not intersect the emitter */
- } else {
- return 1;
- }
-}
-
-static double
-compute_rho
- (const double position[2], /* Position sur la surface réfléchissante */
- const double size[2]) /* Taille de la surface réfléchissante */
-{
- ASSERT(position && size);
- return 0.25
- * (1 - cos(2*PI*position[X]/size[X]))
- * (1 - cos(2*PI*position[Y]/size[Y]));
-}
-
-static void
-compute_grad_rho
- (const double position[2], /* Position sur la surface réfléchissante */
- const double size[2], /* Taille de la surface réfléchissante */
- double out_gradient[2])
-{
- ASSERT(position && size && out_gradient);
- out_gradient[X] = 0.25
- * (((2*PI)/size[X])*sin(2*PI*position[X]/size[X]))
- * (1 - cos(2*PI*position[Y]/size[Y]));
- out_gradient[Y] = 0.25
- * (((2*PI)/size[Y])*sin(2*PI*position[Y]/size[Y]))
- * (1 - cos(2*PI*position[X]/size[X]));
-}
-
-static double
-compute_I
- (const double position[2], /* Position sur la surface émettrice */
- const double size[2]) /* Taille de la surface émettrice */
-{
- ASSERT(position && size);
- return
- (1 - cos(2*PI*position[X]/size[X]))
- * (1 - cos(2*PI*position[Y]/size[Y]));
-}
-
-static void
-compute_grad_I
- (const double position[2], /* Position sur la surface réfléchissante */
- const double size[2], /* Taille de la surface réfléchissante */
- double out_gradient[2])
-{
- ASSERT(position && size && out_gradient);
- out_gradient[X] =
- (((2*PI)/size[X])*sin(2*PI*position[X]/size[X]))
- * (1 - cos(2*PI*position[Y]/size[Y]));
- out_gradient[Y] =
- (((2*PI)/size[Y])*sin(2*PI*position[Y]/size[Y]))
- * (1 - cos(2*PI*position[X]/size[X]));
-}
-
-static void
-decomposition
- (const double chi[3],
- const double normal[3],
- const double omega[3],
- struct deformation_2d* out_vector)
-{
- ASSERT(chi && normal && omega && out_vector);
- ASSERT(d3_is_normalized(normal));
- ASSERT(d3_is_normalized(omega));
-
- out_vector->alpha = d3_dot(chi, normal) / d3_dot(omega, normal);
- out_vector->u[X] = chi[X] - out_vector->alpha*omega[X];
- out_vector->u[Y] = chi[Y] - out_vector->alpha*omega[Y];
- out_vector->u[Z] = chi[Z] - out_vector->alpha*omega[Z];
- out_vector->beta = d3_normalize(out_vector->u, out_vector->u);
-}
-
-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;
-
- struct deformation_2d proj_u_e;
- struct deformation_2d proj_V_e;
- struct deformation_2d proj_V_s;
-
- struct scene scn;
- double emit_e_sz_2d[2];
- double emit_s_sz_2d[2];
-
- 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], pos_emit_e_2d[2];
- double pos_emit_s[3], pos_emit_s_2d[2];
- double pos_recv[3];
-
- double rho;
- double grad_rho[3], grad_rho_2d[2];
- double I;
- double grad_I[3], grad_I_2d[2];
-
- double S;
-
- double w = 0;
- double s = 0;
-
- enum sgs_surface_type surf_emit_s;
-
- res_T res = RES_OK;
- (void)ithread;
-
- setup_scene(sgs->scene.geom, &scn);
-
- /* Sample the sensitivity emissive surface */
- sgs_geometry_sample(sgs->scene.geom, rng, &frag);
- d3_set(pos_emit_s, frag.position);
- d3_set(normal_s, frag.normal);
- surf_emit_s = frag.surface;
-
- /* Sample the cosine weighted sampling of the emissive direction */
- ssp_ran_hemisphere_cos(rng, normal_s, dir_emit_s, NULL);
-
- /* Helper macro used as syntactic sugar */
- #define TRACE_RAY(Org, Dir, StartFrom, Hit) { \
- const double range[2] = {0, DBL_MAX}; \
- sgs_geometry_trace_ray(sgs->scene.geom, (Org), (Dir), range, (StartFrom), (Hit));\
- if(SGS_HIT_NONE(&hit)) { res = RES_BAD_ARG; goto error; } \
- } (void)0
-
- /* Trace the sampled ray and check that if reaches the receiver? */
- TRACE_RAY(pos_emit_s, dir_emit_s, surf_emit_s, &hit);
- if(!hit_receiver(&scn, pos_emit_s, dir_emit_s, &hit, pos_recv)) goto exit;
-
- /* Trace the reflected ray and check that if reaches the source? */
- reflect(dir_spec_s, dir_emit_s, normal_s);
- TRACE_RAY(pos_emit_s, dir_spec_s, surf_emit_s, &hit);
- if(!hit_source(&scn, pos_emit_s, dir_spec_s, &hit, pos_emit_e)) goto exit;
-
- #undef TRACE_RAY
-
- /* Calcul du poids */
-
- d3_set(normal_e, hit.normal);
- d3_minus(dir_spec_e, dir_spec_s);
-
- decomposition(V, normal_e, dir_spec_e, &proj_V_e);
- decomposition(V, normal_s, dir_emit_s, &proj_V_s);
- decomposition(proj_V_s.u, normal_e, dir_spec_e, &proj_u_e);
-
- pos_emit_s_2d[X] = pos_emit_s[X];
- pos_emit_s_2d[Y] = pos_emit_s[Y];
- emit_s_sz_2d[X] = scn.emit_s_sz[X];
- emit_s_sz_2d[Y] = scn.emit_s_sz[Y];
- rho = compute_rho(pos_emit_s_2d, emit_s_sz_2d);
- compute_grad_rho(pos_emit_s_2d, emit_s_sz_2d, grad_rho_2d);
- grad_rho[X] = grad_rho_2d[X];
- grad_rho[Y] = grad_rho_2d[Y];
- grad_rho[Z] = 0;
-
- pos_emit_e_2d[X] = pos_emit_e[Y];
- pos_emit_e_2d[Y] = pos_emit_e[Z];
- emit_e_sz_2d[X] = scn.emit_e_sz[Y];
- emit_e_sz_2d[Y] = scn.emit_e_sz[Z];
- I = compute_I(pos_emit_e_2d, emit_e_sz_2d);
- compute_grad_I(pos_emit_e_2d, emit_e_sz_2d, grad_I_2d);
- grad_I[X] = 0;
- grad_I[Y] = grad_I_2d[X];
- grad_I[Z] = grad_I_2d[Y];
-
- S =
- - proj_V_s.beta * d3_dot(grad_rho, proj_V_s.u) * I
- - rho * proj_V_s.beta * proj_u_e.beta * d3_dot(grad_I, proj_u_e.u)
- + rho * proj_V_e.beta * d3_dot(grad_I, proj_V_e.u);
-
- w = I*scn.emit_s_sz[X]*scn.emit_s_sz[Y]*PI; /* Weight */
- s = S*scn.emit_s_sz[X]*scn.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_sensitivity_translation(struct sgs* sgs)
-{
- struct smc_device_create_args args = SMC_DEVICE_CREATE_ARGS_DEFAULT;
- struct smc_estimator_status status = SMC_ESTIMATOR_STATUS_NULL;
- struct smc_integrator integrator = SMC_INTEGRATOR_NULL;
- struct smc_device* smc = NULL;
- struct smc_estimator* estimator = NULL;
- struct smc_doubleN_context ctx;
- res_T res = RES_OK;
- ASSERT(sgs);
-
- if(sgs_geometry_get_sampling_mask(sgs->scene.geom) != SGS_SURFACE_Z_MAX) {
- sgs_log_err(sgs,
- "Invalid sensitivity emissive surface. Only positive Z face of the box "
- "geometry can be a sensitivity source.\n");
- res = RES_BAD_ARG;
- goto error;
- }
-
- /* Create the Star-MonteCarlo device */
- args.logger = &sgs->logger;
- args.allocator = sgs->allocator;
- args.nthreads_hint = sgs->nthreads;
- args.rng_type = SSP_RNG_MT19937_64;
- res = smc_device_create(&args, &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;
-}
