Base quick start (Python)

To be able to quickly start with a toolkit, it is often advantageous (not only for the impatient users), to look at some code examples first. In this tutorial, we give a short example program which interpolates a bivariate function on a regular sparse grid. Identical versions of the example are given in all languages currently supported by SG++: C++, Python, Java, and MATLAB.

In the example, we create a two-dimensional regular sparse grid of level 3 (with grid points ) using piecewise bilinear basis functions . We then interpolate the function

with

by calculating the coefficients such that for all . This process is called hierarchization in sparse grid contexts; the are called (hierarchical) surpluses. Note that vanishes at the boundary of the domain ; therefore, we don't have to spend sparse grid points on the boundary. Finally, we evaluate the sparse grid function at a point .

At the beginning of the program, we have to import the pysgpp library.

import pysgpp

Before starting, the function , which we want to interpolate, is defined.

f = lambda x0, x1: 16.0 * (x0 - 1.0) * x0 * (x1 - 1.0) * x1

First, we create a two-dimensional grid (type sgpp::base::Grid) with piecewise bilinear basis functions with the help of the factory method sgpp::base::Grid.createLinearGrid().

dim = 2
grid = pysgpp.Grid.createLinearGrid(dim)

Then we obtain the grid's sgpp::base::GridStorage object which allows us, e.g., to access grid points, to obtain the dimensionality (which we print) and the number of grid points.

gridStorage = grid.getStorage()
print("dimensionality:         {}".format(gridStorage.getDimension()))

Now, we use a sgpp::base::GridGenerator to create a regular sparse grid of level 3. Thus, gridStorage.getSize() returns 17, the number of grid points of a two-dimensional regular sparse grid of level 3.

level = 3
grid.getGenerator().regular(level)
print("number of grid points:  {}".format(gridStorage.getSize()))

We create an object of type sgpp::base::DataVector which is essentially a wrapper around a double array. The DataVector is initialized with as many entries as there are grid points. It serves as a coefficient vector for the sparse grid interpolant we want to construct. As the entries of a freshly created DataVector are not initialized, we set them to 0.0. (This is superfluous here as we initialize them in the next few lines anyway.)

alpha = pysgpp.DataVector(gridStorage.getSize())
alpha.setAll(0.0)
print("length of alpha vector: {}".format(len(alpha)))

The for loop iterates over all grid points: For each grid point gp, the corresponding coefficient is set to the function value at the grid point's coordinates which are obtained by getStandardCoordinate(dim). The current coefficient vector is then printed.

for i in range(gridStorage.getSize()):
  gp = gridStorage.getPoint(i)
  alpha[i] = f(gp.getStandardCoordinate(0), gp.getStandardCoordinate(1))

print("alpha before hierarchization: {}".format(alpha))

An object of sgpp::base::OperationHierarchisation is created and used to hierarchize the coefficient vector, which we print.

pysgpp.createOperationHierarchisation(grid).doHierarchisation(alpha)
print("alpha after hierarchization:  {}".format(alpha))

Finally, a second DataVector is created which is used as a point to evaluate the sparse grid function at. An object is obtained which provides an evaluation operation (of type sgpp::base::OperationEvaluation), and the sparse grid interpolant is evaluated at , which is close to (but not exactly at) a grid point.

p = pysgpp.DataVector(dim)
p[0] = 0.52
p[1] = 0.73
opEval = pysgpp.createOperationEval(grid)
print("u(0.52, 0.73) = {}".format(opEval.eval(alpha, p)))

The example results in the following output:

dimensionality:         2
number of grid points:  17
length of alpha vector: 17
alpha before hierarchization: [1, 0.75, 0.75, 0.4375, 0.9375, 0.9375, 0.4375, 0.75, 0.75, 0.4375, 0.9375, 0.9375, 0.4375, 0.5625, 0.5625, 0.5625, 0.5625]
alpha after hierarchization:  [1, 0.25, 0.25, 0.0625, 0.0625, 0.0625, 0.0625, 0.25, 0.25, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625]
u(0.52, 0.73) = 0.7696

It can be clearly seen that the surpluses decay with a factor of 1/4: On the first level, we obtain 1, on the second 1/4, and on the third 1/16 as surpluses.