numerical (empty) → 0.0.0.0
raw patch · 29 files changed
+5898/−0 lines, 29 filesdep +HUnitdep +basedep +ghc-primsetup-changed
Dependencies added: HUnit, base, ghc-prim, hspec, hspec-expectations, monad-ste, numerical, primitive, tagged, transformers, transformers-compat, vector, vector-algorithms
Files
- CHANGELOG.md +1/−0
- LICENSE +30/−0
- README.md +94/−0
- Setup.hs +2/−0
- numerical.cabal +156/−0
- src/Control/NumericalApplicative/Backwards.hs +58/−0
- src/Control/NumericalMonad/State/Strict.hs +246/−0
- src/Numerical/Array.hs +106/−0
- src/Numerical/Array/Address.hs +119/−0
- src/Numerical/Array/Layout.hs +18/−0
- src/Numerical/Array/Layout/Base.hs +554/−0
- src/Numerical/Array/Layout/Builder.hs +310/−0
- src/Numerical/Array/Layout/Dense.hs +914/−0
- src/Numerical/Array/Layout/Sparse.hs +913/−0
- src/Numerical/Array/Locality.hs +81/−0
- src/Numerical/Array/Mutable.hs +411/−0
- src/Numerical/Array/Pure.hs +170/−0
- src/Numerical/Array/Range.hs +74/−0
- src/Numerical/Array/Shape.hs +671/−0
- src/Numerical/Array/Storage.hs +267/−0
- src/Numerical/Data/Vector/HPair.hs +235/−0
- src/Numerical/Data/Vector/Pair.hs +205/−0
- src/Numerical/InternalUtils.hs +14/−0
- src/Numerical/Matrix/Basic.hs +32/−0
- src/Numerical/Nat.hs +107/−0
- src/Numerical/World.hs +18/−0
- tests/Main.hs +32/−0
- tests/NumericalUnit/Layout.hs +7/−0
- tests/NumericalUnit/Shape.hs +53/−0
+ CHANGELOG.md view
@@ -0,0 +1,1 @@+version 0.1.0.0
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2013, Carter Tazio Schonwald and Wellposed Limited++++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 Carter Tazio Schonwald nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ README.md view
@@ -0,0 +1,94 @@+[](http://www.wellposed.com)™++# Currently in Pre alpha release engineering, so nearly ready for experimental consumption+(though please dont publicize yet)++# About Numerical-Core+This is the core Package for Numerical Haskell, a project by Carter Schonwald aka+Wellposed Ltd, and (soon I hope!) other contributors.++Numerical-Core is an open source component of the [Wellposed](http://www.wellposed.com)® Numerical Haskell software suite.++##Build Status++[](http://travis-ci.org/wellposed/numerical)+++#note++this library is **pre alpha release** so not all examples / codes may work as expected.+That said, the current api should be enough to prototype and typecheck algorithms.++++# Numerical Haskell+Numerical Haskell is an effort to bring great numerical computation and data analysis+tooling to haskell, and be the best possible platform for sophisticated efforts in those same domains++## What array Formats are Supported++The initial++## What convention is used for indexing?++When you have an index tuple, just think ``x,y,z`` to keep track of the meaning.+Indexing tuples are written as statically sized lists, eg ``x:*y:*z:*Nil``.+This follows the tradition of x,y,z axes used in plotting. Note well: the underlying memory+order can be row OR column major or other!+++All the computations on these static sized lists get specialized away into+nonrecursive computations at their use sites. So in this special scenario, lists aren't a problem!++++# Contributing+Great! Theres so many awesome ways you could help out. Look at CONTRIBUTING.md for more details.+Right now theres a lot of low hanging fruit in improving test coverage,+and soon there'll be many opportunities on the performance tuning and numerical+algorithms/tooling areas.++## bug reports+see bug.md for how to file a bug report+++# Performance FAQ+1. How do I use Numerical haskell to write fast code thats outstandingly high level !?+ * The leading cause of poor performance in numerical routines (aside from poor choice+ in algorithms) is bad memory locality,+ which has but a single easy cure: ** block recursive algorithms **+ * Yes, you heard me, in compiled languages recursion is pretty cheap outside of the inner+ most loops! It also is a fantastic tool for facilitating good memory locality!+ * I'm totally serious, try out the benchmarks for the various versions of the same routines we+ provide!+2. But, what about fusion?+ * Because of certain aspects of the numerical haskell design, we can't *automagically* use+ the fusion optimization facilities of the underlying array representations such as Vector.++# Community+Many member of the Numerical Haskell community can be found on `#numerical-haskell` on freenode IRC+There is also the [numericalhaskell mailing list](https://groups.google.com/forum/#!forum/numericalhaskell)++# Support+The community provides some basic support through the IRC channel, Mailing list,+and the relevant project [issue trackers](http://github.com/wellposed).++If your support needs can't be resolved though those channels, please do not+hesistate to contact Wellposed (aka Carter) to find out more about our support and+professional services options.++++++++++++++++
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ numerical.cabal view
@@ -0,0 +1,156 @@+cabal-version: 2.2+-- Initial numerics-types.cabal generated by cabal init. For further+-- documentation, see http://haskell.org/cabal/users-guide/++-- The name of the package.+name: numerical++-- The package version. See the Haskell package versioning policy (PVP)+-- for standards guiding when and how versions should be incremented.+-- http://www.haskell.org/haskellwiki/Package_versioning_policy+-- PVP summary: +-+------- breaking API changes+-- | | +----- non-breaking API additions+-- | | | +--- code changes with no API change+version: 0.0.0.0++-- A short (one-line) description of the package.+synopsis: core package for Numerical Haskell project++-- A longer description of the package.+description: the core package for Numerical Haskell. Still immature and incomplete++-- The license under which the package is released.+license: BSD-2-Clause++-- The file containing the license text.+license-file: LICENSE++-- The package author(s).+author: Carter Tazio Schonwald++-- An email address to which users can send suggestions, bug reports, and+-- patches.+maintainer: carter at wellposed dot com++-- A copyright notice.+ copyright: Carter Schonwald++category: Math++build-type: Simple+tested-with: GHC == 8.2.2 , GHC == 8.4.3, GHC == 8.6.2+++++-- Constraint on the version of Cabal needed to build this package.++++extra-source-files:+ README.md+ CHANGELOG.md+++source-repository head+ type: git+ location: http://github.com/wellposed/numerical.git+++library+ -- Modules exported by the library.+ exposed-modules:+ Numerical.Matrix.Basic+ Numerical.Array.Layout+ Numerical.Array.Layout.Base+ Numerical.Array.Layout.Dense+ Numerical.Array.Layout.Sparse+ --Numerical.Array.Layout.Dense.Builder+ Numerical.Array.Layout.Builder+ Numerical.Array+ Numerical.World+ Numerical.InternalUtils+ Control.NumericalMonad.State.Strict+ Control.NumericalApplicative.Backwards+ --Numerical.Array.Operations+ Numerical.Array.Shape+ Numerical.Array.Storage+ Numerical.Nat+ Numerical.Array.Mutable+ Numerical.Array.Pure+ Numerical.Array.Range+ Numerical.Array.Locality+ Numerical.Array.Address+ Numerical.Data.Vector.Pair+ Numerical.Data.Vector.HPair++ ghc-options: -Wall+ -- -ddump-simpl -ddump-to-file+ -- --ghc-option=-ddump-simpl --ghc-option=-ddump-to-file --ghc-options=-dsuppress-all+ -- -O2+ -- Modules included in this library but not exported.+ -- other-modules:++ if impl(ghc >= 8.0) && impl(ghc < 8.2)+ ghc-options: -Wno-redundant-constraints+++ -- LANGUAGE extensions used by modules in this package.+ other-extensions:+ PolyKinds+ BangPatterns+ DataKinds+ TypeFamilies+ DeriveDataTypeable+ TypeOperators+ FlexibleInstances+ FlexibleContexts+ ScopedTypeVariables+++ -- Other library packages from which modules are imported.+ build-depends: base >= 4.10 && < 5++ ,primitive >= 0.5 && < 0.8+ ,vector >= 0.11 && < 0.13+ ,tagged >= 0.7 && < 0.9+ ,monad-ste >= 0.1 && < 0.2+ ,transformers >= 0.4 && < 0.6+ ,transformers-compat >= 0.4 && < 0.6+ ,ghc-prim >=0.2 && <0.6+ ,vector-algorithms >= 0.6.0.1 && < 0.9+ -- ,pqueue >= 1.2 && < 1.3+ -- , quickcheck >=++ -- what version constraints?++++ -- Directories containing source files.+ hs-source-dirs: src++ -- Base language which the package is written in.+ default-language: Haskell2010++Test-suite testsuite+ default-language: Haskell2010+ type: exitcode-stdio-1.0+ build-depends: base+ ,hspec >=2.2 && <2.5+ ,hspec-expectations+ ,HUnit >= 1.2.5 && < 1.7+ ,primitive+ ,vector+ ,tagged+ ,transformers+ ,ghc-prim+ ,numerical+ ghc-options: -threaded+ hs-source-dirs: tests+ main-is: Main.hs+ other-modules:+ -- Only modules which are part of the test suite and not the library+ -- should be listed here. These modules are excluded from the coverage+ -- report because their coverage will be 100%.+ NumericalUnit.Shape+ NumericalUnit.Layout
+ src/Control/NumericalApplicative/Backwards.hs view
@@ -0,0 +1,58 @@+++module Control.NumericalApplicative.Backwards where++import Prelude hiding (foldr, foldr1, foldl, foldl1)+import qualified Control.Applicative as A+import Data.Foldable as F+import Data.Traversable as T++-- | The same functor, but with an 'Applicative' instance that performs+-- actions in the reverse order.+newtype Backwards f a = Backwards { forwards :: f a }++-- | Derived instance.+instance (Functor f) => Functor (Backwards f) where+ fmap f (Backwards a) = Backwards (fmap f a)+ {-# INLINE fmap #-}++-- | Apply @f@-actions in the reverse order.+instance (A.Applicative f) => A.Applicative (Backwards f) where+ pure a = Backwards (A.pure a)+ {-# INLINE pure #-}+ Backwards f <*> Backwards a = Backwards (a <**> f)+ {-# INLINE (<*>) #-}+++-- | Try alternatives in the same order as @f@.+instance (A.Alternative f) => A.Alternative (Backwards f) where+ empty = Backwards A.empty+ Backwards x <|> Backwards y = Backwards (x A.<|> y)++-- | Derived instance.+instance (Foldable f) => Foldable (Backwards f) where+ foldMap f (Backwards t) = foldMap f t+ foldr f z (Backwards t) = foldr f z t+ foldl f z (Backwards t) = foldl f z t+ foldr1 f (Backwards t) = foldl1 f t+ foldl1 f (Backwards t) = foldr1 f t++-- | Derived instance.+instance (Traversable f) => Traversable (Backwards f) where+ traverse f (Backwards t) = fmap Backwards (traverse f t)+ sequenceA (Backwards t) = fmap Backwards (sequenceA t)+ mapM f = A.unwrapMonad . T.traverse (A.WrapMonad . f)+ sequence = T.mapM id+ {-#INLINE traverse #-}+ {-#INLINE sequenceA #-}+ {-#INLINE mapM #-}+ {-#INLINE sequence #-}++(<**>) :: A.Applicative f => f a -> f (a -> b) -> f b+(<**>) = liftA2 (flip ($))+{-# INLINE (<**>) #-}++liftA2 :: A.Applicative f => (a -> b -> c) -> f a -> f b -> f c+liftA2 f a b = f `fmap` a A.<*> b+{-# INLINE liftA2 #-}+
+ src/Control/NumericalMonad/State/Strict.hs view
@@ -0,0 +1,246 @@++module Control.NumericalMonad.State.Strict where++--import Data.Functor.Identity+import Control.Monad.Trans.Class+import Control.Monad.IO.Class+import Control.Applicative+import Control.Monad+import Control.Monad.Fix+++{-++This module is a private copy of the Strict State Monad by Ross Patterson,+patched to unconditionally inline.++Its only purpose is to ensure that certain generic routines in+Numerical.Array.Shape will compositionally unconditionally inline in their use sites++ONLY use if writing generic code in your inner loops++-}+++import Data.Foldable (Foldable(foldMap))+import Data.Traversable (Traversable(traverse))++-- | Identity functor and monad.+newtype Identity a = Identity { runIdentity :: a }++-- ---------------------------------------------------------------------------+-- Identity instances for Functor and Monad++instance Functor Identity where+ fmap f m = Identity (f (runIdentity m))+ {-# INLINE fmap #-}++instance Foldable Identity where+ foldMap f (Identity x) = f x+ {-# INLINE foldMap #-}++instance Traversable Identity where+ traverse f (Identity x) = Identity <$> f x+ {-# INLINE traverse #-}++instance Applicative Identity where+ pure a = Identity a+ {-# INLINE pure #-}+ Identity f <*> Identity x = Identity (f x)+ {-# INLINE (<*>) #-}++instance Monad Identity where+ return a = Identity a+ {-# INLINE return #-}+ m >>= k = k (runIdentity m)+ {-# INLINE (>>=)#-}+instance MonadFix Identity where+ mfix f = Identity (fix (runIdentity . f))+ {-# INLINE mfix #-}+++-- ---------------------------------------------------------------------------+-- | A state monad parameterized by the type @s@ of the state to carry.+--+-- The 'return' function leaves the state unchanged, while @>>=@ uses+-- the final state of the first computation as the initial state of+-- the second.+type State s = StateT s Identity++-- | Construct a state monad computation from a function.+-- (The inverse of 'runState'.)+state :: Monad m+ => (s -> (a, s)) -- ^pure state transformer+ -> StateT s m a -- ^equivalent state-passing computation+state = \f -> StateT (return . f)+{-# INLINE state #-}++-- | Unwrap a state monad computation as a function.+-- (The inverse of 'state'.)+runState :: State s a -- ^state-passing computation to execute+ -> s -- ^initial state+ -> (a, s) -- ^return value and final state+runState = \ m -> runIdentity . runStateT m+{-# INLINE runState#-}+++-- | Evaluate a state computation with the given initial state+-- and return the final value, discarding the final state.+--+-- * @'evalState' m s = 'fst' ('runState' m s)@+evalState :: State s a -- ^state-passing computation to execute+ -> s -- ^initial value+ -> a -- ^return value of the state computation+evalState = \m s -> fst (runState m s)+{-# INLINE evalState #-}++-- | Evaluate a state computation with the given initial state+-- and return the final state, discarding the final value.+--+-- * @'execState' m s = 'snd' ('runState' m s)@+execState :: State s a -- ^state-passing computation to execute+ -> s -- ^initial value+ -> s -- ^final state+execState = \m s -> snd (runState m s)+{-# INLINE execState#-}++-- | Map both the return value and final state of a computation using+-- the given function.+--+-- * @'runState' ('mapState' f m) = f . 'runState' m@+mapState :: ((a, s) -> (b, s)) -> State s a -> State s b+mapState = \ f -> mapStateT (Identity . f . runIdentity)+{-# INLINE mapState #-}++-- | @'withState' f m@ executes action @m@ on a state modified by+-- applying @f@.+--+-- * @'withState' f m = 'modify' f >> m@+withState :: (s -> s) -> State s a -> State s a+withState = \f st -> withStateT f st+{-# INLINE withState #-}+-- ---------------------------------------------------------------------------+-- | A state transformer monad parameterized by:+--+-- * @s@ - The state.+--+-- * @m@ - The inner monad.+--+-- The 'return' function leaves the state unchanged, while @>>=@ uses+-- the final state of the first computation as the initial state of+-- the second.+newtype StateT s m a = StateT { runStateT :: s -> m (a,s) }++-- | Evaluate a state computation with the given initial state+-- and return the final value, discarding the final state.+--+-- * @'evalStateT' m s = 'liftM' 'fst' ('runStateT' m s)@+evalStateT :: (Monad m) => StateT s m a -> s -> m a+evalStateT = \ m s -> do+ (a, _) <- runStateT m s+ return a+{-# INLINE evalStateT #-}++-- | Evaluate a state computation with the given initial state+-- and return the final state, discarding the final value.+--+-- * @'execStateT' m s = 'liftM' 'snd' ('runStateT' m s)@+execStateT :: (Monad m) => StateT s m a -> s -> m s+execStateT = \ m s -> do+ (_, s') <- runStateT m s+ return s'+{-# INLINE execStateT #-}+++-- | Map both the return value and final state of a computation using+-- the given function.+--+-- * @'runStateT' ('mapStateT' f m) = f . 'runStateT' m@+mapStateT :: (m (a, s) -> n (b, s)) -> StateT s m a -> StateT s n b+mapStateT = \ f m -> StateT $ f . runStateT m++-- | @'withStateT' f m@ executes action @m@ on a state modified by+-- applying @f@.+--+-- * @'withStateT' f m = 'modify' f >> m@+withStateT :: (s -> s) -> StateT s m a -> StateT s m a+withStateT = \ f m -> StateT $ runStateT m . f++instance (Functor m) => Functor (StateT s m) where+ fmap = \ f m -> StateT $ \ s ->+ fmap (\ (a, s') -> (f a, s')) $ runStateT m s+ {-# INLINE fmap #-}++instance (Functor m, Monad m) => Applicative (StateT s m) where+ pure = \ a ->return a+ (<*>) = \ a b -> ap a b++instance (Functor m, MonadPlus m) => Alternative (StateT s m) where+ empty = mzero+ {-# INLINE empty #-}+ (<|>) = \ a b -> mplus a b+ {-#INLINE (<|>)#-}++instance (Monad m) => Monad (StateT s m) where+ {-# INLINE return #-}+ return = \ a -> state $ \s -> (a, s)+ {-# INLINE (>>=)#-}+ (>>=) = \m k -> StateT $ \s -> do+ (a, s') <- runStateT m s+ runStateT (k a) s'+ fail str = StateT $ \_ -> fail str++instance (MonadPlus m) => MonadPlus (StateT s m) where+ mzero = StateT $ \_ -> mzero+ {-# INLINE mzero #-}+ mplus = \ m n -> StateT $ \s -> runStateT m s `mplus` runStateT n s+ {-# INLINE mplus #-}+instance (MonadFix m) => MonadFix (StateT s m) where+ mfix = \ f -> StateT $ \s -> mfix $ \ ~(a, _) -> runStateT (f a) s+ {-# INLINE mfix #-}++instance MonadTrans (StateT s) where+ {-#INLINE lift #-}+ lift = \ m -> StateT $ \s -> do+ a <- m+ return (a, s)++instance (MonadIO m) => MonadIO (StateT s m) where+ liftIO = lift . liftIO++-- | Fetch the current value of the state within the monad.+get :: (Monad m) => StateT s m s+get = state $ \s -> (s, s)+{-# INLINE get #-}++-- | @'put' s@ sets the state within the monad to @s@.+put :: (Monad m) => s -> StateT s m ()+put = \s -> state $ \_ -> ((), s)+{-# INLINE put #-}++-- | @'modify' f@ is an action that updates the state to the result of+-- applying @f@ to the current state.+--+-- * @'modify' f = 'get' >>= ('put' . f)@+modify :: (Monad m) => (s -> s) -> StateT s m ()+modify = \f -> state $ \s -> ((), f s)+{-# INLINE modify #-}++-- | A variant of 'modify' in which the computation is strict in the+-- new state.+--+-- * @'modify'' f = 'get' >>= (('$!') 'put' . f)@+modify' :: (Monad m) => (s -> s) -> StateT s m ()+modify' f = do+ s <- get+ put $! f s+{-# INLINE modify' #-}++-- | Get a specific component of the state, using a projection function+-- supplied.+--+-- * @'gets' f = 'liftM' f 'get'@+gets :: (Monad m) => (s -> a) -> StateT s m a+gets = \ f -> state $ \s -> (f s, s)+{-# INLINE gets #-}+
+ src/Numerical/Array.hs view
@@ -0,0 +1,106 @@+{-# LANGUAGE PolyKinds #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE TypeFamilies #-}++module Numerical.Array where++--import Numerical.Array.Shape++--data MArray world rep lay (view:: Locality) sh elem+++++++{-+lets do just IO and not ST for now?+or bite the primstate bullet now?+-}++--data family Array world rep lay (view:: Locality) sh elm+++{-+only Row and Column Major have dense formats that are unique across ALL+possible ranks, not so simple for others. Make different data instances per formats++++++-}+--class Array where++-----------+-- | for now locality is a closed type, may change going forward+-- also should it be in a different module like shape or layout?++++{-+theres several points in the design space of array apis that are nice, but none quite right++Vector is probably the closest+pros:+ which has nice pure vs mutable apis+ simple interface+cons:+ its really designed for Int indexing+ assumes every pure Vector is internally derived from an imperative one+ (this is reflected in where the thaw/freeze)++so there needs to be an Array class, an MArray class,++and the Thawing / Freezing needs to be in a seperate PhasedArray class!+why? Because we can't assume that pure/mutable arrays are the fundamental data type!++-}+++{-+maybe do+ data Locality = Contiguous | Strided++For now lets assume that the concrete (rather than delayed) arrays+have a regular structure when strided. (rather than nonuniform gaps)+-}+{-+rep = storable, unboxed, boxed, delay, etc++lay = row major, column major, morton z, morton w (flipped n),+ --- this ignores symmetry and hermitian being properties as well as packed layouts+ --- also need to have a good sparse story+ --- as currently done, most don't really make sense for != rank-2 arrays,++-- rowMajor is a foldR, columnMajor is a foldL over the shape ices++-- Repa and accelerate use a Snoc List so that Row major fuses well for row major++sh= rank / shape, ie matrix or vector, or some higher tensor thingy+lets borrow from repa/ accelea+++mode= need to have a notion of runnable worlds,+based on "backend" chosen, eg CBlasish, DPH, Repa, LLVM, Free (get the shallow/ deep ast)++view =+ Origin, Slice, and Diced, I might make this a fixed universe for now+ Lets not distinguish whether a contiguous array is the original or derived for now+ doesn't seem to be a meaningful difference and would make type inference crap / not bijective+ Note that this does mean that accidental space leaks may happen++ that suggests (but not for now) having a notion of origin / derived+ that would allow elimiting space leaks. But lets not do that for now++-}++{-++uncheckedReshape :: Array wld rep lay+++-}++
+ src/Numerical/Array/Address.hs view
@@ -0,0 +1,119 @@+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE MultiParamTypeClasses #-}+module Numerical.Array.Address(+ Address(..)+ ,SparseAddress(..)+ ) where++import Data.Data+import Control.Monad (liftM)+import qualified Foreign.Storable as Store+import qualified Data.Vector.Unboxed as UV+import qualified Data.Vector.Generic as GV+import qualified Data.Vector.Generic.Mutable as GMV+import GHC.Generics++-- | 'Address' is the type used for addressing into the underlying memory buffers+-- of numerical arrays, Used for Dense Rank n arrays, and 1dim sparse arrays.+newtype Address = Address Int+ deriving (Eq,Ord,Show,Read,Typeable,Generic,Data,Store.Storable)++-- | 'LogicalAddress' is+-- possibly dead code+newtype LogicalAddress = LogicalAddress Int+ deriving (Eq,Ord,Show,Read,Typeable,Generic,Data,Store.Storable)+-- todo, add unboxed for+++-- | this m+--newtype LogicalExtent+-- sparse address seems to be dead atm+data SparseAddress = SparseAddress {+ outerIndex :: {-# UNPACK #-} !Int+ ,innerIndex :: {-# UNPACK #-} !Int }+ deriving (Eq,Show,Data,Generic,Typeable)++{-+At some point decouple logical and physical address+Logical Address should always be Int64 -- maybe even MORE?!+physical address should be native IntPtr (aka Int)++-}++++instance Num Address where+ {-# INLINE (+) #-}+ (+) (Address a) (Address b) = Address (a+b)+ {-# INLINE (-) #-}+ (-) (Address a) (Address b) = Address (a-b)++ (*) _ _ = error "you cant multiply Addresses"+ negate _ = error "you cant Apply Negate to An Address"+ signum _ = error "error you cant take signum of an Address"+ abs _ = error "error you cant take abs of an Address"+ fromInteger _ = error "you cant use Integer Literals or fromInteger to form an Address"++{-+note that+-}++{-+note that i don't think these vector instances ever matter+-}++newtype instance UV.MVector s Address = MV_Address (UV.MVector s Int)+newtype instance UV.Vector Address = V_Address (UV.Vector Int)++instance UV.Unbox Address where++++instance GMV.MVector UV.MVector Address where+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicOverlaps #-}+ {-# INLINE basicUnsafeNew #-}+ {-# INLINE basicUnsafeReplicate #-}+ {-# INLINE basicUnsafeRead #-}+ {-# INLINE basicUnsafeWrite #-}+ {-# INLINE basicClear #-}+ {-# INLINE basicSet #-}+ {-# INLINE basicUnsafeCopy #-}+ {-# INLINE basicUnsafeGrow #-}+ {-# INLINE basicInitialize #-}+ basicInitialize = \ (MV_Address mva) -> GMV.basicInitialize mva+ basicLength (MV_Address v) = GMV.basicLength v+ basicUnsafeSlice i n (MV_Address v) = MV_Address $ GMV.basicUnsafeSlice i n v+ basicOverlaps (MV_Address v1) (MV_Address v2) = GMV.basicOverlaps v1 v2+ basicUnsafeNew n = MV_Address `liftM` GMV.basicUnsafeNew n+ basicUnsafeReplicate n (Address a) = MV_Address `liftM` GMV.basicUnsafeReplicate n a+ basicUnsafeRead (MV_Address v) i = Address `liftM` GMV.basicUnsafeRead v i+ basicUnsafeWrite (MV_Address v) i (Address a) = GMV.basicUnsafeWrite v i a+ basicClear (MV_Address v) = GMV.basicClear v+ basicSet (MV_Address v) (Address a) = GMV.basicSet v a+ basicUnsafeCopy (MV_Address v1) (MV_Address v2) = GMV.basicUnsafeCopy v1 v2+ basicUnsafeMove (MV_Address v1) (MV_Address v2) = GMV.basicUnsafeMove v1 v2+ basicUnsafeGrow (MV_Address v) n = MV_Address `liftM` GMV.basicUnsafeGrow v n++instance GV.Vector UV.Vector Address where+ {-# INLINE basicUnsafeFreeze #-}+ {-# INLINE basicUnsafeThaw #-}+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicUnsafeIndexM #-}+ {-# INLINE elemseq #-}+ basicUnsafeFreeze (MV_Address v) = V_Address `liftM` GV.basicUnsafeFreeze v+ basicUnsafeThaw (V_Address v) = MV_Address`liftM` GV.basicUnsafeThaw v+ basicLength (V_Address v) = GV.basicLength v+ basicUnsafeSlice i n (V_Address v) = V_Address $ GV.basicUnsafeSlice i n v+ basicUnsafeIndexM (V_Address v) i+ = Address `liftM` GV.basicUnsafeIndexM v i+ basicUnsafeCopy (MV_Address mv) (V_Address v)+ = GV.basicUnsafeCopy mv v+ elemseq _ (Address a) z = GV.elemseq (undefined :: UV.Vector a) a z+
+ src/Numerical/Array/Layout.hs view
@@ -0,0 +1,18 @@+++module Numerical.Array.Layout(+ module Numerical.Array.Layout.Base+ ,module Numerical.Array.Layout.Dense+ ,module Numerical.Array.Layout.Sparse+ ,module Numerical.Array.Address+ ) where+++import Numerical.Array.Layout.Base+import Numerical.Array.Layout.Dense+import Numerical.Array.Layout.Sparse+import Numerical.Array.Address++++
+ src/Numerical/Array/Layout/Base.hs view
@@ -0,0 +1,554 @@+{- | Comments for this modules+++-}++{-# LANGUAGE DataKinds #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE NoImplicitPrelude #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE CPP #-}++{-# LANGUAGE StandaloneDeriving #-}++{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE UndecidableInstances #-}++#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 707+ {-# LANGUAGE AutoDeriveTypeable #-}+#endif+++module Numerical.Array.Layout.Base(+ Layout(..)+ ,DenseLayout(..)+ ,RectilinearLayout(..)+ ,LayoutAddress+ ,LayoutLogicalFormat+ ,Transposed+ ,FormatStorageRep+ ,RectOrientationForm+ ,RectDownRankForm+ ,InnerContigForm+ ,Format+ ,TaggedShape(..)+ ,GDSlice(..) --- right? right?+ ,SMajorOrientation(..)+ ,MajorOrientation(..)+ ,majorCompareRightToLeft+ ,majorCompareLeftToRight+ ,shapeCompareRightToLeft+ ,shapeCompareLeftToRight+ -- * All the various helper types+ ,module Numerical.Array.Storage+ ,module Numerical.Array.Locality+ ,module Numerical.Array.Shape+ ,module Numerical.Array.Range+ ,module Numerical.Array.Address+) where+++import Data.Data++import Numerical.Nat+import Numerical.Array.Address+import Numerical.Array.Locality+import Numerical.Array.Shape+import Numerical.Array.Storage+import Numerical.Array.Range++--import Data.Typeable+#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ < 709+import qualified Control.Applicative as A+import Prelude hiding (foldr,foldr1,foldl1,foldl,map)+import qualified Data.Foldable as F+#elif __GLASGOW_HASKELL__ >= 709+import qualified Control.Applicative as A+import qualified Data.Foldable as F+#endif++#if MIN_VERSION_base(4,8,0)+import Prelude hiding (foldl)+#endif++{-+NB: may need to add some specialization for low rank indexing,+theres 4 choices:+a) INLINE EVERYTHING+b) rewrite rules that take low rank indexing code into specialized versions thereof+c) wait till ghc 7.8.2 to resolve https://ghc.haskell.org/trac/ghc/ticket/8848+ and use SPECIALIZE+d) benchmark and then decide++for now I choose (a), and defer benchmarking variations till everything works :)+++a related concern is the interplay of inlining and specialization+https://ghc.haskell.org/trac/ghc/ticket/5928++-}++++++-- either we need to break ties, or the ties have been broken+majorCompareLeftToRight :: Ordering -> Ordering -> Ordering+majorCompareLeftToRight EQ new = new+majorCompareLeftToRight a _ = a+++majorCompareRightToLeft :: Ordering -> Ordering -> Ordering+majorCompareRightToLeft new EQ = new+majorCompareRightToLeft _ b = b++{-# INLINE shapeCompareLeftToRight #-}+shapeCompareLeftToRight :: (F.Foldable (Shape r),A.Applicative (Shape r), Ord a)+ => Shape r a -> Shape r a -> Ordering+shapeCompareLeftToRight = \ ls rs -> foldl majorCompareLeftToRight EQ $ map2 compare ls rs++{-# INLINE shapeCompareRightToLeft #-}+shapeCompareRightToLeft :: ((F.Foldable (Shape r)),A.Applicative (Shape r), Ord a)+ => Shape r a -> Shape r a -> Ordering+shapeCompareRightToLeft = \ ls rs -> foldl majorCompareRightToLeft EQ $ map2 compare ls rs+++++-- | this is kinda a hack+newtype TaggedShape (form :: *) (rank::Nat) = TaggedShape {unTagShape:: Shape rank Int }+instance Eq (Shape rank Int)=> Eq (TaggedShape f rank) where+ (==) l r = (==) (unTagShape l) (unTagShape r )++instance Show (Shape rank Int) => Show (TaggedShape f rank) where+ show (TaggedShape ix) = "TaggedShape (" ++ show ix ++ " )"++instance forall form rank . (Eq (Shape rank Int),Layout form rank)+ => Ord (TaggedShape form rank) where+ compare left right = basicCompareIndex (Proxy:: Proxy form ) (unTagShape left) (unTagShape right)+++-- | Generalized Dense Slice Projection notation,+-- not sure if it should be defined in this module or elsewhere+-- This provides a type safe interface for the classical+-- general array slice notation.+-- That said, its only useful for dense array formats,+-- at least in general. For formats that aren't "rectilinear dense",+-- this COULD be used as a description format for traversing+-- over various rectilinear subsets of points though?+data GDSlice (from :: Nat) (to :: Nat) :: * where+ GDNil :: GDSlice 'Z 'Z+ GDPick :: Int -> !(GDSlice from to) -> GDSlice ('S from) to+ GDRange :: (Int,Int,Int) {- this is a nonempty interval or error -} -> !(GDSlice from to) -> GDSlice ('S from) ('S to)+ GDAll :: !(GDSlice from to) -> GDSlice ('S from) ('S to)++{-+TODO: for things that++-}+++instance Show (GDSlice 'Z 'Z) where+ show _ = "GDNil"++instance (Show (GDSlice (f) ('S t)),Show (GDSlice f t))=> Show (GDSlice ('S f) ('S t)) where+ show (tup `GDRange` rest) = show tup ++ " `GDRange` (" ++ show rest ++ ")"+ show (GDAll rest) = "GDAll " ++ show rest+ show (ix `GDPick` rest) = show ix ++" `GDPick` " ++ show rest+++instance Show (GDSlice f 'Z)=> Show (GDSlice ('S f) 'Z) where+ show (ix `GDPick` rest) = show ix ++" `GDPick` " ++ show rest+--instance Show (GDSlice f t) where+-- func =++{-+In some (moderately precise sense)++-}+++-- GDRange (from,step,to)+ -- GDAll is just sugar for a special case of GDRange, but maybe its worthwhile sugar?++--computeSlicePlan:: GDSlice from to -> Shape from Int -> Shape from (Either Int (AffineRange Int))+--computeSlicePlan GDNil Nil = Nil+--computeSlicePlan ( ix `GDPick` gdRest )+-- (bd:* shpRest)| ix < bd && ix >= 0 = Left ix :* computeSlicePlan gdRest shpRest+-- | otherwise = error+-- $ "bad indices for computeSlicePlan " ++ show (ix,bd)+--computeSlicePlan ( (strt,step,end) `GDRange` grest) (bd:* shprest)++++data family Format lay (contiguity:: Locality) (rank :: Nat) rep++deriving instance Typeable Format++type family FormatStorageRep ( a:: * ) :: *++type instance FormatStorageRep (Format lay ctg rnk rep)= rep++type family Transposed (form :: *) :: *++type family LayoutAddress (form :: *) :: *++-- TODO / FIXME remove the basic* prefix from all the operations+-- this was done originally because+++-- TODO : should this be pushed into the type class?+-- TODO : should this be pushed into the type class?+-- | every format has a "logical" sibling, that represents the address translation+-- when the underlying buffer layer is contiguous and packed. So it could be claimed+-- that any type that obeys @a~'LayoutLogicalFormat' a@ is one that an be a legal+-- instance of LayoutBuilder?+type family LayoutLogicalFormat (form :: *) :: *++-- | the 'Layout' type class+class Layout form (rank :: Nat) | form -> rank where++ -- | 'basicLogicalShape' gives the extent of the format+ basicLogicalShape :: form -> Shape rank Int++ -- | 'basicLogicalForm' converts a given format into its "contiguous" analogue+ -- this is useful for supporting various address translation manipulation tricks+ -- efficiently. Note that any valid simple format should strive to ensure this is an O(1) operation.+ -- though certain composite 'Layout' instances may provide a slower implementation.+ basicLogicalForm :: (logicalForm ~ LayoutLogicalFormat form ) => form -> logicalForm+++ -- | 'transposedLayout' transposes the format data type+ transposedLayout :: (form ~ Transposed transform,transform~Transposed form)=> form -> transform++ -- | 'basicCompareIndex' lets you compare where two (presumably inbounds)+ -- 'Index' values are in a formats ordering. The logical 'Shape' of the array+ -- is not needed+ basicCompareIndex :: p form-> Shape rank Int ->Shape rank Int -> Ordering++ -- | the (possibly empty) min and max of the valid addresses for a given format.+ -- @minAddress = fmap _RangeMin . rangedFormatAddress@+ -- and @maxAddress = fmap _RangeMax . rangedFormatAddress@+ -- FIXME : This also is a terrible name+ basicAddressRange :: (address ~ LayoutAddress form)=> form -> Maybe (Range address)+ -- FIX ME! this name is crap, i dont like it++ -- | 'basicToAddress' takes an Index, and tries to translate it to an address if its in bounds+ --+ basicToAddress :: (address ~ LayoutAddress form)=>+ form -> Index rank -> Maybe address++ -- | 'basicToIndex' takes an address, and always successfully translates it to+ -- a valid index. Behavior of invalid addresses constructed by a library user+ -- is unspecified.+ basicToIndex ::(address ~ LayoutAddress form)=>+ form -> address -> Index rank++ -- | 'basicNextAddress' takes an address, and tries to compute the next valid+ -- address, or returns Nothing if there is no subsequent valid address.+ basicNextAddress :: (address ~ LayoutAddress form)=>+ form -> address -> Maybe address++ -- | @'basicNextIndex' form ix mbeAddress@ computes the next valid index after+ -- @ix@ if it exists. It takes a @'Maybe' address@ as a hint for where to do the search for the successor.+ -- If the index is in bounds and not the last index, it returns both the index and the associated address.+ basicNextIndex :: (address ~ LayoutAddress form)=>+ form -> Index rank -> Maybe address -> Maybe ( Index rank, address)+++ basicAddressPopCount :: (address ~ LayoutAddress form)=>+ form -> Range address -> Int++ -- | This operation is REALLY unsafe+ -- This should ONLY be used on Formats that are directly+ -- paired with a Buffer or Mutable Buffer (ie a Vector)+ -- This operation being in this class is also kinda a hack+ -- but lets leave it here for now+ basicAddressAsInt :: (address ~ LayoutAddress form)=>+ form -> address -> Int+ basicAddressAsInt =+ \ _ _ ->+ error "called basicAddressAsInt on a Layout thats not meant for this world"++ -- | The semantics of @`basicAffineAddressShift` form addr step@ is that+ -- when step > 0, its equivalent to iteratively computing 'basicNextAddress' @step@ times.+ -- However, the step size can be negative, which means it can+ basicAffineAddressShift :: (address ~ LayoutAddress form) =>+ form -> address -> Int -> Maybe address+++ {-# MINIMAL basicToAddress, basicToIndex, basicNextAddress,basicNextIndex+ ,basicAddressRange,basicLogicalShape,basicCompareIndex+ , transposedLayout, basicAddressPopCount,basicLogicalForm, basicAffineAddressShift #-}+++{- |+these names aren't ideal, but lets punt on bikeshedding till theres >= 2 serious+users+-}+data MajorOrientation = Rowed | Columned | BlockedColumn | BlockedRow+ deriving(Data,Typeable)++data SMajorOrientation (o :: MajorOrientation) where+ SRowed :: SMajorOrientation 'Rowed+ SColumned :: SMajorOrientation 'Columned+ SBlockedRow :: SMajorOrientation 'BlockedRow+ SBlockedColumn :: SMajorOrientation 'BlockedColumn+++-- | Every instance of 'RectilinearLayout' needs to have a corresponding+-- 'RectOrientationForm', 'RectDownRankForm', and 'InnerContigForm'+type family RectOrientationForm form :: MajorOrientation++type family RectDownRankForm form :: *++type family InnerContigForm form :: *++{- | 'RectilinearLayout' is the type class that supports the modle widely+ usable class of slicing operations in Numerical.+ for every instance @'RectilinearLayout' format rank orientation@, a corresponding+ @'RectOrientationForm' form @, @'RectDownRankForm' form@+ and @'InnerContigForm' form@ type family instance should be defined++ The purpose of 'RectilinearLayout' class is to provide++-}+class Layout form rank =>+ RectilinearLayout form (rank :: Nat) (oriented :: MajorOrientation) | form -> rank oriented where++ -- | 'formRectOrientation' provides a runtime mechanism for reflecting+ -- the orientation of the format+ formRectOrientation :: p form -> SMajorOrientation oriented++{-+is array layout always static?+for now lets say yes, cause you can always just existential up the class++-}++ -- | For @'rectlinearShape' form==shp@, we always have that+ -- @'basicLogicalShape' form `weaklyDominates` shp@.+ -- when 'strictlyDominates' holds, that implies that the underlying array format+ -- is a rectilinear layout whose "elements" are tiles of a fixed size array format.+ -- For this initial release and initial set of applicable rectilinear array formats,+ -- the following is always true @'basicLogicalShape' form == basicLogicalShape' form @+ -- Should be @O(1)@ always. Or more precisely @O(rank)@+ rectlinearShape :: form -> Index rank++ unconsOuter:: ('S down ~ rank)=> p form -> Shape rank a -> (a, Shape down a)+ consOuter :: ('S down ~ rank)=> p form -> a -> Shape down a -> Shape rank a++ -- | @'majorAxisSlice' fm (x,y)@ requires that y-x>=1, ie that more than+ -- one sub range wrt the major axis be selected, so that the logical+ -- rank of the selected array stays the same. This operation also preserves+ -- memory locality as applicable.+ -- @O(1)@ / @O(rank)@+ majorAxisSlice :: form -> (Int,Int)-> form+ -- should this be -> Maybe form?+++ -- | @'majorAxixProject' form x@ picks a "row" with respect to the outer most+ -- dimension of the array format. This will+ -- @O(1)@ or @O(rank)@+ majorAxisProject :: (RectilinearLayout downForm subRank oriented,+ rank ~ ('S subRank) , downForm~ RectDownRankForm form) => form -> Int -> downForm++ -- | this is the nonstrided subset of general array slice notation.+ -- Invoke as @'rectilinearSlice' form leastCorner greatestCorner@,+ -- where the least and greatest corners of the sub array are determined+ -- by the 'strictlyDominates' partial order on the bounds of the sub array.+ -- For Dense array formats, this should be @O(1)@ or more precisely @O(rank)@+ -- For the basic Sparse array formats thus far the complexity should be+ -- @O(size of outermost dimension)@, which could be computed by+ -- @fst . unconsOuter [form] . rectilinearShape $ form@+ rectlinearSlice :: (RectilinearLayout icForm rank oriented,icForm~InnerContigForm form )=>form -> Index rank -> Index rank -> icForm -- FIXME, need the range infos????? (icfFOrm, adddress,address)+++{- | 'DenseLayout' only has instances for Dense array formats.+this class will need some sprucing up for the beta,+but its ok for now. NB that 'DenseLayout' is really strictly meant to be used+for optimization purposes, and not meant as a default api+-}+class Layout form rank => DenseLayout form (rank :: Nat) | form -> rank where++++ basicToDenseAddress :: form -> Index rank -> Address++ basicToDenseIndex :: form -> Address -> Index rank++++ basicNextDenseAddress :: form -> Address -> Address+ basicNextDenseAddress = \form shp -> snd+ (basicNextDenseIndex form $ basicToDenseIndex form shp )+ {-# INLINE basicNextDenseAddress #-}++ basicNextDenseIndex :: form -> Index rank ->(Index rank ,Address)+ basicNextDenseIndex = \form shp -> (\ addr ->( basicToDenseIndex form addr, addr) ) $!+ basicNextDenseAddress form $ basicToDenseAddress form shp+ {-# INLINE basicNextDenseIndex #-}+++#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 707+ {-# MINIMAL basicToDenseIndex, basicToDenseAddress,+ (basicNextDenseIndex | basicNextDenseAddress) #-}+#endif++{-+*Numerical.Array.Layout> basicToAddress (FormColumn (2 :* 3 :* 7 :* Nil)) (0:* 2 :* 2 :* Nil)+Address 16+*Numerical.Array.Layout> basicToAddress (FormColumn (2 :* 3 :* 7 :* Nil)) (1:* 0 :* 0 :* Nil)+Address 1+*Numerical.Array.Layout> basicToAddress (FormColumn (2 :* 3 :* 7 :* Nil)) (0:* 0 :* 0 :* Nil)+Address 0+*Numerical.Array.Layout> basicToAddress (FormColumn (2 :* 3 :* 7 :* Nil)) (0:* 1 :* 0 :* Nil)+Address 2+*Numerical.Array.Layout> basicToAddress (FormColumn (2 :* 3 :* 7 :* Nil)) (0:* 0 :* 1 :* Nil)++++-}+++--data Elem ls el where+-- Point :: Elem '[] el+-- (:#) :: a -> Elem ls el -> Elem (a ': ls) el+++{-+ One important invariant about all layouts at all ranks is that for+ any given ints x < y, that the array index for inr++ toIndex shapedLayout (pure x :: Shape rank Int) is strictly less than+ toIndex shapedLayout (pure y :: Shape rank Int).++ more generally++ for rank k tuples,+ xi = x_1 :* ... :* x_k *: Nil and+ yj = y_1 :* ... :* x_k *: Nil+ such that forall \ell, x_\ell < y_\ell+ we have that+ toIndex shapedLayout xi < toIndex shapedLayout yj+++this actually relates to the notion of partial ordering over vectors in convex+geometry!+++so roughly: we have layouts that are dense+we have layouts that can be used as tiles (and are dense)++and we have layouts which can can't be tiled, but can have elements which are tiled++So we have++PrimitiveLayouts++Static Layouts++General Layouts (which are a Top level layout over a static layout)++the Layout class tries to abstract over all three cases+(NB: this only makes sense when the "rank" for the inner+and outer layouts have the same rank!)++-}+++{- Sized is used as a sort of hack to make it easy to express+ the staticly sized layouts. NB, one trade off is that its only+ possible to express "cube" shaped blocks, but on the other+ hand blocking sizes are expressible for every single rank!+-}+--data Sized :: * -> * where+ --(:@) :: Nat -> a -> Sized a+++{-++per se I don't need the StaticLay, PrimLay, Lay constructors, BUT+I really do like how it makes things a teeny bit simpler.. though I may remove them+-}++++--class SimpleDenseLayout lay (rank :: Nat) where+-- type SimpleDenseTranpose lay+-- toIndexSimpleDense :: Shaped rank Int lay -> Shape rank Int -> Int+++--class PrimLayout lay (rank :: Nat) where+-- type TranposedPrim lay+-- toIndexPrim :: Shaped rank Int (PrimLay lay) -> Shape rank Int -> Int+-- fromIndexPrim :: Shaped rank Int (PrimLay lay) -> Int -> Shape rank Int+++{-+for now we will not deal with nested formats, but this will+be a breaking change i plan for later+-}++{-+what is the law for the Layout class?+forall valid formms+toIndex sd (fromIndex sd ix)==ix+fromIndex sd (toIndex sd shp)==shp+-}++{-+if tup1 is strictly less than tup2 (pointwise),+ then any lawful Layout will asign tup1 an index strictly less than that+ asigned to tup2++ transposedLayout . transposedLayout == id++++i treat coordinates as being in x:* y :* z :* Nil, which is Fortran style idexing++in row major we'd have for x:* y :* Nil that X is the inner dimension, and y the outter,+by contrast, in column major, y would be the inner most, and x the outter most.+++++-}+++{- In some respects, the Layout type class is a multidimensional+analogue of the Enum type class in Haskell Prelude,+for Dense / Dense Structured matrix formats+but+ a) requires a witness value, the "Form"+ b) needs to handle multivariate structures+ c) has to deal with structure matrices, like triangular, symmetric, etc+ e) I think every layout should have pure 0 be a valid index, at least for "Dense"+ arrays+ f) transposedLayout . transposedLayout == id+ g)++ Form needs to carry the shape / extent of the matrix++-}+{-++-}++--data View = Origin | Slice+{-+i'm really really hoping to not need a View parameter,+but the nature of the addressing logic needs to change when its a slice+vs a deep copy (for certain classes of arrays that I wish to support very easily)++I will be likely adding this the moment benchmarks validate the distinction++on the+-}
+ src/Numerical/Array/Layout/Builder.hs view
@@ -0,0 +1,310 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE ScopedTypeVariables#-}+-- {-# LANGUAGE DeriveFunctor #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE NoImplicitPrelude #-}++module Numerical.Array.Layout.Builder where++import Control.Monad.Primitive ( PrimMonad, PrimState )+import qualified Data.Vector.Generic as VG+import qualified Data.Vector.Generic.Mutable as VGM+import Numerical.Array.Layout.Base+import Numerical.Array.Layout.Dense as Dense+--import Numerical.Array.Layou.Sparse as Sparse+--import Numerical.Data.Vector.Pair+import Control.Monad.ST (runST)+import Data.Typeable+import qualified Data.Foldable as F+import Data.Traversable as T+import Control.Applicative as A++import Numerical.Data.Vector.Pair+import Numerical.Array.Layout.Sparse+import Data.Vector.Algorithms.Intro as IntroSort+import Data.List (group)+import Numerical.InternalUtils+import Prelude hiding (error)+++data BatchInit v = BatchInit { batchInitSize :: !Int+ ,batchInitKV :: !(Either [v] (IntFun v)) }+ deriving (Typeable)++materializeBatchMV :: (PrimMonad m, VGM.MVector mv a) => BatchInit a -> m (mv (PrimState m) a)+materializeBatchMV (BatchInit size (Left ls )) =+ do+ newMV <- VGM.new size+ _ <- Prelude.mapM (\(ix ,val )-> VGM.unsafeWrite newMV ix val ) (zip [0..] $ take size ls)+ return newMV+materializeBatchMV (BatchInit size (Right (IntFun f) )) =+ do+ newMV <- VGM.new size+ _ <- Prelude.mapM (\ix -> do v <- (f ix) ; VGM.unsafeWrite newMV ix v ) $ take size [0..]+ return newMV++--- not sure if this is EVER useful+newtype AnyMV mv e = AMV (forall s . mv s e )+++instance (Show a)=> Show (BatchInit a) where+ show (BatchInit size (Left ls) ) | size > 100 = "(BatchInit " ++show size +++ "-- only showing the first 100 elements\n"+ ++ "(Left "++(show $ take 100 ls ) ++ "))\n"+ | otherwise ="(BatchInit " ++show size +++ " (Left "++(show ls ) ++ "))\n"+ show (BatchInit size (Right (IntFun f)) ) | size > 100 = "(BatchInit " ++show size +++ "-- only showing the first 100 elements\n"+ ++ "(Left "++(show $ runST (Prelude.mapM f [0..100]) ) ++ "))\n"+ | otherwise ="(BatchInit " ++show size+ ++ "(Left "++(show $ runST (Prelude.mapM f [0,1..size -1]) ) ++ "))\n"+++newtype IntFun a = IntFun (forall m. (PrimMonad m)=> Int -> m a )+-- This may change substantially in a future release, but for now+-- acts like+ deriving (Typeable)++instance Functor IntFun where+ fmap f (IntFun g) = IntFun (\x-> g x >>= (\ y -> return (f y)) )+ {-# INLINE fmap #-}++instance Functor BatchInit where+ {-# INLINE fmap #-}+ fmap = \f bival ->+ case bival of+ (BatchInit size (Left ls))->+ BatchInit size (Left (Prelude.map f ls ))+ (BatchInit size (Right gfun))->+ BatchInit size (Right $ fmap f gfun )++++-- batchInit size should be Word rather than Int for size, but Vector is lame like that+++{-+ChoiceT from monad lib is tempting+as is one of the ListT done right+Bundle from Vector 0.11 and Stream from 0.10 are both alluring too++but all of them make things complicated,+punt for now+++ALso: I may want/need to distinguish sparse vs dense builders+and put them into different classes, punting that for now+-}+++fromListBI :: [a] -> BatchInit a+fromListBI ls = BatchInit (length ls) (Left ls)++fromVectorBI :: VG.Vector v e => v e -> BatchInit e+fromVectorBI v = BatchInit size+ (Right+ (IntFun $+ \i -> if i >= size+ then error $ " out of bounds index on IntFun of size: " ++ show i+ else return $ v VG.! i+ ))+ where+ size = VG.length v++fromMVectorBI :: (VGM.MVector mv e ) => AnyMV mv e -> BatchInit e+fromMVectorBI (AMV v) = BatchInit size+ (Right+ (IntFun $+ \i -> if i >= size+ then error $ " out of bounds index on IntFun of size: " ++ show i+ else v `VGM.read` i+ ))+ where+ size = VGM.length v++++class Layout form (rank::Nat) => LayoutBuilder form (rank::Nat) | form -> rank where++ buildFormatM :: (store~FormatStorageRep form,Buffer store Int ,Buffer store a,PrimMonad m)=>+ Index rank -> proxy form -> a+ -> Maybe (BatchInit (Index rank ,a))+ ->m (form, BufferMut store (PrimState m) a )+++buildFormatPure:: forall store form rank proxy m a. (LayoutBuilder form (rank::Nat)+ ,store~FormatStorageRep form,Buffer store Int ,Buffer store a, Monad m ) =>+ Index rank -> proxy form -> a -> Maybe (BatchInit (Index rank ,a))+ ->m (form, BufferPure store a )+buildFormatPure shape prox defaultValue builder =+ do res@(!_,!_)<-return $! theComputation+ return res+ where+ theComputation :: (form,BufferPure store a )+ !theComputation = runST $+ do (form,buf) <- buildFormatM shape prox defaultValue builder+ pureBuff <- VG.unsafeFreeze buf+ return (form, pureBuff)+{-+this is a funky api for both dense and sparse arrays general builder format.++given the target shape, logical dimensions,a default value (only used for dense arrays)+and the list of manifest values (mostly only used for sparse), build the format+descriptor and the suitably initialized and sized values buffer++this api is only meant for internal use for building new array values+++TODO: compare using a catenable priority heap vs just doing fast sorting.+-}+++{-+the dense instances ignore the builder structure, which does suggest that maybe+there shoudl be a dense builder layout class and a sparse layout class separately+-}++instance LayoutBuilder (Format Direct 'Contiguous ('S 'Z) rep) ('S 'Z) where++ buildFormatM (size:* _) _ defaultValue _ =+ do+ buf<- VGM.replicate size defaultValue+ return (FormatDirectContiguous size,buf)+++-- really wish I didn't have to write the foldable and traversable constraints+-- seems like a code smell?!+instance (F.Foldable (Shape r),T.Traversable (Shape r) ,A.Applicative (Shape r))+ => LayoutBuilder (Format Row 'Contiguous r rep) r where++ buildFormatM ix _ defaultValue _ =+ do+ buf<- VGM.replicate (F.foldl' (*) 0 ix) defaultValue+ return (FormatRowContiguous ix,buf)++instance (F.Foldable (Shape r),T.Traversable (Shape r) ,A.Applicative (Shape r))+ => LayoutBuilder (Format Column 'Contiguous r rep) r where++ buildFormatM ix _ defaultValue _ =+ do+ buf<- VGM.replicate (F.foldl' (*) 0 ix) defaultValue+ return (FormatColumnContiguous ix,buf)++isStrictlyMonotonicV ::(VG.Vector v e)=> (e -> e->Ordering)-> v e -> Maybe Int+isStrictlyMonotonicV cmp v = go 0 (VG.length v)+ where+ go !i !len | i+1 >= len = Nothing+ | (v VG.! i) `lt` (v VG.! (i+1))= go (i+1) len+ | otherwise = Just i++ lt a b = case cmp a b of+ LT -> True+ _ -> False+++instance (Buffer rep Int)=>LayoutBuilder (Format DirectSparse 'Contiguous ('S 'Z) rep ) ('S 'Z) where+++ buildFormatM (size:* _) _ _ Nothing = do+ mvI <- VGM.new 0+ vI <- VG.unsafeFreeze mvI+ mvV <- VGM.new 0+ return $! (FormatDirectSparseContiguous size 0 vI, mvV)++ buildFormatM (size:* _) _ _ (Just builder)= do+ -- need to use let so type inference doesnt totally barf+ mvt@(MVPair (MVLeaf ix) (MVLeaf val)) <- materializeBatchMV $ fmap ( \((ix:*_),v)-> (ix,v)) builder+ -- if i swap to using this instead of ix <- mat.. ; val <- mat..+ --i get CRAZY type errors+ -- could this be a ghc bug?+ --ix <- materializeBatchMV $ fmap fst builtTup+ --val <- materializeBatchMV $ fmap snd builtTup+ _<- IntroSort.sortBy (\x y -> compare (fst x) (fst y)) mvt+ -- this lets me sort a pair of arrays!+ vIx <- VG.unsafeFreeze ix+ optFail <- return $ isStrictlyMonotonicV compare vIx+ --_hoelly+ case optFail of+ Nothing -> return (FormatDirectSparseContiguous size 0 vIx, val)+ Just ixWrong -> error $ "DirectSparse Index duplication at index "++ show (vIx VG.! ixWrong)+++instance (Buffer rep Int) => LayoutBuilder (Format CompressedSparseRow 'Contiguous ('S ('S 'Z)) rep ) ('S ('S 'Z)) where++ buildFormatM (x:* y :* _) _ _ Nothing= do+ mvi <- VGM.new 0+ vi <- VG.unsafeFreeze mvi+ mvval <- VGM.new 0+ return $+ (FormatContiguousCompressedSparseRow+ (FormatContiguousCompressedSparseInternal y x vi vi), mvval )++ buildFormatM (x:* y :* _) proxyFormat _ (Just builder) = do+ mvtup@(MVPair (MVPair (MVLeaf mvectYs) (MVLeaf mvectXs)) (MVLeaf mvectVals))<-+ materializeBatchMV $ fmap (\((xix:* yix :* _),val)-> ((yix,xix),val) ) builder+ _ <- IntroSort.sortBy (\((y1,x1),_) ((y2,x2),_) -> basicCompareIndex proxyFormat (x1:*y1 :* Nil) (x2:*y2:* Nil) )+ mvtup+ vectXs <- unsafeBufferFreeze mvectXs+ vectYs <- unsafeBufferFreeze mvectYs+ --- predicate check here wrt monotonicity of+ --- compute runlength partial sums of where ys go++ -- need to actually check+ yRunsVect <- return $+ VG.replicate (y+1) (0::Int) VG.// computeStarts (computeRunLengths vectYs) 0 y+ --_ <- (error "computeRUnCount") vectYs yRunsMVect+ --yRunsVect <- unsafeBufferFreeze yRunsMVect+ let xyVect = (VPair (VLeaf vectXs) (VLeaf vectYs))+ optFail <- return $+ isStrictlyMonotonicV (\(x1,y1) (x2,y2)->basicCompareIndex proxyFormat (x1:*y1:*Nil) (x2:*y2:*Nil))+ xyVect+ case optFail of+ Nothing -> return $+ (FormatContiguousCompressedSparseRow+ (FormatContiguousCompressedSparseInternal y x vectXs yRunsVect), mvectVals )+ Just i ->+ error $ "illegal duplication in CSR builder (x,y) coordinates "+ ++ show (xyVect VG.! i) ++ " and " ++ show (xyVect VG.! (i+1))+ ++ "starting at position " ++ show i+++computeRunLengths :: (VG.Vector v e, Eq e)=> v e -> [(e,Int)]+computeRunLengths = \y -> fmap (\x ->(head x,length x)) $ group $ VG.toList y++++{-# SPECIALIZE INLINE computeStarts :: [(Int,Int)]->Int->Int ->[(Int,Int)] #-}++computeStarts:: (Enum a, Ord a, Num b )=>[(a,b)]-> a -> a -> [(a,b)]+computeStarts [] start end | start <= end = fmap (\x -> (x ,0)) [start..end]+ | otherwise = error "bad start end arguments to computeStarts"+computeStarts ls start end | start <= end = go start 0 ls+ | otherwise = error "bad start end arguments to computeStarts"+ where+ --go :: a ->b->[(a,b)]-> [(a,b)]+ go !posNext preSum [] | posNext <= end = fmap (\x -> (x,preSum)) [posNext .. end]+ | otherwise = error "impossible go computeStarts "+ go !posNext !preSum gls@((posAt,atSum):rest)+ | posNext < posAt= (posNext,preSum): go (succ posNext) preSum gls+ | posNext == posAt = (posNext,preSum) : go (succ posNext) (preSum + atSum) rest+ | otherwise = error "bad position in prefix stream for computeStarts go, literally unpossible "++++--computeStarts :: (Eq a, Num a)=> [(a,Int)]->Int -> [(a,Int)]+--computeStarts [] len = map (\x -> (x ,0)) [0..len]+--computeStarts ls len = go 0 0 ls+-- where+-- go preSum place [] | place > len = []+-- | place == len = [(place,preSum)]+-- | otherwise = map
+ src/Numerical/Array/Layout/Dense.hs view
@@ -0,0 +1,914 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE NoImplicitPrelude #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE StandaloneDeriving #-}+#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 707+ {-# LANGUAGE AutoDeriveTypeable #-}+#endif++module Numerical.Array.Layout.Dense(+ DenseLayout(..)+ ,Locality(..)+ ,Format(..)+ ,Row+ ,Column+ ,Direct+ ,module Numerical.Array.Layout.Base+ ) where++++import Numerical.Nat+import Control.Applicative+import Numerical.Array.Locality+import Numerical.Array.Layout.Base+import Numerical.Array.Shape as S+import Data.Data(Data,Typeable)+++--import Data.Traversable (Traversable)++import Control.NumericalMonad.State.Strict++import qualified Data.Foldable as F+import Data.Traversable++import Prelude hiding (foldr,foldl,map,scanl,scanr,scanl1,scanr1)+++data Direct++data Row+++data Column+++{-+one important gotcha about shape is that for many formats,+the Shape is the (fmap (+1)) of the largestIndex,+often, but perhaps not always.++-}++{-++need to figure out how to support symmetric and hermitian and triangular+and banded matrices++-}+++--class Layout form rank => DenseLayout form (rank :: Nat) | form -> rank where+ {-+ empty class instances for all the dense Layouts+ -}++type instance LayoutLogicalFormat (Format Direct cont ('S 'Z ) rep )+ = Format Direct 'Contiguous ('S 'Z) rep++-- | @'Format' 'Direct' 'Contiguous' ('S' 'Z')@ is a 1dim array 'Layout' with unit stride+data instance Format Direct 'Contiguous ('S 'Z) rep =+ FormatDirectContiguous {+ logicalShapeDirectContiguous :: {-#UNPACK#-} !Int }+ deriving (Show,Eq,Data)++-- | @'Format' 'Direct' 'Strided' ('S' 'Z')@ is a 1dim array 'Layout' with a regular stride >= 1+data instance Format Direct 'Strided ('S 'Z) rep =+ FormatDirectStrided {+ logicalShapeDirectStrided :: {-#UNPACK#-}!Int+ ,logicalStrideDirectStrided:: {-#UNPACK#-}!Int}+ --deriving (Show,Eq,Data)++type instance LayoutLogicalFormat (Format Row cont n rep )+ = Format Row 'Contiguous n rep++-- | @'Format' 'Row' 'Contiguous' n@ is a rank n Array+data instance Format Row 'Contiguous n rep =+ FormatRowContiguous {+ boundsFormRow :: !(Shape n Int)}+ --deriving (Show,Eq,Data)++data instance Format Row 'Strided n rep =+ FormatRowStrided+ {boundsFormRowStrided:: !(Shape n Int)+ ,strideFormRowStrided:: !(Shape n Int)}+ --deriving (Show,Eq,Data)++data instance Format Row 'InnerContiguous n rep =+ FormatRowInnerContiguous {+ boundsFormRowInnerContig :: !(Shape n Int)+ ,strideFormRowInnerContig:: !(Shape n Int)}+ --deriving (Show,Eq,Data)++type instance LayoutLogicalFormat (Format Column cont n rep )+ = Format Column 'Contiguous n rep++data instance Format Column 'Contiguous n rep =+ FormatColumnContiguous {+ boundsColumnContig :: !(Shape n Int)}+ --deriving (Show,Eq,Data)+--deriving instance (Data (Shape n Int),Typeable n,Typeable rep) =>++data instance Format Column 'InnerContiguous n rep =+ FormatColumnInnerContiguous {+ boundsColumnInnerContig :: !(Shape n Int)+ ,strideFormColumnInnerContig:: !(Shape n Int)+ }++deriving instance Show (Shape n Int) => Show (Format Column 'InnerContiguous n rep)+deriving instance (Data (Shape n Int),Typeable n,Typeable rep) =>Data (Format Column 'InnerContiguous n rep)+ --deriving (Show,Eq,Data)++data instance Format Column 'Strided n rep =+ FormatColumnStrided {+ boundsColumnStrided :: !(Shape n Int)+ ,strideFormColumnStrided:: !(Shape n Int)}+deriving instance Show (Shape n Int) => Show (Format Column 'Strided n rep)+--deriving instance (Eq (Shape n Int)) => Eq (Format Column Strided n rep)+deriving instance (Data (Shape n Int),Typeable n,Typeable rep) => Data (Format Column 'Strided n rep)+ --deriving (Show,Eq,Data)+++type instance Transposed (Format Direct 'Contiguous ('S 'Z) rep) =+ Format Direct 'Contiguous ('S 'Z) rep+type instance Transposed (Format Direct 'Strided ('S 'Z) rep ) =+ Format Direct 'Strided ('S 'Z) rep++type instance Transposed (Format Row 'Contiguous rank rep) =+ Format Column 'Contiguous rank rep+type instance Transposed (Format Row 'InnerContiguous rank rep) =+ Format Column 'InnerContiguous rank rep+type instance Transposed (Format Row 'Strided rank rep) =+ Format Column 'Strided rank rep++type instance Transposed (Format Column 'Contiguous rank rep)=+ Format Row 'Contiguous rank rep+type instance Transposed (Format Column 'InnerContiguous rank rep)=+ Format Row 'InnerContiguous rank rep+type instance Transposed (Format Column 'Strided rank rep)=+ Format Row 'Strided rank rep+++{-+a bunch of routines used to give various Layout operations for+array Formats that have DenseLayout instance+not exported or for human use+-}+{-# INLINE basicAddressRangeGeneric #-}+basicAddressRangeGeneric ::+ (Functor (Shape rank),Applicative (Shape rank),F.Foldable (Shape rank),+ DenseLayout form rank, Address~LayoutAddress form)=> form -> Maybe (Range Address)+basicAddressRangeGeneric = \ form ->+ if (fmap (flip (-) 1)$ basicLogicalShape form) `strictlyDominates` pure 0+ then Just $!+ Range (basicToDenseAddress form $! pure 0)+ (basicToDenseAddress form $!+ fmap (flip (-) 1) $! basicLogicalShape form)+ else Nothing++{-# INLINE basicToAddressDenseGeneric #-}+basicToAddressDenseGeneric :: (Functor (Shape rank),Applicative (Shape rank),F.Foldable (Shape rank),+ DenseLayout form rank,Address~LayoutAddress form) => form -> Shape rank Int -> Maybe Address+basicToAddressDenseGeneric = \ form ix ->+ if (fmap (flip (-) 1)$ basicLogicalShape form) `weaklyDominates` ix+ && ix `weaklyDominates` pure 0+ then Just $ basicToDenseAddress form ix+ else Nothing+{-# INLINE basicToIndexDenseGeneric #-}+basicToIndexDenseGeneric ::+ (Functor (Shape rank),F.Foldable (Shape rank),+ DenseLayout form rank,Address~LayoutAddress form) => form -> Address -> Shape rank Int+basicToIndexDenseGeneric = \form addr ->+ basicToDenseIndex form addr++{-# INLINE basicNextAddressDenseGeneric #-}+basicNextAddressDenseGeneric ::+ (Functor (Shape rank),F.Foldable (Shape rank),+ DenseLayout form rank,Address~LayoutAddress form) => form -> Address-> Maybe Address+basicNextAddressDenseGeneric = \ form addy ->+ case basicAddressRange form of+ Just (Range lo hi ) -> if addy >= lo && addy < hi+ then Just $! basicNextDenseAddress form addy+ else Nothing+ Nothing -> Nothing++{-# INLINE basicNextIndexDenseGeneric #-}+basicNextIndexDenseGeneric :: (Functor (Shape rank),F.Foldable (Shape rank),Applicative (Shape rank),+ DenseLayout form rank,Address~LayoutAddress form) =>+ form -> Shape rank Int -> Maybe Address ->Maybe (Shape rank Int,Address)+basicNextIndexDenseGeneric = \form ix _ ->+ if (fmap (flip (-) 1)$ basicLogicalShape form) `strictlyDominates` ix+ && ix `weaklyDominates` pure 0+ then+ Just $! basicNextDenseIndex form ix+ else+ Nothing++++{- | note that basicAffineAddressShiftGeneric may be suboptimal,+need to investigate what the core looks like+also TODO needs tests+-}+{-# INLINE basicAffineAddressShiftDenseGeneric #-}+basicAffineAddressShiftDenseGeneric :: (DenseLayout form rank+ ,DenseLayout (LayoutLogicalFormat form) rank+ ,Address~ LayoutAddress (LayoutLogicalFormat form))+ => form -> Address -> Int -> Maybe Address+basicAffineAddressShiftDenseGeneric form = \ addy shift ->+ let newForm = basicLogicalForm form in+ do+ nativeIndex <- return $ basicToDenseIndex form addy+ popBaseAddress <- return $ basicToDenseAddress newForm nativeIndex+ rng <- basicAddressRange newForm+ candidateAddress <- return $ popBaseAddress + Address shift+ if (getConst $ rangeMin ( Const) rng) <= candidateAddress+ && candidateAddress <= (getConst $ rangeMax ( Const) rng)+ then return $ basicToDenseAddress form $ basicToDenseIndex newForm candidateAddress+ else Nothing+++++-----+-----+-----+++type instance LayoutAddress (Format Direct 'Contiguous ('S 'Z) rep) = Address+type instance LayoutLogicalFormat (Format Direct 'Contiguous ('S 'Z) rep) = Format Direct 'Contiguous ('S 'Z) rep+instance Layout (Format Direct 'Contiguous ('S 'Z) rep) ('S 'Z) where+++ {-# INLINE basicLogicalShape #-}+ basicLogicalShape = \ x -> (logicalShapeDirectContiguous x) :* Nil++ basicLogicalForm = id++ transposedLayout = id++ {-# INLINE basicCompareIndex #-}+ basicCompareIndex = \ _ (l:* _) (r:* _) -> compare l r++ basicAddressRange = basicAddressRangeGeneric++ basicToAddress = basicToAddressDenseGeneric++ basicToIndex = basicToIndexDenseGeneric+++ basicNextAddress = basicNextAddressDenseGeneric++ basicNextIndex = basicNextIndexDenseGeneric+++ basicAddressPopCount = \ _ (Range (Address lo) (Address hi )) ->+ if hi >= lo then hi - lo+ else error $ "for basicAddressPopCount requires address obey hi >= lo, given: "+ ++ show hi ++ " " ++ show lo+ -- FIX me, add the range error checking+ -- in the style of the Sparse instances+++ basicAddressAsInt = \ _ (Address a) -> a++ basicAffineAddressShift = basicAffineAddressShiftDenseGeneric++ {-# INLINE basicAffineAddressShift #-}+ {-# INLINE basicAddressRange #-}+ {-# INLINE basicToAddress #-}+ {-# INLINE basicToIndex #-}+ {-# INLINE basicNextAddress #-}+ {-# INLINE basicNextIndex #-}+ {-# INLINE basicAddressPopCount #-}++type instance LayoutAddress (Format Direct 'Strided ('S 'Z) rep) = Address++instance Layout (Format Direct 'Strided ('S 'Z) rep) ('S 'Z) where++ {-# INLINE basicLogicalShape #-}+ basicLogicalShape = \x -> (logicalShapeDirectStrided x) :* Nil++ transposedLayout = id++ basicLogicalForm = (\ (n :* Nil ) -> FormatDirectContiguous n) . basicLogicalShape++ {-# INLINE basicCompareIndex #-}+ basicCompareIndex = \ _ (l:* _) (r:* _) -> compare l r++ basicAddressRange = basicAddressRangeGeneric++ basicToAddress = basicToAddressDenseGeneric++ basicToIndex = basicToIndexDenseGeneric++ basicNextAddress = basicNextAddressDenseGeneric++ basicNextIndex = basicNextIndexDenseGeneric++ basicAddressPopCount = \form@(FormatDirectStrided size _ ) (Range loA hiA)->+ let newForm = (FormatDirectContiguous size)+ in+ basicAddressPopCount newForm+ (Range (basicToDenseAddress newForm $ basicToDenseIndex form loA)+ (basicToDenseAddress newForm $ basicToDenseIndex form hiA) )++ basicAddressAsInt = \ _ (Address a) -> a++ basicAffineAddressShift = basicAffineAddressShiftDenseGeneric++ {-# INLINE basicAffineAddressShift #-}+ {-# INLINE basicAddressRange #-}+ {-# INLINE basicToAddress #-}+ {-# INLINE basicToIndex #-}+ {-# INLINE basicNextAddress #-}+ {-# INLINE basicNextIndex #-}+ {-# INLINE basicAddressPopCount #-}+++-- one type family instance for all the rows+type instance LayoutAddress (Format Row locality rank rep) = Address++instance (Applicative (Shape rank), Traversable (Shape rank))+ => Layout (Format Row 'Contiguous rank rep) rank where++ transposedLayout = \(FormatRowContiguous shp) -> FormatColumnContiguous $ reverseShape shp++ {-# INLINE basicLogicalShape #-}+ basicLogicalShape = boundsFormRow++ basicLogicalForm = id++ {-# INLINE basicCompareIndex #-}+ basicCompareIndex = \ _ ls rs -> foldl majorCompareLeftToRight EQ $ S.map2 compare ls rs++ basicAddressPopCount = \ _ (Range (Address lo) (Address hi )) -> hi - lo+ -- FIX me, add the range error checking+ -- in the style of the Sparse instances+ basicAddressRange = basicAddressRangeGeneric++ basicToAddress = basicToAddressDenseGeneric++ basicToIndex = basicToIndexDenseGeneric++ basicNextAddress = basicNextAddressDenseGeneric++ basicNextIndex = basicNextIndexDenseGeneric++ basicAddressAsInt = \ _ (Address a) -> a++ basicAffineAddressShift = basicAffineAddressShiftDenseGeneric++ {-# INLINE basicAffineAddressShift #-}+ {-# INLINE basicAddressRange #-}+ {-# INLINE basicToAddress #-}+ {-# INLINE basicToIndex #-}+ {-# INLINE basicNextAddress #-}+ {-# INLINE basicNextIndex #-}+ {-# INLINE basicAddressPopCount #-}++++instance (Applicative (Shape rank), Traversable (Shape rank))+ => Layout (Format Row 'InnerContiguous rank rep) rank where++ {-# INLINE basicLogicalShape #-}+ basicLogicalShape = boundsFormRowInnerContig++ basicLogicalForm form = FormatRowContiguous $ basicLogicalShape form++ transposedLayout = \(FormatRowInnerContiguous shp stride) ->+ FormatColumnInnerContiguous (reverseShape shp) (reverseShape stride)++ {-# INLINE basicCompareIndex #-}+ basicCompareIndex = \ _ ls rs ->+ foldl majorCompareLeftToRight EQ $ S.map2 compare ls rs++ basicAddressRange = basicAddressRangeGeneric++ basicToAddress = basicToAddressDenseGeneric++ basicToIndex = basicToIndexDenseGeneric++ basicNextAddress = basicNextAddressDenseGeneric++ basicNextIndex = basicNextIndexDenseGeneric++ basicAddressPopCount = \form@(FormatRowInnerContiguous size _) (Range loA hiA)->+ let newForm = (FormatRowContiguous size)+ in+ basicAddressPopCount newForm+ (Range (basicToDenseAddress newForm $ basicToDenseIndex form loA)+ (basicToDenseAddress newForm $ basicToDenseIndex form hiA) )++ basicAddressAsInt = \ _ (Address a) -> a++ basicAffineAddressShift = basicAffineAddressShiftDenseGeneric++ {-# INLINE basicAffineAddressShift #-}+ {-# INLINE basicAddressRange #-}+ {-# INLINE basicToAddress #-}+ {-# INLINE basicToIndex #-}+ {-# INLINE basicNextAddress #-}+ {-# INLINE basicNextIndex #-}+ {-# INLINE basicAddressPopCount #-}++++instance (Applicative (Shape rank),Traversable (Shape rank))+ => Layout (Format Row 'Strided rank rep) rank where++ {-# INLINE basicLogicalShape #-}+ basicLogicalShape = boundsFormRowStrided++ basicLogicalForm form = FormatRowContiguous $ basicLogicalShape form++ transposedLayout = \(FormatRowStrided shp stride) ->+ FormatColumnStrided (reverseShape shp) (reverseShape stride)++ {-# INLINE basicCompareIndex #-}+ basicCompareIndex = \ _ ls rs ->+ foldl majorCompareLeftToRight EQ $ S.map2 compare ls rs++ basicAddressRange = basicAddressRangeGeneric++ basicToAddress = basicToAddressDenseGeneric++ basicToIndex = basicToIndexDenseGeneric++ basicNextAddress = basicNextAddressDenseGeneric++ basicNextIndex = basicNextIndexDenseGeneric++ basicAddressPopCount = \form@(FormatRowStrided size _) (Range loA hiA)->+ let newForm = (FormatRowContiguous size)+ in+ basicAddressPopCount newForm+ (Range (basicToDenseAddress newForm $ basicToDenseIndex form loA)+ (basicToDenseAddress newForm $ basicToDenseIndex form hiA) )++ basicAddressAsInt = \ _ (Address a) -> a++ basicAffineAddressShift = basicAffineAddressShiftDenseGeneric++ {-# INLINE basicAffineAddressShift #-}+ {-# INLINE basicAddressRange #-}+ {-# INLINE basicToAddress #-}+ {-# INLINE basicToIndex #-}+ {-# INLINE basicNextAddress #-}+ {-# INLINE basicNextIndex #-}+ {-# INLINE basicAddressPopCount #-}+++type instance LayoutAddress (Format Column locality rank rep) = Address+instance (Applicative (Shape rank), Traversable (Shape rank))+ => Layout (Format Column 'Contiguous rank rep) rank where++ {-# INLINE basicLogicalShape #-}+ basicLogicalShape = boundsColumnContig++ basicLogicalForm = id++ transposedLayout = \(FormatColumnContiguous shp)-> FormatRowContiguous $ reverseShape shp++ {-# INLINE basicCompareIndex #-}+ basicCompareIndex = \ _ ls rs -> foldr majorCompareRightToLeft EQ $ S.map2 compare ls rs++ basicAddressPopCount = \ _ (Range (Address lo) (Address hi )) ->+ if hi >= lo then hi - lo+ else error $ "for basicAddressPopCount, require address hi >= lo, given: "+ ++ show hi ++ " " ++ show lo+ -- FIX me, add the range error checking+ -- in the style of the Sparse instances+ basicAddressRange = basicAddressRangeGeneric++ basicToAddress = basicToAddressDenseGeneric++ basicToIndex = basicToIndexDenseGeneric++ basicNextAddress = basicNextAddressDenseGeneric++ basicNextIndex = basicNextIndexDenseGeneric++ basicAddressAsInt = \ _ (Address a) -> a+++ basicAffineAddressShift = basicAffineAddressShiftDenseGeneric++ {-# INLINE basicAffineAddressShift #-}+ {-# INLINE basicAddressRange #-}+ {-# INLINE basicToAddress #-}+ {-# INLINE basicToIndex #-}+ {-# INLINE basicNextAddress #-}+ {-# INLINE basicNextIndex #-}+ {-# INLINE basicAddressPopCount #-}+++instance (Applicative (Shape rank), Traversable (Shape rank))+ => Layout (Format Column 'InnerContiguous rank rep) rank where+++ {-# INLINE basicLogicalShape #-}+ basicLogicalShape = boundsColumnInnerContig++ basicLogicalForm form = FormatColumnContiguous $ basicLogicalShape form++ transposedLayout = \(FormatColumnInnerContiguous shp stride)->+ FormatRowInnerContiguous (reverseShape shp) (reverseShape stride)++ {-# INLINE basicCompareIndex #-}+ basicCompareIndex = \ _ ls rs -> foldr majorCompareRightToLeft EQ $ S.map2 compare ls rs++ basicAddressRange = basicAddressRangeGeneric++ basicToAddress = basicToAddressDenseGeneric++ basicToIndex = basicToIndexDenseGeneric++ basicNextAddress= basicNextAddressDenseGeneric++ basicNextIndex= basicNextIndexDenseGeneric++ basicAddressPopCount = \form@(FormatColumnInnerContiguous size _) (Range loA hiA)->+ let newForm = (FormatColumnContiguous size)+ in+ basicAddressPopCount newForm+ (Range (basicToDenseAddress newForm $ basicToDenseIndex form loA)+ (basicToDenseAddress newForm $ basicToDenseIndex form hiA) )++ basicAddressAsInt = \ _ (Address a) -> a+ -- strideRow :: Shape rank Int,++ basicAffineAddressShift = basicAffineAddressShiftDenseGeneric++ {-# INLINE basicAffineAddressShift #-}+ {-# INLINE basicAddressRange #-}+ {-# INLINE basicToAddress #-}+ {-# INLINE basicToIndex #-}+ {-# INLINE basicNextAddress #-}+ {-# INLINE basicNextIndex #-}+ {-# INLINE basicAddressPopCount #-}++instance (Applicative (Shape rank), Traversable (Shape rank))+ => Layout (Format Column 'Strided rank rep) rank where++ {-# INLINE basicLogicalShape #-}+ basicLogicalShape = boundsColumnStrided++ basicLogicalForm form = FormatColumnContiguous $ basicLogicalShape form++ transposedLayout = \(FormatColumnStrided shp stride)->+ FormatRowStrided (reverseShape shp) (reverseShape stride)++ {-# INLINE basicCompareIndex #-}+ basicCompareIndex = \ _ ls rs -> foldr majorCompareRightToLeft EQ $ S.map2 compare ls rs++ basicAddressRange = basicAddressRangeGeneric++ basicToAddress = basicToAddressDenseGeneric++ basicToIndex = basicToIndexDenseGeneric++ basicNextAddress = basicNextAddressDenseGeneric++ basicNextIndex = basicNextIndexDenseGeneric++ basicAddressPopCount = \form@(FormatColumnStrided size _) (Range loA hiA)->+ let newForm = (FormatColumnContiguous size)+ in+ basicAddressPopCount newForm+ (Range (basicToDenseAddress newForm $ basicToDenseIndex form loA)+ (basicToDenseAddress newForm $ basicToDenseIndex form hiA) )++ basicAddressAsInt = \ _ (Address a) -> a++ basicAffineAddressShift = basicAffineAddressShiftDenseGeneric++ {-# INLINE basicAffineAddressShift #-}+ {-# INLINE basicAddressRange #-}+ {-# INLINE basicToAddress #-}+ {-# INLINE basicToIndex #-}+ {-# INLINE basicNextAddress #-}+ {-# INLINE basicNextIndex #-}+ {-# INLINE basicAddressPopCount #-}++----------------------+----------------------+-----+-----+----------------------+----------------------++++---+---+---++{-+these are factored out versions of the+various shared computations in both Row and Column Major+rank n Array format computations++-}++{-# INLINE computeStrideShape #-}+computeStrideShape ::+ ((Int -> State Int Int) -> Shape n Int -> State Int (Shape n Int )) -> Shape n Int -> Shape n Int+computeStrideShape = \trvse shp ->+ flip evalState 1 $+ flip trvse shp $+ -- basically accumulating the product of the+ -- dimensions+ \ val ->+ do accum <- get ;+ put $! (val * accum) ;+ return accum;+++++-----+-----+-----++instance DenseLayout (Format Direct 'Contiguous ('S 'Z) rep) ('S 'Z) where+++ --maxDenseAddress = \ (FormatDirectContiguous ix) -> Address (ix -1)+++ {-#INLINE basicToDenseAddress #-}+ basicToDenseAddress = \ (FormatDirectContiguous _) (j :* _ ) -> Address j++ --basicNextIndex= undefined -- \ _ x -> Just $! x + 1+ --note its unchecked!+ {-# INLINE basicToDenseIndex #-}+ basicToDenseIndex = \ (FormatDirectContiguous _) (Address ix) -> (ix ) :* Nil++ {-# INLINE basicNextDenseAddress #-}+ basicNextDenseAddress = \ _ addr -> addr + 1++++++instance DenseLayout (Format Direct 'Strided ('S 'Z) rep) ('S 'Z) where+++++ {-#INLINE basicToDenseAddress #-}+ basicToDenseAddress = \ (FormatDirectStrided _ strid) (j :* Nil )-> Address (strid * j)++ {-# INLINE basicNextDenseAddress #-}+ basicNextDenseAddress = \ (FormatDirectStrided _ strid) addr -> addr + Address strid++ {-# INLINE basicNextDenseIndex #-}+ basicNextDenseIndex = \ form (i:* Nil ) -> (\ix -> (ix,basicToDenseAddress form ix)) $! (i + 1 :* Nil )+++ {-# INLINE basicToDenseIndex #-}+ basicToDenseIndex = \ (FormatDirectStrided _ stride) (Address ix) -> (ix `div` stride ) :* Nil+++-----+-----+-----++++-- strideRow :: Shape rank Int,+instance (Applicative (Shape rank),F.Foldable (Shape rank), Traversable (Shape rank)) =>+ DenseLayout (Format Row 'Contiguous rank rep) rank where++{-+TODO AUDIT++-}+ {-# INLINE basicToDenseAddress #-}+ --basicToAddress = \rs tup -> let !strider =takePrefix $! S.scanr (*) 1 (boundsFormRow rs)+ basicToDenseAddress = \rs tup ->+ let !strider = computeStrideShape traverse (boundsFormRow rs)+ in Address $! S.foldl' (+) 0 $! map2 (*) strider tup++ {-# INLINE basicNextDenseAddress #-}+ basicNextDenseAddress = \_ addr -> addr + 1++ {-# INLINE basicToDenseIndex #-}+ basicToDenseIndex = \ rs (Address ix) ->+ let !striderShape = computeStrideShape traverse (boundsFormRow rs)++ in+ flip evalState ix $+ flip (S.backwards traverse) striderShape $+ -- want to start from largest stride (which is on the right)+ \ currentStride ->+ do remainderIx <- get ;+ let (!qt,!rm)= quotRem remainderIx currentStride+ put $! rm+ return qt;+++++-----+-----++-- strideRow :: Shape rank Int,+instance (Applicative (Shape rank),F.Foldable (Shape rank), Traversable (Shape rank))+ => DenseLayout (Format Row 'InnerContiguous rank rep) rank where+++ {-# INLINE basicToDenseAddress #-}+ basicToDenseAddress = \rs tup ->+ Address $! S.foldl' (+) 0 $!+ map2 (*) (strideFormRowInnerContig rs ) tup++ {-# INLINE basicNextDenseIndex #-}+ basicNextDenseIndex = \ form@(FormatRowInnerContiguous shape _) ix ->+ --S.map snd $!+ (\index -> (index,basicToDenseAddress form index)) $!+ flip evalState 1 $+ for ((,) <$> ix <*> shape) $+ \(ixv ,shpv )->+ do carry <-get+ let (newCarry,modVal)=divMod (carry + ixv) shpv+ put $! newCarry+ return modVal+++ {-# INLINE basicToDenseIndex #-}+ basicToDenseIndex = \ rs (Address ix) -> flip evalState ix $+ flip ( S.backwards traverse) (strideFormRowInnerContig rs ) $+ \ currentStride ->+ do remainderIx <- get ;+ let (!qt,!rm)= quotRem remainderIx currentStride+ put $! rm+ return qt;++++---+---+-- strideRow :: Shape rank Int,++instance (Applicative (Shape rank),F.Foldable (Shape rank), Traversable (Shape rank))+ => DenseLayout (Format Row 'Strided rank rep) rank where++++ {-# INLINE basicToDenseAddress #-}+ basicToDenseAddress = \rs tup -> Address $!+ S.foldl' (+) 0 $! map2 (*) (strideFormRowStrided rs ) tup++ {-# INLINE basicNextDenseIndex #-}+ basicNextDenseIndex = \ form@(FormatRowStrided shape _) ix ->+ (\index -> (index,basicToDenseAddress form index)) $!+ flip evalState 1 $+ for ((,) <$> ix <*> shape) $+ \(ixv ,shpv )->+ do carry <-get+ let (newCarry,modVal)=divMod (carry + ixv) shpv+ put $! newCarry+ return modVal+++ {-# INLINE basicToDenseIndex #-}+ basicToDenseIndex = \ rs (Address ix) -> flip evalState ix $+ flip (S.backwards traverse ) (strideFormRowStrided rs ) $+ \ currentStride ->+ do remainderIx <- get ;+ let (!qt,!rm)= quotRem remainderIx currentStride+ put $! rm+ return qt;+++++-----+-----+-----+++ -- strideRow :: Shape rank Int,+instance (Applicative (Shape rank),F.Foldable (Shape rank), Traversable (Shape rank))+ => DenseLayout (Format Column 'Contiguous rank rep) rank where++++ {-# INLINE basicToDenseAddress #-}+ basicToDenseAddress = \rs tup ->+ let !strider = computeStrideShape (S.backwards traverse) (boundsColumnContig rs)+ in Address $! S.foldl' (+) 0 $! map2 (*) strider tup++ {-# INLINE basicNextDenseAddress #-}+ basicNextDenseAddress = \_ addr -> addr + 1++ {-# INLINE basicToDenseIndex #-}+ basicToDenseIndex = \ rs (Address ix) ->+ let !striderShape = computeStrideShape (S.backwards traverse) (boundsColumnContig rs)+ in+ flip evalState ix $+ for striderShape $+ \ currentStride ->+ do remainderIx <- get ;+ let (!qt,!rm)= quotRem remainderIx currentStride+ put $! rm+ return qt;++++++ -- strideRow :: Shape rank Int,+instance (Applicative (Shape rank),F.Foldable (Shape rank), Traversable (Shape rank))+ => DenseLayout (Format Column 'InnerContiguous rank rep) rank where+++ {-# INLINE basicToDenseAddress #-}+ basicToDenseAddress = \ form tup -> let !strider = strideFormColumnInnerContig form+ in Address $! foldl' (+) 0 $! map2 (*) strider tup+ {-# INLINE basicNextDenseIndex #-}+ basicNextDenseIndex = \ form@(FormatColumnInnerContiguous shape _) ix ->+ --S.map snd $!+ (\index -> (index,basicToDenseAddress form index)) $!+ flip evalState 1 $+ flip (S.backwards traverse) ((,) <$> ix <*> shape) $+ \(ixv ,shpv )->+ do carry <-get+ let (newCarry,modVal)=divMod (carry + ixv) shpv+ put $! newCarry+ return modVal+++ {-# INLINE basicToDenseIndex #-}+ basicToDenseIndex = \ rs (Address ix) -> flip evalState ix $+ flip S.traverse (strideFormColumnInnerContig rs ) $+ \ currentStride ->+ do remainderIx <- get ;+ let (!qt,!rm)= quotRem remainderIx currentStride+ put $! rm+ return qt;+++++instance (Applicative (Shape rank),F.Foldable (Shape rank), Traversable (Shape rank))+ => DenseLayout (Format Column 'Strided rank rep) rank where++ {-# INLINE basicToDenseAddress #-}+ basicToDenseAddress = \ form tup -> let !strider = strideFormColumnStrided form+ in Address $! foldl' (+) 0 $! map2 (*) strider tup++ {-# INLINE basicNextDenseIndex #-}+ basicNextDenseIndex = \ form@(FormatColumnStrided shape _) ix ->+ --S.map snd $!+ (\index -> (index,basicToDenseAddress form index)) $!+ flip evalState 1 $+ flip (S.backwards traverse) ((,) <$> ix <*> shape) $+ \(ixv ,shpv )->+ do carry <-get+ let (newCarry,modVal)=divMod (carry + ixv) shpv+ put $! newCarry+ return modVal+++ {-# INLINE basicToDenseIndex #-}+ basicToDenseIndex = \ rs (Address ix) -> flip evalState ix $+ flip S.traverse (strideFormColumnStrided rs ) $+ \ currentStride ->+ do remainderIx <- get ;+ let (!qt,!rm)= quotRem remainderIx currentStride+ put $! rm+ return qt;+++++++{-+*Numerical.Array.Layout> basicToAddress (FormColumn (2 :* 3 :* 7 :* Nil)) (0:* 2 :* 2 :* Nil)+Address 16+*Numerical.Array.Layout> basicToAddress (FormColumn (2 :* 3 :* 7 :* Nil)) (1:* 0 :* 0 :* Nil)+Address 1+*Numerical.Array.Layout> basicToAddress (FormColumn (2 :* 3 :* 7 :* Nil)) (0:* 0 :* 0 :* Nil)+Address 0+*Numerical.Array.Layout> basicToAddress (FormColumn (2 :* 3 :* 7 :* Nil)) (0:* 1 :* 0 :* Nil)+Address 2+*Numerical.Array.Layout> basicToAddress (FormColumn (2 :* 3 :* 7 :* Nil)) (0:* 0 :* 1 :* Nil)+-}++
+ src/Numerical/Array/Layout/Sparse.hs view
@@ -0,0 +1,913 @@+{-+the following (currently 5) sparse formats will live here+++DirectSparse 1dim++++one subtlety and a seemingly subtle point will be+that contiguous / inner contiguous sparse arrays+in 2dim (and 1dim) will have an ``inner dimension" shift int.+This is so that slices can be zero copy on *BOTH* the array of values,+and the Format indexing array machinery.++Note that in the 2dim case, it still wont quite be zero copy, because the+offsets into the inner dimension lookup table (not quite the right word)+will have to change when a general slice is used rather than a slice+that acts only on the outermost dimension.+-}++++-- {-# LANGUAGE PolyKinds #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE StandaloneDeriving#-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE NoImplicitPrelude #-}+{-# LANGUAGE UndecidableInstances #-}++#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 707+ {-# LANGUAGE AutoDeriveTypeable #-}+#endif+module Numerical.Array.Layout.Sparse(+ Layout(..)+ ,DirectSparse+ ,CSR+ ,CSC+ ,CompressedSparseRow+ ,CompressedSparseColumn -- FIX ME, re add column support later+ ,Format(FormatDirectSparseContiguous+ ,FormatContiguousCompressedSparseRow+ ,FormatInnerContiguousCompressedSparseRow+ ,FormatContiguousCompressedSparseColumn+ ,FormatInnerContiguousCompressedSparseColumn)+ ,ContiguousCompressedSparseMatrix(..)+ ,InnerContiguousCompressedSparseMatrix(..)+ ,module Numerical.Array.Layout.Base+ ) where++import Data.Data+import Data.Bits (unsafeShiftR)+import Control.Applicative+import Numerical.Array.Layout.Base+--import Numerical.Array.Shape+import Numerical.InternalUtils+import qualified Data.Vector.Generic as V+import Prelude hiding (error )+++data CompressedSparseRow+ deriving (Typeable)++type CSR = CompressedSparseRow++data CompressedSparseColumn+ deriving (Typeable)++type CSC = CompressedSparseColumn++data DirectSparse+ deriving (Typeable)++++data instance Format DirectSparse 'Contiguous ('S 'Z) rep =+ FormatDirectSparseContiguous {+ _logicalShapeDirectSparse:: {-# UNPACK#-} !Int+ ,_logicalBaseIndexShiftDirectSparse::{-# UNPACK#-} !Int+ ,_indexTableDirectSparse :: ! (BufferPure rep Int ) }+++deriving instance Show (BufferPure rep Int ) => Show (Format DirectSparse 'Contiguous ('S 'Z) rep)++++{-+for some listings of the design space of Sparse matrices+as found in other tools,+see < https://software.intel.com/en-us/mkl_11.2_ref >+and then navigate to the section "Sparse Matrix Storage Formats" within+"BLAS and Sparse BLAS Routines"++< http://netlib.org/linalg/html_templates/node90.html > is also pretty readable++theres a subtle detail about the invariants of contiguous vs inner inner contiguous+for CSR and CSC+when I do an inner contiguous / contiguous slice / projection,+what "address shifts" do i need to track to make sure the slices+are zero copy as much as possible++just slicing on the outer dimension doesn't need any row shifts,+but a generalized (a,b) ... (a+x,b+y) selection when a,b!=0 does need a inner+dim shift,++NOTE that translating the inner dimension table's addresses to the corresponding+value buffer's address can require a shift!+This will happen when doing a MajorAxis (outer dimension) slice+the picks out a Suffix of the CSR matrix's rows+++note that there are 2 formulations of CSR (/ CSC) formats++a) 3 array: value, column index, and row start vectors++b) 4 array: value, column index, rowstart, and row end vectors++lets use choice a) for contiguous vectors, and choice b) for+inner contiguous vectors.++In both cases we need to enrich the type with a "buffer shift"+to handle correctly doing lookups on submatrices picked out+by either a major axis slice++-}+++--deriving instance (Show (Shape (S (S Z)) Int), Show (BufferPure rep Int) )+ -- => Show (Format CompressedSparseRow Contiguous (S (S Z)) rep)++--deriving instance (Eq (Shape (S (S Z)) Int), Eq (BufferPure rep Int) )+ -- => Eq (Format CompressedSparseRow Contiguous (S (S Z)) rep)++--deriving instance (Data (Shape (S (S Z)) Int), Data (BufferPure rep Int) )+ --- => Data (Format CompressedSparseRow Contiguous (S (S Z)) rep)++--deriving instance (Typeable (Shape (S (S Z)) Int ), Typeable (BufferPure rep Int) )+ -- => Typeable (Format CompressedSparseRow Contiguous (S (S Z)) rep)+ --deriving (Eq,Data,Typeable)+++{-+NOTE!!!!!+_logicalBaseIndexShiftDirectSparse (and friends)+are so that major axis slices can still use the same buffer,+(needed for both Contiguous and InnerContiguous cases).+So When looking up the Address for a value based upon its+Inner dimension, we need to *SUBTRACT* that shift+to get the correct offset index into the current SLICE.++NB: THIS IS A TERRRIBLE EXPLANATION, FIXMEEEEE++Phrased differently, This address shift is the *Discrepancy/Difference*+between the size of the elided prefix of the Vector and the starting+position of the manifest entries.++(Q: does this ever ever matter, or can i punt that to vector, and only+need this )+++This is kinda a good argument for not punting the Slicing on the raw buffers to+Vector, because it generally makes this a bit more subtle to think about+and someone IS going to implement something wrong this way!+++Another subtle and potentially confusing point is distinguishing between+Affine shifts in the Index Space vs the Address space.++Only the outer dimension lookup table shift is needed in the Contiguous+2dim case, but the 2dim InnerContiguous case is a bit more confusing+because of the potential for a slice along the inner dimension++Rank 1 sparse (like Direct sparse) is only Contiguous,+and either a) doesn't need a shift, or b) only needs an index shift+commensurate matching the leading implicit index of a Major Axis Slice+++theres a BIG corner case in most standard CSR / CSC formats which is+underspecified in most docs about CSC and CSR formats.+Consider Without loss of generality, CSR format+ 1) how are empty rows modeled/signaled?+ 2) if the last row is empty, how is that signaled?++2) The last row is signaled to be be empty by having+ the last entry of _outerDim2InnerDim buffer be set to >=+ length of _innerDimIndex buffer (ie >= 1 + largest index of _innerDimIndex)+1)++note that the outer index table has 1+#rows length, with the last one being the+length of the array++-}++data ContiguousCompressedSparseMatrix rep =+ FormatContiguousCompressedSparseInternal {+ -- does this need the index space shift for outer range slices???+ _outerDimContiguousSparseFormat :: {-# UNPACK #-} !Int+ ,_innerDimContiguousSparseFormat :: {-# UNPACK #-} !Int+ ,_innerDimIndexContiguousSparseFormat :: !(BufferPure rep Int)+ ,_outerDim2InnerDimContiguousSparseFormat:: ! (BufferPure rep Int )+ }+ deriving (Typeable)++deriving instance (Show (BufferPure rep Int))=> Show (ContiguousCompressedSparseMatrix rep)++{-+ outerDim innerDim innerTable outer2InnerStart+-}++++{-+for Row major Compressed Sparse (CSR)+the X dim (columns) are the inner dimension, and Y dim (rows) are the outer dim+-}++++data InnerContiguousCompressedSparseMatrix rep =+ FormatInnerContiguousCompressedSparseInternal {+ _outerDimInnerContiguousSparseFormat :: {-# UNPACK #-} !Int+ ,_innerDimInnerContiguousSparseFormat :: {-# UNPACK #-} !Int+ ,_innerDimIndexShiftInnerContiguousSparseFormat:: {-# UNPACK #-} !Int++ ,_innerDimIndexInnerContiguousSparseFormat :: !(BufferPure rep Int)+ ,_outerDim2InnerDimStartInnerContiguousSparseFormat:: ! (BufferPure rep Int )+ ,_outerDim2InnerDimEndInnerContiguousSparseFormat:: ! (BufferPure rep Int )+ }+ deriving Typeable++deriving instance (Show (BufferPure rep Int))=> Show (InnerContiguousCompressedSparseMatrix rep)+++newtype instance Format CompressedSparseRow 'Contiguous ('S ('S 'Z)) rep =+ FormatContiguousCompressedSparseRow {+ _getFormatContiguousCSR :: (ContiguousCompressedSparseMatrix rep) }++deriving instance Show (ContiguousCompressedSparseMatrix rep)+ => Show (Format CompressedSparseRow 'Contiguous ('S ('S 'Z)) rep)++newtype instance Format CompressedSparseColumn 'Contiguous ('S ('S 'Z)) rep =+ FormatContiguousCompressedSparseColumn {+ _getFormatContiguousCSC :: (ContiguousCompressedSparseMatrix rep) }++deriving instance Show (ContiguousCompressedSparseMatrix rep)+ => Show (Format CompressedSparseColumn 'Contiguous ('S ('S 'Z)) rep)++newtype instance Format CompressedSparseRow 'InnerContiguous ('S ('S 'Z)) rep =+ FormatInnerContiguousCompressedSparseRow {+ _getFormatInnerContiguousCSR :: (InnerContiguousCompressedSparseMatrix rep )+ }+deriving instance Show (InnerContiguousCompressedSparseMatrix rep )+ => Show (Format CompressedSparseRow 'InnerContiguous ('S ('S 'Z)) rep)++newtype instance Format CompressedSparseColumn 'InnerContiguous ('S ('S 'Z)) rep =+ FormatInnerContiguousCompressedSparseColumn {+ _getFormatInnerContiguousCSC :: (InnerContiguousCompressedSparseMatrix rep )+ }++deriving instance Show (InnerContiguousCompressedSparseMatrix rep )+ => Show (Format CompressedSparseColumn 'InnerContiguous ('S ('S 'Z)) rep)++ --deriving (Show,Eq,Data)++{-+ FormatInnerContiguous rowsize columnsize++-}+--newtype instance Format CompressedSparseColumn Contiguous (S (S Z)) rep =+-- FormatContiguousCompressedSparseColumn {+-- _getFormatContiguousCSC :: (ContiguousCompressedSparseMatrix rep)+-- }+ --deriving (Show,Eq,Data)++--newtype instance Format CompressedSparseColumn InnerContiguous (S (S Z)) rep =+-- FormatInnerContiguousCompressedSparseColumn {+-- _getFormatInnerContiguousCSC :: (InnerContiguousCompressedSparseMatrix rep)+-- }+-- --deriving (Show,Eq,Data)++--CSR and CSC go here, and their version of lookups and next address and next index+++++++-- Offset binary search --- cribbed with permission from+-- edward kmett's structured lib++{-+todo: theres some neat micro optimizations that are+possible If I know how indexed structures are paged aligned and what not+eg, when binary search, check both the first and last slot of a page I land on.+Also on >= Nehalem, pages are "paired" so if you land on the lower page, the+upper page is always loaded, etc etc. Not doing these for now.+++also should compare against search strategies defined in+the vector-algorithms package, namely the+galloping ones+-}+++++{-+-- Assuming @l <= h@. Returns @h@ if the predicate is never @True@ over @[l..h)@+-- requires p be a "monotonic" predicate (FFFFFTTTTT)+-}+bsearchUp :: (Int -> Bool) -> Int -> Int -> Int+bsearchUp p = go where+ go l h+ | l == h = l+ | p m = go l m+ | otherwise = go (m+1) h+ where hml = h - l+ m = l + unsafeShiftR hml 1 + unsafeShiftR hml 6+{-# INLINE bsearchUp #-}+{-+ Assuming @l <= h@. Returns @l@ if the predicate is never @True@ over @(l..h]@+ assumes predicate p is monotonic decreasing TTTTTFFFFF+ -}+bsearchDown :: (Int -> Bool) -> Int -> Int -> Int+bsearchDown p = go where+ go l h+ | l == h = l+ | p (m+1) = go (m+1) h+ | otherwise = go l m+ where hml = h - l+ m = l + unsafeShiftR hml 1 + unsafeShiftR hml 6+{-# INLINE bsearchDown #-}++{-+-- Assuming @l <= h@. Returns @h@ if the predicate is never @True@ over @[l..h)@+-- requires p be a "monotonic" predicate (FFFFFTTTTT)+-- does a linear scan on the first constant number of slots+(for now 97 because i had to pick a number thats ~ log MaxInt)+and then falls back to binary search.+Meant to have O(1) average case, O(log n) worst case+-}+basicHybridSearchUp :: (Int -> Bool ) -> Int -> Int -> Int+basicHybridSearchUp p = goCaseMe where+ goCaseMe l h | (h-l <= magicConstant) || p magicConstant+ {- either the range is short, OR+ we know match happens in the first magicConstant size subrange+ -}+ = linearSearchUp p l (min h magicConstant)+ | otherwise = bsearchUp p magicConstant h+{-# INLINE basicHybridSearchUp #-}+++basicHybridSearchDown :: (Int -> Bool)-> Int -> Int -> Int+basicHybridSearchDown p = goCaseMe where+ goCaseMe l h | (h-l <= magicConstant) || p (h- magicConstant)+ {- either the range is short, OR+ we know match happens in the first magicConstant size subrange+ -}+ = linearSearchDown p (max l (h - magicConstant)) h+ | otherwise = bsearchDown p l (h - magicConstant)+{-# INLINE basicHybridSearchDown #-}++{-+i chose 97 because it seemed like a number thats ~ log MaxInt always (within 4x)+And is a range that should stay in L1 cache sizes for most purposes+-}+magicConstant :: Int+magicConstant = 97+++-- Assuming @l <= h@. Returns @h@ if the predicate is never @True@ over @[l..h)@+linearSearchUp :: (Int -> Bool)-> Int -> Int -> Int+linearSearchUp p = go where+ go l h+ | l ==h = l+ | p l = l+ | otherwise = go (l+1) h+{-#INLINE linearSearchUp #-}++-- Assuming @l <= h@. Returns @l@ if the predicate is never @True@ over @(l..h]@+linearSearchDown :: (Int -> Bool)-> Int -> Int -> Int+linearSearchDown p = go where+ go l h+ | l ==h = l+ | p h = h+ | otherwise = go l (h-1)+{-#INLINE linearSearchDown #-}+++++--+-- now assumed each key is unique and ordered+--+-- Assuming @l <= h@. Returns @h@ if the predicate is never @True@ over @[l..h)@++-- should at some point try out a ternary search scheme to have even better+-- cache behavior (and benchmark of course)++searchOrd :: (Int -> Ordering) -> Int -> Int -> Int+searchOrd p = go where+ go l h+ | l == h = l+ | otherwise = case p m of+ LT -> go (m+1) h+ --- entry is less than target, go up!+ EQ -> m+ -- we're there! Finish early+ GT -> go l m+ -- entry is greater than target, go down!+ where hml = h - l+ m = l + unsafeShiftR hml 1 + unsafeShiftR hml 6+{-# INLINE searchOrd #-}++lookupExact :: (Ord k, V.Vector vec k) => vec k -> k -> Maybe Int+lookupExact ks key+ | j <- searchOrd (\i -> compare (ks V.! i) key) 0 (V.length ks - 1)+ , ks V.! j == key = Just $! j+ | otherwise = Nothing+{-# INLINE lookupExact #-}++lookupExactRange :: (Ord k, V.Vector vec k) => vec k -> k -> Int -> Int -> Maybe Int+lookupExactRange ks key lo hi+ | j <- searchOrd (\i -> compare (ks V.! i) key) lo hi+ , ks V.! j == key = Just $! j+ | otherwise = Nothing+{-# INLINE lookupExactRange #-}++--lookupLUB :: (Ord k, V.Vector vec k) => vec k -> k -> Maybe Int+--lookupLUB ks key+-- | j <- search (\i -> compare (ks V.! i) key) 0 (V.length ks - 1)+-- , ks V.! j <= key = Just $! j+-- | otherwise = Nothing+--{-# INLINE lookupLUB #-}++type instance Transposed (Format DirectSparse 'Contiguous ('S 'Z) rep )=+ (Format DirectSparse 'Contiguous ('S 'Z) rep )+++++type instance LayoutAddress (Format DirectSparse 'Contiguous ('S 'Z) rep) = Address+++instance V.Vector (BufferPure rep) Int+ => Layout (Format DirectSparse 'Contiguous ('S 'Z) rep ) ('S 'Z) where++ transposedLayout = id+ -- {-# INLINE transposedLayout #-}++ basicLogicalShape = \ form -> _logicalShapeDirectSparse form :* Nil+ -- {-# INLINE basicLogicalShape #-}++ basicCompareIndex = \ _ (a:* Nil) (b :* Nil) -> compare a b+ -- {-# INLINE basicCompareIndex #-}++ basicAddressRange = \form ->+ case (minAddress form , maxAddress form ) of+ (Just least, Just greatest) -> Just (Range least greatest )+ _ -> Nothing++ where+ minAddress =+ \ (FormatDirectSparseContiguous _ _ lookupTable) ->+ if V.length lookupTable >0 then Just $! Address 0 else Nothing++ maxAddress =+ \ (FormatDirectSparseContiguous _ _ lookupTable) ->+ if (V.length lookupTable >0 )+ then Just $! Address (V.length lookupTable - 1 )+ else Nothing++-- TODO, double check that im doing shift correctly+ {-# INLINE basicToAddress #-}+ basicToAddress =+ \ (FormatDirectSparseContiguous shape indexshift lookupTable) (ix:*_) ->+ if not (ix < shape && ix > 0 ) then Nothing+ else fmap Address $! lookupExact lookupTable (ix + indexshift)++ {-# INLINE basicToIndex #-}+ basicToIndex =+ \ (FormatDirectSparseContiguous _ shift lut) (Address addr) ->+ ((lut V.! addr ) - shift) :* Nil+ {-# INLINE basicAddressAsInt #-}+ basicAddressAsInt = \ _ (Address a) -> a++ {-# INLINE basicNextAddress #-}+ basicNextAddress =+ \ (FormatDirectSparseContiguous _ _ lut) (Address addr) ->+ if addr >= (V.length lut) then Nothing else Just (Address (addr+1))++ -- {-# INLINE basicAddressPopCount #-}+ basicAddressPopCount = \ form (Range loadr@(Address lo) hiadr@(Address hi)) ->+ if not ( lo <= hi ) then+ error $! "basicAddressPopCount was passed a bad Address Range " ++ show loadr ++" " ++ show hiadr+ else+ case basicAddressRange form of+ Nothing -> 0+ Just (Range (Address loBound) (Address hiBound)) ->+ if not $ (loBound<= lo ) && (hi <= hiBound)+ then error $!+ "basicAddressPopCount was passed a bad Address Range: "+ ++show lo++" "++ show hi++"\nwith format Address range"+ ++ show loBound ++ " " ++ show hiBound+ else hi - lo+++{-+ i've said it before, i'll say it again, scanning forward in the index space+ for sparse structures is really weird, :)++ NOTE: also need to remember to do those index space shifts for+ 1dim direct sparse, and test them thoroughly+-}+ -- {-# INLINE basicNextIndex #-}+ basicNextIndex =+ \form@(FormatDirectSparseContiguous size shift lut) (ix:*Nil) mebeAddress ->+ if ix >= size || ix >= (lut V.! (V.length lut -1) - shift ) then Nothing+ -- if ix is out of bounds or the last element, we're done!+ else+ let+ resAddr = Address $! bsearchUp (\lix-> ix < ((lut V.! lix)-shift) )+ 0 (V.length lut )+ in case mebeAddress of+ Nothing -> resAddr `seq` (Just (basicToIndex form resAddr , resAddr))+ -- Q: do i want the Index part of the tuple to be strict or not?+ -- leaving it lazy for now+ -- TODO / FIX / AUDIT ME / NOT SURE+ -- this is the fall back binary search based lookup++ Just (Address adr)->+ -- make sure the address hint is in bounds and+ -- is <= the current position+ if adr >0 && adr < (V.length lut -1) && ix >=((lut V.! adr )-shift)+ then+ -- by construction we know theres at least one applicable index+ -- thats+ let !nextAddr = Address $!+ basicHybridSearchUp+ (\lix-> ix < ((lut V.! lix)-shift ) )+ adr (V.length lut -1)+ in Just (basicToIndex form nextAddr , nextAddr)+ else+ resAddr `seq` (Just (basicToIndex form resAddr , resAddr))+++------------+------------++type instance Transposed (Format CompressedSparseRow 'Contiguous ('S ('S 'Z)) rep )=+ (Format CompressedSparseColumn 'Contiguous ('S ('S 'Z)) rep )+++type instance LayoutAddress (Format CompressedSparseRow 'Contiguous ('S ('S 'Z)) rep ) = SparseAddress++instance (V.Vector (BufferPure rep) Int )+ => Layout (Format CompressedSparseRow 'Contiguous ('S ('S 'Z)) rep ) ('S ('S 'Z)) where++ transposedLayout = \(FormatContiguousCompressedSparseRow repFormat) ->+ (FormatContiguousCompressedSparseColumn repFormat)+ {-# INLINE transposedLayout #-}+++ basicLogicalShape = \ form -> (_innerDimContiguousSparseFormat $ _getFormatContiguousCSR form ) :*+ ( _outerDimContiguousSparseFormat $ _getFormatContiguousCSR form ):* Nil+ -- x_ix :* y_ix+ {-# INLINE basicLogicalShape #-}+++ basicCompareIndex = \ _ as bs -> shapeCompareRightToLeft as bs+ {-# INLINE basicCompareIndex #-}+++ {-# INLINE basicAddressPopCount #-}+ basicAddressPopCount = \ form (Range (SparseAddress _ lo) (SparseAddress _ hi)) ->+ if not ( lo <= hi ) then+ error $! "basicAddressPopCount was passed a bad Address Range " ++ show lo ++" " ++ show hi+ else+ case basicAddressRange form of+ Nothing -> 0+ Just (Range (SparseAddress _ loBound) (SparseAddress _ hiBound)) ->+ if not $ (loBound<= lo ) && (hi <= hiBound)+ then error $!+ "basicAddressPopCount was passed a bad SparseAddress Range: "+ ++show lo++" "++ show hi++"\nwith format SparseAddress range"+ ++ show loBound ++ " " ++ show hiBound+ else hi - lo++ -- {-# INLINE rangedFormatAddress #-}+ basicAddressRange = \ form ->+ case (minAddress form,maxAddress form) of+ (Just least, Just greatest)-> Just (Range least greatest)+ _ -> Nothing++ where+ {-+ probably should deduplicate min/maxAddress+ -}+ minAddress =+ \(FormatContiguousCompressedSparseRow+ (FormatContiguousCompressedSparseInternal y_row_range x_col_range+ columnIndex rowStartIndex)) ->+ if y_row_range < 1 || x_col_range < 1|| (V.length columnIndex < 1)+ then Nothing+ else+ -- the value buffer has the invariant the the end points+ -- of the buffer MUST be valid in bounds values if length buffer > 0+ --SparseAddress $! 0 $! 0++ -- hoisted where into if branch as let so lets could be strict+ let+ !addrShift = columnIndex V.! 0++ -- for now assuming candidateRow is ALWAYS valid+ --- haven't proven this, FIXME+ !candidateRow= {-linearSearchUp-}+ basicHybridSearchUp nonZeroRow 0 (y_row_range-1 )+++ {- FIXME, to get the right complexity+ to linear search on first log #rows + 1 slots, then fall+ back to binary search+ punting for now because this probably wont matter than often++ the solution will be to replace linearSearchUp+ with a hybridSearchUp+ -}+ nonZeroRow =+ \ !row_ix ->+ -- the first row to satisfy this property+ (rowStartIndex V.! (row_ix+1) > rowStartIndex V.! row_ix)+ -- if the start index is >0, already past the min address row!+ || (rowStartIndex V.! row_ix) - addrShift > 0++ --else maxIxP1 > rowStartIndex V.! row_ix+ in Just $! SparseAddress candidateRow $! 0++ maxAddress =+ \(FormatContiguousCompressedSparseRow+ (FormatContiguousCompressedSparseInternal y_row_range x_col_range+ columnIndex rowStartIndex)) ->+ if y_row_range < 1 || x_col_range < 1|| (V.length columnIndex < 1)+ then Nothing+ else+ -- the value buffer has the invariant the the end points+ -- of the buffer MUST be valid in bounds values if length buffer > 0+ --SparseAddress $! 0 $! 0++ -- hoisted where into if branch as let so lets could be strict+ let+ !addrShift = columnIndex V.! 0+ !maxIxP1 = V.length columnIndex++ -- for now assuming candidateRow is ALWAYS valid+ --- haven't proven this, FIXME+ !candidateRow= {-linearSearchDown-}+ basicHybridSearchDown nonZeroRow 0 (y_row_range-1 )+++ {- FIXME, to get the right complexity+ to linear search on last log #rows + 1 slots, then fall+ back to binary search+ punting for now because this probably wont matter than often++ the solution will be to replace linearSearchDown+ with a hybridSearchDown+ -}+ nonZeroRow =+ \ !row_ix ->+ -- the first row to satisfy this property (going down from last row)+ (rowStartIndex V.! (row_ix+1) > rowStartIndex V.! row_ix)+ -- if the start index is >= maxIxP1, havent gone down to max addres yet+ -- if < maxIxp1, we're at or below the max address+ || (rowStartIndex V.! row_ix) - addrShift < maxIxP1++ --else maxIxP1 > rowStartIndex V.! row_ix+ in+ Just $!+ SparseAddress candidateRow $! (V.length columnIndex - 1 )++ -- \ (FormatContiguousCompressedSparseRow+ -- (FormatContiguousCompressedSparseInternal _ y_range+ -- columnIndex _)) ->+ -- SparseAddress (y_range - 1) (V.length columnIndex - 1 )++ {-# INLINE basicAddressAsInt #-}+ basicAddressAsInt = \ _ (SparseAddress _ addr)-> addr++ {-# INLINE basicToIndex #-}+ basicToIndex =+ \ (FormatContiguousCompressedSparseRow+ (FormatContiguousCompressedSparseInternal _ _ columnIndex _))+ (SparseAddress outer inner) ->+ (columnIndex V.! inner ) :* outer :* Nil+ -- outer is the row (y index) and inner is the lookup position for the x index+++{-+theres 3 cases for contiguous next address:+in the middle of a run on a fixed outer dimension,+need to bump the outer dimension, or we're at the end of the entire array++we make the VERY strong assumption that no illegal addresses are ever made!++note that for very very small sparse matrices, the branching will have some+overhead, but in general branch prediction should work out ok.+-}+ {-# INLINE basicNextAddress #-}+ basicNextAddress =+ \ (FormatContiguousCompressedSparseRow+ (FormatContiguousCompressedSparseInternal _ _+ columnIndex rowStartIndex))+ (SparseAddress outer inner) ->+ if inner < (V.length columnIndex -1)+ -- can advance further+ -- && ( outer == (y_row_range-1)+ --- either last row+ || ((inner +1) < (rowStartIndex V.! (outer + 1) - (rowStartIndex V.! 0 )))+ -- or our address is before the next row starts+ -- 3 vector CSR has a +1 slot at the end of the rowStartIndex++ then+ Just (SparseAddress outer (inner+1))+ else+ if inner == (V.length columnIndex -1)+ then Nothing+ else Just (SparseAddress (outer + 1) (inner + 1 ) )+++ -- {-# INLINE basicToAddress #-}+ basicToAddress =+ \ (FormatContiguousCompressedSparseRow+ (FormatContiguousCompressedSparseInternal y_row_range x_col_range+ columnIndex rowStartIndex))+ (ix_x:*ix_y :* _ ) ->+ if not (ix_x >= x_col_range || ix_y >=y_row_range )+ then+ -- slightly different logic when ix_y < range_y-1 vs == range_y-1+ -- because contiguous, don't need the index space shift though!+ let+ shift = (rowStartIndex V.! 0)+ checkIndex i =+ if (columnIndex V.!i) == ix_x+ then Just i+ else Nothing+ in+ (SparseAddress ix_y <$>) $!+ checkIndex =<<+ --- FIXME : need to check+ lookupExactRange columnIndex ix_x+ ((rowStartIndex V.! ix_y) - shift)+ ((rowStartIndex V.! (ix_y+1) ) - shift)++ else (Nothing :: Maybe SparseAddress )+++ -- {-# INLINE basicNextIndex #-}+ {- because nextIndex acts like a range query+ it doesn't make sense for inner loops+ -}+ basicNextIndex =+ \_form@(FormatContiguousCompressedSparseRow+ (FormatContiguousCompressedSparseInternal+ y_row_range x_col_range _columnIndex _rowStartIndex))+ _ix@(innerX :* outerY :*Nil) mebeSparseAddress ->+ if not $ (innerX >=0 && innerX < x_col_range ) && (outerY >= 0 && outerY < y_row_range)+ -- checking if index is inbounds for logical shape+ -- return Nothing if its out of bounds+ -- QUESTION: should it throw an error instead of returning nothing?+ then Nothing+ else+ case mebeSparseAddress of+ Nothing -> error "finish me "+ where+ {- Okay here we check if the proposed current index is manifest, or not+ Is it the right Row To search for the next index,+ Or if We need to search further along. This is the way that+ enables Usage of operations That give a complexity that is O(1)+ in the average/best case and O(log N )in the worst case++ The logical we do is roughly first check If there is an element+ strictly Greater than ix in next we are doing the successor+ That is within that Row And if so we can directly+ do a binary search therein+ -}+ _resRow = error "finish me "++ (Just (SparseAddress _innerix _outerix) )+ -> error "really finish me"+++ --case mebeAddress of+ -- Nothing ->+ -- let+ -- resAddr = Address $! bsearchUp (\lix-> ix < ((lut V.! lix)-shift) )+ -- 0 (V.length lut )+ -- in+ -- resAddr `seq` (Just (basicToIndex form resAddr , resAddr))+ -- -- Q: do i want the Index part of the tuple to be strict or not?+ -- -- leaving it lazy for now+ -- -- TODO / FIX / AUDIT ME / NOT SURE+ -- -- this is the fall back binary search based lookup+ -- Just (Address adr)->+ -- -- make sure the address hint is in bounds and+ -- -- is <= the current position+ -- if adr >0 && adr < (V.length lut -1) && ix >=((lut V.! adr )-shift)+ -- then+ -- -- by construction we know theres at least one applicable index+ -- -- thats+ -- let !nextAddr = Address $!+ -- basicHybridSearchUp+ -- (\lix-> ix < ((lut V.! lix)-shift ) )+ -- adr (V.length lut -1)+ -- in Just (basicToIndex form nextAddr , nextAddr)+ -- else+ -- resAddr `seq` (Just (basicToIndex form resAddr , resAddr))+++++--type instance Transposed (Format CompressedSparseRow InnerContiguous (S (S Z)) rep )=+-- (Format CompressedSparseColumn InnerContiguous (S (S Z)) rep )++--type instance Transposed (Format CompressedSparseColumn InnerContiguous (S (S Z)) rep )=+-- (Format CompressedSparseRow InnerContiguous (S (S Z)) rep )+++--instance Layout (Format CompressedSparseRow InnerContiguous (S (S Z)) rep ) (S (S Z)) where+-- transposedLayout = \(FormatInnerContiguousCompressedSparseRow a b c d e f) ->+-- (FormatInnerContiguousCompressedSparseColumn a b c d e f)+-- {-# INLINE transposedLayout #-}+-- basicFormShape = \ form -> logicalRowShapeInnerContiguousCSR form :*+-- logicalColumnShapeInnerContiguousCSR form :* Nil+-- {-# INLINE basicFormShape #-}+-- basicCompareIndex = \ _ as bs ->shapeCompareRightToLeft as bs+-- {-# INLINE basicCompareIndex#-}++++--instance (V.Vector (BufferPure rep) Int )+-- => SparseLayout (Format CompressedSparseRow InnerContiguous (S (S Z)) rep ) (S (S Z)) where++-- type LayoutAddress (Format CompressedSparseRow+-- InnerContiguous (S (S Z)) rep ) = SparseAddress++-- {-# INLINE minSparseAddress #-}+-- minSparseAddress = \_ -> SparseAddress 0 0++-- {-# INLINE maxSparseAddress#-}+-- maxSparseAddress =+-- \ (FormatInnerContiguousCompressedSparseInternal _ outer_dim_range _+-- innerDimIndex _) ->+-- SparseAddress (outer_dim_range - 1) (V.length innerDimIndex - 1 )+++-- {-#INLINE basicToIndex #-}+-- basicToIndex =+-- \ (FormatInnerContiguousCompressedSparseInternal _ _ _ innerDimIndex _)+-- (SparseAddress outer inner) -> (innerDimIndex V.! inner ) :* outer :* Nil+-- -- outer is the row (y index) and inner is the lookup position for the x index++++--theres 3 cases for contiguous next address:+--in the middle of a run on a fixed outer dimension,+--need to bump the outer dimension, or we're at the end of the entire array++--we make the VERY strong assumption that no illegal addresses are ever made!++--note that for very very small sparse matrices, the branching will have some+--overhead, but in general branch prediction should work out ok.++-- {-# INLINE basicNextAddress #-}+-- basicNextAddress =+-- \ (FormatInnerContiguousCompressedSparseRow+-- (FormatInnerContiguousCompressedSparseInternal _ _ _+-- columnIndex rowstartIndex))+-- (SparseAddress outer inner) ->+-- if not (inner == (V.length columnIndex -1)+-- {- && outer == (y_range-1) -}+-- || (inner +1) == (rowstartIndex V.! (outer + 1)))+-- then+-- Just (SparseAddress outer (inner+1))+-- else+-- if inner == (V.length columnIndex -1)+-- then Nothing+-- else Just (SparseAddress (outer + 1) (inner + 1 ) )++-- -- error "finish me damn it"+-- {-# INLINE basicToSparseAddress #-}+-- basicToSparseAddress =+-- \ (FormatInnerContiguousCompressedSparseRow+-- (FormatInnerContiguousCompressedSparseInternal x_range y_range addrShift+-- columnIndex rowstartIndex))+-- (ix_x:*ix_y :* _ ) ->+-- if not (ix_x >= x_range || ix_y >=y_range )+-- then+-- -- slightly different logic when ix_y < range_y-1 vs == range_y-1+-- -- because contiguous, don't need the index space shift though!+-- SparseAddress ix_y <$>+-- lookupExactRange columnIndex ix_x+-- -- ((rowstartIndex V.! ix_y) - addrShift)+-- (if ix_y < (y_range-1)+-- -- addr shift is for correcting from a major axis slice+-- then (rowstartIndex V.! (ix_y+1) ) - addrShift+-- else V.length columnIndex - 1 )+-- else (Nothing :: Maybe SparseAddress )
+ src/Numerical/Array/Locality.hs view
@@ -0,0 +1,81 @@+++{-# LANGUAGE DataKinds #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE StandaloneDeriving #-}+++module Numerical.Array.Locality(Locality(..),LocalityMax,LocalityMin) where++import Data.Data++data Locality = Contiguous | Strided | InnerContiguous+ deriving (Eq,Show,Read,Typeable,Data)++#if defined(__GLASGOW_HASKELL__) && ( __GLASGOW_HASKELL__ >= 707) && ( __GLASGOW_HASKELL__ < 709)+deriving instance Typeable 'Strided+deriving instance Typeable 'InnerContiguous+deriving instance Typeable 'Contiguous+#endif++#if defined(__GLASGOW_HASKELL__) && ( __GLASGOW_HASKELL__ >= 707)+type family LocalityMax (a :: Locality) (b :: Locality) :: Locality where+ LocalityMax 'Contiguous 'Contiguous = 'Contiguous+ LocalityMax 'Contiguous 'InnerContiguous = 'Contiguous+ LocalityMax 'Contiguous 'Strided = 'Contiguous+ LocalityMax 'InnerContiguous 'Contiguous = 'Contiguous+ LocalityMax 'Strided 'Contiguous = 'Contiguous+ LocalityMax 'InnerContiguous 'InnerContiguous = 'InnerContiguous+ LocalityMax 'InnerContiguous 'Strided = 'InnerContiguous+ LocalityMax 'Strided 'InnerContiguous = 'InnerContiguous+ LocalityMax 'Strided 'Strided = 'Strided+type family LocalityMin (a::Locality) (b ::Locality) :: Locality where+ LocalityMin 'Contiguous 'Contiguous = 'Contiguous+ LocalityMin 'Contiguous 'InnerContiguous = 'InnerContiguous+ LocalityMin 'Contiguous 'Strided = 'Strided+ LocalityMin 'InnerContiguous 'Contiguous = 'InnerContiguous+ LocalityMin 'Strided 'Contiguous = 'Strided+ LocalityMin 'InnerContiguous 'InnerContiguous = 'InnerContiguous+ LocalityMin 'InnerContiguous 'Strided = 'Strided+ LocalityMin 'Strided 'InnerContiguous = 'Strided+ LocalityMin 'Strided 'Strided = 'Strided++#else+type family LocalityMax (a :: Locality) (b :: Locality) :: Locality+type instance LocalityMax a b = LocalityMaxPrivate a b++type family LocalityMaxPrivate (a :: Locality) (b :: Locality) :: Locality+type instance LocalityMaxPrivate 'Contiguous 'Contiguous = 'Contiguous+type instance LocalityMaxPrivate 'Contiguous 'InnerContiguous = 'Contiguous+type instance LocalityMaxPrivate 'Contiguous 'Strided = 'Contiguous+type instance LocalityMaxPrivate 'InnerContiguous 'Contiguous = 'Contiguous+type instance LocalityMaxPrivate 'Strided 'Contiguous = 'Contiguous+type instance LocalityMaxPrivate 'InnerContiguous 'InnerContiguous = 'InnerContiguous+type instance LocalityMaxPrivate 'InnerContiguous 'Strided = 'InnerContiguous+type instance LocalityMaxPrivate 'Strided 'InnerContiguous = 'InnerContiguous+type instance LocalityMaxPrivate 'Strided 'Strided = 'Strided++type family LocalityMin (a::Locality) (b ::Locality) :: Locality+type instance LocalityMin a b = LocalityMinPrivate a b+++type family LocalityMinPrivate (a::Locality) (b ::Locality) :: Locality+type instance LocalityMinPrivate 'Contiguous 'Contiguous = 'Contiguous+type instance LocalityMinPrivate 'Contiguous 'InnerContiguous = 'InnerContiguous+type instance LocalityMinPrivate 'Contiguous 'Strided = 'Strided+type instance LocalityMinPrivate 'InnerContiguous 'Contiguous = 'InnerContiguous+type instance LocalityMinPrivate 'Strided 'Contiguous = 'Strided+type instance LocalityMinPrivate 'InnerContiguous 'InnerContiguous = 'InnerContiguous+type instance LocalityMinPrivate 'InnerContiguous 'Strided = 'Strided+type instance LocalityMinPrivate 'Strided 'InnerContiguous = 'Strided+type instance LocalityMinPrivate 'Strided 'Strided = 'Strided+++#endif+++++
+ src/Numerical/Array/Mutable.hs view
@@ -0,0 +1,411 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE ScopedTypeVariables#-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE FunctionalDependencies #-}+++module Numerical.Array.Mutable(+ MArray(..)+ ,Array(..)+ ,RectilinearArray(..)+ ,DenseArrayBuilder(..)+ ,DenseArray(..)+ ,Boxed+ ,Unboxed+ ,Stored+ --,module Numerical.Array.Layout+ ,module Numerical.Array.Shape+ ) where++import Control.Monad.Primitive ( PrimMonad, PrimState )+--import qualified Numerical.Array.DenseLayout as L+import Numerical.Array.Address+import qualified Numerical.Array.Layout as L+import Numerical.Array.Layout (Layout,Locality(..),LayoutAddress,Format(..),Range(..),AffineRange(..))+import Numerical.Array.Shape+--import Numerical.Nat+--import GHC.Prim(Constraint)+import Numerical.World+--import Numerical.Array.Range+--import Numerical.Array.Storage(Boxed,Unboxed,Stored)+--import Numerical.Array.Locality++import qualified Numerical.Array.Pure as P+import qualified Numerical.Array.Storage as S+import Numerical.Array.Storage (Buffer,Boxed,Unboxed,Stored)+import qualified Data.Vector.Generic as VG+import qualified Data.Vector.Generic.Mutable as VGM++import Control.Monad (liftM)+--import qualified Data.Vector.Storable.Mutable as SM+--import qualified Data.Vector.Unboxed.Mutable as UM+--import qualified Data.Vector.Mutable as BM++{-+For now we're going to just crib the vector style api and Lift it+up into a multi dimensional setting.++the tentative design is to have something like++++you'd think that the following array type is ``right''+but then you'll hit problems supporting+-}++-- data MArray world rep lay (view:: Locality) rank elm where+-- MArray+-- {_marrBuffer :: {-# UNPACK #!(MBuffer world rep elm)+-- ,_marrForm :: {-# UNPACK #- } !(Form lay loc rank)+-- --,_marrShift :: {-# UNPACK #- } !Address+-- }+++-- shift will be zero for most reps i'll ever care about, but in certain cases,+-- might not be. So for now not including it, but might be needed later,+-- though likely in regards to some sparse format of some sort.+--Omitting it for now, but may need to revisit later!+--+--For now any 'Address' shift will need to be via the buffer+--+-- One ssue in the formats is ``logical'' vs ``manifest'' Address.+--+--+--we eedto have 'RepConstraint' be decoupled from the type class instances+-- because we to sometimes have things that are world parametric+--+-- indexing should be oblivious to locality,++++--NB: one important assumption we'll have for now, is that every+++-- dsfdf+--type family RepConstraint world rep el :: Constraint+--type instance MArrayElem++{- | 'MArray' is the generic data family that+-}+data family MArray world rep lay (view::Locality) (rank :: Nat ) st el++data instance MArray Native rep lay locality rank st el =+ MutableNativeArray {+ nativeBuffer :: ! (S.BufferMut rep st el )+ ,nativeFormat :: ! (Format lay locality rank rep)+ }+++-- | Every 'MutableArray' instance has a contiguous version+-- of itself, This contiguous version will ALWAYS have a Builder instance.+type family MutableArrayContiguous (marr :: * -> * -> *) :: * -> * -> *+type instance MutableArrayContiguous (MArray world rep layout locality rank)= MArray world rep layout 'Contiguous rank++-- | Sadly 'ArrMutable' will have to have instances written by hand for now+-- May later migrate the freeze / thaw machinery to Array.Phased, but lets+type family ArrMutable ( arr :: * -> * ) :: * -> * -> *++class P.PureArray (ArrPure marr) rank a => Array marr (rank:: Nat) a | marr -> rank where++ type ArrPure (marr :: * -> * -> * ) :: * -> *++ -- the type of the underlying storage buffer+ --type MutableArrayBuffer marr :: * -> * -> *++ -- really shouldnt appear in end user code, will only+ -- come up in writing new combinators+ -- the abstraction here is a reflection of the need for+ type MArrayAddress (marr :: * -> * -> * ) :: *++ -- | 'basicUnsafeAffineAddressShift' is needed to handle abstracting access in popcount space+ basicUnsafeAffineAddressShift :: (address ~ MArrayAddress marr) => marr st a -> Int -> address -> address+ -- question, should the type be -> address or -> Maybe address++ -- | Unsafely convert a mutable Array to its immutable version without copying.+ -- The mutable Array may not be used after this operation. Assumed O(1) complexity+ basicUnsafeFreeze :: (PrimMonad m, arr ~ ArrPure marr, marr ~ ArrMutable arr)+ => marr (PrimState m) a -> m (arr a)++ -- | Unsafely convert a pure Array to its mutable version without copying.+ -- the pure array may not be used after this operation. Assumed O(1) complexity+ basicUnsafeThaw :: (PrimMonad m, marr ~ ArrMutable arr, arr ~ ArrPure marr )+ => arr a -> m (marr (PrimState m) a)++ -- | gives the shape, a 'rank' length list of the dimensions+ basicShape :: marr st a -> Index rank++ -- | 'basicCardinality' reports the number of manifest addresses/entries are+ -- in the array in a given address sub range.+ -- This is useful for determining when to switch from a recursive algorithm+ -- to a direct algorithm.+ -- Should this be renamed to something like basicPopCount/+ basicCardinality ::(address ~ MArrayAddress marr) => marr st a -> Range address -> Int++ --basicUnsafeRead :: PrimMonad m => marr (PrimState m) a -> Shape rank Int -> m (Maybe a)++ -- | basicMutableSparseIndexToAddres checks if a index is present or not+ -- helpful primitive for authoring codes for (un)structured sparse array format+ basicSparseIndexToAddress :: (address ~ MArrayAddress marr)+ => marr s a -> Index rank -> Maybe address++ -- | 'basicMutableAddressToIndex' assumes you only give it legal manifest addresses+ basicAddressToIndex :: (address ~ MArrayAddress marr) =>marr s a -> address -> Index rank++ -- | return the smallest and largest valid logical address+ basicAddressRange :: (address ~ MArrayAddress marr)=> marr st a -> Maybe (Range address)+++ -- | gives the next valid logical address+ -- undefined on invalid addresses and the greatest valid address.+ -- Note that for invalid addresses in between minAddress and maxAddress,+ -- will return the next valid address.++ basicSparseNextAddress :: (address ~ MArrayAddress marr)=> marr st a -> address -> Maybe address+++ -- I think the case could be made for a basicPreviousAddress opeeration++ -- | gives the next valid array index, the least valid index that is+ -- or+ basicSparseNextIndex ::(address ~ MArrayAddress marr)=>+ marr st a -> Index rank -> Maybe address -> Maybe ( Index rank, address)+++ -- | for a given valid address, @'basicAddressRegion' addr @ will return an AddressInterval+ -- that contains @addr@. This will be a singleton when the "maximal uniform stride interval"+ -- containing @addr@ has strictly less than 3 elements. Otherwise will return an Address range+ -- covering the maximal interval that will have cardinality at least 3.+ basicLocalAffineAddressRegion ::(address ~ MArrayAddress marr)+ => marr st a ->address -> AffineRange address++ -- | this doesn't quite fit in this class, but thats ok, will deal with that later+ basicOverlaps :: marr st a -> marr st a -> Bool++ -- | Reset all elements of the vector to some undefined value, clearing all+ -- references to external objects. This is usually a noop for unboxed+ -- vectors. This method should not be called directly, use 'clear' instead.+ basicClear :: PrimMonad m => marr (PrimState m) a -> m ()++ ---- | Yield the element at the given position. This method should not be+ ---- called directly, use 'unsafeRead' instead.+ basicUnsafeAddressRead :: (PrimMonad m ,address ~ MArrayAddress marr) =>+ marr (PrimState m) a -> address-> m a++ ---- | Replace the element at the given position. This method should not be+ ---- called directly, use 'unsafeAddressWrite' instead.+ basicUnsafeAddressWrite :: (PrimMonad m ,address ~ MArrayAddress marr) =>+ marr (PrimState m) a -> address -> a -> m ()+++ --note the sparsewrite and sparse read are "fused" versions of basicManifestAddress+ -- and address read and write. probably needs to be benchmarked! TODO++ -- | Yield the element at the given position. This method should not be+ -- called directly, use 'unsafeSparseRead' instead.+ basicUnsafeSparseRead :: PrimMonad m => marr (PrimState m) a ->+ Index rank -> m (Maybe a)++ -- Replace the element at the given position. This method should not be+ -- called directly, use 'unsafeWrite' instead.+ -- the following is the type that normal Array indexing,+ -- as folks are used to, lookslike+ -- its wrong+ --basicUnsafeSparseWrite :: PrimMonad m => marr (PrimState m) a ->+ -- Index rank -> m( Maybe (a -> m ()))+-- this might get axed+++instance (Buffer rep el, Layout (Format lay locality rank rep) rank )+ =>Array (MArray Native rep lay locality rank) rank el where++ type ArrPure (MArray Native rep lay locality rank)= P.ImmArray Native rep lay locality rank++ type MArrayAddress (MArray Native rep lay locality rank)= LayoutAddress (Format lay locality rank rep)++ {-# INLINE basicShape #-}+ basicShape = L.basicLogicalShape . nativeFormat++ {-# NOINLINE basicUnsafeFreeze #-}+ basicUnsafeFreeze = \marr -> do+ pureBuffer <- VG.unsafeFreeze $ nativeBuffer marr+ return $ P.ImMutableNativeArray pureBuffer $ nativeFormat marr++ {-# NOINLINE basicUnsafeThaw #-}+ basicUnsafeThaw = \parr -> do+ mutBuffer <- VG.unsafeThaw $ P.nativeBufferPure parr+ return $ MutableNativeArray mutBuffer $ P.nativeFormatPure parr++ {-# INLINE basicSparseIndexToAddress #-}+ basicSparseIndexToAddress = \ marr -> L.basicToAddress (nativeFormat marr)++ {-# INLINE basicAddressToIndex #-}+ basicAddressToIndex = \ marr -> L.basicToIndex (nativeFormat marr)++ {-# INLINE basicSparseNextAddress #-}+ basicSparseNextAddress = \marr -> L.basicNextAddress (nativeFormat marr)++ {-# INLINE basicSparseNextIndex #-}+ basicSparseNextIndex = \marr -> L.basicNextIndex (nativeFormat marr)++ basicOverlaps = \marr1 marr2 -> VGM.overlaps (nativeBuffer marr1) (nativeBuffer marr2)++ basicClear = \marr -> VGM.clear (nativeBuffer marr)++ {-# INLINE basicUnsafeAddressRead #-}+ basicUnsafeAddressRead = \marr addr ->+ VGM.unsafeRead (nativeBuffer marr) (L.basicAddressAsInt (nativeFormat marr) addr)++ {-# INLINE basicUnsafeAddressWrite #-}+ basicUnsafeAddressWrite = \marr addr v->+ VGM.unsafeWrite (nativeBuffer marr) (L.basicAddressAsInt (nativeFormat marr) addr) v++ {-# INLINE basicUnsafeSparseRead #-}+ basicUnsafeSparseRead = \marr ix -> do+ maddr <- return $ basicSparseIndexToAddress marr ix+ maybe (return Nothing) (\addr -> liftM Just $ basicUnsafeAddressRead marr addr ) maddr++ {-# INLINE basicAddressRange #-}+ basicAddressRange = \marr -> L.basicAddressRange (nativeFormat marr)++ basicCardinality = \marr -> L.basicAddressPopCount (nativeFormat marr)+++ basicUnsafeAffineAddressShift = error "carter needs to add this"+ basicLocalAffineAddressRegion = error "crter needs to add this"+{-+++type ArrPure marr :: * -> *++type MArrayAddress marr :: *++basicUnsafeAffineAddressShift :: (address ~ MArrayAddress marr) => marr st a -> Int -> address -> address++basicLocalAffineAddressRegion :: (address ~ MArrayAddress marr) => marr st a -> address -> AffineRange address++-}++++++++class ( Array marr rank a, P.PureDenseArray (ArrPure marr) rank a )=>+ DenseArray marr rank a | marr -> rank where+ -- | for Dense arrays, it is always easy to check if a given index is valid.+ -- this operation better have O(1) complexity or else!+ basicIndexInBounds :: marr st a -> Index rank -> Bool+++ --basicUnsafeAddressDenseRead :: PrimMonad m => marr (PrimState m) a -> Address-> m a++ -- i already have dense address indexing ?+ --basicUnsafeAddressDenseWrite :: PrimMonad m => marr (PrimState m) a -> Address -> a -> m ()++ -- | Yield the element at the given position. This method should not be+ -- called directly, use 'unsafeRead' instead.+ basicUnsafeDenseRead :: PrimMonad m => marr (PrimState m) a -> Index rank -> m a++ -- | Replace the element at the given position. This method should not be+ -- called directly, use 'unsafeWrite' instead.+ basicUnsafeDenseWrite :: PrimMonad m => marr (PrimState m) a -> Index rank -> a -> m ()+++ -- | gives the next valid logical address+ -- undefined on invalid addresses and the greatest valid address.+ -- Note that for invalid addresses in between minAddress and maxAddress,+ -- will return the next valid address.++ basicNextAddress :: marr st a -> Address -> Address+++ -- I think the case could be made for a basicPreviousAddress opeeration++ -- | gives the next valid array index+ -- undefined on invalid indices and the greatest valid index+ basicNextIndex :: marr st a -> Index rank -> Index rank++++++++{-++Mutable (Dense) Array Builder will only have contiguous instances+and only makes sense for dense arrays afaik++BE VERY THOUGHTFUL about what instances you write, or i'll be mad+++not including the general sparse building in the first release,+will include subsequently+-}++--class MutableArray marr (rank:: Nat) a => MutableArrayBuilder marr rank a where+ --basicBuildArray:: Index rank -> b++class DenseArray marr rank a => DenseArrayBuilder marr rank a where+ basicUnsafeNew :: PrimMonad m => Index rank -> m (marr (PrimState m) a)+ basicUnsafeReplicate :: PrimMonad m => Index rank -> a -> m (marr (PrimState m) a)+++++class RectilinearArray marr rank a | marr -> rank where++ -- | @'MutableRectilinearOrientation' marr@ should equal Row or Column for any sane choice+ -- of instance, because every MutableRectilinear instance will have a notion of+ -- what the nominal major axix will be.+ -- The intended use case is side condition constraints like+ -- @'MutableRectilinearOrientation' marr~Row)=> marr -> b @+ -- for operations where majorAxix projections are correct only for Row+ -- major formats. Such as Row based forward/backward substitution (triangular solvers)+ type MutableRectilinearOrientation marr :: *++ type MutableArrayDownRank marr ( st:: * ) a+++ -- | MutableInnerContigArray is the "meet" (minimum) of the locality level of marr and InnerContiguous.+ -- Thus both Contiguous and InnerContiguous are made InnerContiguous, and Strided stays Strided+ -- for now this makes sense to have in the MutableRectilinear class, though that may change.+ -- This could also be thought of as being the GLB (greatest lower bound) on locality+ type MutableInnerContigArray (marr :: * -> * -> *) st a++++ --type MutableArrayBuffer+ --not implementing this .. for now++ -- | @'basicSliceMajorAxis' arr (x,y)@ returns the sub array of the same rank,+ -- with the outermost (ie major axis) dimension of arr restricted to the+ -- (x,y) is an inclusive interval, MUST satisfy x<y , and be a valid+ -- subinterval of the major axis of arr.+ basicMutableSliceMajorAxis :: PrimMonad m => marr (PrimState m) a ->+ (Int,Int)-> m (marr (PrimState m) a)+ --but should it be primmonadic? nah, tis pure!++ -- | semantically, 'basicProjectMajorAxis' arr ix, is the rank reducing version of what+ -- basicSliceMajorAxis arr (ix,ix) would mean _if_ the (ix,ix) tuple was a legal major axis slice+ basicMutableProjectMajorAxis :: PrimMonad m =>marr (PrimState m) a+ -> Int -> m (MutableArrayDownRank marr (PrimState m) a )++ -- | @'basicMutableSlice' arr ix1 ix2@ picks out the (hyper) rectangle in dimension @rank@+ -- where ix1 is the minimal corner and ix2+ basicMutableSlice :: PrimMonad m => marr (PrimState m) a -> Index rank -> Index rank+ -> m (MutableInnerContigArray marr (PrimState m) a )++++
+ src/Numerical/Array/Pure.hs view
@@ -0,0 +1,170 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE ScopedTypeVariables#-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE FunctionalDependencies #-}++module Numerical.Array.Pure where+++--import Numerical.Array.Address+import qualified Numerical.Array.Layout as L++import Numerical.Array.Locality+import Numerical.Array.Shape+import Numerical.Array.Range+import Numerical.Array.Storage as S+import Numerical.World++import qualified Data.Vector.Generic as VG+++{-+a general question that you might ask is "what primops need have a monad constraint"++ie rather than having type a -> b, why are they type Monad m => a -> m b ?++the answer boils down to the following: most array types have+a *PURE* header data structure that can't be mutated,+that contains the Shape, extent, some handle/pointer to the associated underlying+buffer/datastructure. Any (even nominally pure) access to that potentially+mutable buffer should be mediated by a monad.++I further assume that the *structure* and *extent* of this underlying buffer cannot change.++That is, A valid address will always stay valid, even if after some mutation it may+ correspond to a *different* index than it did before.+-}++{-+Fix ME, these names are lame++ImmArray == immutable array++-}+data family ImmArray world rep lay (view::Locality) (rank :: Nat ) el++data instance ImmArray Native rep lay locality rank el =+ ImMutableNativeArray {+ nativeBufferPure :: ! (S.BufferPure rep el )+ ,nativeFormatPure :: ! (L.Format lay locality rank rep)+ }+++class PureArray arr (rank:: Nat) a | arr -> rank where+ type PureArrayAddress (arr :: * -> * ) :: *++ -- | gives the shape, a 'rank' length list of the dimensions+ basicShape :: arr a -> Index rank++ --basicUnsafeRead :: PrimMonad m => marr (PrimState m) a -> Shape rank Int -> m (Maybe a)++ -- | basicMutableSparseIndexToAddres checks if a index is present or not+ -- helpful primitive for authoring codes for (un)structured sparse array format+ -- FIXME : THIS IS A TERRIBLE NAME+ basicSparseIndexToAddress :: ( address ~PureArrayAddress arr) => arr a -> Index rank -> (Maybe address)++ -- |+ basicAddressToIndex :: (address ~PureArrayAddress arr) => arr a -> address -> (Index rank )++ -- | return the Range of valid logical addresses+ basicAddressRange :: (address ~PureArrayAddress arr)=> arr a -> Maybe (Range address)++++ -- | gives the next valid logical address+ -- undefined on invalid addresses and the greatest valid address.+ -- Note that for invalid addresses in between minAddress and maxAddress,+ -- will return the next valid address+ basicNextAddress :: (address ~PureArrayAddress arr)=> arr a -> address -> Maybe address++ -- I think the case could be made for a basicPreviousAddress opeeration++ -- | gives the next valid array index+ -- undefined on invalid indices and the greatest valid index+ basicNextIndex :: (address ~PureArrayAddress arr)=>+ arr a -> Index rank -> Maybe address -> Maybe ( Index rank, address)+++ -- | for a given valid address, @'basicAddressRegion' addr @ will return an AddressInterval+ -- that contains @addr@. This will be a singleton when the "maximal uniform stride interval"+ -- containing @addr@ has strictly less than 3 elements. Otherwise will return an Address range+ -- covering the maximal interval that will have cardinality at least 3.+++ --basicAddressRegion :: (address ~PureArrayAddress arr)=> arr a -> address -> UniformAddressInterval address++ ---- | Yield the element at the given position. This method should not be+ ---- called directly, use 'unsafeRead' instead.+ basicUnsafeAddressRead :: (Monad m , address ~PureArrayAddress arr)=> arr a -> address-> m a++++ -- | Yield the element at the given position. This method should not be+ -- called directly, use 'unsafeSparseRead' instead.+ basicUnsafeSparseRead :: Monad m => arr a -> Index rank -> m (Maybe a)++-- the catch all layout instance++instance (Buffer rep el , L.Layout (L.Format lay locality rank rep) rank)+ =>PureArray (ImmArray Native rep lay locality rank ) rank el where+ type PureArrayAddress (ImmArray Native rep lay locality rank )+ =L.LayoutAddress (L.Format lay locality rank rep)++ {-# INLINE basicShape #-}+ basicShape = L.basicLogicalShape . nativeFormatPure++ {-# INLINE basicSparseIndexToAddress #-}+ basicSparseIndexToAddress= L.basicToAddress . nativeFormatPure++ {-# INLINE basicAddressToIndex #-}+ basicAddressToIndex = L.basicToIndex . nativeFormatPure++ {-# INLINE basicAddressRange #-}+ basicAddressRange = L.basicAddressRange . nativeFormatPure++ {-# INLINE basicNextAddress #-}+ basicNextAddress= L.basicNextAddress . nativeFormatPure++ {-# INLINE basicNextIndex #-}+ basicNextIndex = L.basicNextIndex . nativeFormatPure++ {-# INLINE basicUnsafeSparseRead #-}+ basicUnsafeSparseRead =+ \ arr ix -> case basicSparseIndexToAddress arr ix of+ Nothing -> return Nothing+ (Just addr) -> basicUnsafeAddressRead arr addr >>= ( return . Just)++ {-# INLINE basicUnsafeAddressRead #-}+ basicUnsafeAddressRead =+ \ arr addr ->+ VG.basicUnsafeIndexM (nativeBufferPure arr)+ (L.basicAddressAsInt (nativeFormatPure arr) $ addr)++class PureArray arr rank a => PureDenseArray arr rank a where++ -- | 'basicIndexInBounds' is an O(1) bounds check.+ basicIndexInBounds :: arr a -> Index rank -> Bool++ -- |+ basicUnsafeAddressDenseRead :: (address ~PureArrayAddress arr,Monad m) => arr a -> address-> m a++++ -- | Yield the element at the given position. This method should not be+ -- called directly, use 'unsafeRead' instead.+ basicUnsafeDenseReadM :: Monad m => arr a -> Index rank -> m a++++++
+ src/Numerical/Array/Range.hs view
@@ -0,0 +1,74 @@++{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE DeriveFunctor, DeriveGeneric, DeriveFunctor #-}+{-# LANGUAGE DeriveFoldable,DeriveTraversable #-}+module Numerical.Array.Range (+ Range(..)+ ,AffineRange(..)+ ,HasRange(..)+ ,affineRangeStride) where++import Data.Data+import GHC.Generics+#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ < 709+import Data.Foldable+import Data.Traversable+#endif+++{-+not quite the right module for this notion of range, but lets+fix that later+-}+-- | whenever you are tempted to do a (lo,hi) tuple, use this instead+-- This should perhaps be made lazy, but strict for now.+data Range a =Range {_RangeMin :: !a+ ,_RangeMax :: !a}+ deriving (Eq,Show,Data,Typeable,Generic ,Foldable,Traversable,Functor)++class HasRange r a | r -> a where+ rangeMin:: Functor f => (a -> f a )-> r -> f r+ rangeMax:: Functor f => (a -> f a )-> r -> f r++instance HasRange (Range a) a where+ rangeMax = _rangeMax+ {-# INLINE rangeMax#-}++ rangeMin = _rangeMin+ {-# INLINE rangeMin #-}++instance HasRange (AffineRange a) a where+ rangeMin = _affineRangeMin+ {-# INLINE rangeMin #-}++ rangeMax = _affineRangeMax+ {-# INLINE rangeMax #-}++_rangeMin :: Functor f => (a -> f a)-> Range a -> f (Range a)+_rangeMin = \ fun rec -> fmap (\mup -> rec{_RangeMin= mup}) $ fun (_RangeMin rec )+{-# INLINE _rangeMin#-}++_rangeMax :: Functor f => (a -> f a) -> Range a -> f (Range a)+_rangeMax = \ fun rec -> fmap (\mup -> rec{_RangeMax= mup}) $ fun (_RangeMax rec )+{-# INLINE _rangeMax #-}++-- | this is uniform address interval by any other name+data AffineRange a = AffineRange{_AffineRangeMin :: !a+ ,_AffineRangeStride :: ! Int+ ,_AffineRangeMax :: !a}+ deriving (Eq,Show,Data,Generic,Typeable,Functor,Foldable,Traversable )++_affineRangeMin :: Functor f => (a-> f a) -> AffineRange a -> f (AffineRange a)+_affineRangeMin= \ fun rec -> fmap (\mup -> rec{_AffineRangeMin=mup}) $ fun (_AffineRangeMin rec)+{-# INLINE _affineRangeMin#-}++_affineRangeMax :: Functor f => (a -> f a) -> AffineRange a -> f (AffineRange a)+_affineRangeMax= \ fun rec -> fmap (\mup -> rec{_AffineRangeMax=mup}) $ fun (_AffineRangeMax rec)+{-# INLINE _affineRangeMax #-}++affineRangeStride :: Functor f => (Int -> f Int) -> AffineRange a -> f (AffineRange a)+affineRangeStride = \fun rec -> fmap (\mup -> rec{_AffineRangeStride=mup}) $ fun (_AffineRangeStride rec)+{-# INLINE affineRangeStride #-}
+ src/Numerical/Array/Shape.hs view
@@ -0,0 +1,671 @@+{-# LANGUAGE DataKinds, GADTs, TypeFamilies #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE ExplicitForAll #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE UnboxedTuples #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE StandaloneDeriving #-}+{-# LANGUAGE DeriveFunctor #-}+{-# LANGUAGE NoImplicitPrelude #-}+++module Numerical.Array.Shape(+ -- * Shape+ Shape(..)+ -- * Shape Utilities+ ,foldl+ ,foldr+ ,foldl'+ ,foldl1+ ,foldr1+ ,map+ ,map2+ ,reverseShape+ ,Nat(..)+ ,shapeSize+ ,SNat(..)+ ,weaklyDominates+ ,strictlyDominates+ ,shapeToList+ ,Index+ ,backwards+ -- * Unboxed Vector Morphism+ ,UnBoxedShapeMorphism(..)+ ,unShapeVector+ ,reShapeVector+ ,T.traverse+ --,T.Traversable(..)++ )+ where++--import Data.Data+import Data.Typeable+import Data.Data++import qualified Data.Functor as Fun+import qualified Data.Foldable as F+import qualified Control.Applicative as A+import Control.Monad (liftM)+import Control.Monad.ST (runST)+import qualified Data.Traversable as T++--import Control.NumericalMonad.State.Strict+import Control.NumericalApplicative.Backwards+++import Numerical.Nat+import qualified Data.Monoid as Monoid++import Prelude hiding (map,foldl,foldr,init,scanl,scanr,scanl1,scanr1,foldl1,foldr1)++import qualified Foreign.Storable as Store+import qualified Foreign.Ptr as Ptr++import qualified Data.Vector.Unboxed as UV+import qualified Data.Vector.Unboxed.Mutable as UVM+import qualified Data.Vector.Generic as GV+import qualified Data.Vector.Generic.Mutable as GMV++{-+Shape may get renamed to Index in the near future!++PSA: do not take the INLINE pragmas as a style suggestion.+The only reason for the INLINEs, SPECIALIZE and the+nonrecursive type class definitions of operations+in this module are because shape will be used in the inner loops of+array indexing heavy computations,++-}+++ {-+the concern basically boils down to "will it specialize / inline well"++ -}+++{-+should explore using the Reverse and Backwards transformers in the+Transformers package, but not right now++note also the *Tup operations could be done with a more general State monad+for the tupled accumulation parameter. If theres no perf regression, should+move to using that instead.++-}+++infixr 3 :*++type Index rank = Shape rank Int++data Shape (rank :: Nat) a where+ Nil :: Shape 'Z a+ (:*) :: !(a) -> !(Shape r a ) -> Shape ('S r) a+++deriving instance Typeable Shape+++nilShapeConstrRep :: Constr+nilShapeConstrRep = mkConstr shapeDataTypeRep "Nil" [] Prefix+consShapeConstrRep :: Constr+consShapeConstrRep = mkConstr shapeDataTypeRep ":*" [] Infix++shapeDataTypeRep :: DataType+shapeDataTypeRep = mkDataType "Numerical.Array.Shape.Shape" [nilShapeConstrRep,consShapeConstrRep]+++--deriving instance (Data a, Typeable n ) => Data (Shape n a)+-- gfoldl f z xs = gfoldl f z (shapeToList xs)++-- --gfoldl _ z Nil = z Nil+-- --gfoldl f z (x :* xs) = z (:*) `f` x `f` xs++-- I would like to have (Data (Shape n a)) but that seems tricky++instance (Data a,Typeable 'Z) => Data (Shape 'Z a) where+ gfoldl _ z Nil = z Nil+ gunfold _ z _ = z Nil -- not sure if _ z _ is the right one, but typechecks+ dataTypeOf _ = shapeDataTypeRep+ toConstr _ = nilShapeConstrRep++instance (Data a, Data (Shape n a), Typeable ('S n))=> Data (Shape ('S n) a ) where+ gfoldl k z (a :* b) = (z (:*) `k` a) `k` b+ gunfold k z _ = k (k (z (:*)))+ dataTypeOf _ = shapeDataTypeRep+ toConstr _ = consShapeConstrRep+++-- figure this out!+--look at http://hackage.haskell.org/package/HList-0.3.4.1/docs/src/Data-HList-Data.html+--and https://hackage.haskell.org/package/base-4.3.1.0/docs/Data-Data.html#nilConstr+-- for examples+--instance Data a => Data (Shape Z a) where+ --gfoldl++--deriving instance Data (Shape Z a)+--deriving instance (Data a,Data (Shape n a))=> (Data (Shape (S n) a))++{-+too much work to do data instance with pre 7.8 typeable+-}++++instance Eq (Shape 'Z a) where+ (==) _ _ = True+ {-#INLINE (==)#-}+instance (Eq a,Eq (Shape s a))=> Eq (Shape ('S s) a ) where+ (==) (a:* as) (b:* bs) = (a == b) && (as == bs )+ {-#INLINE (==)#-}+instance Show (Shape 'Z a) where+ show _ = "Nil"++instance (Show a, Show (Shape s a))=> Show (Shape ('S s) a) where+ show (a:* as) = show a ++ " :* " ++ show as++-- at some point also try data model that+-- has layout be dynamically reified, but for now+-- keep it phantom typed for sanity / forcing static dispatch.+-- NB: may need to make it more general at some future point+--data Strided r a lay = Strided { getStrides :: Shape r a }++-- may want to typeclassify this?+++shapeToList :: Shape n a -> [a]+shapeToList Nil = []+shapeToList (a:* as) = a : (shapeToList as )+++{-+the Traversable instance for shape needs both Z and S Z base+cases to interact nicely with the instances defined for+foldable+-}++instance T.Traversable (Shape 'Z) where+ traverse = \ _ Nil -> A.pure Nil+ {-# INLINE traverse #-}+ sequenceA = T.traverse id+ mapM f = A.unwrapMonad . T.traverse (A.WrapMonad . f)+ sequence = T.mapM id+ {-#INLINE sequenceA #-}+ {-#INLINE mapM #-}+ {-#INLINE sequence #-}+++instance T.Traversable (Shape ('S 'Z)) where+ traverse = \ f (a:* as) -> (:*) A.<$> f a A.<*> T.traverse f as+ {-# INLINE traverse #-}+ sequenceA = T.traverse id+ mapM f = A.unwrapMonad . T.traverse (A.WrapMonad . f)+ sequence = T.mapM id+ {-#INLINE sequenceA #-}+ {-#INLINE mapM #-}+ {-#INLINE sequence #-}++instance T.Traversable (Shape ('S n)) => T.Traversable (Shape ('S ('S n))) where+ traverse = \ f (a:* as) -> (:*) A.<$> f a A.<*> T.traverse f as+ {-#INLINE traverse #-}+ sequenceA = T.traverse id+ mapM f = A.unwrapMonad . T.traverse (A.WrapMonad . f)+ sequence = T.mapM id+ {-#INLINE sequenceA #-}+ {-#INLINE mapM #-}+ {-#INLINE sequence #-}++backwards :: (T.Traversable t, A.Applicative f) =>+ ((a -> Backwards f b) -> t a -> Backwards f (t b))+ -> ((a -> f b) -> t a -> f (t b))+backwards= \ traver f container ->+ forwards $ traver (\x -> Backwards $ f x) container+{-#INLINE backwards #-}+++++--instance Fun.Functor (Shape r) where+-- fmap = mapShape+-- {-#INLINE fmap #-}++instance Fun.Functor (Shape 'Z) where+ fmap = \ _ Nil -> Nil+ {-# INLINE fmap #-}++instance (Fun.Functor (Shape r)) => Fun.Functor (Shape ('S r)) where+ fmap = \ f (a :* rest) -> f a :* ( Fun.fmap f rest )+ {-# INLINE fmap #-}++instance A.Applicative (Shape 'Z) where+ pure = \ _ -> Nil+ {-# INLINE pure #-}+ (<*>) = \ _ _ -> Nil+ {-# INLINE (<*>) #-}++instance A.Applicative (Shape r)=> A.Applicative (Shape ('S r)) where+ pure = \ a -> a :* (A.pure a)+ {-# INLINE pure #-}+ (<*>) = \ (f:* fs) (a :* as) -> f a :* ((A.<*>)) fs as+ {-# INLINE (<*>) #-}++{-+only doing Foldable for ranks >= 1 does mean that+we dont get the cute "rank zero arrays are references"+property. But want foldr1 and foldl1 to always succeed++lets try having rank 0 anyways, i'll be happier if i can support it++-}++instance F.Foldable (Shape 'Z) where+ foldl' = \ _ !init _-> init+ foldr' = \ _ !init _ -> init+ foldl = \ _ init _-> init+ foldr = \ _ init _-> init+ foldMap = \ _f _col -> mempty+ {-# INLINE foldMap #-}+ {-# INLINE foldl #-}+ {-# INLINE foldr #-}+ {-# INLINE foldl' #-}+ {-# INLINE foldr' #-}+ foldr1 = \ _ _ -> error "you can't call foldr1 on a rank Z(ero) Shape"+ foldl1 = \_ _ -> error "you can't call foldl1 on a rank Z(ero) Shape"+++instance F.Foldable (Shape ('S 'Z)) where+ foldl' = \ f !init (a:*Nil)-> f init a+ foldr' = \ f !init (a:*Nil)-> f a init+ foldl = \ f init (a:*Nil)-> f init a+ foldr = \ f init (a:*Nil)-> f a init+ foldMap = \ f (a :* Nil ) -> f a+ {-# INLINE foldMap #-}+ {-# INLINE foldl #-}+ {-# INLINE foldr #-}+ {-# INLINE foldl' #-}+ {-# INLINE foldr' #-}+ foldr1 = \ _ (a:* Nil) -> a+ foldl1 = \ _ (a:* Nil) -> a+ {-# INLINE foldl1 #-}+ {-# INLINE foldr1 #-}+instance ( F.Foldable (Shape ('S r)) )=> F.Foldable (Shape ('S ('S r))) where+ foldl' = \ f init (a:* as) -> F.foldl' f (f init a) as+ foldr' = \f !init (a :* as ) -> f a $! F.foldr' f init as+ foldl = \ f init (a:* as) -> F.foldl' f (f init a) as+ foldr = \ f init (a:* as) -> f a $! F.foldr f init as+ foldl1 = \ f (a:* as) -> F.foldl' f a as+ foldr1 = \ f (a :* as) -> F.foldr' f a as+ foldMap = \ f (a :* as ) -> f a Monoid.<> F.foldMap f as+ {-# INLINE foldMap #-}+ {-# INLINE foldl #-}+ {-# INLINE foldr #-}+ {-# INLINE foldl' #-}+ {-# INLINE foldr' #-}+ {-# INLINE foldl1 #-}+ {-# INLINE foldr1 #-}++instance (Semigroup a, A.Applicative (Shape n))=> (Semigroup (Shape n a)) where+ (<>) = \ a b -> A.pure (<>) A.<*> a A.<*> b++instance (Monoid.Monoid a, A.Applicative (Shape n))=> Monoid.Monoid (Shape n a) where+ mempty = A.pure Monoid.mempty+ mappend = \ a b -> A.pure Monoid.mappend A.<*> a A.<*> b+++++{- when you lift a toral order onto vectors, you get+interesting partial order -}++-- | 'weaklyDominates' is the '<=' operator lifted onto a sized vector to+-- induce a partial order relation+weaklyDominates :: (Ord a, A.Applicative (Shape n), F.Foldable (Shape n) )=>+ Shape n a -> Shape n a -> Bool+weaklyDominates = \major minor -> foldl (&&) True $! map2 (>=) major minor+{-# INLINE weaklyDominates #-}++-- | 'strictlyDominates' is the '<' operator lifted onto a sized vector to+-- induce a partial order relation+strictlyDominates :: (Ord a, A.Applicative (Shape n), F.Foldable (Shape n) )=>+ Shape n a -> Shape n a -> Bool++strictlyDominates = \major minor -> foldl (&&) True $! map2 (>) major minor+{-# INLINE strictlyDominates #-}++{-# INLINE reverseShape #-}+reverseShape :: Shape n a -> Shape n a+reverseShape Nil = Nil+reverseShape r@(_ :* Nil)= r+reverseShape (a:* b :* Nil) = b:* a :* Nil+reverseShape (a:* b :* c:* Nil )= c :* b :* a :* Nil+reverseShape (a:* b :* c :* d :* Nil)= d :* c :* b :* a :* Nil+reverseShape list = go SZero Nil list+ where+ go :: SNat n1 -> Shape n1 a-> Shape n2 a -> Shape (n1 + n2) a+ go snat acc Nil = gcastWith (plus_id_r snat) acc+ go snat acc (h :* (t :: Shape n3 a)) =+ gcastWith (plus_succ_r snat (Proxy :: Proxy n3))+ (go (SSucc snat) (h :* acc) t)+++{-+TODO: abstract out all the different unrolled cases i have+++-}+++++{-# INLINE map2 #-}+map2 :: forall a b c r . (A.Applicative (Shape r))=> (a->b ->c) -> (Shape r a) -> (Shape r b) -> (Shape r c )+map2 = \ f shpa shpb -> f A.<$> shpa A.<*> shpb+++{-# INLINE map #-}+map:: forall a b r . (A.Applicative (Shape r))=> (a->b) -> (Shape r a )->( Shape r b)+map = \ f shp -> f A.<$> shp++++{-# INLINE foldr #-}+foldr :: forall a b r . (F.Foldable (Shape r))=> (a->b-> b) -> b -> Shape r a -> b+foldr = \ f init shp -> F.foldr f init shp+++++{-# INLINE foldl #-}+foldl :: forall a b r. (F.Foldable (Shape r))=> (b-> a -> b) -> b -> Shape r a -> b+foldl = \ f init shp -> F.foldl f init shp+++{-# INLINE foldl' #-}+foldl' :: forall a b r . (F.Foldable (Shape r))=> (b-> a -> b) -> b -> Shape r a -> b+foldl' = \ f init shp -> F.foldl' f init shp++{-# INLINE foldr1 #-}+foldr1 :: forall b r . (F.Foldable (Shape ('S r)))=> (b->b-> b) -> Shape ('S r) b -> b+foldr1 = \ f shp -> F.foldr1 f shp+++++{-# INLINE foldl1 #-}+foldl1 :: forall b r. (F.Foldable (Shape ('S r)))=> (b-> b -> b) -> Shape ('S r) b -> b+foldl1 = \ f shp -> F.foldl1 f shp++++++instance Store.Storable a =>Store.Storable (Shape ('S 'Z) a) where+ {-#INLINE sizeOf#-}+ sizeOf = \ _ -> (Store.sizeOf (undefined :: a))+ -- might want to boost the alignment, but ignore for now+ {-# INLINE alignment #-}+ alignment = \ _ -> Store.alignment (undefined :: a )+ {-# INLINE peek #-}+ peek = \ptr -> do res <- Store.peek (Ptr.castPtr ptr) ; return (res :* Nil)+ {-# INLINE poke #-}+ poke = \ptr (a:*_) -> Store.poke (Ptr.castPtr ptr) a+ {-# INLINE pokeElemOff #-}+ {-# INLINE peekElemOff #-}+ peekElemOff = \ ptr off -> Store.peekByteOff ptr (off * Store.sizeOf (undefined :: a ))+ pokeElemOff ptr off val = Store.pokeByteOff ptr (off * Store.sizeOf val) val++ peekByteOff ptr off = Store.peek (ptr `Ptr.plusPtr` off)+ pokeByteOff ptr off = Store.poke (ptr `Ptr.plusPtr` off)+ {-# INLINE peekByteOff #-}+ {-# INLINE pokeByteOff #-}+++instance (Store.Storable a,Store.Storable (Shape ('S n) a)) =>Store.Storable (Shape ('S ('S n)) a) where+ {-#INLINE sizeOf#-}+ sizeOf = \ _ -> Store.sizeOf (undefined :: a) + Store.sizeOf (undefined :: (Shape ('S n) a ))+ -- might want to boost the alignment, but ignore for now+ {-# INLINE alignment #-}+ alignment = \ _ -> Store.alignment (undefined :: a )+ {-# INLINE peek #-}+ peek = \ptr -> do+ a <- Store.peek (Ptr.castPtr ptr) ;+ as <- Store.peek (ptr `Ptr.plusPtr` Store.sizeOf (undefined :: a ))+ return (a:* as)+ {-# INLINE poke #-}+ poke = \ptr (a:*as ) -> do+ Store.poke (Ptr.castPtr ptr) a+ Store.poke (ptr `Ptr.plusPtr` Store.sizeOf (undefined :: a )) as+ {-# INLINE pokeElemOff #-}+ {-# INLINE peekElemOff #-}+ peekElemOff = \ ptr off -> Store.peekByteOff ptr (off * Store.sizeOf (undefined :: (Shape ('S ('S n)) a) ))+ pokeElemOff ptr off val = Store.pokeByteOff ptr (off * Store.sizeOf val) val++ peekByteOff ptr off = Store.peek (ptr `Ptr.plusPtr` off)+ pokeByteOff ptr off = Store.poke (ptr `Ptr.plusPtr` off)+ {-# INLINE peekByteOff #-}+ {-# INLINE pokeByteOff #-}++-- this instance is a bit weird and should never be used+-- but probably legal+instance Store.Storable a =>Store.Storable (Shape 'Z a) where+ {-#INLINE sizeOf#-}+ sizeOf = \ _ -> Store.sizeOf (undefined :: a )+ -- might want to boost the alignment, but ignore for now+ {-# INLINE alignment #-}+ alignment = \ _ -> Store.alignment (undefined :: a )+ {-# INLINE peek #-}+ peek = \ _ -> return Nil+ {-# INLINE poke #-}+ poke = \ _ _-> return ()+ {-# INLINE pokeElemOff #-}+ {-# INLINE peekElemOff #-}+ peekElemOff = \ _ _ -> return Nil+ pokeElemOff = \ _ _ _ -> return ()++ peekByteOff = \ _ _ -> return Nil+ pokeByteOff = \ _ _ _ -> return ()+ {-# INLINE peekByteOff #-}+ {-# INLINE pokeByteOff #-}++{-# INLINE shapeSize #-}+shapeSize :: F.Foldable (Shape n)=>Shape n a -> Int+shapeSize = \ as -> ( F.foldl (\ct _ -> ct +1) 0 as )++unShapeVector ::(UnBoxedShapeMorphism n a, T.Traversable (Shape n), UV.Unbox a) => UV.Vector (Shape n a) -> (Int, Shape n (UV.Vector a))+unShapeVector vs = runST $+ do (l,mvs) <- fmap unShapeMVector $ UV.unsafeThaw vs+ shpvs <- T.traverse UV.unsafeFreeze mvs+ return (l,shpvs)+++reShapeVector::(UnBoxedShapeMorphism n a, T.Traversable (Shape n), UV.Unbox a)=>+ (Int, Shape n (UV.Vector a)) -> UV.Vector (Shape n a)+reShapeVector (l,vs) = runST $+ do mShapeV <- T.traverse UV.unsafeThaw vs+ mvShp <- return $ reShapeMVector (l,mShapeV)+ UV.unsafeFreeze mvShp+++{- THis is a convenience type class so i dont have to export the constructors -}+class (UV.Unbox (Shape n a)) => UnBoxedShapeMorphism n a where+ --unShapeVector :: UV.Vector (Shape n a) -> (Int, Shape n (UV.Vector a))+ --reShapeVector :: (Int, Shape n (UV.Vector a)) -> UV.Vector (Shape n a)++ unShapeMVector :: UVM.MVector s (Shape n a) -> (Int, Shape n (UV.MVector s a))+ reShapeMVector :: (Int, Shape n (UVM.MVector s a)) -> UVM.MVector s (Shape n a)++instance (UV.Unbox a)=> UnBoxedShapeMorphism 'Z a where+ --unShapeVector (V_ShapeZ l)= (l,Nil)+ unShapeMVector (MV_ShapeZ l) = (l,Nil )++ --reShapeVector (l,Nil) = (V_ShapeZ l)+ reShapeMVector (l,Nil ) = (MV_ShapeZ l)+++instance (UV.Unbox a)=> UnBoxedShapeMorphism ('S 'Z) a where+ --unShapeVector (V_ShapeSZ v)= (UV.length v, v :* Nil)++ unShapeMVector (MV_ShapeSZ v) = (UVM.length v,v:* Nil )++ --reShapeVector (l,v :* Nil) = (V_ShapeSZ v)+ reShapeMVector (_,v :* _ ) = (MV_ShapeSZ v)++--UV.V_2+--UVM.MV_2+instance ((UV.Unbox a),UnBoxedShapeMorphism ('S n) a )=> UnBoxedShapeMorphism ('S ('S n)) a where+ --unShapeVector (V_ShapeSSN (UV.V_2 l vhead vtail))= (l, vhead :* snd (unShapeVector vtail) )+ unShapeMVector (MV_ShapeSSN (UVM.MV_2 l vhead vtail)) = (l,vhead:* snd (unShapeMVector vtail ))++ --reShapeVector (l,vh :* vt) = (V_ShapeSSN (UV.V_2 l vh (reShapeVector (l,vt) )))+ reShapeMVector (l,vh :* vt ) = (MV_ShapeSSN (UVM.MV_2 l vh (reShapeMVector (l,vt) )))++newtype instance UV.MVector s (Shape 'Z a) = MV_ShapeZ Int+newtype instance UV.Vector (Shape 'Z a) = V_ShapeZ Int++newtype instance UV.MVector s (Shape ('S 'Z) a) = MV_ShapeSZ (UV.MVector s a)+newtype instance UV.Vector (Shape ('S 'Z) a) = V_ShapeSZ (UV.Vector a)++newtype instance UV.MVector s (Shape ('S ('S n)) a) = MV_ShapeSSN (UV.MVector s (a, Shape ('S n) a) )+newtype instance UV.Vector (Shape ('S ('S n)) a) = V_ShapeSSN (UV.Vector (a, Shape ('S n) a) )+++instance UV.Unbox a => UV.Unbox (Shape 'Z a)+instance UV.Unbox a => UV.Unbox (Shape ('S 'Z) a)+instance (UV.Unbox a,UV.Unbox (Shape ('S n) a) )=> UV.Unbox (Shape ('S ('S n)) a)++++instance UV.Unbox a => GMV.MVector UV.MVector (Shape 'Z a) where+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicOverlaps #-}+ {-# INLINE basicUnsafeNew #-}+ {-# INLINE basicUnsafeRead #-}+ {-# INLINE basicUnsafeWrite #-}+ {-# INLINE basicClear #-}+ {-# INLINE basicSet #-}+ {-# INLINE basicUnsafeCopy #-}+ {-# INLINE basicUnsafeGrow #-}+ {-# INLINE basicInitialize #-}+ basicInitialize = \ (MV_ShapeZ _n) -> return ()+ basicLength = \ (MV_ShapeZ n) -> n+ basicUnsafeSlice = \ _ m (MV_ShapeZ _) -> MV_ShapeZ m+ basicOverlaps = \ _ _ -> False+ basicUnsafeNew = \ n -> return (MV_ShapeZ n)+ basicUnsafeRead = \ (MV_ShapeZ _) _ -> return Nil+ basicUnsafeWrite = \ (MV_ShapeZ _) _ Nil -> return ()+ basicClear = \ _ -> return ()+ basicSet = \ (MV_ShapeZ _) Nil -> return ()+ basicUnsafeCopy = \ (MV_ShapeZ _) (MV_ShapeZ _) -> return ()+ basicUnsafeGrow = \ (MV_ShapeZ n) m -> return $ MV_ShapeZ (n+m)++instance UV.Unbox a => GV.Vector UV.Vector (Shape 'Z a) where+ {-# INLINE basicUnsafeFreeze #-}+ basicUnsafeFreeze = \ (MV_ShapeZ n) -> return $ V_ShapeZ n+ {-# INLINE basicUnsafeThaw #-}+ basicUnsafeThaw = \ (V_ShapeZ n) -> return $ MV_ShapeZ n+ {-# INLINE basicLength #-}+ basicLength = \(V_ShapeZ n) -> n+ {-# INLINE basicUnsafeSlice #-}+ basicUnsafeSlice = \ _ m (V_ShapeZ _) -> V_ShapeZ m+ {-# INLINE basicUnsafeIndexM #-}+ basicUnsafeIndexM = \ (V_ShapeZ _) _ -> return Nil+ {-# INLINE basicUnsafeCopy #-}+ basicUnsafeCopy = \ (MV_ShapeZ _) (V_ShapeZ _) -> return ()+ {-# INLINE elemseq #-}+ elemseq = \ _ -> seq++instance (UV.Unbox a) => GMV.MVector UV.MVector (Shape ('S 'Z) a) where+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicOverlaps #-}+ {-# INLINE basicUnsafeNew #-}+ {-# INLINE basicUnsafeReplicate #-}+ {-# INLINE basicUnsafeRead #-}+ {-# INLINE basicUnsafeWrite #-}+ {-# INLINE basicClear #-}+ {-# INLINE basicSet #-}+ {-# INLINE basicUnsafeCopy #-}+ {-# INLINE basicUnsafeGrow #-}+ {-# INLINE basicInitialize #-}+ basicInitialize = \ (MV_ShapeSZ v) -> GMV.basicInitialize v+ basicLength = \(MV_ShapeSZ v)-> GMV.basicLength v+ basicUnsafeSlice = \ i n (MV_ShapeSZ v) -> MV_ShapeSZ $ GMV.basicUnsafeSlice i n v+ basicOverlaps = \ (MV_ShapeSZ v1) (MV_ShapeSZ v2) -> GMV.basicOverlaps v1 v2+ basicUnsafeNew = \ n -> MV_ShapeSZ `liftM` GMV.basicUnsafeNew n+ basicUnsafeReplicate = \ n (a:*_) -> MV_ShapeSZ `liftM` GMV.basicUnsafeReplicate n a+ basicUnsafeRead = \ (MV_ShapeSZ v) i -> ( :* Nil ) `liftM` GMV.basicUnsafeRead v i+ basicUnsafeWrite = \ (MV_ShapeSZ v) i (a:* _) -> GMV.basicUnsafeWrite v i a+ basicClear = \ (MV_ShapeSZ v) -> GMV.basicClear v+ basicSet = \ (MV_ShapeSZ v) (a:*_) -> GMV.basicSet v a+ basicUnsafeCopy = \ (MV_ShapeSZ v1) (MV_ShapeSZ v2) -> GMV.basicUnsafeCopy v1 v2+ basicUnsafeMove = \ (MV_ShapeSZ v1) (MV_ShapeSZ v2) -> GMV.basicUnsafeMove v1 v2+ basicUnsafeGrow = \ (MV_ShapeSZ v) n -> MV_ShapeSZ `liftM` GMV.basicUnsafeGrow v n++instance ( UV.Unbox a) => GV.Vector UV.Vector (Shape ('S 'Z) a ) where+ {-# INLINE basicUnsafeFreeze #-}+ {-# INLINE basicUnsafeThaw #-}+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicUnsafeIndexM #-}+ {-# INLINE elemseq #-}+ basicUnsafeFreeze = \ (MV_ShapeSZ v) -> V_ShapeSZ `liftM` GV.basicUnsafeFreeze v+ basicUnsafeThaw = \ (V_ShapeSZ v) -> MV_ShapeSZ`liftM` GV.basicUnsafeThaw v+ basicLength = \ (V_ShapeSZ v)-> GV.basicLength v+ basicUnsafeSlice = \ i n (V_ShapeSZ v) -> V_ShapeSZ $ GV.basicUnsafeSlice i n v+ basicUnsafeIndexM = \ (V_ShapeSZ v) i -> ( :* Nil ) `liftM` GV.basicUnsafeIndexM v i+ basicUnsafeCopy = \ (MV_ShapeSZ mv) (V_ShapeSZ v) -> GV.basicUnsafeCopy mv v+ elemseq = \ _ (a:*_) z -> GV.elemseq (undefined :: UV.Vector a) a z+++instance (UV.Unbox a,UV.Unbox (Shape ('S n) a)) => GMV.MVector UV.MVector (Shape ('S ('S n)) a) where+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicOverlaps #-}+ {-# INLINE basicUnsafeNew #-}+ {-# INLINE basicUnsafeReplicate #-}+ {-# INLINE basicUnsafeRead #-}+ {-# INLINE basicUnsafeWrite #-}+ {-# INLINE basicClear #-}+ {-# INLINE basicSet #-}+ {-# INLINE basicUnsafeCopy #-}+ {-# INLINE basicUnsafeGrow #-}+ {-# INLINE basicInitialize #-}+ basicInitialize = \ (MV_ShapeSSN v) -> GMV.basicInitialize v+ basicLength = \ (MV_ShapeSSN v) -> GMV.basicLength v+ basicUnsafeSlice = \ i n (MV_ShapeSSN v) -> MV_ShapeSSN $ GMV.basicUnsafeSlice i n v+ basicOverlaps = \ (MV_ShapeSSN v1) (MV_ShapeSSN v2) -> GMV.basicOverlaps v1 v2+ basicUnsafeNew = \ n -> MV_ShapeSSN `liftM` GMV.basicUnsafeNew n+ basicUnsafeReplicate = \ n (a :* as) -> MV_ShapeSSN `liftM` GMV.basicUnsafeReplicate n (a,as)+ basicUnsafeRead = \ (MV_ShapeSSN v) i -> uncurry (:*) `liftM` GMV.basicUnsafeRead v i+ basicUnsafeWrite = \(MV_ShapeSSN v) i (a :* as ) -> GMV.basicUnsafeWrite v i (a,as)+ basicClear = \ (MV_ShapeSSN v) -> GMV.basicClear v+ basicSet = \ (MV_ShapeSSN v) (a :* as) -> GMV.basicSet v (a,as)+ basicUnsafeCopy = \ (MV_ShapeSSN v1) (MV_ShapeSSN v2) -> GMV.basicUnsafeCopy v1 v2+ basicUnsafeMove = \ (MV_ShapeSSN v1) (MV_ShapeSSN v2) -> GMV.basicUnsafeMove v1 v2+ basicUnsafeGrow = \ (MV_ShapeSSN v) n -> MV_ShapeSSN `liftM` GMV.basicUnsafeGrow v n+++instance (UV.Unbox a,UV.Unbox (Shape ('S n) a)) => GV.Vector UV.Vector (Shape ('S ('S n)) a) where+ {-# INLINE basicUnsafeFreeze #-}+ {-# INLINE basicUnsafeThaw #-}+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicUnsafeIndexM #-}+ {-# INLINE elemseq #-}+ basicUnsafeFreeze = \ (MV_ShapeSSN v) -> V_ShapeSSN `liftM` GV.basicUnsafeFreeze v+ basicUnsafeThaw = \ (V_ShapeSSN v) -> MV_ShapeSSN `liftM` GV.basicUnsafeThaw v+ basicLength = \ (V_ShapeSSN v) -> GV.basicLength v+ basicUnsafeSlice = \ i n (V_ShapeSSN v) -> V_ShapeSSN $ GV.basicUnsafeSlice i n v+ basicUnsafeIndexM = \ (V_ShapeSSN v) i -> uncurry (:*) `liftM` GV.basicUnsafeIndexM v i+ basicUnsafeCopy = \ (MV_ShapeSSN mv) (V_ShapeSSN v) -> GV.basicUnsafeCopy mv v+ elemseq = \ _ (a :* as) z -> GV.elemseq (undefined :: UV.Vector a) a+ $ GV.elemseq (undefined :: UV.Vector (Shape ('S n) a)) as z
+ src/Numerical/Array/Storage.hs view
@@ -0,0 +1,267 @@++{-# LANGUAGE TypeFamilies,FlexibleInstances,MultiParamTypeClasses,FlexibleContexts #-}+{-# LANGUAGE UndecidableInstances,StandaloneDeriving, DeriveDataTypeable #-}+{-# LANGUAGE DeriveFunctor, DeriveFoldable, DeriveTraversable, DeriveGeneric #-}+{-# LANGUAGE CPP #-}+module Numerical.Array.Storage(+ Boxed+ ,Unboxed+ ,Stored+ ,BufferPure(..)+ ,BufferMut(..)+ ,Buffer+ ,MBuffer+ ,unsafeBufferThaw+ ,unsafeBufferFreeze) where+++import Control.Monad.Primitive ( PrimMonad, PrimState )++import qualified Data.Vector.Generic as VG+import qualified Data.Vector.Generic.Mutable as VGM+import qualified Data.Vector as BV+import qualified Data.Vector.Storable as SV+import qualified Data.Vector.Unboxed as UV++--import qualified Data.Functor as F hiding (Functor)+--import qualified Data.Foldable as F hiding (Foldable)+--import qualified Data.Traversable as T hiding (Traversable)+#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ < 709+--import Data.Functor (Functor)+import Data.Foldable (Foldable)+import Data.Traversable (Traversable)+#endif++import Data.Typeable+import Data.Data+import GHC.Generics+++{-+FIXME : should i require that the element type and+mode are both instance of Typeable for Buffers?++-}+++{-+FIX MEEEEE REMINDERS+make the allocators for Storable Buffers do AVX sized alignment+-}++-- | The class instance @'Buffer' mode a@ is a shorthand for saying that a given buffer representation @mode@+-- has a 'VG.Vector' instance for both 'BufferPure' and 'BufferMut'.+class (VG.Vector (BufferPure mode) a, VGM.MVector (BufferMut mode) a)=> Buffer mode a++instance (VG.Vector (BufferPure mode) a, VGM.MVector (BufferMut mode) a)=> Buffer mode a++-- not sure if MBuffer class should exist, fixme. if/when removed, this+class VGM.MVector (BufferMut mode) a=> MBuffer mode a++-- not sure if MBuffer should exist, FIXME+instance VGM.MVector (BufferMut mode) a=> MBuffer mode a++-- | 'Boxed' is the type index for `Buffer`s that use the boxed data structure `Data.Vector.Vector`+-- as the underlying storage representation.+data Boxed+ deriving Typeable++deriving instance Data Boxed++-- | 'Unboxed' is the type index for 'Buffer's that use the unboxed data structure+-- 'Data.Vector.Unboxed.Vector' as the underlying storage representation.+data Unboxed+ deriving Typeable++deriving instance Data Unboxed++-- | 'Stored' is the type index for 'Buffer's that use the 'Foreign.Storable'+-- for values, in pinned byte array buffers, provided by 'Data.Vector.Storable'+data Stored+ deriving Typeable++deriving instance Data Stored++type instance VG.Mutable (BufferPure sort) = BufferMut sort+++data family BufferPure sort elem++deriving instance Typeable BufferPure++newtype instance BufferPure Boxed elem = BoxedBuffer (BV.Vector elem)+ deriving (Show,Data,Generic,Functor,Foldable,Traversable)++++newtype instance BufferPure Unboxed elem = UnboxedBuffer (UV.Vector elem)+ deriving (Show,Data,Generic)+--deriving instance Typeable a => Typeable (BufferPure Unboxed a)++newtype instance BufferPure Stored elem = StorableBuffer (SV.Vector elem)+ deriving (Show,Data,Generic)++data family BufferMut sort st elem+deriving instance Typeable BufferMut+++newtype instance BufferMut Boxed st elem = BoxedBufferMut (BV.MVector st elem)+ --deriving (Show,Data,Generic)+newtype instance BufferMut Unboxed st elem = UnboxedBufferMut (UV.MVector st elem)+ --deriving (Show,Data,Generic)+newtype instance BufferMut Stored st elem = StorableBufferMut (SV.MVector st elem)++-- | 'unsafeBufferFreeze'+unsafeBufferFreeze :: (Buffer rep a,PrimMonad m) => BufferMut rep (PrimState m ) a -> m (BufferPure rep a)+unsafeBufferFreeze = VG.basicUnsafeFreeze++unsafeBufferThaw :: (Buffer rep a,PrimMonad m) => (BufferPure rep a) -> m (BufferMut rep (PrimState m ) a)+unsafeBufferThaw = VG.basicUnsafeThaw++instance (VGM.MVector BV.MVector elem) => VGM.MVector (BufferMut Boxed) elem where+ basicInitialize = \(BoxedBufferMut v) -> VGM.basicInitialize v+ basicLength = \(BoxedBufferMut v) -> VGM.basicLength v+ basicUnsafeSlice =+ \ ix1 ix2 (BoxedBufferMut bv) ->+ BoxedBufferMut $ VGM.basicUnsafeSlice ix1 ix2 bv+ basicOverlaps =+ \ (BoxedBufferMut bv1) (BoxedBufferMut bv2) -> VGM.basicOverlaps bv1 bv2+ basicUnsafeNew = \ size ->+ do+ res<- VGM.basicUnsafeNew size+ return (BoxedBufferMut res)+ basicUnsafeRead= \(BoxedBufferMut bv) ix -> VGM.basicUnsafeRead bv ix+ basicUnsafeWrite = \(BoxedBufferMut bv ) ix val -> VGM.basicUnsafeWrite bv ix val++ {-Q/todo/check fixme, do these other operations need be provided in a pass through way too?+ or will there be no difference in the derived code perf ? -}+-- basicUnsafeClear+-- basicUnsafeSet+-- basicUnsafeCopy+-- basicUnsafeMove+-- basicUnsafeGrow+-- basicUnsafeReplicate+ {-# INLINE basicInitialize #-}+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicOverlaps #-}+ {-# INLINE basicUnsafeNew #-}+ {-# INLINE basicUnsafeRead #-}+ {-# INLINE basicUnsafeWrite #-}++-- {-# INLINE basicUnsafeClear#-}+-- {-# INLINE basicUnsafeSet#-}+-- {-# INLINE basicUnsafeCopy#-}+-- {-# INLINE basicUnsafeMove#-}+-- {-# INLINE basicUnsafeGrow#-}+-- {-# INLINE basicUnsafeReplicate#-}++instance (SV.Storable elem) => VGM.MVector (BufferMut Stored) elem where+ basicInitialize = \(StorableBufferMut v) -> VGM.basicInitialize v+ basicLength = \(StorableBufferMut v) -> VGM.basicLength v+ basicUnsafeSlice =+ \ ix1 ix2 (StorableBufferMut bv) ->+ StorableBufferMut $ VGM.basicUnsafeSlice ix1 ix2 bv+ basicOverlaps =+ \ (StorableBufferMut bv1) (StorableBufferMut bv2) -> VGM.basicOverlaps bv1 bv2+ basicUnsafeNew = \ size ->+ do+ res<- VGM.basicUnsafeNew size+ return (StorableBufferMut res)+ basicUnsafeRead= \(StorableBufferMut bv) ix -> VGM.basicUnsafeRead bv ix+ basicUnsafeWrite = \(StorableBufferMut bv ) ix val -> VGM.basicUnsafeWrite bv ix val+ {-# INLINE basicInitialize #-}+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicOverlaps #-}+ {-# INLINE basicUnsafeNew #-}+ {-# INLINE basicUnsafeRead #-}+ {-# INLINE basicUnsafeWrite #-}+++instance (VGM.MVector UV.MVector elem) => VGM.MVector (BufferMut Unboxed) elem where+ {-# INLINE basicInitialize #-}+ basicInitialize = \(UnboxedBufferMut v) -> VGM.basicInitialize v+ basicLength = \(UnboxedBufferMut v) -> VGM.basicLength v+ basicUnsafeSlice =+ \ ix1 ix2 (UnboxedBufferMut bv) ->+ UnboxedBufferMut $ VGM.basicUnsafeSlice ix1 ix2 bv+ basicOverlaps =+ \ (UnboxedBufferMut bv1) (UnboxedBufferMut bv2) -> VGM.basicOverlaps bv1 bv2+ basicUnsafeNew = \ size ->+ do+ res<- VGM.basicUnsafeNew size+ return (UnboxedBufferMut res)+ basicUnsafeRead= \(UnboxedBufferMut bv) ix -> VGM.basicUnsafeRead bv ix+ basicUnsafeWrite = \(UnboxedBufferMut bv ) ix val -> VGM.basicUnsafeWrite bv ix val++ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicOverlaps #-}+ {-# INLINE basicUnsafeNew #-}+ {-# INLINE basicUnsafeRead #-}+ {-# INLINE basicUnsafeWrite #-}++----+----+instance VG.Vector BV.Vector a => VG.Vector (BufferPure Boxed) a where++ basicUnsafeFreeze =+ \(BoxedBufferMut mv) ->(\ x->return $ BoxedBuffer x) =<< VG.basicUnsafeFreeze mv+ basicUnsafeThaw= \(BoxedBuffer v) ->(\x -> return $ BoxedBufferMut x ) =<< VG.basicUnsafeThaw v+ basicLength = \(BoxedBuffer v) -> VG.basicLength v+ basicUnsafeSlice =+ \ start len (BoxedBuffer v) -> BoxedBuffer $! VG.basicUnsafeSlice start len v+ basicUnsafeIndexM =+ \ (BoxedBuffer v) ix -> VG.basicUnsafeIndexM v ix+ elemseq = \ (BoxedBuffer v) a b -> VG.elemseq v a b+++ {-# INLINE basicUnsafeFreeze #-}+ {-# INLINE basicUnsafeThaw #-}+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicUnsafeIndexM #-}+ {-# INLINE elemseq #-}+++instance (SV.Storable a) => VG.Vector (BufferPure Stored) a where++ basicUnsafeFreeze =+ \(StorableBufferMut mv) -> (\x ->return $StorableBuffer x) =<< VG.basicUnsafeFreeze mv+ basicUnsafeThaw=+ \(StorableBuffer v) -> (\x -> return $ StorableBufferMut x) =<< VG.basicUnsafeThaw v+ basicLength = \(StorableBuffer v) -> VG.basicLength v+ basicUnsafeSlice =+ \ start len (StorableBuffer v) -> StorableBuffer $! VG.basicUnsafeSlice start len v+ basicUnsafeIndexM =+ \ (StorableBuffer v) ix -> VG.basicUnsafeIndexM v ix+ elemseq = \ (StorableBuffer v) a b -> VG.elemseq v a b+++ {-# INLINE basicUnsafeFreeze #-}+ {-# INLINE basicUnsafeThaw #-}+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicUnsafeIndexM #-}+ {-# INLINE elemseq #-}+++instance VG.Vector UV.Vector a => VG.Vector (BufferPure Unboxed) a where++ basicUnsafeFreeze = \(UnboxedBufferMut mv) -> (\x -> return $ UnboxedBuffer x) =<< VG.basicUnsafeFreeze mv+ basicUnsafeThaw= \(UnboxedBuffer v) ->(\x -> return $ UnboxedBufferMut x) =<< VG.basicUnsafeThaw v+ basicLength = \(UnboxedBuffer v) -> VG.basicLength v+ basicUnsafeSlice =+ \ start len (UnboxedBuffer v) -> UnboxedBuffer $! VG.basicUnsafeSlice start len v+ basicUnsafeIndexM =+ \ (UnboxedBuffer v) ix -> VG.basicUnsafeIndexM v ix+ elemseq = \ (UnboxedBuffer v) a b -> VG.elemseq v a b+++ {-# INLINE basicUnsafeFreeze #-}+ {-# INLINE basicUnsafeThaw #-}+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicUnsafeIndexM #-}+ {-# INLINE elemseq #-}
+ src/Numerical/Data/Vector/HPair.hs view
@@ -0,0 +1,235 @@+{- | This module is pretty cool because it gives you a way to talk about+heterogeneous representations for different columns!++might be replaced with an HList of Vectors approach+-}+++{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeFamilyDependencies#-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE PolyKinds #-}+{-# LANGUAGE MultiParamTypeClasses ,FlexibleInstances , FlexibleContexts,UndecidableInstances #-}++module Numerical.Data.Vector.HPair(+ VHProd(..)+ ,vHPair+ ,vUnHPair+ ,MVHProd(..)+ ,HProd(..)+ ,MutableHProdTree+ ,TransformHProdTree+ --,mvUnPair+ --,mvPair+ ) where++import qualified Data.Vector.Generic as V+import qualified Data.Vector.Generic.Mutable as MV++import Control.Monad.Primitive (PrimMonad)++++--type instance V.Mutable (VPair v) = MVPair (V.Mutable v)+++{-+currently primmonad doesn't get its free applicative/functor powers :*(++-}++(<$$$>) :: PrimMonad m => (a->b) -> m a -> m b+(<$$$>) f mv = do v <- mv ; return (f v )+{-# INLINE (<$$$>) #-}++(<***>) :: PrimMonad m => m (a->b) -> m a -> m b+(<***>) mf mv = do f <- mf ; v <- mv ; return (f v)+{-# INLINE (<***>) #-}++{-+probably should just++-}++{-+the names are terrible, fix them later!+HProd , HPair, HUnit, VHPro+-}+++data HProd a where+ HPair :: HProd a-> HProd a -> HProd a+ HUnit :: a -> HProd a++data VHProd (prd:: HProd ( * -> * )) val where+ VHLeaf :: !(v a) -> VHProd ('HUnit v) a+ VHNode :: !(VHProd pra a) -> !(VHProd prb b ) ->VHProd ('HPair pra prb) (a,b)++data MVHProd (prd:: HProd (* -> * -> *) ) (st :: * ) val where+ MVHLeaf :: !(mv st a) -> MVHProd ('HUnit mv) st a+ MVHNode :: !(MVHProd pra st a) -> !(MVHProd prb st b ) -> MVHProd ('HPair pra prb) st (a,b)+++vHPair :: (va a,vb b)->VHProd ('HPair ('HUnit va) ('HUnit vb)) (a,b)+vHPair = \ (va,vb) -> VHNode (VHLeaf va) (VHLeaf vb)+{-# INLINE vHPair #-}++vUnHPair :: VHProd ('HPair ('HUnit va) ('HUnit vb)) (a,b) -> (va a, vb b)+vUnHPair = \ (VHNode (VHLeaf va) (VHLeaf vb))-> (va,vb)+{-# INLINE vUnHPair #-}++type instance V.Mutable (VHProd prod)= MVHProd (MutableHProdTree prod)++type family MutableHProdTree (a :: HProd (* -> *)) = r | r -> a where+ MutableHProdTree ('HUnit v ) = 'HUnit (V.Mutable v)+ MutableHProdTree ('HPair left right) = 'HPair (MutableHProdTree left) (MutableHProdTree right )++type family TransformHProdTree (f :: k-> m) (a :: HProd k) :: HProd m where+ TransformHProdTree f ('HUnit v)= 'HUnit (f v)+ TransformHProdTree f ('HPair left right) = 'HPair (TransformHProdTree f left) (TransformHProdTree f right)+++++--mvPair :: (mv st a,mv st b)->MVPair mv st (a,b)+--mvPair = \ (mva, mvb) -> TheMVPair mva mvb+--{-# INLINE mvPair #-}++--mvUnPair :: MVPair mv st (a,b) -> (mv st a,mv st b)+--mvUnPair = \ (TheMVPair mva mvb)-> (mva,mvb)+--{-# INLINE mvUnPair #-}++instance (MV.MVector (MVHProd (MutableHProdTree ('HPair pa pb )) ) (a,b) ,+ V.Vector (VHProd pa) a, V.Vector (VHProd pb) b)+ => V.Vector (VHProd ('HPair pa pb )) (a,b) where+ {-# INLINE basicUnsafeFreeze #-}+ {-# INLINE basicUnsafeThaw #-}+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicUnsafeIndexM #-}+++ basicUnsafeFreeze = \(MVHNode mva mvb) ->+ VHNode <$$$> V.basicUnsafeFreeze mva <***> V.basicUnsafeFreeze mvb+++ basicUnsafeThaw = \(VHNode va vb) ->+ MVHNode <$$$> V.basicUnsafeThaw va <***> V.basicUnsafeThaw vb+++ basicLength = \(VHNode va _) -> V.basicLength va+++ basicUnsafeSlice = \start len (VHNode va vb) ->+ VHNode (V.basicUnsafeSlice start len va) (V.basicUnsafeSlice start len vb)+++ basicUnsafeIndexM = \(VHNode va vb) ix ->+ do+ a <- V.basicUnsafeIndexM va ix+ b <- V.basicUnsafeIndexM vb ix+ return (a,b)++instance (MV.MVector (MVHProd ('HUnit (V.Mutable v)) ) a ,V.Vector v a)+ => V.Vector (VHProd ('HUnit v)) a where++ {-# INLINE basicUnsafeFreeze #-}+ {-# INLINE basicUnsafeThaw #-}+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicUnsafeIndexM #-}++ basicUnsafeFreeze = \(MVHLeaf mva) ->+ VHLeaf <$$$> V.basicUnsafeFreeze mva+ basicUnsafeThaw = \(VHLeaf va ) ->+ MVHLeaf <$$$> V.basicUnsafeThaw va+ basicLength = \(VHLeaf va ) -> V.basicLength va+ basicUnsafeSlice = \start len (VHLeaf va ) ->+ VHLeaf(V.basicUnsafeSlice start len va)+ basicUnsafeIndexM = \(VHLeaf va) ix -> V.basicUnsafeIndexM va ix+++instance (MV.MVector mv a) => MV.MVector (MVHProd ('HUnit mv )) a where+ basicLength = \ (MVHLeaf mva) -> MV.basicLength mva+ {-# INLINE basicLength #-}++ basicUnsafeSlice = \ start len (MVHLeaf mva )->+ MVHLeaf (MV.basicUnsafeSlice start len mva)+ {-# INLINE basicUnsafeSlice #-}++ basicOverlaps = \ (MVHLeaf mva ) (MVHLeaf mva2 )-> (MV.basicOverlaps mva mva2)+ {-# INLINE basicOverlaps #-}++ basicUnsafeNew =+ \ size ->+ MVHLeaf <$$$> MV.basicUnsafeNew size+ {-# INLINE basicUnsafeNew #-}++ basicUnsafeReplicate =+ \ size a ->+ MVHLeaf <$$$>+ MV.basicUnsafeReplicate size a+ {-# INLINE basicUnsafeReplicate #-}++ basicUnsafeRead = \(MVHLeaf mva ) ix -> MV.basicUnsafeRead mva ix+ {-#INLINE basicUnsafeRead #-}++ basicUnsafeWrite = \ (MVHLeaf mva ) ix a ->+ do+ MV.basicUnsafeWrite mva ix a+ return ()+ {-#INLINE basicUnsafeWrite #-}++ {-#INLINE basicUnsafeGrow #-}+ basicUnsafeGrow = \ (MVHLeaf mva ) growth ->+ MVHLeaf <$$$> MV.basicUnsafeGrow mva growth++++instance (MV.MVector (MVHProd pra) a,MV.MVector (MVHProd prb) b)+ => MV.MVector (MVHProd ('HPair pra prb)) (a,b) where++ basicLength = \ (MVHNode mva _) -> MV.basicLength mva+ {-# INLINE basicLength #-}++ basicUnsafeSlice = \ start len (MVHNode mva mvb )->+ MVHNode (MV.basicUnsafeSlice start len mva) (MV.basicUnsafeSlice start len mvb)+ {-# INLINE basicUnsafeSlice #-}++ basicOverlaps = \ (MVHNode mva mvb) (MVHNode mva2 mvb2)-> (MV.basicOverlaps mva mva2) || (MV.basicOverlaps mvb mvb2)+ {-# INLINE basicOverlaps #-}++ basicUnsafeNew =+ \ size ->+ MVHNode <$$$> MV.basicUnsafeNew size <***> MV.basicUnsafeNew size+ {-# INLINE basicUnsafeNew #-}++ basicUnsafeReplicate =+ \ size (a,b) ->+ MVHNode <$$$>+ MV.basicUnsafeReplicate size a <***>+ MV.basicUnsafeReplicate size b+ {-# INLINE basicUnsafeReplicate #-}++ basicUnsafeRead = \(MVHNode mva mvb) ix ->+ (,) <$$$> MV.basicUnsafeRead mva ix <***> MV.basicUnsafeRead mvb ix++ {-#INLINE basicUnsafeRead #-}++ basicUnsafeWrite = \ (MVHNode mva mvb) ix (a,b) ->+ do+ MV.basicUnsafeWrite mva ix a+ MV.basicUnsafeWrite mvb ix b+ return ()+ {-#INLINE basicUnsafeWrite #-}++ {-#INLINE basicUnsafeGrow #-}+ basicUnsafeGrow = \ (MVHNode mva mvb) growth ->+ MVHNode <$$$> MV.basicUnsafeGrow mva growth <***>+ MV.basicUnsafeGrow mvb growth+++++
+ src/Numerical/Data/Vector/Pair.hs view
@@ -0,0 +1,205 @@++{- | This module is pretty cool because it gives you a way to talk about+open struct of arrays style vectors++might be replaced with an HList of Vectors approach+++-}++++{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE MultiParamTypeClasses ,FlexibleInstances , FlexibleContexts,UndecidableInstances #-}++module Numerical.Data.Vector.Pair(+ VProd(..)+ ,vPair+ ,vUnPair+ ,MVProd(..)+ --,mvUnPair+ ,Prod(..)+ --,mvPair+ ) where++import qualified Data.Vector.Generic as V+import qualified Data.Vector.Generic.Mutable as MV+++--type instance V.Mutable (VPair v) = MVPair (V.Mutable v)+++++data Prod = Pair Prod Prod | Unit+++data family VProd (vect :: * -> * ) (prd:: Prod ) val -- where+data instance VProd v 'Unit a where+ VLeaf :: !(v a) -> VProd v 'Unit a++data instance VProd v ('Pair pra prb ) (a,b) where+ VPair :: !(VProd v pra a) -> !(VProd v prb b ) ->VProd v ('Pair pra prb) (a,b)++data family MVProd (vect :: * -> * -> * ) (prd:: Prod ) (st :: * ) val -- where+data instance MVProd mv 'Unit st a where+ MVLeaf :: !(mv st a) -> MVProd mv 'Unit st a+data instance MVProd mv ('Pair pra prb) st (a,b) where+ MVPair :: !(MVProd mv pra st a) -> !(MVProd mv prb st b ) -> MVProd mv ('Pair pra prb) st (a,b)+++vPair :: (v a,v b)->VProd v ('Pair 'Unit 'Unit) (a,b)+vPair = \ (va,vb) -> VPair (VLeaf va) (VLeaf vb)+{-# INLINE vPair #-}++vUnPair :: VProd v ('Pair 'Unit 'Unit) (a,b) -> (v a, v b)+vUnPair = \ (VPair (VLeaf va) (VLeaf vb))-> (va,vb)+{-# INLINE vUnPair #-}++type instance V.Mutable (VProd vec prod)= MVProd (V.Mutable vec) prod+++--mvPair :: (mv st a,mv st b)->MVPair mv st (a,b)+--mvPair = \ (mva, mvb) -> TheMVPair mva mvb+--{-# INLINE mvPair #-}++--mvUnPair :: MVPair mv st (a,b) -> (mv st a,mv st b)+--mvUnPair = \ (TheMVPair mva mvb)-> (mva,mvb)+--{-# INLINE mvUnPair #-}++instance (MV.MVector (MVProd (V.Mutable v) ('Pair pa pb ) ) (a,b) ,V.Vector (VProd v pa) a,V.Vector (VProd v pb) b)+ => V.Vector (VProd v ('Pair pa pb )) (a,b) where+ {-# INLINE basicUnsafeFreeze #-}+ basicUnsafeFreeze = \(MVPair mva mvb) ->+ VPair <$> V.basicUnsafeFreeze mva <*> V.basicUnsafeFreeze mvb++ {-# INLINE basicUnsafeThaw #-}+ basicUnsafeThaw = \(VPair va vb) ->+ MVPair <$> V.basicUnsafeThaw va <*> V.basicUnsafeThaw vb++ {-# INLINE basicLength #-}+ basicLength = \(VPair va _) -> V.basicLength va++ {-# INLINE basicUnsafeSlice #-}+ basicUnsafeSlice = \start len (VPair va vb) ->+ VPair (V.basicUnsafeSlice start len va) (V.basicUnsafeSlice start len vb)++ {-# INLINE basicUnsafeIndexM #-}+ basicUnsafeIndexM = \(VPair va vb) ix ->+ do+ a <- V.basicUnsafeIndexM va ix+ b <- V.basicUnsafeIndexM vb ix+ return (a,b)++instance (MV.MVector (MVProd (V.Mutable v) 'Unit ) a ,V.Vector v a)+ => V.Vector (VProd v 'Unit) a where++ {-# INLINE basicUnsafeFreeze #-}+ {-# INLINE basicUnsafeThaw #-}+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicUnsafeIndexM #-}++ basicUnsafeFreeze = \(MVLeaf mva) ->+ VLeaf <$> V.basicUnsafeFreeze mva+ basicUnsafeThaw = \(VLeaf va ) ->+ MVLeaf <$> V.basicUnsafeThaw va+ basicLength = \(VLeaf va ) -> V.basicLength va+ basicUnsafeSlice = \start len (VLeaf va ) ->+ VLeaf(V.basicUnsafeSlice start len va)+ basicUnsafeIndexM = \(VLeaf va) ix -> V.basicUnsafeIndexM va ix+++instance (MV.MVector mv a) => MV.MVector (MVProd mv 'Unit) a where+ basicLength = \ (MVLeaf mva) -> MV.basicLength mva+ {-# INLINE basicLength #-}+++ basicInitialize = \ (MVLeaf mva) -> MV.basicInitialize mva+ {-# INLINE basicInitialize #-}++ basicUnsafeSlice = \ start len (MVLeaf mva )->+ MVLeaf (MV.basicUnsafeSlice start len mva)+ {-# INLINE basicUnsafeSlice #-}++ basicOverlaps = \ (MVLeaf mva ) (MVLeaf mva2 )-> (MV.basicOverlaps mva mva2)+ {-# INLINE basicOverlaps #-}++ basicUnsafeNew =+ \ size ->+ MVLeaf <$> MV.basicUnsafeNew size+ {-# INLINE basicUnsafeNew #-}++ basicUnsafeReplicate =+ \ size a ->+ MVLeaf <$>+ MV.basicUnsafeReplicate size a+ {-# INLINE basicUnsafeReplicate #-}++ basicUnsafeRead = \(MVLeaf mva ) ix -> MV.basicUnsafeRead mva ix+ {-#INLINE basicUnsafeRead #-}++ basicUnsafeWrite = \ (MVLeaf mva ) ix a ->+ do+ MV.basicUnsafeWrite mva ix a+ return ()+ {-#INLINE basicUnsafeWrite #-}++ {-#INLINE basicUnsafeGrow #-}+ basicUnsafeGrow = \ (MVLeaf mva ) growth ->+ MVLeaf <$> MV.basicUnsafeGrow mva growth++++instance (MV.MVector (MVProd mv pra) a,MV.MVector (MVProd mv prb) b) => MV.MVector (MVProd mv ('Pair pra prb)) (a,b) where+ basicLength = \ (MVPair mva _) -> MV.basicLength mva+ {-# INLINE basicLength #-}+++ basicInitialize = \ (MVPair mva mvb) ->+ do MV.basicInitialize mva ;+ MV.basicInitialize mvb+ {-# INLINE basicInitialize #-}++ basicUnsafeSlice = \ start len (MVPair mva mvb )->+ MVPair (MV.basicUnsafeSlice start len mva) (MV.basicUnsafeSlice start len mvb)+ {-# INLINE basicUnsafeSlice #-}++ basicOverlaps = \ (MVPair mva mvb) (MVPair mva2 mvb2)-> (MV.basicOverlaps mva mva2) || (MV.basicOverlaps mvb mvb2)+ {-# INLINE basicOverlaps #-}++ basicUnsafeNew =+ \ size ->+ MVPair <$> MV.basicUnsafeNew size <*> MV.basicUnsafeNew size+ {-# INLINE basicUnsafeNew #-}++ basicUnsafeReplicate =+ \ size (a,b) ->+ MVPair <$>+ MV.basicUnsafeReplicate size a <*>+ MV.basicUnsafeReplicate size b+ {-# INLINE basicUnsafeReplicate #-}++ basicUnsafeRead = \(MVPair mva mvb) ix ->+ (,) <$> MV.basicUnsafeRead mva ix <*> MV.basicUnsafeRead mvb ix++ {-# INLINE basicUnsafeRead #-}++ basicUnsafeWrite = \ (MVPair mva mvb) ix (a,b) ->+ do+ MV.basicUnsafeWrite mva ix a+ MV.basicUnsafeWrite mvb ix b+ return ()+ {-#INLINE basicUnsafeWrite #-}++ {-#INLINE basicUnsafeGrow #-}+ basicUnsafeGrow = \ (MVPair mva mvb) growth ->+ MVPair <$> MV.basicUnsafeGrow mva growth <*>+ MV.basicUnsafeGrow mvb growth+++++
+ src/Numerical/InternalUtils.hs view
@@ -0,0 +1,14 @@+{-# LANGUAGE NoImplicitPrelude#-}+module Numerical.InternalUtils(+ error+) where++--import GHC.Stack (errorWithStackTrace,currentCallStack,whoCreated)+import Prelude (error)++++{-+note well: the stack traces only exist+when doing a profiling build in GHC < 7.9/7.10+-}
+ src/Numerical/Matrix/Basic.hs view
@@ -0,0 +1,32 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE ScopedTypeVariables#-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FunctionalDependencies #-}+++module Numerical.Matrix.Basic where+{-+note, the contents of this module will probably be completely relocated elsewhere+at some point++-}+--import Numerical.Array.Mutable as Mut+++--cleverDotProduct :: (Address ~ MArrayAddress mvecL+-- , Address ~ MArrayAddress mvecR+-- ,Array mvecL (S Z) a+-- ,Array mvecR (S Z) a+-- ,Num a)=>++--naiveDotProduct++
+ src/Numerical/Nat.hs view
@@ -0,0 +1,107 @@+{-# LANGUAGE DataKinds, GADTs, TypeFamilies, TypeOperators,+ ConstraintKinds, ScopedTypeVariables, RankNTypes #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE StandaloneDeriving #-}+{-# LANGUAGE DeriveDataTypeable#-}+{-# LANGUAGE CPP #-}++module Numerical.Nat(Nat(..),N0,N1,N2,N3,N4,N5,N6,N7,N8,N9,N10+ ,SNat(..), type (+),plus_id_r,plus_succ_r,gcastWith,Proxy(..),LitNat,U) where+import Data.Typeable+import Data.Data+import qualified GHC.TypeLits as TL+++#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 707+import Data.Type.Equality(gcastWith)+#else+import Data.Proxy+#endif++type LitNat = TL.Nat++data Nat = S !Nat | Z+ deriving (Eq,Show,Read,Typeable,Data)++#if defined(__GLASGOW_HASKELL__) && ( __GLASGOW_HASKELL__ >= 707) && ( __GLASGOW_HASKELL__ < 709)+deriving instance Typeable 'Z+deriving instance Typeable 'S+#endif++{-+use closed type families when available,+need to test that the+-}+++type family U (n:: TL.Nat) :: Nat where+ U 0 = 'Z+ U n = 'S (U (((TL.-)) n 1))+++type family n1 + n2 where+ 'Z + n2 = n2+ ('S n1') + n2 = 'S (n1' + n2)+++-- ghc 7.6 instances++--type family U (n:: (TL.Nat)) :: Nat++---- can't induct, hence crippled+--type instance U n = Z++--type family (n1::Nat) + (n2::Nat) :: Nat+--type instance Z + n2 = n2+--type instance (S n1) + n2 = S (n1 + n2)+--gcastWith :: (a :~: b) -> ((a ~ b) => r) -> r+--gcastWith Refl x = x+--data a :~: b where+-- Refl :: a :~: a++++++-- singleton for Nat+++data SNat :: Nat -> * where+ SZero :: SNat 'Z+ SSucc :: SNat n -> SNat ('S n)++++-- inductive proof of right-identity of ++plus_id_r :: SNat n -> ((n + 'Z) :~: n)+plus_id_r SZero = Refl+plus_id_r (SSucc n) = gcastWith (plus_id_r n) Refl++-- inductive proof of simplification on the rhs of ++plus_succ_r :: SNat n1 -> Proxy n2 -> ((n1 + ('S n2)) :~: ('S (n1 + n2)))+plus_succ_r SZero _ = Refl+plus_succ_r (SSucc n1) proxy_n2 = gcastWith (plus_succ_r n1 proxy_n2) Refl++++type N0 = 'Z++type N1 = 'S N0++type N2 = 'S N1++type N3 = 'S N2++type N4 = 'S N3++type N5 = 'S N4++type N6 = 'S N5++type N7 = 'S N6++type N8 = 'S N7++type N9 = 'S N8++type N10 = 'S N9
+ src/Numerical/World.hs view
@@ -0,0 +1,18 @@++++module Numerical.World where ++++{-| +Every numerical algorithm runs somewhere.++This could be on a CPU, a GPU, ++-}+-- Native is Just Haskell and Cbits, no external Deps+data Native++-- ForeignNative can have foreign lib deps, +data ForeignNative
+ tests/Main.hs view
@@ -0,0 +1,32 @@+module Main where+++import NumericalUnit.Layout+import NumericalUnit.Shape++++import Data.List+import Data.Ord++import Test.Hspec+import Control.Exception (evaluate)++main :: IO ()+main = hspec $ do+ describe "Shape Unit Tests" $ unitTestShape++--main = defaultMain tests++--tests :: Spec+--tests = testGroup "Unit Tests" [unitTestShape] -- , unitTestLayout ]+++--unitTests = testGroup "Unit tests"+-- [ testCase "List comparison (different length)" $+-- [1, 2, 3] `compare` [1,2] @?= GT++-- -- the following test does not hold+-- , testCase "List comparison (same length)" $+-- [1, 2, 3] `compare` [1,2,2] @?= LT+-- ]
+ tests/NumericalUnit/Layout.hs view
@@ -0,0 +1,7 @@+module NumericalUnit.Layout(unitTestLayout) where ++import Test.HUnit++unitTestLayout = [+ + ]
+ tests/NumericalUnit/Shape.hs view
@@ -0,0 +1,53 @@+{-# LANGUAGE DataKinds, GADTs, TypeFamilies #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FunctionalDependencies #-}++module NumericalUnit.Shape(unitTestShape) where+++import Numerical.Array.Shape as S+import qualified Data.Vector.Storable as SV+import qualified Data.Vector.Unboxed as UV+import Prelude as P+import Test.Hspec.Expectations+import Test.Hspec++unitTestShape :: Spec+unitTestShape = describe "unit tests for Shape" $ do+ specify "foldl on shape" $ S.foldl (+) 0 (1:* 2:* 3 :* Nil ) `shouldBe` P.foldl (+) 0 [1,2,3]+ specify "foldr on shape" $ S.foldr (+) 0 (1:* 2:* 3 :* Nil ) `shouldBe` P.foldr (+) 0 [1,2,3]+ specify "foldl1 on shape" $ S.foldl1 (+) (1:* 2:* 3 :* Nil ) `shouldBe` P.foldl1 (+) [1,2,3]+ specify "foldr1 on shape" $ S.foldr1 (+) (1:* 2:* 3 :* Nil ) `shouldBe` P.foldr1 (+) [1,2,3]++ specify "shapeToList on shape" $ S.shapeToList (1:* 2 :* 3 :* Nil) `shouldBe` [1,2,3]++ specify "Show on Nil shape" $ show Nil `shouldBe` "Nil"+ specify "Show on 1:* Nil" $ show (1:* Nil) `shouldBe` "1 :* Nil"++ specify "storable on size 0 shape" $+ do a <- return (svFromList [Nil,Nil :: Shape Z Int]) ; SV.toList a `shouldBe` [Nil,Nil]+ specify "storable on size 1 shape" $+ do a <- return (svFromList [1:*Nil,2:*Nil :: Shape (S Z) Int]) ; SV.toList a `shouldBe` [1:*Nil,2:*Nil]+ specify "storable on size 2 shape" $+ do a <- return (svFromList [3:* 4:* Nil,1:*2:*Nil :: Shape (S (S Z)) Int]) ;+ SV.toList a `shouldBe` [3:* 4:* Nil,1:*2:*Nil]++ specify "unboxed on size 0 shape" $+ do a <- return (uvFromList [Nil,Nil :: Shape Z Int]) ; UV.toList a `shouldBe` [Nil,Nil]+ specify "unboxed on size 1 shape" $+ do a <- return (uvFromList [1:*Nil,2:*Nil :: Shape (S Z) Int]) ; UV.toList a `shouldBe` [1:*Nil,2:*Nil]+ specify "unboxed on size 2 shape" $+ do a <- return (uvFromList [3:* 4:* Nil,1:*2:*Nil :: Shape (S (S Z)) Int]) ;+ UV.toList a `shouldBe` [3:* 4:* Nil,1:*2:*Nil]++ where+ {- The NOINLINE is need to properly check storable /unboxed instances, otherwise fusion removes the allocation! -}+ svFromList :: SV.Storable a => [a] -> SV.Vector a+ svFromList = SV.fromList+ {-# NOINLINE svFromList #-}++ uvFromList :: UV.Unbox a => [a] -> UV.Vector a+ uvFromList = UV.fromList+ {-# NOINLINE uvFromList#-}