fay-0.24.2.0: examples/oscillator.hs

{-# LANGUAGE EmptyDataDecls    #-}

module RingOscillator (main) where

import FFI

-- System parameters.
--
data Params = Params { alpha :: Double
                     , omega :: Double
                     , deven :: Double
                     , dodd :: Double
                     , beta :: Double
                     , sigma :: Double
                     , nosc :: Int } deriving Show

-- Update functions for system parameters used in the event handlers
-- for the parameter select lists.  Not pretty, but there's not really
-- a nicer way to do it.
--
alphaUpd, omegaUpd :: Params -> Double -> Params
devenUpd, doddUpd :: Params -> Double -> Params
betaUpd, sigmaUpd :: Params -> Double -> Params
noscUpd :: Params -> Int -> Params
alphaUpd p val = p { alpha = val }
omegaUpd p val = p { omega = val }
devenUpd p val = p { deven = val }
doddUpd p val = p { dodd = val }
betaUpd p val = p { beta = val }
sigmaUpd p val = p { sigma = val }
noscUpd p val = p { nosc = val }

-- Need this to be able to make a 'Ref Params'.
--

-- Default values for the system parameters (taken from Bridges &
-- Reich 2001).
--
defaultParams = Params { alpha = 1, omega = 1.8, deven = 0.0075, dodd = 0.0125,
                         beta = 1, sigma = 4, nosc = 5 }


-- The system state is represented as: [y, x1, ..., xN, y', x1', ..., xN']
-- where y is the position of the external forcing oscillator, x1..xN
-- are the positions of the oscillators in the ring, and y' and
-- x1'..xN' are the time derivatives.  The initial state has some
-- non-zero displacements to get us going and zero velocities.
--
initialState :: Params -> [Double]
initialState p = replicate (nosc p + 1) 1 ++ replicate (nosc p + 1) 0

-- Default initial state.
--
ic = centre defaultParams (initialState defaultParams)

-- Make data storage for scrolling graph view.
--
makeGraphData :: Int -> Int -> Fay (Double,[Buffer])
makeGraphData no size = do
  bufs <- replicateM no (newBuf size)
  return (0,bufs)


-- Calculate system right hand side, i.e. the equations that determine
-- the time derivative of the system state.
--
rhs :: Params -> [Double] -> [Double]
rhs p state = [dy] ++ dxs ++ [dy'] ++ dxs'
  where n = nosc p

        -- Split system state: Fay doesn't yet allow us to to
        -- something like ((y:xs),(y':xs')) = splitAt (n+1) state, so
        -- we just decompose things manually.
        ss = splitAt (n + 1) state
        sss = fst ss ; sss' = snd ss
        y = head sss ; xs = tail sss
        y' = head sss' ; xs' = tail sss'

        -- Time derivatives of y and the xs are immediately
        -- accessible.
        dy = y' ; dxs = xs'

        -- Forcing system right hand side.
        dy' = 0 - alpha p * (y^2 - 1) * y' - ((omega p)^2::Double) * y

        -- Right hand side for ring oscillators:

        -- Displacement differences.
        xdiffs = (head xs - (xs !! (n - 1))) : zipWith (-) (tail xs) xs

        -- Nonlinear inter-oscillator potential calculation:
        -- V(x) = 1/2 x^2 + 1/4 x^4  =>  V'(x) = x + x^3
        vprime x = x + x^3
        vp = map vprime xdiffs

        -- Differences between adjacent potential values.
        vfacs = zipWith (-) vp (tail vp) ++ [vp !! (n - 1) - head vp]

        -- Damping factors alternate between even and odd numbered
        -- oscillators.
        ds = cycle [dodd p, deven p]

        -- Calculate basic time derivative for each oscillator.
        dxstmp = zipWith3 (\d x' vf -> 0 - d * x' - beta p * vf) ds xs' vfacs

        -- Add forcing term for first oscillator in ring.
        dxs' = (head dxstmp + sigma p * y) : tail dxstmp


-- Take a single fourth order Runge-Kutta step for an autonomous
-- system.  The arguments are the right hand side function, a system
-- state and a time step.
--
rk4 :: ([Double] -> [Double]) -> [Double] -> Double -> [Double]
rk4 f yn h = zipWith5 (\y a b c d -> y+(a+2*b+2*c+d)/6) yn k1 k2 k3 k4
  where k1 = map mh $ f yn
        k2 = map mh $ f (zipWith (+) yn (map half k1))
        k3 = map mh $ f (zipWith (+) yn (map half k2))
        k4 = map mh $ f (zipWith (+) yn k3)
        mh x = h * x
        half x = 0.5 * x


-- Because the *velocities* are damped in our model, but not the
-- displacements, there tends to be a more or less linear drift in the
-- average position of the oscillators over time.  This isn't a
-- problem for displaying the positions on the ring, but it's
-- inconvenient for displaying time traces of the positions as a
-- graph.  To make this a bit easier, we centre the oscillator values
-- for graph display by subtracting the mean displacement for the ring
-- oscillators for each ring oscillator.  The forcing displacement and
-- the time derivative values are not affected by this centring, and
-- in fact none of the time derivative calculations for the system are
-- affected either, since they all depend exclusively on *differences*
-- between displacements, not the absolute values of the
-- displacements.
--
centre :: Params -> [Double] -> [Double]
centre p s = [y] ++ xsc ++ [y'] ++ xs'
  where ss = splitAt (nosc p + 1) s
        sss = fst ss ; sss' = snd ss
        y = head sss ; xs = tail sss
        y' = head sss' ; xs' = tail sss'
        m = (sum xs) / fromIntegral (length xs)
        xsc = map (\x->x-m) xs


-- Display dimensions.
--
ww, wh :: Double                -- Ring canvas dimensions.
ww = 400 ; wh = 400

gww, gwh :: Double              -- Graph canvas dimensions.
gww = 640 ; gwh = 200

-- Other display parameters.
--
dt = 0.025                      -- Integration time step.
framems = 30                    -- Milliseconds between frames.
renderrng = 6.0                 -- Nominal range for graph and ring displays.

-- Sizing for ring view elements.
--
rdotfac, rinner, rdot, fyoff :: Double
rdotfac = 0.0375
rinner = wh / (2 + 7 * rdotfac)
rdot = rdotfac * rinner
fyoff = rinner + 3 * rdot
ctr :: Point
ctr = (0.5 * ww, 5 * rdot + rinner) -- Centre point of ring display.

-- Sizing for graph view elements.
--
gwbuf, ticklen, gwwtime :: Double
gwbuf = 20                      -- Buffer at right end of graph.
ticklen = 10                    -- Graph time tick size.
gwwtime = 10                    -- Time range displayed in graph.
pxpert, pxpersamp :: Double
pxpert = gww / gwwtime          -- Pixels per time unit.
pxpersamp = dt * pxpert         -- Pixels per sample.
nsamp :: Int
nsamp = floor ((gww - gwbuf) / pxpersamp) -- Samples in graph.


-- Main entry point.  Just register an "onload" handler.
--
main :: Fay ()
main = addWindowEventListener "load" run

-- Main function: runs at load time.
--
run :: Event -> Fay Bool
run _ = do
  -- Get DOM elements and canvas contexts.
  [can,graph] <- mapM getElementById ["canvas","graph"]
  [go,stop,reset] <- mapM getElementById ["go","stop","reset"]
  sels <- mapM getElementById ["alpha","omega","deven","dodd","beta","sigma"]
  noscSel <- getElementById "nosc"
  [c,cg] <- mapM (\c -> getContext c "2d") [can,graph]

  -- Set up references to parameter values, current state, a timer
  -- identifier and the data required for rendering the graph view.
  pref <- newRef defaultParams
  xref <- newRef ic
  timerref <- newRef (Nothing :: Maybe Int)
  gd <- makeGraphData (nosc defaultParams) nsamp
  gdataref <- newRef gd

  -- Render initial ring and graph views.
  render c defaultParams ic renderrng
  renderGraph cg pref gdataref (centre defaultParams ic) renderrng

  -- Event listeners for buttons and parameter value selection.
  addEventListener go "click" $
    doGo timerref (animate c cg pref xref gdataref renderrng) framems
  addEventListener stop "click" $ doStop timerref
  addEventListener reset "click" $ doReset c cg timerref pref xref gdataref sels
  forM_ (zip sels [alphaUpd,omegaUpd,devenUpd,doddUpd,betaUpd,sigmaUpd]) $ \(s, pfn) -> do
    addEventListener s "change" $ doParamChange pref pfn
  addEventListener noscSel "change" $
    doNoscChange c cg timerref pref xref gdataref
  return False


-- If the simulation isn't currently running (marked by the timer
-- reference being Nothing), start it by causing the animation
-- function to be called at the appropriate interval.
--
doGo :: Ref (Maybe Int) -> Fay () -> Double -> Event -> Fay Bool
doGo tref anim interval _ = do
  oldtimer <- readRef tref
  case oldtimer of
    Nothing -> do
      timer <- setInterval anim interval
      writeRef tref (Just timer)
    Just _ -> return ()
  return False

-- If the simulation is running (marked by the timer reference having
-- a Just value), stop it.
--
doStop :: Ref (Maybe Int) -> Event -> Fay Bool
doStop tref _ = do
  oldtimer <- readRef tref
  case oldtimer of
    Nothing -> return ()
    Just timer -> do
      clearInterval timer
      writeRef tref Nothing
  return False

-- Stop the simulation and set everything back to the state it was at
-- the beginning, including all parameter values (the default values
-- are the middle ones out of five possibilities).
--
doReset :: Context -> Context -> Ref (Maybe Int) -> Ref Params ->
           Ref [Double] -> Ref (Double,[Buffer]) -> [Element] ->
           Event -> Fay Bool
doReset c cg tref pref xref gdataref sels e = do
  doStop tref e
  writeRef pref defaultParams
  writeRef xref ic
  gd <- makeGraphData (nosc defaultParams) nsamp
  writeRef gdataref gd
  render c defaultParams ic renderrng
  renderGraph cg pref gdataref (centre defaultParams ic) renderrng
  forM_ sels $ \s -> setSelectIndex s 2
  return False

-- Simple parameter change: just get the selected value and update the
-- parameter reference.
--
doParamChange :: Ref Params -> (Params -> Double -> Params) -> Event -> Fay Bool
doParamChange pref updfn e = do
  target <- eventTarget e
  sval <- selectValue target
  p <- readRef pref
  writeRef pref (updfn p (parseDouble sval))
  return False

-- A change in number of oscillators is a bit more complicated.  We
-- stop the simulation and reset everything to start over -- there's
-- no obvious way to derive a state with a different number of
-- oscillators from the current state.
--
doNoscChange :: Context -> Context ->
                Ref (Maybe Int) -> Ref Params ->
                Ref [Double] -> Ref (Double,[Buffer]) -> Event -> Fay Bool
doNoscChange c cg timerref pref xref gdataref e = do
  doStop timerref e
  target <- eventTarget e
  sval <- selectValue target
  let newNosc = parseInt sval
  p <- readRef pref
  let pnew = p { nosc = newNosc }
  let xnew = centre pnew (initialState pnew)
  writeRef pref pnew
  writeRef xref xnew
  gd <- makeGraphData newNosc nsamp
  writeRef gdataref gd
  render c pnew xnew renderrng
  renderGraph cg pref gdataref xnew renderrng
  return False


-- Animation function.  Take a single step of the ODE system and
-- re-render the ring and the graph views.
--
animate :: Context -> Context -> Ref Params -> Ref [Double] ->
           Ref (Double,[Buffer]) -> Double -> Fay ()
animate c cg pref xref gdataref rng = do
  p <- readRef pref
  x <- readRef xref
  let newx = rk4 (rhs p) x dt
  writeRef xref newx
  render c p newx rng
  renderGraph cg pref gdataref (centre p newx) rng


-- Render the ring view.
--
render :: Context -> Params -> [Double] -> Double -> Fay ()
render c p s rng = do
  -- Extract the forcing and the ring oscillator displacements.
  let f = head s
  let xs = take (nosc p) $ tail s

  -- Draw "furniture": ring and baseline for forcing oscillator.
  clearRect c (0,0) (ww, wh)
  beginPath c
  arc c ctr rinner 0 (2*pi)
  let fscale = rinner * 2 * pi / 5
  moveTo c $ offset ctr (-0.5*fscale,-fyoff)
  lineTo c $ offset ctr (0.5*fscale,-fyoff)
  moveTo c $ offset ctr (0, -dtickin)
  lineTo c $ offset ctr (0, -dtickout)
  setLineWidth c 2
  setStrokeStyle c "grey"
  stroke c

  -- Blobs for ring oscillators: first, forced oscillator is picked
  -- out in blue.
  forM_ [0..(nosc p)-1] (\i -> do
    let th0 = fromIntegral i * 2 * pi / fromIntegral (nosc p)
    let th = th0 + (xs !! i) / rng * 2 * pi / fromIntegral (nosc p)
    let dctr = offset ctr (rinner * sin th, -rinner * cos th)
    beginPath c
    arc c dctr rdot 0 (2*pi)
    setFillStyle c (if i == 0 then "blue" else "red")
    fill c)

  -- Blob for forcing oscillator.
  let dctr = offset ctr (f / rng * fscale, -fyoff)
  beginPath c
  arc c dctr rdot 0 (2*pi)
  setFillStyle c "grey"
  fill c
  where dtickin = fyoff - 0.5 * ticklen
        dtickout = fyoff + 0.5 * ticklen


-- Render the graph view.
--
renderGraph :: Context -> Ref Params -> Ref (Double,[Buffer]) ->
               [Double] -> Double -> Fay ()
renderGraph cg pref gdataref x rng = do
  -- Get values for parameters and time and state list.
  p <- readRef pref
  (ts,bufs) <- readRef gdataref

  -- Graph "furniture": axes and time ticks.
  setFont cg "10pt sans-serif"
  setStrokeStyle cg "grey"
  setFillStyle cg "grey"
  setLineWidth cg 2
  clearRect cg (0,0) (gww,gwh)
  beginPath cg
  moveTo cg (0,gwh/2)
  lineTo cg (gww,gwh/2)
  moveTo cg (1,0)
  lineTo cg (1,gwh)
  let ticks = takeWhile (\t -> t <= floor (ts + gwwtime)) [floor ts + 1..]
  let tickxs = map (\t -> (fromIntegral t - ts) * pxpert) ticks
  forM_ (zip ticks tickxs) $ \(t,x) -> do
    moveTo cg (x,gwh/2-ticklen/2)
    lineTo cg (x,gwh/2+ticklen/2)
    let txt = show t
    txtw <- measureText cg txt
    fillText cg txt (x-txtw/2,gwh/2+2*ticklen) Nothing
  stroke cg

  -- Add new samples to circular buffers and draw traces: grey for
  -- forcing, blue for forced ring oscillator, red for the other ring
  -- oscillators.
  setLineWidth cg 1
  newbufs <- forM (zip3 x bufs ("grey":"blue":repeat "red")) $ \(newx,buf,col) -> do
    newbuf <- bufAdd buf newx
    beginPath cg
    setStrokeStyle cg col
    y0 <- bufVal newbuf 0
    moveTo cg (0, gwh/2*(1-y0/renderrng))
    when (bufCurSize buf > 1) $ forM_ [1..bufCurSize buf-1] $ \i -> do
      y <- bufVal newbuf i
      lineTo cg (fromIntegral i * pxpersamp, gwh/2*(1-y/renderrng))
    stroke cg
    return newbuf

  -- Update the graph data reference.
  writeRef gdataref (ts+if bufCurSize (head newbufs) < nsamp then 0 else dt,newbufs)


--------------------------------------------------------------------------------
-- Utilities

zipWith5 :: (a->b->c->d->e->f) -> [a]->[b]->[c]->[d]->[e]->[f]
zipWith5 z (a:as) (b:bs) (c:cs) (d:ds) (e:es) = z a b c d e :
                                                zipWith5 z as bs cs ds es
zipWith5 _ _ _ _ _ _ = []

replicateM :: Int -> Fay a -> Fay [a]
replicateM n x = sequence (replicate n x)

parseDouble :: String -> Double
parseDouble = ffi "parseFloat(%1)"

parseInt :: String -> Int
parseInt = ffi "parseInt(%1)"

--------------------------------------------------------------------------------
-- DOM

data Element
instance Show Element

type Size = (Int,Int)

getElementById :: String -> Fay Element
getElementById = ffi "document['getElementById'](%1)"

data Event

addWindowEventListener :: String -> (Event -> Fay Bool) -> Fay ()
addWindowEventListener = ffi "window['addEventListener'](%1,%2,false)"

addEventListener :: Element -> String -> (Event -> Fay Bool) -> Fay ()
addEventListener = ffi "%1['addEventListener'](%2,%3,false)"

setInterval :: Fay () -> Double -> Fay Int
setInterval = ffi "window['setInterval'](%1,%2)"

clearInterval :: Int -> Fay ()
clearInterval = ffi "window['clearInterval'](%1)"

print :: a -> Fay ()
print = ffi "console.log(%1)"

eventTarget :: Event -> Fay Element
eventTarget = ffi "%1['target']"

selectValue :: Element -> Fay String
selectValue = ffi "%1[%1['selectedIndex']]['value']"

setSelectIndex :: Element -> Int -> Fay ()
setSelectIndex = ffi "%1['selectedIndex']=%2"

--------------------------------------------------------------------------------
-- Ref

-- | A mutable reference like IORef.
data Ref a

-- | Make a new mutable reference.
newRef :: a -> Fay (Ref a)
newRef = ffi "new Fay$$Ref(%1)"

-- | Replace the value in the mutable reference.
writeRef :: Ref a -> a -> Fay ()
writeRef = ffi "Fay$$writeRef(%1,%2)"

-- | Get the referred value from the mutable value.
readRef :: Ref a -> Fay a
readRef = ffi "Fay$$readRef(%1)"

--------------------------------------------------------------------------------
-- Canvas API

data Context
instance Show Context

type Point = (Double,Double)
type Dim = (Double,Double)

getContext :: Element -> String -> Fay Context
getContext = ffi "%1.getContext(%2)"

offset :: Point -> Dim -> Point
offset (x,y) (dx,dy) = (x+dx,y+dy)

-- Basic attributes

setFillStyle :: Context -> String -> Fay ()
setFillStyle = ffi "%1['fillStyle']=%2"

setFont :: Context -> String -> Fay ()
setFont = ffi "%1['font']=%2"

setLineWidth :: Context -> Double -> Fay ()
setLineWidth = ffi "%1['lineWidth']=%2"

setStrokeStyle :: Context -> String -> Fay ()
setStrokeStyle = ffi "%1['strokeStyle']=%2"

-- Path methods

arc :: Context -> Point -> Double -> Double -> Double -> Fay ()
arc c (x,y) r beg end = arc' c x y r beg end True

arcC :: Context -> Point -> Double -> Double -> Double -> Fay ()
arcC c (x,y) r beg end = arc' c x y r beg end False

arc' :: Context -> Double -> Double -> Double -> Double -> Double ->
        Bool -> Fay ()
arc' = ffi "%1['arc'](%2,%3,%4,%5,%6,%7)"

beginPath :: Context -> Fay ()
beginPath = ffi "%1['beginPath']()"

clip :: Context -> Fay ()
clip = ffi "%1['clip']()"

closePath :: Context -> Fay ()
closePath = ffi "%1['closePath']()"

fill :: Context -> Fay ()
fill = ffi "%1['fill']()"

lineTo :: Context -> Point -> Fay ()
lineTo c (x,y) = lineTo' c x y

lineTo' :: Context -> Double -> Double -> Fay ()
lineTo' = ffi "%1['lineTo'](%2,%3)"

moveTo :: Context -> Point -> Fay ()
moveTo c (x,y) = moveTo' c x y

moveTo' :: Context -> Double -> Double -> Fay ()
moveTo' = ffi "%1['moveTo'](%2,%3)"

stroke :: Context -> Fay ()
stroke = ffi "%1['stroke']()"

-- Rectangles

clearRect :: Context -> Point -> Dim -> Fay ()
clearRect c (x,y) (w,h) = clearRect' c x y w h

clearRect' :: Context -> Double -> Double -> Double -> Double -> Fay ()
clearRect' = ffi "%1['clearRect'](%2,%3,%4,%5)"

-- Text

fillText :: Context -> String -> Point -> Maybe Double -> Fay ()
fillText c s (x,y) Nothing = fillText1 c s x y
fillText c s (x,y) (Just mw) = fillText2 c s x y mw

fillText1 :: Context -> String -> Double -> Double -> Fay ()
fillText1 = ffi "%1['fillText'](%2,%3,%4)"

fillText2 :: Context -> String -> Double -> Double -> Double -> Fay ()
fillText2 = ffi "%1['fillText'](%2,%3,%4,%5)"

measureText :: Context -> String -> Fay Double
measureText = ffi "%1['measureText'](%2)['width']"


--------------------------------------------------------------------------------
-- Circular buffers

data Array

data Buffer = Buffer { bufSize :: Int
                     , bufCurSize :: Int
                     , bufNext :: Int
                     , bufArr :: Array }

newBuf :: Int -> Fay Buffer
newBuf size = do
  arr <- newArray size
  return $ Buffer size 0 0 arr

bufAdd :: Buffer -> Double -> Fay Buffer
bufAdd (Buffer sz cursz nxt arr) x = do
  let cursz' = if cursz < sz then cursz + 1 else sz
  setArrayVal arr nxt x
  let nxt' = (nxt + 1) `rem` sz
  return $ Buffer sz cursz' nxt' arr

bufVal :: Buffer -> Int -> Fay Double
bufVal (Buffer sz cursz nxt arr) i = do
  let idx = (if cursz < sz then i else nxt + i) `rem` sz
  arrayVal arr idx >>= return

newArray :: Int -> Fay Array
newArray = ffi "new Array(%1)"

setArrayVal :: Array -> Int -> Double -> Fay ()
setArrayVal = ffi "%1[%2]=%3"

arrayVal :: Array -> Int -> Fay Double
arrayVal = ffi "%1[%2]"