sparseSpACE for interpolation

[OlssonA]

Hi,

I want to use sparseSpACE to benchmark how the SG combination technique performs for interpolation of my 5D function that I cannot evaluate at the boundaries. I have already had positive results with spatially adaptive SGs with hierarchical construction, and I want to see how this method compares.

I tried using the Regression module but I found it difficult to understand how many points were used, and if all of them were used for the grid and not validation etc.

Is there a way to construct an interpolation without using points at the boundaries (with some extrapolation approach instead)? I want to evaluate the interpolant at some data points of my choosing and construct my own error norms, and avoid using any built in error estimators.

Thanks!
/Anton

[obersteiner]

Hi Anton,
Interpolation can be just done using the Integration operation as in the tutorial:

%matplotlib inline
import sparseSpACE
import numpy as np
from sparseSpACE.spatiallyAdaptiveSingleDimension2 import *
from sparseSpACE.Function import *
from sparseSpACE.ErrorCalculator import *

# dimension of the problem
dim = 5

# define integration domain boundaries
a = np.zeros(dim)
b = np.ones(dim)
# toDo use your boundaries here

# define function to be integrated
coefficients = np.array([ float(10**1) * (d+1) for d in range(dim)])
f = GenzProductPeak(midpoint=midpoint,coefficients=coefficients)
# toDo use your function here

# plot function
# f.plot(np.ones(dim)*a,np.ones(dim)*b)

# reference integral solution for calculating errors
# reference_solution = f.getAnalyticSolutionIntegral(a,b)
reference_solution = None

# define error estimator for refinement
errorOperator = ErrorCalculatorSingleDimVolumeGuided()

# define equidistant grid
grid=GlobalTrapezoidalGrid(a=a, b=b, modified_basis=False, boundary=False)

# NEW! define operation which shall be performed in the combination technique
from sparseSpACE.GridOperation import *
operation = Integration(f=f, grid=grid, dim=dim, reference_solution=reference_solution)

# define SingleDim refinement strategy for Spatially Adaptive Combination Technique
adaptiveCombiInstanceSingleDim = SpatiallyAdaptiveSingleDimensions2(np.ones(dim) * a, np.ones(dim) * b, operation=operation)

# refinement stops once this threshold of nodes is exceeded
max_evaluations = 10000

# once error estimate is below tol the refinement stops; I would use max_evaluations if no refernece solution is given
tol = 10**-2

# performing the spatially adaptive refinement with the SingleDim method
adaptiveCombiInstanceSingleDim.performSpatiallyAdaptiv(1, 2, errorOperator, tol = tol, max_evaluations=max_evaluations)

print("Number of points used in refinement:", adaptiveCombiInstanceSingleDim.get_total_num_points())

interpolation_points = np.array([(0,)*dim, (1,)*dim])
# toDo use your interpolation points
print("Interpolating at points", interpolation_points)

 # perform interpolation 
results = adaptiveCombiInstanceSingleDim(interpolation_points)

print("Interpolated values", results)

I switched the boolean flag for the boundary in the grid to off. Therefore no boundary points are included in the grid. Instead currently 0 is assumed at the boundary and the interpolation is done accordingly. It is not yet fully implemented to use the interpolation for the modified basis. If you need this extrapolation towards the boundary, I could look into this.

Best,
Michael

[OlssonA]

Hi,

thank you Michael for the quick response. I do indeed need a modified extrapolating basis as my function does not vanish at the boundary. The linear extrapolation in https://doi.org/10.1016/j.jco.2010.04.001 has been working well for my spatially adaptive SG. If you already use the rescaled hat functions as basis functions, this extension is pretty straightforward.

My function does however vanish at the boundaries of two of the variables, so I will give the current version a try for a 2D slice in these parameters.

Thanks and best,
Anton

[obersteiner]

Hi Anton,

I have implemented the modified extrapolation basis. I needed to restructure the code a bit and find a workaround for some issues. It should work now as expected. You can find the code in the branch modified_basis_interpolation. It is currently only possible to use the modified basis in every dimension or to use 0 boundaries in all dimensions (no mix between both).
If you encounter any issues, please give me an update.

Best,
Michael