diff --git a/Data/Array/Repa/Index.hs b/Data/Array/Repa/Index.hs deleted file mode 100644 --- a/Data/Array/Repa/Index.hs +++ /dev/null @@ -1,171 +0,0 @@ -{-# LANGUAGE TypeOperators, FlexibleInstances, ScopedTypeVariables #-} - --- | Index types. -module Data.Array.Repa.Index - ( - -- * Index types - Z (..) - , (:.) (..) - - -- * Common dimensions. - , DIM0, DIM1, DIM2, DIM3, DIM4, DIM5 - , ix1, ix2, ix3, ix4, ix5) -where -import Data.Array.Repa.Shape -import GHC.Base (quotInt, remInt) - -stage = "Data.Array.Repa.Index" - --- | An index of dimension zero -data Z = Z - deriving (Show, Read, Eq, Ord) - --- | Our index type, used for both shapes and indices. -infixl 3 :. -data tail :. head - = !tail :. !head - deriving (Show, Read, Eq, Ord) - --- Common dimensions -type DIM0 = Z -type DIM1 = DIM0 :. Int -type DIM2 = DIM1 :. Int -type DIM3 = DIM2 :. Int -type DIM4 = DIM3 :. Int -type DIM5 = DIM4 :. Int - - --- | Helper for index construction. --- --- Use this instead of explicit constructors like @(Z :. (x :: Int))@. --- The this is sometimes needed to ensure that 'x' is constrained to --- be in @Int@. -ix1 :: Int -> DIM1 -ix1 x = Z :. x -{-# INLINE ix1 #-} - -ix2 :: Int -> Int -> DIM2 -ix2 y x = Z :. y :. x -{-# INLINE ix2 #-} - -ix3 :: Int -> Int -> Int -> DIM3 -ix3 z y x = Z :. z :. y :. x -{-# INLINE ix3 #-} - -ix4 :: Int -> Int -> Int -> Int -> DIM4 -ix4 a z y x = Z :. a :. z :. y :. x -{-# INLINE ix4 #-} - -ix5 :: Int -> Int -> Int -> Int -> Int -> DIM5 -ix5 b a z y x = Z :. b :. a :. z :. y :. x -{-# INLINE ix5 #-} - - --- Shape ---------------------------------------------------------------------- -instance Shape Z where - {-# INLINE [1] rank #-} - rank _ = 0 - - {-# INLINE [1] zeroDim #-} - zeroDim = Z - - {-# INLINE [1] unitDim #-} - unitDim = Z - - {-# INLINE [1] intersectDim #-} - intersectDim _ _ = Z - - {-# INLINE [1] addDim #-} - addDim _ _ = Z - - {-# INLINE [1] size #-} - size _ = 1 - - {-# INLINE [1] sizeIsValid #-} - sizeIsValid _ = True - - - {-# INLINE [1] toIndex #-} - toIndex _ _ = 0 - - {-# INLINE [1] fromIndex #-} - fromIndex _ _ = Z - - - {-# INLINE [1] inShapeRange #-} - inShapeRange Z Z Z = True - - {-# NOINLINE listOfShape #-} - listOfShape _ = [] - - {-# NOINLINE shapeOfList #-} - shapeOfList [] = Z - shapeOfList _ = error $ stage ++ ".fromList: non-empty list when converting to Z." - - {-# INLINE deepSeq #-} - deepSeq Z x = x - - -instance Shape sh => Shape (sh :. Int) where - {-# INLINE [1] rank #-} - rank (sh :. _) - = rank sh + 1 - - {-# INLINE [1] zeroDim #-} - zeroDim = zeroDim :. 0 - - {-# INLINE [1] unitDim #-} - unitDim = unitDim :. 1 - - {-# INLINE [1] intersectDim #-} - intersectDim (sh1 :. n1) (sh2 :. n2) - = (intersectDim sh1 sh2 :. (min n1 n2)) - - {-# INLINE [1] addDim #-} - addDim (sh1 :. n1) (sh2 :. n2) - = addDim sh1 sh2 :. (n1 + n2) - - {-# INLINE [1] size #-} - size (sh1 :. n) - = size sh1 * n - - {-# INLINE [1] sizeIsValid #-} - sizeIsValid (sh1 :. n) - | size sh1 > 0 - = n <= maxBound `div` size sh1 - - | otherwise - = False - - {-# INLINE [1] toIndex #-} - toIndex (sh1 :. sh2) (sh1' :. sh2') - = toIndex sh1 sh1' * sh2 + sh2' - - {-# INLINE [1] fromIndex #-} - fromIndex (ds :. d) n - = fromIndex ds (n `quotInt` d) :. r - where - -- If we assume that the index is in range, there is no point - -- in computing the remainder for the highest dimension since - -- n < d must hold. This saves one remInt per element access which - -- is quite a big deal. - r | rank ds == 0 = n - | otherwise = n `remInt` d - - {-# INLINE [1] inShapeRange #-} - inShapeRange (zs :. z) (sh1 :. n1) (sh2 :. n2) - = (n2 >= z) && (n2 < n1) && (inShapeRange zs sh1 sh2) - - {-# NOINLINE listOfShape #-} - listOfShape (sh :. n) - = n : listOfShape sh - - {-# NOINLINE shapeOfList #-} - shapeOfList xx - = case xx of - [] -> error $ stage ++ ".toList: empty list when converting to (_ :. Int)" - x:xs -> shapeOfList xs :. x - - {-# INLINE deepSeq #-} - deepSeq (sh :. n) x = deepSeq sh (n `seq` x) - diff --git a/Data/Array/Repa/Shape.hs b/Data/Array/Repa/Shape.hs deleted file mode 100644 --- a/Data/Array/Repa/Shape.hs +++ /dev/null @@ -1,82 +0,0 @@ -{-# LANGUAGE RankNTypes #-} - --- | Class of types that can be used as array shapes and indices. -module Data.Array.Repa.Shape - ( Shape(..) - , inShape - , showShape ) -where - --- Shape ---------------------------------------------------------------------- --- | Class of types that can be used as array shapes and indices. -class Eq sh => Shape sh where - - -- | Get the number of dimensions in a shape. - rank :: sh -> Int - - -- | The shape of an array of size zero, with a particular dimensionality. - zeroDim :: sh - - -- | The shape of an array with size one, with a particular dimensionality. - unitDim :: sh - - -- | Compute the intersection of two shapes. - intersectDim :: sh -> sh -> sh - - -- | Add the coordinates of two shapes componentwise - addDim :: sh -> sh -> sh - - -- | Get the total number of elements in an array with this shape. - size :: sh -> Int - - -- | Check whether this shape is small enough so that its flat - -- indices an be represented as `Int`. If this returns `False` then your - -- array is too big. Mostly used for writing QuickCheck tests. - sizeIsValid :: sh -> Bool - - - -- | Convert an index into its equivalent flat, linear, row-major version. - toIndex :: sh -- ^ Shape of the array. - -> sh -- ^ Index into the array. - -> Int - - -- | Inverse of `toIndex`. - fromIndex - :: sh -- ^ Shape of the array. - -> Int -- ^ Index into linear representation. - -> sh - - -- | Check whether an index is within a given shape. - inShapeRange - :: sh -- ^ Start index for range. - -> sh -- ^ Final index for range. - -> sh -- ^ Index to check for. - -> Bool - - -- | Convert a shape into its list of dimensions. - listOfShape :: sh -> [Int] - - -- | Convert a list of dimensions to a shape - shapeOfList :: [Int] -> sh - - -- | Ensure that a shape is completely evaluated. - infixr 0 `deepSeq` - deepSeq :: sh -> a -> a - - --- | Check whether an index is a part of a given shape. -inShape :: forall sh - . Shape sh - => sh -- ^ Shape of the array. - -> sh -- ^ Index. - -> Bool - -{-# INLINE inShape #-} -inShape sh ix - = inShapeRange zeroDim sh ix - - --- | Nicely format a shape as a string -showShape :: Shape sh => sh -> String -showShape = foldr (\sh str -> str ++ " :. " ++ show sh) "Z" . listOfShape - diff --git a/Data/ExtShape.hs b/Data/ExtShape.hs deleted file mode 100644 --- a/Data/ExtShape.hs +++ /dev/null @@ -1,41 +0,0 @@ -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE TypeOperators #-} - --- | Additional functions on shapes - -module Data.ExtShape where - -import Data.Array.Repa.Index -import Data.Array.Repa.Shape - - - --- | A number of additional operations that are useful together with --- 'PrimitiveArray's. - -class (Eq sh, Shape sh) => ExtShape sh where - - -- | subtract the right coordinates from the left. Does not check if the - -- resulting shape make sense. - - subDim :: sh -> sh -> sh - - -- | Given an index and an extend, return a list of all indices. For - -- @rangeList (Z:.3) (Z:.2)@ this returns @[(Z:.3), (Z:.4), (Z:.5)]@. - - rangeList :: sh -> sh -> [sh] - - - -instance ExtShape Z where - subDim _ _ = Z - {-# INLINE subDim #-} - rangeList _ _ = [Z] - {-# INLINE rangeList #-} - -instance ExtShape sh => ExtShape (sh:.Int) where - subDim (sh1:.n1) (sh2:.n2) = subDim sh1 sh2 :. (n1-n2) - {-# INLINE subDim #-} - rangeList (sh1:.n1) (sh2:.n2) = [sh:.n | sh <- rangeList sh1 sh2, n <- [n1 .. (n1+n2) ] ] - {-# INLINE rangeList #-} - diff --git a/Data/PrimitiveArray.hs b/Data/PrimitiveArray.hs --- a/Data/PrimitiveArray.hs +++ b/Data/PrimitiveArray.hs @@ -1,162 +1,11 @@ -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE MultiParamTypeClasses #-} --- | Vastly extended primitive arrays. Some basic ideas are now modeled after --- the vector package, especially the monadic mutable / pure immutable array --- system. There are eight flavors of arrays among three axes: mutable/pure + --- boxed/unboxed + zero-based/lower-bound. --- --- NOTE all operations in MPrimArrayOps and PrimArrayOps are highly unsafe. No --- bounds-checking is performed at all. - -module Data.PrimitiveArray where - -import Data.Array.Repa.Index -import Data.Array.Repa.Shape -import Data.Primitive.Types -import Data.Primitive -import Control.Monad.ST -import Control.Monad -import Control.Monad.Primitive -import System.IO.Unsafe -import Control.Exception (assert) - -import Data.ExtShape - - - --- | The core set of operations for monadic arrays. - -class (Shape sh, ExtShape sh) => MPrimArrayOps marr sh elm where - - -- | Return the bounds of the array. All bounds are inclusive, as in - -- @[lb..ub]@ - - boundsM :: marr s sh elm -> (sh,sh) - - -- | Given lower and upper bounds and a list of /all/ elements, produce a - -- mutable array. - - fromListM :: PrimMonad m => sh -> sh -> [elm] -> m (marr (PrimState m) sh elm) - - -- | Creates a new array with the given bounds with each element within the - -- array being in an undefined state. - - newM :: PrimMonad m => sh -> sh -> m (marr (PrimState m) sh elm) - - -- | Creates a new array with all elements being equal to 'elm'. - - newWithM :: PrimMonad m => sh -> sh -> elm -> m (marr (PrimState m) sh elm) - - -- | Reads a single element in the array. - - readM :: PrimMonad m => marr (PrimState m) sh elm -> sh -> m elm - - -- | Writes a single element in the array. - - writeM :: PrimMonad m => marr (PrimState m) sh elm -> sh -> elm -> m () - - - --- | Used to connect each immutable array with one mutable array. - -type family MutArray (v :: * -> * -> * ) :: * -> * -> * -> * - - - --- | The core set of functions on immutable arrays. - -class (Shape sh, ExtShape sh, MPrimArrayOps (MutArray arr) sh elm) => PrimArrayOps arr sh elm where - - -- | Returns the bounds of an immutable array, again inclusive bounds: @ [lb..ub] @. - - bounds :: arr sh elm -> (sh,sh) - - -- | Freezes a mutable array an returns its immutable version. This operation - -- is /O(1)/ and both arrays share the same memory. Do not use the mutable - -- array afterwards. - - freeze :: PrimMonad m => MutArray arr (PrimState m) sh elm -> m (arr sh elm) - - -- | Extract a single element from the array. Generally unsafe as not - -- bounds-checking is performed. - - index :: arr sh elm -> sh -> elm - - - --- | Infix index operator. Performs minimal bounds-checking using assert in --- non-optimized code. - -(!) :: PrimArrayOps arr sh elm => arr sh elm -> sh -> elm -(!) arr idx = assert (inBounds arr idx) $ index arr idx -{-# INLINE (!) #-} - --- | Returns true if the index is valid for the array. - -inBoundsM :: MPrimArrayOps marr sh elm => marr s sh elm -> sh -> Bool -inBoundsM marr idx = let (lb,ub) = boundsM marr in inShapeRange lb ub idx -{-# INLINE inBoundsM #-} - --- | Given two arrays with the same dimensionality, their respective starting --- index, and how many steps to go in each dimension (in terms of a dimension --- again), determine if the multidimensional slices have the same value at --- all positions --- --- TODO specialize for DIM1 (and maybe higher dim's) to use memcmp - -sliceEq :: (Eq elm, PrimArrayOps arr sh elm) => arr sh elm -> sh -> arr sh elm -> sh -> sh -> Bool -sliceEq arr1 k1 arr2 k2 xtnd = assert ((inBounds arr1 k1) && (inBounds arr2 k2) && (inBounds arr1 $ k1 `addDim` xtnd) && (inBounds arr2 $ k2 `addDim` xtnd)) $ and res where - res = zipWith (==) xs ys - xs = map (index arr1) $ rangeList k1 xtnd - ys = map (index arr2) $ rangeList k2 xtnd -{-# INLINE sliceEq #-} - --- | Construct a mutable primitive array from a lower and an upper bound, a --- default element, and a list of associations. - -fromAssocsM - :: (PrimMonad m, MPrimArrayOps marr sh elm) - => sh -> sh -> elm -> [(sh,elm)] -> m (marr (PrimState m) sh elm) -fromAssocsM lb ub def xs = do - ma <- newWithM lb ub def - forM_ xs $ \(k,v) -> writeM ma k v - return ma -{-# INLINE fromAssocsM #-} - --- | Return all associations from an array. - -assocs :: PrimArrayOps arr sh elm => arr sh elm -> [(sh,elm)] -assocs arr = map (\k -> (k,index arr k)) $ rangeList lb (ub `subDim` lb) where - (lb,ub) = bounds arr -{-# INLINE assocs #-} - --- | Creates an immutable array from lower and upper bounds and a complete list --- of elements. - -fromList :: PrimArrayOps arr sh elm => sh -> sh -> [elm] -> arr sh elm -fromList lb ub xs = runST $ fromListM lb ub xs >>= freeze -{-# INLINE fromList #-} - --- | Creates an immutable array from lower and upper bounds, a default element, --- and a list of associations. - -fromAssocs :: PrimArrayOps arr sh elm => sh -> sh -> elm -> [(sh,elm)] -> arr sh elm -fromAssocs lb ub def xs = runST $ fromAssocsM lb ub def xs >>= freeze -{-# INLINE fromAssocs #-} - --- | Determines if an index is valid for a given immutable array. - -inBounds :: PrimArrayOps arr sh elm => arr sh elm -> sh -> Bool -inBounds arr idx = let (lb,ub) = bounds arr in inShapeRange lb (ub `addDim` unitDim) idx -{-# INLINE inBounds #-} - --- | Returns all elements of an immutable array as a list. +module Data.PrimitiveArray + ( module Data.PrimitiveArray.Class + , module Data.PrimitiveArray.Dense + , module Data.PrimitiveArray.Index + ) where -toList :: PrimArrayOps arr sh elm => arr sh elm -> [elm] -toList arr = let (lb,ub) = bounds arr in map ((!) arr) $ rangeList lb $ ub `subDim` lb -{-# INLINE toList #-} +import Data.PrimitiveArray.Class +import Data.PrimitiveArray.Dense +import Data.PrimitiveArray.Index diff --git a/Data/PrimitiveArray/Checked.hs b/Data/PrimitiveArray/Checked.hs new file mode 100644 --- /dev/null +++ b/Data/PrimitiveArray/Checked.hs @@ -0,0 +1,29 @@ + +-- | This module exports everything that @Data.PrimitiveArray@ exports, but +-- it will do some bounds-checking on certain operations. +-- +-- Checked are: @(!)@ + +module Data.PrimitiveArray.Checked + ( module Data.PrimitiveArray + , (!) + ) where + +import qualified Data.Vector.Generic as VG + +import Data.PrimitiveArray hiding ((!)) + +-- | Bounds-checked version of indexing. +-- +-- First, we check via @inBounds@, second we check if the linear index is +-- outside of the allocated area. + +--(!) :: PrimArrayOps arr sh elm => arr sh elm -> sh -> elm +(!) arr@(Dense h v) idx + | not (inBounds (upperBound arr) idx) = error $ "(!) / inBounds: out of bounds! " ++ show (h,idx) + | li < 0 || li >= len = error $ "(!) / linearIndex: out of bounds! " ++ show (h,li,len,idx) + | otherwise = unsafeIndex arr idx + where li = linearIndex h idx + len = VG.length v +{-# Inline (!) #-} + diff --git a/Data/PrimitiveArray/Class.hs b/Data/PrimitiveArray/Class.hs new file mode 100644 --- /dev/null +++ b/Data/PrimitiveArray/Class.hs @@ -0,0 +1,285 @@ + +-- | Vastly extended primitive arrays. Some basic ideas are now modeled after the vector package, +-- especially the monadic mutable / pure immutable array system. +-- +-- Note that in general only bulk operations should error out, error handling for index/read/write +-- is too costly. General usage should be to create data structures and run the DP code within an +-- error monad, but keep error handling to high-level operations. + +module Data.PrimitiveArray.Class where + +import Control.Applicative (Applicative, pure, (<$>), (<*>)) +import Control.Exception (assert) +import Control.Monad.Except +import Control.Monad (forM_) +import Control.Monad.Primitive (PrimMonad, liftPrim) +import Control.Monad.ST (runST) +import Data.Proxy +import Data.Vector.Fusion.Util +import Debug.Trace +import GHC.Generics (Generic) +import Prelude as P +import qualified Data.Vector.Fusion.Stream.Monadic as SM +import GHC.Stack +import Data.Kind (Constraint) + +import Data.PrimitiveArray.Index.Class + + + +-- | Mutable version of an array. + +data family MutArr (m :: * -> *) (arr :: *) :: * + +-- | Associate a fill structure with each type of array (dense, sparse, ...). +-- +-- Example: @type instance FillStruc (Sparse w v sh e) = (w sh)@ associates the type @(w sh)@, which +-- is of the same type as the underlying @w@ structure holding index information for a sparse array. + +type family FillStruc arr :: * + + + +-- | The core set of operations for pure and monadic arrays. + +class (Index sh) => PrimArrayOps arr sh elm where + + -- ** Pure operations + + -- | Returns the bounds of an immutable array, again inclusive bounds: @ [lb..ub] @. + upperBound :: arr sh elm -> LimitType sh + + -- | Extract a single element from the array. Generally unsafe as not bounds-checking is + -- performed. + unsafeIndex :: arr sh elm -> sh -> elm + + -- | Index into immutable array, but safe in case @sh@ is not part of the array. + safeIndex :: arr sh elm -> sh -> Maybe elm + + -- | Savely transform the shape space of a table. + transformShape :: Index sh' => (LimitType sh -> LimitType sh') -> arr sh elm -> arr sh' elm + + -- ** Monadic operations + + -- | Return the bounds of the array. All bounds are inclusive, as in @[lb..ub]@. Technically not + -- monadic, but rather working on a monadic array. + upperBoundM :: MutArr m (arr sh elm) -> LimitType sh + + -- | Given lower and upper bounds and a list of /all/ elements, produce a mutable array. + fromListM :: PrimMonad m => LimitType sh -> [elm] -> m (MutArr m (arr sh elm)) + + -- | Creates a new array with the given bounds with each element within the array being in an + -- undefined state. + newM :: PrimMonad m => LimitType sh -> m (MutArr m (arr sh elm)) + + -- | Variant of 'newM' that requires a fill structure. Mostly for special / sparse structures + -- (hence the @S@, also to be interpreted as "safe", since these functions won't fail with sparse + -- structures). + newSM :: (Monad m, PrimMonad m) => LimitType sh -> FillStruc (arr sh elm) -> m (MutArr m (arr sh elm)) + + -- | Creates a new array with all elements being equal to 'elm'. + newWithM :: PrimMonad m => LimitType sh -> elm -> m (MutArr m (arr sh elm)) + + -- | Variant of 'newWithM' + newWithSM :: (Monad m, PrimMonad m) => LimitType sh -> FillStruc (arr sh elm) -> elm -> m (MutArr m (arr sh elm)) + + -- | Reads a single element in the array. + readM :: PrimMonad m => MutArr m (arr sh elm) -> sh -> m elm + + -- | Read from the mutable array, but return @Nothing@ in case @sh@ does not exist. This will + -- allow streaming DP combinators to "jump" over missing elements. + -- + -- Should be used with @Stream.Monadic.mapMaybe@ to get efficient code. + safeReadM :: (Monad m, PrimMonad m) => MutArr m (arr sh elm) -> sh -> m (Maybe elm) + + -- | Writes a single element in the array. + writeM :: PrimMonad m => MutArr m (arr sh elm) -> sh -> elm -> m () + + -- | Write into the mutable array, but if the index @sh@ does not exist, silently continue. + safeWriteM :: (Monad m, PrimMonad m) => MutArr m (arr sh elm) -> sh -> elm -> m () + + -- | Freezes a mutable array an returns its immutable version. This operation is /O(1)/ and both + -- arrays share the same memory. Do not use the mutable array afterwards. + unsafeFreezeM :: PrimMonad m => MutArr m (arr sh elm) -> m (arr sh elm) + + -- | Thaw an immutable array into a mutable one. Both versions share memory. + unsafeThawM :: PrimMonad m => arr sh elm -> m (MutArr m (arr sh elm)) + + +class PrimArrayMap arr sh e e' where + -- -- | Map a function of type @elm -> e@ over the primitive array, returning another primitive array + -- -- of same type and shape but different element. + mapArray :: (e -> e') -> arr sh e -> arr sh e' + + +-- | Sum type of errors that can happen when using primitive arrays. + +data PAErrors + = PAEUpperBound + deriving stock (Eq,Generic) + +instance Show PAErrors where + show (PAEUpperBound) = "Upper bound is too large for @Int@ size!" + + + +-- | Infix index operator. Performs minimal bounds-checking using assert in non-optimized code. +-- +-- @(!)@ is rewritten from phase @[1]@ onwards into an optimized form. Before, it uses a very slow +-- form, that does bounds checking. + +--(!) :: (HasCallStack, PrimArrayOps arr sh elm) => arr sh elm -> sh -> elm +(!) :: (PrimArrayOps arr sh elm) => arr sh elm -> sh -> elm +{-# Inline [1] (!) #-} +{-# Rules "unsafeIndex" [2] (!) = unsafeIndex #-} +(!) = \arr idx -> case safeIndex arr idx of + Nothing -> error $ show (showBound (upperBound arr), showIndex idx) + Just v -> v + + + +-- | Return value at an index that might not exist. + +(!?) :: PrimArrayOps arr sh elm => arr sh elm -> sh -> Maybe elm +{-# Inline (!?) #-} +(!?) = safeIndex + +-- | Returns true if the index is valid for the array. + +inBoundsM :: (Monad m, PrimArrayOps arr sh elm) => MutArr m (arr sh elm) -> sh -> Bool +inBoundsM marr idx = inBounds (upperBoundM marr) idx +{-# INLINE inBoundsM #-} + +-- -- | Given two arrays with the same dimensionality, their respective starting +-- -- index, and how many steps to go in each dimension (in terms of a dimension +-- -- again), determine if the multidimensional slices have the same value at +-- -- all positions +-- -- +-- -- TODO specialize for DIM1 (and maybe higher dim's) to use memcmp +-- +-- sliceEq :: (Eq elm, PrimArrayOps arr sh elm) => arr sh elm -> sh -> arr sh elm -> sh -> sh -> Bool +-- sliceEq arr1 k1 arr2 k2 xtnd = assert ((inBounds arr1 k1) && (inBounds arr2 k2) && (inBounds arr1 $ k1 `addDim` xtnd) && (inBounds arr2 $ k2 `addDim` xtnd)) $ and res where +-- res = zipWith (==) xs ys +-- xs = P.map (unsafeIndex arr1) $ rangeList k1 xtnd +-- ys = P.map (unsafeIndex arr2) $ rangeList k2 xtnd +-- {-# INLINE sliceEq #-} + +-- | Construct a mutable primitive array from a lower and an upper bound, a +-- default element, and a list of associations. + +fromAssocsM + :: (PrimMonad m, PrimArrayOps arr sh elm) + => LimitType sh -> elm -> [(sh,elm)] -> m (MutArr m (arr sh elm)) +fromAssocsM ub def xs = do + ma <- newWithM ub def +-- let s = size ub +-- traceShow (s,length xs) $ when (s < length xs) $ error "bang" + forM_ xs $ \(k,v) -> writeM ma k v + return ma +{-# INLINE fromAssocsM #-} + +-- | Initialize an immutable array but stay within the primitive monad @m@. + +newWithPA + :: (PrimMonad m, PrimArrayOps arr sh elm) + => LimitType sh + -> elm + -> m (arr sh elm) +newWithPA ub def = do + ma ← newWithM ub def + unsafeFreezeM ma +{-# Inlinable newWithPA #-} + +-- | Initialize an immutable array with a fill structure. + +newWithSPA + ∷ (PrimMonad m, PrimArrayOps arr sh elm) + ⇒ LimitType sh + -> FillStruc (arr sh elm) + → elm + → m (arr sh elm) +{-# Inlinable newWithSPA #-} +newWithSPA ub xs def = do + ma ← newWithSM ub xs def + unsafeFreezeM ma + +-- | Safely prepare a primitive array. +-- +-- TODO Check if having a 'MonadError' instance degrades performance. (We +-- should see this once the test with NeedlemanWunsch is under way). + +safeNewWithPA + :: forall m arr sh elm + . (PrimMonad m, MonadError PAErrors m, PrimArrayOps arr sh elm) + => LimitType sh + -> elm + -> m (arr sh elm) +safeNewWithPA ub def = do + case runExcept $ sizeIsValid maxBound [totalSize ub] of + Left (SizeError _) -> throwError PAEUpperBound + Right (CellSize _) -> newWithPA ub def +{-# Inlinable safeNewWithPA #-} + + +-- | Return all associations from an array. + +assocs :: forall arr sh elm . (IndexStream sh, PrimArrayOps arr sh elm) => arr sh elm -> [(sh,elm)] +assocs arr = unId . SM.toList $ assocsS arr +{-# INLINE assocs #-} + +-- | Return all associations from an array. + +assocsS :: forall m arr sh elm . (Monad m, IndexStream sh, PrimArrayOps arr sh elm) => arr sh elm -> SM.Stream m (sh,elm) +assocsS arr = SM.map (\k -> (k,unsafeIndex arr k)) $ streamUp zeroBound' (upperBound arr) +{-# INLINE assocsS #-} + +-- | Creates an immutable array from lower and upper bounds and a complete list +-- of elements. + +fromList :: (PrimArrayOps arr sh elm) => LimitType sh -> [elm] -> arr sh elm +fromList ub xs = runST $ fromListM ub xs >>= unsafeFreezeM +{-# INLINE fromList #-} + +-- | Creates an immutable array from lower and upper bounds, a default element, +-- and a list of associations. + +fromAssocs :: (PrimArrayOps arr sh elm) => LimitType sh -> elm -> [(sh,elm)] -> arr sh elm +fromAssocs ub def xs = runST $ fromAssocsM ub def xs >>= unsafeFreezeM +{-# INLINE fromAssocs #-} + +-- -- | Determines if an index is valid for a given immutable array. +-- +-- inBounds :: PrimArrayOps arr sh elm => arr sh elm -> sh -> Bool +-- inBounds arr idx = let (lb,ub) = bounds arr in inShapeRange lb (ub `addDim` unitDim) idx +-- {-# INLINE inBounds #-} + +-- | Returns all elements of an immutable array as a list. + +toList :: forall arr sh elm . (IndexStream sh, PrimArrayOps arr sh elm) => arr sh elm -> [elm] +toList arr = let ub = upperBound arr in P.map ((!) arr) . unId . SM.toList $ streamUp zeroBound' ub +{-# INLINE toList #-} + + + +{- + +-- * Freeze an inductive stack of tables with a 'Z' at the bottom. + +-- | 'freezeTables' freezes a stack of tables. + +class FreezeTables m t where + type Frozen t :: * + freezeTables :: t -> m (Frozen t) + +instance Applicative m => FreezeTables m Z where + type Frozen Z = Z + freezeTables Z = pure Z + {-# INLINE freezeTables #-} + +instance (Functor m, Applicative m, Monad m, PrimMonad m, FreezeTables m ts, PrimArrayOps arr sh elm) => FreezeTables m (ts:.MutArr m (arr sh elm)) where + type Frozen (ts:.MutArr m (arr sh elm)) = Frozen ts :. arr sh elm + freezeTables (ts:.t) = (:.) <$> freezeTables ts <*> unsafeFreezeM t + {-# INLINE freezeTables #-} + +-} + diff --git a/Data/PrimitiveArray/Dense.hs b/Data/PrimitiveArray/Dense.hs new file mode 100644 --- /dev/null +++ b/Data/PrimitiveArray/Dense.hs @@ -0,0 +1,189 @@ + +-- | Dense primitive arrays where the lower index is zero (or the +-- equivalent of zero for newtypes and enumerations). +-- +-- Actual @write@s to data structures use a more safe @write@ instead of +-- the unsafe @unsafeWrite@. Writes also tend to occur much less in DP +-- algorithms (say, N^2 writes for an N^3 time algorithm -- mostly reads +-- are being executed). +-- +-- TODO consider if we want to force the lower index to be zero, or allow +-- non-zero lower indices. Will have to be considered together with the +-- @Index.Class@ module! +-- +-- TODO while @Unboxed@ is, in princile, @Hashable@, we'd need the +-- corresponding @VU.Vector@ instances ... +-- +-- TODO rename to Dense.Vector, since there are other possibilities to store, +-- without basing on vector. + +module Data.PrimitiveArray.Dense where + +import Control.Lens (makeLenses) +import Control.DeepSeq +import Control.Exception (assert) +import Control.Monad (liftM, forM_, zipWithM_, when) +import Control.Monad.Primitive (PrimState) +import Data.Aeson (ToJSON,FromJSON) +import Data.Binary (Binary) +import Data.Data +import Data.Hashable (Hashable) +import Data.Serialize (Serialize) +import Data.Typeable (Typeable) +import Data.Vector.Binary +import Data.Vector.Generic.Mutable as GM hiding (length) +import Data.Vector.Serialize +import Debug.Trace +import GHC.Generics (Generic) +import qualified Data.Vector as V +import qualified Data.Vector.Fusion.Stream.Monadic as SM +import qualified Data.Vector.Generic as VG +import qualified Data.Vector.Storable as VS +import qualified Data.Vector.Unboxed as VU + +import Data.PrimitiveArray.Class +import Data.PrimitiveArray.Index.Class + + + +data Dense v sh e = Dense { _denseLimit :: !(LimitType sh), _denseV :: !(v e) } +makeLenses ''Dense + +type Unboxed sh e = Dense VU.Vector sh e + +type Storable sh e = Dense VS.Vector sh e + +type Boxed sh e = Dense V.Vector sh e + + + +deriving instance (Eq (LimitType sh), Eq (v e) ) => Eq (Dense v sh e) +deriving instance (Generic (LimitType sh), Generic (v e)) => Generic (Dense v sh e) +deriving instance (Read (LimitType sh), Read (v e) ) => Read (Dense v sh e) +deriving instance (Show (LimitType sh), Show (v e) ) => Show (Dense v sh e) +deriving instance (Functor v) => Functor (Dense v sh) + +deriving instance Typeable (Dense v sh e) + +deriving instance (Data (v e), Data (LimitType sh), Data e, Data sh, Typeable sh, Typeable e, Typeable v) => Data (Dense v sh e) + +instance (Binary (LimitType sh), Binary (v e), Generic (LimitType sh), Generic (v e)) => Binary (Dense v sh e) +instance (Serialize (LimitType sh), Serialize (v e), Generic (LimitType sh), Generic (v e)) => Serialize (Dense v sh e) +instance (ToJSON (LimitType sh), ToJSON (v e), Generic (LimitType sh), Generic (v e)) => ToJSON (Dense v sh e) +instance (FromJSON (LimitType sh), FromJSON (v e), Generic (LimitType sh), Generic (v e)) => FromJSON (Dense v sh e) +instance (Hashable (LimitType sh), Hashable (v e), Generic (LimitType sh), Generic (v e)) => Hashable (Dense v sh e) + +instance (NFData (LimitType sh), NFData (v e)) ⇒ NFData (Dense v sh e) where + rnf (Dense h xs) = rnf h `seq` rnf xs + {-# Inline rnf #-} + + + +data instance MutArr m (Dense v sh e) = MDense !(LimitType sh) !(VG.Mutable v (PrimState m) e) + deriving (Generic,Typeable) + +instance (Show (LimitType sh), Show (VG.Mutable v (PrimState m) e), VG.Mutable v (PrimState m) e ~ mv) ⇒ Show (MutArr m (Dense v sh e)) where + show (MDense sh mv) = show (sh,mv) + +instance (NFData (LimitType sh), NFData (VG.Mutable v (PrimState m) e), VG.Mutable v (PrimState m) e ~ mv) ⇒ NFData (MutArr m (Dense v sh e)) where + rnf (MDense h xs) = rnf h `seq` rnf xs + {-# Inline rnf #-} + +{- +instance + ( Index sh, MutArr m (Dense v sh e) ~ mv + , GM.MVector (VG.Mutable v) e +#if ADPFUSION_DEBUGOUTPUT + , Show sh, Show (LimitType sh), Show e +#endif + ) ⇒ MPrimArrayOps (Dense v) sh e where +-} + +instance + ( Index sh, VG.Vector v e +#if ADPFUSION_DEBUGOUTPUT + , Show sh, Show (LimitType sh), Show e +#endif + ) ⇒ PrimArrayOps (Dense v) sh e where + + -- ** pure operations + + {-# Inline upperBound #-} + upperBound (Dense h _) = h + {-# Inline unsafeFreezeM #-} + unsafeFreezeM (MDense h mba) = Dense h `liftM` VG.unsafeFreeze mba + {-# Inline unsafeThawM #-} + unsafeThawM (Dense h ba) = MDense h `liftM` VG.unsafeThaw ba + {-# Inline unsafeIndex #-} + unsafeIndex (Dense h ba) idx = VG.unsafeIndex ba (linearIndex h idx) + {-# Inline safeIndex #-} + safeIndex (Dense h ba) idx = if inBounds h idx then Just $ unsafeIndex (Dense h ba) idx else Nothing + {-# Inline transformShape #-} + transformShape tr (Dense h ba) = Dense (tr h) ba + + -- ** monadic operations + + {-# Inline upperBoundM #-} + upperBoundM (MDense h _) = h + {-# Inline fromListM #-} + fromListM h xs = do + ma ← newM h + let (MDense _ mba) = ma + -- there need to be at least as many elements, as we want to fill. There could be more, in debug + -- tests, we like to do @[0..]@ and this should not trigger the assert. + SM.zipWithM_ (\k x → assert (length (Prelude.take (size h) xs) == size h) $ unsafeWrite mba k x) (SM.enumFromTo 0 (size h -1)) (SM.fromList xs) + return ma + {-# Inline newM #-} -- TODO was NoInline, check if anything breaks! + newM h = MDense h `liftM` new (size h) + {-# Inline newSM #-} + newSM = error "not implemented, use newM for dense arrays" + {-# Inline newWithM #-} + newWithM h def = do + ma ← newM h + let (MDense _ mba) = ma + GM.set mba def + return ma + {-# Inline newWithSM #-} + newWithSM = error "not implemented, use newWithSM for dense arrays" + {-# Inline readM #-} + readM (MDense h mba) idx = assert (inBounds h idx) $ unsafeRead mba (linearIndex h idx) + {-# Inline safeReadM #-} + safeReadM dense idx = if inBoundsM dense idx then Just <$> readM dense idx else undefined + {-# Inline writeM #-} + writeM (MDense h mba) idx elm = +#if ADPFUSION_DEBUGOUTPUT + (if inBounds h idx then id else traceShow ("writeM", h, idx, elm, size h, linearIndex h idx, inBounds h idx)) +#endif + assert (inBounds h idx) $ unsafeWrite mba (linearIndex h idx) elm + {-# Inline safeWriteM #-} + safeWriteM dense idx elm = when (inBoundsM dense idx) $ writeM dense idx elm + +instance (Index sh, VG.Vector v e, VG.Vector v e') ⇒ PrimArrayMap (Dense v) sh e e' where + {-# Inline mapArray #-} + mapArray f (Dense h xs) = Dense h (VG.map f xs) + + +{- + - + - This stuff tells us how to write efficient generics on large data + - constructors like the Turner and Vienna ctors. + - + +import qualified Data.Generics.TH as T + +data Unboxed sh e = Unboxed !sh !(VU.Vector e) + deriving (Show,Eq,Ord) + +data X e = X (Unboxed DIM1 e) (Unboxed DIM1 e) + deriving (Show,Eq,Ord) + +x :: X Int +x = X z z where z = (Unboxed (Z:.10) (VU.fromList [ 0 .. 10] )) + +pot :: X Int -> X Double +pot = $( T.thmapT (T.mkTs ['f]) [t| X Int |] ) where + f :: Unboxed DIM1 Int -> Unboxed DIM1 Double + f (Unboxed sh xs) = Unboxed sh (VU.map fromIntegral xs) + +-} + diff --git a/Data/PrimitiveArray/HashTable.hs b/Data/PrimitiveArray/HashTable.hs new file mode 100644 --- /dev/null +++ b/Data/PrimitiveArray/HashTable.hs @@ -0,0 +1,64 @@ + +-- | A table representation that internally uses a hashtable from the @hashtables@ library. The +-- implementation is currently a testbed on which idea makes the most sense. +-- +-- In particular, once a hashtable has been created with, say, @newWithPA@, it will be completely +-- void of any entries. To prime the system, call @setValidKeys@ which will setup all keys that are +-- vaild, as well as setup an additional data structure to help with @streamUp@ and @streamDown@. +-- +-- This table does not store default values, since it is assumed that lookups are only done on valid +-- keys, and @ADPfusion@ as the default consumer should have rules "jump over" missing keys. +-- +-- Currently the idea is that any write to an undeclared key will just fail SILENTLY! +-- +-- TODO this also forces rethinking @inBounds@, as this will now depend on the internal structure +-- given via @setValidKeys@. + +module Data.PrimitiveArray.HashTable where + +import Control.Monad.Primitive +import Control.Monad.ST +import Control.Monad.ST.Unsafe +import Data.HashTable.Class as HT +import Data.HashTable.IO as HTIO +import Unsafe.Coerce + +import Data.PrimitiveArray.Class +import Data.PrimitiveArray.Index.Class + + + +data Hashed ht sh e = Hashed + { _hashedUpperBound :: !(LimitType sh) + -- ^ Explicitly store the upper bound. + , _hashedTable :: !(IOHashTable ht sh e) + -- ^ The hashtable to be updated / used. + , _hashedUpDown :: !() + -- ^ Helper structure for the @streamUp@ / @streamDown@ functionality. + -- + -- TODO this should be a recursively constructed hashtable, based on the shape of @sh@. + } + + + +-- | Sets valid keys, working within a primitive Monad. The valid keys should be a hashtable with +-- all correct keys, but values set to something resembling a default value. A good choice will be +-- akin to @mzero@. +-- +-- Internally, this function uses @unsafeCoerce@ to change the @PrimState@ token held by the hash +-- table to @RealWord@, from whatever it is. +-- +-- TODO setup the @hashedUpDown@ part, once it is clear what to do. + +setValidKeys + :: (PrimMonad m, HashTable h) + => LimitType sh + -> h (PrimState m) k v + -> m (Hashed ht sh e) +{-# Inline setValidKeys #-} +setValidKeys ub ks = return $ Hashed + { _hashedUpperBound = ub + , _hashedTable = unsafeCoerce ks + , _hashedUpDown = () + } + diff --git a/Data/PrimitiveArray/Index.hs b/Data/PrimitiveArray/Index.hs new file mode 100644 --- /dev/null +++ b/Data/PrimitiveArray/Index.hs @@ -0,0 +1,29 @@ + +module Data.PrimitiveArray.Index + ( module Data.PrimitiveArray.Index.Class + , module Data.PrimitiveArray.Index.BitSet0 + , module Data.PrimitiveArray.Index.BitSet1 + , module Data.PrimitiveArray.Index.BitSetClasses +-- , module Data.PrimitiveArray.Index.EdgeBoundary + , module Data.PrimitiveArray.Index.Int + , module Data.PrimitiveArray.Index.IOC + , module Data.PrimitiveArray.Index.PhantomInt + , module Data.PrimitiveArray.Index.Point +-- , module Data.PrimitiveArray.Index.Set + , module Data.PrimitiveArray.Index.Subword + , module Data.PrimitiveArray.Index.Unit + ) where + +import Data.PrimitiveArray.Index.Class +--import Data.PrimitiveArray.Index.EdgeBoundary hiding (streamUpMk, streamUpStep, streamDownMk, streamDownStep) +import Data.PrimitiveArray.Index.Int +import Data.PrimitiveArray.Index.IOC +import Data.PrimitiveArray.Index.PhantomInt hiding (streamUpMk, streamUpStep, streamDownMk, streamDownStep) +import Data.PrimitiveArray.Index.Point hiding (streamUpMk, streamUpStep, streamDownMk, streamDownStep) +--import Data.PrimitiveArray.Index.Set hiding (streamUpBsMk, streamUpBsStep, streamDownBsMk, StreamDownBsStep, streamUpBsIMk, streamUpBsIStep, streamDownBsIMk, StreamDownBsIStep, streamUpBsIiMk, streamUpBsIiStep, streamDownBsIiMk, StreamDownBsIiStep) +import Data.PrimitiveArray.Index.BitSet1 hiding (streamUpMk, streamUpStep, streamDownMk, streamDownStep) +import Data.PrimitiveArray.Index.BitSet0 hiding (streamUpMk, streamUpStep, streamDownMk, streamDownStep) +import Data.PrimitiveArray.Index.BitSetClasses +import Data.PrimitiveArray.Index.Subword hiding (streamUpMk, streamUpStep, streamDownMk, streamDownStep) +import Data.PrimitiveArray.Index.Unit + diff --git a/Data/PrimitiveArray/Index/BitSet0.hs b/Data/PrimitiveArray/Index/BitSet0.hs new file mode 100644 --- /dev/null +++ b/Data/PrimitiveArray/Index/BitSet0.hs @@ -0,0 +1,144 @@ + +-- | The most basic bitset structure. Alone, not particularly useful, because +-- two sets @{u,v},{v',w}@ have no way of annotating the connection between the +-- sets. Together with boundaries this yields sets for useful DP algorithms. + +module Data.PrimitiveArray.Index.BitSet0 where + +import Control.DeepSeq (NFData(..)) +import Control.Lens (makeLenses) +import Data.Aeson (FromJSON,ToJSON,FromJSONKey,ToJSONKey) +import Data.Binary (Binary) +import Data.Bits +import Data.Bits.Extras +import Data.Hashable (Hashable) +import Data.Serialize (Serialize) +import Data.Vector.Unboxed.Deriving +import Data.Vector.Unboxed (Unbox(..)) +import Debug.Trace +import GHC.Generics (Generic) +import qualified Data.Vector.Fusion.Stream.Monadic as SM +import Test.QuickCheck + +import Data.Bits.Ordered +import Data.PrimitiveArray.Index.Class +import Data.PrimitiveArray.Index.IOC +import Data.PrimitiveArray.Index.BitSetClasses + + + +-- | Newtype for a bitset. +-- +-- @Int@ integrates better with the rest of the framework. But we should +-- consider moving to @Word@-based indexing, if possible. + +newtype BitSet t = BitSet { _bitSet :: Int } + deriving stock (Eq,Ord,Generic) + deriving newtype (FiniteBits,Ranked,Num,Bits) +makeLenses ''BitSet + +instance FromJSON (BitSet t) +instance FromJSONKey (BitSet t) +instance ToJSON (BitSet t) +instance ToJSONKey (BitSet t) +instance Binary (BitSet t) +instance Serialize (BitSet t) +instance Hashable (BitSet t) + +derivingUnbox "BitSet" + [t| forall t . BitSet t → Int |] + [| \(BitSet s) → s |] + [| BitSet |] + +instance Show (BitSet t) where + show (BitSet s) = "<" ++ (show $ activeBitsL s) ++ ">(" ++ show s ++ ")" + +instance NFData (BitSet t) where + rnf (BitSet s) = rnf s + {-# Inline rnf #-} + +instance Index (BitSet t) where + newtype LimitType (BitSet t) = LtBitSet Int + linearIndex _ (BitSet z) = z + {-# Inline linearIndex #-} + size (LtBitSet pc) = 2 ^ pc -- 2 ^ popCount h - 2 ^ popCount l + 1 + {-# Inline size #-} + inBounds (LtBitSet h) z = popCount z <= h + {-# Inline inBounds #-} + zeroBound = BitSet 0 + {-# Inline zeroBound #-} + zeroBound' = LtBitSet 0 + {-# Inline zeroBound' #-} + totalSize (LtBitSet n) = [2 ^ fromIntegral n] + {-# Inline totalSize #-} + fromLinearIndex _ = BitSet + {-# Inline [0] fromLinearIndex #-} + showBound (LtBitSet b) = ["LtBitSet " ++ show b] + showIndex (BitSet b) = ["BitSet " ++ show b] + +instance SetPredSucc (BitSet t) where + setSucc l h s + | cs > ch = Nothing + | Just s' <- popPermutation ch s = Just s' + | cs >= ch = Nothing + | cs < ch = Just . BitSet $ 2^(cs+1) -1 + where ch = popCount h + cs = popCount s + {-# Inline setSucc #-} + setPred l h s + | cs < cl = Nothing + | Just s' <- popPermutation ch s = Just s' + | cs <= cl = Nothing + | cs > cl = Just . BitSet $ 2^(cs-1) -1 + where cl = popCount l + ch = popCount h + cs = popCount s + {-# Inline setPred #-} + +instance IndexStream z => IndexStream (z:.BitSet I) where + streamUp (ls:..LtBitSet l) (hs:..LtBitSet h) = SM.flatten (streamUpMk l h) (streamUpStep l h) $ streamUp ls hs + streamDown (ls:..LtBitSet l) (hs:..LtBitSet h) = SM.flatten (streamDownMk l h) (streamDownStep l h) $ streamDown ls hs + {-# Inline streamUp #-} + {-# Inline streamDown #-} + +instance IndexStream z ⇒ IndexStream (z:.BitSet O) where + streamUp (ls:..LtBitSet l) (hs:..LtBitSet h) = SM.flatten (streamDownMk l h) (streamDownStep l h) $ streamUp ls hs + streamDown (ls:..LtBitSet l) (hs:..LtBitSet h) = SM.flatten (streamUpMk l h) (streamUpStep l h) $ streamDown ls hs + {-# Inline streamUp #-} + {-# Inline streamDown #-} + +instance IndexStream z ⇒ IndexStream (z:.BitSet C) where + streamUp (ls:..LtBitSet l) (hs:..LtBitSet h) = SM.flatten (streamUpMk l h) (streamUpStep l h) $ streamUp ls hs + streamDown (ls:..LtBitSet l) (hs:..LtBitSet h) = SM.flatten (streamDownMk l h) (streamDownStep l h) $ streamDown ls hs + {-# Inline streamUp #-} + {-# Inline streamDown #-} + +instance IndexStream (Z:.BitSet t) ⇒ IndexStream (BitSet t) where + streamUp l h = SM.map (\(Z:.i) -> i) $ streamUp (ZZ:..l) (ZZ:..h) + {-# Inline streamUp #-} + streamDown l h = SM.map (\(Z:.i) -> i) $ streamDown (ZZ:..l) (ZZ:..h) + {-# Inline streamDown #-} + +streamUpMk ∷ Monad m ⇒ Int → Int → t → m (t, Maybe (BitSet ioc)) +streamUpMk l h z = return (z, if l <= h then Just (BitSet $ 2^l-1) else Nothing) +{-# Inline [0] streamUpMk #-} + +streamUpStep ∷ Monad m ⇒ Int → Int → (t, Maybe (BitSet ioc)) → m (SM.Step (t, Maybe (BitSet ioc)) (t:.BitSet ioc)) +streamUpStep l h (z , Nothing) = return $ SM.Done +streamUpStep l h (z , Just t ) = return $ SM.Yield (z:.t) (z, setSucc (2^l-1) (2^h-1) t) +{-# Inline [0] streamUpStep #-} + +streamDownMk ∷ Monad m ⇒ Int → Int → t → m (t, Maybe (BitSet ioc)) +streamDownMk l h z = return (z, if l <=h then Just (BitSet $ 2^l-1) else Nothing) +{-# Inline [0] streamDownMk #-} + +streamDownStep ∷ Monad m ⇒ Int → Int → (t, Maybe (BitSet ioc)) → m (SM.Step (t, Maybe (BitSet ioc)) (t:.BitSet ioc)) +streamDownStep l h (z , Nothing) = return $ SM.Done +streamDownStep l h (z , Just t ) = return $ SM.Yield (z:.t) (z , setPred (2^l-1) (2^h-1) t) +{-# Inline [0] streamDownStep #-} + +instance Arbitrary (BitSet t) where + arbitrary = BitSet <$> choose (0,2^arbitraryBitSetMax-1) + shrink s = let s' = [ s `clearBit` a | a <- activeBitsL s ] + in s' ++ concatMap shrink s' + diff --git a/Data/PrimitiveArray/Index/BitSet1.hs b/Data/PrimitiveArray/Index/BitSet1.hs new file mode 100644 --- /dev/null +++ b/Data/PrimitiveArray/Index/BitSet1.hs @@ -0,0 +1,172 @@ + +-- | A bitset with one interface. This includes the often-encountered case +-- where @{u,v},{v}@, or sets with a single edge between the old set and a new +-- singleton set are required. Uses are Hamiltonian path problems, and TSP, +-- among others. + +module Data.PrimitiveArray.Index.BitSet1 where + +import Control.DeepSeq (NFData(..)) +import Control.Lens (makeLenses) +import Control.Monad.Except +import Data.Aeson (FromJSON,ToJSON,FromJSONKey,ToJSONKey) +import Data.Binary (Binary) +import Data.Bits +import Data.Bits.Extras +import Data.Hashable (Hashable) +import Data.Serialize (Serialize) +import Data.Vector.Unboxed.Deriving +import Data.Vector.Unboxed (Unbox(..)) +import Debug.Trace +import GHC.Generics (Generic) +import qualified Data.Vector.Fusion.Stream.Monadic as SM +import Test.QuickCheck + +import Data.Bits.Ordered +import Data.PrimitiveArray.Index.BitSet0 (BitSet(..),LimitType(..)) +import Data.PrimitiveArray.Index.BitSetClasses +import Data.PrimitiveArray.Index.Class +import Data.PrimitiveArray.Index.IOC + + + +-- | The bitset with one interface or boundary. + +data BitSet1 i ioc = BitSet1 { _bitset ∷ !(BitSet ioc), _boundary ∷ !(Boundary i ioc) } + deriving (Eq,Ord,Generic,Show) +makeLenses ''BitSet1 + +derivingUnbox "BitSet1" + [t| forall i ioc . BitSet1 i ioc → (Int,Int) |] + [| \ (BitSet1 (BitSet set) (Boundary bnd)) → (set,bnd) |] + [| \ (set,bnd) → BitSet1 (BitSet set) (Boundary bnd) |] + + +-- | +-- +-- NOTE We linearize a bitset as follows: we need @2^number-of-bits * +-- number-of-bits@ elements. The first is due to having a binary set structure. +-- The second is due to pointing to each of those elements as being the +-- boundary. This overcommits on memory since only those bits can be a boundary +-- bits that are actually set. Furthermore, in case no bit is set at all, then +-- there should be no boundary. This is currently rather awkwardly done by +-- restricting enumeration and mapping the 0-set to boundary 0. +-- +-- | TODO The size calculations are off by a factor of two, exactly. Each +-- bitset (say) @00110@ has a mirror image @11001@, whose elements do not have +-- to be indexed. It has to be investigated if a version with exact memory +-- bounds is slower in indexing. + +instance Index (BitSet1 bnd ioc) where + -- This is the number of bits. Meaning that @LtNumBits1 3@ yields @[0,1,2]@. + -- TODO Should we rename this to @NumberOfBits1@? Or have a newtype @NumBits@? + newtype LimitType (BitSet1 bnd ioc) = LtNumBits1 Int + -- Calculate the linear index for a set. Spread out by the possible number of + -- bits to fit the actual boundary results. Add the boundary index. + linearIndex (LtNumBits1 pc) (BitSet1 set (Boundary bnd)) + = linearIndex (LtBitSet pc) set * pc + bnd + {-# Inline linearIndex #-} + size (LtNumBits1 pc) = 2^pc * pc + 1 + {-# Inline size #-} + inBounds (LtNumBits1 pc) (BitSet1 set bnd) = popCount set <= pc && 0 <= bnd && getBoundary bnd <= pc + {-# Inline inBounds #-} + zeroBound = BitSet1 zeroBound zeroBound + {-# Inline zeroBound #-} + zeroBound' = LtNumBits1 0 + {-# Inline zeroBound' #-} + totalSize (LtNumBits1 pc) = + let z = fromIntegral pc + in [z * 2 ^ z] + fromLinearIndex (LtNumBits1 pc) z = error "implement me" + showBound = error "implement me" + showIndex = error "implement me" + +deriving instance Show (LimitType (BitSet1 bnd ioc)) + +instance IndexStream z ⇒ IndexStream (z:.BitSet1 i I) where + streamUp (ls:..LtNumBits1 l) (hs:..LtNumBits1 h) = SM.flatten (streamUpMk l h) (streamUpStep l h) $ streamUp ls hs + streamDown (ls:..LtNumBits1 l) (hs:..LtNumBits1 h) = SM.flatten (streamDownMk l h) (streamDownStep l h) $ streamDown ls hs + {-# Inline streamUp #-} + {-# Inline streamDown #-} + +instance IndexStream z ⇒ IndexStream (z:.BitSet1 i O) where + streamUp (ls:..LtNumBits1 l) (hs:..LtNumBits1 h) = SM.flatten (streamDownMk l h) (streamDownStep l h) $ streamUp ls hs + streamDown (ls:..LtNumBits1 l) (hs:..LtNumBits1 h) = SM.flatten (streamUpMk l h) (streamUpStep l h) $ streamDown ls hs + {-# Inline streamUp #-} + {-# Inline streamDown #-} + +--instance IndexStream z => IndexStream (z:.BS1 i C) where +-- streamUp (ls:..l) (hs:..h) = flatten (streamUpBsIMk l h) (streamUpBsIStep l h) $ streamUp ls hs +-- streamDown (ls:..l) (hs:..h) = flatten (streamDownBsIMk l h) (streamDownBsIStep l h) $ streamDown ls hs +-- {-# Inline streamUp #-} +-- {-# Inline streamDown #-} + +instance IndexStream (Z:.BitSet1 i t) ⇒ IndexStream (BitSet1 i t) where + streamUp l h = SM.map (\(Z:.i) -> i) $ streamUp (ZZ:..l) (ZZ:..h) + {-# Inline streamUp #-} + streamDown l h = SM.map (\(Z:.i) -> i) $ streamDown (ZZ:..l) (ZZ:..h) + {-# Inline streamDown #-} + +streamUpMk ∷ Monad m ⇒ Int → Int → z → m (z, Maybe (BitSet1 c ioc)) +streamUpMk l h z = + let set = BitSet $ 2^l-1 + -- lsbZ set == 0, or no active bits in which case we use 0 + bnd = UndefBoundary + in return (z, if l <= h then Just (BitSet1 set bnd) else Nothing) +{-# Inline [0] streamUpMk #-} + +streamUpStep ∷ Monad m ⇒ Int → Int → (t, Maybe (BitSet1 c ioc)) → m (SM.Step (t, Maybe (BitSet1 c ioc)) (t:.BitSet1 c ioc)) +streamUpStep l h (z, Nothing) = return $ SM.Done +streamUpStep l h (z, Just t ) = return $ SM.Yield (z:.t) (z , setSucc l h t) +{-# Inline [0] streamUpStep #-} + +streamDownMk ∷ Monad m ⇒ Int → Int → z → m (z, Maybe (BitSet1 c ioc)) +streamDownMk l h z = + let set = BitSet $ 2^h-1 + bnd = Boundary 0 -- this is the actual boundary at zero + in return (z, if l <= h then Just (BitSet1 set bnd) else Nothing) +{-# Inline [0] streamDownMk #-} + +streamDownStep ∷ Monad m ⇒ Int → Int → (t, Maybe (BitSet1 c ioc)) → m (SM.Step (t, Maybe (BitSet1 c ioc)) (t:.BitSet1 c ioc)) +streamDownStep l h (z, Nothing) = return $ SM.Done +streamDownStep l h (z, Just t ) = return $ SM.Yield (z:.t) (z , setPred l h t) +{-# Inline [0] streamDownStep #-} + +instance SetPredSucc (BitSet1 t ioc) where + setSucc pcl pch (BitSet1 s (Boundary is)) + | cs > pch = Nothing + | Just is' <- maybeNextActive is s = Just $ BitSet1 s (Boundary is') + | Just s' <- popPermutation pch s = Just $ BitSet1 s' (Boundary $ lsbZ s') + | cs >= pch = Nothing + | cs < pch = let s' = BitSet $ 2^(cs+1)-1 + in Just (BitSet1 s' (Boundary (lsbZ s'))) + where cs = popCount s + {-# Inline setSucc #-} + setPred pcl pch (BitSet1 s (Boundary is)) + | cs < pcl = Nothing + | Just is' <- maybeNextActive is s = Just $ BitSet1 s (Boundary is') + | Just s' <- popPermutation pch s = Just $ BitSet1 s' (Boundary $ lsbZ s') + | cs <= pcl = Nothing + | cs > pcl = let s' = BitSet $ 2^(cs-1)-1 + in Just (BitSet1 s' (Boundary (max 0 $ lsbZ s'))) + where cs = popCount s + {-# Inline setPred #-} + +instance SetPredSucc (FixedMask (BitSet1 t ioc)) where + setPred = error "implement me" + setSucc pcl pch (FixedMask mask bs1) = undefined + +instance Arbitrary (BitSet1 t ioc) where + arbitrary = do + s <- arbitrary + if s==0 + then return (BitSet1 s 0) + else do i <- elements $ activeBitsL s + return (BitSet1 s $ Boundary i) + shrink (BitSet1 s i) = + let s' = [ BitSet1 (s `clearBit` a) i + | a <- activeBitsL s + , Boundary a /= i ] + ++ [ BitSet1 0 0 | popCount s == 1 ] + in s' ++ concatMap shrink s' + diff --git a/Data/PrimitiveArray/Index/BitSetClasses.hs b/Data/PrimitiveArray/Index/BitSetClasses.hs new file mode 100644 --- /dev/null +++ b/Data/PrimitiveArray/Index/BitSetClasses.hs @@ -0,0 +1,175 @@ + +-- | A collection of a number of data types and type classes shared by all +-- bitset variants. + +module Data.PrimitiveArray.Index.BitSetClasses where + +import Control.DeepSeq (NFData(..)) +import Data.Aeson (FromJSON,ToJSON,FromJSONKey,ToJSONKey) +import Data.Binary (Binary) +import Data.Hashable (Hashable) +import Data.Serialize (Serialize) +import Data.Vector.Unboxed.Deriving +import GHC.Generics (Generic) +import qualified Data.Vector.Fusion.Stream.Monadic as SM +import qualified Data.Vector.Unboxed as VU + +import Data.Bits.Ordered +import Data.PrimitiveArray.Index.Class +import Data.PrimitiveArray.Index.IOC + + + +-- * Boundaries, the interface(s) for bitsets. + +-- | Certain sets have an interface, a particular element with special +-- meaning. In this module, certain ``meanings'' are already provided. +-- These include a @First@ element and a @Last@ element. We phantom-type +-- these to reduce programming overhead. + +newtype Boundary boundaryType ioc = Boundary { getBoundary ∷ Int } + deriving stock (Eq,Ord,Generic) + deriving newtype (Num) + +-- | Whenever we can not set the boundary we should have for a set, we use this +-- pattern. All legal boundaries are @>=0@. We also need to set the undefined +-- boundary to @0@, since the @linearIndex@ will use this value to add, which +-- for empty sets would reduce to @0 - UndefBoundary === 0@. + +pattern UndefBoundary ∷ Boundary boundaryType ioc +pattern UndefBoundary = Boundary 0 + +instance Show (Boundary i t) where + show (Boundary i) = "(I:" ++ show i ++ ")" + +derivingUnbox "Boundary" + [t| forall i t . Boundary i t → Int |] + [| \(Boundary i) → i |] + [| Boundary |] + +instance Binary (Boundary i t) +instance Serialize (Boundary i t) +instance ToJSON (Boundary i t) +instance FromJSON (Boundary i t) +instance Hashable (Boundary i t) + +instance NFData (Boundary i t) where + rnf (Boundary i) = rnf i + {-# Inline rnf #-} + +instance Index (Boundary i t) where + newtype LimitType (Boundary i t) = LtBoundary Int + linearIndex _ (Boundary z) = z + {-# INLINE linearIndex #-} + size (LtBoundary h) = h + 1 + {-# INLINE size #-} + inBounds (LtBoundary h) z = 0 <= z && getBoundary z <= h + {-# INLINE inBounds #-} + zeroBound = Boundary 0 + {-# Inline zeroBound #-} + zeroBound' = LtBoundary 0 + {-# Inline zeroBound' #-} + totalSize (LtBoundary n) = [fromIntegral n] + {-# Inline totalSize #-} + fromLinearIndex _ = Boundary + {-# Inline fromLinearIndex #-} + showBound (LtBoundary b) = ["LtBoundary " ++ show b] + showIndex (Boundary b) = ["Boundary " ++ show b] + +instance IndexStream z ⇒ IndexStream (z:.Boundary k I) where + streamUp (ls:..LtBoundary l) (hs:..LtBoundary h) = SM.flatten (streamUpBndMk l h) (streamUpBndStep l h) $ streamUp ls hs + streamDown (ls:..LtBoundary l) (hs:..LtBoundary h) = SM.flatten (streamDownBndMk l h) (streamDownBndStep l h) $ streamDown ls hs + {-# Inline streamUp #-} + {-# Inline streamDown #-} + +instance IndexStream (Z:.Boundary k I) ⇒ IndexStream (Boundary k I) where + streamUp l h = SM.map (\(Z:.i) -> i) $ streamUp (ZZ:..l) (ZZ:..h) + {-# Inline streamUp #-} + streamDown l h = SM.map (\(Z:.i) -> i) $ streamDown (ZZ:..l) (ZZ:..h) + {-# Inline streamDown #-} + +streamUpBndMk l h z = return (z, l) +{-# Inline [0] streamUpBndMk #-} + +streamUpBndStep l h (z , k) + | k > h = return $ SM.Done + | otherwise = return $ SM.Yield (z:.Boundary k) (z, k+1) +{-# Inline [0] streamUpBndStep #-} + +streamDownBndMk l h z = return (z, h) +{-# Inline [0] streamDownBndMk #-} + +streamDownBndStep l h (z , k) + | k < l = return $ SM.Done + | otherwise = return $ SM.Yield (z:.Boundary k) (z,k-1) +{-# Inline [0] streamDownBndStep #-} + +-- | Declare the interface to be the start of a path. + +data First + +-- | Declare the interface to be the end of a path. + +data Last + +-- | Declare the interface to match anything. +-- +-- TODO needed? want to use later in ADPfusion + +data Any + + + +-- * Moving indices within sets. + +-- | Successor and Predecessor for sets. Designed as a class to accomodate +-- sets with interfaces and without interfaces with one function. +-- +-- The functions are not written recursively, as we currently only have +-- three cases, and we do not want to "reset" while generating successors +-- and predecessors. +-- +-- Note that sets have a partial order. Within the group of element with +-- the same @popCount@, we use @popPermutation@ which has the same stepping +-- order for both, @setSucc@ and @setPred@. + +class SetPredSucc s where + -- | Set successor. The first argument is the lower set limit, the second + -- the upper set limit, the third the current set. + setSucc ∷ Int → Int → s → Maybe s + -- | Set predecessor. The first argument is the lower set limit, the + -- second the upper set limit, the third the current set. + setPred ∷ Int → Int → s → Maybe s + +-- | Masks are used quite often for different types of bitsets. We liberate +-- them as a type family. + +type family Mask s ∷ * + +-- | @Fixed@ allows us to fix some or all bits of a bitset, thereby +-- providing @succ/pred@ operations which are only partially free. +-- +-- @f = getFixedMask .&. getFixed@ are the fixed bits. +-- @n = getFixed .&. complement getFixedMask@ are the free bits. +-- @to = complement getFixed@ is the to move mask +-- @n' = popShiftR to n@ yields the population after the move +-- @p = popPermutation undefined n'@ yields the new population permutation +-- @p' = popShiftL to p@ yields the population moved back +-- @final = p' .|. f@ + +data FixedMask t = FixedMask { getMask ∷ (Mask t) , getFixed ∷ !t } + +-- | Assuming a bitset on bits @[0 .. highbit]@, we can apply a mask that +-- stretches out those bits over @[0 .. higherBit]@ with @highbit <= +-- higherBit@. Any active interfaces are correctly set as well. + +class ApplyMask s where + applyMask :: Mask s → s → s + + + +-- | for 'Test.QuickCheck.Arbitrary' + +arbitraryBitSetMax ∷ Int +arbitraryBitSetMax = 6 + diff --git a/Data/PrimitiveArray/Index/Class.hs b/Data/PrimitiveArray/Index/Class.hs new file mode 100644 --- /dev/null +++ b/Data/PrimitiveArray/Index/Class.hs @@ -0,0 +1,344 @@ + +module Data.PrimitiveArray.Index.Class where + +import Control.Applicative +import Control.DeepSeq (NFData(..)) +import Control.Lens hiding (Index, (:>)) +import Control.Monad.Except +import Control.Monad (liftM2) +import Data.Aeson +import Data.Binary +import Data.Data +import Data.Hashable (Hashable) +import Data.Proxy +import Data.Serialize +import Data.Typeable +import Data.Vector.Fusion.Stream.Monadic (Stream) +import Data.Vector.Unboxed.Deriving +import Data.Vector.Unboxed (Unbox(..)) +import GHC.Base (quotRemInt) +import GHC.Generics +import GHC.TypeNats +import qualified Data.Vector.Fusion.Stream.Monadic as SM +import Test.QuickCheck +import Text.Printf +import Data.Type.Equality + + + +infixl 3 :. + +-- | Strict pairs -- as in @repa@. + +data a :. b = !a :. !b + deriving (Eq,Ord,Show,Generic,Data,Typeable) + +derivingUnbox "StrictPair" + [t| forall a b . (Unbox a, Unbox b) => (a:.b) -> (a,b) |] + [| \(a:.b) -> (a, b) |] + [| \(a,b) -> (a:.b) |] + +instance (Binary a, Binary b) => Binary (a:.b) +instance (Serialize a, Serialize b) => Serialize (a:.b) +instance (ToJSON a, ToJSON b) => ToJSON (a:.b) +instance (FromJSON a, FromJSON b) => FromJSON (a:.b) +instance (Hashable a, Hashable b) => Hashable (a:.b) + +instance (ToJSON a , ToJSONKey a, ToJSON b , ToJSONKey b) => ToJSONKey (a:.b) +instance (FromJSON a, FromJSONKey a, FromJSON b, FromJSONKey b) => FromJSONKey (a:.b) + +deriving instance (Read a, Read b) => Read (a:.b) + +instance (NFData a, NFData b) => NFData (a:.b) where + rnf (a:.b) = rnf a `seq` rnf b + {-# Inline rnf #-} + +instance (Arbitrary a, Arbitrary b) => Arbitrary (a :. b) where + arbitrary = liftM2 (:.) arbitrary arbitrary + shrink (a:.b) = [ (a':.b) | a' <- shrink a ] ++ [ (a:.b') | b' <- shrink b ] + +infixr 3 :> + +-- | A different version of strict pairs. Makes for simpler type inference in +-- multi-tape grammars. We use @:>@ when we have special needs, like +-- non-recursive instances on inductives tuples, as used for set indices. +-- +-- This one is @infixr@ so that in @a :> b@ we can have the main type in +-- @a@ and the specializing types in @b@ and then dispatch on @a :> ts@ +-- with @ts@ maybe a chain of @:>@. + +data a :> b = !a :> !b + deriving (Eq,Ord,Show,Generic,Data,Typeable) + +derivingUnbox "StrictIxPair" + [t| forall a b . (Unbox a, Unbox b) => (a:>b) -> (a,b) |] + [| \(a:>b) -> (a, b) |] + [| \(a,b) -> (a:>b) |] + +instance (Binary a, Binary b) => Binary (a:>b) +instance (Serialize a, Serialize b) => Serialize (a:>b) +instance (ToJSON a, ToJSON b) => ToJSON (a:>b) +instance (FromJSON a, FromJSON b) => FromJSON (a:>b) +instance (Hashable a, Hashable b) => Hashable (a:>b) + +deriving instance (Read a, Read b) => Read (a:>b) + +instance (NFData a, NFData b) => NFData (a:>b) where + rnf (a:>b) = rnf a `seq` rnf b + {-# Inline rnf #-} + +--instance (Arbitrary a, Arbitrary b) => Arbitrary (a :> b) where +-- arbitrary = (:>) <$> arbitrary <*> arbitrary +-- shrink (a:>b) = (:>) <$> shrink a <*> shrink b + + + +-- | Base data constructor for multi-dimensional indices. + +data Z = Z + deriving (Eq,Ord,Read,Show,Generic,Data,Typeable,Bounded) + +derivingUnbox "Z" + [t| Z -> () |] + [| const () |] + [| const Z |] + +instance Binary Z +instance Serialize Z +instance ToJSON Z +instance FromJSON Z +instance Hashable Z + +instance Arbitrary Z where + arbitrary = return Z + +instance NFData Z where + rnf Z = () + {-# Inline rnf #-} + + + +-- | Index structures for complex, heterogeneous indexing. Mostly designed for +-- indexing in DP grammars, where the indices work for linear and context-free +-- grammars on one or more tapes, for strings, sets, later on tree structures. + +class Index i where + -- | Data structure encoding the upper limit for each array. + data LimitType i :: * + -- | Given a maximal size, and a current index, calculate + -- the linear index. + linearIndex :: LimitType i -> i -> Int + -- | Given a maximal size and a valid @Int@, return the index. + fromLinearIndex :: LimitType i -> Int -> i + -- | Given the 'LimitType', return the number of cells required for storage. + size :: LimitType i -> Int + -- | Check if an index is within the bounds. + inBounds :: LimitType i -> i -> Bool + -- | A lower bound of @zero@ + zeroBound :: i + -- | A lower bound of @zero@ but for a @LimitType i@. + zeroBound' :: LimitType i + -- | The list of cell sizes for each dimension. its product yields the total + -- size. + totalSize :: LimitType i -> [Integer] + -- | Pretty-print all upper bounds + showBound :: LimitType i -> [String] + -- | Pretty-print all indices + showIndex :: i -> [String] + +-- | Given the maximal number of cells (@Word@, because this is the pointer +-- limit for the machine), and the list of sizes, will check if this is still +-- legal. Consider dividing the @Word@ by the actual memory requirements for +-- each cell, to get better exception handling for too large arrays. +-- +-- One list should be given for each array. + +sizeIsValid :: Monad m => Word -> [[Integer]] -> ExceptT SizeError m CellSize +sizeIsValid maxCells cells = do + let ps = map product cells + s = sum ps + when (fromIntegral maxCells <= s) $ + throwError . SizeError + $ printf "PrimitiveArrays would be larger than maximal cell size. The given limit is %d, but the requested size is %d, with size %s for each array. (Debug hint: %s)" + maxCells s (show ps) (show s) + return . CellSize $ fromIntegral s +{-# Inlinable sizeIsValid #-} + +-- | In case @totalSize@ or variants thereof produce a size that is too big to +-- handle. + +newtype SizeError = SizeError String + deriving (Eq,Ord,Show) + +-- | The total number of cells that are allocated. + +newtype CellSize = CellSize Word + deriving stock (Eq,Ord,Show) + deriving newtype (Num,Bounded,Integral,Real,Enum) + + + +-- | Generate a stream of indices in correct order for dynamic programming. +-- Since the stream generators require @concatMap@ / @flatten@ we have to +-- write more specialized code for @(z:.IX)@ stuff. + +class (Index i) => IndexStream i where + -- | Generate an index stream using 'LimitType's. This prevents having to + -- figure out how the actual limits for complicated index types (like @Set@) + -- would look like, since for @Set@, for example, the @LimitType Set == Int@ + -- provides just the number of bits. + -- + -- This generates an index stream suitable for @forward@ structure filling. + -- The first index is the smallest (or the first indices considered are all + -- equally small in partially ordered sets). Larger indices follow up until + -- the largest one. + streamUp :: Monad m => LimitType i -> LimitType i -> Stream m i + -- | If 'streamUp' generates indices from smallest to largest, then + -- 'streamDown' generates indices from largest to smallest. Outside grammars + -- make implicit use of this. Asking for an axiom in backtracking requests + -- the first element from this stream. + streamDown :: Monad m => LimitType i -> LimitType i -> Stream m i + + + +instance Index Z where + data LimitType Z = ZZ + linearIndex _ _ = 0 + {-# INLINE linearIndex #-} + fromLinearIndex _ _ = Z + {-# Inline fromLinearIndex #-} + size _ = 1 + {-# INLINE size #-} + inBounds _ _ = True + {-# INLINE inBounds #-} + zeroBound = Z + {-# Inline zeroBound #-} + zeroBound' = ZZ + {-# Inline zeroBound' #-} + totalSize ZZ = [1] + {-# Inline [1] totalSize #-} + showBound ZZ = [show ZZ] + showIndex Z = [show Z] + +instance IndexStream Z where + streamUp ZZ ZZ = SM.singleton Z + {-# Inline streamUp #-} + streamDown ZZ ZZ = SM.singleton Z + {-# Inline streamDown #-} + +deriving instance Eq (LimitType Z) +deriving instance Generic (LimitType Z) +deriving instance Read (LimitType Z) +deriving instance Show (LimitType Z) +deriving instance Data (LimitType Z) +deriving instance Typeable (LimitType Z) +deriving instance Bounded (LimitType Z) + +instance (Index zs, Index z) => Index (zs:.z) where + data LimitType (zs:.z) = !(LimitType zs) :.. !(LimitType z) + linearIndex (hs:..h) (zs:.z) = linearIndex hs zs * size h + linearIndex h z + {-# INLINE linearIndex #-} + fromLinearIndex (hs:..h) k = let (l , r) = quotRemInt k (size h) + in fromLinearIndex hs l :. fromLinearIndex h r + {-# Inline fromLinearIndex #-} + size (hs:..h) = size hs * size h + {-# INLINE size #-} + inBounds (hs:..h) (zs:.z) = inBounds hs zs && inBounds h z + {-# INLINE inBounds #-} + zeroBound = zeroBound :. zeroBound + {-# Inline zeroBound #-} + zeroBound' = zeroBound' :.. zeroBound' + {-# Inline zeroBound' #-} + totalSize (hs:..h) = + let tshs = totalSize hs + tsh = totalSize h + in tshs ++ tsh + {-# Inline totalSize #-} + showBound (zs:..z) = showBound zs ++ showBound z + showIndex (zs:.z) = showIndex zs ++ showIndex z + +deriving instance (Eq (LimitType zs) , Eq (LimitType z) ) => Eq (LimitType (zs:.z)) +deriving instance (Generic (LimitType zs), Generic (LimitType z)) => Generic (LimitType (zs:.z)) +deriving instance (Read (LimitType zs) , Read (LimitType z) ) => Read (LimitType (zs:.z)) +deriving instance (Show (LimitType zs) , Show (LimitType z) ) => Show (LimitType (zs:.z)) +deriving instance + ( Data zs, Data (LimitType zs), Typeable zs + , Data z , Data (LimitType z) , Typeable z + ) => Data (LimitType (zs:.z)) +deriving instance (Bounded (LimitType zs), Bounded (LimitType z)) => Bounded (LimitType (zs:.z)) + +--instance (Index zs, Index z) => Index (zs:>z) where +-- type LimitType (zs:>z) = LimitType zs:>LimitType z +-- linearIndex (hs:>h) (zs:>z) = linearIndex hs zs * (size (Proxy :: Proxy z) h) + linearIndex h z +-- {-# INLINE linearIndex #-} +-- size Proxy (ss:>s) = size (Proxy :: Proxy zs) ss * (size (Proxy :: Proxy z) s) +-- {-# INLINE size #-} +-- inBounds (hs:>h) (zs:>z) = inBounds hs zs && inBounds h z +-- {-# INLINE inBounds #-} + + + +-- * Somewhat experimental lens support. +-- +-- The problem here is that tuples are n-ary, while inductive tuples are +-- binary, recursive. + +instance Field1 (Z:.a) (Z:.a') a a' where + {-# Inline _1 #-} + _1 = lens (\(Z:.a) -> a) (\(Z:._) a -> (Z:.a)) + +instance Field1 (Z:.a:.b) (Z:.a':.b) a a' where + {-# Inline _1 #-} + _1 = lens (\(Z:.a:.b) -> a) (\(Z:._:.b) a -> (Z:.a:.b)) + +instance Field1 (Z:.a:.b:.c) (Z:.a':.b:.c) a a' where + {-# Inline _1 #-} + _1 = lens (\(Z:.a:.b:.c) -> a) (\(Z:._:.b:.c) a -> (Z:.a:.b:.c)) + + +instance Field2 (Z:.a:.b) (Z:.a:.b') b b' where + {-# Inline _2 #-} + _2 = lens (\(Z:.a:.b) -> b) (\(Z:.a:._) b -> (Z:.a:.b)) + +instance Field2 (Z:.a:.b:.c) (Z:.a:.b':.c) b b' where + {-# Inline _2 #-} + _2 = lens (\(Z:.a:.b:.c) -> b) (\(Z:.a:._:.c) b -> (Z:.a:.b:.c)) + + +instance Field3 (Z:.a:.b:.c) (Z:.a:.b:.c') c c' where + {-# Inline _3 #-} + _3 = lens (\(Z:.a:.b:.c) -> c) (\(Z:.a:.b:._) c -> (Z:.a:.b:.c)) + + + +-- * Operations for sparsity. + +-- | @manhattan@ turns an index @sh@ into a starting point within 'sparseIndices' of the 'Sparse' +-- data structure. This should reduce the time required to search @sparseIndices@, because +-- @manhattanStart[manhattan sh]@ yields a left bound, while @manhattanStart[manhattan sh +1]@ will +-- yield an excluded right bound. +-- +-- Uses the @Manhattan@ distance. +-- +-- TODO This should probably be moved into the @Index@ module. + +class SparseBucket sh where + -- | The manhattan distance for an index. + manhattan :: LimitType sh -> sh -> Int + -- | The maximal possible manhattan distance. + manhattanMax :: LimitType sh -> Int + +instance SparseBucket Z where + {-# Inline manhattan #-} + manhattan ZZ Z = 0 + {-# Inline manhattanMax #-} + manhattanMax ZZ = 1 + +-- | Manhattan distances add up. + +instance (SparseBucket i, SparseBucket is) => SparseBucket (is:.i) where + {-# Inline manhattan #-} + manhattan (zz:..z) (is:.i) = manhattan zz is + manhattan z i + {-# Inline manhattanMax #-} + manhattanMax (zz:..z) = manhattanMax zz + manhattanMax z + diff --git a/Data/PrimitiveArray/Index/IOC.hs b/Data/PrimitiveArray/Index/IOC.hs new file mode 100644 --- /dev/null +++ b/Data/PrimitiveArray/Index/IOC.hs @@ -0,0 +1,17 @@ + +module Data.PrimitiveArray.Index.IOC where + + + +-- | Phantom type for @Inside@ indices. + +data I + +-- | Phantom type for @Outside@ indices. + +data O + +-- | Phantom type for @Complement@ indices. + +data C + diff --git a/Data/PrimitiveArray/Index/Int.hs b/Data/PrimitiveArray/Index/Int.hs new file mode 100644 --- /dev/null +++ b/Data/PrimitiveArray/Index/Int.hs @@ -0,0 +1,54 @@ + +module Data.PrimitiveArray.Index.Int where + +import qualified Data.Vector.Fusion.Stream.Monadic as SM + +import Data.PrimitiveArray.Index.Class + + + +instance Index Int where + newtype LimitType Int = LtInt Int + linearIndex _ k = k + {-# Inline linearIndex #-} + size (LtInt h) = h+1 + {-# Inline size #-} + inBounds (LtInt h) k = 0 <= k && k <= h + {-# Inline inBounds #-} + zeroBound = 0 + {-# Inline [0] zeroBound #-} + zeroBound' = LtInt 0 + {-# Inline [0] zeroBound' #-} + totalSize (LtInt h) = [fromIntegral $ h+1] + {-# Inline [0] totalSize #-} + fromLinearIndex _ = id + {-# Inline [0] fromLinearIndex #-} + showBound (LtInt b) = ["LtInt " ++ show b] + showIndex i = ["Int " ++ show i] + +deriving instance Show (LimitType Int) + +instance IndexStream z => IndexStream (z:.Int) where + streamUp (ls:.. LtInt l) (hs:.. LtInt h) = SM.flatten mk step $ streamUp ls hs + where mk z = return (z,l) + step (z,k) + | k > h = return $ SM.Done + | otherwise = return $ SM.Yield (z:.k) (z,k+1) + {-# Inline [0] mk #-} + {-# Inline [0] step #-} + {-# Inline streamUp #-} + streamDown (ls:..LtInt l) (hs:..LtInt h) = SM.flatten mk step $ streamDown ls hs + where mk z = return (z,h) + step (z,k) + | k < l = return $ SM.Done + | otherwise = return $ SM.Yield (z:.k) (z,k-1) + {-# Inline [0] mk #-} + {-# Inline [0] step #-} + {-# Inline streamDown #-} + +instance IndexStream Int where + streamUp l h = SM.map (\(Z:.k) -> k) $ streamUp (ZZ:..l) (ZZ:..h) + {-# Inline streamUp #-} + streamDown l h = SM.map (\(Z:.k) -> k) $ streamDown (ZZ:..l) (ZZ:..h) + {-# Inline streamDown #-} + diff --git a/Data/PrimitiveArray/Index/PhantomInt.hs b/Data/PrimitiveArray/Index/PhantomInt.hs new file mode 100644 --- /dev/null +++ b/Data/PrimitiveArray/Index/PhantomInt.hs @@ -0,0 +1,115 @@ + +-- | A linear 0-based int-index with a phantom type. + +module Data.PrimitiveArray.Index.PhantomInt where + +import Control.DeepSeq (NFData(..)) +import Data.Aeson (FromJSON,FromJSONKey,ToJSON,ToJSONKey) +import Data.Binary (Binary) +import Data.Data +import Data.Hashable (Hashable) +import Data.Ix(Ix) +import Data.Serialize (Serialize) +import Data.Typeable +import Data.Vector.Fusion.Stream.Monadic (map,Step(..),flatten) +import Data.Vector.Unboxed.Deriving +import GHC.Generics (Generic) +import Prelude hiding (map) + +import Data.PrimitiveArray.Index.Class +import Data.PrimitiveArray.Index.IOC + + + +-- | A 'PInt' behaves exactly like an @Int@, but has an attached phantom +-- type @p@. In particular, the @Index@ and @IndexStream@ instances are the +-- same as for raw @Int@s. + +newtype PInt (ioc ∷ k) (p ∷ k) = PInt { getPInt :: Int } + deriving stock (Read,Show,Eq,Ord,Generic,Data,Typeable,Ix) + deriving newtype (Real,Num,Enum,Integral) + +pIntI :: Int -> PInt I p +pIntI = PInt +{-# Inline pIntI #-} + +pIntO :: Int -> PInt O p +pIntO = PInt +{-# Inline pIntO #-} + +pIntC :: Int -> PInt C p +pIntC = PInt +{-# Inline pIntC #-} + +derivingUnbox "PInt" + [t| forall t p . PInt t p -> Int |] [| getPInt |] [| PInt |] + +instance Binary (PInt t p) +instance Serialize (PInt t p) +instance FromJSON (PInt t p) +instance FromJSONKey (PInt t p) +instance ToJSON (PInt t p) +instance ToJSONKey (PInt t p) +instance Hashable (PInt t p) +instance NFData (PInt t p) + +instance Index (PInt t p) where + newtype LimitType (PInt t p) = LtPInt Int + linearIndex _ (PInt k) = k + {-# Inline linearIndex #-} + size (LtPInt h) = h+1 + {-# Inline size #-} + inBounds (LtPInt h) (PInt k) = 0 <= k && k <= h + {-# Inline inBounds #-} + fromLinearIndex = error "implement me" + zeroBound = error "implement me" + zeroBound' = error "implement me" + totalSize = error "implement me" + showBound = error "implement me" + showIndex = error "implement me" + +deriving instance Show (LimitType (PInt t p)) +deriving instance Read (LimitType (PInt t p)) +deriving instance Eq (LimitType (PInt t p)) +deriving instance Generic (LimitType (PInt t p)) + +instance IndexStream z => IndexStream (z:.PInt I p) where + streamUp (ls:..LtPInt l) (hs:..LtPInt h) = flatten (streamUpMk l h) (streamUpStep l h) $ streamUp ls hs + streamDown (ls:..LtPInt l) (hs:..LtPInt h) = flatten (streamDownMk l h) (streamDownStep l h) $ streamDown ls hs + {-# Inline streamUp #-} + {-# Inline streamDown #-} + +instance IndexStream z => IndexStream (z:.PInt O p) where + streamUp (ls:..LtPInt l) (hs:..LtPInt h) = flatten (streamDownMk l h) (streamDownStep l h) $ streamUp ls hs + streamDown (ls:..LtPInt l) (hs:..LtPInt h) = flatten (streamUpMk l h) (streamUpStep l h) $ streamDown ls hs + {-# Inline streamUp #-} + {-# Inline streamDown #-} + +instance IndexStream z => IndexStream (z:.PInt C p) where + streamUp (ls:..LtPInt l) (hs:..LtPInt h) = flatten (streamUpMk l h) (streamUpStep l h) $ streamUp ls hs + streamDown (ls:..LtPInt l) (hs:..LtPInt h) = flatten (streamDownMk l h) (streamDownStep l h) $ streamDown ls hs + {-# Inline streamUp #-} + {-# Inline streamDown #-} + +instance IndexStream (Z:.PInt ioc p) => IndexStream (PInt ioc p) where + streamUp l h = map (\(Z:.i) -> i) $ streamUp (ZZ:..l) (ZZ:..h) + {-# INLINE streamUp #-} + streamDown l h = map (\(Z:.i) -> i) $ streamDown (ZZ:..l) (ZZ:..h) + {-# INLINE streamDown #-} + +streamUpMk l h z = return (z,l) +{-# Inline [0] streamUpMk #-} + +streamUpStep l h (z,k) + | k > h = return $ Done + | otherwise = return $ Yield (z:.PInt k) (z,k+1) +{-# Inline [0] streamUpStep #-} + +streamDownMk l h z = return (z,h) +{-# Inline [0] streamDownMk #-} + +streamDownStep l h (z,k) + | k < l = return $ Done + | otherwise = return $ Yield (z:.PInt k) (z,k-1) +{-# Inline [0] streamDownStep #-} + diff --git a/Data/PrimitiveArray/Index/Point.hs b/Data/PrimitiveArray/Index/Point.hs new file mode 100644 --- /dev/null +++ b/Data/PrimitiveArray/Index/Point.hs @@ -0,0 +1,258 @@ + +{-# Language MagicHash #-} + +-- | @Point@ index structures are used for left- and right-linear grammars. +-- Such grammars have at most one syntactic symbol on each r.h.s. of a rule. +-- The syntactic symbol needs to be in an outermost position. + +module Data.PrimitiveArray.Index.Point where + +import Control.Applicative +import Control.DeepSeq (NFData(..)) +import Data.Aeson +import Data.Binary +import Data.Bits +import Data.Bits.Extras (Ranked) +import Data.Hashable (Hashable) +import Data.Serialize +import Data.Vector.Unboxed.Deriving +import Data.Vector.Unboxed (Unbox(..)) +import GHC.Exts +import GHC.Generics (Generic) +import qualified Data.Vector.Fusion.Stream.Monadic as SM +import qualified Data.Vector.Unboxed as VU +import Test.QuickCheck as TQ +import Test.SmallCheck.Series as TS + +import Data.PrimitiveArray.Index.Class +import Data.PrimitiveArray.Index.IOC + + + +-- | A point in a left-linear grammar. The syntactic symbol is in left-most +-- position. + +newtype PointL t = PointL {fromPointL :: Int} + deriving stock (Eq,Ord,Read,Show,Generic) + deriving newtype (Num) + +pointLI :: Int -> PointL I +pointLI = PointL +{-# Inline pointLI #-} + +pointLO :: Int -> PointL O +pointLO = PointL +{-# Inline pointLO #-} + +pointLC :: Int -> PointL C +pointLC = PointL +{-# Inline pointLC #-} + + + +derivingUnbox "PointL" + [t| forall t . PointL t -> Int |] + [| \ (PointL i) -> i |] + [| \ i -> PointL i |] + +instance Binary (PointL t) +instance Serialize (PointL t) +instance FromJSON (PointL t) +instance FromJSONKey (PointL t) +instance ToJSON (PointL t) +instance ToJSONKey (PointL t) +instance Hashable (PointL t) + +instance NFData (PointL t) where + rnf (PointL l) = rnf l + {-# Inline rnf #-} + +instance Index (PointL t) where + newtype LimitType (PointL t) = LtPointL Int + linearIndex _ (PointL z) = z + {-# INLINE linearIndex #-} + fromLinearIndex (LtPointL h) k = (PointL k) + {-# Inline fromLinearIndex #-} + size (LtPointL h) = h + 1 + {-# INLINE size #-} + inBounds (LtPointL h) (PointL x) = 0<=x && x<=h + {-# INLINE inBounds #-} + zeroBound = PointL 0 + {-# Inline [0] zeroBound #-} + zeroBound' = LtPointL 0 + {-# Inline [0] zeroBound' #-} + totalSize (LtPointL h) = [fromIntegral $ h + 1] + {-# Inline [0] totalSize #-} + showBound (LtPointL h) = ["LtPointL " ++ show h] + showIndex (PointL i) = ["PointL " ++ show i] + +deriving instance Eq (LimitType (PointL t)) +deriving instance Generic (LimitType (PointL t)) +deriving instance Read (LimitType (PointL t)) +deriving instance Show (LimitType (PointL t)) + +instance IndexStream z => IndexStream (z:.PointL I) where + streamUp (ls:..LtPointL lf) (hs:..LtPointL ht) = SM.flatten (streamUpMk lf) (streamUpStep PointL ht) $ streamUp ls hs + streamDown (ls:..LtPointL lf) (hs:..LtPointL ht) = SM.flatten (streamDownMk ht) (streamDownStep PointL lf) $ streamDown ls hs + {-# Inline [0] streamUp #-} + {-# Inline [0] streamDown #-} + +instance IndexStream z => IndexStream (z:.PointL O) where + streamUp (ls:..LtPointL lf) (hs:..LtPointL ht) = SM.flatten (streamDownMk ht) (streamDownStep PointL lf) $ streamUp ls hs + streamDown (ls:..LtPointL lf) (hs:..LtPointL ht) = SM.flatten (streamUpMk lf) (streamUpStep PointL ht) $ streamDown ls hs + {-# Inline [0] streamUp #-} + {-# Inline [0] streamDown #-} + +instance IndexStream z => IndexStream (z:.PointL C) where + streamUp (ls:..LtPointL lf) (hs:..LtPointL ht) = SM.flatten (streamUpMk lf) (streamUpStep PointL ht) $ streamUp ls hs + streamDown (ls:..LtPointL lf) (hs:..LtPointL ht) = SM.flatten (streamDownMk ht) (streamDownStep PointL lf) $ streamDown ls hs + {-# Inline [0] streamUp #-} + {-# Inline [0] streamDown #-} + +data SP z = SP !z !Int# + +streamUpMk (I# lf) z = return $ SP z lf +{-# Inline [0] streamUpMk #-} + +streamUpStep wrapper (I# ht) (SP z k) + | 1# <- k ># ht = return $ SM.Done + | otherwise = return $ SM.Yield (z:.wrapper (I# k)) (SP z (k +# 1#)) +{-# Inline [0] streamUpStep #-} + +streamDownMk (I# ht) z = return $ SP z ht +{-# Inline [0] streamDownMk #-} + +streamDownStep wrapper (I# lf) (SP z k) + | 1# <- k <# lf = return $ SM.Done + | otherwise = return $ SM.Yield (z:.wrapper (I# k)) (SP z (k -# 1#)) +{-# Inline [0] streamDownStep #-} + +instance IndexStream (Z:.PointL t) => IndexStream (PointL t) where + streamUp l h = SM.map (\(Z:.i) -> i) $ streamUp (ZZ:..l) (ZZ:..h) + {-# INLINE streamUp #-} + streamDown l h = SM.map (\(Z:.i) -> i) $ streamDown (ZZ:..l) (ZZ:..h) + {-# INLINE streamDown #-} + + +instance Arbitrary (PointL t) where + arbitrary = do + b <- choose (0,100) + return $ PointL b + shrink (PointL j) + | 0 Serial m (PointL t) where + series = PointL . TS.getNonNegative <$> series + + + +-- * @PointR@ + +-- | A point in a right-linear grammars. + +newtype PointR t = PointR {fromPointR :: Int} + deriving stock (Eq,Ord,Read,Show,Generic) + deriving newtype (Num) + + + +derivingUnbox "PointR" + [t| forall t . PointR t -> Int |] + [| \ (PointR i) -> i |] + [| \ i -> PointR i |] + +instance Binary (PointR t) +instance Serialize (PointR t) +instance FromJSON (PointR t) +instance FromJSONKey (PointR t) +instance ToJSON (PointR t) +instance ToJSONKey (PointR t) +instance Hashable (PointR t) + +instance NFData (PointR t) where + rnf (PointR l) = rnf l + {-# Inline rnf #-} + +instance Index (PointR t) where + newtype LimitType (PointR t) = LtPointR Int + linearIndex _ (PointR z) = z + {-# INLINE linearIndex #-} + size (LtPointR h) = h + 1 + {-# INLINE size #-} + inBounds (LtPointR h) (PointR x) = 0<=x && x<=h + {-# INLINE inBounds #-} + zeroBound = PointR 0 + {-# Inline [0] zeroBound #-} + zeroBound' = LtPointR 0 + {-# Inline [0] zeroBound' #-} + totalSize (LtPointR h) = [fromIntegral $ h + 1] + {-# Inline [0] totalSize #-} + fromLinearIndex _ = PointR + {-# Inline [0] fromLinearIndex #-} + showBound (LtPointR b) = ["LtPointR " ++ show b] + showIndex (PointR p) = ["PointR " ++ show p] + +deriving instance Eq (LimitType (PointR t)) +deriving instance Generic (LimitType (PointR t)) +deriving instance Read (LimitType (PointR t)) +deriving instance Show (LimitType (PointR t)) + +instance IndexStream z => IndexStream (z:.PointR I) where + streamUp (ls:..LtPointR lf) (hs:..LtPointR ht) = SM.flatten (streamDownMk ht) (streamDownStep PointR lf) $ streamUp ls hs + streamDown (ls:..LtPointR lf) (hs:..LtPointR ht) = SM.flatten (streamUpMk lf) (streamUpStep PointR ht) $ streamDown ls hs + {-# Inline [0] streamUp #-} + {-# Inline [0] streamDown #-} + +instance IndexStream z => IndexStream (z:.PointR O) where + streamUp (ls:..LtPointR lf) (hs:..LtPointR ht) = SM.flatten (streamUpMk lf) (streamUpStep PointR ht) $ streamUp ls hs + streamDown (ls:..LtPointR lf) (hs:..LtPointR ht) = SM.flatten (streamDownMk ht) (streamDownStep PointR lf) $ streamDown ls hs + {-# Inline [0] streamUp #-} + {-# Inline [0] streamDown #-} + +--instance IndexStream z => IndexStream (z:.PointR C) where +-- streamUp (ls:..LtPointR lf) (hs:..LtPointR ht) = SM.flatten (streamUpMkR lf) (streamUpStepR ht) $ streamUp ls hs +-- streamDown (ls:..LtPointR lf) (hs:..LtPointR ht) = SM.flatten (streamDownMkR ht) (streamDownStepR lf) $ streamDown ls hs +-- {-# Inline [0] streamUp #-} +-- {-# Inline [0] streamDown #-} + +instance IndexStream (Z:.PointR t) => IndexStream (PointR t) where + streamUp l h = SM.map (\(Z:.i) -> i) $ streamUp (ZZ:..l) (ZZ:..h) + {-# INLINE streamUp #-} + streamDown l h = SM.map (\(Z:.i) -> i) $ streamDown (ZZ:..l) (ZZ:..h) + {-# INLINE streamDown #-} + +-- arbitrarily set maximum here to + +arbMaxPointR = 100 + +instance Arbitrary (PointR t) where + arbitrary = do + b <- choose (0,arbMaxPointR) + return $ PointR b + shrink (PointR j) + | j Serial m (PointR t) where +-- series = PointR . TS.getNonNegative <$> series + + + +instance SparseBucket (PointL I) where + {-# Inline manhattan #-} + manhattan (LtPointL u) (PointL i) = i + {-# Inline manhattanMax #-} + manhattanMax (LtPointL u) = u + + +-- | +-- +-- TODO Is this instance correct? Outside indices shrink? + +instance SparseBucket (PointL O) where + {-# Inline manhattan #-} + manhattan (LtPointL u) (PointL i) = u-i + {-# Inline manhattanMax #-} + manhattanMax (LtPointL u) = u + diff --git a/Data/PrimitiveArray/Index/Subword.hs b/Data/PrimitiveArray/Index/Subword.hs new file mode 100644 --- /dev/null +++ b/Data/PrimitiveArray/Index/Subword.hs @@ -0,0 +1,181 @@ + +-- | Index structure for context-free grammars on strings. A @Subword@ captures +-- a pair @(i,j)@ with @i<=j@. + +module Data.PrimitiveArray.Index.Subword where + +import Control.Applicative ((<$>)) +import Control.DeepSeq (NFData(..)) +import Control.Monad (filterM, guard) +import Data.Aeson (FromJSON,FromJSONKey,ToJSON,ToJSONKey) +import Data.Binary (Binary) +import Data.Hashable (Hashable) +import Data.Serialize (Serialize) +import Data.Vector.Fusion.Stream.Monadic (Step(..), map,flatten) +import Data.Vector.Unboxed.Deriving +import GHC.Generics (Generic) +import Prelude hiding (map) +import Test.QuickCheck (Arbitrary(..), choose) +import Test.SmallCheck.Series as TS + +import Math.TriangularNumbers + +import Data.PrimitiveArray.Index.Class +import Data.PrimitiveArray.Index.IOC + + + +-- | A subword wraps a pair of @Int@ indices @i,j@ with @i<=j@. +-- +-- Subwords always yield the upper-triangular part of a rect-angular array. +-- This gives the quite curious effect that @(0,N)@ points to the +-- ``largest'' index, while @(0,0) ... (1,1) ... (k,k) ... (N,N)@ point to +-- the smallest. We do, however, use (0,0) as the smallest as (0,k) gives +-- successively smaller upper triangular parts. + +newtype Subword t = Subword {fromSubword :: (Int:.Int)} + deriving (Eq,Ord,Show,Generic,Read) + +fromSubwordFst :: Subword t -> Int +fromSubwordFst (Subword (i:._)) = i +{-# Inline fromSubwordFst #-} + +fromSubwordSnd :: Subword t -> Int +fromSubwordSnd (Subword (_:.j)) = j +{-# Inline fromSubwordSnd #-} + +derivingUnbox "Subword" + [t| forall t . Subword t -> (Int,Int) |] + [| \ (Subword (i:.j)) -> (i,j) |] + [| \ (i,j) -> Subword (i:.j) |] + +instance Binary (Subword t) +instance Serialize (Subword t) +instance FromJSON (Subword t) +instance FromJSONKey (Subword t) +instance ToJSON (Subword t) +instance ToJSONKey (Subword t) +instance Hashable (Subword t) + +instance NFData (Subword t) where + rnf (Subword (i:.j)) = i `seq` rnf j + {-# Inline rnf #-} + +-- | Create a @Subword t@ where @t@ is inferred. + +subword :: Int -> Int -> Subword t +subword i j = Subword (i:.j) +{-# INLINE subword #-} + +subwordI :: Int -> Int -> Subword I +subwordI i j = Subword (i:.j) +{-# INLINE subwordI #-} + +subwordO :: Int -> Int -> Subword O +subwordO i j = Subword (i:.j) +{-# INLINE subwordO #-} + +subwordC :: Int -> Int -> Subword C +subwordC i j = Subword (i:.j) +{-# INLINE subwordC #-} + + + +instance Index (Subword t) where + newtype LimitType (Subword t) = LtSubword Int + linearIndex (LtSubword n) (Subword (i:.j)) = toLinear n (i,j) + {-# Inline linearIndex #-} + size (LtSubword n) = linearizeUppertri (0,n) + {-# Inline size #-} + inBounds (LtSubword h) (Subword (i:.j)) = 0<=i && i<=j && j<=h + {-# Inline inBounds #-} + zeroBound = subword 0 0 + {-# Inline zeroBound #-} + zeroBound' = LtSubword 0 + {-# Inline zeroBound' #-} + totalSize (LtSubword n) = [fromIntegral (n+1) ^ 2 `div` 2] + {-# Inline totalSize #-} + fromLinearIndex = error "implement me" + showBound = error "implement me" + showIndex = error "implement me" + +deriving instance Eq (LimitType (Subword t)) +deriving instance Generic (LimitType (Subword t)) +deriving instance Read (LimitType (Subword t)) +deriving instance Show (LimitType (Subword t)) + +-- | @Subword I@ (inside) + +instance IndexStream z => IndexStream (z:.Subword I) where + streamUp (ls:..LtSubword l) (hs:..LtSubword h) = flatten (streamUpMk h) (streamUpStep l h) $ streamUp ls hs + streamDown (ls:..LtSubword l) (hs:..LtSubword h) = flatten (streamDownMk l h) (streamDownStep h) $ streamDown ls hs + {-# Inline streamUp #-} + {-# Inline streamDown #-} + +-- | @Subword O@ (outside). +-- +-- Note: @streamUp@ really needs to use @streamDownMk@ / @streamDownStep@ +-- for the right order of indices! + +instance IndexStream z => IndexStream (z:.Subword O) where + streamUp (ls:..LtSubword l) (hs:..LtSubword h) = flatten (streamDownMk l h) (streamDownStep h) $ streamUp ls hs + streamDown (ls:..LtSubword l) (hs:..LtSubword h) = flatten (streamUpMk h) (streamUpStep l h) $ streamDown ls hs + {-# Inline streamUp #-} + {-# Inline streamDown #-} + +-- | @Subword C@ (complement) + +instance IndexStream z => IndexStream (z:.Subword C) where + streamUp (ls:..LtSubword l) (hs:..LtSubword h) = flatten (streamUpMk h) (streamUpStep l h) $ streamUp ls hs + streamDown (ls:..LtSubword l) (hs:..LtSubword h) = flatten (streamDownMk l h) (streamDownStep h) $ streamDown ls hs + {-# Inline streamUp #-} + {-# Inline streamDown #-} + +-- | generic @mk@ for @streamUp@ / @streamDown@ + +streamUpMk h z = return (z,h,h) +{-# Inline [0] streamUpMk #-} + +streamUpStep l h (z,i,j) + | i < l = return $ Done + | j > h = return $ Skip (z,i-1,i-1) + | otherwise = return $ Yield (z:.subword i j) (z,i,j+1) +{-# Inline [0] streamUpStep #-} + +streamDownMk l h z = return (z,l,h) +{-# Inline [0] streamDownMk #-} + +streamDownStep h (z,i,j) + | i > h = return $ Done + | j < i = return $ Skip (z,i+1,h) + | otherwise = return $ Yield (z:.subword i j) (z,i,j-1) +{-# Inline [0] streamDownStep #-} + +instance (IndexStream (Z:.Subword t)) => IndexStream (Subword t) where + streamUp l h = map (\(Z:.i) -> i) $ streamUp (ZZ:..l) (ZZ:..h) + {-# INLINE streamUp #-} + streamDown l h = map (\(Z:.i) -> i) $ streamDown (ZZ:..l) (ZZ:..h) + {-# INLINE streamDown #-} + +instance Arbitrary (Subword t) where + arbitrary = do + a <- choose (0,20) + b <- choose (0,20) + return $ Subword (min a b :. max a b) + shrink (Subword (i:.j)) + | i Serial m (Subword t) where + series = do + i <- TS.getNonNegative <$> series + j <- TS.getNonNegative <$> series + guard $ i<=j + return $ subword i j + {- + let nns :: Series m Int = TS.getNonNegative <$> series + ps <- nns >< nns + let qs = [ subword i j | (i,j) <- ps, i<=j ] + return qs + -} + diff --git a/Data/PrimitiveArray/Index/Unit.hs b/Data/PrimitiveArray/Index/Unit.hs new file mode 100644 --- /dev/null +++ b/Data/PrimitiveArray/Index/Unit.hs @@ -0,0 +1,82 @@ + +-- | Unit indices admit a single element to be memoized. We can't use @()@ +-- because we want to attach phantom types. + +module Data.PrimitiveArray.Index.Unit where + +import Control.Applicative (pure) +import Control.DeepSeq (NFData(..)) +import Data.Aeson (FromJSON,FromJSONKey,ToJSON,ToJSONKey) +import Data.Binary (Binary) +import Data.Hashable (Hashable) +import Data.Serialize (Serialize) +import Data.Vector.Fusion.Stream.Monadic (Step(..), map) +import Data.Vector.Unboxed.Deriving +import GHC.Generics (Generic) +import Prelude hiding (map) +import Test.QuickCheck (Arbitrary(..), choose) + +import Data.PrimitiveArray.Index.Class + + + +data Unit t = Unit + deriving (Eq,Ord,Show,Generic,Read) + +derivingUnbox "Unit" + [t| forall t . Unit t -> () |] + [| \ Unit -> () |] + [| \ () -> Unit |] + +instance Binary (Unit t) +instance Serialize (Unit t) +instance FromJSON (Unit t) +instance FromJSONKey (Unit t) +instance ToJSON (Unit t) +instance ToJSONKey (Unit t) +instance Hashable (Unit t) + +instance NFData (Unit t) where + rnf Unit = () + {-# Inline rnf #-} + +instance Index (Unit t) where + data LimitType (Unit t) = LtUnit + linearIndex _ _ = 0 + {-# Inline linearIndex #-} + size _ = 1 + {-# Inline size #-} + inBounds _ _ = True + {-# Inline inBounds #-} + zeroBound = Unit + {-# Inline zeroBound #-} + zeroBound' = LtUnit + {-# Inline zeroBound' #-} + totalSize LtUnit = return 1 + {-# Inline [0] totalSize #-} + fromLinearIndex _ _ = Unit + {-# Inline fromLinearIndex #-} + showBound _ = ["LtUnit"] + showIndex _ = ["Unit"] + +deriving instance Eq (LimitType (Unit t)) +deriving instance Generic (LimitType (Unit t)) +deriving instance Read (LimitType (Unit t)) +deriving instance Show (LimitType (Unit t)) + +instance IndexStream z => IndexStream (z:.Unit t) where + streamUp (ls:..LtUnit) (hs:..LtUnit) = map (\z -> z:.Unit) $ streamUp ls hs + {-# Inline streamUp #-} + streamDown (ls:..LtUnit) (hs:..LtUnit) = map (\z -> z:.Unit) $ streamDown ls hs + {-# Inline streamDown #-} + +instance (IndexStream (Z:.Unit t)) => IndexStream (Unit t) where + streamUp l h = map (\(Z:.i) -> i) $ streamUp (ZZ:..l) (ZZ:..h) + {-# INLINE streamUp #-} + streamDown l h = map (\(Z:.i) -> i) $ streamDown (ZZ:..l) (ZZ:..h) + {-# INLINE streamDown #-} + +instance Arbitrary (Unit t) where + arbitrary = pure Unit + shrink Unit = [] + diff --git a/Data/PrimitiveArray/Sparse.hs b/Data/PrimitiveArray/Sparse.hs new file mode 100644 --- /dev/null +++ b/Data/PrimitiveArray/Sparse.hs @@ -0,0 +1,7 @@ + +module Data.PrimitiveArray.Sparse + ( module Data.PrimitiveArray.Sparse.IntBinSearch + ) where + +import Data.PrimitiveArray.Sparse.IntBinSearch + diff --git a/Data/PrimitiveArray/Sparse/BinSearch.hs b/Data/PrimitiveArray/Sparse/BinSearch.hs new file mode 100644 --- /dev/null +++ b/Data/PrimitiveArray/Sparse/BinSearch.hs @@ -0,0 +1,233 @@ + +-- | This solution to holding a sparse set of elements for dynamic programming. The underlying +-- representation requires @O (log n)@ access time for each read or write, where @n@ is the number +-- of elements to be stored. It uses an experimental "bucketing" system to provide better lower and +-- upper bounds than otherwise possible. +-- +-- TODO @ADPfusion / FillTyLvl@ handles actually filling the tables. In case all @BigOrder@ tables +-- are dense and of the same dimensional extent, we are fine. However if at least one table is +-- dense, while others are sparse, we will have write to nothing, which should not crash. In case of +-- all-sparse tables for a BigOrder, we have to calculate the union of all indices. This all is +-- currently not happening... +-- +-- TODO require @readMaybe@ and @indexMaybe@ to return @Nothing@ on missing elements. This requires +-- an extension of the @Class@ structure for tables. + +module Data.PrimitiveArray.Sparse.BinSearch where + +import Control.Monad.Primitive (PrimState,PrimMonad) +import Control.Monad.ST (ST) +import Debug.Trace (traceShow) +import qualified Control.Monad.State.Strict as SS +import qualified Data.HashMap.Strict as HMS +import qualified Data.Vector.Algorithms.Intro as VAI +import qualified Data.Vector.Algorithms.Search as VAS +import qualified Data.Vector as V +import qualified Data.Vector.Fusion.Stream.Monadic as SM +import qualified Data.Vector.Generic as VG +import qualified Data.Vector.Generic.Mutable as VGM +import qualified Data.Vector.Storable as VS +import qualified Data.Vector.Unboxed as VU + +import Data.PrimitiveArray.Class +import Data.PrimitiveArray.Index.Class +import Data.PrimitiveArray.Index -- TODO only while @SparseBucket@ is here + + + +-- | This is a sparse matrix, where only a subset of indices have data associated. + +data Sparse w v sh e = Sparse + { sparseUpperBound :: !(LimitType sh) + -- ^ The upper bound for the DP matrix. Not the upper bound of indexes in use, but the theoretical + -- upper bound. + , sparseData :: !(v e) + -- ^ Vector with actually existing data. + , sparseIndices :: !(w sh) + -- ^ The index of each @sh@ is the same as of the corresponding @sparseData@ structure. Indices + -- should be ordered as required by the @streamUp@ function, to facilitate filling @Sparse@ by + -- going from left to right. + , manhattanStart :: !(VU.Vector Int) + -- ^ Provides left/right boundaries into @sparseIndices@ to speed up index search. Should be one + -- larger than the largest index to look up, to always provides a good excluded bound. + } + +type Unboxed w sh e = Sparse w VU.Vector sh e + +type Storable w sh e = Sparse w VS.Vector sh e + +type Boxed w sh e = Sparse w V.Vector sh e + + + +-- | Currently, our mutable variant of sparse matrices will keep indices and manhattan starts +-- immutable as well. + +data instance MutArr m (Sparse w v sh e) = MSparse + { msparseUpperBound :: !(LimitType sh) + , msparseData :: !(VG.Mutable v (PrimState m) e) + , msparseIndices :: !(w sh) -- (VG.Mutable w (PrimState m) sh) + , mmanhattanStart :: !(VU.Vector Int) -- (VU.MVector (PrimState m) Int) + } +-- deriving (Generic,Typeable) + + +type instance FillStruc (Sparse w v sh e) = (w sh) + + + +instance + ( Index sh, SparseBucket sh, Eq sh, Ord sh + , VG.Vector w sh , VG.Vector w (Int,sh), VG.Vector w (Int,(Int,sh)) + , VG.Vector v e +#if ADPFUSION_DEBUGOUTPUT + , Show sh, Show (LimitType sh), Show e +#endif + ) => PrimArrayOps (Sparse w v) sh e where + + -- ** pure operations + + {-# Inline upperBound #-} + upperBound Sparse{..} = sparseUpperBound + {-# Inline unsafeIndex #-} + unsafeIndex Sparse{..} idx = case manhattanIndex sparseUpperBound manhattanStart sparseIndices idx of + Nothing -> error "unsafeIndex of non-existing index" + Just v -> VG.unsafeIndex sparseData v + {-# Inline safeIndex #-} + safeIndex Sparse{..} = fmap (VG.unsafeIndex sparseData) . manhattanIndex sparseUpperBound manhattanStart sparseIndices + + -- ** monadic operations + + {-# Inline unsafeFreezeM #-} + unsafeFreezeM MSparse{..} = do + let sparseUpperBound = msparseUpperBound + sparseIndices = msparseIndices + manhattanStart = mmanhattanStart + sparseData <- VG.unsafeFreeze msparseData + return Sparse{..} + {-# Inline unsafeThawM #-} + unsafeThawM Sparse{..} = do + let msparseUpperBound = sparseUpperBound + msparseIndices = sparseIndices + mmanhattanStart = manhattanStart + msparseData <- VG.unsafeThaw sparseData + return MSparse{..} + {-# Inline upperBoundM #-} + upperBoundM MSparse{..} = msparseUpperBound + {-# Inline newM #-} + newM = error "not implemented, use newSM" + {-# Inline newWithM #-} + newWithM = error "not implemented, use newWithSM" + {-# Inline readM #-} + readM MSparse{..} idx = do + case manhattanIndex msparseUpperBound mmanhattanStart msparseIndices idx of + Nothing -> error "read of non-existing element" + Just v -> VGM.unsafeRead msparseData v + -- | Note that @writeM@ will fail loudly, because we can specialize in @FillTyLvl@ to use + -- non-failing writes. + {-# Inline writeM #-} + writeM MSparse{..} idx elm = do + case manhattanIndex msparseUpperBound mmanhattanStart msparseIndices idx of + Nothing -> error "read of non-existing element" + Just v -> VGM.unsafeWrite msparseData v elm + {-# Inline newSM #-} + newSM h fs' = do + fs <- VG.thaw (VG.map (\i -> (manhattan h i, i)) fs') >>= \v -> VAI.sort v >> VG.unsafeFreeze v + let msparseUpperBound = h + msparseIndices = VG.force $ VG.map snd fs + -- For any manhattan distance not found in the distances, we set to the length of the the + -- @msparseIndices@ vector. Perform reverse-scan to update all manhattan start distances. + go :: VU.MVector s Int -> ST s () + go mv = do + VG.forM_ (VG.reverse $ VG.indexed fs) $ \(i,(mh,_)) -> VGM.write mv mh i + mmanhattanStart = VG.modify go $ VG.replicate (manhattanMax h +1) (VG.length fs) + msparseData <- VGM.new $ VG.length msparseIndices + return $ MSparse {..} + {-# Inline newWithSM #-} + newWithSM h fs' e = do + mv <- newSM h fs' + VGM.set (msparseData mv) e + return mv + {-# Inline safeWriteM #-} + safeWriteM MSparse{..} sh e = case manhattanIndex msparseUpperBound mmanhattanStart msparseIndices sh of + Nothing -> return () + Just v -> VGM.unsafeWrite msparseData v e + {-# Inline safeReadM #-} + safeReadM MSparse{..} sh = case manhattanIndex msparseUpperBound mmanhattanStart msparseIndices sh of + Nothing -> return Nothing + Just v -> Just <$> VGM.unsafeRead msparseData v + -- ** implement me + transformShape = error "implement me" + fromListM = error "implement me" + + + + + +-- * Helper functions. + +-- | Find the index with manhattan helper +-- +-- TODO consider using binary search instead of a linear scan here! +-- e.g.: @k = VAS.binarySearchByBounds (==)@ +-- +-- NOTE running times with 100x100 DP problem "NeedlemanWunsch" +-- full findIndex of sixs: 0,050,000 cells/sec +-- using manhattan buckets, findIndex: 5,000,000 cells/sec +-- using binarySearch on slices: 11,000,000 cells/sec +-- +-- On a 1000x1000 DP NeedlemanWunsch problem, binary search on slices is at 6,500,000 cells/sec. + +manhattanIndex + :: (SparseBucket sh, VG.Vector v sh, Eq sh, Ord sh) + => LimitType sh -> Vector Int -> v sh -> sh -> Maybe Int +{-# Inline manhattanIndex #-} +manhattanIndex ub mstart sixs idx = fmap (+l) . binarySearch idx $ VG.unsafeSlice l (h-l+1) sixs + where + b = manhattan ub idx + -- lower and upper bucket bounds + l = mstart `VU.unsafeIndex` b + h = mstart `VU.unsafeIndex` (b+1) + +binarySearch :: (Eq sh, Ord sh, VG.Vector v sh) => sh -> v sh -> Maybe Int +{-# Inline binarySearch #-} +binarySearch e v = go 0 (VG.length v -1) + where + go :: Int -> Int -> Maybe Int + go !l !r = + let !m = (r+l) `div` 2 + !x = VG.unsafeIndex v m + in if r go l (m-1) + EQ -> Just m + GT -> go (m+1) r + + +-- | Given two index vectors of the same shape, will return the correctly ordered vector of +-- the union of indices. +-- +-- TODO This requires that @Ord (Shape O)@ uses the @Down@ instance of Ord! We need to fix this in +-- the @Index@ modules. +-- +-- TODO Rewrite to allow fusion without intermediate vectors using uncons. This will make it +-- possible to chain applications. @stream@ should be fine for this. + +mergeIndexVectors :: (Eq sh, Ord sh, VG.Vector w sh) => w sh -> w sh -> w sh +{-# Inlinable mergeIndexVectors #-} +mergeIndexVectors xs ys = VG.create $ do + let lxs = VG.length xs + lys = VG.length ys + mv <- VGM.new $ lxs + lys + let go !n !i !j + | i>=lxs && j>=lys = return n + | j>=lys = VG.unsafeIndexM xs i >>= VGM.unsafeWrite mv n >> go (n+1) (i+1) j + | i>=lxs = VG.unsafeIndexM ys j >>= VGM.unsafeWrite mv n >> go (n+1) i (j+1) + | otherwise = do + x <- VG.unsafeIndexM xs i + y <- VG.unsafeIndexM ys j + if | x==y -> VGM.unsafeWrite mv n x >> go (n+1) (i+1) (j+1) + | x< y -> VGM.unsafeWrite mv n x >> go (n+1) (i+1) j + | x> y -> VGM.unsafeWrite mv n y >> go (n+1) i (j+1) + n <- go 0 0 0 + return $ VGM.unsafeTake n mv + diff --git a/Data/PrimitiveArray/Sparse/IntBinSearch.hs b/Data/PrimitiveArray/Sparse/IntBinSearch.hs new file mode 100644 --- /dev/null +++ b/Data/PrimitiveArray/Sparse/IntBinSearch.hs @@ -0,0 +1,286 @@ + +{-# Language MagicHash #-} + +-- | This solution to holding a sparse set of elements for dynamic programming. The underlying +-- representation requires @O (log n)@ access time for each read or write, where @n@ is the number +-- of elements to be stored. It uses an experimental "bucketing" system to provide better lower and +-- upper bounds than otherwise possible. +-- +-- TODO @ADPfusion / FillTyLvl@ handles actually filling the tables. In case all @BigOrder@ tables +-- are dense and of the same dimensional extent, we are fine. However if at least one table is +-- dense, while others are sparse, we will have write to nothing, which should not crash. In case of +-- all-sparse tables for a BigOrder, we have to calculate the union of all indices. This all is +-- currently not happening... +-- +-- This version requires working @fromLinearIndex@ but is potentially faster. + +module Data.PrimitiveArray.Sparse.IntBinSearch where + +import Control.Monad.Primitive (PrimState,PrimMonad) +import Control.Monad.ST (ST) +import Data.Bits.Extras (msb) +import Debug.Trace (traceShow) +import qualified Control.Monad.State.Strict as SS +import qualified Data.HashMap.Strict as HMS +import qualified Data.Vector.Algorithms.Radix as Sort +import qualified Data.Vector.Algorithms.Search as VAS +import qualified Data.Vector as V +import qualified Data.Vector.Fusion.Stream.Monadic as SM +import qualified Data.Vector.Generic as VG +import qualified Data.Vector.Generic.Mutable as VGM +import qualified Data.Vector.Storable as VS +import qualified Data.Vector.Unboxed as VU +import GHC.Exts ( Int(..), Int#(..), (==#), (-#), (/=#), (*#), (+#), (<=#), remInt#, quotInt#, uncheckedIShiftRA#, (<#) ) + +import Data.PrimitiveArray.Class +import Data.PrimitiveArray.Index.Class +import Data.PrimitiveArray.Index -- TODO only while @SparseBucket@ is here + + + +-- | This is a sparse matrix, where only a subset of indices have data associated. + +data Sparse w v sh e = Sparse + { sparseUpperBound :: !(LimitType sh) + -- ^ The upper bound for the DP matrix. Not the upper bound of indexes in use, but the theoretical + -- upper bound. + , sparseData :: !(v e) + -- ^ Vector with actually existing data. + , sparseIndices :: !(VU.Vector Int) + -- ^ Linearly encoded sparse indices + , manhattanStart :: !(VU.Vector Int) + -- ^ Provides left/right boundaries into @sparseIndices@ to speed up index search. Should be one + -- larger than the largest index to look up, to always provides a good excluded bound. + } + +type Unboxed w sh e = Sparse w VU.Vector sh e + +type Storable w sh e = Sparse w VS.Vector sh e + +type Boxed w sh e = Sparse w V.Vector sh e + + + +-- | Currently, our mutable variant of sparse matrices will keep indices and manhattan starts +-- immutable as well. + +data instance MutArr m (Sparse w v sh e) = MSparse + { msparseUpperBound :: !(LimitType sh) + , msparseData :: !(VG.Mutable v (PrimState m) e) + , msparseIndices :: !(VU.Vector Int) + , mmanhattanStart :: !(VU.Vector Int) + } +-- deriving (Generic,Typeable) + + +type instance FillStruc (Sparse w v sh e) = (w sh) + + + +instance + ( Index sh, SparseBucket sh, Eq sh, Ord sh + , VG.Vector w sh , VG.Vector w (Int,sh), VG.Vector w (Int,(Int,sh)), VG.Vector w (Int,Int), VG.Vector w Int + , VG.Vector v e +#if ADPFUSION_DEBUGOUTPUT + , Show sh, Show (LimitType sh), Show e +#endif + ) => PrimArrayOps (Sparse w v) sh e where + + -- ** pure operations + + {-# Inline upperBound #-} + upperBound Sparse{..} = sparseUpperBound + {-# Inline unsafeIndex #-} + unsafeIndex Sparse{..} idx = case manhattanIndex sparseUpperBound manhattanStart sparseIndices idx of + Nothing -> error "unsafeIndex of non-existing index" + Just v -> VG.unsafeIndex sparseData v + {-# Inline safeIndex #-} + safeIndex Sparse{..} = fmap (VG.unsafeIndex sparseData) . manhattanIndex sparseUpperBound manhattanStart sparseIndices + + -- ** monadic operations + + {-# Inline unsafeFreezeM #-} + unsafeFreezeM MSparse{..} = do + let sparseUpperBound = msparseUpperBound + sparseIndices = msparseIndices + manhattanStart = mmanhattanStart + sparseData <- VG.unsafeFreeze msparseData + return Sparse{..} + {-# Inline unsafeThawM #-} + unsafeThawM Sparse{..} = do + let msparseUpperBound = sparseUpperBound + msparseIndices = sparseIndices + mmanhattanStart = manhattanStart + msparseData <- VG.unsafeThaw sparseData + return MSparse{..} + {-# Inline upperBoundM #-} + upperBoundM MSparse{..} = msparseUpperBound + {-# Inline newM #-} + newM = error "not implemented, use newSM" + {-# Inline newWithM #-} + newWithM = error "not implemented, use newWithSM" + {-# Inline readM #-} + readM MSparse{..} idx = do + case manhattanIndex msparseUpperBound mmanhattanStart msparseIndices idx of + Nothing -> error "read of non-existing element" + Just v -> VGM.unsafeRead msparseData v + -- | Note that @writeM@ will fail loudly, because we can specialize in @FillTyLvl@ to use + -- non-failing writes. + {-# Inline writeM #-} + writeM MSparse{..} idx elm = do + case manhattanIndex msparseUpperBound mmanhattanStart msparseIndices idx of + Nothing -> error "read of non-existing element" + Just v -> VGM.unsafeWrite msparseData v elm + {-# Inline newSM #-} + newSM h fs' = do + let msparseUpperBound = h + -- sort sparse indices by (manhattan, linearIndex) + {-# Inline srt #-} + srt x y = let ix = fromLinearIndex h x + iy = fromLinearIndex h y + in compare (manhattan h ix, x) (manhattan h iy, y) + {-# Inline radixsrt #-} + radixsrt i x = let ix = fromLinearIndex h x in Sort.radix i (manhattan h ix, x) + msparseIndices <- do + marr <- VG.thaw (VU.convert $ VG.map (linearIndex h) fs') + Sort.sortBy (Sort.passes (undefined :: (Int,Int))) (Sort.size (undefined :: Int)) radixsrt marr + VG.unsafeFreeze marr + let -- For any manhattan distance not found in the distances, we set to the length of the the + -- @msparseIndices@ vector. Perform reverse-scan to update all manhattan start distances. + go :: VU.MVector s Int -> ST s () + {-# Inline go #-} + go mv = do + VG.forM_ (VG.reverse $ VG.indexed msparseIndices) $ \(i,k) -> let lix = fromLinearIndex h k; mh = manhattan h lix in VGM.write mv mh i + let mmanhattanStart = VU.modify go $ VG.replicate (manhattanMax h +1) (VG.length msparseIndices) + msparseData <- VGM.new $ VG.length msparseIndices + return $ MSparse {..} + {-# Inline newWithSM #-} + newWithSM h fs' e = do + mv <- newSM h fs' + VGM.set (msparseData mv) e + return mv + {-# Inline safeWriteM #-} + safeWriteM MSparse{..} sh e = case manhattanIndex msparseUpperBound mmanhattanStart msparseIndices sh of + Nothing -> return () + Just v -> VGM.unsafeWrite msparseData v e + {-# Inline safeReadM #-} + safeReadM MSparse{..} sh = case manhattanIndex msparseUpperBound mmanhattanStart msparseIndices sh of + Nothing -> return Nothing + Just v -> Just <$> VGM.unsafeRead msparseData v + -- ** implement me + transformShape = error "implement me" + fromListM = error "implement me" + + + +instance (Index sh, VG.Vector v e, VG.Vector v e') ⇒ PrimArrayMap (Sparse w v) sh e e' where + {-# Inline mapArray #-} + mapArray f sparse = sparse{sparseData = VG.map f (sparseData sparse)} + + + + + +-- * Helper functions. + +-- | Find the index with manhattan helper +-- +-- TODO consider using binary search instead of a linear scan here! +-- e.g.: @k = VAS.binarySearchByBounds (==)@ +-- +-- NOTE running times with 100x100 DP problem "NeedlemanWunsch" +-- full findIndex of sixs: 0,050,000 cells/sec +-- using manhattan buckets, findIndex: 5,000,000 cells/sec +-- using binarySearch on slices: 11,000,000 cells/sec +-- +-- On a 1000x1000 DP NeedlemanWunsch problem, binary search on slices is at 6,500,000 cells/sec. + +manhattanIndex + :: (SparseBucket sh, Index sh) + => LimitType sh -> Vector Int -> VU.Vector Int -> sh -> Maybe Int +{-# Inline manhattanIndex #-} +manhattanIndex ub mstart sixs idx = fmap (+l) . binarySearch (linearIndex ub idx) $ VG.unsafeSlice l (h-l) sixs + where + b = manhattan ub idx + -- lower and upper bucket bounds + l = mstart `VU.unsafeIndex` b + h = mstart `VU.unsafeIndex` (b+1) + +binarySearch :: Int -> VU.Vector Int -> Maybe Int +{-# Inline binarySearch #-} +{- +binarySearch k xs = + let r1 = binarySearch1 k xs + r2 = binarySearch2 k xs + in if r1==r2 then r1 else error $ show (k,xs,r1,r2) +-} +{- +-- 1000x1000 at @1000 yields 3,050,000 cells / second +binarySearch (I# e) v = go 0 pp + where + -- largest index to check + (I# r) = VG.length v -1 + -- largest power of two <= (r+1) + pp = (2 ^ (max 0 $ msb $ VG.length v -1)) + -- wrap the actual non-branching worker function + go :: Int -> Int -> Maybe Int + {-# Inline go #-} + go (I# l) (I# p) = let i = I# (go' l p) in if (VG.length v<1 || i<0) then Nothing else Just i + -- @go'@ should be non-branching, and use a minimal number of array reads. + go' :: Int# -> Int# -> Int# + {-# Inline go' #-} + go' l p + -- we are done and will return the proposed position of the last element found or -1 + | 1# <- p ==# 0# = (e ==# x) *# l -# (e /=# x) + | otherwise = let i = go' (l +# (p *# chk *# leq)) (quotInt# p 2#) -- (uncheckedIShiftRA# p 1#) + -- (I# ii) = traceShow (I# l, I# p, I# i2r, I# x, I# leq) (I# i) + in i -- i + where i2r = l +# (p *# chk) -- index to read + (I# x) = VU.unsafeIndex v (I# i2r) + leq = x <=# e + newl = l +# p + chk = newl <=# r +-} +-- +-- 1000x1000 at @1000 yields 6,030,000 cells / second +binarySearch e v = go 0 (VG.length v -1) + where + go :: Int -> Int -> Maybe Int + go !l !r = + let !m = (r+l) `div` 2 + !x = VG.unsafeIndex v m + in if r go l (m-1) + EQ -> Just m + GT -> go (m+1) r +-- + + +-- | Given two index vectors of the same shape, will return the correctly ordered vector of +-- the union of indices. +-- +-- TODO This requires that @Ord (Shape O)@ uses the @Down@ instance of Ord! We need to fix this in +-- the @Index@ modules. +-- +-- TODO Rewrite to allow fusion without intermediate vectors using uncons. This will make it +-- possible to chain applications. @stream@ should be fine for this. + +mergeIndexVectors :: (Eq sh, Ord sh, VG.Vector w sh) => w sh -> w sh -> w sh +{-# Inlinable mergeIndexVectors #-} +mergeIndexVectors xs ys = VG.create $ do + let lxs = VG.length xs + lys = VG.length ys + mv <- VGM.new $ lxs + lys + let go !n !i !j + | i>=lxs && j>=lys = return n + | j>=lys = VG.unsafeIndexM xs i >>= VGM.unsafeWrite mv n >> go (n+1) (i+1) j + | i>=lxs = VG.unsafeIndexM ys j >>= VGM.unsafeWrite mv n >> go (n+1) i (j+1) + | otherwise = do + x <- VG.unsafeIndexM xs i + y <- VG.unsafeIndexM ys j + if | x==y -> VGM.unsafeWrite mv n x >> go (n+1) (i+1) (j+1) + | x< y -> VGM.unsafeWrite mv n x >> go (n+1) (i+1) j + | x> y -> VGM.unsafeWrite mv n y >> go (n+1) i (j+1) + n <- go 0 0 0 + return $ VGM.unsafeTake n mv + diff --git a/Data/PrimitiveArray/Zero.hs b/Data/PrimitiveArray/Zero.hs deleted file mode 100644 --- a/Data/PrimitiveArray/Zero.hs +++ /dev/null @@ -1,71 +0,0 @@ -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE TypeFamilies #-} - --- | Boxed, primitive arrays. A good use-case is to store boxed or unboxed --- vectors. - -module Data.PrimitiveArray.Zero where - -import Control.Monad -import Data.Array.Repa.Index -import Data.Array.Repa.Shape -import Data.Primitive -import Data.Primitive.Array -import Data.Primitive.Types -import Control.Exception (assert) - -import Data.ExtShape -import Data.PrimitiveArray - - - --- | Monadic arrays of primitive type. - -data MArr0 s sh elm = MArr0 !sh !(MutableArray s elm) - --- | Immutable arrays of primitive type. - -data Arr0 sh elm = Arr0 !sh !(Array elm) - - - -type instance MutArray Arr0 = MArr0 - --- NOTE inLb, inUb is including bound, while exUb is excluding upper bound. --- Differentiates between largest included index, first excluded index. - -instance (Shape sh, ExtShape sh) => MPrimArrayOps MArr0 sh elm where - boundsM (MArr0 exUb _) = (zeroDim,exUb `subDim` unitDim) - fromListM inLb inUb xs = do - ma <- newM inLb inUb - let exUb = inUb `addDim` unitDim - let (MArr0 _ mba) = ma - zipWithM_ (\k x -> assert (length xs == size exUb) $ writeArray mba k x) [0.. toIndex exUb inUb] xs - return ma - newM inLb inUb = let exUb = inUb `addDim` unitDim in - unless (inLb == zeroDim) (error "MArr0 lb/=zeroDim") >> - MArr0 exUb `liftM` newArray (size exUb) undefined - newWithM inLb inUb def = do - let exUb = inUb `addDim` unitDim - ma <- newM inLb inUb - let (MArr0 _ mba) = ma - forM_ [0 .. toIndex exUb inUb] $ \k -> writeArray mba k def - return ma - readM (MArr0 exUb mba) idx = assert (inShape exUb idx) $ readArray mba (toIndex exUb idx) - writeM (MArr0 exUb mba) idx elm = assert (inShape exUb idx) $ writeArray mba (toIndex exUb idx) elm - {-# INLINE boundsM #-} - {-# INLINE fromListM #-} - {-# INLINE newM #-} - {-# INLINE newWithM #-} - {-# INLINE readM #-} - {-# INLINE writeM #-} - -instance (Shape sh, ExtShape sh) => PrimArrayOps Arr0 sh elm where - bounds (Arr0 exUb _) = (zeroDim,exUb `subDim` unitDim) - freeze (MArr0 exUb mba) = Arr0 exUb `liftM` unsafeFreezeArray mba - index (Arr0 exUb ba) idx = assert (inShape exUb idx) $ indexArray ba (toIndex exUb idx) - {-# INLINE bounds #-} - {-# INLINE freeze #-} - {-# INLINE index #-} - diff --git a/Data/PrimitiveArray/Zero/Unboxed.hs b/Data/PrimitiveArray/Zero/Unboxed.hs deleted file mode 100644 --- a/Data/PrimitiveArray/Zero/Unboxed.hs +++ /dev/null @@ -1,71 +0,0 @@ -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE FlexibleInstances #-} - --- | Strict, unboxed arrays of primitive type. Uses unboxed vectors internally --- to provide tuple instances. - -module Data.PrimitiveArray.Zero.Unboxed where - -import Control.Monad -import Data.Array.Repa.Index -import Data.Array.Repa.Shape -import Control.Exception (assert) -import Data.Vector.Unboxed as VU hiding (forM_, length, zipWithM_) -import Data.Vector.Unboxed.Mutable as VUM hiding (length) - -import Data.ExtShape -import Data.PrimitiveArray - - - --- | Monadic arrays of primitive type. - -data MArr0 s sh elm = MArr0 !sh !(MVector s elm) - --- | Immutable arrays of primitive type. - -data Arr0 sh elm = Arr0 !sh !(Vector elm) - - - -type instance MutArray Arr0 = MArr0 - --- NOTE inLb, inUb is including bound, while exUb is excluding upper bound. --- Differentiates between largest included index, first excluded index. - -instance (Shape sh, ExtShape sh, VU.Unbox elm) => MPrimArrayOps MArr0 sh elm where - boundsM (MArr0 exUb _) = (zeroDim,exUb `subDim` unitDim) - fromListM inLb inUb xs = do - ma <- newM inLb inUb - let exUb = inUb `addDim` unitDim - let (MArr0 _ mba) = ma - zipWithM_ (\k x -> assert (length xs == size exUb) $ unsafeWrite mba k x) [0.. toIndex exUb inUb] xs - return ma - newM inLb inUb = let exUb = inUb `addDim` unitDim in - unless (inLb == zeroDim) (error "MArr0 lb/=zeroDim") >> - MArr0 exUb `liftM` new (size exUb) - newWithM inLb inUb def = do - let exUb = inUb `addDim` unitDim - ma <- newM inLb inUb - let (MArr0 _ mba) = ma - forM_ [0 .. toIndex exUb inUb] $ \k -> unsafeWrite mba k def - return ma - readM (MArr0 exUb mba) idx = assert (inShape exUb idx) $ unsafeRead mba (toIndex exUb idx) - writeM (MArr0 exUb mba) idx elm = assert (inShape exUb idx) $ unsafeWrite mba (toIndex exUb idx) elm - {-# INLINE boundsM #-} - {-# INLINE fromListM #-} - {-# INLINE newM #-} - {-# INLINE newWithM #-} - {-# INLINE readM #-} - {-# INLINE writeM #-} - -instance (Shape sh, ExtShape sh, VUM.Unbox elm) => PrimArrayOps Arr0 sh elm where - bounds (Arr0 exUb _) = (zeroDim,exUb `subDim` unitDim) - freeze (MArr0 exUb mba) = Arr0 exUb `liftM` unsafeFreeze mba - index (Arr0 exUb ba) idx = assert (inShape exUb idx) $ unsafeIndex ba (toIndex exUb idx) - {-# INLINE bounds #-} - {-# INLINE freeze #-} - {-# INLINE index #-} - diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright Christian Hoener zu Siederdissen 2010 +Copyright Christian Hoener zu Siederdissen 2010-2015 All rights reserved. diff --git a/PrimitiveArray.cabal b/PrimitiveArray.cabal --- a/PrimitiveArray.cabal +++ b/PrimitiveArray.cabal @@ -1,49 +1,216 @@ +Cabal-version: 2.2 Name: PrimitiveArray -Version: 0.4.0.0 -License: BSD3 +Version: 0.10.1.1 +License: BSD-3-Clause License-file: LICENSE -Author: Christian Hoener zu Siederdissen -Maintainer: choener@tbi.univie.ac.at -Copyright: Christian Hoener zu Siederdissen, 2010-2012 -Homepage: http://www.tbi.univie.ac.at/~choener/ +Maintainer: choener@bioinf.uni-leipzig.de +author: Christian Hoener zu Siederdissen, 2011-2021 +copyright: Christian Hoener zu Siederdissen, 2011-2021 +homepage: https://github.com/choener/PrimitiveArray +bug-reports: https://github.com/choener/PrimitiveArray/issues Stability: Experimental Category: Data Build-type: Simple -Cabal-version: >=1.6 -Synopsis: - Efficient multidimensional arrays +tested-with: GHC == 8.8, GHC == 8.10, GHC == 9.0 +Synopsis: Efficient multidimensional arrays Description: - This library provides efficient multidimensional arrays. + . - In general all operations are (highly) unsafe, no - bounds-checking or other sanity-checking is performed. - Operations are aimed toward efficiency as much as possible. - Goals of the library are to have arrays according to three - ideas: immutable/mutable arrays, strict/lazy arrays, - zero-based/lower-bound arrays. Zero-based arrays save one - addition on each access if the lower bound or the array is - always zero. + This library provides efficient multidimensional arrays. Import + @Data.PrimitiveArray@ for indices, lenses, and arrays. . - We have forked two repa modules: Shape and Index. + For + + users, the library also provides the machinary to + fill tables in the correct order required by usual CYK-style + parsers, or regular grammars (used e.g. in alignment + algorithms). This means that unless your grammar requires a + strange order in which parsing is to be performed, it will + mostly "just work". + . + In general operations do not perform bounds-checking or other + sanity-checking and are aimed towards efficiency as much as + possible. Users (like @ADPfusion@) should perform their own + bounds-checking, outside of code that performs "loop-like" + operations. + + +extra-source-files: + README.md + changelog.md + + + +flag debug + description: Enable bounds checking and various other debug operations at the cost of a significant performance penalty. + default: False + manual: True + +flag debugoutput + description: Enable debug output, which spams the screen full of index information + default: False + manual: True + +flag llvm + description: use llvm + default: False + manual: True + +flag debugdump + description: Enable dumping intermediate / core files + default: False + manual: True + +flag dump-core + description: Dump HTML for the core generated by GHC during compilation + default: False + manual: True + + + +common deps + build-depends: base >= 4.7 && < 5.0 + , aeson >= 0.8 + , binary >= 0.7 + , bits >= 0.4 + , cereal >= 0.4 + , cereal-vector >= 0.2 + , containers + , deepseq >= 1.3 + , hashable >= 1.2 + , hashtables >= 1.2 + , lens >= 4.0 + , log-domain >= 0.10 + , mtl >= 2.0 + , primitive >= 0.5.4 + , QuickCheck >= 2.7 + , smallcheck >= 1.1 + , text >= 1.0 + , unordered-containers >= 0.2 + , vector >= 0.11 + , vector-algorithms >= 0.8 + , vector-binary-instances >= 0.2 + , vector-th-unbox >= 0.2 + -- + , DPutils == 0.1.1.* + , OrderedBits == 0.0.2.* + default-extensions: BangPatterns + , CPP + , DataKinds + , DefaultSignatures + , DeriveDataTypeable + , DeriveFunctor + , DeriveGeneric + , DerivingStrategies + , FlexibleContexts + , FlexibleInstances + , FunctionalDependencies + , GADTs + , GeneralizedNewtypeDeriving + , MultiParamTypeClasses + , MultiWayIf + , PatternSynonyms + , PolyKinds + , RankNTypes + , RecordWildCards + , ScopedTypeVariables + , StandaloneDeriving + , TemplateHaskell + , TypeApplications + , TypeFamilies + , TypeOperators + , UndecidableInstances + , UnicodeSyntax + default-language: + Haskell2010 + ghc-options: + -O2 + -funbox-strict-fields + if flag(debug) + cpp-options: -DADPFUSION_CHECKS + ghc-options: -fno-ignore-asserts --disable-optimizations + if flag(debugoutput) + cpp-options: -DADPFUSION_DEBUGOUTPUT + + + Library + import: + deps Exposed-modules: - Data.Array.Repa.Index - Data.Array.Repa.Shape - Data.ExtShape Data.PrimitiveArray - Data.PrimitiveArray.Zero - Data.PrimitiveArray.Zero.Unboxed --- Data.PrimitiveArray.UpperTriangular --- Data.PrimitiveArray.UpperTriangular.Unboxed --- Data.PrimitiveArray.FillTable - Build-depends: - base >= 4 && <5, - primitive == 0.5 , - vector == 0.10 + Data.PrimitiveArray.Checked + Data.PrimitiveArray.Class + Data.PrimitiveArray.Dense + Data.PrimitiveArray.Sparse + Data.PrimitiveArray.Sparse.BinSearch + Data.PrimitiveArray.Sparse.IntBinSearch + Data.PrimitiveArray.HashTable + Data.PrimitiveArray.Index + Data.PrimitiveArray.Index.BitSet0 + Data.PrimitiveArray.Index.BitSet1 + Data.PrimitiveArray.Index.BitSetClasses + Data.PrimitiveArray.Index.Class + Data.PrimitiveArray.Index.Int + Data.PrimitiveArray.Index.IOC + Data.PrimitiveArray.Index.PhantomInt + Data.PrimitiveArray.Index.Point + Data.PrimitiveArray.Index.Subword + Data.PrimitiveArray.Index.Unit + + + +test-suite properties + import: + deps + build-depends: base + , tasty >= 0.11 + , tasty-quickcheck >= 0.8 + , tasty-smallcheck >= 0.8 + , tasty-th >= 0.1 + -- + , PrimitiveArray + type: + exitcode-stdio-1.0 + main-is: + properties.hs + other-modules: + QuickCheck + SmallCheck + Common ghc-options: - -O2 - -funbox-strict-fields + -threaded -rtsopts -with-rtsopts=-N + hs-source-dirs: + tests + build-depends: PrimitiveArray + + +benchmark Lookup + import: + deps + type: + exitcode-stdio-1.0 + main-is: + Lookup.hs + hs-source-dirs: + bench + build-depends: PrimitiveArray + , criterion ^>= 1.5 + if flag(llvm) + ghc-options: + -fllvm + -optlo-O3 + if flag(debugdump) + ghc-options: + -ddump-to-file + -ddump-simpl + -dsuppress-all + if flag(dump-core) + build-depends: dump-core + ghc-options: -fplugin=DumpCore -fplugin-opt DumpCore:core-html + + source-repository head diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,22 @@ +![github action: CI](https://github.com/choener/PrimitiveArray/actions/workflows/ci.yml/badge.svg) +![github action: hackage](https://github.com/choener/PrimitiveArray/actions/workflows/hackage.yml/badge.svg) + +# PrimitiveArray + +[*generalized Algebraic Dynamic Programming Homepage*](http://www.bioinf.uni-leipzig.de/Software/gADP/) + +PrimitiveArray provides operations on multi-dimensional arrays. Internally, the +representation is based on the vector library, while the multi-dimensional +indexing follows repa. + +Primitive arrays are designed to be used together with ADPfusion. + + + +#### Contact + +Christian Hoener zu Siederdissen +Leipzig University, Leipzig, Germany +choener@bioinf.uni-leipzig.de +http://www.bioinf.uni-leipzig.de/~choener/ + diff --git a/bench/Lookup.hs b/bench/Lookup.hs new file mode 100644 --- /dev/null +++ b/bench/Lookup.hs @@ -0,0 +1,50 @@ + +module Main where + +import Criterion.Main + +import Data.PrimitiveArray as PA + + + +go ∷ (Index i) ⇒ Int → Unboxed i Int → i → Int +{-# Inline go #-} +go !c !pa !i = f c i 0 + where f 0 !i !acc = acc + f !k !i !acc = f (k-1) i (acc + pa ! i) + +go1 ∷ Int → Unboxed (Z:.Int) Int → (Z:.Int) → Int +{-# NoInline go1 #-} +go1 = go + +go2 ∷ Int → Unboxed (Z:.Int:.Int) Int → (Z:.Int:.Int) → Int +{-# NoInline go2 #-} +go2 = go + +go3 ∷ Int → Unboxed (Z:.Int:.Int:.Int) Int → (Z:.Int:.Int:.Int) → Int +{-# NoInline go3 #-} +go3 = go + +main ∷ IO () +main = do + let !(pa1 ∷ Unboxed (Z:.Int) Int) = PA.fromAssocs (ZZ:..LtInt 10) 0 [] + let !(pa2 ∷ Unboxed (Z:.Int:.Int) Int) = PA.fromAssocs (ZZ:..LtInt 10:..LtInt 10) 0 [] + let !(pa3 ∷ Unboxed (Z:.Int:.Int:.Int) Int) = PA.fromAssocs (ZZ:..LtInt 10:..LtInt 10:..LtInt 10) 0 [] + defaultMain + [ bgroup "1" + [ bench "10^0" $ whnf (go1 1 pa1) (Z:.5) + , bench "10^3" $ whnf (go1 1000 pa1) (Z:.5) + , bench "10^6" $ whnf (go1 1000000 pa1) (Z:.5) + , bench "10^9" $ whnf (go1 1000000000 pa1) (Z:.5) + ] + , bgroup "2" + [ bench " 1" $ whnf (go2 1 pa2) (Z:.5:.5) + , bench " 1000" $ whnf (go2 1000 pa2) (Z:.5:.5) + , bench "1000000" $ whnf (go2 1000000 pa2) (Z:.5:.5) + ] + , bgroup "3" + [ bench " 1" $ whnf (go3 1 pa3) (Z:.5:.5:.5) + , bench " 1000" $ whnf (go3 1000 pa3) (Z:.5:.5:.5) + , bench "1000000" $ whnf (go3 1000000 pa3) (Z:.5:.5:.5) + ] + ] diff --git a/changelog.md b/changelog.md new file mode 100644 --- /dev/null +++ b/changelog.md @@ -0,0 +1,132 @@ +0.11.1.1 + +- version bump on DPutils + +0.10.1.0 +-------- + +- introduction of @Data.PrimitiveArray.Sparse@ which uses different sparsification options. The + default is @D.P.S.Search@ based on binary search. +- All array operations, pure or mutable are now based on a single, unified class. Mostly because + mutable operations go via a data family anyway. + +0.10.0.0 +-------- + +- Rewrote Data.PrimitiveArray.Dense to accept all vector types using one + interface. This is a breaking change, since @Unboxed@ becomes @Dense + Data.Vector.Unboxed.Vector@, but now @Dense v@ accepts any @v@ as underlying + storage vector. Breaking occurs only at user sites where the actual vector + type needs to be specified. This tends to be very localized. + +0.9.1.1 +------- + +- OrderedBits version bump + +0.9.1.0 +------- + +- Arbitrary instance(s), field lenses that are probably not a good idea (don't use them!) + +0.9.0.0 +------- + +- large-scale changes +- associated data families for bounds + +0.8.1.0 +------- + +- inclusion of Upperlimit data family to simplify declaration of upper limits + +0.8.0.1 +------- + +- PointL delays inlining to phase 0 for table filling. This is part of the + close-to-C optimization effort for linear languages. +- disabling smallcheck until I fix how things are generated + +0.8.0.0 +------- + +- renamed Interface (Iter) to Boundary (Boundary) +- Typeable instances for Dense primitive arrays +- EdgeBoundary index structure +- changes and fixes to quickcheck/smallcheck +- added ScoreMatrix module with simple score and distance matrix structure + (requires log-domain) + +0.7.2.0 +------- + +- JSONKey (To/From) for index types. + +0.7.1.0 +------- + +- minor updates to dependencies +- tasty framework +- Subword/upper triangular indexing provided by DPutils + +0.7.0.1 +------- + +- Data.PrimitiveArray.Checked to capture index out-of-bounds problems + +0.7.0.0 +------- + +- vector <= 0.11 support; including compatibility layer +- redesigned Index structures (for dealing with Inside/Outside/Complement) + +0.6.1.1 +------- + +- Hashable instances for all index structures +- Hashable instances for Unboxed and Boxed arrays. *These require Hashable + instances for vectors, which are not available by default* + +0.6.1.0 +------- + +- OrderedBits < 0.0.1 +- travis.yml update + +0.6.0.0 +------- + +- moved primitive array classes to Data.PrimitiveArray.Class +- added from / to lenses +- Field1 .. Field6 lenses for indices (Z:.a:.b...) (with Z being Field0) + - lens stuff currently commented out; aiming to have an extra package [lens + is fairly heavy] +- FillTables should work now (with PointL, Subword) +- freezing of whole stacks of (Z:.mutarr:.mutarr:. ...) tables +- explicit 'Shape Subword'; this allows for simpler code in a number of places + and is especially useful for CYK-style algorithms that have a + single-dimensional upper-triangular matrix. +- rangeStream of Extshape is new and used by the FillTables module +- Binary, Cereal, Aeson instances for indices and immutable tables +- orphan instances of Binary, Cereal, Aeson for Z, and (:.) +- topmostIndex returns the final index position for CYK-style (bottom to top) + parsing +- removed Data.Array.Repa.Index.Point (we have PointL, PointR in Points.hs) +- added Data.Array.Repa.Index.Set (for sets with an interface, used by + Hamiltonian path problems) +- Data.Array.Repa.Index.Outside is now just a newtype wrapped around other + Index types. We want to be able to say "a Subword, but for Outside + algorithms" +- travis-ci integration + +0.5.4.0 +------- + +- actually implemented PointR + +- added the rather important strictness annotation for mutable arrays in .Zero + +0.5.3.0 +------- + +- fixed vector-th-unbox problem diff --git a/tests/Common.hs b/tests/Common.hs new file mode 100644 --- /dev/null +++ b/tests/Common.hs @@ -0,0 +1,40 @@ + +module Common where + +import Control.Arrow ((&&&), second) +import Data.List (nub, sort, group) +import qualified Data.Map.Strict as M +import qualified Data.Set as S + +import Data.PrimitiveArray.Index.Class + + +-- * generic functions + +-- | Generates a list of, eg, @PointL@s. This are then grouped according to +-- the @linearIndex@. Within each group, there should only be @PointL@s +-- with the same value. + +uniquenessTest ∷ (Ord a, Index a) ⇒ LimitType a → LimitType a → [a] → Bool +uniquenessTest low hgh xs = all allEq ys && all allEq zs + where ys = M.fromListWith S.union . map (second S.singleton) . map (linearIndex hgh &&& id) $ xs + zs = M.fromListWith S.union . map (second S.singleton) . map (id &&& linearIndex hgh) $ xs +{-# Inlineable uniquenessTest #-} +{- +uniquenessTest low xs = all allEq ys && all allEq zs + where hgh = maximum xs + ys = group . sort . map (linearIndex low hgh &&& id) $ xs + zs = group . sort . map (id &&& linearIndex low hgh) $ xs +-} + +-- | are all @xs@ equal to each other + +allEq = (1==) . S.size +{-# Inline allEq #-} + +{- +allEq [] = True +allEq (x:xs) = all (x==) xs +-} + + diff --git a/tests/QuickCheck.hs b/tests/QuickCheck.hs new file mode 100644 --- /dev/null +++ b/tests/QuickCheck.hs @@ -0,0 +1,36 @@ + +module QuickCheck where + +import Data.List (maximumBy) +import Data.Ord (comparing) +import Test.QuickCheck +import Test.Tasty.QuickCheck +import Test.Tasty.TH + +import Data.PrimitiveArray.Index.Class +--import Data.PrimitiveArray.Index.EdgeBoundary +import Data.PrimitiveArray.Index.IOC +import Data.PrimitiveArray.Index.Point +--import Data.PrimitiveArray.Index.Set +--import Data.PrimitiveArray.Index.Subword + +import Common + + + +-- * Uniqueness tests + +-- prop_PointL_I_unique (xs :: [PointL I]) = uniquenessTest (LtPointL 0) (LtPointL $ maximum $ map fromPointL xs) xs + +-- prop_Subword_I_unique (xs :: [Subword I]) = uniquenessTest (subword 0 0) (maximumBy (comparing fromSubwordSnd) xs) xs + +-- prop_EdgeBoundary_I_unique (xs :: [EdgeBoundary I]) = uniquenessTest (0 :-> 0) (maximumBy (comparing fromEdgeBoundarySnd) xs) xs + +-- | TODO check that bitsets produce the correct number of bits when counting + +--prop_BitSet1_First_I_set (numberOfBits ∷ ()) = strm == lst +-- where strm = sort . unId $ streamUp (LtBitSet1 0) (LtBitSet1 0) :: IO [BitSet1 First I] +-- lst = sort [] + +quickcheck_tests = $(testGroupGenerator) + diff --git a/tests/SmallCheck.hs b/tests/SmallCheck.hs new file mode 100644 --- /dev/null +++ b/tests/SmallCheck.hs @@ -0,0 +1,37 @@ + +module SmallCheck where + +import Control.Applicative +import Data.Bits +import Data.List (nub, sort, group, maximumBy) +import Data.Ord (comparing) +import Data.Word (Word) +import Debug.Trace +import Test.SmallCheck +import Test.Tasty +import Test.Tasty.SmallCheck +import Test.Tasty.TH + +import Data.PrimitiveArray.Index.Class +--import Data.PrimitiveArray.Index.EdgeBoundary +import Data.PrimitiveArray.Index.IOC +import Data.PrimitiveArray.Index.Point +--import Data.PrimitiveArray.Index.Set +--import Data.PrimitiveArray.Index.Subword + +import Common + + + +-- * Uniqueness tests. The @xs@ lists are fairly small. + +prop_PointL_I_unique (xs :: [PointL I]) = uniquenessTest (LtPointL 0) (LtPointL $ maximum $ map fromPointL xs) xs + +-- prop_Subword_I_unique (xs :: [Subword I]) = uniquenessTest (subword 0 0) (maximumBy (comparing fromSubwordSnd) xs) xs + +-- prop_EdgeBoundary_I_unique (xs :: [EdgeBoundary I]) = uniquenessTest (0 :-> 0) (maximumBy (comparing fromEdgeBoundarySnd) xs) xs + + + +smallcheck_tests = $(testGroupGenerator) + diff --git a/tests/properties.hs b/tests/properties.hs new file mode 100644 --- /dev/null +++ b/tests/properties.hs @@ -0,0 +1,73 @@ + +module Main where + +import Control.Applicative +import Data.Bits +import Data.List (nub, sort, group) +import Data.Word (Word) +import Test.Tasty +import Test.Tasty.TH +import qualified Test.QuickCheck as QC + +import Data.PrimitiveArray.Index.IOC +import Data.PrimitiveArray.Index.Point +--import Data.PrimitiveArray.Index.Set +import Data.PrimitiveArray.Index.Class + +import QuickCheck +import SmallCheck + + + +-- * Points + +-- | @linearIndex <-> fromLinearIndex@ + +prop_FromLinear_ZP ( x :: PointL I, a') + | ix == frm = True + | otherwise = error $ show (x,a',lt, ix, lin, frm) + where ltx = LtPointL $ QC.getNonNegative a' + fromPointL x + lt = ZZ:..ltx + ix = Z:.x + lin = linearIndex lt ix + frm = fromLinearIndex lt lin + +prop_FromLinear_ZPP ( x :: PointL I, y :: PointL I, a', b') + | ix == frm = True + | otherwise = error $ show (x,y,a',b',lt, ix, lin, frm) + where ltx = LtPointL $ QC.getNonNegative a' + fromPointL x + lty = LtPointL $ QC.getNonNegative b' + fromPointL y + lt = ZZ:..ltx:..lty + ix = Z:.x:.y + lin = linearIndex lt ix + frm = fromLinearIndex lt lin + +-- * Sets + +-- TODO what exactly does the mask fix? Only bits already @1@, or every bit +-- as it is? The mask should actually freeze-fix those bits, where we are +-- set to @1@! + +--prop_Fixed_BitSet_setSucc (u :: Word, Fixed m s :: Fixed (BitSet I)) = traceShow (tgo, tsu) $ tgo == tsu +-- where tgo = go s +-- tsu = (getFixed <$> setSucc (Fixed 0 0) (Fixed 0 h) (Fixed m s)) +-- fb1 = m .&. s -- fixed bits to 1 +-- fb0 = m .&. complement s -- fixed bits to 0 +-- h = bit (fromIntegral $ u `mod` 8) - 1 +-- go x -- continue creating successors, until the mask criterion is met (again). +-- | Nothing <- ssx = Nothing +-- | Just x' <- ssx +-- , fb0 == m .&. complement x' +-- , fb1 == m .&. x' = traceShow ('j',fb0,fb1,m,x,x') $ Just x' +-- | Just x' <- ssx = traceShow ('g',fb0,fb1,m,x,x') $ go x' +-- where ssx = setSucc 0 h x + + + +main :: IO () +main = do + defaultMain $ testGroup "" + [ -- quickcheck_tests +-- , smallcheck_tests + ] +