hmatrix 0.17.0.2 → 0.18.0.0

raw patch · 22 files changed

+575/−261 lines, 22 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• CHANGELOG +13/−0
• THANKS.md +25/−1
• hmatrix.cabal +3/−6
• src/Internal/Algorithms.hs +28/−4
• src/Internal/C/vector-aux.c +49/−0
• src/Internal/Container.hs +3/−5
• src/Internal/Devel.hs +2/−3
• src/Internal/Element.hs +0/−3
• src/Internal/Foreign.hs +0/−102
• src/Internal/LAPACK.hs +28/−26
• src/Internal/Matrix.hs +20/−18
• src/Internal/Modular.hs +8/−8
• src/Internal/Numeric.hs +0/−1
• src/Internal/Random.hs +2/−2
• src/Internal/Sparse.hs +2/−3
• src/Internal/Static.hs +73/−14
• src/Internal/Util.hs +6/−4
• src/Internal/Vector.hs +14/−18
• src/Internal/Vectorized.hs +18/−16
• src/Numeric/LinearAlgebra.hs +5/−3
• src/Numeric/LinearAlgebra/Devel.hs +0/−8
• src/Numeric/LinearAlgebra/Static.hs +276/−16

@@ -1,3 +1,16 @@+0.18.0.0+--------++    * Many new functions and instances in the Static module+    +    * meanCov and gaussianSample use Herm type++    * thinQR, thinRQ+    +    * compactSVDTol++    * unitary changed to normalize, also admits Vector (Complex Double)+   0.17.0.0 -------- 

@@ -192,7 +192,8 @@  - Matt Peddie wrote the interfaces to the interpolation and simulated annealing modules. -- "maxc01" solved uninstallability in FreeBSD and improved urandom+- "maxc01" solved uninstallability in FreeBSD, improved urandom, and fixed a Windows+  link error using rand_s.  - "ntfrgl" added {take,drop}Last{Rows,Columns} and odeSolveVWith with generalized step control function    and fixed link errors related to mod/mod_l.@@ -203,4 +204,27 @@  - Vassil Keremidchiev fixed the cabal options for OpenBlas, fixed several installation   issues, and added support for stack-based build.++- Greg Nwosu fixed arm compilation++- Patrik Jansson changed meanCov and gaussianSample to use Herm type. Fixed stack.yaml.++- Justin Le added NFData instances for Static types, added mapping and outer product+  methods to Domain, and many other functions to the Static module.++- Sidharth Kapur added Normed and numeric instances for several Static types,+fixed the CPP issue in cabal files, and made many other contributions.++- Matt Renaud improved the documentation.++- Joshua Moerman fixed cabal/stack flags for windows.++- Francesco Mazzoli, Niklas Hambüchen, Patrick Chilton, and Andras Slemmer+  discovered a serious and subtle bug in the wrapper helpers causing memory corruption.+  Andras Slemmer fixed the bug. Thank you all.++- Kevin Slagle implemented thinQR and thinRQ, much faster than the original qr,+  and added compactSVDTol.++- "fedeinthemix" suggested a better name and a more general type for unitary. 

@@ -1,5 +1,5 @@ Name:               hmatrix-Version:            0.17.0.2+Version:            0.18.0.0 License:            BSD3 License-file:       LICENSE Author:             Alberto Ruiz@@ -55,7 +55,6 @@                         Internal.Devel                         Internal.Vectorized                         Internal.Matrix-                        Internal.Foreign                         Internal.ST                         Internal.IO                         Internal.Element@@ -79,8 +78,7 @@                         src/Internal/C/vector-aux.c  -    extensions:         ForeignFunctionInterface,-                        CPP+    extensions:         ForeignFunctionInterface      ghc-options:        -Wall                         -fno-warn-missing-signatures@@ -94,7 +92,6 @@     if arch(i386)         cc-options:     -msse2 -    cpp-options:        -DBINARY      if os(OSX)         if flag(openblas)@@ -124,7 +121,7 @@      if os(windows)         if flag(openblas)-            extra-libraries:    libopenblas+            extra-libraries:    libopenblas, libgcc_s_seh-1, libgfortran-3, libquadmath-0         else             extra-libraries:    blas lapack 

@@ -30,6 +30,7 @@ import Internal.Numeric import Data.List(foldl1') import qualified Data.Array as A+import qualified Data.Vector.Storable as Vector import Internal.ST import Internal.Vectorized(range) import Control.DeepSeq@@ -291,9 +292,13 @@  -} compactSVD :: Field t  => Matrix t -> (Matrix t, Vector Double, Matrix t)-compactSVD m = (u', subVector 0 d s, v') where+compactSVD = compactSVDTol 1++-- | @compactSVDTol r@ is similar to 'compactSVD' (for which @r=1@), but uses tolerance @tol=r*g*eps*(max rows cols)@ to distinguish nonzero singular values, where @g@ is the greatest singular value. If @g<r*eps@, then only one singular value is returned.+compactSVDTol :: Field t  => Double -> Matrix t -> (Matrix t, Vector Double, Matrix t)+compactSVDTol r m = (u', subVector 0 d s, v') where     (u,s,v) = thinSVD m-    d = rankSVD (1*eps) m s `max` 1+    d = rankSVD (r*eps) m s `max` 1     u' = takeColumns d u     v' = takeColumns d v @@ -475,9 +480,14 @@ -- | QR factorization. -- -- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular.+-- Note: the current implementation is very slow for large matrices. 'thinQR' is much faster. qr :: Field t => Matrix t -> (Matrix t, Matrix t) qr = {-# SCC "qr" #-} unpackQR . qr' +-- | A version of 'qr' which returns only the @min (rows m) (cols m)@ columns of @q@ and rows of @r@.+thinQR :: Field t => Matrix t -> (Matrix t, Matrix t)+thinQR = {-# SCC "thinQR" #-} thinUnpackQR . qr'+ -- | Compute the QR decomposition of a matrix in compact form. qrRaw :: Field t => Matrix t -> QR t qrRaw m = QR x v@@ -494,9 +504,17 @@ -- | RQ factorization. -- -- If @(r,q) = rq m@ then @m == r \<> q@, where q is unitary and r is upper triangular.+-- Note: the current implementation is very slow for large matrices. 'thinRQ' is much faster. rq :: Field t => Matrix t -> (Matrix t, Matrix t)-rq m =  {-# SCC "rq" #-} (r,q) where-    (q',r') = qr $ trans $ rev1 m+rq = {-# SCC "rq" #-} rqFromQR qr++-- | A version of 'rq' which returns only the @min (rows m) (cols m)@ columns of @r@ and rows of @q@.+thinRQ :: Field t => Matrix t -> (Matrix t, Matrix t)+thinRQ = {-# SCC "thinQR" #-} rqFromQR thinQR++rqFromQR :: Field t => (Matrix t -> (Matrix t, Matrix t)) -> Matrix t -> (Matrix t, Matrix t)+rqFromQR qr0 m = (r,q) where+    (q',r') = qr0 $ trans $ rev1 m     r = rev2 (trans r')     q = rev2 (trans q')     rev1 = flipud . fliprl@@ -723,6 +741,12 @@           vs = zipWith zh [1..mn] cs           hs = zipWith haussholder (toList tau) vs           q = foldl1' mXm hs++thinUnpackQR :: (Field t) => (Matrix t, Vector t) -> (Matrix t, Matrix t)+thinUnpackQR (pq, tau) = (q, r)+  where mn = uncurry min $ size pq+        q = qrgr mn $ QR pq tau+        r = fromRows $ zipWith (\i v -> Vector.replicate i 0 Vector.++ Vector.drop i v) [0..mn-1] (toRows pq)  unpackHess :: (Field t) => (Matrix t -> (Matrix t,Vector t)) -> Matrix t -> (Matrix t, Matrix t) unpackHess hf m

@@ -945,6 +945,8 @@ /* Windows use thread-safe random    See: http://stackoverflow.com/questions/143108/is-windows-rand-s-thread-safe */+#if defined (__APPLE__) || (__FreeBSD__)+ inline double urandom() {     /* the probalility of matching will be theoretically p^3(in fact, it is not)        p is matching probalility of random().@@ -959,6 +961,22 @@     return (double)nrand48(state) / (double)max_random; } +#else++#define _CRT_RAND_S+inline double urandom() {+    unsigned int number;+    errno_t err;+    err = rand_s(&number);+    if (err!=0) {+        printf("something wrong\n");+        return -1;+    }+    return (double)number / (double)UINT_MAX;+}++#endif+ double gaussrand(int *phase, double *pV1, double *pV2, double *pS) { 	double V1=*pV1, V2=*pV2, S=*pS;@@ -985,10 +1003,39 @@  } +#if defined(_WIN32) || defined(WIN32)+ int random_vector(unsigned int seed, int code, DVEC(r)) {     int phase = 0;     double V1,V2,S; +    srand(seed);++    int k;+    switch (code) {+      case 0: { // uniform+        for (k=0; k<rn; k++) {+            rp[k] = urandom();+        }+        OK+      }+      case 1: { // gaussian+        for (k=0; k<rn; k++) {+            rp[k] = gaussrand(&phase,&V1,&V2,&S);+        }+        OK+      }++      default: ERROR(BAD_CODE);+    }+}++#else++int random_vector(unsigned int seed, int code, DVEC(r)) {+    int phase = 0;+    double V1,V2,S;+     srandom(seed);      int k;@@ -1009,6 +1056,8 @@       default: ERROR(BAD_CODE);     } }++#endif  #else 

@@ -28,7 +28,7 @@ import Internal.Matrix import Internal.Element import Internal.Numeric-import Internal.Algorithms(Field,linearSolveSVD)+import Internal.Algorithms(Field,linearSolveSVD,Herm,mTm)  ------------------------------------------------------------------ @@ -206,14 +206,14 @@  , -1.0127225830525157e-2,     3.0373954915729318 ])  -}-meanCov :: Matrix Double -> (Vector Double, Matrix Double)+meanCov :: Matrix Double -> (Vector Double, Herm Double) meanCov x = (med,cov) where     r    = rows x     k    = 1 / fromIntegral r     med  = konst k r `vXm` x     meds = konst 1 r `outer` med     xc   = x `sub` meds-    cov  = scale (recip (fromIntegral (r-1))) (trans xc `mXm` xc)+    cov  = scale (recip (fromIntegral (r-1))) (mTm xc)  -------------------------------------------------------------------------------- @@ -293,5 +293,3 @@     | otherwise = error $ "out of range index in remap"   where     [i',j'] = conformMs [i,j]-    -

@@ -80,8 +80,8 @@   where     type Trans c b     type TransRaw c b-    apply      :: (Trans c b) -> c -> b-    applyRaw   :: (TransRaw c b) -> c -> b+    apply      :: c -> (b -> IO r) -> (Trans c b) -> IO r+    applyRaw   :: c -> (b -> IO r) -> (TransRaw c b) -> IO r     infixl 1 `apply`, `applyRaw`  instance Storable t => TransArray (Vector t)@@ -92,4 +92,3 @@     {-# INLINE apply #-}     applyRaw = avec     {-# INLINE applyRaw #-}-

@@ -3,7 +3,6 @@ {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE MultiParamTypeClasses #-}-{-# LANGUAGE CPP #-}  ----------------------------------------------------------------------------- -- |@@ -35,7 +34,6 @@  ------------------------------------------------------------------- -#ifdef BINARY  import Data.Binary @@ -48,7 +46,6 @@           v <- get           return (reshape c v) -#endif  ------------------------------------------------------------------- 

@@ -1,102 +0,0 @@-{-# LANGUAGE MagicHash, UnboxedTuples #-}--- | FFI and hmatrix helpers.------ Sample usage, to upload a perspective matrix to a shader.------ @ glUniformMatrix4fv 0 1 (fromIntegral gl_TRUE) \`appMatrix\` perspective 0.01 100 (pi\/2) (4\/3) --- @-----module Internal.Foreign -    ( app-    , appVector, appVectorLen-    , appMatrix, appMatrixLen, appMatrixRaw, appMatrixRawLen-    , unsafeMatrixToVector, unsafeMatrixToForeignPtr-    ) where--import Foreign.C.Types(CInt)-import Internal.Vector-import Internal.Matrix-import qualified Data.Vector.Storable as S-import Foreign (Ptr, ForeignPtr, Storable)-import GHC.Base (IO(..), realWorld#)--{-# INLINE unsafeInlinePerformIO #-}--- | If we use unsafePerformIO, it may not get inlined, so in a function that returns IO (which are all safe uses of app* in this module), there would be--- unecessary calls to unsafePerformIO or its internals.-unsafeInlinePerformIO :: IO a -> a-unsafeInlinePerformIO (IO f) = case f realWorld# of-    (# _, x #) -> x--{-# INLINE app #-}--- | Only useful since it is left associated with a precedence of 1, unlike 'Prelude.$', which is right associative.--- e.g.------ @--- someFunction---     \`appMatrixLen\` m---     \`appVectorLen\` v---     \`app\` other---     \`app\` arguments---     \`app\` go here--- @------ One could also write:------ @--- (someFunction ---     \`appMatrixLen\` m---     \`appVectorLen\` v) ---     other ---     arguments ---     (go here)--- @----app :: (a -> b) -> a -> b-app f = f--{-# INLINE appVector #-}-appVector :: Storable a => (Ptr a -> b) -> Vector a -> b-appVector f x = unsafeInlinePerformIO (S.unsafeWith x (return . f))--{-# INLINE appVectorLen #-}-appVectorLen :: Storable a => (CInt -> Ptr a -> b) -> Vector a -> b-appVectorLen f x = unsafeInlinePerformIO (S.unsafeWith x (return . f (fromIntegral (S.length x))))--{-# INLINE appMatrix #-}-appMatrix :: Element a => (Ptr a -> b) -> Matrix a -> b-appMatrix f x = unsafeInlinePerformIO (S.unsafeWith (flatten x) (return . f))--{-# INLINE appMatrixLen #-}-appMatrixLen :: Element a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b-appMatrixLen f x = unsafeInlinePerformIO (S.unsafeWith (flatten x) (return . f r c))-  where-    r = fromIntegral (rows x)-    c = fromIntegral (cols x)--{-# INLINE appMatrixRaw #-}-appMatrixRaw :: Storable a => (Ptr a -> b) -> Matrix a -> b-appMatrixRaw f x = unsafeInlinePerformIO (S.unsafeWith (xdat x) (return . f))--{-# INLINE appMatrixRawLen #-}-appMatrixRawLen :: Element a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b-appMatrixRawLen f x = unsafeInlinePerformIO (S.unsafeWith (xdat x) (return . f r c))-  where-    r = fromIntegral (rows x)-    c = fromIntegral (cols x)--infixl 1 `app`-infixl 1 `appVector`-infixl 1 `appMatrix`-infixl 1 `appMatrixRaw`--{-# INLINE unsafeMatrixToVector #-}--- | This will disregard the order of the matrix, and simply return it as-is. --- If the order of the matrix is RowMajor, this function is identical to 'flatten'.-unsafeMatrixToVector :: Matrix a -> Vector a-unsafeMatrixToVector = xdat--{-# INLINE unsafeMatrixToForeignPtr #-}-unsafeMatrixToForeignPtr :: Storable a => Matrix a -> (ForeignPtr a, Int)-unsafeMatrixToForeignPtr m = S.unsafeToForeignPtr0 (xdat m)-

@@ -18,7 +18,7 @@  import Internal.Devel import Internal.Vector-import Internal.Matrix hiding ((#))+import Internal.Matrix hiding ((#), (#!)) import Internal.Conversion import Internal.Element import Foreign.Ptr(nullPtr)@@ -28,10 +28,13 @@  ----------------------------------------------------------------------------------- -infixl 1 #+infixr 1 # a # b = apply a b {-# INLINE (#) #-} +a #! b = a # b # id+{-# INLINE (#!) #-}+ -----------------------------------------------------------------------------------  type TMMM t = t ::> t ::> t ::> Ok@@ -56,7 +59,7 @@     when (cols a /= rows b) $ error $ "inconsistent dimensions in matrix product "++                                        show (rows a,cols a) ++ " x " ++ show (rows b, cols b)     s <- createMatrix ColumnMajor (rows a) (cols b)-    f (isT a) (isT b) # (tt a) # (tt b) # s #| st+    ((tt a) # (tt b) #! s) (f (isT a) (isT b)) #| st     return s  -- | Matrix product based on BLAS's /dgemm/.@@ -80,7 +83,7 @@     when (cols a /= rows b) $ error $         "inconsistent dimensions in matrix product "++ shSize a ++ " x " ++ shSize b     s <- createMatrix ColumnMajor (rows a) (cols b)-    c_multiplyI m # a # b # s #|"c_multiplyI"+    (a # b #! s) (c_multiplyI m) #|"c_multiplyI"     return s  multiplyL :: Z -> Matrix Z -> Matrix Z -> Matrix Z@@ -88,7 +91,7 @@     when (cols a /= rows b) $ error $         "inconsistent dimensions in matrix product "++ shSize a ++ " x " ++ shSize b     s <- createMatrix ColumnMajor (rows a) (cols b)-    c_multiplyL m # a # b # s #|"c_multiplyL"+    (a # b #! s) (c_multiplyL m) #|"c_multiplyL"     return s  -----------------------------------------------------------------------------@@ -121,7 +124,7 @@     u <- createMatrix ColumnMajor r r     s <- createVector (min r c)     v <- createMatrix ColumnMajor c c-    f # a # u # s # v #| st+    (a # u # s #! v) f #| st     return (u,s,v)   where     r = rows x@@ -149,7 +152,7 @@     u <- createMatrix ColumnMajor r q     s <- createVector q     v <- createMatrix ColumnMajor q c-    f # a # u # s # v #| st+    (a # u # s #! v) f #| st     return (u,s,v)   where     r = rows x@@ -176,7 +179,7 @@ svAux f st x = unsafePerformIO $ do     a <- copy ColumnMajor x     s <- createVector q-    g # a # s #| st+    (a #! s) g #| st     return s   where     r = rows x@@ -197,7 +200,7 @@     a <- copy ColumnMajor x     s <- createVector q     v <- createMatrix ColumnMajor c c-    g # a # s # v #| st+    (a # s #! v) g #| st     return (s,v)   where     r = rows x@@ -218,7 +221,7 @@     a <- copy ColumnMajor x     u <- createMatrix ColumnMajor r r     s <- createVector q-    g # a # u # s #| st+    (a # u #! s) g #| st     return (u,s)   where     r = rows x@@ -237,7 +240,7 @@     a <- copy ColumnMajor m     l <- createVector r     v <- createMatrix ColumnMajor r r-    g # a # l # v #| st+    (a # l #! v) g #| st     return (l,v)   where     r = rows m@@ -252,7 +255,7 @@ eigOnlyAux f st m = unsafePerformIO $ do     a <- copy ColumnMajor m     l <- createVector r-    g # a # l #| st+    (a #! l) g #| st     return l   where     r = rows m@@ -277,7 +280,7 @@     a <- copy ColumnMajor m     l <- createVector r     v <- createMatrix ColumnMajor r r-    g # a # l # v #| "eigR"+    (a # l #! v) g #| "eigR"     return (l,v)   where     r = rows m@@ -305,7 +308,7 @@ eigSHAux f st m = unsafePerformIO $ do     l <- createVector r     v <- copy ColumnMajor m-    f # l # v #| st+    (l #! v) f #| st     return (l,v)   where     r = rows m@@ -356,7 +359,7 @@     | n1==n2 && n1==r = unsafePerformIO . g $ do         a' <- copy ColumnMajor a         s  <- copy ColumnMajor b-        f # a' # s #| st+        (a' #! s) f #| st         return s     | otherwise = error $ st ++ " of nonsquare matrix"   where@@ -387,7 +390,7 @@ linearSolveSQAux2 g f st a b     | n1==n2 && n1==r = unsafePerformIO . g $ do         s <- copy ColumnMajor b-        f # a # s #| st+        (a #! s) f #| st         return s     | otherwise = error $ st ++ " of nonsquare matrix"   where@@ -415,7 +418,7 @@         a' <- copy ColumnMajor a         r  <- createMatrix ColumnMajor (max m n) nrhs         setRect 0 0 b r-        f # a' # r #| st+        (a' #! r) f #| st         return r     | otherwise = error $ "different number of rows in linearSolve ("++st++")"   where@@ -458,7 +461,7 @@  cholAux f st a = do     r <- copy ColumnMajor a-    f # r #| st+    (r # id) f #| st     return r  -- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/.@@ -495,7 +498,7 @@ qrAux f st a = unsafePerformIO $ do     r <- copy ColumnMajor a     tau <- createVector mn-    f # tau # r #| st+    (tau #! r) f #| st     return (r,tau)   where     m = rows a@@ -514,7 +517,7 @@  qrgrAux f st n (a, tau) = unsafePerformIO $ do     res <- copy ColumnMajor (subMatrix (0,0) (rows a,n) a)-    f # (subVector 0 n tau') # res #| st+    ((subVector 0 n tau') #! res) f #| st     return res   where     tau' = vjoin [tau, constantD 0 n]@@ -534,7 +537,7 @@ hessAux f st a = unsafePerformIO $ do     r <- copy ColumnMajor a     tau <- createVector (mn-1)-    f # tau # r #| st+    (tau #! r) f #| st     return (r,tau)   where     m = rows a@@ -556,7 +559,7 @@ schurAux f st a = unsafePerformIO $ do     u <- createMatrix ColumnMajor n n     s <- copy ColumnMajor a-    f # u # s #| st+    (u #! s) f #| st     return (u,s)   where     n = rows a@@ -576,7 +579,7 @@ luAux f st a = unsafePerformIO $ do     lu <- copy ColumnMajor a     piv <- createVector (min n m)-    f # piv # lu #| st+    (piv #! lu) f #| st     return (lu, map (pred.round) (toList piv))   where     n = rows a@@ -598,7 +601,7 @@ lusAux f st a piv b     | n1==n2 && n2==n =unsafePerformIO $ do          x <- copy ColumnMajor b-         f # a # piv' # x #| st+         (a # piv' #! x) f #| st          return x     | otherwise = error st   where@@ -622,7 +625,7 @@ ldlAux f st a = unsafePerformIO $ do     ldl <- copy ColumnMajor a     piv <- createVector (rows a)-    f # piv # ldl #| st+    (piv #! ldl) f #| st     return (ldl, map (pred.round) (toList piv))  -----------------------------------------------------------------------------------@@ -637,4 +640,3 @@ -- | Solve a complex linear system from a precomputed LDL decomposition ('ldlC'), using LAPACK's /zsytrs/. ldlsC :: Matrix (Complex Double) -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double) ldlsC a piv b = lusAux zsytrs "ldlsC" (fmat a) piv b-

@@ -5,10 +5,9 @@ {-# LANGUAGE TypeOperators            #-} {-# LANGUAGE TypeFamilies             #-} {-# LANGUAGE ViewPatterns             #-}+{-# LANGUAGE DeriveGeneric            #-} {-# LANGUAGE ConstrainedClassMethods  #-} -- -- | -- Module      :  Internal.Matrix -- Copyright   :  (c) Alberto Ruiz 2007-15@@ -23,7 +22,7 @@  import Internal.Vector import Internal.Devel-import Internal.Vectorized hiding ((#))+import Internal.Vectorized hiding ((#), (#!)) import Foreign.Marshal.Alloc ( free ) import Foreign.Marshal.Array(newArray) import Foreign.Ptr ( Ptr )@@ -111,15 +110,15 @@  -- C-Haskell matrix adapters {-# INLINE amatr #-}-amatr :: Storable a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b-amatr f x = inlinePerformIO (unsafeWith (xdat x) (return . f r c))+amatr :: Storable a => Matrix a -> (f -> IO r) -> (CInt -> CInt -> Ptr a -> f) -> IO r+amatr x f g = unsafeWith (xdat x) (f . g r c)   where     r  = fi (rows x)     c  = fi (cols x)  {-# INLINE amat #-}-amat :: Storable a => (CInt -> CInt -> CInt -> CInt -> Ptr a -> b) -> Matrix a -> b-amat f x = inlinePerformIO (unsafeWith (xdat x) (return . f r c sr sc))+amat :: Storable a => Matrix a -> (f -> IO r) -> (CInt -> CInt -> CInt -> CInt -> Ptr a -> f) -> IO r+amat x f g = unsafeWith (xdat x) (f . g r c sr sc)   where     r  = fi (rows x)     c  = fi (cols x)@@ -136,10 +135,13 @@     applyRaw = amatr     {-# INLINE applyRaw #-} -infixl 1 #+infixr 1 # a # b = apply a b {-# INLINE (#) #-} +a #! b = a # b # id+{-# INLINE (#!) #-}+ --------------------------------------------------------------------------------  copy ord m = extractR ord m 0 (idxs[0,rows m-1]) 0 (idxs[0,cols m-1])@@ -427,7 +429,8 @@     let nr = if moder == 0 then fromIntegral $ vr@>1 - vr@>0 + 1 else dim vr         nc = if modec == 0 then fromIntegral $ vc@>1 - vc@>0 + 1 else dim vc     r <- createMatrix ord nr nc-    f moder modec # vr # vc # m # r  #|"extract"+    (vr # vc # m #! r) (f moder modec)  #|"extract"+     return r  type Extr x = CInt -> CInt -> CIdxs (CIdxs (OM x (OM x (IO CInt))))@@ -441,7 +444,7 @@  --------------------------------------------------------------- -setRectAux f i j m r = f (fi i) (fi j) # m # r #|"setRect"+setRectAux f i j m r = (m #! r) (f (fi i) (fi j)) #|"setRect"  type SetRect x = I -> I -> x ::> x::> Ok @@ -456,7 +459,7 @@  sortG f v = unsafePerformIO $ do     r <- createVector (dim v)-    f # v # r #|"sortG"+    (v #! r) f #|"sortG"     return r  sortIdxD = sortG c_sort_indexD@@ -483,7 +486,7 @@  compareG f u v = unsafePerformIO $ do     r <- createVector (dim v)-    f # u # v # r #|"compareG"+    (u # v #! r) f #|"compareG"     return r  compareD = compareG c_compareD@@ -500,7 +503,7 @@  selectG f c u v w = unsafePerformIO $ do     r <- createVector (dim v)-    f # c # u # v # w # r #|"selectG"+    (c # u # v # w #! r) f #|"selectG"     return r  selectD = selectG c_selectD@@ -523,7 +526,7 @@  remapG f i j m = unsafePerformIO $ do     r <- createMatrix RowMajor (rows i) (cols i)-    f # i # j # m # r #|"remapG"+    (i # j # m #! r) f #|"remapG"     return r  remapD = remapG c_remapD@@ -546,7 +549,7 @@  rowOpAux f c x i1 i2 j1 j2 m = do     px <- newArray [x]-    f (fi c) px (fi i1) (fi i2) (fi j1) (fi j2) # m #|"rowOp"+    (m # id) (f (fi c) px (fi i1) (fi i2) (fi j1) (fi j2)) #|"rowOp"     free px  type RowOp x = CInt -> Ptr x -> CInt -> CInt -> CInt -> CInt -> x ::> Ok@@ -562,7 +565,7 @@  -------------------------------------------------------------------------------- -gemmg f v m1 m2 m3 = f # v # m1 # m2 # m3 #|"gemmg"+gemmg f v m1 m2 m3 = (v # m1 # m2 #! m3) f #|"gemmg"  type Tgemm x = x :> x ::> x ::> x ::> Ok @@ -590,10 +593,9 @@ saveMatrix name format m = do     cname   <- newCString name     cformat <- newCString format-    c_saveMatrix cname cformat # m #|"saveMatrix"+    (m # id) (c_saveMatrix cname cformat) #|"saveMatrix"     free cname     free cformat     return ()  ---------------------------------------------------------------------------------

@@ -72,7 +72,7 @@   where     compare a b = compare (unMod a) (unMod b) -instance (Integral t, KnownNat m, Integral (Mod m t)) => Real (Mod m t)+instance (Integral t, Real t, KnownNat m, Integral (Mod m t)) => Real (Mod m t)   where     toRational x = toInteger x % 1 @@ -114,7 +114,7 @@   where     show = show . unMod -instance forall n t . (Integral t, KnownNat n) => Num (Mod n t)+instance (Integral t, KnownNat n) => Num (Mod n t)   where     (+) = l2 (\m a b -> (a + b) `mod` (fromIntegral m))     (*) = l2 (\m a b -> (a * b) `mod` (fromIntegral m))@@ -159,11 +159,11 @@         m' = fromIntegral . natVal $ (undefined :: Proxy m)  -instance forall m . KnownNat m => CTrans (Mod m I)-instance forall m . KnownNat m => CTrans (Mod m Z)+instance KnownNat m => CTrans (Mod m I)+instance KnownNat m => CTrans (Mod m Z)  -instance forall m . KnownNat m => Container Vector (Mod m I)+instance KnownNat m => Container Vector (Mod m I)   where     conj' = id     size' = dim@@ -203,7 +203,7 @@     fromZ'   = vmod . fromZ'     toZ'     = toZ' . f2i -instance forall m . KnownNat m => Container Vector (Mod m Z)+instance KnownNat m => Container Vector (Mod m Z)   where     conj' = id     size' = dim@@ -311,7 +311,7 @@  lift2m f a b = liftMatrix vmod (f (f2iM a) (f2iM b)) -instance forall m . KnownNat m => Num (Vector (Mod m I))+instance KnownNat m => Num (Vector (Mod m I))   where     (+) = lift2 (+)     (*) = lift2 (*)@@ -321,7 +321,7 @@     negate = lift1 negate     fromInteger x = fromInt (fromInteger x) -instance forall m . KnownNat m => Num (Vector (Mod m Z))+instance KnownNat m => Num (Vector (Mod m Z))   where     (+) = lift2 (+)     (*) = lift2 (*)

@@ -1,4 +1,3 @@-{-# LANGUAGE CPP #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-}

@@ -31,13 +31,13 @@ gaussianSample :: Seed                -> Int -- ^ number of rows                -> Vector Double -- ^ mean vector-               -> Matrix Double -- ^ covariance matrix+               -> Herm Double   -- ^ covariance matrix                -> Matrix Double -- ^ result gaussianSample seed n med cov = m where     c = dim med     meds = konst' 1 n `outer` med     rs = reshape c $ randomVector seed Gaussian (c * n)-    m = rs `mXm` cholSH cov `add` meds+    m = rs `mXm` chol cov `add` meds  -- | Obtains a matrix whose rows are pseudorandom samples from a multivariate -- uniform distribution.

@@ -144,13 +144,13 @@ gmXv SparseR { gmCSR = CSR{..}, .. } v = unsafePerformIO $ do     dim v /= nCols ~!~ printf "gmXv (CSR): incorrect sizes: (%d,%d) x %d" nRows nCols (dim v)     r <- createVector nRows-    c_smXv # csrVals # csrCols # csrRows # v # r #|"CSRXv"+    (csrVals # csrCols # csrRows # v #! r) c_smXv #|"CSRXv"     return r  gmXv SparseC { gmCSC = CSC{..}, .. } v = unsafePerformIO $ do     dim v /= nCols ~!~ printf "gmXv (CSC): incorrect sizes: (%d,%d) x %d" nRows nCols (dim v)     r <- createVector nRows-    c_smTXv # cscVals # cscRows # cscCols # v # r #|"CSCXv"+    (cscVals # cscRows # cscCols # v #! r) c_smTXv #|"CSCXv"     return r  gmXv Diag{..} v@@ -211,4 +211,3 @@     tr (Diag v n m) = Diag v m n     tr (Dense a n m) = Dense (tr a) m n     tr' = tr-

@@ -1,3 +1,4 @@+{-# LANGUAGE CPP #-} #if __GLASGOW_HASKELL__ >= 708  {-# LANGUAGE DataKinds #-}@@ -11,6 +12,8 @@ {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE TypeOperators #-} {-# LANGUAGE ViewPatterns #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DeriveGeneric #-}  {- | Module      :  Internal.Static@@ -28,18 +31,35 @@ import Numeric.LinearAlgebra hiding (konst,size,R,C) import Internal.Vector as D hiding (R,C) import Internal.ST+import Control.DeepSeq import Data.Proxy(Proxy) import Foreign.Storable(Storable) import Text.Printf +import Data.Binary+import GHC.Generics (Generic)+import Data.Proxy (Proxy(..))+ --------------------------------------------------------------------------------  type ℝ = Double type ℂ = Complex Double  newtype Dim (n :: Nat) t = Dim t-  deriving Show+  deriving (Show, Generic) +instance (KnownNat n, Binary a) => Binary (Dim n a) where+  get = do+    k <- get+    let n = natVal (Proxy :: Proxy n)+    if n == k+      then Dim <$> get+      else fail ("Expected dimension " ++ (show n) ++ ", but found dimension " ++ (show k))++  put (Dim x) = do+    put (natVal (Proxy :: Proxy n))+    put x+ lift1F   :: (c t -> c t)   -> Dim n (c t) -> Dim n (c t)@@ -50,18 +70,22 @@   -> Dim n (c t) -> Dim n (c t) -> Dim n (c t) lift2F f (Dim u) (Dim v) = Dim (f u v) +instance NFData t => NFData (Dim n t) where+    rnf (Dim (force -> !_)) = ()+ --------------------------------------------------------------------------------  newtype R n = R (Dim n (Vector ℝ))-  deriving (Num,Fractional,Floating)+  deriving (Num,Fractional,Floating,Generic,Binary)  newtype C n = C (Dim n (Vector ℂ))-  deriving (Num,Fractional,Floating)+  deriving (Num,Fractional,Floating,Generic)  newtype L m n = L (Dim m (Dim n (Matrix ℝ)))--newtype M m n = M (Dim m (Dim n (Matrix  ℂ)))+  deriving (Generic, Binary) +newtype M m n = M (Dim m (Dim n (Matrix ℂ)))+  deriving (Generic)  mkR :: Vector ℝ -> R n mkR = R . Dim@@ -75,6 +99,18 @@ mkM :: Matrix ℂ -> M m n mkM x = M (Dim (Dim x)) +instance NFData (R n) where+    rnf (R (force -> !_)) = ()++instance NFData (C n) where+    rnf (C (force -> !_)) = ()++instance NFData (L n m) where+    rnf (L (force -> !_)) = ()++instance NFData (M n m) where+    rnf (M (force -> !_)) = ()+ --------------------------------------------------------------------------------  type V n t = Dim n (Vector t)@@ -173,7 +209,7 @@ singleM m = rows m == 1 && cols m == 1  -instance forall n. KnownNat n => Sized ℂ (C n) Vector+instance KnownNat n => Sized ℂ (C n) Vector   where     size _ = fromIntegral . natVal $ (undefined :: Proxy n)     konst x = mkC (LA.scalar x)@@ -189,7 +225,7 @@         r = mkC v :: C n  -instance forall n. KnownNat n => Sized ℝ (R n) Vector+instance KnownNat n => Sized ℝ (R n) Vector   where     size _ = fromIntegral . natVal $ (undefined :: Proxy n)     konst x = mkR (LA.scalar x)@@ -206,7 +242,7 @@   -instance forall m n . (KnownNat m, KnownNat n) => Sized ℝ (L m n) Matrix+instance (KnownNat m, KnownNat n) => Sized ℝ (L m n) Matrix   where     size _ = ((fromIntegral . natVal) (undefined :: Proxy m)              ,(fromIntegral . natVal) (undefined :: Proxy n))@@ -224,7 +260,7 @@         r = mkL x :: L m n  -instance forall m n . (KnownNat m, KnownNat n) => Sized ℂ (M m n) Matrix+instance (KnownNat m, KnownNat n) => Sized ℂ (M m n) Matrix   where     size _ = ((fromIntegral . natVal) (undefined :: Proxy m)              ,(fromIntegral . natVal) (undefined :: Proxy n))@@ -282,7 +318,7 @@  -------------------------------------------------------------------------------- -instance forall n . KnownNat n => Show (R n)+instance KnownNat n => Show (R n)   where     show s@(R (Dim v))       | singleV v = "("++show (v!0)++" :: R "++show d++")"@@ -290,7 +326,7 @@       where         d = size s -instance forall n . KnownNat n => Show (C n)+instance KnownNat n => Show (C n)   where     show s@(C (Dim v))       | singleV v = "("++show (v!0)++" :: C "++show d++")"@@ -298,7 +334,7 @@       where         d = size s -instance forall m n . (KnownNat m, KnownNat n) => Show (L m n)+instance (KnownNat m, KnownNat n) => Show (L m n)   where     show (isDiag -> Just (z,y,(m',n'))) = printf "(diag %s %s :: L %d %d)" (show z) (drop 9 $ show y) m' n'     show s@(L (Dim (Dim x)))@@ -307,7 +343,7 @@       where         (m',n') = size s -instance forall m n . (KnownNat m, KnownNat n) => Show (M m n)+instance (KnownNat m, KnownNat n) => Show (M m n)   where     show (isDiagC -> Just (z,y,(m',n'))) = printf "(diag %s %s :: M %d %d)" (show z) (drop 9 $ show y) m' n'     show s@(M (Dim (Dim x)))@@ -318,7 +354,7 @@  -------------------------------------------------------------------------------- -instance forall n t . (Num (Vector t), Numeric t )=> Num (Dim n (Vector t))+instance (Num (Vector t), Numeric t )=> Num (Dim n (Vector t))   where     (+) = lift2F (+)     (*) = lift2F (*)@@ -482,6 +518,18 @@     (**)  = lift2MD (**)     pi    = M pi +instance Additive (R n) where+    add = (+)++instance Additive (C n) where+    add = (+)++instance (KnownNat m, KnownNat n) => Additive (L m n) where+    add = (+)++instance (KnownNat m, KnownNat n) => Additive (M m n) where+    add = (+)+ --------------------------------------------------------------------------------  @@ -518,6 +566,17 @@         putStr "C " >> putStr (tail . dropWhile (/='x') $ su)  --------------------------------------------------------------------------------++overMatL' :: (KnownNat m, KnownNat n)+          => (LA.Matrix ℝ -> LA.Matrix ℝ) -> L m n -> L m n+overMatL' f = mkL . f . unwrap+{-# INLINE overMatL' #-}++overMatM' :: (KnownNat m, KnownNat n)+          => (LA.Matrix ℂ -> LA.Matrix ℂ) -> M m n -> M m n+overMatM' f = mkM . f . unwrap+{-# INLINE overMatM' #-}+  #else 

@@ -41,7 +41,7 @@     ℕ,ℤ,ℝ,ℂ,iC,     Normed(..), norm_Frob, norm_nuclear,     magnit,-    unitary,+    normalize,     mt,     (~!~),     pairwiseD2,@@ -255,6 +255,7 @@ -- ^ 2-norm of real vector norm = pnorm PNorm2 +-- | p-norm for vectors, operator norm for matrices class Normed a   where     norm_0   :: a -> R@@ -319,10 +320,11 @@     norm_2 = norm_2 . double     norm_Inf = norm_Inf . double -+-- | Frobenius norm (Schatten p-norm with p=2) norm_Frob :: (Normed (Vector t), Element t) => Matrix t -> R norm_Frob = norm_2 . flatten +-- | Sum of singular values (Schatten p-norm with p=1) norm_nuclear :: Field t => Matrix t -> R norm_nuclear = sumElements . singularValues @@ -341,8 +343,8 @@   -- | Obtains a vector in the same direction with 2-norm=1-unitary :: Vector Double -> Vector Double-unitary v = v / scalar (norm v)+normalize :: (Normed (Vector t), Num (Vector t), Field t) => Vector t -> Vector t+normalize v = v / real (scalar (norm_2 v))   -- | trans . inv

@@ -1,4 +1,4 @@-{-# LANGUAGE MagicHash, CPP, UnboxedTuples, BangPatterns, FlexibleContexts #-}+{-# LANGUAGE MagicHash, UnboxedTuples, BangPatterns, FlexibleContexts #-} {-# LANGUAGE TypeSynonymInstances #-}  @@ -39,12 +39,10 @@ import qualified Data.Vector.Storable as Vector import Data.Vector.Storable(Vector, fromList, unsafeToForeignPtr, unsafeFromForeignPtr, unsafeWith) -#ifdef BINARY import Data.Binary import Control.Monad(replicateM) import qualified Data.ByteString.Internal as BS import Data.Vector.Storable.Internal(updPtr)-#endif  type I = CInt type Z = Int64@@ -64,13 +62,13 @@ -- | Number of elements dim :: (Storable t) => Vector t -> Int dim = Vector.length+{-# INLINE dim #-}   -- C-Haskell vector adapter {-# INLINE avec #-}-avec :: Storable a => (CInt -> Ptr a -> b) -> Vector a -> b-avec f v = inlinePerformIO (unsafeWith v (return . f (fromIntegral (Vector.length v))))-infixl 1 `avec`+avec :: Storable a => Vector a -> (f -> IO r) -> ((CInt -> Ptr a -> f) -> IO r)+avec v f g = unsafeWith v $ \ptr -> f (g (fromIntegral (Vector.length v)) ptr)  -- allocates memory for a new vector createVector :: Storable a => Int -> IO (Vector a)@@ -142,6 +140,7 @@                         -> Vector t  -- ^ source                         -> Vector t  -- ^ result subVector = Vector.slice+{-# INLINE subVector #-}   @@ -201,7 +200,7 @@  --------------------------------------------------------------- --- | transforms a complex vector into a real vector with alternating real and imaginary parts +-- | transforms a complex vector into a real vector with alternating real and imaginary parts asReal :: (RealFloat a, Storable a) => Vector (Complex a) -> Vector a asReal v = unsafeFromForeignPtr (castForeignPtr fp) (2*i) (2*n)     where (fp,i,n) = unsafeToForeignPtr v@@ -246,7 +245,7 @@ {-# INLINE zipVectorWith #-}  -- | unzipWith for Vectors-unzipVectorWith :: (Storable (a,b), Storable c, Storable d) +unzipVectorWith :: (Storable (a,b), Storable c, Storable d)                    => ((a,b) -> (c,d)) -> Vector (a,b) -> (Vector c,Vector d) unzipVectorWith f u = unsafePerformIO $ do       let n = dim u@@ -257,7 +256,7 @@               unsafeWith w $ \pw -> do                   let go (-1) = return ()                       go !k   = do z <- peekElemOff pu k-                                   let (x,y) = f z +                                   let (x,y) = f z                                    pokeElemOff      pv k x                                    pokeElemOff      pw k y                                    go (k-1)@@ -305,11 +304,11 @@     return w     where mapVectorM' w' !k !t               | k == t               = do-                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k +                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k                                        y <- f x                                        return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y               | otherwise            = do-                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k +                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k                                        y <- f x                                        _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y                                        mapVectorM' w' (k+1) t@@ -324,7 +323,7 @@                                     x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k                                     f x               | otherwise         = do-                                    x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k +                                    x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k                                     _ <- f x                                     mapVectorM' (k+1) t {-# INLINE mapVectorM_ #-}@@ -338,11 +337,11 @@     return w     where mapVectorM' w' !k !t               | k == t               = do-                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k +                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k                                        y <- f k x                                        return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y               | otherwise            = do-                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k +                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k                                        y <- f k x                                        _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y                                        mapVectorM' w' (k+1) t@@ -357,7 +356,7 @@                                     x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k                                     f k x               | otherwise         = do-                                    x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k +                                    x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k                                     _ <- f k x                                     mapVectorM' (k+1) t {-# INLINE mapVectorWithIndexM_ #-}@@ -380,7 +379,6 @@ --------------------------------------------------------------------------------  -#ifdef BINARY  -- a 64K cache, with a Double taking 13 bytes in Bytestring, -- implies a chunk size of 5041@@ -432,7 +430,6 @@      -- get = fmap bs2v get -#endif   -------------------------------------------------------------------@@ -458,4 +455,3 @@ unzipVector = unzipVectorWith id  --------------------------------------------------------------------

@@ -27,10 +27,13 @@ import System.IO.Unsafe(unsafePerformIO) import Control.Monad(when) -infixl 1 #+infixr 1 # a # b = applyRaw a b {-# INLINE (#) #-} +a #! b = a # b # id+{-# INLINE (#!) #-}+ fromei x = fromIntegral (fromEnum x) :: CInt  data FunCodeV = Sin@@ -103,7 +106,7 @@  sumg f x = unsafePerformIO $ do     r <- createVector 1-    f # x # r #| "sum"+    (x #! r) f #| "sum"     return $ r @> 0  type TVV t = t :> t :> Ok@@ -139,7 +142,7 @@  prodg f x = unsafePerformIO $ do     r <- createVector 1-    f # x # r #| "prod"+    (x #! r) f #| "prod"     return $ r @> 0  @@ -154,24 +157,24 @@  toScalarAux fun code v = unsafePerformIO $ do     r <- createVector 1-    fun (fromei code) # v # r #|"toScalarAux"+    (v #! r) (fun (fromei code)) #|"toScalarAux"     return (r @> 0)  vectorMapAux fun code v = unsafePerformIO $ do     r <- createVector (dim v)-    fun (fromei code) # v # r #|"vectorMapAux"+    (v #! r) (fun (fromei code)) #|"vectorMapAux"     return r  vectorMapValAux fun code val v = unsafePerformIO $ do     r <- createVector (dim v)     pval <- newArray [val]-    fun (fromei code) pval # v # r #|"vectorMapValAux"+    (v #! r) (fun (fromei code) pval) #|"vectorMapValAux"     free pval     return r  vectorZipAux fun code u v = unsafePerformIO $ do     r <- createVector (dim u)-    fun (fromei code) # u # v # r #|"vectorZipAux"+    (u # v #! r) (fun (fromei code)) #|"vectorZipAux"     return r  ---------------------------------------------------------------------@@ -368,7 +371,7 @@              -> Vector Double randomVector seed dist n = unsafePerformIO $ do     r <- createVector n-    c_random_vector (fi seed) ((fi.fromEnum) dist) # r #|"randomVector"+    (r # id) (c_random_vector (fi seed) ((fi.fromEnum) dist)) #|"randomVector"     return r  foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> Double :> Ok@@ -377,7 +380,7 @@  roundVector v = unsafePerformIO $ do     r <- createVector (dim v)-    c_round_vector # v # r #|"roundVector"+    (v #! r) c_round_vector #|"roundVector"     return r  foreign import ccall unsafe "round_vector" c_round_vector :: TVV Double@@ -391,7 +394,7 @@ range :: Int -> Vector I range n = unsafePerformIO $ do     r <- createVector n-    c_range_vector # r #|"range"+    (r # id) c_range_vector #|"range"     return r  foreign import ccall unsafe "range_vector" c_range_vector :: CInt :> Ok@@ -431,7 +434,7 @@  tog f v = unsafePerformIO $ do     r <- createVector (dim v)-    f # v # r #|"tog"+    (v #! r) f #|"tog"     return r  foreign import ccall unsafe "float2double" c_float2double :: Float :> Double :> Ok@@ -450,7 +453,7 @@  stepg f v = unsafePerformIO $ do     r <- createVector (dim v)-    f # v # r #|"step"+    (v #! r) f #|"step"     return r  stepD :: Vector Double -> Vector Double@@ -475,7 +478,7 @@  conjugateAux fun x = unsafePerformIO $ do     v <- createVector (dim x)-    fun # x # v #|"conjugateAux"+    (x #! v) fun #|"conjugateAux"     return v  conjugateQ :: Vector (Complex Float) -> Vector (Complex Float)@@ -493,7 +496,7 @@         let n = dim v         r <- createVector n         let f _ s _ d =  copyArray d s n >> return 0-        f # v # r #|"cloneVector"+        (v #! r) f #|"cloneVector"         return r  --------------------------------------------------------------------------------@@ -501,7 +504,7 @@ constantAux fun x n = unsafePerformIO $ do     v <- createVector n     px <- newArray [x]-    fun px # v #|"constantAux"+    (v # id) (fun px) #|"constantAux"     free px     return v @@ -515,4 +518,3 @@ foreign import ccall unsafe "constantL" cconstantL :: TConst Z  -----------------------------------------------------------------------

@@ -22,7 +22,8 @@      -- * Numeric classes     -- |-    -- The standard numeric classes are defined elementwise:+    -- The standard numeric classes are defined elementwise (commonly referred to+    -- as the Hadamard product or the Schur product):     --     -- >>>  vector [1,2,3] * vector [3,0,-2]     -- fromList [3.0,0.0,-6.0]@@ -116,6 +117,7 @@     svd,     thinSVD,     compactSVD,+    compactSVDTol,     singularValues,     leftSV, rightSV, @@ -125,7 +127,7 @@     geigSH,      -- * QR-    qr, rq, qrRaw, qrgr,+    qr, thinQR, rq, thinRQ, qrRaw, qrgr,      -- * Cholesky     chol, mbChol,@@ -152,7 +154,7 @@     Seed, RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample,      -- * Misc-    meanCov, rowOuters, pairwiseD2, unitary, peps, relativeError, magnit,+    meanCov, rowOuters, pairwiseD2, normalize, peps, relativeError, magnit,     haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv,     iC, sym, mTm, trustSym, unSym,     -- * Auxiliary classes

@@ -12,13 +12,6 @@ --------------------------------------------------------------------------------  module Numeric.LinearAlgebra.Devel(-    -- * FFI helpers-    -- | Sample usage, to upload a perspective matrix to a shader.-    ---    -- @ glUniformMatrix4fv 0 1 (fromIntegral gl_TRUE) \`appMatrix\` perspective 0.01 100 (pi\/2) (4\/3)-    -- @-    module Internal.Foreign,-     -- * FFI tools     -- | See @examples/devel@ in the repository.     @@ -66,7 +59,6 @@  ) where -import Internal.Foreign import Internal.Devel import Internal.ST import Internal.Vector

@@ -42,19 +42,24 @@     blockAt,     matrix,     -- * Complex-    C, M, Her, her, 𝑖,+    ℂ, C, M, Her, her, 𝑖,     -- * Products     (<>),(#>),(<.>),     -- * Linear Systems     linSolve, (<\>),     -- * Factorizations     svd, withCompactSVD, svdTall, svdFlat, Eigen(..),-    withNullspace, qr, chol,+    withNullspace, withOrth, qr, chol,+    -- * Norms+    Normed(..),+    -- * Random arrays+    Seed, RandDist(..),+    randomVector, rand, randn, gaussianSample, uniformSample,     -- * Misc-    mean,+    mean, meanCov,     Disp(..), Domain(..),-    withVector, withMatrix,-    toRows, toColumns,+    withVector, withMatrix, exactLength, exactDims,+    toRows, toColumns, withRows, withColumns,     Sized(..), Diag(..), Sym, sym, mTm, unSym, (<·>) ) where @@ -65,11 +70,16 @@     row,col,vector,matrix,linspace,toRows,toColumns,     (<\>),fromList,takeDiag,svd,eig,eigSH,     eigenvalues,eigenvaluesSH,build,-    qr,size,dot,chol,range,R,C,sym,mTm,unSym)+    qr,size,dot,chol,range,R,C,sym,mTm,unSym,+    randomVector,rand,randn,gaussianSample,uniformSample,meanCov) import qualified Numeric.LinearAlgebra as LA-import Data.Proxy(Proxy)+import qualified Numeric.LinearAlgebra.Devel as LA+import Data.Proxy(Proxy(..)) import Internal.Static import Control.Arrow((***))+import Text.Printf+import Data.Type.Equality ((:~:)(Refl))+import qualified Data.Bifunctor as BF (first)  ud1 :: R n -> Vector ℝ ud1 (R (Dim v)) = v@@ -218,19 +228,14 @@     takeDiag :: m -> d  -instance forall n . (KnownNat n) => Diag (L n n) (R n)-  where-    takeDiag m = mkR (LA.takeDiag (extract m))---instance forall m n . (KnownNat m, KnownNat n, m <= n+1) => Diag (L m n) (R m)+instance KnownNat n => Diag (L n n) (R n)   where-    takeDiag m = mkR (LA.takeDiag (extract m))+    takeDiag x = mkR (LA.takeDiag (extract x))  -instance forall m n . (KnownNat m, KnownNat n, n <= m+1) => Diag (L m n) (R n)+instance KnownNat n => Diag (M n n) (C n)   where-    takeDiag m = mkR (LA.takeDiag (extract m))+    takeDiag x = mkC (LA.takeDiag (extract x))   --------------------------------------------------------------------------------@@ -289,7 +294,21 @@ her m = Her $ (m + LA.tr m)/2  +instance (KnownNat n) => Disp (Sym n)+  where+    disp n (Sym x) = do+        let a = extract x+        let su = LA.dispf n a+        printf "Sym %d" (cols a) >> putStr (dropWhile (/='\n') $ su) +instance (KnownNat n) => Disp (Her n)+  where+    disp n (Her x) = do+        let a = extract x+        let su = LA.dispcf n a+        printf "Her %d" (cols a) >> putStr (dropWhile (/='\n') $ su)++ instance KnownNat n => Eigen (Sym n) (R n) (L n n)   where     eigenvalues (Sym (extract -> m)) =  mkR . LA.eigenvaluesSH . LA.trustSym $ m@@ -319,6 +338,15 @@        Nothing -> error "static/dynamic mismatch"        Just (SomeNat (_ :: Proxy k)) -> f (mkL a :: L n k) +withOrth+    :: forall m n z . (KnownNat m, KnownNat n)+    => L m n+    -> (forall k. (KnownNat k) => L n k -> z)+    -> z+withOrth (LA.orth . extract -> a) f =+    case someNatVal $ fromIntegral $ cols a of+       Nothing -> error "static/dynamic mismatch"+       Just (SomeNat (_ :: Proxy k)) -> f (mkL a :: L n k)  withCompactSVD     :: forall m n z . (KnownNat m, KnownNat n)@@ -365,11 +393,31 @@ toRows :: forall m n . (KnownNat m, KnownNat n) => L m n -> [R n] toRows (LA.toRows . extract -> vs) = map mkR vs +withRows+    :: forall n z . KnownNat n+    => [R n]+    -> (forall m . KnownNat m => L m n -> z)+    -> z+withRows (LA.fromRows . map extract -> m) f =+    case someNatVal $ fromIntegral $ LA.rows m of+       Nothing -> error "static/dynamic mismatch"+       Just (SomeNat (_ :: Proxy m)) -> f (mkL m :: L m n)  toColumns :: forall m n . (KnownNat m, KnownNat n) => L m n -> [R m] toColumns (LA.toColumns . extract -> vs) = map mkR vs +withColumns+    :: forall m z . KnownNat m+    => [R m]+    -> (forall n . KnownNat n => L m n -> z)+    -> z+withColumns (LA.fromColumns . map extract -> m) f =+    case someNatVal $ fromIntegral $ LA.cols m of+       Nothing -> error "static/dynamic mismatch"+       Just (SomeNat (_ :: Proxy n)) -> f (mkL m :: L m n) ++ --------------------------------------------------------------------------------  build@@ -392,6 +440,15 @@        Nothing -> error "static/dynamic mismatch"        Just (SomeNat (_ :: Proxy m)) -> f (mkR v :: R m) +-- | Useful for constraining two dependently typed vectors to match each+-- other in length when they are unknown at compile-time.+exactLength+    :: forall n m . (KnownNat n, KnownNat m)+    => R m+    -> Maybe (R n)+exactLength v = do+    Refl <- sameNat (Proxy :: Proxy n) (Proxy :: Proxy m)+    return $ mkR (unwrap v)  withMatrix     :: forall z@@ -407,6 +464,64 @@                Just (SomeNat (_ :: Proxy n)) ->                   f (mkL a :: L m n) +-- | Useful for constraining two dependently typed matrices to match each+-- other in dimensions when they are unknown at compile-time.+exactDims+    :: forall n m j k . (KnownNat n, KnownNat m, KnownNat j, KnownNat k)+    => L m n+    -> Maybe (L j k)+exactDims m = do+    Refl <- sameNat (Proxy :: Proxy m) (Proxy :: Proxy j)+    Refl <- sameNat (Proxy :: Proxy n) (Proxy :: Proxy k)+    return $ mkL (unwrap m)++randomVector+    :: forall n . KnownNat n+    => Seed+    -> RandDist+    -> R n+randomVector s d = mkR (LA.randomVector s d+                          (fromInteger (natVal (Proxy :: Proxy n)))+                       )++rand+    :: forall m n . (KnownNat m, KnownNat n)+    => IO (L m n)+rand = mkL <$> LA.rand (fromInteger (natVal (Proxy :: Proxy m)))+                       (fromInteger (natVal (Proxy :: Proxy n)))++randn+    :: forall m n . (KnownNat m, KnownNat n)+    => IO (L m n)+randn = mkL <$> LA.randn (fromInteger (natVal (Proxy :: Proxy m)))+                         (fromInteger (natVal (Proxy :: Proxy n)))++gaussianSample+    :: forall m n . (KnownNat m, KnownNat n)+    => Seed+    -> R n+    -> Sym n+    -> L m n+gaussianSample s (extract -> mu) (Sym (extract -> sigma)) =+    mkL $ LA.gaussianSample s (fromInteger (natVal (Proxy :: Proxy m)))+                            mu (LA.trustSym sigma)++uniformSample+    :: forall m n . (KnownNat m, KnownNat n)+    => Seed+    -> R n    -- ^ minimums of each row+    -> R n    -- ^ maximums of each row+    -> L m n+uniformSample s (extract -> mins) (extract -> maxs) =+    mkL $ LA.uniformSample s (fromInteger (natVal (Proxy :: Proxy m)))+                           (zip (LA.toList mins) (LA.toList maxs))++meanCov+    :: forall m n . (KnownNat m, KnownNat n, 1 <= m)+    => L m n+    -> (R n, Sym n)+meanCov (extract -> vs) = mkR *** (Sym . mkL . LA.unSym) $ LA.meanCov vs+ --------------------------------------------------------------------------------  class Domain field vec mat | mat -> vec field, vec -> mat field, field -> mat vec@@ -416,6 +531,15 @@     dot :: forall n . (KnownNat n) => vec n -> vec n -> field     cross :: vec 3 -> vec 3 -> vec 3     diagR ::  forall m n k . (KnownNat m, KnownNat n, KnownNat k) => field -> vec k -> mat m n+    dvmap :: forall n. KnownNat n => (field -> field) -> vec n -> vec n+    dmmap :: forall n m. (KnownNat m, KnownNat n) => (field -> field) -> mat n m -> mat n m+    outer :: forall n m. (KnownNat m, KnownNat n) => vec n -> vec m -> mat n m+    zipWithVector :: forall n. KnownNat n => (field -> field -> field) -> vec n -> vec n -> vec n+    det :: forall n. KnownNat n => mat n n -> field+    invlndet :: forall n. KnownNat n => mat n n -> (mat n n, (field, field))+    expm :: forall n. KnownNat n => mat n n -> mat n n+    sqrtm :: forall n. KnownNat n => mat n n -> mat n n+    inv :: forall n. KnownNat n => mat n n -> mat n n   instance Domain ℝ R L@@ -425,6 +549,15 @@     dot = dotR     cross = crossR     diagR = diagRectR+    dvmap = mapR+    dmmap = mapL+    outer = outerR+    zipWithVector = zipWithR+    det = detL+    invlndet = invlndetL+    expm = expmL+    sqrtm = sqrtmL+    inv = invL  instance Domain ℂ C M   where@@ -433,6 +566,15 @@     dot = dotC     cross = crossC     diagR = diagRectC+    dvmap = mapC+    dmmap = mapM'+    outer = outerC+    zipWithVector = zipWithC+    det = detM+    invlndet = invlndetM+    expm = expmM+    sqrtm = sqrtmM+    inv = invM  -------------------------------------------------------------------------------- @@ -472,6 +614,33 @@     z2 = x!2*y!0-x!0*y!2     z3 = x!0*y!1-x!1*y!0 +outerR :: (KnownNat m, KnownNat n) => R n -> R m -> L n m+outerR (extract -> x) (extract -> y) = mkL (LA.outer x y)++mapR :: KnownNat n => (ℝ -> ℝ) -> R n -> R n+mapR f (unwrap -> v) = mkR (LA.cmap f v)++zipWithR :: KnownNat n => (ℝ -> ℝ -> ℝ) -> R n -> R n -> R n+zipWithR f (extract -> x) (extract -> y) = mkR (LA.zipVectorWith f x y)++mapL :: (KnownNat n, KnownNat m) => (ℝ -> ℝ) -> L n m -> L n m+mapL f = overMatL' (LA.cmap f)++detL :: KnownNat n => Sq n -> ℝ+detL = LA.det . unwrap++invlndetL :: KnownNat n => Sq n -> (L n n, (ℝ, ℝ))+invlndetL = BF.first mkL . LA.invlndet . unwrap++expmL :: KnownNat n => Sq n -> Sq n+expmL = overMatL' LA.expm++sqrtmL :: KnownNat n => Sq n -> Sq n+sqrtmL = overMatL' LA.sqrtm++invL :: KnownNat n => Sq n -> Sq n+invL = overMatL' LA.inv+ --------------------------------------------------------------------------------  mulC :: forall m k n. (KnownNat m, KnownNat k, KnownNat n) => M m k -> M k n -> M m n@@ -510,6 +679,33 @@     z2 = x!2*y!0-x!0*y!2     z3 = x!0*y!1-x!1*y!0 +outerC :: (KnownNat m, KnownNat n) => C n -> C m -> M n m+outerC (extract -> x) (extract -> y) = mkM (LA.outer x y)++mapC :: KnownNat n => (ℂ -> ℂ) -> C n -> C n+mapC f (unwrap -> v) = mkC (LA.cmap f v)++zipWithC :: KnownNat n => (ℂ -> ℂ -> ℂ) -> C n -> C n -> C n+zipWithC f (extract -> x) (extract -> y) = mkC (LA.zipVectorWith f x y)++mapM' :: (KnownNat n, KnownNat m) => (ℂ -> ℂ) -> M n m -> M n m+mapM' f = overMatM' (LA.cmap f)++detM :: KnownNat n => M n n -> ℂ+detM = LA.det . unwrap++invlndetM :: KnownNat n => M n n -> (M n n, (ℂ, ℂ))+invlndetM = BF.first mkM . LA.invlndet . unwrap++expmM :: KnownNat n => M n n -> M n n+expmM = overMatM' LA.expm++sqrtmM :: KnownNat n => M n n -> M n n+sqrtmM = overMatM' LA.sqrtm++invM :: KnownNat n => M n n -> M n n+invM = overMatM' LA.inv+ --------------------------------------------------------------------------------  diagRectR :: forall m n k . (KnownNat m, KnownNat n, KnownNat k) => ℝ -> R k -> L m n@@ -616,3 +812,67 @@   where     checkT _ = test +--------------------------------------------------------------------------------++instance KnownNat n => Normed (R n)+  where+    norm_0 v = norm_0 (extract v)+    norm_1 v = norm_1 (extract v)+    norm_2 v = norm_2 (extract v)+    norm_Inf v = norm_Inf (extract v)++instance (KnownNat m, KnownNat n) => Normed (L m n)+  where+    norm_0 m = norm_0 (extract m)+    norm_1 m = norm_1 (extract m)+    norm_2 m = norm_2 (extract m)+    norm_Inf m = norm_Inf (extract m)++mkSym f = Sym . f . unSym+mkSym2 f x y = Sym (f (unSym x) (unSym y))++instance KnownNat n =>  Num (Sym n)+  where+    (+) = mkSym2 (+)+    (*) = mkSym2 (*)+    (-) = mkSym2 (-)+    abs = mkSym abs+    signum = mkSym signum+    negate = mkSym negate+    fromInteger = Sym . fromInteger++instance KnownNat n => Fractional (Sym n)+  where+    fromRational = Sym . fromRational+    (/) = mkSym2 (/)++instance KnownNat n => Floating (Sym n)+  where+    sin   = mkSym sin+    cos   = mkSym cos+    tan   = mkSym tan+    asin  = mkSym asin+    acos  = mkSym acos+    atan  = mkSym atan+    sinh  = mkSym sinh+    cosh  = mkSym cosh+    tanh  = mkSym tanh+    asinh = mkSym asinh+    acosh = mkSym acosh+    atanh = mkSym atanh+    exp   = mkSym exp+    log   = mkSym log+    sqrt  = mkSym sqrt+    (**)  = mkSym2 (**)+    pi    = Sym pi++instance KnownNat n => Additive (Sym n) where+    add = (+)++instance KnownNat n => Transposable (Sym n) (Sym n) where+    tr  = id+    tr' = id++instance KnownNat n => Transposable (Her n) (Her n) where+    tr          = id+    tr' (Her m) = Her (tr' m)