diff --git a/Dvda.hs b/Dvda.hs --- a/Dvda.hs +++ b/Dvda.hs @@ -6,23 +6,34 @@ -} {-# OPTIONS_GHC -Wall #-} +{-# Language TypeOperators #-} +{-# Language TypeFamilies #-} +{-# Language MultiParamTypeClasses #-} +{-# Language FlexibleInstances #-} module Dvda ( -- * primitives sym , vsym , msym + , svec + , smat , vec , mat -- * operations , scale - , dot +-- , dot , diff + , runDeriv -- * symbolic expression type , Expr + , fullShow + , fullShowNodes -- * construct FunGraphs , FunGraph , makeFunGraph , runFunGraph + , inputs + , outputs , inputs_ , outputs_ , node @@ -32,15 +43,41 @@ , showCollisions , previewGraph -- * compile and link function - , buildHSFunction +-- , buildHSFunction +-- , buildHSFunctionPure +-- , buildHSFunctionFromGraph -- * Heterogenous inputs/outputs , (:*)(..) , Exprs + -- * re-export from repa and hmatrix +-- , Index + , DIM0 + , DIM1 + , DIM2 + , Z(..) + , (:.) + , Matrix + , Vector ) where import Dvda.Expr import Dvda.Graph -import Dvda.HSBuilder +--import Dvda.HSBuilder import Dvda.SymMonad +import Data.Array.Repa ( DIM0, DIM1, DIM2, Z(..), Shape, (:.) ) +import Numeric.LinearAlgebra ( Matrix, Vector ) +-- | Just a nice way to write (Exprs (DIM0 :* DIM1 :* DIM2) Double) +-- | instead of (Expr DIM0 Double :* Expr DIM1 Double :* Expr DIM2 Double) +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 diff --git a/Dvda/BinUn.hs b/Dvda/BinUn.hs --- a/Dvda/BinUn.hs +++ b/Dvda/BinUn.hs @@ -9,6 +9,8 @@ , unaryDeriv , binaryDeriv , isCommutative + , lassoc + , rassoc ) where import Data.Hashable ( Hashable, hash ) @@ -32,14 +34,14 @@ | Cosh | ATanh | ASinh - | ACosh deriving (Eq, Show) + | ACosh deriving (Eq, Show, Enum, Bounded) data BinOp = Add | Sub | Mul | Div | Pow - | LogBase deriving (Eq, Show) + | LogBase deriving (Eq, Show, Enum, Bounded) instance Hashable UnOp where hash Abs = 0 @@ -69,8 +71,8 @@ hash Pow = 22 hash LogBase = 23 -showUnary :: Show a => a -> UnOp -> String -showUnary x Abs = '|': show x ++ "|" +showUnary :: String -> UnOp -> String +showUnary x Abs = '|': x ++ "|" showUnary x Neg = '-':paren x showUnary x Signum = "signum"++paren x showUnary x Exp = "exp"++paren x @@ -139,5 +141,51 @@ isCommutative Pow = False isCommutative LogBase = False -paren :: Show a => a -> String -paren x = "( "++show x++" )" +lassoc :: BinOp -> BinOp -> Bool +lassoc Add Add = True -- a + b + c == (a + b) + c +lassoc Add Sub = True -- a + b - c == (a + b) - c +--lassoc Add Mul = True -- a + b * c == (a + b) * c +--lassoc Add Div = True -- a + b / c == (a + b) / c + +lassoc Sub Add = True -- a - b + c == (a - b) + c +lassoc Sub Sub = True -- a - b - c == (a - b) - c +--lassoc Sub Mul = True -- a - b * c == (a - b) * c +--lassoc Sub Div = True -- a - b / c == (a - b) / c + +lassoc Div Add = True -- a / b + c == (a / b) + c +lassoc Div Sub = True -- a / b - c == (a / b) - c +lassoc Div Mul = True -- a / b * c == (a / b) * c +lassoc Div Div = True -- a / b / c == (a / b) / c + +lassoc Mul Add = True -- a * b + c == (a * b) + c +lassoc Mul Sub = True -- a * b - c == (a * b) - c +lassoc Mul Mul = True -- a * b * c == (a * b) * c +lassoc Mul Div = True -- a * b / c == (a * b) / c + +lassoc _ _ = False + +rassoc :: BinOp -> BinOp -> Bool +--rassoc Add Add = True -- a + b + c == a + (b + c) +--rassoc Add Sub = True -- a + b - c == a + (b - c) +rassoc Add Mul = True -- a + b * c == a + (b * c) +rassoc Add Div = True -- a + b / c == a + (b / c) + +--rassoc Sub Add = True -- a - b + c == a - (b + c) +--rassoc Sub Sub = True -- a - b - c == a - (b - c) +rassoc Sub Mul = True -- a - b * c == a - (b * c) +rassoc Sub Div = True -- a - b / c == a - (b / c) + +--rassoc Div Add = True -- a / b + c == a / (b + c) +--rassoc Div Sub = True -- a / b - c == a / (b - c) +--rassoc Div Mul = True -- a / b * c == a / (b * c) +--rassoc Div Div = True -- a / b / c == a / (b / c) + +--rassoc Mul Add = True -- a * b + c == a * (b + c) +--rassoc Mul Sub = True -- a * b - c == a * (b - c) +--rassoc Mul Mul = True -- a * b * c == a * (b * c) +--rassoc Mul Div = True -- a * b / c == a * (b / c) + +rassoc _ _ = False + +paren :: String -> String +paren x = "( "++ x ++" )" diff --git a/Dvda/CallNative.hs b/Dvda/CallNative.hs new file mode 100644 --- /dev/null +++ b/Dvda/CallNative.hs @@ -0,0 +1,202 @@ +{-# OPTIONS_GHC -Wall #-} +{-# Language TypeOperators #-} +{-# Language TypeFamilies #-} +{-# Language FlexibleContexts #-} +{-# Language FlexibleInstances #-} +{-# Language GADTs #-} + +module Dvda.CallNative ( toNative + , nativeCall + , nativeDiff + , nativeGrad + , nativeJacob + , nativeRun + ) where + +import Data.Hashable ( Hashable ) +import Data.HashMap.Lazy ( HashMap ) +import qualified Data.HashMap.Lazy as HM +import qualified Data.IntMap as IM +import Data.List ( mapAccumL ) +import Data.Maybe ( fromJust, catMaybes ) +import Numeric.LinearAlgebra ( Element, Container ) + +import Dvda +import Dvda.BinUn ( BinOp(Mul), applyBinary, applyUnary ) +import Dvda.Expr ( Expr(..), Const(..), dim ) +import Dvda.Graph ( FunGraph(..), DvdaDim(..), DynamicExpr, fgLookup, fgExprFromKey ) +import Dvda.SymMonad ( rad ) + +class (Hashable (INumT b), Eq (INumT b), Element (INumT b)) => NativeInputs b where + type INumT b + toReplacements :: FunGraph (INumT b) b c -> b -> HashMap (DynamicExpr (INumT b)) (DynamicExpr (INumT b)) + +insToSyms :: DvdaDim sh => FunGraph a b c -> Expr sh a -> Expr sh a -> Maybe (DynamicExpr a, DynamicExpr a) +insToSyms fg e@(ERef _ k) out = fmap (\x -> (makeDynamic x, makeDynamic out)) $ fgExprFromKey (dim e) k fg +insToSyms _ _ _ = Nothing + +instance (DvdaDim sh, Hashable a, Element a, Eq a) => NativeInputs (Expr sh a) where + type INumT (Expr sh a) = a + toReplacements fg@(FunGraph _ _ ins _) xs = HM.fromList $ catMaybes [insToSyms fg ins xs] + +instance (DvdaDim sh, Hashable a, Element a, Eq a) => NativeInputs [Expr sh a] where + type INumT [Expr sh a] = a + toReplacements fg@(FunGraph _ _ ins _) xs = HM.fromList $ catMaybes $ zipWith (insToSyms fg) ins xs + +instance (DvdaDim sh, Hashable a, Element a, Eq a) => NativeInputs [[Expr sh a]] where + type INumT [[Expr sh a]] = a + toReplacements fg@(FunGraph _ _ ins _) xs = + HM.fromList $ catMaybes $ zipWith (insToSyms fg) (concat ins) (concat xs) + +instance (NativeInputs a, NativeInputs b, INumT a ~ INumT b) => NativeInputs (a :* b) where + type INumT (a :* b) = INumT a + toReplacements (FunGraph hm im (in0 :* in1) outs) (x0 :* x1) = HM.union r0 r1 + where + r0 = toReplacements (FunGraph hm im in0 outs) x0 + r1 = toReplacements (FunGraph hm im in1 outs) x1 + +--------------------------------------------------------------------------- +class NativeOutput c where + type ONumT c + traverseOutputs :: (NativeInputs b) + => HashMap (DynamicExpr (ONumT c)) (DynamicExpr (ONumT c)) + -> FunGraph (ONumT c) b c + -> c + -> (FunGraph (ONumT c) b c, c) + +instance (DvdaDim sh, Floating a, Num (Vector a), Container Vector a, Hashable a, Eq a) + => NativeOutput (Expr sh a) where + type ONumT (Expr sh a) = a + traverseOutputs = eval + +instance (DvdaDim sh, Floating a, Num (Vector a), Container Vector a, Hashable a, Eq a) + => NativeOutput [Expr sh a] where + type ONumT [Expr sh a] = a + traverseOutputs = mapAccumL . eval + +instance (DvdaDim sh, Floating a, Num (Vector a), Container Vector a, Hashable a, Eq a) + => NativeOutput [[Expr sh a]] where + type ONumT [[Expr sh a]] = a + traverseOutputs = mapAccumL . mapAccumL . eval + +instance (NativeOutput a, NativeOutput b, ONumT a ~ ONumT b) => NativeOutput (a :* b) where + type ONumT (a :* b) = ONumT a + traverseOutputs replacementMap (FunGraph hm0 im0 ins outs) (x' :* y') = (FunGraph hm2 im2 ins outs, x :* y) + where + err = error "DON'T LOOK AT THESE OUTPUTS YA GOON" + (FunGraph hm1 im1 _ _, x) = traverseOutputs replacementMap (FunGraph hm0 im0 ins err) x' + (FunGraph hm2 im2 _ _, y) = traverseOutputs replacementMap (FunGraph hm1 im1 ins err) y' + + +replace :: (Hashable a, Eq a, Element a, DvdaDim sh) => FunGraph a b c -> Expr sh a -> Expr sh a -> FunGraph a b c +replace fg0@(FunGraph hm0 im0 ins outs) old new = FunGraph hm im ins outs + where + (k, _) = fromJust $ fgLookup old fg0 + hm = HM.insert (makeDynamic new) (k, error "after callNative has happened you can't look at symSets") hm0 + im = IM.insert k (makeDynamic new) im0 + + +eval :: (Hashable a, Eq a, Floating a, Num (Vector a), Container Vector a, DvdaDim sh) + => HashMap (DynamicExpr a) (DynamicExpr a) -> FunGraph a b c -> Expr sh a -> (FunGraph a b c, Expr sh a) +eval _ _ (EDimensionless _) = error "WHO PUT AN EDimensionless IN THIS GRAPH" +eval _ _ (EDeriv _ _) = error "WHO PUT AN EDeriv IN THIS GRAPH" +eval _ _ (EGrad _ _) = error "WHO PUT AN EDeriv IN THIS GRAPH" +eval _ _ (EJacob _ _) = error "WHO PUT AN EJacob IN THIS GRAPH" +eval replacementMap fg expr@(ERef _ k) = eval replacementMap fg (fromJust $ fgExprFromKey (dim expr) k fg) +eval _ fg expr@(EConst _) = (fg, expr) +eval replacementMap fg0 expr@(ESym _ _) = case HM.lookup (makeDynamic expr) replacementMap of + Nothing -> (fg0, expr) + Just replacementExpr' -> (fg1, replacementExpr) + where + replacementExpr = fromDynamic (dim expr) replacementExpr' + fg1 = replace fg0 expr replacementExpr +eval replacementMap fg0 expr@(EUnary op x') = (fg2, newExpr) + where + (fg1, x) = eval replacementMap fg0 x' + newExpr = applyUnary op x + fg2 = replace fg1 expr newExpr +eval replacementMap fg0 expr@(EBinary op x' y') = (fg3, newExpr) + where + (fg1, x) = eval replacementMap fg0 x' + (fg2, y) = eval replacementMap fg1 y' + newExpr = applyBinary op x y + fg3 = replace fg2 expr newExpr +eval replacementMap fg (EScale (EConst (CSingleton _ x)) y) = eval replacementMap fg z + where + z = applyBinary Mul (EConst (CSingleton (dim y) x)) y +eval replacementMap fg0 expr@(EScale x' y') = (fg3, newExpr) + where + (fg1, x) = eval replacementMap fg0 x' + (fg2, y) = eval replacementMap fg1 y' + newExpr = case x of EConst (CSingleton _ c) -> applyBinary Mul (EConst (CSingleton (dim y) c)) y + _ -> EScale x y + fg3 = replace fg2 expr newExpr + +toNative :: (Show a, NativeInputs b, NativeOutput c, a ~ INumT b, a ~ ONumT c) => FunGraph a b c -> b -> c +toNative fg@(FunGraph _ _ _ outs) xs = snd $ traverseOutputs replacementMap fg outs + where + replacementMap = toReplacements fg xs + + +-- | Convenience function for natively computing function This is +-- expected to be very slow. Using code generation instead is +-- recommended +nativeCall :: (Hashable a, Eq a, Show a, Element a, Floating a, Num (Vector a), Container Vector a) + => (Expr Z a -> [Expr Z a]) -> Expr Z a -> [Expr Z a] +nativeCall f = toNative $ runFunGraph $ do + let x = sym "x" + inputs_ x + outputs_ (f x) + +-- | Lift a unary function over @Floating a => a@ to a function over +-- @Floating a => Expr Z a@ +liftNative :: (Hashable a, Eq a, Show a, Element a, Floating a, + Num (Vector a), Container Vector a, + Floating b, b ~ Expr Z a) => (b -> b) -> Expr Z a -> Expr Z a +liftNative f x = case nativeCall (return . f) x of + [] -> error "Function didn't return." + (v:_) -> v + +-- | Evaluate a unary function over @Floating a => a@ using Dvda's +-- internal machinery. The typeclass constraints should make sure the +-- error doesn't happen, but it could anyway. +nativeRun :: (Hashable a, Eq a, Show a, Element a, Floating a, + Num (Vector a), Container Vector a, + Floating b, b ~ Expr Z a) => (b -> b) -> a -> a +nativeRun f x = case liftNative f (EConst (CSingleton Z x)) of + (EConst (CSingleton Z v)) -> v + _ -> error "Function must be unary over class Floating." + + +-- | Convenience function for natively computing jacobian, requires +-- you to pass the number of inputs. This is expected to be very +-- slow. Using code generation instead is recommended +nativeJacob :: (Hashable a, Eq a, Show a, Element a, Floating a, Num (Vector a), Container Vector a) + => Int -> ([Expr Z a] -> [Expr Z a]) -> [Expr Z a] -> [[Expr Z a]] +nativeJacob n f = toNative $ runFunGraph $ do + let xs = map (\k -> sym ("x_"++show k)) [0..(n-1::Int)] + inputs_ xs + ys <- mapM (flip rad xs) (f xs) + outputs_ ys + +-- | Convenience function for natively computing gradient, requires +-- you to pass the number of inputs. This is expected to be very +-- slow. Using code generation instead is recommended +nativeGrad :: (Hashable a, Eq a, Show a, Element a, Floating a, Num (Vector a), Container Vector a) + => Int -> ([Expr Z a] -> Expr Z a) -> [Expr Z a] -> [Expr Z a] +nativeGrad n f = toNative $ runFunGraph $ do + let xs = map (\k -> sym ("x_"++show k)) [0..(n-1::Int)] + inputs_ xs + ys <- rad (f xs) xs + outputs_ ys + +-- | Convenience function for natively computing a derivative. This +-- is expected to be very slow. Using code generation instead is +-- recommended +nativeDiff :: (Hashable a, Eq a, Show a, Element a, Floating a, Num (Vector a), Container Vector a) + => (Expr Z a -> Expr Z a) -> Expr Z a -> Expr Z a +nativeDiff f = toNative $ runFunGraph $ do + let x = sym "x" + inputs_ x + [y] <- rad (f x) [x] + outputs_ y diff --git a/Dvda/Codegen.hs b/Dvda/Codegen.hs new file mode 100644 --- /dev/null +++ b/Dvda/Codegen.hs @@ -0,0 +1,38 @@ +{-# OPTIONS_GHC -Wall #-} + +module Dvda.Codegen ( writeSourceFile + ) where + +import System.Directory +import Control.Monad(when) + +import qualified Dvda.Config as Config + +writeSourceFile :: String -> FilePath -> FilePath -> IO FilePath +writeSourceFile source functionDir sourceName = do + topDir <- Config.dvdaDir + let dir = topDir ++ "/" ++ functionDir + + -- make function directory if it doesn't exist + createDirectoryIfMissing False dir + + -- filenames + let sourcePath = dir ++ "/" ++ sourceName + -- 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 " ++ sourcePath + writeFile sourcePath source + + return sourcePath diff --git a/Dvda/Config.hs b/Dvda/Config.hs --- a/Dvda/Config.hs +++ b/Dvda/Config.hs @@ -4,7 +4,6 @@ module Dvda.Config( -- * directory stuff dvdaDir - , functionDir -- * C syntax , cType , cName @@ -19,12 +18,17 @@ , nameHSSource , nameHSVar , nameHSConst + -- * Octave + , nameOctaveSource + , nameOctaveFunction -- * gcc stuff , gccString , spewGccCall , outputNames -- * ghc stuff , ghcString + -- * symbolic stuff + , simplifyCommutativeOps ) where import System.Directory @@ -62,12 +66,6 @@ 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 @@ -112,5 +110,16 @@ nameHSSource :: String -> String nameHSSource = (++ ".hs") . nameHSModule +nameOctaveSource :: String -> String +nameOctaveSource = (++ ".m") + +nameOctaveFunction :: String -> String +nameOctaveFunction hash = hash + nameHSObject :: String -> String nameHSObject = (++ ".o") . nameHSModule + +-- | whether e.g. x + y == y + x or not +simplifyCommutativeOps :: Bool +simplifyCommutativeOps = True + diff --git a/Dvda/Dot.hs b/Dvda/Dot.hs deleted file mode 100644 --- a/Dvda/Dot.hs +++ /dev/null @@ -1,67 +0,0 @@ -{-# 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 --- a/Dvda/Dual.hs +++ b/Dvda/Dual.hs @@ -5,6 +5,8 @@ {-# Language TypeFamilies #-} module Dvda.Dual ( Dual(..) + , fad + , fad' ) where import Data.Ratio ( numerator, denominator ) @@ -58,3 +60,11 @@ 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) + +-- | Forward derivative propogation. fad' [sin x, 2*x] == [cos x, 2] +fad' :: Num a => (Dual a -> [Dual a]) -> a -> [a] +fad' f x = map dualPerturbation $ f (Dual x 1) + +-- | Forward derivative propogation. fad sin x == cos x +fad :: Num a => (Dual a -> Dual a) -> a -> a +fad f x = dualPerturbation $ f (Dual x 1) diff --git a/Dvda/Examples.hs b/Dvda/Examples.hs --- a/Dvda/Examples.hs +++ b/Dvda/Examples.hs @@ -1,50 +1,69 @@ {-# OPTIONS_GHC -Wall #-} {-# Language TypeOperators #-} -module Dvda.Examples ( exampleFun - , run +module Dvda.Examples ( run , run' , showoff + , bigGraph + , smallGraph + , runCallNative + , composed ) where -import Control.Monad.State (State) -import Data.Array.Repa (DIM0,DIM1,DIM2) +import Data.Array.Repa.Index +import Control.Monad.State import Dvda +import Dvda.Expr +import Dvda.CallNative import Dvda.Graph ( FunGraph(..) ) -exampleFun :: State (FunGraph Double (DIM0 :* DIM1 :* DIM2) (DIM2 :* DIM1 :* DIM0)) () -exampleFun = do - let x = sym "x" +exampleFunGraph :: State (FunGraph + Double (Exprs (DIM0 :* DIM1 :* DIM2) Double) + (Exprs (DIM2 :* DIM1 :* DIM0) Double)) + () +exampleFunGraph = do + let x = sym "x" :: Expr DIM0 Double 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 +-- z2 <- node $ (dot z y)**2 + z2 <- node $ 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" +pureFun :: Exprs (DIM0 :* DIM1 :* DIM2) Double -> Exprs (DIM2 :* DIM1 :* DIM0) Double +pureFun (x :* y :* z) = z1 :* z2 :* z3 + where + z1 = (scale x z)**3 +-- z2 = (dot z y)**2 + z2 = y**2 + z3 = diff ((x*x/2)**x) x + +exampleFunGraph' :: State (FunGraph + Double + (Exprs (DIM0 :* DIM1 :* DIM2) Double) + (Exprs (DIM2 :* DIM1 :* DIM0) Double)) + () +exampleFunGraph' = do + let x = sym "x" :: Expr DIM0 Double 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) + + args = x :* y :* z + + inputs_ args + outputs_ (pureFun args) 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' + let gr@(FunGraph hm im _ _) = runFunGraph exampleFunGraph + (FunGraph hm' im' _ _) = runFunGraph exampleFunGraph' - putStrLn $ funGraphSummary gr + putStrLn $ funGraphSummary' gr putStrLn $ showCollisions gr previewGraph gr putStrLn "\nimperative same as pure+cse?:" @@ -53,43 +72,80 @@ run :: IO () run = do - let gr :: FunGraph Double (DIM0 :* DIM0) (DIM0 :* DIM0) - gr@( FunGraph _ _ _ _) = runFunGraph $ do - let x = sym "x" + let gr@( FunGraph _ _ _ _) = runFunGraph $ do + let x = sym "x" :: Expr DIM0 Double y = sym "y" - z1 = x * y + z1 = x + x / y + 3 z2 = diff z1 x + z3 = diff z1 y inputs_ (x :* y) - outputs_ (z1 :* z2) + outputs_ (z1 :* z2 :* z3) putStrLn $ showCollisions gr - putStrLn "-------------------------------------------" - putStrLn $ funGraphSummary gr - putStrLn "-------------------------------------------" - putStrLn $ funGraphSummary' gr + putStrLn $ fullShowNodes gr + let FunGraph _ _ _ (z:* zx :* zy) = gr + putStrLn $ "\nz: " ++ fullShow gr z + putStrLn $ "dz/dx: " ++ fullShow gr zx + putStrLn $ "dz/dy: " ++ fullShow gr zy 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" +bigGraph :: FunGraph Double + (Exprs (DIM0 :* DIM0 :* DIM0) Double) + (Exprs (DIM0 :* DIM0 :* DIM0 :* DIM0) Double) +bigGraph = makeFunGraph (x' :* y' :* z') (f :* fx :* fy :* fz) + where + x' = sym "x" :: Expr DIM0 Double + 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' - 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' +smallGraph :: FunGraph Double + (Exprs (DIM0 :* DIM0 :* DIM0) Double) + (Exprs (DIM0 :* DIM0) Double) +smallGraph = makeFunGraph (x :* y :* z) (f0 :* f1) + where + x = sym "x" :: Expr DIM0 Double + y = sym "y" + z = sym "z" + f0 = x*y*z + 3 + f1 = 40*f0/x - putStrLn $ showCollisions gr --- putStrLn $ funGraphSummary' gr --- previewGraph gr +runCallNative :: Exprs (Z :* Z) Double +runCallNative = toNative smallGraph (f 1 :* f 2 :* f 3) + where + f = EConst . (CSingleton Z) + +showoff :: IO () +showoff = do + putStrLn $ showCollisions bigGraph + let FunGraph _ _ _ (f :* fx :* fy :* fz) = bigGraph + putStrLn "--------------------------------------------------------------" + putStrLn $ fullShow bigGraph f + putStrLn "--------------------------------------------------------------" + putStrLn $ fullShow bigGraph fx + putStrLn "--------------------------------------------------------------" + putStrLn $ fullShow bigGraph fy + putStrLn "--------------------------------------------------------------" + putStrLn $ fullShow bigGraph fz + putStrLn "--------------------------------------------------------------" +-- putStrLn $ funGraphSummary' bigGraph + previewGraph bigGraph + +composed :: [Expr Z Double] +composed = runDeriv z [t] + where + t = sym "t" + x = symDependent "x" t + y = symDependent "y" x + z = symDependent "z" y diff --git a/Dvda/Expr.hs b/Dvda/Expr.hs --- a/Dvda/Expr.hs +++ b/Dvda/Expr.hs @@ -1,121 +1,214 @@ -{-# OPTIONS_GHC -Wall #-} -{-# Language TypeFamilies #-} -{-# Language MultiParamTypeClasses #-} +{-# Options_ghc -Wall #-} +{-# Language StandaloneDeriving #-} +{-# Language DeriveDataTypeable #-} {-# Language GADTs #-} -{-# Language FlexibleInstances #-} {-# Language FlexibleContexts #-} module Dvda.Expr ( Expr(..) - , FromGExpr + , Const(..) + , Sym(..) , sym + , svec + , smat , vsym , msym , vec , mat , scale - , dot +-- , dot , diff , grad , jacob , hess , dim - , exprOfGExpr + , isVal + , symDependent ) where -import Data.Array.Repa(DIM0,DIM1,DIM2,Z(..),(:.)(..), listOfShape, Shape(rank), shapeOfList) -import qualified Data.Vector.Unboxed as V +import Data.Array.Repa(DIM0,DIM1,DIM2,Z(..),(:.)(..), listOfShape, Shape(shapeOfList), rank ) +import Numeric.LinearAlgebra ( Matrix, Vector, Element ) +import qualified Numeric.LinearAlgebra as LA +import Foreign.Storable ( Storable ) import Data.IntMap ( Key ) +import Data.Hashable ( Hashable, hash, combine ) +import Data.List ( sort ) +import Data.Typeable ( Typeable2 ) -import Dvda.Dot ( Dot(..), dotDims ) -import Dvda.BinUn ( BinOp(..), UnOp(..), showBinary, showUnary ) -import Dvda.GExpr ( GExpr(..) ) -import Dvda.HomoDim ( HomoDim, shapeOfHomo ) +import Dvda.BinUn ( BinOp(..), UnOp(..), showBinary, showUnary, isCommutative, lassoc, rassoc ) +import Dvda.Config ( simplifyCommutativeOps ) +import Dvda.SparseLA ( SparseVec, SparseMat, svFromList, smFromLists ) 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 (EConst (CSingleton sh _)) = sh +dim (EConst (CMat sh _)) = sh +dim (EConst (CVec sh _)) = sh +dim (EConst (CTensor 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" +deriving instance Typeable2 Const +deriving instance Typeable2 Expr +data Const sh a where + CSingleton :: sh -> a -> Const sh a + CVec :: DIM1 -> Vector a -> Const DIM1 a + CMat :: DIM2 -> Matrix a -> Const DIM2 a + CTensor :: sh -> Vector a -> Const sh a + +data Sym = Sym String -- doesn't depend on independent variable, or is an independent variable + | SymDependent String Int Sym -- depends on independent variable, Int specifies the nth derivative + deriving Eq + +instance Show Sym where + show (Sym name) = name + show (SymDependent name k s) = name ++ replicate k '\'' ++ "(" ++ show s ++ ")" + data Expr sh a where - ESym :: sh -> String -> Expr sh a - EConst :: V.Unbox a => sh -> V.Vector a -> Expr sh a + ESym :: sh -> Sym -> Expr sh a + EConst :: Const sh 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 + ERef :: sh -> Key -> Expr sh a EDeriv :: Expr DIM0 a -> Expr DIM0 a -> Expr DIM0 a - EGrad :: Expr DIM0 a -> Expr DIM1 a -> Expr DIM1 a + EGrad :: Expr DIM0 a -> Expr sh a -> Expr sh a EJacob :: Expr DIM1 a -> Expr DIM1 a -> Expr DIM2 a +--------------------------------- show instances ----------------------------- +instance (Shape sh, Show a, Element a) => Show (Const sh a) where + show (CSingleton _ x) = show x + show (CVec sh v) = "CVec " ++ showShapeR sh ++ " " ++ show v + show (CMat sh m) = "CMat " ++ showShapeR sh ++ " " ++ show m + show (CTensor sh v) = "CTensor " ++ showShapeR sh ++ " " ++ show v + +paren :: String -> String +paren x = "("++ x ++")" + +instance (Shape sh, Show a, Element a) => Show (Expr sh a) where + show (ERef sh k) + | rank sh == 0 = "{ref:" ++ show k ++ "}" + | otherwise = "{ref:" ++ show k ++ ",(" ++ showShapeR sh ++ ")}" + show (EDimensionless x) = show x + show (ESym sh s) + | rank sh == 0 = show s + | otherwise = show s++"{"++showShapeR sh++"}" + show (EConst x) = show x + show (EUnary op x) = showUnary (show x) op + show (EBinary op x y) = parenx x (show x) ++ " " ++ showBinary op ++ " " ++ pareny y (show y) + where + parenx (EBinary xop _ _) = if lassoc xop op then id else paren + parenx (EScale _ _) = if lassoc Mul op then id else paren + parenx _ = id + + pareny (EBinary yop _ _) = if rassoc op yop then id else paren + pareny (EScale _ _) = if rassoc op Mul then id else paren + pareny _ = id + show (EScale x y) = parenx x (show x) ++ " " ++ showBinary Mul ++ " " ++ pareny y (show y) + where + parenx (EBinary xop _ _) = if lassoc xop Mul then id else paren + parenx (EScale _ _) = if lassoc Mul Mul then id else paren + parenx _ = id + + pareny (EBinary yop _ _) = if rassoc Mul yop then id else paren + pareny (EScale _ _) = if rassoc Mul Mul then id else paren + pareny _ = id + + 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 ++ ")" + + +--------------------------------- eq instances ------------------------- +instance (Shape sh, Element a, Eq a) => Eq (Const sh a) where + (==) (CSingleton sh0 x0) (CSingleton sh1 x1) = sh0 == sh1 && x0 == x1 + (==) (CVec sh0 v0) (CVec sh1 v1) = sh0 == sh1 && v0 == v1 + (==) (CMat sh0 m0) (CMat sh1 m1) = sh0 == sh1 && (LA.flatten m0) == (LA.flatten m1) + (==) (CTensor sh0 v0) (CTensor sh1 v1) = sh0 == sh1 && v0 == v1 + (==) _ _ = False + +instance (Shape sh, Eq a, Element a) => Eq (Expr sh a) where + (==) (ESym sh0 name0) (ESym sh1 name1) = sh0 == sh1 && name0 == name1 + (==) (EConst c0) (EConst c1) = c0 == c1 + (==) (EDimensionless x0) (EDimensionless x1) = x0 == x1 + (==) (EUnary op0 x0) (EUnary op1 x1) = op0 == op1 && x0 == x1 + (==) (EScale x0 y0) (EScale x1 y1) = x0 == x1 && y0 == y1 + (==) (ERef sh0 k0) (ERef sh1 k1) = sh0 == sh1 && k0 == k1 + (==) (EDeriv x0 y0) (EDeriv x1 y1) = x0 == x1 && y0 == y1 + (==) (EGrad x0 y0) (EGrad x1 y1) = x0 == x1 && y0 == y1 + (==) (EJacob x0 y0) (EJacob x1 y1) = x0 == x1 && y0 == y1 + (==) (EBinary op0 x0 y0) (EBinary op1 x1 y1) = op0 == op1 && commutativeEq + where + commutativeEq + | simplifyCommutativeOps && isCommutative op0 = (x0 == x1 && y0 == y1) || (x0 == y1 && y0 == x1) + | otherwise = x0 == x1 && y0 == y1 + (==) _ _ = False + +------------------------- hashable instances -------------------- +instance (Hashable a, Shape sh, Element a) => Hashable (Const sh a) where + hash (CSingleton sh x) = 24 `combine` hash (listOfShape sh) `combine` hash x + hash (CVec sh v) = LA.foldVector (\x acc -> acc `combine` hash x) (25 `combine` hash (listOfShape sh)) v + hash (CMat sh v) = LA.foldVector (\x acc -> acc `combine` hash x) (26 `combine` hash (listOfShape sh)) (LA.flatten v) + hash (CTensor sh v) = LA.foldVector (\x acc -> acc `combine` hash x) (27 `combine` hash (listOfShape sh)) v + + +instance (Hashable a, Shape sh, Element a) => Hashable (Expr sh a) where + hash (ESym sh name) = 28 `combine` hash (listOfShape sh) `combine` hash name + hash (EConst c) = 29 `combine` hash c + hash (EDimensionless x) = 30 `combine` hash x +-- hash (EBroadcast sh x) = 30 `combine` hash (listOfShape sh) `combine` hash x + hash (EUnary op x) = 31 `combine` hash op `combine` hash x + hash (EBinary op x y) = 32 `combine` hash op `combine` hashx `combine` hashy + where + [hashx,hashy] + | simplifyCommutativeOps && isCommutative op = sort unsorted + | otherwise = unsorted + where + unsorted = [hash x, hash y] + hash (EScale x y) = 33 `combine` hash x `combine` hash y + hash (ERef sh k) = 34 `combine` hash (listOfShape sh) `combine` k + + hash (EDeriv x y) = 35 `combine` hash x `combine` hash y + hash (EGrad x y) = 36 `combine` hash x `combine` hash y + hash (EJacob x y) = 37 `combine` hash x `combine` hash y + +instance Hashable Sym where + hash (Sym name) = 38 `combine` hash name + hash (SymDependent name k s) = 39 `combine` hash name `combine` k `combine` hash s + + +------------------------ symbolic stuff -------------------- isVal :: Eq a => a -> Expr sh a -> Bool isVal x (EDimensionless y) = x == y -isVal x (ESingleton _ y) = x == y +isVal x (EConst (CSingleton _ 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 +makeBinary :: (Eq a, Num (Vector a), LA.Container Vector 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) +makeBinary op f (EDimensionless x) y = makeBinary' op f (EConst (CSingleton (dim y) x)) y +makeBinary op f x (EDimensionless y) = makeBinary' op f x (EConst (CSingleton (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' :: (Eq a, Num (Vector a), LA.Container Vector 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 ++ ")" @@ -126,15 +219,17 @@ 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 +-- 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'' :: (Eq a, Num (Vector a), LA.Container Vector 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 @@ -154,34 +249,59 @@ | isVal 0 x = negate y | isVal 0 y = x | otherwise = makeBinary''' Sub f x y +makeBinary'' op f x y = makeBinary''' op 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) +makeBinary''' :: (Num (Vector a), LA.Container Vector a) => + BinOp -> (a -> a -> a) -> Expr sh a -> Expr sh a -> Expr sh a +-- apply vectorized operations +makeBinary''' Add _ (EConst (CVec sh x)) (EConst (CVec _ y)) = EConst $ CVec sh (x + y) +makeBinary''' Sub _ (EConst (CVec sh x)) (EConst (CVec _ y)) = EConst $ CVec sh (x - y) +makeBinary''' Mul _ (EConst (CVec sh x)) (EConst (CVec _ y)) = EConst $ CVec sh (x * y) +makeBinary''' Div _ (EConst (CVec sh x)) (EConst (CVec _ y)) = EConst $ CVec sh (x / y) +makeBinary''' Add _ (EConst (CMat sh x)) (EConst (CMat _ y)) = EConst $ CMat sh (x + y) +makeBinary''' Sub _ (EConst (CMat sh x)) (EConst (CMat _ y)) = EConst $ CMat sh (x - y) +makeBinary''' Mul _ (EConst (CMat sh x)) (EConst (CMat _ y)) = EConst $ CMat sh (x * y) +makeBinary''' Div _ (EConst (CMat sh x)) (EConst (CMat _ y)) = EConst $ CMat sh (x / y) +makeBinary''' Add _ (EConst (CTensor sh x)) (EConst (CTensor _ y)) = EConst $ CTensor sh (x + y) +makeBinary''' Sub _ (EConst (CTensor sh x)) (EConst (CTensor _ y)) = EConst $ CTensor sh (x - y) +makeBinary''' Mul _ (EConst (CTensor sh x)) (EConst (CTensor _ y)) = EConst $ CTensor sh (x * y) +makeBinary''' Div _ (EConst (CTensor sh x)) (EConst (CTensor _ y)) = EConst $ CTensor sh (x / y) +makeBinary''' _ f (EConst x') (EConst y') = EConst $ czipWith x' y' + where + -- zip like things + czipWith (CSingleton sh x) (CSingleton _ y) = CSingleton sh (f x y) + czipWith (CTensor sh x) (CTensor _ y) = CTensor sh (LA.zipVectorWith f x y) + czipWith (CVec sh x) (CVec _ y) = CVec sh (LA.zipVectorWith f x y) + czipWith (CMat sh x) (CMat _ y) = CMat sh (LA.reshape (LA.cols x) z) + where + z = LA.zipVectorWith f (LA.flatten x) (LA.flatten y) + -- broadcast singletons + czipWith (CSingleton _ x) (CTensor sh y) = CTensor sh (LA.mapVector (f x) y) + czipWith (CSingleton _ x) (CVec sh y) = CVec sh (LA.mapVector (f x) y) + czipWith (CSingleton _ x) (CMat sh y) = CMat sh (LA.mapMatrix (f x) y) + czipWith (CTensor sh x) (CSingleton _ y) = CTensor sh (LA.mapVector (`f` y) x) + czipWith (CVec sh x) (CSingleton _ y) = CVec sh (LA.mapVector (`f` y) x) + czipWith (CMat sh x) (CSingleton _ y) = CMat sh (LA.mapMatrix (`f` y) x) + czipWith _ _ = error "czipWith called on unlike constants" -- | 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 :: Storable a => 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 _ f' (EConst x') = EConst $ cmap f' x' + where + cmap f (CSingleton sh x) = CSingleton sh (f x) + cmap f (CTensor sh x) = CTensor sh (LA.mapVector f x) + cmap f (CVec sh x) = CVec sh (LA.mapVector f x) + cmap f (CMat sh x) = CMat sh (LA.mapMatrix f x) makeUnary op _ x = EUnary op x -instance (Shape sh, Num a, Eq a) => Num (Expr sh a) where +instance (Shape sh, Num a, Eq a, Num (Vector a), LA.Container Vector a) => + Num (Expr sh a) where (*) = makeBinary Mul (*) (+) = makeBinary Add (+) (-) = makeBinary Sub (-) @@ -190,11 +310,13 @@ fromInteger = EDimensionless . fromInteger negate = makeUnary Neg negate -instance (Shape sh, Fractional a, Eq a) => Fractional (Expr sh a) where +instance (Shape sh, Fractional a, Eq a, Num (Vector a), LA.Container Vector a) => + Fractional (Expr sh a) where (/) = makeBinary Div (/) fromRational = EDimensionless . fromRational -instance (Shape sh, Floating a, Eq a) => Floating (Expr sh a) where +instance (Shape sh, Floating a, Eq a, Num (Vector a), LA.Container Vector a) => + Floating (Expr sh a) where pi = EDimensionless pi (**) = makeBinary Pow (**) exp = makeUnary Exp exp @@ -210,45 +332,51 @@ 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 ++ ")" - +------------------------------ convenience functions ------------------------- +-- | symbolic scalar sym :: String -> Expr DIM0 a -sym = ESym Z +sym = (ESym Z) . Sym +symDependent :: String -> Expr DIM0 a -> Expr DIM0 a +symDependent name (ESym _ s) = ESym Z (SymDependent name 0 s) +symDependent _ _ = error "symDependent got non ESym dependency" + +-- | symbolic dense vector vsym :: Int -> String -> Expr DIM1 a -vsym k = ESym (Z :. k) +vsym k = (ESym (Z :. k)) . Sym +-- | symbolic dense matrix msym :: (Int,Int) -> String -> Expr DIM2 a -msym (r,c) = ESym (Z :. r :. c) +msym (r,c) = (ESym (Z :. r :. c)) . Sym -vec :: V.Unbox a => [a] -> Expr DIM1 a -vec xs = EConst (shapeOfList [length xs]) (V.fromList xs) +-- | symbolic dense constant vector +vec :: Storable a => [a] -> Expr DIM1 a +vec xs = EConst $ CVec (shapeOfList [length xs]) (LA.fromList xs) -mat :: V.Unbox a => (Int,Int) -> [a] -> Expr DIM2 a +-- | symbolic dense constant matrix +mat :: Element 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" + | r*c == sum (map length xs) && r == length xs = EConst $ CMat (shapeOfList [c,r]) (LA.fromLists xs) + | otherwise = error $ "bad dims in mat!"++ + "\ngiven (r,c): " ++ show (r,c) ++ + "\nactual (r,c): " ++ show (length xs, map length xs) +-- | symbolic sparse vector +svec :: String -> Int -> SparseVec (Expr DIM0 a) +svec name len = svFromList $ map (\k -> sym $ name ++ "_" ++ show k) [0..len-1] + +-- | symbolic sparse matrix +smat :: String -> (Int,Int) -> SparseMat (Expr DIM0 a) +smat name (rows,cols) = smFromLists allRcs + where + allRcs = map (\row -> map (\col -> (sym $ name ++ "_" ++ show row ++ "_" ++ show col)) [0..cols-1]) [0..rows-1] + + 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 +--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 diff --git a/Dvda/GExpr.hs b/Dvda/GExpr.hs deleted file mode 100644 --- a/Dvda/GExpr.hs +++ /dev/null @@ -1,87 +0,0 @@ -{-# 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 --- a/Dvda/Graph.hs +++ b/Dvda/Graph.hs @@ -1,82 +1,156 @@ {-# OPTIONS_GHC -Wall #-} +{-# Language StandaloneDeriving #-} +{-# Language TypeSynonymInstances #-} +{-# Language FlexibleInstances #-} +{-# Language GADTs #-} +{-# Language RankNTypes #-} module Dvda.Graph ( FunGraph(..) + , DynamicExpr(..) + , DvdaDim(..) , FgNode , SymSet , emptyFunGraph , fgLookup - , fgReverseLookup - , fgGExprFromKey + , fgExprFromKey + , insert , previewGraph , toFGLGraph , collisions , showCollisions , funGraphSummary , funGraphSummary' + , showNodes + , asIfExpr ) where import Data.Graph.Inductive ( Gr, mkGraph ) -import Data.GraphViz ( preview ) +import Data.GraphViz ( Labellable, toLabelValue, preview ) +import Data.GraphViz.Attributes.Complete ( Label ) import Control.Concurrent ( threadDelay ) -import qualified Data.Vector.Unboxed as V( Unbox ) -import Data.Hashable ( Hashable, hash ) import Data.List ( sort ) +import Data.Hashable ( Hashable, hash, combine ) 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 Numeric.LinearAlgebra ( Element ) +import Data.Array.Repa ( Shape, DIM0, DIM1, DIM2 ) +import Control.Monad.State ( State, get, put ) -import Dvda.GExpr ( GExpr(..), getChildren ) +import Dvda.Expr ( Expr(..), Const(..), Sym(..), dim ) -type SymSet a = HS.HashSet (GExpr a) +--------------------- dynamic Expr stuff --------------------------- +data DynamicExpr a = DynamicExpr0 (Expr DIM0 a) + | DynamicExpr1 (Expr DIM1 a) + | DynamicExpr2 (Expr DIM2 a) deriving Show + +asIfExpr :: (forall sh . Expr sh a -> b) -> DynamicExpr a -> b +asIfExpr f (DynamicExpr0 e) = f e +asIfExpr f (DynamicExpr1 e) = f e +asIfExpr f (DynamicExpr2 e) = f e + +instance (Element a, Hashable a) => Hashable (DynamicExpr a) where + hash (DynamicExpr0 expr) = 0 `combine` (hash expr) + hash (DynamicExpr1 expr) = 1 `combine` (hash expr) + hash (DynamicExpr2 expr) = 2 `combine` (hash expr) + +deriving instance (Eq a, Element a) => Eq (DynamicExpr a) + +type SymSet a = HS.HashSet (DynamicExpr 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 + (HM.HashMap (DynamicExpr a) (FgNode a)) -- main lookup + (IM.IntMap (DynamicExpr a)) -- internal for reverse lookup + b + c -- 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 +instance (Hashable a, Hashable b, Hashable c, Element a) => Hashable (FunGraph a b c) where + hash (FunGraph _ im inskeys outskeys) = hash (IM.toList im, inskeys, outskeys) -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 +class Shape sh => DvdaDim sh where + makeDynamic :: Expr sh a -> DynamicExpr a + fromDynamic :: sh -> DynamicExpr a -> Expr sh a -fgGExprFromKey :: (Eq a, Hashable a, V.Unbox a) => Key -> FunGraph a b c -> Maybe (GExpr a) -fgGExprFromKey k (FunGraph _ im _ _) = IM.lookup k im +instance DvdaDim DIM0 where + makeDynamic = DynamicExpr0 + fromDynamic _ (DynamicExpr0 expr) = expr + fromDynamic _ _ = error "DIM0: fromDynamic error" +instance DvdaDim DIM1 where + makeDynamic = DynamicExpr1 + fromDynamic _ (DynamicExpr1 expr) = expr + fromDynamic _ _ = error "DIM1: fromDynamic error" +instance DvdaDim DIM2 where + makeDynamic = DynamicExpr2 + fromDynamic _ (DynamicExpr2 expr) = expr + fromDynamic _ _ = error "DIM2: fromDynamic error" -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 +fgLookup :: (Eq a, Hashable a, Element a, DvdaDim sh) => Expr sh a -> FunGraph a b c -> Maybe (FgNode a) +fgLookup (ERef sh k) fg = fgReverseLookup sh k fg +fgLookup expr (FunGraph hm _ _ _) = HM.lookup (makeDynamic expr) hm + +fgReverseLookup :: (Eq a, Hashable a, Element a, DvdaDim sh) => sh -> Key -> FunGraph a b c -> Maybe (FgNode a) +fgReverseLookup sh k fg = do + expr <- fgExprFromKey sh k fg + fgLookup expr fg + +fgExprFromKey :: DvdaDim sh => sh -> Key -> FunGraph a b c -> Maybe (Expr sh a) +fgExprFromKey sh k (FunGraph _ im _ _) = fmap (fromDynamic sh) (IM.lookup k im) + + +symSet :: (Eq a, Hashable a, Element a, DvdaDim sh) => + FunGraph a b c -> Expr sh a -> HS.HashSet (DynamicExpr a) +symSet fg e@(ESym sh (SymDependent _ _ dep)) = HS.union (HS.singleton (makeDynamic e)) (symSet fg (ESym sh dep)) +symSet _ e@(ESym _ _) = HS.singleton (makeDynamic e) +symSet fg (ERef sh k) = snd $ fromJust $ fgReverseLookup sh k fg +symSet _ (EDimensionless _) = HS.empty +symSet _ (EConst _) = HS.empty +symSet fg (EUnary _ x) = symSet fg x +symSet fg (EBinary _ x y) = (symSet fg x) `HS.union` (symSet fg y) +symSet fg (EScale x y) = (symSet fg x) `HS.union` (symSet fg y) +symSet _ (EDeriv _ _) = error "don't take symSet of EDeriv" +symSet _ (EGrad _ _) = error "don't take symSet of EGrad" +symSet _ (EJacob _ _) = error "don't take symSet of EJacob" + +-- | Try to insert the Expr into the hashmap performing CSE. +-- If the Expr is not yet in the map, insert it and return new key. +-- Otherwise don't insert, just return existing key. +insert :: (Hashable a, Eq a, Element a, DvdaDim sh) => Expr sh a -> State (FunGraph a b c) (Expr sh a) +insert (ERef _ _) = error "don't insert ERef into graph, ya goon" +insert (EConst _) = error "don't insert EConst into graph, ya goon" +insert expr = do + let dexpr = makeDynamic expr + fg@(FunGraph hm im ins outs) <- get + case fgLookup expr fg of + Just (k',_) -> return (ERef (dim expr) k') + Nothing -> do let k = HM.size hm + hm' = HM.insert dexpr (k, symSet fg expr) hm + im' = IM.insert k dexpr im + put (FunGraph hm' im' ins outs) + return (ERef (dim expr) k) + + +funGraphSummary :: (Show a, Element a, Show b, Show c) => FunGraph a b c -> String +funGraphSummary (FunGraph hm _ b c) = + init $ unlines [ "inputs: " ++ show b + , "outputs: " ++ show c , "number of nodes: " ++ show (HM.size hm) - , "graph: " ++ show hm ] +showNodes :: (Show a, Element a) => FunGraph a b c -> String +showNodes (FunGraph _ im _ _) = init $ unlines (map show (IM.toList im)) + -- 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:" +funGraphSummary' :: (Show a, Element a, Show b, Show c) => FunGraph a b c -> String +funGraphSummary' fg@(FunGraph _ im _ _) = + init $ unlines $ [ "graph:" , init $ unlines (map show (IM.toList im)) - , "outputs:" - , init $ unlines (map (show . (\k -> fromJust (IM.lookup k im))) ckeys) - ] + , "" + ] ++ [funGraphSummary fg] -collisions :: (Hashable a, V.Unbox a) => FunGraph a b c -> (Int, Int, Double) +collisions :: (Hashable a, Element 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 @@ -89,26 +163,55 @@ countCollisions n [_] = n countCollisions n [] = n -showCollisions :: (Hashable a, V.Unbox a) => FunGraph a b c -> String +showCollisions :: (Hashable a, Element 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) +emptyFunGraph = FunGraph HM.empty IM.empty inerr outerr where inerr = error "must specify inputs" outerr = error "must specify outputs" -previewGraph :: Show a => FunGraph a b c -> IO () +previewGraph :: (Show a, Element 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 +toFGLGraph :: FunGraph a b c -> Gr (DynamicExpr a) String +toFGLGraph (FunGraph hm _ _ _) = 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 + lnodes = map (\(x,(y,_)) -> (y,x)) $ HM.toList hm +-- lnodes = IM.toList im + ledges = concatMap (\(k,ge) -> map (\ch -> (ch,k,"")) (asIfExpr gc ge)) lnodes + where + gc :: Expr sh a -> [Key] + gc (EBinary _ x y) = gc x ++ gc y + gc (EUnary _ x) = gc x + gc (ERef _ k) = [k] + gc (ESym _ _) = [] + gc (EDimensionless _) = [] + gc (EScale x y) = gc x ++ gc y + gc (EConst _) = [] + gc (EDeriv _ _) = error "don't call getChildren on EDeriv" + gc (EJacob _ _) = error "don't call getChildren on EJacob" + gc (EGrad _ _) = error "don't call getChildren on EGrad" + + +instance (Show a, Element a) => Labellable (DynamicExpr a) where + toLabelValue (DynamicExpr0 e) = tlv e + toLabelValue (DynamicExpr1 e) = tlv e + toLabelValue (DynamicExpr2 e) = tlv e + +tlv :: (Show a, Shape sh, Element a) => Expr sh a -> Data.GraphViz.Attributes.Complete.Label +tlv (EBinary op _ _) = toLabelValue $ show op +tlv (EUnary op _) = toLabelValue $ show op +tlv s@(ESym _ _) = toLabelValue (show s) +tlv (EScale {}) = toLabelValue "scale" +tlv (EConst (CSingleton _ c)) = toLabelValue $ show c +tlv (EConst (CVec _ _)) = toLabelValue "vec" +tlv (EConst (CMat _ _)) = toLabelValue "mat" +tlv (EConst (CTensor _ _)) = toLabelValue "tensor" +tlv _ = error "don't try to preview one of those, ya goon" diff --git a/Dvda/HSBuilder.hs b/Dvda/HSBuilder.hs deleted file mode 100644 --- a/Dvda/HSBuilder.hs +++ /dev/null @@ -1,116 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/Dvda/HSSyntax.hs +++ /dev/null @@ -1,148 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/Dvda/HomoDim.hs +++ /dev/null @@ -1,36 +0,0 @@ -{-# 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/MultipleShooting/MSCoctave.hs b/Dvda/MultipleShooting/MSCoctave.hs new file mode 100644 --- /dev/null +++ b/Dvda/MultipleShooting/MSCoctave.hs @@ -0,0 +1,273 @@ +{-# OPTIONS_GHC -Wall #-} +{-# Language FlexibleContexts #-} +{-# Language TypeFamilies #-} + +module Dvda.MultipleShooting.MSCoctave ( msCoctave + , run + ) where + +import qualified Data.HashMap.Lazy as HM +import qualified Data.HashSet as HS +import Data.List ( zipWith6, transpose, elemIndex ) +import Data.Maybe ( fromJust, catMaybes ) + +import Dvda +import Dvda.Expr ( Expr(..), Const(..), Sym(..) ) +import Dvda.SymMonad ( rad ) +import Dvda.MultipleShooting.MSMonad +import Dvda.MultipleShooting.Types +import Dvda.OctaveSyntax ( toOctaveSource ) +import Dvda.Codegen ( writeSourceFile ) + +{- + min f(x) st: + + c(x) <= 0 + ceq(x) == 0 + A*x <= b + Aeq*x == beq + lb <= x <= ub +-} +--type Integrator a = ([Expr Z a] -> [Expr Z a] -> [Expr Z a] -> [Expr Z a] +-- -> ([Expr Z a] -> [Expr Z a] -> [Expr Z a]) +-- -> Expr Z a -> [Expr Z a]) + +type Integrator a = [Expr Z Double] + -> [Expr Z Double] + -> [Expr Z Double] + -> [Expr Z Double] + -> ([Expr Z Double] + -> [Expr Z Double] -> [Expr Z Double]) + -> Expr Z Double + -> [Expr Z Double] + +msCoctave + :: (Double ~ a) + => State (Step a) b + -> Integrator a + -> Int + -> String + -> FilePath + -> IO () +msCoctave userStep odeError n funDir name = do + _ <- writeSourceFile costSource funDir $ name ++ "_cost.m" + _ <- writeSourceFile constraintsSource funDir $ name ++ "_constraints.m" + _ <- writeSourceFile setupSource funDir $ name ++ "_setup.m" + _ <- writeSourceFile structSource funDir $ name ++ "_struct.m" + _ <- writeSourceFile unstructConstsSource funDir $ name ++ "_unstructConstants.m" + _ <- writeSourceFile unstructSource funDir $ name ++ "_unstruct.m" + _ <- writeSourceFile timeSource funDir $ name ++ "_time.m" + _ <- writeSourceFile outputSource funDir $ name ++ "_outputs.m" + _ <- writeSourceFile plotSource funDir $ name ++ "_plot.m" + _ <- writeSourceFile simSource funDir $ name ++ "_sim.m" + return () + where + steps = map (runOneStep userStep) [0..n-1] + dts = map (fromJust . stepDt) steps + + fromLeft (Left x) = x + fromLeft (Right _) = error "ERROR: fromLeft got Right" + states' = map (fst . unzip . fromJust . fromLeft . stepStates ) steps -- fromJust checked in runOneStep + actions' = map (fst . unzip . fromJust . fromLeft . stepActions) steps -- fromJust checked in runOneStep + + stateNames = map (snd . unzip . fromJust . fromLeft . stepStates ) steps -- fromJust checked in runOneStep + actionNames = map (snd . unzip . fromJust . fromLeft . stepActions) steps -- fromJust checked in runOneStep + + -- ensure that state/action (names) are the same in all steps + states = if all (head stateNames ==) stateNames + then states' + else error "ERROR: different states in different timesteps" + actions = if all (head actionNames ==) actionNames + then actions' + else error "ERROR: different actions in different timesteps" + + params = HS.toList $ foldr HS.union HS.empty (map stepParams steps) + constants = HS.toList $ foldr HS.union HS.empty (map stepConstants steps) + + boundMap = foldr HM.union HM.empty (map stepBounds steps) + + outputMap = foldl (HM.unionWith (++)) HM.empty (map stepOutputs steps) + ------------------------------------------------------------------------------------ + cost = case catMaybes $ map stepCost steps of [] -> error "need to set cost function" + cs -> sum cs + + (ceq, cineq) = foldl f ([],[]) allConstraints + where + f (eqs,ineqs) (Constraint x EQ y) = (eqs ++ [x - y], ineqs) + f (eqs,ineqs) (Constraint x LT y) = (eqs, ineqs ++ [x - y]) + f (eqs,ineqs) (Constraint x GT y) = (eqs, ineqs ++ [y - x]) + + dodeConstraints = map (Constraint (EConst (CSingleton Z 0)) EQ) $ concat $ + zipWith6 odeError (init states) (init actions) (tail states) (tail actions) + (map (execDxdt userStep) [0..]) dts + + allConstraints = dodeConstraints ++ (concatMap stepConstraints steps) ++ periodicConstraints + + periodicConstraints + | HS.size notXU > 0 = error $ "ERROR: can't set periodic constraints for non states/actions:" ++ show notXU + | otherwise = foldl g' [] $ map f' (transpose states ++ transpose actions) + where + pcSets = map stepPeriodic steps + dvSet = HS.fromList (concat states ++ concat actions) + notXU = HS.difference (foldl HS.union HS.empty pcSets) dvSet + pc0 = head pcSets + pcf = last pcSets + + -- match up states/actions by making sure they're in the same state/action list + f' xu = (HS.toList $ HS.filter (`elem` xu) pc0, HS.toList $ HS.filter (`elem` xu) pcf) + g' acc ( [], _) = acc + g' acc ( _, []) = acc + g' acc ([x], [y]) = acc ++ [Constraint x EQ y] + g' _ ( _, _) = error "ERROR: too many matching periodic constraints" + + ------------------------------------------------------------------------------------- + dvs = concat states ++ concat actions ++ params + + costFg = runFunGraph $ do + cost' <- node cost + costGrad <- rad cost' dvs + inputs_ (dvs :* constants) + outputs_ (cost' :* costGrad) + + constraintsFg = runFunGraph $ do + cineqJacob <- mapM (flip rad dvs) cineq + ceqJacob <- mapM (flip rad dvs) ceq + inputs_ (dvs :* constants) + outputs_ (cineq :* ceq :* cineqJacob :* ceqJacob) + + timeFg = runFunGraph $ do + inputs_ (dvs :* constants) + outputs_ $ init $ scanl (+) (EConst (CSingleton Z 0)) dts + + outputFg = runFunGraph $ do + inputs_ (dvs :* constants) + outputs_ $ HM.elems outputMap + + simFg = runFunGraph $ do + let x' = head states + u' = head actions + dxdt' = fromJust $ stepDxdt $ head steps + inputs_ (x' :* u' :* constants) + outputs_ dxdt' + + costSource = toOctaveSource costFg (name ++ "_cost") + constraintsSource = toOctaveSource constraintsFg (name ++ "_constraints") + outputSource = toOctaveSource outputFg (name ++ "_outputs") + timeSource = toOctaveSource timeFg (name ++ "_time") + simSource = toOctaveSource simFg (name ++ "_sim") + + (lbs, ubs, _) = unzip3 $ map getBnd dvs + where + getBnd dv = case HM.lookup dv boundMap of + Nothing -> error $ "please set bounds for " ++ show dv + Just bnd -> bnd + + setupSource = + unlines $ + [ "function [x0, Aineq, bineq, Aeq, beq, lb, ub] = "++ name ++"_setup()" + , "" +-- , "x0 = " ++ show (vectorizeDvs dvsGuess) ++ "';" + , "x0 = zeros(" ++ show (length dvs) ++ ",1);" + , "Aineq = [];" + , "bineq = [];" + , "Aeq = [];" + , "beq = [];" + , "lb = " ++ show lbs ++ "';" + , "ub = " ++ show ubs ++ "';" + ] + + -- take vector of design variables and vector of constants and return nice matlab struct + structSource = + unlines $ + ["function ret = " ++ name ++ "_struct(designVars,constants)" + , "" + , "ret.time = " ++ name ++ "_time(designVars, constants);" + , "outs = " ++ name ++ "_outputs(designVars, constants);" + , concat $ zipWith (\name' k -> "ret." ++name'++ " = outs("++show k++",:);\n") (HM.keys outputMap) [(1::Int)..] + ] ++ + toStruct dvs "designVars" (map show params) (map (\x -> [x]) params) ++ + toStruct constants "constants" (map show constants) (map (\x -> [x]) constants) + where + dvsToIdx dvs' = fromJust . flip HM.lookup (HM.fromList (zip dvs' [(1::Int)..])) + toStruct dvs' nm = zipWith (\name' vars -> "ret." ++ name' ++ " = " ++ nm ++ "(" ++ show (map (dvsToIdx dvs') vars) ++ ");\n") + + + -- take nice matlab structs and return vectors of design variables and constants + unstructSource = + unlines $ + [ "function dvs = " ++ name ++ "_unstruct(dvStruct)\n" + , "dvs = zeros(" ++ show (length dvs) ++ ", 1);" + , "" + , concatMap fromParam params + , concat $ zipWith fromXUS (head stateNames) (transpose states) + , concat $ zipWith fromXUS (head actionNames) (transpose actions) + ] + where + fromParam e@(ESym _ (Sym nm)) = + "dvs(" ++ show (1 + (fromJust $ e `elemIndex` dvs)) ++ ") = dvStruct." ++ nm ++ ";\n" + fromParam _ = error "param not ESym" + + fromXU nm e k = + "dvs(" ++ show (1 + (fromJust $ e `elemIndex` dvs)) ++ ") = dvStruct." ++ nm ++ "(" ++ show k ++ ");\n" + fromXUS name' xs = (concat $ zipWith (fromXU name') xs [(1::Int)..]) ++ "\n" + + unstructConstsSource = + unlines $ + [ "function constants = " ++ name ++ "_unstructConstants(constStruct)\n" + , "constants = zeros(" ++ show (length constants) ++ ", 1);" + , "" + , concatMap fromConst constants + ] + where + fromConst e@(ESym _ (Sym nm)) = + "constants(" ++ show (1 + (fromJust $ e `elemIndex` constants)) ++ ") = constStruct." ++ nm ++ ";\n" + fromConst _ = error "const not ESym" + + plotSource = + unlines $ + [ "function " ++ name ++ "_plot(designVars, constants)\n" + , "x = " ++ name ++ "_struct(designVars, constants);\n" + , init $ unlines $ zipWith f (HM.keys outputMap) [(1::Int)..] + ] + where + rows = ceiling $ sqrt $ (fromIntegral ::Int -> Double) $ HM.size outputMap + cols = (HM.size outputMap `div` rows) + 1 + f name' k = unlines $ + [ "subplot(" ++ show rows ++ "," ++ show cols ++ ","++show k++");" + , "plot( x.time, x." ++ name' ++ " );" + , "xlabel('time');" + , "ylabel('" ++ name'' ++ "');" + , "title('" ++ name'' ++ "');" + ] + where + name'' = foldl (\acc x -> if x == '_' then acc ++ "\\_" else acc ++ [x]) "" name' + + +spring :: State (Step Double) () +spring = do + [x, v] <- setStates ["x","v"] + [u] <- setActions ["u"] + [k, b] <- addConstants ["k", "b"] + setDxdt [v, -k*x - b*v + u] + setDt 0.1 + let cost = 2*x*x + 3*v*v + 10*u*u + setCost cost + addOutput cost "cost" + + setBound x (5,5) (TIMESTEP 0) + setBound v (0,0) (TIMESTEP 0) + + setBound x (-5,5) ALWAYS + setBound v (-10,10) ALWAYS + setBound u (-200, 100) ALWAYS + + setBound v (0,0) (TIMESTEP (n'-1)) + + setPeriodic x + +n' :: Int +n' = 20 + +run :: IO () +run = msCoctave spring simpsonsRuleError' n' "../Documents/MATLAB/" "cartpole" +--run = msCoctave spring eulerError' n' "../Documents/MATLAB/" "cartpole" diff --git a/Dvda/MultipleShooting/MSMonad.hs b/Dvda/MultipleShooting/MSMonad.hs new file mode 100644 --- /dev/null +++ b/Dvda/MultipleShooting/MSMonad.hs @@ -0,0 +1,215 @@ +{-# OPTIONS_GHC -Wall #-} +{-# Language FlexibleContexts #-} + +module Dvda.MultipleShooting.MSMonad ( State + , setStates + , setActions + , addParam + , addParams + , addConstant + , addConstants + , setDxdt + , setCost + , setDt + , addOutput + , getTimeStep + , setPeriodic + , addConstraint + , setBound + , runOneStep + , execDxdt + ) where + +import Data.Array.Repa ( Z(..) ) +import Data.Hashable ( Hashable ) +import qualified Data.HashMap.Lazy as HM +import qualified Data.HashSet as HS +import Data.List ( nub, sort ) --, union ) +import Data.Maybe ( isJust, isNothing ) +import Control.Monad ( when, zipWithM_ ) +import Control.Monad.State ( State ) +import qualified Control.Monad.State as State +--import Debug.Trace ( trace ) +--import Numeric.LinearAlgebra ( Element ) +import Text.Printf ( printf ) + +import Dvda ( sym ) +import Dvda.Expr ( Expr(..) ) +import Dvda.MultipleShooting.Types + +failDuplicates :: [String] -> [String] +failDuplicates names + | length names == length (nub names) = names + | otherwise = error $ "ERROR: saw duplicate names in: " ++ show (sort names) + +checkOctaveName :: String -> String +checkOctaveName name + | any (`elem` "\"'~!@#$%^&*()+`-=[]{}\\|;:,.<>/?") name = + error $ "ERROR: addOutput saw illegal octave variable character in string: \"" ++ name ++ "\"" + | otherwise = name + +setStates :: [String] -> State (Step a) [Expr Z a] +setStates names' = do + step <- State.get + case stepStates step of (Right states) -> return states + (Left (Just _)) -> error "states already set, don't call setStates twice" + (Left Nothing) -> do + let names = failDuplicates (map checkOctaveName names') + syms = map (sym . (++ "_" ++ show (stepIdx step))) (failDuplicates names) + State.put $ step {stepStates = Left (Just (zip syms names))} + zipWithM_ addOutput syms names + return syms + +setActions :: [String] -> State (Step a) [Expr Z a] +setActions names' = do + step <- State.get + case stepActions step of (Right actions) -> return actions + (Left (Just _)) -> error "actions already set, don't call setActions twice" + (Left Nothing) -> do + let names = failDuplicates (map checkOctaveName names') + syms = map (sym . (++ "_" ++ show (stepIdx step))) (failDuplicates names) + State.put $ step {stepActions = Left (Just (zip syms names))} + zipWithM_ addOutput syms names + return syms + +addParam :: (Eq (Expr Z a), Hashable (Expr Z a)) => String -> State (Step a) (Expr Z a) +addParam name = do + [blah] <- addParams [name] + return blah + +addConstant :: (Eq (Expr Z a), Hashable (Expr Z a)) => String -> State (Step a) (Expr Z a) +addConstant name = do + [blah] <- addParams [name] + return blah + +addParams :: (Eq (Expr Z a), Hashable (Expr Z a)) => [String] -> State (Step a) [Expr Z a] +addParams names = do + step <- State.get + let syms = map (sym . checkOctaveName) names + params0 = stepParams step + State.put $ step {stepParams = HS.union params0 (HS.fromList syms)} + return syms + +addConstants :: (Eq (Expr Z a), Hashable (Expr Z a)) => [String] -> State (Step a) [Expr Z a] +addConstants names = do + step <- State.get + let syms = map (sym . checkOctaveName) names + constants0 = stepConstants step + State.put $ step {stepConstants = HS.union constants0 (HS.fromList syms)} + return syms + +addOutput :: Expr Z a -> String -> State (Step a) () +addOutput var name = do + step <- State.get + let hm = stepOutputs step + err = error $ "ERROR: already have an output with name: \"" ++ name ++ "\"" + State.put $ step {stepOutputs = HM.insertWith err (checkOctaveName name) [var] hm} + +setDt :: Expr Z a -> State (Step a) () +setDt expr = do + step <- State.get + when (isJust (stepDt step)) $ error "dt already set, don't call setDt twice" + State.put $ step {stepDt = Just expr} + +getTimeStep :: State (Step a) Int +getTimeStep = do + step <- State.get + return (stepIdx step) + +setPeriodic :: (Eq (Expr Z a), Hashable (Expr Z a)) => Expr Z a -> State (Step a) () +setPeriodic var = do + step <- State.get + State.put $ step {stepPeriodic = HS.insert var (stepPeriodic step)} + +------------------------------------------- + +setDxdt :: [Expr Z a] -> State (Step a) () +setDxdt vars = do + step <- State.get + when (isJust (stepDxdt step)) $ error "dxdt already set, don't call setDxdt twice" + State.put $ step {stepDxdt = Just vars} + +setCost :: Expr Z a -> State (Step a) () +setCost var = do + step <- State.get + when (isJust (stepCost step)) $ error "cost already set, don't call setCost twice" + State.put $ step {stepCost = Just var} + +setBound :: (Show a, Eq a, Show (Expr Z a), Eq (Expr Z a), Hashable (Expr Z a)) + => Expr Z a -> (a, a) -> BCTime -> State (Step a) () +setBound var@(ESym _ _) (lb, ub) bnd = do + step <- State.get + let k = stepIdx step + newbnd = (lb,ub,bnd) + oldBounds = stepBounds step + + err old = error $ printf "ERROR: setBound called twice on %s (old bound: %s, new bound: %s)" (show var) (show old) (show newbnd) + + let putNewBnd = case bnd of + (TIMESTEP j) -> if j /= k + then Nothing + else case (HM.lookup var (stepBounds step)) of + Just oldbnd@(_, _, TIMESTEP _) -> err oldbnd + _ -> Just newbnd + ALWAYS -> case (HM.lookup var (stepBounds step)) of + Just oldbnd@(_,_,ALWAYS) -> err oldbnd + Just (_,_,TIMESTEP _) -> Nothing + Nothing -> Just newbnd + + when (isJust putNewBnd) $ + State.put $ step {stepBounds = HM.insert var newbnd oldBounds} +setBound _ _ _ = do + -- if execDxdt has put the x/u, they won't be symbolic - ignore them + step <- State.get + case stepStates step of + Left _ -> error "WARNING - setBound called on non-design variable, use addConstraint instead" + _ -> return () + + +addConstraint :: Expr Z a -> Ordering -> Expr Z a -> State (Step a) () +addConstraint x ordering y = + State.state (\step -> ((), step {stepConstraints = (stepConstraints step) ++ [Constraint x ordering y]})) + + +runOneStep :: State (Step a) b -> Int -> Step a +runOneStep userStep k + | isNothing (stepDxdt ret) = error "ERROR: need to set dxdt" + | isNothing (stepDt ret) = error "ERROR: need to set timestep dt" + | otherwise = stateErr `seq` actionErr `seq` ret + where + stateErr = case stepStates ret of Left Nothing -> error "ERROR: need to set states" + _ -> () + actionErr = case stepActions ret of Left Nothing -> error "ERROR: need to set actions" + _ -> () + ret = State.execState userStep $ Step { stepStates = Left Nothing + , stepActions = Left Nothing + , stepDxdt = Nothing + , stepCost = Nothing + , stepDt = Nothing + , stepBounds = HM.empty + , stepConstraints = [] + , stepParams = HS.empty + , stepConstants = HS.empty + , stepIdx = k + , stepOutputs = HM.empty + , stepPeriodic = HS.empty + } + +execDxdt :: Num (Expr Z a) => State (Step a) b -> Int -> [Expr Z a] -> [Expr Z a] -> [Expr Z a] +execDxdt userStep k x u = case stepDxdt $ State.execState userStep step0 of + Nothing -> error "ERROR: need to set dxdt" + Just dxdt -> dxdt + where + step0 = Step { stepStates = Right x + , stepActions = Right u + , stepDxdt = Nothing + , stepDt = Nothing + , stepCost = Nothing + , stepBounds = HM.empty + , stepConstraints = [] + , stepParams = HS.empty + , stepConstants = HS.empty + , stepIdx = k + , stepOutputs = HM.empty + , stepPeriodic = HS.empty + } diff --git a/Dvda/MultipleShooting/MultipleShooting.hs b/Dvda/MultipleShooting/MultipleShooting.hs new file mode 100644 --- /dev/null +++ b/Dvda/MultipleShooting/MultipleShooting.hs @@ -0,0 +1,166 @@ +{-# OPTIONS_GHC -Wall #-} +{-# Language FlexibleContexts #-} + +module Dvda.MultipleShooting.MultipleShooting ( Cost(..) + , MultipleShooting(..) + , Constraint'(..) + , DesignVars(..) + , multipleShooting + , simpleSystem + , boundEqs + , boundEq + , boundInterval + , boundIntervals + , ltZero + , replaceFinalCost + , vectorizeDvs + , dvIdx + , numDvs + , interpolateInitialGuess + , simpsonsRuleError + , eulerError + ) where + +import Text.Printf ( printf ) +import Data.List ( elemIndex, zipWith6 ) +import Data.Maybe ( isJust, fromJust ) +import Data.Array.Repa ( Z(..) ) +import Debug.Trace ( trace ) + +import Dvda ( svec ) +import Dvda.Expr ( Expr(..), Sym(..) ) +import Dvda.SparseLA ( SparseVec, svCats, svZeros, svSize, sparseListFromSv, svFromList ) +import Dvda.MultipleShooting.Types + +data Cost a = Cost (SparseVec (Expr Z a) -> SparseVec (Expr Z a) -> Expr Z a) (Int,Int) + +data Constraint' a = Constraint' Ordering (SparseVec (Expr Z a)) (SparseVec (Expr Z a)) deriving Show + +data System a = System { sysOdes :: [Ode a] + , sysCosts :: [Cost a] + , sysDts :: [Expr Z a] + } + +data DesignVars a = DesignVars { dvStates :: [SparseVec a] + , dvActions :: [SparseVec a] + , dvParams :: [a] + } + +data MultipleShooting a = MultipleShooting { msSystem :: System a + , msDesignVars :: DesignVars (Expr Z a) + , msDodeConstraints :: [Constraint' a] + , msConstants :: [Expr Z a] + , msObjFun :: Expr Z a + } + + + +data Bound a = Bound { boundVar :: Expr Z a + , boundL :: a + , boundU :: a + } + +instance Show a => Show (Bound a) where + show bound = show (boundL bound) ++ " <= " ++ name ++ " <= " ++ show (boundU bound) + where + name = safeGetSymNameFromExpr (boundVar bound) + +safeGetSymNameFromExpr :: Expr sh a -> String +safeGetSymNameFromExpr (ESym _ (Sym name)) = name +safeGetSymNameFromExpr _ = trace "Warning - Bound has non-symbolic value" "{NOT A DESIGN VARIABLE}" + +numDvs :: MultipleShooting a -> Int +numDvs = length . vectorizeDvs . msDesignVars + +dvIdx :: Eq (Expr Z a) => MultipleShooting a -> Expr Z a -> Int +dvIdx ms val + | isJust idx = fromJust idx + | otherwise = error $ "Error - idxOfDvs fail" + where + idx = elemIndex val (vectorizeDvs $ msDesignVars ms) + +vectorizeDvs :: DesignVars a -> [a] +vectorizeDvs (DesignVars {dvStates = states, dvActions = actions, dvParams = params}) = + sparseListFromSv (svCats [svCats states, svCats actions]) ++ params + +--vectorizedIndices :: Multipleshooting a -> DesignVars Int +--vectorizedIndices ms +-- | any ((/=) (head odeDims)) (tail odeDims) = error "vectorizedIndices got ODE dimension mismatch" +-- | otherwise = DesignVars +-- where +-- odeDims = map (\(Ode _ _ _ d) -> d) $ sysOdes (msSystem ms) + + +boundEq :: Eq (Expr Z a) => Expr Z a -> a -> Bound a +boundEq x val = Bound { boundL = val + , boundU = val + , boundVar = x + } + +boundEqs :: Eq (Expr Z a) => SparseVec (Expr Z a) -> SparseVec a -> [Bound a] +boundEqs xs vals = zipWith boundEq (sparseListFromSv xs) (sparseListFromSv vals) + +boundInterval :: Eq (Expr Z a) => Expr Z a -> (a, a) -> Bound a +boundInterval x (lb, ub) = Bound { boundL = lb + , boundU = ub + , boundVar = x + } + +boundIntervals :: Eq (Expr Z a) => SparseVec (Expr Z a) -> [(a,a)] -> [Bound a] +boundIntervals xs bnds = zipWith boundInterval (sparseListFromSv xs) bnds + + +multipleShooting :: Fractional (Expr Z a) => System a -> [Expr Z a] -> [Expr Z a] + -> (SparseVec (Expr Z a) -> SparseVec (Expr Z a) -> SparseVec (Expr Z a) -> SparseVec (Expr Z a) -> Ode a -> Expr Z a -> SparseVec (Expr Z a)) + -> MultipleShooting a +multipleShooting sys params constants odeError + | dimensionsMatch = MultipleShooting { msSystem = sys + , msDesignVars = DesignVars { dvStates = states + , dvActions = actions + , dvParams = params + } + , msDodeConstraints = dodeConstraints + , msConstants = constants + , msObjFun = objFun + } + | otherwise = error $ printf "Error in multipleShooting: lengths of odes (%d), costs (%d), dts (%d) are not consistent" nOdes nCosts nDts + where + dimensionsMatch = (nOdes == nDts) && (nCosts == nOdes + 1) && (and $ zipWith (==) odeDims costDims) + + odeDims = map (\(Ode _ d) -> d) (sysOdes sys) + costDims = map (\(Cost _ d) -> d) (sysCosts sys) + + nOdes = length (sysOdes sys) + nCosts = length (sysCosts sys) + nDts = length (sysDts sys) + + states = zipWith (\(nx,_) k -> svec ("x_"++show k) nx) costDims [0..nCosts-1] + actions = zipWith (\(_,nu) k -> svec ("u_"++show k) nu) costDims [0..nCosts-1] + + dodeConstraints = map eqZero $ zipWith6 odeError (init states) (init actions) (tail states) (tail actions) + (sysOdes sys) (sysDts sys) + + objFun = sum $ zipWith3 (\(Cost cost _) x u -> cost x u) (sysCosts sys) states actions + +simpleSystem :: Ode a -> Cost a -> Expr Z a -> Int -> System a +simpleSystem ode cost dt n = System { sysOdes = replicate (n-1) ode + , sysCosts = replicate n cost + , sysDts = replicate (n-1) dt + } + +replaceFinalCost :: Cost a -> System a -> System a +replaceFinalCost cost sysIn = sysIn {sysCosts = init (sysCosts sysIn) ++ [cost]} + +eqZero :: SparseVec (Expr Z a) -> Constraint' a +eqZero g = Constraint' EQ g (svZeros $ svSize g) + +ltZero :: Fractional a => SparseVec (Expr Z a) -> Constraint' a +ltZero g = Constraint' LT g (svZeros $ svSize g) + +interpolateInitialGuess :: Fractional a => SparseVec a -> SparseVec a -> Int -> [SparseVec a] +interpolateInitialGuess x0 xf n' = map (combine x0 xf) alphas + where + n = fromIntegral n' + alphas = map (/ (n-1)) $ map fromIntegral [0..n'-1] + combine v0 vf alpha = + svFromList $ zipWith (\x y -> alpha*x + (1-alpha)*y) (sparseListFromSv v0) (sparseListFromSv vf) diff --git a/Dvda/MultipleShooting/Types.hs b/Dvda/MultipleShooting/Types.hs new file mode 100644 --- /dev/null +++ b/Dvda/MultipleShooting/Types.hs @@ -0,0 +1,88 @@ +{-# OPTIONS_GHC -Wall #-} +{-# Language FlexibleContexts #-} + +module Dvda.MultipleShooting.Types ( Step(..) + , Constraint(..) + , Ode(..) + , BCTime(..) + , eulerError + , simpsonsRuleError + , eulerError' + , simpsonsRuleError' + ) where + +import Data.HashMap.Lazy ( HashMap ) +import Data.HashSet ( HashSet ) + +import Dvda ( Z ) +import Dvda.Expr ( Expr(..) ) +import Dvda.SparseLA + +data BCTime = ALWAYS | TIMESTEP Int deriving (Show, Eq) + +data Constraint a = Constraint (Expr Z a) Ordering (Expr Z a) deriving Show + +data Step a = Step { stepStates :: Either (Maybe [(Expr Z a, String)]) [Expr Z a] + , stepActions :: Either (Maybe [(Expr Z a, String)]) [Expr Z a] + , stepParams :: HashSet (Expr Z a) + , stepConstants :: HashSet (Expr Z a) + , stepDxdt :: Maybe [Expr Z a] + , stepCost :: Maybe (Expr Z a) + , stepDt :: Maybe (Expr Z a) + , stepBounds :: HashMap (Expr Z a) (a,a, BCTime) + , stepConstraints :: [Constraint a] + , stepIdx :: Int + , stepOutputs :: HashMap String [Expr Z a] + , stepPeriodic :: HashSet (Expr Z a) + } + +data Ode a = Ode (SparseVec (Expr Z a) -> SparseVec (Expr Z a) -> SparseVec (Expr Z a)) (Int,Int) + +wrapOdeError :: Fractional (Expr Z a) + => (SparseVec (Expr Z a) -> SparseVec (Expr Z a) -> SparseVec (Expr Z a) -> SparseVec (Expr Z a) -> Ode a -> Expr Z a -> SparseVec (Expr Z a)) + -> [Expr Z a] -> [Expr Z a] -> [Expr Z a] -> [Expr Z a] + -> ([Expr Z a] -> [Expr Z a] -> [Expr Z a]) + -> Expr Z a + -> [Expr Z a] +wrapOdeError odeError xk uk xkp1 ukp1 dxdt dt = + denseListFromSv $ odeError xk' uk' xkp1' ukp1' (Ode dxdt' (error "FUUUUCK")) dt + where + xk' = svFromList xk + xkp1' = svFromList xkp1 + uk' = svFromList uk + ukp1' = svFromList ukp1 + dxdt' x u = svFromList $ dxdt (denseListFromSv x) (denseListFromSv u) + +eulerError' :: Fractional (Expr Z a) + => [Expr Z a] -> [Expr Z a] -> [Expr Z a] -> [Expr Z a] + -> ([Expr Z a] -> [Expr Z a] -> [Expr Z a]) + -> Expr Z a + -> [Expr Z a] +eulerError' = wrapOdeError eulerError + +simpsonsRuleError' :: Fractional (Expr Z a) + => [Expr Z a] -> [Expr Z a] -> [Expr Z a] -> [Expr Z a] + -> ([Expr Z a] -> [Expr Z a] -> [Expr Z a]) + -> Expr Z a + -> [Expr Z a] +simpsonsRuleError' = wrapOdeError simpsonsRuleError + +eulerError :: Fractional (Expr Z a) => SparseVec (Expr Z a) -> SparseVec (Expr Z a) -> SparseVec (Expr Z a) -> SparseVec (Expr Z a) -> Ode a -> Expr Z a -> SparseVec (Expr Z a) +eulerError xk uk xkp1 _ (Ode ode _) dt = xkp1 - (xk + svScale dt f0) + where + f0 = ode xk uk + +simpsonsRuleError :: Fractional (Expr Z a) => SparseVec (Expr Z a) -> SparseVec (Expr Z a) -> SparseVec (Expr Z a) -> SparseVec (Expr Z a) -> Ode a -> Expr Z a -> SparseVec (Expr Z a) +simpsonsRuleError xk uk xkp1 ukp1 (Ode ode _) dt = xkp1 - xk - (svScale (dt/6.0) (f0 + fourFm + f1)) + where + f0 = ode xk uk + f1 = ode xkp1 ukp1 + + um = svScale 0.5 (uk + ukp1) + xm = xm' - xm'' + where + xm' = svScale 0.5 (xk + xkp1) + xm'' = svScale (0.125 * dt) (f1 - f0) + + fm = ode xm um + fourFm = svScale 4 fm diff --git a/Dvda/OctaveSyntax.hs b/Dvda/OctaveSyntax.hs new file mode 100644 --- /dev/null +++ b/Dvda/OctaveSyntax.hs @@ -0,0 +1,165 @@ +{-# OPTIONS_GHC -Wall #-} +{-# Language GADTs #-} +{-# Language FlexibleInstances #-} +{-# Language TypeOperators #-} + +module Dvda.OctaveSyntax ( GenOctave + , toOctaveSource + ) where + +import Data.Maybe ( fromJust ) +import Data.List ( intersperse ) +import Data.IntMap ( Key ) +import qualified Data.IntMap as IM +import Numeric.LinearAlgebra ( Element ) +import Text.Printf + +import Dvda ( DIM0 ) +import Dvda.Expr ( Expr(..), Const(..), isVal ) +import Dvda.Graph ( FunGraph(..), DynamicExpr, asIfExpr ) +import Dvda.BinUn ( BinOp(..), UnOp(..) ) +import Dvda.SymMonad ( (:*)(..) ) +import qualified Dvda.Config as Config + +class GenOctave a where + numObjects :: a -> Int + writeOutputs :: a -> Int -> String + writeInputs :: a -> Int -> (String, IM.IntMap String) + +instance GenOctave (Expr DIM0 Double) where + numObjects _ = 1 + writeOutputs e outputK = printf "%% output %d\noutput%d = %s; %% Expr DIM0 Double\n" outputK outputK (writeExpr e) + writeInputs e@(ERef _ k) inputK = (printf "%% input %d\n%s\n" inputK decl, IM.singleton k decl) + where + decl = printf "%s = x%d;" (writeExpr e) inputK + writeInputs e inputK = error $ "input " ++ show inputK ++ " is non-symbolic: " ++ show e + + +instance GenOctave [Expr DIM0 Double] where + numObjects _ = 1 + writeOutputs [] outputK = printf "%% output %d\noutput%d = []; %% [Expr DIM0 Double]\n" outputK outputK; + writeOutputs exprs outputK = + unlines $ + (printf "%% output %d\noutput%d = zeros(%d,1); %% [Expr DIM0 Double]" outputK outputK (length exprs)): + zipWith f [(1::Int)..] exprs + where + f outIdx e = printf "%soutput%d(%d) = %s;" maybeComment outputK outIdx (writeExpr e) + where + maybeComment + | isVal 0 e = "% " + | otherwise = "" + writeInputs exprs inputK = ((printf "%% input %d\n" inputK) ++ unlines (map snd keyDecls), IM.fromList keyDecls) + where + keyDecls = zipWith f [(1::Int)..] exprs + f outIdx e@(ERef _ k) = (k, printf "%s = x%d(%d);" (writeExpr e) inputK outIdx) + f outIdx e = error $ "input " ++ show inputK ++ ", " ++ show outIdx ++" is non-symbolic: " ++ show e + +instance GenOctave [[Expr DIM0 Double]] where + numObjects _ = 1 + writeOutputs [] outputK = printf "%% output %d\noutput%d = []; %% [[Expr DIM0 Double]]\n" outputK outputK; + writeOutputs exprs outputK = + unlines $ + (printf "output%d = zeros(%d,%d); %% [[Expr DIM0 Double]]" outputK (length exprs) (length (head exprs))): + zipWith f [(r,c) | r <- [1..length exprs], c <- [1..(length (head exprs))]] (concat exprs) + where + f (rowIdx,colIdx) e = printf "%soutput%d(%d,%d) = %s;" maybeComment outputK rowIdx colIdx (writeExpr e) + where + maybeComment + | isVal 0 e = "% " + | otherwise = "" + writeInputs exprs inputK = ((printf "% input %d\n" inputK) ++ unlines (map snd keyDecls), IM.fromList keyDecls) + where + keyDecls = zipWith f [(r,c) | r <- [1..length exprs], c <- [1..(length (head exprs))]] (concat exprs) + f (rowIdx,colIdx) e@(ERef _ k) = (k, printf "%s = x%d(%d,%d);" (writeExpr e) inputK rowIdx colIdx) + f outIdx e = error $ "input " ++ show inputK ++ ", " ++ show outIdx ++" is non-symbolic: " ++ show e + +instance (GenOctave a, GenOctave b) => GenOctave (a :* b) where + numObjects (x :* y) = numObjects x + numObjects y + writeOutputs (x :* y) outputK = writeOutputs x outputK ++ "\n" ++ writeOutputs y (outputK + numObjects x) + writeInputs (x :* y) inputK = (headerX ++ '\n' : headerY, IM.unionWith err imx imy) + where + (headerX,imx) = writeInputs x inputK + (headerY,imy) = writeInputs y (inputK + numObjects x) + err = error "writeInputs (x :* y) got inputs from two sources" + +-- assign a scalar +sassign :: Key -> String +sassign k = Config.nameHSVar k ++ " = " + +octaveBinary :: BinOp -> String +octaveBinary Add = "+" +octaveBinary Sub = "-" +octaveBinary Mul = "*" +octaveBinary Div = "/" +octaveBinary Pow = "^" +octaveBinary LogBase = "error('no logBase here lol')" + +octaveUnary :: UnOp -> String +octaveUnary Abs = "abs" +octaveUnary Neg = "-" +octaveUnary Signum = "sign" +octaveUnary Exp = "exp" +octaveUnary Sqrt = "sqrt" +octaveUnary Log = "log" +octaveUnary Sin = "sin" +octaveUnary Cos = "cos" +octaveUnary Tan = "tan" +octaveUnary ASin = "asin" +octaveUnary ACos = "acos" +octaveUnary ATan = "atan" +octaveUnary Sinh = "sinh" +octaveUnary Cosh = "cosh" +octaveUnary Tanh = "tanh" +octaveUnary ASinh = "asinh" +octaveUnary ATanh = "atanh" +octaveUnary ACosh = "acosh" + +writeExpr :: (Show a, Element a) => Expr sh a -> String +writeExpr (ERef _ k) = Config.nameHSVar k +writeExpr (EBinary op x y) = writeExpr x ++ " " ++ octaveBinary op ++ " " ++ writeExpr y +writeExpr (EUnary op x) = octaveUnary op ++ "( " ++ writeExpr x ++ " )" +writeExpr (EScale x y) = "LA.scale " ++ writeExpr x ++ " " ++ writeExpr y +writeExpr (EConst (CSingleton _ x)) = show x +writeExpr (EConst (CVec _ x)) = show x -- Config.nameHSConst k +writeExpr (EConst (CMat _ x)) = show x -- Config.nameHSConst k +writeExpr (EConst (CTensor _ x)) = show x -- Config.nameHSConst k +writeExpr (ESym _ _) = error "ESym shouldn't be handled here" +writeExpr (EDimensionless _) = error "EDimensionless shouldn't be handled here" +writeExpr (EJacob _ _) = error "EJacob shouldn't be handled here" +writeExpr (EDeriv _ _) = error "EDeriv shouldn't be handled here" +writeExpr (EGrad _ _) = error "EGrad shouldn't be handled here" + +writeAssignment :: (Show a, Element a) => IM.IntMap String -> (Key, DynamicExpr a) -> (String, String) +writeAssignment inputMap (k, dexpr) + | asIfExpr isSym dexpr = (fromJust $ IM.lookup k inputMap, drop 13 (show dexpr)) + | otherwise = (sassign k ++ asIfExpr writeExpr dexpr ++ ";", drop 13 (show dexpr)) + where + isSym (ESym _ _) = True + isSym _ = False + + +toOctaveSource :: (Show a, Element a, GenOctave b, GenOctave c) => + FunGraph a b c -> String -> String +toOctaveSource (FunGraph _ im inputs outputs) funName = + unlines $ + [ "function [" ++ outputHeader ++ "] = " ++ funName ++ "(" ++ inputHeader ++ ")" + , "" + , "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% inputs: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" + , unlines $ map ('%' :) $ lines inputDecls + , "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% body: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" + , unlines $ (zipWith3 (\d s c -> d ++ s ++ "% " ++ c) decls extraSpaces comments) + , "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% outputs: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" + , writeOutputs outputs 0 + , "end" + ] + where + header prefix num = concat $ intersperse "," $ map (\k -> prefix ++ show k) [0..(num - 1)] + inputHeader = header "x" (numObjects inputs) + outputHeader = header "output" (numObjects outputs) + (decls, comments) = unzip $ map (writeAssignment inputMap) (IM.toList im) + + (inputDecls, inputMap) = writeInputs inputs 0 + + lengths = map length decls + longestDecl = maximum lengths + extraSpaces = map (\n -> replicate (longestDecl - n + 4) ' ') lengths diff --git a/Dvda/SparseLA.hs b/Dvda/SparseLA.hs new file mode 100644 --- /dev/null +++ b/Dvda/SparseLA.hs @@ -0,0 +1,245 @@ +{-# OPTIONS_GHC -Wall #-} + +module Dvda.SparseLA ( SparseVec + , SparseMat + , svFromList + , smFromLists + , svFromSparseList + , smFromSparseList + , denseListFromSv + , sparseListFromSv + , svZeros + , smZeros + , svSize + , smSize + , svMap + , smMap + , svBinary + , smBinary + , svAdd + , svSub + , svMul + , smAdd + , smSub + , smMul + , svScale + , smScale + , getRow + , getCol + , svCat + , svCats + , sVV + , sMV + ) where + +import Data.List ( foldl' ) +import Data.Maybe ( fromJust, fromMaybe ) --, isNothing ) +--import qualified Data.Traversable as T +import Data.IntMap ( IntMap ) +import qualified Data.IntMap as IM + +-- map from row to (map from col to value) +data SparseMat a = SparseMat (Int,Int) (IntMap (IntMap a)) + +instance Show a => Show (SparseMat a) where + show (SparseMat rowsCols xs) = "SparseMat " ++ show vals ++ " " ++ show rowsCols + where + vals = concatMap f (IM.toList xs) + f (row,m) = map g (IM.toList m) + where + g (col, val) = ((row, col), val) + +instance Num a => Num (SparseMat a) where + x + y = fromJust $ smAdd x y + x - y = fromJust $ smSub x y + x * y = fromJust $ smMul x y + abs = smMap abs + signum = smMap signum + fromInteger = error "fromInteger not declared for Num SparseMat" + +-- puts zeroes where there aren't entries +denseListFromSv :: Num a => SparseVec a -> [a] +denseListFromSv v@(SparseVec _ im) = IM.elems $ IM.union im (IM.fromList $ zip [0..n-1] (repeat 0)) + where + n = svSize v + +sparseListFromSv :: SparseVec a -> [a] +sparseListFromSv (SparseVec _ im) = IM.elems im + +svZeros :: Int -> SparseVec a +svZeros n = SparseVec n IM.empty + +smZeros :: (Int, Int) -> SparseMat a +smZeros rowsCols = SparseMat rowsCols IM.empty + +smSize :: SparseMat a -> (Int,Int) +smSize (SparseMat rowsCols _) = rowsCols + +smMap :: (a -> b) -> SparseMat a -> SparseMat b +smMap f (SparseMat sh maps) = SparseMat sh (IM.map (IM.map f) maps) + +smFromLists :: [[a]] -> SparseMat a +smFromLists blah = smFromSparseList sparseList (rows, cols) + where + rows = length blah + cols = length (head blah) + sparseList = concat $ zipWith (\row xs -> zipWith (\col x -> ((row,col),x)) [0..] xs) [0..] blah + +smFromSparseList :: [((Int,Int),a)] -> (Int,Int) -> SparseMat a +smFromSparseList xs' rowsCols = SparseMat rowsCols (foldr f IM.empty xs') + where + f ((row,col), val) = IM.insertWith g row (IM.singleton col val) + where + g = IM.union +-- g = IM.unionWith (error $ "smFromList got 2 values for entry: "++show (row,col)) + +---- more efficient using mergeWithKey, but needs containers 0.5 so wait till ghc 7.6 :( +-- smBinary :: (a -> b -> c) -> (IntMap a -> IntMap c) -> (IntMap b -> IntMap c) +-- -> SparseMat a -> SparseMat b -> Maybe (SparseMat c) +-- smBinary fBoth fLeft fRight (SparseMat shx xs) (SparseMat shy ys) +-- | shx /= shy = Nothing +-- | isNothing merged = Nothing +-- | otherwise = Just $ SparseMat shx (fromJust merged) +-- where +-- merged = T.sequence $ IM.mergeWithKey f (IM.map (Just . fLeft)) (IM.map (Just . fRight)) xs ys +-- where +-- cols = Repa.shapeOfList [head $ Repa.listOfShape shx] +-- f _ x y = case svBinary fBoth fLeft fRight (SparseVec cols x) (SparseVec cols y) of +-- Just (SparseVec _ im) -> Just (Just im) +-- Nothing -> Just Nothing + +smBinary :: (a -> a -> a) -> (IntMap a -> IntMap a) -> (IntMap a -> IntMap a) + -> SparseMat a -> SparseMat a -> Maybe (SparseMat a) +smBinary fBoth fLeft fRight (SparseMat shx@(_,cols) xs) (SparseMat shy ys) + | shx /= shy = Nothing + | otherwise = Just $ SparseMat shx merged + where + merged = IM.unionWith f (IM.map fLeft xs) (IM.map fRight ys) + where + f x y = case svBinary fBoth fLeft fRight (SparseVec cols x) (SparseVec cols y) of + Just (SparseVec _ im) -> im + Nothing -> error "goons everywhere" + +-------------------------------------------------------------------------------------- +data SparseVec a = SparseVec Int (IntMap a) + +svSize :: SparseVec a -> Int +svSize (SparseVec sh _) = sh + +instance Show a => Show (SparseVec a) where + show sv@(SparseVec _ xs) = "SparseVec " ++ show vals ++ " " ++ show rows + where + rows = svSize sv + vals = IM.toList xs + +instance Num a => Num (SparseVec a) where + x + y = fromJust $ svAdd x y + x - y = fromJust $ svSub x y + x * y = fromJust $ svMul x y + abs = svMap abs + signum = svMap signum + fromInteger = error "fromInteger not declared for Num SparseVec" + +svFromList :: [a] -> SparseVec a +svFromList xs = svFromSparseList (zip [0..] xs) (length xs) + +svFromSparseList :: [(Int,a)] -> Int -> SparseVec a +svFromSparseList xs rows = SparseVec rows (IM.fromList xs) + +svMap :: (a -> b) -> SparseVec a -> SparseVec b +svMap f (SparseVec sh maps) = SparseVec sh (IM.map f maps) + +svBinary :: (a -> b -> c) -> (IntMap a -> IntMap c) -> (IntMap b -> IntMap c) + -> SparseVec a -> SparseVec b -> Maybe (SparseVec c) +svBinary fBoth fLeft fRight (SparseVec shx xs) (SparseVec shy ys) + | shx /= shy = Nothing + | otherwise = Just $ SparseVec shx merged + where + -- more efficient using mergeWithKey, but needs containers 0.5 so wait till ghc 7.6 :( +-- merged = IM.mergeWithKey (\_ x y -> Just (fBoth x y)) fLeft fRight xs ys + merged = IM.unionWithKey f (fLeft xs) (fRight ys) + where + f k _ _ = fBoth (fromJust $ IM.lookup k xs) (fromJust $ IM.lookup k ys) + +--------------------------------------------------------------------------- +svAdd :: Num a => SparseVec a -> SparseVec a -> Maybe (SparseVec a) +svAdd = svBinary (+) id id + +svSub :: Num a => SparseVec a -> SparseVec a -> Maybe (SparseVec a) +svSub = svBinary (-) id (IM.map negate) + +svMul :: Num a => SparseVec a -> SparseVec a -> Maybe (SparseVec a) +svMul = svBinary (*) (\_ -> IM.empty) (\_ -> IM.empty) + + +smAdd :: Num a => SparseMat a -> SparseMat a -> Maybe (SparseMat a) +smAdd = smBinary (+) id id + +smSub :: Num a => SparseMat a -> SparseMat a -> Maybe (SparseMat a) +smSub = smBinary (-) id (IM.map negate) + +smMul :: Num a => SparseMat a -> SparseMat a -> Maybe (SparseMat a) +smMul = smBinary (*) (\_ -> IM.empty) (\_ -> IM.empty) + +-------------------------------------------------------------------------- + +svScale :: Num a => a -> SparseVec a -> SparseVec a +svScale x (SparseVec sh xs) = SparseVec sh (IM.map (x *) xs) + +smScale :: Num a => a -> SparseMat a -> SparseMat a +smScale x (SparseMat sh xs) = SparseMat sh (IM.map (IM.map (x *)) xs) + + +-------------------------------------------------------------------------- +getRow :: Int -> SparseMat a -> SparseVec a +getRow row sm@(SparseMat (_,cols) xs) + | row >= (\(rows,_) -> rows) (smSize sm) = + error $ "getRow saw out of bounds index " ++ show row ++ " for matrix size " ++ show (smSize sm) + | otherwise = SparseVec cols out + where + out = fromMaybe IM.empty (IM.lookup row xs) + +getCol :: Int -> SparseMat a -> SparseVec a +getCol col sm@(SparseMat (rows,_) xs) + | col >= (\(_,cols) -> cols) (smSize sm) = + error $ "getCol saw out of bounds index " ++ show col ++ " for matrix size " ++ show (smSize sm) + | otherwise = SparseVec rows out + where + out = IM.mapMaybe (IM.lookup col) xs + +--------------------------------------------------------------------------- +sVV :: Num a => SparseVec a -> SparseVec a -> Maybe a +sVV x y = fmap (\(SparseVec _ xs) -> sum (IM.elems xs)) (svMul x y) + +sMV :: Num a => SparseMat a -> SparseVec a -> Maybe (SparseVec a) +sMV (SparseMat (mrows,mcols) ms) vec@(SparseVec vsize _) + | mcols /= vsize = Nothing + | otherwise = Just $ SparseVec mrows out + where + out = IM.mapMaybe f ms + where + f im = sVV (SparseVec mcols im) vec + +--------------------------------------------------------------------------- +svCat :: SparseVec a -> SparseVec a -> SparseVec a +svCat svx@(SparseVec _ xs) svy@(SparseVec _ ys) = SparseVec (shx + shy) (IM.union xs newYs) + where + shx = svSize svx + shy = svSize svy + newYs = IM.fromList $ map (\(k,x) -> (k+shx, x)) $ IM.toList ys + +svCats :: [SparseVec a] -> SparseVec a +svCats [] = SparseVec 0 IM.empty +svCats (xs0:xs) = foldl' svCat xs0 xs + +--mx' :: SparseMat Double +--mx' = smFromList [((0,0), 10), ((0,2), 20), ((1,0), 30)] (2,3) +-- +--my' :: SparseMat Double +--my' = smFromList [((0,0), 1), ((0,1), 7)] (2,3) +-- +--x' :: SparseVec Int +--x' = svFromList [(0,10), (1, 20)] 4 +-- +--y' :: SparseVec Int +--y' = svFromList [(0,7), (3, 30)] 4 diff --git a/Dvda/SymMonad.hs b/Dvda/SymMonad.hs --- a/Dvda/SymMonad.hs +++ b/Dvda/SymMonad.hs @@ -1,16 +1,14 @@ {-# OPTIONS_GHC -Wall #-} -{-# Language GADTs #-} -{-# Language FlexibleContexts #-} {-# Language TypeOperators #-} {-# Language TypeFamilies #-} -{-# Language MultiParamTypeClasses #-} {-# Language FlexibleInstances #-} +{-# Language FlexibleContexts #-} +{-# Language GADTs #-} +{-# Language DoAndIfThenElse #-} module Dvda.SymMonad ( (:*)(..) - , HList(..) - , Exprs + , MkFunGraph(..) , node - , node' , inputs , inputs_ , outputs @@ -19,306 +17,326 @@ , runFunGraph , rad , getSensitivities + , recover + , fullShow + , fullShowNodes + , runDeriv ) 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 Control.Monad ( foldM, liftM ) +import Control.Monad.State ( State, get, put, runState ) +import Data.Array.Repa ( DIM0, DIM1, DIM2, 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 Numeric.LinearAlgebra ( Element, Vector, Matrix ) +import qualified Numeric.LinearAlgebra as LA +-- import Debug.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 ) +import Dvda.BinUn ( applyUnary, applyBinary ) +import Dvda.Graph ( FunGraph(..), DynamicExpr(..), DvdaDim(..), insert, emptyFunGraph, fgLookup, fgExprFromKey ) +import Dvda.Expr ( Expr(..), Const(..), Sym(..), dim ) --- | 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 +---- | take all sub expressions of an Expr and turn them into nodes +---- return an Expr that is just a ref +node :: (Hashable a, Eq a, Floating a, Num (Vector a), LA.Container Vector a, DvdaDim sh) => + Expr sh a -> State (FunGraph a b c) (Expr sh a) +node (EDimensionless _) = error "don't put EDimensionless in graph, ya goon" +node (EJacob _ _) = error "can't do node EJacob yet" +node e@(ERef _ _) = return e +node e@(EConst _) = return e +node e@(ESym _ (SymDependent _ _ dep)) = do + _ <- node (ESym Z dep) + insert e +node e@(ESym _ _) = insert e +node (EUnary op x') = do x <- node x' - arg <- node arg' - outs <- rad x [arg] - node' (head outs) -node' (EGrad x' arg') = do + insert $ EUnary op x +node (EBinary op x' y') = do x <- node x' - arg <- node arg' + y <- node y' + insert $ EBinary op x y +node (EScale x' y') = do + x <- node x' + y <- node y' + insert $ EScale x y +node (EDeriv x_ arg_) = do + x <- node x_ + arg <- node arg_ outs <- rad x [arg] - node' (head outs) + 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 +-- gradient of expression w.r.t. list of args +rad :: (Eq a, Floating a, Num (Vector a), Hashable a, LA.Container Vector a, DvdaDim sh0, DvdaDim sh) => + Expr sh0 a -> [Expr sh a] -> State (FunGraph a b c) [Expr sh a] +rad expr' args' = do + expr <- node expr' + args'' <- mapM node args' 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 - + let args = map (\(ERef sh k) -> fromJust $ fgExprFromKey sh k fg) args'' + argSet = HS.fromList (map makeDynamic args) -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) + sensitivities <- getSensitivities argSet expr (EConst (CSingleton (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 + + let getSens arg = case HM.lookup (makeDynamic arg) sensitivities of + Just sens -> node $ fromDynamic (dim arg) sens +-- Nothing -> trace "WARNING: taking deriviative df/dx where f is not a function of x" $ +-- return $ EConst (CSingleton (dim arg) 0) + Nothing -> return $ EConst (CSingleton (dim arg) 0) + mapM getSens args --- 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 +-- | combine two (DynamicExpr a, DynamicExpr a) hashmaps +-- if there is a conflict, add the two sensitivities together +unionWithPlus :: (Hashable a, Eq a, Num (Vector a), LA.Container Vector a, Floating a) => + HM.HashMap (DynamicExpr a) (DynamicExpr a) -> HM.HashMap (DynamicExpr a) (DynamicExpr a) + -> State (FunGraph a b c) (HM.HashMap (DynamicExpr a) (DynamicExpr a)) +unionWithPlus xs ys = foldM addCommon union0 commonDExprs where -- the gexprs that occur in both maps - commonGExprs = HM.keys $ HM.intersection xs ys + commonDExprs = 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) - + addCommon hm commonDExpr = do + let xsens = fromJust $ HM.lookup commonDExpr xs + ysens = fromJust $ HM.lookup commonDExpr ys + xysens <- case (xsens,ysens) of + (DynamicExpr0 x, DynamicExpr0 y) -> do + ret <- node (x + y) + return (makeDynamic ret) + (DynamicExpr1 x, DynamicExpr1 y) -> do + ret <- node (x + y) + return (makeDynamic ret) + (DynamicExpr2 x, DynamicExpr2 y) -> do + ret <- node (x + y) + return (makeDynamic ret) + (_, _) -> error "unionWithPlus got different dimensions" + return (HM.insert commonDExpr xysens hm) -lookupSymSet :: (Eq a, Hashable a, Unbox a) => Key -> StateT (FunGraph a b c) Identity (HS.HashSet (GExpr a)) -lookupSymSet k = do + +lookupSymSet :: (Eq a, Hashable a, Element a, DvdaDim sh) => + Expr sh a -> State (FunGraph a b c) (Maybe (HS.HashSet (DynamicExpr a))) +lookupSymSet expr = do fg <- get - let (_,symSet) = fromJust $ fgReverseLookup k fg - return symSet + case fgLookup expr fg of Just (_,symSet) -> return (Just symSet) + Nothing -> return Nothing +getSensitivities :: (Eq a, Floating a, Num (Vector a), Hashable a, LA.Container Vector a, DvdaDim sh) => + HS.HashSet (DynamicExpr a) -> Expr sh a -> Expr sh a + -> State (FunGraph a b c) (HM.HashMap (DynamicExpr a) (DynamicExpr a)) +getSensitivities _ (EGrad _ _) _ = error "don't call getSensitivities on EGrad" +getSensitivities _ (EJacob _ _) _ = error "don't call getSensitivities on EJacob" +getSensitivities _ (EDeriv _ _) _ = error "don't call getSensitivities on EDeriv" +getSensitivities _ (EScale _ _) _ = error "cant' do getSensitivities on EScale yet (needs EinSum?)" +getSensitivities _ (EDimensionless _) _ = return HM.empty +getSensitivities _ (EConst _) _ = return HM.empty +getSensitivities args (ERef sh k) sens = do + fg <- get + let expr = fromJust $ fgExprFromKey sh k fg + getSensitivities args expr sens +getSensitivities args primal@(ESym sh (SymDependent name k dep')) sens = do + let dprimal = makeDynamic primal + primalMap = + if HS.member dprimal args + then HM.fromList [(dprimal, makeDynamic sens)] + -- don't backprop if there aren't any interesting symbols farther in the tree + else HM.empty -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 + dep = ESym sh dep' + + depSymSet <- liftM fromJust $ lookupSymSet dep + + let commonSyms = HS.intersection args depSymSet + + dependentMap <- case HS.size commonSyms of + 0 -> return HM.empty + _ -> getSensitivities commonSyms dep (sens*primal') + where + primal' = ESym sh (SymDependent name (k+1) dep') + + return $ HM.union primalMap dependentMap +getSensitivities args primal@(ESym _ _) sens = do + let dprimal = makeDynamic primal + if HS.member dprimal args + then return $ HM.fromList [(dprimal, makeDynamic sens)] + -- 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 + +getSensitivities args (EUnary op g) sens = do + symSetG <- liftM fromJust $ lookupSymSet g 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 + let dfdg = dualPerturbation $ applyUnary op (Dual g 1) + getSensitivities args g (sens*dfdg) +getSensitivities args (EBinary op g h) sens = do + symSetG <- lookupSymSet g + symSetH <- lookupSymSet h - 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) + let 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) + gsens <- case liftM HS.size (liftM (HS.intersection args) symSetG) of + Nothing -> return HM.empty + Just 0 -> return HM.empty + _ -> getSensitivities args g (sens*dfdg) + hsens <- case liftM HS.size (liftM (HS.intersection args) symSetH) of + Nothing -> return HM.empty + Just 0 -> return HM.empty + _ -> getSensitivities args h (sens*dfdh) unionWithPlus gsens hsens +--getSensitivities args (EScale g h) sens = do +-- symSetG <- lookupSymSet g +-- symSetH <- lookupSymSet h +-- +-- fg <- get +-- let dfdg = h +-- dfdh = g +-- +-- gsens <- case liftM HS.size (liftM (HS.intersection args) symSetG) of +-- Nothing -> return HM.empty +-- Just 0 -> return HM.empty +-- _ -> getSensitivities args g (sens*dfdg) +-- hsens <- case liftM HS.size (liftM (HS.intersection args) symSetH) of +-- Nothing -> return HM.empty +-- 0 -> return HM.empty +-- _ -> getSensitivities args h (sens*dfdh) +-- unionWithPlus gsens hsens +--getSensitivities _ (EDeriv _ _) _ = error "don't call getSensitivities on EDeriv" +--getSensitivities _ (EGrad _ _) _ = error "don't call getSensitivities on EGrad" +--getSensitivities _ (EJacob _ _) _ = error "don't call getSensitivities on EJacob" - ---------------------- heterogenous inputs/outputs ------------------ data a :* b = a :* b deriving Show infixr 6 :* -class HList a where + +---------------------------------- input/output class --------------------------------------------- +class MkFunGraph 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 + type GenT a + mkNodes :: a -> State (FunGraph (NumT a) b c) a -instance (HList a, HList b, NumT a ~ NumT b) => HList (a :* b) where +instance (Hashable a, Eq a, Floating a, Num (Vector a), LA.Container Vector a) => + MkFunGraph (Expr DIM0 a) where + type NumT (Expr DIM0 a) = a + type GenT (Expr DIM0 a) = a + mkNodes = node + +instance (Hashable a, Eq a, Floating a, Num (Vector a), LA.Container Vector a) => + MkFunGraph (Expr DIM1 a) where + type NumT (Expr DIM1 a) = a + type GenT (Expr DIM1 a) = Vector a + mkNodes = node + +instance (Hashable a, Eq a, Floating a, Num (Vector a), LA.Container Vector a) => + MkFunGraph (Expr DIM2 a) where + type NumT (Expr DIM2 a) = a + type GenT (Expr DIM2 a) = Matrix a + mkNodes = node + +instance (Hashable a, Eq a, Floating a, Num (Vector a), LA.Container Vector a, MkFunGraph (Expr sh a), DvdaDim sh) => + MkFunGraph [Expr sh a] where + type NumT [Expr sh a] = a + type GenT [Expr sh a] = [GenT (Expr sh a)] + mkNodes = mapM node + +instance (Hashable a, Eq a, Floating a, Num (Vector a), LA.Container Vector a, MkFunGraph (Expr sh a), DvdaDim sh) => + MkFunGraph [[Expr sh a]] where + type NumT [[Expr sh a]] = a + type GenT [[Expr sh a]] = [[GenT (Expr sh a)]] + mkNodes = mapM (mapM node) + +--instance (Show a, MkFunGraph a) => MkFunGraph [a] where +-- type NumT [a] = NumT a +-- type GenT [a] = [GenT a] +-- type KeyT [a] = [KeyT a] +-- mkNodes xs = do +-- (x',kxs) <- mapM mkNodes xs >>= (return . unzip) +-- return (x', concat kxs) + +instance (MkFunGraph a, MkFunGraph b, NumT a ~ NumT b) => MkFunGraph (a :* b) where type NumT (a :* b) = NumT a - type DimT (a :* b) = DimT a :* DimT b + type GenT (a :* b) = GenT a :* GenT 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 + x' <- mkNodes x + y' <- mkNodes y + return (x' :* 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 +inputs :: MkFunGraph b => b -> State (FunGraph (NumT b) b c) b +inputs exprs_ = do + exprs <- mkNodes exprs_ FunGraph hm im _ outs <- get - put (FunGraph hm im (getHDim exprs, keys) outs) - return exprs' + put $ FunGraph hm im exprs outs + return exprs -outputs :: HList c => c -> StateT (FunGraph (NumT c) b (DimT c)) Identity c -outputs exprs = do - (exprs',keys) <- mkNodes exprs +outputs :: MkFunGraph c => c -> State (FunGraph (NumT c) b c) c +outputs exprs_ = do + exprs <- mkNodes exprs_ FunGraph hm im ins _ <- get - put (FunGraph hm im ins (getHDim exprs,keys)) - return exprs' + put $ FunGraph hm im ins exprs + return exprs -inputs_ :: HList b => b -> StateT (FunGraph (NumT b) (DimT b) c) Identity () +inputs_ :: MkFunGraph b => b -> State (FunGraph (NumT b) b c) () inputs_ exprs = do _ <- inputs exprs return () -outputs_ :: HList c => c -> StateT (FunGraph (NumT c) b (DimT c)) Identity () +outputs_ :: MkFunGraph c => c -> State (FunGraph (NumT c) b c) () 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 +------------------ utility function ----------------- +runFunGraph :: State (FunGraph a b c) 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 :: (MkFunGraph b, MkFunGraph c, NumT b ~ NumT c) => + b -> c -> FunGraph (NumT b) b c makeFunGraph ins outs = runFunGraph $ do inputs_ ins outputs_ outs + +-- | Show an Expr, looking up all ERefs +fullShow :: (Show a, Element a, DvdaDim sh) => FunGraph a b c -> Expr sh a -> String +fullShow fg = show . (recover fg) + +fullShowNodes :: (Show a, Element a) => FunGraph a b c -> String +fullShowNodes fg@(FunGraph _ im _ _) = + init $ unlines $ map (\(a,b) -> show a ++ ": " ++ (fullShow fg) (fromDynamic Z b)) (IM.toList im) + +-- | Take a FunGraph and an expression and traverse the expression. +-- . +-- Each time an ERef is found, look it up in the FunGraph and continue traversal +recover :: DvdaDim sh => FunGraph a b c -> Expr sh a -> Expr sh a +recover fg (ERef sh k) = recover fg (fromJust $ fgExprFromKey sh k fg) +recover _ e@(EDimensionless _) = e +recover _ e@(ESym _ _) = e +recover _ e@(EConst _) = e +recover fg (EUnary op x) = EUnary op (recover fg x) +recover fg (EBinary op x y) = EBinary op (recover fg x) (recover fg y) +recover fg (EDeriv x y) = EDeriv (recover fg x) (recover fg y) +recover fg (EGrad x y) = EGrad (recover fg x) (recover fg y) +recover fg (EJacob x y) = EJacob (recover fg x) (recover fg y) +recover fg (EScale x y) = EScale (recover fg x) (recover fg y) + +-- | "Pure" gradient which which runs rad and then calls recover to substitute values for ERefs +runDeriv :: (Eq a, Floating a, Num (Vector a), Hashable a, LA.Container Vector a, DvdaDim sh) + => Expr sh a -> [Expr sh a] -> [Expr sh a] +runDeriv expr args = map (recover fg) deda + where + (deda, fg) = runState (rad expr args) emptyFunGraph diff --git a/dvda.cabal b/dvda.cabal --- a/dvda.cabal +++ b/dvda.cabal @@ -1,14 +1,14 @@ Name: dvda -Version: 0.1.1 +Version: 0.2.0 License: BSD3 License-file: LICENSE Author: Greg Horn Maintainer: gregmainland@gmail.edu Stability: Experimental Category: Numerical, Math -Build-type: Simple +Build-type: Custom Synopsis: Efficient automatic differentiation -Cabal-version: >= 1.6 +Cabal-version: >= 1.8 Description: { dvda == DVDA Verifiably Differentiates Algorithmically . @@ -46,16 +46,22 @@ Library Exposed-modules: Dvda Dvda.BinUn + Dvda.CallNative + Dvda.Codegen Dvda.Config - Dvda.Dot +-- Dvda.Dot Dvda.Dual Dvda.Examples Dvda.Expr - Dvda.GExpr Dvda.Graph - Dvda.HSBuilder - Dvda.HSSyntax - Dvda.HomoDim +-- Dvda.HSBuilder +-- Dvda.HSSyntax + Dvda.MultipleShooting.MSCoctave + Dvda.MultipleShooting.MSMonad + Dvda.MultipleShooting.MultipleShooting + Dvda.MultipleShooting.Types + Dvda.OctaveSyntax + Dvda.SparseLA Dvda.SymMonad -- Dvda.CFunction -- Dvda.Codegen.CBuilder @@ -67,27 +73,42 @@ Build-depends: base >= 4 && < 5, hashable >= 1.1 && < 1.2, - vector >= 0.9 && < 0.10, - repa >= 3.1 && < 3.2, + repa >= 3.2 && < 3.3, 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 +-- process >= 1.1 && < 1.2, +-- text >= 0.11 && < 0.12, +-- plugins >= 1.5 && < 1.6, +-- deepseq >= 1.3 && < 1.4, + hmatrix >= 0.14 && < 0.15 + +-- latc >= 0.1 && < 0.2 -- unix -- text, --- QuickCheck, Ghc-options: -Wall -- Ghc-options: -O2 -Wall -threaded GHC-Prof-Options: -prof -fprof-auto + +flag test + description: Build test program. + default: False + +Test-suite test + type: exitcode-stdio-1.0 + hs-source-dirs: test, . + main-is: Test.hs + build-depends: base, + QuickCheck == 2.4.*, + ad, + test-framework-quickcheck2, + test-framework + ghc-options: -Wall -- Executable stressTest -- if flag(stressTest) diff --git a/test/Test.hs b/test/Test.hs new file mode 100644 --- /dev/null +++ b/test/Test.hs @@ -0,0 +1,146 @@ +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# OPTIONS_GHC -Wall #-} + +import Test.Framework (defaultMain, testGroup, Test(..)) +import Test.Framework.Providers.QuickCheck2 (testProperty) + +import Test.QuickCheck + +import Data.Array.Repa(DIM0, DIM1, DIM2, Z(..), Shape, shapeOfList) +import Data.Hashable ( Hashable ) + +import Numeric.LinearAlgebra ( Matrix, Vector, Element, fromList, fromLists, Container, (><) ) +import qualified Numeric.LinearAlgebra as LA +import Foreign.Storable (Storable) + +import Control.Monad + +import Dvda.BinUn +import Dvda.Expr + +import Dvda.CallNative +import Dvda.Dual + +import qualified Numeric.AD as AD + +-- Arbitrary math operations + +instance Arbitrary UnOp where + arbitrary = oneof $ map return [minBound..maxBound] + +instance Arbitrary BinOp where + arbitrary = oneof $ map return [minBound..maxBound] + +-- Arbitrary constants + +class Shape sh => ArbSingleton sh where + arbCSingleton :: Arbitrary a => Gen (Const sh a) + +instance ArbSingleton Z where + arbCSingleton = (liftM (CSingleton Z)) $ arbitrary + +cvec xs = CVec (shapeOfList [length xs]) (fromList xs) + +arbCVector :: (Storable a, Arbitrary a) => Gen (Const DIM1 a) +arbCVector = liftM cvec $ listOf1 $ arbitrary + +cmat (r,c) xs + | r*c == sum (map length xs) && r == length xs = CMat (shapeOfList [c,r]) (fromLists xs) + | otherwise = error $ "bad dims in mat!"++ + "\ngiven (r,c): " ++ show (r,c) ++ + "\nactual (r,c): " ++ show (length xs, map length xs) + +arbCMatrix :: (Element a, Arbitrary a) => Gen (Const DIM2 a) +arbCMatrix = do r <- choose (0, 100) + c <- choose (0, 100) + l <- vectorOf (r*c) arbitrary + return $ cmat (r, c) l + +instance Arbitrary a => Arbitrary (Const Z a) where + arbitrary = arbCSingleton + +instance (Arbitrary a, Storable a) => Arbitrary (Const DIM1 a) where + arbitrary = arbCVector + +instance (Arbitrary a, Element a) => Arbitrary (Const DIM2 a) where + arbitrary = arbCMatrix + +-- Arbitrary expressions + +class Shape sh => ArbExpr sh where + arbExpr :: Arbitrary a => Gen (Expr sh a) + +arbConst :: (Shape sh, Arbitrary (Const sh a), Arbitrary a) => Gen (Expr sh a) +arbConst = liftM EConst arbitrary + +arbUnary :: (Shape sh, Arbitrary (Expr sh a)) => Gen (Expr sh a) +arbUnary = liftM2 EUnary arbitrary arbitrary + +arbBinary :: (Shape sh, Arbitrary (Expr sh a)) => Gen (Expr sh a) +arbBinary = liftM3 EBinary arbitrary arbitrary arbitrary + +instance ArbExpr Z where + arbExpr = oneof [arbConst, arbUnary, arbBinary] + +instance Arbitrary a => Arbitrary (Expr Z a) where + arbitrary = arbExpr + +-- Arbitrary numerical functions + +binary :: Floating a => [a -> a -> a] +binary = [(*), (+), (-), (/)] + +unary :: Floating a => [a -> a] +unary = [abs, negate, signum, exp, sqrt, log, sin, cos, tan, asin, acos, atan, tanh, sinh, cosh, atanh, asinh, acosh] + +-- We have to do this by hand because of all kinds of stupid type +-- shit. Otherwise we'd map over unary. + +cosProp, sinProp, tanProp, coshProp, sinhProp, tanhProp, expProp, signumProp, sqrtProp, negateProp, absProp :: (Eq a, Floating a, Num (Vector a), Show a, Hashable a, Container Vector a) => a -> Bool +cosProp x = cos x == nativeRun cos x +sinProp x = sin x == nativeRun sin x +tanProp x = tan x == nativeRun tan x +coshProp x = cosh x == nativeRun cosh x +sinhProp x = sinh x == nativeRun sinh x +tanhProp x = tanh x == nativeRun tanh x +expProp x = exp x == nativeRun exp x +signumProp x = signum x == nativeRun signum x +sqrtProp x = sqrt x == nativeRun sqrt x +negateProp x = negate x == nativeRun negate x +absProp x = abs x == nativeRun abs x + +props :: [Double -> Bool] +props = [ cosProp, sinProp, tanProp + , coshProp, sinhProp, tanhProp + , expProp, signumProp, sqrtProp + , negateProp, absProp + ] + +propNames :: [String] +propNames = [ "cos", "sin", "tan", "cosh", "sinh", "tanh", "exp", "signum", "sqrt", "negate", "abs"] + +acosProp, asinProp, atanProp :: (Floating a, Num (Vector a), Ord a, Show a, Hashable a, Container Vector a) => a -> Property + +acosProp x = x >= 0 && x <= 1 ==> acos x == nativeRun acos x +asinProp x = x >= 0 && x <= 1 ==> asin x == nativeRun asin x +atanProp x = x >= 0 && x <= 1 ==> atan x == nativeRun atan x + +prop'Names :: [String] +prop'Names = ["acos", "asin", "atan"] + +props' :: [Double -> Property] +props' = [ acosProp, asinProp, atanProp ] + +mkUnaryTest (n, t) = testProperty ("unary_" ++ n) t + +tests, uts :: [Test] +uts = (map mkUnaryTest (zip propNames props)) ++ (map mkUnaryTest (zip prop'Names props')) + +tests = [ + testGroup "Unary functions 1" uts + ] + +main :: IO () +main = defaultMain tests +