SG++
predictiveANOVARefinement.py

This example can be found under base/examples/predictiveANOVARefinement.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.createLinearGrid(dim)
19 HashGridStorage = grid.getStorage()
20 print "dimensionality: {}".format(dim)
21
22 # create regular grid, level 3
23 level = 3
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*10)+x1
31
32 # create coefficient vectors
33 alpha = DataVector(HashGridStorage.getSize())
34 print "length of alpha vector: {}".format(alpha.getSize())
35
36
37 #dataPoints
38
39 rows = 100
40 cols = 100
41
42 dataSet = DataMatrix(rows*cols,dim)
43 vals = DataVector(rows*cols)
44
45 # Create a "List" of points where the error should be calculated.
46 # This represents a regular 2d grid with a step size of 1 / rows and 1 / cols.
47 for i in xrange(rows):
48  for j in xrange(cols):
49  #xcoord
50  dataSet.set(i*cols+j,0,i*1.0/rows)
51  #ycoord
52  dataSet.set(i*cols+j,1,j*1.0/cols)
53  vals[i*cols+j] = f(i*1.0/rows,j*1.0/cols)
54
55 def calculateError(dataSet,f,grid,alpha,error):
56  print "calculating error"
57  #traverse dataSet
58  vec = DataVector(2)
59  opEval = createOperationEval(grid)
60  for i in xrange(dataSet.getNrows()):
61  dataSet.getRow(i,vec)
62  error[i] = pow(f(dataSet.get(i,0),dataSet.get(i,1))-opEval.eval(alpha,vec),2)
63  return error
64
65 #store old files
66 xCoordsOld = []
67 yCoordsOld = []
68
69 opEval = createOperationEval(grid)
70 for i in xrange(HashGridStorage.getSize()):
71  gridPointCoordinates = DataVector(dim)
72  HashGridStorage.getPoint(i).getStandardCoordinates(gridPointCoordinates)
73  xCoordsOld.append(gridPointCoordinates[0])
74  yCoordsOld.append(gridPointCoordinates[1])
75
76 # now refine adaptively 20 times
77 for refnum in range(20):
78  # set function values in alpha
79  for i in xrange(HashGridStorage.getSize()):
80  gp = HashGridStorage.getPoint(i)
81  alpha[i] = f(gp.getStandardCoordinate(0), gp.getStandardCoordinate(1))
82
83  # hierarchize
84  createOperationHierarchisation(grid).doHierarchisation(alpha)
85
86  #initialize plotter
87  plotter.hold(True)
88
89
90  xCoordinates = []
91  yCoordinates = []
92
93  #print all points
94
95  opEval = createOperationEval(grid)
96
97  for i in xrange(HashGridStorage.getSize()):
98  gridPointCoordinates = DataVector(dim)
99  HashGridStorage.getPoint(i).getStandardCoordinates(gridPointCoordinates)
100  xCoordinates.append(gridPointCoordinates[0])
101  yCoordinates.append(gridPointCoordinates[1])
102  bla = opEval.eval(alpha,gridPointCoordinates)
103
104  plotter.scatter(xCoordinates, yCoordinates, c='b')
105  plotter.scatter(xCoordsOld, yCoordsOld, c='r')
106  xCoordsOld = xCoordinates
107  yCoordsOld = yCoordinates
108
109  #show plot
110
111  plotter.hold(False)
112  plotter.show()
113
114  #calculate squared offset
115  errorVector = DataVector(dataSet.getNrows())
116  calculateError(dataSet, f, grid, alpha, errorVector)
117
118  #refinement stuff
119  refinement = HashRefinement()
120  decorator = PredictiveANOVARefinement(refinement)
121  # refine a single grid point each time
122  print "Error over all = %s" % errorVector.sum()
123  indicator = PredictiveRefinementIndicator(grid,dataSet,errorVector,1)
124  decorator.free_refine(HashGridStorage,indicator)
125
126  print "Refinement step %d, new grid size: %d" % (refnum+1, HashGridStorage.getSize())
127
128  #
129  #plot grid
130  #
131
132
133  # extend alpha vector (new entries uninitialized)
134  alpha.resize(HashGridStorage.getSize())