Note
Go to the end to download the full example code.
Iterated Functions SystemΒΆ
This examples show how to generate fractal sets using iterated functions systems. See https://en.wikipedia.org/wiki/Iterated_function_system for an introduction.
import openturns as ot
import openturns.viewer as viewer
from matplotlib import pyplot as plt
# sphinx_gallery_thumbnail_number = 4
import math as m
Tree traversal algorithm (the chaos game)
def drawIFS(f_i, skip=100, iterations=1000, batch_size=1, name="IFS", color="blue"):
# Any set of initial points should work in theory
initialPoints = ot.Normal(2).getSample(batch_size)
# Compute the contraction factor of each function
all_r = [m.sqrt(abs(f[1].computeDeterminant())) for f in f_i]
# Find the box counting dimension, ie the value s such that r_1^s+...+r_n^s-1=0
equation = "-1.0+" + "+".join([str(r) + "^s" for r in all_r])
dim = len(f_i)
fs = ot.SymbolicFunction("s", equation)
# tweak search bounds
xMin, xMax = 0.0, -m.log(dim) / m.log(max(all_r))
fMax = fs([xMax])[0]
eps = ot.SpecFunc.Precision**0.5
if abs(fMax) < eps:
xMax += eps
s = ot.Brent().solve(fs, 0.0, xMin, xMax)
# Add a small perturbation to sample even the degenerated transforms
probabilities = [r**s + 1e-2 for r in all_r]
# Build the sampling distribution
support = [[i] for i in range(dim)]
choice = ot.UserDefined(support, probabilities)
currentPoints = initialPoints
points = ot.Sample(0, 2)
# Convert the f_i into LinearEvaluation to benefit from the evaluation over
# a Sample
phi_i = [ot.LinearEvaluation([0.0] * 2, f[0], f[1]) for f in f_i]
# Burning phase
for i in range(skip):
index = int(round(choice.getRealization()[0]))
currentPoints = phi_i[index](currentPoints)
# Iteration phase
for i in range(iterations):
index = int(round(choice.getRealization()[0]))
currentPoints = phi_i[index](currentPoints)
points.add(currentPoints)
# Draw the IFS
graph = ot.Graph()
graph.setTitle(name)
graph.setXTitle("x")
graph.setYTitle("y")
graph.setGrid(True)
cloud = ot.Cloud(points)
cloud.setColor(color)
cloud.setPointStyle("dot")
graph.add(cloud)
return graph, s
Definition of some IFS
Spiral
rho1 = 0.9
theta1 = 137.5 * m.pi / 180.0
f1 = [
[0.0] * 2,
ot.SquareMatrix(
2,
[
rho1 * m.cos(theta1),
-rho1 * m.sin(theta1),
rho1 * m.sin(theta1),
rho1 * m.cos(theta1),
],
),
]
rho2 = 0.15
f2 = [[1.0, 0.0], rho2 * ot.IdentityMatrix(2)]
f_i = [f1, f2]
graph, s = drawIFS(
f_i, skip=100, iterations=100000, batch_size=1, name="Spiral", color="blue"
)
print("Box counting dimension=%.3f" % s)
view = viewer.View(graph)
Box counting dimension=1.146
Fern
f1 = [[0.0] * 2, ot.SquareMatrix(2, [0.0, 0.0, 0.0, 0.16])]
f2 = [[0.0, 1.6], ot.SquareMatrix(2, [0.85, 0.04, -0.04, 0.85])]
f3 = [[0.0, 1.6], ot.SquareMatrix(2, [0.2, -0.26, 0.23, 0.22])]
f4 = [[0.0, 0.44], ot.SquareMatrix(2, [-0.15, 0.28, 0.26, 0.24])]
f_i = [f1, f2, f3, f4]
graph, s = drawIFS(
f_i, skip=100, iterations=100000, batch_size=1, name="Fern", color="green"
)
print("Box counting dimension=%.3f" % s)
# sphinx_gallery_thumbnail_number = 2
view = viewer.View(graph)
Box counting dimension=1.834
Dragon
f1 = [[0.0, 0.0], ot.SquareMatrix(2, [0.5, -0.5, 0.5, 0.5])]
f2 = [[1.0, 0.0], ot.SquareMatrix(2, [-0.5, -0.5, 0.5, -0.5])]
f_i = [f1, f2]
graph, s = drawIFS(
f_i, skip=100, iterations=100000, batch_size=1, name="Dragon", color="red"
)
print("Box counting dimension=%.3f" % s)
view = viewer.View(graph)
Box counting dimension=2.000
Sierpinski triangle
f1 = [[0.0, 0.0], ot.SquareMatrix(2, [0.5, 0.0, 0.0, 0.5])]
f2 = [[0.5, 0.0], ot.SquareMatrix(2, [0.5, 0.0, 0.0, 0.5])]
f3 = [[0.25, m.sqrt(3.0) / 4.0], ot.SquareMatrix(2, [0.5, 0.0, 0.0, 0.5])]
f_i = [f1, f2, f3]
graph, s = drawIFS(
f_i,
skip=100,
iterations=100000,
batch_size=1,
name="Sierpinski's triangle",
color="magenta",
)
print("Box counting dimension=%.3f" % s)
view = viewer.View(graph)
plt.show()
Box counting dimension=1.585