class GSL::ScatterInterp

A class for interpolation from scattered multidimensional data using radial basis functions.

E.g.

x = GSL::Vector.alloc([0,3,6])
y = GSL::Vector.alloc([0,1,2])
z = GSL::Vector.alloc([0,5,10])
normalising_radius = GSL::Vector.alloc([3,1])

int = GSL::ScatterInterp.alloc(:linear, [x,y,z], false, normalising_radius)
puts int.eval(4.5, 1.7)

Constants

THIN_PLATE_SPLINES

Public Class Methods

alloc(func, datavecs, norm, r0=1.0) click to toggle source

Create a new interpolation class, for interpolating a function of one or more variables from a scattered dataset. datavecs is an array of vectors; the last vector should be the values of the function to be interpolated from, the other vectors are the coordinates or parameters corresponding to those values. Norm is a boolean; should the normalised functions be used? Default false. Func is the particular basis function to be used: can be

  • :linear

  • :cubic

  • :thin_plate_splines

  • :multiquadratic

  • :inverse_multiquadratic

r0 is a normalising radius which should be on the order of the scale length of the variation. If the scale length differs for each direction, an array of values can be passed; the length of this array should be one less than the number of datavecs, i.e. it should be the number of coordinates.

# File lib/gsl_extras.rb, line 115
def self.alloc(func, datavecs, norm, r0=1.0)
        new(func, datavecs, norm, r0)
end
new(func, datavecs, norm=false, r0=1.0) click to toggle source
# File lib/gsl_extras.rb, line 121
        def initialize(func, datavecs, norm=false, r0=1.0)
#               p datavecs
                 @norm = norm
                @npoints = datavecs[0].size
                datavecs.map!{|v| v.to_a}
                datavecs.each{|vec| raise ArgumentError.new("Datavectors must all have the same size ") unless vec.size == @npoints}
                @func = func
                data = datavecs.pop
                @gridpoints = GSL::Matrix.alloc(*datavecs).transpose
#                       puts @gridpoints.shape
                @dim = datavecs.size
                @r0 = r0
                if @r0.kind_of? Numeric
                        v = GSL::Vector.alloc(@dim)
                        v.set_all(@r0)
                        @r0 = v
                end
                m = GSL::Matrix.alloc(@npoints, @npoints)
                for i in 0...@npoints
                        for j in 0...@npoints
#                                       ep i, j
#                                       if true or i>= j
#                                       p @gridpoints.row(i), @gridpoints.row(j) if i == j
                                m[i,j] = function(@gridpoints.row(i), @gridpoints.row(j)) #(radius(@gridpoints.row(i), @gridpoints.row(j)))
#                                       else
#                                               m[i,j] = 0.0
#                                       end
                        end
#                               data[i] = data[i] * m.get_row(i).sum if norm
                end
#                       ep m
                @weights = m.LU_solve(GSL::Vector.alloc(data))
#                       ep @weights
        end