[SciPy-Dev] Mass 1D interpolation, advice

[josef.pktd at gmail.com]

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

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