commit c465bdbb86537b914653ccaa7d7073179f4c5c7b
parent 0c4e3e0eb13663e6fcf47211c27096c176701e43
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Sun, 10 Mar 2024 11:11:18 +0100
Update H table
In addition to the function value, record the first, second and third
derivatives of the tabulated x-value for later use in interpolation.
Return an error when a numerical uncertainty is detected on the
evaluation of H, i.e. when the function is not increasing. In addition,
do not register x arguments whose value is equal to the previous
tabulated x, in order to avoid stair stepping. Thus, if successful,
tabulation guarantees that H(x) is a strictly increasing normalized
function.
Diffstat:
3 files changed, 55 insertions(+), 17 deletions(-)
diff --git a/src/swf_H.c b/src/swf_H.c
@@ -163,6 +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 norm = 0;
res_T res = RES_OK;
@@ -172,8 +173,8 @@ 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 number of arguments x to be tabulated. Note that there is one
- * more argument than the number of intervals, hence the +1 */
+ /* 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;
@@ -183,14 +184,29 @@ swf_H3d_tabulate
/* Setup the tabulation */
FOR_EACH(i, 0, nitems) {
- items[i].x = MMIN(args->x_min + (double)i*args->delta_x, args->x_max);
- items[i].f_x = swf_H3d_eval(items[i].x);
+ /* 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);
+
+ 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; }
+
+ /* 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;
+ }
}
- ASSERT(items[nitems-1].x == args->x_max);
+
+ CHK(j && darray_item_resize(&tab->items, j) == RES_OK);
+ nitems = j;
/* Normalize the tabulation */
- norm = items[nitems - 1].f_x;
- FOR_EACH(i, 0, nitems) items[i].f_x /= norm;
+ norm = items[nitems - 1].fx;
+ FOR_EACH(i, 0, nitems) items[i].fx /= norm;
exit:
if(out_tab) *out_tab = tab;
diff --git a/src/swf_tabulation.c b/src/swf_tabulation.c
@@ -32,13 +32,26 @@ release_tabulation(ref_T* ref)
}
static FINLINE int
-cmp_item(const void* a, const void* b)
+cmp_fx(const void* a, const void* b)
{
const double f_a = *((const double*)a);
- const double f_b = ((const struct item*)b)->f_x;
+ const double f_b = ((const struct item*)b)->fx;
return f_a < f_b ? -1 : (f_a > f_b ? 1 : 0);
}
+static INLINE double
+linear_interpolation
+ (const double y,
+ const struct item* a,
+ const struct item* b)
+{
+ double u, x;
+ ASSERT(a && b);
+ u = (y - a->fx) / (b->fx - a->fx);
+ x = u*(b->x - a->x) + a->x;
+ return x;
+}
+
/*******************************************************************************
* Exported symbols
******************************************************************************/
@@ -63,7 +76,6 @@ swf_tabulation_inverse(const struct swf_tabulation* tab, const double y)
{
const struct item* items = NULL;
const struct item* find = NULL;
- double u = 0; /* Linear interpolation parameter */
double x = 0; /* Argument corresponding to input value y */
size_t nitems = 0;
ASSERT(tab && y >= 0 && y < 1);
@@ -72,17 +84,24 @@ swf_tabulation_inverse(const struct swf_tabulation* tab, const double y)
nitems = darray_item_size_get(&tab->items);
ASSERT(nitems);
- find = search_lower_bound(&y, items, nitems, sizeof(*items), cmp_item);
- ASSERT(find && find->f_x >= y);
+ find = search_lower_bound(&y, items, nitems, sizeof(*items), cmp_fx);
+
+ if(y == 0) return 0;
+
+ /* Input y is not in the tabulated range: returns the nearest x */
+ if(!find) return items[nitems-1].x;
+ if(find == items) return find->x;
+
+ ASSERT(find->fx >= y);
/* The input y correspond exactly to a tabulated argument */
- if(find->f_x == y) return find->x;
+ if(find->fx == y) return find->x;
/* Linear interpolation of tabulated arguments containing the argument
* corresponding to the input value */
ASSERT(find > items && find[-1].f_x < y);
- u = (y - find[-1].f_x) / (find[0].f_x - find[-1].f_x);
- x = u*(find[0].x - find[-1].x) + find[-1].x;
+ x = linear_interpolation(y, &find[-1], &find[0]);
+
return x;
}
diff --git a/src/swf_tabulation.h b/src/swf_tabulation.h
@@ -21,9 +21,12 @@
struct item {
double x; /* Function argument */
- double f_x; /* Value of the function */
+ double fx; /* Value of the function */
+ double dfx; /* Value of the derivative of the function */
+ double d2fx; /* Value of the second derivative of the function */
+ double d3fx; /* Value of the third derivative of the function */
};
-#define ITEM_NULL__ {0, 0}
+#define ITEM_NULL__ {0, 0, 0, 0, 0}
static const struct item ITEM_NULL = ITEM_NULL__;
/* Declare the dynamic array of items */
