text:
{-# LANGUAGE CPP                 #-}
{-# LANGUAGE DeriveDataTypeable  #-}
{-# LANGUAGE DeriveGeneric       #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications    #-}
-- | Tests for forms.
module FormTest where

#include "ghc_modules.h"

-- base
import           Control.Applicative          (Alternative (..))
import           Data.Char                    (toUpper)
import           Data.Complex
import           Data.Data
import qualified Data.Fixed                   as Fixed
import           Data.Functor.Compose
import           Data.Functor.Const
import           Data.Functor.Identity
import qualified Data.Functor.Product         as Product
import qualified Data.Functor.Sum             as Sum
import           Data.Int
import           Data.List                    (isPrefixOf, isSubsequenceOf)
import           Data.List.NonEmpty           (NonEmpty (..))
import           Data.Monoid
import           Data.Ratio
import qualified Data.Semigroup               as Semigroup
import           Data.Version
import           Data.Word
import           GHC.Generics                 (Generic (..))
import           Numeric.Natural
import           Text.Show.Functions          ()

#if !MIN_VERSION_ghc(8,4,0)
import           Data.Monoid                  ((<>))
#endif

-- binary
import           Data.Binary                  (decode, encode)

-- deepseq
import           Control.DeepSeq

-- ghc
import           GHC_Data_FastString          (fsLit, unpackFS)
import           GHC_Data_StringBuffer        (stringToStringBuffer)
import           GHC_Driver_Session           (HasDynFlags (..))
import           GHC_Types_SrcLoc             (GenLocated (..), SrcSpan (..),
                                               noSrcSpan)
import           GHC_Driver_Ppr               (showPpr)

#if MIN_VERSION_ghc(9,0,0)
import           GHC_Types_SrcLoc             (UnhelpfulSpanReason (..))
#endif

-- transformers
import           Control.Monad.Trans.State

-- hspec
import           Test.Hspec

-- QuickCheck
import           Test.QuickCheck

-- finkel-kernel
import           Language.Finkel.Fnk
import           Language.Finkel.Form
import           Language.Finkel.Homoiconic
import qualified Language.Finkel.Homoiconic   as Homoiconic
import           Language.Finkel.Lexer
import           Language.Finkel.Reader
import           Language.Finkel.SpecialForms

-- Internal
import           Orphan                       ()

formTests :: Spec
formTests = do
  mapM_ readShow
    [ "foo", "#'a", "12345", "6.789", "0.001"
    , "(foo bar buzz)"
    , "(#'a #'\\SP #'\\ \"bcd\")"
    , "[#'\\BEL #'\\BS #'\\FF #'\\LF #'\\CR #'\\HT #'\\VT]"
    , "[()]"
    , "(1 -2 345 6.789 0.001)" ]

  readUnicodeStringProp
  readShowFormProp

  dataInstanceTests
  qFunctionTests

  fracTest 1.23
  fracTest (-1.23)
  fracTest 0
  fracTest 1e-9

  showTest
  pprTest
  functorTest "(a \"foo\" \\x [True False])"
  applicativeTest
  monadTest
  foldableTest
  traversableTest
  binaryTest

  eqTest "(a \"bcd\" \\e [f g] (h i))"
  eqPropTest

  locationTest Nothing "foo"
  locationTest (Just "locationTest") "foo"
  asLocOfTest

  lengthTest 3 "(a b c)"
  lengthTest 5 "(a (b (c)) d e)"
  lengthTest 1 "()"
  lengthTest 8 "[a (b (c d e) [f g]) h]"
  lengthTest 1 "foo"

  homoiconicTests
  rnfTest
  listTest
  numTest
  fractionalTest
  monoidTest
  alternativeTest

  fromCodeTest
  fromCodeErrorTest

  dataToCodeTest
  genericHomoiconicTest

  unquoteSpliceTest

readShow :: String -> Spec
readShow str =
  describe ("read and show `" ++ str ++ "'") $
    it "should match the input" $
      show (parseE str) `shouldBe` str

readUnicodeStringProp :: Spec
readUnicodeStringProp =
  describe "read and show unicode string" $
    it "should return itself" $
      property (\uni ->
                  let str = getUnicodeString uni
                  in  parseE (show str) == toCode (aString NoSourceText str))

readShowFormProp :: Spec
readShowFormProp =
  describe "read and show form property" $
    it "should match the input" $
      property (\form ->
                  form == form && parseE (show form) `eqForm` form)

dataInstanceTests :: Spec
dataInstanceTests = do
  let gfoldl_self atom =
        gfoldl (\mb_f x -> case mb_f of
                   Just f -> pure (f x)
                   _ -> error "should not happen")
               return atom
      t_gfoldl_self x = gfoldl_self x `shouldBe` Just x
      t_show_constr x y = show (toConstr x) `shouldBe` y
  describe "Data instance for Atom" $ do
    let aunit = AUnit
        asym = ASymbol (fsLit "foo")
        achar = AChar NoSourceText 'a'
        astr = AString NoSourceText (fsLit "string")
        aint = AInteger (mkIntegralLit (42 :: Int))
        afrac = aFractional (1.23 :: Double)
    it "should return Just self with simple gfoldl" $ do
      t_gfoldl_self aunit
      t_gfoldl_self asym
      t_gfoldl_self achar
      t_gfoldl_self astr
      t_gfoldl_self aint
      t_gfoldl_self afrac
    it "should show itself with toConstr" $ do
      t_show_constr aunit "AUnit"
      t_show_constr asym "ASymbol"
      t_show_constr achar "AChar"
      t_show_constr astr "AString"
      t_show_constr aint "AInteger"
      t_show_constr afrac "AFractional"
    it "should return AUnit with simple gunfold" $ do
      gunfold (const Nothing) Just (toConstr AUnit) `shouldBe` Just AUnit
    it "should return AUnit constr" $ do
      let dtype = dataTypeOf AUnit
          cnstr = toConstr AUnit
      readConstr dtype "AUnit" `shouldBe` Just cnstr
  describe "Data instance for Form" $ do
    let fatom = Atom AUnit
        qc c = qChar c "" 0 0 0 0
        flist = List [qc 'a', qc 'b']
        fhslist = HsList [qc 'a', qc 'b']
        ftend :: Form Atom
        ftend = TEnd
    it "should return Just self with simple gfoldl" $ do
      t_gfoldl_self fatom
      t_gfoldl_self flist
      t_gfoldl_self fhslist
      t_gfoldl_self ftend
    it "should show itself with toConstr" $ do
      t_show_constr fatom "Atom"
      t_show_constr flist "List"
      t_show_constr fhslist "HsList"
      t_show_constr ftend "TEnd"
    it "should return TEnd with simple gunfold" $
      gunfold (const Nothing) Just (toConstr ftend) `shouldBe` Just ftend
    it "should return Atom constr" $ do
      let dtype = dataTypeOf fatom
          cnstr = toConstr fatom
      readConstr dtype "Atom" `shouldBe` Just cnstr
    it "should return same result from dataCast1 and gcast1" $ do
      (dataCast1 [TEnd] :: Maybe [Form Atom]) `shouldBe` Just [TEnd]
  describe "Data instance for LForm" $ do
    let qc = qChar 'x' "" 0 0 0 0
        d1, d2, d3 :: Data a => a
        d1 = fromConstr (toConstr noSrcSpan)
        d2 = fromConstrB (fromConstr (toConstr AUnit))
                         (toConstr (Atom AUnit))
        d3 = evalState m 0
          where
            m = fromConstrM act (toConstr (L noSrcSpan (Atom AUnit)))
            act :: Data d => State Int d
            act = do
              i <- get
              modify succ
              case i of
                0 -> return d1
                1 -> return d2
                _ -> error ("index " ++ show i ++ " for LForm")
        d4 = fromConstrB d3 (toConstr qu)
        gc1 :: Data a => Maybe [LForm a]
        gc1 =
          let a :: Data a => a
              a = fromConstr (toConstr AUnit)
          in  dataCast1 [LForm (L noSrcSpan (Atom a))]
        qu = qUnit "" 0 0 0 0
    it "should return Just self with simple gfoldl" $ do
      t_gfoldl_self qc
    it "should show itself with toConstr" $ do
      t_show_constr qc "LForm"
    it "should return LForm constr" $ do
      let dtype = dataTypeOf qc
          cnstr = toConstr qc
      readConstr dtype "LForm" `shouldBe` Just cnstr
    it "should construct qUnit from constructors" $ do
      d4 `shouldBe` qu
    it "should return qUnit from dataCast1" $
      gc1 `shouldBe` Just [qu]

qFunctionTests :: Spec
qFunctionTests = do
  describe "qSymbol function" $
    it "should equal to quoted symbol" $
      qSymbol "foo" "" 0 0 0 0 `shouldBe` toCode (ASymbol (fsLit "foo"))
  describe "qChar function" $
    it "should equal to quoted char" $
      qChar 'x' "" 0 0 0 0 `shouldBe` toCode 'x'
  describe "qString function" $
    it "should equal to quoted string" $
      qString "foo" "" 0 0 0 0 `shouldBe` toCode "foo"
  describe "qInteger function" $
    it "should equal to quoted integer" $
      qInteger 42 "" 0 0 0 0 `shouldBe` toCode (42 :: Integer)
  describe "qFractional function" $
    it "should equal to quoted fractional" $
      qFractional (1.23 :: Double) "" 0 0 0 0
        `shouldBe` toCode (1.23 :: Double)
  describe "qUnit function" $
    it "should equal to quoted unit" $
      qUnit "" 0 0 0 0 `shouldBe` toCode ()

  let qc x = qChar x "" 0 0 0 0
  describe "qList function" $
    it "should equal to quoted codes" $
      let xs = [qc 'a', qc 'b']
      in qList xs "" 0 0 0 0 `shouldBe` toCode (List xs)
  describe "qHsList function" $
    it "should equal to quoted haskell list" $
      let xs = [qc 'a', qc 'b']
      in  qHsList xs "" 0 0 0 0 `shouldBe` toCode (HsList xs)

fracTest :: Double -> Spec
fracTest x =
  describe ("read and show a fractional number `" ++ show x ++ "'") $
    it "should match the input" $
       show (aFractional x) `shouldBe` show x

showTest :: Spec
showTest =
  describe "showing TEnd" $
    it "should be \"TEnd\"" $
      show (TEnd :: Form Atom) `shouldBe` "TEnd"

functorTest :: String -> Spec
functorTest str = do
  describe ("Functor instance of Code `" ++ str ++ "'") $
    it "should obey the Functor law" $
       let c = parseE str
       in  fmap id c `shouldBe` c
  describe "fmap to TEnd" $
    it "should be TEnd" $ do
      let te :: Form Atom
          te = TEnd
          f :: Atom -> Atom
          f _ = AUnit
      fmap f te `shouldBe` te

applicativeTest :: Spec
applicativeTest = do
  let atom_a = AChar NoSourceText 'a'
      char_a = toCode atom_a
      al2 = qList [char_a, char_a] "" 0 0 0 0
      ahl2 = qHsList [char_a, char_a] "" 0 0 0 0
      unit = toCode ()
      f1 a b = (a,b)
      a_pair = lf (Atom (atom_a, atom_a))
      a_pair_ls = mk_a_pairs List 2 -- lf (List [a_pair, a_pair])
      a_pair_hls = mk_a_pairs HsList 2 -- lf (HsList [a_pair, a_pair])
      mk_a_pairs f n = lf (f (replicate n a_pair))
      lf = LForm . genSrc
      tend :: Form Atom
      tend = TEnd

  describe "pure" $
    it "should result in atom" $
      toCode (pure AUnit :: Form Atom) `shouldBe` unit

  describe "<*>" $ do
    it "should apply f1 to atom and atom" $
      f1 <$> char_a <*> char_a `shouldBe` a_pair
    it "should apply f1 to atom and list" $
      f1 <$> char_a <*> al2 `shouldBe` a_pair_ls
    it "should apply f1 to atom and hslist" $
      f1 <$> char_a <*> ahl2 `shouldBe` a_pair_hls
    it "should aply f1 to list and atom" $
      f1 <$> al2 <*> char_a `shouldBe` a_pair_ls
    it "should apply f1 to list and list" $
      f1 <$> al2 <*> al2 `shouldBe` mk_a_pairs List 4
    it "should apply f1 to list and hslist" $
      f1 <$> al2 <*> ahl2 `shouldBe` mk_a_pairs List 4
    it "should aply f1 to list and atom" $
      f1 <$> ahl2 <*> char_a `shouldBe` a_pair_hls
    it "should apply f1 to list and list" $
      f1 <$> ahl2 <*> al2 `shouldBe` mk_a_pairs HsList 4
    it "should apply f1 to list and list" $
      f1 <$> ahl2 <*> ahl2 `shouldBe` mk_a_pairs HsList 4
    it "should result in TEnd" $ do
      f1 <$> char_a <*> lf tend `shouldBe` lf TEnd
      f1 <$> lf tend <*> char_a `shouldBe` lf TEnd

monadTest :: Spec
monadTest = do
  let f1 x = case x of
               AChar st c -> AChar st (toUpper c)
               _          -> x
      qh x = qHsList x "" 0 0 0 0
      ql x = qList x "" 0 0 0 0
  describe "bind" $ do
    it "should apply f1 to atom" $
      do {x <- toCode 'x'; return (f1 x)} `shouldBe` toCode 'X'
    it "should apply f1 to list" $
      do {x <- ql [toCode 'x', toCode 'x']; return (f1 x)} `shouldBe`
         ql [toCode 'X', toCode 'X']
    it "should apply f1 to hslist" $
      do {x <- qh [toCode 'x',toCode 'x']; return (f1 x)} `shouldBe`
         qh [toCode 'X', toCode 'X']
    it "should apply f1 to TEnd" $
      do {x <- toCode (TEnd :: Form Atom); return (f1 x)} `shouldBe`
         toCode (TEnd :: Form Atom)

foldableTest :: Spec
foldableTest = do
  let fsum = foldr (\x acc -> case x of
                                AInteger il -> acc + il_value il
                                _           -> acc)
                   0
  describe "taking sum of 1 to 10 with foldr" $ do
    let str1 = "(1 2 3 4 5 6 7 8 9 10)"
    it ("should be 55 for " ++ str1) $
      fsum (parseE str1) `shouldBe` 55
    it "should be 0 for TEnd" $ do
#if MIN_VERSION_ghc(9,0,0)
      let sp = UnhelpfulSpan (UnhelpfulOther (fsLit "<foldableTest>"))
#else
      let sp = UnhelpfulSpan (fsLit "<foldableTest>")
#endif
      fsum (LForm (L sp TEnd)) `shouldBe` 0
  describe "length of nil" $
    it "should be 0" $
      length nil `shouldBe` 0

traversableTest :: Spec
traversableTest = do
  let f atom = case atom of
                 ASymbol sym | '$' : _ <- unpackFS sym
                             -> return (aSymbol "_")
                 _ -> return atom
  describe "replacing symbol with mapM" $ do
    let str1 = "(a $b c [$d e] [f g $h] ($i $j))"
        str2 = "(a  _ c [ _ e] [f g  _] ( _  _))"
    it "should replace $.* with _" $ do
      let form1 = parseE str1
          form2 = parseE str2
      mapM f form1 `shouldBe` Just form2
  describe "traversing TEnd" $
    it "should be Just TEnd" $
      mapM f TEnd `shouldBe` Just TEnd

binaryTest :: Spec
binaryTest = do
  describe "Instances of Get and Put from binary" $ do
    it "should return the original value" $ do
      let to_from_bin :: Code -> Bool
          to_from_bin x = decode (encode x) == x
      property to_from_bin

eqTest :: String -> Spec
eqTest str =
  describe "parsing same string twice" $
   it "should result in equal codes" $
     let c1 = parseE str
         c2 = parseE str
     in  c1 `shouldBe` c2

eqPropTest :: Spec
eqPropTest = do
  describe "Eq instance for LForm" $
    it "should ignore location information" $ do
      let g :: Code -> Bool
          g x@(LForm (L _ body)) = x == LForm (L sp body)
#if MIN_VERSION_ghc(9,0,0)
          sp = UnhelpfulSpan (UnhelpfulOther (fsLit "<eqPropTest>"))
#else
          sp = UnhelpfulSpan (fsLit "<eqPropTest>")
#endif
      property g
  describe "comparing TEnd with TEnd" $
    it "should be True" $
      (TEnd :: Form Atom) `shouldBe` (TEnd :: Form Atom)

locationTest :: Maybe FilePath -> String -> Spec
locationTest mb_path str =
  describe ("location of `" ++ str ++ "'") $
    it (case mb_path of
          Just path -> "should contain `" ++ path ++ "'"
          Nothing   -> "should be unhelpful") $ do
       let c = parseE' mb_path str
           l = showLoc c
       case mb_path of
         Just path -> (path `isPrefixOf` l) `shouldBe` True
         Nothing   -> ("anon" `isSubsequenceOf` l) `shouldBe` True

asLocOfTest :: Spec
asLocOfTest =
  describe ("apply asLocOf to code") $ do
    it "should return the value of first arg" $
      property (\ a b -> asLocOf a b == a)
    it "should use the location of second arg" $
      property (\ a b -> case (asLocOf a b, b) of
                   (LForm (L l1 _), LForm (L l2 _)) -> l1 == l2)

lengthTest :: Int -> String -> Spec
lengthTest n str =
  describe ("length of " ++ str) $
   it ("should be " ++ show n) $
     length (parseE str) `shouldBe` n

homoiconicTests :: Spec
homoiconicTests = do
  let t x = describe ("to/from code " ++ show x) $
              it "shoult match the input" $
                 case fromCode (toCode x) of
                   Just y  -> y `shouldBe` x
                   Nothing -> error ("got Nothing with " ++ show x)
  t (aIntegral (42 :: Int))
  t ()
  t 'x'
  t "string"
  t (42 :: Int)
  t (42 :: Int8)
  t (42 :: Int16)
  t (42 :: Int32)
  t (42 :: Int64)
  t (42 :: Integer)
  t (42 :: Word)
  t (42 :: Word8)
  t (42 :: Word16)
  t (42 :: Word32)
  t (42 :: Word64)
  t (0.123456789 :: Double)
  t (1.234 :: Float)
  t ([1,2,3] :: [Int])
  t (Fixed.MkFixed 2 :: Fixed.Pico)
  t (Identity 'a')
  t (1 :+ 2 :: Complex Int)
  t (Compose (Just (Just 'a')))
  t (Const True :: Const Bool Char)
  t (Product.Pair (Just 'a') (Just 'b'))
  t [Sum.InL (Just 'a'), Sum.InR (Right 'b'), Sum.InR (Left "foo")]
  t ('a' :| ['b', 'c', 'd'])
  t (All False)
  t (Alt (Just True))
  t (Any False)
  t (Dual 'x')
  t (First (Just 'a'))
  t (Last (Just 'a'))
  t (Product (42 :: Int))
  t (Sum (42 :: Int))
  t (Proxy :: Proxy ())
  t (Version [1,2,3] ["foo", "bar"])
  t (1 % 3 :: Rational)
  t (Semigroup.Arg 'x' False)
  t (Semigroup.First 'a')
  t (Semigroup.Last 'z')
  t (Semigroup.Max (42 :: Int))
  t (Semigroup.Min (42 :: Int))
#if !MIN_VERSION_ghc(9,0,0)
  t (Semigroup.Option (Just "foo"))
#endif
  t (Semigroup.WrapMonoid True)
  t (42 :: Natural)
  t (Atom (aIntegral (42 :: Int)))
  t (LForm (genSrc (Atom (aIntegral (42 :: Int)))))
  t [True, False]
  t [EQ, LT, GT]
  t (Just (42 :: Int))
  t [Right True, Left "foo"]
  t (Just 'x', [Right False, Left "foo"])
  t (Just 'x', [Right False, Left "foo"], EQ)
  t (Just 'x', [Right False, Left "foo"], EQ, 42::Int)
  t (Just 'x', [Right False, Left "foo"], EQ, 42::Int ,False)
  t (Just 'x', [Right False, Left "foo"], EQ, 42::Int
    ,False, Just [Right (Just EQ), Left (3.1 :: Double)])

rnfTest :: Spec
rnfTest = do
  describe "rnf of arbitrary form" $
    it "should return ()" $
       property (rnf :: Code -> ())
  describe "rnf of TEnd" $
    it "should return ()" $
      rnf (TEnd :: Form Atom) `shouldBe` ()

listTest :: Spec
listTest =
  describe "list from arbitrary form applied to arbitrary function" $
    it "should be a list" $ do
      let f :: (Code -> Code) -> Code -> Bool
          f g form = isListL (toListL (g form))
      property f

cInt :: Int -> Code
cInt = toCode

cDouble :: Double -> Code
cDouble = toCode

numTest :: Spec
numTest =
  describe "Num instance for Code" $ do
    it "should evaluate +" $ do
      cInt 2 + cInt 40 `shouldBe` 42
      cInt 2 + cDouble 40 `shouldBe` 42.0
      cDouble 2 + cInt 40 `shouldBe` 42.0
      cDouble 2.0 + cDouble 40.0 `shouldBe` 42.0
    it "should evaluate * for AInteger" $ do
      cInt 6 * cInt 7  `shouldBe` 42
      cDouble 6 * cDouble 7 `shouldBe` 42.0
    it "should evalue - for AInteger" $ do
      cInt 50 - cInt 8 `shouldBe` 42
      cDouble 50 - cDouble 8 `shouldBe` 42.0
    it "should evaluate abs" $ do
      abs (cInt (-42)) `shouldBe` 42
      abs (cDouble (-42)) `shouldBe` 42.0
    it "should evaluate signum" $ do
      signum (cInt (-42)) `shouldBe` -1
      signum (cDouble (-42)) `shouldBe` -1.0
    it "should evaluate fromInteger" $
      fromInteger 42 `shouldBe` cInt 42
    it "should result to nil with invalid values" $ do
      3 + toCode 'a' `shouldBe` nil
      signum (toCode 'a') `shouldBe` nil

fractionalTest :: Spec
fractionalTest =
  describe "Fractional instance for Code" $ do
    it "should evaluate /" $ do
      cInt 84 / cInt 2 `shouldBe` 42.0
      cDouble 84.0 / cDouble 2.0 `shouldBe` 42.0
    it "should evaluate recip" $ do
      recip (cInt 4) `shouldBe` 0.25
      recip (cDouble 4.0) `shouldBe` 0.25

monoidTest :: Spec
monoidTest = do
  let a = toCode 'a'
      b = toCode 'b'
      lst = toCode (List [a,b])
      hslst = toCode (HsList [a,b])
  describe "Atom <> XXX" $ do
    it "should result in '(a b)" $
      a <> b `shouldBe` lst
    it "should result in '(a a b)" $
      a <> lst `shouldBe` toCode (List [a,a,b])
    it "should result in '(a a b)" $
      a <> hslst `shouldBe` toCode (List [a,a,b])
  describe "List <> XXX" $ do
    it "should result in '(a b a)" $
      lst <> a `shouldBe` toCode (List [a,b,a])
    it "should result in '(a b a b)" $
      lst <> lst `shouldBe` toCode (List [a,b,a,b])
    it "should result in '(a b a b)" $
      lst <> hslst `shouldBe` toCode (List [a,b,a,b])
  describe "HsList <> XXX" $ do
    it "should result in '(a b a)" $
      hslst <> a `shouldBe` toCode (List [a,b,a])
    it "should result in '(a b a b)" $
      hslst <> lst `shouldBe` toCode (List [a,b,a,b])
    it "should result in '(a b a b)" $
      hslst <> hslst `shouldBe` toCode (List [a,b,a,b])
  describe "TEnd" $ do
    it "should result in non-TEnd" $ do
      let at = Atom True
      at <> TEnd `shouldBe` at
      TEnd <> at `shouldBe` at
  describe "mempty" $ do
    it "should be empty list" $
      mempty `shouldBe` (List [] :: Form Atom)
    it "should be nil" $
      mempty `shouldBe` nil

alternativeTest :: Spec
alternativeTest =
  describe "Alternative" $ do
    describe "empty" $ do
      it "should be nil for Code" $
        empty `shouldBe` nil
      it "should be empty form for Form" $
        empty `shouldBe` (mempty :: Form Atom)
    describe "<|>" $ do
      let a = toCode 'a'
          b = toCode 'b'
          c = toCode 'c'
      it "should append elements for Code" $
        let ab = toCode (List [a,b])
            cb = toCode (List [c,b])
        in  ab <|> cb `shouldBe` toCode (List [a,b,c,b])
      it "should append elements for Form Atom" $
        let ab = List [a,b]
            cb = List [c,b]
        in  ab <|> cb `shouldBe` List [a,b,c,b]

data Foo = Foo deriving (Eq, Show)

instance Homoiconic Foo where
  toCode foo = toCode (aSymbol (show foo))
  parseCode _ = Homoiconic.Failure "Foo"

fromCodeTest :: Spec
fromCodeTest = do
  describe "default toCode implementation" $
    it "should return Nothing" $
      (fromCode nil :: Maybe Foo) `shouldBe` Nothing

  describe "getting Nothing from fromCode" $ do
    let ng title x = it title $ x `shouldBe` Nothing
        q x = qSymbol x "" 0 0 0 0
        foo = q "foo"
        foo1 = qList [q "Foo", q "a"] "" 0 0 0 0
        foo2 = qList [q "Foo", q "a", q "b"] "" 0 0 0 0

    it "should result to Nothing with explicit Nothing" $ do
      fromCode (q "Nothing") `shouldBe` (Just Nothing :: Maybe (Maybe ()))
      toCode (Nothing :: Maybe ()) `shouldBe` q "Nothing"

    ng "()" (fromCode foo :: Maybe ())
    ng "Int" (fromCode foo :: Maybe Int)
    ng "Char" (fromCode nil :: Maybe Char)
    ng "[Int]" (fromCode nil :: Maybe [Int])
    ng "[Char]" (fromCode foo :: Maybe [Char])
    ng "Bool" (fromCode foo :: Maybe Bool)
    ng "Ordering" (fromCode foo :: Maybe Ordering)
    ng "Maybe ()" (fromCode foo1 :: Maybe (Maybe ()))
    ng "Either String ()" (fromCode foo1 :: Maybe (Either String ()))
    ng "(,)" (fromCode foo :: Maybe ((), ()))
    ng "(,,)" (fromCode foo :: Maybe ((), (), ()))
    ng "(,,,)" (fromCode foo :: Maybe ((), (), (), ()))
    ng "(,,,,)" (fromCode foo :: Maybe ((), (), (), (), ()))
    ng "(,,,,,)" (fromCode foo :: Maybe ((), (), (), (), (), ()))
    ng "Data.Functor.Sum.Sum Maybe (Either String) ()"
       (fromCode foo1 :: Maybe (Sum.Sum Maybe (Either String) ()))
    ng "Data.Functor.Product.Product Maybe (Either String) ()"
       (fromCode foo2 :: Maybe (Product.Product Maybe (Either String) ()))
    ng "Sum ()" (fromCode foo1 :: Maybe (Sum ()))
    ng "Proxy ()" (fromCode foo :: Maybe (Proxy ()))
    ng "Atom" (fromCode foo1 :: Maybe Atom)

fromCodeErrorTest :: Spec
fromCodeErrorTest = do
  let err :: Homoiconic a => a -> String -> Spec
      err a ty = it ("should show error message containing " ++ ty) $ do
        case parseCode x `asTypeOf` Homoiconic.Success a of
          Homoiconic.Failure e -> e `shouldSatisfy` isSubsequenceOf ty
          _           -> expectationFailure "parser unexpectedly succeeded"
      x = toCode (List [ toCode (aSymbol "foo")
                       , toCode (aSymbol "x") ])
  describe "Failing to parse Code" $ do
    err (1 :: Integer) "integral"
    err (1 :: Double) "fractional"
    err () "()"
    err 'a' "Char"
    err "string" "String"
    err False "Bool"
    err EQ "Ordering"
    err (Just False) "Maybe"
    err (Right True `asTypeOf` Left "foo") "Either"
    err ((),()) ","
    err ((),(),()) "(,,)"
    err ((),(),(),()) "(,,,)"
    err ((),(),(),(),()) "(,,,,)"
    err ((),(),(),(),(),()) "(,,,,,)"
    err (Sum.InL (Just True) `asTypeOf` Sum.InR (Just False)) "Sum"
    err (All True) "All"
    err (Proxy :: Proxy Int) "Proxy"
    err (aSymbol "foo") "Atom"

data D1 = D1a | D1b | D1c
  deriving (Eq, Show, Data, Typeable, Generic)

instance Homoiconic D1
instance Arbitrary D1 where
  arbitrary = oneof (map pure [D1a, D1b, D1c])

data D2 a = D2a a | D2b a a
  deriving (Eq, Show, Data, Typeable, Generic)

instance Homoiconic a => Homoiconic (D2 a)
instance Arbitrary a => Arbitrary (D2 a) where
  arbitrary = oneof [D2a <$> arbitrary, D2b <$> arbitrary <*> arbitrary]

data D3 a b = D3a Int a b
  deriving (Eq, Show, Data, Typeable, Generic)

instance (Homoiconic a, Homoiconic b) => Homoiconic (D3 a b)
instance (Arbitrary a, Arbitrary b) => Arbitrary (D3 a b) where
  arbitrary = D3a <$> arbitrary <*> arbitrary <*> arbitrary

data D4 a = D4a (a, a, a, a)
  deriving (Eq, Show, Data, Typeable, Generic)

instance Homoiconic a => Homoiconic (D4 a)
instance Arbitrary a => Arbitrary (D4 a) where
  arbitrary =
    D4a <$> ((,,,) <$> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary)

dataToCodeTest :: Spec
dataToCodeTest = do
  let s = toCode . aSymbol
  describe "D1 to Code" $
    it "should match D1 values" $ do
      let e1a = s "D1a"
          e1b = s "D1b"
          e1c = s "D1c"
      dataToCode D1a `shouldBe` e1a
      dataToCode D1b `shouldBe` e1b
      dataToCode D1c `shouldBe` e1c
  describe "D2 to Code" $
    it "should match `D2 Char' values"$ do
      let e2a = toCode (List [s "D2a", toCode 'x'])
          e2b = toCode (List [s "D2b", toCode 'y', toCode 'z'])
      dataToCode (D2a 'x') `shouldBe` e2a
      dataToCode (D2b 'y' 'z') `shouldBe` e2b
  describe "D2 with D1 to Code" $
    it "should match `D2 D1' values" $ do
      let e2b = toCode (List [s "D2b", s "D1a", s "D1b"])
      dataToCode (D2b D1a D1b) `shouldBe` e2b
  describe "D2 with Double to Code" $
    it "should match `D2 Double' values" $ do
      let e2c = toCode (List [s "D2a", toCode (1.23 :: Double)])
      dataToCode (D2a (1.23 :: Double)) `shouldBe` e2c
  describe "D3 to Code" $
    it "should match `D3' value" $ do
      let e3 = toCode (List [ s "D3a"
                            , toCode (42 :: Int)
                            , toCode False
                            , toCode 'a' ])
      dataToCode (D3a 42 False 'a') `shouldBe` e3
  describe "D4 to Code" $
    it "should match `D4' value" $ do
       let e4 = toCode (List [ s "D4a"
                             , toCode (List [ s ","
                                            , toCode 'w'
                                            , toCode 'x'
                                            , toCode 'y'
                                            , toCode 'z'])])
       dataToCode (D4a ('w', 'x', 'y', 'z')) `shouldBe` e4

genericHomoiconicTest :: Spec
genericHomoiconicTest = do
  describe "Result data type for parsing" $ do
    let s1, s2, f1 :: Homoiconic.Result Bool
        s1 = Homoiconic.Success True
        s2 = Homoiconic.Success False
        f1 = Homoiconic.Failure "foo"
    it "should equal to itself" $ do
      s1 `shouldBe` s1
      s1 `shouldNotBe` s2
      f1 `shouldNotBe` s2
    it "should show its conteints" $ do
      show s1 `shouldBe` "Success True"
      show f1 `shouldBe` "Failure \"foo\""
    let isFailure (Homoiconic.Failure _) = True
        isFailure _                      = False
    it "should fail on binding Failure" $
      ((Homoiconic.Failure "X_X" >>= pure) :: Homoiconic.Result Bool)
        `shouldSatisfy` isFailure

  describe "parseCode failure" $ do
    let d2a_sym = toCode $ aSymbol "D2a"
        d3a_sym = toCode $ aSymbol "D3a"
        unit_sym = toCode ()
        fail_with substring r = case r of
          Homoiconic.Failure msg -> substring `isSubsequenceOf` msg
          _                      -> False

    it "should fail on non list Code" $ do
      let r = parseCode @(D2 ()) d2a_sym
      r `shouldSatisfy` fail_with "Not a list"
    it "should fail on unexpected field constructor" $ do
      let r = parseCode @(D3 (D2 Bool) Bool) x
          x = toCode (List [d3a_sym, unit_sym])
      r `shouldSatisfy` fail_with "()"
    it "should fail on unexpected leftover values" $ do
      let r = parseCode @(D2 ()) x
          x = toCode (List [d2a_sym, unit_sym, unit_sym])
      r `shouldSatisfy` fail_with "leftover"
    it "should fail on invalid constructor (no arg)" $ do
      let r = parseCode @FormTest.D1 (toCode (aSymbol "foo"))
      r `shouldSatisfy` fail_with "foo"
    it "should fail on invalid constructor (with args)" $ do
      let r = parseCode @(D2 Bool) (toCode (List [toCode (aSymbol "foo")]))
      r `shouldSatisfy` fail_with "foo"

  describe "genericFromCode" $ do
    it "should result in Code value" $ do
      let d2 :: D2 Int
          d2 = D2a 42
      genericFromCode (toCode d2) `shouldBe` Just d2
    it "should fail with invalid value" $ do
      genericFromCode (toCode ()) `shouldBe` (Nothing :: Maybe (D2 Int))

  describe "To/from code via Generic" $ do
    it "should return itself" $
      let g_to_from :: D4 (D3 (D2 Int) D1) -> Bool
          g_to_from x = Just x == fromCode (toCode x)
      in  property g_to_from
    it "should return itself with nested Maybe data type" $
      let g_to_from :: D3 (Maybe (Maybe (Maybe Int))) (Maybe (Maybe Bool))
                    -> Bool
          g_to_from x = Just x == fromCode (toCode x)
      in  property g_to_from

  let s = toCode . aSymbol
      noParse :: (Eq a, Show a, Homoiconic a) => Maybe a -> IO ()
      noParse = (`shouldBe` Nothing)

  describe "Malformed codes via Generic" $ do
    it "result to Nothing with invalid constructor" $
      noParse (fromCode (s "Foo") :: Maybe D1)
    it "result to Nothing with arity mismatch" $
      noParse
        (fromCode (toCode (List [s "D2a", toCode False, toCode False]))
          :: Maybe (D2 Bool))
    it "result to Nothing with invalid field value" $
      noParse
        (fromCode (toCode (List [s "D2a", s "FOO"])) :: Maybe (D2 Bool))
    it "result to Nothing with non list" $
      noParse (fromCode (toCode "FOO") :: Maybe (D2 Bool))
    it "result to Nothing with invalid form" $
      noParse
        (fromCode (toCode (List [s "D2a", toCode (List [s "Just", s "___"])]))
           :: Maybe (D2 (Maybe Bool)))

unquoteSpliceTest :: Spec
unquoteSpliceTest =
  describe "unquote splicing List" $
    it "should return list contents" $
      property
        (\form ->
           (isListL form || isHsListL form ||
            isStringL form || isUnitL form)
           ==> (0 <= length (unquoteSplice form)))

eqForm :: Code -> Code -> Bool
eqForm a b =
  case (unCode a, unCode b) of
    -- Ignoring rounding error for fractional literals.
    (Atom (AFractional x), Atom (AFractional y))
      -> abs (fl_value x - fl_value y) <= toRational epsilon

    -- Recursively compare with `eqForm' for 'List' and 'HsList'.
    (List [], List []) -> True
    (List (x:xs), List (y:ys)) ->
      let ql z = qList z "" 0 0 0 0
      in  eqForm x y && eqForm (ql xs) (ql ys)

    (HsList [], HsList []) -> True
    (HsList (x:xs), HsList (y:ys)) ->
      let qh z = qHsList z "" 0 0 0 0
      in  eqForm x y && eqForm (qh xs) (qh ys)

    -- Treating empty 'List' and Atom symbol 'nil' as same value.
    (Atom (ASymbol sym), List []) | sym == fsLit "nil" -> True
    (List [], Atom (ASymbol sym)) | sym == fsLit "nil" -> True

    -- Using '==' for other Atom values.
    (Atom x, Atom y) -> x == y

    _ -> False
  where
    epsilon :: Double
    epsilon = 1e-7

pprTest :: Spec
pprTest =
  describe "ppr" $
    it "should return expected String" $ do
      let str0 = "(1 2.34 () #'a \"foo\" [foo bar buzz])"
      str1 <- runPpr str0
      str1 `shouldBe` str0

runPpr :: String -> IO String
runPpr str =
  runFnk (do dflags <- getDynFlags
             return (showPpr dflags (parseE str)))
         defaultFnkEnv

parseE :: String -> Code
parseE = parseE' Nothing

parseE' :: Maybe FilePath -> String -> Code
parseE' mb_path str =
  let inp = stringToStringBuffer str
  in  case runSP sexpr mb_path inp of
        Right (expr, _) -> expr
        Left err        -> error (show err)

isListL :: Code -> Bool
isListL (LForm (L _ (List _))) = True
isListL _                      = False

isHsListL :: Code -> Bool
isHsListL (LForm (L _ (HsList _))) = True
isHsListL _                        = False

isStringL :: Code -> Bool
isStringL (LForm (L _ (Atom (AString _ _)))) = True
isStringL _                                  = False

isUnitL :: Code -> Bool
isUnitL (LForm (L _ (Atom AUnit))) = True
isUnitL _                          = False

max_stars_repo_name: finkel-lang/finkel

max_issues_repo_name: finkel-lang/finkel

max_forks_repo_name: finkel-lang/finkel


text:
module Day01 (solveDay01) where

import           Data.Char       (isDigit)
import           Data.Complex
import           Data.Foldable   (foldl')
import           Data.List.Split (splitOn)
import           Data.Map.Strict (Map)
import qualified Data.Map.Strict as M

type Point   = Complex Double
type Instr   = String
type Command = (Instr, Point)

origin :: Point
origin = 0 :+ 0

commands :: [String] -> [Command]
commands = map parseCommand
  where
    parseCommand s = (dir, dist)
      where dir  = take 1 s
            dist = (read (takeWhile isDigit (drop 1 s)) :: Double) :+ 0

-- √-1
i :: Complex Double
i = 0 :+ 1

step :: (Point, Point) -> Command -> (Point, Point)
step (heading, point) ("R", dist) = (h', (h' * dist) + point)
  where h' = heading * conjugate i
step (heading, point) ("L", dist) = (h', (h' * dist) + point)
  where h' = heading * i
step _ _                          = error "Unrecognized command"

step' :: (Point, [Point]) -> Command -> (Point, [Point])
step' (heading, points) (cmd, dist)
  | cmd == "R" = let h' = heading * conjugate i
                 in (h', points ++ explodePoint (h' * dist))
  | cmd == "L" = let h' = heading * i
                 in (h', points ++ explodePoint (h' * dist))
  | otherwise  = error "Unrecognized command in `steps`"

processCommands' :: [Command] -> [Point]
processCommands' = snd . foldl' step' (north, [])
  where north = 0 :+ 1

processCommands :: [Command] -> Point
processCommands = snd . foldl' step (north, origin)
  where north = 0 :+ 1

calcAns :: Point -> Double
calcAns (x :+ y) = abs x + abs y

explodePoint :: Point -> [Point]
explodePoint (x :+ y)
  | x == 0  = replicate (abs . floor $ y) (0 :+ signum y)
  | y == 0  = replicate (abs . floor $ x) (signum x :+ 0)
  | otherwise = error "Invalid point provided to `explodePoint`"

walk :: [Point] -> Point
walk ps = go ps origin M.empty
  where
    go :: [Point] -> Point -> Map (Double, Double) Bool -> Point
    go []     _ _ = error "Did not encounter any places twice"
    go (s:ss) p m = case M.lookup (x, y) m of
                        Just _  -> p'
                        Nothing -> go ss p' m'
      where p'@(x :+ y) = p + s
            m'          = M.insert (x, y) True m

solveDay01 :: FilePath -> IO ()
solveDay01 pathToInput = do
  putStrLn "Solution for day one: "
  ls <- readFile pathToInput
  let cleanInput = splitOn ", " ls
  -- Part 1
  print . calcAns . processCommands . commands $ cleanInput
  -- Part 2
  print . calcAns . walk . processCommands' . commands $ cleanInput
  putStrLn ""

ext: hs

max_stars_repo_name: justanotherdot/advent-linguist

max_issues_repo_name: justanotherdot/advent-linguist

max_forks_repo_name: justanotherdot/advent-linguist


text:
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE MultiWayIf   #-}
{-# LANGUAGE RankNTypes   #-}
{-# LANGUAGE ViewPatterns #-}
module Main where

import Control.Lens                       hiding (inside, transform)
import Control.Monad.Extra
import Control.Monad.ST
import Data.Array.Repa
import Data.Complex
import Data.Matrix                        (Matrix)
import Data.STRef
import GHC.Base                           (until)
import Graphics.Gloss.Interface.Pure.Game hiding (Vector, shift)
import Graphics.Gloss.Raster.Array
import Relude

import qualified Data.ByteString as B
import qualified Data.Matrix     as Mtx

data World = World { _v   :: Matrix Float
                   , _mtx :: Matrix Float
                   } deriving (Eq,Show)
makeLenses ''World

initialWorld :: World
initialWorld = World { _v = Mtx.fromLists [[-2],[-1]] -- FIXME no hardcode
                     , _mtx = Mtx.fromLists [[1/250,0],[0,1/250]]
                     }

type Vector a = Matrix a

data Settings = Settings { _windowSize :: (Int, Int)
                         , _offset     :: (Int, Int)
                         , _pointSize  :: (Int,Int)
                         , _threshold  :: forall a. Num a => a
                         , _fps        :: Int
                         , _iteration  :: Int
                         }
makeLenses ''Settings

settings :: Settings
settings = Settings { _windowSize = (750,500)
                    , _offset     = (250,250)
                    , _pointSize  = (1,1)
                    , _threshold  = 2
                    , _fps        = 60
                    , _iteration  = 100
                    }

window :: Display
window = InWindow "fractol" (settings^.windowSize) (settings^.offset)

toTup :: DIM2 -> (Int, Int)
toTup (Z :. y :. x) = (x,y)

toVec :: (Int, Int) -> Vector Float
toVec (!x, !y) = seq y $ seq x $ on f fromIntegral x y
  where
    f !x1 !x2 = Mtx.fromLists [[x1],[x2]]
{-# INLINE toVec #-}

transform :: World -> Vector Float -> Vector Float
transform !w !p = w^.mtx * p + w^.v
{-# INLINE transform #-}

mandel :: World -> Array D DIM2 Color
mandel !w = fromFunction (ix2 500 750) render
  where
    render (iter . transform w . toVec . toTup -> (len, n)) =
      let diff = len - settings^.threshold
          m = 2 * n
      in if diff >= 0
           then rgb' (0.9 * diff) (0.6 * diff) (0.1 * diff)
           else black

toComplex :: Vector Float -> Complex Float
toComplex !v = let [!x,!y] = Mtx.toList v
          in x :+ y
{-# INLINE toComplex #-}

iter :: Vector Float -> (Float,Int)
iter (toComplex -> !c) = let (!zn, !len, !n) = until condition endo start
                        in  (len, n)
  where
    start = (0,0,0)
    endo (!zn, !len, !n) = let new = zn * zn + c
                          in (new, magnitude new, succ n)
    {-# INLINE endo #-}
    condition (_, !len, !n) = len >= settings^.threshold
                                || n > settings^.iteration
    {-# INLINE condition #-}

{-# INLINE iter #-}

cast :: (Float, Float) -> (Int, Int)
cast (!x,!y) = on (,) floor x y
{-# INLINE cast #-}

shift :: Vector Float -> Vector Float
shift = (+ Mtx.fromLists [[375],[250]])
{-# INLINE shift #-}

eventHandler :: Event -> World -> World
eventHandler (EventKey (Char 'p') Up _ (shift . toVec . cast -> p)) w =
  let truePoint = transform w p
  in w & v -~ Mtx.scaleMatrix 0.1 truePoint
       & mtx %~ (Mtx.scaleMatrix 1.1 (Mtx.identity 2) *)
eventHandler (EventKey (Char 'n') Up _ (shift . toVec . cast -> p)) w =
  let truePoint = transform w p
  in w & v +~ Mtx.scaleMatrix 0.1 truePoint
       & mtx %~ (Mtx.scaleMatrix 0.9 (Mtx.identity 2) *)
eventHandler _ w = w

main :: IO ()
main = do
  let fractol = playArray window (settings^.pointSize) (settings^.fps)
  fractol initialWorld mandel eventHandler (const id)

max_stars_repo_name: krendelhoff/fractol

max_issues_repo_name: krendelhoff/fractol

max_forks_repo_name: krendelhoff/fractol


text:
-- |
-- Module    : GRN.Density
-- Copyright : (c) 2011 Jason Knight
-- License   : BSD3
--
-- Maintainer  : jason@jasonknight.us
-- Stability   : experimental
-- Portability : portable
--
-- Simulates many different random values in the unknown Kmap entries and
-- estimates the density of this distribution for each gene.
--
-- Then we will convolve this with the 'measurement kernel' to provide a non
-- dirac prior for every gene (even deterministic ones).
--

module GRN.Density where

import GRN.StateTransition
import GRN.Render
import GRN.Sparse
import GRN.Types

import qualified Data.Map as M
import Data.List
import Control.Monad
import Text.Printf
import System.Console.ParseArgs

import qualified Data.Vector as V
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Generic as G

import Control.Parallel.Strategies

{-import Numeric.FFT.Vector.Unitary-}

import Data.Graph.Analysis.Algorithms.Common (componentsOf)

import Statistics.Sample.KernelDensity
import Statistics.Sample
import Numeric.MathFunctions.Constants
import Statistics.Distribution hiding (stdDev)
import Statistics.Distribution.Normal
import Graphics.Gnuplot.Simple hiding (Points)

step = 0.01 -- Step size
nl = 101 :: Int -- Number of poitns
smoothingVar = 0.001 -- Variance of smoothing gaussian

points = U.enumFromStepN 0.00 step nl

normSmooth mu = G.map (density $ normalDistr mu smoothingVar) points
stdNorm = normSmooth 0.5

myKernelEst :: U.Vector Double -> U.Vector Double
myKernelEst sample = fst $ kde_ nl 0.0 1.0 sample

myEst :: U.Vector Double -> U.Vector Double
myEst xv  | stdDev xv < 0.001 = normalize.normSmooth $ (G.sum xv / (fromIntegral.G.length $ xv))
          | otherwise = normalize.clamp.convolve stdNorm.myKernelEst $ xv

normalize xv = let scale = 1.0 / (step * G.sum xv) in G.map (*scale) xv
clamp xv = G.map (\x->if x<0 then 0 else x) xv

-- Assume that both xs and ys are the same size
convolve xs ys = error "FFTW out of date." --G.slice start l uncentered
  {-where-}
    {-xsF = run dftR2C (xs G.++ zeros)-}
    {-ysF = run dftR2C (ys G.++ zeros)-}
    {-l = G.length xs-}
    {-start = floor $ (fromIntegral l)/2-}
    {-zeros = U.replicate l 0.0-}
    {-mult = G.zipWith (*) xsF ysF-}
    {-uncentered = run dftC2R mult-}

calcEsts :: GeneMC -> AttractorDensity
calcEsts = M.map myEst

simRuns :: Args String -> ParseData -> IO () --U.Vector SSA
simRuns args p = do
  let
    titles = map (\x->defaultStyle {lineSpec = CustomStyle [LineTitle x]} ) (M.keys p)
    titleAll = zip titles
    prepForPlot = titleAll.map (G.toList.G.zip points.snd).M.toList
    ks = buildKmaps p
    dm = getRequiredArg args "submode" -- Number of simruns

    vals = concatMap (\(Kmap _ _ x)->M.elems x) (M.elems ks)
    unknowns = [C x | C x <- vals] ++ [X | X <- vals]

  print unknowns

  when (dm == "m") $ do
    let
      allssas = calcSSAsDOK args p ks
      allests = calcEsts allssas
    --print allssas
    plotPathsStyle [XRange (0.0,1.0)] (prepForPlot allests)

  when (dm == "g") $ do
    let
      grComps = componentsOf $ simulate args $ kmapToStateGraph ks
      grNComps = map normalizeGraph grComps
      allests = map (\(d,g)->(d, g, calcEsts.calcSSAs args p ks $ g)) grNComps

    mapM_ (plotPathsStyle [XRange (0.0,1.0)] . prepForPlot . (\(_,_,x)->x)) allests
    --when gen $ drawStateGraph snd.G.head $ simVec) args

max_stars_repo_name: binarybana/grn-pathways

max_issues_repo_name: binarybana/grn-pathways

max_forks_repo_name: binarybana/grn-pathways


text:
-- {-# DeriveDataTypeable #-}
module Main where

  import System.Environment
  import System.Random
  import Numeric.LinearAlgebra (Matrix, assoc, Z, toLists)
  import Control.Monad (join)
  import Data.Char (toLower)

  import System.Console.CmdArgs.Implicit

  import ODMatrix.Sim.BusState

  data SimParams = Params {
      size :: Int,
      capacity :: Int,
      seed :: Int,
      cases :: Int,
      lambdas :: Bool
  } deriving (Show, Data, Typeable)

  params = Params {
    size = 10 &= help "Size of the route",
    capacity = 50 &= help "Capacity of the unit",
    seed = 0 &= help "Seed for the simulation",
    cases = 10 &= help "Number of cases to simulate",
    lambdas = False &= help "Indicates if the lists of lambdas will be provided manually"
  }

  getLambdas :: RandomGen g => Bool -> g -> Int -> Int -> IO ([Double],[Double])
  getLambdas False g n k = do -- Random genetated
    let (bg,ag) = split g
        bs = take n $ randomRs (0,fromIntegral k) bg
        as = take n $ randomRs (0,fromIntegral k) ag
    return (bs, as)
  getLambdas True _ n _ = do  -- Manually provided
    input <- getContents
    let ls = lines input
        bs = take n ls
        as = take n $ drop n ls
    return (map read bs, map read as)

  main :: IO ()
  main = do
    -- ODMatrix Simulation
    putStrLn "{"

    args <- cmdArgs params

    let (Params n k seed cases manual_lambdas) = args
        (g1, g2) = split $ mkStdGen seed

    putStrLn "\"params\": {"
    putStrLn $ "\"size\":" ++ (show n) ++ ","
    putStrLn $ "\"capacity\":" ++ (show k) ++ ","
    putStrLn $ "\"seed\":" ++ (show seed) ++ ","
    putStrLn $ "\"cases\":" ++ (show cases) ++ ","
    putStrLn $ "\"manual_lambdas\":" ++ (map toLower $ show manual_lambdas)
    putStrLn "},"

    -- putStrLn $  if manual_lambdas
    --             then "# Reading lambdas..."
    --             else "# Generating randoms lambdas..."

    (bs,as) <- getLambdas manual_lambdas g1 n k

    putStrLn "\"lambdas\" : {"

    -- Boarding lambdas
    putStr "\"boards\": "
    putStr . show $ bs
    putStrLn ","

    -- Alighting lambdas
    putStr "\"alights\": "
    putStrLn . show $ as

    putStrLn "},"

    -- Generating origin destination matrices

    putStrLn "\"odms\": {"
    let ms = map (assoc (n,n) 0) (busSims k bs as g2) :: [Matrix Z]
        ms' = map (\(i,m) -> "\"" ++ (show i) ++ "\": " ++ (show m) ++ ",\n") .
              zip [0..] .
              take cases . map (join . toLists) $ ms

    putStrLn . init . init . join $ ms'
    putStrLn "}"

    putStrLn "}"

max_stars_repo_name: renecura/odmatrix-sim

max_issues_repo_name: renecura/odmatrix-sim

max_forks_repo_name: renecura/odmatrix-sim


text:
module STCR2Z2T0S0PointSet where

import           Data.Array.Repa         as R
import           Data.Binary             (decodeFile)
import           Data.Complex
import           Data.List               as L
import           DFT.Plan
import           FokkerPlanck.MonteCarlo
import           FokkerPlanck.Pinwheel
import           Image.IO
import           STC.PowerMethod
import           System.Directory
import           System.Environment
import           System.FilePath
import           Types
import Text.Printf

main = do
  args@(numPointStr:numOrientationStr:numScaleStr:thetaSigmaStr:scaleSigmaStr:maxScaleStr:taoStr:numTrailStr:maxTrailStr:theta0FreqsStr:thetaFreqsStr:scale0FreqsStr:scaleFreqsStr:histFilePath:numIterationStr:writeSourceFlagStr:rStr:numThreadStr:_) <-
    getArgs
  print args
  let numPoint = read numPointStr :: Int
      numOrientation = read numOrientationStr :: Int
      numScale = read numScaleStr :: Int
      thetaSigma = read thetaSigmaStr :: Double
      scaleSigma = read scaleSigmaStr :: Double
      maxScale = read maxScaleStr :: Double
      tao = read taoStr :: Double
      numTrail = read numTrailStr :: Int
      maxTrail = read maxTrailStr :: Int
      theta0Freq = read theta0FreqsStr :: Double
      theta0Freqs = [-theta0Freq .. theta0Freq]
      thetaFreq = read thetaFreqsStr :: Double
      thetaFreqs = [-thetaFreq .. thetaFreq]
      scale0Freq = read scale0FreqsStr :: Double
      scaleFreq = read scaleFreqsStr :: Double
      scale0Freqs = [-scale0Freq .. scale0Freq]
      scaleFreqs = [-scaleFreq .. scaleFreq]
      numIteration = read numIterationStr :: Int
      writeSourceFlag = read writeSourceFlagStr :: Bool
      numThread = read numThreadStr :: Int
      folderPath = "output/test/STCR2Z2T0S0PointSet"
  createDirectoryIfMissing True folderPath
  flag <- doesFileExist histFilePath
  radialArr <-
    if flag
      then R.map magnitude . getNormalizedHistogramArr <$>
           decodeFile histFilePath
      else do
        putStrLn "Couldn't find a Green's function data. Start simulation..."
        solveMonteCarloR2Z2T0S0Radial
          numThread
          numTrail
          maxTrail
          numPoint
          numPoint
          thetaSigma
          scaleSigma
          maxScale
          tao
          theta0Freqs
          thetaFreqs
          scale0Freqs
          scaleFreqs
          histFilePath
          (emptyHistogram
             [ (round . sqrt . fromIntegral $ 2 * (div numPoint 2) ^ 2)
             , L.length scale0Freqs
             , L.length theta0Freqs
             , L.length scaleFreqs
             , L.length thetaFreqs
             ]
             0)
  arrR2Z2T0S0 <-
    computeUnboxedP $
    computeR2Z2T0S0ArrayRadial
      (cutoff 8 radialArr)
      numPoint
      numPoint
      1
      maxScale
      thetaFreqs
      scaleFreqs
      theta0Freqs
      scale0Freqs
  plan <- makeR2Z2T0S0Plan emptyPlan arrR2Z2T0S0
  let n = 30
      m = round $ (fromIntegral n) * (sqrt 2) / 2
      a = 20
      b = 7
      c = 50
      a' = round $ (fromIntegral a) * (sqrt 2) / 2
      b' = round $ (fromIntegral b) * (sqrt 2) / 2
      c' = round $ (fromIntegral c) * (sqrt 2) / 2
      r = read rStr :: Double
      numTheta = 8
      deltaTheta = (2 * pi) / numTheta
      xs = -- [R2S1RPPoint (i,0,0,1) | i <- [-30,-20..30]]
        -- ((L.map
        --     (\(i, j) -> R2S1RPPoint (round i, round j, 0, 1))
        --     [ (r * cos (k * deltaTheta) + 0, r * sin (k * deltaTheta) + 0)
        --     | k <- [0 .. numTheta-1]
        --     ]) -- L.++
         -- (L.map
         --    (\(i, j) -> R2S1RPPoint (round i, round j, 0, 1))
         --    [ ((r + 2) * cos (k * deltaTheta), (r + 2) * sin (k * deltaTheta))
         --    | k <- [0 .. numTheta - 1]
         --    ]) L.++
         -- (L.map
         --    (\(i, j) -> R2S1RPPoint (round i, round j, 0, 1))
         --    [ ( (r + 5) * cos (k * deltaTheta)
         --      , (r + 5) * sin (k * deltaTheta))
         --    | k <- [0 .. numTheta - 1]
         --    ])
               -- L.++ [R2S1RPPoint (3, 5, 0, 1)]
         -- )
               ([R2S1RPPoint (i, i, 0, 0) | i <- [a',a' + b' .. c']] L.++
                [R2S1RPPoint (i, -i, 0, 0) | i <- [-a',-(a' + b') .. -c']] L.++
                [R2S1RPPoint (i, i, 0, 0) | i <- [-a',-(a' + b') .. -c']] L.++
                [R2S1RPPoint (i, -i, 0, 0) | i <- [a',a' + b' .. c']] L.++
                [R2S1RPPoint (i, 0, 0, 0) | i <- [a,a + b .. c]] L.++
                [R2S1RPPoint (0, i, 0, 0) | i <- [-a,-(a + b) .. -c]] L.++
                [R2S1RPPoint (i, 0, 0, 0) | i <- [-a,-(a + b) .. -c]] L.++
                [R2S1RPPoint (0, i, 0, 0) | i <- [a,a + b .. c]] -- L.++
                -- [R2S1RPPoint (3, 5, 0, 1)]
               )
                -- [ r2s1rppoint (n, 0, 0, 1)
                -- , R2S1RPPoint (0, n, 0, 1)
                -- , R2S1RPPoint (-n, 0, 0, 1)
                -- , R2S1RPPoint (0, -n, 0, 1)
                -- , R2S1RPPoint (m, m, 0, 1)
                -- , R2S1RPPoint (-m, m, 0, 1)
                -- , R2S1RPPoint (m, -m, 0, 1)
                -- , R2S1RPPoint (-m, -m, 0, 1)
                -- ]
  let bias = computeBiasR2T0S0 numPoint numPoint theta0Freqs scale0Freqs xs
      eigenVec =
        computeInitialEigenVectorR2T0S0
          numPoint
          numPoint
          theta0Freqs
          scale0Freqs
          thetaFreqs
          scaleFreqs
          xs
  powerMethodR2Z2T0S0Reversal
    plan
    folderPath
    numPoint
    numPoint
    numOrientation
    thetaFreqs
    theta0Freqs
    numScale
    scaleFreqs
    scale0Freqs
    arrR2Z2T0S0
    numIteration
    writeSourceFlag
    -- (printf "_%d" (round r :: Int))
    (printf
       "_%d_%d_%d_%d_%d_%.2f_%.2f"
       numPoint
       (round thetaFreq :: Int)
       (round scaleFreq :: Int)
       (round maxScale :: Int)
       (round tao :: Int)
       thetaSigma
       scaleSigma)
    0.5
    0.001
    bias
    eigenVec

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
module Prelude
  ( module P
  ) where

import Data.Tuple as P (fst,snd)
import Data.String as P (String)
import GHC.Show as P (Show(..))
import Data.Int as P (Int)
import Data.Bool as P (Bool(..),not)
import Data.Monoid as P (Monoid(..))
import Data.Semigroup as P (Semigroup(..))
import Data.Complex as P (Complex(..))
import Data.Functor as P (Functor(..))
import Data.Semiring as P (Semiring(..), Ring(..), (+), (*), (-))
import Control.Monad as P (Monad(..), (=<<))
import Control.Applicative as P (Applicative(..))
import Data.Eq as P (Eq(..))
import Data.Ord as P (Ord(..))
import GHC.Base as P (Double)
import Data.Function as P (($),(.),(&),id,flip)
import Data.Map as P (Map)
import GHC.Float as P (Floating(..))
import GHC.Real as P (fromIntegral, (/))

max_stars_repo_name: chessai/qprob

max_issues_repo_name: chessai/qprob

max_forks_repo_name: chessai/qprob


text:
{-# LANGUAGE DataKinds                            #-}
{-# LANGUAGE FlexibleContexts                     #-}
{-# LANGUAGE FlexibleInstances                    #-}
{-# LANGUAGE MultiParamTypeClasses                #-}
{-# LANGUAGE PartialTypeSignatures                #-}
{-# LANGUAGE TupleSections                        #-}
{-# LANGUAGE TypeApplications                     #-}
{-# LANGUAGE TypeOperators                        #-}
{-# LANGUAGE TypeSynonymInstances                 #-}
{-# LANGUAGE UndecidableInstances                 #-}
{-# OPTIONS_GHC -fno-warn-orphans                 #-}
{-# OPTIONS_GHC -fno-warn-partial-type-signatures #-}

import           Backprop.Learn
import           Control.Monad.Trans.Maybe
import           Data.Bitraversable
import           Data.Default
import           Data.IDX
import           Data.Traversable
import           Data.Tuple
import           Numeric.LinearAlgebra.Static.Backprop
import           Numeric.Opto
import           Numeric.Opto.Run.Simple
import           System.Environment
import           System.FilePath
import           Text.Printf
import qualified Data.Vector.Generic                   as VG
import qualified Numeric.LinearAlgebra                 as HM
import qualified Numeric.LinearAlgebra.Static          as H
import qualified System.Random.MWC                     as MWC

mnistNet :: Model _ _ (R 784) (R 10)
mnistNet = fca     @300 softMax
        <~ dropout      0.25
        <~ fca     @784 logistic

main :: IO ()
main = MWC.withSystemRandom $ \g -> do
    datadir:_  <- getArgs
    Just train <- loadMNIST (datadir </> "train-images-idx3-ubyte")
                            (datadir </> "train-labels-idx1-ubyte")
    Just test  <- loadMNIST (datadir </> "t10k-images-idx3-ubyte")
                            (datadir </> "t10k-labels-idx1-ubyte")
    putStrLn "Loaded data."
    net0 <- initParamNormal mnistNet 0.2 g

    let so = def { soTestSet   = Just test
                 , soEvaluate  = runTest
                 , soSkipSamps = 2500
                 }
        runTest chnk net = printf "Error: %.2f%%" ((1 - score) * 100)
          where
            score = testModelAll maxIxTest mnistNet (TJust net) chnk

    simpleRunner so train SOSingle def net0
      (adam def $ modelGradStoch crossEntropy noReg mnistNet g)
      g

loadMNIST
    :: FilePath
    -> FilePath
    -> IO (Maybe [(R 784, R 10)])
loadMNIST fpI fpL = runMaybeT $ do
    i <- MaybeT          $ decodeIDXFile       fpI
    l <- MaybeT          $ decodeIDXLabelsFile fpL
    d <- MaybeT . return $ labeledIntData l i
    MaybeT . return $ for d (bitraverse mkImage mkLabel . swap)
  where
    mkImage = H.create . VG.convert . VG.map (\i -> fromIntegral i / 255)
    mkLabel n = H.create $ HM.build 10 (\i -> if round i == n then 1 else 0)

max_stars_repo_name: mstksg/backprop-learn

max_issues_repo_name: mstksg/backprop-learn

max_forks_repo_name: mstksg/backprop-learn


text:
-- |
-- Copyright  : (c) 2012 Daniël de Kok
-- License    : BSD3
--
-- Maintainer : Daniël de Kok <me@danieldk.eu>
-- Stability  : experimental
--
-- This module provides functionality to perform approximate randomization
-- tests (Noreen, 1989).

{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE DoAndIfThenElse #-}
{-# LANGUAGE FlexibleContexts #-}

module Statistics.Test.ApproxRand (
  -- * Description
  -- $description

  -- * Examples
  -- $examples

  -- * Data types
  TestOptions(..),
  TestResult(..),
  Significance(..),
  RandWithError,

  -- * Approximate randomization tests
  approxRandTest,
  approxRandStats,

  approxRandPairTest,
  approxRandPairStats,

  -- * Test statistics
  TestStatistic,
  differenceMean,
  meanDifference,
  varianceRatio
) where

import           Prelude hiding ((++))
import           Control.Monad (liftM, replicateM, when)
import           Control.Monad.Error (ErrorT)
import           Control.Monad.Error.Class (throwError)
import           Control.Monad.Mersenne.Random (R(..), Rand(..), getBool)
import           Control.Monad.ST (runST)
import           Control.Monad.Trans.Class (lift)
import           Data.Vector.Generic ((++))
import qualified Data.Vector.Generic as VG
import qualified Data.Vector.Generic.Mutable as GM
import           Data.Word (Word)
import           Statistics.Sample (variance)
import           Statistics.Test.Types (TestType(..))
import           Statistics.Types
import           System.Random.Mersenne.Pure64 (PureMT, randomInt, randomWord)

-- $description
--
-- Approximate randomization tests rely on a simple premise: given a test
-- statistic, if the null-hypothesis (the samples do not differ) is true,
-- we can randomly swap values between samples without an (extreme) impact
-- on the test statistic. Otherwise, the null-hypothesis must be rejected.
--
-- The test works by generating a given number of sample shuffles and computing
-- the test statistic for each shuffle. If /r/ is the number of shuffled
-- samples where the test statistic is at least as high as the test statistic
-- applied on the original samples; and /N/ the number of shuffles, then
-- the null-hypothesis is rejected iff /(r + 1):(N + 1) < p-value/ (for
-- one-sided tests).
--
-- Two kinds of test are supported:
--
-- * /Paired sample/ ('approxRandPairTest'): values from samples are shuffled
--   pair-wise. This requires the samples to have an equal length.
--
-- * /Unpaired sample/ ('approxRandTest'): values from samples are shuffled
--   among both samples. Consequently the i-th element of one sample does not
--   bear a relationship with the i-th element of the other sample. The
--   shuffled samples retain the sizes of the original samples.
--
-- Both tests can be performed as a one-tailed or two-tailed test.

-- $examples
-- Both unpaired and paired sample tests use the 'Rand' monad to obtain
-- random numbers. We can obtain a pseudo-random number generator that
-- is seeded using the system clock using the
-- 'System.Random.Mersenne.Pure64.newPureMT' function (please refer to
-- the documentation of 'System.Random.Mersenne.Pure64' for more
-- information):
--
-- > prng <- newPureMT
--
-- Suppose that we have the samples 's1' and 's2'. We could now perform
-- a Two-Tailed randomization test with 10,000 shuffles and the mean
-- difference as the test statistic, by running 'approxRandTest' in the 'Rand'
-- monad (at the /p = 0.01/ level):
--
-- > evalRandom (approxRandTest TwoTailed meanDifference 10000 0.01 s1 s2) prng
--
-- It is also possible to obtain the test statistics of the shuffled samples
-- directly (e.g. to inspect the distribution of test statistics) using the
-- 'approxRandStats'/'approxRandPiarStats' functions:
--
-- > evalRandom (approxRandStats meanDifference 10000 0.01 s1 s2) prng

-- | Computations with random numbers that can fail.
type RandWithError a = ErrorT String Rand a

-- | Options for randomization tests
--
data TestOptions = TestOptions {
    toTestType      :: TestType,      -- ^ Type of test ('OneTailed' or 'TwoTailed')
    toTestStatistic :: TestStatistic, -- ^ Test statistic
    toIterations    :: Int,           -- ^ Number of shuffled samples to create
    toPValue        :: Double         -- ^ he p-value at which to test (e.g. 0.05)
}

-- |
-- The result of hypothesis testing.
data TestResult = TestResult {
    trSignificance     :: Significance, -- ^ Significance
    trStat             :: Double,       -- ^ Test statistic for the samples
    trRandomizedStats  :: Sample        -- ^ Test statistics for the
                                        --   randomized samples
  } deriving (Eq, Ord, Show)

-- |
-- Significance.
data Significance =
    Significant    Double -- ^ The null hypothesis should be rejected
  | NotSignificant Double -- ^ Data is compatible with the null hypothesis
  deriving (Eq, Ord, Show)

-- |
-- Apply a pair-wise approximate randomization test.
--
-- In pair-wise approximate randomization tests the data points at a given
-- index are swapped between samples with a probability of 0.5. Since
-- swapping is pairwise, the samples should have the same length.
approxRandPairTest ::
     TestOptions              -- ^ Options for the test
  -> Sample                   -- ^ First sample
  -> Sample                   -- ^ Second sample
  -> RandWithError TestResult -- ^ The test result
approxRandPairTest (TestOptions testType stat n pTest) s1 s2 = do
  stats <- approxRandPairStats stat n s1 s2
  let tOrig = stat s1 s2
  let sig = significance testType pTest n $ countExtremes tOrig stats
  return $ TestResult sig tOrig stats

-- |
-- Apply an approximate randomization test.
--
-- In approximate randomization tests, the values of two samples are
-- shuffled among those samples. A test statistic is calculated for
-- the original samples and the shuffled samples, to detect whether the
-- difference of the samples is extreme or not.
approxRandTest ::
     TestOptions     -- ^ Options for the test
  -> Sample          -- ^ First sample
  -> Sample          -- ^ Second sample
  -> Rand TestResult -- ^ The test result
approxRandTest (TestOptions testType stat n pTest) s1 s2 = do
  stats <- approxRandStats stat n s1 s2
  let tOrig = stat s1 s2
  let sig = significance testType pTest n $ countExtremes tOrig stats
  return $ TestResult sig tOrig stats

-- | Determine the significance.
significance ::
     TestType     -- ^ Type of test ('OneTailed' or 'TwoTailed')
  -> Double       -- ^ The p-value at which to test (e.g. 0.05)
  -> Int          -- ^ Number of sample shuffles
  -> (Int, Int)   -- ^ Extreme statistic counts
  -> Significance -- ^ The test result
significance TwoTailed pTest n =
  significant (pTest / 2) . pValue n . uncurry min
significance OneTailed pTest n =
  significant pTest . pValue n . snd

-- | Wrap a p-value in a 'TestResult'.
significant ::
     Double       -- ^ The p-value at which to test
  -> Double       -- ^ The p-value
  -> Significance -- ^ The test result
significant pTail p =
  if p < pTail then Significant p else NotSignificant p

-- | Calculate a p-value
pValue ::
     Int    -- ^ Number of extreme outcomes
  -> Int    -- ^ Number of shuffles
  -> Double -- ^ The p-value
pValue n r = (fromIntegral r + 1) / (fromIntegral n + 1)

-- |
-- Count extreme test statistic values. If the test statistic value of the
-- original sample is in the right tail, we want to count values equal to
-- or larger than that value. If the value is in the left tail, we want to
-- count value smaller than or equal to that value. Since we do not know
-- the tail (yet), we count both.
--
-- Note: we can determine the tail by (1) averaging the test statistics of
-- the randomized samples, or (2) taking the smaller of the two counts.
countExtremes ::
     Double     -- ^ Test statistic value of the original samples
  -> Sample     -- ^ Test statistic values of the randomized samples.
  -> (Int, Int) -- ^ Count of left- and right-tail extremes.
countExtremes tOrig =
  VG.foldl' count (0, 0)
  where
    count (left, right) tPerm =
      let !newLeft = if tPerm <= tOrig then succ left else left in
      let !newRight = if tPerm >= tOrig then succ right else right in
        (newLeft, newRight)

-- |
-- Generate a given number of pairwise shuffled samples, and calculate
-- the test statistic for each shuffle.
--
-- Since the data points at a given index are swapped (with a probability of
-- 0.5), the samples should have the same length.
approxRandPairStats ::
     TestStatistic          -- ^ Test statistic
  -> Int                    -- ^ Number of shuffled samples to create
  -> Sample                 -- ^ First sample
  -> Sample                 -- ^ Second sample
  -> RandWithError Sample   -- ^ The statistics of the shuffles
approxRandPairStats stat n s1 s2 = do
  when (VG.length s1 /= VG.length s2) $
    throwError "Cannot calculate pairwise statistic: samples have different sizes"
  lift $ liftM VG.fromList $ replicateM n $
    uncurry stat `liftM` shuffleVectorsPairwise s1 s2

-- |
-- Generate a given number of shuffled samples, and calculate the test
-- statistic for each shuffle.
--
-- This function does not require the samples to have an equal length.
approxRandStats ::
     TestStatistic -- ^ Test statistic
  -> Int           -- ^ Number of shuffled samples to create
  -> Sample        -- ^ First sample
  -> Sample        -- ^ Second sample
  -> Rand Sample   -- ^ The statistics of the shuffles
approxRandStats stat n s1 s2 =
  liftM VG.fromList $ replicateM n $ uncurry stat `liftM` shuffleVectors s1 s2

-- | Pair-wise shuffle of two vectors.
shuffleVectorsPairwise :: (VG.Vector v a, VG.Vector v Bool) =>
  v a -> v a -> Rand (v a, v a)
shuffleVectorsPairwise vec1 vec2 = do
  randomVec <- randomVector (VG.length vec1)
  let pv1 = VG.zipWith3 permute vec1 vec2 randomVec
  let pv2 = VG.zipWith3 permute vec2 vec1 randomVec
  return (pv1, pv2)
  where
    permute val1 val2 coin =
      if coin then val1 else val2

randomVector :: VG.Vector v Bool => Int -> Rand (v Bool)
randomVector len =
  VG.replicateM len getBool

-- Shuffle values amongst two vectors, keeping the original vector lengths.
shuffleVectors :: VG.Vector v a => v a -> v a -> Rand (v a, v a)
shuffleVectors v1 v2 = do
  shuffledVectors <- shuffleVector $ v1 ++ v2
  return (VG.slice 0 (VG.length v1) shuffledVectors,
    VG.slice (VG.length v1) (VG.length v2) shuffledVectors)

-- Fisher-Yates shuffle in the Rand monad
shuffleVector :: VG.Vector v a => v a -> Rand (v a)
shuffleVector v =
  Rand $ \s -> case shuffleVector' s v of (sv, s') -> R sv s'

-- Fisher-Yates shuffle
shuffleVector' :: VG.Vector v a => PureMT -> v a -> (v a, PureMT)
shuffleVector' gen v = runST $ do
  let maxIdx = VG.length v - 1
  vm   <- VG.thaw v
  gen' <- swaps vm 0 maxIdx gen
  vmf  <- VG.unsafeFreeze vm
  return (vmf, gen')
  where
    swaps vm idx maxIdx gen'
      | idx < maxIdx = do
          let (newIdx, gen'') = randomIntR gen' (idx, maxIdx)
          GM.unsafeSwap vm idx newIdx
          swaps vm (idx + 1) maxIdx gen''
      | otherwise = return gen'

-- |
-- A test stastic calculates the difference between two samples. See
-- 'meanDifference' and 'varianceRatio' for examples.
type TestStatistic = Sample -> Sample -> Double

-- |
-- Calculates the difference mean of two samples (/mean(s1 - s2)/). When the
-- two samples do not have an equal length, the trailing elements of the
-- longer vector are ignored.
differenceMean :: TestStatistic
differenceMean v1 v2 =
  VG.sum (subVector v1 v2) / fromIntegral (VG.length v1)

-- | Calculates the mean difference of two samples (/mean(s1) - mean(s2)/).
meanDifference :: TestStatistic
meanDifference s1 s2 =
  mean s1 - mean s2

-- | Calculate the mean of a sample.
mean :: Sample -> Double
mean = do
  t <- VG.sum
  l <- VG.length
  return $ t / fromIntegral l

-- | Calculate the ratio of sample variances (/var(s1) : var(s2)/).
varianceRatio :: TestStatistic
varianceRatio v1 v2 =
  variance v1 / variance v2

-- | Subtract two vectors.
subVector :: (VG.Vector v n, Num n) => v n -> v n -> v n
subVector = VG.zipWith (-)

subIIW :: Int -> Int -> Word
subIIW a b = fromIntegral a - fromIntegral b
{-# INLINE subIIW #-}

addIWI :: Int -> Word -> Int
addIWI a b = a + fromIntegral b
{-# INLINE addIWI #-}

-- | Generate Int numbers within a range
randomIntR :: PureMT -> (Int, Int) -> (Int, PureMT)
randomIntR gen (a, b)
  | n == 0    = randomInt gen
  | otherwise = loop gen
  where
    (a', b') = if a < b then (a, b) else (b, a)
    -- Number of different Ints that should be generated
    n = 1 + subIIW b' a'
    -- The total range of Word can hold x complete n ranges
    x = maxBound `div` n
    -- Pick from a range the is dividable by n without remainders
    s = x * n
    loop gen'
      | r >= s    = loop gen'' -- r is outside the range, discard it...
      | otherwise = (addIWI a' (r `div` x), gen'')
      where
        (!r, !gen'') = randomWord gen'
{-# INLINE randomIntR #-}

max_stars_repo_name: danieldk/approx-rand-test

max_issues_repo_name: danieldk/approx-rand-test

max_forks_repo_name: danieldk/approx-rand-test


text:
{-# LANGUAGE ScopedTypeVariables #-}

module AI.Learning.NeuralNetwork
    (
      -- * Representation
      NeuralNetwork
      -- * Prediction and Training
    , nnPredict
    , nnTrain
    , nnTrainIO
    ) where

import Control.Monad.Random hiding (fromList)
import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Util
import System.IO.Unsafe

import AI.Learning.Core
import AI.Util.Matrix

-----------------------
-- NN Representation --
----------------------

-- |Representation for a single-hidden-layer neural network. If the network
--  has K input nodes, H hidden nodes and L output layers then the dimensions
--  of the matrices theta0 and theta1 are
--
--  * size Layer1 = (K+1) x H
--  * size Layer2 = (H+1) x L
--
--  the (+1)s account for the addition of bias nodes in the input layer and
--  the hidden layer. Therefore the total number of parameters is
--  H * (K + L) + H + L
data NeuralNetwork = NN (Matrix Double) (Matrix Double)

-- |Three-tuple describing the shape of a network: (input, hidden, output).
type NNShape = (Int,Int,Int)

instance Show NeuralNetwork where
    show (NN theta0 theta1) = "Neural Net:\n\n" ++ dispf 3 theta0 ++
                                           "\n" ++ dispf 3 theta1

fromVector :: NNShape -> Vector Double -> NeuralNetwork
fromVector (k,h,l) vec = NN theta0 theta1
    where theta0 = reshape h $ takeVector ((k + 1) * h) vec
          theta1 = reshape l $ dropVector ((k + 1) * h) vec

toVector :: Matrix Double -> Matrix Double -> Vector Double
toVector theta0 theta1 = join [flatten theta0, flatten theta1]

----------------------
-- NN Train/Predict --
----------------------

-- |Make predictions using a neural network.
nnPredict :: NeuralNetwork -> Matrix Double -> Matrix Double
nnPredict nn x = h where (_,_,h) = nnForwardProp nn x

-- |Train a neural network from a training set. Note that you must supply
--  an initial vector of weights. Supplying initial weights all equal to
--  zero will generally give poor performance.
nnTrain :: Int              -- number of hidden neurons
        -> Matrix Double    -- y
        -> Matrix Double    -- x
        -> Double           -- lambda
        -> Vector Double    -- initial weights
        -> NeuralNetwork
nnTrain h y x lambda initialVec = fromVector shape vec
    where shape = (cols x, h, cols y)
          vec   = minimizeS cost grad initialVec
          cost  = fst . nnCostGradient shape y x lambda
          grad  = snd . nnCostGradient shape y x lambda

-- |Train a neural network from a training set.
--
--  Note that this is implemented as an IO action, because it randomly sets the
--  weights in the network before performing numerical optimization (this is to
--  break the symmetry that exists in the network when all weights are zero).
--  As a result, it is possible to get different results when training the
--  network multiple times with the same data.
nnTrainIO :: Int              -- number of hidden neurons
          -> Matrix Double    -- y
          -> Matrix Double    -- x
          -> Double           -- lambda
          -> IO NeuralNetwork
nnTrainIO h y x lambda = nnTrain h y x lambda `fmap` initialVec (cols x, h, cols y)

-- |Choose initial random weights for a neural network.
initialVec :: NNShape -> IO (Vector Double)
initialVec (k,h,l) = do
    let len = h * (k + l) + h + l
    xs <- getRandomRs (0.0, 0.01)
    return . fromList $ take len xs

------------------------------
-- Forward/Back Propagation --
------------------------------

-- |Perform forward propagation through a neural network, returning the matrices
--  created in the process.
nnForwardProp :: NeuralNetwork                                 -- neural net
              -> Matrix Double                                 -- design matrix (x)
              -> (Matrix Double, Matrix Double, Matrix Double) -- results of fwd prop
nnForwardProp (NN theta0 theta1) x = (a0,a1,a2)
    where a0 = addOnes $ x
          a1 = addOnes $ sigmoid (a0 <> theta0)
          a2 = sigmoid (a1 <> theta1)

-- |Perform backward propagiation through a neural network. You must supply the
--  target values and the results of forward propagation for each layer, and
--  the function returns the gradient matrices for the neural network.
nnBackProp :: NeuralNetwork                                 -- neural net
           -> Matrix Double                                 -- target (y)
           -> (Matrix Double, Matrix Double, Matrix Double) -- results of fwd prop
           -> (Matrix Double, Matrix Double)                -- gradient (delta0, delta1)
nnBackProp (NN _ theta1) y (a0,a1,a2) = (dropColumns 1 delta0, delta1)
    where
        d2     = a2 - y
        d1     = (d2 <> trans theta1) * a1 * (1 - a1)
        delta0 = trans a0 <> d1
        delta1 = trans a1 <> d2

-- |Perform back and forward propagation through a neural network, returning the
--  final predictions (variable /a2/) and the gradient matrices (variables
--  /delta0/ and /delta1/) produced.
nnFwdBackProp :: NeuralNetwork -> Matrix Double -> Matrix Double -> (Matrix Double, Matrix Double, Matrix Double)
nnFwdBackProp nn@(NN theta0 theta1) y x = (a2, delta0, delta1)
    where
        (a0,a1,a2)      = nnForwardProp nn x
        (delta0,delta1) = nnBackProp nn y (a0,a1,a2)

-- |Compute the penalty function and gradient vector for a neural network
--  given a training set.
nnCostGradient :: NNShape                   -- (K,H,L)
               -> Matrix Double             -- targets (y)
               -> Matrix Double             -- design matrix (x)
               -> Double                    -- regularization parameter (lambda)
               -> Vector Double             -- neural network
               -> (Double, Vector Double)   -- (cost, gradient)
nnCostGradient shape y x lambda vec = (cost, grad)
    where
        m = fromIntegral (rows x)
        nn@(NN theta0 theta1) = fromVector shape vec
        (h, delta0, delta1)   = nnFwdBackProp nn y x

        cost  = (cost1 + cost2) / m
        cost1 = negate $ sumMatrix $ y * log h + (1-y) * log (1-h)
        cost2 = lambda/2 * (normMatrix theta0 + normMatrix theta1)

        grad  = (1/m) `scale` (grad1 + grad2)
        grad1 = toVector delta0 delta1
        grad2 = lambda `scale` toVector (insertNils theta0) (insertNils theta1)

        normMatrix m = sumMatrix $ (dropRows 1 m) ^ 2
        insertNils m = vertcat [0, dropRows 1 m]

nnCost shape y x lambda = fst . nnCostGradient shape y x lambda
nnGrad shape y x lambda = snd . nnCostGradient shape y x lambda

-- |Use central differencing to compute an approximation to the gradient
--  vector for a neural network. This is mainly used for checking the
--  implementation of backprop.
nnGradApprox :: NNShape -> Matrix Double -> Matrix Double -> Double -> Vector Double -> Vector Double
nnGradApprox shape y x lambda vec = fromList $ g `map` [0..n-1]
    where
        h = 1e-6
        n = dim vec
        f v = nnCost shape y x lambda v
        g i = (f (vec + e i) - f (vec - e i)) / (2*h)
        e i = fromList $ replicate i 0 ++ [h] ++ replicate (n-i-1) 0

-------------
-- Testing --
-------------

testNN :: NeuralNetwork
testNN = NN t0 t1
    where t0 = fromLists [[11.9934, -5.1396], [-7.7162, 10.1512], [-7.668, 10.1835]]
          t1 = fromLists [[-16.8806], [10.0445], [8.7476]]

testFwdProp :: IO ()
testFwdProp = do
    putStrLn "***\nCompare forward propagation to the MATLAB implementation.\n"
    let theta0 = fromLists [[0], [-10]]
        theta1 = fromLists [[5],[-10]]
        nn = NN theta0 theta1
        x = fromLists [[-1],[0],[1]]
    -- sigmoid [10, 0, -10] = [1, 0.5, 0]
    -- sigmoid [-5, 0, 5]  = [0.0, 0.5, 1]
    let y = nnPredict nn x
    putStrLn "Predictions (should be roughly 0.0, 0.5, 1.0)"
    disp 2 y

testBackProp :: IO ()
testBackProp = do
    putStrLn "***\nCompare back propagation to the MATLAB implementation.\n"
    let nn = testNN
        x  = fromLists [[0, 0], [0, 1], [1, 0], [1, 1]]
        y  = fromLists [[0], [1], [1], [0]]
        (h, delta0, delta1) = nnFwdBackProp nn y x
    putStrLn "Predictions (should be 0.0011, 0.8487, 0.8476, 0.0004)"
    disp 4 h
    putStrLn "Delta 0 (should be -0.0401, -0.0171, -0.0205, -0.0088, -0.0195, -0.0084)"
    disp 4 delta0
    putStrLn "Delta 1 (should be -0.3022, -0.2985, -0.3013)"
    disp 4 delta1

testCostGradient :: IO ()
testCostGradient = do
    putStrLn "***\nTest cost/gradient against the MATLAB implementation.\n"
    let nn@(NN theta0 theta1) = testNN
        x = fromLists [[0, 0], [0, 1], [1, 0], [1, 1]]
        y = fromLists [[0], [1], [1], [0]]
        lambda = 1e-4
        (cost, grad) = nnCostGradient (2,2,1) y x lambda (toVector theta0 theta1)
    putStrLn "Cost (should be around 0.0890)"
    print cost
    putStrLn "Gradient (should be -0.0100, -0.0043, -0.0053, -0.0019, -0.0051, -0.0019, -0.0755, -0.0744, -0.0751)"
    disp 4 (asColumn grad)

test :: Int -> Double -> IO ()
test n lambda = do
    putStrLn "***\nLearning XOR function.\n"
    x <- rand n 2
    e <- fmap (0.01*) (rand n 1)
    let y = mapMatrix (\x -> if x > 0.5 then 1.0 else 0.0) (xor x)
    nn <- nnTrainIO 4 y x lambda
    let ypred = nnPredict nn x
    -- Show predictions
    putStrLn "Predictions:"
    disp 2 $ takeRows 10 $ horzcat [x, y, ypred]
    -- Show neural net
    print nn
    -- Final test; should approximately compute xor function
    let xx = fromLists [[0,0],[0,1],[1,0],[1,1]]
        yy = nnPredict nn xx
    putStrLn "Exclusive or:"
    disp 2 $ horzcat [xx,yy]

xor :: Matrix Double -> Matrix Double
xor x = let [u,v] = toColumns x in asColumn (u + v - 2 * u * v)

max_stars_repo_name: cagix/aima-haskell


text:
module Widget.Total
  ( Total()
  , mkTotal
  ) where

import Import
import File
import Format.SplitTime
import Format.Percentile

import Statistics.Distribution.Empirical

data Total = Total FileFormat (Zipper EmpiricalDistribution) (Zipper NominalDiffTime) Window

instance Widget Total where
  update                 = updateTotal
  handle                 = handleTotal
  redraw                 = redrawTotal
  window (Total _ _ _ w) = w

mkTotal :: FileFormat -> [EmpiricalDistribution] -> Window -> Total
mkTotal f p = Total f (listToZipper p) (listToZipper [0])

updateTotal :: NominalDiffTime -> Total -> Curses Total
updateTotal d (Total f c z w) = updateWindow w $ do
  let nz = replace (curs z + d) z
      nt = Total f c nz w
      str = show $ splitTime (sum $ zipperToList nz)

  (rows, columns) <- windowSize

  moveCursor 1 (columns - 14)
  drawString $
    if null $ left z
      then "----"
      else if null $ right c
             then show $ percentile (curs c) (sum $ zipperToList z)
             else show $ percentile (curs $ prev c) (sum $ left z)

  moveCursor 1 (columns - (fromIntegral $ length str))
  insertString str

  return nt

handleTotal :: TimerAction -> Total -> Curses Total
handleTotal Reset   (Total f c _ w) = return $ mkTotal f (zipperToList c) w
handleTotal Advance (Total f c z w) = return $ Total f (next c) (push 0 z) w
handleTotal Reverse (Total f c z w) =
  return $
    if null $ left z
      then Total f (prev c) z w
      else
        let cz = curs z
            z' = trunc $ prev z
            nz = replace (curs z' + cz) z'
        in Total f (prev c) nz w
handleTotal Skip    (Total f c z w) = return $ Total f (next c) (push 0 z) w
handleTotal _ t = return t

redrawTotal :: Total -> Curses ()
redrawTotal t@(Total f c z w) = updateWindow w $ do
  (rows, columns) <- windowSize

  erase

  -- Top horizontal line
  -- This is the only thing in row 0
  moveCursor 0 0
  drawLineH Nothing columns

  -- Current total
  -- Tag
  moveCursor 1 (columns - 22)
  drawString "Total:"

  -- Percentile
  moveCursor 1 (columns - 14)
  insertString $
    if null $ left z
      then "----"
      else if null $ right c
             then show $ percentile (curs c) (sum $ zipperToList z)
             else show $ percentile (curs $ prev c) (sum $ left z)

  -- Current cumulative time over all splits
  do
     let str = show $ splitTime $ sum $ zipperToList z

     moveCursor 1 (columns - (fromIntegral $ length str))
     insertString str

  -- The following are consistent internal to a run
  let target = if null $ right c
                 then curs c
                 else last $ right c

  -- Personal Best
  -- Tag
  moveCursor 2 (columns - 30)
  drawString "Personal Best:"

  do
    let pb  = personalBest f
        str = denoteZero $ splitTime pb
        p   = percentile target pb

    -- Percentile
    moveCursor 2 (columns - 14)
    insertString $ ignoreHigh p

    -- Current personal best
    moveCursor 2 (columns - (fromIntegral $ length str))
    insertString str

  -- Draw Sum of Bests
  -- Tag
  moveCursor 3 (columns - 29)
  drawString "Sum of Bests:"

  do
    let sob = sumOfBests f
        str = denoteZero $ splitTime $ sob
        p   = percentile target sob

    -- Percentile
    moveCursor 3 (columns - 14)

    if (sob > 0)
      then insertString $ ignoreHigh p
      else insertString "----"

    -- Current sum of bests
    moveCursor 3 (columns - (fromIntegral $ length str))
    insertString str

  -- Various vertical lines
  moveCursor 0 (columns - 15)
  drawGlyph glyphPlus

  moveCursor 1 (columns - 15)
  drawLineV Nothing 3

  moveCursor 0 (columns - 10)
  drawGlyph glyphPlus

  moveCursor 1 (columns - 10)
  drawLineV Nothing 3

max_stars_repo_name: argiopetech/timer

max_issues_repo_name: argiopetech/timer

max_forks_repo_name: argiopetech/timer


text:
{-# LANGUAGE ViewPatterns #-}

module Edge where

import Codec.Image.DevIL
import Data.Array.Unboxed
import Math.Probably.MCMC
import Math.Probably.Sampler
import Math.Probably.FoldingStats
import Control.Monad.State.Strict
import System.Cmd
import System.Environment
import Data.Array.IO
import System.IO

import Data.Array.Unboxed
import Numeric.LinearAlgebra hiding (find)
import qualified Math.Probably.PDF as PDF
import Data.List hiding (map)
import Data.Maybe
import Control.Applicative
import Data.Ord
import qualified Data.Vector as V
import qualified Graphics.Rendering.OpenGL.GLU.Tessellation as GLU
import Graphics.Rendering.OpenGL

import CVUtils
import qualified Graphics.UI.GLUT as GLUT

loadPoints :: IO [(R,R)]
loadPoints = fmap read $ readFile "edge.dat"

complexPolygon :: [(R,R)] -> GLU.ComplexPolygon GLfloat
complexPolygon points
  = let p2v (x,y) = Vertex3 (realToFrac x) (realToFrac y) 0
    in GLU.ComplexPolygon
         [GLU.ComplexContour $map (\v-> GLU.AnnotatedVertex (p2v v) 0) points]

noOpCombiner _newVertex _weightedProperties = 0.0 ::GLfloat

triangulateEdge :: [(R,R)] -> IO [[Vector R]]
triangulateEdge pts = do
   fmap getTriangles $ GLU.triangulate
            GLU.TessWindingPositive 0 (Normal3 0 0 0) noOpCombiner
            $ complexPolygon pts

type Triangles = [[Vector R]]

getTriangles (GLU.Triangulation tris) = map unTri tris where
   unTri
    (GLU.Triangle (GLU.AnnotatedVertex (Vertex3 x1 y1 _) _)
                  (GLU.AnnotatedVertex (Vertex3 x2 y2 _) _)
                  (GLU.AnnotatedVertex (Vertex3 x3 y3 _) _))
       = [fromList [realToFrac x1,realToFrac y1],
          fromList [realToFrac x2,realToFrac y2],
          fromList [realToFrac x3,realToFrac y3]]

--http://frame3dd.svn.sourceforge.net/viewvc/frame3dd/trunk/src/microstran/vec3.c?view=markup
cross a b = fromList [(a@>1)*(b@>2)-(a@>2)*(b@>1)
                     ,(a@>2)*(b@>0)-(a@>0)*(b@>2)
                     ,(a@>0)*(b@>1)-(a@>1)*(b@>0)]

crossZ x y = (x@>0)*(y@>1)-(x@>1)*(y@>0)

crossZp a b p = ((b@>0)-(a@>0))*((p@>1)-(a@>1))-((b@>1)-(a@>1))*((p@>0)-(a@>0))

--http://www.blackpawn.com/texts/pointinpoly/default.html

sameSide :: (Vector R,Vector R,Vector R,Vector R) -> Bool
sameSide(p1,p2, a,b)
 = let cp1 = crossZ (b-a) (p1-a)
       cp2 = crossZ (b-a) (p2-a)
    in cp1 * cp2 >= 0

sameSide' :: (Vector R,Vector R,Vector R,Vector R) -> Bool
sameSide'(p1,p2, a,b)
 = let cp1 = crossZp a b p1
       cp2 = crossZp a b p2
    in cp1 * cp2 >= 0

pointInTriangle :: Vector R -> [Vector R] -> Bool
pointInTriangle p [a,b,c]
    = sameSide'(p,a, b,c) && sameSide'(p,b, a,c)
        && sameSide'(p,c, a,b)

markBg :: Image -> IO Image
markBg im = do
    (progName,args) <-  GLUT.getArgsAndInitialize
    pts <- loadPoints
    tris <- triangulateEdge pts
 --   print tris
    print $ head tris
    --print $ cross (head tris!!0) (head tris!!1)

--    print $ pointInTriangle (fromList [100,100]) $ head tris
    mutIm <- thaw im
    ((loy,lox,_),(hiy,hix,_)) <- getBounds (mutIm::MImage)
    forM_ [(y,x,0) | x<- [lox..hix],
                     y<- [loy..hiy]] $ \ix@(y,x,_)-> do
      when (any (pointInTriangle (fromList [realToFrac x, realToFrac y])) tris) $ writeArray mutIm ix 255
    freeze (mutIm::MImage)
--    return im

loadTriangles :: IO [[Vector R]]
loadTriangles = do
    (progName,args) <-  GLUT.getArgsAndInitialize
    pts <- loadPoints
    triangulateEdge pts

writeVisMat :: String -> BitImage -> Image -> IO ()
writeVisMat nm vm im = do
   mutIm <- thaw im
   ((loy,lox,_),(hiy,hix,_)) <- getBounds (mutIm::MImage)
   forM_ [(y,x,0) | x<- [lox..hix],
                     y<- [loy..hiy]] $ \ix@(y,x,_)-> do
     when (not $ readBitImage vm x y) $ writeArray mutIm ix 255

   newIm <- freeze (mutIm::MImage)
   writeImage nm newIm
   return ()

max_stars_repo_name: glutamate/cvutils

max_issues_repo_name: glutamate/cvutils

max_forks_repo_name: glutamate/cvutils


text:
{-# LANGUAGE
  Rank2Types,
  TypeFamilies,
  TupleSections
 #-}

-- Import all models under maintenance.
-- Models not imported here will not be compiled
-- when invoking `stack bench`.
import qualified BetaBin
import qualified Dice
import qualified DPmixture
import qualified Gamma
import qualified HMM
import Plotting

import Numeric.LogDomain
import Control.Monad.Bayes.Class
import Control.Monad.Bayes.Simple
import Control.Monad.Bayes.Weighted
import Control.Monad.Bayes.Sampler
import Control.Monad.Bayes.Enumerator
import Control.Monad.Bayes.Inference
import Control.Monad.Bayes.Population

import System.IO
import Control.Arrow (first, second)
import Control.Monad (unless, when)
import Control.Monad.IO.Class (MonadIO, liftIO)

import Statistics.Sample
import qualified Data.Vector as Vector
import Graphics.Rendering.Chart.Easy
import Graphics.Rendering.Chart.Backend.Cairo
import Options.Applicative
import System.Directory

opts :: ParserInfo Bool
opts = flip info fullDesc $ switch
  ( long "trial"
 <> help "Run a quick version of benchmarks to check that all is working correctly")

tryCache :: (MonadIO m, Read a, Show a) => FilePath -> m a -> m a
tryCache filepath fresh = do
  exists <- liftIO $ doesFileExist filepath
  if exists then
    liftIO $ fmap read (readFile filepath)
  else do
    value <- fresh
    liftIO $ writeFile filepath (show value)
    return value

main = do
  -- make sure `putStrLn` prints to console immediately
  hSetBuffering stdout LineBuffering

  trial <- execParser opts
  when trial $ putStrLn "Trial run"

  sampleIO hmmBenchmark

meanVar :: (MonadDist m, CustomReal m ~ Double) => Int -> m Double -> m (Double, Double)
meanVar n d = do
  xs <- Vector.replicateM n d
  return $ meanVariance xs

smcParams :: [Int]
smcParams = [10,20,50,100,200,500,1000]

ns :: [Int]
--ns = [10,20,50,100,200,500,1000]
ns = [10,20..1000]

smcParamsDouble :: [Double]
smcParamsDouble = map fromIntegral smcParams

smcResults :: (MonadDist m, CustomReal m ~ Double) => [m (Vector.Vector Double)]
smcResults = map (\p -> Vector.replicateM 10 $ fmap HMM.hmmKL $ explicitPopulation $ smcMultinomial (length HMM.values) p HMM.hmm) smcParams

hmmBenchmark :: SamplerIO ()
hmmBenchmark = do
  liftIO $ putStrLn "running HMM benchmark"

  isSamples <- fmap (drop 5000) $ explicitPopulation $ importance 10000 HMM.hmm
  let isRes = map (\n -> HMM.hmmKL $ take n isSamples) ns
  -- mhSamples <- fmap (drop 5000) $ traceMH 10000 HMM.hmm
  -- let mhRes = map (\n -> HMM.hmmKL $ take n $ map (,1) mhSamples) ns
  mhPriorSamples <- fmap (drop 5000) $ mhPrior HMM.hmm 10000
  let mhPriorRes = map (\n -> HMM.hmmKL $ take n $ map (,1) mhPriorSamples) ns
  pimhSamples <- pimh (length HMM.values) 100 1000 HMM.hmm
  let pimhRes = map (\n -> HMM.hmmKL $ take n $ map (,1) pimhSamples) ns
  liftIO $ toFile (fo_format .~ PDF $ def) "anytime.pdf" $ do
    layout_title .= "HMM"
    anytimePlot "#samples" "KL" ns [
      ("IS", isRes),
      -- ("MHtrace", mhRes),
      ("MHprior", mhPriorRes),
      ("PIMH", pimhRes)]

  smcRes <- sequence smcResults
  liftIO $ toFile (fo_format .~ PDF $ def) "smc.pdf" $ do
    layout_title .= "HMM"
    oneShotPlot "#particles" "KL" [
      ("SMC", zip smcParamsDouble (map (errBars 2) smcRes))]

max_stars_repo_name: eligottlieb/monad-bayes

max_issues_repo_name: eligottlieb/monad-bayes

max_forks_repo_name: eligottlieb/monad-bayes