interp.c

   1 /* A full-fledged interpolating lookup table implementation. */
   2
   3 #include <stdlib.h> /* malloc, free */
   4 #include <stdio.h>  /* printf       */
   5 #include <math.h>   /* sin          */
   6 #include <assert.h> /* assert       */
   7 #include "tavl.h"   /* tavl_*       */
   8 #include "interp.h"
   9
  10 typedef struct evaluated_func_struct {
  11   double x;
  12   double f_of_x;
  13 } evaluated_func_t;
  14
  15 /* Comparison function for pointers to |evaluated_func_t|s.
  16    |param| is not used. */
  17 /* adjusted from compare_ints in test.c */
  18 int
  19 compare_evaluated_func_ts (const void *pa, const void *pb, void *param)
  20 {
  21   const evaluated_func_t *a = pa;
  22   const evaluated_func_t *b = pb;
  23
  24   if (a->x < b->x)
  25     return -1;
  26   else if (a->x > b->x)
  27     return +1;
  28   else
  29     return 0;
  30 }
  31
  32 void *allocate_evaluated_func_t (double x, double f_of_x)
  33 {
  34   evaluated_func_t *eft=NULL;
  35   eft = (evaluated_func_t *)malloc(sizeof(evaluated_func_t));
  36   assert(eft != NULL);
  37   eft->x = x;
  38   eft->f_of_x = f_of_x;
  39   return eft;
  40 }
  41
  42 void free_evaluated_func_t (void *tavl_item, void *tavl_param)
  43 {
  44   free(tavl_item);
  45 }
  46
  47
  48 interp_table_t *interp_table_allocate (interp_fn_t *fn, interp_tol_fn_t *tol)
  49 {
  50   interp_table_t *itable=NULL;
  51   itable = (interp_table_t *)malloc(sizeof(interp_table_t));
  52   assert(itable != NULL);
  53   itable->table = (void *)tavl_create(&compare_evaluated_func_ts, NULL, NULL);
  54   itable->fn = fn;
  55   itable->tol = tol;
  56   return itable;
  57 }
  58
  59 void interp_table_free (interp_table_t *itable)
  60 {
  61   tavl_destroy((struct tavl_table *)itable->table, &free_evaluated_func_t);
  62   free(itable);
  63 }
  64
  65 double interp_table_eval (interp_table_t *itable, void *params, double x)
  66 {
  67   struct tavl_table *table = (struct tavl_table *)itable->table;
  68   const evaluated_func_t *a, *b, *c;
  69   double y, weight;
  70
  71   /* attempt to look up the x variable */
  72   a = tavl_bracket(table, &x, (const void **)&b, (const void **)&c);
  73   if (a != NULL) /* direct hit */
  74     y = a->f_of_x;
  75   else {
  76     if (c == NULL || b == NULL) { /* out of interpolation range, run an actual eval */
  77       y = (*itable->fn)(x, params);
  78       tavl_insert(table, allocate_evaluated_func_t (x,y) );
  79     } else {
  80       //printf("bracketed x=%g with %g, %g\n", x, b->x, c->x);
  81       /* should we interpolate? */
  82       if ( (*itable->tol)(b->x, b->f_of_x, c->x, c->f_of_x) == 0 ) {
  83         /* interpolate */
  84         weight = (x - b->x) / (c->x - b->x);
  85         y = b->f_of_x + (c->f_of_x - b->f_of_x)*weight;
  86         //printf("interpolated f(%g) ~= %g\n", x, y);
  87         //fprintf(stderr, "%g\t%g\n", x, y);
  88       } else { /* tolerance exceeded, run an actual eval */
  89         y = (*itable->fn)(x, params);
  90         tavl_insert(table, allocate_evaluated_func_t (x,y) );
  91       }
  92     }
  93   }
  94   return y;
  95 }