SG++
predictiveRefinement.py

This example can be found under base/examples/predictiveRefinement.py.

1 # Copyright (C) 2008-today The SG++ project
2 # This file is part of the SG++ project. For conditions of distribution and
3 # use, please see the copyright notice provided with SG++ or at
4 # sgpp.sparsegrids.org
5
6 #!/usr/bin/python
7
8 # import modules
9 import sys
10 import math
11
12 from pysgpp import *
13 import matplotlib.pyplot as plotter
14 from mpl_toolkits.mplot3d import Axes3D
15
16 # create a two-dimensional piecewise bi-linear grid
17 dim = 2
18 grid = Grid.createModLinearGrid(dim)
19 HashGridStorage = grid.getStorage()
20 print "dimensionality: {}".format(dim)
21
22 # create regular grid, level 3
23 level = 1
24 gridGen = grid.getGenerator()
25 gridGen.regular(level)
26 print "number of initial grid points: {}".format(HashGridStorage.getSize())
27
28 # definition of function to interpolate - nonsymmetric(!)
29 #f = lambda x0, x1: 16.0 * (x0-1)*x0 * (x1-1)*x1-x1
30 f = lambda x0, x1: math.sin(x0*math.pi)
31
32 # create coefficient vectors
33 alpha = DataVector(HashGridStorage.getSize())
34 print "length of alpha vector: {}".format(alpha.getSize())
35
36 #dataPoints
37
38 rows = 100
39 cols = 100
40
41 dataSet = DataMatrix(rows*cols,dim)
42 vals = DataVector(rows*cols)
43
44 # Create a "List" of points where the error should be calculated.
45 # This represents a regular 2d grid with a step size of 1 / rows and 1 / cols.
46 for i in xrange(rows):
47  for j in xrange(cols):
48  #xcoord
49  dataSet.set(i*cols+j,0,i*1.0/rows)
50  #ycoord
51  dataSet.set(i*cols+j,1,j*1.0/cols)
52  vals[i*cols+j] = f(i*1.0/rows,j*1.0/cols)
53
54 def calculateError(dataSet,f,grid,alpha,error):
55  print "calculating error"
56  #traverse dataSet
57  vec = DataVector(2)
58  opEval = createOperationEval(grid)
59  for i in xrange(dataSet.getNrows()):
60  dataSet.getRow(i,vec)
61  error[i] = pow(f(dataSet.get(i,0),dataSet.get(i,1))-opEval.eval(alpha,vec),2)
62  return error
63
64 #store old files
65 xCoordsOld = []
66 yCoordsOld = []
67
68 opEval = createOperationEval(grid)
69 for i in xrange(HashGridStorage.getSize()):
70  gridPointCoordinates = DataVector(dim)
71  HashGridStorage.getPoint(i).getStandardCoordinates(gridPointCoordinates)
72  xCoordsOld.append(gridPointCoordinates[0])
73  yCoordsOld.append(gridPointCoordinates[1])
74
75 # now refine adaptively 20 times
76 for refnum in range(20):
77  # set function values in alpha
78  for i in xrange(HashGridStorage.getSize()):
79  gp = HashGridStorage.getPoint(i)
80  alpha[i] = f(gp.getStandardCoordinate(0), gp.getStandardCoordinate(1))
81
82  # hierarchize
83  createOperationHierarchisation(grid).doHierarchisation(alpha)
84
85  #initialize plotter
86  plotter.hold(True)
87
88
89  xCoordinates = []
90  yCoordinates = []
91
92  #print all points
93
94  opEval = createOperationEval(grid)
95
96  for i in xrange(HashGridStorage.getSize()):
97  gridPointCoordinates = DataVector(dim)
98  HashGridStorage.getPoint(i).getStandardCoordinates(gridPointCoordinates)
99  xCoordinates.append(gridPointCoordinates[0])
100  yCoordinates.append(gridPointCoordinates[1])
101  bla = opEval.eval(alpha,gridPointCoordinates)
102
103  plotter.scatter(xCoordinates, yCoordinates, c='b')
104  plotter.scatter(xCoordsOld, yCoordsOld, c='r')
105  xCoordsOld = xCoordinates
106  yCoordsOld = yCoordinates
107
108  #show plot
109
110  plotter.hold(False)
111  plotter.show()
112
113  #calculate squared offset
114  errorVector = DataVector(dataSet.getNrows())
115  calculateError(dataSet, f, grid, alpha, errorVector)
116
117  #refinement stuff
118  refinement = HashRefinement()
119  decorator = PredictiveRefinement(refinement)
120  # refine a single grid point each time
121  print "Error over all = %s" % errorVector.sum()
122  indicator = PredictiveRefinementIndicator(grid,dataSet,errorVector,1)
123  decorator.free_refine(HashGridStorage,indicator)
124
125  print "Refinement step %d, new grid size: %d" % (refnum+1, HashGridStorage.getSize())
126
127  #
128  #plot grid
129  #
130
131
132  # extend alpha vector (new entries uninitialized)
133  alpha.resize(HashGridStorage.getSize())