diff --git a/FFT/src/DFT.hs b/FFT/src/DFT.hs new file mode 100644 --- /dev/null +++ b/FFT/src/DFT.hs @@ -0,0 +1,61 @@ +{-# 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) + diff --git a/FFT/src/FFT.hs b/FFT/src/FFT.hs new file mode 100644 --- /dev/null +++ b/FFT/src/FFT.hs @@ -0,0 +1,136 @@ +{-# 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' + diff --git a/FFT/src/Main.hs b/FFT/src/Main.hs new file mode 100644 --- /dev/null +++ b/FFT/src/Main.hs @@ -0,0 +1,194 @@ +{-# 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 Read a real valued input vector from this file." + , " -vector-real-step Use a real valued step function for input." + , " -ppm Use a PPM file for input." + , "" + , "OUTPUT:" + , " -out-vector-magnitude Write the magnitute of the transformed vector to file." + , " -out-ppm-magnitude Write transformed matrix to a ppm file." + , " -out-ppm-clip ... 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) + diff --git a/FFT/src/Roots.hs b/FFT/src/Roots.hs new file mode 100644 --- /dev/null +++ b/FFT/src/Roots.hs @@ -0,0 +1,43 @@ +{-# 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 diff --git a/FFT/src/StrictComplex.hs b/FFT/src/StrictComplex.hs new file mode 100644 --- /dev/null +++ b/FFT/src/StrictComplex.hs @@ -0,0 +1,32 @@ +{-# 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) + + diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,13 @@ +Copyright (c) 2010 The DPH Team + + 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: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. diff --git a/Laplace/src/Main.hs b/Laplace/src/Main.hs new file mode 100644 --- /dev/null +++ b/Laplace/src/Main.hs @@ -0,0 +1,91 @@ +{-# 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" +-- +import Solver +import Data.Array.Repa as A +import ColorRamp +import PPM +import System.Environment + +main :: IO () +main + = do args <- getArgs + case args of + [steps, fileInput, fileOutput] + -> laplace (read steps) fileInput fileOutput + + _ -> do + putStr usage + return () + + +-- | Command line usage information. +usage :: String +usage = unlines + [ "Usage: laplace " + , "" + , " 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." + , "" + , " Format of input file:" + , " Boundary values are indicated in greyscale," + , " ie from the list [(x, x, x) | x <- [0 .. 255]]" + , " Non-boundary values are indicated in blue," + , " ie (0, 0, 255)" + , " Any other pixel value is an error." + , "" + ] + + +-- | Solve it. +laplace :: Int -- ^ Number of iterations to use. + -> FilePath -- ^ Input file. + -> FilePath -- ^ Output file + -> IO () + +laplace steps fileInput fileOutput + = do + -- Load up the file containing boundary conditions. + (matBoundMask, matBoundValue) + <- readPPMAsMatrix2 loadPixel fileInput + + -- Use the boundary condition values as the initial matrix. + let matInitial = matBoundValue + + -- Run the solver. + let matFinal = solveLaplace + steps + matBoundMask + matBoundValue + matInitial + + matFinal `deepSeqArray` return () + + -- Write out the matrix as a colorised PPM image + writeMatrixAsNormalisedPPM + fileOutput + (rampColorHotToCold 0.0 1.0) + matFinal + + +-- | Extract boundary mask and value from a pixel value. +loadPixel :: Int -> Int -> Int -> (Double, Double) +loadPixel r g b + -- A non-boundary value. + | r == 0 && g == 0 && b == 255 + = (1, 0) + + -- A boundary value. + | (r == g) && (r == b) + = (0, (fromIntegral r) / 255) + + | otherwise + = error $ "Unhandled pixel value in input " ++ show (r, g, b) + + diff --git a/Laplace/src/Solver.hs b/Laplace/src/Solver.hs new file mode 100644 --- /dev/null +++ b/Laplace/src/Solver.hs @@ -0,0 +1,82 @@ +{-# LANGUAGE BangPatterns #-} +module Solver + (solveLaplace) +where +import Data.Array.Repa as A +import qualified Data.Array.Repa.Shape as S +import Data.Array.Parallel.Unlifted (Elt) +import qualified Data.Array.Parallel.Unlifted as U + + +-- | Solver for the Laplace equation. +solveLaplace + :: Int -- ^ Number of iterations to use. + -> Array DIM2 Double -- ^ Boundary value mask. + -> Array DIM2 Double -- ^ Boundary values. + -> Array DIM2 Double -- ^ Initial state. + -> Array DIM2 Double + +{-# INLINE solveLaplace #-} +solveLaplace steps arrBoundMask@Manifest{} arrBoundValue@Manifest{} arrInit + = go steps arrInit + where go i arr@Manifest{} + | i == 0 = arr + | otherwise = go (i - 1) + $ force + $ applyBoundary arrBoundMask arrBoundValue + $ relaxLaplace arr + + +-- | Perform matrix relaxation for the Laplace equation, +-- using a stencil function. +-- +-- Computation fn is +-- u'(i,j) = (u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1)) / 4 +-- +relaxLaplace + :: Array DIM2 Double + -> Array DIM2 Double + +{-# INLINE relaxLaplace #-} +relaxLaplace !arr + = traverse arr id elemFn + where + _ :. height :. width + = extent arr + + {-# INLINE elemFn #-} + elemFn get d@(sh :. i :. j) + = if isBorder i j + then get d + else (get (sh :. (i-1) :. j) + + get (sh :. i :. (j-1)) + + get (sh :. (i+1) :. j) + + get (sh :. i :. (j+1))) / 4 + + -- Check if this element is on the border of the matrix. + -- If so we can't apply the stencil because we don't have all the neighbours. + {-# INLINE isBorder #-} + isBorder i j + = (i == 0) || (i >= width - 1) + || (j == 0) || (j >= height - 1) + + +-- | Apply the boundary conditions to this matrix. +-- The mask matrix has 0 in places where boundary conditions hold +-- and 1 otherwise. +-- +-- The value matrix has the boundary condition value in places where it holds, +-- and 0 otherwise. +-- +applyBoundary + :: Array DIM2 Double -- ^ Boundary condition mask. + -> Array DIM2 Double -- ^ Boundary condition values. + -> Array DIM2 Double -- ^ Initial matrix. + -> Array DIM2 Double -- ^ Matrix with boundary conditions applied. + +{-# INLINE applyBoundary #-} +applyBoundary arrBoundMask arrBoundValue arr + = A.zipWith (+) arrBoundValue + $ A.zipWith (*) arrBoundMask arr + + diff --git a/MMult/src/Main.hs b/MMult/src/Main.hs new file mode 100644 --- /dev/null +++ b/MMult/src/Main.hs @@ -0,0 +1,102 @@ +{-# LANGUAGE PatternGuards #-} + +import Solver +import Data.Array.Repa +import Data.Maybe +import Matrix +import System.Environment +import Control.Monad + +-- Arg Parsing ------------------------------------------------------------------------------------ +data Arg + = ArgSolver String + | ArgMatrixRandom Int Int + | ArgMatrixFile FilePath + | ArgOutFile FilePath + deriving Show + +isArgMatrix arg + = case arg of + ArgMatrixRandom{} -> True + ArgMatrixFile{} -> True + _ -> False + +parseArgs [] = [] +parseArgs (flag:xx) + | "-file" <- flag + , f:rest <- xx + = ArgMatrixFile f : parseArgs rest + + | "-out" <- flag + , f:rest <- xx + = ArgOutFile f : parseArgs rest + + | "-random" <- flag + , x:y:rest <- xx + = ArgMatrixRandom (read x) (read y) : parseArgs rest + + | otherwise + = error $ "bad arg " ++ flag ++ "\n" + +printHelp + = putStr + $ unlines + [ "Usage: mmult [args..]" + , "" + , " -random Use a random matrix of this size." + , " -file Read a matrix from this file." + , " -out Write resulting matrix to this file." + , "" + , " Format of matrix file:" + , " MATRIX" + , " " + , " " + , "" ] + + +-- | Get a matrix from a file, or generate a random one. +getMatrix :: Arg -> IO (Array DIM2 Double) +getMatrix arg + = case arg of + ArgMatrixFile fileName + -> readMatrixFromTextFile fileName + + ArgMatrixRandom height width + -> genRandomMatrix (Z :. height :. width) + + +-- Main ------------------------------------------------------------------------------------------- +main :: IO () +main + = do args <- liftM parseArgs $ getArgs + main' args + +main' args + | [argMat1, argMat2] <- filter isArgMatrix args + , mArgOut <- listToMaybe [s | ArgOutFile s <- args] + = do + -- Get matrices from files, + -- or generate random ones we were asked to. + mat1 <- getMatrix argMat1 + mat2 <- getMatrix argMat2 + + mat1 + `deepSeqArray` mat2 + `deepSeqArray` return () + + -- Run the solver. + let matResult = mmMult mat1 mat2 + + matResult `deepSeqArray` return () + + -- Write the output to file if requested. + case mArgOut of + Nothing -> return () + Just fileOut -> writeMatrixAsTextFile matResult fileOut + + | otherwise + = printHelp + + + + diff --git a/MMult/src/Solver.hs b/MMult/src/Solver.hs new file mode 100644 --- /dev/null +++ b/MMult/src/Solver.hs @@ -0,0 +1,14 @@ +{-# 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)) diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/lib/ColorRamp.hs b/lib/ColorRamp.hs new file mode 100644 --- /dev/null +++ b/lib/ColorRamp.hs @@ -0,0 +1,47 @@ + +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 + diff --git a/lib/Matrix.hs b/lib/Matrix.hs new file mode 100644 --- /dev/null +++ b/lib/Matrix.hs @@ -0,0 +1,86 @@ + +-- | 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 + + + diff --git a/lib/PPM.hs b/lib/PPM.hs new file mode 100644 --- /dev/null +++ b/lib/PPM.hs @@ -0,0 +1,171 @@ +{-# 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) + diff --git a/lib/Vector.hs b/lib/Vector.hs new file mode 100644 --- /dev/null +++ b/lib/Vector.hs @@ -0,0 +1,137 @@ + +-- | 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 + + + + diff --git a/repa-examples.cabal b/repa-examples.cabal new file mode 100644 --- /dev/null +++ b/repa-examples.cabal @@ -0,0 +1,59 @@ +Name: repa-examples +Version: 1.0.0.0 +License: MIT +License-file: LICENSE +Author: The DPH Team +Maintainer: Ben Lippmeier +Build-Type: Simple +Cabal-Version: >=1.6 +Stability: experimental +Category: Data Structures +Homepage: http://trac.haskell.org/repa +Bug-reports: http://trac.haskell.org/repa/newticket +Description: + Examples using the Repa array library. + +Synopsis: + Examples using the Repa array library. + + +Executable repa-laplace + Build-depends: + base >= 3 && < 5, + random >= 1.0 && < 1.1, + dph-prim-par >= 0.4 && < 0.5, + repa >= 1.0 && < 1.1 + + Main-is: Laplace/src/Main.hs + other-modules: Solver Matrix Vector ColorRamp PPM + hs-source-dirs: Laplace/src . lib + ghc-options: -Odph -threaded + + +Executable repa-mmult + Build-depends: + base >= 3 && < 5, + random >= 1.0 && < 1.1, + dph-prim-par >= 0.4 && < 0.5, + repa >= 1.0 && < 1.1 + + Main-is: MMult/src/Main.hs + other-modules: Solver Matrix Vector ColorRamp PPM + hs-source-dirs: MMult/src . lib + ghc-options: -Odph -threaded + + +Executable repa-fft + 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 + + Main-is: FFT/src/Main.hs + other-modules: FFT DFT Roots StrictComplex Matrix Vector ColorRamp PPM + hs-source-dirs: FFT/src . lib + ghc-options: -Odph -threaded + +