1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
|
Compiling to an executable file with the @-O2@ optimization level is recomended.
For example: @ghc -o 'mandelbrot_iter' -O2 MandelbrotIter.lhs ; ./mandelbrot_iter@
> {-# LANGUAGE PackageImports #-}
> import Control.Arrow
> import Control.Monad.Trans.Iter
> import "mtl" Control.Monad.Reader
> import "mtl" Control.Monad.List
> import "mtl" Control.Monad.Identity
> import Control.Monad.IO.Class
> import Data.Complex
> import Graphics.HGL (runGraphics, Window, withPen,
> line, RGB (RGB), RedrawMode (Unbuffered, DoubleBuffered), openWindowEx,
> drawInWindow, mkPen, Style (Solid))
Some fractals can be defined by infinite sequences of complex numbers. For example,
to render the <https://en.wikipedia.org/wiki/Mandelbrot_set Mandelbrot set>,
the following sequence is generated for each point @c@ in the complex plane:
@
z₀ = c
z₁ = z₀² + c
z₂ = z₁² + c
…
@
If, after some iterations, |z_i| ≥ 2, the point is not in the set. We
can compute if a point is not in the Mandelbrot set this way:
@
escaped :: Complex Double -> Int
escaped c = loop 0 0 where
loop z n = if (magnitude z) >= 2 then n
else loop (z*z + c) (n+1)
@
If @c@ is not in the Mandelbrot set, we get the number of iterations required to
prove that fact. But, if @c@ is in the mandelbrot set, 'escaped' will
run forever.
We can use the 'Iter' monad to delimit this effect. By applying
'delay' before the recursive call, we decompose the computation into
terminating steps.
> escaped :: Complex Double -> Iter Int
> escaped c = loop 0 0 where
> loop z n = if (magnitude z) >= 2 then return n
> else delay $ loop (z*z + c) (n+1)
>
If we draw each point on a canvas after it escapes, we can get a _negative_
image of the Mandelbrot set. Drawing pixels is a side-effect, so it
should happen inside the IO monad. Also, we want to have an
environment to store the size of the canvas, and the target window.
By using 'IterT', we can add all these behaviours to our non-terminating
computation.
> data Canvas = Canvas { width :: Int, height :: Int, window :: Window }
>
> type FractalM a = IterT (ReaderT Canvas IO) a
Any simple, non-terminating computation can be lifted into a richer environment.
> escaped' :: Complex Double -> IterT (ReaderT Canvas IO) Int
> escaped' = liftIter . escaped
Then, to draw a point, we can just retrieve the number of iterations until it
finishes, and draw it. The color will depend on the number of iterations.
> mandelbrotPoint :: (Int, Int) -> FractalM ()
> mandelbrotPoint p = do
> c <- scale p
> n <- escaped' c
> let color = if (even n) then RGB 0 0 255 -- Blue
> else RGB 0 0 127 -- Darker blue
> drawPoint color p
The pixels on the screen don't match the region in the complex plane where the
fractal is; we need to map them first. The region we are interested in is
Im z = [-1,1], Re z = [-2,1].
> scale :: (Int, Int) -> FractalM (Complex Double)
> scale (xi,yi) = do
> (w,h) <- asks $ (fromIntegral . width) &&& (fromIntegral . height)
> let (x,y) = (fromIntegral xi, fromIntegral yi)
> let im = (-y + h / 2 ) / (h/2)
> let re = ( x - w * 2 / 3 ) / (h/2)
> return $ re :+ im
Drawing a point is equivalent to drawing a line of length one.
> drawPoint :: RGB -> (Int,Int) -> FractalM ()
> drawPoint color p@(x,y) = do
> w <- asks window
> let point = line (x,y) (x+1, y+1)
> liftIO $ drawInWindow w $ mkPen Solid 1 color (flip withPen point)
We may want to draw more than one point. However, if we just sequence the computations
monadically, the first point that is not a member of the set will block the whole
process. We need advance all the points at the same pace, by interleaving the
computations.
> drawMandelbrot :: FractalM ()
> drawMandelbrot = do
> (w,h) <- asks $ width &&& height
> let ps = [mandelbrotPoint (x,y) | x <- [0 .. (w-1)], y <- [0 .. (h-1)]]
> interleave_ ps
To run this computation, we can just use @retract@, which will run indefinitely:
> runFractalM :: Canvas -> FractalM a -> IO a
> runFractalM canvas = flip runReaderT canvas . retract
Or, we can trade non-termination for getting an incomplete result,
by cutting off after a certain number of steps.
> runFractalM' :: Integer -> Canvas -> FractalM a -> IO (Maybe a)
> runFractalM' n canvas = flip runReaderT canvas . retract . cutoff n
Thanks to the 'IterT' transformer, we can separate timeout concerns from
computational concerns.
> main :: IO ()
> main = do
> let windowWidth = 800
> let windowHeight = 480
> runGraphics $ do
> w <- openWindowEx "Mandelbrot" Nothing (windowWidth, windowHeight) DoubleBuffered (Just 1)
> let canvas = Canvas windowWidth windowHeight w
> runFractalM' 100 canvas drawMandelbrot
> putStrLn $ "Fin"
|
