text:
{-# LANGUAGE ScopedTypeVariables #-}

{-|
Module      : Numeric.LinearAlgebra.Static.Util
Description : Functions to manipulate vectors and matrices
Copyright   : (c) Ryan Orendorff, 2020
License     : BSD3
Stability   : experimental
-}
module Numeric.LinearAlgebra.Static.Util
    ( atM
    , atV
    , foldR
    , foldR1
    )
where

import           GHC.TypeLits
import           Numeric.LinearAlgebra.Static
import qualified Numeric.LinearAlgebra         as LA
import           Data.Fin
import qualified Data.Vector.Storable          as V

---------------------------------------------------------------
-- Get data from vector/matrix with compile time bound check --
---------------------------------------------------------------

-- | Get a value from a vector at a given coordinate/index with guarantee
-- that the element index is within the length of the vector at compile
-- time.
atV
    :: forall n
     . (KnownNat n)
    => R n
    -> Fin n
    -> Double
atV v p = extract v LA.! (fromFin p)

-- | Get element from matrix at a given coordinate/index with guarantee that
-- the element index is within the shape of the matrix at compile time.
atM
    :: forall m n
     . (KnownNat m, KnownNat n)
    => L m n
    -> Fin m
    -> Fin n
    -> Double
atM mat m n = extract mat `LA.atIndex` (fromFin m, fromFin n)

-- | Fold a vector.
foldR :: KnownNat n => (b -> Double -> b) -> b -> R n -> b
foldR f initial x = V.foldl' f initial (extract x)

-- | Fold a vector without an initial element.
foldR1 :: KnownNat n => (Double -> Double -> Double) -> R n -> Double
foldR1 f x = V.foldl1' f (extract x)

max_stars_repo_name: ryanorendorff/convex

max_issues_repo_name: ryanorendorff/convex

max_forks_repo_name: ryanorendorff/convex


text:
{-# LANGUAGE AllowAmbiguousTypes, FlexibleContexts, MagicHash #-}

-- | Utility functions for type-safe PPCA
-- This module contains mostly redefined functions for DimMatrix wrapper
-- around the REPA matrix
module GPLVM.TypeSafe.Util
  ( cholM
  , deleteColumnsM
  , deleteRowsM
  , detSM
  , diagM
  , eigSHM
  , emptyM
  , extendX
  , extendY
  , foldAllSM
  , getColumn
  , getFirstColumns
  , getRow
  , invSM
  , makeLessThenNat
  , mapDiagonalM
  , mapMM
  , meanColumnM
  , meanColumnWithNanM
  , meanCovSM
  , mulM
  , pinvSM
  , solveM
  , substractMeanM
  , sumAllSM
  , sumListMatricesM
  , toListM
  , trace2SM
  , transposeM
  , withMat
  , withVec
  , vecToMat
  , (+^^)
  , (^++^)
  , (-^^)
  , (*^^)
  , (^!^)
  , Decisions
  ) where

import Universum hiding (All, Any, Nat, One, Vector,
                         map, toList, transpose,
                         (%~), (++))

import GPLVM.Types
import GPLVM.TypeSafe.Types
import GPLVM.Util

import Data.Array.Repa
import Numeric.LinearAlgebra.Repa hiding (Matrix, Vector, diag)

import Data.Singletons.Decide (Decision(..))
import Data.Type.Natural hiding (Z)

-- | Creates new Matrix with dimensions in its type and passes it to continuation.
-- Also creates SinI instances for dimensions
withMat
  :: Matrix D Double
  -> (forall (x :: Nat) (y :: Nat). (SingI y, SingI x) => DimMatrix D x y Double -> k)
  -> k
withMat m f =
    let (Z  :.  y  :.  x) = extent m
    in
    case toSing (intToNat y) of
      SomeSing (sy :: Sing m) -> withSingI sy $
        case toSing (intToNat x) of
          SomeSing (sx :: Sing n) -> withSingI sx $ f (DimMatrix @D @m @n m)

-- | Vector cousin of withMat
withVec
  :: Vector D a
  -> (forall (n :: Nat). SingI n => DimVector D n a -> k)
  -> k
withVec vec f =
  case toSing (intToNat n) of
    SomeSing (sy :: Sing m) -> withSingI sy $ f (DimVector @D @m vec)
  where
    (Z :. n) = extent vec

vecToMat
  :: forall (n :: Nat) a.
  SingI n
  => DimVector D n a
  -> DimMatrix D n One a
vecToMat (DimVector vec) =
  DimMatrix $
  fromFunction (Z :. dim :. 1) (\(Z :. a :. b) -> index vec (Z :. a))
  where
    dim = natToInt $ demote @n :: Int

-- | Creates type-level number which is less then given Nat type
makeLessThenNat :: forall (max :: Nat). (SingI max) => Int -> LessThen max
makeLessThenNat i =
  case toSing (intToNat i) of
    (SomeSing (tt:: Sing (uu :: Nat))) ->
          case (tt %< (Sing :: Sing max)) of
            (Proved LS) -> withSingI tt $ Less (Proxy @uu)
            (Disproved _) -> error "the first number is not less then given type-level number"

-- | Creates empty column with desired number of rows
emptyM
  :: forall y x. (SingI y, SingI x) => (x ~ Zero)
  => DimMatrix D y x Double
emptyM =
  let x = fromIntegral $ toNatural (sing :: Sing x)
      y = fromIntegral $ toNatural (sing :: Sing y)
  in DimMatrix (delay  $ fromListUnboxed (Z  :.  y  :.  x) []) :: DimMatrix D y x Double

-- | Returns the column with given index

getColumn
  :: forall n x y. (((n <= x) ~ 'True), SingI n)
  => DimMatrix D y x Double
  -> DimMatrix D y One Double
getColumn (DimMatrix m) =
  let i = (fromIntegral $ toNatural (Sing :: Sing n)) :: Int
  in DimMatrix $ extend (Any   :.  (1 :: Int)) $ slice m (Any :.  i)

-- | Replicates the given column
extendX :: forall n y. (SingI n)
        => DimMatrix D y One Double
        -> DimMatrix D y n Double
extendX (DimMatrix m) =
  let n = natToInt $ fromSing (Sing :: Sing n) :: Int
  in DimMatrix $ (extend (Any :. n) $ slice m (Any   :.  (0 :: Int)))

-- | Replicates the given row
extendY :: forall n x. (SingI n)
        => DimMatrix D One x Double
        -> DimMatrix D n x Double
extendY (DimMatrix m) =
  let n = natToInt $ fromSing (Sing :: Sing n) :: Int
  in DimMatrix $ (extend (Any :. n :. All) $ slice m (Any   :.  (0 :: Int) :. All))

-- | Returns the row with given index
getRow
  :: forall n x y. (((n <= y) ~ 'True), SingI n)
  => DimMatrix D y x Double
  -> DimMatrix D One x Double
getRow (DimMatrix m) =
  let j = (fromIntegral $ toNatural (Sing :: Sing n)) :: Int
  in DimMatrix $ extend (Any   :.  (1 :: Int)  :.  All ) $ slice m (Any   :.  (j :: Int)  :.  All )

-- | Column of the mean values along every dimension
meanColumnM
  :: forall y x. DimMatrix D y x Double
  -> DimMatrix D x One Double
meanColumnM (DimMatrix m) = DimMatrix $ extend (Any   :.  (1 :: Int)) $ meanColumn m

-- | Similar to "meanColumnM" but it does not take into account NaN values
meanColumnWithNanM
  :: forall y x. DimMatrix D y x Double
  -> DimMatrix D y One Double
meanColumnWithNanM (DimMatrix m) = DimMatrix $ extend (Any   :.  (1 :: Int)) $ meanColumnWithNan m

-- todo: make sure that length of deleting indexes is equal to toDel
-- of change the list to list with length in its type
deleteRowsM
  :: forall toDel y1 x1. DimMatrix D y1 x1 Double
  -> [Int]
  -> DimMatrix D (y1 - toDel) x1 Double
deleteRowsM (DimMatrix m) indexes = DimMatrix $
  deleteRows indexes m

deleteColumnsM
  :: forall toDel y1 x1. DimMatrix D y1 x1 Double
  -> [Int]
  -> DimMatrix D y1 (x1 - toDel) Double
deleteColumnsM (DimMatrix m) indexes = DimMatrix $
  deleteColumns indexes m

getFirstColumns
  :: forall toSave y1 x1. (SingI toSave, (toSave <= x1) ~ 'True)
  => DimMatrix D y1 x1 Double
  -> DimMatrix D y1 toSave Double
getFirstColumns (DimMatrix m@(ADelayed (Z :. y :. x) f)) =
  let desiredColumns = natToInt $ fromSing (Sing :: Sing toSave)
  in DimMatrix $ ADelayed (Z :. y :. desiredColumns) f

sumListMatricesM
  :: [DimMatrix D y x Double]
  -> DimMatrix D y x Double
sumListMatricesM [] = error "Empty matrices list"
sumListMatricesM mss = foldl1 (\ms m -> ms +^^ m) mss

-- == Numeric functions from REPA of Numeric.LinearAlgebra for DimMatrix ==

mulM
  :: forall y1 x1 y2 x2.
  (x1 ~ y2)
  => DimMatrix D y1 x1 Double
  -> DimMatrix D y2 x2 Double
  -> DimMatrix D y1 x2 Double
mulM (DimMatrix m1) (DimMatrix m2) = DimMatrix $ delay  $ m1 `mulS` m2

transposeM
  :: DimMatrix D y x Double
  -> DimMatrix D x y Double
transposeM (DimMatrix m) = DimMatrix $ transpose m

mapMM :: (Double -> Double)
  -> DimMatrix D y x Double
  -> DimMatrix D y x Double
mapMM f (DimMatrix m) =  DimMatrix $ map f m

mapDiagonalM :: (Double -> Double)
  -> DimMatrix D y x Double
  -> DimMatrix D y x Double
mapDiagonalM f (DimMatrix m) = DimMatrix $ mapDiagonal f m

invSM
  :: (y ~ x)
  => DimMatrix D y x Double
  -> DimMatrix D y x Double
invSM (DimMatrix m) = DimMatrix $ delay  $ invS m

pinvSM
  :: DimMatrix D y x Double
  -> DimMatrix D x y Double
pinvSM (DimMatrix m) = DimMatrix $ delay $ pinvS m

substractMeanM
  :: DimMatrix D y x Double
  -> DimMatrix D y x Double
substractMeanM (DimMatrix m) = DimMatrix $ substractMean m

(+^^) :: forall y1 x1 y2 x2.
  ( x1 ~ x2
  , y1 ~ y2
  )
  => DimMatrix D y1 x1 Double
  -> DimMatrix D y2 x2 Double
  -> DimMatrix D y2 x2 Double
(+^^) (DimMatrix m1) (DimMatrix m2) = DimMatrix $ m1 +^ m2

(^++^) :: forall y1 x1 y2 x2.
  ( y1 ~ y2
  )
  => DimMatrix D y1 x1 Double
  -> DimMatrix D y2 x2 Double
  -> DimMatrix D y2 (x1 :+: x2) Double
(^++^) (DimMatrix m1) (DimMatrix m2) = DimMatrix $ m1 ++ m2

(-^^) :: forall y1 x1 y2 x2.
  ( x1 ~ x2
  , y1 ~ y2
  )
  => DimMatrix D y1 x1 Double
  -> DimMatrix D y2 x2 Double
  -> DimMatrix D y2 x2 Double
(-^^) (DimMatrix m1) (DimMatrix m2) = DimMatrix $ m1 -^ m2

(*^^) :: forall y1 x1 y2 x2.
  ( x1 ~ x2
  , y1 ~ y2
  )
  => DimMatrix D y1 x1 Double
  -> DimMatrix D y2 x2 Double
  -> DimMatrix D y2 x2 Double
(*^^) (DimMatrix m1) (DimMatrix m2) = DimMatrix $ m1 *^ m2

(^!^) :: forall y x. DimMatrix D y x Double
      -> (LessThen y, LessThen x)
      -> Double
(^!^) (DimMatrix m) ((Less (Proxy :: Proxy y1)), (Less (Proxy :: Proxy x1))) =
  let y1 = natToInt $ fromSing (Sing :: Sing y1)
      x1 = natToInt $ fromSing (Sing :: Sing x1)
  in m ! (Z :. y1 :. x1)

cholM
  :: DimMatrix D y x Double
  -> DimMatrix D y x Double
cholM (DimMatrix m) = DimMatrix $ delay  $ chol $ trustSym $ computeS m

meanCovSM
  :: DimMatrix D y x Double
  -> (DimMatrix D One x Double, DimMatrix D x x Double)
meanCovSM (DimMatrix m) =
  let (meanVec, covMatrix) = meanCovS m
  in ((DimMatrix $ extend (Any  :. (1 :: Int) :. All) meanVec), DimMatrix $ delay $ unSym covMatrix)

solveM
  :: forall x1 y1 x2 y2.
  ( y1 ~ y2
  , x1 ~ y2
  )
  => DimMatrix D y1 x1 Double
  -> DimMatrix D y2 x2 Double
  -> DimMatrix D y1 x2 Double
solveM (DimMatrix m1) (DimMatrix m2) = DimMatrix $ delay $ m1 `solveS` m2

sumAllSM
  :: DimMatrix D y x Double
  -> Double
sumAllSM (DimMatrix m) = sumAllS m

foldAllSM
  :: (Double -> Double -> Double)
  -> Double
  -> DimMatrix D y x Double
  -> Double
foldAllSM f initValue (DimMatrix m) = foldAllS f initValue m

detSM
  :: DimMatrix D y x Double
  -> Double
detSM (DimMatrix m) = detS m

eigSHM
  :: (x ~ y, y1 ~ x) => DimMatrix D y x Double
  -> (DimMatrix D y1 One Double, DimMatrix D y x Double)
eigSHM (DimMatrix m) =
  let hermM = trustSymS m
      (eigenVals, eigenVecs) = eigSH hermM
  in ((DimMatrix $ extend (Any   :.  (1 :: Int)) eigenVals), DimMatrix $ delay eigenVecs)

diagM
  :: forall x y. (y ~ x)
  => DimMatrix D y x Double
  -> DimMatrix D y One Double
diagM (DimMatrix m) = DimMatrix $ diag m

trace2SM
  :: DimMatrix D y x Double
  -> Double
trace2SM (DimMatrix m) = trace2S $ computeS m

toListM
  :: DimMatrix D x1 y1 Double
  -> [Double]
toListM (DimMatrix m) = toList m

-- | this type family is introduced as a syntax-sugar for multiple decisions
type family Decisions k :: Type

type instance Decisions (a, b, c) = (Decision a, Decision b, Decision c)

ext: hs

max_stars_repo_name: serokell/GPLVMHaskell

max_issues_repo_name: serokell/GPLVMHaskell

max_forks_repo_name: serokell/GPLVMHaskell


text:
{-# LANGUAGE FlexibleInstances #-}

module School.Types.FloatEq
( FloatEq(..)
, compareDouble
) where

import Numeric.LinearAlgebra (Element)
import Numeric.LinearAlgebra.Data (Matrix, Vector, cols, toList, toLists)
import School.Train.TrainState (HandlerStore(..))
import School.Unit.UnitActivation (UnitActivation(..))
import School.Unit.UnitGradient (UnitGradient(..))
import School.Unit.UnitParams (UnitParams(..))

class FloatEq a where
  (~=) :: a -> a -> Bool
  (~/) :: a -> a -> Bool
  x ~= y = not (x ~/ y)
  x ~/ y = not (x ~= y)

instance FloatEq Int where
  (~=) = (==)

instance (FloatEq a) => FloatEq (Maybe a) where
  Nothing ~= Nothing = True
  (Just x) ~= (Just y) = x ~= y
  _ ~= _ = False

compareDouble :: Double -> Double -> Double -> Bool
compareDouble prec d1 d2 = if abs d1 < prec
  then abs d2 < prec
  else abs ((d2 - d1) / d1) < prec

doubleEq :: Double -> Double -> Bool
doubleEq = compareDouble 5e-11

instance FloatEq Double where
  (~=) = doubleEq

instance (FloatEq a) => FloatEq [a] where
  l ~= m = (length l == length m)
        && and (zipWith (~=) l m)

mToList :: (Element a) => Matrix a -> [a]
mToList = concat . toLists

instance (Element a, FloatEq a) => FloatEq (Matrix a) where
  m1 ~= m2 = (cols m1 == cols m2)
          && (mToList m1 ~= mToList m2)

instance (Element a, FloatEq a) => FloatEq (Vector a) where
  v1 ~= v2 = toList v1 ~= toList v2

instance (FloatEq b) => FloatEq (Either a b) where
  (Left _) ~= (Left _) = True
  (Right res1) ~= (Right res2) = res1 ~= res2
  _ ~= _ = False

instance (Element a, FloatEq a) => FloatEq (UnitActivation a) where
  (BatchActivation m1) ~= (BatchActivation m2) = m1 ~= m2
  _ ~= _ = False

instance (Element a, FloatEq a) => FloatEq (UnitGradient a) where
  (BatchGradient g1) ~= (BatchGradient g2) = g1 ~= g2
  _ ~= _ = False

instance (Element a, FloatEq a) => FloatEq (UnitParams a) where
  AffineParams { affineBias = b1, affineWeights = w1 } ~= AffineParams { affineBias = b2, affineWeights = w2 } =
    b1 ~= b2 && w1 ~= w2
  EmptyParams ~= EmptyParams = True
  _ ~= _ = False

instance (FloatEq a) => FloatEq (HandlerStore a) where
  (CostList l1) ~= (CostList l2) = l1 ~= l2
  NoStore ~= NoStore = True
  _ ~= _ = False

max_stars_repo_name: jfulseca/School

max_issues_repo_name: jfulseca/School

max_forks_repo_name: jfulseca/School


text:
-- Kolmogorov-Smirnov tests for distribution
--
-- Note that it's not most powerful test for normality.
module MWC.KS (
  tests
  ) where

import qualified Data.Vector.Unboxed as U

import Statistics.Test.KolmogorovSmirnov

import Statistics.Types
import Statistics.Distribution
import Statistics.Distribution.Binomial
import Statistics.Distribution.Exponential
import Statistics.Distribution.Gamma
import Statistics.Distribution.Normal
import Statistics.Distribution.Uniform
import Statistics.Distribution.Beta

import qualified System.Random.MWC               as MWC
import qualified System.Random.MWC.Distributions as MWC

import Test.HUnit hiding (Test)
import Test.Framework
import Test.Framework.Providers.HUnit

tests :: MWC.GenIO -> Test.Framework.Test
tests g = testGroup "Kolmogorov-Smirnov"
  [ testCase "standard"           $ testKS standard           MWC.standard    g
  , testCase "normal m=1 s=2"     $ testKS (normalDistr 1 2) (MWC.normal 1 2) g
    -- Gamma distribution
  , testCase "gamma  k=1   θ=1"   $ testKS (gammaDistr  1   1  ) (MWC.gamma  1   1  ) g
  , testCase "gamma  k=0.3 θ=0.4" $ testKS (gammaDistr  0.3 0.4) (MWC.gamma  0.3 0.4) g
  , testCase "gamma  k=0.3 θ=3"   $ testKS (gammaDistr  0.3 3  ) (MWC.gamma  0.3 3  ) g
  , testCase "gamma  k=3   θ=0.4" $ testKS (gammaDistr  3   0.4) (MWC.gamma  3   0.4) g
  , testCase "gamma  k=3   θ=3"   $ testKS (gammaDistr  3   3  ) (MWC.gamma  3   3  ) g
    -- Uniform
  , testCase "uniform -2 .. 3"    $ testKS (uniformDistr (-2) 3) (MWC.uniformR (-2,3)) g
    -- Exponential
  , testCase "exponential l=1"    $ testKS (exponential 1)       (MWC.exponential 1) g
  , testCase "exponential l=3"    $ testKS (exponential 3)       (MWC.exponential 3) g
    -- Beta
  , testCase "beta a=0.3,b=0.5"    $ testKS (betaDistr 0.3 0.5)       (MWC.beta 0.3 0.5) g
  , testCase "beta a=0.1,b=0.8"    $ testKS (betaDistr 0.3 0.5)       (MWC.beta 0.3 0.5) g
  , testCase "beta a=0.8,b=0.1"    $ testKS (betaDistr 0.3 0.5)       (MWC.beta 0.3 0.5) g
  ]

testKS :: (Distribution d) => d -> (MWC.GenIO -> IO Double) -> MWC.GenIO -> IO ()
testKS distr generator g = do
  sample <- U.replicateM 1000 (generator g)
  case kolmogorovSmirnovTest distr sample of
    Nothing -> assertFailure "Not enough data to make test"
    Just t
      | isSignificant (mkPValue 0.01) t == Significant
        -> assertFailure ("KS test failed: " ++ show t)
      | otherwise
        -> return ()

max_stars_repo_name: Shimuuar/mwc-random

max_issues_repo_name: Shimuuar/mwc-random

max_forks_repo_name: Shimuuar/mwc-random


text:
-- |
-- This module reexports most of the definitions from the \"base\" package,
-- which are meant to be imported unqualified.
--
-- This is adapted from the base-prelude package by Nikita Volkov.
module Prelude.Base
  ( module E
  ) where

-- Reexports
-------------------------

import Control.Applicative as E
import Control.Arrow as E
import Control.Concurrent as E
import Control.Exception as E
import Control.Monad as E
import Control.Monad.Fix as E hiding (fix)
import Control.Monad.ST as E
import Data.Bits as E
import Data.Bool as E
import Data.Char as E
import Data.Complex as E
import Data.Data as E
import Data.Dynamic as E
import Data.Either as E
import Data.Fixed as E
import Data.Foldable as E
import Data.Functor as E
import Data.Function as E
import Data.Int as E
import Data.IORef as E
import Data.Ix as E
import Data.List as E
import Data.Maybe as E
import Data.Monoid as E
import Data.Ord as E
import Data.Ratio as E
import Data.STRef as E
import Data.String as E
import Data.Traversable as E
import Data.Tuple as E
import Data.Unique as E
import Data.Version as E
import Data.Word as E
import Debug.Trace as E
import Foreign.Storable as E (Storable)
import GHC.Conc as E hiding (withMVar, threadWaitWriteSTM, threadWaitWrite, threadWaitReadSTM, threadWaitRead)
import GHC.Generics as E (Generic)
import GHC.IO.Exception as E (IOError,ioError,userError)
import Prelude as E
import System.Environment as E
import System.Exit as E
import System.IO as E (Handle, hClose)
import System.IO.Error as E
import System.Mem as E
import System.Mem.StableName as E
import System.Timeout as E
import Text.ParserCombinators.ReadP as E (ReadP, readP_to_S, readS_to_P)
import Text.ParserCombinators.ReadPrec as E (ReadPrec, readPrec_to_P, readP_to_Prec, readPrec_to_S, readS_to_Prec)
import Text.Read as E (Read(..), readMaybe, readEither)
import Unsafe.Coerce as E

max_stars_repo_name: jackblades/lens-prelude

max_issues_repo_name: jackblades/lens-prelude

max_forks_repo_name: jackblades/lens-prelude


text:
{-# LANGUAGE BangPatterns          #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE KindSignatures        #-}
{-# LANGUAGE LambdaCase            #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE DuplicateRecordFields #-}
{-# LANGUAGE TypeApplications      #-}

module Main where

import Control.Monad.Primitive

import qualified Data.Vector.Sized as VS
--import Data.Vector.Algorithms.Intro as VAI

--import Neldermead
import Data.Vector.GDP.Sorted
import ObjRSimplex
import Numeric.LinearAlgebra.Static

import qualified GDP

simplexEx1 :: VSimplex R 2
simplexEx1 = VS.fromTuple (v1, v2, v3)
  where
    (v1 :: R 2) = vector [0.0, 3.0]
    (v2 :: R 2) = vector [1.0, 4.0]
    (v3 :: R 2) = vector [2.0, 0.0]

createObjSimplex :: (ObjRF n GDP.~~ fObjName)
  -> VSimplex R n
  -> ObjRSimplex fObjName n
createObjSimplex f = ObjRSimplex . VS.map (makeFObjRVtx f)

sortObjSimplex :: forall fObjName n m. PrimMonad m
  => (CmpObjRVtx fObjName n GDP.~~ RVtxObjCmp fObjName)
  -> ObjRSimplex fObjName n
  -> m (SortedObjRSimplex fObjName n)
sortObjSimplex cmp smplx = do
  v <- VS.thaw $ vtxs (smplx :: ObjRSimplex fObjName n)
  SortedObjRSimplex <$> sortByAndFreeze cmp v

main :: IO ()
main = do
  GDP.name norm_1 $ \f -> do
    let smplx = createObjSimplex f simplexEx1
        c = getRVtxObjCmp f
    (SortedObjRSimplex v) <- sortObjSimplex c smplx
    putStrLn . show . VS.map (GDP.the . getObjRVtx) $ GDP.the v

max_stars_repo_name: brnzhg/howod

max_issues_repo_name: brnzhg/howod

max_forks_repo_name: brnzhg/howod


text:
{-# LANGUAGE LambdaCase        #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE RecordWildCards   #-}

module Lib
    ( createBasicDataLayerIO
    , createConfig
    -- * library functions
    , getZendeskResponses
    , getAttachmentsFromComment
    , inspectAttachment
    , ExtractLogFileFunc
    , ExtractErrorCodeFunc
    , processTicket
    , processTicketSafe
    , processTickets
    , fetchTickets
    , showStatistics
    , filterAnalyzedTickets
    , exportZendeskDataToLocalDB
    -- * Optional
    , fetchTicket
    , fetchTicketComments
    , fetchAgents
    , fetchAndShowTickets
    , fetchAndShowTicketsFrom
    , inspectLocalZipAttachment
    , processTicketsFromTime
    ) where

import           Universum

import           UnliftIO.Async (mapConcurrently)

import           Control.Exception.Safe (Handler (..), catches)
import           Data.Attoparsec.Text.Lazy (eitherResult, parse)
import qualified Data.ByteString.Lazy as BS
import           Data.List (nub)
import           Data.Text (stripEnd)

import           System.Directory (createDirectoryIfMissing)
import           System.IO (BufferMode (..), hSetBuffering)

import           Configuration (defaultConfig)
import           DataSource (App, Attachment (..), AttachmentContent (..), Comment (..),
                             Config (..), DBLayer (..), DataLayer (..), DeletedTicket (..),
                             ExportFromTime (..), IOLayer (..), TicketId (..), TicketInfo (..),
                             TicketStatus (..), TicketTag (..), TicketTags (..), User (..),
                             UserId (..), ZendeskResponse (..), asksDBLayer, asksIOLayer,
                             basicDataLayer, connPoolDBLayer, createProdConnectionPool,
                             renderTicketStatus, runApp)

import           Exceptions (ClassifierExceptions (..), ProcessTicketExceptions (..),
                             ZipFileExceptions (..))
import           Http.Layer (basicHTTPNetworkLayer)
import           Http.Queue (createSchedulerConfig, runScheduler)

import           LogAnalysis.Classifier (extractErrorCodes, extractIssuesFromLogs,
                                         prettyFormatAnalysis, prettyFormatLogReadError,
                                         prettyFormatNoIssues, prettyFormatNoLogs)
import           LogAnalysis.Exceptions (LogAnalysisException (..))
import           LogAnalysis.KnowledgeCSVParser (parseKnowLedgeBase)
import           LogAnalysis.Types (Analysis, ErrorCode (..), Knowledge, LogFile, renderErrorCode,
                                    setupAnalysis)
import           Statistics (showStatistics)
import           System.Environment (getEnv)
import           Util (extractLogsFromZip)

------------------------------------------------------------
-- Functions
------------------------------------------------------------

-- | Create configuration.
createConfig :: IO Config
createConfig = do
    putTextLn "Starting Zendesk classifier!"

    createDirectoryIfMissing True "logs"
    hSetBuffering stdout NoBuffering

    putTextLn "Welcome to Zendesk classifier!"

    filePath    <- getEnv "LCPATH"

    token       <- readFile $ filePath <> "/token"       -- Zendesk token
    assignFile  <- readFile $ filePath <> "/assign_to"   -- Select assignee
    knowledges  <- setupKnowledgebaseEnv $ filePath <> "/knowledge.csv"
    assignTo    <- case readEither assignFile of
        Right agentid -> return agentid
        Left  _       -> throwM FailedToReadAssignFile

    connPool    <- createProdConnectionPool

    pure defaultConfig
        { cfgToken         = stripEnd token
        , cfgAssignTo      = assignTo
        , cfgAgentId       = UserId (fromIntegral assignTo)
        , cfgKnowledgebase = knowledges
        , cfgDBLayer       = connPoolDBLayer connPool
        }

  where
    -- Read CSV file and setup knowledge base
    setupKnowledgebaseEnv :: FilePath -> IO [Knowledge]
    setupKnowledgebaseEnv path = do
        kfile <- toLText <$> readFile path
        let kb = parse parseKnowLedgeBase kfile
        case eitherResult kb of
            Left _   -> throwM FailedToParseKnowledgebase
            Right ks -> return ks

-- | Create basic data layer. Do this ONCE!
createBasicDataLayerIO :: Config -> IO (DataLayer App)
createBasicDataLayerIO config = runApp createBasicDataLayer config
  where
    -- | Create a basic @DataLayer@.
    createBasicDataLayer :: App (DataLayer App)
    createBasicDataLayer = do
        -- for now, we limit the rate-limit to 700/minute, which is lower then
        -- 750/minute, but takes into consideration that other things might be called
        -- at the same time as well.
        shedConfig  <- createSchedulerConfig 700

        -- Run the sheduler which resets the number of requests to 0 every minute.
        _           <- runScheduler shedConfig

        -- create the data layer
        pure $ basicDataLayer (basicHTTPNetworkLayer shedConfig)

-- | The function for exporting Zendesk data to our local DB so
-- we can have faster analysis and runs.
-- We expect that the local database exists and has the correct schema.
exportZendeskDataToLocalDB
    :: DataLayer App
    -> ExportFromTime
    -> App [TicketInfo]
exportZendeskDataToLocalDB dataLayer exportFromTime = do

    deleteAllData               <- asksDBLayer dlDeleteAllData

    notDeletedExportedTickets   <- fetchTicketsExportedFromTime dataLayer exportFromTime

    -- Map concurrently if required.
    allTicketData               <-  mapConcurrently fetchTicketData notDeletedExportedTickets

    -- Clear the data.
    deleteAllData

    -- Save all the data.
    _           <- mapM saveTicketDataToLocalDB allTicketData

    -- TODO(ks): We currently don't download the attachment since that
    -- would take a ton of time and would require probably up to 100Gb of space!
    pure notDeletedExportedTickets
  where
    -- | Get ticket comments and create a sum type out of them.
    fetchTicketData :: TicketInfo -> App (TicketInfo,[Comment])
    fetchTicketData ticket = do

        let getTicketComments       = zlGetTicketComments dataLayer

        -- First fetch comments, FAIL-FAST
        ticketComments              <- getTicketComments $ tiId ticket

        pure (ticket, ticketComments)

-- | Fetch all tickets changed from a specific date.
fetchTicketsExportedFromTime :: DataLayer App -> ExportFromTime -> App [TicketInfo]
fetchTicketsExportedFromTime dataLayer exportFromTime = do

    let listDeletedTickets      = zlListDeletedTickets dataLayer
    let exportTickets           = zlExportTickets dataLayer

    -- Yeah, not sure why this happens. Very weird.
    exportedTickets             <- nub <$> exportTickets exportFromTime

    let filteredTicketIds = filterAnalyzedTickets exportedTickets
    let sortedTicketIds   = sortBy compare filteredTicketIds

    -- Yeah, there is a strange behaviour when we export tickets,
    -- we might get deleted tickets as well.
    deletedTickets              <- listDeletedTickets

    pure . sort $ filter (not . isDeletedTicket deletedTickets) sortedTicketIds
  where
    -- | Filter for deleted tickets.
    isDeletedTicket :: [DeletedTicket] -> TicketInfo -> Bool
    isDeletedTicket dTickets TicketInfo{..} =
        tiId `elem` map dtId dTickets

-- | We are missing an intermediate step here to test and separate.
saveTicketDataToLocalDB :: (TicketInfo,[Comment]) -> App ()
saveTicketDataToLocalDB (ticket, ticketComments) = do

    insertTicketInfo            <- asksDBLayer dlInsertTicketInfo
    insertTicketComments        <- asksDBLayer dlInsertTicketComments
    insertCommentAttachments    <- asksDBLayer dlInsertCommentAttachments

    let ticketId = tiId ticket

    insertTicketInfo ticket

    -- For all the ticket comments
    void $ forM ticketComments $ \ticketComment -> do
        insertTicketComments ticketId ticketComment

        -- For all the ticket comment attachments
        forM (cAttachments ticketComment) $ \ticketAttachment ->
            insertCommentAttachments ticketComment ticketAttachment

    pure ()

-- | 'processTicket' with exception handling
processTicketSafe :: DataLayer App -> TicketId -> App ()
processTicketSafe dataLayer tId = catch (void $ processTicket dataLayer tId)
    -- Print and log any exceptions related to process ticket
    -- TODO(ks): Remove IO from here, return the error.
    (\(e :: ProcessTicketExceptions) -> do
        printText <- asksIOLayer iolPrintText
        printText $ "Error on TicketId " <> show (getTicketId tId) <> ":" <> show e
        -- TODO(hs): Implement concurrent logging
        appendF <- asksIOLayer iolAppendFile
        appendF "./logs/errors.log" (show e <> "\n"))

-- | Process ticket with given 'TicketId'
processTicket :: HasCallStack => DataLayer App -> TicketId -> App ZendeskResponse
processTicket dataLayer tId = do

    -- We see 3 HTTP calls here.
    let getTicketInfo       = zlGetTicketInfo dataLayer
    let getTicketComments   = zlGetTicketComments dataLayer
    let postTicketComment   = zlPostTicketComment dataLayer

    mTicketInfo         <- getTicketInfo tId
    comments            <- getTicketComments tId

    let attachments     = getAttachmentsFromComment comments
    case mTicketInfo of
        Nothing         -> throwM $ TicketInfoNotFound tId
        Just ticketInfo -> do
            zendeskResponse <- getZendeskResponses dataLayer comments attachments ticketInfo

            -- post ticket comment, Carl said this is ok
            postTicketComment ticketInfo zendeskResponse

            pure zendeskResponse

-- | When we want to process all possible tickets.
processTickets :: HasCallStack => DataLayer App -> App [ZendeskResponse]
processTickets dataLayer = do

    printText <- asksIOLayer iolPrintText
    printText "Classifier will start processing tickets"

    -- Fetching all the tickets that needs to be processed
    allTickets <- fetchTickets dataLayer

    printText $ "Number of tickets to be analyzed: " <> show (length allTickets)

    zendeskResponses <- mapM (processTicket dataLayer . tiId) allTickets

    putTextLn "All the tickets has been processed."

    pure zendeskResponses

-- | Fetch all tickets that needs to be analyzed
fetchTickets :: DataLayer App -> App [TicketInfo]
fetchTickets dataLayer = do
    allTickets <- zlListToBeAnalysedTickets dataLayer

    let filteredTickets = sortBy compare $ filterAnalyzedTickets allTickets
    -- Any ticket that has a "to_be_analysed" tag
    pure filteredTickets

-- | Fetch a single ticket with given 'TicketId'
fetchTicket :: DataLayer App -> TicketId -> App (Maybe TicketInfo)
fetchTicket dataLayer ticketId = do
    let getTicketInfo       = zlGetTicketInfo dataLayer
    getTicketInfo ticketId

-- | Fetch comments with given 'TicketId'
fetchTicketComments :: DataLayer App -> TicketId -> App [Comment]
fetchTicketComments dataLayer ticketId = do
    let getTicketComments   = zlGetTicketComments dataLayer
    getTicketComments ticketId

-- | A pure function for fetching 'Attachment' from 'Comment'
getAttachmentsFromComment :: [Comment] -> [Attachment]
getAttachmentsFromComment comments = do
    -- Filter tickets without logs
    let commentsWithAttachments :: [Comment]
        commentsWithAttachments = filter commentHasAttachment comments

    -- Filter out ticket without logs
    let attachments :: [Attachment]
        attachments = concatMap cAttachments commentsWithAttachments

    -- Filter out non-logs
    filter isAttachmentZip attachments

  where
    commentHasAttachment :: Comment -> Bool
    commentHasAttachment comment = length (cAttachments comment) > 0

    -- For readability
    isAttachmentZip :: Attachment -> Bool
    isAttachmentZip attachment = aContentType attachment `elem`
        ["application/zip", "application/x-zip-compressed"]

-- | Inspects the comment, attachment, ticket info and create 'ZendeskResponse'
getZendeskResponses :: DataLayer App -> [Comment] -> [Attachment] -> TicketInfo -> App ZendeskResponse
getZendeskResponses dataLayer comments attachments ticketInfo
    | not (null attachments) = inspectAttachments dataLayer ticketInfo attachments
    | not (null comments)    = responseNoLogs ticketInfo
    | otherwise              = throwM $ CommentAndAttachmentNotFound (tiId ticketInfo)
    -- No attachment, no comments means something is wrong with ticket itself
  where
    -- Create 'ZendeskResponse' stating no logs were found on the ticket
    responseNoLogs :: TicketInfo -> App ZendeskResponse
    responseNoLogs TicketInfo{..} = do
        Config {..} <- ask
        pure ZendeskResponse
                { zrTicketId = tiId
                , zrComment  = prettyFormatNoLogs
                , zrTags     = TicketTags [renderTicketStatus NoLogAttached]
                , zrIsPublic = cfgIsCommentPublic
                }

-- | Inspect the latest attachment
inspectAttachments :: DataLayer App -> TicketInfo -> [Attachment] -> App ZendeskResponse
inspectAttachments dataLayer ticketInfo attachments = do

    config          <- ask
    let getAttachment   = zlGetAttachment dataLayer

    lastAttach      <- handleMaybe . safeHead . reverse . sort $ attachments
    att             <- handleMaybe =<< getAttachment lastAttach
    inspectAttachment config ticketInfo att extractLogsFromZip extractIssuesFromLogs
  where
    handleMaybe :: Maybe a -> App a
    handleMaybe Nothing  = throwM $ AttachmentNotFound (tiId ticketInfo)
    handleMaybe (Just a) = return a

-- This is to enable testing on 'inspectAttachment'
type ExtractLogFileFunc = Int -> LByteString -> Either ZipFileExceptions [LogFile]
type ExtractErrorCodeFunc m = [LogFile] -> Analysis -> m Analysis

-- | Given number of file of inspect, knowledgebase and attachment,
-- analyze the logs and return the results.
inspectAttachment :: (MonadCatch m)
                  => Config
                  -> TicketInfo
                  -> AttachmentContent
                  -> ExtractLogFileFunc
                  -> ExtractErrorCodeFunc m
                  -> m ZendeskResponse
inspectAttachment Config{..} ticketInfo attachment extractLogFileFunc extractIssueFunc =
    flip catches [Handler handleZipFileException, Handler handleLogAnalysisException] $ do

        let analysisEnv = setupAnalysis cfgKnowledgebase
        -- let logFiles    = either throwM id $ extractLogFileFunc cfgNumOfLogsToAnalyze
        --                 $ getAttachmentContent attachment
        logFiles <- either throwM pure $
           extractLogFileFunc cfgNumOfLogsToAnalyze $ getAttachmentContent attachment

        -- Perform analysis
        analysisResult <- extractIssueFunc logFiles analysisEnv
        let errorCodes  = extractErrorCodes analysisResult
        let commentRes  = prettyFormatAnalysis analysisResult ticketInfo

        pure $ mkZendeskResponse commentRes errorCodes
   where
     handleZipFileException :: (Monad m) => ZipFileExceptions -> m ZendeskResponse
     handleZipFileException _ =
        pure $ mkZendeskErrorResponse (prettyFormatLogReadError ticketInfo) SentLogCorrupted

     handleLogAnalysisException :: (MonadThrow m) => LogAnalysisException -> m ZendeskResponse
     handleLogAnalysisException = \case
        LogReadException ->
            pure $ mkZendeskErrorResponse (prettyFormatLogReadError ticketInfo) SentLogCorrupted

        NoKnownIssueFound ->
            pure $ mkZendeskErrorResponse (prettyFormatNoIssues ticketInfo) NoKnownIssue

        (JSONDecodeFailure errorText) -> throwM $ JSONDecodeFailure errorText

     mkZendeskResponse :: Text -> [Text] -> ZendeskResponse
     mkZendeskResponse comment errorCodes = ZendeskResponse
        { zrTicketId = tiId ticketInfo
        , zrComment  = comment
        , zrTags     = TicketTags errorCodes
        , zrIsPublic = cfgIsCommentPublic
        }

     mkZendeskErrorResponse :: Text -> ErrorCode -> ZendeskResponse
     mkZendeskErrorResponse comment errorCode = mkZendeskResponse comment [renderErrorCode errorCode]

-- | Filter tickets
filterAnalyzedTickets :: [TicketInfo] -> [TicketInfo]
filterAnalyzedTickets ticketsInfo =
    filter ticketsFilter ticketsInfo
  where
    ticketsFilter :: TicketInfo -> Bool
    ticketsFilter ticketInfo =
           isTicketOpen ticketInfo
        && isTicketBlacklisted ticketInfo
        && isTicketInGoguenTestnet ticketInfo

    unsolvedTicketStatus :: [TicketStatus]
    unsolvedTicketStatus = map TicketStatus ["new", "open", "hold", "pending"]

    isTicketOpen :: TicketInfo -> Bool
    isTicketOpen TicketInfo{..} = tiStatus `elem` unsolvedTicketStatus

    -- | If we have a ticket we are having issues with...
    isTicketBlacklisted :: TicketInfo -> Bool
    isTicketBlacklisted TicketInfo{..} = tiId `notElem` [TicketId 9377, TicketId 10815, TicketId 15066, TicketId 30849]

    isTicketInGoguenTestnet :: TicketInfo -> Bool
    isTicketInGoguenTestnet TicketInfo{..} = "goguen_testnets" `notElem` getTicketTags tiTags

------------------------------------------------------------
-- Utility functions
------------------------------------------------------------

fetchAgents :: DataLayer App -> App [User]
fetchAgents dataLayer = do
    Config{..}          <- ask

    let listAdminAgents = zlListAdminAgents dataLayer
    printText           <- asksIOLayer iolPrintText

    printText "Fetching Zendesk agents"

    agents              <- listAdminAgents

    mapM_ print agents
    pure agents

-- | Fetch tickets that need to be analyzed and print them on console
fetchAndShowTickets :: DataLayer App -> App ()
fetchAndShowTickets dataLayer = do
    allTickets <- fetchTickets dataLayer

    showTickets allTickets

-- | Fetch and show tickets from a specific time.
fetchAndShowTicketsFrom :: DataLayer App -> ExportFromTime -> App ()
fetchAndShowTicketsFrom dataLayer exportFromTime = do
    allTickets <- fetchTicketsExportedFromTime dataLayer exportFromTime

    showTickets allTickets

-- | Display all the TicketId and its tags
showTickets :: [TicketInfo] -> App ()
showTickets tickets = do
    putTextLn $ "There are " <> show (length tickets) <> " tickets."

    forM_ tickets $ \ticket -> do
        let ticketId   = getTicketId $ tiId ticket
        let ticketTags = (sort . getTicketTags . tiTags) ticket
        putTextLn $ "TicketId: " <> show ticketId <> ", Tags: " <> show ticketTags

-- | Inspection of the local zip.
-- This function prints out the analysis result on the console.
-- Can apply same refactoring as 'inspectAttachment'
inspectLocalZipAttachment :: FilePath -> App ()
inspectLocalZipAttachment filePath =
    flip catches [ Handler $ exceptionHandler @ZipFileExceptions
                 , Handler $ exceptionHandler @LogAnalysisException] $ do
        config          <- ask
        printText       <- asksIOLayer iolPrintText

        -- Read the zip file
        fileContent     <- liftIO $ BS.readFile filePath
        logFiles        <- either throwM pure $ extractLogsFromZip 100 fileContent

        let analysisEnv = setupAnalysis $ cfgKnowledgebase config
        analysisResult  <- extractIssuesFromLogs logFiles analysisEnv

        let errorCodes = extractErrorCodes analysisResult

        printText "Analysis result:"
        void $ mapM (printText . show) analysisResult

        printText "Error codes:"
        void $ mapM printText errorCodes
  where
    exceptionHandler :: (Exception e) => e -> App ()
    exceptionHandler = print

-- | When we want to process all tickets from a specific time onwards.
-- Run in parallel.
processTicketsFromTime :: DataLayer App -> ExportFromTime -> App ()
processTicketsFromTime dataLayer exportFromTime = do

    allTickets          <- fetchTicketsExportedFromTime dataLayer exportFromTime

    putTextLn $ "There are " <> show (length allTickets) <> " tickets."

    mZendeskResponses   <- mapConcurrently (processTicket dataLayer . tiId) allTickets

    -- This is single-threaded.
    mapM_ processSingleTicket mZendeskResponses

    putTextLn "All the tickets has been processed."
  where
    -- | Process a single ticket after they were analyzed.
    processSingleTicket :: ZendeskResponse -> App ()
    processSingleTicket zendeskResponse = do

        -- We first fetch the function from the configuration
        printText           <- asksIOLayer iolPrintText
        appendF             <- asksIOLayer iolAppendFile -- We need to remove this.

        let ticketId = getTicketId $ zrTicketId zendeskResponse

        -- Printing id before inspecting the ticket so that when the process stops by the
        -- corrupted log file, we know which id to blacklist.
        printText $ "Analyzing ticket id: " <> show ticketId

        let tags = getTicketTags $ zrTags zendeskResponse
        forM_ tags $ \tag -> do
            let formattedTicketIdAndTag = show ticketId <> " " <> tag
            printText formattedTicketIdAndTag
            appendF "logs/analysis-result.log" (formattedTicketIdAndTag <> "\n")

max_stars_repo_name: input-output-hk/log-classifier

max_issues_repo_name: input-output-hk/log-classifier

max_forks_repo_name: input-output-hk/log-classifier


text:
{-# LANGUAGE BangPatterns        #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeOperators       #-}
{-# LANGUAGE ViewPatterns        #-}
module FokkerPlanck.GPUKernel where

import           Data.Array.Accelerate              as A
import           Data.Array.Accelerate.Data.Complex as A
import qualified Data.Complex                       as C
import Debug.Trace

{-# INLINE moveParticle #-}
moveParticle ::
     (A.RealFloat a, A.Num a, Elt a-- , FromIntegral Int a
     )
  => Exp (a, a, a, a)
  -> Exp (a, a, a, a)
moveParticle (unlift -> (phi, rho, theta, r)) =
  let !x = theta - phi
      !cosX = A.cos x
      !newPhi = phi + (A.atan2 (r * A.sin x) (rho + r * cosX))
      !newRho = A.sqrt $ (rho * rho + r * r + 2 * r * rho * cosX)
      -- !xx = A.round $ newRho * A.cos newPhi :: Exp Int
      -- !y = A.round $ newRho * A.sin newPhi :: Exp Int
      -- !newRho' = A.sqrt . A.fromIntegral $ xx * xx + y * y
      -- !newPhi' = A.atan2 (A.fromIntegral y) (A.fromIntegral xx)
  in lift (newPhi, newRho, theta, r)

{-# INLINE normalizationTerm #-}
normalizationTerm ::
     (A.Floating a, A.Num a, A.RealFloat a, A.Elt (Complex a))
  => Exp a
  -> Exp a
  -> Exp (Complex a)
normalizationTerm halfLogPeriod freq =
  (lift $ 1 A.:+ ((-A.pi) * freq / halfLogPeriod)) /
  (A.exp (lift $ halfLogPeriod A.:+ (-A.pi * freq)) -
   A.exp (lift $ (-halfLogPeriod) A.:+ (A.pi * freq)))

-- {-# INLINE coefficient #-}
-- coefficient ::
--      forall a. (A.Floating a, A.Num a, A.RealFloat a, A.Elt (Complex a), Prelude.Fractional a)
--   => Exp a
--   -> Exp a
--   -> Exp a
--   -> Exp a
--   -> Exp a
--   -> Exp (a, a, a, a)
--   -> Exp (A.Complex a)
-- coefficient halfLogPeriod rFreq thetaFreq rhoFreq phiFreq particle -- (unlift -> (phi, rho, theta, r))
--  =
--   let (phi, rho, theta, r) = unlift particle :: (Exp a, Exp a, Exp a, Exp a)
--      -- (normalizationTerm halfLogPeriod rhoFreq) *
--      -- (normalizationTerm halfLogPeriod rFreq) *
--      -- (lift (A.cos (phiFreq * phi + thetaFreq * (theta - phi)) A.:+ 0)) *
--                            -- (A.cis $
--                            --  (-2 * A.pi) * (rhoFreq * rho + rFreq * (r - rho)) /
--                            --  (A.exp halfLogPeriod))

--   in (lift $
--       ((A.exp $ (-0.5) * (rho + r)) *
--        (A.cos (phiFreq * phi  + thetaFreq * (theta - phi)))
--       ) A.:+
--       0) *
--      (A.cis $ (-1) * ((rhoFreq * rho + rFreq * (r - rho)
--                       ) -- * (2 * A.pi) / (A.exp halfLogPeriod)
--                       -- + phiFreq * phi  + thetaFreq * (theta - phi)
--                      ))
--      -- (lift $ ((A.cos (phiFreq * phi + thetaFreq * (theta - phi) / 2))) A.:+ 0) *
--      -- (A.cis $ (-2 * A.pi) * (rhoFreq * rho + rFreq * (r - rho)) / (A.exp halfLogPeriod))
--      -- ((lift $ rho A.:+ 0) A.** (lift $ (-0.5) A.:+ (rFreq - rhoFreq))) *
--      -- ((lift $ r A.:+ 0) A.** (lift $ (-0.5) A.:+ (-rFreq)))
--      -- (A.cis $ (-A.pi) * (rhoFreq * rho + rFreq * (r - rho)
--      --                        ) / (halfLogPeriod))
--      -- (A.cis $
--      --  (-A.pi) * (rhoFreq * (A.log rho) + rFreq * (A.log (r / rho))) /
--      --  halfLogPeriod)
--   -- in ((lift $ rho A.:+ 0) A.**
--   --     (lift $ 0 A.:+ (-A.pi) * (rhoFreq - rFreq) / halfLogPeriod)) *
--   --    ((lift $ r A.:+ 0) A.** (lift $ 0 A.:+ (-A.pi) * rFreq / halfLogPeriod)) *
--   --    (A.cis $ (-phiFreq) * phi - thetaFreq * (theta - phi))

coefficient ::
     forall a. (A.Floating a, A.Num a, A.RealFloat a, A.Elt (Complex a), Prelude.Fractional a)
  => Exp a
  -> Exp a
  -> Exp a
  -> Exp a
  -> Exp a
  -> Exp (a, a, a, a)
  -> Exp (A.Complex a)
coefficient halfLogPeriod rFreq thetaFreq rhoFreq phiFreq particle =
  let (phi, rho, theta, r) = unlift particle :: (Exp a, Exp a, Exp a, Exp a)
  in (lift $
      ((A.exp $ (-0.5) * (rho + r)) *
       (A.cos (phiFreq * phi + thetaFreq * (theta - phi)))) A.:+
      0) *
     (A.cis $ (-1) * (rhoFreq * rho + rFreq * (r + rho)))

-- {-# INLINE gpuKernel #-}
gpuKernel ::
     forall a. (A.Eq a, A.Floating a, A.Num a, A.RealFloat a, A.Elt (Complex a), A.FromIntegral Int a, Prelude.Fractional a)
  => Exp a
  -> Exp a
  -> Exp (A.Complex a)
  -> Acc (A.Vector (a, a, a, a))
  -> Acc (A.Vector (a, a, a, a))
  -> Acc (A.Vector (A.Complex a))
gpuKernel !maxScaleExp !halfLogPeriodExp !deltaLogRhoComplexExp freqArr particles =
  let -- delta = constant 0.01
      movedParticleArr =
        compute .
        A.map
          (\particle ->
             let (phi, rho, theta, r) =
                   unlift particle :: (Exp a, Exp a, Exp a, Exp a)
                 logRho = A.log rho
                   -- (A.fromIntegral $
                   --  (A.round (((A.log rho) + halfLogPeriodExp) / delta) :: Exp Int)) *
                   -- delta -
                   -- halfLogPeriodExp
                 -- newRho =
                 --   (A.fromIntegral $ (A.round (rho / delta) :: Exp Int)) * delta :: Exp a
             in lift (phi, logRho, theta, (A.log r) :: Exp a)) -- .
        -- afst .
        -- A.filter
        --   (\particle ->
        --      let (_, rho, _, _) =
        --            unlift particle :: (Exp a, Exp a, Exp a, Exp a)
        --      in rho A.> (A.constant 0) ) -- .
        -- A.map moveParticle
        $
        particles
  in A.map
       (\(unlift -> (rFreq, thetaFreq, rhoFreq, phiFreq))
          -- (lift $ delta A.:+ 0) /
          -- (lift $ (8 * A.pi * A.pi * halfLogPeriodExp) A.:+ 0) *
         ->
          (sfoldl
             (\s particle ->
                s +
                (coefficient
                   halfLogPeriodExp
                   rFreq
                   thetaFreq
                   rhoFreq
                   phiFreq
                   particle))
             0
             (constant Z)
             movedParticleArr) * deltaLogRhoComplexExp
          -- (lift $
          --  (16 * A.pi * A.pi * halfLogPeriodExp * halfLogPeriodExp) A.:+ 0)
             -- movedParticleArr
        ) $
     freqArr

{-# INLINE convolveKernel #-}
convolveKernel ::
     (A.Num a, A.RealFloat a, A.Elt (Complex a))
  => Acc (Array DIM4 (Complex a))
  -> Acc (Array DIM4 (Complex a))
  -> Acc (Scalar Int)
  -> Acc (Scalar Int)
  -> Acc (Array DIM4 (Complex a))
  -> Acc (Vector (Complex a))
convolveKernel coefficients harmonics thetaIdx rIdx input =
  let (Z :. numRhoFreq :. numPhiFreq :. cols :. rows) =
        unlift . shape $ input :: (Z :. Exp Int :. Exp Int :. Exp Int :. Exp Int)
      coefficientsArr =
        A.replicate (lift (Z :. All :. All :. cols :. rows)) .
        slice coefficients $
        (lift (Z :. (the rIdx) :. (the thetaIdx) :. All :. All))
      harmonicsArr =
        backpermute
          (shape input)
          (\(unlift -> Z :. rho :. phi :. col :. row :: (Z :. Exp Int :. Exp Int :. Exp Int :. Exp Int)) ->
             lift (Z :. (rho + the rIdx) :. (phi + the thetaIdx) :. col :. row))
          harmonics
  in flatten .
     A.sum .
     A.sum .
     backpermute
       (lift (Z :. cols :. rows :. numRhoFreq :. numPhiFreq))
       (\(unlift -> Z :. col :. row :. rho :. phi :: (Z :. Exp Int :. Exp Int :. Exp Int :. Exp Int)) ->
          lift (Z :. rho :. phi :. col :. row)) .
     A.zipWith (*) harmonicsArr . A.zipWith (*) coefficientsArr $
     input

{-# INLINE coefficient' #-}
coefficient' ::
     forall a.
     ( A.Floating a
     , A.Num a
     , A.RealFloat a
     , A.Elt (Complex a)
     , Prelude.Fractional a
     )
  => Exp a
  -> Exp a
  -> Exp a
  -> Exp a
  -> Exp a
  -> Exp (a, a, a, a, a)
  -> Exp (A.Complex a)
coefficient' sigma rFreq thetaFreq rhoFreq phiFreq particle =
  let (phi, rho, theta, r, v) =
        unlift particle :: (Exp a, Exp a, Exp a, Exp a, Exp a)
  in (lift $ (v * (A.exp $ (sigma - 1) * (rho + r))) A.:+ 0) *
     (A.cis $
      (-1) *
      (rhoFreq * rho + rFreq * (r - rho) + phiFreq * phi +
       thetaFreq * (theta - phi)))

gpuKernel' ::
     forall a.
     ( A.Eq a
     , A.Floating a
     , A.Num a
     , A.RealFloat a
     , A.Elt (Complex a)
     , A.FromIntegral Int a
     , Prelude.Fractional a
     )
  => Exp a
  -> Acc (A.Vector (a, a, a, a))
  -> Acc (A.Vector (a, a, a, a, a))
  -> Acc (A.Vector (A.Complex a))
gpuKernel' sigma freqArr xs =
  A.map
    (\(unlift -> (rFreq, thetaFreq, rhoFreq, phiFreq)) ->
       A.sfoldl
         (\s particle ->
            s + (coefficient' sigma rFreq thetaFreq rhoFreq phiFreq particle))
         0
         (constant Z)
         xs)
    freqArr

{-# INLINE coefficient'' #-}
coefficient'' ::
     forall a.
     ( A.Floating a
     , A.Num a
     , A.RealFloat a
     , A.Elt (Complex a)
     , Prelude.Fractional a
     )
  => Exp a
  -> Exp a
  -> Exp a
  -> Exp a
  -> Exp a
  -> Exp a
  -> Exp (a, a, a, a, a)
  -> Exp (A.Complex a)
coefficient'' sigma period rFreq thetaFreq rhoFreq phiFreq particle =
  let (phi, rho, theta, r, v) =
        unlift particle :: (Exp a, Exp a, Exp a, Exp a, Exp a)
   in lift
        ((v * A.exp ((sigma - 1) * rho) *
          A.cos (phiFreq * phi + thetaFreq * (theta - phi))
         ) :+
         0) *
      A.cis ((-1) * (2 * A.pi / period * (rhoFreq * rho + rFreq * (r - rho))))

gpuKernel'' ::
     forall a.
     ( A.Eq a
     , A.Floating a
     , A.Num a
     , A.RealFloat a
     , A.Elt (Complex a)
     , A.FromIntegral Int a
     , Prelude.Fractional a
     )
  => Exp a
  -> Exp a
  -> Acc (A.Vector (a, a, a, a))
  -> Acc (A.Vector (a, a, a, a, a))
  -> Acc (A.Vector (A.Complex a))
gpuKernel'' sigma period freqArr xs =
  A.map
    (\(unlift -> (rFreq, thetaFreq, rhoFreq, phiFreq)) ->
       A.sfoldl
         (\s particle ->
            s +
            coefficient'' sigma period rFreq thetaFreq rhoFreq phiFreq particle)
         0
         (constant Z)
         xs)
    freqArr

{-# INLINE pinwheelAcc #-}
pinwheelAcc ::
     forall a. (A.Floating a, A.Num a, A.RealFloat a, A.Elt (Complex a), Prelude.Fractional a)
  => Exp a
  -> Exp a
  -> Exp a
  -> Exp a
  -> Exp (a, a, a, a, A.Complex a)
  -> Exp (A.Complex a)
pinwheelAcc rFreq thetaFreq rhoFreq phiFreq particle =
  let (phi, rho, theta, r, v) =
        unlift particle :: (Exp a, Exp a, Exp a, Exp a, Exp (A.Complex a))
  in v * (lift $ ((A.exp $ (-0.5) * (rho + r))) A.:+ 0) *
     (A.cis $
      (-1) *
      (rhoFreq * rho + rFreq * (r - rho) + phiFreq * phi +
       thetaFreq * (theta - phi)))

pinwheelCoefficientsAcc ::
     forall a.
     ( A.Eq a
     , A.Floating a
     , A.Num a
     , A.RealFloat a
     , A.Elt (Complex a)
     , A.FromIntegral Int a
     , Prelude.Fractional a
     )
  => Acc (A.Vector (a, a, a, a))
  -> Acc (A.Vector (a, a, a, a, A.Complex a))
  -> Acc (A.Vector (A.Complex a))
pinwheelCoefficientsAcc freqArr xs =
  A.map
    (\(unlift -> (rFreq, thetaFreq, rhoFreq, phiFreq)) ->
       A.sfoldl
         (\s particle ->
            s + (pinwheelAcc rFreq thetaFreq rhoFreq phiFreq particle))
         0
         (constant Z)
         xs)
    freqArr

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-# LANGUAGE RecordWildCards #-}

import Data.Complex (realPart, imagPart)
import Diagrams.Backend.CmdLine (mainRender)
import Diagrams.Backend.Rasterific.CmdLine ()

import Args (parseArgs, Args(..))
import ExponentialSum (expSums)
import Plot (display)

main :: IO ()
main = do
  Args{..} <- parseArgs
  mainRender options $ display $ select start end $ getPoints coeffs
  where
    getPoints = map toCoords . expSums . map fromIntegral
    select st en = take (en - st) . drop st
    toCoords z = (realPart z, imagPart z)

max_stars_repo_name: mihaimaruseac/esotd

max_issues_repo_name: mihaimaruseac/esotd

max_forks_repo_name: mihaimaruseac/esotd


text:
{-# LANGUAGE AllowAmbiguousTypes, DataKinds, RankNTypes, TypeFamilies,
             TypeOperators #-}
{-# OPTIONS_GHC -Wno-missing-export-lists #-}
{-# OPTIONS_GHC -fconstraint-solver-iterations=16 #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise #-}
-- | Shaped tensor-based implementation of Convolutional Neural Network
-- for classification of MNIST digits. Sports 2 hidden layers.
module HordeAd.Tool.MnistCnnShaped where

import Prelude

import qualified Data.Array.DynamicS as OT
import           Data.Array.Internal (valueOf)
import qualified Data.Array.Shape
import qualified Data.Array.ShapedS as OS
import           Data.Proxy (Proxy)
import qualified Data.Vector.Generic as V
import           GHC.TypeLits (KnownNat, type (+), type (<=), type Div)
import qualified Numeric.LinearAlgebra as HM

-- until stylish-haskell accepts NoStarIsType
import qualified GHC.TypeLits

import HordeAd.Core.DualNumber
import HordeAd.Core.Engine
import HordeAd.Core.PairOfVectors (DualNumberVariables, varS)
import HordeAd.Tool.MnistData

convMnistLayerS
  :: forall kheight_minus_1 kwidth_minus_1 out_channels
            in_height in_width in_channels batch_size d r m.
     ( KnownNat kheight_minus_1, KnownNat kwidth_minus_1, KnownNat out_channels
     , KnownNat in_height, KnownNat in_width
     , KnownNat in_channels, KnownNat batch_size
     , 1 <= kheight_minus_1
     , 1 <= kwidth_minus_1  -- wrongly reported as redundant
     , DualMonad d r m )
  => DualNumber d (OS.Array '[ out_channels, in_channels
                             , kheight_minus_1 + 1, kwidth_minus_1 + 1 ] r)
  -> DualNumber d (OS.Array '[batch_size, in_channels, in_height, in_width] r)
  -> DualNumber d (OS.Array '[out_channels] r)
  -> m (DualNumber d (OS.Array '[ batch_size, out_channels
                                , (in_height + kheight_minus_1) `Div` 2
                                , (in_width + kwidth_minus_1) `Div` 2 ] r))
convMnistLayerS ker x bias = do
  let yConv = conv24 ker x
      replicateBias
        :: DualNumber d (OS.Array '[] r)
        -> DualNumber d (OS.Array '[ in_height + kheight_minus_1
                                   , in_width + kwidth_minus_1 ] r)
      replicateBias = konstS . fromS0
      biasStretched = ravelFromListS
                      $ replicate (valueOf @batch_size)
                      $ mapS replicateBias bias
        -- TODO: this is weakly typed; add and use replicateS instead
        -- or broadcastS or stretchS, possibly with transposeS?
  yRelu <- reluAct $ yConv + biasStretched
  maxPool24 @1 @2 yRelu

convMnistTwoS
  :: forall kheight_minus_1 kwidth_minus_1 num_hidden out_channels
            in_height in_width in_channels batch_size d r m.
     ( KnownNat kheight_minus_1, KnownNat kwidth_minus_1
     , KnownNat num_hidden, KnownNat out_channels
     , KnownNat in_height, KnownNat in_width, KnownNat batch_size
     , in_channels ~ 1
     , 1 <= kheight_minus_1
     , 1 <= kwidth_minus_1
     , DualMonad d r m )
  => OS.Array '[batch_size, in_channels, in_height, in_width] r
  -- All below is the type of all paramters of this nn. The same is reflected
  -- in the length function below and read from variables further down.
  -> DualNumber d (OS.Array '[ out_channels, in_channels
                             , kheight_minus_1 + 1, kwidth_minus_1 + 1 ] r)
  -> DualNumber d (OS.Array '[out_channels] r)
  -> DualNumber d (OS.Array '[ out_channels, out_channels
                             , kheight_minus_1 + 1, kwidth_minus_1 + 1 ] r)
  -> DualNumber d (OS.Array '[out_channels] r)
  -> DualNumber d (OS.Array '[ num_hidden
                             , out_channels
                                 GHC.TypeLits.*
                                   (((in_height + kheight_minus_1) `Div` 2
                                     + kheight_minus_1) `Div` 2)
                                 GHC.TypeLits.*
                                   (((in_width + kwidth_minus_1) `Div` 2
                                     + kwidth_minus_1) `Div` 2)
                             ] r)
  -> DualNumber d (OS.Array '[num_hidden] r)
  -> DualNumber d (OS.Array '[SizeMnistLabel, num_hidden] r)
  -> DualNumber d (OS.Array '[SizeMnistLabel] r)
  -> m (DualNumber d (OS.Array '[SizeMnistLabel, batch_size] r))
convMnistTwoS x ker1 bias1 ker2 bias2
              weigthsDense biasesDense weigthsReadout biasesReadout = do
  t1 <- convMnistLayerS ker1 (constant x) bias1
  t2 <- convMnistLayerS ker2 t1 bias2
  let m1 = mapS reshapeS t2
      m2 = transpose2S m1
      denseLayer = weigthsDense <>$ m2 + asColumnS biasesDense
  denseRelu <- reluAct denseLayer
  returnLet $ weigthsReadout <>$ denseRelu + asColumnS biasesReadout

convMnistLenS
  :: forall kheight_minus_1 kwidth_minus_1 num_hidden out_channels
            in_height in_width.
     ( KnownNat kheight_minus_1, KnownNat kwidth_minus_1
     , KnownNat num_hidden, KnownNat out_channels
     , KnownNat in_height, KnownNat in_width )
  => Proxy kheight_minus_1
  -> Proxy kwidth_minus_1
  -> Proxy num_hidden
  -> Proxy out_channels
  -> Proxy in_height
  -> Proxy in_width
  -> (Int, [Int], [(Int, Int)], [OT.ShapeL])
convMnistLenS _ _ _ _ _ _ =
  ( 0
  , []
  , []
  , [ Data.Array.Shape.shapeT @'[ out_channels, 1
                                , kheight_minus_1 + 1, kwidth_minus_1 + 1 ]
    , Data.Array.Shape.shapeT @'[out_channels]
    , Data.Array.Shape.shapeT @'[ out_channels, out_channels
                                , kheight_minus_1 + 1, kwidth_minus_1 + 1 ]
    , Data.Array.Shape.shapeT @'[out_channels]
    , Data.Array.Shape.shapeT @'[ num_hidden
                                , out_channels
                                    GHC.TypeLits.*
                                      ((in_height + kheight_minus_1) `Div` 2
                                       + kheight_minus_1) `Div` 2
                                    GHC.TypeLits.*
                                      ((in_width + kwidth_minus_1) `Div` 2
                                       + kheight_minus_1) `Div` 2
                                ]
    , Data.Array.Shape.shapeT @'[num_hidden]
    , Data.Array.Shape.shapeT @'[SizeMnistLabel, num_hidden]
    , Data.Array.Shape.shapeT @'[SizeMnistLabel]
    ]
  )

convMnistS
  :: forall kheight_minus_1 kwidth_minus_1 num_hidden out_channels
            in_height in_width batch_size d r m.
     ( KnownNat kheight_minus_1, KnownNat kwidth_minus_1
     , KnownNat num_hidden, KnownNat out_channels
     , KnownNat in_height, KnownNat in_width, KnownNat batch_size
     , 1 <= kheight_minus_1
     , 1 <= kwidth_minus_1
     , DualMonad d r m )
  => OS.Array '[batch_size, 1, in_height, in_width] r
  -> DualNumberVariables d r
  -> m (DualNumber d (OS.Array '[SizeMnistLabel, batch_size] r))
convMnistS x variables = do
  let ker1 = varS variables 0
      bias1 = varS variables 1
      ker2 = varS variables 2
      bias2 = varS variables 3
      weigthsDense = varS variables 4
      biasesDense = varS variables 5
      weigthsReadout = varS variables 6
      biasesReadout = varS variables 7
  convMnistTwoS @kheight_minus_1 @kwidth_minus_1 @num_hidden @out_channels
                x ker1 bias1 ker2 bias2
                weigthsDense biasesDense weigthsReadout biasesReadout

convMnistLossFusedS
  :: forall kheight_minus_1 kwidth_minus_1 num_hidden out_channels
            in_height in_width batch_size d r m.
     ( KnownNat kheight_minus_1, KnownNat kwidth_minus_1
     , KnownNat num_hidden, KnownNat out_channels
     , KnownNat in_height, KnownNat in_width, KnownNat batch_size
     , 1 <= kheight_minus_1
     , 1 <= kwidth_minus_1
     , DualMonad d r m )
  => Proxy kheight_minus_1
  -> Proxy kwidth_minus_1
  -> Proxy num_hidden
  -> Proxy out_channels
  -> ( OS.Array '[batch_size, in_height, in_width] r
     , OS.Array '[batch_size, SizeMnistLabel] r )
  -> DualNumberVariables d r
  -> m (DualNumber d r)
convMnistLossFusedS _ _ _ _ (glyphS, labelS) variables = do
  let xs :: OS.Array '[batch_size, 1, in_height, in_width] r
      xs = OS.reshape glyphS
  result <- convMnistS @kheight_minus_1 @kwidth_minus_1
                       @num_hidden @out_channels
                       xs variables
  let targets2 = HM.tr $ HM.reshape (valueOf @SizeMnistLabel)
                       $ OS.toVector labelS
  vec <- lossSoftMaxCrossEntropyL targets2 (fromS2 result)
  returnLet $ scale (recip $ fromIntegral (valueOf @batch_size :: Int))
            $ sumElements0 vec

-- For simplicity, testing is performed in mini-batches of 1.
-- See RNN for testing done in batches.
convMnistTestS
  :: forall kheight_minus_1 kwidth_minus_1 num_hidden out_channels
            in_height in_width r.
     ( KnownNat kheight_minus_1, KnownNat kwidth_minus_1
     , KnownNat num_hidden, KnownNat out_channels
     , KnownNat in_height, KnownNat in_width
     , 1 <= kheight_minus_1
     , 1 <= kwidth_minus_1
     , IsScalar 'DModeGradient r )
  => Proxy r
  -> Proxy kheight_minus_1
  -> Proxy kwidth_minus_1
  -> Proxy num_hidden
  -> Proxy out_channels
  -> [( OS.Array '[in_height, in_width] r
      , OS.Array '[SizeMnistLabel] r )]
  -> Domains r
  -> r
convMnistTestS _ _ _ _ _ inputs parameters =
  let matchesLabels :: ( OS.Array '[in_height, in_width] r
                       , OS.Array '[SizeMnistLabel] r )
                    -> Bool
      matchesLabels (glyph, label) =
        let tx :: OS.Array '[1, 1, in_height, in_width] r
            tx = OS.reshape glyph
            nn :: DualNumberVariables 'DModeGradient r
               -> DualMonadValue r (DualNumber 'DModeGradient (OS.Array '[SizeMnistLabel, 1] r))
            nn = convMnistS @kheight_minus_1 @kwidth_minus_1
                            @num_hidden @out_channels
                            tx
            value = primalValue nn parameters
        in V.maxIndex (OS.toVector value) == V.maxIndex (OS.toVector label)
  in fromIntegral (length (filter matchesLabels inputs))
     / fromIntegral (length inputs)

max_stars_repo_name: Mikolaj/horde-ad

max_issues_repo_name: Mikolaj/horde-ad

max_forks_repo_name: Mikolaj/horde-ad


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE TemplateHaskell   #-}
module DBPnet.Correlation
    ( cisCorMat
    , transCorMat
    ) where

import           Bio.Data.Bed
import           Conduit
import Lens.Micro
import           Control.Monad
import qualified Data.ByteString.Char8       as B
import           Data.List.Split
import qualified Data.Matrix.Unboxed         as MU
import qualified Data.Matrix.Unboxed.Mutable as MUM
import qualified Data.Vector                 as V
import qualified Data.Vector.Unboxed         as U
import           Statistics.Correlation      (spearman)

import           DBPnet.Utils                (readLoops)

cisCorMat :: [(String, FilePath)]
          -> IO ([B.ByteString], MU.Matrix Double)
cisCorMat dat = do
    counts <- fmap V.fromList $ forM dat $ \(i, fl) -> do
        beds <- readBed fl
        return (i, U.fromList $ map (^.bdgValue) beds)

    let l = V.length counts
        header = map B.pack $ V.toList $ fst . V.unzip $ counts
        mat :: MU.Matrix Double
        mat = MU.create $ do
            m <- MUM.new (l,l)
            forM_ [0..l-1] $ \i -> do
                let (_, tfA) = counts V.! i
                forM_ [0..l-1] $ \j -> case () of
                    _ | j < i -> MUM.read m (j,i) >>= MUM.write m (i,j)
                      | otherwise -> do
                          let (_, tfB) = counts V.! j
                              vec = U.filter (\(a,b) -> a /= 0 || b /= 0) $ U.zip tfA tfB
                              cor | U.length vec < 500 = 0
                                  | otherwise = spearman vec
                          MUM.write m (i,j) cor
            return m
    return (header, mat)

transCorMat :: [(String, FilePath)] -> FilePath
            -> IO ([B.ByteString], MU.Matrix Double)
transCorMat dat hic = do
    loops <- concatMap (\(a,b) -> [a,b]) <$> readLoops hic

    counts <- fmap V.fromList $ forM dat $ \(i,fl) -> do
        beds <- readBed fl :: IO [BEDGraph]
        let rs = runIdentity $ runConduit $ yieldMany loops .|
                intersectBedWith fn beds .| sinkList
            fn _ [] = 0
            fn _ x = maximum $ map (^.bdgValue) x
        return (i, U.fromList $ map (\[a,b] -> (a,b)) $ chunksOf 2 rs)

    let l = V.length counts
        header = map B.pack $ V.toList $ fst . V.unzip $ counts
        mat :: MU.Matrix Double
        mat = MU.create $ do
            m <- MUM.new (l,l)
            forM_ [0..l-1] $ \i -> do
                let (_, tfA) = counts V.! i
                forM_ [0..l-1] $ \j -> case () of
                    _ | j < i -> MUM.read m (j,i) >>= MUM.write m (i,j)
                      | otherwise -> do
                          let (_, tfB) = counts V.! j
                              vec = U.filter (\(a,b) -> a /= 0 || b /= 0) $ U.zipWith f tfA tfB
                              cor | U.length vec < 500 = 0
                                  | otherwise = spearman vec
                          MUM.write m (i,j) cor
            return m
    return (header, mat)
  where
    f (x1,x2) (y1,y2) | x1 > x2 = (x1, y2)
                      | otherwise = (x2, y1)

max_stars_repo_name: kaizhang/DBPnet

max_issues_repo_name: kaizhang/DBPnet

max_forks_repo_name: kaizhang/DBPnet


text:
module AI.Learning.Bootstrap where

import Control.Monad.Random
import Foreign.Storable (Storable)
import Numeric.LinearAlgebra
import qualified Data.List as L

import AI.Util.Matrix
import AI.Util.Util

---------------
-- Bootstrap --
---------------

-- |Generate a bootstrap sample of size @sz@.
genBootstrapSample :: RandomGen g => Int -> Rand g [Int]
genBootstrapSample sz = go sz []
    where go 0 accum = return accum
          go n accum = do
            i <- getRandomR (0,sz-1)
            go (n - 1) (i:accum)

sampleVector :: (Storable a, RandomGen g) => Vector a -> Rand g (Vector a)
sampleVector v = do
    idx <- genBootstrapSample (dim v)
    return (v `subRefVec` idx)

sampleMatrixRows :: (Element a, RandomGen g) => Matrix a -> Rand g (Matrix a)
sampleMatrixRows m = do
    idx <- genBootstrapSample (rows m)
    return $ m `subRefRows` idx

sampleMatrixCols :: (Element a, RandomGen g) => Matrix a -> Rand g (Matrix a)
sampleMatrixCols m = do
    idx <- genBootstrapSample (cols m)
    return $ m `subRefCols` idx

-- |Generate a bootstrap sample of a statistic from a data set of type /a/.
bootStrapResample :: RandomGen g =>
                     (a -> Rand g a)    -- Sampling function
                  -> (a -> b)           -- Statistic to be resampled
                  -> Int                -- Number of resamples
                  -> a                  -- Data
                  -> Rand g [b]
bootStrapResample sample func nSamples x = go [] nSamples
    where
        go samples 0 = return samples
        go samples n = do
            x' <- sample x
            go (func x' : samples) (n-1)

max_stars_repo_name: cagix/aima-haskell


text:
module NN.NeuralNetworkTest
  (
    testsNeuralNetwork
  ) where

import Test.HUnit
import AI.HNN.FF.Network
import Numeric.LinearAlgebra.HMatrix hiding (corr)

import NUtil
import NN.NeuralNetwork
import qualified Data.Map.Strict as Map

nnVars :: Vars
nnVars = Map.fromList [("numberInputs", 2)
                      ,("numberHidden", 2)
                      ,("numberOutputs", 1)
                      ,("timesToTrain", 3)
                      ,("learningRate", 0.8)
                      ,("sigmoidOrTanh", 0)]

samples :: Samples Double
samples = [ fromList [0, 0] --> fromList [0]
          , fromList [0, 1] --> fromList [1]
          , fromList [1, 0] --> fromList [1]
          , fromList [1, 1] --> fromList [0]
          ]

expectedOutputVals :: [Double]
expectedOutputVals = [0,1,1,0]

testsNeuralNetwork :: Test
testsNeuralNetwork = TestList [ TestLabel "label" testsLabel]

testsLabel :: Test
testsLabel = TestList [ TestLabel "test1" testCreate
                      , TestLabel "test2" testGet
                      , TestLabel "test3" testTrainV]

testCreate :: Test
testCreate = TestCase (do
                          net <- create nnVars :: IO Net
                          let smartNet = trainNTimes 1500 0.8 tanh tanh' (_network net) samples
--                              outputValues = map (get smartNet . fst) samples
                              outputValues = convert1 $ map (output smartNet tanh . fst) samples
                          assertEqual "testCreate"
                            True
                            (floatVectorCompare outputValues expectedOutputVals 0.2))

testGet :: Test
testGet = TestCase (do
                       net <- create nnVars :: IO Net
                       let smartNet = trainNTimes 1000 0.8 tanh tanh' (_network net) samples
                           outputValues = convert1 $ map (get (Net smartNet nnVars) . fst) samples
                       assertEqual "testGet"
                         True
                         (floatVectorCompare outputValues expectedOutputVals 0.2))

testTrainV :: Test
testTrainV = TestCase (do
                         net <- create nnVars :: IO Net
                         let smartNet = repeatNet trainV net samples 1000
                             outputValues = convert1 $ map (get smartNet . fst) samples
                         assertEqual "testTrainV"
                           True
                           (floatVectorCompare outputValues expectedOutputVals 0.2))

convert1 :: [Vector Double] -> [Double]
convert1 = map (! 0)

floatVectorCompare :: [Double] -> [Double] -> Double -> Bool
floatVectorCompare a b tol = foldr (&&) True $ zipWith (\x y -> abs (x - y) < tol ) a b

repeatNet :: (n -> c -> n) -> n -> c -> Int -> n
repeatNet _ net _ 0 = net
repeatNet f net c n = let net' = f net c in repeatNet f net' c (n-1)

max_stars_repo_name: eklinkhammer/neural-algorithms

max_issues_repo_name: eklinkhammer/neural-algorithms

max_forks_repo_name: eklinkhammer/neural-algorithms


text:
module Main where

import Numeric.NLOPT
import Numeric.LinearAlgebra (fromList, toList)

main :: IO ()
main = do
    let objf x = let [x1, x2, x3] = toList x
                 in 3 * x1 * x1 + 2 * x1 * x2 + x1 * x3
                    + 2.5 * x2 * x2 + 2 * x2 * x3 + 2 * x3 * x3
                    - 8 * x1 - 3 * x2 - 3 * x3

        objfGrad x = let [x1, x2, x3] = toList x
                     in fromList [ 6 * x1 + 2 * x2 + x3 - 8
                                 , 2 * x1 + 5 * x2 + 2 * x3 - 3
                                 , x1 + 2 * x2 + 4 * x3 - 3 ]

        objfD x = (objf x, objfGrad x)

        stop = ObjectiveRelativeTolerance 1e-9 :| []

        c1f x = let [x1, _, x3] = toList x in x1 + x3 - 3
        c1fGrad _ = fromList [1, 0, 1]
        c1fD x = (c1f x, c1fGrad x)
        c1 = EqualityConstraint (Scalar c1fD) 1e-6

        c2f x = let [_, x2, x3] = toList x in x2 + x3
        c2fGrad _ = fromList [0, 1, 1]
        c2fD x = (c2f x, c2fGrad x)
        c2 = EqualityConstraint (Scalar c2fD) 1e-6

        algorithm = SLSQP objfD [] [] [c1, c2]
        problem = LocalProblem 3 stop algorithm

        x0 = fromList [0, 0, 0]

        solC = minimizeLocal problem x0

    print solC

    putStrLn "-- Done."

max_stars_repo_name: cbpark/hmatrix-nlopt-test

max_issues_repo_name: cbpark/hmatrix-nlopt-test

max_forks_repo_name: cbpark/hmatrix-nlopt-test


text:
{-# LANGUAGE BangPatterns        #-}
{-# LANGUAGE CPP                 #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TupleSections       #-}
{-# LANGUAGE TypeFamilies        #-}
{-# LANGUAGE TypeOperators       #-}

import           Control.Monad.Random
import           Control.Monad.Trans.Except

import qualified Data.ByteString              as B
import           Data.Char                    (isUpper, toLower, toUpper)
import           Data.List                    (foldl')
import qualified Data.Map                     as M
import           Data.Maybe                   (fromMaybe)
import           Data.Serialize
import           Data.Singletons.Prelude
import           Data.Vector                  (Vector)
import qualified Data.Vector                  as V
import           GHC.TypeLits
import           Numeric.LinearAlgebra.Static (konst)
import           Options.Applicative

import           Grenade
import           Grenade.Recurrent
import           Grenade.Utils.OneHot

import           System.IO.Unsafe             (unsafeInterleaveIO)

-- The defininition for our natural language recurrent network.
-- This network is able to learn and generate simple words in
-- about an hour.
--
-- Grab the input from
-- https://raw.githubusercontent.com/karpathy/char-rnn/master/data/tinyshakespeare/input.txt
--
-- This is a first class recurrent net.
--
-- The F and R types are tagging types to ensure that the runner and
-- creation function know how to treat the layers.
--
-- As an example, here's a short sequence generated.
--
-- > KING RICHARD III:
-- > And as the heaven her his words, we the son, I show sand stape but the lament to shall were the sons with a strend

type FF = FeedForward
type R = Recurrent

-- The definition of our network
type Shakespeare = RecurrentNetwork '[ R (LSTM 40 80), R (LSTM 80 40), FF (FullyConnected 40 40), FF Logit]
                                    '[ 'D1 40, 'D1 80, 'D1 40, 'D1 40, 'D1 40 ]

-- The definition of the "sideways" input, which the network is fed recurrently.
type Shakespearian = RecurrentInputs  '[ R (LSTM 40 80), R (LSTM 80 40), FF (FullyConnected 40 40), FF Logit]

randomNet :: IO Shakespeare
randomNet = randomRecurrent

-- | Load the data files and prepare a map of characters to a compressed int representation.
loadShakespeare :: FilePath -> ExceptT String IO (Vector Int, M.Map Char Int, Vector Char)
loadShakespeare path = do
  contents     <- lift $ readFile path
  let annotated = annotateCapitals contents
  (m,cs)       <- ExceptT . return $ hotMap (Proxy :: Proxy 40) annotated
  hot          <- ExceptT . return . note "Couldn't generate hot values" $ traverse (`M.lookup` m) annotated
  return (V.fromList hot, m, cs)

trainSlice :: Optimizer o -> Shakespeare -> Shakespearian -> Vector Int -> Int -> Int -> (Shakespeare, Shakespearian)
trainSlice !opt !net !recIns input offset size =
  let e = fmap (x . oneHot) . V.toList $ V.slice offset size input
  in case reverse e of
    (o : l : xs) ->
      let examples = reverse $ (l, Just o) : ((,Nothing) <$> xs)
      in  trainRecurrent opt net recIns examples
    _ -> error "Not enough input"
    where
      x = fromMaybe (error "Hot variable didn't fit.")

runShakespeare :: ShakespeareOpts -> ExceptT String IO ()
runShakespeare opts = do
  (shakespeare, oneHotMap, oneHotDictionary) <- loadShakespeare $ trainingFile opts
  (net0, i0) <- lift $
    case loadPath opts of
      Just loadFile -> netLoad loadFile
      Nothing       -> (,0) <$> randomNet

  (trained, bestInput) <- lift $ foldM (\(!net, !io) size -> do
    xs <- take (iterations opts `div` 10) <$> getRandomRs (0, length shakespeare - size - 1)
    let (!trained, !bestInput) = foldl' (\(!n, !i) offset -> trainSlice (rate opts) n i shakespeare offset size) (net, io) xs
    results <- take 1000 <$> generateParagraph trained bestInput (temperature opts) oneHotMap oneHotDictionary ( S1D $ konst 0)
    putStrLn ("TRAINING STEP WITH SIZE: " ++ show size)
    putStrLn (unAnnotateCapitals results)
    return (trained, bestInput)
    ) (net0, i0) $ replicate 10 (sequenceSize opts)

  case savePath opts of
    Just saveFile -> lift . B.writeFile saveFile $ runPut (put trained >> put bestInput)
    Nothing       -> return ()

generateParagraph :: forall layers shapes n a. (Last shapes ~ 'D1 n, Head shapes ~ 'D1 n, KnownNat n, Ord a)
  => RecurrentNetwork layers shapes
  -> RecurrentInputs layers
  -> RealNum
  -> M.Map a Int
  -> Vector a
  -> S ('D1 n)
  -> IO [a]
generateParagraph n s temp hotmap hotdict =
  go s
    where
  go x y =
    do let (_, ns, o) = runRecurrent n x y
       un            <- sample temp hotdict o
       Just re       <- return $ makeHot hotmap un
       rest          <- unsafeInterleaveIO $ go ns re
       return (un : rest)

data ShakespeareOpts = ShakespeareOpts {
    trainingFile :: FilePath
  , iterations   :: Int
  , rate         :: Optimizer 'SGD
  , sequenceSize :: Int
  , temperature  :: RealNum
  , loadPath     :: Maybe FilePath
  , savePath     :: Maybe FilePath
  }

shakespeare' :: Parser ShakespeareOpts
shakespeare' = ShakespeareOpts <$> argument str (metavar "TRAIN")
                               <*> option auto (long "examples" <> short 'e' <> value 1000000)
                               <*> (OptSGD
                                    <$> option auto (long "train_rate" <> short 'r' <> value 0.01)
                                    <*> option auto (long "momentum" <> value 0.95)
                                    <*> option auto (long "l2" <> value 0.000001)
                                    )
                               <*> option auto (long "sequence-length" <> short 's' <> value 50)
                               <*> option auto (long "temperature" <> short 't' <> value 0.4)
                               <*> optional (strOption (long "load"))
                               <*> optional (strOption (long "save"))

main :: IO ()
main = do
    shopts <- execParser (info (shakespeare' <**> helper) idm)
    res <- runExceptT $ runShakespeare shopts
    case res of
      Right () -> pure ()
      Left err -> putStrLn err

netLoad :: FilePath -> IO (Shakespeare, Shakespearian)
netLoad modelPath = do
  modelData <- B.readFile modelPath
  either fail return $ runGet get modelData

-- Replace capitals with an annotation and the lower case letter
-- http://fastml.com/one-weird-trick-for-training-char-rnns/
annotateCapitals :: String -> String
annotateCapitals (x : rest)
    | isUpper x
    = '^' : toLower x : annotateCapitals rest
    | otherwise
    = x : annotateCapitals rest
annotateCapitals []
    = []

unAnnotateCapitals :: String -> String
unAnnotateCapitals ('^' : x : rest)
    = toUpper x : unAnnotateCapitals rest
unAnnotateCapitals (x : rest)
    =  x : unAnnotateCapitals rest
unAnnotateCapitals []
    = []

-- | Tag the 'Nothing' value of a 'Maybe'
note :: a -> Maybe b -> Either a b
note a = maybe (Left a) Right

max_stars_repo_name: schnecki/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
{-# LANGUAGE BangPatterns          #-}
{-# LANGUAGE DeriveDataTypeable    #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE TemplateHaskell       #-}
{-# LANGUAGE TypeFamilies          #-}
-- |
-- Monoids for calculating various statistics in constant space. This
-- module contains algorithms that should be generally avoided unless
-- there's specific reason to use them.
module Statistics.Monoid.Extra (
    -- * Mean
    WelfordMean(..)
  , asWelfordMean
  , MeanKahan(..)
  , asMeanKahan
    -- $references
  ) where

import Data.Semigroup               (Semigroup(..))
import Data.Monoid                  (Monoid(..))
import Data.Data                    (Typeable,Data)
import Data.Vector.Unboxed          (Unbox)
import Data.Vector.Unboxed.Deriving (derivingUnbox)
import Numeric.Sum
import GHC.Generics                 (Generic)

import Statistics.Monoid.Class
-- COMPAT
import qualified Data.Vector.Generic         -- Needed for GHC7.4
import qualified Data.Vector.Generic.Mutable -- Needed for GHC7.4

----------------------------------------------------------------
-- Mean
----------------------------------------------------------------

-- | Incremental calculation of mean. Sum of elements is calculated
--   using compensated Kahan summation. It's provided only for sake of
--   completeness. 'Statistics.Monoid.Numeric.KBNSum' should be used
--   instead.
data MeanKahan = MeanKahan !Int !KahanSum
             deriving (Show,Eq,Typeable,Data,Generic)

asMeanKahan :: MeanKahan -> MeanKahan
asMeanKahan = id

instance Semigroup MeanKahan where
  (<>) = mappend
  {-# INLINE (<>) #-}

instance Monoid MeanKahan where
  mempty = MeanKahan 0 mempty
  MeanKahan 0  _  `mappend` m               = m
  m               `mappend` MeanKahan 0  _  = m
  MeanKahan n1 s1 `mappend` MeanKahan n2 s2 = MeanKahan (n1+n2) (s1 `mappend` s2)

instance Real a => StatMonoid MeanKahan a where
  addValue (MeanKahan n m) x = MeanKahan (n+1) (addValue m x)

instance CalcCount MeanKahan where
  calcCount (MeanKahan n _) = n
instance CalcMean MeanKahan where
  calcMean (MeanKahan 0 _) = Nothing
  calcMean (MeanKahan n s) = Just (kahan s / fromIntegral n)

-- | Incremental calculation of mean. Note that this algorithm doesn't
--   offer better numeric precision than plain summation. Its only
--   advantage is protection against double overflow:
--
--   > λ> calcMean $ asMeanKBN     $ reduceSample (replicate 100 1e308)
--   > Just NaN
--   > λ> calcMean $ asWelfordMean $ reduceSample (replicate 100 1e308)
--   > Just 1.0e308
--
--   Unless this feature is needed 'Statistics.Monoid.Numeric.KBNSum'
--   should be used. Algorithm is due to Welford [Welford1962]
data WelfordMean = WelfordMean !Int    -- Number of entries
                               !Double -- Current mean
  deriving (Show,Eq,Typeable,Data,Generic)

-- | Type restricted 'id'
asWelfordMean :: WelfordMean -> WelfordMean
asWelfordMean = id

instance Semigroup WelfordMean where
  (<>) = mappend
  {-# INLINE (<>) #-}

instance Monoid WelfordMean where
  mempty = WelfordMean 0 0
  mappend (WelfordMean 0 _) m = m
  mappend m (WelfordMean 0 _) = m
  mappend (WelfordMean n x) (WelfordMean k y)
    = WelfordMean (n + k) ((x*n' + y*k') / (n' + k'))
    where
      n' = fromIntegral n
      k' = fromIntegral k
  {-# INLINE mempty  #-}
  {-# INLINE mappend #-}

-- | \[ s_n = s_{n-1} + \frac{x_n - s_{n-1}}{n} \]
instance Real a => StatMonoid WelfordMean a where
  addValue (WelfordMean n m) !x
    = WelfordMean n' (m + (realToFrac x - m) / fromIntegral n')
    where
      n' = n+1
  {-# INLINE addValue #-}

instance CalcCount WelfordMean where
  calcCount (WelfordMean n _) = n
instance CalcMean WelfordMean where
  calcMean (WelfordMean 0 _) = Nothing
  calcMean (WelfordMean _ m) = Just m

----------------------------------------------------------------
-- Unboxed instances
----------------------------------------------------------------

derivingUnbox "MeanKahan"
  [t| MeanKahan -> (Int,Double,Double) |]
  [| \(MeanKahan a (KahanSum b c)) -> (a,b,c)   |]
  [| \(a,b,c) -> MeanKahan a (KahanSum b c) |]

derivingUnbox "WelfordMean"
  [t| WelfordMean -> (Int,Double) |]
  [| \(WelfordMean a b) -> (a,b)  |]
  [| \(a,b) -> WelfordMean a b    |]

-- $references
--
-- * [Welford1962] Welford, B.P. (1962) Note on a method for
--   calculating corrected sums of squares and
--   products. /Technometrics/
--   4(3):419-420. <http://www.jstor.org/stable/1266577>


text:
{-# LANGUAGE ConstraintKinds #-}
module Spectrogram(
     spectrogram
    ) where

import Signal
import Windows
import qualified Data.Vector.Unboxed as U
import Data.Vector.Unboxed((!))
import Data.Complex
import Data.Bits
import Common
import Transform
import Graphics.PDF
import Plot
import Control.Applicative((<$>))
import Data.List(foldl1')

import Debug.Trace

debug a = trace (show a) a

_spectrum :: (Sample a, Num a)
          => Int
          -> Double
          -> U.Vector a
          -> U.Vector Double
_spectrum n t d =
    let l = 1 `shiftL` n :: Int
        complexd = U.map (:+ 0) . U.map toDouble $ d
        m (x :+ y) =
            let x' = toDouble x
                y' = toDouble y
            in
            t*(x'*x' + y'*y') / fromIntegral l
    in
    U.map m . fft  $ complexd

prologSpect l b wi hi Nothing _ = return ()
prologSpect l b wi hi (Just r) _ = do
    withNewContext $ do
        applyMatrix $ translate (l :+ b)
        drawXObject r

genPicture :: ((Double, Double) -> Double)
           -> Int
           -> Int
           -> CoordinateMapping Double Double
           -> PDF (Maybe (PDFReference RawImage))
genPicture value w h (plotToPixel, pixelToPlot) = Just <$> do
    runPixmap w h $ do
        let drawValue p = do
            let pl = pixelToPlot p
                v = log (1 + value pl) / log 2.0
                col = Rgb v 0 0
            setColor col 1.0
            pixel p

        mapM_ drawValue ( [ fromIntegral col :+ fromIntegral row | col <- [0..w-1], row <- [0..h-1]])

-- We display only up to sampling frequency / 2 because the spectrum is even.
-- So, it allow a zooming a the frequency that matters instead of displaying redudant information
spectrogram :: (Sample a, Show a)
            => Sampled Time a -- ^ Signal
            -> Time -- ^ Duration
            -> (Int -> Int -> a -> a) -- ^ Window
            -> Int -- ^ Overlap
            -> StyledSignal Double Double
spectrogram signal duration window overlap =
    let samplingF = rate signal
        winSize :: Int
        winExp = 8
        winSize = 1 `shiftL` winExp
        freqResolution = (getF samplingF) / fromIntegral winSize
        frames = frameWithWinAndOverlap winSize overlap window signal
        nbFrames :: Int
        nbFrames = (floor (getT duration * getF (rate frames)))
        listOfSpectra = map (U.slice 0 (winSize `shiftR` 1)) .
                        map (_spectrum winExp (1.0 / getF samplingF)) . takeS nbFrames $ (getSignal frames)
        theSpectrum = U.concat listOfSpectra

        thePeak = U.maximum theSpectrum
        energyLocation = foldl1' (U.zipWith max)  listOfSpectra
        (lowest,_) = U.span (<= thePeak * 0.1) energyLocation
        (highest,_) = U.span (<= thePeak * 0.1) (U.reverse energyLocation)
        startFrequency = fromIntegral (U.length lowest) / fromIntegral winSize * getF samplingF
        stopFrequency = fromIntegral (U.length energyLocation - U.length highest) / fromIntegral winSize * getF samplingF
        normalizedSpectrum = U.map (/ thePeak) theSpectrum
        totalElements = U.length normalizedSpectrum
        value (i,j) =
            let pos = floor i*(winSize `shiftR` 1) + floor (j / freqResolution)
            in
            if pos >= totalElements
                then
                    normalizedSpectrum ! (totalElements - 1)
                else
                    normalizedSpectrum ! pos
        left = leftMargin (defaultPlotStyle)
        bottom = bottomMargin (defaultPlotStyle)
        style = defaultPlotStyle { title = Just "Spectrogram"
                                 , horizontalBounds = Just (0, fromIntegral nbFrames)
                                 , verticalBounds = Just (startFrequency,stopFrequency)
                                 , prolog = prologSpect left bottom
                                 , prologRsrc = genPicture value
                                 , horizontalLabel = Just "frame"
                                 , verticalLabel = Just "Hz"
                                 }
        l=[] :: [[Double]]
        theFrames = [] :: [Double]
    in
    signalsWithStyle False theFrames l style

max_stars_repo_name: alpheccar/Signal

max_issues_repo_name: cpehle/Signal

max_forks_repo_name: cpehle/Signal


text:
-----------------------------------------------------------------------------
-- |
-- Module      :  FlashADC
-- Copyright   :  José Edil Guimarães de Medeiros
-- License     :  BSD-style (see the file LICENSE)
--
-- Maintainer  :  j.edil@ene.unb.br
-- Stability   :  experimental
-- Portability :  portable
--
-- This module shows the model of a flash architecture ADC.
--
-- =Running the demo
--
-- To run the demo you need the ForSyDe installed in your
-- environment. For plotting, you need the
-- <http://gnuplot.sourceforge.net Gnuplot package>.
--
-- To calculate 10000 samples for R=1e3 Omhs, C=1e-6 F, T=1e-6run in @ghci@:
--
-- >>> simulate 4
--
-- To plot the response for a sine wave input, run in @ghci@:
--
-- >>> plotOutput 4 0.01 1000
-----------------------------------------------------------------------------

module ForSyDe.Shallow.Example.Synchronous.FlashADC where

import ForSyDe.Shallow
import Data.Complex

type Resistance = Double
type Voltage = Double
type Code = Integer
type Bits = Integer
type Time = Double
type Frequency = Double

-- | 'flashADC' is the top level module.
flashADC :: [Resistance]                -- ^ Resistance values
          -> Signal Voltage             -- ^ Input signal
          -> Signal Code                -- ^ Output signal
flashADC resistors input = decoder $ compNetwork input $ resNetwork resistors

-- | 'decoder' takes the thermometer code and outputs integers.
decoder :: [Signal Bits]                -- ^ Bit inputs
        -> Signal Code                  -- ^ Output signal
decoder = foldl1 (zipWithSY (+))

-- | 'compNetwork' implements the comparator array.
compNetwork :: Signal Voltage           -- ^ ADC input signal
            -> [Voltage]                -- ^ Voltage thresholds
            -> [Signal Bits]            -- ^ Bit outputs
compNetwork input = zipWith (\i v -> mapSY (comparator v) i) (repeat input)

-- | 'comparator' is the one bit quantizer.
comparator :: Voltage                   -- ^ (+) input
           -> Voltage                   -- ^ (-) input
           -> Bits                      -- ^ Output
comparator i v
  | v <= i = 0
  | otherwise = 1

-- | 'resNetwork' implements the resistor voltage scaling network.
resNetwork :: [Resistance]              -- ^ Resistor values
           -> [Voltage]                 -- ^ Threshold voltages
resNetwork resistors = init $ tail $ scanl (\v r -> v + vdd * r / (sumR)) 0 resistors
  where vdd = 1
        sumR = sum resistors

-- | 'simulate' takes the system parameters and runs the simulation
-- with a sine wave input.
simulate :: Int                         -- ^ ADC resolution
         -> Int                         -- ^ Number of samples
         -> Signal Code                 -- ^ ADC output
simulate res n = flashADC resistors input
  where resistors = replicate (2^res) 1
        input = sineWave' 0.5 0.5 50 n

-- | 'plotFFT' runs the simulation with a sine wave input
plotFFT :: Int                  -- ^ ADC resolution
        -> Int                  -- ^ FFT depth (power of 2)
        -> IO String            -- ^ plot
plotFFT res n = plotCT' (toRational 1.0) [(g_out, "g")]
  where g_out = d2aConverter DAhold (toRational 1.0) $ signal g
        a = fromSignal $ simulate res n
        b = (:+) <$> (map ((\x -> x/2^res - 0.5).fromIntegral) a)
        c = zipWith ($) b $ repeat 0.0
        d = vector c
        e = fromVector $ fft n d
        f = map ((*20).(logBase 10).(\x -> 2*x/(fromIntegral n)) . magnitude) e
        g = (take (n `div` 2) f)

-- | 'plotOutput' uses the CTLib plot capabilities to plot the
-- output. In a later version, a plotter to Synchornous signals will be
-- developed.
plotOutput :: Int                       -- ^ Resolution
           -> Double                    -- ^ Discretization timestep
           -> Int                       -- ^ Number of samples
           -> IO String                 -- ^ plot
plotOutput res t n = plotCT' (toRational t) [(output, "output")]
  where output = d2aConverter DAhold (toRational t) adcSignal
        adcSignal = signal $ map (toRational) $ fromSignal $
                    simulate res n

-- | 'sineWave' is an auxiliary function that creates a sine wave signal for
-- simulation purposes.
sineWave' :: Voltage                     -- ^ Amplitude
          -> Voltage                     -- ^ Offset
          -> Frequency                   -- ^ Frequency
          -> Int                         -- ^ Number of samples
          -> Signal Time                 -- ^ Output signal
sineWave' amp offset freq n = signal sineW
  where sineW = map (\t -> (amp/2) * sin(2*pi*freq*t) + offset) grid
        grid = linspace 0 1 n

-- | 'linspace' is an auxiliary function that creates an uniform
-- spaced sampling grid for simulation purposes.
linspace :: Time                        -- ^ Start time
         -> Time                        -- ^ Stop time
         -> Int                         -- ^ Number of samples
         -> [Time]                      -- ^ Sampling grid
linspace start stop n = map (scale) nodes
  where n' = fromIntegral n
        nodes = map (fromIntegral) [0..(n-1)]
        scale t = t * length / (n' - 1) + start
        length = stop - start

max_stars_repo_name: Rojods/forsyde-shallow-examples

max_issues_repo_name: Rojods/forsyde-shallow-examples

max_forks_repo_name: Rojods/forsyde-shallow-examples