packages feed

sparse 0.5 → 0.6

raw patch · 16 files changed

+1449/−268 lines, 16 filesdep +QuickCheckdep +arraydep +bytestringdep −monadplusdep −natsdep −pointeddep ~basedep ~semigroupsbuild-type:Customsetup-changednew-uploader

Dependencies added: QuickCheck, array, bytestring, containers, contravariant, criterion, deepseq, directory, doctest, filepath, hybrid-vectors, lens, linear, mtl, primitive, simple-reflect, sparse, test-framework, test-framework-quickcheck2, test-framework-th, transformers, vector, vector-algorithms

Dependencies removed: monadplus, nats, pointed

Dependency ranges changed: base, semigroups

Files

+ .ghci view
@@ -0,0 +1,1 @@+:set -isrc -idist/build/autogen -optP-include -optPdist/build/autogen/cabal_macros.h
+ .gitignore view
@@ -0,0 +1,16 @@+dist+docs+wiki+TAGS+tags+wip+stats+.DS_Store+.*.swp+.*.swo+*.o+*.hi+*~+*#+.cabal-sandbox+cabal.sandbox.config
+ .travis.yml view
@@ -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}"
+ .vim.custom view
@@ -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 <F2> :exec ":!hasktags -x -c --ignore src"<CR><CR>++" 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"
− COPYING 
@@ -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 <organization> 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 <COPYRIGHT HOLDER> 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.-
+ LICENSE view
@@ -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.
Setup.lhs view
@@ -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}
+ benchmarks/mm.hs view
@@ -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
sparse.cabal view
@@ -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 <ekmett@gmail.com>+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 <hans@hanshoglund.se>-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. <https://www.fpcomplete.com/user/edwardk/revisiting-matrix-multiplication-part-1>+  .+  2. <https://www.fpcomplete.com/user/edwardk/revisiting-matrix-multiplication-part-2>+  .+  3. <https://www.fpcomplete.com/user/edwardk/revisiting-matrix-multiplication-part-3>+  .+  4. <https://www.fpcomplete.com/user/edwardk/revisiting-matrix-multiplication-part-4>+  .+  5. <https://www.fpcomplete.com/user/edwardk/revisiting-matrix-multiplication-part-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++
− src/Data/Sparse.hs 
@@ -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
+ src/Sparse/Matrix.hs view
@@ -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 <ekmett@gmail.com>+-- 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 #-}
+ src/Sparse/Matrix/Internal/Fusion.hs view
@@ -0,0 +1,122 @@+{-# LANGUAGE BangPatterns #-}+-----------------------------------------------------------------------------+-- |+-- Copyright   :  (C) 2013 Edward Kmett+-- License     :  BSD-style (see the file LICENSE)+-- Maintainer  :  Edward Kmett <ekmett@gmail.com>+-- 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 #-}
+ src/Sparse/Matrix/Internal/Heap.hs view
@@ -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 <ekmett@gmail.com>+-- Stability   :  experimental+-- Portability :  non-portable+--+-- Bootstrapped catenable non-empty pairing heaps as described in+--+-- <https://www.fpcomplete.com/user/edwardk/revisiting-matrix-multiplication-part-5>+-----------------------------------------------------------------------------+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 #-}
+ src/Sparse/Matrix/Internal/Key.hs view
@@ -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 <ekmett@gmail.com>+-- 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 <http://en.wikipedia.org/wiki/Z-order_curve> for more information+-- about Morton order.+--+-- How to perform the comparison without interleaving is described in+--+-- <https://www.fpcomplete.com/user/edwardk/revisiting-matrix-multiplication-part-2>+--+----------------------------------------------------------------------------+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 #-}
+ tests/doctests.hsc view
@@ -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 <ekmett@gmail.com>+-- 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
+ tests/properties.hs view
@@ -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