diff --git a/Dvda.hs b/Dvda.hs new file mode 100644 --- /dev/null +++ b/Dvda.hs @@ -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 + + diff --git a/Dvda/BinUn.hs b/Dvda/BinUn.hs new file mode 100644 --- /dev/null +++ b/Dvda/BinUn.hs @@ -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++" )" diff --git a/Dvda/Config.hs b/Dvda/Config.hs new file mode 100644 --- /dev/null +++ b/Dvda/Config.hs @@ -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 diff --git a/Dvda/Dot.hs b/Dvda/Dot.hs new file mode 100644 --- /dev/null +++ b/Dvda/Dot.hs @@ -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 + diff --git a/Dvda/Dual.hs b/Dvda/Dual.hs new file mode 100644 --- /dev/null +++ b/Dvda/Dual.hs @@ -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) diff --git a/Dvda/Examples.hs b/Dvda/Examples.hs new file mode 100644 --- /dev/null +++ b/Dvda/Examples.hs @@ -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 diff --git a/Dvda/Expr.hs b/Dvda/Expr.hs new file mode 100644 --- /dev/null +++ b/Dvda/Expr.hs @@ -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 diff --git a/Dvda/GExpr.hs b/Dvda/GExpr.hs new file mode 100644 --- /dev/null +++ b/Dvda/GExpr.hs @@ -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 _ _) = [] diff --git a/Dvda/Graph.hs b/Dvda/Graph.hs new file mode 100644 --- /dev/null +++ b/Dvda/Graph.hs @@ -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 diff --git a/Dvda/HSBuilder.hs b/Dvda/HSBuilder.hs new file mode 100644 --- /dev/null +++ b/Dvda/HSBuilder.hs @@ -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) diff --git a/Dvda/HSSyntax.hs b/Dvda/HSSyntax.hs new file mode 100644 --- /dev/null +++ b/Dvda/HSSyntax.hs @@ -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) )" diff --git a/Dvda/HomoDim.hs b/Dvda/HomoDim.hs new file mode 100644 --- /dev/null +++ b/Dvda/HomoDim.hs @@ -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" diff --git a/Dvda/SymMonad.hs b/Dvda/SymMonad.hs new file mode 100644 --- /dev/null +++ b/Dvda/SymMonad.hs @@ -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 diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/dvda.cabal b/dvda.cabal new file mode 100644 --- /dev/null +++ b/dvda.cabal @@ -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