eegkit / py_tsg Goto Github PK

# EXAMPLE 1:
# make a classical Smolyak grid using Clenshaw-Curtis quadrature
# integrate the function f(x,y) = exp( -x^2 ) * cos( y )
# the exact answer is: 2.513723354063905e+00

import scipy, scipy.integrate, itertools
import pyublas
import _py_tsg as tsg

grid = tsg.TSG()
def fn1(x):		return scipy.exp(-x[0]*x[0])*scipy.cos(x[1])
dimension = 2
outputs = 0
level = 7
grid.make_global_grid( dimension, outputs, level, tsg.TypeDepth.type_level, tsg.TypeOneDRule.rule_clenshawcurtis, scipy.array([], dtype=int), 0, 0 )
points = grid.get_points()
weights = grid.get_weights()
sum = scipy.sum( [w*fn1(x) for (x,w) in zip(points, weights)] )
print("\nExample 1")
print("grid has: %d points" % grid.get_num_points())
print("integral is: %.17f" % sum)
print("error: %.17f" % scipy.fabs( sum - 2.513723354063905e+00 ))

# make a Clenshaw-Curtis rule that interpolates exactly polynomials of order up to 10
# integrate the function f(x,y) = exp( -x^2 ) * cos( y )
# exact values is: 6.990131267703512e-01	

import scipy, scipy.integrate, itertools
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

import pyublas
import _py_tsg as tsg

def apply_fn_to_grid(gridList, fn):
	xy_list = list(itertools.product(*gridList))
	xlists = []	
	[xlists.append( [xy[i] for xy in xy_list] ) for i in range(len(gridList))]
	f_list = [fn(xy) for xy in xy_list]
	return (xlists, f_list)

def plot_gridpoints_2(grid):
	assert(grid.get_num_dimensions() == 2)
	points = grid.get_points()
	fig = plt.figure()
	plt.scatter(points[:,0], points[:,1], s=3)
	print("%d gridpoints" % len(points))

def scatter_grid_2(grid, n_gridpoints=20):
	assert(grid.get_num_dimensions() == 2)
	(low, high) = grid.get_transform_AB()	
	grid1 = scipy.linspace(low[0], high[0], n_gridpoints)
	grid2 = scipy.linspace(low[1], high[1], n_gridpoints)
	(xlists, f_list) = apply_fn_to_grid([grid1, grid2], lambda x: grid.evaluate(scipy.array(x)))
	fig = plt.figure()
	ax = Axes3D(fig)	
	ax.scatter3D(xlists[0], xlists[1], f_list, s=3)	
	return ax

grid = tsg.TSG()
def fn1(x):		return scipy.exp(-x[0]*x[0])*scipy.cos(x[1])
dimension = 2
outputs = 1
precision = 10
grid.make_global_grid( dimension, outputs, precision, tsg.TypeDepth.type_basis, tsg.TypeOneDRule.rule_clenshawcurtis, scipy.array([], dtype=int), 0, 0 )
points = grid.get_needed_points()
values = fn1(points.T)
grid.load_needed_points(values)
plot_gridpoints_2(grid)
scatter_grid_2(grid, 50)