diff --git a/Dvda.hs b/Dvda.hs --- a/Dvda.hs +++ b/Dvda.hs @@ -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 diff --git a/Dvda/AD.hs b/Dvda/AD.hs new file mode 100644 --- /dev/null +++ b/Dvda/AD.hs @@ -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 diff --git a/Dvda/BinUn.hs b/Dvda/BinUn.hs deleted file mode 100644 --- a/Dvda/BinUn.hs +++ /dev/null @@ -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 ++" )" diff --git a/Dvda/CGen.hs b/Dvda/CGen.hs new file mode 100644 --- /dev/null +++ b/Dvda/CGen.hs @@ -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 \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 diff --git a/Dvda/CSE.hs b/Dvda/CSE.hs new file mode 100644 --- /dev/null +++ b/Dvda/CSE.hs @@ -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) + diff --git a/Dvda/CallNative.hs b/Dvda/CallNative.hs deleted file mode 100644 --- a/Dvda/CallNative.hs +++ /dev/null @@ -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 diff --git a/Dvda/Codegen.hs b/Dvda/Codegen.hs deleted file mode 100644 --- a/Dvda/Codegen.hs +++ /dev/null @@ -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 diff --git a/Dvda/Codegen/Gcc.hs b/Dvda/Codegen/Gcc.hs new file mode 100644 --- /dev/null +++ b/Dvda/Codegen/Gcc.hs @@ -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 diff --git a/Dvda/Codegen/WriteFile.hs b/Dvda/Codegen/WriteFile.hs new file mode 100644 --- /dev/null +++ b/Dvda/Codegen/WriteFile.hs @@ -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 diff --git a/Dvda/Config.hs b/Dvda/Config.hs deleted file mode 100644 --- a/Dvda/Config.hs +++ /dev/null @@ -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 - diff --git a/Dvda/Examples.hs b/Dvda/Examples.hs --- a/Dvda/Examples.hs +++ b/Dvda/Examples.hs @@ -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 diff --git a/Dvda/Expr.hs b/Dvda/Expr.hs --- a/Dvda/Expr.hs +++ b/Dvda/Expr.hs @@ -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 diff --git a/Dvda/FunGraph.hs b/Dvda/FunGraph.hs new file mode 100644 --- /dev/null +++ b/Dvda/FunGraph.hs @@ -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) diff --git a/Dvda/Graph.hs b/Dvda/Graph.hs deleted file mode 100644 --- a/Dvda/Graph.hs +++ /dev/null @@ -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" diff --git a/Dvda/MultipleShooting/CoctaveTemplates.hs b/Dvda/MultipleShooting/CoctaveTemplates.hs new file mode 100644 --- /dev/null +++ b/Dvda/MultipleShooting/CoctaveTemplates.hs @@ -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 ++ "\"" diff --git a/Dvda/MultipleShooting/MSCoctave.hs b/Dvda/MultipleShooting/MSCoctave.hs --- a/Dvda/MultipleShooting/MSCoctave.hs +++ b/Dvda/MultipleShooting/MSCoctave.hs @@ -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" diff --git a/Dvda/MultipleShooting/MSMonad.hs b/Dvda/MultipleShooting/MSMonad.hs --- a/Dvda/MultipleShooting/MSMonad.hs +++ b/Dvda/MultipleShooting/MSMonad.hs @@ -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 - } diff --git a/Dvda/MultipleShooting/MultipleShooting.hs b/Dvda/MultipleShooting/MultipleShooting.hs deleted file mode 100644 --- a/Dvda/MultipleShooting/MultipleShooting.hs +++ /dev/null @@ -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) diff --git a/Dvda/MultipleShooting/Types.hs b/Dvda/MultipleShooting/Types.hs --- a/Dvda/MultipleShooting/Types.hs +++ b/Dvda/MultipleShooting/Types.hs @@ -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 diff --git a/Dvda/OctaveSyntax.hs b/Dvda/OctaveSyntax.hs deleted file mode 100644 --- a/Dvda/OctaveSyntax.hs +++ /dev/null @@ -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 diff --git a/Dvda/Reify.hs b/Dvda/Reify.hs new file mode 100644 --- /dev/null +++ b/Dvda/Reify.hs @@ -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) diff --git a/Dvda/ReifyGraph.hs b/Dvda/ReifyGraph.hs new file mode 100644 --- /dev/null +++ b/Dvda/ReifyGraph.hs @@ -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] diff --git a/Dvda/SymMonad.hs b/Dvda/SymMonad.hs deleted file mode 100644 --- a/Dvda/SymMonad.hs +++ /dev/null @@ -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 diff --git a/Dvda/Vis.hs b/Dvda/Vis.hs new file mode 100644 --- /dev/null +++ b/Dvda/Vis.hs @@ -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" diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,5 @@ -Copyright (c) 2011, Greg Horn +Copyright (c) 2011-2012 Greg Horn +Copyright (c) 2009 Andy Gill All rights reserved. diff --git a/TestMain.hs b/TestMain.hs new file mode 100644 --- /dev/null +++ b/TestMain.hs @@ -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] diff --git a/dvda.cabal b/dvda.cabal --- a/dvda.cabal +++ b/dvda.cabal @@ -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 diff --git a/test/Test.hs b/test/Test.hs deleted file mode 100644 --- a/test/Test.hs +++ /dev/null @@ -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 -