dvda 0.2.2 → 0.3
raw patch · 28 files changed
+2038/−2539 lines, 28 filesdep +dvdadep +file-locationdep +hashtablesdep −hmatrixdep −repa
Dependencies added: dvda, file-location, hashtables, process
Dependencies removed: hmatrix, repa
Files
- Dvda.hs +6/−57
- Dvda/AD.hs +74/−0
- Dvda/BinUn.hs +0/−191
- Dvda/CGen.hs +295/−0
- Dvda/CSE.hs +153/−0
- Dvda/CallNative.hs +0/−202
- Dvda/Codegen.hs +0/−38
- Dvda/Codegen/Gcc.hs +32/−0
- Dvda/Codegen/WriteFile.hs +32/−0
- Dvda/Config.hs +0/−125
- Dvda/Examples.hs +41/−138
- Dvda/Expr.hs +562/−326
- Dvda/FunGraph.hs +152/−0
- Dvda/Graph.hs +0/−217
- Dvda/MultipleShooting/CoctaveTemplates.hs +122/−0
- Dvda/MultipleShooting/MSCoctave.hs +217/−207
- Dvda/MultipleShooting/MSMonad.hs +61/−121
- Dvda/MultipleShooting/MultipleShooting.hs +0/−166
- Dvda/MultipleShooting/Types.hs +35/−33
- Dvda/OctaveSyntax.hs +0/−165
- Dvda/Reify.hs +88/−0
- Dvda/ReifyGraph.hs +16/−0
- Dvda/SymMonad.hs +0/−342
- Dvda/Vis.hs +69/−0
- LICENSE +2/−1
- TestMain.hs +18/−0
- dvda.cabal +63/−64
- test/Test.hs +0/−146
Dvda.hs view
@@ -6,43 +6,18 @@ -} {-# OPTIONS_GHC -Wall #-}-{-# Language TypeOperators #-}-{-# Language TypeFamilies #-}-{-# Language MultiParamTypeClasses #-}-{-# Language FlexibleInstances #-} module Dvda ( -- * primitives sym , symDependent , symDependentN- , vsym- , msym- , svec- , smat- , vec- , mat -- * operations- , scale--- , dot- , diff- , runDeriv+ , rad -- * symbolic expression type , Expr- , fullShow- , fullShowNodes -- * construct FunGraphs- , FunGraph- , makeFunGraph- , runFunGraph- , inputs- , outputs- , inputs_- , outputs_- , node+ , toFunGraph -- * show/summarize FunGraphs- , funGraphSummary- , funGraphSummary'- , showCollisions , previewGraph -- * compile and link function -- , buildHSFunction@@ -50,36 +25,10 @@ -- , 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.AD ( rad )+import Dvda.Expr ( Expr, sym, symDependent, symDependentN )+import Dvda.FunGraph ( toFunGraph, (:*)(..) )+import Dvda.Vis ( previewGraph ) --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/AD.hs view
@@ -0,0 +1,74 @@+{-# OPTIONS_GHC -Wall #-}++module Dvda.AD ( backprop+ , rad+ ) where++import Data.Hashable ( Hashable )++import Dvda.Dual hiding ( fad, fad' )+import Dvda.Expr+import Dvda.HashMap ( HashMap )+import qualified Dvda.HashMap as HM++--fad :: Num a => (Dual a -> [Dual a]) -> a -> [a]+--fad f x = map dualPerturbation $ f (Dual x 1)++bpBinary :: (Eq a, Num a)+ => Expr a -> Expr a -> Expr a+ -> (Dual (Expr a) -> Dual (Expr a) -> Dual (Expr a))+ -> [(Expr a, Expr a)]+bpBinary sens g h binop = gsens ++ hsens+ where+ dfdg = dualPerturbation $ binop (Dual g 1) (Dual h 0)+ dfdh = dualPerturbation $ binop (Dual g 0) (Dual h 1)+ gsens = backpropNode (sens*dfdg) g+ hsens = backpropNode (sens*dfdh) h++bpUnary :: (Eq a, Num a)+ => Expr a -> Expr a+ -> (Dual (Expr a) -> Dual (Expr a))+ -> [(Expr a, Expr a)]+bpUnary sens g unop = backpropNode (sens*dfdg) g+ where+ dfdg = dualPerturbation $ unop (Dual g 1)++backpropNode :: (Eq a, Num a) => Expr a -> Expr a -> [(Expr a, Expr a)]+backpropNode sens e@(ESym (SymDependent name k dep_)) = (e,sens):(backpropNode (sens*primal') dep)+ where+ primal' = ESym (SymDependent name (k+1) dep_)+ dep = ESym dep_+backpropNode sens e@(ESym (Sym _)) = [(e,sens)]+backpropNode _ (EConst _) = []+backpropNode _ (ENum (FromInteger _)) = []+backpropNode _ (EFractional (FromRational _)) = []+backpropNode sens (ENum (Mul x y)) = bpBinary sens x y (*)+backpropNode sens (ENum (Add x y)) = bpBinary sens x y (+)+backpropNode sens (ENum (Sub x y)) = bpBinary sens x y (-)+backpropNode sens (ENum (Abs x)) = bpUnary sens x abs+backpropNode sens (ENum (Negate x)) = bpUnary sens x negate+backpropNode sens (ENum (Signum x)) = bpUnary sens x signum+backpropNode sens (EFractional (Div x y)) = bpBinary sens x y (/)+backpropNode sens (EFloating (Pow x y)) = bpBinary sens x y (**)+backpropNode sens (EFloating (LogBase x y)) = bpBinary sens x y logBase+backpropNode sens (EFloating (Exp x)) = bpUnary sens x exp+backpropNode sens (EFloating (Log x)) = bpUnary sens x log+backpropNode sens (EFloating (Sin x)) = bpUnary sens x sin+backpropNode sens (EFloating (Cos x)) = bpUnary sens x cos+backpropNode sens (EFloating (ASin x)) = bpUnary sens x asin+backpropNode sens (EFloating (ATan x)) = bpUnary sens x atan+backpropNode sens (EFloating (ACos x)) = bpUnary sens x acos+backpropNode sens (EFloating (Sinh x)) = bpUnary sens x sinh+backpropNode sens (EFloating (Cosh x)) = bpUnary sens x cosh+backpropNode sens (EFloating (Tanh x)) = bpUnary sens x tanh+backpropNode sens (EFloating (ASinh x)) = bpUnary sens x asinh+backpropNode sens (EFloating (ATanh x)) = bpUnary sens x atanh+backpropNode sens (EFloating (ACosh x)) = bpUnary sens x acosh++backprop :: (Num a, Eq a, Hashable a) => Expr a -> HashMap (Expr a) (Expr a)+backprop x = HM.fromListWith (+) (backpropNode 1 x)++rad :: (Num a, Eq a, Hashable a) => Expr a -> [Expr a] -> [Expr a]+rad x args = map (\arg -> HM.lookupDefault 0 arg sensitivities) args+ where+ sensitivities = backprop x
− Dvda/BinUn.hs
@@ -1,191 +0,0 @@-{-# OPTIONS_GHC -Wall #-}--module Dvda.BinUn ( BinOp(..)- , UnOp(..)- , showBinary- , showUnary- , applyUnary- , applyBinary- , unaryDeriv- , binaryDeriv- , isCommutative- , lassoc- , rassoc- ) where--import Data.Hashable ( Hashable, hash )--import Dvda.Dual ( Dual(..), dualPerturbation )--data UnOp = Abs- | Neg- | Signum- | Exp- | Sqrt- | Log- | Sin- | Cos- | Tan- | ASin- | ACos- | ATan- | Tanh- | Sinh- | Cosh- | ATanh- | ASinh- | ACosh deriving (Eq, Show, Enum, Bounded)--data BinOp = Add- | Sub- | Mul- | Div- | Pow- | LogBase deriving (Eq, Show, Enum, Bounded)--instance Hashable UnOp where- hash Abs = 0- hash Neg = 1- hash Signum = 2- hash Exp = 3- hash Sqrt = 4- hash Log = 5- hash Sin = 6- hash Cos = 7- hash Tan = 8- hash ASin = 9- hash ACos = 10- hash ATan = 11- hash Tanh = 12- hash Sinh = 13- hash Cosh = 14- hash ATanh = 15- hash ASinh = 16- hash ACosh = 17--instance Hashable BinOp where- hash Add = 18- hash Sub = 19- hash Mul = 20- hash Div = 21- hash Pow = 22- hash LogBase = 23- -showUnary :: String -> UnOp -> String-showUnary x Abs = '|': x ++ "|"-showUnary x Neg = '-':paren x-showUnary x Signum = "signum"++paren x-showUnary x Exp = "exp"++paren x-showUnary x Sqrt = "sqrt"++paren x-showUnary x Log = "log"++paren x-showUnary x Sin = "sin"++paren x-showUnary x Cos = "cos"++paren x-showUnary x Tan = "tan"++paren x-showUnary x ASin = "asin"++paren x-showUnary x ACos = "acos"++paren x-showUnary x ATan = "atan"++paren x-showUnary x Sinh = "sinh"++paren x-showUnary x Cosh = "cosh"++paren x-showUnary x Tanh = "tanh"++paren x-showUnary x ASinh = "asinh"++paren x-showUnary x ATanh = "atanh"++paren x-showUnary x ACosh = "acosh"++paren x--applyUnary :: Floating a => UnOp -> a -> a-applyUnary Abs = abs-applyUnary Neg = negate-applyUnary Signum = signum-applyUnary Exp = exp-applyUnary Sqrt = sqrt-applyUnary Log = log-applyUnary Sin = sin-applyUnary Cos = cos-applyUnary Tan = tan-applyUnary ASin = asin-applyUnary ACos = acos-applyUnary ATan = atan-applyUnary Sinh = sinh-applyUnary Cosh = cosh-applyUnary Tanh = tanh-applyUnary ASinh = asinh-applyUnary ATanh = atanh-applyUnary ACosh = acosh--applyBinary :: Floating a => BinOp -> a -> a -> a-applyBinary Add = (+)-applyBinary Sub = (-)-applyBinary Mul = (*)-applyBinary Div = (/)-applyBinary Pow = (**)-applyBinary LogBase = logBase--unaryDeriv :: Floating a => UnOp -> (a,a) -> a-unaryDeriv op (x,x') = dualPerturbation $ applyUnary op (Dual x x')--binaryDeriv :: Floating a => BinOp -> (a,a) -> (a,a) -> a-binaryDeriv op (x,x') (y,y') = dualPerturbation $ applyBinary op (Dual x x') (Dual y y')--showBinary :: BinOp -> String-showBinary Add = "+"-showBinary Sub = "-"-showBinary Mul = "*"-showBinary Div = "/"-showBinary Pow = "**"-showBinary LogBase = "`logbase`"--isCommutative :: BinOp -> Bool-isCommutative Add = True-isCommutative Sub = False-isCommutative Mul = True-isCommutative Div = False-isCommutative Pow = False-isCommutative LogBase = False--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/CGen.hs view
@@ -0,0 +1,295 @@+{-# OPTIONS_GHC -Wall #-}+{-# Language TemplateHaskell #-}+{-# Language TypeFamilies #-}+{-# Language FlexibleContexts #-}++module Dvda.CGen ( showC+ , showMex+ , MatrixStorageOrder(..)+ ) where+++import Data.Hashable ( Hashable )+import Data.List ( intercalate )+import FileLocation ( err )+import Text.Printf ( printf )++import Dvda.Expr ( GExpr(..), Floatings(..), Nums(..), Fractionals(..) )+import Dvda.FunGraph ( FunGraph, MVS(..), topSort, fgInputs, fgOutputs, fgLookupGExpr )+import Dvda.HashMap ( HashMap )+import qualified Dvda.HashMap as HM++data MatrixStorageOrder = RowMajor | ColMajor++-- | take a list of pair of inputs to indices which reference them+-- create a hashmap from GSyms to strings which hold the declaration+makeInputMap :: (Eq a, Hashable a, Show a)+ => MatrixStorageOrder -> [MVS (GExpr a Int)] -> HashMap (GExpr a Int) String+makeInputMap matStorageOrder ins = HM.fromList $ concat $ zipWith writeInput [(0::Int)..] ins+ where+ writeInput inputK (Sca g) = [(g, printf "*input%d; /* %s */" inputK (show g))]+ writeInput inputK (Vec gs) = zipWith f [(0::Int)..] gs+ where+ f inIdx g = (g, printf "input%d[%d]; /* %s */" inputK inIdx (show g))+ writeInput inputK (Mat gs)+ | any ((ncols /=) . length) gs =+ error $ "writeInputs [[GraphRef]] matrix got inconsistent column dimensions: "++ show (map length gs)+ | otherwise = zipWith f [(r,c) | r <- [0..(nrows-1)], c <- [0..(ncols-1)]] (concat gs)+ where+ nrows = length gs+ ncols = if nrows == 0 then 0 else length (head gs)+ f (rowIdx,colIdx) g = (g,printf "input%d[%d][%d]; /* %s */" inputK fstIdx sndIdx (show g))+ where+ (fstIdx,sndIdx) = case matStorageOrder of RowMajor -> (rowIdx,colIdx)+ ColMajor -> (colIdx,rowIdx)++writeInputPrototypes :: MatrixStorageOrder -> [MVS a] -> [String]+writeInputPrototypes matStorageOrder ins = concat $ zipWith inputPrototype [(0::Int)..] ins+ where+ inputPrototype inputK (Sca _) = ["const double * input" ++ show inputK]+ inputPrototype inputK (Vec gs) = ["const double input" ++ show inputK ++ "[" ++ show (length gs) ++ "]"]+ inputPrototype inputK (Mat gs)+ | any ((ncols /=) . length) gs =+ error $ "writeInputs [[GraphRef]] matrix got inconsistent column dimensions: "++ show (map length gs)+ | otherwise = ["const double input" ++ show inputK ++ "[" ++ show fstIdx ++ "][" ++ show sndIdx ++ "]"]+ where+ nrows = length gs+ ncols = if nrows == 0 then 0 else length (head gs)+ (fstIdx,sndIdx) = case matStorageOrder of RowMajor -> (nrows,ncols)+ ColMajor -> (ncols,nrows)++writeOutputs :: MatrixStorageOrder -> [MVS Int] -> ([String], [String])+writeOutputs matStorageOrder ins = (concatMap fst dcs, concatMap snd dcs)+ where+ dcs :: [([String],[String])]+ dcs = zipWith writeOutput ins [0..]++ writeOutput :: MVS Int -> Int -> ([String], [String])+ writeOutput (Sca gref) outputK = (decls, prototype)+ where+ decls = [printf "/* output %d */" outputK, printf "(*output%d) = %s;" outputK (nameNode gref)]+ prototype = ["double * const output" ++ show outputK]+ writeOutput (Vec grefs) outputK = (decls, prototype)+ where+ prototype = ["double output" ++ show outputK ++ "[" ++ show (length grefs) ++ "]"]+ decls = (printf "/* output %d */" outputK):+ zipWith f [(0::Int)..] grefs+ where+ f outIdx gref = printf "output%d[%d] = %s;" outputK outIdx (nameNode gref)+ writeOutput (Mat grefs) outputK+ | any ((ncols /=) . length) grefs =+ error $ "writeOutputs [[GraphRef]] matrix got inconsistent column dimensions: "++ show (map length grefs)+ | otherwise = (decls, prototype)+ where+ nrows = length grefs+ ncols = if nrows == 0 then 0 else length (head grefs)+ prototype = ["double output" ++ show outputK ++ "[" ++ show fstIdx ++ "][" ++ show sndIdx ++ "]"]+ where+ (fstIdx,sndIdx) = case matStorageOrder of RowMajor -> (nrows,ncols)+ ColMajor -> (ncols,nrows)+ decls = (printf "/* output %d */" outputK):+ zipWith f [(r,c) | r <- [0..(nrows-1)], c <- [0..(ncols-1)]] (concat grefs)+ where+ f (rowIdx,colIdx) gref = printf "output%d[%d][%d] = %s;" outputK fstIdx sndIdx (nameNode gref)+ where+ (fstIdx,sndIdx) = case matStorageOrder of RowMajor -> (rowIdx,colIdx)+ ColMajor -> (colIdx,rowIdx)+++createMxOutputs :: [MVS Int] -> [String]+createMxOutputs xs = concat $ zipWith createMxOutput xs [0..]+ where+ createMxOutput :: MVS Int -> Int -> [String]+ createMxOutput (Sca _) outputK =+ [ " if ( " ++ show outputK ++ " < nlhs ) {"+ , " plhs[" ++ show outputK ++ "] = mxCreateDoubleScalar( 0 );"+ , " outputs[" ++ show outputK ++ "] = mxGetPr( plhs[" ++ show outputK ++ "] );"+ , " } else"+ , " outputs[" ++ show outputK ++ "] = (double*)malloc( sizeof(double) );"+ ]+ createMxOutput (Vec grefs) outputK =+ [ " if ( " ++ show outputK ++ " < nlhs ) {"+ , " plhs[" ++ show outputK ++ "] = mxCreateDoubleMatrix( " ++ show (length grefs) ++ ", 1, mxREAL );"+ , " outputs[" ++ show outputK ++ "] = mxGetPr( plhs[" ++ show outputK ++ "] );"+ , " } else"+ , " outputs[" ++ show outputK ++ "] = (double*)malloc( " ++ show (length grefs) ++ "*sizeof(double) );"+ ]+ createMxOutput (Mat grefs) outputK =+ [ " if ( " ++ show outputK ++ " < nlhs ) {"+ , " plhs[" ++ show outputK ++ "] = mxCreateDoubleMatrix( " ++ show nrows++ ", " ++ show ncols ++ ", mxREAL );"+ , " outputs[" ++ show outputK ++ "] = mxGetPr( plhs[" ++ show outputK ++ "] );"+ , " } else"+ , " outputs[" ++ show outputK ++ "] = (double*)malloc( " ++ show (nrows*ncols) ++ "*sizeof(double) );"+ ]+ where+ nrows = length grefs+ ncols = if nrows == 0 then 0 else length (head grefs)+++checkMxInputDims :: MVS a -> String -> Int -> [String]+checkMxInputDims (Sca _) functionName inputK =+ [ " if ( 1 != mxGetM( prhs[" ++ show inputK ++ "] ) || 1 != mxGetN( prhs[" ++ show inputK ++ "] ) ) {"+ , " char errMsg[200];"+ , " sprintf(errMsg,"+ , " \"mex function '" ++ functionName ++ "' got incorrect dimensions for input " ++ show (1+inputK) ++ "\\n\""+ , " \"expected dimensions: (1, 1) but got (%zu, %zu)\","+ , " mxGetM( prhs[" ++ show inputK ++ "] ),"+ , " mxGetN( prhs[" ++ show inputK ++ "] ) );"+ , " mexErrMsgTxt(errMsg);"+ , " }"+ ]+checkMxInputDims (Vec grefs) functionName inputK =+ [ " if ( !( " ++ show nrows ++ " == mxGetM( prhs[" ++ show inputK ++ "] ) && 1 == mxGetN( prhs[" ++ show inputK ++ "] ) ) && !( " ++ show nrows ++ " == mxGetN( prhs[" ++ show inputK ++ "] ) && 1 == mxGetM( prhs[" ++ show inputK ++ "] ) ) ) {"+ , " char errMsg[200];"+ , " sprintf(errMsg,"+ , " \"mex function '" ++ functionName ++ "' got incorrect dimensions for input " ++ show (1+inputK) ++ "\\n\""+ , " \"expected dimensions: (" ++ show nrows ++ ", 1) or (1, " ++ show nrows ++ ") but got (%zu, %zu)\","+ , " mxGetM( prhs[" ++ show inputK ++ "] ),"+ , " mxGetN( prhs[" ++ show inputK ++ "] ) );"+ , " mexErrMsgTxt(errMsg);"+ , " }"+ ]+ where+ nrows = length grefs+checkMxInputDims (Mat grefs) functionName inputK =+ [ " if ( " ++ show nrows ++ " != mxGetM( prhs[" ++ show inputK ++ "] ) || " ++ show ncols ++ " != mxGetN( prhs[" ++ show inputK ++ "] ) ) {"+ , " char errMsg[200];"+ , " sprintf(errMsg,"+ , " \"mex function '" ++ functionName ++ "' got incorrect dimensions for input " ++ show (1+inputK) ++ "\\n\""+ , " \"expected dimensions: (" ++ show nrows ++ ", " ++ show ncols ++ ") but got (%zu, %zu)\","+ , " mxGetM( prhs[" ++ show inputK ++ "] ),"+ , " mxGetN( prhs[" ++ show inputK ++ "] ) );"+ , " mexErrMsgTxt(errMsg);"+ , " }"+ ]+ where+ nrows = length grefs+ ncols = if nrows == 0 then 0 else length (head grefs)+++-- | Turns a FunGraph into a string containing C code+showC :: (Eq a, Show a, Hashable a) => MatrixStorageOrder -> String -> FunGraph a -> String+showC matStorageOrder functionName fg = txt+ where+ inPrototypes = writeInputPrototypes matStorageOrder (fgInputs fg)+ (outDecls, outPrototypes) = writeOutputs matStorageOrder (fgOutputs fg)+ inputMap = makeInputMap matStorageOrder (fgInputs fg)+ mainDecls = let f k = case fgLookupGExpr fg k of+ Just v -> cAssignment inputMap k v+ Nothing -> error $ "couldn't find node " ++ show k ++ " in fungraph :("+ in map f $ reverse $ topSort fg+ + body = unlines $ map (" "++) $+ mainDecls ++ [""] +++ outDecls+ + txt = "#include <math.h>\n\n" +++ "void " ++ functionName ++ " ( " ++ (intercalate ", " (inPrototypes++outPrototypes)) ++ " )\n{\n" +++ body ++ "}\n"++nameNode :: Int -> String+nameNode k = "v_" ++ show k++cAssignment :: (Eq a, Hashable a, Show a) => HashMap (GExpr a Int) String -> Int -> GExpr a Int -> String+cAssignment inputMap k g@(GSym _) = case HM.lookup g inputMap of+ Nothing -> error $ "cAssignment: couldn't find " ++ show g ++ " in the input map"+ Just str -> "const double " ++ nameNode k ++ " = " ++ str+cAssignment inputMap k gexpr = "const double " ++ nameNode k ++ " = " ++ toCOp gexpr ++ ";"+ where+ bin :: Int -> Int -> String -> String+ bin x y op = nameNode x ++ " " ++ op ++ " " ++ nameNode y+ + un :: Int -> String -> String+ un x op = op ++ "( " ++ nameNode x ++ " )"++ asTypeOfG :: a -> GExpr a b -> a+ asTypeOfG x _ = x+ + toCOp (GSym _) = $(err "This should be impossible")+ toCOp (GConst c) = show c+ toCOp (GNum (Mul x y)) = bin x y "*"+ toCOp (GNum (Add x y)) = bin x y "+"+ toCOp (GNum (Sub x y)) = bin x y "-"+ toCOp (GNum (Negate x)) = un x "-"+ toCOp (GNum (Abs x)) = un x "abs"+ toCOp (GNum (Signum x)) = un x "sign"+ toCOp (GNum (FromInteger x)) = show x+ toCOp (GFractional (Div x y)) = bin x y "/"+ toCOp (GFractional (FromRational x)) = show (fromRational x `asTypeOfG` gexpr)+ toCOp (GFloating (Pow x y)) = "pow( " ++ nameNode x ++ ", " ++ nameNode y ++ " )"+ toCOp (GFloating (LogBase x y)) = "log( " ++ nameNode y ++ ") / log( " ++ nameNode x ++ " )"+ toCOp (GFloating (Exp x)) = un x "exp"+ toCOp (GFloating (Log x)) = un x "log"+ toCOp (GFloating (Sin x)) = un x "sin"+ toCOp (GFloating (Cos x)) = un x "cos"+ toCOp (GFloating (ASin x)) = un x "asin"+ toCOp (GFloating (ATan x)) = un x "atan"+ toCOp (GFloating (ACos x)) = un x "acos"+ toCOp (GFloating (Sinh x)) = un x "sinh"+ toCOp (GFloating (Cosh x)) = un x "cosh"+ toCOp (GFloating (Tanh x)) = un x "tanh"+ toCOp (GFloating (ASinh _)) = error "C generation doesn't support ASinh"+ toCOp (GFloating (ATanh _)) = error "C generation doesn't support ATanh"+ toCOp (GFloating (ACosh _)) = error "C generation doesn't support ACosh"+++showMex :: (Eq a, Show a, Hashable a) => String -> FunGraph a -> String+showMex functionName fg = cText ++ "\n\n\n" ++ mexFun functionName (fgInputs fg) (fgOutputs fg)+ where+ cText = showC ColMajor functionName fg -- matlab is column major >_<++mexFun :: String -> [MVS a] -> [MVS Int] -> String+mexFun functionName ins outs =+ unlines $+ [ "#include \"mex.h\""+ , []+ , "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])"+ , "{"+ , " /* check number of inputs */"+ , " if ( " ++ show nrhs ++ " != nrhs ) {"+ , " char errMsg[200];"+ , " sprintf(errMsg,"+ , " \"mex function '" ++ functionName ++ "' given incorrect number of inputs\\n\""+ , " \"expected: " ++ show nrhs ++ " but got %d\","+ , " nrhs);"+ , " mexErrMsgTxt(errMsg);"+ , " }"+ , []+ , " /* check the dimensions of the input arrays */"+ ] ++ concat (zipWith (\x -> checkMxInputDims x functionName) ins [0..]) +++ [ []+ , " /* check number of outputs */"+ , " if ( " ++ show nlhs ++ " < nlhs ) {"+ , " char errMsg[200];"+ , " sprintf(errMsg,"+ , " \"mex function '" ++ functionName ++ "' saw too many outputs\\n\""+ , " \"expected <= " ++ show nlhs ++ " but got %d\","+ , " nlhs);"+ , " mexErrMsgTxt(errMsg);"+ , " }"+ , []+ , " /* create the output arrays, if no output is provided by user create a dummy output */"+ , " double * outputs[" ++ show nlhs ++ "];"+ ] ++ createMxOutputs outs ++ -- e.g.: plhs[0] = mxCreateDoubleMatrix(1,ncols,mxREAL);+ [ []+ , " /* call the c function */"+ , " " ++ functionName ++ "( " ++ intercalate ", " (inputPtrs ++ outputPtrs) ++ " );"+ , []+ , " /* free the unused dummy outputs */"+ , " int k;"+ , " for ( k = " ++ show (nlhs - 1) ++ "; nlhs <= k; k-- )"+ , " free( outputs[k] );"+ , "}"+ ]+ where+ nlhs = length outs+ nrhs = length ins+ inputPtrs = zipWith (\i k -> cast i "const " ++ "mxGetPr(prhs[" ++ show k ++ "])") ins [(0::Int)..]+ outputPtrs = zipWith (\o k -> cast o "" ++ "(outputs[" ++ show k ++ "])") outs [(0::Int)..]++ cast :: MVS a -> String -> String+ cast (Sca _) _ = ""+ cast (Vec _) _ = ""+ cast (Mat xs) cnst = "(" ++ cnst ++ "double (*)[" ++ show nrows ++ "])" -- column major order+ where+ nrows = length xs
+ Dvda/CSE.hs view
@@ -0,0 +1,153 @@+{-# OPTIONS_GHC -Wall #-}++module Dvda.CSE ( cse+ ) where++import Control.Monad.ST ( ST, runST )+import Data.Foldable ( toList )+import Data.Hashable ( Hashable )+import Data.IntMap ( IntMap )+import qualified Data.IntMap as IM+import Data.Tuple ( swap )++import Dvda.Expr ( GExpr(..), Floatings(..), Fractionals(..), Nums(..) )+import Dvda.FunGraph++import qualified Data.HashTable.Class as HT+import qualified Data.HashTable.ST.Cuckoo as C+type HashTable s v k = C.HashTable s v k++cse :: (Eq a, Hashable a) => FunGraph a -> FunGraph a+cse fg = nodelistToFunGraph (map swap htList) (fgInputs fg) outputIndices+ where+ (htList, im) = cse' (fgLookupGExpr fg) (fgOutputs fg)+ -- since the fgInputs are all symbolic (GSym _) there is no need for mapping old inputs to new inputs+ outputIndices = let+ oldIndexToNewIndex k = case IM.lookup k im of+ Just k' -> k'+ Nothing -> error $+ "CSE error, in mapping old output indices to new, found an old one which was missing from" +++ "the old --> new Int mapping"+ in map (fmap oldIndexToNewIndex) (fgOutputs fg)++cse' ::+ (Eq a, Hashable a)+ => (Int -> Maybe (GExpr a Int))+ -> [MVS Int]+ -> ([(GExpr a Int, Int)], IntMap Int)+cse' lookupFun outputIndices = runST $ do+ ht <- HT.new+ let -- folding function+ f (im,n) [] = return (im,n)+ f (im0,n0) (k:ks) = do+ (_,im,n) <- insertOldNode k lookupFun ht im0 n0+ f (im,n) ks+ -- outputs+ (oldToNewIdx,_) <- f (IM.empty,0) (concatMap toList outputIndices)+ htList <- HT.toList ht+ return (htList, oldToNewIdx)++ +---- | take in an Int that represents a node in the original graph+---- see if that int has been inserted in the new graph+insertOldNode ::+ (Eq a, Hashable a)+ => Int -- ^ Int to be inserted+ -> (Int -> Maybe (GExpr a Int)) -- ^ function to lookup old GExpr from old Int reference+ -> HashTable s (GExpr a Int) Int -- ^ hashmap of new GExprs to their new Int references+ -> IntMap Int -- ^ intmap of old int reference to new int references+ -> Int -- ^ next free index+ -> ST s (Int, IntMap Int, Int)+insertOldNode kOld lookupOldGExpr ht oldNodeToNewNode0 nextFreeInt0 =+ case IM.lookup kOld oldNodeToNewNode0 of+ -- if the int has already been inserted in the new graph, return it+ Just k -> return (k, oldNodeToNewNode0, nextFreeInt0)+ -- if the int has not yet been inserted, then insert it+ -- get the old GExpr to which this node corresponds+ Nothing -> case lookupOldGExpr kOld of+ Nothing -> error $ "in CSE, insertOldNode got an old key \"" ++ show kOld +++ "\" with was not found in the old graph"+ -- insert this old GExpr+ Just oldGExpr -> do+ (k, oldNodeToNewNode1, nextFreeInt1) <- insertOldGExpr oldGExpr lookupOldGExpr ht oldNodeToNewNode0 nextFreeInt0+ return (k, IM.insert kOld k oldNodeToNewNode1, nextFreeInt1)++insertOldGExpr ::+ (Eq a, Hashable a)+ => GExpr a Int -- ^ GExpr to be inserted+ -> (Int -> Maybe (GExpr a Int)) -- ^ function to lookup old GExpr from old Int reference+ -> HashTable s (GExpr a Int) Int -- ^ hashmap of new GExprs to their new Int references+ -> IntMap Int -- ^ intmap of old int reference to new int references+ -> Int -- ^ next free index+ -> ST s (Int, IntMap Int, Int)++insertOldGExpr g@(GSym _) = \_ -> cseInsert g+insertOldGExpr g@(GConst _) = \_ -> cseInsert g+insertOldGExpr g@(GNum (FromInteger _)) = \_ -> cseInsert g+insertOldGExpr g@(GFractional (FromRational _)) = \_ -> cseInsert g++insertOldGExpr (GNum (Mul x y)) = insertOldGExprBinary GNum Mul x y+insertOldGExpr (GNum (Add x y)) = insertOldGExprBinary GNum Add x y+insertOldGExpr (GNum (Sub x y)) = insertOldGExprBinary GNum Sub x y+insertOldGExpr (GFractional (Div x y)) = insertOldGExprBinary GFractional Div x y+insertOldGExpr (GFloating (Pow x y)) = insertOldGExprBinary GFloating Pow x y+insertOldGExpr (GFloating (LogBase x y)) = insertOldGExprBinary GFloating LogBase x y+ +insertOldGExpr (GNum (Negate x)) = insertOldGExprUnary GNum Negate x+insertOldGExpr (GNum (Abs x)) = insertOldGExprUnary GNum Abs x+insertOldGExpr (GNum (Signum x)) = insertOldGExprUnary GNum Signum x+insertOldGExpr (GFloating (Exp x)) = insertOldGExprUnary GFloating Exp x+insertOldGExpr (GFloating (Log x)) = insertOldGExprUnary GFloating Log x+insertOldGExpr (GFloating (Sin x)) = insertOldGExprUnary GFloating Sin x+insertOldGExpr (GFloating (Cos x)) = insertOldGExprUnary GFloating Cos x+insertOldGExpr (GFloating (ASin x)) = insertOldGExprUnary GFloating ASin x+insertOldGExpr (GFloating (ATan x)) = insertOldGExprUnary GFloating ATan x+insertOldGExpr (GFloating (ACos x)) = insertOldGExprUnary GFloating ACos x+insertOldGExpr (GFloating (Sinh x)) = insertOldGExprUnary GFloating Sinh x+insertOldGExpr (GFloating (Cosh x)) = insertOldGExprUnary GFloating Cosh x+insertOldGExpr (GFloating (Tanh x)) = insertOldGExprUnary GFloating Tanh x+insertOldGExpr (GFloating (ASinh x)) = insertOldGExprUnary GFloating ASinh x+insertOldGExpr (GFloating (ATanh x)) = insertOldGExprUnary GFloating ATanh x+insertOldGExpr (GFloating (ACosh x)) = insertOldGExprUnary GFloating ACosh x++insertOldGExprBinary ::+ (Eq a, Hashable a)+ => (f -> GExpr a Int)+ -> (Int -> Int -> f)+ -> Int -> Int+ -> (Int -> Maybe (GExpr a Int)) -- ^ function to lookup old GExpr from old Int reference+ -> HashTable s (GExpr a Int) Int -- ^ hashmap of new GExprs to their new Int references+ -> IntMap Int -- ^ intmap of old int reference to new int references+ -> Int -- ^ next free index+ -> ST s (Int, IntMap Int, Int)+insertOldGExprBinary gnum mul kxOld kyOld lookupOldGExpr ht oldNodeToNewNode0 nextFreeInt0 = do+ (kx, oldNodeToNewNode1,nextFreeInt1) <- insertOldNode kxOld lookupOldGExpr ht oldNodeToNewNode0 nextFreeInt0+ (ky, oldNodeToNewNode2,nextFreeInt2) <- insertOldNode kyOld lookupOldGExpr ht oldNodeToNewNode1 nextFreeInt1+ let newGExpr = gnum (mul kx ky)+ cseInsert newGExpr ht oldNodeToNewNode2 nextFreeInt2++insertOldGExprUnary ::+ (Eq a, Hashable a)+ => (f -> GExpr a Int)+ -> (Int -> f)+ -> Int+ -> (Int -> Maybe (GExpr a Int)) -- ^ function to lookup old GExpr from old Int reference+ -> HashTable s (GExpr a Int) Int -- ^ hashmap of new GExprs to their new Int references+ -> IntMap Int -- ^ intmap of old int reference to new int references+ -> Int -- ^ next free index+ -> ST s (Int, IntMap Int, Int)+insertOldGExprUnary gnum neg kxOld lookupOldGExpr ht oldNodeToNewNode0 nextFreeInt0 = do+ (kx, oldNodeToNewNode1,nextFreeInt1) <- insertOldNode kxOld lookupOldGExpr ht oldNodeToNewNode0 nextFreeInt0+ let newGExpr = gnum (neg kx)+ cseInsert newGExpr ht oldNodeToNewNode1 nextFreeInt1++cseInsert :: (Eq a, Hashable a) => GExpr a Int -> HashTable s (GExpr a Int) Int -> IntMap Int -> Int+ -> ST s (Int, IntMap Int, Int)+cseInsert gexpr ht oldNodeToNewNode0 nextFreeInt0 = do+ lu <- HT.lookup ht gexpr+ case lu of+ Just k -> return (k, oldNodeToNewNode0, nextFreeInt0)+ Nothing -> do+ HT.insert ht gexpr nextFreeInt0+ return (nextFreeInt0, oldNodeToNewNode0, nextFreeInt0+1)+
− Dvda/CallNative.hs
@@ -1,202 +0,0 @@-{-# 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 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.HashMap ( HashMap )-import qualified Dvda.HashMap as HM-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
@@ -1,38 +0,0 @@-{-# 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/Codegen/Gcc.hs view
@@ -0,0 +1,32 @@+{-# OPTIONS_GHC -Wall #-}++module Dvda.Codegen.Gcc ( compileWithGcc+ ) where++import System.Process(runCommand, waitForProcess)+import System.Exit(ExitCode(ExitSuccess))+import Control.Monad(when)++-- | whether to print the gcc call when generating code+spewGccCall :: Bool+spewGccCall = True++-- | take in source file and object, return string suitible for calling to compile+gccString :: FilePath -> FilePath -> String+gccString src obj = "gcc -O2 -std=gnu99 -fPIC -shared -Wall -Wextra -Werror " ++ src ++ " -o " ++ obj++-- | take in name of source and future object, compile object+compileWithGcc :: FilePath -> FilePath -> IO ()+compileWithGcc srcname objname = do+ -- compile new object+ let compileString = gccString srcname objname++ -- print compilation string+ when spewGccCall $ putStrLn compileString+ + -- run compilation string+ p <- runCommand compileString+ + -- check for errors+ exitCode <- waitForProcess p+ when (exitCode /= ExitSuccess) $ error $ "failed compiling " ++ srcname
+ Dvda/Codegen/WriteFile.hs view
@@ -0,0 +1,32 @@+{-# OPTIONS_GHC -Wall #-}++module Dvda.Codegen.WriteFile ( writeSourceFile+ ) where++import System.Directory+import Control.Monad ( when )++---- | return directory to use for temp files+---- | create this directory and print message if it doesn't exist+--dvdaDir :: IO FilePath+--dvdaDir = do+-- dir <- getAppUserDataDirectory "dvda"++writeSourceFile :: String -> FilePath -> FilePath -> IO FilePath+writeSourceFile source functionDir sourceName = do+ -- make function directory if it doesn't exist+ createDirectoryIfMissing False functionDir+ + -- filenames+ let sourcePath = functionDir ++ "/" ++ sourceName+ + -- if the source already exists, make sure it matches the old source+ srcExists <- doesFileExist sourcePath+ when srcExists $ do+ putStrLn $ "file \"" ++ sourcePath ++ "\" already exists, overwriting"+ + -- write source+ putStrLn $ "writing " ++ sourcePath+ writeFile sourcePath source++ return sourcePath
− Dvda/Config.hs
@@ -1,125 +0,0 @@--- Config.hs--{-# OPTIONS_GHC -Wall #-}--module Dvda.Config( -- * directory stuff- dvdaDir- -- * C syntax- , cType- , cName- , nameCSource- , nameCInclude- , nameCObject- , nameCFunction- -- * Haskell syntax- , nameHSObject- , nameHSModule- , nameHSFunction- , nameHSSource- , nameHSVar- , nameHSConst- -- * Octave- , nameOctaveSource- , nameOctaveFunction- -- * gcc stuff- , gccString- , spewGccCall- , outputNames- -- * ghc stuff- , ghcString- -- * symbolic stuff- , simplifyCommutativeOps- ) where--import System.Directory-import Control.Monad(unless)---- | what symbolic variable names to use when generating a function-outputNames :: [String]-outputNames = map (\x -> "out"++show x) [(0::Integer)..]---- | whether to print the gcc call when generating code-spewGccCall :: Bool-spewGccCall = True---- | return directory to use for temp files--- | create this directory and print message if it doesn't exist-dvdaDir :: IO FilePath-dvdaDir = do- dir <- getAppUserDataDirectory "dvda"- - -- print message if creating directory- exist <- doesDirectoryExist dir- unless exist $ putStrLn $ "creating directory \""++dir++"\" for codegen source/objects"-- -- make the directory if missing- createDirectoryIfMissing True dir- - return dir----- | take in source file and object, return string suitible for calling to compile-gccString :: FilePath -> FilePath -> String-gccString src obj = "gcc -O2 -std=gnu99 -fPIC -shared -Wall -Wextra -Werror " ++ src ++ " -o " ++ obj--ghcString :: FilePath -> FilePath -> String-ghcString src obj = "ghc -c " ++ src ++ " -o " ++ obj----- c syntax--- | type to use when generating c code-cType :: String-cType = "double"---- | name convention for c variables-cName :: Int -> String-cName k- | k < 0 = error "cName got negative index" - | otherwise = 't':show k --nameCSource :: String -> String-nameCSource hash = nameCFunction hash ++ ".c"--nameCInclude :: String -> String-nameCInclude hash = nameCFunction hash ++ ".h"--nameCObject :: String -> String-nameCObject hash = "c_" ++ nameCFunction hash ++ ".o"--nameCFunction :: String -> String-nameCFunction hash = "call_" ++ hash---- haskell syntax-nameHSVar :: Int -> String-nameHSVar k- | k < 0 = error "nameHSVar got negative index" - | otherwise = 'v':show k---- haskell syntax-nameHSConst :: Int -> String-nameHSConst k- | k < 0 = error "nameHSConst got negative index" - | otherwise = 'c':show k--nameHSFunction :: String -> String-nameHSFunction hash = "call_" ++ hash--nameHSModule :: String -> String-nameHSModule hash = "Call_" ++ hash--nameHSSource :: String -> String-nameHSSource = (++ ".hs") . nameHSModule--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/Examples.hs view
@@ -1,151 +1,54 @@ {-# OPTIONS_GHC -Wall #-}-{-# Language TypeOperators #-} -module Dvda.Examples ( run- , run'- , showoff- , bigGraph- , smallGraph- , runCallNative- , composed+module Dvda.Examples ( doCse+ , showFg+ , cgen+ , mexgen ) where -import Data.Array.Repa.Index-import Control.Monad.State--import Dvda import Dvda.Expr-import Dvda.CallNative-import Dvda.Graph ( FunGraph(..) )--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 $ y**2- z3 <- node $ diff ((x*x/2)**x) x- - outputs_ (z1 :* z2 :* z3)--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"- - args = x :* y :* z- - inputs_ args- outputs_ (pureFun args)--run' :: IO ()-run' = do- let gr@(FunGraph hm im _ _) = runFunGraph exampleFunGraph- (FunGraph hm' im' _ _) = runFunGraph exampleFunGraph'- - putStrLn $ funGraphSummary' gr- putStrLn $ showCollisions gr- previewGraph gr- putStrLn "\nimperative same as pure+cse?:"- print $ hm == hm'- print $ im == im'--run :: IO ()-run = do- let gr@( FunGraph _ _ _ _) = runFunGraph $ do- let x = sym "x" :: Expr DIM0 Double- y = sym "y"- z1 = x + x / y + 3- z2 = diff z1 x- z3 = diff z1 y-- inputs_ (x :* y)- outputs_ (z1 :* z2 :* z3)-- putStrLn $ showCollisions 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--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'+import Dvda.FunGraph+import Dvda.CGen+import Dvda.Vis ( previewGraph )+import Dvda.CSE ( cse )+import Dvda.AD ( rad ) -smallGraph :: FunGraph Double- (Exprs (DIM0 :* DIM0 :* DIM0) Double)- (Exprs (DIM0 :* DIM0) Double)-smallGraph = makeFunGraph (x :* y :* z) (f0 :* f1)+-- a random function to use in different examples+someFunGraph :: IO (FunGraph Double)+someFunGraph = toFunGraph inputs outputs where- x = sym "x" :: Expr DIM0 Double+ x = sym "x" :: Expr Double y = sym "y" z = sym "z"+ w = sym "w"+ w1 = sym "w1"+ w2 = sym "w2"+ w3 = sym "w3"+ f0 = x*y + z + w1 + w2+ f2 = f0 * w2/w3+ + f1 = [f0/2, f0*y, w, 0.0, 0]+ boo = x - f0 = x*y*z + 3- f1 = 40*f0/x+ inputs = boo :* [y]:*[[z]] :* [w3,w1,w2,w]+ outputs = f0:*f1:*f2:*[[f0*f0]]:*(rad f2 [x,y,z,w,w1,w2,w3]) -runCallNative :: Exprs (Z :* Z) Double-runCallNative = toNative smallGraph (f 1 :* f 2 :* f 3)- where- f = EConst . (CSingleton Z)+-- | do cse on a fungraph and count nodes+doCse :: IO ()+doCse = do+ fg' <- someFunGraph+ putStrLn $ "fungraph has " ++ show (countNodes fg') ++ " nodes"+ let fg = cse fg'+ putStrLn $ "fungraph has " ++ show (countNodes fg) ++ " nodes after cse" -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+-- | show a fungraph+showFg :: IO ()+showFg = someFunGraph >>= previewGraph -composed :: [Expr Z Double]-composed = runDeriv z [t]- where- t = sym "t"- x = symDependent "x" t- y = symDependent "y" x- z = symDependent "z" y+-- | c code generation+cgen :: IO ()+cgen = fmap (showC RowMajor "foo") someFunGraph >>= putStrLn++-- | mex function generation+mexgen :: IO ()+mexgen = fmap (showMex "foo") someFunGraph >>= putStrLn
Dvda/Expr.hs view
@@ -1,71 +1,44 @@-{-# Options_ghc -Wall #-}+{-# OPTIONS_GHC -Wall #-}+{-# Language GADTs #-}+{-# Language TemplateHaskell #-}+{-# Language TypeFamilies #-} {-# Language StandaloneDeriving #-} {-# Language DeriveDataTypeable #-}-{-# Language GADTs #-}+{-# Language FlexibleInstances #-} {-# Language FlexibleContexts #-} module Dvda.Expr ( Expr(..)- , Const(..)+ , GExpr(..)+ , Nums(..)+ , Fractionals(..)+ , Floatings(..) , Sym(..)- , RefHash(..)- , sym- , svec- , smat- , vsym- , msym- , vec- , mat- , scale--- , dot- , diff- , grad- , jacob- , hess- , dim , isVal+ , sym , symDependent , symDependentN+ , const'+ , getParents+ , extractLinearPart+ , getConst+ , substitute+ , sketchySubstitute ) where -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 Control.Applicative ( (<$>), (<*>), pure )+import Data.Data ( Data, Typeable, Typeable1, Typeable2 ) import Data.Hashable ( Hashable, hash, combine )-import Data.List ( sort )-import Data.Typeable ( Typeable2 ) -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+--import Test.QuickCheck -- ( Arbitrary(..) ) -dim :: Expr sh a -> sh-dim (ESym 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 (EUnary _ x) = dim x-dim (EBinary _ x1 _) = dim x1-dim (EScale _ y) = 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))+import qualified Dvda.HashMap as HM+import Dvda.Reify ( MuRef(..) ) -deriving instance Typeable2 Const-deriving instance Typeable2 Expr+commutativeMul :: Bool+commutativeMul = True -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+commutativeAdd :: Bool+commutativeAdd = True 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@@ -75,328 +48,591 @@ show (Sym name) = name show (SymDependent name k s) = name ++ replicate k '\'' ++ "(" ++ show s ++ ")" -data RefHash = RefHash Int deriving (Eq, Show)+data Expr a where+ ESym :: Sym -> Expr a+ EConst :: a -> Expr a+ ENum :: Num a => Nums (Expr a) -> Expr a+ EFractional :: Fractional a => Fractionals (Expr a) -> Expr a+ EFloating :: Floating a => Floatings (Expr a) -> Expr a -data Expr sh a where- ESym :: sh -> Sym -> Expr sh a- EConst :: Const sh a -> Expr sh a- EDimensionless :: 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- ERef :: sh -> RefHash -> Key -> Expr sh a+data Nums a = Mul a a+ | Add a a+ | Sub a a+ | Negate a+ | Abs a+ | Signum a+ | FromInteger Integer - EDeriv :: Expr DIM0 a -> Expr DIM0 a -> Expr DIM0 a- EGrad :: Expr DIM0 a -> Expr sh a -> Expr sh a- EJacob :: Expr DIM1 a -> Expr DIM1 a -> Expr DIM2 a+data Fractionals a = Div a a+ | FromRational Rational deriving Eq ---------------------------------- 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+data Floatings a = Pow a a+ | LogBase a a+ | Exp a+ | Log a+ | Sin a+ | Cos a+ | ASin a+ | ATan a+ | ACos a+ | Sinh a+ | Cosh a+ | Tanh a+ | ASinh a+ | ATanh a+ | ACosh a deriving Eq -paren :: String -> String-paren x = "("++ x ++")"+deriving instance Data Sym+deriving instance Data a => Data (Nums a)+deriving instance Data a => Data (Fractionals a)+deriving instance Data a => Data (Floatings a)+deriving instance (Data a, Floating a) => Data (Expr a)+deriving instance (Data a, Data b, Floating a) => Data (GExpr a b) -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+deriving instance Typeable Sym+deriving instance Typeable1 Nums+deriving instance Typeable1 Fractionals+deriving instance Typeable1 Floatings+deriving instance Typeable1 Expr+deriving instance Typeable2 GExpr - 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+----------------------- Show instances -------------------------+showsInfixBinary :: (Show a, Show b) => Int -> Int -> String -> a -> b -> ShowS+showsInfixBinary d prec op u v = showParen (d > prec) $+ showsPrec prec u .+ showString op .+ showsPrec prec v - 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 ++ ")"+showsUnary :: Show a => Int -> Int -> String -> a -> ShowS+showsUnary d prec op u = showParen (d > prec) $+ showString op .+ showsPrec prec u +instance Show a => Show (Nums a) where+ showsPrec d (Mul x y) = showsInfixBinary d 7 " * " x y+ showsPrec d (Add x y) = showsInfixBinary d 6 " + " x y+ showsPrec d (Sub x y) = showsInfixBinary d 6 " - " x y+ showsPrec d (Negate x) = showsUnary d 7 "-" x+ showsPrec d (Abs x) = showsUnary d 10 "abs" x+ showsPrec d (Signum x) = showsUnary d 10 "signum" x+ showsPrec _ (FromInteger k) = showString (show k) ---------------------------------- 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+instance Show a => Show (Fractionals a) where+ showsPrec d (Div x y) = showsInfixBinary d 7 " / " x y+ showsPrec _ (FromRational r) = showString $ show (fromRational r :: Double)++instance Show a => Show (Floatings a) where+ showsPrec d (Pow x y) = showsInfixBinary d 8 " ** " x y+ showsPrec d (LogBase x y) = showParen (d > 10) $ showString $ "logBase(" ++ show x ++ ", " ++ show y ++ ")"+ showsPrec d (Exp x) = showsUnary d 10 "exp" x+ showsPrec d (Log x) = showsUnary d 10 "log" x+ showsPrec d (Sin x) = showsUnary d 10 "sin" x+ showsPrec d (Cos x) = showsUnary d 10 "cos" x+ showsPrec d (ASin x) = showsUnary d 10 "asin" x+ showsPrec d (ATan x) = showsUnary d 10 "atan" x+ showsPrec d (ACos x) = showsUnary d 10 "acos" x+ showsPrec d (Sinh x) = showsUnary d 10 "sinh" x+ showsPrec d (Cosh x) = showsUnary d 10 "cosh" x+ showsPrec d (Tanh x) = showsUnary d 10 "tanh" x+ showsPrec d (ASinh x) = showsUnary d 10 "asinh" x+ showsPrec d (ATanh x) = showsUnary d 10 "atanh" x+ showsPrec d (ACosh x) = showsUnary d 10 "acosh" x++instance Show a => Show (Expr a) where+ showsPrec _ (ESym s) = showString (show s)+ showsPrec _ (EConst x) = showString (show x)+ showsPrec d (ENum x) = showsPrec d x+ showsPrec d (EFractional x) = showsPrec d x+ showsPrec d (EFloating x) = showsPrec d x+++----------------------- Eq instances -------------------------+instance Eq a => Eq (Expr a) where+ (==) (ESym x) (ESym y) = x == y+ (==) (EConst x) (EConst y) = x == y+ (==) (ENum x) (ENum y) = x == y+ (==) (EFractional x) (EFractional y) = x == y+ (==) (EFloating x) (EFloating y) = x == y (==) _ _ = False++instance Eq a => Eq (Nums a) where+ (Mul x0 y0) == (Mul x1 y1) = if commutativeMul+ then (x0 == x1 && y0 == y1) || (x0 == y1 && x1 == y0)+ else x0 == x1 && y0 == y1+ (Add x0 y0) == (Add x1 y1) = if commutativeAdd+ then (x0 == x1 && y0 == y1) || (x0 == y1 && x1 == y0)+ else x0 == x1 && y0 == y1+ (Sub x0 y0) == (Sub x1 y1) = x0 == x1 && y0 == y1+ (Negate x) == (Negate y) = x == y+ (Abs x) == (Abs y) = x == y+ (Signum x) == (Signum y) = x == y+ (FromInteger x) == (FromInteger y) = x == y+ _ == _ = 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 h0 k0) (ERef sh1 h1 k1) = sh0 == sh1 && h0 == h1 && 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++----------------------------- hashable instances --------------------------+instance Hashable Sym where+ hash (Sym name) = hash "Sym" `combine` hash name+ hash (SymDependent name k s) = hash ("SymDependent", name, k, s)++instance Hashable a => Hashable (Nums a) where+ hash (Mul x y) = hash "Mul" `combine` hx `combine` hy where- commutativeEq- | simplifyCommutativeOps && isCommutative op0 = (x0 == x1 && y0 == y1) || (x0 == y1 && y0 == x1)- | otherwise = x0 == x1 && y0 == y1- (==) _ _ = False+ hx' = hash x+ hy' = hash y+ (hx, hy)+ | commutativeMul = (min hx' hy', max hx' hy')+ | otherwise = (hx', hy')+ hash (Add x y) = hash "Add" `combine` hx `combine` hy+ where+ hx' = hash x+ hy' = hash y+ (hx, hy)+ | commutativeAdd = (min hx' hy', max hx' hy')+ | otherwise = (hx', hy')+ hash (Sub x y) = hash "Sub" `combine` hash x `combine` hash y+ hash (Negate x) = hash "Negate" `combine` hash x+ hash (Abs x) = hash "Abs" `combine` hash x+ hash (Signum x) = hash "Signum" `combine` hash x+ hash (FromInteger x) = hash "FromInteger" `combine` hash x -------------------------- 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 => Hashable (Fractionals a) where+ hash (Div x y) = hash "Div" `combine` hash x `combine` hash y+ hash (FromRational x) = hash "FromRational" `combine` hash x +instance Hashable a => Hashable (Floatings a) where+ hash (Pow x y) = hash "Pow" `combine` hash x `combine` hash y+ hash (LogBase x y) = hash "LogBase" `combine` hash x `combine` hash y+ hash (Exp x) = hash "Exp" `combine` hash x+ hash (Log x) = hash "Log" `combine` hash x+ hash (Sin x) = hash "Sin" `combine` hash x+ hash (Cos x) = hash "Cos" `combine` hash x+ hash (ASin x) = hash "ASin" `combine` hash x+ hash (ATan x) = hash "ATan" `combine` hash x+ hash (ACos x) = hash "ACos" `combine` hash x+ hash (Sinh x) = hash "Sinh" `combine` hash x+ hash (Cosh x) = hash "Cosh" `combine` hash x+ hash (Tanh x) = hash "Tanh" `combine` hash x+ hash (ASinh x) = hash "ASinh" `combine` hash x+ hash (ATanh x) = hash "ATanh" `combine` hash x+ hash (ACosh x) = hash "ACosh" `combine` hash x -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 _ (RefHash h) _) = h+instance Hashable a => Hashable (Expr a) where+ hash (ESym name) = hash "ESym" `combine` hash name+ hash (EConst x) = hash "EConst" `combine` hash x+ hash (ENum x) = hash "ENum" `combine` hash x+ hash (EFractional x) = hash "EFractional" `combine` hash x+ hash (EFloating x) = hash "EFloating" `combine` hash x - 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+--deriving instance Enum a => Enum (Nums a)+--deriving instance Bounded a => Bounded (Nums a) -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+--deriving instance Enum a => Enum (Fractionals a)+--deriving instance Bounded a => Bounded (Fractionals a) +--deriving instance Enum a => Enum (Floatings a)+--deriving instance Bounded a => Bounded (Floatings a) ------------------------- symbolic stuff ---------------------isVal :: Eq a => a -> Expr sh a -> Bool-isVal x (EDimensionless y) = x == y-isVal x (EConst (CSingleton _ y)) = x == y-isVal _ _ = False+instance (Num a, Eq a) => Num (Expr a) where+ (*) (EConst x) (EConst y) = EConst (x*y)+ (*) (ENum (FromInteger kx)) (ENum (FromInteger ky)) = ENum $ FromInteger (kx * ky)+ (*) (EFractional (FromRational rx)) (EFractional (FromRational ry)) = EFractional $ FromRational (rx * ry)+ (*) (EConst x) (ENum (FromInteger ky)) = EConst $ x * fromInteger ky+ (*) (ENum (FromInteger kx)) (EConst y) = EConst $ fromInteger kx * y+ (*) (EConst x) (EFractional (FromRational ry)) = EConst $ x * fromRational ry+ (*) (EFractional (FromRational rx)) (EConst y) = EConst $ fromRational rx * y+ (*) (ENum (FromInteger kx)) (EFractional (FromRational ry)) = EFractional $ FromRational (fromInteger kx * ry)+ (*) (EFractional (FromRational rx)) (ENum (FromInteger ky)) = EFractional $ FromRational (rx * fromInteger ky)+ (*) x y+ | isVal 0 x || isVal 0 y = 0+ | isVal 1 x = y+ | isVal 1 y = x+ | otherwise = ENum $ Mul x y --- | first layer of binary simplification: infer dimension of EDimensionless if possible-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 (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+ (+) (EConst x) (EConst y) = EConst (x+y)+ (+) (ENum (FromInteger kx)) (ENum (FromInteger ky)) = ENum $ FromInteger (kx + ky)+ (+) (EFractional (FromRational rx)) (EFractional (FromRational ry)) = EFractional $ FromRational (rx + ry)+ (+) (EConst x) (ENum (FromInteger ky)) = EConst $ x + fromInteger ky+ (+) (ENum (FromInteger kx)) (EConst y) = EConst $ fromInteger kx + y+ (+) (EConst x) (EFractional (FromRational ry)) = EConst $ x + fromRational ry+ (+) (EFractional (FromRational rx)) (EConst y) = EConst $ fromRational rx + y+ (+) (ENum (FromInteger kx)) (EFractional (FromRational ry)) = EFractional $ FromRational (fromInteger kx + ry)+ (+) (EFractional (FromRational rx)) (ENum (FromInteger ky)) = EFractional $ FromRational (rx + fromInteger ky)+ (+) x y+ | isVal 0 x = y+ | isVal 0 y = x+ | otherwise = ENum $ Add x y + (-) (EConst x) (EConst y) = EConst (x-y)+ (-) (ENum (FromInteger kx)) (ENum (FromInteger ky)) = ENum $ FromInteger (kx - ky)+ (-) (EFractional (FromRational rx)) (EFractional (FromRational ry)) = EFractional $ FromRational (rx - ry)+ (-) (EConst x) (ENum (FromInteger ky)) = EConst $ x - fromInteger ky+ (-) (ENum (FromInteger kx)) (EConst y) = EConst $ fromInteger kx - y+ (-) (EConst x) (EFractional (FromRational ry)) = EConst $ x - fromRational ry+ (-) (EFractional (FromRational rx)) (EConst y) = EConst $ fromRational rx - y+ (-) (ENum (FromInteger kx)) (EFractional (FromRational ry)) = EFractional $ FromRational (fromInteger kx - ry)+ (-) (EFractional (FromRational rx)) (ENum (FromInteger ky)) = EFractional $ FromRational (rx - fromInteger ky)+ (-) x y+ | isVal 0 x = negate y+ | isVal 0 y = x+ | otherwise = ENum $ Sub x y --- | second layer of binary simplification: check dimensions-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 ++ ")"- where- shx = dim x- shy = dim y- sdx = showShapeR shx- sdy = showShapeR shy- sop = show op+ abs (EConst x) = EConst (abs x)+ abs (ENum (FromInteger k)) = ENum (FromInteger (abs k))+ abs (EFractional (FromRational r)) = EFractional (FromRational (abs r))+ abs x = ENum $ Abs x + negate (EConst x) = EConst (negate x)+ negate (ENum (FromInteger k)) = ENum (FromInteger (negate k))+ negate (EFractional (FromRational r)) = EFractional (FromRational (negate r))+ negate x = ENum $ Negate x --- | 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'' :: (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- | isVal 1 x = y- | isVal 1 y = x- | otherwise = makeBinary''' Mul f x y-makeBinary'' Add f x y- | isVal 0 x = y- | isVal 0 y = x- | otherwise = makeBinary''' Add f x y-makeBinary'' Div f x y- | isVal 0 y = error "divide by zero"- | isVal 1 y = x- | isVal 0 x = x- | otherwise = makeBinary''' Div f x y-makeBinary'' Sub f x y- | isVal 0 x = negate y- | isVal 0 y = x- | otherwise = makeBinary''' Sub f x y-makeBinary'' op f x y = makeBinary''' op f x y+ signum (EConst x) = EConst (signum x)+ signum (ENum (FromInteger k)) = ENum (FromInteger (signum k))+ signum (EFractional (FromRational r)) = EFractional (FromRational (signum r))+ signum x = ENum $ Signum x + fromInteger = ENum . FromInteger --- | fourth layer of binary simplification: make reasonable simplifications-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+instance (Fractional a, Eq a) => Fractional (Expr a) where+ (/) (EConst x) (EConst y) = EConst (x/y)+-- (/) (ENum (FromInteger kx)) (ENum (FromInteger ky)) = ENum $ FromInteger (kx / ky)+ (/) (EFractional (FromRational rx)) (EFractional (FromRational ry)) = EFractional $ FromRational (rx / ry)+ (/) (EConst x) (ENum (FromInteger ky)) = EConst $ x / fromInteger ky+ (/) (ENum (FromInteger kx)) (EConst y) = EConst $ fromInteger kx / y+ (/) (EConst x) (EFractional (FromRational ry)) = EConst $ x / fromRational ry+ (/) (EFractional (FromRational rx)) (EConst y) = EConst $ fromRational rx / y+ (/) (ENum (FromInteger kx)) (EFractional (FromRational ry)) = EFractional $ FromRational (fromInteger kx / ry)+ (/) (EFractional (FromRational rx)) (ENum (FromInteger ky)) = EFractional $ FromRational (rx / fromInteger ky)+ (/) x y+ | isVal 0 y = error "Fractional (Expr a) divide by zero"+ | isVal 0 x = 0+ | isVal 1 y = x+ | otherwise = EFractional $ Div x y + fromRational = EFractional . FromRational --- | apply unary operations on constants-makeUnary :: Storable a => UnOp -> (a -> a) -> Expr sh a -> Expr sh a-makeUnary _ f (EDimensionless x) = EDimensionless (f x)-makeUnary _ f' (EConst x') = EConst $ cmap f' x'+instance (Floating a, Eq a) => Floating (Expr a) where+ pi = EConst pi+ x ** y = EFloating $ Pow x y+ logBase x y = EFloating $ LogBase x y+ exp = applyFloatingUn ( exp, Exp)+ log = applyFloatingUn ( log, Log)+ sin = applyFloatingUn ( sin, Sin)+ cos = applyFloatingUn ( cos, Cos)+ asin = applyFloatingUn ( asin, ASin)+ atan = applyFloatingUn ( atan, ATan)+ acos = applyFloatingUn ( acos, ACos)+ sinh = applyFloatingUn ( sinh, Sinh)+ cosh = applyFloatingUn ( cosh, Cosh)+ tanh = applyFloatingUn ( tanh, Tanh)+ asinh = applyFloatingUn (asinh, ASinh)+ atanh = applyFloatingUn (atanh, ATanh)+ acosh = applyFloatingUn (acosh, ACosh)++applyFloatingUn :: Floating a => (t -> a, Expr t -> Floatings (Expr a)) -> Expr t -> Expr a+applyFloatingUn (f,_) (EConst x) = EConst (f x)+applyFloatingUn (f,_) (ENum (FromInteger x)) = EConst (f $ fromInteger x)+applyFloatingUn (f,_) (EFractional (FromRational x)) = EConst (f $ fromRational x)+applyFloatingUn (_,f) x = EFloating (f x)++---------------------------------- GExprs --------------------------------+data GExpr a b where+ GSym :: Sym -> GExpr a b+ GConst :: a -> GExpr a b+ GNum :: Num a => Nums b -> GExpr a b+ GFractional :: Fractional a => Fractionals b -> GExpr a b+ GFloating :: Floating a => Floatings b -> GExpr a b++-- you might use this to use Expr's nice Show instance+gexprToExpr :: (b -> Expr a) -> GExpr a b -> Expr a+gexprToExpr _ (GSym s@(Sym _)) = ESym s+gexprToExpr _ (GSym sd@(SymDependent _ _ _)) = ESym sd+gexprToExpr _ (GConst c) = EConst c+gexprToExpr f (GNum (Mul x y)) = ENum (Mul (f x) (f y))+gexprToExpr f (GNum (Add x y)) = ENum (Add (f x) (f y))+gexprToExpr f (GNum (Sub x y)) = ENum (Sub (f x) (f y))+gexprToExpr f (GNum (Negate x)) = ENum (Negate (f x))+gexprToExpr f (GNum (Abs x)) = ENum (Abs (f x))+gexprToExpr f (GNum (Signum x)) = ENum (Signum (f x))+gexprToExpr _ (GNum (FromInteger x)) = ENum (FromInteger x)+gexprToExpr f (GFractional (Div x y)) = EFractional (Div (f x) (f y))+gexprToExpr _ (GFractional (FromRational x)) = EFractional (FromRational x)+gexprToExpr f (GFloating (Pow x y)) = EFloating (Pow (f x) (f y))+gexprToExpr f (GFloating (LogBase x y)) = EFloating (LogBase (f x) (f y))+gexprToExpr f (GFloating (Exp x)) = EFloating (Exp (f x))+gexprToExpr f (GFloating (Log x)) = EFloating (Log (f x))+gexprToExpr f (GFloating (Sin x)) = EFloating (Sin (f x))+gexprToExpr f (GFloating (Cos x)) = EFloating (Cos (f x))+gexprToExpr f (GFloating (ASin x)) = EFloating (ASin (f x))+gexprToExpr f (GFloating (ATan x)) = EFloating (ATan (f x))+gexprToExpr f (GFloating (ACos x)) = EFloating (ACos (f x))+gexprToExpr f (GFloating (Sinh x)) = EFloating (Sinh (f x))+gexprToExpr f (GFloating (Cosh x)) = EFloating (Cosh (f x))+gexprToExpr f (GFloating (Tanh x)) = EFloating (Tanh (f x))+gexprToExpr f (GFloating (ASinh x)) = EFloating (ASinh (f x))+gexprToExpr f (GFloating (ATanh x)) = EFloating (ATanh (f x))+gexprToExpr f (GFloating (ACosh x)) = EFloating (ACosh (f x))++getParents :: GExpr a b -> [b]+getParents (GSym _) = []+getParents (GConst _) = []+getParents (GNum (Mul x y)) = [x,y]+getParents (GNum (Add x y)) = [x,y]+getParents (GNum (Sub x y)) = [x,y]+getParents (GNum (Negate x)) = [x]+getParents (GNum (Abs x)) = [x]+getParents (GNum (Signum x)) = [x]+getParents (GNum (FromInteger _)) = []+getParents (GFractional (Div x y)) = [x,y]+getParents (GFractional (FromRational _)) = []+getParents (GFloating (Pow x y)) = [x,y]+getParents (GFloating (LogBase x y)) = [x,y]+getParents (GFloating (Exp x)) = [x]+getParents (GFloating (Log x)) = [x]+getParents (GFloating (Sin x)) = [x]+getParents (GFloating (Cos x)) = [x]+getParents (GFloating (ASin x)) = [x]+getParents (GFloating (ATan x)) = [x]+getParents (GFloating (ACos x)) = [x]+getParents (GFloating (Sinh x)) = [x]+getParents (GFloating (Cosh x)) = [x]+getParents (GFloating (Tanh x)) = [x]+getParents (GFloating (ASinh x)) = [x]+getParents (GFloating (ATanh x)) = [x]+getParents (GFloating (ACosh x)) = [x]++instance (Show a, Show b) => Show (GExpr a b) where+ show = show . (gexprToExpr (\x -> ESym (Sym ("{" ++ show x ++ "}"))))+ +deriving instance (Eq a, Eq b) => Eq (GExpr a b)++instance (Hashable a, Hashable b) => Hashable (GExpr a b) where+ hash (GSym name) = hash "GSym" `combine` hash name+ hash (GConst x) = hash "GConst" `combine` hash x+ hash (GNum x) = hash "GNum" `combine` hash x+ hash (GFractional x) = hash "GFractional" `combine` hash x+ hash (GFloating x) = hash "GFloating" `combine` hash x++instance MuRef (Expr a) where+ type DeRef (Expr a) = GExpr a+ mapDeRef _ (ESym name) = pure (GSym name)+ mapDeRef _ (EConst c) = pure (GConst c)+ mapDeRef f (ENum (Mul x y)) = GNum <$> (Mul <$> (f x) <*> (f y))+ mapDeRef f (ENum (Add x y)) = GNum <$> (Add <$> (f x) <*> (f y))+ mapDeRef f (ENum (Sub x y)) = GNum <$> (Sub <$> (f x) <*> (f y))+ mapDeRef f (ENum (Negate x)) = GNum <$> (Negate <$> (f x))+ mapDeRef f (ENum (Abs x)) = GNum <$> (Negate <$> (f x))+ mapDeRef f (ENum (Signum x)) = GNum <$> (Signum <$> (f x))+ mapDeRef _ (ENum (FromInteger k)) = pure $ GNum (FromInteger k)++ mapDeRef f (EFractional (Div x y)) = GFractional <$> (Div <$> (f x) <*> (f y))+ mapDeRef _ (EFractional (FromRational x)) = pure $ GFractional (FromRational x)++ mapDeRef f (EFloating (Pow x y)) = GFloating <$> (Pow <$> (f x) <*> (f y))+ mapDeRef f (EFloating (LogBase x y)) = GFloating <$> (LogBase <$> (f x) <*> (f y))+ mapDeRef f (EFloating (Exp x)) = GFloating <$> (Exp <$> (f x))+ mapDeRef f (EFloating (Log x)) = GFloating <$> (Log <$> (f x))+ mapDeRef f (EFloating (Sin x)) = GFloating <$> (Sin <$> (f x))+ mapDeRef f (EFloating (Cos x)) = GFloating <$> (Cos <$> (f x))+ mapDeRef f (EFloating (ASin x)) = GFloating <$> (ASin <$> (f x))+ mapDeRef f (EFloating (ATan x)) = GFloating <$> (ATan <$> (f x))+ mapDeRef f (EFloating (ACos x)) = GFloating <$> (ACos <$> (f x))+ mapDeRef f (EFloating (Sinh x)) = GFloating <$> (Sinh <$> (f x))+ mapDeRef f (EFloating (Cosh x)) = GFloating <$> (Cosh <$> (f x))+ mapDeRef f (EFloating (Tanh x)) = GFloating <$> (Tanh <$> (f x))+ mapDeRef f (EFloating (ASinh x)) = GFloating <$> (ASinh <$> (f x))+ mapDeRef f (EFloating (ATanh x)) = GFloating <$> (ATanh <$> (f x))+ mapDeRef f (EFloating (ACosh x)) = GFloating <$> (ACosh <$> (f x))++substitute :: (Eq a, Hashable a, Show a) => Expr a -> [(Expr a, Expr a)] -> Expr a+substitute expr subList+ | nonSymInputs /= [] = error $ "substitute got non-ESym input: " ++ show nonSymInputs+ | otherwise = subs expr 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+ isSym (ESym _) = True+ isSym _ = False+ nonSymInputs = filter (not . isSym . fst) subList+ lookup' e = let hm = HM.fromList subList in+ HM.lookupDefault e e hm+ + subs e@(ESym _) = lookup' e+ subs e@(EConst _) = e+ subs e@(ENum (FromInteger _)) = e+ subs e@(EFractional (FromRational _)) = e+ subs (ENum (Mul x y)) = (subs x) * (subs y)+ subs (ENum (Add x y)) = (subs x) + (subs y)+ subs (ENum (Sub x y)) = (subs x) - (subs y)+ subs (ENum (Negate x)) = negate (subs x)+ subs (ENum (Abs x)) = abs (subs x)+ subs (ENum (Signum x)) = signum (subs x)+ + subs (EFractional (Div x y)) = (subs x) / (subs y)+ + subs (EFloating (Pow x y)) = (subs x) ** (subs y)+ subs (EFloating (LogBase x y)) = logBase (subs x) (subs y)+ subs (EFloating (Exp x)) = exp (subs x)+ subs (EFloating (Log x)) = log (subs x)+ subs (EFloating (Sin x)) = sin (subs x)+ subs (EFloating (Cos x)) = cos (subs x)+ subs (EFloating (ASin x)) = asin (subs x)+ subs (EFloating (ATan x)) = atan (subs x)+ subs (EFloating (ACos x)) = acos (subs x)+ subs (EFloating (Sinh x)) = sinh (subs x)+ subs (EFloating (Cosh x)) = cosh (subs x)+ subs (EFloating (Tanh x)) = tanh (subs x)+ subs (EFloating (ASinh x)) = asinh (subs x)+ subs (EFloating (ATanh x)) = atanh (subs x)+ subs (EFloating (ACosh x)) = acosh (subs x) -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 (-)- abs = makeUnary Abs abs- signum = makeUnary Signum signum- fromInteger = EDimensionless . fromInteger- negate = makeUnary Neg negate+-- | this substitute is sketchy because it doesn't perform simplifications that are often assumed to be done+sketchySubstitute :: (Eq a, Hashable a, Show a) => Expr a -> [(Expr a, Expr a)] -> Expr a+sketchySubstitute expr subList+ | nonSymInputs /= [] = error $ "substitute got non-ESym input: " ++ show nonSymInputs+ | otherwise = subs expr+ where+ isSym (ESym _) = True+ isSym _ = False+ nonSymInputs = filter (not . isSym . fst) subList+ lookup' e = let hm = HM.fromList subList in+ HM.lookupDefault e e hm+ + subs e@(ESym _) = lookup' e+ subs e@(EConst _) = e+ subs e@(ENum (FromInteger _)) = e+ subs e@(EFractional (FromRational _)) = e+ subs (ENum (Mul x y)) = ENum (Mul (subs x) (subs y))+ subs (ENum (Add x y)) = ENum (Add (subs x) (subs y))+ subs (ENum (Sub x y)) = ENum (Sub (subs x) (subs y))+ subs (ENum (Negate x)) = ENum (Negate (subs x))+ subs (ENum (Abs x)) = ENum (Negate (subs x))+ subs (ENum (Signum x)) = ENum (Signum (subs x))+ + subs (EFractional (Div x y)) = EFractional (Div (subs x) (subs y))+ + subs (EFloating (Pow x y)) = EFloating (Pow (subs x) (subs y))+ subs (EFloating (LogBase x y)) = EFloating (LogBase (subs x) (subs y))+ subs (EFloating (Exp x)) = EFloating (Exp (subs x))+ subs (EFloating (Log x)) = EFloating (Log (subs x))+ subs (EFloating (Sin x)) = EFloating (Sin (subs x))+ subs (EFloating (Cos x)) = EFloating (Cos (subs x))+ subs (EFloating (ASin x)) = EFloating (ASin (subs x))+ subs (EFloating (ATan x)) = EFloating (ATan (subs x))+ subs (EFloating (ACos x)) = EFloating (ACos (subs x))+ subs (EFloating (Sinh x)) = EFloating (Sinh (subs x))+ subs (EFloating (Cosh x)) = EFloating (Cosh (subs x))+ subs (EFloating (Tanh x)) = EFloating (Tanh (subs x))+ subs (EFloating (ASinh x)) = EFloating (ASinh (subs x))+ subs (EFloating (ATanh x)) = EFloating (ATanh (subs x))+ subs (EFloating (ACosh x)) = EFloating (ACosh (subs x)) -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, Num (Vector a), LA.Container Vector a) =>- Floating (Expr sh a) where- pi = EDimensionless pi- (**) = makeBinary Pow (**)- exp = makeUnary Exp exp- log = makeUnary Log log- sin = makeUnary Sin sin- cos = makeUnary Cos cos- asin = makeUnary ASin asin- atan = makeUnary ATan atan- acos = makeUnary ACos acos- sinh = makeUnary Sinh sinh- cosh = makeUnary Cosh cosh- asinh = error "no instance for asinh"- atanh = error "no instance for atanh"- acosh = error "no instance for acosh"+---------------------------------- utility functions ------------------------------- ------------------------------- convenience functions ------------------------- -- | symbolic scalar-sym :: String -> Expr DIM0 a-sym = (ESym Z) . Sym+sym :: String -> Expr a+sym = ESym . Sym -- | Symbolic scalar which is a function of some independent variable, like time. -- . -- This lets you do d(f(g(t)))/dt == f'(g(t))*g'(t)-symDependent :: String -> Expr DIM0 a -> Expr DIM0 a+symDependent :: String -> Expr a -> Expr a symDependent name s = symDependentN name s 0 -- | same as symDependent but it can start as the Nth derivative-symDependentN :: String -> Expr DIM0 a -> Int -> Expr DIM0 a-symDependentN name (ESym _ s) n = ESym Z (SymDependent name n s)+symDependentN :: String -> Expr a -> Int -> Expr a+symDependentN name (ESym s) n = ESym (SymDependent name n s) symDependentN _ _ _ = error "symDependent got non ESym dependency" --- | symbolic dense vector-vsym :: Int -> String -> Expr DIM1 a-vsym k = (ESym (Z :. k)) . Sym --- | symbolic dense matrix-msym :: (Int,Int) -> String -> Expr DIM2 a-msym (r,c) = (ESym (Z :. r :. c)) . Sym---- | symbolic dense constant vector-vec :: Storable a => [a] -> Expr DIM1 a-vec xs = EConst $ CVec (shapeOfList [length xs]) (LA.fromList xs)+const' :: a -> Expr a+const' = EConst --- | symbolic dense constant matrix-mat :: Element a => (Int,Int) -> [[a]] -> Expr DIM2 a-mat (r,c) xs - | 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)+-- | Checks to see if an Expr is equal to a value+isVal :: Eq a => a -> Expr a -> Bool+isVal v (EConst c) = v == c+isVal v (ENum (FromInteger k)) = v == fromInteger k+isVal v (EFractional (FromRational r)) = v == fromRational r+isVal _ _ = False --- | symbolic sparse vector-svec :: String -> Int -> SparseVec (Expr DIM0 a)-svec name len = svFromList $ map (\k -> sym $ name ++ "_" ++ show k) [0..len-1]+-- | if the expression is a constant, a fromInteger, or a fromRational, return the constant part+-- otherwise return nothing+getConst :: Expr a -> Maybe a+getConst (EConst x) = Just x+getConst (ENum (FromInteger k)) = Just (fromInteger k)+getConst (EFractional (FromRational r)) = Just (fromRational r)+getConst _ = Nothing --- | symbolic sparse matrix-smat :: String -> (Int,Int) -> SparseMat (Expr DIM0 a)-smat name (rows,cols) = smFromLists allRcs+-- | Separate nonlinear and linear parts of an expression+-- @extractLinearPart (fNonLin(x)+a*x) x == (fNonLin(x), a)+extractLinearPart :: (Num a, Eq a, Show a) => Expr a -> Expr a -> (Expr a, a)+extractLinearPart e@(EConst _) _ = (e,0)+extractLinearPart e@(ENum (FromInteger _)) _ = (e,0)+extractLinearPart e@(EFractional (FromRational _)) _ = (e,0)+extractLinearPart e@(ESym _) arg+ | e == arg = (0,1)+ | otherwise = (e,0)+extractLinearPart (ENum (Add x y)) arg = (xNonlin+yNonlin, xLin+yLin) where- allRcs = map (\row -> map (\col -> (sym $ name ++ "_" ++ show row ++ "_" ++ show col)) [0..cols-1]) [0..rows-1]+ (xNonlin,xLin) = extractLinearPart x arg+ (yNonlin,yLin) = extractLinearPart y arg+extractLinearPart (ENum (Sub x y)) arg = (xNonlin-yNonlin, xLin-yLin)+ where+ (xNonlin,xLin) = extractLinearPart x arg+ (yNonlin,yLin) = extractLinearPart y arg+extractLinearPart (ENum (Negate x)) arg = (-xNonlin, -xLin)+ where+ (xNonlin,xLin) = extractLinearPart x arg+extractLinearPart e@(ENum (Mul x y)) arg = case (getConst x, getConst y) of+ (Nothing,Nothing) -> (e,0)+ (Just cx, Nothing) -> let (yNl,yL) = extractLinearPart y arg in ((EConst cx)*yNl,cx*yL)+ (Nothing, Just cy) -> let (xNl,xL) = extractLinearPart x arg in (xNl*(EConst cy),xL*cy)+ _ -> error $ "extractLinearPart got ENum (Mul x y) where x and y are both constants\n"+++ "x: " ++ show x ++ "\ny: " ++ show y+extractLinearPart e@(EFractional (Div x y)) arg = case getConst y of+ Nothing -> (e,0)+ Just cy -> let (xNl,xL) = extractLinearPart x arg in (xNl/(EConst cy),xL/cy)+extractLinearPart e@(ENum (Abs _)) _ = (e,0)+extractLinearPart e@(ENum (Signum _)) _ = (e,0)+extractLinearPart e@(EFloating (Pow _ _)) _ = (e,0)+extractLinearPart e@(EFloating (LogBase _ _)) _ = (e,0)+extractLinearPart e@(EFloating (Exp _)) _ = (e,0)+extractLinearPart e@(EFloating (Log _)) _ = (e,0)+extractLinearPart e@(EFloating (Sin _)) _ = (e,0)+extractLinearPart e@(EFloating (Cos _)) _ = (e,0)+extractLinearPart e@(EFloating (ASin _)) _ = (e,0)+extractLinearPart e@(EFloating (ATan _)) _ = (e,0)+extractLinearPart e@(EFloating (ACos _)) _ = (e,0)+extractLinearPart e@(EFloating (Sinh _)) _ = (e,0)+extractLinearPart e@(EFloating (Cosh _)) _ = (e,0)+extractLinearPart e@(EFloating (Tanh _)) _ = (e,0)+extractLinearPart e@(EFloating (ASinh _)) _ = (e,0)+extractLinearPart e@(EFloating (ATanh _)) _ = (e,0)+extractLinearPart e@(EFloating (ACosh _)) _ = (e,0) -scale :: Expr DIM0 a -> Expr sh a -> Expr sh a-scale = EScale----dot :: (Dot sh1 sh2, DotT sh1 sh2 ~ sh) => Expr sh1 a -> Expr sh2 a -> Expr sh a---dot = EDot--diff :: Expr DIM0 a -> Expr DIM0 a -> Expr DIM0 a-diff = EDeriv--grad :: Expr DIM0 a -> Expr DIM1 a -> Expr DIM1 a-grad = EGrad--jacob :: Expr DIM1 a -> Expr DIM1 a -> Expr DIM2 a-jacob = EJacob+------------------------------- arbitrary instances --------------------------+--instance Arbitrary a => Arbitrary (Nums a) where+-- arbitrary = liftM ENums arbitrary+--data Nums a = Mul a a+-- | Add a a+-- | Sub a a+-- | Negate a+-- | Abs a+-- | Signum a+-- | FromInteger Integer+ +--instance Arbitrary a => Arbitrary (Expr a) where+-- arbitrary = oneof [arbConst, arbUnary, arbBinary]+--+--arbConst :: Arbitrary a => Gen (Expr a)+--arbConst = liftM EConst arbitrary+--+--arbUnary :: Arbitrary (Expr a) => Gen (Expr a)+--arbUnary = liftM2 EUnary arbitrary arbitrary+--+--arbBinary :: Arbitrary (Expr a) => Gen (Expr a)+--arbBinary = liftM3 EBinary arbitrary arbitrary arbitrary -hess :: Expr DIM0 a -> Expr DIM1 a -> Expr DIM2 a-hess expr args = jacob (grad expr args) args+--arbNum :: (Num a, Arbitrary (Expr a)) => Gen (Expr a)+--arbNum = liftM ENum arbitrary -- arbitrary arbitrary
+ Dvda/FunGraph.hs view
@@ -0,0 +1,152 @@+{-# OPTIONS_GHC -Wall #-}+{-# Language TypeOperators #-}+{-# Language TypeFamilies #-}+{-# Language FlexibleInstances #-}++module Dvda.FunGraph ( FunGraph+ , ToFunGraph+ , NumT+ , (:*)(..)+ , MVS(..)+ , toFunGraph+ , countNodes+ , fgInputs+ , fgOutputs+ , fgLookupGExpr+ , fgReified+ , topSort+-- , fgGraph+ , nodelistToFunGraph+ ) where++import Control.Applicative+import Data.Foldable ( Foldable )+import qualified Data.Foldable as F+import qualified Data.Graph as Graph+import Data.Hashable ( Hashable )+import qualified Data.HashSet as HS+import Data.Traversable ( Traversable )+import qualified Data.Traversable as T++import Dvda.Expr+import Dvda.Reify ( ReifyGraph(..), reifyGraphs )++data FunGraph a = FunGraph { fgGraph :: Graph.Graph+ , fgInputs :: [MVS (GExpr a Int)]+ , fgOutputs :: [MVS Int]+ , fgReified :: [(Int, GExpr a Int)]+ , fgLookupGExpr :: Int -> Maybe (GExpr a Int)+ , fgVertexFromKey :: Int -> Maybe Int+ , fgNodeFromVertex :: Int -> (GExpr a Int, Int, [Int])+ }++instance Show a => Show (FunGraph a) where+ show fg = "FunGraph\ninputs:\n" ++ show (fgInputs fg) ++ "\noutputs:\n" ++ show (fgOutputs fg) ++ "\ngraph:\n" ++ show (fgGraph fg)++---- | matrix or vector or scalar+data MVS a = Mat [[a]] | Vec [a] | Sca a deriving Show++instance Functor MVS where+ fmap f (Sca x) = Sca (f x)+ fmap f (Vec xs) = Vec (map f xs)+ fmap f (Mat xs) = Mat (map (map f) xs)++instance Foldable MVS where+ foldr f x0 (Sca x) = foldr f x0 [x]+ foldr f x0 (Vec xs) = foldr f x0 xs+ foldr f x0 (Mat xs) = foldr f x0 (concat xs)++instance Traversable MVS where+ traverse f (Sca x) = Sca <$> f x+ traverse f (Vec xs) = Vec <$> T.traverse f xs+ traverse f (Mat xs) = Mat <$> T.traverse (T.traverse f) xs++class ToFunGraph a where+ type NumT a+ toMVSList :: a -> [MVS (Expr (NumT a))]+instance ToFunGraph (Expr a) where+ type NumT (Expr a) = a+ toMVSList x = [Sca x]+instance ToFunGraph [Expr a] where+ type NumT [Expr a] = NumT (Expr a)+ toMVSList x = [Vec x]+instance ToFunGraph [[Expr a]] where+ type NumT [[Expr a]] = NumT [Expr a]+ toMVSList x = [Mat x]++data a :* b = a :* b deriving Show+infixr 6 :*+instance (ToFunGraph a, ToFunGraph b, NumT a ~ NumT b) => ToFunGraph (a :* b) where+ type NumT (a :* b) = NumT a+ toMVSList (x :* y) = toMVSList x ++ toMVSList y++-- | find any symbols which are parents of outputs, but are not supplied by the user+detectMissingInputs :: (Eq a, Hashable a, Show a) => [MVS (Expr a)] -> [(Int,GExpr a Int)] -> [GExpr a Int]+detectMissingInputs exprs gr = HS.toList $ HS.difference allGraphInputs allUserInputs+ where+ allUserInputs = let f (ESym name) acc = (GSym name):acc+ f _ e = error $ "detectMissingInputs given non-ESym input \"" ++ show e ++ "\""+ in HS.fromList $ foldr f [] (concatMap F.toList exprs)++ allGraphInputs = let f (_,(GSym name)) acc = (GSym name):acc+ f _ acc = acc+ in HS.fromList $ foldr f [] gr++-- | if the same input symbol (like ESym "x") is given at two different places throw an exception+findConflictingInputs :: (Eq a, Hashable a, Show a) => [MVS (Expr a)] -> [Expr a]+findConflictingInputs exprs = HS.toList redundant+ where+ redundant = snd $ foldl f (HS.empty, HS.empty) (concatMap F.toList exprs)+ where+ f (knownExprs, redundantExprs) expr@(ESym _)+ | HS.member expr knownExprs = (knownExprs, HS.insert expr redundantExprs)+ | otherwise = (HS.insert expr knownExprs, redundantExprs)+ f _ e = error $ "findConflictingInputs saw non-ESym input \"" ++ show e ++ "\""+++-- | Take inputs and outputs which are of classes ToFunGraph (heterogenous lists of @Expr a@)+-- and traverse the outputs reifying all expressions and creating a hashmap of StableNames (stable pointers).+-- Once the hashmap is created, lookup the provided inputs and return a FunGraph which contains an+-- expression graph, input/output indices, and other useful functions. StableNames is non-deterministic+-- so this function may return graphs with more or fewer CSE's eliminated.+-- If CSE is then performed on the graph, the result is deterministic.+toFunGraph :: (Eq a, Hashable a, Show a, ToFunGraph b, ToFunGraph c, NumT b ~ a, NumT c ~ a)+ => b -> c -> IO (FunGraph a)+toFunGraph inputs outputs = mvsToFunGraph (toMVSList inputs) (toMVSList outputs)++mvsToFunGraph :: (Eq a, Hashable a, Show a) => [MVS (Expr a)] -> [MVS (Expr a)] -> IO (FunGraph a)+mvsToFunGraph inputMVSExprs outputMVSExprs = do+ -- reify the outputs+ (ReifyGraph rgr, outputMVSIndices) <- reifyGraphs outputMVSExprs+ let fg = nodelistToFunGraph rgr inputMVSGExprs outputMVSIndices+ inputMVSGExprs = map (fmap f) inputMVSExprs+ where+ f (ESym name) = (GSym name)+ f x = error $ "ERROR: mvsToFunGraph given non-ESym input \"" ++ show x ++ "\""+ return $ case (detectMissingInputs inputMVSExprs rgr, findConflictingInputs inputMVSExprs) of+ ([],[]) -> fg+ (xs,[]) -> error $ "mvsToFunGraph found inputs that were not provided by the user: " ++ show xs+ ( _,xs) -> error $ "mvsToFunGraph found idential inputs set more than once: " ++ show xs++nodelistToFunGraph :: [(Int,GExpr a Int)] -> [MVS (GExpr a Int)] -> [MVS Int] -> FunGraph a+nodelistToFunGraph rgr inputMVSIndices outputMVSIndices =+ FunGraph { fgGraph = gr+ , fgInputs = inputMVSIndices+ , fgOutputs = outputMVSIndices+ , fgLookupGExpr = lookupG+ , fgReified = rgr+ , fgVertexFromKey = lookupKey+ , fgNodeFromVertex = lookupVertex+ }+ where+ -- make sure all the inputs are symbolic, and find their indices in the Expr graph+ (gr, lookupVertex, lookupKey) = Graph.graphFromEdges $ map (\(k,gexpr) -> (gexpr, k, getParents gexpr)) rgr+ lookupG k = (\(g,_,_) -> g) <$> lookupVertex <$> lookupKey k+++---------------------------------- utilities -----------------------------+countNodes :: FunGraph a -> Int+countNodes = length . Graph.vertices . fgGraph++topSort :: FunGraph a -> [Int]+topSort fg = map ((\(_,k,_) -> k) . (fgNodeFromVertex fg)) $ Graph.topSort (fgGraph fg)
− Dvda/Graph.hs
@@ -1,217 +0,0 @@-{-# OPTIONS_GHC -Wall #-}-{-# Language StandaloneDeriving #-}-{-# Language TypeSynonymInstances #-}-{-# Language FlexibleInstances #-}-{-# Language GADTs #-}-{-# Language RankNTypes #-}--module Dvda.Graph ( FunGraph(..)- , DynamicExpr(..)- , DvdaDim(..)- , FgNode- , SymSet- , emptyFunGraph- , fgLookup- , fgExprFromKey- , insert- , previewGraph- , toFGLGraph- , collisions- , showCollisions- , funGraphSummary- , funGraphSummary'- , showNodes- , asIfExpr- ) where--import Data.Graph.Inductive ( Gr, mkGraph )-import Data.GraphViz ( Labellable, toLabelValue, preview )-import Data.GraphViz.Attributes.Complete ( Label )-import Control.Concurrent ( threadDelay )-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.IntMap as IM-import Numeric.LinearAlgebra ( Element )-import Data.Array.Repa ( Shape, DIM0, DIM1, DIM2 )-import Control.Monad.State ( State, get, put )--import Dvda.Expr ( Expr(..), Const(..), Sym(..), RefHash(..), dim )-import qualified Dvda.HashMap as HM----------------------- 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 (DynamicExpr a) (FgNode a)) -- main lookup- (IM.IntMap (DynamicExpr a)) -- internal for reverse lookup- b- c -- deriving Show- -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)--class Shape sh => DvdaDim sh where- makeDynamic :: Expr sh a -> DynamicExpr a- fromDynamic :: sh -> DynamicExpr a -> Expr sh a--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"--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) (RefHash (hash 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) (RefHash (hash 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)- ]--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, 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))- , ""- ] ++ [funGraphSummary fg]--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- numTotal = length allHashes- numCollisions = countCollisions 0 allHashes- where- countCollisions n (x:y:ys)- | x == y = countCollisions (n+1) (y:ys)- | otherwise = countCollisions n (y:ys)- countCollisions n [_] = n- countCollisions n [] = n--showCollisions :: (Hashable a, 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 outerr- where- inerr = error "must specify inputs"- outerr = error "must specify outputs"---previewGraph :: (Show a, Element a) => FunGraph a b c -> IO ()-previewGraph fungraph = do- preview $ toFGLGraph fungraph- threadDelay 10000--toFGLGraph :: FunGraph a b c -> Gr (DynamicExpr a) String-toFGLGraph (FunGraph hm _ _ _) = mkGraph lnodes ledges- where- 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/MultipleShooting/CoctaveTemplates.hs view
@@ -0,0 +1,122 @@+{-# OPTIONS_GHC -Wall #-}++module Dvda.MultipleShooting.CoctaveTemplates ( writeMexAll+ , writeSetupSource+ , writeUnstructConsts+ , writeToStruct+ , writeUnstruct+ , writePlot+ )where++import Data.Maybe ( fromMaybe )+import Data.Hashable ( Hashable )+import Data.List ( elemIndex, transpose )++import Dvda.Expr ( Expr(..), Sym(..) )+import Dvda.HashMap ( HashMap )+import qualified Dvda.HashMap as HM++writeMexAll :: String -> String+writeMexAll name = unlines $ map f ["time", "outputs", "sim", "cost", "constraints"]+ where+ f x = "tic\nfprintf('mexing " ++ file ++ "... ')\n"++"mex " ++ file ++ "\nt1 = toc;\nfprintf('finished in %.2f seconds\\n', t1)"+ where+ file = name ++ "_" ++ x ++ ".c"+++writeSetupSource :: Show a => String -> [Expr a] -> [a] -> [a] -> String+writeSetupSource name dvs lbs ubs =+ unlines $+ [ "function [x0, Aineq, bineq, Aeq, beq, lb, ub] = "++ name ++"_setup()"+ , ""+ , "x0 = zeros(" ++ show (length dvs) ++ ",1);"+ , "Aineq = [];"+ , "bineq = [];"+ , "Aeq = [];"+ , "beq = [];"+ , "lb = " ++ show lbs ++ "';"+ , "ub = " ++ show ubs ++ "';"+ ]+++-- take nice matlab structs and return vector of design constants+writeUnstructConsts :: Eq a => String -> [Expr a] -> String+writeUnstructConsts name constants =+ unlines $+ [ "function constants = " ++ name ++ "_unstructConstants(constStruct)\n"+ , "constants = zeros(" ++ show (length constants) ++ ", 1);"+ , ""+ , concatMap fromConst constants+ ]+ where+ readName e = case e of+ ESym (Sym nm) -> nm+ _ -> error "const not ESym Sym"+ fromConst e = "constants(" ++ show (1 + (fromJustErr "fromConst error" $ e `elemIndex` constants)) ++ ") = constStruct." ++ readName e ++ ";\n"+++---- take vector of design variables and vector of constants and return nice matlab struct+writeToStruct :: (Eq a, Show a, Hashable a)+ => String -> [Expr a] -> [Expr a] -> [Expr a] -> HashMap String [Expr a] -> String+writeToStruct name dvs params constants outputMap =+ 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' = (fromJustErr "toStruct error") . (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 vector of design variables+writeUnstruct :: (Eq a, Show a)+ => String+ -> [Expr a] -> [Expr a]+ -> [Expr a] -> [[Expr a]]+ -> [Expr a] -> [[Expr a]]+ -> String+writeUnstruct name dvs params states allStates actions allActions =+ unlines $+ [ "function dvs = " ++ name ++ "_unstruct(dvStruct)\n"+ , "dvs = zeros(" ++ show (length dvs) ++ ", 1);"+ , ""+ , concatMap fromParam params+ , concat $ zipWith fromXU states (transpose allStates)+ , concat $ zipWith fromXU actions (transpose allActions)+ ]+ where+ dvIdx e = fromMaybe (error $ "dvIdx error - " ++ show e ++ " is not a design variable")+ (e `elemIndex` dvs)+ fromParam e = "dvs(" ++ show (1 + dvIdx e) ++ ") = dvStruct." ++ show e ++ ";\n"+ fromXU e es =+ "dvs(" ++ show (map ((1 +) . dvIdx) es) ++ ") = dvStruct." ++ show e ++ ";\n"++writePlot :: String -> HashMap String [Expr a] -> String+writePlot name outputMap =+ 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'+++fromJustErr :: String -> Maybe a -> a+fromJustErr _ (Just x) = x+fromJustErr message Nothing = error $ "fromJustErr got Nothing, message: \"" ++ message ++ "\""
Dvda/MultipleShooting/MSCoctave.hs view
@@ -1,23 +1,26 @@ {-# OPTIONS_GHC -Wall #-}-{-# Language FlexibleContexts #-}-{-# Language TypeFamilies #-} module Dvda.MultipleShooting.MSCoctave ( msCoctave , run ) where +import qualified Control.Monad.State as State+import Data.Hashable ( Hashable ) import qualified Data.HashSet as HS-import Data.List ( zipWith6, transpose, elemIndex )-import Data.Maybe ( fromJust, catMaybes )+import Data.List ( zipWith6 )+import Data.Maybe ( fromMaybe ) -import Dvda-import Dvda.Codegen ( writeSourceFile )-import Dvda.Expr ( Expr(..), Const(..), Sym(..) )+import Dvda.AD ( rad )+import Dvda.CGen ( showMex )+import Dvda.CSE ( cse )+import Dvda.Codegen.WriteFile ( writeSourceFile )+import Dvda.Expr ( Expr(..), sym, substitute )+import Dvda.FunGraph ( (:*)(..), toFunGraph, countNodes )+import Dvda.HashMap ( HashMap ) import qualified Dvda.HashMap as HM+import Dvda.MultipleShooting.CoctaveTemplates import Dvda.MultipleShooting.MSMonad import Dvda.MultipleShooting.Types-import Dvda.OctaveSyntax ( toOctaveSource )-import Dvda.SymMonad ( rad ) {- min f(x) st:@@ -28,246 +31,253 @@ 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]+type Integrator a = [Expr Double]+ -> [Expr Double]+ -> [Expr Double]+ -> [Expr Double]+ -> ([Expr Double]+ -> [Expr Double] -> [Expr Double])+ -> Expr Double+ -> [Expr 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 ()+-- take user provided bounds and make sure they're complete+-- return functions which will lookup bounds on given state/action @ timestep, and given param+setupBounds :: (Eq a, Hashable a, Show a)+ => [(Expr a, (a,a, BCTime))]+ -> Int+ -> (Expr a -> Int -> (a,a), Expr a -> (a,a))+setupBounds userBounds nSteps = (lookupAll, lookupParam) 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)+ lookupAll x k+ | k >= nSteps = error "don't ask for bounds at timestep >= number of total timesteps"+ | otherwise = case HM.lookup (x,k) specificTimestepBounds of+ Just bnd -> bnd+ Nothing -> case HM.lookup x everyTimestepBounds of+ Just bnd -> bnd+ Nothing -> error $ "need to set bounds for \"" ++ show x ++ "\" at timestep " ++ show k - boundMap = foldr HM.union HM.empty (map stepBounds steps)+ lookupParam x = case HM.lookup x everyTimestepBounds of+ Just bnd -> bnd+ Nothing -> error $ "need to set bounds for \"" ++ show x ++ "\"" - 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+ -- bounds set at only one timestep+-- everyTimestepBounds :: HashMap (Expr a) (a,a)+ everyTimestepBounds = let+ everyTS (e,(lb,ub,ALWAYS)) = [(e,(lb,ub))]+ everyTS _ = []+ f (e,lbub) hm =+ if HM.member e hm+ then error $ "you set bounds twice for \"" ++ show e ++ "\""+ else HM.insert e lbub hm+ in foldr f HM.empty $ concatMap everyTS userBounds - (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])+ -- bounds set at specific timestep+-- specificTimestepBounds :: HashMap (Expr a, Int) (a,a)+ specificTimestepBounds = let+ specificTS (e,(lb,ub,TIMESTEP k)) = [((e,k),(lb,ub))]+ specificTS _ = []+ f (e,lbub) hm =+ if HM.member e hm+ then error $ "you set bounds twice for \"" ++ show e ++ "\""+ else HM.insert e lbub hm+ in foldr f HM.empty $ concatMap specificTS userBounds - dodeConstraints = map (Constraint (EConst (CSingleton Z 0)) EQ) $ concat $- zipWith6 odeError (init states) (init actions) (tail states) (tail actions)- (map (execDxdt userStep) [0..]) dts+vectorizeDvs :: [[a]] -> [[a]] -> [a] -> [a]+vectorizeDvs allStates allActions params = concat allStates ++ concat allActions ++ params - allConstraints = dodeConstraints ++ (concatMap stepConstraints steps) ++ periodicConstraints+msCoctave ::+ State (Step Double) b+ -> Integrator Double+ -> Int+ -> String+ -> FilePath+ -> IO ()+msCoctave userStep' odeError n funDir name = do+ let step = State.execState userStep' $+ Step { stepStates = Nothing+ , stepActions = Nothing+ , stepDxdt = Nothing+ , stepDt = Nothing+ , stepLagrangeTerm = Nothing+ , stepMayerTerm = Nothing+ , stepBounds = []+ , stepConstraints = []+ , stepParams = HS.empty+ , stepConstants = HS.empty+ , stepOutputs = HM.empty+ , stepPeriodic = HS.empty+ }+ getWithErr :: String -> (Step Double -> Maybe c) -> c+ getWithErr fieldName f = case f step of+ Nothing -> error $ "need to set " ++ fieldName+ Just ret -> ret - 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+ actions = getWithErr "actions" stepActions+ dt = getWithErr "dt" stepDt+ (states,outputs,dxdt,lagrangeState) = let+ states' = getWithErr "states" stepStates+ dxdt' = getWithErr "dxdt" stepDxdt+ outputs' = stepOutputs step+ in+ case stepLagrangeTerm step of+ Nothing -> (states',outputs',dxdt',Nothing)+ Just (lagrangeTerm,(lb,ub)) ->+ ( states' ++ [lagrangeState']+ , HM.union outputs' $ HM.fromList+ [(lagrangeStateName, lagrangeState'), (lagrangeTermName, lagrangeTerm)]+ , dxdt'++[lagrangeTerm]+ , Just (lagrangeState',(lb,ub)) )+ where+ lagrangeState' = sym lagrangeStateName+ + params = HS.toList (stepParams step)+ constants = HS.toList (stepConstants step) - -- 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"+ allStates = [[sym $ show x ++ "__" ++ show k | x <- states] | k <- [0..(n-1)]]+ allActions = [[sym $ show u ++ "__" ++ show k | u <- actions] | k <- [0..(n-1)]]+ dvs = vectorizeDvs allStates allActions params - -------------------------------------------------------------------------------------- dvs = concat states ++ concat actions ++ params+ outputMap :: HashMap String [Expr Double]+ outputMap = HM.map f outputs+ where+ f output = zipWith (subStatesActions output) allStates allActions - costFg = runFunGraph $ do- cost' <- node cost- costGrad <- rad cost' dvs- inputs_ (dvs :* constants)- outputs_ (cost' :* costGrad)+ subStatesActions f x u = substitute f (zip states x ++ zip actions u) - constraintsFg = runFunGraph $ do- cineqJacob <- mapM (flip rad dvs) cineq- ceqJacob <- mapM (flip rad dvs) ceq- inputs_ (dvs :* constants)- outputs_ (cineq :* ceq :* cineqJacob :* ceqJacob)+ subAllTimesteps :: Expr Double -> [Expr Double]+ subAllTimesteps something = zipWith (subStatesActions something) allStates allActions - timeFg = runFunGraph $ do- inputs_ (dvs :* constants)- outputs_ $ init $ scanl (+) (EConst (CSingleton Z 0)) dts+ (lbs,ubs) = unzip $ vectorizeDvs stateBounds actionBounds paramBounds+ where+ (getAllBounds,getParamBounds) = setupBounds bounds n+ stateBounds = [[getAllBounds x k | x <- states ] | k <- [0..(n-1)]]+ actionBounds = [[getAllBounds u k | u <- actions] | k <- [0..(n-1)]]+ paramBounds = [getParamBounds p | p <- params] - outputFg = runFunGraph $ do- inputs_ (dvs :* constants)- outputs_ $ HM.elems outputMap+ bounds = stepBounds step ++ lagrangeBound+ where+ lagrangeBound = case lagrangeState of+ Nothing -> []+ Just (ls,(lb,ub)) -> [(ls,(0,0,TIMESTEP 0)),(ls, (lb, ub, ALWAYS))] - simFg = runFunGraph $ do- let x' = head states- u' = head actions- dxdt' = fromJust $ stepDxdt $ head steps- inputs_ (x' :* u' :* constants)- outputs_ dxdt'+ cost = subStatesActions finalCost (last allStates) (last allActions)+ where+ finalCost = case (stepMayerTerm step, lagrangeState) of+ (Just mc, Nothing) -> mc+ (Nothing, Just (ls,_)) -> ls+ (Just mc, Just (ls,_)) -> mc + ls+ (Nothing,Nothing) -> error "need to set cost function" - 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")+ (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])+ + execDxdt x u = map (flip substitute (zip states x ++ zip actions u)) dxdt - (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+ dodeConstraints = map (Constraint 0 EQ) $ concat $+ zipWith6 odeError (init allStates) (init allActions) (tail allStates) (tail allActions)+ (repeat execDxdt) (repeat dt) - 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 ++ "';"- ]+ allConstraints = dodeConstraints ++ (concatMap (g . (fmap subAllTimesteps)) (stepConstraints step)) ++ periodicConstraints+ where+ g (Constraint [] _ _) = []+ g (Constraint _ _ []) = []+ g (Constraint (x:xs) ord (y:ys)) = Constraint x ord y : g (Constraint xs ord ys)+ + periodicConstraints = map lookup' $ HS.toList (stepPeriodic step)+ where+ lookup' x = fromMaybe (error $ "couldn't find periodic thing \"" ++ show x ++ "\" in hashmap")+ $ HM.lookup x xuMap+ xuMap = HM.fromList $ zip states (zipWith setEqual (head allStates) (last allStates )) +++ zip actions (zipWith setEqual (head allActions) (last allActions))+ where+ setEqual x y = Constraint x EQ y - -- 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")+ (costSource,costFg0,costFg) <- do+ let costGrad = rad cost dvs+ fg0 <- toFunGraph (dvs :* constants) (cost :* costGrad)+ let fg = cse fg0+ return (showMex (name ++ "_cost") fg, fg0, fg)+ + (constraintsSource,constraintsFg0,constraintsFg) <- do+ let cineqJacob = map (flip rad dvs) cineq+ ceqJacob = map (flip rad dvs) ceq+ fg0 <- toFunGraph (dvs :* constants) (cineq :* ceq :* cineqJacob :* ceqJacob)+ let fg = cse fg0+ return (showMex (name ++ "_constraints") fg, fg0, fg) + (timeSource,timeFg) <- do+ fg <- toFunGraph (dvs :* constants) (take n $ scanl (+) 0 (repeat dt))+ return (showMex (name ++ "_time") fg, fg) - -- 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"+ (outputSource,outputFg) <- do+ fg <- toFunGraph (dvs :* constants) (HM.elems outputMap)+ return (showMex (name ++ "_outputs") fg, fg) - 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"+ (simSource,simFg) <- do+ fg <- toFunGraph (states :* actions :* params :* constants) dxdt+ return (showMex (name ++ "_sim") fg, fg)+ + let setupSource = writeSetupSource name dvs lbs ubs+ mexAllSource = writeMexAll name+ unstructConstsSource = writeUnstructConsts name constants+ structSource = writeToStruct name dvs params constants outputMap+ unstructSource = writeUnstruct name dvs params states allStates actions allActions+ plotSource = writePlot name outputMap - 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"+ _ <- writeSourceFile mexAllSource funDir $ name ++ "_mex_all.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 plotSource funDir $ name ++ "_plot.m" - 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'+ _ <- writeSourceFile timeSource funDir $ name ++ "_time.c"+ _ <- writeSourceFile outputSource funDir $ name ++ "_outputs.c"+ _ <- writeSourceFile simSource funDir $ name ++ "_sim.c"+ _ <- writeSourceFile costSource funDir $ name ++ "_cost.c"+ _ <- writeSourceFile constraintsSource funDir $ name ++ "_constraints.c" + putStrLn $ "nodes in time: " ++ show (countNodes timeFg)+ putStrLn $ "nodes in output: " ++ show (countNodes outputFg)+ putStrLn $ "nodes in sim: " ++ show (countNodes simFg)+ putStrLn $ "nodes in cost: " ++ show (countNodes costFg) +++ " (" ++ show (countNodes costFg0) ++ " before CSE)"+ putStrLn $ "nodes in constraints: " ++ show (countNodes constraintsFg) +++ " (" ++ show (countNodes constraintsFg0) ++ " before CSE)"+ spring :: State (Step Double) () spring = do [x, v] <- setStates ["x","v"]- [u] <- setActions ["u"]+ [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"+ setDxdt [v, -k*x - b*v + u]+ setDt (tEnd/((fromIntegral n')-1)) + setLagrangeTerm cost (-1,2000)+ 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 u (-200, 200) ALWAYS setBound v (0,0) (TIMESTEP (n'-1)) setPeriodic x+ setPeriodic u +tEnd :: Expr Double+tEnd = 1.5+ n' :: Int-n' = 20+n' = 18 run :: IO ()-run = msCoctave spring simpsonsRuleError' n' "../Documents/MATLAB/" "cartpole"---run = msCoctave spring eulerError' n' "../Documents/MATLAB/" "cartpole"+run = msCoctave spring simpsonsRuleError' n' "../Documents/MATLAB/" "spring"+--run = msCoctave spring eulerError' n' "../Documents/MATLAB/" "spring"
Dvda/MultipleShooting/MSMonad.hs view
@@ -1,5 +1,4 @@ {-# OPTIONS_GHC -Wall #-}-{-# Language FlexibleContexts #-} module Dvda.MultipleShooting.MSMonad ( State , setStates@@ -9,34 +8,35 @@ , addConstant , addConstants , setDxdt- , setCost+ , setLagrangeTerm+ , setMayerTerm , setDt , addOutput- , getTimeStep , setPeriodic , addConstraint , setBound- , runOneStep- , execDxdt+ , lagrangeStateName+ , lagrangeTermName ) where -import Data.Array.Repa ( Z(..) ) import Data.Hashable ( Hashable ) import qualified Data.HashSet as HS-import Data.List ( nub, sort ) --, union )-import Data.Maybe ( isJust, isNothing )+import Data.List ( nub, sort )+import Data.Maybe ( isJust, fromMaybe )+import Data.Monoid ( mappend ) 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 qualified Dvda.HashMap as HM++import Dvda.Expr ( Expr(..), sym ) import Dvda.MultipleShooting.Types +lagrangeStateName,lagrangeTermName :: String+lagrangeStateName = "lagrangeState"+lagrangeTermName = "lagrangeTerm"+ failDuplicates :: [String] -> [String] failDuplicates names | length names == length (nub names) = names@@ -44,45 +44,48 @@ checkOctaveName :: String -> String checkOctaveName name- | any (`elem` "\"'~!@#$%^&*()+`-=[]{}\\|;:,.<>/?") name =- error $ "ERROR: addOutput saw illegal octave variable character in string: \"" ++ name ++ "\""+ | any (`elem` badChars) name =+ error $ "ERROR: saw illegal octave variable character in string: \"" ++ name +++ "\", illegal characters: " ++ badChars+ | name == lagrangeStateName = error "don't call your variable \"" ++ lagrangeStateName ++ "\", it's reserved"+ | name == lagrangeTermName = error "don't call your variable \"" ++ lagrangeTermName ++ "\", it's reserved" | otherwise = name+ where+ badChars = "\"'~!@#$%^&*()+`-=[]{}\\|;:,.<>/?" -setStates :: [String] -> State (Step a) [Expr Z a]+setStates :: [String] -> State (Step a) [Expr 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+ case stepStates step of Just _ -> error "states already set, don't call setStates twice"+ 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))}+ syms = map sym (failDuplicates names)+ State.put $ step {stepStates = Just syms} zipWithM_ addOutput syms names return syms -setActions :: [String] -> State (Step a) [Expr Z a]+setActions :: [String] -> State (Step a) [Expr 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+ case stepActions step of Just _ -> error "actions already set, don't call setActions twice"+ 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))}+ syms = map sym (failDuplicates names)+ State.put $ step {stepActions = Just syms} zipWithM_ addOutput syms names return syms -addParam :: (Eq (Expr Z a), Hashable (Expr Z a)) => String -> State (Step a) (Expr Z a)+addParam :: (Eq a, Hashable a) => String -> State (Step a) (Expr 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 :: (Eq a, Hashable a) => String -> State (Step a) (Expr a) addConstant name = do- [blah] <- addParams [name]+ [blah] <- addConstants [name] return blah -addParams :: (Eq (Expr Z a), Hashable (Expr Z a)) => [String] -> State (Step a) [Expr Z a]+addParams :: (Eq a, Hashable a) => [String] -> State (Step a) [Expr a] addParams names = do step <- State.get let syms = map (sym . checkOctaveName) names@@ -90,7 +93,7 @@ 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 :: (Eq a, Hashable a) => [String] -> State (Step a) [Expr a] addConstants names = do step <- State.get let syms = map (sym . checkOctaveName) names@@ -98,118 +101,55 @@ State.put $ step {stepConstants = HS.union constants0 (HS.fromList syms)} return syms -addOutput :: Expr Z a -> String -> State (Step a) ()+addOutput :: Expr 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}+ State.put $ step {stepOutputs = HM.insertWith err (checkOctaveName name) var hm} -setDt :: Expr Z a -> State (Step a) ()+setDt :: Expr 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 :: (Eq a, Hashable a, Show a) => Expr a -> State (Step a) () setPeriodic var = do step <- State.get- State.put $ step {stepPeriodic = HS.insert var (stepPeriodic step)}- + let newPeriodic+ | var `HS.member` (stepPeriodic step) = error $ "you called setPeriodic twice on \"" ++ show var ++ "\""+ | not (var `elem` (fromMaybe [] (mappend (stepStates step) (stepActions step)))) =+ error $ "you can only make states or actions periodic, you can't make \"" ++ show var ++ "\" periodic"+ | otherwise = HS.insert var (stepPeriodic step)+ State.put $ step {stepPeriodic = newPeriodic} ------------------------------------------- -setDxdt :: [Expr Z a] -> State (Step a) ()+setDxdt :: [Expr 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+setLagrangeTerm :: Expr a -> (a,a) -> State (Step a) ()+setLagrangeTerm var (lb,ub) = 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)+ when (isJust (stepLagrangeTerm step)) $ error "Lagrange term already set, don't call setLagrangeTerm twice"+ State.put $ step {stepLagrangeTerm = Just (var,(lb,ub))} - 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+setMayerTerm :: Expr a -> State (Step a) ()+setMayerTerm var = do+ step <- State.get+ when (isJust (stepMayerTerm step)) $ error "Mayer term already set, don't call setMayerTerm twice"+ State.put $ step {stepMayerTerm = Just var} - 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+setBound :: (Show a, Eq a, Hashable a)+ => Expr a -> (a, a) -> BCTime -> State (Step a) ()+setBound var@(ESym _) (lb, ub) bctime = do step <- State.get- case stepStates step of- Left _ -> error "WARNING - setBound called on non-design variable, use addConstraint instead"- _ -> return ()-+ State.put $ step {stepBounds = (var, (lb,ub,bctime)):(stepBounds step)}+setBound _ _ _ = error "WARNING - setBound called on non-design variable, try addConstraint instead" -addConstraint :: Expr Z a -> Ordering -> Expr Z a -> State (Step a) ()+addConstraint :: Expr a -> Ordering -> Expr 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
@@ -1,166 +0,0 @@-{-# 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
@@ -13,37 +13,39 @@ import Data.HashSet ( HashSet ) -import Dvda ( Z ) import Dvda.Expr ( Expr(..) ) import Dvda.HashMap ( HashMap ) 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 Constraint a = Constraint a Ordering a deriving Show+instance Functor Constraint where+ fmap f (Constraint x ordering y) = Constraint (f x) ordering (f y)+ -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 Step a = Step { stepStates :: Maybe [Expr a]+ , stepActions :: Maybe [Expr a]+ , stepParams :: HashSet (Expr a)+ , stepConstants :: HashSet (Expr a)+ , stepDxdt :: Maybe [Expr a]+ , stepLagrangeTerm :: Maybe (Expr a, (a,a))+ , stepMayerTerm :: Maybe (Expr a)+ , stepDt :: Maybe (Expr a)+ , stepBounds :: [(Expr a, (a,a, BCTime))]+ , stepConstraints :: [Constraint (Expr a)]+ , stepOutputs :: HashMap String (Expr a)+ , stepPeriodic :: HashSet (Expr a) } -data Ode a = Ode (SparseVec (Expr Z a) -> SparseVec (Expr Z a) -> SparseVec (Expr Z a)) (Int,Int)+data Ode a = Ode (SparseVec (Expr a) -> SparseVec (Expr a) -> SparseVec (Expr 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 :: Fractional (Expr a)+ => (SparseVec (Expr a) -> SparseVec (Expr a) -> SparseVec (Expr a) -> SparseVec (Expr a) -> Ode a -> Expr a -> SparseVec (Expr a))+ -> [Expr a] -> [Expr a] -> [Expr a] -> [Expr a]+ -> ([Expr a] -> [Expr a] -> [Expr a])+ -> Expr a+ -> [Expr a] wrapOdeError odeError xk uk xkp1 ukp1 dxdt dt = denseListFromSv $ odeError xk' uk' xkp1' ukp1' (Ode dxdt' (error "FUUUUCK")) dt where@@ -53,26 +55,26 @@ 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' :: Fractional (Expr a)+ => [Expr a] -> [Expr a] -> [Expr a] -> [Expr a]+ -> ([Expr a] -> [Expr a] -> [Expr a])+ -> Expr a+ -> [Expr 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' :: Fractional (Expr a)+ => [Expr a] -> [Expr a] -> [Expr a] -> [Expr a]+ -> ([Expr a] -> [Expr a] -> [Expr a])+ -> Expr a+ -> [Expr 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 :: Fractional (Expr a) => SparseVec (Expr a) -> SparseVec (Expr a) -> SparseVec (Expr a) -> SparseVec (Expr a) -> Ode a -> Expr a -> SparseVec (Expr 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 :: Fractional (Expr a) => SparseVec (Expr a) -> SparseVec (Expr a) -> SparseVec (Expr a) -> SparseVec (Expr a) -> Ode a -> Expr a -> SparseVec (Expr a) simpsonsRuleError xk uk xkp1 ukp1 (Ode ode _) dt = xkp1 - xk - (svScale (dt/6.0) (f0 + fourFm + f1)) where f0 = ode xk uk
− Dvda/OctaveSyntax.hs
@@ -1,165 +0,0 @@-{-# 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/Reify.hs view
@@ -0,0 +1,88 @@+{-# OPTIONS_GHC -Wall #-}+{-# Language RankNTypes #-}+{-# Language TemplateHaskell #-}+{-# Language TypeFamilies #-}++-- this file is a modified version from Andy Gill's data-reify package++module Dvda.Reify ( MuRef(..)+ , ReifyGraph(..)+ , reifyGraphs+ ) where++import Control.Concurrent.MVar ( newMVar, takeMVar, putMVar, MVar, readMVar )+import Control.Applicative ( Applicative )+import Data.Hashable ( Hashable, hash )+import Data.Traversable ( Traversable )+import qualified Data.Traversable as T+import System.Mem.StableName ( StableName, makeStableName, hashStableName )+import Unsafe.Coerce ( unsafeCoerce )++import Dvda.ReifyGraph ( ReifyGraph(..) )++import qualified Data.HashTable.IO as H+type HashTable k v = H.CuckooHashTable k v++class MuRef a where+ type DeRef a :: * -> *+ mapDeRef :: Applicative f+ => (forall b . (MuRef b, DeRef a ~ DeRef b) => b -> f u)+ -> a+ -> f (DeRef a u)++-- | 'reifyGraph' takes a data structure that admits 'MuRef', and returns a 'ReifyGraph' that contains+-- the dereferenced nodes, with their children as 'Int' rather than recursive values.+reifyGraphs :: (MuRef s, Traversable t) => [t s] -> IO (ReifyGraph (DeRef s), [t Int])+reifyGraphs m = do+ stableNameMap <- H.new >>= newMVar+ graph <- newMVar []+ uVar <- newMVar 0+ roots <- mapM (T.mapM (findNodes stableNameMap graph uVar)) m+ pairs <- readMVar graph+ return (ReifyGraph pairs, roots)++findNodes :: MuRef s+ => MVar (HashTable DynStableName Int)+ -> MVar [(Int,DeRef s Int)]+ -> MVar Int+ -> s+ -> IO Int+findNodes stableNameMap graph uVar j | j `seq` True = do+ st <- makeDynStableName j+ tab <- takeMVar stableNameMap+ amIHere <- H.lookup tab st+ case amIHere of+ -- if the j's StableName is already in the table, return the element+ Just var -> do putMVar stableNameMap tab+ return var+ -- if j's StableName is not yet in the table, recursively call findNodes+ Nothing -> do var <- newUnique uVar+ H.insert tab st var+ putMVar stableNameMap tab+ res <- mapDeRef (findNodes stableNameMap graph uVar) j+ tab' <- takeMVar graph+ putMVar graph $ (var,res) : tab'+ return var+findNodes _ _ _ _ = error "findNodes: strictness seq function failed to return True"++newUnique :: MVar Int -> IO Int+newUnique var = do+ v <- takeMVar var+ let v' = succ v+ putMVar var v'+ return v'+ +-- Stable names that not use phantom types.+-- As suggested by Ganesh Sittampalam.+data DynStableName = DynStableName (StableName ())++instance Hashable DynStableName where+ hash (DynStableName sn) = hashStableName sn+ +instance Eq DynStableName where+ (DynStableName sn1) == (DynStableName sn2) = sn1 == sn2++makeDynStableName :: a -> IO DynStableName+makeDynStableName a = do+ st <- makeStableName a+ return $ DynStableName (unsafeCoerce st)
+ Dvda/ReifyGraph.hs view
@@ -0,0 +1,16 @@+{-# OPTIONS_GHC -Wall #-}+{-# Language FlexibleContexts #-}+{-# Language UndecidableInstances #-}++-- this file is a modified version from Andy Gill's data-reify package++module Dvda.ReifyGraph ( ReifyGraph(..)+ ) where++data ReifyGraph e = ReifyGraph [(Unique,e Unique)]++type Unique = Int++-- | If 'e' is s Functor, and 'e' is 'Show'-able, then we can 'Show' a 'Graph'.+instance (Show (e Int)) => Show (ReifyGraph e) where+ show (ReifyGraph netlist) = show [ (u,e) | (u,e) <- netlist]
− Dvda/SymMonad.hs
@@ -1,342 +0,0 @@-{-# OPTIONS_GHC -Wall #-}-{-# Language TypeOperators #-}-{-# Language TypeFamilies #-}-{-# Language FlexibleInstances #-}-{-# Language FlexibleContexts #-}-{-# Language GADTs #-}-{-# Language DoAndIfThenElse #-}--module Dvda.SymMonad ( (:*)(..)- , MkFunGraph(..)- , node- , inputs- , inputs_- , outputs- , outputs_- , makeFunGraph- , runFunGraph- , rad- , getSensitivities- , recover- , fullShow- , fullShowNodes- , runDeriv- ) where--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.Maybe ( fromJust )-import qualified Data.HashSet as HS-import qualified Data.IntMap as IM-import Numeric.LinearAlgebra ( Element, Vector, Matrix )-import qualified Numeric.LinearAlgebra as LA--- import Debug.Trace--import Dvda.Dual ( Dual(..), dualPerturbation )-import Dvda.BinUn ( applyUnary, applyBinary )-import Dvda.Graph ( FunGraph(..), DynamicExpr(..), DvdaDim(..), insert, emptyFunGraph, fgLookup, fgExprFromKey )-import Dvda.Expr ( Expr(..), Const(..), Sym(..), dim )-import qualified Dvda.HashMap as HM------ | 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'- insert $ EUnary op x-node (EBinary op x' y') = do- x <- node x'- 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 (EGrad x_ arg_) = do- x <- node x_- arg <- node arg_- outs <- rad x [arg]- node (head outs)----- 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 args = map (\(ERef sh _ k) -> fromJust $ fgExprFromKey sh k fg) args''- argSet = HS.fromList (map makeDynamic args)-- sensitivities <- getSensitivities argSet expr (EConst (CSingleton (dim expr) 1))- -- order inputs requested by user- - 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 (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- commonDExprs = HM.keys $ HM.intersection xs ys- -- the initial union that needs conflicts fixed- union0 = xs `HM.union` ys- 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, Element a, DvdaDim sh) =>- Expr sh a -> State (FunGraph a b c) (Maybe (HS.HashSet (DynamicExpr a)))-lookupSymSet expr = do- fg <- get- 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-- 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 (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- 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- - let dfdg = dualPerturbation $ applyBinary op (Dual g 1) (Dual h 0)- dfdh = dualPerturbation $ applyBinary op (Dual g 0) (Dual h 1)- - 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 :*------------------------------------- input/output class ----------------------------------------------class MkFunGraph a where- type NumT a- type GenT a- mkNodes :: a -> State (FunGraph (NumT a) b c) a--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 GenT (a :* b) = GenT a :* GenT b- mkNodes (x :* y) = do- x' <- mkNodes x- y' <- mkNodes y- return (x' :* y')--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 exprs outs- return 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 exprs- return exprs--inputs_ :: MkFunGraph b => b -> State (FunGraph (NumT b) b c) ()-inputs_ exprs = do- _ <- inputs exprs- return ()--outputs_ :: MkFunGraph c => c -> State (FunGraph (NumT c) b c) ()-outputs_ exprs = do- _ <- outputs exprs- return ()-------------------- utility function ------------------runFunGraph :: State (FunGraph a b c) d -> FunGraph a b c-runFunGraph f = snd $ runState f emptyFunGraph--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/Vis.hs view
@@ -0,0 +1,69 @@+{-# OPTIONS_GHC -Wall #-}+{-# Language TypeOperators #-}+{-# Language TypeFamilies #-}+{-# Language FlexibleInstances #-}++module Dvda.Vis ( previewGraph+ ) where++import Control.Concurrent ( threadDelay )+import Data.GraphViz ( Labellable, toLabelValue, preview )+import Data.GraphViz.Attributes.Complete ( Label )+import qualified Data.Graph.Inductive as FGL++import Dvda.Expr+import Dvda.FunGraph++previewGraph :: Show a => FunGraph a -> IO ()+previewGraph fg = do+ preview $ toFGLGraph fg+ threadDelay 10000++toFGLGraph :: FunGraph a -> FGL.Gr (FGLNode a) (FGLEdge a)+toFGLGraph fg = FGL.mkGraph fglNodes fglEdges+ where+ fglNodes = map (\(k,gexpr) -> (k, FGLNode (k, gexpr))) $ fgReified fg+ fglEdges = concatMap nodeToEdges $ fgReified fg+ where+ nodeToEdges (k,gexpr) = map (\p -> (p,k,FGLEdge (p,k,gexpr))) (getParents gexpr)++data FGLNode a = FGLNode (Int, GExpr a Int)+data FGLEdge a = FGLEdge (Int, Int, GExpr a Int)+instance Eq (FGLEdge a) where+ (==) (FGLEdge (p0,k0,_)) (FGLEdge (p1,k1,_)) = (==) (p0,k0) (p1,k1)+instance Ord (FGLEdge a) where+ compare (FGLEdge (p0,k0,_)) (FGLEdge (p1,k1,_)) = compare (p0,k0) (p1,k1)++instance Labellable (FGLEdge a) where+ toLabelValue (FGLEdge (p,k,_)) = toLabelValue $ show p ++ " --> " ++ show k++tlv :: Int -> String -> Label+tlv k s = toLabelValue $ show k ++ ": " ++ s++instance Show a => Labellable (FGLNode a) where+ toLabelValue (FGLNode (k, (GSym s))) = tlv k (show s)+ toLabelValue (FGLNode (k, (GConst c))) = tlv k (show c)+ toLabelValue (FGLNode (k, (GNum (Mul _ _)))) = tlv k "*"+ toLabelValue (FGLNode (k, (GNum (Add _ _)))) = tlv k "+"+ toLabelValue (FGLNode (k, (GNum (Sub _ _)))) = tlv k "-"+ toLabelValue (FGLNode (k, (GNum (Negate _)))) = tlv k "-"+ toLabelValue (FGLNode (k, (GNum (Abs _)))) = tlv k "abs"+ toLabelValue (FGLNode (k, (GNum (Signum _)))) = tlv k "signum"+ toLabelValue (FGLNode (k, (GNum (FromInteger x)))) = tlv k (show x)+ toLabelValue (FGLNode (k, (GFractional (Div _ _)))) = tlv k "/"+ toLabelValue (FGLNode (k, (GFractional (FromRational x)))) = tlv k (show (fromRational x :: Double))+ toLabelValue (FGLNode (k, (GFloating (Pow _ _)))) = tlv k "**"+ toLabelValue (FGLNode (k, (GFloating (LogBase _ _)))) = tlv k "logBase"+ toLabelValue (FGLNode (k, (GFloating (Exp _)))) = tlv k "exp"+ toLabelValue (FGLNode (k, (GFloating (Log _)))) = tlv k "log"+ toLabelValue (FGLNode (k, (GFloating (Sin _)))) = tlv k "sin"+ toLabelValue (FGLNode (k, (GFloating (Cos _)))) = tlv k "cos"+ toLabelValue (FGLNode (k, (GFloating (ASin _)))) = tlv k "asin"+ toLabelValue (FGLNode (k, (GFloating (ATan _)))) = tlv k "atan"+ toLabelValue (FGLNode (k, (GFloating (ACos _)))) = tlv k "acos"+ toLabelValue (FGLNode (k, (GFloating (Sinh _)))) = tlv k "sinh"+ toLabelValue (FGLNode (k, (GFloating (Cosh _)))) = tlv k "cosh"+ toLabelValue (FGLNode (k, (GFloating (Tanh _)))) = tlv k "tanh"+ toLabelValue (FGLNode (k, (GFloating (ASinh _)))) = tlv k "asinh"+ toLabelValue (FGLNode (k, (GFloating (ATanh _)))) = tlv k "atanh"+ toLabelValue (FGLNode (k, (GFloating (ACosh _)))) = tlv k "acosh"
LICENSE view
@@ -1,4 +1,5 @@-Copyright (c) 2011, Greg Horn+Copyright (c) 2011-2012 Greg Horn+Copyright (c) 2009 Andy Gill All rights reserved.
+ TestMain.hs view
@@ -0,0 +1,18 @@+{-# OPTIONS_GHC -Wall #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts #-}++module Main where++import Test.Framework ( defaultMain ) ++import Dvda.Tests.Unary ( unaryTests )++-- 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]+main :: IO ()+main = defaultMain [unaryTests]
dvda.cabal view
@@ -1,37 +1,36 @@ Name: dvda-Version: 0.2.2+Version: 0.3 License: BSD3 License-file: LICENSE Author: Greg Horn Maintainer: gregmainland@gmail.edu+Copyright: (c) 2011 - 2012 Greg Horn+ (c) 2009 Andy Gill Stability: Experimental Category: Numerical, Math Build-type: Custom-Synopsis: Efficient automatic differentiation+Synopsis: Efficient automatic differentiation and code generation Cabal-version: >= 1.8 Description: { dvda == DVDA Verifiably Differentiates Algorithmically .-This library provides a symbolic type `Dvda.Expr` which is+This library provides a symbolic scalar type `Dvda.Expr` which is manipulated mathematically through its Num\/Fractional\/Floating instances.-Expr can be a scalar, vector, or matrix. Binary operations (adding\/multiplying\/etc)-are all elementwise. .-Matrix/vector/scalar safety is enforced at compile time-.--Efficient derivatives can be computed. Internally reverse automatic differentiation-is performed including efficient common subexpression elimination.+Automatic differentiation can be performed with `Dvda.AD`. Expressions can be turned into+computational graphs (@FunGraph@s) using toFunGraph. This uses unsafe reification for performance reasons,+and explicit common subexpression elimination using hashing can be performed using `Dvda.CSE` .-Function graphs can be JIT compiled into efficient functions using "buildHSFunction".-This is the intended way to use this library.+@FunGraph@s can be converted to C code and MATLAB mex functions. In the future there will be JIT compilation+so you can call these functions efficiently from Haskell. . Pretty graphviz plots! .-If the runtime JIT stuff works in terminal ghci but not emacs haskell-mode, you may need to add-`(setenv "PATH" (concatenate 'string (getenv "PATH") ":/usr/local/bin"))` to your .emacs file-.-To get started look in `Dvda.Examples` or CompileTest.hs in the github repo+To get started check out the source for `Dvda.Examples`+--If the runtime JIT stuff works in terminal ghci but not emacs haskell-mode, you may need to add+--`(setenv "PATH" (concatenate 'string (getenv "PATH") ":/usr/local/bin"))` to your .emacs file+--.+--To get started look in `Dvda.Examples` or CompileTest.hs in the github repo } source-repository head@@ -45,54 +44,54 @@ Library Exposed-modules: Dvda- Dvda.BinUn- Dvda.CallNative- Dvda.Codegen- Dvda.Config--- Dvda.Dot Dvda.Dual- Dvda.Examples+-- Dvda.OldExamples+-- Dvda.OctaveSyntax+-- Dvda.Tests.Function+-- Dvda.Tests.Unary+ Dvda.SparseLA++ Dvda.AD+ Dvda.CGen+ -- Dvda.Codegen.CPlugins+ Dvda.Codegen.Gcc+ Dvda.Codegen.WriteFile+ Dvda.CSE Dvda.Expr- Dvda.Graph--- Dvda.HSBuilder--- Dvda.HSSyntax+ Dvda.Examples+ Dvda.FunGraph+ Dvda.MultipleShooting.CoctaveTemplates Dvda.MultipleShooting.MSCoctave Dvda.MultipleShooting.MSMonad- Dvda.MultipleShooting.MultipleShooting Dvda.MultipleShooting.Types- Dvda.OctaveSyntax- Dvda.SparseLA- Dvda.SymMonad--- Dvda.CFunction--- Dvda.Codegen.CBuilder--- Dvda.Codegen.CCallWrapper--- Dvda.Codegen.CSyntax--- Dvda.Codegen.Utils+ Dvda.Reify+ Dvda.ReifyGraph+ Dvda.Vis Other-modules: Dvda.HashMap Build-depends: base >= 4 && < 5,+ file-location >= 0.4.4 && < 0.5, hashable >= 1.1 && < 1.2,- repa >= 3.2 && < 3.3, containers >= 0.4 && < 0.5, unordered-containers >= 0.2 && < 0.3,+ hashtables >= 1.0.1.6 && < 1.1, 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,+ QuickCheck == 2.4.*,+ test-framework-quickcheck2,+ test-framework,+ 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, - Ghc-options: -Wall--- Ghc-options: -O2 -Wall -threaded- GHC-Prof-Options: -prof -fprof-auto+ Ghc-options: -Wall -O2+ GHC-Prof-Options: -Wall -O2 -prof -fprof-auto -fprof-cafs -rtsopts+ GHC-Shared-Options: -fPIC flag test@@ -100,24 +99,24 @@ 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)--- Buildable: True--- else--- Buildable: False--- --- Main-Is: StressTest.hs--- --- Ghc-Options: -O2--- --- GHC-Prof-Options: -prof -fprof-auto+ type: exitcode-stdio-1.0+ hs-source-dirs: .+ main-is: TestMain.hs+ build-depends: base,+ dvda,+ file-location >= 0.4.4 && < 0.5,+ hashable >= 1.1 && < 1.2,+ hashtables >= 1.0.1.6 && < 1.1,+ 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,+ QuickCheck == 2.4.*,+ process >= 1.1 && < 1.2,+-- directory >= 1.1 && < 1.2,+ ad,+ test-framework-quickcheck2,+ test-framework+ ghc-options: -Wall
− test/Test.hs
@@ -1,146 +0,0 @@-{-# 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-