Skip to content

Commit

Permalink
WIP
Browse files Browse the repository at this point in the history
  • Loading branch information
@qnikst
qnikst committed Aug 9, 2017
1 parent 3863bc8 commit 57c6870
Show file tree
Hide file tree
Showing 2 changed files with 74 additions and 173 deletions.
29 changes: 17 additions & 12 deletions src/Math/Integrators/RK.hs
View file
@@ -1,4 +1,4 @@
{-# LANGUAGE QuasiQuotes, FlexibleContexts #-}
{-# LANGUAGE TemplateHaskell, FlexibleContexts #-}
{-# OPTIONS_GHC -Wwarn #-}
-- | Runge-Kutta module
-- TODO: add description and history notes
Expand All @@ -17,23 +17,28 @@ module Math.Integrators.RK

-- import Linear

-- import Math.Integrators.RK.Template
import Math.Integrators.RK.Template
-- import Math.Integrators.RK.Types
-- import Math.Integrators.Internal
-- import Math.Integrators.Implicit

-- | Runge-Kutta
-- @
-- 0 |
-- 0.5 | 0.5
-- 0.5 | 0 & 0.5
-- 1 | 0 & 0 & 1
-- - - + - - - - - - - -
-- | 1/6 & 2/6 & 2/6 & 1/6
-- @
rk45 = $(generateRK [[0::Double]
,[0.5, 0.5]
,[0.5, 0 , 0.5]
,[1 , 0 , 0 , 1]
, [ 1/6,2/6,2/6,1/6]
])

{-
rk45 :: (VectorSpace a, Floating (Scalar a)) => (Double -> a -> a) -> Integrator (Double,a)
rk45 = [qrk|
0 |
0.5 | 0.5
0.5 | 0 & 0.5
1 | 0 & 0 & 1
- - + - - - - - - - -
| 1/6 & 2/6 & 2/6 & 1/6
|]
rk46 :: (VectorSpace a, Floating (Scalar a)) => (Double -> a -> a) -> Integrator (Double,a)
rk46 = [qrk|
0 |
Expand Down
218 changes: 57 additions & 161 deletions src/Math/Integrators/RK/Template.hs
View file
@@ -1,169 +1,65 @@
{-# LANGUAGE TemplateHaskell #-}
{-# OPTIONS_GHC -Wwarn #-} -- this module will be removed in future versions
module Math.Integrators.RK.Template
where

import Data.Maybe

import Language.Haskell.TH.Quote
import Language.Haskell.TH

import Control.Monad

import Math.Integrators.RK.Types
import Math.Integrators.RK.Internal
import Math.Integrators.RK.Parser


qrk :: QuasiQuoter
qrk = QuasiQuoter {quoteExp = x}
-- Let \((b_i)\),\((a_{ij},(i,j=1,\ldots,s)\)\) be real numbers and let
-- \((c=\sum_{j=1}a_{ij})\). An $s$-stage Runge-Kutta method is given by
--
-- \(\(k_i = f(t_0+ c_ih, y_0+h\sum_{j=1}^s a_{ij} k_j), i=1,\ldots,s\)\)
-- \(\(y_1 = y_0 + h\sum_{i=1}^s b_ik_i\)\)
--
-- Te coefficients are usually displayed as follows:
--
-- c_1 | a_11 .. a_{1s}
-- ... | ... ...
-- c_s | a_s1 ... a_{ss}
-- ----+----------------
-- | b_1 ... b_s
--
-- This module provide a template haskell based approach for method generation
-- all elements takes a table as a list:
--
-- [[c1, a_11 ... a_1s]
-- , ...
-- ,[c_s, a_s1 ... a_ss]
-- , [ b_1 ... b_s]
-- ]
--
module Math.Integrators.RK.Template
where
x s = rk $! readMatrixTable s

-----------------------------------------------------------
-- List of helpers
jv = Just . VarE
jv' = Just . varE
jld = Just . LitE. {-DoublePrimL-} RationalL . toRational
ld = LitE . RationalL . toRational
plus = VarE $! mkName "+"
vplus = VarE $! mkName "^+^"
vmult = VarE $! mkName "*^"
mult = VarE $! mkName "*"
ld' = litE . RationalL . toRational
plus' = varE $! mkName "+"
vplus'= varE $! mkName "^+^"
vmult'= varE $! mkName "*^"
mult' = varE $! mkName "*"
vN s = varE (mkName s)

foldOp op = foldl1 (\x y -> infixE (Just x) op (Just y))

realToFracN = varE (mkName "realToFrac")
zeroVN = varE (mkName "zeroV")
f = mkName "f"
t = mkName "t"
h = mkName "h"
y = mkName "y"
tpy = mkName "tpy"

rk :: [MExp] -> Q Exp
rk mExp = do
let (ab,_:c:[]) = break (==Delimeter) mExp
lenA = length ab
kn <- forM [1..lenA] (\_ -> newName "k")
let kvv = zip kn ab
ks' <- forM kvv $ \(k,r) -> do
t <- rkt1 r kn
return $ ValD (VarP k) (NormalB t) []
y' <- rkt2 c kn
if isExplicit mExp
then return $ LamE [VarP f, VarP h, TupP [VarP t,VarP y]] $ LetE ks' y'
else irk mExp
where
rkt1 (Row (Just c,ls)) ks = do
let ft = infixE (jv' t) plus' (Just $ infixE (Just $ ld' c) mult' (jv' h))
st = if null ls
then (varE y)
else
infixE (jv' y) vplus'
(Just $ infixE (Just $ appE realToFracN (varE h)) vmult'
(Just $ foldOp vplus' $
zipWith (\k l -> infixE (Just $ appE realToFracN (ld' l)) vmult' (jv' k))
ks
ls
)
)
appE (appE (varE f) ft) st
rkt2 (Row (_,ls)) ks = do
tupE [ infixE (jv' t) plus' (jv' h)
, infixE (jv' y) vplus'
(Just $ infixE (Just $ appE realToFracN (varE h)) vmult'
(Just $ foldOp vplus' $
zipWith (\k l -> infixE (Just $ appE realToFracN (ld' l)) vmult' (jv' k)) ks ls
)
)
]


test = [Row (Just 1,[2,3]),Row (Just 4,[5,6]),Delimeter, Row (Nothing, [7,8])]

irk mExp = do
fpoint' <- fpoint mExp
lamE [varP tpy, varP f, varP h, tupP [varP t,varP y]] $
caseE (varE tpy) [match (conP (mkName "Math.Integrators.RK.Types.FixedPoint") [varP $ mkName "breakRule"])
(normalB (fpointRun mExp))
fpoint'
,match (conP (mkName "Math.Integrators.RK.Types.NewtonIteration") []) (normalB (varE $ mkName "undefined")) []
]

fpointRun mExp = do
let (ab,_:(Row (_,ls)):[]) = break (==Delimeter) mExp
lenA = length ab
zs <- forM [1..lenA] (\_ -> newName "z")
letE [valD (tupP $ map varP zs)
(normalB $
appE
(appE
(appE
(varE $! mkName "Math.Integrators.Implicit.fixedPointSolver")
(varE $! mkName "method")
)
(varE $ mkName "breakRule")
)
(tupE $ replicate lenA zeroVN) {- TODO: give avaliability to user -}
)
[]]
(appE (varE (mkName "solution")) (tupE $ map varE zs))
import Language.Haskell.TH.Syntax
import Language.Haskell.TH

fpoint mExp = do
let (ab,_:(Row (_,ls)):[]) = break (==Delimeter) mExp
lenA = length ab
zs <- forM [1..lenA] (\_ -> newName "z")
zs' <- forM [1..lenA] (const $ newName "z'")
newVar :: String -> Q (Name, ExpQ, PatQ)
newVar s = do
z <- newName s
return (z,varE z, varP z)

-- | Generate Runge-Kutta method based on the table.
generateRK :: Lift a => [[a]] -> ExpQ
generateRK matrix = do
(_,fE,fP) <- newVar "f"
(_,hE,hP) <- newVar "h"
(_,y0E, y0P) <- newVar "y0"
let makeKas = zipWith go [(0::Int)..] ks
where
go i name = funD name [clause [] (normalB (k i)) []]
k i = [| $fE (t0 + $(mkC (cs !! i)) * $hE) ($y0E + $hE * $(sumOp (as !! i))) |]
y1 = [| $y0E + $hE * $(sumOp bs) |]
[| \ $fP $hP $y0P -> $(letE makeKas [| $y1 |] )
|]
where
-- split the matrix into parameters
as = map (drop 1) $ init matrix
bs = last matrix
cs = concat $ take 1 $ init matrix
s = length matrix - 1
ks = map (\i -> mkName ("k" ++ show i)) [0..(s-1)]
-- Generate $k_i$ according to formula above
mkC c = [| c |]
sumOp aa = internal 0 aa
where internal i [x] = [| x * $(varE $ ks !! i)|]
internal i (x:xs) = [| x * $(varE $ ks !! i) + $(internal (i+1) xs) |]
internal _ [] = [| 0 |]

return $
[ funD (mkName "method") [clause [tupP $ map varP zs] (normalB $ letE (map (topRow zs) $! zip zs' ab) (tupE $ map varE zs')) [] ]
, funD (mkName "solution") [clause [tupP $ map varP zs] (normalB $ solutionRow zs ls) []]
]
where
topRow zs (x,(Row (Just c,ls))) =
valD (varP x)
(normalB $ infixE (Just $ appE realToFracN (varE h))
vmult'
(Just $ foldOp vplus' $
zipWith (\z l -> infixE (Just $ appE realToFracN (ld' l))
vmult'
(Just $ appE (appE (varE f)
(infixE (jv' t)
plus'
(Just $ infixE (Just $ appE realToFracN $ varE h) mult' (Just $ ld' c))
)
)
(infixE (jv' y) vplus' (jv' z))
)
)
zs
ls
)
)
[]
topRow _ _ = error "not a row"
solutionRow zs ls =
(tupE [infixE (jv' t) plus' (jv' h)
,infixE (jv' y)
vplus'
(Just $ infixE (Just $ appE realToFracN $ varE h)
vmult'
(Just $ foldOp vplus'
(zipWith (\z b -> infixE (Just $ appE realToFracN
(ld' b))
vmult'
(jv' z))
zs
ls
)
)
)
]
)

0 comments on commit 57c6870

Please sign in to comment.