[SciPy-Dev] Multilinear interpolation

[Pablo Winant]

On 25/02/2013 14:32, Pauli Virtanen wrote:
> Pablo Winant <pablo.winant <at> gmail.com> writes:
> [clip]
>> So I made an attempt using tempita in the case of regularly spaced
>> grids. For now, it is in a temporary repository in github:
>>
>> https://github.com/albop/python_interpolation
> [clip]
>
> Some comments:
>
> - "#pragma omp parallel for" statement probably doesn't do
>    anything: Cython interprets it as a comment.

yes, my mistake, I intended to write "with nogil, parallel():" instead
> - The code assumes C-contiguous arrays, but does not check for it.
right
> - Do you actually need code generation here? The templated loops
>    can be replaced by N-dim untemplated loops. The question is
>    probably mostly about speed.
I wonder too. The motivation was indeed speed but there may be a better 
way.

> - The inner evaluation loop should be inside `with nogil`
>
> [clip]
>> Before, I do the same for irregularly spaced grid, I have some questions:
>> - I included the templates as Python strings in the .pyx.in file . Was
>> there a better way ?
> As the script is executable by Python, the file extension should be .py.
> You can either put the template itself to a separate file (.pyx.in)
> --- note that you can put Python code inside tempita templates.
> Alternatively, rename the script to `generate_multilinear.py`.
>
> Note that tempita templates can contain executable blocks, so you can
> for instance put some logic there, and use the template syntax for the
> loops.
OK, as the python code is minimal, I will try to include everything in 
the template.

>> - I was not sure about how to deal with single/double precision in the
>> cython code.
> You can use fused types instead of templating; the end result is the same,
> but Cython takes care of picking the correct routine.
>
I tried it but I didn't see how to do it. I would like to do define a 
function like

def function(ndarray[fused_type, ndim=2]):
     ...
     cdef fused_type internal_variable
     ...
     return something


but it doesn't work. From this thread 
[https://groups.google.com/forum/?fromgroups=#!topic/cython-users/GwvZiAdRJpM] 
I understand I should do :

ctypedef fused  myarray:
     ndarray[type1, ndim=2]
     ndarray[type2, ndim=2]

def function(myarray):
     ...
     cdef (?) internal_variable

but then how do I define scalars of the same dtype as the array arguments ?