I recently started learning Haskell and as my first project I decided to port a particle simulation I had written in C.
The simulation is pretty simple. We have a cubic box with particles (spheres) inside it. Then we choose to either move a particle or to change the volume of the box and check for any collisions between the particles. If there are any collisions we keep the old configuration.
To speed up collision detection I have also implemented a 'cell linked list' which is in principle a uniform grid with a minimum size of a cell equal to the particle diameter. In the C code I have implemented this using arrays while in the Haskell code I use a Vector of Lists.
You can find my full code here, but here are some of the most important snippets,
The main Simulation Loop, where type Eval a = ReaderT (Env a) (StateT (SimState a) IO)
runSimulation :: Int -> Eval Double ()
runSimulation steps = do
Env npart ps <-; ask
forM_ [1..steps] (\i -> do
forM_ [1..npart+1] (\_ -> do
SimState config _ variables <-; get
selection <-; liftIO $ getRandomR (0, npart)
if selection < npart
then do
rdx <-; liftIO $ getRandomRs (0.0, _dx variables) >>= \r -> return $ take 3 r
rn <-; liftIO $ getRandomR (0, npart - 1)
lift $ moveParticle (vec3fromList rdx) rn
else do
let dv = _dv variables
rdv <-; liftIO $ getRandomR (-dv, dv)
racc <-; liftIO $ getRandomR (0.0, 1.0)
let vol = boxVolume $ _box config
acc = exp $ (-ps) * rdv + fromIntegral npart * log ((vol + rdv) / vol)
when (racc < acc) $ changeVolume rdv
)
when (i `mod` 100 == 0) $ do
state@(SimState _ _ variables) <-; get
let accVol = fromIntegral (_nVol variables) / 100.0
accMov = fromIntegral (_nMov variables) / (fromIntegral npart * 100)
olddx = _dx variables
olddv = _dv variables
newdx = (*) olddx (if accMov > 0.3 && olddx < 0.45 then 1.04 else 0.94)
newdv = (*) olddv (if accVol > 0.15 then 1.04 else 0.94)
liftIO $ print i
liftIO $ printf "dx = %.6f, dv = %.6f, nMov = %d, nVol = %d\n" newdx newdv (_nMov variables) (_nVol variables)
put $! state{_vars = variables{_dx = newdx, _dv = newdv, _nMov = 0, _nVol = 0}}
printConfig $ "Data/test" ++ show i ++ ".dat"
Here I'm effectively randomly choosing to either displace a particle or change the volume and every 100 steps I print useful information.
These two moves moveParticle, changeVolume are listed below,
moveParticle :: (RealFrac a, Floating a) => Vec3 a -> Int -> StateIO a ()
moveParticle dx nPart = do
SimState config@(Configuration box particles) cll variables <-; get
let particle = particles V.! nPart
particle' = Particle $ applyBC $ _position particle .+. dx
cllIdx = cellIndex (_nCells cll) (_position particle')
isCollision = any (checkCollision box) [(particle', particles V.! pId) | pId <-; neighbours cll cllIdx, pId /= nPart]
unless isCollision $ do --Accept move if there is no collision
let cllIdxOld = cellIndex (_nCells cll) (_position particle)
cll' = if cllIdx == cllIdxOld then cll else cllInsertAt (cllRemoveAt cll cllIdxOld nPart) cllIdx nPart
particles' = modifyVectorElement nPart particle' particles
put $! SimState config{_particles = particles'} cll' variables{_nMov = _nMov variables + 1}
changeVolume :: (Floating a, RealFrac a) => a -> Eval a ()
changeVolume dv = do
SimState config@(Configuration box particles) cll variables <-; get
Env npart _ <-; ask
let v = boxVolume box
box' = let scalef = ((v + dv) / v)**(1.0 / 3.0) in scaleBox scalef box
isCollision = any (checkCollision box') combinations
{-- A list of all the particles pairs that need to be checked for collision.
For each particle, the particles in the neighbouring cells are used. --}
combinations = do
pId <-; [0..(npart - 1)]
let particle = particles V.! pId
cllIdx = cellIndex (_nCells cll) (_position particle)
pId' <-; neighbours cll cllIdx
guard (pId /= pId')
return (particle, particles V.! pId')
unless isCollision $ do --Accept move if there is no collision
let new_cell_size = zipWith ((/) . fromIntegral) (_nCells cll) box'
new_n = map truncate box'
config' = config{_box = box'}
recreate = any (< 1.0) new_cell_size || any (> 2.0) new_cell_size || any (> 2) (zipWith ((abs .) . (-)) new_n (_nCells cll))
cll' = if recreate then generateCellList config' else cll
put $! SimState config' cll' variables{_nVol = _nVol variables + 1}
I would really appreciate a code review, especially since I just started learning Haskell. The code right now feels a bit messy. As you may have also noticed, I have given the accessor functions an underscore in their name with the intention to use Lenses later on.
I'm especially interested in any performance improvements. My Haskell implementation is up to 10 times slower than the C code which is unacceptable for a simulation. I should also mention that I've tried converting the Cell type to a Data.Seq and the Box type to a Vector but didn't see any performance improvements.
You can find the profiling report here and below I have also posted the space profiling graphs. I don't quite understand from these what is making the Haskell implementation so slow.
UPDATE: I have made most of the changes Peter suggested. The largest boost in performance came indeed from using the ST Monad to modify a Vector and from substituting the MonadRandom with the mersenne-random package.
The code I was using to modify Vectors is the following:
modifyVectorElement :: Int -> a -> V.Vector a -> V.Vector a
modifyVectorElement !i !x !v = V.modify (\v' -> write v' i x) v
, which I change to:
modifyVectorElement :: Int -> a -> Vector a -> Vector a
modifyVectorElement !i !x !v = runST $ unsafeThaw v >>= (\mv -> write mv i x >> unsafeFreeze mv)
Including all the other changes Peter suggested, brought the time for a 1000 iterations from ~30s to ~3.7s, which is about 25-35% slower than the C equivalent code. This is very satisfactory, but I'm sure more optimizations may be possible.
I'm still surprised by how big a difference optimization makes in Haskell compared to imperative languages, where beyond having a good algorithm, the performance gains from optimizations are much smaller.
For anyone that may be interested, I have updated the code on GitHub with the latest changes.
