text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE MagicHash #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE ViewPatterns #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE DefaultSignatures #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
-----------------------------------------------------------------------------
-- |
-- Copyright   :  (C) 2013-2015 Edward Kmett
-- License     :  BSD-style (see the file LICENSE)
-- Maintainer  :  Edward Kmett <ekmett@gmail.com>
-- Stability   :  experimental
-- Portability :  non-portable
--
-- This module provides a choice of a "best" vector type for a given type
-- that is as unboxed as possible.
-----------------------------------------------------------------------------
module Data.Vector.Array

  ( Arrayed(..)
  , Array
  , MArray
  -- * Internals
  , V_Complex(V_Complex)
  , MV_Complex(MV_Complex)
  , V_Pair(V_Pair)
  , MV_Pair(MV_Pair)
  , prefetchArray0#, prefetchArray1#, prefetchArray2#, prefetchArray3#
  , prefetchSmallArray0#, prefetchSmallArray1#, prefetchSmallArray2#, prefetchSmallArray3#
  , prefetchMutableArray0#, prefetchMutableArray1#, prefetchMutableArray2#, prefetchMutableArray3#
  , prefetchSmallMutableArray0#, prefetchSmallMutableArray1#, prefetchSmallMutableArray2#, prefetchSmallMutableArray3#
  , prefetchPrim0#, prefetchPrim1#, prefetchPrim2#, prefetchPrim3#
  , prefetchMutablePrim0#, prefetchMutablePrim1#, prefetchMutablePrim2#, prefetchMutablePrim3#
  ) where

import Control.Monad
import Data.Complex
import Data.Int
import Data.Primitive hiding (Array)
import qualified Data.Primitive as Primitive
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Generic.Mutable as GM
import qualified Data.Vector.Primitive as P
import qualified Data.Vector.Primitive.Mutable as PM
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Unboxed.Base as UB
import qualified Data.Vector.Fusion.Stream as Stream
import qualified Data.Vector as B
import qualified Data.Vector.Mutable as BM
import Data.Word
import GHC.Prim
import GHC.Types
import Text.Read

#include "ghcautoconf.h"
#include "DerivedConstants.h"

#if SIZEOF_VOID_P == 4
#define W 4#
#elif SIZEOF_VOID_P == 8
#define W 8#
#else
#error unknown pointer size
#endif

#if OFFSET_StgMutArrPtrs_size != SIZEOF_VOID_P
#error did the card marking machinery change?
#endif

-- | A vector of product-like data types that know how to store themselves in a Vector optimally,
-- maximizing the level of unboxing provided, but not guaranteeing to unbox it all.

type Array a = Arr a a
type MArray s a = G.Mutable (Arr a) s a

prefetchArray0#, prefetchArray1#, prefetchArray2#, prefetchArray3# :: Array# a -> Int# -> State# d -> State# d
prefetchArray0# a i s = unsafeCoerce# prefetchByteArray0# a ((i *# W) +# W) s
prefetchArray1# a i s = unsafeCoerce# prefetchByteArray1# a ((i *# W) +# W) s
prefetchArray2# a i s = unsafeCoerce# prefetchByteArray2# a ((i *# W) +# W) s
prefetchArray3# a i s = unsafeCoerce# prefetchByteArray3# a ((i *# W) +# W) s

prefetchMutableArray0#, prefetchMutableArray1#, prefetchMutableArray2#, prefetchMutableArray3# :: MutableArray# d a -> Int# -> State# d -> State# d
prefetchMutableArray0# a i s = unsafeCoerce# prefetchMutableByteArray0# a ((i *# W) +# W) s
prefetchMutableArray1# a i s = unsafeCoerce# prefetchMutableByteArray1# a ((i *# W) +# W) s
prefetchMutableArray2# a i s = unsafeCoerce# prefetchMutableByteArray2# a ((i *# W) +# W) s
prefetchMutableArray3# a i s = unsafeCoerce# prefetchMutableByteArray3# a ((i *# W) +# W) s

prefetchSmallArray0#, prefetchSmallArray1#, prefetchSmallArray2#, prefetchSmallArray3# :: SmallArray# a -> Int# -> State# d -> State# d
prefetchSmallArray0# a i s = unsafeCoerce# prefetchByteArray0# a (i *# W) s
prefetchSmallArray1# a i s = unsafeCoerce# prefetchByteArray1# a (i *# W) s
prefetchSmallArray2# a i s = unsafeCoerce# prefetchByteArray2# a (i *# W) s
prefetchSmallArray3# a i s = unsafeCoerce# prefetchByteArray3# a (i *# W) s

prefetchSmallMutableArray0#, prefetchSmallMutableArray1#, prefetchSmallMutableArray2#, prefetchSmallMutableArray3# :: SmallMutableArray# d a -> Int# -> State# d -> State# d
prefetchSmallMutableArray0# a i s = unsafeCoerce# prefetchMutableByteArray0# a (i *# W) s
prefetchSmallMutableArray1# a i s = unsafeCoerce# prefetchMutableByteArray1# a (i *# W) s
prefetchSmallMutableArray2# a i s = unsafeCoerce# prefetchMutableByteArray2# a (i *# W) s
prefetchSmallMutableArray3# a i s = unsafeCoerce# prefetchMutableByteArray3# a (i *# W) s

class (G.Vector (Arr a) a, Monoid (Arr a a)) => Arrayed a where
  type Arr a :: * -> *
  type Arr a = B.Vector

  prefetchArr0# :: Arr a a -> Int# -> State# s -> State# s
  default prefetchArr0# :: (Arr a ~ B.Vector) => Arr a a -> Int# -> State# s -> State# s
  prefetchArr0# (coerceVector -> BVector (I# o) _ (Primitive.Array arr)) i s = prefetchArray0# arr (o +# i) s

  prefetchArr1# :: Arr a a -> Int# -> State# s -> State# s
  default prefetchArr1# :: (Arr a ~ B.Vector) => Arr a a -> Int# -> State# s -> State# s
  prefetchArr1# (coerceVector -> BVector (I# o) _ (Primitive.Array arr)) i s = prefetchArray1# arr (o +# i) s

  prefetchArr2# :: Arr a a -> Int# -> State# s -> State# s
  default prefetchArr2# :: (Arr a ~ B.Vector) => Arr a a -> Int# -> State# s -> State# s
  prefetchArr2# (coerceVector -> BVector (I# o) _ (Primitive.Array arr)) i s = prefetchArray2# arr (o +# i) s

  prefetchArr3# :: Arr a a -> Int# -> State# s -> State# s
  default prefetchArr3# :: (Arr a ~ B.Vector) => Arr a a -> Int# -> State# s -> State# s
  prefetchArr3# (coerceVector -> BVector (I# o) _ (Primitive.Array arr)) i s = prefetchArray3# arr (o +# i) s

  prefetchMutableArr0# :: G.Mutable (Arr a) s a -> Int# -> State# s -> State# s
  default prefetchMutableArr0# :: (G.Mutable (Arr a) ~ BM.MVector) => G.Mutable (Arr a) s a -> Int# -> State# s -> State# s
  prefetchMutableArr0# (BM.MVector (I# o) _ (MutableArray arr)) i s = prefetchMutableArray0# arr (o +# i) s

  prefetchMutableArr1# :: G.Mutable (Arr a) s a -> Int# -> State# s -> State# s
  default prefetchMutableArr1# :: (G.Mutable (Arr a) ~ BM.MVector) => G.Mutable (Arr a) s a -> Int# -> State# s -> State# s
  prefetchMutableArr1# (BM.MVector (I# o) _ (MutableArray arr)) i s = prefetchMutableArray1# arr (o +# i) s

  prefetchMutableArr2# :: G.Mutable (Arr a) s a -> Int# -> State# s -> State# s
  default prefetchMutableArr2# :: (G.Mutable (Arr a) ~ BM.MVector) => G.Mutable (Arr a) s a -> Int# -> State# s -> State# s
  prefetchMutableArr2# (BM.MVector (I# o) _ (MutableArray arr)) i s = prefetchMutableArray2# arr (o +# i) s

  prefetchMutableArr3# :: G.Mutable (Arr a) s a -> Int# -> State# s -> State# s
  default prefetchMutableArr3# :: (G.Mutable (Arr a) ~ BM.MVector) => G.Mutable (Arr a) s a -> Int# -> State# s -> State# s
  prefetchMutableArr3# (BM.MVector (I# o) _ (MutableArray arr)) i s = prefetchMutableArray3# arr (o +# i) s

prefetchPrim0#, prefetchPrim1#, prefetchPrim2#, prefetchPrim3# :: forall s a. Prim a => P.Vector a -> Int# -> State# s -> State# s
prefetchPrim0# (P.Vector (I# o) _ (ByteArray arr)) i s = prefetchByteArray0# arr ((o +# i) *# sizeOf# (undefined :: a)) s
prefetchPrim1# (P.Vector (I# o) _ (ByteArray arr)) i s = prefetchByteArray1# arr ((o +# i) *# sizeOf# (undefined :: a)) s
prefetchPrim2# (P.Vector (I# o) _ (ByteArray arr)) i s = prefetchByteArray2# arr ((o +# i) *# sizeOf# (undefined :: a)) s
prefetchPrim3# (P.Vector (I# o) _ (ByteArray arr)) i s = prefetchByteArray3# arr ((o +# i) *# sizeOf# (undefined :: a)) s

prefetchMutablePrim0#, prefetchMutablePrim1#, prefetchMutablePrim2#, prefetchMutablePrim3# :: forall s a. Prim a => PM.MVector s a -> Int# -> State# s -> State# s
prefetchMutablePrim0# (PM.MVector (I# o) _ (MutableByteArray arr)) i s = prefetchMutableByteArray0# arr ((o +# i) *# sizeOf# (undefined :: a)) s
prefetchMutablePrim1# (PM.MVector (I# o) _ (MutableByteArray arr)) i s = prefetchMutableByteArray1# arr ((o +# i) *# sizeOf# (undefined :: a)) s
prefetchMutablePrim2# (PM.MVector (I# o) _ (MutableByteArray arr)) i s = prefetchMutableByteArray2# arr ((o +# i) *# sizeOf# (undefined :: a)) s
prefetchMutablePrim3# (PM.MVector (I# o) _ (MutableByteArray arr)) i s = prefetchMutableByteArray3# arr ((o +# i) *# sizeOf# (undefined :: a)) s

-- * Unboxed vectors

instance Arrayed () where
  type Arr () = U.Vector
  prefetchArr0# _ _ s = s
  prefetchArr1# _ _ s = s
  prefetchArr2# _ _ s = s
  prefetchArr3# _ _ s = s
  prefetchMutableArr0# _ _ s = s
  prefetchMutableArr1# _ _ s = s
  prefetchMutableArr2# _ _ s = s
  prefetchMutableArr3# _ _ s = s

instance Arrayed Double where
  type Arr Double = U.Vector
  prefetchArr0# (UB.V_Double v) i s = prefetchPrim0# v i s
  prefetchArr1# (UB.V_Double v) i s = prefetchPrim1# v i s
  prefetchArr2# (UB.V_Double v) i s = prefetchPrim2# v i s
  prefetchArr3# (UB.V_Double v) i s = prefetchPrim3# v i s
  prefetchMutableArr0# (UB.MV_Double v) i s = prefetchMutablePrim0# v i s
  prefetchMutableArr1# (UB.MV_Double v) i s = prefetchMutablePrim1# v i s
  prefetchMutableArr2# (UB.MV_Double v) i s = prefetchMutablePrim2# v i s
  prefetchMutableArr3# (UB.MV_Double v) i s = prefetchMutablePrim3# v i s

instance Arrayed Float where
  type Arr Float = U.Vector
  prefetchArr0# (UB.V_Float v) i s = prefetchPrim0# v i s
  prefetchArr1# (UB.V_Float v) i s = prefetchPrim1# v i s
  prefetchArr2# (UB.V_Float v) i s = prefetchPrim2# v i s
  prefetchArr3# (UB.V_Float v) i s = prefetchPrim3# v i s
  prefetchMutableArr0# (UB.MV_Float v) i s = prefetchMutablePrim0# v i s
  prefetchMutableArr1# (UB.MV_Float v) i s = prefetchMutablePrim1# v i s
  prefetchMutableArr2# (UB.MV_Float v) i s = prefetchMutablePrim2# v i s
  prefetchMutableArr3# (UB.MV_Float v) i s = prefetchMutablePrim3# v i s

instance Arrayed Int where
  type Arr Int = U.Vector
  prefetchArr0# (UB.V_Int v) i s = prefetchPrim0# v i s
  prefetchArr1# (UB.V_Int v) i s = prefetchPrim1# v i s
  prefetchArr2# (UB.V_Int v) i s = prefetchPrim2# v i s
  prefetchArr3# (UB.V_Int v) i s = prefetchPrim3# v i s
  prefetchMutableArr0# (UB.MV_Int v) i s = prefetchMutablePrim0# v i s
  prefetchMutableArr1# (UB.MV_Int v) i s = prefetchMutablePrim1# v i s
  prefetchMutableArr2# (UB.MV_Int v) i s = prefetchMutablePrim2# v i s
  prefetchMutableArr3# (UB.MV_Int v) i s = prefetchMutablePrim3# v i s

instance Arrayed Int8 where
  type Arr Int8 = U.Vector
  prefetchArr0# (UB.V_Int8 v) i s = prefetchPrim0# v i s
  prefetchArr1# (UB.V_Int8 v) i s = prefetchPrim1# v i s
  prefetchArr2# (UB.V_Int8 v) i s = prefetchPrim2# v i s
  prefetchArr3# (UB.V_Int8 v) i s = prefetchPrim3# v i s
  prefetchMutableArr0# (UB.MV_Int8 v) i s = prefetchMutablePrim0# v i s
  prefetchMutableArr1# (UB.MV_Int8 v) i s = prefetchMutablePrim1# v i s
  prefetchMutableArr2# (UB.MV_Int8 v) i s = prefetchMutablePrim2# v i s
  prefetchMutableArr3# (UB.MV_Int8 v) i s = prefetchMutablePrim3# v i s

instance Arrayed Int16 where
  type Arr Int16 = U.Vector
  prefetchArr0# (UB.V_Int16 v) i s = prefetchPrim0# v i s
  prefetchArr1# (UB.V_Int16 v) i s = prefetchPrim1# v i s
  prefetchArr2# (UB.V_Int16 v) i s = prefetchPrim2# v i s
  prefetchArr3# (UB.V_Int16 v) i s = prefetchPrim3# v i s
  prefetchMutableArr0# (UB.MV_Int16 v) i s = prefetchMutablePrim0# v i s
  prefetchMutableArr1# (UB.MV_Int16 v) i s = prefetchMutablePrim1# v i s
  prefetchMutableArr2# (UB.MV_Int16 v) i s = prefetchMutablePrim2# v i s
  prefetchMutableArr3# (UB.MV_Int16 v) i s = prefetchMutablePrim3# v i s

instance Arrayed Int32 where
  type Arr Int32 = U.Vector
  prefetchArr0# (UB.V_Int32 v) i s = prefetchPrim0# v i s
  prefetchArr1# (UB.V_Int32 v) i s = prefetchPrim1# v i s
  prefetchArr2# (UB.V_Int32 v) i s = prefetchPrim2# v i s
  prefetchArr3# (UB.V_Int32 v) i s = prefetchPrim3# v i s
  prefetchMutableArr0# (UB.MV_Int32 v) i s = prefetchMutablePrim0# v i s
  prefetchMutableArr1# (UB.MV_Int32 v) i s = prefetchMutablePrim1# v i s
  prefetchMutableArr2# (UB.MV_Int32 v) i s = prefetchMutablePrim2# v i s
  prefetchMutableArr3# (UB.MV_Int32 v) i s = prefetchMutablePrim3# v i s

instance Arrayed Int64 where
  type Arr Int64 = U.Vector
  prefetchArr0# (UB.V_Int64 v) i s = prefetchPrim0# v i s
  prefetchArr1# (UB.V_Int64 v) i s = prefetchPrim1# v i s
  prefetchArr2# (UB.V_Int64 v) i s = prefetchPrim2# v i s
  prefetchArr3# (UB.V_Int64 v) i s = prefetchPrim3# v i s
  prefetchMutableArr0# (UB.MV_Int64 v) i s = prefetchMutablePrim0# v i s
  prefetchMutableArr1# (UB.MV_Int64 v) i s = prefetchMutablePrim1# v i s
  prefetchMutableArr2# (UB.MV_Int64 v) i s = prefetchMutablePrim2# v i s
  prefetchMutableArr3# (UB.MV_Int64 v) i s = prefetchMutablePrim3# v i s

instance Arrayed Word where
  type Arr Word = U.Vector
  prefetchArr0# (UB.V_Word v) i s = prefetchPrim0# v i s
  prefetchArr1# (UB.V_Word v) i s = prefetchPrim1# v i s
  prefetchArr2# (UB.V_Word v) i s = prefetchPrim2# v i s
  prefetchArr3# (UB.V_Word v) i s = prefetchPrim3# v i s
  prefetchMutableArr0# (UB.MV_Word v) i s = prefetchMutablePrim0# v i s
  prefetchMutableArr1# (UB.MV_Word v) i s = prefetchMutablePrim1# v i s
  prefetchMutableArr2# (UB.MV_Word v) i s = prefetchMutablePrim2# v i s
  prefetchMutableArr3# (UB.MV_Word v) i s = prefetchMutablePrim3# v i s

instance Arrayed Word8 where
  type Arr Word8 = U.Vector
  prefetchArr0# (UB.V_Word8 v) i s = prefetchPrim0# v i s
  prefetchArr1# (UB.V_Word8 v) i s = prefetchPrim1# v i s
  prefetchArr2# (UB.V_Word8 v) i s = prefetchPrim2# v i s
  prefetchArr3# (UB.V_Word8 v) i s = prefetchPrim3# v i s
  prefetchMutableArr0# (UB.MV_Word8 v) i s = prefetchMutablePrim0# v i s
  prefetchMutableArr1# (UB.MV_Word8 v) i s = prefetchMutablePrim1# v i s
  prefetchMutableArr2# (UB.MV_Word8 v) i s = prefetchMutablePrim2# v i s
  prefetchMutableArr3# (UB.MV_Word8 v) i s = prefetchMutablePrim3# v i s

instance Arrayed Word16 where
  type Arr Word16 = U.Vector
  prefetchArr0# (UB.V_Word16 v) i s = prefetchPrim0# v i s
  prefetchArr1# (UB.V_Word16 v) i s = prefetchPrim1# v i s
  prefetchArr2# (UB.V_Word16 v) i s = prefetchPrim2# v i s
  prefetchArr3# (UB.V_Word16 v) i s = prefetchPrim3# v i s
  prefetchMutableArr0# (UB.MV_Word16 v) i s = prefetchMutablePrim0# v i s
  prefetchMutableArr1# (UB.MV_Word16 v) i s = prefetchMutablePrim1# v i s
  prefetchMutableArr2# (UB.MV_Word16 v) i s = prefetchMutablePrim2# v i s
  prefetchMutableArr3# (UB.MV_Word16 v) i s = prefetchMutablePrim3# v i s

instance Arrayed Word32 where
  type Arr Word32 = U.Vector
  prefetchArr0# (UB.V_Word32 v) i s = prefetchPrim0# v i s
  prefetchArr1# (UB.V_Word32 v) i s = prefetchPrim1# v i s
  prefetchArr2# (UB.V_Word32 v) i s = prefetchPrim2# v i s
  prefetchArr3# (UB.V_Word32 v) i s = prefetchPrim3# v i s
  prefetchMutableArr0# (UB.MV_Word32 v) i s = prefetchMutablePrim0# v i s
  prefetchMutableArr1# (UB.MV_Word32 v) i s = prefetchMutablePrim1# v i s
  prefetchMutableArr2# (UB.MV_Word32 v) i s = prefetchMutablePrim2# v i s
  prefetchMutableArr3# (UB.MV_Word32 v) i s = prefetchMutablePrim3# v i s

instance Arrayed Word64 where
  type Arr Word64 = U.Vector
  prefetchArr0# (UB.V_Word64 v) i s = prefetchPrim0# v i s
  prefetchArr1# (UB.V_Word64 v) i s = prefetchPrim1# v i s
  prefetchArr2# (UB.V_Word64 v) i s = prefetchPrim2# v i s
  prefetchArr3# (UB.V_Word64 v) i s = prefetchPrim3# v i s
  prefetchMutableArr0# (UB.MV_Word64 v) i s = prefetchMutablePrim0# v i s
  prefetchMutableArr1# (UB.MV_Word64 v) i s = prefetchMutablePrim1# v i s
  prefetchMutableArr2# (UB.MV_Word64 v) i s = prefetchMutablePrim2# v i s
  prefetchMutableArr3# (UB.MV_Word64 v) i s = prefetchMutablePrim3# v i s

-- * Boxed vectors

instance Arrayed Integer
instance Arrayed [a]
instance Arrayed (Maybe a)
instance Arrayed (Either a b)
instance Arrayed (IO a)

-- * Pairs are boxed or unboxed based on their components

#ifndef HLINT
data MV_Pair :: * -> * -> * where
  MV_Pair:: {-# UNPACK #-} !Int -> !(G.Mutable (Arr a) s a) -> !(G.Mutable (Arr b) s b) -> MV_Pair s (a, b)

data V_Pair :: * -> * where
  V_Pair :: {-# UNPACK #-} !Int -> !(Array a) -> !(Array b) -> V_Pair (a, b)
#endif

type instance G.Mutable V_Pair = MV_Pair

instance (Arrayed a, Arrayed b) => GM.MVector MV_Pair (a, b) where
  {-# INLINE basicLength #-}
  {-# INLINE basicUnsafeSlice #-}
  {-# INLINE basicOverlaps #-}
  {-# INLINE basicUnsafeNew #-}
  {-# INLINE basicUnsafeReplicate #-}
  {-# INLINE basicUnsafeRead #-}
  {-# INLINE basicUnsafeWrite #-}
  {-# INLINE basicClear #-}
  {-# INLINE basicSet #-}
  {-# INLINE basicUnsafeCopy #-}
  {-# INLINE basicUnsafeGrow #-}
  basicLength (MV_Pair l _ _) = l
  basicUnsafeSlice i n (MV_Pair _ u v)                   = MV_Pair n (GM.basicUnsafeSlice i n u) (GM.basicUnsafeSlice i n v)
  basicOverlaps (MV_Pair _ u1 v1) (MV_Pair _ u2 v2)      = GM.basicOverlaps u1 u2 || GM.basicOverlaps v1 v2
  basicUnsafeNew n                                       = liftM2 (MV_Pair n) (GM.basicUnsafeNew n) (GM.basicUnsafeNew n)
  basicUnsafeReplicate n (x, y)                          = liftM2 (MV_Pair n) (GM.basicUnsafeReplicate n x) (GM.basicUnsafeReplicate n y)
  basicUnsafeRead (MV_Pair _ u v) i                      = liftM2 (,) (GM.basicUnsafeRead u i) (GM.basicUnsafeRead v i)
  basicUnsafeWrite (MV_Pair _ u v) i (x, y)              = GM.basicUnsafeWrite u i x >> GM.basicUnsafeWrite v i y
  basicClear (MV_Pair _ u v)                             = GM.basicClear u >> GM.basicClear v
  basicSet (MV_Pair _ u v) (x, y)                        = GM.basicSet u x >> GM.basicSet v y
  basicUnsafeCopy (MV_Pair _ u1 v1) (MV_Pair _ u2 v2)    = GM.basicUnsafeCopy u1 u2 >> GM.basicUnsafeCopy v1 v2
  basicUnsafeMove (MV_Pair _ u1 v1) (MV_Pair _ u2 v2)    = GM.basicUnsafeMove u1 u2 >> GM.basicUnsafeMove v1 v2
  basicUnsafeGrow (MV_Pair _ u v) n                      = liftM2 (MV_Pair n) (GM.basicUnsafeGrow u n) (GM.basicUnsafeGrow v n)

instance (Arrayed a, Arrayed b) => G.Vector V_Pair (a, b) where
  {-# INLINE basicLength #-}
  {-# INLINE basicUnsafeFreeze #-}
  {-# INLINE basicUnsafeThaw #-}
  {-# INLINE basicUnsafeSlice #-}
  {-# INLINE basicUnsafeIndexM #-}
  {-# INLINE elemseq #-}
  basicLength (V_Pair v _ _) = v
  basicUnsafeFreeze (MV_Pair n u v)                = liftM2 (V_Pair n) (G.basicUnsafeFreeze u) (G.basicUnsafeFreeze v)
  basicUnsafeThaw (V_Pair n u v)                   = liftM2 (MV_Pair n) (G.basicUnsafeThaw u) (G.basicUnsafeThaw v)
  basicUnsafeSlice i n (V_Pair _ u v)              = V_Pair n (G.basicUnsafeSlice i n u) (G.basicUnsafeSlice i n v)
  basicUnsafeIndexM (V_Pair _ u v) i               = liftM2 (,) (G.basicUnsafeIndexM u i) (G.basicUnsafeIndexM v i)
  basicUnsafeCopy (MV_Pair _ mu mv) (V_Pair _ u v) = G.basicUnsafeCopy mu u >> G.basicUnsafeCopy mv v
  elemseq _ (x, y) z  = G.elemseq (undefined :: Array a) x
                       $ G.elemseq (undefined :: Array b) y z

instance (Arrayed a, Arrayed b, Show a, Show b, c ~ (a, b)) => Show (V_Pair c) where
  showsPrec = G.showsPrec

instance (Arrayed a, Arrayed b, Read a, Read b, c ~ (a, b)) => Read (V_Pair c) where
  readPrec = G.readPrec
  readListPrec = readListPrecDefault

instance (Arrayed a, Arrayed b, Eq a, Eq b, c ~ (a, b)) => Eq (V_Pair c) where
  xs == ys = Stream.eq (G.stream xs) (G.stream ys)
  {-# INLINE (==) #-}

instance (Arrayed a, Arrayed b, c ~ (a, b)) => Monoid (V_Pair c) where
  mappend = (G.++)
  {-# INLINE mappend #-}
  mempty = G.empty
  {-# INLINE mempty #-}
  mconcat = G.concat
  {-# INLINE mconcat #-}

instance (Arrayed a, Arrayed b) => Arrayed (a, b) where
  type Arr (a, b) = V_Pair
  prefetchArr0# (V_Pair _ u v) i s = prefetchArr0# v i (prefetchArr0# u i s)
  prefetchArr1# (V_Pair _ u v) i s = prefetchArr1# v i (prefetchArr1# u i s)
  prefetchArr2# (V_Pair _ u v) i s = prefetchArr2# v i (prefetchArr2# u i s)
  prefetchArr3# (V_Pair _ u v) i s = prefetchArr3# v i (prefetchArr3# u i s)
  prefetchMutableArr0# (MV_Pair _ u v) i s = prefetchMutableArr0# v i (prefetchMutableArr0# u i s)
  prefetchMutableArr1# (MV_Pair _ u v) i s = prefetchMutableArr1# v i (prefetchMutableArr1# u i s)
  prefetchMutableArr2# (MV_Pair _ u v) i s = prefetchMutableArr2# v i (prefetchMutableArr2# u i s)
  prefetchMutableArr3# (MV_Pair _ u v) i s = prefetchMutableArr3# v i (prefetchMutableArr3# u i s)

-- * Complex numbers are boxed or unboxed based on their components

#ifndef HLINT
data MV_Complex :: * -> * -> * where
  MV_Complex :: {-# UNPACK #-} !Int -> !(G.Mutable (Arr a) s a) -> !(G.Mutable (Arr a) s a) -> MV_Complex s (Complex a)

data V_Complex :: * -> * where
  V_Complex :: {-# UNPACK #-} !Int -> !(Array a) -> !(Array a) -> V_Complex (Complex a)
#endif

type instance G.Mutable V_Complex = MV_Complex

instance (Arrayed a, RealFloat a) => GM.MVector MV_Complex (Complex a) where
  {-# INLINE basicLength #-}
  {-# INLINE basicUnsafeSlice #-}
  {-# INLINE basicOverlaps #-}
  {-# INLINE basicUnsafeNew #-}
  {-# INLINE basicUnsafeReplicate #-}
  {-# INLINE basicUnsafeRead #-}
  {-# INLINE basicUnsafeWrite #-}
  {-# INLINE basicClear #-}
  {-# INLINE basicSet #-}
  {-# INLINE basicUnsafeCopy #-}
  {-# INLINE basicUnsafeGrow #-}
  basicLength (MV_Complex l _ _) = l
  basicUnsafeSlice i n (MV_Complex _ u v)                   = MV_Complex n (GM.basicUnsafeSlice i n u) (GM.basicUnsafeSlice i n v)
  basicOverlaps (MV_Complex _ u1 v1) (MV_Complex _ u2 v2)   = GM.basicOverlaps u1 u2 || GM.basicOverlaps v1 v2
  basicUnsafeNew n                                          = liftM2 (MV_Complex n) (GM.basicUnsafeNew n) (GM.basicUnsafeNew n)
  basicUnsafeReplicate n (x :+ y)                           = liftM2 (MV_Complex n) (GM.basicUnsafeReplicate n x) (GM.basicUnsafeReplicate n y)
  basicUnsafeRead (MV_Complex _ u v) i                      = liftM2 (:+) (GM.basicUnsafeRead u i) (GM.basicUnsafeRead v i)
  basicUnsafeWrite (MV_Complex _ u v) i (x :+ y)            = GM.basicUnsafeWrite u i x >> GM.basicUnsafeWrite v i y
  basicClear (MV_Complex _ u v)                             = GM.basicClear u >> GM.basicClear v
  basicSet (MV_Complex _ u v) (x :+ y)                      = GM.basicSet u x >> GM.basicSet v y
  basicUnsafeCopy (MV_Complex _ u1 v1) (MV_Complex _ u2 v2) = GM.basicUnsafeCopy u1 u2 >> GM.basicUnsafeCopy v1 v2
  basicUnsafeMove (MV_Complex _ u1 v1) (MV_Complex _ u2 v2) = GM.basicUnsafeMove u1 u2 >> GM.basicUnsafeMove v1 v2
  basicUnsafeGrow (MV_Complex _ u v) n                      = liftM2 (MV_Complex n) (GM.basicUnsafeGrow u n) (GM.basicUnsafeGrow v n)

instance (Arrayed a, RealFloat a) => G.Vector V_Complex (Complex a) where
  {-# INLINE basicLength #-}
  {-# INLINE basicUnsafeFreeze #-}
  {-# INLINE basicUnsafeThaw #-}
  {-# INLINE basicUnsafeSlice #-}
  {-# INLINE basicUnsafeIndexM #-}
  {-# INLINE elemseq #-}
  basicLength (V_Complex v _ _) = v
  basicUnsafeFreeze (MV_Complex n u v)                   = liftM2 (V_Complex n) (G.basicUnsafeFreeze u) (G.basicUnsafeFreeze v)
  basicUnsafeThaw (V_Complex n u v)                      = liftM2 (MV_Complex n) (G.basicUnsafeThaw u) (G.basicUnsafeThaw v)
  basicUnsafeSlice i n (V_Complex _ u v)                 = V_Complex n (G.basicUnsafeSlice i n u) (G.basicUnsafeSlice i n v)
  basicUnsafeIndexM (V_Complex _ u v) i                  = liftM2 (:+) (G.basicUnsafeIndexM u i) (G.basicUnsafeIndexM v i)
  basicUnsafeCopy (MV_Complex _ mu mv) (V_Complex _ u v) = G.basicUnsafeCopy mu u >> G.basicUnsafeCopy mv v
  elemseq _ (x :+ y) z = G.elemseq (undefined :: Arr a a) x
                       $ G.elemseq (undefined :: Arr a a) y z

instance (Arrayed a, RealFloat a, Show a, b ~ Complex a) => Show (V_Complex b) where
  showsPrec = G.showsPrec

instance (Arrayed a, RealFloat a, Read a, b ~ Complex a) => Read (V_Complex b) where
  readPrec = G.readPrec
  readListPrec = readListPrecDefault

instance (Arrayed a, RealFloat a, Eq a, b ~ Complex a) => Eq (V_Complex b) where
  xs == ys = Stream.eq (G.stream xs) (G.stream ys)
  {-# INLINE (==) #-}

instance (Arrayed a, RealFloat a, b ~ Complex a) => Monoid (V_Complex b) where
  mappend = (G.++)
  {-# INLINE mappend #-}
  mempty = G.empty
  {-# INLINE mempty #-}
  mconcat = G.concat
  {-# INLINE mconcat #-}

instance (Arrayed a, RealFloat a) => Arrayed (Complex a) where
  type Arr (Complex a) = V_Complex
  prefetchArr0# (V_Complex _ u v) i s = prefetchArr0# v i (prefetchArr0# u i s)
  prefetchArr1# (V_Complex _ u v) i s = prefetchArr1# v i (prefetchArr1# u i s)
  prefetchArr2# (V_Complex _ u v) i s = prefetchArr2# v i (prefetchArr2# u i s)
  prefetchArr3# (V_Complex _ u v) i s = prefetchArr3# v i (prefetchArr3# u i s)
  prefetchMutableArr0# (MV_Complex _ u v) i s = prefetchMutableArr0# v i (prefetchMutableArr0# u i s)
  prefetchMutableArr1# (MV_Complex _ u v) i s = prefetchMutableArr1# v i (prefetchMutableArr1# u i s)
  prefetchMutableArr2# (MV_Complex _ u v) i s = prefetchMutableArr2# v i (prefetchMutableArr2# u i s)
  prefetchMutableArr3# (MV_Complex _ u v) i s = prefetchMutableArr3# v i (prefetchMutableArr3# u i s)

-- LAME!
data BVector a
  = BVector {-# UNPACK #-}!Int
            {-# UNPACK #-}!Int
            {-# UNPACK #-}!(Primitive.Array a)

coerceVector :: B.Vector a -> BVector a
coerceVector = unsafeCoerce#

ext: hs

max_stars_repo_name: jmatsushita/structures

max_issues_repo_name: jmatsushita/structures

max_forks_repo_name: jmatsushita/structures


text:
{-# LANGUAGE DeriveDataTypeable #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE MultiParamTypeClasses #-}

module Scheme.Types where

import Control.Exception (Exception, throw)
import Control.Monad.Except (ExceptT)
import Data.Array (Array, elems)
import qualified Data.ByteString as BS
import Data.Complex (Complex ((:+)), imagPart, realPart)
import Data.IORef (IORef)
import Data.Ratio (denominator, numerator)
import Data.Text (Text)
import qualified Data.Text as T
import Data.Typeable (Typeable)
import System.IO (Handle)
import Text.Pretty.Simple (pPrint, pPrintLightBg)
import Text.Show.Functions ()

type Env = IORef [(Text, IORef SchemeVal)]

type SchemeResult = Either SchemeError

type IOSchemeResult = ExceptT SchemeError IO

extractValue :: Either a p -> p
extractValue (Right val) = val
extractValue _ = error "Called with Left"

data Number
  = Integer Integer
  | Real Double
  | Rational Rational
  | Complex (Complex Double)
  deriving (Show, Typeable, Eq)

showNum :: Number -> Text
showNum (Integer i) = T.pack $ show i
showNum (Real d) = T.pack $ show d
showNum (Rational r) = T.pack $ show (numerator r) <> "/" <> show (denominator r)
showNum (Complex p) = T.pack $ show (realPart p) <> "+" <> show (imagPart p) <> "i"

instance Ord Number where
  compare (Integer x) (Integer y) = compare x y
  compare (Integer x) (Real y) = compare (fromInteger x) y
  compare (Integer x) (Rational y) = compare (fromInteger x) y
  compare (Real x) (Integer y) = compare x (fromInteger y)
  compare (Rational x) (Integer y) = compare x (fromInteger y)
  compare (Real x) (Real y) = compare x y
  compare (Real x) (Rational y) = compare (toRational x) y
  compare (Rational x) (Real y) = compare x (toRational y)
  compare (Rational x) (Rational y) = compare x y
  compare _ _ = throw $ InvalidOperation "Cannot compare complex numbers"

instance Num Number where
  fromInteger x = Integer x

  (Integer x) + (Integer y) = Integer (x + y)
  (Integer x) + (Real y) = Real (fromInteger x + y)
  (Integer x) + (Rational y) = Rational (fromInteger x + y)
  (Integer x) + (Complex y) = Complex (fromInteger x + y)
  (Real x) + (Integer y) = Real (x + fromInteger y)
  (Rational x) + (Integer y) = Rational (x + fromInteger y)
  (Complex x) + (Integer y) = Complex (x + fromInteger y)
  (Real x) + (Real y) = Real (x + y)
  (Real x) + (Rational y) = Rational (toRational x + y)
  (Real x) + (Complex y) = Complex ((x :+ 0.0) + y)
  (Rational x) + (Real y) = Rational (x + toRational y)
  (Complex x) + (Real y) = Complex (x + (y :+ 0.0))
  (Rational x) + (Rational y) = Rational (x + y)
  (Complex x) + (Complex y) = Complex (x + y)
  _ + _ = throw $ InvalidOperation "Cannot add rationals and complex numbers"

  (Integer x) * (Integer y) = Integer (x * y)
  (Integer x) * (Real y) = Real (fromInteger x * y)
  (Integer x) * (Rational y) = Rational (fromInteger x * y)
  (Integer x) * (Complex y) = Complex (fromInteger x * y)
  (Real x) * (Integer y) = Real (x * fromInteger y)
  (Rational x) * (Integer y) = Rational (x * fromInteger y)
  (Complex x) * (Integer y) = Complex (x * fromInteger y)
  (Real x) * (Real y) = Real (x * y)
  (Real x) * (Rational y) = Rational (toRational x * y)
  (Real x) * (Complex y) = Complex ((x :+ 0.0) * y)
  (Rational x) * (Real y) = Rational (x * toRational y)
  (Complex x) * (Real y) = Complex (x * (y :+ 0.0))
  (Rational x) * (Rational y) = Rational (x * y)
  (Complex x) * (Complex y) = Complex (x * y)
  _ * _ = throw $ InvalidOperation "Cannot multiply rationals and complex numbers"

  (Integer x) - (Integer y) = Integer (x - y)
  (Integer x) - (Real y) = Real (fromInteger x - y)
  (Integer x) - (Rational y) = Rational (fromInteger x - y)
  (Integer x) - (Complex y) = Complex (fromInteger x - y)
  (Real x) - (Integer y) = Real (x - fromInteger y)
  (Rational x) - (Integer y) = Rational (x - fromInteger y)
  (Complex x) - (Integer y) = Complex (x - fromInteger y)
  (Real x) - (Real y) = Real (x - y)
  (Real x) - (Rational y) = Rational (toRational x - y)
  (Real x) - (Complex y) = Complex ((x :+ 0.0) - y)
  (Rational x) - (Real y) = Rational (x - toRational y)
  (Complex x) - (Real y) = Complex (x - (y :+ 0.0))
  (Rational x) - (Rational y) = Rational (x - y)
  (Complex x) - (Complex y) = Complex (x - y)
  _ - _ = throw $ InvalidOperation "Cannot subtract rationals and complex numbers"

  negate (Integer x) = Integer (negate x)
  negate (Real x) = Real (negate x)
  negate (Rational x) = Rational (negate x)
  negate (Complex x) = Complex (negate x)

  abs (Integer x) = Integer (abs x)
  abs (Real x) = Real (abs x)
  abs (Rational x) = Rational (abs x)
  abs (Complex x) = Complex (abs x)

  signum (Integer x) = Integer (signum x)
  signum (Real x) = Real (signum x)
  signum (Rational x) = Rational (signum x)
  signum (Complex x) = Complex (signum x)

instance Fractional Number where
  fromRational x = Rational x

  (Integer x) / (Integer y) = Rational (fromInteger x / fromInteger y)
  (Integer x) / (Real y) = Real (fromInteger x / y)
  (Integer x) / (Rational y) = Rational (fromInteger x / y)
  (Integer x) / (Complex y) = Complex (fromInteger x / y)
  (Real x) / (Integer y) = Real (x / fromInteger y)
  (Rational x) / (Integer y) = Rational (x / fromInteger y)
  (Complex x) / (Integer y) = Complex (x / fromInteger y)
  (Real x) / (Real y) = Real (x / y)
  (Real x) / (Rational y) = Rational (toRational x / y)
  (Real x) / (Complex y) = Complex ((x :+ 0.0) / y)
  (Rational x) / (Real y) = Rational (x / toRational y)
  (Complex x) / (Real y) = Complex (x / (y :+ 0.0))
  (Rational x) / (Rational y) = Rational (x / y)
  (Complex x) / (Complex y) = Complex (x / y)
  _ / _ = throw $ InvalidOperation "Cannot divide rationals and complex numbers"

  recip (Integer x) = Real (recip $ fromInteger x)
  recip (Real x) = Real (recip x)
  recip (Rational x) = Rational (recip x)
  recip (Complex x) = Complex (recip x)

data Fn = Fn {macro :: Bool, params :: [Text], vararg :: Maybe Text, body :: [SchemeVal], closure :: Env}

instance Show Fn where
  show _ = "<fn>"

data SchemeVal
  = List [SchemeVal]
  | PairList [SchemeVal] SchemeVal
  | Vector (Array Int SchemeVal)
  | Bytevector BS.ByteString
  | String Text
  | Character Char
  | Symbol Text
  | Boolean Bool
  | Number Number
  | Nil
  | Primitive ([SchemeVal] -> SchemeResult SchemeVal)
  | Fun Fn
  | IOFun ([SchemeVal] -> IOSchemeResult SchemeVal)
  | Port Handle
  deriving (Show, Typeable)

instance Eq SchemeVal where
  (==) (List x) (List y) = length x == length y && all (uncurry (==)) (zip x y)
  (==) (PairList xs x) (PairList ys y) = List (xs ++ [x]) == List (ys ++ [y])
  (==) (Vector x) (Vector y) = x == y
  (==) (Bytevector x) (Bytevector y) = x == y
  (==) (String x) (String y) = x == y
  (==) (Character x) (Character y) = x == y
  (==) (Symbol x) (Symbol y) = x == y
  (==) (Boolean x) (Boolean y) = x == y
  (==) (Number x) (Number y) = x == y
  (==) _ _ = False

unwordVals :: [SchemeVal] -> Text
unwordVals xs = T.unwords $ showVal <$> xs

showVal :: SchemeVal -> Text
showVal (List (Symbol "quote" : xs)) = "'" <> unwordVals xs
showVal (List (Symbol "quasiquote" : xs)) = "`" <> unwordVals xs
showVal (List (Symbol "unquote" : xs)) = "," <> unwordVals xs
showVal (List (Symbol "unquote-splicing" : xs)) = ",@" <> unwordVals xs
showVal (List contents) = "(" <> unwordVals contents <> ")"
showVal (PairList contents cdr) = "(" <> unwordVals contents <> " . " <> showVal cdr <> ")"
showVal (Vector vec) = T.pack $ "#(" <> unwords (map show $ elems vec) <> ")"
showVal (Bytevector vec) = T.pack $ "#u8(" <> unwords (map show $ BS.unpack vec) <> ")"
showVal (String s) = "\"" <> s <> "\""
showVal (Character a) = T.pack $ "#\\" <> [a]
showVal (Symbol s) = s
showVal (Boolean True) = "#t"
showVal (Boolean False) = "#f"
showVal (Number n) = showNum n
showVal Nil = "nil"
showVal Primitive {} = "<prim>"
showVal Fun {} = "<fun>"
showVal IOFun {} = "<io>"
showVal (Port _) = "<port>"

dumpAST :: SchemeVal -> IO ()
dumpAST = dumpAST' True

dumpAST' :: Bool -> SchemeVal -> IO ()
dumpAST' True = pPrint
dumpAST' False = pPrintLightBg

data SchemeError
  = Generic Text
  | ArgumentLengthMismatch Int [SchemeVal]
  | TypeMismatch Text SchemeVal
  | UnboundSymbol Text
  | ParserError Text
  | NotFunction SchemeVal
  | InvalidOperation Text
  | ReservedName Text
  | EmptyList
  deriving (Show)

instance Exception SchemeError

showError :: SchemeError -> Text
showError (TypeMismatch err vap) = "Invalid type: expected " <> err <> ", but found " <> T.pack (show vap)
showError (Generic err) = "Unexpectec error: " <> err
showError (UnboundSymbol sym) = "Unbound symbol: " <> sym
showError (ArgumentLengthMismatch ex act) = "Expected " <> T.pack (show ex) <> " but found " <> T.pack (show $ length act)
showError (ParserError err) = "Parsing error, could not parse input: " <> err
showError (NotFunction err) = "Attempt at calling " <> showVal err <> " as a function"
showError (InvalidOperation err) = "Invalid: " <> err
showError (ReservedName err) = "Cannot define : " <> err <> " because it is a reserved name"
showError EmptyList = "Empty list is not allowed"


text:
{-# OPTIONS_GHC -fno-cse #-}
{-# OPTIONS_GHC -fno-full-laziness #-}
{-# LANGUAGE ExistentialQuantification #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE BangPatterns #-}
module Trace(
      clearTrace
    , traceNames
    , CanBeTraced
    , traceValues
    , trace
    , traceSample
    , forceSignal
    ) where

import Common(HasDoubleRepresentation(..))
import System.IO.Unsafe
import qualified Data.Map as M
import Data.IORef
import Control.Applicative((<$>))
import Signal(Signal,mapS,fromListS,takeS)
import Fixed(Resolution(..))
import Statistics.Sample.KernelDensity
import qualified Data.Vector.Unboxed as U
import Viewer
import Plot
import Graphics.PDF
import Displayable
import Internal
import Control.DeepSeq
import Text.Printf
import qualified Graphics.PDF as PDF

import qualified Debug.Trace as T
debug a = T.trace (show a) a

data CanBeTraced = forall a. (HasDoubleRepresentation a,Resolution a) => CBT !a

instance HasDoubleRepresentation CanBeTraced where
    toDouble (CBT a) = toDouble a
    fromDouble = error "Can't create a CanBeTraced from a double"

instance Resolution CanBeTraced where
    smallestValue (CBT a) = CBT (smallestValue a)
    maxValue (CBT a) = CBT (maxValue a)
    minValue (CBT a) = CBT (minValue a)
    signedFormat (CBT a) = signedFormat a
    bitWidth (CBT a) = bitWidth a

{-# NOINLINE globalTrace #-}
globalTrace  :: IORef (M.Map String [CanBeTraced])
globalTrace = unsafePerformIO $ newIORef (M.empty)

clearTrace :: IO ()
clearTrace = writeIORef globalTrace M.empty

traceNames :: IO [String]
traceNames = M.keys <$> (readIORef  globalTrace)

getTrace :: String -> IO (Maybe [CanBeTraced])
getTrace name = M.lookup name <$> (readIORef globalTrace)

removeZeroAndLog :: [Double] -> [Double] -> ([Double],[Double])
removeZeroAndLog a b = unzip . map (\(t,x) -> (t, log x)) . filter ((> 0) . snd ) . zip a $ b

traceValues :: String -> IO ()
traceValues s = do
    let statStyle = LabelStyle 10 Centered PDF.S
    r <- getTrace s
    case r of
        Nothing -> return ()
        Just l -> do
            let h = head l
                sm = toDouble (smallestValue h)
                mi = toDouble (minValue h)
                ma = toDouble (maxValue h)
                triplePartition ta tb !la !lb !lc [] = (la,lb,lc)
                triplePartition ta tb !la !lb !lc (a:h) | ta a = triplePartition ta tb (a:la) lb lc h
                                                        | tb a = triplePartition ta tb la (a:lb) lc h
                                                        | otherwise = triplePartition ta tb la lb (a:lc) h
                (negatives, nulls, positives) = triplePartition (<0) (==0) [] [] [] (map toDouble l)
                positiveSat | null positives = 0 :: Double
                            | otherwise = 100*fromIntegral (length (filter (== ma) positives)) / fromIntegral (length positives)
                negativeSat | null negatives = 0 :: Double
                            | otherwise = 100*fromIntegral (length (filter (== mi) negatives)) / fromIntegral (length negatives)
                nullValues | null nulls = 0 :: Double
                           | otherwise = 100*fromIntegral (length nulls) / fromIntegral (length l)
                drawHist bmi bma l h label theTitle = do
                            let targetLine wi hi (ptF,_) = do
                                    withNewContext $ do
                                        strokeColor $ Rgb 1.0 0.8 0.8
                                        let start = (log bma) - (log bmi)
                                            (xa :+ ya) = ptF ((log bmi),(-start))
                                            (xb :+ yb) = ptF ((log bma),0)
                                        stroke $ Line  xa ya xb yb
                                tickv mav mbv =
                                    map (\t -> (fromIntegral t)/(fromIntegral 10)*(mbv-mav) + mav) ([0..10] :: [Int])
                                tickh _ _ =
                                    let mav = log bmi
                                        mbv = log bma
                                    in
                                    map (\t -> (fromIntegral t)/(fromIntegral 10)*(mbv-mav) + mav) ([0..10] :: [Int])
                                plotStyle =  (defaultPlotStyle { title = theTitle
                                                               , verticalLabel = Just "Log density"
                                                               , horizontalLabel = Just $ "Log amplitude (" ++ label ++ ")"
                                                               , horizontalTickValues = tickh
                                                               , verticalTickValues = tickv
                                                               , epilog = targetLine
                                                               , axis = False
                                                               , horizontalBounds = Just (log bmi, log bma)
                                                               })
                            if null l
                                then do
                                    (\_ -> \_ -> return (R Nothing Nothing), \_ -> \_ -> \_ -> return ())
                                else do
                                    let (t,p) = kde_ 32 (log bmi) (log bma) (U.fromList (map log l))
                                        lt' = U.toList t
                                        lp' = let templ = U.toList p
                                                  s = sum templ
                                              in
                                              map (/ s) templ
                                        (lt,lp) = removeZeroAndLog lt' lp'
                                    let r = (log bma - log bmi) / fromIntegral (length lp)
                                    drawing $ discreteSignalsWithStyle (length lt) plotStyle [ AS (fromListS 0 lp)]
            let (wi,he) = (\(x,y) -> (fromIntegral x, fromIntegral y)) $ viewerSize
            if (signedFormat h)
                then do
                    let posD = drawHist sm ma positives (he / 2 - 10) "pos" (Just "Log density of log amplitude")
                        negD = drawHist sm (-mi) (map negate negatives) (he / 2 - 10) "neg" Nothing
                        pict = Pair posD negD $ \pos -> \neg -> do
                            withNewContext $ do
                                applyMatrix $ translate (0 :+ (he/2 + 10))
                                pos
                            withNewContext $ do
                                neg
                            drawStringLabel statStyle ((printf "negative max: %2.3f %%" negativeSat) :: String)
                                             (leftMargin defaultPlotStyle + (wi / 4.0))  (he / 2 - 10) 100 20
                            drawStringLabel statStyle ((printf "nulls: %2.3f %%" nullValues) :: String)
                                             (leftMargin defaultPlotStyle + (2.0*wi/4.0)) (he/2 - 10) 100 20
                            drawStringLabel statStyle ((printf "positive max: %2.3f %%" positiveSat) :: String)
                                             (leftMargin defaultPlotStyle + (3.0*wi/4.0)) (he/2 - 10) 100 20
                    display pict
                else do
                        let dD = drawHist sm ma positives (he - 40) "unsigned" (Just "Log density of log amplitude")
                            pict = Mono dD $ \d -> do
                                withNewContext $ do
                                    applyMatrix $ translate (0 :+ 40)
                                    d
                                drawStringLabel statStyle ((printf "nulls: %2.3f %%" nullValues) :: String)
                                                    (leftMargin defaultPlotStyle + wi / 3.0) 10 100 20
                                drawStringLabel statStyle ((printf "positive max: %2.3f %%" positiveSat) :: String)
                                                    (leftMargin defaultPlotStyle + (2*wi/3.0)) 10 100 20
                        display pict

{-# NOINLINE traceSample #-}
traceSample :: (HasDoubleRepresentation a, Resolution a) => String -> a -> a
traceSample s a = unsafePerformIO $ do
    nv <- M.insertWith' (++) s [CBT a] <$> (readIORef globalTrace)
    writeIORef globalTrace nv
    return a

{-# NOINLINE trace #-}
trace :: (HasDoubleRepresentation a, Resolution a)
      => String
      -> Signal a
      -> Signal a
trace s = mapS (traceSample s)

{-# NOINLINE myTake #-}
myTake :: NFData a => Int -> [a] -> [a]
myTake 0 l = []
myTake i (h:l) = h:myTake (i-1) l
myTake _ [] = []

{-# NOINLINE forceSignal #-}
forceSignal :: (NFData a, Show a)
            => Int
            -> Signal a
            -> IO ()
forceSignal i (Signal s) = do
    let l = myTake i s
    l `deepseq` return ()

max_stars_repo_name: alpheccar/Signal

max_issues_repo_name: cpehle/Signal

max_forks_repo_name: cpehle/Signal


text:
{-# LANGUAGE ScopedTypeVariables, FlexibleContexts, BangPatterns #-}
module Statistics.Sampler.Slice (
  slice,
  newSlicerState,
  adaptOff,
  SlicerState()
) where

import           Random.CRI
import qualified Statistics.Distribution.Random.Exponential as D
import qualified Statistics.Distribution.Random.Uniform as D
import           Numeric.MathFunctions.Constants (m_epsilon)
import           Control.Monad (when)

import           Prelude hiding (max)

-- | slice sampling as described in
-- Radford M. Neal, "Slice sampling",
-- Ann. Statist. Volume 31, Number 3 (2003), 705-767.
-- http://projecteuclid.org/euclid.aos/1056562461

data SlicerState  = SlicerState {
    lower    :: !Double,             -- ^ lower bound of distribution
    upper    :: !Double,             -- ^ upper bound of distribution
    width    :: !Double,             -- ^ width of step out size (approximate scale parameter)
    steps    :: !Int,                -- ^ maximum number of step outs
    adapt    :: !Bool,               -- ^ adapt phase underway
    sumdiff  :: !Double,             -- ^ store sumdiff for adaption phase
    iter     :: !Int                 -- ^ number of iterations
    } deriving Show

newSlicerState :: Double    -- ^ lower bound
               -> Double    -- ^ upper bound
               -> Double    -- ^ initial width
               -> SlicerState
newSlicerState l u w = SlicerState {
                          lower   = l,
                          upper   = u,
                          width   = w,
                          steps   = 10,
                          adapt   = True,
                          iter    = 1,
                          sumdiff = 0.0
                        }

adaptOff :: SlicerState -> SlicerState
adaptOff st = st { adapt = False }

slice :: (Source m g Double) =>
         SlicerState
      -> (Double -> Double)       -- ^ x -> log(f x) where f is proportional to probibility density
      -> Double                   -- ^ current value
      -> g m                      -- ^ a random number generator
      -> m (SlicerState, Double)  -- ^ return slicer state and new sample value
slice st g x0 rng =
  do
    let !g0 = g x0

    when (isInfinite g0) $
        error $ "Infinite value found in slice sampler: " ++ show x0 ++ " -> " ++ show g0

    when (isNaN g0) $
        error $ "NaN found in slice sampler: " ++ show x0 ++ " -> " ++ show g0

    -- 1. define slice
    e <- D.exponential rng
    let !z = g0 - e

    -- 2. find interval
    u <- D.uniform rng
    let !l = x0 - width st * u
        !r = l + width st

    v :: Double <- D.uniform rng
    let !j = floor (fromIntegral (steps st) * v)
        !k = (steps st - 1) - j

    let !left = calc_left j l
        calc_left n l'
            | l' < lower st  = lower st
            | n == 0         = l'
            | z >= g l'      = l'
            | otherwise      = calc_left  (n-1) (l' - width st)

    let !right = calc_right k r
        calc_right n r'
            | r' > upper st  = upper st
            | n == 0         = r'
            | z >= g r'      = r'
            | otherwise      = calc_right (n-1) (r' + width st)

    -- 3. loop until accept (guaranteed)
    let sample left' right' =
          do
            u' <- D.uniform rng
            let !x = left' + u' * (right' - left')
            if z - m_epsilon <= g x
              then return x -- accept
              else
                if x < x0
                  then sample x     right'
                  else sample left' x

    !x1 <- sample left right

    let !st' | adapt st  = adaptSlicer x0 x1 st
             | otherwise = st

    return (st', x1)

adaptSlicer :: Double -> Double -> SlicerState -> SlicerState
adaptSlicer old new oldst = newst
  where
    !iterf    = fromIntegral (iter oldst)
    !sumdiff' = sumdiff oldst +  iterf * abs (new - old)
    !newst    = oldst {
                    sumdiff = sumdiff',
                    iter    = iter oldst + 1,
                    width   = if iter oldst > 50
                                 then 2 * sumdiff' / (iterf * (iterf - 1))
                                 else width oldst
                    }

max_stars_repo_name: finlay/random-dist

max_issues_repo_name: finlay/random-dist

max_forks_repo_name: finlay/random-dist


text:
{-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}
-- |
-- Module    : Statistics.Distribution.Hypergeometric
-- Copyright : (c) 2009 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- The Hypergeometric distribution.  This is the discrete probability
-- distribution that measures the probability of /k/ successes in /l/
-- trials, without replacement, from a finite population.
--
-- The parameters of the distribution describe /k/ elements chosen
-- from a population of /l/, with /m/ elements of one type, and
-- /l/-/m/ of the other (all are positive integers).

module Statistics.Distribution.Hypergeometric
    (
      HypergeometricDistribution
    -- * Constructors
    , hypergeometric
    -- ** Accessors
    , hdM
    , hdL
    , hdK
    ) where

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

data HypergeometricDistribution = HD {
      hdM :: {-# UNPACK #-} !Int
    , hdL :: {-# UNPACK #-} !Int
    , hdK :: {-# UNPACK #-} !Int
    } deriving (Eq, Read, Show, Typeable, Data, Generic)

instance FromJSON HypergeometricDistribution
instance ToJSON HypergeometricDistribution

instance Binary HypergeometricDistribution where
    get = HD <$> get <*> get <*> get
    put (HD x y z) = put x >> put y >> put z

instance D.Distribution HypergeometricDistribution where
    cumulative = cumulative

instance D.DiscreteDistr HypergeometricDistribution where
    probability = probability

instance D.Mean HypergeometricDistribution where
    mean = mean

instance D.Variance HypergeometricDistribution where
    variance = variance

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

instance D.MaybeVariance HypergeometricDistribution where
    maybeStdDev   = Just . D.stdDev
    maybeVariance = Just . D.variance

instance D.Entropy HypergeometricDistribution where
  entropy = directEntropy

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

variance :: HypergeometricDistribution -> Double
variance (HD m l k) = (k' * ml) * (1 - ml) * (l' - k') / (l' - 1)
  where m' = fromIntegral m
        l' = fromIntegral l
        k' = fromIntegral k
        ml = m' / l'

mean :: HypergeometricDistribution -> Double
mean (HD m l k) = fromIntegral k * fromIntegral m / fromIntegral l

directEntropy :: HypergeometricDistribution -> Double
directEntropy d@(HD m _ _) =
    negate . sum $
  takeWhile (< negate m_epsilon) $
  dropWhile (not . (< negate m_epsilon)) $
  [ let x = probability d n in x * log x | n <- [0..m]]

hypergeometric :: Int               -- ^ /m/
               -> Int               -- ^ /l/
               -> Int               -- ^ /k/
               -> HypergeometricDistribution
hypergeometric m l k
  | not (l > 0)            = error $ msg ++ "l must be positive"
  | not (m >= 0 && m <= l) = error $ msg ++ "m must lie in [0,l] range"
  | not (k > 0 && k <= l)  = error $ msg ++ "k must lie in (0,l] range"
  | otherwise = HD m l k
    where
      msg = "Statistics.Distribution.Hypergeometric.hypergeometric: "

-- Naive implementation
probability :: HypergeometricDistribution -> Int -> Double
probability (HD mi li ki) n
  | n < max 0 (mi+ki-li) || n > min mi ki = 0
  | otherwise =
      choose mi n * choose (li - mi) (ki - n) / choose li ki

cumulative :: HypergeometricDistribution -> Double -> Double
cumulative d@(HD mi li ki) x
  | isNaN x      = error "Statistics.Distribution.Hypergeometric.cumulative: NaN argument"
  | isInfinite x = if x > 0 then 1 else 0
  | n <  minN    = 0
  | n >= maxN    = 1
  | otherwise    = D.sumProbabilities d minN n
  where
    n    = floor x
    minN = max 0 (mi+ki-li)
    maxN = min mi ki

max_stars_repo_name: StefanHubner/statistics

max_issues_repo_name: StefanHubner/statistics

max_forks_repo_name: StefanHubner/statistics


text:
{-|
Module      : MachineLearning.Model.Measure
Description : TIR expression data structures
Copyright   : (c) Fabricio Olivetti de Franca, 2022
License     : GPL-3
Maintainer  : fabricio.olivetti@gmail.com
Stability   : experimental
Portability : POSIX

Performance measures for Regression and Classification.
-}
module MachineLearning.Model.Measure
  ( Measure(..)
  , toMeasure
  , measureAll
  , _rmse
  )
  where

import Data.Semigroup (Sum(..))

import qualified Numeric.LinearAlgebra       as LA
import qualified Numeric.Morpheus.Statistics as Stat
import qualified Data.Vector.Storable        as V

type Vector = LA.Vector Double
-- * Performance measures

-- | A performance measure has a string name and a function that
-- takes a vector of the true values, a vector of predict values
-- and returns a `Double`.
data Measure = Measure { _name :: String
                       , _fun  :: Vector -> Vector -> Double -- ^ true values -> predicted values -> measure
                       }

instance Show Measure where
  show (Measure n _) = n

-- | Mean for a vector of doubles
mean :: Vector -> Double
mean xs = V.sum xs / fromIntegral (V.length xs)
{-# INLINE mean #-}

-- | Variance for a vector of doubles
var :: Vector -> Double
var xs = sum' / fromIntegral (V.length xs)
  where
    mu   = mean xs
    sum' = V.foldl (\s x -> s + (x-mu)^(2 :: Int)) 0 xs
{-# INLINE var #-}

-- | generic mean error measure
meanError :: (Vector -> Vector) -- ^ a function to be applied to the error terms (abs, square,...)
          -> Vector             -- ^ target values
          -> Vector             -- ^ fitted values
          -> Double
meanError op ys ysHat = mean $ op $ ysHat - ys
{-# INLINE meanError #-}

-- * Common error measures for regression:
-- MSE, MAE, RMSE, NMSE, r^2

-- | Mean Squared Error
mse :: Vector -> Vector -> Double
--mse           = meanError (^(2 :: Int))
mse ys ysHat = mean $ (ysHat - ys) ^(2 :: Int)
{-# INLINE mse #-}

-- | Mean Absolute Error
mae :: Vector -> Vector -> Double
mae ys ysHat = mean $ abs (ysHat - ys) -- meanError abs
{-# INLINE mae #-}

-- | Normalized Mean Squared Error
nmse :: Vector -> Vector -> Double
nmse ys ysHat = mse ysHat ys / var ys
{-# INLINE nmse #-}

-- | Root of the Mean Squared Error
rmse :: Vector -> Vector -> Double
rmse ys ysHat = sqrt $ mse ysHat ys
{-# INLINE rmse #-}

-- | negate R^2 - minimization metric
rSq :: Vector -> Vector -> Double
rSq ys ysHat = negate (1 - r/t)
  where
    ym      = Stat.mean ys
    t       = sumOfSq $ V.map (\yi -> yi - ym) ys
    r       = sumOfSq $ ys - ysHat
    sumOfSq = V.foldl (\s di -> s + di^(2 :: Int)) 0
{-# INLINE rSq #-}

-- * Regression measures
_rmse, _mae, _nmse, _r2 :: Measure
_rmse = Measure "RMSE" rmse
_mae  = Measure "MAE" mae
_nmse = Measure "NMSE" nmse
_r2   = Measure "R^2" rSq

-- * Classification measures
_accuracy,_recall,_precision,_f1,_logloss :: Measure
_accuracy  = Measure "Accuracy" accuracy
_recall    = Measure "Recall" recall
_precision = Measure "Precision" precision
_f1        = Measure "F1" f1
_logloss   = Measure "Log-Loss" logloss

-- | Accuracy: ratio of correct classification
accuracy :: Vector -> Vector -> Double
accuracy ys ysHat = -equals/tot
  where
    ys'    = map round $ LA.toList ys
    ysHat' = map round $ LA.toList ysHat
    (Sum equals, Sum tot) = foldMap cmp $ zip ysHat' ys'
    cmp :: (Integer, Integer) -> (Sum Double, Sum Double)
    cmp (yH, y)
      | yH == y   = (Sum 1, Sum 1)
      | otherwise = (Sum 0, Sum 1)

-- | Precision: ratio of correct positive classification
precision :: Vector -> Vector -> Double
precision ys ysHat = equals/tot
  where
    ys'    = map round $ LA.toList ys
    ysHat' = map round $ LA.toList ysHat
    (Sum equals, Sum tot) = foldMap cmp $ zip ysHat' ys'
    cmp :: (Integer, Integer) -> (Sum Double, Sum Double)
    cmp (1, 1)  = (Sum 1, Sum 1)
    cmp (1, 0)  = (Sum 0, Sum 1)
    cmp (_, _) = (Sum 0, Sum 0)

-- | Recall: ratio of retrieval of positive labels
recall :: Vector -> Vector -> Double
recall ys ysHat = equals/tot
  where
    ys'    = map round $ LA.toList ys
    ysHat' = map round $ LA.toList ysHat
    (Sum equals, Sum tot) = foldMap cmp $ zip ysHat' ys'

    cmp :: (Integer, Integer) -> (Sum Double, Sum Double)
    cmp (1, 1)  = (Sum 1, Sum 1)
    cmp (0, 1)  = (Sum 0, Sum 1)
    cmp (_, _) = (Sum 0, Sum 0)

-- | Harmonic average between Precision and Recall
f1 :: Vector -> Vector -> Double
f1 ys ysHat = 2*prec*rec/(prec+rec)
  where
    prec = precision ysHat ys
    rec  = recall ysHat ys

-- | LogLoss of a classifier that returns a probability.
logloss :: Vector -> Vector -> Double
logloss ys ysHat = mean $ -(ys * log ysHat' + (1 - ys)*log(1 - ysHat'))
  where
    ysHat' = LA.cmap (min (1.0 - 1e-15) . max 1e-15) ysHat

-- | List of all measures
measureAll :: [Measure]
measureAll = [_rmse, _mae, _nmse, _r2
             , _accuracy, _recall, _precision, _f1, _logloss
             ]

-- | Read a string into a measure
toMeasure :: String -> Measure
toMeasure input
  | null cmp  = error ("Invalid measure: " ++ input)
  | otherwise = (snd.head) cmp
  where
    cmp                       = filter fst $ map isThis measureAll
    isThis m@(Measure name _) = (name == input, m)

max_stars_repo_name: folivetti/tir

max_issues_repo_name: folivetti/tir

max_forks_repo_name: folivetti/tir


text:
{-# LANGUAGE ConstraintKinds       #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE KindSignatures        #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE PatternSynonyms       #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE TupleSections         #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeInType            #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE UndecidableInstances  #-}
{-# LANGUAGE ViewPatterns          #-}

module Backprop.Learn.Model.Types (
  -- * Model type
    ModelFunc, ModelFuncStoch
  , Model(..)
  , modelD
  -- * Specialized Models
  -- ** Stateless
  , ModelFuncStateless, ModelFuncStochStateless
  , ModelStateless, pattern ModelStateless, runLearnStateless, runLearnStochStateless
  , modelStatelessD
  -- ** Stateless and Parameterless
  , BFunc, BFuncStoch
  , Func, pattern Func, runFunc, runFuncStoch
  , funcD
  -- * Manipulating models as functions
  , ModelFuncM, withModelFunc0, withModelFunc, withModelFunc2
  -- * Utility
  , PMaybe(..), TMaybe, fromPJust, Learnables, Learnable
  , Regularize(..), Initialize(..)
  -- , HKD
  ) where

import           Backprop.Learn.Initialize
import           Backprop.Learn.Regularize
import           Control.Category
import           Control.DeepSeq
import           Control.Monad.Primitive
import           Control.Monad.Trans.Reader
import           Data.Functor.Identity
import           Data.Kind
import           Data.Type.Functor.Product
import           Data.Type.Tuple
import           Data.Vinyl
import           GHC.TypeNats
import           Numeric.Backprop
import           Numeric.Opto.Update
import           Prelude hiding               ((.))
import qualified Data.Binary                  as Bi
import qualified Numeric.LinearAlgebra.Static as H
import qualified System.Random.MWC            as MWC

type ModelFunc p s a b = forall z. Reifies z W
    => PMaybe (BVar z) p
    -> BVar z a
    -> PMaybe (BVar z) s
    -> (BVar z b, PMaybe (BVar z) s)

type ModelFuncStoch p s a b = forall m z. (PrimMonad m, Reifies z W)
    => MWC.Gen (PrimState m)
    -> PMaybe (BVar z) p
    -> BVar z a
    -> PMaybe (BVar z) s
    -> m (BVar z b, PMaybe (BVar z) s)

-- | General parameterized model with potential state
data Model :: Maybe Type -> Maybe Type -> Type -> Type -> Type where
    Model :: { runLearn      :: ModelFunc      p s a b
             , runLearnStoch :: ModelFuncStoch p s a b
             } -> Model p s a b

-- | Construct a deterministic model, with no stochastic component.
modelD :: ModelFunc p s a b -> Model p s a b
modelD f = Model { runLearn      = f
                 , runLearnStoch = \_ p x -> pure . f p x
                 }

type ModelFuncStateless p a b = forall z. Reifies z W
    => PMaybe (BVar z) p
    -> BVar z a
    -> BVar z b

type ModelFuncStochStateless p a b = forall m z. (PrimMonad m, Reifies z W)
    => MWC.Gen (PrimState m)
    -> PMaybe (BVar z) p
    -> BVar z a
    -> m (BVar z b)

newtype MFS  p a b = MFS  (ModelFuncStateless      p a b)
newtype MFSS p a b = MFSS (ModelFuncStochStateless p a b)

-- | Parameterized model with no state
type ModelStateless p = Model p 'Nothing

-- | Construct a 'ModelStateless'.
pattern ModelStateless
    :: ModelFuncStateless p a b
    -> ModelFuncStochStateless p a b
    -> ModelStateless p a b
pattern ModelStateless { runLearnStateless, runLearnStochStateless } <-
    (\m -> (_runLearnStateless m, _runLearnStochStateless m)
       -> (MFS runLearnStateless, MFSS runLearnStochStateless)
    )
  where
    ModelStateless l ls = Model
      { runLearn      = \  p x s -> (l p x, s)
      , runLearnStoch = \g p x s -> (, s) <$> ls g p x
      }

-- | Construct a deterministic stateless model, with no stochastic
-- component.
modelStatelessD :: ModelFuncStateless p a b -> ModelStateless p a b
modelStatelessD f = ModelStateless
    { runLearnStateless      = f
    , runLearnStochStateless = \_ p -> pure . f p
    }

_runLearnStateless :: Model p 'Nothing a b -> MFS p a b
_runLearnStateless md = MFS $ \p x -> fst $ runLearn md p x PNothing

_runLearnStochStateless :: Model p 'Nothing a b -> MFSS p a b
_runLearnStochStateless md = MFSS $ \g p x -> fst <$> runLearnStoch md g p x PNothing

type BFunc a b = forall z. Reifies z W
    => BVar z a
    -> BVar z b
type BFuncStoch a b = forall m z. (PrimMonad m, Reifies z W)
    => MWC.Gen (PrimState m)
    -> BVar z a
    -> m (BVar z b)

newtype BF  a b = BF  (BFunc      a b)
newtype BFS a b = BFS (BFuncStoch a b)

type Func = Model 'Nothing 'Nothing

-- | Unparameterized model with no state
pattern Func :: BFunc a b -> BFuncStoch a b -> Func a b
pattern Func { runFunc, runFuncStoch } <-
    (\m -> (_runFunc m, _runFuncStoch m)
       -> (BF runFunc, BFS runFuncStoch)
    )
  where
    Func r rs = ModelStateless
      { runLearnStateless      = const r
      , runLearnStochStateless = const . rs
      }

-- | Construct an deterministic unparameterized stateless model, with no
-- stochastic component.
funcD :: BFunc a b -> Func a b
funcD f = Func { runFunc      = f
               , runFuncStoch = const (pure . f)
               }

_runFunc :: Model 'Nothing 'Nothing a b -> BF a b
_runFunc md = BF $ runLearnStateless md PNothing

_runFuncStoch :: Model 'Nothing 'Nothing a b -> BFS a b
_runFuncStoch md = BFS $ flip (runLearnStochStateless md) PNothing

-- | Share parameter and sequence state
instance Category (Model p s) where
    id    = modelD $ \_ x s -> (x, s)
    f . g = Model
      { runLearn      = \p x s0 -> let (y, s1) = runLearn g p x s0
                                   in  runLearn f p y s1
      , runLearnStoch = \gen p x s0 -> do
          (y, s1) <- runLearnStoch g gen p x s0
          runLearnStoch f gen p y s1
      }

type ModelFuncM m p s a b = forall z. Reifies z W
    => PMaybe (BVar z) p
    -> BVar z a
    -> PMaybe (BVar z) s
    -> m (BVar z b, PMaybe (BVar z) s)

withModelFunc0
    :: (forall m. Monad m => ModelFuncM m p s a b)
    -> Model p s a b
withModelFunc0 f = Model
    { runLearn      = \p x -> runIdentity . f p x
    , runLearnStoch = \g p x -> flip runReaderT g . f p x
    }

withModelFunc
    :: (forall m. Monad m => ModelFuncM m p s a b -> ModelFuncM m q t c d)
    -> Model p s a b
    -> Model q t c d
withModelFunc f m = Model
    { runLearn      = \p x -> runIdentity . f (runLearnModelFuncM m) p x
    , runLearnStoch = \g p x -> flip runReaderT g . f (runLearnStochModelFuncM m) p x
    }

withModelFunc2
    :: (forall m. Monad m => ModelFuncM m p s a b -> ModelFuncM m q t c d -> ModelFuncM m r u e f)
    -> Model p s a b
    -> Model q t c d
    -> Model r u e f
withModelFunc2 f m n = Model
    { runLearn      = \p x ->
            runIdentity
          . f (runLearnModelFuncM m) (runLearnModelFuncM n) p x
    , runLearnStoch = \g p x ->
            flip runReaderT g
          . f (runLearnStochModelFuncM m) (runLearnStochModelFuncM n) p x
    }

runLearnModelFuncM
    :: Model p s a b
    -> ModelFuncM Identity p s a b
runLearnModelFuncM m p x = Identity . runLearn m p x

runLearnStochModelFuncM
    :: PrimMonad m
    => Model p s a b
    -> ModelFuncM (ReaderT (MWC.Gen (PrimState m)) m) p s a b
runLearnStochModelFuncM m p x s = ReaderT $ \g -> runLearnStoch m g p x s

-- | Combination of common constraints for type-level lists.
type Learnables as = ( RecApplicative as
                     , ReifyConstraint Backprop TF as
                     , RMap as
                     , RApply as
                     )

-- | Helpful utility class for ensuring that a parameter or state type is
-- Learnable
class ( Initialize a, Regularize a, Bi.Binary a, NFData a, LinearInPlace m Double a )
        => Learnable m a

instance (Initialize a, Regularize a, Bi.Binary a, NFData a, LinearInPlace m Double a)
    => Learnable m (TF a)
instance (Initialize a, Initialize b, Regularize a, Regularize b, Bi.Binary a, Bi.Binary b, NFData a, NFData b, LinearInPlace m Double a, LinearInPlace m Double b)
    => Learnable m (a :# b)

instance (RPureConstrained Initialize as, RPureConstrained Regularize as, RPureConstrained Bi.Binary as, RPureConstrained NFData as, ReifyConstraint Backprop TF as, RMap as, RApply as, RFoldMap as, RecordToList as, LinearInPlace m Double (Rec TF as))
        => Learnable m (Rec TF as)
instance (KnownNat n, PrimMonad m) => Learnable m (H.R n)
instance (KnownNat n, KnownNat o, PrimMonad m) => Learnable m (H.L n o)

max_stars_repo_name: mstksg/backprop-learn

max_issues_repo_name: mstksg/backprop-learn

max_forks_repo_name: mstksg/backprop-learn


text:
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE DeriveGeneric #-}
{-# OPTIONS_GHC -funbox-strict-fields #-}

module Warden.Anomaly.Data (
    FeatureVector(..)
  , Features(..)
  , featuresFromList
  , Dimensionality(..)
  , dimensionality
  , fromDataRow
  , toMatrix
  ) where

import           Data.ByteString (ByteString)
import           Data.List.NonEmpty (NonEmpty(..))
import qualified Data.List.NonEmpty as NE
import qualified Data.Vector as V
import qualified Data.Vector.Unboxed as VU

import           Numeric.LinearAlgebra (R)
import qualified Numeric.LinearAlgebra as HM

import           P

import           Warden.Data.Numeric
import           Warden.Data.Row
import           Warden.Row (parseProbablyNumeric)

newtype FeatureVector =
  FeatureVector {
    unFeatureVector :: VU.Vector Double
  } deriving (Eq, Show)

numFeatures :: FeatureVector -> Dimensionality
numFeatures = Dimensionality . VU.length . unFeatureVector

-- | Contains one or more feature vectors. All feature vectors have the same
-- dimensionality.
newtype Features =
  Features {
    unFeatures :: V.Vector FeatureVector
  } deriving (Eq, Show)

newtype Dimensionality =
  Dimensionality {
    unDimensionality :: Int
  } deriving (Eq, Show)

dimensionality :: Features -> Dimensionality
dimensionality (Features fvs) =
  numFeatures $ V.head fvs

featuresFromList :: NonEmpty FeatureVector -> Maybe' Features
featuresFromList (v:|vs) =
  let
    dim = numFeatures v
    bads = filter (\x -> numFeatures x /= dim) vs
  in
  case bads of
    [] ->
      Just' . Features . V.fromList . NE.toList $ v :| vs
    _ ->
      Nothing'

getNumericField :: V.Vector ByteString -> Int -> Maybe' Double
getNumericField fs ix = do
  f <- strictMaybe $ fs V.!? ix
  case parseProbablyNumeric f of
    NoNumericField ->
      Nothing'
    MNumericField (NumericField x) ->
      pure x

fromDataRow :: [Int] -> Row -> Maybe' FeatureVector
fromDataRow nixes (SVFields fs) = do
  nums <- mapM (getNumericField fs) nixes
  pure . FeatureVector $ VU.fromList nums
fromDataRow _ (RowFailure _) =
  Nothing'

toMatrix :: [FeatureVector] -> HM.Matrix R
toMatrix =
  HM.fromRows . fmap (HM.vector . VU.toList . unFeatureVector)

max_stars_repo_name: ambiata/warden

max_issues_repo_name: ambiata/warden

max_forks_repo_name: ambiata/warden


text:
{-# LANGUAGE NamedFieldPuns, TemplateHaskell #-}

module School.Train.Test.GradientDescent
( gradientDescentTest ) where

import Conduit ((.|), liftIO, runConduit, sinkNull, yield, yieldMany)
import Data.Either (isLeft)
import Numeric.LinearAlgebra ((><))
import School.TestUtils (addClasses, assertRight, empty, isSorted,
                         randomAffineParams, randomNNInts, randomMatrix,
                         unitCorrect, weight1)
import School.FileIO.AppIO (runAppIO)
import School.FileIO.FileType (FileType, toExtension)
import School.FileIO.FileHeader (FileHeader(..))
import School.FileIO.MatrixSink (matrixIntSink)
import School.Train.AppTrain (AppTrain)
import School.Train.GradientDescent
import School.Train.SimpleDescentUpdate (simpleDescentUpdate)
import School.Train.StoppingCondition (maxIterations)
import School.Train.IterationHandler (storeCost)
import School.Train.TrainState (TrainState(..), HandlerStore(..), def)
import School.Types.PingPong (pingPongSingleton)
import School.Types.DataType (DataType(..))
import School.Types.Error (Error)
import School.Types.FloatEq ((~=))
import School.Unit.CostFunction (CostFunction)
import School.Unit.Affine (affine)
import School.Unit.MultiNoulli (multiNoulli)
import School.Unit.RecLin (recLin)
import School.Unit.UnitParams (UnitParams(..))
import System.Directory (removeFile)
import Test.Tasty (TestTree)
import Test.Tasty.QuickCheck hiding ((><))
import Test.Tasty.TH
import Test.QuickCheck.Monadic (assert, monadicIO)

prop_no_units :: Property
prop_no_units = monadicIO $ do
  result <- liftIO $ gradientDescent (yield empty)
                                     []
                                     weight1
                                     simpleDescentUpdate
                                     mempty
                                     mempty
                                     def
  assert $ isLeft result

prop_iterations :: Positive Int -> Positive Int -> Positive Int -> Property
prop_iterations (Positive n) (Positive b) (Positive f) = monadicIO $ do
  input <- liftIO $ randomMatrix b f
  let source =  yieldMany . repeat $ input
  result <- liftIO $ gradientDescent source
                                    [recLin]
                                    weight1
                                    simpleDescentUpdate
                                    (maxIterations n)
                                    mempty
                                    def
  assertRight ((== n) . iterationCount) result

prop_cost_decline :: Positive Int -> Positive Int -> Positive Int -> Property
prop_cost_decline (Positive b) (Positive f) (Positive o) = monadicIO $ do
  input <- liftIO $ randomMatrix b f
  let source =  yieldMany . repeat $ input
  paramList <- liftIO $ pingPongSingleton <$> randomAffineParams f o
  let initState = def { handlerStore = CostList []
                                  , learningRate = 1e-2
                                  , paramList
                                  }
  result <- liftIO $ gradientDescent source
                                     [affine]
                                     weight1
                                     simpleDescentUpdate
                                     (maxIterations 5)
                                     storeCost
                                     initState
  assertRight ((\(CostList c) -> isSorted c) . handlerStore)
              result

multiSingle :: CostFunction Double (AppTrain a)
multiSingle = multiNoulli Nothing Nothing

prop_multinoulli_single_file :: Positive Int -> Positive Int -> Property
prop_multinoulli_single_file (Positive c) (Positive b) = monadicIO $ do
  activation <- liftIO $ randomMatrix b c
  classes <- liftIO $ randomNNInts (c - 1) b
  let input = addClasses classes activation
  let source =  yieldMany . repeat $ input
  let initState = def { handlerStore = CostList []
                      , paramList = pingPongSingleton EmptyParams
                      }
  result <- liftIO $ gradientDescent source
                                     [unitCorrect classes]
                                     multiSingle
                                     simpleDescentUpdate
                                     (maxIterations 5)
                                     storeCost
                                     initState
  let check = CostList $ replicate 5 (-1)
  assertRight ((~= check) . handlerStore) result

multiTwo :: [Int]
         -> Int
         -> FileType
         -> IO ( Either Error ( FilePath
                              , CostFunction Double (AppTrain a)))
multiTwo classes iterations fType = do
  let fName = "test." ++ (toExtension fType)
  let cols = 1
  let nClasses = length classes
  let rows = iterations * nClasses
  let header = FileHeader { dataType = INT32B, cols, rows }
  let matrix = (rows >< cols) . concat . (replicate iterations)
             . (map fromIntegral) $ classes
  writeRes <- runAppIO . runConduit $ yield matrix
                                   .| matrixIntSink fType header fName
                                   .| sinkNull
  let multi = multiNoulli (Just fName) (Just $ header { rows = nClasses })
  either (return . Left)
         (return . pure . Right $ (fName, multi))
         writeRes

prop_multinoulli_two_files :: Positive Int -> Positive Int -> FileType -> Property
prop_multinoulli_two_files (Positive c) (Positive b) fType = monadicIO $ do
  activation <- liftIO $ randomMatrix b c
  classes <- liftIO $ randomNNInts (c - 1) b
  let input = addClasses classes activation
  let source =  yieldMany . repeat $ input
  let initState = def { handlerStore = CostList []
                      , paramList = pingPongSingleton EmptyParams
                      }
  let iterations = 5
  multiResult <- liftIO $ multiTwo classes iterations fType
  either (\e -> do
            liftIO . putStrLn $ "ERROR " ++ e
            assert False)
         (\(fName, multi) -> do
            result <- liftIO $ gradientDescent source
                                               [unitCorrect classes]
                                               multi
                                               simpleDescentUpdate
                                               (maxIterations iterations)
                                               storeCost
                                               initState
            let check = CostList $ replicate iterations (-1)
            liftIO $ removeFile fName
            assertRight ((~= check) . handlerStore) result
          )
          multiResult

gradientDescentTest :: TestTree
gradientDescentTest = $(testGroupGenerator)

max_stars_repo_name: jfulseca/School

max_issues_repo_name: jfulseca/School

max_forks_repo_name: jfulseca/School


text:
{-# LANGUAGE Strict     #-}
{-# LANGUAGE StrictData #-}
module FokkerPlanck.Analytic
  ( R2S1RP(..)
  , computePji
  , computePjiCorner
  , computePjiCorner'
  ) where

import           Data.Complex
import           Data.List               as L
import           Numeric.GSL.Polynomials
import           Text.Printf

data R2S1RP = R2S1RP
  { xR2S1RP, yR2S1RP, thetaR2S1RP, gammaR2S1RP :: Double
  } deriving (Show)

instance Eq R2S1RP where
  (==) (R2S1RP x1 y1 theta1 _) (R2S1RP x2 y2 theta2 _) =
    (x1 == x2) && (y1 == y2)

{-# INLINE computeCoefficients #-}
computeCoefficients :: R2S1RP -> R2S1RP -> (Double, Double, Double)
computeCoefficients (R2S1RP x_i y_i theta_i _) (R2S1RP x_j y_j theta_j gamma) =
  let x = x_j - x_i
      y = y_j - y_i
      a = (2 + cos (theta_j - theta_i)) / 3
      b =
        (x * (cos theta_j + cos theta_i) + y * (sin theta_j + sin theta_i)) /
        gamma
      c = (x ^ 2 + y ^ 2) / (gamma ^ 2)
   in (a, b, c)

{-# INLINE solveCubicEquation #-}
solveCubicEquation :: Double -> Double -> Double -> Double -> [Complex Double]
solveCubicEquation a b c d =
  let x0 = b ^ 2 - 3 * a * c
      x1 = 2 * b ^ 3 - 9 * a * b * c + 27 * a ^ 2 * d
      x2 = x1 ^ 2 - 4 * x0 ^ 3
      y =
        if x2 >= 0
          then ((x1 + sqrt x2) / 2) ** (1 / 3) :+ 0
          else ((x1 :+ sqrt (abs x2)) / (2 :+ 0)) ** (1 / 3)
      z = (-0.5) :+ 0.8660254037844386 -- (sqrt 3) / 2
  in [ ((b :+ 0) + z ^ k * y + (x0 :+ 0) / (z ^ k * y)) / ((-3) * a :+ 0)
     | k <- [0 .. 2]
     ]

{-# INLINE isRealPositive #-}
isRealPositive :: Complex Double -> Bool
isRealPositive (a :+ b) = (a > 0) && (abs b < 1e-10)

{-# INLINE computePjit #-}
computePjit ::
     Double -> Double -> (Double, Double, Double) -> Double -> Double
computePjit sigma tau (a, b, c) t =
  3 * exp ((-6) * (a * (t ^ 2) - b * t + c) / (sigma * (t ^ 3))) *
  exp (-t / tau) /
  sqrt ((pi * sigma) ^ 3 * (t ^ 7) / 2)

computePji :: Double -> Double -> R2S1RP -> R2S1RP -> Double
computePji sigma tau x_i x_j =
  if x_i == x_j
    then 0
    else let coef@(a, b, c) = computeCoefficients x_i x_j
             roots' =
               polySolve [9 * c / sigma, -6 * b / sigma, 3 * a / sigma, -7 / 4]
             roots = L.filter isRealPositive roots'
             (pOpt, tOpt) =
               L.maximumBy (\x y -> compare (fst x) (fst y)) .
               L.map
                 (\root ->
                    let realRoot = realPart root
                     in (computePjit sigma tau coef realRoot, realRoot)) $
               roots
             f =
               sqrt
                 (2 * pi * (tOpt ^ 5) /
                  (12 * (3 * c - b * tOpt) / sigma + 7 * (tOpt ^ 3) / 2))
          in f * pOpt

{-# INLINE findIntersectionPoint #-}
findIntersectionPoint :: Double -> R2S1RP -> Maybe Double
findIntersectionPoint delta (R2S1RP x y theta _)
  | theta == pi || theta == 0 || z < delta = Nothing
  | otherwise = Just z
  where
    z = x - y / tan theta

{-# INLINE rotate #-}
rotate :: (Double,Double) -> Double -> (Double,Double)
rotate (x, y) theta =
  (x * cos theta - y * sin theta, x * sin theta + y * cos theta)

{-# INLINE computePjiCorner #-}
computePjiCorner :: Double -> Double -> Double -> [Double] -> R2S1RP -> Double
computePjiCorner delta sigma tau orientations point@(R2S1RP x y theta gamma) =
  case findIntersectionPoint delta point of
    Nothing -> 0
    Just z ->
      let p1 = computePji sigma tau (R2S1RP 0 0 0 gamma) (R2S1RP z 0 0 gamma)
                      -- (newX, newY) = rotate (x - z, y) (-theta)
                      -- p2 =
                      --   computePji
                      --     sigma
                      --     tau
                      --     (R2S1RP 0 0 0 gamma)
                      --     (R2S1RP newX newY 0 gamma)
          p2 = computePji sigma tau (R2S1RP z 0 theta gamma) point
                      -- p1 = L.foldl' (\s ori -> s + computePji sigma tau
                      --   (R2S1RP 0 0 0 gamma) (R2S1RP z 0 ori gamma)) 0
                      --   orientations p2 = L.foldl' (\s ori -> s +
                      --   computePji sigma tau (R2S1RP z 0 ori gamma)
                      --   point) 0 orientations
       in p1 * p2

{-# INLINE computePjiCorner' #-}
computePjiCorner' :: Double -> Double -> Double -> Double -> [Double] -> R2S1RP -> Double
computePjiCorner' delta sigma tau threshold orientations point@(R2S1RP x y theta gamma) =
  case findIntersectionPoint delta point of
    Nothing -> 0
    Just z ->
      let p1 = computePji sigma tau (R2S1RP 0 0 0 gamma) (R2S1RP z 0 0 gamma)
       in if p1 < threshold
            then 0
            else let p11 =
                       L.foldl'
                         (\s ori ->
                            s +
                            computePji
                              sigma
                              tau
                              (R2S1RP 0 0 0 gamma)
                              (R2S1RP z 0 ori gamma))
                         0
                         orientations
                                               -- p2 = computePji sigma tau (R2S1RP z 0 theta gamma) point
                     p2 =
                       L.foldl'
                         (\s ori ->
                            s +
                            computePji
                              sigma
                              tau
                              (R2S1RP x y (theta + pi) gamma)
                              (R2S1RP z 0 ori gamma)
                                                -- point
                          )
                         0
                         orientations
                                       -- p3 =
                                       --   L.foldl'
                                       --     (\s ori ->
                                       --        s +
                                       --        computePji
                                       --          sigma
                                       --          tau
                                       --          (R2S1RP 0 0 0 gamma)
                                       --          (R2S1RP z 0 ori gamma) *
                                       --        computePji
                                       --          sigma
                                       --          tau
                                       --          (R2S1RP z 0 theta gamma)
                                       --          point)
                                       --     0
                                       --     orientations
                  in p11 * p2

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field