8 files changed, 375 insertions, 0 deletions

diff --git a/examples/experiments/Static.hs b/examples/experiments/Static.hs
new file mode 100644
index 0000000..c0cfcce
--- /dev/null
+++ b/examples/experiments/Static.hs
@@ -0,0 +1,123 @@
+{-# OPTIONS_GHC -fglasgow-exts -fth -fallow-overlapping-instances -fallow-undecidable-instances #-}
+
+module Static where
+
+import Language.Haskell.TH
+import Numeric.LinearAlgebra
+import Foreign
+import Language.Haskell.TH.Syntax
+import Data.Packed.Internal(Vector(..),Matrix(..))
+
+instance Lift Double where
+  lift x = return (LitE (RationalL (toRational x)))
+
+instance Lift (Ptr Double) where
+    lift p = [e| p |]
+
+instance Lift (ForeignPtr Double) where
+    lift p = [e| p |]
+
+instance (Lift a, Storable a, Lift (Ptr a), Lift (ForeignPtr a)) => Lift (Vector a ) where
+    lift (V n fp p) = [e| V $(lift n) $(lift fp) $(lift p) |]
+
+instance (Lift (Vector a)) => Lift (Matrix a) where
+    lift (MC r c v vt) = [e| MC $(lift r) $(lift c) $(lift v) $(lift vt) |]
+    lift (MF r c v vt) = [e| MF $(lift r) $(lift c) $(lift v) $(lift vt) |]
+
+tdim :: Int -> ExpQ
+tdim 0 = [| Z |]
+tdim n = [| S $(tdim (n-1)) |]
+
+
+data Z = Z deriving Show
+data S a = S a deriving Show
+
+class Dim a
+
+instance Dim Z
+instance Dim a => Dim (S a)
+
+class Sum a b c | a b -> c -- , a c -> b, b c -> a
+
+instance Sum Z a a
+instance Sum a Z a
+instance Sum a b c => Sum a (S b) (S c)
+
+newtype SVec d t = SVec (Vector t) deriving Show
+newtype SMat r c t = SMat (Matrix t) deriving Show
+
+createl :: d -> [Double] -> SVec d Double
+createl d l = SVec (fromList l)
+
+createv :: Storable t => d -> Vector t -> SVec d t
+createv d v = SVec v
+
+--vec'' v = [|createv ($(tdim (dim v))) v|]
+
+vec' :: [Double] -> ExpQ
+vec' d = [| createl ($(tdim (length d))) d |]
+
+
+createm :: (Dim r, Dim c) => r -> c -> (Matrix Double) -> SMat r c Double
+createm _ _ m = SMat m
+
+createml :: (Dim r, Dim c) => r -> c -> Int -> Int -> [Double] -> SMat r c Double
+createml _ _ r c l = SMat ((r><c) l)
+
+mat :: Int -> Int -> [Double] -> ExpQ
+mat r c l = [| createml ($(tdim r)) ($(tdim c)) r c l |]
+
+vec :: [Double] -> ExpQ
+vec d = mat (length d) 1 d
+
+
+mat' :: Matrix Double -> ExpQ
+mat' m = [| createm ($(tdim (rows m))) ($(tdim (cols m))) m |]
+
+covec :: [Double] -> ExpQ
+covec d = mat 1 (length d) d
+
+scalar :: SMat (S Z) (S Z) Double -> Double
+scalar (SMat m) = flatten m @> 0
+
+v = fromList [1..5] :: Vector Double
+l = [1,1.5..5::Double]
+
+k = [11..30::Int]
+
+rawv (SVec v) = v
+raw (SMat m) = m
+
+liftStatic :: (Matrix a -> Matrix b -> Matrix c) -> SMat dr dc a -> SMat dr dc b -> SMat dr dc c
+liftStatic f a b = SMat (f (raw a) (raw b))
+
+a |+| b = liftStatic (+) a b
+
+prod :: SMat r k Double -> SMat k c Double -> SMat r c Double
+prod a b = SMat (raw a <> raw b)
+
+strans :: SMat r c Double -> SMat c r Double
+strans = SMat . trans . raw
+
+sdot a b = scalar (prod a b)
+
+jv :: (Field t, Sum r1 r2 r3) => SMat r1 c t -> SMat r2 c t -> SMat r3 c t
+jv a b = SMat ((raw a) <-> (raw b))
+
+-- curiously, we cannot easily fold jv because the matrics are not of the same type.
+
+jh a b = strans (jv (strans a) (strans b))
+
+
+homog :: (Field t) => SMat r c t -> SMat (S r) c t
+homog m = SMat (raw m <-> constant 1 (cols (raw m)))
+
+inhomog :: (Linear Vector t) => SMat (S (S r)) c t -> SMat r c t
+inhomog (SMat m) = SMat (sm <> d) where
+    sm = takeRows r' m
+    d = diag $ 1 / (flatten $ dropRows r' m)
+    r' = rows m -1
+
+
+ht t vs = inhomog (t `prod` homog vs)
+
diff --git a/examples/experiments/error.hs b/examples/experiments/error.hs
new file mode 100644
index 0000000..16305dc
--- /dev/null
+++ b/examples/experiments/error.hs
@@ -0,0 +1,16 @@
+import Numeric.GSL
+import Prelude hiding (catch)
+import Control.Exception
+
+test x = catch
+       (print x)
+       (\e -> putStrLn $ "captured ["++ show e++"]")
+
+main = do
+    setErrorHandlerOff
+
+    test $ log_e (-1)
+    test $ 5 + (fst.exp_e) 1000
+    test $ bessel_zero_Jnu_e (-0.3) 2
+
+    putStrLn "Bye"
\ No newline at end of file
diff --git a/examples/experiments/listlike.hs b/examples/experiments/listlike.hs
new file mode 100644
index 0000000..43216cb
--- /dev/null
+++ b/examples/experiments/listlike.hs
@@ -0,0 +1,32 @@
+{-# OPTIONS_GHC -fglasgow-exts #-}
+
+import qualified Data.ListLike as LL
+import Numeric.LinearAlgebra
+import Data.Monoid
+import Data.Packed.Internal.Vector
+import Foreign
+
+instance (Storable a) => Monoid (Vector a) where
+    mempty = V { dim = 0, fptr = undefined, ptr = undefined }
+    mappend a b = mconcat [a,b]
+    mconcat = j . filter ((>0).dim)
+        where j [] = mempty
+              j l  = join l
+
+instance Storable a => LL.FoldableLL (Vector a) a where
+    foldl f x v = foldl f x (toList v)
+    foldr f x v = foldr f x (toList v)
+
+instance Storable a => LL.ListLike (Vector a) a where
+    singleton a = fromList [a]
+    head a = a @> 0
+    tail a | dim a == 1 = mempty
+           | otherwise  = subVector 1 (dim a -1) a
+    genericLength = fromIntegral.dim
+
+
+v k = fromList [1..k] :: Vector Double
+
+f k = k+(3::Double)
+
+main = print $ (LL.map f [1..5] :: Vector Double)
\ No newline at end of file
diff --git a/examples/experiments/parallel.hs b/examples/experiments/parallel.hs
new file mode 100644
index 0000000..7256fb6
--- /dev/null
+++ b/examples/experiments/parallel.hs
@@ -0,0 +1,28 @@
+import System(getArgs)
+import Numeric.LinearAlgebra
+import Control.Parallel.Strategies
+import System.Time
+
+inParallel = parMap rwhnf id
+
+parMul x y = fromBlocks [inParallel[x <> y1, x <> y2]]
+    where p = cols y `div` 2
+          (y1,y2) = splitColumnsAt p y
+
+main = do
+    n <- (read . head) `fmap` getArgs
+    let m = ident n :: Matrix Double
+    time $ print $ vectorMax $ takeDiag $ parMul m m
+    time $ print $ vectorMax $ takeDiag $ m <> m
+
+a = (2><3) [1..6::Double]
+b = (3><4) [1..12::Double]
+
+splitRowsAt p m    = (takeRows p m, dropRows p m)
+splitColumnsAt p m = (takeColumns p m, dropColumns p m)
+
+time act = do
+    t0 <- getClockTime
+    act
+    t1 <- getClockTime
+    print $ tdSec $ normalizeTimeDiff $ diffClockTimes t1 t0
diff --git a/examples/experiments/speed.hs b/examples/experiments/speed.hs
new file mode 100644
index 0000000..22d7220
--- /dev/null
+++ b/examples/experiments/speed.hs
@@ -0,0 +1,90 @@
+{- speed tests
+
+GNU-Octave (see speed.m in this folder):
+
+./speed.m
+  -0.017877255967426   0.000000000000000  -0.999840189089781
+   0.000000000000000   1.000000000000000   0.000000000000000
+   0.999840189089763   0.000000000000000  -0.017877255967417
+9.69 seconds
+
+Mathematica:
+
+rot[a_]:={{ Cos[a], 0, Sin[a]},
+          { 0,      1,      0},
+          { -Sin[a],0,Cos[a]}}//N
+
+test := Timing[
+    n = 100000;
+    angles = Range[0.,n-1]/(n-1);
+    Fold[Dot,IdentityMatrix[3],rot/@angles] // MatrixForm
+]
+
+2.08013 Second
+     {{\(-0.017877255967432837`\), "0.`", \(-0.9998401890898042`\)},
+      {"0.`", "1.`", "0.`"},
+      {"0.9998401890898042`", "0.`", \(-0.017877255967432837`\)}}
+
+$ ghc --make -O speed
+
+$ ./speed 5 100000 1
+(3><3)
+ [ -1.7877255967425523e-2, 0.0,    -0.9998401890897632
+ ,                    0.0, 1.0,                    0.0
+ ,      0.999840189089781, 0.0, -1.7877255967416586e-2 ]
+0.33 CPU seconds
+
+cos 50000 = -0.0178772559665563
+sin 50000 = -0.999840189089790
+
+-}
+
+import Numeric.LinearAlgebra
+import System
+import Data.List(foldl1')
+import System.CPUTime
+import Text.Printf
+import Debug.Trace
+
+debug x = trace (show x) x
+
+timing act = do
+    t0 <- getCPUTime
+    act
+    t1 <- getCPUTime
+    printf "%.2f CPU seconds\n" $ (fromIntegral ((t1 - t0) `div` (10^10)) / 100 :: Double) :: IO ()
+
+op a b = trans $ (trans a) <> (trans b)
+
+op2 a b = trans $ (trans a) + (trans b)
+
+rot' :: Double -> Matrix Double
+rot' a = matrix [[ c,0,s],
+                 [ 0,1,0],
+                 [-s,0,c]]
+    where c = cos a
+          s = sin a
+          matrix = fromLists
+
+rot :: Double -> Matrix Double
+rot a = (3><3) [ c,0,s
+               , 0,1,0
+               ,-s,0,c ]
+    where c = cos a
+          s = sin a
+
+fun n r = foldl1' (<>) (map r angles)
+    where angles = toList $ linspace n (0,1)
+
+main = do
+    args <- getArgs
+    let [p,n,d] = map read args
+    let ms = replicate n ((ident d :: Matrix Double))
+    let mz = replicate n (diag (constant (0::Double) d))
+    timing $ case p of
+        0 -> print $ foldl1' (<>) ms
+        1 -> print $ foldl1' (<>) (map trans ms)
+        2 -> print $ foldl1' op ms
+        3 -> print $ foldl1' op2 mz
+        4 -> print $ fun n rot'
+        5 -> print $ fun n rot
diff --git a/examples/experiments/speed.m b/examples/experiments/speed.m
new file mode 100644
index 0000000..2f41665
--- /dev/null
+++ b/examples/experiments/speed.m
@@ -0,0 +1,21 @@
+#! /usr/bin/octave -qf
+1; % measuring Octave computing times
+
+function r = rot(a)
+    c = cos(a);
+    s = sin(a);
+    r = [ c , 0, s;
+          0,  1, 0;
+         -s,  0, c];
+end
+
+t0=time();
+x = linspace(0,1,1E5);
+ac = eye(3);
+for a = x
+    ac = ac*rot(a);
+end
+
+format long
+disp(ac);
+printf("%.2f seconds\n",time()-t0)
diff --git a/examples/experiments/testmnist.m b/examples/experiments/testmnist.m
new file mode 100644
index 0000000..38625a7
--- /dev/null
+++ b/examples/experiments/testmnist.m
@@ -0,0 +1,29 @@
+#! /usr/bin/octave -qf
+% measuring Octave computing times
+
+t0=time();
+load mnist.txt
+disp("load");
+disp(time()-t0)
+
+
+x = mnist(:,1:784);
+d = mnist(:,785);
+
+
+t0=time();
+xc = x - repmat(mean(x),rows(x),1);
+disp("x - repmat(mean(x),rows(x),1)");
+disp(time()-t0)
+
+t0=time();
+mc = (xc'*xc)/rows(x);
+disp("(xc'*xc)/rows(x)");
+disp(time()-t0)
+
+t0=time();
+[v,l]=eig(mc);
+disp("eig");
+disp(time()-t0)
+
+disp(flipud(diag(l))(1:10));
\ No newline at end of file
diff --git a/examples/experiments/useStatic.hs b/examples/experiments/useStatic.hs
new file mode 100644
index 0000000..619af8f
--- /dev/null
+++ b/examples/experiments/useStatic.hs
@@ -0,0 +1,36 @@
+{-# OPTIONS -fno-monomorphism-restriction #-}
+
+import Static
+import Numeric.LinearAlgebra
+
+
+x = $(vec [1,2])
+
+y = $(vec [5,7])
+
+z a = vec [a,a]
+
+w = $(vec [1,2,3])
+
+cx = $(covec [1,2,3])
+
+
+t3 = $(tdim 3)
+
+crm33 = createml t3 t3 3 3
+
+rot a = crm33 [a,0,0,0,a,0,0,0,1]
+
+--q = x |+| y |+| $(z 5)
+
+m = $(mat 2 3 [1..6])
+
+n = $(mat 3 5 [1..15])
+
+infixl 7 <*>
+(<*>) = prod
+
+r1 = m <*> n
+r2 = strans (strans n <*> strans m)
+
+--r' = prod n m