static dimensions, cont.

2 files changed, 75 insertions, 81 deletions

diff --git a/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs b/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs
index 5e2ea84..50372f1 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs
@@ -3,11 +3,14 @@
 
 module Numeric.LinearAlgebra.Util.CG(
     cgSolve, cgSolve',
-    CGMat, CGState(..), R, V
+    CGState(..), R, V
 ) where
 
 import Data.Packed.Numeric
+import Numeric.Sparse
 import Numeric.Vector()
+import Numeric.LinearAlgebra.Algorithms(linearSolveLS, relativeError, NormType(..))
+import Control.Arrow((***))
 
 {-
 import Util.Misc(debug, debugMat)
@@ -51,7 +54,7 @@ cg sym at a (CGState p r r2 x _) = CGState p' r' r'2 x' rdx
     rdx = norm2 dx / max 1 (norm2 x)
 
 conjugrad
-  :: (Transposable m, Contraction m V V)
+  :: (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
 
@@ -82,27 +85,88 @@ takeUntil q xs = a++ take 1 b
   where
     (a,b) = break q xs
 
-class (Transposable m, Contraction m V V) => CGMat m
-
 cgSolve
-  :: CGMat m
-  => Bool     -- ^ is symmetric
-  -> m        -- ^ coefficient matrix
+  :: Bool          -- ^ is symmetric
+  -> GMatrix       -- ^ coefficient matrix
   -> Vector Double -- ^ right-hand side
-  -> Vector Double        -- ^ solution
+  -> Vector Double -- ^ solution
 cgSolve sym a b  = cgx $ last $ cgSolve' sym 1E-4 1E-3 n a b 0
   where
     n = max 10 (round $ sqrt (fromIntegral (dim b) :: Double))
 
 cgSolve'
-  :: CGMat m
-  => Bool      -- ^ symmetric
+  :: Bool      -- ^ symmetric
   -> R         -- ^ relative tolerance for the residual (e.g. 1E-4)
   -> R         -- ^ relative tolerance for δx (e.g. 1E-3)
   -> Int       -- ^ maximum number of iterations
-  -> m         -- ^ coefficient matrix
+  -> GMatrix   -- ^ coefficient matrix
   -> V         -- ^ initial solution
   -> V         -- ^ right-hand side
   -> [CGState] -- ^ solution
 cgSolve' sym er es n a b x = take n $ conjugrad sym a b x er es
 
+
+--------------------------------------------------------------------------------
+
+instance Testable GMatrix
+  where
+    checkT _ = (ok,info)
+      where
+        sma = convo2 20 3
+        x1 = vect [1..20]
+        x2 = vect [1..40]
+        sm = mkSparse sma
+        dm = toDense sma
+
+        s1 = sm !#> x1
+        d1 = dm #> x1
+
+        s2 = tr sm !#> x2
+        d2 = tr dm #> x2
+
+        sdia = mkDiagR 40 20 (vect [1..10])
+        s3 =    sdia !#> x1
+        s4 = tr sdia !#> x2
+        ddia = diagRect 0 (vect [1..10])  40 20
+        d3 = ddia #> x1
+        d4 = tr ddia #> x2
+
+        v = testb 40
+        s5 = cgSolve False sm v
+        d5 = denseSolve dm v
+
+        info = do
+            print sm
+            disp (toDense sma)
+            print s1; print d1
+            print s2; print d2
+            print s3; print d3
+            print s4; print d4
+            print s5; print d5
+            print $ relativeError Infinity s5 d5
+
+        ok = s1==d1
+          && s2==d2
+          && s3==d3
+          && s4==d4
+          && relativeError Infinity s5 d5 < 1E-10
+
+        disp = putStr . dispf 2
+
+        vect = fromList :: [Double] -> Vector Double
+
+        convomat :: Int -> Int -> AssocMatrix
+        convomat n k = [ ((i,j `mod` n),1) | i<-[0..n-1], j <- [i..i+k-1]]
+
+        convo2 :: Int -> Int -> AssocMatrix
+        convo2 n k = m1 ++ m2
+          where
+            m1 = convomat n k
+            m2 = map (((+n) *** id) *** id) m1
+            
+        testb n = vect $ take n $ cycle ([0..10]++[9,8..1])
+        
+        denseSolve a = flatten . linearSolveLS a . asColumn
+
+        -- mkDiag v = mkDiagR (dim v) (dim v) v
+
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util/Static.hs b/packages/base/src/Numeric/LinearAlgebra/Util/Static.hs
deleted file mode 100644
index a3f8eb0..0000000
--- a/packages/base/src/Numeric/LinearAlgebra/Util/Static.hs
+++ /dev/null
@@ -1,70 +0,0 @@
-{-# LANGUAGE DataKinds #-}
-{-# LANGUAGE KindSignatures #-}
-{-# LANGUAGE GeneralizedNewtypeDeriving #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE FlexibleContexts #-}
-{-# LANGUAGE ScopedTypeVariables #-}
-{-# LANGUAGE EmptyDataDecls #-}
-{-# LANGUAGE Rank2Types #-}
-{-# LANGUAGE FlexibleInstances #-}
-{-# LANGUAGE TypeOperators #-}
-
-module Numeric.LinearAlgebra.Util.Static(
-    Static (ddata),
-    R,
-    vect0, sScalar, vect2, vect3, (&)
-) where
-
-
-import GHC.TypeLits
-import Data.Packed.Numeric
-import Numeric.Vector()
-import Numeric.LinearAlgebra.Util(Numeric,ℝ)
-
-lift1F :: (Vector t -> Vector t) -> Static n (Vector t) -> Static n (Vector t)
-lift1F f (Static v) = Static (f v)
-
-lift2F :: (Vector t -> Vector t -> Vector t) -> Static n (Vector t) -> Static n (Vector t) -> Static n (Vector t)
-lift2F f (Static u) (Static v) = Static (f u v)
-
-newtype Static (n :: Nat) t = Static { ddata :: t } deriving Show
-
-type R n = Static n (Vector ℝ)
-
-
-infixl 4 &
-(&) :: R n -> ℝ -> R (n+1)
-Static v & x = Static (vjoin [v, scalar x])
-
-vect0 :: R 0
-vect0 = Static (fromList[])
-
-sScalar :: ℝ -> R 1
-sScalar = Static . scalar
-
-
-vect2 :: ℝ -> ℝ -> R 2
-vect2 x1 x2 = Static (fromList [x1,x2])
-
-vect3 :: ℝ -> ℝ -> ℝ -> R 3
-vect3 x1 x2 x3 = Static (fromList [x1,x2,x3])
-
-
-
-
-
-
-instance forall n t . (KnownNat n, Num (Vector t), Numeric t )=> Num (Static n (Vector t))
-  where
-    (+) = lift2F add
-    (*) = lift2F mul
-    (-) = lift2F sub
-    abs = lift1F abs
-    signum = lift1F signum
-    negate = lift1F (scale (-1))
-    fromInteger x = Static (konst (fromInteger x) d)
-      where
-        d = fromIntegral . natVal $ (undefined :: Proxy n)
-
-data Proxy :: Nat -> *
-