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}