text:
{-# LANGUAGE AllowAmbiguousTypes  #-}
{-# LANGUAGE UndecidableInstances #-}
{-# OPTIONS_GHC -fno-warn-deprecations #-}

module Numeric.Tensile.Tensor.Internal where

import Data.Bits
import Data.ByteString (ByteString())
import Data.Complex
import Data.Int
import Data.Kind
import Data.Proxy
import Data.Word
import Data.Vector.Storable (Storable(..))
import Unsafe.Coerce (unsafeCoerce)

import Numeric.Tensile.Dimensions

import Data.Vector.Sized (Vector)
import TensorFlow.Types
import System.IO.Unsafe (unsafePerformIO)

import qualified Data.Finite as F
import qualified Data.Vector as V
import qualified Data.Vector.Sized as N
import qualified Data.Vector.Storable as S
import qualified Data.Vector.Storable.Mutable as M
import qualified TensorFlow.Tensor as T
import qualified TensorFlow.Build as T
import qualified TensorFlow.GenOps.Core as O
import qualified TensorFlow.Ops as O
import qualified TensorFlow.Session as Ss

import Control.Monad.IO.Class (liftIO)
import Control.Monad.Trans.Reader

{-
fromTensor :: Elt e => Tensor d e -> [e]
fromTensor t = unsafePerformIO $ Ss.runSession $ do
  a <- T.render $ unTensor t
  b <- Ss.run a
  return $ V.toList b

:: OneOf '[Int32, Int64] tidx
=> Tensor v'1 Bool
input

-> Tensor v'2 tidx
reduction_indices

-> Tensor Build Bool

transpose
  :: Elt e
  => Permutable d d'
  => Dims d -> Perm (Rank d) -> Tensor d e -> Tensor d' e
transpose d (Perm p) (Tensor t) = Tensor $ (flip O.transpose) (O.vector w) t
  where v = [1.. fromIntegral $ rank d] :: [IVal]
        w = P.permuteList p v

-}
-- | Generally specialized as in 'TVal' or 'IVal'.
type Elt e = OneOf '[Int32, Int64, Word32, Word64, Float, Double, Complex Float, Complex Double] e

-- TODO: move to application / test stanza
--type Elt = TensorType -- TODO deal with Variant and exclude ResourceHandle
type TVal = Float
type IVal = Int32
type BVal = Bool

type Dynamic d = Reader (Dims d)

newtype Tensor (d :: [Nat]) (e :: Type) = Tensor { unTensor :: T.Tensor T.Build e }

-- | A real or complex-valued tensor of shape 'd'.
type T d = Tensor d TVal

-- | An integer or non-negative integer-valued tensor of shape 'd'.
type I d = Tensor d IVal

-- | A boolean-valued tensor of shape 'd'.
type B d = Tensor d BVal

--TODO update Show instance
instance (Elt e, Show e) => Show (Tensor d e) where show _ = "No Show instance."

--TODO implement
instance (KnownDims d, Elt e, Eq e) => Eq (Tensor d e) where
    (==) ta tb = all $ ta `eq` tb
      where all :: B d -> Bool
            all (Tensor b) = Prelude.head $ toList $ Tensor $ O.all b (O.vector ([] :: [IVal]))

instance (KnownDims d, Num e, Elt e) => Num (Tensor d e) where
    (+) (Tensor v1) (Tensor v2) = Tensor $ (+) v1 v2
    {-# INLINE (+) #-}
    (-) (Tensor v1) (Tensor v2) = Tensor $ (-) v1 v2
    {-# INLINE (-) #-}
    (*) (Tensor v1) (Tensor v2) = Tensor $ (*) v1 v2
    {-# INLINE (*) #-}
    negate (Tensor v1) = Tensor $ negate v1
    {-# INLINE negate #-}
    abs (Tensor v1) = Tensor $ abs v1
    {-# INLINE abs #-}
    signum (Tensor v1) = Tensor $ signum v1
    {-# INLINE signum #-}
    fromInteger = constant (dims @d) . fromInteger
    {-# INLINE fromInteger #-}

{-
instance (KnownDim (Size d), Fractional e, Elt e) => Fractional (Tensor d e)  where
    (/) = _lift2 (/)
    {-# INLINE (/) #-}
    recip = _lift recip
    {-# INLINE recip #-}
    fromRational = _replicate (fromIntegral . dimVal $ (dim :: Dim (Size d))) . fromRational
    {-# INLINE fromRational #-}

instance (KnownDim (Size d), Floating e, Elt e) => Floating (Tensor d e) where
    pi = Tensor $ S.singleton pi  --TODO make this dim safe
    {-# INLINE pi #-}
    exp = _lift exp
    {-# INLINE exp #-}
    log = _lift log
    {-# INLINE log #-}
    sqrt = _lift sqrt
    {-# INLINE sqrt #-}
    sin = _lift sin
    {-# INLINE sin #-}
    cos = _lift cos
    {-# INLINE cos #-}
    tan = _lift tan
    {-# INLINE tan #-}
    asin = _lift asin
    {-# INLINE asin #-}
    acos = _lift acos
    {-# INLINE acos #-}
    atan = _lift atan
    {-# INLINE atan #-}
    sinh = _lift sinh
    {-# INLINE sinh #-}
    cosh = _lift cosh
    {-# INLINE cosh #-}
    tanh = _lift tanh
    {-# INLINE tanh #-}
    (**) = _lift2 (**)
    {-# INLINE (**) #-}
    logBase = _lift2 logBase
    {-# INLINE logBase #-}
    asinh = _lift asinh
    {-# INLINE asinh #-}
    acosh = _lift acosh
    {-# INLINE acosh #-}
    atanh = _lift atanh
    {-# INLINE atanh #-}

instance (KnownDim (Size d), Elt e, Eq e, Bits e, Num e) => Bits (Tensor d e) where
    (.&.) = _lift2 (.&.)
    {-# INLINE (.&.) #-}
    (.|.) = _lift2 (.|.)
    {-# INLINE (.|.) #-}
    xor = _lift2 xor
    {-# INLINE xor #-}
    complement = _lift complement
    shift t i = _lift (flip shift i) t
    rotate t i = _lift (flip rotate i) t
    bit = _replicate (fromIntegral . dimVal $ (dim :: Dim (Size d))) . bit
    testBit = testBitDefault
    bitSizeMaybe _ = bitSizeMaybe @e undefined
    bitSize _ = bitSize @e undefined
    isSigned _ = isSigned @e undefined
    popCount = popCountDefault
-}

{-

instance (KnownDim (Size d), Elt e, Real e) => Real (Tensor d e) where
  toRational = undefined --TODO find a reasonable sum-based implementation or scrap the typeclass

instance (KnownDim (Size d), Elt e, Enum e) => Enum (Tensor d e) where
  toEnum = Tensor . S.replicate (fromIntegral . dimVal $ (dim :: Dim (Size d))) . toEnum
  {-# INLINE toEnum #-}
  fromEnum = undefined --TODO find a reasonable sum-based implementation or scrap the typeclass

instance (KnownDim (Size d), Elt e, Integral e) => Integral (Tensor d e) where
  quot (Tensor a) (Tensor b) = _lift2 quot a b
  rem (Tensor a) (Tensor b) = _lift2 rem a b
  div (Tensor a) (Tensor b) = _lift2 div a b
  mod (Tensor a) (Tensor b) = _lift2 mod a b
  quotRem ta tb = (quot ta tb, rem ta tb)
  divMod ta tb = (div ta tb, mod ta tb)
  toInteger _ = undefined --TODO find a reasonable sum-based implementation or scrap the typeclass

-}

eq :: Elt e => Eq e => Tensor d e -> Tensor d e -> Tensor d BVal
eq (Tensor a) (Tensor b) = Tensor $ O.equal a b

neq :: Elt e => Eq e => Tensor d e -> Tensor d e -> Tensor d BVal
neq (Tensor a) (Tensor b) = Tensor $ O.notEqual a b

lt :: Elt e => Ord e => Tensor d e -> Tensor d e -> Tensor d BVal
--lt :: Ord TVal => T d -> T d -> B d
lt (Tensor a) (Tensor b) = Tensor $ O.less a b

lte :: Elt e => Ord e => Tensor d e -> Tensor d e -> Tensor d BVal
lte (Tensor a) (Tensor b) = Tensor $ O.lessEqual a b

gt :: Elt e => Ord e => Tensor d e -> Tensor d e -> Tensor d BVal
gt (Tensor a) (Tensor b) = Tensor $ O.greater a b

gte :: Elt e => Ord e => Tensor d e -> Tensor d e -> Tensor d BVal
gte (Tensor a) (Tensor b) = Tensor $ O.greaterEqual a b

maximum :: Elt e => Ord e => Tensor d e -> Tensor d e -> Tensor d e
maximum (Tensor a) (Tensor b) = Tensor $ O.maximum a b

minimum :: Elt e => Ord e => Tensor d e -> Tensor d e -> Tensor d e
minimum (Tensor a) (Tensor b) = Tensor $ O.minimum a b

toList :: TensorDataType V.Vector a => Tensor d a -> [a]
toList t = unsafePerformIO $ Ss.runSession $ do
  a <- T.render $ unTensor t
  b <- Ss.run a
  return $ V.toList b

fromList :: forall d e. Elt e => Dims d -> [e] -> Maybe (Tensor d e)
fromList d l = if size d == (fromIntegral $ length l) then Just $ f l else Nothing
  where f = Tensor . O.constant (toShape d)

fromList'' :: forall d e. Elt e => Dynamic d ([e] -> Maybe (Tensor d e))
fromList'' = reader fromList

fromList' :: forall d e. Elt e => KnownDims d => [e] -> Maybe (Tensor d e)
fromList' = withDyn @d fromList''

-- reflectDims @d $ \d -> runReader fromList'' d $ l

liftDyn :: (Dims d -> r) -> Dynamic d r
liftDyn f = reader f

withDyn :: forall d r. KnownDims d => Dynamic d r -> r
withDyn dyn = reflectDims @d $ runReader dyn

{-
 - http://hackage.haskell.org/package/contravariant-1.5/docs/Data-Functor-Contravariant.html
 -
=> f (Dims d)

TODO:
 - do lowerPerm and withDyn have anything in common? (e.g. operator constraints 'lowering' to Predicate (Dims d)). both are functors?
 - type-level Predicate for '[Nat] ? Equivalence '[Nat] ~ ('[Nat],'[Nat]) ~> Bool   Could also include >1 req?
 -
 -
withDyn :: forall d r. KnownDims d => Predicate (Dims d) -> Dynamic d r -> Maybe r

withDyn2
  :: forall x y r. (KnownDims x, KnownDims y)
  => Predicate (Dims x,)
  => Predicate (Dims y)
  -> Dynamic d r -> Maybe r

divided :: Divisible f => f a -> f b -> f (a, b) Source#
divided = divide id

Predicate (Tensor d e)
Equivalence (Tensor d e)
Predicate (Dims d)
Equivalence (Dims d)

divided :: Divisible f => f a -> f b -> f (a, b)

note in general we have

Dims x -> Bool
Dims x -> Dims y -> Bool

type Dynamic d = Reader (Dims d)

fromList :: forall d e. Elt e => KnownDims d => [e] -> Maybe (Tensor d e)
fromList v = if size (dims @d) == (fromIntegral $ length v) then Just $ f v else Nothing
  where f = Tensor . O.constant (toShape (dims @d))

i = 0 :+ 1 :+ 2 :+ S :: Idxs '[1,2,3]
j = 0 :+ 1 :+ 1 :+ S :: Idxs '[1,2,3]

-- fromVector :: Elt e => Dims d -> Vector e -> Maybe (Tensor d e)
-- fromVector d v = undefined
--
-- toVector :: Elt e => Tensor d e -> Vector e
-- toVector = unTensor
--
fromSizedVector :: Elt e => Vector (Size d) e -> Tensor d e
fromSizedVector = Tensor . S.convert . N.fromSized

toSizedVector :: Elt e => Tensor d e -> Vector (Size d) e
toSizedVector = coerce . S.convert . unTensor
  where coerce :: S.Vector e -> N.Vector n e
        coerce = unsafeCoerce

fill :: Elt e => Dims d -> (Idxs d -> e) -> Tensor d e
fill d f = Tensor $ S.create $ do
  mv <- M.new $ fromIntegral $ product $ listDims d
  let act ix = M.write mv (minorToMajor d ix) $ f ix
  forMIdxs_ d act
  return mv

-- constant :: TensorType a => Shape -> [a] -> Tensor Build a

constant :: Elt e => Dims d -> e -> Tensor d e
constant d t = fill d $ const t

-}

constant :: Elt e => Dims d -> e -> Tensor d e
constant d = Tensor . O.constant (toShape d) . pure

toShape :: Dims d -> Shape
toShape = Shape . fmap fromIntegral . listDims

toShape' :: Idxs d -> Shape
toShape' = Shape . fmap fromIntegral . listIdxs

toShapeTensor :: Idxs d -> I '[Rank d]
toShapeTensor i = Tensor . O.constant (toShape' i) $ l
  where l = fmap fromIntegral . listIdxs $ i

-- TODO use more efficient density list approach
fill :: Elt e => Dims d -> (Idxs d -> e) -> Tensor d e
fill d k = Tensor . O.constant (toShape d) $ Prelude.reverse l
  where l = foldIdxs d (\i x -> k i : x) []

{-

fromIntegral' :: Int -> Int32
fromIntegral' = fromIntegral

toList $ fill (dims @'[2,4]) (fromIntegral' . fromEnum)

t1 = fill (dims @'[2,4]) $ const 1 . fromEnum
t2 = fill (dims @'[4,2]) $ const 1 . fromEnum

t2 = fill (dims @'[2,4]) ((+1) . fromIntegral' . fromEnum)

f :: Elt e => Tensor '[2,4] e -> Tensor '[4,2] e
f = transpose (dims @'[2,4]) (reversal (dims @'[2,4]))

t1' = f t1

foo :: Dims d -> (Idxs d -> TVal) -> IO ()
foo d f = printTensor $ fill d f

> d = dims @'[2,4]
> foo d (minorToMajor d)
[7,6,5,4,3,2,1,0]
> foo d (transposeIdxs $ majorToMinor d)
[7,5,3,1,6,4,2,0]

pred_sum_idxs :: forall ds. Dims ds -> Bool
pred_sum_idxs ds = foldIdxs ds (\_ a -> a + 1) 0 == (product . listDims $ ds)

-}

{-

Splits a tensor into sub tensors.

If num_or_size_splits is an integer type, then value is split along dimension axis into num_split smaller tensors. Requires that num_split evenly divides value.shape[axis].

If num_or_size_splits is not an integer type, it is presumed to be a Tensor size_splits, then splits value into len(size_splits) pieces. The shape of the i-th piece has the same size as the value except along dimension axis where the size is size_splits[i].

slice:
t = tf.constant([[[1, 1, 1], [2, 2, 2]],
                 [[3, 3, 3], [4, 4, 4]],
                 [[5, 5, 5], [6, 6, 6]]]) 3 x 2 x 3
tf.slice(t, [1, 0, 0], [1, 1, 3])  # [[[3, 3, 3]]]
tf.slice(t, [1, 0, 0], [1, 2, 3])  # [[[3, 3, 3],
                                   #   [4, 4, 4]]]
tf.slice(t, [1, 0, 0], [2, 1, 3])  # [[[3, 3, 3]],
                                   #  [[5, 5, 5]]]

-- want d' < d,
slice :: Idxs d -> Dims d' -> Tensor d e -> Tensor d' e

slice :: (TensorType t, OneOf '[Int32, Int64] index)	=> Tensor v'1 t	-> Tensor v'2 index	-> Tensor v'3 index	-> Tensor Build t

Idxs d -> I '[Rank d]

i = toEnum 6 :: Idxs '[3,2,3]
j = toEnum 9 :: Idxs '[3,2,3]

diffIdxs (dims @'[3,2,3]) i j

stepIdx (dims @'[3,2,3]) 3 i
i = 1 :+ 0 :+ 0 :+ S :: Idxs '[3,2,3]

-}

{-

pack
  :: Elt e
  => Vector n (Tensor (x ++ y) e) -> Tensor (x ++ n :+ y) e
pack = undefined

pack0
  :: forall d e n. Elt e
  => KnownDims d
  => KnownDim n
  => Vector n (Tensor d e) -> Tensor (n :+ d) e
pack0 v = Tensor res
  where d = dims @d
        n = fromIntegral $ dimVal $ dim @n
        size = product $ listDims d
        res = S.create $ do
          mv <- M.new $ fromIntegral $ size * n
          flip N.imapM_ v $ \i t ->
            let i' = idxToWord . idxFromFinite $ i
                off = fromIntegral $ i' * size
                v' = unTensor t
                act ix = M.write mv (off + fromEnum ix) $ v' S.! (fromEnum ix) -- could use a tensor op instead here
            in forMIdxs_ d act
          return mv

-- TODO use http://hackage.haskell.org/package/reflection-2.1.4/docs/Data-Reflection.html#t:reifyNat
unpack0
  :: forall d e n. Elt e
  => KnownDims d
  => KnownNat n
  => Tensor (n :+ d) e -> Vector n (Tensor d e)
unpack0 t = N.generate f
  where d = dims @d
        size = fromIntegral $ product $ listDims d
        f i = fill d $ \ix ->
          let i' = fromIntegral $ F.getFinite i
              off = i' * size
              v = unTensor t
          in v S.! (off + fromEnum ix)
-}

{-

see example http://jeankossaifi.com/blog/unfolding.html

t :: Vector 4 (Tensor '[2,2] Word)
t = N.generate $ \f ->
  let d = dims @'[2,2]
      i' = idxToWord . idxFromFinite $ f
  in fill d (const i')

t' :: Tensor '[4,2,2] Word
t' = pack0 t

t'' :: Vector 4 (Tensor '[2,2] Word)
t'' = unpack0 t'

generate :: forall n a. KnownNat n => (Finite n -> a) -> Vector n a

t :: Data.Vector.Sized.Vector 4 (Tensor '[2,2] Word)
t = generate $ \f ->
  let d = dims @'[2,2]
      i' = idxToWord . idxFromFinite $ f
  in fill d (const i')

t' :: Tensor '[4,2,2] Word
t' = pack0 t
N.imapM_ :: Monad m => (Finite n -> a -> m b) -> Vector n a -> m ()

overDim_ :: Monad m
         => Dims ds -- ^ Shape of a space
         -> (Idxs ds -> Int -> m ()) -- ^ Function to call on each dimension
         -> Int -- ^ Initial offset
         -> Int -- ^ Offset step
         -> m ()

stack ::  Dim n -> Vector n (Tensor (x ++ y) e) -> Tensor (x +: n :+ y) e
unstack :: (KnownDims x, Elt e) => Dim n -> Tensor (x +: n :+ y) e -> Vector n (Tensor (x ++ y) e)
-}

max_stars_repo_name: cmk/tensile

max_issues_repo_name: cmk/tensile

max_forks_repo_name: cmk/tensile


text:
-- ch1 ex.1.23

import Statistics.Distribution.Normal
import qualified Statistics.Distribution as S

normalD :: NormalDistribution
normalD = normalDistr mean stdD
  where
    mean = 75
    stdD = 20

main :: IO ()
main = do
  -- ratio greater than 95 sec:
  putStrLn $ "the ratio of longer than 95 second: "
          ++ show (S.complCumulative normalD 95)

  putStrLn $ "Between 35 sec and 115 sec:         "
          ++ show ((S.complCumulative normalD 35) - (S.complCumulative normalD 115))

  putStrLn $ "more than 2 min:                    "
          ++ show (S.complCumulative normalD 120)

ext: hs

max_stars_repo_name: raydsameshima/Stat

max_issues_repo_name: raydsameshima/Stat

max_forks_repo_name: raydsameshima/Stat


text:
-----------------------------------------------------------------------------
-- |
-- Module     : Numeric.LinearAlgebra.Packed
-- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : experimental
--
-- Packed matrices.
--
module Numeric.LinearAlgebra.Packed (
    -- * Immutable packed matrices
    Packed,
    dim,

    -- * Read-only packed matrices
    RPacked(..),

    -- * Conversions between vectors and packed matrices
    fromVector,
    toVector,

    -- * Mutable interface
    module Numeric.LinearAlgebra.Packed.ST,

    -- * Hermitian views
    module Numeric.LinearAlgebra.Packed.Herm,

    -- * Triangular views
    module Numeric.LinearAlgebra.Packed.Tri,

    -- * Cholesky factorizations
    module Numeric.LinearAlgebra.Packed.Cholesky,

    -- * Basic multivariate statistics
    module Numeric.LinearAlgebra.Packed.Statistics,

    ) where

import Numeric.LinearAlgebra.Packed.Base
import Numeric.LinearAlgebra.Packed.ST hiding ( RPacked(..) )
import Numeric.LinearAlgebra.Packed.Herm
import Numeric.LinearAlgebra.Packed.Tri
import Numeric.LinearAlgebra.Packed.Cholesky
import Numeric.LinearAlgebra.Packed.Statistics

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies, MagicHash,
  FlexibleInstances, FlexibleContexts, UnboxedTuples, DeriveDataTypeable, CPP,
  BangPatterns #-}
#ifdef __GLASGOW_HASKELL__
#if __GLASGOW_HASKELL__ < 610
{-# OPTIONS_GHC -frewrite-rules #-}
#else
{-# OPTIONS_GHC -fenable-rewrite-rules #-}
#endif
#endif
-----------------------------------------------------------------------------
-- |
-- Module      : Data.Array.CArray.Base
-- Copyright   : (c) 2001 The University of Glasgow
--               (c) 2008 Jed Brown
-- License     : BSD-style
--
-- Maintainer  : jed@59A2.org
-- Stability   : experimental
-- Portability : non-portable
--
-- This module provides both the immutable 'CArray' and mutable 'IOCArray'.  The
-- underlying storage is exactly the same - pinned memory on the GC'd heap.
-- Elements are stored according to the class 'Storable'.  You can obtain a
-- pointer to the array contents to manipulate elements from languages like C.
--
-- 'CArray' is 16-byte aligned by default.  If you create a 'CArray' with
-- 'unsafeForeignPtrToCArray' then it may not be aligned.  This will be an issue
-- if you intend to use SIMD instructions.
--
-- 'CArray' is similar to 'Data.Array.Unboxed.UArray' but slower if you stay
-- within Haskell.  'CArray' can handle more types and can be used by external
-- libraries.
--
-- 'IOCArray' is equivalent to 'Data.Array.Storable.StorableArray' and similar
-- to 'Data.Array.IO.IOUArray' but slower.  'IOCArray' has O(1) versions of
-- 'unsafeFreeze' and 'unsafeThaw' when converting to/from 'CArray'.
-----------------------------------------------------------------------------

module Data.Array.CArray.Base where

import Control.Applicative
import Control.Monad
import Data.Ix.Shapable
import Data.Array.Base
import Data.Array.MArray
import Data.Array.IArray        ()
import qualified Data.ByteString.Internal as S
import Data.Binary
import Data.Complex
import Data.List
import System.IO.Unsafe         (unsafePerformIO)
import Foreign.Storable
import Foreign.ForeignPtr
import Foreign.Ptr
import Foreign.Marshal.Alloc
import Foreign.Marshal.Array    (copyArray)
import Data.Word                ()

import Data.Generics            (Data(..), Typeable(..))
import GHC.Ptr                  (Ptr(..))
import GHC.ForeignPtr           (ForeignPtr(..), mallocPlainForeignPtrBytes)

-- | The immutable array type.
data CArray i e = CArray !i !i Int !(ForeignPtr e)
    deriving (Data, Typeable)

-- | Absolutely equivalent representation, but used for the mutable interface.
data IOCArray i e = IOCArray !i !i Int !(ForeignPtr e)
    deriving (Data, Typeable)

instance Storable e => MArray IOCArray e IO where
    getBounds (IOCArray l u _ _) = return (l,u)

    getNumElements (IOCArray _ _ n _) = return n

    newArray (l,u) e0 = do
        fp <- mallocForeignPtrArrayAligned n
        withForeignPtr fp $ \a ->
            sequence_ [pokeElemOff a i e0 | i <- [0 .. n - 1]]
        return (IOCArray l u n fp)
        where n = rangeSize (l,u)

    unsafeNewArray_ (l,u) = do
        let n = rangeSize (l,u)
        fp <- mallocForeignPtrArrayAligned n
        return (IOCArray l u n fp)

    newArray_ = unsafeNewArray_

    unsafeRead (IOCArray _ _ _ fp) i =
        withForeignPtr fp $ \a -> peekElemOff a i

    unsafeWrite (IOCArray _ _ _ fp) i e =
        withForeignPtr fp $ \a -> pokeElemOff a i e

-- | The pointer to the array contents is obtained by 'withCArray'.
-- The idea is similar to 'ForeignPtr' (used internally here).
-- The pointer should be used only during execution of the 'IO' action
-- retured by the function passed as argument to 'withCArray'.
withCArray :: CArray i e -> (Ptr e -> IO a) -> IO a
withCArray (CArray _ _ _ fp) f = withForeignPtr fp f

withIOCArray :: IOCArray i e -> (Ptr e -> IO a) -> IO a
withIOCArray (IOCArray _ _ _ fp) f = withForeignPtr fp f

-- | If you want to use it afterwards, ensure that you
-- 'touchCArray' after the last use of the pointer,
-- so the array is not freed too early.
touchIOCArray :: IOCArray i e -> IO ()
touchIOCArray (IOCArray _ _ _ fp) = touchForeignPtr fp

-- | /O(1)/ Construct a 'CArray' from an arbitrary 'ForeignPtr'.  It is
-- the caller's responsibility to ensure that the 'ForeignPtr' points to
-- an area of memory sufficient for the specified bounds.
unsafeForeignPtrToCArray
   :: Ix i => ForeignPtr e -> (i,i) -> IO (CArray i e)
unsafeForeignPtrToCArray p (l,u) =
   return (CArray l u (rangeSize (l,u)) p)

-- | /O(1)/ Construct a 'CArray' from an arbitrary 'ForeignPtr'.  It is
-- the caller's responsibility to ensure that the 'ForeignPtr' points to
-- an area of memory sufficient for the specified bounds.
unsafeForeignPtrToIOCArray
   :: Ix i => ForeignPtr e -> (i,i) -> IO (IOCArray i e)
unsafeForeignPtrToIOCArray p (l,u) =
   return (IOCArray l u (rangeSize (l,u)) p)

-- | /O(1)/ Extract ForeignPtr from a CArray.
toForeignPtr :: CArray i e -> (Int, ForeignPtr e)
toForeignPtr (CArray _ _ n fp) = (n, fp)

-- | /O(1)/ Turn a CArray into a ByteString.  Unsafe because it uses
-- 'castForeignPtr' and thus is not platform independent.
unsafeCArrayToByteString :: (Storable e) => CArray i e -> S.ByteString
unsafeCArrayToByteString (CArray _ _ l fp) = go undefined fp
    where go :: (Storable e) => e -> ForeignPtr e -> S.ByteString
          go dummy fp' = S.fromForeignPtr (castForeignPtr fp') 0 (l * sizeOf dummy)

-- | /O(1)/ Turn a ByteString into a CArray.  Unsafe because it uses
-- 'castForeignPtr' and thus is not platform independent.  Returns 'Nothing' if
-- the range specified is larger than the size of the ByteString or the start of
-- the ByteString does not fulfil the alignment requirement of the resulting
-- CArray (as specified by the Storable instance).
unsafeByteStringToCArray :: (Ix i, Storable e, IArray CArray e)
                            => (i,i) -> S.ByteString -> Maybe (CArray i e)
unsafeByteStringToCArray lu bs = go undefined lu
    where go :: (Ix i, Storable e, IArray CArray e) => e -> (i,i) -> Maybe (CArray i e)
          go dummy (l,u) | safe = Just (CArray l u n fp)
                         | otherwise = Nothing
              where n = rangeSize (l,u)
                    !((ForeignPtr addr contents), off, len) = S.toForeignPtr bs
                    !p@(Ptr addr') = Ptr addr `plusPtr` off
                    fp = ForeignPtr addr' contents
                    safe = sizeOf dummy * n <= len && p == p `alignPtr` alignment dummy

copy :: (Ix i, Storable e) => CArray i e -> IO (CArray i e)
copy ain@(CArray l u n _) =
    createCArray (l,u) $ \op ->
        withCArray ain $ \ip ->
            copyArray op ip n

freezeIOCArray :: (Ix i, Storable e) => IOCArray i e -> IO (CArray i e)
freezeIOCArray = unsafeFreezeIOCArray >=> copy

unsafeFreezeIOCArray :: (Ix i) => IOCArray i e -> IO (CArray i e)
unsafeFreezeIOCArray (IOCArray l u n fp) = return (CArray l u n fp)

thawIOCArray :: (Ix i, Storable e) => CArray i e -> IO (IOCArray i e)
thawIOCArray = copy >=> unsafeThawIOCArray

unsafeThawIOCArray :: (Ix i) => CArray i e -> IO (IOCArray i e)
unsafeThawIOCArray (CArray l u n fp) = return (IOCArray l u n fp)

-- Since we can remove the (Storable e) restriction for these, the rules are
-- compact and general.
{-# RULES
"unsafeFreeze/IOCArray" unsafeFreeze = unsafeFreezeIOCArray
"unsafeThaw/IOCArray" unsafeThaw = unsafeThawIOCArray
  #-}

-- Since we can't parameterize the rules with the (Storable e) constraint, we
-- have to specialize manually.  This is unfortunate since it is less general.
{-# RULES
"freeze/IOCArray/Int"    freeze      = freezeIOCArray :: (Ix i) => IOCArray i Int -> IO (CArray i Int)
"freeze/IOCArray/Float"  freeze      = freezeIOCArray :: (Ix i) => IOCArray i Float -> IO (CArray i Float)
"freeze/IOCArray/Double" freeze      = freezeIOCArray :: (Ix i) => IOCArray i Double -> IO (CArray i Double)
"thaw/IOCArray/Int"      thaw        = thawIOCArray   :: (Ix i) => CArray i Int -> IO (IOCArray i Int)
"thaw/IOCArray/Float"    thaw        = thawIOCArray   :: (Ix i) => CArray i Float -> IO (IOCArray i Float)
"thaw/IOCArray/Double"   thaw        = thawIOCArray   :: (Ix i) => CArray i Double -> IO (IOCArray i Double)
  #-}

instance Storable e => IArray CArray e where
    {-# INLINE bounds #-}
    bounds (CArray l u _ _) = (l,u)
    {-# INLINE numElements #-}
    numElements (CArray _ _ n _) = n
    {-# NOINLINE unsafeArray #-}
    unsafeArray lu ies = unsafePerformIO $ unsafeArrayCArray lu ies (zeroElem (undefined :: e))
    {-# INLINE unsafeAt #-}
    unsafeAt (CArray _ _ _ fp) i = S.inlinePerformIO $
                                   withForeignPtr fp $ \a -> peekElemOff a i
    {-# NOINLINE unsafeReplace #-}
    unsafeReplace arr ies = unsafePerformIO $ unsafeReplaceCArray arr ies
    {-# NOINLINE unsafeAccum #-}
    unsafeAccum f arr ies = unsafePerformIO $ unsafeAccumCArray f arr ies
    {-# NOINLINE unsafeAccumArray #-}
    unsafeAccumArray f e0 lu ies = unsafePerformIO $ unsafeAccumArrayCArray f e0 lu ies

-- | Hackish way to get the zero element for a Storable type.
{-# NOINLINE zeroElem #-}
zeroElem :: Storable a => a -> a
zeroElem u = unsafePerformIO $ do
    allocaBytes n $ \p -> do
        sequence_ [pokeByteOff p off (0 :: Word8) | off <- [0 .. n - 1] ]
        peek (castPtr p)
    where n = sizeOf u

{-# INLINE unsafeArrayCArray #-}
unsafeArrayCArray :: (MArray IOCArray e IO, Storable e, Ix i)
                  => (i,i) -> [(Int, e)] -> e -> IO (CArray i e)
unsafeArrayCArray lu ies default_elem = do
        marr <- newArray lu default_elem
        sequence_ [unsafeWrite marr i e | (i, e) <- ies]
        unsafeFreezeIOCArray marr

{-# INLINE unsafeReplaceCArray #-}
unsafeReplaceCArray :: (MArray IOCArray e IO, Storable e, Ix i)
                       => CArray i e -> [(Int, e)] -> IO (CArray i e)
unsafeReplaceCArray arr ies = do
    marr <- thawIOCArray arr
    sequence_ [unsafeWrite marr i e | (i, e) <- ies]
    unsafeFreezeIOCArray marr

{-# INLINE unsafeAccumCArray #-}
unsafeAccumCArray :: (MArray IOCArray e IO, Storable e, Ix i)
                            => (e -> e' -> e) -> CArray i e -> [(Int, e')]
                                              -> IO (CArray i e)
unsafeAccumCArray f arr ies = do
    marr <- thawIOCArray arr
    sequence_ [do
        old <- unsafeRead marr i
        unsafeWrite marr i (f old new)
        | (i, new) <- ies]
    unsafeFreezeIOCArray marr

{-# INLINE unsafeAccumArrayCArray #-}
unsafeAccumArrayCArray :: (MArray IOCArray e IO, Storable e, Ix i)
                          => (e -> e' -> e) -> e -> (i,i) -> [(Int, e')]
                                            -> IO (CArray i e)
unsafeAccumArrayCArray f e0 lu ies = do
    marr <- newArray lu e0
    sequence_ [do
        old <- unsafeRead marr i
        unsafeWrite marr i (f old new)
        | (i, new) <- ies]
    unsafeFreezeIOCArray marr

{-# INLINE eqCArray #-}
eqCArray :: (IArray CArray e, Ix i, Eq e)
                   => CArray i e -> CArray i e -> Bool
eqCArray arr1@(CArray l1 u1 n1 _) arr2@(CArray l2 u2 n2 _) =
    if n1 == 0 then n2 == 0 else
        l1 == l2 && u1 == u2 &&
           and [unsafeAt arr1 i == unsafeAt arr2 i | i <- [0 .. n1 - 1]]

{-# INLINE cmpCArray #-}
cmpCArray :: (IArray CArray e, Ix i, Ord e)
                    => CArray i e -> CArray i e -> Ordering
cmpCArray arr1 arr2 = compare (assocs arr1) (assocs arr2)

{-# INLINE cmpIntCArray #-}
cmpIntCArray :: (IArray CArray e, Ord e, Storable e)
                       => CArray Int e -> CArray Int e -> Ordering
cmpIntCArray arr1@(CArray l1 u1 n1 _) arr2@(CArray l2 u2 n2 _) =
    if n1 == 0 then if n2 == 0 then EQ else LT else
    if n2 == 0 then GT else
    case compare l1 l2 of
        EQ    -> foldr cmp (compare u1 u2) [0 .. (n1 `min` n2) - 1]
        other -> other
    where
    cmp i rest = case compare (unsafeAt arr1 i) (unsafeAt arr2 i) of
        EQ    -> rest
        other -> other

-- {-# RULES "cmpCArray/Int" cmpCArray = cmpIntCArray #-}

instance (Ix ix, Eq e, IArray CArray e) => Eq (CArray ix e) where
    (==) = eqCArray

instance (Ix ix, Ord e, IArray CArray e) => Ord (CArray ix e) where
    compare = cmpCArray

instance (Ix ix, Show ix, Show e, IArray CArray e) => Show (CArray ix e) where
    showsPrec = showsIArray

--
-- General purpose array operations which happen to be very fast for CArray.
--

-- | O(1) reshape an array.  The number of elements in the new shape must not
-- exceed the number in the old shape.  The elements are in C-style ordering.
reshape :: (Ix i, Ix j) => (j,j) -> CArray i e -> CArray j e
reshape (l',u') (CArray _ _ n fp) | n' > n = error "reshape: new size too large"
                                  | otherwise = CArray l' u' n' fp
    where n' = rangeSize (l', u')

-- | O(1) make a rank 1 array from an arbitrary shape.
-- It has the property that 'reshape (0, size a - 1) a == flatten a'.
flatten :: Ix i => CArray i e -> CArray Int e
flatten (CArray _ _ n fp) = CArray 0 (n - 1) n fp

--
-- None of the following are specific to CArray.  Some could have slightly
-- faster versions specialized to CArray.  In general, slicing is expensive
-- because the slice is not contiguous in memory, so must be copied.  There are
-- many specialized versions.
--

-- | Generic slice and map.  This takes the new range, the inverse map on
-- indices, and function to produce the next element.  It is the most general
-- operation in its class.
ixmapWithIndP :: (Ix i, Ix i', IArray a e, IArray a' e')
                 => (i',i') -> (i' -> i) -> (i -> e -> i' -> e') -> a i e -> a' i' e'
ixmapWithIndP lu f g arr = listArray lu
                           [ let i = f i' in g i (arr ! i) i' | i' <- range lu ]

-- | Less polymorphic version.
ixmapWithInd :: (Ix i, Ix i', IArray a e, IArray a e')
                => (i',i') -> (i' -> i) -> (i -> e -> i' -> e') -> a i e -> a i' e'
ixmapWithInd = ixmapWithIndP

-- | Perform an operation on the elements, independent of their location.
ixmapWithP  :: (Ix i, Ix i', IArray a e, IArray a' e')
               => (i',i') -> (i' -> i) -> (e -> e') -> a i e -> a' i' e'
ixmapWithP lu f g arr = listArray lu
                        [ g (arr ! f i') | i' <- range lu ]

-- | Less polymorphic version.
ixmapWith  :: (Ix i, Ix i', IArray a e, IArray a e')
              => (i',i') -> (i' -> i) -> (e -> e') -> a i e -> a i' e'
ixmapWith = ixmapWithP

-- | More polymorphic version of 'ixmap'.
ixmapP :: (Ix i, Ix i', IArray a e, IArray a' e)
          => (i',i') -> (i' -> i) -> a i e -> a' i' e
ixmapP lu f arr = ixmapWithP lu f id arr

-- | More friendly sub-arrays with element mapping.
sliceStrideWithP :: (Ix i, Shapable i, Ix i', IArray a e, IArray a' e')
                    => (i',i') -> (i,i,i) -> (e -> e') -> a i e -> a' i' e'
sliceStrideWithP lu (start,next,end) f arr
    | all (inRange (bounds arr)) [start,next,end] = listArray lu es
    | otherwise = error "sliceStrideWith: out of bounds"
    where is = offsetShapeFromThenTo (shape arr) (index' start) (index' next) (index' end)
          es = map (f . (unsafeAt arr)) is
          index' = indexes arr

-- | Less polymorphic version.
sliceStrideWith :: (Ix i, Shapable i, Ix i', IArray a e, IArray a e')
                   => (i',i') -> (i,i,i) -> (e -> e') -> a i e -> a i' e'
sliceStrideWith = sliceStrideWithP

-- | Strided sub-array without element mapping.
sliceStrideP :: (Ix i, Shapable i, Ix i', IArray a e, IArray a' e)
                => (i',i') -> (i,i,i) -> a i e -> a' i' e
sliceStrideP lu sne = sliceStrideWithP lu sne id

-- | Less polymorphic version.
sliceStride :: (Ix i, Shapable i, Ix i', IArray a e)
               => (i',i') -> (i,i,i) -> a i e -> a i' e
sliceStride = sliceStrideP

-- | Contiguous sub-array with element mapping.
sliceWithP :: (Ix i, Shapable i, Ix i', IArray a e, IArray a' e')
              => (i',i') -> (i,i) -> (e -> e') -> a i e -> a' i' e'
sliceWithP lu (start,end) f arr
    | all (inRange (bounds arr)) [start,end] = listArray lu es
    | otherwise = error "sliceWith: out of bounds"
    where is = offsetShapeFromTo (shape arr) (index' start) (index' end)
          es = map (f . (unsafeAt arr)) is
          index' = indexes arr

-- | Less polymorphic version.
sliceWith :: (Ix i, Shapable i, Ix i', IArray a e, IArray a e')
             => (i',i') -> (i,i) -> (e -> e') -> a i e -> a i' e'
sliceWith = sliceWithP

-- | Contiguous sub-array without element mapping.
sliceP :: (Ix i, Shapable i, Ix i', IArray a e, IArray a' e)
          => (i',i') -> (i,i) -> a i e -> a' i' e
sliceP lu se = sliceWithP lu se id

-- | Less polymorphic version.
slice ::  (Ix i, Shapable i, Ix i', IArray a e)
          => (i',i') -> (i,i) -> a i e -> a i' e
slice = sliceP

-- | In-place map on CArray.  Note that this is /IN PLACE/ so you should not
-- retain any reference to the original.  It flagrantly breaks referential
-- transparency!
{-# INLINE mapCArrayInPlace #-}
mapCArrayInPlace :: (Ix i, IArray CArray e, Storable e) => (e -> e) -> CArray i e -> CArray i e
mapCArrayInPlace f a = S.inlinePerformIO $ do
    withCArray a $ \p ->
        forM_ [0 .. size a - 1] $ \i ->
            peekElemOff p i >>= pokeElemOff p i . f
    return a

-----------------------------------------
-- These are meant to be internal only
indexes :: (Ix i, Shapable i, IArray a e) => a i e -> i -> [Int]
indexes a i = map pred $ (sShape . fst . bounds) a i

offsetShapeFromThenTo :: [Int] -> [Int] -> [Int] -> [Int] -> [Int]
offsetShapeFromThenTo s a b c = foldr (liftA2 (+)) [0] (ilists stride a b c)
    where ilists = zipWith4 (\s' a' b' c' -> map (*s') $ enumFromThenTo a' b' c')
          stride = shapeToStride s

offsetShapeFromTo :: [Int] -> [Int] -> [Int] -> [Int]
offsetShapeFromTo = offsetShapeFromTo' id

offsetShapeFromTo' :: ([[Int]] -> [[Int]]) -> [Int] -> [Int] -> [Int] -> [Int]
offsetShapeFromTo' f s a b = foldr (liftA2 (+)) [0] (f $ ilists stride a b)
    where ilists = zipWith3 (\s' a' b' -> map (*s') $ enumFromTo a' b')
          stride = shapeToStride s

offsets :: (Ix a, Shapable a) => (a, a) -> a -> [Int]
offsets lu i = reverse . osets (index lu i) . reverse . scanl1 (*) . uncurry sShape $ lu
    where osets 0 [] = []
          osets i' (b:bs) = r : osets d bs
              where (d,r) = i' `divMod` b
          osets _ _ = error "osets"
-----------------------------------------

-- | p-norm on the array taken as a vector
normp :: (Ix i, RealFloat e', Abs e e', IArray a e) => e' -> a i e -> e'
normp p a | 1 <= p && not (isInfinite p) = (** (1/p)) $ foldl' (\z e -> z + (abs_ e) ** p) 0 (elems a)
          | otherwise = error "normp: p < 1"

-- | 2-norm on the array taken as a vector (Frobenius norm for matrices)
norm2 :: (Ix i, Floating e', Abs e e', IArray a e) => a i e -> e'
norm2 a = sqrt $ foldl' (\z e -> z + abs_ e ^ (2 :: Int)) 0 (elems a)

-- | Sup norm on the array taken as a vector
normSup :: (Ix i, Num e', Ord e', Abs e e', IArray a e) => a i e -> e'
normSup a = foldl' (\z e -> z `max` abs_ e) 0 (elems a)

-- | Polymorphic version of amap.
liftArrayP :: (Ix i, IArray a e, IArray a1 e1)
              => (e -> e1) -> a i e -> a1 i e1
liftArrayP f a = listArray (bounds a) (map f (elems a))

-- | Equivalent to amap.  Here for consistency only.
liftArray :: (Ix i, IArray a e, IArray a e1)
              => (e -> e1) -> a i e -> a i e1
liftArray = liftArrayP

-- | Polymorphic 2-array lift.
liftArray2P :: (Ix i, IArray a e, IArray a1 e1, IArray a2 e2)
              => (e -> e1 -> e2) -> a i e -> a1 i e1 -> a2 i e2
liftArray2P f a b | aBounds == bounds b =
                     listArray aBounds (zipWith f (elems a) (elems b))
                 | otherwise = error "liftArray2: array bounds must match"
    where aBounds = bounds a

-- | Less polymorphic version.
liftArray2 :: (Ix i, IArray a e, IArray a e1, IArray a e2)
              => (e -> e1 -> e2) -> a i e -> a i e1 -> a i e2
liftArray2 = liftArray2P

-- | Polymorphic 3-array lift.
liftArray3P :: (Ix i, IArray a e, IArray a1 e1, IArray a2 e2, IArray a3 e3)
               => (e -> e1 -> e2 -> e3) -> a i e -> a1 i e1 -> a2 i e2 -> a3 i e3
liftArray3P f a b c | aBounds == bounds b && aBounds == bounds c =
                       listArray aBounds (zipWith3 f (elems a) (elems b) (elems c))
                   | otherwise = error "liftArray2: array bounds must match"
    where aBounds = bounds a

-- | Less polymorphic version.
liftArray3 :: (Ix i, IArray a e, IArray a e1, IArray a e2, IArray a e3)
              => (e -> e1 -> e2 -> e3) -> a i e -> a i e1 -> a i e2 -> a i e3
liftArray3 = liftArray3P

-- | Hack so that norms have a sensible type.
class Abs a b | a -> b where
    abs_ :: a -> b
instance Abs (Complex Double) Double where
    abs_ = magnitude
instance Abs (Complex Float) Float where
    abs_ = magnitude
instance Abs Double Double where
    abs_ = abs
instance Abs Float Float where
    abs_ = abs

-- | Allocate an array which is 16-byte aligned.  Essential for SIMD instructions.
mallocForeignPtrArrayAligned :: Storable a => Int -> IO (ForeignPtr a)
mallocForeignPtrArrayAligned n = doMalloc undefined
  where
    doMalloc :: Storable b => b -> IO (ForeignPtr b)
    doMalloc dummy = mallocForeignPtrBytesAligned (n * sizeOf dummy)

-- | Allocate memory which is 16-byte aligned.  This is essential for SIMD
-- instructions.  We know that mallocPlainForeignPtrBytes will give word-aligned
-- memory, so we pad enough to be able to return the desired amount of memory
-- after aligning our pointer.
mallocForeignPtrBytesAligned :: Int -> IO (ForeignPtr a)
mallocForeignPtrBytesAligned n = do
    (ForeignPtr addr contents) <- mallocPlainForeignPtrBytes (n + pad)
    let !(Ptr addr') = alignPtr (Ptr addr) 16
    return (ForeignPtr addr' contents)
    where pad = 16 - sizeOf (undefined :: Word)

-- | Make a new CArray with an IO action.
createCArray :: (Ix i, Storable e) => (i,i) -> (Ptr e -> IO ()) -> IO (CArray i e)
createCArray lu f = do
    fp <- mallocForeignPtrArrayAligned (rangeSize lu)
    withForeignPtr fp f
    unsafeForeignPtrToCArray fp lu

unsafeCreateCArray :: (Ix i, Storable e) => (i,i) -> (Ptr e -> IO ()) -> CArray i e
unsafeCreateCArray lu =  unsafePerformIO . createCArray lu

instance (Ix i, Binary i, Binary e, Storable e) => Binary (CArray i e) where
    put a = do
        put (bounds a)
        mapM_ put (elems a)
    get = do
        lu <- get
        es <- replicateM (rangeSize lu) get
        return $ listArray lu es

max_stars_repo_name: BeautifulDestinations/carray

max_issues_repo_name: BeautifulDestinations/carray

max_forks_repo_name: BeautifulDestinations/carray


text:
-- Copyright 2017 Elkin Cruz.
-- Copyright 2017 James Bowen.
--
-- Licensed under the Apache License, Version 2.0 (the "License");
-- you may not use this file except in compliance with the License.
-- You may obtain a copy of the License at
--
--     http://www.apache.org/licenses/LICENSE-2.0
--
-- Unless required by applicable law or agreed to in writing, software
-- distributed under the License is distributed on an "AS IS" BASIS,
-- WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-- See the License for the specific language governing permissions and
-- limitations under the License.

{-# LANGUAGE RankNTypes           #-}
{-# LANGUAGE KindSignatures       #-}
{-# LANGUAGE DataKinds            #-}
{-# LANGUAGE TypeOperators        #-}
{-# LANGUAGE TypeFamilies         #-}
{-# LANGUAGE ScopedTypeVariables  #-}
{-# LANGUAGE UndecidableInstances #-}
{-# OPTIONS_GHC -Wno-orphans #-}

module TensorFlow.DepTyped.Ops (
  constant,
  placeholder,
  add,
  mul,
  matMul,
  argMax,
  softmax,
  scalar,
  oneHot,
  oneHot_,
  reduceMean,
  softmaxCrossEntropyWithLogits,
  equal,
  truncatedNormal,
  relu,
  sub,
  cast,
  square,
  reshape,
  sigmoid,
  shape
) where

import           GHC.TypeLits (Nat, Symbol, KnownNat, natVal)
import           Data.Proxy (Proxy(Proxy))
import           Data.Kind (Type)

import           Data.Vector.Sized (Vector, toList)
import           Data.Int (Int8, Int16, Int32, Int64)
import           Data.Word (Word8, Word16)
import           Data.Complex (Complex)
import           Data.ByteString (ByteString)

import qualified TensorFlow.Ops as TF (constant, add, matMul, placeholder, argMax, scalar,
                                       softmax, oneHot, reduceMean, softmaxCrossEntropyWithLogits,
                                       equal, truncatedNormal, vector, mul, relu, sub, cast, abs, sign,
                                       neg, reshape, shape)
import qualified TensorFlow.GenOps.Core as TF (square, sigmoid)
import qualified TensorFlow.Types as TF (TensorType, Shape(Shape), type OneOf)
--import qualified TensorFlow.Tensor as TF (Tensor)
import           TensorFlow.Tensor (Value)
import           TensorFlow.Build (Build, MonadBuild)

import           TensorFlow.DepTyped.Base (KnownNatList(natListVal), ShapeProduct, UnionPlaceholder,
                                           BroadcastShapes, RemoveAxisFromShape, AddAxisToEndShape,
                                           KnownNatListLength)
import           TensorFlow.DepTyped.Tensor (Tensor(Tensor), Placeholder)

instance ( TF.TensorType a
         , Num a
         , shps ~ '[]
         , v ~ Build
         , TF.OneOf '[ Double, Float, Int32, Int64
                     , Complex Float, Complex Double] a) => Num (Tensor s shps v a) where
    (Tensor t1) + (Tensor t2) = Tensor (t1 `TF.add` t2)
    (Tensor t1) * (Tensor t2) = Tensor (t1 `TF.mul` t2)
    (Tensor t1) - (Tensor t2) = Tensor (t1 `TF.sub` t2)
    abs (Tensor t)    = Tensor (TF.abs t)
    fromInteger       = Tensor . TF.scalar . fromInteger
    signum (Tensor t) = Tensor (TF.sign t)
    negate (Tensor t) = Tensor (TF.neg t)

constant :: forall (s :: [Nat]) (n :: Nat) a.
            (TF.TensorType a, ShapeProduct s ~ n, KnownNatList s)
         => Vector n a
         -> Tensor s '[] Build a
constant v = Tensor $ TF.constant shape' (toList v)
  where shape' = TF.Shape . fmap fromInteger $ natListVal (Proxy :: Proxy s)

placeholder :: forall (name :: Symbol) (shape :: [Nat]) a m.
               (TF.TensorType a, KnownNatList shape, MonadBuild m)
            => m (Placeholder name shape a)
placeholder = Tensor <$> TF.placeholder shape'
  where shape' = TF.Shape . fmap fromInteger $ natListVal (Proxy :: Proxy shape)

add :: forall (shape1 :: [Nat]) (shape2 :: [Nat]) (phs1 :: [(Symbol, [Nat], Type)]) (phs2 :: [(Symbol, [Nat], Type)]) v1 v2 a.
       TF.OneOf '[(Complex Double), (Complex Float), ByteString, Int16, Int32, Int64, Int8, Word16, Word8, Double, Float] a
    => Tensor shape1 phs1 v1 a -> Tensor shape2 phs2 v2 a -> Tensor (BroadcastShapes shape1 shape2) (UnionPlaceholder phs1 phs2) Build a
add (Tensor t1) (Tensor t2) = Tensor (t1 `TF.add` t2)

sub :: forall (shape1 :: [Nat]) (shape2 :: [Nat]) (phs1 :: [(Symbol, [Nat], Type)]) (phs2 :: [(Symbol, [Nat], Type)]) v1 v2 a.
       TF.OneOf '[(Complex Double), (Complex Float), Int32, Int64, Word16, Double, Float] a
    => Tensor shape1 phs1 v1 a -> Tensor shape2 phs2 v2 a -> Tensor (BroadcastShapes shape1 shape2) (UnionPlaceholder phs1 phs2) Build a
sub (Tensor t1) (Tensor t2) = Tensor (t1 `TF.sub` t2)

mul :: TF.OneOf '[(Complex Double), (Complex Float), Int16, Int32, Int64, Int8, Word16, Word8, Double, Float] a
    => Tensor shape1 phs1 v1 a -> Tensor shape2 phs2 v2 a -> Tensor (BroadcastShapes shape1 shape2) (UnionPlaceholder phs1 phs2) Build a
mul (Tensor t1) (Tensor t2) = Tensor (t1 `TF.mul` t2)

matMul :: forall (shapeout::[Nat]) (n::Nat) a (i::Nat) (o::Nat) (p::[(Symbol, [Nat], Type)]) (q::[(Symbol, [Nat], Type)]) v'1 v'2.
          (TF.OneOf '[(Complex Double), (Complex Float), Int32, Word16, Double, Float] a,
           shapeout ~ '[i,o])
       => Tensor '[i,n] p v'1 a -> Tensor '[n,o] q v'2 a -> Tensor '[i,o] (UnionPlaceholder p q) Build a
matMul (Tensor t1) (Tensor t2) = Tensor (t1 `TF.matMul` t2)

argMax :: forall (n::Nat) (output_shape::[Nat]) (shape::[Nat]) (phs::[(Symbol,[Nat],Type)]) t output_type v.
          (KnownNat n,
           output_shape ~ RemoveAxisFromShape n shape,
           TF.OneOf '[(Complex Double), (Complex Float), Int16, Int32, Int64, Int8, Word16, Word8, Double, Float] t,
           TF.OneOf '[Int32, Int64] output_type)
       => Proxy n
       -> Tensor shape phs v t
       -> Tensor output_shape phs Build output_type
argMax p (Tensor t) = Tensor $ TF.argMax t (TF.scalar (fromInteger $ natVal p :: Int64))

softmax :: (TF.OneOf '[Word16, Double, Float] t, KnownNat batchSize, KnownNat outs)
        => Tensor '[batchSize, outs] phs v t
        -> Tensor '[batchSize, outs] phs Build t
softmax (Tensor t) = Tensor $ TF.softmax t

scalar :: TF.TensorType a => a -> Tensor '[1] '[] Build a
scalar = Tensor . TF.scalar

oneHot :: (TF.TensorType t,
           KnownNat n,
           TF.OneOf '[Int32, Int64, Word8] tI,
           output_shape ~ AddAxisToEndShape shape n)
       => Proxy n
       -> Tensor '[1] '[] v'2 t
       -> Tensor '[1] '[] v'3 t
       -> Tensor shape phs v'1 tI
       -> Tensor output_shape phs Build t
oneHot p (Tensor t1) (Tensor t2) (Tensor tinput) = Tensor $ TF.oneHot tinput (TF.scalar (fromInteger $ natVal p :: Int32)) t1 t2

-- TODO(helq): search how to use oneHot with both tensors and "numbers" entered by hand,
-- see usage in non-dependent-typed version of example tensorflow mnist
oneHot_ :: (TF.TensorType t,
           KnownNat n,
           TF.OneOf '[Int32, Int64, Word8] tI,
           output_shape ~ AddAxisToEndShape shape n)
        => Proxy n
        -> t
        -> t
        -> Tensor shape phs v'1 tI
        -> Tensor output_shape phs Build t
oneHot_ p t1 t2 = oneHot p (scalar t1) (scalar t2)

reduceMean :: forall (shape::[Nat]) (phs::[(Symbol,[Nat],Type)]) a v.
              TF.OneOf '[Double, Float, Complex Float, Complex Double] a
           => Tensor shape phs v a
           -> Tensor '[1] phs Build a
reduceMean (Tensor t) = Tensor $ TF.reduceMean t

softmaxCrossEntropyWithLogits :: (TF.OneOf '[Word16, Double, Float] a,
                                  addphs ~ UnionPlaceholder phs1 phs2)
                              => Tensor '[batchSize, n] phs1 v'1 a
                              -> Tensor '[batchSize, n] phs2 v'2 a
                              -> (Tensor '[batchSize] addphs Build a, Tensor '[batchSize, n] addphs Build a)
softmaxCrossEntropyWithLogits (Tensor feats) (Tensor labels) = (Tensor loss, Tensor backprop)
  where (loss, backprop) = TF.softmaxCrossEntropyWithLogits feats labels

equal :: (TF.OneOf '[(Complex Double), (Complex Float), Bool, ByteString, Int16, Int32, Int64, Int8, Word16, Word8, Double, Float] a)
      => Tensor shape1 phs1 v'1 a -> Tensor shape2 phs2 v'2 a -> Tensor (BroadcastShapes shape1 shape2) (UnionPlaceholder phs1 phs2) Build Bool
equal (Tensor t1) (Tensor t2) = Tensor (t1 `TF.equal` t2)

truncatedNormal :: forall (shape::[Nat]) a m.
                   (MonadBuild m,
                    TF.OneOf '[Word16, Double, Float] a,
                    KnownNatList shape)
                => m (Tensor shape '[] Value a)
truncatedNormal = Tensor <$> TF.truncatedNormal (TF.vector shape_)
  where shape_ = fmap fromInteger $ natListVal (Proxy :: Proxy shape)

relu :: TF.OneOf '[Int16, Int32, Int64, Int8, Word16, Word8, Double, Float] a
     => Tensor shape phs v a -> Tensor shape phs Build a
relu (Tensor t) = Tensor $ TF.relu t

cast :: (TF.TensorType srcT, TF.TensorType dstT)
     => Tensor shape phs v srcT -> Tensor shape phs Build dstT
cast (Tensor t) = Tensor $ TF.cast t

square :: TF.OneOf '[(Complex Double), (Complex Float), Int32, Int64, Word16, Double, Float] a
    => Tensor shape phs v a -> Tensor shape phs Build a
square (Tensor t) = Tensor (TF.square t)

reshape :: forall t shape1 shape2 phs v.
           (TF.TensorType t, ShapeProduct shape1 ~ ShapeProduct shape2, KnownNatList shape2)
        => Tensor shape1 phs v t
        -> Tensor shape2 phs Build t
reshape (Tensor t) = Tensor $ TF.reshape t $ TF.vector $ map
            (fromInteger :: Integer -> Int32)
            (natListVal (Proxy :: Proxy shape2))

sigmoid :: forall (phs::[(Symbol,[Nat],Type)]) a v s.
           TF.OneOf '[Complex Double, Complex Float, Word16, Double, Float] a
        =>   Tensor s phs v a -> Tensor s phs Build a
sigmoid (Tensor t) = Tensor (TF.sigmoid t)

shape :: forall t shape phs v.
         (TF.TensorType t, KnownNatList shape) -- , TF.OneOf '[Int32, Int64] out_type
       => Tensor shape phs v t
       -> Tensor '[KnownNatListLength shape] phs Build Int32
shape (Tensor t) = Tensor (TF.shape t)

max_issues_repo_name: vilunov/tensorflow-haskell-deptyped


text:
{-# LANGUAGE RankNTypes                               #-}
{-# LANGUAGE RecordWildCards                          #-}
{-# LANGUAGE ScopedTypeVariables                      #-}
{-# LANGUAGE TypeApplications                         #-}
{-# LANGUAGE TypeFamilies                             #-}
{-# LANGUAGE TypeInType                               #-}
{-# LANGUAGE TypeOperators                            #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}

module Tests.Util (
    groupTree
  , CloseEnough(..)
  , generateData
  , withSize
  , dot
  ) where

import           Data.Complex
import           Data.Proxy
import           Data.Type.Equality
import           GHC.TypeLits.Compare
import           GHC.TypeNats
import           Hedgehog
import           Hedgehog.Internal.Property
import           Numeric.Natural
import           Statistics.Transform
import           Test.Tasty
import           Test.Tasty.Hedgehog
import qualified Data.Vector.Sized          as V
import qualified Hedgehog.Gen               as Gen
import qualified Hedgehog.Range             as Range

groupTree :: Group -> TestTree
groupTree Group{..} = testGroup (unGroupName groupName)
                                (map (uncurry go) groupProperties)
  where
    go :: PropertyName -> Property -> TestTree
    go n = testProperty (mkName (unPropertyName n))
    mkName = map deUnderscore . drop (length @[] @Char "prop_")
    deUnderscore '_' = ' '
    deUnderscore c   = c

newtype CloseEnough n = CE { getCE :: V.Vector n Double }
  deriving Show

instance KnownNat n => Eq (CloseEnough n) where
    CE x == CE y = ((d `dot` d) / sqrt ((x `dot` x) * (y `dot` y))) < 0.0001
      where
        d = V.zipWith (-) x y

dot :: Num a => V.Vector n a -> V.Vector n a -> a
dot a b = sum $ V.zipWith (*) a b

withSize
    :: Monad m
    => Range Natural
    -> (forall n. (KnownNat n, 1 <= 2^n) => Proxy n -> PropertyT m a)
    -> PropertyT m a
withSize r f = do
    n <- forAll $ Gen.integral r
    case someNatVal n of
      SomeNat (p :: Proxy n) -> do
        LE Refl <- pure $ Proxy @1 %<=? Proxy @(2^n)
        f p

generateData
    :: KnownNat n
    => Gen (V.Vector (2^n) Double)
generateData = fmap (fmap realPart . ifftSized) . V.generateM $ \i ->
    let i' = recip . (+ 1) . fromIntegral $ i
    in  mkPolar <$> Gen.double (Range.exponentialFloat (i' / 10) i')
                <*> Gen.double (Range.constant (-pi) pi)

ifftSized
    :: V.Vector (2^n) (Complex Double)
    -> V.Vector (2^n) (Complex Double)
ifftSized = V.withVectorUnsafe ifft

max_stars_repo_name: mstksg/emd

max_issues_repo_name: mstksg/emd

max_forks_repo_name: mstksg/emd


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

-----------------------------------------------------------------------------
-- |
-- Module      :  Graphics.Rendering.Plot.Render.Plot.Data
-- Copyright   :  (c) A. V. H. McPhail 2010, 2013
-- License     :  BSD3
--
-- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
-- Stability   :  provisional
-- Portability :  portable
--
-- Rendering 'Figure's
--
-----------------------------------------------------------------------------

module Graphics.Rendering.Plot.Render.Plot.Data (
                                       -- * Rendering
                                       renderData
                                       ) where

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

import Data.List(partition)
--import Prelude.Unicode

--import Foreign.Storable
--import Foreign.Ptr

--import Data.Packed.Vector
--import Data.Packed.Matrix
--import Data.Packed()
import Numeric.LinearAlgebra

import qualified Data.Array.IArray as A
--import qualified Data.Array.MArray as M
import qualified Data.Array.Base as B

import Data.Word

import Data.Maybe
import qualified Graphics.Rendering.Cairo as C
import qualified Graphics.Rendering.Cairo.Matrix as CM

import Control.Monad.Reader
import Control.Monad.State
import Control.Monad.Maybe

import Graphics.Rendering.Plot.Types

import Graphics.Rendering.Plot.Render.Types
import Graphics.Rendering.Plot.Render.Plot.Format
import Graphics.Rendering.Plot.Render.Plot.Glyph
--import Graphics.Rendering.Plot.Render.Plot.Annotation
import Graphics.Rendering.Plot.Defaults

import Prelude hiding(min,max,abs)
import qualified Prelude

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

findMinIdx, findMaxIdx :: Vector Double -> Double -> Int -> Int -> Int
findMinIdx v x n max
    | n >= max       = error "findMinIdx: data not in range"
    | v @> n >= x    = n
    | otherwise     = findMinIdx v x (n+1) max

findMaxIdx v x n min
    | n < 0          = error "findMaxIdx: data not in range"
    | v @> n <= x     = n
    | otherwise      = findMaxIdx v x (n-1) min

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

greySurfaceFromMatrix :: C.SurfaceData Int Word8 -> Surface -> Int -> Int -> Int -> IO ()
greySurfaceFromMatrix s m stride r' c' = do
  let r = rows m
      c = cols m
  let fm = flatten m
      mx = maxElement m
      mn = minElement m
  mapM_ (\ri -> mapM_ (\(rj,ci) -> do
     let mi = ((rj `div` r')*c) + (ci `div` c')
     let e = round . (* 255) . (/ (mx-mn)) . (\x -> x - mn) $ (fm @> mi)
     let si = (rj*stride) + ci
     B.unsafeWrite s si e) $ zip (repeat ri) [0..((c*c')-1)]) [0..((r*r')-1)]

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

zeroToOne :: Double -> Double
zeroToOne x
    | x == 0.0  = 1.0
    | otherwise = x

renderData :: Ranges -> BarSetting -> SampleData -> DataSeries -> Render ()
renderData _ _ _ (DS_Surf m) = do
  (BoundingBox x y w h) <- get
  let r = rows m
      c = cols m
  cairo $ do
    C.save
    --C.setAntialias C.AntialiasNone
    let r'' = Prelude.min 4 ((round h) `div` r)
        c'' = Prelude.min 4 ((round w) `div` c)
        r' = if r'' < 1 then 1 else r''
        c' = if c'' < 1 then 1 else c''
    s <- liftIO $ C.createImageSurface C.FormatA8 (c*c') (r*r')
    p <- liftIO $ C.imageSurfaceGetPixels s
    C.surfaceFlush s
    stride <- liftIO $ C.imageSurfaceGetStride s
    liftIO $ greySurfaceFromMatrix p m stride r' c'
    C.surfaceMarkDirty s
    C.setSourceSurface s x y
    pa <- C.getSource
    pm <- liftIO $ C.patternGetMatrix pa
    let pm' = CM.scale ((fromIntegral (c*c'))/w) ((fromIntegral (r*r'))/h) pm
    liftIO $ C.patternSetMatrix pa pm'
    --C.patternSetFilter pa C.FilterBest
    C.rectangle x y w h --(fromIntegral c) (fromIntegral r)
    C.paint
    C.stroke
    C.restore
    return ()

renderData r bc sd ds = do
  let aos = case ds of
              (DS_Y         os') -> zip (repeat (AbsFunction id)) (A.elems os')
              (DS_1toN abs' os') -> zip (repeat abs')             (A.elems os')
              (DS_1to1 aos')     -> A.elems aos'
              _                  -> error "renderData: DataSeries not handled"
  let (los,ups) = partition (\(_,DecSeries o _) -> isLower o) aos
  (BoundingBox x y w h) <- get
  let (xsc,xmin',xmax') = getRanges XAxis Lower r
  let (xmin,xmax) = if xsc == Log then (logBase 10 $ zeroToOne xmin',logBase 10 $ zeroToOne xmax') else (xmin',xmax')
  let xscale = w/(xmax-xmin)
  cairo $ C.save
  let (yscl,yminl',ymaxl') = getRanges YAxis Lower r
  let (yminl,ymaxl) = if yscl == Log then (logBase 10 $ zeroToOne yminl',logBase 10 $ zeroToOne ymaxl') else (yminl',ymaxl')
  let yscalel = h/(ymaxl-yminl)
  -- transform to data coordinates
  cairo $ do
    C.translate x (y+h)
    --C.scale xscale yscalel
    C.translate (-xmin*xscale) (yminl*yscalel)
    flipVertical
  los' <- configureBars xsc yscl xmin xmax xscale yscalel bc los
  mapM_ (renderSeries xsc yscl xmin xmax xscale yscalel sd) los'
  cairo $ C.restore
  when (not $ null ups)
           (do
             cairo $ C.save
             let (yscu,yminu',ymaxu') = getRanges YAxis Upper r
             let (yminu,ymaxu) = if yscu == Log then (logBase 10 $ zeroToOne yminu',logBase 10 $ zeroToOne ymaxu') else (yminu',ymaxu')
             let yscaleu = h/(ymaxu-yminu)
             -- transform to data coordinates
             cairo $ do
               C.translate x (y+h)
               --C.scale xscale yscaleu
               C.translate (-xmin*xscale) (yminu*yscaleu)
               flipVertical
             ups' <- configureBars xsc yscu xmin xmax xscale yscaleu bc ups
             mapM_ (renderSeries xsc yscu xmin xmax xscale yscaleu sd) ups'
             cairo $ C.restore)
             -- could filter annotations as well
  return ()

logSeries :: Scale -> Vector Double -> Vector Double
logSeries Log a = logBase 10 $ mapVector zeroToOne a
logSeries _   a = a

midpoints :: (Num (Vector t), Container Vector t)
            => (t1, Vector t) -> (t1, Vector t)
midpoints(mi,v) = let v' = subVector 1 (dim v - 1) v
                      w' = subVector 0 (dim v - 1) v
                  in (mi,(v'+w')/2.0)

logSeriesMinMax :: Scale -> (Vector Double,Vector Double) -> (Vector Double,Vector Double)
logSeriesMinMax Log    (v,w) = (logSeries Log v,logSeries Log w)
logSeriesMinMax Linear x     = x

getBar :: (Integer,(Abscissae,DecoratedSeries))
        -> Maybe (Integer,(Abscissae,DecoratedSeries,BarType))
getBar (ix,(as,DecSeries os ds)) = let d = decorationGetBarType ds
   in case d of
        Just d' -> Just (ix,(as,DecSeries os ds,d'))
        Nothing -> Nothing

getBarWidth :: (Integer,(Abscissae,DecoratedSeries,BarType)) -> Render Double
getBarWidth (_,(_,_,(ColourBar _)))                  = do
  (BarOptions w _ _) <- asks (_baroptions . _renderoptions)
  return w
getBarWidth (_,(_,_,(TypeBar (BarOptions w _ _) _))) = return w

replace :: Integer -> a -> [a] -> [a]
replace n x xs = let (pre,post) = splitAt (fromIntegral n) xs
                 in pre ++ x:(tail post)

shiftAbscissa :: (Integer,(Abscissae,DecoratedSeries,BarType)) -> Double
              -> (Integer,(Abscissae,DecoratedSeries))
shiftAbscissa (j,(AbsFunction f,ds,_)) s  = (j,(AbsFunction ((+) s . f),ds))
shiftAbscissa (j,(AbsPoints mi t,ds,_)) s = (j,(AbsPoints mi (addConstant s t),ds))

replaceBars :: [(Integer,(Abscissae,DecoratedSeries))]
            -> [(Abscissae,DecoratedSeries)]
            -> [(Abscissae,DecoratedSeries)]
replaceBars [] as = as
replaceBars ((j,ds):dss) as = replaceBars dss $ replace j ds as

scanStacked :: Vector Double
            -> (Integer,(Abscissae,DecoratedSeries,BarType))
            -> Vector Double
scanStacked v (_,(_,DecSeries (OrdPoints _ o _) _,_)) = v + (getOrdData o)
scanStacked _ _                                       = error "Data.hs:scanStacked: unreachable code"

convertBarToCandle :: (Vector Double,Vector Double)
                   -> OrdSeries
                   -> OrdSeries
convertBarToCandle (v,w) os =
    case os of
      Plain _                 -> MinMax (v,w) Nothing
      Error _ (Right (el,eh)) -> MinMax (v,w) (Just (el,eh))
      _                       -> error "convertBarToCandle: unreachable code"

mkCandlesFromBars :: (Vector Double,Vector Double)
                  -> (Integer,(Abscissae,DecoratedSeries,BarType))
                  -> (Integer,(Abscissae,DecoratedSeries))
mkCandlesFromBars (v,w) (j,(a,DecSeries (OrdPoints ax o mb_l) _,bt)) =
  (j,(a,DecSeries (OrdPoints ax ordSeries mb_l) (DecCand bt)))
      where
        ordSeries = convertBarToCandle (v,w) o
mkCandlesFromBars _     _                                            = error "Data.hs:mkCandlesFromBars: unreachable code"

getOrdData' :: DecoratedSeries -> Series
getOrdData' (DecSeries (OrdPoints _ os _) _) = getOrdData os
getOrdData' _                                = error "Data.hs:getOrdData': unreachable code"

configureBars :: Scale -> Scale
             -> Double -> Double -> Double -> Double
             -> BarSetting
             -> [(Abscissae,DecoratedSeries)]
             -> Render [(Abscissae,DecoratedSeries)]
configureBars _ _ _ _ xscale _ bs aos = do
   let bars = mapMaybe getBar $ zip [0..] aos
   case bs of
     BarNone   -> return aos
     BarSpread -> do
        widths <- mapM getBarWidth bars
        let half_total = (sum widths) / 2
        let shifts' = init $ scanl (+) 0 widths
            shifts = map ((/ xscale) . flip (-) half_total) shifts'
        let shifted = zipWith shiftAbscissa bars shifts
        let aos' = replaceBars shifted aos
        return aos'
     BarStack  -> do
        let od = (getOrdData' . (\(_,b,_) -> b) . snd . head) bars
        let ln = dim $ od
        let zero = constant 0 ln
        let pairs = pair $ scanl scanStacked zero bars
        let candles = zipWith mkCandlesFromBars pairs bars
        let aos' = replaceBars candles aos
        return aos'
            where pair [] = []
                  pair [_] = error "Data.hs:configureBars:pair: unreachable code"
                  pair [x,y] = [(x,y)]
                  pair (x:y:xs) = (x,y) : pair (y:xs)

renderSeries :: Scale -> Scale
             -> Double -> Double -> Double -> Double -> SampleData
             -> (Abscissae,DecoratedSeries) -> Render ()
renderSeries xsc ysc xmin xmax xscale yscale sd (abs,(DecSeries o d)) = do
  dat  <- case o of
     (OrdFunction _ f _)            -> do
        (BoundingBox _ _ w _) <- get
        let t = logSeries xsc $ linspace (round w) (xmin,xmax)
        return $ Left $ Left ((True,t),logSeries ysc $ mapVector f t)
     (OrdPoints _ (Plain o') _)     -> do
        let t = case abs of
                  AbsFunction f    ->
                    if isHist d
                    then (True,mapVector f $ fromList [0.0..(fromIntegral $ dim o')])
                    else (True,mapVector f $ fromList [1.0..(fromIntegral $ dim o')])
                  AbsPoints mi t'  -> (mi,t')
        return $ Left $ Left ((fst t,logSeries xsc $ snd t),logSeries ysc $ o')
     (OrdPoints _ (Error o' (Left e)) _) -> do
        let t = case abs of
                  AbsFunction f    ->
                    if isHist d
                    then (True,mapVector f $ fromList [0.0..(fromIntegral $ dim o')])
                    else (True,mapVector f $ fromList [1.0..(fromIntegral $ dim o')])
                  AbsPoints mi t'  -> (mi,t')
        let t' = (fst t,logSeries xsc $ snd t)
        return $ Left $ Right $ Left ((t',logSeries ysc $ o'),(t',logSeries ysc $ e))
     (OrdPoints _ (Error o' (Right (l,h))) _) -> do
        let t = case abs of
                  AbsFunction f    ->
                    if isHist d
                    then (True,mapVector f $ fromList [0.0..(fromIntegral $ dim o')])
                    else (True,mapVector f $ fromList [1.0..(fromIntegral $ dim o')])
                  AbsPoints mi t'  -> (mi,t')
        let t' = (fst t,logSeries xsc $ snd t)
        return $ Left $ Right $ Right ((t',logSeries ysc $ o'),(t',logSeries ysc $ l),(t',logSeries ysc $ h))
     (OrdPoints _ (MinMax o' Nothing) _) -> do
        let t = case abs of
                  AbsFunction f    -> (True,mapVector f $ fromList [1.0..(fromIntegral $ dim $ fst o')])
                  AbsPoints mi t'  -> (mi,t')
        let t' = (fst t,logSeries xsc $ snd t)
        return $ Right $ Left (t',logSeriesMinMax ysc $ o')
     (OrdPoints _ (MinMax o' (Just (l,h))) _) -> do
        let t = case abs of
                  AbsFunction f    -> (True,mapVector f $ fromList [1.0..(fromIntegral $ dim l)])
                  AbsPoints mi t'  -> (mi,t')
        let t' = (fst t,logSeries xsc $ snd t)
        return $ Right $ Right ((t',logSeriesMinMax ysc o'),(t',(logSeries ysc l,logSeries ysc h)))
  case d of
    (DecLine lt)   -> do
        formatLineSeries lt
        case dat of
          Left (Left (t',y')) -> do
            renderSamples xscale yscale xmin xmax sd Nothing renderLineSample endLineSample t' y'
          _  -> error "Data.hs renderSeries: cannot have error bars with line type"
    (DecPoint pt)  -> do
        (pz,g) <- formatPointSeries pt
        case dat of
          Left (Left (t',y')) -> do
            renderSamples xscale yscale xmin xmax sd Nothing (renderPointSample pz g) endPointSample t' y'
          Left (Right (Left ((t',y'),(_,e')))) -> do
            renderSamples xscale yscale xmin xmax sd Nothing (renderPointSample pz g) endPointSample t' y'
            renderSamples xscale yscale xmin xmax sd Nothing (renderPointSample pz Bot) endPointSample t' (y'-e')
            renderSamples xscale yscale xmin xmax sd Nothing (renderPointSample pz Top) endPointSample t' (y'+e')
          Left (Right (Right ((t',y'),(_,l),(_,h)))) -> do
            renderSamples xscale yscale xmin xmax sd Nothing (renderPointSample pz g) endPointSample t' y'
            renderSamples xscale yscale xmin xmax sd Nothing (renderPointSample pz Bot) endPointSample t' l
            renderSamples xscale yscale xmin xmax sd Nothing (renderPointSample pz Top) endPointSample t' h
          _  -> error "Data.hs renderSeries: cannot have MinMax data series with point type"
    (DecLinPt lt pt) -> do
       formatLineSeries lt
       (pz,g) <- formatPointSeries pt
       case dat of
         Left (Left (t',y')) -> do
           renderSamples xscale yscale xmin xmax sd Nothing renderLineSample endLineSample t' y'
           renderSamples xscale yscale xmin xmax sd Nothing (renderPointSample pz g) endPointSample t' y'
         Left (Right (Left ((t',y'),(_,e')))) -> do
           renderSamples xscale yscale xmin xmax sd Nothing renderLineSample endLineSample t' y'
           renderSamples xscale yscale xmin xmax sd Nothing (renderPointSample pz g) endPointSample t' y'
           renderSamples xscale yscale xmin xmax sd Nothing (renderPointSample pz Bot) endPointSample t' (y'-e')
           renderSamples xscale yscale xmin xmax sd Nothing (renderPointSample pz Top) endPointSample t' (y'+e')
--               error "Data.hs renderSeries: cannot have single error value with line-points type"
         Left (Right (Right ((t',y'),(_,l),(_,h)))) -> do
           renderSamples xscale yscale xmin xmax sd Nothing renderLineSample endLineSample t' y'
           renderSamples xscale yscale xmin xmax sd Nothing (renderPointSample pz g) endPointSample t' y'
           renderSamples xscale yscale xmin xmax sd Nothing (renderPointSample pz Bot) endPointSample t' l
           renderSamples xscale yscale xmin xmax sd Nothing (renderPointSample pz Top) endPointSample t' h
         _  -> error "Data.hs renderSeries: cannot have MinMax data series with line-point type"
    (DecImpulse lt) -> do
       formatLineSeries lt
       case dat of
         Left (Left (t',y')) -> do
           renderSamples xscale yscale xmin xmax sd Nothing renderImpulseSample endImpulseSample t' y'
         _  -> error "Data.hs renderSeries: cannot have error bars with impulse type"
    (DecStep lt) -> do
        formatLineSeries lt
        case dat of
          Left (Left (t',y')) -> do
            renderSamples xscale yscale xmin xmax sd Nothing renderStepSample endStepSample t' y'
          _  -> error "Data.hs renderSeries: cannot have error bars with step type"
    (DecArea lt) -> do
       formatLineSeries lt
       case dat of
         Left (Left (t',y')) -> do
           let ln = dim $ snd t'
               xmin_ix = findMinIdx (snd t') xmin 0 (ln-1)
               x0 = (snd t') @> xmin_ix
               y0 = y' @> xmin_ix
           renderSamples xscale yscale xmin xmax sd Nothing renderAreaSample (endAreaSample x0 y0) t' y'
         _  -> error "Data.hs renderSeries: cannot have error bars with area type"
    (DecBar bt)   -> do
       (bw,bc,c) <- formatBarSeries bt
       (gw,_) <- formatPointSeries defaultPointType
       case dat of
         Left (Left (t',y')) -> do
           renderSamples xscale yscale xmin xmax sd Nothing (renderBarSample bw bc c) endBarSample t' y'
         Left (Right (Left ((t',y'),(_,e')))) -> do
           renderSamples xscale yscale xmin xmax sd Nothing (renderBarSample bw bc c) endBarSample t' y'
           renderSamples xscale yscale xmin xmax sd Nothing (renderPointSampleUpDown gw) endPointSample t' e'
         Left (Right (Right ((t',y'),(_,l'),(_,h')))) -> do
           renderSamples xscale yscale xmin xmax sd Nothing (renderBarSample bw bc c) endBarSample t' y'
           renderSamples xscale yscale xmin xmax sd Nothing (renderPointSample gw Bot) endPointSample t' l'
           renderSamples xscale yscale xmin xmax sd Nothing (renderPointSample gw Top) endPointSample t' h'
         _  -> error "Data.hs renderSeries: cannot have MinMax data series with bar type"
    (DecHist bt)  -> do
       (bw,bc,c) <- formatBarSeries bt
       (gw,_) <- formatPointSeries defaultPointType
       case dat of
         Left (Left (t',y')) -> do
           let ln = dim $ snd $ t'
               xmin_ix = findMinIdx (snd t') xmin 0 (ln-1)
               rest (m,v) = (m,subVector 1 (dim v - 1) v)
               x0 = (snd t') @> xmin_ix
               y0 = 0
           renderSamples xscale yscale xmin xmax sd (Just $ C.moveTo x0 y0) (renderHistSample bw bc c) endHistSample (rest t') y'
         Left (Right (Left ((t',y'),(_,e')))) -> do
           let ln = dim $ snd $ t'
               xmin_ix = findMinIdx (snd t') xmin 0 (ln-1)
               rest (m,v) = (m,subVector 1 (dim v - 1) v)
               x0 = (snd t') @> xmin_ix
               y0 = 0
           renderSamples xscale yscale xmin xmax sd (Just $ C.moveTo x0 y0) (renderHistSample bw bc c) endHistSample (rest t') y'
           renderSamples xscale yscale xmin xmax sd Nothing (renderPointSampleUpDown gw) endPointSample (midpoints t') e'
         Left (Right (Right ((t',y'),(_,l'),(_,h')))) -> do
           let ln = dim $ snd $ t'
               xmin_ix = findMinIdx (snd t') xmin 0 (ln-1)
               rest (m,v) = (m,subVector 1 (dim v - 1) v)
               x0 = (snd t') @> xmin_ix
               y0 = 0
           renderSamples xscale yscale xmin xmax sd (Just $ C.moveTo x0 y0) (renderHistSample bw bc c) endHistSample (rest t') y'
           renderSamples xscale yscale xmin xmax sd Nothing (renderPointSample gw Bot) endPointSample (midpoints t') l'
           renderSamples xscale yscale xmin xmax sd Nothing (renderPointSample gw Top) endPointSample (midpoints t') h'
         _  -> error "Data.hs renderSeries: cannot have MinMax data series with histogram type"
    (DecCand bt)  → do
       (bw,bc,c) <- formatBarSeries bt
       case dat of
         Left _ -> error "Candles series requires two data series (MinMax series)"
         Right (Left (t',y')) -> do
           renderMinMaxSamples xscale yscale xmin xmax sd Nothing (renderCandleSample bw bc c) endCandleSample t' y'
         Right (Right ((t',y'),(_,e'))) -> do
           renderMinMaxSamples xscale yscale xmin xmax sd Nothing (renderWhiskerSample bw bc c False) endWhiskerSample t' e'
           renderMinMaxSamples xscale yscale xmin xmax sd Nothing (renderCandleSample bw bc c) endCandleSample t' y'
    (DecWhisk bt)  → do
      (bw,bc,c) <- formatBarSeries bt
      case dat of
        Left _ -> error "Candles series requires two data series (MinMax series)"
        Right (Left (t',y')) -> do
          renderMinMaxSamples xscale yscale xmin xmax sd Nothing (renderCandleSample bw bc c) endCandleSample t' y'
        Right (Right ((t',y'),(_,e'))) -> do
          renderMinMaxSamples xscale yscale xmin xmax sd Nothing (renderWhiskerSample bw bc c True) endWhiskerSample t' e'
          renderMinMaxSamples xscale yscale xmin xmax sd Nothing (renderCandleSample bw bc c) endCandleSample t' y'
  return ()

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

renderSamples :: Double -> Double
              -> Double -> Double
              -> SampleData
              -> Maybe (C.Render ())
              -> (Double -> Double -> Double -> Double -> C.Render ()) -> (Double -> Double -> C.Render ())
              -> (Bool,Vector Double) -> Vector Double -> Render ()
renderSamples xscale yscale xmin xmax sd s f e (mono,t) y = do
  (BoundingBox _ _ w _) <- get
  let ln = dim t
      (xmin_ix,xmax_ix,num_pts) = if mono
                                  then (findMinIdx t xmin 0 (ln-1)
                                       ,findMaxIdx t xmax (ln-1) 0
                                       ,xmax_ix - xmin_ix + 1)
                                  else (0,ln-1,ln)
      diff'' = floor $ (fromIntegral num_pts)/w
      diff' = if diff'' <= 1 then 1 else diff''
      diff = if mono then if sd then diff' else 1 else 1
  cairo $ do
    case s of
      Nothing -> C.moveTo ((t @> xmin_ix)*xscale) ((y @> xmin_ix)*yscale)
      Just s' -> s'
    _ <- runMaybeT $ do
            mapVectorWithIndexM_ (\j y' -> do
               when (j >= xmin_ix && j `mod` diff == 0)
                        (renderSample j xmax_ix t (f xscale yscale) y')
               return ()) y
    (e xscale yscale)

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

renderMinMaxSamples :: Double -> Double
              -> Double -> Double
              -> SampleData
              -> Maybe (C.Render ())
              -> (Double -> Double -> Double -> (Double,Double) -> C.Render ()) -> (Double -> Double -> C.Render ())
              -> (Bool,Vector Double) -> (Vector Double,Vector Double) -> Render ()
renderMinMaxSamples xscale yscale xmin xmax sd s f e (mono,t) y = do
  (BoundingBox _ _ w _) <- get
  let ln = dim t
      (xmin_ix,xmax_ix,num_pts) = if mono
                                  then (findMinIdx t xmin 0 (ln-1)
                                       ,findMaxIdx t xmax (ln-1) 0
                                       ,xmax_ix - xmin_ix + 1)
                                  else (0,ln-1,ln)
      diff'' = floor $ (fromIntegral num_pts)/w
      diff' = if diff'' <= 1 then 1 else diff''
      diff = if mono then if sd then diff' else 1 else 1
  cairo $ do
    case s of
      Nothing -> C.moveTo ((t @> xmin_ix)*xscale) (((fst $ y) @> xmin_ix)*yscale)
      Just s' -> s'
    _ <- runMaybeT $ mapVectorWithIndexM_ (\j t' -> do
        when (j >= xmin_ix && j `mod` diff == 0)
           (renderMinMaxSample j xmax_ix t' (f xscale yscale) (e xscale yscale) y)
        return ()) t
    return ()

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

renderSample :: Int -> Int -> Vector Double
             -> (Double -> Double -> C.Render ())
             -> Double -> MaybeT C.Render ()
renderSample ix xmax_ix t f y
    | ix >= xmax_ix            = do
                                lift $ f (t @> ix) y
                                fail "end of bounded area"
    | otherwise               = do
                                lift $ f (t @> ix) y

renderMinMaxSample :: Int -> Int -> Double
             -> (Double -> (Double,Double) -> C.Render ()) -> C.Render ()
             -> (Vector Double,Vector Double) -> MaybeT C.Render ()
renderMinMaxSample ix xmax_ix t f e (yl,yu)
    | ix >= xmax_ix            = do
                                lift $ do
                                       f t (yl @> ix,yu @> ix)
                                       e
                                fail "end of bounded area"
    | otherwise               = do
                                lift $ f t (yl @> ix,yu @> ix)

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

renderLineSample :: Double → Double -> Double -> Double -> C.Render ()
renderLineSample xscale yscale x y = C.lineTo (x*xscale) (y*yscale)

endLineSample :: Double -> Double -> C.Render ()
endLineSample _ _ = C.stroke

renderPointSample :: LineWidth -> Glyph -> Double -> Double -> Double -> Double -> C.Render ()
renderPointSample pz g xscale yscale x y = do
  C.moveTo (x*xscale) (y*yscale)
  renderGlyph pz g

renderPointSampleUpDown :: LineWidth -> Double -> Double -> Double -> Double -> C.Render ()
renderPointSampleUpDown pz xscale yscale x y = do
  C.moveTo (x*xscale) (y*yscale)
  renderGlyph pz (if y < 0 then Bot else Top)

endPointSample :: Double -> Double -> C.Render ()
endPointSample _ _ = return ()

renderImpulseSample :: Double -> Double -> Double -> Double -> C.Render ()
renderImpulseSample xscale yscale x y = do
  let x' = x*xscale
  C.moveTo x' 0
  C.lineTo x (y*yscale)
  C.stroke

endImpulseSample :: Double -> Double -> C.Render ()
endImpulseSample _ _ = return ()

renderStepSample :: Double -> Double -> Double -> Double -> C.Render ()
renderStepSample xscale yscale x y = do
  (x',_) <- C.getCurrentPoint
  let y' = y*yscale
  C.lineTo x'         y'
  C.lineTo (x*xscale) y'

endStepSample :: Double -> Double -> C.Render ()
endStepSample _ _ = C.stroke

renderAreaSample :: Double -> Double -> Double -> Double -> C.Render ()
renderAreaSample xscale yscale x y = C.lineTo (x*xscale) (y*yscale)

endAreaSample :: Double -> Double -> Double -> Double -> C.Render ()
endAreaSample x0 _ xscale _ = do
  (x',_) <- C.getCurrentPoint
  let x0' = x0*xscale
  C.lineTo x'  0
  C.lineTo x0' 0
  -- C.lineTo x0 y0
  C.closePath
  C.fill
  C.stroke

renderBarSample :: Width -> Color -> Color -> Double -> Double -> Double -> Double -> C.Render ()
renderBarSample bw c bc xscale yscale x y = do
  setColour bc
  C.rectangle ((x*xscale)-bw/2) 0 bw (y*yscale)
  C.strokePreserve
  setColour c
  C.fill

endBarSample :: Double -> Double -> C.Render ()
endBarSample _ _ = return ()

renderHistSample :: Width -> Color -> Color -> Double -> Double -> Double -> Double -> C.Render ()
renderHistSample _ c bc xscale yscale x y = do
  (x',_) <- C.getCurrentPoint
  let x'' = x*xscale
  C.stroke
  setColour bc
  C.rectangle x' 0 (x''-x') (y*yscale)
  C.strokePreserve
  setColour c
  C.fill
  C.moveTo x'' 0

endHistSample :: Double -> Double -> C.Render ()
endHistSample _ _ = return ()

renderCandleSample :: Width -> Color -> Color -> Double -> Double -> Double -> (Double,Double) -> C.Render ()
renderCandleSample bw c bc xscale yscale x (yl,yu) = do
  setColour bc
  let (yl',yu') = (yl*yscale,yu*yscale)
  C.rectangle ((x*xscale)-bw/2) yl' bw (yu'-yl')
  C.strokePreserve
  if (yl < yu)
     then do
       setColour c
     else do
       C.setSourceRGBA 1 1 1 1
  C.fill

endCandleSample :: Double -> Double -> C.Render ()
endCandleSample _ _ = return ()

renderWhiskerSample :: Width -> Color -> Color → Bool -> Double -> Double -> Double -> (Double,Double) -> C.Render ()
renderWhiskerSample bw _ bc whiskers xscale yscale x (yl,yu) = do
  setColour bc
  let (x',yl',yu') = (x*xscale,yl*yscale,yu*yscale)
  C.moveTo x' yl'
  C.lineTo x' yu'
  if whiskers
     then do
       C.moveTo (x'-bw/2) yu'
       C.lineTo (x'+bw/2) yu'
       C.moveTo (x'-bw/2) yl'
       C.lineTo (x'+bw/2) yl'
     else return ()
  C.stroke

endWhiskerSample :: Double -> Double -> C.Render ()
endWhiskerSample _ _ = return ()

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

max_stars_repo_name: tsurucapital/plot

max_issues_repo_name: tsurucapital/plot

max_forks_repo_name: tsurucapital/plot


text:
module Proginduction where

import MHMonad
import Distr
import Data.Monoid
import Statistics.Distribution
import Statistics.Distribution.Normal (normalDistr)
import Graphics.Rendering.Chart
import Graphics.Rendering.Chart.Backend.Diagrams
import Graphics.Rendering.Chart.State
import Data.Colour
import Data.Colour.Names
import Data.Default.Class
import Control.Lens

-- A grammar for a simple expression language
data BinOp = Add | Multiply
data UnOp = Neg
data Expr = BinOp BinOp Expr Expr
          | UnOp UnOp Expr
          | Number Double
          | Var
          | Cond BoolExpr Expr Expr
data BoolExpr = Leq Expr Expr

-- Simple evaluation function for expressions
eval :: Expr -> Double -> Double
eval (BinOp Add e1 e2) x = (eval e1 x) + (eval e2 x)
eval (BinOp Multiply e1 e2) x = (eval e1 x) * (eval e2 x)
eval (UnOp Neg e ) x = - (eval e x)
eval (Number r) x = r
eval (Cond (Leq e1 e2) e3 e4) x = if (eval e1 x) <= (eval e2 x) then (eval e3 x) else (eval e4 x)
eval Var x = x

-- a probabilistic context free grammar
pcfgbinop :: Meas BinOp
pcfgbinop = do n <- categorical [0.5,0.5]
               case n of 0 -> return Add
                         1 -> return Multiply

pcfgexpr :: Meas Expr
pcfgexpr = do n <- categorical [0.2,0.1,0.3,0.39,0.01]
              case n of 0 -> do {o <- pcfgbinop ; e1 <- pcfgexpr ; e2 <- pcfgexpr ; return $ BinOp o e1 e2 }
                        1 -> do {e <- pcfgexpr ; return $ UnOp Neg e }
                        2 -> do {r <- normal 0 3 ; return $ Number r }
                        3 -> return $ Var
                        4 -> do {e1 <- return Var ; r <- normal 0 4 ; e3 <- pcfgexpr ; e4 <- pcfgexpr ; return $ Cond (Leq e1 (Number r)) e3 e4 }

-- length of an expression
len :: Expr -> Integer
len (BinOp _ e1 e2) = len e1 + len e2 + 1
len (UnOp _ e1) = len e1 + 1
len (Number r) = 1
len Var = 1
len (Cond (Leq e1 e2) e3 e4) = len e1 + len e2 + len e3 + len e4 + 1

-- a random expression is drawn from the context free grammar,
-- but then we weight to get a more reasonable length
randexpr :: Meas (Double -> Double)
randexpr = do { e <- pcfgexpr ;
                score ((normalPdf 25 5) (fromIntegral $ len e)) ;
                return $ \x -> eval e x }

-- functions for plotting a histogram
totalweight :: [(Integer,Double)] -> Double
totalweight rxs = sum $ map (\(_,r) -> r) rxs

toHistogram :: [(Integer,Double)] -> Integer -> Integer -> [Double]
toHistogram rxs from to = if from > to then [] else
  ((sum $ map (\(x,r) -> if x == from then r else 0) rxs) / (totalweight rxs)) : (toHistogram rxs (from + 1) to)

chart :: [Double] -> Renderable ()
chart rs = toRenderable layout
 where
  layout =
        layout_title .~ "Sample Bars" ++ btitle
      $ layout_title_style . font_size .~ 10
      $ layout_x_axis . laxis_generate .~ autoIndexAxis alabels
      $ layout_y_axis . laxis_override .~ axisGridHide
      $ layout_left_axis_visibility . axis_show_ticks .~ False
      $ layout_plots .~ [ plotBars bars2 ]
      $ def :: Layout PlotIndex Double

  bars2 = plot_bars_titles .~ [""]
      $ plot_bars_values .~ addIndexes (map (\x -> [x]) rs) -- [45,30],[30,20],[70,25]]
--      $ plot_bars_style .~ BarsClustered
      $ plot_bars_spacing .~ BarsFixGap 30 5
      $ plot_bars_item_styles .~ map mkstyle (cycle defaultColorSeq)
      $ def

  alabels = map show [0 .. ((length rs)-1)]

  btitle = ""
  bstyle = Just (solidLine 1.0 $ opaque black)
  mkstyle c = (solidFillStyle c, bstyle)

-- plot a histogram of lengths when the weighted pcfg grammar is used
test1 =
  let test =
        do
          e <- pcfgexpr
          -- score ((density $ normalDistr 25 5) (fromIntegral $ len e))
          return $ len e
  in
  do
    rxs' <- weightedsamples test
    let rxs = map (\(x,r) -> (x,r)) rxs'
    _ <- renderableToFile def "histogram.svg" (chart (toHistogram (take 100000 rxs) 0 40))
    return ()

sample_fun f =
  [ (x, f x) | x <- [(-0.25),(-0.25+0.1)..6.2]]

plot_funs :: String -> [(Double,Double)] -> [Double -> Double] -> IO ()
plot_funs filename dataset funs =
  let graphs  = map sample_fun funs                 in
  let my_lines  = plot_lines_style . line_color .~ blue `withOpacity` 0.1
                $ plot_lines_values .~ graphs $ def in
  let my_dots = plot_points_style .~ filledCircles 4 (opaque black)
              $ plot_points_values .~ dataset
              $ def in
  let my_layout = layout_plots .~ [toPlot my_lines , toPlot my_dots]
                $ layout_x_axis .
                  laxis_generate .~ scaledAxis def (0,6)
                $ layout_y_axis . laxis_generate .~ scaledAxis def (-2,10)
                $ def in
  let graphic =  toRenderable my_layout in
  do
     putStr ("Generating " ++ filename ++ "...")
     renderableToFile def filename graphic;
     putStrLn (" Done!")
     return ()

dataset :: [(Double, Double)]
dataset = [(0,0.6), (1, 0.7), (2,1.2), (3,3.2), (4,6.8), (5, 8.2), (6,8.4)]

fit prior dataset = do f <- prior
                       mapM (\(x,y) -> score $ (density $ normalDistr (f x) 0.3) y) dataset
                       return f

every n xs = case drop (n-1) xs of
              (y:ys) -> y : every n ys
              [] -> []

-- Plot some fitted expressions
test3 =
  do
    fs' <- mh $ fit randexpr dataset
    let fs = map (\(f,_)->f) $ take 10 $ every 1000 $ drop 1000000 $ fs'
    plot_funs "prog-ind.svg" dataset fs
    return ()

-- Histogram of the posterior length
test4 =
  do
    let posteriorLength =
          do e <- pcfgexpr
             mapM (\(x,y) -> score $ (density $ normalDistr (eval e x) 0.3) y) dataset
             return (len e)
    ls <- mh $ posteriorLength
    let rxs = map (\(x,r) -> (x,getProduct r)) ls
    _ <- renderableToFile def "histogram.svg" (chart (toHistogram (take 100000 $ drop 100000 rxs) 0 40))
    return ()

max_stars_repo_name: PiotrJander/oplss-staton

max_issues_repo_name: PiotrJander/oplss-staton

max_forks_repo_name: PiotrJander/oplss-staton


text:
{-# LANGUAGE OverloadedStrings, ExistentialQuantification, FlexibleContexts, TypeApplications, RecordWildCards #-}

module Main where
--monad-bayes
import Control.Monad
import Control.Monad.Bayes.Sampler
import Control.Monad.Bayes.Traced
import Control.Monad.Bayes.Weighted
import Control.Monad.Bayes.Class
--hvega
import qualified Graphics.Vega.VegaLite as VL
import qualified Graphics.Vega.Tutorials.VegaLite as VT
--data
import qualified Data.Histogram as DH
import Data.Text (Text, pack)
import Data.List (sort, partition)
import Data.Vector.Unboxed.Base (Unbox)
import qualified Data.Vector.Unboxed as Vec
--numeric
import Numeric.Log
import Numeric.LinearAlgebra (Matrix, Vector, toList, toInt, vector, matrix, dot, (#>), (!), cmap, scalar, size)
--src
import BNN
import PlotVL
--csv
import Control.Exception (IOException)
import qualified Control.Exception as Exception
import qualified Data.Foldable as Foldable
import Control.Applicative
import qualified Data.ByteString.Lazy as BL
import Data.Csv
import qualified Data.Vector as V

--Data Type
data Item =
  Item
    { idNumber :: Double,
      date :: String,
      new_cases :: Double
    }
  deriving (Eq, Show)

--Decoding
instance FromNamedRecord Item where
  parseNamedRecord m =
    Item
      <$> m .: "idNumber"
      <*> m .: "date"
      <*> m .: "new_cases"

decodeItems :: BL.ByteString -> Either String (Vector Item)
decodeItems = fmap snd . Data.Csv.decodeByName

decodeItemsFromFile :: FilePath -> IO (Either String (Vector Item))
decodeItemsFromFile filePath = catchShowIO (BL.readFile filePath) >>= return . either Left decodeItems

catchShowIO :: IO a -> IO (Either String a)
catchShowIO action = fmap Right action `Exception.catch` handleIOException
  where
    handleIOException :: IOException -> IO (Either String a)
    handleIOException = return . Left . show

fillingData :: BL.ByteString ->  [Double] -> [Data]
fillingData csvData noise = case decodeByName csvData of
  Left err -> putStrLn err
  Right (_,v) -> [Data (idNumber v) ((new_cases v) + n) | n <- noise]

main :: IO ()
main = do
    -- putStrLn "Enter the number of samples to train: "
    -- input1 <- getLine
    -- let nsamples = (read input1 :: Int)

    putStrLn "Enter the number of nodes: "
    input3 <- getLine
    let nnodes = (read input3 :: Int)

    putStrLn "Enter the number of predicted samples: "
    input4 <- getLine
    let npredictive = (read input4 :: Int)

    let nsamples = 365
    noise <- sampleIOfixed $ replicateM nsamples $ normal 0.0 0.5

    -- case decodeByName csvData of
    --   Left err -> putStrLn err
    --   Right (_, v) -> decodeItems

    csvData <- BL.readFile "/home/leduin/Desktop/10th Semester/Bayesian Statistics/Project/owid-covid-data-ecuador.csv"

    let observations = fillingData csvData noise

    -- let observations = case decodeByName csvData of
    --   Left err -> putStrLn err
    --   Right (_,v) -> [Data (idNumber v) ((new_cases v) + n) | n <- noise]

-- A SIMPLE NEURAL NETWORK

    --[ Data x (0.5 * x - 2 + n)
    --let observations = take nsamples
    --         [ Data x (2 * sin (x) + 1 + n)
    --         | (x, n) <- zip (map ((/(fromIntegral nsamples)) . fromIntegral ) [0, 10 ..]) noise
    --         ]

        --nnodes = 3
        --mkSampler = prior . mh 60000
        mkSampler = prior . mh npredictive
    predicted <-
        sampleIOfixed $ mkSampler $ predDist $
        postNN (priorNN nnodes) observations

    let hist = histo2D (0, 10, 10) (-10, 20, 10)
                ((\(nn, d) -> (xValue d, yValue d)) <$> predicted)
        cents = Vec.toList $ DH.binsCenters $ DH.bins hist
        val = Vec.toList $ DH.histData hist

    VL.toHtmlFile "final.html" $
        plot
            (600, 300)
            (L [imagePlot "xModel" "yModel" "zModel", scatterBlue "xObs" "yObs" (0, 10) (-10, 10)])
            ( Cols
                [ ("xModel", VL.Numbers (fst <$> cents)),
                ("yModel", VL.Numbers (snd <$> cents)),
                ("zModel", VL.Numbers val),
                ("xObs", VL.Numbers (xValue <$> observations)),
                ("yObs", VL.Numbers (yValue <$> observations))
                ]
            )

    putStrLn "Success! Check html images in directory"

max_stars_repo_name: leduin/Bayesian-Neural-Network

max_issues_repo_name: leduin/Bayesian-Neural-Network

max_forks_repo_name: leduin/Bayesian-Neural-Network


text:
{-# LANGUAGE TemplateHaskell #-}
module NewTTRS.Law where

import Control.Lens
import Data.Array.Unboxed
import Data.List
import Data.Monoid
import Statistics.Distribution (complCumulative, cumulative)
import Statistics.Distribution.Normal (normalDistr)

import NewTTRS.Outcome

data Law = Law { lawRaw :: !(UArray Int Double)
               , lawMean, lawStddev :: !Double }

data LawUpdate = LawUpdate
  { playerLaw , opponentLaw :: Law
  , updateOutcome :: Outcome
  }

-- | Law assigned to unrated players
defaultLaw :: Law
defaultLaw = normalLaw 1800 450

-- | The list of discrete scores characterized by a law
omega :: [Int]
omega = [0,10..3600]

-- | The parameter used to characterize upset probability.
alpha :: Double
alpha = 0.0148540595817432

-- | Generate a normalized law from a list of probabilities
-- which correspond to the elements of 'omega'.
lawFromList :: [Double] -> Law
lawFromList xs = Law (listArray (0,360) normalized) mean (sqrt variance)
  where
  normalized = fmap (/ sum xs) xs
  mean     = sum [ fromIntegral p       * x | (p,x) <- zip omega normalized ]
  variance = sum [ fromIntegral (p * p) * x | (p,x) <- zip omega normalized ]
           - mean * mean

-- | Find the probability that the score is in the range of
-- [i-5,i+5] given a law.
lawAt :: Law -> Int -> Double
lawAt law i = lawRaw law ! (i `div` 10)
{-# INLINE lawAt #-}

-- | Probability of an upset given the difference in two ratings.
upsetProbability :: Int -> Double
upsetProbability d = recip (1 + exp (alpha * fromIntegral d))

-- | Perform a bayesian inference given a player's law, the opponent's law
-- and the outcome of the matches between the two.
lawUpdate ::
  Law {- ^ Player's law to be updated -} ->
  Law {- ^ Opponent's law -} ->
  Outcome {- ^ Player's outcome against opponent -} ->
  Law
lawUpdate lp lq outcome
    = lawFromList
      [ sum [ upsetProbability (q - p) ^ w
            * upsetProbability (p - q) ^ l
            * lawAt lq q
            | q <- omega ]
      * lawAt lp p
      | p <- omega
      ]
  where
  w = view outcomeWins outcome
  l = view outcomeLosses outcome

-- | Perform the inverse computation of 'lawUpdate'
lawUnupdate ::
  Law {- ^ Player's law to be updated -} ->
  Law {- ^ Opponent's law -} ->
  Outcome {- ^ Player's outcome against opponent -} ->
  Law
lawUnupdate lp lq outcome
    = lawFromList
      [ lawAt lp p
      / sum [ upsetProbability (q - p) ^ w
            * upsetProbability (p - q) ^ l
            * lawAt lq q
            | q <- omega ]
      | p <- omega
      ]
  where
  w = view outcomeWins outcome
  l = view outcomeLosses outcome

-- | Chance that a player with the first law will defeat a player with the second law
chanceToWin :: Law -> Law -> Double
chanceToWin lp lq =
  sum [ upsetProbability (q - p)
      * lawAt lp p
      * lawAt lq q
      | q <- omega
      , p <- omega
      ]

-- | Generate a discretized normal law given a mean and standard deviation.
normalLaw ::
   Double {- ^ mean -} ->
   Double {- ^ standard deviation -} ->
   Law
normalLaw mean stddev = lawFromList $ mkNormal 0 3600 mean stddev

-- | Generate a discrete approximation of a normal distribution.
mkNormal ::
  Int    {- ^ lower bound -} ->
  Int    {- ^ upper bound -} ->
  Double {- ^ mean -} ->
  Double {- ^ standard deviation -} ->
  [Double]
mkNormal lo hi mean stddev =
    [ c (lo + 5)  ]
 ++ [ c (p+10) - c p | p <- [lo+5, lo+15 .. hi - 15]]
 ++ [ c' (hi - 5) ]
  where
  distr = normalDistr mean stddev
  c     = cumulative      distr . fromIntegral
  c'    = complCumulative distr . fromIntegral

-- | Degrade a law given a certain number of days. When days is
-- less than 1, no degrading is done.
timeEffect :: Int -> Law -> Law
timeEffect days law
  | days <= 0 = law
  | otherwise = lawFromList
              $ sum [ lawAt law x * timeAt y | x <- omega, y <- [-3600,-3590.. -x]]
              : [ sum [ lawAt law x * timeAt (r - x) | x <- omega ]
                | r <- omega \\ [0,3600] ]
             ++ [ sum [lawAt law x * timeAt y | x <- omega, y <- [3600-x,3610-x..3600]] ]
  where
  timeArray :: UArray Int Double
  timeArray = listArray (-360,360) $ mkNormal (-3600) 3600 0 (70 * sqrt (fromIntegral days / 365))
  timeAt y = timeArray ! (y `div` 10)

-- | Compute a metric of the law used to order players
lawScore :: Law -> Double
lawScore law = lawMean law - 2 * lawStddev law

-- | Return the raw law discrete approximations
lawElems :: Law -> [Double]
lawElems = elems . lawRaw

max_stars_repo_name: glguy/newttrs

max_issues_repo_name: glguy/newttrs

max_forks_repo_name: glguy/newttrs


text:
module Statistics.GModeling.Models.FixedTreeLDA
  (
  ) where

import Statistics.GModeling.DSL

data LDATreeLabels =
    Alpha1 | Alpha2 | Beta
  | TopicDepth | TopicPath | Topic | Doc
  | Symbols | Symbol

ldaTree :: Network LDATreeLabels
ldaTree =
  [
    Only Alpha1         :-> TopicDepth
  , Only Alpha2         :-> TopicPath
  , Only Beta           :-> Symbols
  , (TopicPath :@ Doc)  :-> Topic
  , (TopicDepth :@ Doc) :-> Topic
  , (Symbols :@ Topic)  :-> Symbol
  ]

max_stars_repo_name: nanonaren/gmodeling

max_issues_repo_name: nanonaren/gmodeling

max_forks_repo_name: nanonaren/gmodeling