diff --git a/DPutils.cabal b/DPutils.cabal new file mode 100644 --- /dev/null +++ b/DPutils.cabal @@ -0,0 +1,102 @@ +Name: DPutils +Version: 0.0.0.2 +License: BSD3 +License-file: LICENSE +Maintainer: choener@bioinf.uni-leipzig.de +author: Christian Hoener zu Siederdissen, 2016 +copyright: Christian Hoener zu Siederdissen, 2016 +homepage: https://github.com/choener/DPutils +bug-reports: https://github.com/choener/DPutils/issues +Stability: Experimental +Category: Data +Build-type: Simple +Cabal-version: >=1.10.0 +tested-with: GHC == 7.8.4, GHC == 7.10.3, GHC == 8.0.1 +Synopsis: utilities for DP +Description: + Small set of utility functions + . + + + +extra-source-files: + README.md + changelog.md + + + +Library + Exposed-modules: + Data.Paired.Common + Data.Paired.Foldable + Data.Paired.Vector + Math.TriangularNumbers + Pipes.Parallel + build-depends: base >= 4.7 && < 5.0 + , containers + , kan-extensions >= 4.0 && < 6.0 + , parallel >= 3.0 && < 4.0 + , pipes >= 4.0 && < 4.3 + , QuickCheck >= 2.7 + , vector >= 0.10 && < 0.12 + default-extensions: BangPatterns + , CPP + , ScopedTypeVariables + , FlexibleContexts + default-language: + Haskell2010 + ghc-options: + -O2 + -funbox-strict-fields + + + +test-suite properties + type: + exitcode-stdio-1.0 + main-is: + properties.hs + ghc-options: + -O2 -threaded -rtsopts -with-rtsopts=-N + hs-source-dirs: + tests + default-language: + Haskell2010 + default-extensions: CPP + , TemplateHaskell + build-depends: base + , containers + , DPutils + , kan-extensions + , parallel + , pipes + , QuickCheck + , tasty >= 0.11 + , tasty-quickcheck >= 0.8 + , tasty-th >= 0.1 + , vector + + + +benchmark benchmark + type: + exitcode-stdio-1.0 + build-depends: base + , criterion >= 1.1 + , DPutils + , vector + hs-source-dirs: + tests + main-is: + benchmark.hs + default-language: + Haskell2010 + ghc-options: + -O2 + + + +source-repository head + type: git + location: git://github.com/choener/DPutils + diff --git a/Data/Paired/Common.hs b/Data/Paired/Common.hs new file mode 100644 --- /dev/null +++ b/Data/Paired/Common.hs @@ -0,0 +1,29 @@ + +module Data.Paired.Common where + + + +-- | Shall we combine elements on the main diagonal as well? +-- +-- If we choose @NoDiag@, we deal with upper triangular matrices that are +-- effectively one element smaller. + +data OnDiag = OnDiag | NoDiag + deriving (Eq) + +-- | Select only a subset of the possible enumerations. + +data Enumerate + -- | Enumerate all elements + = All + -- | Enumerate from a value and at most @N@ elements + | FromN Int Int + deriving (Eq) + +-- | If the size of the input is known before-hand or not. + +data SizeHint + = UnknownSize + | KnownSize Int + deriving (Eq) + diff --git a/Data/Paired/Foldable.hs b/Data/Paired/Foldable.hs new file mode 100644 --- /dev/null +++ b/Data/Paired/Foldable.hs @@ -0,0 +1,121 @@ + +-- | Efficient enumeration of subsets of triangular elements. Given a list +-- @[1..n]@ we want to enumerate a subset @[(i,j)]@ of ordered pairs in +-- such a way that we only have to hold the elements necessary for this +-- subset in memory. + +module Data.Paired.Foldable where + +import Data.IntMap as IM +import Data.Foldable as F +import Data.List as L +import Control.Arrow ((***)) +import Data.Vector as V +import Data.Vector.Generic as VG +import Debug.Trace (traceShow) +import Text.Printf + +import Data.Paired.Common +import Math.TriangularNumbers + + + +-- | Generalized upper triangular elements. Given a list of elements +-- @[e_1,...,e_k]@, we want to return pairs @(e_i,e_j)@ such that we have +-- all ordered pairs with @i SizeHint + -- ^ If the size of @t a@ is known beforehand, give the appropriate + -- @KnownSize n@, otherwise give @UnknownSize@. Using @UnknownSize@ will + -- force the complete spine of @t a@. + -> OnDiag + -- ^ The enumeration will include the pairs on the main diagonal with + -- @OnDiag@, meaning @(i,i)@ will be included for all @i@. Otherwise, + -- @NoDiag@ will exclude these elements. + -> Enumerate + -- ^ Either enumerate @All@ elements or enumerate the @s@ elements + -- starting at @k@ with @FromN k s@. + -> t a + -- ^ The foldable data structure to enumerate over. + -> Either String (IntMap a, Int, [((Int,Int),(a,a))]) + -- ^ If there is any error then return @Left errorMsg@. Otherwise we have + -- @Right (imap, numElems, list)@. The @imap@ structure holds the subset + -- of elements with which we actually generate elements. @numElems@ is + -- the total number of elements that will be generated. This is + -- calculated without touch @list@. Finally, @list@ is the lazy list of + -- elements to be generated. +upperTri sz d e xs + | szLen /= readLen = Left $ printf "Expected SizeHint %d elements, but processed only %d elements!" szLen readLen + | otherwise = Right (imp, numElems, ys) + where ys = case e of {All -> id ; FromN _ s -> L.take s} + . L.unfoldr go $ initEnum e d + -- how many elements we will emit depends on enumeration and on + -- diagonal element counting + numElems + | All <- e = allSize + | FromN s k <- e = if s+k > allSize then max 0 (allSize - s) else k + -- The length of the input. With a given size hint, @xs :: t a@ + -- will only be touched once. +#if MIN_VERSION_base(4,8,0) + szLen = case sz of { UnknownSize -> F.length xs ; KnownSize z -> z } +#else + szLen = case sz of { UnknownSize -> L.length . F.toList $ xs ; KnownSize z -> z } +#endif + szLn' = case d of { OnDiag -> szLen - 1 ; NoDiag -> szLen - 2 } + -- Construct an intmap @imp@ of all elements in the accepted range. + -- At the same time, return the length or size of the foldable + -- container we gave as input. @xs@ is touched only once and can + -- be efficiently consumed. + (!imp,!readLen) = F.foldl' (\(!i,!l) x -> (if inRange l then IM.insert l x i else i,l+1)) (IM.empty, 0) xs + allSize = szLen * (szLen + if d == OnDiag then 1 else -1) `div` 2 + -- we need three ranges. @cMin@ and @cMax@ are the range for the + -- slow-moving first element in the tuple. @rMin@ and @rMax@ are + -- the first and last element of the range starting at @cMin@ (we + -- can actually start at @cMax@ but it doesn't matter). + -- Finally, @lMin@ and @lMax@ are the range to the left of @cMin@. + (lMin,lMax,cMin,cMax,rMin,rMax) = case e of + All -> (0, szLen-1, 0, szLen-1, 0, szLen-1) + FromN s k -> + let (cmin,rmin) = fromLinear szLn' s + (cmax,_ ) = fromLinear szLn' (s+k) + rmax = rmin+k -- if this is @>= len@ we are safe anyway. + lmin = if rmin+k >= szLen then 0 else cmin + lmax = if rmin+k >= szLen then lmin + toLinear szLn' (cmin+1,cmin+1+rmin+k-szLn') else cmax + in (lmin, lmax, cmin, cmax, rmin, rmax) + -- Determine if an element at linear index @z@ is in the range to + -- be consumed. + inRange z = lMin <= z && z <= lMax + || cMin <= z && z <= cMax + || rMin <= z && z <= rMax + -- index into the generated vector @xs@ when generating elements + -- via @go@ + go (k,l) + | k >= szLen = Nothing + | l >= szLen = go (k+1,k+1 + if d == OnDiag then 0 else 1) + | otherwise = Just (((k,l),(imp IM.! k, imp IM.! l)), (k,l+1)) + -- Initialize the enumeration at the correct pair @(i,j)@. From + -- then on we can @take@ the correct number of elements, or stream + -- all of them. + initEnum All OnDiag = (0,0) + initEnum All NoDiag = (0,1) + initEnum (FromN s k) OnDiag + | s >= allSize = (szLen,szLen) + | otherwise = fromLinear szLn' s + initEnum (FromN s k) NoDiag + | s >= allSize = (szLen,szLen) + | otherwise = id *** (+1) $ fromLinear szLn' s + diff --git a/Data/Paired/Vector.hs b/Data/Paired/Vector.hs new file mode 100644 --- /dev/null +++ b/Data/Paired/Vector.hs @@ -0,0 +1,46 @@ + +module Data.Paired.Vector + ( module Data.Paired.Vector + , module Data.Paired.Common + ) where + +import Data.Vector.Generic as VG + +import Data.Paired.Common + + + +-- | Upper triangular elements. + +upperTriVG + :: (Vector v a, Vector w (a,a)) + => OnDiag + -> v a + -> (Int, w (a,a)) +upperTriVG d as = (z, unfoldrN z go (0,if d == OnDiag then 0 else 1)) + where la = VG.length as + z = la * (la + if d == OnDiag then 1 else 0) `div` 2 + go (k,l) + | k >= la = Nothing + | l >= la = go (k+1,k+1 + if d == OnDiag then 0 else 1) + | otherwise = Just ((as `VG.unsafeIndex` k, as `VG.unsafeIndex` l), (k,l+1)) +{-# Inline upperTriVG #-} + +-- | Outer pairing of all @as@ with all @bs@. This one is quasi-trivial, +-- but here for completeness. + +rectangularVG + :: (Vector va a, Vector vb b, Vector w (a,b)) + => va a + -> vb b + -> (Int, w (a,b)) +rectangularVG as bs = (z, unfoldrN z go (0,0)) + where la = VG.length as + lb = VG.length bs + z = la * lb + go (k,l) + | k >= la = Nothing + | l >= lb = go (k+1,0) + | otherwise = Just ((as `VG.unsafeIndex` k, bs `VG.unsafeIndex` l), (k,l+1)) +{-# Inline rectangularVG #-} + diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright Christian Hoener zu Siederdissen 2016 + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Christian Hoener zu Siederdissen nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/Math/TriangularNumbers.hs b/Math/TriangularNumbers.hs new file mode 100644 --- /dev/null +++ b/Math/TriangularNumbers.hs @@ -0,0 +1,78 @@ + +-- | Triangular numbers and various helper functions. +-- +-- Main use is for linearization of triangular array indexing. +-- +-- Triangular numbers: +-- @ +-- T_n = Σ_{k=1)^n k = 1 + 2 + 3 + ... + n = +-- +-- n * (n+1) / 2 = (n+1) `choose` 2 +-- @ +-- +-- + +module Math.TriangularNumbers where + + + +-- | Triangular numbers. +-- +-- https://oeis.org/A000217 + +triangularNumber :: Int -> Int +triangularNumber x = (x * (x+1)) `quot` 2 +{-# INLINE triangularNumber #-} + +-- | Size of an upper triangle starting at 'i' and ending at 'j'. "(0,N)" what +-- be the normal thing to use. + +linearizeUppertri :: (Int,Int) -> Int +linearizeUppertri (i,j) = triangularNumber $ j-i+1 +{-# INLINE linearizeUppertri #-} + +-- | Subword indexing. Given the longest subword and the current subword, +-- calculate a linear index "[0,..]". "(l,n)" in this case means "l"ower bound, +-- length "n". And "(i,j)" is the normal index. +-- +-- @ +-- 0 1 2 3 <- j = ... +-- +-- 0 1 2 3 i=0 +-- _ 4 5 6 i=1 +-- _ _ 7 8 i=2 +-- 9 i=3 +-- +-- i=2, j=3 -> (4+1) * i - tri i + j +-- +-- _ +-- _ _ the triangular number to subtract. +-- @ + +toLinear :: Int -> (Int,Int) -> Int +toLinear n (i,j) = adr n (i,j) + where + adr n (i,j) = (n+1)*i - triangularNumber i + j + {-# Inline adr #-} +{-# INLINE toLinear #-} + + + +-- | Linear index to paired. +-- +-- We have indices in @[0,N]@, and linear index @k@. +-- +-- @ +-- (N+1)*i - (i*(i+1)/2) + j == K +-- @ + +fromLinear :: Int -> Int -> (Int,Int) +fromLinear n' k' = (i,j) + where ll = (2*n+1) / 2 + rr = sqrt $ ((2*(n+1)+1) / 2)^2 - 2*k + n = fromIntegral n' + k = fromIntegral k' + i = floor $ ll - rr + 1 + j = k' - toLinear n' (i,0) +{-# Inline fromLinear #-} + diff --git a/Pipes/Parallel.hs b/Pipes/Parallel.hs new file mode 100644 --- /dev/null +++ b/Pipes/Parallel.hs @@ -0,0 +1,60 @@ + +-- | Pipes that introduce parallelism on different levels. + +module Pipes.Parallel where + +import Control.Monad.Codensity (lowerCodensity) +import Control.Monad (replicateM) +import Control.Parallel.Strategies (Strategy, parMap) +import Pipes + + + +-- | Evaluates chunks of pipes elements in parallel with a pure function. + +pipePar + :: (Monad m) + => Int + -- ^ number of elements to evaluate in parallel + -> Strategy b + -- ^ with which strategy + -> (a -> b) + -- ^ function to be mapped in parallel + -> Pipe a b m () +pipePar n strat f = pipeParBA n strat f (\as -> return ((),as)) (\() bs -> return bs) +{- + where + go = do + xs <- lowerCodensity . replicateM n $ lift await + let ys = parMap strat f xs + lowerCodensity $ mapM_ (lift . yield) ys + go +-} + +-- | Evaluates chunks of pipes elements in parallel with a pure function. +-- Before and after each parallel step, a monadic function is run. This +-- allows generation of certain statistics or information during runs. + +pipeParBA + :: (Monad m) + => Int + -- ^ number of elements to evaluate in parallel + -> Strategy b + -- ^ with which strategy + -> (a -> b) + -- ^ pure function to run in parallel + -> ([a] -> m (x,[a])) + -- ^ function to run before + -> (x -> [b] -> m [b]) + -- ^ function to run after + -> Pipe a b m () +pipeParBA n strat f bef aft = go + where + go = do + as' <- lowerCodensity . replicateM n $ lift await + (x,as) <- lift $ bef as' + let bs' = parMap strat f as + bs <- lift $ aft x bs' + lowerCodensity $ mapM_ (lift . yield) bs + go + diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,26 @@ +[![Build Status](https://travis-ci.org/choener/DPutils.svg?branch=master)](https://travis-ci.org/choener/DPutils) + +# DPutils + +Small set of utility functions. Currently centered around dynamic programming. + +*Math.TriangularNumbers* provides indexing into upper triangular tables. With +back and forth between the index pair and the linear index. + +*Data.Paired.Vector* provided rectangular and upper-triangular pairing of +elements from a vector. + +*Data.Paired.Foldable* is a more powerful generalization of such pairing for +any foldable container. We try to only retain elements that will be needed for +the pairing, while others are being filtered out. + +*Pipes.Parallel* provides some simple tools for parallelisation of tasks with +the pipes eco-system. + +#### Contact + +Christian Hoener zu Siederdissen +Leipzig University, Leipzig, Germany +choener@bioinf.uni-leipzig.de +http://www.bioinf.uni-leipzig.de/~choener/ + diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/changelog.md b/changelog.md new file mode 100644 --- /dev/null +++ b/changelog.md @@ -0,0 +1,14 @@ +0.0.0.2 +------- + +- Math.TriangularNumbers collects functions for triangular numbers and indexing +- Data.Paired.Foldable for upper triangular pairings [work in progress] +- renamed modules into new Paired hierarchy +- helper functions for parallelisation within a pipes pipe + +0.0.0.1 +------- + +- initial checkin +- upper triangular and rectangular combinations + diff --git a/tests/benchmark.hs b/tests/benchmark.hs new file mode 100644 --- /dev/null +++ b/tests/benchmark.hs @@ -0,0 +1,30 @@ + +module Main where + +import Criterion.Main +import Data.Vector as V + +import Math.TriangularNumbers + + + +-- | Run + +benchToLinear :: Int -> Int -> Int +benchToLinear i j = V.sum $ V.map (\z -> toLinear z (i,j)) $ V.enumFromN 10 1000 +{-# NoInline benchToLinear #-} + +benchFromLinear :: Int -> Int +benchFromLinear k = V.sum $ V.map (\z -> let (i,j) = fromLinear z k in i+j) $ V.enumFromN 10 1000 +{-# NoInline benchFromLinear #-} + +-- + +main :: IO () +main = defaultMain + [ bench "toLinear" $ nf (toLinear 10) (4,7) + , bench "fromLinear" $ nf (fromLinear 10) 41 + , bench "benchToLinear" $ nf (benchToLinear 4) 7 + , bench "benchFromLinear" $ nf benchFromLinear 41 + ] + diff --git a/tests/properties.hs b/tests/properties.hs new file mode 100644 --- /dev/null +++ b/tests/properties.hs @@ -0,0 +1,144 @@ + +module Main where + +import Data.List as L +import Data.Map.Strict as M +import Data.Tuple (swap) +import Data.Vector as V +import Debug.Trace +import Test.QuickCheck +import Test.Tasty +import Test.Tasty.QuickCheck as QC +import Test.Tasty.TH + +import Data.Paired.Vector as DPV +import Data.Paired.Foldable as DPF +import Math.TriangularNumbers + + + +-- * Data.Paired.Vector + +-- | + +prop_vector_upperTri_On :: NonNegative Int -> Bool +prop_vector_upperTri_On (NonNegative k) = V.toList vs == ls + where vs = snd $ upperTriVG OnDiag v + ls = [ (a,b) + | as@(a:_) <- L.init . L.tails $ V.toList v + , b <- as + ] + v = V.enumFromTo 0 k + +-- | + +prop_vector_upperTri_No :: NonNegative Int -> Bool +prop_vector_upperTri_No (NonNegative k) = V.toList vs == ls + where vs = snd $ upperTriVG NoDiag v + ls = [ (a,b) + | (a:as) <- L.init . L.tails $ V.toList v + , b <- as + ] + v = V.enumFromTo 0 k + +-- | + +prop_vector_rectangular :: NonNegative Int -> NonNegative Int -> Bool +prop_vector_rectangular (NonNegative k) (NonNegative l) = V.toList vs == ls + where vs = snd $ rectangularVG as bs + ls = [ (a,b) + | a <- V.toList as + , b <- V.toList bs + ] + as = V.enumFromTo 0 k + bs = V.enumFromTo 0 l + + + +-- * Data.Paired.Foldable + +-- | Generalized upper triangular elements. We want to enumerate all +-- elements, including those on the main diagonal. + +prop_foldable_upperTri_On_All :: (NonNegative Int, Bool) -> Bool +prop_foldable_upperTri_On_All (NonNegative n, b) + | chk = True + | otherwise = traceShow (ls,vs) False + where Right (_,_,vs) = DPF.upperTri (if b then KnownSize n else UnknownSize) OnDiag All xs + ls = [ ((a,b),(a,b)) + | as@(a:_) <- L.init . L.tails $ xs + , b <- as + ] + xs = [ 0 .. n-1 ] + chk = vs == ls + +-- | Only a subset of elements, starting at @k@ (counting from 0) and +-- taking @s@ elements. + +prop_foldable_upperTri_On_FromN :: (NonNegative Int, NonNegative Int, NonNegative Int, Bool) -> Bool +prop_foldable_upperTri_On_FromN (NonNegative n, NonNegative k, NonNegative s, b) + | chk = True + | otherwise = traceShow (ls,vs) False + where Right (_,_,vs) = DPF.upperTri (if b then KnownSize n else UnknownSize) OnDiag (FromN k s) xs + ls = L.take s + . L.drop k + $ [ ((a,b),(a,b)) + | as@(a:_) <- L.init . L.tails $ xs + , b <- as + ] + xs = [ 0 .. n-1 ] + chk = vs == ls + +prop_foldable_upperTri_No_All :: (NonNegative Int, Bool) -> Bool +prop_foldable_upperTri_No_All (NonNegative n, b) + | chk = True + | otherwise = traceShow (ls,vs) False + where Right (_,_,vs) = DPF.upperTri (if b then KnownSize n else UnknownSize) NoDiag All xs + ls = [ ((a,b),(a,b)) + | (a:as) <- L.init . L.tails $ xs + , b <- as + ] + xs = [ 0 .. n-1 ] + chk = vs == ls + +prop_foldable_upperTri_No_FromN :: (NonNegative Int, NonNegative Int, NonNegative Int, Bool) -> Bool +prop_foldable_upperTri_No_FromN (NonNegative n, NonNegative k, NonNegative s, b) + | chk = True + | otherwise = traceShow (ls,vs) False + where Right (_,_,vs) = DPF.upperTri (if b then KnownSize n else UnknownSize) NoDiag (FromN k s) xs + ls = L.take s + . L.drop k + $ [ ((a,b),(a,b)) + | (a:as) <- L.init . L.tails $ xs + , b <- as + ] + xs = [ 0 .. n-1 ] + chk = vs == ls + + + +-- * Math.TriangularNumbers + +-- | Test that each index pair @(i,j)@ is assigned a unique linear index +-- @k@ given @0 <= i <= j <= n@. + +prop_uniqueLinear :: NonNegative Int -> Bool +prop_uniqueLinear (NonNegative n) = M.null $ M.filter ((/=1) . L.length) mp + where mp = M.fromListWith (L.++) [ (toLinear n (i,j), [(i,j)]) | i <- [0..n], j <- [i..n] ] + +-- | Back and forth translation between paired and linear indices is +-- unique. + +prop_BackForth :: NonNegative Int -> Bool +prop_BackForth (NonNegative n) = L.and xs + where mb = M.fromList ls + mf = M.fromList $ L.map swap ls + ls = [ (toLinear n (i,j), (i,j)) | i <- [0..n], j <- [i..n] ] + xs = [ (mb M.! k == (i,j)) && (mf M.! (i,j) == k) && fromLinear n k == (i,j) + | (k,(i,j)) <- ls ] + +-- + +main :: IO () +main = $(defaultMainGenerator) +