Update Repa libraries (fixes Blur for now)
authorDavid Terei <davidterei@gmail.com>
Mon, 16 Jan 2012 09:13:48 +0000 (01:13 -0800)
committerDavid Terei <davidterei@gmail.com>
Mon, 16 Jan 2012 09:14:58 +0000 (01:14 -0800)
80 files changed:
fibon/Repa/Blur/Fibon/Instance.hs [deleted file]
fibon/Repa/Blur/Fibon/data/ref/input/lena.bmp [deleted file]
fibon/Repa/Blur/Fibon/data/ref/output/out.expected.bmp [deleted file]
fibon/Repa/Blur/Fibon/data/test/input/step20.bmp [deleted file]
fibon/Repa/Blur/Fibon/data/test/output/out.expected.bmp [deleted file]
fibon/Repa/Blur/Makefile
fibon/Repa/Blur/src/Main.hs
fibon/Repa/_RepaLib/bmp/Codec/BMP.hs
fibon/Repa/_RepaLib/bmp/Codec/BMP/Base.hs
fibon/Repa/_RepaLib/bmp/Codec/BMP/BitmapInfo.hs
fibon/Repa/_RepaLib/bmp/Codec/BMP/BitmapInfoV3.hs
fibon/Repa/_RepaLib/bmp/Codec/BMP/BitmapInfoV4.hs
fibon/Repa/_RepaLib/bmp/Codec/BMP/BitmapInfoV5.hs
fibon/Repa/_RepaLib/bmp/Codec/BMP/Compression.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/bmp/Codec/BMP/Error.hs
fibon/Repa/_RepaLib/bmp/Codec/BMP/FileHeader.hs
fibon/Repa/_RepaLib/bmp/Codec/BMP/Pack.hs
fibon/Repa/_RepaLib/bmp/Codec/BMP/Unpack.hs
fibon/Repa/_RepaLib/bmp/bmp.cabal
fibon/Repa/_RepaLib/quickcheck/LICENSE
fibon/Repa/_RepaLib/quickcheck/QuickCheck.cabal
fibon/Repa/_RepaLib/quickcheck/README
fibon/Repa/_RepaLib/quickcheck/Test/QuickCheck.hs
fibon/Repa/_RepaLib/quickcheck/Test/QuickCheck/All.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/quickcheck/Test/QuickCheck/Arbitrary.hs
fibon/Repa/_RepaLib/quickcheck/Test/QuickCheck/Exception.hs
fibon/Repa/_RepaLib/quickcheck/Test/QuickCheck/Function.hs
fibon/Repa/_RepaLib/quickcheck/Test/QuickCheck/Gen.hs
fibon/Repa/_RepaLib/quickcheck/Test/QuickCheck/Modifiers.hs
fibon/Repa/_RepaLib/quickcheck/Test/QuickCheck/Monadic.hs
fibon/Repa/_RepaLib/quickcheck/Test/QuickCheck/Poly.hs
fibon/Repa/_RepaLib/quickcheck/Test/QuickCheck/Property.hs
fibon/Repa/_RepaLib/quickcheck/Test/QuickCheck/State.hs
fibon/Repa/_RepaLib/quickcheck/Test/QuickCheck/Test.hs
fibon/Repa/_RepaLib/quickcheck/Test/QuickCheck/Text.hs
fibon/Repa/_RepaLib/repa-algorithms/Data/Array/Repa/Algorithms/Complex.hs
fibon/Repa/_RepaLib/repa-algorithms/Data/Array/Repa/Algorithms/Convolve.hs
fibon/Repa/_RepaLib/repa-algorithms/Data/Array/Repa/Algorithms/DFT.hs
fibon/Repa/_RepaLib/repa-algorithms/Data/Array/Repa/Algorithms/FFT.hs
fibon/Repa/_RepaLib/repa-algorithms/Data/Array/Repa/Algorithms/Iterate.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/repa-algorithms/Data/Array/Repa/Algorithms/Matrix.hs
fibon/Repa/_RepaLib/repa-algorithms/Data/Array/Repa/Algorithms/Randomish.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/repa-algorithms/repa-algorithms.cabal
fibon/Repa/_RepaLib/repa-bytestring/Data/Array/Repa/ByteString.hs
fibon/Repa/_RepaLib/repa-bytestring/repa-bytestring.cabal
fibon/Repa/_RepaLib/repa-io/Data/Array/Repa/IO/BMP.hs
fibon/Repa/_RepaLib/repa-io/Data/Array/Repa/IO/Binary.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/repa-io/Data/Array/Repa/IO/Matrix.hs
fibon/Repa/_RepaLib/repa-io/Data/Array/Repa/IO/Timing.hs
fibon/Repa/_RepaLib/repa-io/Data/Array/Repa/IO/Vector.hs
fibon/Repa/_RepaLib/repa-io/repa-io.cabal
fibon/Repa/_RepaLib/repa/Data/Array/Repa.hs
fibon/Repa/_RepaLib/repa/Data/Array/Repa/Arbitrary.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/repa/Data/Array/Repa/Index.hs
fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/Base.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/Elt.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/EvalBlockwise.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/EvalChunked.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/EvalCursored.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/EvalReduction.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/Forcing.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/Gang.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/Select.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/repa/Data/Array/Repa/Operators/IndexSpace.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/repa/Data/Array/Repa/Operators/Interleave.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/repa/Data/Array/Repa/Operators/Mapping.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/repa/Data/Array/Repa/Operators/Modify.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/repa/Data/Array/Repa/Operators/Reduction.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/repa/Data/Array/Repa/Operators/Select.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/repa/Data/Array/Repa/Operators/Traverse.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/repa/Data/Array/Repa/Properties.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/repa/Data/Array/Repa/QuickCheck.hs [deleted file]
fibon/Repa/_RepaLib/repa/Data/Array/Repa/Shape.hs
fibon/Repa/_RepaLib/repa/Data/Array/Repa/Slice.hs
fibon/Repa/_RepaLib/repa/Data/Array/Repa/Specialised/Dim2.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/repa/Data/Array/Repa/Stencil.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/repa/Data/Array/Repa/Stencil/Base.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/repa/Data/Array/Repa/Stencil/Template.hs [new file with mode: 0644]
fibon/Repa/_RepaLib/repa/Test.hs [deleted file]
fibon/Repa/_RepaLib/repa/repa.cabal

diff --git a/fibon/Repa/Blur/Fibon/Instance.hs b/fibon/Repa/Blur/Fibon/Instance.hs
deleted file mode 100644 (file)
index 2a546b7..0000000
+++ /dev/null
@@ -1,26 +0,0 @@
-{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
-module Fibon.Benchmarks.Repa.Blur.Fibon.Instance(
-  mkInstance
-)
-where
-import Fibon.BenchmarkInstance
-
-sharedConfig = BenchmarkInstance {
-    flagConfig = FlagConfig {
-        configureFlags = ["--ghc-option=-threaded"]
-      , buildFlags     = []
-      , runFlags       = []
-      }
-    , stdinInput     = Nothing
-    , output         = [(OutputFile "out.bmp", Diff "out.expected.bmp")]
-    , exeName        = "repa-blur"
-  }
-flgCfg = flagConfig sharedConfig
-
-mkInstance Test = sharedConfig {
-        flagConfig = flgCfg {runFlags = words "1 step20.bmp out.bmp"}
-    }
-mkInstance Ref  = sharedConfig {
-        flagConfig = flgCfg {runFlags = words "12 lena.bmp out.bmp"}
-    }
-
diff --git a/fibon/Repa/Blur/Fibon/data/ref/input/lena.bmp b/fibon/Repa/Blur/Fibon/data/ref/input/lena.bmp
deleted file mode 100644 (file)
index 53d59e9..0000000
Binary files a/fibon/Repa/Blur/Fibon/data/ref/input/lena.bmp and /dev/null differ
diff --git a/fibon/Repa/Blur/Fibon/data/ref/output/out.expected.bmp b/fibon/Repa/Blur/Fibon/data/ref/output/out.expected.bmp
deleted file mode 100644 (file)
index 71fe990..0000000
Binary files a/fibon/Repa/Blur/Fibon/data/ref/output/out.expected.bmp and /dev/null differ
diff --git a/fibon/Repa/Blur/Fibon/data/test/input/step20.bmp b/fibon/Repa/Blur/Fibon/data/test/input/step20.bmp
deleted file mode 100644 (file)
index 0c4ca29..0000000
Binary files a/fibon/Repa/Blur/Fibon/data/test/input/step20.bmp and /dev/null differ
diff --git a/fibon/Repa/Blur/Fibon/data/test/output/out.expected.bmp b/fibon/Repa/Blur/Fibon/data/test/output/out.expected.bmp
deleted file mode 100644 (file)
index 26a2796..0000000
Binary files a/fibon/Repa/Blur/Fibon/data/test/output/out.expected.bmp and /dev/null differ
index 2c1bd54..ca3cd78 100644 (file)
@@ -1,36 +1,76 @@
 TOP = ../../..
 include $(TOP)/mk/boilerplate.mk
-SRCS = ../_RepaLib/bmp/Codec/BMP/CIEXYZ.hs \
-       ../_RepaLib/bmp/Codec/BMP/Error.hs \
+SRCS = ../_RepaLib/bmp/Codec/BMP/Base.hs \
+       ../_RepaLib/bmp/Codec/BMP/BitmapInfo.hs \
        ../_RepaLib/bmp/Codec/BMP/BitmapInfoV3.hs \
-       ../_RepaLib/bmp/Codec/BMP/FileHeader.hs \
        ../_RepaLib/bmp/Codec/BMP/BitmapInfoV4.hs \
        ../_RepaLib/bmp/Codec/BMP/BitmapInfoV5.hs \
-       ../_RepaLib/bmp/Codec/BMP/BitmapInfo.hs \
-       ../_RepaLib/bmp/Codec/BMP/Base.hs \
-       ../_RepaLib/bmp/Codec/BMP/Unpack.hs \
-       ../_RepaLib/bmp/Codec/BMP/Pack.hs \
+       ../_RepaLib/bmp/Codec/BMP/CIEXYZ.hs \
+       ../_RepaLib/bmp/Codec/BMP/Compression.hs \
+       ../_RepaLib/bmp/Codec/BMP/Error.hs \
+       ../_RepaLib/bmp/Codec/BMP/FileHeader.hs \
        ../_RepaLib/bmp/Codec/BMP.hs \
+       ../_RepaLib/bmp/Codec/BMP/Pack.hs \
+       ../_RepaLib/bmp/Codec/BMP/Unpack.hs \
+       ../_RepaLib/quickcheck/Test/QuickCheck/All.hs \
+       ../_RepaLib/quickcheck/Test/QuickCheck/Arbitrary.hs \
        ../_RepaLib/quickcheck/Test/QuickCheck/Exception.hs \
-       ../_RepaLib/quickcheck/Test/QuickCheck/Text.hs \
-       ../_RepaLib/quickcheck/Test/QuickCheck/State.hs \
+       ../_RepaLib/quickcheck/Test/QuickCheck/Function.hs \
        ../_RepaLib/quickcheck/Test/QuickCheck/Gen.hs \
-       ../_RepaLib/quickcheck/Test/QuickCheck/Arbitrary.hs \
+       ../_RepaLib/quickcheck/Test/QuickCheck.hs \
        ../_RepaLib/quickcheck/Test/QuickCheck/Modifiers.hs \
+       ../_RepaLib/quickcheck/Test/QuickCheck/Monadic.hs \
+       ../_RepaLib/quickcheck/Test/QuickCheck/Poly.hs \
        ../_RepaLib/quickcheck/Test/QuickCheck/Property.hs \
+       ../_RepaLib/quickcheck/Test/QuickCheck/State.hs \
        ../_RepaLib/quickcheck/Test/QuickCheck/Test.hs \
-       ../_RepaLib/quickcheck/Test/QuickCheck.hs \
-       ../_RepaLib/repa/Data/Array/Repa/QuickCheck.hs \
-       ../_RepaLib/repa/Data/Array/Repa/Shape.hs \
+       ../_RepaLib/quickcheck/Test/QuickCheck/Text.hs \
+       ../_RepaLib/repa-algorithms/Data/Array/Repa/Algorithms/Complex.hs \
+       ../_RepaLib/repa-algorithms/Data/Array/Repa/Algorithms/Convolve.hs \
+       ../_RepaLib/repa-algorithms/Data/Array/Repa/Algorithms/DFT/Center.hs \
+       ../_RepaLib/repa-algorithms/Data/Array/Repa/Algorithms/DFT.hs \
+       ../_RepaLib/repa-algorithms/Data/Array/Repa/Algorithms/DFT/Roots.hs \
+       ../_RepaLib/repa-algorithms/Data/Array/Repa/Algorithms/FFT.hs \
+       ../_RepaLib/repa-algorithms/Data/Array/Repa/Algorithms/Iterate.hs \
+       ../_RepaLib/repa-algorithms/Data/Array/Repa/Algorithms/Matrix.hs \
+       ../_RepaLib/repa-algorithms/Data/Array/Repa/Algorithms/Randomish.hs \
+       ../_RepaLib/repa-bytestring/Data/Array/Repa/ByteString.hs \
+       ../_RepaLib/repa/Data/Array/Repa/Arbitrary.hs \
+       ../_RepaLib/repa/Data/Array/Repa.hs \
        ../_RepaLib/repa/Data/Array/Repa/Index.hs \
+       ../_RepaLib/repa/Data/Array/Repa/Internals/Base.hs \
+       ../_RepaLib/repa/Data/Array/Repa/Internals/Elt.hs \
+       ../_RepaLib/repa/Data/Array/Repa/Internals/EvalBlockwise.hs \
+       ../_RepaLib/repa/Data/Array/Repa/Internals/EvalChunked.hs \
+       ../_RepaLib/repa/Data/Array/Repa/Internals/EvalCursored.hs \
+       ../_RepaLib/repa/Data/Array/Repa/Internals/EvalReduction.hs \
+       ../_RepaLib/repa/Data/Array/Repa/Internals/Forcing.hs \
+       ../_RepaLib/repa/Data/Array/Repa/Internals/Gang.hs \
+       ../_RepaLib/repa/Data/Array/Repa/Internals/Select.hs \
+       ../_RepaLib/repa/Data/Array/Repa/Operators/IndexSpace.hs \
+       ../_RepaLib/repa/Data/Array/Repa/Operators/Interleave.hs \
+       ../_RepaLib/repa/Data/Array/Repa/Operators/Mapping.hs \
+       ../_RepaLib/repa/Data/Array/Repa/Operators/Modify.hs \
+       ../_RepaLib/repa/Data/Array/Repa/Operators/Reduction.hs \
+       ../_RepaLib/repa/Data/Array/Repa/Operators/Select.hs \
+       ../_RepaLib/repa/Data/Array/Repa/Operators/Traverse.hs \
+       ../_RepaLib/repa/Data/Array/Repa/Properties.hs \
+       ../_RepaLib/repa/Data/Array/Repa/Shape.hs \
        ../_RepaLib/repa/Data/Array/Repa/Slice.hs \
-       ../_RepaLib/repa/Data/Array/Repa.hs \
-       ../_RepaLib/repa-bytestring/Data/Array/Repa/ByteString.hs \
-       ../_RepaLib/repa-algorithms/Data/Array/Repa/Algorithms/Convolve.hs \
-       ../_RepaLib/repa-io/Data/Array/Repa/IO/Timing.hs \
+       ../_RepaLib/repa/Data/Array/Repa/Specialised/Dim2.hs \
+       ../_RepaLib/repa/Data/Array/Repa/Stencil/Base.hs \
+       ../_RepaLib/repa/Data/Array/Repa/Stencil.hs \
+       ../_RepaLib/repa/Data/Array/Repa/Stencil/Template.hs \
+       ../_RepaLib/repa-io/Data/Array/Repa/IO/Binary.hs \
        ../_RepaLib/repa-io/Data/Array/Repa/IO/BMP.hs \
+       ../_RepaLib/repa-io/Data/Array/Repa/IO/ColorRamp.hs \
+       ../_RepaLib/repa-io/Data/Array/Repa/IO/Internals/Text.hs \
+       ../_RepaLib/repa-io/Data/Array/Repa/IO/Matrix.hs \
+       ../_RepaLib/repa-io/Data/Array/Repa/IO/Timing.hs \
+       ../_RepaLib/repa-io/Data/Array/Repa/IO/Vector.hs \
        src/Main.hs
 PROG_ARGS += 12 lena.bmp out.bmp
 HC_OPTS += -threaded -i. -i../_RepaLib/repa -i../_RepaLib/repa-algorithms -i../_RepaLib/repa-io -i../_RepaLib/bmp -i../_RepaLib/repa-bytestring -i../_RepaLib/quickcheck -package base -package binary -package bytestring -package dph-base -package dph-prim-par -package dph-prim-seq -package extensible-exceptions -package ghc -package mtl -package old-time -package random -package vector
 CLEAN_FILES += out.bmp
 include $(TOP)/mk/target.mk
+
index fcca60b..a7ad0ce 100644 (file)
@@ -1,13 +1,15 @@
-{-# LANGUAGE PackageImports, BangPatterns #-}
+{-# LANGUAGE PackageImports, BangPatterns, TemplateHaskell, QuasiQuotes #-}
 {-# OPTIONS -Wall -fno-warn-missing-signatures -fno-warn-incomplete-patterns #-}
 
 import Data.List
 import Control.Monad
 import System.Environment
+import Data.Word
 import Data.Array.Repa.IO.BMP
 import Data.Array.Repa.IO.Timing
-import Data.Array.Repa.Algorithms.Convolve
+import Data.Array.Repa.Algorithms.Iterate
 import Data.Array.Repa                         as A
+import Data.Array.Repa.Stencil         as A
 import Prelude                         as P
 
 main 
@@ -17,53 +19,71 @@ main
         _                              -> usage
 
 usage  = putStr $ unlines
-       [ "repa-blur <iterations> <fileIn.bmp> <fileOut.bmp>" ]
+       [ "repa-blur <iterations::Int> <fileIn.bmp> <fileOut.bmp>" ]
        
--- TODO: component-wise is filthy.
---       do this with a DIM3 array.
 run iterations fileIn fileOut
- = do  (arrRed, arrGreen, arrBlue)
-               <- liftM (either (error . show) id) 
-               $  readComponentsFromBMP fileIn
+ = do  comps   <- liftM (either (error . show) id) 
+               $  readComponentsListFromBMP fileIn
 
-       arrRed `deepSeqArray` arrGreen `deepSeqArray` arrBlue `deepSeqArray` return ()
-
-       ((arrRed', arrGreen', arrBlue'), tElapsed)
-               <- time $ let result    = blurComponents iterations arrRed arrGreen arrBlue
-                         in  result `seq` return result
+       comps `deepSeqArrays` return ()
+       
+       (comps', _)
+        <- time $ let  comps'  = P.map (process iterations) comps
+                  in   comps' `deepSeqArrays` return comps'
        
        -- putStr $ prettyTime tElapsed
-       putStrLn "Done"
+        putStrLn "Done"
                        
-       writeComponentsToBMP fileOut arrRed' arrGreen' arrBlue'
+       writeComponentsListToBMP fileOut comps'
 
+{-# NOINLINE process #-}
+process        :: Int -> Array DIM2 Word8 -> Array DIM2 Word8
+process iterations = demote . blur iterations . promote
 
+       
+{-# NOINLINE promote #-}
+promote        :: Array DIM2 Word8 -> Array DIM2 Double
+promote arr@(Array _ [Region RangeAll (GenManifest _)])
+ = arr `deepSeqArray` force
+ $ A.map ffs arr
 
--- | Blur all the components of an image.
-blurComponents iterations arrRed arrGreen arrBlue
- = let process arr     
-               = force 
-               $ A.map truncate 
-               $ blurs iterations
-               $ force
-               $ A.map fromIntegral 
-               $ arr
+ where {-# INLINE ffs #-}
+       ffs     :: Word8 -> Double
+       ffs x   =  fromIntegral (fromIntegral x :: Int)
 
-       [arrRed', arrGreen', arrBlue']
-               = P.map process [arrRed, arrGreen, arrBlue]
-   in  arrRed' `deepSeqArray` arrGreen' `deepSeqArray` arrBlue' `deepSeqArray`
-         (arrRed', arrGreen', arrBlue')
 
+{-# NOINLINE demote #-}
+demote :: Array DIM2 Double -> Array DIM2 Word8
+demote arr@(Array _ [Region RangeAll (GenManifest _)])
+ = arr `deepSeqArray` force
+ $ A.map ffs arr
 
--- | Run several iterations of blurring.
-blurs  :: Int -> Array DIM2 Double -> Array DIM2 Double
-blurs 0 arr    = arr
-blurs n arr    = blurs (n - 1) (force $ blur arr)
+ where {-# INLINE ffs #-}
+       ffs     :: Double -> Word8
+       ffs x   =  fromIntegral (truncate x :: Int)
 
 
--- | Run a single iteration of blurring.
 {-# NOINLINE blur #-}
+blur   :: Int -> Array DIM2 Double -> Array DIM2 Double
+blur !iterations arr@(Array _ [Region RangeAll (GenManifest _)])
+       = arr `deepSeqArray` force2
+       $ iterateBlockwise' iterations arr
+       $ A.map (/ 159)
+       . mapStencil2 BoundClamp
+         [stencil2|    2  4  5  4  2
+                       4  9 12  9  4
+                       5 12 15 12  5
+                       4  9 12  9  4
+                       2  4  5  4  2 |]
+                       
+
+{- version using convolveOut
+-- | Run several iterations of blurring.
+blur   :: Int -> Array DIM2 Double -> Array DIM2 Double
+blur 0 arr     = arr
+blur n arr     = blurs (n - 1) (force $ blur arr)
+
+-- | Run a single iteration of blurring.
 blur :: Array DIM2 Double -> Array DIM2 Double
 blur input@Manifest{}
  = convolveOut outClamp kernel input
@@ -75,5 +95,5 @@ blur input@Manifest{}
                           5.0, 12.0, 15.0, 12.0, 5.0,
                           4.0,  9.0, 12.0,  9.0, 4.0,
                           2.0,  4.0,  5.0,  4.0, 2.0]
-       
+-}     
        
index 77d67c3..6624e3a 100644 (file)
@@ -2,7 +2,8 @@
 
 -- | Reading and writing uncompressed BMP files.
 --
---   Reading works for both uncompressed 24bit RGB WindowsV3 and 32bit RGBA WindowsV4 formats.
+--   Reading works for both uncompressed 24bit RGB and 32bit RGBA
+--   WindowsV3, WindowsV4 and WindowsV5 formats.
 -- 
 --   Writing is limited to the uncompressed 24bit RGB WindowsV3 format.
 --
 --  >    let (width, height) = bmpDimensions bmp
 --  >    ... 
 --      
+-- Release Notes:
 --
+--  >  * bmp 1.2.0
+--  >    Accept files with zero padding on the end of the file.
+--  >    Accept RGBA files with V3 headers.
+--
+--  >  * bmp 1.1.2   
+--  >    Accept files with the image size field set to zero.
 --
 module Codec.BMP
        ( BMP             (..)
@@ -77,14 +85,12 @@ hGetBMP h
                = BSL.splitAt (fromIntegral sizeOfFileHeader) buf
        
        if (fromIntegral $ BSL.length bufFileHeader) /= sizeOfFileHeader
-        then   return $ Left ErrorReadOfFileHeaderFailed
+        then   return $ Left ErrorFileHeaderTruncated
         else   hGetBMP2 bufRest (decode bufFileHeader)
        
                
 hGetBMP2 buf fileHeader
  -- Check the magic before doing anything else.
- --    If the specified file is not a BMP file then we'd prefer to get 
- --    this error than a `ReadOfImageHeaderFailed`.
  | fileHeaderType fileHeader /= bmpMagic
  = return $ Left $ ErrorBadMagic (fileHeaderType fileHeader)
        
@@ -96,48 +102,59 @@ hGetBMP2 buf fileHeader
        -- split off the image header
        let (bufImageHeader, bufRest)
                = BSL.splitAt (fromIntegral sizeHeader) buf
+        
+        -- How much non-header data is present in the file.
+        -- For uncompressed data without a colour table, the remaining data should
+        -- be the image, but there may also be padding bytes on the end.
+        let physicalBufferSize
+                = (fromIntegral $ BSL.length bufRest) :: Word32
 
        if (fromIntegral $ BSL.length bufImageHeader) /= sizeHeader
-        then   return $ Left ErrorReadOfImageHeaderFailed
-        else   hGetBMP3 fileHeader bufImageHeader bufRest
+        then   return $ Left ErrorImageHeaderTruncated
+        else   hGetBMP3 fileHeader bufImageHeader bufRest physicalBufferSize
 
                        
-hGetBMP3 fileHeader bufImageHeader bufRest
+hGetBMP3 fileHeader bufImageHeader bufRest physicalBufferSize
        | BSL.length bufImageHeader == 40 
        = do    let info        = decode bufImageHeader
-               case checkBitmapInfoV3 info of
+               case checkBitmapInfoV3 info physicalBufferSize of
                 Just err       -> return $ Left err
-                Nothing        -> hGetBMP4 fileHeader (InfoV3 info) bufRest
-                                       (fromIntegral $ dib3ImageSize info)
+                Nothing
+                  | Just imageSize      <- imageSizeFromBitmapInfoV3 info
+                  -> hGetBMP4 fileHeader (InfoV3 info) bufRest imageSize
+                  | otherwise
+                  -> return $ Left $ ErrorInternalErrorPleaseReport
 
        | BSL.length bufImageHeader == 108
        = do    let info        = decode bufImageHeader
-               case checkBitmapInfoV4 info of
+               case checkBitmapInfoV4 info physicalBufferSize of
                 Just err       -> return $ Left err
-                Nothing        -> hGetBMP4 fileHeader (InfoV4 info) bufRest
-                                       (fromIntegral 
-                                               $ dib3ImageSize 
-                                               $ dib4InfoV3 info)
+                Nothing        
+                  | Just imageSize      <- imageSizeFromBitmapInfoV4 info
+                  -> hGetBMP4 fileHeader (InfoV4 info) bufRest imageSize
+                  | otherwise
+                  -> return $ Left $ ErrorInternalErrorPleaseReport
                
        | BSL.length bufImageHeader == 124
        = do    let info        = decode bufImageHeader
-               case checkBitmapInfoV5 info of
+               case checkBitmapInfoV5 info physicalBufferSize of
                 Just err       -> return $ Left err
-                Nothing        -> hGetBMP4 fileHeader (InfoV5 info) bufRest
-                                       (fromIntegral
-                                               $ dib3ImageSize 
-                                               $ dib4InfoV3
-                                               $ dib5InfoV4 info)
+                Nothing        
+                  | Just imageSize      <- imageSizeFromBitmapInfoV5 info
+                  -> hGetBMP4 fileHeader (InfoV5 info) bufRest imageSize
+                  | otherwise
+                  -> return $ Left $ ErrorInternalErrorPleaseReport
                
        | otherwise
-       = return 
-               $ Left 
-               $ ErrorUnhandledBitmapHeaderSize (fromIntegral $ BSL.length bufImageHeader)
+       = return $ Left 
+                $ ErrorUnhandledBitmapHeaderSize 
+                 $ fromIntegral $ BSL.length bufImageHeader
 
 
 hGetBMP4 fileHeader imageHeader bufImage (sizeImage :: Int)
- = if (fromIntegral $ BSL.length bufImage) /= sizeImage
-        then return $ Left ErrorReadOfImageDataFailed
+ = let  bufLen  = fromIntegral $ BSL.length bufImage
+   in   if bufLen < sizeImage
+        then return $ Left $ ErrorImageDataTruncated sizeImage bufLen
         else return 
                $ Right $ BMP 
                { bmpFileHeader         = fileHeader
@@ -152,6 +169,7 @@ writeBMP fileName bmp
  = do  h       <- openBinaryFile fileName WriteMode
        hPutBMP h bmp
        hFlush h
+       hClose h
 
 
 -- | Put a BMP image to a file handle.
index e31c360..7ada6e3 100644 (file)
@@ -1,11 +1,12 @@
 {-# OPTIONS_HADDOCK hide #-}
 module Codec.BMP.Base
-       (BMP    (..))
+       ( BMP   (..))
 where
 import Codec.BMP.FileHeader
 import Codec.BMP.BitmapInfo
 import Data.ByteString
 
+
 -- | A BMP image.
 --     For an uncompressed image, the image data contains triples of BGR component values.
 --     Each line may also have zero pad values on the end, to bring them up to a multiple
@@ -18,6 +19,3 @@ data BMP
        deriving Show
 
 
-
-
-
index a539b35..50ca1c7 100644 (file)
@@ -18,6 +18,7 @@ data BitmapInfo
        | InfoV5 BitmapInfoV5
        deriving (Show)
 
+
 instance Binary BitmapInfo where
  get
   = do size    <- lookAhead getWord32le 
index 9d03d9e..10bc9dd 100644 (file)
@@ -4,13 +4,16 @@ module Codec.BMP.BitmapInfoV3
        ( BitmapInfoV3  (..)
        , Compression (..)
        , sizeOfBitmapInfoV3
-       , checkBitmapInfoV3)
+       , checkBitmapInfoV3
+        , imageSizeFromBitmapInfoV3)
 where
 import Codec.BMP.Error
+import Codec.BMP.Compression
 import Data.Binary
 import Data.Binary.Get 
 import Data.Binary.Put
 
+
 -- | Device Independent Bitmap (DIB) header for Windows V3.
 data BitmapInfoV3
        = BitmapInfoV3                  
@@ -33,6 +36,10 @@ data BitmapInfoV3
        , dib3Compression       :: Compression
 
          -- | (+20) Size of raw image data.
+         --   Some encoders set this to zero, so we need to calculate it based on the
+         --   overall file size.
+         -- 
+         --   If it is non-zero then we check it matches the file size - header size.
        , dib3ImageSize         :: Word32
 
          -- | (+24) Prefered resolution in pixels per meter, along the X axis.
@@ -49,16 +56,6 @@ data BitmapInfoV3
        }
        deriving (Show)
 
-data Compression
-       = CompressionRGB
-       | CompressionRLE8
-       | CompressionRLE4
-       | CompressionBitFields
-       | CompressionJPEG
-       | CompressionPNG
-       | CompressionUnknown Word32
-       deriving (Show, Eq)
-
 
 -- | Size of `BitmapInfoV3` header (in bytes)
 sizeOfBitmapInfoV3 :: Int
@@ -105,53 +102,70 @@ instance Binary BitmapInfoV3 where
        putWord32le     $ dib3ColorsUsed header
        putWord32le     $ dib3ColorsImportant header
        
-       
-instance Binary Compression where
- get
-  = do c       <- getWord32le
-       case c of
-        0      -> return $ CompressionRGB
-        1      -> return $ CompressionRLE8
-        2      -> return $ CompressionRLE4
-        3      -> return $ CompressionBitFields
-        4      -> return $ CompressionJPEG
-        5      -> return $ CompressionPNG
-        _      -> return $ CompressionUnknown c
-       
- put c
-  = case c of
-       CompressionRGB          -> putWord32le 0
-       CompressionRLE8         -> putWord32le 1
-       CompressionRLE4         -> putWord32le 2
-       CompressionBitFields    -> putWord32le 3
-       CompressionJPEG         -> putWord32le 4
-       CompressionPNG          -> putWord32le 5
-       CompressionUnknown x    -> putWord32le x
-       
-       
--- | Check headers for problems and unsupported features.       
---     With a V3 header we only support the uncompressed 24bit RGB format.
-checkBitmapInfoV3 :: BitmapInfoV3 ->  Maybe Error
-checkBitmapInfoV3 header
                
+-- | Check headers for problems and unsupported features.       
+checkBitmapInfoV3 :: BitmapInfoV3 -> Word32 -> Maybe Error
+checkBitmapInfoV3 header physicalBufferSize
+
+        -- We only handle a single color plane.
        | dib3Planes header /= 1
-       = Just  $ ErrorUnhandledPlanesCount 
-               $ fromIntegral $ dib3Planes header
-       
-       | dib3ImageSize header == 0
-       = Just  $ ErrorZeroImageSize
+       = Just  $ ErrorUnhandledPlanesCount $ dib3Planes header
        
-       | dib3ImageSize header `mod` dib3Height header /= 0
-       = Just  $ ErrorLacksWholeNumberOfLines
-
+        -- We only handle 24 and 32 bit images.
        | dib3BitCount header /= 24
-       = Just  $ ErrorUnhandledColorDepth
-               $ fromIntegral $ dib3BitCount header
+        , dib3BitCount header /= 32
+       = Just  $ ErrorUnhandledColorDepth $ dib3BitCount header
+
+        -- If the image size field in the header is non-zero, 
+        -- then it must be the same as the physical size of the image buffer.
+        | headerImageSize               <- dib3ImageSize header
+        , headerImageSize /= 0
+        , fromIntegral headerImageSize /= physicalBufferSize
+        = Just  $ ErrorImagePhysicalSizeMismatch
+                        headerImageSize physicalBufferSize
+
+        -- Check that the physical buffer contains enough image data.
+        -- It may contain more, as some encoders put padding bytes
+        -- on the end.
+        | Just calculatedImageSize      <- imageSizeFromBitmapInfoV3 header
+        , fromIntegral physicalBufferSize < calculatedImageSize
+        = Just  $ ErrorImageDataTruncated 
+                        calculatedImageSize
+                        (fromIntegral physicalBufferSize)
+
+        -- We only handle uncompresssed images.
+        | dib3Compression header /= CompressionRGB
+        = Just  $ ErrorUnhandledCompressionMode (dib3Compression header)
 
-       | dib3Compression header /= CompressionRGB
-       = Just  $ ErrorUnhandledCompressionMode
-       
        | otherwise
        = Nothing
        
 
+-- | Compute the size of the image data from the header.
+--
+--   * We can't just use the 'dib3ImageSize' field because some encoders
+--     set this to zero.
+--
+--   * We also can't use the physical size of the data in the file because
+--     some encoders add zero padding bytes on the end.   
+--
+imageSizeFromBitmapInfoV3 :: BitmapInfoV3 -> Maybe Int
+imageSizeFromBitmapInfoV3 header
+        | dib3BitCount    header == 32
+        , dib3Planes      header == 1
+        , dib3Compression header == CompressionRGB
+        = Just $ fromIntegral (dib3Width header * dib3Height header * 4)
+
+        | dib3BitCount    header == 24
+        , dib3Planes      header == 1
+        , dib3Compression header == CompressionRGB
+        = let   imageBytesPerLine = dib3Width header * 3
+                tailBytesPerLine  = imageBytesPerLine `mod` 4
+                padBytesPerLine   = if tailBytesPerLine > 0
+                                        then 4 - tailBytesPerLine
+                                        else 0
+          in    Just $ fromIntegral (dib3Height header * imageBytesPerLine + padBytesPerLine)
+
+        | otherwise
+        = Nothing
+
index 4a564b6..e6a969f 100644 (file)
@@ -4,7 +4,8 @@ module Codec.BMP.BitmapInfoV4
        ( BitmapInfoV4  (..)
        , CIEXYZ        (..)
        , sizeOfBitmapInfoV4
-       , checkBitmapInfoV4)
+       , checkBitmapInfoV4
+        , imageSizeFromBitmapInfoV4)
 where
 import Codec.BMP.Error
 import Codec.BMP.CIEXYZ
@@ -92,27 +93,44 @@ instance Binary BitmapInfoV4 where
 --     With a V4 header we support both the uncompressed 24bit RGB format,
 --     and the uncompressed 32bit RGBA format.
 --
-checkBitmapInfoV4 :: BitmapInfoV4 ->  Maybe Error
-checkBitmapInfoV4 headerV4
+checkBitmapInfoV4 :: BitmapInfoV4 -> Word32 -> Maybe Error
+checkBitmapInfoV4 headerV4 physicalBufferSize
                
+        -- We only handle a single color plane.
        | dib3Planes headerV3 /= 1
-       = Just  $ ErrorUnhandledPlanesCount 
-               $ fromIntegral $ dib3Planes headerV3
-       
-       | dib3ImageSize headerV3 == 0
-       = Just  $ ErrorZeroImageSize
-       
-       | dib3ImageSize headerV3 `mod` dib3Height headerV3 /= 0
-       = Just  $ ErrorLacksWholeNumberOfLines
+       = Just  $ ErrorUnhandledPlanesCount $ dib3Planes headerV3
+
+        -- We only handle 24 and 32 bit images.
+        | dib3BitCount headerV3 /= 24
+        , dib3BitCount headerV3 /= 32
+        = Just  $ ErrorUnhandledColorDepth $ dib3BitCount headerV3
+
+        -- If the image size field in the header is non-zero, 
+        -- then it must be the same as physical size of the image buffer.
+        | headerImageSize               <- dib3ImageSize headerV3
+        , headerImageSize /= 0
+        , fromIntegral headerImageSize /= physicalBufferSize
+        = Just  $ ErrorImagePhysicalSizeMismatch
+                        headerImageSize physicalBufferSize
+
+        -- Check that the physical buffer contains enough image data.
+        -- It may contain more, as some encoders put padding bytes
+        -- on the end.
+        | Just calculatedImageSize      <- imageSizeFromBitmapInfoV4 headerV4
+        , fromIntegral physicalBufferSize < calculatedImageSize
+        = Just  $ ErrorImageDataTruncated 
+                        calculatedImageSize
+                        (fromIntegral physicalBufferSize)
+
 
        -- Check for valid compression modes ----
 
-       -- uncompressed 24bit RGB
-       | dib3BitCount    headerV3 == 24 
-       , dib3Compression headerV3 == CompressionRGB
-       = Nothing
+        -- uncompressed 32bit RGBA stated as CompressionRGB.
+        | dib3BitCount    headerV3 == 32
+        , dib3Compression headerV3 == CompressionRGB
+        = Nothing
        
-       -- uncompressed 32bit RGBA
+       -- uncompressed 32bit RGBA stated as CompressionBitFields.
        | dib3BitCount    headerV3 == 32
        , dib3Compression headerV3 == CompressionBitFields
        , dib4AlphaMask   headerV4 == 0xff000000
@@ -120,11 +138,53 @@ checkBitmapInfoV4 headerV4
        , dib4GreenMask   headerV4 == 0x0000ff00
        , dib4BlueMask    headerV4 == 0x000000ff
        = Nothing
+
+        -- uncompressed 24bit RGB
+        | dib3BitCount    headerV3 == 24 
+        , dib3Compression headerV3 == CompressionRGB
+        = Nothing
        
        -- Some unsupported compression mode ----
        | otherwise
-       = Just $ ErrorUnhandledCompressionMode
+       = Just $ ErrorUnhandledCompressionMode (dib3Compression headerV3)
        
        where   headerV3 = dib4InfoV3 headerV4
-       
-       
\ No newline at end of file
+
+
+-- | Compute the size of the image data from the header.
+--
+--   * We can't just use the 'dib3ImageSize' field because some encoders
+--     set this to zero.
+--
+--   * We also can't use the physical size of  the data in the file because
+--     some encoders add zero padding bytes on the end.   
+imageSizeFromBitmapInfoV4 :: BitmapInfoV4 -> Maybe Int
+imageSizeFromBitmapInfoV4 headerV4
+        | dib3BitCount    headerV3 == 32
+        , dib3Planes      headerV3 == 1
+        , dib3Compression headerV3 == CompressionRGB
+        = Just $ fromIntegral (dib3Width headerV3 * dib3Height headerV3 * 4)
+
+        | dib3BitCount    headerV3 == 32
+        , dib3Planes      headerV3 == 1
+        , dib3Compression headerV3 == CompressionBitFields
+        , dib4AlphaMask   headerV4 == 0xff000000
+        , dib4RedMask     headerV4 == 0x00ff0000
+        , dib4GreenMask   headerV4 == 0x0000ff00
+        , dib4BlueMask    headerV4 == 0x000000ff
+        = Just $ fromIntegral (dib3Width headerV3 * dib3Height headerV3 * 4)        
+
+        | dib3BitCount    headerV3 == 24
+        , dib3Planes      headerV3 == 1
+        , dib3Compression headerV3 == CompressionRGB
+        = let   imageBytesPerLine = dib3Width headerV3 * 3
+                tailBytesPerLine  = imageBytesPerLine `mod` 4
+                padBytesPerLine   = if tailBytesPerLine > 0
+                                        then 4 - tailBytesPerLine
+                                        else 0
+          in    Just $ fromIntegral (dib3Height headerV3 * imageBytesPerLine + padBytesPerLine)
+
+        | otherwise
+        = Nothing
+
+        where   headerV3 = dib4InfoV3 headerV4
index 35cc7b7..40aca39 100644 (file)
@@ -2,7 +2,8 @@
 module Codec.BMP.BitmapInfoV5
        ( BitmapInfoV5  (..)
        , sizeOfBitmapInfoV5
-       , checkBitmapInfoV5)
+       , checkBitmapInfoV5
+        , imageSizeFromBitmapInfoV5)
 where
 import Codec.BMP.Error
 import Codec.BMP.BitmapInfoV4
@@ -61,8 +62,14 @@ instance Binary BitmapInfoV5 where
        
 -- | Check headers for problems and unsupported features.       
 --     The V5 header doesn't give us any more useful info than the V4 one.
-checkBitmapInfoV5 :: BitmapInfoV5 ->  Maybe Error
-checkBitmapInfoV5 header
-       = checkBitmapInfoV4 $ dib5InfoV4 header
+checkBitmapInfoV5 :: BitmapInfoV5 -> Word32 -> Maybe Error
+checkBitmapInfoV5 header expectedImageSize
+       = checkBitmapInfoV4 (dib5InfoV4 header) expectedImageSize
+
+
+-- | Compute the size of the image data from the header.
+imageSizeFromBitmapInfoV5 :: BitmapInfoV5 -> Maybe Int
+imageSizeFromBitmapInfoV5 
+        = imageSizeFromBitmapInfoV4 . dib5InfoV4
 
 
diff --git a/fibon/Repa/_RepaLib/bmp/Codec/BMP/Compression.hs b/fibon/Repa/_RepaLib/bmp/Codec/BMP/Compression.hs
new file mode 100644 (file)
index 0000000..955745b
--- /dev/null
@@ -0,0 +1,43 @@
+{-# OPTIONS_HADDOCK hide #-}
+module Codec.BMP.Compression
+        (Compression(..))
+where
+import Data.Word
+import Data.Binary
+import Data.Binary.Get  
+import Data.Binary.Put
+
+
+data Compression
+        = CompressionRGB
+        | CompressionRLE8
+        | CompressionRLE4
+        | CompressionBitFields
+        | CompressionJPEG
+        | CompressionPNG
+        | CompressionUnknown Word32
+        deriving (Show, Eq)
+
+
+instance Binary Compression where
+ get
+  = do  c       <- getWord32le
+        case c of
+         0      -> return $ CompressionRGB
+         1      -> return $ CompressionRLE8
+         2      -> return $ CompressionRLE4
+         3      -> return $ CompressionBitFields
+         4      -> return $ CompressionJPEG
+         5      -> return $ CompressionPNG
+         _      -> return $ CompressionUnknown c
+        
+ put c
+  = case c of
+        CompressionRGB          -> putWord32le 0
+        CompressionRLE8         -> putWord32le 1
+        CompressionRLE4         -> putWord32le 2
+        CompressionBitFields    -> putWord32le 3
+        CompressionJPEG         -> putWord32le 4
+        CompressionPNG          -> putWord32le 5
+        CompressionUnknown x    -> putWord32le x
+        
index ab79110..7768b92 100644 (file)
@@ -2,23 +2,59 @@
 module Codec.BMP.Error
        (Error(..))
 where
+import Codec.BMP.Compression
 import Data.Word
 
 -- | Things that can go wrong when loading a BMP file.
 data Error
-       = ErrorReadOfFileHeaderFailed
-       | ErrorReadOfImageHeaderFailed
-       | ErrorReadOfImageDataFailed
-       | ErrorBadMagic                         Word16
+
+        -- | Magic number was not at the start of the file, 
+        --   so this probably isn't a BMP file.
+        = ErrorBadMagic                         Word16
+
+        -- | File is too short to contain a file header.
+       | ErrorFileHeaderTruncated
+
+        -- | File is too short to contain an image header.
+       | ErrorImageHeaderTruncated
+
+        -- | File is too short to contain the image data.
+       | ErrorImageDataTruncated
+        { errorBytesNeeded      :: Int
+        , errorBytesAvailable   :: Int }
+
+        -- | Reserved fields should be zero.
        | ErrorReservedFieldNotZero
-       | ErrorDodgyFileHeaderFieldOffset       Int
-       | ErrorDodgyFileHeaderFieldFileSize     Int
-       | ErrorFileIsTruncated
-       | ErrorUnhandledBitmapHeaderSize        Int
-       | ErrorUnhandledPlanesCount             Int
-       | ErrorUnhandledColorDepth              Int
+
+        -- | The offset to the image data from the file header doesn't
+        --   point anywhere sensible.
+       | ErrorDodgyFileHeaderFieldOffset
+        { errorFileHeaderOffset :: Word32 }
+
+        -- | We handle V3 V4 and V5 image headers, but the size of 
+        --   the header indicates it has some other format.
+       | ErrorUnhandledBitmapHeaderSize
+        { errorBitmapHeaderSize :: Word32 }
+
+        -- | We only handle single color planes.
+       | ErrorUnhandledPlanesCount
+        { errorPlanesCount      :: Word16 }
+
+        -- | We only handle 24 and 32 bit images.
+       | ErrorUnhandledColorDepth
+        { errorColorDepth       :: Word16 }
+
+        -- | We only handle uncompressed images.
        | ErrorUnhandledCompressionMode
-       | ErrorZeroImageSize
-       | ErrorLacksWholeNumberOfLines
+        { errorCompression      :: Compression}
+
+        -- | Mismatch between the image size stated in the header
+        --   and that which is calculuated from the other fields.
+        | ErrorImagePhysicalSizeMismatch 
+        { errorImageSizeFromHeader  :: Word32
+        , errorImageSizeOfBuffer    :: Word32 }
+
+        -- | Something went wrong in the library.
+        | ErrorInternalErrorPleaseReport
        deriving (Eq, Show)
 
index 4bb4d2f..5515cf9 100644 (file)
@@ -70,10 +70,11 @@ checkFileHeader header
        | fileHeaderType header /= bmpMagic
        = Just  $ ErrorBadMagic (fileHeaderType header)
 
-       | fileHeaderFileSize header 
-               < fromIntegral (sizeOfFileHeader + sizeOfBitmapInfoV3)
-       = Just  $ ErrorDodgyFileHeaderFieldFileSize 
-               $ fromIntegral $ fileHeaderFileSize header
+        | fileHeaderFileSize header < fromIntegral sizeOfFileHeader
+        = Just  $ ErrorFileHeaderTruncated
+
+        | fileHeaderFileSize header < fromIntegral (sizeOfFileHeader + sizeOfBitmapInfoV3)
+        = Just  $ ErrorImageHeaderTruncated
 
        | fileHeaderReserved1 header /= 0
        = Just  $ ErrorReservedFieldNotZero
@@ -86,4 +87,4 @@ checkFileHeader header
                $ fromIntegral $ fileHeaderOffset header
 
        | otherwise
-       = Nothing
\ No newline at end of file
+       = Nothing
index a41afff..63861bb 100644 (file)
@@ -65,7 +65,7 @@ packRGBA32ToBMP width height str
                
        errs    = catMaybes             
                        [ checkFileHeader   fileHeader
-                       , checkBitmapInfoV3 bitmapInfoV3 ]
+                       , checkBitmapInfoV3 bitmapInfoV3 (fromIntegral $ BS.length imageData)]
                
    in  case errs of
         [] -> BMP 
index 52f69e7..329b572 100644 (file)
@@ -26,7 +26,7 @@ unpackBMPToRGBA32 bmp
    in  case bitCount of
         24     -> packRGB24ToRGBA32 width height (bmpRawImageData bmp)
         32     -> packRGB32ToRGBA32 width height (bmpRawImageData bmp)
-        _      -> error "Codec.BMP.unpackBMPToRGBA32: unhandled bitcount"
+        _      -> error "Codec.BMP.unpackBMPToRGBA32: unhandled bitcount."
 
 
 -- | Unpack raw, uncompressed 24 bit BMP image data to a string of RGBA component values.
@@ -41,13 +41,16 @@ packRGB24ToRGBA32 width height str
  = let bytesPerLine    = BS.length str `div` height
        padPerLine      = bytesPerLine - width * 3
        sizeDest        = width * height * 4
-   in  if height * (width * 3 + padPerLine) /= BS.length str
-        then error "Codec.BMP.unpackRGB24ToRGBA32: given image dimensions don't match input data."
+
+        -- We allow padding bytes on the end of the image data.
+   in  if BS.length str < height * (width * 3 + padPerLine)
+        then error "Codec.BMP.unpackRGB24ToRGBA32: image data is truncated."
         else unsafePerformIO
                        $ allocaBytes sizeDest      $ \bufDest -> 
                  BS.unsafeUseAsCString str $ \bufSrc  ->
                   do   packRGB24ToRGBA32' width height padPerLine (castPtr bufSrc) (castPtr bufDest)
                        packCStringLen (bufDest, sizeDest)
+
                
 -- We're doing this via Ptrs because we don't want to take the
 -- overhead of doing the bounds checks in ByteString.index.
@@ -89,8 +92,8 @@ packRGB32ToRGBA32
                
 packRGB32ToRGBA32 width height str
   = let sizeDest = height * width * 4
-    in  if sizeDest /= BS.length str
-        then error "Codec.BMP.unpackRGB24ToRGBA32: given image dimensions don't match input data."
+    in  if  BS.length str < sizeDest
+        then error "Codec.BMP.packRGB24ToRGBA32: image data is truncated."
         else unsafePerformIO
                        $ allocaBytes sizeDest      $ \bufDest -> 
                  BS.unsafeUseAsCString str $ \bufSrc  ->
index e2cd54f..250d2d1 100644 (file)
@@ -1,5 +1,5 @@
 Name:                bmp
-Version:             1.1.1.1
+Version:             1.2.0.2
 License:             MIT
 License-file:        LICENSE
 Author:              Ben Lippmeier
@@ -8,7 +8,7 @@ Build-Type:          Simple
 Cabal-Version:       >=1.6
 Stability:           stable
 Category:            Codec
-Homepage:            http://code.haskell.org/~benl/code/bmp-head
+Homepage:            http://code.ouroborus.net/bmp
 Bug-reports:         bmp@ouroborus.net
 Description:
        Read and write uncompressed BMP image files. 100% robust Haskell implementation.
@@ -16,7 +16,7 @@ Description:
 Synopsis:
         Read and write uncompressed BMP image files.
 
-Tested-with: GHC == 6.12.1
+Tested-with: GHC == 7.2
 
 Library
   Build-Depends: 
@@ -25,13 +25,14 @@ Library
         binary               >= 0.5.0.2 && < 0.6.0.0
 
   ghc-options:
-        -Wall -fno-warn-missing-signatures -fno-warn-missing-fields
+        -Wall -fno-warn-missing-signatures
 
   Exposed-modules:
         Codec.BMP
 
   Other-modules:
         Codec.BMP.Base
+        Codec.BMP.Compression
         Codec.BMP.BitmapInfo
         Codec.BMP.BitmapInfoV3
         Codec.BMP.BitmapInfoV4
@@ -41,4 +42,4 @@ Library
         Codec.BMP.FileHeader
         Codec.BMP.Pack
         Codec.BMP.Unpack
-       
\ No newline at end of file
+       
index 5e89f2c..607e610 100644 (file)
@@ -1,5 +1,6 @@
-Copyright (c) 2000-2006, Koen Claessen
-Copyright (c) 2006, Bjorn Bringert
+Copyright (c) 2000-2011, Koen Claessen
+Copyright (c) 2006-2008, Björn Bringert
+Copyright (c) 2009-2011, Nick Smallbone
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without 
index 1bad34a..9ad6a05 100644 (file)
@@ -1,13 +1,15 @@
 Name: QuickCheck
-Version: 2.3.0.2
-Cabal-Version: >= 1.2
+Version: 2.4.2
+Cabal-Version: >= 1.6
 Build-type: Simple
 License: BSD3
 License-file: LICENSE
 Extra-source-files: README
-Copyright: Koen Claessen <koen@chalmers.se>
+Copyright: 2000-2011 Koen Claessen, 2006-2008 Björn Bringert, 2009-2011 Nick Smallbone
 Author: Koen Claessen <koen@chalmers.se>
 Maintainer: QuickCheck developers <quickcheck@projects.haskell.org>
+Bug-reports: mailto:quickcheck@projects.haskell.org
+Tested-with: GHC >=6.8, Hugs
 Homepage: http://code.haskell.org/QuickCheck
 Category:          Testing
 Synopsis:          Automatic testing of Haskell programs
@@ -25,29 +27,46 @@ Description:
        the distribution of test data, and define test
        data generators.
 
-flag splitBase
+source-repository head
+  type:     darcs
+  location: http://code.haskell.org/QuickCheck/devel
+
+source-repository this
+  type:     darcs
+  location: http://code.haskell.org/QuickCheck/stable
+  tag:      2.4.2
+
+flag base3
   Description: Choose the new smaller, split-up base package.
 
-flag extensibleExceptions
-  Description: Choose the even newer, even smaller, split-up base package.
+flag base4
+  Description: Choose the even newer base package with extensible exceptions.
+
+flag templateHaskell
+  Description: Build Test.QuickCheck.All, which uses Template Haskell.
 
 library
-  Build-depends: mtl
-  if flag(extensibleExceptions)
+  -- Choose which library versions to use.
+  if flag(base4)
     Build-depends: base >= 4 && < 5, random
   else
-    if flag(splitBase)
+    if flag(base3)
       Build-depends: base >= 3 && < 4, random
     else
       Build-depends: base < 3
+
+  -- On old versions of GHC use the ghc package to catch ctrl-C.  
   if impl(ghc >= 6.7) && impl(ghc < 6.13)
       Build-depends: ghc
+
+  -- We want to use extensible-exceptions even if linking against base-3.
   if impl(ghc >= 6.9)
     Build-depends: extensible-exceptions
+
+  -- Modules that are always built.
   Exposed-Modules:
     Test.QuickCheck,
     Test.QuickCheck.Arbitrary,
-    Test.QuickCheck.Function,
     Test.QuickCheck.Gen,
     Test.QuickCheck.Monadic,
     Test.QuickCheck.Modifiers,
@@ -56,6 +75,25 @@ library
     Test.QuickCheck.Text,
     Test.QuickCheck.Poly,
     Test.QuickCheck.State
+
+  -- Choose which optional features to build based on which compiler
+  -- we're using. It would be nice to use flags for this but Cabal's
+  -- dependency resolution isn't good enough.
+  if impl(ghc)
+    Exposed-Modules: Test.QuickCheck.Function
+    if flag(templateHaskell) && impl(ghc >= 6.12)
+      Build-depends: template-haskell >= 2.4
+      Exposed-Modules: Test.QuickCheck.All
+    -- GHC < 7.0 doesn't cope with multiple LANGUAGE pragmas in the same
+    -- file, I think...
+    if impl(ghc < 7)
+      Extensions: GeneralizedNewtypeDeriving, MultiParamTypeClasses, Rank2Types, TypeOperators
+  else
+    -- If your Haskell compiler can cope without some of these, please
+    -- send a message to the QuickCheck mailing list!
+    cpp-options: -DNO_TIMEOUT -DNO_NEWTYPE_DERIVING
+    if !impl(hugs)
+      cpp-options: -DNO_ST_MONAD -DNO_MULTI_PARAM_TYPE_CLASSES
   Other-Modules:
     Test.QuickCheck.Exception
   GHC-options:
index c991eb8..aa28bcf 100644 (file)
@@ -1,4 +1,4 @@
-This is QuickCheck 2.1.1, a library for random testing of program properties.
+This is QuickCheck 2.4.1.1, a library for random testing of program properties.
 
 === Installation ===
 
@@ -12,6 +12,10 @@ $ runghc Setup.lhs configure
 $ runghc Setup.lhs build
 $ runghc Setup.lhs install
 
+If you get errors about Template Haskell, try
+
+$ cabal install -f-templateHaskell
+
 === Bugs ===
 
 Please report bugs to the QuickCheck mailing list at
index 0950630..298790f 100644 (file)
@@ -1,3 +1,4 @@
+{-# LANGUAGE CPP #-}
 module Test.QuickCheck
   ( 
     -- * Running tests
@@ -7,6 +8,12 @@ module Test.QuickCheck
   , quickCheckWith
   , quickCheckWithResult
   , quickCheckResult
+    -- ** Running tests verbosely
+  , verboseCheck
+  , verboseCheckWith
+  , verboseCheckWithResult
+  , verboseCheckResult
+  , verbose
     
     -- * Random generation
   , Gen
@@ -31,7 +38,7 @@ module Test.QuickCheck
   , sample
   , sample'
 
-    -- * Arbitrary and CoArbitrary classes.
+    -- * Arbitrary and CoArbitrary classes
   , Arbitrary(..)
   , CoArbitrary(..)
   
@@ -62,7 +69,9 @@ module Test.QuickCheck
   , NonNegative(..)
   , Smart(..)
   , Shrink2(..)
+#ifndef NO_MULTI_PARAM_TYPE_CLASSES
   , Shrinking(..)
+#endif
   , ShrinkState(..)
 
     -- * Properties
@@ -73,9 +82,15 @@ module Test.QuickCheck
   , (==>)
   , forAll
   , forAllShrink
+    -- *** Experimental combinators for conjunction and disjunction
   , (.&.)
+  , (.&&.)
+  , conjoin
+  , (.||.)
+  , disjoin
     -- *** Handling failure
   , whenFail
+  , printTestCase
   , whenFail'
   , expectFailure
   , within
diff --git a/fibon/Repa/_RepaLib/quickcheck/Test/QuickCheck/All.hs b/fibon/Repa/_RepaLib/quickcheck/Test/QuickCheck/All.hs
new file mode 100644 (file)
index 0000000..f2849b7
--- /dev/null
@@ -0,0 +1,121 @@
+{-# LANGUAGE TemplateHaskell, Rank2Types #-}
+-- | Experimental features using Template Haskell.
+-- You need to have a @{-\# LANGUAGE TemplateHaskell \#-}@ pragma in
+-- your module for any of these to work.
+module Test.QuickCheck.All(
+  -- ** Testing all properties in a module.
+  quickCheckAll,
+  verboseCheckAll,
+  forAllProperties,
+  -- ** Testing polymorphic properties.
+  polyQuickCheck,
+  polyVerboseCheck,
+  mono) where
+
+import Language.Haskell.TH
+import Test.QuickCheck.Property hiding (Result)
+import Test.QuickCheck.Test
+import Data.Char
+import Data.List
+import Control.Monad
+
+-- | Test a polymorphic property, defaulting all type variables to 'Integer'.
+--
+-- Invoke as @$('polyQuickCheck' 'prop)@, where @prop@ is a property.
+-- Note that just evaluating @'quickCheck' prop@ in GHCi will seem to
+-- work, but will silently default all type variables to @()@!
+polyQuickCheck :: Name -> ExpQ
+polyQuickCheck x = [| quickCheck $(mono x) |]
+
+-- | Test a polymorphic property, defaulting all type variables to 'Integer'.
+-- This is just a convenience function that combines 'polyQuickCheck' and 'verbose'.
+polyVerboseCheck :: Name -> ExpQ
+polyVerboseCheck x = [| verboseCheck $(mono x) |]
+
+type Error = forall a. String -> a
+
+-- | Monomorphise an arbitrary name by defaulting all type variables to 'Integer'.
+--
+-- For example, if @f@ has type @'Ord' a => [a] -> [a]@
+-- then @$('mono' 'f)@ has type @['Integer'] -> ['Integer']@.
+mono :: Name -> ExpQ
+mono t = do
+  ty0 <- fmap infoType (reify t)
+  let err msg = error $ msg ++ ": " ++ pprint ty0
+  (polys, ctx, ty) <- deconstructType err ty0
+  case polys of
+    [] -> return (VarE t)
+    _ -> do
+      integer <- [t| Integer |]
+      ty' <- monomorphise err integer ty
+      return (SigE (VarE t) ty')
+
+infoType :: Info -> Type
+infoType (ClassOpI _ ty _ _) = ty
+infoType (DataConI _ ty _ _) = ty
+infoType (VarI _ ty _ _) = ty
+
+deconstructType :: Error -> Type -> Q ([Name], Cxt, Type)
+deconstructType err ty0@(ForallT xs ctx ty) = do
+  let plain (PlainTV _) = True
+      plain _ = False
+  unless (all plain xs) $ err "Higher-kinded type variables in type"
+  return (map (\(PlainTV x) -> x) xs, ctx, ty)
+deconstructType _ ty = return ([], [], ty)
+
+monomorphise :: Error -> Type -> Type -> TypeQ
+monomorphise err mono ty@(VarT n) = return mono
+monomorphise err mono (AppT t1 t2) = liftM2 AppT (monomorphise err mono t1) (monomorphise err mono t2)
+monomorphise err mono ty@(ForallT _ _ _) = err $ "Higher-ranked type"
+monomorphise err mono ty = return ty
+
+-- | Test all properties in the current module, using a custom
+-- 'quickCheck' function. The same caveats as with 'quickCheckAll'
+-- apply.
+-- 
+-- @$'forAllProperties'@ has type @('Property' -> 'IO' 'Result') -> 'IO' 'Bool'@.
+-- An example invocation is @$'forAllProperties' 'quickCheckResult'@,
+-- which does the same thing as @$'quickCheckAll'@.
+forAllProperties :: Q Exp -- :: (Property -> IO Result) -> IO Bool
+forAllProperties = do
+  Loc { loc_filename = filename } <- location
+  when (filename == "<interactive>") $ error "don't run this interactively"
+  ls <- runIO (fmap lines (readFile filename))
+  let prefixes = map (takeWhile (\c -> isAlphaNum c || c == '_') . dropWhile (\c -> isSpace c || c == '>')) ls
+      idents = nubBy (\x y -> snd x == snd y) (filter (("prop_" `isPrefixOf`) . snd) (zip [1..] prefixes))
+      quickCheckOne :: (Int, String) -> Q [Exp]
+      quickCheckOne (l, x) = do
+        exists <- return False `recover` (reify (mkName x) >> return True)
+        if exists then sequence [ [| ($(stringE $ x ++ " on " ++ filename ++ ":" ++ show l),
+                                     property $(mono (mkName x))) |] ]
+         else return []
+  [| runQuickCheckAll $(fmap (ListE . concat) (mapM quickCheckOne idents)) |]
+
+-- | Test all properties in the current module.
+-- The name of the property must begin with @prop_@.
+-- Polymorphic properties will be defaulted to 'Integer'.
+-- Returns 'True' if all tests succeeded, 'False' otherwise.
+--
+-- Using 'quickCheckAll' interactively doesn't work.
+-- Instead, add a definition to your module along the lines of
+-- 
+-- > runTests = $quickCheckAll
+-- 
+-- and then execute @runTests@.
+quickCheckAll :: Q Exp
+quickCheckAll = [| $(forAllProperties) quickCheckResult |]
+
+-- | Test all properties in the current module.
+-- This is just a convenience function that combines 'quickCheckAll' and 'verbose'.
+verboseCheckAll :: Q Exp
+verboseCheckAll = [| $(forAllProperties) verboseCheckResult |]
+
+runQuickCheckAll :: [(String, Property)] -> (Property -> IO Result) -> IO Bool
+runQuickCheckAll ps qc =
+  fmap and . forM ps $ \(xs, p) -> do
+    putStrLn $ "=== " ++ xs ++ " ==="
+    r <- qc p
+    return $ case r of
+      Success { } -> True
+      Failure { } -> False
+      NoExpectedFailure { } -> False
index c2965f8..3309ed8 100644 (file)
@@ -1,6 +1,6 @@
 module Test.QuickCheck.Arbitrary
   ( 
-  -- * Arbitrary and CoArbitrary classes.
+  -- * Arbitrary and CoArbitrary classes
     Arbitrary(..)
   , CoArbitrary(..)
   
@@ -127,12 +127,23 @@ instance Arbitrary a => Arbitrary [a] where
   shrink xs = shrinkList shrink xs
 
 shrinkList :: (a -> [a]) -> [a] -> [[a]]
-shrinkList shr xs0 = removeChunks xs0 ++ shrinkOne xs0
+shrinkList shr xs = concat [ removes k n xs | k <- takeWhile (>0) (iterate (`div`2) n) ]
+                 ++ shrinkOne xs
  where
+  n = length xs
+
   shrinkOne []     = []
   shrinkOne (x:xs) = [ x':xs | x'  <- shr x ]
                   ++ [ x:xs' | xs' <- shrinkOne xs ] 
 
+  removes k n xs
+    | k > n     = []
+    | null xs2  = [[]]
+    | otherwise = xs2 : map (xs1 ++) (removes k (n-k) xs2)
+   where
+    xs1 = take k xs
+    xs2 = drop k xs
+
 {-
   -- "standard" definition for lists:
   shrink []     = []
@@ -141,26 +152,6 @@ shrinkList shr xs0 = removeChunks xs0 ++ shrinkOne xs0
                ++ [ x':xs | x'  <- shrink x ]
 -}
 
-removeChunks :: [a] -> [[a]]
-removeChunks xs0 = remC (length xs0) xs0
-   where
-    remC 0 _  = []
-    remC 1 _  = [[]]
-    remC n xs = xs1
-              : xs2
-              : ( [ xs1' ++ xs2 | xs1' <- remC n1 xs1, not (null xs1') ]
-            `ilv` [ xs1 ++ xs2' | xs2' <- remC n2 xs2, not (null xs2') ]
-                )
-     where
-      n1  = n `div` 2
-      xs1 = take n1 xs
-      n2  = n - n1
-      xs2 = drop n1 xs
-
-      []     `ilv` bs     = bs
-      as     `ilv` []     = as
-      (a:as) `ilv` (b:bs) = a : b : (as `ilv` bs)
-
 instance (Integral a, Arbitrary a) => Arbitrary (Ratio a) where
   arbitrary = arbitrarySizedFractional
   shrink    = shrinkRealFrac
@@ -264,12 +255,12 @@ instance Arbitrary Char where
            ++ [' ','\n']
    where
     a <. b  = stamp a < stamp b
-    stamp a = ( not (isLower a)
+    stamp a = ( (not (isLower a)
               , not (isUpper a)
-              , not (isDigit a)
-              , not (a==' ')
+              , not (isDigit a))
+              , (not (a==' ')
               , not (isSpace a)
-              , a
+              , a)
               )
     
 instance Arbitrary Float where
index 548f933..44e2aaf 100644 (file)
@@ -1,3 +1,6 @@
+-- Hide away the nasty implementation-specific ways of catching
+-- exceptions behind a nice API. The main trouble is catching ctrl-C.
+
 {-# LANGUAGE CPP #-}
 module Test.QuickCheck.Exception where
 
@@ -13,7 +16,7 @@ module Test.QuickCheck.Exception where
 #endif
 #endif
 
-#if defined OLD_EXCEPTIONS
+#if defined(OLD_EXCEPTIONS)
 import Control.Exception(evaluate, try, Exception(..))
 #else
 import Control.Exception.Extensible(evaluate, try, SomeException(SomeException)
index 7c730bd..3ca4912 100644 (file)
@@ -1,25 +1,42 @@
 {-# LANGUAGE TypeOperators, GADTs #-}
+-- | Generation of random shrinkable, showable functions.
+-- Not really documented at the moment!
+--
+-- Example of use:
+-- 
+-- >>> :{
+-- >>> let prop :: Fun String Integer -> Bool
+-- >>>     prop (Fun _ f) = f "monkey" == f "banana" || f "banana" == f "elephant"
+-- >>> :}
+-- >>> quickCheck prop
+-- *** Failed! Falsifiable (after 3 tests and 134 shrinks):     
+-- {"elephant"->1, "monkey"->1, _->0}
+--
+-- To generate random values of type @'Fun' a b@,
+-- you must have an instance @'Function' a@.
+-- If your type has a 'Show' instance, you can use 'functionShow' to write the instance; otherwise,
+-- use 'functionMap' to give a bijection between your type and a type that is already an instance of 'Function'.
+-- See the @'Function' [a]@ instance for an example of the latter.
 module Test.QuickCheck.Function
   ( Fun(..)
   , apply
   , (:->)
-  , FunArbitrary(..)
-  , funArbitraryMap
-  , funArbitraryShow
+  , Function(..)
+  , functionMap
+  , functionShow
   )
  where
 
 --------------------------------------------------------------------------
 -- imports
 
-import Test.QuickCheck.Gen
 import Test.QuickCheck.Arbitrary
-import Test.QuickCheck.Property
 import Test.QuickCheck.Poly
-import Test.QuickCheck.Modifiers
 
 import Data.Char
 import Data.Word
+import Data.List( intersperse )
+import Data.Maybe( fromJust )
 
 --------------------------------------------------------------------------
 -- concrete functions
@@ -42,17 +59,17 @@ instance Functor ((:->) a) where
   fmap f (Map g h p) = Map g h (fmap f p)
 
 instance (Show a, Show b) => Show (a:->b) where
-  -- only use this on finite functions
-  show p =
-    "{" ++ (case table p of
-             []        -> ""
-             (_,c):xcs -> concat [ show x ++ "->" ++ show c ++ ","
-                                 | (x,c) <- xcs
-                                 ]
-                       ++ "_->" ++ show c)
-        ++ "}"
-   where
-    xcs = table p
+  show p = showFunction p Nothing
+
+-- only use this on finite functions
+showFunction :: (Show a, Show b) => (a :-> b) -> Maybe b -> String
+showFunction p md =
+  "{" ++ concat (intersperse ", " ( [ show x ++ "->" ++ show c
+                                    | (x,c) <- table p
+                                    ]
+                                 ++ [ "_->" ++ show d
+                                    | Just d <- [md]
+                                    ] )) ++ "}"
 
 -- turning a concrete function into an abstract function (with a default result)
 abstract :: (a :-> c) -> c -> (a -> c)
@@ -73,54 +90,71 @@ table Nil         = []
 table (Table xys) = xys
 table (Map _ h p) = [ (h x, c) | (x,c) <- table p ]
 
---------------------------------------------------------------------------
--- FunArbitrary
+-- finding a default result
 
-class FunArbitrary a where
-  funArbitrary :: Arbitrary c => Gen (a :-> c)
+-- breadth-first search (can't use depth-first search!)
+data Steps a = Step (Steps a) | Fail | Result a
 
-instance (FunArbitrary a, Arbitrary c) => Arbitrary (a :-> c) where
-  arbitrary = funArbitrary
-  shrink    = shrinkFun shrink
+(#) :: Steps a -> Steps a -> Steps a
+Result x # t        = Result x
+s        # Result y = Result y
+Fail     # t        = t
+s        # Fail     = s
+Step s   # Step t   = Step (s # t)
+
+instance Monad Steps where
+  Step s   >>= k = Step (s >>= k)
+  Result x >>= k = k x
+  Fail     >>= k = Fail
+
+  return x = Result x
+
+run :: Steps a -> Maybe a
+run (Step s)   = run s
+run (Result x) = Just x
+run Fail       = Nothing
+
+defaultSteps :: (a :-> c) -> Steps c
+defaultSteps (Pair p)          = defaultSteps p >>= defaultSteps
+defaultSteps (p :+: q)         = Step (defaultSteps p # defaultSteps q)
+defaultSteps (Unit c)          = Result c
+defaultSteps (Table ((_,y):_)) = Result y
+defaultSteps (Map _ _ p)       = defaultSteps p
+defaultSteps _                 = Fail
 
--- basic instances: pairs, sums, units
+defaultResult :: (a :-> c) -> Maybe c
+defaultResult = run . defaultSteps
 
-instance (FunArbitrary a, FunArbitrary b) => FunArbitrary (a,b) where
-  funArbitrary =
-    do p <- funArbitrary
-       return (Pair p)
+--------------------------------------------------------------------------
+-- Function
+
+class Function a where
+  function :: (a->b) -> (a:->b)
+
+-- basic instances
+  
+instance Function () where
+  function f = Unit (f ())
 
-instance (FunArbitrary a, FunArbitrary b) => FunArbitrary (Either a b) where
-  funArbitrary =
-    do p <- funArbitrary
-       q <- funArbitrary
-       return (p :+: q)
+instance Function Word8 where
+  function f = Table [(x,f x) | x <- [0..255]]
 
-instance FunArbitrary () where
-  funArbitrary =
-    do c <- arbitrary
-       return (Unit c)
+instance (Function a, Function b) => Function (a,b) where
+  function f = Pair (function `fmap` function (curry f))
 
-instance FunArbitrary Word8 where
-  funArbitrary =
-    do xys <- sequence [ do y <- arbitrary
-                            return (x,y)
-                       | x <- [0..255]
-                       ]
-       return (Table xys)
+instance (Function a, Function b) => Function (Either a b) where
+  function f = function (f . Left) :+: function (f . Right)
 
--- other instances (using Map)
+-- other instances
 
-funArbitraryMap :: (FunArbitrary a, Arbitrary c) => (b -> a) -> (a -> b) -> Gen (b :-> c)
-funArbitraryMap g h =
-  do p <- funArbitrary
-     return (Map g h p)
+functionMap :: Function b => (a->b) -> (b->a) -> (a->c) -> (a:->c)
+functionMap g h f = Map g h (function (\b -> f (h b)))
 
-funArbitraryShow :: (Show a, Read a, Arbitrary c) => Gen (a :-> c)
-funArbitraryShow = funArbitraryMap show read
+functionShow :: (Show a, Read a) => (a->c) -> (a:->c)
+functionShow f = functionMap show read f
 
-instance FunArbitrary a => FunArbitrary [a] where
-  funArbitrary = funArbitraryMap g h
+instance Function a => Function [a] where
+  function = functionMap g h
    where
     g []     = Left ()
     g (x:xs) = Right (x,xs)
@@ -128,8 +162,8 @@ instance FunArbitrary a => FunArbitrary [a] where
     h (Left _)       = []
     h (Right (x,xs)) = x:xs
 
-instance FunArbitrary a => FunArbitrary (Maybe a) where
-  funArbitrary = funArbitraryMap g h
+instance Function a => Function (Maybe a) where
+  function = functionMap g h
    where
     g Nothing  = Left ()
     g (Just x) = Right x
@@ -137,8 +171,8 @@ instance FunArbitrary a => FunArbitrary (Maybe a) where
     h (Left _)  = Nothing
     h (Right x) = Just x
 
-instance FunArbitrary Bool where
-  funArbitrary = funArbitraryMap g h
+instance Function Bool where
+  function = functionMap g h
    where
     g False = Left ()
     g True  = Right ()
@@ -146,8 +180,8 @@ instance FunArbitrary Bool where
     h (Left _)  = False
     h (Right _) = True
 
-instance FunArbitrary Integer where
-  funArbitrary = funArbitraryMap gInteger hInteger
+instance Function Integer where
+  function = functionMap gInteger hInteger
    where
     gInteger n | n < 0     = Left (gNatural (abs n - 1))
                | otherwise = Right (gNatural n)
@@ -161,34 +195,40 @@ instance FunArbitrary Integer where
     hNatural []     = 0
     hNatural (w:ws) = fromIntegral w + 256 * hNatural ws
 
-instance FunArbitrary Int where
-  funArbitrary = funArbitraryMap fromIntegral fromInteger
+instance Function Int where
+  function = functionMap fromIntegral fromInteger
 
-instance FunArbitrary Char where
-  funArbitrary = funArbitraryMap ord' chr'
+instance Function Char where
+  function = functionMap ord' chr'
    where
     ord' c = fromIntegral (ord c) :: Word8
     chr' n = chr (fromIntegral n)
 
 -- poly instances
 
-instance FunArbitrary A where
-  funArbitrary = funArbitraryMap unA A
+instance Function A where
+  function = functionMap unA A
 
-instance FunArbitrary B where
-  funArbitrary = funArbitraryMap unB B
+instance Function B where
+  function = functionMap unB B
 
-instance FunArbitrary C where
-  funArbitrary = funArbitraryMap unC C
+instance Function C where
+  function = functionMap unC C
 
-instance FunArbitrary OrdA where
-  funArbitrary = funArbitraryMap unOrdA OrdA
+instance Function OrdA where
+  function = functionMap unOrdA OrdA
 
-instance FunArbitrary OrdB where
-  funArbitrary = funArbitraryMap unOrdB OrdB
+instance Function OrdB where
+  function = functionMap unOrdB OrdB
 
-instance FunArbitrary OrdC where
-  funArbitrary = funArbitraryMap unOrdC OrdC
+instance Function OrdC where
+  function = functionMap unOrdC OrdC
+
+-- instance Arbitrary
+
+instance (Function a, CoArbitrary a, Arbitrary b) => Arbitrary (a:->b) where
+  arbitrary = function `fmap` arbitrary
+  shrink    = shrinkFun shrink
 
 --------------------------------------------------------------------------
 -- shrinking
@@ -203,8 +243,8 @@ shrinkFun shr (Pair p) =
 shrinkFun shr (p :+: q) =
   [ p .+. Nil | not (isNil q) ] ++
   [ Nil .+. q | not (isNil p) ] ++
-  [ p' .+. q  | p' <- shrinkFun shr p ] ++
-  [ p  .+. q' | q' <- shrinkFun shr q ]
+  [ p  .+. q' | q' <- shrinkFun shr q ] ++
+  [ p' .+. q  | p' <- shrinkFun shr p ]
  where
   isNil :: (a :-> b) -> Bool
   isNil Nil = True
@@ -237,22 +277,26 @@ shrinkFun shr (Map g h p) =
 --------------------------------------------------------------------------
 -- the Fun modifier
 
-data Fun a b = Fun (a :-> b) (a -> b)
+data Fun a b = Fun (a :-> b, b) (a -> b)
 
-fun :: (a :-> b) -> Fun a b
-fun p = Fun p (abstract p (snd (head (table p))))
+mkFun :: (a :-> b) -> b -> Fun a b
+mkFun p d = Fun (p,d) (abstract p d)
 
 apply :: Fun a b -> (a -> b)
 apply (Fun _ f) = f
 
 instance (Show a, Show b) => Show (Fun a b) where
-  show (Fun p _) = show p
+  show (Fun (p,d) _) = showFunction p (Just d)
 
-instance (FunArbitrary a, Arbitrary b) => Arbitrary (Fun a b) where
-  arbitrary = fun `fmap` arbitrary
+instance (Function a, CoArbitrary a, Arbitrary b) => Arbitrary (Fun a b) where
+  arbitrary =
+    do p <- arbitrary
+       return (mkFun p (fromJust (defaultResult p)))
 
-  shrink (Fun p _) =
-    [ fun p' | p' <- shrink p, _:_ <- [table p'] ]
+  shrink (Fun (p,d) _) =
+       [ mkFun p' d' | p' <- shrink p, Just d' <- [defaultResult p'] ]
+    ++ [ mkFun p' d  | p' <- shrink p ]
+    ++ [ mkFun p d'  | d' <- shrink d ]
 
 --------------------------------------------------------------------------
 -- the end.
index ce4b9b5..518e9c5 100644 (file)
@@ -1,14 +1,15 @@
+-- | Test case generation.
 module Test.QuickCheck.Gen where
 
 --------------------------------------------------------------------------
 -- imports
 
 import System.Random
-  ( RandomGen(..)
-  , Random(..)
+  ( Random
   , StdGen
-  , newStdGen
+  , randomR
   , split
+  , newStdGen
   )
 
 import Control.Monad
@@ -20,14 +21,6 @@ import Control.Applicative
   ( Applicative(..)
   )
 
-import Control.Monad.Reader()
-  -- needed for "instance Monad (a ->)"
-  
-  -- 2005-09-16:
-  -- GHC gives a warning for this. I reported this as a bug. /Koen
-
--- * Test case generation
-
 --------------------------------------------------------------------------
 -- ** Generator type
 
index 1c2fdd7..0f4ca57 100644 (file)
@@ -1,4 +1,40 @@
-{-# LANGUAGE MultiParamTypeClasses, GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE CPP #-}
+#ifndef NO_MULTI_PARAM_TYPE_CLASSES
+{-# LANGUAGE MultiParamTypeClasses #-}
+#endif
+#ifndef NO_NEWTYPE_DERIVING
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+#endif
+-- | Modifiers for test data.
+--
+-- These types do things such as restricting the kind of test data that can be generated.
+-- They can be pattern-matched on in properties as a stylistic
+-- alternative to using explicit quantification.
+-- 
+-- Examples:
+-- 
+-- @
+-- -- Functions cannot be shown (but see "Test.QuickCheck.Function")
+-- prop_TakeDropWhile ('Blind' p) (xs :: ['A']) =
+--   takeWhile p xs ++ dropWhile p xs == xs
+-- @
+--
+-- @
+-- prop_TakeDrop ('NonNegative' n) (xs :: ['A']) =
+--   take n xs ++ drop n xs == xs
+-- @
+--
+-- @
+-- -- cycle does not work for empty lists
+-- prop_Cycle ('NonNegative' n) ('NonEmpty' (xs :: ['A'])) =
+--   take n (cycle xs) == take n (xs ++ cycle xs)
+-- @
+--
+-- @
+-- -- Instead of 'forAll' 'orderedList'
+-- prop_Sort ('Ordered' (xs :: ['OrdA'])) =
+--   sort xs == xs
+-- @
 module Test.QuickCheck.Modifiers
   (
   -- ** Type-level modifiers for changing generator behavior
@@ -27,32 +63,13 @@ import Data.List
   )
 
 --------------------------------------------------------------------------
--- ** arbitrary modifiers
-
--- These datatypes are mainly here to *pattern match* on in properties.
--- This is a stylistic alternative to using explicit quantification.
--- In other words, they should not be replaced by type synonyms, and their
--- constructors should be exported.
-
--- Examples:
-{-
-prop_TakeDropWhile (Blind p) (xs :: [A]) =           -- because functions cannot be shown
-  takeWhile p xs ++ dropWhile p xs == xs
-
-prop_TakeDrop (NonNegative n) (xs :: [A]) =          -- (BTW, also works for negative n)
-  take n xs ++ drop n xs == xs
-
-prop_Cycle (NonNegative n) (NonEmpty (xs :: [A])) =  -- cycle does not work for empty lists
-  take n (cycle xs) == take n (xs ++ cycle xs)
-
-prop_Sort (Ordered (xs :: [OrdA])) =                 -- instead of "forAll orderedList"
-  sort xs == xs
--}
-
---------------------------------------------------------------------------
 -- | @Blind x@: as x, but x does not have to be in the 'Show' class.
 newtype Blind a = Blind a
- deriving ( Eq, Ord, Num, Integral, Real, Enum )
+ deriving ( Eq, Ord
+#ifndef NO_NEWTYPE_DERIVING
+          , Num, Integral, Real, Enum
+#endif
+          )
 
 instance Show (Blind a) where
   show _ = "(*)"
@@ -65,7 +82,11 @@ instance Arbitrary a => Arbitrary (Blind a) where
 --------------------------------------------------------------------------
 -- | @Fixed x@: as x, but will not be shrunk.
 newtype Fixed a = Fixed a
- deriving ( Eq, Ord, Num, Integral, Real, Enum, Show, Read )
+ deriving ( Eq, Ord, Show, Read
+#ifndef NO_NEWTYPE_DERIVING
+          , Num, Integral, Real, Enum
+#endif
+          )
 
 instance Arbitrary a => Arbitrary (Fixed a) where
   arbitrary = Fixed `fmap` arbitrary
@@ -103,8 +124,11 @@ instance Arbitrary a => Arbitrary (NonEmptyList a) where
 --------------------------------------------------------------------------
 -- | @Positive x@: guarantees that @x \> 0@.
 newtype Positive a = Positive a
- deriving ( Eq, Ord, Num, Integral, Real, Enum, Show, Read )
-
+ deriving ( Eq, Ord, Show, Read
+#ifndef NO_NEWTYPE_DERIVING
+          , Num, Integral, Real, Enum
+#endif
+          )
 instance (Num a, Ord a, Arbitrary a) => Arbitrary (Positive a) where
   arbitrary =
     (Positive . abs) `fmap` (arbitrary `suchThat` (/= 0))
@@ -118,7 +142,11 @@ instance (Num a, Ord a, Arbitrary a) => Arbitrary (Positive a) where
 --------------------------------------------------------------------------
 -- | @NonZero x@: guarantees that @x \/= 0@.
 newtype NonZero a = NonZero a
- deriving ( Eq, Ord, Num, Integral, Real, Enum, Show, Read )
+ deriving ( Eq, Ord, Show, Read
+#ifndef NO_NEWTYPE_DERIVING
+          , Num, Integral, Real, Enum
+#endif
+          )
 
 instance (Num a, Ord a, Arbitrary a) => Arbitrary (NonZero a) where
   arbitrary = fmap NonZero $ arbitrary `suchThat` (/= 0)
@@ -128,14 +156,18 @@ instance (Num a, Ord a, Arbitrary a) => Arbitrary (NonZero a) where
 --------------------------------------------------------------------------
 -- | @NonNegative x@: guarantees that @x \>= 0@.
 newtype NonNegative a = NonNegative a
- deriving ( Eq, Ord, Num, Integral, Real, Enum, Show, Read )
+ deriving ( Eq, Ord, Show, Read
+#ifndef NO_NEWTYPE_DERIVING
+          , Num, Integral, Real, Enum
+#endif
+          )
 
 instance (Num a, Ord a, Arbitrary a) => Arbitrary (NonNegative a) where
   arbitrary =
     frequency
       -- why is this distrbution like this?
       [ (5, (NonNegative . abs) `fmap` arbitrary)
-      , (1, return 0)
+      , (1, return (NonNegative 0))
       ]
 
   shrink (NonNegative x) =
@@ -147,7 +179,11 @@ instance (Num a, Ord a, Arbitrary a) => Arbitrary (NonNegative a) where
 --------------------------------------------------------------------------
 -- | @Shrink2 x@: allows 2 shrinking steps at the same time when shrinking x
 newtype Shrink2 a = Shrink2 a
- deriving ( Eq, Ord, Num, Integral, Real, Enum, Show, Read )
+ deriving ( Eq, Ord, Show, Read
+#ifndef NO_NEWTYPE_DERIVING
+          , Num, Integral, Real, Enum
+#endif
+          )
 
 instance Arbitrary a => Arbitrary (Shrink2 a) where
   arbitrary =
@@ -206,6 +242,7 @@ instance Arbitrary a => Arbitrary (Smart a) where
     -- == take k ys ++ drop k ys
     -- == ys
 
+#ifndef NO_MULTI_PARAM_TYPE_CLASSES
 --------------------------------------------------------------------------
 -- | @Shrinking _ x@: allows for maintaining a state during shrinking.
 data Shrinking s a =
@@ -228,5 +265,7 @@ instance (Arbitrary a, ShrinkState s a) => Arbitrary (Shrinking s a) where
     | (x',s') <- shrinkState x s
     ]
 
+#endif /* NO_MULTI_PARAM_TYPE_CLASSES */
+
 --------------------------------------------------------------------------
 -- the end.
index e862a1c..098c6b9 100644 (file)
@@ -1,5 +1,10 @@
+{-# LANGUAGE CPP #-}
+#ifndef NO_ST_MONAD
 {-# LANGUAGE Rank2Types #-}
+#endif
 -- | Allows testing of monadic values.
+-- See the paper \"Testing Monadic Code with QuickCheck\":
+-- <http://www.cse.chalmers.se/~rjmh/Papers/QuickCheckST.ps>.
 module Test.QuickCheck.Monadic where
 
 --------------------------------------------------------------------------
@@ -73,13 +78,15 @@ monadic' :: Monad m => PropertyM m a -> Gen (m Property)
 monadic' (MkPropertyM m) = m (const (return (return (property True))))
 
 monadicIO :: PropertyM IO a -> Property
-monadicIO = monadic property
+monadicIO = monadic morallyDubiousIOProperty
 
+#ifndef NO_ST_MONAD
 monadicST :: (forall s. PropertyM (ST s) a) -> Property
 monadicST m = property (runSTGen (monadic' m))
 
 runSTGen :: (forall s. Gen (ST s a)) -> Gen a
 runSTGen g = MkGen $ \r n -> runST (unGen g r n)
+#endif
 
 --------------------------------------------------------------------------
 -- the end.
index 41a390c..a66a75f 100644 (file)
@@ -1,3 +1,15 @@
+{-# LANGUAGE CPP #-}
+#ifndef NO_NEWTYPE_DERIVING
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+#endif
+-- | Types to help with testing polymorphic properties.
+--
+-- Types 'A', 'B' and 'C' are @newtype@ wrappers around 'Integer' that
+-- implement 'Eq', 'Show', 'Arbitrary' and 'CoArbitrary'. Types
+-- 'OrdA', 'OrdB' and 'OrdC' also implement 'Ord' and 'Num'.
+--
+-- See also "Test.QuickCheck.All" for an experimental way of testing
+-- polymorphic properties.
 module Test.QuickCheck.Poly
   ( A(..), B(..), C(..)
   , OrdA(..), OrdB(..), OrdC(..)
@@ -63,7 +75,11 @@ instance CoArbitrary C where
 -- OrdA
 
 newtype OrdA = OrdA{ unOrdA :: Integer }
-  deriving ( Eq, Ord )
+  deriving ( Eq, Ord
+#ifndef NO_NEWTYPE_DERIVING
+           , Num
+#endif
+           )
 
 instance Show OrdA where
   showsPrec n (OrdA x) = showsPrec n x
@@ -78,7 +94,11 @@ instance CoArbitrary OrdA where
 -- OrdB
 
 newtype OrdB = OrdB{ unOrdB :: Integer }
-  deriving ( Eq, Ord )
+  deriving ( Eq, Ord
+#ifndef NO_NEWTYPE_DERIVING
+           , Num
+#endif
+           )
 
 instance Show OrdB where
   showsPrec n (OrdB x) = showsPrec n x
@@ -93,7 +113,11 @@ instance CoArbitrary OrdB where
 -- OrdC
 
 newtype OrdC = OrdC{ unOrdC :: Integer }
-  deriving ( Eq, Ord )
+  deriving ( Eq, Ord
+#ifndef NO_NEWTYPE_DERIVING
+           , Num
+#endif
+           )
 
 instance Show OrdC where
   showsPrec n (OrdC x) = showsPrec n x
index c0a381d..0ef4da7 100644 (file)
@@ -1,4 +1,4 @@
-{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE CPP #-}
 module Test.QuickCheck.Property where
 
 --------------------------------------------------------------------------
@@ -10,26 +10,47 @@ import Test.QuickCheck.Text( showErr, putLine )
 import Test.QuickCheck.Exception
 import Test.QuickCheck.State
 
-import Control.Concurrent
-  ( forkIO
-  , threadDelay
-  , killThread
-  , newEmptyMVar
-  , takeMVar
-  , putMVar
-  )
-
-import System.IO
-  ( hFlush
-  , stdout
-  )
+#ifndef NO_TIMEOUT
+import System.Timeout(timeout)
+#endif
+import Data.Maybe
 
 --------------------------------------------------------------------------
--- fixeties
+-- fixities
 
 infixr 0 ==>
 infixr 1 .&.
--- infixr 1 .&&.
+infixr 1 .&&.
+infixr 1 .||.
+
+-- The story for exception handling:
+--
+-- To avoid insanity, we have rules about which terms can throw
+-- exceptions when we evaluate them:
+--   * A rose tree must evaluate to WHNF without throwing an exception
+--   * The 'ok' component of a Result must evaluate to Just True or
+--     Just False or Nothing rather than raise an exception
+--   * IORose _ must never throw an exception when executed
+--
+-- Both rose trees and Results may loop when we evaluate them, though,
+-- so we have to be careful not to force them unnecessarily.
+--
+-- We also have to be careful when we use fmap or >>= in the Rose
+-- monad that the function we supply is total, or else use
+-- protectResults afterwards to install exception handlers. The
+-- mapResult function on Properties installs an exception handler for
+-- us, though.
+--
+-- Of course, the user is free to write "error "ha ha" :: Result" if
+-- they feel like it. We have to make sure that any user-supplied Rose
+-- Results or Results get wrapped in exception handlers, which we do by:
+--   * Making the 'property' function install an exception handler
+--     round its argument. This function always gets called in the
+--     right places, because all our Property-accepting functions are
+--     actually polymorphic over the Testable class so they have to
+--     call 'property'.
+--   * Installing an exception handler round a Result before we put it
+--     in a rose tree (the only place Results can end up).
 
 --------------------------------------------------------------------------
 -- * Property and Testable types
@@ -47,69 +68,89 @@ instance Testable Bool where
   property = property . liftBool
 
 instance Testable Result where
-  property = return . MkProp . return . return
+  property = return . MkProp . protectResults . return
 
 instance Testable Prop where
-  property = return . protectProp
+  property (MkProp r) = return . MkProp . ioRose . return $ r
 
 instance Testable prop => Testable (Gen prop) where
   property mp = do p <- mp; property p
 
-instance Testable prop => Testable (IO prop) where
-  property = fmap (MkProp . IORose . fmap unProp) . promote . fmap property
+-- | Do I/O inside a property. This can obviously lead to unrepeatable
+-- testcases, so use with care.
+morallyDubiousIOProperty :: Testable prop => IO prop -> Property
+morallyDubiousIOProperty = fmap (MkProp . ioRose . fmap unProp) . promote . fmap property
 
 instance (Arbitrary a, Show a, Testable prop) => Testable (a -> prop) where
   property f = forAllShrink arbitrary shrink f
 
+-- ** Exception handling
+protect :: (AnException -> a) -> IO a -> IO a
+protect f x = either f id `fmap` tryEvaluateIO x
+
 --------------------------------------------------------------------------
 -- ** Type Prop
 
--- is this the right level to be abstract at?
-
-newtype Prop = MkProp{ unProp :: Rose (IO Result) }
-
-protectProp :: Prop -> Prop
-protectProp (MkProp r) =
-  MkProp . IORose $ do
-    (x, rs) <- unpackRose r
-    return (MkRose x rs)
+newtype Prop = MkProp{ unProp :: Rose Result }
 
 -- ** type Rose
 
--- We never allow a rose tree to be _|_. This makes avoiding
--- exceptions easier.
--- This relies on the fact that the 'property' function never returns _|_.
 data Rose a = MkRose a [Rose a] | IORose (IO (Rose a))
-
-join :: Rose (Rose a) -> Rose a
-join (IORose rs) = IORose (fmap join rs)
-join (MkRose (IORose rm) rs) = IORose $ do r <- rm; return (join (MkRose r rs))
-join (MkRose (MkRose x ts) tts) =
+-- Only use IORose if you know that the argument is not going to throw an exception!
+-- Otherwise, try ioRose.
+ioRose :: IO (Rose Result) -> Rose Result
+ioRose = IORose . protectRose
+
+joinRose :: Rose (Rose a) -> Rose a
+joinRose (IORose rs) = IORose (fmap joinRose rs)
+joinRose (MkRose (IORose rm) rs) = IORose $ do r <- rm; return (joinRose (MkRose r rs))
+joinRose (MkRose (MkRose x ts) tts) =
   -- first shrinks outer quantification; makes most sense
-  MkRose x (map join tts ++ ts)
+  MkRose x (map joinRose tts ++ ts)
   -- first shrinks inner quantification
-  --MkRose x (ts ++ map join tts)
+  --MkRose x (ts ++ map joinRose tts)
 
 instance Functor Rose where
-  fmap f (IORose rs) = IORose (fmap (fmap f) rs)
+  -- f must be total
+  fmap f (IORose rs)   = IORose (fmap (fmap f) rs)
   fmap f (MkRose x rs) = MkRose (f x) [ fmap f r | r <- rs ]
 
 instance Monad Rose where
   return x = MkRose x []
-  m >>= k  = join (fmap k m)
-
-unpackRose :: Rose (IO Result) -> IO (IO Result, [Rose (IO Result)])
-unpackRose rose = either (\e -> (return (exception "Exception" e), [])) id
-                  `fmap` tryEvaluateIO (unpack rose)
-  where unpack (MkRose x xs) = return (x, xs)
-        unpack (IORose m) = m >>= unpack
+  -- k must be total
+  m >>= k  = joinRose (fmap k m)
+
+-- Execute the "IORose" bits of a rose tree, returning a tree
+-- constructed by MkRose.
+reduceRose :: Rose Result -> IO (Rose Result)
+reduceRose r@(MkRose _ _) = return r
+reduceRose (IORose m) = m >>= reduceRose
+
+-- Apply a function to the outermost MkRose constructor of a rose tree.
+-- The function must be total!
+onRose :: (a -> [Rose a] -> Rose a) -> Rose a -> Rose a
+onRose f (MkRose x rs) = f x rs
+onRose f (IORose m) = IORose (fmap (onRose f) m)
+
+-- Wrap a rose tree in an exception handler.
+protectRose :: IO (Rose Result) -> IO (Rose Result)
+protectRose = protect (return . exception "Exception")
+
+-- Wrap all the Results in a rose tree in exception handlers.
+protectResults :: Rose Result -> Rose Result
+protectResults = onRose $ \x rs ->
+  IORose $ do
+    y <- protectResult (return x)
+    return (MkRose y (map protectResults rs))
 
 -- ** Result type
 
 -- | Different kinds of callbacks
 data Callback
-  = PostTest (State -> Result -> IO ())         -- ^ Called just after a test
-  | PostFinalFailure (State -> Result -> IO ()) -- ^ Called with the final failing test-case
+  = PostTest CallbackKind (State -> Result -> IO ())         -- ^ Called just after a test
+  | PostFinalFailure CallbackKind (State -> Result -> IO ()) -- ^ Called with the final failing test-case
+data CallbackKind = Counterexample    -- ^ Affected by the 'verbose' combinator
+                  | NotCounterexample -- ^ Not affected by the 'verbose' combinator
 
 -- | The result of a single test.
 data Result
@@ -138,8 +179,7 @@ exception msg err = failed{ reason = msg ++ ": '" ++ showErr err ++ "'",
                             interrupted = isInterrupt err }
 
 protectResult :: IO Result -> IO Result
-protectResult m = either (exception "Exception") id `fmap` tryEvaluateIO (fmap force m)
-  where force res = ok res == Just False `seq` res
+protectResult = protect (exception "Exception")
 
 succeeded :: Result 
 succeeded = result{ ok = Just True }
@@ -153,28 +193,19 @@ rejected = result{ ok = Nothing }
 --------------------------------------------------------------------------
 -- ** Lifting and mapping functions
 
-liftBool :: Bool -> Property
-liftBool b = liftResult $
-  result
-  { ok     = Just b
-  , reason = if b then "" else "Falsifiable"
-  }
-
-liftResult :: Result -> Property
-liftResult r = liftIOResult (return r)
-
-liftIOResult :: IO Result -> Property
-liftIOResult m = property (MkProp (return m))
+liftBool :: Bool -> Result
+liftBool True = succeeded
+liftBool False = failed { reason = "Falsifiable" }
 
 mapResult :: Testable prop => (Result -> Result) -> prop -> Property
-mapResult f = mapIOResult (fmap f)
+mapResult f = mapRoseResult (protectResults . fmap f)
 
-mapIOResult :: Testable prop => (IO Result -> IO Result) -> prop -> Property
-mapIOResult f = mapRoseIOResult (fmap (f . protectResult))
+mapTotalResult :: Testable prop => (Result -> Result) -> prop -> Property
+mapTotalResult f = mapRoseResult (fmap f)
 
--- f here has to be total.
-mapRoseIOResult :: Testable prop => (Rose (IO Result) -> Rose (IO Result)) -> prop -> Property
-mapRoseIOResult f = mapProp (\(MkProp t) -> MkProp (f t))
+-- f here mustn't throw an exception (rose tree invariant).
+mapRoseResult :: Testable prop => (Rose Result -> Rose Result) -> prop -> Property
+mapRoseResult f = mapProp (\(MkProp t) -> MkProp (f t))
 
 mapProp :: Testable prop => (Prop -> Prop) -> prop -> Property
 mapProp f = fmap f . property
@@ -193,31 +224,29 @@ shrinking :: Testable prop =>
              (a -> [a])  -- ^ 'shrink'-like function.
           -> a           -- ^ The original argument
           -> (a -> prop) -> Property
-shrinking shrinker x0 pf = fmap (MkProp . join . fmap unProp) (promote (props x0))
+shrinking shrinker x0 pf = fmap (MkProp . joinRose . fmap unProp) (promote (props x0))
  where
   props x =
     MkRose (property (pf x)) [ props x' | x' <- shrinker x ]
 
 -- | Disables shrinking for a property altogether.
 noShrinking :: Testable prop => prop -> Property
-noShrinking = mapRoseIOResult f
-  where f (MkRose mres _ts) = MkRose mres []
-        f (IORose rm) = IORose (fmap f rm)
+noShrinking = mapRoseResult (onRose (\res _ -> MkRose res []))
 
 -- | Adds a callback
 callback :: Testable prop => Callback -> prop -> Property
-callback cb = mapResult (\res -> res{ callbacks = cb : callbacks res })
+callback cb = mapTotalResult (\res -> res{ callbacks = cb : callbacks res })
 
--- | Prints a message to the terminal after the last failure of a property.
-whenFailPrint :: Testable prop => String -> prop -> Property
-whenFailPrint s =
-  callback $ PostFinalFailure $ \st _res ->
+-- | Prints a message to the terminal as part of the counterexample.
+printTestCase :: Testable prop => String -> prop -> Property
+printTestCase s =
+  callback $ PostFinalFailure Counterexample $ \st _res ->
     putLine (terminal st) s
 
 -- | Performs an 'IO' action after the last failure of a property.
 whenFail :: Testable prop => IO () -> prop -> Property
 whenFail m =
-  callback $ PostFinalFailure $ \_st _res ->
+  callback $ PostFinalFailure NotCounterexample $ \_st _res ->
     m
 
 -- | Performs an 'IO' action every time a property fails. Thus,
@@ -225,14 +254,26 @@ whenFail m =
 -- failures along the way.
 whenFail' :: Testable prop => IO () -> prop -> Property
 whenFail' m =
-  callback $ PostTest $ \_st res ->
+  callback $ PostTest NotCounterexample $ \_st res ->
     if ok res == Just False
       then m
       else return ()
 
+-- | Prints out the generated testcase every time the property is tested,
+-- like 'verboseCheck' from QuickCheck 1.
+-- Only variables quantified over /inside/ the 'verbose' are printed.
+verbose :: Testable prop => prop -> Property
+verbose = mapResult (\res -> res { callbacks = newCallbacks (callbacks res) ++ callbacks res })
+  where newCallbacks cbs =
+          PostTest Counterexample (\st res -> putLine (terminal st) (status res ++ ":")):
+          [ PostTest Counterexample f | PostFinalFailure Counterexample f <- cbs ]
+        status MkResult{ok = Just True} = "Passed"
+        status MkResult{ok = Just False} = "Failed"
+        status MkResult{ok = Nothing} = "Skipped (precondition false)"
+
 -- | Modifies a property so that it is expected to fail for some test cases.
 expectFailure :: Testable prop => prop -> Property
-expectFailure = mapResult (\res -> res{ expect = False })
+expectFailure = mapTotalResult (\res -> res{ expect = False })
 
 -- | Attaches a label to a property. This is used for reporting
 -- test case distribution.
@@ -259,10 +300,10 @@ cover :: Testable prop =>
       -> Int    -- ^ The required percentage (0-100) of test cases.
       -> String -- ^ Label for the test case class.
       -> prop -> Property
-cover b n s = mapResult $ \res ->
-        case b of
-         True  -> res{ stamp  = (s,n) : stamp res }
-         False -> res
+cover True n s = n `seq` s `listSeq` (mapTotalResult $ \res -> res { stamp = (s,n) : stamp res })
+  where [] `listSeq` z = z
+        (x:xs) `listSeq` z = x `seq` xs `listSeq` z
+cover False _ _ = property
 
 -- | Implication for properties: The resulting property holds if
 -- the first argument is 'False', or if the given property holds.
@@ -273,39 +314,20 @@ True  ==> p = property p
 -- | Considers a property failed if it does not complete within
 -- the given number of microseconds.
 within :: Testable prop => Int -> prop -> Property
-within n = mapIOResult race
- where
-  race ior =
-    do put "Race starts ..."
-       resV <- newEmptyMVar
-       
-       let waitAndFail =
-             do put "Waiting ..."
-                threadDelay n
-                put "Done waiting!"
-                putMVar resV (failed {reason = "Time out"})
-           
-           evalProp =
-             do put "Evaluating Result ..."
-                res <- protectResult ior
-                put "Evaluating OK ..."
-                putMVar resV res
-       
-       pid1  <- forkIO evalProp
-       pid2  <- forkIO waitAndFail
-
-       put "Blocking ..."
-       res <- takeMVar resV
-       put "Killing threads ..."
-       killThread pid1
-       killThread pid2
-       put ("Got Result: " ++ show (ok res))
-       return res
-         
-
-  put s | True      = do return ()
-        | otherwise = do putStrLn s
-                         hFlush stdout
+within n = mapRoseResult f
+  -- We rely on the fact that the property will catch the timeout
+  -- exception and turn it into a failed test case.
+  where
+    f rose = ioRose $ do
+      let m `orError` x = fmap (fromMaybe (error x)) m
+      MkRose res roses <- timeout n (reduceRose rose) `orError`
+                          "within: timeout exception not caught in Rose Result"
+      res' <- timeout n (protectResult (return res)) `orError`
+              "within: timeout exception not caught in Result"
+      return (MkRose res' (map f roses))
+#ifdef NO_TIMEOUT
+    timeout _ = fmap Just
+#endif
 
 -- | Explicit universal quantification: uses an explicitly given
 -- test case generator.
@@ -313,8 +335,7 @@ forAll :: (Show a, Testable prop)
        => Gen a -> (a -> prop) -> Property
 forAll gen pf =
   gen >>= \x ->
-    whenFailPrint (show x) $
-      property (pf x)
+    printTestCase (show x) (pf x)
 
 -- | Like 'forAll', but tries to shrink the argument for failing test cases.
 forAllShrink :: (Show a, Testable prop)
@@ -322,21 +343,89 @@ forAllShrink :: (Show a, Testable prop)
 forAllShrink gen shrinker pf =
   gen >>= \x ->
     shrinking shrinker x $ \x' ->
-      whenFailPrint (show x') $
-        property (pf x')
+      printTestCase (show x') (pf x')
 
+-- | Nondeterministic choice: 'p1' '.&.' 'p2' picks randomly one of
+-- 'p1' and 'p2' to test. If you test the property 100 times it
+-- makes 100 random choices.
 (.&.) :: (Testable prop1, Testable prop2) => prop1 -> prop2 -> Property
 p1 .&. p2 =
   arbitrary >>= \b ->
-    whenFailPrint (if b then "LHS" else "RHS") $
+    printTestCase (if b then "LHS" else "RHS") $
       if b then property p1 else property p2
 
-{-
--- TODO
-
+-- | Conjunction: 'p1' '.&&.' 'p2' passes if both 'p1' and 'p2' pass.
 (.&&.) :: (Testable prop1, Testable prop2) => prop1 -> prop2 -> Property
-p1 .&&. p2 = error "not implemented yet"
--}
+p1 .&&. p2 = conjoin [property p1, property p2]
+
+-- | Take the conjunction of several properties.
+conjoin :: Testable prop => [prop] -> Property
+conjoin ps = 
+  do roses <- mapM (fmap unProp . property) ps
+     return (MkProp (conj [] roses))
+ where
+  conj cbs [] =
+    MkRose succeeded{callbacks = cbs} []
+
+  conj cbs (p : ps) = IORose $ do
+    rose@(MkRose result _) <- reduceRose p
+    case ok result of
+      _ | not (expect result) ->
+        return (return failed { reason = "expectFailure may not occur inside a conjunction" })
+      Just True -> return (conj (cbs ++ callbacks result) ps)
+      Just False -> return rose
+      Nothing -> do
+        rose2@(MkRose result2 _) <- reduceRose (conj (cbs ++ callbacks result) ps)
+        return $
+          -- Nasty work to make sure we use the right callbacks
+          case ok result2 of
+            Just True -> MkRose (result2 { ok = Nothing }) []
+            Just False -> rose2
+            Nothing -> rose2
+
+-- | Disjunction: 'p1' '.||.' 'p2' passes unless 'p1' and 'p2' simultaneously fail.
+(.||.) :: (Testable prop1, Testable prop2) => prop1 -> prop2 -> Property
+p1 .||. p2 = disjoin [property p1, property p2]
+
+-- | Take the disjunction of several properties.
+disjoin :: Testable prop => [prop] -> Property
+disjoin ps = 
+  do roses <- mapM (fmap unProp . property) ps
+     return (MkProp (foldr disj (MkRose failed []) roses))
+ where
+  disj :: Rose Result -> Rose Result -> Rose Result
+  disj p q =
+    do result1 <- p
+       case ok result1 of
+         _ | not (expect result1) -> return expectFailureError
+         Just True -> return result1
+         Just False -> do
+           result2 <- q
+           return $
+             if expect result2 then
+               case ok result2 of
+                 Just True -> result2
+                 Just False -> result1 >>> result2
+                 Nothing -> result2
+             else expectFailureError
+         Nothing -> do
+           result2 <- q
+           return (case ok result2 of
+                     _ | not (expect result2) -> expectFailureError
+                     Just True -> result2
+                     _ -> result1)
+
+  expectFailureError = failed { reason = "expectFailure may not occur inside a disjunction" }
+  result1 >>> result2 | not (expect result1 && expect result2) = expectFailureError
+  result1 >>> result2 =
+    result2
+    { reason      = if null (reason result2) then reason result1 else reason result2
+    , interrupted = interrupted result1 || interrupted result2
+    , stamp       = stamp result1 ++ stamp result2
+    , callbacks   = callbacks result1 ++
+                    [PostFinalFailure Counterexample $ \st _res -> putLine (terminal st) ""] ++
+                    callbacks result2
+    }
 
 --------------------------------------------------------------------------
 -- the end.
index 237ecc9..4ad9779 100644 (file)
@@ -7,7 +7,7 @@ import System.Random( StdGen )
 -- State
 
 -- | State represents QuickCheck's internal state while testing a property.
--- The state is made visible to callback functions.
+-- The state is made visible to callback functions.
 data State
   = MkState
   -- static
index d580c8a..d43d322 100644 (file)
@@ -9,13 +9,11 @@ import qualified Test.QuickCheck.Property as P
 import Test.QuickCheck.Text
 import Test.QuickCheck.State
 import Test.QuickCheck.Exception
-import Data.IORef
 
 import System.Random
-  ( RandomGen(..)
+  ( split
   , newStdGen
   , StdGen
-  , split
   )
 
 import Data.Char
@@ -112,7 +110,7 @@ quickCheckWithResult a p =
                  , maxDiscardedTests = maxDiscard a
                  , computeSize       = case replay a of
                                          Nothing    -> computeSize'
-                                         Just (_,s) -> \_ _ -> s
+                                         Just (_,s) -> computeSize' `at0` s
                  , numSuccessTests   = 0
                  , numDiscardedTests = 0
                  , collected         = []
@@ -130,6 +128,28 @@ quickCheckWithResult a p =
           | otherwise =
             (n `mod` maxSize a) * maxSize a `div` (maxSuccess a `mod` maxSize a) + d `div` 10
         n `roundTo` m = (n `div` m) * m
+        at0 f s 0 0 = s
+        at0 f s n d = f n d
+
+-- | Tests a property and prints the results and all test cases generated to 'stdout'.
+-- This is just a convenience function that means the same as 'quickCheck' '.' 'verbose'.
+verboseCheck :: Testable prop => prop -> IO ()
+verboseCheck p = quickCheck (verbose p)
+
+-- | Tests a property, using test arguments, and prints the results and all test cases generated to 'stdout'.
+-- This is just a convenience function that combines 'quickCheckWith' and 'verbose'.
+verboseCheckWith :: Testable prop => Args -> prop -> IO ()
+verboseCheckWith args p = quickCheckWith args (verbose p)
+
+-- | Tests a property, produces a test result, and prints the results and all test cases generated to 'stdout'.
+-- This is just a convenience function that combines 'quickCheckResult' and 'verbose'.
+verboseCheckResult :: Testable prop => prop -> IO Result
+verboseCheckResult p = quickCheckResult (verbose p)
+
+-- | Tests a property, using test arguments, produces a test result, and prints the results and all test cases generated to 'stdout'.
+-- This is just a convenience function that combines 'quickCheckWithResult' and 'verbose'.
+verboseCheckWithResult :: Testable prop => Args -> prop -> IO Result
+verboseCheckWithResult a p = quickCheckWithResult a (verbose p)
 
 --------------------------------------------------------------------------
 -- main test loop
@@ -195,34 +215,27 @@ runATest st f =
        ++ ")"
         )
      let size = computeSize st (numSuccessTests st) (numDiscardedTests st)
-     (mres, ts) <- unpackRose (unProp (f rnd1 size))
-     res <- mres
+     MkRose res ts <- protectRose (reduceRose (unProp (f rnd1 size)))
      callbackPostTest st res
      
-     case ok res of
-       Just True -> -- successful test
+     case res of
+       MkResult{ok = Just True, stamp = stamp, expect = expect} -> -- successful test
          do test st{ numSuccessTests = numSuccessTests st + 1
                    , randomSeed      = rnd2
-                   , collected       = stamp res : collected st
-                   , expectedFailure = expect res
+                   , collected       = stamp : collected st
+                   , expectedFailure = expect
                    } f
        
-       Nothing -> -- discarded test
+       MkResult{ok = Nothing, expect = expect} -> -- discarded test
          do test st{ numDiscardedTests = numDiscardedTests st + 1
                    , randomSeed        = rnd2
-                   , expectedFailure   = expect res
+                   , expectedFailure   = expect
                    } f
          
-       Just False -> -- failed test
+       MkResult{ok = Just False} -> -- failed test
          do if expect res
               then putPart (terminal st) (bold "*** Failed! ")
               else putPart (terminal st) "+++ OK, failed as expected. "
-            putTemp (terminal st)
-              ( short 30 (P.reason res)
-             ++ " (after "
-             ++ number (numSuccessTests st+1) "test"
-             ++ ")..."
-              )
             numShrinks <- foundFailure st res ts
             theOutput <- terminalOutput (terminal st)
             if not (expect res) then
@@ -300,40 +313,39 @@ success st =
 --------------------------------------------------------------------------
 -- main shrinking loop
 
-foundFailure :: State -> P.Result -> [Rose (IO P.Result)] -> IO Int
+foundFailure :: State -> P.Result -> [Rose P.Result] -> IO Int
 foundFailure st res ts =
   do localMin st{ numTryShrinks = 0 } res ts
 
-localMin :: State -> P.Result -> [Rose (IO P.Result)] -> IO Int
+localMin :: State -> P.Result -> [Rose P.Result] -> IO Int
 localMin st res _ | P.interrupted res = localMinFound st res
 localMin st res ts = do
+  putTemp (terminal st)
+    ( short 26 (P.reason res)
+   ++ " (after " ++ number (numSuccessTests st+1) "test"
+   ++ concat [ " and "
+            ++ show (numSuccessShrinks st)
+            ++ concat [ "." ++ show (numTryShrinks st) | numTryShrinks st > 0 ]
+            ++ " shrink"
+            ++ (if numSuccessShrinks st == 1
+                && numTryShrinks st == 0
+                then "" else "s")
+             | numSuccessShrinks st > 0 || numTryShrinks st > 0
+             ]
+   ++ ")..."
+    )
   r <- tryEvaluate ts
   case r of
     Left err ->
       localMinFound st
-         (exception "Exception while generating shrink-list" err)
+         (exception "Exception while generating shrink-list" err) { callbacks = callbacks res }
     Right ts' -> localMin' st res ts'
 
-localMin' :: State -> P.Result -> [Rose (IO P.Result)] -> IO Int
+localMin' :: State -> P.Result -> [Rose P.Result] -> IO Int
 localMin' st res [] = localMinFound st res
 localMin' st res (t:ts) =
   do -- CALLBACK before_test
-    (mres', ts') <- unpackRose t
-    res' <- mres'
-    putTemp (terminal st)
-      ( short 35 (P.reason res)
-     ++ " (after " ++ number (numSuccessTests st+1) "test"
-     ++ concat [ " and "
-              ++ show (numSuccessShrinks st)
-              ++ concat [ "." ++ show (numTryShrinks st) | numTryShrinks st > 0 ]
-              ++ " shrink"
-              ++ (if numSuccessShrinks st == 1
-                  && numTryShrinks st == 0
-                  then "" else "s")
-               | numSuccessShrinks st > 0 || numTryShrinks st > 0
-               ]
-     ++ ")..."
-      )
+    MkRose res' ts' <- protectRose (reduceRose t)
     callbackPostTest st res'
     if ok res' == Just False
       then foundFailure st{ numSuccessShrinks = numSuccessShrinks st + 1 } res' ts'
@@ -357,11 +369,21 @@ localMinFound st res =
 
 callbackPostTest :: State -> P.Result -> IO ()
 callbackPostTest st res =
-  sequence_ [ f st res | PostTest f <- callbacks res ]
+  sequence_ [ safely st (f st res) | PostTest _ f <- callbacks res ]
 
 callbackPostFinalFailure :: State -> P.Result -> IO ()
 callbackPostFinalFailure st res =
-  sequence_ [ f st res | PostFinalFailure f <- callbacks res ]
+  sequence_ [ safely st (f st res) | PostFinalFailure _ f <- callbacks res ]
+
+safely :: State -> IO () -> IO ()
+safely st x = do
+  r <- tryEvaluateIO x
+  case r of
+    Left e ->
+      putLine (terminal st)
+        ("*** Exception in callback: " ++ show e)
+    Right x ->
+      return x
 
 --------------------------------------------------------------------------
 -- the end.
index d609e28..99fdb85 100644 (file)
@@ -40,7 +40,7 @@ newtype Str = MkStr String
 instance Show Str where
   show (MkStr s) = s
 
-ranges :: Integral a => a -> a -> Str
+ranges :: (Show a, Integral a) => a -> a -> Str
 ranges k n = MkStr (show n' ++ " -- " ++ show (n'+k-1))
  where
   n' = k * (n `div` k)
index a280694..a3dabd4 100644 (file)
@@ -1,4 +1,4 @@
-{-# LANGUAGE TypeOperators, TypeSynonymInstances #-}
+{-# LANGUAGE TypeOperators, TypeSynonymInstances, FlexibleInstances #-}
 
 -- | Strict complex doubles.
 module Data.Array.Repa.Algorithms.Complex
@@ -7,7 +7,8 @@ module Data.Array.Repa.Algorithms.Complex
        , arg)
 where
 
--- | Strict complex doubles.
+
+-- | Complex doubles.
 type Complex 
        = (Double, Double)
 
index ea5ecb9..41b3389 100644 (file)
@@ -1,5 +1,11 @@
 {-# LANGUAGE BangPatterns, PackageImports #-}
 {-# OPTIONS -Wall -fno-warn-missing-signatures -fno-warn-incomplete-patterns #-}
+
+-- | Old support for stencil based convolutions. 
+--
+--   NOTE: This is slated to be merged with the new Stencil support in the next version
+--         of Repa. We'll still expose the `convolve` function though.
+--
 module Data.Array.Repa.Algorithms.Convolve
        ( convolve
 
@@ -9,10 +15,9 @@ module Data.Array.Repa.Algorithms.Convolve
        , convolveOut )
 where
 import Data.Array.Repa                                         as A
-import qualified Data.Array.Repa.Shape                 as S
 import qualified Data.Vector.Unboxed                   as V
+import qualified Data.Array.Repa.Shape                 as S
 import Prelude                                         as P
-import "dph-prim-par" Data.Array.Parallel.Unlifted     (Elt)
 
 
 -- Plain Convolve ---------------------------------------------------------------------------------
@@ -27,11 +32,11 @@ convolve
 
 {-# INLINE convolve #-}
 convolve makeOut
-       kernel@(Manifest       (_ :. krnHeight :. krnWidth) krnVec)
-        image@(Manifest imgSh@(_ :. imgHeight :. imgWidth) imgVec)
+       kernel@(Array       (_ :. krnHeight :. krnWidth) [Region RangeAll (GenManifest krnVec)])
+        image@(Array imgSh@(_ :. imgHeight :. imgWidth) [Region RangeAll (GenManifest imgVec)])
 
  = kernel `deepSeqArray` image `deepSeqArray` 
-   force $ traverse image id update
+   force $ unsafeTraverse image id update
  where 
        !krnHeight2     = krnHeight `div` 2
        !krnWidth2      = krnWidth  `div` 2
@@ -117,11 +122,11 @@ convolveOut
 
 {-# INLINE convolveOut #-}
 convolveOut getOut
-       kernel@(Manifest krnSh@(_ :. krnHeight :. krnWidth) _)
-        image@(Manifest imgSh@(_ :. imgHeight :. imgWidth) _)
+       kernel@(Array krnSh@(_ :. krnHeight :. krnWidth) _)
+        image@(Array imgSh@(_ :. imgHeight :. imgWidth) _)
 
  = kernel `deepSeqArray` image `deepSeqArray` 
-   force $ traverse image id stencil
+   force $ unsafeTraverse image id stencil
  where 
        !krnHeight2     = krnHeight `div` 2
        !krnWidth2      = krnWidth  `div` 2
@@ -155,7 +160,7 @@ convolveOut getOut
                 | otherwise
                 = let  !ix@(sh :. y :. x)      = S.fromIndex krnSh count
                        !ix'                    = sh :. y + jkrnHeight' :. x + ikrnWidth'
-                       !here                   = kernel !: ix * (get' ix')
+                       !here                   = kernel `unsafeIndex` ix * (get' ix')
                   in   integrate (count + 1) (acc + here)
 
           in   integrate 0 0
index 0be1662..37bf0b4 100644 (file)
@@ -22,6 +22,7 @@ where
 import Data.Array.Repa.Algorithms.DFT.Roots
 import Data.Array.Repa.Algorithms.Complex
 import Data.Array.Repa                         as A
+import Prelude                                 as P
 
 -- | Compute the DFT along the low order dimension of an array.
 dft    :: forall sh
@@ -62,8 +63,8 @@ dftWithRoots rofu arr
        | _ :. rLen     <- extent rofu
        , _ :. vLen     <- extent arr
        , rLen /= vLen
-       = error $  "dftWithRoots: length of vector (" ++ show vLen ++ ")"
-               ++ " does not match the length of the roots (" ++ show rLen ++ ")"
+       = error $    "dftWithRoots: length of vector (" P.++ show vLen P.++ ")"
+               P.++ " does not match the length of the roots (" P.++ show rLen P.++ ")"
 
        | otherwise
        = traverse arr id (\_ k -> dftWithRootsSingle rofu arr k)
@@ -84,8 +85,8 @@ dftWithRootsSingle rofu arrX (_ :. k)
        | _ :. rLen     <- extent rofu
        , _ :. vLen     <- extent arrX
        , rLen /= vLen
-       = error $  "dftWithRootsSingle: length of vector (" ++ show vLen ++ ")"
-               ++ " does not match the length of the roots (" ++ show rLen ++ ")"
+       = error $    "dftWithRootsSingle: length of vector (" P.++ show vLen P.++ ")"
+               P.++ " does not match the length of the roots (" P.++ show rLen P.++ ")"
 
        | otherwise
        = let   sh@(_ :. len)   = extent arrX
@@ -93,7 +94,7 @@ dftWithRootsSingle rofu arrX (_ :. k)
                -- All the roots we need to multiply with.
                wroots          = fromFunction sh elemFn
                elemFn (sh' :. n) 
-                       = rofu !: (sh' :. (k * n) `mod` len)
+                       = rofu ! (sh' :. (k * n) `mod` len)
 
          in  A.sumAll $ A.zipWith (*) arrX wroots
 
index 8c2c40c..bed9f1f 100644 (file)
@@ -1,18 +1,22 @@
 {-# LANGUAGE TypeOperators, PatternGuards, RankNTypes, ScopedTypeVariables, BangPatterns, FlexibleContexts #-}
 {-# OPTIONS -fno-warn-incomplete-patterns #-}
 
--- | Fast computation of Discrete Fourier Transforms using the Cooley-Tuckey algorithm.
+-- | Fast computation of Discrete Fourier Transforms using the Cooley-Tuckey algorithm. 
+--   Time complexity is O(n log n) in the size of the input. 
 --
---   Time complexity is O(n log n) in the size of the input.
+--   This uses a naive divide-and-conquer algorithm, the absolute performance is about
+--   50x slower than FFTW in estimate mode.
 --
 module Data.Array.Repa.Algorithms.FFT
        ( Mode(..)
+       , isPowerOfTwo
        , fft3d
        , fft2d
        , fft1d)
 where
 import Data.Array.Repa.Algorithms.Complex
 import Data.Array.Repa                         as A
+import Prelude                                  as P
 
 data Mode
        = Forward
@@ -29,23 +33,42 @@ signOfMode mode
        Inverse         ->   1
 
 
+{-# INLINE isPowerOfTwo #-}
+-- | Check if an `Int` is a power of two.
+isPowerOfTwo :: Int -> Bool
+isPowerOfTwo n
+       | 0     <- n            = True
+       | 2     <- n            = True
+       | n `mod` 2 == 0        = isPowerOfTwo (n `div` 2)
+       | otherwise             = False
+
+
+
+
 -- 3D Transform -----------------------------------------------------------------------------------
--- | Compute the DFT of a 3d array.
+-- | Compute the DFT of a 3d array. Array dimensions must be powers of two else `error`.
 fft3d  :: Mode
        -> Array DIM3 Complex
        -> Array DIM3 Complex
 
 fft3d mode arr
  = let _ :. depth :. height :. width   = extent arr
-       sign    = signOfMode mode
-       scale   = fromIntegral (depth * width * height) 
+       !sign   = signOfMode mode
+       !scale  = fromIntegral (depth * width * height) 
                
-   in  arr `deepSeqArray` 
-       case mode of
-               Forward -> fftTrans3d sign $ fftTrans3d sign $ fftTrans3d sign arr
-               Reverse -> fftTrans3d sign $ fftTrans3d sign $ fftTrans3d sign arr
-               Inverse -> force $ A.map (/ scale) 
-                               $ fftTrans3d sign $ fftTrans3d sign $ fftTrans3d sign arr
+   in  if not (isPowerOfTwo depth && isPowerOfTwo height && isPowerOfTwo width)
+        then error $ unlines
+               [ "Data.Array.Repa.Algorithms.FFT: fft3d"
+               , "  Array dimensions must be powers of two,"
+               , "  but the provided array is " 
+                       P.++ show height P.++ "x" P.++ show width P.++ "x" P.++ show depth ]
+                  
+        else arr `deepSeqArray` 
+               case mode of
+                       Forward -> fftTrans3d sign $ fftTrans3d sign $ fftTrans3d sign arr
+                       Reverse -> fftTrans3d sign $ fftTrans3d sign $ fftTrans3d sign arr
+                       Inverse -> force $ A.map (/ scale) 
+                                       $ fftTrans3d sign $ fftTrans3d sign $ fftTrans3d sign arr
 
 fftTrans3d 
        :: Double
@@ -69,7 +92,7 @@ rotate3d arr
 
 
 -- Matrix Transform -------------------------------------------------------------------------------
--- | Compute the DFT of a matrix.
+-- | Compute the DFT of a matrix. Array dimensions must be powers of two else `error`.
 fft2d  :: Mode
        -> Array DIM2 Complex
        -> Array DIM2 Complex
@@ -79,11 +102,17 @@ fft2d mode arr
        sign    = signOfMode mode
        scale   = fromIntegral (width * height) 
                
-   in  arr `deepSeqArray` 
-       case mode of
-               Forward -> fftTrans2d sign $ fftTrans2d sign arr
-               Reverse -> fftTrans2d sign $ fftTrans2d sign arr
-               Inverse -> force $ A.map (/ scale) $ fftTrans2d sign $ fftTrans2d sign arr
+   in  if not (isPowerOfTwo height && isPowerOfTwo width)
+        then error $ unlines
+               [ "Data.Array.Repa.Algorithms.FFT: fft2d"
+               , "  Array dimensions must be powers of two,"
+               , "  but the provided array is " P.++ show height P.++ "x" P.++ show width ]
+        
+        else arr `deepSeqArray` 
+               case mode of
+                       Forward -> fftTrans2d sign $ fftTrans2d sign arr
+                       Reverse -> fftTrans2d sign $ fftTrans2d sign arr
+                       Inverse -> force $ A.map (/ scale) $ fftTrans2d sign $ fftTrans2d sign arr
 
 fftTrans2d 
        :: Double
@@ -98,7 +127,7 @@ fftTrans2d sign arr'
 
 
 -- Vector Transform -------------------------------------------------------------------------------
--- | Compute the DFT of a vector.
+-- | Compute the DFT of a vector. Array dimensions must be powers of two else `error`.
 fft1d  :: Mode 
        -> Array DIM1 Complex 
        -> Array DIM1 Complex
@@ -108,11 +137,17 @@ fft1d mode arr
        sign    = signOfMode mode
        scale   = fromIntegral len
        
-   in  arr `deepSeqArray`
-       case mode of
-               Forward -> fftTrans1d sign arr
-               Reverse -> fftTrans1d sign arr
-               Inverse -> force $ A.map (/ scale) $ fftTrans1d sign arr
+   in  if not $ isPowerOfTwo len
+        then error $ unlines 
+                [ "Data.Array.Repa.Algorithms.FFT: fft1d"
+               , "  Array dimensions must be powers of two, "
+                , "  but the provided array is " P.++ show len ]
+             
+        else arr `deepSeqArray`
+               case mode of
+                       Forward -> fftTrans1d sign arr
+                       Reverse -> fftTrans1d sign arr
+                       Inverse -> force $ A.map (/ scale) $ fftTrans1d sign arr
 
 fftTrans1d
        :: Double 
@@ -141,14 +176,15 @@ fft !sign !sh !lenVec !vec
 
         where  swivel (sh' :. ix)
                 = case ix of
-                       0       -> vec !: (sh' :. offset) + vec !: (sh' :. (offset + stride))
-                       1       -> vec !: (sh' :. offset) - vec !: (sh' :. (offset + stride))
+                       0       -> (vec `unsafeIndex` (sh' :. offset)) + (vec `unsafeIndex` (sh' :. (offset + stride)))
+                       1       -> (vec `unsafeIndex` (sh' :. offset)) - (vec `unsafeIndex` (sh' :. (offset + stride)))
 
                {-# INLINE combine #-}
-               combine !len' evens@Manifest{} odds@Manifest{}
+               combine !len'   evens@(Array _ [Region RangeAll GenManifest{}]) 
+                                odds@(Array _ [Region RangeAll GenManifest{}])
                 = evens `deepSeqArray` odds `deepSeqArray`
-                  let  odds'   = traverse odds id (\get ix@(_ :. k) -> twiddle sign k len' * get ix) 
-                  in   force   $ (evens +^ odds') +:+ (evens -^ odds')
+                  let  odds'   = unsafeTraverse odds id (\get ix@(_ :. k) -> twiddle sign k len' * get ix) 
+                  in   force   $ (evens +^ odds') A.++ (evens -^ odds')
 
 
 -- Compute a twiddle factor.
@@ -163,4 +199,3 @@ twiddle sign k' n'
        where   k       = fromIntegral k'
                n       = fromIntegral n'
       
-
diff --git a/fibon/Repa/_RepaLib/repa-algorithms/Data/Array/Repa/Algorithms/Iterate.hs b/fibon/Repa/_RepaLib/repa-algorithms/Data/Array/Repa/Algorithms/Iterate.hs
new file mode 100644 (file)
index 0000000..3853cc9
--- /dev/null
@@ -0,0 +1,45 @@
+{-# LANGUAGE BangPatterns, PackageImports #-}
+{-# OPTIONS -Wall -fno-warn-missing-signatures -fno-warn-incomplete-patterns #-}
+module Data.Array.Repa.Algorithms.Iterate
+       (iterateBlockwise, iterateBlockwise')
+where
+import Data.Array.Repa
+
+-- | Iterate array transformation function a fixed number of times, applying `force2` between
+--   each iteration.
+iterateBlockwise
+       :: (Elt a, Num a)
+       => Int                                  -- ^ Number of iterations to run for.
+       -> (Array DIM2 a -> Array DIM2 a)       -- ^ Fn to step the array.
+       -> Array DIM2 a                         -- ^ Initial array value.
+       -> Array DIM2 a
+       
+{-# INLINE iterateBlockwise #-}
+iterateBlockwise steps f arrInit@(Array shInit [Region RangeAll (GenManifest vecInit)])
+ = arrInit `deepSeqArray`
+  goSolve steps shInit vecInit
+
+ where -- NOTE: We manually unpack the current array into its shape and vector to
+       --       stop GHC from unboxing the vector again for every loop. deepSeqing
+       --       the arrays at the start of solveLaplace makes the unboxings happen
+       --       at that point in the corresponding core code.
+       goSolve !i !shCurrent !vecCurrent
+        = let  !arrCurrent     = fromVector shCurrent vecCurrent
+          in   if i == 0 
+                then arrCurrent
+                else let arrNew@(Array _ [Region RangeAll (GenManifest _)]) 
+                               = force2 $ f arrCurrent
+                     in  goSolve (i - 1) (extent arrNew) (toVector arrNew)
+
+
+-- | As above, but with the parameters flipped.
+iterateBlockwise'
+       :: (Elt a, Num a)
+       => Int                                  -- ^ Number of iterations to run for.
+       -> Array DIM2 a                         -- ^ Initial array value.
+       -> (Array DIM2 a -> Array DIM2 a)       -- ^ Fn to step the array.
+       -> Array DIM2 a
+
+{-# INLINE iterateBlockwise' #-}
+iterateBlockwise' steps arr fn
+       = iterateBlockwise steps fn arr
\ No newline at end of file
index 66ff52a..7f230ab 100644 (file)
@@ -13,35 +13,29 @@ module Data.Array.Repa.Algorithms.Matrix
        (multiplyMM)
 where
 import Data.Array.Repa as A
-import "dph-prim-par" Data.Array.Parallel.Unlifted                     (Elt)
 
 -- | Matrix-matrix multiply.
 multiplyMM
-       :: (Num a, Elt a)
-       => Array DIM2 a
-       -> Array DIM2 a
-       -> Array DIM2 a
+       :: Array DIM2 Double
+       -> Array DIM2 Double
+       -> Array DIM2 Double
 
--- NOTE: We need this INLINE so multiplyMM gets specialised for the appropriate types.
-{-# INLINE multiplyMM #-}
-multiplyMM arr1 arr2
- = multiplyMM' (force arr1) (force arr2)
-
- where multiplyMM' arr@Manifest{} brr@Manifest{}
-       = A.force $ A.sum (A.zipWith (*) arrRepl brrRepl)
-       where
-        trr            = force $ transpose2D brr
-        arrRepl        = A.replicate (Z :. All   :. colsB :. All) arr
-        brrRepl        = A.replicate (Z :. rowsA :. All   :. All) trr
-        (Z :. _     :. rowsA) = extent arr
-        (Z :. colsB :. _    ) = extent brr
+{-# NOINLINE multiplyMM #-}
+multiplyMM arr@(Array _ [Region RangeAll (GenManifest _)])
+          brr@(Array _ [Region RangeAll (GenManifest _)])
+ = [arr, brr] `deepSeqArrays`
+   A.force $ A.sum (A.zipWith (*) arrRepl brrRepl)
+ where trr@(Array _ [Region RangeAll (GenManifest _)])
+                       = force $ transpose2D brr
+       arrRepl         = trr `deepSeqArray` A.extend (Z :. All   :. colsB :. All) arr
+       brrRepl         = trr `deepSeqArray` A.extend (Z :. rowsA :. All   :. All) trr
+       (Z :. _     :. rowsA) = extent arr
+       (Z :. colsB :. _    ) = extent brr
        
 
 transpose2D :: Elt e => Array DIM2 e -> Array DIM2 e
 {-# INLINE transpose2D #-}
 transpose2D arr
  = backpermute new_extent swap arr
- where
-       swap (Z :. i :. j)      = Z :. j :. i
+ where swap (Z :. i :. j)      = Z :. j :. i
        new_extent              = swap (extent arr)
-
diff --git a/fibon/Repa/_RepaLib/repa-algorithms/Data/Array/Repa/Algorithms/Randomish.hs b/fibon/Repa/_RepaLib/repa-algorithms/Data/Array/Repa/Algorithms/Randomish.hs
new file mode 100644 (file)
index 0000000..ac0b592
--- /dev/null
@@ -0,0 +1,115 @@
+{-# LANGUAGE BangPatterns #-}
+
+module Data.Array.Repa.Algorithms.Randomish
+       ( randomishIntArray
+       , randomishIntVector
+       , randomishDoubleArray
+       , randomishDoubleVector)
+where
+import Data.Word
+import Data.Vector.Unboxed                     (Vector)
+import Data.Array.Repa                         as R
+import qualified Data.Vector.Unboxed.Mutable   as MV
+import qualified Data.Vector.Unboxed           as V
+import qualified Data.Vector.Generic           as G
+
+
+-- | Use the ''minimal standard'' Lehmer generator to quickly generate some random
+--   numbers with reasonable statistical properties. By ''reasonable'' we mean good
+--   enough for games and test data, but not cryptography or anything where the
+--   quality of the randomness really matters. 
+--
+--   By nature of the algorithm, the maximum value in the output is clipped
+--   to (valMin + 2^31 - 1)
+-- 
+--   From ''Random Number Generators: Good ones are hard to find''
+--   Stephen K. Park and Keith W. Miller.
+--   Communications of the ACM, Oct 1988, Volume 31, Number 10.
+--
+randomishIntArray
+       :: Shape sh
+       => sh                   -- ^ Shape of array
+       -> Int                  -- ^ Minumum value in output.
+       -> Int                  -- ^ Maximum value in output.
+       -> Int                  -- ^ Random seed.       
+       -> Array sh Int         -- ^ Array of randomish numbers.
+
+randomishIntArray !sh !valMin !valMax !seed
+       = fromVector sh $ randomishIntVector (R.size sh) valMin valMax seed
+
+
+randomishIntVector 
+       :: Int                  -- ^ Length of vector.
+       -> Int                  -- ^ Minumum value in output.
+       -> Int                  -- ^ Maximum value in output.
+       -> Int                  -- ^ Random seed.       
+       -> Vector Int           -- ^ Vector of randomish numbers.
+
+randomishIntVector !len !valMin' !valMax' !seed'
+ = let -- a magic number
+       -- (don't change it, the randomness depends on this specific number).
+       multiplier :: Word64
+       multiplier = 16807
+
+       -- a merzenne prime
+       -- (don't change it, the randomness depends on this specific number).
+       modulus :: Word64
+       modulus = 2^(31 :: Integer) - 1
+
+       -- if the seed is 0 all the numbers in the sequence are also 0.
+       seed    
+        | seed' == 0   = 1
+        | otherwise    = seed'
+
+       !valMin = fromIntegral valMin'
+       !valMax = fromIntegral valMax' + 1
+       !range  = valMax - valMin
+
+       {-# INLINE f #-}
+       f x             = multiplier * x `mod` modulus
+ in G.create 
+     $ do      
+       vec             <- MV.new len
+
+       let go !ix !x 
+               | ix == len     = return ()
+               | otherwise
+               = do    let x'  = f x
+                       MV.write vec ix $ fromIntegral $ (x `mod` range) + valMin
+                       go (ix + 1) x'
+
+       go 0 (f $ f $ f $ fromIntegral seed)
+       return vec
+
+
+-- | Generate some randomish doubles with terrible statistical properties.
+--   This just takes randomish ints then scales them, so there's not much randomness in low-order bits.
+randomishDoubleArray
+       :: Shape sh
+       => sh                   -- ^ Shape of array
+       -> Double               -- ^ Minumum value in output.
+       -> Double               -- ^ Maximum value in output.
+       -> Int                  -- ^ Random seed.       
+       -> Array sh Double      -- ^ Array of randomish numbers.
+
+randomishDoubleArray !sh !valMin !valMax !seed
+       = fromVector sh $ randomishDoubleVector (R.size sh) valMin valMax seed
+
+
+-- | Generate some randomish doubles with terrible statistical properties.
+--   This just takes randmish ints then scales them, so there's not much randomness in low-order bits.
+randomishDoubleVector
+       :: Int                  -- ^ Length of vector
+       -> Double               -- ^ Minimum value in output
+       -> Double               -- ^ Maximum value in output
+       -> Int                  -- ^ Random seed.
+       -> Vector Double        -- ^ Vector of randomish doubles.
+
+randomishDoubleVector !len !valMin !valMax !seed
+ = let range   = valMax - valMin
+
+       mx      = 2^(30 :: Integer) - 1
+       mxf     = fromIntegral mx
+       ints    = randomishIntVector len 0 mx seed
+       
+   in  V.map (\n -> valMin + (fromIntegral n / mxf) * range) ints
index 6f96364..af2e912 100644 (file)
@@ -1,5 +1,5 @@
 Name:                repa-algorithms
-Version:             1.1.1.0
+Version:             2.2.0.1
 License:             BSD3
 License-file:        LICENSE
 Author:              The DPH Team
@@ -8,24 +8,19 @@ 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
+Homepage:            http://repa.ouroborus.net
+Bug-reports:         repa@ouroborus.net
 Description:
-        NOTE: You must use the GHC head branch > 6.13.20100309 to get decent performance.
         Reusable algorithms using the Repa array library.
 
 Synopsis:
         Algorithms using the Repa array library.
 
-Tested-with: GHC == 6.13.20100309, GHC == 6.12.1
-
 Library
   Build-Depends: 
         base                 == 4.*,
-        dph-base             == 0.5.*,
-        dph-prim-par         == 0.5.*,
-        vector               == 0.7.*,
-        repa                 == 1.1.*
+        vector               == 0.9.*,
+        repa                 == 2.2.*
 
   ghc-options:
         -Wall -fno-warn-missing-signatures
@@ -38,9 +33,11 @@ Library
 
   Exposed-modules:
         Data.Array.Repa.Algorithms.Complex
+        Data.Array.Repa.Algorithms.Randomish
         Data.Array.Repa.Algorithms.DFT
         Data.Array.Repa.Algorithms.DFT.Roots
         Data.Array.Repa.Algorithms.DFT.Center
         Data.Array.Repa.Algorithms.FFT
         Data.Array.Repa.Algorithms.Matrix
         Data.Array.Repa.Algorithms.Convolve
+        Data.Array.Repa.Algorithms.Iterate
index 5245279..e02fc81 100644 (file)
@@ -3,7 +3,9 @@
 -- | Conversions between Repa Arrays and ByteStrings.
 module Data.Array.Repa.ByteString
        ( toByteString
-       , fromByteString)
+       , fromByteString
+       , copyFromPtrWord8
+       , copyToPtrWord8)
 where
 import Data.Word
 import Data.Array.Repa
@@ -11,29 +13,31 @@ import Foreign.Marshal.Alloc
 import Foreign.Ptr
 import Foreign.Storable
 import System.IO.Unsafe
+import qualified Data.Vector.Unboxed   as V
 import qualified Data.ByteString       as BS
 import Data.ByteString                 (ByteString)
-import qualified "dph-prim-par" Data.Array.Parallel.Unlifted as U
 
 
--- | Convert an `Array` to a (strict) `ByteString`.
+-- | Copy an `Array` to a (strict) `ByteString`.
 toByteString 
        :: Shape sh
        => Array sh Word8
        -> ByteString
 
+{-# NOINLINE toByteString #-}
 toByteString arr
- = unsafePerformIO
- $ allocaBytes (size $ extent arr)     $ \(bufDest :: Ptr Word8) ->
-   let         uarr    = toUArray arr
-       len     = size $ extent arr
+ =  withManifest' (force arr) $ \arr' 
+ -> unsafePerformIO
+ $ allocaBytes (size $ extent arr')    $ \(bufDest :: Ptr Word8) ->
+   let         vec     = toVector arr'
+       len     = size $ extent arr'
 
        copy offset
         | offset >= len
         = return ()
 
         | otherwise
-        = do   pokeByteOff bufDest offset (uarr U.!: offset)
+        = do   pokeByteOff bufDest offset (vec `V.unsafeIndex` offset)
                copy (offset + 1)
 
     in do
@@ -41,8 +45,8 @@ toByteString arr
        BS.packCStringLen (castPtr bufDest, len)
 
 
--- | Convert a (strict) `ByteString` to an `Array`.
---     The given array size must match the length of the `ByteString`, else `error`.
+-- | Copy a (strict) `ByteString` to a new `Array`.
+--     The given array extent must match the length of the `ByteString`, else `error`.
 fromByteString 
        :: Shape sh
        => sh
@@ -58,3 +62,40 @@ fromByteString sh str
        | otherwise
        = fromFunction sh (\ix -> str `BS.index` toIndex sh ix)
 
+
+-- Ptr utils ------------------------------------------------------------------
+-- | Copy some data from somewhere into a new `Array`.
+copyFromPtrWord8 
+       :: Shape sh
+       => sh
+       -> Ptr Word8
+       -> IO (Array sh Word8)
+
+{-# INLINE copyFromPtrWord8 #-}        
+copyFromPtrWord8 sh ptr
+ = do  return  $ fromFunction sh (\ix -> unsafePerformIO (peekElemOff ptr (toIndex sh ix)))
+
+
+-- | Copy array data somewhere.s
+copyToPtrWord8 
+       :: Shape sh
+       => Ptr Word8
+       -> Array sh Word8
+       -> IO ()
+       
+{-# INLINE copyToPtrWord8 #-}
+copyToPtrWord8 ptr arr
+ = let vec     = toVector arr
+       len     = size $ extent arr
+       
+       copy offset
+        | offset >= len
+        = return ()
+       
+        | otherwise
+        = do   pokeByteOff ptr offset (vec `V.unsafeIndex` offset)
+               copy (offset + 1)
+
+   in do
+       copy 0
+       return ()
index d614120..7fdec4f 100644 (file)
@@ -1,5 +1,5 @@
 Name:                repa-bytestring
-Version:             1.1.1.0
+Version:             2.2.0.1
 License:             BSD3
 License-file:        LICENSE
 Author:              The DPH Team
@@ -8,23 +8,20 @@ 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
+Homepage:            http://repa.ouroborus.net
+Bug-reports:         repa@ouroborus.net
 Description:
-        NOTE: You must use the GHC head branch > 6.13.20100309 to get decent performance.
         Conversions between Repa Arrays and ByteStrings.
 
 Synopsis:
         Conversions between Repa Arrays and ByteStrings.
 
-Tested-with: GHC == 6.13.20100309, GHC == 6.12.1
-
 Library
   Build-Depends: 
         base                 == 4.*,
-        dph-prim-par         == 0.5.*,
-        repa                 == 1.1.*,
-        bytestring           == 0.9.1.*
+        repa                 == 2.2.*,
+        vector               == 0.9.*,
+        bytestring           == 0.9.*
 
   ghc-options:
         -Odph -Wall -fno-warn-missing-signatures
index c4afce2..fe89663 100644 (file)
@@ -1,15 +1,18 @@
-{-# LANGUAGE PackageImports #-} 
+{-# LANGUAGE PackageImports, PatternGuards, ExplicitForAll  #-} 
 
 -- | Reading and writing arrays as uncompressed 24 and 32 bit Windows BMP files.
 module Data.Array.Repa.IO.BMP
        ( readImageFromBMP
        , readComponentsFromBMP
+       , readComponentsListFromBMP
        , readMatrixFromGreyscaleBMP
+
+       -- Writing.
        , writeImageToBMP
        , writeComponentsToBMP
+       , writeComponentsListToBMP
        , writeMatrixToGreyscaleBMP)
 where
-import qualified "dph-prim-par" Data.Array.Parallel.Unlifted   as U
 import Data.Array.Repa                         as A
 import Data.Array.Repa.ByteString               as A
 import Prelude                                 as P
@@ -35,13 +38,28 @@ readMatrixFromGreyscaleBMP filePath
        case eComps of 
         Left err       -> return $ Left err
         Right (arrRed, arrGreen, arrBlue)
-         -> let arr    = force 
+         -> let arr    = force2 
                        $ A.fromFunction (extent arrRed)
-                          (\ix -> sqrt ( (fromIntegral (arrRed   !: ix) / 255) ^ (2 :: Int)
-                                       + (fromIntegral (arrGreen !: ix) / 255) ^ (2 :: Int)
-                                       + (fromIntegral (arrBlue  !: ix) / 255) ^ (2 :: Int)))
+                          (\ix -> sqrt ( (fromIntegral (arrRed   ! ix) / 255) ^ (2 :: Int)
+                                       + (fromIntegral (arrGreen ! ix) / 255) ^ (2 :: Int)
+                                       + (fromIntegral (arrBlue  ! ix) / 255) ^ (2 :: Int)))
             in arr `deepSeqArray` return (Right arr)
-               
+
+
+-- | Like `readComponentsFromBMP`, but return the components as a list.
+readComponentsListFromBMP
+       :: FilePath
+       -> IO (Either Error [Array DIM2 Word8])
+
+readComponentsListFromBMP filePath
+ = do  eComps  <- readComponentsFromBMP filePath
+       case eComps of
+        Left err
+         -> return $ Left err
+
+        Right (arrRed, arrGreen, arrBlue)      
+         -> return $ Right [arrRed, arrGreen, arrBlue]
+
 
 -- | Read RGB components from a BMP file.
 --     Returns arrays of red, green and blue components, all with the same extent.
@@ -66,18 +84,18 @@ readComponentsFromBMP' bmp
        shapeFn _       = Z :. height :. width
 
        arrRed  
-        = traverse arr shapeFn
+        = force2 $ traverse arr shapeFn
                (\get (sh :. x) -> get (sh :. (x * 4)))
 
        arrGreen
-        = traverse arr shapeFn
+        = force2 $ traverse arr shapeFn
                (\get (sh :. x) -> get (sh :. (x * 4 + 1)))
 
        arrBlue
-        = traverse arr shapeFn
+        = force2 $ traverse arr shapeFn
                (\get (sh :. x) -> get (sh :. (x * 4 + 2)))
        
-   in  (arrRed, arrGreen, arrBlue)
+   in  [arrRed, arrGreen, arrBlue] `deepSeqArrays` (arrRed, arrGreen, arrBlue)
 
 
 -- | Read a RGBA image from a BMP file.
@@ -109,20 +127,36 @@ readImageFromBMP' bmp
 --     Negative values are discarded. Positive values are normalised to the maximum 
 --     value in the matrix and used as greyscale pixels.
 writeMatrixToGreyscaleBMP 
-       :: FilePath
-       -> Array DIM2 Double
+       :: forall a. (Num a, Elt a, Fractional a, RealFrac a)
+       => FilePath
+       -> Array DIM2 a
        -> IO ()
 
-{-# NOINLINE writeMatrixToGreyscaleBMP #-}
+{-# NOINLINE   writeMatrixToGreyscaleBMP #-}
+{-# SPECIALISE writeMatrixToGreyscaleBMP :: FilePath -> Array DIM2 Float  -> IO () #-}
+{-# SPECIALISE writeMatrixToGreyscaleBMP :: FilePath -> Array DIM2 Double -> IO () #-}
 writeMatrixToGreyscaleBMP fileName arr
  = let arrNorm         = normalisePositive01 arr
-
-       scale :: Double -> Word8
        scale x         = fromIntegral (truncate (x * 255) :: Int)
-
        arrWord8        = A.map scale arrNorm
    in  writeComponentsToBMP fileName arrWord8 arrWord8 arrWord8
                
+               
+-- | Like `writeComponentsToBMP` but take the components as a list.
+--   The list must have 3 arrays, for the red, green blue components
+--   respectively, else `error`.
+writeComponentsListToBMP
+       :: FilePath 
+       -> [Array DIM2 Word8]
+       -> IO ()
+
+writeComponentsListToBMP filePath comps
+       | [red, green, blue]    <- comps
+       = writeComponentsToBMP filePath red green blue
+       
+       | otherwise
+       = error "Data.Array.Repa.IO.BMP.writeComponentsListToBMP: wrong number of components"
+       
 
 -- | Write RGB components to a BMP file.
 --     All arrays must have the same extent, else `error`.
@@ -167,7 +201,7 @@ writeImageToBMP fileName arrImage
        = error "Data.Array.Repa.IO.BMP: lowest order dimension must be 4"
 
        | otherwise
-       = let   bmp     = packRGBA32ToBMP height width 
+       = let   bmp     = packRGBA32ToBMP width height
                        $ A.toByteString arrImage
          in    writeBMP fileName bmp
        
@@ -179,7 +213,7 @@ writeImageToBMP fileName arrImage
 -- | Normalise a matrix to to [0 .. 1], discarding negative values.
 --     If the maximum value is 0 then return the array unchanged.
 normalisePositive01
-       :: (Shape sh, U.Elt a, Fractional a, Ord a)
+       :: (Shape sh, Elt a, Fractional a, Ord a)
        => Array sh a
        -> Array sh a
 
diff --git a/fibon/Repa/_RepaLib/repa-io/Data/Array/Repa/IO/Binary.hs b/fibon/Repa/_RepaLib/repa-io/Data/Array/Repa/IO/Binary.hs
new file mode 100644 (file)
index 0000000..924aa8d
--- /dev/null
@@ -0,0 +1,83 @@
+{-# LANGUAGE FlexibleInstances, ScopedTypeVariables #-}
+{-# OPTIONS -fno-warn-orphans #-}
+
+-- | Reading and writing Repa arrays as binary files.
+module Data.Array.Repa.IO.Binary
+        ( readArrayFromStorableFile
+        , writeArrayToStorableFile)
+where
+import Foreign.Storable
+import Foreign.Ptr
+import Foreign.ForeignPtr
+import Foreign.Marshal.Alloc
+import System.IO
+import Data.Array.Repa                  as R
+import Prelude                          as P
+import Control.Monad
+
+
+-- | Read an array from a file.
+--   Data appears in host byte order.
+--   If the file size does match the provided shape then `error`.
+readArrayFromStorableFile 
+        :: forall a sh
+        .  (Shape sh, Storable a, Elt a)
+        => FilePath 
+        -> sh
+        -> IO (Array sh a)
+
+readArrayFromStorableFile filePath sh
+ = do
+        -- Determine number of bytes per element.
+        let (fake   :: Array sh a)      = R.fromList R.zeroDim []
+        let (bytes1 :: Integer)         = fromIntegral $ sizeOf (fake R.! R.zeroDim)
+
+        -- Determine how many elements the whole file will give us.
+        h :: Handle     <- openBinaryFile filePath ReadMode
+        bytesTotal      <- hFileSize h
+
+        let lenTotal      =  bytesTotal `div` bytes1
+        let bytesExpected =  bytes1 * lenTotal
+        
+        when (bytesTotal /= bytesExpected)
+         $ error $ unlines
+                ["Data.Array.Repa.IO.Binary.readArrayFromStorableFile: not a whole number of elements in file"
+                , "element length = " P.++ show bytes1
+                , "file size      = " P.++ show bytesTotal
+                , "slack space    = " P.++ show (bytesTotal `mod` bytes1) ]
+         
+        let bytesTotal' = fromIntegral bytesTotal
+        buf :: Ptr a    <- mallocBytes bytesTotal' 
+        bytesRead       <- hGetBuf h buf bytesTotal'
+        when (bytesTotal' /= bytesRead)
+         $ error "Data.Array.Repa.IO.Binary.readArrayFromStorableFile: read failed"
+
+        hClose h
+        fptr     <- newForeignPtr finalizerFree buf        
+
+        -- Converting the foreign ptr like this means that the array
+        -- elements are used directly from the buffer, and not copied.
+        let arr  =  R.unsafeFromForeignPtr sh fptr
+
+        return   $  arr `asTypeOf` fake
+
+
+-- | Write an array to a file.
+--   Data appears in host byte order.
+writeArrayToStorableFile
+        :: forall sh a 
+        .  (Shape sh, Storable a, Elt a)
+        => FilePath 
+        -> Array sh a
+        -> IO ()
+
+writeArrayToStorableFile filePath arr
+ = do   let bytes1      = sizeOf (arr R.! R.zeroDim)
+        let bytesTotal  = bytes1 * (R.size $ R.extent arr)
+        
+        buf :: Ptr a    <- mallocBytes bytesTotal
+        R.forceWith (pokeElemOff buf) arr
+        
+        h <- openBinaryFile filePath WriteMode
+        hPutBuf h buf bytesTotal 
+        hClose h
index 6e3291a..31add7b 100644 (file)
@@ -19,15 +19,18 @@ import System.IO
 import Data.List                               as L
 import Data.Array.Repa                         as A
 import Prelude                                 as P
-import qualified "dph-prim-par" Data.Array.Parallel.Unlifted   as U
 
 
+-------------------------------------------------------------------------------
 -- | Read a matrix from a text file.
+--   WARNING: This is implemented fairly naively, just using `Strings` 
+--   under the covers. It will be slow for large data files.
+-- 
+--   It also doesn't do graceful error handling.
+--   If the file has the wrong format you'll get a confusing `error`.
 --
---   WARNING: This doesn't do graceful error handling. If the file has the wrong format
---   you'll get a confusing `error`.
 readMatrixFromTextFile
-       :: (U.Elt a, Num a, Read a)
+       :: (Elt a, Num a, Read a)
        => FilePath
        -> IO (Array DIM2 a)    
 
@@ -47,7 +50,7 @@ readMatrixFromTextFile fileName
 
 -- | Write a matrix as a text file.
 writeMatrixToTextFile 
-       :: (U.Elt a, Show a)
+       :: (Elt a, Show a)
        => FilePath
        -> Array DIM2 a
        -> IO ()
@@ -60,7 +63,7 @@ writeMatrixToTextFile fileName arr
        let Z :. width :. height        
                = extent arr
 
-       hPutStrLn file $ show width ++ " " ++ show height
+       hPutStrLn file $ show width P.++ " " P.++ show height
                
        hWriteValues file $ toList arr
        hClose file
index 94c4661..afbc8b3 100644 (file)
@@ -1,5 +1,9 @@
 module Data.Array.Repa.IO.Timing
-       (time, showTime, prettyTime)
+       ( Time
+       , milliseconds, cpuTime, wallTime
+       , time, minus, plus
+       , showTime
+       , prettyTime)
 where
 import GHC.Exts        (traceEvent)
 import System.CPUTime
@@ -7,6 +11,7 @@ import System.Time
 
 
 -- Time -----------------------------------------------------------------------
+-- | Abstract representation of process time.
 data Time 
        = Time 
        { cpu_time  :: Integer
@@ -17,11 +22,18 @@ zipT :: (Integer -> Integer -> Integer) -> Time -> Time -> Time
 zipT f (Time cpu1 wall1) (Time cpu2 wall2) 
        = Time (f cpu1 cpu2) (f wall1 wall2)
 
+-- | Subtract second time from the first.
 minus :: Time -> Time -> Time
 minus = zipT (-)
 
 
+-- | Add two times.
+plus :: Time -> Time -> Time
+plus  = zipT (+)
+
+
 -- TimeUnit -------------------------------------------------------------------
+-- | Conversion 
 type TimeUnit 
        = Integer -> Integer
 
@@ -50,7 +62,7 @@ showTime t = (show $ wallTime milliseconds t)
           ++ "/"
           ++ (show $ cpuTime  milliseconds t)
 
--- | Pretty print the times.
+-- | Pretty print the times, in milliseconds.
 prettyTime :: Time -> String
 prettyTime t
        = unlines
@@ -58,6 +70,11 @@ prettyTime t
        , "cpuTimeMS       = " ++ (show $ cpuTime  milliseconds t) ]
 
 -- Timing benchmarks ----------------------------------------------------------
+
+-- | Time some IO action.
+--   Make sure to deepseq the result before returning it from the action. If you
+--   don't do this then there's a good chance that you'll just pass a suspension
+--   out of the action, and the computation time will be zero.
 time :: IO a -> IO (a, Time)
 {-# NOINLINE time #-}
 time p = do
index d8630b2..4062558 100644 (file)
@@ -20,14 +20,13 @@ import Prelude                                      as P
 import System.IO
 import Control.Monad
 import Data.Char
-import qualified "dph-prim-par" Data.Array.Parallel.Unlifted   as U
 
 
 -- | Read a vector from a text file.
 --   WARNING: This doesn't do graceful error handling. If the file has the wrong format
 --   you'll get a confusing `error`.
 readVectorFromTextFile
-       :: (U.Elt a, Num a, Read a)
+       :: (Elt a, Num a, Read a)
        => FilePath
        -> IO (Array DIM1 a)    
 
@@ -55,7 +54,7 @@ readInt str
        
 -- | Write a vector as a text file.
 writeVectorToTextFile 
-       :: (U.Elt a, Show a)
+       :: (Elt a, Show a)
        => Array DIM1 a
        -> FilePath
        -> IO ()
index 915df32..93cc7a6 100644 (file)
@@ -1,5 +1,5 @@
 Name:                repa-io
-Version:             1.1.1.0
+Version:             2.2.0.1
 License:             BSD3
 License-file:        LICENSE
 Author:              The DPH Team
@@ -8,25 +8,23 @@ 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
+Homepage:            http://repa.ouroborus.net
+Bug-reports:         repa@ouroborus.net
 Description:
-        NOTE: You must use the GHC head branch > 6.13.20100309 to get decent performance.
         Read and write Repa arrays in various formats.
 
 Synopsis:
         Read and write Repa arrays in various formats.
 
-Tested-with: GHC == 6.13.20100309, GHC == 6.12.1
-
 Library
   Build-Depends: 
         base                 == 4.*,
-        dph-prim-par         == 0.5.*,
-        repa                 == 1.1.*,
-        repa-bytestring      == 1.1.*,
+        repa                 == 2.2.*,
+        repa-bytestring      == 2.2.*,
         bmp                  == 1.1.*,
-        old-time             == 1.0.*
+        old-time             == 1.0.*,
+        vector               == 0.9.*,
+        binary               == 0.5.*
 
   ghc-options:
         -Odph -Wall -fno-warn-missing-signatures
@@ -36,6 +34,7 @@ Library
         Data.Array.Repa.IO.BMP
         Data.Array.Repa.IO.Vector
         Data.Array.Repa.IO.Matrix
+        Data.Array.Repa.IO.Binary
         Data.Array.Repa.IO.Timing
       
   Other-modules:
index 5d50ec8..2b46963 100644 (file)
 {-# LANGUAGE PatternGuards, PackageImports, ScopedTypeVariables, RankNTypes #-}
-{-# LANGUAGE TypeOperators, FlexibleContexts, NoMonomorphismRestriction #-}
+{-# LANGUAGE TypeOperators, FlexibleContexts, NoMonomorphismRestriction, FlexibleInstances, UndecidableInstances #-}
+{-# OPTIONS -fno-warn-orphans #-}
 
 -- | See the repa-examples package for examples.
---   
---   More information is also at <http://trac.haskell.org/repa>
--- 
---   NOTE:     To get decent performance you must use GHC head branch > 6.13.20100309.
 --
---   WARNING:  Most of the functions that operate on indices don't perform bounds checks.
---             Doing these checks would interfere with code optimisation and reduce performance.               
---             Indexing outside arrays, or failing to meet the stated obligations will
---             likely cause heap corruption.
+--   More information at <http://repa.ouroborus.net>.
 --
---   
+--   There is a draft tutorial at <http://www.haskell.org/haskellwiki/Numeric_Haskell:_A_Repa_Tutorial>
+--
+-- @Release Notes:
+--  For 2.2.0.1:
+--   * Added unsafeFromForeignPtr, which helps use foreign source
+--     arrays without intermediate copying.
+--   * Added forceWith and forceWith2, which can be used to force
+--     arrays into foreign result buffers without intermediate copying.
+--
+--  For 2.1.0.1:
+--   * The fold and foldAll functions now run in parallel and require the
+--     starting element to be neutral with respect to the reduction operator.
+--                                   -- thanks to Trevor McDonell
+--   * Added (\/\/) update function.   -- thanks to Trevor McDonell
+--   * Dropped unneeded Elt constraints from traverse functions.
+-- @
 module Data.Array.Repa
        ( module Data.Array.Repa.Shape
        , module Data.Array.Repa.Index
        , module Data.Array.Repa.Slice
-       
-       , Array (..)
 
-        -- * Constructors
-       , fromUArray
+       -- from Data.Array.Repa.Internals.Elt -----------------------
+       , Elt(..)
+
+       -- from Data.Array.Repa.Internals.Base ----------------------
+       , Array(..)
+       , Region(..)
+       , Range(..)
+       , Rect(..)
+       , Generator(..)
+       , deepSeqArray, deepSeqArrays
+       , singleton,    toScalar
+       , extent,       delay
+
+       --
+       , withManifest, withManifest'
+
+       -- * Indexing
+       , (!),  index
+       , (!?), safeIndex
+       , unsafeIndex
+
+       -- * Construction
        , fromFunction
-       , unit
-
-        -- * Projections
-       , extent
-       , delay
-       , toUArray
-       , index, (!:)
-       , toScalar
-
-        -- * Basic Operations
-       , force
-       , deepSeqArray
-       
-        -- * Conversion
+       , fromVector
        , fromList
+       , unsafeFromForeignPtr
+
+       -- from Data.Array.Repa.Interlals.Forcing -------------------
+       -- * Forcing
+       , force,  forceWith
+       , force2, forceWith2
+       , toVector
        , toList
 
-        -- * Index space transformations
+       -- from Data.Array.Repa.Operators.IndexSpace ----------------
+       -- * Index space transformations
        , reshape
-       , append, (+:+)
+       , append, (++)
        , transpose
-       , replicate
+       , extend
        , slice
        , backpermute
        , backpermuteDft
 
-         -- * Structure preserving operations
+       -- from Data.Array.Repa.Operators.Mapping -------------------
+        -- * Structure preserving operations
        , map
        , zipWith
        , (+^), (-^), (*^), (/^)
 
-        -- * Reductions
+        -- from Data.Array.Repa.Operations.Modify -------------------
+        -- * Bulk updates
+        , (//)
+
+       -- from Data.Array.Repa.Operators.Reduction -----------------
+       -- * Reductions
        , fold, foldAll
        , sum,  sumAll
-       
-        -- * Generic traversal
+
+       -- from Data.Array.Repa.Operators.Traverse ------------------
+       -- * Generic Traversal
        , traverse
        , traverse2
        , traverse3
        , traverse4
-               
-        -- * Interleaving
+       , unsafeTraverse
+       , unsafeTraverse2
+       , unsafeTraverse3
+       , unsafeTraverse4
+
+       -- from Data.Array.Repa.Operators.Interleave ----------------
+       -- * Interleaving
        , interleave2
        , interleave3
        , interleave4
-               
-        -- * Testing
-       , arbitrarySmallArray
-       , props_DataArrayRepa)
+
+       -- from Data.Array.Repa.Operators.Select --------------------
+       -- * Selection
+       , select)
+
 where
 import Data.Array.Repa.Index
 import Data.Array.Repa.Slice
 import Data.Array.Repa.Shape
-import Data.Array.Repa.QuickCheck
+import Data.Array.Repa.Internals.Elt
+import Data.Array.Repa.Internals.Base
+import Data.Array.Repa.Internals.Forcing
+import Data.Array.Repa.Operators.Traverse
+import Data.Array.Repa.Operators.IndexSpace
+import Data.Array.Repa.Operators.Interleave
+import Data.Array.Repa.Operators.Mapping
+import Data.Array.Repa.Operators.Modify
+import Data.Array.Repa.Operators.Reduction
+import Data.Array.Repa.Operators.Select
 import qualified Data.Array.Repa.Shape as S
 
-import "dph-prim-par" Data.Array.Parallel.Unlifted             (Elt)
-import qualified "dph-prim-par" Data.Array.Parallel.Unlifted   as U
-import qualified "dph-prim-seq" Data.Array.Parallel.Unlifted.Sequential.Vector as USeq
-
-import Test.QuickCheck
-import Prelude                         hiding (sum, map, zipWith, replicate)   
+import Prelude                         hiding (sum, map, zipWith, (++))
 import qualified Prelude               as P
 
 stage  = "Data.Array.Repa"
-       
--- | Possibly delayed arrays.
-data Array sh a
-       = -- | An array represented as some concrete unboxed data.
-         Manifest !sh !(U.Array a)
-
-          -- | An array represented as a function that computes each element.
-       | Delayed  !sh !(sh -> a)
-
--- Constructors ----------------------------------------------------------------------------------
-
--- | Create a `Manifest` array from an unboxed `U.Array`. 
---     The elements are in row-major order.
-fromUArray
-       :: Shape sh
-       => sh
-       -> U.Array a
-       -> Array sh a
-
-{-# INLINE fromUArray #-}
-fromUArray sh uarr
-       = sh   `S.deepSeq` 
-         uarr `seq`
-         Manifest sh uarr
-
-
--- | Create a `Delayed` array from a function.
-fromFunction 
-       :: Shape sh
-       => sh
-       -> (sh -> a)
-       -> Array sh a
-       
-{-# INLINE fromFunction #-}
-fromFunction sh fnElems
-       = sh `S.deepSeq` Delayed sh fnElems
-
-
--- | Wrap a scalar into a singleton array.
-unit :: Elt a => a -> Array Z a
-{-# INLINE unit #-}
-unit   = Delayed Z . const
-
-
--- Projections ------------------------------------------------------------------------------------
-
--- | Take the extent of an array.
-extent :: Array sh a -> sh
-{-# INLINE extent #-}
-extent arr
- = case arr of
-       Manifest sh _   -> sh
-       Delayed  sh _   -> sh
-
--- | Unpack an array into delayed form.
-delay  :: (Shape sh, Elt a) 
-       => Array sh a 
-       -> (sh, sh -> a)
-
-{-# INLINE delay #-}   
-delay arr
- = case arr of
-       Delayed  sh fn          -> (sh, fn)
-       Manifest sh uarr        -> (sh, \i -> uarr U.!: S.toIndex sh i)
-
-
--- | Convert an array to an unboxed `U.Array`, forcing it if required.
---     The elements come out in row-major order.
-toUArray 
-       :: (Shape sh, Elt a)
-       => Array sh a 
-       -> U.Array a
-{-# INLINE toUArray #-}
-toUArray arr
- = case force arr of
-       Manifest _ uarr -> uarr
-       _               -> error $ stage ++ ".toList: force failed"
-
-
--- | Get an indexed element from an array.
---
---   OBLIGATION: The index must be within the array. 
---
---     @inRange zeroDim (shape arr) ix == True@
---
-index, (!:)
-       :: forall sh a
-       .  (Shape sh, Elt a)
-       => Array sh a
-       -> sh 
-       -> a
-
-{-# INLINE index #-}
-index arr ix
- = case arr of
-       Delayed  _  fn          -> fn ix
-       Manifest sh uarr        -> uarr U.!: (S.toIndex sh ix)
-
-{-# INLINE (!:) #-}
-(!:) arr ix = index arr ix
-
-
--- | Take the scalar value from a singleton array.
-toScalar :: Elt a => Array Z a -> a
-{-# INLINE toScalar #-}
-toScalar arr
- = case arr of
-       Delayed  _ fn           -> fn Z
-       Manifest _ uarr         -> uarr U.!: 0
-
-
--- Basic Operations -------------------------------------------------------------------------------
-
--- | Force an array, so that it becomes `Manifest`.
-force  :: (Shape sh, Elt a)
-       => Array sh a -> Array sh a
-       
-{-# INLINE force #-}
-force arr
- = Manifest sh' uarr'
- where (sh', uarr')
-       = case arr of
-               Manifest sh uarr
-                -> sh `S.deepSeq` uarr `seq` (sh, uarr)
-
-               Delayed sh fn
-                -> let uarr    =  U.map (fn . S.fromIndex sh) 
-                               $! U.enumFromTo (0 :: Int) (S.size sh - 1)
-                    in sh `S.deepSeq` uarr `seq` (sh, uarr)
-
-                               
--- | Ensure an array's structure is fully evaluated.
---     This evaluates the extent and outer constructor, but does not `force` the elements.
-infixr 0 `deepSeqArray`
-deepSeqArray 
-       :: Shape sh
-       => Array sh a 
-       -> b -> b
-
-{-# INLINE deepSeqArray #-}
-deepSeqArray arr x 
- = case arr of
-       Delayed  sh _           -> sh `S.deepSeq` x
-       Manifest sh uarr        -> sh `S.deepSeq` uarr `seq` x
-
-
--- Conversion -------------------------------------------------------------------------------------
--- | Convert a list to an array.
---     The length of the list must be exactly the `size` of the extent given, else `error`.
---
-fromList
-       :: (Shape sh, Elt a)
-       => sh
-       -> [a]
-       -> Array sh a
-       
-{-# INLINE fromList #-}
-fromList sh xx
-       | U.length uarr /= S.size sh
-       = error $ unlines
-               [ stage ++ ".fromList: size of array shape does not match size of list"
-               , "        size of shape = " ++ (show $ S.size sh)      ++ "\n"
-               , "        size of list  = " ++ (show $ U.length uarr)  ++ "\n" ]
-       
-       | otherwise
-       = Manifest sh uarr
-
-       where   uarr    = U.fromList xx
-       
-       
--- | Convert an array to a list.
-toList         :: (Shape sh, Elt a)
-       => Array sh a
-       -> [a]
-
-{-# INLINE toList #-}
-toList arr
- = case force arr of
-       Manifest _ uarr -> U.toList uarr
-       _               -> error $ stage ++ ".toList: force failed"
 
 
 -- Instances --------------------------------------------------------------------------------------
 -- Show
 instance (Shape sh, Elt a, Show a) => Show (Array sh a) where
-       show arr = show $ toList arr
+        show arr =
+          let shape = showShape (extent arr)
+              elems = show      (toList arr)
+          in
+          "Array (" P.++ shape P.++ ") " P.++ elems
+
 
 -- Eq
 instance (Shape sh, Elt a, Eq a) => Eq (Array sh a) where
 
        {-# INLINE (==) #-}
-       (==) arr1  arr2 
-               = toScalar 
-               $ fold (&&) True 
-               $ reshape (Z :. (S.size $ extent arr1)) 
+       (==) arr1  arr2
+               = foldAll (&&) True
+               $ reshape (Z :. (S.size $ extent arr1))
                $ zipWith (==) arr1 arr2
-               
+
        {-# INLINE (/=) #-}
-       (/=) a1 a2 
+       (/=) a1 a2
                = not $ (==) a1 a2
 
 -- Num
@@ -311,501 +174,36 @@ instance (Shape sh, Elt a, Num a) => Num (Array sh a) where
        signum          = map signum
 
        {-# INLINE fromInteger #-}
-       fromInteger n    = Delayed failShape (\_ -> fromInteger n) 
-        where failShape = error $ stage ++ ".fromInteger: Constructed array has no shape."
+       fromInteger n    = fromFunction failShape (\_ -> fromInteger n)
+        where failShape = error $ stage P.++ ".fromInteger: Constructed array has no shape."
 
 
--- Index space transformations --------------------------------------------------------------------
--- | Impose a new shape on the elements of an array.
---     The new extent must be the same size as the original, else `error`.
---
-reshape        :: (Shape sh, Shape sh', Elt a) 
-       => sh'
-       -> Array sh a
-       -> Array sh' a
-
-{-# INLINE reshape #-}
-reshape sh' arr
-       | not $ S.size sh' == S.size (extent arr)
-       = error $ stage ++ ".reshape: reshaped array will not match size of the original"
-
-       | otherwise
-       = case arr of
-               Manifest _  uarr -> Manifest sh' uarr
-               Delayed  sh f    -> Delayed sh' (f . fromIndex sh . toIndex sh')
-       
-
--- | Append two arrays.
---
---   OBLIGATION: The higher dimensions of both arrays must have the same extent.
---
---   @tail (listOfShape (shape arr1)) == tail (listOfShape (shape arr2))@
---
-append, (+:+)  
+-- | Force an array before passing it to a function.
+withManifest
        :: (Shape sh, Elt a)
-       => Array (sh :. Int) a
-       -> Array (sh :. Int) a
-       -> Array (sh :. Int) a
-
-{-# INLINE append #-}
-append arr1 arr2 
- = traverse2 arr1 arr2 fnExtent fnElem
- where
-       (_ :. n)        = extent arr1
-
-       fnExtent (sh :. i) (_  :. j) 
-               = sh :. (i + j)
-
-       fnElem f1 f2 (sh :. i)
-               | i < n         = f1 (sh :. i)
-               | otherwise     = f2 (sh :. (i - n))
-
-{-# INLINE (+:+) #-}
-(+:+) arr1 arr2 = append arr1 arr2
-
-
--- | Transpose the lowest two dimensions of an array. 
---     Transposing an array twice yields the original.
-transpose 
-       :: (Shape sh, Elt a) 
-       => Array (sh :. Int :. Int) a
-       -> Array (sh :. Int :. Int) a
-
-{-# INLINE transpose #-}
-transpose arr 
- = traverse arr
-       (\(sh :. m :. n)        -> (sh :. n :.m))
-       (\f -> \(sh :. i :. j)  -> f (sh :. j :. i))
-
-
--- | Replicate an array, according to a given slice specification.
-replicate
-       :: ( Slice sl
-          , Shape (FullShape sl)
-          , Shape (SliceShape sl)
-          , Elt e)
-       => sl
-       -> Array (SliceShape sl) e
-       -> Array (FullShape sl) e
-
-{-# INLINE replicate #-}
-replicate sl arr
-       = backpermute 
-               (fullOfSlice sl (extent arr)) 
-               (sliceOfFull sl)
-               arr
-
--- | Take a slice from an array, according to a given specification.
-slice  :: ( Slice sl
-          , Shape (FullShape sl)
-          , Shape (SliceShape sl)
-          , Elt e)
-       => Array (FullShape sl) e
-       -> sl
-       -> Array (SliceShape sl) e
-
-{-# INLINE slice #-}
-slice arr sl
-       = backpermute 
-               (sliceOfFull sl (extent arr))
-               (fullOfSlice sl)
-               arr
-
-
--- | Backwards permutation of an array's elements.
---     The result array has the same extent as the original.
-backpermute
-       :: forall sh sh' a
-       .  (Shape sh, Shape sh', Elt a) 
-       => sh'                          -- ^ Extent of result array.
-       -> (sh' -> sh)                  -- ^ Function mapping each index in the result array
-                                       --      to an index of the source array.
-       -> Array sh a                   -- ^ Source array.
-       -> Array sh' a
-
-{-# INLINE backpermute #-}
-backpermute newExtent perm arr
-       = traverse arr (const newExtent) (. perm) 
-       
-
--- | Default backwards permutation of an array's elements.
---     If the function returns `Nothing` then the value at that index is taken
---     from the default array (@arrDft@)
-backpermuteDft
-       :: forall sh sh' a
-       .  (Shape sh, Shape sh', Elt a) 
-       => Array sh' a                  -- ^ Default values (@arrDft@)
-       -> (sh' -> Maybe sh)            -- ^ Function mapping each index in the result array
-                                       --      to an index in the source array.
-       -> Array sh  a                  -- ^ Source array.
-       -> Array sh' a
-
-{-# INLINE backpermuteDft #-}
-backpermuteDft arrDft fnIndex arrSrc
-       = Delayed (extent arrDft) fnElem
-       where   fnElem ix       
-                = case fnIndex ix of
-                       Just ix'        -> arrSrc !: ix'
-                       Nothing         -> arrDft !: ix
-                               
-
--- Structure Preserving Operations ----------------------------------------------------------------
--- | Apply a worker function to each element of an array, 
---     yielding a new array with the same extent.
-map    :: (Shape sh, Elt a, Elt b) 
-       => (a -> b)
-       -> Array sh a
-       -> Array sh b
-
-{-# INLINE map #-}
-map f arr
-       = traverse arr id (f .)
-       
-
--- | Combine two arrays, element-wise, with a binary operator.
---     If the extent of the two array arguments differ, 
---     then the resulting array's extent is their intersection.
-zipWith :: (Shape sh, Elt a, Elt b, Elt c) 
-       => (a -> b -> c) 
-       -> Array sh a
-       -> Array sh b
-       -> Array sh c
-
-{-# INLINE zipWith #-}
-zipWith f arr1 arr2
-       = arr1 `deepSeqArray` 
-         arr2 `deepSeqArray`
-         Delayed       (S.intersectDim (extent arr1) (extent arr2))
-                       (\ix -> f (arr1 !: ix) (arr2 !: ix))
-
-{-# INLINE (+^) #-}
-(+^)   = zipWith (+)
-
-{-# INLINE (-^) #-}
-(-^)   = zipWith (-)
-
-{-# INLINE (*^) #-}
-(*^)   = zipWith (*)
-
-{-# INLINE (/^) #-}
-(/^)   = zipWith (/)
-
-
--- Reductions -------------------------------------------------------------------------------------
-
--- IMPORTANT: 
---     These reductions use the sequential version of foldU, mapU and enumFromToU.
---     If we use parallel versions then we'll end up with nested parallelism
---     and the gang will abort at runtime.
-
--- | Fold the innermost dimension of an array.
---     Combine this with `transpose` to fold any other dimension.
-fold   :: (Shape sh, Elt a)
-       => (a -> a -> a)
-       -> a 
-       -> Array (sh :. Int) a
-       -> Array sh a
-
-{-# INLINE fold #-}
-fold f x arr
- = x `seq` arr `deepSeqArray` 
-   let sh' :. n        = extent arr
-       elemFn i        = USeq.fold f x
-                       $ USeq.map
-                               (\ix -> arr !: (i :. ix)) 
-                               (USeq.enumFromTo 0 (n - 1))
-   in  Delayed sh' elemFn
-
-
--- | Fold all the elements of an array.
-foldAll :: (Shape sh, Elt a)
-       => (a -> a -> a)
-       -> a
-       -> Array sh a
-       -> a
-       
-{-# INLINE foldAll #-}
-foldAll f x arr
-       = USeq.fold f x
-       $ USeq.map ((arr !:) . (S.fromIndex (extent arr)))
-       $ USeq.enumFromTo
-               0
-               ((S.size $ extent arr) - 1)
-
-
-
--- | Sum the innermost dimension of an array.
-sum    :: (Shape sh, Elt a, Num a)
-       => Array (sh :. Int) a
-       -> Array sh a
-
-{-# INLINE sum #-}
-sum arr        = fold (+) 0 arr
-
-
--- | Sum all the elements of an array.
-sumAll :: (Shape sh, Elt a, Num a)
-       => Array sh a
-       -> a
-
-{-# INLINE sumAll #-}
-sumAll arr
-       = USeq.fold (+) 0
-       $ USeq.map ((arr !:) . (S.fromIndex (extent arr)))
-       $ USeq.enumFromTo
-               0
-               ((S.size $ extent arr) - 1)
-
-
--- Generic Traversal -----------------------------------------------------------------------------
--- | Unstructured traversal.
-traverse
-       :: forall sh sh' a b
-       .  (Shape sh, Shape sh', Elt a)
-       => Array sh a                           -- ^ Source array.
-       -> (sh  -> sh')                         -- ^ Function to produce the extent of the result.
-       -> ((sh -> a) -> sh' -> b)              -- ^ Function to produce elements of the result. 
-                                               --   It is passed a lookup function to get elements of the source.
-       -> Array sh' b
-       
-{-# INLINE traverse #-}
-traverse arr transExtent newElem
-       = arr `deepSeqArray` 
-          Delayed (transExtent sh) (newElem f)
-       where   (sh, f) = delay arr
-
-
--- | Unstructured traversal over two arrays at once.
-traverse2
-       :: forall sh sh' sh'' a b c
-       .  ( Shape sh, Shape sh', Shape sh''
-          , Elt a,    Elt b,     Elt c)
-        => Array sh a                          -- ^ First source array.
-       -> Array sh' b                          -- ^ Second source array.
-        -> (sh -> sh' -> sh'')                 -- ^ Function to produce the extent of the result.
-        -> ((sh -> a) -> (sh' -> b)            
-                      -> (sh'' -> c))          -- ^ Function to produce elements of the result.
-                                               --   It is passed lookup functions to get elements of the 
-                                               --   source arrays.
-        -> Array sh'' c 
-
-{-# INLINE traverse2 #-}
-traverse2 arrA arrB transExtent newElem
-       = arrA `deepSeqArray` arrB `deepSeqArray`
-         Delayed 
-               (transExtent (extent arrA) (extent arrB)) 
-               (newElem     ((!:) arrA) ((!:) arrB))
-
-
--- | Unstructured traversal over three arrays at once.
-traverse3
-       :: forall sh1 sh2 sh3 sh4
-                 a   b   c   d 
-       .  ( Shape sh1, Shape sh2, Shape sh3, Shape sh4
-          , Elt a,     Elt b,     Elt c,     Elt d)
-        => Array sh1 a                 
-       -> Array sh2 b                  
-       -> Array sh3 c                  
-        -> (sh1 -> sh2 -> sh3 -> sh4)  
-        -> (  (sh1 -> a) -> (sh2 -> b) 
-           -> (sh3 -> c)
-           ->  sh4 -> d )              
-        -> Array sh4 d
-
-{-# INLINE traverse3 #-}
-traverse3 arrA arrB arrC transExtent newElem
-       = arrA `deepSeqArray` arrB `deepSeqArray` arrC `deepSeqArray`
-         Delayed 
-               (transExtent (extent arrA) (extent arrB) (extent arrC)) 
-               (newElem     (arrA !:) (arrB !:) (arrC !:))
-
-
--- | Unstructured traversal over four arrays at once.
-traverse4
-       :: forall sh1 sh2 sh3 sh4 sh5 
-                 a   b   c   d   e
-       .  ( Shape sh1, Shape sh2, Shape sh3, Shape sh4, Shape sh5
-          , Elt a,     Elt b,     Elt c,     Elt d,     Elt e)
-        => Array sh1 a                         
-       -> Array sh2 b                  
-       -> Array sh3 c                  
-       -> Array sh4 d                          
-        -> (sh1 -> sh2 -> sh3 -> sh4 -> sh5 )  
-        -> (  (sh1 -> a) -> (sh2 -> b) 
-           -> (sh3 -> c) -> (sh4 -> d)
-           ->  sh5 -> e )              
-        -> Array sh5 e 
-
-{-# INLINE traverse4 #-}
-traverse4 arrA arrB arrC arrD transExtent newElem
-       = arrA `deepSeqArray` arrB `deepSeqArray` arrC `deepSeqArray` arrD `deepSeqArray` 
-         Delayed 
-               (transExtent (extent arrA) (extent arrB) (extent arrC) (extent arrD)) 
-               (newElem     (arrA !:) (arrB !:) (arrC !:) (arrD !:))
-
-
--- Interleaving -----------------------------------------------------------------------------------
--- | Interleave the elments of two arrays. 
---   All the input arrays must have the same extent, else `error`.
---   The lowest dimenion of the result array is twice the size of the inputs.
---
--- @
---  interleave2 a1 a2   b1 b2  =>  a1 b1 a2 b2
---              a3 a4   b3 b4      a3 b3 a4 b4
--- @
---
-interleave2
-       :: (Shape sh, Elt a)
-       => Array (sh :. Int) a
-       -> Array (sh :. Int) a
-       -> Array (sh :. Int) a
-       
-{-# INLINE interleave2 #-}
-interleave2 arr1 arr2
- = arr1 `deepSeqArray` arr2 `deepSeqArray`
-   traverse2 arr1 arr2 shapeFn elemFn
- where
-       shapeFn dim1 dim2
-        | dim1 == dim2
-        , sh :. len    <- dim1
-        = sh :. (len * 2)
-       
-        | otherwise
-        = error "Data.Array.Repa.interleave2: arrays must have same extent"
-               
-       elemFn get1 get2 (sh :. ix)
-        = case ix `mod` 3 of
-               0       -> get1 (sh :. ix `div` 2)
-               1       -> get2 (sh :. ix `div` 2)
-               _       -> error "Data.Array.Repa.interleave2: this never happens :-P"
-
-
--- | Interleave the elments of three arrays. 
-interleave3
-       :: (Shape sh, Elt a)
-       => Array (sh :. Int) a
-       -> Array (sh :. Int) a
-       -> Array (sh :. Int) a
-       -> Array (sh :. Int) a
-       
-{-# INLINE interleave3 #-}
-interleave3 arr1 arr2 arr3
- = arr1 `deepSeqArray` arr2 `deepSeqArray` arr3 `deepSeqArray`
-   traverse3 arr1 arr2 arr3 shapeFn elemFn
- where
-       shapeFn dim1 dim2 dim3
-        | dim1 == dim2
-        , dim1 == dim3
-        , sh :. len    <- dim1
-        = sh :. (len * 3)
-       
-        | otherwise
-        = error "Data.Array.Repa.interleave3: arrays must have same extent"
-               
-       elemFn get1 get2 get3 (sh :. ix)
-        = case ix `mod` 3 of
-               0       -> get1 (sh :. ix `div` 3)
-               1       -> get2 (sh :. ix `div` 3)
-               2       -> get3 (sh :. ix `div` 3)
-               _       -> error "Data.Array.Repa.interleave3: this never happens :-P"
-
-
--- | Interleave the elments of four arrays. 
-interleave4
+       => (Array sh a -> b) -> Array sh a -> b
+
+{-# INLINE withManifest #-}
+withManifest f arr
+ = case arr of
+       Array sh [Region RangeAll (GenManifest vec)]
+         -> vec `seq` f (Array sh [Region RangeAll (GenManifest vec)])
+
+       _ -> f (force arr)
+
+
+-- | Force an array before passing it to a function.
+withManifest'
        :: (Shape sh, Elt a)
-       => Array (sh :. Int) a
-       -> Array (sh :. Int) a
-       -> Array (sh :. Int) a
-       -> Array (sh :. Int) a
-       -> Array (sh :. Int) a
-       
-{-# INLINE interleave4 #-}
-interleave4 arr1 arr2 arr3 arr4
- = arr1 `deepSeqArray` arr2 `deepSeqArray` arr3 `deepSeqArray` arr4 `deepSeqArray`
-   traverse4 arr1 arr2 arr3 arr4 shapeFn elemFn
- where
-       shapeFn dim1 dim2 dim3 dim4
-        | dim1 == dim2
-        , dim1 == dim3
-        , dim1 == dim4
-        , sh :. len    <- dim1
-        = sh :. (len * 4)
-       
-        | otherwise
-        = error "Data.Array.Repa.interleave4: arrays must have same extent"
-               
-       elemFn get1 get2 get3 get4 (sh :. ix)
-        = case ix `mod` 4 of
-               0       -> get1 (sh :. ix `div` 4)
-               1       -> get2 (sh :. ix `div` 4)
-               2       -> get3 (sh :. ix `div` 4)
-               3       -> get4 (sh :. ix `div` 4)
-               _       -> error "Data.Array.Repa.interleave4: this never happens :-P"
-
-
--- Arbitrary --------------------------------------------------------------------------------------
--- | Create an arbitrary small array, restricting the size of each of the dimensions to some value.
-arbitrarySmallArray 
-       :: (Shape sh, Elt a, Arbitrary sh, Arbitrary a)
-       => Int
-       -> Gen (Array (sh :. Int) a)
-
-arbitrarySmallArray maxDim
- = do  sh      <- arbitrarySmallShape maxDim
-       xx      <- arbitraryListOfLength (S.size sh)
-       return  $ fromList sh xx
-
-
-
--- Properties -------------------------------------------------------------------------------------
-
--- | QuickCheck properties for this module and its children.
-props_DataArrayRepa :: [(String, Property)]
-props_DataArrayRepa
- =  props_DataArrayRepaIndex
- ++ [(stage ++ "." ++ name, test) | (name, test)
-    <- [ ("id_force/DIM5",                     property prop_id_force_DIM5)
-       , ("id_toScalarUnit",                   property prop_id_toScalarUnit)
-       , ("id_toListFromList/DIM3",            property prop_id_toListFromList_DIM3) 
-       , ("id_transpose/DIM4",                 property prop_id_transpose_DIM4)
-       , ("reshapeTransposeSize/DIM3",         property prop_reshapeTranspose_DIM3)
-       , ("appendIsAppend/DIM3",               property prop_appendIsAppend_DIM3)
-       , ("sumAllIsSum/DIM3",                  property prop_sumAllIsSum_DIM3) ]]
-       
-
--- The Eq instance uses fold and zipWith.
-prop_id_force_DIM5
- =     forAll (arbitrarySmallArray 10)                 $ \(arr :: Array DIM5 Int) ->
-       arr == force arr
-       
-prop_id_toScalarUnit (x :: Int)
- =     toScalar (unit x) == x
-
--- Conversions ------------------------
-prop_id_toListFromList_DIM3
- =     forAll (arbitrarySmallShape 10)                 $ \(sh :: DIM3) ->
-       forAll (arbitraryListOfLength (S.size sh))      $ \(xx :: [Int]) ->
-       toList (fromList sh xx) == xx
-
--- Index Space Transforms -------------
-prop_id_transpose_DIM4
- =     forAll (arbitrarySmallArray 20)                 $ \(arr :: Array DIM3 Int) ->
-       transpose (transpose arr) == arr
-
--- A reshaped array has the same size and sum as the original
-prop_reshapeTranspose_DIM3
- =     forAll (arbitrarySmallArray 20)                 $ \(arr :: Array DIM3 Int) ->
-   let arr'    = transpose arr
-       sh'     = extent arr'
-   in  (S.size $ extent arr) == S.size (extent (reshape sh' arr))
-     && (sumAll arr          == sumAll arr')
-
-prop_appendIsAppend_DIM3
- =     forAll (arbitrarySmallArray 20)                 $ \(arr1 :: Array DIM3 Int) ->
-       sumAll (append arr1 arr1) == (2 * sumAll arr1)
-
--- Reductions --------------------------
-prop_sumAllIsSum_DIM3
- =     forAll (arbitrarySmallShape 100)                $ \(sh :: DIM2) ->
-       forAll (arbitraryListOfLength (S.size sh))      $ \(xx :: [Int]) -> 
-       sumAll (fromList sh xx) == P.sum xx
+       => Array sh a -> (Array sh a -> b) -> b
+
+{-# INLINE withManifest' #-}
+withManifest' arr f
+ = case arr of
+       Array sh [Region RangeAll (GenManifest vec)]
+        -> vec `seq` f (Array sh [Region RangeAll (GenManifest vec)])
+
+       _ -> f (force arr)
+
+
+
diff --git a/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Arbitrary.hs b/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Arbitrary.hs
new file mode 100644 (file)
index 0000000..4f59dfb
--- /dev/null
@@ -0,0 +1,99 @@
+{-# LANGUAGE TypeOperators, FlexibleInstances #-}
+{-# OPTIONS -fno-warn-orphans #-}
+
+-- Utils to help with testing. Not exported.
+module Data.Array.Repa.Arbitrary
+       ( arbitraryShape
+       , arbitrarySmallShape
+       , arbitraryListOfLength
+       , arbitrarySmallArray)
+where
+import Data.Array.Repa.Internals.Elt
+import Data.Array.Repa.Internals.Base
+import Data.Array.Repa.Index
+import Data.Array.Repa.Shape   as S
+import Control.Monad
+import Test.QuickCheck
+
+
+-- Arbitrary --------------------------------------------------------------------------------------
+instance Arbitrary Z where
+       arbitrary       = return Z
+
+-- | Generate an arbitrary index, which may have 0's for some components.
+instance (Shape sh, Arbitrary sh) => Arbitrary (sh :. Int)  where
+       arbitrary
+        = do   sh1             <- arbitrary
+               let sh1Unit     = if size sh1 == 0 then unitDim else sh1
+
+               -- Make sure not to create an index so big that we get
+               --      integer overflow when converting it to the linear form.
+               n               <- liftM abs $ arbitrary
+               let nMax        = maxBound `div` (size sh1Unit)
+               let nMaxed      = n `mod` nMax
+
+               return  $ sh1 :. nMaxed
+
+-- | Generate an aribrary shape that does not have 0's for any component.
+arbitraryShape
+       :: (Shape sh, Arbitrary sh)
+       => Gen (sh :. Int)
+
+arbitraryShape
+ = do  sh1             <- arbitrary
+       let sh1Unit     = if size sh1 == 0 then unitDim else sh1
+
+       -- Make sure not to create an index so big that we get
+       --      integer overflow when converting it to the linear form.
+       n               <- liftM abs $ arbitrary
+       let nMax        = maxBound `div` size sh1Unit
+       let nMaxed      = n `mod` nMax
+       let nClamped    = if nMaxed == 0 then 1 else nMaxed
+
+       return $ sh1Unit :. nClamped
+
+
+-- | Generate an arbitrary shape where each dimension is more than zero,
+--     but less than a specific value.
+arbitrarySmallShape
+       :: (Shape sh, Arbitrary sh)
+       => Int
+       -> Gen (sh :. Int)
+
+arbitrarySmallShape maxDim
+ = do  sh              <- arbitraryShape
+       let dims        = listOfShape sh
+
+       let clamp x
+               = case x `mod` maxDim of
+                       0       -> 1
+                       n       -> n
+
+       return  $ if True
+                       then shapeOfList $ map clamp dims
+                       else sh
+
+
+arbitraryListOfLength
+       :: Arbitrary a
+       => Int -> Gen [a]
+
+arbitraryListOfLength n
+       | n == 0                = return []
+       | otherwise
+       = do    i       <- arbitrary
+               rest    <- arbitraryListOfLength (n - 1)
+               return  $ i : rest
+
+-- | Create an arbitrary small array, restricting the size of each of the
+--   dimensions to some value.
+arbitrarySmallArray
+       :: (Shape sh, Elt a, Arbitrary sh, Arbitrary a)
+       => Int
+       -> Gen (Array (sh :. Int) a)
+
+arbitrarySmallArray maxDim
+ = do  sh      <- arbitrarySmallShape maxDim
+       xx      <- arbitraryListOfLength (S.size sh)
+       return  $ fromList sh xx
+
index 170b969..149319c 100644 (file)
@@ -2,7 +2,7 @@
 
 -- | Index types.
 module Data.Array.Repa.Index
-       ( 
+       (
        -- * Index types
          Z     (..)
        , (:.)  (..)
@@ -13,16 +13,9 @@ module Data.Array.Repa.Index
        , DIM2
        , DIM3
        , DIM4
-       , DIM5 
-       
-       -- * Testing
-       , arbitraryShape
-       , arbitrarySmallShape
-       , props_DataArrayRepaIndex)
+       , DIM5)
 where
 import Data.Array.Repa.Shape
-import Test.QuickCheck
-import Control.Monad
 import GHC.Base                (quotInt, remInt)
 
 stage  = "Data.Array.Repa.Index"
@@ -60,6 +53,9 @@ instance Shape Z where
        {-# INLINE intersectDim #-}
        intersectDim _ _        = Z
 
+       {-# INLINE addDim #-}
+       addDim _ _              = Z
+
        {-# INLINE size #-}
        size _                  = 1
 
@@ -74,8 +70,8 @@ instance Shape Z where
        fromIndex _ _           = Z
 
 
-       {-# INLINE inRange #-}
-       inRange Z Z Z           = True
+       {-# INLINE inShapeRange #-}
+       inShapeRange Z Z Z      = True
 
        listOfShape _           = []
        shapeOfList []          = Z
@@ -84,7 +80,7 @@ instance Shape Z where
        {-# INLINE deepSeq #-}
        deepSeq Z x             = x
 
-       
+
 instance Shape sh => Shape (sh :. Int) where
        {-# INLINE rank #-}
        rank   (sh  :. _)
@@ -97,9 +93,13 @@ instance Shape sh => Shape (sh :. Int) where
        unitDim = unitDim :. 1
 
        {-# INLINE intersectDim #-}
-       intersectDim (sh1 :. n1) (sh2 :. n2) 
+       intersectDim (sh1 :. n1) (sh2 :. n2)
                = (intersectDim sh1 sh2 :. (min n1 n2))
 
+       {-# INLINE addDim #-}
+       addDim (sh1 :. n1) (sh2 :. n2)
+               = addDim sh1 sh2 :. (n1 + n2)
+
        {-# INLINE size #-}
        size  (sh1 :. n)
                = size sh1 * n
@@ -108,16 +108,16 @@ instance Shape sh => Shape (sh :. Int) where
        sizeIsValid (sh1 :. n)
                | size sh1 > 0
                = n <= maxBound `div` size sh1
-               
+
                | otherwise
                = False
-               
+
        {-# INLINE toIndex #-}
-       toIndex (sh1 :. sh2) (sh1' :. sh2') 
+       toIndex (sh1 :. sh2) (sh1' :. sh2')
                = toIndex sh1 sh1' * sh2 + sh2'
 
        {-# INLINE fromIndex #-}
-       fromIndex (ds :. d) n 
+       fromIndex (ds :. d) n
                = fromIndex ds (n `quotInt` d) :. r
                where
                -- If we assume that the index is in range, there is no point
@@ -127,9 +127,9 @@ instance Shape sh => Shape (sh :. Int) where
                r       | rank ds == 0  = n
                        | otherwise     = n `remInt` d
 
-       {-# INLINE inRange #-}
-       inRange (zs :. z) (sh1 :. n1) (sh2 :. n2) 
-               = (n2 >= z) && (n2 < n1) && (inRange zs sh1 sh2)
+       {-# INLINE inShapeRange #-}
+       inShapeRange (zs :. z) (sh1 :. n1) (sh2 :. n2)
+               = (n2 >= z) && (n2 < n1) && (inShapeRange zs sh1 sh2)
 
 
                listOfShape (sh :. n)
@@ -138,98 +138,8 @@ instance Shape sh => Shape (sh :. Int) where
        shapeOfList xx
         = case xx of
                []      -> error $ stage ++ ".toList: empty list when converting to  (_ :. Int)"
-               x:xs    -> shapeOfList xs :. x                  
+               x:xs    -> shapeOfList xs :. x
 
-       {-# INLINE deepSeq #-} 
+       {-# INLINE deepSeq #-}
        deepSeq (sh :. n) x = deepSeq sh (n `seq` x)
 
-
-
-
--- Arbitrary --------------------------------------------------------------------------------------
-instance Arbitrary Z where
-       arbitrary       = return Z
-
--- | Generate an arbitrary index, which may have 0's for some components.
-instance (Shape sh, Arbitrary sh) => Arbitrary (sh :. Int)  where
-       arbitrary 
-        = do   sh1             <- arbitrary
-               let sh1Unit     = if size sh1 == 0 then unitDim else sh1
-               
-               -- Make sure not to create an index so big that we get
-               --      integer overflow when converting it to the linear form.
-               n               <- liftM abs $ arbitrary
-               let nMax        = maxBound `div` (size sh1Unit)
-               let nMaxed      = n `mod` nMax
-               
-               return  $ sh1 :. nMaxed 
-
--- | Generate an aribrary shape that does not have 0's for any component.
-arbitraryShape 
-       :: (Shape sh, Arbitrary sh) 
-       => Gen (sh :. Int)
-
-arbitraryShape 
- = do  sh1             <- arbitrary
-       let sh1Unit     = if size sh1 == 0 then unitDim else sh1
-
-       -- Make sure not to create an index so big that we get
-       --      integer overflow when converting it to the linear form.
-       n               <- liftM abs $ arbitrary
-       let nMax        = maxBound `div` size sh1Unit
-       let nMaxed      = n `mod` nMax
-       let nClamped    = if nMaxed == 0 then 1 else nMaxed
-       
-       return $ sh1Unit :. nClamped
-       
-       
--- | Generate an arbitrary shape where each dimension is more than zero, 
---     but less than a specific value.
-arbitrarySmallShape 
-       :: (Shape sh, Arbitrary sh)
-       => Int
-       -> Gen (sh :. Int)
-
-arbitrarySmallShape maxDim
- = do  sh              <- arbitraryShape
-       let dims        = listOfShape sh
-
-       let clamp x
-               = case x `mod` maxDim of
-                       0       -> 1
-                       n       -> n
-                                               
-       return  $ if True 
-                       then shapeOfList $ map clamp dims
-                       else sh
-
-
-genInShape2 :: DIM2 -> Gen DIM2
-genInShape2 (Z :. yMax :. xMax)
- = do  y       <- liftM (`mod` yMax) $ arbitrary
-       x       <- liftM (`mod` xMax) $ arbitrary
-       return  $ Z :. y :. x
-
-
--- Properties -------------------------------------------------------------------------------------
--- | QuickCheck properties for this module.
-props_DataArrayRepaIndex :: [(String, Property)]
-props_DataArrayRepaIndex
-  = [(stage ++ "." ++ name, test) | (name, test)
-     <-        [ ("toIndexFromIndex/DIM1",     property prop_toIndexFromIndex_DIM1) 
-       , ("toIndexFromIndex/DIM2",     property prop_toIndexFromIndex_DIM2) ]]
-
-prop_toIndexFromIndex_DIM1 sh ix
-       =   (sizeIsValid sh)
-       ==> (inShape sh ix)
-       ==> fromIndex sh (toIndex sh ix) == ix
-       where   _types  = ( sh :: DIM1
-                         , ix :: DIM1)
-
-prop_toIndexFromIndex_DIM2
- =     forAll arbitraryShape   $ \(sh :: DIM2) ->
-       forAll (genInShape2 sh) $ \(ix :: DIM2) ->
-       fromIndex sh (toIndex sh ix) == ix
-
-
-
diff --git a/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/Base.hs b/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/Base.hs
new file mode 100644 (file)
index 0000000..5f160ef
--- /dev/null
@@ -0,0 +1,410 @@
+{-# LANGUAGE ExplicitForAll, TypeOperators, FlexibleInstances, UndecidableInstances, BangPatterns,
+             ExistentialQuantification #-}
+module Data.Array.Repa.Internals.Base
+       ( Array (..)
+       , Region(..)
+       , Range (..)
+       , Rect  (..)
+       , Generator(..)
+       , deepSeqArray, deepSeqArrays
+
+       , singleton, toScalar
+       , extent,    delay
+
+       -- * Predicates
+       , inRange
+
+       -- * Indexing
+       , (!),  index
+       , (!?), safeIndex
+       , unsafeIndex
+
+       -- * Construction
+       , fromFunction
+       , fromVector
+       , fromList
+       , unsafeFromForeignPtr)
+where
+import Data.Array.Repa.Index
+import Data.Array.Repa.Internals.Elt
+import Data.Array.Repa.Shape                   as S
+import qualified Data.Vector.Unboxed           as V
+import Data.Vector.Unboxed                     (Vector)
+import Foreign.ForeignPtr
+import Foreign.Storable
+import System.IO.Unsafe
+
+stage  = "Data.Array.Repa.Array"
+
+-- Array ----------------------------------------------------------------------
+-- | Repa arrays.
+data Array sh a
+       = Array
+       { -- | The entire extent of the array.
+         arrayExtent           :: sh
+
+         -- | Arrays can be partitioned into several regions.
+       , arrayRegions          :: [Region sh a] }
+
+
+-- | Defines the values in a region of the array.
+data Region sh a
+       = Region
+       { -- | The range of elements this region applies to.
+         regionRange           :: Range sh
+
+         -- | How to compute the array elements in this region.
+       , regionGenerator       :: Generator sh a }
+
+
+-- | Represents a range of elements in the array.
+data Range sh
+         -- | Covers the entire array.
+       = RangeAll
+
+         -- | The union of a possibly disjoint set of rectangles.
+       | RangeRects
+       { rangeMatch    :: sh -> Bool
+       , rangeRects    :: [Rect sh] }
+
+
+-- | A rectangle\/cube of arbitrary dimension.
+--   The indices are of the minimum and maximim elements to fill.
+data Rect sh
+       = Rect sh sh
+
+-- | Generates array elements for a particular region in the array.
+data Generator sh a
+       -- | Elements are already computed and sitting in this vector.
+       = GenManifest (Vector a)
+       --   NOTE: Don't make the vector field strict. If you do then deepSeqing arrays
+       --         outside of loops won't cause the unboxings to be floated out.
+
+       -- | Elements can be computed using these cursor functions.
+       | forall cursor
+       . GenCursor
+       { -- | Make a cursor to a particular element.
+         genMakeCursor         :: sh -> cursor
+
+         -- | Shift the cursor by an offset, to get to another element.
+       , genShiftCursor        :: sh -> cursor -> cursor
+
+         -- | Load\/compute the element at the given cursor.
+       , genLoadElem           :: cursor -> a }
+
+
+-- DeepSeqs -------------------------------------------------------------------
+-- | Ensure the structure for an array is fully evaluated.
+--   As we are in a lazy language, applying the @force@ function to a delayed array doesn't
+--   actually compute it at that point. Rather, Haskell builds a suspension representing the
+--   appliction of the @force@ function to that array. Use @deepSeqArray@ to ensure the array
+--   is actually computed at a particular point in the program.
+infixr 0 `deepSeqArray`
+deepSeqArray :: Shape sh => Array sh a -> b -> b
+{-# INLINE deepSeqArray #-}
+deepSeqArray (Array ex rgns) x
+       = ex `S.deepSeq` rgns `deepSeqRegions` x
+
+-- | Like `deepSeqArray` but seqs all the arrays in a list.
+--   This is specialised up to lists of 4 arrays. Using more in the list will break fusion.
+infixr 0 `deepSeqArrays`
+deepSeqArrays :: Shape sh => [Array sh a] -> b -> b
+{-# INLINE deepSeqArrays #-}
+deepSeqArrays as y
+ = case as of
+       []              -> y
+       [a]             -> a  `deepSeqArray` y
+       [a1, a2]        -> a1 `deepSeqArray` a2 `deepSeqArray` y
+       [a1, a2, a3]    -> a1 `deepSeqArray` a2 `deepSeqArray` a3 `deepSeqArray` y
+       [a1, a2, a3, a4]-> a1 `deepSeqArray` a2 `deepSeqArray` a3 `deepSeqArray` a4 `deepSeqArray` y
+       _               -> deepSeqArrays' as y
+
+deepSeqArrays' as' y
+ = case as' of
+       []      -> y
+       x : xs  -> x `deepSeqArray` xs `deepSeqArrays` y
+
+-- | Ensure the structure for a region is fully evaluated.
+infixr 0 `deepSeqRegion`
+deepSeqRegion :: Shape sh => Region sh a -> b -> b
+{-# INLINE deepSeqRegion #-}
+deepSeqRegion (Region range gen) x
+       = range `deepSeqRange` gen `deepSeqGen` x
+
+
+-- | Ensure the structure for some regions are fully evaluated.
+infixr 0 `deepSeqRegions`
+deepSeqRegions :: Shape sh => [Region sh a] -> b -> b
+{-# INLINE deepSeqRegions #-}
+deepSeqRegions rs y
+ = case rs of
+       []              -> y
+       [r]             -> r  `deepSeqRegion`  y
+       [r1, r2]        -> r1 `deepSeqRegion` r2 `deepSeqRegion` y
+       rs'             -> deepSeqRegions' rs' y
+
+deepSeqRegions' rs' y
+ = case rs' of
+       []      -> y
+       x : xs  -> x `deepSeqRegion` xs `deepSeqRegions'` y
+
+
+-- | Ensure a range is fully evaluated.
+infixr 0 `deepSeqRange`
+deepSeqRange :: Shape sh => Range sh -> b -> b
+{-# INLINE deepSeqRange #-}
+deepSeqRange range x
+ = case range of
+       RangeAll                -> x
+       RangeRects f rects      -> f `seq` rects `seq` x
+
+
+-- | Ensure a Generator's structure is fully evaluated.
+infixr 0 `deepSeqGen`
+deepSeqGen :: Shape sh => Generator sh a -> b -> b
+{-# INLINE deepSeqGen #-}
+deepSeqGen gen x
+ = case gen of
+       GenManifest vec         -> vec `seq` x
+       GenCursor{}             -> x
+
+
+-- Predicates -------------------------------------------------------------------------------------
+inRange :: Shape sh => Range sh -> sh -> Bool
+{-# INLINE inRange #-}
+inRange RangeAll _             = True
+inRange (RangeRects fn _) ix   = fn ix
+
+
+-- Singletons -------------------------------------------------------------------------------------
+-- | Wrap a scalar into a singleton array.
+singleton :: Elt a => a -> Array Z a
+{-# INLINE singleton #-}
+singleton      = fromFunction Z . const
+
+-- | Take the scalar value from a singleton array.
+toScalar :: Elt a => Array Z a -> a
+{-# INLINE toScalar #-}
+toScalar arr   = arr ! Z
+
+
+-- Projections ------------------------------------------------------------------------------------
+-- | Take the extent of an array.
+extent :: Array sh a -> sh
+{-# INLINE extent #-}
+extent arr     = arrayExtent arr
+
+
+-- | Unpack an array into delayed form.
+delay  :: (Shape sh, Elt a)
+       => Array sh a
+       -> (sh, sh -> a)
+
+{-# INLINE delay #-}
+delay arr@(Array sh _)
+       = (sh, (arr !))
+
+
+-- Indexing ---------------------------------------------------------------------------------------
+-- | Get an indexed element from an array.
+--   This uses the same level of bounds checking as your Data.Vector installation.
+(!), index
+       :: forall sh a
+       .  (Shape sh, Elt a)
+       => Array sh a
+       -> sh
+       -> a
+
+{-# INLINE (!) #-}
+(!) arr ix = index arr ix
+
+{-# INLINE index #-}
+index arr ix
+ = case arr of
+       Array _ []
+        -> zero
+
+       Array sh [Region _ gen1]
+        -> indexGen sh gen1 ix
+
+       Array sh [Region r1 gen1, Region _ gen2]
+        | inRange r1 ix        -> indexGen sh gen1 ix
+        | otherwise            -> indexGen sh gen2 ix
+
+       _ -> index' arr ix
+
+
+ where {-# INLINE indexGen #-}
+       indexGen sh gen ix'
+        = case gen of
+               GenManifest vec
+                -> vec V.! (S.toIndex sh ix')
+
+               GenCursor makeCursor _ loadElem
+                -> loadElem $ makeCursor ix'
+
+       index' (Array sh (Region range gen : rs)) ix'
+        | inRange range ix     = indexGen sh gen ix'
+        | otherwise            = index' (Array sh rs) ix'
+
+        index' (Array _ []) _
+        = zero
+
+
+
+-- | Get an indexed element from an array.
+--   If the element is out of range then `Nothing`.
+(!?), safeIndex
+       :: forall sh a
+       .  (Shape sh, Elt a)
+       => Array sh a
+       -> sh
+       -> Maybe a
+
+{-# INLINE (!?) #-}
+(!?) arr ix = safeIndex arr ix
+
+
+{-# INLINE safeIndex #-}
+safeIndex arr ix
+ = case arr of
+       Array _ []
+        -> Nothing
+
+       Array sh [Region _ gen1]
+        -> indexGen sh gen1 ix
+
+       Array sh [Region r1 gen1, Region r2 gen2]
+        | inRange r1 ix        -> indexGen sh gen1 ix
+        | inRange r2 ix        -> indexGen sh gen2 ix
+        | otherwise            -> Nothing
+
+       _ -> index' arr ix
+
+
+ where {-# INLINE indexGen #-}
+       indexGen sh gen ix'
+        = case gen of
+               GenManifest vec
+                -> vec V.!? (S.toIndex sh ix')
+
+               GenCursor makeCursor _ loadElem
+                -> Just (loadElem $ makeCursor ix')
+
+       index' (Array sh (Region range gen : rs)) ix'
+        | inRange range ix     = indexGen sh gen ix'
+        | otherwise            = index' (Array sh rs) ix'
+
+        index' (Array _ []) _
+        = Nothing
+
+
+-- | Get an indexed element from an array, without bounds checking.
+--   This assumes that the regions in the array give full coverage.
+--   An array with no regions gets zero for every element.
+unsafeIndex
+       :: forall sh a
+       .  (Shape sh, Elt a)
+       => Array sh a
+       -> sh
+       -> a
+
+{-# INLINE unsafeIndex #-}
+unsafeIndex arr ix
+ = case arr of
+       Array _ []
+        -> zero
+
+       Array sh [Region _ gen1]
+        -> unsafeIndexGen sh gen1 ix
+
+       Array sh [Region r1 gen1, Region _ gen2]
+        | inRange r1 ix        -> unsafeIndexGen sh gen1 ix
+        | otherwise            -> unsafeIndexGen sh gen2 ix
+
+       _ -> unsafeIndex' arr ix
+
+ where {-# INLINE unsafeIndexGen #-}
+       unsafeIndexGen sh gen ix'
+        = case gen of
+               GenManifest vec
+                -> vec `V.unsafeIndex` (S.toIndex sh ix')
+
+               GenCursor makeCursor _ loadElem
+                -> loadElem $ makeCursor ix'
+
+       unsafeIndex' (Array sh (Region range gen : rs)) ix'
+        | inRange range ix     = unsafeIndexGen sh gen ix'
+        | otherwise            = unsafeIndex' (Array sh rs) ix'
+
+        unsafeIndex' (Array _ []) _
+        = zero
+
+
+-- Conversions ------------------------------------------------------------------------------------
+-- | Create a `Delayed` array from a function.
+fromFunction
+       :: Shape sh
+       => sh
+       -> (sh -> a)
+       -> Array sh a
+
+{-# INLINE fromFunction #-}
+fromFunction sh fnElems
+       = sh `S.deepSeq`
+         Array sh [Region
+                       RangeAll
+                       (GenCursor id addDim fnElems)]
+
+-- | Create a `Manifest` array from an unboxed `Vector`.
+--     The elements are in row-major order.
+fromVector
+       :: Shape sh
+       => sh
+       -> Vector a
+       -> Array sh a
+
+{-# INLINE fromVector #-}
+fromVector sh vec
+       = sh  `S.deepSeq` vec `seq`
+         Array sh [Region RangeAll (GenManifest vec)]
+
+
+-- | Convert a list to an array.
+--     The length of the list must be exactly the `size` of the extent given, else `error`.
+fromList
+       :: (Shape sh, Elt a)
+       => sh
+       -> [a]
+       -> Array sh a
+
+{-# INLINE fromList #-}
+fromList sh xx
+       | V.length vec /= S.size sh
+       = error $ unlines
+               [ stage ++ ".fromList: size of array shape does not match size of list"
+               , "        size of shape = " ++ (show $ S.size sh)      ++ "\n"
+               , "        size of list  = " ++ (show $ V.length vec)   ++ "\n" ]
+
+       | otherwise
+       = Array sh [Region RangeAll (GenManifest vec)]
+
+       where   vec     = V.fromList xx
+
+
+-- | Convert a `Ptr` to an `Array`. 
+--   The data is used directly, and not copied.
+--   You promise not to modify the pointed-to data any further.
+--   
+unsafeFromForeignPtr
+        :: (Shape sh, Elt a, Storable a)
+        => sh
+        -> ForeignPtr a   
+        -> Array sh a
+
+unsafeFromForeignPtr sh fptr
+ = fromFunction sh 
+        (\ix -> unsafePerformIO 
+             $  withForeignPtr fptr
+                        (\ptr -> peekElemOff ptr $ toIndex sh ix))
+
diff --git a/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/Elt.hs b/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/Elt.hs
new file mode 100644 (file)
index 0000000..722f6fd
--- /dev/null
@@ -0,0 +1,283 @@
+-- | Values that can be stored in Repa Arrays.
+{-# LANGUAGE MagicHash, UnboxedTuples, TypeSynonymInstances, FlexibleInstances #-}
+module Data.Array.Repa.Internals.Elt
+       (Elt (..))
+where
+import GHC.Prim
+import GHC.Exts
+import GHC.Types
+import GHC.Word
+import GHC.Int
+import Data.Vector.Unboxed
+
+
+-- Note that the touch# function is special because we can pass it boxed or unboxed
+-- values. The argument type has kind ?, not just * or #.
+
+-- | Element types that can be stored in Repa arrays.
+--   Repa uses `Data.Vector.Unboxed` to store the actual data. The implementation
+--   of this library is based on type families and picks an efficient, specialised
+--   representation for every element type. In particular, unboxed vectors of pairs
+--   are represented as pairs of unboxed vectors.
+class (Show a, Unbox a)        => Elt a where
+
+       -- | We use this to prevent bindings from being floated inappropriatey.
+       --   Doing a `seq` sometimes isn't enough, because the GHC simplifier can
+       --   erase these, and/or still move around the bindings.
+       touch :: a -> IO ()
+
+       -- | Generic zero value, helpful for debugging.
+       zero  :: a
+
+       -- | Generic one value, helpful for debugging.
+       one   :: a
+
+
+-- Bool -----------------------------------------------------------------------
+instance Elt Bool where
+ {-# INLINE touch #-}
+ touch b
+  = IO (\state -> case touch# b state of
+                       state' -> (# state', () #))
+
+ {-# INLINE zero #-}
+ zero = False
+
+ {-# INLINE one #-}
+ one  = True
+
+
+-- Floating -------------------------------------------------------------------
+instance Elt Float where
+ {-# INLINE touch #-}
+ touch (F# f)
+  = IO (\state -> case touch# f state of
+                       state' -> (# state', () #))
+
+ {-# INLINE zero #-}
+ zero = 0
+
+ {-# INLINE one #-}
+ one = 1
+
+
+instance Elt Double where
+ {-# INLINE touch #-}
+ touch (D# d)
+  = IO (\state -> case touch# d state of
+                       state' -> (# state', () #))
+
+ {-# INLINE zero #-}
+ zero = 0
+
+ {-# INLINE one #-}
+ one = 1
+
+
+-- Int ------------------------------------------------------------------------
+instance Elt Int where
+ {-# INLINE touch #-}
+ touch (I# i)
+  = IO (\state -> case touch# i state of
+                       state' -> (# state', () #))
+
+ {-# INLINE zero #-}
+ zero = 0
+
+ {-# INLINE one #-}
+ one = 1
+
+instance Elt Int8 where
+ {-# INLINE touch #-}
+ touch (I8# w)
+  = IO (\state -> case touch# w state of
+                       state' -> (# state', () #))
+
+ {-# INLINE zero #-}
+ zero = 0
+
+ {-# INLINE one #-}
+ one = 1
+
+
+instance Elt Int16 where
+ {-# INLINE touch #-}
+ touch (I16# w)
+  = IO (\state -> case touch# w state of
+                       state' -> (# state', () #))
+
+ {-# INLINE zero #-}
+ zero = 0
+
+ {-# INLINE one #-}
+ one = 1
+
+
+instance Elt Int32 where
+ {-# INLINE touch #-}
+ touch (I32# w)
+  = IO (\state -> case touch# w state of
+                       state' -> (# state', () #))
+
+ {-# INLINE zero #-}
+ zero = 0
+
+ {-# INLINE one #-}
+ one = 1
+
+
+instance Elt Int64 where
+ {-# INLINE touch #-}
+ touch (I64# w)
+  = IO (\state -> case touch# w state of
+                       state' -> (# state', () #))
+
+ {-# INLINE zero #-}
+ zero = 0
+
+ {-# INLINE one #-}
+ one = 1
+
+
+-- Word -----------------------------------------------------------------------
+instance Elt Word where
+ {-# INLINE touch #-}
+ touch (W# i)
+  = IO (\state -> case touch# i state of
+                       state' -> (# state', () #))
+
+ {-# INLINE zero #-}
+ zero = 0
+
+ {-# INLINE one #-}
+ one = 1
+
+
+instance Elt Word8 where
+ {-# INLINE touch #-}
+ touch (W8# w)
+  = IO (\state -> case touch# w state of
+                       state' -> (# state', () #))
+
+ {-# INLINE zero #-}
+ zero = 0
+
+ {-# INLINE one #-}
+ one = 1
+
+
+instance Elt Word16 where
+ {-# INLINE touch #-}
+ touch (W16# w)
+  = IO (\state -> case touch# w state of
+                       state' -> (# state', () #))
+
+ {-# INLINE zero #-}
+ zero = 0
+
+ {-# INLINE one #-}
+ one = 1
+
+
+instance Elt Word32 where
+ {-# INLINE touch #-}
+ touch (W32# w)
+  = IO (\state -> case touch# w state of
+                       state' -> (# state', () #))
+
+ {-# INLINE zero #-}
+ zero = 0
+
+ {-# INLINE one #-}
+ one = 1
+
+
+instance Elt Word64 where
+ {-# INLINE touch #-}
+ touch (W64# w)
+  = IO (\state -> case touch# w state of
+                       state' -> (# state', () #))
+
+ {-# INLINE zero #-}
+ zero = 0
+
+ {-# INLINE one #-}
+ one = 1
+
+
+-- Tuple ----------------------------------------------------------------------
+instance (Elt a, Elt b) => Elt (a, b) where
+ {-# INLINE touch #-}
+ touch (a, b)
+  = do touch a
+       touch b
+
+ {-# INLINE zero #-}
+ zero = (zero, zero)
+
+ {-# INLINE one #-}
+ one =  (one, one)
+
+
+instance (Elt a, Elt b, Elt c) => Elt (a, b, c) where
+ {-# INLINE touch #-}
+ touch (a, b, c)
+  = do touch a
+       touch b
+       touch c
+
+ {-# INLINE zero #-}
+ zero = (zero, zero, zero)
+
+ {-# INLINE one #-}
+ one =  (one, one, one)
+
+
+instance (Elt a, Elt b, Elt c, Elt d) => Elt (a, b, c, d) where
+ {-# INLINE touch #-}
+ touch (a, b, c, d)
+  = do touch a
+       touch b
+       touch c
+       touch d
+
+ {-# INLINE zero #-}
+ zero = (zero, zero, zero, zero)
+
+ {-# INLINE one #-}
+ one =  (one, one, one, one)
+
+
+instance (Elt a, Elt b, Elt c, Elt d, Elt e) => Elt (a, b, c, d, e) where
+ {-# INLINE touch #-}
+ touch (a, b, c, d, e)
+  = do touch a
+       touch b
+       touch c
+       touch d
+       touch e
+
+ {-# INLINE zero #-}
+ zero = (zero, zero, zero, zero, zero)
+
+ {-# INLINE one #-}
+ one =  (one, one, one, one, one)
+
+
+instance (Elt a, Elt b, Elt c, Elt d, Elt e, Elt f) => Elt (a, b, c, d, e, f) where
+ {-# INLINE touch #-}
+ touch (a, b, c, d, e, f)
+  = do touch a
+       touch b
+       touch c
+       touch d
+       touch e
+       touch f
+
+ {-# INLINE zero #-}
+ zero = (zero, zero, zero, zero, zero, zero)
+
+ {-# INLINE one #-}
+ one =  (one, one, one, one, one, one)
+
+
diff --git a/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/EvalBlockwise.hs b/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/EvalBlockwise.hs
new file mode 100644 (file)
index 0000000..98e1291
--- /dev/null
@@ -0,0 +1,153 @@
+{-# LANGUAGE BangPatterns #-}
+
+-- | Old non-cursored, blockwise filling functions.
+--   NOTE: this isn't currently used.
+module Data.Array.Repa.Internals.EvalBlockwise
+       ( fillVectorBlockwiseP
+       , fillVectorBlock
+       , fillVectorBlockP)
+where
+import Data.Array.Repa.Internals.Elt
+import Data.Array.Repa.Internals.Gang
+import Data.Vector.Unboxed.Mutable             as VM
+import GHC.Base                                        (remInt, quotInt)
+import Prelude                                 as P
+
+
+-- Blockwise filling ------------------------------------------------------------------------------
+fillVectorBlockwiseP
+       :: Elt a
+       => IOVector a           -- ^ vector to write elements into
+       -> (Int -> a)           -- ^ fn to evaluate an element at the given index
+       -> Int                  -- ^ width of image.
+       -> IO ()
+
+{-# INLINE [0] fillVectorBlockwiseP #-}
+fillVectorBlockwiseP !vec !getElemFVBP !imageWidth
+ =     gangIO theGang fillBlock
+
+ where !threads        = gangSize theGang
+       !vecLen         = VM.length vec
+       !imageHeight    = vecLen `div` imageWidth
+       !colChunkLen    = imageWidth `quotInt` threads
+       !colChunkSlack  = imageWidth `remInt`  threads
+
+
+       {-# INLINE colIx #-}
+       colIx !ix
+        | ix < colChunkSlack   = ix * (colChunkLen + 1)
+        | otherwise            = ix * colChunkLen + colChunkSlack
+
+
+       -- just give one column to each thread
+       {-# INLINE fillBlock #-}
+       fillBlock :: Int -> IO ()
+       fillBlock !ix
+        = let  !x0     = colIx ix
+               !x1     = colIx (ix + 1)
+               !y0     = 0
+               !y1     = imageHeight
+          in   fillVectorBlock vec getElemFVBP imageWidth x0 y0 x1 y1
+
+
+-- Block filling ----------------------------------------------------------------------------------
+-- | Fill a block in a 2D image, in parallel.
+--   Coordinates given are of the filled edges of the block.
+--   We divide the block into columns, and give one column to each thread.
+fillVectorBlockP
+       :: Elt a
+       => IOVector a           -- ^ vector to write elements into
+       -> (Int -> a)           -- ^ fn to evaluate an element at the given index.
+       -> Int                  -- ^ width of whole image
+       -> Int                  -- ^ x0 lower left corner of block to fill
+       -> Int                  -- ^ y0 (low x and y value)
+       -> Int                  -- ^ x1 upper right corner of block
+       -> Int                  -- ^ y1 (high x and y value, last index to fill)
+       -> IO ()
+
+{-# INLINE [0] fillVectorBlockP #-}
+fillVectorBlockP !vec !getElem !imageWidth !x0 !y0 !x1 !y1
+ =     gangIO theGang fillBlock
+ where !threads        = gangSize theGang
+       !blockWidth     = x1 - x0 + 1
+
+       -- All columns have at least this many pixels.
+       !colChunkLen    = blockWidth `quotInt` threads
+
+       -- Extra pixels that we have to divide between some of the threads.
+       !colChunkSlack  = blockWidth `remInt` threads
+
+       -- Get the starting pixel of a column in the image.
+       {-# INLINE colIx #-}
+       colIx !ix
+        | ix < colChunkSlack   = x0 + ix * (colChunkLen + 1)
+        | otherwise            = x0 + ix * colChunkLen + colChunkSlack
+
+       -- Give one column to each thread
+       {-# INLINE fillBlock #-}
+       fillBlock :: Int -> IO ()
+       fillBlock !ix
+        = let  !x0'    = colIx ix
+               !x1'    = colIx (ix + 1) - 1
+               !y0'    = y0
+               !y1'    = y1
+          in   fillVectorBlock vec getElem imageWidth x0' y0' x1' y1'
+
+
+-- | Fill a block in a 2D image.
+--   Coordinates given are of the filled edges of the block.
+fillVectorBlock
+       :: Elt a
+       => IOVector a           -- ^ vector to write elements into.
+       -> (Int -> a)           -- ^ fn to evaluate an element at the given index.
+       -> Int                  -- ^ width of whole image
+       -> Int                  -- ^ x0 lower left corner of block to fill
+       -> Int                  -- ^ y0 (low x and y value)
+       -> Int                  -- ^ x1 upper right corner of block
+       -> Int                  -- ^ y1 (high x and y value, last index to fill)
+       -> IO ()
+
+{-# INLINE [0] fillVectorBlock #-}
+fillVectorBlock !vec !getElemFVB !imageWidth !x0 !y0 !x1 !y1
+ = do  -- putStrLn $ "fillVectorBlock: " P.++ show (x0, y0, x1, y1)
+       fillBlock ixStart (ixStart + (x1 - x0))
+ where
+       -- offset from end of one line to the start of the next.
+       !ixStart        = x0 + y0 * imageWidth
+       !ixFinal        = x1 + y1 * imageWidth
+
+       {-# INLINE fillBlock #-}
+       fillBlock !ixLineStart !ixLineEnd
+        | ixLineStart > ixFinal        = return ()
+        | otherwise
+        = do   fillLine4 ixLineStart
+               fillBlock (ixLineStart + imageWidth) (ixLineEnd + imageWidth)
+
+        where  {-# INLINE fillLine4 #-}
+               fillLine4 !ix
+                | ix + 4 > ixLineEnd   = fillLine1 ix
+                | otherwise
+                = do
+                       let d0          = getElemFVB (ix + 0)
+                       let d1          = getElemFVB (ix + 1)
+                       let d2          = getElemFVB (ix + 2)
+                       let d3          = getElemFVB (ix + 3)
+
+                       touch d0
+                       touch d1
+                       touch d2
+                       touch d3
+
+                       VM.unsafeWrite vec (ix + 0) d0
+                       VM.unsafeWrite vec (ix + 1) d1
+                       VM.unsafeWrite vec (ix + 2) d2
+                       VM.unsafeWrite vec (ix + 3) d3
+                       fillLine4 (ix + 4)
+
+               {-# INLINE fillLine1 #-}
+               fillLine1 !ix
+                | ix > ixLineEnd       = return ()
+                | otherwise
+                = do   VM.unsafeWrite vec ix (getElemFVB ix)
+                       fillLine1 (ix + 1)
+
diff --git a/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/EvalChunked.hs b/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/EvalChunked.hs
new file mode 100644 (file)
index 0000000..4f70653
--- /dev/null
@@ -0,0 +1,62 @@
+-- | Evaluate a vector by breaking it up into linear chunks and filling each chunk
+--   in parallel.
+{-# LANGUAGE BangPatterns #-}
+module Data.Array.Repa.Internals.EvalChunked
+       ( fillChunkedS
+       , fillChunkedP)
+where
+import Data.Array.Repa.Internals.Gang
+import GHC.Base                                        (remInt, quotInt)
+import Prelude                                 as P
+
+
+-- | Fill something sequentially.
+fillChunkedS
+       :: Int                  -- ^ Number of elements
+       -> (Int -> a -> IO ())  -- ^ Update function to write into result buffer
+       -> (Int -> a)           -- ^ Fn to get the value at a given index.
+       -> IO ()
+
+{-# INLINE [0] fillChunkedS #-}
+fillChunkedS !len !write !getElem
+ = fill 0
+ where fill !ix
+        | ix >= len    = return ()
+        | otherwise
+        = do   write ix (getElem ix)
+               fill (ix + 1)
+
+
+-- | Fill something in parallel.
+fillChunkedP
+        :: Int                  -- ^ Number of elements
+       -> (Int -> a -> IO ())  -- ^ Update function to write into result buffer
+       -> (Int -> a)           -- ^ Fn to get the value at a given index.
+       -> IO ()
+
+{-# INLINE [0] fillChunkedP #-}
+fillChunkedP !len !write !getElem
+ =     gangIO theGang
+        $  \thread -> fill (splitIx thread) (splitIx (thread + 1))
+
+ where
+       -- Decide now to split the work across the threads.
+       -- If the length of the vector doesn't divide evenly among the threads,
+       -- then the first few get an extra element.
+       !threads        = gangSize theGang
+       !chunkLen       = len `quotInt` threads
+       !chunkLeftover  = len `remInt`  threads
+
+       {-# INLINE splitIx #-}
+       splitIx thread
+        | thread < chunkLeftover = thread * (chunkLen + 1)
+        | otherwise              = thread * chunkLen  + chunkLeftover
+
+       -- Evaluate the elements of a single chunk.
+       {-# INLINE fill #-}
+       fill !ix !end
+        | ix >= end            = return ()
+        | otherwise
+        = do   write ix (getElem ix)
+               fill (ix + 1) end
+
diff --git a/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/EvalCursored.hs b/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/EvalCursored.hs
new file mode 100644 (file)
index 0000000..74ca6ab
--- /dev/null
@@ -0,0 +1,136 @@
+
+{-# LANGUAGE BangPatterns, UnboxedTuples #-}
+module Data.Array.Repa.Internals.EvalCursored
+       ( fillCursoredBlock2P
+       , fillCursoredBlock2 )
+where
+import Data.Array.Repa.Index
+import Data.Array.Repa.Internals.Elt
+import Data.Array.Repa.Internals.Gang
+import GHC.Base                                        (remInt, quotInt)
+import Prelude                                 as P
+
+
+-- Block filling ----------------------------------------------------------------------------------
+-- | Fill a block in a 2D image, in parallel.
+--   Coordinates given are of the filled edges of the block.
+--   We divide the block into columns, and give one column to each thread.
+fillCursoredBlock2P
+       :: Elt a
+       => (Int -> a -> IO ())          -- ^ Update function to write into result buffer
+       -> (DIM2   -> cursor)           -- ^ make a cursor to a particular element
+       -> (DIM2   -> cursor -> cursor) -- ^ shift the cursor by an offset
+       -> (cursor -> a)                -- ^ fn to evaluate an element at the given index.
+       -> Int                  -- ^ width of whole image
+       -> Int                  -- ^ x0 lower left corner of block to fill
+       -> Int                  -- ^ y0 (low x and y value)
+       -> Int                  -- ^ x1 upper right corner of block to fill
+       -> Int                  -- ^ y1 (high x and y value, index of last elem to fill)
+       -> IO ()
+
+{-# INLINE [0] fillCursoredBlock2P #-}
+fillCursoredBlock2P
+       !write
+       !makeCursorFCB !shiftCursorFCB !getElemFCB
+       !imageWidth !x0 !y0 !x1 !y1
+ =     gangIO theGang fillBlock
+ where !threads        = gangSize theGang
+       !blockWidth     = x1 - x0 + 1
+
+       -- All columns have at least this many pixels.
+       !colChunkLen    = blockWidth `quotInt` threads
+
+       -- Extra pixels that we have to divide between some of the threads.
+       !colChunkSlack  = blockWidth `remInt` threads
+
+       -- Get the starting pixel of a column in the image.
+       {-# INLINE colIx #-}
+       colIx !ix
+        | ix < colChunkSlack   = x0 + ix * (colChunkLen + 1)
+        | otherwise            = x0 + ix * colChunkLen + colChunkSlack
+
+       -- Give one column to each thread
+       {-# INLINE fillBlock #-}
+       fillBlock :: Int -> IO ()
+       fillBlock !ix
+        = let  !x0'    = colIx ix
+               !x1'    = colIx (ix + 1) - 1
+               !y0'    = y0
+               !y1'    = y1
+          in   fillCursoredBlock2
+                       write
+                       makeCursorFCB shiftCursorFCB getElemFCB
+                       imageWidth x0' y0' x1' y1'
+
+
+-- | Fill a block in a 2D image.
+--   Coordinates given are of the filled edges of the block.
+fillCursoredBlock2
+       :: Elt a
+       => (Int -> a -> IO ())          -- ^ Update function to write into result buffer
+       -> (DIM2   -> cursor)           -- ^ make a cursor to a particular element
+       -> (DIM2   -> cursor -> cursor) -- ^ shift the cursor by an offset
+       -> (cursor -> a)                -- ^ fn to evaluate an element at the given index.
+       -> Int                          -- ^ width of whole image
+       -> Int                          -- ^ x0 lower left corner of block to fill
+       -> Int                          -- ^ y0 (low x and y value)
+       -> Int                          -- ^ x1 upper right corner of block to fill
+       -> Int                          -- ^ y1 (high x and y value, index of last elem to fill)
+       -> IO ()
+
+{-# INLINE [0] fillCursoredBlock2 #-}
+fillCursoredBlock2
+       !write
+       !makeCursor !shiftCursor !getElem
+       !imageWidth !x0 !y0 !x1 !y1
+
+ = fillBlock y0
+
+ where {-# INLINE fillBlock #-}
+       fillBlock !y
+        | y > y1       = return ()
+        | otherwise
+        = do   fillLine4 x0
+               fillBlock (y + 1)
+
+        where  {-# INLINE fillLine4 #-}
+               fillLine4 !x
+                | x + 4 > x1           = fillLine1 x
+                | otherwise
+                = do   -- Compute each source cursor based on the previous one so that
+                       -- the variable live ranges in the generated code are shorter.
+                       let srcCur0     = makeCursor  (Z :. y :. x)
+                       let srcCur1     = shiftCursor (Z :. 0 :. 1) srcCur0
+                       let srcCur2     = shiftCursor (Z :. 0 :. 1) srcCur1
+                       let srcCur3     = shiftCursor (Z :. 0 :. 1) srcCur2
+
+                       -- Get the result value for each cursor.
+                       let val0        = getElem srcCur0
+                       let val1        = getElem srcCur1
+                       let val2        = getElem srcCur2
+                       let val3        = getElem srcCur3
+
+                       -- Ensure that we've computed each of the result values before we
+                       -- write into the array. If the backend code generator can't tell
+                       -- our destination array doesn't alias with the source then writing
+                       -- to it can prevent the sharing of intermediate computations.
+                       touch val0
+                       touch val1
+                       touch val2
+                       touch val3
+
+                       -- Compute cursor into destination array.
+                       let !dstCur0    = x + y * imageWidth
+                       write (dstCur0)     val0
+                       write (dstCur0 + 1) val1
+                       write (dstCur0 + 2) val2
+                       write (dstCur0 + 3) val3
+                       fillLine4 (x + 4)
+
+               {-# INLINE fillLine1 #-}
+               fillLine1 !x
+                | x > x1               = return ()
+                | otherwise
+                = do   write (x + y * imageWidth) (getElem $ makeCursor (Z :. y :. x))
+                       fillLine1 (x + 1)
+
diff --git a/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/EvalReduction.hs b/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/EvalReduction.hs
new file mode 100644 (file)
index 0000000..ddb53cb
--- /dev/null
@@ -0,0 +1,121 @@
+{-# LANGUAGE BangPatterns #-}
+module Data.Array.Repa.Internals.EvalReduction 
+        ( foldS,    foldP
+        , foldAllS, foldAllP)
+where
+import Data.Array.Repa.Internals.Elt
+import Data.Array.Repa.Internals.Gang
+import qualified Data.Vector.Unboxed            as V
+import qualified Data.Vector.Unboxed.Mutable    as M
+import GHC.Base                                 ( quotInt, divInt )
+
+
+-- | Sequential reduction of a multidimensional array along the innermost dimension.
+foldS :: Elt a
+      => M.IOVector a           -- ^ vector to write elements into
+      -> (Int -> a)             -- ^ function to get an element from the given index
+      -> (a -> a -> a)          -- ^ binary associative combination function
+      -> a                      -- ^ starting value (typically an identity)
+      -> Int                    -- ^ inner dimension (length to fold over)
+      -> IO ()
+{-# INLINE foldS #-}
+foldS vec !f !c !r !n = iter 0 0
+  where
+    !end = M.length vec
+
+    {-# INLINE iter #-}
+    iter !sh !sz | sh >= end = return ()
+                 | otherwise =
+                     let !next = sz + n
+                     in  M.unsafeWrite vec sh (reduce f c r sz next) >> iter (sh+1) next
+
+
+-- | Parallel reduction of a multidimensional array along the innermost dimension.
+--   Each output value is computed by a single thread, with the output values
+--   distributed evenly amongst the available threads.
+foldP :: Elt a
+      => M.IOVector a           -- ^ vector to write elements into
+      -> (Int -> a)             -- ^ function to get an element from the given index
+      -> (a -> a -> a)          -- ^ binary associative combination operator 
+      -> a                      -- ^ starting value. Must be neutral with respect
+                                -- ^ to the operator. eg @0 + a = a@.
+      -> Int                    -- ^ inner dimension (length to fold over)
+      -> IO ()
+{-# INLINE foldP #-}
+foldP vec !f !c !r !n
+  = gangIO theGang
+  $ \tid -> fill (split tid) (split (tid+1))
+  where
+    !threads  = gangSize theGang
+    !len      = M.length vec
+    !step     = (len + threads - 1) `quotInt` threads
+
+    {-# INLINE split #-}
+    split !ix = len `min` (ix * step)
+
+    {-# INLINE fill #-}
+    fill !start !end = iter start (start * n)
+      where
+        {-# INLINE iter #-}
+        iter !sh !sz | sh >= end = return ()
+                     | otherwise =
+                         let !next = sz + n
+                         in  M.unsafeWrite vec sh (reduce f c r sz next) >> iter (sh+1) next
+
+
+-- | Sequential reduction of all the elements in an array.
+foldAllS :: Elt a
+         => (Int -> a)          -- ^ function to get an element from the given index
+         -> (a -> a -> a)       -- ^ binary associative combining function
+         -> a                   -- ^ starting value
+         -> Int                 -- ^ number of elements
+         -> IO a
+{-# INLINE foldAllS #-}
+foldAllS !f !c !r !len = return $! reduce f c r 0 len
+
+
+-- | Parallel tree reduction of an array to a single value. Each thread takes an
+--   equally sized chunk of the data and computes a partial sum. The main thread
+--   then reduces the array of partial sums to the final result.
+--
+--   We don't require that the initial value be a neutral element, so each thread
+--   computes a fold1 on its chunk of the data, and the seed element is only
+--   applied in the final reduction step.
+--
+foldAllP :: Elt a
+         => (Int -> a)          -- ^ function to get an element from the given index
+         -> (a -> a -> a)       -- ^ binary associative combining function
+         -> a                   -- ^ starting value
+         -> Int                 -- ^ number of elements
+         -> IO a
+{-# INLINE foldAllP #-}
+foldAllP !f !c !r !len
+  | len == 0    = return r
+  | otherwise   = do
+      mvec <- M.unsafeNew chunks
+      gangIO theGang $ \tid -> fill mvec tid (split tid) (split (tid+1))
+      vec  <- V.unsafeFreeze mvec
+      return $! V.foldl' c r vec
+  where
+    !threads    = gangSize theGang
+    !step       = (len + threads - 1) `quotInt` threads
+    chunks      = ((len + step - 1) `divInt` step) `min` threads
+
+    {-# INLINE split #-}
+    split !ix   = len `min` (ix * step)
+
+    {-# INLINE fill #-}
+    fill !mvec !tid !start !end
+      | start >= end = return ()
+      | otherwise    = M.unsafeWrite mvec tid (reduce f c (f start) (start+1) end)
+
+
+-- | Sequentially reduce values between the given indices
+{-# INLINE reduce #-}
+reduce :: (Int -> a) -> (a -> a -> a) -> a -> Int -> Int -> a
+reduce !f !c !r !start !end = iter start r
+  where
+    {-# INLINE iter #-}
+    iter !i !z | i >= end  = z
+               | otherwise = iter (i+1) (f i `c` z)
+
diff --git a/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/Forcing.hs b/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/Forcing.hs
new file mode 100644 (file)
index 0000000..5d27ad3
--- /dev/null
@@ -0,0 +1,215 @@
+{-# LANGUAGE BangPatterns #-}
+module Data.Array.Repa.Internals.Forcing
+       ( toVector
+       , toList
+       , force,  forceWith
+       , force2, forceWith2)
+where
+import Data.Array.Repa.Internals.EvalChunked
+import Data.Array.Repa.Internals.EvalCursored
+import Data.Array.Repa.Internals.Elt
+import Data.Array.Repa.Internals.Base
+import Data.Array.Repa.Index
+import Data.Array.Repa.Shape                   as S
+import qualified Data.Vector.Unboxed           as V
+import qualified Data.Vector.Unboxed.Mutable   as VM
+import Data.Vector.Unboxed                     (Vector)
+import System.IO.Unsafe
+
+stage  = "Data.Array.Repa.Internals.Forcing"
+
+
+-- Conversions that also force the array ----------------------------------------------------------
+-- | Convert an array to an unboxed `Data.Vector`, forcing it if required.
+--     The elements come out in row-major order.
+toVector
+       :: (Shape sh, Elt a)
+       => Array sh a
+       -> Vector a
+{-# INLINE toVector #-}
+toVector arr
+ = case force arr of
+       Array _ [Region _ (GenManifest vec)]    -> vec
+       _       -> error $ stage ++ ".toVector: force failed"
+
+
+-- | Convert an array to a list, forcing it if required.
+toList         :: (Shape sh, Elt a)
+       => Array sh a
+       -> [a]
+
+{-# INLINE toList #-}
+toList arr
+ = V.toList $ toVector arr
+
+
+-- Forcing ----------------------------------------------------------------------------------------
+-- | Force an array, so that it becomes `Manifest`.
+--   The array is split into linear chunks and each chunk evaluated in parallel.
+force  :: (Shape sh, Elt a)
+       => Array sh a -> Array sh a
+
+{-# INLINE [2] force #-}
+force arr
+ = unsafePerformIO
+ $ do  (sh, vec)       <- forceIO arr
+       return $ sh `seq` vec `seq`
+                Array sh [Region RangeAll (GenManifest vec)]
+
+ where forceIO arr'
+        = case arr' of
+               -- Don't force an already forced array.
+               Array sh [Region RangeAll (GenManifest vec)]
+                ->     return (sh, vec)
+
+               Array sh _
+                -> do  mvec    <- VM.unsafeNew (S.size sh)
+                        forceWith (VM.unsafeWrite mvec) arr'
+                       vec     <- V.unsafeFreeze mvec
+                       return  (sh, vec)
+
+
+-- | Force an array, passing elements to the provided update function.
+--   Provide something like @(Foreign.Ptr.pokeElemOff ptr)@ to write elements into a buffer.
+--   The array is split into linear chunks and each chunk is evaluated in parallel.
+forceWith
+        :: (Shape sh, Elt a)
+        => (Int -> a -> IO ())
+        -> Array sh a
+        -> IO ()
+
+{-# INLINE [2] forceWith #-}        
+forceWith !update arr@(Array sh _)
+        = fillChunkedP  
+                (S.size sh)
+               update
+               (\ix -> arr `unsafeIndex` fromIndex sh ix)
+
+
+-- | Force an array, so that it becomes `Manifest`.
+--   This forcing function is specialised for DIM2 arrays, and does blockwise filling.
+force2 :: Elt a => Array DIM2 a -> Array DIM2 a
+{-# INLINE [2] force2 #-}
+force2 arr
+ = unsafePerformIO
+ $ do  (sh, vec)       <- forceIO2 arr
+       return $ sh `seq` vec `seq`
+                Array sh [Region RangeAll (GenManifest vec)]
+
+ where forceIO2 arr'
+        = arr' `deepSeqArray`
+          case arr' of
+               -- Don't force an already forced array.
+               Array sh [Region RangeAll (GenManifest vec)]
+                ->     return (sh, vec)
+
+               -- Create a vector to hold the new array and load in the regions.
+               Array sh _
+                -> do  mvec    <- VM.new (S.size sh)
+                        forceWith2 (VM.unsafeWrite mvec) arr'
+                        vec     <- V.unsafeFreeze mvec
+                        return (sh, vec)
+
+
+-- | Force an array, passing elements to the provided update function.
+--   Provide something like @(Foreign.Ptr.pokeElemOff ptr)@ to write elements into a buffer.
+--   This forcing function is specialised for DIM2 arrays, and does blockwise filling.
+forceWith2
+        :: Elt a
+        => (Int -> a -> IO ())
+        -> Array DIM2 a
+        -> IO ()
+
+{-# INLINE [2] forceWith2 #-}
+forceWith2 !write arr
+ = arr `deepSeqArray`
+   case arr of
+       -- If the array is already manifest then copy it into the buffer.
+       -- We don't need a particular traversal order just for a copy.
+       Array _ [Region RangeAll (GenManifest _)]
+        -> forceWith write arr
+
+       -- NOTE We must specialise this for common numbers of regions so that
+       --      we get fusion for them. If we just have the last case (arbitrary
+       --      region list) then the worker won't fuse with the filling /
+       --      evaluation code.
+       Array sh [r1]
+        -> do  fillRegion2P write sh r1
+
+       Array sh [r1, r2]
+        -> do  fillRegion2P write sh r1
+               fillRegion2P write sh r2
+
+       Array sh regions
+        -> do  mapM_ (fillRegion2P write sh) regions
+
+
+-- FillRegion2P -----------------------------------------------------------------------------------
+-- | Fill an array region into a vector.
+--   This is specialised for DIM2 regions.
+--   The region is evaluated in parallel in a blockwise manner, where each block is
+--   evaluated independently and in a separate thread. For delayed or cursored regions
+--   access their source elements from the local neighbourhood, this specialised version
+--   should given better cache performance than plain `fillRegionP`.
+--
+fillRegion2P
+       :: Elt a
+       => (Int -> a -> IO ())  -- ^ Update function to write into result buffer
+       -> DIM2                 -- ^ Extent of entire array.
+       -> Region DIM2 a        -- ^ Region to fill.
+       -> IO ()
+
+{-# INLINE [1] fillRegion2P #-}
+fillRegion2P write sh@(_ :. height :. width) (Region range gen)
+ = write `seq` height `seq` width `seq`
+   case range of
+       RangeAll
+        -> fillRect2 write sh gen
+               (Rect   (Z :. 0          :. 0)
+                       (Z :. height - 1 :. width - 1))
+
+       RangeRects _ [r1]
+        -> do  fillRect2 write sh gen r1
+
+       RangeRects _ [r1, r2]
+        -> do  fillRect2 write sh gen r1
+               fillRect2 write sh gen r2
+
+       RangeRects _ [r1, r2, r3]
+        -> do  fillRect2 write sh gen r1
+               fillRect2 write sh gen r2
+               fillRect2 write sh gen r3
+
+       RangeRects _ [r1, r2, r3, r4]
+        -> do  fillRect2 write sh gen r1
+               fillRect2 write sh gen r2
+               fillRect2 write sh gen r3
+               fillRect2 write sh gen r4
+
+       RangeRects _ rects
+        -> mapM_ (fillRect2 write sh gen) rects
+
+
+-- | Fill a rectangle in a vector.
+fillRect2
+       :: Elt a
+       => (Int -> a -> IO ())  -- ^ Update function to write into result buffer
+       -> DIM2                 -- ^ Extent of entire array.
+       -> Generator DIM2 a     -- ^ Generator for array elements.
+       -> Rect DIM2            -- ^ Rectangle to fill.
+       -> IO ()
+
+{-# INLINE fillRect2 #-}
+fillRect2 write sh@(_ :. _ :. width) gen (Rect (Z :. y0 :. x0) (Z :. y1 :. x1))
+ = write `seq` width `seq` y0 `seq` x0 `seq` y1 `seq` x1 `seq`
+   case gen of
+       GenManifest vec
+        -> fillCursoredBlock2P write
+               id addDim (\ix -> vec `V.unsafeIndex` toIndex sh ix)
+               width x0 y0 x1 y1
+
+       -- Cursor based arrays.
+       GenCursor makeCursor shiftCursor loadElem
+         -> fillCursoredBlock2P write
+               makeCursor shiftCursor loadElem
+               width x0 y0 x1 y1
diff --git a/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/Gang.hs b/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/Gang.hs
new file mode 100644 (file)
index 0000000..e1ffd23
--- /dev/null
@@ -0,0 +1,249 @@
+{-# LANGUAGE CPP #-}
+
+-- | Gang Primitives.
+--   Based on DPH code by Roman Leshchinskiy
+--
+--   Gang primitives.
+--
+#define TRACE_GANG 0
+
+module Data.Array.Repa.Internals.Gang
+       ( Gang, seqGang, forkGang, gangSize, gangIO, gangST, traceGang, traceGangST
+       , theGang)
+where
+import GHC.IO
+import GHC.ST
+import GHC.Conc                  (forkOn)
+
+import Control.Concurrent.MVar
+import Control.Exception         (assert)
+
+import Control.Monad             (zipWithM, zipWithM_)
+import GHC.Conc                        (numCapabilities)
+import System.IO
+
+#if TRACE_GANG
+import GHC.Exts                  (traceEvent)
+import System.Time ( ClockTime(..), getClockTime )
+#endif
+
+-- TheGang ----------------------------------------------------------------------------------------
+-- | The gang is shared by all computations.
+theGang :: Gang
+{-# NOINLINE theGang #-}
+theGang = unsafePerformIO $ forkGang numCapabilities
+
+
+-- Requests ---------------------------------------------------------------------------------------
+-- | The 'Req' type encapsulates work requests for individual members of a gang.
+data Req
+       -- | Instruct the worker to run the given action then signal it's done
+       --   by writing to the MVar.
+       = ReqDo        (Int -> IO ()) (MVar ())
+
+       -- | Tell the worker that we're shutting the gang down. The worker should
+        --   signal that it's received the request down by writing to the MVar
+        --   before returning to its caller (forkGang)
+       | ReqShutdown  (MVar ())
+
+
+-- | Create a new request for the given action.
+newReq :: (Int -> IO ()) -> IO Req
+newReq p
+ = do  mv      <- newEmptyMVar
+       return  $ ReqDo p mv
+
+
+-- | Block until a thread request has been executed.
+--   NOTE: only one thread can wait for the request.
+waitReq :: Req -> IO ()
+waitReq req
+ = case req of
+       ReqDo     _ varDone     -> takeMVar varDone
+       ReqShutdown varDone     -> takeMVar varDone
+
+
+-- Gang ------------------------------------------------------------------------------------------
+-- | A 'Gang' is a group of threads which execute arbitrary work requests.
+--   To get the gang to do work, write Req-uest values to its MVars
+data Gang
+       = Gang !Int           -- Number of 'Gang' threads
+               [MVar Req]     -- One 'MVar' per thread
+               (MVar Bool)    -- Indicates whether the 'Gang' is busy
+
+
+instance Show Gang where
+  showsPrec p (Gang n _ _)
+       = showString "<<"
+        . showsPrec p n
+        . showString " threads>>"
+
+
+-- | A sequential gang has no threads.
+seqGang :: Gang -> Gang
+seqGang (Gang n _ mv) = Gang n [] mv
+
+
+-- | The worker thread of a 'Gang'.
+--   The threads blocks on the MVar waiting for a work request.
+gangWorker :: Int -> MVar Req -> IO ()
+gangWorker threadId varReq
+ = do  traceGang $ "Worker " ++ show threadId ++ " waiting for request."
+       req     <- takeMVar varReq
+
+       case req of
+        ReqDo action varDone
+         -> do traceGang $ "Worker " ++ show threadId ++ " begin"
+               start   <- getGangTime
+               action threadId
+               end     <- getGangTime
+               traceGang $ "Worker " ++ show threadId ++ " end (" ++ diffTime start end ++ ")"
+
+               putMVar varDone ()
+               gangWorker threadId varReq
+
+        ReqShutdown varDone
+         -> do traceGang $ "Worker " ++ show threadId ++ " shutting down."
+               putMVar varDone ()
+
+
+-- | Finaliser for worker threads.
+--   We want to shutdown the corresponding thread when it's MVar becomes unreachable.
+--     Without this Repa programs can complain about "Blocked indefinitely on an MVar"
+--     because worker threads are still blocked on the request MVars when the program ends.
+--     Whether the finalizer is called or not is very racey. It happens about 1 in 10 runs
+--     when for the repa-edgedetect benchmark, and less often with the others.
+--
+--   We're relying on the comment in System.Mem.Weak that says
+--    "If there are no other threads to run, the runtime system will check for runnable
+--     finalizers before declaring the system to be deadlocked."
+--
+--   If we were creating and destroying the gang cleanly we wouldn't need this, but theGang
+--     is created with a top-level unsafePerformIO. Hacks beget hacks beget hacks...
+--
+finaliseWorker :: MVar Req -> IO ()
+finaliseWorker varReq
+ = do  varDone <- newEmptyMVar
+       putMVar varReq (ReqShutdown varDone)
+       takeMVar varDone
+       return ()
+
+
+-- | Fork a 'Gang' with the given number of threads (at least 1).
+forkGang :: Int -> IO Gang
+forkGang n
+ = assert (n > 0)
+ $ do
+       -- Create the vars we'll use to issue work requests.
+       mvs     <- sequence . replicate n $ newEmptyMVar
+
+       -- Add finalisers so we can shut the workers down cleanly if they become unreachable.
+       mapM_ (\var -> addMVarFinalizer var (finaliseWorker var)) mvs
+
+       -- Create all the worker threads
+       zipWithM_ forkOn [0..]
+               $ zipWith gangWorker [0 .. n-1] mvs
+
+       -- The gang is currently idle.
+       busy    <- newMVar False
+
+       return $ Gang n mvs busy
+
+
+-- | The number of threads in the 'Gang'.
+gangSize :: Gang -> Int
+gangSize (Gang n _ _) = n
+
+
+-- | Issue work requests for the 'Gang' and wait until they have been executed.
+--   If the gang is already busy then just run the action in the
+--   requesting thread.
+--
+--   TODO: We might want to print a configurable warning that this is happening.
+--
+gangIO :: Gang
+       -> (Int -> IO ())
+       -> IO ()
+
+{-# NOINLINE gangIO #-}
+gangIO (Gang n mvs busy) p
+ = do  traceGang   "gangIO: issuing work requests (SEQ_IF_GANG_BUSY)"
+       b <- swapMVar busy True
+
+       traceGang $ "gangIO: gang is currently " ++ (if b then "busy" else "idle")
+       if b
+        then do
+               hPutStr stderr
+                $ unlines      [ "Data.Array.Repa: Performing nested parallel computation sequentially."
+                               , "  You've probably called the 'force' function while another instance was"
+                               , "  already running. This can happen if the second version was suspended due"
+                               , "  to lazy evaluation. Use 'deepSeqArray' to ensure that each array is fully"
+                               , "  evaluated before you 'force' the next one."
+                               , "" ]
+
+               mapM_ p [0 .. n-1]
+
+        else do
+               parIO n mvs p
+               _ <- swapMVar busy False
+               return ()
+
+
+-- | Issue some requests to the worker threads and wait for them to complete.
+parIO  :: Int                  -- ^ Number of threads in the gang.
+       -> [MVar Req]           -- ^ Request vars for worker threads.
+       -> (Int -> IO ())       -- ^ Action to run in all the workers, it's given the ix of
+                               --   the particular worker thread it's running on.
+       -> IO ()
+
+parIO n mvs p
+ = do  traceGang "parIO: begin"
+
+       start   <- getGangTime
+       reqs    <- sequence . replicate n $ newReq p
+
+       traceGang "parIO: issuing requests"
+       _ <- zipWithM putMVar mvs reqs
+
+       traceGang "parIO: waiting for requests to complete"
+       mapM_ waitReq reqs
+       end     <- getGangTime
+
+       traceGang $ "parIO: end " ++ diffTime start end
+
+
+-- | Same as 'gangIO' but in the 'ST' monad.
+gangST :: Gang -> (Int -> ST s ()) -> ST s ()
+gangST g p = unsafeIOToST . gangIO g $ unsafeSTToIO . p
+
+
+-- Tracing ----------------------------------------------------------------------------------------
+#if TRACE_GANG
+getGangTime :: IO Integer
+getGangTime
+ = do  TOD sec pico <- getClockTime
+       return (pico + sec * 1000000000000)
+
+diffTime :: Integer -> Integer -> String
+diffTime x y = show (y-x)
+
+traceGang :: String -> IO ()
+traceGang s
+ = do  t <- getGangTime
+       traceEvent $ show t ++ " @ " ++ s
+
+#else
+getGangTime :: IO ()
+getGangTime = return ()
+
+diffTime :: () -> () -> String
+diffTime _ _ = ""
+
+traceGang :: String -> IO ()
+traceGang _ = return ()
+
+#endif
+
+traceGangST :: String -> ST s ()
+traceGangST s = unsafeIOToST (traceGang s)
+
diff --git a/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/Select.hs b/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Internals/Select.hs
new file mode 100644 (file)
index 0000000..1da46df
--- /dev/null
@@ -0,0 +1,118 @@
+{-# LANGUAGE BangPatterns, ExplicitForAll, ScopedTypeVariables, PatternGuards #-}
+module Data.Array.Repa.Internals.Select
+       (selectChunkedS, selectChunkedP)
+where
+import Data.Array.Repa.Internals.Gang
+import Data.Array.Repa.Shape
+import Data.Vector.Unboxed                     as V
+import Data.Vector.Unboxed.Mutable             as VM
+import GHC.Base                                        (remInt, quotInt)
+import Prelude                                 as P
+import Control.Monad                           as P
+import Data.IORef
+
+-- | Select indices matching a predicate
+selectChunkedS
+       :: (Shape sh, Unbox a)
+       => (sh -> Bool)         -- ^ See if this predicate matches.
+       -> (sh -> a)            -- ^  .. and apply fn to the matching index
+       -> IOVector a           -- ^  .. then write the result into the vector.
+       -> sh                   -- ^ Extent of indices to apply to predicate.
+       -> IO Int               -- ^ Number of elements written to destination array.
+
+{-# INLINE selectChunkedS #-}
+selectChunkedS match produce !vDst !shSize
+ = fill 0 0
+ where lenSrc  = size shSize
+       lenDst  = VM.length vDst
+
+       fill !nSrc !nDst
+        | nSrc >= lenSrc       = return nDst
+        | nDst >= lenDst       = return nDst
+
+        | ixSrc        <- fromIndex shSize nSrc
+        , match ixSrc
+        = do   VM.unsafeWrite vDst nDst (produce ixSrc)
+               fill (nSrc + 1) (nDst + 1)
+
+        | otherwise
+        =      fill (nSrc + 1) nDst
+
+
+-- | Select indices matching a predicate, in parallel.
+--   The array is chunked up, with one chunk being given to each thread.
+--   The number of elements in the result array depends on how many threads
+--   you're running the program with.
+selectChunkedP
+       :: forall a
+       .  Unbox a
+       => (Int -> Bool)        -- ^ See if this predicate matches.
+       -> (Int -> a)           --   .. and apply fn to the matching index
+       -> Int                  -- Extent of indices to apply to predicate.
+       -> IO [IOVector a]      -- Chunks containing array elements.
+
+{-# INLINE selectChunkedP #-}
+selectChunkedP !match !produce !len
+ = do
+       -- Make IORefs that the threads will write their result chunks to.
+       -- We start with a chunk size proportial to the number of threads we have,
+       -- but the threads themselves can grow the chunks if they run out of space.
+       refs    <- P.replicateM threads
+               $ do    vec     <- VM.new $ len `div` threads
+                       newIORef vec
+
+       -- Fire off a thread to fill each chunk.
+       gangIO theGang
+        $ \thread -> makeChunk (refs !! thread)
+                       (splitIx thread)
+                       (splitIx (thread + 1) - 1)
+
+       -- Read the result chunks back from the IORefs.
+       -- If a thread had to grow a chunk, then these might not be the same ones
+       -- we created back in the first step.
+       P.mapM readIORef refs
+
+ where -- See how many threads we have available.
+       !threads        = gangSize theGang
+       !chunkLen       = len `quotInt` threads
+       !chunkLeftover  = len `remInt`  threads
+
+
+       -- Decide where to split the source array.
+       {-# INLINE splitIx #-}
+       splitIx thread
+        | thread < chunkLeftover = thread * (chunkLen + 1)
+        | otherwise              = thread * chunkLen  + chunkLeftover
+
+
+       -- Fill the given chunk with elements selected from this range of indices.
+       makeChunk :: IORef (IOVector a) -> Int -> Int -> IO ()
+       makeChunk !ref !ixSrc !ixSrcEnd
+        = do   vecDst  <- VM.new (len `div` threads)
+               vecDst' <- fillChunk ixSrc ixSrcEnd vecDst 0 (VM.length vecDst - 1)
+               writeIORef ref vecDst'
+
+
+       -- The main filling loop.
+       fillChunk :: Int -> Int -> IOVector a -> Int -> Int -> IO (IOVector a)
+       fillChunk !ixSrc !ixSrcEnd !vecDst !ixDst !ixDstEnd
+         -- If we've finished selecting elements, then slice the vector down
+         -- so it doesn't have any empty space at the end.
+        | ixSrc >= ixSrcEnd
+        =      return  $ VM.slice 0 ixDst vecDst
+
+        -- If we've run out of space in the chunk then grow it some more.
+        | ixDst >= ixDstEnd
+        = do   let ixDstEnd'   = VM.length vecDst * 2 - 1
+               vecDst'         <- VM.grow vecDst (ixDstEnd + 1)
+               fillChunk (ixSrc + 1) ixSrcEnd vecDst' (ixDst + 1) ixDstEnd'
+
+        -- We've got a maching element, so add it to the chunk.
+        | match ixSrc
+        = do   VM.unsafeWrite vecDst ixDst (produce ixSrc)
+               fillChunk (ixSrc + 1) ixSrcEnd vecDst (ixDst + 1)  ixDstEnd
+
+        -- The element doesnt match, so keep going.
+        | otherwise
+        =      fillChunk (ixSrc + 1) ixSrcEnd vecDst ixDst ixDstEnd
+
diff --git a/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Operators/IndexSpace.hs b/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Operators/IndexSpace.hs
new file mode 100644 (file)
index 0000000..9d22002
--- /dev/null
@@ -0,0 +1,165 @@
+{-# OPTIONS_HADDOCK hide #-}
+{-# LANGUAGE TypeOperators, ExplicitForAll, FlexibleContexts #-}
+
+module Data.Array.Repa.Operators.IndexSpace
+       ( reshape
+       , append, (++)
+       , transpose
+       , extend
+       , slice
+       , backpermute
+       , backpermuteDft)
+where
+import Data.Array.Repa.Index
+import Data.Array.Repa.Slice
+import Data.Array.Repa.Internals.Elt
+import Data.Array.Repa.Internals.Base
+import Data.Array.Repa.Operators.Traverse
+import Data.Array.Repa.Shape           as S
+import Prelude                         hiding ((++))
+import qualified Prelude               as P
+
+stage  = "Data.Array.Repa.Operators.IndexSpace"
+
+-- Index space transformations --------------------------------------------------------------------
+-- | Impose a new shape on the elements of an array.
+--   The new extent must be the same size as the original, else `error`.
+--
+--   TODO: This only works for arrays with a single region.
+--
+reshape        :: (Shape sh, Shape sh', Elt a)
+       => sh'
+       -> Array sh a
+       -> Array sh' a
+
+{-# INLINE reshape #-}
+reshape sh' arr
+       | not $ S.size sh' == S.size (extent arr)
+       = error $ stage P.++ ".reshape: reshaped array will not match size of the original"
+
+reshape sh' (Array sh [Region RangeAll gen])
+ = Array sh' [Region RangeAll gen']
+ where gen' = case gen of
+               GenManifest vec
+                -> GenManifest vec
+
+               GenCursor makeCursor _ loadElem
+                -> GenCursor
+                       id
+                       addDim
+                       (loadElem . makeCursor . fromIndex sh . toIndex sh')
+
+reshape _ _
+       = error $ stage P.++ ".reshape: can't reshape a partitioned array"
+
+
+-- | Append two arrays.
+--
+append, (++)
+       :: (Shape sh, Elt a)
+       => Array (sh :. Int) a
+       -> Array (sh :. Int) a
+       -> Array (sh :. Int) a
+
+{-# INLINE append #-}
+append arr1 arr2
+ = unsafeTraverse2 arr1 arr2 fnExtent fnElem
+ where
+       (_ :. n)        = extent arr1
+
+       fnExtent (sh :. i) (_  :. j)
+               = sh :. (i + j)
+
+       fnElem f1 f2 (sh :. i)
+               | i < n         = f1 (sh :. i)
+               | otherwise     = f2 (sh :. (i - n))
+
+{-# INLINE (++) #-}
+(++) arr1 arr2 = append arr1 arr2
+
+
+-- | Transpose the lowest two dimensions of an array.
+--     Transposing an array twice yields the original.
+transpose
+       :: (Shape sh, Elt a)
+       => Array (sh :. Int :. Int) a
+       -> Array (sh :. Int :. Int) a
+
+{-# INLINE transpose #-}
+transpose arr
+ = unsafeTraverse arr
+       (\(sh :. m :. n)        -> (sh :. n :.m))
+       (\f -> \(sh :. i :. j)  -> f (sh :. j :. i))
+
+
+-- | Extend an array, according to a given slice specification.
+--   (used to be called replicate).
+extend
+       :: ( Slice sl
+          , Shape (FullShape sl)
+          , Shape (SliceShape sl)
+          , Elt e)
+       => sl
+       -> Array (SliceShape sl) e
+       -> Array (FullShape sl) e
+
+{-# INLINE extend #-}
+extend sl arr
+       = backpermute
+               (fullOfSlice sl (extent arr))
+               (sliceOfFull sl)
+               arr
+
+-- | Take a slice from an array, according to a given specification.
+slice  :: ( Slice sl
+          , Shape (FullShape sl)
+          , Shape (SliceShape sl)
+          , Elt e)
+       => Array (FullShape sl) e
+       -> sl
+       -> Array (SliceShape sl) e
+
+{-# INLINE slice #-}
+slice arr sl
+       = backpermute
+               (sliceOfFull sl (extent arr))
+               (fullOfSlice sl)
+               arr
+
+
+-- | Backwards permutation of an array's elements.
+--     The result array has the same extent as the original.
+backpermute
+       :: forall sh sh' a
+       .  (Shape sh, Shape sh', Elt a)
+       => sh'                          -- ^ Extent of result array.
+       -> (sh' -> sh)                  -- ^ Function mapping each index in the result array
+                                       --      to an index of the source array.
+       -> Array sh a                   -- ^ Source array.
+       -> Array sh' a
+
+{-# INLINE backpermute #-}
+backpermute newExtent perm arr
+       = traverse arr (const newExtent) (. perm)
+
+
+-- | Default backwards permutation of an array's elements.
+--     If the function returns `Nothing` then the value at that index is taken
+--     from the default array (@arrDft@)
+backpermuteDft
+       :: forall sh sh' a
+       .  (Shape sh, Shape sh', Elt a)
+       => Array sh' a                  -- ^ Default values (@arrDft@)
+       -> (sh' -> Maybe sh)            -- ^ Function mapping each index in the result array
+                                       --      to an index in the source array.
+       -> Array sh  a                  -- ^ Source array.
+       -> Array sh' a
+
+{-# INLINE backpermuteDft #-}
+backpermuteDft arrDft fnIndex arrSrc
+       = fromFunction (extent arrDft) fnElem
+       where   fnElem ix
+                = case fnIndex ix of
+                       Just ix'        -> arrSrc ! ix'
+                       Nothing         -> arrDft ! ix
+
diff --git a/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Operators/Interleave.hs b/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Operators/Interleave.hs
new file mode 100644 (file)
index 0000000..f67c96d
--- /dev/null
@@ -0,0 +1,111 @@
+{-# OPTIONS_HADDOCK hide #-}
+{-# LANGUAGE TypeOperators, PatternGuards #-}
+
+module Data.Array.Repa.Operators.Interleave
+       ( interleave2
+       , interleave3
+       , interleave4)
+where
+import Data.Array.Repa.Index
+import Data.Array.Repa.Internals.Elt
+import Data.Array.Repa.Internals.Base
+import Data.Array.Repa.Operators.Traverse
+import Data.Array.Repa.Shape                   as S
+
+-- | Interleave the elements of two arrays.
+--   All the input arrays must have the same extent, else `error`.
+--   The lowest dimension of the result array is twice the size of the inputs.
+--
+-- @
+--  interleave2 a1 a2   b1 b2  =>  a1 b1 a2 b2
+--              a3 a4   b3 b4      a3 b3 a4 b4
+-- @
+--
+interleave2
+       :: (Shape sh, Elt a)
+       => Array (sh :. Int) a
+       -> Array (sh :. Int) a
+       -> Array (sh :. Int) a
+
+{-# INLINE interleave2 #-}
+interleave2 arr1 arr2
+ = arr1 `deepSeqArray` arr2 `deepSeqArray`
+   unsafeTraverse2 arr1 arr2 shapeFn elemFn
+ where
+       shapeFn dim1 dim2
+        | dim1 == dim2
+        , sh :. len    <- dim1
+        = sh :. (len * 2)
+
+        | otherwise
+        = error "Data.Array.Repa.interleave2: arrays must have same extent"
+
+       elemFn get1 get2 (sh :. ix)
+        = case ix `mod` 3 of
+               0       -> get1 (sh :. ix `div` 2)
+               1       -> get2 (sh :. ix `div` 2)
+               _       -> error "Data.Array.Repa.interleave2: this never happens :-P"
+
+
+-- | Interleave the elements of three arrays.
+interleave3
+       :: (Shape sh, Elt a)
+       => Array (sh :. Int) a
+       -> Array (sh :. Int) a
+       -> Array (sh :. Int) a
+       -> Array (sh :. Int) a
+
+{-# INLINE interleave3 #-}
+interleave3 arr1 arr2 arr3
+ = arr1 `deepSeqArray` arr2 `deepSeqArray` arr3 `deepSeqArray`
+   unsafeTraverse3 arr1 arr2 arr3 shapeFn elemFn
+ where
+       shapeFn dim1 dim2 dim3
+        | dim1 == dim2
+        , dim1 == dim3
+        , sh :. len    <- dim1
+        = sh :. (len * 3)
+
+        | otherwise
+        = error "Data.Array.Repa.interleave3: arrays must have same extent"
+
+       elemFn get1 get2 get3 (sh :. ix)
+        = case ix `mod` 3 of
+               0       -> get1 (sh :. ix `div` 3)
+               1       -> get2 (sh :. ix `div` 3)
+               2       -> get3 (sh :. ix `div` 3)
+               _       -> error "Data.Array.Repa.interleave3: this never happens :-P"
+
+
+-- | Interleave the elements of four arrays.
+interleave4
+       :: (Shape sh, Elt a)
+       => Array (sh :. Int) a
+       -> Array (sh :. Int) a
+       -> Array (sh :. Int) a
+       -> Array (sh :. Int) a
+       -> Array (sh :. Int) a
+
+{-# INLINE interleave4 #-}
+interleave4 arr1 arr2 arr3 arr4
+ = arr1 `deepSeqArray` arr2 `deepSeqArray` arr3 `deepSeqArray` arr4 `deepSeqArray`
+   unsafeTraverse4 arr1 arr2 arr3 arr4 shapeFn elemFn
+ where
+       shapeFn dim1 dim2 dim3 dim4
+        | dim1 == dim2
+        , dim1 == dim3
+        , dim1 == dim4
+        , sh :. len    <- dim1
+        = sh :. (len * 4)
+
+        | otherwise
+        = error "Data.Array.Repa.interleave4: arrays must have same extent"
+
+       elemFn get1 get2 get3 get4 (sh :. ix)
+        = case ix `mod` 4 of
+               0       -> get1 (sh :. ix `div` 4)
+               1       -> get2 (sh :. ix `div` 4)
+               2       -> get3 (sh :. ix `div` 4)
+               3       -> get4 (sh :. ix `div` 4)
+               _       -> error "Data.Array.Repa.interleave4: this never happens :-P"
+
diff --git a/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Operators/Mapping.hs b/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Operators/Mapping.hs
new file mode 100644 (file)
index 0000000..1b68901
--- /dev/null
@@ -0,0 +1,109 @@
+{-# OPTIONS_HADDOCK hide #-}
+{-# LANGUAGE NoMonomorphismRestriction, PatternGuards #-}
+
+module Data.Array.Repa.Operators.Mapping
+       ( map
+       , zipWith
+       , (+^)
+       , (-^)
+       , (*^)
+       , (/^))
+where
+import Data.Array.Repa.Internals.Elt
+import Data.Array.Repa.Internals.Base
+import Data.Array.Repa.Shape           as S
+import qualified Data.Vector.Unboxed   as V
+import qualified Prelude               as P
+import Prelude                         (($), (.), (+), (*), (+), (/), (-))
+
+-- | Apply a worker function to each element of an array, yielding a new array with the same extent.
+--
+--   This is specialised for arrays of up to four regions, using more breaks fusion.
+--
+map    :: (Shape sh, Elt a, Elt b)
+       => (a -> b)
+       -> Array sh a
+       -> Array sh b
+
+{-# INLINE map #-}
+map f (Array sh regions)
+ = Array sh (mapRegions regions)
+
+ where {-# INLINE mapRegions #-}
+       mapRegions rs
+        = case rs of
+               []               -> []
+               [r]              -> [mapRegion r]
+               [r1, r2]         -> [mapRegion r1, mapRegion r2]
+               [r1, r2, r3]     -> [mapRegion r1, mapRegion r2, mapRegion r3]
+               [r1, r2, r3, r4] -> [mapRegion r1, mapRegion r2, mapRegion r3, mapRegion r4]
+               _                -> mapRegions' rs
+
+       mapRegions' rs
+        = case rs of
+               []               -> []
+               (r : rs')        -> mapRegion r : mapRegions' rs'
+
+       {-# INLINE mapRegion #-}
+       mapRegion (Region range gen)
+        = Region range (mapGen gen)
+
+       {-# INLINE mapGen #-}
+       mapGen gen
+        = case gen of
+               GenManifest vec
+                -> GenCursor
+                       P.id
+                       addDim
+                       (\ix -> f $ V.unsafeIndex vec $ S.toIndex sh ix)
+
+               GenCursor makeCursor shiftCursor loadElem
+                -> GenCursor makeCursor shiftCursor (f . loadElem)
+
+
+-- | Combine two arrays, element-wise, with a binary operator.
+--     If the extent of the two array arguments differ,
+--     then the resulting array's extent is their intersection.
+--
+zipWith :: (Shape sh, Elt a, Elt b, Elt c)
+       => (a -> b -> c)
+       -> Array sh a
+       -> Array sh b
+       -> Array sh c
+
+{-# INLINE zipWith #-}
+zipWith f arr1 arr2
+       | Array sh2 [_] <- arr1
+       , Array sh1 [ Region g21 (GenCursor make21 _ load21)
+                   , Region g22 (GenCursor make22 _ load22)] <- arr2
+
+       = let   {-# INLINE load21' #-}
+               load21' ix      = f (arr1 `unsafeIndex` ix) (load21 $ make21 ix)
+
+               {-# INLINE load22' #-}
+               load22' ix      = f (arr1 `unsafeIndex` ix) (load22 $ make22 ix)
+
+         in    Array (S.intersectDim sh1 sh2)
+                     [ Region g21 (GenCursor P.id addDim load21')
+                     , Region g22 (GenCursor P.id addDim load22') ]
+
+       | P.otherwise
+       = let   {-# INLINE getElem' #-}
+               getElem' ix     = f (arr1 `unsafeIndex` ix) (arr2 `unsafeIndex` ix)
+         in    fromFunction
+                       (S.intersectDim (extent arr1) (extent arr2))
+                       getElem'
+
+
+{-# INLINE (+^) #-}
+(+^)   = zipWith (+)
+
+{-# INLINE (-^) #-}
+(-^)   = zipWith (-)
+
+{-# INLINE (*^) #-}
+(*^)   = zipWith (*)
+
+{-# INLINE (/^) #-}
+(/^)   = zipWith (/)
+
diff --git a/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Operators/Modify.hs b/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Operators/Modify.hs
new file mode 100644 (file)
index 0000000..c0334f4
--- /dev/null
@@ -0,0 +1,53 @@
+{-# OPTIONS_HADDOCK hide #-}
+
+module Data.Array.Repa.Operators.Modify 
+        ( -- * Bulk updates
+         (//))
+where
+import Data.Array.Repa.Shape
+import Data.Array.Repa.Internals.Elt
+import Data.Array.Repa.Internals.Base
+
+{-
+stage :: String
+stage = "Data.Array.Repa.Operators.Modify"
+-}
+
+-- Bulk updates ----------------------------------------------------------------
+-- | For each pair @(sh, a)@ from the list of index/value pairs, replace the
+-- element at position @sh@ by @a@.
+--
+-- > update <5,9,2,7> [(2,1),(0,3),(2,8)] = <3,9,8,7>
+--
+{-# INLINE (//) #-}
+(//) :: (Shape sh, Elt a) => Array sh a -> [(sh,a)] -> Array sh a
+(//) arr us 
+        = fromFunction
+                (extent arr) 
+                (\sh -> case lookup sh us of
+                            Just a  -> a
+                            Nothing -> index arr sh)
+
+{-
+-- For each pair @(sh, a)@ from the array of index/value pairs, replace the
+-- element at position @sh@ by @a@.
+--
+-- > update <5,9,2,7> <(2,1),(0,3),(2,8)> = <3,9,8,7>
+--
+{-# INLINE update #-}
+update :: Shape sh
+       => Array sh a            -- ^ initial array
+       -> Array sh (sh, a)      -- ^ array of shape/value pairs
+       -> Array sh a
+update _arr _us = error $ stage ++ ".update: not defined yet"
+
+
+-- Same as 'update', but without bounds checks
+--
+{-# INLINE unsafeUpdate #-}
+unsafeUpdate :: Shape sh
+             => Array sh a
+             -> Array sh (sh, a)
+             -> Array sh a
+unsafeUpdate _arr _us = error $ stage ++ ".unsafeUpdate: not defined yet"
+-}
diff --git a/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Operators/Reduction.hs b/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Operators/Reduction.hs
new file mode 100644 (file)
index 0000000..a2db9b6
--- /dev/null
@@ -0,0 +1,82 @@
+{-# OPTIONS_HADDOCK hide #-}
+{-# LANGUAGE BangPatterns, ExplicitForAll, TypeOperators #-}
+
+module Data.Array.Repa.Operators.Reduction
+       ( fold, foldAll
+       , sum,  sumAll)
+where
+import Data.Array.Repa.Index
+import Data.Array.Repa.Internals.Elt
+import Data.Array.Repa.Internals.Base
+import Data.Array.Repa.Shape                   as S
+import qualified Data.Vector.Unboxed           as V
+import qualified Data.Vector.Unboxed.Mutable    as M
+import Prelude                                 hiding (sum)
+
+import Data.Array.Repa.Internals.EvalReduction
+import System.IO.Unsafe
+
+
+-- | Reduction of the innermost dimension of an arbitrary rank array. The first
+--   argument needs to be an /associative/ operator. The starting element must
+--   be neutral with respect to the operator, for example @0@ is neutral with
+--   respect to @(+)@ as @0 + a = a@. These restrictions are required to support
+--   parallel evaluation, as the starting element may be used multiple
+--   times depending on the number of threads.
+
+--   Combine this with `transpose` to fold any other dimension.
+fold   :: (Shape sh, Elt a)
+       => (a -> a -> a)
+       -> a
+       -> Array (sh :. Int) a
+       -> Array sh a
+{-# INLINE [1] fold #-}
+fold f z arr 
+ = let  sh@(sz :. n) = extent arr
+   in   case rank sh of
+           -- specialise rank-1 arrays, else one thread does all the work. We can't
+           -- match against the shape constructor, otherwise type error: (sz ~ Z)
+           --
+           1 -> let !x = V.singleton $ foldAll f z arr
+                in  Array sz [Region RangeAll (GenManifest x)]
+
+           _ -> unsafePerformIO 
+              $ do mvec   <- M.unsafeNew (S.size sz)
+                   foldP mvec (\ix -> arr `unsafeIndex` fromIndex sh ix) f z n
+                   !vec   <- V.unsafeFreeze mvec
+                   return $ Array sz [Region RangeAll (GenManifest vec)]
+
+
+-- | Reduction of an array of arbitrary rank to a single scalar value. The first
+--   argument needs to be an /associative/ operator. The starting element must
+--   be neutral with respect to the operator, for example @0@ is neutral with
+--   respect to @(+)@ as @0 + a = a@. These restrictions are required to support
+--   parallel evaluation, as the starting element may be used multiple
+--   times depending on the number of threads.
+foldAll :: (Shape sh, Elt a)
+       => (a -> a -> a)
+       -> a
+       -> Array sh a
+       -> a
+{-# INLINE [1] foldAll #-}
+foldAll f z arr 
+ = let  sh = extent arr
+        n  = size sh
+   in   unsafePerformIO $ foldAllP (\ix -> arr `unsafeIndex` fromIndex sh ix) f z n
+
+
+-- | Sum the innermost dimension of an array.
+sum    :: (Shape sh, Elt a, Num a)
+       => Array (sh :. Int) a
+       -> Array sh a
+{-# INLINE sum #-}
+sum arr        = fold (+) 0 arr
+
+
+-- | Sum all the elements of an array.
+sumAll :: (Shape sh, Elt a, Num a)
+       => Array sh a
+       -> a
+{-# INLINE sumAll #-}
+sumAll arr = foldAll (+) 0 arr
+
diff --git a/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Operators/Select.hs b/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Operators/Select.hs
new file mode 100644 (file)
index 0000000..a7948c2
--- /dev/null
@@ -0,0 +1,44 @@
+{-# OPTIONS_HADDOCK hide #-}
+{-# LANGUAGE BangPatterns #-}
+
+module Data.Array.Repa.Operators.Select
+       (select)
+where
+import Data.Array.Repa.Index
+import Data.Array.Repa.Internals.Elt
+import Data.Array.Repa.Internals.Base
+import Data.Array.Repa.Internals.Select
+import qualified Data.Vector.Unboxed           as V
+import System.IO.Unsafe
+
+
+-- | Produce an array by applying a predicate to a range of integers.
+--   If the predicate matches, then use the second function to generate
+--   the element.
+--
+--   This is a low-level function helpful for writing filtering operations on arrays.
+--   Use the integer as the index into the array you're filtering.
+--
+select :: Elt a
+       => (Int -> Bool)        -- ^ If the Int matches this predicate,
+       -> (Int -> a)           -- ^  ... then pass it to this fn to produce a value
+       -> Int                  -- ^ Range between 0 and this maximum.
+       -> Array DIM1 a         -- ^ Array containing produced values.
+
+{-# INLINE select #-}
+select match produce len
+ = unsafePerformIO
+ $ do  (sh, vec)       <- selectIO
+       return $ sh `seq` vec `seq`
+                Array sh [Region RangeAll (GenManifest vec)]
+
+ where {-# INLINE selectIO #-}
+       selectIO
+        = do   vecs            <- selectChunkedP match produce len
+               vecs'           <- mapM V.unsafeFreeze vecs
+
+               -- TODO: avoid copy.
+               let result      = V.concat vecs'
+
+               return  (Z :. V.length result, result)
+
diff --git a/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Operators/Traverse.hs b/fibon/Repa/_RepaLib/repa/Data/Array/Repa/Operators/Traverse.hs
new file mode 100644 (file)
index 0000000..725c850
--- /dev/null
@@ -0,0 +1,126 @@
+{-# OPTIONS_HADDOCK hide #-}
+{-# LANGUAGE ExplicitForAll #-}
+
+module Data.Array.Repa.Operators.Traverse
+       ( traverse,  unsafeTraverse
+       , traverse2, unsafeTraverse2
+       , traverse3, unsafeTraverse3
+       , traverse4, unsafeTraverse4)
+where
+import Data.Array.Repa.Internals.Elt
+import Data.Array.Repa.Internals.Base
+import Data.Array.Repa.Shape   as S
+
+-- Generic Traversal -----------------------------------------------------------------------------
+-- | Unstructured traversal.
+traverse
+       :: forall sh sh' a b
+       .  (Shape sh, Shape sh', Elt a)
+       => Array sh a                           -- ^ Source array.
+       -> (sh  -> sh')                         -- ^ Function to produce the extent of the result.
+       -> ((sh -> a) -> sh' -> b)              -- ^ Function to produce elements of the result.
+                                               --   It is passed a lookup function to get elements of the source.
+       -> Array sh' b
+
+{-# INLINE traverse #-}
+traverse arr transExtent newElem
+       = arr `deepSeqArray`
+          fromFunction (transExtent (extent arr)) (newElem (arr !))
+
+
+{-# INLINE unsafeTraverse #-}
+unsafeTraverse arr transExtent newElem
+       = arr `deepSeqArray`
+         fromFunction (transExtent (extent arr)) (newElem (unsafeIndex arr))
+
+
+-- | Unstructured traversal over two arrays at once.
+traverse2, unsafeTraverse2
+       :: forall sh sh' sh'' a b c
+       .  ( Shape sh, Shape sh', Shape sh''
+          , Elt a,    Elt b)
+        => Array sh a                          -- ^ First source array.
+       -> Array sh' b                          -- ^ Second source array.
+        -> (sh -> sh' -> sh'')                 -- ^ Function to produce the extent of the result.
+        -> ((sh -> a) -> (sh' -> b)
+                      -> (sh'' -> c))          -- ^ Function to produce elements of the result.
+                                               --   It is passed lookup functions to get elements of the
+                                               --   source arrays.
+        -> Array sh'' c
+
+{-# INLINE traverse2 #-}
+traverse2 arrA arrB transExtent newElem
+       = arrA `deepSeqArray` arrB `deepSeqArray`
+         fromFunction
+               (transExtent (extent arrA) (extent arrB))
+               (newElem     (arrA !) (arrB !))
+
+{-# INLINE unsafeTraverse2 #-}
+unsafeTraverse2 arrA arrB transExtent newElem
+       = arrA `deepSeqArray` arrB `deepSeqArray`
+         fromFunction
+               (transExtent (extent arrA) (extent arrB))
+               (newElem     (unsafeIndex arrA) (unsafeIndex arrB))
+
+
+-- | Unstructured traversal over three arrays at once.
+traverse3, unsafeTraverse3
+       :: forall sh1 sh2 sh3 sh4
+                 a   b   c   d
+       .  ( Shape sh1, Shape sh2, Shape sh3, Shape sh4
+          , Elt a,     Elt b,     Elt c)
+        => Array sh1 a
+       -> Array sh2 b
+       -> Array sh3 c
+        -> (sh1 -> sh2 -> sh3 -> sh4)
+        -> (  (sh1 -> a) -> (sh2 -> b)
+           -> (sh3 -> c)
+           ->  sh4 -> d )
+        -> Array sh4 d
+
+{-# INLINE traverse3 #-}
+traverse3 arrA arrB arrC transExtent newElem
+       = arrA `deepSeqArray` arrB `deepSeqArray` arrC `deepSeqArray`
+         fromFunction
+               (transExtent (extent arrA) (extent arrB) (extent arrC))
+               (newElem     (arrA !) (arrB !) (arrC !))
+
+{-# INLINE unsafeTraverse3 #-}
+unsafeTraverse3 arrA arrB arrC transExtent newElem
+       = arrA `deepSeqArray` arrB `deepSeqArray` arrC `deepSeqArray`
+         fromFunction
+               (transExtent (extent arrA) (extent arrB) (extent arrC))
+               (newElem     (unsafeIndex arrA) (unsafeIndex arrB) (unsafeIndex arrC))
+
+
+-- | Unstructured traversal over four arrays at once.
+traverse4, unsafeTraverse4
+       :: forall sh1 sh2 sh3 sh4 sh5
+                 a   b   c   d   e
+       .  ( Shape sh1, Shape sh2, Shape sh3, Shape sh4, Shape sh5
+          , Elt a,     Elt b,     Elt c,     Elt d)
+        => Array sh1 a
+       -> Array sh2 b
+       -> Array sh3 c
+       -> Array sh4 d
+        -> (sh1 -> sh2 -> sh3 -> sh4 -> sh5 )
+        -> (  (sh1 -> a) -> (sh2 -> b)
+           -> (sh3 -> c) -> (sh4 -> d)
+           ->  sh5 -> e )
+        -> Array sh5 e
+
+{-# INLINE traverse4 #-}
+traverse4 arrA arrB arrC arrD transExtent newElem
+       = arrA `deepSeqArray` arrB `deepSeqArray` arrC `deepSeqArray` arrD `deepSeqArray`
+         fromFunction
+               (transExtent (extent arrA) (extent arrB) (extent arrC) (extent arrD))
+               (newElem     (arrA !) (arrB !) (arrC !) (arrD !))
+
+
+{-# INLINE unsafeTraverse4 #-}
+unsafeTraverse4 arrA arrB arrC arrD transExtent newElem
+       = arrA `deepSeqArray` arrB `deepSeqArray` arrC `deepSeqArray` arrD `deepSeqArray`
+         fromFunction