accelerate-examples 0.12.1.0 → 0.13.0.0
raw patch · 76 files changed
+4810/−1423 lines, 76 filesdep +accelerate-fftdep +cerealdep +repadep ~QuickCheckdep ~acceleratedep ~accelerate-cudanew-component:exe:accelerate-cannynew-component:exe:accelerate-hashcatnew-component:exe:accelerate-nbodynew-component:exe:accelerate-smoothlifenew-uploader
Dependencies added: accelerate-fft, cereal, repa, repa-io
Dependency ranges changed: QuickCheck, accelerate, accelerate-cuda, accelerate-io, array, attoparsec, bmp, bytestring, bytestring-lexing, cmdargs, criterion, cuda, deepseq, directory, fclabels, filepath, gloss, hashtables, mtl, mwc-random, pgm, pretty, test-framework, test-framework-quickcheck2, vector, vector-algorithms
Files
- accelerate-examples.cabal +219/−86
- examples/canny/src-acc/Canny.hs +209/−0
- examples/canny/src-acc/Config.hs +72/−0
- examples/canny/src-acc/Main.hs +80/−0
- examples/canny/src-acc/Wildfire.hs +102/−0
- examples/crystal/Config.hs +17/−50
- examples/crystal/Main.hs +10/−15
- examples/fluid/src-acc/Config.hs +309/−0
- examples/fluid/src-acc/Event.hs +83/−0
- examples/fluid/src-acc/Fluid.hs +172/−0
- examples/fluid/src-acc/Main.hs +91/−0
- examples/fluid/src-acc/Type.hs +67/−0
- examples/fluid/src-acc/World.hs +122/−0
- examples/fluid/src/Config.hs +0/−168
- examples/fluid/src/Event.hs +0/−81
- examples/fluid/src/Fluid.hs +0/−169
- examples/fluid/src/Main.hs +0/−84
- examples/fluid/src/Type.hs +0/−67
- examples/fluid/src/World.hs +0/−198
- examples/hashcat/Config.hs +79/−0
- examples/hashcat/Digest.hs +97/−0
- examples/hashcat/MD5.hs +164/−0
- examples/hashcat/Main.hs +80/−0
- examples/mandelbrot/Config.hs +32/−60
- examples/mandelbrot/Main.hs +51/−135
- examples/mandelbrot/Mandel.hs +195/−0
- examples/mandelbrot/World.hs +174/−0
- examples/n-body/Common/Body.hs +180/−0
- examples/n-body/Common/Dump.hs +52/−0
- examples/n-body/Common/Tree.hs +7/−0
- examples/n-body/Common/Type.hs +39/−0
- examples/n-body/Common/Util.hs +76/−0
- examples/n-body/Common/World.hs +63/−0
- examples/n-body/Config.hs +163/−0
- examples/n-body/Gloss/Draw.hs +113/−0
- examples/n-body/Gloss/Event.hs +27/−0
- examples/n-body/Gloss/Simulate.hs +48/−0
- examples/n-body/Main.hs +87/−0
- examples/n-body/Random/Position.hs +44/−0
- examples/n-body/Solver/BarnsHut.hs +12/−0
- examples/n-body/Solver/Naive.hs +23/−0
- examples/quickcheck/Arbitrary/Array.hs +2/−2
- examples/quickcheck/Arbitrary/Shape.hs +2/−2
- examples/quickcheck/Config.hs +9/−8
- examples/quickcheck/Main.hs +1/−1
- examples/quickcheck/Test/Base.hs +10/−13
- examples/quickcheck/Test/IndexSpace.hs +4/−13
- examples/quickcheck/Test/Mapping.hs +3/−3
- examples/quickcheck/Test/PrefixSum.hs +107/−20
- examples/quickcheck/Test/Reduction.hs +5/−11
- examples/smoothlife/Config.hs +232/−0
- examples/smoothlife/Gloss/Draw.hs +72/−0
- examples/smoothlife/Gloss/Event.hs +11/−0
- examples/smoothlife/Gloss/Simulate.hs +35/−0
- examples/smoothlife/Main.hs +69/−0
- examples/smoothlife/Random/Splat.hs +56/−0
- examples/smoothlife/SmoothLife.hs +160/−0
- examples/tests/Benchmark.hs +7/−1
- examples/tests/Config.hs +7/−1
- examples/tests/Main.hs +1/−0
- examples/tests/PGM.hs +32/−0
- examples/tests/Test.hs +15/−19
- examples/tests/image-processing/Canny.hs +314/−81
- examples/tests/image-processing/IntegralImage.hs +13/−22
- examples/tests/image-processing/PGM.hs +0/−32
- examples/tests/io/BlockCopy.hs +1/−1
- examples/tests/io/VectorCopy.hs +36/−45
- examples/tests/primitives/Backpermute.hs +3/−11
- examples/tests/primitives/FoldSeg.hs +1/−1
- examples/tests/primitives/Permute.hs +1/−1
- examples/tests/primitives/ScanSeg.hs +1/−1
- examples/tests/simple/DotP.hs +2/−1
- examples/tests/simple/Filter.hs +3/−20
- lib/ParseArgs.hs +127/−0
- lib/Random/Array.hs +90/−0
- lib/Util.hs +19/−0
accelerate-examples.cabal view
@@ -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:+ <http://hackage.haskell.org/package/accelerate>+ .+ License: BSD3 License-file: LICENSE Author: The Accelerate Team Maintainer: Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+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
+ examples/canny/src-acc/Canny.hs view
@@ -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+
+ examples/canny/src-acc/Config.hs view
@@ -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 = []+
+ examples/canny/src-acc/Main.hs view
@@ -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+ ]+ ]+
+ examples/canny/src-acc/Wildfire.hs view
@@ -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+
examples/crystal/Config.hs view
@@ -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 = []
examples/crystal/Main.hs view
@@ -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
+ examples/fluid/src-acc/Config.hs view
@@ -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)+
+ examples/fluid/src-acc/Event.hs view
@@ -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)+
+ examples/fluid/src-acc/Fluid.hs view
@@ -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)+
+ examples/fluid/src-acc/Main.hs view
@@ -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+
+ examples/fluid/src-acc/Type.hs view
@@ -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)+
+ examples/fluid/src-acc/World.hs view
@@ -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)+
− examples/fluid/src/Config.hs
@@ -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"- ]-
− examples/fluid/src/Event.hs
@@ -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)-
− examples/fluid/src/Fluid.hs
@@ -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)-
− examples/fluid/src/Main.hs
@@ -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-
− examples/fluid/src/Type.hs
@@ -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)-
− examples/fluid/src/World.hs
@@ -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)-
+ examples/hashcat/Config.hs view
@@ -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 = []+
+ examples/hashcat/Digest.hs view
@@ -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))+
+ examples/hashcat/MD5.hs view
@@ -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+
+ examples/hashcat/Main.hs view
@@ -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 (:"+
examples/mandelbrot/Config.hs view
@@ -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)"+ ]
examples/mandelbrot/Main.hs view
@@ -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
+ examples/mandelbrot/Mandel.hs view
@@ -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+--}
+ examples/mandelbrot/World.hs view
@@ -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')+
+ examples/n-body/Common/Body.hs view
@@ -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+
+ examples/n-body/Common/Dump.hs view
@@ -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+
+ examples/n-body/Common/Tree.hs view
@@ -0,0 +1,7 @@+--+-- A Barns-Hut tree+--++module Common.Tree+ where+
+ examples/n-body/Common/Type.hs view
@@ -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)+
+ examples/n-body/Common/Util.hs view
@@ -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)+
+ examples/n-body/Common/World.hs view
@@ -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' }+
+ examples/n-body/Config.hs view
@@ -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 = []+
+ examples/n-body/Gloss/Draw.hs view
@@ -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+
+ examples/n-body/Gloss/Event.hs view
@@ -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+
+ examples/n-body/Gloss/Simulate.hs view
@@ -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+
+ examples/n-body/Main.hs view
@@ -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+
+ examples/n-body/Random/Position.hs view
@@ -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+
+ examples/n-body/Solver/BarnsHut.hs view
@@ -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!"+
+ examples/n-body/Solver/Naive.hs view
@@ -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+
examples/quickcheck/Arbitrary/Array.hs view
@@ -1,6 +1,6 @@+{-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE FlexibleInstances #-}-{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE TypeOperators #-} module Arbitrary.Array where
examples/quickcheck/Arbitrary/Shape.hs view
@@ -1,6 +1,6 @@+{-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE FlexibleInstances #-}-{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE TypeOperators #-} module Arbitrary.Shape where
examples/quickcheck/Config.hs view
@@ -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
examples/quickcheck/Main.hs view
@@ -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.
examples/quickcheck/Test/Base.hs view
@@ -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
examples/quickcheck/Test/IndexSpace.hs view
@@ -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))
examples/quickcheck/Test/Mapping.hs view
@@ -1,7 +1,7 @@-{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE RankNTypes #-} {-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE TypeOperators #-}-{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE TypeOperators #-} module Test.Mapping where
examples/quickcheck/Test/PrefixSum.hs view
@@ -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 )
examples/quickcheck/Test/Reduction.hs view
@@ -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)
+ examples/smoothlife/Config.hs view
@@ -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])+ ]+
+ examples/smoothlife/Gloss/Draw.hs view
@@ -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+
+ examples/smoothlife/Gloss/Event.hs view
@@ -0,0 +1,11 @@++module Gloss.Event+ where++import Gloss.Simulate+import Graphics.Gloss.Interface.Pure.Game+++react :: Event -> Simulate -> Simulate+react _ = id+
+ examples/smoothlife/Gloss/Simulate.hs view
@@ -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+ }+
+ examples/smoothlife/Main.hs view
@@ -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+
+ examples/smoothlife/Random/Splat.hs view
@@ -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+
+ examples/smoothlife/SmoothLife.hs view
@@ -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' ))
examples/tests/Benchmark.hs view
@@ -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
examples/tests/Config.hs view
@@ -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"
examples/tests/Main.hs view
@@ -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
+ examples/tests/PGM.hs view
@@ -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+
examples/tests/Test.hs view
@@ -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)
examples/tests/image-processing/Canny.hs view
@@ -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 #-}
examples/tests/image-processing/IntegralImage.hs view
@@ -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)
− examples/tests/image-processing/PGM.hs
@@ -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-
examples/tests/io/BlockCopy.hs view
@@ -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 ()
examples/tests/io/VectorCopy.hs view
@@ -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+ ]+
examples/tests/primitives/Backpermute.hs view
@@ -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 =
examples/tests/primitives/FoldSeg.hs view
@@ -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
examples/tests/primitives/Permute.hs view
@@ -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 =
examples/tests/primitives/ScanSeg.hs view
@@ -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 --
examples/tests/simple/DotP.hs view
@@ -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
examples/tests/simple/Filter.hs view
@@ -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
+ lib/ParseArgs.hs view
@@ -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)+
+ lib/Random/Array.hs view
@@ -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+
+ lib/Util.hs view
@@ -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+