packages feed

permutation (empty) → 0.1

raw patch · 6 files changed

+557/−0 lines, 6 filesdep +basedep +containersbuild-type:Customsetup-changed

Dependencies added: base, containers

Files

+ Data/Permutation.hs view
@@ -0,0 +1,265 @@+{-# LANGUAGE CPP #-}+{-# OPTIONS_GHC -fglasgow-exts #-}+-----------------------------------------------------------------------------+-- |+-- Module     : Data.Permutation+-- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>+-- License    : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability  : experimental+--++module Data.Permutation (+    -- * The Permutation type+    Permutation,++    -- * Creating permutations+    permutation,+    identity,+    inverse,+    +    -- * Permutation properties+    size,+    apply,++    -- * Applying permutations+    applyWith,+    invertWith,++    -- * Converstion to/from other types+    -- ** @ForeignPtr@s+    fromForeignPtr,+    toForeignPtr,+    +    -- ** Lists+    toList,+    fromList,++    -- * Unsafe operations+    withPermutationPtr,+    unsafePermutation,+    unsafeApply,+    +    ) where+        +import Control.Monad         ( foldM, liftM )+import Data.IntSet           ( IntSet )+import qualified Data.IntSet as IntSet+import Foreign               ( Ptr, ForeignPtr, mallocForeignPtrArray, +                               withForeignPtr, pokeArray, peekArray, +                               advancePtr, peek, peekElemOff, pokeElemOff ) +import System.IO.Unsafe      ( unsafePerformIO )++#if defined(__GLASGOW_HASKELL__)+import GHC.Base                 ( realWorld# )+import GHC.IOBase               ( IO(IO) )+#endif++inlinePerformIO :: IO a -> a+#if defined(__GLASGOW_HASKELL__)+inlinePerformIO (IO m) = case m realWorld# of (# _, r #) -> r+#else+inlinePerformIO = unsafePerformIO+#endif+{-# INLINE inlinePerformIO #-}+++-- | Represents a permutation of the integers @[0..n)@.        +data Permutation =+        Perm {-# UNPACK #-} !Int             +             {-# UNPACK #-} !(ForeignPtr Int)+             {-# UNPACK #-} !Int++-- | Get the size, raw data array and offset+toForeignPtr :: Permutation -> (Int, ForeignPtr Int, Int)+toForeignPtr (Perm n f o) = (n, f, o)++-- | Convert size and raw array to a permutation.  No validation is+-- performed on the arguments.+fromForeignPtr :: Int -> ForeignPtr Int -> Int -> Permutation+fromForeignPtr = Perm++-- | Get the size of the permutation.+size :: Permutation -> Int+size (Perm n _ _) = n+{-# INLINE size #-}++-- | Perform an operation, given a pointer to the start of the+-- permutation data+withPermutationPtr :: Permutation -> (Ptr Int -> IO a) -> IO a+withPermutationPtr (Perm _ fptr off) f =+    withForeignPtr fptr $ \ptr ->+        f (ptr `advancePtr` off)++-- | Apply a permutation to an integer.  The integer must be in the range+--   @[0..n)@.+apply :: Permutation -> Int -> Int+apply p@(Perm n _ _) i +    | i < 0 || i >= n =+        error $ +            "applyPerm: Tried to apply permutation of size `" ++ show n +++            "' to the value `" ++ show i ++ "'."+    | otherwise =+        unsafeApply p i+{-# INLINE apply #-}++-- | Same as 'apply' but does not range-check the argument.+unsafeApply :: Permutation -> Int -> Int+unsafeApply p i =+    inlinePerformIO $ do+        withPermutationPtr p $ flip peekElemOff i+{-# INLINE unsafeApply #-}++-- | Create a permutation from a list of values.  The list must be of length+--   @n@ and contain each integer in @{0, 1, ..., (n-1) }@ exactly once.+--   The permutation that is returned will send the integer @i@ to its index+--   in the list.+permutation :: Int -> [Int] -> Permutation+permutation n is = +    let p = unsafePermutation n is+    in case isValid p of+           False -> error $ "Not a valid permutation."+           True  -> p++-- | Same as 'permutation', but does not check that the inputs are valid.+unsafePermutation :: Int -> [Int] -> Permutation+unsafePermutation n is = +    unsafePerformIO $ do+        fptr <- mallocForeignPtrArray n+        withForeignPtr fptr $ \ptr -> pokeArray ptr is+        return $ fromForeignPtr n fptr 0+{-# NOINLINE unsafePermutation #-}+++fromList :: [Int] -> Permutation+fromList is = permutation (length is) is++toList :: Permutation -> [Int]+toList p = unsafePerformIO $ +               withPermutationPtr p $ peekArray (size p)++-- | Create an identity permutation of the given size.+identity :: Int -> Permutation+identity n =+    unsafePerformIO $ do+        fptr <- mallocForeignPtrArray n+        withForeignPtr fptr $ \ptr -> pokeArray ptr [0..(n-1)]+        return $ fromForeignPtr n fptr 0+++-- | Get the inverse of a permutation.+inverse :: Permutation -> Permutation+inverse p =+    let n = size p+    in+        unsafePerformIO $ do+            fptr <- mallocForeignPtrArray n+            withForeignPtr fptr $ \ptr -> do+                pokeArray ptr [0..(n-1)]+                invertWith (swap ptr) p+            return $ fromForeignPtr n fptr 0+    where+        swap :: Ptr Int -> Int -> Int -> IO ()+        swap ptr i j = do+            x <- peekElemOff ptr i+            y <- peekElemOff ptr j+            pokeElemOff ptr i y+            pokeElemOff ptr j x+{-# NOINLINE inverse #-}+++++isValid :: Permutation -> Bool+isValid p@(Perm n _ _) =+    unsafePerformIO $+        withPermutationPtr p $ \ptr -> do+            liftM and $+                mapM (\i -> peekElemOff ptr i +                            >>= \x -> isValidI ptr x i)+                     [0..(n-1)]+             +    where+        isValidI :: Ptr Int -> Int -> Int -> IO Bool+        isValidI ptr x i =+            liftM and $ +                sequence [ inRange x, isUnique x ptr i ]+    +        inRange :: Int -> IO Bool+        inRange x =+            return $ x >= 0 && x < n+        +        isUnique :: Int -> Ptr Int -> Int -> IO Bool+        isUnique x ptr' n'+            | n' == 0 =+                return True+            | otherwise = do+                x' <- peek ptr'+                if x' == x +                    then return False+                    else isUnique x (ptr' `advancePtr` 1) (n'-1)+    +-- | @applyWith swap perm@ applies the permutation as a sequence of swaps.  After +-- this function is applied, @OUT[i] = IN[P[i]]@+applyWith :: (Monad m) => (Int -> Int -> m ()) -> Permutation -> m ()+applyWith swap p =+    let n = size p+    in foldM (flip $ doCycle swap) IntSet.empty [0..(n-1)] >> return ()+    +    where+        doCycle :: (Monad m) => +            (Int -> Int -> m ()) -> Int -> IntSet -> m (IntSet)+        doCycle swp i visited = +            if i `IntSet.member` visited +                then return visited+                else let visited' = IntSet.insert i visited+                         next     = unsafeApply p i+                     in doCycle' swp i i next visited'+            +        doCycle' :: (Monad m) => +            (Int -> Int -> m ()) -> Int -> Int -> Int -> IntSet -> m (IntSet)+        doCycle' swp start cur next visited+            | next == start =+                return visited +            | otherwise = +                let visited' = IntSet.insert next visited+                    next'    = unsafeApply p next+                in do+                    swp cur next+                    doCycle' swp start next next' visited'++-- | @invertWith swap p@ applies the inverse of the permutation as a +-- sequence of swaps.  After this function is applied, @OUT[P[i]] = IN[i]@+invertWith :: (Monad m) => (Int -> Int -> m ()) -> Permutation -> m ()+invertWith swap p =+    let n = size p+    in foldM (flip $ doCycle swap) IntSet.empty [0..(n-1)] >> return ()+    +    where+        doCycle :: Monad m => +            (Int -> Int -> m ()) -> Int -> IntSet -> m (IntSet)+        doCycle swp i visited = +            if i `IntSet.member` visited +                then return visited+                else let visited' = IntSet.insert i visited+                         cur      = unsafeApply p i+                     in doCycle' swp i cur visited'+            +        doCycle' :: Monad m => (Int -> Int -> m ()) -> Int -> Int -> IntSet -> m (IntSet)+        doCycle' swp start cur visited+            | cur == start =+                return visited +            | otherwise = +                let visited' = IntSet.insert cur visited+                    cur'     = unsafeApply p cur+                in do+                    swp start cur+                    doCycle' swp start cur' visited'+++instance Show Permutation where+    show p = "permutation "     ++ show (size p) ++ " " ++ show (toList p)+    +instance Eq Permutation where+    (==) p q = (size p == size q) && (toList p == toList q)+    
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) Patrick Perry <patperry@stanford.edu> 2008++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:++1. Redistributions of source code must retain the above copyright+   notice, this list of conditions and the following disclaimer.++2. 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.++3. Neither the name of the author nor the names of his contributors+   may be used to endorse or promote products derived from this software+   without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``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 AUTHORS 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.
+ Setup.lhs view
@@ -0,0 +1,14 @@+#!/usr/bin/env runhaskell+> import Distribution.Simple+> import System.Cmd+> import System.Exit ( ExitCode(..) )+>+> testing _ _ _ _ = do+>     err <- system "make -C tests"+>     system "make -s -C tests clean"+>     if err /= ExitSuccess+>         then ioError $ userError $ "failed"+>         else return ()+>+> main = defaultMainWithHooks defaultUserHooks+>        {runTests=testing}
+ permutation.cabal view
@@ -0,0 +1,26 @@+name:            permutation+version:         0.1+homepage:        http://stat.stanford.edu/~patperry/code/permutation+synopsis:        Permutations and combinations+description:+    A library for representing and applying permutations.+    .+category:        Math+license:         BSD3+license-file:    LICENSE+copyright:       (c) 2008. Patrick Perry <patperry@stanford.edu>+author:          Patrick Perry+maintainer:      Patrick Perry <patperry@stanford.edu>+cabal-version: >= 1.2.0+build-type:      Custom+tested-with:     GHC ==6.8.2++extra-source-files:  tests/Unit.hs+                     tests/Makefile++library+    exposed-modules: Data.Permutation++    build-depends:   base, containers++    ghc-options:     -Wall
+ tests/Makefile view
@@ -0,0 +1,8 @@+all:+	ghc -O -i. -i.. Unit.hs --make+	./Unit++clean:+	rm -f *~ Unit *.hi *.o \+		../Data/Permutation.hi ../Data/Permutation.o \+		../System/Random/Permutation.o
+ tests/Unit.hs view
@@ -0,0 +1,214 @@++import Control.Monad ( forM_ )+import Data.List ( sortBy )+import Data.Ord  ( comparing )+import Foreign+import System.Random+import Test.HUnit++import Data.Permutation++++swapElems :: Storable a => Ptr a -> Int -> Int -> IO ()+swapElems ptr i j =  do+    x <- peekElemOff ptr i+    y <- peekElemOff ptr j+    pokeElemOff ptr i y+    pokeElemOff ptr j x++newRandomArray :: (Random a, Storable a) => Int -> IO (Ptr a)+newRandomArray n = do+    xs <- sequence $ replicate n randomIO+    newArray xs++newArrayCopy :: Storable a => Int -> Ptr a -> IO (Ptr a)+newArrayCopy n x = do+    y <- mallocArray n+    copyArray y x n+    return y++newRandomPerm :: Int -> IO (Permutation)+newRandomPerm n = do+    xs <- sequence $ replicate n randomIO :: IO [Int]+    let ixs  = zip [0..] xs+        ixs' = sortBy (comparing snd) ixs+        is   = (fst . unzip) ixs'+        p    = permutation n is+    return p+    ++testIdentity :: Test+testIdentity = TestCase $ do+    let n = 30+        p = identity n+    +    x <- newRandomArray n :: IO (Ptr Double)+    y <- newArrayCopy n x+    applyWith (swapElems y) p+    +    xs <- peekArray n x+    ys <- peekArray n y+    +    assertEqual ""+        xs+        ys+        +        +testSmallCycle :: Test+testSmallCycle = TestCase $ do+    let n = 3+        p = permutation 3  [ 2, 0, 1 ]++    x <- newArray [ 6, 7, 8 ] :: IO (Ptr Int)+    applyWith (swapElems x) p+    +    xs <- peekArray n x+    assertEqual ""+        [ 8, 6, 7 ]        +        xs++testSmallCycleInv :: Test+testSmallCycleInv = TestCase $ do+    let n = 3+        p = permutation 3  [ 2, 0, 1 ]++    x <- newArray [ 6, 7, 8 ] :: IO (Ptr Int)+    invertWith (swapElems x) p+    +    xs <- peekArray n x+    assertEqual ""+        [ 7, 8, 6 ]+        xs+++testSingleCycleInv :: Test+testSingleCycleInv = TestCase $ do+    let n = 5+        p = permutation n [ 3, 0, 1, 4, 2 ]+        +    x <- newArray [ 0.237,  0.382,  0.818, -0.413,  0.037 ] :: IO (Ptr Double)+    invertWith (swapElems x) p+    +    xs <- peekArray n x+    assertEqual ""+        [ 0.382, 0.818, 0.037, 0.237, -0.413 ]+        xs++testSingleCycle :: Test+testSingleCycle = TestCase $ do+    let n = 5+        p = permutation n [ 3, 0, 1, 4, 2 ]++    x <- newArray [ 0.237,  0.382,  0.818, -0.413,  0.037 ] :: IO (Ptr Double)+    applyWith (swapElems x) p+    +    xs <- peekArray n x+    assertEqual ""+        [ -0.413, 0.237, 0.382, 0.037, 0.818 ]+        xs+        +testMultiCycleInv :: Test+testMultiCycleInv = TestCase $ do+    let n = 5+        p = permutation n [ 2, 3, 0, 4, 1 ]+        +    x <- newArray [ 0.237,  0.382,  0.818, -0.413,  0.037 ] :: IO (Ptr Double)+    invertWith (swapElems x) p+    +    xs <- peekArray n x+    assertEqual ""+        [ 0.818, 0.037, 0.237, 0.382, -0.413 ]+        xs++testMultiCycle :: Test+testMultiCycle = TestCase $ do+    let n = 5+        p = permutation n [ 2, 3, 0, 4, 1 ]+    x <- newArray [ 0.237,  0.382,  0.818, -0.413,  0.037 ] :: IO (Ptr Double)+    applyWith (swapElems x) p+    +    xs <- peekArray n x+    assertEqual ""+        [ 0.818, -0.413, 0.237, 0.037, 0.382 ]+        xs++testApply :: Int -> Test+testApply n = TestCase $ do+    x <- newRandomArray n :: IO (Ptr Double)+    y <- newArrayCopy n x+    p <- newRandomPerm n+    applyWith  (swapElems y) p+    +    xs <- peekArray n x+    ys <- peekArray n y++    forM_ [0..(n-1)] $ \i ->+        assertEqual ("position " ++ show i) (xs !! (apply p i)) (ys !! i)++testInvert :: Int -> Test+testInvert n = TestCase $ do+    x <- newRandomArray n :: IO (Ptr Double)+    y <- newArrayCopy n x+    p <- newRandomPerm n+    invertWith  (swapElems y) p+    +    xs <- peekArray n x+    ys <- peekArray n y++    forM_ [0..(n-1)] $ \i ->+        assertEqual ("position " ++ show i) (xs !! i) (ys !! (apply p i))+++testApplyThenInv :: Int -> Test+testApplyThenInv n = TestCase $ do+    x <- newRandomArray n :: IO (Ptr Double)+    y <- newArrayCopy n x+    p <- newRandomPerm n+    applyWith  (swapElems y) p+    invertWith (swapElems y) p+    +    xs <- peekArray n x+    ys <- peekArray n y+    +    assertEqual "" xs ys++testInvThenApply :: Int -> Test+testInvThenApply n = TestCase $ do+    x <- newRandomArray n :: IO (Ptr Double)+    y <- newArrayCopy n x+    p <- newRandomPerm n+    invertWith (swapElems y) p+    applyWith  (swapElems y) p+    +    xs <- peekArray n x+    ys <- peekArray n y+    +    assertEqual "" xs ys++++ +smallTests = [ TestLabel "identity"             testIdentity+             , TestLabel "small cycle"          testSmallCycle+             , TestLabel "small cycle inverse"  testSmallCycleInv+             , TestLabel "single cycle"         testSingleCycle+             , TestLabel "single cycle inverse" testSingleCycleInv+             , TestLabel "multi cycle"          testMultiCycle+             , TestLabel "multi cycle inverse"  testMultiCycleInv+             ]++ +main =+    let ns      = [ 0, 1, 2, 3, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048 ]+        aTests  = map (\n -> TestLabel ("apply (" ++ show n ++ ")") +                                       (testApply n)) ns+        iTests  = map (\n -> TestLabel ("apply (" ++ show n ++ ")") +                                       (testInvert n)) ns+        aiTests = map (\n -> TestLabel ("apply then inverse (" ++ show n ++ ")") +                                       (testApplyThenInv n)) ns+        iaTests = map (\n -> TestLabel ("inverse then apply (" ++ show n ++ ")") +                                       (testInvThenApply n)) ns+        tests   = TestList $ smallTests ++ aTests ++ aiTests ++ iaTests+    in runTestTT tests+