SG++-Doxygen-Documentation
test_Rosenblatt.py

This example can be found under datadriven/examples/test_Rosenblatt.py.

1 import pysgpp
2 import numpy as np
3 import matplotlib.pyplot as plt
4 import math
5 from pysgpp.extensions.datadriven.uq.plot.plot2d import plotSG2d
6 
7 pi = math.pi
8 
9 class interpolation_function(object):
10  def __init__(self, d, f):
11  self.f = f
12  self.d = d
13  self.grid = pysgpp.Grid.createBsplineClenshawCurtisGrid(d, 3)
14  self.gridStorage = self.grid.getStorage()
15  try :
16  self.hierarch = pysgpp.createOperationHierarchisation(self.grid)
17  except :
18  self.hierarch = pysgpp.createOperationMultipleHierarchisation(self.grid)
19  self.opeval = pysgpp.createOperationEvalNaive(self.grid)
20  self.alpha = pysgpp.DataVector(self.gridStorage.getSize())
21 
22  def create_interpolation(self, grid_lvl):
23  self.gridStorage.clear()
24  self.grid.getGenerator().regular(grid_lvl)
25  self.alpha = pysgpp.DataVector(self.gridStorage.getSize())
26  self.min_f = float('inf')
27  for i in range(self.gridStorage.getSize()):
28  gp = self.gridStorage.getPoint(i)
29  x = [self.gridStorage.getCoordinate(gp, j) for j in range(self.d)]
30  self.alpha[i] = self.f(x)
31  if self.alpha[i] < self.min_f:
32  self.min_f = self.alpha[i]
33  self.min_x = x
34  self.hierarch.doHierarchisation(self.alpha)
35 
36  def __call__(self, x):
37  if (self.d == 1 and not isinstance(x, list)):
38  x = [x]
39  return self.opeval.eval(self.alpha, pysgpp.DataVector(x))
40 
41 # ---------------------------------------------
42 
43 def distrib(x):
44  if isinstance(x, list):
45  x = x[0]
46  return (np.sin(2.5*pi*x - (pi - pi / 4.)) + 1. / 2**0.5) / 0.527044
47 
48 def parabola(x):
49  res = 1.
50  for i in range(len(x)):
51  res *= x[i] * (1. - x[i]) * 4.;
52  return res
53 
54 def const_0(x):
55  return 0.0
56 
57 def eval_rosenblatt1d(sg_pdf, xs):
58  op = pysgpp.createOperationRosenblattTransformation1D(sg_pdf.grid)
59  ys = []
60  for i,x in enumerate(xs):
61  print("---------------{}---------------".format(i))
62  ys.append(op.doTransformation1D(sg_pdf.alpha, x))
63  return ys
64  # ---------------------------------------------
65 
66 def eval_rosenblattdd(sg_pdf, xs):
67  op = pysgpp.createOperationRosenblattTransformation(sg_pdf.grid)
68  X, Y = np.meshgrid(xs, xs)
69  input_points = pysgpp.DataMatrix(len(xs), sg_pdf.d)
70  output_points = pysgpp.DataMatrix(len(xs), sg_pdf.d)
71  for i in range(len(xs)):
72  for j in range(sg_pdf.d):
73  input_points.set(i, j , 0.5)
74  op.doTransformation(sg_pdf.alpha, input_points, output_points)
75  print(output_points)
76 
77  # ---------------------------------------------
78 
79 
80 def eval_inverse_rosenblatt1d(sg_pdf, xs):
81  op = pysgpp.createOperationInverseRosenblattTransformation1D(sg_pdf.grid)
82  ys = []
83  print(op.doTransformation1D(sg_pdf.alpha, 0.999942))
84  for i,x in enumerate(xs):
85  # print("---------------{}---------------".format(i))
86  ys.append(op.doTransformation1D(sg_pdf.alpha, x))
87  return ys
88 
89 xs = np.arange(0., 1.01, 0.01)
90 l_max = 2
91 d = 2
92 interpolation = interpolation_function(d, parabola)
93 interpolation.create_interpolation(l_max)
94 
95 # test()
96 # grid_points = np.arange(0, 1.01, 2**-l_max)
97 # ys = [interpolation(x) for x in xs]
98 
99 # plotSG2d(interpolation.grid, interpolation.alpha)
100 # ys = eval_inverse_rosenblatt1d(interpolation, xs)
101 # ys = eval_rosenblatt1d(interpolation, xs)
102 # print(ys)
103 eval_rosenblattdd(interpolation, xs)
104 # plt.plot(xs, ys)
105 # plt.scatter(grid_points, np.zeros_like(grid_points))
106 # plt.pcolormesh(X, Y, Z)
107 # plt.show()