SG++-Doxygen-Documentation
LTwoDotTest.py

This example can be found under pde/examples/LTwoDotTest.py.

1 import pysgpp
2 from pysgpp.extensions.datadriven.uq.operations.sparse_grid import getBasis
3 import numpy as np
4 import random
5 import matplotlib.pyplot as plt
6 
7 def test_base():
8  n = 100000
9  b = pysgpp.SLinearModifiedBase()
10  # print(sum([b.eval(0,0,x) for x in np.linspace(0, 1, n)]) / n)
11  for l in range(1,3,1):
12  print(l)
13  for i in range(1, 2**l, 1):
14  s = 1
15  s = sum([b.eval(l,i,x) for x in np.linspace(0, 1, n)]) / n
16  print(("Level:",l , "index:", i))
17  epsilon = abs(s - b.getIntegral(l,i))
18  if(epsilon > 10**-6):
19  print("s: {} res: {}".format(s, b.getIntegral(l,i)))
20  print("error: {}".format(epsilon))
21 
22 def test_LTwoDot(grid, l):
23  res = 10000
24  b = grid.getBasis();
25  grid.getGenerator().regular(l)
26  gridStorage = grid.getStorage()
27  size = gridStorage.getSize()
28  # print(size)
29  m = pysgpp.DataMatrix(size, size)
30  opMatrix = pysgpp.createOperationLTwoDotExplicit(m, grid)
31  # print m
32  # m_comp = pysgpp.DataMatrix(size, size)
33  for i in range(gridStorage.getSize()):
34  for j in range(i, gridStorage.getSize()):
35  gpi = gridStorage.getPoint(i)
36  gpj = gridStorage.getPoint(j)
37  sol = 1
38  # print "--------"
39  # print "i:{} j:{}".format(i, j)
40  for k in range(d):
41  lik = gpi.getLevel(k)
42  iik = gpi.getIndex(k)
43  ljk = gpj.getLevel(k)
44  ijk = gpj.getIndex(k)
45  # print "i l,i: {},{} j l,i: {},{}".format(lik, iik, ljk, ijk)
46  xs = np.linspace(0, 1, res)
47  tmp = sum([b.eval(lik, iik, x)*b.eval(ljk, ijk, x) for x in xs]) / res
48  sol *= tmp
49  # print("lik:{} iik:{} ljk:{} ijk:{} k:{} tmp: {}".format(lik, iik, ljk, ijk, k,tmp))
50  # print(sol)
51  error = abs(m.get(i,j) - sol)
52  # print error
53  if(error >= 10**-4):
54  print("i:{} j:{} error: {}".format(i, j, error))
55  print("iik:{} lik:{} ijk:{} ljk:{} error: {}".format(iik, lik, ijk, ljk, error))
56  print("is:{} should:{}".format(m.get(i,j), sol))
57 
58 def test_LTwoDotImplicit(grid, l):
59  grid.getGenerator().regular(l)
60  gridStorage = grid.getStorage()
61  size = gridStorage.getSize()
62  # print(size)
63  m = pysgpp.DataMatrix(size, size)
64  opExplicit = pysgpp.createOperationLTwoDotExplicit(m, grid)
65  op = pysgpp.createOperationLTwoDotProduct(grid)
66  alpha = pysgpp.DataVector(size)
67  resultExplicit = pysgpp.DataVector(size)
68  result = pysgpp.DataVector(size)
69  for i in range(size):
70  alpha[i] = 1
71  opExplicit.mult(alpha, resultExplicit)
72  op.mult(alpha, result)
73  for i in range(size):
74  if result[i] != resultExplicit[i]:
75  print("Error result entry {} differs".format(i))
76 
77  if abs(result[i] - resultExplicit[i]) > 1e-16:
78  # print result[i] - resultExplicit[i]
79  print("result:{}".format(result[i]))
80  print("resultExplicit:{}".format(resultExplicit[i]))
81 
82 def test_laplace(grid, lmax):
83  resolution = 100000
84  grid.getGenerator().regular(lmax)
85  gridStorage = grid.getStorage()
86  size = gridStorage.getSize()
87  b = getBasis(grid)
88  op = pysgpp.createOperationLaplace(grid)
89  alpha = pysgpp.DataVector(size)
90  result = pysgpp.DataVector(size)
91 
92  for point_i in range(size):
93  for point_j in range(size):
94  gp_i = gridStorage.getPoint(point_i)
95  gp_j = gridStorage.getPoint(point_j)
96  print("--------")
97  for i in range(0, size):
98  alpha[i] = 0
99  alpha[point_i] = 1
100  op.mult(alpha, result)
101  xs = np.linspace(0, 1, resolution)
102  approx = sum([b.evalDx(gp_i.getLevel(0), gp_i.getIndex(0), x) * b.evalDx(gp_j.getLevel(0), gp_j.getIndex(0), x) for x in xs]) / resolution
103  print("i,j: {},{} result: {} approx:{}".format(point_i, point_j, result[point_j], approx))
104  if(abs(result.get(point_j) - approx) > 1e-1):
105  print("--------")
106  print("points: {},{} ".format(point_i, point_j))
107  print("approx:{}".format(approx))
108  print("result:{}".format(result.get(point_j)))
109  # print result
110  print("--------")
111 
112 def test_poly_evaldx():
113  l = 3
114  i = 1
115  x = 0.12
116  eps = 0.0001
117  b = pysgpp.SPolyModifiedClenshawCurtisBase(3)
118  tang = b.evalDx(l, i, x)
119  sec = (b.eval(l, i, x + eps) - b.eval(l, i, x - eps)) / (2*eps)
120  print("evalDx:{}".format(tang))
121  print("sekante:{}".format(sec))
122  print("evals: {} {}".format( b.eval(l, i, x - eps), b.eval(l, i, x + eps) ))
123 
124 def plot_evaldx():
125  l = 3
126  i = 1
127  xs = np.linspace(0, 1, 50)
128  b = pysgpp.SPolyClenshawCurtisBase(3)
129  plt.plot(xs, [b.evalDx(l, i, x) for x in xs])
130  plt.show()
131 
132 def plot_evaldx_prod(grid, lmax, i, j):
133  grid.getGenerator().regular(lmax)
134  gridStorage = grid.getStorage()
135  b = getBasis(grid)
136  gp_i = gridStorage.getPoint(i)
137  gp_j = gridStorage.getPoint(j)
138  lik = gp_i.getLevel(0)
139  ljk = gp_j.getLevel(0)
140  iik = gp_i.getIndex(0)
141  ijk = gp_j.getIndex(0)
142  xs = np.linspace(0, 1, 100000)
143 
144  plt.plot(xs, [b.evalDx(lik, iik, x) * b.evalDx(ljk, ijk, x) for x in xs])
145  plt.show()
146 
147 
148 def test_laplace2(grid, lmax):
149  # grid.getGenerator().regular(lmax)
150  gridStorage = grid.getStorage()
151  size = gridStorage.getSize()
152  m = pysgpp.DataMatrix(size, size)
153  op = pysgpp.createOperationLaplaceExplicit(m, grid)
154  print(m)
155 
156 
157 
158 def test_firstMoment(grid, lmax):
159  grid.getGenerator().regular(lmax)
160  resolution = 100000
161  gridStorage = grid.getStorage()
162  b = grid.getBasis()
163  op = pysgpp.createOperationFirstMoment(grid)
164  alpha = pysgpp.DataVector(grid.getSize(), 1.0)
165  bounds = pysgpp.DataMatrix(1, 2, 0.0)
166  bounds.set(0, 1, 1.0)
167  res = 0.0
168  for i in range(grid.getSize()):
169  lev = gridStorage.getPoint(i).getLevel(0)
170  ind = gridStorage.getPoint(i).getIndex(0)
171  temp_res = 0.0
172  for c in range(resolution):
173  x = float(c) / resolution
174  temp_res += x * b.eval(lev, ind, x)
175  res += alpha.get(i) * temp_res / resolution
176  print("--FirstMoment--")
177  print(res)
178  print(op.doQuadrature(alpha, bounds))
179  print (res - op.doQuadrature(alpha, bounds))
180 
181 def test_secondtMoment(grid, lmax):
182  gridStorage = grid.getStorage()
183  gridStorage.clear()
184  grid.getGenerator().regular(lmax)
185  resolution = 1000000
186  b = grid.getBasis()
187  op = pysgpp.createOperationSecondMoment(grid)
188  alpha = pysgpp.DataVector(grid.getSize(), 1.0)
189  bounds = pysgpp.DataMatrix(1, 2, 0.0)
190  bounds.set(0, 1, 1.0)
191  res = 0.0
192  for i in range(grid.getSize()):
193  lev = gridStorage.getPoint(i).getLevel(0)
194  ind = gridStorage.getPoint(i).getIndex(0)
195  temp_res = 0.0
196  for c in range(resolution):
197  x = float(c) / resolution
198  temp_res += x * x *b.eval(lev, ind, x)
199  res += alpha.get(i) * temp_res / resolution
200 
201  print("--SecondMoment--")
202  print(res)
203  print(op.doQuadrature(alpha, bounds))
204  print (res - op.doQuadrature(alpha, bounds))
205 
206 # test_poly_evaldx()
207 # plot_evaldx()
208 # test_base()
209 d = 1
210 l = 3
211 grid = pysgpp.Grid.createModLinearGrid(d)
212 # plot_evaldx_prod(grid, 4, 1, 4)
213 # test_laplace(grid, l)
214 # test_laplace2(grid, l)
215 # test_LTwoDot(grid, l)
216 # test_LTwoDotImplicit(grid, l)
217 test_firstMoment(grid, l)
218 test_secondtMoment(grid, l)