diff --git a/Hungarian-Munkres.cabal b/Hungarian-Munkres.cabal --- a/Hungarian-Munkres.cabal +++ b/Hungarian-Munkres.cabal @@ -2,21 +2,23 @@ -- documentation, see http://haskell.org/cabal/users-guide/ name: Hungarian-Munkres -version: 0.1.3 +version: 0.1.4 synopsis: A Linear Sum Assignment Problem (LSAP) solver description: This library provide a Haskell binding to the libhungarian, a solver for Linear Sum Assignment Problem (LSAP) implemented in C language. It uses Hungarian algorithm , and runs - in O(n^3) time. + in O(n^3) time. This implementation is efficient. Benchmarks + versus pure haskell implementation are included (run + "cabal bench"). license: GPL-3 license-file: LICENSE author: Kai Zhang maintainer: Kai Zhang copyright: (c) 2014 Kai Zhang -category: Algorithm +category: Algorithms build-type: Simple -extra-source-files: cbits/hungarian.h +extra-source-files: include/hungarian.h cabal-version: >=1.10 library @@ -26,8 +28,22 @@ build-depends: base >=4.7 && <4.8 hs-source-dirs: src c-sources: cbits/hungarian.c + include-dirs: include + includes: hungarian.h + default-language: Haskell2010 +benchmark bench + type: exitcode-stdio-1.0 + hs-source-dirs: benchmarks + main-is: bench.hs default-language: Haskell2010 + build-depends: + base >=4.7 && <4.8 + , array + , Hungarian-Munkres + , Munkres + , random + , criterion test-suite tests type: exitcode-stdio-1.0 diff --git a/benchmarks/bench.hs b/benchmarks/bench.hs new file mode 100644 --- /dev/null +++ b/benchmarks/bench.hs @@ -0,0 +1,27 @@ +import Algorithms.Hungarian +import Criterion.Main +import Data.Algorithm.Munkres +import Data.Array.Unboxed +import System.Random + +sample :: [Double] +sample = take 4000 $ randomRs (-100, 100) (mkStdGen 4) + +sample' :: [Double] +sample' = take 40000 $ randomRs (-1000, 1000) (mkStdGen 24) + +sampleArray :: UArray (Int, Int) Double +sampleArray = listArray ((1,1), (40,50)) sample + +sampleArray' :: UArray (Int, Int) Double +sampleArray' = listArray ((1,1), (200,200)) sample' + +main :: IO () +main = defaultMain + [ bgroup "Hungarian Algorithm Benchmarks" + [ bench "C version (40 * 50)" $ nf (\x -> hungarian x 40 50) sample + , bench "Pure Haskell (40 * 50)" $ nf hungarianMethodDouble sampleArray + , bench "C version (200 * 200)" $ nf (\x -> hungarian x 200 200) sample' + , bench "Pure Haskell (200 * 200)" $ nf hungarianMethodDouble sampleArray' + ] + ] diff --git a/cbits/hungarian.c b/cbits/hungarian.c --- a/cbits/hungarian.c +++ b/cbits/hungarian.c @@ -33,23 +33,32 @@ #define hungarian_test_alloc(X) do {if ((void *)(X) == NULL) fprintf(stderr, "Out of memory in %s, (%s, line %d).\n", __FUNCTION__, __FILE__, __LINE__); } while (0) -double hungarian(double* data, int* result, int rows, int cols) { - int i,j; +double hungarian(double* data, int rows, int cols, size_t *from, size_t *to) { + size_t i, j, k, n; hungarian_problem_t p; int matrix_size = hungarian_init(&p, data , rows, cols, HUNGARIAN_MODE_MINIMIZE_COST) ; double cost = hungarian_solve(&p); - for(i=0;i + +typedef struct { + int num_rows; + int num_cols; + double** cost; + int** assignment; +} hungarian_problem_t; + +double hungarian(double* data, int rows, int cols, size_t *from, size_t *to); + +/** This method initialize the hungarian_problem structure and init + * the cost matrices (missing lines or columns are filled with 0). + * It returns the size of the quadratic(!) assignment matrix. **/ +int hungarian_init(hungarian_problem_t* p, + double* cost_matrix, + int rows, + int cols, + int mode); + +/** Free the memory allocated by init. **/ +void hungarian_free(hungarian_problem_t* p); + +/** This method computes the optimal assignment. **/ +double hungarian_solve(hungarian_problem_t* p); + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/src/Algorithms/Hungarian.hs b/src/Algorithms/Hungarian.hs --- a/src/Algorithms/Hungarian.hs +++ b/src/Algorithms/Hungarian.hs @@ -1,4 +1,5 @@ {-# LANGUAGE ForeignFunctionInterface #-} +{-# LANGUAGE EmptyDataDecls #-} module Algorithms.Hungarian ( hungarian @@ -6,33 +7,35 @@ ) where import Data.List -import Foreign.Ptr +import Foreign import Foreign.C -import Foreign.Marshal.Array import System.IO.Unsafe foreign import ccall "hungarian" - c_hungarian :: Ptr CDouble -> Ptr CInt -> CInt -> CInt -> IO Double + c_hungarian :: Ptr CDouble -> CInt -> CInt -> Ptr CSize -> Ptr CSize -> IO Double -- | solve the LSAP by hungarian algorithm, return assignment and score hungarian :: [Double] -- ^ row majored flat matrix -> Int -- ^ number of rows -> Int -- ^ number of columns -> ([(Int, Int)], Double) -hungarian costMatrix nrows ncols = unsafePerformIO $ do +hungarian costMatrix rows cols = unsafePerformIO $ do withArray (map realToFrac costMatrix) $ \input -> - allocaArray (nrows * ncols) $ \output -> do - cost <- c_hungarian input output (fromIntegral nrows) (fromIntegral ncols) - results <- peekArray (nrows * ncols) output - return (getAssign results, realToFrac cost) + allocaArray n $ \from -> allocaArray n $ \to -> do + cost <- c_hungarian input (fromIntegral rows) (fromIntegral cols) + from to + froms <- peekArray n from + tos <- peekArray n to + return (zipWith f froms tos, realToFrac cost) where - getAssign :: [CInt] -> [(Int, Int)] - getAssign = snd . foldl' step (0,[]) - step (i,assign) x | x == 0 = (i+1, assign) - | otherwise = (i+1, (i `div` ncols, i `mod` ncols) : assign) + f x y = (fromIntegral x, fromIntegral y) + n = min rows cols {-# INLINE hungarian #-} -- | solve the LSAP by hungarian algorithm, return score only hungarianScore :: [Double] -> Int -> Int -> Double -hungarianScore costMatrix nrows = snd . hungarian costMatrix nrows +hungarianScore costMatrix rows cols = unsafePerformIO $ do + withArray (map realToFrac costMatrix) $ \input -> do + fmap realToFrac $ c_hungarian input (fromIntegral rows) + (fromIntegral cols) nullPtr nullPtr {-# INLINE hungarianScore #-}