diff --git a/accelerate-examples.cabal b/accelerate-examples.cabal --- a/accelerate-examples.cabal +++ b/accelerate-examples.cabal @@ -1,34 +1,29 @@ Name: accelerate-examples -Version: 0.12.1.0 +Version: 0.13.0.0 +Cabal-version: >=1.6 +Tested-with: GHC >= 7.6 +Build-type: Simple + +Synopsis: Examples using the Accelerate library +Description: + This package demonstrates a number of computation kernels and applications + showcasing the /Accelerate/ language and associated backend implementations. + It is also used for performance and regression testing. + . + Refer to the main /Accelerate/ package for more information: + + . + License: BSD3 License-file: LICENSE Author: The Accelerate Team Maintainer: Trevor L. McDonell +Homepage: https://github.com/AccelerateHS/accelerate-examples +Bug-reports: https://github.com/AccelerateHS/accelerate/issues + Category: Compilers/Interpreters, Concurrency, Data, Parallelism Stability: Experimental -Build-type: Simple -Cabal-version: >=1.6 -Homepage: http://www.cse.unsw.edu.au/~chak/project/accelerate/ -Synopsis: Examples using the Accelerate library -Description: This package defines a number of example programs for verification and - performance testing of the Accelerate language and associated backend - implementations. By default the package attempts to build with all available - backends, but this might not be desirable or even possible depending on your - available hardware. To disable a specific component, install with the extra - cabal flag(s): - . - . > cabal install accelerate-examples -f-opt - . - Where the available options are: - . - * cuda: A parallel CUDA backend for NVIDIA GPUs - . - * opencl: A parallel OpenCL backend - . - * io: Extra tests for reading and writing arrays in various formats - . - Flag gui Description: Enable gloss-based GUIs, where applicable. If not enabled, the application always runs in benchmark mode. @@ -40,11 +35,7 @@ Flag opencl Description: Enable the OpenCL parallel backend - Default: True - -Flag io - Description: Provide access to the block copy I/O functionality - Default: True + Default: False -- Run some basic randomised tests of the primitive library operations, and unit @@ -69,22 +60,23 @@ -fno-warn-orphans -fno-full-laziness -fno-excess-precision + if impl(ghc >= 7.0) ghc-options: -rtsopts if flag(cuda) CPP-options: -DACCELERATE_CUDA_BACKEND - build-depends: accelerate-cuda >= 0.12, - cuda >= 0.4 + build-depends: accelerate-cuda == 0.13.*, + cuda >= 0.5 && < 0.6 if flag(opencl) CPP-options: -DACCELERATE_OPENCL_BACKEND build-depends: accelerate-opencl == 0.1.* - build-depends: accelerate >= 0.12, - QuickCheck >= 2.0, - test-framework >= 0.5, - test-framework-quickcheck2 >= 0.2, + build-depends: accelerate == 0.13.*, + QuickCheck >= 2.0 && < 2.6, + test-framework >= 0.5 && < 0.9, + test-framework-quickcheck2 >= 0.2 && < 0.4, random @@ -99,11 +91,11 @@ examples/tests/primitives examples/tests/simple - other-modules: Benchmark, Config, Random, Test, Util, Validate, + other-modules: Benchmark, Config, PGM, Random, Test, Util, Validate, -- tests: -- image-processing - Canny, IntegralImage, PGM, + Canny, IntegralImage, -- io BlockCopy, VectorCopy, @@ -120,49 +112,51 @@ c-sources: examples/tests/io/fill_with_values.cpp extra-libraries: stdc++ - ghc-options: -Wall -O2 + ghc-options: -Wall -O2 -threaded -fno-excess-precision if impl(ghc >= 7.0) ghc-options: -rtsopts + if impl(ghc >= 7.6) + ghc-options: -fsimpl-tick-factor=130 + if flag(cuda) CPP-options: -DACCELERATE_CUDA_BACKEND - build-depends: accelerate-cuda >= 0.12 + build-depends: accelerate-cuda == 0.13.* if flag(opencl) CPP-options: -DACCELERATE_OPENCL_BACKEND build-depends: accelerate-opencl == 0.1.* - if flag(io) - CPP-options: -DACCELERATE_IO - build-depends: accelerate-io >= 0.12 - - build-depends: accelerate >= 0.12, - array >= 0.3, - attoparsec >= 0.10, + build-depends: accelerate == 0.13.*, + accelerate-io == 0.13.*, base == 4.*, - bytestring >= 0.9, - bytestring-lexing >= 0.2, - cmdargs >= 0.6, - criterion >= 0.5, - deepseq >= 1.1, - directory >= 1.0, - filepath >= 1.0, - hashtables >= 1.0.1, - mtl >= 1.1, - mwc-random >= 0.8, - pgm >= 0.1, - pretty >= 1.0, - vector >= 0.7, - vector-algorithms >= 0.4 + array >= 0.3 && < 0.5, + attoparsec >= 0.10 && < 0.11, + bytestring >= 0.9 && < 0.11, + bytestring-lexing >= 0.2 && < 0.5, + cmdargs >= 0.6 && < 0.11, + criterion >= 0.5 && < 0.9, + deepseq >= 1.1 && < 1.4, + directory >= 1.0 && < 1.3, + filepath >= 1.0 && < 1.4, + hashtables >= 1.0 && < 1.2, + mtl >= 1.1 && < 2.2, + mwc-random >= 0.8 && < 0.13, + pgm >= 0.1 && < 0.2, + pretty >= 1.0 && < 1.2, + repa >= 3.1 && < 3.3, + vector >= 0.7 && < 0.11, + vector-algorithms >= 0.4 && < 0.6 -- A quasicrystal demo as the sum of waves in a plane --- +-- -------------------------------------------------- + Executable accelerate-crystal - hs-source-dirs: examples/crystal + hs-source-dirs: examples/crystal lib Main-is: Main.hs - other-modules: Config - ghc-options: -O2 -Wall + other-modules: Config ParseArgs + ghc-options: -O2 -Wall -threaded if impl(ghc >= 7.0) ghc-options: -rtsopts @@ -171,25 +165,57 @@ if flag(cuda) CPP-options: -DACCELERATE_CUDA_BACKEND - build-depends: accelerate-cuda >= 0.12 + build-depends: accelerate-cuda == 0.13.* if flag(opencl) CPP-options: -DACCELERATE_OPENCL_BACKEND build-depends: accelerate-opencl == 0.1.* - build-depends: accelerate >= 0.12, + build-depends: accelerate == 0.13.*, base == 4.*, - criterion >= 0.5, - fclabels >= 1.0, - gloss >= 1.5 + criterion >= 0.5 && < 0.9, + fclabels >= 1.0 && < 1.2, + gloss >= 1.7 && < 1.8 + +-- Image edge detection +-- -------------------- + +Executable accelerate-canny + hs-source-dirs: examples/canny/src-acc lib + Main-is: Main.hs + other-modules: Config ParseArgs Canny Wildfire + ghc-options: -Wall -threaded -Odph -fllvm -optlo-O3 + + if impl(ghc >= 7.0) + ghc-options: -rtsopts + + if flag(cuda) + CPP-options: -DACCELERATE_CUDA_BACKEND + build-depends: accelerate-cuda == 0.13.* + + if flag(opencl) + CPP-options: -DACCELERATE_OPENCL_BACKEND + build-depends: accelerate-opencl == 0.1.* + + build-depends: accelerate == 0.13.*, + accelerate-io == 0.13.*, + base == 4.*, + criterion >= 0.5 && < 0.9, + fclabels >= 1.0 && < 1.2, + repa >= 3.1 && < 3.3, + repa-io >= 3.1 && < 3.3, + vector >= 0.7 && < 0.11 + + -- A simple mandelbrot generator +-- ----------------------------- Executable accelerate-mandelbrot - hs-source-dirs: examples/mandelbrot + hs-source-dirs: examples/mandelbrot lib Main-is: Main.hs - other-modules: Config - ghc-options: -O2 -Wall + other-modules: Config World Mandel ParseArgs + ghc-options: -O2 -Wall -threaded if impl(ghc >= 7.0) ghc-options: -rtsopts @@ -198,46 +224,153 @@ if flag(cuda) CPP-options: -DACCELERATE_CUDA_BACKEND - build-depends: accelerate-cuda >= 0.12 + build-depends: accelerate-cuda == 0.13.* if flag(opencl) CPP-options: -DACCELERATE_OPENCL_BACKEND build-depends: accelerate-opencl == 0.1.* - build-depends: accelerate >= 0.12, + build-depends: accelerate == 0.13.*, base == 4.*, - criterion >= 0.5, - fclabels >= 1.0, - gloss >= 1.5 + criterion >= 0.5 && < 0.9, + fclabels >= 1.0 && < 1.2, + gloss >= 1.7 && < 1.8 -- A stable fluid simulation --- +-- ------------------------- + Executable accelerate-fluid Main-is: Main.hs - hs-source-dirs: examples/fluid/src + hs-source-dirs: examples/fluid/src-acc lib other-modules: Config + ParseArgs Event Fluid Type World - ghc-options: -O2 -Wall + ghc-options: -O2 -Wall -threaded + if impl(ghc >= 7.0) + ghc-options: -rtsopts if flag(gui) CPP-options: -DACCELERATE_ENABLE_GUI if flag(cuda) cpp-options: -DACCELERATE_CUDA_BACKEND - build-depends: accelerate-cuda >= 0.12 + build-depends: accelerate-cuda == 0.13.* - Build-depends: accelerate >= 0.12, - accelerate-io >= 0.12, + Build-depends: accelerate == 0.13.*, + accelerate-io == 0.13.*, base == 4.*, - bmp >= 1.2, - criterion >= 0.5, - fclabels >= 1.0, - gloss >= 1.7 + bmp >= 1.2 && < 1.3, + criterion >= 0.5 && < 0.9, + fclabels >= 1.0 && < 1.2, + gloss >= 1.7 && < 1.8 + + +-- Simulation of gravitational attraction between solid particles +-- -------------------------------------------------------------- + +Executable accelerate-nbody + Main-is: Main.hs + hs-source-dirs: examples/n-body lib + other-modules: Config + ParseArgs + Common.Body + Common.Dump + Common.Tree + Common.Type + Common.Util + Common.World + Gloss.Draw + Gloss.Event + Gloss.Simulate + Random.Array + Random.Position + Solver.BarnsHut + Solver.Naive + + ghc-options: -O2 -Wall -threaded + if impl(ghc >= 7.0) + ghc-options: -rtsopts + + if flag(gui) + cpp-options: -DACCELERATE_ENABLE_GUI + + if flag(cuda) + cpp-options: -DACCELERATE_CUDA_BACKEND + build-depends: accelerate-cuda == 0.13.* + + Build-depends: accelerate == 0.13.*, + base == 4.*, + criterion >= 0.5 && < 0.9, + fclabels >= 1.0 && < 1.2, + gloss >= 1.7 && < 1.8 + + +-- A celular automata +-- ------------------ + +Executable accelerate-smoothlife + Main-is: Main.hs + hs-source-dirs: examples/smoothlife lib + other-modules: Config + ParseArgs + SmoothLife + Gloss.Draw + Gloss.Event + Gloss.Simulate + Random.Array + Random.Splat + + ghc-options: -O2 -Wall -threaded + if impl(ghc >= 7.0) + ghc-options: -rtsopts + + if flag(gui) + cpp-options: -DACCELERATE_ENABLE_GUI + + if flag(cuda) + cpp-options: -DACCELERATE_CUDA_BACKEND + build-depends: accelerate-cuda == 0.13.* + + Build-depends: accelerate == 0.13.*, + accelerate-fft == 0.13.*, + base == 4.*, + criterion >= 0.5 && < 0.9, + fclabels >= 1.0 && < 1.2, + gloss >= 1.7 && < 1.8 + + +-- A password recovery tool +-- ------------------------ + +Executable accelerate-hashcat + Main-is: Main.hs + hs-source-dirs: examples/hashcat lib + other-modules: Config + ParseArgs + Digest + MD5 + Util + + ghc-options: -O2 -Wall -threaded + if impl(ghc >= 7.0) + ghc-options: -rtsopts + + if flag(cuda) + cpp-options: -DACCELERATE_CUDA_BACKEND + build-depends: accelerate-cuda == 0.13.* + + Build-depends: accelerate == 0.13.*, + base == 4.*, + bytestring >= 0.9 && < 0.11, + bytestring-lexing >= 0.2 && < 0.5, + cereal >= 0.3 && < 0.4, + criterion >= 0.5 && < 0.9, + fclabels >= 1.0 && < 1.2 source-repository head diff --git a/examples/canny/src-acc/Canny.hs b/examples/canny/src-acc/Canny.hs new file mode 100644 --- /dev/null +++ b/examples/canny/src-acc/Canny.hs @@ -0,0 +1,209 @@ +{-# LANGUAGE ViewPatterns #-} +-- An implementation of the Canny edge detection algorithm +-- +-- J. F. Canny, "A Computational Approach to Edge Detection" in _Pattern +-- Analysis and Machine Intelligence_, (6), 1986. +-- +-- This module defines the first seven data-parallel sections of the algorithm. +-- The last phase, which uses a recursive algorithm to "connect" pixels that +-- form the output image, is defined in the module Wildfire and implemented +-- using Repa. +-- + +module Canny where + +import Prelude as P + +import Data.Array.Accelerate as A +import Data.Array.Accelerate.IO as A + + +-- Canny algorithm ------------------------------------------------------------- + +canny :: Float -> Float -> Acc (Image RGBA) -> (Acc (Image Float), Acc (Vector Int)) +canny (constant -> low) (constant -> high) + = stage1 + . nonMaximumSuppression low high + . gradientMagDir low + . gaussianY + . gaussianX + . toGreyscale + where + stage1 x = (x, selectStrong x) + + +-- Accelerate component -------------------------------------------------------- + +type RGBA = Word32 +type Image a = Array DIM2 a + +type Stencil5x1 a = (Stencil3 a, Stencil5 a, Stencil3 a) +type Stencil1x5 a = (Stencil3 a, Stencil3 a, Stencil3 a, Stencil3 a, Stencil3 a) + +-- Classification of the output pixel +data Orient = Undef | PosD | Vert | NegD | Horiz +data Edge = None | Weak | Strong + +orient :: Orient -> Int +orient Undef = 0 +orient PosD = 64 +orient Vert = 128 +orient NegD = 192 +orient Horiz = 255 + +orient' :: Orient -> Exp Int +orient' = constant . orient + +edge :: Edge -> Float +edge None = 0 +edge Weak = 0.5 +edge Strong = 1.0 + +edge' :: Edge -> Exp Float +edge' = constant . edge + +convolve5x1 :: (Elt a, IsNum a) => [Exp a] -> Stencil5x1 a -> Exp a +convolve5x1 kernel (_, (a,b,c,d,e), _) + = P.sum $ P.zipWith (*) kernel [a,b,c,d,e] + +convolve1x5 :: (Elt a, IsNum a) => [Exp a] -> Stencil1x5 a -> Exp a +convolve1x5 kernel ((_,a,_), (_,b,_), (_,c,_), (_,d,_), (_,e,_)) + = P.sum $ P.zipWith (*) kernel [a,b,c,d,e] + + +-- RGB to Greyscale conversion, in the range [0,255] +-- +toGreyscale :: Acc (Image RGBA) -> Acc (Image Float) +toGreyscale = A.map (\rgba -> 255 * luminanceOfRGBA32 rgba) + + +-- Separable Gaussian blur in the x- and y-directions +-- +gaussianX :: Acc (Image Float) -> Acc (Image Float) +gaussianX = stencil (convolve5x1 gaussian) Clamp + where + gaussian = [ 1, 4, 6, 4, 1 ] + +gaussianY :: Acc (Image Float) -> Acc (Image Float) +gaussianY = stencil (convolve1x5 gaussian) Clamp + where + gaussian = P.map (/256) [ 1, 4, 6, 4, 1 ] + + +-- Gradients in the x- and y- directions +-- +gradientX :: Acc (Image Float) -> Acc (Image Float) +gradientX = stencil grad Clamp + where + grad :: Stencil3x3 Float -> Exp Float + grad ((u, _, x) + ,(v, _, y) + ,(w, _, z)) = x + (2*y) + z - u - (2*v) - w + +gradientY :: Acc (Image Float) -> Acc (Image Float) +gradientY = stencil grad Clamp + where + grad :: Stencil3x3 Float -> Exp Float + grad ((x, y, z) + ,(_, _, _) + ,(u, v, w)) = x + (2*y) + z - u - (2*v) - w + + +-- Classify the magnitude and orientation of the image gradient. +-- +-- Because accelerate supports generalised stencil functions, not just +-- convolutions, we can combine the x- and y- sobel operators and save some +-- memory bandwidth. +-- +gradientMagDir + :: Exp Float + -> Acc (Image Float) + -> Acc (Array DIM2 (Float,Int)) +gradientMagDir low = stencil magdir Clamp + where + magdir :: Stencil3x3 Float -> Exp (Float,Int) + magdir ((v0, v1, v2) + ,(v3, _, v4) + ,(v5, v6, v7)) = + let + -- Image gradients + dx = v2 + (2*v4) + v7 - v0 - (2*v3) - v5 + dy = v0 + (2*v1) + v2 - v5 - (2*v6) - v7 + + -- Magnitude + mag = sqrt (dx * dx + dy + dy) + + -- Direction + -- + -- Determine the angle of the vector and rotate it around a bit to + -- make the segments easier to classify + theta = atan2 dy dx + alpha = (theta - (pi/8)) * (4/pi) + + -- Normalise the angle to between [0..8) + norm = alpha + 8 * A.fromIntegral (boolToInt (alpha <=* 0)) + + -- Try to avoid doing explicit tests, to avoid warp divergence + undef = abs dx <=* low &&* abs dy <=* low + dir = boolToInt (A.not undef) * ((64 * (1 + A.floor norm `mod` 4)) `A.min` 255) + in + lift (mag, dir) + + +-- Non-maximum suppression classifies pixels that are the local maximum along +-- the direction of the image gradient as either strong or weak edges. All other +-- pixels are not considered edges at all. +-- +-- The image intensity is in the range [0,1] +-- +nonMaximumSuppression + :: Exp Float + -> Exp Float + -> Acc (Image (Float,Int)) + -> Acc (Image Float) +nonMaximumSuppression low high magdir = + generate (shape magdir) $ \ix -> + let -- The input parameters + -- + (mag, dir) = unlift (magdir ! ix) + Z :. h :. w = unlift (shape magdir) + Z :. y :. x = unlift ix + + -- Determine the points that lie either side of this point along to the + -- direction of the image gradient. + -- + -- The direction coding: + -- + -- 192 128 64 + -- | + -- 255 --- --- + -- + offsetx = dir >* orient' Vert ? (-1, dir <* orient' Vert ? (1, 0)) + offsety = dir <* orient' Horiz ? (-1, 0) + + (fwd, _) = unlift $ magdir ! lift (clamp (Z :. y+offsety :. x+offsetx)) :: (Exp Float, Exp Int) + (rev, _) = unlift $ magdir ! lift (clamp (Z :. y-offsety :. x-offsetx)) :: (Exp Float, Exp Int) + + clamp (Z:.u:.v) = Z :. 0 `A.max` u `A.min` (h-1) :. 0 `A.max` v `A.min` (w-1) + + -- Try to avoid doing explicit tests to avoid warp divergence. + -- + none = dir ==* orient' Undef ||* mag <* low ||* mag <* fwd ||* mag <* rev + strong = mag >=* high + in + A.fromIntegral (boolToInt (A.not none) * (1 + boolToInt strong)) * 0.5 + + +-- Extract the linear indices of the strong edges +-- +selectStrong + :: Acc (Image Float) + -> Acc (Array DIM1 Int) +selectStrong img = + let strong = A.map (\x -> boolToInt (x ==* edge' Strong)) (flatten img) + (targetIdx, len) = A.scanl' (+) 0 strong + indices = A.enumFromN (index1 $ size img) 0 + zeros = A.fill (index1 $ the len) 0 + in + A.permute const zeros (\ix -> strong!ix ==* 0 ? (ignore, index1 $ targetIdx!ix)) indices + diff --git a/examples/canny/src-acc/Config.hs b/examples/canny/src-acc/Config.hs new file mode 100644 --- /dev/null +++ b/examples/canny/src-acc/Config.hs @@ -0,0 +1,72 @@ +{-# LANGUAGE TemplateHaskell #-} + +module Config where + +import ParseArgs +import Data.Label + +data Config + = Config + { + -- Standard options + _configBackend :: Backend + , _configHelp :: Bool + , _configBenchmark :: Bool + + -- How to execute the program + , _configThreshLow :: Float + , _configThreshHigh :: Float + + } + +$(mkLabels [''Config]) + +defaults :: Config +defaults = Config + { + _configBackend = maxBound + , _configHelp = False + , _configBenchmark = False + + , _configThreshLow = 50 + , _configThreshHigh = 100 + } + + +options :: [OptDescr (Config -> Config)] +options = + [ Option ['w'] ["threshold-low"] + (ReqArg (set configThreshLow . read) "FLOAT") + (describe configThreshLow "threshold value for weak edges") + + , Option ['s'] ["threshold-high"] + (ReqArg (set configThreshHigh . read) "FLOAT") + (describe configThreshHigh "threshold value for strong edges") + + , Option [] ["benchmark"] + (NoArg (set configBenchmark True)) + (describe configBenchmark "benchmark instead of displaying animation") + + , Option ['h', '?'] ["help"] + (NoArg (set configHelp True)) + "show this help message" + ] + where + describe f msg + = msg ++ " (" ++ show (get f defaults) ++ ")" + + +-- | Process the command line options +-- + +header :: [String] +header = + [ "accelerate-canny (c) [2007..2013] The Accelerate Team" + , "" + , "Usage: accelerate-canny [OPTIONS] fileIn.bmp fileOut.bmp" + , "" + ] + +footer :: [String] +footer = [] + diff --git a/examples/canny/src-acc/Main.hs b/examples/canny/src-acc/Main.hs new file mode 100644 --- /dev/null +++ b/examples/canny/src-acc/Main.hs @@ -0,0 +1,80 @@ +{-# LANGUAGE TypeOperators #-} + +import Canny +import Config +import Wildfire +import ParseArgs + +import Prelude as P +import Data.Label +import Criterion.Main ( defaultMainWith, bgroup, bench, whnf ) +import System.Exit +import System.Environment + +import Data.Array.Accelerate as A +import qualified Data.Array.Accelerate.IO as A +import qualified Data.Array.Repa.IO.BMP as R +import qualified Data.Array.Repa.Repr.Unboxed as R + + +-- Main routine ---------------------------------------------------------------- + +main :: IO () +main + = do argv <- getArgs + (conf, cconf, nops) <- parseArgs configHelp configBackend options defaults header footer argv + (fileIn, fileOut) <- case nops of + (i:o:_) -> return (i,o) + _ -> parseArgs configHelp configBackend options defaults header footer ("--help":argv) + >> exitSuccess + + -- Read in the image file + img <- either (error . show) id `fmap` A.readImageFromBMP fileIn + + -- Set up the algorithm parameters + let threshLow = get configThreshLow conf + threshHigh = get configThreshHigh conf + backend = get configBackend conf + + -- Set up the partial results so that we can benchmark individual + -- kernel stages. + low = constant threshLow + high = constant threshHigh + + grey' = run backend $ toGreyscale (use img) + blurX' = run backend $ gaussianX (use grey') + blurred' = run backend $ gaussianY (use blurX') + magdir' = run backend $ gradientMagDir low (use blurred') + suppress' = run backend $ nonMaximumSuppression low high (use magdir') + + if P.not (get configBenchmark conf) + then do + -- Connect the strong and weak edges of the image using Repa, and + -- write the final image to file + -- + let (image, strong) = run backend $ A.lift (canny threshLow threshHigh (use img)) + edges <- wildfire (A.toRepa image) (A.toRepa strong) + R.writeImageToBMP fileOut (R.zip3 edges edges edges) + + else do + -- Run each of the individual kernel stages through criterion, as + -- well as the end-to-end step process. + -- + withArgs (P.drop 2 nops) $ defaultMainWith cconf (return ()) + [ bgroup "kernels" + [ bench "greyscale" $ whnf (run1 backend toGreyscale) img + , bench "blur-x" $ whnf (run1 backend gaussianX) grey' + , bench "blur-y" $ whnf (run1 backend gaussianY) blurX' + , bench "grad-x" $ whnf (run1 backend gradientX) blurred' + , bench "grad-y" $ whnf (run1 backend gradientY) blurred' + , bench "mag-orient" $ whnf (run1 backend (gradientMagDir low)) blurred' + , bench "suppress" $ whnf (run1 backend (nonMaximumSuppression low high)) magdir' + , bench "select" $ whnf (run1 backend selectStrong) suppress' + ] + + , bgroup "canny" + [ bench "run" $ whnf (run backend . (P.snd . canny threshLow threshHigh)) (use img) + , bench "run1" $ whnf (run1 backend (P.snd . canny threshLow threshHigh)) img + ] + ] + diff --git a/examples/canny/src-acc/Wildfire.hs b/examples/canny/src-acc/Wildfire.hs new file mode 100644 --- /dev/null +++ b/examples/canny/src-acc/Wildfire.hs @@ -0,0 +1,102 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE TypeOperators #-} +-- This module defines the final phase of the Canny algorithm, a recursive +-- algorithm to "connect" pixels of the output image. +-- +-- This function originally written by Ben Lippmeier for Repa. +-- +-- NOTE: for best performance this needs to be compiled with the following GHC options: +-- -fllvm -optlo-O3 -Odph -fno-liberate-case +-- -funfolding-use-threshold100 -funfolding-keeness-factor100 +-- + +module Wildfire where + +import Canny + +import Data.Word +import Data.Array.Accelerate.IO ( A ) +import Data.Array.Repa.Repr.Unboxed ( U ) +import qualified Data.Array.Repa as R +import qualified Data.Vector.Unboxed as V +import qualified Data.Vector.Unboxed.Mutable as VM + + +-- Repa component -------------------------------------------------------------- + +-- | Trace out strong edges in the final image. +-- Also trace out weak edges that are connected to strong edges. +-- +wildfire + :: R.Array A R.DIM2 Float -- ^ Image with strong and weak edges set. + -> R.Array A R.DIM1 Int -- ^ Array containing flat indices of strong edges. + -> IO (R.Array U R.DIM2 Word8) +{-# NOINLINE wildfire #-} +wildfire img arrStrong + = do (sh, vec) <- wildfireIO + return $ sh `seq` vec `seq` R.fromUnboxed sh vec + + where lenImg = R.size $ R.extent img + lenStrong = R.size $ R.extent arrStrong + shImg = R.extent img + + wildfireIO + = do -- Stack of image indices we still need to consider. + vStrong <- R.toUnboxed `fmap` R.computeUnboxedP (R.delay arrStrong) + vStrong' <- V.thaw vStrong + vStack <- VM.grow vStrong' (lenImg - lenStrong) + + -- Burn in new edges. + vImg <- VM.unsafeNew lenImg + VM.set vImg 0 + burn vImg vStack lenStrong + vImg' <- V.unsafeFreeze vImg + return (R.extent img, vImg') + + + burn :: VM.IOVector Word8 -> VM.IOVector Int -> Int -> IO () + burn !vImg !vStack !top + | top == 0 + = return () + + | otherwise + = do let !top' = top - 1 + n <- VM.unsafeRead vStack top' + let (R.Z R.:. y R.:. x) = R.fromIndex (R.extent img) n + + let {-# INLINE push #-} + push ix t = + if R.inShape shImg ix + then pushWeak vImg vStack ix t + else return t + + VM.write vImg n 255 + >> push (R.Z R.:. y - 1 R.:. x - 1) top' + >>= push (R.Z R.:. y - 1 R.:. x ) + >>= push (R.Z R.:. y - 1 R.:. x + 1) + + >>= push (R.Z R.:. y R.:. x - 1) + >>= push (R.Z R.:. y R.:. x + 1) + + >>= push (R.Z R.:. y + 1 R.:. x - 1) + >>= push (R.Z R.:. y + 1 R.:. x ) + >>= push (R.Z R.:. y + 1 R.:. x + 1) + + >>= burn vImg vStack + + -- If this ix is weak in the source then set it to strong in the + -- result and push the ix onto the stack. + {-# INLINE pushWeak #-} + pushWeak vImg vStack ix top + = do let n = R.toIndex (R.extent img) ix + xDst <- VM.unsafeRead vImg n + let xSrc = img `R.unsafeIndex` ix + + if xDst == 0 + && xSrc == edge Weak + then do + VM.unsafeWrite vStack top (R.toIndex (R.extent img) ix) + return (top + 1) + + else return top + diff --git a/examples/crystal/Config.hs b/examples/crystal/Config.hs --- a/examples/crystal/Config.hs +++ b/examples/crystal/Config.hs @@ -1,27 +1,11 @@ -{-# LANGUAGE CPP, TemplateHaskell, PatternGuards #-} - -module Config ( - - Options, optBackend, optSize, optZoom, optScale, optDegree, optBench, - processArgs, run +{-# LANGUAGE CPP #-} +{-# LANGUAGE TemplateHaskell #-} -) where +module Config where +import ParseArgs import Data.Label -import System.Exit -import System.Console.GetOpt -import Data.Array.Accelerate ( Arrays, Acc ) -import qualified Data.Array.Accelerate.Interpreter as Interp -#ifdef ACCELERATE_CUDA_BACKEND -import qualified Data.Array.Accelerate.CUDA as CUDA -#endif -data Backend = Interpreter -#ifdef ACCELERATE_CUDA_BACKEND - | CUDA -#endif - deriving (Bounded, Show) - data Options = Options { _optBackend :: Backend @@ -36,8 +20,8 @@ $(mkLabels [''Options]) -defaultOptions :: Options -defaultOptions = Options +defaults :: Options +defaults = Options { _optBackend = maxBound , _optSize = 200 , _optZoom = 3 @@ -52,21 +36,9 @@ } -run :: (Arrays a, Arrays b) => Options -> (Acc a -> Acc b) -> a -> b -run opts f = case _optBackend opts of - Interpreter -> head . Interp.stream f . return -#ifdef ACCELERATE_CUDA_BACKEND - CUDA -> CUDA.run1 f -#endif - - options :: [OptDescr (Options -> Options)] options = - [ Option [] ["interpreter"] (NoArg (set optBackend Interpreter)) "reference implementation (sequential)" -#ifdef ACCELERATE_CUDA_BACKEND - , Option [] ["cuda"] (NoArg (set optBackend CUDA)) "implementation for NVIDIA GPUs (parallel)" -#endif - , Option [] ["size"] (ReqArg (set optSize . read) "INT") "visualisation size (200)" + [ Option [] ["size"] (ReqArg (set optSize . read) "INT") "visualisation size (200)" , Option [] ["zoom"] (ReqArg (set optZoom . read) "INT") "pixel replication factor (3)" , Option [] ["scale"] (ReqArg (set optScale . read) "FLOAT") "feature size of visualisation (30)" , Option [] ["degree"] (ReqArg (set optDegree . read) "INT") "number of waves to sum for each point (5)" @@ -75,19 +47,14 @@ ] -processArgs :: [String] -> IO (Options, [String]) -processArgs argv = - case getOpt' Permute options argv of - (o,_,n,[]) -> case foldl (flip id) defaultOptions o of - opts | False <- get optHelp opts -> return (opts, n) - opts | True <- get optBench opts -> return (opts, "--help":n) - _ -> putStrLn (helpMsg []) >> exitSuccess - (_,_,_,err) -> error (helpMsg err) - where - helpMsg err = concat err ++ usageInfo header options - header = unlines - [ "accelerate-crystal (c) [2011..2012] The Accelerate Team" - , "" - , "Usage: accelerate-crystal [OPTIONS]" - ] +header :: [String] +header = + [ "accelerate-crystal (c) [2011..2013] The Accelerate Team" + , "" + , "Usage: accelerate-crystal [OPTIONS]" + , "" + ] + +footer :: [String] +footer = [] diff --git a/examples/crystal/Main.hs b/examples/crystal/Main.hs --- a/examples/crystal/Main.hs +++ b/examples/crystal/Main.hs @@ -9,13 +9,14 @@ module Main where import Config +import ParseArgs import Data.Word import Data.Label import Foreign.Ptr import Control.Monad import Control.Exception -import Criterion.Main ( defaultMain, bench, whnf ) +import Criterion.Main ( defaultMainWith, bench, whnf ) import Foreign.ForeignPtr import System.Environment import System.IO.Unsafe @@ -130,9 +131,9 @@ -- Rendering ------------------------------------------------------------------ --- | Compute a single frame of the animation as a Gloss picture. ---frame :: Size -> Scale -> Zoom -> Degree -> Float -> G.Picture +-- | Compute a single frame of the animation as a Gloss picture. +-- frame :: Render -> Size -> Zoom -> Float -> G.Picture frame render size zoom time = G.scale zoom' zoom' pic where @@ -167,33 +168,27 @@ -- Main ----------------------------------------------------------------------- main :: IO () main - = do (config, nops) <- processArgs =<< getArgs + = do argv <- getArgs + (config, crit, nops) <- parseArgs optHelp optBackend options defaults header footer argv let size = get optSize config zoom = get optZoom config scale = get optScale config degree = get optDegree config - render = run config $ makeImage size scale degree + backend = get optBackend config + render = run1 backend $ makeImage size scale degree force arr = A.indexArray arr (Z:.0:.0) `seq` arr void . evaluate $ render (A.fromList Z [0]) if get optBench config - then withArgs nops $ defaultMain + then withArgs nops $ defaultMainWith crit (return ()) [ bench "crystal" $ whnf (force . render) (A.fromList Z [1.0]) ] #ifndef ACCELERATE_ENABLE_GUI else return () - -#elif MIN_VERSION_gloss(1,6,0) +#else else G.animate (G.InWindow "Quasicrystals" (size * zoom, size * zoom) (10, 10)) - (G.black) - (frame render size zoom) -#else - else G.animateInWindow - "Quasicrystals" - (size * zoom, size * zoom) - (10, 10) (G.black) (frame render size zoom) #endif diff --git a/examples/fluid/src-acc/Config.hs b/examples/fluid/src-acc/Config.hs new file mode 100644 --- /dev/null +++ b/examples/fluid/src-acc/Config.hs @@ -0,0 +1,309 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE PatternGuards #-} +{-# LANGUAGE TemplateHaskell #-} +-- +-- Configuration parameters +-- + +module Config where + +import Type +import ParseArgs + +import Data.Label +import Control.Monad +import Prelude as P +import Data.Array.Accelerate as A +import Data.Array.Accelerate.IO as A + + +data Initial a + = FromFile FilePath + | FromFunction (Backend -> Int -> Int -> a) + + +data Config = Config + { + -- simulation + _viscosity :: !Float + , _diffusion :: !Float + , _timestep :: !Float + , _inputDensity :: !Float + , _inputVelocity :: !Float + , _simulationSteps :: !Int + , _simulationWidth :: !Int + , _simulationHeight :: !Int + + -- visualisation + , _displayScale :: !Int + , _displayFramerate :: !Int + , _initialDensity :: DensityField + , _initialVelocity :: VelocityField + + -- misc + , _optBackend :: !Backend + , _optBench :: !Bool + , _optHelp :: !Bool + + -- extra options to specify initial conditions for command parsing + , _setupDensity :: Initial DensityField + , _setupVelocity :: Initial VelocityField + } + +$(mkLabels [''Config]) + + +defaults :: Config +defaults = Config + { _viscosity = 0 + , _diffusion = 0 + , _timestep = 0.1 + , _inputDensity = 100 + , _inputVelocity = 20 + , _simulationSteps = 40 + , _simulationWidth = 100 + , _simulationHeight = 100 + + , _displayScale = 5 + , _displayFramerate = 25 + , _initialDensity = error "initial density??" + , _initialVelocity = error "initial velocity??" + + , _optBackend = maxBound + , _optBench = False + , _optHelp = False + + , _setupDensity = FromFunction makeField_empty + , _setupVelocity = FromFunction makeField_empty + } + + +-- | The set of available command-line options +-- +options :: [OptDescr (Config -> Config)] +options = + -- Simulation options + [ Option [] ["viscosity"] + (ReqArg (parse viscosity) "FLOAT") + (describe viscosity "viscosity for velocity damping") + + , Option [] ["diffusion"] + (ReqArg (parse diffusion) "FLOAT") + (describe diffusion "diffusion rate for mass dispersion") + + , Option [] ["delta"] + (ReqArg (parse timestep) "FLOAT") + (describe timestep "simulation time between each frame") + + , Option [] ["density"] + (ReqArg (parse inputDensity) "FLOAT") + (describe inputDensity "magnitude of user input density") + + , Option [] ["velocity"] + (ReqArg (parse inputVelocity) "FLOAT") + (describe inputVelocity "magnitude of user input velocity") + + , Option [] ["iterations"] + (ReqArg (parse simulationSteps) "INT") + (describe simulationSteps "number of iterations of the linear solver") + + , Option [] ["width"] + (ReqArg (parse simulationWidth) "INT") + (describe simulationWidth "grid width of simulation") + + , Option [] ["height"] + (ReqArg (parse simulationHeight) "INT") + (describe simulationHeight "grid height of simulation") + + -- Display options + , Option [] ["scale"] + (ReqArg (parse displayScale) "INT") + (describe displayScale "feature size of visualisation") + + , Option [] ["framerate"] + (ReqArg (parse displayFramerate) "INT") + (describe displayFramerate "frame rate for visualisation") + + -- Initial conditions + , Option [] ["bmp-density"] + (ReqArg (set setupDensity . FromFile) "FILE.bmp") + "file for initial fluid density" + + , Option [] ["bmp-velocity"] + (ReqArg (set setupVelocity . FromFile) "FILE.bmp") + "file for initial fluid velocity" + + , Option [] ["init-checks"] + (NoArg init_checks) + "initial density field with zero velocity field" + + , Option [] ["init-man"] + (NoArg init_man) + "initial density field with swirling velocity" + + , Option [] ["init-elk"] + (NoArg init_elk) + "initial density field with swirling velocity" + + -- Miscellaneous + , Option [] ["benchmark"] + (NoArg (set optBench True)) + (describe optBench "benchmark instead of displaying animation") + + , Option ['h','?'] ["help"] + (NoArg (set optHelp True)) + (describe optHelp "show help message") + ] + where + parse f x = set f (read x) + describe f msg = msg ++ " (" ++ show (get f defaults) ++ ")" + + init_checks = set setupDensity (FromFunction makeDensity_checks) + . set setupVelocity (FromFunction makeField_empty) + + init_man = set setupDensity (FromFunction makeDensity_checks) + . set setupVelocity (FromFunction makeVelocity_man) + + init_elk = set setupDensity (FromFunction makeDensity_checks) + . set setupVelocity (FromFunction makeVelocity_elk) + + +header :: [String] +header = + [ "accelerate-fluid (c) [2011..2013] The Accelerate Team" + , "" + , "Usage: accelerate-fluid [OPTIONS]" + , "" + ] + +footer :: [String] +footer = + [ "" + , "Runtime usage:" + , " click add density sources to the image" + , " shift-click add velocity sources" + , " r reset the image" + , " d toggle display of density field" + , " v toggle display of velocity field lines" + ] + + +-- Initial conditions +-- ------------------ + +initialiseConfig :: Config -> IO Config +initialiseConfig conf = do + let backend = get optBackend conf + width = get simulationWidth conf + height = get simulationHeight conf + + dens <- case get setupDensity conf of + FromFile fn -> loadDensity_bmp backend fn width height + FromFunction f -> return (f backend width height) + + velo <- case get setupVelocity conf of + FromFile fn -> loadVelocity_bmp backend fn width height + FromFunction f -> return (f backend width height) + + return . set initialDensity dens + . set initialVelocity velo + $ conf + + +makeField_empty :: FieldElt e => Backend -> Int -> Int -> Field e +makeField_empty backend width height + = run backend + $ A.fill (constant (Z:.height:.width)) (constant zero) + + +makeDensity_checks :: Backend -> Int -> Int -> DensityField +makeDensity_checks backend width height + = let width' = constant $ P.fromIntegral width + height' = constant $ P.fromIntegral height + yc = constant $ P.fromIntegral (height `div` 2) + xc = constant $ P.fromIntegral (width `div` 2) + + checks ix + = let Z :. y :. x = unlift ix + x' = A.fromIntegral x + y' = A.fromIntegral y + tx = 10 * (x' - xc) / width' + ty = 10 * (y' - yc) / height' + xk1 = abs tx >* 3*pi/2 ? (0 , cos tx) + yk1 = abs ty >* 3*pi/2 ? (0 , cos ty) + d1 = xk1 * yk1 + in + 0 `A.max` d1 + in + run backend $ A.generate (constant (Z:.height:.width)) checks + + +makeVelocity_man :: Backend -> Int -> Int -> VelocityField +makeVelocity_man backend width height + = let width' = constant $ P.fromIntegral width + height' = constant $ P.fromIntegral height + yc = constant $ P.fromIntegral (height `div` 2) + xc = constant $ P.fromIntegral (width `div` 2) + + man ix + = let Z :. y :. x = unlift ix + x' = A.fromIntegral x + y' = A.fromIntegral y + xk2 = cos (19 * (x' - xc) / width') + yk2 = cos (17 * (y' - yc) / height') + d2 = xk2 * yk2 / 5 + in + lift (constant 0, d2) + in + run backend $ A.generate (constant (Z:.height:.width)) man + + +makeVelocity_elk :: Backend -> Int -> Int -> VelocityField +makeVelocity_elk backend width height + = let width' = constant $ P.fromIntegral width + height' = constant $ P.fromIntegral height + yc = constant $ P.fromIntegral (height `div` 2) + xc = constant $ P.fromIntegral (width `div` 2) + + elk ix + = let Z :. y :. x = unlift ix + x' = A.fromIntegral x + y' = A.fromIntegral y + tx = 12 * (x' - xc) / width' + ty = 12 * (y' - yc) / height' + xk2 = cos tx + yk2 = -cos ty + d2 = xk2 * yk2 / 5 + in + lift (constant 0, d2) + in + run backend $ A.generate (constant (Z:.height:.width)) elk + + +loadDensity_bmp :: Backend -> FilePath -> Int -> Int -> IO DensityField +loadDensity_bmp backend filepath width height + = do arr <- either (error . show) id `fmap` readImageFromBMP filepath + let Z:.h:.w = arrayShape arr + + when (w /= width || h /= height) + $ error "accelerate-fluid: density-bmp does not match width x height" + + return . run backend $ A.map luminanceOfRGBA32 (use arr) + + +loadVelocity_bmp :: Backend -> FilePath -> Int -> Int -> IO VelocityField +loadVelocity_bmp backend filepath width height + = do arr <- either (error . show) id `fmap` readImageFromBMP filepath + let Z:.h:.w = arrayShape arr + + when (w /= width || h /= height) + $ error "accelerate-fluid: velocity-bmp does not match width x height" + + let conv rgb = + let (r,g,_,_) = unlift (unpackRGBA32 rgb) :: (Exp Word8, Exp Word8, Exp Word8, Exp Word8) + r' = A.fromIntegral (-128 + A.fromIntegral r :: Exp Int) + g' = A.fromIntegral (-128 + A.fromIntegral g :: Exp Int) + in lift (r' * 0.0001, g' * 0.0001) + + return . run backend $ A.map conv (use arr) + diff --git a/examples/fluid/src-acc/Event.hs b/examples/fluid/src-acc/Event.hs new file mode 100644 --- /dev/null +++ b/examples/fluid/src-acc/Event.hs @@ -0,0 +1,83 @@ +{-# LANGUAGE PatternGuards #-} +-- +-- Event handling +-- + +module Event where + +import Config +import World +import Data.Label +import Graphics.Gloss.Interface.Pure.Game +import Data.Array.Accelerate ( Z(..), (:.)(..) ) + + +-- Event locations are returned as window coordinates, where the origin is in +-- the centre of the window and increases to the right and up. If the simulation +-- size is (100,100) with scale factor of 4, then the event coordinates are +-- returned in the range [-200,200]. +-- +react :: Config -> Event -> World -> World +react opt event world = + case event of + EventKey (Char c) s m _ -> keyboard c s m + EventKey (MouseButton LeftButton) s m uv -> mouse m s (coord uv) + EventMotion uv -> motion (coord uv) + _ -> world + where + -- Inject a new density source when the left button is clicked. + -- + -- If the shift key is held, remember the location and add a new velocity + -- source between the old and new points as the mouse moves. + -- + mouse key button xy + | Up <- shift key, Down <- button + = world { currentSource = Density xy + , densitySource = addDensity xy } + + | Down <- shift key, Down <- button + = world { currentSource = Velocity xy } + + | Down <- shift key, Up <- button + = case currentSource world of + Velocity x0y0 -> world { currentSource = None + , velocitySource = addVelocity x0y0 xy } + _ -> world { currentSource = None} + + | otherwise + = world { currentSource = None } + + -- Handle key presses + -- + keyboard 'r' Down _ = initialise opt + keyboard 'd' Down _ = world { displayDensity = not (displayDensity world) } + keyboard 'v' Down _ = world { displayVelocity = not (displayVelocity world) } + keyboard _ _ _ = world + + -- As the mouse moves, keep inserting density sources, or adding source + -- velocities + -- + motion xy = + case currentSource world of + Density _ -> world { currentSource = Density xy + , densitySource = addDensity xy } + Velocity x0y0 -> world { currentSource = Velocity xy + , velocitySource = addVelocity x0y0 xy } + _ -> world + -- + addDensity (x,y) = (Z:.y:.x, density) : densitySource world + + addVelocity (x0,y0) (x1,y1) = let u = fromIntegral (x1-x0) + v = fromIntegral (y1-y0) + in (Z:.y0:.x0, (u * velocity, v * velocity)) : velocitySource world + -- + density = get inputDensity opt + velocity = get inputVelocity opt + zoom = fromIntegral $ get displayScale opt + width = fromIntegral $ get simulationWidth opt + height = fromIntegral $ get simulationHeight opt + scaleX = width / (width * zoom + 1) + scaleY = height / (height * zoom + 1) + coord (u,v) = ( truncate $ u * scaleX + width /2 + , truncate $ v * scaleY + height/2) + diff --git a/examples/fluid/src-acc/Fluid.hs b/examples/fluid/src-acc/Fluid.hs new file mode 100644 --- /dev/null +++ b/examples/fluid/src-acc/Fluid.hs @@ -0,0 +1,172 @@ +-- +-- Fluid simulation +-- + +module Fluid ( + + Simulation, fluid + +) where + +import Type +import Prelude as P +import Data.Array.Accelerate as A + + +type Simulation + = Acc ( DensitySource -- locations to add density sources + , VelocitySource -- locations to add velocity sources + , DensityField -- the current density field + , VelocityField ) -- the current velocity field + -> Acc ( DensityField, VelocityField ) + + +-- A fluid simulation +-- +fluid :: Int -> Timestep -> Viscosity -> Diffusion -> Simulation +fluid steps dt dp dn inputs = + let (ds, vs, df, vf) = A.unlift inputs + vf' = velocity steps dt dp vs vf + df' = density steps dt dn ds vf' df + in + A.lift (df', vf') + +-- The velocity over a timestep evolves due to three causes: +-- 1. the addition of forces +-- 2. viscous diffusion +-- 3. self-advection +-- +velocity + :: Int + -> Timestep + -> Viscosity + -> Acc VelocitySource + -> Acc VelocityField + -> Acc VelocityField +velocity steps dt dp vs + = project steps + . (\vf' -> advect dt vf' vf') + . project steps + . diffuse steps dt dp + . inject vs + + +-- Ensure the velocity field conserves mass +-- +project :: Int -> Acc VelocityField -> Acc VelocityField +project steps vf = A.stencil2 poisson (A.Constant zero) vf (A.Constant zero) p + where + grad = A.stencil divF (A.Constant zero) vf + p1 = A.stencil2 pF (A.Constant zero) grad (A.Constant zero) + p = foldl1 (.) (P.replicate steps p1) grad + + poisson :: A.Stencil3x3 Velocity -> A.Stencil3x3 Float -> Exp Velocity + poisson (_,(_,uv,_),_) ((_,t,_), (l,_,r), (_,b,_)) = uv .-. 0.5 .*. A.lift (r-l, b-t) + + divF :: A.Stencil3x3 Velocity -> Exp Float + divF ((_,t,_), (l,_,r), (_,b,_)) = -0.5 * (A.fst r - A.fst l + A.snd b - A.snd t) + + pF :: A.Stencil3x3 Float -> A.Stencil3x3 Float -> Exp Float + pF (_,(_,x,_),_) ((_,t,_), (l,_,r), (_,b,_)) = 0.25 * (x + l + t + r + b) + + +-- The density over a timestep evolves due to three causes: +-- 1. the addition of source particles +-- 2. self-diffusion +-- 3. motion through the velocity field +-- +density + :: Int + -> Timestep + -> Diffusion + -> Acc DensitySource + -> Acc VelocityField + -> Acc DensityField + -> Acc DensityField +density steps dt dn ds vf + = advect dt vf + . diffuse steps dt dn + . inject ds + + +-- Inject sources into the field +-- +-- TLM: sources should be a vector of (index, value) pairs, but no fusion means +-- that extracting the components for permute (via unzip) is extra work. +-- +inject + :: FieldElt e + => Acc (Vector Index, Vector e) + -> Acc (Field e) + -> Acc (Field e) +inject source field = + let (is, ps) = A.unlift source + in A.size ps ==* 0 + ?| ( field, A.permute (.+.) field (is A.!) ps ) + + +diffuse + :: FieldElt e + => Int + -> Timestep + -> Diffusion + -> Acc (Field e) + -> Acc (Field e) +diffuse steps dt dn df0 = + a ==* 0 + ?| ( df0 , foldl1 (.) (P.replicate steps diffuse1) df0 ) + where + a = A.constant dt * A.constant dn * (A.fromIntegral (A.size df0)) + c = 1 + 4*a + + diffuse1 df = A.stencil2 relax (A.Constant zero) df0 (A.Constant zero) df + + relax :: FieldElt e => A.Stencil3x3 e -> A.Stencil3x3 e -> Exp e + relax (_,(_,x0,_),_) ((_,t,_), (l,_,r), (_,b,_)) = (x0 .+. a .*. (l.+.t.+.r.+.b)) ./. c + + +advect + :: FieldElt e + => Timestep + -> Acc VelocityField + -> Acc (Field e) + -> Acc (Field e) +advect dt vf df = A.generate sh backtrace + where + sh = A.shape vf + Z :. h :. w = A.unlift sh + width = A.fromIntegral w + height = A.fromIntegral h + + backtrace ix = s0.*.(t0.*.d00 .+. t1.*.d10) .+. s1.*.(t0.*.d01 .+. t1.*.d11) + where + Z:.j:.i = A.unlift ix + (u, v) = A.unlift (vf A.! ix) + + -- backtrack densities based on velocity field + clamp z = A.max (-0.5) . A.min (z + 0.5) + x = width `clamp` (A.fromIntegral i - A.constant dt * width * u) + y = height `clamp` (A.fromIntegral j - A.constant dt * height * v) + + -- discrete locations surrounding point + i0 = A.truncate (x + 1) - 1 + j0 = A.truncate (y + 1) - 1 + i1 = i0 + 1 + j1 = j0 + 1 + + -- weighting based on location between the discrete points + s1 = x - A.fromIntegral i0 + t1 = y - A.fromIntegral j0 + s0 = 1 - s1 + t0 = 1 - t1 + + -- read the density values surrounding the calculated advection point + get ix'@(Z :. j' :. i') + = (j' <* 0 ||* i' <* 0 ||* j' >=* h ||* i' >=* w) + ? (A.constant zero, df A.! A.lift ix') + + d00 = get (Z :. j0 :. i0) + d10 = get (Z :. j1 :. i0) + d01 = get (Z :. j0 :. i1) + d11 = get (Z :. j1 :. i1) + diff --git a/examples/fluid/src-acc/Main.hs b/examples/fluid/src-acc/Main.hs new file mode 100644 --- /dev/null +++ b/examples/fluid/src-acc/Main.hs @@ -0,0 +1,91 @@ +{-# LANGUAGE CPP #-} +-- +-- A stable fluid simulation +-- +-- Jos Stam, "Real-time Fluid Dynamics for Games" +-- + +module Main where + +import Config +import ParseArgs +import World +import Fluid +import Event +import Data.Label +import Criterion.Main ( defaultMainWith, bench, whnf ) +import Control.Exception +import System.Environment +import Graphics.Gloss.Interface.IO.Game + +import Prelude as P +import Data.Array.Accelerate as A + + +main :: IO () +main = do + argv <- getArgs + (c,crit,noms) <- parseArgs optHelp optBackend options defaults header footer argv + opt <- initialiseConfig c + + let -- configuration parameters + -- + width = get simulationWidth opt * get displayScale opt + height = get simulationHeight opt * get displayScale opt + steps = get simulationSteps opt + fps = get displayFramerate opt + dp = get viscosity opt + dn = get diffusion opt + dt = get timestep opt + backend = get optBackend opt + + -- Prepare user-input density and velocity sources + -- + sources s = let (ix, ss) = P.unzip s + sh = Z :. length ix + in ( A.fromList sh ix, A.fromList sh ss ) + + -- for benchmarking + -- + force w = + indexArray (densityField w) (Z:.0:.0) `seq` + indexArray (velocityField w) (Z:.0:.0) `seq` w + + -- Prepare to execute the next step of the simulation. + -- + -- Critically, we use the run1 execution form to ensure we bypass all + -- front-end conversion phases. + -- + simulate world = + let step = run1 backend (fluid steps dt dp dn) + ds = sources (densitySource world) + vs = sources (velocitySource world) + (df', vf') = step ( ds, vs, densityField world, velocityField world ) + in + force $ world { densityField = df', velocityField = vf' + , densitySource = [], velocitySource = [] } + + -- warming up... + -- + initialWorld <- evaluate (initialise opt) + _ <- evaluate (simulate initialWorld) + +#ifndef ACCELERATE_ENABLE_GUI + if True +#else + if get optBench opt +#endif + -- benchmark + then withArgs noms $ defaultMainWith crit (return ()) + [ bench "fluid" $ whnf simulate initialWorld ] + + -- simulate + else playIO + (InWindow "accelerate-fluid" (width, height) (10, 20)) + black -- background colour + fps -- display framerate + initialWorld -- initial state of the simulation + (render opt) -- render world state into a picture + (\e -> return . react opt e) -- handle user events + (\_ -> return . simulate) -- one step of the simulation + diff --git a/examples/fluid/src-acc/Type.hs b/examples/fluid/src-acc/Type.hs new file mode 100644 --- /dev/null +++ b/examples/fluid/src-acc/Type.hs @@ -0,0 +1,67 @@ +{-# LANGUAGE FlexibleInstances #-} +-- +-- Types used throughout the simulation +-- + +module Type ( + + Timestep, Viscosity, Diffusion, Index, Density, Velocity, + Field, FieldElt(..), DensityField, VelocityField, DensitySource, VelocitySource, + + RGBA, Image + +) where + +import Data.Word +import Data.Array.Accelerate + +type Timestep = Float +type Viscosity = Float +type Diffusion = Float +type Index = DIM2 +type Density = Float +type Velocity = (Float, Float) + +type Field a = Array DIM2 a +type DensityField = Field Density +type VelocityField = Field Velocity + +type Source a = (Vector Index, Vector a) +type DensitySource = Source Density +type VelocitySource = Source Velocity + +type RGBA = Word32 +type Image a = Array DIM2 a + + +infixl 6 .+. +infixl 6 .-. +infixl 7 .*. +infixl 7 ./. + +class Elt e => FieldElt e where + zero :: e + (.+.) :: Exp e -> Exp e -> Exp e + (.-.) :: Exp e -> Exp e -> Exp e + (.*.) :: Exp Float -> Exp e -> Exp e + (./.) :: Exp e -> Exp Float -> Exp e + +instance FieldElt Density where + zero = 0 + (.+.) = (+) + (.-.) = (-) + (.*.) = (*) + (./.) = (/) + +instance FieldElt Velocity where + zero = (0, 0) + (.+.) = app2 (+) + (.-.) = app2 (-) + c .*. xy = let (x,y) = unlift xy in lift (c*x, c*y) + xy ./. c = let (x,y) = unlift xy in lift (x/c, y/c) + +app2 :: Elt e => (Exp e -> Exp e -> Exp e) -> Exp (e,e) -> Exp (e,e) -> Exp (e,e) +app2 f xu yv = let (x,u) = unlift xu + (y,v) = unlift yv + in lift (f x y, f u v) + diff --git a/examples/fluid/src-acc/World.hs b/examples/fluid/src-acc/World.hs new file mode 100644 --- /dev/null +++ b/examples/fluid/src-acc/World.hs @@ -0,0 +1,122 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE PatternGuards #-} +-- +-- Displaying the world state +-- + +module World (World(..), Source(..), initialise, render) where + +import Type +import Config + +import Data.Word +import Data.Label +import Graphics.Gloss.Interface.IO.Game +import Foreign.Ptr +import Foreign.ForeignPtr +import Foreign.Storable +import Foreign.Marshal.Alloc +import Data.Array.Accelerate ( Z(..), (:.)(..) ) +import Data.Array.Accelerate.Array.Sugar ( Array(..) ) + +import qualified Data.Array.Accelerate as A +import qualified Data.Array.Accelerate.Array.Data as A + + +data World = World + { + -- current state of the simulation + densityField :: !DensityField + , velocityField :: !VelocityField + + -- user input + , densitySource :: [(Index, Density)] + , velocitySource :: [(Index, Velocity)] + , currentSource :: Source + , displayDensity :: Bool + , displayVelocity :: Bool + } + deriving Show + +data Source = Density (Int, Int) + | Velocity (Int, Int) + | None + deriving (Eq, Show) + + +-- Initialisation -------------------------------------------------------------- +-- -------------- -- + +initialise :: Config -> World +initialise opt = + World + { densityField = get initialDensity opt + , velocityField = get initialVelocity opt + , densitySource = [] + , velocitySource = [] + , currentSource = None + , displayDensity = True + , displayVelocity = False + } + + +-- Rendering ------------------------------------------------------------------- +-- --------- -- + +render :: Config -> World -> IO Picture +render opt world = do + den <- if displayDensity world + then renderDensity $ densityField world + else return blank + + vel <- if displayVelocity world + then renderVelocity $ velocityField world + else return blank + -- + return $ Scale zoom zoom $ Pictures [ den, vel ] + where + zoom = fromIntegral $ get displayScale opt + + +renderDensity :: DensityField -> IO Picture +renderDensity df@(Array _ ad) = do + dst <- mallocBytes (n*4) + fill 0 src dst + fptr <- newForeignPtr finalizerFree dst + return $ bitmapOfForeignPtr w h fptr False + where + ((),src) = A.ptrsOfArrayData ad + Z:.h:.w = A.arrayShape df + n = h*w + colour !f = let c = 0 `max` f `min` 1 + in floatToWord8 (255*c) + + fill !i !s !d | i >= n = return () + | otherwise = do c <- colour `fmap` peek s + poke d 0xFF -- A + pokeByteOff d 1 c -- B + pokeByteOff d 2 c -- G + pokeByteOff d 3 c -- R + fill (i+1) (plusPtr s 4) (plusPtr d 4) + + +renderVelocity :: VelocityField -> IO Picture +renderVelocity vf + = return + $ Translate (fromIntegral $ -w `div` 2) (fromIntegral $ -h `div` 2) + $ Pictures [ field (x,y) | y <- [2,7..h], x <- [2,7..w] ] + where + Z:.h:.w = A.arrayShape vf + field (x0,y0) = + let x = fromIntegral x0 + y = fromIntegral y0 + (u,v) = A.indexArray vf (Z:.y0:.x0) + in Color red $ Line [ (x,y), (x+ 25*u, y+ 25*v) ] + + +-- Float to Word8 conversion because the one in the GHC libraries doesn't have +-- enough specialisations and goes via Integer. +{-# INLINE floatToWord8 #-} +floatToWord8 :: Float -> Word8 +floatToWord8 f = fromIntegral (truncate f :: Int) + diff --git a/examples/fluid/src/Config.hs b/examples/fluid/src/Config.hs deleted file mode 100644 --- a/examples/fluid/src/Config.hs +++ /dev/null @@ -1,168 +0,0 @@ -{-# LANGUAGE CPP, TemplateHaskell, PatternGuards #-} --- --- Configuration parameters --- - -module Config ( - - Options, - viscosity, diffusion, densityBMP, velocityBMP, simulationWidth, - simulationHeight, displayScale, displayFramerate, optBench, - - processArgs, run, run1 - -) where - -import Data.Label -import System.Exit -import System.Console.GetOpt -import Data.Array.Accelerate ( Acc, Arrays ) -import qualified Data.Array.Accelerate.Interpreter as I -#ifdef ACCELERATE_CUDA_BACKEND -import qualified Data.Array.Accelerate.CUDA as CUDA -#endif - -data Backend - = Interpreter -#ifdef ACCELERATE_CUDA_BACKEND - | CUDA -#endif - deriving (Show, Bounded) - - -data Options = Options - { - -- simulation - _viscosity :: Float - , _diffusion :: Float - - -- visualisation - , _densityBMP :: Maybe FilePath - , _velocityBMP :: Maybe FilePath - , _simulationWidth :: Int - , _simulationHeight :: Int - , _displayScale :: Int - , _displayFramerate :: Int - - -- misc - , _optBackend :: Backend - , _optBench :: Bool - , _optHelp :: Bool - } - deriving Show - -$(mkLabels [''Options]) - -defaultOptions :: Options -defaultOptions = Options - { _viscosity = 0 - , _diffusion = 0 - - , _densityBMP = Nothing - , _velocityBMP = Nothing - , _simulationWidth = 250 - , _simulationHeight = 250 - , _displayScale = 2 - , _displayFramerate = 20 - - , _optBackend = maxBound -#ifdef ACCELERATE_ENABLE_GUI - , _optBench = False -#else - , _optBench = True -#endif - , _optHelp = False - } - --- Execute an Accelerate expression using the selected backend --- -run :: Arrays a => Options -> Acc a -> a -run opts = case _optBackend opts of - Interpreter -> I.run -#ifdef ACCELERATE_CUDA_BACKEND - CUDA -> CUDA.run -#endif - -run1 :: (Arrays a, Arrays b) => Options -> (Acc a -> Acc b) -> a -> b -run1 opts f = case _optBackend opts of - Interpreter -> head . I.stream f . return -#ifdef ACCELERATE_CUDA_BACKEND - CUDA -> CUDA.run1 f -#endif - - --- Process command line arguments --- -options :: [OptDescr (Options -> Options)] -options = - [ Option [] ["interpreter"] (NoArg (set optBackend Interpreter)) - "reference implementation (sequential)" - -#ifdef ACCELERATE_CUDA_BACKEND - , Option [] ["cuda"] (NoArg (set optBackend CUDA)) - "implementation for NVIDIA GPUs (parallel)" - -#endif - , Option [] ["viscosity"] (ReqArg (set viscosity . read) "FLOAT") - $ "viscosity for velocity dampening " ++ def viscosity - - , Option [] ["diffusion"] (ReqArg (set diffusion . read) "FLOAT") - $ "diffusion rate for mass dispersion " ++ def diffusion - - , Option [] ["width"] (ReqArg (set simulationWidth . read) "INT") - $ "grid width for simulation " ++ def simulationWidth - - , Option [] ["height"] (ReqArg (set simulationHeight . read) "INT") - $ "grid height for simulation " ++ def simulationHeight - - , Option [] ["scale"] (ReqArg (set displayScale . read) "INT") - $ "feature size of visualisation " ++ def displayScale - - , Option [] ["framerate"] (ReqArg (set displayFramerate . read) "INT") - $ "frame rate for visualisation " ++ def displayFramerate - - , Option [] ["bmp-density"] (ReqArg (set densityBMP . Just) "FILE.bmp") - "file for initial fluid density" - - , Option [] ["bmp-velocity"] (ReqArg (set velocityBMP . Just) "FILE.bmp") - "file for initial fluid velocity" - - -- - , Option [] ["benchmark"] (NoArg (set optBench True)) - "benchmark instead of displaying animation (False)" - - , Option "h?" ["help"] (NoArg (set optHelp True)) - "show help message" - ] - where - parens s = "(" ++ s ++ ")" - def f = parens (show (get f defaultOptions)) - - -processArgs :: [String] -> IO (Options, [String]) -processArgs argv = - case getOpt' Permute options argv of - (o,_,n,[]) -> case foldl (flip id) defaultOptions o of - opts | False <- get optHelp opts -> return (opts, n) - opts | True <- get optBench opts -> return (opts, "--help":n) - _ -> putStrLn (helpMsg []) >> exitSuccess - (_,_,_,err) -> error (helpMsg err) - where - helpMsg err = concat err ++ usageInfo header options ++ footer - header = unlines - [ "accelerate-fluid (c) 2011 Trevor L. McDonell" - , "" - , "Usage: accelerate-fluid [OPTIONS]" - ] - - footer = unlines - [ "" - , "Runtime usage:" - , "" - , " click add density sources to the image" - , " shift-click add velocity sources" - , " r reset the image" - , " d toggle display of density field" - , " v toggle display of velocity field lines" - ] - diff --git a/examples/fluid/src/Event.hs b/examples/fluid/src/Event.hs deleted file mode 100644 --- a/examples/fluid/src/Event.hs +++ /dev/null @@ -1,81 +0,0 @@ -{-# LANGUAGE PatternGuards #-} --- --- Event handling --- - -module Event where - -import Config -import World -import Data.Label -import Graphics.Gloss.Interface.Pure.Game -import Data.Array.Accelerate ( Z(..), (:.)(..) ) - - --- Event locations are returned as window coordinates, where the origin is in --- the centre of the window and increases to the right and up. If the simulation --- size is (100,100) with scale factor of 4, then the event coordinates are --- returned in the range [-200,200]. --- -react :: Options -> Event -> World -> IO World -react opt event world = - case event of - EventKey (Char c) s m _ -> keyboard c s m - EventKey (MouseButton LeftButton) s m uv -> return $ mouse m s (coord uv) - EventMotion uv -> return $ motion (coord uv) - _ -> return $ world - where - -- Inject a new density source when the left button is clicked. - -- - -- If the shift key is held, remember the location and add a new velocity - -- source between the old and new points as the mouse moves. - -- - mouse key button xy - | Up <- shift key, Down <- button - = world { currentSource = Density xy - , densitySource = addDensity xy } - - | Down <- shift key, Down <- button - = world { currentSource = Velocity xy } - - | Down <- shift key, Up <- button - = case currentSource world of - Velocity x0y0 -> world { currentSource = None - , velocitySource = addVelocity x0y0 xy } - _ -> world { currentSource = None} - - | otherwise - = world { currentSource = None } - - -- Handle key presses - -- - keyboard 'r' Down _ = initialise opt - keyboard 'd' Down _ = return $ world { displayDensity = not (displayDensity world) } - keyboard 'v' Down _ = return $ world { displayVelocity = not (displayVelocity world) } - keyboard _ _ _ = return world - - -- As the mouse moves, keep inserting density sources, or adding source - -- velocities - -- - motion xy = - case currentSource world of - Density _ -> world { currentSource = Density xy - , densitySource = addDensity xy } - Velocity x0y0 -> world { currentSource = Velocity xy - , velocitySource = addVelocity x0y0 xy } - _ -> world - -- - addDensity (x,y) = (Z:.y:.x, 1) : densitySource world - - addVelocity (x0,y0) (x1,y1) = let u = fromIntegral (x1-x0) - v = fromIntegral (y1-y0) - in (Z:.y0:.x0, (u * width, v * height)) : velocitySource world - -- - zoom = fromIntegral $ get displayScale opt - width = fromIntegral $ get simulationWidth opt - height = fromIntegral $ get simulationHeight opt - scaleX = width / (width * zoom + 1) - scaleY = height / (height * zoom + 1) - coord (u,v) = ( truncate $ u * scaleX + width /2 - , truncate $ v * scaleY + height/2) - diff --git a/examples/fluid/src/Fluid.hs b/examples/fluid/src/Fluid.hs deleted file mode 100644 --- a/examples/fluid/src/Fluid.hs +++ /dev/null @@ -1,169 +0,0 @@ --- --- Fluid simulation --- - -module Fluid ( - - Simulation, fluid - -) where - -import Type - -import Data.Array.Accelerate ( Z(..), (:.)(..), Exp, Acc, Scalar, Vector, (?|), (==*) ) -import qualified Data.Array.Accelerate as A - - -type Simulation - = Acc ( Scalar Timestep -- time to evolve the simulation - , DensitySource -- locations to add density sources - , VelocitySource -- locations to add velocity sources - , DensityField -- the current density field - , VelocityField ) -- the current velocity field - -> Acc ( DensityField, VelocityField ) - - --- A fluid simulation --- -fluid :: Viscosity -> Diffusion -> Simulation -fluid dp dn inputs = - let (dt, ds, vs, df, vf) = A.unlift inputs - df' = density dn dt ds vf df - vf' = velocity dp dt vs vf - in - A.lift (df', vf') - --- The velocity over a timestep evolves due to three causes: --- 1. the addition of forces --- 2. viscous diffusion --- 3. self-advection --- -velocity - :: Viscosity - -> Acc (Scalar Timestep) - -> Acc VelocitySource - -> Acc VelocityField - -> Acc VelocityField -velocity dp dt vs vf0 - = project - . advect dt vf0 - . project - . diffuse dp dt - $ inject vs vf0 - - --- Ensure the velocity field conserves mass --- -project :: Acc VelocityField -> Acc VelocityField -project vf = A.stencil2 poisson A.Mirror vf A.Mirror p - where - steps = 20 - grad = A.stencil divF A.Mirror vf - p1 = A.stencil2 pF (A.Constant 0) grad A.Mirror - p = foldl1 (.) (replicate steps p1) grad - - poisson :: A.Stencil3x3 Velocity -> A.Stencil3x3 Float -> Exp Velocity - poisson (_,(_,uv,_),_) ((_,t,_), (l,_,r), (_,b,_)) = uv .-. 0.5 .*. A.lift (r-l, t-b) - - divF :: A.Stencil3x3 Velocity -> Exp Float - divF ((_,t,_), (l,_,r), (_,b,_)) = -0.5 * (A.fst r - A.fst l + A.snd t - A.snd b) - - pF :: A.Stencil3x3 Float -> A.Stencil3x3 Float -> Exp Float - pF (_,(_,x,_),_) ((_,t,_), (l,_,r), (_,b,_)) = 0.25 * (x + l + t + r + b) - - - --- The density over a timestep evolves due to three causes: --- 1. the addition of source particles --- 2. self-diffusion --- 3. motion through the velocity field --- -density - :: Diffusion - -> Acc (Scalar Timestep) - -> Acc DensitySource - -> Acc VelocityField - -> Acc DensityField - -> Acc DensityField -density dn dt ds vf - = advect dt vf - . diffuse dn dt - . inject ds - - --- Inject sources into the field --- --- TLM: sources should be a vector of (index, value) pairs, but no fusion means --- that extracting the components for permute (via unzip) is extra work. --- -inject - :: FieldElt e - => Acc (Vector Index, Vector e) - -> Acc (Field e) - -> Acc (Field e) -inject source field = - let (is, ps) = A.unlift source - in A.size ps ==* 0 - ?| ( field, A.permute (.+.) field (is A.!) ps ) - - -diffuse - :: FieldElt e - => Diffusion - -> Acc (Scalar Timestep) - -> Acc (Field e) - -> Acc (Field e) -diffuse dn dt df0 = - a ==* 0 - ?| ( df0 , foldl1 (.) (replicate steps diffuse1) df0 ) - where - steps = 20 - a = A.the dt * A.constant dn * (A.fromIntegral (A.size df0)) - c = 1 + 4*a - - diffuse1 df = A.stencil2 relax (A.Constant zero) df0 A.Mirror df - - relax :: FieldElt e => A.Stencil3x3 e -> A.Stencil3x3 e -> Exp e - relax (_,(_,x0,_),_) ((_,t,_), (l,_,r), (_,b,_)) = (x0 .+. a .*. (l.+.t.+.r.+.b)) ./. c - - -advect - :: FieldElt e - => Acc (Scalar Timestep) - -> Acc VelocityField - -> Acc (Field e) - -> Acc (Field e) -advect dt' vf df = A.generate sh backtrace - where - dt = A.the dt' - sh = A.shape vf - Z :. h :. w = A.unlift sh - - backtrace ix = s0.*.(t0.*.d00 .+. t1.*.d10) .+. s1.*.(t0.*.d01 .+. t1.*.d11) - where - Z:.j:.i = A.unlift ix - (u, v) = A.unlift (vf A.! ix) - - -- backtrack densities based on velocity field - clamp z = A.max 0.5 . A.min (A.fromIntegral z - 1.5) - x = w `clamp` (A.fromIntegral i - dt * u) - y = h `clamp` (A.fromIntegral j - dt * v) - - -- discrete locations surrounding point - i0 = A.truncate x - j0 = A.truncate y - i1 = i0 + 1 - j1 = j0 + 1 - - -- weighting based on location between the discrete points - s1 = x - A.fromIntegral i0 - t1 = y - A.fromIntegral j0 - s0 = 1 - s1 - t0 = 1 - t1 - - -- read the density values surrounding the calculated advection point - d00 = df A.! A.lift (Z :. j0 :. i0) - d10 = df A.! A.lift (Z :. j1 :. i0) - d01 = df A.! A.lift (Z :. j0 :. i1) - d11 = df A.! A.lift (Z :. j1 :. i1) - diff --git a/examples/fluid/src/Main.hs b/examples/fluid/src/Main.hs deleted file mode 100644 --- a/examples/fluid/src/Main.hs +++ /dev/null @@ -1,84 +0,0 @@ -{-# LANGUAGE CPP #-} --- --- A stable fluid simulation --- --- Jos Stam, "Real-time Fluid Dynamics for Games" --- - -module Main where - -import Config -import World -import Fluid -import Event -import Data.Label -import Criterion.Main -import System.Environment -import Graphics.Gloss.Interface.IO.Game - -import Prelude as P -import Data.Array.Accelerate as A - - -main :: IO () -main = do - (opt,noms) <- processArgs =<< getArgs - let -- configuration parameters - -- - width = get simulationWidth opt * get displayScale opt - height = get simulationHeight opt * get displayScale opt - fps = get displayFramerate opt - dp = get viscosity opt - dn = get diffusion opt - - -- Prepare user-input density and velocity sources - -- - sources s = let (ix, ss) = P.unzip s - sh = Z :. length ix - in ( A.fromList sh ix, A.fromList sh ss ) - - -- for benchmarking - -- - force world = - indexArray (densityField world) (Z:.0:.0) `seq` - indexArray (velocityField world) (Z:.0:.0) `seq` () - - -- Prepare to execute the next step of the simulation. - -- - -- Critically, we use the run1 execution form to ensure we bypass all - -- front-end conversion phases. - -- - simulate timestep world = - let step = run1 opt (fluid dp dn) - dt = A.fromList Z [timestep] - ds = sources (densitySource world) - vs = sources (velocitySource world) - (df', vf') = step ( dt, ds, vs, densityField world, velocityField world ) - in - return $ world { densityField = df', velocityField = vf' - , densitySource = [], velocitySource = [] } - - -- warming up... - initialWorld <- initialise opt - _ <- simulate 0.1 initialWorld - - if get optBench opt - -- benchmark - then withArgs noms $ defaultMain - [ bench "fluid" $ whnfIO (force `fmap` simulate 1.0 initialWorld) ] - -#ifndef ACCELERATE_ENABLE_GUI - else return () - -#else - -- simulate - else playIO - (InWindow "accelerate-fluid" (width, height) (10, 20)) - black -- background colour - fps -- display framerate - initialWorld -- initial state of the simulation - (render opt) -- render world state into a picture - (react opt) -- handle user events - simulate -- one step of the simulation -#endif - diff --git a/examples/fluid/src/Type.hs b/examples/fluid/src/Type.hs deleted file mode 100644 --- a/examples/fluid/src/Type.hs +++ /dev/null @@ -1,67 +0,0 @@ -{-# LANGUAGE FlexibleInstances #-} --- --- Types used throughout the simulation --- - -module Type ( - - Timestep, Viscosity, Diffusion, Index, Density, Velocity, - Field, FieldElt(..), DensityField, VelocityField, DensitySource, VelocitySource, - - RGBA, Image - -) where - -import Data.Word -import Data.Array.Accelerate - -type Timestep = Float -type Viscosity = Float -type Diffusion = Float -type Index = DIM2 -type Density = Float -type Velocity = (Float, Float) - -type Field a = Array DIM2 a -type DensityField = Field Density -type VelocityField = Field Velocity - -type Source a = (Vector Index, Vector a) -type DensitySource = Source Density -type VelocitySource = Source Velocity - -type RGBA = Word32 -type Image a = Array DIM2 a - - -infixl 6 .+. -infixl 6 .-. -infixl 7 .*. -infixl 7 ./. - -class Elt e => FieldElt e where - zero :: e - (.+.) :: Exp e -> Exp e -> Exp e - (.-.) :: Exp e -> Exp e -> Exp e - (.*.) :: Exp Float -> Exp e -> Exp e - (./.) :: Exp e -> Exp Float -> Exp e - -instance FieldElt Density where - zero = 0 - (.+.) = (+) - (.-.) = (-) - (.*.) = (*) - (./.) = (/) - -instance FieldElt Velocity where - zero = (0, 0) - (.+.) = app2 (+) - (.-.) = app2 (-) - c .*. xy = let (x,y) = unlift xy in lift (c*x, c*y) - xy ./. c = let (x,y) = unlift xy in lift (x/c, y/c) - -app2 :: Elt e => (Exp e -> Exp e -> Exp e) -> Exp (e,e) -> Exp (e,e) -> Exp (e,e) -app2 f xu yv = let (x,u) = unlift xu - (y,v) = unlift yv - in lift (f x y, f u v) - diff --git a/examples/fluid/src/World.hs b/examples/fluid/src/World.hs deleted file mode 100644 --- a/examples/fluid/src/World.hs +++ /dev/null @@ -1,198 +0,0 @@ -{-# LANGUAGE BangPatterns #-} -{-# LANGUAGE PatternGuards #-} --- --- Displaying the world state --- - -module World (World(..), Source(..), initialise, render) where - -import Type -import Config - -import Codec.BMP -import Data.Bits -import Data.Int -import Data.Word -import Data.Label -import Control.Monad -import Graphics.Gloss.Interface.IO.Game -import Foreign.Ptr -import Foreign.ForeignPtr -import Foreign.Storable -import Foreign.Marshal.Alloc -import Data.Array.Accelerate ( Z(..), (:.)(..), Exp ) -import Data.Array.Accelerate.Array.Sugar ( Array(..) ) - -import qualified Data.Array.Accelerate.Array.Data as A -import qualified Data.Array.Accelerate.IO as A -import qualified Data.Array.Accelerate as A - - -data World = World - { - -- current state of the simulation - densityField :: !DensityField - , velocityField :: !VelocityField - - -- user input - , densitySource :: [(Index, Density)] - , velocitySource :: [(Index, Velocity)] - , currentSource :: Source - , displayDensity :: Bool - , displayVelocity :: Bool - } - deriving Show - -data Source = Density (Int, Int) - | Velocity (Int, Int) - | None - deriving (Eq, Show) - - --- Initialisation -------------------------------------------------------------- --- -------------- -- - -initialise :: Options -> IO World -initialise opt = - let width = get simulationWidth opt - height = get simulationHeight opt - in do - density <- initialDensity opt width height - velocity <- initialVelocity opt width height - - return $ World - { densityField = density - , velocityField = velocity - , densitySource = [] - , velocitySource = [] - , currentSource = None - , displayDensity = True - , displayVelocity = False - } - - -initialDensity :: Options -> Int -> Int -> IO DensityField -initialDensity opt width height - -- Load the file, and use the luminance value as scalar density - -- - | Just file <- get densityBMP opt - = do - arr <- readImageFromBMP file - let (Z:.h:.w) = A.arrayShape arr - - when (w /= width || h /= height) - $ error "fluid: density-bmp does not match width x height" - - return . run opt $ A.map densityOfRGBA (A.use arr) - - -- No density file given, just set the field to zero - -- - | otherwise - = return $ A.fromList (Z :. height :. width) (repeat 0) - - -initialVelocity :: Options -> Int -> Int -> IO VelocityField -initialVelocity opt width height - -- Load the file, and use the red and green channels for x- and y- velocity - -- values respectively - -- - | Just file <- get velocityBMP opt - = do - arr <- readImageFromBMP file - let (Z:.h:.w) = A.arrayShape arr - - when (w /= width || h /= height) - $ error "fluid: velocity-bmp does not match width x height" - - return . run opt $ A.map velocityOfRGBA (A.use arr) - - -- No density file given, just set to zero - -- - | otherwise - = return $ A.fromList (Z :. height :. width) (repeat (0,0)) - - -readImageFromBMP :: FilePath -> IO (Image RGBA) -readImageFromBMP file = do - bmp <- either (error . show) id `fmap` readBMP file - let (w,h) = bmpDimensions bmp - -- - A.fromByteString (Z :. h :. w) ((), unpackBMPToRGBA32 bmp) - - -densityOfRGBA :: Exp RGBA -> Exp Density -densityOfRGBA rgba = - let b = (0.11 / 255) * A.fromIntegral ((rgba `div` 0x100) .&. 0xFF) - g = (0.59 / 255) * A.fromIntegral ((rgba `div` 0x10000) .&. 0xFF) - r = (0.3 / 255) * A.fromIntegral ((rgba `div` 0x1000000) .&. 0xFF) - in - r + g + b - -velocityOfRGBA :: Exp RGBA -> Exp Velocity -velocityOfRGBA rgba = - let g = A.fromIntegral (-128 + A.fromIntegral ((rgba `div` 0x10000) .&. 0xFF) :: Exp Int32) - r = A.fromIntegral (-128 + A.fromIntegral ((rgba `div` 0x1000000) .&. 0xFF) :: Exp Int32) - in - A.lift (r * 0.001, g * 0.001) - - --- Rendering ------------------------------------------------------------------- --- --------- -- - -render :: Options -> World -> IO Picture -render opt world = do - den <- if displayDensity world - then renderDensity $ densityField world - else return blank - - vel <- if displayVelocity world - then renderVelocity $ velocityField world - else return blank - -- - return $ Scale zoom zoom $ Pictures [ den, vel ] - where - zoom = fromIntegral $ get displayScale opt - - -renderDensity :: DensityField -> IO Picture -renderDensity df@(Array _ ad) = do - dst <- mallocBytes (n*4) - fill 0 src dst - fptr <- newForeignPtr finalizerFree dst - return $ bitmapOfForeignPtr w h fptr False - where - ((),src) = A.ptrsOfArrayData ad - Z:.h:.w = A.arrayShape df - n = h*w - colour !f = let c = 0 `max` f `min` 1 - in floatToWord8 (255*c) - - fill !i !s !d | i >= n = return () - | otherwise = do c <- colour `fmap` peek s - poke d 0xFF -- A - pokeByteOff d 1 c -- B - pokeByteOff d 2 c -- G - pokeByteOff d 3 c -- R - fill (i+1) (plusPtr s 4) (plusPtr d 4) - - -renderVelocity :: VelocityField -> IO Picture -renderVelocity vf - = return - $ Translate (fromIntegral $ -w `div` 2) (fromIntegral $ -h `div` 2) - $ Pictures [ field (x,y) | y <- [2,7..h], x <- [2,7..w] ] - where - Z:.h:.w = A.arrayShape vf - field (x0,y0) = - let x = fromIntegral x0 - y = fromIntegral y0 - (u,v) = A.indexArray vf (Z:.y0:.x0) - in Color red $ Line [ (x,y), (x+u, y+v) ] - - --- Float to Word8 conversion because the one in the GHC libraries doesn't have --- enough specialisations and goes via Integer. -{-# INLINE floatToWord8 #-} -floatToWord8 :: Float -> Word8 -floatToWord8 f = fromIntegral (truncate f :: Int) - diff --git a/examples/hashcat/Config.hs b/examples/hashcat/Config.hs new file mode 100644 --- /dev/null +++ b/examples/hashcat/Config.hs @@ -0,0 +1,79 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE PatternGuards #-} +{-# LANGUAGE TemplateHaskell #-} + +module Config where + +import ParseArgs +import Data.Label + + +data Config + = Config + { + -- How to execute the search + _configBackend :: Backend + + -- Input data + , _configStrings :: [String] + , _configDict :: FilePath + , _configMaxWords :: Maybe Int + , _configSkipWords :: Int + + , _configHelp :: Bool + } + deriving Show + +$(mkLabels [''Config]) + +defaults :: Config +defaults = Config + { + _configBackend = maxBound + , _configStrings = [] + , _configDict = [] + , _configMaxWords = Nothing + , _configSkipWords = 0 + , _configHelp = False + } + + +-- | The set of available command-line options +-- +options :: [OptDescr (Config -> Config)] +options = + [ Option ['s'] [] + (ReqArg (modify configStrings . (:)) "STRING") + "Lookup the plain text of a given checksum" + + , Option ['d'] ["dictionary"] + (ReqArg (set configDict) "FILE") + "Plain text word list to search against" + + , Option ['j'] ["skip-words"] + (ReqArg (set configSkipWords . read) "INT") + "Skip this many entries from the start of the word list" + + , Option ['n'] ["max-words"] + (ReqArg (set configMaxWords . Just . read) "INT") + "Use at most this many words from the list" + + , Option ['h', '?'] ["help"] + (NoArg (set configHelp True)) + "show this help message" + ] + + +-- | Process the command line options +-- +header :: [String] +header = + [ "accelerate-hashcat (c) [2013] The Accelerate Team" + , "" + , "Usage: accelerate-hashcat -d dictionary [OPTIONS] [file ...]" + , "" + ] + +footer :: [String] +footer = [] + diff --git a/examples/hashcat/Digest.hs b/examples/hashcat/Digest.hs new file mode 100644 --- /dev/null +++ b/examples/hashcat/Digest.hs @@ -0,0 +1,97 @@ +{-# LANGUAGE BangPatterns #-} + +module Digest (readDict, extract) + where + +import MD5 +import Config + +import Prelude +import Data.Label +import Control.Monad +import Control.Applicative +import Data.ByteString.Internal ( w2c ) +import qualified Data.Serialize as S +import qualified Data.ByteString as S +import qualified Data.ByteString.Lazy.Char8 as L + +import Data.Array.Accelerate.Array.Data as A +import Data.Array.Accelerate.Array.Sugar as A +import qualified Data.Array.Accelerate as A + + +-- MD5 block sizes +-- +blockSizeBits, blockSizeBytes, blockSizeWords :: Int +blockSizeBits = 512 +blockSizeBytes = blockSizeBits `div` 8 +blockSizeWords = blockSizeBytes `div` 4 + + +-- Create an MD5 block from the given message. This appends The '1' bit to the +-- message, pads the block with zeros until the length in bits is 448, and +-- finally appends the length in bits (mod 2^64). +-- +md5block :: L.ByteString -> S.ByteString +md5block msg = do + let + len = fromIntegral (L.length msg) + lenBits = 8 * fromIntegral len + lenZeroPad + | len + 1 <= blockSizeBytes - 8 = (blockSizeBytes - 8) - (len + 1) + | otherwise = (2 * blockSizeBytes - 8) - (len + 1) + -- + S.runPut $! do + S.putLazyByteString msg + S.putWord8 0x80 + mapM_ S.putWord8 (replicate lenZeroPad 0) + S.putWord64le lenBits + + +-- Create a dictionary of blocks ready to digest from the given bytestring. This +-- reads one entry per line. Because we only do a single MD5 chunk, we discard +-- any entries with (length > blockSizeBytes - 8 = 55) characters. Because of +-- this restriction, we also need to a precise count of how many entries will be +-- in the final array, since it is stored column-major and hence there is no +-- easy way to truncate it. +-- +readDict :: Config -> FilePath -> IO Dictionary +readDict c fp = do + entries <- length . chunk <$> L.readFile fp + blocks <- map md5block . chunk <$> L.readFile fp + + let sh = Z :. blockSizeWords :. entries + (adata,_) = runArrayData $ do + arr <- newArrayData (size sh) + + let go !_ [] = return () + go !n (b:bs) = do + foldM_ (\i w -> do unsafeWriteArrayData arr (toIndex sh (Z:.i:.n)) (fromElt w) + return (i+1)) 0 (bytes b) + go (n+1) bs + + go 0 blocks + return (arr, undefined) + + adata `seq` return $ Array (fromElt sh) adata + + where + chunk = maybe id take (get configMaxWords c) + . drop (get configSkipWords c) + . filter (\w -> fromIntegral (L.length w) < blockSizeBytes - 8) + . L.lines + + bytes = either error id + . S.runGet (replicateM blockSizeWords S.getWord32le) + +-- Extract a word from the dictionary at a given index +-- +extract :: Dictionary -> Int -> L.ByteString +extract dict i + | i > n = error "extract: index too large" + | otherwise = L.takeWhile (/= w2c 0x80) bytes + where + Z :. _ :. n = A.arrayShape dict + bytes = S.runPutLazy $ + forM_ [0 .. blockSizeWords-1] $ \c -> S.putWord32le (dict `A.indexArray` (Z:.c:.i)) + diff --git a/examples/hashcat/MD5.hs b/examples/hashcat/MD5.hs new file mode 100644 --- /dev/null +++ b/examples/hashcat/MD5.hs @@ -0,0 +1,164 @@ +{-# LANGUAGE ViewPatterns #-} + +module MD5 ( + + Dictionary, + hashcat, readMD5, showMD5 + +) where + + +import Prelude as P +import Numeric +import Control.Monad +import Control.Applicative +import Data.Maybe +import Data.List ( foldl' ) +import Data.Bits hiding ( rotateL ) +import qualified Data.Serialize as S +import qualified Data.ByteString as B +import qualified Data.ByteString.Lazy as L +import Data.ByteString.Lex.Integral ( readHexadecimal ) + +import Data.Array.Accelerate as A + + +-- Generate an MD5 hash for every word in the dictionary, and if an entry +-- matches the given unknown md5, returns the index into the dictionary of said +-- match. If not found, this returns (-1). +-- +hashcat :: Acc Dictionary + -> Acc (Scalar MD5) + -> Acc (Scalar Int) +hashcat dict passwd + = reshape (constant Z) + $ permute const res (\ix -> crypt A.! ix `cmp` the passwd ? (constant (Z:.0), ignore)) + (enumFromN (index1 n) 0) + where + n = A.size crypt + res = use $ fromList (Z:.1) [-1] :: Acc (Vector Int) + crypt = md5 dict + + cmp x y = let (x1,x2,x3,x4) = unlift x + (y1,y2,y3,y4) = unlift y + in x1 ==* y1 &&* x2 ==* y2 &&* x3 ==* y3 &&* x4 ==* y4 + + +-- An MD5 round processes 512 bits of the input, as 16 x 32-bit values. We +-- require that the input dictionary of words to hash is of dimension +-- (Z :. 16 :. n), which will then output n MD5 hashes values. +-- +-- NOTE: All values are unsigned 32-bit integers and wrap modulo 2^32 when +-- calculating. All values are little endian. +-- +-- TODO: How to apply salt to passwords? +-- +type MD5 = (Word32, Word32, Word32, Word32) +type Dictionary = Array DIM2 Word32 + +md5 :: Acc Dictionary -> Acc (Vector MD5) +md5 dict + = let n = A.snd . unindex2 $ A.shape dict + in + A.generate (index1 n) (md5Round dict) + + +-- A single round of MD5. For inputs longer than 512-bits we would simply +-- reapply this until we run out of input data to hash. +-- +md5Round :: Acc Dictionary -> Exp DIM1 -> Exp MD5 +md5Round word (unindex1 -> ix) + = lift + $ foldl step (a0,b0,c0,d0) [0..64] + where + step (a,b,c,d) i + | i < 16 = shfl ((b .&. c) .|. ((complement b) .&. d)) + | i < 32 = shfl ((b .&. d) .|. (c .&. (complement d))) + | i < 48 = shfl ((b `xor` c `xor` d)) + | i < 64 = shfl (c `xor` (b .|. (complement d))) + | otherwise = (a+a0,b+b0,c+c0,d+d0) + where + shfl f = (d, b + ((a + f + k i + get i) `rotateL` r i), b, c) + + get :: Int -> Exp Word32 + get i + | i < 16 = word A.! index2 (constant i) ix + | i < 32 = word A.! index2 (constant ((5*i + 1) `rem` 16)) ix + | i < 48 = word A.! index2 (constant ((3*i + 5) `rem` 16)) ix + | otherwise = word A.! index2 (constant ((7*i) `rem` 16)) ix + + -- Initial values. For a multi-round implementation we would initialise the + -- context with these values, and then after applying the round update the + -- context with the addition of the new result. Since we only ever apply one + -- round, we begin with these values and add them again at the end, which + -- seems odd but is correct. + -- + a0, b0, c0, d0 :: Exp Word32 + a0 = 0x67452301 + b0 = 0xefcdab89 + c0 = 0x98badcfe + d0 = 0x10325476 + + -- The binary integer part of sines and cosines (in radians) as constants + -- + -- K[i] := floor(abs(sin(i + 1)) × (2 pow 32)) + -- + k :: Int -> Exp Word32 + k i = constant (ks P.!! i) + where + ks = [ 0xd76aa478, 0xe8c7b756, 0x242070db, 0xc1bdceee + , 0xf57c0faf, 0x4787c62a, 0xa8304613, 0xfd469501 + , 0x698098d8, 0x8b44f7af, 0xffff5bb1, 0x895cd7be + , 0x6b901122, 0xfd987193, 0xa679438e, 0x49b40821 + , 0xf61e2562, 0xc040b340, 0x265e5a51, 0xe9b6c7aa + , 0xd62f105d, 0x02441453, 0xd8a1e681, 0xe7d3fbc8 + , 0x21e1cde6, 0xc33707d6, 0xf4d50d87, 0x455a14ed + , 0xa9e3e905, 0xfcefa3f8, 0x676f02d9, 0x8d2a4c8a + , 0xfffa3942, 0x8771f681, 0x6d9d6122, 0xfde5380c + , 0xa4beea44, 0x4bdecfa9, 0xf6bb4b60, 0xbebfbc70 + , 0x289b7ec6, 0xeaa127fa, 0xd4ef3085, 0x04881d05 + , 0xd9d4d039, 0xe6db99e5, 0x1fa27cf8, 0xc4ac5665 + , 0xf4292244, 0x432aff97, 0xab9423a7, 0xfc93a039 + , 0x655b59c3, 0x8f0ccc92, 0xffeff47d, 0x85845dd1 + , 0x6fa87e4f, 0xfe2ce6e0, 0xa3014314, 0x4e0811a1 + , 0xf7537e82, 0xbd3af235, 0x2ad7d2bb, 0xeb86d391 ] + + -- Specifies the per-round shift amounts + -- + r :: Int -> Exp Int + r i = constant (rs P.!! i) + where + rs = [ 7, 12, 17, 22, 7, 12, 17, 22, 7, 12, 17, 22, 7, 12, 17, 22 + , 5, 9, 14, 20, 5, 9, 14, 20, 5, 9, 14, 20, 5, 9, 14, 20 + , 4, 11, 16, 23, 4, 11, 16, 23, 4, 11, 16, 23, 4, 11, 16, 23 + , 6, 10, 15, 21, 6, 10, 15, 21, 6, 10, 15, 21, 6, 10, 15, 21 ] + + +-- Display the result in big endian hexadecimal format, which is consistent with +-- other MD5 hash implementations. +-- +showMD5 :: MD5 -> String +showMD5 (a,b,c,d) = + let bs = S.runPutLazy $ S.putWord32be d >> S.putWord32be c >> S.putWord32be b >> S.putWord32be a + showsHex str w = let x = showHex w str + in if length x < length str + 2 + then '0':x + else x + in + foldl' showsHex "" (L.unpack bs) + + +-- Read an input hexadecimal string into our MD5 representation. This string is +-- the same as that displayed by the standard md5 utility, and also produced by +-- 'showMD5'. +-- +readMD5 :: L.ByteString -> MD5 +readMD5 = + let get = (,,,) <$> readHex32be <*> readHex32be <*> readHex32be <*> readHex32be + + readHex32be = do + s <- B.concat . P.reverse <$> replicateM 4 (S.getBytes 2) + return . P.fst $ fromMaybe (error "readHex32be: parse failure") (readHexadecimal s) + in + either error id . S.runGetLazy get + diff --git a/examples/hashcat/Main.hs b/examples/hashcat/Main.hs new file mode 100644 --- /dev/null +++ b/examples/hashcat/Main.hs @@ -0,0 +1,80 @@ + +module Main where + +import MD5 +import Digest +import Config +import ParseArgs +import Util + +import Data.Label +import Text.Printf +import Control.Monad +import Control.Applicative +import Criterion.Measurement +import System.IO +import System.Environment +import Data.Array.Accelerate ( Z(..), (:.)(..) ) +import qualified Data.Array.Accelerate as A +import qualified Data.ByteString.Lazy.Char8 as L + + +main :: IO () +main = do + argv <- getArgs + (conf, _cconf, files) <- parseArgs configHelp configBackend options defaults header footer argv + + -- Read the plain text word lists. This creates a vector of MD5 chunks ready + -- for hashing. + -- + putStr "Loading dictionary... " >> hFlush stdout + (tdict, dict) <- time $ readDict conf (get configDict conf) + + let (Z :. _ :. entries) = A.arrayShape dict + putStrLn $ printf "%d words in %s" entries (secs tdict) + + -- Attempt to recover one hash at a time by comparing it to entries in the + -- database. This rehashes the entire word list every time, rather than + -- pre-computing the hashes and performing a lookup. That approach, known as a + -- rainbow table, while much faster when multiple iterations of the hashing + -- function are applied, but is defeated by salting passwords. This is true + -- even if the salt is known, so long as it is unique for each password. + -- + let backend = get configBackend conf + + recover hash = + let abcd = readMD5 hash + idx = run1 backend (hashcat (A.use dict)) (A.fromList Z [abcd]) + -- + in case idx `A.indexArray` Z of + -1 -> Nothing + n -> Just (extract dict n) + + recoverAll :: [L.ByteString] -> IO (Int,Int) + recoverAll = + foldM (\(i,n) h -> maybe (return (i,n+1)) (\t -> showText h t >> return (i+1,n+1)) (recover h)) (0,0) + + showText hash text = do + L.putStr hash >> putStr ": " >> L.putStrLn text + + -- Queue up all the message digests to process + -- + digests <- concat . (map L.pack (get configStrings conf) :) + <$> mapM (\f -> L.lines `fmap` L.readFile f) files + + -- Run the lookup for each unknown hash against the given wordlists. + -- + putStrLn "Beginning recovery..." + (trec, (r, t)) <- time (recoverAll digests) + + -- And print a summary of results + -- + let percent = fromIntegral r / fromIntegral t * 100.0 :: Double + persec = fromIntegral (t * entries) / trec + putStrLn $ printf "\nRecovered %d/%d (%.2f %%) digests in %s, %s" + r t percent + (showFFloatSIBase (Just 2) 1000 trec "s") + (showFFloatSIBase (Just 2) 1000 persec "Hash/sec") + + when (r == t) $ putStrLn "All hashes recovered (:" + diff --git a/examples/mandelbrot/Config.hs b/examples/mandelbrot/Config.hs --- a/examples/mandelbrot/Config.hs +++ b/examples/mandelbrot/Config.hs @@ -1,34 +1,18 @@ {-# LANGUAGE CPP #-} -{-# LANGUAGE PatternGuards #-} {-# LANGUAGE TemplateHaskell #-} -module Config ( - - Options, optBackend, optSize, optLimit, optBench, - processArgs, run, run1 - -) where +module Config where +import ParseArgs import Data.Label -import System.Exit -import System.Console.GetOpt -import Data.Array.Accelerate ( Arrays, Acc ) -import qualified Data.Array.Accelerate.Interpreter as Interp -#ifdef ACCELERATE_CUDA_BACKEND -import qualified Data.Array.Accelerate.CUDA as CUDA -#endif -data Backend = Interpreter -#ifdef ACCELERATE_CUDA_BACKEND - | CUDA -#endif - deriving (Bounded, Show) - data Options = Options { _optBackend :: Backend - , _optSize :: Int + , _optWidth :: Int + , _optHeight :: Int , _optLimit :: Int + , _optFramerate :: Int , _optBench :: Bool , _optHelp :: Bool } @@ -36,11 +20,13 @@ $(mkLabels [''Options]) -defaultOptions :: Options -defaultOptions = Options +defaults :: Options +defaults = Options { _optBackend = maxBound - , _optSize = 512 + , _optWidth = 800 + , _optHeight = 600 , _optLimit = 255 + , _optFramerate = 25 #ifdef ACCELERATE_ENABLE_GUI , _optBench = False #else @@ -49,48 +35,34 @@ , _optHelp = False } - -run :: Arrays a => Options -> Acc a -> a -run opts = case _optBackend opts of - Interpreter -> Interp.run -#ifdef ACCELERATE_CUDA_BACKEND - CUDA -> CUDA.run -#endif - -run1 :: (Arrays a, Arrays b) => Options -> (Acc a -> Acc b) -> a -> b -run1 opts f = case _optBackend opts of - Interpreter -> head . Interp.stream f . return -#ifdef ACCELERATE_CUDA_BACKEND - CUDA -> CUDA.run1 f -#endif - - options :: [OptDescr (Options -> Options)] options = - [ Option [] ["interpreter"] (NoArg (set optBackend Interpreter)) "reference implementation (sequential)" -#ifdef ACCELERATE_CUDA_BACKEND - , Option [] ["cuda"] (NoArg (set optBackend CUDA)) "implementation for NVIDIA GPUs (parallel)" -#endif - , Option [] ["size"] (ReqArg (set optSize . read) "INT") "visualisation size (512)" + [ Option [] ["width"] (ReqArg (set optWidth . read) "INT") "visualisation width (800)" + , Option [] ["height"] (ReqArg (set optHeight . read) "INT") "visualisation height (600)" , Option [] ["limit"] (ReqArg (set optLimit . read) "INT") "iteration limit for escape (255)" + , Option [] ["framerate"] (ReqArg (set optFramerate . read) "INT")"visualisation framerate (10)" + , Option [] ["static"] (NoArg (set optFramerate 0)) "do not animate the image" , Option [] ["benchmark"] (NoArg (set optBench True)) "benchmark instead of displaying animation (False)" , Option "h?" ["help"] (NoArg (set optHelp True)) "show help message" ] -processArgs :: [String] -> IO (Options, [String]) -processArgs argv = - case getOpt' Permute options argv of - (o,_,n,[]) -> case foldl (flip id) defaultOptions o of - opts | False <- get optHelp opts -> return (opts, n) - opts | True <- get optBench opts -> return (opts, "--help":n) - _ -> putStrLn (helpMsg []) >> exitSuccess - (_,_,_,err) -> error (helpMsg err) - where - helpMsg err = concat err ++ usageInfo header options - header = unlines - [ "accelerate-mandelbrot (c) [2011..2012] The Accelerate Team" - , "" - , "Usage: accelerate-mandelbrot [OPTIONS]" - ] +header :: [String] +header = + [ "accelerate-mandelbrot (c) [2011..2013] The Accelerate Team" + , "" + , "Usage: accelerate-mandelbrot [OPTIONS]" + , "" + ] + +footer :: [String] +footer = + [ "" + , "Runtime usage:" + , " arrows translate display" + , " z ; zoom in" + , " x q zoom out" + , " f single precision calculations" + , " d double precision calculations (if supported)" + ] diff --git a/examples/mandelbrot/Main.hs b/examples/mandelbrot/Main.hs --- a/examples/mandelbrot/Main.hs +++ b/examples/mandelbrot/Main.hs @@ -1,125 +1,34 @@ {-# LANGUAGE CPP #-} -- --- A Mandelbrot set generator. Submitted by Simon Marlow as part of Issue #49. +-- A Mandelbrot set generator. +-- Originally submitted by Simon Marlow as part of Issue #49. -- - +import World import Config +import ParseArgs + import Data.Label import Control.Monad -import Control.Exception import Foreign.Ptr import Foreign.ForeignPtr import System.IO.Unsafe -import System.Environment -import Criterion.Main ( defaultMain, bench, whnf ) +import System.Environment ( getArgs, withArgs ) +import Criterion.Main ( defaultMainWith, bench, whnf ) import Data.Array.Accelerate.Array.Data ( ptrsOfArrayData ) import Data.Array.Accelerate.Array.Sugar ( Array(..) ) import Prelude as P -import Data.Array.Accelerate as A hiding ( size ) +import Data.Array.Accelerate as A import qualified Graphics.Gloss as G --- Mandelbrot Set -------------------------------------------------------------- - -type F = Float -- Double not supported on all devices - -type Complex = (F,F) -type ComplexPlane = Array DIM2 Complex - -mandelbrot :: F -> F -> F -> F -> Int -> Int -> Int -> Acc (Array DIM2 (F,F,Int)) -mandelbrot x y x' y' screenX screenY depth - = foldl (flip ($)) zs0 (P.take depth (repeat go)) - where - cs = genPlane x y x' y' screenX screenY - zs0 = mkinit cs - - go :: Acc (Array DIM2 (F,F,Int)) -> Acc (Array DIM2 (F,F,Int)) - go = A.zipWith iter cs - - -genPlane :: F -> F - -> F -> F - -> Int - -> Int - -> Acc ComplexPlane -genPlane lowx lowy highx highy viewx viewy - = generate (constant (Z:.viewy:.viewx)) - (\ix -> let pr = unindex2 ix - x = A.fromIntegral (A.fst pr) - y = A.fromIntegral (A.snd pr) - in - lift ( elowx + (x * exsize) / eviewx - , elowy + (y * eysize) / eviewy)) - where - elowx, elowy, exsize, eysize, eviewx, eviewy :: Exp F - - elowx = constant lowx - elowy = constant lowy - - exsize = constant (highx - lowx) - eysize = constant (highy - lowy) - - eviewx = constant (P.fromIntegral viewx) - eviewy = constant (P.fromIntegral viewy) - - -next :: Exp Complex -> Exp Complex -> Exp Complex -next c z = c `plus` (z `times` z) - - -plus :: Exp Complex -> Exp Complex -> Exp Complex -plus = lift2 f - where f :: (Exp F, Exp F) -> (Exp F, Exp F) -> (Exp F, Exp F) - f (x1,y1) (x2,y2) = (x1+x2,y1+y2) - -times :: Exp Complex -> Exp Complex -> Exp Complex -times = lift2 f - where f :: (Exp F, Exp F) -> (Exp F, Exp F) -> (Exp F, Exp F) - f (x,y) (x',y') = (x*x'-y*y', x*y'+y*x') - -dot :: Exp Complex -> Exp F -dot = lift1 f - where f :: (Exp F, Exp F) -> Exp F - f (x,y) = x*x + y*y - - -iter :: Exp Complex -> Exp (F,F,Int) -> Exp (F,F,Int) -iter c z = f (unlift z) - where - f :: (Exp F, Exp F, Exp Int) -> Exp (F,F,Int) - f (x,y,i) = - (dot z' >* 4.0) ? ( lift (x,y,i) - , lift (A.fst z', A.snd z', i+1) ) - where z' = A.curry (next c) x y - - -mkinit :: Acc ComplexPlane -> Acc (Array DIM2 (F,F,Int)) -mkinit cs = A.map (lift1 f) cs - where f :: (Exp F, Exp F) -> (Exp F, Exp F, Exp Int) - f (x,y) = (x,y,0) - - --- Rendering ------------------------------------------------------------------- - -type RGBA = Word32 - -prettyRGBA :: Exp Int -> Exp (F, F, Int) -> Exp RGBA -prettyRGBA lIMIT s' = r + g + b + a - where - (_, _, s) = unlift s' :: (Exp F, Exp F, Exp Int) - t = A.fromIntegral $ ((lIMIT - s) * 255) `quot` lIMIT - r = (t `mod` 128 + 64) * 0x1000000 - g = (t * 2 `mod` 128 + 64) * 0x10000 - b = (t * 3 `mod` 256 ) * 0x100 - a = 0xFF - +-- Main ------------------------------------------------------------------------ -makePicture :: Options -> Int -> Acc (Array DIM2 (F, F, Int)) -> G.Picture -makePicture opt limit zs = pic +makePicture :: World -> G.Picture +makePicture world = pic where - arrPixels = run opt $ A.map (prettyRGBA (constant limit)) zs + arrPixels = renderWorld world (Z:.h:.w) = arrayShape arrPixels {-# NOINLINE rawData #-} @@ -128,47 +37,54 @@ in unsafePerformIO $ newForeignPtr_ (castPtr ptr) - pic = G.bitmapOfForeignPtr h w rawData False + pic = G.bitmapOfForeignPtr w h rawData False --- Main ------------------------------------------------------------------------ - main :: IO () main - = do (config, nops) <- processArgs =<< getArgs - let size = get optSize config - limit = get optLimit config - -- - x = -0.25 -- should get this from command line as well - y = -1.0 - x' = 0.0 - y' = -0.75 - -- - force arr = indexArray arr (Z:.0:.0) `seq` arr - image = makePicture config limit - $ mandelbrot x y x' y' size size limit + = do + argv <- getArgs + (conf, cconf, rest) <- parseArgs optHelp optBackend options defaults header footer argv - void $ evaluate image + let world = initialWorld conf view + fps = get optFramerate conf + width = get optWidth conf + height = get optHeight conf - if get optBench config - then withArgs nops $ defaultMain - [ bench "mandelbrot" $ - whnf (force . run config . mandelbrot x y x' y' size size) limit ] + -- Centre coordinates: Re(c) = -0.7; Im(c) = 0 + -- View width: 3.0769 + -- + view = (-2.23, -1.15, 0.83, 1.15) -#ifndef ACCELERATE_ENABLE_GUI - else return () + force arr = indexArray arr (Z:.0:.0) `seq` arr -#elif MIN_VERSION_gloss(1,6,0) - else G.display - (G.InWindow "Mandelbrot" (size, size) (10, 10)) + mandel + | get optBench conf + = withArgs rest $ defaultMainWith cconf (return ()) + [ bench "mandelbrot" $ whnf (force . renderWorld) world ] +#ifdef ACCELERATE_ENABLE_GUI + | fps == 0 + = G.display + (G.InWindow "Mandelbrot" (width, height) (10, 10)) G.black - image -#else - else G.displayInWindow - "Mandelbrot" - (size, size) - (10, 10) + (makePicture world) + + | fps > 0 + = G.play + (G.InWindow "Mandelbrot" (width, height) (10, 10)) G.black - image + fps + world + makePicture + (react conf) + (const refocus) #endif + | otherwise + = return () + + + unless (P.null rest) $ + putStrLn $ "Warning: unrecognized options: " ++ show rest + + mandel diff --git a/examples/mandelbrot/Mandel.hs b/examples/mandelbrot/Mandel.hs new file mode 100644 --- /dev/null +++ b/examples/mandelbrot/Mandel.hs @@ -0,0 +1,195 @@ +{-# LANGUAGE ScopedTypeVariables #-} +-- +-- A Mandelbrot set generator. +-- Originally submitted by Simon Marlow as part of Issue #49. +-- +module Mandel ( + + -- Types + View, Render, Bitmap, + + -- Pretty pictures + mandelbrot, prettyRGBA, + +) where + +import Prelude as P +import Data.Array.Accelerate as A +import Data.Array.Accelerate.IO as A + +-- Types ----------------------------------------------------------------------- + +-- Current view into the complex plane +type View a = (a, a, a, a) + +-- Complex numbers +type Complex a = (a, a) +type ComplexPlane a = Array DIM2 (Complex a) + +-- Image data +type Bitmap = Array DIM2 RGBA32 + +-- Action to render a frame +type Render a = Scalar (View a) -> Bitmap + + +-- Mandelbrot Set -------------------------------------------------------------- + +-- Compute the mandelbrot as repeated application of the recurrence relation: +-- +-- Z_{n+1} = c + Z_n^2 +-- +-- This is applied until the iteration limit is reached, regardless of whether +-- or not the pixel has diverged. This implies wasted work, but no thread +-- divergence. +-- +mandelbrot + :: forall a. (Elt a, IsFloating a) + => Int + -> Int + -> Int + -> Acc (Scalar (View a)) + -> Acc (Array DIM2 (Complex a,Int)) +mandelbrot screenX screenY depth view + = foldr ($) zs0 (P.take depth (repeat go)) + where + cs = genPlane screenX screenY view + zs0 = mkinit cs + + go :: Acc (Array DIM2 (Complex a, Int)) -> Acc (Array DIM2 (Complex a, Int)) + go = A.zipWith iter cs + + +genPlane :: (Elt a, IsFloating a) => Int -> Int -> Acc (Scalar (View a)) -> Acc (ComplexPlane a) +genPlane screenX screenY view + = generate (constant (Z:.screenY:.screenX)) + (\ix -> let pr = unindex2 ix + x = A.fromIntegral (A.snd pr :: Exp Int) + y = A.fromIntegral (A.fst pr :: Exp Int) + in + lift ( xmin + (x * sizex) / viewx + , ymin + (y * sizey) / viewy)) + where + (xmin,ymin,xmax,ymax) = unlift (the view) + + sizex = xmax - xmin + sizey = ymax - ymin + + viewx = constant (P.fromIntegral screenX) + viewy = constant (P.fromIntegral screenY) + + +next :: (Elt a, IsNum a) => Exp (Complex a) -> Exp (Complex a) -> Exp (Complex a) +next c z = c `plus` (z `times` z) + + +plus :: (Elt a, IsNum a) => Exp (Complex a) -> Exp (Complex a) -> Exp (Complex a) +plus = lift2 f + where f :: (Elt a, IsNum a) => (Exp a, Exp a) -> (Exp a, Exp a) -> (Exp a, Exp a) + f (x1,y1) (x2,y2) = (x1+x2, y1+y2) + +times :: forall a. (Elt a, IsNum a) => Exp (Complex a) -> Exp (Complex a) -> Exp (Complex a) +times = lift2 f + where f :: (Elt a, IsNum a) => (Exp a, Exp a) -> (Exp a, Exp a) -> (Exp a, Exp a) + f (x,y) (x',y') = (x*x'-y*y', x*y'+y*x') + +dot :: forall a. (Elt a, IsNum a) => Exp (Complex a) -> Exp a +dot = lift1 f + where f :: (Elt a, IsNum a) => (Exp a, Exp a) -> Exp a + f (x,y) = x*x + y*y + + +iter :: forall a. (Elt a, IsNum a) => Exp (Complex a) -> Exp (Complex a, Int) -> Exp (Complex a, Int) +iter c zi = f (A.fst zi) (A.snd zi) + where + f :: Exp (Complex a) -> Exp Int -> Exp (Complex a, Int) + f z i = + let z' = next c z + in (dot z' >* 4) ? ( zi , lift (z', i+1) ) + + +mkinit :: Elt a => Acc (ComplexPlane a) -> Acc (Array DIM2 (Complex a, Int)) +mkinit cs = A.zip cs (A.fill (A.shape cs) 0) + + +-- Rendering ------------------------------------------------------------------- + +prettyRGBA :: forall a. (Elt a, IsFloating a) => Exp Int -> Exp (Complex a, Int) -> Exp RGBA32 +prettyRGBA lIMIT s = + let cmax = A.fromIntegral lIMIT + c = A.fromIntegral (A.snd s) + in + c ==* cmax ? ( 0xFF000000, escapeToColour (cmax - c) ) + +-- Directly convert the iteration count on escape to a colour. The base set +-- (x,y,z) yields a dark background with light highlights. +-- +escapeToColour :: Exp Int -> Exp RGBA32 +escapeToColour m = constant 0xFFFFFFFF - (packRGBA32 $ lift (x,y,z,w)) + where + x = constant 0 + w = A.fromIntegral (3 * m) + z = A.fromIntegral (5 * m) + y = A.fromIntegral (7 * m) + + +{-- +-- A simple colour scheme +-- +prettyRGBA :: Elt a => Exp Int -> Exp (Complex a, Int) -> Exp RGBA32 +prettyRGBA lIMIT s' = r + g + b + a + where + s = A.snd s' + t = A.fromIntegral $ ((lIMIT - s) * 255) `quot` lIMIT + r = (t `rem` 128 + 64) * 0x1000000 + g = (t * 2 `rem` 128 + 64) * 0x10000 + b = (t * 3 `rem` 256 ) * 0x100 + a = 0xFF +--} +{-- +prettyRGBA :: forall a. (Elt a, IsFloating a) => Exp Int -> Exp (Complex a, Int) -> Exp RGBA32 +prettyRGBA lIMIT s = + let cmax = A.fromIntegral lIMIT :: Exp a + c = A.fromIntegral (A.snd s) + in + c >* 0.98 * cmax ? ( 0xFF000000, rampColourHotToCold 0 cmax c ) + +-- Standard Hot-to-Cold hypsometric colour ramp. Colour sequence is +-- Red, Yellow, Green, Cyan, Blue +-- +rampColourHotToCold + :: (Elt a, IsFloating a) + => Exp a -- ^ minimum value of the range + -> Exp a -- ^ maximum value of the range + -> Exp a -- ^ data value + -> Exp RGBA32 +rampColourHotToCold vmin vmax vNotNorm + = let v = vmin `A.max` vNotNorm `A.min` vmax + dv = vmax - vmin + -- + result = v <* vmin + 0.28 * dv + ? ( lift ( constant 0.0 + , 4 * (v-vmin) / dv + , constant 1.0 + , constant 1.0 ) + + , v <* vmin + 0.5 * dv + ? ( lift ( constant 0.0 + , constant 1.0 + , 1 + 4 * (vmin + 0.25 * dv - v) / dv + , constant 1.0 ) + + , v <* vmin + 0.75 * dv + ? ( lift ( 4 * (v - vmin - 0.5 * dv) / dv + , constant 1.0 + , constant 0.0 + , constant 1.0 ) + + , lift ( constant 1.0 + , 1 + 4 * (vmin + 0.75 * dv - v) / dv + , constant 0.0 + , constant 1.0 ) + ))) + in + rgba32OfFloat result +--} diff --git a/examples/mandelbrot/World.hs b/examples/mandelbrot/World.hs new file mode 100644 --- /dev/null +++ b/examples/mandelbrot/World.hs @@ -0,0 +1,174 @@ +{-# LANGUAGE GADTs #-} +{-# LANGUAGE TypeOperators #-} + +module World ( + + -- Types + World, + + -- Updating the World state + renderWorld, initialWorld, refocus, react + +) where + +import Mandel +import Config +import ParseArgs + +import Data.Char +import Data.Label +import Data.Array.Accelerate as A +import Graphics.Gloss.Interface.Pure.Game hiding ( translate, scale ) + + +-- World state +-- ----------- + +data Zoom = In | Out +data Move = Fwd | Rev + +data Precision = Float | Double + +data World where + World :: (Elt a, RealFloat a) + => View a + -> Render a + -> Maybe Zoom + -> Maybe Move -- horizontal movement + -> Maybe Move -- vertical movement + -> World + + +-- Render the picture +-- +renderWorld :: World -> Bitmap +renderWorld (World view render _ _ _) = render $ A.fromList Z [view] + + +-- Initialise the World state +-- +initialWorld :: Options -> View Float -> World +initialWorld config view + = setPrecisionOfWorld Float config + $ World view undefined Nothing Nothing Nothing + + +-- Reset the rendering routines to compute with the specified precision +-- +setPrecisionOfWorld :: Precision -> Options -> World -> World +setPrecisionOfWorld f config (World p _ z h v) + = let + width = get optWidth config + height = get optHeight config + limit = get optLimit config + backend = get optBackend config + + render :: (Elt a, IsFloating a) => Render a + render = run1 backend + $ A.map (prettyRGBA (constant limit)) + . mandelbrot width height limit + + in case f of + Float -> World (convertView p :: View Float) render z h v + Double -> World (convertView p :: View Double) render z h v + + +-- Event handling +-- -------------- + +-- Refocus the viewport by adjusting the limits of the x- and y- range of the +-- display, based on the current key state. +-- +refocus :: World -> World +refocus = move . zoom + where + -- translate the display + -- + move :: World -> World + move world@(World viewport r z h v) + = World (translate (dy,dx) viewport) r z h v + where + dx = case get horizontal world of + Nothing -> 0 + Just Fwd -> 0.025 + Just Rev -> -0.025 + + dy = case get vertical world of + Nothing -> 0 + Just Fwd -> 0.025 + Just Rev -> -0.025 + + translate (j,i) (x,y,x',y') = + let sizex = x' - x + sizey = y' - y + in (x+i*sizex, y+j*sizey, x'+i*sizex, y'+j*sizey) + + -- zoom the display in or out + -- + zoom :: World -> World + zoom world@(World viewport r z h v) + = World (scale s viewport) r z h v + where + s = case get zooming world of + Nothing -> 1 + Just In -> 0.975 + Just Out -> 1.025 + + scale alpha (x,y,x',y') = + let dx = sizex * alpha / 2 + dy = sizey * alpha / 2 + sizex = x' - x + sizey = y' - y + midx = x + sizex / 2 + midy = y + sizey / 2 + in (midx - dx, midy - dy, midx + dx, midy + dy) + + +-- Event locations are returned as window coordinates, where the origin is in +-- the centre of the window and increases to the right and up. If the simulation +-- size is (100,100) with scale factor of 4, then the event coordinates are +-- returned in the range [-200,200]. +-- +react :: Options -> Event -> World -> World +react opt event world + = case event of + EventKey (Char c) s _ _ -> char (toLower c) s world + EventKey (SpecialKey c) s _ _ -> special c s world + _ -> world + where + char ';' = toggle zooming In + char 'z' = toggle zooming In + char 'q' = toggle zooming Out + char 'x' = toggle zooming Out + char 'd' = precision Double + char 'f' = precision Float + char _ = const id + + special KeyUp = toggle vertical Fwd + special KeyDown = toggle vertical Rev + special KeyRight = toggle horizontal Fwd + special KeyLeft = toggle horizontal Rev + special _ = const id + + toggle f x Down = set f (Just x) + toggle f _ Up = set f Nothing + + precision f Down = setPrecisionOfWorld f opt + precision _ _ = id + + +-- Miscellaneous +-- ------------- + +zooming :: World :-> Maybe Zoom +zooming = lens (\(World _ _ z _ _) -> z) (\z (World p r _ h v) -> World p r z h v) + +horizontal :: World :-> Maybe Move +horizontal = lens (\(World _ _ _ h _) -> h) (\h (World p r z _ v) -> World p r z h v) + +vertical :: World :-> Maybe Move +vertical = lens (\(World _ _ _ _ v) -> v) (\v (World p r z h _) -> World p r z h v) + +convertView :: (Real a, Fractional b) => View a -> View b +convertView (x,y,x',y') = (realToFrac x, realToFrac y, realToFrac x', realToFrac y') + diff --git a/examples/n-body/Common/Body.hs b/examples/n-body/Common/Body.hs new file mode 100644 --- /dev/null +++ b/examples/n-body/Common/Body.hs @@ -0,0 +1,180 @@ +-- +-- Simulation bodies with mass +-- + +module Common.Body ( + + -- * Types + Velocity, Accel, PointMass, Body, + + -- * Calculations + accel, advanceBody, + + -- ** Getters + pointMassOfBody, velocityOfBody, accelOfBody, positionOfPointMass, + massOfPointMass, + + -- ** Setters + unitBody, setMassOfBody, setAccelOfBody, setStartVelOfBody, + +) where + +import Common.Type +import Common.Util + +import Data.Array.Accelerate as A + + +-- Acceleration ---------------------------------------------------------------- +-- +-- | Calculate the acceleration on a point due to some other point as an inverse +-- separation-squared relation. The force on a body i caused by its +-- gravitational attraction to a body j is given by: +-- +-- m_i m_j r_ij +-- f_ij = G --------- . ------ +-- |r_ij|^2 |r_ij| +-- +-- The total force on a body F is given by this interaction to all other +-- bodies: +-- +-- F_i = Sum_{i/=j} f_ij +-- +-- m_j r_ij +-- = G m_i . Sum_{j/=i} -------- +-- |r_ij|^3 +-- +-- As the bodies approach each other, the force between them grows without +-- bound, which is an undesirable situation for numerical integration. Since +-- collisions between bodies are precluded, a softening factor (epsilon^2 > 0) +-- is added: +-- +-- m_j r_ij +-- F_i = G m_i . Sum ----------------------------- +-- ( |r_ij|^2 + epsilon^2) ^ 3/2 +-- +-- Note that the condition (i /= j) is no longer required, because (f_ii = 0) +-- when (epsilon^2 > 0). The softening factor models the interaction between +-- two Plummer point masses: bodies that behave as if they were spherical +-- galaxies (and thus may pass through each other). +-- +-- To integrate over time, we need the acceleration (a_i = F_i / m_i), and so +-- the above can be simplified by removing m_i from the RHS. This function +-- computes the component of the Sum for two bodies i and j. +-- +accel :: Exp R -- ^ Smoothing parameter + -> Exp Body -- ^ The point being accelerated + -> Exp Body -- ^ Neighbouring point + -> Exp Accel + +accel epsilon bodyi bodyj = s *. r + where + pmi = pointMassOfBody bodyi + pmj = pointMassOfBody bodyj + mj = massOfPointMass pmj + + r = positionOfPointMass pmj .-. positionOfPointMass pmi + rsqr = dot r r + epsilon * epsilon + invr = 1 / sqrt rsqr + invr3 = invr * invr * invr + + s = mj * invr3 + + +-- Body ------------------------------------------------------------------------ +-- + +-- | Make a stationary Body of unit mass +-- +unitBody :: Exp (Vec R) -> Exp Body +unitBody pos = lift (pointmass, vec 0, vec 0) + where + pointmass = lift (pos, constant 1) :: Exp PointMass + + +-- | Take the Velocity of a Body +-- +velocityOfBody :: Exp Body -> Exp Velocity +velocityOfBody body = vel + where + (_, vel, _) = unlift body :: (Exp PointMass, Exp Velocity, Exp Accel) + + +-- | Take the Acceleration of a Body +-- +accelOfBody :: Exp Body -> Exp Accel +accelOfBody body = acc + where + (_, _, acc) = unlift body :: (Exp PointMass, Exp Velocity, Exp Accel) + + +-- | Take the PointMass of a Body +-- +pointMassOfBody :: Exp Body -> Exp PointMass +pointMassOfBody body = mp + where + (mp, _, _) = unlift body :: (Exp PointMass, Exp Velocity, Exp Accel) + + +-- | Take the position or mass of a PointMass +-- +positionOfPointMass :: Exp PointMass -> Exp Position +positionOfPointMass = A.fst + +massOfPointMass :: Exp PointMass -> Exp Mass +massOfPointMass = A.snd + + +-- | Set the mass of a Body. +-- +setMassOfBody :: Exp Mass -> Exp Body -> Exp Body +setMassOfBody mass body = lift (pointmass, vel, acc) + where + vel = velocityOfBody body + acc = accelOfBody body + pos = positionOfPointMass (pointMassOfBody body) + pointmass = lift (pos, mass) :: Exp PointMass + + +-- | Set the acceleration of a Body. +-- +setAccelOfBody :: Exp Accel -> Exp Body -> Exp Body +setAccelOfBody acc body = lift (pm, vel, acc) + where + pm = pointMassOfBody body + vel = velocityOfBody body + + +-- | Set the starting velocity of a Body. +-- It is set to rotate around the origin, with the speed proportional +-- to the sqrt of the distance from it. This seems to make nice simulations. +-- +setStartVelOfBody :: Exp R -> Exp Body -> Exp Body +setStartVelOfBody startVel body = lift (pm, vel'', acc) + where + pm = pointMassOfBody body + acc = accelOfBody body + pos = positionOfPointMass pm + + pos' = normalise pos + vel' = lift (y', -x', z') + vel'' = (sqrt (magnitude pos) * startVel) *. vel' + + (x',y',z') = unlift pos' :: Vec (Exp R) + + +-- | Advance a body forwards in time. +-- +advanceBody :: Exp Time -> Exp Body -> Exp Body +advanceBody time body = lift ( pm', vel', acc ) + where + pm = pointMassOfBody body + pos = positionOfPointMass pm + vel = velocityOfBody body + acc = accelOfBody body + mass = massOfPointMass pm + + pm' = lift (pos', mass) :: Exp PointMass + pos' = pos .+. time *. vel + vel' = vel .+. time *. acc + diff --git a/examples/n-body/Common/Dump.hs b/examples/n-body/Common/Dump.hs new file mode 100644 --- /dev/null +++ b/examples/n-body/Common/Dump.hs @@ -0,0 +1,52 @@ + +module Common.Dump ( dumpWorld, printWorld ) + where + +import Common.Body +import Common.World + +import System.IO +import Data.List +import Text.Printf +import qualified Data.Array.Accelerate as A + + +-- | Dump the bodies in a world to a file. +-- +dumpWorld :: World -> FilePath -> IO () +dumpWorld world filePath + = withFile filePath WriteMode + $ \h -> hWriteWorld h world + +-- | Print the bodies in a world to stdout +-- +printWorld :: World -> IO () +printWorld world + = hWriteWorld stdout world + + +-- | Write a world to a handle +-- +hWriteWorld :: Handle -> World -> IO () +hWriteWorld h world + = mapM_ (hWriteBody h) + $ A.toList + $ worldBodies world + +-- | Write a single body to a handle +-- +hWriteBody :: Handle -> Body -> IO () +hWriteBody h (((px, py), mass), (vx, vy), (ax, ay)) + = hPutStrLn h + $ concat + $ ( (padRc 8 ' ' $ show mass) + : " " + : (intersperse " " $ map (\f -> printf "%15.8f" f) [ px, py, vx, vy, ax, ay ])) + + +-- | Right justify a doc, padding with a given character. +-- +padRc :: Int -> Char -> String -> String +padRc n c str + = replicate (n - length str) c ++ str + diff --git a/examples/n-body/Common/Tree.hs b/examples/n-body/Common/Tree.hs new file mode 100644 --- /dev/null +++ b/examples/n-body/Common/Tree.hs @@ -0,0 +1,7 @@ +-- +-- A Barns-Hut tree +-- + +module Common.Tree + where + diff --git a/examples/n-body/Common/Type.hs b/examples/n-body/Common/Type.hs new file mode 100644 --- /dev/null +++ b/examples/n-body/Common/Type.hs @@ -0,0 +1,39 @@ + +module Common.Type where + + +-- Types ----------------------------------------------------------------------- +-- We're using tuples instead of ADTs and defining Elt instances +-- + +-- | Not all compute devices support double precision +-- +type R = Float + +-- | A data point in space +-- +type Vec a = (a, a, a) + +-- | Units of time +-- +type Time = R + +-- | The velocity of a point. +-- +type Velocity = Vec R + +-- | The acceleration of a point. +-- +type Accel = Vec R + +-- | A point in 2D space with its mass. +-- +type Mass = R +type Position = Vec R +type PointMass = (Position, Mass) + +-- | Bodies consist of a Position and Mass, but also carry their velocity and +-- acceleration between steps of the simulation. +-- +type Body = (PointMass, Velocity, Accel) + diff --git a/examples/n-body/Common/Util.hs b/examples/n-body/Common/Util.hs new file mode 100644 --- /dev/null +++ b/examples/n-body/Common/Util.hs @@ -0,0 +1,76 @@ + +module Common.Util ( + + magnitude, dot, normalise, vec, + + (.+.), (.-.), (.*.), (+.), (-.), (*.), + +) where + +import Common.Type +import Data.Array.Accelerate as A + +-- | The magnitude of a vector. +-- +magnitude :: (Elt a, IsFloating a) => Exp (Vec a) -> Exp a +magnitude v = sqrt (dot v v) + + +-- | Dot product of a vector +-- +dot :: (Elt a, IsNum a) => Exp (Vec a) -> Exp (Vec a) -> Exp a +dot v1 v2 + = let (x1,y1,z1) = unlift v1 + (x2,y2,z2) = unlift v2 + in + x1 * x2 + y1 * y2 + z1 * z2 + + +-- | Normalise a vector, so it has a magnitude of 1. +-- +normalise :: (Elt a, IsFloating a) => Exp (Vec a) -> Exp (Vec a) +normalise v = (1 / magnitude v) *. v + +-- | Replicate a value into a vector +-- +vec :: Elt a => Exp a -> Exp (Vec a) +vec x = lift (x,x,x) + +-- | Basic arithmetic component-wise +-- +infixl 7 .*. +infixl 6 .+. +infixl 6 .-. + +(.+.), (.-.), (.*.) :: (Elt a, IsNum a) => Exp (Vec a) -> Exp (Vec a) -> Exp (Vec a) +(.+.) = vzipWith (+) +(.-.) = vzipWith (-) +(.*.) = vzipWith (*) + +-- | Apply a scalar value component-wise to each element of the vector +-- +infixl 7 *. +infixl 6 +. +infixl 6 -. + +(+.), (-.), (*.) :: (Elt a, IsNum a) => Exp a -> Exp (Vec a) -> Exp (Vec a) +(+.) c = vmap (c+) +(-.) c = vmap (c-) +(*.) c = vmap (c*) + +-- | Arithmetic lifted to our vector type. As far as possible, want to gloss +-- over whether we are calculating in 2D or 3D. +-- +vmap :: (Elt a, Elt b) => (Exp a -> Exp b) -> Exp (Vec a) -> Exp (Vec b) +vmap f v + = let (x1,y1,z1) = unlift v + in + lift (f x1, f y1, f z1) + +vzipWith :: (Elt a, Elt b, Elt c) => (Exp a -> Exp b -> Exp c) -> Exp (Vec a) -> Exp (Vec b) -> Exp (Vec c) +vzipWith f v1 v2 + = let (x1,y1,z1) = unlift v1 + (x2,y2,z2) = unlift v2 + in + lift (f x1 x2, f y1 y2, f z1 z2) + diff --git a/examples/n-body/Common/World.hs b/examples/n-body/Common/World.hs new file mode 100644 --- /dev/null +++ b/examples/n-body/Common/World.hs @@ -0,0 +1,63 @@ + +module Common.World ( + + World(..), advanceBodies, advanceWorld + +) where + +import Common.Type +import Common.Util +import Common.Body + +import Data.Array.Accelerate as A + + +data World + = World + { + worldBodies :: !(Vector Body) -- ^ Bodies in the simulation + , worldSteps :: {-# UNPACK #-} !Int -- ^ Number of steps taken in the simulation so far + , worldTime :: {-# UNPACK #-} !Time -- ^ Current simulation time + } + + +-- | Move bodies under the influence of acceleration +-- +advanceBodies + :: (Acc (Vector Body) -> Acc (Vector Accel)) -- ^ Function to compute accelerations at each point + -> Acc (Scalar Time) -- ^ Time step + -> Acc (Vector Body) -- ^ Bodies + -> Acc (Vector Body) +advanceBodies calcAccels timeStep bodies + = let + -- Calculate the accelerations on each body. + accels = calcAccels bodies + + -- Apply the accelerations to the bodies and advance them + advance b a = let m = massOfPointMass (pointMassOfBody b) + a' = m *. a + in advanceBody (the timeStep) (setAccelOfBody a' b) + in + A.zipWith advance bodies accels + + +-- | Advance a cluster of bodies forward in time +-- +advanceWorld + :: (Scalar Time -> Vector Body -> Vector Body) -- ^ Function to update body positions + -> Time + -> World + -> World +advanceWorld advance timeStep world + = let + -- Update the bodies + bodies' = advance (fromList Z [timeStep]) (worldBodies world) + + -- Update the world + steps' = worldSteps world + 1 + time' = worldTime world + timeStep + + in world { worldBodies = bodies' + , worldSteps = steps' + , worldTime = time' } + diff --git a/examples/n-body/Config.hs b/examples/n-body/Config.hs new file mode 100644 --- /dev/null +++ b/examples/n-body/Config.hs @@ -0,0 +1,163 @@ +{-# LANGUAGE TemplateHaskell #-} + +module Config where + +import ParseArgs +import Common.Type + +import Data.Char +import Data.List +import Data.Label + + +data Solver = Naive | BarnsHut + deriving (Enum, Bounded, Show) + + +data Config + = Config + { + -- How to execute the simulation + _configBackend :: Backend + , _configSolver :: Solver + + -- How to present the output + , _configWindowSize :: Int + , _configShouldDrawTree :: Bool + , _configRate :: Int + + -- System setup + , _configBodyCount :: Int + , _configBodyMass :: R + , _configTimeStep :: R + , _configEpsilon :: R + + -- Initial conditions + , _configStartDiscSize :: R + , _configStartSpeed :: R + + -- Terminating conditions + , _configMaxSteps :: Maybe Int + , _configBenchmark :: Bool + , _configHelp :: Bool + + -- Dump final particle locations to file + , _configDumpFinal :: Maybe FilePath + } + deriving Show + +$(mkLabels [''Config]) + + +defaults :: Config +defaults = Config + { + _configBackend = maxBound + , _configSolver = Naive -- no barns-hut yet! + + , _configWindowSize = 1000 + , _configShouldDrawTree = False -- no barns-hut yet! + , _configRate = 30 + + , _configBodyCount = 1000 + , _configBodyMass = 40 + , _configTimeStep = 1 + , _configEpsilon = 50 + + , _configStartDiscSize = 500 + , _configStartSpeed = 1 + + , _configMaxSteps = Nothing + , _configBenchmark = False + , _configHelp = False + + , _configDumpFinal = Nothing + } + + +-- | The set of available command-line options +-- +options :: [OptDescr (Config -> Config)] +options = + [ Option ['s'] ["solver"] + (ReqArg (set configSolver . solver) "ALGORITHM") + ("solver to use, one of: " ++ intercalate ", " (map show [minBound .. maxBound :: Solver])) + + , Option [] ["size"] + (ReqArg (set configWindowSize . read) "INT") + (describe configWindowSize "visualisation size") + + , Option [] ["framerate"] + (ReqArg (set configRate . read) "INT") + (describe configRate "visualisation frame rate") + + , Option [] ["draw-tree"] + (NoArg (set configShouldDrawTree True)) + "draw the Barns-Hut quad tree" + + , Option ['n'] ["bodies"] + (ReqArg (set configBodyCount . read) "INT") + (describe configBodyCount "number of bodies in the simulation") + + , Option [] ["mass"] + (ReqArg (set configBodyMass . read) "FLOAT") + (describe configBodyMass "mass of each body") + + , Option [] ["timestep"] + (ReqArg (set configTimeStep . read) "FLOAT") + (describe configTimeStep "time step between simulation states") + + , Option [] ["epsilon"] + (ReqArg (set configEpsilon . read) "FLOAT") + (describe configEpsilon "smoothing parameter") + + , Option [] ["disc"] + (ReqArg (set configStartDiscSize . read) "FLOAT") + (describe configStartDiscSize "initial size of particle disc") + + , Option [] ["speed"] + (ReqArg (set configStartSpeed . read) "FLOAT") + (describe configStartSpeed "initial rotation speed of the disc") + + , Option [] ["max-steps"] + (ReqArg (set configMaxSteps . read) "INT") + (describe configMaxSteps "exit simulation after this many steps") + + , Option [] ["benchmark"] + (NoArg (set configBenchmark True)) + (describe configBenchmark "benchmark instead of displaying animation") + + , Option [] ["dump-final"] + (ReqArg (set configDumpFinal . Just) "FILE") + "dump final body positions to file" + + , Option ['h', '?'] ["help"] + (NoArg (set configHelp True)) + "show this help message" + ] + where + solver algorithm + | a `elem` ["n", "naive"] = Naive + | a `elem` ["bh", "barnshut", "barns-hut"] = BarnsHut + | otherwise = error $ "Unknown solver method: " ++ algorithm + where + a = map toLower algorithm + + describe f msg + = msg ++ " (" ++ show (get f defaults) ++ ")" + + +-- | Process the command line options +-- + +header :: [String] +header = + [ "accelerate-nbody (c) [2012..2013] The Accelerate Team" + , "" + , "Usage: accelerate-nbody [OPTIONS]" + , "" + ] + +footer :: [String] +footer = [] + diff --git a/examples/n-body/Gloss/Draw.hs b/examples/n-body/Gloss/Draw.hs new file mode 100644 --- /dev/null +++ b/examples/n-body/Gloss/Draw.hs @@ -0,0 +1,113 @@ + +-- +-- Drawing the world as a gloss picture. +-- +module Gloss.Draw ( draw ) + where + +import Config +import Common.Type +import Common.World +import Gloss.Simulate + +import Data.Label +import Graphics.Gloss +import qualified Data.Array.Accelerate as A + + +-- | Draw the simulation, optionally showing the Barnes-Hut tree. +-- +draw :: Config -> Simulate -> Picture +draw conf universe + = let + shouldDrawTree = get simulateDrawTree universe + world = get simulateWorld universe + + picPoints = Color (makeColor 1 1 1 0.4) + $ Pictures + $ map (drawBody conf) + $ A.toList + $ worldBodies world + + picTree = Blank +-- picTree = drawBHTree +-- $ L.buildTree +-- $ map massPointOfBody +-- $ V.toList +-- $ worldBodies world + + in Pictures [ if shouldDrawTree + then Color (makeColor 0.5 1.0 0.5 0.2) $ picTree + else Blank + + , picPoints ] + + +{-- +-- | Draw a list version Barnes-Hut tree. +drawBHTree :: L.BHTree -> Picture +drawBHTree bht + = drawBHTree' 0 bht + +drawBHTree' depth bht + = let + -- The bounding box + L.Box left down right up = L.bhTreeBox bht + [left', down', right', up'] = map realToFrac [left, down, right, up] + + picCell = lineLoop [(left', down'), (left', up'), (right', up'), (right', down')] + + + -- Draw a circle with an area equal to the mass of the centroid. + centroidX = realToFrac $ L.bhTreeCenterX bht + centroidY = realToFrac $ L.bhTreeCenterY bht + + centroidMass = L.bhTreeMass bht + circleRadius = realToFrac $ sqrt (centroidMass / pi) + + midX = (left' + right') / 2 + midY = (up' + down') / 2 + + picCentroid + | _:_ <- L.bhTreeBranch bht + , depth >= 1 + = Color (makeColor 0.5 0.5 1.0 0.4) + $ Pictures + [ Line [(midX, midY), (centroidX, centroidY)] + , Translate centroidX centroidY + $ ThickCircle + (circleRadius * 4 / 2) + (circleRadius * 4) ] + + | otherwise + = Blank + + -- The complete picture for this cell. + picHere = Pictures [picCentroid, picCell] + + -- Pictures of children. + picSubs = map (drawBHTree' (depth + 1)) + $ L.bhTreeBranch bht + + in Pictures (picHere : picSubs) +--} + + +-- | Draw a single body. Set the size of the body depending on it's mass, in +-- five size categories. +-- +drawBody :: Config -> Body -> Picture +drawBody conf ((position, mass), _, _) + = let sizeMax = get configBodyMass conf / 5 + size = 1 `max` mass / sizeMax + in + drawPoint position (size + 1) + + +-- | Draw a point using a filled circle. +-- +drawPoint :: Position -> R -> Picture +drawPoint (x, y, _) size + = Translate (realToFrac x) (realToFrac y) + $ ThickCircle (size / 2) size + diff --git a/examples/n-body/Gloss/Event.hs b/examples/n-body/Gloss/Event.hs new file mode 100644 --- /dev/null +++ b/examples/n-body/Gloss/Event.hs @@ -0,0 +1,27 @@ + +module Gloss.Event ( react ) + where + +import Gloss.Simulate + +import Data.Char +import Data.Label +import Graphics.Gloss.Interface.Pure.Game + + +-- React to user input events +-- +react :: Event -> Simulate -> Simulate +react event universe + = let + char ' ' Down = modify simulatePause not + char 't' Down = modify simulateDrawTree not + char _ _ = id + + special _ _ = id + + in case event of + EventKey (Char c) s _ _ -> char (toLower c) s universe + EventKey (SpecialKey c) s _ _ -> special c s universe + _ -> universe + diff --git a/examples/n-body/Gloss/Simulate.hs b/examples/n-body/Gloss/Simulate.hs new file mode 100644 --- /dev/null +++ b/examples/n-body/Gloss/Simulate.hs @@ -0,0 +1,48 @@ +{-# LANGUAGE TemplateHaskell #-} + +module Gloss.Simulate ( + + Simulate(..), simulateWorld, simulateDrawTree, simulatePause, + initialise, simulate, + +) where + + +import Common.Type +import Common.World + +import Data.Label + + +data Simulate = Simulate + { + _simulateWorld :: {-# UNPACK #-} !World + + -- parameters to control the visualisation + -- + , _simulateDrawTree :: !Bool + , _simulatePause :: !Bool + } + + +-- | Populate the universe with a set of bodies +-- +initialise :: World -> Simulate +initialise world + = Simulate + { + _simulateWorld = world + , _simulateDrawTree = False + , _simulatePause = False + } + +$(mkLabels [''Simulate]) + + +-- | Advance the simulation state +-- +simulate :: (Time -> World -> World) -> Float -> Simulate -> Simulate +simulate advance timeStep state + | _simulatePause state = state + | otherwise = modify simulateWorld (advance (realToFrac timeStep)) state + diff --git a/examples/n-body/Main.hs b/examples/n-body/Main.hs new file mode 100644 --- /dev/null +++ b/examples/n-body/Main.hs @@ -0,0 +1,87 @@ +-- +-- An N-Body simulation +-- + +-- friends +import Config +import ParseArgs +import Common.Body +import Common.World +import Gloss.Draw +import Gloss.Event +import Gloss.Simulate +import Random.Array +import Random.Position +import qualified Solver.Naive as Naive +import qualified Solver.BarnsHut as BarnsHut + +import Data.Array.Accelerate as A hiding ( size ) + +-- system +import Prelude as P +import Data.Label +import System.Environment +import Criterion.Main ( defaultMainWith, bench, whnf ) +import Graphics.Gloss.Interface.Pure.Game + + +main :: IO () +main + = do argv <- getArgs + (conf, cconf, nops) <- parseArgs configHelp configBackend options defaults header footer argv + + let solver = case get configSolver conf of + Naive -> Naive.calcAccels + BarnsHut -> BarnsHut.calcAccels + + n = get configBodyCount conf + size = get configWindowSize conf + fps = get configRate conf + epsilon = get configEpsilon conf + mass = get configBodyMass conf + radius = get configStartDiscSize conf + backend = get configBackend conf + + -- Generate random particle positions in a disc layout centred at + -- the origin. Start the system rotating with particle speed + -- proportional to distance from the origin + -- + positions = randomArray (cloud (size,size) radius) (Z :. n) + masses = randomArray (uniformR (1, mass)) (Z :. n) + + bodies = run backend + $ A.map (setStartVelOfBody . constant $ get configStartSpeed conf) + $ A.zipWith setMassOfBody (A.use masses) + $ A.map unitBody + $ A.use positions + + -- The initial simulation state + -- + universe = initialise world + world = World { worldBodies = bodies + , worldSteps = 0 + , worldTime = 0 } + + -- Advancing the simulation + -- + advance = advanceWorld step + step = P.curry + $ run1 backend + $ A.uncurry + $ advanceBodies (solver $ constant epsilon) + + -- Forward unto dawn + -- + if get configBenchmark conf + then withArgs nops $ defaultMainWith cconf (return ()) + [ bench "n-body" $ whnf (advance 0.1) world ] + + else play + (InWindow "N-Body" (size, size) (10, 10)) -- window size & position + black -- background colour + fps -- number of simulation steps per second + universe -- initial world + (draw conf) -- fn to convert a world into a picture + react -- fn to handle input events + (simulate advance) -- fn to advance the world + diff --git a/examples/n-body/Random/Position.hs b/examples/n-body/Random/Position.hs new file mode 100644 --- /dev/null +++ b/examples/n-body/Random/Position.hs @@ -0,0 +1,44 @@ +{-# LANGUAGE ViewPatterns #-} + +-- +-- Generators for various distributions of particle positions +-- +module Random.Position + where + +import Common.Type + +import Control.Monad.ST ( ST ) +import System.Random.MWC ( GenST, uniformR ) +import Data.Array.Accelerate.Array.Sugar as A + + +-- | Points distributed as a disc +-- +disc :: Position -> R -> sh -> GenST s -> ST s Position +disc (originX, originY, originZ) radiusMax _ix gen + = do radius <- uniformR (0,radiusMax) gen + theta <- uniformR (0, pi) gen + phi <- uniformR (0, 2*pi) gen + + return ( originX + radius * sin theta * cos phi + , originY + radius * sin theta * sin phi + , originZ + radius * cos theta ) + + +-- | A point cloud with areas of high and low density +-- +cloud :: Shape sh => (Int,Int) -> R -> sh -> GenST s -> ST s Position +cloud (fromIntegral -> sizeX, fromIntegral -> sizeY) radiusMax ix gen + = let + blob (sx,sy,sz) r + = disc (sx * sizeX, sy * sizeY, sz * (sizeX `min` sizeY)) + (radiusMax * r) + + in case A.size ix `mod` 5 of + 0 -> blob ( 0.25, 0.25, 0.25) 1.00 ix gen + 1 -> blob (-0.10, 0.10, 0.10) 0.60 ix gen + 2 -> blob (-0.05, 0.30,-0.30) 0.35 ix gen + 3 -> blob (-0.20,-0.12,-0.12) 0.45 ix gen + _ -> blob ( 0.15,-0.10, 0.20) 0.75 ix gen + diff --git a/examples/n-body/Solver/BarnsHut.hs b/examples/n-body/Solver/BarnsHut.hs new file mode 100644 --- /dev/null +++ b/examples/n-body/Solver/BarnsHut.hs @@ -0,0 +1,12 @@ + +module Solver.BarnsHut + where + +import Common.Type +import Data.Array.Accelerate as A + +-- | Calculate accelerations on the particles using the Barns-Hut algorithm +-- +calcAccels :: Exp R -> Acc (Vector Body) -> Acc (Vector Accel) +calcAccels = error "BarnsHut.calcAccels: not implemented yet!" + diff --git a/examples/n-body/Solver/Naive.hs b/examples/n-body/Solver/Naive.hs new file mode 100644 --- /dev/null +++ b/examples/n-body/Solver/Naive.hs @@ -0,0 +1,23 @@ + +module Solver.Naive + where + +import Common.Type +import Common.Body +import Common.Util + +import Data.Array.Accelerate as A + + +-- | Calculate accelerations on these particles in a naïve O(n^2) way +-- +calcAccels :: Exp R -> Acc (Vector Body) -> Acc (Vector Accel) +calcAccels epsilon bodies + = let n = A.size bodies + + cols = A.replicate (lift $ Z :. n :. All) bodies + rows = A.replicate (lift $ Z :. All :. n) bodies + + in + A.fold (.+.) (vec 0) $ A.zipWith (accel epsilon) rows cols + diff --git a/examples/quickcheck/Arbitrary/Array.hs b/examples/quickcheck/Arbitrary/Array.hs --- a/examples/quickcheck/Arbitrary/Array.hs +++ b/examples/quickcheck/Arbitrary/Array.hs @@ -1,6 +1,6 @@ +{-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE TypeOperators #-} module Arbitrary.Array where diff --git a/examples/quickcheck/Arbitrary/Shape.hs b/examples/quickcheck/Arbitrary/Shape.hs --- a/examples/quickcheck/Arbitrary/Shape.hs +++ b/examples/quickcheck/Arbitrary/Shape.hs @@ -1,6 +1,6 @@ +{-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE TypeOperators #-} module Arbitrary.Shape where diff --git a/examples/quickcheck/Config.hs b/examples/quickcheck/Config.hs --- a/examples/quickcheck/Config.hs +++ b/examples/quickcheck/Config.hs @@ -1,12 +1,12 @@ -{-# LANGUAGE TupleSections #-} +{-# LANGUAGE CPP #-} {-# LANGUAGE TemplateHaskell #-} -{-# LANGUAGE CPP #-} +{-# LANGUAGE TupleSections #-} module Config ( -- options & test configuration Options, - processArgs, + parseArgs, optBackend, double, float, int64, int32, int16, int8, -- running tests @@ -25,6 +25,7 @@ import qualified Data.Array.Accelerate.Interpreter as Interp #ifdef ACCELERATE_CUDA_BACKEND +import qualified Foreign.CUDA.Analysis as CUDA import qualified Foreign.CUDA.Driver as CUDA import qualified Data.Array.Accelerate.CUDA as CUDA #endif @@ -57,8 +58,8 @@ , _float = True , _int64 = True , _int32 = True - , _int16 = True - , _int8 = True + , _int16 = False + , _int8 = False } @@ -129,12 +130,12 @@ n <- CUDA.count devs <- mapM CUDA.device [0 .. n-1] props <- mapM CUDA.props devs - return $! opts { _double = any (\dev -> CUDA.computeCapability dev >= 1.3) props } + return $! opts { _double = any (\dev -> CUDA.computeCapability dev >= CUDA.Compute 1 3) props } #endif -processArgs :: [String] -> IO (Options, RunnerOptions) -processArgs argv = do +parseArgs :: [String] -> IO (Options, RunnerOptions) +parseArgs argv = do args <- interpretArgs argv (options, runner) <- case args of Left msg -> error msg diff --git a/examples/quickcheck/Main.hs b/examples/quickcheck/Main.hs --- a/examples/quickcheck/Main.hs +++ b/examples/quickcheck/Main.hs @@ -15,7 +15,7 @@ main = do -- process command line args, and print a brief usage message -- - (options, runner) <- processArgs =<< getArgs + (options, runner) <- parseArgs =<< getArgs -- the default execution order uses some knowledge of what functionality is -- required for each operation in the CUDA backend. diff --git a/examples/quickcheck/Test/Base.hs b/examples/quickcheck/Test/Base.hs --- a/examples/quickcheck/Test/Base.hs +++ b/examples/quickcheck/Test/Base.hs @@ -1,6 +1,6 @@ {-# LANGUAGE DefaultSignatures #-} {-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE TypeOperators #-} module Test.Base where @@ -58,16 +58,6 @@ epsilonRel = 0.001 epsilonAbs = 0.00001 - --- some "missing" instances --- -instance Eq Z where - Z == Z = True - -instance (Eq sh, Shape sh) => Eq (sh :. Int) where - (sh1 :. n1) == (sh2 :. n2) = n1 == n2 && sh1 == sh2 - (sh1 :. n1) /= (sh2 :. n2) = n1 /= n2 || sh1 /= sh2 - instance (Eq e, Eq sh, Shape sh) => Eq (Array sh e) where a1 == a2 = arrayShape a1 == arrayShape a2 && toList a1 == toList a2 @@ -116,6 +106,13 @@ splitEvery :: Int -> [a] -> [[a]] splitEvery _ [] = cycle [[]] -splitEvery n xs = let (h,t) = splitAt n xs - in h : splitEvery n t +splitEvery n xs = + let (h,t) = splitAt n xs + in h : splitEvery n t + +splitPlaces :: Integral i => [i] -> [a] -> [[a]] +splitPlaces [] _ = [] +splitPlaces (i:is) vs = + let (h,t) = splitAt (P.fromIntegral i) vs + in h : splitPlaces is t diff --git a/examples/quickcheck/Test/IndexSpace.hs b/examples/quickcheck/Test/IndexSpace.hs --- a/examples/quickcheck/Test/IndexSpace.hs +++ b/examples/quickcheck/Test/IndexSpace.hs @@ -1,5 +1,5 @@ -{-# LANGUAGE TypeOperators #-} {-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeOperators #-} module Test.IndexSpace where @@ -100,21 +100,12 @@ -- Reverse a vector -- - reverseRef xs = fromList (arrayShape xs) (reverse $ toList xs) - reverseAcc xs = - let xs' = use xs - n = unindex1 $ shape xs' - in - backpermute (shape xs') (\ix -> index1 $ n - unindex1 ix - 1) xs' + reverseAcc xs = Acc.reverse (use xs) + reverseRef xs = fromList (arrayShape xs) (P.reverse $ toList xs) -- Transpose a 2D matrix -- - transposeAcc xs = - let xs' = use xs - swap = lift1 $ \(Z:.x:.y) -> Z:.y:.x :: Z :. Exp Int :. Exp Int - in - backpermute (swap (shape xs')) swap xs' - + transposeAcc xs = Acc.transpose (use xs) transposeRef xs = let swap (Z:.x:.y) = Z :. y :. x in newArray (swap (arrayShape xs)) (\ix -> indexArray xs (swap ix)) diff --git a/examples/quickcheck/Test/Mapping.hs b/examples/quickcheck/Test/Mapping.hs --- a/examples/quickcheck/Test/Mapping.hs +++ b/examples/quickcheck/Test/Mapping.hs @@ -1,7 +1,7 @@ -{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE RankNTypes #-} {-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE TypeOperators #-} module Test.Mapping where diff --git a/examples/quickcheck/Test/PrefixSum.hs b/examples/quickcheck/Test/PrefixSum.hs --- a/examples/quickcheck/Test/PrefixSum.hs +++ b/examples/quickcheck/Test/PrefixSum.hs @@ -1,10 +1,10 @@ -{-# LANGUAGE RankNTypes #-} -{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE RankNTypes #-} {-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeOperators #-} module Test.PrefixSum where -import Prelude as P hiding ( sum ) +import Prelude as P import Test.QuickCheck import Data.Label import Data.Maybe @@ -14,7 +14,7 @@ import Config import Test.Base -import Arbitrary.Array () +import Arbitrary.Array import Data.Array.Accelerate as Acc @@ -40,32 +40,66 @@ testElt ok _ | P.not (get ok opt) = Nothing | otherwise = Just $ testGroup (show (typeOf (undefined :: e))) - [ testProperty "scanl" (test_scanl :: Vector e -> Property) - , testProperty "scanl'" (test_scanl' :: Vector e -> Property) - , testProperty "scanl1" (test_scanl1 :: Vector e -> Property) - , testProperty "scanr" (test_scanr :: Vector e -> Property) - , testProperty "scanr'" (test_scanr' :: Vector e -> Property) - , testProperty "scanr1" (test_scanr1 :: Vector e -> Property) + [ testProperty "scanl" (test_scanl :: Vector e -> Property) + , testProperty "scanl'" (test_scanl' :: Vector e -> Property) + , testProperty "scanl1" (test_scanl1 :: Vector e -> Property) + , testProperty "scanr" (test_scanr :: Vector e -> Property) + , testProperty "scanr'" (test_scanr' :: Vector e -> Property) + , testProperty "scanr1" (test_scanr1 :: Vector e -> Property) + -- + , testProperty "scanl1Seg" (test_scanl1seg :: Vector e -> Property) + , testProperty "scanr1Seg" (test_scanr1seg :: Vector e -> Property) + , testProperty "scanlSeg" (test_scanlseg :: Vector e -> Property) + , testProperty "scanrSeg" (test_scanrseg :: Vector e -> Property) + , testProperty "scanl'Seg" (test_scanl'seg :: Vector e -> Property) + , testProperty "scanr'Seg" (test_scanr'seg :: Vector e -> Property) ] -- left scan -- - test_scanl xs = run opt (Acc.scanl (+) 0 (use xs)) .==. scanlRef (+) 0 xs - test_scanl1 xs = run opt (Acc.scanl1 Acc.min (use xs)) .==. scanl1Ref P.min xs - test_scanl' xs = - let (vec, sum) = Acc.scanl' (+) 0 (use xs) - in (run opt vec, run opt sum) .==. scanl'Ref (+) 0 xs + test_scanl xs = run opt (Acc.scanl (+) 0 (use xs)) .==. scanlRef (+) 0 xs + test_scanl1 xs = run opt (Acc.scanl1 Acc.min (use xs)) .==. scanl1Ref P.min xs + test_scanl' xs = run opt (Acc.lift $ Acc.scanl' (+) 0 (use xs)) .==. scanl'Ref (+) 0 xs -- right scan -- - test_scanr xs = run opt (Acc.scanr (+) 0 (use xs)) .==. scanrRef (+) 0 xs - test_scanr1 xs = run opt (Acc.scanr1 Acc.max (use xs)) .==. scanr1Ref P.max xs - test_scanr' xs = - let (vec, sum) = Acc.scanr' (+) 0 (use xs) - in (run opt vec, run opt sum) .==. scanr'Ref (+) 0 xs + test_scanr xs = run opt (Acc.scanr (+) 0 (use xs)) .==. scanrRef (+) 0 xs + test_scanr1 xs = run opt (Acc.scanr1 Acc.max (use xs)) .==. scanr1Ref P.max xs + test_scanr' xs = run opt (Acc.lift $ Acc.scanr' (+) 0 (use xs)) .==. scanr'Ref (+) 0 xs + -- segmented left/right scan + -- + test_scanl1seg elt = + forAll arbitrarySegments1 $ \(seg :: Vector Int32) -> + forAll (arbitrarySegmentedArray seg) $ \xs -> + run opt (Acc.scanl1Seg (+) (use xs) (use seg)) .==. scanl1SegRef (+) (xs `asTypeOf` elt) seg + test_scanr1seg elt = + forAll arbitrarySegments1 $ \(seg :: Vector Int32) -> + forAll (arbitrarySegmentedArray seg) $ \xs -> + run opt (Acc.scanr1Seg (+) (use xs) (use seg)) .==. scanr1SegRef (+) (xs `asTypeOf` elt) seg + test_scanlseg elt = + forAll arbitrarySegments $ \(seg :: Vector Int32) -> + forAll (arbitrarySegmentedArray seg) $ \xs -> + run opt (Acc.scanlSeg (+) 0 (use xs) (use seg)) .==. scanlSegRef (+) 0 (xs `asTypeOf` elt) seg + + test_scanrseg elt = + forAll arbitrarySegments $ \(seg :: Vector Int32) -> + forAll (arbitrarySegmentedArray seg) $ \xs -> + run opt (Acc.scanrSeg (+) 0 (use xs) (use seg)) .==. scanrSegRef (+) 0 (xs `asTypeOf` elt) seg + + test_scanl'seg elt = + forAll arbitrarySegments $ \(seg :: Vector Int32) -> + forAll (arbitrarySegmentedArray seg) $ \xs -> + run opt (lift $ Acc.scanl'Seg (+) 0 (use xs) (use seg)) .==. scanl'SegRef (+) 0 (xs `asTypeOf` elt) seg + + test_scanr'seg elt = + forAll arbitrarySegments $ \(seg :: Vector Int32) -> + forAll (arbitrarySegmentedArray seg) $ \xs -> + run opt (lift $ Acc.scanr'Seg (+) 0 (use xs) (use seg)) .==. scanr'SegRef (+) 0 (xs `asTypeOf` elt) seg + + -- Reference implementation -- ------------------------ @@ -102,4 +136,57 @@ = Acc.fromList (arrayShape vec) . P.scanr1 f . Acc.toList $ vec + + +-- segmented operations +-- +scanlSegRef :: (Elt e, Integral i) => (e -> e -> e) -> e -> Vector e -> Vector i -> Vector e +scanlSegRef f z vec seg = + let seg' = toList seg + vec' = toList vec + n = P.sum $ P.map (\x -> P.fromIntegral x + 1) seg' + in fromList (Z :. n) $ + concat [ P.scanl f z v | v <- splitPlaces seg' vec' ] + +scanl1SegRef :: (Elt e, Integral i) => (e -> e -> e) -> Vector e -> Vector i -> Vector e +scanl1SegRef f vec seg = + let seg' = toList seg + vec' = toList vec + n = P.sum $ P.map P.fromIntegral seg' + in fromList (Z :. n) $ + concat [ P.scanl1 f v | v <- splitPlaces seg' vec' ] + +scanl'SegRef :: (Elt e, Integral i) => (e -> e -> e) -> e -> Vector e -> Vector i -> (Vector e, Vector e) +scanl'SegRef f z vec seg = + let seg' = toList seg + vec' = toList vec + scanl'_ v = let res = P.scanl f z v in (P.init res, P.last res) + (scans, sums) = P.unzip [ scanl'_ v | v <- splitPlaces seg' vec'] + in ( fromList (arrayShape vec) (concat scans) + , fromList (arrayShape seg) sums ) + +scanrSegRef :: (Elt e, Integral i) => (e -> e -> e) -> e -> Vector e -> Vector i -> Vector e +scanrSegRef f z vec seg = + let seg' = toList seg + vec' = toList vec + n = P.sum $ P.map (\x -> P.fromIntegral x + 1) seg' + in fromList (Z :. n) $ + concat [ P.scanr f z v | v <- splitPlaces seg' vec' ] + +scanr1SegRef :: (Elt e, Integral i) => (e -> e -> e) -> Vector e -> Vector i -> Vector e +scanr1SegRef f vec seg = + let seg' = toList seg + vec' = toList vec + n = P.sum $ P.map P.fromIntegral seg' + in fromList (Z :. n) $ + concat [ P.scanr1 f v | v <- splitPlaces seg' vec' ] + +scanr'SegRef :: (Elt e, Integral i) => (e -> e -> e) -> e -> Vector e -> Vector i -> (Vector e, Vector e) +scanr'SegRef f z vec seg = + let seg' = toList seg + vec' = toList vec + scanr'_ v = let res = P.scanr f z v in (P.tail res, P.head res) + (scans, sums) = P.unzip [ scanr'_ v | v <- splitPlaces seg' vec'] + in ( fromList (arrayShape vec) (concat scans) + , fromList (arrayShape seg) sums ) diff --git a/examples/quickcheck/Test/Reduction.hs b/examples/quickcheck/Test/Reduction.hs --- a/examples/quickcheck/Test/Reduction.hs +++ b/examples/quickcheck/Test/Reduction.hs @@ -1,6 +1,6 @@ -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeOperators #-} module Test.Reduction where @@ -16,7 +16,7 @@ import Config import Test.Base import Arbitrary.Array -import Data.Array.Accelerate as Acc +import Data.Array.Accelerate as Acc hiding (indexHead, indexTail) import Data.Array.Accelerate.Array.Sugar as Sugar @@ -207,7 +207,7 @@ (sh :. n) = arrayShape arr (Z :. sz) = arrayShape seg seg' = toList seg - foldseg xs = P.map (foldl' f z) (split seg' xs) + foldseg xs = P.map (foldl' f z) (splitPlaces seg' xs) fold1SegRef :: (Shape sh, Elt e, Elt i, Integral i) => (e -> e -> e) -> Array (sh :. Int) e -> Segments i -> Array (sh :. Int) e fold1SegRef f arr seg = fromList (sh :. sz) $ concat [ foldseg sub | sub <- splitEvery n (toList arr) ] @@ -215,11 +215,5 @@ (sh :. n) = arrayShape arr (Z :. sz) = arrayShape seg seg' = toList seg - foldseg xs = P.map (foldl1' f) (split seg' xs) - -split :: Integral i => [i] -> [a] -> [[a]] -split [] _ = [] -split (i:is) vs = - let (h,t) = splitAt (P.fromIntegral i) vs - in h : split is t + foldseg xs = P.map (foldl1' f) (splitPlaces seg' xs) diff --git a/examples/smoothlife/Config.hs b/examples/smoothlife/Config.hs new file mode 100644 --- /dev/null +++ b/examples/smoothlife/Config.hs @@ -0,0 +1,232 @@ +{-# LANGUAGE TemplateHaskell #-} + +module Config + where + +import Prelude as P hiding ((.), id, fst, snd) +import qualified Prelude as P + +import ParseArgs +import Data.Label +import Control.Category + +import Data.Array.Accelerate ( Array, DIM2 ) +import Data.Array.Accelerate.Math.Complex ( Complex ) + + +-- | Types +-- +type R = Float +type C = Complex R +type Matrix a = Array DIM2 a + +data SigmoidFunction + = Hard + | Smooth + | Atan + | Atancos + | Overshoot + | Linear + | Hermite + | Sin + deriving (Eq, Show, Read, Enum, Bounded) + + +data ColourScheme + = RedBlack + | WhiteBlack + | BlackWhite + | BrownGreen + | GoldBrown + | Rainbow1 + | Rainbow2 + | Rainbow3 + deriving (Eq, Show, Read, Enum, Bounded) + + +data Config = Config + { + -- How to execute the simulation + _configBackend :: Backend + + -- How to present the output + , _configWindowSize :: Int + , _configWindowZoom :: Int + , _configFramerate :: Int + , _configTimestep :: Float + , _configColourScheme :: ColourScheme + + -- Initial conditions + , _configRim :: R -- b + , _configDiscRadius :: (R,R) -- (ri, ra) + , _configBirthInterval :: (R,R) -- (b1, b2) + , _configDeathInterval :: (R,R) -- (d1, d2) + , _configStep :: (R,R) -- (alpha_n, alpha_m) + + -- Timestep and sigmoid functions + , _configTimestepMode :: Int + , _configSigmode :: Int + , _configSigtype :: SigmoidFunction + , _configMixtype :: SigmoidFunction + + -- Terminating conditions + , _configMaxSteps :: Maybe Int + , _configBenchmark :: Bool + , _configHelp :: Bool + } + +$(mkLabels [''Config]) + + +defaults :: Config +defaults = Config + { + _configBackend = maxBound + + , _configWindowSize = 256 + , _configWindowZoom = 3 + , _configFramerate = 30 + , _configTimestep = 0.1 + , _configColourScheme = WhiteBlack + + -- generic smooth glider + --, _configRim = 1 + --, _configDiscRadius = (4, 12) + --, _configBirthInterval = (0.278, 0.365) + --, _configDeathInterval = (0.267, 0.445) + --, _configStep = (0.028, 0.147) + + --, _configTimestepMode = 0 + --, _configSigmode = 2 + --, _configSigtype = Smooth + --, _configMixtype = Smooth + + -- SmoothLifeL + , _configRim = 1 + , _configDiscRadius = (10.0/3.0, 10) + , _configBirthInterval = (0.257, 0.336) + , _configDeathInterval = (0.365, 0.549) + , _configStep = (0.028, 0.147) + + , _configTimestepMode = 1 + , _configSigmode = 2 + , _configSigtype = Smooth + , _configMixtype = Smooth + + , _configMaxSteps = Nothing + , _configBenchmark = False + , _configHelp = False + } + + +options :: [OptDescr (Config -> Config)] +options = + [ Option [] ["size"] + (ReqArg (set configWindowSize . read) "INT") + (describe configWindowSize "visualisation size") + + , Option [] ["zoom"] + (ReqArg (set configWindowZoom . read) "INT") + (describe configWindowZoom "visualisation pixel replication factor") + + , Option [] ["framerate"] + (ReqArg (set configFramerate . read) "INT") + (describe configFramerate "visualisation frame rate") + + , Option [] ["timestep"] + (ReqArg (set configTimestep . read) "FLOAT") + (describe configTimestep "simulation timestep") + + , Option [] ["colours"] + (ReqArg (set configColourScheme . read) "COLOURS") + (describe configColourScheme "visualisation colour scheme") + + , Option [] ["b", "rim"] + (ReqArg (set configRim . read) "FLOAT") + (describe configRim "anti-aliasing zone around the rim") + + , Option [] ["ri", "inner-radius"] + (ReqArg (set (fst . configDiscRadius) . read) "FLOAT") + (describe (fst . configDiscRadius) "inner radius") + + , Option [] ["ra", "outer-radius"] + (ReqArg (set (snd . configDiscRadius) . read) "FLOAT") + (describe (snd . configDiscRadius) "outer radius") + + , Option [] ["b1", "birth-low"] + (ReqArg (set (fst . configBirthInterval) . read) "FLOAT") + (describe (fst . configBirthInterval) "lower birth interval") + + , Option [] ["b2", "birth-high"] + (ReqArg (set (snd . configBirthInterval) . read) "FLOAT") + (describe (snd . configBirthInterval) "upper birth interval") + + , Option [] ["d1", "death-low"] + (ReqArg (set (fst . configDeathInterval) . read) "FLOAT") + (describe (fst . configDeathInterval) "lower death interval") + + , Option [] ["d2", "death-high"] + (ReqArg (set (snd . configDeathInterval) . read) "FLOAT") + (describe (snd . configDeathInterval) "upper death interval") + + , Option [] ["sn", "alpha_n"] + (ReqArg (set (fst . configStep) . read) "FLOAT") + (describe (fst . configStep) "lower step interval") + + , Option [] ["sm", "alpha_m"] + (ReqArg (set (snd . configStep) . read) "FLOAT") + (describe (snd . configStep) "upper step interval") + + , Option [] ["tsm", "timestep-mode"] + (ReqArg (set configTimestepMode . read) "INT") + (describe configTimestepMode "timestep mode 0..2") + + , Option [] ["sigmode"] + (ReqArg (set configSigmode . read) "INT") + (describe configSigmode "sigmode 1..4") + + , Option [] ["sigtype"] + (ReqArg (set configSigtype . read) "SIGMOID") + (describe configSigtype "sigtype") + + , Option [] ["mixtype"] + (ReqArg (set configMixtype . read) "SIGMOID") + (describe configMixtype "mixtype") + + , Option [] ["max-steps"] + (ReqArg (set configMaxSteps . read) "INT") + (describe configMaxSteps "exit simulation after this many steps") + + , Option [] ["benchmark"] + (NoArg (set configBenchmark True)) + (describe configBenchmark "benchmark instead of displaying animation") + + , Option "h?" ["help"] + (NoArg (set configHelp True)) + "show this help message" + ] + where + describe f msg + = msg ++ " (" ++ show (get f defaults) ++ ")" + + fst = lens P.fst (\a (_,b) -> (a,b)) + snd = lens P.snd (\b (a,_) -> (a,b)) + +header :: [String] +header = + [ "accelerate-smoothlife (c) [2012..2013] The Accelerate Team" + , "" + , "Usage: accelerate-smoothlife [OPTIONS]" + , "" + ] + +footer :: [String] +footer = + [ "" + , "Available sigmoid functions:" + , " " ++ unwords (map show [minBound .. maxBound :: SigmoidFunction]) + , "" + , "Available colour schemes:" + , " " ++ unwords (map show [minBound .. maxBound :: ColourScheme]) + ] + diff --git a/examples/smoothlife/Gloss/Draw.hs b/examples/smoothlife/Gloss/Draw.hs new file mode 100644 --- /dev/null +++ b/examples/smoothlife/Gloss/Draw.hs @@ -0,0 +1,72 @@ + +module Gloss.Draw + where + +import Config + +import Prelude as P +import Data.Label +import Graphics.Gloss +import Data.Array.Accelerate as A hiding ( size ) +import Data.Array.Accelerate.IO as A +import Foreign.Ptr +import Foreign.ForeignPtr +import System.IO.Unsafe + +import Data.Array.Accelerate.Array.Data ( ptrsOfArrayData ) +import Data.Array.Accelerate.Array.Sugar ( Array(..) ) + + +colourise :: ColourScheme -> Acc (Matrix R) -> Acc (Matrix RGBA32) +colourise scheme = A.map (rgba32OfFloat . colour scheme) + where + phase = 0.01 + alpha = constant 1 + + fract x = x - A.fromIntegral (A.floor x :: Exp Int) + mix x y a = x*(1-a) + y*a + + colour RedBlack f = lift (f, constant 0, constant 0, alpha) + colour WhiteBlack f = lift (f, f, f, alpha) + colour BlackWhite f = lift (x, x, x, alpha) where x = 1-f + colour BrownGreen f = lift (mix 0.5 0.5 f, mix 0.3 0.75 f, mix 0 1 f, alpha) + colour GoldBrown f = + let ssf = sqrt (sqrt f) + in lift ( mix 0.5 (mix 1 0.3 f) ssf + , mix 0.3 (mix 0.95 0.2 f) ssf + , constant 0 + , alpha) + colour Rainbow1 f = rainbow (fract phase * 6) f + colour Rainbow2 f = rainbow (6 * sqrt (sqrt (1-f))) (sqrt (sqrt f)) + colour Rainbow3 f = rainbow' (sqrt (sqrt (1-f))) (sqrt (sqrt f)) + + rainbow p x + = p >* 0 &&* p <* 1 + ? (lift (x, x*p, constant 0, alpha) + , p <* 2 ? (lift (x*(2-p), x, constant 0, alpha) + , p <* 3 ? (lift (constant 0, x, x*(p-2), alpha) + , p <* 4 ? (lift (constant 0, x*(4-p), x, alpha) + , p <* 5 ? (lift (x*(p-4), constant 0, x, alpha) + , p <* 6 ? (lift (x, constant 0, x*(6-p), alpha) + , {-else-} (constant (0,0,0,0)))))))) + + rainbow' p x + = let c y = 0.5 * sin (1.7 * cos (pi * (p + y / 3.0 + phase))) + 0.5 + in lift (x * c 0, x * c 1, x * c 2, alpha) + + +draw :: Config -> Matrix RGBA32 -> Picture +draw conf (Array _ adata) = scale zoom zoom pic + where + zoom = P.fromIntegral + $ get configWindowZoom conf + size = get configWindowSize conf + + rawData = let ((), ptr) = ptrsOfArrayData adata + in unsafePerformIO $ newForeignPtr_ (castPtr ptr) + + pic = bitmapOfForeignPtr + size size -- raw image size + rawData -- the image data + False -- don't cache in texture memory + diff --git a/examples/smoothlife/Gloss/Event.hs b/examples/smoothlife/Gloss/Event.hs new file mode 100644 --- /dev/null +++ b/examples/smoothlife/Gloss/Event.hs @@ -0,0 +1,11 @@ + +module Gloss.Event + where + +import Gloss.Simulate +import Graphics.Gloss.Interface.Pure.Game + + +react :: Event -> Simulate -> Simulate +react _ = id + diff --git a/examples/smoothlife/Gloss/Simulate.hs b/examples/smoothlife/Gloss/Simulate.hs new file mode 100644 --- /dev/null +++ b/examples/smoothlife/Gloss/Simulate.hs @@ -0,0 +1,35 @@ +{-# LANGUAGE TemplateHaskell #-} + +module Gloss.Simulate ( + + Simulate(..), + +) where + +import Config +import Data.Label + + +data Mode = Hard | Linear | Hermite | Sine | Smooth | Atan | Atancos | Overshoot + deriving (Eq, Show) + +data Simulate = Simulate + { + -- Parameters to control the simulation + _sigmoidMode :: Mode + , _mixingMode :: Mode + } + +$(mkLabels [''Simulate]) + + +-- Initial simulation state +-- +initialise :: Config -> Simulate +initialise _ + = Simulate + { + _sigmoidMode = Smooth + , _mixingMode = Smooth + } + diff --git a/examples/smoothlife/Main.hs b/examples/smoothlife/Main.hs new file mode 100644 --- /dev/null +++ b/examples/smoothlife/Main.hs @@ -0,0 +1,69 @@ +{-# LANGUAGE CPP #-} +-- +-- A cellular automata simulation +-- + +-- friends +import Config +import ParseArgs +import SmoothLife +import Gloss.Draw +import Random.Array +import Random.Splat + +-- system +import Prelude as P +import Data.Array.Accelerate as A +import Data.Label +import Control.Exception +import System.Environment +import Graphics.Gloss +import Criterion.Main ( defaultMainWith, bench, whnf ) + + +main :: IO () +main + = do argv <- getArgs + (conf, cconf, nops) <- parseArgs configHelp configBackend options defaults header footer argv + + let -- visualisation configuration + n = get configWindowSize conf + zoom = get configWindowZoom conf + fps = get configFramerate conf + (ra,rb) = get configDiscRadius conf + + dish = Z :. n :. n + width = n * zoom + height = n * zoom + + backend = get configBackend conf + scheme = get configColourScheme conf + + render = draw conf + . run1 backend (colourise scheme) + advance = run1 backend (smoothlife conf) + + -- initialise with patches of random data + dots <- randomCircles dish ra rb + agar <- randomArrayWithSystemRandom (splat dots) dish + + world <- evaluate (advance agar) + + -- Rise minions! + if get configBenchmark conf + then withArgs nops $ defaultMainWith cconf (return ()) + [ bench "smoothlife" $ whnf advance world ] + +#ifndef ACCELERATE_ENABLE_GUI + else return () +#else + else play + (InWindow "Smooth Life" (width, height) (10, 20)) + black + fps + world + render + (\_ -> id) + (\_ -> advance) +#endif + diff --git a/examples/smoothlife/Random/Splat.hs b/examples/smoothlife/Random/Splat.hs new file mode 100644 --- /dev/null +++ b/examples/smoothlife/Random/Splat.hs @@ -0,0 +1,56 @@ + +module Random.Splat + where + +import Config + +import Control.Monad +import Control.Monad.ST +import System.Random.MWC +import Data.Array.Accelerate ( DIM2, Z(..), (:.)(..) ) + + + +-- Generate a random point on the grid by checking which circles it overlaps. +-- +splat :: [(DIM2, R, Bool)] -> DIM2 -> GenST s -> ST s R +splat circles ix _ = + let overlaps = filter (inside ix) circles + negative (_,_,s) = s + in + return $ case overlaps of + [] -> 0 + xs -> if negative (last xs) then 0 else 1 + + +-- Generate some random circles with centre position somewhere inside the grid +-- and some maximum radius. The boolean represents whether life should be +-- created or destroyed inside this region. +-- +randomCircles :: DIM2 -> R -> R -> IO [(DIM2, R, Bool)] +randomCircles (Z :. height :. width) radiusMin radiusMax = + withSystemRandom . asGenIO $ \gen -> + let + circle = do + a <- uniform gen + r <- uniformR (radiusMin, radiusMax) gen + y <- uniformR (0, height-1) gen + x <- uniformR (0, width-1) gen + return $ (Z:.y:.x, r, a) + + n = 1 `max` round ( fromIntegral (width * height) / (radiusMax * radiusMax) ) + in + replicateM n circle + + +-- Test if a point is inside the given circle, given by centre coordinate and +-- radius. +-- +inside :: DIM2 -> (DIM2, R, Bool) -> Bool +inside (Z:.y:.x) (Z:.cy:.cx, r, _) = + let + dx = fromIntegral $ cx - x + dy = fromIntegral $ cy - y + in + dx*dx + dy*dx <= r*r + diff --git a/examples/smoothlife/SmoothLife.hs b/examples/smoothlife/SmoothLife.hs new file mode 100644 --- /dev/null +++ b/examples/smoothlife/SmoothLife.hs @@ -0,0 +1,160 @@ +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeOperators #-} +-- +-- A cellular automata simulation over a smooth domain. +-- +-- http://arxiv.org/abs/1111.1567 +-- http://sourceforge.net/projects/smoothlife/ +-- + +module SmoothLife + where + +import Config +import ParseArgs + +import Prelude as P +import Data.Label +import Data.Array.Accelerate as A hiding ( size ) +import Data.Array.Accelerate.Math.FFT +import Data.Array.Accelerate.Math.DFT.Centre +import Data.Array.Accelerate.Math.Complex + +import Data.Array.Accelerate.CUDA.Foreign + +-- Smooth life +-- ~~~~~~~~~~~ + +smoothlife + :: Config + -> Acc (Matrix R) + -> Acc (Matrix R) +smoothlife conf aa + = aa'' + where + -- A simulation step + -- + aaf = fft2D' Forward size size (complex aa) + nf = A.zipWith (*) aaf (use krf') + mf = A.zipWith (*) aaf (use kdf') + n = A.map (\x -> real x / kflr'') (fft2D' Inverse size size nf) + m = A.map (\x -> real x / kfld'') (fft2D' Inverse size size mf) + aa' = snm conf sn sm b1 b2 d1 d2 n m + aa'' = clamp $ A.zipWith timestepMode aa' aa + + -- simulation parameters + -- + b = get1 configRim conf + (ri,ra) = get2 configDiscRadius conf + (b1,b2) = get2 configBirthInterval conf + (d1,d2) = get2 configDeathInterval conf + (sn,sm) = get2 configStep conf -- aka. alpha_n alpha_m + dt = get1 configTimestep conf + + timestepMode f g = timestepModes f g P.!! get configTimestepMode conf + + size = get configWindowSize conf + sh = constant (Z:.size:.size) + + -- initial state + -- + kflr = A.sum kr + kfld = A.sum kd + krf = fft2D' Forward size size (shift2D (complex kr)) + kdf = fft2D' Forward size size (shift2D (complex kd)) + + kd = A.generate sh (\ix -> 1 - linear (radius ix) ri b) + kr = A.generate sh (\ix -> let r = radius ix + in linear r ri b * (1 - linear r ra b)) + + kflr'' = constant (kflr' `A.indexArray` Z) + kfld'' = constant (kfld' `A.indexArray` Z) + (kflr', kfld', krf', kdf') + = run (get configBackend conf) $ lift (kflr, kfld, krf, kdf) + + -- Auxiliary + -- + get1 f c = constant $ get f c + get2 f c = let (x,y) = get f c in (constant x, constant y) + + complex = A.map (\x -> lift (x, constant 0)) + + radius ix = + let Z:.y':.x' = unlift ix :: Z :. Exp Int :. Exp Int + x = A.fromIntegral $ x' - constant (size `div` 2) + y = A.fromIntegral $ y' - constant (size `div` 2) + in + sqrt (x*x + y*y) + + linear x l u + = x <* l-u/2 ? ( 0.0 + , x >* l+u/2 ? ( 1.0 + , (x - l + u / 2) / u )) + + clamp = A.map + (\x -> A.min (A.max x 0.0) 1.0) + + timestepModes f g + = [ f + , g + dt*(2.0*f-1.0) + , g + dt*(f-g) ] + + +-- Equation 6: s(n,m) +-- ~~~~~~~~~~~~~~~~~~ + +-- Also a few additional modes as discovered from the source code +-- +snm :: Config -> Exp R -> Exp R -> Exp R -> Exp R -> Exp R -> Exp R -> Acc (Matrix R) -> Acc (Matrix R) -> Acc (Matrix R) +snm conf sn sm b1 b2 d1 d2 + = A.zipWith sigmode + where + sigtype = getSigmoidFunction (get configSigtype conf) + mixtype = getSigmoidFunction (get configMixtype conf) + sigmode n m = sigmodes n m P.!! (get configSigmode conf - 1) + + sigmoid_ab :: Exp R -> Exp R -> Exp R -> Exp R + sigmoid_ab x a b + = sigtype x a sn * (1.0 - sigtype x b sn) + + sigmoid_mix :: Exp R -> Exp R -> Exp R -> Exp R + sigmoid_mix x y m + = x * (1 - mixtype m 0.5 sm) + y * mixtype m 0.5 sm + + mix :: Exp R -> Exp R -> Exp R -> Exp R + mix x y a = x * (1 - a) + y*a + + -- available sigmodes + sigmodes n m + = [ mix (sigmoid_ab n b1 b2) + (sigmoid_ab n d1 d2) m + , sigmoid_mix (sigmoid_ab n b1 b2) + (sigmoid_ab n d1 d2) m + , sigmoid_ab n (mix b1 d1 m) + (mix b2 d2 m) + , sigmoid_ab n (sigmoid_mix b1 d1 m) + (sigmoid_mix b2 d2 m)] + + +getSigmoidFunction :: SigmoidFunction -> Exp R -> Exp R -> Exp R -> Exp R +getSigmoidFunction f x a ea + = let + -- __expf is CUDA's faster but less precise version of exp. + cexp = foreignExp (cudaExp "math_functions.h __expf") exp + in + case f of + Hard -> x >=* a ? (1, 0) + Smooth -> 1.0/(1.0+cexp(-(x-a)*4.0/ea)) + Atan -> atan ((x-a) * pi/ea) / pi + 0.5 + Atancos -> 0.5 * (0.5 * atan ((x-a) / ea) / pi * cos ((x-a) * 1.4) * 1.1 + 1.0) + Overshoot -> 0.5 + (1.0 / (1.0 + exp (-(x-a)*4.0/ea)) - 0.5) * (1.0 + exp(-(x-a)*(x-a)/ea/ea)) + Linear -> bounded (\x' a' ea' -> (x'-a')/ea' + 0.5) x a ea + Hermite -> bounded (\x' a' ea' -> let v = (x' - (a'-ea'/2.0))/ea' in v * v * (3.0-2.0*v)) x a ea + Sin -> bounded (\x' a' ea' -> sin (pi * (x'-a')/ea') * 0.5 + 0.5) x a ea + + where + bounded :: (Exp R -> Exp R -> Exp R -> Exp R) -> Exp R -> Exp R -> Exp R -> Exp R + bounded f' x' a' ea' + = x' <* a'-ea'/2.0 ? ( 0.0 + , x' >* a'+ea'/2.0 ? ( 1.0 + , f' x' a' ea' )) diff --git a/examples/tests/Benchmark.hs b/examples/tests/Benchmark.hs --- a/examples/tests/Benchmark.hs +++ b/examples/tests/Benchmark.hs @@ -1,3 +1,4 @@ +{-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleContexts #-} {-# OPTIONS_GHC -fno-warn-orphans #-} @@ -9,7 +10,6 @@ import Data.List import Data.Monoid import Criterion -import Criterion.IO import Criterion.Config import Criterion.Main import Criterion.Monad @@ -21,6 +21,12 @@ import System.IO import System.Directory import System.Environment + +#if MIN_VERSION_criterion(0,8,0) +import Criterion.IO.Printf +#else +import Criterion.IO +#endif instance (Ix dim, IArray UArray e) => NFData (UArray dim e) where diff --git a/examples/tests/Config.hs b/examples/tests/Config.hs --- a/examples/tests/Config.hs +++ b/examples/tests/Config.hs @@ -79,6 +79,7 @@ -- criterion hooks , cfgPerformGC :: Bool , cfgConfidence :: Maybe Double + , cfgSamples :: Maybe Int , cfgResamples :: Maybe Int , cfgSummaryFile :: Maybe FilePath @@ -161,9 +162,14 @@ , cfgResamples = def &= explicit - &= name "s" &= name "resamples" &= help "Number of bootstrap resamples to perform" + + , cfgSamples = def + &= explicit + &= name "s" + &= name "samples" + &= help "Number of samples to collect" , cfgSummaryFile = def &= name "u" diff --git a/examples/tests/Main.hs b/examples/tests/Main.hs --- a/examples/tests/Main.hs +++ b/examples/tests/Main.hs @@ -49,6 +49,7 @@ unless (verbose == Quiet) $ putStrLn "" let args = [ maybe "" (\ci -> "--ci=" ++ show ci) (cfgConfidence cfg) , maybe "" (\r -> "--resamples=" ++ show r) (cfgResamples cfg) + , maybe "" (\r -> "--samples=" ++ show r) (cfgSamples cfg) , maybe "" (\f -> "--summary=" ++ f) (cfgSummaryFile cfg) , if cfgPerformGC cfg then "-g" else "-G" , case verbose of diff --git a/examples/tests/PGM.hs b/examples/tests/PGM.hs new file mode 100644 --- /dev/null +++ b/examples/tests/PGM.hs @@ -0,0 +1,32 @@ +-- +-- Load a PGM file. MacOS X users might find the following quicklook plugin +-- useful for viewing PGM files: +-- +-- http://code.google.com/p/quicklook-pfm/ +-- + +module PGM where + +import Control.Applicative +import Graphics.Pgm + +import Prelude as P +import Data.Array.Accelerate as Acc +import Data.Array.Unboxed hiding (Array) +import qualified Data.ByteString as B + + +-- Read an 8-bit PGM file, and marshal to an Accelerate array as floating-point +-- data in the range [0,1]. +-- +readPGM :: FilePath -> IO (Array DIM2 Float) +readPGM fp = do + img <- either (error . show) head . pgmsToArrays <$> B.readFile fp :: IO (UArray (Int,Int) Word8) + return . fromIArray $ amap (\x -> P.fromIntegral x / 255) img + + +writePGM :: FilePath -> Array DIM2 Float -> IO () +writePGM fp img = + let arr = toIArray img :: UArray (Int,Int) Float + in arrayToFile fp $ amap (\x -> P.round (255 * x)) arr + diff --git a/examples/tests/Test.hs b/examples/tests/Test.hs --- a/examples/tests/Test.hs +++ b/examples/tests/Test.hs @@ -31,10 +31,8 @@ import qualified Radix import qualified SliceExamples -#ifdef ACCELERATE_IO import qualified BlockCopy -import qualified VectorCopy () -- FIXME: Need to add tests from here -#endif +import qualified VectorCopy import qualified Canny import qualified IntegralImage @@ -139,17 +137,16 @@ , mkTest "permute-hist" "histogram" $ Permute.run "histogram" n , mkTest "backpermute-reverse" "reverse a vector" $ Backpermute.run "reverse" n , mkTest "backpermute-transpose" "transpose a matrix" $ Backpermute.run2d "transpose" n - , mkTest "init" "vector init" $ Vector.run "init" n - , mkTest "tail" "vector tail" $ Vector.run "tail" n - , mkTest "take" "vector take" $ Vector.run "take" n - , mkTest "drop" "vector drop" $ Vector.run "drop" n - , mkTest "slit" "vector slit" $ Vector.run "slit" n - , mkTest "gather" "backpermute via index mapping vector" $ Gather.run "gather" n - , mkTest "gather-if" "cond. backpermute via index mapping vector" $ Gather.run "gather-if" n - , mkTest "scatter" "permute via index mapping vector" $ Scatter.run "scatter" n - , mkTest "scatter-if" "cond. permute via index mapping vector" $ Scatter.run "scatter-if" n + , mkTest "init" "vector init" $ Vector.run "init" n + , mkTest "tail" "vector tail" $ Vector.run "tail" n + , mkTest "take" "vector take" $ Vector.run "take" n + , mkTest "drop" "vector drop" $ Vector.run "drop" n + , mkTest "slit" "vector slit" $ Vector.run "slit" n + , mkTest "gather" "backpermute via index mapping vector" $ Gather.run "gather" n + , mkTest "gather-if" "conditional backwards index mapping" $ Gather.run "gather-if" n + , mkTest "scatter" "permute via index mapping vector" $ Scatter.run "scatter" n + , mkTest "scatter-if" "conditional permutation via index mapping" $ Scatter.run "scatter-if" n - -- simple examples , mkTest "sasum" "sum of absolute values" $ SASUM.run n , mkTest "saxpy" "scalar alpha*x + y" $ SAXPY.run n @@ -159,19 +156,18 @@ , mkTest "black-scholes" "Black-Scholes option pricing" $ BlackScholes.run n , mkTest "radixsort" "radix sort" $ Radix.run n -#ifdef ACCELERATE_IO -- Array IO - , mkIO "io" "array IO test" $ BlockCopy.run -#endif + , mkIO "io" "array ptr copy test" $ BlockCopy.run + , mkIO "io" "array vector copy test" $ VectorCopy.run -- image processing - , mkNoRef "canny" "canny edge detection" $ Canny.run img - , mkNoRef "integral-image" "image integral (2D scan)" $ IntegralImage.run img + , mkIO "canny" "canny edge detection" $ Canny.canny (backend cfg) img 50 100 + , mkNoRef "integral-image" "image integral (2D scan)" $ IntegralImage.run img -- slices , mkTest "slices" "replicate (Z:.2:.All:.All)" $ SliceExamples.run1 , mkTest "slices" "replicate (Z:.All:.2:.All)" $ SliceExamples.run2 , mkTest "slices" "replicate (Z:.All:.All:.2)" $ SliceExamples.run3 - , mkTest "slices" "replicate (Any:.2)" $ SliceExamples.run4 + , mkTest "slices" "replicate (Any:.2)" $ SliceExamples.run4 , mkTest "slices" "replicate (Z:.2:.2:.2)" $ SliceExamples.run5 -- , mkIO "sharing-recovery" "simple" $ return (show SharingRecovery.simple) diff --git a/examples/tests/image-processing/Canny.hs b/examples/tests/image-processing/Canny.hs --- a/examples/tests/image-processing/Canny.hs +++ b/examples/tests/image-processing/Canny.hs @@ -1,117 +1,350 @@ +{-# LANGUAGE Rank2Types #-} +{-# LANGUAGE BangPatterns #-} {-# LANGUAGE TypeOperators #-} -module Canny where +module Canny (canny) where -import PGM +import Prelude as P -import Data.Array.Accelerate hiding (zipWith, unindex2) -import qualified Data.Array.Accelerate as Acc +import Data.Array.Accelerate as A +import Data.Array.Accelerate.IO as A +import Data.Array.Repa.Repr.Unboxed ( U ) +import qualified Data.Array.Repa as R +import qualified Data.Vector.Unboxed as V +import qualified Data.Vector.Unboxed.Mutable as VM +-- import qualified Data.Array.Repa.IO.BMP as R +-- import qualified Data.Array.Repa.Repr.Unboxed as R -type Image a = Array DIM2 a -type Stencil7x1 a = (Stencil3 a, Stencil7 a, Stencil3 a) -type Stencil1x7 a = (Stencil3 a, Stencil3 a, Stencil3 a, Stencil3 a, Stencil3 a, Stencil3 a, Stencil3 a) +canny :: (forall a. Arrays a => Acc a -> a) + -> FilePath + -> Float + -> Float + -> IO (R.Array U R.DIM2 Word8) +canny run fileIn threshLow threshHigh = do -convolve7x1 :: (Elt a, IsNum a) => [Exp a] -> Stencil7x1 a -> Exp a -convolve7x1 kernel (_, (a,b,c,d,e,f,g), _) = - sum $ zipWith (*) kernel [a,b,c,d,e,f,g] + -- Read in the image BMP file + -- + img <- either (error . show) id `fmap` readImageFromBMP fileIn -convolve1x7 :: (Elt a, IsNum a) => [Exp a] -> Stencil1x7 a -> Exp a -convolve1x7 kernel ((_,a,_), (_,b,_), (_,c,_), (_,d,_), (_,e,_), (_,f,_), (_,g,_)) = - sum $ zipWith (*) kernel [a,b,c,d,e,f,g] + -- Setup and run the Accelerate kernel stages, which identify strong and weak + -- edges in the image. + -- + let low = unit (constant threshLow) + high = unit (constant threshHigh) + grey = toGreyscale + blurred = gaussianY . gaussianX . grey + magdir = gradientMagDir low . blurred + suppress = nonMaximumSuppression low high . magdir --- Gaussian smoothing + (image,strong) = run $ stage1 (use img) + stage1 x = let img' = suppress x + in lift $ (img', selectStrong img') + + -- Now use Repa to trace out weak edges connected to strong edges + -- + edges <- wildfire (A.toRepa image) (A.toRepa strong) + +-- R.writeImageToBMP "edges.bmp" (R.zip3 edges edges edges) + return edges + +{-- + -- Benchmark each (accelerate) stage. The "kernels" group is intended to be + -- run under nvvp, whereas the canny group is for end-to-end benchmarks. + -- + (args,rest) <- break (== "--") `fmap` getArgs + let opts = if null rest then [] else P.tail rest + force2 a = A.indexArray a (Z:.0:.0 :: DIM2) `seq` () + force1 a = A.indexArray a (Z:.0 :: DIM1) `seq` () + + -- Need to force partial results so that we benchmark individual stages + -- + grey' = run $ toGreyscale (use img) + blurX' = run $ gaussianX (use grey') + blurred' = run $ gaussianY (use blurX') + gradX' = run $ gradientX (use blurred') + gradY' = run $ gradientY (use blurred') + magdir' = run $ gradientMagDir low (use blurred') + suppress' = run $ nonMaximumSuppression low high (use magdir') + + withArgs opts $ defaultMain + [ bgroup "kernels" + [ bench "greyscale" $ whnf (force2 . run . grey) (use img) + , bench "blur-x" $ whnf (force2 . run . gaussianX) (use grey') + , bench "blur-y" $ whnf (force2 . run . gaussianY) (use blurX') + , bench "grad-x" $ whnf (force2 . run . gradientX) (use blurred') + , bench "grad-y" $ whnf (force2 . run . gradientY) (use blurred') + , bench "mag-orient" $ whnf (force2 . run . gradientMagDir low) (use blurred') + , bench "suppress" $ whnf (force2 . run . nonMaximumSuppression low high) (use magdir') + , bench "strong" $ whnf (force1 . run . selectStrong) (use suppress') + ] + + , bgroup "canny" + [ bench "run" $ whnf (force1 . P.snd . run . stage1) (use img) + , bench "run1" $ whnf (force1 . P.snd . run1 stage1) img + ] + ] +--} + +-- Accelerate component -------------------------------------------------------- + +type RGBA = Word32 +type Image a = Array DIM2 a + +type Stencil5x1 a = (Stencil3 a, Stencil5 a, Stencil3 a) +type Stencil1x5 a = (Stencil3 a, Stencil3 a, Stencil3 a, Stencil3 a, Stencil3 a) + +-- Classification of the output pixel +data Orient = Undef | PosD | Vert | NegD | Horiz +data Edge = None | Weak | Strong + +orient :: Orient -> Int +orient Undef = 0 +orient PosD = 64 +orient Vert = 128 +orient NegD = 192 +orient Horiz = 255 + +orient' :: Orient -> Exp Int +orient' = constant . orient + +edge :: Edge -> Float +edge None = 0 +edge Weak = 0.5 +edge Strong = 1.0 + +edge' :: Edge -> Exp Float +edge' = constant . edge + +convolve5x1 :: (Elt a, IsNum a) => [Exp a] -> Stencil5x1 a -> Exp a +convolve5x1 kernel (_, (a,b,c,d,e), _) + = P.foldl1 (+) + $ P.zipWith (*) kernel [a,b,c,d,e] + +convolve1x5 :: (Elt a, IsNum a) => [Exp a] -> Stencil1x5 a -> Exp a +convolve1x5 kernel ((_,a,_), (_,b,_), (_,c,_), (_,d,_), (_,e,_)) + = P.foldl1 (+) + $ P.zipWith (*) kernel [a,b,c,d,e] + + +-- RGB to Greyscale conversion, in the range [0,255] -- -gaussian :: (Elt a, IsFloating a) => [Exp a] -gaussian = [ 0.00442012927963 - , 0.05384819825462 - , 0.24133088157513 - , 0.39788735772974 - , 0.24133088157513 - , 0.05384819825462 - , 0.00442012927963 ] +toGreyscale :: Acc (Image RGBA) -> Acc (Image Float) +toGreyscale = A.map (\rgba -> 255 * luminanceOfRGBA32 rgba) -gaussianX :: (Elt a, IsFloating a) => Acc (Image a) -> Acc (Image a) -gaussianX = stencil (convolve7x1 gaussian) (Constant 0) -gaussianY :: (Elt a, IsFloating a) => Acc (Image a) -> Acc (Image a) -gaussianY = stencil (convolve1x7 gaussian) (Constant 0) +-- Separable Gaussian blur in the x- and y-directions +-- +gaussianX :: Acc (Image Float) -> Acc (Image Float) +gaussianX = stencil (convolve5x1 gaussian) Clamp + where + gaussian = [ 1, 4, 6, 4, 1 ] +gaussianY :: Acc (Image Float) -> Acc (Image Float) +gaussianY = stencil (convolve1x5 gaussian) Clamp + where + gaussian = P.map (/256) [ 1, 4, 6, 4, 1 ] --- Gaussian derivative and gradient quantisation + +-- Gradients in the x- and y- directions -- -gaussian' :: (Elt a, IsFloating a) => [Exp a] -gaussian' = [ 0.02121662054222 - , 0.17231423441479 - , 0.38612941052022 - , 0.0 - ,-0.38612941052022 - ,-0.17231423441479 - ,-0.02121662054222 ] +gradientX :: Acc (Image Float) -> Acc (Image Float) +gradientX = stencil grad Clamp + where + grad :: Stencil3x3 Float -> Exp Float + grad ((u, _, x) + ,(v, _, y) + ,(w, _, z)) = x + (2*y) + z - u - (2*v) - w -gradientX :: (Elt a, IsFloating a) => Acc (Image a) -> Acc (Image a) -gradientX = stencil (convolve7x1 gaussian') (Constant 0) +gradientY :: Acc (Image Float) -> Acc (Image Float) +gradientY = stencil grad Clamp + where + grad :: Stencil3x3 Float -> Exp Float + grad ((x, y, z) + ,(_, _, _) + ,(u, v, w)) = x + (2*y) + z - u - (2*v) - w -gradientY :: (Elt a, IsFloating a) => Acc (Image a) -> Acc (Image a) -gradientY = stencil (convolve1x7 gaussian') (Constant 0) -gradientMagnitude :: (Elt a, IsFloating a) => Acc (Image a) -> Acc (Image a) -> Acc (Image a) -gradientMagnitude = Acc.zipWith magdir +-- Classify the magnitude and orientation of the image gradient. +-- +-- Because accelerate supports generalised stencil functions, not just +-- convolutions, we can combine the x- and y- sobel operators and save some +-- memory bandwidth. +-- +gradientMagDir + :: Acc (Scalar Float) + -> Acc (Image Float) + -> Acc (Array DIM2 (Float,Int)) +gradientMagDir threshLow = stencil magdir Clamp where - magdir dx dy = let mag = sqrt (dx*dx + dy*dy) - -- dir = atan2 dy dx - in mag -- lift (mag, dir) + magdir :: Stencil3x3 Float -> Exp (Float,Int) + magdir ((v0, v1, v2) + ,(v3, _, v4) + ,(v5, v6, v7)) = + let + -- Image gradients + dx = v2 + (2*v4) + v7 - v0 - (2*v3) - v5 + dy = v0 + (2*v1) + v2 - v5 - (2*v6) - v7 + -- Magnitude + mag = sqrt (dx * dx + dy + dy) --- Non-maximum suppression + -- Direction + -- + -- Determine the angle of the vector and rotate it around a bit to + -- make the segments easier to classify + theta = atan2 dy dx + alpha = (theta - (pi/8)) * (4/pi) + + -- Normalise the angle to between [0..8) + norm = alpha + 8 * A.fromIntegral (boolToInt (alpha <=* 0)) + + -- Try to avoid doing explicit tests, to avoid warp divergence + low = the threshLow + undef = abs dx <=* low &&* abs dy <=* low + dir = boolToInt (A.not undef) * ((64 * (1 + A.floor norm `mod` 4)) `A.min` 255) + in + lift (mag, dir) + + +-- Non-maximum suppression classifies pixels that are the local maximum along +-- the direction of the image gradient as either strong or weak edges. All other +-- pixels are not considered edges at all. -- +-- The image intensity is in the range [0,1] +-- nonMaximumSuppression - :: (Elt a, IsFloating a) - => Exp a - -> Acc (Image a) - -> Acc (Image a) - -> Acc (Image a) - -> Acc (Image a) -nonMaximumSuppression threshold gradX gradY gradM = - generate (shape gradX) $ \ix -> - let dx = gradX ! ix - dy = gradY ! ix - mag = gradM ! ix - alpha = 1.3065629648763766 -- 0.5 / sin (pi / 8.0) - offsetx = Acc.round (alpha * dx / mag) - offsety = Acc.round (alpha * dy / mag) + :: Acc (Scalar Float) + -> Acc (Scalar Float) + -> Acc (Image (Float,Int)) + -> Acc (Image Float) +nonMaximumSuppression threshLow threshHigh magdir = + generate (shape magdir) $ \ix -> + let -- The input parameters -- - (m,n) = unindex2 (shape gradX) - (x,y) = unindex2 ix - fwd = gradM ! clamp (x+offsetx, y+offsety) - rev = gradM ! clamp (x-offsetx, y-offsety) + low = the threshLow + high = the threshHigh + (mag, dir) = unlift (magdir ! ix) + Z :. h :. w = unlift (shape magdir) + Z :. y :. x = unlift ix + + -- Determine the points that lie either side of this point along to the + -- direction of the image gradient. -- - unindex2 uv = let Z:.u:.v = unlift uv in (u,v) - clamp (u,v) = lift (Z:. 0 `Acc.max` u `Acc.min` (m-1) :. 0 `Acc.max` v `Acc.min` (n-1)) + -- The direction coding: + -- + -- 192 128 64 + -- | + -- 255 --- --- + -- + offsetx = dir >* orient' Vert ? (-1, dir <* orient' Vert ? (1, 0)) + offsety = dir <* orient' Horiz ? (-1, 0) + + (fwd, _) = unlift $ magdir ! lift (clamp (Z :. y+offsety :. x+offsetx)) :: (Exp Float, Exp Int) + (rev, _) = unlift $ magdir ! lift (clamp (Z :. y-offsety :. x-offsetx)) :: (Exp Float, Exp Int) + + clamp (Z:.u:.v) = Z :. 0 `A.max` u `A.min` (h-1) :. 0 `A.max` v `A.min` (w-1) + + -- Try to avoid doing explicit tests to avoid warp divergence. + -- + none = dir ==* orient' Undef ||* mag <* low ||* mag <* fwd ||* mag <* rev + strong = mag >=* high in - (mag <* threshold ||* fwd >* mag ||* rev >* mag) ? (0, 1) + A.fromIntegral (boolToInt (A.not none) * (1 + boolToInt strong)) * 0.5 --- Canny edge detection +-- Extract the linear indices of the strong edges -- -canny :: (Elt a, IsFloating a) => Image a -> Acc (Image a) -canny img = - let smooth = gaussianX . gaussianY $ use img - gradX = gradientX smooth - gradY = gradientY smooth - gradMag = gradientMagnitude gradX gradY +selectStrong + :: Acc (Image Float) + -> Acc (Array DIM1 Int) +selectStrong img = + let strong = A.map (\x -> boolToInt (x ==* edge' Strong)) (flatten img) + (targetIdx, len) = A.scanl' (+) 0 strong + indices = A.enumFromN (index1 $ size img) 0 + zeros = A.fill (index1 $ the len) 0 in - nonMaximumSuppression 0.1 gradX gradY gradMag + A.permute const zeros (\ix -> strong!ix ==* 0 ? (ignore, index1 $ targetIdx!ix)) indices --- Main --- ---- +-- Repa component -------------------------------------------------------------- --- TLM: should compare to a pre-saved reference image -run :: FilePath -> IO (() -> Acc (Array DIM2 Float)) -run file = do - pgm <- readPGM file - return (\() -> canny pgm) +-- | Trace out strong edges in the final image. +-- Also trace out weak edges that are connected to strong edges. +-- +wildfire + :: R.Array A R.DIM2 Float -- ^ Image with strong and weak edges set. + -> R.Array A R.DIM1 Int -- ^ Array containing flat indices of strong edges. + -> IO (R.Array U R.DIM2 Word8) + +wildfire img arrStrong + = do (sh, vec) <- wildfireIO + return $ sh `seq` vec `seq` R.fromUnboxed sh vec + + where lenImg = R.size $ R.extent img + lenStrong = R.size $ R.extent arrStrong + shImg = R.extent img + + wildfireIO + = do -- Stack of image indices we still need to consider. + vStrong <- R.toUnboxed `fmap` R.computeUnboxedP (R.delay arrStrong) + vStrong' <- V.thaw vStrong + vStack <- VM.grow vStrong' (lenImg - lenStrong) + + -- Burn in new edges. + vImg <- VM.unsafeNew lenImg + VM.set vImg 0 + burn vImg vStack lenStrong + vImg' <- V.unsafeFreeze vImg + return (R.extent img, vImg') + + + burn :: VM.IOVector Word8 -> VM.IOVector Int -> Int -> IO () + burn !vImg !vStack !top + | top == 0 + = return () + + | otherwise + = do let !top' = top - 1 + n <- VM.unsafeRead vStack top' + let (R.Z R.:. y R.:. x) = R.fromIndex (R.extent img) n + + let {-# INLINE push #-} + push ix t = + if R.inShape shImg ix + then pushWeak vImg vStack ix t + else return t + + VM.write vImg n 255 + >> push (R.Z R.:. y - 1 R.:. x - 1) top' + >>= push (R.Z R.:. y - 1 R.:. x ) + >>= push (R.Z R.:. y - 1 R.:. x + 1) + + >>= push (R.Z R.:. y R.:. x - 1) + >>= push (R.Z R.:. y R.:. x + 1) + + >>= push (R.Z R.:. y + 1 R.:. x - 1) + >>= push (R.Z R.:. y + 1 R.:. x ) + >>= push (R.Z R.:. y + 1 R.:. x + 1) + + >>= burn vImg vStack + + -- If this ix is weak in the source then set it to strong in the + -- result and push the ix onto the stack. + {-# INLINE pushWeak #-} + pushWeak vImg vStack ix top + = do let n = R.toIndex (R.extent img) ix + xDst <- VM.unsafeRead vImg n + let xSrc = img `R.unsafeIndex` ix + + if xDst == 0 + && xSrc == edge Weak + then do + VM.unsafeWrite vStack top (R.toIndex (R.extent img) ix) + return (top + 1) + + else return top +{-# NOINLINE wildfire #-} diff --git a/examples/tests/image-processing/IntegralImage.hs b/examples/tests/image-processing/IntegralImage.hs --- a/examples/tests/image-processing/IntegralImage.hs +++ b/examples/tests/image-processing/IntegralImage.hs @@ -2,48 +2,39 @@ module IntegralImage where -import PGM - -import Data.Array.Accelerate as Acc +import Data.Array.Accelerate as A +import Data.Array.Accelerate.IO as A -- |The value of each element in an integral image is the sum of all input elements -- above and to the left, inclusive. It is calculated by performing an inclusive/post -- scan from left-to-right then top-to-bottom. -- -integralImage :: (Elt a, IsNum a) => Array DIM2 a -> Acc (Array DIM2 a) +integralImage :: (Elt a, IsFloating a) => Array DIM2 Word32 -> Acc (Array DIM2 a) integralImage img = sumTable where -- scan rows - rowArr = reshape (lift $ Z:.(w * h)) arr - rowSegs = Acc.replicate (lift $ Z:.h) $ unit w - rowSum = reshape (lift (Z:.w:.h)) $ Acc.scanl1Seg (+) rowArr rowSegs + rowArr = reshape (index1 (w * h)) $ arr + rowSegs = A.replicate (index1 h) $ unit w + rowSum = reshape (index2 w h) $ scanl1Seg (+) rowArr rowSegs -- scan cols - colArr = reshape (lift $ Z:.(h * w)) $ transpose2D rowSum - colSegs = Acc.replicate (lift $ Z:.w) $ unit h - colSum = reshape (lift (Z:.h:.w)) $ Acc.scanl1Seg (+) colArr colSegs + colArr = reshape (index1 (h * w)) $ transpose rowSum + colSegs = A.replicate (index1 w) $ unit h + colSum = reshape (index2 h w) $ A.scanl1Seg (+) colArr colSegs -- transpose back - sumTable = transpose2D colSum + sumTable = transpose colSum -- - arr = use img + arr = A.map luminanceOfRGBA32 (use img) Z:.w:.h = unlift $ shape arr --- |Simple 2D matrix transpose. --- -transpose2D :: Elt a => Acc (Array DIM2 a) -> Acc (Array DIM2 a) -transpose2D arr = backpermute (swap $ shape arr) swap arr - where - swap = lift1 $ \(Z:.x:.y) -> Z:.y:.x :: Z :. Exp Int :. Exp Int - - -- Run integralImage over the input PGM -- run :: FilePath -> IO (() -> Acc (Array DIM2 Float)) run file = do - pgm <- readPGM file - return (\() -> integralImage pgm) + bmp <- either (error . show) id `fmap` readImageFromBMP file + return (\() -> integralImage bmp) diff --git a/examples/tests/image-processing/PGM.hs b/examples/tests/image-processing/PGM.hs deleted file mode 100644 --- a/examples/tests/image-processing/PGM.hs +++ /dev/null @@ -1,32 +0,0 @@ --- --- Load a PGM file. MacOS X users might find the following quicklook plugin --- useful for viewing PGM files: --- --- http://code.google.com/p/quicklook-pfm/ --- - -module PGM where - -import Control.Applicative -import Graphics.Pgm - -import Prelude as P -import Data.Array.Accelerate as Acc -import Data.Array.Unboxed hiding (Array) -import qualified Data.ByteString as B - - --- Read an 8-bit PGM file, and marshal to an Accelerate array as floating-point --- data in the range [0,1]. --- -readPGM :: FilePath -> IO (Array DIM2 Float) -readPGM fp = do - img <- either (error . show) head . pgmsToArrays <$> B.readFile fp :: IO (UArray (Int,Int) Word8) - return . fromIArray $ amap (\x -> P.fromIntegral x / 255) img - - -writePGM :: FilePath -> Array DIM2 Float -> IO () -writePGM fp img = - let arr = toIArray img :: UArray (Int,Int) Float - in arrayToFile fp $ amap (\x -> P.round (255 * x)) arr - diff --git a/examples/tests/io/BlockCopy.hs b/examples/tests/io/BlockCopy.hs --- a/examples/tests/io/BlockCopy.hs +++ b/examples/tests/io/BlockCopy.hs @@ -17,7 +17,7 @@ assertEqual preface expected actual = unless (actual == expected) (throw $ AssertionFailed msg) where - msg = (if null preface then "" else preface ++ "\n") ++ + msg = (if P.null preface then "" else preface ++ "\n") ++ "expected: " ++ show expected ++ "\n but got: " ++ show actual run :: IO () diff --git a/examples/tests/io/VectorCopy.hs b/examples/tests/io/VectorCopy.hs --- a/examples/tests/io/VectorCopy.hs +++ b/examples/tests/io/VectorCopy.hs @@ -1,74 +1,65 @@ -{-# LANGUAGE TemplateHaskell, TypeFamilies #-} -module VectorCopy where +module VectorCopy (run) where -import Data.Array.Accelerate hiding (fromList) -import Data.Array.Accelerate.Array.Sugar (EltRepr) -import Data.Array.Accelerate.IO +-- friends +import Data.Array.Accelerate +import Data.Array.Accelerate.IO ( toVectors, fromVectors ) -import Data.Vector.Storable +-- standard library +import Test.QuickCheck -import Foreign -import Test.QuickCheck -import Test.QuickCheck.All -import Test.QuickCheck.Monadic +-- Print expected/received message on inequality +-- +infix 4 .==. +(.==.) :: (Eq a, Show a) => a -> a -> Property +(.==.) ans ref = printTestCase message (ref == ans) + where + message = unlines ["*** Expected:", show ref + ,"*** Received:", show ans ] -roundtrip :: ( Arbitrary a - , Eq a - , Elt a - , Storable a - , BlockPtrs (EltRepr a) ~ ((), Ptr a) ) - => [a] -> Property -roundtrip xs = monadicIO $ do - let xsv = fromList xs - xsarr <- run $ fromVectorIO xsv - xsv' <- run $ toVectorIO xsarr - assert (xsv == xsv') -unsaferoundtrip :: ( Arbitrary a - , Eq a - , Elt a - , Storable a - , BlockPtrs (EltRepr a) ~ ((), Ptr a) ) - => [a] -> Bool -unsaferoundtrip xs = xsv == (toVector (fromVector xsv)) - where xsv = fromList xs +roundtrip :: (Arbitrary a, Eq a, Elt a) + => [a] + -> Property +roundtrip xs = + let n = length xs + sh = Z:.n + arr = fromList sh xs + in + xs .==. toList (fromVectors sh (toVectors arr)) + prop_Int8_roundtrip :: [Int8] -> Property prop_Int8_roundtrip = roundtrip -prop_Int8_unsaferoundtrip :: [Int8] -> Bool -prop_Int8_unsaferoundtrip = unsaferoundtrip prop_Int16_roundtrip :: [Int16] -> Property prop_Int16_roundtrip = roundtrip -prop_Int16_unsaferoundtrip :: [Int16] -> Bool -prop_Int16_unsaferoundtrip = unsaferoundtrip prop_Int32_roundtrip :: [Int32] -> Property prop_Int32_roundtrip = roundtrip -prop_Int32_unsaferoundtrip :: [Int32] -> Bool -prop_Int32_unsaferoundtrip = unsaferoundtrip prop_Int64_roundtrip :: [Int64] -> Property prop_Int64_roundtrip = roundtrip -prop_Int64_unsaferoundtrip :: [Int64] -> Bool -prop_Int64_unsaferoundtrip = unsaferoundtrip prop_Int_roundtrip :: [Int] -> Property prop_Int_roundtrip = roundtrip -prop_Int_unsaferoundtrip :: [Int] -> Bool -prop_Int_unsaferoundtrip = unsaferoundtrip prop_Float_roundtrip :: [Float] -> Property prop_Float_roundtrip = roundtrip -prop_Float_unsaferoundtrip :: [Float] -> Bool -prop_Float_unsaferoundtrip = unsaferoundtrip prop_Double_roundtrip :: [Double] -> Property prop_Double_roundtrip = roundtrip -prop_Double_unsaferoundtrip :: [Double] -> Bool -prop_Double_unsaferoundtrip = unsaferoundtrip -test :: IO Bool -test = $quickCheckAll + +run :: IO () +run = mapM_ quickCheck + [ property prop_Int8_roundtrip + , property prop_Int16_roundtrip + , property prop_Int32_roundtrip + , property prop_Int64_roundtrip + , property prop_Int_roundtrip + , property prop_Float_roundtrip + , property prop_Double_roundtrip + ] + diff --git a/examples/tests/primitives/Backpermute.hs b/examples/tests/primitives/Backpermute.hs --- a/examples/tests/primitives/Backpermute.hs +++ b/examples/tests/primitives/Backpermute.hs @@ -17,22 +17,14 @@ -- ----- reverseAcc :: Vector Float -> Acc (Vector Float) -reverseAcc xs = - let xs' = use xs - len = unindex1 (shape xs') - in - backpermute (shape xs') (\ix -> index1 $ len - (unindex1 ix) - 1) xs' +reverseAcc xs = Acc.reverse (use xs) reverseRef :: UArray Int Float -> UArray Int Float -reverseRef xs = listArray (bounds xs) (reverse (elems xs)) +reverseRef xs = listArray (bounds xs) (P.reverse (elems xs)) transposeAcc :: Acc.Array DIM2 Float -> Acc (Acc.Array DIM2 Float) -transposeAcc mat = - let mat' = use mat - swap = lift1 $ \(Z:.x:.y) -> Z:.y:.x :: Z:.Exp Int:.Exp Int - in - backpermute (swap $ shape mat') swap mat' +transposeAcc mat = Acc.transpose (use mat) transposeRef :: UArray (Int,Int) Float -> UArray (Int,Int) Float transposeRef mat = diff --git a/examples/tests/primitives/FoldSeg.hs b/examples/tests/primitives/FoldSeg.hs --- a/examples/tests/primitives/FoldSeg.hs +++ b/examples/tests/primitives/FoldSeg.hs @@ -49,7 +49,7 @@ -- generate elements -- - let x = P.fromIntegral $ sum (elems seg) + let x = P.fromIntegral $ P.sum (elems seg) vec <- randomUArrayR (-1,1) gen x vec' <- convertUArray vec diff --git a/examples/tests/primitives/Permute.hs b/examples/tests/primitives/Permute.hs --- a/examples/tests/primitives/Permute.hs +++ b/examples/tests/primitives/Permute.hs @@ -18,7 +18,7 @@ zeros = generate (constant (Z:. n-m)) (const 0) ones = generate (shape vec') (const 1) in - permute (+) zeros (\ix -> index1 $ Acc.floor (vec' Acc.! ix)) ones + permute (+) zeros (\ix -> index1 (Acc.floor (vec' Acc.! ix) :: Exp Int)) ones histogramRef :: (Int,Int) -> UArray Int Float -> UArray Int Int32 histogramRef (m,n) vec = diff --git a/examples/tests/primitives/ScanSeg.hs b/examples/tests/primitives/ScanSeg.hs --- a/examples/tests/primitives/ScanSeg.hs +++ b/examples/tests/primitives/ScanSeg.hs @@ -48,7 +48,7 @@ -- generate elements -- - let ne = P.fromIntegral $ sum (elems seg) + let ne = P.fromIntegral $ P.sum (elems seg) vec <- randomUArrayR (-1,1) gen ne vec' <- convertUArray vec -- diff --git a/examples/tests/simple/DotP.hs b/examples/tests/simple/DotP.hs --- a/examples/tests/simple/DotP.hs +++ b/examples/tests/simple/DotP.hs @@ -7,6 +7,7 @@ import System.Random.MWC import Data.Array.Unboxed import Data.Array.Accelerate as Acc +import Prelude as P -- Dot product @@ -23,7 +24,7 @@ -> UArray Int Float -> UArray () Float dotpRef xs ys - = listArray ((), ()) [sum [x * y | x <- elems xs | y <- elems ys]] + = listArray ((), ()) [P.sum [x * y | x <- elems xs | y <- elems ys]] -- Main diff --git a/examples/tests/simple/Filter.hs b/examples/tests/simple/Filter.hs --- a/examples/tests/simple/Filter.hs +++ b/examples/tests/simple/Filter.hs @@ -11,28 +11,11 @@ -- Filter -- ------ -filterAcc :: Elt a - => (Exp a -> Exp Bool) - -> Vector a - -> Acc (Vector a) -filterAcc p vec - = let arr = Acc.use vec - flags = Acc.map (boolToInt . p) arr - (targetIdx, len) = Acc.scanl' (+) 0 flags - arr' = Acc.backpermute (index1 $ the len) id arr - in - Acc.permute const arr' (\ix -> flags!ix ==* 0 ? (ignore, index1 $ targetIdx!ix)) arr - -- FIXME: This is abusing 'permute' in that the first two arguments are - -- only justified because we know the permutation function will - -- write to each location in the target exactly once. - -- Instead, we should have a primitive that directly encodes the - -- compaction pattern of the permutation function. +filterAcc :: Elt a => (Exp a -> Exp Bool) -> Vector a -> Acc (Vector a) +filterAcc p vec = Acc.filter p (use vec) -filterRef :: IArray UArray e - => (e -> Bool) - -> UArray Int e - -> UArray Int e +filterRef :: IArray UArray e => (e -> Bool) -> UArray Int e -> UArray Int e filterRef p xs = let xs' = Prelude.filter p (elems xs) in diff --git a/lib/ParseArgs.hs b/lib/ParseArgs.hs new file mode 100644 --- /dev/null +++ b/lib/ParseArgs.hs @@ -0,0 +1,127 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE PatternGuards #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE ViewPatterns #-} + +module ParseArgs ( + + Backend(..), parseArgs, run, run1, + module System.Console.GetOpt, + +) where + +import Data.Char +import Data.List +import Data.Label +import System.Exit +import System.Console.GetOpt +import qualified Criterion.Main as Criterion +import qualified Criterion.Config as Criterion + +import Data.Array.Accelerate ( Arrays, Acc ) +import qualified Data.Array.Accelerate.Interpreter as Interp +#ifdef ACCELERATE_CUDA_BACKEND +import qualified Data.Array.Accelerate.CUDA as CUDA +#endif + + +-- | Execute Accelerate expressions +-- +run :: Arrays a => Backend -> Acc a -> a +run Interpreter = Interp.run +#ifdef ACCELERATE_CUDA_BACKEND +run CUDA = CUDA.run +#endif + + +run1 :: (Arrays a, Arrays b) => Backend -> (Acc a -> Acc b) -> a -> b +run1 Interpreter f = Interp.run1 f +#ifdef ACCELERATE_CUDA_BACKEND +run1 CUDA f = CUDA.run1 f +#endif + + +-- | The set of backends available to execute the program +-- +data Backend = Interpreter +#ifdef ACCELERATE_CUDA_BACKEND + | CUDA +#endif + deriving (Eq, Bounded, Show) + +availableBackends :: (f :-> Backend) -> [OptDescr (f -> f)] +availableBackends backend = + [ Option [] ["interpreter"] + (NoArg (set backend Interpreter)) + "reference implementation (sequential)" + +#ifdef ACCELERATE_CUDA_BACKEND + , Option [] ["cuda"] + (NoArg (set backend CUDA)) + "implementation for NVIDIA GPUs (parallel)" +#endif + ] + + +-- | Complete the options set by appending a description of the available +-- execution backends. +-- +withBackends :: (f :-> Backend) -> [OptDescr (f -> f)] -> [OptDescr (f -> f)] +withBackends backend xs = availableBackends backend ++ xs + + +-- | Create the help text including a list of the available (and selected) +-- Accelerate backends. +-- +fancyHeader :: (config :-> Backend) -> config -> [String] -> [String] -> String +fancyHeader backend opts header footer = unlines (header ++ body ++ footer) + where + active this = if this == map toLower (show (get backend opts)) then "*" else "" + (ss,bs,ds) = unzip3 $ map (\(b,d) -> (active b, b, d)) $ concatMap extract (availableBackends backend) + table = zipWith3 paste (sameLen ss) (sameLen bs) ds + paste x y z = " " ++ x ++ " " ++ y ++ " " ++ z + sameLen xs = flushLeft ((maximum . map length) xs) xs + flushLeft n xs = [ take n (x ++ repeat ' ') | x <- xs ] + -- + extract (Option _ los _ descr) = + let losFmt = intercalate ", " los + in case lines descr of + [] -> [(losFmt, "")] + (x:xs) -> (losFmt, x) : [ ("",x') | x' <- xs ] + -- + body = "Available backends:" : table + + +-- | Process the command line arguments and return a tuple consisting of the +-- options structure, options for Criterion, and a list of unrecognised and +-- non-options. +-- +-- We drop any command line arguments following a "--". +-- +parseArgs :: (config :-> Bool) -- ^ access a help flag from the options structure + -> (config :-> Backend) -- ^ access the chosen backend from the options structure + -> [OptDescr (config -> config)] -- ^ the option descriptions + -> config -- ^ default option set + -> [String] -- ^ header text + -> [String] -- ^ footer text + -> [String] -- ^ command line arguments + -> IO (config, Criterion.Config, [String]) +parseArgs help backend (withBackends backend -> options) config header footer (takeWhile (/= "--") -> argv) = + let + helpMsg err = concat err + ++ usageInfo (unlines header) options + ++ usageInfo "\nGeneric criterion options:" Criterion.defaultOptions + + in case getOpt' Permute options argv of + (o,n,u,[]) -> do + + -- pass unrecognised options to criterion + (cconf, n') <- Criterion.parseArgs Criterion.defaultConfig Criterion.defaultOptions u + + case foldr id config o of + conf | False <- get help conf + -> putStrLn (fancyHeader backend conf header footer) >> return (conf, cconf, n ++ n') + _ -> putStrLn (helpMsg []) >> exitSuccess + + (_,_,_,err) -> error (helpMsg err) + diff --git a/lib/Random/Array.hs b/lib/Random/Array.hs new file mode 100644 --- /dev/null +++ b/lib/Random/Array.hs @@ -0,0 +1,90 @@ +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE TypeOperators #-} + +module Random.Array ( + + (:~>), + uniform, uniformR, + randomArray, randomArrayWithSeed, randomArrayWithSystemRandom, + +) where + +import Control.Monad.ST +import System.Random.MWC hiding ( uniform, uniformR ) +import qualified System.Random.MWC as R + +import Data.Array.Accelerate as A +import Data.Array.Accelerate.Array.Data as A +import Data.Array.Accelerate.Array.Sugar as Sugar + + +-- | A PNRG from indices to variates +-- +type sh :~> e = forall s. sh -> GenST s -> ST s e + + +-- | Uniformly distributed random variates. +-- +uniform :: (Shape sh, Elt e, Variate e) => sh :~> e +uniform _ = R.uniform + +-- | Uniformly distributed random variates in a given range. +-- +uniformR :: (Shape sh, Elt e, Variate e) => (e, e) -> sh :~> e +uniformR bounds _ = R.uniformR bounds + + +-- | Generate an array of random values using the supplied generator function. +-- The generator for variates is initialised with a fixed seed. +-- +randomArray :: (Shape sh, Elt e) => sh :~> e -> sh -> Array sh e +randomArray f sh + = let + (adata, _) = runArrayData $ do + gen <- create + arr <- runRandomArray f sh gen + return (arr, undefined) + in + adata `seq` Array (fromElt sh) adata + + +-- | Generate an array of random values using a supplied generator function and +-- seed value. +-- +randomArrayWithSeed :: (Shape sh, Elt e) => Seed -> sh :~> e -> sh -> Array sh e +randomArrayWithSeed seed f sh + = let + (adata, _) = runArrayData $ do + gen <- restore seed + arr <- runRandomArray f sh gen + return (arr, undefined) + in + adata `seq` Array (fromElt sh) adata + + +-- | Generate an array of random values using a supplied generator function, +-- initialised with the system's source of pseudo-random numbers. +-- +-- TODO: find a way to do this directly, without going via save/restore. +-- +randomArrayWithSystemRandom :: forall sh e. (Shape sh, Elt e) => sh :~> e -> sh -> IO (Array sh e) +randomArrayWithSystemRandom f sh + = do + seed <- withSystemRandom (asGenIO save) + return $! randomArrayWithSeed seed f sh + + +-- Common function to create a mutable array and fill it with random values +-- +runRandomArray :: (Shape sh, Elt e) => sh :~> e -> sh -> GenST s -> ST s (MutableArrayData s (EltRepr e)) +runRandomArray f sh gen + = let + n = Sugar.size sh + in do + arr <- newArrayData n + let write ix = unsafeWriteArrayData arr (Sugar.toIndex sh ix) + . fromElt =<< f ix gen + + iter sh write (>>) (return ()) + return arr + diff --git a/lib/Util.hs b/lib/Util.hs new file mode 100644 --- /dev/null +++ b/lib/Util.hs @@ -0,0 +1,19 @@ + +module Util where + +import Numeric +import Data.List + + +-- Utilities +-- +showFFloatSIBase :: RealFloat a => Maybe Int -> a -> a -> ShowS +showFFloatSIBase p b n + = showString + . nubBy (\x y -> x == ' ' && y == ' ') + $ showFFloat p n' [ ' ', si_unit ] + where + n' = n / (b ^^ (pow-4)) + pow = max 0 . min 8 . (+) 4 . floor $ logBase b n + si_unit = "pnµm kMGT" !! pow +