added a lot of comments

haskell/flame.lhs

@@ -1,63 +1,118 @@
+ # Flame fractal in Haskell
+
+First computing the flame fractal is difficult to do in parallel.
+You can take a look at [electricsheep.org](http://electricsheep.org)
+to view some examples.
+
+Here is literate Haskell.
+
+To start your program in Haskell you have to first declare all you imports.
+I personally find the number of import annoying.
+Not really a language issue but more a usage issue.
+For example, instead of the actual
+
+
 > module Main where
-> import Data.Hashable
 > import Data.HashMap as Dict -- cabal install HashMap
+> import Data.Hashable
 > import Data.Maybe   as Maybe
 > import Data.Word (Word8)
-> 
+> -- I need to write picture files
-> import Codec.Picture -- cabal install juicyPixels FTW
+> -- I also prefer to declare my own Pixel data type
+> import Codec.Picture hiding (Pixel) -- cabal install juicyPixels FTW
 > import Control.Monad
-> import Control.Monad.State
 > import System.Environment (getArgs)
-> 
+
+I would have preferred a more concise alternative such as:
+
+~~~
+import Maybe, Word, Picture, System
+import HashMap as Dict
+~~~
+
+Not a big deal thought.
+
+Now, instead of using common types like `(Int,Int)` for point,
+I prefer to use the power of Haskell types to help me discover errors.
+Therefore I create a type for each of the element I need.
+
+I will need a type for 2D points, Colors, extended colors (added color with the number)
+Furthermore for efficiency reason I'll use "low level" Haskell.
+I will replace `data Point = P Int Int` by an unboxed strict variant.
+As it is a GHC optimization it is far more verbose to declare.
+We tell GHC to unbox our type using the `{-# UNPACK #-}` comment before each field.
+Furthermore to make each field _strict_ we use a `!` before the type.
+
 > -- Data types
-> --  Global argument passed to most functions
+> --  Global state
-> data Global = Global { imgWidth  :: Int
+> data Global = Global {
->                      , imgHeight :: Int
+>		  filename  :: String
->                      , nbPoints  :: Int }
+>       , imgWidth  :: Int
-> 
+>       , imgHeight :: Int
+>       , nbPoints  :: Int }
+>
 > -- Real Points
 > data Point = P {-# UNPACK #-} !Float
 >                {-# UNPACK #-} !Float
-> 
+>
-> data Color = Color {-# UNPACK #-} !Word8
+
->                    {-# UNPACK #-} !Word8
+A Pixel, just a position not the color.
->                    {-# UNPACK #-} !Word8
+It is mostly like a Point but with integer.
-> 
+This way I won't mess up between the two representation (on screen) and
-> data ExtColor = ExtColor {-# UNPACK #-} !Int
+into the mathematical plane.
->                          {-# UNPACK #-} !Int
+
->                          {-# UNPACK #-} !Int
+> data Pixel = Pixel {-# UNPACK #-} !Int
->                          {-# UNPACK #-} !Int
+>                    {-# UNPACK #-} !Int
-> 
+>                    deriving (Eq,Ord)
-> data YPixel = YPixel {-# UNPACK #-} !Int
+
->                      {-# UNPACK #-} !Int
+
->                      deriving (Eq,Ord)
+I need to talk a bit about the color data structure I use.
-> 
+Instead of simple RGB each color coded on 8bit.
-> instance Hashable YPixel where hashWithSalt n (YPixel x y) = hashWithSalt n (x,y)
+I need something a bit more complex.
-> 
+
-> type YMap = Map YPixel ExtColor
+Each time a point will be lighten I'll give it a color.
-> 
+But some point can be lighten multiple times and in different ways.
-> addExtColor (ExtColor r g b n) (ExtColor r' g' b' n') =
+Then each time I light the same point I add the new light color to this point.
->   ExtColor (r+r') (g+g') (b+b') (n+n')
+But in order to retrieve the right brightness,
-> 
+I remember how many times the point was lighten up.
-> cmap f (Color r g b) = Color (f r) (f g) (f b)
+Thus this strange structure for color.
-> ecmap f (ExtColor r g b n) = ExtColor (f r) (f g) (f b) n
+Note that for drawing the image I'll use a more standard `PixelRGB8` which is
-> 
+simply three `Word8` values.
-> gammaCorrection :: Float -> Color -> Color
+
-> gammaCorrection gamma = cmap (round . (**(1/gamma)) . fromIntegral)
+> data Color = Color {-# UNPACK #-} !Int
-> 
+>                    {-# UNPACK #-} !Int
-> colorToPixelRGB8 :: Color -> PixelRGB8
+>                    {-# UNPACK #-} !Int
-> colorToPixelRGB8 (Color r g b) = PixelRGB8 r g b
+>                    {-# UNPACK #-} !Int
-> 
+
-> colorFromExt :: ExtColor -> Color
+My hash key will be pixel, I then need to make my Pixel data type
-> colorFromExt (ExtColor r g b n) = Color (fromIntegral $ div r n)
+an instance of Hashable.
->                                         (fromIntegral $ div g n)
+
->                                         (fromIntegral $ div b n)
+> instance Hashable Pixel where
-> -- Basic functions
+>   hashWithSalt n (Pixel x y) = hashWithSalt n (x,y)
-> neg x = 0-x
+
-> 
+Now I can use `Map` with `Pixel`s as key:
-> rgb :: Word8 -> Word8 -> Word8 -> Color
+
-> rgb r g b = Color r g b
+> type YMap = Map Pixel Color
-> -- Colors (theme is solarized)
+
+I want to be able to add two `Color` and
+to translate to `PixelRGB8` colors which are used to make the image.
+
+
+> addColor (Color r g b n) (Color r' g' b' n') =
+>   Color (r+r') (g+g') (b+b') (n+n')
+>
+> colorFromExt :: Color -> PixelRGB8
+> colorFromExt (Color r g b n) = PixelRGB8 (fromIntegral $ div r n)
+>                                             (fromIntegral $ div g n)
+>                                             (fromIntegral $ div b n)
+
+In most of my project I use the solarized theme.
+I generally use only a small part of these colors.
+But it is simply generally better to use a basic scheme than to use hard colors.
+
+> -- Colors from the theme solarized
+> rgb :: Int -> Int -> Int -> Color
+> rgb r g b = Color r g b 1
 > black   = rgb 0   0  0
 > base03  = rgb 0   43  54
 > base02  = rgb 7   54  66
@@ -75,12 +130,13 @@
 > blue    = rgb 38  139 210
 > cyan    = rgb 42  161 152
 > green   = rgb 133 153 0
-> 
+
-> extend :: Color -> ExtColor
+> -- very basic change of representation between point and pixel
-> extend (Color r g b) = ExtColor (fromIntegral r) (fromIntegral g) (fromIntegral b) 1
+> pixelFromPoint (P x y) = Pixel (round x) (round y)
-> 
+
-> pixelFromPoint (P x y) = YPixel (round x) (round y)
+Next I needed some pseudo random number generation.
-> 
+I don't need real good random number generation inside the Random for now.
+
 > -- PSEUDO RANDOM NUMBER GENERATION
 > -- !!!!!!!! DONT WORK ON 32 BITS Architecture !!!!!!!
 > nextint n = (a*n + c) `rem` m
@@ -89,28 +145,86 @@
 >       c = 1
 >       m = 2^32
 > -- generate a random sequence of length k starting with some seed
-> randlist seed n = take n $ iterate nextint seed
+> randlist seed = iterate nextint seed
 > -- END OF PSEUDO RANDOM NUMBER GENERATION
-> 
+
-> 
+ ## The Flame Set
-> {-
+
->  - Flame Set
+Let $F_i$ be a finite family of functions.
->  -
+And let consider the set $S$ which is the minimal non trivial set which is
->  - S = U_{i} F_i(S)
+the fixed point of the union of these function.
->  -
+
->  - F_i being transformations
+More precisely, let F_i be transformations.
->  - General form:
+Then consider the set S such that for all i F_i(S)=S
->  -   F = affine .  linearcomp [variation] . affine
+Or equivalently :
->  -   affine is a linear function (x,y) -> (ax+by+c,dx+ey+f)
+
->  -   variation is some kind of function with some contraction properties
+$$ S = U_{i} F_i(S) $$
->         ex: (x,y) -> (x,y), (sin x, sin y), etc...
+
->  -   linearcomp [f] is a linear composition of functions: (x,y) -> Sum vi*f(x,y)
+Consider the set is non trivial (understand non empty), then There is at least
-> -}
+one point in it.
-> 
+Finding the exact set is a difficult task.
+But finding an approximation can be done this way:
+
+ Let S_0 = {x} where x is a random point in the unit square.
+ Let S_1 = x_1 = F_i(x)   for a random i
+ Let S_2 = x_2 = F_j(x_1) for a random j
+ ...
+ Let S_n = x_n = F_k(x_{n-1}) for a random k
+ ...
+
+Each S_n will be closer to S.
+At each step you add another point to S_i.
+Also to remove bad initialization we generally don't consider the 20th firsts
+steps. And we return only {x_21,....,x_n}.
+
+In order to find only interesting elements we much choose our F_i family wisely.
+Here are standard form of the F_i:
+
+    F_i = affine . linear_combination [variations] . affine
+
+affine are function of the following form:
+
+    affine (x,y) = (ax+by+c , dx+ey+f)
+
+It correspond to a  composition of translation, rotation and scaling.
+
+linear_combination are function of the following form:
+
+    linear_combination [v_i] = sum p_iv_i
+
+affine have 6 parameters,
+linear_combination [x] have length [x] parameters.
+
+And we can use many different fonctions for the variations.
+
+Example of variations:
+
+- (x,y) → (x,y)
+- (x,y) → (sin x,sin y)
+- (x,y) → (x/r^2,y/r^2)
+- (x,y) → (x sin r^2 - y cos r^2, x cos r^2 + y sin r^2)
+- (x,y) → ((x-y)(x+y)/r,2xy/r)
+
+Wich are coded here:
+
+> -- Some variations
+> vs :: [Point -> Point]
+> vs = [ \ (P x y) -> P x y
+>      , \ (P x y) -> P (sin x) (sin y)
+>      , \ (P x y) -> let r2 = x*x+y*y in P (x/r2) (y/r2)
+>      , \ (P x y) -> let r2 = x*x+y*y in P (x*(sin r2) - y*(cos r2)) (x*(cos r2) + y * (sin r2))
+>      , \ (P x y) -> let r = sqrt (x^2+y^2) in P ((x - y)*(x + y)/r) (2*x*y/r)
+>      ]
+
+To define affine function a standard usage is to use matrices.
+
 > data Matrice = M Float Float Float Float Float Float
 > aff :: Matrice -> Point -> Point
 > aff (M a b c d e f) (P x y) = P (a*x + b*y + c) (d*x + e*y +f)
-> 
+
+If you use the identity variation,
+the following functions generate the sierpinsky set.
+
 > -- Some affine functions to generate the sierpinsky set
 > -- Equivalent to
 > -- sierp = [ \(x,y)->(x/2,y/2)
@@ -127,56 +241,68 @@
 >           0.5  0.0  0.0
 >           0.0  0.5  0.5
 >         ]
+
+Here are the functions for the fern functions.
+
 > fern :: [ Point -> Point ]
 > fern = [ aff $ M
 >           0.0  0.0  0.0
 >           0.0  0.16 0.0
 >         , aff $ M
 >           0.85        0.04  0.0
->           (neg 0.04)  0.85  1.6
+>           (- 0.04)  0.85  1.6
 >         , aff $ M
->           0.2  (neg 0.26)  0.0
+>           0.2  (- 0.26)  0.0
 >           0.23  0.22       1.6
 >         , aff $ M
->           (neg 0.15)  0.28  0.0
+>           (- 0.15)  0.28  0.0
 >           0.26        0.24  0.44
 >         ]
-> 
+
-> -- Some variations
+Also in order to zoom on all points we generally add a final transformation
-> vs :: [Point -> Point]
+which is applied to all points. It helps zoom on the fractal for example.
-> vs = [ \ (P x y) -> P x y
+
->      , \ (P x y) -> P (sin x) (sin y)
->      , \ (P x y) -> let r2 = x*x+y*y in P (x/r2) (y/r2)
->      , \ (P x y) -> let r2 = x*x+y*y in P (x*(sin r2) - y*(cos r2)) (x*(cos r2) + y * (sin r2))
->      , \ (P x y) -> let r = sqrt (x^2+y^2) in P ((x - y)*(x + y)/r) (2*x*y/r)
->      ]
-> 
-> -- Some final functions
-> fs :: [((Int,ExtColor),Point -> Point)]
-> fs = [ ((  1,extend    red),(vs !! 0) . (fern !! 0))
->      , (( 86,extend  green),(vs !! 0) . (fern !! 1))
->      , (( 95,extend   blue),(vs !! 0) . (fern !! 2))
->      , ((100,extend yellow),(vs !! 0) . (fern !! 3))
->      ]
-> 
 > -- Transformation functions
 > -- translate
 > trans :: (Float,Float) -> Point -> Point
 > trans (tx,ty) = aff $ M 1 0 tx 0 1 ty
 > -- rotate
 > rot :: Float -> Point -> Point
-> rot phi = aff $ M (cos phi) (sin phi) 0.0 (neg (sin phi)) (cos phi) 0.0
+> rot phi = aff $ M (cos phi) (sin phi) 0.0 (- (sin phi)) (cos phi) 0.0
 > -- zoom
 > zoom :: Float -> Point -> Point
 > zoom z = aff $ M z 0 0 0 z 0
-> 
+
+As the final function goal is to help the final rendering position,
+it seems natural to add the size of the view as parameter.
+
 > -- The final transformation to transform the final result (zoom,rotate,translate)
 > final :: Int -> Point -> Point
-> final width = trans (w/2,w/2) . zoom (w/10)  . rot (neg pi)
+> final width = trans (w/2,w/2) . zoom (w/10)  . rot (- pi)
 >               where w = fromIntegral width
-> 
+
-> sierpset :: Int -> Point -> [Int] -> YMap -> YMap
+And now the F_i functions.
-> sierpset w startpoint rands tmpres =
+As we can see, it is not only a list of functions.
+But we add informations to each function:
+- a probability to be used,
+- a color.
+
+> -- F_i
+> fs :: [((Int, Color), Point -> Point)]
+> fs = [ ((  1,   red), (vs !! 0) . (fern !! 0))
+>      , (( 86, green), (vs !! 0) . (fern !! 1))
+>      , (( 95,  blue), (vs !! 0) . (fern !! 2))
+>      , ((100,yellow), (vs !! 0) . (fern !! 3))
+>      ]
+
+Up until now it was only a verbose babbling.
+Here is the heart of our program.
+Where the interesting stuff is going on.
+For now, this is a rather naive implementation.
+Naive in the sense that it doesn't use Monad as helper.
+
+> flameset :: Int -> Point -> [Int] -> YMap -> YMap
+> flameset w startpoint rands tmpres =
 >   if rands == []
 >     then tmpres
 >     else
@@ -193,39 +319,48 @@
 >         -- Search the old color
 >         oldvalue = Dict.lookup savepoint tmpres
 >         -- Set the new color.
->         newvalue = addExtColor col (Maybe.fromMaybe (extend black) oldvalue)
+>         newvalue = addColor col (Maybe.fromMaybe black oldvalue)
 >         -- update the dict
 >         newtmpres = Dict.insert savepoint newvalue tmpres
 >       in
->         sierpset w newpoint (tail rands) newtmpres
+>         flameset w newpoint (tail rands) newtmpres
-> 
+
-> sierpinsky :: Int -> Int -> YMap
+The flame function is just a call to the flameset function with initial values.
-> sierpinsky w n = sierpset w (P 0.13 0.47) (randlist 0 n) Dict.empty
+Clearly there is something to be done here.
-> 
+
+> flame :: Int -> Int -> YMap
+> flame w n = flameset w (P 0.13 0.47) (take n $ randlist 0) Dict.empty
+
+A function to read the command line arguments.
+
 > initGlobalParams args =
->   Global { imgWidth  = read (args !! 0)
+>   Global { filename  = args !! 0
->          , imgHeight = read (args !! 1)
+>		   , imgWidth  = read (args !! 1)
->          , nbPoints  = read (args !! 2) }
+>          , imgHeight = read (args !! 2)
-> 
+>          , nbPoints  = read (args !! 3) }
+
+The functions needed to transform the dictionary as a picture file.
+
 > imageFromDict :: YMap -> Int -> Int -> Image PixelRGB8
 > imageFromDict dict width height = generateImage colorOfPoint width height
 >   where
 >     colorOfPoint :: Int -> Int -> PixelRGB8
->     colorOfPoint x y = colorToPixelRGB8 $ colorFromExt $
+>     colorOfPoint x y = colorFromExt $
->                        fromMaybe (extend base03)
+>							fromMaybe base03 -- background color
->                                  (Dict.lookup (YPixel x y) dict)
+>                        	          (Dict.lookup (Pixel x y) dict)
-> 
+>
-> writeImage :: Int -> Int -> Int -> YMap -> IO ()
+> writeImage :: String -> Int -> Int -> Int -> YMap -> IO ()
-> writeImage w h n dict = writePng "flame.png" $ imageFromDict dict w h
+> writeImage filename w h n dict = writePng filename $ imageFromDict dict w h
-> 
+>
 > main :: IO ()
 > main = do
 >   args <- getArgs
->   if (length args<3)
+>   if (length args<4)
->     then print $ "Usage flame w h n"
+>     then print $ "Usage flame ficname w h n"
 >     else do
 >       env  <- return (initGlobalParams args)
+>       fic <- return (filename env)
 >       w <- return (imgWidth env)
 >       h <- return (imgHeight env)
 >       n <- return (nbPoints env)
->       writeImage w h n (sierpinsky w n)
+>       writeImage fic w h n (flame w n)