diff --git a/.ghci b/.ghci new file mode 100644 --- /dev/null +++ b/.ghci @@ -0,0 +1,1 @@ +:set -isrc -idist/build/autogen -optP-include -optPdist/build/autogen/cabal_macros.h diff --git a/.gitignore b/.gitignore new file mode 100644 --- /dev/null +++ b/.gitignore @@ -0,0 +1,16 @@ +dist +docs +wiki +TAGS +tags +wip +stats +.DS_Store +.*.swp +.*.swo +*.o +*.hi +*~ +*# +.cabal-sandbox +cabal.sandbox.config diff --git a/.travis.yml b/.travis.yml new file mode 100644 --- /dev/null +++ b/.travis.yml @@ -0,0 +1,8 @@ +language: haskell +notifications: + irc: + channels: + - "irc.freenode.org#haskell-lens" + skip_join: true + template: + - "\x0313sparse\x03/\x0306%{branch}\x03 \x0314%{commit}\x03 %{build_url} %{message}" diff --git a/.vim.custom b/.vim.custom new file mode 100644 --- /dev/null +++ b/.vim.custom @@ -0,0 +1,31 @@ +" Add the following to your .vimrc to automatically load this on startup + +" if filereadable(".vim.custom") +" so .vim.custom +" endif + +function StripTrailingWhitespace() + let myline=line(".") + let mycolumn = col(".") + silent %s/ *$// + call cursor(myline, mycolumn) +endfunction + +" enable syntax highlighting +syntax on + +" search for the tags file anywhere between here and / +set tags=TAGS;/ + +" highlight tabs and trailing spaces +set listchars=tab:‗‗,trail:‗ +set list + +" f2 runs hasktags +map :exec ":!hasktags -x -c --ignore src" + +" strip trailing whitespace before saving +" au BufWritePre *.hs,*.markdown silent! cal StripTrailingWhitespace() + +" rebuild hasktags after saving +au BufWritePost *.hs silent! :exec ":!hasktags -x -c --ignore src" diff --git a/COPYING b/COPYING deleted file mode 100644 --- a/COPYING +++ /dev/null @@ -1,26 +0,0 @@ - -Copyright (c) 2013, Hans Höglund -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 the nor the - names of its 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 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/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,26 @@ +Copyright 2013 Edward Kmett + +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. + +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. diff --git a/Setup.lhs b/Setup.lhs --- a/Setup.lhs +++ b/Setup.lhs @@ -1,4 +1,55 @@ -#! /usr/bin/env runhaskell +#!/usr/bin/runhaskell +\begin{code} +{-# OPTIONS_GHC -Wall #-} +module Main (main) where -> import Distribution.Simple -> main = defaultMain +import Data.List ( nub ) +import Data.Version ( showVersion ) +import Distribution.Package ( PackageName(PackageName), Package, PackageId, InstalledPackageId, packageVersion, packageName ) +import Distribution.PackageDescription ( PackageDescription(), TestSuite(..) ) +import Distribution.Simple ( defaultMainWithHooks, UserHooks(..), simpleUserHooks ) +import Distribution.Simple.Utils ( rewriteFile, createDirectoryIfMissingVerbose, copyFiles ) +import Distribution.Simple.BuildPaths ( autogenModulesDir ) +import Distribution.Simple.Setup ( BuildFlags(buildVerbosity), Flag(..), fromFlag, HaddockFlags(haddockDistPref)) +import Distribution.Simple.LocalBuildInfo ( withLibLBI, withTestLBI, LocalBuildInfo(), ComponentLocalBuildInfo(componentPackageDeps) ) +import Distribution.Text ( display ) +import Distribution.Verbosity ( Verbosity, normal ) +import System.FilePath ( () ) + +main :: IO () +main = defaultMainWithHooks simpleUserHooks + { buildHook = \pkg lbi hooks flags -> do + generateBuildModule (fromFlag (buildVerbosity flags)) pkg lbi + buildHook simpleUserHooks pkg lbi hooks flags + , postHaddock = \args flags pkg lbi -> do + copyFiles normal (haddockOutputDir flags pkg) [] + postHaddock simpleUserHooks args flags pkg lbi + } + +haddockOutputDir :: Package p => HaddockFlags -> p -> FilePath +haddockOutputDir flags pkg = destDir where + baseDir = case haddockDistPref flags of + NoFlag -> "." + Flag x -> x + destDir = baseDir "doc" "html" display (packageName pkg) + +generateBuildModule :: Verbosity -> PackageDescription -> LocalBuildInfo -> IO () +generateBuildModule verbosity pkg lbi = do + let dir = autogenModulesDir lbi + createDirectoryIfMissingVerbose verbosity True dir + withLibLBI pkg lbi $ \_ libcfg -> do + withTestLBI pkg lbi $ \suite suitecfg -> do + rewriteFile (dir "Build_" ++ testName suite ++ ".hs") $ unlines + [ "module Build_" ++ testName suite ++ " where" + , "deps :: [String]" + , "deps = " ++ (show $ formatdeps (testDeps libcfg suitecfg)) + ] + where + formatdeps = map (formatone . snd) + formatone p = case packageName p of + PackageName n -> n ++ "-" ++ showVersion (packageVersion p) + +testDeps :: ComponentLocalBuildInfo -> ComponentLocalBuildInfo -> [(InstalledPackageId, PackageId)] +testDeps xs ys = nub $ componentPackageDeps xs ++ componentPackageDeps ys + +\end{code} diff --git a/benchmarks/mm.hs b/benchmarks/mm.hs new file mode 100644 --- /dev/null +++ b/benchmarks/mm.hs @@ -0,0 +1,49 @@ +import Control.Applicative +import Control.DeepSeq +import Criterion.Main +import Data.Array.Unboxed as A +import Data.Vector.Generic as G +import Data.Vector.Unboxed as U +import Sparse.Matrix as M +import Sparse.Matrix.Heap as Heap + +instance NFData (UArray i e) + +main :: IO () +main = defaultMain + [ bench "naive I_32" $ nf (\x -> mmul x x) $ array ((0,0),(31,31)) $ [ ((i, j), if i == j then 1 else 0) | i <- [0..31], j <- [0..31] ] + , bench "I_32 new" $ nf (\x -> x * x) (ident 32 :: Mat U.Vector Int) + , bench "I_64 new" $ nf (\x -> x * x) (ident 64 :: Mat U.Vector Int) + , bench "I_128 new" $ nf (\x -> x * x) (ident 128 :: Mat U.Vector Int) + -- , bench "I_256" $ nf (\x -> x * x) (ident 256 :: Mat U.Vector Int) + -- , bench "I_512" $ nf (\x -> x * x) (ident 1024 :: Mat U.Vector Int) + -- , bench "I_1024" $ nf (\x -> x * x) (ident 1024 :: Mat U.Vector Int) + , bench "naive 32x32" $ nf (\x -> mmul x x) $ listArray ((0,0),(31,31)) $ Prelude.replicate (32*32) 1 + , bench "32x32 Int" $ nf (\x -> x * x) blockInt + , bench "32x32 ()" $ nf (\x -> multiplyWith const (Heap.streamHeapWith const) x x) blockUnit + , bench "naive 128x128" $ nf (\x -> mmul x x) $ listArray ((0,0),(127,127)) $ Prelude.replicate (128*128) 1 + , bench "128x128 Int" $ nf (\x -> x * x) blockInt128 + , bench "128x128 ()" $ nf (\x -> multiplyWith const (Heap.streamHeapWith const) x x) blockUnit128 + ] + +blockInt :: Mat U.Vector Int +blockInt = M.fromList $ Prelude.zip (Key <$> [0..31] <*> [0..31]) (repeat 1) + +blockInt128 :: Mat U.Vector Int +blockInt128 = M.fromList $ Prelude.zip (Key <$> [0..127] <*> [0..127]) (repeat 1) + +blockUnit :: Mat U.Vector () +blockUnit = M.fromList $ Prelude.zip (Key <$> [0..31] <*> [0..31]) (repeat ()) + +blockUnit128 :: Mat U.Vector () +blockUnit128 = M.fromList $ Prelude.zip (Key <$> [0..127] <*> [0..127]) (repeat ()) + +mmul :: UArray (Int,Int) Int -> UArray (Int,Int) Int -> UArray (Int,Int) Int +mmul x y = accumArray (+) 0 ((i0,k0),(i1,k1)) $ do + i <- range (i0,i1) + j <- range (max j0 j0',min j1 j1') + k <- range (k0,k1) + return ((i,k),x A.!(i,j) * y A.!(j,k)) + where + ((i0,j0),(i1,j1)) = bounds x + ((j0',k0),(j1',k1)) = bounds y diff --git a/sparse.cabal b/sparse.cabal --- a/sparse.cabal +++ b/sparse.cabal @@ -1,32 +1,153 @@ +name: sparse +category: Data, Vector +version: 0.6 +license: BSD3 +cabal-version: >= 1.8 +license-file: LICENSE +author: Edward A. Kmett +maintainer: Edward A. Kmett +stability: experimental +homepage: http://github.com/ekmett/sparse +bug-reports: http://github.com/ekmett/sparse/issues +copyright: Copyright (C) 2013 Edward A. Kmett +build-type: Custom +synopsis: A playground of sparse linear algebra primitives using Morton ordering -name: sparse -version: 0.5 -cabal-version: >= 1.10 -author: Hans Hoglund -maintainer: Hans Hoglund -license: BSD3 -license-file: COPYING -synopsis: Lightweight parsing library based on partial functions. -category: -tested-with: GHC -build-type: Simple +extra-source-files: + .ghci + .travis.yml + .gitignore + .vim.custom description: - To be written. + A playground of sparse linear algebra primitives using Morton ordering + . + The design of this library is describe in the following articles on FP Complete's School of Haskell. + . + 1. + . + 2. + . + 3. + . + 4. + . + 5. source-repository head - type: git - location: git://github.com/hanshoglund/sparse.git + type: git + location: git://github.com/ekmett/sparse.git +-- Build the properties test if we're building tests +flag test-properties + default: True + manual: True + +-- You can disable the doctests test suite with -f-test-doctests +flag test-doctests + default: True + manual: True + +flag optimize + default: False + manual: True + +flag llvm + default: False + manual: True + library + build-depends: + base >= 4 && < 5, + contravariant >= 0.4.2 && < 1, + deepseq >= 1.1 && < 1.4, + hybrid-vectors >= 0.1 && < 1, + lens >= 3.9 && < 4, + primitive >= 0.5 && < 0.6, + transformers >= 0.3 && < 0.4, + vector >= 0.10 && < 0.11, + vector-algorithms >= 0.5 && < 0.6 + + hs-source-dirs: src + + exposed-modules: + Sparse.Matrix + Sparse.Matrix.Internal.Fusion + Sparse.Matrix.Internal.Heap + Sparse.Matrix.Internal.Key + + ghc-options: -Wall + + if flag(optimize) + ghc-options: -O2 + else + ghc-options: -O0 + + if flag(llvm) + ghc-options: -fllvm + + if impl(ghc<6.13) + Ghc-Options: -finline-if-enough-args -fno-method-sharing + +test-suite properties + type: exitcode-stdio-1.0 + main-is: properties.hs + hs-source-dirs: tests + ghc-options: -w -threaded -rtsopts -with-rtsopts=-N -fsimpl-tick-factor=400 + if !flag(test-properties) + buildable: False + else build-depends: - base >= 4 && < 5, - nats, - semigroups, - pointed, - monadplus - - hs-source-dirs: src - default-language: Haskell2010 - exposed-modules: - Data.Sparse + base, + containers, + hybrid-vectors, + lens, + linear >= 1.2 && < 2, + QuickCheck >= 2.5, + sparse, + test-framework >= 0.6, + test-framework-quickcheck2 >= 0.3, + test-framework-th >= 0.2, + transformers, + vector + +-- Verify the results of the examples +test-suite doctests + type: exitcode-stdio-1.0 + main-is: doctests.hs + ghc-options: -Wall -threaded + hs-source-dirs: tests + + if !flag(test-doctests) + buildable: False + else + build-depends: + base, + bytestring, + containers, + directory >= 1.0, + deepseq, + doctest >= 0.9.1, + filepath, + mtl, + semigroups >= 0.9, + simple-reflect >= 0.3.1 + + if impl(ghc<7.6.1) + ghc-options: -Werror + +-- matrix-matrix multiplication +benchmark mm + type: exitcode-stdio-1.0 + main-is: mm.hs + ghc-options: -Wall -O2 -threaded -fdicts-cheap -funbox-strict-fields -fsimpl-tick-factor=400000 + hs-source-dirs: benchmarks + build-depends: + array, + base, + criterion, + deepseq, + sparse, + vector + + diff --git a/src/Data/Sparse.hs b/src/Data/Sparse.hs deleted file mode 100644 --- a/src/Data/Sparse.hs +++ /dev/null @@ -1,215 +0,0 @@ - -{-# LANGUAGE GeneralizedNewtypeDeriving, - OverloadedStrings, - TypeOperators, - DeriveFunctor, - DeriveFoldable, - FlexibleInstances - #-} - - -------------------------------------------------------------------------------------- --- | --- Copyright : (c) Hans Hoglund 2012 --- --- License : BSD-style --- --- Maintainer : hans@hanshoglund.se --- Stability : experimental --- Portability : non-portable (GNTD, DeriveFunctor, OverloadedStrings) --- --- Lightweight parsing library based on partial functions. --- -------------------------------------------------------------------------------------- - -module Data.Sparse ( - -- * Sparse - SparseT, - Sparse, - asSparse, - - -- * Running - runSparseT, - runSparseT', - runSparse, - runSparse', - - -- * Primitives - headP, - splitP, - - -- * Basic parsers - char, - charIs, - string, - stringIs, - - -- * Combinators - optionally, - optionallyMaybe, - Data.Sparse.optional, - between, - skipMany1, - skipMany, - many1, - sepBy, - sepBy1, - sepEndBy1, - sepEndBy, - endBy1, - endBy, - count -) where - -import Data.String -import Data.Semigroup -import Data.Foldable(Foldable) -import Control.Applicative -import Control.Monad.Plus - --- TODO -instance Semigroup (Partial a b) where (<>) = mplus - - - - - - -newtype a ?-> b = PartialP { getPartialP :: a -> Maybe (a, b) } - -instance Functor ((?->) r) where - fmap f (PartialP g) = PartialP (fmap (fmap f) . g) - -instance Monad ((?->) r) where - return x = PartialP (\a -> Just (a, x)) - PartialP f >>= k = PartialP $ \r -> (f r >>= \(r1, x) -> getPartialP (k x) r1) - -instance MonadPlus ((?->) r) where - mzero = PartialP (const Nothing) - PartialP f `mplus` PartialP g = PartialP $ \x -> f x `mplus` g x - -instance Applicative ((?->) r) where - pure = return - (<*>) = ap - -instance Alternative ((?->) r) where - empty = mzero - (<|>) = mplus - -instance Semigroup ((?->) a b) where - (<>) = mplus - -instance Monoid ((?->) a b) where - mempty = mzero - mappend = mplus - - ----------- - -newtype SparseT a b = SparseT { getSparseT :: a ?-> b } - deriving (Semigroup, Monoid, Functor, Applicative, Alternative, Monad, MonadPlus) - -instance IsString (SparseT String String) where - fromString = string - -type Sparse = SparseT String - -runSparseT :: SparseT a b -> a -> Maybe b -runSparseT = fmap (fmap snd) . runSparseT' - -runSparseT' :: SparseT a b -> a -> Maybe (a, b) -runSparseT' = getPartialP . getSparseT - -runSparse :: Sparse a -> String -> Maybe a -runSparse = runSparseT - -runSparse' :: Sparse a -> String -> Maybe (String, a) -runSparse' = runSparseT' - ----------- - --- | Consumes one input element. --- --- Fails if the predicate fails, or if there is no more input. --- -headP = SparseT . PartialP . headP' - --- | Consume one or more input elements. --- --- The given function receives the /entire/ remaining input, and must return --- the number of consumed elements. --- --- Fails if the predicate return 0 or less, or if there is no more input. --- -splitP = SparseT . PartialP . splitP' - -headP' :: (a -> Bool) -> [a] -> Maybe ([a], a) -headP' p [] = Nothing -headP' p (x:xs) = if not (p x) then Nothing else Just (xs, x) - -splitP' :: ([a] -> Int) -> [a] -> Maybe ([a], [a]) -splitP' p [] = Nothing -splitP' p ys = let n = p ys in if n < 1 then Nothing else Just (drop n ys, take n ys) - ----------- - -char :: Char -> Sparse Char -char c = charIs (== c) - -charIs :: (Char -> Bool) -> Sparse Char -charIs p = headP p - -string :: String -> Sparse String -string s = stringIs (length s) (== s) - -stringIs :: Int -> (String -> Bool) -> Sparse String -stringIs n p = splitP (\xs -> if p (take n xs) then n else 0) - -asSparse = id -asSparse :: Sparse a -> Sparse a - ----------- - -optionally x p = p <|> return x -optionallyMaybe p = optionally Nothing (liftM Just p) -optional p = do{ p; return ()} <|> return () -between open close p - = do{ open; x <- p; close; return x } -skipMany1 p = do{ p; skipMany p } -skipMany p = scan - where - scan = do{ p; scan } <|> return () -many1 p = do{ x <- p; xs <- many p; return (x:xs) } -sepBy p sep = sepBy1 p sep <|> return [] -sepBy1 p sep = do{ x <- p - ; xs <- many (sep >> p) - ; return (x:xs) - } -sepEndBy1 p sep = do{ x <- p - ; do{ sep - ; xs <- sepEndBy p sep - ; return (x:xs) - } - <|> return [x] - } -sepEndBy p sep = sepEndBy1 p sep <|> return [] -endBy1 p sep = many1 (do{ x <- p; sep; return x }) -endBy p sep = many (do{ x <- p; sep; return x }) -count n p | n <= 0 = return [] - | otherwise = sequence (replicate n p) - ----------- - - --- test :: Sparse [String] -test = asSparse $ string "hans" >> many1 (string ";") - - - - -single x = [x] -list z f xs = case xs of - [] -> z - ys -> f ys - -[a,b,c,d,e,f,g,x,y,z,m,n,o,p,q,r] = undefined diff --git a/src/Sparse/Matrix.hs b/src/Sparse/Matrix.hs new file mode 100644 --- /dev/null +++ b/src/Sparse/Matrix.hs @@ -0,0 +1,410 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE DefaultSignatures #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE KindSignatures #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE PatternGuards #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeSynonymInstances #-} +{-# LANGUAGE UndecidableInstances #-} + +----------------------------------------------------------------------------- +-- | +-- Copyright : (C) 2013 Edward Kmett +-- License : BSD-style (see the file LICENSE) +-- Maintainer : Edward Kmett +-- Stability : experimental +-- Portability : non-portable +-- +-- Sparse Matrices in Morton order +-- +---------------------------------------------------------------------------- + +module Sparse.Matrix + ( + -- * Sparse Matrices + Mat(..) + -- * Keys + , Key(..) + -- * Construction + , Sparse.Matrix.fromList + , Sparse.Matrix.singleton + , transpose + , ident + , empty + -- * Consumption + , size + , null + -- * Distinguishable Zero + , Eq0(..) + -- * Customization + , addWith + , multiplyWith + -- * Lenses + , _Mat, keys, values + ) where + +import Control.Applicative hiding (empty) +import Control.Arrow +import Control.DeepSeq +import Control.Lens +import Data.Bits +import Data.Complex +import Data.Foldable +import Data.Function (on) +import qualified Data.Vector as V +import qualified Data.Vector.Algorithms.Insertion as Sort +import qualified Data.Vector.Generic as G +import qualified Data.Vector.Hybrid as H +import qualified Data.Vector.Hybrid.Internal as H +import qualified Data.Vector.Unboxed as U +import Data.Vector.Fusion.Stream (Stream, sized) +import Data.Vector.Fusion.Stream.Size +import Data.Word +import Prelude hiding (head, last, null) +import Sparse.Matrix.Internal.Fusion as Fusion +import Sparse.Matrix.Internal.Key +import Sparse.Matrix.Internal.Heap as Heap hiding (head) +import Text.Read + +-- import Debug.Trace +-- import Numeric.Lens + +-- * Distinguishable Zero + +class Num a => Eq0 a where + -- | Return whether or not the element is 0. + -- + -- It may be okay to never return 'True', but you won't be + -- able to thin spurious zeroes introduced into your matrix. + -- + isZero :: a -> Bool +#ifndef HLINT + default isZero :: (Num a, Eq a) => a -> Bool + isZero = (0 ==) + {-# INLINE isZero #-} +#endif + + -- | Remove results that are equal to zero from a simpler function. + -- + -- When used with @addWith@ or @multiplyWith@'s additive argument + -- this can help retain the sparsity of the matrix. + nonZero :: (x -> y -> a) -> x -> y -> Maybe a + nonZero f a b = case f a b of + c | isZero c -> Nothing + | otherwise -> Just c + {-# INLINE nonZero #-} + + -- | + -- Add two matrices. By default this assumes 'isZero' can + -- possibly return 'True' after an addition. For some + -- ring-like structures, this doesn't hold. There you can + -- use: + -- + -- @ + -- 'addMats' = 'addWith' ('+') + -- @ + -- + -- By default this will use + -- + -- @ + -- 'addMats' = 'addWith0' '$' 'nonZero' ('+') + -- @ + addMats :: G.Vector v a => Mat v a -> Mat v a -> Mat v a + addMats = addWith0 $ nonZero (+) + {-# INLINE addMats #-} + + -- | Convert from a 'Heap' to a 'Stream'. + -- + -- If addition of non-zero valus in your ring-like structure + -- cannot yield zero, then you can use + -- + -- @ + -- 'addHeap' = 'Heap.streamHeapWith' ('+') + -- @ + -- + -- instead of the default definition: + -- + -- @ + -- 'addHeap' = 'Heap.streamHeapWith0' '$' 'nonZero' ('+') + -- @ + addHeap :: Maybe (Heap a) -> Stream (Key, a) + addHeap = Heap.streamHeapWith0 $ nonZero (+) + +instance Eq0 Int +instance Eq0 Word +instance Eq0 Integer +instance Eq0 Float +instance Eq0 Double +instance (RealFloat a, Eq0 a) => Eq0 (Complex a) where + isZero (a :+ b) = isZero a && isZero b + {-# INLINE isZero #-} + +-- * Sparse Matrices + +-- invariant: all vectors are the same length +data Mat v a = Mat {-# UNPACK #-} !Int !(U.Vector Word) !(U.Vector Word) !(v a) + deriving (Eq,Ord) + +instance (G.Vector v a, Show a) => Show (Mat v a) where + showsPrec d m = G.showsPrec d (m^._Mat) + +instance (G.Vector v a, Read a) => Read (Mat v a) where + readPrec = (_Mat #) <$> G.readPrec + +instance NFData (v a) => NFData (Mat v a) where + rnf (Mat _ xs ys vs) = rnf xs `seq` rnf ys `seq` rnf vs `seq` () + +-- | bundle up the matrix in a form suitable for vector-algorithms +_Mat :: Iso (Mat u a) (Mat v b) (H.Vector U.Vector u (Key, a)) (H.Vector U.Vector v (Key, b)) +_Mat = iso (\(Mat n xs ys vs) -> H.V (V_Key n xs ys) vs) + (\(H.V (V_Key n xs ys) vs) -> Mat n xs ys vs) +{-# INLINE _Mat #-} + +-- | Access the keys of a matrix +keys :: Lens' (Mat v a) (U.Vector Key) +keys f (Mat n xs ys vs) = f (V_Key n xs ys) <&> \ (V_Key n' xs' ys') -> Mat n' xs' ys' vs +{-# INLINE keys #-} + +-- | Access the keys of a matrix +values :: Lens (Mat u a) (Mat v b) (u a) (v b) +values f (Mat n xs ys vs) = Mat n xs ys <$> f vs +{-# INLINE values #-} + +instance Functor v => Functor (Mat v) where + fmap = over (values.mapped) + {-# INLINE fmap #-} + +instance Foldable v => Foldable (Mat v) where + foldMap = foldMapOf (values.folded) + {-# INLINE foldMap #-} + +instance Traversable v => Traversable (Mat v) where + traverse = values.traverse + {-# INLINE traverse #-} + +type instance IxValue (Mat v a) = a +type instance Index (Mat v a) = Key + +-- traverse a Vector +eachV :: (Applicative f, G.Vector v a, G.Vector v b) => (a -> f b) -> v a -> f (v b) +eachV f v = G.fromListN (G.length v) <$> traverse f (G.toList v) + +instance (Applicative f, G.Vector v a, G.Vector v b) => Each f (Mat v a) (Mat v b) a b where + each f = _Mat $ eachV $ \(k,v) -> (,) k <$> indexed f k v + {-# INLINE each #-} + +instance (Functor f, Contravariant f, G.Vector v a) => Contains f (Mat v a) where + contains = containsIx + +instance (Applicative f, G.Vector v a) => Ixed f (Mat v a) where + ix ij@(Key i j) f m@(Mat n xs ys vs) + | Just i' <- xs U.!? l, i == i' + , Just j' <- ys U.!? l, j == j' = indexed f ij (vs G.! l) <&> \v -> Mat n xs ys (vs G.// [(l,v)]) + | otherwise = pure m + where l = search (\k -> Key (xs U.! k) (ys U.! k) >= ij) 0 n + {-# INLINE ix #-} + +instance (G.Vector v a, Num a, Eq0 a) => Eq0 (Mat v a) where + isZero (Mat n _ _ _) = n == 0 + {-# INLINE isZero #-} + +-- * Construction + +-- | Build a sparse matrix. +fromList :: G.Vector v a => [(Key, a)] -> Mat v a +fromList xs = _Mat # H.modify (Sort.sortBy (compare `on` fst)) (H.fromList xs) +{-# INLINABLE fromList #-} + +-- | Transpose a matrix +transpose :: G.Vector v a => Mat v a -> Mat v a +transpose xs = xs & _Mat %~ H.modify (Sort.sortBy (compare `on` fst)) . H.map (first swap) +{-# INLINE transpose #-} + +-- | @singleton@ makes a matrix with a singleton value at a given location +singleton :: G.Vector v a => Key -> a -> Mat v a +singleton k v = _Mat # H.singleton (k,v) +{-# INLINE singleton #-} + +-- | @ident n@ makes an @n@ x @n@ identity matrix +-- +-- >>> ident 4 :: Mat U.Vector Int +-- fromList [(Key 0 0,1),(Key 1 1,1),(Key 2 2,1),(Key 3 3,1)] +ident :: (G.Vector v a, Num a) => Int -> Mat v a +ident w = Mat w (U.generate w fromIntegral) (U.generate w fromIntegral) (G.replicate w 1) +{-# INLINE ident #-} + +-- | The empty matrix +-- +-- >>> empty :: Mat U.Vector Int +-- fromList [] +empty :: G.Vector v a => Mat v a +empty = Mat 0 U.empty U.empty G.empty +{-# INLINE empty #-} + +-- * Consumption + +-- | Count the number of non-zero entries in the matrix +-- +-- >>> size (ident 4 :: Mat U.Vector Int) +-- 4 +size :: Mat v a -> Int +size (Mat n _ _ _) = n +{-# INLINE size #-} + +-- | +-- >>> null (empty :: Mat U.Vector Int) +-- True +null :: Mat v a -> Bool +null (Mat n _ _ _) = n == 0 +{-# INLINE null #-} + +instance (G.Vector v a, Num a, Eq0 a) => Num (Mat v a) where + {-# SPECIALIZE instance (Num a, Eq0 a) => Num (Mat V.Vector a) #-} + {-# SPECIALIZE instance Num (Mat U.Vector Int) #-} + {-# SPECIALIZE instance Num (Mat U.Vector Double) #-} + {-# SPECIALIZE instance Num (Mat U.Vector (Complex Double)) #-} + abs = over each abs + {-# INLINE abs #-} + signum = over each signum + {-# INLINE signum #-} + negate = over each negate + {-# INLINE negate #-} + fromInteger 0 = empty + fromInteger _ = error "Mat: fromInteger n" + {-# INLINE fromInteger #-} + (+) = addMats + {-# INLINE (+) #-} + (-) = addWith0 $ nonZero (-) + {-# INLINE (-) #-} + (*) = multiplyWith (*) addHeap + {-# INLINEABLE (*) #-} + +-- * Utilities + +-- | assuming @l <= h@. Returns @h@ if the predicate is never @True@ over @[l..h)@ +search :: (Int -> Bool) -> Int -> Int -> Int +search p = go where + go l h + | l == h = l + | p m = go l m + | otherwise = go (m+1) h + where m = l + div (h-l) 2 +{-# INLINE search #-} + +split1 :: G.Vector v a => Word -> Word -> Mat v a -> (Mat v a, Mat v a) +split1 ai bi (Mat n xs ys vs) = (m0,m1) + where + !aibi = xor ai bi + !k = search (\l -> xor (xs U.! l) bi `lts` aibi) 0 n + (xs0,xs1) = U.splitAt k xs + (ys0,ys1) = U.splitAt k ys + (vs0,vs1) = G.splitAt k vs + !m0 = Mat k xs0 ys0 vs0 + !m1 = Mat (n-k) xs1 ys1 vs1 +{-# INLINE split1 #-} + +split2 :: G.Vector v a => Word -> Word -> Mat v a -> (Mat v a, Mat v a) +split2 aj bj (Mat n xs ys vs) = (m0,m1) + where + !ajbj = xor aj bj + !k = search (\l -> xor (ys U.! l) bj `lts` ajbj) 0 n + (xs0,xs1) = U.splitAt k xs + (ys0,ys1) = U.splitAt k ys + (vs0,vs1) = G.splitAt k vs + !m0 = Mat k xs0 ys0 vs0 + !m1 = Mat (n-k) xs1 ys1 vs1 +{-# INLINE split2 #-} + +-- | Merge two matrices where the indices coincide into a new matrix. This provides for generalized +-- addition, but where the summation of two non-zero entries is necessarily non-zero. +addWith :: G.Vector v a => (a -> a -> a) -> Mat v a -> Mat v a -> Mat v a +addWith f xs ys = _Mat # G.unstream (mergeStreamsWith f (G.stream (xs^._Mat)) (G.stream (ys^._Mat))) +{-# INLINE addWith #-} + +-- | Merge two matrices where the indices coincide into a new matrix. This provides for generalized +-- addition. Return 'Nothing' for zero. +addWith0 :: G.Vector v a => (a -> a -> Maybe a) -> Mat v a -> Mat v a -> Mat v a +addWith0 f xs ys = _Mat # G.unstream (mergeStreamsWith0 f (G.stream (xs^._Mat)) (G.stream (ys^._Mat))) +{-# INLINE addWith0 #-} + +-- | Multiply two matrices using the specified multiplication and addition operation. +multiplyWith :: G.Vector v a => (a -> a -> a) -> (Maybe (Heap a) -> Stream (Key, a)) -> Mat v a -> Mat v a -> Mat v a +{-# INLINEABLE multiplyWith #-} +multiplyWith times make x0 y0 = case compare (size x0) 1 of + LT -> empty + EQ | size y0 == 1 -> _Mat # (G.unstream $ hint $ make $ go11 (lo x0) (head x0) (lo y0) (head y0)) + | otherwise -> _Mat # (G.unstream $ hint $ make $ go12 (lo x0) (head x0) (lo y0) y0 (hi y0)) + GT -> case compare (size y0) 1 of + LT -> empty + EQ -> _Mat # (G.unstream $ hint $ make $ go21 (lo x0) x0 (hi x0) (lo y0) (head y0)) + GT -> _Mat # (G.unstream $ hint $ make $ go22 (lo x0) x0 (hi x0) (lo y0) y0 (hi y0)) + where + hint x = sized x $ Max (size x0 * size y0) + go11 (Key i j) a (Key j' k) b + | j == j' = Just $ Heap.singleton (Key i k) (times a b) + | otherwise = Nothing + + -- internal cases in go22 + go22L0 xa x ya y yb + | size x == 1 = go12 xa (head x) ya y yb + | otherwise = go22 xa x (hi x) ya y yb + {-# INLINE go22L0 #-} + + go22L1 x xb ya y yb + | size x == 1 = go12 xb (head x) ya y yb + | otherwise = go22 (lo x) x xb ya y yb + {-# INLINE go22L1 #-} + + go22R0 xa x xb ya y + | size y == 1 = go21 xa x xb ya (head y) + | otherwise = go22 xa x xb ya y (hi y) + {-# INLINE go22R0 #-} + + go22R1 xa x xb y yb + | size y == 1 = go21 xa x xb yb (head y) + | otherwise = go22 xa x xb (lo y) y yb + {-# INLINE go22R1 #-} + + -- x and y have at least 2 non-zero elements each + go22 xa@(Key xai xaj) x xb@(Key xbi xbj) ya@(Key yaj yak) y yb@(Key ybj ybk) + | gts (xor xaj yaj) (xiyj .|. ykxj) = Nothing + | ges xiyj ykxj + = if ges xi yj then case split1 xai xbi x of (m0,m1) -> go22L0 xa m0 ya y yb `mfby` go22L1 m1 xb ya y yb -- we can split on i, fby + else case split1 yaj ybj y of (m0,m1) -> go22R0 xa x xb ya m0 `madd` go22R1 xa x xb m1 yb -- we split on j, mix + | ges yk xj = case split2 yak ybk y of (m0,m1) -> go22R0 xa x xb ya m0 `mfby` go22R1 xa x xb m1 yb -- we can split on k, fby + | otherwise = case split2 xaj xbj x of (m0,m1) -> go22L0 xa m0 ya y yb `madd` go22L1 m1 xb ya y yb -- we split on j, mix + where + xi = xor xai xbi + xj = xor xaj xbj + yj = xor yaj ybj + yk = xor yak ybk + xiyj = xi .|. yj + ykxj = yk .|. xj + + go21 _ mx _ yb b = Heap.timesSingleton times (G.stream (mx^._Mat)) yb b -- linear scan. use tree and fast rejects? + go12 xa a _ my _ = Heap.singletonTimes times xa a (G.stream (my^._Mat)) + + madd Nothing xs = xs + madd xs Nothing = xs + madd (Just x) (Just y) = Just (mix x y) + {-# INLINE madd #-} + + mfby Nothing xs = xs + mfby xs Nothing = xs + mfby (Just x) (Just y) = Just (fby x y) + {-# INLINE mfby #-} + + lo (Mat _ xs ys _) = Key (U.head xs) (U.head ys) + {-# INLINE lo #-} + + hi (Mat _ xs ys _) = Key (U.last xs) (U.last ys) + {-# INLINE hi #-} + + head :: G.Vector v a => Mat v a -> a + head (Mat _ _ _ vs) = G.head vs + {-# INLINE head #-} diff --git a/src/Sparse/Matrix/Internal/Fusion.hs b/src/Sparse/Matrix/Internal/Fusion.hs new file mode 100644 --- /dev/null +++ b/src/Sparse/Matrix/Internal/Fusion.hs @@ -0,0 +1,122 @@ +{-# LANGUAGE BangPatterns #-} +----------------------------------------------------------------------------- +-- | +-- Copyright : (C) 2013 Edward Kmett +-- License : BSD-style (see the file LICENSE) +-- Maintainer : Edward Kmett +-- Stability : experimental +-- Portability : non-portable +-- +-- Matrix stream fusion internals +-- +----------------------------------------------------------------------------- +module Sparse.Matrix.Internal.Fusion + ( mergeStreamsWith, mergeStreamsWith0 + ) where + +import Data.Vector.Fusion.Stream.Monadic (Step(..), Stream(..)) +import Data.Vector.Fusion.Stream.Size +import Sparse.Matrix.Internal.Key + +-- | The state for 'Stream' fusion that is used by 'mergeStreamsWith'. +-- +-- This form permits cancellative addition. +data MergeState sa sb i a + = MergeL sa sb i a + | MergeR sa sb i a + | MergeLeftEnded sb + | MergeRightEnded sa + | MergeStart sa sb + +-- | This is the internal stream fusion combinator used to merge streams for addition. +-- +-- This form permits cancellative addition. +mergeStreamsWith0 :: Monad m => (a -> a -> Maybe a) -> Stream m (Key, a) -> Stream m (Key, a) -> Stream m (Key, a) +mergeStreamsWith0 f (Stream stepa sa0 na) (Stream stepb sb0 nb) + = Stream step (MergeStart sa0 sb0) (toMax na + toMax nb) where + step (MergeStart sa sb) = do + r <- stepa sa + return $ case r of + Yield (i, a) sa' -> Skip (MergeL sa' sb i a) + Skip sa' -> Skip (MergeStart sa' sb) + Done -> Skip (MergeLeftEnded sb) + step (MergeL sa sb i a) = do + r <- stepb sb + return $ case r of + Yield (j, b) sb' -> case compare i j of + LT -> Yield (i, a) (MergeR sa sb' j b) + EQ -> case f a b of + Just c -> Yield (i, c) (MergeStart sa sb') + Nothing -> Skip (MergeStart sa sb') + GT -> Yield (j, b) (MergeL sa sb' i a) + Skip sb' -> Skip (MergeL sa sb' i a) + Done -> Yield (i, a) (MergeRightEnded sa) + step (MergeR sa sb j b) = do + r <- stepa sa + return $ case r of + Yield (i, a) sa' -> case compare i j of + LT -> Yield (i, a) (MergeR sa' sb j b) + EQ -> case f a b of + Just c -> Yield (i, c) (MergeStart sa' sb) + Nothing -> Skip (MergeStart sa' sb) + GT -> Yield (j, b) (MergeL sa' sb i a) + Skip sa' -> Skip (MergeR sa' sb j b) + Done -> Yield (j, b) (MergeLeftEnded sb) + step (MergeLeftEnded sb) = do + r <- stepb sb + return $ case r of + Yield (j, b) sb' -> Yield (j, b) (MergeLeftEnded sb') + Skip sb' -> Skip (MergeLeftEnded sb') + Done -> Done + step (MergeRightEnded sa) = do + r <- stepa sa + return $ case r of + Yield (i, a) sa' -> Yield (i, a) (MergeRightEnded sa') + Skip sa' -> Skip (MergeRightEnded sa') + Done -> Done + {-# INLINE [0] step #-} +{-# INLINE [1] mergeStreamsWith0 #-} + + +-- | This is the internal stream fusion combinator used to merge streams for addition. +mergeStreamsWith :: Monad m => (a -> a -> a) -> Stream m (Key, a) -> Stream m (Key, a) -> Stream m (Key, a) +mergeStreamsWith f (Stream stepa sa0 na) (Stream stepb sb0 nb) + = Stream step (MergeStart sa0 sb0) (toMax na + toMax nb) where + step (MergeStart sa sb) = do + r <- stepa sa + return $ case r of + Yield (i, a) sa' -> Skip (MergeL sa' sb i a) + Skip sa' -> Skip (MergeStart sa' sb) + Done -> Skip (MergeLeftEnded sb) + step (MergeL sa sb i a) = do + r <- stepb sb + return $ case r of + Yield (j, b) sb' -> case compare i j of + LT -> Yield (i, a) (MergeR sa sb' j b) + EQ -> Yield (i, f a b) (MergeStart sa sb') + GT -> Yield (j, b) (MergeL sa sb' i a) + Skip sb' -> Skip (MergeL sa sb' i a) + Done -> Yield (i, a) (MergeRightEnded sa) + step (MergeR sa sb j b) = do + r <- stepa sa + return $ case r of + Yield (i, a) sa' -> case compare i j of + LT -> Yield (i, a) (MergeR sa' sb j b) + EQ -> Yield (i, f a b) (MergeStart sa' sb) + GT -> Yield (j, b) (MergeL sa' sb i a) + Skip sa' -> Skip (MergeR sa' sb j b) + Done -> Yield (j, b) (MergeLeftEnded sb) + step (MergeLeftEnded sb) = do + r <- stepb sb + return $ case r of + Yield (j, b) sb' -> Yield (j, b) (MergeLeftEnded sb') + Skip sb' -> Skip (MergeLeftEnded sb') + Done -> Done + step (MergeRightEnded sa) = do + r <- stepa sa + return $ case r of + Yield (i, a) sa' -> Yield (i, a) (MergeRightEnded sa') + Skip sa' -> Skip (MergeRightEnded sa') + Done -> Done + {-# INLINE [0] step #-} +{-# INLINE [1] mergeStreamsWith #-} diff --git a/src/Sparse/Matrix/Internal/Heap.hs b/src/Sparse/Matrix/Internal/Heap.hs new file mode 100644 --- /dev/null +++ b/src/Sparse/Matrix/Internal/Heap.hs @@ -0,0 +1,211 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE PatternGuards #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE TypeFamilies #-} +----------------------------------------------------------------------------- +-- | +-- Copyright : (C) 2013 Edward Kmett +-- License : BSD-style (see the file LICENSE) +-- Maintainer : Edward Kmett +-- Stability : experimental +-- Portability : non-portable +-- +-- Bootstrapped catenable non-empty pairing heaps as described in +-- +-- +----------------------------------------------------------------------------- +module Sparse.Matrix.Internal.Heap + ( Heap(..) + , fby + , mix + , singleton + , head + , tail + , fromList + , fromAscList + , streamHeapWith + , streamHeapWith0 + , timesSingleton + , singletonTimes + ) where + +import Control.Applicative +import Control.Lens +import Data.Foldable +import Data.Monoid +import Data.Vector.Fusion.Stream.Monadic hiding (singleton, fromList, head, tail) +import Data.Vector.Fusion.Stream.Size +import Data.Vector.Fusion.Util +import Sparse.Matrix.Internal.Key +import Prelude hiding (head, tail) + +-- | Bootstrapped _catenable_ non-empty pairing heaps +data Heap a = Heap {-# UNPACK #-} !Key a [Heap a] [Heap a] [Heap a] + deriving (Show,Read) + +-- | Append two heaps where we know every key in the first occurs before every key in the second +-- +-- >>> head $ singleton (Key 1 1) 1 `fby` singleton (Key 2 2) 2 +-- (Key 1 1,1) +fby :: Heap a -> Heap a -> Heap a +fby (Heap i a as ls rs) r = Heap i a as ls (r:rs) + +-- | Interleave two heaps making a new 'Heap' +-- +-- >>> head $ singleton (Key 1 1) 1 `mix` singleton (Key 2 2) 2 +-- (Key 1 1,1) +mix :: Heap a -> Heap a -> Heap a +mix x@(Heap i a as al ar) y@(Heap j b bs bl br) + | i <= j = Heap i a (y:pops as al ar) [] [] + | otherwise = Heap j b (x:pops bs bl br) [] [] + +-- | +-- >>> head $ singleton (Key 1 1) 1 +-- (Key 1 1,1) +head :: Heap a -> (Key, a) +head (Heap i a _ _ _) = (i, a) + +-- | +-- >>> tail $ singleton (Key 1 1) 1 +-- Nothing +tail :: Heap a -> Maybe (Heap a) +tail (Heap _ _ xs fs rs) = pop xs fs rs + +-- | +-- >>> singleton (Key 1 1) 1 +-- Heap (Key 1 1) 1 [] [] [] +singleton :: Key -> a -> Heap a +singleton k v = Heap k v [] [] [] + +-- | Build a 'Heap' from a jumbled up list of elements. +fromList :: [(Key,a)] -> Heap a +fromList ((k0,v0):xs) = Prelude.foldr (\(k,v) r -> mix (singleton k v) r) (singleton k0 v0) xs +fromList [] = error "empty Heap" + +-- | Build a 'Heap' from an list of elements that must be in strictly ascending Morton order. +fromAscList :: [(Key,a)] -> Heap a +fromAscList ((k0,v0):xs) = Prelude.foldr (\(k,v) r -> fby (singleton k v) r) (singleton k0 v0) xs +fromAscList [] = error "empty Heap" + +-- * Internals + +fbys :: Heap a -> [Heap a] -> [Heap a] -> Heap a +fbys (Heap i a as [] []) ls' rs' = Heap i a as ls' rs' +fbys (Heap i a as ls []) ls' rs' = Heap i a as ls $ rs' <> reverse ls' +fbys (Heap i a as ls rs) ls' rs' = Heap i a as ls $ rs' <> reverse ls' <> rs + +pops :: [Heap a] -> [Heap a] -> [Heap a] -> [Heap a] +pops xs [] [] = xs +pops (x:xs) ls rs = [fbys (Prelude.foldl mix x xs) ls rs] +pops [] (l:ls) rs = [fbys l ls rs] +pops [] [] rs = case reverse rs of + f:fs -> [fbys f fs []] + _ -> [] -- caught above by the 'go as [] []' case + +pop :: [Heap a] -> [Heap a] -> [Heap a] -> Maybe (Heap a) +pop (x:xs) ls rs = Just $ fbys (Prelude.foldl mix x xs) ls rs +pop [] (l:ls) rs = Just $ fbys l ls rs +pop [] [] rs = case reverse rs of + f:fs -> Just (fbys f fs []) + [] -> Nothing + +-- * Instances + +instance Functor Heap where + fmap f (Heap k a xs ls rs) = Heap k (f a) (fmap f <$> xs) (fmap f <$> ls) (fmap f <$> rs) + +instance FunctorWithIndex Key Heap where + imap f (Heap k a xs ls rs) = Heap k (f k a) (imap f <$> xs) (imap f <$> ls) (imap f <$> rs) + +instance Foldable Heap where + foldMap f = go where + go (Heap _ a xs ls rs) = case pop xs ls rs of + Nothing -> f a + Just h -> f a `mappend` go h + {-# INLINE foldMap #-} + +instance FoldableWithIndex Key Heap where + ifoldMap f = go where + go (Heap i a xs ls rs) = case pop xs ls rs of + Nothing -> f i a + Just h -> f i a `mappend` go h + {-# INLINE ifoldMap #-} + +instance Traversable Heap where + traverse f xs = fromAscList <$> traverse (traverse f) (itoList xs) + {-# INLINE traverse #-} + +instance TraversableWithIndex Key Heap where + itraverse f xs = fromAscList <$> traverse (\(k,v) -> (,) k <$> f k v) (itoList xs) + {-# INLINE itraverse #-} + +data HeapState a + = Start !(Heap a) + | Ready {-# UNPACK #-} !Key a !(Heap a) + | Final {-# UNPACK #-} !Key a + | Finished + +-- | Convert a 'Heap' into a 'Stream' folding together values with identical keys using the supplied +-- addition operator. +streamHeapWith :: Monad m => (a -> a -> a) -> Maybe (Heap a) -> Stream m (Key, a) +streamHeapWith f h0 = Stream step (maybe Finished Start h0) Unknown where + step (Start (Heap i a xs ls rs)) = return $ Skip $ maybe (Final i a) (Ready i a) $ pop xs ls rs + step (Ready i a (Heap j b xs ls rs)) = return $ case compare i j of + LT -> Yield (i, a) $ maybe (Final j b) (Ready j b) $ pop xs ls rs + EQ | c <- f a b -> Skip $ maybe (Final i c) (Ready i c) $ pop xs ls rs + GT -> Yield (j, b) $ maybe (Final i a) (Ready i a) $ pop xs ls rs + step (Final i a) = return $ Yield (i,a) Finished + step Finished = return Done + {-# INLINE [1] step #-} +{-# INLINE [0] streamHeapWith #-} + +-- | Convert a 'Heap' into a 'Stream' folding together values with identical keys using the supplied +-- addition operator that is allowed to return a sparse 0, by returning 'Nothing'. +streamHeapWith0 :: Monad m => (a -> a -> Maybe a) -> Maybe (Heap a) -> Stream m (Key, a) +streamHeapWith0 f h0 = Stream step (maybe Finished Start h0) Unknown where + step (Start (Heap i a xs ls rs)) = return $ Skip $ maybe (Final i a) (Ready i a) $ pop xs ls rs + step (Ready i a (Heap j b xs ls rs)) = return $ case compare i j of + LT -> Yield (i, a) $ maybe (Final j b) (Ready j b) $ pop xs ls rs + EQ -> case f a b of + Nothing -> Skip $ maybe Finished Start $ pop xs ls rs + Just c -> Skip $ maybe (Final i c) (Ready i c) $ pop xs ls rs + GT -> Yield (j, b) $ maybe (Final i a) (Ready i a) $ pop xs ls rs + step (Final i a) = return $ Yield (i,a) Finished + step Finished = return Done + {-# INLINE [1] step #-} +{-# INLINE [0] streamHeapWith0 #-} + +-- | This is an internal 'Heap' fusion combinator used to multiply on the right by a singleton 'Key'/value pair. +timesSingleton :: (a -> b -> c) -> Stream Id (Key, a) -> Key -> b -> Maybe (Heap c) +timesSingleton f (Stream stepa sa0 _) (Key j k) b = start sa0 where + start sa = case unId (stepa sa) of + Yield (Key i j', a) sa' + | j == j' -> Just $ run (singleton (Key i k) (f a b)) sa' + | otherwise -> start sa' + Skip sa' -> start sa' + Done -> Nothing + run h sa = case unId (stepa sa) of + Yield (Key i j', a) sa' + | j == j' -> run (h `mix` singleton (Key i k) (f a b)) sa' + | otherwise -> run h sa' + Skip sa' -> run h sa' + Done -> h +{-# INLINE timesSingleton #-} + +-- | This is an internal 'Heap' fusion combinator used to multiply on the right by a singleton 'Key'/value pair. +singletonTimes :: (a -> b -> c) -> Key -> a -> Stream Id (Key, b) -> Maybe (Heap c) +singletonTimes f (Key i j) a (Stream stepb sb0 _) = start sb0 where + start sb = case unId (stepb sb) of + Yield (Key j' k, b) sb' + | j == j' -> Just $ run (singleton (Key i k) (f a b)) sb' + | otherwise -> start sb' + Skip sb' -> start sb' + Done -> Nothing + run h sb = case unId (stepb sb) of + Yield (Key j' k, b) sb' + | j == j' -> run (h `mix` singleton (Key i k) (f a b)) sb' + | otherwise -> run h sb' + Skip sb' -> run h sb' + Done -> h +{-# INLINE singletonTimes #-} diff --git a/src/Sparse/Matrix/Internal/Key.hs b/src/Sparse/Matrix/Internal/Key.hs new file mode 100644 --- /dev/null +++ b/src/Sparse/Matrix/Internal/Key.hs @@ -0,0 +1,253 @@ +{-# LANGUAGE KindSignatures #-} +{-# LANGUAGE DefaultSignatures #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeSynonymInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE GADTs #-} +----------------------------------------------------------------------------- +-- | +-- Copyright : (C) 2013 Edward Kmett +-- License : BSD-style (see the file LICENSE) +-- Maintainer : Edward Kmett +-- Stability : experimental +-- Portability : non-portable +-- +-- Keys in Morton order +-- +-- This module provides combinators for shuffling together the bits of two +-- key components to get a key that is based on their interleaved bits. +-- +-- See for more information +-- about Morton order. +-- +-- How to perform the comparison without interleaving is described in +-- +-- +-- +---------------------------------------------------------------------------- +module Sparse.Matrix.Internal.Key + ( + -- * Keys in Morton order + Key(..) + , swap + -- * Most significant bit comparisons + , compares + , lts, les, eqs, nes, ges, gts + -- * Unboxed vector constructors + , U.MVector(..) + , U.Vector(..) + ) where + +import Data.Bits +import Control.Monad +import Control.Lens +import qualified Data.Vector.Generic as G +import qualified Data.Vector.Generic.Mutable as GM +import qualified Data.Vector.Unboxed as U +import Data.Word + +-- * Morton Order + +-- | @Key i j@ logically orders the keys as if the bits of the keys @i@ and @j@ +-- were interleaved. This is equivalent to storing the keys in \"Morton Order\". +-- +-- >>> Key 100 200 ^. _1 +-- 100 +-- +-- >>> Key 100 200 ^. _2 +-- 200 +data Key = Key {-# UNPACK #-} !Word {-# UNPACK #-} !Word + deriving (Show, Read, Eq) + +instance Ord Key where + Key a b `compare` Key c d + | xor a c `lts` xor b d = compare b d + | otherwise = compare a c + +instance (a ~ Word, b ~ Word) => Field1 Key Key a b where + _1 f (Key i j) = indexed f (0 :: Int) i <&> (Key ?? j) + {-# INLINE _1 #-} + +instance (a ~ Word, b ~ Word) => Field2 Key Key a b where + _2 f (Key i j) = indexed f (1 :: Int) j <&> Key i + {-# INLINE _2 #-} + +instance U.Unbox Key + +data instance U.MVector s Key = MV_Key {-# UNPACK #-} !Int !(U.MVector s Word) !(U.MVector s Word) +data instance U.Vector Key = V_Key {-# UNPACK #-} !Int !(U.Vector Word) !(U.Vector Word) + +instance GM.MVector U.MVector Key where + {-# INLINE basicLength #-} + {-# INLINE basicUnsafeSlice #-} + {-# INLINE basicOverlaps #-} + {-# INLINE basicUnsafeNew #-} + {-# INLINE basicUnsafeReplicate #-} + {-# INLINE basicUnsafeRead #-} + {-# INLINE basicUnsafeWrite #-} + {-# INLINE basicClear #-} + {-# INLINE basicSet #-} + {-# INLINE basicUnsafeCopy #-} + {-# INLINE basicUnsafeGrow #-} + basicLength (MV_Key l _ _) = l + basicUnsafeSlice i n (MV_Key _ u v) = MV_Key n (GM.basicUnsafeSlice i n u) (GM.basicUnsafeSlice i n v) + basicOverlaps (MV_Key _ u1 v1) (MV_Key _ u2 v2) = GM.basicOverlaps u1 u2 || GM.basicOverlaps v1 v2 + basicUnsafeNew n = liftM2 (MV_Key n) (GM.basicUnsafeNew n) (GM.basicUnsafeNew n) + basicUnsafeReplicate n (Key x y) = liftM2 (MV_Key n) (GM.basicUnsafeReplicate n x) (GM.basicUnsafeReplicate n y) + basicUnsafeRead (MV_Key _ u v) i = liftM2 Key (GM.basicUnsafeRead u i) (GM.basicUnsafeRead v i) + basicUnsafeWrite (MV_Key _ u v) i (Key x y) = GM.basicUnsafeWrite u i x >> GM.basicUnsafeWrite v i y + basicClear (MV_Key _ u v) = GM.basicClear u >> GM.basicClear v + basicSet (MV_Key _ u v) (Key x y) = GM.basicSet u x >> GM.basicSet v y + basicUnsafeCopy (MV_Key _ u1 v1) (MV_Key _ u2 v2) = GM.basicUnsafeCopy u1 u2 >> GM.basicUnsafeCopy v1 v2 + basicUnsafeMove (MV_Key _ u1 v1) (MV_Key _ u2 v2) = GM.basicUnsafeMove u1 u2 >> GM.basicUnsafeMove v1 v2 + basicUnsafeGrow (MV_Key _ u v) n = liftM2 (MV_Key n) (GM.basicUnsafeGrow u n) (GM.basicUnsafeGrow v n) + +instance G.Vector U.Vector Key where + {-# INLINE basicLength #-} + {-# INLINE basicUnsafeFreeze #-} + {-# INLINE basicUnsafeThaw #-} + {-# INLINE basicUnsafeSlice #-} + {-# INLINE basicUnsafeIndexM #-} + {-# INLINE elemseq #-} + basicLength (V_Key v _ _) = v + basicUnsafeFreeze (MV_Key n u v) = liftM2 (V_Key n) (G.basicUnsafeFreeze u) (G.basicUnsafeFreeze v) + basicUnsafeThaw (V_Key n u v) = liftM2 (MV_Key n) (G.basicUnsafeThaw u) (G.basicUnsafeThaw v) + basicUnsafeSlice i n (V_Key _ u v) = V_Key n (G.basicUnsafeSlice i n u) (G.basicUnsafeSlice i n v) + basicUnsafeIndexM (V_Key _ u v) i = liftM2 Key (G.basicUnsafeIndexM u i) (G.basicUnsafeIndexM v i) + basicUnsafeCopy (MV_Key _ mu mv) (V_Key _ u v) = G.basicUnsafeCopy mu u >> G.basicUnsafeCopy mv v + elemseq _ (Key x y) z = G.elemseq (undefined :: U.Vector Word) x + $ G.elemseq (undefined :: U.Vector Word) y z + +-- | Swaps the key components around +-- +-- >>> swap (Key 100 200) +-- Key 200 100 +swap :: Key -> Key +swap (Key i j) = Key j i +{-# INLINE swap #-} + +-- | compare the position of the most significant bit of two words +-- +-- >>> compares 4 7 +-- EQ +-- +-- >>> compares 7 9 +-- LT +-- +-- >>> compares 9 7 +-- GT +compares :: Word -> Word -> Ordering +compares a b = case compare a b of + LT | a < xor a b -> LT + GT | b < xor a b -> GT + _ -> EQ +{-# INLINE compares #-} + +-- | @'lts' a b@ returns 'True' when the position of the most significant bit of @a@ is less than the position of the most signficant bit of @b@. +-- +-- >>> lts 4 10 +-- True +-- +-- >>> lts 4 7 +-- False +-- +-- >>> lts 7 8 +-- True +lts :: Word -> Word -> Bool +lts a b = a < b && a < xor a b +{-# INLINE lts #-} + +-- | @'les' a b@ returns 'True' when the position of the most significant bit of @a@ is less than or equal to the position of the most signficant bit of @b@. +-- +-- >>> les 4 10 +-- True +-- +-- >>> les 4 7 +-- True +-- +-- >>> les 7 4 +-- True +-- +-- >>> les 10 4 +-- False +les :: Word -> Word -> Bool +les a b = a <= b || xor a b <= b +{-# INLINE les #-} + +-- | @'eqs' a b@ returns 'True' when the position of the most significant bit of @a@ is equal to the position of the most signficant bit of @b@. +-- +-- >>> eqs 4 7 +-- True +-- +-- >>> eqs 4 8 +-- False +-- +-- >>> eqs 7 4 +-- True +-- +-- >>> eqs 8 4 +-- False +eqs :: Word -> Word -> Bool +eqs a b = case compare a b of + LT -> a >= xor a b + GT -> b >= xor a b + EQ -> True +{-# INLINE eqs #-} + +-- | @'nes' a b@ returns 'True' when the position of the most significant bit of @a@ is not equal to the position of the most signficant bit of @b@. +-- +-- >>> nes 4 7 +-- False +-- +-- >>> nes 4 8 +-- True +-- +-- >>> nes 7 4 +-- False +-- +-- >>> nes 8 4 +-- True +nes :: Word -> Word -> Bool +nes a b = case compare a b of + LT -> a < xor a b + GT -> b < xor a b + EQ -> False +{-# INLINE nes #-} + +-- | @'gts' a b@ returns 'True' when the position of the most significant bit of @a@ is greater than to the position of the most signficant bit of @b@. +-- +-- >>> gts 4 10 +-- False +-- +-- >>> gts 4 7 +-- False +-- +-- >>> gts 7 4 +-- False +-- +-- >>> gts 10 4 +-- True +gts :: Word -> Word -> Bool +gts a b = a > b && xor a b > b +{-# INLINE gts #-} + +-- | @'gts' a b@ returns 'True' when the position of the most significant bit of @a@ is greater than or equal to the position of the most signficant bit of @b@. +-- +-- >>> ges 4 10 +-- False +-- +-- >>> ges 4 7 +-- True +-- +-- >>> ges 7 4 +-- True +-- +-- >>> ges 10 4 +-- True +ges :: Word -> Word -> Bool +ges a b = a >= b || a >= xor a b +{-# INLINE ges #-} diff --git a/tests/doctests.hsc b/tests/doctests.hsc new file mode 100644 --- /dev/null +++ b/tests/doctests.hsc @@ -0,0 +1,73 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE ForeignFunctionInterface #-} +----------------------------------------------------------------------------- +-- | +-- Module : Main (doctests) +-- Copyright : (C) 2012-13 Edward Kmett +-- License : BSD-style (see the file LICENSE) +-- Maintainer : Edward Kmett +-- Stability : provisional +-- Portability : portable +-- +-- This module provides doctests for a project based on the actual versions +-- of the packages it was built with. It requires a corresponding Setup.lhs +-- to be added to the project +----------------------------------------------------------------------------- +module Main where + +import Build_doctests (deps) +import Control.Applicative +import Control.Monad +import Data.List +import System.Directory +import System.FilePath +import Test.DocTest + +##if defined(mingw32_HOST_OS) +##if defined(i386_HOST_ARCH) +##define USE_CP +import Control.Applicative +import Control.Exception +import Foreign.C.Types +foreign import stdcall "windows.h SetConsoleCP" c_SetConsoleCP :: CUInt -> IO Bool +foreign import stdcall "windows.h GetConsoleCP" c_GetConsoleCP :: IO CUInt +##elif defined(x86_64_HOST_ARCH) +##define USE_CP +import Control.Applicative +import Control.Exception +import Foreign.C.Types +foreign import ccall "windows.h SetConsoleCP" c_SetConsoleCP :: CUInt -> IO Bool +foreign import ccall "windows.h GetConsoleCP" c_GetConsoleCP :: IO CUInt +##endif +##endif + +-- | Run in a modified codepage where we can print UTF-8 values on Windows. +withUnicode :: IO a -> IO a +##ifdef USE_CP +withUnicode m = do + cp <- c_GetConsoleCP + (c_SetConsoleCP 65001 >> m) `finally` c_SetConsoleCP cp +##else +withUnicode m = m +##endif + +main :: IO () +main = withUnicode $ getSources >>= \sources -> doctest $ + "-isrc" + : "-idist/build/autogen" + : "-optP-include" + : "-optPdist/build/autogen/cabal_macros.h" + : "-hide-all-packages" + : map ("-package="++) deps ++ sources + +getSources :: IO [FilePath] +getSources = filter (isSuffixOf ".hs") <$> go "src" + where + go dir = do + (dirs, files) <- getFilesAndDirectories dir + (files ++) . concat <$> mapM go dirs + +getFilesAndDirectories :: FilePath -> IO ([FilePath], [FilePath]) +getFilesAndDirectories dir = do + c <- map (dir ) . filter (`notElem` ["..", "."]) <$> getDirectoryContents dir + (,) <$> filterM doesDirectoryExist c <*> filterM doesFileExist c diff --git a/tests/properties.hs b/tests/properties.hs new file mode 100644 --- /dev/null +++ b/tests/properties.hs @@ -0,0 +1,50 @@ +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE PatternGuards #-} +module Main where + +import Control.Applicative +import Control.Monad (guard) +import Control.Lens +import Data.List (nub) +import Data.Map as M +import Data.Maybe (fromMaybe) +import Data.Vector as B +import Data.Vector.Unboxed as U +import Data.Vector.Hybrid as H +import Data.Vector.Generic as G +import Data.Word +import Instances +import Sparse.Matrix as SM +import Test.Framework.TH +import Test.Framework.Providers.QuickCheck2 +import Test.QuickCheck +import Test.QuickCheck.Function +import Linear + +-- model for matrix multiplication +type Linear a = Map Word (Map Word a) + +nonEmpty :: Lens' (Maybe (Map k Int)) (Map k Int) +nonEmpty f m = f (fromMaybe M.empty m) <&> \ m -> m <$ guard (not (M.null m)) + +-- | matrix multiplication in linear will leave empty maps inside the outer map in sparse multiplication +sane :: Linear Int -> Linear Int +sane = M.filter (not . M.null) + +toLinear :: Mat U.Vector Int -> Linear Int +toLinear = sane . H.foldr (\(k,v) r -> r & at (k^._1) . nonEmpty . at (k^._2) ?~ v) M.empty . view _Mat + +fromLinear :: Linear Int -> Mat U.Vector Int +fromLinear m = SM.fromList $ do + (i, n) <- M.toList m + (j, a) <- M.toList n + return (Key i j, a) + +prop_to_from x = toLinear (fromLinear x) == sane x +prop_from_to x = fromLinear (toLinear x) == x + +prop_model :: Mat U.Vector Int -> Mat U.Vector Int -> Gen Prop +prop_model x y | z <- x * y, z' <- fromLinear (toLinear x !*! toLinear y) + = label (show z Prelude.++ " == " Prelude.++ show z') (z == z') + +main = $defaultMainGenerator