text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE Trustworthy #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE ExplicitNamespaces #-}
{-# OPTIONS_GHC -fno-warn-unused-imports #-}

module Protolude (
  -- * Base functions
  module Base,
  identity,
  pass,
#if !MIN_VERSION_base(4,8,0)
  (&),
  scanl',
#endif
  -- * Function functions
  module Function,
  applyN,
  -- * List functions
  module List,
  map,
  uncons,
  unsnoc,
  -- * Data Structures
  module DataStructures,
  -- * Show functions
  module Show,
  show,
  print,
  -- * Bool functions
  module Bool,
  -- * Monad functions
  module Monad,
  liftIO1,
  liftIO2,
  -- * Functor functions
  module Functor,
  -- * Either functions
  module Either,
  -- * Applicative functions
  module Applicative,
  guarded,
  guardedA,
  -- * String conversion
  module ConvertText,
  -- * Debug functions
  module Debug,

  -- * Panic functions
  module Panic,
  -- * Exception functions
  module Exception,
  Protolude.throwIO,
  Protolude.throwTo,
  -- * Semiring functions
  module Semiring,

  -- * String functions
  module String,
  -- * Safe functions
  module Safe,
  -- * Eq functions
  module Eq,
  -- * Ord functions
  module Ord,
  -- * Traversable functions
  module Traversable,
  -- * Foldable functions
  module Foldable,
  -- * Semigroup functions
#if MIN_VERSION_base(4,9,0)
  module Semigroup,
#endif
  -- * Monoid functions
  module Monoid,
  -- * Bifunctor functions
  module Bifunctor,
  -- * Bifunctor functions
  module Hashable,

  -- * Deepseq functions
  module DeepSeq,

  -- * Tuple functions
  module Tuple,

  module Typeable,

#if MIN_VERSION_base(4,7,0)
  -- * Typelevel programming
  module Typelevel,
#endif

  -- * Monads
  module Fail,
  module State,
  module Reader,
  module Except,
  module Trans,
  module ST,
  module STM,

  -- * Integers
  module Int,
  module Bits,

  -- * Complex functions
  module Complex,

  -- * Char functions
  module Char,

  -- * Maybe functions
  module Maybe,

  -- * Generics functions
  module Generics,

  -- * ByteString functions
  module ByteString,
  LByteString,

  -- * Text functions
  module Text,
  LText,

  -- * Read functions
  module Read,
  readMaybe,
  readEither,

  -- * System functions
  module System,
  die,

  -- * Concurrency functions
  module Concurrency,

  -- * Foreign functions
  module Foreign,
) where

-- Protolude module exports.
import Protolude.Debug as Debug
import Protolude.List as List
import Protolude.Show as Show
import Protolude.Bool as Bool
import Protolude.Monad as Monad
import Protolude.Functor as Functor
import Protolude.Either as Either
import Protolude.Applicative as Applicative
import Protolude.ConvertText as ConvertText
import Protolude.Panic as Panic
import Protolude.Exceptions as Exception
import Protolude.Semiring as Semiring
import qualified Protolude.Conv as Conv

import Protolude.Base as Base hiding (
    putStr           -- Overriden by Show.putStr
  , putStrLn         -- Overriden by Show.putStrLn
  , print            -- Overriden by Protolude.print
  , show             -- Overriden by Protolude.show
  , showFloat        -- Custom Show instances deprecated.
  , showList         -- Custom Show instances deprecated.
  , showSigned       -- Custom Show instances deprecated.
  , showSignedFloat  -- Custom Show instances deprecated.
  , showsPrec        -- Custom Show instances deprecated.
  )
import qualified Protolude.Base as PBase

-- Used for 'show', not exported.
import Data.String (String)
import Data.String as String (IsString)

-- Maybe'ized version of partial functions
import Protolude.Safe as Safe (
    headMay
  , headDef
  , initMay
  , initDef
  , initSafe
  , tailMay
  , tailDef
  , tailSafe
  , lastDef
  , lastMay
  , foldr1May
  , foldl1May
  , foldl1May'
  , maximumMay
  , minimumMay
  , maximumDef
  , minimumDef
  , atMay
  , atDef
  )

-- Applicatives
import Control.Applicative as Applicative (
    Applicative(..)
  , Alternative(..)
  , Const(Const,getConst)
  , ZipList(ZipList,getZipList)
  , (<**>)
  , liftA
  , liftA2
  , liftA3
  , optional
  )

-- Base typeclasses
import Data.Eq as Eq (
    Eq(..)
  )
import Data.Ord as Ord (
    Ord(..)
  , Ordering(LT,EQ,GT)
  , Down(Down)
  , comparing
  )
import Data.Traversable as Traversable
import Data.Foldable as Foldable (
    Foldable,
    fold,
    foldMap,
    foldr,
    foldr',
    foldl,
    foldl',
    toList,
#if MIN_VERSION_base(4,8,0)
    null,
    length,
#endif
    elem,
    maximum,
    minimum,
    foldrM,
    foldlM,
    traverse_,
    for_,
    mapM_,
    forM_,
    sequence_,
    sequenceA_,
    asum,
    msum,
    concat,
    concatMap,
    and,
    or,
    any,
    all,
    maximumBy,
    minimumBy,
    notElem,
    find,
  )
import Data.Functor.Identity as Functor (
    Identity(Identity, runIdentity)
  )

#if MIN_VERSION_base(4,9,0)
import Data.List.NonEmpty as List (
    NonEmpty((:|))
  , nonEmpty
  )
import Data.Semigroup as Semigroup (
    Semigroup(sconcat, stimes)
  , WrappedMonoid
  , diff
  , cycle1
  , stimesMonoid
  , stimesIdempotent
  , stimesIdempotentMonoid
  , mtimesDefault
  )
#endif

#if !MIN_VERSION_base(4,15,0)
import Data.Semigroup as Semigroup (
  , Option(..)
  , option
  )
#endif

import Data.Monoid as Monoid

#if !MIN_VERSION_base(4,8,0)
import Protolude.Bifunctor as Bifunctor (Bifunctor(bimap, first, second))
#else
import Data.Bifunctor as Bifunctor (Bifunctor(bimap, first, second))
#endif

-- Deepseq
import Control.DeepSeq as DeepSeq (
    NFData(..)
  , ($!!)
  , deepseq
  , force
  )

-- Data structures
import Data.Tuple as Tuple (
    fst
  , snd
  , curry
  , uncurry
  , swap
  )

import Data.List as List (
    splitAt
  , break
  , intercalate
  , isPrefixOf
  , drop
  , filter
  , reverse
  , replicate
  , take
  , sortBy
  , sort
  , intersperse
  , transpose
  , subsequences
  , permutations
  , scanl
#if MIN_VERSION_base(4,8,0)
  , scanl'
#endif
  , scanr
  , iterate
  , repeat
  , cycle
  , unfoldr
  , takeWhile
  , dropWhile
  , group
  , inits
  , tails
  , zipWith
  , zip
  , unzip
  , genericLength
  , genericTake
  , genericDrop
  , genericSplitAt
  , genericReplicate
  )

#if !MIN_VERSION_base(4,8,0)
-- These imports are required for the scanl' rewrite rules
import GHC.Exts (build)
import Data.List (tail)
#endif

-- Hashing
import Data.Hashable as Hashable (
    Hashable
  , hash
  , hashWithSalt
  , hashUsing
  )

import Data.Map as DataStructures (Map)
import Data.Set as DataStructures (Set)
import Data.Sequence as DataStructures (Seq)
import Data.IntMap as DataStructures (IntMap)
import Data.IntSet as DataStructures (IntSet)

import Data.Typeable as Typeable (
    TypeRep
  , Typeable
  , typeOf
  , cast
  , gcast
#if MIN_VERSION_base(4,7,0)
  , typeRep
  , eqT
#endif
  )

#if MIN_VERSION_base(4,7,0)
import Data.Proxy as Typelevel (
    Proxy(..)
  )

import Data.Type.Coercion as Typelevel (
    Coercion(..)
  , coerceWith
  , repr
  )

import Data.Type.Equality as Typelevel (
    (:~:)(..)
  , type (==)
  , sym
  , trans
  , castWith
  , gcastWith
  )

#endif

#if MIN_VERSION_base(4,8,0)
import Data.Void as Typelevel (
    Void
  , absurd
  , vacuous
  )
#endif

import Control.Monad.Fail as Fail (
    MonadFail
  )

-- Monad transformers
import Control.Monad.State as State (
    MonadState
  , State
  , StateT(StateT)
  , put
  , get
  , gets
  , modify
  , state
  , withState

  , runState
  , execState
  , evalState

  , runStateT
  , execStateT
  , evalStateT
  )

import Control.Monad.Reader as Reader (
    MonadReader
  , Reader
  , ReaderT(ReaderT)
  , ask
  , asks
  , local
  , reader
  , runReader
  , runReaderT
  )

import Control.Monad.Trans.Except as Except (
    throwE
  , catchE
  )

import Control.Monad.Except as Except (
    MonadError
  , Except
  , ExceptT(ExceptT)
  , throwError
  , catchError
  , runExcept
  , runExceptT
  , mapExcept
  , mapExceptT
  , withExcept
  , withExceptT
  )

import Control.Monad.Trans as Trans (
    MonadIO
  , lift
  , liftIO
  )

-- Base types
import Data.Int as Int (
    Int
  , Int8
  , Int16
  , Int32
  , Int64
  )
import Data.Bits as Bits (
  Bits,
  (.&.),
  (.|.),
  xor,
  complement,
  shift,
  rotate,
#if MIN_VERSION_base(4,7,0)
  zeroBits,
#endif
  bit,
  setBit,
  clearBit,
  complementBit,
  testBit,
#if MIN_VERSION_base(4,7,0)
  bitSizeMaybe,
#endif
  bitSize,
  isSigned,
  shiftL,
  shiftR,
  rotateL,
  rotateR,
  popCount,
#if MIN_VERSION_base(4,7,0)
  FiniteBits,
  finiteBitSize,
  bitDefault,
  testBitDefault,
  popCountDefault,
#endif
#if MIN_VERSION_base(4,8,0)
  toIntegralSized,
  countLeadingZeros,
  countTrailingZeros,
#endif
  )
import Data.Word as Bits (
    Word
  , Word8
  , Word16
  , Word32
  , Word64
#if MIN_VERSION_base(4,7,0)
  , byteSwap16
  , byteSwap32
  , byteSwap64
#endif
  )

import Data.Either as Either (
    Either(Left,Right)
  , either
  , lefts
  , rights
  , partitionEithers
#if MIN_VERSION_base(4,7,0)
  , isLeft
  , isRight
#endif
  )

import Data.Complex as Complex (
    Complex((:+))
  , realPart
  , imagPart
  , mkPolar
  , cis
  , polar
  , magnitude
  , phase
  , conjugate
  )
import Data.Char as Char (
    Char
  , ord
  , chr
  , digitToInt
  , intToDigit
  , toUpper
  , toLower
  , toTitle
  , isAscii
  , isLetter
  , isDigit
  , isHexDigit
  , isPrint
  , isAlpha
  , isAlphaNum
  , isUpper
  , isLower
  , isSpace
  , isControl
  )
import Data.Bool as Bool (
  Bool(True, False),
  (&&),
  (||),
  not,
  otherwise
  )
import Data.Maybe as Maybe (
    Maybe(Nothing, Just)
  , maybe
  , isJust
  , isNothing
  , fromMaybe
  , listToMaybe
  , maybeToList
  , catMaybes
  , mapMaybe
  )

import Data.Function as Function (
    const
  , (.)
  , ($)
  , flip
  , fix
  , on
#if MIN_VERSION_base(4,8,0)
  , (&)
#endif
  )

-- Genericss
import GHC.Generics as Generics (
    Generic(..)
  , Generic1
  , Rep
  , K1(..)
  , M1(..)
  , U1(..)
  , V1
  , D1
  , C1
  , S1
  , (:+:)(..)
  , (:*:)(..)
  , (:.:)(..)
  , Rec0
  , Constructor(..)
  , Datatype(..)
  , Selector(..)
  , Fixity(..)
  , Associativity(..)
#if ( __GLASGOW_HASKELL__ >= 800 )
  , Meta(..)
  , FixityI(..)
  , URec
#endif
  )

-- ByteString
import qualified Data.ByteString.Lazy
import Data.ByteString as ByteString (ByteString)

-- Text
import Data.Text as Text (
    Text
  , lines
  , words
  , unlines
  , unwords
  )
import qualified Data.Text.Lazy

import Data.Text.IO as Text (
    getLine
  , getContents
  , interact
  , readFile
  , writeFile
  , appendFile
  )

import Data.Text.Lazy as Text (
    toStrict
  , fromStrict
  )

import Data.Text.Encoding as Text (
    encodeUtf8
  , decodeUtf8
  , decodeUtf8'
  , decodeUtf8With
  )

import Data.Text.Encoding.Error as Text (
    OnDecodeError
  , OnError
  , UnicodeException
  , lenientDecode
  , strictDecode
  , ignore
  , replace
  )

-- IO
import System.Environment as System (getArgs)
import qualified System.Exit
import System.Exit as System (
    ExitCode(..)
  , exitWith
  , exitFailure
  , exitSuccess
  )
import System.IO as System (
    Handle
  , FilePath
  , IOMode(..)
  , stdin
  , stdout
  , stderr
  , withFile
  , openFile
  )

-- ST
import Control.Monad.ST as ST (
    ST
  , runST
  , fixST
  )

-- Concurrency and Parallelism
import Control.Exception as Exception (
    Exception,
    toException,
    fromException,
#if MIN_VERSION_base(4,8,0)
    displayException,
#endif
    SomeException(SomeException)
  , IOException
  , ArithException(
    Overflow,
    Underflow,
    LossOfPrecision,
    DivideByZero,
    Denormal,
    RatioZeroDenominator
    )
  , ArrayException(IndexOutOfBounds, UndefinedElement)
  , AssertionFailed(AssertionFailed)
#if MIN_VERSION_base(4,7,0)
  , SomeAsyncException(SomeAsyncException)
  , asyncExceptionToException
  , asyncExceptionFromException
#endif
  , AsyncException(StackOverflow, HeapOverflow, ThreadKilled, UserInterrupt)
  , NonTermination(NonTermination)
  , NestedAtomically(NestedAtomically)
  , BlockedIndefinitelyOnMVar(BlockedIndefinitelyOnMVar)
  , BlockedIndefinitelyOnSTM(BlockedIndefinitelyOnSTM)
#if MIN_VERSION_base(4,8,0)
  , AllocationLimitExceeded(AllocationLimitExceeded)
#endif
#if MIN_VERSION_base(4,10,0)
  , CompactionFailed(CompactionFailed)
#endif
  , Deadlock(Deadlock)
  , NoMethodError(NoMethodError)
  , PatternMatchFail(PatternMatchFail)
  , RecConError(RecConError)
  , RecSelError(RecSelError)
  , RecUpdError(RecUpdError)
#if MIN_VERSION_base(4,9,0)
  , ErrorCall(ErrorCall, ErrorCallWithLocation)
#else
  , ErrorCall(ErrorCall)
#endif
#if MIN_VERSION_base(4,9,0)
  , TypeError(TypeError)
#endif
  , ioError
  , catch
  , catches
  , Handler(Handler)
  , catchJust
  , handle
  , handleJust
  , try
  , tryJust
  , evaluate
  , mapException
  , mask
  , mask_
  , uninterruptibleMask
  , uninterruptibleMask_
  , MaskingState(..)
  , getMaskingState
#if MIN_VERSION_base(4,9,0)
  , interruptible
#endif
  , allowInterrupt
  , bracket
  , bracket_
  , bracketOnError
  , finally
  , onException
  )
import qualified Control.Exception as PException

import Control.Monad.STM as STM (
    STM
  , atomically
#if !(MIN_VERSION_stm(2,5,0))
  , always
  , alwaysSucceeds
#endif
  , retry
  , orElse
  , check
  , throwSTM
  , catchSTM
  )

import Control.Concurrent.MVar as Concurrency (
    MVar
  , newEmptyMVar
  , newMVar
  , takeMVar
  , putMVar
  , readMVar
  , swapMVar
  , tryTakeMVar
  , tryPutMVar
  , isEmptyMVar
  , withMVar
#if MIN_VERSION_base(4,7,0)
  , withMVarMasked
#endif
  , modifyMVar_
  , modifyMVar
  , modifyMVarMasked_
  , modifyMVarMasked
#if MIN_VERSION_base(4,7,0)
  , tryReadMVar
  , mkWeakMVar
#endif
  , addMVarFinalizer
  )
import Control.Concurrent.Chan as Concurrency (
    Chan
  , newChan
  , writeChan
  , readChan
  , dupChan
  , getChanContents
  , writeList2Chan
  )
import Control.Concurrent.QSem as Concurrency (
    QSem
  , newQSem
  , waitQSem
  , signalQSem
  )
import Control.Concurrent.QSemN as Concurrency (
    QSemN
  , newQSemN
  , waitQSemN
  , signalQSemN
  )
import Control.Concurrent as Concurrency (
    ThreadId
  , forkIO
  , forkFinally
  , forkIOWithUnmask
  , killThread
  , forkOn
  , forkOnWithUnmask
  , getNumCapabilities
  , setNumCapabilities
  , threadCapability
  , yield
  , threadDelay
  , threadWaitRead
  , threadWaitWrite
#if MIN_VERSION_base(4,7,0)
  , threadWaitReadSTM
  , threadWaitWriteSTM
#endif
  , rtsSupportsBoundThreads
  , forkOS
#if MIN_VERSION_base(4,9,0)
  , forkOSWithUnmask
#endif
  , isCurrentThreadBound
  , runInBoundThread
  , runInUnboundThread
  , mkWeakThreadId
  , myThreadId
  )
import Control.Concurrent.Async as Concurrency (
    Async(..)
  , Concurrently(..)
  , async
  , asyncBound
  , asyncOn
  , withAsync
  , withAsyncBound
  , withAsyncOn
  , wait
  , poll
  , waitCatch
  , cancel
  , cancelWith
  , asyncThreadId
  , waitAny
  , waitAnyCatch
  , waitAnyCancel
  , waitAnyCatchCancel
  , waitEither
  , waitEitherCatch
  , waitEitherCancel
  , waitEitherCatchCancel
  , waitEither_
  , waitBoth
  , link
  , link2
  , race
  , race_
  , concurrently
  )

import Foreign.Ptr as Foreign (IntPtr, WordPtr)
import Foreign.Storable as Foreign (Storable)
import Foreign.StablePtr as Foreign (StablePtr)

-- Read instances hiding unsafe builtins (read)
import qualified Text.Read as Read
import Text.Read as Read (
    Read
  , reads
  )

-- Type synonymss for lazy texts
type LText = Data.Text.Lazy.Text
type LByteString = Data.ByteString.Lazy.ByteString

#if !MIN_VERSION_base(4,8,0)
infixl 1 &

(&) :: a -> (a -> b) -> b
x & f = f x
#endif

-- | The identity function, returns the give value unchanged.
identity :: a -> a
identity x = x

map :: Functor.Functor f => (a -> b) -> f a -> f b
map = Functor.fmap

uncons :: [a] -> Maybe (a, [a])
uncons [] = Nothing
uncons (x:xs) = Just (x, xs)

unsnoc :: [x] -> Maybe ([x],x)
unsnoc = Foldable.foldr go Nothing
  where
    go x mxs = Just (case mxs of
       Nothing -> ([], x)
       Just (xs, e) -> (x:xs, e))

-- | Apply a function n times to a given value
applyN :: Int -> (a -> a) -> a -> a
applyN n f = Foldable.foldr (.) identity (List.replicate n f)

-- | Parse a string using the 'Read' instance.
-- Succeeds if there is exactly one valid result.
--
-- >>> readMaybe ("123" :: Text) :: Maybe Int
-- Just 123
--
-- >>> readMaybe ("hello" :: Text) :: Maybe Int
-- Nothing
readMaybe :: (Read b, Conv.StringConv a String) => a -> Maybe b
readMaybe = Read.readMaybe . Conv.toS

-- | Parse a string using the 'Read' instance.
-- Succeeds if there is exactly one valid result.
-- A 'Left' value indicates a parse error.
--
-- >>> readEither "123" :: Either Text Int
-- Right 123
--
-- >>> readEither "hello" :: Either Text Int
-- Left "Prelude.read: no parse"
readEither :: (Read a, Conv.StringConv String e, Conv.StringConv e String) => e -> Either e a
readEither = first Conv.toS . Read.readEither . Conv.toS

-- | The print function outputs a value of any printable type to the standard
-- output device. Printable types are those that are instances of class Show;
-- print converts values to strings for output using the show operation and adds
-- a newline.
print :: (Trans.MonadIO m, PBase.Show a) => a -> m ()
print = liftIO . PBase.print

-- | Lifted throwIO
throwIO :: (Trans.MonadIO m, Exception e) => e -> m a
throwIO = liftIO . PException.throwIO

-- | Lifted throwTo
throwTo :: (Trans.MonadIO m, Exception e) => ThreadId -> e -> m ()
throwTo tid e = liftIO (PException.throwTo tid e)

-- | Do nothing returning unit inside applicative.
pass :: Applicative f => f ()
pass = pure ()

guarded :: (Alternative f) => (a -> Bool) -> a -> f a
guarded p x = Bool.bool empty (pure x) (p x)

guardedA :: (Functor.Functor f, Alternative t) => (a -> f Bool) -> a -> f (t a)
guardedA p x = Bool.bool empty (pure x) `Functor.fmap` p x

-- | Lift an 'IO' operation with 1 argument into another monad
liftIO1 :: MonadIO m => (a -> IO b) -> a -> m b
liftIO1 = (.) liftIO

-- | Lift an 'IO' operation with 2 arguments into another monad
liftIO2 :: MonadIO m => (a -> b -> IO c) -> a -> b -> m c
liftIO2 = ((.).(.)) liftIO

show :: (Show a, Conv.StringConv String b) => a -> b
show x = Conv.toS (PBase.show x)
{-# SPECIALIZE show :: Show  a => a -> Text  #-}
{-# SPECIALIZE show :: Show  a => a -> LText  #-}
{-# SPECIALIZE show :: Show  a => a -> String  #-}

#if MIN_VERSION_base(4,8,0)
-- | Terminate main process with failure
die :: Text -> IO a
die err = System.Exit.die (ConvertText.toS err)
#else
-- | Terminate main process with failure
die :: Text -> IO a
die err = hPutStrLn stderr err >> exitFailure
#endif

#if !MIN_VERSION_base(4,8,0)
-- This is a literal copy of the implementation in GHC.List in base-4.10.1.0.

-- | A strictly accumulating version of 'scanl'
{-# NOINLINE [1] scanl' #-}
scanl'           :: (b -> a -> b) -> b -> [a] -> [b]
scanl' = scanlGo'
  where
    scanlGo'           :: (b -> a -> b) -> b -> [a] -> [b]
    scanlGo' f !q ls    = q : (case ls of
                            []   -> []
                            x:xs -> scanlGo' f (f q x) xs)

{-# RULES
"scanl'"  [~1] forall f a bs . scanl' f a bs =
  build (\c n -> a `c` foldr (scanlFB' f c) (flipSeqScanl' n) bs a)
"scanlList'" [1] forall f a bs .
    foldr (scanlFB' f (:)) (flipSeqScanl' []) bs a = tail (scanl' f a bs)
 #-}

{-# INLINE [0] scanlFB' #-}
scanlFB' :: (b -> a -> b) -> (b -> c -> c) -> a -> (b -> c) -> b -> c
scanlFB' f c = \b g -> \x -> let !b' = f x b in b' `c` g b'

{-# INLINE [0] flipSeqScanl' #-}
flipSeqScanl' :: a -> b -> a
flipSeqScanl' a !_b = a
#endif


text:
{-# LANGUAGE UnicodeSyntax #-}
{-# LANGUAGE BangPatterns #-}

{-|
Module      : Statistics.FastBayes.Linear
Description : Bayesian linear regression via maximum marginal likelihood.
Copyright   : (c) Melinae, 2014
                  Chad Scherrer, 2014
License     : MIT
Maintainer  : chad.scherrer@gmail.com
Stability   : experimental
Portability : POSIX

This module gives an implementation of Bayesian linear regression, with the scale of the prior chosen by marginal likelihood.

The inputs for a Bayesian linear model are identical to those of a classical linear model, except that in addition to a design matrix and response, we must also specify a prior distribution on the weights and the noise. This leaves us with an open question of how these should be specified.

In his book /Pattern Recognition and Machine Learning/, Christopher Bishop provides details for an approach that simplifies the situation significantly, and allows for much faster inference. The structure of the linear model allows us to integrate the posterior over the weights, resulting in the /marginal likelihood/, expressed as a function of the prior precision and noise precision. This, in turn, can be easily optimized.
-}

module Statistics.FastBayes.Linear
  ( test
  , Fit
  , fit0
  , fit1
  , priorPrecision
  , noisePrecision
  , effectiveNumParameters
  , logEvidence
  , mapWeights
  , intercept
  , hessian
  , features
  , Predict
  , predict
  , studentize
  , studentizeM
  , inputs
  , fitted
  , fittedSD
  )
  where

import qualified Data.Vector.Storable as V
import Data.Packed.Matrix
import Numeric.LinearAlgebra

data Fit a = Fit
  { priorPrecision         :: Double             -- ^The precision (inverse variance) of the prior distribution, determined by maximizing the marginal likelihood
  , noisePrecision         :: Double             -- ^The precision (inverse variance) of the noise
  , effectiveNumParameters :: Double             -- ^The effective number of parameters in the model
  , logEvidence            :: Double             -- ^The log of the evidence, which is useful for model comparison (different features, same response)
  , mapWeights             :: Vector Double      -- ^The MAP (maximum a posteriori) values for the parameter weights
  , intercept              :: Double             -- ^The y-intercept
  , hessian                :: Matrix Double      -- ^The Hessian (matrix of second derivatives) for the posterior distribution
  , features               :: a → Vector Double  -- ^Extract features from a data point
  }

-- instance Show (Fit a) where

-- |@fit0 lim x y@ fits a Bayesian linear model to a design matrix @x@ and response vector @y@. This is an iterative algorithm, resulting in a sequence (list) of (α,β) values. Here α is the prior precision, and β is the noise precision. The @lim@ function passed in is used to specify how the limit of this sequence should be computed.
fit0 ::
  ([(Double, Double)] → (Double, Double)) -- ^How to take the limit of the (α,β) sequence. A simple approach is, /e.g./, @fit (!! 1000) x y@
  → Bool                                  -- ^Whether to studentize (translate and re-scale) the features
  → (a → Vector Double)                   -- ^The features
  → [a]                                   -- ^The input data
  → Vector Double                         -- ^The response vector
  → Fit a
fit0 lim student f xs y = Fit α β γ logEv m 0 h f'
  where
  (x, f') | not student = (fromRows $ map f xs,      f)
          | student     = (                 x0, fs . f)
            where
            (x0, fs) = studentizeM . fromRows . map f $ xs
  n = rows x
  p = cols x
  α0 = 1.0
  β0 = 1.0

  -- A vector of the eigenvalues of xtx
  (_,sqrtEigs,_) = compactSVD x
  eigs = V.map square sqrtEigs

  xtx = trans x <> x
  xty = trans x <> y
  getHessian a b = diag (V.replicate p a) + scale b xtx

  m = scale β $ h <\> xty

  go :: Double → Double → [(Double, Double)]
  go a0 b0 = (a0, b0) : go a b
    where
    h0 = getHessian a0 b0
    m0 = scale b0 $ h0 <\> xty
    c = V.sum $ V.map (\x' → x' / (a0 + x')) eigs
    a = c / (m0 <.> m0)
    b = recip $ (normSq $ y - x <> m0) / (fromIntegral n - c)

  γ = V.sum $ V.map (\x' → x' / (α + x')) eigs

  h = getHessian α β
  logEv = 0.5 *
    ( fromIntegral p * log α
    + fromIntegral n * log β
    - (β * normSq (y - x <> m) + α * (m <.> m))
    - logDetH
    - fromIntegral n * log (2*pi)
    )
    where
    (_,(logDetH, _)) = invlndet h

  (α, β) = lim $ go α0 β0

fit1 ::
  ([(Double, Double)] → (Double, Double)) -- ^How to take the limit of the (α,β) sequence. A simple approach is, /e.g./, @fit (!! 1000) x y@
  → Bool                                  -- ^Whether to studentize (translate and re-scale) the features
  → (a → Vector Double)                   -- ^The features
  → [a]                                   -- ^The input data
  → Vector Double                         -- ^The response vector
  → Fit a
fit1 lim student f xs y = f0 {intercept = μ}
  where
  μ = mean y
  n = V.length y
  f0 = fit0 lim student f xs (y - konst μ n)

square :: Double → Double
square x = x * x

normSq :: Vector Double → Double
normSq x = x <.> x

mean :: Vector Double -> Double
mean v = V.sum v / fromIntegral (V.length v)

studentizeM :: Matrix Double -> (Matrix Double, Vector Double → Vector Double)
studentizeM m = (x, f)
  where
  (xs, fs) = unzip . map studentize . toColumns $ m
  x = fromColumns xs
  f = fromList . zipWith ($) fs . toList

studentize :: Vector Double -> (Vector Double, Double → Double)
studentize v = (scale (1/σ) v0, f)
  where
  n   = V.length v
  μ   = mean v
  v0  = v - konst μ n
  σ   = sqrt $ normSq v0 / fromIntegral n
  f x = (x - μ) / σ

test = myFit 11
  where
  x = [1,1.1..pi]
  y = fromList $ map sin x
  f n x = fromList . map (x **) $ [1..n]
  myFit n = fit1 (!!100) True (f n) x y

data Predict a = Predict
  { inputs          :: [a]
  , fitted          :: [Double]
  , fittedSD        :: [Double]
  }

predict :: Fit a -> [a] -> Predict a
predict myFit xs = Predict xs ys (toList . V.map sqrt $ konst (1/β) n + takeDiag cov)
  where
  n = rows x
  α  = priorPrecision myFit
  β  = noisePrecision myFit
  w  = mapWeights     myFit
  y0 = intercept      myFit
  h  = hessian        myFit
  f  = features       myFit
  x  = fromRows $ map f xs
  cov = x <> (h <\> trans x)
  y = x <> w + konst y0 n
  ys = toList y

ext: hs

max_stars_repo_name: cscherrer/fastbayes

max_issues_repo_name: cscherrer/fastbayes

max_forks_repo_name: cscherrer/fastbayes


text:
module Malliavin where

import Data.Random.Distribution.Normal
import Numeric.LinearAlgebra.HMatrix as H
import Data.Ratio

data L2Map a = L2Map {eval :: a -> Double,
                      norm_l2 :: a -> Double}
l2_1 = L2Map {eval = (\t -> 1), norm_l2 = id}
l2id = L2Map {eval = id, norm_l2 = (\x -> x ^ 3 / 3)}

data L1Map a = L1Map {norm_l1 :: a -> Double}
l1_1 = L1Map {eval = (\t -> 1), norm_l1 = id}
l1id = L1Map {eval = id, norm_l1 = (\x -> x ^ 2 / 2)}

-- | 'hermite' is an infinite list of Hermite polynomials for given x
hermite :: (Enum a, Num a) => a -> [a]
hermite x = s
  where s@(_:ts) = 1 : x : zipWith3 (\hn2 hn1 n1 -> x * hn1 - n1 * hn2) s ts [1..]

-- | 'gaussianProcess' samples from Gaussian process
gaussianProcess :: Seed                   -- random state
                -> Int                    -- number of samples m
                -> Int                    -- number of dimensions n
                -> ((Int, Int) -> Double) -- function that maps indices of the matrix into dot products of its elements
                -> [Vector Double]        -- m n-th dimensional samples of Gaussian process
gaussianProcess seed 0 _ _ = []
gaussianProcess seed m 0 _ = replicate m $ vector []
gaussianProcess seed m n dotProducts = toRows $ gaussianSample seed m mean cov_matrix
  where mean = n |> repeat 0
        cov_matrix = H.sym $ (n><n) $ map (\i -> dotProducts (quot i n, rem i n)) [0..]

-- | 'second' function takes a list and gives each second element of it
second (x:y:xs) = y : second xs
second _ = []

-- | 'coinTossExpansion' is a Wiener chaos expansion for coin-toss r.v. to n-th element
coinTossExpansion :: Int    -- number of elements in the sum
                  -> Double -- gaussian random variable
                  -> Double -- the sum
coinTossExpansion n xi = sum (take n $ 0.5 : zipWith (*) fn (second $ hermite xi))
  where fn = 1.0 / (sqrt $ 2 * pi) : zipWith (\fn1 k -> -fn1 * k / (k ^ 2 + 3 * k + 2)) fn [1, 3..]

-- | 'coinTossSequence' is a coin-toss sequence of given size
coinTossSequence :: Seed  -- random state ω
                 -> Int   -- size of resulting sequence m
                 -> [Int] -- coin-toss sequence of size m
coinTossSequence seed m = map (round.coinTossExpansion 100) (toList nvec)
  where nvec = gaussianProcess seed m 1 (\(i,j) -> 1) !! 0

-- | 'mapOverInterval' map function f over the interval [0, T]
mapOverInterval :: (Fractional a) => Int      -- n, size of the output list
                                  -> a        -- T, end of the time interval
                                  -> (a -> b) -- f, function that maps from fractional numbers to some abstract space
                                  -> [b]      -- list of values f(t), t \in [0, T]
mapOverInterval n endT fn = [fn $ (endT * fromIntegral i) / fromIntegral (n - 1) | i <- [0..(n-1)]]

type ItoIntegral = Seed         -- ω, random state
                -> Int          -- n, sample size
                -> Double       -- T, end of the time interval
                -> L2Map Double -- h, L2-function
                -> [Double]     -- list of values sampled from Ito integral

-- | 'itoIntegral'' trajectory of Ito integral on the time interval [0, T]
itoIntegral' :: ItoIntegral
itoIntegral' seed n endT h = 0 : (toList $ gp !! 0)
  where gp = gaussianProcess seed 1 (n-1) (\(i, j) -> norm_l2 h $ fromIntegral (min i j + 1) * t)
        t = endT / fromIntegral n

-- | 'itoIntegral' faster implementation of itoIntegral' function
itoIntegral :: ItoIntegral
itoIntegral seed 0 _    _ = []
itoIntegral seed n endT h = scanl (+) 0 increments
  where increments = toList $ sigmas * gaussianVector
        gaussianVector = flatten $ gaussianSample seed (n-1) (vector [0]) (H.sym $ matrix 1 [1])
        sigmas = fromList $ zipWith (\x y -> sqrt(x-y)) (tail hnorms) hnorms
        hnorms = mapOverInterval n endT $ norm_l2 h

-- | 'brownianMotion' trajectory of Brownian motion a.k.a. Wiener process on the time interval [0, T]
brownianMotion :: Seed     -- ω, random state
               -> Int      -- n, sample size
               -> Double   -- T, end of the time interval
               -> [Double] -- list of values sampled from Brownian motion
brownianMotion seed n endT = itoIntegral seed n endT l2_1

type DiffusionProcess = Seed         -- ω, random state
                     -> Int          -- n, sample size
                     -> Double       -- T, end of the time interval
                     -> L1Map Double -- function \mu with L1 norm
                     -> L2Map Double -- function \sigma with L2 norm
                     -> [Double]     -- list of values sampled from stochastic diffusion process X(t)

-- | 'diffusionProcess' stochastic process, defined by SDE: dX(t) = \mu(t)dt + \sigma(t)dB(t)
diffusionProcess :: DiffusionProcess
diffusionProcess seed n endT mu sigma = zipWith (+) trend noise
  where trend = mapOverInterval n endT $ norm_l1 mu
        noise = itoIntegral seed n endT sigma

-- | 'multipleItoIntegral' multiple m-th Ito integral of m-th tensor power of h
multipleItoIntegral :: Seed         -- ω, random state
                    -> Int          -- m, tensor power
                    -> Int          -- n, sample size
                    -> Double       -- T, end of the time interval
                    -> L2Map Double -- function h with L2 norm
                    -> [Double]     -- list of values sampled from multiple m-th Wiener-Ito integral of h^(⊗m)
multipleItoIntegral seed m n endT h = zipWith fn itoTrajectory hnorm
  where itoTrajectory = itoIntegral seed n endT h
        hnorm = mapOverInterval n endT $ norm_l2 h
        fn :: Double -> Double -> Double
        fn _ 0 = 0
        fn x y = (y ^ m) * (hermite (x / y) !! m)

multipleBrownianMotion :: Seed     -- ω, random state
               -> Int      -- m, tensor power
               -> Int      -- n, sample size
               -> Double   -- T, end of the time interval
               -> [Double] -- list of values sampled from Brownian motion
multipleBrownianMotion seed m n endT = multipleItoIntegral seed m n endT l2_1

-- | 'totalVariation' calculates total variation of the sample process X(t)
totalVariation :: [Double] -- process trajectory X(t)
               -> Double   -- V(X(t))
totalVariation x@(_:tx) = sum $ zipWith (\x xm1 -> abs $ x - xm1) tx x

-- | 'quadraticVariation' calculates quadratic variation of the sample process X(t)
quadraticVariation :: [Double] -- process trajectory X(t)
                   -> Double   -- [X(t)]
quadraticVariation x@(_:tx) = sum $ zipWith (\x xm1 -> (x - xm1) ^ 2) tx x

max_stars_repo_name: vishalbelsare/Malliavin

max_issues_repo_name: vishalbelsare/Malliavin

max_forks_repo_name: vishalbelsare/Malliavin


text:
{-# LANGUAGE CPP                 #-}
{-# LANGUAGE DefaultSignatures   #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE FlexibleInstances   #-}
{-# LANGUAGE OverloadedStrings   #-}
{-# LANGUAGE RankNTypes          #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeOperators       #-}
-- | Module containing 'Homoiconic' and 'FromCode' type classes and its instance
-- declarations.
module Language.Finkel.Homoiconic
  ( -- * Homoiconic class
    Homoiconic(..)
  , fromCode
  , Result(..)

    -- * Generic functions
  , genericToCode
  , genericFromCode
  , genericParseCode

    -- * Generic classes
  , GToCode(..)
  , GParseCode(..)

    -- * Data.Data function
  , dataToCode
  ) where

#include "ghc_modules.h"

-- base
import           Control.Applicative             (Alternative (..))
import           Data.Complex                    (Complex (..))
import           Data.Data
import           Data.Fixed                      (Fixed (..))
import           Data.Functor.Compose            (Compose (..))
import           Data.Functor.Const              (Const (..))
import           Data.Functor.Identity           (Identity (..))
import           Data.Int                        (Int16, Int32, Int64, Int8)
import           Data.List.NonEmpty              (NonEmpty (..))
import           Data.Monoid                     (All (..), Alt (..), Any (..),
                                                  Dual (..), First (..),
                                                  Last (..), Product (..),
                                                  Sum (..))
import           Data.Ratio                      (Ratio, denominator, numerator,
                                                  (%))
import           Data.Version                    (Version (..))
import           Data.Word                       (Word16, Word32, Word64, Word8)
import           GHC.Generics                    (C, Constructor (..), D,
                                                  Generic (..), K1 (..),
                                                  M1 (..), S, U1 (..), V1,
                                                  (:*:) (..), (:+:) (..))
import           Numeric.Natural                 (Natural)

import qualified Data.Functor.Product            as Product
import qualified Data.Functor.Sum                as Sum
import qualified Data.Semigroup                  as Semigroup

#if !MIN_VERSION_ghc(8,8,0)
import           Control.Monad.Fail              (MonadFail (..))
import           Prelude                         hiding (fail)
#endif

-- ghc
import           GHC_Data_FastString             (FastString, unpackFS)
import           GHC_Types_SrcLoc                (GenLocated (..), SrcSpan,
                                                  getLoc)

-- Internal
import           Language.Finkel.Form

-- -------------------------------------------------------------------
--
-- Homoiconic type class
--
-- -------------------------------------------------------------------

-- | Class for handling Haskell value as 'Code'.
--
-- Instance of 'Homoiconic' should satisfy the law:
--
-- @
-- 'parseCode' ('toCode' x) ≡ 'Success' x
-- @
--
-- The function 'listToCode' and 'parseHsListCode' are used when handling
-- Haskell list values specially (e.g., 'Char'). These functions have default
-- implementations, which simply applies 'toCode' to elements of the argument
-- list, and which parses elements of 'HsList', respectively.
--
-- One can implement 'Homoiconic' instance with 'GHC.Generics.Generic', e.g.:
--
-- @
-- {-# LANGUAGE DeriveGeneric #-}
--
-- data MyData
--   = MyInt Int
--   | MyChar Char
--   deriving (Generic)
--
-- instance Homoiconic MyData
-- @
--
-- Sample snippet using above @MyData@:
--
-- >>> toCode (MyInt 42)
-- (MyInt 42)
-- >>> fromCode (toCode (MyChar 'a')) :: Maybe MyData
-- Just (MyChar 'a')
---
class Homoiconic a where
  -- | Convert Haskell value to 'Code'.
  toCode :: a -> Code
  {-# INLINE toCode #-}

  default toCode :: (Generic a, GToCode (Rep a)) => a -> Code
  toCode = genericToCode

  -- | Convert list of Haskell values to 'Code'.
  listToCode :: [a] -> Code
  listToCode xs =
     let xs' = map toCode xs
         l = getLoc (mkLocatedForm xs')
     in  LForm (L l (HsList xs'))
  {-# INLINE listToCode #-}

  -- | Convert 'Code' to Haskell value, or 'Failure' if the code could
  -- not be converted.
  parseCode :: Code -> Result a
  {-# INLINE parseCode #-}

  default parseCode :: (Generic a, GParseCode (Rep a)) => Code -> Result a
  parseCode = genericParseCode

  -- | Convert 'Code' to list of Haskell values, or 'Failure' if the code
  -- could not be converted.
  parseHsListCode :: Code -> Result [a]
  parseHsListCode xs =
    case unCode xs of
      HsList as -> mapM parseCode as
      _         -> fail "got non HsList value"
  {-# INLINE parseHsListCode #-}

-- | Like 'parseCode', but the result wrapped with 'Maybe' instead of 'Result'.
fromCode :: Homoiconic a => Code -> Maybe a
fromCode code = case parseCode code of
  Success a -> Just a
  _         -> Nothing

-- -------------------------------------------------------------------
--
-- Instances of Homoiconic
--
-- -------------------------------------------------------------------

--
-- Prelude
--

instance Homoiconic () where
  toCode _ = LForm (genSrc (Atom AUnit))
  parseCode a =
    case unCode a of
      Atom AUnit -> pure ()
      _          -> failedToParse "()"

instance Homoiconic Char where
  toCode = LForm . genSrc . Atom . AChar NoSourceText
  listToCode = LForm . genSrc . Atom . aString NoSourceText
  parseCode a =
    case unCode a of
      Atom (AChar _ x) -> pure x
      _                -> failedToParse "Char"
  parseHsListCode a = case unCode a of
                      Atom (AString _ s) -> pure (unpackFS s)
                      _                  -> failedToParse "String"

instance Homoiconic Int where
  toCode = integralToCode
  parseCode = integralFromCode

instance Homoiconic Word where
  toCode = integralToCode
  parseCode = integralFromCode

instance Homoiconic Integer where
  toCode = integralToCode
  parseCode = integralFromCode

instance Homoiconic Float where
  toCode = realFracToCode
  parseCode = fractionalFromCode

instance Homoiconic Double where
  toCode = realFracToCode
  parseCode = fractionalFromCode

instance Homoiconic a => Homoiconic [a] where
  toCode = listToCode
  parseCode = parseHsListCode

instance Homoiconic Bool where
  toCode = showAsSymbolCode
  parseCode a =
    case unCode a of
      Atom (ASymbol sym) | sym == "True"  -> pure True
                         | sym == "False" -> pure False
      _                                   -> failedToParse "Bool"

instance Homoiconic Ordering where
  toCode = showAsSymbolCode
  parseCode a =
    case unCode a of
      Atom (ASymbol sym) | sym == "EQ" -> pure EQ
                         | sym == "LT" -> pure LT
                         | sym == "GT" -> pure GT
      _                                -> failedToParse "Ordering"

instance Homoiconic a => Homoiconic (Maybe a) where
  toCode a =
    case a of
      Nothing -> toCode (aSymbol "Nothing")
      Just x  -> toCode1 "Just" x
  parseCode a =
    case unCode a of
      Atom (ASymbol "Nothing")                      -> pure Nothing
      List [LForm (L _ (Atom (ASymbol "Just"))), x] -> pure <$> parseCode x
      _                                             -> failedToParse "Maybe"

instance (Homoiconic a, Homoiconic b) => Homoiconic (Either a b) where
  toCode a =
    case a of
      Right x -> toCode1 "Right" x
      Left x  -> toCode1 "Left" x
  parseCode a =
    case unCode a of
      List [LForm (L _ (Atom (ASymbol x))), y]
        | x == "Right" -> fmap Right (parseCode y)
        | x == "Left"  -> fmap Left (parseCode y)
      _                -> failedToParse "Either"

instance (Homoiconic a, Homoiconic b) => Homoiconic (a, b) where
  toCode (a1, a2) = toCode2 "," a1 a2
  parseCode = parseCode2 "," (,)

instance (Homoiconic a, Homoiconic b, Homoiconic c)
         => Homoiconic (a, b, c) where
  toCode (a1, a2, a3) =
    toCode (List [symbolCode ",", toCode a1, toCode a2, toCode a3])
  parseCode a =
    case unCode a of
      List [LForm (L _ (Atom (ASymbol ","))), a1, a2, a3]
        -> (,,) <$> parseCode a1 <*> parseCode a2 <*> parseCode a3
      _ -> failedToParse "(,,)"

instance (Homoiconic a, Homoiconic b, Homoiconic c, Homoiconic d)
         => Homoiconic (a, b, c, d) where
  toCode (a1, a2, a3, a4) =
    toCode (List [ symbolCode ",", toCode a1, toCode a2, toCode a3
                 , toCode a4])
  parseCode a =
    case unCode a of
      List [LForm (L _ (Atom (ASymbol ","))), a1, a2, a3, a4]
        -> (,,,) <$>
           parseCode a1 <*> parseCode a2 <*> parseCode a3 <*> parseCode a4
      _ -> failedToParse "(,,,)"

instance (Homoiconic a, Homoiconic b, Homoiconic c, Homoiconic d, Homoiconic e)
         => Homoiconic (a, b, c, d, e) where
  toCode (a1, a2, a3, a4, a5) =
    toCode (List [ symbolCode ",", toCode a1, toCode a2, toCode a3
                 , toCode a4, toCode a5])
  parseCode a =
    case unCode a of
      List [LForm (L _ (Atom (ASymbol ","))), a1, a2, a3, a4, a5]
        -> (,,,,) <$>
           parseCode a1 <*> parseCode a2 <*> parseCode a3 <*>
           parseCode a4 <*> parseCode a5
      _ -> failedToParse "(,,,,)"

instance (Homoiconic a, Homoiconic b, Homoiconic c, Homoiconic d, Homoiconic e, Homoiconic f)
         => Homoiconic (a, b, c, d, e, f) where
  toCode (a1, a2, a3, a4, a5, a6) =
    toCode (List [ symbolCode ",", toCode a1, toCode a2, toCode a3
                 , toCode a4, toCode a5, toCode a6])
  parseCode a =
    case unCode a of
      List [LForm (L _ (Atom (ASymbol ","))), a1, a2, a3, a4, a5, a6]
        -> (,,,,,) <$>
           parseCode a1 <*> parseCode a2 <*> parseCode a3 <*>
           parseCode a4 <*> parseCode a5 <*> parseCode a6
      _ -> failedToParse "(,,,,,)"

--
-- Data.Complex
--

instance Homoiconic a => Homoiconic (Complex a) where
  toCode (a :+ b) = toCode2 ":+" a b
  parseCode = parseCode2 ":+" (:+)

--
-- Data.Fixed
--

instance Homoiconic (Fixed a) where
  toCode (MkFixed a) = toCode1 "MkFixed" a
  parseCode = parseCode1 "MkFixed" MkFixed

--
-- Data.Functor.Compose

instance Homoiconic (f (g a)) => Homoiconic (Compose f g a) where
  toCode (Compose a) = toCode1 "Compose" a
  parseCode = parseCode1 "Compose" Compose

--
-- Data.Functor.Const
--

instance Homoiconic a => Homoiconic (Const a b) where
  toCode (Const a) = toCode1 "Const" a
  parseCode = parseCode1 "Const" Const

--
-- Data.Functor.Identity
--

instance Homoiconic a=> Homoiconic (Identity a) where
  toCode (Identity a) = toCode1 "Identity" a
  parseCode = parseCode1 "Identity" Identity

--
-- Data.Functor.Product
--

instance (Homoiconic (f a), Homoiconic (g a))
         => Homoiconic (Product.Product f g a) where
  toCode (Product.Pair a b) = toCode2 "Pair" a b
  parseCode = parseCode2 "Pair" Product.Pair

--
-- Data.Functor.Sum
--

instance (Homoiconic (f a), Homoiconic (g a)) => Homoiconic (Sum.Sum f g a) where
  toCode a =
    case a of
      Sum.InL x -> toCode1 "InL" x
      Sum.InR x -> toCode1 "InR" x
  parseCode a =
    case unCode a of
      List [LForm (L _ (Atom (ASymbol tag))), b]
        | tag == "InL" -> Sum.InL <$> parseCode b
        | tag == "InR" -> Sum.InR <$> parseCode b
      _ -> failedToParse "Sum"

--
-- Data.Int
--

instance Homoiconic Int8 where
  toCode = integralToCode
  parseCode = integralFromCode

instance Homoiconic Int16 where
  toCode = integralToCode
  parseCode = integralFromCode

instance Homoiconic Int32 where
  toCode = integralToCode
  parseCode = integralFromCode

instance Homoiconic Int64 where
  toCode = integralToCode
  parseCode = integralFromCode

--
-- Data.List.NonEmpty
--

instance Homoiconic a => Homoiconic (NonEmpty a) where
  toCode (a :| as) = toCode2 ":|" a as
  parseCode = parseCode2 ":|" (:|)

--
-- Data.Monoid
--

instance Homoiconic All where
  toCode (All a) = toCode1 "All" a
  parseCode = parseCode1 "All" All

instance Homoiconic (f a) => Homoiconic (Alt f a) where
  toCode (Alt a) = toCode1 "Alt" a
  parseCode = parseCode1 "Alt" Alt

instance Homoiconic Any where
  toCode (Any a) = toCode1 "Any" a
  parseCode = parseCode1 "Any" Any

instance Homoiconic a => Homoiconic (Dual a) where
  toCode (Dual a) = toCode1 "Dual" a
  parseCode = parseCode1 "Dual" Dual

instance Homoiconic a => Homoiconic (First a) where
  toCode (First a) = toCode1 "First" a
  parseCode = parseCode1 "First" First

instance Homoiconic a => Homoiconic (Last a) where
  toCode (Last a) = toCode1 "Last" a
  parseCode = parseCode1 "Last" Last

instance Homoiconic a => Homoiconic (Product a) where
  toCode (Product a) = toCode1 "Product" a
  parseCode = parseCode1 "Product" Product

instance Homoiconic a => Homoiconic (Sum a) where
  toCode (Sum a) = toCode1 "Sum" a
  parseCode = parseCode1 "Sum" Sum

--
-- Data.Proxy
--

instance Homoiconic a => Homoiconic (Proxy a) where
  toCode _ = symbolCode "Proxy"
  parseCode a = case unCode a of
                 Atom (ASymbol "Proxy") -> pure Proxy
                 _                      -> failedToParse "Proxy"

--
-- Data.Version
--

instance Homoiconic Version where
  toCode (Version b t) = toCode2 "Version" b t
  parseCode = parseCode2 "Version" Version

--
-- Data.Ratio
--

instance (Integral a, Homoiconic a) => Homoiconic (Ratio a) where
  toCode a =
    let n = toCode (numerator a)
        d = toCode (denominator a)
    in toCode (List [symbolCode ":%", n, d])
  parseCode = parseCode2 ":%" (%)

--
-- Data.Semigroup
--

instance (Homoiconic a, Homoiconic b) => Homoiconic (Semigroup.Arg a b) where
  toCode (Semigroup.Arg a b) = toCode2 "Arg" a b
  parseCode = parseCode2 "Arg" Semigroup.Arg

instance Homoiconic a => Homoiconic (Semigroup.First a) where
  toCode (Semigroup.First a) = toCode1 "First" a
  parseCode = parseCode1 "First" Semigroup.First

instance Homoiconic a => Homoiconic (Semigroup.Last a) where
  toCode (Semigroup.Last a) = toCode1 "Last" a
  parseCode = parseCode1 "Last" Semigroup.Last

instance Homoiconic a => Homoiconic (Semigroup.Max a) where
  toCode (Semigroup.Max a) = toCode1 "Max" a
  parseCode = parseCode1 "Max" Semigroup.Max

instance Homoiconic a => Homoiconic (Semigroup.Min a) where
  toCode (Semigroup.Min a) = toCode1 "Min" a
  parseCode = parseCode1 "Min" Semigroup.Min

#if !MIN_VERSION_ghc(9,0,0)
instance Homoiconic a => Homoiconic (Semigroup.Option a) where
  toCode (Semigroup.Option a) = toCode1 "Option" a
  parseCode = parseCode1 "Option" Semigroup.Option
#endif

instance Homoiconic a => Homoiconic (Semigroup.WrappedMonoid a) where
  toCode (Semigroup.WrapMonoid a) = toCode1 "WrapMonoid" a
  parseCode = parseCode1 "WrapMonoid" Semigroup.WrapMonoid

--
-- Data.Word
--

instance Homoiconic Word8 where
  toCode = integralToCode
  parseCode = integralFromCode

instance Homoiconic Word16 where
  toCode = integralToCode
  parseCode = integralFromCode

instance Homoiconic Word32 where
  toCode = integralToCode
  parseCode = integralFromCode

instance Homoiconic Word64 where
  toCode = integralToCode
  parseCode = integralFromCode

--
-- Numeric.Natural
--

instance Homoiconic Natural where
  toCode = integralToCode
  parseCode = integralFromCode

--
-- Language.Finkel.Form
--

instance Homoiconic Atom where
  toCode = LForm . genSrc . Atom
  parseCode a =
    case unCode a of
      Atom x -> pure x
      _      -> failedToParse "Atom"

instance Homoiconic (Form Atom) where
  toCode = LForm . genSrc
  parseCode = pure . unCode

instance Homoiconic (LForm Atom) where
  toCode = id
  parseCode = pure

-- -------------------------------------------------------------------
--
-- Generic toCode
--
-- -------------------------------------------------------------------

-- | Generic variant of 'toCode'.
genericToCode :: (Generic a, GToCode (Rep a)) => a -> Code
genericToCode = unCodeArgs . gToCode . from
{-# INLINABLE genericToCode #-}

-- | To distinguish arguments of constructor from non-argument.
data CodeArgs
  = NonArg Code
  | Args [Code]

unCodeArgs :: CodeArgs -> Code
unCodeArgs ca = case ca of
  NonArg c -> c
  Args cs  -> toCode (List cs)
{-# INLINABLE unCodeArgs #-}

instance Semigroup.Semigroup CodeArgs where
  Args xs <> Args ys   = Args (xs Semigroup.<> ys)
  Args xs <> NonArg y  = Args (xs Semigroup.<> [y])
  NonArg x <> Args ys  = Args (x : ys)
  NonArg x <> NonArg y = Args [x, y]
  {-# INLINE (<>) #-}

-- | For making 'Code' with 'Generic' instances.
class GToCode f where
  gToCode :: f a -> CodeArgs

instance GToCode V1 where
  gToCode _ = NonArg undefined
  {-# INLINE gToCode #-}

instance GToCode U1 where
  gToCode U1 = NonArg nil
  {-# INLINE gToCode #-}

instance (GToCode f, GToCode g) => GToCode (f :+: g) where
  gToCode lr = case lr of
    L1 x -> gToCode x
    R1 x -> gToCode x
  {-# INLINE gToCode #-}

instance (GToCode f, GToCode g) => GToCode (f :*: g) where
  gToCode (f :*: g) = gToCode f Semigroup.<> gToCode g
  {-# INLINE gToCode #-}

instance Homoiconic c => GToCode (K1 i c) where
  gToCode (K1 x) = NonArg (toCode x)
  {-# INLINE gToCode #-}

instance GToCode f => GToCode (M1 D c f) where
  gToCode (M1 x) = gToCode x
  {-# INLINE gToCode #-}

instance (Constructor c, GToCode f) => GToCode (M1 C c f) where
  gToCode m1@(M1 x) =
    let constr = toCode (aSymbol (conName m1))
    in  case gToCode x of
      NonArg c -> if null c
        then NonArg constr
        else NonArg (toCode (List [constr, c]))
      Args cs -> NonArg (toCode (List (constr : cs)))
  {-# INLINE gToCode #-}

instance GToCode f => GToCode (M1 S c f) where
  gToCode (M1 x) = gToCode x
  {-# INLINE gToCode #-}

-- -------------------------------------------------------------------
--
-- Generic FromCode
--
-- -------------------------------------------------------------------

-- | Generic variant of 'fromCode'.
genericFromCode :: (Generic a, GParseCode (Rep a)) => Code -> Maybe a
genericFromCode x = case genericParseCode x of
  Success a -> Just a
  _         -> Nothing
{-# INLINABLE genericFromCode #-}

-- | Generic function to get result value from 'Code'.
genericParseCode :: (Generic a, GParseCode (Rep a)) => Code -> Result a
genericParseCode =
  let f a xs = if null xs
                 then pure (to a)
                 else fail "Unexpected leftover"
  in  runCodeP gParseCode fail f
{-# INLINABLE genericParseCode #-}

-- | For getting value from 'Code' with 'Generic' instances.
class GParseCode f where
  gParseCode :: CodeP (f a)

instance GParseCode V1 where
  gParseCode = pure undefined
  {-# INLINE gParseCode #-}

instance GParseCode U1 where
  gParseCode = pure U1
  {-# INLINE gParseCode #-}

instance (GParseCode f, GParseCode g) => GParseCode (f :+: g) where
  gParseCode = fmap L1 gParseCode <|> fmap R1 gParseCode
  {-# INLINE gParseCode #-}

instance (GParseCode f, GParseCode g) => GParseCode (f :*: g) where
  gParseCode = (:*:) <$> gParseCode <*> gParseCode
  {-# INLINE gParseCode #-}

instance Homoiconic c => GParseCode (K1 i c) where
  gParseCode =
    unconsP (\l c cs ->
              case parseCode c of
                Success a -> contP (K1 a) (LForm (L l (List cs)))
                _         -> failP ("Unexpected: " ++ show c))
  {-# INLINE gParseCode #-}

instance GParseCode f => GParseCode (M1 D c f) where
  gParseCode = fmap M1 gParseCode
  {-# INLINE gParseCode #-}

instance {-# OVERLAPPABLE #-} Constructor c => GParseCode (M1 C c U1) where
  gParseCode =
    let c1 :: M1 C c U1 a
        c1 = undefined
    in  eqP (toCode (aSymbol (conName c1))) *> fmap M1 gParseCode
  {-# INLINE gParseCode #-}

instance {-# OVERLAPPABLE #-} (Constructor c, GParseCode f)
  => GParseCode (M1 C c f) where
  gParseCode =
    let c1 :: M1 C c f a
        c1 = undefined
    in  eqCarP (toCode (aSymbol (conName c1))) *> fmap M1 gParseCode
  {-# INLINE gParseCode #-}

instance GParseCode f => GParseCode (M1 S c f) where
  gParseCode = fmap M1 gParseCode
  {-# INLINE gParseCode #-}

-- -------------------------------------------------------------------
--
-- Code parser for GParseCode
--
-- -------------------------------------------------------------------

-- | Dedicated data type to hold parsed result of 'Code'.
--
-- Using dedicated data type when parsing 'Code' data type for 'parseCode'. This
-- data type is intentionally not defined as an instance of 'Homoiconic', so
-- that the user defined data types can tell the parse error from explicit
-- failure constructor of the target type, e,g, 'Nothing' for 'Maybe', 'Left'
-- for 'Either', ... etc.
data Result a
  = Success a
  | Failure String
  deriving (Eq, Show)

instance Functor Result where
  fmap f r = case r of
    Success a -> Success (f a)
    Failure e -> Failure e
  {-# INLINE fmap #-}

instance Applicative Result where
  pure a = Success a
  {-# INLINE pure #-}
  f <*> m = f >>= flip fmap m
  {-# INLINE (<*>) #-}

instance Monad Result where
  m >>= k = case m of
    Success a -> k a
    Failure e -> Failure e
  {-# INLINE (>>=) #-}

instance MonadFail Result where
  fail = Failure
  {-# INLINE fail #-}

failedToParse :: String -> Result a
failedToParse ty = Failure ("Failed to parse " ++ ty)
{-# INLINABLE failedToParse #-}

-- | Simple parser for 'Code'.
newtype CodeP a =
  CodeP {runCodeP :: forall r. (String -> r) -- On failure
                  -> (a -> Code -> r)        -- On success
                  -> Code                    -- Input
                  -> r}

instance Functor CodeP where
  fmap f p = CodeP (\err go -> runCodeP p err (go . f))
  {-# INLINE fmap #-}

instance Applicative CodeP where
  pure a = CodeP (\_ go -> go a)
  {-# INLINE pure #-}

  f <*> p = f >>= flip fmap p
  {-# INLINE (<*>) #-}

instance Monad CodeP where
  m >>= k = CodeP (\err go -> runCodeP m err (\a -> runCodeP (k a) err go))
  {-# INLINE (>>=) #-}

instance Alternative CodeP where
  empty = failP "Alternative.empty"
  {-# INLINE empty #-}

  p1 <|> p2 =
    CodeP (\err go cs ->
            runCodeP p1 (\_ -> runCodeP p2 err go cs) go cs)
  {-# INLINE (<|>) #-}

failP :: String -> CodeP a
failP msg = CodeP (\err _ _ -> err msg)
{-# INLINABLE failP #-}

contP :: a -> Code -> CodeP a
contP a cs = CodeP (\_ go _ -> go a cs)
{-# INLINEABLE contP #-}

unconsP :: (SrcSpan -> Code -> [Code] -> CodeP a) -> CodeP a
unconsP f =
  CodeP (\err go cs ->
    case cs of
      LForm (L l (List (x : xs))) -> runCodeP (f l x xs) err go cs
      _                           -> err "Not a list")
{-# INLINEABLE unconsP #-}

eqP :: Code -> CodeP ()
eqP x =
  CodeP (\err go cs ->
    if cs == x
      then go () nil
      else err ("eqP: unexpected " ++ show cs))
{-# INLINABLE eqP #-}

eqCarP :: Code -> CodeP ()
eqCarP x =
  unconsP (\l c cs ->
             if x == c
               then contP () (LForm (L l (List cs)))
               else failP ("eqCarP: unexpected " ++ show c))
{-# INLINABLE eqCarP #-}

-- -------------------------------------------------------------------
--
-- Data to Code
--
-- -------------------------------------------------------------------

dataToCode :: Data d => d -> Code
dataToCode x =
  let constr = toConstr x
      isTupleStr cs = case cs of
                        '(':cs1 -> go cs1
                        _       -> False
                        where go xs = case xs of
                                        ',':xs' -> go xs'
                                        [')']   -> True
                                        _       -> False
      cstr = case showConstr constr of
               str | isTupleStr str -> ","
                   | otherwise  -> str
      hd = toCode (aSymbol cstr)
  in  case constrRep constr of
         IntConstr n   -> toCode (aIntegral n)
         FloatConstr f -> toCode (aFractional (fromRational f :: Double))
         CharConstr c  -> toCode c
         _             ->
           case gmapQ dataToCode x of
             [] -> hd
             _  -> toCode (List (hd:gmapQ dataToCode x))

-- -------------------------------------------------------------------
--
-- Auxiliary
--
-- -------------------------------------------------------------------

realFracToCode :: (Real a, Show a) => a -> Code
realFracToCode = LForm . genSrc . Atom . aFractional
{-# INLINABLE realFracToCode #-}

fractionalFromCode :: Fractional a => Code -> Result a
fractionalFromCode a =
  case unCode a of
    Atom (AFractional x) -> pure (fromRational (fl_value x))
    _                    -> failedToParse "fractional"
{-# INLINABLE fractionalFromCode #-}

symbolCode :: String -> Code
symbolCode = LForm . genSrc . Atom . aSymbol
{-# INLINABLE symbolCode #-}

showAsSymbolCode :: Show a => a -> Code
showAsSymbolCode = symbolCode . show
{-# INLINABLE showAsSymbolCode #-}

integralToCode :: Integral a => a -> Code
integralToCode = LForm . genSrc . Atom . aIntegral
{-# INLINABLE integralToCode #-}

integralFromCode :: Integral a => Code -> Result a
integralFromCode a =
  case unCode a of
    Atom (AInteger n) -> pure (fromIntegral (il_value n))
    _                 -> failedToParse "integral"

toCode1 :: Homoiconic a => FastString -> a -> Code
toCode1 tag arg1 =
  toCode (List [LForm (genSrc (Atom (ASymbol tag))), toCode arg1])
{-# INLINABLE toCode1 #-}

toCode2 :: (Homoiconic a, Homoiconic b) => FastString -> a -> b -> Code
toCode2 tag arg1 arg2 =
  toCode (List [ LForm (genSrc (Atom (ASymbol tag)))
               , toCode arg1, toCode arg2 ])
{-# INLINABLE toCode2 #-}

parseCode1 :: (Homoiconic a) => FastString -> (a -> h) -> Code -> Result h
parseCode1 tag f a =
  case unCode a of
    List [LForm (L _ (Atom (ASymbol tag'))), x]
      | tag == tag' -> f <$> parseCode x
    _               -> failedToParse (unpackFS tag)
{-# INLINABLE parseCode1 #-}

parseCode2 :: (Homoiconic a, Homoiconic b)
          => FastString -> (a -> b -> h) -> Code -> Result h
parseCode2 tag f a =
  case unCode a of
    List [LForm (L _ (Atom (ASymbol tag'))), x, y]
      | tag == tag' -> f <$> parseCode x <*> parseCode y
    _               -> failedToParse (unpackFS tag)
{-# INLINABLE parseCode2 #-}

max_stars_repo_name: finkel-lang/finkel

max_issues_repo_name: finkel-lang/finkel

max_forks_repo_name: finkel-lang/finkel


text:
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE DeriveDataTypeable #-}
{-# LANGUAGE DeriveGeneric #-}

-- HdpH implementaion of a Mandelbrot set creation
--
-- Based largely on code by Rob Stewart <R.Stewart@hw.ac.uk> who adapted it
-- from the monad-par implementation.
--
-- Author: Blair Archbald
-- Date: 2/3/2015
-----------------------------------------------------------------------------

module Main where

import Control.DeepSeq hiding (force)
import Control.Exception (evaluate)
import Control.Monad (when)
import Control.Parallel.HdpH
import Control.Parallel.HdpH.Strategies
import qualified Control.Monad.Par as MonadPar
import qualified Control.Monad.Par.Combinator as MonadPar.C
import qualified Control.Parallel.HdpH as HdpH (declareStatic)
import qualified Control.Parallel.HdpH.Strategies as Strategies (declareStatic)

import Data.Binary (Binary)
import Data.Complex
import Data.Functor ((<$>))
import Data.List (stripPrefix)
import Data.List.Split (splitOn)
import Data.Maybe (isJust,fromJust)
import Data.Monoid (mconcat)
import Data.Serialize
import Data.Typeable (Typeable)
import Data.Vector.Cereal () -- To get the generic deriving for cereal and vectors.
import qualified Data.Vector.Unboxed as V

import GHC.Generics (Generic)

import System.Clock
import System.Environment (getArgs)
import System.IO (stdout, stderr, hSetBuffering, BufferMode(..))
import System.Random (mkStdGen, setStdGen)

import Prelude

type MSTime = Double

data VecTree = Leaf (V.Vector Int)
             | MkNode VecTree VecTree
             deriving (Eq,Show,Generic,Typeable)
instance Serialize VecTree
instance NFData VecTree
instance ToClosure VecTree
  where locToClosure = $(here)

instance ToClosure (Int, Int)
  where locToClosure = $(here)

----------------------------------------------------------------------------
-- sequential mandel function

mandel :: Int -> Complex Double -> Int
mandel max_depth c = loop 0 0
  where
   loop i !z
    | i == max_depth       = i
    | magnitude z >= 2.0   = i
    | otherwise            = loop (i+1) (z*z + c)

checkSum :: VecTree -> Int
checkSum (Leaf vec) = V.foldl (+) 0 vec
checkSum (MkNode v1 v2) = checkSum v1 + checkSum v2

-----------------------------------------------
-- Mandel using monad-par
--  From: https://github.com/simonmar/monad-par/blob/master/examples/src/mandel.hs

monadparRunMandel :: Int -> Int -> Int -> Int -> MonadPar.Par VecTree
monadparRunMandel = monadparRunMandel' (-2) (-2) 2 2

--This is where the skeleton comes from.
monadparRunMandel' :: Double -> Double -- (minX, MinY)
                    -> Double -> Double -- (maxX, maxY)
                    -> Int -> Int -- (winX, winY)
                    -> Int -- Depth
                    -> Int -- Threshold
                    -> MonadPar.Par VecTree
monadparRunMandel' minX minY maxX maxY winX winY max_depth threshold =
  MonadPar.C.parMapReduceRangeThresh threshold (MonadPar.C.InclusiveRange 0 (winY-1))
     (\y ->
       do
          let vec = V.generate winX (\x -> mandelStep y x)
          seq (vec V.! 0) $ return (Leaf vec))
     (\ a b -> return$ MkNode a b)
     (Leaf V.empty)
  where
    mandelStep i j = mandel max_depth z
        where z = ((fromIntegral j * r_scale) / fromIntegral winY + minY) :+
                  ((fromIntegral i * c_scale) / fromIntegral winX + minX)
    r_scale  =  maxY - minY  :: Double
    c_scale =   maxX - minX  :: Double

-----------------------------------------------
-- Mandel using the generic HpdH D&C skeleton

hdphDandCMandel :: Int -> Int -> Int -> Int -> Par VecTree
hdphDandCMandel = hdphDandCMandel' (-2) (-2) 2 2

hdphDandCMandel' :: Double -> Double -- (minX, MinY)
                 -> Double -> Double -- (maxX, maxY)
                 -> Int -> Int -- (winX, winY)
                 -> Int -- Depth
                 -> Int -- Threshold
                 -> Par VecTree
hdphDandCMandel' minX minY maxX maxY winX winY maxDepth threshold = do
  res <- parDivideAndConquer
          $(mkClosure [|dc_trivial threshold|])
          $(mkClosure [|dc_decompose|])
          $(mkClosure [|dc_combine|])
          $(mkClosure [|dc_algorithm (minX, minY, maxX, maxY, winX, winY, maxDepth)|])
          (toClosure (0,winY-1))
  return $ unClosure res

dc_trivial :: Int -> Thunk (Closure (Int, Int) -> Bool)
dc_trivial threshold = Thunk $ \bnds -> let (min,max) = unClosure bnds
                                        in   max - min <= threshold

dc_decompose :: Closure (Int, Int) -> [Closure (Int, Int)]
dc_decompose bnds = let (min, max) = unClosure bnds
                        mid = min + (max - min) `quot` 2
                     in [toClosure (min, mid), toClosure (mid+1, max)]

-- This could prove to be a bottleneck and reason to use Rob's Skeleton
-- (TODO: compare performance between approaches)
dc_combine :: Closure a -> [Closure VecTree] -> Closure VecTree
dc_combine _ ts = toClosure $ foldl1 MkNode $ map unClosure ts

dc_algorithm :: (Double, Double, -- (minX, minY)
                 Double, Double, -- (maxX, maxY)
                 Int, Int, -- (winX, winY)
                 Int) -- maxDepth
              -> Thunk (Closure (Int,Int) -> Par (Closure VecTree))
dc_algorithm (minX, minY, maxX, maxY, winX, winY, maxDepth) = Thunk $ \bnds ->
  let (min,max) = unClosure bnds
      v = foldl go V.empty [min..max]
  in  do vec <- force v
         return $ toClosure (Leaf vec)
  where
    go a y = a V.++ V.generate winX (\x -> mandelStep y x)
    mandelStep i j = mandel maxDepth (calcZ i j)
    calcZ i j = ((fromIntegral j * r_scale) / fromIntegral winY + minY) :+
                ((fromIntegral i * c_scale) / fromIntegral winX + minX)
    r_scale =  maxY - minY  :: Double
    c_scale =  maxX - minX  :: Double

---------------------------------------------------------------------------
-- Timed Version
hdphDandCMandelTimed :: Int -> Int -> Int -> Int -> Par VecTree
hdphDandCMandelTimed = hdphDandCMandelTimed' (-2) (-2) 2 2

hdphDandCMandelTimed' :: Double -> Double -- (minX, MinY)
                 -> Double -> Double -- (maxX, maxY)
                 -> Int -> Int -- (winX, winY)
                 -> Int -- Depth
                 -> Int -- Threshold
                 -> Par VecTree
hdphDandCMandelTimed' minX minY maxX maxY winX winY maxDepth threshold = do
  res <- parDivideAndConquer
          $(mkClosure [|dc_trivial threshold|])
          $(mkClosure [|dc_decompose|])
          $(mkClosure [|dc_combine|])
          $(mkClosure [|dc_algorithm_timed (minX, minY, maxX, maxY, winX, winY, maxDepth)|])
          (toClosure (0,winY-1))
  return $ unClosure res

dc_algorithm_timed :: (Double, Double, -- (minX, minY)
                       Double, Double, -- (maxX, maxY)
                       Int, Int, -- (winX, winY)
                       Int) -- maxDepth
                   -> Thunk (Closure (Int,Int) -> Par (Closure VecTree))
dc_algorithm_timed (minX, minY, maxX, maxY, winX, winY, maxDepth) = Thunk $ \bnds ->
  let (min,max) = unClosure bnds
      v = foldl go V.empty [min..max]
  in  do s <- io $ getTime Monotonic
         vec <- force v
         e <- io $ getTime Monotonic
         io $ putStrLn $ show min ++ "," ++ show max ++ ":" ++ show (timeDiffMSecs s e)
         return $ toClosure (Leaf vec)
  where
    go a y = a V.++ V.generate winX (\x -> mandelStep y x)
    mandelStep i j = mandel maxDepth (calcZ i j)
    calcZ i j = ((fromIntegral j * r_scale) / fromIntegral winY + minY) :+
                ((fromIntegral i * c_scale) / fromIntegral winX + minX)
    r_scale =  maxY - minY  :: Double
    c_scale =  maxX - minX  :: Double
-----------------------------------------------------------------------------
-- initialisation, argument processing and 'main'

timeDiffMSecs :: TimeSpec -> TimeSpec -> MSTime
timeDiffMSecs (TimeSpec s1 n1) (TimeSpec s2 n2) = fromIntegral (t2 - t1)
                                                          /
                                                  fromIntegral (10 ^ 6)
  where t1 = (fromIntegral s1 * 10 ^ 9) + fromIntegral n1
        t2 = (fromIntegral s2 * 10 ^ 9) + fromIntegral n2

timeIOMs :: IO a -> IO (a, MSTime)
timeIOMs action = do s  <- getTime Monotonic
                     x  <- action
                     e  <- getTime Monotonic
                     return (x, timeDiffMSecs s e)

-----------------------------------------------------
-- Static Closures

$(return []) -- Bring the types into scope so that reify works.
declareStatic :: StaticDecl
declareStatic =
  mconcat
    [HdpH.declareStatic,         -- declare Static deserialisers
     Strategies.declareStatic,   -- from imported modules
     declare (staticToClosure :: StaticToClosure VecTree),
     declare (staticToClosure :: StaticToClosure (Int,Int)),

     declare $(static 'dc_trivial),
     declare $(static 'dc_decompose),
     declare $(static 'dc_combine),
     declare $(static 'dc_algorithm),

     declare $(static 'dc_algorithm_timed),

     declare $(static 'mandel)]
---------------------------------------------------------------------------

-- parse (optional) arguments in this order:
-- * version to run
-- * X value for Mandel
-- * Y value for Mandel
-- * Depth value for Mandel
-- * Threshold for Mandel
parseArgs :: [String] -> (Int, Int, Int, Int, Int, Int)
parseArgs []     = (defVers, defX, defY, defDepth,defThreshold,defExpected)
parseArgs (s:ss) =
  let go :: Int -> [String] -> (Int, Int, Int, Int,Int, Int)
      go v []              = (v, defX, defY, defDepth,defThreshold,defExpected)
      go v [s1]            = (v, defX, read s1,  defDepth,defThreshold,defExpected)
      go v [s1,s2]         = (v, read s1,  read s2,  defDepth,defThreshold,defExpected)
      go v [s1,s2,s3]      = (v, read s1,  read s2,  read s3, defThreshold,defExpected)
      go v [s1,s2,s3,s4] = (v, read s1,  read s2,  read s3, read s4,defExpected)
      go v (s1:s2:s3:s4:s5:_) = (v, read s1,  read s2,  read s3, read s4, read s5)
  in case stripPrefix "v" s of
       Just s' -> go (read s') ss
       Nothing -> go defVers (s:ss)

-- defaults from Simon Marlow from monad-par example
defVers, defX, defY, defDepth, defThreshold, defExpected :: Int
defVers      = 0
defX         = 1024
defY         = 1024
defDepth     = 256
defThreshold = 1
defExpected  = 0

main :: IO ()
main = do
  hSetBuffering stdout LineBuffering
  hSetBuffering stderr LineBuffering
  register Main.declareStatic

  hdphCfg    <- flip updateConf defaultRTSConf =<< getArgs
  (conf, as) <- case hdphCfg of
    Left err -> error $ "Could not initialise HdpH config: " ++ show err
    Right a  -> return a

  let (version, valX, valY, valDepth,valThreshold,expected) = parseArgs as
  case version of
      1 -> do (p, t) <- timeIOMs $ evaluate =<< return (MonadPar.runPar
                (monadparRunMandel valX valY valDepth valThreshold))
              printOutput (Just p,t)
      2 -> do res <- timeIOMs $ evaluate =<< runParIO conf
                                (hdphDandCMandel valX valY valDepth valThreshold)
              printOutput res
      3 -> do res <- timeIOMs $ evaluate =<< runParIO conf
                                (hdphDandCMandelTimed valX valY valDepth valThreshold)
              printOutput res
      _ -> return ()

printOutput :: (Maybe VecTree, MSTime) -> IO ()
printOutput (p,t) = case p of
                     Just pixels -> do
                       putStrLn $ "CHECKSUM: " ++ show (checkSum pixels)
                       putStrLn $ "RUNTIME: "  ++ show t ++ "ms"
                     Nothing -> return ()

max_stars_repo_name: BlairArchibald/HdpH-applications

max_issues_repo_name: BlairArchibald/HdpH-applications

max_forks_repo_name: BlairArchibald/HdpH-applications


text:
{-| This module defines the 'Produce' typeclass, used for generating random
    values for testing in StrictCheck.

    'Produce' is a strict generalization of "Test.QuickCheck"'s 'Arbitrary'
    typeclass. Paired with 'Consume' (a generalization of 'CoArbitrary') it can
    be used to create random non-strict functions, whose strictness behavior is
    dependent on the values given to them.
-}

module Test.StrictCheck.Produce
  ( Produce(..)
  -- * Tools for writing 'Produce' instances
  , recur
  , build
  -- * Producing non-strict functions
  , returning
  , variadic
  -- * Integration with "Test.QuickCheck"'s @Arbitrary@
  , Lazy(..)
  , freely
  -- * Abstract types representing input to a function
  , Input
  , Inputs
  -- * The traversal distribution for processing @Input@s
  , draws
  ) where

import Test.QuickCheck hiding (variant)
import Test.QuickCheck.Gen.Unsafe

import Test.StrictCheck.Internal.Inputs
import Test.StrictCheck.Consume
import Test.StrictCheck.Curry

import Generics.SOP
import Data.Complex
import Data.Monoid ((<>))

-------------------------------------------------------
-- The user interface for creating Produce instances --
-------------------------------------------------------

-- TODO: parameterize over destruction pattern?

-- | Produce an arbitrary value of type @b@, such that destructing that value
-- incrementally evaluates some input to a function.
--
-- Writing instances of @Produce@ is very similar to writing instances of
-- QuickCheck's 'Arbitrary'. The distinction: when making a recursive call to
-- produce a subfield of a structure, __always__ use 'build' or 'recur', and
-- __never__ a direct call to 'produce' itself. This ensures that the input can
-- potentially be demanded at any step of evaluation of the produced value.
--
-- If, in the course of generating a value of type @b@, you need to generate a
-- random value of some other type, which is /not/ going to be a subpart of the
-- resultant @b@ (e.g. a length or depth), use a direct call to @arbitrary@ or
-- some other generator which does not consume input.
--
-- An example instance of @Produce@:
--
-- > data D a
-- >   = X a
-- >   | Y [Int]
-- >
-- > instance Produce a => Produce (D a) where
-- >   produce =
-- >     oneof [ fmap X recur
-- >           , fmap Y recur
-- >           ]
class Produce b where
  produce :: (?inputs::Inputs) => Gen b

theInputs :: (?inputs::Inputs) => [Input]
theInputs = destruct ?inputs

-- | Given an input-consuming producer, wrap it in an outer layer of input
-- consumption, so that this consumption can be interleaved when the producer is
-- called recursively to generate a subfield of a larger produced datatype.
build :: (?inputs::Inputs) => ((?inputs::Inputs) => Gen a) -> Gen a
build gen = do
  (v, is') <- draws theInputs
  vary v $ let ?inputs = Inputs is' in gen

-- | Destruct some inputs to generate an output. This function handles the
-- interleaving of input destruction with output construction. When producing a
-- data type, it should be called to produce each subfield -- *not* produce
-- itself.
recur :: (Produce a, ?inputs::Inputs) => Gen a
recur = build produce

---------------------------------------
-- How to make random lazy functions --
---------------------------------------

-- NOTE: This instance must be defined in this module, as it has to break the
-- abstraction of the Inputs type. No other instance needs to break this.
-- Incidentally, it also must break Gen's abstraction barrier, because it needs
-- to use promote to make a function.

instance (Consume a, Produce b) => Produce (a -> b) where
  produce = returning produce

-- | Create an input-consuming producer of input-consuming functions, given an
-- input-consuming producer for results of that function.
returning
  :: (Consume a, ?inputs::Inputs)
  => ((?inputs::Inputs) => Gen b)
  -> Gen (a -> b)
returning out =
  promote $ \a ->
    let ?inputs = Inputs (consume a : theInputs)
    in build out

-- | Create an input-consuming producer of input-consuming functions, of any
-- arity. This will usually be used in conjuntion with type application, to
-- specify the type(s) of the argument(s) to the function.
variadic ::
  forall args result.
  (All Consume args, Curry args, ?inputs::Inputs)
  => ((?inputs::Inputs) => Gen result)
  -> Gen (args ⋯-> result)
variadic out =
  fmap (curryAll @args @_ @(NP I)) . promote $ \args ->
    let ?inputs =
          Inputs . (++ theInputs) $
            hcollapse $ hcliftA (Proxy @Consume) (K . consume . unI) args
    in build out

-------------------------------------------------------------------------
-- Random destruction of the original input, as transformed into Input --
-------------------------------------------------------------------------

-- | Destruct a random subpart of the given 'Input's, returning the 'Variant'
-- corresponding to the combined information harvested during this process, and
-- the remaining "leaves" of the inputs yet to be destructed
--
-- To maximize the likelihood that different random consumption paths through
-- the same value will diverge (desirable when generating functions with
-- interesting strictness), @draws@ destructs the forest of @Input@s as a
-- depth-first random traversal with a budget sampled from a geometric
-- distribution with expectation 1.
draws :: [Input] -> Gen (Variant, [Input])
draws inputs = go [inputs]
  where
    -- Mutually recursive:
    go, inwardFrom :: [[Input]] -> Gen (Variant, [Input])

    go levels =
      oneof                               -- 50% choice between:
        [ return (mempty, concat levels)  -- stop consuming input, or
        , inwardFrom levels ]             -- keep consuming input

    inwardFrom levels =
      case levels of
        [            ] -> return mempty         -- if no more input: stop
        [  ] : outside -> inwardFrom outside    -- if nothing here: backtrack
        here : outside -> do                    -- if something here: go deeper
          (Input v inside, here') <- pick here
          vary v $ do
            (entropy, levels') <- go (inside : here' : outside)  -- back to 'go'
            return (v <> entropy, levels')

    -- Pick a random list element and the remaining list
    pick :: [a] -> Gen (a, [a])
    pick as = do
      index <- choose (0, length as - 1)
      let (before, picked : after) = splitAt index as
      return (picked, before ++ after)

---------------------------------------------
-- Integration with QuickCheck's Arbitrary --
---------------------------------------------

-- | We hook into QuickCheck's existing Arbitrary infrastructure by using
-- a newtype to differentiate our special way of generating things.
newtype Lazy a
  = Lazy { runLazy :: a }

instance Produce a => Arbitrary (Lazy a) where
  arbitrary = Lazy <$> freely produce

-- | Actually produce an output, given an input-consuming producer. If a
-- function is to be produced, it will be almost-certainly non-strict.
freely :: ((?inputs::Inputs) => Gen a) -> Gen a
freely p = let ?inputs = Inputs [] in p

---------------
-- Instances --
---------------

instance Produce ()       where produce = arbitrary
instance Produce Bool     where produce = arbitrary
instance Produce Ordering where produce = arbitrary

instance Produce Char     where produce = arbitrary
instance Produce Word     where produce = arbitrary
instance Produce Int      where produce = arbitrary
instance Produce Double   where produce = arbitrary
instance Produce Float    where produce = arbitrary
instance Produce Rational where produce = arbitrary
instance Produce Integer  where produce = arbitrary

instance (Arbitrary a, RealFloat a) => Produce (Complex a) where
  produce = arbitrary

instance Produce a => Produce (Maybe a) where
  produce =
    oneof [ return Nothing
          , Just <$> recur
          ]

instance (Produce a, Produce b) => Produce (Either a b) where
  produce =
    oneof [ Left <$> recur
          , Right <$> recur
          ]

instance (Produce a) => Produce [a] where
  produce =
    frequency [ (1, return [])
              , (1, (:) <$> recur
                        <*> recur)
              ]

max_stars_repo_name: janpath/StrictCheck

max_issues_repo_name: janpath/StrictCheck

max_forks_repo_name: janpath/StrictCheck


text:
module Observable.MCMC.MALA where

import Control.Monad
import Control.Monad.Primitive
import Control.Monad.State.Strict
import Observable.Core
import Observable.Types
import Data.HashMap.Strict (HashMap)
import qualified Data.HashMap.Strict as HashMap
import Data.Vector.Unboxed (Vector)
import qualified Data.Vector.Unboxed as V
import Statistics.Distribution
import Statistics.Distribution.Normal

sphereGauss :: Vector Double -> Vector Double -> Double -> Double
sphereGauss xs mu sd = product $ zipWith density normalDists xsAsList where
  xsAsList = V.toList xs
  muAsList = V.toList mu
  normalDists = map (`normalDistr` sd) muAsList

(.*) :: Double -> Vector Double -> Vector Double
z .* xs = V.map (* z) xs

(.+) :: Vector Double -> Vector Double -> Vector Double
xs .+ ys = V.zipWith (+) xs ys

localMean :: Target Double -> Vector Double -> Double -> Vector Double
localMean target position e = position .+ scaledGradient position where
  grad = handleGradient (gradient target)
  scaledGradient p = (0.5 * e * e) .* grad p

perturb
  :: PrimMonad m
  => Target Double
  -> Vector Double
  -> Double
  -> Observable m (Vector Double)
perturb target position e = do
  zs <- V.replicateM (V.length position) standardNormal
  return $ localMean target position e .+ (e .* zs)

mala e = do
  Chain current target _ store <- get
  let step = getStepSize e store
  proposal <- lift $ perturb target current step
  zc       <- lift unit

  let cMean    = localMean target current step
      pMean    = localMean target proposal step
      next     = nextState target (current, cMean) (proposal, pMean) step zc
      newStore = updateStepSize step store

  put $ Chain next target (logObjective target next) newStore
  return next

getStepSize :: Maybe Double -> OptionalStore -> Double
getStepSize (Just step) _ = step
getStepSize Nothing store = step where
  (ODouble step) = HashMap.lookupDefault (ODouble 1.0) MALA store

updateStepSize :: Double -> OptionalStore -> OptionalStore
updateStepSize step = HashMap.insert MALA (ODouble step)

nextState
  :: Target Double
  -> (Vector Double, Vector Double)
  -> (Vector Double, Vector Double)
  -> Double
  -> Double
  -> Vector Double
nextState target (current, cMean) (proposal, pMean) e z
    | z < acceptProb = proposal
    | otherwise      = current
  where
    ratio = acceptRatio target (current, cMean) (proposal, pMean) e
    acceptProb | isNaN ratio = 0
               | otherwise   = ratio

acceptRatio
  :: Target Double
  -> (Vector Double, Vector Double)
  -> (Vector Double, Vector Double)
  -> Double
  -> Double
acceptRatio target (current, cMean) (proposal, pMean) e = exp . min 0 $
    logObjective target proposal + log (sphereGauss current pMean e)
  - logObjective target current  - log (sphereGauss proposal cMean e)

max_stars_repo_name: jtobin/deprecated-observable

max_issues_repo_name: jtobin/deprecated-observable

max_forks_repo_name: jtobin/deprecated-observable


text:
{-|
Module: MachineLearning.NeuralNetwork.TanhActivation.
Description: Tanh Activation.
Copyright: (c) Alexander Ignatyev, 2017
License: BSD-3
Stability: experimental
Portability: POSIX

Tanh Activation.
-}

module MachineLearning.NeuralNetwork.TanhActivation
(
  tanh
  , gradient
)

where

import qualified Numeric.LinearAlgebra as LA
import MachineLearning.Types (Matrix)

tanhGradient :: Matrix -> Matrix
tanhGradient x = 1 - tanhx*tanhx
  where tanhx = tanh x

gradient :: Matrix -> Matrix -> Matrix
gradient x dx = dx * (tanhGradient x)

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
{-# LANGUAGE ScopedTypeVariables, TypeFamilies #-}

module Fuml.Supervised where

import qualified Fuml.Base.LinearRegression as LR
import Fuml.Core
import Numeric.LinearAlgebra
import Control.Monad.Identity

newtype Beta = Beta { unBeta :: Vector Double }

ols :: Supervisor Identity Double Beta Double
ols = Supervisor $ \_ theData ->
              let beta = LR.ols theData
              in return $ linPredict beta

type Weight = Double

olsWeighted :: Supervisor Identity (Double, Weight) Beta Double
olsWeighted
   = Supervisor $ \_ theDataW ->
              let wvs = fromList $ map (\(_,(_,w)) -> w) theDataW
                  theData = map (\(xs,(y,_)) -> (xs,y)) theDataW
                  beta = LR.wols wvs theData
              in return $ linPredict beta

ridge :: Double -> Supervisor Identity Double Beta Double
ridge gammac = Supervisor $ \_ xys ->
  let npars = size $ fst $ head xys
      gamma = diag $ konst gammac npars
      beta = LR.ridge gamma xys
  in return $ linPredict beta

linPredict :: Vector Double -> Predict Beta Double
linPredict beta = Predict (Beta beta) (`dot` beta)

max_stars_repo_name: ekalosak/open

max_issues_repo_name: ekalosak/open

max_forks_repo_name: ekalosak/open


text:
{-# LANGUAGE RecordWildCards #-}
-- Chi square tests for random generators
module MWC.ChiSquare (
  tests
  ) where

import Control.Applicative
import Control.Monad

import Data.Typeable
import Data.Word
import Data.List (find)
import qualified Data.Vector.Unboxed              as U
import qualified Data.Vector.Unboxed.Mutable      as M
import qualified System.Random.MWC                as MWC
import qualified System.Random.MWC.Distributions  as MWC
import qualified System.Random.MWC.CondensedTable as MWC

import Statistics.Types
import Statistics.Test.ChiSquared
import Statistics.Distribution
import Statistics.Distribution.Poisson
import Statistics.Distribution.Binomial
import Statistics.Distribution.Geometric

import Test.HUnit hiding (Test)
import Test.Framework
import Test.Framework.Providers.HUnit

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

tests :: MWC.GenIO -> Test.Framework.Test
tests g = testGroup "Chi squared tests"
    -- Word8 tests
  [ uniformRTest (0,255 :: Word8) g
  , uniformRTest (0,254 :: Word8) g
  , uniformRTest (0,129 :: Word8) g
  , uniformRTest (0,126 :: Word8) g
  , uniformRTest (0,10  :: Word8) g
    -- * Tables
  , ctableTest   [1] g
  , ctableTest   [0.5,  0.5] g
  , ctableTest   [0.25, 0.25, 0.25, 0.25] g
  , ctableTest   [0.25, 0.5,  0.25] g
  , ctableTest   [1/3 , 1/3, 1/3] g
  , ctableTest   [0.1,  0.9] g
  , ctableTest   (replicate 10 0.1) g
    -- ** Poisson
  , poissonTest 0.2  g
  , poissonTest 1.32 g
  , poissonTest 6.8  g
  , poissonTest 100  g
    -- ** Binomial
  , binomialTest 4   0.5 g
  , binomialTest 10  0.1 g
  , binomialTest 10  0.6 g
  , binomialTest 10  0.8 g
  , binomialTest 100 0.3 g
    -- ** Geometric
  , geometricTest 0.1 g
  , geometricTest 0.5 g
  , geometricTest 0.9 g
  ]

----------------------------------------------------------------
-- | RNG and corresonding distribution
data Generator = Generator {
    generator    :: MWC.GenIO -> IO Int
  , probabilites :: U.Vector Double
  }

-- | Apply chi square test for a distribution
sampleTest :: String              -- ^ Name of test
           -> Generator           -- ^ Generator to test
           -> Int                 -- ^ N of events
           -> MWC.GenIO           -- ^ PRNG state
           -> Test.Framework.Test
sampleTest nm (Generator{..}) n g = testCase nm $ do
  let size = U.length $ probabilites
  h <- histogram (generator g) size n
  let w      = U.map (* fromIntegral n) probabilites
      Just t = chi2test 0 $ U.zip h w
  case isSignificant (mkPValue 0.01) t of
    Significant    -> assertFailure ("Significant: " ++ show t)
    NotSignificant -> return ()
{-# INLINE sampleTest #-}

-- | Fill histogram using supplied generator
histogram :: IO Int             -- ^ Random generator
          -> Int                -- ^ N of outcomes
          -> Int                -- ^ N of events
          -> IO (U.Vector Int)
histogram gen size n = do
  arr <- M.replicate size 0
  replicateM_ n $ do i <- gen
                     when (i < size) $ M.write arr i . (+1) =<< M.read arr i
  U.unsafeFreeze arr
{-# INLINE histogram #-}

-- | Test uniformR
uniformRTest :: (MWC.Variate a, Typeable a, Show a, Integral a) => (a,a) -> MWC.GenIO -> Test.Framework.Test
uniformRTest (a,b)
  = sampleTest ("uniformR: " ++ show (a,b) ++ " :: " ++ show (typeOf a)) gen (10^5)
  where
    n   = fromIntegral b - fromIntegral a + 1
    gen = Generator { generator    = \g -> fromIntegral . subtract a <$> MWC.uniformR (a,b) g
                    , probabilites = U.replicate n (1 / fromIntegral n)
                    }
{-# INLINE uniformRTest #-}

-- | Test for condensed tables
ctableTest :: [Double] -> MWC.GenIO -> Test.Framework.Test
ctableTest ps
  = sampleTest ("condensedTable: " ++ show ps) gen (10^4)
  where
    gen = Generator
          { generator    = MWC.genFromTable $ MWC.tableFromProbabilities $ U.fromList $ zip [0..] ps
          , probabilites = U.fromList ps
          }

-- | Test for condensed table for poissson distribution
poissonTest :: Double -> MWC.GenIO -> Test.Framework.Test
poissonTest lam
  = sampleTest ("poissonTest: " ++ show lam) gen (10^4)
  where
    pois      = poisson lam
    Just nMax = find (\n -> probability pois n < 2**(-33)) [floor lam ..]
    gen = Generator
          { generator    = MWC.genFromTable (MWC.tablePoisson lam)
          , probabilites = U.generate nMax (probability pois)
          }

-- | Test for condensed table for binomial distribution
binomialTest :: Int -> Double -> MWC.GenIO -> Test.Framework.Test
binomialTest n p
  = sampleTest ("binomialTest: " ++ show p ++ " " ++ show n) gen (10^4)
  where
    binom = binomial n p
    gen   = Generator
            { generator    = MWC.genFromTable (MWC.tableBinomial n p)
            , probabilites = U.generate (n+1) (probability binom)
            }

-- | Test for geometric distribution
geometricTest :: Double -> MWC.GenIO -> Test.Framework.Test
geometricTest gd
  = sampleTest ("geometricTest: " ++ show gd) gen (10^4)
  where
    n   = 1000
    gen = Generator
          { generator    = MWC.geometric1 gd
          , probabilites = U.generate (n+1) (probability $ geometric gd)
          }

max_stars_repo_name: bos/mwc-random

max_issues_repo_name: bos/mwc-random

max_forks_repo_name: bos/mwc-random


text:
-- ex 1.13 (1.2, 1.17)

import qualified Data.Vector as V

import Statistics.Sample.Histogram (histogram)
import Statistics.Sample as S (mean, stdDev, range)

import Statistics.Distribution
import Statistics.Distribution.Normal
import Statistics.Distribution.Binomial

-- Average wind speeds for 45 us cities:
samples :: V.Vector Double
samples = V.fromList
  [ 8.9,12.4, 8.6,11.3, 9.2, 8.8,35.1, 6.2, 7.0
  , 7.1,11.8,10.7, 7.6, 9.1, 9.2, 8.2, 9.0, 8.7
  , 9.1,10.9,10.3, 9.6, 7.8,11.5, 9.3, 7.9, 8.8
  , 8.8,12.7, 8.4, 7.8, 5.7,10.5,10.5, 9.6, 8.9
  ,10.2,10.3, 7.7,10.6, 8.3, 8.8, 9.5, 8.8, 9.4
  ]

-- 1.2 histogram
hist :: (V.Vector Double, V.Vector Int)
hist = histogram 12 samples

main :: IO ()
main = do
  let (lowerBs, hs) = hist
  putStrLn $ "The lower bounds: " ++ show lowerBs
  putStrLn $ "The histogram: " ++ show hs

  let av = S.mean samples
      sd = S.stdDev samples
  putStrLn $ "The average value is: " ++ show av ++ " and the standard deviation is: " ++ show sd

  let f k =
       V.length . V.filter (\x -> (av-k*sd)<x && x< (av+k*sd)) $ samples
  putStrLn $ "Our sample size is: " ++ show (V.length samples)
  putStrLn $ "Within one sigma:   " ++ show (f 1)
  putStrLn $ "Within two sigma:   " ++ show (f 2)
  putStrLn $ "Within three sigma: " ++ show (f 3)

  putStrLn $ "range/4 is:           " ++ show (S.range samples / 4)
  putStrLn $ "The average value is: " ++ show av

max_stars_repo_name: raydsameshima/Stat

max_issues_repo_name: raydsameshima/Stat

max_forks_repo_name: raydsameshima/Stat


text:
-- |
-- Module      : Jeopardy.Graphics
-- Description : Rendering functionality for the Jeopardy game
-- Copyright   : (c) Jonatan H Sundqvist, 2015
-- License     : MIT
-- Maintainer  : Jonatan H Sundqvist
-- Stability   : experimental|stable
-- Portability : POSIX (not sure)
--

-- Created July 7 2015

-- TODO | - 'App' type with configurations, event listeners, etc.
--        -

-- SPEC | -
--        -

module Jeopardy.Graphics where

---------------------------------------------------------------------------------------------------
-- We'll need these
---------------------------------------------------------------------------------------------------
import qualified Graphics.Rendering.Cairo as Cairo --

import Control.Monad (when, forM_) --
import Data.Complex

import qualified Southpaw.Picasso.Palette as Palette --
import Southpaw.Picasso.RenderUtils                  --

import Jeopardy.Core --

---------------------------------------------------------------------------------------------------
-- Types
---------------------------------------------------------------------------------------------------
type Point = Complex Double

---------------------------------------------------------------------------------------------------
-- Data
---------------------------------------------------------------------------------------------------
π = pi
τ = 2*π

---------------------------------------------------------------------------------------------------
-- Functions
---------------------------------------------------------------------------------------------------
-- Game specific rendering ------------------------------------------------------------------------
-- |
renderGame :: Game -> Cairo.Render ()
renderGame game = do
	--
	let (dx, padx) = (42, 5)
	mapM_ (\ (n, cat) -> renderCategory ((50+n*(dx+padx)*2):+50) (dx:+dx) 4 cat) . zip [1..] $ _board game

-- |
renderAnswerTile :: Point -> Double -> String -> Cairo.Render ()
renderAnswerTile (cx:+cy) radius answer = do
	Cairo.arc cx cy radius 0 τ
	Palette.choose Palette.blue
	Cairo.fill

	Cairo.arc cx cy (radius*0.86) 0 τ
	Palette.choose Palette.cadetblue
	Cairo.fill

	Palette.choose Palette.black
	Cairo.setFontSize 18
	renderCentredText (cx:+cy) answer

-- |
renderCategoryHeading :: Point -> Point -> String -> Cairo.Render ()
renderCategoryHeading (cx:+cy) (dx:+dy) title = do
	Cairo.rectangle (cx-dx) (cy-dy) (dx*2) (dy*2)
	Palette.choose Palette.blue
	Cairo.fill

	Cairo.rectangle (cx-(dx*0.86)) (cy-(dy*0.86)) (dx*2*0.86) (dy*2*0.86)
	Palette.choose Palette.cadetblue
	Cairo.fill

	Palette.choose Palette.black
	Cairo.setFontSize 18
	renderCentredText (cx:+cy) title

-- |
-- TODO: Add layout option parameters
renderCategory :: Point -> Point -> Double -> Category -> Cairo.Render ()
renderCategory (cx:+cy) (dx:+dy) pady (Category { _questions=q, _title=t }) = do
	renderCategoryHeading (cx:+cy) (dx:+(dy*0.8)) t
	mapM_ (\ (n, q) -> renderAnswerTile (cx:+(cy+n*(dy+pady)*2)) dx . show $ _value q) . zip [1..] $ q

-- General rendering utilities --------------------------------------------------------------------
-- |
renderCircle :: Double -> Point -> Cairo.Render ()
renderCircle radius (cx:+cy) = Cairo.arc cx cy radius 0 τ

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

max_stars_repo_name: SwiftsNamesake/Leopardy

max_issues_repo_name: SwiftsNamesake/Leopardy

max_forks_repo_name: SwiftsNamesake/Leopardy


text:
module Main where

-- third party imports
import Criterion.Main
import Data.Complex
import Data.Tuple
import Graphics.Image
import System.Environment

-- local imports
import CommonFunctions
import JuliaSet

main :: IO ()
main = do
  args <- getArgs
  if length args == 0 then do
    usage
  else do
    if length args < 5 then do
      visualizeJuliaBench
    else do
      let outputFilename = args !! 4
      -- putStrLn "visualizing julia set..."
      runWithCliArgs args
      -- putStrLn ("Done! Image output to " ++ outputFilename)

visualizeJuliaBench = defaultMain [
    bgroup "visualize julia set" [
      -- bench "f1 ((-0.4) :+ 0.65) VU" $ whnf makeImageVU g,
      -- bench "f1 ((-0.4) :+ 0.65) VS" $ whnf makeImageVS g,
      bench "f1 ((-0.4) :+ 0.65) RSU" $ whnf makeImageRSU g,
      bench "f1 ((-0.4) :+ 0.65) RPU" $ whnf makeImageRPU g,
      bench "f1 ((-0.4) :+ 0.65) RSS" $ whnf makeImageRSS g,
      bench "f1 ((-0.4) :+ 0.65) RPS" $ whnf makeImageRPS g
    ]
  ]
  where
    default_func = f1 ((-0.4) :+ 0.65)
    g = pixelToJuliaSetValue default_func 1.5 width height 100
    -- makeImageVU = makeImageR VU (width, height)
    -- makeImageVS = makeImageR VS (width, height)
    makeImageRSU = makeImageR RSU (width, height)
    makeImageRPU = makeImageR RPU (width, height)
    makeImageRSS = makeImageR RSS (width, height)
    makeImageRPS = makeImageR RPS (width, height)
    width = 10000
    height = 10000

-- takes string cli args and returns processed args for visualizeJuliaSet
runWithCliArgs :: [String] -> IO()
runWithCliArgs args
  -- | (length args) < 5 = do
  --     usage
  --     putStrLn "using default values " ++ (show (
  --       default_func, 1.5, 1000, 1000, 100, "images/output.png"))
  --     (default_func, 1.5, 1000, 1000, 100, "images/output.png")
  | (length args) < 6 = visualizeJuliaSet default_func r width height
                          maxIter outputFilename "RPU"
  | (length args) < 7 = visualizeJuliaSet default_func r width height
                          maxIter outputFilename arrayType
  | (length args) < 8 = visualizeJuliaSet (parseFunctionParams func_num [])
                          r width height maxIter outputFilename arrayType
  | (length args) < 9 = visualizeJuliaSet
                          (parseFunctionParams func_num params)
                          r width height maxIter outputFilename arrayType
  where
    default_func = f1 ((-0.4) :+ 0.65)
    r = read (args !! 0) :: Double
    width = read (args !! 1) :: Int
    height = read (args !! 2) :: Int
    maxIter = read (args !! 3) :: Int
    outputFilename = args !! 4
    arrayType = args !! 5
    func_num = read (args!!6) :: Int
    params = Prelude.map (read::String -> Complex Double)
                         (slice 7 (length args) args)

slice :: Int -> Int -> [a] -> [a]
slice from to xs = take (to - from + 1) (drop from xs)

parseFunctionParams :: Int
                    -> [Complex Double]
                    -> (Complex Double -> Complex Double)
parseFunctionParams func_num params
    | func_num == 1 && params == [] = f1 ((-0.4) :+ 0.65)
    | func_num == 1 = f1 (params!!0)

usage = do
  let prog = "stack run"
  -- prog <- getProgName
  putStrLn ("Usages: ")
  putStrLn ("    " ++ prog)
  putStrLn ("    " ++ prog ++ " escape_radius width height "
            ++ "max_iter output_filename [arrType]")
  putStrLn ("              [func_num [constant1 [constant2] ... ]]")
  putStrLn ("    " ++ prog ++ " -- --output benchmark.html")
  putStrLn ""
  putStrLn "usage 1 (no argments) displays this help"
  putStrLn ("usage 2 visualizes the julia set created with the " ++
           "given parameters")
  putStrLn "usage 3 benchmarks the program"
  putStrLn ""
  putStrLn "notes on usage 2:"
  putStrLn ("number of constants must correspond to the number of constants "
    ++ "function func_num takes")
  putStrLn ("example: " ++ prog ++ " 1.5 1000 1000 100 "
    ++ "\"images/output.png\"")
  putStrLn ("example: " ++ prog ++ " 1.5 1000 1000 100 "
    ++ "\"images/output.png\" RPU 1 \"(0.285 :+ 0)\"")