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