MandelBrot Set ด้วย Haskell

จากโจทย์ใน codenone ที่ Implement ด้วย Ruby ไปแล้ว
คราวนี้มาลอง implement ด้วย haskell บ้าง

เริ่มด้วย algorithm การคำนวณค่า
ถ้าเป็น ruby เขียนแบบนี้

def is_mandelbrot(x,y) 
  x0 = x
  y0 = y
  x2 = x*x
  y2 = y*y
 
  iter = 0
  maxiter = 30
 
  while ( x2 + y2 < 4  &&  iter < maxiter ) 
    y = 2*x*y + y0
    x = x2 - y2 + x0
    x2 = x*x    
    y2 = y*y
    iter += 1
  end
 
  if iter == maxiter || iter == 0
    0.0
  else
    iter.to_f * 100 / maxiter
  end
  
end

ใน imperative language จะมีพวก while loop ให้ใช้
แต่ใน haskell ไม่มี loop แบบนี้
ดังนั้นเราต้องแปลงให้เป็น tail-recursive แทน

max_iter = 30
is_mandel x y = is_mandel' x y max_iter x y
is_mandel' x y cnt x0 y0 
    | cnt == 0 = 0.0
    | exceed = (max_iter - cnt) / max_iter
    | otherwise = is_mandel' (x2-y2+x0) (2*x*y+y0) (cnt-1) x0 y0 
    where x2 = x * x
          y2 = y * y
          exceed = x2 * y2 > 4

ขั้นถัดไปก็คือคำนวณว่าจะ render ภาพอย่างไร
โดยในโจทย์ เขาจะให้เรากำหนด coordinate ของมุมบนซ้าย กับมุมล่างขวา
จากนั้นก็ให้เราคำนวณรูป โดยกำหนดว่า resolution ของการ plot ต้องมีอย่างน้อย 200 pixel,
ถ้าเขียนด้วย python จะเขียนดังนี้ (code นี้คุณสุกรีเป็นคนเขียน)

def mandelbrot(ul,lr,callback,step=(200,200),maxiter=255,limit=2):
    xs,ys = step
    x0,y0 = ul.real,ul.imag
    x1,y1 = lr.real,lr.imag
    xd,yd = (x1-x0)/xs,(y1-y0)/ys
    i = 0
    y = y0
    for i in range(ys):
        x = x0
        for j in range(xs):
            color = mandelbrot_point(x,y,maxiter,limit)
            callback(i,j,x,y,color)
        y += yd

จะเห็นว่ามี for loop ซ้อนกัน 2 ชั้น
haskell ก็ไม่มี for loop ให้ใช้เหมือนเดิม
ดังนั้น เราจะเปลียนแนวคิดไปเป็น List แทน
เริ่มโดยหาจุดที่เป็นไปได้บน แกน X ก่อน
โดยต้องการ list หน้าตาแบบนี้
[(-2,0),(-1.98,1),(-1.96,2),....]
ตัวแรกใน tuple คือ ค่าที่จะนำไปคำนวณสมการ
ส่วนตัวที่สองใน tuple ก็คือค่าที่ใช้ plot จุดบนแกน x

loopX startX endX = take w $ iterate (\(x,x') -> (x + stepX, x' + 1)) (startX,0)
    where stepX = (endX - startX) / width
          w = floor width

take คือการดึงค่าออกจาก list ตามจำนวนที่ระบุ เช่น
take 3 [1,2,3,4,5] ก็จะได้ [1,2,3]
ส่วน iterate มันอธิบายยาก ลองดูตัวอย่าง
iterate (\x -> x + 1) 1 ได้ผลลัพท์คือ [1,2,3,4,... จน infinity]
iterate (\x -> x + 2) 1 ได้ผลลัพท์คือ [1,3,5,7,...]
ดังนั้น code ข้างบน ก็อ่านได้ว่า
"ดึง element ออกมาจาก list ตามจำนวน pixel ที่เราจะ plot"
ซึ่งก็คือ 200 element

จากนั้นก็หาจุดที่เป็นไปได้บน แกน Y

loopY startY endY = take h $ iterate (\(y,y') -> (y - stepY, y' + 1)) (startY,0)
    where stepY = (startY - endY) / height
          h = floor height

นำ list ของแกน X และ แกน Y
มาทำ cartesial product โดยใช้ List Comprehension
โดย sx คือ x ที่มุมบนซ้าย
ex คือ x ที่มุมล่างขวา
sy คือ y ที่มุมบนซ้าย
ey คือ y ที่มุมล่างขวา

mandelset sx sy ex ey = [(x',y',is_mandel x y) | (x,x') <- (loopX sx ex), (y,y') <- (loopY sy ey)]

ถึงขั้นนี้ เราก็ได้ผลลัพท์ออกมาเป็น list ของ tuple แบบนี้
[(x coordinate,y coordinate, ค่าที่ได้จากการคำนวณว่าอยู่ใน set ของ mandelbrot )]
ซึ่งสามารถนำไป plot ได้แล้ว

ขั้นตอนที่ยากอีกขั้นของ Haskell ก็คือ จะเลือกใช้ Graphic library อันไหนดี
เพราะมันมีให้เลือกเยอะ,
แต่ที่เหมือนกันหมดทุกเจ้าก็คือ
ถ้าอยากรู้ว่าทำงานอย่างไร ก็ให้อ่าน source code ของ example เอา

ผมตกลงเลือกใช้ Graphics.UI.GLUT
(หลังจากเลือกไปตัวหนึ่งแล้ว แต่มัน run ช้ามาก,
กับอีกตัวหนึ่งที่ install ไม่ผ่าน)

ปัญหาของโปรแกรมในส่วนของ UI ก็คือ State
เนื่องจาก Haskell มันเป็น pure functional
การจัดการกับ state ก็เลยต้องใช้ Monad เข้ามาช่วย
ลองดูการ define state
โดยผมจะเก็บค่า มุมบนซ้ายกับมุมล่างขวา
ซึ่งค่านี้จะเปลี่ยนไปเรื่อยๆ เมื่อ user กด zoom

data State = State { sx, sy, ex, ey :: IORef Float }

makeState :: IO State
makeState = do 
  sx <- newIORef (-2.0)
  sy <- newIORef 2.0
  ex <- newIORef 2.0
  ey <- newIORef (-2.0)
  return $ State { sx = sx, sy = sy, ex = ex, ey = ey }

เวลาจะ get ค่า ก็ทำแบบนี้
sx <- get (sx state)
ส่วนเวลาจะ set ค่า ก็ต้องทำแบบนี้
writeIORef (sx state) val

code ที่เหลือเป็น code ที่เกี่ยวกับ open-gl แล้ว
มีหน้าตาประมาณนี้ (ไม่ได้แสดงส่วนที่เกี่ยวกับการ zoom)

display :: State -> DisplayCallback
display state = do
   clear [ ColorBuffer ]
   sx <- get (sx state)
   sy <- get (sy state)
   ex <- get (ex state)
   ey <- get (ey state)
   renderPrimitive Points $ do
            mapM_ plot $ mandelset sx sy ex ey
   swapBuffers
   
reshape :: ReshapeCallback
reshape size@(Size w h) = do
   viewport $= (Position 0 0, size)
   matrixMode $= Projection
   loadIdentity
   ortho2D 0 (fromIntegral w) 0 (fromIntegral h) 
   scale 1 (-1) (1::GLfloat)
   translate (Vector3 0 (-400) (0 :: GLfloat))


keyboardMouse :: State -> KeyboardMouseCallback
keyboardMouse state (MouseButton b) Down _ (Position x y) =  case b of
      LeftButton -> zoom state x y
      _ -> return ()
keyboardMouse _ _ _ _ _ = return ()

main :: IO ()
main = do
  (progName, _args) <- getArgsAndInitialize
  initialDisplayMode $= [ DoubleBuffered, RGBMode ]
  initialWindowSize $= Size (floor width) (floor height)
  createWindow progName
  myInit
  state <- makeState 
  displayCallback $= display state
  reshapeCallback $= Just reshape
  keyboardMouseCallback $= Just (keyboardMouse state)
  mainLoop