repa-examples 1.0.0.0 → 1.1.0.0
raw patch · 16 files changed
+302/−993 lines, 16 filesdep +repa-algorithmsdep +repa-iodep ~basedep ~repanew-component:exe:repa-fft2dnew-component:exe:repa-fft2d-highpass
Dependencies added: repa-algorithms, repa-io
Dependency ranges changed: base, repa
Files
- FFT/src/DFT.hs +0/−61
- FFT/src/FFT.hs +0/−136
- FFT/src/FFT2d/Main.hs +52/−0
- FFT/src/HighPass/Main.hs +73/−0
- FFT/src/Main.hs +0/−194
- FFT/src/Roots.hs +0/−43
- FFT/src/StrictComplex.hs +0/−32
- LICENSE +22/−11
- Laplace/src/Main.hs +83/−26
- MMult/src/Main.hs +29/−7
- MMult/src/Solver.hs +0/−14
- lib/ColorRamp.hs +0/−47
- lib/Matrix.hs +0/−86
- lib/PPM.hs +0/−171
- lib/Vector.hs +0/−137
- repa-examples.cabal +43/−28
− FFT/src/DFT.hs
@@ -1,61 +0,0 @@-{-# LANGUAGE TypeOperators #-}--module DFT - ( dft- , idft- , dftWithRoots)-where-import Data.Array.Repa as A-import Data.Ratio-import StrictComplex-import Roots---- | Compute the (non-fast) Discrete Fourier Transform of a vector.-dft :: Shape sh- => Array (sh :. Int) Complex- -> Array (sh :. Int) Complex--dft v- = let rofu = calcRofu (extent v)- in force $ dftWithRoots rofu v----- | Compute the inverse (non-fast) Discrete Fourier Transform of a vector.-idft :: Shape sh- => Array (sh :. Int) Complex- -> Array (sh :. Int) Complex--idft v- = let _ :. len = extent v- scale = fromIntegral len :*: 0- rofu = calcInverseRofu (extent v)- in force $ A.map (/ scale) $ dftWithRoots rofu v----- | Generic function for computation of forward or inverse Discrete Fourier Transforms.-dftWithRoots- :: Shape sh- => Array (sh :. Int) Complex -- ^ Roots of unity for this vector length.- -> Array (sh :. Int) Complex -- ^ Input vector.- -> Array (sh :. Int) Complex--dftWithRoots rofu arr- = traverse arr id (\_ k -> dftK rofu arr k)- ---- | Compute one value of the DFT.-dftK :: Shape sh- => Array (sh :. Int) Complex -- ^ Roots of unity for this vector length.- -> Array (sh :. Int) Complex -- ^ Input vector.- -> (sh :. Int) -- ^ Index of the value we want.- -> Complex--dftK rofu arrX (_ :. k)- = A.sumAll $ A.zipWith (*) arrX wroots- where sh@(_ :. len) = extent arrX-- -- All the roots we need to multiply with.- wroots = fromFunction sh elemFn- elemFn (sh :. n) - = rofu !: (sh :. (k * n) `mod` len)-
− FFT/src/FFT.hs
@@ -1,136 +0,0 @@-{-# LANGUAGE TypeOperators, PatternGuards #-}---- | Computation of Fast Fourier Transforms using the Cooley-Tuckey algorithm.-module FFT - ( fft- , ifft- , fft2d- , fftWithRoots )-where-import Data.Array.Repa as A-import Data.Ratio-import StrictComplex-import Roots---- Vector Transform ---------------------------------------------------------------------------------- | Compute the (fast) Discrete Fourier Transform of a vector.-fft :: Shape sh- => Array (sh :. Int) Complex- -> Array (sh :. Int) Complex--fft v- = let rofu = calcRofu (extent v)- in force $ fftWithRoots rofu v----- | Compute the (fast) Inverse Discrete Fourier Transform of a vector.-ifft :: Shape sh- => Array (sh :. Int) Complex- -> Array (sh :. Int) Complex--ifft v- = let _ :. len = extent v- scale = fromIntegral len :*: 0- rofu = calcInverseRofu (extent v)- in force $ A.map (/ scale) $ fftWithRoots rofu v----- Matrix Transform ---------------------------------------------------------------------------------- | Compute the (fast) Discrete Fourier Transform of a square matrix.-fft2d :: Array DIM2 Complex- -> Array DIM2 Complex--fft2d arr- = let rofu = calcRofu (extent arr)- fftTrans = transpose . fftWithRoots rofu- in fftTrans $ fftTrans arr----- Cube Transform ------------------------------------------------------------------------------------ | Compute the (fast) Discrete Fourier Transform of a 3d array.-fft3d :: Array DIM3 Complex- -> Array DIM3 Complex--fft3d arr- = let rofu = calcRofu (extent arr)-- transpose3 arr- = traverse arr - (\(Z :. k :. l :. m) -> (Z :. l :. m :. k)) - (\f (Z :. l :. m :. k) -> f (Z :. k :. l :. m)) -- fftTrans = transpose3 . fftWithRoots rofu- - in fftTrans $ fftTrans $ fftTrans arr-- --- Worker -------------------------------------------------------------------------------------------- | Generic function for computation of forward or inverse Discrete Fourier Transforms.--- The length of the roots vector must be the same as the values vector.--- The length of these vectors must be a power of two.-fftWithRoots - :: Shape sh- => Array (sh :. Int) Complex -- ^ Roots of unity.- -> Array (sh :. Int) Complex -- ^ Input values.- -> Array (sh :. Int) Complex--fftWithRoots rofu v- | not $ (denominator $ toRational (logBase 2 $ fromIntegral vLen)) == 1- = error $ "fft: vector length of " ++ show vLen ++ " is not a power of 2"- - | rLen /= vLen- = error $ "fft: length of vector is not the length of the roots"- - | otherwise- = fftWithRoots' rofu v-- where _ :. rLen = extent rofu- _ :. vLen = extent v---fftWithRoots'- :: Shape sh- => Array (sh :. Int) Complex- -> Array (sh :. Int) Complex- -> Array (sh :. Int) Complex--{-# INLINE fftWithRoots' #-}-fftWithRoots' rofu v- = case extent v of- _ :. 2 -> fft_two rofu v- dim -> fft_split rofu v dim--{-# INLINE fft_two #-}-fft_two rofu v- = let vFn' vFn (sh :. 0) = vFn (sh :. 0) + vFn (sh :. 1)- vFn' vFn (sh :. 1) = vFn (sh :. 0) - vFn (sh :. 1)- in traverse v id vFn'- -{-# INLINE fft_split #-}-fft_split rofu v vLen- = let fft_lr = force $ fftWithRoots' (splitRofu rofu) (splitVector v)-- fft_l = traverse2 fft_lr rofu - (\(sh :. 2 :. n) _ -> sh :. n)- (\f r (sh :. i) -> f (sh :. 0 :. i) + r (sh :. i) * f (sh :. 1 :. i))-- fft_r = traverse2 fft_lr rofu - (\(sh :. 2 :. n) _ -> sh :. n)- (\f r (sh :. i) -> f (sh :. 0 :. i) - r (sh :. i) * f (sh :. 1 :. i))-- in fft_l +:+ fft_r--{-# INLINE splitRofu #-}-splitRofu rofu- = traverse rofu- (\(rSh :. rLen) -> rSh :. (2::Int) :. (rLen `div` 2))- (\rFn (sh :. _ :. i) -> rFn (sh :. 2*i))--{-# INLINE splitVector #-}-splitVector v - = let vFn' vFn (sh :. 0 :. i) = vFn (sh :. 2*i)- vFn' vFn (sh :. 1 :. i) = vFn (sh :. 2*i+1)- in traverse v- (\(vSh :. vLen) -> vSh :. 2 :. (vLen `div` 2)) - vFn'-
+ FFT/src/FFT2d/Main.hs view
@@ -0,0 +1,52 @@++-- | Perform the 2D FFT on a BMP image.+import Data.Array.Repa.Algorithms.FFT+import Data.Array.Repa.Algorithms.DFT.Center+import Data.Array.Repa.Algorithms.Complex+import Data.Array.Repa.IO.BMP+import Data.Array.Repa as A+import System.Environment+import Control.Monad++main :: IO ()+main + = do args <- getArgs+ case args of+ [fileIn, clipMag, fileMag, filePhase]+ -> mainWithArgs fileIn (read clipMag) fileMag filePhase++ _ -> putStr $ unlines+ [ "Usage: repa-fft2d <fileIn.bmp> <clip mag :: Int> <fileOutMag.bmp> <fileOutPhase.bmp>"+ , ""+ , " The output magnitude has a high dynamic range. We need to clip it otherwise"+ , " most of the pixels in the output BMP will be black. Start with a value equal"+ , " to about the width of the image (eg 512)"+ , "" ]+ + +mainWithArgs fileIn clipMag fileMag filePhase+ = do + -- Load in the matrix.+ arrReal <- liftM (either (\e -> error $ show e) id)+ $ readMatrixFromGreyscaleBMP fileIn+ let arrComplex = force $ A.map (\r -> r :*: 0) arrReal+ + -- Apply the centering transform so that the output has the zero+ -- frequency in the middle of the image.+ let arrCentered = centerMatrix arrComplex+ + -- Do the 2d transform.+ let arrFreq = fft2d arrCentered+ + -- Write out the magnitude of the transformed array, + -- clipping it at the given value.+ let clip m = if m >= clipMag then clipMag else m+ let arrMag = A.map (clip . mag) arrFreq+ writeMatrixToGreyscaleBMP fileMag arrMag++ -- Write out the phase of the transofmed array, + -- scaling it to make full use of the 8 bit greyscale.+ let scaledArg x = (arg x + pi) * (255 / (2 * pi))+ let arrPhase = A.map scaledArg arrFreq+ writeMatrixToGreyscaleBMP filePhase arrPhase+
+ FFT/src/HighPass/Main.hs view
@@ -0,0 +1,73 @@+{-# LANGUAGE PatternGuards #-}++-- | Perform high pass filtering on a BMP image.+import Data.Array.Repa.Algorithms.FFT+import Data.Array.Repa.Algorithms.DFT.Center+import Data.Array.Repa.Algorithms.Complex+import Data.Array.Repa.IO.BMP+import Data.Array.Repa as A+import System.Environment+import Control.Monad++main :: IO ()+main + = do args <- getArgs+ case args of+ [cutoff, fileIn, fileOut]+ -> mainWithArgs (read cutoff) fileIn fileOut++ _ -> putStr $ unlines+ [ "Usage: repa-fft-highpass <cutoff::Int> <fileIn.bmp> <fileOut.bmp>"+ , "" ]+ + +mainWithArgs cutoff fileIn fileOut+ = do + -- Load in the matrix.+ (arrRed, arrGreen, arrBlue)+ <- liftM (either (\e -> error $ show e) id)+ $ readComponentsFromBMP fileIn+ + -- Do the transform on each component individually+ let arrRed' = transform cutoff arrRed+ let arrGreen' = transform cutoff arrGreen+ let arrBlue' = transform cutoff arrBlue+ + -- Write it back to file.+ writeComponentsToBMP fileOut+ arrRed' arrGreen' arrBlue' + ++transform cutoff arrReal+ = let arrComplex = force $ A.map (\r -> (fromIntegral r) :*: 0) arrReal+ + -- Do the 2d transform.+ arrCentered = centerMatrix arrComplex+ arrFreq = fft2d arrCentered++ -- Zap out the low frequency components.+ _ :. height :. width = extent arrFreq+ centerX = width `div` 2+ centerY = height `div` 2+ + {-# INLINE highpass #-}+ highpass get ix@(_ :. y :. x)+ | x > centerX + cutoff+ || x < centerX - cutoff+ || y > centerY + cutoff+ || y < centerY - cutoff+ = get ix+ + | otherwise+ = 0+ + arrFilt = traverse arrFreq id highpass++ -- Do the inverse transform to get back to image space.+ arrInv = ifft2d arrFilt+ + -- Write out the magnitude of the transformed array, + arrMag = A.map (truncate . mag) arrInv++ in arrMag+
− FFT/src/Main.hs
@@ -1,194 +0,0 @@-{-# LANGUAGE ParallelListComp, PatternGuards, ScopedTypeVariables #-}--module Main where--import FFT-import DFT-import Roots-import StrictComplex-import Vector-import Matrix-import PPM-import ColorRamp--import Data.Array.Repa as A-import Data.List as L-import Data.Maybe -import Prelude as P-import Control.Monad-import System.Environment----- Arg Parsing -------------------------------------------------------------------------------------data Arg- -- | Use DFT instead of FFT- = ArgDFT -- -- | Compute the inverse transform.- | ArgInverse- - -- | Use a vector step function for input.- | ArgVectorRealStep - Int -- length of 'on' part- Int -- total length of vector-- -- | Read the a real valued vector from this file.- | ArgVectorReal- FilePath-- -- | Read a PPM file as input.- | ArgPPM - FilePath-- -- | Write the result to this file.- | ArgOutMagnitude- FilePath-- -- | Write the magnitude of the transformed matrix as a PPM file.- | ArgOutPPMMagnitude- FilePath- - -- | Clip tranformed values to this level.- | ArgOutPPMClip- Double- - deriving (Show, Eq)---parseArgs [] = []-parseArgs (flag:xx)- | "-dft" <- flag- = ArgDFT : parseArgs xx-- | "-inverse" <- flag- = ArgInverse : parseArgs xx-- | "-vector-real-step" <- flag- , onlen:len:rest <- xx- = ArgVectorRealStep (read onlen) (read len) : parseArgs rest-- | "-vector-real" <- flag- , fileName:rest <- xx- = ArgVectorReal fileName : parseArgs rest- - | "-ppm" <- flag- , fileName:rest <- xx- = ArgPPM fileName : parseArgs rest-- | "-ppm-clip" <- flag- , num:rest <- xx- = ArgOutPPMClip (read num) : parseArgs rest- - | "-out-magnitude" <- flag- , fileName:rest <- xx- = ArgOutMagnitude fileName : parseArgs rest- - | "-out-ppm-magnitude" <- flag- , fileName:rest <- xx- = ArgOutPPMMagnitude fileName : parseArgs rest- - | "-out-ppm-clip" <- flag- , level:rest <- xx- = ArgOutPPMClip (read level) : parseArgs rest- - | otherwise - = error $ "bad arg " ++ flag ++ "\n"---help = unlines- [ "Usage: fft [args..]"- , ""- , " -dft Use (slow) DFT instead of (fast) FFT."- , " -inverse Compute the inverse transform."- , ""- , "INPUT:"- , " -vector-real <filename> Read a real valued input vector from this file."- , " -vector-real-step <onlen::Int> <len::Int> Use a real valued step function for input."- , " -ppm <filename> Use a PPM file for input."- , ""- , "OUTPUT:"- , " -out-vector-magnitude <file-name> Write the magnitute of the transformed vector to file."- , " -out-ppm-magnitude <file-name> Write transformed matrix to a ppm file."- , " -out-ppm-clip <val::Int> ... while clipping transformed values to this level." - , ""- , "NOTE:" - , " - For the fast algorithm, the length of the input vector/dimensions of the array"- , " must be powers of two."- , ""- , " - When using PPM input files, pixels are converted to grey-scale and then used as"- , " the real values of the input matrix." ]- --- Main --------------------------------------------------------------------------------------------main :: IO ()-main - = do args <- liftM parseArgs $ getArgs- - -- Decide which algorithm to use- let alg | elem ArgDFT args - , elem ArgInverse args = idft-- | elem ArgDFT args = dft- - | elem ArgInverse args = ifft- - | otherwise = fft-- main' args alg---main' args alg-- -- | Transform a real-valued step function - | [(onLength, vecLength)] <- [(ol, vl) | ArgVectorRealStep ol vl <- args]- = let offLength = vecLength - onLength- step_real = P.replicate onLength 1 ++ P.replicate offLength 0- step = P.map (:*: 0) step_real- - arr = A.fromList (Z :. vecLength) step- arrT = alg arr- in outVector args arrT- - -- | Transform some vector from a file- | [fileName] <- [f | ArgVectorReal f <- args]- = do arr_real :: Array DIM1 Double <- readVectorFromTextFile fileName - let arr = A.map (\r -> r :*: 0) arr_real- let arrT = alg arr- outVector args arrT- - -- | Transform a PPM file.- | [fileName] <- [f | ArgPPM f <- args]- = do let loadPixel r g b = sqrt (fromIntegral r^2 + fromIntegral g^2 + fromIntegral b^2)- arr_double <- readPPMAsMatrix loadPixel fileName- let arr_real = (A.map (\r -> r :*: 0) arr_double) :: Array DIM2 Complex- let arrT = fft2d arr_real-- outPPM args arrT- - -- Not sure what you mean...- | otherwise- = putStr help- - -outVector args vec- | [fileName] <- [f | ArgOutMagnitude f <- args ]- = writeVectorAsTextFile (A.map mag vec) fileName- - | otherwise- = return ()--outPPM :: [Arg] -> Array DIM2 Complex -> IO ()-outPPM args arr- | [fileName] <- [f | ArgOutPPMMagnitude f <- args ]- , mClipLevel <- listToMaybe [l | ArgOutPPMClip l <- args]- = do let arr_mag = A.map mag arr- let arr_clipped = maybe arr_mag- (\level -> A.map (\x -> if x > level then level else x) arr_mag)- mClipLevel- - writeMatrixAsNormalisedPPM fileName - pixelGrey arr_clipped- - -pixelColor x = (x, x, x)-pixelGrey x = (x, x, x)-
− FFT/src/Roots.hs
@@ -1,43 +0,0 @@-{-# LANGUAGE TypeOperators #-}--module Roots- ( calcRofu- , calcInverseRofu)-where-import Data.Array.Repa-import StrictComplex----- Roots of Unity ------------------------------------------------------------------------------------- | Fill a vector with roots of unity (Rofu)-calcRofu - :: Shape sh- => (sh :. Int) -- ^ Length of resulting vector.- -> Array (sh :. Int) Complex--calcRofu sh@(_ :. n) - = force $ fromFunction sh f- where- f :: Shape sh => (sh :. Int) -> Complex- f (_ :. i) = (cos (2 * pi * (fromIntegral i) / len))- :*: (- sin (2 * pi * (fromIntegral i) / len))-- len = fromIntegral n----- | Fill a vector with roots of unity (Rofu)--- for the inverse transform.-calcInverseRofu- :: Shape sh- => (sh :. Int) -- ^ Length of resulting vector.- -> Array (sh :. Int) Complex--calcInverseRofu sh@(_ :. n) - = force $ fromFunction sh f- where- f :: Shape sh => (sh :. Int) -> Complex- f (_ :. i) = (cos (2 * pi * (fromIntegral i) / len))- :*: (sin (2 * pi * (fromIntegral i) / len))-- len = fromIntegral n
− FFT/src/StrictComplex.hs
@@ -1,32 +0,0 @@-{-# LANGUAGE TypeOperators, TypeSynonymInstances #-}--module StrictComplex- ( Complex(..)- , mag- , (:*:)(..))-where-import Data.Array.Parallel.Base ((:*:)(..))---- | Strict complex doubles.-type Complex - = Double :*: Double--instance Num Complex where- (r :*: i) + (r' :*: i') = r+r' :*: i+i'- (r :*: i) - (r' :*: i') = r-r' :*: i-i'- (r :*: i) * (r' :*: i') = r*r' - i*i' :*: r*i' + r'*i- fromInteger n = fromInteger n :*: 0.0--instance Fractional Complex where- (a :*: b) / (c :*: d) - = let den = c^2 + d^2- re = (a * c + b * d) / den- im = (b * c - a * d) / den- in re :*: im- - --- | Take the magnitude of a complex number.-mag :: Complex -> Double-mag (r :*: i) = sqrt (r * r + i * i)--
LICENSE view
@@ -1,13 +1,24 @@-Copyright (c) 2010 The DPH Team+Copyright (c) 2010, University of New South Wales.+All rights reserved. - Permission is hereby granted, free of charge, to any person- obtaining a copy of this software and associated documentation- files (the "Software"), to deal in the Software without- restriction, including without limitation the rights to use,- copy, modify, merge, publish, distribute, sublicense, and/or sell- copies of the Software, and to permit persons to whom the- Software is furnished to do so, subject to the following- condition:+Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:+ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.+ * Redistributions in binary form must reproduce the above copyright+ notice, this list of conditions and the following disclaimer in the+ documentation and/or other materials provided with the distribution.+ * Neither the name of the University of New South Wales nor the+ names of its contributors may be used to endorse or promote products+ derived from this software without specific prior written permission. - The above copyright notice and this permission notice shall be- included in all copies or substantial portions of the Software.+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ''AS IS'' AND ANY+EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE+DISCLAIMED. IN NO EVENT SHALL COPYRIGHT HOLDERS BE LIABLE FOR ANY+DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES+(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND+ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
Laplace/src/Main.hs view
@@ -1,16 +1,16 @@ {-# LANGUAGE BangPatterns #-} -- | Solver for the Laplace equation--- Writes a PPM file of the completed solution.------ You can use the ImageMagick convert program to make a png--- with "convert out.ppm out.png"+-- You supply a BMP files specifying the boundary conditions.+-- The output is written back to another BMP file. -- import Solver import Data.Array.Repa as A-import ColorRamp-import PPM+import Data.Array.Repa.IO.BMP +import Data.Array.Repa.IO.ColorRamp import System.Environment+import Data.Word+import Control.Monad main :: IO () main @@ -27,11 +27,11 @@ -- | Command line usage information. usage :: String usage = unlines- [ "Usage: laplace <iterations> <input.ppm> <output.ppm>"+ [ "Usage: laplace <iterations> <input.bmp> <output.bmp>" , "" , " iterations :: Int Number of iterations to use in the solver."- , " input.ppm :: FileName ASCII 8 bit RGB PPM file for initial and boundary values."- , " output.ppm :: FileName PPM file to write output to."+ , " input.bmp :: FileName Uncompressed RGB24 or RGBA32 BMP file for initial and boundary values."+ , " output.bmp :: FileName BMP file to write output to." , "" , " Format of input file:" , " Boundary values are indicated in greyscale,"@@ -52,38 +52,95 @@ laplace steps fileInput fileOutput = do -- Load up the file containing boundary conditions.- (matBoundMask, matBoundValue) - <- readPPMAsMatrix2 loadPixel fileInput+ arrImage <- liftM (either (error . show) id)+ $ readImageFromBMP fileInput + let arrBoundValue = force $ slurpDoublesFromImage slurpBoundValue arrImage+ let arrBoundMask = force $ slurpDoublesFromImage slurpBoundMask arrImage+ -- Use the boundary condition values as the initial matrix.- let matInitial = matBoundValue+ let arrInitial = arrBoundValue -- Run the solver.- let matFinal = solveLaplace+ let arrFinal = solveLaplace steps- matBoundMask- matBoundValue- matInitial+ arrBoundMask+ arrBoundValue+ arrInitial - matFinal `deepSeqArray` return ()+ arrFinal `deepSeqArray` return () - -- Write out the matrix as a colorised PPM image - writeMatrixAsNormalisedPPM+ -- Make the result image+ let arrImageOut = makeImageFromDoubles (rampColorHotToCold 0.0 1.0) arrFinal++ -- Write out the image to file. + writeImageToBMP fileOutput- (rampColorHotToCold 0.0 1.0)- matFinal+ arrImageOut --- | Extract boundary mask and value from a pixel value.-loadPixel :: Int -> Int -> Int -> (Double, Double)-loadPixel r g b++slurpDoublesFromImage+ :: (Word8 -> Word8 -> Word8 -> Double)+ -> Array DIM3 Word8+ -> Array DIM2 Double+ +slurpDoublesFromImage mkDouble arrBound+ = traverse arrBound+ (\(Z :. height :. width :. _) + -> Z :. height :. width)++ (\get (Z :. y :. x)+ -> mkDouble+ (get (Z :. y :. x :. 0))+ (get (Z :. y :. x :. 1))+ (get (Z :. y :. x :. 2)))+++makeImageFromDoubles+ :: (Double -> (Double, Double, Double))+ -> Array DIM2 Double+ -> Array DIM3 Word8+ +makeImageFromDoubles fnColor arrDoubles+ = traverse arrDoubles+ (\(Z :. height :. width)+ -> Z :. height :. width :. 4)+ + (\get (Z :. y :. x :. c)+ -> let (r, g, b) = rampColorHotToCold 0 1 (get (Z :. y :. x))+ in case c of+ 0 -> truncate (r * 255)+ 1 -> truncate (g * 255)+ 2 -> truncate (b * 255)+ 3 -> 0)+++-- | Extract the boundary value from a RGB triple.+slurpBoundValue :: Word8 -> Word8 -> Word8 -> Double+slurpBoundValue r g b -- A non-boundary value. | r == 0 && g == 0 && b == 255 - = (1, 0)+ = 0 -- A boundary value. | (r == g) && (r == b) - = (0, (fromIntegral r) / 255)+ = fromIntegral r / 255+ + | otherwise+ = error $ "Unhandled pixel value in input " ++ show (r, g, b)+++-- | Extract boundary mask from a RGB triple.+slurpBoundMask :: Word8 -> Word8 -> Word8 -> Double+slurpBoundMask r g b+ -- A non-boundary value.+ | r == 0 && g == 0 && b == 255 + = 1++ -- A boundary value.+ | (r == g) && (r == b) + = 0 | otherwise = error $ "Unhandled pixel value in input " ++ show (r, g, b)
MMult/src/Main.hs view
@@ -1,11 +1,13 @@ {-# LANGUAGE PatternGuards #-} -import Solver-import Data.Array.Repa +import Data.Array.Repa as A+import Data.Array.Repa.IO.Matrix+import Data.Array.Repa.Algorithms.Matrix import Data.Maybe-import Matrix import System.Environment import Control.Monad+import System.Random+import qualified Data.Array.Parallel.Unlifted as U -- Arg Parsing ------------------------------------------------------------------------------------ data Arg@@ -64,6 +66,28 @@ ArgMatrixRandom height width -> genRandomMatrix (Z :. height :. width) ++-- Random -----------------------------------------------------------------------------------------+-- | Generate a random(ish) matrix.+genRandomMatrix + :: DIM2 + -> IO (Array DIM2 Double)++genRandomMatrix sh+ = do uarr <- genRandomUArray (A.size sh)+ return $ fromUArray sh uarr++-- | Generate a random(ish) UArray of doubles.+-- The std random function is too slow to generate really big vectors+-- with. Instead, we generate a short random vector and repeat that.+genRandomUArray :: Int -> IO (U.Array Double)+genRandomUArray n + = do let k = 1000+ rg <- newStdGen+ let randvec = U.randomRs k (-100, 100) rg+ let vec = U.map (\i -> randvec U.!: (i `mod` k)) (U.enumFromTo 0 (n-1))+ return vec+ -- Main ------------------------------------------------------------------------------------------- main :: IO ()@@ -85,18 +109,16 @@ `deepSeqArray` return () -- Run the solver.- let matResult = mmMult mat1 mat2+ let matResult = multiplyMM mat1 mat2 matResult `deepSeqArray` return () -- Write the output to file if requested. case mArgOut of Nothing -> return ()- Just fileOut -> writeMatrixAsTextFile matResult fileOut+ Just fileOut -> writeMatrixToTextFile fileOut matResult | otherwise = printHelp --
− MMult/src/Solver.hs
@@ -1,14 +0,0 @@-{-# LANGUAGE TypeOperators #-}--module Solver where-import Data.Array.Repa--mmMult :: Array DIM2 Double- -> Array DIM2 Double- -> Array DIM2 Double--mmMult arr1@Manifest{} arr2@Manifest{}- = fold (+) 0 - $ traverse2 arr1 (force $ transpose arr2)- (\(sh :. m1 :. n1) -> \(_ :. n2 :. m2) -> (sh :. m1 :. n2 :. n1))- (\f1 -> \f2 -> \(sh :. i :. j :. k) -> f1 (sh :. i :. k) * f2 (sh :. j :. k))
− lib/ColorRamp.hs
@@ -1,47 +0,0 @@--module ColorRamp- (rampColorHotToCold)-where----- Color Ramps -------------------------------------------------------------------------------------- | Standard Hot -> Cold hypsometric color ramp.--- Sequence is red, yellow, green, cyan, blue.-rampColorHotToCold - :: (Ord a, Floating a) - => a - -> a - -> a - -> (a, a, a)- -{-# INLINE rampColorHotToCold #-}-rampColorHotToCold vmin vmax vNotNorm- = let - v | vNotNorm < vmin = vmin- | vNotNorm > vmax = vmax- | otherwise = vNotNorm- - dv = vmax - vmin -- result | v < vmin + 0.25 * dv- = ( 0- , 4 * (v - vmin) / dv- , 1.0)- - | v < vmin + 0.5 * dv- = ( 0- , 1.0- , 1 + 4 * (vmin + 0.25 * dv - v) / dv)- - | v < vmin + 0.75 * dv- = ( 4 * (v - vmin - 0.5 * dv) / dv- , 1.0- , 0.0)- - | otherwise- = ( 1.0- , 1 + 4 * (vmin + 0.75 * dv - v) / dv- , 0)- - in result-
− lib/Matrix.hs
@@ -1,86 +0,0 @@---- | Reading and writing matricies as ASCII files.--- We use ASCII so we can generate and check simple test data by hand,--- and we don't want to fool around with byte order issues.------ Matrix file format is like:------ MATRIX -- header--- 100 100 -- width and height--- 1.23 1.56 1.23 ... -- data, separated by whitespace--- ....------ TODO: Merge this with PPM.hs--- TODO: Merge this with Vector.hs. --- We should really have fns that read and write arrays of arbitrary dimension.----module Matrix- ( readMatrixFromTextFile- , writeMatrixAsTextFile- , genRandomMatrix)-where-import Data.List as L-import Data.Array.Repa as A-import qualified Data.Array.Parallel.Unlifted as U-import Prelude as P-import System.IO-import Control.Monad-import Data.Char-import System.Random-import Vector----- Reading ------------------------------------------------------------------------------------------- | Read a matrix from a text file.-readMatrixFromTextFile- :: (U.Elt a, Num a, Read a)- => FilePath -- ^ File name of matrix file.- -> IO (Array DIM2 a) --readMatrixFromTextFile fileName- = do handle <- openFile fileName ReadMode- - "MATRIX" <- hGetLine handle- [width, height] <- liftM (P.map read . words) $ hGetLine handle- str <- hGetContents handle- let vals = readValues str-- let dims = Z :. width :. height- let mat = fromList dims vals-- return mat------ | Write a matrix as a text file.-writeMatrixAsTextFile - :: (U.Elt a, Show a)- => Array DIM2 a -- ^ Matrix to write.- -> FilePath -- ^ File name of output file.- -> IO ()--writeMatrixAsTextFile arr fileName- = do file <- openFile fileName WriteMode -- hPutStrLn file "MATRIX"-- let Z :. width :. height - = extent arr-- hPutStrLn file $ show width ++ " " ++ show height- - hWriteValues file $ toList arr- hClose file----- | Generate a random(ish) matrix.-genRandomMatrix - :: DIM2 - -> IO (Array DIM2 Double)--genRandomMatrix sh- = do uarr <- genRandomUArray (A.size sh)- return $ fromUArray sh uarr---
− lib/PPM.hs
@@ -1,171 +0,0 @@-{-# LANGUAGE ScopedTypeVariables, BangPatterns #-}---- | Writing out matricies as PPM image files.-module PPM- ( writeMatrixAsNormalisedPPM- , writeMatrixAsPPM- , readPPMAsMatrix- , readPPMAsMatrix2)-where-import qualified Data.Array.Parallel.Unlifted as U-import Data.List as L-import Data.Array.Repa as A-import Prelude as P-import System.IO-import Control.Monad-import Data.Char----- Write --------------------------------------------------------------------------------------------- | Convert a matrix to a PPM image,--- while normalising it to the maximum value present in the matrix.-writeMatrixAsNormalisedPPM - :: FilePath -- ^ Filename of output file.- -> (Double -> (Double, Double, Double)) -- ^ Function for producing colors from data values,- -> Array DIM2 Double -- ^ Matrix of values (need not be normalised).- -> IO ()--writeMatrixAsNormalisedPPM fileName colorFn arr- = let -- Use the maximum elem in the array as the white value.- vals = U.toList $ toUArray arr- maxVal = maximum vals-- -- Normalise the array to the range [0..1] for display.- arrNorm = A.map (/ maxVal) arr-- in writeMatrixAsPPM fileName colorFn arrNorm----- | Convert a matrix to a PPM image.--- Matrix elements should be normalised to [0..1]-writeMatrixAsPPM - :: FilePath- -> (Double -> (Double, Double, Double)) -- ^ Function for producing colors from data values.- -> Array DIM2 Double -- ^ Matrix of values, normalised to [0..1]- -> IO ()--writeMatrixAsPPM fileName colorFn arr- = let - -- Break flat array data into individual rows- Z :. width :. height - = extent arr-- in do- file <- openFile fileName WriteMode- hPutStrLn file $ "P3"- hPutStrLn file $ show width ++ " " ++ show height- hPutStrLn file $ "255"- - hWritePixels file colorFn $ toList arr- hClose file----- | Write out pixel values to a file.-hWritePixels - :: Handle - -> (Double -> (Double, Double, Double)) -- ^ Function for producing colors from data values.- -> [Double] -- ^ Data values.- -> IO ()--hWritePixels h colorFn !xx- = go xx- where- go [] = return ()- go (x:xs)- = do let (r, g, b) = colorFn x- hPutStr h $ showInt $ truncate (r * 255)- hPutStr h $ " "- hPutStr h $ showInt $ truncate (g * 255)- hPutStr h $ " "- hPutStr h $ showInt $ truncate (b * 255)- hPutStr h $ "\n"- go xs- -showInt :: Int -> String-showInt i = show i----- Read --------------------------------------------------------------------------------------------readPPMAsMatrix - :: (Int -> Int -> Int -> Double) -- ^ Function for producing array values from RGB pixel values- -> FilePath -- ^ File name of ppm file.- -> IO (Array DIM2 Double) -- ^ Loaded matrix.- -readPPMAsMatrix pointFn fileName- = do file <- openFile fileName ReadMode- - "P3" <- hGetLine file- [width, height] <- liftM (P.map read . words) $ hGetLine file- _maxVal <- hGetLine file- vals <- loadPixels pointFn file-- let dims = Z :. width :. height- let mat = fromList dims vals-- return mat----- | Read the values of two separate matricies that are encoded in a single PPM file.-readPPMAsMatrix2 - :: (Int -> Int -> Int -> (Double, Double)) - -- ^ Function for producing array values from RGB pixel values- -> FilePath -- ^ File name of ppm file.- -> IO ( Array DIM2 Double -- ^ Loaded matrix.- , Array DIM2 Double)- -readPPMAsMatrix2 pointFn fileName- = do file <- openFile fileName ReadMode- - "P3" <- hGetLine file- [width, height] <- liftM (P.map read . words) $ hGetLine file- _maxVal <- hGetLine file- vals <- loadPixels pointFn file-- let dims = Z :. width :. height- let mat1 = fromList dims $ P.map fst vals- let mat2 = fromList dims $ P.map snd vals-- return (mat1, mat2)- ---- | Load list of array data from the data part of a PPM file.-loadPixels - :: (Int -> Int -> Int -> a) -- ^ Function for producing array values from RGB pixel values.- -> Handle -- ^ Handle of file to load from.- -> IO [a]- -loadPixels pointFn handle- = do str <- hGetContents handle- let ints = readInts str- let vals = convertLine pointFn ints- return vals----- | Convert groups of RGB components to values.-convertLine - :: (Int -> Int -> Int -> a) -- ^ Function for producing array values from RGB pixel values.- -> [Int] -- ^ Ints read from the PPM file.- -> [a] -- ^ Output values.--convertLine fn vs- = case vs of- [] -> []- r : g : b : rest -> fn r g b : convertLine fn rest- _ -> error "PPM.convertLine: bad pixel"- --- | Read a string containing ints separated by whitespace. -readInts :: String -> [Int]-readInts cs = readInts' [] cs- where readInts' _ [] = []- readInts' acc (c : rest)- | isSpace c- = if null acc - then readInts' [] rest- else read (reverse acc) : readInts' [] rest-- | isDigit c- = readInts' (c : acc) rest-- | otherwise- = error $ "unexpected char in PPM file " ++ show (ord c)-
− lib/Vector.hs
@@ -1,137 +0,0 @@---- | Reading and writing vectors as ASCII files.--- We use ASCII so we can generate and check simple test data by hand,--- and we don't want to fool around with byte order issues.------ Vector file format is like:------ VECTOR -- header--- 100 -- length of vector--- 1.23 1.56 1.23 ... -- data, separated by whitespace--- ....-------module Vector- ( readVectorFromTextFile- , writeVectorAsTextFile- , genRandomVector-- , readValues- , hWriteValues- , genRandomUArray)-where-import Data.List as L-import Data.Array.Repa as A-import qualified Data.Array.Parallel.Unlifted as U-import Prelude as P-import System.IO-import Control.Monad-import Data.Char-import System.Random----- | Read a vector from a text file.-readVectorFromTextFile- :: (U.Elt a, Num a, Read a)- => FilePath -- ^ File name of vector file.- -> IO (Array DIM1 a) --readVectorFromTextFile fileName- = do handle <- openFile fileName ReadMode- - "VECTOR" <- hGetLine handle- [len] <- liftM (P.map readInt . words) $ hGetLine handle- str <- hGetContents handle- let vals = readValues str-- let dims = Z :. len- let vec = fromList dims vals-- return vec--readInt :: String -> Int-readInt str- | and $ P.map isDigit str- = read str- - --- | Write a vector as a text file.-writeVectorAsTextFile - :: (U.Elt a, Show a)- => Array DIM1 a -- ^ Vector to write.- -> FilePath -- ^ File name of output file.- -> IO ()--writeVectorAsTextFile arr fileName- = do file <- openFile fileName WriteMode -- hPutStrLn file "VECTOR"-- let Z :. len- = extent arr-- hPutStrLn file $ show len- hWriteValues file $ toList arr- hClose file- ---- | Generate a random(ish) vector.-genRandomVector - :: DIM2 - -> IO (Array DIM2 Double)--genRandomVector sh- = do uarr <- genRandomUArray (A.size sh)- return $ fromUArray sh uarr- ----- Stuff shared with Matrix module ---------------------------------------------------------------- | Write out values to a file.-hWriteValues- :: Show a- => Handle - -> [a] -- ^ Data values.- -> IO ()--hWriteValues handle xx- = go xx- where- go [] = return ()- go (x:xs)- = do hPutStr handle $ show x- hPutStr handle $ "\n"- go xs----- | Read a string containing ints separated by whitespace. -readValues :: (Num a, Read a) => String -> [a]-readValues cs = readValues' [] cs- where readValues' _ [] = []- readValues' acc (c : rest)- | isSpace c- = if null acc - then readValues' [] rest- else read (reverse acc) : readValues' [] rest-- | isDigit c || c == '.' || c == 'e' || c == '-'- = readValues' (c : acc) rest-- | otherwise- = error $ "unexpected char in Matrix file " ++ show (ord c)----- | Generate a random(ish) UArray of doubles.--- The std random function is too slow to generate really big vectors--- with. Instead, we generate a short random vector and repeat that.-genRandomUArray :: Int -> IO (U.Array Double)-genRandomUArray n - = do let k = 1000- rg <- newStdGen- let randvec = U.randomRs k (-100, 100) rg- let vec = U.map (\i -> randvec U.!: (i `mod` k)) (U.enumFromTo 0 (n-1))- return vec----
repa-examples.cabal view
@@ -1,6 +1,6 @@ Name: repa-examples-Version: 1.0.0.0-License: MIT+Version: 1.1.0.0+License: BSD3 License-file: LICENSE Author: The DPH Team Maintainer: Ben Lippmeier <benl@ouroborus.net>@@ -16,44 +16,59 @@ Synopsis: Examples using the Repa array library. --Executable repa-laplace+Executable repa-fft2d Build-depends: - base >= 3 && < 5, - random >= 1.0 && < 1.1,- dph-prim-par >= 0.4 && < 0.5,- repa >= 1.0 && < 1.1+ base == 4.*,+ dph-prim-par == 0.4.*,+ dph-base == 0.4.*,+ repa == 1.1.*,+ repa-algorithms == 1.1.*,+ repa-io == 1.1.* - Main-is: Laplace/src/Main.hs- other-modules: Solver Matrix Vector ColorRamp PPM- hs-source-dirs: Laplace/src . lib+ Main-is: FFT/src/FFT2d/Main.hs+ hs-source-dirs: FFT/src . ghc-options: -Odph -threaded -Executable repa-mmult+Executable repa-fft2d-highpass Build-depends: - base >= 3 && < 5, - random >= 1.0 && < 1.1,- dph-prim-par >= 0.4 && < 0.5,- repa >= 1.0 && < 1.1+ base == 4.*,+ dph-prim-par == 0.4.*,+ dph-base == 0.4.*,+ repa == 1.1.*,+ repa-algorithms == 1.1.*,+ repa-io == 1.1.* - Main-is: MMult/src/Main.hs- other-modules: Solver Matrix Vector ColorRamp PPM- hs-source-dirs: MMult/src . lib+ Main-is: FFT/src/HighPass/Main.hs+ hs-source-dirs: FFT/src . ghc-options: -Odph -threaded -Executable repa-fft+Executable repa-laplace Build-depends: - base >= 3 && < 5, - random >= 1.0 && < 1.1,- dph-prim-par >= 0.4 && < 0.5,- dph-base >= 0.4 && < 0.5,- repa >= 1.0 && < 1.1+ base == 4.*,+ dph-prim-par == 0.4.*,+ dph-base == 0.4.*,+ repa == 1.1.*,+ repa-io == 1.1.* - Main-is: FFT/src/Main.hs- other-modules: FFT DFT Roots StrictComplex Matrix Vector ColorRamp PPM- hs-source-dirs: FFT/src . lib+ Main-is: Laplace/src/Main.hs+ other-modules: Solver+ hs-source-dirs: Laplace/src . ghc-options: -Odph -threaded ++Executable repa-mmult+ Build-depends: + base == 4.*,+ dph-prim-par == 0.4.*,+ dph-base == 0.4.*,+ repa == 1.1.*,+ repa-io == 1.1.*,+ repa-algorithms == 1.1.*,+ random == 1.0.*++ Main-is: MMult/src/Main.hs+ hs-source-dirs: MMult/src .+ ghc-options: -Odph -threaded