text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE TypeOperators #-}

{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise #-}

{-|
Module      : Numeric.LinearAlgebra.MatrixFree.Static
Description : Compose matrix-free linear maps
Copyright   : (c) Ryan Orendorff, 2020
License     : BSD3
Stability   : experimental

Enables the user to compose matrix-free `LinearMap` objects in the same way
that one would compose regular matrices.
-}
module Numeric.LinearAlgebra.MatrixFree.Combinators
    ( (===)
    , (|||)
    , blockMatrix
    , blockDiag
    , LinearMapList(..)
    , applyLinearMapList
    )
where

import           GHC.TypeLits
import           Numeric.LinearAlgebra          ( Transposable(..) )
import           Numeric.LinearAlgebra.Static   ( (#)
                                                , split
                                                , R
                                                )
import           Numeric.LinearAlgebra.MatrixFree
                                                ( LinearMap(..), (##>))
import           Numeric.LinearAlgebra.Static.RList (Sum, RList(..))

-- Stacking matrices either horizontally or vertically are linked through a
-- transpose, so the underlying functions are shared between the two
-- stacking cases. Only in a switched order.
stackLinearMaps :: (KnownNat m, KnownNat n)
                => LinearMap m p -> LinearMap n p -> R p -> R (m + n)
stackLinearMaps la lb v = la ##> v # lb ##> v

addLinearMaps :: (KnownNat n, KnownNat p)
              => LinearMap m n -> LinearMap m p -> R (n + p) -> R m
addLinearMaps la lb v = let (v1, v2) = split v in la ##> v1 + lb ##> v2

-- | Vertically concatenate two `LinearMap`s \(\begin{bmatrix}A\\B\end{bmatrix}\)
(===) :: (KnownNat m, KnownNat n, KnownNat p)
      => LinearMap m n -> LinearMap p n -> LinearMap (m + p) n
(===) la lb = LinearMap (stackLinearMaps la lb) (addLinearMaps (tr la) (tr lb))

infixl 2 ===

-- | Horizontally concatenate two `LinearMap`s \(\begin{bmatrix}A&B\end{bmatrix}\)
(|||) :: (KnownNat m, KnownNat n, KnownNat p)
      => LinearMap m n -> LinearMap m p -> LinearMap m (n + p)
(|||) la lb = tr (tr la === tr lb)

infixl 3 |||

-- | Generate a block matrix from four `LinearMap`s
-- \[
-- \begin{bmatrix}
--   A & B \\
--   C & D
-- \end{bmatrix}
-- \]
blockMatrix
    :: (KnownNat m, KnownNat n, KnownNat p, KnownNat q)
    => LinearMap m n             -- ^ Upper left matrix
    -> LinearMap m p             -- ^ Upper right matrix
    -> LinearMap q n             -- ^ Lower left matrix
    -> LinearMap q p             -- ^ Lower right matrix
    -> LinearMap (m + q) (n + p) -- ^ Resulting block matrix
blockMatrix a b c d = a ||| b === c ||| d

-- | A non-empty list of "LinearMap" matrices of different sizes.
data LinearMapList (ms :: [Nat]) (ns :: [Nat]) where
    LOne :: (KnownNat m, KnownNat n) => LinearMap m n -> LinearMapList '[m] '[n]
    (:-:) :: (KnownNat m, KnownNat n) => LinearMap m n -> LinearMapList ms ns
                                      -> LinearMapList (m ': ms) (n ': ns)

-- | Concatenate a "LinearMapList" of "LinearMap"s as a block diagonal.
-- Note that this has poor performance; it runs `split` on the
-- input vector O(n) times, where n is the number of entries in the list.
blockDiag :: LinearMapList ms ns -> LinearMap (Sum ms) (Sum ns)
blockDiag ls = case blockDiagWorker ls of
  LWorker lsW -> lsW

data LWorker m n where
  LWorker :: (KnownNat m, KnownNat n) => LinearMap m n -> LWorker m n

blockDiagWorker :: LinearMapList ms ns -> LWorker (Sum ms) (Sum ns)
blockDiagWorker (LOne l) = LWorker l
blockDiagWorker (l :-: ls) = case blockDiagWorker ls of
  LWorker lsW -> LWorker (diag2 l lsW)

-- Combine two `LinearMap`s on a diagonal without having to use
-- konst 0, which would just be a waste of processing time (sum a list and
-- multiply it by zero)
diag2 :: (KnownNat m, KnownNat n, KnownNat p, KnownNat q)
      => LinearMap m n
      -> LinearMap p q
      -> LinearMap (m + p) (n + q)
diag2 la lb = LinearMap (f la lb) (f (tr la) (tr lb))
  where
    f :: (KnownNat m', KnownNat n', KnownNat p', KnownNat q')
      => LinearMap m' n' -> LinearMap p' q' -> R (n' + q') -> R (m' + p')
    f l1 l2 v = let (v1, v2) = split v in
      (l1 ##> v1) # (l2 ##> v2)

-- | Apply a list of "LinearMap"s over a list of static vectors
applyLinearMapList :: LinearMapList ms ns -> RList ns -> RList ms
applyLinearMapList (LOne _) (_ ::: _) = error "This case is impossible"
applyLinearMapList (_ :-: _) (ROne _) = error "This case is impossible"
applyLinearMapList (LOne l) (ROne r) = ROne (l ##> r)
applyLinearMapList (l :-: ls) (r ::: rs) = (l ##> r) ::: applyLinearMapList ls rs

max_stars_repo_name: ryanorendorff/convex

max_issues_repo_name: ryanorendorff/convex

max_forks_repo_name: ryanorendorff/convex


text:
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE PolyKinds           #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TemplateHaskell     #-}
{-# LANGUAGE TypeFamilies        #-}

{-# OPTIONS_GHC -Wno-orphans -Wno-deprecations #-}
-- {-# OPTIONS_GHC -ddump-splices -dsuppress-module-prefixes #-}

-- | This module is a modified version of Generics.SOP.Instances.
-- mostly the same.
module Staged.GHC.Generics.Instances () where

import Staged.GHC.Generics.TH

import Control.Exception
import Data.Char
import Data.Complex
import Data.Data
import Data.Fixed
import Data.Functor.Compose
import qualified Data.Functor.Const
import Data.Functor.Identity
import Data.Functor.Product
import Data.Functor.Sum
import Data.List.NonEmpty
import qualified Data.Monoid
import Data.Ord
import qualified Data.Semigroup
import Data.Version
import Data.Void
import Foreign.C.Error
import Foreign.C.Types
import GHC.ByteOrder
import GHC.Conc
import GHC.ExecutionStack
import GHC.Exts
-- import GHC.Events -- platform-specific, omitted
import GHC.Fingerprint
import GHC.Float
import qualified GHC.Generics
import GHC.IO.Buffer
import GHC.IO.Device
import GHC.IO.Encoding
import GHC.IO.Encoding.Failure
import GHC.IO.Exception
import GHC.IO.Handle
import GHC.RTS.Flags
import qualified GHC.Stack
import GHC.StaticPtr
import GHC.Stats
import System.Console.GetOpt
import System.IO
import Text.Printf
import Text.Read.Lex

-- Types from the Prelude:

-- there is manual instance for this
-- deriveGeneric ''Bool

deriveGeneric ''Ordering
deriveGeneric ''Maybe
deriveGeneric ''Either
deriveGeneric ''()
deriveGeneric ''(,)              -- 2
deriveGeneric ''(,,)
deriveGeneric ''(,,,)
deriveGeneric ''(,,,,)           -- 5
deriveGeneric ''(,,,,,)
-- deriveGeneric ''(,,,,,,)
-- deriveGeneric ''(,,,,,,,)
-- deriveGeneric ''(,,,,,,,,)
-- deriveGeneric ''(,,,,,,,,,)      -- 10
-- deriveGeneric ''(,,,,,,,,,,)
-- deriveGeneric ''(,,,,,,,,,,,)
-- deriveGeneric ''(,,,,,,,,,,,,)
-- deriveGeneric ''(,,,,,,,,,,,,,)
-- deriveGeneric ''(,,,,,,,,,,,,,,) -- 15
-- deriveGeneric ''(,,,,,,,,,,,,,,,)
-- deriveGeneric ''(,,,,,,,,,,,,,,,,)
-- deriveGeneric ''(,,,,,,,,,,,,,,,,,)
-- deriveGeneric ''(,,,,,,,,,,,,,,,,,,)
-- deriveGeneric ''(,,,,,,,,,,,,,,,,,,,) -- 20
-- deriveGeneric ''(,,,,,,,,,,,,,,,,,,,,)
-- deriveGeneric ''(,,,,,,,,,,,,,,,,,,,,,)
-- deriveGeneric ''(,,,,,,,,,,,,,,,,,,,,,,)
-- deriveGeneric ''(,,,,,,,,,,,,,,,,,,,,,,,)
-- deriveGeneric ''(,,,,,,,,,,,,,,,,,,,,,,,,) -- 25
-- deriveGeneric ''(,,,,,,,,,,,,,,,,,,,,,,,,,)
-- deriveGeneric ''(,,,,,,,,,,,,,,,,,,,,,,,,,,)
-- deriveGeneric ''(,,,,,,,,,,,,,,,,,,,,,,,,,,,)
-- deriveGeneric ''(,,,,,,,,,,,,,,,,,,,,,,,,,,,,)
-- deriveGeneric ''(,,,,,,,,,,,,,,,,,,,,,,,,,,,,,) -- 30
-- deriveGeneric ''[]

deriveGeneric1 ''Maybe
deriveGeneric1 ''Either
deriveGeneric1 ''(,)              -- 2
deriveGeneric1 ''(,,)
deriveGeneric1 ''(,,,)
deriveGeneric1 ''(,,,,)           -- 5
deriveGeneric1 ''(,,,,,)
deriveGeneric1 ''(,,,,,,)
deriveGeneric1 ''(,,,,,,,)
deriveGeneric1 ''(,,,,,,,,)
deriveGeneric1 ''(,,,,,,,,,)      -- 10
deriveGeneric1 ''(,,,,,,,,,,)
deriveGeneric1 ''(,,,,,,,,,,,)
deriveGeneric1 ''(,,,,,,,,,,,,)
deriveGeneric1 ''(,,,,,,,,,,,,,)
deriveGeneric1 ''(,,,,,,,,,,,,,,) -- 15
deriveGeneric1 ''(,,,,,,,,,,,,,,,)
deriveGeneric1 ''(,,,,,,,,,,,,,,,,)
deriveGeneric1 ''(,,,,,,,,,,,,,,,,,)
deriveGeneric1 ''(,,,,,,,,,,,,,,,,,,)
deriveGeneric1 ''(,,,,,,,,,,,,,,,,,,,) -- 20
deriveGeneric1 ''(,,,,,,,,,,,,,,,,,,,,)
deriveGeneric1 ''(,,,,,,,,,,,,,,,,,,,,,)
deriveGeneric1 ''(,,,,,,,,,,,,,,,,,,,,,,)
deriveGeneric1 ''(,,,,,,,,,,,,,,,,,,,,,,,)
deriveGeneric1 ''(,,,,,,,,,,,,,,,,,,,,,,,,) -- 25
deriveGeneric1 ''(,,,,,,,,,,,,,,,,,,,,,,,,,)
deriveGeneric1 ''(,,,,,,,,,,,,,,,,,,,,,,,,,,)
deriveGeneric1 ''(,,,,,,,,,,,,,,,,,,,,,,,,,,,)
deriveGeneric1 ''(,,,,,,,,,,,,,,,,,,,,,,,,,,,,)
deriveGeneric1 ''(,,,,,,,,,,,,,,,,,,,,,,,,,,,,,) -- 30
deriveGeneric1 ''[]

-- Other types from base:

-- From Control.Exception:
deriveGeneric ''IOException
deriveGeneric ''ArithException
deriveGeneric ''ArrayException
deriveGeneric ''AssertionFailed
deriveGeneric ''AsyncException
deriveGeneric ''NonTermination
deriveGeneric ''NestedAtomically
deriveGeneric ''BlockedIndefinitelyOnMVar
deriveGeneric ''BlockedIndefinitelyOnSTM
deriveGeneric ''AllocationLimitExceeded
deriveGeneric ''Deadlock
deriveGeneric ''NoMethodError
deriveGeneric ''PatternMatchFail
deriveGeneric ''RecConError
deriveGeneric ''RecSelError
deriveGeneric ''RecUpdError
deriveGeneric ''ErrorCall
deriveGeneric ''TypeError
deriveGeneric ''MaskingState

-- From Data.Char:
deriveGeneric ''GeneralCategory

-- From Data.Complex:
deriveGeneric ''Complex

-- From Data.Data:
deriveGeneric ''DataRep
deriveGeneric ''Fixity
deriveGeneric ''ConstrRep

-- From Data.Fixed:
deriveGeneric ''Fixed
deriveGeneric ''E0
deriveGeneric ''E1
deriveGeneric ''E2
deriveGeneric ''E3
deriveGeneric ''E6
deriveGeneric ''E9
deriveGeneric ''E12

-- From Data.Functor.Compose
deriveGeneric ''Compose
deriveGeneric1 ''Compose

-- From Data.Functor.Const
deriveGeneric ''Data.Functor.Const.Const
deriveGeneric1 ''Data.Functor.Const.Const

-- From Data.Functor.Identity
deriveGeneric ''Identity
deriveGeneric1 ''Identity

-- From Data.Functor.Product
deriveGeneric ''Product
deriveGeneric1 ''Product

-- From Data.Functor.Sum
deriveGeneric ''Sum
deriveGeneric1 ''Sum

-- From Data.List.NonEmpty
deriveGeneric ''NonEmpty
deriveGeneric1 ''NonEmpty

-- From Data.Monoid:
deriveGeneric ''Data.Monoid.Dual
deriveGeneric ''Data.Monoid.Endo
deriveGeneric ''Data.Monoid.All
deriveGeneric ''Data.Monoid.Any
deriveGeneric ''Data.Monoid.Sum
deriveGeneric ''Data.Monoid.Product
deriveGeneric ''Data.Monoid.First
deriveGeneric ''Data.Monoid.Last
deriveGeneric ''Data.Monoid.Alt

deriveGeneric1 ''Data.Monoid.Dual
-- deriveGeneric1 ''Data.Monoid.Endo
deriveGeneric1 ''Data.Monoid.Sum
deriveGeneric1 ''Data.Monoid.Product
deriveGeneric1 ''Data.Monoid.First
deriveGeneric1 ''Data.Monoid.Last
deriveGeneric1 ''Data.Monoid.Alt

-- From Data.Ord:
deriveGeneric ''Down
deriveGeneric1 ''Down

-- From Data.Proxy:
deriveGeneric ''Proxy
deriveGeneric1 ''Proxy

-- From Data.Semigroup:
deriveGeneric ''Data.Semigroup.Min
deriveGeneric ''Data.Semigroup.Max
deriveGeneric ''Data.Semigroup.First
deriveGeneric ''Data.Semigroup.Last
deriveGeneric ''Data.Semigroup.WrappedMonoid
deriveGeneric ''Data.Semigroup.Option
deriveGeneric ''Data.Semigroup.Arg

deriveGeneric1 ''Data.Semigroup.Min
deriveGeneric1 ''Data.Semigroup.Max
deriveGeneric1 ''Data.Semigroup.First
deriveGeneric1 ''Data.Semigroup.Last
deriveGeneric1 ''Data.Semigroup.WrappedMonoid
deriveGeneric1 ''Data.Semigroup.Option
deriveGeneric1 ''Data.Semigroup.Arg

-- From Data.Version:
deriveGeneric ''Version

-- From Data.Void:
deriveGeneric ''Void

-- From Foreign.C.Error:
deriveGeneric ''Errno

-- From Foreign.C.Types:
deriveGeneric ''CChar
deriveGeneric ''CSChar
deriveGeneric ''CUChar
deriveGeneric ''CShort
deriveGeneric ''CUShort
deriveGeneric ''CInt
deriveGeneric ''CUInt
deriveGeneric ''CLong
deriveGeneric ''CULong
deriveGeneric ''CPtrdiff
deriveGeneric ''CSize
deriveGeneric ''CWchar
deriveGeneric ''CSigAtomic
deriveGeneric ''CLLong
deriveGeneric ''CULLong
deriveGeneric ''CIntPtr
deriveGeneric ''CUIntPtr
deriveGeneric ''CIntMax
deriveGeneric ''CUIntMax
deriveGeneric ''CClock
deriveGeneric ''CTime
deriveGeneric ''CUSeconds
deriveGeneric ''CSUSeconds
deriveGeneric ''CFloat
deriveGeneric ''CDouble

-- From GHC.ByteOrder:
deriveGeneric ''ByteOrder

-- From GHC.Conc:
deriveGeneric ''ThreadStatus
deriveGeneric ''BlockReason

-- From GHC.ExecutionStack:
deriveGeneric ''Location
deriveGeneric ''SrcLoc

-- From GHC.Exts:
deriveGeneric ''RuntimeRep
deriveGeneric ''VecCount
deriveGeneric ''VecElem

-- From GHC.Generics:
deriveGeneric ''GHC.Generics.K1
deriveGeneric ''GHC.Generics.U1
deriveGeneric ''GHC.Generics.V1
deriveGeneric ''GHC.Generics.Par1
deriveGeneric ''GHC.Generics.M1
deriveGeneric ''GHC.Generics.R
deriveGeneric ''GHC.Generics.S
deriveGeneric ''GHC.Generics.D
deriveGeneric ''GHC.Generics.C
deriveGeneric ''(GHC.Generics.:*:)
deriveGeneric ''(GHC.Generics.:+:)
deriveGeneric ''(GHC.Generics.:.:)
deriveGeneric ''GHC.Generics.Associativity
deriveGeneric ''GHC.Generics.DecidedStrictness
deriveGeneric ''GHC.Generics.SourceStrictness
deriveGeneric ''GHC.Generics.SourceUnpackedness
deriveGeneric ''GHC.Generics.Fixity

-- From GHC.IO.Buffer:
deriveGeneric ''Buffer
deriveGeneric ''BufferState

-- From GHC.IO.Device:
deriveGeneric ''IODeviceType

-- From GHC.IO.Encoding:
deriveGeneric ''BufferCodec
deriveGeneric ''CodingProgress

-- From GHC.IO.Encoding.Failure:
deriveGeneric ''CodingFailureMode

-- From GHC.Fingerprint
deriveGeneric ''Fingerprint

-- From GHC.Float
deriveGeneric ''FFFormat

-- From GHC.IO.Exception:
deriveGeneric ''FixIOException
deriveGeneric ''IOErrorType

-- From GHC.IO.Handle:
deriveGeneric ''HandlePosn
deriveGeneric ''LockMode

-- From GHC.RTS.Flags:
deriveGeneric ''RTSFlags
deriveGeneric ''GiveGCStats
deriveGeneric ''GCFlags
deriveGeneric ''ConcFlags
deriveGeneric ''MiscFlags
deriveGeneric ''DebugFlags
deriveGeneric ''DoCostCentres
deriveGeneric ''CCFlags
deriveGeneric ''DoHeapProfile
deriveGeneric ''ProfFlags
deriveGeneric ''DoTrace
deriveGeneric ''TraceFlags
deriveGeneric ''TickyFlags
deriveGeneric ''ParFlags

-- From GHC.Stack:
deriveGeneric ''GHC.Stack.SrcLoc
deriveGeneric ''GHC.Stack.CallStack

-- From GHC.StaticPtr:
deriveGeneric ''StaticPtrInfo

-- From GHC.Stats:
deriveGeneric ''RTSStats
deriveGeneric ''GCDetails

-- From System.Console.GetOpt:

deriveGeneric ''ArgOrder
deriveGeneric ''OptDescr
deriveGeneric ''ArgDescr

-- From System.Exit:

deriveGeneric ''ExitCode

-- From System.IO:

deriveGeneric ''IOMode
deriveGeneric ''BufferMode
deriveGeneric ''SeekMode
deriveGeneric ''Newline
deriveGeneric ''NewlineMode

-- From Text.Printf:

deriveGeneric ''FieldFormat
deriveGeneric ''FormatAdjustment
deriveGeneric ''FormatSign
deriveGeneric ''FormatParse

-- From Text.Read.Lex:

deriveGeneric ''Lexeme
deriveGeneric ''Number

-- Abstract / primitive datatypes (we don't derive Generic for these):
--
-- Ratio
-- Integer
-- ThreadId
-- Chan
-- MVar
-- QSem
-- QSemN
-- DataType
-- Dynamic
-- IORef
-- TypeRep
-- TyCon
-- TypeRepKey
-- KProxy -- not abstract, but intended for kind-level use
-- STRef
-- Unique
-- ForeignPtr
-- CFile
-- CFpos
-- CJmpBuf
-- Pool
-- Ptr
-- FunPtr
-- IntPtr
-- WordPtr
-- StablePtr
-- Char
-- Double
-- Float
-- Int
-- Int8
-- Int16
-- Int32
-- Int64
-- Word
-- Word8
-- Word16
-- Word32
-- Word64
-- IO
-- ST
-- (->)
-- RealWorld
-- Handle
-- HandlePosn
-- TextEncoding
-- StableName
-- Weak
-- ReadP
-- ReadPrec
-- STM
-- TVar
-- Natural
-- Event
-- EventManager
-- CostCentre
-- CostCentreStack
--
-- Datatypes we cannot currently handle:
--
-- SomeException
-- SomeAsyncException
-- Handler
-- Coercion
-- (:~:)

ext: hs

max_stars_repo_name: phadej/staged-gg

max_issues_repo_name: phadej/staged-gg

max_forks_repo_name: phadej/staged-gg


text:
-- |
-- Module      : Control.Monad.Bayes.Class
-- Description : Types for probabilistic modelling
-- Copyright   : (c) Adam Scibior, 2015-2020
-- License     : MIT
-- Maintainer  : leonhard.markert@tweag.io
-- Stability   : experimental
-- Portability : GHC
--
-- This module defines 'MonadInfer', which can be used to represent a simple model
-- like the following:
--
-- @
-- import Control.Monad (when)
-- import Control.Monad.Bayes.Class
--
-- model :: MonadInfer m => m Bool
-- model = do
--   rain <- bernoulli 0.3
--   sprinkler <-
--     bernoulli $
--     if rain
--       then 0.1
--       else 0.4
--   let wetProb =
--     case (rain, sprinkler) of
--       (True,  True)  -> 0.98
--       (True,  False) -> 0.80
--       (False, True)  -> 0.90
--       (False, False) -> 0.00
--   score wetProb
--   return rain
-- @
module Control.Monad.Bayes.Class
  ( MonadSample,
    random,
    uniform,
    normal,
    gamma,
    beta,
    bernoulli,
    categorical,
    logCategorical,
    uniformD,
    geometric,
    poisson,
    dirichlet,
    MonadCond,
    score,
    factor,
    condition,
    MonadInfer,
    normalPdf,
  )
where

import Control.Monad (when)
import Control.Monad.Trans.Class
import Control.Monad.Trans.Cont
import Control.Monad.Trans.Identity
import Control.Monad.Trans.List
import Control.Monad.Trans.Maybe
import Control.Monad.Trans.RWS hiding (tell)
import Control.Monad.Trans.Reader
import Control.Monad.Trans.State
import Control.Monad.Trans.Writer
import qualified Data.Vector as V
import Data.Vector.Generic as VG
import Numeric.Log
import Statistics.Distribution
import Statistics.Distribution.Beta (betaDistr)
import Statistics.Distribution.Gamma (gammaDistr)
import Statistics.Distribution.Geometric (geometric0)
import Statistics.Distribution.Normal (normalDistr)
import qualified Statistics.Distribution.Poisson as Poisson
import Statistics.Distribution.Uniform (uniformDistr)

-- | Monads that can draw random variables.
class Monad m => MonadSample m where
  -- | Draw from a uniform distribution.
  random ::
    -- | \(\sim \mathcal{U}(0, 1)\)
    m Double

  -- | Draw from a uniform distribution.
  uniform ::
    -- | lower bound a
    Double ->
    -- | upper bound b
    Double ->
    -- | \(\sim \mathcal{U}(a, b)\).
    m Double
  uniform a b = draw (uniformDistr a b)

  -- | Draw from a normal distribution.
  normal ::
    -- | mean μ
    Double ->
    -- | standard deviation σ
    Double ->
    -- | \(\sim \mathcal{N}(\mu, \sigma^2)\)
    m Double
  normal m s = draw (normalDistr m s)

  -- | Draw from a gamma distribution.
  gamma ::
    -- | shape k
    Double ->
    -- | scale θ
    Double ->
    -- | \(\sim \Gamma(k, \theta)\)
    m Double
  gamma shape scale = draw (gammaDistr shape scale)

  -- | Draw from a beta distribution.
  beta ::
    -- | shape α
    Double ->
    -- | shape β
    Double ->
    -- | \(\sim \mathrm{Beta}(\alpha, \beta)\)
    m Double
  beta a b = draw (betaDistr a b)

  -- | Draw from a Bernoulli distribution.
  bernoulli ::
    -- | probability p
    Double ->
    -- | \(\sim \mathrm{B}(1, p)\)
    m Bool
  bernoulli p = fmap (< p) random

  -- | Draw from a categorical distribution.
  categorical ::
    Vector v Double =>
    -- | event probabilities
    v Double ->
    -- | outcome category
    m Int
  categorical ps = fromPMF (ps !)

  -- | Draw from a categorical distribution in the log domain.
  logCategorical ::
    (Vector v (Log Double), Vector v Double) =>
    -- | event probabilities
    v (Log Double) ->
    -- | outcome category
    m Int
  logCategorical = categorical . VG.map (exp . ln)

  -- | Draw from a discrete uniform distribution.
  uniformD ::
    -- | observable outcomes @xs@
    [a] ->
    -- | \(\sim \mathcal{U}\{\mathrm{xs}\}\)
    m a
  uniformD xs = do
    let n = Prelude.length xs
    i <- categorical $ V.replicate n (1 / fromIntegral n)
    return (xs !! i)

  -- | Draw from a geometric distribution.
  geometric ::
    -- | success rate p
    Double ->
    -- | \(\sim\) number of failed Bernoulli trials with success probability p before first success
    m Int
  geometric = discrete . geometric0

  -- | Draw from a Poisson distribution.
  poisson ::
    -- | parameter λ
    Double ->
    -- | \(\sim \mathrm{Pois}(\lambda)\)
    m Int
  poisson = discrete . Poisson.poisson

  -- | Draw from a Dirichlet distribution.
  dirichlet ::
    Vector v Double =>
    -- | concentration parameters @as@
    v Double ->
    -- | \(\sim \mathrm{Dir}(\mathrm{as})\)
    m (v Double)
  dirichlet as = do
    xs <- VG.mapM (`gamma` 1) as
    let s = VG.sum xs
    let ys = VG.map (/ s) xs
    return ys

-- | Draw from a continuous distribution using the inverse cumulative density
-- function.
draw :: (ContDistr d, MonadSample m) => d -> m Double
draw d = fmap (quantile d) random

-- | Draw from a discrete distribution using a sequence of draws from
-- Bernoulli.
fromPMF :: MonadSample m => (Int -> Double) -> m Int
fromPMF p = f 0 1
  where
    f i r = do
      when (r < 0) $ error "fromPMF: total PMF above 1"
      let q = p i
      when (q < 0 || q > 1) $ error "fromPMF: invalid probability value"
      b <- bernoulli (q / r)
      if b then pure i else f (i + 1) (r - q)

-- | Draw from a discrete distributions using the probability mass function.
discrete :: (DiscreteDistr d, MonadSample m) => d -> m Int
discrete = fromPMF . probability

-- | Monads that can score different execution paths.
class Monad m => MonadCond m where
  -- | Record a likelihood.
  score ::
    -- | likelihood of the execution path
    Log Double ->
    m ()

-- | Synonym for 'score'.
factor ::
  MonadCond m =>
  -- | likelihood of the execution path
  Log Double ->
  m ()
factor = score

-- | Hard conditioning.
condition :: MonadCond m => Bool -> m ()
condition b = score $ if b then 1 else 0

-- | Monads that support both sampling and scoring.
class (MonadSample m, MonadCond m) => MonadInfer m

-- | Probability density function of the normal distribution.
normalPdf ::
  -- | mean μ
  Double ->
  -- | standard deviation σ
  Double ->
  -- | sample x
  Double ->
  -- | relative likelihood of observing sample x in \(\mathcal{N}(\mu, \sigma^2)\)
  Log Double
normalPdf mu sigma x = Exp $ logDensity (normalDistr mu sigma) x

----------------------------------------------------------------------------
-- Instances that lift probabilistic effects to standard tranformers.

instance MonadSample m => MonadSample (IdentityT m) where
  random = lift random
  bernoulli = lift . bernoulli

instance MonadCond m => MonadCond (IdentityT m) where
  score = lift . score

instance MonadInfer m => MonadInfer (IdentityT m)

instance MonadSample m => MonadSample (MaybeT m) where
  random = lift random

instance MonadCond m => MonadCond (MaybeT m) where
  score = lift . score

instance MonadInfer m => MonadInfer (MaybeT m)

instance MonadSample m => MonadSample (ReaderT r m) where
  random = lift random
  bernoulli = lift . bernoulli

instance MonadCond m => MonadCond (ReaderT r m) where
  score = lift . score

instance MonadInfer m => MonadInfer (ReaderT r m)

instance (Monoid w, MonadSample m) => MonadSample (WriterT w m) where
  random = lift random
  bernoulli = lift . bernoulli
  categorical = lift . categorical

instance (Monoid w, MonadCond m) => MonadCond (WriterT w m) where
  score = lift . score

instance (Monoid w, MonadInfer m) => MonadInfer (WriterT w m)

instance MonadSample m => MonadSample (StateT s m) where
  random = lift random
  bernoulli = lift . bernoulli
  categorical = lift . categorical

instance MonadCond m => MonadCond (StateT s m) where
  score = lift . score

instance MonadInfer m => MonadInfer (StateT s m)

instance (MonadSample m, Monoid w) => MonadSample (RWST r w s m) where
  random = lift random

instance (MonadCond m, Monoid w) => MonadCond (RWST r w s m) where
  score = lift . score

instance (MonadInfer m, Monoid w) => MonadInfer (RWST r w s m)

instance MonadSample m => MonadSample (ListT m) where
  random = lift random
  bernoulli = lift . bernoulli
  categorical = lift . categorical

instance MonadCond m => MonadCond (ListT m) where
  score = lift . score

instance MonadInfer m => MonadInfer (ListT m)

instance MonadSample m => MonadSample (ContT r m) where
  random = lift random

instance MonadCond m => MonadCond (ContT r m) where
  score = lift . score

instance MonadInfer m => MonadInfer (ContT r m)

max_stars_repo_name: saeedhadikhanloo/monad-bayes

max_issues_repo_name: saeedhadikhanloo/monad-bayes

max_forks_repo_name: saeedhadikhanloo/monad-bayes


text:
{-# LANGUAGE ForeignFunctionInterface #-}
{-|
Module      : Grenade.Layers.Internal.Activations
Description : Fast activation functions that call efficient C function
Maintainer  : Theo Charalambous
License     : BSD2
Stability   : experimental
-}
module Grenade.Layers.Internal.Activations
  ( relu
  , relu1d
  , leakyRelu
  , leakyRelu1d
  ) where

import qualified Data.Vector.Storable        as U (unsafeFromForeignPtr0,
                                                   unsafeToForeignPtr0)

import           Foreign                     (mallocForeignPtrArray, withForeignPtr)
import           Foreign.Ptr                 (Ptr)
import           Numeric.LinearAlgebra       (Matrix, Vector, flatten)
import qualified Numeric.LinearAlgebra.Devel as U
import           System.IO.Unsafe            (unsafePerformIO)

import           Grenade.Types

-- | Apply ReLU on a matrix
relu :: Int -> Int -> Int -> Matrix RealNum -> Matrix RealNum
relu channels rows cols = leakyRelu channels rows cols 0
{-# INLINE relu #-}

-- | Apply ReLU on a vector
relu1d :: Int -> Vector RealNum -> Vector RealNum
relu1d size = leakyRelu1d size 0
{-# INLINE relu1d #-}

-- | Apply LeakyReLU on a matrix
leakyRelu :: Int -> Int -> Int -> RealNum -> Matrix RealNum -> Matrix RealNum
leakyRelu channels rows cols alpha m
  = let vec = flatten m
        out = leakyRelu1d (channels * rows * cols) alpha vec
    in  U.matrixFromVector U.RowMajor (rows * channels) cols out
{-# INLINE leakyRelu #-}

-- | Apply LeakyReLU on a vector
leakyRelu1d :: Int -> RealNum -> Vector RealNum -> Vector RealNum
leakyRelu1d size alpha vec
  = unsafePerformIO $ do
      outPtr        <- mallocForeignPtrArray size
      let (inPtr, _) = U.unsafeToForeignPtr0 vec

      withForeignPtr inPtr $ \inPtr' ->
        withForeignPtr outPtr $ \outPtr' ->
          leaky_relu_forward inPtr' size alpha outPtr'

      return $ U.unsafeFromForeignPtr0 outPtr size
{-# INLINE leakyRelu1d #-}

foreign import ccall unsafe
   leaky_relu_forward
      :: Ptr RealNum -> Int -> RealNum -> Ptr RealNum -> IO ()

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
{-# LANGUAGE CPP                   #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE OverloadedStrings     #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE PolyKinds             #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE TypeApplications      #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE AllowAmbiguousTypes   #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE UndecidableInstances  #-}
{-# LANGUAGE ConstraintKinds #-}
module Frames.VegaLite.Regression
  (
    Frame2DRegressionScatterFit (..)
  , frameScatterWithFit
  , keyedLayeredFrameScatterWithFit
  , scatterWithFit
  , FitToPlot (..)
  , regressionCoefficientPlot
  , regressionCoefficientPlotMany
  ) where

import qualified Frames.VegaLite.Utils as FV
import qualified Frames.Transform as FT
import qualified Frames.Regression as FR
import qualified Math.Regression.Regression as RE

import qualified Control.Foldl          as FL
import qualified Data.Text              as T
import qualified Data.Vinyl             as V
import qualified Data.Vinyl.Functor     as V
import qualified Data.Vinyl.TypeLevel   as V
import qualified Frames                 as F
import qualified Frames.Melt            as F
import qualified Graphics.Vega.VegaLite as GV
import qualified Data.List              as List
import Text.Printf (printf)

import qualified Statistics.Types               as S

#if MIN_VERSION_hvega(0,4,0)
gvTitle :: Text -> GV.PropertySpec
gvTitle x = GV.title x []
#else
gvTitle :: Text -> (GV.VLProperty, GV.VLSpec)
gvTitle = GV.title
#endif

-- | Plot regression coefficents with error bars
-- | Flex version handles a foldable of results so we can, e.g.,
-- | 1. Compare across time or other variable in the data
-- | 2. Compare different fitting methods for the same data

-- TODO: Fix replicated y-axis ticks.  There since we need the difference between them (trailing "'"s) to provide y-offset
-- to the different results.  Otherwise they would overlap in y.  Maybe fix by switching to layers and only having Y-axis labels
-- on 0th?

-- TODO: Make this typed?  Using a regression result with typed parameters?  Power-to-weight?

regressionCoefficientPlot :: T.Text -> [T.Text] -> RE.RegressionResult Double -> S.CL Double -> GV.VegaLite
regressionCoefficientPlot title names r cl = regressionCoefficientPlotFlex False id title names (V.Identity ("",r)) cl

regressionCoefficientPlotMany :: Foldable f
                              => (k -> T.Text) -> T.Text -> [T.Text] -> f (k, RE.RegressionResult Double) -> S.CL Double -> GV.VegaLite
regressionCoefficientPlotMany = regressionCoefficientPlotFlex True

regressionCoefficientPlotFlex :: Foldable f
                              => Bool -> (k -> Text) -> T.Text -> [T.Text] -> f (k, RE.RegressionResult Double) -> S.CL Double -> GV.VegaLite
regressionCoefficientPlotFlex haveLegend printKey title names results cl =
  let toRow m (RE.NamedEstimate n e eci _) = [("Parameter",GV.Str (n <> T.replicate m "'")), ("Estimate",GV.Number e), ("Confidence",GV.Number eci)]
      addKey k l = ("Key", GV.Str $ printKey k) : l
      dataRowFold = FL.Fold (\(l,n) (k,regRes) -> (l ++ fmap (flip GV.dataRow [] . addKey k . toRow n) (RE.namedEstimates names regRes cl), n+1)) ([],0) (GV.dataFromRows [] . List.concat . fst)
      dat = FL.fold dataRowFold results
      xLabel = "Estimate (with " <> T.pack (printf "%2.0f" (100 * S.confidenceLevel cl)) <> "% confidence error bars)"
      estimateXEnc = GV.position GV.X [GV.PName "Estimate", GV.PmType GV.Quantitative, GV.PAxis [GV.AxTitle xLabel]]
      estimateYEnc = GV.position GV.Y [GV.PName "Parameter", GV.PmType GV.Ordinal]
      handleLegend l = if haveLegend then l else GV.MLegend [] : l
      estimateColorEnc = GV.color $ handleLegend  [GV.MName "Key", GV.MmType GV.Nominal]
      estimateEnc = estimateXEnc . estimateYEnc . estimateColorEnc
      estLoCalc = "datum.Estimate - datum.Confidence/2"
      estHiCalc = "datum.Estimate + datum.Confidence/2"
      calcEstConf = GV.calculateAs estLoCalc "estLo" . GV.calculateAs estHiCalc "estHi"
      estConfLoEnc = GV.position GV.X [GV.PName "estLo", GV.PmType GV.Quantitative, GV.PAxis [GV.AxTitle ""]]
      estConfHiEnc = GV.position GV.X2 [GV.PName "estHi", GV.PmType GV.Quantitative, GV.PAxis [GV.AxTitle ""]]
      estConfEnc = estConfLoEnc . estConfHiEnc . estimateYEnc . estimateColorEnc
      estSpec = GV.asSpec [(GV.encoding . estimateEnc) [], GV.mark GV.Point []]
      confSpec = GV.asSpec [(GV.encoding . estConfEnc) [], GV.mark GV.Rule []]
      configuration = GV.configure
        . GV.configuration (GV.ViewStyle [GV.ViewContinuousWidth 800, GV.ViewContinuousHeight 400]) . GV.configuration (GV.PaddingStyle $ GV.PSize 50)
      vl = GV.toVegaLite
        [
          gvTitle title
        , (GV.transform . calcEstConf) []
        , GV.layer [estSpec, confSpec]
        , dat
        , configuration []
        ]
  in vl

--
{-
type family WeightElemOf (rs :: [(Symbol, Type)]) (w :: (Symbol, Type)) :: Constraint where
  WeightElemOf _ FR.Unweighted = ()
  WeightElemOf rs w  = F.ElemOf rs w
-}

frameRegressionError :: forall y wc as w rs. FR.FrameRegressionResult y wc as w rs -> Error rs
frameRegressionError (FR.FrameUnweightedRegressionResult _) = const 0
frameRegressionError (FR.FrameWeightedRegressionResult wf _) = (\r -> 1/(wf $ F.rgetField @w r))

type ScatterFitConstraints x y = ( F.ColumnHeaders '[x]
                                 , F.ColumnHeaders '[y]
                                 , FV.ToVLDataValue (F.ElField x)
                                 , FV.ToVLDataValue (F.ElField y)
                                 , V.KnownField y
                                 , Real (V.Snd y)
                                 )

type ScatterFitC1 rs x y = ( F.ElemOf (rs V.++ [YError,YFit,YFitError]) x
                           , F.ElemOf (rs V.++ [YError,YFit,YFitError]) y
                           , F.ElemOf (rs V.++ [YError,YFit,YFitError]) YError
                           , F.ElemOf (rs V.++ [YError,YFit,YFitError]) YFit
                           , F.ElemOf (rs V.++ [YError,YFit,YFitError]) YFitError)

class Frame2DRegressionScatterFit rs a where
  regressionResultToError :: a -> Error rs
  regressionResultToFit :: Maybe T.Text -> a -> S.CL Double -> FitToPlot rs
--  scatterPlot :: (Foldable f, Functor f) => T.Text -> Maybe T.Text -> a -> Double -> f (F.Record rs) -> GV.VegaLite
  scatterPlotSpec :: (Foldable f, Functor f) => Maybe T.Text -> a -> S.CL Double -> f (F.Record rs) -> GV.VLSpec

frameScatterWithFit :: (Frame2DRegressionScatterFit rs a, Foldable f, Functor f)
                    => T.Text -> Maybe T.Text -> a -> S.CL Double -> f (F.Record rs) -> GV.VegaLite
frameScatterWithFit title fitNameM frr cl frame =
  let configuration = GV.configure
        . GV.configuration (GV.ViewStyle [GV.ViewContinuousWidth 800, GV.ViewContinuousHeight 400]) . GV.configuration (GV.PaddingStyle $ GV.PSize 50)
      swfSpec = GV.layer [scatterPlotSpec fitNameM frr cl frame]
  in GV.toVegaLite [configuration [], swfSpec, gvTitle title]

keyedLayeredFrameScatterWithFit :: (Frame2DRegressionScatterFit rs a, Foldable f, Functor f, Foldable g, Functor g)
  => T.Text -> (k -> T.Text) -> g (k, a) -> S.CL Double -> f (F.Record rs) -> GV.VegaLite
keyedLayeredFrameScatterWithFit title keyText keyedFits cl dat =
  let toSpec (k, a) = scatterPlotSpec (Just $ keyText k) a cl dat
      specs = FL.fold FL.list (fmap toSpec keyedFits)
      configuration = GV.configure
        . GV.configuration (GV.ViewStyle [GV.ViewContinuousWidth 800, GV.ViewContinuousHeight 400]) . GV.configuration (GV.PaddingStyle $ GV.PSize 50)
  in GV.toVegaLite [configuration [], GV.layer specs, gvTitle title]

instance ( V.KnownField x
         , Real (V.Snd x)
         , ScatterFitConstraints x y
         , ScatterFitC1 rs x y
         , F.ElemOf (rs V.++ [YError,YFit,YFitError]) x
         , F.ElemOf (rs V.++ [YError,YFit,YFitError]) y
         , F.ElemOf rs x
         , F.ElemOf rs y
         , V.KnownField w
         , rs F.⊆ rs
         ) => Frame2DRegressionScatterFit rs (FR.FrameRegressionResult y 'True '[x] w rs) where
  regressionResultToError = frameRegressionError
  regressionResultToFit fitNameM frr cl =
    let label = fromMaybe "fit" fitNameM
        dof = RE.degreesOfFreedom (FR.regressionResult frr)
        predF = RE.predictFromEstimateAtConfidence dof frr cl
    in FitToPlot label predF

--  scatterPlot title fitNameM frr ci frame =
--    scatterWithFit @x @y @rs @rs (regressionResultToError frr) (regressionResultToFit fitNameM frr ci) Nothing title frame

  scatterPlotSpec fitNameM frr cl frame =
      scatterWithFitSpec @x @y @rs @rs (regressionResultToError frr) (regressionResultToFit fitNameM frr cl) Nothing frame

-- in this case we have y = a(x1) + b(x2) so we plot (y/x1) vs (a + b(x1/x2))
type YOverX1 = "y_over_x1" F.:-> Double
type X2OverX1 = "x2_over_x1" F.:-> Double

instance ( F.ElemOf rs x1
         , F.ElemOf rs x2
         , F.ElemOf rs y
         , V.KnownField y
         , V.KnownField x1
         , V.KnownField x2
         , V.KnownField w
         , Real (V.Snd x1)
         , Real (V.Snd x2)
         , Real (V.Snd y)
         , F.ElemOf [x1,x2] x1
         , F.ElemOf [x1,x2] x2
         , F.ElemOf (rs V.++ [YOverX1, X2OverX1]) x1
         , F.ElemOf (rs V.++ [YOverX1, X2OverX1]) x2
         , F.ElemOf (rs V.++ [YOverX1, X2OverX1]) y
         , F.ElemOf (rs V.++ [YOverX1, X2OverX1]) X2OverX1
         , F.ElemOf (rs V.++ [YOverX1, X2OverX1]) YOverX1
         , rs F.⊆ (rs V.++  [YOverX1, X2OverX1])
         , ScatterFitC1 (rs V.++ [YOverX1, X2OverX1]) YOverX1 X2OverX1
         ) => Frame2DRegressionScatterFit rs (FR.FrameRegressionResult y 'False '[x1,x2] w rs) where
  regressionResultToError frr =
    let x1 = realToFrac . F.rgetField @x1
    in (\r -> frameRegressionError frr r/x1 r)

  regressionResultToFit fitNameM frr cl =
    let label = fromMaybe "fit" fitNameM
        dof = RE.degreesOfFreedom (FR.regressionResult frr)
        predF r =
          let x1 = realToFrac $ F.rgetField @x1 r
              (y,dy) = RE.predictFromEstimateAtConfidence dof frr cl r
          in (y/x1, dy/x1)
    in FitToPlot label predF

  scatterPlotSpec fitNameM frr cl frame =
    let mut :: F.Record rs -> F.Record [YOverX1, X2OverX1]
        mut r =
          let y = F.rgetField @y r
              x1 = F.rgetField @x1 r
              x2 = F.rgetField @x2 r
          in (realToFrac y/realToFrac x1) F.&: (realToFrac x2/realToFrac x1) F.&: V.RNil
        mutData = fmap (FT.mutate mut) frame
        yName = FV.colName @y
        x1Name = FV.colName @x1
        x2Name = FV.colName @x2
        xLabel = x2Name <> "/" <> x1Name
        yLabel = yName <> "/" <> x1Name
    in scatterWithFitSpec @X2OverX1 @YOverX1 @(rs V.++ [YOverX1,X2OverX1]) @rs (regressionResultToError @rs frr) (regressionResultToFit @rs fitNameM frr cl) (Just (xLabel,yLabel)) mutData

{-
frame2DRegressionScatter :: forall x y rs a f. (Foldable f, Frame2DRegressionScatterFit rs x y a)
  => T.Text -> Maybe T.Text -> a -> Double -> f (F.Record rs) -> GV.VegaLite
frame2DRegressionScatter title fitNameM frr ci frame =
-}

--
type YError = "yError" F.:-> Double
type YFit = "yFit" F.:-> Double
type YFitError = "yFitError" F.:->Double

type Error rs = (F.Record rs -> Double)
data FitToPlot rs = FitToPlot { fitLabel :: T.Text, fitFunction :: F.Record rs -> (Double, Double) }

-- | 2D Scatter of Data with calculated error and fit function.
-- | Use TypeApplications to specify x and y columns for scatter and then provide calculated error and fit.
-- | Since both calculations use the record itself as the domain, you can put the error and fit into the record and just use
-- | field selection.  But this allows more flexibility and doesn't require adding things to the input frame record just for the plot.

scatterWithFit :: forall x y rs as f. ( ScatterFitConstraints x y
                                      , as F.⊆ rs
                                      , ScatterFitC1 rs x y
                                      , Foldable f)
               => Error as -> FitToPlot as -> Maybe (T.Text, T.Text) -> Text -> f (F.Record rs) -> GV.VegaLite
scatterWithFit err fit axisLabelsM title frame =
  let configuration = GV.configure
        . GV.configuration (GV.ViewStyle [GV.ViewContinuousWidth 800, GV.ViewContinuousHeight 400]) . GV.configuration (GV.PaddingStyle $ GV.PSize 50)
      swfSpec = GV.specification $ scatterWithFitSpec @x @y @rs @as @f err fit axisLabelsM frame
  in GV.toVegaLite [configuration [], swfSpec, gvTitle title]

scatterWithFitSpec :: forall x y rs as f. ( ScatterFitConstraints x y
                                          , as F.⊆ rs
                                          , ScatterFitC1 rs x y
                                          , Foldable f)
               => Error as -> FitToPlot as -> Maybe (T.Text, T.Text) -> f (F.Record rs) -> GV.VLSpec
scatterWithFitSpec err fit axisLabelsM frame =
  let mut :: F.Record rs -> F.Record '[YError, YFit, YFitError]
      mut r =
        let a = F.rcast r
            (f,fe) = fitFunction fit a
        in err a F.&: f F.&: fe F.&: V.RNil
      vegaDat = FV.recordsToVLData (F.rcast @[x,y,YError,YFit,YFitError] . FT.mutate mut) frame
  in scatterWithFitSpec' @x @y @YError @YFit @YFitError axisLabelsM (fitLabel fit) vegaDat

{-
scatterWithFit' :: forall x y ye fy fye. ( F.ColumnHeaders '[x]
                                         , F.ColumnHeaders '[y]
                                         , F.ColumnHeaders '[ye]
                                         , F.ColumnHeaders '[fy]
                                         , F.ColumnHeaders '[fye])
  => Maybe (T.Text, T.Text) -> Text -> Text -> GV.Data -> GV.VegaLite
scatterWithFit' axisLabelsM title fitLbl dat =
  let configuration = GV.configure
        . GV.configuration (GV.ViewStyle [GV.ViewContinuousWidth 800, GV.ViewContinuousHeight 400]) . GV.configuration (GV.PaddingStyle $ GV.PSize 50)
      swfSpec = GV.specification $ scatterWithFitSpec' @x @y @ye @fy @fye axisLabelsM fitLbl dat
  in GV.toVegaLite [configuration [], swfSpec, gvTitle title]
-}
-- TODO: Add xErrors as well, in scatter and in fit
scatterWithFitSpec' :: forall x y ye fy fye. ( F.ColumnHeaders '[x]
                                             , F.ColumnHeaders '[y]
                                             , F.ColumnHeaders '[ye]
                                             , F.ColumnHeaders '[fy]
                                             , F.ColumnHeaders '[fye])
  => Maybe (T.Text, T.Text) -> Text -> GV.Data -> GV.VLSpec
scatterWithFitSpec' axisLabelsM fitLbl dat =
-- create 4 new cols so we can use rules/areas for errors
  let yLoCalc yName yErrName = "datum." <> yName <> " - (datum." <> yErrName <> ")/2"
      yHiCalc yName yErrName = "datum." <> yName <> " + (datum." <> yErrName <> ")/2"
      calcs = GV.calculateAs (yLoCalc (FV.colName @y) (FV.colName @ye)) "yLo"
              . GV.calculateAs (yHiCalc (FV.colName @y) (FV.colName @ye)) "yHi"
              . GV.calculateAs (yLoCalc (FV.colName @fy) (FV.colName @fye)) "fyLo"
              . GV.calculateAs (yHiCalc (FV.colName @fy) (FV.colName @fye)) "fyHi"
              . GV.calculateAs ("\"" <> fitLbl <> "\"")  "fitLabel"
      xLabel = maybe (FV.colName @x) fst axisLabelsM
      yLabel = maybe (FV.colName @y) snd axisLabelsM
      xEnc = GV.position GV.X [FV.pName @x, GV.PmType GV.Quantitative, GV.PAxis [GV.AxTitle xLabel]]
      yEnc = GV.position GV.Y [FV.pName @y, GV.PmType GV.Quantitative, GV.PAxis [GV.AxTitle yLabel]]
      yLEnc = GV.position GV.Y [GV.PName "yLo", GV.PmType GV.Quantitative]
      yHEnc = GV.position GV.Y2 [GV.PName "yHi", GV.PmType GV.Quantitative]
--      yErrorEnc = GV.position GV.YError [FV.pName @ye, GV.PmType GV.Quantitative]
      yFitEnc t = GV.position GV.Y [FV.pName @fy, GV.PmType GV.Quantitative, GV.PAxis [GV.AxTitle t]]
      yFitLEnc = GV.position GV.Y [GV.PName "fyLo", GV.PmType GV.Quantitative]
      yFitHEnc = GV.position GV.Y2 [GV.PName "fyHi", GV.PmType GV.Quantitative]
      colorEnc = GV.color [GV.MName "fitLabel", GV.MmType GV.Nominal]
  --    yFitErrorEnc = GV.position GV.YError [FV.pName @fye, GV.PmType GV.Quantitative]
      scatterEnc = xEnc . yEnc
      scatterBarEnc = xEnc . yLEnc . yHEnc
      fitEnc t = xEnc . yFitEnc t . colorEnc
      fitBandEnc = xEnc . yFitLEnc . yFitHEnc . colorEnc

      selectScalesS = GV.select ("scalesS" <> fitLbl) GV.Interval [GV.BindScales]
      selectScalesSE = GV.select("scalesSE" <> fitLbl) GV.Interval [GV.BindScales]
      selectScalesF = GV.select ("scalesF" <> fitLbl) GV.Interval [GV.BindScales]
      selectScalesFE = GV.select ("scalesFE" <> fitLbl) GV.Interval [GV.BindScales]
      scatterSpec = GV.asSpec
        [
          (GV.encoding . scatterEnc) []
        , GV.mark GV.Point []
        , (GV.selection . selectScalesS) []
        ]
      scatterBarSpec = GV.asSpec
        [
          (GV.encoding . scatterBarEnc) []
        , GV.mark GV.Rule []
        , (GV.selection . selectScalesSE) []
        ]
      fitLineSpec = GV.asSpec
        [
          (GV.encoding . fitEnc fitLbl) []
        , GV.mark GV.Line []
        , (GV.selection . selectScalesF) []
        ]
      fitBandSpec = GV.asSpec
        [
          (GV.encoding . fitBandEnc) []
        , GV.mark GV.Area [GV.MOpacity 0.5, GV.MFillOpacity 0.3]
        , (GV.selection . selectScalesFE) []
        ]
      layers = GV.layer [ scatterSpec, scatterBarSpec, fitLineSpec, fitBandSpec]
      configuration = GV.configure
        . GV.configuration (GV.ViewStyle [GV.ViewContinuousWidth 800, GV.ViewContinuousHeight 400]) . GV.configuration (GV.PaddingStyle $ GV.PSize 50)
      spec =
        GV.asSpec
        [
          (GV.transform . calcs) []
        , layers
        , dat
        , configuration []
        ]
  in spec

max_stars_repo_name: teto/Frames-utils

max_issues_repo_name: teto/Frames-utils

max_forks_repo_name: teto/Frames-utils


text:
module Main where

import Network.Layer
import Network.Network
import Network.Neuron
import Network.Trainer

import Numeric.LinearAlgebra
import System.Random

main :: IO ()
main = do

  g <- newStdGen
  let l   = LayerDefinition sigmoidNeuron 2 connectFully
      l'  = LayerDefinition sigmoidNeuron 2 connectFully
      l'' = LayerDefinition sigmoidNeuron 1 connectFully

  let n = createNetwork normals g [l, l', l'']

  let t = BackpropTrainer 3 quadraticCost quadraticCost'

  let dat = [(fromList [0, 1], fromList [1]),
             (fromList [1, 1], fromList [0]),
             (fromList [1, 0], fromList [1]),
             (fromList [0, 0], fromList [0])]

  let n' = trainNTimes n t online dat 10000

  putStrLn "==> XOR predictions: "
  putStr "(0,0): "
  print $ predict (fromList [0, 0]) n'
  putStr "(1,0): "
  print $ predict (fromList [1, 0]) n'
  putStr "(0,1): "
  print $ predict (fromList [0, 1]) n'
  putStr "(1,1): "
  print $ predict (fromList [1, 1]) n'

  putStrLn "\n==> Trained network: "
  print $ map layerToShowable $ layers n'

{-
  saveNetwork "xor.ann" n'

  putStrLn "==> Network saved and reloaded: "
  n'' <- loadNetwork "xor.ann" [l, l', l'']

  print $ predict (fromList [0, 0]) n''
  print $ predict (fromList [1, 0]) n''
  print $ predict (fromList [0, 1]) n''
  print $ predict (fromList [1, 1]) n''
-}

max_stars_repo_name: IMOKURI/nn-sample

max_issues_repo_name: IMOKURI/nn-sample

max_forks_repo_name: IMOKURI/nn-sample


text:
{-# LANGUAGE UnicodeSyntax, DataKinds, TypeOperators, KindSignatures,
             TypeInType, GADTs, MultiParamTypeClasses, FunctionalDependencies,
             TypeFamilies, AllowAmbiguousTypes, FlexibleInstances,
             UndecidableInstances, InstanceSigs, TypeApplications,
             ScopedTypeVariables, EmptyCase, FlexibleContexts,
             RankNTypes, LambdaCase
#-}

-- Implementation of density matrix interpretation of quantum computation
module LinTrans where

import Prelim

import qualified Data.Complex as C
import Data.Complex (Complex(..))
import Control.Monad.State.Lazy
import Data.Singletons
import Data.Maybe (fromJust)
--import Data.Tuple hiding (swap)
import Unsafe.Coerce
import Data.Constraint (Dict(..))
import Data.List ( (\\) )
import qualified Debug.Trace as Tr

-- I don't like the show instance for ℂ
newtype ℂ = ℂ (Complex Double)

liftC :: (Complex Double -> Complex Double) -> ℂ -> ℂ
liftC f (ℂ c) = ℂ $ f c

i :: ℂ
i = ℂ $ 0 :+ 1
conjugate = liftC C.conjugate
instance Num ℂ where
  ℂ m + ℂ n = ℂ $ m + n
  ℂ m - ℂ n = ℂ $ m - n
  ℂ m * ℂ n = ℂ $ m * n
  abs = liftC abs
  signum = liftC signum
  fromInteger = ℂ . fromInteger
instance Fractional ℂ where
  fromRational = ℂ . fromRational
  ℂ m / ℂ n    = ℂ $ m / n
instance Floating ℂ where
  pi    = ℂ $ pi
  exp   = liftC exp
  log   = liftC log
  sin   = liftC sin
  cos   = liftC cos
  asin  = liftC asin
  acos  = liftC acos
  atan  = liftC atan
  sinh  = liftC sinh
  cosh  = liftC cosh
  asinh = liftC asinh
  acosh = liftC acosh
  atanh = liftC atanh
instance Show ℂ where
  show (ℂ (α :+ β)) =   if β == 0 then show α
                        else if α == 0 then show β ++ "i"
                        else show α ++ " + " ++ show β ++ "i"

-- A matrix of dimension m × n is a function from
-- pairs of nats bounded by m and n respectively.
type Matrix m n = (BNat m,BNat n) -> ℂ
type Squared n = Matrix (Two `RaiseTo` n) (Two `RaiseTo` n)

matrix :: forall m n. (SingI m, SingI n) => [ℂ] -> Matrix m n
matrix ls (i,j) = ls !! (n * fromIntegral i + fromIntegral j)
  where
    n = toInt (sing :: Sing n)

--------------
-- Printing --
--------------
rows :: forall m n. (SingI m, SingI n) => Matrix m n -> [[ℂ]]
rows mat = fmap f $ allBNat (sing :: Sing m)
  where
    f i = fmap (\j -> mat (i,j)) $ allBNat (sing :: Sing n)

showRows :: (SingI m, SingI n) => Matrix m n -> [String]
showRows mat = show <$> rows mat

instance (SingI m, SingI n) => Show (Matrix m n) where
  show mat = unlines $ showRows mat

---------------------------------------
-- Primitive matrices and operatoins --
---------------------------------------

ident :: Matrix m m
ident (i,j) = if i==j then 1 else 0

transpose :: Matrix m n -> Matrix n m
transpose mat (i,j) = conjugate $ mat (j,i)

{-
-- if i : BNat m then raise m n = i : BNat (S m)
-- so pruneRows mat is just the first m rows of mat
pruneRows :: Sing m -> Matrix (S m) n -> Matrix m n
pruneRows m mat (i,j) = mat (raise m i,j)

pruneCols :: Sing n -> Matrix m (S n) -> Matrix m n
pruneCols n mat (i,j) = mat (i,raise n j)
-}

sumTo :: Sing n -> BNat n -> (BNat n -> ℂ) -> ℂ
sumTo _ (BNat SZ)     _ = 0
sumTo n (BNat (SS m)) f = f m' + sumTo n m' f
  where
    m' = case succLtTrans m n of Dict -> BNat m -- m < S m < n

dot :: Sing n -> Matrix (S Z) n -> Matrix n (S Z) -> ℂ
dot n v1 v2 = sumTo n (maxBNat n) (\x -> v1 (0,x) * v2(x,0))

row :: BNat m -> Matrix m n -> Matrix (S Z) n
row i mat (_,j) = mat (i,j)

col :: BNat n -> Matrix m n -> Matrix m (S Z)
col j mat (i,_) = mat (i,j)

plus :: Matrix m n -> Matrix m n -> Matrix m n
plus mat1 mat2 (i,j) = mat1(i,j) + mat2(i,j)

mult :: forall m n p. SingI n => Matrix m n -> Matrix n p -> Matrix m p
mult mat1 mat2 (i,j) =
--  Tr.trace ("Composing for intermediate " ++ show (sing :: Sing n)) $
    dot sing (row i mat1) (col j mat2)

kron :: forall m1 m2 n1 n2. (SingI m1, SingI n1, SingI m2, SingI n2)
     => Matrix m1 n1 -> Matrix m2 n2 -> Matrix (m1 `Times` m2) (n1 `Times` n2)
kron mat1 mat2 (BNat i,BNat j) =
--    Tr.trace ("Computing ⊗ for " ++ show (m1, n1) ++ "×" ++ show (m2, n2)) $
    mat1 (bNat i1, bNat j1) * mat2 (bNat i2,bNat j2)
  where
    i1 = i `divSNat` m2
    i2 = i `modSNat` m2
    j1 = j `divSNat` n2
    j2 = j `modSNat` n2
    m2 = (sing :: Sing m2)
    n2 = (sing :: Sing n2)

trace :: forall m. SingI m => Matrix m m -> ℂ
trace mat = foldr f 0 ls
  where
    ls = allBNat (sing :: Sing m)
    f i x = mat (i,i) + x

----------------------
-- Quantum Matrices --
----------------------

square :: forall n. SingI n => [ℂ] -> Squared n
square ls = withSingI (two `raiseToSNat` (sing :: Sing n)) $ matrix ls

density0 :: Matrix Two Two
density0 = matrix [1,0
                  ,0,0]

density1 :: Matrix Two Two
density1 = matrix [0,0
                  ,0,1]

hadamard = square @One $ [1/sqrt 2, 1/sqrt 2
                         ,1/sqrt 2, -1/sqrt 2]

pauliX = square @One $ [0,1,1,0]
pauliY = square @One $ [0,-i,i,0]
pauliZ = square @One $ [1,0,0,-1]

cnotD  = square @Two $ [1,0,0,0
                       ,0,1,0,0
                       ,0,0,0,1
                       ,0,0,1,0]

newD :: Bool -> Squared One
newD True  = density1
newD False = density0

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

----------------------
-- Density Matrices --
----------------------

data Density where
  Density :: Sing n -> Squared n -> Density
instance Show Density where
  show (Density n m) = withSingI (two `raiseToSNat` n) $ show m

size :: forall m. SingI m => Matrix m m -> Int
size _ = toInt (sing :: Sing m)
logsizeD :: Density -> Int
logsizeD (Density n mat) = toInt n
sizeD :: Density -> Int
sizeD ρ = 2 ^ logsizeD ρ

transposeD :: Density -> Density
transposeD (Density m ρ) = Density m $ transpose ρ

multD :: Density -> Density -> Density
multD = undefined

kronD :: Density -> Density -> Density
kronD (Density m mat1) (Density n mat2) =
      withSingI (raiseToSNat two n) $
      withSingI (raiseToSNat two m) $
      withSingI (raiseToSNat two (plusSNat m n)) $
      case raiseToPlus two m n of {Dict ->
        Density (m `plusSNat` n) $ mat1 `kron` mat2
      }

identD :: Int -> Density
identD n = case fromIntegral @_ @(SomeSing Nat) n of SomeSing n' -> Density n' ident

nil :: Int -> Density
nil n = case fromIntegral n of
          SomeSing n' -> withSingI (two `raiseToSNat` n') $
                         Density n' . matrix $ fromIntegral <$> repeat 0

------------------
-- Permutations --
------------------

-- encode and decode modulo a key k
encode :: Int -> [Int] -> Int
encode k [] = 0
encode k (i : ls) = i + k * encode k ls

decode :: Int -> Int -> [Int]
decode k i | i == 0 = []
decode k i | otherwise = i `mod` k : decode k (i `div` k)

fromAssocList :: [(Int,Int)] -> Int -> Int
fromAssocList []           = id
fromAssocList ((a,b) : ls) = assocFun a b . fromAssocList ls
  where
    assocFun a b i = if i == a then b else i

-- Check if two lists are equal up to the permutation specified, i.e.
-- if permuting ls1 by f is equal to ls2
isPermutation :: Eq a => (Int -> Int) -> [a] -> [a] -> Bool
isPermutation f ls1 ls2 = all (\i -> ls1 !! i == ls2 !! f i) [0..len]
  where
    len = length ls1

swapFun :: forall n. SingI n => (Int -> Int) -> Squared n
swapFun f (i,j) = if isPermutation f (decode n' i') (decode n' j') then 1 else 0
  where
    n' = toInt (sing :: Sing n)
    i' = toInt i
    j' = toInt j

-- swap takes a list of qubit/bit variables and produces a matrix that permutes
-- those qubits.
-- i.e. swap ls |φ_0 ⋯ φ_n⟩ = |φ_f(0) ⋯ ρ_f(n)⟩
-- where f = fromAssocList ls
swap :: forall n. SingI n => [Int] -> Squared n
swap ls = swapFun @n (fromAssocList $ zip [0..] ls)

-------------------
-- Density Monad --
-------------------

-- A density monad is a nondeterminism state monad:
-- Density -> [Density]
-- An element op of type DensityMonad corresponds to the superoperator
-- \ρ -> ∑ (op ρ)
type DensityMonad = StateT Density []

newM :: Bool -> DensityMonad Int
newM b = do
    ρ <- get
    put $ ρ `kronD` (Density one $ newD b)
    return $ sizeD ρ

-- runQ applies the superoperator to the identity density matrix of size 1.
runQ :: DensityMonad a -> [(a,Density)]
runQ m = runStateT m (Density SZ ident)

-- getDensity combines the result of runQ into a single density matrix
getDensity :: DensityMonad a -> Density
getDensity m = foldr f (nil n) ls
  where
    ls = runQ m
    n = logsizeD . snd $ head ls
    f (_,ρ) ρ0 = ρ + ρ0

-----------------
-- Application --
-----------------

superM :: forall m n. Matrix (Two `RaiseTo` m) (Two `RaiseTo` n) -> DensityMonad ()
superM mat = do Density _ ρ ← get
                put . Density _ $ mat `mult` ρ `mult` transpose mat
  where
    matD = Density (sing :: Sing m) mat

applyUnitaryM :: forall m. SingI m => Squared m -> [Int] -> DensityMonad ()
applyUnitaryM u ls = do superM @m σ -- moves the relevant qubits to the front of the density matrix
                        superM @m u -- applies operation
                        superM @m $ transpose σ -- moves the qubits back
  where
    σ = swap @m ls

applyMatrixM :: Matrix (Two `RaiseTo` m) (Two `RaiseTo` n)
             -> [Int] -> DensityMonad ()
applyMatrixM mat ls = undefined

measM :: Int -> DensityMonad Bool
measM i = do
    ρ <- get
    (b,ρ') <- lift [ (False, applyMatrixD @One density0 [i] ρ)
                   , (True,  applyMatrixD @One density1 [i] ρ) ]
    put ρ'
    return b

------------------
-- Num instance --
------------------

instance (SingI m, SingI n) => Num (Matrix m n) where
  (+) mat1 mat2 (i,j) = mat1(i,j) + mat2(i,j)
  (-) mat1 mat2 (i,j) = mat1(i,j) - mat2(i,j)
  (*) mat1 mat2 (i,j) = mat1(i,j) * mat2(i,j)
  abs mat (i,j) = abs $ mat(i,j)
  signum mat = undefined -- trace mat
  fromInteger n = matrix . repeat $ fromIntegral n
instance Num Density where
  Density n1 mat1 + Density n2 mat2 =
    withSingI (two `raiseToSNat` n1) $
    withSingI (two `raiseToSNat` n2) $
    case eqSNat n1 n2 of
      Left Dict -> Density n1 $ mat1 + mat2
      Right _   -> error $ "Cannot add mismatched matrices "
                            ++ show mat1 ++ " and " ++ show mat2
  _ - _ = undefined
  _ * _ = undefined
  abs _ = undefined
  signum = undefined
  fromInteger n = undefined

instance Show (DensityMonad a) where
  show = show . getDensity

-- Debugging help
m0 = square @Two [1,0,0,0
                 ,0,1,0,0
                 ,0,0,0,1
                 ,0,0,1,0]
rho = Density three $ \case
  (0,0) -> -0.5
  (0,2) ->  0.5
  (3,0) -> -0.5
  (3,2) ->  0.5
  (_,_) -> 0

{-
-- produces an association list
listToAssoc :: forall n. SingI n => [Int] -> [(BNat n, BNat n)]
listToAssoc ls = (fromIntegral <$> [0..n]) `zip` (fromIntegral <$> pad ls)
  where
    n = maxBNat (sing :: Sing n)

pad :: [Int] -> [Int]
pad ls = ls ++ ([0..] \\ ls)

mkPermutation :: [(BNat n,BNat n)] -> BNat n -> BNat n
mkPermutation ls i =
  case lookup i ls of
    Just i' -> i'
    Nothing -> error $ "Permutation " ++ show ls
                        ++ " not well defined at index " ++ show i

permute :: forall m. SingI m => [(BNat m, BNat m)] -> Matrix m m -> Matrix m m
permute ls mat (i,j) = mat (mkPermutation ls i, mkPermutation ls j)

expand :: forall m n. (SingI m, SingI n)
       => [(BNat (m `Times` n), BNat (m `Times` n))]
       -> Matrix m m -> Matrix (m `Times` n) (m `Times` n)
expand ls mat = Tr.trace ("Expanding..." ++ show mat) $
                withSingI (timesSNat (sing :: Sing m) (sing :: Sing n)) $
                  permute ls $ mat `kron` ident @n

applyMatrix :: forall m p. (SingI m, SingI p, (p < m) ~ 'False)
            => Squared m -> [Int] -> Squared p -> Squared p
applyMatrix mat ls ρ =
  Tr.trace ("Applying " ++ show ls) $
  case raiseToPlus two m minus of {Dict ->
  case plusMinus m p of {Dict ->
  withSingI (raiseToSNat two m) $
  withSingI (raiseToSNat two minus) $
  withSingI (raiseToSNat two p) $
    let mat1 = expand @(Two `RaiseTo` m) @(Two `RaiseTo` (p `Minus` m))
                       (swap <$> ls') (transpose mat)
        mat2 = expand @(Two `RaiseTo` m) @(Two `RaiseTo` (p `Minus` m)) ls' mat
        res  = Tr.trace "Middle of application" $
               mat1 `mult` ρ `mult` mat2
    in mat1 `seq` mat2 `seq` res `seq` Tr.trace "End of application" res
  }}
  where
    ls' :: [(BNat (Two `RaiseTo` p), BNat (Two `RaiseTo` p))]
    ls' = withSingI (raiseToSNat two p) $ listToAssoc ls
    m = (sing :: Sing m)
    p = (sing :: Sing p)
    minus :: Sing (p `Minus` m)
    minus = minusSNat p m

applyMatrixD :: forall m. SingI m => Squared m -> [Int] -> Density -> Density
applyMatrixD mat perm (Density (p :: Sing p) ρ) =
  Tr.trace ("Applying " ++ mat' ++ " to qubits " ++ show perm ++ " of \n" ++ show (Density p ρ)) $
  case compareSNat p (sing :: Sing m) of
    -- p < m
    Left Dict -> error "Cannot apply matrix to state: matrix too large"
    -- m <= p
    Right Dict ->
      let res = withSingI p $ Density p $ applyMatrix @m @p mat perm ρ
      in  res `seq` Tr.trace "Done!" res
  where
    mat' = withSingI (two `raiseToSNat` (sing :: Sing m)) $ show mat

applyUnitaryM :: forall m. SingI m => Squared m -> [Int] -> DensityMonad ()
applyUnitaryM mat ls = get >>= \ρ -> put $ applyMatrixD @m mat ls ρ

measM :: Int -> DensityMonad Bool
measM i = do
    ρ <- get
    (b,ρ') <- lift [ (False, applyMatrixD @One density0 [i] ρ)
                   , (True,  applyMatrixD @One density1 [i] ρ) ]
    put ρ'
    return b

-}

max_stars_repo_name: jpaykin/LNLHaskell

max_issues_repo_name: jpaykin/LNLHaskell

max_forks_repo_name: jpaykin/LNLHaskell


text:
{- |
Module          : $Header$
Description     : This is the parser file for the Scheme compiler.
Copyright       : (c) Michael Buchel
License         : BSD3

Maintainer      : abuchel2@uwo.ca
Stability       : experimental
Portability     : portable

This is the parser, where I am using Parsec for simplicity.
| -}
module Parser where

import Abstract
import Control.Monad
import Control.Monad.Except
import Data.Array
import Data.Complex
import Data.Ratio
import Error
import Evaluator
import Numbers
import Text.ParserCombinators.Parsec

-- | Symbol function which identifies all the symbols allowed in Scheme.
symbol :: Parser Char
symbol = oneOf "!$%&|*+-/:<=>?@^_~"

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

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

-- | Parses an atom.
parseAtom :: Parser LispVal
parseAtom = do
	first <- letter <|> symbol
	rest <- many (letter <|> digit <|> symbol)
	let atom = first:rest
	return $ Atom atom

-- | Parses a boolean value.
parseBool :: Parser LispVal
parseBool = do
	char '#'
	(char 't' >> return (Bool True)) <|> (char 'f' >> return (Bool False))

-- | Parses decimals without the identifier string.
parseDecimal1 :: Parser LispVal
parseDecimal1 = many1 digit >>= (return . Number . read)

-- | Parses decimals with the identifier string.
parseDecimal2 :: Parser LispVal
parseDecimal2 = do
	try $ string "#d"
	many1 digit >>= (return . Number . read)

-- | Parses hexadecimals.
parseHex :: Parser LispVal
parseHex = do
	try $ string "#x"
	many1 hexDigit >>= (return . Number . hex2dig)

-- | Parses octal numbers.
parseOct :: Parser LispVal
parseOct = do
	try $ string "#o"
	many1 octDigit >>= (return . Number . oct2dig)

-- | Parses binary numbers.
parseBin :: Parser LispVal
parseBin = do
	try $ string "#b"
	many1 (oneOf "10") >>= (return . Number . bin2dig)

-- | Parses a number.
parseNumber :: Parser LispVal
parseNumber = parseDecimal1 <|> parseDecimal2 <|> parseHex <|> parseOct <|> parseBin

-- | Parses a floating point number.
parseFloat :: Parser LispVal
parseFloat = do
	x <- many1 digit
	char '.'
	y <- many1 digit
	return $ Float $ getFloat x y

-- | Parses a ratio value.
parseRatio :: Parser LispVal
parseRatio = do
	x <- many1 digit
	spaces >> char '/' >> spaces
	y <- many1 digit
	return $ Ratio ((read x) % (read y))

-- | Parses a complex number.
parseComplex :: Parser LispVal
parseComplex = do
	x <- (try parseFloat <|> parseDecimal1)
	spaces >> char '+' >> spaces
	y <- (try parseFloat <|> parseDecimal1)
	oneOf "ij"
	return $ Complex (toDouble x :+ toDouble y)

-- | Helper list parser.
parseListHelp :: Parser LispVal
parseListHelp = do
	char '.' >> (skipMany1 space)
	x <- parseExpr
	spaces >> char ')'
	return x

-- | Parses an individual list.
parseList :: Parser LispVal
parseList = do
	char '(' >> spaces
	headList <- parseExpr `sepEndBy` (skipMany1 space)
	(parseListHelp >>= return . (DottedList headList)) <|> (spaces >> char ')' >> return (List headList))

-- | For Scheme syntax sugar.
parseQuoted :: Parser LispVal
parseQuoted = do
	char '\''
	x <- parseExpr
	return $ List [Atom "quote", x]

-- | For Scheme syntax sugar.
parseQuasiQuoted :: Parser LispVal
parseQuasiQuoted = do
	char '`'
	x <- parseExpr
	return $ List [Atom "quasiquote", x]

-- | For Scheme syntax sugar.
parseUnQuoted :: Parser LispVal
parseUnQuoted = do
	char ','
	x <- parseExpr
	return $ List [Atom "unquote", x]

-- | Parses a character value.
parseChar :: Parser LispVal
parseChar = do
	try $ string "#\\"
	value <- try (string "newline" <|> string "space") <|> do {x <- anyChar; notFollowedBy alphaNum; return [x]}
	return $ Character $ case value of
				  "space" -> ' '
				  "newline" -> '\n'
				  otherwise -> (value !! 0)

-- | Parses a vector.
parseVector :: Parser LispVal
parseVector = do
	arrayVals <- sepBy parseExpr spaces
	return $ Vector (listArray (0, (length arrayVals - 1)) arrayVals)

-- | Helper function for parsing vectors.
parseVectors :: Parser LispVal
parseVectors = do
	string "#(" >> spaces
	x <- parseVector
	spaces >> char ')'
	return x

-- | Helper parsing function which tries to parse the expression given to it.
parseExpr :: Parser LispVal
parseExpr = parseAtom
	<|> parseString
	<|> parseQuoted
	<|> parseQuasiQuoted
	<|> parseUnQuoted
	<|> try parseComplex
	<|> try parseFloat
	<|> try parseRatio
	<|> try parseNumber
	<|> try parseBool
	<|> try parseChar
	<|> try parseVectors
	<|> parseList

-- | A function to read the expression from a string.
readExpr :: String -- ^ Input string
	-> ThrowsError LispVal
readExpr input = case parse (spaces >> parseExpr) "parser" input of
		      Left err -> throwError $ Parser err
		      Right val -> return val

-- | Find the end of an expression.
findExpr :: String -- ^ Input string.
	-> Integer -- ^ Previous integer.
	-> Integer -- ^ Nesting of expression.
	-> Integer
findExpr [] x _ = x
findExpr (x : xs) n y = case x of
			     '(' -> findExpr xs (n + 1) (y + 1)
			     ')' -> if (y - 1) == 0 then (n + 1) else findExpr xs (n + 1) (y - 1)
			     otherwise -> findExpr xs (n + 1) y

max_stars_repo_name: mbuchel/compiler

max_issues_repo_name: mbuchel/compiler

max_forks_repo_name: mbuchel/compiler


text:
{- |
Raw, unnormalized versions of the transforms in @fftw@.

Note that the forwards and backwards transforms of this module are not actually
inverses.  For example, @run idft (run dft v) /= v@ in general.

For more information on the individual transforms, see
<http://www.fftw.org/fftw3_doc/What-FFTW-Really-Computes.html>.
-}
module Numeric.FFT.Vector.Unnormalized(
                    -- * Creating and executing 'Plan's
                    run,
                    plan,
                    execute,
                    -- * Complex-to-complex transforms
                    dft,
                    idft,
                    -- * Real-to-complex transforms
                    dftR2C,
                    dftC2R,
                    -- * Real-to-real transforms
                    -- $dct_size
                    -- ** Discrete cosine transforms
                    dct1,
                    dct2,
                    dct3,
                    dct4,
                    -- ** Discrete sine transforms
                    dst1,
                    dst2,
                    dst3,
                    dst4,
                    ) where

import Numeric.FFT.Vector.Base
import Foreign
import Foreign.C
import Data.Complex

#include <fftw3.h>

-- | Whether the complex fft is forwards or backwards.
type CDirection = CInt

-- | The type of the cosine or sine transform.
type CKind = (#type fftw_r2r_kind)

foreign import ccall unsafe fftw_plan_dft_1d
    :: CInt -> Ptr (Complex Double) -> Ptr (Complex Double) -> CDirection
        -> CFlags -> IO (Ptr CPlan)

foreign import ccall unsafe fftw_plan_dft_r2c_1d
    :: CInt -> Ptr Double -> Ptr (Complex Double) -> CFlags -> IO (Ptr CPlan)

foreign import ccall unsafe fftw_plan_dft_c2r_1d
    :: CInt -> Ptr (Complex Double) -> Ptr Double -> CFlags -> IO (Ptr CPlan)

foreign import ccall unsafe fftw_plan_r2r_1d
    :: CInt -> Ptr Double -> Ptr Double -> CKind -> CFlags -> IO (Ptr CPlan)

dft1D :: CDirection -> Transform (Complex Double) (Complex Double)
dft1D d = Transform {
            inputSize = id,
            outputSize = id,
            creationSizeFromInput = id,
            makePlan = \n a b -> withPlanner . fftw_plan_dft_1d n a b d,
            normalization = const id
            }

-- | A forward discrete Fourier transform.  The output and input sizes are the same (@n@).
--
-- @y_k = sum_(j=0)^(n-1) x_j e^(-2pi i j k/n)@
dft :: Transform (Complex Double) (Complex Double)
dft = dft1D (#const FFTW_FORWARD)

-- | A backward discrete Fourier transform.  The output and input sizes are the same (@n@).
--
-- @y_k = sum_(j=0)^(n-1) x_j e^(2pi i j k/n)@
idft :: Transform (Complex Double) (Complex Double)
idft = dft1D (#const FFTW_BACKWARD)

-- | A forward discrete Fourier transform with real data.  If the input size is @n@,
-- the output size will be @n \`div\` 2 + 1@.
dftR2C :: Transform Double (Complex Double)
dftR2C = Transform {
            inputSize = id,
            outputSize = \n -> n `div` 2 + 1,
            creationSizeFromInput = id,
            makePlan = \n a b -> withPlanner . fftw_plan_dft_r2c_1d n a b,
            normalization = const id
        }

-- | A backward discrete Fourier transform which produces real data.
--
-- This 'Transform' behaves differently than the others:
--
--  - Calling @plan dftC2R n@ creates a 'Plan' whose /output/ size is @n@, and whose
--    /input/ size is @n \`div\` 2 + 1@.
--
--  - If @length v == n@, then @length (run dftC2R v) == 2*(n-1)@.
dftC2R :: Transform (Complex Double) Double
dftC2R = Transform {
            inputSize = \n -> n `div` 2 + 1,
            outputSize = id,
            creationSizeFromInput = \n -> 2 * (n-1),
            makePlan = \n a b -> withPlanner . fftw_plan_dft_c2r_1d n a b,
            normalization = const id
        }

r2rTransform :: CKind -> Transform Double Double
r2rTransform kind = Transform {
                    inputSize = id,
                    outputSize = id,
                    creationSizeFromInput = id,
                    makePlan = \n a b -> withPlanner . fftw_plan_r2r_1d n a b kind,
                    normalization = const id
                }

-- $dct_size
-- The real-even (DCT) and real-odd (DST) transforms.  The input and output sizes
-- are the same (@n@).

-- | A type-1 discrete cosine transform.
--
-- @y_k = x_0 + (-1)^k x_(n-1) + 2 sum_(j=1)^(n-2) x_j cos(pi j k\/(n-1))@
dct1 :: Transform Double Double
dct1 = r2rTransform (#const  FFTW_REDFT00)

-- | A type-2 discrete cosine transform.
--
-- @y_k = 2 sum_(j=0)^(n-1) x_j cos(pi(j+1\/2)k\/n)@
dct2 :: Transform Double Double
dct2 = r2rTransform (#const  FFTW_REDFT10)

-- | A type-3 discrete cosine transform.
--
-- @y_k = x_0 + 2 sum_(j=1)^(n-1) x_j cos(pi j(k+1\/2)\/n)@
dct3 :: Transform Double Double
dct3 = r2rTransform (#const  FFTW_REDFT01)

-- | A type-4 discrete cosine transform.
--
-- @y_k = 2 sum_(j=0)^(n-1) x_j cos(pi(j+1\/2)(k+1\/2)\/n)@
dct4 :: Transform Double Double
dct4 = r2rTransform (#const  FFTW_REDFT11)

-- | A type-1 discrete sine transform.
--
-- @y_k = 2 sum_(j=0)^(n-1) x_j sin(pi(j+1)(k+1)\/(n+1))@
dst1 :: Transform Double Double
dst1 = r2rTransform (#const  FFTW_RODFT00)

-- | A type-2 discrete sine transform.
--
-- @y_k = 2 sum_(j=0)^(n-1) x_j sin(pi(j+1\/2)(k+1)\/n)@
dst2 :: Transform Double Double
dst2 = r2rTransform (#const  FFTW_RODFT10)

-- | A type-3 discrete sine transform.
--
-- @y_k = (-1)^k x_(n-1) + 2 sum_(j=0)^(n-2) x_j sin(pi(j+1)(k+1\/2)/n)@
dst3 :: Transform Double Double
dst3 = r2rTransform (#const  FFTW_RODFT01)

-- | A type-4 discrete sine transform.
--
-- @y_k = sum_(j=0)^(n-1) x_j sin(pi(j+1\/2)(k+1\/2)\/n)@
dst4 :: Transform Double Double
dst4 = r2rTransform (#const FFTW_RODFT11)

ext: hsc

max_stars_repo_name: TravisWhitaker/vector-fftw

max_issues_repo_name: TravisWhitaker/vector-fftw

max_forks_repo_name: TravisWhitaker/vector-fftw


text:
module Parser
  (readExpr) where

import Control.Monad
import Text.ParserCombinators.Parsec hiding (spaces)
import Data.Ratio
import Numeric
import Data.Complex

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

readExpr :: String -> String
readExpr input = case parse parseExpr "lisp" input of
                   Left err -> "No match. " ++ show err
                   Right val -> "Found Val " ++ show val

spaces :: Parser()
spaces = skipMany1 space

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

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

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

parseCharacter :: Parser LispVal
parseCharacter = do
  string "#\\"
  x <- many $ anyChar
  return $ case x of
             "space" -> Character ' '
             "newline" -> Character '\n'
             _ -> Character (x !! 0)

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

parseAtom :: Parser LispVal
parseAtom = do first <- letter <|> symbol
               rest <- many (letter <|> digit <|> symbol)
               return $ Atom $ first:rest

parseList :: Parser LispVal
parseList = liftM List $ sepBy parseExpr spaces

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

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

parseNumber :: Parser LispVal
parseNumber = parseDecimal <|> parseHex <|> parseOct <|> parseBin

parseDecimal :: Parser LispVal
parseDecimal = many1 digit >>= (return . Number . read)

parseDecimalAlt :: Parser LispVal
parseDecimalAlt = do try $ string "#d"
                     x <- many1 digit
                     return $ Number (read x)

parseHex :: Parser LispVal
parseHex = do try $ string "#x"
              x <- many1 hexDigit
              return $ Number (hex2dig x)

parseOct :: Parser LispVal
parseOct = do try $ string "#o"
              x <- many1 octDigit
              return $ Number (oct2dig x)

parseBin :: Parser LispVal
parseBin = do try $ string "#b"
              x <- many1 (oneOf "10")
              return $ Number (bin2dig x)

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

parseFloat :: Parser LispVal
parseFloat = do x <- many1 digit
                char '.'
                y <- many1 digit
                return $ Float (fst.head $ readFloat (x ++ "." ++ y))

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

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

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

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

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

parseExpr :: Parser LispVal
parseExpr = parseAtom
         <|> parseString
         <|> try parseComplex
         <|> try parseRatio
         <|> try parseFloat
         <|> try parseNumber
         <|> try parseBool
         <|> try parseCharacter
         <|> parseQuoted
         <|> parseBackquote
         <|> parseCommaList
         <|> do char '('
                x <- try parseList <|> parseDottedList
                char ')'
                return x


text:
{-# LANGUAGE FlexibleContexts #-}
-- |
-- Module    : Statistics.Sample
-- Copyright : (c) 2008 Don Stewart, 2009 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- Commonly used sample statistics, also known as descriptive
-- statistics.

module Statistics.Sample
    (
    -- * Statistics of location
      mean

    -- ** Two-pass functions (numerically robust)
    -- $robust
    , variance
    , varianceUnbiased
    , stdDev

    -- * References
    -- $references
    ) where

import Statistics.Sample.Internal (robustSumVar, sum)
import qualified Data.Vector as V
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed as U

-- Operator ^ will be overriden
import Prelude hiding ((^), sum)

-- | /O(n)/ Arithmetic mean.  This uses Kahan-Babuška-Neumaier
-- summation, so is more accurate than 'welfordMean' unless the input
-- values are very large.
mean :: (G.Vector v Double) => v Double -> Double
mean xs = sum xs / fromIntegral (G.length xs)
{-# SPECIALIZE mean :: U.Vector Double -> Double #-}
{-# SPECIALIZE mean :: V.Vector Double -> Double #-}

-- $variance
--
-- The variance—and hence the standard deviation—of a
-- sample of fewer than two elements are both defined to be zero.

-- $robust
--
-- These functions use the compensated summation algorithm of Chan et
-- al. for numerical robustness, but require two passes over the
-- sample data as a result.
--
-- Because of the need for two passes, these functions are /not/
-- subject to stream fusion.

-- | Maximum likelihood estimate of a sample's variance.  Also known
-- as the population variance, where the denominator is /n/.
variance :: (G.Vector v Double) => v Double -> Double
variance samp
    | n > 1     = robustSumVar (mean samp) samp / fromIntegral n
    | otherwise = 0
    where
      n = G.length samp
{-# SPECIALIZE variance :: U.Vector Double -> Double #-}
{-# SPECIALIZE variance :: V.Vector Double -> Double #-}

-- | Unbiased estimate of a sample's variance.  Also known as the
-- sample variance, where the denominator is /n/-1.
varianceUnbiased :: (G.Vector v Double) => v Double -> Double
varianceUnbiased samp
    | n > 1     = robustSumVar (mean samp) samp / fromIntegral (n-1)
    | otherwise = 0
    where
      n = G.length samp
{-# SPECIALIZE varianceUnbiased :: U.Vector Double -> Double #-}
{-# SPECIALIZE varianceUnbiased :: V.Vector Double -> Double #-}

-- | Standard deviation.  This is simply the square root of the
-- unbiased estimate of the variance.
stdDev :: (G.Vector v Double) => v Double -> Double
stdDev = sqrt . varianceUnbiased
{-# SPECIALIZE stdDev :: U.Vector Double -> Double #-}
{-# SPECIALIZE stdDev :: V.Vector Double -> Double #-}

-- $cancellation
--
-- The functions prefixed with the name @fast@ below perform a single
-- pass over the sample data using Knuth's algorithm. They usually
-- work well, but see below for caveats. These functions are subject
-- to array fusion.
--
-- /Note/: in cases where most sample data is close to the sample's
-- mean, Knuth's algorithm gives inaccurate results due to
-- catastrophic cancellation.

-- $references
--
-- * Chan, T. F.; Golub, G.H.; LeVeque, R.J. (1979) Updating formulae
--   and a pairwise algorithm for computing sample
--   variances. Technical Report STAN-CS-79-773, Department of
--   Computer Science, Stanford
--   University. <ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/79/773/CS-TR-79-773.pdf>
--
-- * Knuth, D.E. (1998) The art of computer programming, volume 2:
--   seminumerical algorithms, 3rd ed., p. 232.
--
-- * Welford, B.P. (1962) Note on a method for calculating corrected
--   sums of squares and products. /Technometrics/
--   4(3):419–420. <http://www.jstor.org/stable/1266577>
--
-- * West, D.H.D. (1979) Updating mean and variance estimates: an
--   improved method. /Communications of the ACM/
--   22(9):532–535. <http://doi.acm.org/10.1145/359146.359153>

max_stars_repo_name: runeksvendsen/hs-gauge

max_issues_repo_name: runeksvendsen/hs-gauge

max_forks_repo_name: runeksvendsen/hs-gauge


text:
{-# OPTIONS_GHC -Wall #-}

-- Dependant on cabal packages: sdl2, wave.

module Main where

import qualified Control.Concurrent as C
--import qualified Control.Monad as M
import qualified Data.Vector.Storable.Mutable as V
import qualified Data.Set as S
import Foreign.ForeignPtr as P

import qualified SDL
import qualified SDL.Audio as A

import qualified Codec.Audio.Wave as W

import qualified System.IO as IO

import qualified Statistics.Sample as St

micSpec :: IO.Handle -> A.OpenDeviceSpec
micSpec h = A.OpenDeviceSpec {A.openDeviceFreq = A.Mandate 48000
                             ,A.openDeviceFormat = A.Mandate A.Signed16BitNativeAudio
                             ,A.openDeviceChannels = A.Mandate A.Mono
                             ,A.openDeviceSamples = 4096
                             ,A.openDeviceCallback = \_ (V.MVector size ptr) -> P.withForeignPtr ptr (\p -> IO.hPutBuf h p size)
                             ,A.openDeviceUsage = A.ForCapture
                             ,A.openDeviceName = Nothing}

waveSpec :: W.Wave
waveSpec = W.Wave {W.waveFileFormat = W.WaveVanilla
                  , W.waveSampleRate = 48000
                  , W.waveSampleFormat = W.SampleFormatPcmInt 16
                  , W.waveChannelMask = S.singleton W.SpeakerFrontCenter
                  , W.waveDataOffset = 0
                  , W.waveDataSize = 0
                  , W.waveSamplesTotal = 0
                  , W.waveOtherChunks = []}

record :: IO.Handle -> IO ()
record h = do
  SDL.initialize [SDL.InitAudio]
  (dev, _) <- A.openAudioDevice $ micSpec h
  A.setAudioDevicePlaybackState dev A.Play
--  _ <- M.forever (C.threadDelay maxBound)
  _ <- C.threadDelay 10000000
  return ()

main :: IO ()
main =  W.writeWaveFile "mic.rec" waveSpec record

max_stars_repo_name: MiroslavVitkov/micid

max_issues_repo_name: MiroslavVitkov/micid

max_forks_repo_name: MiroslavVitkov/micid


text:
{-# LANGUAGE BangPatterns #-}
-- import System.Environment
import Data.Array.Repa as Repa
import Data.Array.Repa.Eval as Eval
import Data.Array.Repa.IO.BMP as BMP
import Data.Array.Repa.Algorithms.Pixel as Pixel
import Data.Vector.Unboxed.Base as Unboxed

import Control.Monad.Identity

import GHC.Word
import Data.Complex
import Data.List
import Julia
import Rendering
import RealFunctions
import ComplexFunctions

width :: (Num a) => a
width = 600

height :: (Num a) => a
height = 600

pixelRenderer :: (RealFloat a) => Int -> Int -> a
pixelRenderer x y = smooth samples -- avg' (Prelude.map (f) samples)
    where samples = sample 128 (magnitudeLimit 2) dragon start
          start   = afold (remapped x y)
          afold   = (1.5*(cis(pi/(4)))*).(foldPhase (pi/6))
          f z     = cos (atan2 (imagPart z) (realPart z)) --(cos . abs . phase) z
          f1 x z  = 1 + (z*x)
          f2 x z  = 1 - (z*x)

paramRenderer :: (RealFloat a) => a -> Int -> Int -> a
paramRenderer progress x y = smooth samples -- progress between 0 and 100%
    where samples = sample 128 (magnitudeLimit 2) (dragon) start
          start   = ((cis n *).(foldPhase (pi/6))) (remapped x y)
          n       = progress * 2 * pi

remapped :: (RealFloat a) => Int -> Int -> Complex a
remapped x y = (remap (fromIntegral x :+ fromIntegral y) width height)

avg :: (Num a, Fractional a) => [a] -> a
avg l = (sum l) / ((fromIntegral . length) l)

avg' :: (Num a, Fractional a) => [a] -> a
avg' = uncurry (/) . foldl' (\(s,n) x -> (s+x, n+1)) (0, 0)

-- 2D shaped array
shape :: DIM2
shape = Z :. (height :: Int) :. (width :: Int)

canvas :: Array U DIM2 Double
canvas = fromListUnboxed shape (take (fromIntegral (size shape)) (repeat 0))

render :: (Source r a, RealFloat a, Elt a) => Array r DIM2 a -> Array D DIM2 a
render !input = traverse input id (\_ (Z :. y :. x) -> pixelRenderer x y)

-- render a frame with progress between 0% and 100%
renderFrame :: (Source r a, RealFloat a, Elt a) => Array r DIM2 a -> a -> Array D DIM2 a
renderFrame !input progress = traverse input id (\_ (Z :. y :. x) -> paramRenderer progress x y)

convert :: (Shape sh) => Array U sh Double
     -> Array D sh (GHC.Word.Word8, GHC.Word.Word8, GHC.Word.Word8)
convert !array = Repa.map (rgb8OfGreyDouble) array

getMinimum :: (Ord b, Num b, Source r b, Elt b, Unbox b, Shape sh) =>
    Array r sh b -> IO (b)
getMinimum !array = foldAllP (min) 100000 array

getMaximum :: (Ord b, Num b, Source r b, Elt b, Unbox b, Shape sh) =>
    Array r sh b -> IO (b)
getMaximum !array = foldAllP (max) 0 array

correctExpo :: (Source r b, Unbox b, Elt b, Floating b, Shape sh, RealFrac b) =>
    Array r sh b -> b -> b -> IO (Array U sh b)
correctExpo !array !mini !maxi = computeP $ Repa.map ((map' mini maxi 0.0 1.0)) array
--
compute :: IO ()
compute = do
    img <- computeP (render canvas) :: IO (Array U DIM2 Double) -- super intensive
    mini <- getMinimum img
    maxi <- getMaximum img
    corrected <- correctExpo img mini maxi :: IO (Array U DIM2 Double)
    image <- computeP (convert corrected)
    writeImageToBMP "out.bmp" image

computeFrame :: Double -> Int -> IO ()
computeFrame progress frameNum = do
    img <- computeP (renderFrame canvas progress) :: IO (Array U DIM2 Double) -- super intensive
    mini <- getMinimum img
    maxi <- getMaximum img
    corrected <- correctExpo img mini maxi :: IO (Array U DIM2 Double)
    image <- computeP (convert corrected)
    let filename = "gifs/" Data.List.++ (show frameNum) Data.List.++ ".bmp"
    writeImageToBMP filename image

main :: IO ()
main = compute -- foldl (>>) (return ()) (Data.List.zipWith3 ($) (repeat computeFrame) progresses [1..])
    where progresses = takeWhile (<= 1.0) $ iterate (+0.005) 0.0

max_stars_repo_name: Vetii/Fractal

max_issues_repo_name: Vetii/Fractal

max_forks_repo_name: Vetii/Fractal


text:
-- |
-- Copyright: (c) 2020-2021 Tom Westerhout
-- SPDX-License-Identifier: BSD-3-Clause
-- Maintainer: Tom Westerhout <14264576+twesterhout@users.noreply.github.com>
module Numeric.PRIMME.Dense
  ( -- * Treat dense matrices as linear operators
    primmeFromDense,
  )
where

import Control.Monad.ST (RealWorld)
import Data.Complex
import Data.Proxy
import Data.Vector.Storable (MVector, Vector)
import qualified Data.Vector.Storable as V
import qualified Data.Vector.Storable.Mutable as MV
import Foreign.C.Types (CChar, CInt)
import Foreign.Marshal.Utils
import Foreign.Ptr (Ptr)
import Numeric.PRIMME.Types

type BlasInt = CInt

type BlasHemmType a =
  -- | SIDE
  Ptr CChar ->
  -- | UPLO
  Ptr CChar ->
  -- | M
  Ptr BlasInt ->
  -- | N
  Ptr BlasInt ->
  -- | ALPHA
  Ptr a ->
  -- | A
  Ptr a ->
  -- | LDA
  Ptr BlasInt ->
  -- | B
  Ptr a ->
  -- | LDB
  Ptr BlasInt ->
  -- | BETA
  Ptr a ->
  -- | C
  Ptr a ->
  -- | LDC
  Ptr BlasInt ->
  IO ()

foreign import ccall unsafe "ssymm_" ssymm_ :: BlasHemmType Float

foreign import ccall unsafe "dsymm_" dsymm_ :: BlasHemmType Double

foreign import ccall unsafe "chemm_" chemm_ :: BlasHemmType (Complex Float)

foreign import ccall unsafe "zhemm_" zhemm_ :: BlasHemmType (Complex Double)

hemm :: BlasDatatype a => Int -> Int -> a -> Vector a -> Int -> Vector a -> Int -> a -> MVector RealWorld a -> Int -> IO ()
hemm m n α a aStride b bStride β c cStride = do
  with (fromIntegral (fromEnum 'L') :: CChar) $ \side' ->
    with (fromIntegral (fromEnum 'U') :: CChar) $ \uplo' -> do
      with (fromIntegral m) $ \m' ->
        with (fromIntegral n) $ \n' ->
          with (fromIntegral aStride) $ \aStride' ->
            with (fromIntegral bStride) $ \bStride' ->
              with (fromIntegral cStride) $ \cStride' ->
                with α $ \α' ->
                  with β $ \β' ->
                    V.unsafeWith a $ \aPtr ->
                      V.unsafeWith b $ \bPtr ->
                        MV.unsafeWith c $ \cPtr ->
                          hemm' side' uplo' m' n' α' aPtr aStride' bPtr bStride' β' cPtr cStride'
  where
    hemm' :: forall b. BlasDatatype b => BlasHemmType b
    hemm' = case blasTag (Proxy :: Proxy b) of
      FloatTag -> ssymm_
      DoubleTag -> dsymm_
      ComplexFloatTag -> chemm_
      ComplexDoubleTag -> zhemm_

blockHemm :: BlasDatatype a => a -> Block a -> Block a -> a -> MBlock RealWorld a -> IO ()
blockHemm α (Block (m, n) aStride a) (Block (n', k) bStride b) β (MBlock (m', k') cStride c)
  | m /= n || m /= m' || n /= n' || k /= k' =
    error $
      "dimension mismatch: " <> show (m, n) <> " x " <> show (n', k) <> " = " <> show (m', k')
  | otherwise = hemm m' k' α a aStride b bStride β c cStride

-- | Treat a dense symmetric or Hermitian matrix as a operator. Internally we
-- call
-- [@?symm@](https://www.netlib.org/lapack/explore-html/db/dc9/group__single__blas__level3_ga8e8391a9873114d97e2b63e39fe83b2e.html#ga8e8391a9873114d97e2b63e39fe83b2e)
-- or
-- [@?hemm@](https://www.netlib.org/lapack/explore-html/db/def/group__complex__blas__level3_gad2d1853a142397404eae974b6574ece3.html#gad2d1853a142397404eae974b6574ece3)
-- BLAS functions so it is really important that the matrix is indeed Hermitian.
primmeFromDense ::
  BlasDatatype a =>
  -- | Matrix in column-major order
  Block a ->
  -- | Hermitian linear operator
  PrimmeOperator a
primmeFromDense a
  | isHermitian a = \b c -> blockHemm 1 a b 0 c
  | otherwise = error "expected a Hermitian matrix"

-- | Determine whether a matrix is Hermitian (or symmetric when @a@ is real).
--
-- /Note:/ this function returns 'False' for non-square matrices. No exceptions
-- are thrown.
isHermitian :: forall a. BlasDatatype a => Block a -> Bool
isHermitian (Block (n, n') stride v)
  | n == n' =
    -- Iterate in column-major order here and over the lower part of the matrix
    loop 0 (n - 1) $ \j ->
      loop (j + 1) n $ \i ->
        check i j
  | otherwise = False
  where
    access !i !j = v V.! (i + stride * j)
    check !i !j = access i j `unsafeEq` conj (access j i)
    loop !i !high f
      | i < high =
        if f i
          then loop (i + 1) high f
          else False
      | otherwise = True
    -- Yes, we really want exact binary comparison of floating point types
    unsafeEq :: a -> a -> Bool
    unsafeEq = case blasTag (Proxy :: Proxy a) of
      FloatTag -> (==)
      DoubleTag -> (==)
      ComplexFloatTag -> (==)
      ComplexDoubleTag -> (==)
    -- Return complex conjugate of a number.
    !conj = case blasTag (Proxy :: Proxy a) of
      FloatTag -> id
      DoubleTag -> id
      ComplexFloatTag -> Data.Complex.conjugate
      ComplexDoubleTag -> Data.Complex.conjugate

max_stars_repo_name: twesterhout/primme-hs

max_issues_repo_name: twesterhout/primme-hs

max_forks_repo_name: twesterhout/primme-hs


text:
{-# LANGUAGE NamedFieldPuns #-}

module School.App.MultiLayerPerceptron
( MLPOptions(..)
, multiLayerPerceptron
) where

import Numeric.LinearAlgebra (R)
import School.FileIO.FileHeader (FileHeader(..))
import School.FileIO.Source (source)
import School.Train.GradientDescent (gradientDescent)
import School.Train.IterationHandler (logCost)
import School.Train.SimpleDescentUpdate (simpleDescentUpdate)
import School.Train.StoppingCondition (maxIterations)
import School.Train.TrainState (TrainState(..), def)
import School.Types.FloatEq ((~=))
import School.Types.PingPong (PingPong, toPingPong)
import School.Types.DataType (DataType(..))
import School.Unit.Affine (affine)
import School.Unit.RecLin (recLin)
import School.Unit.MultiNoulli (multiNoulli)
import School.Unit.Unit (Unit)
import School.Unit.UnitParams (UnitParams(..))
import School.Unit.WeightDecay (weightDecay)
import School.Utils.Random (randAffineParams)
import System.FilePath (FilePath)
import System.Random (StdGen, mkStdGen)

data MLPOptions =
  MLPOptions { hiddenDims :: [Int]
             , initRate :: R
             , maxIter :: Int
             , nBatches :: Int
             , nClasses :: Int
             , nFeatures :: Int
             , randomSeed :: Int
             , weightDecayCoeff :: R
             , writeCost :: Maybe FilePath
             }

heLimits :: Int -> (R, R)
heLimits d = (0, sqrt (2 / fromIntegral d))

getUnits :: [Int]
         -> StdGen
         -> ([Unit R], [UnitParams R])
         -> ([Unit R], [UnitParams R])
getUnits (dim1:dim2:dim3:dims) gen (units, params) =
  getUnits (dim2:dim3:dims) gen' acc where
    addUnits = [affine, recLin]
    limits = heLimits dim1
    (affineParams, gen') =
      randAffineParams dim1 dim2 limits gen
    addParams = [affineParams, EmptyParams]
    acc = (units ++ addUnits, params ++ addParams)
getUnits (dim1:dim2:_) gen (units, params) =
  (units ++ [addUnit], params ++ [addParam]) where
    addUnit = affine
    limits = heLimits dim1
    (addParam, _) =
      randAffineParams dim1 dim2 limits gen
getUnits [_] _ acc = acc
getUnits [] _ acc = acc

setupUnits :: Int
           -> Int
           -> [Int]
           -> Int
           -> Either String
                     ([Unit R], PingPong (UnitParams R))
setupUnits nFeatures nClasses hiddenDims seed = do
  let allDims = nFeatures:hiddenDims ++ [nClasses]
  let gen = mkStdGen seed
  let (units, paramList) = getUnits allDims
                                    gen
                                    ([], [])
  params <- toPingPong paramList
  return (units, params)

multiLayerPerceptron :: MLPOptions
                     -> FilePath
                     -> IO (Either String
                           (Maybe [R], [UnitParams R]))
multiLayerPerceptron MLPOptions { nClasses = nC
                                , nBatches = nB
                                , nFeatures = nF
                                , hiddenDims
                                , initRate = learningRate
                                , maxIter
                                , randomSeed
                                , weightDecayCoeff
                                , writeCost
                                }
                      path = do
  let header = FileHeader { cols = nF
                            , rows = nB
                            , dataType = DBL64B
                            }
  let setup = setupUnits nF nC hiddenDims randomSeed
  case setup of
    Left e -> return . Left $ e
    Right (units, paramList) -> do
      let multi = multiNoulli Nothing Nothing
      let cost = if weightDecayCoeff ~= 0
                   then multi
                   else mappend multi (weightDecay weightDecayCoeff)
      let update = simpleDescentUpdate
      let condition = maxIterations maxIter
      let handler = case writeCost of
                      Nothing -> mempty
                      Just p -> logCost p
      let initState = def { paramList
                          , learningRate
                          }
      endState <- gradientDescent (source header path)
                                  units
                                  cost
                                  update
                                  condition
                                  handler
                                  initState
      return $ Left ""

max_stars_repo_name: jfulseca/School

max_issues_repo_name: jfulseca/School

max_forks_repo_name: jfulseca/School


text:
module Main (main) where

import Criterion.Types
import qualified Data.Aeson as Aeson
import qualified Data.Vector as V
import Properties
import Statistics.Types (estimateFromErr, mkCL)
import Test.Tasty (defaultMain, testGroup)
import Test.Tasty.HUnit (testCase)
import Test.HUnit

r1 :: Report
r1 = Report 0 "" [] v1 s1 (Outliers 0 0 0 0 0) []
 where
  m1 = Measured 4.613000783137977e-05 3.500000000000378e-05 31432 1 0 0 0 0.0 0.0 0.0 0.0
  v1 = V.fromList [m1]
  est1 = estimateFromErr 0.0 (0.0, 0.0) (mkCL 0.0)
  s1 = SampleAnalysis [] 0.0 est1 est1 (OutlierVariance Unaffected "" 0.0)

m2 :: Measured
m2 = Measured {measTime = 1.1438998626545072e-5
              , measCpuTime = 1.2000000001677336e-5
              , measCycles = 6208
              , measIters = 1

              , measAllocated = minBound
              , measNumGcs = minBound
              , measBytesCopied = minBound

              , measMutatorWallSeconds = -1/0
              , measMutatorCpuSeconds = -1/0
              , measGcWallSeconds = -1/0
              , measGcCpuSeconds = -1/0}

main :: IO ()
main = defaultMain $ testGroup "Tests"
       [ Properties.tests
       , testCase "json-roundtrip1"
           (assertEqual "round trip simple Measured"
              (Right m2) (Aeson.eitherDecode (Aeson.encode m2)))
       , testCase "json-roundtrip2"
           (assertEqual "round trip simple Report"
              (Right r1) (Aeson.eitherDecode (Aeson.encode r1)))
       ]

max_stars_repo_name: bgamari/criterion

max_issues_repo_name: bgamari/criterion

max_forks_repo_name: bgamari/criterion