Monday, March 6, 2017

butterfly curve

Temple Fay discovered this delightfully suggestive curve in 1989. It can be defined either parametrically or as a polar equation; I use the former method here.
{-# LANGUAGE NoMonomorphismRestriction #-}

import Data.Bits ((.|.), shiftL)
import Control.Arrow ((***), (&&&))
import Graphics.UI.SDL as SDL
import Graphics.UI.SDL.Primitives (pixel)

(xres, yres, zz) = (1050, 1050, round *** round)

main = withInit [InitVideo] $ do
  w <- setVideoMode xres yres 32 [NoFrame]
  enableEvent SDLMouseMotion False
  setCaption "Butterfly Curve" "Butterfly Curve"
  fillRect w (Just $ Rect 0 0 xres yres) $ Pixel 0
  plot w $ center $ scale curve
  loop w []
 where
  scale  = map ((150 *) *** (* 150))
  center = map (((xres / 2) +) *** (+ ((yres + 190) / 2)))

curve = map (f &&& (negate . g)) ts
 where
  ts  = [-999, -998.99.. 999]
  f t = sin t * (e ** cos t - 2 * cos (4*t) - h (t / 12))
  g t = cos t * (e ** cos t - 2 * cos (4*t) - h (t / 12))
  h   = foldr1 (.) $ replicate 5 sin
  e   = exp 1

plot w = mapM_ (f . zz)
 where
  f (x,y) = pixel w (fromIntegral x) (fromIntegral y) $ Pixel rgb
   where
    rgb   = rgb'  .|. shiftL (255 `div` (max x y `div` min x y)) 8
    rgb'  = rgb'' .|. shiftL (255 `div` (xres `div` x)) 16
    rgb'' = 0xFF  .|. shiftL (255 `div` (yres `div` y)) 24

loop w ps = do
  delay 128
  event <- pollEvent
  case event of
   KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
   _                              -> loop w ps

Friday, February 24, 2017

bézier curves

De Casteljau's algorithm is a fast, numerically stable way to rasterize Bézier curves (I was disappointed to see Wikipedia already gives succinct Haskell for it!). Clicking appends nodes to the curve's defining polygon; you can also drag any node to alter the curve.
{-# LANGUAGE NoMonomorphismRestriction #-}

import Control.Monad (when)
import Control.Arrow ((***))
import Data.List (findIndex)
import Graphics.UI.SDL as SDL hiding (init)
import Graphics.UI.SDL.Primitives (circle, line)

dim = 400

main = withInit [InitVideo] $ do
  w <- setVideoMode dim dim 32 []
  enableEvent SDLMouseMotion False
  setCaption "Bézier Curves" "Bézier Curves"
  loop w []

plot w ps = do
  fillRect w (Just $ Rect 0 0 dim dim) $ Pixel 0xFF222255
  mapM_ (f 2 0xFFFFFFFF . zz) [head ps, last ps]
  when b $ mapM_ (f 3 0x888888FF . zz) controls
  when b $ mapM_ (f 2 0xBBBBBBFF . zz) controls
 where
  f r c (x,y)   = circle w x y r $ Pixel c
  (b, controls) = (length ps > 2, tail $ init ps)

limn w [_] = SDL.flip w
limn w ((a,b):(x,y):ps) = do
  line w a b x y $ Pixel 0xFFFFFFFF
  limn w $ (x,y) : ps

loop w ps = do
  delay 128
  event <- pollEvent
  case event of
   KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
   MouseButtonDown x y _          -> click x y
   _                              -> loop w ps
 where
  click x y = let p = rr (x,y) in
   case findIndex ((10 >) . dist p) ps of
    Just i  -> drag w i ps
    Nothing -> do
      let ps' = p : ps
      plot w ps' >> SDL.flip w
      when (length ps' > 2) $ render w ps'
      loop w ps'

drag w i ps = do
  delay 16
  (x,y,_) <- getMouseState
  event   <- pollEvent
  let ps'  = swap i ps $ rr (x,y)
  plot w ps'
  when (length ps' > 2) $ render w ps'
  case event of
   MouseButtonUp x y _ -> loop w ps'
   _                   -> drag w i ps'

render w ps = limn w $ map zz curve
 where
  curve = map (casteljau ps) [0, 0.001.. 1]

casteljau [p] t = p
casteljau ps  t = casteljau ps' t
 where
  ps'             = zipWith (g t) ps $ tail ps
  g t (a,b) (c,d) = (f t a c, f t b d)
  f t a b         = (1 - t) * a + t * b

swap 0 ps p = p : tail ps
swap i ps p = take i ps ++ p : drop (i+1) ps

dist (a,b) (c,d) = sqrt $ (a-c)^2 + (b-d)^2

rr = fromIntegral *** fromIntegral
zz = round *** round

Wednesday, February 22, 2017

convex hulls

This code demonstrates the Graham scan, an O(n log n) method for finding the convex hull (smallest enclosing polygon) of a planar point set. It's a great example of exploiting order: it works by sorting the set by angle about a known extreme point, allowing the hull points to be found in linear time.

A variation on this algorithm, noted by A.M Andrew, sorts the set lexicographically and finds the upper and lower hull chains separately. This 'monotone chain' technique is often preferred, since it's easier to do robustly.

I've also written a couple JavaScript examples, the Graham scan and a slower (but constant-space and adaptable to higher dimensions) method called the Jarvis march.

{-# LANGUAGE NoMonomorphismRestriction #-}

import Control.Arrow ((***))
import Data.List (maximumBy, delete, sort, sortBy, unfoldr)
import Data.Ord (comparing)
import Graphics.UI.SDL as SDL
import Graphics.UI.SDL.Primitives (filledCircle, line)
import System.Random.Mersenne.Pure64 (newPureMT, randomDouble)

res = 250

main = withInit [InitVideo] $ do
  w  <- setVideoMode res res 32 []
  ps <- randPoints
  enableEvent SDLMouseMotion False
  setCaption "Graham Scan" "Graham Scan"
  fillRect w (Just $ Rect 0 0 res res) $ Pixel 0
  limn w $ map (round *** round) $ hull ps
  plot w ps
  pause

plot w ps = do
  mapM_ (f . (round *** round)) ps
  SDL.flip w
 where
  f (x,y) = filledCircle w x y 1 $ Pixel 0xFFFFFFFF

limn w ps = f $ ps ++ [head ps]
 where
  f [_] = return ()
  f ((a,b):(x,y):ps) = do
    line w a b x y $ Pixel 0xFF0000FF
    f $ (x,y) : ps

pause = do
  delay 128
  e <- pollEvent
  case e of
   KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
   _                              -> pause

hull qs = go (drop 2 ps) $ reverse $ take 2 ps
 where
  o  = bottomRightP qs
  ps = o : sortBy (ccw o) (delete o qs)

go [] qs = qs

go (p:ps) s@(a:b:qs)
  | ccw a b p /= GT = go (p:ps) $ b:qs
  | otherwise       = go ps $ p:s

ccw (ax, ay) (bx, by) (cx, cy)
  | d < 0 = LT
  | d > 0 = GT
  | True  = EQ
 where
  d = (bx - ax) * (cy - ay) - (by - ay) * (cx - ax)

bottomRightP = maximumBy (comparing snd) . sort

randPoints = fmap f newPureMT
 where
  f = uncurry zip . splitAt 20 . g
  g = map (* res) . unfoldr (Just . randomDouble)

Wednesday, February 15, 2017

more Π obscurantism

Here's an illustration of an approach to pi pointed out by Kevin Brown. Define f(n) as the nearest greater or equal multiple of n-1, then of n-2, etc (yielding OEIS sequence 2491). Then, inverting a result found by Duane Broline and Daniel Loeb, pi = n2 / f(n).

But as you can see from the comment, the series converges very slowly!

main = print $ e**2 / f e e 1
 where
  e = 900000  -- yields 3.1416003...

f n 1 _ = n
f n k l = f n' (k - 1) $ n' / k
 where
  n' = head $ dropWhile (< n) $ map (* k) [l..]

Friday, June 3, 2016

segment intersection search

Here's a simple method I learned from Gareth Rees for finding the intersection (or parallel / collinear status) of two line segments; Rees credits Ronald Goldman, though I'd imagine the technique goes further back.
import Control.Arrow ((***))
import Control.Monad (liftM2)
import Data.Int (Int16)
import Data.List (unfoldr)
import Graphics.UI.SDL as SDL
import Graphics.UI.SDL.Primitives (line, circle)
import System.Random.Mersenne.Pure64 (pureMT, randomInt)

(xres, yres) = (1600, 900)
(red, white) = (0xFF0000FF, 0xFFFFFFFF)

main = withInit [InitVideo] $ do
  w <- setVideoMode xres yres 32 []
  enableEvent SDLMouseMotion False
  setCaption "Line Intersection" "Line Intersection"
  run w . take 15 $ filter f segments
 where
  f pq = dist pq < 1000
 
run w ss = do
  delay 256
  drawSegs w ss
  let ss' = map (vectorize . cst) ss
  drawPoints w $ liftM2 crossing ss' ss'
  SDL.flip w
  e <- pollEvent
  case e of
   KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
   _                              -> run w ss
  
drawSegs w = mapM_ f
 where
  f ((a,b), (c,d)) = line w a b c d $ Pixel red

drawPoints w = mapM_ (f . (round *** round))
 where
  f (x,y) = circle w x y 5 $ Pixel white 

-------------------------------------------------------

segments = h $ zip (f xres 13) $ f yres 23
 where
  f m      = map (fromIntegral . (`mod` m)) . g . pureMT
  g        = unfoldr $ fmap randomInt . Just
  h (p:ps) = (p, head ps) : h (tail ps)

dist (p,q) = sqrt $ (a-c)^2 + (b-d)^2
 where
  [a,b,c,d] = map f [fst p, snd p, fst q, snd q]
  f         = fromIntegral

crossing (p,r) (q,s)
  | f    = p .+. (t .^. r)  -- exactly one intersection
  | True = (-5,-5)          -- parallel, colinear, etc.
 where
  t   = (q .-. p) .*. s / (r .*. s)
  u   = (q .-. p) .*. r / (r .*. s)
  f   = (r .*. s) /= 0 && h t && h u
  h n = 0 <= n && n <= 1

(a,b) .+. (c,d) = (a+c, b+d)
(a,b) .-. (c,d) = (a-c, b-d)
n     .^. (x,y) = (n*x, n*y)   -- vector scaling
(a,b) .*. (c,d) =  a*d - b*c   -- 2D cross product

vectorize ((a,b), (c,d)) = ((a,b), (c-a, d-b))

cst (p,q) = (f p, f q)
 where
  f = fromIntegral *** fromIntegral

Monday, March 28, 2016

today in oblique approaches

If your standard library offers complex numbers and you're not in any hurry, you can't ask for much simpler (or more obscure) ways to compute pi than this.
import Data.Complex (magnitude, Complex((:+)))
 
pi = length $ takeWhile (< 2) $ map magnitude zs
 where
  zs = iterate (\z -> z^2 + ((-0.75) :+ 0.0000001)) 0

Wednesday, February 17, 2016

ikeda map

Continuing the theme of strange attractors, here's the well-known one embedded in the Ikeda map. The thumbnail at right shows the central 'vortex' of the attractor, and links to a larger viewport.

In these images, I've plotted the real and imaginary components along the x and y axes respectively. But the more popular way to visualize this attractor adds an extra parameter to the system and is expressed in trigonometric functions. Such adaptation of the code below yields these results.

{-# LANGUAGE NoMonomorphismRestriction #-}
 
import Data.Complex
import Graphics.UI.SDL as SDL
  
(xRes, yRes) = (1366, 768)
[a, b, k, p] = [0.85, 0.9, 0.4, 7.7]
 
main = withInit [InitVideo] $ do
  w <- setVideoMode xRes yRes 32 [NoFrame]
  s <- createRGBSurface [] 1 1 32 0 0 0 0
  fillRect s Nothing $ Pixel 0xFFFFFF
  enableEvent SDLMouseMotion False
  setCaption "Ikeda" "Ikeda"
  render w s $ ikeda $ 0 :+ 0
  SDL.flip w
  run w
 
run w = do
  e <- pollEvent
  delay 64
  case e of
   KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
   _                              -> run w
 
render w s = mapM_ draw
 where
  rect (x,y) = Just $ Rect (round x) (round y) 1 1
  g (x,y)    = (x + xRes / 24, y + yRes / 2.25)
  draw z     = blitSurface s Nothing w $ rect $ g p
   where
    p = (1050 * realPart z, 1050 * imagPart z)
 
ikeda = take 1500000 . filter g . iterate f
 where
  f z      = a + b * z * exp(i * (k-p) * (1 + abs z^2))
  i        = 0 :+ 1
  g (r:+i) = -35 < r && r < 35 && -18 < i && i < 18

Saturday, February 13, 2016

hénon attractor

French astronomer Michel Hénon reported on this chaotic dynamical system in 1976.
{-# LANGUAGE NoMonomorphismRestriction #-}

import Graphics.UI.SDL as SDL
 
(xRes, yRes) = (1366, 768)
(a, b)       = (1.4, 0.3)

main = withInit [InitVideo] $ do
  w <- setVideoMode xRes yRes 32 [NoFrame]
  s <- createRGBSurface [] 1 1 32 0 0 0 0
  fillRect s Nothing $ Pixel 0xFFFFFF
  enableEvent SDLMouseMotion False
  setCaption "Hénon" "Hénon"
  render w s henon
  SDL.flip w
  run w

run w = do
  e <- pollEvent
  delay 64
  case e of
   KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
   _                              -> run w

render w s = mapM_ draw
 where
  rect (x,y) = Just $ Rect (round x) (round y) 1 1
  g (x,y)    = (500*x + xRes / 2, 900*y + yRes / 2)
  draw (x,y) = blitSurface s Nothing w $ rect $ g (x,y)

henon = take 99999 $ filter g $ iterate f (0,0)
 where
  f (x,y) = (1 - a*x^2 + y, b*x)
  g (x,y) = -1.5 < x && x < 1.5 && -0.45 < y && y < 0.45

Tuesday, January 26, 2016

langton's ant

For a round of code golf, I wrote this spare implementation of Chris Langton's remarkably simple universal computer. If you want amenities like pause, random starting pattern or even quit, check out this more complete version.
import Control.Arrow ((***))
import Control.Monad (when, join)
import Data.Set (insert, delete, member, empty)
import Graphics.UI.SDL as SDL

(xres, yres, sq) = (1366, 768, 4)

main = withInit [InitVideo] $ do
  w <- setVideoMode xres yres 32 [NoFrame]
  enableEvent SDLMouseMotion False
  setCaption "Langton's Ant" "Langton's Ant"
  run w [1..] p (0,1) empty
 where
  p = (xres `div` 2 `div` sq, yres `div` 2 `div` sq)

run w (n:ns) p v ps = do
  when (n `mod` 7 == 0) render
  run w ns (move p $ g v) (g v) $ f p ps
 where
  b          = member p ps
  (f,g)      = if b then (delete, fl) else (insert, fr)
  fr (x,y)   = if x == 0 then (-y,x) else (y,x)
  fl (x,y)   = if x == 0 then (y,x)  else (y,-x)
  move (x,y) = (x+) *** (+y)
  render     = do
    fillRect w (Just $ Rect 0 0 xres yres) $ Pixel 0
    mapM_ (draw w . join (***) (* sq)) ps
    SDL.flip w
 
draw w p = f p =<< g [SWSurface] sq sq 32 0 0 0 0
 where
  rect x y  = Just $ Rect x y sq sq
  g         = createRGBSurface
  f (x,y) s = do fillRect s (rect 0 0) $ Pixel 0xFFFFFF
                 blitSurface s (rect 0 0) w $ rect x y

Wednesday, December 23, 2015

connett circles

Like Barry Martin's 'Hopalong' fractal, this dynamical system from John Connett was first published in Scientific American in 1986. However, it's not actually a fractal: when zooming in (by clicking two points to specify a rectangle), viewers will see the pattern is not infinitely detailed. Instead, peacock-like images appear.
{-# LANGUAGE NoMonomorphismRestriction #-}

import Data.Complex (magnitude, Complex((:+)))
import Graphics.UI.SDL as SDL
 
(xres, yres, cast) = (768, 768, fromIntegral)
 
main = withInit [InitVideo] $ do
  w <- setVideoMode xres yres 32 [NoFrame]
  enableEvent SDLMouseMotion False
  setCaption "Connett Circles" "Connett Circles"
  zoom w (-2000, -2000) (2000, 2000)
 
zoom w (x,y) (s,t) = do
  let r = area x s y t
  render w r
  run w r (0,0)

run w (xs,ys) p@(x,y) = do
  delay 32
  e <- pollEvent
  case e of
   KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
   MouseButtonUp i j _            -> click $ f i j
   _                              -> run w (xs,ys) p
 where
  f i j = (xs !! cast i, ys !! cast j)
  click q
    | p == (0,0) = run w (xs,ys) q
    | otherwise  = zoom w p q

render w r = draw w . concat . f $ set r
 where
  f = zipWith zip $ map (zip [0..] . repeat) [0..]

draw w r = do
  s <- createRGBSurface [] 1 1 32 0 0 0 0
  mapM_ (draw s) r
  SDL.flip w
 where
  rect (x,y)   = Just $ Rect x y 1 1
  draw s (p,n) = do fillRect s Nothing $ Pixel $ rgb n
                    blitSurface s Nothing w $ rect p

------------------------------------------------------

rgb n 
  | m `mod` 7 == 0 = 0xFFFFFF
  | m `mod` 6 == 0 = 0x00FFFF
  | m `mod` 5 == 0 = 0xFFFF00
  | m `mod` 4 == 0 = 0xFF0000
  | m `mod` 3 == 0 = 0xFF00
  | m `mod` 2 == 0 = 0xFF
  | otherwise      = 0
 where
  f x y = round $ magnitude $ (cast x^2) :+ (cast y^2)
  m     = abs n

set (xs, ys) = [[f x y | x <- xs] | y <- ys]
 where
  f x y = round $ magnitude $ (x^2) :+ (y^2)

area x x' y y' = ([t, t+a.. s], [u, u-b.. v])
 where
  (a,b)    = ((s - t) / dx, (u - v) / dy)
  (s,t)    = (max x x', min x x')
  (u,v)    = (max y y', min y y')
  (dx, dy) = (cast xres, cast yres)

martin attractor

This pattern generator, discovered by Barry Martin, was nicknamed 'Hopalong' when Scientific American introduced it in their September '86 issue.

Clicking the window adjusts the viewport position; there is also an alternate version with color and animation

Update: For a certain Rubyist friend, I wrote another lazily-evaluated, colored and animated implementation in Ruby.

import Control.Arrow ((***))
import Graphics.UI.SDL as SDL

(xres, yres) = (850, 850)
(a,b,c)      = (5, 1, 20)

main = withInit [InitVideo] $ do
  w <- setVideoMode xres yres 32 [NoFrame]
  enableEvent SDLMouseMotion False
  setCaption "Hopalong Pattern" "Barry Martin"
  render w p ps
  loop w p ps
 where
  p  = (xres `div` 2, yres `div` 2)
  ps = take 7000000 points

loop w p ps = do
  delay 32
  e <- pollEvent
  case e of
   KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
   MouseButtonUp x y _            -> click $ f x y
   _                              -> loop w p ps
 where
  click p = render w p ps >> loop w p ps
  f x y   = (fromIntegral x, fromIntegral y)

render w (x,y) ps = do
  s <- createRGBSurface [SWSurface] 1 1 32 0 0 0 0
  fillRect w (Just $ Rect 0 0 xres yres) (Pixel 0)
  mapM_ (draw s 0xFFFFFF . rect . center) ps
  SDL.flip w
 where
  center     = ((xres - x) +) *** (+ (yres - y))
  rect (x,y) = Just $ Rect x y 1 1
  draw s c p = do fillRect s Nothing $ Pixel c
                  blitSurface s Nothing w p

points = map (round *** round) $ iterate f (0,0)
 where
  f (x,y) = (g x y, a - x)
  g x y   = y - (signum x * (sqrt . abs $ b*x - c))

Tuesday, December 22, 2015

lyapunov fractals

It took me some experimentation to figure out how to color this derivation of the logistic map; I'm still not quite sure how the hues should scale as you zoom. But the bi-tonal method shown below works well enough to produce the image at right.
{-# LANGUAGE NoMonomorphismRestriction #-}

import Control.Monad (liftM2)
import Data.Bits (shiftL)
import Graphics.UI.SDL as SDL
 
(xres, yres, cast) = (400, 400, fromIntegral)

main = withInit [InitVideo] $ do
  w <- setVideoMode xres yres 32 [NoFrame]
  enableEvent SDLMouseMotion False
  setCaption "Lyapunov" "Lyapunov"
  render w 1024 (3.2, 3.2) (3.8, 3.8)

render w i (x,y) (s,t) = do
  draw w (zip ps $ area x s y t) i
  run w
 where
  ps = liftM2 (,) [0.. xres] [0.. yres]

draw w cs i = do
  s <- createRGBSurface [] 1 1 32 0 0 0 0
  mapM_ (draw s) cs
  SDL.flip w
 where
  draw s (p,q) = do 
    let rect (x,y) = Just $ Rect x y 1 1
    fillRect s Nothing $ rgb $ lyapunov i q
    blitSurface s Nothing w $ rect p

run w = do
  delay 32
  e <- pollEvent
  case e of
   KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
   _                              -> run w

-----------------------------------------------------

rgb n | n < 0 = Pixel $ f n
      | True  = Pixel $ f n `shiftL` 16
 where
  f x = round (min 255 $ abs n * 255)

lyapunov i p = (1 / cast i) * sum rs
 where
  rs = filter (not . isInfinite) $ f xs
  f  = map (log . abs) . zipWith ($) (str p)
  xs = map ((1 -) . (* 2)) $ verhulst i p

verhulst i p = take i $ go 0.5 $ str p
 where
  go x (f:fs) = x : go (f (x * (1-x))) fs

str (x,y) = map (*) $ concat $ repeat [x,x,y,x,y]

area x x' y y' = liftM2 (,) [s, s+a.. t] [u, u+b.. v]
 where
  (a,b) = ((t - s) / m, (v - u) / n)
  (s,t) = (min x x', max x x')
  (u,v) = (min y y', max y y')
  (m,n) = (cast xres, cast yres)

Tuesday, December 15, 2015

the mandelbrot set

Benoit's famous fractal. Click any two points to specify a zoom rectangle.
import Data.Complex (magnitude, Complex((:+)))
import Graphics.UI.SDL as SDL
 
(xres, yres) = (900, 750)

main = withInit [InitVideo] $ do
  w <- setVideoMode xres yres 32 [NoFrame]
  enableEvent SDLMouseMotion False
  setCaption "Mandelbrot" "Mandelbrot"
  zoom w 1 (-2, 1) (0.5, -1)

zoom w i (x,y) (s,t) = do
  let r = area x s y t
  render w r i
  run w (i + 1) r (0,0)

render w r i = draw w . concat . f $ set i r
 where
  f = zipWith zip $ map (zip [0..] . repeat) [0..]

draw w cs = do
  s <- createRGBSurface [] 1 1 32 0 0 0 0
  mapM_ (draw s) cs
  SDL.flip w
 where
  rgb n        = div (floor $ min n 2 * 512) 2
  rect (x,y)   = Just $ Rect x y 1 1
  draw s (p,n) = do fillRect s Nothing $ Pixel $ rgb n
                    blitSurface s Nothing w $ rect p

run w i (xs,ys) p@(x,y) = do
  e <- pollEvent
  case e of
   KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
   MouseButtonUp i j _            -> click $ f i j
   _                              -> run w i (xs,ys) p
 where
  f i j = (xs !! fromIntegral i, ys !! fromIntegral j)
  click q
    | p == (0,0) = run w i (xs,ys) q
    | otherwise  = zoom w i p q

------------------------------------------------------

area x x' y y' = ([t, t+a.. s], [u, u-b.. v])
 where
  (a,b)    = ((s - t) / dx, (u - v) / dy)
  (s,t)    = (max x x', min x x')
  (u,v)    = (max y y', min y y')
  (dx, dy) = (fromIntegral xres, fromIntegral yres)

set i (xs, ys) = [[f x y | x <- xs] | y <- ys]
 where
  f x = magnitude . g . (x :+)
  g a = iterate ((+ a) . (^2)) 0 !! (i * 20)

Monday, November 2, 2015

a lava lamp

This 'lava lamp' is actually a simple cyclic cellular automaton. The rule is courtesy of Jason Rampe.
import Control.Arrow ((***))
import Control.Monad (join, liftM2)
import Data.List (unfoldr)
import Data.Maybe (fromJust)
import qualified Data.HashMap.Strict as M
import Graphics.UI.SDL as SDL
import System.Random.Mersenne.Pure64 (newPureMT, randomInt)
 
(xres, yres, unit) = (160, 90, 10)
 
main = withInit [InitVideo] $ do
  w <- setVideoMode (f xres) (f yres) 32 [NoFrame]
  enableEvent SDLMouseMotion False
  setCaption "Cyclic CA" "Cyclic CA"
  pause w =<< randPattern
 where
  f = (* unit)
 
pause w cs = do
  delay 128
  e <- pollEvent
  case e of
   KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
   KeyUp (Keysym SDLK_SPACE _ _)  -> run w cs
   _                              -> pause w cs

run w cs = do
  drawCells w cs
  SDL.flip w
  e <- pollEvent
  case e of
   KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
   KeyUp (Keysym SDLK_SPACE  _ _) -> pause w cs
   _                              -> run w cs'
 where
  cs' = M.mapWithKey (next cs) cs

drawCells w cs = do
  s <- createRGBSurface [] unit unit 32 0 0 0 0
  mapM_ (draw s) $ M.toList cs
 where
  rgb n        = [0xFF0000, 0xAA0000, 0x770000] !! n
  rect (x,y)   = Just $ Rect x y unit unit
  scale        = rect . join (***) (* unit)
  draw s (p,n) = do fillRect s Nothing $ Pixel $ rgb n
                    blitSurface s Nothing w $ scale p

randPattern = fmap (M.fromList . zip ps) ns
 where
  g  = map (`mod` 3) . unfoldr (Just . randomInt)
  ns = (return . g) =<< newPureMT
  ps = liftM2 (,) [0.. xres - 1] [0.. yres - 1]

next cs p n = if m >= 10 then n' else n
 where
  m    = length $ filter (== n') area
  area = map (fromJust . (`M.lookup` cs)) $ moore p
  n'   = (n + 1) `mod` 3

moore (x,y) = filter p ps
 where
  ps      = liftM2 (,) [x-2.. x+2] [y-2.. y+2]
  p (x,y) = x >= 0 && x < xres && y >= 0 && y < yres

Wednesday, August 19, 2015

random undirected graphs

This little program generates random, undirected graphs without loops or isolated vertices - though the graph is not necessarily connected.
import Data.List (nub, unfoldr)
import Control.Arrow ((***))
import Control.Monad (join, replicateM)
import Graphics.UI.SDL as SDL
import Graphics.UI.SDL.Primitives (line)
import System.Random.Shuffle (shuffle')
import System.Random.Mersenne.Pure64
       (newPureMT, randomInt)

[xres, yres, unit] = map (1 *) [512, 512, 5]

main = withInit [InitVideo] $ do
  window <- setVideoMode xres yres 32 [NoFrame]
  enableEvent SDLMouseMotion False
  setCaption "Random Graphs" "Random Graphs"
  render window =<< genGraph
  run

run = do
  delay 128
  e <- pollEvent
  case e of
   KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
   _                              -> run

render w ls = do
  s <- createRGBSurface [SWSurface] unit unit 32 0 0 0 0
  fillRect s Nothing $ Pixel 0xFFFFFF
  fillRect w (Just $ Rect 0 0 xres yres) (Pixel 0)
  mapM_ (draw w s) ls
  SDL.flip w

draw w s ((a,b), (x,y)) = do
  blitSurface s Nothing w $ rect (a,b)
  blitSurface s Nothing w $ rect (x,y)
  line w a' b' x' y' $ Pixel 0xFF0000FF
 where
  rect (x,y)       = Just $ Rect x y unit unit
  [a', b', x', y'] = map ((+2) . fromIntegral) [a,b,x,y]

genGraph = do
  gs <- replicateM 2 newPureMT
  ps <- return . nub . take 15 $ rands gs
  qs <- permute ps
  return $ filter p $ zip ps qs
 where
  permute xs       = fmap (shuffle' xs $ length xs) newPureMT
  p ((a,b), (x,y)) = (a,b) /= (y,x) && (a,b) /= (x,y)

rands [g, g'] = map (f . h) $ zip xs ys
 where
  xs    = map ((8+) . (`mod` x)) $ unfoldr (Just . randomInt) g
  ys    = map ((8+) . (`mod` y)) $ unfoldr (Just . randomInt) g'
  (x,y) = (xres - 16, yres - 16)
  (f,h) = (join (***) (* unit), join (***) (`div` unit))