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