マンデルブロ集合

ｐｂｍファイルに書き出してマンデルブロ集合を描きます。（画像サイズは１０００ｘ１０００に固定しています。）

	--mandel.hs
	import Data.Array.Repa
	import Data.Complex
	import Control.DeepSeq
	import Control.Parallel.Strategies
	import System.Environment
	getXY :: IO ((Double,Double),Double)
	getXY = getArgs >>= (\[x,y,xx,_] -> return ((read x,read y),read xx))
	getRep :: IO Int
	getRep = getArgs >>= (\ [_,_,_,rep] -> return $ read rep)
	mkRowInD :: Double -> Double -> Int -> [Double]
	mkRowInD a gap n = take n [x | x <- [a,(a+gap)..]]
	mkRowInRepa :: [Double] -> Array U DIM1 Double
	mkRowInRepa l = fromListUnboxed (Z :. (length l)::DIM1) l
	mkImgParts y gap n = mkRowInRepa ds
	where ds = mkRowInD y gap n
	initC :: (RealFloat a) => Complex a
	initC = 0:+0
	isMandel :: Double -> Double -> Int -> Char
	isMandel a b n = isMandel' initC z n where z = a :+ b
	isMandel' :: (Eq a, Num a, RealFloat a1) => Complex a1 -> Complex a1 -> a -> Char
	isMandel' z c n
	| n==0 = '1'
	| realPart newC >= 2 = '0'
	| imagPart newC >= 2 = '0'
	| otherwise = isMandel' newC c (n-1)
	where newC = z*z+c
	pMap = Prelude.map
	isMandels a0 d n r =
	take n $ pMap (\x -> let d0 = a0 ! (Z:.x)
	in isMandel d0 d r) [0..]
	lineAll a y gap n r =
	pMap (\x -> let d = (a' x) ! (Z:.x)in isMandels a d n r)
	[0..(n-1)]
	`using` parListChunk 100 rdeepseq
	where a' x = mkImgParts y gap n
	pix = 1000 :: Int
	main = do
	((x,y),xx) <- getXY
	rep <- getRep
	let gap = abs (x - xx) / fromIntegral pix
	let r = mkRowInD x gap pix
	let a0 = mkRowInRepa r
	putStrLn "P1 "
	print pix
	print pix
	print $ unwords $ lineAll a0 y (-gap) pix rep

view raw mandel.hs hosted with ❤ by GitHub

---

コンパイルします。

$ ghc -O2 -threaded -rtsopts --make -XFlexibleContexts -eventlog mandel.hs
$ ./mande -2 2 2 3000 +RTS -N2 -l | display

displayコマンドにパイプします。
動作を確認。
$ threadscope mamdel.eventlog

殆どの時間２Coreで動いています。

---

Data.Array.Repa.IO.BMP

モジュールには

readImageFromBMP :: FilePath -> IO (Either Error (Array U DIM2 (Word8, Word8, Word8)))

writeImageToBMP :: FilePath -> Array U DIM2 (Word8, Word8, Word8) -> IO ()

という関数がありBMPファイルの読み書きができる。


画像サイズもコマンドラインから与えます。

	import Data.Array.Repa.IO.BMP
	import Data.Array.Repa as R
	import Data.Complex
	import GHC.Word
	import System.Environment
	bk,wt :: (GHC.Word.Word8, GHC.Word.Word8, GHC.Word.Word8)
	bk = (0,0,0)
	wt = (255,255,255)
	getXY :: IO ((Double,Double),Double)
	getXY = getArgs >>= (\[x,y,xx,_,_] -> return ((read x,read y),read xx))
	getRep :: IO Int
	getRep = getArgs >>= (\ [_,_,_,_,rep] -> return $ read rep)
	getPix = getArgs >>= (\ [_,_,_,p,_] -> return $ read p)
	mkRowInD :: Double -> Double -> Int -> [Double]
	mkRowInD a gap n = Prelude.take n [x | x <- [a,(a+gap)..]]
	mkRowInRepa :: [Double] -> Array U DIM1 Double
	mkRowInRepa l = fromListUnboxed (Z :. (Prelude.length l)::DIM1) l
	mkRealPart x gap n = mkRowInRepa ds
	where ds = mkRowInD x gap n
	mkImgPart y gap n = mkRowInRepa ds
	where ds = mkRowInD y (-gap) n
	initC :: (RealFloat a) => Complex a
	initC = 0:+0
	isMandel a b n = isMandel' initC z n where z = a :+ b
	isMandel' z c n
	| n == 0 = True
	| realPart newC >= 2 = False
	| imagPart newC >= 2 = False
	| otherwise = isMandel' newC c (n-1)
	where newC = z*z+c
	lineAll x y r = R.map (\(a,b) ->
	if isMandel (x ! (Z :. b))
	(y ! (Z :. a))
	r
	then wt
	else bk)
	cros pix = [(x,y) | x <- [0..(pix-1)],y <- [0..(pix-1)]]
	main = do
	((x,y),xx) <- getXY
	pix <- getPix
	rep <- getRep
	let gap = abs (x - xx) / fromIntegral pix
	let ra = mkRealPart x gap pix
	let ia = mkImgPart y gap pix
	let vr = fromListUnboxed (Z :. (pix*pix)) $ cros pix :: Array U DIM1 (Int,Int)
	let bs = lineAll ra ia rep vr
	computeP $ reshape (Z :. pix :. pix) bs
	>>=
	(writeImageToBMP "a.bmp")

view raw mandel2.hs hosted with ❤ by GitHub

$ ./mandel2 -2 2 2 1000 3000 +RTS -N2 -l

1次元の配列に対し　map し、　reshape すればよいので、コード自体はわかりやすい。
index のTable (Int のタプルの配列)を作るのに時間がかかっているようなので、
配列を連結するやり方も試みた。

	import Data.Array.Repa.IO.BMP
	import Data.Array.Repa as R
	import Data.Complex
	import Data.Vector.Unboxed
	import Data.List(foldl1')
	import GHC.Word
	import System.Environment
	bk,wt :: (GHC.Word.Word8, GHC.Word.Word8, GHC.Word.Word8)
	bk = (0,0,0)
	wt = (255,255,255)
	getXY :: IO ((Double,Double),Double)
	getXY = getArgs >>= (\[x,y,xx,_,_] -> return ((read x,read y),read xx))
	getRep :: IO Int
	getRep = getArgs >>= (\ [_,_,_,_,rep] -> return $ read rep)
	getPix = getArgs >>= (\ [_,_,_,p,_] -> return $ read p)
	mkRowInD :: Double -> Double -> Int -> [Double]
	mkRowInD a gap n = Prelude.take n [x | x <- [a,(a+gap)..]]
	mkRowInRepa :: [Double] -> Array U DIM1 Double
	mkRowInRepa l = fromListUnboxed (Z :. (Prelude.length l)::DIM1) l
	mkRealPart x gap n = mkRowInRepa ds
	where ds = mkRowInD x gap n
	mkImgPart y gap n = mkRowInRepa ds
	where ds = mkRowInD y (-gap) n
	initC :: (RealFloat a) => Complex a
	initC = 0:+0
	isMandel a b n = isMandel' initC z n where z = a :+ b
	isMandel' z c n
	| n == 0 = True
	| realPart newC >= 2 = False
	| imagPart newC >= 2 = False
	| otherwise = isMandel' newC c (n-1)
	where newC = z*z+c
	oneLine is xd yd y pix r = R.map (\x -> isMandel (xd R.! (Z :. x)) (yd R.! (Z :. y)) r) is
	lineAll xs ys pix r =
	Prelude.map (\x -> oneLine is xs ys x pix r) [0.. (pix-1)]
	where is = fromUnboxed (Z :. (pix::Int)) (enumFromN (0 :: Int) pix)
	main = do
	((x,y),xx) <- getXY
	pix <- getPix
	rep <- getRep
	let gap = abs (x - xx) / fromIntegral pix
	let ra = mkRealPart x gap pix
	let ia = mkImgPart y gap pix
	let o = lineAll ra ia pix rep
	let pp = Data.List.foldl1' (\x y -> x R.++ y) o
	let ppp = R.map (\x -> if x then wt else bk) pp
	computeP $ reshape (Z :. pix :. pix) ppp
	>>=
	(writeImageToBMP "a.bmp")

view raw mandel3.hs hosted with ❤ by GitHub

---

また　 mandel2.hsにおいてIndexをIntではなくWord16に置き換えてみたが効果は微妙だった。

---

世の中にはGPUを使ってマンデルブロー集合を極めて高速に描画するプログラムもあるらしい。