text:
{-# LANGUAGE GADTs #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE DeriveFunctor #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RecursiveDo #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE FunctionalDependencies #-}

module Dynamical.Sim.Internal where

import Control.Applicative
import Control.DeepSeq (force, NFData(..))
import Control.Monad
import Control.Monad.Fix (MonadFix(mfix))
import Control.Monad.State (StateT, State, get, put, runState, modify, execState, evalState)
import Control.Monad.Trans (lift)
import Control.Monad.Writer (Writer, tell, WriterT, runWriterT)
import Control.Parallel (pseq)
import qualified Control.Parallel.Strategies as Par
import Data.Complex (Complex((:+)))
import Data.Default.Class (Default(def))
import Data.Fixed (Fixed, mod')
import Data.IntMap.Strict (IntMap)
import qualified Data.IntMap.Strict as IntMap
import Data.IntSet (IntSet)
import qualified Data.IntSet as Set
import Data.List
import Data.Map (Map)
import qualified Data.Map as Map
import Data.Maybe
import Data.Monoid (Monoid(), Any(Any), (<>), mempty, mappend)
import Data.Proxy (Proxy(Proxy))
import Data.Scientific (Scientific)
import Data.StableMemo (memo)
import Data.Time (getCurrentTime, diffUTCTime, UTCTime)
import qualified Data.Vector as V
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Generic.Mutable as GM
import qualified Data.Vector.Mutable as VM
import qualified Data.Vector.Unboxed as UV
import qualified Data.Vector.Unboxed.Mutable as UM
import qualified Debug.Trace as Dbg
import GHC.Generics (Generic)
import qualified GHC.Prim as Prim
import GHC.TypeLits (Symbol, symbolVal, KnownSymbol)
import qualified Graphics.Rendering.Chart.Backend.Cairo as Chart
import qualified Graphics.Rendering.Chart.Easy as Chart
import Linear (V1(..), V2(..), angle, V3(..), V4(..), norm, _x, _y,_z,_w, qd,Additive(..),signorm,(*^),(^/), perp)
import System.IO.Unsafe (unsafePerformIO, unsafeInterleaveIO)
import System.Mem.StableName (hashStableName,eqStableName,StableName,makeStableName)
import Text.Printf (printf, PrintfArg)
import Unsafe.Coerce (unsafeCoerce)

-------------------------------------------------------------------
-- Signal
-------------------------------------------------------------------

-- | A `Signal t a` represents a time varying value of type `a`, where
-- time is measured using a `t`.
data Signal t a where
    SPure   :: a -> Signal t a
    SAp     :: Signal t (a -> b) -> Signal t a -> Signal t b
    SAdd    :: Num a => Signal t a -> Signal t a -> Signal t a
    SMul    :: Num a => Signal t a -> Signal t a -> Signal t a
    SDiv    :: Fractional a => Signal t a -> Signal t a -> Signal t a
    SMap    :: (a -> b) -> Signal t a -> Signal t b
    SInt    :: Splat t a => Int -> (a -> b) -> Signal t b
    SFn     :: Int -> (t -> a) -> Signal t a
    SSwitch :: Signal t a -> Event t (Sim t (Signal t a)) -> Signal t a
    SShare  :: Signal t a -> Signal t a

instance Monoid a => Monoid (Signal t a) where
    mempty = pure mempty
    mappend = liftA2 mappend

instance Functor (Signal t) where
    fmap f (SMap g s) = SMap (f . g) s
    fmap f (SInt i g) = SInt i (f . g)
    fmap f (SPure a) = SPure (f a)
    -- fmap f (SAp g a) = SAp (fmap (f.) g) a
    fmap f (SFn i g) = SFn i (f . g)
    -- fmap f (SSwitch s e) = SSwitch (fmap f s) (fmap (fmap (fmap f)) e)
    fmap f s = SMap f s

instance Applicative (Signal t) where
    pure = SPure
    (SPure f) <*> (SPure a) = pure (f a)
    (SPure f) <*> s = fmap f s
    a <*> b = SAp a b

dbgSig :: Signal t a -> String
dbgSig (SPure _) = "SPure"
dbgSig (SAp a b) = "SAp (" ++ dbgSig a ++ ") (" ++ dbgSig b ++ ")"
dbgSig (SMap _ s) = "SMap (" ++ dbgSig s ++ ")"
dbgSig (SInt i _) = "SInt " ++ show i
dbgSig (SFn i _) = "SFn " ++ show i
dbgSig (SSwitch s e) = "SSwitch (" ++ dbgSig s ++ ")"
dbgSig (SShare s) = "SShare (" ++ dbgSig s ++ ")"
dbgSig (SAdd a b) = "SAdd (" ++ dbgSig a ++ ") (" ++ dbgSig b ++ ")"
dbgSig (SDiv a b) = "SDiv (" ++ dbgSig a ++ ") (" ++ dbgSig b ++ ")"
dbgSig (SMul a b) = "SMul (" ++ dbgSig a ++ ") (" ++ dbgSig b ++ ")"

instance Num a => Num (Signal t a) where
    fromInteger = pure . fromInteger
    (SPure a) + (SPure b) = pure $ a + b
    a + b = SAdd a b
    (SPure a) * (SPure b) = pure $ a * b
    a * b = SMul a b
    abs = fmap abs
    signum = fmap signum
    negate = fmap negate

instance Fractional a => Fractional (Signal t a) where
    (SPure a) / (SPure b) = pure $ a / b
    a / b = SDiv a b
    recip = fmap recip
    fromRational = pure . fromRational

instance Floating a => Floating (Signal t a) where
    pi = pure pi
    exp = fmap exp
    log = fmap log
    sqrt = fmap sqrt
    (**) = liftA2 (**)
    logBase = liftA2 (**)
    sin = fmap sin
    cos = fmap cos
    tan = fmap tan
    asin = fmap asin
    acos = fmap acos
    atan = fmap atan
    sinh = fmap sinh
    cosh = fmap cosh
    tanh = fmap tanh
    asinh = fmap asinh
    acosh = fmap acosh
    atanh = fmap atanh

-------------------------------------------------------------------
-- Event
-------------------------------------------------------------------

-- | An `Event t a` represents a series of instantaneous events.
data Event t a where
    ERoot  :: (Num a, Ord a) => Signal t a -> Event t ()
    ETag   :: Event t a -> Signal t b -> Event t (a,b)
    EMap   :: (a -> b) -> Event t a -> Event t b

evalEvent :: Time t => Network t o -> Event t a -> a
evalEvent n (ERoot s) = ()
evalEvent n (EMap f e) = f $ evalEvent n e
evalEvent n (ETag e s) = (evalEvent n e, evalSignal n s)

instance Functor (Event t) where
    fmap = EMap

-------------------------------------------------------------------
-- Sim
-------------------------------------------------------------------

-- | The `Sim t` monad is used for constructing certain `Signal` values.
newtype Sim t a = Sim {unSim :: forall o. State (Network t o) a}

instance Num a => Num (Sim t a) where
    fromInteger = pure . fromInteger
    (+) = liftA2 (+)
    (*) = liftA2 (*)
    abs = fmap abs
    signum = fmap signum
    negate = fmap negate

instance Functor (Sim t) where
    fmap f (Sim s) = Sim (fmap f s)

instance Applicative (Sim t) where
    pure a = Sim (pure a)
    (Sim f) <*> (Sim a) = Sim (f <*> a)

instance Monad (Sim t) where
    (Sim a) >>= f = Sim $ do
        r <- a
        let (Sim k) = f r
        k

instance MonadFix (Sim t) where
    mfix f = Sim $ mfix (unSim . f)

-------------------------------------------------------------------
-- Network
-------------------------------------------------------------------

-- | This type family is used for determining the best way to store
-- a vector of values. It should evaluate to something that is an
-- instance of @Data.Vector.Generic.Vector@.
type family NetStoreVec a :: * -> *
type instance NetStoreVec Double = UV.Vector
type instance NetStoreVec Float = UV.Vector
type instance NetStoreVec Word = UV.Vector
type instance NetStoreVec Int = UV.Vector
type instance NetStoreVec (Fixed a) = V.Vector
type instance NetStoreVec Scientific = V.Vector

type NetStore t = NetStoreVec t t
type NetStoreMut t = G.Mutable (NetStoreVec t) t

type Time t = (NFData (NetStoreVec t t), Show t, G.Vector (NetStoreVec t) t, Show (NetStore t))

-- | A `Network` represents a compiled simulation.
data Network t o = Network
    { netIntState :: !(NetStore t)
    , netIntDeriv :: !(IntMap (Signal t (NetStore t)))
    , netFnTime   :: !(NetStore t)
    , netRoot     :: (Signal t o)  -- TODO: Make this a container of signals?
    } deriving (Generic)

-- TODO: Rewrite

-- | Cleans up a `Network` after a switching event, removing signals
-- that are no longer needed.
gc :: forall v t o. Time t => Network t o -> Network t o
gc n@Network{..} =
    let
        (intAll, intReachable, fnReachable) = findNames netRoot

        findNames :: Signal t a -> (IntSet, IntSet, IntSet)
        findNames s = execState (go s)  mempty

        go :: forall a. Signal t a -> State (IntSet, IntSet, IntSet) ()
        go (SInt ix (f :: x -> a)) = do
            (i,i',f) <- get
            when (not $ Set.member ix i) $ do
                let
                    pt = Proxy :: Proxy t
                    pa = Proxy :: Proxy x
                    len = splen pt pa - 1
                    names = Set.fromList [ix..ix + len]
                put (i `Set.union` names, Set.insert ix i', f)
                go (netIntDeriv IntMap.! ix)
        go (SAp a b) = go a >> go b
        go (SMap _ s) = go s
        go (SAdd a b) = go a >> go b
        go (SMul a b) = go a >> go b
        go (SDiv a b) = go a >> go b
        go (SPure _) = return ()
        go (SFn ix _) = modify $ mappend (mempty, mempty, Set.singleton ix)
        go (SSwitch s e) = go s >> goE e
        go (SShare s) = go s

        goE :: Event t a -> State (IntSet, IntSet, IntSet) ()
        goE (ERoot s) = go s
        goE (EMap _ e) = goE e
        goE (ETag e s) = goE e >> go s

        intMap  = UV.fromList (Set.toAscList intReachable)
        intMapAll = UV.fromList (Set.toAscList intAll)
        intMap' = IntMap.fromList $ zip (Set.toAscList intReachable) [0..]

        fnMap   = UV.fromList (Set.toAscList fnReachable)
        fnMap' = IntMap.fromList $ zip (Set.toAscList fnReachable) [0..]

        rename :: Signal t a -> Signal t a
        rename (SAp a b) = SAp (rename a) (rename b)
        rename (SMap f a) = SMap f (rename a)
        rename s@(SPure _) = s
        rename (SInt ix f) = SInt (intMap' IntMap.! ix) f
        rename (SFn ix f) = SFn (fnMap' IntMap.! ix) f
        rename (SSwitch s e) = SSwitch (rename s) (renameE e)
        rename (SShare s) = SShare (rename s)

        renameE :: Event t a -> Event t a
        renameE (ERoot s) = ERoot (rename s)
        renameE (EMap f e) = EMap f (renameE e)
        renameE (ETag e s) = ETag (renameE e) (rename s)

        intCnt = UV.length intMap
        fnCnt  = UV.length fnMap

        intState' = G.generate intCnt $ \ix -> netIntState G.! (intMapAll UV.! ix)

        restrictKeys m s = IntMap.filterWithKey (\k _ -> k `Set.member` s) m
        intDeriv' = fmap rename $ restrictKeys netIntDeriv intReachable

        fnTime' = G.generate fnCnt $ \ix -> netFnTime G.! (fnMap UV.! ix)
        root' = rename netRoot

        -- intDeriv' = execWriter $ getDeriv root'

    in n
        { netIntState = intState'
        , netIntDeriv = IntMap.mapKeys (\k -> intMap' IntMap.! k) intDeriv'
        , netFnTime = fnTime'
        , netRoot = root'
        }

eventOccured :: Time t => Network t o -> Network t o -> Event t a -> Maybe a
eventOccured old new (EMap f e) = f <$> eventOccured old new e
eventOccured old new (ETag e s) = flip (,) (evalSignal new s) <$> eventOccured old new e
eventOccured old new (ERoot s)
    | oldS < 0 && 0 <= newS = Just ()
    | newS <= 0 && 0 < oldS = Just ()
    | otherwise = Nothing
        where
            oldS = evalSignal old s
            newS = evalSignal new s

runSwitches' :: forall v t o. Time t => Network t o -> Network t o -> (Bool, Network t o)
runSwitches' old new =
    let
        go :: Signal t a -> WriterT Any (Sim t) (Signal t a)
        go s@(SPure a) = return s
        go (SMap f s) = go s >>= pure . SMap f
        go (SAdd a b) = do
            a' <- go a
            b' <- go b
            return $ SAdd a' b'
        go (SMul a b) = do
            a' <- go a
            b' <- go b
            return $ SMul a' b'
        go (SDiv a b) = do
            a' <- go a
            b' <- go b
            return $ SDiv a' b'
        go (SAp f a) = do
            f' <- go f
            a' <- go a
            return $ SAp f' a'
        go s@(SInt _ _) = return s
        go s@(SFn ix f) = return s
        go (SSwitch s e) = case eventOccured old new e of
            Just v -> tell (Any True) >> lift v
            Nothing -> do
                s' <- go s
                e' <- goE e
                return $ SSwitch s' e'
        go (SShare s) = go s >>= pure . SShare

        goE :: Event t a -> WriterT Any (Sim t) (Event t a)
        goE (ERoot s) = go s >>= pure . ERoot
        goE (EMap f e) = goE e >>= pure . EMap f
        goE (ETag e s) = do
            e' <- goE e
            s' <- go s
            return (ETag e' s')

        (((root', derivs'), Any changed), net') = addSim new $ runWriterT $ do
            newRoot <- go $ netRoot new
            newDerivs <- mapM go $ netIntDeriv new
            return (newRoot, newDerivs)
    in
        if changed
        then (changed, gc $ net'
            { netRoot = root'
            , netIntDeriv = derivs'
            })
        else (changed, new)

runSwitches old new = snd $ runSwitches' old new
anyEvent old new = fst $ runSwitches' old new

newDoubleSim :: Sim Double (Signal Double a) -> Network Double a
newDoubleSim = newSim

newSim :: Time t => Sim t (Signal t a) -> Network t a
newSim (Sim s) =
    let
        (r,n) = runState s Network
            { netIntState = G.empty
            , netIntDeriv = IntMap.empty
            , netFnTime  = G.empty
            , netRoot = r
            }
    in n

addSim :: Time t => Network t o -> Sim t a -> (a, Network t o)
addSim n (Sim s) = runState s n

-- TODO: Below are three implementations of evalSignal. Each with different
-- sharing characteristics.
--  * evalSignal' does explicit sharing
--  * evalSignal'' does no sharing
--  * evalSignal''' does implicit sharing
-- We need to benchmark and determine which is best.
-- Once we determine which solution to go for, we should integrate it
-- more with the rest of the code such that sharing caries over to
-- events.
--
-- It's probably worth redefining `share` to `id` when not using sharing
-- or using implicit sharing.
evalSignal'
    :: forall t o a. (Time t)
    => Network t o -> Signal t a -> State (IntMap (StableName Prim.Any, Any)) a
evalSignal' c me@(SShare s) = do
    m <- get
    let
        name = unsafePerformIO $ makeStableName me
        nameh = hashStableName name
        mv = do
            (n,v) <- IntMap.lookup nameh m
            if eqStableName name n
            then return $ unsafeCoerce v
            else Nothing
    case mv of
        Just v -> return v
        Nothing -> do
            v' <- evalSignal' c s
            put $ IntMap.insert nameh (unsafeCoerce name, unsafeCoerce v') m
            return v'
evalSignal' c (SPure a) = return a
evalSignal' c (SAdd a b) = do
    a' <- evalSignal' c a
    b' <- evalSignal' c b
    return $ a' + b'
evalSignal' c (SMul a b) = do
    a' <- evalSignal' c a
    b' <- evalSignal' c b
    return $ a' * b'
evalSignal' c (SDiv a b) = do
    a' <- evalSignal' c a
    b' <- evalSignal' c b
    return $ a' / b'
evalSignal' c (SMap f s) = f <$> evalSignal' c s
evalSignal' c (SAp f a) = do
    f' <- evalSignal' c f
    a' <- evalSignal' c a
    return $ f' a'
evalSignal' c (SFn ix f) = return $ f $ netFnTime c G.! ix
evalSignal' c (SInt ix (f :: x -> a)) = do
    let
        len = splen (Proxy :: Proxy t) (Proxy :: Proxy x)
    return $ f $ unsplat $ G.slice ix len (netIntState c)
evalSignal' c (SSwitch s _) = evalSignal' c s

evalSignal'' :: forall t o a. Time t => Network t o -> Signal t a -> a
evalSignal'' c = go
    where
        go :: forall b. Signal t b -> b
        go (SPure a) = a `seq` a
        go (SMap f s) = let x = f $! go s in x `seq` x
        go (SAdd a b) = let x = go a + go b in x `seq` x
        go (SMul a b) = let x = go a * go b in x `seq` x
        go (SDiv a b) = let x = go a / go b in x `seq` x
        go (SAp f a) =
            let
                f' = go f
                a' = go a
                x = f' a'
            in x `seq` x
        go (SFn ix f) = let x = f $! netFnTime c G.! ix in x `seq` x
        go (SInt ix (f :: x -> b)) = let x = f $! unsplat $ G.slice ix (splen (Proxy :: Proxy t) (Proxy :: Proxy x)) (netIntState c) in x `seq` x
        go (SSwitch s _) = go s
        go (SShare s) = go s

evalSignal''' :: forall t o a. Time t => Network t o -> Signal t a -> a
evalSignal''' c = go
    where
        go :: Signal t c -> c
        go = memo gogo

        gogo :: forall b. Signal t b -> b
        gogo (SPure a) = a
        gogo (SMap f s) = f $ go s
        gogo (SAdd a b) = go a + go b
        gogo (SMul a b) = go a * go b
        gogo (SDiv a b) = go a / go b
        gogo (SAp f a) = go f $ go a
        gogo (SFn ix f) = f $ netFnTime c G.! ix
        gogo (SInt ix (f :: x -> b)) = f $ unsplat $ G.slice ix (splen (Proxy :: Proxy t) (Proxy :: Proxy x)) (netIntState c)
        gogo (SSwitch s _) = go s
        gogo (SShare s) = go s

evalSignal :: forall t o a. Time t => Network t o -> Signal t a -> a
--evalSignal n s = evalState (evalSignal' n s) mempty
evalSignal = evalSignal''
--evalSignal = evalSignal'''
--evalSignal n = memo (evalSignal'' n)
--{-# NOINLINE evalSignal #-}

evalRoot :: Time t => Network t o -> o
evalRoot c = evalSignal c (netRoot c)

deltaNet :: (Num t, Time t) => Network t o -> t -> NetStore t -> Network t o
deltaNet n@Network{..} h dv = deltaNetState (deltaNetTime n h) dv

deltaNetState :: (Num t, Time t) => Network t o -> NetStore t -> Network t o
deltaNetState n@Network{..} dv = n
    { netIntState = G.zipWith (+) netIntState dv
    }

deltaNetTime :: (Num t, Time t) => Network t o -> t -> Network t o
deltaNetTime n@Network{..} h = n
    { netFnTime  = G.map (+h) netFnTime
    }

derivs :: (Num t, Time t) => Network t o -> t -> NetStore t -> NetStore t
derivs n h dv =
    let
        n' = deltaNet n h dv
    in derivsNow n'

derivsNow :: (Time t) => Network t o -> NetStore t
derivsNow n = G.concat $ parMapChunk (evalSignal n . snd) $ IntMap.toAscList $ netIntDeriv n

-- TODO: Why is this burning so much CPU and not giving me /that/ much of
-- a speed up. Changing the chunk size doesn't seem to make that much
-- difference.
--
-- 1 -N1  (100% cpu)
-- real	0m14.623s
-- user	0m14.458s
-- sys	0m0.149s
--
-- All the rest 580% cpu
-- 1 -N
-- real	0m10.852s
-- user	0m39.608s
-- sys	0m22.108s
--
-- 20 -N
-- real	0m10.932s
-- user	0m35.629s
-- sys	0m23.757s
--
-- 200 -N
-- real	0m11.231s
-- user	0m37.745s
-- sys	0m24.182s
parMapChunk f as =
    let
        a = map f as
    in a `Par.using` Par.parListChunk 50 (Par.rdeepseq)

-------------------------------------------------------------------
-- Splat
-------------------------------------------------------------------

class Splat t a where
    splen :: Proxy t -> Proxy a -> Int
    splat :: Time t => a -> NetStore t
    unsplat :: Time t => NetStore t -> a

instance Splat t t where
    splen _ _ = 1
    splat t = G.singleton t
    unsplat v = v G.! 0

instance Splat t (t,t) where
    splen t _ = 2
    splat (a,b) = splat a G.++ splat b
    unsplat v = (v G.! 0, v G.! 1)

instance Splat t (t,t,t) where
    splen t _ = 3
    splat (a,b,c) = G.fromList [a,b,c]
    unsplat v = (v G.! 0, v G.! 1, v G.! 2)

instance Splat t (t,t,t,t) where
    splen t _ = 4
    splat (a,b,c,d) = G.fromList [a,b,c,d]
    unsplat v = (v G.! 0, v G.! 1, v G.! 2, v G.! 3)

instance Splat t (Complex t) where
    splen t _ = 2
    splat (a :+ b) = G.fromList [a,b]
    unsplat v = (v G.! 0) :+ (v G.! 1)

instance Splat t (V1 t) where
    splen t _ = splen t (Proxy :: Proxy t)
    splat (V1 a) = splat a
    unsplat = V1 . unsplat

instance Splat t (V2 t) where
    splen t _ = splen t (Proxy :: Proxy (t,t))
    splat (V2 a b) = splat (a,b)
    unsplat v = let (a,b) = unsplat v in V2 a b

instance Splat t (V3 t) where
    splen t _ = splen t (Proxy :: Proxy (t,t,t))
    splat (V3 a b c) = splat (a,b,c)
    unsplat v = let (a,b,c) = unsplat v in V3 a b c

instance Splat t (V4 t) where
    splen t _ = splen t (Proxy :: Proxy (t,t,t,t))
    splat (V4 a b c d) = splat (a,b,c,d)
    unsplat v = let (a,b,c,d) = unsplat v in V4 a b c d

-------------------------------------------------------------------
-- API
-------------------------------------------------------------------

-- | Integrate the input `Signal` with respect to time.
integral :: (Splat t a, Time t) => a -> Signal t a -> Sim t (Signal t a)
integral i s = Sim $ do
    st <- get
    let
        sp = splat i
        ix = G.length (netIntState st)
    put st
        { netIntState = netIntState st G.++ sp
        , netIntDeriv = IntMap.insert ix (fmap splat s) (netIntDeriv st)
        }
    return $ SInt ix id

-- | Create a `Signal` that is a pure function of time. Each `timeFn'` has
-- a local concept of time, and this time starts from the provided time
-- value.
timeFn' :: Time t => t -> (t -> a) -> Sim t (Signal t a)
timeFn' t f = Sim $ do
    st <- get
    let
        ix = G.length (netFnTime st)
    put st
        { netFnTime = netFnTime st `G.snoc` t
        }
    return $ SFn ix f

-- | As `timeFn'` but local time starts at 0.
timeFn :: (Time t, Num t) => (t -> a) -> Sim t (Signal t a)
timeFn = timeFn' 0

-- | Start out as the input `Signal`, but when the `Event` occurs, become
-- the signal that it carries, and remain as that signal forever. Future
-- `Event`s from this event source will be ignored. See `switch` for
-- alternate behavior.
become :: Signal t a -> Event t (Sim t (Signal t a)) -> Signal t a
become = SSwitch

-- | Start out as the input `Signal` and everytime the `Event` occurs, change
-- to the new `Signal` carried in the `Event`. This can easily cause
-- a space leak if used recursively. Look at `become` for a safer solution.
switch :: Signal t a -> Event t (Sim t (Signal t a)) -> Signal t a
switch s e = become s (fmap (\s' -> s' >>= \s'' -> pure (switch s'' e)) e)

-- | Useful combination of `become` and `fmap`
becomeOn :: Signal t a -> Event t b -> (b -> Sim t (Signal t a)) -> Signal t a
becomeOn s e f = become s (fmap f e)

-- | Turns the input `Signal t a` into a `Signal t (Maybe a)`, resulting
-- in `Nothing` as soon as the `Event` fires, and remaining `Nothing` forever
-- after that point.
becomeNothingOn :: Signal t a -> Event t b -> Signal t (Maybe a)
becomeNothingOn s e = becomeOn (fmap Just s) e $ \_ -> return $ pure Nothing

-- | Memoize the result of evaluating this Signal so that repeated
-- uses don't have to re-evaluate everything. Integrations are
-- already shared by default, so avoid sharing them directly.
--
-- TODO: Depending on the signal evaluation strategy (see comments in code)
-- sharing may implicitly happen for all signals, or not at all.
share :: Signal t a -> Signal t a
share = SShare

-- | Emit an event every time the input `Signal` crosses 0.
root :: (Ord t, Num t) => Signal t t -> Event t ()
root = ERoot

-- | Tag the given event with the value of the signal when it occurs.
tag :: Event t a -> Signal t b -> Event t (a,b)
tag  = ETag

-------------------------------------------------------------------
-- Simulation
-------------------------------------------------------------------

data SimResult t o = SimResult
    { stepNumber :: !Integer
    , globalTime :: !t
    , result     :: !o
    } deriving (Eq, Ord, Show, Functor)

asTuple :: SimResult t o -> (t,o)
asTuple (SimResult _ t o) = (t,o)

-- | Run the `Integrator` over the given `Signal` outputting a list of
-- time stamp and `Network`s.
simulate'
    :: forall v t o. (Num t, Time t)
    => Integrator t o -> Sim t (Signal t o) -> [SimResult t (Network t o)]
simulate' integrator s =
    let
        n = newSim s
    in unfoldr (\(cnt,i,n,t) ->
        let
            (n',r,i') =  runIntegrator i n
            t' = t + r
        in
            cnt `seq` t `seq` n  `seq` (Just $!
                ( SimResult cnt t n   -- Previous step
                , (succ cnt, i', n', t')
                ))
        ) (0,integrator,n,0)

-- | Run the `Integrator` over the given `Signal`, outputting a list of
-- time stamp and signal values.
simulate :: (Num t, Time t) => Integrator t o -> Sim t (Signal t o) -> [SimResult t o]
simulate i s = fmap evalRoot <$> simulate' i s

-- | Causes progress messages to be printed every 100 steps. Uses lazy IO,
-- so be careful.
--  TODO: Rejig things to be more streamy.
trace :: (Show t, Show o, PrintfArg t) => [SimResult t o] -> IO [SimResult t o]
trace [] = return []
trace (s@(SimResult c t o):rest)
    | 0 == c `mod` 100 = do
        printf "S: %5d  T: %5.3v  V: " c t
        putStrLn $ show o
        r <- unsafeInterleaveIO $ trace rest
        return $ s : r
    | otherwise = do
        r <- unsafeInterleaveIO $ trace rest
        return $ s : r

-- | `simulate'` specialised to `Double`
simulateDouble'
    :: Integrator Double o
    -> Sim Double (Signal Double o)
    -> [SimResult Double (Network Double o)]
simulateDouble' = simulate'

-- | `simulate` specialised to `Double`
simulateDouble
    :: Integrator Double o
    -> Sim Double (Signal Double o)
    -> [SimResult Double o]
simulateDouble = simulate

-- | `simulate` which stops as soon as the `Signal` value becomes `Nothing`
simulateJust :: (Num t, Time t) => Integrator t (Maybe o) -> Sim t (Signal t (Maybe o)) -> [SimResult t o]
simulateJust i s = fmap (\(SimResult n t (Just x)) -> SimResult n t x) $ takeWhile (isJust . result) $ simulate i s

-- | `simulateJust` specialised to `Double`
simulateJustDouble :: Integrator Double (Maybe o) -> Sim Double (Signal Double (Maybe o)) -> [SimResult Double o]
simulateJustDouble = simulateJust

-- | Simulate until the specified time is reached. The termination time is
-- marked by an `Event`, meaning that adaptive integration methods which
-- narrow down on event locations will be do so to get accurate finishing
-- times.
simulateUntil :: (Ord t, Num t, Time t) => t -> Integrator t (Maybe o) -> Sim t (Signal t o) -> [SimResult t o]
simulateUntil duration i s = simulateJust i $ do
    t <- timer duration
    s' <- s
    return $ s' `becomeNothingOn` t

-- | `simulateUntil` specialised to `Double`
simulateUntilDouble :: Double -> Integrator Double (Maybe o) -> Sim Double (Signal Double o) -> [SimResult Double o]
simulateUntilDouble = simulateUntil

-- | Run a simulation in real time, printing the result values to the terminal.
runRk4RealTime :: forall a. Show a => Sim Double (Signal Double a) -> IO ()
runRk4RealTime s =
    let
        n = newSim s
        go :: UTCTime -> Network Double a -> IO ()
        go prev n = do
            now <- getCurrentTime
            let
                h = realToFrac $ diffUTCTime now prev
                (n',_,_) = runIntegrator (rk4 h) n
            print $ evalRoot n'
            go now n'
    in do
        now <- getCurrentTime
        go now n

-- | Same as `runRk4RealTime` but stops when a `Nothing` value is produced.
runRk4RealTimeJust :: forall a. Show a => Sim Double (Signal Double (Maybe a)) -> IO ()
runRk4RealTimeJust s =
    let
        n = newSim s
        go :: UTCTime -> Network Double (Maybe a) -> IO ()
        go prev n = do
            now <- getCurrentTime
            let
                h = realToFrac $ diffUTCTime now prev
                (n',_,_) = runIntegrator (rk4 h) n
            case evalRoot n' of
                Nothing -> return ()
                Just x -> print x >> go now n'
    in do
        now <- getCurrentTime
        go now n

sample :: (Ord t, Real t) => t -> [(t,a)] -> [(t,a)]
sample freq = go 500
    where
        go _ [] = []
        go n ((t,r):xs)
            | t' < n = (t,r) : go t' xs
            | otherwise = go t' xs
            where
                t' = t `mod'` freq

-------------------------------------------------------------------
-- Integrators
-------------------------------------------------------------------

newtype Integrator t o = Integrator
    { runIntegrator :: Network t o -> (Network t o, t, Integrator t o)
    } deriving Generic

instance NFData (Integrator t o)

eulerSimple :: (Num t, Time t) => t -> Integrator t o
eulerSimple h = integrator
    where
        integrator = Integrator go
        go n =
            let
                ds = derivsNow n
                n' = deltaNet n h (G.map (*h) ds)
            in (runSwitches n n', h, integrator)

-- | Same as `euler` except that it bisects around `Event`s until it is within
-- the provided tolerance.
eulerBisect :: (Fractional t, Ord t, Time t) => t -> t -> Integrator t o
eulerBisect tol h = integrator
    where
        integrator = Integrator $ go h
        go h' n =
            let
                ds = derivsNow n
                n' = deltaNet n h' (G.map (*h') ds)
                (e,n'') = runSwitches' n n'
                res
                    -- TODO: Don't emit nearby points, jump strait to it.
                    | e && h' > tol = go (h'/2) n
                    | otherwise = (n'', h', integrator)
            in res

-- | Configuration for adaptive methods.
--
-- Unchecked assumptions: 0 < `minStep` < `eventTolerance` < `maxStep`
data AdaptiveConfig t = AdaptiveConfig
    { minStep :: Maybe t
    , maxStep :: Maybe t
    , tolerance :: t
    , eventTolerance :: t
    }

instance Fractional t => Default (AdaptiveConfig t) where
    def = AdaptiveConfig
        { minStep = Nothing
        , maxStep = Nothing
        , tolerance = 0.001
        , eventTolerance = 0.001
        }

-- If maxStep < minStep, we choose maxStep always.
adaptiveClamp :: Ord t => AdaptiveConfig t -> t -> t
adaptiveClamp AdaptiveConfig{..} t =
    let
        t' = maybe t (max t) minStep
    in
        maybe t' (min t') maxStep

clamp :: Ord a => a -> a -> a -> a
clamp l h t
    | t < l = l
    | t > h = h
    | otherwise = t

euler :: (Ord t, Fractional t, Time t, Show t) => AdaptiveConfig t -> Integrator t o
euler conf@AdaptiveConfig{..} = Integrator $ go 0.1
    where
        go h n =
            let
                eulStep h n = deltaNet n h (G.map (*h) (derivsNow n))

                n10  = eulStep h n
                n1_2 = eulStep (h/2) n
                n11  = eulStep (h/2) n1_2
                tau  = G.zipWith (-) (netIntState n11) (netIntState n10)
                tau' = G.maximum $ G.map abs tau
                n12  = deltaNetState n11 tau

                h' = adaptiveClamp conf $ 0.9 * h * clamp 0.3 3 (tolerance / tau')

                (e,n') = runSwitches' n n12

                res
                    | e && h > eventTolerance && h > fromMaybe 0 minStep
                        = go (adaptiveClamp conf (h/2)) n
                    | tau' > tolerance && h' > fromMaybe 0 minStep = go h' n
                    | otherwise = (n', h, Integrator $ go h')
            in res

rk4 :: (Fractional t, Time t) => t -> Integrator t o
rk4 h = integrator
    where
        integrator = Integrator go
        go n =
            let
                h2 = h / 2
                h6 = h / 6
                k1 = derivsNow n
                k2 = derivs n h2 (G.map (*h2) k1)
                k3 = derivs n h2 (G.map (*h2) k2)
                k4 = derivs n h  (G.map (*h)  k3)
                n' = deltaNet n h (G.map (* h6)
                    $ G.zipWith4 (\a b c d -> a + 2 * b + 2 * c + d) k1 k2 k3 k4)
            in
                (runSwitches n n', h, integrator)

-------------------------------------------------------------------
-- Visualisation
-------------------------------------------------------------------

class Plot a where
    plot :: [SimResult Double a] -> Chart.EC (Chart.Layout Double Double) ()

newtype a ::: (name :: Symbol) = Named a
    deriving (Read,Enum,Real,RealFrac,RealFloat,Eq,Ord,Show,Num,Fractional,Floating,Chart.PlotValue)

infixr 0 :::

instance Plot Double where
    plot = Chart.plot . Chart.points "" . fmap asTuple

instance (KnownSymbol name, Chart.PlotValue a) => Plot (a ::: name) where
    plot = Chart.plot . Chart.points name . fmap (fmap Chart.toValue) . fmap asTuple
        where
            name = symbolVal (Proxy :: Proxy name)

instance (Plot a, Plot b) => Plot (a,b) where
    plot as = do
        plot (fmap fst <$> as)
        plot (fmap snd <$> as)

instance (Plot a, Plot b, Plot c) => Plot (a,b,c) where
    plot as = do
        plot (fmap (\(a,_,_) -> a) <$> as)
        plot (fmap (\(_,b,_) -> b) <$> as)
        plot (fmap (\(_,_,c) -> c) <$> as)

instance Plot a => Plot (V2 a) where
    plot as = do
        plot (fmap (Chart.^. _x) <$> as)
        plot (fmap (Chart.^. _y) <$> as)

instance Plot a => Plot (V4 a) where
    plot as = do
        plot (fmap (Chart.^. _x) <$> as)
        plot (fmap (Chart.^. _y) <$> as)
        plot (fmap (Chart.^. _z) <$> as)
        plot (fmap (Chart.^. _w) <$> as)

plotSimUntil
    :: Plot o
    => Double
    -> Integrator Double (Maybe o)
    -> Sim Double (Signal Double o)
    -> Chart.EC (Chart.Layout Double Double) ()
plotSimUntil t i s = plot (simulateUntil t i s)

plotSimUntilToFile n title t i s = do
    r <- trace $ simulateUntil t i s
    Chart.toFile Chart.def n $ do
        Chart.layout_title Chart..= title
        plot r

class PlotPara a where
    plotPara :: [SimResult Double a] -> Chart.EC (Chart.Layout Double Double) ()

instance PlotPara (V2 Double) where
    plotPara = Chart.plot . Chart.line "" . (:[]) . fmap (\(SimResult n t (V2 a b)) -> (a,b))

instance (KnownSymbol name) => PlotPara (V2 Double ::: name) where
    plotPara = Chart.plot . Chart.line name . (:[]) . fmap (\(SimResult n t (Named (V2 a b))) -> (a,b))
        where
            name = symbolVal (Proxy :: Proxy name)

instance (PlotPara a, PlotPara b) => PlotPara (a,b) where
    plotPara as = do
        plotPara (fmap fst <$> as)
        plotPara (fmap snd <$> as)

instance (PlotPara a, PlotPara b, PlotPara c) => PlotPara (a,b,c) where
    plotPara as = do
        plotPara (fmap (\(a,_,_) -> a) <$> as)
        plotPara (fmap (\(_,b,_) -> b) <$> as)
        plotPara (fmap (\(_,_,c) -> c) <$> as)

instance PlotPara (Map String (V2 Double)) where
    plotPara [] = return ()
    plotPara as@(SimResult _ _ a : _) = do
        let
            ks = Map.keys a
            kk = if length ks > 10 then \_ -> "" else id
        forM_ a $ \(V2 x y) -> do
            Chart.plot $ Chart.points "" [(x,y)]
        forM_ ks $ \k -> do
            Chart.plot $ Chart.line (kk k) . (:[]) $ fmap (\(SimResult n t m) ->
                let
                    V2 a b = m Map.! k
                in (a,b)
                ) as

plotParaSimUntil
    :: PlotPara o
    => Double
    -> Integrator Double (Maybe o)
    -> Sim Double (Signal Double o)
    -> Chart.EC (Chart.Layout Double Double) ()
plotParaSimUntil t i s = plotPara $ simulateUntil t i s

plotParaSimUntilToFile n title t i s = do
    r <- trace $ simulateUntil t i s
    Chart.toFile Chart.def n $ do
        Chart.layout_title Chart..= title
        plotPara r

-------------------------------------------------------------------
-- Utility
-------------------------------------------------------------------

-- | Applicative tupples
(-:) :: Applicative f => f a -> f b -> f (a,b)
(-:) = liftA2 (,)
infixr 5 -:

-- | Applicative `Named` tuples.
(-::) :: Applicative f => f a -> f b -> f (a ::: an, b ::: bn)
a -:: b = (,) <$> (Named <$> a) <*> (Named <$> b)
infixr 5 -::

threeNames :: Applicative f => f a -> f b -> f c -> f (a ::: an, b ::: nm, c ::: cn)
threeNames a b c = (,,) <$> (Named <$> a) <*> (Named <$> b) <*> (Named <$> c)

-------------------------------------------------------------------
-- Examples/Test
-------------------------------------------------------------------

-- Ex1 demonstrates basic integration with recursion.
ex1 :: (Time t, Fractional t) => Sim t (Signal t (t, t))
ex1 = do
    rec a <- integral 0.5 a
    b <- integral 0 1
    return $ a -: b

-- Ex2 shows how to use `become` to change a singal to a new one upon
-- an event.
ex2 :: forall t. (Fractional t, Ord t, Time t) => Sim t (Signal t t)
ex2 = do
    t <- integral 5 (-1)
    let
        e = root t
        e' = fmap (const $ pure 3) e
    return $ become 1 e'

-- Ex3 show the use of `switch` to control a signal.
ex3 :: (Time t, Floating t, Ord t) => Sim t (Signal t (t,t))
ex3 = do
    t <- timeFn sin
    let
        e = tag (root t) (signum t)
        e' = fmap (\(_,v) -> return $ pure v) e
    return $ (,) <$> t <*> switch 0 e'

-- Ex4 can be used to demonstrate the accuracy of various integrators.
ex4
    :: (Floating t, Time t)
    => Sim t (Signal t
        ( t ::: "sin(t)"
        , t ::: "∫∫∫∫sin(t)"
        ))
ex4 = do
    s <- timeFn sin
    s' <- integral (-1) s >>= integral 0 >>= integral 1 >>= integral 0
    return $ s -:: s'

-- Ex5 demonstrates using `Maybe` to terminate a simulation with
-- `simulateJust`.
--
-- Note: The use of an `Event` rather than `fmap`ing over the signal
-- causes the root finding algorithm of many integrators to kick in allowing
-- for very precise finishing times.
ex5 :: (Floating t, Time t, Ord t) => Sim t (Signal t (Maybe (V4 t)))
ex5 = do
    i <- (integral 0 $ pure (V4 0 1 2 (-0.5)) )
    let
        s = fmap (\i -> 5 - norm i) i
        e = root s
    return $ i `becomeNothingOn` e

-- Ex6 demonstrates the use of `share` to reduce the number
-- of times a signal is evaluated. The trace will show
-- "Calc A" once, but "Calc B" twice when using an explicit
-- sharing eval strategy.
ex6 :: (Time t, Floating t) => Sim t (Signal t (V4 t))
ex6 = do
    let
        a = share $ fmap (Dbg.trace "Calc A") $ 5 + 2
        b = fmap (Dbg.trace "Calc B") $ 5 + 2
    return $ V4 <$> a <*> a <*> b <*> b

ex7 :: (Ord t, Floating t, Time t) => Sim t (Signal t t)
ex7 = do
    s <- timeFn sin
    let
        e = root s
        s1 = becomeOn 0 e (\_ -> pure $ signum s)
    integral 0 s1

timer :: (Num t, Time t, Ord t) => t -> Sim t (Event t ())
timer duration = timeFn (\t -> duration - t) >>= pure . root

stepAndHold :: (Num t, Ord t, Time t) => [(t,Sim t (Signal t a))] -> Sim t (Signal t a) -> Sim t (Signal t a)
stepAndHold [] a = a
stepAndHold ((t,a):ls) last = do
    e <- timer t
    s <- a
    return $ s `becomeOn` e $ \_ -> stepAndHold ls last

ex8 :: (Ord t, Floating t, Time t) => Sim t (Signal t (Maybe (t ::: "Deriv",t ::: "Integral")))
ex8 = do
    t <- timer 6
    d <- stepAndHold [(0.5,1),(0.2, -2),(1,0)] (timeFn sin)
    i <- integral 0 d
    let r = d -:: i
    return $ r `becomeNothingOn` t

ex9
    :: forall t. (Floating t, Time t)
    => Sim t (Signal t
        ( V2 t ::: "Sun"
        , V2 t ::: "Earth"
        , V2 t ::: "Moon"
        ))
ex9 = do
    rec
        let
            degRad :: t -> t
            degRad t = t / 180 * pi

            pos2init = 149.6e9 *^ angle $ degRad 3.201e2
            pos3init = pos2init + 384.4e6 *^ angle (degRad 2.0644e2)

        pos1 <- integral (V2 0 0) vel1
        pos2 <- integral pos2init vel2
        pos3 <- integral pos3init vel3

        let
            vel2init = 29766.42101876582 *^ signorm (perp pos2init)
            vel3init = vel2init + 1023.005 *^ (signorm (perp pos3init))
        vel1 <- integral (V2 0 0) acc1
        vel2 <- integral vel2init acc2
        vel3 <- integral vel3init acc3

        let
            g = 6.674e-11
            m1 = 1.9891e30
            m2 = 5.972e24
            m3 = 7.347e22
            -- f = g * m1 * m2 / r^2 = m1 * a
            dir :: Signal t (V2 t) -> Signal t (V2 t) -> Signal t t -> Signal t (V2 t)
            dir a b c =
                let
                    d = (^-^) <$> b <*> a
                    u = signorm <$> d
                in (*^) <$> c <*> u
            f1 = sum
                [ dir pos1 pos2 $ g * m1 * m2 / (qd <$> pos1 <*> pos2)
                , dir pos1 pos3 $ g * m1 * m3 / (qd <$> pos1 <*> pos3)
                ]
            f2 = sum
                [ dir pos2 pos1 $ g * m2 * m1 / (qd <$> pos2 <*> pos1)
                , dir pos2 pos3 $ g * m2 * m3 / (qd <$> pos2 <*> pos3)
                ]
            f3 = sum
                [ dir pos3 pos1 $ g * m3 * m1 / (qd <$> pos3 <*> pos1)
                , dir pos3 pos2 $ g * m3 * m2 / (qd <$> pos3 <*> pos2)
                ]
            acc1 = (^/) <$> f1 <*> m1
            acc2 = (^/) <$> f2 <*> m2
            acc3 = (^/) <$> f3 <*> m3
    return $ threeNames pos1 pos2 pos3

max_stars_repo_name: luke-clifton/dynamical

max_issues_repo_name: luke-clifton/dynamical

max_forks_repo_name: luke-clifton/dynamical


text:
{-# OPTIONS_GHC -fno-warn-type-defaults #-}
{-# LANGUAGE DataKinds, FlexibleContexts, FlexibleInstances, GADTs        #-}
{-# LANGUAGE MultiParamTypeClasses, NoMonomorphismRestriction             #-}
{-# LANGUAGE OverloadedStrings, PolyKinds, TypeApplications, TypeFamilies #-}
{-# LANGUAGE UndecidableInstances                                         #-}
{-# OPTIONS_GHC -fno-warn-type-defaults -fno-warn-orphans -Wall #-}
module Main where
import Cases

import           Algebra.Bridge.Singular
import           Algebra.Field.Prime
import           Algebra.Prelude.Core
import           Control.Monad                   (void)
import qualified Data.Text                       as T
import qualified Data.Vector.Unboxed             as V
import           Prelude                         (read)
import           Statistics.Resampling
import           Statistics.Resampling.Bootstrap
import           Statistics.Types
import qualified System.Random.MWC               as Rand

benchIdeal :: IsSingularPolynomial poly => Text -> Ideal poly -> IO Double
benchIdeal fun i = fmap ((/1000) . read . T.unpack) $ singular $
  toProg fun i

toProg :: (SingularOrder (Arity poly) (MOrder poly), SingularCoeff (Coefficient poly), IsOrderedPolynomial poly) => Text -> Ideal poly -> Text
toProg fun i =
  let expr = funE fun [idealE' i]
  in prettySingular $ do
    directC "system(\"--ticks-per-sec\", 1000000)"
    void $ ringC "R" expr
    declOnlyC IdealT "G"
    directC "timer=1"
    directC "int t=rtimer"
    letC "G" expr
    directC "print(rtimer-t)"
    directC "exit"

analyse :: (IsSingularPolynomial poly) => String -> Text -> Ideal poly -> IO ()
analyse lab fun ideal = do
  gen <- Rand.create
  i2Gr <- V.replicateM 50 $ benchIdeal fun ideal
  res  <- resample gen [Mean, StdDev] 1000 i2Gr
  let [Estimate mn (ConfInt lbmn ubmn _) ,Estimate dv _]
        = bootstrapBCA cl95 i2Gr res
  putStrLn lab
  mapM_ (putStrLn . ('\t':))
    ["Mean:\t" ++ show mn ++ "(ms)"
    ,"MeanLB:\t" ++ show lbmn
    ,"MeanUB:\t" ++ show ubmn
    ,"StdDev:\t" ++ show dv
    ]

runTestCases :: (IsMonomialOrder n o, KnownNat n)
             => String -> Ideal (OrderedPolynomial Rational o n) -> IO ()
runTestCases lab i = do
  analyse (lab ++ " (Q, Lex, Sing(gr))") "groebner" $ fmap (changeOrder Lex) i
  analyse (lab ++ " (Q, Lex, Sing(sba))") "sba" $ fmap (changeOrder Lex) i
  analyse (lab ++ " (Q, Grevlex, Sing(gr))") "groebner" $ fmap (changeOrder Grevlex) i
  analyse (lab ++ " (Q, Grevlex, Sing(sba))") "sba" $ fmap (changeOrder Grevlex) i
  analyse (lab ++ " (F_65521, Lex, Sing(gr))") "groebner" $ fmap (mapCoeff ratToF . changeOrder Lex) i
  analyse (lab ++ " (F_65521, Lex, Sing(sba))") "sba" $ fmap (mapCoeff ratToF . changeOrder Lex) i
  analyse (lab ++ " (F_65521, Grevlex, Sing(gr))") "groebner" $ fmap (mapCoeff ratToF . changeOrder Grevlex) i
  analyse (lab ++ " (F_65521, Grevlex, Sing(sba))") "sba" $ fmap (mapCoeff ratToF . changeOrder Grevlex) i

main :: IO ()
main = do
  runTestCases "Cyclic-4" $ cyclic (sing :: Sing 4)
  runTestCases "Cyclic-5" $ cyclic (sing :: Sing 5)
  runTestCases "Cyclic-6" $ cyclic (sing :: Sing 6)
  runTestCases "Katsura-5" $ katsura (sing :: Sing 5)
  runTestCases "Katsura-6" $ katsura (sing :: Sing 6)
  runTestCases "Katsura-7" $ katsura (sing :: Sing 7)

ratToF :: Rational -> F 65521
ratToF = modRat'

ext: hs


text:
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications    #-}
-- |
-- Module      : Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Level2
-- Copyright   : [2017..2020] Trevor L. McDonell
-- License     : BSD3
--
-- Maintainer  : Trevor L. McDonell <trevor.mcdonell@gmail.com>
-- Stability   : experimental
-- Portability : non-portable (GHC extensions)
--

module Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Level2
  where

import Data.Complex
import Data.Array.Accelerate.Representation.Array
import Data.Array.Accelerate.Representation.Shape
import Data.Array.Accelerate.Sugar.Elt

import Data.Array.Accelerate.LLVM.PTX.Foreign
import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Base
import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Context
import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Level3
import Data.Array.Accelerate.Numeric.LinearAlgebra.Type

import Foreign.Marshal                                              ( with )
import qualified Foreign.CUDA.Ptr                                   as CUDA
import qualified Foreign.CUDA.BLAS                                  as BLAS

import Control.Monad.Reader

-- NOTE: cuBLAS requires matrices to be stored in column-major order
-- (Fortran-style), but Accelerate uses C-style arrays in row-major order.
--
-- If the operation is N or T, we can just swap the operation. For
-- conjugate-transpose (H) operations (on complex valued arguments), since there
-- is no conjugate-no-transpose operation, we implement that via 'gemm', which
-- I assume is more efficient than ?geam followed by ?gemv.
--
gemv :: NumericR s e
     -> Transpose
     -> ForeignAcc ((((((), Scalar e), Matrix e), Vector e)) -> Vector e)
gemv eR opA = ForeignAcc "ptx.gemv" (gemv' eR opA)

gemv' :: NumericR s e
      -> Transpose
      -> ((((), Scalar e), Matrix e), Vector e)
      -> Par PTX (Future (Vector e))
gemv' NumericRcomplex32 H = as_gemm NumericRcomplex32 H
gemv' NumericRcomplex64 H = as_gemm NumericRcomplex64 H
gemv' nR                t = as_gemv nR t

as_gemm
    :: NumericR s e
    -> Transpose
    -> ((((), Scalar e), Matrix e), Vector e)
    -> Par PTX (Future (Vector e))
as_gemm nR opA ((((), alpha), matA), Array sh adata) = do
  let matB = Array (sh,1) adata
  --
  future <- new
  result <- gemm' nR opA N ((((), alpha), matA), matB)
  fork $ do Array (sh',_) vecy <- get result
            put future (Array sh' vecy)
  return future

as_gemv
    :: NumericR s e
    -> Transpose
    -> ((((), Scalar e), Matrix e), Vector e)
    -> Par PTX (Future (Vector e))
as_gemv nR opA ((((), alpha), matA), vecx) = do
  let
      (((), rowsA), colsA) = shape matA

      sizeY   = case opA of
                  N -> rowsA
                  _ -> colsA

      opA'    = encodeTranspose
              $ case opA of
                  N -> T
                  _ -> N

      aR      = ArrayR dim1 eR
      eR      = case nR of
                  NumericRfloat32   -> eltR @Float
                  NumericRfloat64   -> eltR @Double
                  NumericRcomplex32 -> eltR @(Complex Float)
                  NumericRcomplex64 -> eltR @(Complex Double)
  --
  future  <- new
  stream  <- asks ptxStream
  vecy    <- allocateRemote aR ((), sizeY)
  alpha'  <- indexRemote eR alpha 0
  ()      <- liftPar $ do
    withArray nR matA stream   $ \ptr_A -> do
     withArray nR vecx stream  $ \ptr_x -> do
      withArray nR vecy stream $ \ptr_y -> do
       withBLAS                $ \hdl   -> do
         case nR of
           NumericRfloat32 -> liftIO $
            with alpha' $ \ptr_alpha ->
             with 0     $ \ptr_beta  ->
               BLAS.sgemv hdl opA' colsA rowsA ptr_alpha ptr_A colsA ptr_x 1 ptr_beta ptr_y 1

           NumericRfloat64 -> liftIO $
            with alpha' $ \ptr_alpha ->
             with 0     $ \ptr_beta  ->
               BLAS.dgemv hdl opA' colsA rowsA ptr_alpha ptr_A colsA ptr_x 1 ptr_beta ptr_y 1

           NumericRcomplex32 -> liftIO $
            withV2 nR alpha' $ \ptr_alpha ->
             with 0          $ \ptr_beta  ->
               BLAS.cgemv hdl opA' colsA rowsA ptr_alpha (CUDA.castDevPtr ptr_A) colsA (CUDA.castDevPtr ptr_x) 1 ptr_beta (CUDA.castDevPtr ptr_y)  1

           NumericRcomplex64 -> liftIO $
            withV2 nR alpha' $ \ptr_alpha ->
             with 0          $ \ptr_beta  ->
               BLAS.zgemv hdl opA' colsA rowsA ptr_alpha (CUDA.castDevPtr ptr_A) colsA (CUDA.castDevPtr ptr_x) 1 ptr_beta (CUDA.castDevPtr ptr_y)  1
  --
  put future vecy
  return future

max_stars_repo_name: statusfailed/accelerate-blas

max_issues_repo_name: statusfailed/accelerate-blas

max_forks_repo_name: statusfailed/accelerate-blas


text:
import Statistics.Distribution
import Statistics.Distribution.Normal
import Statistics.Matrix
import Statistics.MLE
import System.Random.MWC
import Control.Applicative
import Control.Monad
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed as U

main :: IO ()
main = do
	let n = 1000 :: Int
	let k = 3 :: Int
	let b0 = 0.0
	let b1 = U.fromList $ take k (enumFrom 1.0)
	let ratio b = G.map (\bk -> bk / (b G.! 0)) b

	gen <- create
	e <- fromList n 1 <$> replicateM n (genContinous (normalDistr 0 1) gen)
	x <- fromList n k <$> replicateM (n * k) (genContinous (normalDistr 0 1) gen)
	let y = G.zipWith3 (\x y z -> if x+y+z > 0 then 1.0 else 0.0 :: Double ) (G.replicate n b0) (x `multiplyV` b1) (toVector e)
	putStrLn $ "True values: " ++ show (ratio b1) ++ " Estimated values: " ++ show (ratio (probit x y))

max_stars_repo_name: StefanHubner/statistics

max_issues_repo_name: StefanHubner/statistics

max_forks_repo_name: StefanHubner/statistics


text:
module D03
    ( Claim (..), mkClaims, claimsIntersectionSize, isolatedClaim
    ) where

import Data.List
import qualified Data.HashMap.Strict as HM
import Data.Complex
import Text.Regex

data Claim = Claim { identifier :: Int,
                     topLeft :: Complex Int,
                     bottomRight :: Complex Int }
  deriving (Show, Eq)

type World = HM.HashMap (Complex Int) [Claim]

isolatedClaim :: String -> Int
isolatedClaim input = identifier $ head $ claims \\ badClaims
  where
    badClaims = concat $ HM.elems $ HM.filter (\cs -> (length cs) > 1) $ buildWorld claims
    claims = mkClaims input

claimsIntersectionSize :: String -> Int
claimsIntersectionSize = length . filter (> 1) . map length . HM.elems . buildWorld . mkClaims

buildWorld :: [Claim] -> World
buildWorld = foldl addClaimToWorld HM.empty

addClaimToWorld :: World -> Claim -> World
addClaimToWorld world claim = foldl addPoint world [x:+y | x<-[(realPart p1)..(realPart p2)], y<-[(imagPart p1)..(imagPart p2)]]
  where
    p1 = topLeft claim
    p2 = bottomRight claim
    addPoint w p = HM.insertWith (++) p [claim] w

mkClaims :: String -> [Claim]
mkClaims = map mkClaim . lines

mkClaim :: String -> Claim
mkClaim def = Claim (params !! 0) (left :+ top) ((left+width-1) :+ (top+height-1))
  where
    top = params !! 2
    left = params !! 1
    width = params !! 3
    height = params !! 4
    params = extract $ matchRegexAll (mkRegex "#([0-9]+) @ ([0-9]+),([0-9]+): ([0-9]+)x([0-9]+)") def
    extract (Just (_,_,_,subs)) = map read subs


text:
{-# LANGUAGE NoMonomorphismRestriction, StandaloneDeriving, UndecidableInstances #-}
module L where
import Prelude hiding ((.),(++),flip)
import Control.Applicative
import Control.Arrow
import Control.Arrow.Operations
import Control.Monad
import Control.Monad.Cont
import Control.Monad.Error
import Control.Monad.Fix
import Control.Monad.Identity
import Control.Monad.Instances
import Control.Monad.Logic
import Control.Monad.RWS
import Control.Monad.Reader
import Control.Monad.ST (ST, runST, fixST)
import Control.Monad.State
import Control.Monad.Writer
import Control.Parallel
import Control.Parallel.Strategies
import Data.Array
import Data.Bits
import Data.Bool
import Data.Char
import Data.Complex
import Data.Dynamic
import Data.Either
import Data.Eq
import Data.Fixed
import Data.Function hiding ((.),flip)
import Data.Generics hiding (GT)
import Data.Graph
import Data.Int
import Data.Ix
import Data.List hiding ((++),map)
import Data.Maybe
import Data.Monoid
import Data.Number.BigFloat
import Data.Number.CReal
import Data.Number.Dif
import Data.Number.Interval
import Data.Number.Symbolic
import Data.Ord
import Data.Ratio
import Data.STRef
import Data.Tree
import Data.Tuple
import Data.Typeable
import Data.Word
import Numeric
import ShowQ
import System.Random
import Text.PrettyPrint.HughesPJ hiding (empty)
import Text.Printf
import Text.Regex.Posix
import qualified Control.Arrow.Transformer as AT
import qualified Control.Arrow.Transformer.All as AT
import qualified Data.ByteString as BS
import qualified Data.ByteString.Char8 as BSC
import qualified Data.ByteString.Lazy as BSL
import qualified Data.ByteString.Lazy.Char8 as BSLC
import qualified Data.Foldable
import qualified Data.Generics
import qualified Data.IntMap as IM
import qualified Data.IntSet as IS
import qualified Data.Map as M
import qualified Data.Sequence as Seq
import qualified Data.Set as S
import qualified Data.Traversable
import qualified Data.MemoCombinators as Memo
import ShowIO

import SimpleReflect hiding (var)
import Math.OEIS
import Control.Monad.Random

import Data.VectorSpace hiding (Sum, getSum)
import Data.NumInstances
import Data.LinearMap

import qualified Test.IOSpec as IOSpec

import Data.Fixed

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

flip :: Functor f => f (a -> b) -> a -> f b
flip f x = fmap ($ x) f

(++) :: (Monoid m) => m -> m -> m
(++) = mappend
infixr 5 ++

asTypeIn :: a -> (a -> b) -> a
a `asTypeIn` f = a where _ = f a
infixl 0 `asTypeIn`

f `asAppliedTo` x = f `asTypeIn` \f -> f x
infixl 0 `asAppliedTo`

describeSequence = fmap description . lookupSequence

newtype Mu f = In (f (Mu f))
out (In x) = x
deriving instance Show (f (Mu f)) => Show (Mu f)

newtype Rec a = InR { outR :: Rec a -> a }

data Bin a = Tip | Branch a (Bin a) (Bin a)
  deriving (Eq, Ord, Show, Read)
-- Fun Stuff

cake = ["One 18.25 ounce package chocolate cake mix.",
        "One can prepared coconut pecan frosting.",
        "Three slash four cup vegetable oil.",
        "Four large eggs.",
        "One cup semi-sweet chocolate chips.",
        "Three slash four cups butter or margarine.",
        "One and two third cups granulated sugar.",
        "Two cups all-purpose flour.",
        "Don't forget garnishes such as:",
        "Fish shaped crackers.",
        "Fish shaped candies.",
        "Fish shaped solid waste.",
        "Fish shaped dirt.",
        "Fish shaped ethylbenzene.",
        "Pull and peel licorice.",
        "Fish shaped organic compounds and sediment shaped sediment.",
        "Candy coated peanut butter pieces. Shaped like fish.",
        "One cup lemon juice.",
        "Alpha resins.",
        "Unsaturated polyester resin.",
        "Fiberglass surface resins.",
        "And volatile malted milk impoundments.",
        "Nine large egg yolks.",
        "Twelve medium geosynthetic membranes.",
        "One cup granulated sugar.",
        "An entry called 'how to kill someone with your bare hands.'",
        "Two cups rhubarb, sliced.",
        "Two slash three cups granulated rhubarb.",
        "One tablespoon all-purpose rhubarb.",
        "One teaspoon grated orange rhubarb.",
        "Three tablespoons rhubarb, on fire.",
        "One large rhubarb.",
        "One cross borehole electro-magnetic imaging rhubarb.",
        "Two tablespoons rhubarb juice.",
        "Adjustable aluminum head positioner.",
        "Slaughter electric needle injector.",
        "Cordless electric needle injector.",
        "Injector needle driver.",
        "Injector needle gun.",
        "Cranial caps.",
        "And it contains proven preservatives, deep penetration agents, and gas and odor control chemicals.",
        "That will deodorize and preserve putrid tissue."]

{-# LINE 1 "<local>" #-}
(****) = (\f n v -> (iterate f v) !! n)
toChurch = ((appEndo . Data.Foldable.foldMap Endo) .) . replicate
fromChuch f = f (+1) 0
evalCont :: Cont o o -> o; evalCont = (`runCont` id)

max_stars_repo_name: jwiegley/lambdabot


text:
{-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  Statistics.Distribution.Beta
-- Copyright   :  (C) 2012 Edward Kmett,
-- License     :  BSD-style (see the file LICENSE)
--
-- Maintainer  :  Edward Kmett <ekmett@gmail.com>
-- Stability   :  provisional
-- Portability :  DeriveDataTypeable
--
----------------------------------------------------------------------------
module Statistics.Distribution.Beta
  ( BetaDistribution
    -- * Constructor
  , betaDistr
  , improperBetaDistr
    -- * Accessors
  , bdAlpha
  , bdBeta
  ) where

import Data.Aeson (FromJSON, ToJSON)
import Data.Binary (Binary)
import Data.Data (Data, Typeable)
import GHC.Generics (Generic)
import Numeric.SpecFunctions (
  incompleteBeta, invIncompleteBeta, logBeta, digamma)
import Numeric.MathFunctions.Constants (m_NaN)
import qualified Statistics.Distribution as D
import Data.Binary (put, get)
import Control.Applicative ((<$>), (<*>))

-- | The beta distribution
data BetaDistribution = BD
 { bdAlpha :: {-# UNPACK #-} !Double
   -- ^ Alpha shape parameter
 , bdBeta  :: {-# UNPACK #-} !Double
   -- ^ Beta shape parameter
 } deriving (Eq, Read, Show, Typeable, Data, Generic)

instance FromJSON BetaDistribution
instance ToJSON BetaDistribution

instance Binary BetaDistribution where
    put (BD x y) = put x >> put y
    get = BD <$> get <*> get

-- | Create beta distribution. Both shape parameters must be positive.
betaDistr :: Double             -- ^ Shape parameter alpha
          -> Double             -- ^ Shape parameter beta
          -> BetaDistribution
betaDistr a b
  | a > 0 && b > 0 = improperBetaDistr a b
  | otherwise      =
      error $  "Statistics.Distribution.Beta.betaDistr: "
            ++ "shape parameters must be positive. Got a = "
            ++ show a
            ++ " b = "
            ++ show b

-- | Create beta distribution. This construtor doesn't check parameters.
improperBetaDistr :: Double             -- ^ Shape parameter alpha
                  -> Double             -- ^ Shape parameter beta
                  -> BetaDistribution
improperBetaDistr = BD

instance D.Distribution BetaDistribution where
  cumulative (BD a b) x
    | x <= 0    = 0
    | x >= 1    = 1
    | otherwise = incompleteBeta a b x

instance D.Mean BetaDistribution where
  mean (BD a b) = a / (a + b)

instance D.MaybeMean BetaDistribution where
  maybeMean = Just . D.mean

instance D.Variance BetaDistribution where
  variance (BD a b) = a*b / (apb*apb*(apb+1))
    where apb = a + b

instance D.MaybeVariance BetaDistribution where
  maybeVariance = Just . D.variance

instance D.Entropy BetaDistribution where
  entropy (BD a b) =
    logBeta a b
    - (a-1) * digamma a
    - (b-1) * digamma b
    + (a+b-2) * digamma (a+b)

instance D.MaybeEntropy BetaDistribution where
  maybeEntropy = Just . D.entropy

instance D.ContDistr BetaDistribution where
  density (BD a b) x
   | a <= 0 || b <= 0 = m_NaN
   | x <= 0 = 0
   | x >= 1 = 0
   | otherwise = exp $ (a-1)*log x + (b-1)*log (1-x) - logBeta a b

  quantile (BD a b) p
    | p == 0         = 0
    | p == 1         = 1
    | p > 0 && p < 1 = invIncompleteBeta a b p
    | otherwise      =
        error $ "Statistics.Distribution.Gamma.quantile: p must be in [0,1] range. Got: "++show p

instance D.ContGen BetaDistribution where
  genContVar = D.genContinous


text:
{-# LANGUAGE DeriveGeneric             #-}
{-# LANGUAGE ExistentialQuantification #-}
{-# LANGUAGE FlexibleInstances         #-}
{-# LANGUAGE InstanceSigs              #-}
{-# LANGUAGE LambdaCase                #-}
{-# LANGUAGE MultiParamTypeClasses     #-}
{-# LANGUAGE PolyKinds                 #-}
{-# LANGUAGE StandaloneDeriving        #-}
{-# LANGUAGE TypeFamilies              #-}
{-# LANGUAGE TypeInType                #-}

module Learn.Neural.Layer.FullyConnected (
    FullyConnected
  ) where

import           Data.Kind
import           Data.Singletons.Prelude
import           Data.Singletons.TypeLits
import           GHC.Generics                   (Generic)
import           GHC.Generics.Numeric
import           Learn.Neural.Layer
import           Numeric.BLAS
import           Numeric.Backprop
import           Statistics.Distribution
import           Statistics.Distribution.Normal
import qualified Generics.SOP                   as SOP

data FullyConnected :: Type

instance (Num (b '[o,i]), Num (b '[o])) => Num (CParam FullyConnected b '[i] '[o]) where
    FCP w1 b1 + FCP w2 b2 = FCP (w1 + w2) (b1 + b2)
    FCP w1 b1 - FCP w2 b2 = FCP (w1 - w2) (b1 - b2)
    FCP w1 b1 * FCP w2 b2 = FCP (w1 * w2) (b1 * b2)
    negate (FCP w b) = FCP (negate w) (negate b)
    signum (FCP w b) = FCP (signum w) (signum b)
    abs    (FCP w b) = FCP (abs    w) (abs    b)
    fromInteger x = FCP (fromInteger x) (fromInteger x)

instance (Fractional (b '[o,i]), Fractional (b '[o])) => Fractional (CParam FullyConnected b '[i] '[o]) where
    FCP w1 b1 / FCP w2 b2 = FCP (w1 / w2) (b1 / b2)
    recip (FCP w b)       = FCP (recip w) (recip b)
    fromRational x        = FCP (fromRational x) (fromRational x)

instance (Floating (b '[o,i]), Floating (b '[o])) => Floating (CParam FullyConnected b '[i] '[o]) where
    sqrt (FCP w b)       = FCP (sqrt w) (sqrt b)

instance Num (CState FullyConnected b '[i] '[o]) where
    _ + _         = FCS
    _ * _         = FCS
    _ - _         = FCS
    negate _      = FCS
    abs    _      = FCS
    signum _      = FCS
    fromInteger _ = FCS

instance Fractional (CState FullyConnected b '[i] '[o]) where
    _ / _          = FCS
    recip _        = FCS
    fromRational _ = FCS

instance Floating (CState FullyConnected b '[i] '[o]) where
    sqrt _ = FCS

deriving instance Generic (CParam FullyConnected b '[i] '[o])
instance SOP.Generic (CParam FullyConnected b '[i] '[o])

instance (BLAS b, KnownNat i, KnownNat o, Floating (b '[o,i]), Floating (b '[o]))
        => Component FullyConnected b '[i] '[o] where
    data CParam  FullyConnected b '[i] '[o] =
            FCP { _fcWeights :: !(b '[o,i])
                , _fcBiases  :: !(b '[o])
                }
    data CState  FullyConnected b '[i] '[o] = FCS
    type CConstr FullyConnected b '[i] '[o] = Num (b '[o,i])
    data CConf   FullyConnected b '[i] '[o] = forall d. ContGen d => FCC d

    componentOp = componentOpDefault

    initParam = \case
      i `SCons` SNil -> \case
        so@(o `SCons` SNil) -> \(FCC d) g -> do
          w <- genA (o `SCons` (i `SCons` SNil)) $ \_ ->
            realToFrac <$> genContVar d g
          b <- genA so $ \_ ->
            realToFrac <$> genContVar d g
          return $ FCP w b
        _ -> error "inaccessible."
      _ -> error "inaccessible."

    initState _ _ _ _ = return FCS

    defConf = FCC (normalDistr 0 0.01)

instance (BLAS b, KnownNat i, KnownNat o, Floating (b '[o,i]), Floating (b '[o]))
        => ComponentFF FullyConnected b '[i] '[o] where
    componentOpFF = bpOp . withInps $ \(x :< p :< Ø) -> do
        w :< b :< Ø <- gTuple #<~ p
        y <- matVecOp ~$ (w :< x :< Ø)
        z <- (+.)     ~$ (y :< b :< Ø)
        return . only $ z

instance (BLAS b, KnownNat i, KnownNat o, Floating (b '[o,i]), Floating (b '[o]))
        => ComponentLayer r FullyConnected b '[i] '[o] where
    componentRunMode = RMIsFF

max_stars_repo_name: mstksg/backprop-learn

max_issues_repo_name: mstksg/backprop-learn

max_forks_repo_name: mstksg/backprop-learn


text:
--------------------------------------------------------------------------------
--
--  Copyright (c) 2010 - 2013 Tad Doxsee
--  MIT License
--
--  Author: Tad Doxsee
--
--------------------------------------------------------------------------------

{-# LANGUAGE CPP                      #-}
{-# LANGUAGE ForeignFunctionInterface #-}
{-# LANGUAGE EmptyDataDecls           #-}

module Numeric.LinearAlgebra.CHOLMOD.CholmodXFace where

-- base
import Foreign.Ptr
import Foreign.ForeignPtr
import Foreign.C.Types
import Foreign.Marshal.Array (withArray)

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

-- hcholmod
import Numeric.LinearAlgebra.CHOLMOD.CholmodXFaceLow

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

data Matrix a = Matrix {fPtr :: ForeignPtr a}
-- a can be types Dense, Sparse, or Triplet

class MatrixC a where
   getNRow  :: Matrix a -> IO CSize
   getNCol  :: Matrix a -> IO CSize
   getNZMax :: Matrix a -> IO CSize
   getX     :: Matrix a -> IO (V.MVector s CDouble)
   mtxFree  :: ForeignPtr Common -> Matrix a -> IO CInt
--   finalizer :: FinalizerEnvPtr Common a

instance MatrixC Triplet where
    getNRow  m = withForeignPtr (fPtr m) triplet_get_nrow
    getNCol  m = withForeignPtr (fPtr m) triplet_get_ncol
    getNZMax m = withForeignPtr (fPtr m) triplet_get_nzmax
    getX     m = tripletGetX m

    mtxFree cfp  m = withForeignPtr (fPtr m) $ \mp -> do
                       withForeignPtr cfp $ \cp -> do
                         triplet_free mp cp

--    finalizer  = triplet_free_ptr
--    free c   m = to fix: complete free functions

instance MatrixC Dense where
    getNRow  m = withForeignPtr (fPtr m) dense_get_nrow
    getNCol  m = withForeignPtr (fPtr m) dense_get_ncol
    getNZMax m = withForeignPtr (fPtr m) dense_get_nzmax
    getX     m = denseGetX m

    mtxFree cfp  m = withForeignPtr (fPtr m) $ \mp -> do
                       withForeignPtr cfp $ \cp -> do
                         dense_free mp cp
--    finalizer  = dense_free_ptr

instance MatrixC Sparse where
    getNRow  m = withForeignPtr (fPtr m) sparse_get_nrow
    getNCol  m = withForeignPtr (fPtr m) sparse_get_ncol
    getNZMax m = withForeignPtr (fPtr m) sparse_get_nzmax
    getX     m = sparseGetX m

    mtxFree cfp  m = withForeignPtr (fPtr m) $ \mp -> do
                       withForeignPtr cfp $ \cp -> do
                         sparse_free mp cp

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

allocTriplet :: CSize -> CSize -> CSize ->
                SType -> XType -> ForeignPtr Common ->
                         IO (Matrix Triplet)

allocTriplet nRow nCol nZMax sType xType cfp =
    withForeignPtr cfp $ \cp -> do
      (triplet_allocate nRow nCol nZMax sType xType cp)
      >>= newForeignPtr_
      >>= (return . Matrix)

allocCommon :: IO (ForeignPtr Common)
allocCommon = newForeignPtr common_finish_and_free_ptr =<< common_allocate

startC :: ForeignPtr Common -> IO ()
startC cfp = withForeignPtr cfp start

tripletGetRowIndices :: Matrix Triplet -> IO (V.IOVector CInt)
tripletGetRowIndices m = do
  ip <- withForeignPtr (fPtr m) (triplet_get_row_indices)     :: IO (Ptr CInt)
  ipp <- newForeignPtr_ ip                              :: IO (ForeignPtr CInt)
  nZMax <- getNZMax m                                         :: IO CSize
  return $ V.unsafeFromForeignPtr ipp 0 (fromIntegral nZMax)

tripletGetColIndices :: Matrix Triplet -> IO (V.IOVector CInt)
tripletGetColIndices m = do
  jp <- withForeignPtr (fPtr m) (triplet_get_col_indices) :: IO (Ptr CInt)
  jpp <- newForeignPtr_ jp                               :: IO (ForeignPtr CInt)
  nZMax <- getNZMax m                                    :: IO CSize
  return $ V.unsafeFromForeignPtr jpp 0 (fromIntegral nZMax)

tripletGetX :: Matrix Triplet -> IO (V.MVector s CDouble)
tripletGetX m = do
  xp <- withForeignPtr (fPtr m) triplet_get_x  :: IO (Ptr CDouble)
  xpp <- newForeignPtr_ xp                     :: IO (ForeignPtr CDouble)
  nZMax <- getNZMax m                          :: IO CSize
  return $ V.unsafeFromForeignPtr xpp 0 (fromIntegral nZMax)

tripletSetNNZ :: Matrix Triplet -> CSize -> IO ()
tripletSetNNZ m nnz = withForeignPtr (fPtr m) $ \mp -> (triplet_set_nnz mp nnz)

tripletGetNNZ :: Matrix Triplet -> IO CSize
tripletGetNNZ m = withForeignPtr (fPtr m) triplet_get_nnz

tripletToSparse :: Matrix Triplet -> ForeignPtr Common -> IO (Matrix Sparse)
tripletToSparse t cfp = do
    nnz <- tripletGetNNZ t
    withForeignPtr cfp $ \cp -> do
      withForeignPtr (fPtr t) $ \tp -> do
         sp <- triplet_to_sparse tp nnz cp :: IO (Ptr Sparse)
         sfp <- newForeignPtr_ sp :: IO (ForeignPtr Sparse)
         return (Matrix sfp)

zeros :: CSize -> CSize -> XType -> ForeignPtr Common -> IO (Matrix Dense)
zeros nRow nCol xt cfp =
  withForeignPtr cfp $ \cp -> do
    op <- zerosL nRow nCol xt cp :: IO (Ptr Dense)
    ofp<- newForeignPtr_ op :: IO (ForeignPtr Dense)
    return (Matrix ofp)

ones :: CSize -> CSize -> XType -> ForeignPtr Common -> IO (Matrix Dense)
ones nRow nCol xt cfp =
  withForeignPtr cfp $ \cp -> do
    op <- onesL nRow nCol xt cp :: IO (Ptr Dense)
    ofp<- newForeignPtr_ op :: IO (ForeignPtr Dense)
    return (Matrix ofp)

eye :: CSize -> CSize -> XType -> ForeignPtr Common -> IO (Matrix Dense)
eye nRow nCol xt cfp =
  withForeignPtr cfp $ \cp -> do
    op <- eyeL nRow nCol xt cp :: IO (Ptr Dense)
    ofp<- newForeignPtr_ op :: IO (ForeignPtr Dense)
    return (Matrix ofp)

analyze :: Matrix Sparse -> ForeignPtr Common -> IO (ForeignPtr Factor)
analyze m cfp = withForeignPtr (fPtr m) $ \mp -> do
                   withForeignPtr cfp $ \cp -> do
                     lp <- analyzeL mp cp :: IO (Ptr Factor)
                     newForeignPtrEnv factor_free_ptr cp lp

factorize :: Matrix Sparse -> ForeignPtr Factor -> ForeignPtr Common -> IO ()
factorize m ffp cfp =
    withForeignPtr (fPtr m) $ \mp -> do
      withForeignPtr ffp $ \fp -> do
        withForeignPtr cfp $ \cp ->  (factorizeL mp fp cp)

solve ::  SystemType
      -> ForeignPtr Factor
      -> Matrix Dense
      -> ForeignPtr Common
      -> IO (Matrix Dense)

solve st ffp m cfp =
    withForeignPtr ffp $ \fp -> do
      withForeignPtr (fPtr m) $ \mp -> do
        withForeignPtr cfp $ \cp -> do
          xp <- solveL st fp mp cp
          xfp <- newForeignPtr_ xp
          return (Matrix xfp)

denseGetX :: Matrix Dense -> IO (V.MVector s CDouble)
denseGetX m = do

  nZMax <- getNZMax m                       :: IO CSize
  xp <- withForeignPtr (fPtr m) dense_get_x :: IO (Ptr CDouble)
  xpp <- newForeignPtr_ xp                  :: IO (ForeignPtr CDouble)
  return $ V.unsafeFromForeignPtr xpp 0 (fromIntegral nZMax)

denseCopy :: Matrix Dense -> ForeignPtr Common -> IO (Matrix Dense)
denseCopy m cfp =
    withForeignPtr (fPtr m) $ \mp -> do
      withForeignPtr cfp $ \cp -> do
        m2p <- dense_copy mp cp
        m2fp <- newForeignPtr_ m2p
        return (Matrix m2fp)

sdMult :: Matrix Sparse  -- ^ A, sparse matrix to multiply
       -> DoTranspose    -- ^
       -> [CDouble]      -- ^ alpha[2], scale factor for A
       -> [CDouble]      -- ^ beta[2],  scale factor for Y
       -> Matrix Dense   -- ^ X, dense matrix to multiply
       -> Matrix Dense   -- ^ Y, resulting dense matrix
       -> ForeignPtr Common
       -> IO CInt

sdMult s doT alphaL betaL x y cfp =
    withForeignPtr (fPtr s) $ \sp -> do
      withArray alphaL $ \alpha -> do
        withArray betaL $ \beta -> do
          withForeignPtr (fPtr x) $ \xp -> do
            withForeignPtr (fPtr y) $ \yp -> do
              withForeignPtr cfp $ \cp -> do
                sdmult sp doT alpha beta xp yp cp

denseNorm :: Matrix Dense -> NormType -> ForeignPtr Common -> IO CDouble

denseNorm m nt cfp =
    withForeignPtr (fPtr m) $ \mp -> do
      withForeignPtr cfp $ \cp -> do
        dense_norm mp nt cp

sparseToDense :: Matrix Sparse -> ForeignPtr Common -> IO (Matrix Dense)
sparseToDense m cfp =
    withForeignPtr (fPtr m) $ \mp -> do
      withForeignPtr cfp $ \cp -> do
        m2p <- sparse_to_dense mp cp
        m2fp <- newForeignPtr_ m2p
        return (Matrix m2fp)

sparseGetX :: Matrix Sparse -> IO (V.MVector s CDouble)
sparseGetX m = do
  xp <- withForeignPtr (fPtr m) sparse_get_x   :: IO (Ptr CDouble)
  xpp <- newForeignPtr_ xp                     :: IO (ForeignPtr CDouble)
  nZMax <- getNZMax m                          :: IO CSize
  return $ V.unsafeFromForeignPtr xpp 0 (fromIntegral nZMax)

max_stars_repo_name: tdox/hcholmod

max_issues_repo_name: tdox/hcholmod

max_forks_repo_name: tdox/hcholmod


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE ExistentialQuantification #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE PolyKinds #-}

module CSimpleSpec where

import Commons
import Control.Applicative (liftA2)
import Control.Concurrent
import Control.Monad (replicateM_, unless, when)
import Data.Array
import Data.Complex (Complex (..))
import qualified Data.IntMap.Strict as IM
import Data.List (intercalate, sort, tails)
import Data.List.Split (splitOn)
import qualified Data.Map.Strict as Map
import Data.Maybe (fromJust, mapMaybe)
import qualified Data.String.Interpolate as I
import qualified Data.Text as T
import qualified Data.Text.IO as TIO
import Debug.Trace (traceShowId)
import GHC.IO.Exception (ExitCode (..))
import HashedExpression.Codegen
import HashedExpression.Codegen.CSimple
import HashedExpression.Embed.FFTW
import HashedExpression.Internal.Expression hiding (InnerProductSpaceOp (..))
import HashedExpression.Internal.Inner
import HashedExpression.Internal.Node
import HashedExpression.Internal.Normalize (normalize)
import HashedExpression.Internal.Utils
import HashedExpression.Interp
import HashedExpression.Prettify (showExp, showExpDebug)
import HashedExpression.Value
import System.FilePath
import System.Process (readProcess, readProcessWithExitCode)
import Test.Hspec
import Test.QuickCheck
import Var
import Prelude hiding ((!!))

-- | Compute the inner offset:
-- e.g localOffset [3, 4, 5] [2, 0, 1] = 2 * (4 * 5) + 0 * (5) + 1
localOffset :: [Int] -> [Int] -> Int
localOffset shape indices
  | length shape == length indices =
    sum . zipWith (*) indices . map product . tail . tails $ shape
  | otherwise = error $ "shape and indices are not compatible" ++ show shape ++ show indices

-- | Generate a fully working C program that compute the expression and
-- print out the result, mostly used for testing
singleExpressionCProgram ::
  (DimensionType d, NumType et) => ValMaps -> Expression d et -> Code
singleExpressionCProgram valMaps expr =
  [ "#include <math.h>", --
    "#include <stdio.h>",
    "#include <stdlib.h>",
    if containsFTNode $ exMap expr
      then T.pack $ fftUtils
      else "",
    "int main()" --
  ]
    ++ scoped (initMemory ++ assignVals ++ codes ++ printValue ++ releaseMemory)
  where
    (mp, n) = unwrap expr
    bound = product (retrieveShape n mp)
    et = retrieveElementType n mp
    codeGen = initCodegen (CodegenInit mp []) CSimpleConfig
    [i, j, k, nooffset] = ["i", "j", "k", "0"]
    initMemory = [[I.i|double *ptr = malloc(sizeof(double) * #{cMemSize codeGen});|]]
    -- assign value to variables
    assignVals = assigningValues codeGen valMaps
    -- codes to compute
    codes = evaluating codeGen [n]
    -- print the value of expression
    printValue
      | et == R = for i bound [[I.i|printf("%f ", #{(!!) codeGen n i});|]]
      | et == C =
        for i bound [[I.i|printf("%f ", #{(!!) codeGen n i});|]]
          ++ [[I.i|printf("\\n");|]]
          ++ for i bound [[I.i|printf("%f ", #{imAt codeGen n i});|]]
    releaseMemory = ["free(ptr);"]
    -------------------------------------------------------------------------------
    assigningValues :: CSimpleCodegen -> ValMaps -> Code
    assigningValues CSimpleCodegen {..} valMaps = concatMap codesForVar vars
      where
        [i, j, k, nooffset] = ["i", "j", "k", "0"]
        vars :: [(Int, String)]
        vars =
          let toVar nId
                | Var varName <- retrieveNode nId mp = Just (nId, varName)
                | otherwise = Nothing
           in mapMaybe toVar . IM.keys $ mp
        codesForVar :: (Int, String) -> Code
        codesForVar (n, varName) =
          case Map.lookup varName valMaps of
            Just (VScalar val) -> [[I.i|#{n !! nooffset} = #{val};|]]
            Just (V1D array1d) ->
              let assignIndex id = [[I.i|#{n !! (showT id)} = #{array1d ! id};|]]
               in concatMap assignIndex $ indices array1d
            Just (V2D array2d) ->
              let shape = retrieveShape n mp
                  assignIndex (id1, id2) =
                    [[I.i|#{n !! showT (localOffset shape [id1, id2])} = #{array2d ! (id1, id2)};|]]
               in concatMap assignIndex $ indices array2d
            Just (V3D array3d) ->
              let shape = retrieveShape n mp
                  assignIndex (id1, id2, id3) =
                    [[I.i|#{n !! showT (localOffset shape [id1, id2, id3])} = #{array3d ! (id1, id2, id2)};|]]
               in concatMap assignIndex $ indices array3d
            _ -> []

hasFFTW :: IO Bool
hasFFTW = do
  fileName <- generate $ vectorOf 10 $ elements ['A' .. 'Z']
  let cFilePath = "C" </> fileName <.> "c"
      executablePath = "C" </> fileName
      codes = ["#include <fftw3.h>", "int main() {}"]
      program = T.intercalate "\n" . map T.pack $ codes
  TIO.writeFile cFilePath program
  (exitCode, output, _) <- readProcessWithExitCode "gcc" [cFilePath, "-o", executablePath, "-lm", "-lfftw3"] ""
  readProcessWithExitCode "rm" [cFilePath] ""
  readProcessWithExitCode "rm" [executablePath] ""
  return $ exitCode == ExitSuccess

-- |
evaluateCodeC :: (DimensionType d, NumType et) => Bool -> Expression d et -> ValMaps -> IO (ExitCode, String)
evaluateCodeC withFT exp valMaps = do
  readProcessWithExitCode "mkdir" ["C"] ""
  fileName <- generate $ vectorOf 10 $ elements ['A' .. 'Z']
  let cFilePath = "C" </> fileName <.> "c"
      executablePath = "C" </> fileName
      program = singleExpressionCProgram valMaps exp
      libs
        | withFT = ["-lm", "-lfftw3"]
        | otherwise = ["-lm"]
  TIO.writeFile cFilePath (T.intercalate "\n" program)
  readProcess "gcc" ([cFilePath, "-o", executablePath] ++ libs) ""
  (exitCode, output, _) <- readProcessWithExitCode executablePath [] ""
  readProcess "rm" [cFilePath] ""
  readProcess "rm" [executablePath] ""
  return (exitCode, output)

-- | Parse output of the C program
readR :: String -> [Double]
readR = map read . filter (not . null) . splitOn " "

-- |
readC :: String -> ([Double], [Double])
readC str = (readR rePart, readR imPart)
  where
    [rePart, imPart] = splitOn "\n" str

-- |
prop_CEqualInterpScalarR :: SuiteScalarR -> Expectation
prop_CEqualInterpScalarR (Suite exp valMaps) =
  if containsFTNode $ exMap exp
    then do
      hasFTLib <- hasFFTW
      when hasFTLib $ proceed True
    else proceed False
  where
    proceed withFT = do
      let normalizedExp = normalize exp
      (exitCode, outputSimple) <- evaluateCodeC withFT normalizedExp valMaps
      let resultNormalize = read . head . splitOn " " $ outputSimple
      let resultInterpNormalize = eval valMaps normalizedExp
      resultNormalize `shouldApprox` resultInterpNormalize

-- |
prop_CEqualInterpScalarC :: SuiteScalarC -> Expectation
prop_CEqualInterpScalarC (Suite exp valMaps) =
  if containsFTNode $ exMap exp
    then do
      hasFTLib <- hasFFTW
      when hasFTLib $ proceed True
    else proceed False
  where
    proceed withFT = do
      let normalizedExp = normalize exp
      (exitCode, outputCodeC) <- evaluateCodeC withFT normalizedExp valMaps
      let ([im], [re]) = readC outputCodeC
      let resultNormalize = im :+ re
      let resultInterpNormalize = eval valMaps normalizedExp
      resultNormalize `shouldApprox` resultInterpNormalize

-- |
prop_CEqualInterpOneR :: SuiteOneR -> Expectation
prop_CEqualInterpOneR (Suite exp valMaps) =
  if containsFTNode $ exMap exp
    then do
      hasFTLib <- hasFFTW
      when hasFTLib $ proceed True
    else proceed False
  where
    proceed withFT = do
      let normalizedExp = normalize exp
      (exitCode, outputSimple) <- evaluateCodeC withFT normalizedExp valMaps
      let resultNormalize = listArray (0, defaultDim1D - 1) $ readR outputSimple
      let resultInterpNormalize = eval valMaps normalizedExp
      resultNormalize `shouldApprox` resultInterpNormalize

-- |
prop_CEqualInterpOneC :: SuiteOneC -> Expectation
prop_CEqualInterpOneC (Suite exp valMaps) =
  if containsFTNode $ exMap exp
    then do
      hasFTLib <- hasFFTW
      when hasFTLib $ proceed True
    else proceed False
  where
    proceed withFT = do
      let normalizedExp = normalize exp
      (exitCode, outputCodeC) <- evaluateCodeC withFT normalizedExp valMaps
      let (re, im) = readC outputCodeC
          resultNormalize = listArray (0, defaultDim1D - 1) $ zipWith (:+) re im
          resultInterpNormalize = eval valMaps normalizedExp
      resultNormalize `shouldApprox` resultInterpNormalize

-- |
prop_CEqualInterpTwoR :: SuiteTwoR -> Expectation
prop_CEqualInterpTwoR (Suite exp valMaps) =
  if containsFTNode $ exMap exp
    then do
      hasFTLib <- hasFFTW
      when hasFTLib $ proceed True
    else proceed False
  where
    proceed withFT = do
      let normalizedExp = normalize exp
      (exitCode, outputSimple) <- evaluateCodeC withFT normalizedExp valMaps
      let resultNormalize = listArray ((0, 0), (default1stDim2D - 1, default2ndDim2D - 1)) $ readR outputSimple
      let resultInterpNormalize = eval valMaps normalizedExp
      resultNormalize `shouldApprox` resultInterpNormalize

-- |
prop_CEqualInterpTwoC :: SuiteTwoC -> Expectation
prop_CEqualInterpTwoC (Suite exp valMaps) =
  if containsFTNode $ exMap exp
    then do
      hasFTLib <- hasFFTW
      when hasFTLib $ proceed True
    else proceed False
  where
    proceed withFT = do
      let normalizedExp = normalize exp
      (exitCode, outputCodeC) <- evaluateCodeC withFT (normalize exp) valMaps
      let (re, im) = readC outputCodeC
      let resultNormalize = listArray ((0, 0), (default1stDim2D - 1, default2ndDim2D - 1)) $ zipWith (:+) re im
      let resultInterpNormalize = eval valMaps normalizedExp
      resultNormalize `shouldApprox` resultInterpNormalize

-- | Spec
spec :: Spec
spec =
  describe "C simple spec" $ do
    specify
      "Evaluate hash interp should equal to C code evaluation (Expression Scalar R)"
      $ property prop_CEqualInterpScalarR
    specify
      "Evaluate hash interp should equal to C code evaluation (Expression Scalar C)"
      $ property prop_CEqualInterpScalarC
    specify
      "Evaluate hash interp should equal to C code evaluation (Expression One R)"
      $ property prop_CEqualInterpOneR
    specify
      "Evaluate hash interp should equal to C code evaluation (Expression One C)"
      $ property prop_CEqualInterpOneC
    specify
      "Evaluate hash interp should equal to C code evaluation (Expression Two R)"
      $ property prop_CEqualInterpTwoR
    specify
      "Evaluate hash interp should equal to C code evaluation (Expression Two C)"
      $ property prop_CEqualInterpTwoC

max_stars_repo_name: Turboscient/HashedExpression

max_issues_repo_name: Turboscient/HashedExpression

max_forks_repo_name: Turboscient/HashedExpression


text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE NoMonomorphismRestriction #-}

{-# OPTIONS -Wall #-}
-- | some functions for working with uniformly-sampled data
module Statistics.Iteratee.Uniform (
  someRollingFunction
, movingAverage
) where

import Statistics.Iteratee.Compat
import Statistics.Iteratee.Sample

import Control.Monad.Identity
import Data.Iteratee as I
import Data.ListLike (ListLike)

#if MIN_VERSION_iteratee(0,9,0)
#else
import qualified Data.ListLike as LL
#endif

#if MIN_VERSION_iteratee(0,9,0)
roll'
    :: (Monad m, ListLike s el)
    => Int  -- ^ length of chunk (t)
    -> Int  -- ^ amount to consume (d)
    -> Iteratee s m [s]
roll' = roll
#else
roll'
    :: (Monad m, Nullable s, ListLike s el)
    => Int  -- ^ length of chunk (t)
    -> Int  -- ^ amount to consume (d)
    -> Iteratee s m [s]
roll' t d
  | t > d  = liftI (go LL.empty)
  | otherwise = error "Iteratee.roll: (t <= d).  Reverse the args?"
    where
        go prev (Chunk vec) =
                let withPrev = prev `LL.append` vec
                in if LL.length withPrev > t
                    then idone [LL.take t withPrev] (Chunk $ LL.drop d withPrev)
                    else liftI (go withPrev)
        go prev e = idone [prev] e
#endif

someRollingFunction
    :: (Monad m, ListLikey s el)
    => Int
    -> (s -> summary)
    -> Enumeratee s [summary] m a
someRollingFunction count mkSummary =
    convStream (roll' count 1)
    ><> mapStream mkSummary
{-# INLINABLE someRollingFunction #-}

movingAverage
    :: (Fractional el, Monad m, ListLikey s el)
    => Int
    -> Enumeratee s [el] m a
movingAverage n = someRollingFunction n chunkMean
  where
    chunkMean = runIdentity . (run <=< flip enumPure1Chunk mean)
{-# INLINABLE movingAverage #-}

max_stars_repo_name: JohnLato/iter-stats

max_issues_repo_name: JohnLato/iter-stats

max_forks_repo_name: JohnLato/iter-stats


text:
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FlexibleContexts #-}

module Trade.Statistics.YieldStatistics where

import Data.Time.Clock (NominalDiffTime)

import qualified Data.List as List
import qualified Data.Vector as Vec

import Data.Ord (comparing)

import qualified Statistics.Sample as Sample

import Trade.Type.DeltaSignal (DeltaSignal(..))
import qualified Trade.Type.DeltaSignal.Algorithm as DSA
import Trade.Type.Yield (LogYield(..), logYield2yield)

import qualified Trade.Report.Table as Table

import Trade.Report.ToReport (ToReport, toReport)

import Trade.Statistics.Statistics (Statistics(..), formatStat) -- , formatYield, formatStat)
import Trade.Statistics.Algorithm (mean, stdDev, skewness, kurtosis)

import Trade.Report.Pretty (Pretty, pretty)

data YieldStatistics ohlc = YieldStatistics {
  count :: !Int
  , maximumYield :: LogYield ohlc
  , minimumYield :: LogYield ohlc
  , maximumDuration :: LogYield ohlc
  , minimumDuration :: LogYield ohlc
  , meanYield :: Statistics NominalDiffTime Double
  , stdDevYield :: Statistics NominalDiffTime Double
  , skewnessYield :: Statistics Double Double
  , kurtosisYield :: Statistics Double Double
  }

yieldStatistics :: [LogYield ohlc] -> Maybe (YieldStatistics ohlc)
yieldStatistics [] = Nothing
yieldStatistics ys = Just $
  let (dts, zs) = unzip (map (\(LogYield dt y) -> (dt, y)) ys)
      dtsVec = Vec.fromList dts
      zsVec = Vec.fromList zs
  in YieldStatistics {
    count = length ys
    , maximumYield = List.maximumBy (comparing logYield) ys
    , minimumYield = List.minimumBy (comparing logYield) ys
    , maximumDuration = List.maximumBy (comparing logDuration) ys
    , minimumDuration = List.minimumBy (comparing logDuration) ys
    , meanYield = Statistics (mean dtsVec) (exp (mean zsVec))
    , stdDevYield = Statistics (stdDev dtsVec) (exp (stdDev zsVec))
    , skewnessYield = Statistics (skewness dtsVec) (Sample.skewness zsVec)
    , kurtosisYield = Statistics (kurtosis dtsVec) (Sample.kurtosis zsVec)
    }

yieldStatistics2table ::
  (Pretty ohlc) =>
  YieldStatistics ohlc -> [[String]]
yieldStatistics2table ys =
  [ "No. of trades" : [pretty (count ys)]
  , []
  , [ "", "Yield", "Duration" ]
  , []
  , "Maximum yield trade" : [pretty (logYield2yield (maximumYield ys))]
  , "Minimum yield trade" : [pretty (logYield2yield (minimumYield ys))]
  , []
  , "Maximum duration trade" : [pretty (logYield2yield (maximumDuration ys))]
  , "Minimum duration trade" : [pretty (logYield2yield (minimumDuration ys))]
  , []
  , "Mean" : formatStat (meanYield ys)
  , "Standard dev." : formatStat (stdDevYield ys)
  , "Skewness (log yield)" : formatStat (skewnessYield ys)
  , "Kurtosis (log yield)" : formatStat (kurtosisYield ys)
  ]

instance (Pretty ohlc) => ToReport (YieldStatistics ohlc) where
  toReport = Table.table . yieldStatistics2table

toYieldStatistics ::
  (Functor f) =>
  f [DeltaSignal ohlc] -> f (Maybe (YieldStatistics ohlc))
toYieldStatistics = fmap (yieldStatistics . map DSA.yield)

max_stars_repo_name: fphh/trade

max_issues_repo_name: fphh/trade

max_forks_repo_name: fphh/trade


text:
{-# LANGUAGE FlexibleContexts #-}
module PlotGreensFunctionR2Z2T0S0 where

import           Control.Monad.Parallel  as MP
import           Data.Array.Repa         as R
import           Data.Binary             (decodeFile)
import           Data.Complex
import           Data.List               as L
import           FokkerPlanck.MonteCarlo
import           FokkerPlanck.Pinwheel
import           Image.IO
import           System.Directory
import           System.Environment
import           System.FilePath
import           Types
import           Utils.Array
import Data.Vector.Unboxed as VU

main = do
  args@(numPointStr:numOrientationStr:thetaSigmaStr:scaleSigmaStr:maxScaleStr:taoStr:lenStr:initStr:numTrailStr:maxTrailStr:theta0FreqsStr:thetaFreqsStr:scale0FreqsStr:scaleFreqsStr:histFileName:alphaStr:numThreadStr:_) <-
    getArgs
  print args
  let numPoint = read numPointStr :: Int
      numOrientation = read numOrientationStr :: Int
      thetaSigma = read thetaSigmaStr :: Double
      scaleSigma = read scaleSigmaStr :: Double
      maxScale = read maxScaleStr :: Double
      tao = read taoStr :: Double
      len = read lenStr :: Int
      init = read initStr :: (Double, Double, Double, Double, Double, Double)
      numTrail = read numTrailStr :: Int
      maxTrail = read maxTrailStr :: Int
      theta0Freq = read theta0FreqsStr :: Double
      thetaFreq = read thetaFreqsStr :: Double
      scale0Freq = read scale0FreqsStr :: Double
      scaleFreq = read scaleFreqsStr :: Double
      theta0Freqs = [-theta0Freq .. -theta0Freq]
      thetaFreqs = [thetaFreq .. thetaFreq]
      scale0Freqs = [(scale0Freq - 1) .. scale0Freq]
      scaleFreqs = [-scaleFreq .. scaleFreq]
      alpha = read alphaStr :: Double
      numThread = read numThreadStr :: Int
      folderPath = "output/app/PlotGreensFunctionR2Z2T0S0"
      histPath = folderPath </> histFileName
      radialHistPath = folderPath </> ("Radial_" L.++ histFileName)
  createDirectoryIfMissing True (folderPath </> "GreensR2Z2T0S0")
  flag <- doesFileExist histPath
  hist <-
    if flag
      then decodeFile histPath
      else return $
           emptyHistogram
             [ numPoint
             , numPoint
             , L.length scale0Freqs
             , L.length theta0Freqs
             , L.length scaleFreqs
             , L.length thetaFreqs
             ]
             0
  arr <-
    solveMonteCarloR2Z2T0S0
      numThread
      numTrail
      maxTrail
      numPoint
      numPoint
      thetaSigma
      scaleSigma
      maxScale
      tao
      theta0Freqs
      thetaFreqs
      scale0Freqs
      scaleFreqs
      histPath
      hist
  let arr4d =
        R.slice
          arr
          (Z :. (L.length theta0Freqs - 1) :. (L.length scale0Freqs - 1) :. All :.
           All :.
           All :.
           All)
      arr4dMax = VU.maximum . VU.map magnitude . toUnboxed . computeS $ arr4d
      -- arr4dNormalized = R.map (\x -> x / (arr4dMax :+ 0) * (255 :+ 0))
  -- MP.mapM_
  --   (\(i, j) ->
  --      plotImageRepaComplex
  --        (folderPath </> "GreensR2Z2T0S0" </> show (i + 1) L.++ "_" L.++
  --         show (j + 1) L.++
  --         ".png") .
  --      ImageRepa 8 .
  --      computeS . R.extend (Z :. (1 :: Int) :. All :. All) . R.slice arr4dNormalized $
  --      (Z :. i :. j :. All :. All))
  --   [ (i, j)
  --   | i <- [0 .. (L.length thetaFreqs) - 1]
  --   , j <- [0 .. (L.length scaleFreqs) - 1]
  --   ]
  flag <- doesFileExist radialHistPath
  radialArr <-
    if flag
      then R.map magnitude . getNormalizedHistogramArr <$>
           decodeFile radialHistPath
      else do
        putStrLn "Couldn't find a Green's function data. Start simulation..."
        solveMonteCarloR2Z2T0S0Radial
          numThread
          numTrail
          maxTrail
          numPoint
          numPoint
          thetaSigma
          scaleSigma
          maxScale
          tao
          theta0Freqs
          thetaFreqs
          scale0Freqs
          scaleFreqs
          radialHistPath
          (emptyHistogram
             [ if maxScale >= 2
                 then round maxScale + 1
                 else (round . sqrt . fromIntegral $ 2 * (div numPoint 2) ^ 2)
             , L.length scale0Freqs
             , L.length theta0Freqs
             , L.length scaleFreqs
             , L.length thetaFreqs
             ]
             0)
  arr' <-
    computeUnboxedP $
    computeR2Z2T0S0ArrayRadial
      (PinwheelHollow0 1)
      radialArr
      numPoint
      numPoint
      1
      maxScale
      thetaFreqs
      scaleFreqs
      theta0Freqs
      scale0Freqs
  let arr4d' =
        R.slice
          arr'
          (Z :. All :. All :. (L.length thetaFreqs - 1) :.
           (L.length scaleFreqs - 1) :.
           All :.
           All)
  MP.mapM_
    (\(i, j) ->
       plotImageRepaComplex
         (folderPath </> "GreensR2Z2T0S0" </> "Radial_" L.++ show (i + 1) L.++
          "_" L.++
          show (j + 1) L.++
          ".png") .
       ImageRepa 8 .
       computeS . R.extend (Z :. (1 :: Int) :. All :. All) . R.slice arr4d' $
       (Z :. i :. j :. All :. All))
    [ (i, j)
    | i <- [0 .. (L.length thetaFreqs) - 1]
    , j <- [0 .. (L.length scaleFreqs) - 1]
    ]

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
module Repl.Parser (parseExpr) where

import Data.Array
import Data.Char
import Data.Complex
import Data.Maybe
import Data.Ratio
import Numeric
import Text.ParserCombinators.Parsec

import Core.Language

parseExpr :: Parser SchemeVal
parseExpr = tryAll parsers

parsers = [parseAtom,
           parseAtomEnclosed,
           parseBoolean,
           parseChar,
           parseHexChar,
           parseSpecialChar,
           parseString,
           parseComplex,
           parseReal,
           parseRational,
           parseInteger,
           parseList,
           parseImproperList,
           parseVector,
           parseBytevector,
           parseQuote,
           parseUnquote,
           parseQuasiquote]

-- PARSERS ---------------------------------------------------------------------

parseAtom :: Parser SchemeVal
parseAtom = do
  str <- many1 $ letter <|> digit <|> symbol
  if isDigit (head str) then fail "digit not allowed" else return $ Atom str

parseAtomEnclosed :: Parser SchemeVal
parseAtomEnclosed = do
  char '|'
  str <- many1 $ letter <|> digit <|> symbol <|> hexChar <|> space
  char '|'
  return $ Atom str

parseBoolean :: Parser SchemeVal
parseBoolean = do
  char '#'
  c <- oneOf "tf"
  return $ Boolean $ c == 't'

parseChar :: Parser SchemeVal
parseChar = do
  string "#\\"
  x <- anyChar
  notFollowedBy alphaNum
  return $ Char x

parseHexChar :: Parser SchemeVal
parseHexChar = fmap (Char . chr . fromHex) (string "#\\x" >> many1 hexDigit)

parseSpecialChar :: Parser SchemeVal
parseSpecialChar = do
  string "#\\"
  str <- tryAll $ map (string . fst) specialChars
  return $ Char $ fromJust $ lookup str specialChars

parseString :: Parser SchemeVal
parseString = do
  char '"'
  str <- many $ tryAll [escChar, hexChar, noneOf "\"\\"]
  char '"'
  return $ String str

parseInteger :: Parser SchemeVal
parseInteger = Integer <$> tryAll [int, hex, dec, oct, bin]

parseRational :: Parser SchemeVal
parseRational = do
  a <- many1 digit
  char '/'
  b <- many1 digit
  return $ Rational $ read a % read b

parseReal :: Parser SchemeVal
parseReal = do
  a <- many1 digit
  char '.'
  b <- many1 digit
  return $ Real $ fst . head . readFloat $ a ++ "." ++ b

parseComplex :: Parser SchemeVal
parseComplex = do
  x <- tryAll [parseReal, parseInteger]
  char '+'
  y <- tryAll [parseReal, parseInteger]
  char 'i'
  return $ Complex $ fromNumber x :+ fromNumber y

parseList :: Parser SchemeVal
parseList = do
  char '('
  xs <- sepBy parseExpr spaces
  char ')'
  return $ List xs

parseImproperList :: Parser SchemeVal
parseImproperList = do
  char '('
  xs <- endBy parseExpr spaces
  char '.'
  x <- spaces >> parseExpr
  char ')'
  return $ ImproperList xs x

parseVector :: Parser SchemeVal
parseVector = do
  string "#("
  xs <- sepBy parseExpr spaces
  char ')'
  return $ Vector $ list2array xs

parseBytevector :: Parser SchemeVal
parseBytevector = do
  string "#u8("
  xs <- sepBy byte spaces
  char ')'
  return $ Bytevector $ list2array xs

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

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

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

-- UTILS -----------------------------------------------------------------------

tryAll = foldl1 (<|>) . map try

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

hexChar = do
  string "\\x"
  x <- many1 hexDigit
  char ';'
  return $ chr $ fromHex x

escChar = do
  char '\\'
  x <- oneOf "\\\"|abnrt"
  return $ case x of
    '\\' -> x
    '"'  -> x
    '|'  -> '|'
    'a'  -> '\a'
    'b'  -> '\b'
    'n'  -> '\n'
    'r'  -> '\r'
    't'  -> '\t'

specialChars = [("alarm", '\a'), ("backspace", '\b'), ("delete", '\DEL'),
                ("escape", '\ESC'), ("newline", '\n'), ("null", '\NUL'),
                ("return", '\r'), ("space", ' '), ("tab", '\t')]

int = read <$> many1 digit

hex = fmap fromHex $ string "#x" >> many1 hexDigit
dec = fmap fromDec $ string "#d" >> many1 decDigit
oct = fmap fromOct $ string "#o" >> many1 octDigit
bin = fmap fromBin $ string "#b" >> many1 binDigit

fromHex = fst . head . readHex
fromDec = read
fromOct = fst . head . readOct
fromBin = foldl1 (\ a b -> 2 * a + b) . map digitToInt

decDigit = digit

binDigit = oneOf "01"

fromNumber :: SchemeVal -> Double
fromNumber (Real x) = x
fromNumber (Integer x) = fromIntegral x

list2array xs = listArray (0, length xs - 1) xs

byte = do
  x <- read <$> many1 digit
  if 0 <= x && x < 256 then return x else fail $ show x ++ " is not a byte"

max_stars_repo_name: calintat/scheme

max_issues_repo_name: calintat/scheme

max_forks_repo_name: calintat/scheme