use parBuffer
[nofib.git] / parallel / mandel / Mandel.lhs
1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 %%%     This is for when it is used in the report      %%%
3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4
5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
6 %%%         When used as a standalone document         %%%
7 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
8 \begin{onlystandalone}
9 \begin{code}
10 {-# LANGUAGE BangPatterns #-}
11 module Mandel where
12 import Complex -- 1.3
13 import PortablePixmap
14 import Control.Parallel
15 import Control.Parallel.Strategies (using)
16 -- import qualified NewStrategies as NS
17 default ()
18 \end{code}
19 \end{onlystandalone}
20
21         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
22 The Mandelbrot set is the set of complex numbers for which the values
23 of the polynomial
24 \begin{equation}
25 \label{poly}
26 f_c(z) = z^2 + c
27 \end{equation}
28 are connected\cite{falconer}. A graphical representation of this set
29 can be rendered by plotting for different points in the complex plane
30 an intensity that is proportional to either the rate of divergence of
31 the above polynomial, or a constant if the polynomial is found to converge.
32
33 We present a \DPHaskell{} implementation of the mandelbrot set in a bottom up
34 manner that can be decomposed into the following stages :
35
36 \begin{enumerate}
37 \item   First we generate an infinite list of complex numbers caused by
38         applying the mandelbrot polynomial to a single complex number
39         in the complex plane.
40 \item   Next we walk over the list, determining the rate of divergence of
41         the list of complex numbers.
42 \item   We parallelize the above stages, by mapping the composition
43         of the preceeding functions over the complex plain.
44 \item   Lastly we describe how the complex plain is initialised, and the
45         results are graphically rendered.
46 \end{enumerate}
47
48         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
49 \section{Generating the list}
50
51 \begin{code}
52 mandel::(Num a) => a -> [a]
53 mandel c = infiniteMandel
54            where
55                 infiniteMandel = c : (map (\z -> z*z +c) infiniteMandel)
56 \end{code}
57
58 The definition of @mandel@ uses the numerically overloaded functions
59 @*@ and @+@; operationally it generates an infinite list of numbers caused by
60 itearting a number (i.e an object that has a type that belongs to class @Num@),
61 with the function @\z -> z*z +c@. For example the evaluation of ``@mandel 2@''
62 would result in the list of integers :
63
64 \begin{pseudocode}
65 [2, 6, 38, 1446, 2090918, 4371938082726, 19113842599189892819591078,
66 365338978906606237729724396156395693696687137202086, ^C{Interrupted!}
67 \end{pseudocode}
68
69 , whereas if the function were applied to the complex number\footnote{complex
70 numbers are defined in Haskells prelude as the algebraic data type
71 @data (RealFloat a) => Complex a = a :+ a@} (1.0 :+ 2.0), then the functions
72 behaviour would be equivalent of the mandelbrot polynomial(\ref{poly}) in
73 which \hbox{$c = 1.0 + i2.0$} :
74
75 \begin{pseudocode}
76 [(1.0 :+ 2.0), (-2.0 :+ 6.0), (-31.0 :+ -22.0), (478.0 :+ 1366.0),
77 (-1637471.0 :+ 1305898.0), ^C{Interrupted!}
78 \end{pseudocode}
79
80         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
81 \section{Walking the list}
82
83 ``@whenDiverge@'' determines the rate of divergence of a list of complex
84 numbers. If diveregence cannot be proved within a fixed number of
85 iteration of the mandelbrot polynomial, convergence is assumed.
86 This process is encoded in \DPHaskell{} by @take@ing a @limit@ of
87 complex numbers off the infinite stream generated by the mandel function.
88 This {\em finite} list is then traversed, each element of which is checked
89 for divergence.
90
91 \begin{code}
92 whenDiverge::  Int -> Double -> Complex Double -> Int
93 whenDiverge limit radius c
94   = walkIt (take limit (mandel c))
95   where
96      walkIt []     = 0                                   -- Converged
97      walkIt (x:xs) | diverge x radius  = 0               -- Diverged
98                    | otherwise           = 1 + walkIt xs -- Keep walking
99 \end{code}
100
101 \begin{equation}
102 \|x + iy\|_{2} = \sqrt{x^2 + y^2}
103 \end{equation}
104
105 We assume that divergence occurs if the norm ($\|x\|_{2}$) of a
106 complex number exceeds the radius of the region to be rendered.
107 The predicate @diverge@ encodes the divergence check, where @magnitude@
108 is the prelude function that calculates the euclidean norm
109 ($\|x\|_{2}$) of a complex number.
110
111 \begin{code}
112 -- VERY IMPORTANT FUNCTION: sits in inner loop
113
114 diverge::Complex Double -> Double -> Bool
115 diverge cmplx radius =  magnitude cmplx > radius
116 \end{code}
117
118         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
119 \section{Parallelising the mandelbrot set}
120
121 The mandelbrot set can be trivially parallelised by ``mapping''
122 the @whenDiverge@ function over a complex plain of values.
123
124
125 \begin{code}
126 parallelMandel:: [[Complex Double]] -> Int -> Double -> [Int]
127 parallelMandel mat limit radius
128    = concat $ 
129 -- NewStrategies version:
130 --       NS.parListBuffer 50 (NS.seqList id) $
131 --       map (map (whenDiverge limit radius)) mat
132
133 -- NewStrategies version:
134 --       NS.parListBufferRev 50 (NS.seqList id) $
135 --       map (map (whenDiverge limit radius)) mat
136
137 -- lazyParList version:
138 --                  lazyParList 50
139 --                              [ let l = map (whenDiverge limit radius) xs
140 --                                in seqList l `pseq` l
141 --                              | xs <- mat ]
142
143 --   lazyParList1 version:
144                  parBuffer 70
145                              [ let l = map (whenDiverge limit radius) xs
146                                in seqList l `pseq` l
147                              | xs <- mat ]
148
149 --   = lazyParListChunk 100 100 $ map (whenDiverge limit radius) mat
150 --   = lazyParMap 512 (whenDiverge limit radius) mat
151
152 parBuffer :: Int -> [a] -> [a]
153 parBuffer n xs = return xs (start n xs)
154   where
155     return (x:xs) (y:ys) = y `par` (x : return xs ys) 
156     return xs [] = xs
157
158     start !n [] = []
159     start 0 ys = ys
160     start !n (y:ys) = y `par` start (n-1) ys
161
162 -- parListN :: Int -> [a] -> [a]
163 -- parListN 0  xs     = xs 
164 -- parListN !n []     = []
165 -- parListN !n (x:xs) = x `par` parListN (n-1) xs
166
167 lazyParList :: Int -> [a] -> [a]
168 lazyParList !n xs = go xs (parListN n xs)
169    where 
170          go []     _ys    = []
171          go (x:xs) []     = x : xs
172          go (x:xs) (y:ys) = y `par` (x : go xs ys)
173
174 lazyParList1 :: Int -> [a] -> [a]
175 lazyParList1 !n xs = go xs (parListN1 n xs [])
176    where 
177          go []     _ys    = []
178          go (x:xs) []     = x : xs
179          go (x:xs) (y:ys) = y `par` (x : go xs ys)
180
181 -- parMap :: (a -> b) -> [a] -> [b]
182 -- parMap f [] = []
183 -- parMap f (x:xs) = let fx = f x; fxs = parMap f xs in fx `par` fxs `pseq` fx:fxs
184
185 parList :: [a] -> ()
186 parList [] = ()
187 parList (x:xs) = x `par` parList xs
188
189 -- parListN version 1: leads to fights as all capabilities try to
190 -- steal the early sparks, and the main thread gets blocked.
191 parListN :: Int -> [a] -> [a]
192 parListN 0  xs     = xs 
193 parListN !n []     = []
194 parListN !n (x:xs) = x `par` parListN (n-1) xs
195
196 -- like parListN, but starts the sparks in reverse order
197 parListN1 :: Int -> [a] -> [a] -> [a]
198 parListN1 0  xs     ys = parList ys `pseq` xs
199 parListN1 !n []     ys = parList ys `pseq` []
200 parListN1 !n (x:xs) ys = parListN1 (n-1) xs (x:ys)
201
202 seqList :: [a] -> ()
203 seqList [] = ()
204 seqList (x:xs) = x `pseq` seqList xs
205 -- 
206 -- parListChunk :: Int -> [a] -> ()
207 -- parListChunk n [] = ()
208 -- parListChunk n xs = let (ys,zs) = splitAt n xs in 
209 --                     seqList ys `par` parListChunk n zs
210
211 -- parListChunkWHNF :: Int -> [a] -> [a]
212 -- parListChunkWHNF n
213 --   = concat
214 --   . (`using` parList)
215 --   . map (`using` seqList)
216 --   . chunk n
217
218 -- chunk n [] = []
219 -- chunk n xs = as : chunk n bs where (as,bs) = splitAt n xs
220 -- 
221 -- lazyParList :: Int -> [a] -> [a]
222 -- lazyParList !n xs = go xs (parListN' n xs [])
223 --   where 
224 --         go []     _ys    = []
225 --         go (x:xs) []     = x : xs
226 --         go (x:xs) (y:ys) = y `par` (x : go xs ys)
227
228 -- lazyParListChunk :: Int -> Int -> [a] -> [a]
229 -- lazyParListChunk !n !size xs = go chunks seqchunks (parListN n seqchunks)
230 --   where 
231 --         chunks = chunkList size xs
232 --         seqchunks = map seqList chunks
233 -- 
234 --         go :: [[a]] -> [()] -> [()] -> [a]
235 --         go []     _ _ys         = []
236 --         go (x:xs) _ []          = concat (x:xs)
237 --         go (x:xs) (y:ys) (z:zs) = z `par` y `pseq` (x ++ go xs ys zs)
238 -- 
239 -- chunkList :: Int -> [a] -> [[a]]
240 -- chunkList !n [] = []
241 -- chunkList !n xs = chunk : chunkList n rest
242 --   where (chunk,rest) = splitAt n xs
243 \end{code}
244
245 \section{Initialisation of data and graphical rendering.}
246
247 We render the mandelbrot set by producing a PixMap file that
248 provides a portable representation of a 32bit RGB picture. The format of
249 this file is outside the scope of this paper, but we provide an interface
250 to this representation by defining an abstract data type @Pixmap@ that
251 has the functions @createPixmap@ (signature shown below), and @show@
252 defined (overloaded function from class @Text@).
253
254 \begin{pseudocode}
255 createPixmap :: Integer ->              -- The width of the Pixmap
256                 Integer ->              -- The height of the Pixmap
257                 Int ->                  -- The depth of of the RGB colour
258                 [(Int,Int,Int)] ->      -- A matrix of the RGB colours
259                 PixMap
260 \end{pseudocode}
261
262 @mandelset@ controls the evaluation and rendering of the mandelbrot set. The
263 function requires that a ``viewport'' on the complex plane is defined,
264 such that the region bounded by the viewport will be rendered in a ``Window''
265 represented by the pixmap.
266 \begin{code}
267 mandelset::Double ->                    -- Minimum X viewport
268            Double ->                    -- Minimum Y viewport
269            Double ->                    -- Maximum X viewport
270            Double ->                    -- maximum Y viewport
271            Integer ->                   -- Window width
272            Integer ->                   -- Window height
273            Int ->                       -- Window depth
274            PixMap                       -- result pixmap
275 mandelset x y x' y' screenX screenY lIMIT
276    = createPixmap screenX screenY lIMIT (map prettyRGB result)
277    where
278 \end{code}
279
280 @windowToViewport@ is a function that is defined within the scope of the
281 mandelset function (therefore arguments to mandelset will be in scope). The
282 function represents the relationship between the pixels in a window,
283 and points on a complex plain.
284
285 \begin{code}
286       windowToViewport s t
287            = ((x + (((coerce s) * (x' - x)) / (fromInteger screenX))) :+
288               (y + (((coerce t) * (y' - y)) / (fromInteger screenY))))
289
290       coerce::Integer -> Double
291       coerce  s   = encodeFloat (toInteger s) 0
292 \end{code}
293
294 The complex plain is initialised, and the mandelbrot set is calculated.
295 \begin{code}
296       result = parallelMandel
297                   [[windowToViewport s t | s<-[1..screenX]]
298                    | t <- [1..screenY]]
299                   lIMIT
300                   ((max (x'-x) (y'-y)) / 2.0)
301
302       prettyRGB::Int -> (Int,Int,Int)
303       prettyRGB s = let t = (lIMIT - s) in (s,t,t)
304 \end{code}