text:
{-# LANGUAGE CPP
           , GADTs
           , KindSignatures
           , TypeOperators
           , TypeFamilies
           , EmptyCase
           , DataKinds
           , PolyKinds
           , ExistentialQuantification
           , FlexibleContexts
           , OverloadedStrings
           #-}

{-# OPTIONS_GHC -Wall -fwarn-tabs #-}

module Language.Hakaru.Sample where

import           Numeric.SpecFunctions            (logFactorial)
import qualified Data.Number.LogFloat             as LF
import qualified Math.Combinatorics.Exact.Binomial as EB
-- import qualified Numeric.Integration.TanhSinh     as TS
import qualified System.Random.MWC                as MWC
import qualified System.Random.MWC.CondensedTable as MWC
import qualified System.Random.MWC.Distributions  as MWCD

import qualified Data.Vector                      as V
import           Data.STRef
import           Data.Sequence (Seq)
import qualified Data.Foldable                    as F
import qualified Data.List.NonEmpty               as L
import           Data.List.NonEmpty               (NonEmpty(..))
import           Data.Maybe                       (fromMaybe)

#if __GLASGOW_HASKELL__ < 710
import           Control.Applicative   (Applicative(..), (<$>))
#endif
import           Control.Monad
import           Control.Monad.ST
import           Control.Monad.Identity
import           Control.Monad.Trans.Maybe
import           Control.Monad.State.Strict
import qualified Data.IntMap                      as IM

import Data.Number.Nat     (fromNat)
import Data.Number.Natural (fromNatural, fromNonNegativeRational, Natural, unsafeNatural)
import Language.Hakaru.Types.DataKind
import Language.Hakaru.Types.Coercion
import Language.Hakaru.Types.Sing
import Language.Hakaru.Types.HClasses
import Language.Hakaru.Syntax.IClasses
import Language.Hakaru.Syntax.TypeOf
import Language.Hakaru.Syntax.Value
import Language.Hakaru.Syntax.Reducer
import Language.Hakaru.Syntax.Datum
import Language.Hakaru.Syntax.DatumCase
import Language.Hakaru.Syntax.AST
import Language.Hakaru.Syntax.ABT

data EAssoc =
    forall a. EAssoc {-# UNPACK #-} !(Variable a) !(Value a)

newtype Env = Env (IM.IntMap EAssoc)

emptyEnv :: Env
emptyEnv = Env IM.empty

updateEnv :: EAssoc -> Env -> Env
updateEnv v@(EAssoc x _) (Env xs) =
    Env $ IM.insert (fromNat $ varID x) v xs

updateEnvs
    :: List1 Variable xs
    -> List1 Value xs
    -> Env
    -> Env
updateEnvs Nil1         Nil1         env = env
updateEnvs (Cons1 x xs) (Cons1 y ys) env =
    updateEnvs xs ys (updateEnv (EAssoc x y) env)

lookupVar :: Variable a -> Env -> Maybe (Value a)
lookupVar x (Env env) = do
    EAssoc x' e' <- IM.lookup (fromNat $ varID x) env
    Refl         <- varEq x x'
    return e'

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

-- Makes use of Atkinson's algorithm as described in:
-- Monte Carlo Statistical Methods pg. 55
--
-- Further discussion at:
-- http://www.johndcook.com/blog/2010/06/14/generating-poisson-random-values/
poisson_rng :: Double -> MWC.GenIO -> IO Int
poisson_rng lambda g' = make_poisson g'
    where
    smu   = sqrt lambda
    b     = 0.931 + 2.53*smu
    a     = -0.059 + 0.02483*b
    vr    = 0.9277 - 3.6224/(b - 2)
    arep  = 1.1239 + 1.1368/(b - 3.4)
    lnlam = log lambda

    make_poisson :: MWC.GenIO -> IO Int
    make_poisson g = do
        u <- MWC.uniformR (-0.5,0.5) g
        v <- MWC.uniformR (0,1) g
        let us = 0.5 - abs u
            k = floor $ (2*a / us + b)*u + lambda + 0.43
        case () of
            () | us >= 0.07 && v <= vr -> return k
            () | k < 0                 -> make_poisson g
            () | us <= 0.013 && v > us -> make_poisson g
            () | accept_region us v k  -> return k
            _                          -> make_poisson g

    accept_region :: Double -> Double -> Int -> Bool
    accept_region us v k =
        log (v * arep / (a/(us*us)+b))
        <=
        -lambda + fromIntegral k * lnlam - logFactorial k

normalize :: [Value 'HProb] -> (LF.LogFloat, Double, [Double])
normalize []          = (0, 0, [])
normalize [(VProb x)] = (x, 1, [1])
normalize xs          = (m, y, ys)
    where
    xs' = map (\(VProb x) -> x) xs
    m   = maximum xs'
    ys  = [ LF.fromLogFloat (x/m) | x <- xs' ]
    y   = sum ys

normalizeVector
    :: Value ('HArray 'HProb) -> (LF.LogFloat, Double, V.Vector Double)
normalizeVector (VArray xs) =
    let xs' = V.map (\(VProb x) -> x) xs in
    case V.length xs of
    0 -> (0, 0, V.empty)
    1 -> (V.unsafeHead xs', 1, V.singleton 1)
    _ ->
        let m   = V.maximum xs'
            ys  = V.map (\x -> LF.fromLogFloat (x/m)) xs'
            y   = V.sum ys
        in (m, y, ys)

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

runEvaluate
    :: (ABT Term abt)
    => abt '[] a
    -> Value a
runEvaluate prog = evaluate prog emptyEnv

evaluate
    :: (ABT Term abt)
    => abt '[] a
    -> Env
    -> Value a
evaluate e env = caseVarSyn e (evaluateVar env) (flip evaluateTerm env)

evaluateVar :: Env -> Variable a -> Value a
evaluateVar env v =
    case lookupVar v env of
    Nothing -> error "variable not found!"
    Just a  -> a

evaluateTerm
    :: (ABT Term abt)
    => Term abt a
    -> Env
    -> Value a
evaluateTerm t env =
    case t of
    o :$          es -> evaluateSCon    o es    env
    NaryOp_  o    es -> evaluateNaryOp  o es    env
    Literal_ v       -> evaluateLiteral v
    Empty_   _       -> evaluateEmpty
    Array_   n    es -> evaluateArray   n es    env
    ArrayLiteral_ es -> VArray . V.fromList $ map (flip evaluate env) es
    Bucket b e    rs -> evaluateBucket  b e  rs env
    Datum_   d       -> evaluateDatum   d       env
    Case_    o    es -> evaluateCase    o es    env
    Superpose_    es -> evaluateSuperpose es    env
    Reject_  _       -> VMeasure $ \_ _ -> return Nothing

evaluateSCon
    :: (ABT Term abt)
    => SCon args a
    -> SArgs abt args
    -> Env
    -> Value a
evaluateSCon Lam_ (e1 :* End) env =
    caseBind e1 $ \x e1' ->
        VLam $ \v -> evaluate e1' (updateEnv (EAssoc x v) env)
evaluateSCon App_ (e1 :* e2 :* End) env =
    case evaluate e1 env of
    VLam f -> f (evaluate e2 env)
evaluateSCon Let_ (e1 :* e2 :* End) env =
    let v = evaluate e1 env
    in caseBind e2 $ \x e2' ->
        evaluate e2' (updateEnv (EAssoc x v) env)
evaluateSCon (CoerceTo_   c) (e1 :* End) env =
    coerceTo c $ evaluate e1 env
evaluateSCon (UnsafeFrom_ c) (e1 :* End) env =
    coerceFrom c $ evaluate e1 env
evaluateSCon (PrimOp_ o)     es env = evaluatePrimOp    o es env
evaluateSCon (ArrayOp_ o)    es env = evaluateArrayOp   o es env
evaluateSCon (MeasureOp_  m) es env = evaluateMeasureOp m es env
evaluateSCon Dirac           (e1 :* End) env =
    VMeasure $ \p _ -> return $ Just (evaluate e1 env, p)
evaluateSCon MBind (e1 :* e2 :* End) env =
    case evaluate e1 env of
    VMeasure m1 -> VMeasure $ \ p g -> do
        x <- m1 p g
        case x of
            Nothing -> return Nothing
            Just (a, p') ->
                caseBind e2 $ \x' e2' ->
                    case evaluate e2' (updateEnv (EAssoc x' a) env) of
                    VMeasure y -> y p' g

evaluateSCon Plate (n :* e2 :* End) env =
    case evaluate n env of
    VNat n' -> caseBind e2 $ \x e' ->
        VMeasure $ \(VProb p) g -> runMaybeT $ do
            (v', ps) <- fmap V.unzip . V.mapM (performMaybe g) $
                V.generate (fromInteger $ fromNatural n') $ \v ->
                    evaluate e' $
                    updateEnv (EAssoc x . VNat $ intToNatural v) env
            return
                ( VArray v'
                , VProb $ p * V.product (V.map (\(VProb y) -> y) ps)
                )
    where
    performMaybe
        :: MWC.GenIO
        -> Value ('HMeasure a)
        -> MaybeT IO (Value a, Value 'HProb)
    performMaybe g (VMeasure m) = MaybeT $ m (VProb 1) g

evaluateSCon Chain (n :* s :* e :* End) env =
    case (evaluate n env, evaluate s env) of
    (VNat n', start) ->
        caseBind e $ \x e' ->
            let s' = VLam $ \v -> evaluate e' (updateEnv (EAssoc x v) env) in
            VMeasure (\(VProb p) g -> runMaybeT $ do
                (evaluates, sout) <- runStateT (replicateM (unsafeInt n') $ convert g s') start
                let (v', ps) = unzip evaluates
                    bodyType :: Sing ('HMeasure (HPair a b)) -> Sing ('HArray a)
                    bodyType = SArray . fst . sUnPair . sUnMeasure
                return
                    ( VDatum $ dPair_ (bodyType $ caseBind e (const typeOf)) (typeOf s)
                        (VArray . V.fromList $ v') sout
                    , VProb $ p * product (map (\(VProb y) -> y) ps)
                    ))
    where
    convert
        :: MWC.GenIO
        -> Value (s ':-> 'HMeasure (HPair a s))
        -> StateT (Value s) (MaybeT IO) (Value a, Value 'HProb)
    convert g (VLam f) = StateT $ \s' ->
        case f s' of
        VMeasure f' -> do
            (as'', p') <- MaybeT (f' (VProb 1) g)
            let (a, s'') = unPair as''
            return ((a, p'), s'')

    unPair :: Value (HPair a b) -> (Value a, Value b)
    unPair (VDatum (Datum "pair" _typ
        (Inl (Et (Konst a)
            (Et (Konst b) Done))))) = (a, b)
    unPair x = case x of {}

evaluateSCon (Summate hd hs) (e1 :* e2 :* e3 :* End) env =
    case (evaluate e1 env, evaluate e2 env) of
    (lo, hi) ->
        caseBind e3 $ \x e3' ->
            foldl (\t i ->
                   evalOp (Sum  hs) t $
                     evaluate e3' (updateEnv (EAssoc x i) env))
                  (identityElement $ Sum hs)
                  (enumFromUntilValue hd lo hi)

evaluateSCon (Product hd hs) (e1 :* e2 :* e3 :* End) env =
    case (evaluate e1 env, evaluate e2 env) of
    (lo, hi) ->
        caseBind e3 $ \x e3' ->
            foldl (\t i ->
                   evalOp (Prod hs) t $
                     evaluate e3' (updateEnv (EAssoc x i) env))
                  (identityElement $ Prod hs)
                  (enumFromUntilValue hd lo hi)

evaluateSCon s _ _ = error $ "TODO: evaluateSCon{" ++ show s ++ "}"

evaluatePrimOp
    ::  ( ABT Term abt, typs ~ UnLCs args, args ~ LCs typs)
    => PrimOp typs a
    -> SArgs abt args
    -> Env
    -> Value a
evaluatePrimOp Not (e1 :* End) env =
    case evaluate e1 env of
      VDatum a -> if a == dTrue
                  then VDatum dFalse
                  else VDatum dTrue

evaluatePrimOp Pi  End         _   = VProb . LF.logFloat $ pi
evaluatePrimOp Cos (e1 :* End) env =
    case evaluate e1 env of
      VReal v1 -> VReal . cos $ v1

evaluatePrimOp Sin (e1 :* End) env =
    case evaluate e1 env of
      VReal v1 -> VReal . sin $ v1

evaluatePrimOp Tan (e1 :* End) env =
    case evaluate e1 env of
      VReal v1 -> VReal . tan $ v1

evaluatePrimOp RealPow (e1 :* e2 :* End) env =
    case (evaluate e1 env, evaluate e2 env) of
      (VProb v1, VReal v2) -> VProb $ LF.pow v1 v2

evaluatePrimOp Choose (e1 :* e2 :* End) env =
    case (evaluate e1 env, evaluate e2 env) of
      (VNat v1, VNat v2) -> VNat $ EB.choose v1 v2

evaluatePrimOp Exp (e1 :* End) env =
    case evaluate e1 env of
      VReal v1 -> VProb . LF.logToLogFloat $ v1

evaluatePrimOp Log (e1 :* End) env =
    case evaluate e1 env of
      VProb v1 -> VReal . LF.logFromLogFloat $ v1

evaluatePrimOp (Infinity h) End _ =
    case h of
      HIntegrable_Nat  -> error "Can not evaluate infinity for natural numbers"
      HIntegrable_Prob -> VProb $ LF.logFloat LF.infinity

evaluatePrimOp (Equal _) (e1 :* e2 :* End) env = (VDatum . dBool) $ evaluate e1 env == evaluate e2 env

evaluatePrimOp (Less _) (e1 :* e2 :* End) env =
    case (evaluate e1 env, evaluate e2 env) of
    (VNat  v1, VNat  v2) -> VDatum $ if v1 < v2 then dTrue else dFalse
    (VInt  v1, VInt  v2) -> VDatum $ if v1 < v2 then dTrue else dFalse
    (VProb v1, VProb v2) -> VDatum $ if v1 < v2 then dTrue else dFalse
    (VReal v1, VReal v2) -> VDatum $ if v1 < v2 then dTrue else dFalse
    _                    -> error "TODO: evaluatePrimOp{Less}"
evaluatePrimOp (NatPow _) (e1 :* e2 :* End) env =
    case evaluate e2 env of
    VNat  v2 ->
        let v2' = fromNatural v2 in
        case evaluate e1 env of
          VNat  v1 -> VNat  (v1 ^ v2')
          VInt  v1 -> VInt  (v1 ^ v2')
          VProb v1 -> VProb (v1 ^ v2')
          VReal v1 -> VReal (v1 ^ v2')
          _        -> error "NatPow should always return some kind of number"
evaluatePrimOp (Negate _) (e1 :* End) env =
    case evaluate e1 env of
    VInt  v -> VInt  (negate v)
    VReal v -> VReal (negate v)
    v       -> case v of {}
evaluatePrimOp (Abs   _) (e1 :* End) env =
    case evaluate e1 env of
    VInt  v -> VNat  . unsafeNatural   $ abs v
    VReal v -> VProb . LF.logFloat $ abs v
    v       -> case v of {}
evaluatePrimOp (Recip _) (e1 :* End) env =
    case evaluate e1 env of
    VProb v -> VProb (recip v)
    VReal v -> VReal (recip v)
    v       -> case v of {}
evaluatePrimOp (NatRoot _) (e1 :* e2 :* End) env =
    case (evaluate e1 env, evaluate e2 env) of
    (VProb v1, VNat v2) -> VProb $ LF.pow v1 (recip . fromIntegral $ v2)
    v                   -> case v of {}

evaluatePrimOp (Floor) (e1 :* End) env =
    case (evaluate e1 env) of
    VProb v1 -> VNat (floor (LF.fromLogFloat v1))

evaluatePrimOp prim _ _ =
    error ("TODO: evaluatePrimOp{" ++ show prim ++ "}")

evaluateArrayOp
    :: ( ABT Term abt
       , typs ~ UnLCs args
       , args ~ LCs typs)
    => ArrayOp typs a
    -> SArgs abt args
    -> Env
    -> Value a
evaluateArrayOp (Index _) = \(e1 :* e2 :* End) env ->
    case (evaluate e1 env, evaluate e2 env) of
    (VArray v, VNat n) -> v V.! unsafeInt n

evaluateArrayOp (Size _) = \(e1 :* End) env ->
    case evaluate e1 env of
    VArray v -> VNat . intToNatural $ V.length v

evaluateArrayOp (Reduce _) = \(e1 :* e2 :* e3 :* End) env ->
    case ( evaluate e1 env
         , evaluate e2 env
         , evaluate e3 env) of
    (f, a, VArray v) -> V.foldl' (lam2 f) a v

evaluateMeasureOp
    :: ( ABT Term abt
       , typs ~ UnLCs args
       , args ~ LCs typs)
    => MeasureOp typs a
    -> SArgs abt args
    -> Env
    -> Value ('HMeasure a)

evaluateMeasureOp Lebesgue = \(e1 :* e2 :* End) env ->
  case (evaluate e1 env, evaluate e2 env) of
    (VReal v1, VReal v2) | v1 < v2 ->
      VMeasure $ \(VProb p) g ->
        case (isInfinite v1, isInfinite v2) of
          (False, False) -> do
            x <- MWC.uniformR (v1, v2) g
            return $ Just (VReal $ x,
                           VProb $ p * LF.logFloat (v2 - v1))
          (False, True) -> do
            u <- MWC.uniform g
            let l = log u
            let n = -l
            return $ Just (VReal $ v1 + n,
                           VProb $ p * LF.logToLogFloat n)
          (True, False) -> do
            u <- MWC.uniform g
            let l = log u
            let n = -l
            return $ Just (VReal $ v2 - n,
                           VProb $ p * LF.logToLogFloat n)
          (True, True) -> do
            (u,b) <- MWC.uniform g
            let l = log u
            let n = -l
            return $ Just (VReal $ if b then n else l,
                           VProb $ p * 2 * LF.logToLogFloat n)
    (VReal _, VReal _) -> error "Lebesgue with length 0 or flipped endpoints"

evaluateMeasureOp Counting = \End _ ->
    VMeasure $ \(VProb p) g -> do
        let success = LF.logToLogFloat (-3 :: Double)
        let pow x y = LF.logToLogFloat (LF.logFromLogFloat x *
                                       (fromIntegral y :: Double))
        u' <- MWCD.geometric0 (LF.fromLogFloat success) g
        let u = toInteger u'
        b <- MWC.uniform g
        return $ Just
            ( VInt  $ if b then -1-u else u
            , VProb $ p * 2 / pow (1-success) u / success)

evaluateMeasureOp Categorical = \(e1 :* End) env ->
    VMeasure $ \p g -> do
        let (_,y,ys) = normalizeVector (evaluate e1 env)
        if not (y > (0::Double)) -- TODO: why not use @y <= 0@ ??
        then error "Categorical needs positive weights"
        else do
            u <- MWC.uniformR (0, y) g
            return $ Just
                ( VNat
                . intToNatural
                . fromMaybe 0
                . V.findIndex (u <=)
                . V.scanl1' (+)
                $ ys
                , p)

evaluateMeasureOp Uniform = \(e1 :* e2 :* End) env ->
    case (evaluate e1 env, evaluate e2 env) of
    (VReal v1, VReal v2) -> VMeasure $ \p g -> do
        x <- MWC.uniformR (v1, v2) g
        return $ Just (VReal x, p)

evaluateMeasureOp Normal = \(e1 :* e2 :* End) env ->
    case (evaluate e1 env, evaluate e2 env) of
    (VReal v1, VProb v2) -> VMeasure $ \ p g -> do
        x <- MWCD.normal v1 (LF.fromLogFloat v2) g
        return $ Just (VReal x, p)

evaluateMeasureOp Poisson = \(e1 :* End) env ->
    case evaluate e1 env of
    VProb v1 -> VMeasure $ \ p g -> do
        x <- MWC.genFromTable (MWC.tablePoisson (LF.fromLogFloat v1)) g
        return $ Just (VNat $ intToNatural x, p)

evaluateMeasureOp Gamma = \(e1 :* e2 :* End) env ->
    case (evaluate e1 env, evaluate e2 env) of
    (VProb v1, VProb v2) -> VMeasure $ \ p g -> do
        x <- MWCD.gamma (LF.fromLogFloat v1) (LF.fromLogFloat v2) g
        return $ Just (VProb $ LF.logFloat x, p)

evaluateMeasureOp Beta = \(e1 :* e2 :* End) env ->
    case (evaluate e1 env, evaluate e2 env) of
    (VProb v1, VProb v2) -> VMeasure $ \ p g -> do
        x <- MWCD.beta (LF.fromLogFloat v1) (LF.fromLogFloat v2) g
        return $ Just (VProb $ LF.logFloat x, p)

evaluateNaryOp
    :: (ABT Term abt)
    => NaryOp a -> Seq (abt '[] a) -> Env -> Value a
evaluateNaryOp s es =
    F.foldr (evalOp s) (identityElement s) . mapEvaluate es

identityElement :: NaryOp a -> Value a
identityElement And                   = VDatum dTrue
identityElement (Sum HSemiring_Nat)   = VNat  0
identityElement (Sum HSemiring_Int)   = VInt  0
identityElement (Sum HSemiring_Prob)  = VProb 0
identityElement (Sum HSemiring_Real)  = VReal 0
identityElement (Prod HSemiring_Nat)  = VNat  1
identityElement (Prod HSemiring_Int)  = VInt  1
identityElement (Prod HSemiring_Prob) = VProb 1
identityElement (Prod HSemiring_Real) = VReal 1
identityElement (Max  HOrd_Prob)      = VProb 0
identityElement (Max  HOrd_Real)      = VReal LF.negativeInfinity
identityElement (Min  HOrd_Prob)      = VProb (LF.logFloat LF.infinity)
identityElement (Min  HOrd_Real)      = VReal LF.infinity
identityElement _                     = error "Missing identity elements?"

evalOp
    :: NaryOp a -> Value a -> Value a -> Value a
evalOp And (VDatum a) (VDatum b)
    | a == dTrue && b == dTrue = VDatum dTrue
    | otherwise = VDatum dFalse
evalOp (Sum  HSemiring_Nat)  (VNat  a) (VNat  b) = VNat  (a + b)
evalOp (Sum  HSemiring_Int)  (VInt  a) (VInt  b) = VInt  (a + b)
evalOp (Sum  HSemiring_Prob) (VProb a) (VProb b) = VProb (a + b)
evalOp (Sum  HSemiring_Real) (VReal a) (VReal b) = VReal (a + b)
evalOp (Prod HSemiring_Nat)  (VNat  a) (VNat  b) = VNat  (a * b)
evalOp (Prod HSemiring_Int)  (VInt  a) (VInt  b) = VInt  (a * b)
evalOp (Prod HSemiring_Prob) (VProb a) (VProb b) = VProb (a * b)
evalOp (Prod HSemiring_Real) (VReal a) (VReal b) = VReal (a * b)
evalOp (Max  HOrd_Prob)      (VProb a) (VProb b) = VProb (max a b)
evalOp (Max  HOrd_Real)      (VReal a) (VReal b) = VReal (max a b)
evalOp (Min  HOrd_Prob)      (VProb a) (VProb b) = VProb (min a b)
evalOp (Min  HOrd_Real)      (VReal a) (VReal b) = VReal (min a b)

evalOp op                    _          _        =
    error ("TODO: evalOp{" ++ show op ++ "}")

mapEvaluate
    :: (ABT Term abt)
    => Seq (abt '[] a) -> Env -> Seq (Value a)
mapEvaluate es env = fmap (flip evaluate env) es

evaluateLiteral :: Literal a -> Value a
evaluateLiteral (LNat  n) = VNat  . fromInteger $ fromNatural n -- TODO: catch overflow errors
evaluateLiteral (LInt  n) = VInt  $ fromInteger n -- TODO: catch overflow errors
evaluateLiteral (LProb n) = VProb . fromRational $ fromNonNegativeRational n
evaluateLiteral (LReal n) = VReal $ fromRational n

evaluateEmpty :: Value ('HArray a)
evaluateEmpty = VArray V.empty

evaluateArray
    :: (ABT Term abt)
    => (abt '[] 'HNat)
    -> (abt '[ 'HNat ] a)
    -> Env
    -> Value ('HArray a)
evaluateArray n e env =
    case evaluate n env of
    VNat n' -> caseBind e $ \x e' ->
        VArray $ V.generate (unsafeInt n') $ \v ->
            let v' = VNat $ intToNatural v in
            evaluate e' (updateEnv (EAssoc x v') env)

evaluateBucket
    :: (ABT Term abt)
    => abt '[] 'HNat
    -> abt '[] 'HNat
    -> Reducer abt '[] a
    -> Env
    -> Value a
evaluateBucket b e rs env =
    case (evaluate b env, evaluate e env) of
      (VNat b', VNat e') -> runST $ do
          s' <- init Nil1 rs env
          mapM_ (\i -> accum (VNat i) Nil1 rs s' env) [b' .. e' - 1]
          done s'
    where init :: (ABT Term abt)
               => List1 Value xs
               -> Reducer abt xs a
               -> Env
               -> ST s (VReducer s a)
          init ix (Red_Fanout r1 r2)    env  =
              VRed_Pair (type_ r1) (type_ r2) <$> init ix r1 env <*> init ix r2 env
          init ix (Red_Index  n  _  mr) env' =
              let (vars, n') = caseBinds n in
              case evaluate n' (updateEnvs vars ix env') of
                VNat n'' -> VRed_Array <$> V.generateM (fromIntegral n'')
                            (\bb -> init (Cons1 (vnat bb) ix) mr env')
          init ix (Red_Split _ r1 r2)   env' =
              VRed_Pair (type_ r1) (type_ r2) <$> init ix r1 env <*> init ix r2 env'
          init _  Red_Nop               _    = return VRed_Unit
          init _  (Red_Add h _) _ = VRed_Num <$> newSTRef (identityElement (Sum h))

          type_ = typeOfReducer

          vnat :: Int -> Value 'HNat
          vnat  = VNat . fromIntegral

          accum :: (ABT Term abt)
                => Value 'HNat
                -> List1 Value xs
                -> Reducer abt xs a
                -> VReducer s a
                -> Env
                -> ST s ()
          accum n ix (Red_Fanout r1 r2)   (VRed_Pair _ _ v1 v2) env' =
              accum n ix r1 v1 env >> accum n ix r2 v2 env'
          accum n ix (Red_Index n' a1 r2) (VRed_Array v)          env' =
              caseBind a1 $ \i a1' ->
              let (vars, a1'') = caseBinds a1'
                  VNat ov = evaluate a1''
                            (updateEnv (EAssoc i n) (updateEnvs vars ix env'))
                  ov' = fromIntegral ov in
              accum n (Cons1 (VNat ov) ix) r2 (v V.! ov') env
          accum n ix (Red_Split bb r1 r2) (VRed_Pair _ _ v1 v2) env' =
              caseBind bb $ \i b' ->
                  let (vars, b'') = caseBinds b' in
                  case evaluate b''
                       (updateEnv (EAssoc i n) (updateEnvs vars ix env')) of
                  VDatum bb -> if bb == dTrue then
                                   accum n ix r1 v1 env'
                               else
                                   accum n ix r2 v2 env'
          accum n ix (Red_Add h ee) (VRed_Num s) env' =
              caseBind ee $ \i e' ->
                  let (vars, e'') = caseBinds e'
                      v = evaluate e''
                          (updateEnv (EAssoc i n) (updateEnvs vars ix env')) in
                  modifySTRef' s (evalOp (Sum h) v)
          accum _ _ Red_Nop _ _ = return ()
          accum _ _ _ _ _ = error "Some impossible combinations happened?"

          done :: VReducer s a -> ST s (Value a)
          done (VRed_Num s)            = readSTRef s
          done VRed_Unit               = return (VDatum dUnit)
          done (VRed_Pair s1 s2 v1 v2) = do
            v1' <- done v1
            v2' <- done v2
            return (VDatum $ dPair_ s1 s2 v1' v2')
          done (VRed_Array v)          = VArray <$> V.sequence (V.map done v)

evaluateDatum
    :: (ABT Term abt)
    => Datum (abt '[]) (HData' a)
    -> Env
    -> Value (HData' a)
evaluateDatum d env = VDatum (fmap11 (flip evaluate env) d)

evaluateCase
    :: forall abt a b
    .  (ABT Term abt)
    => abt '[] a
    -> [Branch a abt b]
    -> Env
    -> Value b
evaluateCase o es env =
    case runIdentity $ matchBranches evaluateDatum' (evaluate o env) es of
    Just (Matched rho b) ->
        evaluate b (extendFromMatch (fromAssocs rho) env)
    _ -> error "Missing cases in match expression"
    where
    extendFromMatch :: [Assoc Value] -> Env -> Env
    extendFromMatch []                env' = env'
    extendFromMatch (Assoc x v : xvs) env' =
        extendFromMatch xvs (updateEnv (EAssoc x v) env')

    evaluateDatum' :: DatumEvaluator Value Identity
    evaluateDatum' = return . Just . getVDatum

    getVDatum :: Value (HData' a) -> Datum Value (HData' a)
    getVDatum (VDatum a) = a

evaluateSuperpose
    :: (ABT Term abt)
    => NonEmpty (abt '[] 'HProb, abt '[] ('HMeasure a))
    -> Env
    -> Value ('HMeasure a)
evaluateSuperpose ((q, m) :| []) env =
    case evaluate m env of
    VMeasure m' ->
        let VProb q' = evaluate q env
        in  VMeasure (\(VProb p) g -> m' (VProb $ p * q') g)

evaluateSuperpose pms@((_, m) :| _) env =
    case evaluate m env of
    VMeasure m' ->
        let pms'     = L.toList pms
            weights  = map ((flip evaluate env) . fst) pms'
            (x,y,ys) = normalize weights
        in VMeasure $ \(VProb p) g ->
            if not (y > (0::Double)) then return Nothing else do
            u <- MWC.uniformR (0, y) g
            case [ m1 | (v,(_,m1)) <- zip (scanl1 (+) ys) pms', u <= v ] of
                m2 : _ ->
                    case evaluate m2 env of
                    VMeasure m2' -> m2' (VProb $ p * x * LF.logFloat y) g
                []     -> m' (VProb $ p * x * LF.logFloat y) g

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

-- Useful 'short-hand'
intToNatural :: Int -> Natural
intToNatural = unsafeNatural . toInteger

unsafeInt :: Natural -> Int
unsafeInt = fromInteger . fromNatural
----------------------------------------------------------- fin.

max_stars_repo_name: vmchale/hakaru

max_issues_repo_name: vmchale/hakaru

max_forks_repo_name: vmchale/hakaru


text:
{-# LANGUAGE Haskell2010
    , TypeFamilies
    , FlexibleContexts
    , Trustworthy
    , StandaloneDeriving
    , DeriveDataTypeable
    , CPP
 #-}
{-# OPTIONS -Wall -fno-warn-name-shadowing #-}
{-# OPTIONS -cpp  -pgmPcpphs  -optP--cpp #-}

module Numeric.Vector (

    Vector

) where

import Data.Ratio
import Data.Complex
import Data.Int

import Data.Array.IArray
import Data.Array.Unboxed

import Data.Typeable

import qualified Data.Array.Unsafe as U

data family Vector e

#if defined(__GLASGOW_HASKELL__) && (__GLASGOW_HASKELL__ >= 707)
deriving instance Typeable Vector
#else
deriving instance Typeable1 Vector
#endif

data instance Vector Int
    = IntVector !Int (UArray Int Int)

data instance Vector Float
    = FloatVector !Int (UArray Int Float)

data instance Vector Double
    = DoubleVector !Int (UArray Int Double)

data instance Vector Integer
    = IntegerVector !Int (Array Int Integer)

data instance Vector (Ratio a)
    = RatioVector !Int (Array Int (Ratio a))

data instance Vector (Complex a)
    = ComplexVector !Int (Array Int (Complex a))

-- <.>
dotProd = undefined

-- ><
vectorProd = undefined

ext: hs

max_stars_repo_name: phadej/bed-and-breakfast

max_issues_repo_name: phadej/bed-and-breakfast

max_forks_repo_name: phadej/bed-and-breakfast


text:
{-# LANGUAGE FlexibleContexts #-}
module Statistics.Distribution.Random.Gamma (
   gamma
) where

{-
 -  Shape parameter a >= 1.  Algorithm GD in:
 -
 -  Ahrens, J.H. and Dieter, U. (1982).
 -  Generating gamma variates by a modified
 -  rejection technique.
 -  Comm. ACM, 25, 47-54.
 -
 -
 -  Shape parameter 0 < a < 1. Algorithm GS in:
 -
 -  Ahrens, J.H. and Dieter, U. (1974).
 -  Computer methods for sampling from gamma, beta,
 -  poisson and binomial distributions.
 -  Computing, 12, 223-246.
 -}

import qualified Statistics.Distribution.Random.Exponential as D
import qualified Statistics.Distribution.Random.Normal as D
import qualified Statistics.Distribution.Random.Uniform as D
import           Data.List (foldl')
import           Data.Number.LogFloat (expm1)
import           Data.Word
import           Random.CRI

sqrt32, exp_m1 :: Double
sqrt32 = sqrt 32
{-# NOINLINE sqrt32 #-}
exp_m1 = exp (-1)
{-# NOINLINE exp_m1 #-}

q1, q2, q3, q4, q5, q6, q7 :: Double
q1 =  0.04166669
q2 =  0.02083148
q3 =  0.00801191
q4 =  0.00144121
q5 = -7.388e-5
q6 =  2.4511e-4
q7 =  2.424e-4
qs :: [Double]
qs = [q7, q6, q5, q4, q3, q2, q1]
{-# INLINE qs #-}

a1, a2, a3, a4, a5, a6, a7 :: Double
a1 =  0.3333333
a2 = -0.250003
a3 =  0.2000062
a4 = -0.1662921
a5 =  0.1423657
a6 = -0.1367177
a7 =  0.1233795
as :: [Double]
as = [a7, a6, a5, a4, a3, a2, a1]
{-# INLINE as #-}

horner :: [Double] -> Double -> Double
horner q r = foldl' (\ a b -> (a + b) * r) 0 q
{-# INLINE horner #-}

gamma :: (Source m g Double, Source m g Word32) => Double -> Double -> g m -> m Double
gamma shape scale rng
    | shape == 0    =  return 0
    | shape <  1    =  gammaGS shape scale rng
    | otherwise     =  gammaGD shape scale rng

gammaGS :: (Source m g Double, Source m g Word32) => Double -> Double -> g m -> m Double
gammaGS shape scale rng =
    let e  = 1 + exp_m1 * shape
        go = do
            ru <- D.uniform     rng
            re <- D.exponential rng
            let p                  = e * ru
                x      | p >= 1    = - log ((e - p) / shape)
                       | otherwise = exp (log p / shape)
                accept | p >= 1    = re >= (1 - shape) * log x
                       | otherwise = re >= x
            if accept then return (scale * x) else go
    in go

gammaGD :: (Source m g Double, Source m g Word32) => Double -> Double -> g m -> m Double
gammaGD shape scale rng =

    let s2 = shape - 0.5
        s  = sqrt s2
        d  = sqrt32 - s * 12
        r  = 1 / shape
        q0 = horner qs r

        (b, si, c)
            | shape <=  3.686 = (0.463 + s + 0.178 * s2, 1.235,            0.195 / s - 0.079 + 0.16 * s)
            | shape <= 13.022 = (1.654 + 0.0076 * s2,    1.68 / s + 0.275, 0.062 / s + 0.024)
            | otherwise       = (1.77,                   0.75,             0.1515 / s)

        calc_q t
            | abs v <= 0.25   = q0 + 0.5 * t * t * horner as v
            | otherwise       = q0 - s * t + 0.25 * t * t + (s2 + s2) * log (1 + v)
            where v = t / (s + s)

        choose_t = do
            e  <- D.exponential rng
            u' <- D.uniform rng
            let uu = u' + u' - 1
            let tt = if uu < 0 then b - si * e else b + si * e
            if tt >= -0.71874483771719
                then do
                    let qq = calc_q tt
                    if qq > 0 && c * abs uu <= expm1 qq * exp (e - 0.5 * tt * tt)
                        then return tt
                        else choose_t
                else choose_t

    in do
        t <- D.normal rng
        let x = s + 0.5 * t
            ret_val = scale * x * x
        if t >= 0
            then return ret_val
            else do

                u <- D.uniform rng
                if d * u <= t * t * t
                    then return ret_val
                    else do

                        let q = calc_q t
                        if x > 0 && log (1 - u) <= q
                            then return ret_val
                            else do
                                ttt <- choose_t
                                let xx = s + 0.5 * ttt
                                return $ scale * xx * xx

max_stars_repo_name: finlay/random-dist

max_issues_repo_name: finlay/random-dist

max_forks_repo_name: finlay/random-dist


text:
{-|
  Copyright   :  (C) 2019, Google Inc.
  License     :  BSD2 (see the file LICENSE)
  Maintainer  :  Christiaan Baaij <christiaan.baaij@gmail.com>
-}
{-# LANGUAGE TemplateHaskell #-}

{-# OPTIONS_GHC -Wno-orphans #-}
-- {-# OPTIONS_GHC -ddump-splices #-}

module Clash.Class.AutoReg.Instances where

import           Clash.Class.AutoReg.Internal
import           Clash.CPP                           (maxTupleSize)

import           Data.Complex (Complex)
import           Data.Ord (Down)
import           Data.Ratio (Ratio)

deriveAutoReg ''Complex
deriveAutoReg ''Down
deriveAutoReg ''Ratio

-- | __N.B.__: The documentation only shows instances up to /3/-tuples. By
-- default, instances up to and including /12/-tuples will exist. If the flag
-- @large-tuples@ is set instances up to the GHC imposed limit will exist. The
-- GHC imposed limit is either 62 or 64 depending on the GHC version.
deriveAutoRegTuples [2..maxTupleSize]

max_stars_repo_name: oliverbunting/clash-compiler

max_issues_repo_name: oliverbunting/clash-compiler

max_forks_repo_name: oliverbunting/clash-compiler


text:
{-# LANGUAGE ConstraintKinds       #-}
{-# LANGUAGE DerivingVia           #-}
{-# Language FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# Language MultiParamTypeClasses #-}
{-# LANGUAGE StandaloneDeriving    #-}
{-# LANGUAGE UndecidableInstances  #-}
module Geometry.Intersections
  ( Intersections (..)
  , intersections, isIntersecting
  ) where

import Data.Complex
import Data.List (nub)

import Geometry.Base
import Geometry.Line
import Geometry.Circle
import Geometry.Polygon
import Geometry.Plot
import Geometry.Decorations

--------------------------------------------------------------------------------
-- | Class provides `intersections'` function returning a list (possible empty)
-- of all potential intersection points (co-dimension 1).
class (Manifold a, Manifold b) => Intersections a b where
  intersections' :: a -> b -> [Cmp]

-- | Filters intersection points, keeping ones belonging to given manifolds.
intersections :: (Figure a, Figure b, Intersections a b)
  => a -> b -> [Cmp]
intersections a b
  | isTrivial a = filter (isContaining b) [cmp $ refPoint a]
  | isTrivial b = filter (isContaining a) [cmp $ refPoint b]
  | otherwise =
    nub $
    filter (isContaining a) $
    filter (isContaining b) $
    intersections' a b

-- | Returns `True` if tho curves have intersection points.
isIntersecting :: (Figure a, Figure b, Intersections a b) => a -> b -> Bool
isIntersecting a b = not . null $ intersections a b

instance (Figure a, Figure b, Intersections a b) =>
  Intersections a (Maybe b) where
  intersections' a = maybe [] (intersections' a)

instance (Figure a, Figure b, Intersections b a) =>
  Intersections (Maybe a) b where
  intersections' = flip intersections'

instance {-# OVERLAPPING #-}
  (Figure a, Figure b, Intersections b a) =>
  Intersections (Maybe a) (Maybe b) where
  intersections' a = maybe [] (intersections' a)

instance (Figure a, Figure b, Intersections a b) =>
         Intersections a (Decorated b) where
  intersections' a = intersections a . fromDecorated

instance (Figure a, Figure b, Intersections b a) =>
         Intersections (Decorated a) b where
  intersections' = flip intersections'

instance {-# OVERLAPPING #-}
  (Figure a, Figure b, Intersections a b) =>
  Intersections (Decorated a) (Decorated b) where
  intersections' = intersections' . fromDecorated

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

intersectionCC c1 c2 | d == 0 || b < 0 = []
                     | b == 0 = [x]
                     | b > 0 = [x, conjugate x]
  where
    r1 = radius c1
    r2 = radius c2
    d = magnitude (center c1 - center c2)
    a = r1**2-r2**2+d**2
    b = -(r2-r1-d)*(r2-r1+d)*(r2+r1-d)*(r2+r1+d)
    x = scale (1/(2*d)) $ a :+ sqrt b

intersectionCL c l | b > r = []
                   | b == r = [0:+b]
                   | b < r = [a:+b, (-a):+b]
  where
    b = getY (refPoint l)
    r = radius c
    a = sqrt (r**2 - b**2)

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

instance Intersections Line Line where
  intersections' l1 l2
    | l2 `isContaining` refPoint l1 = [refPoint l1]
    | l1 `isContaining` refPoint l2 = [refPoint l2]
    | otherwise = lineIntersection l1 l2

instance Intersections Line Circle  where
  intersections' l cir =
    invt <$> intersectionCL (t cir) (t l)
    where
      t :: Trans x => x -> x
      t = rotate (-a) . translate' (negate (center cir))
      invt = translate' (center cir) . rotate a
      a = angle l

instance Intersections Line Polyline where
  intersections' l = foldMap (intersections l) . segments . closePolyline

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

instance Intersections Circle Line where
  intersections' = flip intersections'

instance Intersections Circle Circle where
  intersections' cir1 cir2 =
    invt <$> intersectionCC (t cir1) (t cir2)
    where
      t = rotate (-a) . translate' (negate c1)
      invt = translate' c1 . rotate a
      c1 = center cir1
      c2 = center cir2
      a = azimuth c1 c2

instance Intersections Circle Polyline where
  intersections' c = foldMap (intersections c) . segments . closePolyline

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

instance Intersections Polyline Line where
  intersections' = flip intersections'

instance Intersections Polyline Circle where
  intersections' = flip intersections'

instance Intersections Polyline Polyline where
  intersections' p = foldMap (intersections p) . segments . closePolyline

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

deriving via Line instance
  (Manifold a, Intersections a Line) => Intersections a Ray

instance Intersections Ray Line where
  intersections' = intersections' . asLine

instance Intersections Ray Polyline where
  intersections' = intersections' . asLine

instance Intersections Ray Circle where
  intersections' = intersections' . asLine

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

deriving via Line instance
  (Manifold a, Intersections a Line) => Intersections a Segment

instance Intersections Segment Line where
  intersections' = intersections' . asLine

instance Intersections Segment Polyline where
  intersections' = intersections' . asLine

instance Intersections Segment Circle where
  intersections' = intersections' . asLine

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

deriving via Polyline instance
  (Manifold a, Intersections a Polyline) => Intersections a Polygon

instance Intersections Polygon Line  where
  intersections' = intersections' . asPolyline

instance Intersections Polygon Polyline  where
  intersections' = intersections' . asPolyline

instance Intersections Polygon Circle  where
  intersections' = intersections' . asPolyline

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

deriving via Polygon instance
  (Manifold a, Intersections a Polyline) => Intersections a Triangle

instance Intersections Triangle Line  where
  intersections' = intersections' . asPolyline

instance Intersections Triangle Polyline  where
  intersections' = intersections' . asPolyline

instance Intersections Triangle Circle  where
  intersections' = intersections' . asPolyline

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

deriving via Polygon instance
  (Manifold a, Intersections a Polyline) => Intersections a Rectangle

instance Intersections Rectangle Line  where
  intersections' = intersections' . asPolyline

instance Intersections Rectangle Polyline  where
  intersections' = intersections' . asPolyline

instance Intersections Rectangle Circle  where
  intersections' = intersections' . asPolyline

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

instance Pnt a => Intersections (Plot a) Line where
  intersections' p l = intersectionsP p l $ intersections (asPolyline p) l

instance Pnt a => Intersections (Plot a) Circle where
  intersections' p c = intersectionsP p c $ intersections (asPolyline p) c

instance Pnt a => Intersections (Plot a) Polyline where
  intersections' p pl = intersectionsP p pl $ intersections (asPolyline p) pl

instance (Pnt a, Manifold b, Intersections (Plot a) b) => Intersections b (Plot a) where
  intersections' = flip intersections'

instance {-# OVERLAPPING #-}
  (Pnt a, Pnt b) => Intersections (Plot a) (Plot b) where
  intersections' p1 p2 = intersectionsP p1 p2 $
                         intersections (asPolyline p1) (asPolyline p2)

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

deriving via (Plot a) instance
  (Manifold b, Affine a, Trans a, Intersections b (Plot a)) => Intersections b (ClosedPlot a)

instance Pnt a => Intersections (ClosedPlot a) Line  where
  intersections' p l = intersectionsP p l $ intersections (asPolyline p) l

instance Pnt a => Intersections (ClosedPlot a) Polyline  where
  intersections' p pl = intersectionsP p pl $ intersections (asPolyline p) pl

instance Pnt a => Intersections (ClosedPlot a) Circle  where
  intersections' p c = intersectionsP p c $ intersections (asPolyline p) c

max_stars_repo_name: MethaHardworker/geometry

max_issues_repo_name: MethaHardworker/geometry

max_forks_repo_name: MethaHardworker/geometry


text:
{-# LANGUAGE FlexibleContexts, TypeFamilies #-}
module Main where

import Data.Complex (Complex(..), magnitude)
import qualified Diagrams.Prelude as D

import qualified Diagrams.Backend.Rasterific.CmdLine as CmdLine

main :: IO ()
main = CmdLine.mainWith (picture 32 2.0)

type C = Complex Double

-- calc
mandelbrotSet :: Traversable t =>
                 Int -> Double -> t C -> t Bool
mandelbrotSet n inf lst = fmap (testMandelbrot n inf) lst

testMandelbrot :: Int -> Double -> C -> Bool
testMandelbrot n inf c = isDiverged (mandelbrot c) n inf

isDiverged :: (C -> C) -> Int -> Double -> Bool
isDiverged fn n inf = filter (\z -> magnitude z > inf) seq == []
  where seq = take n $ iterate fn 0

mandelbrot :: C -> C -> C
mandelbrot c z = z*z + c

-- draw
toPixel :: Bool -> Int
toPixel True = 0
toPixel False = 1

toSquare n =
  D.opacity (fromIntegral n)
  $ D.fc D.black
  $ D.lw D.medium
  $ D.square 1

side n v0 v1 =
   let sv = (v1 - v0) / fromIntegral n
   in  [v0, (v0 + sv) .. v1]

grid :: Int -> Double -> Double -> [[C]]
grid n v0 v1 = map (\y -> map (:+ y) s) s
  where s = side n v0 v1

image n inf = map (map (toSquare . toPixel . (testMandelbrot n inf))) $ grid 64 (-2) 2

picture :: Int -> Double -> D.Diagram CmdLine.B
picture n inf = D.bgFrame 3 D.pink $ (D.vcat . map D.hcat $ image n inf)

max_stars_repo_name: cormojs/simple-mandelbrot

max_issues_repo_name: cormojs/simple-mandelbrot

max_forks_repo_name: cormojs/simple-mandelbrot


text:
{-|
Module: MachineLearning.PCA
Description: Principal Component Analysis.
Copyright: (c) Alexander Ignatyev, 2017-2018.
License: BSD-3
Stability: experimental
Portability: POSIX

Principal Component Analysis (PCA) - dimensionality reduction algorithm.
It is mostly used to speed up supervising learning (Regression, Classification, etc) and visualization of data.
-}

module MachineLearning.PCA
(
  getDimReducer
  , getDimReducer_rv
)

where

import Prelude hiding ((<>))
import Data.Maybe (fromMaybe)
import qualified Data.Vector.Storable as V
import qualified Numeric.LinearAlgebra as LA
import Numeric.LinearAlgebra ((<>), (??))
import qualified MachineLearning as ML

import MachineLearning.Types (R, Vector, Matrix)

-- | Computes "covariance matrix".
covarianceMatrix :: Matrix -> Matrix
covarianceMatrix x = ((LA.tr x) <> x) / (fromIntegral $ LA.rows x)

-- | Compute eigenvectors (matrix U) and singular values (matrix S) of the given covariance matrix.
pca :: Matrix -> (Matrix, Vector)
pca x = (u, s)
  where sigma = covarianceMatrix x
        (u, s, _) = LA.svd sigma

-- | Gets dimensionality reduction function, retained variance (0..1) and reduced X
-- for given matrix X and number of dimensions to retain.
getDimReducer :: Matrix -> Int -> (Matrix -> Matrix, R, Matrix)
getDimReducer x k = (reducer, retainedVariance, reducer xNorm)
  where muSigma = ML.meanStddev x
        xNorm = ML.featureNormalization muSigma x
        (u, s) = pca xNorm
        u' = u ?? (LA.All, LA.Take k)
        reducer xx = (ML.featureNormalization muSigma xx) <> u'
        retainedVariance = (V.sum $ V.slice 0 k s) / (V.sum s)

-- | Gets dimensionality reduction function, retained number of dimensions and reduced X
-- for given matrix X and variance to retain (0..1].
getDimReducer_rv :: Matrix -> R -> (Matrix -> Matrix, Int, Matrix)
getDimReducer_rv x rv = (reducer, k, reducer xNorm)
  where muSigma = ML.meanStddev x
        xNorm = ML.featureNormalization muSigma x
        (u, s) = pca xNorm
        sum_s = V.sum s
        variances = (V.scanl' (+) 0 s) / (LA.scalar sum_s)
        k = fromMaybe ((V.length s) - 1) $ V.findIndex (\v -> v >= rv) variances
        u' = u ?? (LA.All, LA.Take k)
        reducer xx = (ML.featureNormalization muSigma xx) <> u'

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE TupleSections #-}
{-# LANGUAGE OverloadedLabels #-}
{-# OPTIONS_GHC -Wall #-}

-- | Quantile chart patterns.
module Chart.Quantile
  (
    -- * quantile chart patterns
    quantileChart,
    histChart,

    -- * testing
    qua,
    qLines,
  )
where

import Chart
import Optics.Core
import Data.Foldable
import Data.Maybe
import Data.Text (Text)
import Prelude hiding (abs)
import Statistics.Distribution.Normal
import Statistics.Distribution
import Data.Bool

-- * charts

-- | Chart template for quantile data.
quantileChart ::
  [Text] ->
  [LineStyle] ->
  [AxisOptions] ->
  [[Double]] ->
  ChartSvg
quantileChart names ls as xs =
  mempty & #charts .~ chart0 & #hudOptions .~ hudOptions'
  where
    hudOptions' :: HudOptions
    hudOptions' =
      defaultHudOptions
        & ( #legends
              .~
                [(12, defaultLegendOptions
                    & #textStyle % #size .~ 0.1
                    & #vgap .~ 0.05
                    & #innerPad .~ 0.2
                    & #place .~ PlaceRight
                    & #content .~ zip names (fmap (\l -> LineChart l [[Point 0 0, Point 1 1]]) ls)
                )]
          )
        & set #axes ((5,) <$> as)

    chart0 = unnamed $
      zipWith (\s d -> LineChart s [d])
        ls
        (zipWith Point [0 ..] <$> xs)

-- | histogram chart
histChart ::
  Range Double ->
  Int ->
  [Double] ->
  ChartSvg
histChart r g xs =
  barChart defaultBarOptions barData'
  where
    barData' = BarData [freqs] xs'' []
    hcuts = grid OuterPos r g
    h = fill hcuts xs
    counts =
      (\(Rect _ _ _ w) -> w)
        <$> makeRects (IncludeOvers (width r / fromIntegral g)) h
    freqs = (/sum counts) <$> counts
    xs' =
      (\(Rect x x' _ _) -> (x + x') / 2)
        <$> makeRects (IncludeOvers (width r / fromIntegral g)) h
    xs'' = ["unders"] <> take (length xs' - 2) (drop 1 (comma (Just 2) <$> xs')) <> ["overs"]

qua :: Double -> Double -> Double -> Double -> Double
qua u s t p
  | t <= 0 = u*t
  | otherwise = Statistics.Distribution.quantile (normalDistr (u*t) (s*sqrt t)) p

qLines :: Double -> Int -> Colour -> [LineStyle]
qLines s n c = (\x -> defaultLineStyle & #color .~ x & #size .~ s) <$> cqs n c

cqs :: Int -> Colour -> [Colour]
cqs n c = fmap (\x -> mix x c (greyed c)) xs
  where
    xs = (\t -> fmap (\x -> fromIntegral (abs (x-(t `div` 2)-bool 0 1 (x>(t `div`2) && 1==t `mod` 2))) / fromIntegral (t `div` 2)) [0..t]) (n - 1)


text:
{-|
Module: MachineLearning
Description: Machine Learning
Copyright: (c) Alexander Ignatyev, 2016
License: BSD-3
Stability: experimental
Portability: POSIX

-}

module MachineLearning
(
  addBiasDimension
  , removeBiasDimension
  , meanStddev
  , featureNormalization
  , mapFeatures
  , splitToXY
)

where

import MachineLearning.Types (Vector, Matrix)
import qualified Numeric.LinearAlgebra as LA
import Numeric.LinearAlgebra((|||), (??))
import qualified Numeric.Morpheus.Statistics as Stat

import Control.Monad (replicateM, mfilter, MonadPlus)
import Data.List (sort, group, foldl')
import qualified Data.Vector as V

-- | Add biad dimension to the future matrix
addBiasDimension :: Matrix -> Matrix
addBiasDimension x = 1 ||| x

-- | Remove biad dimension
removeBiasDimension :: Matrix -> Matrix
removeBiasDimension x = x ?? (LA.All, LA.Drop 1)

-- | Caclulates mean and stddev values of every feature.
-- Takes feature matrix X, returns pair of vectors of means and stddevs.
meanStddev x =
  let means = Stat.columnMean x
      stddevs = Stat.columnStddev_m means x
  in (LA.asRow means, LA.asRow stddevs)

featureNormalization (means, stddevs) x = (x - means) / stddevs

-- | Maps the features into all polynomial terms of X up to the degree-th power
mapFeatures :: Int -> Matrix -> Matrix
mapFeatures 1 x = x
mapFeatures degree x = LA.fromColumns $ cols ++ (foldl' (\l d -> (terms d) ++ l) [] [degree, degree-1 .. 2])
  where cols = LA.toColumns x
        vv = V.fromList cols
        ncols = V.length vv
        makeTerm :: [(Int, Int)] -> Vector
        makeTerm = foldl' (\c (index, power) -> c * (vv V.! index) ^ power) 1
        terms :: Int -> [Vector]
        terms d = foldl' (\l x -> (makeTerm x) : l) [] $ polynomialTerms d [ncols-1, ncols-2 .. 0]

polynomialTerms :: Ord a => Int -> [a] -> [[(a, Int)]]
polynomialTerms degree terms =
  map (\x -> map (\y -> (head y, length y)) $ group x)
  $ combinationsWithReplacement degree terms

combinationsWithReplacement :: (MonadPlus m, Ord a) => Int -> m a -> m [a]
combinationsWithReplacement sample objects = mfilter (\a -> sort a == a) $ replicateM sample objects

-- | Splits data matrix to features matrix X and vector of outputs y
splitToXY m =
  let x = m ?? (LA.All, LA.DropLast 1)
      y = LA.flatten $ m ?? (LA.All, LA.TakeLast 1)
  in (x, y)


text:
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE DeriveAnyClass        #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE Strict                #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}
{-|
Module      : Grenade.Layers.Relu
Description : Rectifying linear unit layer
Copyright   : (c) Huw Campbell, 2016-2017
License     : BSD2
Stability   : experimental
-}
module Grenade.Layers.Relu (
    Relu (..)
  , SpecRelu (..)
  , specRelu1D
  , specRelu2D
  , specRelu3D
  , relu
  ) where

import           Control.Parallel.Strategies
import           Data.Constraint                (Dict (..))
import           Data.Reflection                (reifyNat)
import           Data.Serialize
import           Data.Singletons
import           GHC.Generics                   (Generic)
import           GHC.TypeLits
import qualified Numeric.LinearAlgebra.Static   as LAS
import           Unsafe.Coerce                  (unsafeCoerce)

import           Grenade.Core
import           Grenade.Dynamic
import           Grenade.Dynamic.Internal.Build
import           Grenade.Utils.Conversion
import           Grenade.Utils.Vector

-- | A rectifying linear unit.
--   A layer which can act between any shape of the same dimension, acting as a
--   diode on every neuron individually.
data Relu = Relu
  deriving (Generic, NFData, Show)

instance UpdateLayer Relu where
  type Gradient Relu = ()
  runUpdate _ _ _ = Relu

instance RandomLayer Relu where
  createRandomWith _ _ = return Relu

instance Serialize Relu where
  put _ = return ()
  get = return Relu

-- Run forward 1D
runForward1D :: (KnownNat i) => Relu -> S ('D1 i) -> (Tape Relu ('D1 i) ('D1 i), S ('D1 i))
runForward1D _ (S1D y) = (S1D y, S1D (relu y))
    where
      relu = LAS.dvmap (\a -> if a <= 0 then 0 else a)
runForward1D _ (S1DV y) =
  (S1DV y, S1DV (mapVectorInPlace relu y)) -- y will be replaced, which however still leads to a correct implementation!
  -- force (mapVectorInPlace relu y) `seq` (S1DV y, S1DV y) -- y will be replaced, which however still leads to a correct implementation!
  -- (S1DV y, S1DV $ mapVector relu y) -- y will be replaced, which however still leads to a correct implementation!
    where
      relu a = if a <= 0 then 0 else a

reluDifZip :: (Ord a, Floating a) => a -> a -> a
reluDifZip a g = if a <= 0 then 0 else g

-- Run backward 1D
runBackward1D :: (KnownNat i) => Relu -> S ('D1 i) -> S ('D1 i) -> (Gradient Relu, S ('D1 i))
runBackward1D _ (S1D y) (S1D dEdy) = ((), S1D (relu' y * dEdy))
    where
      relu' = LAS.dvmap (\a -> if a <= 0 then 0 else 1)
runBackward1D _ (S1DV y) (S1DV dEdy) = ((), zipWithVectorInPlaceSnd reluDifZip y dEdy `seq` S1DV dEdy)
runBackward1D l y dEdy = runBackward1D l y (toLayerShape y dEdy)

instance (KnownNat i) => Layer Relu ('D1 i) ('D1 i) where
  type Tape Relu ('D1 i) ('D1 i) = S ('D1 i)
  runForwards = runForward1D
  runBackwards = runBackward1D

runBackward2D :: (KnownNat i, KnownNat j) => Relu -> S ('D2 i j) -> S ('D2 i j) -> (Gradient Relu, S ('D2 i j))
runBackward2D _ (S2D y) (S2D dEdy) = ((), S2D (relu' y * dEdy))
    where
      relu' = LAS.dmmap (\a -> if a <= 0 then 0 else 1)
runBackward2D _ (S2DV y) (S2DV dEdy) = ((), zipWithVectorInPlaceSnd reluDifZip y dEdy `seq` S2DV dEdy)
runBackward2D l y dEdy = runBackward2D l y (toLayerShape y dEdy)

instance (KnownNat i, KnownNat j) => Layer Relu ('D2 i j) ('D2 i j) where
  type Tape Relu ('D2 i j) ('D2 i j) = S ('D2 i j)

  runForwards _ (S2D y) = (S2D y, S2D (relu y))
    where
      relu = LAS.dmmap (\a -> if a <= 0 then 0 else a)
  runForwards _ (S2DV y) = (S2DV y, S2DV $ mapVectorInPlace relu y ) -- This replaces the y vector, but that doesn't harm the update
    where
      relu a = if a <= 0 then 0 else a
  runBackwards = runBackward2D

instance (KnownNat i, KnownNat j, KnownNat k) => Layer Relu ('D3 i j k) ('D3 i j k) where

  type Tape Relu ('D3 i j k) ('D3 i j k) = S ('D3 i j k)

  runForwards _ (S3D y) = (S3D y, S3D (relu y))
    where
      relu = LAS.dmmap (\a -> if a <= 0 then 0 else a)
  runBackwards _ (S3D y) (S3D dEdy) = ((), S3D (relu' y * dEdy))
    where
      relu' = LAS.dmmap (\a -> if a <= 0 then 0 else 1)

-------------------- DynamicNetwork instance --------------------

instance FromDynamicLayer Relu where
  fromDynamicLayer inp _ _ = SpecNetLayer $ SpecRelu (tripleFromSomeShape inp)

instance ToDynamicLayer SpecRelu where
  toDynamicLayer _ _ (SpecRelu (rows, cols, depth)) =
     reifyNat rows $ \(_ :: (KnownNat rows) => Proxy rows) ->
     reifyNat cols $ \(_ :: (KnownNat cols) => Proxy cols) ->
     reifyNat depth $ \(_ :: (KnownNat depth) => Proxy depth) ->
     case (rows, cols, depth) of
         (_, 1, 1)    -> return $ SpecLayer Relu (sing :: Sing ('D1 rows)) (sing :: Sing ('D1 rows))
         (_, _, 1) -> return $ SpecLayer Relu (sing :: Sing ('D2 rows cols)) (sing :: Sing ('D2 rows cols))
         _    -> case (unsafeCoerce (Dict :: Dict()) :: Dict (KnownNat (rows GHC.TypeLits.* depth))) of
           Dict -> return $ SpecLayer Relu (sing :: Sing ('D3 rows cols depth)) (sing :: Sing ('D3 rows cols depth))

-- | Create a specification for a elu layer.
specRelu1D :: Integer -> SpecNet
specRelu1D i = specRelu3D (i, 1, 1)

-- | Create a specification for a elu layer.
specRelu2D :: (Integer, Integer) -> SpecNet
specRelu2D (i, j) = specRelu3D (i, j, 1)

-- | Create a specification for a elu layer.
specRelu3D :: (Integer, Integer, Integer) -> SpecNet
specRelu3D = SpecNetLayer . SpecRelu

-- | Add a Relu layer to your build.
relu :: BuildM ()
relu = buildGetLastLayerOut >>= buildAddSpec . SpecNetLayer . SpecRelu

-------------------- GNum instances --------------------

instance GNum Relu where
  _ |* Relu = Relu
  _ |+ Relu = Relu
  sumG _ = Relu

max_stars_repo_name: schnecki/grenade

max_issues_repo_name: schnecki/grenade

max_forks_repo_name: schnecki/grenade


text:
module Numeric.Layer.Sigmoid where
import Numeric.Layer
import Dependent.Size
import Numeric.Vector.Sized
import qualified Numeric.Vector.Sized as Sized
import GHC.TypeLits(KnownNat,Nat)

newtype Sigmoid (n :: Nat) = Sigmoid ()

instance Show (Sigmoid n) where
  showsPrec d (Sigmoid _) = showString "Sigmoid"

sigmoid :: Floating x => x -> x
sigmoid x = 1 / (1 + exp (-x))

sigmoid' :: Num a => a -> a
sigmoid' y = (1 - y) * y

instance KnownNat n => Layer (Sigmoid n) where
  type Inputs (Sigmoid n) = ZZ ::. n
  type Outputs (Sigmoid n) = ZZ ::. n
  type Tape (Sigmoid n) = SizedArray (ZZ ::. n)
  type Gradient (Sigmoid n) = ()
  forward _ x = (y, tape)
    where
      y = Sized.map sigmoid x
      tape = y
  backward _ tape dy = (gradient, dx)
    where
      gradient = ()
      dx = Sized.zipWith (*) dy (Sized.map sigmoid' tape)
  applyGradient layer _ = layer

max_stars_repo_name: mixed-signals/mixed-signals

max_issues_repo_name: mixed-signals/mixed-signals

max_forks_repo_name: mixed-signals/mixed-signals


text:
module Main where

import Linear (
  V2(V2), V3(V3), (^+^), (^-^), (*^))
import Data.Complex (
  Complex((:+)))
import Codec.Picture (
  Image, Pixel8, generateImage, writePng)

data Edge2 a = Edge2 (V2 a) (V2 a)
  deriving (Show, Eq, Ord)

type ScalarField2 a = V2 a -> a

renderE0 :: Floating a => Ord a => ScalarField2 a
renderE0 (V2 x1 x2)
  | x1 < 0 = 0
  | x1 > 1 = 0
  | otherwise = x2 * exp (negate (sqr x2))

sqr :: Num a => a -> a
sqr x = x * x

cotransformScalarField2 :: (V2 a -> V2 a) -> ScalarField2 a -> ScalarField2 a
cotransformScalarField2 g f = \x -> f (g x)

translate2 :: Num a => V2 a -> V2 a -> V2 a
translate2 = (^+^)

rotate2 :: RealFloat a => Complex a -> V2 a -> V2 a
rotate2 r (V2 x1 x2) = do
  let y1 :+ y2 = r * (x1 :+ x2)
  V2 y1 y2

scale22 :: Num a => a -> V2 a -> V2 a
scale22 s (V2 x1 x2) = V2 x1 (s * x2)

renderEdge2 :: RealFloat a => Ord a => Edge2 a -> ScalarField2 a
renderEdge2 (Edge2 a b) = do
  let w = 0.1
  let V2 c1 c2 = b ^-^ a
  let r = c1 :+ c2
  cotransformScalarField2 (
    scale22 (recip w) . rotate2 (recip r) . translate2 (negate a)) renderE0

-- Write image

resolution :: Int
resolution = 600

imageToCameraSpace :: V2 Int -> V2 Double
imageToCameraSpace (V2 i1 i2) = do
  let x1 = fromIntegral i1 / fromIntegral resolution * 2 - 1
  let x2 = 1 - fromIntegral i2 / fromIntegral resolution * 2
  V2 x1 x2

discreteGrey :: Double -> Pixel8
discreteGrey a
  | a < -1 = 0
  | a > 1 = 255
  | otherwise = round ((0.5 + 0.5 * a) * 255)

manifestImage :: ScalarField2 Double -> Image Pixel8
manifestImage rendering = do
  let imageFunction i1 i2 = discreteGrey (rendering (imageToCameraSpace (V2 i1 i2)))
  generateImage imageFunction resolution resolution

main :: IO ()
main = do

  let testImage = renderEdge2 (Edge2 (V2 (-0.8) 0.3) (V2 0.4 0.8))

  writePng "out/testImage.png" (manifestImage testImage)

max_stars_repo_name: phischu/stereo-vision

max_issues_repo_name: phischu/stereo-vision

max_forks_repo_name: phischu/stereo-vision


text:
module Parser (readExpr, readExprList) where

import Types

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

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

spaces :: Parser ()
spaces = skipMany1 space

escapedChar :: Parser String
escapedChar = do backslash <- char '\\'
                 escaped <- oneOf "nrt\"\\"
                 return [backslash, escaped]

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

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

parseNumber :: Parser LispVal
parseNumber = fmap (Number . read) (many1 digit)

parseHashPrefix :: Parser LispVal
parseHashPrefix = char '#' >> (parseOctal
                               <|> parseDec
                               <|> parseHex
                               <|> parseBin
                               <|> parseChar
                               <|> parseBool)

-- Simple types parsing
parseNumberHelper prefix filter reader = do char prefix
                                            number <- many1 filter
                                            return $ (Number . reader) number

parseOctal :: Parser LispVal
parseOctal = parseNumberHelper 'o' octDigit (fst . head . readOct)

parseHex :: Parser LispVal
parseHex = parseNumberHelper 'h' hexDigit (fst . head . readHex)

parseDec :: Parser LispVal
parseDec = parseNumberHelper 'd' digit read

parseBin :: Parser LispVal
parseBin = parseNumberHelper 'b' (oneOf "01") binToDec
  where binToDec = foldl (\acc x -> acc * 2 + (if x == '0' then 0 else 1)) 0

parseChar :: Parser LispVal
parseChar = fmap Character (char '\\' >> anyChar)

parseBool :: Parser LispVal
parseBool = fmap (Bool . (== 't')) (oneOf "ft")

parseFloat :: Parser LispVal
parseFloat = do whole <- many1 digit
                _ <- char '.'
                rest <- many1 digit
                return $ Float $ read (whole ++ "." ++ rest)

parseComplex :: Parser LispVal
parseComplex = do real <- parseNumber <|> parseFloat
                  _ <- char '+'
                  imag <- parseNumber <|> parseFloat
                  _ <- char 'i'
                  return $ Complex (toDouble real :+ toDouble imag)
              where toDouble (Float d) = d
                    toDouble (Number n) = fromInteger n
                    toDouble n = error $ "Expected a number, received: " ++ show n

parseRational :: Parser LispVal
parseRational = do nom <- parseNumber
                   _ <- char '/'
                   den <- parseNumber
                   return $ Rational (toNum nom % toNum den)
                where toNum (Number n) = n
                      toNum n = error $ "Expected a number, received: " ++ show n

-- Recursive type parsing
parseList :: Parser LispVal
parseList = fmap List $ sepBy parseExp spaces

parseDottedList :: Parser LispVal
parseDottedList = do car <- endBy parseExp spaces
                     cdr <- char '.' >> spaces >> parseExp
                     return $ DottedList car cdr

parseLists :: Parser LispVal
parseLists = do _ <- char '('
                x <- try parseList <|> parseDottedList
                _ <- char ')'
                return x

parseQuoted :: Parser LispVal
parseQuoted = char '\'' >> parseExp >>= (\e -> return $ List [Atom "quote", e])

parseExp :: Parser LispVal
parseExp = parseHashPrefix
           <|> parseAtom
           <|> parseString
           <|> try parseFloat
           <|> try parseComplex
           <|> try parseRational
           <|> parseNumber
           <|> parseQuoted
           <|> parseLists

readOrThrow :: Parser a -> String -> ThrowsError a
readOrThrow parser input = case parse parser "lisp" input of
  Left  err -> throwError $ Parser err
  Right val -> return val

readExpr = readOrThrow parseExp
readExprList = readOrThrow (endBy parseExp spaces)

max_stars_repo_name: kimsnj/writeAscheme

max_issues_repo_name: kimsnj/writeAscheme

max_forks_repo_name: kimsnj/writeAscheme