SG++-Doxygen-Documentation

simple code that provides convergence plots for various analytic models

1 from argparse import ArgumentParser
2 from pysgpp.extensions.datadriven.uq.parameters.ParameterBuilder import ParameterBuilder
3 from pysgpp.extensions.datadriven.uq.plot.colors import insert_legend
4 from pysgpp.extensions.datadriven.uq.plot.plot1d import plotFunction1d
5 from pysgpp.pysgpp_swig import DataVector, CombigridOperation,\
6  CombigridMultiOperation, CombigridTensorOperation
7 import pysgpp
8 
9 import matplotlib.pyplot as plt
10 import numpy as np
11 from pysgpp.extensions.datadriven.uq.dists import Uniform, Beta, TLognormal
12 
13 
14 def expModel(x, params):
15  return np.exp(-x[0] * x[1])
16 
17 
18 def arctanModel(x, params):
19  return np.arctan(50.0 * (x[0] - .35)) + np.pi / 2.0 + 4.0 * x[1] ** 3 + np.exp(x[0] * x[1] - 1.0)
20 
21 
22 def buildAtanParams(dist_type):
23  dist = generateDistribution(dist_type, [0, 1])
24 
25  parameterBuilder = ParameterBuilder()
26  up = parameterBuilder.defineUncertainParameters()
27 
28  up.new().isCalled("x1").withDistribution(dist)
29  up.new().isCalled("x2").withDistribution(dist)
30 
31  return parameterBuilder.andGetResult()
32 
33 
34 def boreholeModel(x, params):
35  z = params.getJointTransformation().unitToProbabilistic(x)
36 
37  num = 2 * np.pi * z[2] * (z[3] - z[5])
38  den = np.log(z[1] / z[0]) * (1 + (2 * z[6] * z[2]) / (np.log(z[1] / z[0]) * (z[0] ** 2) * z[7]) + z[2] / z[4])
39  return num / den
40 
41 
42 def buildBoreholeParams(dist_type):
43  boreholeLimits = [[0.05, 0.15],
44  [100, 50000],
45  [63070, 115600],
46  [990, 1110],
47  [63.1, 116],
48  [700, 820],
49  [1120, 1680],
50  [9855, 12045]]
51  boreholeParamNames = ["r_w", "r", "T_u", "H_u", "T_l", "H_l", "L", "K_w"]
52 
53  parameterBuilder = ParameterBuilder()
54  up = parameterBuilder.defineUncertainParameters()
55 
56  for k in range(8):
57  xlim = boreholeLimits[k]
58  dist = generateDistribution(dist_type, xlim)
59  up.new().isCalled(boreholeParamNames[k]).withDistribution(dist)
60 
61  return parameterBuilder.andGetResult()
62 
63 
64 def sin2Model(x, params):
65  return np.sin(x ** 2)
66 
67 
68 def buildSin2Params(dist_type):
69  dist = generateDistribution(dist_type, alpha=0.01)
70 
71  parameterBuilder = ParameterBuilder()
72  up = parameterBuilder.defineUncertainParameters()
73  up.new().isCalled("x1").withDistribution(dist)
74 
75  return parameterBuilder.andGetResult()
76 
77 
78 def generateDistribution(dist_type, xlim=None, alpha=None):
79  if dist_type == "uniform":
80  return Uniform(xlim[0], xlim[1])
81  elif dist_type == "beta":
82  return Beta(5, 4, xlim[0], xlim[1] - xlim[0])
83  elif dist_type == "lognormal":
84  return TLognormal.by_alpha(1e-12, np.exp(-1), alpha=alpha)
85  else:
86  raise AttributeError("dist type '%s' unknown" % dist_type)
87 
88 
89 def buildModel(name, dist_type):
90  if name == "arctan":
91  params = buildAtanParams(dist_type)
92  model = arctanModel
93  elif name == "borehole":
94  params = buildBoreholeParams(dist_type)
95  model = boreholeModel
96  elif name == "exp":
97  params = buildAtanParams(dist_type)
98  model = expModel
99  elif name == "sin2":
100  params = buildSin2Params(dist_type)
101  model = sin2Model
102  else:
103  raise AttributeError("model unknown")
104 
105  return model, params
106 
107 
108 def buildSparseGrid(gridType, basisType, degree=5, growthFactor=2, orthogonal_basis=None):
109  if orthogonal_basis is None:
110  if gridType == "ClenshawCurtis" and basisType == "bspline":
111  return CombigridMultiOperation.createExpClenshawCurtisBsplineInterpolation(numDims, func, degree)
112  elif gridType == "UniformBoundary" and basisType == "bspline":
113  return CombigridMultiOperation.createExpUniformBoundaryBsplineInterpolation(numDims, func, degree)
114  elif gridType == "Leja" and basisType == "poly":
115  return CombigridMultiOperation.createExpLejaPolynomialInterpolation(numDims, func)
116  elif gridType == "L2Leja" and basisType == "poly":
117  return CombigridMultiOperation.createExpL2LejaPolynomialInterpolation(numDims, func)
118  elif gridType == "ClenshawCurtis" and basisType == "poly":
119  return CombigridMultiOperation.createExpClenshawCurtisPolynomialInterpolation(numDims, func)
120  else:
121  raise AttributeError("not supported")
122  else:
123  if gridType == "Leja" and basisType == "poly":
124  return CombigridTensorOperation.createExpLejaPolynomialInterpolation(orthogonal_basis, numDims, func)
125  elif gridType == "L2Leja" and basisType == "poly":
126  return CombigridTensorOperation.createExpL2LejaPolynomialInterpolation(orthogonal_basis, numDims, func)
127  elif gridType == "ClenshawCurtis" and basisType == "poly":
128  return CombigridTensorOperation.createExpClenshawCurtisPolynomialInterpolation(orthogonal_basis, numDims, func)
129  else:
130  raise AttributeError("not supported")
131 
132 
133 def buildLevelManager(name):
134  if name == "regular":
135  return pysgpp.RegularLevelManager()
136  elif name == "averaging":
137  return pysgpp.AveragingLevelManager()
138  elif name == "weightedRatio":
139  return pysgpp.WeightedRatioLevelManager()
140  elif name == "variance":
141  return pysgpp.AveragingLevelManager()
142  else:
143  raise AttributeError("level manager '%s' not supported" % name)
144 
145 
146 class RefinementWrapper(object):
147 
148  def __init__(self,
149  gridType,
150  basisType,
151  degree,
152  growthFactor,
153  levelManagerType,
154  distType,
155  samples,
156  params):
157  self.operation = buildSparseGrid(gridType,
158  basisType,
159  degree,
160  growthFactor)
161 
162  # set evaluation points
163  self.numDims = samples.shape[1]
164  self.samples_mat = pysgpp.DataMatrix(samples)
165  self.samples_mat.transpose()
166  self.operation.setParameters(self.samples_mat)
167 
168  self.levelManagerType = levelManagerType
169  self.levelManager = buildLevelManager(levelManagerType)
170  if levelManagerType == "variance":
171  self.orthogonal_basis = params.getUnivariateOrthogonalPolynomials(dtype="orthogonal")[0]
172  self.tensor_operation = buildSparseGrid(gridType,
173  basisType,
174  degree,
175  growthFactor,
176  orthogonal_basis=self.orthogonal_basis)
177  print( self.tensor_operation.getLevelManager() )
178  self.tensor_operation.getLevelManager().addRegularLevels(1)
179  self.tensor_operation.setLevelManager(self.levelManager)
180  else:
181  self.operation.setLevelManager(self.levelManager)
182 
183  def evaluate(self, n):
184  if self.levelManagerType == "regular":
185  self.operation.getLevelManager().addRegularLevels(n)
186  elif self.levelManagerType == "variance":
187  numPoints = self.tensor_operation.getLevelManager().numGridPoints()
188  self.tensor_operation.getLevelManager().addLevelsAdaptive(np.max([n - numPoints, 0]))
189 
190 # tensorResult = self.tensor_operation.getResult()
191 # print "Total function evaluations: %i" % self.tensor_operation.numGridPoints()
192 # print "E(u) = %g" % tensorResult.get(pysgpp.IndexVector(self.numDims, 0)).getValue()
193 # print "Var(u) = %g" % tensorResult.norm() ** 2
194 
195  levelStructure = self.tensor_operation.getLevelManager().getLevelStructure()
196  self.operation.getLevelManager().addLevelsFromStructure(levelStructure)
197  else: # adaptive
198  numPoints = self.operation.getLevelManager().numGridPoints()
199  self.operation.getLevelManager().addLevelsAdaptive(np.max([n - numPoints, 0]))
200 
201  return self.operation.getResult().array()
202 
203 
204 if __name__ == "__main__":
205  # parse the input arguments
206  parser = ArgumentParser(description='Get a program and run it with input')
207  parser.add_argument('--version', action='version', version='%(prog)s 1.0')
208  parser.add_argument('--model', default="arctan", type=str, help="define true model")
209  parser.add_argument('--degree', default=5, type=int, help="polynomial degree of B-splines")
210  parser.add_argument('--minLevel', default=0, type=int, help="minimum level of regular grids")
211  parser.add_argument('--maxLevel', default=4, type=int, help="maximum level of regular grids")
212  parser.add_argument('--maxNumGridPoints', default=200, type=int, help="maximum number of grid points")
213  parser.add_argument('--growthFactor', default=2, type=int, help="Leja growth factor")
214  parser.add_argument('--levelManager', default="variance", type=str, help="define level manager")
215  parser.add_argument('--dist', default="beta", type=str, help="define marginal distribution")
216  args = parser.parse_args()
217 
218  model, params = buildModel(args.model, args.dist)
219 
220  # We have to wrap f in a pysgpp.MultiFunction object.
221  func = pysgpp.multiFunc(lambda x: model(x, params))
222  numDims = params.getStochasticDim()
223 
224  # compute reference values
225  n = 10000
226  x = params.getIndependentJointDistribution().rvs(n)
227  y = np.array([model(xi, params) for xi in x])
228 
229  results = {}
230  for gridType, levelManagerType, basisType in [
231  ("ClenshawCurtis", "variance", "poly"),
232  ("ClenshawCurtis", "averaging", "poly"),
233  ("UniformBoundary", "averaging", "bspline"),
234  ("UniformBoundary", "regular", "bspline"),
235  ("L2Leja", "variance", "poly"),
236  # ("Leja", "variance", "poly"),
237  ("ClenshawCurtis", "regular", "poly")]:
238 
239  refinement_wrapper = RefinementWrapper(gridType,
240  basisType,
241  args.degree,
242  args.growthFactor,
243  levelManagerType,
244  distType=args.dist,
245  samples=x,
246  params=params)
247  numGridPoints = np.array([])
248  l2errors = np.array([])
249  adaptive = levelManagerType != "regular"
250 
251  if adaptive:
252  n_sequence = np.unique(np.logspace(0, np.log10(args.maxNumGridPoints), dtype="int"))
253  else:
254  n_sequence = list(range(args.minLevel, args.maxLevel + 1))
255 
256  for i, n in enumerate(n_sequence):
257  # refine the grid
258  n_grid_points = refinement_wrapper.operation.numGridPoints()
259 
260  # compute L2 error
261  y_surrogate = refinement_wrapper.evaluate(n)
262 
263  if i == 0 or numGridPoints[-1] < n_grid_points:
264  print(( gridType, basisType, levelManagerType, n, n_grid_points ))
265 
266  l2error = np.sqrt(np.mean((y - y_surrogate) ** 2))
267  l2errors = np.append(l2errors, l2error)
268  numGridPoints = np.append(numGridPoints, refinement_wrapper.operation.numGridPoints())
269 
270  results[basisType, levelManagerType, gridType] = numGridPoints, l2errors
271 
272  print( "E(u) ~ %g" % np.mean(y) )
273  print( "V(u) ~ %g" % np.var(y) )
274 
275  fig = plt.figure()
276  ax = plt.subplot(111)
277 
278  for (basisType, levelManagerType, gridType), (numGridPoints, l2errors) in list(results.items()):
279  plt.loglog(numGridPoints, l2errors, label="%s %s %s" % (gridType, levelManagerType, basisType))
280 
281  plt.xlabel("number of grid points")
282  plt.ylabel("L2 error")
283  plt.tight_layout()
284  insert_legend(fig, loc="right", ncol=2)
285  plt.show()