[SciPy-Dev] Mass 1D interpolation, advice

[josef.pktd at gmail.com]

On Wed, Apr 15, 2015 at 11:24 AM,  <josef.pktd at gmail.com> wrote:
> On Wed, Apr 15, 2015 at 11:04 AM,  <josef.pktd at gmail.com> wrote:
>> On Wed, Apr 15, 2015 at 10:33 AM, Jonathan Stickel <jjstickel at gmail.com> wrote:
>>> On 4/14/15 16:12 , Matthew Brett wrote:
>>>> Hi,
>>>>
>>>> In brain imaging (my world) - we need to do 1D interpolation of large
>>>> datasets, specifically, we have data in the order of (X, Y) = (200,
>>>> 4000), where we are interpolating at the same ~200 X locations for
>>>> each of the 4000 rows.
>>>>
>>>> We need to do this with extrapolation at the boundaries, so I think we
>>>> can't use scipy.interpolate.interp1d which is otherwise well suited
>>>> for the task.
>>>>
>>>> scipy.interpolate.InterpolatedUnivariateSpline does deal with
>>>> extrapolation, but only works with vectors.
>>>>
>>>> I think, in order to make speed reasonable, I will need to loop over
>>>> the rows in compiled code.
>>>>
>>>> My questions are:
>>>>
>>>> a) how easy would it be to try and link some external Cython code
>>>> against the fitpack code in scipy?
>>>> b) if that is not easy, what can I do that will also be of most use to
>>>> other scipy users?
>>>>
>>>> Cheers,
>>>>
>>>> Matthew
>>>
>>> I am not sure if I understand your problem statement. Most interpolation
>>> problems are posed as:  given data (x,y), determine yp at xp, where xp
>>> are different from x. This is fine for cases where you have a single y
>>> and possibly multiple sets xp. I had case awhile ago where I had a
>>> single xp and multiple y (at the same x). In this case, the problem
>>> statement is:  given known locations x and interpolant locations xp,
>>> what are the values for yp (corresponding to xp) given y (corresponding
>>> to xp). To solve this problem efficiently, I wrote the code I posted here:
>>>
>>> http://scipy-central.org/item/21/1/simple-piecewise-polynomial-interpolation
>>
>> nice,
>>
>> 3 questions
>>
>> Can you construct P = Xp*inv(X)*H  directly from the block diagonal
>> parts instead of going through sparse block-diagonal matrices?
>>
>> Why do you need monotonically increasing x? I would have thought it
>> works for any x
>>
>> (I'm not clear about this)
>>
>> My guess is that replacing inv by pinv would be more stable. But,
>> would pinv also be able to handle the case when there are some
>> identical x (by implicitly taking an average)?
>
> Something is wrong in the way I think about this.
> (I might be thinking about a different use case where I was recently
> using Kernel Ridge Regression.)
>
> Sorry for any noise


The projection matrix for interpolation would need to have shape
(number of new points to be interpolated, number of sample points)

Josef

>
> Josef
>
>>
>> Josef
>>
>>>
>>> Although it is not true spline interpolation (rather a simpler local
>>> polynomial interpolation), it is reasonably fast and accurate.
>>> Extrapolation is allowed. However, I later found that
>>> scipy.interpolate.splrep/splev are so fast that they could be used
>>> instead with equivalent speed for my problem. By comparison,
>>> scipy.interpolate.interp1d is quite slow.
>>>
>>> HTH,
>>> Jonathan
>>> _______________________________________________
>>> SciPy-Dev mailing list
>>> SciPy-Dev at scipy.org
>>> http://mail.scipy.org/mailman/listinfo/scipy-dev