1""" 2Generate low discrepancy sequences 3================================== 4""" 5# %% 6# In this examples we are going to expose the available low discrepancy sequences in order to approximate some integrals. 7# 8# The following low-discrepancy sequences are available: 9# 10# - Sobol 11# - Faure 12# - Halton 13# - reverse Halton 14# - Haselgrove 15# 16# To illustrate these sequences we generate their first 1024 points and compare with the sequence obtained from the pseudo random generator (Merse Twister) as the latter has a higher discrepancy. 17 18# %% 19from __future__ import print_function 20import openturns as ot 21import math as m 22import openturns.viewer as viewer 23from matplotlib import pylab as plt 24ot.Log.Show(ot.Log.NONE) 25 26# %% 27# 1. Sobol sequence 28dimension = 2 29size = 1024 30sequence = ot.SobolSequence(dimension) 31sample = sequence.generate(size) 32graph = ot.Graph("Sobol", "", "", True, "") 33cloud = ot.Cloud(sample) 34graph.add(cloud) 35view = viewer.View(graph) 36 37# %% 38# 2. Halton sequence 39dimension = 2 40sequence = ot.HaltonSequence(dimension) 41sample = sequence.generate(size) 42graph = ot.Graph("Halton", "", "", True, "") 43cloud = ot.Cloud(sample) 44graph.add(cloud) 45view = viewer.View(graph) 46 47# %% 48# 3. Halton sequence in high dimension: bad filling in upper dimensions 49dimension = 20 50sequence = ot.HaltonSequence(dimension) 51sample = sequence.generate(size).getMarginal([dimension-2, dimension-1]) 52graph = ot.Graph("Halton (" + str(dimension - 2) + "," + str(dimension-1) + ")", 53 "dim " + str(dimension-2), "dim " + str(dimension-1), True, "") 54cloud = ot.Cloud(sample) 55graph.add(cloud) 56view = viewer.View(graph) 57 58# %% 59# 4. Scrambled Halton sequence in high dimension 60dimension = 20 61sequence = ot.HaltonSequence(dimension) 62sequence.setScrambling("RANDOM") 63sample = sequence.generate(size).getMarginal([dimension-2, dimension-1]) 64graph = ot.Graph("Halton (" + str(dimension - 2) + "," + str(dimension-1) + ")", 65 "dim " + str(dimension-2), "dim " + str(dimension-1), True, "") 66cloud = ot.Cloud(sample) 67graph.add(cloud) 68view = viewer.View(graph) 69 70# %% 71# 5. Reverse Halton sequence 72dimension = 2 73sequence = ot.ReverseHaltonSequence(dimension) 74sample = sequence.generate(size) 75print('discrepancy=', 76 ot.LowDiscrepancySequenceImplementation.ComputeStarDiscrepancy(sample)) 77graph = ot.Graph("Reverse Halton", "", "", True, "") 78cloud = ot.Cloud(sample) 79graph.add(cloud) 80view = viewer.View(graph) 81 82# %% 83# 6. Haselgrove sequence 84dimension = 2 85sequence = ot.HaselgroveSequence(dimension) 86sample = sequence.generate(size) 87graph = ot.Graph("Haselgrove", "", "", True, "") 88cloud = ot.Cloud(sample) 89graph.add(cloud) 90view = viewer.View(graph) 91 92# %% 93# Compare with uniform random sequence 94distribution = ot.ComposedDistribution([ot.Uniform(0.0, 1.0)]*2) 95sample = distribution.getSample(size) 96print('discrepancy=', 97 ot.LowDiscrepancySequenceImplementation.ComputeStarDiscrepancy(sample)) 98graph = ot.Graph("Mersenne Twister", "", "", True, "") 99cloud = ot.Cloud(sample) 100graph.add(cloud) 101view = viewer.View(graph) 102plt.show() 103