text:
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE TypeFamilies #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}
{-# OPTIONS_GHC -fno-warn-unused-do-bind #-}

module Network.NTP where

import           Control.Applicative ((<$>))
import           Control.Concurrent (forkIO)
import           Control.Concurrent.STM.TChan
import           Control.Exception (IOException, bracket, handle)
import           Control.Monad.IO.Class
import           Control.Monad.Loops (unfoldM)
import           Control.Monad.STM (atomically)
import           Control.Monad.State
import           Data.Bits (shiftR, shiftL, (.|.), (.&.))
import qualified Data.ByteString as B
import           Data.IORef
import           Data.Int (Int64)
import           Data.Maybe (listToMaybe)
import           Data.Serialize (encode, decode)
import qualified Data.Time.Clock as T
import qualified Data.Vector.Generic as Generic
import qualified Data.Vector.Generic.Mutable as MGeneric
import qualified Data.Vector.Unboxed as U
import           Data.Word (Word32, Word64)
import           Network.Socket hiding (send, sendTo, recv, recvFrom)
import           Network.Socket.ByteString
import qualified Statistics.Function as S
import qualified Statistics.Sample as S

import           Network.NTP.Types hiding (Server)
import           System.Counter
import           Text.PrefixSI

--import           Debug.Trace
--import           Text.Printf (printf)

------------------------------------------------------------------------
-- NTP Monad

newtype NTP a = NTP { runNTP :: StateT NTPData IO a }
    deriving (Functor, Monad, MonadIO, MonadState NTPData)

data NTPData = NTPData {
      ntpSocket   :: Socket
    , ntpIncoming :: TChan AddrPacket
    , ntpRunning  :: IORef Bool
    }

data AddrPacket = AddrPacket ClockIndex SockAddr Packet
type Logger = String -> IO ()

withNTP :: Logger -> NTP a -> IO a
withNTP writeLog =
    fmap fst . withSocketsDo . bracket create destroy . runStateT . runNTP
  where
    create = do
        ntpSocket   <- initSocket
        ntpIncoming <- newTChanIO
        ntpRunning  <- newIORef True

        let enqueue = atomically . writeTChan ntpIncoming
            receiveLoop = do
              pkt <- receive ntpSocket
              ok  <- readIORef ntpRunning
              if ok then either writeLog enqueue pkt >> receiveLoop
                    else return ()

        -- start receive packet loop
        forkIO receiveLoop

        return NTPData{..}

    destroy NTPData{..} = do
        writeIORef ntpRunning False
        sClose ntpSocket

    initSocket = do
        sock <- socket AF_INET Datagram defaultProtocol
        bindSocket sock (SockAddrInet 0 0)
        return sock

------------------------------------------------------------------------
-- Clock

type ClockIndex = Int64
type ClockDiff  = Int64
type ClockDiffD = Double
type Seconds    = Double

data Clock = Clock {
      clockTime0      :: Time
    , clockIndex0     :: ClockIndex
    , clockFrequency  :: Double
    , clockPrecision  :: Word64
    }

-- | Initializes a 'Clock' based on the high precision counter.
initCounterClock :: Logger -> IO Clock
initCounterClock logger = do
    logger "Analyzing high precision counter..."

    counter <- analyzeCounter

    let clockFrequency = fromIntegral (cntFrequency counter)
        clockPrecision = cntPrecision counter
        freq = clockFrequency
        prec = fromIntegral clockPrecision / freq

    logger ("Frequency = " ++ showFreq freq)
    logger ("Precision = " ++ showTime prec)

    clockTime0  <- fromUTCTime <$> T.getCurrentTime
    clockIndex0 <- getCurrentIndex

    return Clock{..}

-- | Gets the current time according to the high res counter.
getCurrentIndex :: IO ClockIndex
getCurrentIndex = fromIntegral <$> readCounter

-- | Gets the current time according to the 'Clock'.
getCurrentTime :: Clock -> IO Time
getCurrentTime clock = clockTime clock <$> getCurrentIndex

-- | Gets the time at the given index.
clockTime :: Clock -> ClockIndex -> Time
clockTime Clock{..} index = clockTime0 `add` fromSeconds off
  where
    off  = diff / freq
    -- TODO: Unsigned subtract breaks when index is earlier than clockIndex0.
    --       This is the case for servers which have not responded for a while
    --       as we continually move clockIndex0 forward using adjustOrigin.
    diff = fromIntegral (index - clockIndex0)
    freq = clockFrequency

-- | Gets the index at the given time.
clockIndex :: Clock -> Time -> ClockIndex
clockIndex Clock{..} time = clockIndex0 + round (diff * freq)
  where
    diff = toSeconds (time `sub` clockTime0)
    freq = clockFrequency

-- | Gets a clock difference in seconds
fromDiff :: Clock -> ClockDiff -> Seconds
fromDiff Clock{..} d = fromIntegral d / clockFrequency

adjustFrequency :: Double -- ^ the drift rate in seconds/second
                -> Clock -> Clock
adjustFrequency adj c = c
    { clockFrequency = (1 - adj) * clockFrequency c }

adjustOffset :: Seconds -- ^ the offset in seconds
             -> Clock -> Clock
adjustOffset adj c = c
    { clockTime0 = clockTime0 c `add` fromSeconds adj }

adjustOrigin :: ClockIndex -- ^ the new clockIndex0
             -> Clock -> Clock
adjustOrigin newIndex0 c = c
    { clockTime0  = clockTime c newIndex0
    , clockIndex0 = newIndex0 }

setFrequency :: Double -> Clock -> Clock
setFrequency freq c = c { clockFrequency = freq }

------------------------------------------------------------------------
-- Server

data Server = Server {
      svrHostName     :: HostName
    , svrSockAddr     :: SockAddr
    , svrRawSamples   :: U.Vector Sample
    , svrMinRoundtrip :: ClockDiff
    , svrBaseError    :: ClockDiff
    , svrStratum      :: Int
    , svrRefId        :: B.ByteString
    , svrRefAddr      :: HostAddress
    , svrClock        :: Clock
    }

svrName :: Server -> String
svrName svr | host /= addr = host ++ " (" ++ addr ++ ")"
            | otherwise    = host
  where
    host = svrHostName svr
    addr = (takeWhile (/= ':') . show . svrSockAddr) svr

svrAddr :: Server -> HostAddress
svrAddr Server{..} = case svrSockAddr of
    (SockAddrInet _ addr) -> addr
    _                     -> error "svrAddr: unexpected non-IPv4 address"

-- | Resolves a list of IP addresses registered for the specified
-- hostname and creates 'Server' instances for each of them.
resolveServer :: Clock -> HostName -> IO (Maybe Server)
resolveServer clock host = ioErrorToNothing $
    listToMaybe . map (mkServer . addrAddress) <$> getHostAddrInfo
  where
    mkServer addr = Server host addr U.empty 0 0 0 B.empty 0 clock

    getHostAddrInfo = getAddrInfo hints (Just host) (Just "ntp")

    hints = Just defaultHints { addrFamily     = AF_INET
                              , addrSocketType = Datagram }

    ioErrorToNothing = handle (\(_ :: IOException) -> return Nothing)

transmit :: Server -> NTP ()
transmit Server{..} = do
    NTPData{..} <- get
    liftIO $ do
        now <- fromIntegral <$> getCurrentIndex
        let msg = (encode . requestMsg . Time) now
        sendAllTo ntpSocket msg svrSockAddr

receive :: Socket -> IO (Either String AddrPacket)
receive sock = handleIOErrors $ do
    (bs, addr) <- recvFrom sock 128
    t          <- getCurrentIndex
    return $ case decode bs of
        Left err  -> Left err
        Right pkt -> Right (AddrPacket t addr pkt)
  where
    handleIOErrors = handle (\(e :: IOException) -> return (Left (show e)))

updateServers :: [Server] -> NTP [(Server, Bool)]
updateServers svrs = do
    NTPData{..} <- get
    liftIO $ do
        pkts <- unfoldM (tryReadChan ntpIncoming)
        return $ map (update pkts) (zip svrs (repeat False))
  where
    update :: [AddrPacket] -> (Server, Bool) -> (Server, Bool)
    update = fconcat . map go
      where
        go (AddrPacket t addr pkt) (svr, up)
          | addr == svrSockAddr svr = (updateServer svr t pkt, True)
          | otherwise               = (svr, up)

    fconcat :: [a -> a] -> a -> a
    fconcat = foldl (.) id

    tryReadChan chan = atomically $ do
        empty <- isEmptyTChan chan
        if empty then return Nothing
                 else Just <$> readTChan chan

updateServer :: Server -> ClockIndex -> Packet -> Server
updateServer svr t4 p@Packet{..} =
    adjustClock svr
    { svrRawSamples   = samples
    , svrMinRoundtrip = round minRoundtrip
    , svrBaseError    = round (3 * stdDevRoundtrip)
    , svrStratum      = fromIntegral ntpStratum
    , svrRefId        = decodeRefId p
    , svrRefAddr      = byteSwap32 ntpRefId
    }
  where
    newSample = (fromIntegral (unTime ntpT1), ntpT2, ntpT3, t4)

    -- force the evaluation of the last sample in the list
    -- to avoid building up lots of unevaluated thunks
    samples = U.take maxSamples (newSample `U.cons` svrRawSamples svr)

    minRoundtrip    = U.head bestRoundtrips
    stdDevRoundtrip = S.stdDev bestRoundtrips

    -- Trim roundtrips so that outliers due to network congestion are
    -- not included in our quality weights.
    roundtrips     = U.map (fromIntegral . roundtrip) samples
    bestRoundtrips = (U.take half . S.partialSort half) roundtrips
    half           = (U.length roundtrips + 1) `div` 2

byteSwap32 :: Word32 -> Word32
byteSwap32 n = a .|. b .|. c .|. d
  where
    a = n `shiftR` 24
    b = (n `shiftL` 8) .&. 0x00ff0000
    c = (n `shiftR` 8) .&. 0x0000ff00
    d = n `shiftL` 24

maxSamples :: Int
maxSamples = max phaseSamples freqSamples

freqSamples :: Int
freqSamples = round (1000 * samplesPerSecond)

phaseSamples :: Int
phaseSamples = round (50 * samplesPerSecond)

samplesPerSecond :: Double
samplesPerSecond = 0.5

adjustClock :: Server -> Server
adjustClock svr@Server{..} = svr { svrClock = adjust svrClock }
  where
    adjust = if isNaN phase then id else adjustOffset phase
           . if isNaN freq  then id else adjustFrequency freq
           . adjustOrigin earliestTime

    earliestTime = time1 (U.last svrRawSamples)

    phase  = S.meanWeighted (U.take phaseSamples weightedOffsets)

    (_, freq) = linearRegression (U.take freqSamples times)
                                 (U.take freqSamples weightedOffsets)

    weightedOffsets = U.zip offsets weights
    times   = U.map (fromDiff svrClock . (\x -> time4 x - earliestTime)) svrRawSamples
    offsets = U.map (toSeconds . offset svrClock) svrRawSamples
    weights = U.map (quality svrClock svr) svrRawSamples

------------------------------------------------------------------------
-- Sample Type

type Sample = (ClockIndex, Time, Time, ClockIndex)

time1 :: Sample -> ClockIndex
time1 (x,_,_,_) = x

time2 :: Sample -> Time
time2 (_,x,_,_) = x

time3 :: Sample -> Time
time3 (_,_,x,_) = x

time4 :: Sample -> ClockIndex
time4 (_,_,_,x) = x

deriving instance MGeneric.MVector U.MVector Time
deriving instance Generic.Vector U.Vector Time
deriving instance U.Unbox Time

------------------------------------------------------------------------
-- Sample Functions

offset :: Clock -> Sample -> Duration
offset clock sample = remoteTime sample `sub` localTime clock sample

localTime :: Clock -> Sample -> Time
localTime clock (t1,_,_,t4) = clockTime clock (t1 + (t4 - t1) `div` 2)

remoteTime :: Sample -> Time
remoteTime (_,t2,t3,_) = mid t2 t3

networkDelay :: Clock -> Sample -> Duration
networkDelay clock sample = total - server
  where
    total  = fromSeconds (fromDiff clock (roundtrip sample))
    server = serverDelay sample

roundtrip :: Sample -> ClockDiff
roundtrip (t1,_,_,t4) = t4 - t1

serverDelay :: Sample -> Duration
serverDelay (_,t2,t3,_) = t3 `sub` t2

quality :: Clock -> Server -> Sample -> Double
quality clock svr s = exp (-x*x)
  where
    x = sampleError / baseError

    sampleError = currentError clock svr s
    baseError   = fromDiff clock (svrBaseError svr)

initialError :: Clock -> Server -> Sample -> Seconds
initialError clock Server{..} s = fromDiff clock (roundtrip s - svrMinRoundtrip)

-- TODO could be more accurate if we calculated estimated frequency error
currentError :: Clock -> Server -> Sample -> Seconds
currentError clock svr s = initial + drift * age
  where
    initial = initialError clock svr s

    age = fromDiff clock (currentTime - sampleTime)
    currentTime = time4 (U.head (svrRawSamples svr))
    sampleTime  = time4 s

    drift = 1 / 10000000 -- 0.1 ppm

------------------------------------------------------------------------
-- Vector/Statistics Utils

-- | Maps a functions over a list and turns it in to a vector.
vector :: (a -> Double) -> [a] -> U.Vector Double
vector f = U.fromList . map f

-- | Simple linear regression between 2 samples.
--   Takes two vectors Y={yi} and X={xi} and returns
--   (alpha, beta, r*r) such that Y = alpha + beta*X
linearRegression :: S.Sample -> S.WeightedSample -> (Double, Double)
linearRegression xs ys = (alpha, beta)
  where
    !c     = U.sum (U.zipWith (*) (U.map (subtract mx) xs) (U.map (subtract my . fst) ys)) / (n-1)
    !r     = c / (sx * sy)
    !mx    = S.mean xs
    !my    = S.meanWeighted ys
    !sx    = S.stdDev xs
    !sy    = sqrt (S.varianceWeighted ys)
    !n     = fromIntegral (U.length xs)
    !beta  = r * sy / sx
    !alpha = my - beta * mx
{-# INLINE linearRegression #-}

max_stars_repo_name: jacobstanley/ntpmon

max_issues_repo_name: jacobstanley/ntpmon

max_forks_repo_name: jacobstanley/ntpmon


text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}
-- |
-- Module    : Statistics.Types
-- Copyright : (c) 2009 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- Data types common used in statistics
module Statistics.Types
    ( -- * Confidence level
      CL
      -- ** Accessors
    , confidenceLevel
    , significanceLevel
      -- ** Constructors
    , mkCL
    , mkCLE
    , mkCLFromSignificance
    , mkCLFromSignificanceE
      -- ** Constants and conversion to nσ
    , cl90
    , cl95
    , cl99
      -- *** Normal approximation
    , nSigma
    , nSigma1
    , getNSigma
    , getNSigma1
      -- * p-value
    , PValue
      -- ** Accessors
    , pValue
      -- ** Constructors
    , mkPValue
    , mkPValueE
      -- * Estimates and upper/lower limits
    , Estimate(..)
    , NormalErr(..)
    , ConfInt(..)
    , UpperLimit(..)
    , LowerLimit(..)
      -- ** Constructors
    , estimateNormErr
    , (±)
    , estimateFromInterval
    , estimateFromErr
      -- ** Accessors
    , confidenceInterval
    , asymErrors
    , Scale(..)
      -- * Other
    , Sample
    , WeightedSample
    , Weights
    ) where

import Control.Monad                ((<=<), liftM2, liftM3)
import Control.DeepSeq              (NFData(..))
import Data.Aeson                   (FromJSON(..), ToJSON)
import Data.Binary                  (Binary(..))
import Data.Data                    (Data,Typeable)
import Data.Maybe                   (fromMaybe)
import Data.Vector.Unboxed          (Unbox)
import Data.Vector.Unboxed.Deriving (derivingUnbox)
import GHC.Generics                 (Generic)

#if __GLASGOW_HASKELL__ == 704
import qualified Data.Vector.Generic
import qualified Data.Vector.Generic.Mutable
#endif

import Statistics.Internal
import Statistics.Types.Internal
import Statistics.Distribution
import Statistics.Distribution.Normal

----------------------------------------------------------------
-- Data type for confidence level
----------------------------------------------------------------

-- |
-- Confidence level. In context of confidence intervals it's
-- probability of said interval covering true value of measured
-- value. In context of statistical tests it's @1-α@ where α is
-- significance of test.
--
-- Since confidence level are usually close to 1 they are stored as
-- @1-CL@ internally. There are two smart constructors for @CL@:
-- 'mkCL' and 'mkCLFromSignificance' (and corresponding variant
-- returning @Maybe@). First creates @CL@ from confidence level and
-- second from @1 - CL@ or significance level.
--
-- >>> cl95
-- mkCLFromSignificance 0.05
--
-- Prior to 0.14 confidence levels were passed to function as plain
-- @Doubles@. Use 'mkCL' to convert them to @CL@.
newtype CL a = CL a
               deriving (Eq, Typeable, Data, Generic)

instance Show a => Show (CL a) where
  showsPrec n (CL p) = defaultShow1 "mkCLFromSignificance" p n
instance (Num a, Ord a, Read a) => Read (CL a) where
  readPrec = defaultReadPrecM1 "mkCLFromSignificance" mkCLFromSignificanceE

instance (Binary a, Num a, Ord a) => Binary (CL a) where
  put (CL p) = put p
  get        = maybe (fail errMkCL) return . mkCLFromSignificanceE =<< get

instance (ToJSON a)                 => ToJSON   (CL a)
instance (FromJSON a, Num a, Ord a) => FromJSON (CL a) where
  parseJSON = maybe (fail errMkCL) return . mkCLFromSignificanceE <=< parseJSON

instance NFData   a => NFData   (CL a) where
  rnf (CL a) = rnf a

-- |
-- >>> cl95 > cl90
-- True
instance Ord a => Ord (CL a) where
  CL a <  CL b = a >  b
  CL a <= CL b = a >= b
  CL a >  CL b = a <  b
  CL a >= CL b = a <= b
  max (CL a) (CL b) = CL (min a b)
  min (CL a) (CL b) = CL (max a b)

-- | Create confidence level from probability β or probability
--   confidence interval contain true value of estimate. Will throw
--   exception if parameter is out of [0,1] range
--
-- >>> mkCL 0.95    -- same as cl95
-- mkCLFromSignificance 0.05
mkCL :: (Ord a, Num a) => a -> CL a
mkCL
  = fromMaybe (error "Statistics.Types.mkCL: probability is out if [0,1] range")
  . mkCLE

-- | Same as 'mkCL' but returns @Nothing@ instead of error if
--   parameter is out of [0,1] range
--
-- >>> mkCLE 0.95    -- same as cl95
-- Just (mkCLFromSignificance 0.05)
mkCLE :: (Ord a, Num a) => a -> Maybe (CL a)
mkCLE p
  | p >= 0 && p <= 1 = Just $ CL (1 - p)
  | otherwise        = Nothing

-- | Create confidence level from probability α or probability that
--   confidence interval does not contain true value of estimate. Will
--   throw exception if parameter is out of [0,1] range
--
-- >>> mkCLFromSignificance 0.05    -- same as cl95
-- mkCLFromSignificance 0.05
mkCLFromSignificance :: (Ord a, Num a) => a -> CL a
mkCLFromSignificance = fromMaybe (error errMkCL) . mkCLFromSignificanceE

-- | Same as 'mkCLFromSignificance' but returns @Nothing@ instead of error if
--   parameter is out of [0,1] range
--
-- >>> mkCLFromSignificanceE 0.05    -- same as cl95
-- Just (mkCLFromSignificance 0.05)
mkCLFromSignificanceE :: (Ord a, Num a) => a -> Maybe (CL a)
mkCLFromSignificanceE p
  | p >= 0 && p <= 1 = Just $ CL p
  | otherwise        = Nothing

errMkCL :: String
errMkCL = "Statistics.Types.mkPValCL: probability is out if [0,1] range"

-- | Get confidence level. This function is subject to rounding
--   errors. If @1 - CL@ is needed use 'significanceLevel' instead
confidenceLevel :: (Num a) => CL a -> a
confidenceLevel (CL p) = 1 - p

-- | Get significance level.
significanceLevel :: CL a -> a
significanceLevel (CL p) = p

-- | 90% confidence level
cl90 :: Fractional a => CL a
cl90 = CL 0.10

-- | 95% confidence level
cl95 :: Fractional a => CL a
cl95 = CL 0.05

-- | 99% confidence level
cl99 :: Fractional a => CL a
cl99 = CL 0.01

----------------------------------------------------------------
-- Data type for p-value
----------------------------------------------------------------

-- | Newtype wrapper for p-value.
newtype PValue a = PValue a
               deriving (Eq,Ord, Typeable, Data, Generic)

instance Show a => Show (PValue a) where
  showsPrec n (PValue p) = defaultShow1 "mkPValue" p n
instance (Num a, Ord a, Read a) => Read (PValue a) where
  readPrec = defaultReadPrecM1 "mkPValue" mkPValueE

instance (Binary a, Num a, Ord a) => Binary (PValue a) where
  put (PValue p) = put p
  get            = maybe (fail errMkPValue) return . mkPValueE =<< get

instance (ToJSON a)                 => ToJSON   (PValue a)
instance (FromJSON a, Num a, Ord a) => FromJSON (PValue a) where
  parseJSON = maybe (fail errMkPValue) return . mkPValueE <=< parseJSON

instance NFData a => NFData (PValue a) where
  rnf (PValue a) = rnf a

-- | Construct PValue. Throws error if argument is out of [0,1] range.
--
mkPValue :: (Ord a, Num a) => a -> PValue a
mkPValue = fromMaybe (error errMkPValue) . mkPValueE

-- | Construct PValue. Returns @Nothing@ if argument is out of [0,1] range.
mkPValueE :: (Ord a, Num a) => a -> Maybe (PValue a)
mkPValueE p
  | p >= 0 && p <= 1 = Just $ PValue p
  | otherwise        = Nothing

-- | Get p-value
pValue :: PValue a -> a
pValue (PValue p) = p

-- | P-value expressed in sigma. This is convention widely used in
--   experimental physics. N sigma confidence level corresponds to
--   probability within N sigma of normal distribution.
--
--   Note that this correspondence is for normal distribution. Other
--   distribution will have different dependency. Also experimental
--   distribution usually only approximately normal (especially at
--   extreme tails).
nSigma :: Double -> PValue Double
nSigma n
  | n > 0     = PValue $ 2 * cumulative standard (-n)
  | otherwise = error "Statistics.Extra.Error.nSigma: non-positive number of sigma"

-- | P-value expressed in sigma for one-tail hypothesis. This correspond to
--   probability of obtaining value less than @N·σ@.
nSigma1 :: Double -> PValue Double
nSigma1 n
  | n > 0     = PValue $ cumulative standard (-n)
  | otherwise = error "Statistics.Extra.Error.nSigma1: non-positive number of sigma"

-- | Express confidence level in sigmas
getNSigma :: PValue Double -> Double
getNSigma (PValue p) = negate $ quantile standard (p / 2)

-- | Express confidence level in sigmas for one-tailed hypothesis.
getNSigma1 :: PValue Double -> Double
getNSigma1 (PValue p) = negate $ quantile standard p

errMkPValue :: String
errMkPValue = "Statistics.Types.mkPValue: probability is out if [0,1] range"

----------------------------------------------------------------
-- Point estimates
----------------------------------------------------------------

-- |
-- A point estimate and its confidence interval. It's parametrized by
-- both error type @e@ and value type @a@. This module provides two
-- types of error: 'NormalErr' for normally distributed errors and
-- 'ConfInt' for error with normal distribution. See their
-- documentation for more details.
--
-- For example @144 ± 5@ (assuming normality) could be expressed as
--
-- > Estimate { estPoint = 144
-- >          , estError = NormalErr 5
-- >          }
--
-- Or if we want to express @144 + 6 - 4@ at CL95 we could write:
--
-- > Estimate { estPoint = 144
-- >          , estError = ConfInt
-- >                       { confIntLDX = 4
-- >                       , confIntUDX = 6
-- >                       , confIntCL  = cl95
-- >                       }
--
-- Prior to statistics 0.14 @Estimate@ data type used following definition:
--
-- > data Estimate = Estimate {
-- >      estPoint           :: {-# UNPACK #-} !Double
-- >    , estLowerBound      :: {-# UNPACK #-} !Double
-- >    , estUpperBound      :: {-# UNPACK #-} !Double
-- >    , estConfidenceLevel :: {-# UNPACK #-} !Double
-- >    }
--
-- Now type @Estimate ConfInt Double@ should be used instead. Function
-- 'estimateFromInterval' allow to easily construct estimate from same inputs.
data Estimate e a = Estimate
    { estPoint           :: !a
      -- ^ Point estimate.
    , estError           :: !(e a)
      -- ^ Confidence interval for estimate.
    } deriving (Eq, Read, Show, Generic
#if __GLASGOW_HASKELL__ >= 708
               , Typeable, Data
#endif
               )

instance (Binary   (e a), Binary   a) => Binary   (Estimate e a) where
  get = liftM2 Estimate get get
  put (Estimate ep ee) = put ep >> put ee
instance (FromJSON (e a), FromJSON a) => FromJSON (Estimate e a)
instance (ToJSON   (e a), ToJSON   a) => ToJSON   (Estimate e a)
instance (NFData   (e a), NFData   a) => NFData   (Estimate e a) where
    rnf (Estimate x dx) = rnf x `seq` rnf dx

-- |
-- Normal errors. They are stored as 1σ errors which corresponds to
-- 68.8% CL. Since we can recalculate them to any confidence level if
-- needed we don't store it.
newtype NormalErr a = NormalErr
  { normalError :: a
  }
  deriving (Eq, Read, Show, Typeable, Data, Generic)

instance Binary   a => Binary   (NormalErr a) where
  get = fmap NormalErr get
  put = put . normalError
instance FromJSON a => FromJSON (NormalErr a)
instance ToJSON   a => ToJSON   (NormalErr a)
instance NFData   a => NFData   (NormalErr a) where
    rnf (NormalErr x) = rnf x

-- | Confidence interval. It assumes that confidence interval forms
--   single interval and isn't set of disjoint intervals.
data ConfInt a = ConfInt
  { confIntLDX :: !a
    -- ^ Lower error estimate, or distance between point estimate and
    --   lower bound of confidence interval.
  , confIntUDX :: !a
    -- ^ Upper error estimate, or distance between point estimate and
    --   upper bound of confidence interval.
  , confIntCL  :: !(CL Double)
    -- ^ Confidence level corresponding to given confidence interval.
  }
  deriving (Read,Show,Eq,Typeable,Data,Generic)

instance Binary   a => Binary   (ConfInt a) where
  get = liftM3 ConfInt get get get
  put (ConfInt l u cl) = put l >> put u >> put cl
instance FromJSON a => FromJSON (ConfInt a)
instance ToJSON   a => ToJSON   (ConfInt a)
instance NFData   a => NFData   (ConfInt a) where
    rnf (ConfInt x y _) = rnf x `seq` rnf y

----------------------------------------
-- Constructors

-- | Create estimate with normal errors
estimateNormErr :: a            -- ^ Point estimate
                -> a            -- ^ 1σ error
                -> Estimate NormalErr a
estimateNormErr x dx = Estimate x (NormalErr dx)

-- | Synonym for 'estimateNormErr'
(±) :: a      -- ^ Point estimate
    -> a      -- ^ 1σ error
    -> Estimate NormalErr a
(±) = estimateNormErr

-- | Create estimate with asymmetric error.
estimateFromErr
  :: a                     -- ^ Central estimate
  -> (a,a)                 -- ^ Lower and upper errors. Both should be
                           --   positive but it's not checked.
  -> CL Double             -- ^ Confidence level for interval
  -> Estimate ConfInt a
estimateFromErr x (ldx,udx) cl = Estimate x (ConfInt ldx udx cl)

-- | Create estimate with asymmetric error.
estimateFromInterval
  :: Num a
  => a                     -- ^ Point estimate. Should lie within
                           --   interval but it's not checked.
  -> (a,a)                 -- ^ Lower and upper bounds of interval
  -> CL Double             -- ^ Confidence level for interval
  -> Estimate ConfInt a
estimateFromInterval x (lx,ux) cl
  = Estimate x (ConfInt (x-lx) (ux-x) cl)

----------------------------------------
-- Accessors

-- | Get confidence interval
confidenceInterval :: Num a => Estimate ConfInt a -> (a,a)
confidenceInterval (Estimate x (ConfInt ldx udx _))
  = (x - ldx, x + udx)

-- | Get asymmetric errors
asymErrors :: Estimate ConfInt a -> (a,a)
asymErrors (Estimate _ (ConfInt ldx udx _)) = (ldx,udx)

-- | Data types which could be multiplied by constant.
class Scale e where
  scale :: (Ord a, Num a) => a -> e a -> e a

instance Scale NormalErr where
  scale a (NormalErr e) = NormalErr (abs a * e)

instance Scale ConfInt where
  scale a (ConfInt l u cl) | a >= 0    = ConfInt  (a*l)  (a*u) cl
                           | otherwise = ConfInt (-a*u) (-a*l) cl

instance Scale e => Scale (Estimate e) where
  scale a (Estimate x dx) = Estimate (a*x) (scale a dx)

----------------------------------------------------------------
-- Upper/lower limit
----------------------------------------------------------------

-- | Upper limit. They are usually given for small non-negative values
--   when it's not possible detect difference from zero.
data UpperLimit a = UpperLimit
    { upperLimit        :: !a
      -- ^ Upper limit
    , ulConfidenceLevel :: !(CL Double)
      -- ^ Confidence level for which limit was calculated
    } deriving (Eq, Read, Show, Typeable, Data, Generic)

instance Binary   a => Binary   (UpperLimit a) where
  get = liftM2 UpperLimit get get
  put (UpperLimit l cl) = put l >> put cl
instance FromJSON a => FromJSON (UpperLimit a)
instance ToJSON   a => ToJSON   (UpperLimit a)
instance NFData   a => NFData   (UpperLimit a) where
    rnf (UpperLimit x cl) = rnf x `seq` rnf cl

-- | Lower limit. They are usually given for large quantities when
--   it's not possible to measure them. For example: proton half-life
data LowerLimit a = LowerLimit {
    lowerLimit        :: !a
    -- ^ Lower limit
  , llConfidenceLevel :: !(CL Double)
    -- ^ Confidence level for which limit was calculated
  } deriving (Eq, Read, Show, Typeable, Data, Generic)

instance Binary   a => Binary   (LowerLimit a) where
  get = liftM2 LowerLimit get get
  put (LowerLimit l cl) = put l >> put cl
instance FromJSON a => FromJSON (LowerLimit a)
instance ToJSON   a => ToJSON   (LowerLimit a)
instance NFData   a => NFData   (LowerLimit a) where
    rnf (LowerLimit x cl) = rnf x `seq` rnf cl

----------------------------------------------------------------
-- Deriving unbox instances
----------------------------------------------------------------

derivingUnbox "CL"
  [t| forall a. Unbox a => CL a -> a |]
  [| \(CL a) -> a |]
  [| CL           |]

derivingUnbox "PValue"
  [t| forall a. Unbox a => PValue a -> a |]
  [| \(PValue a) -> a |]
  [| PValue           |]

derivingUnbox "Estimate"
  [t| forall a e. (Unbox a, Unbox (e a)) => Estimate e a -> (a, e a) |]
  [| \(Estimate x dx) -> (x,dx) |]
  [| \(x,dx) -> (Estimate x dx) |]

derivingUnbox "NormalErr"
  [t| forall a. Unbox a => NormalErr a -> a |]
  [| \(NormalErr a) -> a |]
  [| NormalErr           |]

derivingUnbox "ConfInt"
  [t| forall a. Unbox a => ConfInt a -> (a, a, CL Double) |]
  [| \(ConfInt a b c) -> (a,b,c) |]
  [| \(a,b,c) -> ConfInt a b c   |]

derivingUnbox "UpperLimit"
  [t| forall a. Unbox a => UpperLimit a -> (a, CL Double) |]
  [| \(UpperLimit a b) -> (a,b) |]
  [| \(a,b) -> UpperLimit a b   |]

derivingUnbox "LowerLimit"
  [t| forall a. Unbox a => LowerLimit a -> (a, CL Double) |]
  [| \(LowerLimit a b) -> (a,b) |]
  [| \(a,b) -> LowerLimit a b   |]

ext: hs


text:
-----------------------------------------------------------------------------
-- |
-- Module      :  ForSyDe.Atom.Skel.Vector.DSP
-- Copyright   :  (c) George Ungureanu, KTH/EECS/ESY 2019-2020
-- License     :  BSD-style (see the file LICENSE)
--
-- Maintainer  :  ugeorge@kth.se
-- Stability   :  experimental
-- Portability :  portable
--
-- This module exports a set of 'Vector' patterns commonly used in signal processing
-- designs.
-----------------------------------------------------------------------------
module ForSyDe.Atom.Skel.Vector.DSP where

import Data.Complex
import qualified Data.Number.FixedFunctions as F
import ForSyDe.Atom.MoC as MoC
import ForSyDe.Atom.Skel.Vector as V hiding (duals, unduals)
import ForSyDe.Atom.Skel.Vector.Matrix as M
import ForSyDe.Atom ((><))

-- | Return the Taylor window. The Taylor window allows for a selectable sidelobe
-- suppression with a minimum broadening. This window is commonly used in radar
-- processing [1].  Code inspired from a <https://github.com/scipy/scipy/pull/8032 pull request for SciPy>.
--
-- Reference:
--
-- #ref1# [1] W. Carrara, R. Goodman, and R. Majewski "Spotlight Synthetic Aperture
-- Radar: Signal Processing Algorithms" Pages 512-513, July 1995.
taylor :: (Eq a, Floating a)
       => Int
       -- ^ number of points in the output window.
       -> Int
       -- ^ Number of nearly constant level sidelobes adjacent to the mainlobe
       -> a
       -- ^ Desired peak sidelobe level in decibels (db) relative to the mainlobe
       -> Vector a
       -- ^ The window, with the center value normalized to one (the value one appears
       -- only if the number of samples is odd).
taylor n nbar level = V.farm11 (*scale) w
  where
    -- explicit conversions to floating
    bN   = fromIntegral n
    nBar = fromIntegral nbar
    ma   = map fromIntegral [1 .. nbar-1]
    -- calculate intermediate values
    b  = 10**((-level) / 20)
    a  = log(b + sqrt(b**2 - 1)) / pi
    s2 = nBar ** 2 / (a**2 + (nBar - 0.5)**2)
    -- functions for calculating coefficients
    fmcalc m  = let numer = (-1)**(m+1) * product[1.0 - m**2/s2/(a**2 + (j - 0.5)**2) | j <- ma]
                    denom = 2 * product[1 - m**2/j**2 | j <- ma, j /= m]
                in numer / denom
    ccalc m x = cos(2 * pi * x * (m - bN/2 + 1/2) / bN)
    wcalc m   = 2 * dotvv (vector $ map fmcalc ma) (vector $ map (ccalc m) ma) + 1
    -- calculate window coefficients
    w  = vector $ map (wcalc . fromIntegral) [0..n-1]
    -- normalize (Note that this is not described in the original text [1])
    scale = 1 / wcalc (bN - 1) / 2

-- (*^*) :: (RealFrac a, RealFrac b) => a -> b -> a
-- x *^* y = realToFrac $ realToFrac x ** realToFrac y

-- | Returns a Taylor window with default arguments: 4 sidelobes, and peak sidelobe level of -30dB
taylor' n = taylor 4 (-30)

-- | Calculates the dot product between two vectors.
dotvv :: Num a => Vector a -> Vector a -> a
dotvv a b
  | V.length a == V.length b = V.reduce (+) (V.farm21 (*) a b)
  | otherwise                = error "Vector sizes must match"

-- | Calculates the dot product between a vector and a matrix
dotvm :: Num a => Vector a -> Vector (Vector a) -> Vector a
dotvm = dotvm' (+) (*)

-- | Higher order version of 'dotvv'. Applies dot product on vectors of /structured/
-- /types/, e.g. 'Signal's. See 'dotvv'.
dotvv' :: Num a
     => ((a -> a -> a) -> f a -> f a -> f a)
     -- ^ higher-order function that can wrap the (*) operation.
     -> Vector (f a) -> Vector (f a) -> f a
dotvv' wrap a b
  | V.length a == V.length b = V.reduce (wrap (+)) $ V.farm21 (wrap (*)) a b
  | otherwise                = error "Vector sizes must match"

-- | Higher order version of 'dotvm'. Implements the template for a dot operation
-- between a vector and a matrix.
--
-- >>> let mA = vector [vector[1,-1,1],  vector[1,-1,-1],  vector[1,1,-1],  vector[1,1,1]]
-- >>> let y  = vector[1,0,0,0]
-- >>> dotVecMat (+) (*) mA y
-- <1,1,1,1>
dotvm' :: (b -> b -> b)     -- ^ kernel function for a row/column reduction, e.g. @(+)@ for dot product
       -> (a -> b -> b)     -- ^ binary operation for pair-wise elements, e.g. @(*)@ for dot product
       -> Vector a          -- ^ /length/ = @xa@
       -> Vector (Vector b) -- ^ /size/ = @(xa,ya)@
       -> Vector b          -- ^ /length/ = @ya@
dotvm' f g vs = V.reduce (V.farm21 f) . V.farm21 (\x -> V.farm11 (g x)) vs

-- | Higher order pattern implementing the template for a dot operation between a
-- vector and a matrix.
--
-- >>> let mA = vector [vector[1,-1,1],  vector[1,-1,-1],  vector[1,1,-1],  vector[1,1,1]]
-- >>> let y  = vector[1,0,0,0]
-- >>> dotVecMat (+) (*) mA y
-- <1,1,1,1>
dotmv' :: (a -> a -> a)
       -- ^ kernel function for a row/column reduction, e.g. @(+)@ for dot product
       -> (b -> a -> a)
       -- ^ binary operation for pair-wise elements, e.g. @(*)@ for dot product
       -> Vector (Vector b) -- ^ /size/ = @(xa,ya)@
       -> Vector a      -- ^ /length/ = @xa@
       -> Vector a      -- ^ /length/ = @xa@
dotmv' f g mA y = V.farm11 (\x -> V.reduce f $ V.farm21 g x y) mA

-- | Compute a Hanning window.
--
-- Inspired from <https://hackage.haskell.org/package/sdr-0.1.0.6/src/hs_sources/SDR/FilterDesign.hs>
hanning :: (Floating n)
        => Int -- ^ The length of the window
        -> Vector n
hanning size = V.farm11 func $ V.vector [0..size-1]
  where
    func idx = let i = fromIntegral idx
                   n = fromIntegral size
               in 0.5 * (1 - cos((2 * pi * i) / (n - 1)))

-- | Compute a Hamming window.
--
-- Inspired from <https://hackage.haskell.org/package/sdr-0.1.0.6/src/hs_sources/SDR/FilterDesign.hs>
hamming :: (Floating n)
        => Int -- ^ The length of the window
        -> Vector n
hamming size = V.farm11 func $ V.vector [0..size-1]
  where
    func idx = let i = fromIntegral idx
                   n = fromIntegral size
               in 0.54 - 0.46 * cos((2 * pi * i) / (n - 1))

-- | Compute a Blackman window.
--
-- Inspired from <https://hackage.haskell.org/package/sdr-0.1.0.6/src/hs_sources/SDR/FilterDesign.hs>
blackman :: (Floating n)
        => Int -- ^ The length of the window
        -> Vector n
blackman size = V.farm11 func $ V.vector [0..size-1]
  where
    func idx = let i = fromIntegral idx
                   n = fromIntegral size
               in 0.42 - 0.5 * cos((2 * pi * i) / (n - 1)) + 0.08 * cos((4 * pi * i) / (n - 1))

-- | Moving average filter (FIR) on numbers, applied in reverse order (more
-- optimized).
--
-- >>> let v = vector [0,0,0,0,0,1]
-- >>> let c = vector [1,2,3]
-- >>> fir c
-- <0,0,0,3,2,1>
fir :: Num a
    => Vector a  -- ^ vector of coefficients
    -> Vector a  -- ^ input vector of numbers; /size/ = @n@
    -> Vector a  -- ^ output vector of numbers; /size/ = @n@
fir coefs = V.reverse . V.farm11 applyFilter . tails . V.reverse
  where
    applyFilter = V.reduce (+) . V.farm21 (*) coefs

-- | Higher order version of 'fir'. Can create a systolic FIR process network. See
-- <ForSyDe-Atom.html#ungureanu20a [Ungureanu20a]> for a discussion on the relation
-- between 'fir' and 'fir''.
--
-- >>> let c = vector [1,2,3]
-- >>> let s = SY.signal [1,0,0,0,0,0,0,0]
-- >>> fir' (SY.comb21 (+)) (\c -> SY.comb11 (*c)) (SY.delay 0) c s
-- {3,2,1,0,0,0,0,0}
fir' :: (a -> a -> a)  -- ^ process/operation replacing '+'
     -> (c -> a -> a)  -- ^ process/operation replacing '*'
     -> (a -> a)       -- ^ delay process
     -> Vector c       -- ^ vector of coefficients
     -> a              -- ^ input signal/structure
     -> a              -- ^ output signal/structure
fir' plus times delay coefs =
  -- V.reduce plus . V.farm21 (\c -> times c) coefs . V.recur (V.fanoutn n delay <: id)
  V.reduce plus . V.farm21 times coefs . V.recuri (V.fanoutn n delay)
  where n = V.length coefs - 1

-- | generates the "twiddle" coefficients for a FFT network.
twiddles :: Floating a => Int -> V.Vector (Complex a)
twiddles bN = (bitrev . V.take (bN `div` 2)) (V.farm11 bW $ vector [0..])
  where bW x = (cis . negate) (-2 * pi * fromIntegral x / fromIntegral bN)

-- | Radix-2 decimation-in-frequecy Fast Fourier Transform, applied on numbers.  For
-- the difference between DIT- and DIF-FFT, see
-- <https://www.slideshare.net/chappidi_saritha/decimation-in-time-and-frequency>
fft :: RealFloat a
    => Int -- ^ number of FFT stages (@== log_2@ of /length/ of input vector)
    -> V.Vector (Complex a) -> V.Vector (Complex a)
fft k vs | n == 2^k = V.reverse $ bitrev $ (stage `V.pipe1` V.iterate k (*2) 2) vs
         | otherwise = error $ "n=" ++ show n ++ "; k=" ++ show k
  where
    stage   w = V.concat . V.farm21 segment (twiddles n) . V.group w
    segment t = (><) unduals . (><) (V.farm22 (butterfly t)) . duals
    n         = V.length vs        -- length of input
    -------------------------------------------------
    butterfly w x0 x1 = let t = w * x1 in (x0 + t, x0 - t) -- kernel function

-- | Higher order version of 'fft'. Takes the "butterfly" function as argument.
--
-- > butterfly w x0 x1 = let t = w * x1 in (x0 + t, x0 - t)
-- > fft' butterfly === fft
fft' :: Floating a
     => (Complex a -> a -> a -> (a, a)) -- ^ "butterfly" function.
     -> Int  -- ^ number of FFT stages (@== log_2@ of /length/ of input vector)
     -> V.Vector a -> V.Vector a
fft' butterfly k vs | n == 2^k = bitrev $ (stage `V.pipe1` (V.iterate k (*2) 2)) vs
  where
    stage   w = V.concat . V.farm21 segment (twiddles n) . V.group w
    segment t = (><) unduals . (><) (V.farm22 (butterfly t)) . duals
    n         = V.length vs        -- length of input

-- | splits a vector in two equal parts.
--
-- >>> duals $ vector [1,2,3,4,5,6,7]
-- (<1,2,3>,<4,5,6>)
duals    :: Vector a -> (Vector a, Vector a)
duals v  = (V.take k v, V.drop k v)
  where k = V.length v `div` 2

-- | concatenates a previously split vector. See also 'duals'
unduals  :: Vector a -> Vector a -> Vector a
unduals x y =  x <++> y

-- | performs a bit-reverse permutation.
--
-- <<fig/skel-vector-comm-bitrev.png>>
-- <<fig/skel-vector-comm-bitrev-net.png>>
--
-- >>> bitrev $ vector ["000","001","010","011","100","101","110","111"]
-- <"111","011","101","001","110","010","100","000">
bitrev (x:>Null) = V.unit x
bitrev xs        = bitrev (evens xs) <++> bitrev (odds xs)

-- duals   v = let k = length v `div` 2
--             in  S.farm22 (,) (take k v) (drop k v)

-- unduals = (<++>)

-- fft n = V.vector . map (fmap realToFrac) . Sh.fromVector . ShDSP.fft (2^n) . Sh.vector . map (fmap realToFrac) . V.fromVector

max_stars_repo_name: forsyde/forsyde-atom

max_issues_repo_name: forsyde/forsyde-atom

max_forks_repo_name: forsyde/forsyde-atom


text:
{-# language TemplateHaskell #-}
{-# language FlexibleInstances #-}
{-# language GADTs #-}
{-# language TypeOperators #-}

{-# OPTIONS_GHC -fno-warn-orphans #-}

module Data.Eq.StableName.Instances ( ) where

import Control.Applicative (liftA2, ZipList)
import Control.Concurrent.Chan (Chan)
import Control.Concurrent.MVar (MVar)
import Control.Exception
import Data.Char (GeneralCategory)
import Data.Complex (Complex)
import Data.Data (TyCon, Fixity, Constr, Proxy, ConstrRep, DataRep)
import Data.Eq.StableName.Class (StableNameEq, (==$), ($==))
import Data.Eq.StableName.TH (deriveStableNameEq)
import Data.Fixed (Fixed)
import Data.Function (on)
import Data.Functor.Compose (Compose(Compose), getCompose)
import Data.Functor.Const (Const)
import Data.Functor.Identity (Identity)
import Data.Functor.Product (Product(Pair))
import Data.Functor.Sum (Sum(InL, InR))
import Data.IORef (IORef)
import Data.Int (Int8, Int16, Int32, Int64)
import Data.List.NonEmpty (NonEmpty)
import Data.Ord (Down)
import Data.Ratio (Ratio)
import Data.STRef (STRef)
import Data.Type.Coercion (Coercion)
import Data.Type.Equality ((:~:), (:~~:))
import Data.Unique (Unique)
import Data.Version (Version)
import Data.Void (Void)
import Data.Word (Word8, Word16, Word32, Word64)
import Foreign.C.Error (Errno)
import Foreign.C.Types
import Foreign.Ptr (IntPtr, WordPtr, Ptr, FunPtr)
import GHC.Conc (BlockReason, ThreadId, TVar, ThreadStatus)
import GHC.Event (Lifetime, TimeoutKey, FdKey, Event)
import GHC.ExecutionStack (SrcLoc)
import GHC.Exts (SpecConstrAnnotation)
import GHC.Fingerprint (Fingerprint)
import GHC.ForeignPtr (ForeignPtr)
import GHC.Generics
import GHC.IO.Buffer (BufferState)
import GHC.IO.Device (SeekMode, IODeviceType)
import GHC.IO.Encoding.Types (CodingProgress)
import GHC.IO.Handle (Newline, Handle, NewlineMode, BufferMode, HandlePosn)
import GHC.Stable (StablePtr)
import GHC.TypeLits (SomeSymbol)
import GHC.TypeNats (SomeNat)
import GHC.Types (Module, TrName)
import Numeric.Natural (Natural)
import System.Exit (ExitCode)
import System.IO (IOMode)
import System.IO.Error (IOErrorType)
import System.Mem.StableName (StableName)
import System.Posix.Types
import Text.Read.Lex (Number, Lexeme)
import Type.Reflection (SomeTypeRep, TypeRep)
import qualified Data.Monoid
import qualified Data.Semigroup

instance StableNameEq Bool                        where (==$) = eq
instance StableNameEq Ordering                    where (==$) = eq
instance StableNameEq Char                        where (==$) = eq
instance StableNameEq Float                       where (==$) = eq
instance StableNameEq Double                      where (==$) = eq
instance StableNameEq Int                         where (==$) = eq
instance StableNameEq Int8                        where (==$) = eq
instance StableNameEq Int16                       where (==$) = eq
instance StableNameEq Int32                       where (==$) = eq
instance StableNameEq Int64                       where (==$) = eq
instance StableNameEq Integer                     where (==$) = eq
instance StableNameEq Natural                     where (==$) = eq
instance StableNameEq Word                        where (==$) = eq
instance StableNameEq Word8                       where (==$) = eq
instance StableNameEq Word16                      where (==$) = eq
instance StableNameEq Word32                      where (==$) = eq
instance StableNameEq Word64                      where (==$) = eq
instance StableNameEq SomeTypeRep                 where (==$) = eq
instance StableNameEq TyCon                       where (==$) = eq
instance StableNameEq Module                      where (==$) = eq
instance StableNameEq TrName                      where (==$) = eq
instance StableNameEq GeneralCategory             where (==$) = eq
instance StableNameEq IOMode                      where (==$) = eq
instance StableNameEq IntPtr                      where (==$) = eq
instance StableNameEq WordPtr                     where (==$) = eq
instance StableNameEq CUIntMax                    where (==$) = eq
instance StableNameEq CIntMax                     where (==$) = eq
instance StableNameEq CUIntPtr                    where (==$) = eq
instance StableNameEq CIntPtr                     where (==$) = eq
instance StableNameEq CSUSeconds                  where (==$) = eq
instance StableNameEq CUSeconds                   where (==$) = eq
instance StableNameEq CTime                       where (==$) = eq
instance StableNameEq CClock                      where (==$) = eq
instance StableNameEq CSigAtomic                  where (==$) = eq
instance StableNameEq CWchar                      where (==$) = eq
instance StableNameEq CSize                       where (==$) = eq
instance StableNameEq CPtrdiff                    where (==$) = eq
instance StableNameEq CDouble                     where (==$) = eq
instance StableNameEq CFloat                      where (==$) = eq
instance StableNameEq CBool                       where (==$) = eq
instance StableNameEq CULLong                     where (==$) = eq
instance StableNameEq CLLong                      where (==$) = eq
instance StableNameEq CULong                      where (==$) = eq
instance StableNameEq CLong                       where (==$) = eq
instance StableNameEq CUInt                       where (==$) = eq
instance StableNameEq CInt                        where (==$) = eq
instance StableNameEq CUShort                     where (==$) = eq
instance StableNameEq CShort                      where (==$) = eq
instance StableNameEq CUChar                      where (==$) = eq
instance StableNameEq CSChar                      where (==$) = eq
instance StableNameEq CChar                       where (==$) = eq
instance StableNameEq SomeNat                     where (==$) = eq
instance StableNameEq SomeSymbol                  where (==$) = eq
instance StableNameEq DecidedStrictness           where (==$) = eq
instance StableNameEq SourceStrictness            where (==$) = eq
instance StableNameEq SourceUnpackedness          where (==$) = eq
instance StableNameEq Associativity               where (==$) = eq
instance StableNameEq Data.Monoid.Any             where (==$) = eq
instance StableNameEq Data.Monoid.All             where (==$) = eq
instance StableNameEq ArithException              where (==$) = eq
instance StableNameEq IOException                 where (==$) = eq
instance StableNameEq MaskingState                where (==$) = eq
instance StableNameEq BufferState                 where (==$) = eq
instance StableNameEq CodingProgress              where (==$) = eq
instance StableNameEq SeekMode                    where (==$) = eq
instance StableNameEq IODeviceType                where (==$) = eq
instance StableNameEq Newline                     where (==$) = eq
instance StableNameEq Handle                      where (==$) = eq
instance StableNameEq IOErrorType                 where (==$) = eq
instance StableNameEq AsyncException              where (==$) = eq
instance StableNameEq ArrayException              where (==$) = eq
instance StableNameEq Errno                       where (==$) = eq
instance StableNameEq Fd                          where (==$) = eq
instance StableNameEq CTimer                      where (==$) = eq
instance StableNameEq CKey                        where (==$) = eq
instance StableNameEq CId                         where (==$) = eq
instance StableNameEq CFsFilCnt                   where (==$) = eq
instance StableNameEq CFsBlkCnt                   where (==$) = eq
instance StableNameEq CClockId                    where (==$) = eq
instance StableNameEq CBlkCnt                     where (==$) = eq
instance StableNameEq CBlkSize                    where (==$) = eq
instance StableNameEq CRLim                       where (==$) = eq
instance StableNameEq CTcflag                     where (==$) = eq
instance StableNameEq CSpeed                      where (==$) = eq
instance StableNameEq CCc                         where (==$) = eq
instance StableNameEq CUid                        where (==$) = eq
instance StableNameEq CNlink                      where (==$) = eq
instance StableNameEq CGid                        where (==$) = eq
instance StableNameEq CSsize                      where (==$) = eq
instance StableNameEq CPid                        where (==$) = eq
instance StableNameEq COff                        where (==$) = eq
instance StableNameEq CMode                       where (==$) = eq
instance StableNameEq CIno                        where (==$) = eq
instance StableNameEq CDev                        where (==$) = eq
instance StableNameEq Lifetime                    where (==$) = eq
instance StableNameEq Event                       where (==$) = eq
instance StableNameEq BlockReason                 where (==$) = eq
instance StableNameEq ThreadId                    where (==$) = eq
instance StableNameEq TimeoutKey                  where (==$) = eq
instance StableNameEq Unique                      where (==$) = eq
instance StableNameEq Data.Data.Fixity            where (==$) = eq
instance StableNameEq Constr                      where (==$) = eq
instance StableNameEq SpecConstrAnnotation        where (==$) = eq
instance StableNameEq Void                        where (==$) = eq
instance StableNameEq (StablePtr a)               where (==$) = eq
instance StableNameEq (Ptr a)                     where (==$) = eq
instance StableNameEq (FunPtr a)                  where (==$) = eq
instance StableNameEq (MVar a)                    where (==$) = eq
instance StableNameEq (IORef a)                   where (==$) = eq
instance StableNameEq (ForeignPtr a)              where (==$) = eq
instance StableNameEq (TVar a)                    where (==$) = eq
instance StableNameEq (StableName a)              where (==$) = eq
instance StableNameEq (Fixed a)                   where (==$) = eq
instance StableNameEq (V1 p)                      where (==$) = eq
instance StableNameEq (U1 p)                      where (==$) = eq
instance StableNameEq (Type.Reflection.TypeRep a) where (==$) = eq
instance StableNameEq (STRef s a)                 where (==$) = eq
instance StableNameEq (Proxy a)                   where (==$) = eq
instance StableNameEq FdKey                       where (==$) = eq
instance StableNameEq (URec Word p)               where (==$) = eq
instance StableNameEq (URec Int p)                where (==$) = eq
instance StableNameEq (URec Float p)              where (==$) = eq
instance StableNameEq (URec Double p)             where (==$) = eq
instance StableNameEq (URec Char p)               where (==$) = eq
instance StableNameEq (URec (Ptr ()) p)           where (==$) = eq
instance StableNameEq (a :~: b)                   where (==$) = eq
instance StableNameEq (a :~~: b)                  where (==$) = eq
instance StableNameEq (Coercion a b)              where (==$) = eq

eq :: (Eq a) => a -> a -> IO Bool
eq x y = pure (x == y)

instance StableNameEq (f p) => StableNameEq (Rec1 f p) where
    (==$) = ($==) `on` unRec1

instance StableNameEq c => StableNameEq (K1 i c p) where
    (==$) = ($==) `on` unK1

instance StableNameEq (f p) => StableNameEq (M1 i c f p) where
    (==$) = ($==) `on` unM1

instance StableNameEq (f (g p)) => StableNameEq ((:.:) f g p) where
    (==$) = ($==) `on` unComp1

instance (StableNameEq (f p), StableNameEq (g p)) => StableNameEq ((:+:) f g p) where
    L1 f1 ==$ L1 f2 = f1 $== f2
    R1 g1 ==$ R1 g2 = g1 $== g2
    _     ==$ _     = pure False

instance (StableNameEq (f p), StableNameEq (g p)) => StableNameEq ((:*:) f g p) where
    (f1 :*: g1) ==$ (f2 :*: g2) = liftA2 (&&) (f1 $== f2) (g1 $== g2)

instance (StableNameEq (f a), StableNameEq (g a)) => StableNameEq (Data.Functor.Sum.Sum f g a) where
    Data.Functor.Sum.InL f1 ==$ Data.Functor.Sum.InL f2 = f1 $== f2
    Data.Functor.Sum.InR g1 ==$ Data.Functor.Sum.InR g2 = g1 $== g2
    _                       ==$ _      = pure False

instance (StableNameEq (f a), StableNameEq (g a)) => StableNameEq (Product f g a) where
    Pair f1 g1 ==$ Pair f2 g2 = liftA2 (&&) (f1 $== f2) (g1 $== g2)

instance StableNameEq (f (g a)) => StableNameEq (Compose f g a) where
    (==$) = ($==) `on` getCompose

instance (StableNameEq a) => StableNameEq (Data.Semigroup.Arg a b) where
    Data.Semigroup.Arg a _ ==$ Data.Semigroup.Arg b _ = a $== b

concat <$> traverse deriveStableNameEq
  [ ''()
  , ''(,)
  , ''(,,)
  , ''(,,,)
  , ''(,,,,)
  , ''(,,,,,)
  , ''(,,,,,,)
  , ''(,,,,,,,)
  , ''(,,,,,,,,)
  , ''(,,,,,,,,,)
  , ''(,,,,,,,,,,)
  , ''(,,,,,,,,,,,)
  , ''(,,,,,,,,,,,,)
  , ''(,,,,,,,,,,,,,)
  , ''(,,,,,,,,,,,,,,)
  , ''[]
  , ''Maybe
  , ''SrcLoc
  , ''Number
  , ''Lexeme
  , ''Fingerprint
  , ''GHC.Generics.Fixity
  , ''ErrorCall
  , ''NewlineMode
  , ''BufferMode
  , ''ExitCode
  , ''ThreadStatus
  , ''HandlePosn
  , ''Version
  , ''ConstrRep
  , ''DataRep
  , ''Ratio
  , ''Par1
  , ''Down
  , ''Data.Monoid.Last
  , ''Data.Monoid.First
  , ''Data.Monoid.Product
  , ''Data.Monoid.Sum
  , ''Data.Monoid.Dual
  , ''Identity
  , ''ZipList
  , ''Chan
  , ''NonEmpty
  , ''Data.Semigroup.Option
  , ''Data.Semigroup.WrappedMonoid
  , ''Data.Semigroup.Last
  , ''Data.Semigroup.First
  , ''Data.Semigroup.Max
  , ''Data.Semigroup.Min
  , ''Complex
  , ''Either
  , ''Const
  ]

max_stars_repo_name: basvandijk/stable-name-eq

max_issues_repo_name: basvandijk/stable-name-eq

max_forks_repo_name: basvandijk/stable-name-eq


text:
module Perceptron where

import Prelude hiding ((<>))
import Numeric.LinearAlgebra
import qualified Data.Monoid
import System.IO.Unsafe

orInputs :: Matrix R
orInputs = (4><2)
  [0, 0,
   0, 1,
   1, 0,
   1, 1]

orTargets :: Matrix R
orTargets = col [0, 1, 1, 1]

orInitialWeights :: Matrix R
orInitialWeights = initialWeights (cols orInputs) 1

-- todo: set the initial weights correctly
initialWeights :: Int -> Int -> Matrix R
initialWeights r c = unsafePerformIO $ rand (r + 1) c

addBiasNodes :: Matrix R -> Matrix R
addBiasNodes m = let rs = rows m
                     bias = (rs><1) $ repeat (-1)
                  in m ||| bias

feed :: Matrix R -> Vector R -> Vector R
feed weights = flatten . fwd weights . asColumn

fwd :: Matrix R -> Matrix R -> Matrix R
fwd weights inputs = let withBias = addBiasNodes inputs
                      in step $ withBias <> weights

update :: Double -> Matrix R -> Matrix R -> Matrix R -> Matrix R
update learningRate weights targets inputs  =
  weights - scale learningRate (tr withBias <> (activations - targets))
    where withBias = addBiasNodes inputs
          activations = fwd weights inputs

max_stars_repo_name: mfandl/hml

max_issues_repo_name: mfandl/hml

max_forks_repo_name: mfandl/hml


text:
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  Graphics.Rendering.Plot.Render.Types
-- Copyright   :  (c) A. V. H. McPhail 2010
-- License     :  BSD3
--
-- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
-- Stability   :  provisional
-- Portability :  portable
--
-- Rendering 'Figure's
--
-----------------------------------------------------------------------------

module Graphics.Rendering.Plot.Render.Types where

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

--import Data.Either

--import Data.Packed.Vector
--import Numeric.LinearAlgebra.Linear

--import Data.Word

import Data.Maybe

import Data.Colour.SRGB
import Data.Colour.Names

--import qualified Data.Array.IArray as A

import qualified Graphics.Rendering.Cairo as C
import qualified Graphics.Rendering.Cairo.Matrix as CM
import qualified Graphics.Rendering.Pango as P

import Control.Monad.Reader
import Control.Monad.State
--import Control.Monad.Trans

import Graphics.Rendering.Plot.Types
import Graphics.Rendering.Plot.Defaults

--import Graphics.Rendering.Plot.Figure.Text

--import qualified Text.Printf as Printf

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

-----------------------------------------------------------------------------
{-
newtype Render a = FR { runRender :: StateT BoundingBox C.Render a }
    deriving(Monad, MonadState BoundingBox, MonadTrans (StateT BoundingBox))
-}

data RenderEnv = RenderEnv {
                  _pangocontext    :: P.PangoContext
                  , _renderoptions :: Options
                 }

newtype BoundedT m a = BT { runRender :: ReaderT RenderEnv (StateT BoundingBox m) a }
    deriving(Monad, MonadState BoundingBox, MonadReader RenderEnv)

instance MonadTrans BoundedT where
    lift m = BT $ lift $ lift m

type Render = BoundedT C.Render

evalRender :: Render a -> RenderEnv -> BoundingBox -> C.Render a
evalRender m r = evalStateT (runReaderT (runRender m) r)

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

cairo :: C.Render a -> Render a
cairo = lift

pango :: IO a -> C.Render a
pango = liftIO

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

bbX, bbY, bbW, bbH :: Render Double
bbX = gets _bbX
bbY = gets _bbY
bbW = gets _bbW
bbH = gets _bbH

bbLeftWidth, bbCentreWidth, bbRightWidth, bbBottomHeight, bbCentreHeight, bbTopHeight :: Render Double
bbLeftWidth    = gets $ \(BoundingBox x _ _ _) -> x
bbCentreWidth  = gets $ \(BoundingBox x _ w _) -> x + w / 2
bbRightWidth   = gets $ \(BoundingBox x _ w _) -> x + w
bbBottomHeight = gets $ \(BoundingBox _ y _ h) -> y + h
bbCentreHeight = gets $ \(BoundingBox _ y _ h) -> y + h / 2
bbTopHeight    = gets $ \(BoundingBox _ y _ _) -> y

bbShiftLeft, bbShiftRight, bbLowerTop, bbRaiseBottom :: Double -> Render ()
bbShiftLeft   n = modify $ \(BoundingBox x y w h) ->
                            BoundingBox (x+n) y     (w-n) h
bbShiftRight  n = modify $ \(BoundingBox x y w h) ->
                            BoundingBox x     y     (w-n) h
bbLowerTop    n = modify $ \(BoundingBox x y w h) ->
                            BoundingBox x     (y+n) w     (h-n)
bbRaiseBottom n = modify $ \(BoundingBox x y w h) ->
                            BoundingBox x     y     w     (h-n)

applyPads :: Padding -> Render ()
applyPads (Padding l r b t) = modify (\(BoundingBox x y w h) ->
                                      BoundingBox (x+l) (y+t) (w-l-r) (h-t-b))

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

clipBoundary :: Render ()
clipBoundary = do
  (BoundingBox x y w h) <- get
  cairo $ do
    C.rectangle x y w h
    C.clip

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

-- | output file type
data OutputType = PNG | PS | PDF | SVG

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

setColour :: Color -> C.Render ()
setColour c = let (RGB r g b) = toSRGB c
              in C.setSourceRGBA r g b 1 -- no transparent colours

setDashes :: [Dash] -> C.Render ()
setDashes [] = C.setDash [] 0
setDashes xs = do
  let xs' = map (\d -> case d of { Dot -> 1 ; Dash -> 3 }) xs
  C.setDash xs' 0

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

getDefaultTextOptions :: P.PangoContext -> IO TextOptions
getDefaultTextOptions pc = do
  fd <- P.contextGetFontDescription pc
  getTextOptionsFD fd

getTextOptionsFD :: P.FontDescription -> IO TextOptions
getTextOptionsFD fd = do
  ff' <- P.fontDescriptionGetFamily fd
  fs' <- P.fontDescriptionGetStyle fd
  fv' <- P.fontDescriptionGetVariant fd
  fw' <- P.fontDescriptionGetWeight fd
  fc' <- P.fontDescriptionGetStretch fd
  fz' <- P.fontDescriptionGetSize fd
  let ff = fromMaybe defaultFontFamily ff'
      fs = fromMaybe defaultFontStyle fs'
      fv = fromMaybe defaultFontVariant fv'
      fw = fromMaybe defaultFontWeight fw'
      fc = fromMaybe defaultFontStretch fc'
      fz = fromMaybe defaultFontSize fz'
  return $ TextOptions (FontOptions ff fs fv fw fc) fz black

setTextOptions :: TextOptions -> P.PangoLayout -> C.Render ()
setTextOptions to lo = do
  fd' <- pango $ P.layoutGetFontDescription lo
  fd <- case fd' of
         Nothing   -> pango $ P.fontDescriptionNew
         Just fd'' -> return fd''
  setTextOptionsFD to fd
  pango $ P.layoutSetFontDescription lo (Just fd)

setTextOptionsFD :: TextOptions -> P.FontDescription -> C.Render ()
setTextOptionsFD (TextOptions (FontOptions ff fs fv fw fc) fz c) fd = do
  pango $ do
    P.fontDescriptionSetFamily fd ff
    P.fontDescriptionSetStyle fd fs
    P.fontDescriptionSetVariant fd fv
    P.fontDescriptionSetWeight fd fw
    P.fontDescriptionSetStretch fd fc
    P.fontDescriptionSetSize fd fz
  setColour c

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

textPad :: Double
textPad = 2

data TextXAlign = TLeft | Centre | TRight
data TextYAlign = TBottom | Middle | TTop

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

setLineOptions :: LineOptions -> C.Render ()
setLineOptions (LineOptions ds lw) = do
  setDashes ds
  C.setLineWidth lw

setLineStyle :: LineType -> C.Render ()
setLineStyle NoLine          = return ()
setLineStyle (ColourLine c)  = setColour c
setLineStyle (TypeLine lo c) = do
  setLineOptions lo
  setColour c

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

setPointOptions :: PointOptions -> C.Render ()
setPointOptions (PointOptions pz c) = do
  setColour c
  C.scale pz pz

setPointStyle :: PointType -> C.Render Glyph
setPointStyle (FullPoint po g) = do
  setPointOptions po
  return g

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

flipVerticalMatrix :: CM.Matrix
flipVerticalMatrix = CM.Matrix 1 0 0 (-1) 0 0

flipVertical :: C.Render ()
flipVertical = C.transform flipVerticalMatrix

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

max_stars_repo_name: tsurucapital/plot

max_issues_repo_name: tsurucapital/plot

max_forks_repo_name: tsurucapital/plot


text:
module Statistics.CSV
    ( txCntInChainMemPoolToCSV
    , focusToCSV
    ) where

import           Control.Monad.Random (MonadRandom (..), evalRandT)
import qualified Data.Text as T
import           System.IO (hPutStrLn)
import           System.Random (mkStdGen)

import           Pos.Infra.Util.JsonLog.Events (JLMemPool (..), MemPoolModifyReason (..))
import           Statistics.Focus (Focus (..))
import           Types
import           Universum

txCntInChainMemPoolToCSV :: FilePath
                         -> Double
                         -> [(NodeId, Timestamp, Int)]
                         -> [(NodeId, Timestamp, JLMemPool)]
                         -> IO ()
txCntInChainMemPoolToCSV f sp txCnt mp =
    flip evalRandT (mkStdGen 918273) $ liftIO $ withFile f WriteMode $ \h -> do
        hPutStrLn h "time,txCount,txType,node"
        for_ txCnt $ \(n, ts, cnt) -> csvLine h "written" n ts (fromIntegral cnt)
        for_ mp $ \(n, ts, p@JLMemPool{..}) ->
            whenM (inSample jlmReason) $ do
                csvLine h (toTxType "Wait" p) n ts jlmWait
                csvLine h (toTxType "Modify" p) n ts jlmModify
                csvLine h (toTxType "SizeAfter" p) n ts (fromIntegral jlmSizeAfter)
  where
    csvLine :: MonadIO m => Handle -> String -> NodeId -> Timestamp -> Integer -> m ()
    csvLine h txType node time txCount = liftIO $ hPutStrLn h $ show (fromIntegral time :: Integer)
        ++ "," ++ show txCount ++ "," ++ txType ++ "," ++ toString node

    draw :: MonadRandom m => m Bool
    draw = (<= sp) <$> getRandomR (0, 1)

    inSample :: MonadRandom m => MemPoolModifyReason -> m Bool
    inSample ProcessTransaction = draw
    inSample _                      = return True

    toTxType :: String -> JLMemPool -> String
    toTxType s JLMemPool{..} =
        let reason = case jlmReason of
                ApplyBlock             -> "ApplyBlock"
                ApplyBlockWithRollback -> "ApplyBlockWithRollback"
                ProcessTransaction     -> "ProcessTransaction"
        in  "mp_" ++ reason ++ "_" ++ s

focusToCSV :: FilePath -> [(Timestamp, NodeId, Focus)] -> IO ()
focusToCSV f xs = withFile f WriteMode $ \h -> do
    hPutStrLn h "time,delta_first_seconds,delta_step_seconds,node,type,block/error"
    case xs of
        []                -> return ()
        ((ts0, _, _) : _) -> foldM_ (step h ts0) ts0 xs
  where
    step :: Handle -> Timestamp -> Timestamp -> (Timestamp, NodeId, Focus) -> IO Timestamp
    step h ts0 ts (ts', n, y) = do
        let dt0 = fromIntegral (ts' - ts0) / 1000000 :: Double
            dt  = fromIntegral (ts' - ts ) / 1000000 :: Double
        case y of
            Received me         -> csvLine h ts' dt0 dt n "received" $ maybe "" show me
            InCreatedBlock hash -> csvLine h ts' dt0 dt n "created" $ T.take 6 hash
            InAdoptedBlock hash -> csvLine h ts' dt0 dt n "adopted" $ T.take 6 hash
        return ts'

    csvLine :: Handle -> Timestamp -> Double -> Double -> NodeId -> String -> BlockHash -> IO ()
    csvLine h ts dt0 dt node t he =
        hPutStrLn h $ show ts ++ "," ++ show dt0 ++ "," ++ show dt ++ "," ++ show node ++ "," ++ t
            ++ "," ++ toString he

max_stars_repo_name: serokell/cardano-sl

max_issues_repo_name: serokell/cardano-sl

max_forks_repo_name: serokell/cardano-sl


text:
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE TupleSections #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE TypeSynonymInstances #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE CPP #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE DefaultSignatures #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE UndecidableInstances #-}

{-# OPTIONS_GHC -Wall #-}

#include "ConCat/Ops.inc"

-- | Commutative monoid intended to be used with a multiplicative monoid

module ConCat.Additive where

import Prelude hiding (zipWith)
import Data.Monoid (Monoid(..), Sum(..), Product(..))
import Data.Semigroup (Semigroup(..))
import Data.Complex hiding (magnitude)
import Data.Ratio
import Foreign.C.Types (CSChar, CInt, CShort, CLong, CLLong, CIntMax, CFloat, CDouble)
import GHC.Generics (U1(..),Par1(..),(:*:)(..),(:.:)(..))
import GHC.TypeLits (KnownNat)

import Data.Key(Zip(..))
import Data.Pointed
import Data.Functor.Rep (Representable(..))
import Data.Vector.Sized (Vector)
import Data.Finite.Internal

import ConCat.Misc
import ConCat.Rep (HasRep(abst),inAbst,inAbst2)
import qualified ConCat.Rep
import ConCat.Orphans ()

-- | Commutative monoid intended to be used with a multiplicative monoid
class Additive a where
  zero  :: a
  infixl 6 ^+^
  (^+^) :: a -> a -> a
  default zero :: (Pointed h, Additive b, a ~ h b) => a
  zero = pointNI zero
  default (^+^) :: (Zip h, Additive b, a ~ h b) => Binop a
  (^+^) = zipWithNI (^+^)
  {-# INLINE zero #-}
  {-# INLINE (^+^) #-}

-- These definitions and the corresponding Catify rewrites in AltCat prevent the point and zipWith methods from getting inlined too soon.
-- See 2018-04-09 notes.
pointNI :: Pointed h => a -> h a
pointNI = point
{-# INLINE [0] pointNI #-}

zipWithNI :: Zip h => (a -> b -> c) -> h a -> h b -> h c
zipWithNI = zipWith
{-# INLINE [0] zipWithNI #-}

instance Additive () where
  zero = ()
  () ^+^ () = ()

#define ScalarType(t) \
  instance Additive (t) where { zero = 0 ; (^+^) = (+) }

ScalarType(Int)
ScalarType(Integer)
ScalarType(Float)
ScalarType(Double)
ScalarType(CSChar)
ScalarType(CInt)
ScalarType(CShort)
ScalarType(CLong)
ScalarType(CLLong)
ScalarType(CIntMax)
ScalarType(CDouble)
ScalarType(CFloat)

instance Integral a => Additive (Ratio a) where
  zero  = 0
  (^+^) = (+)

instance (RealFloat v, Additive v) => Additive (Complex v) where
  zero  = zero :+ zero
  (^+^) = (+)

-- The 'RealFloat' constraint is unfortunate here. It's due to a
-- questionable decision to place 'RealFloat' into the definition of the
-- 'Complex' /type/, rather than in functions and instances as needed.

instance (Additive u,Additive v) => Additive (u,v) where
  zero              = (zero,zero)
  (u,v) ^+^ (u',v') = (u^+^u',v^+^v')

instance (Additive u,Additive v,Additive w)
    => Additive (u,v,w) where
  zero                   = (zero,zero,zero)
  (u,v,w) ^+^ (u',v',w') = (u^+^u',v^+^v',w^+^w')

instance (Additive u,Additive v,Additive w,Additive x)
    => Additive (u,v,w,x) where
  zero                        = (zero,zero,zero,zero)
  (u,v,w,x) ^+^ (u',v',w',x') = (u^+^u',v^+^v',w^+^w',x^+^x')

type AddF f = (Pointed f, Zip f)

instance KnownNat n => Additive (Finite n) where
  zero = 0
  (^+^) = (+)

#if 1

#define AdditiveFunctor(f) instance (AddF (f), Additive v) => Additive ((f) v)

AdditiveFunctor((->) a)
AdditiveFunctor(Sum)
AdditiveFunctor(Product)
AdditiveFunctor(U1)
AdditiveFunctor(Par1)
AdditiveFunctor(f :*: g)
AdditiveFunctor(g :.: f)

#else

instance Additive v => Additive (a -> v)
instance Additive v => Additive (Sum     v)
instance Additive v => Additive (Product v)

instance Additive v => Additive (U1 v)
instance Additive v => Additive (Par1 v)
instance (Additive v, AddF f, AddF g) => Additive ((f :*: g) v)
instance (Additive v, AddF f, AddF g) => Additive ((g :.: f) v)

#endif

-- instance (Eq i, Additive v) => Additive (Arr i v) where
--   zero = point zero
--   as ^+^ bs = fmap (uncurry (^+^)) (zipC (as,bs))

-- TODO: Define and use zipWithC (^+^) as bs.

instance (Additive v, KnownNat n) => Additive (Vector n v)

-- Maybe is handled like the Maybe-of-Sum monoid
instance Additive a => Additive (Maybe a) where
  zero                = Nothing
  Nothing ^+^ b'      = b'
  a' ^+^ Nothing      = a'
  Just a' ^+^ Just b' = Just (a' ^+^ b')

-- -- Memo tries
-- instance (HasTrie u, Additive v) => Additive (u :->: v) where
--   zero  = pure   zero
--   (^+^) = liftA2 (^+^)

-- Experiment
instance Additive Bool where
  -- zero = undefined
  -- _ ^+^ _ = undefined
  zero = False
  (^+^) = (||)
  {-# INLINE zero #-}
  {-# INLINE (^+^) #-}

{--------------------------------------------------------------------
    Monoid wrapper
--------------------------------------------------------------------}

-- | Monoid under group addition.  Alternative to the @Sum@ in
-- "Data.Monoid", which uses 'Num' instead of 'Additive'.
newtype Add a = Add { getAdd :: a }
  deriving (Eq, Ord, Read, Show, Bounded)

instance HasRep (Add a) where
  type Rep (Add a) = a
  abst = Add
  repr = getAdd

instance Functor Add where fmap = inAbst

instance Applicative Add where
  pure  = abst
  (<*>) = inAbst2 ($)

instance Additive a => Semigroup (Add a) where
  (<>) = inAbst2 (^+^)

instance Additive a => Monoid (Add a) where
  mempty  = abst zero
  mappend = (<>)

instance Additive a => Additive (Add a) where
  zero  = mempty
  (^+^) = mappend

-- sumA' :: (Foldable h, Additive a) => h a -> a
-- sumA' = getAdd . foldMap Add

-- Enables translation of sumA to jamPF in AltCat.
type SummableF h = (Representable h, Eq (Rep h), Zip h, Pointed h, Foldable h)

class    SummableF h => Summable h
instance SummableF h => Summable h

-- The constraint ‘Representable h’
--   is no smaller than the instance head
-- (Use UndecidableInstances to permit this)

sumA :: (
  -- Summable h, Additive a
  Foldable h, Additive a
  ) => h a -> a
sumA = getAdd . foldMap Add
{-# OPINLINE sumA #-}

max_stars_repo_name: kenranunderscore/concat

max_issues_repo_name: kenranunderscore/concat

max_forks_repo_name: kenranunderscore/concat


text:
{-# LANGUAGE FlexibleContexts #-}
-- |
-- Module      : Data.Array.Accelerate.Math.Qtfd
-- Copyright   : [2017] Rinat Stryungis
-- License     : BSD3
--
-- Maintainer  : Rinat Stryungis <lazybonesxp@gmail.com>
-- Stability   : experimental
-- Portability : non-portable (GHC extensions)
--
-- Computation of quadratic time-frequency distributions using the accelerate-fft library.
--
-- This module uses the accelerate-fft library. And the base implementation of fft
-- uses a naive divide-and-conquer fft implementation
-- whose absolute performance is appalling. It also requires that you know on
-- the Haskell side the size of the data being transformed, and that this is
-- a power-of-two in each dimension.
--
-- For performance, compile accelerate-fft against the foreign library bindings (using any
-- number of '-fllvm-ptx', and '-fllvm-cpu' for the accelerate-llvm-ptx, and
-- accelerate-llvm-native backends, respectively), which have none of the above
-- restrictions.
-- Both of this flags are enabled by default.

module Data.Array.Accelerate.Math.Qtfd where

import Data.Array.Accelerate as A
import qualified Data.Array.Accelerate.Data.Complex as ADC
import qualified Data.Array.Accelerate.Math.BornJordan as BJ
import qualified Data.Array.Accelerate.Math.ChoiWilliams as CW
import qualified Data.Array.Accelerate.Math.FFT as AMF
import qualified Data.Array.Accelerate.Math.ModifiedB as EMB
import qualified Data.Array.Accelerate.Math.PseudoWigner as P
import Data.Array.Accelerate.Math.Wigner'

import Data.Complex
import Data.List
import Math.Gamma
-- | Wigner-ville distribution. It takes 1D array of complex floating numbers and returns 2D array of real numbers.
--  Columns of result array represents time and rows - frequency. Frequency range is from 0 to n/4, where n is a sampling frequency.

wignerVille :: (A.RealFloat e, A.IsFloating e, A.FromIntegral Int e, Elt e, Elt (ADC.Complex e), AMF.Numeric e)
  => Acc (Array DIM1 (ADC.Complex e))  -- ^ Data array
  -> Acc (Array DIM2 e)
wignerVille arr =
  let times = A.enumFromN (A.index1 leng) 0 :: Acc (Array DIM1 Int)
      leng = A.length arr
      taumx = taumaxs times
      lims = limits taumx
  in A.map ADC.real . A.transpose . AMF.fft AMF.Forward $ createMatrix arr taumx lims

wignerVilleNew :: (A.RealFloat e, A.IsFloating e, A.FromIntegral Int e, Elt e, Elt (ADC.Complex e), AMF.Numeric e)
  => Acc (Array DIM1 (ADC.Complex e))  -- ^ Data array
  -> Acc (Array DIM2 e)
wignerVilleNew arr =
  let times = A.enumFromN (A.index1 leng) 0 :: Acc (Array DIM1 Int)
      leng = A.length arr
      taumx = taumaxs times
      lims = limits taumx
  in A.map ADC.real . A.transpose . AMF.fft AMF.Forward $ createMatrix arr taumx lims

-- | Pseudo Wigner-ville distribution.
-- It takes 1D array of complex floating numbers, window and returns 2D array of real numbers.
-- Columns of result array represents time and rows - frequency. Frequency range is from 0 to n/4, where n is a sampling frequency.

pWignerVille :: (A.RealFloat e, A.IsFloating e, A.FromIntegral Int e, Elt e, Elt (ADC.Complex e), AMF.Numeric e)
  => Acc (Array DIM1 e)                -- ^ Smoothing window. Length of it must be odd.
  -> Acc (Array DIM1 (ADC.Complex e))  -- ^ Data array
  -> Acc (Array DIM2 e)
pWignerVille window arr =
  let times = A.enumFromN (A.index1 leng) 0 :: Acc (Array DIM1 Int)
      leng = A.length arr
      taumx = P.taumaxs times window
      lims = P.limits taumx
  in A.map ADC.real $ A.transpose $ AMF.fft AMF.Forward $ P.createMatrix arr window taumx lims

-- | Choi-Williams distribution. It takes 1D array of complex floating numbers,
-- and returns 2D array of real numbers.
-- Columns of result array represents time and rows - frequency. Frequency range is from 0 to n/4, where n is a sampling frequency.

choiWilliams :: (A.RealFloat e, A.IsFloating e, A.FromIntegral Int e, Elt e, Elt (ADC.Complex e), AMF.Numeric e)
  => A.Exp e                           -- ^ sigma
  -> Maybe (Acc (Array DIM1 e), Acc (Array DIM1 e))
  -> A.Exp Int                         -- ^ Smoothing window over mu, must be odd and symmetrical
  -> A.Exp Int                         -- ^ Smoothing window over tau, must be odd and symmetrical
  -> Bool                              -- ^ If normalise
  -> Acc (Array DIM1 (ADC.Complex e))  -- ^ Data array
  -> Acc (Array DIM2 e)
choiWilliams sigma mWindowArrays uWindow nWindow normalise arr =
  let times = A.enumFromN (A.index1 leng) 0 :: Acc (Array DIM1 Int)
      leng = A.length arr
      taumx = taumaxs times
      lims = limits taumx
  in A.transpose $ A.map ((*4) . ADC.real) $
    AMF.fft AMF.Forward $
      CW.summedOverMu arr taumx lims sigma mWindowArrays uWindow nWindow normalise

bornJordan :: (A.RealFloat e, A.IsFloating e, A.FromIntegral Int e, Elt e, Elt (ADC.Complex e), AMF.Numeric e)
  => A.Exp e                           -- ^ sigma
  -> Maybe (Acc (Array DIM1 e), Acc (Array DIM1 e))
  -> A.Exp Int                         -- ^ Smoothing window over mu, must be odd and symmetrical
  -> A.Exp Int                         -- ^ Smoothing window over tau, must be odd and symmetrical
  -> Acc (Array DIM1 (ADC.Complex e))  -- ^ Data array
  -> Acc (Array DIM2 e)
bornJordan sigma mWindowArrays uWindow nWindow arr =
  let times = A.enumFromN (A.index1 leng) 0 :: Acc (Array DIM1 Int)
      leng = A.length arr
      taumx = taumaxs times
      lims = limits taumx
  in A.transpose $ A.map ((*2) . ADC.real) $
    AMF.fft AMF.Forward $
      BJ.summedOverMu arr taumx lims sigma mWindowArrays uWindow nWindow

eModifiedB :: (A.RealFloat e, A.IsFloating e, A.FromIntegral Int e, Elt e, Elt (ADC.Complex e), AMF.Numeric e)
  => A.Exp e                           -- ^ alpha
  -> (Acc (Array DIM1 e))              -- ^ gammas
  -> Maybe (Acc (Array DIM1 e), Acc (Array DIM1 e))
  -> A.Exp Int
  -> A.Exp Int
  -> Bool
  -> Acc (Array DIM1 (ADC.Complex e))  -- ^ Data array
  -> Acc (Array DIM2 e)
eModifiedB alpha gammas mWindowArrays uWindow nWindow normalise arr =
  let times = A.enumFromN (A.index1 leng) 0 :: Acc (Array DIM1 Int)
      leng = A.length arr
      taumx = taumaxs times
      lims = limits taumx
  in A.transpose $ A.map (*2) $ A.map ADC.real $ AMF.fft AMF.Forward $ A.transpose $ EMB.summedOverMu arr taumx lims alpha gammas mWindowArrays uWindow nWindow normalise

makeGammaArray2 :: Double -> Int -> [Double]
makeGammaArray2 beta wLength =
  let dwLength = Prelude.fromIntegral wLength :: Double
      v = [(-0.5), ((-0.5) + 1.0/dwLength)..0.5]
      s = Prelude.map (\x -> beta :+ pi*x) v
      f = Prelude.map (\x -> abs $ ((magnitude $ gamma x) Prelude.^ 2)/((gamma beta) Prelude.^ 2)) s
      n = Prelude.floor $ (Prelude.fromIntegral $ Prelude.length f)/2.0
      frst = Data.List.take (n-1) f
      secnd  = Data.List.drop n f
  in secnd Prelude.++ frst

max_stars_repo_name: Haskell-mouse/wigner-ville-accelerate

max_issues_repo_name: Haskell-mouse/wigner-ville-accelerate

max_forks_repo_name: Haskell-mouse/wigner-ville-accelerate


text:
{-# OPTIONS_GHC -fno-warn-missing-signatures #-}

module Main where

import Wyas.Parser
import Wyas.Evaluator
import Wyas.Types
import Wyas.Error

import Test.Tasty
import Test.Tasty.HUnit

import Data.Ratio ((%))
import Data.Complex (Complex((:+)))
import Data.Vector (fromList)

testTypeParsing = testGroup "type parsing tests"
  [ testCase "parse decimal point number"
        $ (extractValue . readExpr) "3"
        @?= Number 3

  , testCase "parse floating point number"
        $ (extractValue . readExpr) "42.0"
        @?= Float 42.0

  , testCase "parse rational number"
        $ (extractValue . readExpr) "3/4"
        @?= (Ratio $ 3 % 4)

  , testCase "parse complex number (decimal)"
        $ (extractValue . readExpr) "3+1i"
        @?= (Complex $ 3.0 :+ 1.0)

  , testCase "parse complex number (floating)"
        $ (extractValue . readExpr) "3.0+1.0i"
        @?= (Complex $ 3.0 :+ 1.0)

  , testCase "parse number (binary prefix)"
        $ (extractValue . readExpr) "#b1000"
        @?= Number 8

  , testCase "parse number (hex prefix)"
        $ (extractValue . readExpr) "#xFF"
        @?= Number 255

  , testCase "parse number (oct prefix)"
        $ (extractValue . readExpr) "#o7"
        @?= Number 7

  , testCase "parse char"
        $ (extractValue . readExpr) "#\\a"
        @?= Character 'a'

  , testCase "parse space char"
        $ (extractValue . readExpr) "#\\ "
        @?= Character ' '

  , testCase "parse special newline char"
        $ (extractValue . readExpr) "#\\newline"
        @?= Character '\n'

  , testCase "parse special space char"
        $ (extractValue . readExpr) "#\\space"
        @?= Character ' '

  , testCase "parse string"
        $ (extractValue . readExpr) "\"test 123\""
        @?= String "test 123"

  , testCase "parse escaped quote in string"
        $ (extractValue . readExpr) "\"\\\"\""
        @?= String "\""

  , testCase "parse quote in string"
        $ (extractValue . readExpr) "\"test \\\"123\\\" test\""
        @?= String "test \"123\" test"

  , testCase "parse backslash in string"
        $ (extractValue . readExpr) "\"te\\\\st\""
        @?= String "te\\st"

  , testCase "parse newline in string"
        $ (extractValue . readExpr) "\"te\nst\""
        @?= String "te\nst"

  , testCase "parse atom"
        $ (extractValue . readExpr) "test_func123"
        @?= Atom "test_func123"

  , testCase "parse bool (true)"
        $ (extractValue . readExpr) "#t"
        @?= Bool True

  , testCase "parse bool (false)"
        $ (extractValue . readExpr) "#f"
        @?= Bool False
  ]

testListParsing = testGroup "list parsing tests"
  [ testCase "parse simple list"
        $ (extractValue . readExpr) "(a b)"
        @?= List [Atom "a", Atom "b"]

  , testCase "parse nested list"
        $ (extractValue . readExpr) "(a (b) c)"
        @?= List [Atom "a",
                  List [Atom "b"],
                  Atom "c"]

  , testCase "parse dotted list"
        $ (extractValue . readExpr) "(a (b . c) d)"
        @?= List [Atom "a",
                  DottedList [(Atom "b")] (Atom "c"),
                  Atom "d"]

  , testCase "parse quoted dotted"
        $ (extractValue . readExpr) "(a '(b (c . d)) e)"
        @?= List [Atom "a",
                  List [Atom "quote",
                        List [Atom "b",
                              DottedList [(Atom "c")] (Atom "d")]],
                  Atom "e"]

  , testCase "parse quasiquoted list"
        $ (extractValue . readExpr) "`(a (b c))"
        @?= List [Atom "quasiquote",
                  List [Atom "a",
                        List [Atom "b", Atom "c"]]]

  , testCase "parse quasiquoted list (unquote)"
        $ (extractValue . readExpr) "`(a ,(b c) d)"
        @?= List [Atom "quasiquote",
                  List [Atom "a",
                        List [Atom "unquote",
                              List [Atom "b",
                                    Atom "c"]],
                        Atom "d"]]
  , testCase "parse vector"
        $ (extractValue . readExpr) "#(1 a 3)"
        @?= Vector (fromList [Number 1, Atom "a", Number 3])
  ]

testEvaluation = testGroup "evaluation tests"
  [ testCase "simple addition"
        $ extractValue (readExpr "(+ 2 2)" >>= eval)
        @?= Number 4
  , testCase "negative number"
        $ extractValue (readExpr "(+ 2 (-4 1))" >>= eval)
        @?= Number 2
  , testCase "nested calls"
        $ extractValue (readExpr "(+ 2 (- 4 1))" >>= eval)
        @?= Number 5
  , testCase "chained calls"
        $ extractValue (readExpr "(- (+ 4 6 3) 3 5 2)" >>= eval)
        @?= Number 3
  , testCase "check symbol?"
        $ extractValue (readExpr "(symbol? test)" >>= eval)
        @?= Bool True
  , testCase "check number?"
        $ extractValue (readExpr "(number? 3)" >>= eval)
        @?= Bool True
  , testCase "check string?"
        $ extractValue (readExpr "(string? \"test\")" >>= eval)
        @?= Bool True
  , testCase "symbol to string"
        $ extractValue (readExpr "(symbol->string test)" >>= eval)
        @?= String "test"
  ]

testBinop = testGroup "evaluate binops"
  [ testCase "less than (numbers)"
        $ extractValue (readExpr "(< 2 3)" >>= eval)
        @?= Bool True
  , testCase "greater than (numbers)"
        $ extractValue (readExpr "(> 2 3)" >>= eval)
        @?= Bool False
  , testCase "greater or equal than (numbers)"
        $ extractValue (readExpr "(>= 3 3)" >>= eval)
        @?= Bool True
  , testCase "string equality"
        $ extractValue (readExpr "(string=? \"test\" \"test\")" >>= eval)
        @?= Bool True
  , testCase "less than (strings)"
        $ extractValue (readExpr "(string<? \"abc\" \"bba\")" >>= eval)
        @?= Bool True
  , testCase "equal on lists"
        $ extractValue (readExpr "(equal? \'(1 \"2\") \'(1 2))" >>= eval)
        @?= Bool True
  ]

tests :: TestTree
tests = testGroup "Parser Tests" [ testTypeParsing
                                 , testListParsing
                                 , testEvaluation
                                 , testBinop
                                 ]

main = defaultMain tests

max_stars_repo_name: grtlr/wyas

max_issues_repo_name: grtlr/wyas

max_forks_repo_name: grtlr/wyas