summaryrefslogtreecommitdiff
path: root/examples/experiments/Static.hs
diff options
context:
space:
mode:
Diffstat (limited to 'examples/experiments/Static.hs')
-rw-r--r--examples/experiments/Static.hs123
1 files changed, 123 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 @@
1{-# OPTIONS_GHC -fglasgow-exts -fth -fallow-overlapping-instances -fallow-undecidable-instances #-}
2
3module Static where
4
5import Language.Haskell.TH
6import Numeric.LinearAlgebra
7import Foreign
8import Language.Haskell.TH.Syntax
9import Data.Packed.Internal(Vector(..),Matrix(..))
10
11instance Lift Double where
12 lift x = return (LitE (RationalL (toRational x)))
13
14instance Lift (Ptr Double) where
15 lift p = [e| p |]
16
17instance Lift (ForeignPtr Double) where
18 lift p = [e| p |]
19
20instance (Lift a, Storable a, Lift (Ptr a), Lift (ForeignPtr a)) => Lift (Vector a ) where
21 lift (V n fp p) = [e| V $(lift n) $(lift fp) $(lift p) |]
22
23instance (Lift (Vector a)) => Lift (Matrix a) where
24 lift (MC r c v vt) = [e| MC $(lift r) $(lift c) $(lift v) $(lift vt) |]
25 lift (MF r c v vt) = [e| MF $(lift r) $(lift c) $(lift v) $(lift vt) |]
26
27tdim :: Int -> ExpQ
28tdim 0 = [| Z |]
29tdim n = [| S $(tdim (n-1)) |]
30
31
32data Z = Z deriving Show
33data S a = S a deriving Show
34
35class Dim a
36
37instance Dim Z
38instance Dim a => Dim (S a)
39
40class Sum a b c | a b -> c -- , a c -> b, b c -> a
41
42instance Sum Z a a
43instance Sum a Z a
44instance Sum a b c => Sum a (S b) (S c)
45
46newtype SVec d t = SVec (Vector t) deriving Show
47newtype SMat r c t = SMat (Matrix t) deriving Show
48
49createl :: d -> [Double] -> SVec d Double
50createl d l = SVec (fromList l)
51
52createv :: Storable t => d -> Vector t -> SVec d t
53createv d v = SVec v
54
55--vec'' v = [|createv ($(tdim (dim v))) v|]
56
57vec' :: [Double] -> ExpQ
58vec' d = [| createl ($(tdim (length d))) d |]
59
60
61createm :: (Dim r, Dim c) => r -> c -> (Matrix Double) -> SMat r c Double
62createm _ _ m = SMat m
63
64createml :: (Dim r, Dim c) => r -> c -> Int -> Int -> [Double] -> SMat r c Double
65createml _ _ r c l = SMat ((r><c) l)
66
67mat :: Int -> Int -> [Double] -> ExpQ
68mat r c l = [| createml ($(tdim r)) ($(tdim c)) r c l |]
69
70vec :: [Double] -> ExpQ
71vec d = mat (length d) 1 d
72
73
74mat' :: Matrix Double -> ExpQ
75mat' m = [| createm ($(tdim (rows m))) ($(tdim (cols m))) m |]
76
77covec :: [Double] -> ExpQ
78covec d = mat 1 (length d) d
79
80scalar :: SMat (S Z) (S Z) Double -> Double
81scalar (SMat m) = flatten m @> 0
82
83v = fromList [1..5] :: Vector Double
84l = [1,1.5..5::Double]
85
86k = [11..30::Int]
87
88rawv (SVec v) = v
89raw (SMat m) = m
90
91liftStatic :: (Matrix a -> Matrix b -> Matrix c) -> SMat dr dc a -> SMat dr dc b -> SMat dr dc c
92liftStatic f a b = SMat (f (raw a) (raw b))
93
94a |+| b = liftStatic (+) a b
95
96prod :: SMat r k Double -> SMat k c Double -> SMat r c Double
97prod a b = SMat (raw a <> raw b)
98
99strans :: SMat r c Double -> SMat c r Double
100strans = SMat . trans . raw
101
102sdot a b = scalar (prod a b)
103
104jv :: (Field t, Sum r1 r2 r3) => SMat r1 c t -> SMat r2 c t -> SMat r3 c t
105jv a b = SMat ((raw a) <-> (raw b))
106
107-- curiously, we cannot easily fold jv because the matrics are not of the same type.
108
109jh a b = strans (jv (strans a) (strans b))
110
111
112homog :: (Field t) => SMat r c t -> SMat (S r) c t
113homog m = SMat (raw m <-> constant 1 (cols (raw m)))
114
115inhomog :: (Linear Vector t) => SMat (S (S r)) c t -> SMat r c t
116inhomog (SMat m) = SMat (sm <> d) where
117 sm = takeRows r' m
118 d = diag $ 1 / (flatten $ dropRows r' m)
119 r' = rows m -1
120
121
122ht t vs = inhomog (t `prod` homog vs)
123