README.md 5.48 KB
Newer Older
Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
1 2
# librinterpolate

3
librinterpolate: a library to perform linear interpolation on a constant gridded dataset.
Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
4
(c) Robert Izzard 2005-2019
5

6
Usage:
7 8 9

The table should be organised (in memory) like this

Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
10
```
11 12 13
p1[0] p2[0] p3[0] ... d1[0] d2[0] d3[0] d4[0] ...
p1[1] p2[1] p3[1] ... d1[1] d2[1] d3[1] d4[1] ...
p1[2] p2[2] p3[2] ... d1[2] d2[2] d3[2] d4[2] ...
Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
14 15
```

16 17 18 19

so there are n parameters (p1, p2, p3...) and d data items
per data line (d1, d2, d3, d4 ...). There are l lines of data.

20
The parameters should be ordered from low to high values.
Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
21
The parameters should be on a constant* grid (which does NOT
22 23 24 25 26 27
need to be regular).

What does this mean?

This is a good data table:

Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
28
```
29 30 31 32 33 34 35 36 37 38 39
0.1 -100 10 ...data...
0.1 -100 25 ...data...
0.1 -100 30 ...data...
0.1  -50 10 ...data...
0.1  -50 25 ...data...
0.1  -50 30 ...data...
0.3 -100 10 ...data...
0.3 -100 25 ...data...
0.3 -100 30 ...data...
0.3  -50 10 ...data...
0.3  -50 25 ...data...
40
0.3  -50 30 ...data...
41 42 43 44 45
0.9 -100 10 ...data...
0.9 -100 25 ...data...
0.9 -100 30 ...data...
0.9  -50 10 ...data...
0.9  -50 25 ...data...
46
0.9  -50 30 ...data...
Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
47 48
```

49 50 51 52 53 54 55 56 57

In the above case, the parameters have values:
p0 : 0.1,0.3,0.9
p1 : -100, -50
p2 : 10,25,30

The parameter "hypercube" then is the 3D-hypercube of
323 = 18 data lines.

Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
58
Note that the points on the cube are constant but not regularly spaced,
59 60 61 62 63 64
e.g. p0 has spacing (0.3-0.1)=0.2 and (0.9-0.3)=0.6, which are different,
BUT e.g. the spacing for p1 (-100 - -50 = -50) is the same, whatever the
value of p0. The same is true for p3 which always has spacings 15 and 5
from (25-10) and (30-25) respectively.

Note that the table is assumed to be sorted from SMALLEST
Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
65 66
to LARGEST parameter values. It is also assumed to be regular and fully filled.
So no missing data please.
67

Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
68
To interpolate data, n parameters are passed into this
69
routine in the array x. The result of the interpolation is put
70
into the array r (of size d).
71

Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
72 73
If you enable RINTERPOLATE_CACHE (set on by default) then results are cached to
avoid slowing the code too much.  This means
74 75 76 77 78 79
that the interpolate routine checks your input parameters x against
the last few sets of parameters passed in. If you have used these x
recently, the result is simply returned. This saves extra interpolation
and is often faster.  This is only true in some cases of course - if your
x are always different you waste time checking the cache. This is why
the cache_hint variable exists: if this is false then the cache is skipped.
Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
80
Of course only you know if you are likely to call the interpolate routine
81 82
repeated with the same values... I cannot possibly know this in advance!

Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
83
The interpolation process involves finding the lines of the data table
84 85 86 87 88
which span each parameter x. This makes a hypercube of length 2^n (e.g.
in the above it is 8, for simple 1D linear interpolation it would be the
two spanning values). Linear interpolation is then done in the largest
dimension, above this is the third parameter (p2), then the second (p1)
and finally on the remaining variable (p0). This would reduce the table
89
from 2^3 lines to 2^2 to 2^1 to (finally) 2^0 i.e. one line which is the
90 91
interpolation result.

Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
92 93
To find the spanning lines a binary search is performed. This code was 
originally donated by Evert Glebbeek. See e.g.
94
http://en.wikipedia.org/wiki/Binary_search_algorithm
95
and note the comment "Although the basic idea of binary search is comparatively
Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
96
straightforward, the details can be surprisingly tricky... " haha!
97 98 99 100 101 102

Each table has its own unique table_id number. This is just to allow
us to set up caches (one per table) and save arrays such as varcount and
steps (see below) so they only have to be calculated once.
Since binary_c2 these table_id numbers have been allocated automatically,
you do not have to set one yourself, but this does assume that the tables
Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
103
are each at fixed memory locations - so please do not move the data around once you start to use it.
104

Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
105
Example code is given in test_rinterpolate.c and reproduced here
106

Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
107 108 109
```
    /* Number of parameters */
    const rinterpolate_counter_t N = 2;
110

Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
111 112
     /* Number of data */
    const rinterpolate_counter_t D = 3;
113

Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
114 115
    /* length of each line (in doubles) i.e. N+D */
    const rinterpolate_counter_t ND = N + D;
116

Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
117 118 119
    /* total number of lines */
    const rinterpolate_counter_t L = 100;

120
    /* make rinterpolate data (for cache etc.) */
Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
121
    struct rinterpolate_data_t * rinterpolate_data = NULL;
122
    rinterpolate_counter_t status = rinterpolate_alloc_dataspace(&rinterpolate_data);
Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
123 124

    /* data table : it is up to you to set the data in here*/
125 126
    rinterpolate_float_t table[ND*L];

Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
127
    /* ... set the data ... */
128 129


Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
130
    /*
131 132
     * Arrays for the interpolation location and
     * interpolated data. You need to set
Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
133 134 135
     * the interpolation location, x.
     */
    rinterpolate_float_t x[N],r[D];
136

Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
137
    /* ... set the x array ... */
138

Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
139
    /* choose whether to cache (0=no, 1=yes) */
140 141
    int usecache = 1;

Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
142 143 144 145 146 147 148 149 150 151 152
    /* do the interpolation */
    rinterpolate(table,
                 rinterpolate_data,
                 N,
                 D,
                 L,
                 x,
                 r,
                 usecache);

    /* the array r contains the result */
153 154


Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
155
    /* ... rest of code ... */
156 157


Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
158 159
    /* free memory on exit */
    rinterpolate_free_data(rinterpolate_data);
Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
160
    free(rinterpolate_data);
161

Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
162 163
```

164 165


166 167
I have optimized this as best I can, please let me know if
you can squeeze any more speed out of the function.
168
I am sorry this has made most of the function unreadable! The
169 170
comments should help, but you will need to know some tricks...

Izzard, Robert Dr (Physics)'s avatar
Izzard, Robert Dr (Physics) committed
171
(c) Robert Izzard, 2005-2020, please send bug fixes!