dvda (empty) → 0.1
raw patch · 16 files changed
+1749/−0 lines, 16 filesdep +basedep +containersdep +deepseqsetup-changed
Dependencies added: base, containers, deepseq, directory, fgl, graphviz, hashable, mtl, plugins, process, repa, text, transformers, unordered-containers, vector
Files
- Dvda.hs +46/−0
- Dvda/BinUn.hs +143/−0
- Dvda/Config.hs +116/−0
- Dvda/Dot.hs +67/−0
- Dvda/Dual.hs +60/−0
- Dvda/Examples.hs +95/−0
- Dvda/Expr.hs +263/−0
- Dvda/GExpr.hs +87/−0
- Dvda/Graph.hs +114/−0
- Dvda/HSBuilder.hs +116/−0
- Dvda/HSSyntax.hs +148/−0
- Dvda/HomoDim.hs +36/−0
- Dvda/SymMonad.hs +324/−0
- LICENSE +30/−0
- Setup.hs +2/−0
- dvda.cabal +102/−0
+ Dvda.hs view
@@ -0,0 +1,46 @@+{- |+ Module : Dvda+ Description : Top level module++ This is the top level module which exports the API+ -}++{-# OPTIONS_GHC -Wall #-}++module Dvda ( -- * primitives+ sym+ , vsym+ , msym+ , vec+ , mat+ -- * operations+ , scale+ , dot+ , diff+ -- * symbolic expression type+ , Expr+ -- * construct FunGraphs+ , FunGraph+ , makeFunGraph+ , runFunGraph+ , inputs_+ , outputs_+ , node+ -- * show/summarize FunGraphs+ , funGraphSummary+ , funGraphSummary'+ , showCollisions+ , previewGraph+ -- * compile and link function+ , buildHSFunction+ -- * Heterogenous inputs/outputs+ , (:*)(..)+ , Exprs+ ) where++import Dvda.Expr+import Dvda.Graph+import Dvda.HSBuilder+import Dvda.SymMonad++
+ Dvda/BinUn.hs view
@@ -0,0 +1,143 @@+{-# OPTIONS_GHC -Wall #-}++module Dvda.BinUn ( BinOp(..)+ , UnOp(..)+ , showBinary+ , showUnary+ , applyUnary+ , applyBinary+ , unaryDeriv+ , binaryDeriv+ , isCommutative+ ) where++import Data.Hashable ( Hashable, hash )++import Dvda.Dual ( Dual(..), dualPerturbation )++data UnOp = Abs+ | Neg+ | Signum+ | Exp+ | Sqrt+ | Log+ | Sin+ | Cos+ | Tan+ | ASin+ | ACos+ | ATan+ | Tanh+ | Sinh+ | Cosh+ | ATanh+ | ASinh+ | ACosh deriving (Eq, Show)++data BinOp = Add+ | Sub+ | Mul+ | Div+ | Pow+ | LogBase deriving (Eq, Show)++instance Hashable UnOp where+ hash Abs = 0+ hash Neg = 1+ hash Signum = 2+ hash Exp = 3+ hash Sqrt = 4+ hash Log = 5+ hash Sin = 6+ hash Cos = 7+ hash Tan = 8+ hash ASin = 9+ hash ACos = 10+ hash ATan = 11+ hash Tanh = 12+ hash Sinh = 13+ hash Cosh = 14+ hash ATanh = 15+ hash ASinh = 16+ hash ACosh = 17++instance Hashable BinOp where+ hash Add = 18+ hash Sub = 19+ hash Mul = 20+ hash Div = 21+ hash Pow = 22+ hash LogBase = 23+ +showUnary :: Show a => a -> UnOp -> String+showUnary x Abs = '|': show x ++ "|"+showUnary x Neg = '-':paren x+showUnary x Signum = "signum"++paren x+showUnary x Exp = "exp"++paren x+showUnary x Sqrt = "sqrt"++paren x+showUnary x Log = "log"++paren x+showUnary x Sin = "sin"++paren x+showUnary x Cos = "cos"++paren x+showUnary x Tan = "tan"++paren x+showUnary x ASin = "asin"++paren x+showUnary x ACos = "acos"++paren x+showUnary x ATan = "atan"++paren x+showUnary x Sinh = "sinh"++paren x+showUnary x Cosh = "cosh"++paren x+showUnary x Tanh = "tanh"++paren x+showUnary x ASinh = "asinh"++paren x+showUnary x ATanh = "atanh"++paren x+showUnary x ACosh = "acosh"++paren x++applyUnary :: Floating a => UnOp -> a -> a+applyUnary Abs = abs+applyUnary Neg = negate+applyUnary Signum = signum+applyUnary Exp = exp+applyUnary Sqrt = sqrt+applyUnary Log = log+applyUnary Sin = sin+applyUnary Cos = cos+applyUnary Tan = tan+applyUnary ASin = asin+applyUnary ACos = acos+applyUnary ATan = atan+applyUnary Sinh = sinh+applyUnary Cosh = cosh+applyUnary Tanh = tanh+applyUnary ASinh = asinh+applyUnary ATanh = atanh+applyUnary ACosh = acosh++applyBinary :: Floating a => BinOp -> a -> a -> a+applyBinary Add = (+)+applyBinary Sub = (-)+applyBinary Mul = (*)+applyBinary Div = (/)+applyBinary Pow = (**)+applyBinary LogBase = logBase++unaryDeriv :: Floating a => UnOp -> (a,a) -> a+unaryDeriv op (x,x') = dualPerturbation $ applyUnary op (Dual x x')++binaryDeriv :: Floating a => BinOp -> (a,a) -> (a,a) -> a+binaryDeriv op (x,x') (y,y') = dualPerturbation $ applyBinary op (Dual x x') (Dual y y')++showBinary :: BinOp -> String+showBinary Add = "+"+showBinary Sub = "-"+showBinary Mul = "*"+showBinary Div = "/"+showBinary Pow = "**"+showBinary LogBase = "`logbase`"++isCommutative :: BinOp -> Bool+isCommutative Add = True+isCommutative Sub = False+isCommutative Mul = True+isCommutative Div = False+isCommutative Pow = False+isCommutative LogBase = False++paren :: Show a => a -> String+paren x = "( "++show x++" )"
+ Dvda/Config.hs view
@@ -0,0 +1,116 @@+-- Config.hs++{-# OPTIONS_GHC -Wall #-}++module Dvda.Config( -- * directory stuff+ dvdaDir+ , functionDir+ -- * C syntax+ , cType+ , cName+ , nameCSource+ , nameCInclude+ , nameCObject+ , nameCFunction+ -- * Haskell syntax+ , nameHSObject+ , nameHSModule+ , nameHSFunction+ , nameHSSource+ , nameHSVar+ , nameHSConst+ -- * gcc stuff+ , gccString+ , spewGccCall+ , outputNames+ -- * ghc stuff+ , ghcString+ ) where++import System.Directory+import Control.Monad(unless)++-- | what symbolic variable names to use when generating a function+outputNames :: [String]+outputNames = map (\x -> "out"++show x) [(0::Integer)..]++-- | whether to print the gcc call when generating code+spewGccCall :: Bool+spewGccCall = True++-- | return directory to use for temp files+-- | create this directory and print message if it doesn't exist+dvdaDir :: IO FilePath+dvdaDir = do+ dir <- getAppUserDataDirectory "dvda"+ + -- print message if creating directory+ exist <- doesDirectoryExist dir+ unless exist $ putStrLn $ "creating directory \""++dir++"\" for codegen source/objects"++ -- make the directory if missing+ createDirectoryIfMissing True dir+ + return dir+++-- | take in source file and object, return string suitible for calling to compile+gccString :: FilePath -> FilePath -> String+gccString src obj = "gcc -O2 -std=gnu99 -fPIC -shared -Wall -Wextra -Werror " ++ src ++ " -o " ++ obj++ghcString :: FilePath -> FilePath -> String+ghcString src obj = "ghc -c " ++ src ++ " -o " ++ obj+++functionDir :: String -> IO FilePath+functionDir hash = do+ -- dvda directory+ topDir <- dvdaDir+ return (topDir ++ "/" ++ nameCFunction hash)++-- c syntax+-- | type to use when generating c code+cType :: String+cType = "double"++-- | name convention for c variables+cName :: Int -> String+cName k+ | k < 0 = error "cName got negative index" + | otherwise = 't':show k ++nameCSource :: String -> String+nameCSource hash = nameCFunction hash ++ ".c"++nameCInclude :: String -> String+nameCInclude hash = nameCFunction hash ++ ".h"++nameCObject :: String -> String+nameCObject hash = "c_" ++ nameCFunction hash ++ ".o"++nameCFunction :: String -> String+nameCFunction hash = "call_" ++ hash++-- haskell syntax+nameHSVar :: Int -> String+nameHSVar k+ | k < 0 = error "nameHSVar got negative index" + | otherwise = 'v':show k++-- haskell syntax+nameHSConst :: Int -> String+nameHSConst k+ | k < 0 = error "nameHSConst got negative index" + | otherwise = 'c':show k++nameHSFunction :: String -> String+nameHSFunction hash = "call_" ++ hash++nameHSModule :: String -> String+nameHSModule hash = "Call_" ++ hash++nameHSSource :: String -> String+nameHSSource = (++ ".hs") . nameHSModule++nameHSObject :: String -> String+nameHSObject = (++ ".o") . nameHSModule
+ Dvda/Dot.hs view
@@ -0,0 +1,67 @@+{-# OPTIONS_GHC -Wall #-}+{-# Language TypeFamilies #-}+{-# Language MultiParamTypeClasses #-}+{-# Language FlexibleContexts #-}+{-# Language FlexibleInstances #-}++module Dvda.Dot ( Dot(..)+ ) where++import Data.Array.Repa(DIM0,DIM1,DIM2,Z(..),(:.)(..), listOfShape, Shape, shapeOfList)++import Dvda.HomoDim ( HomoDim(..), homoOfShape ) +++class (Shape sh1, Shape sh2, Shape (DotT sh1 sh2)) => Dot sh1 sh2 where+ type DotT sh1 sh2+ dotDims :: sh1 -> sh2 -> DotT sh1 sh2++instance Dot HomoDim HomoDim where+ type DotT HomoDim HomoDim = HomoDim+ dotDims (HomoDim x@[_,_]) (HomoDim y@[_,_]) =+ homoOfShape $ dotDims (shapeOfList x :: DIM2) (shapeOfList y :: DIM2)+ dotDims (HomoDim x@[_,_]) (HomoDim y@[_]) =+ homoOfShape $ dotDims (shapeOfList x :: DIM2) (shapeOfList y :: DIM1)+ dotDims (HomoDim x@[_]) (HomoDim y@[_,_]) =+ homoOfShape $ dotDims (shapeOfList x :: DIM1) (shapeOfList y :: DIM2)+ dotDims (HomoDim x@[_]) (HomoDim y@[_]) =+ homoOfShape $ dotDims (shapeOfList x :: DIM1) (shapeOfList y :: DIM1)+ dotDims x y = error $ "dotDims HomoDim not instanced for " ++ show x ++ " " ++ show y+ ++instance Dot DIM2 DIM2 where -- matrix-matrix+ type DotT DIM2 DIM2 = DIM2+ dotDims sh1 sh2 + | c1 == r2 = Z :. r1 :. c2+ | otherwise = error $ "MM dimension mismatch: " ++ show sh1' ++ ", " ++ show sh2'+ where+ sh1'@[r1,c1] = reverse $ listOfShape sh1+ sh2'@[r2,c2] = reverse $ listOfShape sh2+ +instance Dot DIM1 DIM1 where -- vector-vector+ type DotT DIM1 DIM1 = DIM0+ dotDims sh1 sh2 + | r1 == r2 = Z+ | otherwise = error $ "VV dimension mismatch: " ++ show sh1' ++ ", " ++ show sh2'+ where+ sh1'@[r1] = listOfShape sh1+ sh2'@[r2] = listOfShape sh2++instance Dot DIM2 DIM1 where -- matrix-vector+ type DotT DIM2 DIM1 = DIM1+ dotDims sh1 sh2 + | c1 == r2 = Z :. r1+ | otherwise = error $ "MV dimension mismatch: " ++ show sh1' ++ ", " ++ show sh2'+ where+ sh1'@[r1,c1] = reverse $ listOfShape sh1+ sh2'@[r2] = reverse $ listOfShape sh2++instance Dot DIM1 DIM2 where -- vector-matrix+ type DotT DIM1 DIM2 = DIM1+ dotDims sh1 sh2 + | c1 == r2 = Z :. c2+ | otherwise = error $ "VM dimension mismatch: " ++ show sh1' ++ ", " ++ show sh2'+ where+ sh1'@[c1] = reverse $ listOfShape sh1+ sh2'@[r2,c2] = reverse $ listOfShape sh2+
+ Dvda/Dual.hs view
@@ -0,0 +1,60 @@+{-# OPTIONS_GHC -Wall #-}+{-# Language GADTs #-}+{-# Language FlexibleContexts #-}+{-# Language TypeOperators #-}+{-# Language TypeFamilies #-}++module Dvda.Dual ( Dual(..)+ ) where++import Data.Ratio ( numerator, denominator )++data Dual a = Dual { dualPrimal :: a+ , dualPerturbation :: a+ } deriving (Show, Eq)++instance Num a => Num (Dual a) where+ (Dual x x') * (Dual y y') = Dual (x * y) (x*y' + x'*y)+ (Dual x x') + (Dual y y') = Dual (x + y) (x' + y')+ (Dual x x') - (Dual y y') = Dual (x - y) (x' - y')+ negate (Dual x x') = Dual (-x) (-x')+ abs (Dual x x') = Dual (abs x) (signum x * x')+ signum (Dual x _) = Dual (signum x) 0 -- technically this should be a dirac delta+ fromInteger x = Dual (fromInteger x) 0+ +instance Fractional a => Fractional (Dual a) where+ (Dual x x') / (Dual y y') = Dual (x/y) (x'/y - x/(y*y)*y')+ fromRational x = num/den+ where+ num = fromIntegral $ numerator x+ den = fromIntegral $ denominator x++instance Floating a => Floating (Dual a) where+ pi = Dual pi 0+ + exp (Dual x x') = Dual (exp x) (exp x *x')+ sqrt (Dual x x') = Dual (sqrt x) (x'/(2*sqrt x))+ log (Dual x x') = Dual (log x) (x'/x)+ + (Dual x x')**(Dual y y') = Dual (x**y) $ ( x'*y + x*y'*log x ) * x**(y-1)++ logBase (Dual b b') (Dual e e') = Dual primal pert'+ where+ primal = logBase b e+ pert' = (e'/e - primal*b'/b) / log b+ + sin (Dual x x') = Dual (sin x) $ cos x * x'+ cos (Dual x x') = Dual (cos x) $ -(sin x)*x'+ tan (Dual x x') = Dual (tan x) $ x'/(cos x * cos x)+ + asin (Dual x x') = Dual (asin x) $ x' / sqrt (1 - x*x)+ acos (Dual x x') = Dual (acos x) $ -x' / sqrt (1 - x*x)+ atan (Dual x x') = Dual (atan x) $ x' / (1 + x*x)+ + sinh (Dual x x') = Dual (sinh x) $ cosh x * x'+ cosh (Dual x x') = Dual (cosh x) $ sinh x * x'+ tanh (Dual x x') = Dual (tanh x) $ x'/(cosh x * cosh x)+ + asinh (Dual x x') = Dual (asinh x) $ x'/ sqrt (1 + x*x)+ acosh (Dual x x') = Dual (acosh x) $ x'/( sqrt (x - 1) * sqrt (x + 1) )+ atanh (Dual x x') = Dual (atanh x) $ x'/(1 - x*x)
+ Dvda/Examples.hs view
@@ -0,0 +1,95 @@+{-# OPTIONS_GHC -Wall #-}+{-# Language TypeOperators #-}++module Dvda.Examples ( exampleFun+ , run+ , run'+ , showoff+ ) where++import Control.Monad.State (State)+import Data.Array.Repa (DIM0,DIM1,DIM2)++import Dvda+import Dvda.Graph ( FunGraph(..) )++exampleFun :: State (FunGraph Double (DIM0 :* DIM1 :* DIM2) (DIM2 :* DIM1 :* DIM0)) ()+exampleFun = do+ let x = sym "x"+ y = vsym 5 "y"+ z = msym (3,5) "Z"+ inputs_ (x :* y :* z)+ + z1 <- node $ (scale x z)**3+ z2 <- node $ (dot z y)**2+ z3 <- node $ diff ((x*x/2)**x) x+ + outputs_ (z1 :* z2 :* z3)++exampleFun' :: State (FunGraph Double (DIM0 :* DIM1 :* DIM2) (DIM2 :* DIM1 :* DIM0)) ()+exampleFun' = do+ let x = sym "x"+ y = vsym 5 "y"+ z = msym (3,5) "Z"+ z1 = (scale x z)**3+ z2 = (dot z y)**2+ z3 = diff ((x*x/2)**x) x++ inputs_ (x :* y :* z)+ outputs_ (z1 :* z2 :* z3)++run' :: IO ()+run' = do+ let gr :: FunGraph Double (DIM0 :* DIM1 :* DIM2) (DIM2 :* DIM1 :* DIM0)+ gr@(FunGraph hm im _ _) = runFunGraph exampleFun+ (FunGraph hm' im' _ _) = runFunGraph exampleFun'+ + putStrLn $ funGraphSummary gr+ putStrLn $ showCollisions gr+ previewGraph gr+ putStrLn "\nimperative same as pure+cse?:"+ print $ hm == hm'+ print $ im == im'++run :: IO ()+run = do+ let gr :: FunGraph Double (DIM0 :* DIM0) (DIM0 :* DIM0)+ gr@( FunGraph _ _ _ _) = runFunGraph $ do+ let x = sym "x"+ y = sym "y"+ z1 = x * y+ z2 = diff z1 x++ inputs_ (x :* y)+ outputs_ (z1 :* z2)++ putStrLn $ showCollisions gr+ putStrLn "-------------------------------------------"+ putStrLn $ funGraphSummary gr+ putStrLn "-------------------------------------------"+ putStrLn $ funGraphSummary' gr+ previewGraph gr++showoff :: IO ()+showoff = do+ let gr :: FunGraph Double (DIM0 :* DIM0 :* DIM0) (DIM0 :* DIM0 :* DIM0 :* DIM0)+ gr = makeFunGraph (x' :* y' :* z') (f :* fx :* fy :* fz)+ where+ x' = sym "x"+ y' = sym "y"+ z' = sym "z"++ f0 x y z = (z + x*y)*log(cos x / tanh y)**(z/exp y)+ fx0 = f0 (f0 x' y' z') (f0 z' y' x') (f0 y' x' z')+ fy0 = f0 (f0 z' x' y') (f0 x' z' y') (f0 z' z' y')+ fz0 = f0 (f0 x' y' z') (f0 x' y' x') (f0 y' x' y')+ f = f0 fx0 fy0 fz0+ + fx = diff f x'+ fy = diff f y'+ fz = diff f z'+++ putStrLn $ showCollisions gr+-- putStrLn $ funGraphSummary' gr+-- previewGraph gr
+ Dvda/Expr.hs view
@@ -0,0 +1,263 @@+{-# OPTIONS_GHC -Wall #-}+{-# Language TypeFamilies #-}+{-# Language MultiParamTypeClasses #-}+{-# Language GADTs #-}+{-# Language FlexibleInstances #-}+{-# Language FlexibleContexts #-}++module Dvda.Expr ( Expr(..)+ , FromGExpr+ , sym+ , vsym+ , msym+ , vec+ , mat+ , scale+ , dot+ , diff+ , grad+ , jacob+ , hess+ , dim+ , exprOfGExpr+ ) where++import Data.Array.Repa(DIM0,DIM1,DIM2,Z(..),(:.)(..), listOfShape, Shape(rank), shapeOfList)+import qualified Data.Vector.Unboxed as V+import Data.IntMap ( Key )++import Dvda.Dot ( Dot(..), dotDims )+import Dvda.BinUn ( BinOp(..), UnOp(..), showBinary, showUnary )+import Dvda.GExpr ( GExpr(..) )+import Dvda.HomoDim ( HomoDim, shapeOfHomo )++showShapeR :: Shape sh => sh -> String+showShapeR = show . reverse . listOfShape++class Shape sh => FromGExpr sh where+ fromMM :: HomoDim -> HomoDim -> Key -> Key -> Expr sh a+ fromMV :: HomoDim -> HomoDim -> Key -> Key -> Expr sh a+ fromVM :: HomoDim -> HomoDim -> Key -> Key -> Expr sh a+ fromVV :: HomoDim -> HomoDim -> Key -> Key -> Expr sh a+ fromMM shx shy = error $ "sorry, no fromMM instance for: " ++ show shx ++ ", " ++ show shy+ fromMV shx shy = error $ "sorry, no fromMV instance for: " ++ show shx ++ ", " ++ show shy+ fromVM shx shy = error $ "sorry, no fromVM instance for: " ++ show shx ++ ", " ++ show shy+ fromVV shx shy = error $ "sorry, no fromVV instance for: " ++ show shx ++ ", " ++ show shy++instance FromGExpr DIM2 where+ fromMM shx shy kx ky = EDot (ERef (shapeOfHomo shx :: DIM2) kx) (ERef (shapeOfHomo shy :: DIM2) ky)++instance FromGExpr DIM1 where+ fromMV shx shy kx ky = EDot (ERef (shapeOfHomo shx :: DIM2) kx) (ERef (shapeOfHomo shy :: DIM1) ky)+ fromVM shx shy kx ky = EDot (ERef (shapeOfHomo shx :: DIM1) kx) (ERef (shapeOfHomo shy :: DIM2) ky)++instance FromGExpr DIM0 where+ fromVV shx shy kx ky = EDot (ERef (shapeOfHomo shx :: DIM1) kx) (ERef (shapeOfHomo shy :: DIM1) ky)++dim :: Expr sh a -> sh+dim (ESym sh _) = sh+dim (EConst sh _) = sh+dim (EDimensionless _) = error "EDimensionless doesn't have a dimension, ya goon"+dim (ESingleton sh _) = sh+dim (EUnary _ x) = dim x+dim (EBinary _ x1 _) = dim x1+dim (EScale _ y) = dim y+dim (EDot x y) = dotDims (dim x) (dim y)+dim (ERef sh _) = sh+dim (EDeriv _ _) = Z+dim (EGrad _ args) = dim args+dim (EJacob x args) = Z :. head (listOfShape (dim x)) :. head (listOfShape (dim args))++exprOfGExpr :: (Shape sh, V.Unbox a, FromGExpr sh) => GExpr a -> Expr sh a+exprOfGExpr (GBinary sh' op kx ky) = EBinary op (ERef sh kx) (ERef sh ky)+ where+ sh = shapeOfHomo sh'+exprOfGExpr (GUnary sh op kx) = EUnary op (ERef (shapeOfHomo sh) kx)+exprOfGExpr (GSym sh name) = ESym (shapeOfHomo sh) name+exprOfGExpr (GSingleton sh a) = ESingleton (shapeOfHomo sh) a+exprOfGExpr (GScale sh kx ky) = EScale (ERef Z kx) (ERef (shapeOfHomo sh) ky)+exprOfGExpr (GConst sh v) = EConst (shapeOfHomo sh) v+exprOfGExpr (GDot shx shy kx ky) = case (rank shx, rank shy) of+ (2,2) -> fromMM shx shy kx ky+ (2,1) -> fromMV shx shy kx ky+ (1,2) -> fromVM shx shy kx ky+ (1,1) -> fromVV shx shy kx ky+ nm -> error $ "can't convert GDot of rank: " ++ show nm ++ " to Expr"++data Expr sh a where+ ESym :: sh -> String -> Expr sh a+ EConst :: V.Unbox a => sh -> V.Vector a -> Expr sh a+ EDimensionless :: a -> Expr sh a+ ESingleton :: sh -> a -> Expr sh a+ EUnary :: UnOp -> Expr sh a -> Expr sh a+ EBinary :: BinOp -> Expr sh a -> Expr sh a -> Expr sh a+ EScale :: Expr DIM0 a -> Expr sh a -> Expr sh a+ EDot :: Dot sh1 sh2 => Expr sh1 a -> Expr sh2 a -> Expr (DotT sh1 sh2) a+ ERef :: sh -> Int -> Expr sh a++ EDeriv :: Expr DIM0 a -> Expr DIM0 a -> Expr DIM0 a+ EGrad :: Expr DIM0 a -> Expr DIM1 a -> Expr DIM1 a+ EJacob :: Expr DIM1 a -> Expr DIM1 a -> Expr DIM2 a++isVal :: Eq a => a -> Expr sh a -> Bool+isVal x (EDimensionless y) = x == y+isVal x (ESingleton _ y) = x == y+isVal _ _ = False++-- | first layer of binary simplification: infer dimension of EDimensionless if possible+makeBinary :: (Num a, Eq a, Shape sh) => BinOp -> (a -> a -> a) -> Expr sh a -> Expr sh a -> Expr sh a+-- | can't infer dimension, just apply operation+makeBinary _ f (EDimensionless x) (EDimensionless y) = EDimensionless (f x y)+-- | infer dimension, then call makeBinary' for further simplification+makeBinary op f (EDimensionless x) y = makeBinary' op f (ESingleton (dim y) x) y+makeBinary op f x (EDimensionless y) = makeBinary' op f x (ESingleton (dim x) y)+-- | dimension inferred, call makeBinary'+makeBinary op f x y = makeBinary' op f x y++-- | second layer of binary simplification: check dimensions+makeBinary' :: (Num a, Eq a, Shape sh) => BinOp -> (a -> a -> a) -> Expr sh a -> Expr sh a -> Expr sh a+makeBinary' op f x y+ | shx == shy = makeBinary'' op f x y+ | otherwise = error $ "Binary op \""++ sop ++"\" dimension mismatch ya goon (" ++ sdx ++ ", " ++ sdy ++ ")"+ where+ shx = dim x+ shy = dim y+ sdx = showShapeR shx+ sdy = showShapeR shy+ sop = show op++-- | third layer of binary simplification: 0*x == x*0 == 0+-- | 1*x == x*1 == x+-- | 0+x == x+0 == x+-- | x/0 == error+-- | x/1 == x+-- | 0/x == 0+-- | x - 0 == 0+-- | 0 - x == neg x+makeBinary'' :: (Num a, Eq a, Shape sh) => BinOp -> (a -> a -> a) -> Expr sh a -> Expr sh a -> Expr sh a+makeBinary'' Mul f x y+ | isVal 0 x = x+ | isVal 0 y = y+ | isVal 1 x = y+ | isVal 1 y = x+ | otherwise = makeBinary''' Mul f x y+makeBinary'' Add f x y+ | isVal 0 x = y+ | isVal 0 y = x+ | otherwise = makeBinary''' Add f x y+makeBinary'' Div f x y+ | isVal 0 y = error "divide by zero"+ | isVal 1 y = x+ | isVal 0 x = x+ | otherwise = makeBinary''' Div f x y+makeBinary'' Sub f x y+ | isVal 0 x = negate y+ | isVal 0 y = x+ | otherwise = makeBinary''' Sub f x y++-- | apply operation to constant vectors+makeBinary'' _ f (EConst sh x) (EConst _ y) = EConst sh (V.zipWith f x y)+-- | broadcast constant operations+makeBinary'' _ f (ESingleton _ x) (EConst sh y) = EConst sh (V.map (f x) y)+makeBinary'' _ f (EConst sh x) (ESingleton _ y) = EConst sh (V.map (`f` y) x)+-- | otherwise make symbolic binary+makeBinary'' op _ x y = EBinary op x y++-- | fourth layer of binary simplification: make reasonable simplifications+makeBinary''' :: Shape sh => BinOp -> (a -> a -> a) -> Expr sh a -> Expr sh a -> Expr sh a+-- | apply operation to constant vectors+makeBinary''' _ f (EConst sh x) (EConst _ y) = EConst sh (V.zipWith f x y)+-- | broadcast constant operations+makeBinary''' _ f (ESingleton _ x) (EConst sh y) = EConst sh (V.map (f x) y)+makeBinary''' _ f (EConst sh x) (ESingleton _ y) = EConst sh (V.map (`f` y) x)+-- | otherwise make symbolic binary+makeBinary''' op _ x y = EBinary op x y+++-- | apply unary operations on constants+makeUnary :: Shape sh => UnOp -> (a -> a) -> Expr sh a -> Expr sh a+makeUnary _ f (EDimensionless x) = EDimensionless (f x)+makeUnary _ f (ESingleton sh x) = ESingleton sh (f x)+makeUnary _ f (EConst sh x) = EConst sh (V.map f x)+makeUnary op _ x = EUnary op x++instance (Shape sh, Num a, Eq a) => Num (Expr sh a) where+ (*) = makeBinary Mul (*)+ (+) = makeBinary Add (+)+ (-) = makeBinary Sub (-)+ abs = makeUnary Abs abs+ signum = makeUnary Signum signum+ fromInteger = EDimensionless . fromInteger+ negate = makeUnary Neg negate++instance (Shape sh, Fractional a, Eq a) => Fractional (Expr sh a) where+ (/) = makeBinary Div (/)+ fromRational = EDimensionless . fromRational++instance (Shape sh, Floating a, Eq a) => Floating (Expr sh a) where+ pi = EDimensionless pi+ (**) = makeBinary Pow (**)+ exp = makeUnary Exp exp+ log = makeUnary Log log+ sin = makeUnary Sin sin+ cos = makeUnary Cos cos+ asin = makeUnary ASin asin+ atan = makeUnary ATan atan+ acos = makeUnary ACos acos+ sinh = makeUnary Sinh sinh+ cosh = makeUnary Cosh cosh+ asinh = error "no instance for asinh"+ atanh = error "no instance for atanh"+ acosh = error "no instance for acosh"++paren :: Show a => a -> String+paren x = "( "++show x++" )"++instance (Shape sh, Show a) => Show (Expr sh a) where+ show (ESingleton _ x) = show x+ show (EDimensionless x) = show x+ show (ESym sh name) = name++"{"++showShapeR sh++"}"+ show (EConst sh x) = "{" ++ showShapeR sh ++ ", "++show (V.toList x)++"}" + show (EUnary op x) = showUnary x op+ show (EBinary op x y) = paren x ++ showBinary op ++ paren y+ show (EScale s x) = paren s ++ "*" ++ paren x+ show (EDot _ _) = "EDot ?? ??"+ show (ERef sh k) = "{ref:" ++ showShapeR sh ++ ":" ++ show k ++ "}"+ show (EDeriv x y) = "deriv(" ++ show x ++ ", " ++ show y ++ ")"+ show (EGrad x y) = "grad(" ++ show x ++ ", " ++ show y ++ ")"+ show (EJacob x y) = "jacob(" ++ show x ++ ", " ++ show y ++ ")"++sym :: String -> Expr DIM0 a+sym = ESym Z++vsym :: Int -> String -> Expr DIM1 a+vsym k = ESym (Z :. k)++msym :: (Int,Int) -> String -> Expr DIM2 a+msym (r,c) = ESym (Z :. r :. c)++vec :: V.Unbox a => [a] -> Expr DIM1 a+vec xs = EConst (shapeOfList [length xs]) (V.fromList xs)++mat :: V.Unbox a => (Int,Int) -> [a] -> Expr DIM2 a+mat (r,c) xs + | r*c == length xs = EConst (shapeOfList [c,r]) (V.fromList xs)+ | otherwise = error "bad dims in mat"++scale :: Expr DIM0 a -> Expr sh a -> Expr sh a+scale = EScale++dot :: (Dot sh1 sh2, DotT sh1 sh2 ~ sh) => Expr sh1 a -> Expr sh2 a -> Expr sh a+dot = EDot++diff :: Expr DIM0 a -> Expr DIM0 a -> Expr DIM0 a+diff = EDeriv++grad :: Expr DIM0 a -> Expr DIM1 a -> Expr DIM1 a+grad = EGrad++jacob :: Expr DIM1 a -> Expr DIM1 a -> Expr DIM2 a+jacob = EJacob++hess :: Expr DIM0 a -> Expr DIM1 a -> Expr DIM2 a+hess expr args = jacob (grad expr args) args
+ Dvda/GExpr.hs view
@@ -0,0 +1,87 @@+{-# OPTIONS_GHC -Wall #-}++module Dvda.GExpr ( GExpr(..)+ , getChildren+ , gdim+ ) where++import Data.IntMap ( Key )+import Data.Hashable ( Hashable, hash, combine )+import Data.GraphViz ( Labellable, toLabelValue )+import qualified Data.Vector.Unboxed as V++import Dvda.BinUn ( BinOp, UnOp, isCommutative )+import Dvda.HomoDim ( HomoDim(..) )+import Dvda.Dot ( dotDims )++simplifyCommutativeOps :: Bool+simplifyCommutativeOps = True++data GExpr a = GBinary HomoDim BinOp Key Key+ | GUnary HomoDim UnOp Key+ | GSym HomoDim String+ | GSingleton HomoDim a+ | GScale HomoDim Key Key+ | GDot HomoDim HomoDim Key Key+ | GConst HomoDim (V.Vector a) deriving Show++instance (Eq a, V.Unbox a) => Eq (GExpr a) where + (==) (GBinary shx opx x0 x1) (GBinary shy opy y0 y1) = and [opx == opy, shx == shy, commutativeConditions]+ where+ commutativeConditions = if simplifyCommutativeOps && isCommutative opx+ then (and [x0 == y0, x1 == y1]) || (and [x0 == y1, x1 == y0])+ else (and [x0 == y0, x1 == y1])+ (==) (GUnary shx opx x) (GUnary shy opy y) = and [shx == shy, opx == opy, x == y]+ (==) (GSym shx namex) (GSym shy namey) = and [shx == shy, namex == namey]+ (==) (GSingleton shx x) (GSingleton shy y) = and [shx == shy, x == y]+ (==) (GScale shx x0 x1) (GScale shy y0 y1) = and [shx == shy, x0 == y0, x1 == y1]+ (==) (GDot shx0 shx1 x0 x1) (GDot shy0 shy1 y0 y1) = and [shx0 == shy0, shx1==shy1, x0 == y0, x1 == y1]+ (==) (GConst shx x) (GConst shy y) = and [shx == shy, x == y]+ (==) _ _ = False++gdim :: GExpr a -> HomoDim+gdim (GBinary sh _ _ _) = sh+gdim (GUnary sh _ _) = sh+gdim (GSym sh _) = sh+gdim (GSingleton sh _) = sh+gdim (GScale sh _ _) = sh+gdim (GDot shx shy _ _) = dotDims shx shy+gdim (GConst sh _) = sh++instance (V.Unbox a, Hashable a) => Hashable (GExpr a) where+ hash (GBinary _ op k1 k2) = 24 `combine` hash op `combine` hk1 `combine` hk2+ where+ -- if the binary operator is commutative then always put the lesser hash first+ -- so that e.g. x*y and y*x are not computed twice+ (hk1, hk2)+ | simplifyCommutativeOps && isCommutative op && hk2' < hk1' = (hk2', hk1')+ | otherwise = (hk1', hk2')+ hk1' = hash k1+ hk2' = hash k2+ hash (GUnary _ op k) = 25 `combine` hash op `combine` hash k+ hash (GSym sh name) = 26 `combine` hash sh `combine` hash name+ hash (GSingleton sh x) = 27 `combine` hash sh `combine` hash x+ hash (GScale _ k1 k2) = 28 `combine` hash k1 `combine` hash k2+ hash (GDot _ _ k1 k2) = 29 `combine` hash k1 `combine` hash k2+ hash (GConst sh v) = V.foldl (\acc x -> acc `combine` hash x) (30 `combine` hash sh) v+++instance Show a => Labellable (GExpr a) where+ toLabelValue (GBinary _ op _ _) = toLabelValue $ show op+ toLabelValue (GUnary _ op _) = toLabelValue $ show op+ toLabelValue (GSym (HomoDim []) name) = toLabelValue name+ toLabelValue (GSym (HomoDim sh) name) = toLabelValue $ name ++ "{" ++ (tail . init . show . reverse) sh ++ "}"+ toLabelValue (GSingleton _ x) = toLabelValue $ show x+ toLabelValue (GScale {}) = toLabelValue "scale"+ toLabelValue (GDot {}) = toLabelValue "dot"+ toLabelValue (GConst {}) = toLabelValue "const"+++getChildren :: GExpr a -> [Int]+getChildren (GBinary _ _ k1 k2) = [k1,k2]+getChildren (GUnary _ _ k) = [k]+getChildren (GSym _ _ ) = []+getChildren (GSingleton _ _) = []+getChildren (GScale _ k1 k2) = [k1,k2]+getChildren (GDot _ _ k1 k2) = [k1,k2]+getChildren (GConst _ _) = []
+ Dvda/Graph.hs view
@@ -0,0 +1,114 @@+{-# OPTIONS_GHC -Wall #-}++module Dvda.Graph ( FunGraph(..)+ , FgNode+ , SymSet+ , emptyFunGraph+ , fgLookup+ , fgReverseLookup+ , fgGExprFromKey+ , previewGraph+ , toFGLGraph+ , collisions+ , showCollisions+ , funGraphSummary+ , funGraphSummary'+ ) where++import Data.Graph.Inductive ( Gr, mkGraph )+import Data.GraphViz ( preview )+import Control.Concurrent ( threadDelay )+import qualified Data.Vector.Unboxed as V( Unbox )+import Data.Hashable ( Hashable, hash )+import Data.List ( sort )+import Data.Maybe ( fromJust )+import Data.IntMap ( Key )+import qualified Data.HashSet as HS+import qualified Data.HashMap.Strict as HM+import qualified Data.IntMap as IM++import Dvda.GExpr ( GExpr(..), getChildren )++type SymSet a = HS.HashSet (GExpr a)+type FgNode a = (Key, SymSet a)++data FunGraph a b c = FunGraph+ (HM.HashMap (GExpr a) (FgNode a)) -- main lookup+ (IM.IntMap (GExpr a)) -- internal for reverse lookup+ (b,[Key])+ (c,[Key]) -- deriving Show+ +instance (Hashable a, V.Unbox a) => Hashable (FunGraph a b c) where+ hash (FunGraph _ im (_, inskeys) (_, outskeys)) = hash (IM.toList im, inskeys, outskeys)+ +fgLookup :: (Eq a, Hashable a, V.Unbox a) => GExpr a -> FunGraph a b c -> Maybe (FgNode a)+fgLookup gexpr (FunGraph hm _ _ _) = HM.lookup gexpr hm++fgReverseLookup :: (Eq a, Hashable a, V.Unbox a) => Key -> FunGraph a b c -> Maybe (FgNode a)+fgReverseLookup k fg = do+ gexpr <- fgGExprFromKey k fg+ fgLookup gexpr fg++fgGExprFromKey :: (Eq a, Hashable a, V.Unbox a) => Key -> FunGraph a b c -> Maybe (GExpr a)+fgGExprFromKey k (FunGraph _ im _ _) = IM.lookup k im++funGraphSummary :: (Show a, V.Unbox a, Show b, Show c) => FunGraph a b c -> String+funGraphSummary (FunGraph hm _ (b,bkeys) (c,ckeys)) =+ init $ unlines [ "input dims: " ++ show b+ , "input nodes:" ++ show bkeys+ , "output dims: " ++ show c+ , "output nodes:" ++ show ckeys+ , "number of nodes: " ++ show (HM.size hm)+ , "graph: " ++ show hm+ ]++-- more extensive+funGraphSummary' :: (Show a, V.Unbox a, Show b, Show c) => FunGraph a b c -> String+funGraphSummary' (FunGraph hm im (b,bkeys) (c,ckeys)) =+ init $ unlines [ "input dims: " ++ show b+ , "input nodes:" ++ show bkeys+ , "output dims: " ++ show c+ , "output nodes:" ++ show ckeys+ , "number of nodes: " ++ show (HM.size hm)+ , "graph:" + , init $ unlines (map show (IM.toList im))+ , "outputs:"+ , init $ unlines (map (show . (\k -> fromJust (IM.lookup k im))) ckeys)+ ]++collisions :: (Hashable a, V.Unbox a) => FunGraph a b c -> (Int, Int, Double)+collisions (FunGraph gr _ _ _) = (numCollisions, numTotal, fromIntegral numCollisions / fromIntegral numTotal)+ where+ allHashes = sort $ map (hash . fst) $ HM.toList gr+ numTotal = length allHashes+ numCollisions = countCollisions 0 allHashes+ where+ countCollisions n (x:y:ys)+ | x == y = countCollisions (n+1) (y:ys)+ | otherwise = countCollisions n (y:ys)+ countCollisions n [_] = n+ countCollisions n [] = n++showCollisions :: (Hashable a, V.Unbox a) => FunGraph a b c -> String+showCollisions gr = show numCollisions ++ '/' : show numTotal ++ " collisions ("++show (100*frac)++" %)"+ where+ (numCollisions, numTotal, frac) = collisions gr++emptyFunGraph :: FunGraph a b c+emptyFunGraph = FunGraph HM.empty IM.empty (inerr,inerr) (outerr,outerr)+ where+ inerr = error "must specify inputs"+ outerr = error "must specify outputs"+++previewGraph :: Show a => FunGraph a b c -> IO ()+previewGraph fungraph = do+ preview $ toFGLGraph fungraph+ threadDelay 10000++toFGLGraph :: FunGraph a b c -> Gr (GExpr a) String+toFGLGraph (FunGraph gexprs _ _ _) = mkGraph lnodes ledges+ where+ lnodes = map (\(x,(y,_)) -> (y,x)) $ HM.toList gexprs+-- lnodes = IM.toList gexprs+ ledges = concatMap (\(k,ge) -> map (\ch -> (ch,k,"")) (getChildren ge)) lnodes
+ Dvda/HSBuilder.hs view
@@ -0,0 +1,116 @@+{-# OPTIONS_GHC -Wall #-}+{-# Language GADTs #-}+{-# Language FlexibleContexts #-}+{-# Language TypeOperators #-}+{-# Language TypeFamilies #-}++module Dvda.HSBuilder ( buildHSFunction+ ) where++import qualified Data.Hashable as H+import qualified Data.Vector.Unboxed as V+import System.Directory+import Control.Monad(when)+import qualified System.Plugins.Make as Make+import qualified System.Plugins.Load as Load+--import System.Process( runCommand, waitForProcess )+--import System.Exit( ExitCode(ExitSuccess) )++import Dvda.HSSyntax ( writeHSSource )+import Dvda.Graph ( FunGraph(..) )+import Dvda.SymMonad ( Exprs )+import qualified Dvda.Config as Config+++-- | make source functions+buildHSFunction :: (H.Hashable a, Show a, V.Unbox a, Show b, Show c) =>+ FunGraph a b c -> IO (Exprs b Double -> Exprs c Double)+buildHSFunction fg = do+ -- source and hash+ let hash = show $ abs $ H.hash fg+ source = writeHSSource fg hash ++ -- function directory+ dir <- Config.functionDir hash+ + -- make function directory if it doesn't exist+ createDirectoryIfMissing False dir+ + -- filenames+ let sourcePath = dir ++ "/" ++ Config.nameHSSource hash+ -- objectPath = dir ++ "/" ++ Config.nameHSObject hash+ + -- if the source already exists, make sure it matches the old source+ srcExists <- doesFileExist sourcePath+ when srcExists $ do+ oldSrc <- readFile sourcePath+ when (source /= oldSrc) $ putStrLn $+ "====================================================\n" ++ + "WARNING: Hash not unique or source code has been edited\n"++ + "If you have not edited the auto-generated code, please let me\n" +++ "know that Hash collisions are a problem at gregmainland@gmail.com\n" +++ "====================================================\n\n"+ + -- write source+ putStrLn "writing source"+ writeFile sourcePath source++ -- compile code+ putStrLn "building source"+ objpath <- makeWithPlugins sourcePath+-- objpath <- makeWithProcess sourcePath objectPath+ + -- load object+ putStrLn "loading object"+ loadWithPlugins objpath hash+++makeWithPlugins :: FilePath -> IO FilePath+makeWithPlugins sourcePath = do + status <- Make.make sourcePath [] -- ["-v3"]+ + case status of (Make.MakeSuccess _ path) -> do putStrLn "Success!"+ return path+ (Make.MakeFailure code) -> do mapM_ putStrLn code+ error "Make Failure"+ ++loadWithPlugins :: FilePath -> String -> IO a+loadWithPlugins objpath hash = do+ status' <- Load.load_ objpath [] (Config.nameHSFunction hash)+ case status' of (Load.LoadFailure codes) -> do mapM_ putStrLn codes+ error "Load Failure"+ (Load.LoadSuccess _ fun) -> do putStrLn "load success!"+ return fun++-- -- | take in name of source and future object, compile object+-- makeWithProcess :: FilePath -> FilePath -> IO FilePath+-- makeWithProcess srcname objname = do+-- -- compile new object+-- let compileString = Config.ghcString srcname objname+-- displayString = Config.ghcString (shortName srcname) (shortName objname)+-- +-- -- print compilation string+-- when Config.spewGccCall $ putStrLn displayString+-- +-- -- run compilation string+-- p <- runCommand compileString+-- +-- -- check for errors+-- exitCode <- waitForProcess p+-- when (exitCode /= ExitSuccess) $ error $ "failed compiling " ++ srcname+-- +-- return objname+-- +-- +-- -- | shorten path name for display purposes+-- shortName :: String -> String+-- shortName full+-- | length name <= maxN = name ++ extension+-- | otherwise = take firstN name ++ "..." ++ drop (length name - lastN) name ++ extension+-- where+-- firstN = 20+-- lastN = 10+-- maxN = firstN + lastN+-- +-- (name, extension) = break (== '.') $ reverse $ takeWhile (/= '/') (reverse full)
+ Dvda/HSSyntax.hs view
@@ -0,0 +1,148 @@+{-# OPTIONS_GHC -Wall #-}++module Dvda.HSSyntax ( writeHSSource+ ) where++import Data.IntMap ( Key )+import Data.List ( intersperse )+import qualified Data.Vector.Unboxed as V+import qualified Data.IntMap as IM+import qualified Data.Text.Lazy as T++import Dvda.GExpr ( GExpr(..) )+import Dvda.Graph ( FunGraph(..) )+import Dvda.BinUn ( BinOp(..), UnOp(..) )+import qualified Dvda.Config as Config+++-- assign a scalar+sassign :: Key -> String+sassign k = Config.nameHSVar k ++ " = "++hBinary :: BinOp -> String+hBinary Add = "(+)"+hBinary Sub = "(-)"+hBinary Mul = "(*)"+hBinary Div = "(/)"+hBinary Pow = "(**)"+hBinary LogBase = "logBase"++hUnary :: UnOp -> String+hUnary Abs = "abs"+hUnary Neg = "negate"+hUnary Signum = "signum"+hUnary Exp = "exp"+hUnary Sqrt = "sqrt"+hUnary Log = "log"+hUnary Sin = "sin"+hUnary Cos = "cos"+hUnary Tan = "tan"+hUnary ASin = "asin"+hUnary ACos = "acos"+hUnary ATan = "atan"+hUnary Sinh = "sinh"+hUnary Cosh = "cosh"+hUnary Tanh = "tanh"+hUnary ASinh = "asinh"+hUnary ATanh = "atanh"+hUnary ACosh = "acosh"++pretty :: (Show a, V.Unbox a) => (Int, GExpr a) -> String+pretty (_, (GBinary _ op kx ky)) = hBinary op ++ " " ++ Config.nameHSVar kx ++ " " ++ Config.nameHSVar ky+pretty (_, (GUnary _ op kx)) = hUnary op ++ " " ++ Config.nameHSVar kx+pretty (_, (GSingleton _ x)) = show x+pretty (_, (GScale _ kx ky)) = "scale " ++ Config.nameHSVar kx ++ " " ++ Config.nameHSVar ky+pretty (_, (GDot _ _ kx ky)) = "dot " ++ Config.nameHSVar kx ++ " " ++ Config.nameHSVar ky+--pretty (k, (GConst _ vec)) = Config.nameHSConst k+pretty (_, (GConst _ vec)) = show vec -- Config.nameHSConst k+pretty (_, (GSym _ _)) = error "GSym shouldn't be handled here"++writeAssignment :: (Show a, V.Unbox a) => (Key, GExpr a) -> String+writeAssignment (k, gexpr@(GSym _ _)) = "-- " ++ Config.nameHSVar k ++ ": " ++ show gexpr+writeAssignment (k, gexpr) = sassign k ++ pretty (k,gexpr) ++ " -- " ++ show gexpr++writeHSSource :: (V.Unbox a, Show a, Show b, Show c) => FunGraph a b c -> String -> String+writeHSSource (FunGraph _ im (insT,ins) (outsT,outs)) hash =+ init $ unlines $+ [ "-- {-# OPTIONS_GHC -Wall #-}"+ , "{-# Language GADTs #-}"+ , "{-# Language FlexibleContexts #-}"+ , "{-# Language TypeOperators #-}"+ , "{-# Language TypeFamilies #-}"+ , ""+ , "module " ++ Config.nameHSModule hash ++ " ( " ++ Config.nameHSFunction hash ++ " ) where"+ , ""+ , "import Data.Array.Repa"+ , "import Data.Vector.Unboxed as V"+ , "import Dvda"+ , ""+-- , "-- constants:"+-- , constants+-- , ""+-- , Config.nameHSFunction hash ++ " :: Floating a => " +-- , spaces ++ rewriteType (show insT) ++ " -> " +-- , spaces ++ rewriteType (show outsT)+ , Config.nameHSFunction hash ++ " :: " ++ rewriteType (show insT) ++ " ->"+ , spaces ++ rewriteType (show outsT)+ , Config.nameHSFunction hash ++ " ( " ++ inputs ++ " ) = " ++ outputs+ , " where"+ , init $ unlines $ map (" " ++) body + ]+ where+ spaces = replicate ((length (Config.nameHSFunction hash)) + 4) ' '+ inputs = concat $ intersperse " :* " (map Config.nameHSVar ins)+ outputs = concat $ intersperse " :* " (map Config.nameHSVar outs)+ body = map writeAssignment (IM.toList im)+++intercalate :: String -> [String] -> String+intercalate _ [] = []+intercalate _ [x] = x+intercalate int (x:xs) = (x++int) ++ intercalate int xs++rewriteType :: String -> String+rewriteType typeString = final+ where+ text = T.pack typeString+ -- "Z :* ((Z :. 5) :* ((Z :. 3) :. 5))"+ + cleaned = T.filter (\x -> not (elem x "() ")) text+ -- "Z:*Z:.5:*Z:.3:.5"+ + grouped :: [T.Text]+ grouped = T.splitOn (T.pack ":*") cleaned+ -- ["Z", "Z:.5", "Z:.3:.5"]+ + + grouped' :: [[T.Text]]+ grouped' = map (T.splitOn (T.pack ":.")) grouped+ -- [["Z"], ["Z","5"], ["Z","3","5"]]++ counted :: [Int]+ counted = map (\x -> length x - 1) grouped'+ -- [0, 1, 2]++ addExpr = map (\x -> "(Expr DIM" ++ show x ++ " Double)") counted+ -- ["(Expr DIM0 Double)", "(Expr DIM1 Double)", "(Expr DIM2 Double)"]+ + final = "( " ++ (intercalate " :* " addExpr) ++ " )"+ -- "( (Expr DIM0 Double) :* (Expr DIM1 Double) :* (Expr DIM2 Double) )"+++-- rewriteType :: String -> String+-- rewriteType typeString = final+-- where+-- text = T.pack typeString+-- -- "Z :* ((Z :. 5) :* ((Z :. 3) :. 5))"+-- +-- cleaned = T.filter (\x -> not (elem x "() ")) text+-- -- "Z:*Z:.5:*Z:.3:.5"+-- +-- grouped = T.splitOn (T.pack ":*") cleaned+-- -- ["Z", "Z:.5", "Z:.3:.5"]+-- +-- addExpr = map (\x -> T.append "(Expr (" (T.append x ") a)")) grouped+-- -- ["(Expr (Z) a)", "(Expr (Z:.5) a)", "(Expr (Z:.3:.5) a)"]+-- +-- final = "( " ++ T.unpack (T.intercalate " :* " addExpr) ++ " )"+-- -- "( (Expr (Z) a) :* (Expr (Z:.5) a) :* (Expr (Z:.3:.5) a) )"
+ Dvda/HomoDim.hs view
@@ -0,0 +1,36 @@+{-# OPTIONS_GHC -Wall #-}++module Dvda.HomoDim ( HomoDim(..)+ , homoOfShape+ , shapeOfHomo+ ) where++import Control.DeepSeq+import Data.Array.Repa ( Shape(..) )+import Data.Hashable ( Hashable, hash )++newtype HomoDim = HomoDim [Int] deriving (Eq, Show)++homoOfShape :: Shape sh => sh -> HomoDim+homoOfShape = shapeOfList . listOfShape++shapeOfHomo :: Shape sh => HomoDim -> sh+shapeOfHomo = shapeOfList . listOfShape++instance Hashable HomoDim where+ hash (HomoDim xs) = hash xs++instance Shape HomoDim where+ listOfShape (HomoDim xs) = xs+ shapeOfList = HomoDim+ deepSeq xs y = listOfShape xs `deepseq` y+ rank = length . listOfShape+ size = product . listOfShape+ zeroDim = shapeOfList []+ addDim x y = shapeOfList $ zipWith (+) (listOfShape x) (listOfShape y)+ unitDim = error "need to finish instancing Shape HomoDim"+ intersectDim = error "need to finish instancing Shape HomoDim"+ sizeIsValid = error "need to finish instancing Shape HomoDim"+ toIndex = error "need to finish instancing Shape HomoDim"+ fromIndex = error "need to finish instancing Shape HomoDim"+ inShapeRange = error "need to finish instancing Shape HomoDim"
+ Dvda/SymMonad.hs view
@@ -0,0 +1,324 @@+{-# OPTIONS_GHC -Wall #-}+{-# Language GADTs #-}+{-# Language FlexibleContexts #-}+{-# Language TypeOperators #-}+{-# Language TypeFamilies #-}+{-# Language MultiParamTypeClasses #-}+{-# Language FlexibleInstances #-}++module Dvda.SymMonad ( (:*)(..)+ , HList(..)+ , Exprs+ , node+ , node'+ , inputs+ , inputs_+ , outputs+ , outputs_+ , makeFunGraph+ , runFunGraph+ , rad+ , getSensitivities+ ) where++import Control.Monad ( foldM, zipWithM )+import Control.Monad.State ( MonadState, StateT, get, put, liftM, runState )+import Data.Functor.Identity ( Identity )+import Data.Array.Repa ( Shape, Z, (:.) )+import Data.Hashable ( Hashable )+import Data.Vector.Unboxed ( Unbox )+import Data.Maybe ( fromJust )+import qualified Data.HashMap.Strict as HM+import qualified Data.HashSet as HS+import qualified Data.IntMap as IM+import Data.IntMap ( Key )+import Debug.Trace ( trace )++import Dvda.Dual ( Dual(..), dualPerturbation )+import Dvda.BinUn ( BinOp(..), applyUnary, applyBinary )+import Dvda.Graph ( FunGraph(..), emptyFunGraph, fgReverseLookup, fgGExprFromKey )+import Dvda.GExpr ( GExpr(..), gdim )+import Dvda.Expr ( Expr(..), FromGExpr, dim, exprOfGExpr )+import Dvda.HomoDim ( homoOfShape )++-- | take all sub expressions of an Expr and turn them into nodes+-- return an Expr that is just a ref+node :: (Shape sh, Hashable a, Unbox a, Floating a, Eq a) => Expr sh a -> StateT (FunGraph a b c) Identity (Expr sh a)+node expr = liftM (ERef (dim expr)) (node' expr)+ +node' :: (Shape sh, Hashable a, Unbox a, Floating a, Eq a) => Expr sh a -> StateT (FunGraph a b c) Identity Key+node' (EDimensionless _) = error "don't put EDimensionless in graph, ya goon"+node' (ERef _ k) = return k+node' (ESym sh name) = insert $ GSym (homoOfShape sh) name+node' (EConst sh x) = insert $ GConst (homoOfShape sh) x+node' (ESingleton sh x) = insert $ GSingleton (homoOfShape sh) x+node' (EUnary op x) = do+ xk <- node' x+ insert $ GUnary (homoOfShape $ dim x) op xk+node' (EBinary op x y) = do+ xk <- node' x+ yk <- node' y+ insert $ GBinary (homoOfShape $ dim x) op xk yk+node' (EScale x y) = do+ xk <- node' x+ yk <- node' y+ insert $ GScale (homoOfShape $ dim y) xk yk+node' (EDot x y) = do+ xk <- node' x+ yk <- node' y+ let shx = homoOfShape $ dim x+ shy = homoOfShape $ dim y+ insert $ GDot shx shy xk yk+node' (EDeriv x' arg') = do+ x <- node x'+ arg <- node arg'+ outs <- rad x [arg]+ node' (head outs)+node' (EGrad x' arg') = do+ x <- node x'+ arg <- node arg'+ outs <- rad x [arg]+ node' (head outs)+++-- | Try to insert the GExpr into the hashmap performing CSE.+-- If the GExpr is not yet in the map, insert it and return new key.+-- Otherwise don't insert, just return existing key.+insert :: (Hashable a, Unbox a, Floating a, Eq a) => GExpr a -> StateT (FunGraph a b c) Identity Key+insert gexpr = do+ fg <- get+ let symSet (GSym _ _) = HS.singleton gexpr+ symSet (GSingleton _ _) = HS.empty+ symSet (GConst _ _) = HS.empty+ symSet (GUnary _ _ k) = snd $ fromJust $ fgReverseLookup k fg+ symSet (GBinary _ _ xk yk) = symMapX `HS.union` symMapY+ where+ (_,symMapX) = fromJust $ fgReverseLookup xk fg+ (_,symMapY) = fromJust $ fgReverseLookup yk fg+ symSet (GScale _ xk yk) = symMapX `HS.union` symMapY+ where+ (_,symMapX) = fromJust $ fgReverseLookup xk fg+ (_,symMapY) = fromJust $ fgReverseLookup yk fg+ symSet (GDot _ _ xk yk) = symMapX `HS.union` symMapY+ where+ (_,symMapX) = fromJust $ fgReverseLookup xk fg+ (_,symMapY) = fromJust $ fgReverseLookup yk fg++ (FunGraph hm im ins outs) <- get+ case HM.lookup gexpr hm of+ Just (k',_) -> return k'+ Nothing -> do let k = HM.size hm+ hm' = HM.insert gexpr (k,symSet gexpr) hm+ im' = IM.insert k gexpr im+ put (FunGraph hm' im' ins outs)+ return k+++gexprOfExpr :: (Eq a, Floating a, Hashable a, Unbox a, Shape sh, FromGExpr sh) =>+ Expr sh a -> StateT (FunGraph a b c) Identity (GExpr a)+gexprOfExpr expr = do+ k <- node' expr+ fg <- get+ return (fromJust $ fgGExprFromKey k fg)+ +-- gradient of expression w.r.t. list of args+rad :: (Eq a, Hashable a, Unbox a, Floating a, Shape sh, FromGExpr sh, Shape sh0, FromGExpr sh0) => + Expr sh0 a -> [Expr sh a] -> StateT (FunGraph a b c) Identity [Expr sh a]+rad expr_ args_ = do+ expr <- gexprOfExpr expr_+ args <- mapM gexprOfExpr args_+ let argSet = HS.fromList args+ sensitivities <- getSensitivities argSet expr (ESingleton (dim expr_) 1)+ -- order inputs requested by user+ let getSens x argDim = case HM.lookup x sensitivities of+ Just sens -> return sens+ Nothing -> trace "WARNING: taking deriviative df/dx where f is not a function of x (inserting 0 in graph)" $+ node' (ESingleton argDim 0)+ argDims = map dim args_+ orderedSensitivities <- zipWithM getSens args argDims+ return $ zipWith ERef argDims orderedSensitivities+++-- combine two (GExpr, Key) hashmaps+-- if there is a conflict, add the two GExprs together+unionWithPlus :: (Eq a, Floating a, Hashable a, Unbox a) =>+ HM.HashMap (GExpr a) Key -> HM.HashMap (GExpr a) Key ->+ StateT (FunGraph a b c) Identity (HM.HashMap (GExpr a) Key)+unionWithPlus xs ys = foldM addCommon union0 commonGExprs+ where+ -- the gexprs that occur in both maps+ commonGExprs = HM.keys $ HM.intersection xs ys+ -- the initial union that needs conflicts fixed+ union0 = xs `HM.union` ys+ addCommon hm commonGExpr = do+ let xsensk = fromJust $ HM.lookup commonGExpr xs+ ysensk = fromJust $ HM.lookup commonGExpr ys+ k <- insert $ GBinary (gdim commonGExpr) Add xsensk ysensk+ return (HM.insert commonGExpr k hm)+ ++lookupSymSet :: (Eq a, Hashable a, Unbox a) => Key -> StateT (FunGraph a b c) Identity (HS.HashSet (GExpr a))+lookupSymSet k = do+ fg <- get+ let (_,symSet) = fromJust $ fgReverseLookup k fg+ return symSet+++getSensitivities :: (Eq a, Hashable a, Unbox a, Floating a, Shape sh, FromGExpr sh) => + HS.HashSet (GExpr a) -> GExpr a -> Expr sh a ->+ StateT (FunGraph a b c) Identity (HM.HashMap (GExpr a) Key)+getSensitivities _ (GSingleton _ _) _ = return HM.empty+getSensitivities _ (GConst _ _) _ = return HM.empty+getSensitivities args primal@(GSym _ _) sens = if HS.member primal args then do+ k <- node' sens+ return $ HM.fromList [(primal, k)]+ -- don't backprop if there aren't any interesting symbols farther in the tree+ else return HM.empty+getSensitivities args (GUnary _ op gk) sens = do+ symSetG <- lookupSymSet gk+ case HS.size (HS.intersection args symSetG) of+ -- don't backprop if there aren't any interesting symbols farther in the tree+ 0 -> return HM.empty+ _ -> do+ fg <- get+ let g' = fromJust $ fgGExprFromKey gk fg+ g = exprOfGExpr g'+ dfdg = dualPerturbation $ applyUnary op (Dual g 1)+ getSensitivities args g' (sens*dfdg)+getSensitivities args (GBinary _ op gk hk) sens = do+ symSetG <- lookupSymSet gk+ symSetH <- lookupSymSet hk+ + fg <- get+ let g' = fromJust $ fgGExprFromKey gk fg+ h' = fromJust $ fgGExprFromKey hk fg+ g = exprOfGExpr g'+ h = exprOfGExpr h'+ dfdg = dualPerturbation $ applyBinary op (Dual g 1) (Dual h 0)+ dfdh = dualPerturbation $ applyBinary op (Dual g 0) (Dual h 1)+ + gsens <- case HS.size (HS.intersection args symSetG) of+ 0 -> return HM.empty+ _ -> getSensitivities args g' (sens*dfdg)+ hsens <- case HS.size (HS.intersection args symSetH) of+ 0 -> return HM.empty+ _ -> getSensitivities args h' (sens*dfdh)+ unionWithPlus gsens hsens+getSensitivities args (GDot _ _ gk hk) sens = do+ symSetG <- lookupSymSet gk+ symSetH <- lookupSymSet hk+ + fg <- get+ let g' = fromJust $ fgGExprFromKey gk fg+ h' = fromJust $ fgGExprFromKey hk fg+ g = exprOfGExpr g'+ h = exprOfGExpr h'+ dfdg = h+ dfdh = g+ + gsens <- case HS.size (HS.intersection args symSetG) of+ 0 -> return HM.empty+ _ -> getSensitivities args g' (sens*dfdg)+ hsens <- case HS.size (HS.intersection args symSetH) of+ 0 -> return HM.empty+ _ -> getSensitivities args h' (sens*dfdh)+ unionWithPlus gsens hsens+getSensitivities args (GScale _ gk hk) sens = do+ symSetG <- lookupSymSet gk+ symSetH <- lookupSymSet hk+ + fg <- get+ let g' = fromJust $ fgGExprFromKey gk fg+ h' = fromJust $ fgGExprFromKey hk fg+ g = exprOfGExpr g'+ h = exprOfGExpr h'+ dfdg = h+ dfdh = g+ + gsens <- case HS.size (HS.intersection args symSetG) of+ 0 -> return HM.empty+ _ -> getSensitivities args g' (sens*dfdg)+ hsens <- case HS.size (HS.intersection args symSetH) of+ 0 -> return HM.empty+ _ -> getSensitivities args h' (sens*dfdh)+ unionWithPlus gsens hsens+ ++++---------------------- heterogenous inputs/outputs ------------------+data a :* b = a :* b deriving Show+infixr 6 :*++class HList a where+ type NumT a+ type DimT a+-- mkNodes :: (NumT a ~ b) => a -> State (FunGraph b c d) (a,[Key])+ mkNodes :: a -> StateT (FunGraph (NumT a) b c) Identity (a,[Key])+ getHDim :: a -> DimT a++instance (HList a, HList b, NumT a ~ NumT b) => HList (a :* b) where+ type NumT (a :* b) = NumT a+ type DimT (a :* b) = DimT a :* DimT b+ mkNodes (x :* y) = do+ (exs,kxs) <- mkNodes x+ (eys,kys) <- mkNodes y+ return (exs :* eys, kxs++kys)+ getHDim (x :* y) = getHDim x :* getHDim y++instance (Shape sh, Hashable a, Unbox a, Eq a, Floating a) => HList (Expr sh a) where+ type NumT (Expr sh a) = a+ type DimT (Expr sh a) = sh+ mkNodes expr = do+ expr'@(ERef _ k) <- node expr+ return (expr', [k])+ getHDim = dim+ +inputs :: HList b => b -> StateT (FunGraph (NumT b) (DimT b) c) Identity b+inputs exprs = do+ (exprs', keys) <- mkNodes exprs+ FunGraph hm im _ outs <- get+ put (FunGraph hm im (getHDim exprs, keys) outs)+ return exprs'++outputs :: HList c => c -> StateT (FunGraph (NumT c) b (DimT c)) Identity c+outputs exprs = do+ (exprs',keys) <- mkNodes exprs+ FunGraph hm im ins _ <- get+ put (FunGraph hm im ins (getHDim exprs,keys))+ return exprs'++inputs_ :: HList b => b -> StateT (FunGraph (NumT b) (DimT b) c) Identity ()+inputs_ exprs = do+ _ <- inputs exprs+ return ()++outputs_ :: HList c => c -> StateT (FunGraph (NumT c) b (DimT c)) Identity ()+outputs_ exprs = do+ _ <- outputs exprs+ return ()++--------------------------------------------------------------+class ExprList sh a where+ type Exprs sh a+ +instance (ExprList sh0 a, ExprList sh1 a) => ExprList (sh0 :* sh1) a where+ type Exprs (sh0 :* sh1) a = (Exprs sh0 a) :* (Exprs sh1 a)+ +instance ExprList Z a where+ type Exprs Z a = Expr Z a++instance Shape sh => ExprList (sh :. Int) a where+ type Exprs (sh :. Int) a = Expr (sh :. Int) a+++---------------- utility function -----------------+runFunGraph :: StateT (FunGraph a b c) Identity d -> FunGraph a b c+runFunGraph f = snd $ runState f emptyFunGraph++--makeFunGraph :: (HList c, HList b, NumT b ~ NumT c, NumT b ~ a, Eq a, Floating a, Hashable a, Unbox a) =>+makeFunGraph :: (HList c, HList b, NumT b ~ NumT c, NumT b ~ a) =>+ b -> c -> FunGraph a (DimT b) (DimT c)+makeFunGraph ins outs = runFunGraph $ do+ inputs_ ins+ outputs_ outs
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2011, Greg Horn++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ dvda.cabal view
@@ -0,0 +1,102 @@+Name: dvda+Version: 0.1+License: BSD3+License-file: LICENSE+Author: Greg Horn+Maintainer: gregmainland@gmail.edu+Stability: Experimental+Category: Math+Build-type: Simple+Synopsis: Efficient automatic differentiation+Cabal-version: >= 1.6+Description: {+DVDA Verifiably Differentiates Algorithmically+.+This library provides a symbolic type `Dvda.Expr` which is+manipulated mathematically through its Num\/Fractional\/Floating instances.+Expr can be a scalar, vector, or matrix. Binary operations (adding\/multiplying\/etc)+are all elementwise.+.+Matrix/vector/scalar compatability is enforced at compile time+.++Efficient derivatives can be computed. Internally reverse automatic differentiation+is performed including efficient common subexpression elimination.+.+Function graphs can be JIT compiled into efficienty functions using "buildHSFunction".+This is the intended way to use this library.+.+Pretty graphviz plots!+.+If the runtime JIT stuff works in terminal ghci but not emacs haskell-mode, you may need to add+`(setenv "PATH" (concatenate 'string (getenv "PATH") ":/usr/local/bin"))` to your .emacs file+.+To get started look in `Dvda.Examples` or CompileTest.hs in the github repo+}++source-repository head+ type: git+ location: git://github.com/ghorn/dvda.git+-- tag: ++Flag stressTest+ Description: Build a profilable hard executable+ Default: False++Library+ Exposed-modules: Dvda+ Dvda.BinUn+ Dvda.Config+ Dvda.Dot+ Dvda.Dual+ Dvda.Examples+ Dvda.Expr+ Dvda.GExpr+ Dvda.Graph+ Dvda.HSBuilder+ Dvda.HSSyntax+ Dvda.HomoDim+ Dvda.SymMonad+-- Dvda.CFunction+-- Dvda.Codegen.CBuilder+-- Dvda.Codegen.CCallWrapper+-- Dvda.Codegen.CSyntax+-- Dvda.Codegen.Utils++ Other-modules: ++ Build-depends: base >= 4 && < 5,+ hashable >= 1.1 && < 1.2,+ vector >= 0.9 && < 0.10,+ repa >= 3.1 && < 3.2,+ containers >= 0.4 && < 0.5,+ unordered-containers >= 0.2 && < 0.3,+ graphviz >= 2999.12 && < 2999.13,+ fgl >= 5.4 && < 5.5,+ mtl >= 2.0 && < 2.1,+ directory >= 1.1 && < 1.2,+ process >= 1.1 && < 1.2,+ text >= 0.11 && < 0.12,+ transformers >= 0.2 && < 0.3,+ plugins >= 1.5 && < 1.6,+ deepseq >= 1.3 && < 1.4+-- unix+-- text,+-- QuickCheck,++ Ghc-options: -Wall+-- Ghc-options: -O2 -Wall -threaded+ GHC-Prof-Options: -prof -fprof-auto+++-- Executable stressTest+-- if flag(stressTest)+-- Buildable: True+-- else+-- Buildable: False+-- +-- Main-Is: StressTest.hs+-- +-- Ghc-Options: -O2+-- +-- GHC-Prof-Options: -prof -fprof-auto