[SciPy-user] Fwd: scipy.interpolate.Rbf() MemoryError for more than 4000 (x, y) pairs ?

josef.pktd@gmai... josef.pktd@gmai...
Fri Apr 3 11:50:07 CDT 2009


send again, previous message bounced


---------- Forwarded message ----------
From:  <josef.pktd@gmail.com>
Date: Fri, Apr 3, 2009 at 12:36 PM
Subject: Re: [SciPy-user] scipy.interpolate.Rbf() MemoryError for more
than 4000 (x, y) pairs ?
To: SciPy Users List <scipy-user@scipy.org>


On Fri, Apr 3, 2009 at 12:22 PM,  <josef.pktd@gmail.com> wrote:
> On Fri, Apr 3, 2009 at 11:40 AM, Jim Vickroy <Jim.Vickroy@noaa.gov> wrote:
>> Hello Group,
>>
>> I'm a radial basis functions novice so I may be making a mistake; however,
>> it seems that scipy.interpolate.Rbf() can handle a maximum of ~4000 (x,y)
>> points.  Is that correct?
>>
>> Here is a sample output from the attached script (rbf-demo.py) for 5000
>> (x,y) points:
>>
>> <output>
>> DOS>rbf-demo.py
>> Python version: 2.5.4 (r254:67916, Dec 23 2008, 15:10:54) [MSC v.1310 32 bit
>> (Intel)]
>> numpy version:  1.2.1
>> scipy version:  0.7.0
>> multiquadric         in-painting required 200 seconds for 5000 points
>> Traceback (most recent call last):
>> File "rbf-demo.py", line 28, in <module>
>>  rbf   = scipy.interpolate.Rbf(x, y, z, function=function)
>> File "C:\Python25\lib\site-packages\scipy\interpolate\rbf.py", line 132, in
>> __init__
>>  r = self._call_norm(self.xi, self.xi)
>> File "C:\Python25\lib\site-packages\scipy\interpolate\rbf.py", line 147, in
>> _call_norm
>>  return self.norm(x1, x2)
>> File "C:\Python25\lib\site-packages\scipy\interpolate\rbf.py", line 100, in
>> _euclidean_norm
>>  return sqrt( ((x1 - x2)**2).sum(axis=0) )
>> MemoryError
>>
>> DOS>
>> </output>
>>
>> Thanks,
>> -- jv
>>
>> P.S.
>> I'm considering using Rbf to replace a horizontal band of bad-pixel values
>> in 512x512 images obtained from a damaged detector.
>>
>
>
> Besides the memory error, I don't think rbf can do this, rbf uses all
> points to do the interpolation, the kernel is number of observation
> squared. In some examples, that I tried, if the number of points is to
> large, then you don't get a smooth solution, the interpolation has too
> much variation (like fitting a very high order polynomial).
>
> From my experimentation with rbf, it only works well with a small
> number of points, so the best is to run rbf on a local neighborhood of
> the missing points.
>
> For 1d, I have the script somewhere but I never tried the 2d version.
>
> Josef
>

here is the script, it produces 12 graphs to illustrate this

Josef
-------------- next part --------------
'''script to try out scipy.interpolate.rbf

rbf solves a system of equations with the number of equations and variables
equal to the number of node points. This produces reasonable results for up
to 10 or 20 points, but for 100 points the results are very jagged.

Scott Sinclair recommended using moving windows, which I implemented for the
one dimensional case. For the higher dimensional case, an efficient neighbor
search could be done with the scipy.distance functions.

RBF kernel regression ?

Josef
'''

import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate


def show_rbf_example(basefn,n):
    '''example for global rbf

    Parameters
    ----------

    basefn: function
        sample points are created with basefn plus noise
    n: integer
        number of sample points to create

    returns
    -------
    
    None
    output: one plots

    Notes
    -----

    If the number of sample points is larger than 20 (in the examples),
    then the interpolation is very jagged.
    '''
    xr = np.linspace(0,3,n)
    yr = basefn(xr)+0.1*np.random.randn(xr.size)
    plt.figure()
    plt.plot(xr,yr,'o')

    xi = np.linspace(0,3,n*5)
    rbfi = interpolate.Rbf(xr, yr) #use default 'multiquadric'::
    #yi = rbfi(xi)
    plt.plot(xi,rbfi(xi),'r.-')
    rbfig = interpolate.Rbf(xr, yr, function='gaussian')
    plt.plot(xi,rbfig(xi),'g.-')
    plt.title('interpolation with global rbf, n = %d' % n + \
              '\nred: multiquadric, green: gaussian')



def windowedrbf(xi,x,d, hw=2):
    '''interpolate 1D points using local rbf

    Parameters
    ----------

    xi: 1D array
        points for which interpolated values are calculated
    x: 1D array
        coordinates of node points
    d: 1D array
        values, hight at coordinates x
    hw: integer (default = 2)
        half window size, number of points to include in the rbf
        in each direction.
        The number of included points in the rbf is 2*hw for points
        that are not nodes and 2*hw+1 for node points

    '''

    ind = np.argsort(x)
    x = x[ind]
    y = d[ind]
    ny = len(y)

    # find neighborhoods
    indx = np.argsort(xi)
    revindx = np.argsort(indx)  # inverse sort index 
    xs = xi[indx]
    intervals = np.searchsorted(x,xs)
    neigh = {}
    for ii,iv in enumerate(intervals):
        neigh.setdefault(iv,[]).append(ii)

    # interpolate for points in each neighborhood
    di = []       
    for k in neigh:
        sl = slice(max(0,k-hw),min(ny,k+hw),1)
        #print sl
        #print x[sl], y[sl]
        
        # get rbf for neighborhood
        rbfi = interpolate.Rbf(x[sl], y[sl])
        xiindlocal = neigh[k]
        #print xiindlocal
        xilocal = xs[xiindlocal]
        di.extend(zip(xilocal,rbfi(xilocal)))
    xdi = np.array(di)
    return xdi[revindx,:], neigh, intervals

def show_rbfw_example(basefn, n, gap=True, extrapol=False):
    '''example to plot comparison between global and local rbf

    Parameters
    ----------

    basefn: function
        sample points are created with basefn plus noise
    n: integer
        number of sample points to create
    gap: boolean (default: True)
        if true, create a 10% gap starting in middle of sample
    extrapol: boolean
        if true, extrapolate 5 points (1 node point distance)

    returns
    -------
    
    None
    output: two plots

    '''
    xr = np.linspace(0,3,n)
    yr = basefn(xr)+0.1*np.random.randn(xr.size)
    
    if gap:
        rm = slice(n/2,np.int(n*(1/2.0+0.1)))
        xr = np.delete(xr,rm)
        yr = np.delete(yr,rm)
    if extrapol:
        extrap = 3/float(n)  # extrapolate 5 points (1 nodepoint distance)
    else:
        extrap = 0
    xi = 3 + extrap - np.linspace(0,3+extrap,5*n)
    res, resn, iv = windowedrbf(xi,xr,yr)
    assert np.any(res[:,0]==xi), 'original order of xi array not preserved'
##    plt.figure()
##    plt.plot(res[:,0],res[:,1],'r.',xr,yr,'o')
##    plt.title('interprolation with local rbf, n = %d' % n)
    plt.figure()
    plt.plot(xr,yr,'o',res[:,0],res[:,1],'r.-')
    plt.title('interprolation with local rbf, n = %d' % n)
    plt.figure()
    rbfi = interpolate.Rbf(xr, yr) #use default 'multiquadric'::
    plt.plot(xi,rbfi(xi),'r.-',xr,yr,'o',)
    plt.title('interpolation with global rbf, n = %d' % n)

#basefn = lambda (x): np.sin(x**2)
def basefn1(x):
    return np.sin(x**2)
def basefn2(x):
    return x**2

# global rbf works ok up to 20 points
##for basefn in [basefn1, basefn2]:
##    for n in [20]:#[3,10,20,100]:
##        show_rbf_example(basefn,n)

for basefn in [basefn1, basefn2]:
    for n in [10,20,100]:
        show_rbfw_example(basefn,n)

#plt.show()
-------------- next part --------------
A non-text attachment was scrubbed...
Name: rbf100.png
Type: image/png
Size: 60550 bytes
Desc: not available
Url : http://mail.scipy.org/pipermail/scipy-user/attachments/20090403/a03ce14a/attachment-0001.png 


More information about the SciPy-user mailing list