Class: GSL::ScatterInterp

Inherits:
Object
  • Object
show all
Defined in:
lib/gsl_extras.rb

Constant Summary collapse

THIN_PLATE_SPLINES = 
:thin_plate_splines

Class Method Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(func, datavecs, norm = false, r0 = 1.0) ⇒ ScatterInterp

Returns a new instance of ScatterInterp.



121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
 
# 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

Class Method Details

.alloc(func, datavecs, norm, r0 = 1.0) ⇒ Object

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.



115
116
117
 
# File 'lib/gsl_extras.rb', line 115

def self.alloc(func, datavecs, norm, r0=1.0)
	new(func, datavecs, norm, r0)
end

Instance Method Details

#eval(*pars) ⇒ Object

Return the interpolated value for the given parameters.



209
210
211
212
213
214
215
216
217
 
# File 'lib/gsl_extras.rb', line 209

def eval(*pars)
	raise ArgumentError("wrong number of points") if pars.size != @dim
# 			p vals
	pars = GSL::Vector.alloc(pars)
	return @npoints.times.inject(0.0) do |sum, i|
# 			sum + function(radius(vals, @gridpoints.row(i)))*@weights[i]
		sum + function(pars, @gridpoints.row(i))*@weights[i]
	end
end

#function(vec1, vec2) ⇒ Object

Return the value of the interpolation kernel for the separation between the two given vectors. If linear was chosen this will just be the normalised distance between the two points.



176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
 
# File 'lib/gsl_extras.rb', line 176

def function(vec1, vec2)
	case @func
	when :linear
		return normalized_radius(vec1, vec2)
	when :cubic_alt
		return normalized_radius(vec1, vec2)**(1.5)
	when :thin_plate_splines
		return 0.0 if radius(vec1, vec2) == 0.0
		return normalized_radius(vec1, vec2)**2.0 * Math.log(normalized_radius(vec1, vec2))
	when :thin_plate_splines_alt
		rnorm = ((@r0.prod.abs)**(2.0/@r0.size)*(((vec1-vec2).square / @r0.square).sum + 1.0))
		return rnorm * Math.log(rnorm)
	when :multiquadratic
# 			return Math.sqrt(radius(vec1, vec2)**2 + Math.sqrt(@r0.square.sum))
		(@r0.prod.abs)**(1.0/@r0.size)*Math.sqrt(((vec1-vec2).square / @r0.square).sum + 1.0)
	when :inverse_multiquadratic
		1.0 / ((@r0.prod.abs)**(1.0/@r0.size)*Math.sqrt(((vec1-vec2).square / @r0.square).sum + 1.0))
	when :cubic
		((@r0.prod.abs)**(2.0/@r0.size)*(((vec1-vec2).square / @r0.square).sum + 1.0))**(1.5)
# 			invs = ((vec1-vec2).square + @r0.square).sqrt**(-1)
# 			invs.sum
# 				p @ro
# 			return 1.0 / Math.sqrt(radius(vec1, vec2)**2 + Math.sqrt(@r0.square.sum))
# 		when :inverse_multiquadratic
# # 				p @ro
# 			return 1.0 / Math.sqrt(radius(vec1, vec2)**2 + Math.sqrt(@r0.square.sum))
	else
		raise ArgumentError.new("Bad radial basis function: #{@func}")
	end
end

#gaussian_smooth_eval(*vals, sigma_vec) ⇒ Object

Evaluate the function,



221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
 
# File 'lib/gsl_extras.rb', line 221

def gaussian_smooth_eval(*vals, sigma_vec)
	npix = 7
	raise "npix must be odd" if npix%2==0
	case vals.size
	when 2
# 			delt0 = 3.0*0.999999*sigma_vec[0] / (npix-1)
# 			delt1 = 3.0*0.999999*sigma_vec[1] / (npix-1)
# 			sig3 = 3.0*sigma
		vals0 = GSL::Vector.linspace(vals[0] - 3.0* sigma_vec[0], vals[0] + 3.0* sigma_vec[0], npix)
		vals1 = GSL::Vector.linspace(vals[1] - 3.0* sigma_vec[1], vals[1] + 3.0* sigma_vec[1], npix)
		mat = GSL::Matrix.alloc(vals0.size, vals1.size)
		for i in 0...vals0.size
			for j in 0...vals1.size
				mat[i,j] = eval(vals0[i], vals1[j])
			end
		end
		mat.gaussian_smooth(*sigma_vec.to_a)
		cent = (npix - 1) / 2
		return mat[cent, cent]
		
	else
		raise 'Not supported for this number of dimensions yet'
	end
end

#normalized_radius(vec1, vec2) ⇒ Object

:nodoc:



163
164
165
166
167
168
169
170
171
172
 
# File 'lib/gsl_extras.rb', line 163

def normalized_radius(vec1, vec2)
	#Math.sqrt(((vec1 - vec2) / @r0).square.sum)
	case @r0
	when Numeric
		Math.sqrt(((vec1 - vec2) / @r0).square.sum)
	else 
		Math.sqrt(((vec1/@r0.to_gslv - vec2/@r0.to_gslv)).square.sum)
	end

end

#radius(vec1, vec2) ⇒ Object

:nodoc:



157
158
159
 
# File 'lib/gsl_extras.rb', line 157

def radius(vec1, vec2) 	# :nodoc: 
	Math.sqrt((vec1 -  vec2).square.sum)
end

#to_contour(grid_size = @npoints) ⇒ Object

Create a GSL::Contour object for making contours of the interpolated function. Only works for functions of 2 variables.



248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
 
# File 'lib/gsl_extras.rb', line 248

def to_contour(grid_size=@npoints)
	m = Matrix.alloc(grid_size, grid_size)
	raise TypeError("Must be 3d data") unless @gridpoints.shape[1] == 2
# 			p @gridpoints.shape
# 			p @gridpoints
	xmax, xmin = @gridpoints.col(0).max, @gridpoints.col(0).min
	ymax, ymin = @gridpoints.col(1).max, @gridpoints.col(1).min
	p 'x', xmax, 'y', ymax
	xvec = GSL::Vector.alloc((0...grid_size).to_a) * (xmax - xmin) / (grid_size - 1.0).to_f + xmin
	yvec = GSL::Vector.alloc((0...grid_size).to_a) * (ymax - ymin) / (grid_size - 1.0).to_f + ymin	
	p 'x', xvec.max, 'y', yvec.max

	for i in 0...grid_size
		for j in 0...grid_size
# 					p xvec[i], yvec[j]
			m[i,j] = eval(xvec[i], yvec[j])
		end
	end
	p 'm', m.max
	Contour.alloc(xvec, yvec, m)
end