dvda 0.1.1 → 0.2.0
raw patch · 23 files changed
+2457/−943 lines, 23 filesdep +QuickCheckdep +addep +hmatrixdep −deepseqdep −pluginsdep −processdep ~basedep ~repabuild-type:Custom
Dependencies added: QuickCheck, ad, hmatrix, test-framework, test-framework-quickcheck2
Dependencies removed: deepseq, plugins, process, text, transformers, vector
Dependency ranges changed: base, repa
Files
- Dvda.hs +40/−3
- Dvda/BinUn.hs +54/−6
- Dvda/CallNative.hs +202/−0
- Dvda/Codegen.hs +38/−0
- Dvda/Config.hs +16/−7
- Dvda/Dot.hs +0/−67
- Dvda/Dual.hs +10/−0
- Dvda/Examples.hs +106/−50
- Dvda/Expr.hs +244/−116
- Dvda/GExpr.hs +0/−87
- Dvda/Graph.hs +152/−49
- Dvda/HSBuilder.hs +0/−116
- Dvda/HSSyntax.hs +0/−148
- Dvda/HomoDim.hs +0/−36
- Dvda/MultipleShooting/MSCoctave.hs +273/−0
- Dvda/MultipleShooting/MSMonad.hs +215/−0
- Dvda/MultipleShooting/MultipleShooting.hs +166/−0
- Dvda/MultipleShooting/Types.hs +88/−0
- Dvda/OctaveSyntax.hs +165/−0
- Dvda/SparseLA.hs +245/−0
- Dvda/SymMonad.hs +260/−242
- dvda.cabal +37/−16
- test/Test.hs +146/−0
Dvda.hs view
@@ -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
Dvda/BinUn.hs view
@@ -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 ++" )"
+ Dvda/CallNative.hs view
@@ -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
+ Dvda/Codegen.hs view
@@ -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
Dvda/Config.hs view
@@ -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+
− Dvda/Dot.hs
@@ -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-
Dvda/Dual.hs view
@@ -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)
Dvda/Examples.hs view
@@ -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
Dvda/Expr.hs view
@@ -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
− Dvda/GExpr.hs
@@ -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 _ _) = []
Dvda/Graph.hs view
@@ -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"
− Dvda/HSBuilder.hs
@@ -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)
− Dvda/HSSyntax.hs
@@ -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) )"
− Dvda/HomoDim.hs
@@ -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"
+ Dvda/MultipleShooting/MSCoctave.hs view
@@ -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"
+ Dvda/MultipleShooting/MSMonad.hs view
@@ -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+ }
+ Dvda/MultipleShooting/MultipleShooting.hs view
@@ -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)
+ Dvda/MultipleShooting/Types.hs view
@@ -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
+ Dvda/OctaveSyntax.hs view
@@ -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
+ Dvda/SparseLA.hs view
@@ -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
Dvda/SymMonad.hs view
@@ -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
dvda.cabal view
@@ -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)
+ test/Test.hs view
@@ -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+