Hi David, I'm moving this thread over to gsl-discuss from help-gsl.
I've been taking a look at your interp2d code since I'm currently doing
some work which needs 2d interpolation. I like your coding style and
think this should be moved into the repository at some point.
Just a few initial thoughts:
1. Your INDEX_2D macro for indexing the z array appears to use
column-major ordering, even though the comment in the .h file says
row-major. IE: it is currently defined as:
#define INDEX_2D(xi, yi, xsize, ysize) ((yi) * (xsize) + (xi))
which would store each column contiguously in memory instead of each
row. The GSL gsl_matrix structure uses row-major ordering (see
gsl_matrix_double.h):
283 return m->data[i * m->tda + j] ;
and so if a user wants to store their 2D grid in a gsl_matrix, its not
straightforward to just pass m->data as the z array argument to your
functions; they'd have to transpose first.
2. It would be really great if you could document your library with
texinfo so when we fold the code into the repository the docs are ready
to go. You can look at some other extensions for examples (ie see the
doc/alf.texi file in the ALF extension)
Patrick
On 08/16/2013 05:09 PM, David Zaslavsky wrote:
> Hi there,
>
> For quite some time I've been working on a 2D interpolation library
> compatible with the GSL. I've mentioned it a couple times on this list.
> With some recent work it's gotten to the point where I consider it
> ready for a first release. How would I go about getting it listed in the
> "Extensions" section on the GSL web page?
>
> Also if anyone would like to try compiling and running the test suite,
> I'd appreciate knowing about any bugs that pop up. I've only tried it on
> Linux with GCC 4.6.3.
> https://github.com/diazona/interp2d
>
> :) David
>
[-- Attachment #1: Type: text/plain, Size: 4334 bytes --] I just noticed this message seems not to have made it out to the list when I originally sent it on Friday. So I'm trying again. ----- Hi Patrick, Thanks for the feedback. First let me say that texinfo documentation will not be a problem. I'm happy to add it, but it's just a matter of time - I have a lot of other things going on so it will probably take on the order of some weeks for me to find time to do that. On the indexing issue: interp2d treats the x index as the column coordinate and the y index as the row coordinate, so that xsize is the number of columns and ysize is the number of rows. This is done for compatibility with Cartesian coordinates: most people envision x as the horizontal coordinate and y as the vertical one. Given this interpretation, the current implementation of INDEX_2D does use row-major order. If I'm not mistaken, it's still possible to use a GSL matrix's data array in interp2d without transposing if you pass the y index as i and the x index as j to gsl_matrix_get and friends. gsl_matrix_alloc(ysize, xsize) gsl_matrix_get(m, yi, xi) etc. So I would say the real issue is not the implementation, but the variable naming. I think it's worth putting in some thought on how best to resolve this without confusing people. For compatibility with gsl_matrix, it might be better to ditch the names xi, yi, xsize, ysize in favor of j, i, size2, size1 respectively. But I must admit I find the x and y naming convention much more intuitive. I guess something else I could do is switch the order of the x and y arguments in the interp2d functions, to put yi, ysize, etc. first. That would make it better correspond with the standard notation in which the first index labels the row of the matrix and the second index labels the column. But then again it would cause some headaches with people who are already using the code. (I'm not sure how many users the library has.) Also, I should mention that I don't consider interp2d to be feature-complete, nor am I in any particular rush to make it complete. There are some features I'm interested in adding, including additional interpolation schemes, more interpolation methods without boundary checks (for extrapolation), etc. How would this move affect my ability to continue to add features? :) David On 01/03/2014 02:00 PM, Patrick Alken wrote: > Hi David, I'm moving this thread over to gsl-discuss from help-gsl. > > I've been taking a look at your interp2d code since I'm currently doing > some work which needs 2d interpolation. I like your coding style and > think this should be moved into the repository at some point. > > Just a few initial thoughts: > > 1. Your INDEX_2D macro for indexing the z array appears to use > column-major ordering, even though the comment in the .h file says > row-major. IE: it is currently defined as: > > #define INDEX_2D(xi, yi, xsize, ysize) ((yi) * (xsize) + (xi)) > > which would store each column contiguously in memory instead of each > row. The GSL gsl_matrix structure uses row-major ordering (see > gsl_matrix_double.h): > > 283 return m->data[i * m->tda + j] ; > > and so if a user wants to store their 2D grid in a gsl_matrix, its not > straightforward to just pass m->data as the z array argument to your > functions; they'd have to transpose first. > > 2. It would be really great if you could document your library with > texinfo so when we fold the code into the repository the docs are ready > to go. You can look at some other extensions for examples (ie see the > doc/alf.texi file in the ALF extension) > > Patrick > > On 08/16/2013 05:09 PM, David Zaslavsky wrote: >> Hi there, >> >> For quite some time I've been working on a 2D interpolation library >> compatible with the GSL. I've mentioned it a couple times on this list. >> With some recent work it's gotten to the point where I consider it >> ready for a first release. How would I go about getting it listed in the >> "Extensions" section on the GSL web page? >> >> Also if anyone would like to try compiling and running the test suite, >> I'd appreciate knowing about any bugs that pop up. I've only tried it on >> Linux with GCC 4.6.3. >> https://github.com/diazona/interp2d >> >> :) David >> > [-- Attachment #2: OpenPGP digital signature --] [-- Type: application/pgp-signature, Size: 295 bytes --]
On 01/05/2014 01:31 PM, David Zaslavsky wrote: > I just noticed this message seems not to have made it out to the list > when I originally sent it on Friday. So I'm trying again. > > ----- > > Hi Patrick, > > Thanks for the feedback. First let me say that texinfo documentation > will not be a problem. I'm happy to add it, but it's just a matter of > time - I have a lot of other things going on so it will probably take on > the order of some weeks for me to find time to do that. > > > On the indexing issue: interp2d treats the x index as the column > coordinate and the y index as the row coordinate, so that xsize is the > number of columns and ysize is the number of rows. This is done for > compatibility with Cartesian coordinates: most people envision x as the > horizontal coordinate and y as the vertical one. Given this > interpretation, the current implementation of INDEX_2D does use > row-major order. This makes sense, we can probably leave it the way it is, we'll just need to make sure the documentation is very clear on how to construct the z array so users aren't confused with the indexing. > > > Also, I should mention that I don't consider interp2d to be > feature-complete, nor am I in any particular rush to make it complete. > There are some features I'm interested in adding, including additional > interpolation schemes, more interpolation methods without boundary > checks (for extrapolation), etc. How would this move affect my ability > to continue to add features? It already has bilinear and bicubic implemented. This is more than enough to make it useful to people so I think we can put it in now and continue adding features. I saw from some comments that you used some code from ALGLIB - can we assume you're using the GPL v2 version of ALGLIB? In that case there shouldn't be any problem putting it in GSL under GPL v3. The next steps would be to get the code ready for inclusion in GSL. This includes: 1. renaming external functions to have the gsl_ prefix and make sure all other functions are declared static 2. documentation 3. tests 4. I suggest checking out the 'master' branch of gsl and making your own branch 'interp2d' or similar. Once your branch is ready we can discuss how to merge it - you could either post it to somewhere like github or if you'd like, I don't see any problem with giving you access to savannah git so you can continue to add features, etc. Patrick
[-- Attachment #1: Type: text/plain, Size: 2314 bytes --] Hi Patrick, On 01/06/2014 11:57 AM, Patrick Alken wrote: > On 01/05/2014 01:31 PM, David Zaslavsky wrote: >> On the indexing issue: interp2d treats the x index as the column >> coordinate and the y index as the row coordinate, so that xsize is the >> number of columns and ysize is the number of rows. This is done for >> compatibility with Cartesian coordinates: most people envision x as the >> horizontal coordinate and y as the vertical one. Given this >> interpretation, the current implementation of INDEX_2D does use >> row-major order. > > This makes sense, we can probably leave it the way it is, we'll just > need to make sure the documentation is very clear on how to construct > the z array so users aren't confused with the indexing. Sure, that would be a reasonable option, but I'd like to do a bit of research (i.e. ask a few more people for opinions) to see what solution will be the least confusing. After all, if there is a change to be made, this is the time to do it, not later after the code is incorporated into GSL. I think backward compatibility will not be too much of a problem because people will have to make other changes anyway to transition from the interp2d_* functions to the corresponding gsl_* functions. >> Also, I should mention that I don't consider interp2d to be >> feature-complete, nor am I in any particular rush to make it complete. >> There are some features I'm interested in adding, including additional >> interpolation schemes, more interpolation methods without boundary >> checks (for extrapolation), etc. How would this move affect my ability >> to continue to add features? > > It already has bilinear and bicubic implemented. This is more than > enough to make it useful to people so I think we can put it in now and > continue adding features. > > I saw from some comments that you used some code from ALGLIB - can we > assume you're using the GPL v2 version of ALGLIB? In that case there > shouldn't be any problem putting it in GSL under GPL v3. The code was adapted from ALGLIB by Thomas Beutlich, and I'm pretty sure he told me that he used the GPL v2 version. I'll double-check that, though. When I have some time I'll fork the GSL repository and work on adding the interp2d code in a new branch. :) David [-- Attachment #2: OpenPGP digital signature --] [-- Type: application/pgp-signature, Size: 295 bytes --]
[-- Attachment #1: Type: text/plain, Size: 2391 bytes --] Hi Patrick, Sorry for the delayed response; I've been quite busy this week. On 01/06/2014 11:57 AM, Patrick Alken wrote: > On 01/05/2014 01:31 PM, David Zaslavsky wrote: >> On the indexing issue: interp2d treats the x index as the column >> coordinate and the y index as the row coordinate, so that xsize is the >> number of columns and ysize is the number of rows. This is done for >> compatibility with Cartesian coordinates: most people envision x as the >> horizontal coordinate and y as the vertical one. Given this >> interpretation, the current implementation of INDEX_2D does use >> row-major order. > > This makes sense, we can probably leave it the way it is, we'll just > need to make sure the documentation is very clear on how to construct > the z array so users aren't confused with the indexing. Sure, that would be a reasonable option, but I'd like to do a bit of research (i.e. ask a few more people for opinions) to see what solution will be the least confusing. After all, if there is a change to be made, this is the time to do it, not later after the code is incorporated into GSL. I think backward compatibility will not be too much of a problem because people will have to make other changes anyway to transition from the interp2d_* functions to the corresponding gsl_* functions. >> Also, I should mention that I don't consider interp2d to be >> feature-complete, nor am I in any particular rush to make it complete. >> There are some features I'm interested in adding, including additional >> interpolation schemes, more interpolation methods without boundary >> checks (for extrapolation), etc. How would this move affect my ability >> to continue to add features? > > It already has bilinear and bicubic implemented. This is more than > enough to make it useful to people so I think we can put it in now and > continue adding features. > > I saw from some comments that you used some code from ALGLIB - can we > assume you're using the GPL v2 version of ALGLIB? In that case there > shouldn't be any problem putting it in GSL under GPL v3. The code was adapted from ALGLIB by Thomas Beutlich, and I'm pretty sure he told me that he used the GPL v2 version. I'll double-check that, though. When I have some time I'll fork the GSL repository and work on adding the interp2d code in a new branch. :) David [-- Attachment #2: OpenPGP digital signature --] [-- Type: application/pgp-signature, Size: 295 bytes --]
> Sure, that would be a reasonable option, but I'd like to do a bit of
> research (i.e. ask a few more people for opinions) to see what solution
> will be the least confusing. After all, if there is a change to be made,
> this is the time to do it, not later after the code is incorporated into
> GSL. I think backward compatibility will not be too much of a problem
> because people will have to make other changes anyway to transition from
> the interp2d_* functions to the corresponding gsl_* functions.
>
I think 2 additional functions could be added to help solve this issue
and hide the indexing
from the user:
interp2d_grid_set(double zarr[], size_t i, size_t j, double z)
{
/* set point (i,j) of grid 'zarr' to value 'z' */
}
and
interp2d_grid_fill(double zarr[], int (*grid_func)(double x, double y,
void *params))
{
/* fill entire grid 'zarr' using callback function 'grid_func'; this
would be useful
* in cases where the user has a continuous/analytic function which
is expensive to call, but
* they'd like to make a grid and then interpolate from it
*/
}
Hi all,
I consider the addition of these routines to GSL as a major breakthrough.
Now I have three questions about generalization:
1. can these 2D interpolation routines be extended to scattered (i.e. non-
rectangular) data?
2. can these 2D interpolation routines be extended to n dimensions?
3. both
Thank you for your hard work!
--
jeremy
On Friday 10 January 2014 09:37:04 Patrick Alken wrote:
> > Sure, that would be a reasonable option, but I'd like to do a bit of
> > research (i.e. ask a few more people for opinions) to see what solution
> > will be the least confusing. After all, if there is a change to be made,
> > this is the time to do it, not later after the code is incorporated into
> > GSL. I think backward compatibility will not be too much of a problem
> > because people will have to make other changes anyway to transition from
> > the interp2d_* functions to the corresponding gsl_* functions.
>
> I think 2 additional functions could be added to help solve this issue
> and hide the indexing
> from the user:
>
> interp2d_grid_set(double zarr[], size_t i, size_t j, double z)
> {
> /* set point (i,j) of grid 'zarr' to value 'z' */
> }
>
> and
>
> interp2d_grid_fill(double zarr[], int (*grid_func)(double x, double y,
> void *params))
> {
> /* fill entire grid 'zarr' using callback function 'grid_func'; this
> would be useful
> * in cases where the user has a continuous/analytic function which
> is expensive to call, but
> * they'd like to make a grid and then interpolate from it
> */
> }
On 01/10/2014 05:27 PM, jeremy theler wrote: > Hi all, > > I consider the addition of these routines to GSL as a major breakthrough. > Now I have three questions about generalization: > > 1. can these 2D interpolation routines be extended to scattered (i.e. non- > rectangular) data? Not currently, but this can always be done with least-squares models. I'm not aware of any available software that interpolates on nonuniform grids - do you know of any? > > 2. can these 2D interpolation routines be extended to n dimensions? Once the 2D code is imported I'd like to make at least a 3D linear interpolater, which should be straightforward.
[-- Attachment #1: Type: text/plain, Size: 2469 bytes --] Hi Jeremy, The algorithms used by interp2d require a rectangular grid. Doing interpolation without that grid would require an entirely different algorithm, so that would be a separate project. The algorithms can be generalized to higher-dimensional interpolation, but it becomes exponentially more complicated as the number of dimensions increases. For linear interpolation this probably wouldn't be too bad, but other interpolation schemes would be quite a pain to code in a higher number of dimensions. I believe it would be possible to make an interpolation routine that automatically generalizes to any number of dimensions (I think Mathematica does this, for example), but that would require pretty much rewriting the entire code base. :) David On 01/10/2014 07:27 PM, jeremy theler wrote: > Hi all, > > I consider the addition of these routines to GSL as a major breakthrough. > Now I have three questions about generalization: > > 1. can these 2D interpolation routines be extended to scattered (i.e. non- > rectangular) data? > > 2. can these 2D interpolation routines be extended to n dimensions? > > 3. both > > > Thank you for your hard work! > > -- > jeremy > > > On Friday 10 January 2014 09:37:04 Patrick Alken wrote: >>> Sure, that would be a reasonable option, but I'd like to do a bit of >>> research (i.e. ask a few more people for opinions) to see what solution >>> will be the least confusing. After all, if there is a change to be made, >>> this is the time to do it, not later after the code is incorporated into >>> GSL. I think backward compatibility will not be too much of a problem >>> because people will have to make other changes anyway to transition from >>> the interp2d_* functions to the corresponding gsl_* functions. >> I think 2 additional functions could be added to help solve this issue >> and hide the indexing >> from the user: >> >> interp2d_grid_set(double zarr[], size_t i, size_t j, double z) >> { >> /* set point (i,j) of grid 'zarr' to value 'z' */ >> } >> >> and >> >> interp2d_grid_fill(double zarr[], int (*grid_func)(double x, double y, >> void *params)) >> { >> /* fill entire grid 'zarr' using callback function 'grid_func'; this >> would be useful >> * in cases where the user has a continuous/analytic function which >> is expensive to call, but >> * they'd like to make a grid and then interpolate from it >> */ >> } [-- Attachment #2: OpenPGP digital signature --] [-- Type: application/pgp-signature, Size: 295 bytes --]
On Friday 10 January 2014 18:55:37 Patrick Alken wrote: > On 01/10/2014 05:27 PM, jeremy theler wrote: > > Hi all, > > > > I consider the addition of these routines to GSL as a major breakthrough. > > > > Now I have three questions about generalization: > > 1. can these 2D interpolation routines be extended to scattered (i.e. > > non- > > > > rectangular) data? > > Not currently, but this can always be done with least-squares models. > I'm not aware of any available software that interpolates on nonuniform > grids - do you know of any? I think I saw a python numerical library (whose name I forgot) that performed some kind of scattered-data interpolation. I think it used a weigthed average of first n neighbors. What I do in these cases is a nearest-neighbors search using a k-dimensional tree to sort the data. Some simple 2D examples can be seen in http://www.talador.com.ar/jeremy/wasora/realbook/._realbook010.html For two dimensional functions such as f(x,y) I once tried to do a 3-nearest neighbors search, fit a plane as z = a*x + b*y + c an then evaluate the plane z at the desired (x,y) location, but that approach had more problems than advantages. Actually, wasora is a piece of software I am writing that acts as a high-level interface to GSL for most of its features (except multidim interpolation for now, and ODE/DAE systems for which it uses SUNDIALS). Some examples of the GSL manual are rewritten as wasora inputs in http://www.talador.com.ar/jeremy/wasora/realbook/._realbook015.html If someone wants to try, comments are welcome. > > 2. can these 2D interpolation routines be extended to n dimensions? > > Once the 2D code is imported I'd like to make at least a 3D linear > interpolater, which should be straightforward. Good. Are n dimensions straightforward also? -- jeremy
> Good. Are n dimensions straightforward also?
I'll put it on the todo list, but feel free to make one yourself and
contribute it if you'd like