text:
{-# LANGUAGE ExistentialQuantification #-}

module Lib
    ( libMain
    ) where

import Control.Monad
import Control.Monad.Except
import Data.Array
import Data.Complex
import Data.Ratio
import Data.Void
import Numeric
import System.Environment
import Text.Megaparsec
import Text.Megaparsec.Char
import qualified Text.Megaparsec.Char.Lexer as L

-- Types
-------------------------------------------------------------

data LispVal = Atom String
             | List [LispVal]
             | DottedList [LispVal] LispVal
             | Vector (Array Int LispVal)
             | Number Integer
             | String String
             | Bool Bool
             | Character Char
             | Float Double
             | Ratio Rational
             | Complex (Complex Double)

data LispError = NumArgs Integer [LispVal]
               | TypeMismatch String LispVal
               | Parser String -- we'll just pass the ouput from errorBundlePretty here
               | BadSpecialForm String LispVal
               | NotFunction String String
               | UnboundVar String String
               | Default String

type ThrowsError = Either LispError

-- Parser
-------------------------------------------------------------

-- Numeric conversion helpers
oct2dig x = fst $ head (readOct x)
hex2dig x = fst $ head (readHex x)
bin2dig  = bin2dig' 0
bin2dig' digint "" = digint
bin2dig' digint (x:xs) = let old = 2 * digint + (if x == '0' then 0 else 1) in
                         bin2dig' old xs

toDouble :: LispVal -> Double
toDouble (Float f)  = realToFrac f
toDouble (Number n) = fromIntegral n

type Parser = Parsec Void String

symbol :: Parser Char
symbol = oneOf "!$%&|*+-/:<=>?@^_~"

spaces :: Parser ()
spaces = space1

escapedChars :: Parser Char
escapedChars = do
    char '\\'
    x <- oneOf "\\\"nrt"
    return $ case x of
        '\\' -> x
        '"'  -> x
        'n'  -> '\n'
        'r'  -> '\r'
        't'  -> '\t'
        x    -> x

parseString :: Parser LispVal
parseString = do
    char '"'
    x <- many $ escapedChars <|> noneOf "\"\\"
    char '"'
    return $ String x

-- parses characters specified like #\a #\A #\space #\newline
parseCharacter :: Parser LispVal
parseCharacter = do
    try $ string "#\\"
    value <- try (string "newline" <|> string "space")
            <|> do { x <- L.charLiteral; notFollowedBy alphaNumChar; return [x] }
    return $ Character $ case value of
        "space"     -> ' '
        "newline"   -> '\n'
        _           -> head value

parseAtom :: Parser LispVal
parseAtom = do
    first <- letterChar <|> symbol
    rest <- many (letterChar <|> digitChar <|> symbol)
    let atom = first:rest
    return $ Atom atom

parseBool :: Parser LispVal
parseBool = do
    char '#'
    (Bool True <$ char 't') <|> (Bool False <$ char 'f')

parseDecimal1 :: Parser LispVal
parseDecimal1 = Number . read <$> some digitChar

parseDecimal2 :: Parser LispVal
parseDecimal2 = do
    try $ string "#x"
    Number . read <$> some digitChar

parseHex :: Parser LispVal
parseHex = do
    try $ string "#x"
    Number . hex2dig <$> some hexDigitChar

parseOct :: Parser LispVal
parseOct = do
    try $ string "#o"
    Number . oct2dig <$> some octDigitChar

parseBin :: Parser LispVal
parseBin = do
    try $ string "#b"
    Number . bin2dig <$> some binDigitChar

parseNumber :: Parser LispVal
parseNumber = parseDecimal1 <|> parseDecimal2 <|> parseHex <|> parseOct <|> parseBin

parseFloat :: Parser LispVal
parseFloat = Float <$> L.float

parseRatio :: Parser LispVal
parseRatio = do
    x <- some digitChar
    char '/'
    y <- some digitChar
    return $ Ratio (read x % read y)

parseComplex :: Parser LispVal
parseComplex = do
    x <- try parseFloat <|> try parseDecimal1
    char '+'
    y <- try parseFloat <|> try parseDecimal1
    char 'i'
    return $ Complex (toDouble x :+ toDouble y)

parseList :: Parser LispVal
parseList = List <$> sepBy parseExpr spaces

parseDottedList :: Parser LispVal
parseDottedList = do
    head <- endBy parseExpr spaces
    tail <- char '.' >> spaces >> parseExpr
    return $ DottedList head tail

parseQuoted :: Parser LispVal
parseQuoted = do
    char '\''
    x <- parseExpr
    return $ List [Atom "quote", x]

parseQuasiQuoted :: Parser LispVal
parseQuasiQuoted = do
    char '`'
    x <- parseExpr
    return $ List [Atom "quasiquote", x]

parseUnQuote :: Parser LispVal
parseUnQuote = do
    char ','
    x <- parseExpr
    return $ List [Atom "unquote", x]

parseUnQuoteSplicing :: Parser LispVal
parseUnQuoteSplicing = do
    char ','
    char '@'
    x <- parseExpr
    return $ List [Atom "unquote-splicing", x]

parseVector :: Parser LispVal
parseVector = do
    arrayValues <- sepBy parseExpr spaces
    return $ Vector (listArray (0, length arrayValues - 1) arrayValues)

parseExpr :: Parser LispVal
parseExpr = parseAtom
         <|> parseString
         <|> try parseComplex
         <|> try parseFloat
         <|> try parseRatio
         <|> try parseNumber
         <|> try parseBool
         <|> try parseCharacter
         <|> parseQuoted
         <|> parseQuasiQuoted
         <|> parseUnQuote
         <|> parseUnQuoteSplicing
         <|> try (do
             string "#("
             x <- parseVector
             char ')'
             return x)
         <|> do
                char '('
                x <- try parseList <|> parseDottedList
                char ')'
                return x

readExpr :: String -> ThrowsError LispVal
readExpr input = case parse parseExpr "lisp" input of
    Left err -> throwError $ Parser $ errorBundlePretty err
    Right val -> return val

-- Eval
-------------------------------------------------------------

unwordsList :: [LispVal] -> String
unwordsList = unwords . map showVal

showVal :: LispVal -> String
showVal (String contents) = "\"" ++ contents ++ "\""
showVal (Atom name) = name
showVal (Number contents) = show contents
showVal (Bool True) = "#t"
showVal (Bool False) = "#f"
showVal (List contents) = "(" ++ unwordsList contents ++ ")"
showVal (DottedList head tail) = "(" ++ unwordsList head ++ " . " ++ showVal tail ++ ")"
showVal (Float contents) = show contents
showVal (Character char) = show char
showVal (Ratio contents) = show contents
showVal (Complex contents) = show contents
--TODO: Implement the rest so this isn't required
-- showVal _ = error "showVal not implemented for type"

instance Show LispVal where show = showVal

showError :: LispError -> String
showError (UnboundVar message varname)  = message ++ ": " ++ varname
showError (BadSpecialForm message form) = message ++ ": " ++ show form
showError (NotFunction message func)    = message ++ ": " ++ show func
showError (NumArgs expected found)      = "Expected " ++ show expected ++ " args; found values " ++ unwordsList found
showError (TypeMismatch expected found) = "Invalid type: expected " ++ expected ++ ", found " ++ show found
showError (Parser parseErr)             = "Parse error at " ++ show parseErr
showError (Default message)             = "Error: " ++ message

instance Show LispError where show = showError

trapError action = catchError action (return . show)

-- Non exhaustiveness here is intentional, it's a programmer error to apply extractValue to a Left val
extractValue :: ThrowsError a -> a
extractValue (Right val) = val

unaryOp :: (LispVal -> LispVal) -> [LispVal] -> ThrowsError LispVal
unaryOp f [] = throwError $ NumArgs 1 []
unaryOp f [v] = return $ f v
unaryOp _ args = throwError $ NumArgs 1 args

--TODODB: Should implement numberp for Complex
symbolp, numberp, stringp, boolp, listp :: LispVal -> LispVal
symbolp (Atom _)   = Bool True
symbolp _          = Bool False
numberp (Number _) = Bool True
numberp _          = Bool False
stringp (String _) = Bool True
stringp _          = Bool False
boolp   (Bool _)   = Bool True
boolp   _          = Bool False
listp   (List _)   = Bool True
listp   (DottedList _ _) = Bool False
listp   _          = Bool False

unpackNum :: LispVal -> ThrowsError Integer
unpackNum (Number n) = return n
unpackNum (String n) = let parsed = reads n in
                           if null parsed
                             then throwError $ TypeMismatch "number" $ String n
                             else return $ fst $ head parsed
unpackNum (List [n]) = unpackNum n
unpackNum notNum     = throwError $ TypeMismatch "number" notNum

unpackStr :: LispVal -> ThrowsError String
unpackStr (String s) = return s
unpackStr (Number s) = return $ show s
unpackStr (Bool s)   = return $ show s
unpackStr notString  = throwError $ TypeMismatch "string" notString

unpackBool :: LispVal -> ThrowsError Bool
unpackBool (Bool b) = return b
unpackBool notBool  = throwError $ TypeMismatch "boolean" notBool

numericBinop :: (Integer -> Integer -> Integer) -> [LispVal] -> ThrowsError LispVal
numericBinop op []              = throwError $ NumArgs 2 []
numericBinop op singleVal@[_]   = throwError $ NumArgs 2 singleVal
numericBinop op params          = mapM unpackNum params >>= return . Number . foldl1 op

boolBinop :: (LispVal -> ThrowsError a) -> (a -> a -> Bool) -> [LispVal] -> ThrowsError LispVal
boolBinop unpacker op args = if length args /= 2
                             then throwError $ NumArgs 2 args
                             else do left <- unpacker $ head args
                                     right <- unpacker $ args !! 1
                                     return $ Bool $ left `op` right

numBoolBinop  = boolBinop unpackNum
strBoolBinop  = boolBinop unpackStr
boolBoolBinop = boolBinop unpackBool

primitives :: [(String, [LispVal] -> ThrowsError LispVal)]
primitives = [  ("+", numericBinop (+)),
                ("-", numericBinop (-)),
                ("*", numericBinop (*)),
                ("/", numericBinop div),
                ("mod", numericBinop mod),
                ("quotient", numericBinop quot),
                ("remainder", numericBinop rem),
                ("symbol?", unaryOp symbolp),
                ("string?", unaryOp stringp),
                ("number?", unaryOp numberp),
                ("bool?", unaryOp boolp),
                ("list?", unaryOp listp),
                ("=", numBoolBinop (==)),
                ("<", numBoolBinop (<)),
                (">", numBoolBinop (>)),
                ("/=", numBoolBinop (/=)),
                (">=", numBoolBinop (>=)),
                ("<=", numBoolBinop (<=)),
                ("&&", boolBoolBinop (&&)),
                ("||", boolBoolBinop (||)),
                ("string=?", strBoolBinop (==)),
                ("string<?", strBoolBinop (<)),
                ("string>?", strBoolBinop (>)),
                ("string<=?", strBoolBinop (<=)),
                ("string>=?", strBoolBinop (>=)),
                ("car", car),
                ("cdr", cdr),
                ("cons", cons),
                ("eq?", eqv),
                ("eqv?", eqv),
                ("equal?", equal)]

-- lookup primitive, if not present, error
apply :: String -> [LispVal] -> ThrowsError LispVal
apply func args = maybe (throwError $ NotFunction "Unrecognized primitive function args" func)
                  ($ args)
                  (lookup func primitives)

eval :: LispVal -> ThrowsError LispVal
eval val@(String _) = return val
eval val@(Number _) = return val
eval val@(Bool _)   = return val
eval val@(Float _)  = return val
eval val@(Complex _) = return val
eval (List [Atom "quote", val]) = return val
eval (List [Atom "if", pred, conseq, alt]) = do
    result <- eval pred
    case result of
        Bool False -> eval alt
        Bool True -> eval conseq
        _ -> throwError $ TypeMismatch "bool" pred
eval (List (Atom func : args)) = mapM eval args >>= apply func
eval badForm = throwError $ BadSpecialForm "Unrecognized special form" badForm

-- List primitives

{-
(car '(a b c)) = a
(car '(a)) = a
(car '(a b . c)) = a
(car 'a) = error – not a list
(car 'a 'b) = error – car only takes one argument
-}
car :: [LispVal] -> ThrowsError LispVal
car [List (x : xs)]         = return x
car [DottedList (x : xs) _] = return x
car [badArg]                = throwError $ TypeMismatch "pair" badArg
car badArgList              = throwError $ NumArgs 1 badArgList

{-
(cdr '(a b c)) = (b c)
(cdr '(a b)) = (b)
(cdr '(a)) = NIL
(cdr '(a . b)) = b
(cdr '(a b . c)) = (b . c)
(cdr 'a) = error – not a list
(cdr 'a 'b) = error – too many arguments
-}
cdr :: [LispVal] -> ThrowsError LispVal
cdr [List (x : xs)]         = return $ List xs
cdr [DottedList [_] x]      = return x
cdr [DottedList (_ : xs) x] = return $ DottedList xs x
cdr [badArg]                = throwError $ TypeMismatch "pair" badArg
cdr badArgList              = throwError $ NumArgs 1 badArgList

cons :: [LispVal] -> ThrowsError LispVal
cons [x1, List []]              = return $ List [x1]
cons [x, List xs]               = return $ List $ x : xs
cons [x, DottedList xs xlast]   = return $ DottedList (x:xs) xlast
cons [x1, x2]                   = return $ DottedList [x1] x2
cons badArgList                 = throwError $ NumArgs 2 badArgList

-- There are 3 levels of equivalence predicates in scheme. eq? eqv? and equal?
-- eq? and eqv? as basically the same: Two items are the same if they print the same
-- We can write one function for both, and register just the one function for both symbols
eqv :: [LispVal] -> ThrowsError LispVal
eqv [Bool arg1, Bool arg2]             = return $ Bool $ arg1 == arg2
eqv [Number arg1, Number arg2]         = return $ Bool $ arg1 == arg2
eqv [String arg1, String arg2]         = return $ Bool $ arg1 == arg2
eqv [Atom arg1, Atom arg2]             = return $ Bool $ arg1 == arg2
eqv [DottedList xs x, DottedList ys y] = eqv [List $ xs ++ [x], List $ ys ++ [y]]
eqv [List arg1, List arg2]             = return $ Bool $ (length arg1 == length arg2) &&
                                                             all eqvPair (zip arg1 arg2)
     where eqvPair (x1, x2) = case eqv [x1, x2] of
                                Left err -> False
                                Right (Bool val) -> val
eqv [_, _]                                 = return $ Bool False
eqv badArgList                             = throwError $ NumArgs 2 badArgList

-- Now we want an equal? function that ignores differences in the type tags.
-- e.g. we want to keep weak typing and let that work across different type tags.
--      (equal? 2 "2") = #t
--
-- Existential types to the rescue.
-- We need a data type that can hold any function from a `LispVal -> something`
-- (provided that something supports equality)
data Unpacker = forall a. Eq a => AnyUnpacker (LispVal -> ThrowsError a)

-- "For any type that is an instance of Eq, you can define an Unpacker that takes a function
-- from LispVal to that type, and may throw an error."
-- All we have to do is wrap our functions with the AnyUnpacker constructor,
-- then we can create a list of Unpackers that does what we want.

-- Take an Unpacker and determine if two LispVals are equal after unpack
unpackEquals :: LispVal -> LispVal -> Unpacker -> ThrowsError Bool
unpackEquals arg1 arg2 (AnyUnpacker unpacker) =
             do unpacked1 <- unpacker arg1
                unpacked2 <- unpacker arg2
                return $ unpacked1 == unpacked2
        `catchError` const (return False)

equal :: [LispVal] -> ThrowsError LispVal
equal [arg1, arg2] = do
    let unpackers = [AnyUnpacker unpackNum, AnyUnpacker unpackStr, AnyUnpacker unpackBool]
    primitiveEquals <- liftM or $ mapM (unpackEquals arg1 arg2) unpackers
    eqvEquals <- eqv [arg1, arg2]
    return $ Bool (primitiveEquals || let (Bool x) = eqvEquals in x)
equal badArgList = throwError $ NumArgs 2 badArgList

libMain :: IO ()
libMain = do
    args <- getArgs
    let evaled = fmap show $ readExpr (args !! 0) >>= eval
    putStrLn $ extractValue $ trapError evaled


text:
module Mnist (Image(..), normalisedData, readImages, readLabels, toMatrix)
  where

import qualified Data.ByteString.Lazy as BL
import Data.Binary.Get
import Data.Word
import qualified Data.List.Split as S
import Numeric.LinearAlgebra

data Image = Image {
      iRows :: Int
    , iColumns :: Int
    , iPixels :: [Word8]
    } deriving (Eq, Show)

toMatrix :: Image -> Matrix Double
toMatrix image = (r><c) p :: Matrix Double
  where r = iRows image
        c = iColumns image
        p = map fromIntegral (iPixels image)

{-
toColumnVector :: Image -> Matrix Double
toColumnVector i = (r><1) q :: Matrix Double
  where r = Mnist.rows i * Mnist.columns i
        p = map fromIntegral (pixels i)
        q = map normalise p
-}

normalisedData :: Image -> [Double]
normalisedData image = map normalisePixel (iPixels image)

--normalisedData :: Image -> [Double]
--normalisedData i = map (/m) x
--    where x = map normalisePixel (pixels i)
--          m = sqrt( sum (zipWith (*) x x))

normalisePixel :: Word8 -> Double
normalisePixel p = (fromIntegral p) / 255.0

-- MNIST label file format
--
-- [offset] [type]          [value]          [description]
-- 0000     32 bit integer  0x00000801(2049) magic number (MSB first)
-- 0004     32 bit integer  10000            number of items
-- 0008     unsigned byte   ??               label
-- 0009     unsigned byte   ??               label
-- ........
-- xxxx     unsigned byte   ??               label
--
-- The labels values are 0 to 9.

deserialiseLabels :: Get (Word32, Word32, [Word8])
deserialiseLabels = do
  magicNumber <- getWord32be
  count <- getWord32be
  labelData <- getRemainingLazyByteString
  let labels = BL.unpack labelData
  return (magicNumber, count, labels)

readLabels :: FilePath -> IO [Int]
readLabels filename = do
  content <- BL.readFile filename
  let (_, _, labels) = runGet deserialiseLabels content
  return (map fromIntegral labels)

-- MNIST Image file format
--
-- [offset] [type]          [value]          [description]
-- 0000     32 bit integer  0x00000803(2051) magic number
-- 0004     32 bit integer  ??               number of images
-- 0008     32 bit integer  28               number of rows
-- 0012     32 bit integer  28               number of columns
-- 0016     unsigned byte   ??               pixel
-- 0017     unsigned byte   ??               pixel
-- ........
-- xxxx     unsigned byte   ??               pixel
--
-- Pixels are organized row-wise. Pixel values are 0 to 255. 0 means background (white), 255
-- means foreground (black).

deserialiseHeader :: Get (Word32, Word32, Word32, Word32, [[Word8]])
deserialiseHeader = do
  magicNumber <- getWord32be
  imageCount <- getWord32be
  r <- getWord32be
  c <- getWord32be
  packedData <- getRemainingLazyByteString
  let len = fromIntegral (r * c)
  let unpackedData = S.chunksOf len (BL.unpack packedData)
  return (magicNumber, imageCount, r, c, unpackedData)

readImages :: FilePath -> IO [Image]
readImages filename = do
  content <- BL.readFile filename
  let (_, _, r, c, unpackedData) = runGet deserialiseHeader content
  return (map (Image (fromIntegral r) (fromIntegral c)) unpackedData)

ext: hs

max_stars_repo_name: mhwombat/backprop-example

max_issues_repo_name: mhwombat/backprop-example

max_forks_repo_name: mhwombat/backprop-example


text:
{-# LANGUAGE TypeApplications #-}

module SSSFNProps where

import Data.VectorSpace
import qualified GHC.Exts as GHC
import Numeric.LinearAlgebra hiding ((#>), (===), inv)
import SSSFN
import Test.QuickCheck

--------------------------------------------------------------------------------

infix 4 ~~~

(~~~) ::
  (GHC.IsList as, Show as, a ~ GHC.Item as, RealFloat a, Show a) =>
  as ->
  as ->
  Property
xs ~~~ ys =
  let xs' = GHC.toList xs
      ys' = GHC.toList ys
      sc = 1 `max` maximum (map abs xs') `max` maximum (map abs ys')
      di = maximum (zipWith absdiff xs' ys')
   in counterexample
        ( show xs
            ++ " ~~~ "
            ++ show ys
            ++ " (difference "
            ++ show di
            ++ ", scale "
            ++ show sc
            ++ ")"
        )
        (length xs' == length ys' && di <= eps * sc)
  where
    absdiff x y = abs (x - y)
    eps = 100 * peps

--------------------------------------------------------------------------------

prop_Grid_add_assoc :: Grid Double -> Grid Double -> Grid Double -> Property
prop_Grid_add_assoc x y z = (x ^+^ y) ^+^ z ~~~ x ^+^ (y ^+^ z)

prop_Grid_add_lunit :: Grid Double -> Property
prop_Grid_add_lunit x = zeroV ^+^ x ~~~ x

prop_Grid_add_runit :: Grid Double -> Property
prop_Grid_add_runit x = x ~~~ x ^+^ zeroV

prop_Grid_add_comm :: Grid Double -> Grid Double -> Property
prop_Grid_add_comm x y = x ^+^ y ~~~ y ^+^ x

prop_Grid_add_linv :: Grid Double -> Property
prop_Grid_add_linv x = negateV x ^+^ x ~~~ zeroV

prop_Grid_add_rinv :: Grid Double -> Property
prop_Grid_add_rinv x = x ^+^ negateV x ~~~ zeroV

prop_Grid_add_scale_dist :: Double -> Grid Double -> Grid Double -> Property
prop_Grid_add_scale_dist a x y = a *^ (x ^+^ y) ~~~ a *^ x ^+^ a *^ y

prop_Grid_scale_assoc :: Double -> Double -> Grid Double -> Property
prop_Grid_scale_assoc a b x = a *^ (b *^ x) ~~~ (a * b) *^ x

prop_Grid_scale_lunit :: Grid Double -> Property
prop_Grid_scale_lunit x = 1 *^ x ~~~ x

prop_Grid_scale_szero :: Grid Double -> Property
prop_Grid_scale_szero x = 0 *^ x ~~~ zeroV

prop_Grid_scale_zero :: Double -> Property
prop_Grid_scale_zero a = a *^ (zeroV :: Grid Double) ~~~ zeroV

prop_Grid_scale_comm :: Double -> Double -> Grid Double -> Property
prop_Grid_scale_comm a b x = a *^ (b *^ x) ~~~ b *^ (a *^ x)

prop_Grid_scale_Grid_add_dist :: Double -> Double -> Grid Double -> Property
prop_Grid_scale_Grid_add_dist a b x = (a + b) *^ x ~~~ a *^ x ^+^ b *^ x

--------------------------------------------------------------------------------

prop_Op_add_assoc :: Op Double -> Op Double -> Op Double -> Property
prop_Op_add_assoc x y z = (x + y) + z ~~~ x + (y + z)

prop_Op_add_lunit :: Op Double -> Property
prop_Op_add_lunit x = 0 + x ~~~ x

prop_Op_add_runit :: Op Double -> Property
prop_Op_add_runit x = x ~~~ x + 0

prop_Op_add_comm :: Op Double -> Op Double -> Property
prop_Op_add_comm x y = x + y ~~~ y + x

prop_Op_add_linv :: Op Double -> Property
prop_Op_add_linv x = negate x + x ~~~ 0

prop_Op_add_rinv :: Op Double -> Property
prop_Op_add_rinv x = x + negate x ~~~ 0

prop_Op_add_scale_dist :: Double -> Op Double -> Op Double -> Property
prop_Op_add_scale_dist a x y = a *^ (x + y) ~~~ a *^ x + a *^ y

prop_Op_mul_assoc :: Op Double -> Op Double -> Op Double -> Property
prop_Op_mul_assoc x y z = (x * y) * z ~~~ x * (y * z)

prop_Op_mul_lunit :: Op Double -> Property
prop_Op_mul_lunit x = 1 * x ~~~ x

prop_Op_mul_runit :: Op Double -> Property
prop_Op_mul_runit x = x ~~~ x * 1

prop_Op_mul_scale_assoc :: Double -> Op Double -> Op Double -> Property
prop_Op_mul_scale_assoc a x y = a *^ (x * y) ~~~ (a *^ x) * y

prop_Op_scale_assoc :: Double -> Double -> Op Double -> Property
prop_Op_scale_assoc a b x = a *^ (b *^ x) ~~~ (a * b) *^ x

prop_Op_scale_lunit :: Op Double -> Property
prop_Op_scale_lunit x = 1 *^ x ~~~ x

prop_Op_scale_lzero :: Op Double -> Property
prop_Op_scale_lzero x = 0 *^ x ~~~ 0

prop_Op_scale_rzero :: Double -> Property
prop_Op_scale_rzero a = a *^ (0 :: Op Double) ~~~ 0

prop_Op_scale_comm :: Double -> Double -> Op Double -> Property
prop_Op_scale_comm a b x = a *^ (b *^ x) ~~~ b *^ (a *^ x)

prop_Op_scale_add_dist :: Double -> Double -> Op Double -> Property
prop_Op_scale_add_dist a b x = (a + b) *^ x ~~~ (a *^ x) + (b *^ x)

max_stars_repo_name: eschnett/sssfn

max_issues_repo_name: eschnett/sssfn

max_forks_repo_name: eschnett/sssfn


text:
import Test.HUnit
import FFT.DFTData
import FFT.Orig
import FFT.Samples
import System.Environment
import Data.Complex

testDftInnerMap = TestCase (assertEqual "DFTInnerMap" expected actual)
    where
        expected = (sum $ dftInnerMap samples512)
        actual = (sum $ dft samples512)
        samples512 = samples 1 512

testDftFoldMap =  TestCase (assertEqual "DFTFoldMap" expected actual)
    where
        expected = (sum $ dftFoldMap samples512)
        actual = (sum $ dft samples512)
        samples512 = samples 1 512

testDftOuterMap = TestCase (assertEqual "DFTOuterMap" expected actual)
    where
        expected = (sum $ dftOuterMap samples512)
        actual = (sum $ dft samples512)
        samples512 = samples 1 512

testDftInnerMapReduce = TestCase (assertEqual "DFTInnerMapReduce" expected actual)
    where
        expected = (sum $ dft samples512)
        actual = (sum $ dft samples512)
        samples512 = samples 1 512

tests = TestList [TestLabel "testDftOuterMap" testDftOuterMap,
                    TestLabel "testDftInnerMap" testDftInnerMap,
                    TestLabel "testDftFoldMap" testDftFoldMap,
                    TestLabel "testDftInnerMapReduce" testDftInnerMapReduce
                    ]

main = runTestTT tests

max_stars_repo_name: BNJHope/parallel-fourier-transforms

max_issues_repo_name: BNJHope/parallel-fourier-transforms

max_forks_repo_name: BNJHope/parallel-fourier-transforms


text:
-- Building blocks for making fully connected neural networks (FCNs).
--
-- Original author: David Banas <capn.freako@gmail.com>
-- Original date:   January 18, 2018
--
-- Copyright (c) 2018 David Banas; all rights reserved World wide.

{-# OPTIONS_GHC -Wall #-}
{-# OPTIONS_GHC -Wno-unused-top-binds #-}

{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE ExplicitForAll #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE TypeOperators #-}

{-|
Module      : Haskell_ML.FCN
Description : Allows: creation, training, running, saving, and loading,
              of multi-layer, fully connected neural networks.
Copyright   : (c) David Banas, 2018
License     : BSD-3
Maintainer  : capn.freako@gmail.com
Stability   : experimental
Portability : ?
-}
module Haskell_ML.FCN
  ( FCNet(), TrainEvo(..)
  , randNet, runNet, runNet', netTest, hiddenStruct
  , getWeights, getBiases
  , trainNTimes, classificationAccuracy
  ) where

import Control.Monad.Random
import Data.Binary
import Data.List
import Data.Singletons.Prelude
import Data.Singletons.TypeLits
import Data.Vector.Storable (toList, maxIndex)
import GHC.Generics (Generic)
import Numeric.LinearAlgebra.Static hiding (mean)  -- ^ @R@ & @L@ come from here.
-- import Internal.Static (lift1F)  -- Doesn't work; module is hidden.

import Haskell_ML.Util  hiding (classificationAccuracy, maxIndex, getWeights, getBiases)

-- | A fully connected, multi-layer network with fixed input/output
-- widths, but variable (and existentially hidden!) internal structure.
data FCNet :: Nat -> Nat -> * where
  FCNet :: Network i hs o -> FCNet i o

-- | Returns a value of type `FCNet`, filled with random weights
-- ready for training, tucked inside the appropriate Monad, which must
-- be an instance of `MonadRandom` . (IO is such an instance.)
--
-- The input/output widths are determined by the compiler automatically,
-- via type inferencing.
--
-- The internal structure of the network is determined by the list of
-- integers passed in. Each integer in the list indicates the output
-- width of one hidden layer, with the first entry in the list
-- corresponding to the hidden layer nearest to the input layer.
randNet :: (KnownNat i, KnownNat o, MonadRandom m)
        => [Integer]
        -> m (FCNet i o)
randNet hs = withSomeSing hs (fmap FCNet . randNetwork')

-- | Data type for holding training evolution data.
data TrainEvo = TrainEvo
  { accs  :: [Double]                   -- ^ training accuracies
  , diffs :: [([[Double]],[[Double]])]  -- ^ differences of weights/biases, by layer
  }

-- | Train a network on several epochs of the training data, keeping
-- track of accuracy and weight/bias changes per layer, after each.
trainNTimes :: (KnownNat i, KnownNat o)
            => Int           -- ^ Number of epochs
            -> Double        -- ^ learning rate
            -> FCNet i o     -- ^ the network to be trained
            -> [(R i, R o)]  -- ^ the training pairs
            -> (FCNet i o, TrainEvo)
trainNTimes = trainNTimes' [] []

trainNTimes' :: (KnownNat i, KnownNat o)
             => [Double]                    -- accuracies
             -> [([[Double]], [[Double]])]  -- weight/bias differences
             -> Int -> Double -> FCNet i o -> [(R i, R o)] -> (FCNet i o, TrainEvo)
trainNTimes' accs diffs 0 _    net _   = (net, TrainEvo accs diffs)
trainNTimes' accs diffs n rate net prs = trainNTimes' (accs ++ [acc]) (diffs ++ [diff]) (n-1) rate net' prs
  where net'  = trainNet rate net prs
        acc   = classificationAccuracy res ref
        res   = runNet net' $ map fst prs
        ref   = map snd prs
        diff  = ( zipWith (zipWith (-)) (getWeights net') (getWeights net)
                , zipWith (zipWith (-)) (getBiases  net') (getBiases  net) )

-- | Calculate the classification accuracy, given:
--
--   - a list of results vectors, and
--   - a list of reference vectors.
classificationAccuracy :: (KnownNat n) => [R n] -> [R n] -> Double
classificationAccuracy us vs = mean $ cmpr us vs
  where cmpr :: (KnownNat n) => [R n] -> [R n] -> [Double]
        cmpr xs ys = for (zipWith maxComp xs ys) $ \case
                       True  -> 1.0
                       False -> 0.0

        maxComp :: (KnownNat n) => R n -> R n -> Bool
        maxComp u v = maxIndex (extract u) == maxIndex (extract v)

-- | Run a network on a list of inputs.
runNet :: (KnownNat i, KnownNat o)
       => FCNet i o  -- ^ the network to run
       -> [R i]      -- ^ the list of inputs
       -> [R o]      -- ^ the list of outputs
runNet (FCNet n) = map (runNetwork n)

-- | Run a network on a list of inputs,
-- enforcing 4-bit precision for activation outputs.
runNet' :: (KnownNat i, KnownNat o)
       => FCNet i o  -- ^ the network to run
       -> [R i]      -- ^ the list of inputs
       -> [R o]      -- ^ the list of outputs
runNet' (FCNet n) = map (runNetwork' n)

-- | `Binary` instance definition for `FCNet`.
--
-- With this definition, the user of our library is able to use standard
-- `put` and `get` calls, to serialize his created/trained network for
-- future use. And we don't need to provide auxilliary `saveNet` and
-- `loadNet` functions in the API.
instance (KnownNat i, KnownNat o) => Binary (FCNet i o) where
    put = putFCNet
    get = getFCNet

-- | Basic sanity test of our code, taken from Justin's repository.
--
-- Printed output should contain two offset solid circles.
netTest :: MonadRandom m => Double -> Int -> m String
netTest rate n = do
    inps <- replicateM n $ do
      s <- getRandom
      return $ randomVector s Uniform * 2 - 1
    let outs = flip map inps $ \v ->
                 if v `inCircle` (fromRational 0.33, 0.33)
                      || v `inCircle` (fromRational (-0.33), 0.33)
                   then fromRational 1
                   else fromRational 0
    net0 :: Network 2 '[16, 8] 1 <- randNetwork
    let trained = sgd rate (zip inps outs) net0

        outMat = [ [ render (norm_2 (runNetwork trained (vector [x / 25 - 1,y / 10 - 1])))
                   | x <- [0..50] ]
                 | y <- [0..20] ]

        render r | r <= 0.2  = ' '
                 | r <= 0.4  = '.'
                 | r <= 0.6  = '-'
                 | r <= 0.8  = '='
                 | otherwise = '#'

    return $ unlines outMat
  where
    inCircle :: KnownNat n => R n -> (R n, Double) -> Bool
    v `inCircle` (o, r) = norm_2 (v - o) <= r

-- | Returns a list of integers corresponding to the widths of the hidden
-- layers of a `FCNet`.
hiddenStruct :: FCNet i o -> [Integer]
hiddenStruct (FCNet net) = hiddenStruct' net

hiddenStruct' :: Network i hs o -> [Integer]
hiddenStruct' = \case
    W _    -> []
    _ :&~ (n' :: Network h hs' o)
           -> natVal (Proxy @h)
            : hiddenStruct' n'

-- | Returns a list of lists of Doubles, each containing the weights of
-- one layer of the network.
getWeights :: (KnownNat i, KnownNat o) => FCNet i o -> [[Double]]
getWeights (FCNet net) = getWeights' net

getWeights' :: (KnownNat i, KnownNat o) => Network i hs o -> [[Double]]
getWeights' (W Layer{..})       = [concatMap (toList . extract) (toRows nodes)]
getWeights' (Layer{..} :&~ net) = concatMap (toList . extract) (toRows nodes) : getWeights' net

-- | Returns a list of lists of Doubles, each containing the biases of
-- one layer of the network.
getBiases :: (KnownNat i, KnownNat o) => FCNet i o -> [[Double]]
getBiases (FCNet net) = getBiases' net

getBiases' :: (KnownNat i, KnownNat o) => Network i hs o -> [[Double]]
getBiases' (W Layer{..})       = [toList $ extract biases]
getBiases' (Layer{..} :&~ net) = toList (extract biases) : getBiases' net

-----------------------------------------------------------------------
-- All following functions are for internal library use only!
-- They are not exported through the API.
-----------------------------------------------------------------------

-- A single network layer mapping an input of width `i` to an output of
-- width `o`, via simple matrix/vector mult.
data Layer i o = Layer { biases :: !(R o)
                       , nodes  :: !(L o i)
                       }
  deriving (Show, Generic)

instance (KnownNat i, KnownNat o) => Binary (Layer i o)

-- Generates a value of type `Layer i o`, filled with uniformly
-- distributed random values, tucked inside the appropriate Monad, which
-- must be an instance of `MonadRandom`.
--
-- Note: normally distributed values were tried (See commented code, below.)
--       and found not to perform as well as uniformly distributed values.
--
--       This could be due to my use of `logistic`, as opposed to `relu`,
--       as my activation function; I'm not sure.
--       (See code and comments near end of this file.)
randLayer :: forall m i o. (MonadRandom m, KnownNat i, KnownNat o)
          => m (Layer i o)
randLayer = do
  s1 :: Int <- getRandom
  s2 :: Int <- getRandom
  let b = randomVector  s1 Uniform * 2 - 1  -- Need to change to Gaussian equivalent.
      n = uniformSample s2 (-1) 1           -- Need to change to Gaussian equivalent.
  -- PLEASE, KEEP THIS COMMENTED OUT CODE!
  -- let m = eye
  --     b = randomVector s2 Gaussian
  --     n = gaussianSample s1 (takeDiag m) (sym m)
  return $ Layer b n

-- This is the network structure that `FCNet i o` wraps, hiding its
-- internal structure existentially, outside of the library.
data Network :: Nat -> [Nat] -> Nat -> * where
  W     :: !(Layer i o)
        -> Network i '[] o

  (:&~) :: KnownNat h
        => !(Layer i h)
        -> !(Network h hs o)
        -> Network i (h ': hs) o

infixr 5 :&~

-- Generates a value of type `Network i hs o`
-- filled with random weights, ready to begin training.
--
-- Note: `hs` is determined explicitly, via the first argument, while
--       `i` and `o` are determined implicitly, via type inference.
randNetwork :: forall m i hs o. (MonadRandom m, KnownNat i, SingI hs, KnownNat o)
            => m (Network i hs o)
randNetwork = randNetwork' sing

randNetwork' :: forall m i hs o. (MonadRandom m, KnownNat i, KnownNat o)
             => Sing hs -> m (Network i hs o)
randNetwork' = \case
  SNil            -> W     <$> randLayer
  SNat `SCons` ss -> (:&~) <$> randLayer <*> randNetwork' ss

-- Binary instance definition for `Network i hs o`.
putNet :: (KnownNat i, KnownNat o)
       => Network i hs o
       -> Put
putNet = \case
    W w    -> put w
    w :&~ n -> put w *> putNet n

getNet :: forall i hs o. (KnownNat i, KnownNat o)
       => Sing hs
       -> Get (Network i hs o)
getNet = \case
    SNil            -> W    <$> get
    SNat `SCons` ss -> (:&~) <$> get <*> getNet ss

instance (KnownNat i, SingI hs, KnownNat o) => Binary (Network i hs o) where
    put = putNet
    get = getNet sing

putFCNet :: (KnownNat i, KnownNat o)
         => FCNet i o
         -> Put
putFCNet (FCNet net) = do
  put (hiddenStruct' net)
  putNet net

getFCNet :: (KnownNat i, KnownNat o)
         => Get (FCNet i o)
getFCNet = do
  hs <- get
  withSomeSing hs (fmap FCNet . getNet)

runLayer :: (KnownNat i, KnownNat o)
         => Layer i o
         -> R i
         -> R o
runLayer (Layer b n) v = b + n #> v

runNetwork :: (KnownNat i, KnownNat o)
           => Network i hs o
           -> R i
           -> R o
runNetwork = \case
  W w        -> \(!v) -> logistic (runLayer w v)
  (w :&~ n') -> \(!v) -> let v' = logistic (runLayer w v)
                         in runNetwork n' v'

-- 4-bit quantization experiment
runNetwork' :: (KnownNat i, KnownNat o)
            => Network i hs o
            -> R i
            -> R o
runNetwork' = \case
  -- W w -> \(!v) ->
  W Layer{..} -> \(!v) ->
    let vq = quantizeV v
        wq = Layer (quantizeV biases) (quantizeM nodes)
     in quantizeV $ logistic (runLayer wq vq)
  -- (w :&~ n') -> \(!v) ->
  (Layer{..} :&~ n') -> \(!v) ->
    let v' = quantizeV $ logistic (runLayer wq vq)
        vq = quantizeV v
        wq = Layer (quantizeV biases) (quantizeM nodes)
     in runNetwork' n' v'

quantizeV :: KnownNat n => R n -> R n
quantizeV v = dvmap (fromIntegral . floor) (abs v' * 15 + 0.5) / 15
  where v' = dvmap (max 0 . min 1) v

quantizeM :: (KnownNat i, KnownNat o) => L o i -> L o i
quantizeM m = dmmap (fromIntegral . floor) (abs m' * 15 + 0.5) / 15
  where m' = dmmap (max 0 . min 1) m

-- Trains a value of type `FCNet i o`, using the supplied list of
-- training pairs (i.e. - matched input/output vectors).
trainNet :: (KnownNat i, KnownNat o)
         => Double        -- learning rate
         -> FCNet i o     -- the network to be trained
         -> [(R i, R o)]  -- the training pairs
         -> FCNet i o     -- the trained network
trainNet rate (FCNet net) trn_prs = FCNet $ sgd rate trn_prs net

-- Train a network of type `Network i hs o` using a list of training
-- pairs and the Stochastic Gradient Descent (SGD) approach.
sgd :: forall i hs o. (KnownNat i, KnownNat o)
    => Double           -- learning rate
    -> [(R i, R o)]     -- training pairs
    -> Network i hs o   -- network to train
    -> Network i hs o   -- trained network
sgd rate trn_prs net = foldl' (sgdStep rate) net trn_prs

-- Train a network of type `Network i hs o` using a single training pair.
--
-- This code was taken directly from Justin Le's public GitHub archive:
-- https://github.com/mstksg/inCode/blob/43adae31b5689a95be83a72866600033fcf52b50/code-samples/dependent-haskell/NetworkTyped.hs#L77
-- and modified only slightly.
sgdStep :: forall i hs o. (KnownNat i, KnownNat o)
         => Double           -- learning rate
         -> Network i hs o   -- network to train
         -> (R i, R o)       -- training pair
         -> Network i hs o   -- trained network
sgdStep rate net trn_pr = fst $ go x0 net
  where
    x0     = fst trn_pr
    target = snd trn_pr
    go  :: forall j js. KnownNat j
        => R j              -- input vector
        -> Network j js o   -- network to train
        -> (Network j js o, R j)
    go !x (W w@(Layer wB wN))
        = let y    = runLayer w x
              o    = logistic y
              -- the gradient (how much y affects the error)
              --   (logistic' is the derivative of logistic)
              dEdy = logistic' y * (o - target)
              -- new bias weights and node weights
              wB'  = wB - konst rate * dEdy
              wN'  = wN - konst rate * (dEdy `outer` x)
              w'   = Layer wB' wN'
              -- bundle of derivatives for next step
              dWs  = tr wN #> dEdy
          in  (W w', dWs)
    -- Handle the inner layers.
    go !x (w@(Layer wB wN) :&~ n)
        = let y          = runLayer w x
              o          = logistic y
              -- get dWs', bundle of derivatives from rest of the net
              (n', dWs') = go o n
              -- the gradient (how much y affects the error)
              dEdy       = logistic' y * dWs'
              -- new bias weights and node weights
              wB'  = wB - konst rate * dEdy
              wN'  = wN - konst rate * (dEdy `outer` x)
              w'   = Layer wB' wN'
              -- bundle of derivatives for next step
              dWs  = tr wN #> dEdy
          in  (w' :&~ n', dWs)

-- Doesn't work, because the "constructors of R are not in scope."
-- What am I to do, here?!
-- Orphan `Ord` instance, for R n.
-- deriving instance (KnownNat n) => Ord (R n)

-- | Normalize a vector to a probability vector, via softmax.
-- softMax :: (KnownNat n)
--         => R n  -- ^ vector to be normalized
--         -> R n
-- softMax v = exp v / norm_0 v

-- Rectified Linear Unit
-- relu :: (KnownNat n)
--      => R n
--      -> R n
-- relu = max 0

-- relu' :: (KnownNat n)
--       => R n
--       -> R n
-- relu' v = if v > 0 then 1
--                    else 0

-- Logistic non-linear activation function.
logistic :: Floating a => a -> a
logistic x = 1 / (1 + exp (-x))

logistic' :: Floating a => a -> a
logistic' x = logix * (1 - logix)
  where
    logix = logistic x


text:
-- | A tool for measuring cycle time of tasks
--
-- Cycle time is the difference between creation date and completion date.
--
module Main (main) where

import RIO

import Asana.Api
import Asana.Api.Gid (Gid, gidToText)
import Asana.Api.Project (Project(pCreatedAt, pGid, pName), getProjects)
import Asana.App
  (AppM, appExt, loadAppWith, parseBugProjectId, parseYear, runApp)
import qualified Data.Csv as Csv
import Data.Foldable (maximum, minimum)
import Data.List (intercalate, nub)
import qualified Data.Map.Strict as Map
import qualified Data.Vector as V
import RIO.ByteString (writeFile)
import RIO.ByteString.Lazy (toStrict)
import RIO.Text (isPrefixOf, unpack)
import RIO.Time (NominalDiffTime, diffUTCTime, toGregorian, utctDay)
import Statistics.Quantile (median, medianUnbiased)
import Statistics.Sample (kurtosis, mean, skewness, stdDev)
import Statistics.Sample.Histogram (histogram)
import System.IO.Temp (emptySystemTempFile)
import Text.Printf (printf)

data AppExt = AppExt
  { appYear :: Integer
  , appBugProject :: Gid
  }

main :: IO ()
main = do
  app <- loadAppWith $ AppExt <$> parseYear <*> parseBugProjectId
  runApp app $ do
    logDebug "Fetch projects"
    year <- asks $ appYear . appExt
    projects <- filter (approvedProject year) <$> getProjects

    taskGids <- fetchRelevantTaskGids projects

    tasks <- pooledForConcurrentlyN maxRequests taskGids getTask
    logDebug "Write task CSV"
    writeTaskCsv tasks

    let
      cycleTimes = zscoreFilter . V.fromList . fmap realToFrac $ mapMaybe
        mayCycleTime
        tasks

    when (null cycleTimes) $ logWarn "No cycle time to compute"

    logDebug "Write histogram CSV"
    writeHistogramCsv cycleTimes

    logDebug "Display Cycle Time"
    let
      fmt = intercalate
        "\n"
        [ "Cycle Time"
        , "- mean:     %d days"
        , "- mediam:   %d days"
        , "- skewness: %.3f"
        , "- kurtosis: %.3f"
        , "- stdDev:   %d days"
        , "- max:      %d days"
        , "- min:      %d days"
        ]
    logInfo . fromString $ printf
      fmt
      (toDays $ mean cycleTimes)
      (toDays $ median medianUnbiased cycleTimes)
      (skewness cycleTimes)
      (kurtosis cycleTimes)
      (toDays $ stdDev cycleTimes)
      (toDays $ maximum cycleTimes)
      (toDays $ minimum cycleTimes)

zscoreFilter :: Vector Double -> Vector Double
zscoreFilter p =
  let
    av = mean p
    sdev = stdDev p
    zs x = (x - av) / sdev
  in V.filter ((<= 3) . zs) p

toDays :: Double -> Int
toDays = floor . (/ 86400)

approvedProject :: Integer -> Project -> Bool
approvedProject compareYear project = year == compareYear && isPrefixOf
  "Iteration"
  (pName project)
  where (year, _, _) = toGregorian $ utctDay (pCreatedAt project)

mayCycleTime :: Task -> Maybe NominalDiffTime
mayCycleTime task = do
  completedAt <- tCompletedAt task
  -- FIXME: Sometimes this is true
  guard $ completedAt > tCreatedAt task
  Just . diffUTCTime completedAt $ tCreatedAt task

writeHistogramCsv :: Vector Double -> AppM ext ()
writeHistogramCsv cycleTimes = do
  filePath <- liftIO $ emptySystemTempFile ".csv"
  let histogramCsv = uncurry V.zip $ histogram @_ @_ @Double 100 cycleTimes
  writeFile filePath . toStrict . Csv.encode $ V.toList histogramCsv
  logInfo . fromString $ "Histogram written to " <> filePath

writeTaskCsv :: [Task] -> AppM ext ()
writeTaskCsv tasks = do
  filePath <- liftIO $ emptySystemTempFile ".csv"
  let
    toRecord task = Map.fromList
      [ ("name" :: String, unpack $ tName task)
      , ("task ID" :: String, unpack $ taskUrl task)
      , ( "cycle time"
        , maybe "" (show . realToFrac @_ @Double) $ mayCycleTime task
        )
      ]
  writeFile filePath
    . toStrict
    . Csv.encodeByName (V.fromList ["name", "task ID", "cycle time"])
    $ toRecord
    <$> tasks
  logInfo . fromString $ "Tasks written to " <> filePath

fetchRelevantTaskGids :: Traversable t => t Project -> AppM AppExt [Gid]
fetchRelevantTaskGids projects = do
  taskGids <-
    fmap (nub . concat)
    . pooledForConcurrentlyN maxRequests projects
    $ \project -> do
        logDebug . fromString $ "Project tasks: " <> show
          (gidToText $ pGid project)
        fmap nGid <$> getProjectTasks (pGid project) AllTasks
  bugProjectGid <- asks $ appBugProject . appExt
  bugTaskGids <- fmap nGid <$> getProjectTasks bugProjectGid AllTasks

  pure $ filter (`notElem` bugTaskGids) taskGids

max_stars_repo_name: freckle/asana

max_issues_repo_name: freckle/asana

max_forks_repo_name: freckle/asana


text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE DeriveDataTypeable #-}
#if __GLASGOW_HASKELL__ >= 786
{-# LANGUAGE PatternSynonyms #-}
{-# LANGUAGE ViewPatterns #-}
#endif

module Data.Complex.Polar (
#if __GLASGOW_HASKELL__ >= 800
    Polar((:<)),
#elif __GLASGOW_HASKELL >= 786
    Polar,
    pattern (:<),
#else
    Polar,
#endif
    fromPolar,
    fromComplex,
    realPart,
    imagPart,
    conjugate,
    mkPolar,
    cis,
    polar,
    magnitude,
    phase
) where

import           Data.Complex (Complex(..))
import qualified Data.Complex as C

import Data.Typeable
#ifdef __GLASGOW_HASKELL__
import Data.Data (Data)
#endif

#ifdef __HUGS__
import Hugs.Prelude(Num(fromInt), Fractional(fromDouble))
#endif

import Text.Read (parens)
import Text.ParserCombinators.ReadPrec (prec, step)
import Text.Read.Lex (Lexeme (Ident))

infix 6 :<<

-- -----------------------------------------------------------------------------
-- The Polar type

-- | Complex numbers are an algebraic type.
--
-- For a complex number @z@, @'abs' z@ is a number with the magnitude of @z@,
-- but oriented in the positive real direction, whereas @'signum' z@
-- has the phase of @z@, but unit magnitude.
data (RealFloat a) => Polar a = !a :<< !a -- ^ forms a complex number from its magnitude
                                          --   and its phase in radians.
#if __GLASGOW_HASKELL__
    deriving (Eq, Data, Typeable)
#else
    deriving Eq
#endif

#if __GLASGOW_HASKELL__ >= 786
infix 6 :<

-- | Smart constructor that canonicalizes the magnitude and phase.
pattern r :< theta <-
  ( \(r :<< theta) -> Just (r, theta) ->
      Just (r, theta)
    )
  where
    r :< theta = mkPolar r theta
#endif

instance (RealFloat a, Show a) => Show (Polar a) where
    {-# SPECIALISE instance Show (Polar Float) #-}
    {-# SPECIALISE instance Show (Polar Double) #-}
    showsPrec d (r :<< theta) =
      showParen (d >= 11) (showString "mkPolar "
                           . showsPrec 11 r
                           . showString " "
                           . showsPrec 11 theta)

instance (RealFloat a, Read a) => Read (Polar a) where
    {-# SPECIALISE instance Read (Polar Float) #-}
    {-# SPECIALISE instance Read (Polar Double) #-}
    readsPrec d = readParen (d > 10)
                        (\s -> do ("mkPolar", s2) <- lex s
                                  (r, s3) <- readsPrec 11 s2
                                  (theta, s4) <- readsPrec 11 s3
                                  return (mkPolar r theta, s4))

-- | Wrap phase back in interval @(-'pi', 'pi']@.
wrap :: (RealFloat a) => a -> a
{-# SPECIALISE wrap :: Float  -> Float   #-}
{-# SPECIALISE wrap :: Double -> Double  #-}
wrap phi | phi <= (-pi) = wrap (phi+2*pi)
wrap phi | phi > pi     = wrap (phi-2*pi)
wrap phi                = phi

-- | Convert to rectangular form.
fromPolar :: (RealFloat a) => Polar a -> Complex a
{-# INLINE fromPolar #-}
fromPolar p = realPart p :+ imagPart p

-- | Convert to polar form.
fromComplex :: (RealFloat a) => Complex a -> Polar a
{-# INLINE fromComplex #-}
fromComplex = uncurry mkPolar_ . C.polar

-- | Extracts the real part of a complex number.
realPart :: (RealFloat a) => Polar a -> a
realPart (r :<< theta) = r * cos theta

-- | Extracts the imaginary part of a complex number.
imagPart :: (RealFloat a) => Polar a -> a
imagPart (r :<< theta) = r * sin theta

-- | The conjugate of a complex number.
conjugate :: (RealFloat a) => Polar a -> Polar a
{-# SPECIALISE conjugate :: Polar Float  -> Polar Float #-}
{-# SPECIALISE conjugate :: Polar Double -> Polar Double #-}
conjugate (r :<< theta) = mkPolar r (negate theta)

-- | Form a complex number from polar components of magnitude and phase.
-- The magnitude and phase are expected to be in canonical form.
mkPolar_ :: RealFloat a => a -> a -> Polar a
{-# SPECIALISE mkPolar_ :: Float  -> Float  -> Polar Float  #-}
{-# SPECIALISE mkPolar_ :: Double -> Double -> Polar Double #-}
mkPolar_ r theta = r :<< theta

-- | Form a complex number from polar components of magnitude and phase.
mkPolar :: RealFloat a => a -> a -> Polar a
{-# SPECIALISE mkPolar :: Float  -> Float  -> Polar Float  #-}
{-# SPECIALISE mkPolar :: Double -> Double -> Polar Double #-}
mkPolar r theta | r == 0 = 0 :<< 0
mkPolar r theta | r < 0 = mkPolar (- r) (theta + pi)
mkPolar r theta = mkPolar_ r (wrap theta)

-- | @'cis' t@ is a complex value with magnitude @1@
-- and phase @t@ (modulo @2*'pi'@).
cis :: (RealFloat a) => a -> Polar a
{-# SPECIALISE cis :: Float  -> Polar Float  #-}
{-# SPECIALISE cis :: Double -> Polar Double #-}
cis theta = mkPolar 1 theta

-- | The function 'polar' takes a complex number and
-- returns a (magnitude, phase) pair in canonical form:
-- the magnitude is nonnegative, and the phase in the range @(-'pi', 'pi']@;
-- if the magnitude is zero, then so is the phase.
polar :: (RealFloat a) => Polar a -> (a,a)
{-# SPECIALISE polar :: Polar Double -> (Double,Double) #-}
{-# SPECIALISE polar :: Polar Float  -> (Float,Float)   #-}
polar (r :<< theta) = (r,theta)

-- | The nonnegative magnitude of a complex number.
magnitude :: (RealFloat a) => Polar a -> a
{-# SPECIALISE magnitude :: Polar Float  -> Float  #-}
{-# SPECIALISE magnitude :: Polar Double -> Double #-}
magnitude (r :<< _) = r

-- | The phase of a complex number, in the range @(-'pi', 'pi']@.
-- If the magnitude is zero, then so is the phase.
phase :: (RealFloat a) => Polar a -> a
{-# SPECIALISE phase :: Polar Float  -> Float  #-}
{-# SPECIALISE phase :: Polar Double -> Double #-}
phase (_ :<< theta) = theta

instance (RealFloat a) => Num (Polar a) where
    {-# SPECIALISE instance Num (Polar Float)  #-}
    {-# SPECIALISE instance Num (Polar Double) #-}
    z + z'           = fromComplex (fromPolar z + fromPolar z')
    z - z'           = fromComplex (fromPolar z - fromPolar z')
    z * z'           = mkPolar  (magnitude z * magnitude z') (phase z + phase z')
    negate z         = mkPolar  (negate (magnitude z)) (phase z)
    abs z            = mkPolar_ (magnitude z) 0
    signum (0 :<< _) = 0
    signum z         = mkPolar_ 1 (phase z)
    fromInteger      = flip mkPolar 0 . fromInteger

instance (RealFloat a) => Fractional (Polar a) where
    {-# SPECIALISE instance Fractional (Polar Float)  #-}
    {-# SPECIALISE instance Fractional (Polar Double) #-}
    z / z'          = mkPolar (magnitude z / magnitude z') (phase z - phase z')
    fromRational r  = mkPolar (fromRational r) 0

instance  (RealFloat a) => Floating (Polar a)  where
    {-# SPECIALISE instance Floating (Polar Float) #-}
    {-# SPECIALISE instance Floating (Polar Double) #-}
    pi             = mkPolar_ pi 0
    exp (r :<< theta) = mkPolar (exp (r * cos theta)) (r * sin theta)
    log (r :<< theta) = fromComplex (log r :+ theta)

    -- sqrt (0 :<< _)       =  0
    -- sqrt z@(r :<< theta) =  u :+ (if y < 0 then -v else v)
    --                           where (u,v) = if x < 0 then (v',u') else (u',v')
    --                                 v'    = abs y / (u'*2)
    --                                 u'    = sqrt ((magnitude z + abs x) / 2)
    sqrt = fromComplex.sqrt.fromPolar

    -- sin (x:+y)     =  sin x * cosh y :+ cos x * sinh y
    sin = fromComplex.sin.fromPolar

    -- cos (x:+y)     =  cos x * cosh y :+ (- sin x * sinh y)
    cos = fromComplex.cos.fromPolar

    -- tan (x:+y)     =  (sinx*coshy:+cosx*sinhy)/(cosx*coshy:+(-sinx*sinhy))
    --                   where sinx  = sin x
    --                         cosx  = cos x
    --                         sinhy = sinh y
    --                         coshy = cosh y
    tan = fromComplex.tan.fromPolar

    -- sinh (x:+y)    =  cos y * sinh x :+ sin  y * cosh x
    sinh = fromComplex.sinh.fromPolar

    -- cosh (x:+y)    =  cos y * cosh x :+ sin y * sinh x
    cosh = fromComplex.cosh.fromPolar

    -- tanh (x:+y)    =  (cosy*sinhx:+siny*coshx)/(cosy*coshx:+siny*sinhx)
    --                   where siny  = sin y
    --                         cosy  = cos y
    --                         sinhx = sinh x
    --                         coshx = cosh x
    tanh = fromComplex.tanh.fromPolar

    -- asin z@(x:+y)  =  y':+(-x')
    --                   where  (x':+y') = log (((-y):+x) + sqrt (1 - z*z))
    asin = fromComplex.asin.fromPolar

    -- acos z         =  y'':+(-x'')
    --                   where (x'':+y'') = log (z + ((-y'):+x'))
    --                         (x':+y')   = sqrt (1 - z*z)
    acos = fromComplex.acos.fromPolar

    -- atan z@(x:+y)  =  y':+(-x')
    --                   where (x':+y') = log (((1-y):+x) / sqrt (1+z*z))
    atan = fromComplex.atan.fromPolar

    asinh z        =  log (z + sqrt (1+z*z))
    acosh z        =  log (z + (z+1) * sqrt ((z-1)/(z+1)))
    atanh z        =  log ((1+z) / sqrt (1-z*z))

max_stars_repo_name: kaoskorobase/polar

max_issues_repo_name: kaoskorobase/polar

max_forks_repo_name: kaoskorobase/polar


text:
{-# LANGUAGE FlexibleContexts #-}
module InterpolationAlgorithms
  (riemannNormCoord,estimateGrad)
  where

import           Control.Monad
import           Control.Monad.ST
import           Numeric.LinearAlgebra
import           Numeric.LinearAlgebra.Devel

consInterpMat :: Int
              -> Int
              -> [Vector Double]
              -> Matrix Double
consInterpMat rowlen collen vecs = runST $ do
  m <- unsafeThawMatrix mat
  mapM_ (\(i, v) -> write' m i $ flatOuter v) vecs'
  unsafeFreezeMatrix m
  where
    flatOuter vec' = vjoin [1, vec', flatten . join outer $ vec'] :: Vector Double
    collen' = 1 + collen + (collen * collen)
    mat = konst 0 (rowlen, collen') :: Matrix Double
    vecs' = zip [0 :: Int ..] vecs
    write' m i v' = setMatrix m i 0 (asRow v')

estimateGrad :: Int
             -> Int
             -> Matrix Double
             -> Matrix Double
estimateGrad nums gradlen pcaMat = mat
  where
    mat = consInterpMat nums gradlen pcaMat'
    pcaMat' = map (subVector 1 gradlen . (<# getPCAEigVec pcaMat)) matRows
    matRows = toRows pcaMat

riemannNormCoord :: Matrix Double -> Matrix Double -> [Vector Double]
riemannNormCoord res mat = toColumns leastsquaresSol
  where
    leastsquaresSol = cmap (* (-0.5)) $ (res <\> mat) ?? (Pos $ idxs [1, 2], All)

getPCAEigVec :: Matrix Double -> Matrix Double
getPCAEigVec = snd . eigSH . trustSym . snd . meanCov

max_stars_repo_name: emmanueldenloye/manifoldRNC

max_issues_repo_name: emmanueldenloye/manifoldRNC

max_forks_repo_name: emmanueldenloye/manifoldRNC


text:
{-# LANGUAGE DeriveFunctor, GeneralizedNewtypeDeriving #-}
{-# LANGUAGE DataKinds, TypeFamilies #-}
{-# LANGUAGE FlexibleContexts, FlexibleInstances #-}

module DelayedSamplingPure where

import Prelude hiding ((<>))

import Control.Exception.Base (assert)
import Control.Monad (forM_, when)

import Data.List (delete)
import Data.Maybe (isJust)
import Data.Proxy
import Data.IORef

import Debug.Trace

import GHC.TypeLits

import Numeric.AD
import Numeric.AD.Rank1.Tower (Tower)

import Numeric.LinearAlgebra.Static hiding (M)

import Unsafe.Coerce (unsafeCoerce)

import Distributions
import MVDistributions
import Inference

data Node' b where
  NMarginalized :: Maybe (CDistr b c, Node' c) -> [InitTree b] -> Node' b

data RootNode a where
  RootNode :: MDistr a -> Node' a -> RootNode a

data InitTree b where
  InitTree :: CDistr b c -> [InitTree c] -> InitTree b

data MPath where
  Root :: MPath
  Child :: MPath -> MPath

data IPath = IPath [Int]

data AnyPath where
  AMPath :: MPath -> AnyPath
  AIPath :: IPath -> AnyPath

data MarginalT = MGaussianT | MBetaT | MBernoulliT | forall (n :: Nat). MMVGaussianT (Proxy n)

data SMarginalT (m :: MarginalT) where
  SMGaussianT :: SMarginalT MGaussianT
  SMBetaT :: SMarginalT MBetaT
  SMBernoulliT :: SMarginalT MBernoulliT
  SMMVGaussianT :: forall (n :: Nat) (p :: Proxy n). SMarginalT (MMVGaussianT p)

type family MType (m :: MarginalT) where
  MType MGaussianT = Double
  MType MBetaT = Double
  MType MBernoulliT = Bool
  MType (MMVGaussianT (p :: Proxy n)) = R n

data MDistr (m :: MarginalT) where
  MGaussian :: !Double -> !Double -> MDistr MGaussianT
  MBeta :: !Double -> !Double -> MDistr MBetaT
  MBernoulli :: !Double -> MDistr MBernoulliT
  MMVGaussian :: forall (n :: Nat) (p :: Proxy n). KnownNat n => !(R n) -> !(Sym n) -> MDistr (MMVGaussianT p)

--deriving instance Eq (MDistr m)
deriving instance Show (MDistr m)

data CDistr (m :: MarginalT) (m' :: MarginalT) where
  AffineMeanGaussian :: !Double -> !Double -> !Double -> CDistr MGaussianT MGaussianT
  CBernoulli :: CDistr MBetaT MBernoulliT
  MVAffineMeanGaussian :: forall (n :: Nat) (m :: Nat) (pn :: Proxy n) (pm :: Proxy m).
    (KnownNat n, KnownNat m) => (L m n) -> (R m) -> (Sym m) -> CDistr (MMVGaussianT pn) (MMVGaussianT pm)

gaussianConditioning :: Double -> Double -> Double -> Double -> (Double, Double)
gaussianConditioning mu var obs obsvar = (mu', var')
  where
  ivar = recip var
  iobsvar = recip obsvar
  inf = ivar + iobsvar
  var' = recip inf
  mu' = (ivar * mu + iobsvar * obs) / inf

gaussianMeanGaussian :: Double -> CDistr MGaussianT MGaussianT
gaussianMeanGaussian = AffineMeanGaussian 1 0

mdistrToDistr :: MDistr a -> Distr (MType a)
mdistrToDistr (MGaussian mu var) = normal mu var
mdistrToDistr (MMVGaussian mu var) = mvNormal mu var
mdistrToDistr (MBeta a b) = beta a b
mdistrToDistr (MBernoulli p) = bernoulli p

cdistrToDistr :: CDistr m m' -> MType m -> MDistr m'
cdistrToDistr (AffineMeanGaussian m b obsvar) mu = MGaussian (m * mu + b) obsvar
cdistrToDistr (MVAffineMeanGaussian f b obsvar) mu = MMVGaussian (f #> mu + b) obsvar
cdistrToDistr CBernoulli p = MBernoulli p

makeMarginal :: MDistr a -> CDistr a b -> MDistr b
makeMarginal (MGaussian mu var) (AffineMeanGaussian m b obsvar) =
  MGaussian (m * mu + b) (m^2 * var + obsvar)
makeMarginal (MMVGaussian mu var) (MVAffineMeanGaussian f b obsvar) =
  MMVGaussian (f #> mu + b) (conjugate f var + obsvar)
makeMarginal (MBeta a b) CBernoulli = MBernoulli (a / (a + b))
makeMarginal (MBernoulli _) _ = error "impossible"

makeConditional :: MDistr a -> CDistr a b -> MType b -> MDistr a
makeConditional (MGaussian mu var) (AffineMeanGaussian m b obsvar) obs =
  MGaussian mu' var'
  where (mu', var') = gaussianConditioning mu var ((obs - b) / m) (obsvar / m^2)
makeConditional (MMVGaussian mu var) (MVAffineMeanGaussian f b obsvar) obs =
  MMVGaussian mu' var'
  where (mu', var') = fst (kalmanUpdate obs f obsvar (mu, var))
makeConditional (MBeta a b) CBernoulli tf =
  if tf then MBeta (a + 1) b else MBeta a (b + 1)
makeConditional (MBernoulli _) _ _ = error "impossible"

data State b where
  Initialized :: State b
  Marginalized :: MDistr b -> State b

deriving instance Show (State MGaussianT)
deriving instance Show (State MBetaT)

isMarginalized :: State b -> Bool
isMarginalized (Marginalized _) = True
isMarginalized _ = False

data DSDistr a b where
   UDistr :: !(MDistr b) -> DSDistr a b
   CDistr :: !(IORef (Node z a)) -> !(CDistr a b) -> DSDistr a b

typeOfMDistr :: MDistr b -> SMarginalT b
typeOfMDistr (MGaussian _ _) = SMGaussianT
typeOfMDistr (MBeta _ _) = SMBetaT
typeOfMDistr (MBernoulli _) = SMBernoulliT
typeOfMDistr (MMVGaussian _ _) = SMMVGaussianT

typeOfCDistr :: CDistr m m' -> SMarginalT m'
typeOfCDistr (AffineMeanGaussian _ _ _) = SMGaussianT
typeOfCDistr CBernoulli = SMBernoulliT
typeOfCDistr (MVAffineMeanGaussian _ _ _) = SMMVGaussianT

typeOfDSDistr :: DSDistr a b -> SMarginalT b
typeOfDSDistr (UDistr d) = typeOfMDistr d
typeOfDSDistr (CDistr _ c) = typeOfCDistr c

data Node a b = Node
  { name :: !String
  , children :: [RefNodeFrom b]
  , state :: !(State b)
  , distr :: !(DSDistr a b)
  }
  | RealizedNode (MType b)

data Path e (a :: MarginalT) (b :: MarginalT) where
  PNil :: Path e a a
  PSnoc :: e b c -> Path e a b -> Path e a c

data PathTo e b where
  PathTo :: !(Path e a b) -> PathTo e b

data RefNodeTo b where
  RefNodeTo :: !(IORef (Node a b)) -> RefNodeTo b

data RefNodeFrom a where
  RefNodeFrom :: !(IORef (Node a b)) -> RefNodeFrom a

instance Eq (RefNodeFrom a) where
  RefNodeFrom x == RefNodeFrom y = x == unsafeCoerce y

instance Eq (RefNodeTo a) where
  RefNodeTo x == RefNodeTo y = x == unsafeCoerce y

-- initialize without parent node
constant' :: MDistr a -> RootNode a
constant' d = RootNode d (NMarginalized Nothing [])

-- mModifyChildren :: ([InitTree a] -> [InitTree a]) -> MPath -> Node' x -> Node' x
-- mModifyChildren f Root (NMarginalized c t) = NMarginalized c (f t)
-- mModifyChildren f (Child p) (NMarginalized (Just (cdistr, c)) t) = NMarginalized (Just (cdistr, mModifyChildren f p c)) t

-- modifyChildren :: ([InitTree a] -> [InitTree a]) -> AnyPath -> Node' x -> Node' x
-- modifyChildren f (AMPath p) = mModifyChildren f p

-- -- initialize with parent node
-- newConditional' :: AnyPath -> CDistr a b -> Node' x -> (Node' x, AnyPath)
-- newConditional' par cdistr = do
--   childRef <- newIORef child
--   modifyIORef' par $ \n -> n { children = (RefNodeFrom childRef) : children n }
--   --childRef' <- mkWeakIORef childRef $ putStrLn ("deleted " ++ str)
--   pure childRef
--   where
--   child = Node
--     { name = str
--     , children = []
--     , state = Initialized
--     , distr = CDistr par cdistr }

max_stars_repo_name: psg-mit/probzelus-haskell

max_issues_repo_name: psg-mit/probzelus-haskell

max_forks_repo_name: psg-mit/probzelus-haskell


text:
-- Copyright (c) Facebook, Inc. and its affiliates.

{-# LANGUAGE CPP #-}
module UnitTests (main) where

import Test.HUnit
import TestRunner

import Util.ByteString
import Util.Aeson
import qualified Util.Text as UT
import Data.Aeson.Types
import Data.Aeson (decode)
import Data.ByteString (packCString, packCStringLen, ByteString)
import Data.Complex
import qualified Data.Vector as V
import Numeric.FFT.Vector.Unitary as FFTW

import Foreign.Storable
import Foreign.C.String
import Foreign.C.Types
import Foreign.Marshal.Array (advancePtr)
import Foreign.Ptr (nullPtr, Ptr)

import TextShow

intToByteStringTest :: Test
intToByteStringTest = TestCase $ do
  assertEqual "min int" (intToByteString minBound) "-9223372036854775808"
  assertEqual "max int" (intToByteString maxBound) "9223372036854775807"
  assertEqual "zero" (intToByteString 0) "0"

parseValueStrictTest :: Test
parseValueStrictTest = TestCase $ do
  assertEqual "false" (Right $ Bool False) $ parseValueStrict' "false"
  assertEqual "0.1" (Right $ Number 0.1) $ parseValueStrict' "0.1"
  assertEqual "consume all input" (Left "endOfInput") $
    parseValueStrict' "falsee"

#ifdef FACEBOOK
-- | We've patched the official aeson library to handle NaN. See #5183005
-- This test ensures that we notice if we forget to patch on upgrade.
parseNaN :: Test
parseNaN = TestCase $
  case decode "[NaN]" of
    Just [a :: Double] ->
      assertBool "NaN" $ isNaN a
    _ -> assertFailure "couldn't parse NaN"
#endif

-- | Check that we can still call vector-fftw after T22849884
useFFTW :: Test
useFFTW = TestCase $
  assertEqual "useFFTW" [17,5,2,1,1,1,1,1,2,5] $
    applyFFT [1..10]
  where
    applyFFT :: [Int] -> [Int]
    applyFFT series =
      map (floor . magnitude)
      . V.toList
      . FFTW.execute (FFTW.plan FFTW.dft $ length series)
      . V.fromList
      $ map (doubleToComplex . fromIntegral) series
      where
        doubleToComplex :: Double -> Complex Double
        doubleToComplex d = d :+ 0.0

withCStringLenTest :: Test
withCStringLenTest = TestCase $
  UT.withCStringLen "Hello World" $ \(ptr, _) -> do
    val <- peek ptr
    assertEqual "first character" val (castCharToCChar 'H')

useByteStringsAsCStringsTest :: Test
useByteStringsAsCStringsTest = TestCase $
  useByteStringsAsCStrings ["hello", "world"] $ \ptr -> do
    hello <- packCString =<< peek ptr
    world <- packCString =<< peek (advancePtr ptr 1)
    term <- peek $ advancePtr ptr 2
    assertEqual "first word" hello "hello"
    assertEqual "second word" world "world"
    assertEqual "terminator" term nullPtr

-- This should probably be a QuickCheck Property
bsListAsCStrLenArrTest :: [ByteString] -> Test
bsListAsCStrLenArrTest bs = TestCase $
  bsListAsCStrLenArr bs $ \strPtr lenPtr len -> do
    let packBS :: CSize -> IO ByteString
        packBS i = packCStringLen =<< (,)
                   <$> peek (getPtr strPtr)
                   <*> (fromIntegral <$> peek (getPtr lenPtr))
            where getPtr :: Storable a => Ptr a -> Ptr a
                  getPtr = (`advancePtr` fromIntegral i)
    got <- mapM packBS [0..len - 1]
    assertEqual "same list" got bs

data MyEnum = MyA | MyB
  deriving (Bounded, Enum, Eq, Show)

instance UT.TextRead MyEnum
instance TextShow MyEnum where
  showb MyA = "MyA"
  showb MyB = "MyB"

readTextTest :: Test
readTextTest = TestCase $ do
  assertEqual "myA" (Just MyA) (UT.readText "MyA")
  assertEqual "Nothing" (Nothing :: Maybe MyEnum) (UT.readText "MyC")

insertCommasAndAndTest :: Test
insertCommasAndAndTest = TestCase $ do
  assertEqual "null" "" $ UT.insertCommasAndAnd []
  assertEqual "singleton" "foo" $ UT.insertCommasAndAnd ["foo"]
  assertEqual "pair" "foo and bar" $ UT.insertCommasAndAnd ["foo", "bar"]
  assertEqual "triple" "foo, bar, and baz" $
    UT.insertCommasAndAnd ["foo", "bar", "baz"]
  assertEqual "quad" "foo, bar, baz, and quux" $
    UT.insertCommasAndAnd ["foo", "bar", "baz", "quux"]
  assertEqual "10" "1, 2, 3, 4, 5, 6, 7, 8, 9, and 10" $
    UT.insertCommasAndAnd [ showt (i::Int) | i <- [1..10] ]

main :: IO ()
main = testRunner $ TestList
  [ TestLabel "intToByteString" intToByteStringTest
  , TestLabel "parseValueStrictTest" parseValueStrictTest
#ifdef FACEBOOK
  , TestLabel "parseNaN" parseNaN
#endif
  , TestLabel "useFFTW" useFFTW
  , TestLabel "withCStringLen" withCStringLenTest
  , TestLabel "useByteStringsAsCStrings" useByteStringsAsCStringsTest
  , TestLabel "bsListAsCStrLenArr" (bsListAsCStrLenArrTest ["hello", "world"])
  , TestLabel "readText" readTextTest
  , TestLabel "insertCommasAndAnd" insertCommasAndAndTest
  ]

max_stars_repo_name: nhawkes/hsthrift

max_issues_repo_name: nhawkes/hsthrift

max_forks_repo_name: nhawkes/hsthrift


text:
--------------------------------------------------------------------------------
--
--  Copyright (c) 2010 - 2013 Tad Doxsee
--  All rights reserved.
--
--  Author: Tad Doxsee
--
--------------------------------------------------------------------------------

{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE  FlexibleInstances    #-}

-- * base
import Foreign           (ForeignPtr, Storable)
import Foreign.C.Types   (CDouble, CInt, CSize)

-- * vector
import qualified Data.Vector.Storable.Mutable as V

-- * cholmod
import Numeric.LinearAlgebra.CHOLMOD.CholmodXFaceLow
import Numeric.LinearAlgebra.CHOLMOD.CholmodXFace

--------------------------------------------------------------------------------

main :: IO ()
main = do

  c <- allocCommon :: IO (ForeignPtr Common)
  startC c

  let n = 5                  :: CSize
      nz = (n+1) * n `div` 2 :: CSize
      nn = fromIntegral n    :: CSize

  at <-allocTriplet n n nz stSquareSymmetricLower xtReal c
       :: IO (Matrix Triplet)

  atNRow  <- getNRow  at :: IO CSize
  atNCol  <- getNCol  at :: IO CSize
  atNZMax <- getNZMax at :: IO CSize

  putStrLn ""
  putStrLn $ "at: (" ++ show atNRow
                     ++ " x " ++ show atNCol ++ ")"
  putStrLn $ "nzmax: " ++ show atNZMax

  iv <- tripletGetRowIndices at :: IO (V.IOVector CInt)
  jv <- tripletGetColIndices at :: IO (V.IOVector CInt)
  xt  <- tripletGetX at

  let nni = fromIntegral nn

  let ij      = [(i, j) | i <- [0 .. nni-1], j <- [0 .. i]] :: [(CInt, CInt)]
      (ia,ja) = unzip ij                                    :: ([CInt], [CInt])

      -- the elements of the matrix are (lower half)
      --
      --  11.0
      --  21.0  22.0
      --  31.0  32.0  33.0
      --  41.0  42.0  43.0  44.0
      --  51.0  52.0  53.0  54.0  55.0

      xp = [11, 21, 22, 31, 32, 33, 41, 42, 43, 44, 51, 52, 53, 54, 55]
           :: [CDouble]

  writev iv ia
  writev jv ja
  writev xt xp

  let nnz = fromIntegral $ Prelude.length xp :: CSize

  tripletSetNNZ at nnz

  as <- tripletToSparse at c :: IO (Matrix Sparse)
  b  <- ones n 1 xtReal c    :: IO (Matrix Dense)
  l  <- analyze as c         :: IO (ForeignPtr Factor)
  factorize as l c
  x  <- solve stA l b c      :: IO (Matrix Dense)

  xNRow <- getNRow x
  xNCol <- getNCol x

  xv <- getX x :: IO (V.IOVector CDouble)

  putStrLn ""
  putStrLn $ "xv: (" ++ show xNRow
                     ++ " x " ++ show xNCol ++ ")"

  xl <- readv xv
  mapM_ (putStrLn . show) xl

  r <- denseCopy b c

  let oneL = [1, 0] :: [CDouble]
      m1L = [-1, 0] :: [CDouble]

  _ <- sdMult as noTranspose m1L oneL x r c

  norm <- denseNorm r infNorm c :: IO CDouble

  putStrLn $ "norm: " ++ show norm

  putStrLn "done"

writev :: (Storable a) => V.IOVector a -> [a] -> IO ()

writev v xs =
  sequence_ [V.write v i x | (i, x) <- zip [0 .. (Prelude.length xs - 1)] xs]

readv :: (Storable a) => V.IOVector a -> IO [a]
readv v = sequence [V.read v i | i <- [0 .. (V.length v) - 1]]

max_stars_repo_name: tdox/hcholmod

max_issues_repo_name: tdox/hcholmod

max_forks_repo_name: tdox/hcholmod


text:
-- Copyright 2016 TensorFlow authors.
--
-- Licensed under the Apache License, Version 2.0 (the "License");
-- you may not use this file except in compliance with the License.
-- You may obtain a copy of the License at
--
--     http://www.apache.org/licenses/LICENSE-2.0
--
-- Unless required by applicable law or agreed to in writing, software
-- distributed under the License is distributed on an "AS IS" BASIS,
-- WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-- See the License for the specific language governing permissions and
-- limitations under the License.

-- | This module contains definitions for some built-in TensorFlow operations.
--
-- Note that certain, "stateful" ops like 'variable' and 'assign' return a
-- 'Build' action (e.g., @Build (Tensor Ref a)@ instead of a pure value; the
-- returned 'Tensor's are always rendered in the current 'Build' context.  This
-- approach helps us avoid problems with inlining or common subexpression
-- elimination, by writing
--
-- > do
-- >     v <- variable []
-- >     w <- assign v 3
-- >     render $ w * w
--
-- instead of
--
-- > let
-- >    v = variable []
-- >    w = assign v 3
-- > in w * w
--
-- since the latter could be reasonably transformed by the compiler into (or
-- vice versa)
--
-- > let
-- >    v = variable []
-- >    w = assign v 3
-- >    w' = assign v 3
-- > in w * w'
--
-- Ops should return a 'Build' action if their original 'OpDef' marks them as
-- stateful, or if they take any Refs as input.  (This mirrors the rules that
-- TensorFlow uses to avoid common subexpression elimination.)
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE OverloadedLists #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE UndecidableInstances #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}

module TensorFlow.Ops
    ( CoreOps.add
    , CoreOps.abs
    , CoreOps.addN
    , CoreOps.argMax
    , CoreOps.assign
    , CoreOps.broadcastGradientArgs
    , CoreOps.cast
    , CoreOps.concat
    , constant
    , CoreOps.equal
    , expandDims
    , initializedVariable
    , zeroInitializedVariable
    , CoreOps.fill
    , CoreOps.oneHot
    , CoreOps.matMul
    , matTranspose
    , CoreOps.mean
    , CoreOps.mul
    , CoreOps.neg
    , CoreOps.pack
    , placeholder
    , CoreOps.range
    , reducedShape
    , CoreOps.relu
    , CoreOps.reluGrad
    , CoreOps.reshape
    , restore
    , restoreFromName
    , save
    , scalar
    , shape
    , CoreOps.sign
    , CoreOps.size
    , CoreOps.softmax
    , CoreOps.softmaxCrossEntropyWithLogits
    , CoreOps.sparseToDense
    , CoreOps.sub
    , CoreOps.sum
    , CoreOps.transpose
    , truncatedNormal
    , CoreOps.variable
    , vector
    , zeros
    , CoreOps.zerosLike
    , scalarize
    ) where

import Data.ByteString (ByteString)
import Data.Complex (Complex)
import Data.Int (Int32, Int64)
import Prelude hiding (abs, sum, concat)
import Data.ProtoLens (def)
import Data.Text.Encoding (encodeUtf8)
import Lens.Family2 ((.~), (&))
import Text.Printf (printf)
import Proto.Tensorflow.Core.Framework.Tensor
    ( TensorProto
    , dtype
    , tensorShape
    )
import qualified Proto.Tensorflow.Core.Framework.TensorShape
  as TensorShape
import TensorFlow.Build
import TensorFlow.BuildOp
import TensorFlow.ControlFlow (group)
import TensorFlow.Output (unNodeName)
import TensorFlow.Tensor
import TensorFlow.Types

import qualified TensorFlow.GenOps.Core as CoreOps

import qualified Prelude (abs)

-- TODO: Look into hs-boot refactoring to allow mutually recursive imports.
-- | Must be defined as an orphan because of the dependency order between Ops
-- and Tensor.
--
-- The indirect constraint "v ~ Value" helps disambiguate types, for example in
-- "neg 1 :: Tensor Value Float", it helps find the type of the subexpression
-- "1".
instance ( TensorType a
         , Num a
         , v ~ Value
         , OneOf '[ Double, Float, Int32, Int64
                  , Complex Float, Complex Double] a) => Num (Tensor v a) where
    (+) = CoreOps.add
    (*) = CoreOps.mul
    (-) = CoreOps.sub
    abs = CoreOps.abs
    fromInteger = scalar . fromInteger
    signum = CoreOps.sign
    negate = CoreOps.neg

matTranspose :: forall a v . TensorType a
             => Tensor v a -> Tensor Value a
matTranspose = flip CoreOps.transpose (vector [1, 0 :: Int32])

placeholder :: forall a . TensorType a => Shape -> Build (Tensor Value a)
placeholder shape' =
    buildOp $ opDef "Placeholder"
            & opAttr "dtype" .~ tensorType (undefined :: a)
            & opAttr "shape" .~ shape'

-- | Creates a variable initialized to the given value.
-- Initialization happens next time session runs.
initializedVariable :: forall a . TensorType a
                    => Tensor Value a -> Build (Tensor Ref a)
initializedVariable initializer = do
    v <- CoreOps.variable []  -- The shape is not known initially.
    (i :: Tensor Ref a) <-
        buildOp (opDef "Assign"
                 & opAttr "T" .~ tensorType (undefined :: a)
                 & opAttr "use_locking" .~ True
                 & opAttr "validate_shape" .~ False
                 )
        v initializer
    addInitializer =<< group i
    return v

-- | Creates a zero-initialized variable with the given shape.
zeroInitializedVariable
  :: (TensorType a, Num a) =>
     TensorFlow.Types.Shape -> Build (Tensor TensorFlow.Tensor.Ref a)
zeroInitializedVariable = initializedVariable . zeros

-- TODO: Support heterogeneous list of tensors.
save :: forall a v . TensorType a
        => ByteString     -- ^ File path.
        -> [Tensor v a]  -- ^ Tensors to save.
        -> Build ControlNode
save path xs = do
    let toByteStringTensor = scalar . encodeUtf8 . unNodeName
    names <- mapM (fmap toByteStringTensor . renderNodeName) xs
    let types = replicate (length xs) (tensorType (undefined :: a))
    let saveOp = buildOp $ opDef "Save"
                         & opAttr "T" .~ types
    saveOp (scalar path) (CoreOps.pack names) xs

-- | Restore a tensor's value from a checkpoint file.
--
-- This version allows restoring from a checkpoint file that uses a different
-- tensor name than the variable.
restoreFromName :: forall a . TensorType a
                => ByteString    -- ^ File path.
                -> ByteString    -- ^ Tensor name override.
                -> Tensor Ref a  -- ^ Tensor to restore.
                -> Build ControlNode
restoreFromName path name x = do
    let restoreOp = buildOp $ opDef "Restore"
                            & opAttr "dt" .~ tensorType (undefined :: a)
    group =<< CoreOps.assign x
                (restoreOp (scalar path) (scalar name) :: Tensor Value a)

-- | Restore a tensor's value from a checkpoint file.
restore :: forall a . TensorType a
        => ByteString    -- ^ File path.
        -> Tensor Ref a  -- ^ Tensor to restore.
        -> Build ControlNode
restore path x = do
    name <- encodeUtf8 . unNodeName <$> renderNodeName x
    restoreFromName path name x

-- | Create a constant tensor.
--
-- The values should be in row major order, e.g.,
--
--   element 0:   index (0, ..., 0)
--   element 1:   index (0, ..., 1)
--   ...
constant :: forall a . TensorType a => Shape -> [a] -> Tensor Value a
constant (Shape shape') values
    | invalidLength = error invalidLengthMsg
    | otherwise = buildOp $ opDef "Const"
                          & opAttr "value" .~ typedNode
                          & opAttr "dtype" .~ nodeType
  where
    invalidLength = product shape' /= fromIntegral (length values)
    invalidLengthMsg = printf "invalid tensor length: expected %d got %d"
                              (product shape')
                              (length values)
    nodeType = tensorType (undefined :: a)
    typedNode :: TensorProto
    typedNode = def
                & dtype .~ nodeType
                & tensorShape.TensorShape.dim .~
                      [def & TensorShape.size .~ x | x <- shape']
                & tensorVal .~ values

-- | Reshape a N-D tensor down to a scalar.
--
-- See `TensorFlow.GenOps.Core.reshape`.
scalarize :: (TensorType a) => Tensor v a -> Tensor Value a
scalarize t = CoreOps.reshape t (vector scalarShape)
    where
        scalarShape = [] :: [Int32]

-- | Create a constant vector.
vector :: TensorType a => [a] -> Tensor Value a
vector xs = constant [fromIntegral $ length xs] xs

-- | Create a constant scalar.
scalar :: forall a . TensorType a => a -> Tensor Value a
scalar x = constant [] [x]

-- Random tensor from the unit normal distribution with bounded values.
truncatedNormal :: forall a v . TensorType a
                => Tensor v Int64  -- ^ Shape.
                -> Build (Tensor Value a)
truncatedNormal = buildOp $ opDef "TruncatedNormal"
                          & opAttr "dtype" .~ tensorType (undefined :: a)
                          & opAttr "T" .~ tensorType (undefined :: Int64)

zeros :: forall a . (Num a, TensorType a) => Shape -> Tensor Value a
zeros (Shape shape') = CoreOps.fill (vector $ map fromIntegral shape') (scalar 0)

shape :: (TensorType t) => Tensor v1 t -> Tensor Value Int32
shape = CoreOps.shape

expandDims :: (TensorType t) => Tensor v1 t -> Tensor v2 Int32 -> Tensor Value t
expandDims = CoreOps.expandDims

-- | Helper function for reduction ops (translation of math_ops.reduced_shape).
reducedShape :: (OneOf '[ Int32, Int64 ] t1, OneOf '[ Int32, Int64 ] t2) =>
                Tensor v1 t1 -> Tensor v2 t2 -> Tensor Value Int32
reducedShape inputShape axes =
    let inputShape32 = toInt32 inputShape         -- [2, 3, 5, 7]
        axes32 = toInt32 axes                     -- [1, 2]
        toInt32 x = CoreOps.cast x :: Tensor Value Int32
        inputRank = CoreOps.size inputShape32     -- 4
        axesMod = (axes32 + inputRank) `CoreOps.mod` inputRank
        axesShape = shape axesMod                 -- [2]
    in CoreOps.dynamicStitch                      -- [2, 1, 1, 7]
         [CoreOps.range 0 inputRank 1,            -- [0, 1, 2, 3]
           axesMod]                               -- [1, 2]
         [inputShape32,                           -- [2, 3, 5, 7]
           CoreOps.fill axesShape 1]              -- [1, 1]

max_stars_repo_name: blackgnezdo/tensorflow-haskell

max_issues_repo_name: blackgnezdo/tensorflow-haskell

max_forks_repo_name: blackgnezdo/tensorflow-haskell


text:
module Statistics (
    tests_Statistics
    ) where

import Control.Monad( replicateM )
import Data.AEq
import Data.List( foldl' )
import Debug.Trace
import Test.Framework
import Test.Framework.Providers.QuickCheck2
import Test.QuickCheck hiding ( vector )
import qualified Test.QuickCheck as QC

import Numeric.LinearAlgebra
import qualified Numeric.LinearAlgebra.Packed as P
import qualified Numeric.LinearAlgebra.Matrix as M
import qualified Numeric.LinearAlgebra.Vector as V

import Test.QuickCheck.LinearAlgebra( VectorList(..), WeightedVectorList(..) )
import qualified Test.QuickCheck.LinearAlgebra as Test

import Typed

testAEq actual expected =
    if actual ~== expected
        then True
        else trace ("expected: " ++ show expected ++ "\nactual: " ++ show actual)
                   False

tests_Statistics = testGroup "Statistics"
    [ testPropertyD "sum" prop_sum
    , testPropertyD "weightedSum" prop_weightedSum
    , testPropertyD "mean" prop_mean
    , testPropertyD "weightedMean" prop_weightedMean
    , testPropertyD "weightedMean (equal weights)" prop_weightedMean_eqw
    , testPropertyD "Matrix cov" prop_Matrix_cov
    , testPropertyD "Matrix covWithMean" prop_Matrix_covWithMean
    , testPropertyD "Matrix weightedCov (equal weights)" prop_Matrix_weightedCov_eqw
    , testPropertyD "Matrix weightedCov" prop_Matrix_weightedCov
    , testPropertyD "Matrix weightedCovWithMean" prop_Matrix_weightedCovWithMean
    , testPropertyD "Packed cov" prop_Packed_cov
    , testPropertyD "Packed covWithMean" prop_Packed_covWithMean
    , testPropertyD "Packed weightedCov (equal weights)" prop_Packed_weightedCov_eqw
    , testPropertyD "Packed weightedCov" prop_Packed_weightedCov
    , testPropertyD "Packed weightedCovWithMean" prop_Packed_weightedCovWithMean
    ]

prop_sum t (VectorList p xs) =
    V.sum p xs ~== foldl' V.add (V.constant p 0) xs
  where
    _ = typed t (head xs)

prop_weightedSum t (WeightedVectorList p wxs) =
    V.weightedSum p wxs
            ~== V.sum p (map (uncurry V.scale) wxs)
  where
    n = length wxs
    _ = typed t (snd $ head wxs)

prop_mean t (VectorList p xs) =
    V.mean p xs ~== V.scale (1/n) (V.sum p xs)
  where
    n = fromIntegral $ max (length xs) 1
    _ = typed t (head xs)

prop_weightedMean_eqw t (VectorList p xs) = let
    wxs = zip (repeat 1) xs
    in V.weightedMean p wxs ~== V.mean p xs
  where
    _ = typed t (head xs)

prop_weightedMean t (WeightedVectorList p wxs) = let
        w_sum = (sum . fst . unzip) wxs
        in V.weightedMean p wxs
            ~== if w_sum == 0 then V.constant p 0
                              else V.scale (1/w_sum) (V.weightedSum p wxs)
  where
    n = length wxs
    _ = typed t (snd $ head wxs)

prop_Matrix_cov t (VectorList p xs) =
    forAll (Test.vector p) $ \z ->
    forAll (elements [ UnbiasedCov, MLCov ]) $ \method -> let
        xbar = V.mean p xs
        ys = [ V.sub x xbar | x <- xs ]
        scale = case method of { UnbiasedCov -> 1/(n-1) ; MLCov -> 1/n }
        cov' = foldl' (flip $ \y -> M.rank1Update scale y y)
                      (M.constant (p,p) 0)
                      ys
        cov = M.cov p method xs

        in M.hermMulVector cov z ~== M.mulVector NoTrans cov' z
  where
    n = fromIntegral $ length xs
    _ = typed t $ head xs

prop_Matrix_covWithMean t (VectorList p xs) =
    forAll (Test.vector p) $ \z ->
    forAll (elements [ UnbiasedCov, MLCov ]) $ \method -> let
        xbar = V.mean p xs
        cov' = M.covWithMean xbar method xs
        cov = M.cov p method xs
        in M.hermMulVector cov z ~== M.hermMulVector cov' z
  where
    n = fromIntegral $ length xs
    _ = typed t $ head xs

prop_Matrix_weightedCov_eqw t (VectorList p xs) =
    forAll (Test.vector p) $ \z ->
    forAll (elements [ UnbiasedCov, MLCov ]) $ \method -> let
        wxs = zip (repeat 1) xs
        cov = M.weightedCov p method wxs
        cov' = M.cov p method xs
        in M.hermMulVector cov z ~== M.hermMulVector cov' z
  where
    n = fromIntegral $ length xs
    _ = typed t $ head xs

prop_Matrix_weightedCov t (WeightedVectorList p wxs) =
    forAll (Test.vector p) $ \z ->
    forAll (elements [ UnbiasedCov, MLCov ]) $ \method -> let
        (ws,xs) = unzip wxs

        w_sum = sum ws
        ws' = [ w / w_sum | w <- ws ]
        w2_sum = sum [ w*w | w <- ws' ]
        scale = case method of { UnbiasedCov -> 1/(1-w2_sum) ; MLCov -> 1 }

        xbar = V.weightedMean p wxs
        wys = [ (w, V.sub x xbar) | (w,x) <- zip ws' xs ]
        cov' = if w_sum == 0
                    then M.constant (p,p) 0
                    else foldl' (flip $ \(w,y) -> M.rank1Update (scale*w) y y)
                                (M.constant (p,p) 0)
                                wys

        cov = M.weightedCov p method wxs

        in M.hermMulVector cov z ~== M.mulVector NoTrans cov' z
  where
    n = fromIntegral $ length wxs
    _ = typed t $ snd $ head wxs

prop_Matrix_weightedCovWithMean t (WeightedVectorList p wxs) =
    forAll (Test.vector p) $ \z ->
    forAll (elements [ UnbiasedCov, MLCov ]) $ \method -> let
        xbar = V.weightedMean p wxs
        cov' = M.weightedCov p method wxs
        cov = M.weightedCovWithMean xbar method wxs
        in M.hermMulVector cov z ~== M.hermMulVector cov' z
  where
    n = fromIntegral $ length wxs
    _ = typed t $ snd $ head wxs

prop_Packed_cov t (VectorList p xs) =
    forAll (Test.vector p) $ \z ->
    forAll (elements [ UnbiasedCov, MLCov ]) $ \method -> let
        cov' = M.cov p method xs
        cov = P.cov p method xs
        in P.hermMulVector cov z ~== M.hermMulVector cov' z
  where
    n = fromIntegral $ length xs
    _ = typed t $ head xs

prop_Packed_covWithMean t (VectorList p xs) =
    forAll (Test.vector p) $ \z ->
    forAll (elements [ UnbiasedCov, MLCov ]) $ \method -> let
        xbar = V.mean p xs
        cov = P.covWithMean xbar method xs
        cov' = P.cov p method xs
        in P.hermMulVector cov z ~== P.hermMulVector cov' z
  where
    n = fromIntegral $ length xs
    _ = typed t $ head xs

prop_Packed_weightedCov_eqw t (VectorList p xs) =
    forAll (Test.vector p) $ \z ->
    forAll (elements [ UnbiasedCov, MLCov ]) $ \method -> let
        wxs = zip (repeat 1) xs
        cov = P.weightedCov p method wxs
        cov' = P.cov p method xs
        in P.hermMulVector cov z ~== P.hermMulVector cov' z
  where
    n = fromIntegral $ length xs
    _ = typed t $ head xs

prop_Packed_weightedCov t (WeightedVectorList p wxs) =
    forAll (Test.vector p) $ \z ->
    forAll (elements [ UnbiasedCov, MLCov ]) $ \method -> let
        cov' = M.weightedCov p method wxs
        cov = P.weightedCov p method wxs
        in P.hermMulVector cov z ~== M.hermMulVector cov' z
  where
    n = fromIntegral $ length wxs
    _ = typed t $ snd $ head wxs

prop_Packed_weightedCovWithMean t (WeightedVectorList p wxs) =
    forAll (Test.vector p) $ \z ->
    forAll (elements [ UnbiasedCov, MLCov ]) $ \method -> let
        xbar = V.weightedMean p wxs
        cov' = P.weightedCov p method wxs
        cov = P.weightedCovWithMean xbar method wxs
        in P.hermMulVector cov z ~== P.hermMulVector cov' z
  where
    n = fromIntegral $ length wxs
    _ = typed t $ snd $ head wxs

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
{-# LANGUAGE DeriveAnyClass    #-}
{-# LANGUAGE DeriveGeneric     #-}
{-# LANGUAGE GADTs             #-}
{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE TypeFamilies      #-}
-- Ops.hs ---
--
-- Filename: Ops.hs
-- Description:
-- Author: Manuel Schneckenreither
-- Maintainer:
-- Created: Wed Nov  1 17:40:44 2017 (+0100)
-- Version:
-- Package-Requires: ()
-- Last-Updated:
--           By:
--     Update #: 155
-- URL:
-- Doc URL:
-- Keywords:
-- Compatibility:
--
--

-- Commentary:
--
--
--
--

-- Change Log:
--
--
--
--
--
--
--

-- Code:

module SimSim.Order.Ops
  ( generateOrdersContDueDateDistr
  , generateOrdersDiscDueDateDistr
  , generateOrdersFixedDueDateSlack
  , generateOrdersUniform
  ) where

import           ClassyPrelude
import           Control.DeepSeq
import           Data.Serialize
import           GHC.Generics
import           Statistics.Distribution
import           Statistics.Distribution.DiscreteUniform
import           Statistics.Distribution.Uniform
import           System.Random.MWC

import           SimSim.Order.Type
import           SimSim.ProductType
import           SimSim.Simulation.Type
import           SimSim.Time

type MinimumTime = Time

generateOrdersContDueDateDistr ::
     (ContGen interarrivalTime, DiscreteGen productType, ContGen dueDateSlack) => SimSim -> interarrivalTime -> productType -> (MinimumTime, dueDateSlack) -> IO [Order]
generateOrdersContDueDateDistr sim dInterarrival dProductType (minTime, dueDateSlack) = do
  let repeatUntilTime acc = do
        interArr <- genContVar dInterarrival gen
        if acc + interArr < periodLen
          then ((acc + interArr) :) <$> repeatUntilTime (acc + interArr)
          else return []
  arrivals <- map ((+currentTime) . timeFromDouble) <$> repeatUntilTime 0
  productTypes <- map (Product . (+ 1) . (`mod` length pts)) <$> mapM (\_ -> genDiscreteVar dProductType gen) arrivals
  dueDates <- mapM (\a -> (\x -> timeFromDouble x + a + minTime) <$> genContVar dueDateSlack gen) arrivals
  return $ zipWith3 newOrder productTypes arrivals dueDates
  where
    currentTime = simCurrentTime sim
    gen = simRandGenDemand $ simInternal sim
    periodLen = timeToDouble $ simPeriodLength sim
    pts = productTypes sim

generateOrdersDiscDueDateDistr ::
     (ContGen interarrivalTime, DiscreteGen productType, DiscreteGen dueDateSlack) => SimSim -> interarrivalTime -> productType -> (MinimumTime, dueDateSlack) -> IO [Order]
generateOrdersDiscDueDateDistr sim dInterarrival dProductType (minTime, dueDateSlack) = do
  let repeatUntilTime acc = do
        interArr <- genContVar dInterarrival gen
        if acc + interArr < periodLen
          then ((acc + interArr) :) <$> repeatUntilTime (acc + interArr)
          else return []
  arrivals <- map ((+currentTime) . timeFromDouble) <$> repeatUntilTime 0
  productTypes <- map (Product . (+ 1) . (`mod` length pts)) <$> mapM (\_ -> genDiscreteVar dProductType gen) arrivals
  dueDates <- mapM (\_ -> (\x -> x + minTime + currentTime) . Time . fromIntegral <$> genDiscreteVar dueDateSlack gen) arrivals
  return $ zipWith3 newOrder productTypes arrivals dueDates
  where
    currentTime = simCurrentTime sim
    gen = simRandGenDemand $ simInternal sim
    periodLen = timeToDouble $ simPeriodLength sim
    pts = productTypes sim

generateOrdersFixedDueDateSlack :: (ContGen interarrivalTime, ContGen productType) => SimSim -> interarrivalTime -> productType -> Time -> IO [Order]
generateOrdersFixedDueDateSlack sim dInterarrival dProductType dueDateSlack = do
  let repeatUntilTime acc = do
        interArr <- genContVar dInterarrival gen
        if acc + interArr < periodLen
          then ((acc + interArr) :) <$> repeatUntilTime (acc + interArr)
          else return []
  arrivals <- map ((+currentTimeMinusPeriod) . timeFromDouble) <$> repeatUntilTime 0
  -- print currentTimeMinusPeriod
  -- print $ map timeToDouble arrivals
  -- productTypes <- map (Product . (+ 1) . (`mod` length pts)) <$> mapM (\_ -> genDiscreteVar dProductType gen) arrivals
  productTypes <- map (Product . (+ 1) . (`mod` length pts) . round) <$> mapM (\_ -> genContVar dProductType gen) arrivals
  let dueDates = replicate (length arrivals) (currentTime + dueDateSlack)
  return $ zipWith3 newOrder productTypes arrivals dueDates
  where
    currentTimeMinusPeriod = currentTime - timeFromDouble periodLen
    currentTime = simCurrentTime sim
    gen = simRandGenDemand $ simInternal sim
    periodLen = timeToDouble $ simPeriodLength sim
    pts = productTypes sim

generateOrdersUniform :: SimSim -> Int -> Int -> Time -> IO [Order]
generateOrdersUniform sim minOrders maxOrders dueDateSlack = do
  nrOrders <- round <$> genContVar (uniformDistr (fromIntegral minOrders-0.4999) (fromIntegral maxOrders+0.4999)) gen
  let arrivals = replicate nrOrders currentTime
  -- productTypes <- map Product <$> mapM (\_ -> genDiscreteVar (discreteUniform (length pts)) gen) arrivals
  productTypes <- map (Product . round) <$> mapM (\_ -> genContVar (uniformDistr 0.5001 (fromIntegral (length pts)+0.4999)) gen) arrivals
  let dueDates = replicate (length arrivals) (currentTime + dueDateSlack)
  return $ zipWith3 newOrder productTypes arrivals dueDates
  where
    currentTime = simCurrentTime sim
    gen = simRandGenDemand $ simInternal sim
    pts = productTypes sim

--
-- Ops.hs ends here

max_stars_repo_name: schnecki/simsim

max_issues_repo_name: schnecki/simsim

max_forks_repo_name: schnecki/simsim