fix static blocks, GMatrix

6 files changed, 55 insertions, 43 deletions

diff --git a/packages/base/src/Numeric/HMatrix.hs b/packages/base/src/Numeric/HMatrix.hs
index 1c70ef6..7f27fd4 100644
--- a/packages/base/src/Numeric/HMatrix.hs
+++ b/packages/base/src/Numeric/HMatrix.hs
@@ -41,7 +41,7 @@ module Numeric.HMatrix (
     -- ** dot
     (<·>),
     -- ** matrix-vector
-     (#>),(!#>),
+    (#>), (!#>),
     -- ** matrix-matrix
      (<>),
     -- | The matrix x matrix product is also implemented in the "Data.Monoid" instance, where
@@ -135,7 +135,7 @@ module Numeric.HMatrix (
     -- * Misc
     meanCov, peps, relativeError, haussholder, optimiseMult, udot,
     -- * Auxiliary classes
-    Element, Container, Product, Contraction(..), Numeric, LSDiv,
+    Element, Container, Product, Numeric, LSDiv,
     Complexable, RealElement,
     RealOf, ComplexOf, SingleOf, DoubleOf,
     IndexOf,
diff --git a/packages/base/src/Numeric/LinearAlgebra/Data.hs b/packages/base/src/Numeric/LinearAlgebra/Data.hs
index 33a2c9a..20f3b81 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Data.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Data.hs
@@ -68,7 +68,7 @@ module Numeric.LinearAlgebra.Data(
 
     module Data.Complex,
 
-    Vector, Matrix, GMatrix, CSR(..), mkCSR
+    Vector, Matrix, GMatrix, nRows, nCols
 
 ) where
 
diff --git a/packages/base/src/Numeric/LinearAlgebra/Devel.hs b/packages/base/src/Numeric/LinearAlgebra/Devel.hs
index ca9e53a..fce8b71 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Devel.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Devel.hs
@@ -49,10 +49,15 @@ module Numeric.LinearAlgebra.Devel(
     mapMatrixWithIndex, mapMatrixWithIndexM, mapMatrixWithIndexM_,
     liftMatrix, liftMatrix2, liftMatrix2Auto,
 
+    -- * Misc
+    CSR(..), fromCSR, mkCSR,
+    GMatrix(..)
+
 ) where
 
 import Data.Packed.Foreign
 import Data.Packed.Development
 import Data.Packed.ST
 import Data.Packed
+import Numeric.Sparse
 
diff --git a/packages/base/src/Numeric/LinearAlgebra/Real.hs b/packages/base/src/Numeric/LinearAlgebra/Real.hs
index 1e8b544..5634031 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Real.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Real.hs
@@ -29,26 +29,26 @@ module Numeric.LinearAlgebra.Real(
     vec2, vec3, vec4, 𝕧, (&),
     -- * Matrix
     L, Sq,
-    𝕞,
-    (#),(¦),(——),
-     Konst(..),
-     eye,
-     diagR, diag,
-     blockAt,
+    row, col, (¦),(——),
+    Konst(..),
+    eye,
+    diagR, diag,
+    blockAt,
     -- * Products
-     (<>),(#>),(<·>),
-     -- * Pretty printing
-     Disp(..),
-     -- * Misc
-     Dim, unDim,
+    (<>),(#>),(<·>),
+    -- * Pretty printing
+    Disp(..),
+    -- * Misc
+    Dim, unDim,
     module Numeric.HMatrix
 ) where
 
 
 import GHC.TypeLits
-import Numeric.HMatrix hiding ((<>),(#>),(<·>),Konst(..),diag, disp,(¦),(——))
+import Numeric.HMatrix hiding ((<>),(#>),(<·>),Konst(..),diag, disp,(¦),(——),row,col)
 import qualified Numeric.HMatrix as LA
 import Data.Packed.ST
+import Data.Proxy(Proxy)
 
 newtype Dim (n :: Nat) t = Dim t
   deriving Show
@@ -56,7 +56,7 @@ newtype Dim (n :: Nat) t = Dim t
 unDim :: Dim n t -> t
 unDim (Dim x) = x
 
-data Proxy :: Nat -> *
+-- data Proxy :: Nat -> *
 
 
 lift1F
@@ -223,7 +223,7 @@ instance Disp (R n)
             else putStr "Dim " >> putStr (tail . dropWhile (/='x') $ su)
 
 --------------------------------------------------------------------------------
-
+{-
 infixl 3 #
 (#) :: L r c -> R c -> L (r+1) c
 Dim (Dim m) # Dim v = Dim (Dim (m LA.—— asRow v))
@@ -233,14 +233,31 @@ Dim (Dim m) # Dim v = Dim (Dim (m LA.—— asRow v))
 𝕞  = Dim (Dim (LA.konst 0 (0,d)))
   where
     d = fromIntegral . natVal $ (undefined :: Proxy n)
+-}
+
+row :: R n -> L 1 n
+row (Dim v) = Dim (Dim (asRow v))
+
+col :: R n -> L n 1
+col = tr . row
 
 infixl 3 ¦
-(¦) :: L r c1 -> L r c2 -> L r (c1+c2)
-Dim (Dim a) ¦ Dim (Dim b) = Dim (Dim (a LA.¦ b))
+(¦) :: (KnownNat r, KnownNat c1, KnownNat c2) => L r c1 -> L r c2 -> L r (c1+c2)
+a ¦ b = rjoin (expk a) (expk b)
+  where
+    Dim (Dim a') `rjoin` Dim (Dim b') = Dim (Dim (a' LA.¦ b'))
 
 infixl 2 ——
-(——) :: L r1 c -> L r2 c -> L (r1+r2) c
-Dim (Dim a) —— Dim (Dim b) = Dim (Dim (a LA.—— b))
+(——) :: (KnownNat r1, KnownNat r2, KnownNat c) => L r1 c -> L r2 c -> L (r1+r2) c
+a —— b = cjoin (expk a) (expk b)
+  where
+    Dim (Dim a') `cjoin` Dim (Dim b') = Dim (Dim (a' LA.—— b'))
+
+expk :: (KnownNat n, KnownNat m) => L m n -> L m n
+expk x | singleton x = konst (d2 x `atIndex` (0,0))
+       | otherwise = x
+  where
+    singleton (d2 -> m) = rows m == 1 && cols m == 1
 
 
 {-
@@ -338,10 +355,4 @@ instance (KnownNat n', KnownNat m') => Testable (L n' m')
   where
     checkT _ = test
 
-{-
-do (snd test)
-fst test
--}
-
-
 
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs b/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs
index f821b57..b82c74f 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs
@@ -41,22 +41,21 @@ cg sym at a (CGState p r r2 x _) = CGState p' r' r'2 x' rdx
     ap1 = a p
     ap  | sym       = ap1
         | otherwise = at ap1
-    pap | sym       = p ◇ ap1
+    pap | sym       = p <·> ap1
         | otherwise = norm2 ap1 ** 2
     alpha = r2 / pap
     dx = scale alpha p
     x' = x + dx
     r' = r - scale alpha ap
-    r'2 = r' ◇ r'
+    r'2 = r' <·> r'
     beta = r'2 / r2
     p' = r' + scale beta p
 
     rdx = norm2 dx / max 1 (norm2 x)
 
 conjugrad
-  :: (Transposable m mt, Contraction m V V, Contraction mt V V)
-  => Bool -> m -> V -> V -> R -> R -> [CGState]
-conjugrad sym a b = solveG (tr a ◇) (a ◇) (cg sym) b
+  :: Bool -> GMatrix -> V -> V -> R -> R -> [CGState]
+conjugrad sym a b = solveG (tr a !#>) (a !#>) (cg sym) b
 
 solveG
     :: (V -> V) -> (V -> V)
@@ -72,9 +71,9 @@ solveG mat ma meth rawb x0' ϵb ϵx
     b = mat rawb
     x0  = if x0' == 0 then konst 0 (dim b) else x0'
     r0  = b - a x0
-    r20 = r0 ◇ r0
+    r20 = r0 <·> r0
     p0  = r0
-    nb2 = b ◇ b
+    nb2 = b <·> b
     ok CGState {..}
         =  cgr2 <nb2*ϵb**2
         || cgdx < ϵx
@@ -115,7 +114,7 @@ instance Testable GMatrix
         sma = convo2 20 3
         x1 = vect [1..20]
         x2 = vect [1..40]
-        sm = (mkSparse . mkCSR) sma
+        sm = mkSparse sma
         dm = toDense sma
 
         s1 = sm !#> x1
diff --git a/packages/base/src/Numeric/Sparse.hs b/packages/base/src/Numeric/Sparse.hs
index 3c19c93..1b8a7b3 100644
--- a/packages/base/src/Numeric/Sparse.hs
+++ b/packages/base/src/Numeric/Sparse.hs
@@ -3,7 +3,7 @@
 {-# LANGUAGE FlexibleInstances #-}
 
 module Numeric.Sparse(
-    GMatrix, CSR(..), mkCSR,
+    GMatrix(..), CSR(..), mkCSR, fromCSR,
     mkSparse, mkDiagR, mkDense,
     AssocMatrix,
     toDense,
@@ -95,9 +95,11 @@ mkDense m = Dense{..}
     nRows = rows m
     nCols = cols m
 
+mkSparse :: AssocMatrix -> GMatrix
+mkSparse = fromCSR . mkCSR
 
-mkSparse :: CSR -> GMatrix
-mkSparse csr = SparseR {..}
+fromCSR :: CSR -> GMatrix
+fromCSR csr = SparseR {..}
   where
     gmCSR @ CSR {..} = csr
     nRows = csrNRows
@@ -149,11 +151,6 @@ infixr 8 !#>
 (!#>) :: GMatrix -> Vector Double -> Vector Double
 (!#>) = gmXv
 
-
-instance Contraction GMatrix (Vector Double) (Vector Double)
-  where
-    contraction = gmXv
-
 --------------------------------------------------------------------------------
 
 foreign import ccall unsafe "smXv"