commit 60b0db63ea3bc9d669eb28e3312c20fdc9d7cba7
parent 83bcb658fe11bb59afa6f4572067fceb185362b7
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 11 Mar 2024 15:09:42 +0100
Tabulation update for function H
From now on, the smaller the x, the smaller the tabulation step. This
ensures that the function is finely tabulated where estimation of H is
difficult (i.e. close to 0).
The test is updated to check the relative error on x over the whole
range. We check that the relative error is less than 1e-5 or 1e-6,
except for Hx values outside [1e-9, 1-1e-9].
Diffstat:
| M | src/swf.h | | | 12 | ++++++++++-- |
| M | src/swf_H.c | | | 58 | ++++++++++++++++++++++++++++++---------------------------- |
| M | src/test_swf_H3d.c | | | 66 | +++++++++++++++++++++++++++++++++--------------------------------- |
3 files changed, 73 insertions(+), 63 deletions(-)
diff --git a/src/swf.h b/src/swf.h
@@ -42,10 +42,18 @@ struct mem_allocator;
struct swf_H_tabulate_args {
double x_min;
double x_max;
- double delta_x;
+
+ /* The step member variable is used to calculate the x arguments to be
+ * tabulated. Each x to be tabulated is calculated by adding its product with
+ * the step value to the previous one. The delta x is therefore relative to
+ * each x value: the smaller the x, the smaller the tabulation step. And this
+ * is what we're looking for, since the function is more difficult to estimate
+ * near 0 */
+ double step; /* x_next = x + x*step: */
+
struct mem_allocator* allocator; /* NULL <=> use default allocator */
};
-#define SWF_H_TABULATE_ARGS_DEFAULT__ {7.5e-3, 3.8, 1e-5, NULL}
+#define SWF_H_TABULATE_ARGS_DEFAULT__ {7.5e-3, 3.8, 1e-3, NULL}
static const struct swf_H_tabulate_args SWF_H_TABULATE_ARGS_DEFAULT =
SWF_H_TABULATE_ARGS_DEFAULT__;
diff --git a/src/swf_H.c b/src/swf_H.c
@@ -38,11 +38,11 @@ check_swf_H_tabulate_args(const struct swf_H_tabulate_args* args)
return RES_BAD_ARG;
/* Delta X cannot be null */
- if(args->delta_x <= 0)
+ if(args->step <= 0)
return RES_BAD_ARG;
/* Delta X cannot be greater than the X range */
- if(args->delta_x > args->x_max - args->x_min)
+ if(args->step > args->x_max - args->x_min)
return RES_BAD_ARG;
return RES_OK;
@@ -163,8 +163,7 @@ swf_H3d_tabulate
struct item* items = NULL;
size_t nitems = 0;
size_t i = 0;
- size_t j = 0;
- double x_range = 0;
+ double x = 0;
double norm = 0;
res_T res = RES_OK;
@@ -173,36 +172,39 @@ swf_H3d_tabulate
if((res = check_swf_H_tabulate_args(args)) != RES_OK) goto error;
if((res = tabulation_create(args->allocator, &tab)) != RES_OK) goto error;
- /* Calculate the maximum number of arguments x to be tabulated. Note that
- * there is one more argument than the number of intervals, hence the +1 */
- x_range = args->x_max - args->x_min;
- nitems = (size_t)ceil(x_range / args->delta_x) + 1;
-
- /* Tabulation memory space allocation */
- if((res = darray_item_resize(&tab->items, nitems)) != RES_OK) goto error;
- items = darray_item_data_get(&tab->items);
-
/* Setup the tabulation */
- FOR_EACH(i, 0, nitems) {
+ i = 0;
+ for(x = args->x_min; x < args->x_max; x += x*args->step) {
+ const struct item* pitem = darray_item_cdata_get(&tab->items) + i - 1;
+ struct item item = ITEM_NULL;
+
/* Evaluate H(x) and its derivatives up to order 3 */
- items[j].x = MMIN(args->x_min + (double)i*args->delta_x, args->x_max);
- items[j].fx = H3d(items[j].x, &items[j].dfx, &items[j].d2fx, &items[j].d3fx);
+ item.x = MMIN(x, args->x_max);
+ item.fx = H3d(item.x, &item.dfx, &item.d2fx, &item.d3fx);
+
+ if(i > 0) {
+ /* Detect oscillations due to numerical problems.
+ * These oscillations occur close to 0 */
+ if(item.fx < pitem->fx) {
+ res = RES_BAD_OP;
+ goto error;
+ }
+ /* Do not tabulate argument whose value is (numerically) equal to the
+ * previous one: it would artificially add staircase steps */
+ if(item.fx == pitem->fx) {
+ continue;
+ }
+ }
- if(j == 0) {
- ++j;
- } else {
- /* Detect oscillations due to numerical problems.
- * These oscillations occur close to 0 */
- if(items[j].fx < items[j-1].fx) { res = RES_BAD_OP; goto error; }
+ ++i;
- /* Reject tabulation whose value is (numerically) equal to the previous
- * one: it would artificially add staircase steps */
- j += items[j].fx > items[j-1].fx;
- }
+ res = darray_item_push_back(&tab->items, &item);
+ if(res != RES_OK) goto error;
}
- CHK(j && darray_item_resize(&tab->items, j) == RES_OK);
- nitems = j;
+ ASSERT(darray_item_size_get(&tab->items) == i);
+ items = darray_item_data_get(&tab->items);
+ nitems = darray_item_size_get(&tab->items);
/* Normalize the tabulation */
norm = items[nitems - 1].fx;
diff --git a/src/test_swf_H3d.c b/src/test_swf_H3d.c
@@ -24,12 +24,6 @@
/*******************************************************************************
* Helper functions
******************************************************************************/
-static double
-rand_canonic(void)
-{
- return (double)rand() / ((double)RAND_MAX+1);
-}
-
static INLINE void
check_tabulation_creation(void)
{
@@ -55,13 +49,13 @@ check_tabulation_creation(void)
args.x_min = SWF_H_TABULATE_ARGS_DEFAULT.x_min;
args.x_max = SWF_H_TABULATE_ARGS_DEFAULT.x_max;
- args.delta_x = args.x_max - args.x_min + 1.e-3;
+ args.step = args.x_max - args.x_min + 1.e-3;
CHK(swf_H3d_tabulate(&args, &tab) == RES_BAD_ARG);
- args.delta_x = 0;
+ args.step = 0;
CHK(swf_H3d_tabulate(&args, &tab) == RES_BAD_ARG);
- args.delta_x = 1.0e-3;
+ args.step = 1.0e-2;
CHK(swf_H3d_tabulate(&args, &tab) == RES_OK);
CHK(swf_tabulation_ref_put(tab) == RES_OK);
}
@@ -69,37 +63,43 @@ check_tabulation_creation(void)
static void
check_tabulation_inversion(void)
{
- const double delta_Hx = 1.0e-6;
+ const double nsteps = 10;
struct swf_H_tabulate_args args = SWF_H_TABULATE_ARGS_DEFAULT;
struct swf_tabulation* tab = NULL;
- double errl = 0;
- double errq = 0;
- size_t n = 0;
+ double pHx = 0;
+ double x0 = 0;
size_t i = 0;
CHK(swf_H3d_tabulate(&args, &tab) == RES_OK);
- /*errl = swf_tabulation_max_x_error(tab, SWF_LINEAR);
- errq = swf_tabulation_max_x_error(tab, SWF_QUADRATIC);
- CHK(errq < errl);*/
-
- n = (size_t)ceil(1.0/delta_Hx);
-
- FOR_EACH(i, 0, n) {
- const double Hx = (double)i*delta_Hx + rand_canonic()*delta_Hx;
- const double x_ref = swf_H3d_inverse(Hx);
- const double xl = swf_tabulation_inverse(tab, SWF_LINEAR, Hx);
- const double xq = swf_tabulation_inverse(tab, SWF_QUADRATIC, Hx);
-
- if(!x_ref) {
- CHK(x_ref == xl);
- } else {
- errl = fabs((x_ref - xl) / x_ref);
- errq = fabs((x_ref - xq) / x_ref);
- /*printf("%e %e %e %e %e\n", x_ref, xq, xl, errq, errl);*/
- CHK(errl < 1.e-4);
- CHK(errq < 1.e-6);
+ for(x0 = args.x_min; x0 < args.x_max; x0 += x0*args.step) {
+ const double x1 = x0 + x0*args.step;
+ const double delta_x = (x1 - x0) / (double)nsteps;
+
+ for(i = 1; i < nsteps; ++i) {
+ const double x = x0 + (double)i*delta_x;
+ const double Hx = swf_H3d_eval(x);
+ const double xl = swf_tabulation_inverse(tab, SWF_LINEAR, Hx);
+ const double xq = swf_tabulation_inverse(tab, SWF_QUADRATIC, Hx);
+ const double errl = fabs((x - xl) / x);
+ const double errq = fabs((x - xq) / x);
+
+ /* Do not check the value of x if the corresponding H is indistinguishable
+ * from an H calculated from a previous x. In other words, the H's are
+ * numerically equal and any of the corresponding x-values is a valid
+ * inversion result, numerically speaking. */
+ if(Hx == pHx) continue;
+
+ printf("%e %.20e %e %e\n", x, Hx, errq, errl);
+ if(1e-9 < Hx && Hx < (1.0 - 1e-9)) {
+ CHK(errl < 1.e-5);
+ CHK(errq < 1.e-6);
+ } else {
+ CHK(errl < 1.e-1);
+ CHK(errq < 1.e-1);
+ }
+ pHx = Hx;
}
}
